简介：这份压缩包提供了一套完整的EMD-KPCA-LSTM多维时间序列预测程序，面向光伏功率预测、回归预测等场景，适合具备Matlab和深度学习基础的研究者、赛题选手或相关专业学生。案例使用北半球光伏功率数据，以太阳辐射度、气温、气压、大气湿度为输入特征，光伏功率为输出，先经经验模态分解将各环境序列分解为多个本征模态函数，降低非平稳性；再通过核主成分分析提取关键影响因子，消除相关性与冗余，降低输入维度；最后利用长短期记忆网络进行动态时间建模，实现功率预测，整套流程逻辑完整，可直接运行复现实验。资源包共24个文件，大小约2.64MB，其中包含8个m格式源程序，覆盖EMD分解、KPCA降维、LSTM训练与误差计算等环节；另有3个mat数据文件保存中间与最终结果，1个xlsx原始数据表，1个txt说明文件和1个pdf文档，方便理解算法细节与操作流程；同时附有10张png运行结果截图，可直观对比预测效果。目前已有182人学习使用，适合用作论文实验、毕业设计或课程项目的参考实现。
1. EMD-KPCA-LSTM 在光伏功率预测里到底解决了什么
做光伏功率预测时，我最初把太阳辐射度、气温、气压、大气湿度四个环境特征直接拼进一个 LSTM 网络，训练轮数不少，验证集误差却一直下不去。原因在于气象序列非平稳性太强，几个特征之间又存在复杂耦合，网络低层很难同时消化这些信息。146 号资源包给出的解法，是把“经验模态分解、核主成分分析、LSTM 网络”三段式串成一条完整预测链路：先用 EMD 把每个输入特征按时间尺度拆成多个分量，再用 KPCA 压缩出关键因子，最后交给 LSTM 做动态时序建模。源码是 MATLAB 程序，跑通后可以平移到负荷预测、风电预测这类多维时间序列预测任务。适合想完整复现或改造一套预测流程的算法从业者。
2. 先做 EMD 分解：非平稳环境量被拆成多尺度 IMF
2.1 为什么环境量要先拆开再进网络
光伏功率的波动来源是环境量，而不是功率本身。太阳辐射度序列里有分钟级云层遮挡造成的快速抖动，也有小时级太阳高度角变化带来的趋势；气温、气压、湿度同样有各自的变化周期。把这些不同频段的信号揉在一起直接作为 LSTM 输入，网络既要去学高频抖动，又要在同一组权重里表达低频趋势，任务太重。
EMD 的本质是按信号自身的时间特征做筛分，不需要预设基函数，也不依赖平稳性假设，这一点比小波变换更省心。分解之后，每个固有模态函数（IMF）只覆盖一个相对窄的频段：高频 IMF 反映快速波动，低频 IMF 和残余项则反映趋势。四个环境特征各自分解，得到的 IMF 数量还不一样，这在后面拼接时是关键坑点，稍后展开。
需要留意的是，这套流程分解的是输入特征，不是预测目标。光伏功率原始序列保留下来作为回归目标，中间不经过分解重构，链路更短，也避免了多个预测分量叠加产生的误差累积。
2.2 用 emdtest.m 跑一遍逐特征分解
资源包里的 emd.m 是核心工具，emdtest.m 可以看作它的使用示例。我一般这样组织调用流程：
%% 读取北半球光伏数据，前四列是环境特征，第五列是功率 T = readtable('北半球光伏数据.xlsx'); feat = table2array(T(:, 1:4)); % 太阳辐射度、气温、气压、大气湿度 power = table2array(T(:, 5)); % 光伏功率，作为预测目标 %% 逐特征做 EMD 分解，每个特征得到一组 IMF imfCell = cell(size(feat, 2), 1); residCell = cell(size(feat, 2), 1); for i = 1:size(feat, 2) [imfCell{i}, residCell{i}] = emd(feat(:, i)); fprintf('特征%d: 分解出 %d 个IMF\n', i, size(imfCell{i}, 2)); end这段代码把四个环境特征分别送进 emd.m，返回结果按特征存放：imfCell{i} 是一个矩阵，行数与原始序列长度一致，每一列是一个 IMF 或残余项。分解顺序默认按频率从高到低排列，第一列通常是波动最剧烈的高频分量，最后一列接近单调趋势。
EMD 内部做的是经典筛分：找极值点、用三次样条插值构造上下包络、取包络均值、原信号减均值得到候选分量，再反复迭代直到满足 IMF 条件。停止条件一般由筛分门限控制，默认门限在 0.2 左右基本够用。门限调小会让分解更充分，但容易出现虚假振荡分量；门限调大则分解粗糙，低频趋势拆不干净。光伏序列本身噪声大，我建议先按默认门限跑，看 IMF 数量再决定要不要收紧。
原始 xlsx 前几行偶尔会混入空值或异常标记。读进来之后先sum(isnan(feat))检查一遍，有缺失就按前后均值补掉，否则 EMD 在插值包络时会把空值带进结果。
2.3 IMF 数量对齐：拼接前必须统一列数
四个特征分解后，IMF 个数很难一致。辐射度序列波动复杂，可能拆出 7 个分量；气压相对平稳，可能只有 3 到 4 个。如果直接横向拼接，列数对不上，矩阵根本拼不起来。即使强行用 cell 存放，后续 KPCA 也需要一个规则矩阵。
我习惯先统计最小 IMF 数，按这个数量统一截断：
%% 统一 IMF 数量：截断到最小个数，避免补零列干扰核距离 imfCounts = cellfun(@(c) size(c, 2), imfCell); minCount = min(imfCounts); fprintf('IMF数量范围: %d ~ %d，统一取 %d\n', min(imfCounts), max(imfCounts), minCount); %% 拼接成 KPCA 的输入 featureBlock = []; for i = 1:numel(imfCell) featureBlock = [featureBlock, imfCell{i}(:, 1:minCount)]; end这里有两个选择：截断到最小值，或者补零到最大值。我选截断而不是补零，原因是补进去的全零列在后续 RBF 核计算里会产生虚假的距离贡献，等于凭空造了一批相似样本，KPCA 主成分会被带偏。截断会丢掉一些高频细节，但保留下来的 IMF 已经覆盖了各特征的主要变化尺度，代价可控。
拼接后的 featureBlock 行数还是原始样本数，列数变成 4 乘 minCount。以 4 个特征、每个 5 个 IMF 为例，这里就是 20 维。这 20 维存在明显的相关性和冗余，接下来交给 KPCA 处理。源码包里的 emd_data.mat 就对应这个阶段的时间结果，第一次跑完分解后，建议把 featureBlock 保存成同样的 mat 文件，后续调参可以直接复用。
3. KPCA 降维：核矩阵中心化与主成分取舍
3.1 为什么 PCA 不够用
EMD 拆出来的 IMF 之间不一定只有线性关系。辐射度的高频分量与气温的某段趋势可能在非线性运算下互相纠缠，PCA 只能找到线性投影方向，对这种耦合无能为力。KPCA 先通过核函数把样本映射到高维特征空间，在这个空间里做主成分分析，相当于用非线性方式提取主方向。
资源包中 kernel.m 实现的就是径向基核函数，distanceMatrix.m 负责计算样本间的欧氏距离平方。RBF 核有一个关键参数 gamma，控制核宽度，gamma 越大，核函数衰减越快，模型越容易盯着局部样本，过拟合风险也越高；gamma 太小，所有样本的核值都趋于一致，KPCA 退化成近似均值投影。光伏数据特征尺度差异大，gamma 先从 1 除以输入维数这个量级起步，比如 20 维输入就从 0.05 开始试，再根据主成分贡献率调整。
3.2 核矩阵中心化：kPCA.m 里最容易被忽视的一步
KPCA 和普通 PCA 最大的区别在于，核矩阵拿到手之后不能直接做特征分解，必须先中心化。中心化用的均值是整个核矩阵的全局均值，而不是逐列均值。漏掉这一步，主成分方向会被核空间的平移量污染，降维后的 score 分布完全走样。
%% 核矩阵构建，对应源码包 kernel.m + distanceMatrix.m D = distanceMatrix(X, X); % n×n，两两样本的欧氏距离平方 K = exp(-gamma * D); % RBF 核矩阵 %% 核矩阵中心化，这是 KPCA 的核心步骤 n = size(X, 1); oneN = ones(n) / n; Kc = K - oneN * K - K * oneN + oneN * K * oneN; %% 特征分解，取前 nPC 个主成分 [V, D] = eig((Kc + Kc') / 2); % 对称化处理，避免数值误差 lam = diag(D); [lam, idx] = sort(lam, 'descend'); V = V(:, idx); %% 特征向量归一化：score 才是真正的 KPCA 主成分 Vn = V(:, 1:nPC) ./ sqrt(lam(1:nPC))'; score = Kc * Vn; % n×nPC，降维后的特征序列代码里(Kc + Kc')/2是保障特征分解数值稳定的小技巧，矩阵理论上对称，实际计算时浮点误差会破坏对称性。特征向量除以特征值的平方根，这一步对应数学上的归一化约束，少了它，score 会混入尺度偏差，LSTM 训练收敛速度会被拖慢。
对于预测任务，测试集映射时有一个应用陷阱：测试样本的核矩阵也必须用训练集的全局均值做中心化，不能用测试核矩阵自身的均值。训练时保存 oneN 和 Kc 的均值信息，预测阶段加载后复用，否则主成分空间整体平移，预测结果会系统性偏移。这个问题的具体现象和解法放在后面排错章节。
3.3 主成分数量：按累计贡献率定，而不是拍脑袋
nPC 取多少，直接决定 LSTM 输入维度。取少了丢信息，取多了又把 KPCA 的降维意义抹掉。我按累计特征值贡献率来定，而不是固定取某个数字：
cumRatio = cumsum(lam) / sum(lam); nPC = find(cumRatio >= 0.95, 1, 'first'); fprintf('达到95%%累计贡献率需要取 %d 个主成分\n', nPC);0.95 是一个常用阈值，在时间序列场景里可以先从 0.95 开始，如果后续 LSTM 验证集误差偏高，再往上调整到 0.98。有一点需要说明：KPCA 的特征值占比并不完全等价于原始方差解释率，它描述的是映射后特征空间里的能量集中程度。在模拟项目 X 上跑，20 维左右的 IMF 特征经过 KPCA，通常能压缩到 10 到 13 个主成分，降维比例接近一半。这个压缩效果是有实际意义的，输入维度降下来之后，LSTM 的网络规模可以小一截，训练速度提升明显。
主成分个数确定后，score 矩阵替代原先的 featureBlock，成为 LSTM 的输入。资源包里的 kpcaTest.m 我理解就是为了验证这一步效果而存在的小脚本，跑一下可以看到主成分的分布情况。
4. LSTM 动态建模：从主成分序列到功率预测
4.1 时间窗构造与数据集切分顺序
KPCA 降维后得到的是 score 序列，每个时刻一个向量。LSTM 不能只吃单时刻向量，它需要看一段历史才能学习时间依赖。这里要构造滑动窗口样本：取过去 window 个时刻的主成分作为输入，取窗口末端时刻的功率作为标签。
%% 切分序列：先切数据，再构造窗口，顺序不能反 trainN = round(length(power) * 0.8); scoreTrain = score(1:trainN, :); powerTrain = power(1:trainN, :); %% 构造时间窗样本 window = 24; % 输入过去24个时刻 numSamples = trainN - window; XTrain = cell(numSamples, 1); YTrain = zeros(numSamples, 1); for t = 1:numSamples XTrain{t} = scoreTrain(t:t+window-1, :)'; % numFeat × window YTrain(t) = powerTrain(t + window - 1); end这段代码有两个关键点。第一，XTrain 是 cell 数组，每个元素是一个 numFeat 乘 window 的矩阵，这是 MATLAB 序列输入的标准格式，网络才能正确识别时间维度；第二，切分训练集和验证集必须在构造窗口之前完成，如果先构造窗口再随机划分样本，验证集里可能混入时间上与训练集交叠的样本，等于把未来信息泄漏给模型。
window 的取值会影响训练效果。取 24 是假设数据按小时采样，能覆盖一天的周期；取 48 可以看更长的跨日规律，但样本数量会减少。光伏数据夜里功率为 0，窗口跨到夜间时会引入大量零值片段，这部分对梯度贡献不大，但会让 LSTM 花更多轮次去拟合零值区间。源码里提到的剔除零值逻辑，在实际跑的时候可以结合数据分布决定是剔除还是保留。
提示：时间窗构造必须在切完训练集和验证集之后做，顺序反了，验证集指标会虚高，真实预测能力反而差。
4.2 网络结构与训练参数配置
EMD_KPCA_LSTM.m 是资源包的主脚本，网络搭建使用 MATLAB 深度学习工具箱。我在复现时会这样组织：
numFeatures = size(scoreTrain, 2); layers = [ sequenceInputLayer(numFeatures) lstmLayer(64, 'OutputMode', 'last') % 输出最后一个时刻的状态 dropoutLayer(0.2) % 抑制过拟合 fullyConnectedLayer(1) regressionLayer ]; options = trainingOptions('adam', ... 'InitialLearnRate', 1e-3, ... 'MiniBatchSize', 32, ... 'MaxEpochs', 150, ... 'Shuffle', 'never', ... 'ValidationData', {XVal, YVal}, ... 'ValidationFrequency', 10, ... 'Verbose', false); net = trainNetwork(XTrain, YTrain, layers, options);这里 lstmLayer 的 OutputMode 设成 last，因为我们做的是单步回归，只需要序列最后一个时刻的隐藏状态来预测窗口末端的功率。如果要做多步预测，可以改成 sequence 输出再接全连接层，但资源包的思路是单步滚动预测，我先按这个写。
训练参数里最容易被忽略的是 Shuffle。图像分类任务里随机打乱训练样本是常规操作，但时间序列不行。一旦打乱，模型看到的梯度顺序就失去了时间连续性，LSTM 学到的动态模式会被破坏。ValidationData 单独准备，不能从 XTrain 里抽。
| 参数 | 推荐取值 | 作用与风险 |
|---|---|---|
| InitialLearnRate | 1e-3 | 过大训练震荡，过小收敛极慢 |
| MiniBatchSize | 32 | 太大削弱序列动态性，太小梯度噪声大 |
| MaxEpochs | 150 | 结合早停策略，防止过拟合 |
| Shuffle | never | 时序数据必须保持顺序 |
| ValidationFrequency | 10 | 频次高能早发现发散，但拖慢训练 |
dropout 层设 0.2，适合光伏这种噪声偏大的数据。如果想要更强的表达能力，可以再加一层 lstmLayer，把第一个 LSTM 的 OutputMode 改成 sequence，再接第二个 LSTM。层数最多两层就够，样本量不大时三层以上基本是过拟合。
4.3 误差计算与结果判读
资源包的 calc_error.m 负责统计预测误差。我按通用做法实现 RMSE、MAE、MAPE 三个指标：
function err = calc_error(yTrue, yPred) err.RMSE = sqrt(mean((yTrue - yPred).^2)); err.MAE = mean(abs(yTrue - yPred)); err.MAPE = mean(abs((yTrue - yPred) ./ (yTrue + eps))) * 100; end注意 MAPE 的计算里加了 eps 防止零值除零。光伏功率夜间为 0，如果把 24 小时数据全拿去算 MAPE，深夜的零值会直接把指标拉爆。实际项目中我只统计太阳辐射度大于阈值的时段，比如辐射度大于 50 W/m² 的样本，这样算出来的 MAPE 才有横向可比性。这也是为什么源码里专门处理零值的原因。
预测曲线画出来后，重点看两个位置：清晨和傍晚，功率快速爬升和下降时预测是否跟得上；午间多云时段，功率剧烈抖动时误差是否集中。如果只有夜间拟合好，白天误差大，说明模型没学到辐射度和功率的真实映射关系，要回头检查 KPCA 的主成分是否保留了辐射度相关的高频 IMF。
5. 跑这套流程的五个坑：现象、原因与解决
5.1 IMF 数量不一致导致特征拼错位
现象：四个特征分别 EMD 分解后，imfCell 里列数不一样，强行拼矩阵时报维度错误，或者拼出来的特征序列在预测时出现周期性错位。
原因：每个特征的非平稳程度不同，分解出的 IMF 数量天然不一致，直接横向拼接默认了各列是同一分量，实际却是不同时间尺度混在一起。
解决：先统计 imfCounts，统一截断到最小值后再拼接。截断优于补零，补零列会在 RBF 核里产生虚假距离。对齐后把拼接矩阵存成 emd_data.mat，后面每次重跑都从这个中间结果开始。
5.2 测试集核矩阵中心化用错
现象：训练集误差一路下降，测试集预测曲线整体偏离真实值，误差比训练集大一个量级，看起来像 LSTM 过拟合，但调 dropout 和学习率都没用。
原因：预测阶段按测试集自己的核矩阵均值做了中心化，相当于换了坐标系。KPCA 主成分空间是基于训练集均值定义的，测试集用自身均值映射后，score 整体发生了平移。
解决：训练阶段把核矩阵中心化使用的 oneN 和均值信息保存到 KPCA_data.mat，预测阶段加载同一个 oneN，对测试核矩阵做同样的中心化操作。我习惯把整个 kPCA 模型封装成一个结构体，里面包含 gamma、特征向量、特征值和训练集均值，映射新数据时全部从结构体读取。
5.3 随机种子不固定，实验结果难以复现
现象：同样一段代码，连续跑三次，验证集 RMSE 差别能到百分之五到百分之十，让人分不清是模型差异还是随机波动。
原因：LSTM 权重初始化依赖随机数，GPU 上的浮点计算也有不确定性。没有固定种子，每次训练等于从不同起点出发。
解决：训练脚本开头加上rng(0)固定全局随机种子，并把训练日志和验证集指标一起保存。从那以后我每次做模型对比都会把种子固定在同一个值，否则两个网络结构之间的对比毫无意义。
5.4 洗牌导致的时间泄漏
现象：验证集 MAPE 低到反常，甚至比训练集还低，但真实场景下模型表现很差。
原因：构造完窗口样本后直接randperm打乱再切训练集和验证集，验证集样本的时间窗可能落在训练集样本的区间中间，模型在验证阶段偷看到了未来的数据。
解决：先按时间顺序把原始序列切分到训练集和验证集，再各自构造窗口样本，并且训练时 Shuffle 设成 never。这样的验证集才是真正意义上的未来时间段。
5.5 夜间 0 值主导损失
现象：MAPE 数值异常大，误差散点图在功率低值区堆了一大片。
原因：光伏功率夜间为 0，零值样本占总样本比例很高，回归损失被这些零值主导，LSTM 花大量容量去拟合零点，白天有效预测被挤占。
解决：有几种处理方式，我建议先计算日照时段占比，然后把训练目标限定在辐射度大于某个阈值的样本上，只统计这部分误差；或者像源码提示的那样，剔除零值样本再训练。注意剔除后样本的连续性会断裂，窗口构造要以剔除后的时间序列为准。
6. 把分解结果和核矩阵落盘缓存：调参效率的关键
6.1 先存中间结果，再调网络
EMD 分解和核矩阵构建是整个流程里计算量最大的两段，前者是迭代筛分，后者要算 n 乘 n 的距离矩阵，样本量上千之后开销很明显。但这两步的结果和网络参数无关，换 LSTM 层数、换学习率、换时间窗长度时，它们不需要重新计算。
我习惯把流程拆成两段跑。第一阶段只做分解和降维，结果存盘：
%% 第一次运行生成缓存，之后直接加载 if ~exist('emd_cache.mat', 'file') imfCell = decomposeAll(feat); % 四个特征逐列EMD save('emd_cache.mat', 'imfCell'); else load('emd_cache.mat', 'imfCell'); end if ~exist('kpca_cache.mat', 'file') K = exp(-gamma * distanceMatrix(featureBlock, featureBlock)); save('kpca_cache.mat', 'K', 'gamma'); else load('kpca_cache.mat', 'K', 'gamma'); end第二阶段专心调 LSTM，每次只从缓存里加载降维后的 score。在模拟项目 X 上，这个改动把每次重训的预处理时间从几分钟压到几秒，层数对比实验跑起来顺畅很多。资源包里本身带了 emd_data.mat 和 KPCA_data.mat 两个中间文件，说明原作者也是这个思路，中间结果落盘是标准操作。
这个技巧配合固定随机种子，整个 E-K-L 流程就变成：分解一次、降维一次、网络随便调。从那以后我每次做这类时间序列流程，都会强制把分解结果、核矩阵、随机种子存成一份带时间戳的 mat 文件，再进入调参循环，重跑时先看缓存有没有命中。这份习惯让我少走了很多重复计算的路，希望帮到你。
本文还有配套的精品资源，点击获取
