以下是使用Python实现EMD-SSA-VMD-LSTM混合模型进行时间序列预测的完整代码，结合经验模态分解（EMD）、奇异谱分析（SSA）、变分模态分解（VMD）与LSTM深度学习模型。该方案适用于复杂非平稳信号的预测任务，代码包含数据生成、多级分解、模型构建和结果可视化。

---

一、环境配置与依赖

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from PyEMD import EMD
from vmdpy import VMD
from scipy.linalg import hankel, svd
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_absolute_error, mean_squared_error
import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader

# 设置随机种子保证可复现性
np.random.seed(42)
torch.manual_seed(42)

---

二、数据生成（示例数据）

def generate_complex_signal(length=1000):
    """生成含多成分的非平稳信号"""
    t = np.linspace(0, 10, length)
    
    # 趋势项
    trend = 0.02 * t**2 + 0.1 * t
    
    # 周期成分
    seasonal1 = 1.5 * np.sin(2 * np.pi * 0.8 * t)
    seasonal2 = 0.8 * np.cos(2 * np.pi * 2.5 * t)
    
    # 脉冲噪声
    impulse = np.zeros(length)
    impulse[np.random.choice(length, 20)] = np.random.uniform(-3, 3, 20)
    
    # 高斯噪声
    noise = 0.3 * np.random.randn(length)
    
    return trend + seasonal1 + seasonal2 + impulse + noise

# 生成数据并可视化
data = generate_complex_signal()
plt.figure(figsize=(12,4))
plt.plot(data, color='darkblue')
plt.title("Generated Non-stationary Signal")
plt.show()

---

三、多级信号分解

1. 经验模态分解（EMD）

def emd_decomposition(signal):
    emd = EMD()
    imfs = emd(signal)
    return imfs

imfs_emd = emd_decomposition(data)
print(f"EMD分解得到 {imfs_emd.shape[0]} 个IMF分量")

2. 奇异谱分析（SSA）

def ssa_decomposition(signal, window=30, rank=3):
    """奇异谱分析核心函数"""
    # 构建轨迹矩阵
    L = window
    K = len(signal) - L + 1
    X = hankel(signal[:L], signal[L-1:])
    
    # 奇异值分解
    U, S, VT = svd(X, full_matrices=False)
    
    # 选择主成分重构
    X_rank = (U[:, :rank] * S[:rank]) @ VT[:rank, :]
    
    # 对角平均化
    reconstructed = np.zeros_like(signal)
    for i in range(len(signal)):
        X_diag = np.diagonal(X_rank, offset=-(L-1-i))
        reconstructed[i] = X_diag.mean() if X_diag.size > 0 else 0
    return reconstructed

# 对每个EMD-IMF执行SSA分解
components_ssa = []
for imf in imfs_emd:
    ssa_comp = ssa_decomposition(imf, window=30, rank=3)
    components_ssa.append(ssa_comp)

3. 变分模态分解（VMD）

def vmd_decomposition(signal, alpha=2000, K=4):
    u, _, _ = VMD(signal, alpha=alpha, tau=0, K=K, DC=0, init=1, tol=1e-7)
    return u

# 对SSA结果进行VMD分解
final_components = []
for comp in components_ssa:
    vmd_comps = vmd_decomposition(comp, K=2)
    final_components.extend(vmd_comps)

# 合并所有分量
all_components = np.vstack(final_components)
print(f"总分解分量数: {all_components.shape[0]}")

---

四、数据预处理

1. 归一化处理

scalers = []
scaled_components = []
for comp in all_components:
    scaler = MinMaxScaler(feature_range=(-1, 1))
    scaled = scaler.fit_transform(comp.reshape(-1, 1)).flatten()
    scaled_components.append(scaled)
    scalers.append(scaler)

scaled_components = np.array(scaled_components)

2. 数据集构建

class HybridDataset(Dataset):
    def __init__(self, components, lookback=60, horizon=1):
        self.components = components
        self.lookback = lookback
        self.horizon = horizon
        
    def __len__(self):
        return self.components.shape[1] - self.lookback - self.horizon + 1
    
    def __getitem__(self, idx):
        x = self.components[:, idx:idx+self.lookback].T  # (lookback, n_components)
        y = self.components[:, idx+self.lookback:idx+self.lookback+self.horizon].T
        return torch.FloatTensor(x), torch.FloatTensor(y)

lookback = 60
horizon = 1
dataset = HybridDataset(scaled_components, lookback, horizon)
dataloader = DataLoader(dataset, batch_size=32, shuffle=True)

---

五、混合LSTM模型

1. 模型架构

class MultiScaleLSTM(nn.Module):
    def __init__(self, input_size, hidden_size=128, output_size=1):
        super().__init__()
        
        # 特征提取层
        self.lstm1 = nn.LSTM(input_size, hidden_size, batch_first=True)
        self.dropout1 = nn.Dropout(0.3)
        
        # 时序预测层
        self.lstm2 = nn.LSTM(hidden_size, hidden_size//2, batch_first=True)
        self.dropout2 = nn.Dropout(0.2)
        
        # 输出层
        self.fc = nn.Linear(hidden_size//2, output_size)
        
    def forward(self, x):
        # 输入形状: (batch_size, seq_len, input_size)
        out, (h, c) = self.lstm1(x)
        out = self.dropout1(out)
        
        out, _ = self.lstm2(out)
        out = self.dropout2(out[:, -1, :])
        
        return self.fc(out)

2. 模型训练

model = MultiScaleLSTM(input_size=scaled_components.shape[0])
criterion = nn.MSELoss()
optimizer = torch.optim.AdamW(model.parameters(), lr=0.001, weight_decay=1e-4)

# 训练循环
for epoch in range(100):
    total_loss = 0
    for x, y in dataloader:
        optimizer.zero_grad()
        pred = model(x)
        loss = criterion(pred, y.squeeze())
        loss.backward()
        optimizer.step()
        total_loss += loss.item()
    
    print(f"Epoch {epoch+1}/100 | Loss: {total_loss/len(dataloader):.4f}")

---

六、预测与结果重构

1. 多步预测

def recursive_forecast(model, initial_seq, steps=50):
    current_seq = initial_seq.clone()
    predictions = []
    
    for _ in range(steps):
        with torch.no_grad():
            pred = model(current_seq.unsqueeze(0))
            predictions.append(pred.numpy()[0][0])
            
            # 更新输入序列
            current_seq = torch.cat([current_seq[1:], pred.unsqueeze(0)])
    
    return np.array(predictions)

# 获取初始序列
test_input = scaled_components[:, -lookback:]
test_input = torch.FloatTensor(test_input.T)  # (lookback, n_components)

# 执行预测
pred_steps = 50
prediction = recursive_forecast(model, test_input, pred_steps)

2. 结果反归一化

# 重构所有分量预测
pred_components = []
for i in range(len(scalers)):
    pred_scaled = prediction * 0  # 初始化
    pred_scaled[i::len(scalers)] = prediction  # 分量位置插值
    pred_components.append(scalers[i].inverse_transform(pred_scaled.reshape(-1, 1)))

# 合成最终结果
final_pred = np.sum(pred_components, axis=0).flatten()

# 获取真实值
true_values = data[-pred_steps:]

---

七、性能评估与可视化

1. 评估指标

mae = mean_absolute_error(true_values, final_pred)
rmse = np.sqrt(mean_squared_error(true_values, final_pred))
print(f"MAE: {mae:.4f}")
print(f"RMSE: {rmse:.4f}")

2. 结果可视化

plt.figure(figsize=(12,6))
plt.plot(true_values, label='True', marker='o', linestyle='--')
plt.plot(final_pred, label='Predicted', marker='x', linewidth=2)
plt.fill_between(range(len(final_pred)), 
                 final_pred - 1.96*rmse, 
                 final_pred + 1.96*rmse, 
                 alpha=0.2, color='orange')
plt.title("EMD-SSA-VMD-LSTM Multi-step Prediction")
plt.legend()
plt.grid(True)
plt.show()

---

八、完整数据流说明

步骤	技术实现	数学表达
信号生成	合成趋势项+周期项+噪声	$$x(t)=\sum _{i=1}^n{a}_i{f}_i(t)+ϵ(t)$$
EMD分解	自适应分解非平稳信号	$$x(t)=\sum _{k=1}^K{c}_k(t)+r(t)$$
SSA分解	轨迹矩阵SVD分解	$$\mathbf{X}={\mathbf{U}\mathbf{\Sigma }\mathbf{V}}^T$$
VMD分解	变分模态优化分解	min⁡{uk},{ωk}∑k∥∂t[uk(t)e−jωkt]∥22\min_{\{u_k\},\{\omega_k\}} \sum_k \|\partial_t[u_k(t)e^{-j\omega_k t}]\|_2^2{uk​},{ωk​}min​k∑​∥∂t​[uk​(t)e−jωk​t]∥22​
特征融合	多分量时序对齐	$${\mathbf{X}}_{\text{stack}}=[{\mathbf{C}}_1^T;{\mathbf{C}}_2^T;\dots ;{\mathbf{C}}_n^T]$$
LSTM建模	门控机制时序建模	$$f_t=\sigma (W_f\cdot [h_{t-1},x_t]+b_f)$$
结果重构	逆归一化加权求和	$$\hat{y}=\sum _{k=1}^K{\text{scaler}}_k^{-1}({\hat{c}}_k)$$

---

九、参数调优建议

参数	优化策略	典型值范围
EMD最大IMF数	根据信号复杂度调整	5-10
SSA窗口长度	取1/3周期长度	20-50
VMD模态数(K)	频谱分析确定	3-6
LSTM隐藏层	防止过拟合	64-256
学习率	余弦退火调整	1e-4~1e-3
输入序列长度	覆盖主要周期	60-120

---

十、扩展方向

1. 自适应分解

   # 自动确定VMD的K值
   from vmdpy import VMD
   def auto_vmd(signal, max_K=8):
       for K in range(3, max_K+1):
           u, _, _ = VMD(signal, alpha=2000, K=K)
           if np.any(np.isnan(u)):
               return K-1
       return max_K
   

2. 概率预测

   # 修改输出层为分位数回归
   self.fc = nn.Linear(hidden_size//2, 3)  # 输出3个分位数
   

3. 在线学习

   # 增量训练机制
   def online_update(model, new_data):
       model.train()
       optimizer.zero_grad()
       outputs = model(new_data)
       loss = criterion(outputs, targets)
       loss.backward()
       optimizer.step()
   

---

源码说明

1. 数据兼容性

   • 支持CSV输入：修改generate_complex_signal()为pd.read_csv()
   • 多变量扩展：调整输入维度为(n_features, seq_len)

2. 性能优化

   • 启用CUDA加速：model.to('cuda')
   • 使用混合精度训练：scaler = torch.cuda.amp.GradScaler()

3. 工业级部署

   # 模型保存与加载
   torch.save(model.state_dict(), 'multiscale_lstm.pth')
   model.load_state_dict(torch.load('multiscale_lstm.pth'))
   

---

该方案通过三级分解（EMD-SSA-VMD）充分提取信号多尺度特征，结合深度LSTM建模复杂时序依赖，在非平稳信号预测中展现出显著优势。实际应用时需根据数据特性调整分解参数与模型结构，并通过误差分析持续优化。