简介这是一套基于Python实现CEEMDAN-ISOS-VMD-GRU-ARIMA组合时间序列预测的完整工程面向需要完成课程设计或毕业设计的高校学生以及想掌握多阶段分解预测模型的研究者。压缩包内包含一个主程序文件和两个CSV格式样例数据集资源整体只有52KB环境配置好后可直接运行省去整理数据的麻烦。代码在AnacondaPyCharmTensorFlow环境下编写每一步都配有保姆级注释从CEEMDAN自适应分解到ISOS-VMD降噪再到GRU与ARIMA的融合预测逻辑链路完整适合基础薄弱的学习者逐行研读。程序采用参数化设计分解层数、迭代次数、模型维度等超参数均可自由修改并能打印中间结果方便对比不同参数下的预测效果。目前已有387人参与学习下载适用于计算机、电子信息、数学等相关专业的实训作业、期末大作业及毕业设计也可作为进一步开展组合优化预测研究的参考起点。1. 时间序列预测不是换模型是先拆干净CEEMDAN-ISOS-VMD-GRU-ARIMA 解决什么问题先看一个几乎所有做时间序列预测的人都会遇到的场景拿一份电力负荷或者风速数据训练集上 GRU 拟合得很好一到测试段就整体滞后一拍换成 ARIMA又对突变段毫无反应。问题不在模型不够新在输入序列本身——非平稳、多尺度、带噪声单一模型根本吃不下。CEEMDAN-ISOS-VMD-GRU-ARIMA 这套长名字要解决的就是拆解和分工CEEMDAN 先做一轮自适应分解ISOS 再帮 VMD 找到最优的分解参数最后高频分量交给 GRU、低频趋势交给 ARIMA各自预测再叠回去。这套方案适合两类人。一类是已经在用 Python 做数据分析与可视化想把手头风电功率、负荷或量化交易策略里的价格序列预测精度再往上提一截的分析师另一类是刚做完 Python 入门想找一个能落地的完整预测工程方向的开发者。它不算新但把五段技术串成一条可复现流水线细节全在参数和边界处理上。下面按我自己的实现路径展开代码可以直接照着改。2. 算法链路与选型逻辑为什么分解要套两轮预测要 GRU 配 ARIMA2.1 CEEMDAN 打底解决 EEMD 的残留噪声与模态混叠先理解为什么要做第一轮分解。原始序列里同时混着长期趋势、周期波动和随机噪声GRU 这类模型虽然能拟合非线性但面对这种多尺度叠加的信号容易把噪声也学进参数里。经验模态分解EMD能把序列按频率从高到低拆成若干本征模态函数但 EMD 有模态混叠问题同一个模态里可能同时出现不同频率的成分。EEMD 通过加白噪声来解决分解结果却会残留一部分噪声重构误差不干净。CEEMDAN 的做法是在每个分解阶段加入自适应白噪声并且在每一层只求一次残余分解完备性比 EEMD 好重构误差也小得多。在 PyEMD 库里调用它不需要自己写迭代逻辑通常我只需要关注三个参数集成次数 trials、噪声幅值 epsilon 和极值检测方式。trials 越大分解越稳定但越慢一般在 50 到 200 之间epsilon 太大会让分量失真太小则退化回 EMD。做 CEEMDAN 之前还有一个前提序列最好已经标准化否则幅值差异会把分解带偏。2.2 ISOS 给 VMD 找参数K 和 alpha 为什么不能拍脑袋CEEMDAN 拆完之后某些分量仍然很复杂比如第一个高频 IMF 可能还是包含多个频带或者残差项里还有一段明显的中长期波动。这时候就要上第二层分解 VMD。VMD 和 EMD 思路不同它把信号分解问题放到变分框架里求解能控制每个模态的中心频率和带宽但代价是必须预设两个关键参数模态数 K 和惩罚因子 alpha。K 选小了欠分解两个频率成分挤在一个模态里K 选大了模态被拆碎出现假分量。alpha 控制模态带宽alpha 太大每个模态都被压得很窄容易把真实成分切成好几段alpha 太小模态之间互相重叠。手动试 K 和 alpha 不是不行但每换一组数据就得重新试网格搜索又要跑几十次 VMD一次 VMD 在长序列上可能要几十秒非常熬人。ISOS 我按改进的麻雀搜索算法来落地把 K 和 alpha 当作二维变量去搜索目标函数用最小包络熵。麻雀搜索收敛快、实现简单种群和迭代次数都不用设很大跑一轮也就十几分钟比网格搜索省事得多。2.3 GRU 和 ARIMA 的分工逻辑高频给深度网络趋势给统计模型分解做完预测器也要分开。CEEMDAN 拆出来的第一个 IMF 通常是高频分量波动密集、非线性强ARIMA 这类线性模型拟合这种分量几乎必炸GRU 却擅长从这种短时波动里抓依赖关系它比 LSTM 少一个门参数更少在单个分量这种小样本上反而不容易过拟合训练也快。低频趋势项和残差项则平稳得多用 ARIMA 更合适可解释性强而且对轻微趋势用一阶差分就能解决。所以整套链路的逻辑是CEEMDAN 负责把序列从复杂变简单ISOS 负责解决 VMD 的调参问题VMD 负责把剩下的复杂分量拆到可用GRU 和 ARIMA 按频率分工各自预测最后把预测值相加再做反归一化回到原始尺度。每一步解决的都是上一环节暴露出来的具体问题没有哪一环是凑数。3. 环境、数据与 CEEMDAN 第一次分解依赖安装和最小实现3.1 Python 环境准备与三个核心库的选型开始之前先把依赖装齐。我常用的 Python 版本是 3.9 到 3.11装包命令固定这几行pip install numpy pandas scikit-learn matplotlib pip install tensorflow statsmodels pip install PyEMD vmdpy如果网络条件一般可以在 pip 后面加-i指定国内镜像源常用的清华、阿里源都行。PyEMD 提供 CEEMDANvmdpy 提供 VMDstatsmodels 提供 ARIMAtensorflow 提供 GRU。注意 PyEMD 的 CEEMDAN 接口在不同小版本里略有差异装完之后先打印一下函数签名再往下写不然一运行就报参数错很浪费感情。3.2 数据读取、标准化与训练集划分数据格式我统一用两列的 CSV一列时间一列数值。读取和标准化的代码要特别小心一个坑标准化器的 fit 只能用训练段不能用全序列。否则测试段的极值会被提前看见这就是信息泄漏离线指标会虚高。import numpy as np import pandas as pd from sklearn.preprocessing import MinMaxScaler # 读入数据ds 是时间列y 是待预测序列 df pd.read_csv(load.csv, parse_dates[ds], index_colds) raw df[y].values.astype(float) # 缺失值先线性插值不要直接删 raw pd.Series(raw).interpolate(limit_directionboth).values # 按时间顺序切分前 80% 训练后 20% 测试 n_train int(len(raw) * 0.8) train_raw, test_raw raw[:n_train], raw[n_train:] # 标准化只用训练段统计量 scaler MinMaxScaler(feature_range(0, 1)) scaler.fit(train_raw.reshape(-1, 1)) scaled scaler.transform(raw.reshape(-1, 1)).ravel() print(scaled length:, len(scaled), train length:, n_train)这里有几个参数值得解释。MinMaxScaler 的 feature_range 我固定用 0 到 1GRU 的激活函数对输入尺度敏感0 到 1 比 -1 到 1 在 MSE 损失下更容易收敛。缺失值处理用插值而不是填充因为填充会引入一段平坦的伪信号分解时会产生虚假模态。train_raw 只用于 fitscaled 是对全序列的变换测试段数值落在训练段范围内这是可接受的部署口径。3.3 CEEMDAN 分解的最小 Python 代码与分量初筛标准化做完就可以做第一次分解。CEEMDAN 在 PyEMD 里的调用非常简单难点全在参数和后续分量筛选from PyEMD import CEEMDAN # trials 是集成次数epsilon 是噪声幅值parabol 是抛物线样条插值 ceemdan CEEMDAN(trials100, epsilon0.005, extrema_detectionparabol) imfs ceemdan(scaled) print(IMF shape:, imfs.shape)返回的 imfs 是二维数组行数是分量个数列数等于序列长度。最后一行是残差项一般是一个单调趋势或一个幅度很小的平稳项。参数上我解释一下trials 设 100 是稳定性和耗时的折中序列超过一万点可以降到 50epsilon 设 0.005 是我大部分数据上的起点如果分解出的第一个 IMF 振幅明显偏大或出现上下包络不对称就把 epsilon 调到 0.001 到 0.01 之间再试。extrema_detection 用 parabol 是因为它在连续信号上的包络拟合比默认的 simple 更平滑模态混叠更少。分解完不要急着训练先画一张分量图肉眼检查有没有模态混叠import matplotlib.pyplot as plt plt.figure(figsize(12, 2 * imfs.shape[0])) for i in range(imfs.shape[0]): plt.subplot(imfs.shape[0], 1, i 1) plt.plot(imfs[i], linewidth0.8) plt.ylabel(fIMF{i}) plt.savefig(ceemdan_imfs.png, dpi150)这一步能看出很多问题如果相邻两个 IMF 频率几乎一样说明 CEEMDAN 没拆干净如果某个 IMF 在某个时间段突然幅度归零可能是边界效应。分量图存下来之后下一个环节才决定要对哪些分量做 VMD。4. ISOS 优化 VMD 与 GRU-ARIMA 分量预测核心代码与参数4.1 包络熵与 ISOS 搜索 K、alpha 的 Python 实现VMD 的适应度函数我常用最小包络熵。包络熵的含义是信号经希尔伯特变换得到包络包络归一化后求香农熵。模态混叠少、频带干净时包络形态简单熵值低混叠严重时包络起伏杂乱熵值高。代码很短from scipy.signal import hilbert def envelope_entropy(imf): env np.abs(hilbert(imf)) p env / (np.sum(env) 1e-12) return -np.sum(p * np.log(p 1e-12))加 1e-12 只是防止 log 零出 NaN不影响结果。有了适应度函数再把 VMD 封装一下from vmdpy import VMD def vmd_envelope_entropy(series, K, alpha, tau0): # VMD 参数: f, alpha, tau, K, DC, init, tol u, _, _ VMD(series, alpha, tau, int(K), 0, 1, 1e-7) entropies [envelope_entropy(u[i]) for i in range(u.shape[0])] return np.mean(entropies)ISOS 我按改进麻雀搜索实现麻雀搜索里发现者负责全局探索、加入者跟随发现者、警戒者负责局部扰动。完整论文里的位置更新公式比较长实际工程我会做简化保留精英个体发现者按指数递减步长搜索加入者向最优位置靠拢并带随机扰动警戒者做小幅度随机游走。这样实现短、收敛也够用。import numpy as np class ISOS: def __init__(self, objective, dim2, lbNone, ubNone, pop12, max_iter15, pd_ratio0.2, sd_ratio0.1): self.objective objective self.dim dim self.lb np.array(lb) self.ub np.array(ub) self.pop pop self.max_iter max_iter self.pd_num max(2, int(pop * pd_ratio)) self.sd_num max(1, int(pop * sd_ratio)) self.X np.random.uniform(lb, ub, size(pop, dim)) self.f np.array([objective(x) for x in self.X]) best_idx int(np.argmin(self.f)) self.best_x self.X[best_idx].copy() self.best_f self.f[best_idx] def clip(self, x): return np.clip(x, self.lb, self.ub) def fit(self): for _ in range(self.max_iter): idx np.argsort(self.f) X self.X[idx] worst X[-1].copy() R2 np.random.rand() # 发现者 for i in range(self.pd_num): if R2 0.8: new X[i] * np.exp(-i / (np.random.rand() * self.max_iter 1e-6)) else: new X[i] np.random.randn(self.dim) * 0.1 X[i] self.clip(new) # 加入者 for i in range(self.pd_num, self.pop): if i self.pop / 2: new worst np.random.randn(self.dim) * 0.1 else: A (np.random.rand(self.dim) * 2 - 1) * 0.1 new X[i] np.abs(X[i] - X[0]) A X[i] self.clip(new) # 警戒者 for _ in range(self.sd_num): j np.random.randint(self.pop) X[j] self.clip(X[j] 0.01 * np.random.randn(self.dim)) self.X X self.f np.array([self.objective(x) for x in self.X]) cur int(np.argmin(self.f)) if self.f[cur] self.best_f: self.best_f self.f[cur] self.best_x self.X[cur].copy() return self.best_x, self.best_f调用时把 K 的搜索范围设为 3 到 10alpha 设为 100 到 3000。K 必须取整数所以在目标函数里 round 一下lb [3, 100] ub [10, 3000] def obj(x): return vmd_envelope_entropy(target_imf, round(x[0]), x[1]) optimizer ISOS(obj, lblb, ubub, pop10, max_iter15) best_k, best_alpha optimizer.fit() print(fbest K{best_k:.0f}, alpha{best_alpha:.1f}, entropy{optimizer.best_f:.4f})两个值得注意的点。第一目标函数里的 VMD 是在 CEEMDAN 某个分量上做的不要在整条原始序列上直接搜那样搜出来的参数会被其他频带干扰。第二ISOS 每次运行有随机性最好连续跑两次取熵值最小的那一组避免某次初始化太差陷入局部最优。4.2 用最优参数对 CEEMDAN 分量做 VMD 二次分解参数搜出来之后对选中的分量做最终分解u, _, _ VMD(target_imf, best_alpha, 0, round(best_k), 0, 1, 1e-7) print(VMD sub-modals shape:, u.shape)这里我一般只对两个地方做 VMD一个是 CEEMDAN 的第一个高频 IMF如果它的包络熵依然明显高于其他 IMF另一个是残差项里如果存在明显的中长期波动。原因很简单所有分量都套 VMD 会导致分量总数膨胀到十几个甚至二十几个每个分量都要单独训练一个模型计算量成倍增加预测误差却不会跟着降。筛选哪些 CEEMDAN 分量值得保留我按方差贡献率来var_ratio np.array([np.var(imfs[i]) for i in range(imfs.shape[0])]) var_ratio var_ratio / var_ratio.sum() keep var_ratio 0.01 print(keep components:, np.where(keep)[0])这个 0.01 的阈值是经验值数据噪声占比高时可以提高到 0.02但要保证保留分量的累计方差贡献在 95% 以上。被过滤掉的低能量分量通常属于纯噪声预测它对最终结果贡献极低反而会把误差引入叠加结果。4.3 GRU 建模高频分量滑窗构造与训练参数拿到分量清单后高频分量统一交给 GRU。先构造滑窗这一步直接决定模型能不能学到跨步依赖def make_windows(data, n_steps): X, y [], [] for i in range(n_steps, len(data)): X.append(data[i - n_steps:i]) y.append(data[i]) return np.array(X).reshape(-1, n_steps, 1), np.array(y) n_steps 24 X_train, y_train make_windows(train_part, n_steps) X_test, y_test make_windows(test_part, n_steps)滑窗长度的选择如果数据有明显周期n_steps 至少取两个周期长度没有明显周期就用 24 到 48 起步。窗口太短GRU 只能看到局部抖动预测会变成滞后一拍的复制窗口太长训练样本变少高频分量本来数据量就不大过拟合风险上升。GRU 网络我保持最小可用结构from tensorflow.keras.models import Sequential from tensorflow.keras.layers import GRU, Dense, Dropout def build_gru(n_steps, units64): model Sequential() model.add(GRU(units, activationtanh, input_shape(n_steps, 1))) model.add(Dropout(0.2)) model.add(Dense(1)) model.compile(optimizeradam, lossmse, metrics[mae]) return model model build_gru(n_steps) model.fit(X_train, y_train, epochs60, batch_size32, validation_split0.1, verbose0)units 设 64 是我在几千点数据量下的默认值分量样本少就降到 32样本多可以试 128。Dropout 放 0.2 防止高频噪声被死记。训练时 watch 验证集 loss如果验证 loss 前 10 轮不降就降学习率或加 epochs不要无脑堆层数。GRU 在分量预测上用单层就够堆两层在这个任务里收益很微弱。测试段预测用滚动方式把上一个预测值当作下一步输入这和实际部署一致def rolling_predict(model, init_window, n_forecast, n_steps): preds [] current init_window.copy() for _ in range(n_forecast): p model.predict(current.reshape(1, n_steps, 1), verbose0)[0, 0] preds.append(p) current np.roll(current, -1) current[-1] p return np.array(preds)注意 init_window 要取自训练段末尾的 n_steps 个真实值不能取测试段开头的值否则就是拿答案当输入。4.4 ARIMA 拟合低频趋势与整体叠加重构低频分量和残差项交给 ARIMA。建模前先做 ADF 平稳性检验非平稳就先差一阶from statsmodels.tsa.stattools import adfuller from statsmodels.tsa.arima.model import ARIMA adf_stat adfuller(train_part) print(ADF p-value:, adf_stat[1])p 值大于 0.05 就在 ARIMA 里设 d1否则 d0。选 p 和 q 我一般用 AIC 小范围搜索阶数控制在 5 以内import itertools best_aic, best_order float(inf), None for p, q in itertools.product(range(3), range(3)): try: res ARIMA(train_part, order(p, 1, q)).fit() except Exception: continue if res.aic best_aic: best_aic, best_order res.aic, (p, 1, q) print(best ARIMA order:, best_order)然后用最优阶数做测试段预测。测试段不长就直接 forecast测试段超过 50 步建议滚动 re-fit虽然慢但精度高def arima_forecast(train_part, n_forecast, order): res ARIMA(train_part, orderorder).fit() return res.forecast(stepsn_forecast)最后一环是叠加。把所有分量的预测值按时间点相加得到归一化尺度的预测序列再走一遍 inverse_transform 回到原始单位pred_scaled np.zeros(n_forecast) for pred in component_preds: pred_scaled pred pred_raw scaler.inverse_transform(pred_scaled.reshape(-1, 1)).ravel()到这里整套预测流水线就通了。注意 component_preds 里的每条预测序列长度必须一致如果某个模型输出短一截叠加时 numpy 会广播出错报错信息还不直观。我习惯在叠加前断言所有长度相等省得排半天。5. 避坑记录与常见问题排查五个会让整套模型翻车的细节5.1 现象离线指标很好一上线就崩用全序列做了 MinMaxScaler 的 fit或者在整段序列上做了 CEEMDAN 和 VMD再切训练测试段测试段的统计信息已经被分解器看到了。CEEMDAN 和 VMD 都是全局方法整段分解时边界处的模态会受到未来数据影响这属于信息泄漏。解决方法是标准化器只 fit 训练段分解尽量在训练段上做如果必须对全序列分解做分步预测就要接受这种边界效应并在最后评估时说明误差口径。更稳妥的做法是滚动分解每次只对训练段加最近一段历史窗口做预测窗口长度固定。5.2 现象VMD 拆出的子模态长得几乎一样K 和 alpha 没配对好。K 太大把真实成分切碎alpha 太大又把每个模态压得太窄结果几个子模态共享同一段频带。判断方法是画 VMD 子模态的频谱中心频率重叠就是混叠。解决方法是把适应度从纯包络熵改成包络熵加模态中心频率间隔惩罚或者缩小 ISOS 的搜索范围。我实际使用时会增加一个约束任意两个模态中心频率的差不能小于某个阈值否则给目标函数加一个大的惩罚项逼 ISOS 远离这种参数组合。5.3 现象GRU 训练 loss 不降预测整体滞后一拍三个常见原因输入没标准化、滑窗太短、学习率偏高。标准化只做训练段 fit 后再 transform滑窗小于两个周期会让模型只学会上一时刻的复制。还有一个隐蔽问题高频分量里如果还有残余噪声GRU 会去拟合噪声训练 loss 降到一定程度就不再下降。解决方法是把该分量重新做一次 VMD 或者对分量做轻平滑而不是改网络结构。验证滞后最简单的办法是画预测值和真实值的对比图滞后一拍几乎都是同一时间套准地偏差一位。5.4 现象ARIMA 对非平稳分量硬建模预测值直接发散直接把 CEEMDAN 的第一个 IMF 喂给 ARIMAADF 检验没做差分数也没设ARIMA 在非平稳序列上预测几步之后就开始指数式发散。解决方法是只对低频和残差分量用 ARIMA建模前必须跑 adfullerp 值大于 0.05 就 d 至少为 1。还有一种情况是残差项里还带着微弱的周期成分ARIMA 的 p 和 q 阶数压不住这时可以换 SARIMAX 加季节项。统计模型的容错率低参数要按检验结果定不能一套 (2,1,2) 打天下。5.5 现象每个分量都预测得不错叠加之后反而变差分量的个别预测误差在叠加时没有抵消而是同向累加。高频分量的 GRU 在突变点附近容易整体低估几个低频分量又把趋势略微高估叠加后系统偏差被放大。解决方法是先看每个分量的独立误差找出偏差方向一致的分量组对它们做加权而不是简单求和或者加一层融合模型把所有分量预测值作为特征用线性回归学一个最佳权重。这个融合层很轻量但能把消融对比里的分数救回来不少。6. 验证这套模型真的变强了消融对比与误差口径6.1 消融对比表想说服别人或者说服自己这套长链条值得上最直接的办法是做一组消融实验。同一份数据、同一个切分和标准化口径依次跑下面几挡方案RMSEMAEMAPE只用 GRU用你的数据跑用你的数据跑用你的数据跑CEEMDAN GRUCEEMDAN VMD GRUCEEMDAN ISOS-VMD GRUCEEMDAN ISOS-VMD GRU ARIMA每一挡之间只加一个环节这样才能看出 ISOS 搜出来的参数和 ARIMA 接替低频分量到底贡献了多少。三组指标都在原始尺度上算不要在归一化尺度里比较因为 inverse_transform 之后误差量级会被放大但这也是真实部署看到的口径。6.2 验证时的两个习惯第一个习惯是固定误差口径。测试段预测必须用滚动方式上一步预测值作为下一步输入不允许用真实值喂给 GRU 和 ARIMA。很多离线指标虚高就是因为预测时偷偷用了真实滞后值。第二个习惯是每个分量单独画预测图。我现在每次跑完这套流水线都会先看一眼高频 IMF 的 GRU 预测有没有滞后、残差项的 ARIMA 预测有没有发散确认每个环节都合理之后再去看叠加结果。这样做的原因是叠加后的总误差会掩盖分量的内部问题而内部问题才是下一次迭代真正要改的地方。这套流程跑顺之后最大的体会是模型的创新不在网络结构而在把信号拆到什么程度、每个分量用什么模型去接。先看分解再谈预测最后用消融对比守住底线方向就不会跑偏。希望帮到你。本文还有配套的精品资源点击获取