CEEMDAN-ISOS-VMD-GRU-ARIMA：非平稳时间序列预测全链路拆解

简介：这份资源面向计算机、电子信息工程、数学等专业的大学生及算法初学者，提供一套完整的CEEMDAN-ISOS-VMD-GRU-ARIMA时间序列预测实现方案，可用于课程设计、期末大作业或毕业设计。资源包共3个文件，包含2个CSV数据文件与1个Python源码文件，压缩包约52KB，数据文件用于模型训练与验证，源码文件则串联起信号分解、组合预测与结果评估的完整流程。代码基于Anaconda、PyCharm与TensorFlow环境编写，采用参数化编程思路，参数修改方便，且配有保姆级注释，几乎一行一注释，便于新手理解每一步的算法逻辑与数据流向。目前已有388人学习下载，读者可从中获得从数据预处理、分解策略到GRU与ARIMA组合建模的完整实现路径，并借助清晰注释快速上手调试与二次开发。

1. CEEMDAN-ISOS-VMD-GRU-ARIMA：一条把非平稳序列拆到骨头里的预测链路

风电功率、光伏出力、负荷曲线这类时间序列，直接丢进 GRU 或 ARIMA，MAPE 经常卡在 8% 到 15% 下不去。不是模型不行，是序列本身太"脏"——趋势、周期、突变、噪声全叠在一起，单一模型学不动。这套 CEEMDAN-ISOS-VMD-GRU-ARIMA 的思路就是先把序列按频率拆开，再用不同模型分别吃不同分量，最后重构。CEEMDAN 负责自适应分解出 IMF 分量，ISOS 做分量复杂度评估和重构，VMD 对高频分量二次分解，GRU 抓非线性残差，ARIMA 兜底线性趋势。适合做电力负荷、风速、光伏、销量这类强非平稳预测的从业者，也适合想把信号分解和深度学习串起来的新手。下面按"为什么这么拆 → 每步怎么落地 → 参数怎么调 → 哪里会翻车"的顺序讲透。

2. 分解链路的设计逻辑：为什么不是直接上 GRU

2.1 CEEMDAN 相比 EMD、EEMD 解决了什么

EMD 的模态混叠是老问题：一个 IMF 里混着不同时间尺度的成分，后续建模根本没法对应。EEMD 靠加高斯白噪声再做集合平均缓解，但噪声残留和重构误差又冒出来。CEEMDAN（Complete Ensemble EMD with Adaptive Noise）在每个分解阶段加入自适应白噪声，并且把每阶残差作为下一阶的输入，最终重构误差几乎为零。

落地时核心参数只有三个：

常见做法是用PyEMD库的 CEEMDAN 实现，或者自己按 Torres 2011 的论文写。我一般直接用 PyEMD，省得在筛选停止准则上反复调。

from PyEMD import CEEMDAN import numpy as np # signal: 一维 numpy 数组，长度建议 >= 1000 ceemdan = CEEMDAN(trials=200) # trials 对应 NR IMFs = ceemdan(signal) # 返回 shape=(n_imfs, len(signal)) residue = signal - IMFs.sum(axis=0) print("重构误差 max:", np.abs(residue).max())

trials就是噪声实现次数，200 是精度和耗时的平衡点。IMFs最后一行通常是趋势项，前面几行是高频到低频。重构误差如果大于 1e-10，说明 trials 太小或者信号里有 NaN，先查数据再做分解。

2.2 ISOS 做分量重构：把相似分量合并降维

CEEMDAN 出来十几个 IMF 是常态，逐个建模计算量爆炸，而且相邻 IMF 频率接近，分开建模反而引入冗余。ISOS（Improved Symbiotic Organisms Search）在这里的作用是：以样本熵或排列熵为相似度指标，把熵值接近的 IMF 聚成一组，组内叠加成一个新分量。

为什么用 ISOS 而不是直接 K-means？因为分量个数不确定，K-means 要预设 K。ISOS 是群智能优化，把"分组方案"编码成生物体位置，适应度函数用组内熵方差最小化，自动搜出最优分组数。

import antropy as ant import numpy as np def entropy_features(IMFs): # 对每个 IMF 算排列熵，作为 ISOS 的输入特征 return np.array([ant.perm_entropy(imf, normalize=True) for imf in IMFs]) def fitness(groups, entropies): # groups: list of list, 每个子列表是 IMF 索引 score = 0.0 for g in groups: if len(g) == 0: return 1e9 score += np.var(entropies[g]) return score

perm_entropy的order默认 3，序列短可以降到 2。适应度里组内方差越小说明同组分量越相似。ISOS 的种群规模和迭代次数一般设 30 和 100，再大收益递减。分组完成后，每组 IMF 直接相加得到重构分量，分量数通常从 12~15 降到 4~6。

2.3 VMD 二次分解高频分量：把突变拎出来

重构后频率最高的那个分量往往还混着突变和噪声，直接喂 GRU 会让模型去拟合噪声。VMD（Variational Mode Decomposition）把信号分解成 K 个窄带模态，对高频分量再拆一次，能把突变单独分离。

VMD 两个关键参数：K模态数和alpha带宽约束。K 太小欠分解，太大会过分解出虚假模态。常用做法是先设 K=3~5，看中心频率有没有重叠，重叠就减 K。

from vmdpy import VMD # 对高频重构分量 high_freq 做二次分解 alpha = 2000 # 带宽约束，越大带宽越窄 tau = 0 # 噪声容限，无先验噪声设 0 K = 4 # 模态数 DC = 0 # 不含直流分量 init = 1 # 中心频率初始化 tol = 1e-7 u, u_hat, omega = VMD(high_freq, alpha, tau, K, DC, init, tol) # u shape=(K, len), 每行是一个模态

alpha=2000对负荷类信号比较稳，风速信号可以降到 1000。分解完看omega最后一行的中心频率，如果两个模态频率差小于 0.01，说明 K 设大了。VMD 对端点效应敏感，序列首尾各截掉 5% 再分解，或者做镜像延拓。

2.4 GRU 与 ARIMA 的分工：线性和非线性各管一段

重构后的分量按熵值分两类：低熵分量规律性强，用 ARIMA 建模；高熵分量非线性强，用 GRU。这不是拍脑袋，是 ARIMA 对线性自相关结构建模效率高，GRU 对非线性映射强，混着用反而互相干扰。

GRU 输入用滑动窗口构造，窗口长度一般取 24（日周期）或 168（周周期）。ARIMA 的 (p,d,q) 用 AIC 定阶，d 用 ADF 检验确定。两部分预测结果按分量重构的逆运算叠加，得到最终预测。

提示：GRU 和 ARIMA 的预测误差量级要统一，建议都在归一化空间里建模，最后再反归一化，否则叠加时高幅值分量会淹没低幅值分量。

3. 从原始序列到预测结果：完整落地步骤

3.1 数据准备与平稳性检验

拿到序列先做三件事：缺失值插补、异常值剔除、平稳性检验。缺失值用线性插值，异常值用 3σ 或 IQR，别用均值填充，会抹掉突变信息。

import pandas as pd import numpy as np from statsmodels.tsa.stattools import adfuller df = pd.read_csv("load.csv", parse_dates=["timestamp"]) series = df["value"].interpolate(method="linear").values # 3σ 异常值替换为插值 mu, sigma = series.mean(), series.std() outliers = np.abs(series - mu) > 3 * sigma series[outliers] = np.nan series = pd.Series(series).interpolate().values # ADF 检验 adf_stat, p_value, _, _, crit, _ = adfuller(series, autolag="AIC") print(f"ADF={adf_stat:.4f}, p={p_value:.4f}") if p_value > 0.05: print("非平稳，ARIMA 需要差分")

autolag="AIC"自动选滞后阶。p 值大于 0.05 说明非平稳，ARIMA 的 d 至少取 1。做完差分再检验一次，直到平稳。注意 CEEMDAN 分解前不需要序列平稳，分解本身就能处理非平稳，但 ARIMA 分支必须平稳。

3.2 CEEMDAN 分解与分量导出

from PyEMD import CEEMDAN ceemdan = CEEMDAN(trials=200) IMFs = ceemdan(series) np.save("imfs.npy", IMFs) # 检查分量数和重构误差 print("IMF 数量:", IMFs.shape[0]) recon = IMFs.sum(axis=0) print("最大重构误差:", np.abs(recon - series).max())

trials=200对长度 5000 以内的序列够用。如果分量数超过 15，先别急着调参，检查序列里有没有阶跃或断点，CEEMDAN 对突变会过度分解。导出后每个 IMF 单独存一列，方便后续按索引分组。

3.3 ISOS 分组重构

import antropy as ant import numpy as np IMFs = np.load("imfs.npy") entropies = np.array([ant.perm_entropy(imf, normalize=True) for imf in IMFs]) # ISOS 简化实现：种群初始化 + 共生阶段 def isos_group(entropies, pop_size=30, max_iter=100): n = len(entropies) # 每个个体是一个分组边界向量，排序后切分 pop = np.random.rand(pop_size, n) best_pos, best_fit = None, 1e9 for _ in range(max_iter): for i in range(pop_size): order = np.argsort(pop[i]) # 按熵值排序后均分 4 组 groups = np.array_split(order, 4) fit = sum(np.var(entropies[g]) for g in groups) if fit < best_fit: best_fit, best_pos = fit, groups # 共生更新（简化：向最优靠拢） for i in range(pop_size): pop[i] += 0.1 * (pop[np.argmin([sum(np.var(entropies[g]) for g in np.array_split(np.argsort(p), 4)) for p in pop])] - pop[i]) return best_pos groups = isos_group(entropies) reconstructed = np.array([IMFs[g].sum(axis=0) for g in groups]) np.save("reconstructed.npy", reconstructed) print("重构后分量数:", reconstructed.shape[0])

这里把分组问题简化成排序切分，实际 ISOS 要做互利、偏利、寄生三个阶段。pop_size=30、max_iter=100对 15 个 IMF 足够。分组数不固定时，可以在适应度里加惩罚项控制组数。重构后分量按熵值从低到高排列，低熵走 ARIMA，高熵走 GRU。

3.4 VMD 二次分解高频分量

from vmdpy import VMD import numpy as np reconstructed = np.load("reconstructed.npy") high_freq = reconstructed[-1] # 熵值最高的分量 u, u_hat, omega = VMD(high_freq, alpha=2000, tau=0, K=4, DC=0, init=1, tol=1e-7) np.save("vmd_modes.npy", u) print("中心频率:", omega[-1])

K=4是起点，看omega[-1]里相邻频率差，小于 0.01 就减 K。VMD 分解出的模态按频率从低到高排列，低频模态可以并回 ARIMA 分支，高频模态进 GRU。端点效应处理：分解前对high_freq做镜像延拓，分解后截掉延拓部分。

3.5 GRU 建模与训练

import torch import torch.nn as nn import numpy as np class GRUNet(nn.Module): def __init__(self, input_size=1, hidden_size=64, num_layers=2, output_size=1): super().__init__() self.gru = nn.GRU(input_size, hidden_size, num_layers, batch_first=True) self.fc = nn.Linear(hidden_size, output_size) def forward(self, x): out, _ = self.gru(x) return self.fc(out[:, -1, :]) def make_windows(data, window=24): X, y = [], [] for i in range(len(data) - window): X.append(data[i:i+window]) y.append(data[i+window]) return np.array(X), np.array(y) # 对每个高频模态分别建模 modes = np.load("vmd_modes.npy") preds = [] for mode in modes: X, y = make_windows(mode, window=24) X = torch.tensor(X, dtype=torch.float32).unsqueeze(-1) y = torch.tensor(y, dtype=torch.float32).unsqueeze(-1) model = GRUNet() opt = torch.optim.Adam(model.parameters(), lr=1e-3) loss_fn = nn.MSELoss() for epoch in range(200): model.train() opt.zero_grad() loss = loss_fn(model(X), y) loss.backward() opt.step() model.eval() with torch.no_grad(): pred = model(X).squeeze().numpy() preds.append(pred)

hidden_size=64、num_layers=2对大多数分量够用，分量长可以加到 128。window=24对应日周期，周周期数据改 168。学习率 1e-3 配 Adam，200 epoch 后看 loss 曲线，震荡就降到 5e-4。每个模态单独训练，别把所有模态堆一起，否则模型会偏向高幅值模态。

3.6 ARIMA 建模与叠加

from statsmodels.tsa.arima.model import ARIMA import numpy as np reconstructed = np.load("reconstructed.npy") low_freq = reconstructed[:-1] # 低熵分量走 ARIMA arima_preds = [] for comp in low_freq: model = ARIMA(comp, order=(2, 1, 2)) fit = model.fit() pred = fit.forecast(steps=len(comp) - 24) arima_preds.append(pred) # 叠加：ARIMA 预测 + GRU 预测 + VMD 模态预测 final = np.sum(arima_preds, axis=0) + np.sum(preds, axis=0)

order=(2,1,2)是起点，用 AIC 网格搜 (0~3, d, 0~3)。forecast的 steps 要和 GRU 预测长度对齐。叠加时注意各分支预测起点一致，差一个时间步结果就全错。反归一化如果做了归一化，在这一步之后做。

4. 避坑与排查：这套链路最容易翻车的 5 个地方

4.1 CEEMDAN 分量数爆炸，跑到 20 个以上

现象：分解出来 20 多个 IMF，后面几个幅值接近零。原因：序列里有阶跃或异常值，CEEMDAN 把突变当成独立模态反复拆。解决：分解前做异常值检测，阶跃点做平滑或分段处理；trials从 200 降到 100，噪声标准差Nstd从 0.2 降到 0.1。

4.2 ISOS 分组结果每次跑都不一样

现象：同一份数据跑两次，分组数差 2 个。原因：ISOS 是随机初始化种群，没设随机种子。解决：np.random.seed(42)固定种子；种群规模从 30 加到 50，迭代从 100 加到 200，让搜索更稳定。如果还抖，改用熵值阈值直接切分，牺牲一点最优性换可复现。

4.3 VMD 分解出虚假模态，中心频率重叠

现象：omega[-1]里两个模态频率差小于 0.005。原因：K 设大了，或者alpha太小导致带宽过宽。解决：K 从 4 降到 3，alpha从 2000 加到 3000。判断标准：分解后每个模态的频谱应该只有一个主峰，多峰就是过分解。

4.4 GRU 预测滞后一个时间步

现象：预测曲线整体右移，跟真实值错位。原因：滑动窗口构造时y的索引偏了，或者归一化用了全局均值导致信息泄漏。解决：检查make_windows里y.append(data[i+window])是否正确；归一化用训练集均值方差，别用全序列。滞后一步在时序预测里是血泪经验，MAPE 能差 3 个点。

4.5 ARIMA 和 GRU 叠加后误差反而变大

现象：单分支 MAPE 8%，叠加后 12%。原因：两个分支预测误差量级不匹配，或者叠加时没对齐时间索引。解决：所有分支在归一化空间预测，叠加后再反归一化；检查forecast的 steps 和 GRU 预测长度是否一致；如果还差，给每个分支按验证集误差加权，别直接等权相加。

5. 进阶技巧：用滚动预测和误差加权把 MAPE 再压 2 个点

这套链路跑通后，单次预测 MAPE 大概在 5%~8%。想再往下压，两个技巧最实用。

第一个是滚动预测。别一次性预测未来 24 步，改成每预测一步就把真实值喂回去，重新构造窗口。GRU 分支用model(X)单步推理，ARIMA 用fit.append()更新。代价是计算量翻 24 倍，但 MAPE 能降 1~2 个点。代码上把make_windows改成在线版本：

def rolling_forecast(model, init_window, steps): window = list(init_window) preds = [] for _ in range(steps): x = torch.tensor(window[-24:], dtype=torch.float32).view(1, 24, 1) with torch.no_grad(): p = model(x).item() preds.append(p) window.append(p) # 用预测值滚动，有真实值就替换成真实值 return np.array(preds)

有真实值时把window.append(p)换成window.append(true_value)，误差不会累积。

第二个是分支误差加权。验证集上分别算 ARIMA 和 GRU 的 MAPE，权重取倒数归一化：

叠加时final = w_arima * arima_pred + w_gru * gru_pred。权重每月用最近数据重算一次，别固定死。我一般还会留一个 10% 的验证集做早停，GRU 训练时loss连续 20 epoch 不降就停，省得过拟合。

最后一个习惯：每次改参数前先存一版预测结果和 MAPE，改完对比。这套链路参数多，不记录的话跑两周就忘了哪版最好。希望帮到你。

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