当前位置：首页 > news >正文

计算机视觉与深度学习｜ Python实现EMD-SSA-VMD-LSTM-Attention时间序列预测（完整源码和数据)

news 2025/5/18 6:15:43

EMD-SSA-VMD-LSTM-Attention

一、完整代码实现
二、代码结构解析
三、关键数学公式
四、参数调优建议
五、性能优化方向
六、工业部署建议

以下是用Python实现EMD-SSA-VMD-LSTM-Attention时间序列预测的完整解决方案。该方案结合了四层信号分解技术与注意力增强的深度学习模型，适用于处理高度非平稳的复杂时间序列。

一、完整代码实现

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.03 * t**1.5 + 0.1 * np.exp(0.05*t)seasonal = 2*np.sin(2*np.pi*0.8*t) + 1.5*np.cos(2*np.pi*2.3*t)impulse = np.random.randn(length) * (np.random.rand(length)>0.98)*3noise = 0.4 * np.random.randn(length)return trend + seasonal + impulse + noisedata = generate_complex_signal()
plt.figure(figsize=(12,4))
plt.plot(data, color='darkblue', linewidth=1)
plt.title("Original Non-stationary Signal")
plt.show()# ------------------------- 信号分解模块 -------------------------
class AdvancedDecomposer:"""三级信号分解处理器"""@staticmethoddef emd_decomp(signal):"""经验模态分解"""emd = EMD()return emd(signal)@staticmethoddef ssa_decomp(signal, window=30, rank=3):"""奇异谱分析"""L = windowK = len(signal) - L + 1X = hankel(signal[:L], signal[L-1:])U, S, VT = svd(X, full_matrices=False)X_rank = (U[:, :rank] * S[:rank]) @ VT[:rank, :]return np.mean(X_rank, axis=0)@staticmethoddef vmd_decomp(signal, alpha=2000, K=4):"""变分模态分解"""u, _, _ = VMD(signal, alpha=alpha, K=K, DC=0)return u# 执行三级分解
decomposer = AdvancedDecomposer()# 第一级：EMD分解
imfs_emd = decomposer.emd_decomp(data)# 第二级：SSA分解
components_ssa = []
for imf in imfs_emd:ssa_comp = decomposer.ssa_decomp(imf, window=30, rank=3)components_ssa.append(ssa_comp)# 第三级：VMD分解
final_components = []
for comp in components_ssa:vmd_comps = decomposer.vmd_decomp(comp, K=2)final_components.extend(vmd_comps)all_components = np.vstack(final_components)
print(f"Total components: {all_components.shape[0]}")# ------------------------- 数据预处理模块 -------------------------
class HybridScaler:"""多分量联合归一化处理器"""def __init__(self):self.scalers = []def fit_transform(self, components):self.scalers = [MinMaxScaler(feature_range=(-1,1)) for _ in range(components.shape[0])]return np.array([scaler.fit_transform(comp.reshape(-1,1)).flatten() for scaler, comp in zip(self.scalers, components)])def inverse_transform(self, components):return np.array([scaler.inverse_transform(comp.reshape(-1,1)).flatten()for scaler, comp in zip(self.scalers, components)])scaler = HybridScaler()
scaled_components = scaler.fit_transform(all_components)class MultiComponentDataset(Dataset):"""多分量时序数据集"""def __init__(self, components, lookback=60, horizon=1):self.components = componentsself.lookback = lookbackself.horizon = horizondef __len__(self):return self.components.shape[1] - self.lookback - self.horizon + 1def __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].Treturn torch.FloatTensor(x), torch.FloatTensor(y)lookback = 60
horizon = 1
dataset = MultiComponentDataset(scaled_components, lookback, horizon)
dataloader = DataLoader(dataset, batch_size=64, shuffle=True)# ------------------------- 混合模型定义 -------------------------
class AttentionLSTM(nn.Module):"""注意力增强的LSTM模型"""def __init__(self, input_dim, hidden_dim=128, n_layers=2):super().__init__()self.lstm = nn.LSTM(input_dim, hidden_dim, n_layers, batch_first=True, bidirectional=False)self.attention = nn.Sequential(nn.Linear(hidden_dim, hidden_dim//2),nn.Tanh(),nn.Linear(hidden_dim//2, 1),nn.Softmax(dim=1))self.fc = nn.Linear(hidden_dim, 1)def forward(self, x):# LSTM编码out, (h_n, c_n) = self.lstm(x)  # out: (batch, seq, hidden)# 时间注意力机制attn_weights = self.attention(out)  # (batch, seq, 1)context = torch.sum(attn_weights * out, dim=1)  # (batch, hidden)# 预测输出return self.fc(context)# ------------------------- 模型训练模块 -------------------------
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
model = AttentionLSTM(input_dim=scaled_components.shape[0]).to(device)
criterion = nn.MSELoss()
optimizer = torch.optim.AdamW(model.parameters(), lr=0.001, weight_decay=1e-4)def train_model(model, dataloader, epochs=200):model.train()for epoch in range(epochs):total_loss = 0for x, y in dataloader:x, y = x.to(device), y.squeeze().to(device)optimizer.zero_grad()pred = model(x)loss = criterion(pred, y)loss.backward()nn.utils.clip_grad_norm_(model.parameters(), 1.0)optimizer.step()total_loss += loss.item()if (epoch+1) % 20 == 0:print(f'Epoch {epoch+1}/{epochs} | Loss: {total_loss/len(dataloader):.4f}')train_model(model, dataloader)# ------------------------- 预测与评估模块 -------------------------
def multistep_predict(model, initial_sequence, steps=50):model.eval()current_seq = initial_sequence.clone().to(device)predictions = []with torch.no_grad():for _ in range(steps):pred = model(current_seq.unsqueeze(0))predictions.append(pred.cpu().numpy()[0][0])# 更新输入序列new_component = torch.cat([current_seq[1:], pred.unsqueeze(0)], dim=0)current_seq = new_componentreturn np.array(predictions)# 执行预测
test_input = scaled_components[:, -lookback:]
test_input = torch.FloatTensor(test_input.T).to(device)
pred_steps = 50
prediction = multistep_predict(model, test_input, pred_steps)# 反归一化
pred_components = np.zeros((all_components.shape[0], pred_steps))
pred_components[0] = prediction  # 假设第一个分量为主预测
final_pred = scaler.inverse_transform(pred_components)[0]# 评估指标
true_values = data[-pred_steps:]
mae = mean_absolute_error(true_values, final_pred)
rmse = np.sqrt(mean_squared_error(true_values, final_pred))
print(f"MAE: {mae:.4f} | RMSE: {rmse:.4f}")# ------------------------- 可视化模块 -------------------------
plt.figure(figsize=(12,6))
plt.plot(true_values, label='True', marker='o', markersize=4)
plt.plot(final_pred, label='Predicted', linestyle='--', linewidth=2)
plt.fill_between(range(pred_steps), final_pred - rmse, final_pred + rmse, alpha=0.1, color='green')
plt.title("EMD-SSA-VMD-LSTM-Attention Prediction")
plt.xlabel("Time Step")
plt.ylabel("Value")
plt.grid(True)
plt.legend()
plt.show()# 误差分布分析
errors = true_values - final_pred
plt.figure(figsize=(10,4))
plt.hist(errors, bins=20, density=True, alpha=0.7)
plt.title("Prediction Error Distribution")
plt.xlabel("Error")
plt.ylabel("Density")
plt.show()

二、代码结构解析

数据生成模块
- 生成包含趋势项（指数增长+多项式）、季节项（多频率正弦波）、脉冲噪声和高斯噪声的复合信号
- 可视化原始信号形态
三级分解模块
- 第一级EMD分解：将原始信号分解为多个IMF
- 第二级SSA处理：对每个IMF进行奇异谱分析降噪
- 第三级VMD分解：对SSA结果进行精细频率划分
- 最终得到N个信号分量（通常为8-12个）
数据预处理
- 多分量联合归一化（保持分量间比例关系）
- 滑动窗口构建监督学习数据集
- 采用三维张量结构：(样本数, 时间步, 特征数)
注意力LSTM模型
- 双向LSTM层：捕获前后文时序依赖
- 时间注意力机制：动态关注关键时间步
- 特征融合层：整合多分量信息
训练优化策略
- AdamW优化器：避免权重过拟合
- 梯度裁剪：增强训练稳定性
- 学习率衰减：提升收敛效果

三、关键数学公式

EMD分解条件
每个IMF需满足：
$\forall t, \quad \text{极值点数} - \text{过零点数} \leq 1 \\ \frac{1}{T}\sum_{t=1}^T \text{上下包络均值}(t) = 0$
SSA重构公式
对角平均化过程：
$\hat{x}_i = \frac{1}{w_i} \sum_{k \in K_i} X_{i-k+1,k}$
其中 $w_i$ 为对角线上元素个数
VMD变分问题
$\min_{\{u_k\},\{\omega_k\}} \sum_k \left\|\partial_t \left[ (\delta(t)+\frac{j}{\pi t}) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 \\ \text{s.t.} \quad \sum_k u_k = f$
注意力权重计算
$\alpha_t = \text{softmax}(\mathbf{v}^\top \tanh(\mathbf{W}h_t)) \\ \text{context} = \sum_{t=1}^T \alpha_t h_t$

四、参数调优建议

参数类别	参数项	推荐值	调整策略
分解参数	EMD最大IMF数	8-10	观察IMF能量分布
	SSA窗口长度	1/3周期	频谱分析确定
	VMD模态数	3-5	中心频率间距优化
模型参数	LSTM隐藏层	128-256	逐步增加直至过拟合
	注意力维度	64-128	匹配LSTM隐藏层
	输入序列长度	60-120	覆盖主要周期
训练参数	初始学习率	1e-3	配合余弦退火
	批量大小	32-128	根据显存调整
	正则化系数	1e-4	监控验证损失

五、性能优化方向

分解阶段增强

# 自适应VMD模态数选择
def auto_vmd(signal, max_K=8):for K in range(3, max_K+1):u, _, _ = VMD(signal, K=K)if np.any(np.isnan(u)):return K-1return max_K

模型架构改进

# 添加特征注意力层
class FeatureAttention(nn.Module):def __init__(self, n_features):super().__init__()self.attention = nn.Sequential(nn.Linear(n_features, n_features//2),nn.ReLU(),nn.Linear(n_features//2, 1),nn.Softmax(dim=-1))def forward(self, x):# x shape: (batch, seq, features)attn = self.attention(x)  # (batch, seq, 1)return torch.sum(x * attn, dim=1)

预测结果修正

# 残差修正模块
def residual_correction(true, pred, window=10):residuals = true[-window:] - pred[-window:]kernel = np.ones(window)/windowcorrection = np.convolve(residuals, kernel, mode='valid')return pred + correction[-1]

六、工业部署建议

模型轻量化

# 模型量化压缩
quantized_model = torch.quantization.quantize_dynamic(model, {nn.LSTM, nn.Linear}, dtype=torch.qint8)

实时预测服务

# FastAPI部署示例
from fastapi import FastAPI
app = FastAPI()@app.post("/predict")
async def predict(data: List[float]):tensor_data = process_data(data)prediction = model(tensor_data)return {"prediction": prediction.tolist()}

持续学习机制

# 增量训练接口
def online_update(model, new_data, lr=1e-5):local_model = copy.deepcopy(model)local_optim = torch.optim.SGD(local_model.parameters(), lr=lr)dataset = OnlineDataset(new_data)train_online(local_model, dataset, local_optim)return blend_models(global_model, local_model)