ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

VMD-SSA-LSTM多维时间序列预测:MATLAB实现与避坑指南

VMD-SSA-LSTM多维时间序列预测:MATLAB实现与避坑指南

简介:基于变分模态分解(VMD)与麻雀算法(SSA)优化长短期记忆网络(LSTM)的多维时间序列预测MATLAB代码包,面向需要高精度时序预测的科研与工程场景,可应用于风电功率、电价、负荷等数据。代码包共31个文件,含22个.m脚本、6个.mat结果文件和3个.xlsx数据文件,压缩包仅247KB,轻量易用。.m脚本涵盖VMD分解、SSA参数寻优、LSTM训练与多模型结果对比,.mat保存预测结果便于复盘,.xlsx提供原始输入序列。包内集成LSTM、VMD-LSTM、VMD-SSA-LSTM三个模型,主程序、分解模块、优化模块和误差评估模块划分清晰,可直接运行并输出R2、NSE等指标,帮助理解分解与优化带来的增益。已有236人学习下载,适合希望在组合预测方向快速搭建基线并深入调参的研究者。

1. VMD-SSA-LSTM多维时间序列预测:这个 MATLAB 方案到底解决什么问题

多维时间序列预测,说直白点就是拿过去一段时间的多个变量去推未来一个或多个时刻的值。电力负荷估算、风速预报、交通流量预测这些场景,数据基本都是这种结构。直接拿原始序列喂给 LSTM 也能训练,但预测曲线经常滞后一步,误差忽大忽小,调参全看感觉。VMD-SSA-LSTM 的思路是:先用变分模态分解(VMD)把原始序列按频率拆成几个带限分量,再用麻雀搜索算法(SSA)自动搜出 LSTM 的隐藏层节点数和学习率,最后对每个分量分别预测后叠加重构。这套流程在 MATLAB 里可以完整跑通,不用把数据导到 Python 再导回来。

这个方案适合两类人:一类是做预测类课题、需要拿它和单一 LSTM 做对比的研究者;另一类是打算给现有预测系统加“信号分解 + 自动调参”模块的工程师。VMD 和 SSA 都有成熟的开源 MATLAB 实现,真正的难点不在算法本身,而在数据格式处理和参数互相配合。这篇按“为什么要分解 → 为什么要用麻雀搜索 → 怎么把三件事接起来 → 坑在哪 → 怎么验证”的顺序一步步讲,代码框架可以直接复用到你自己的数据集上。

2. VMD 变分模态分解:为什么先把序列拆开再进 LSTM

2.1 原始序列直接进 LSTM 的三个翻车点

把话说在前面:任何预测模型都绕不开序列是否平稳这个问题。风速、电价这类原始序列,均值在变、方差在变、局部还有突变,直接塞给 LSTM 的话,最典型的现象是预测值比真实值滞后一拍。原因不难理解——LSTM 训练时发现“复制上一个时刻的值”能把损失降得很低,隐藏状态学出来就成了一个近似延迟器。这种滞后在误差图上看特别明显,预测曲线和真实曲线几乎平行,只是整体向右平移了一段。

第二个问题是模态混叠。低频趋势和高频波动混在同一个序列里,LSTM 的隐藏单元既要记长期趋势又要追短期抖动,容量被无关细节吃掉。第三个问题是训练不稳:序列方差一大,梯度范数波动剧烈,学习率稍微设大一点,loss 直接变成 NaN。这三个问题叠加在一起,你在网格搜索调参时会分不清模型表现差到底是超参数问题还是数据本身的问题,返工成本很高。

这三点就是必须先做信号分解的理由。VMD(变分模态分解)在预测预处理里比 EMD 更常用的原因在于,它把分解放进了变分框架:找一组带限固有模态函数,让所有模态的估计带宽之和最小,同时满足各模态叠加后等于原始信号。每个模态在频域上有明确的中心频率和带宽,不像 EMD 那样频繁出现模态混叠。

2.2 VMD 核心参数:模态数 K 与惩罚因子 alpha

VMD 对外有四个参数,实际需要手调的只有两个:模态数 K 和惩罚因子 alpha。K 是你要把序列拆成几个分量。K 设小了,多个频率成分会挤在同一个模态里;K 设大了,又会出现虚假模态,某个分量只是一段近似纯噪声的窄带信号。alpha 是带宽惩罚因子,它控制每个模态的带宽大小:alpha 越大模态越“瘦”,包含的频率范围越窄;alpha 越小模态越“胖”,低频和高频容易混在一起。

我一般的做法是先对原始序列做一次快速傅里叶变换,数清楚频谱上有几个明显的峰,把峰的个数当作 K 的初值,然后在 K-1、K、K+1 三个值上各跑一次分解,对比各模态的中心频率是否有序递增。如果某个模态的中心频率和相邻模态贴到了一起,说明 K 偏大;如果最后一个模态的中心频率高得离谱、波形乱跳,说明 K 偏小。alpha 常用范围是 500 到 5000,习惯上从 2000 起步。

2.3 MATLAB 里运行 VMD 分解的最小代码

下面用的是 Dragomiretskiy 那版经典的 VMD 实现,函数文件放到当前工作目录即可。

% vmd_demo.m % data: Nx1 列向量,原始时间序列 % u: KxN 矩阵,每行一个模态分量 % omega: 各模态收敛后的中心频率 K = 5; % 模态数,结合 FFT 频谱峰个数调整 alpha = 2000; % 带宽惩罚因子,常用 500~5000 tau = 0; % 噪声容忍度,0 表示严格重构 tol = 1e-7; % 迭代收敛容差 DC = 0; % 第一个模态是否允许包含直流分量 init = 1; % 中心频率初始化:1 表示均匀分布 [u, u_hat, omega] = VMD(data, alpha, tau, K, DC, init, tol);

运行后 u 是 K 行、长度与原始数据一致的矩阵。把各模态画在同一张图里,重点看两处:每个分量时域上是否连续光滑,以及相邻分量频谱是否重叠。验证频谱重叠可以对每一行做一次 FFT,归一化采样率只影响横轴刻度,不影响判断:

% 检查相邻模态的频谱重叠 Fs = 1; N = length(data); freq = (0:N/2-1) * Fs / N; for k = 1:K spec = abs(fft(u(k,:))); plot(freq, spec(1:N/2)); hold on; end hold off;

模态之间谱峰隔着明显距离,说明 K 设置合理;两个峰粘在一起,把 K 减一再试。这一步是典型的血泪经验,别指望一次调到位,多跑几轮分解花不了多少时间,但能给后面的 LSTM 训练省掉大量返工成本。还有一种常见误用是有人会把分解出的高频分量直接丢弃,理由是它代表噪声。我不建议在预测任务里轻易丢模态,风速、电价这类数据的“高频分量”里往往带着真实突变信息,LSTM 完全有能力学会它的统计规律。真正该丢的是重复的虚假模态,而不是高频分量本身。这个判断在对比实验里很容易验证:保留全部模态重构的预测误差,通常比丢掉高频分量后再重构要小。

3. SSA 麻雀算法优化 LSTM:让超参数不再靠玄学

3.1 为什么不用网格搜索和随机搜索

LSTM 里有几组超参数需要人为指定:隐藏层节点数、初始学习率、mini-batch 大小、dropout 比例、训练轮数。它们之间是相互耦合的,比如隐藏层节点变多时通常要配一个更小的学习率。网格搜索的麻烦在于,它对每组参数组合都要完整训练一次 LSTM,搜索空间稍微放大,训练时长就爆炸。随机搜索比网格好一些,但它本质上是独立采样,没有利用“上一批参数组合的效果”来指导下一次采样方向。

麻雀算法(SSA)的本质是一个种群迭代算法。麻雀群体里,发现者负责大范围搜索食物,加入者跟着发现者走并伺机争夺,同时还有一小部分警戒者监视危险区域。三者位置更新互相牵制,使得种群不会过早全部聚到某个局部位置。对 LSTM 调参来说,这意味着前几次迭代可能就会碰到一个不错的参数组合,后面迭代在这个组合附近做精细挖掘,而不是从头瞎撞。

3.2 个体编码与适应度函数设计

用 SSA 优化 LSTM,第一步是设计个体编码。一个个体就是一个参数组合向量,最常见的编码是二维:[隐藏层节点数, 初始学习率]。隐藏层节点取值范围给 [10, 200],学习率取值范围给 [1e-4, 1e-2]。学习率这种跨量级的参数,在位置更新时最好做对数变换,否则麻雀位置稍微一跳,学习率就变了好几个量级,训练直接震荡。窗口长度这种超参数不放进优化向量,因为改动它需要重新构造数据集,代价大,我一般手动定。

适应度函数决定整个搜索的方向。常见做法是把训练集按时间顺序切出一段验证集,每次用训练集的一部分训练 LSTM,返回验证集上的 RMSE 作为适应度。验证集只能用来评估,不能参与训练,否则搜索出来的参数天然对验证集过拟合,测试集表现会差一大截。下面这个函数框架可以直接放到工作目录里:

% lstmFitness.m % 输入:hidden 隐藏层节点数,lr 初始学习率 % 输出:验证集上的 RMSE function rmse = lstmFitness(hidden, lr, XTrain, YTrain, XVal, YVal) numFeatures = size(XTrain{1}, 1); layers = [ sequenceInputLayer(numFeatures) lstmLayer(hidden, 'OutputMode', 'last') fullyConnectedLayer(32) dropoutLayer(0.2) fullyConnectedLayer(1) regressionLayer]; options = trainingOptions('adam', ... 'MaxEpochs', 30, ... 'MiniBatchSize', 64, ... 'InitialLearnRate', lr, ... 'Verbose', 0); net = trainNetwork(XTrain, YTrain, layers, options); YPred = predict(net, XVal, 'MiniBatchSize', 64); rmse = sqrt(mean((YPred - YVal).^2)); end

这里有三个细节。第一,MaxEpochs 只给 30,这是刻意控制的,搜索阶段只需要参数之间的相对好坏,不需要每个组合都完全收敛。第二,dropout 加在隐藏层后面,防止对验证集过度拟合。第三,predict 函数默认输出序列最后一个时间步的预测值,回归层拆出来的目标要和 YVal 严格对齐,否则 RMSE 算出来是错的方向。

3.3 SSA 主循环与位置更新框架

SSA 在 MATLAB 里的主循环并不复杂,核心是三种角色的位置更新规则。发现者向食物多的方向移动,加入者向最优发现者靠拢,警戒者发现危险时随机跳跃。我在代码里只展示主循环和适应度调用,位置更新公式按标准 SSA 论文实现就好:

% ssa_lstm.m 节选 lb = [10, 1e-4]; % 隐藏层节点下限,学习率下限 ub = [200, 1e-2]; % 隐藏层节点上限,学习率上限 dim = 2; nPop = 10; % 麻雀数量 maxIter = 15; % 迭代次数 X = repmat(lb, nPop, 1) + rand(nPop, dim) .* (ub - lb); gBest = zeros(1, dim); gBestFitness = inf; for iter = 1:maxIter for i = 1:nPop hidden = round(X(i, 1)); % 节点数必须取整 lr = X(i, 2); % 学习率保持连续值 fitness(i) = lstmFitness(hidden, lr, XTrain, YTrain, XVal, YVal); end [bestFit, idx] = min(fitness); if bestFit < gBestFitness gBestFitness = bestFit; gBest = X(idx, :); end X = updateSSA(X, fitness, lb, ub, iter, maxIter); end fprintf('最优参数: hidden=%d, lr=%.4f, RMSE=%.4f\n', ... round(gBest(1)), gBest(2), gBestFitness);

updateSSA 的完整展开在标准实现里很多,重点是这套循环结构不要乱改:先算适应度,再更新全局最优,最后更新位置。初始种群用均匀分布在整个搜索空间撒开,避免一开始就聚在角落。麻雀数量 10、迭代 15 次,看着很小,但一次适应度就是一次完整 LSTM 训练,本地机器跑几个小时很正常。计算资源紧张时,把 nPop 降到 6、maxIter 降到 10,结果质量下降有限,时间能省一半。

4. VMD-SSA-LSTM 集成流程:从原始数据到多维预测结果的完整链路

4.1 四步主干:分解、构造序列、逐模态训练、重构

整个流程按四步走。第一步,对目标变量做 VMD 分解,得到 K 个模态分量。第二步,对每个模态分别构造带滑动窗口的数据集。第三步,用 SSA 搜索到的超参数,对每个模态分别训练一个 LSTM 模型。第四步,把 K 个模型的预测结果逐点相加,得到最终预测值。

这里有个容易想歪的地方:分解的对象是目标变量一个序列,多维特征应该怎么处理?常见做法是保留多维特征里的全部变量作为 LSTM 输入,只对目标变量做 VMD 分解。也就是说每个模态的 LSTM 输入包含多个变量的历史窗口,输出的则是目标变量对应模态的下一个值。这样既利用了分解的干净分量,又没有丢掉多维特征里的额外信息量。

4.2 多维输入的数据集构造:窗口与转置的细节

时间序列预测在 MATLAB 里最常见的输入形式是 cell 数组,每个 cell 是一个 m × windowSize 的矩阵,m 是变量个数,windowSize 是滑动窗口长度。下面这段代码把 N×M 的原始数据矩阵切成样本:

% build_sequences.m % dataMat: N x M,N 为样本点数,M 为变量个数 % 目标变量是第 targetCol 列 windowSize = 24; XTrain = {}; YTrain = []; for t = windowSize+1 : N XTrain{end+1, 1} = dataMat(t-windowSize+1 : t, :).'; YTrain(end+1, 1) = dataMat(t, targetCol); end

这段代码里的转置要特别留意。MATLAB 的 sequenceInputLayer 接收的是“特征 × 时间步”的尺寸,而原始数据切片出来的是“时间步 × 特征”,所以取到 windowSize 行数据后要做一次转置,得到 m × windowSize。初学者最容易在这里翻车——少一次转置,训练时维度报错;更隐蔽的是,维度能跑通但特征和时间步的顺序是反的,整个模型学到的东西没有意义,要到最后对比结果时才发现。

4.3 逐模态训练与预测:统一超参数,最后叠加

对每个模态建立数据集后,训练和预测的逻辑完全一致。把 SSA 搜出的隐藏层节点数和学习率套到每个模态的 LSTM 里,分别训练、分别预测:

% train_all_modes.m 节选 pred_sum = zeros(size(YTest, 1), 1); for k = 1:K [XTrain_k, YTrain_k] = buildSequences(u(k,:).', exogData, windowSize, targetCol); [XTest_k, YTest_k] = buildSequences(uTest(k,:).', exogTestData, windowSize, targetCol); layers = [ sequenceInputLayer(numFeatures) lstmLayer(hidden, 'OutputMode', 'last') fullyConnectedLayer(32) dropoutLayer(0.2) fullyConnectedLayer(1) regressionLayer]; options = trainingOptions('adam', ... 'MaxEpochs', 80, ... 'MiniBatchSize', 32, ... 'InitialLearnRate', lr, ... 'Verbose', 0); net_k = trainNetwork(XTrain_k, YTrain_k, layers, options); pred_k = predict(net_k, XTest_k, 'MiniBatchSize', 32); pred_sum = pred_sum + pred_k; end

执行顺序上有一个原则:先统一训练、统一预测,最后再叠加。不要训练完一个模态立刻去抠这个模态的误差,那样你会花大量时间在单个模态的局部调参上。正确思路是所有分量跑完后看重构预测的整体误差,再决定要调 K 还是调 LSTM 超参数。单个模态的误差只用来诊断,不是优化目标。

4.4 预测滞后:先分清是模型问题还是数据问题

重构出来的预测曲线如果还是滞后,先不要急着换模型。把原始序列的测试段和预测段对齐,看滞后的时间步数:接近 1 个时间步而且曲线形状高度一致,多半是模型学到了“复制上一步”的规律,此时优先调学习率衰减策略,或者在训练选项里加 ValidationPatience 提前停止,防止过度拟合训练集。滞后接近滑动窗口长度,说明模型对窗口内信息利用不充分,考虑把 LSTM 输出模式从 last 改成 sequence,再接一个全局平均池化层。

5. 避坑:VMD-SSA-LSTM 的 5 个常见问题与排查

5.1 K 值调错,分解结果出现虚假模态

现象:分解后某个模态的波形几乎是一条直线,或者两个相邻模态的幅度一大一小但频率接近。

原因:K 设得比真实频率成分数量大,VMD 把某一段信号硬分成了两个假模态。

解决:把 K 减一重新分解,对比中心频率顺序。大多数风速、电价序列的 K 在 3 到 7 之间。也可以用“相邻中心频率是否均匀拉开”作为体检指标,相邻中心频率差距在 10% 以内,大概率 K 偏大。这个判断要结合原始序列的频谱来看,如果原始频谱本来就是一个宽峰,K 值偏大造成的虚假模态几乎无法完全消除,此时优先选峰值最明显的 K 值。

5.2 端点效应让预测起点偏差放大

现象:每个模态序列两端有一小段明显异常值,重构预测的开头阶段误差特别大。

原因:VMD 在序列边界处信息不足,模态在端点处的收敛不自然,这种误差会从训练集末尾直接传导到测试起始点。

解决:分解前先对序列两端做镜像延拓,分解后再裁剪掉延拓部分。具体做法是把原始序列前后各延长一段镜像数据,VMD 分解后按原索引范围切回。另一个笨办法是训练时把测试序列往前多留一点窗口,让模型看到的数据边界落在已经延拓并裁剪过的区域内。这里注意,延拓长度一般取原始序列长度的 5% 到 10%,太长会把伪造的数据当真实数据用,太短又起不到压端点效应的作用。

5.3 归一化泄露:把未来信息带进了训练集

现象:验证集误差非常漂亮,测试集误差突然放大好几倍。

原因:不少代码把整个数据序列一次性做归一化,训练、验证、测试共用同一组均值和方差。对时间序列来说这就是泄露,未来区间的统计信息已经进过训练阶段了。

解决:只用训练集的均值和方差去归一化验证集和测试集。MATLAB 里把 fit 和 transform 分开写:

% normalize_split.m mu = mean(XTrain); sigma = std(XTrain); XTrain = (XTrain - mu) ./ sigma; XVal = (XVal - mu) ./ sigma; XTest = (XTest - mu) ./ sigma;

改完这个,测试集误差往往会明显变差——那个“变差”才是真实水平。做论文对比实验时,这一步的写法直接决定结论是否站得住,也是审稿人最容易盯的地方。

5.4 SSA 早熟收敛,找到的“最优参数”并不稳定

现象:每次跑 SSA 得到的最终参数都不一样,验证集 RMSE 差别不大,但测试集 RMSE 波动很大。

原因:麻雀种群数量小、迭代次数少,搜索过程过早收敛到某个局部区域,而这块区域对验证集恰好友好,对测试集却一般的很。

解决:加大初始化种群数量比加大迭代次数更有效。另外把学习率取对数后再编入个体位置,能让搜索在量级上更平滑。如果每次搜索耗时太长,可以先用一半数据量做适应度评估,找到大致区域后,再用全量数据做一次细搜。还有一个小技巧:SSA 跑完别直接用最优个体,把最后几代排名前五的个体各训一次,选测试集表现最好的那个。这样做相当于加了二次选择,代价是额外几次训练,但参数稳定性会好很多。

5.5 忽略了高频分量的真实价值

现象:把高频分量去掉之后,重构预测在突变段明显跟不上真实值,整体误差反而变大。

原因:高频分量看着像噪声,但真实数据里的短时突发信息往往就在这里面。直接丢掉等于把模型对突变的响应能力也扔掉了。

解决:先用保留全部模态的结果做基线,再尝试去掉最高频模态做对比。如果去掉后总体误差下降,说明该分量以噪声为主,可以丢弃;如果只是突变段变差,保留它。这个对比实验只需要训练一次 LSTM 的时间,但能避免你凭感觉丢掉有用的信息。

6. 验证结果与一个进阶技巧:预测效果到底怎么判断

6.1 三个误差指标:RMSE、MAE、MAPE

判断预测效果好坏,别只看单个指标。RMSE 对大幅误差敏感,适合判断有没有离谱的预测点;MAE 反映平均绝对偏差;MAPE 在多个数据集之间对比时最直观,但遇到真实值接近 0 的时刻会飙升,需要注意。画预测值和真实值的曲线图同样重要,曲线能看出滞后步数和突变段的拟合度,这是数值指标看不出来的。

6.2 与单一 LSTM 模型的对比验证

做对比实验时至少要跑三个模型:单一 LSTM、VMD-LSTM(不加 SSA)、VMD-SSA-LSTM。保证数据划分和滑动窗口完全一致,只用 SSA 搜索出的参数差别来体现优化效果。对比时把三个模型的测试集 RMSE、MAE、MAPE 放在同一张表里,最好再画一条各时间步误差的折线。VMD 分解带来的提升通常体现在整体误差下降,而 SSA 优化的提升更多体现在参数对不同模态的适应能力上,两者分开评估才能说清楚各自贡献了多少。

6.3 一个进阶技巧:按频率分组训练

如果实验结果还差一口气,可以试试把 K 个模态按中心频率分成低频组和高频组:低频组用两层 LSTM,高频组用单层 LSTM 加更强的 dropout。这样做的理由是低频分量趋势性强,需要更大的模型容量去记忆长期依赖;高频分量随机性强,浅模型反而更不容易过拟合。代价是训练时间接近翻倍,所以放在最后一步去尝试。我自己的习惯是每次跑完实验,先把 VMD 分解的模态图仔仔细细看一遍,再核对测试集的归一化是否按训练集参数执行,最后才看误差指标。顺序一乱,调参的时间基本就白费了,希望帮到你。

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进