今年年初有个做设备状态监测的朋友拿了一组振动数据来找我抱怨说ARMA预测结果完全没法看问我能不能救一下。我扫了一眼那组曲线明显的趋势项叠着两三个周期的波动再加一堆毛刺噪声——典型的不平稳序列。这种情况下直接上ARMA说难听点就是让模型去猜一个它根本没学过、也不符合假设的东西。后来我把序列先用小波分解拆开提取出一层低频趋势分量和好几层细节分量对每一层分别做ARMA建模预测最后把各分量预测结果相加。这么一折腾预测精度和稳定性都上来了。这篇文章就把这套基于小波分解和ARMA预测的流程完整写出来包括能直接复制的Matlab代码、每一处参数选择的理由以及我在实际项目里踩过的坑和调参心得。如果你也在处理非平稳、含噪声的时间序列手里有Matlab但不清楚怎么把这些工具串起来这篇应该能帮到你。1. 这套组合拳到底解决了什么问题1.1 ARMA的局限平稳性是这个模型的命门先说清楚为什么单独用ARMA会翻车。ARMA自回归滑动平均模型的核心假设是序列来自一个平稳随机过程也就是说序列的均值、方差、自相关结构都不随时间的推移而改变。平稳性不是教科书上拿来考试的抽象概念它直接决定了ARMA的数学推导能不能成立。如果序列里有趋势项均值就在变化如果有周期性波动自相关结构就不稳定如果噪声方差忽大忽小ARMA估计出来的参数置信区间基本没什么意义。大多数人拿到数据后的第一反应是直接拟合ARMA试试结果就是拟合阶段看着还行一到预测就飞掉。为什么因为ARMA本质上是利用序列的线性相关结构做外推它只能捕捉那种相对稳定的短期依赖。你给它一个叠了趋势、周期、噪声的复杂序列它会把趋势当成均值的一部分、把噪声当成有用信号来学。预测步数一拉长误差迅速扩散。有人会说那先用差分去掉趋势、再对差分序列做ARMA不就行了差分确实能处理一部分非平稳趋势但周期性成分它不是线性趋势差分一次根本消不掉多差几次又容易把有用信息差掉。还有人会用STL或X11做季节性分解这方法在规则周期场景里挺好用但面对多个周期叠加、周期还不固定的数据就力不从心。1.2 小波分解的本质把信号拆成一堆慢件和快件小波分解的核心思想可以类比成拆表一块机械手表有秒针、分针、时针它们的运动速度差异很大但共同驱动着表盘。小波分解就是按频率快慢把原始信号拆成几个正交分量低频部分保留趋势和大尺度波动高频部分保留周期细节和噪声。关键在于这些分量之间没有信息冗余而且满足完美重构条件——把各分量相加能完全还原原始信号。用数学语言说小波分解是把信号投影到一组由小波基函数伸缩和平移构成的子空间里。不同尺度对应不同频带因此天然适合处理那种低频趋势多周期叠加高频噪声的混合信号。分解之后每一层分量都相对平稳、频带更窄ARMA模型在这样干净得多的信号上拟合捕捉到的才是真正有意义的短期相关结构。这就是先分解、后预测、再重构框架能提高预测精度的根本原因你不是在同一个模型里强行处理一堆性质完全不同的成分而是把复杂问题拆成一堆简单问题各个击破。1.3 适用范围什么数据值得用小波预处理小波分解不是万能的它有明确的适用边界。我自己的经验是至少满足以下两个条件才值得用这套方案序列长度要够最好有几百个点以上。分解层数太少效果有限太多又会有边界效应和数据缩短的问题。长度不够的时候小波分解的价值会被边缘误差吞掉。信号有明显的多频带结构也就是趋势项、周期项、噪声项叠加得明显。如果你手里的序列基本就是白噪声或者本身已经非常平稳那直接ARMA就够加小波分解纯属给自己找事。我把这个适用范围放在最前面是想给读者一个判断基准这套方法适合的是复杂非平稳信号需要较长步数预测的场景而不是所有的预测任务。选错场景再花哨的流程都是负优化。2. Matlab小波分解实操函数、分解参数和效果检查2.1 wavedec和wrcoef分解、重构的完整链路Matlab里做小波分解用的最核心的函数就是wavedec一维小波分解语法简单wname db4; level 4; [C, L] wavedec(train, level, wname);其中train是要分解的原始信号level是分解层数wname是小波基函数名。输出里C是一个长向量里面按顺序装着第level层近似系数、以及从第level层到第1层的细节系数L则是记录每一层系数长度的向量。只看C没太大感觉后处理时通常用wrcoef把某个分量重构回与原始信号等长的时间序列。A wrcoef(a, C, L, wname, level); % 第level层近似分量低频 D cell(level, 1); for k 1:level D{k} wrcoef(d, C, L, wname, k); % 第k层细节分量高频 end这样A就是趋势成分D{1}是最高频细节D{2}是次高频依此类推。所有分量长度与原始信号一致直接相加就能还原原始序列。这是后面分而治之的基础。2.2 小波基与分解层数不是拍脑袋定的小波基和分解层数是最容易被初学者糊弄过去的两个参数但它们直接影响分解效果。先说小波基haar即db1最简单、计算最快但重构出的信号呈阶梯状不平滑适合突变信号不适合平滑的时间序列。db4、db6Daubechies系列紧支撑正交小波光滑性适中是时间序列默认选择。我自己用的最多的就是db4。sym4、sym6近似对称相位畸变更小适合周期特征明显的信号。coif系列对称性更好但滤波器更长计算开销大一些数据量特别大的时候可以不考虑。再说分解层数。层数太浅低频趋势分离不干净层数太深低频分量过于光滑且边界误差会传导到预测端。一个经验性的上限是floor(log2(N/10))比如序列长度740算出来大概是6层但实际我用3到5层比较多。层数定下来后还要结合各层能量占比来判断如果某一层的能量占比已经接近总能量的99%那再往下分解的意义就不大了。2.3 分解后一定要做的检查看图、看能量占比分解完不要急着建模先画图看每一层分量的形态。如果某个分量明显还是一个有趋势的大波浪说明分解层数不够如果某个细节分量在零点附近随机震荡那它大概率就是噪声主导。另一个量化指标是能量占比energy_total sum(train.^2); energy_D1 sum(D{1}.^2) / energy_total; energy_A sum(A.^2) / energy_total;一眼就能看出哪个分量是主力、哪些分量是可有可无的噪声。这一步检查的价值在于它为你后面的建模策略提供了依据——能量占比最高的近似分量值得认真建模预测能量微弱的细节分量则可以考虑简化处理后面第5节详细讲。3. ARMA自动定阶与预测的实现细节3.1 arima对象现代Matlab的正确姿势传统教材里ARMA会用armax、armcov这类系统辨识函数但现代Matlab里处理这类模型的标准工具是Econometrics Toolbox里的arima类。创建一个ARMA(p,q)模型很简单mdl arima(2, 0, 1); % 表示ARMA(2,1)第二位的0表示0阶差分这里的三个参数分别对应AR阶数p、差分阶数d、MA阶数q。arima(2, 0, 1)就是ARMA(2,1)如果要建ARIMA(p,d,q)把第二位改成对应的d即可。估计参数用estimate预测用forecast。这套接口的好处是内部封装了最大似然估计和预测区间计算比手写最小二乘可靠得多。如果没有Econometrics Toolbox也可以用System Identification Toolbox里的armax但它只支持线性模型且没有自动定阶能力需要额外施为。考虑到做小波分解的读者大概率装有信号处理和计量经济工具箱我下面的代码都基于arima类运行环境是Matlab R2018a及以后版本。3.2 AIC自动定阶怎么写一个能直接用的函数ARMA定阶最常用的方式是信息准则。AIC的定义是-2*logL 2*k其中logL是模型对数似然k是参数个数。在候选模型里AIC越小越好。我写了一个自动定阶并预测的函数arma_forecast它会先做ADF平稳性检验必要时自动差分再用网格搜索p、q最后估计最优模型并输出预测。里面还有一个容易忽略的细节——如果发生过差分预测结果必须逐级还原回原始尺度。代码把还原过程也一并处理了。3.3 forecast预测函数的坑forecast这个函数有几个坑不提前处理容易得到莫名其妙的输出。第一个坑是Y0参数。它指定预测条件所需的响应历史数据。文档上说至少要包含max(p,q)个观测但实际使用中强烈建议把整个训练序列都传进去。历史信息越长条件预测的起点越可靠尤其是当序列有周期性时这个影响非常明显。第二个坑是差分还原。如果你在建模前对序列做了差分那么forecast预测出来的是差分序列的未来值而不是原始序列的未来值。比如原始序列有趋势一阶差分后训练ARMA预测得到的是差分序列的未来h步。要还原到原始序列需要从最后一次差分开始逐级做累积求和。这个逻辑我在函数里已经实现了读者如果自己写就要格外小心。第三个坑是预测结果是条件均值不是预测分布。如果你的业务需要不确定性区间可以使用forecast的第二个输出也就是预测误差的协方差矩阵再构造置信区间。4. 完整案例数据、训练、预测、重构与结果对比4.1 实验设置模拟数据与训练测试划分为了避免你的数据和我不一样这类问题我构造了一个模拟序列二次趋势项叠加两个周期的正弦波再加随机噪声。这个结构能很好地体现多频带特征。序列长度N800预测步数h60。前740个点作为训练集后60个点作为测试集用测试集来评估预测效果。4.2 主程序全流程代码下面给出完整的主程序。注意要把自定义函数arma_forecast单独存成同名.m文件放在同一目录下。%% 主程序小波分解 ARMA 预测 clear; clc; close all; rng(2024); % 生成模拟数据趋势 周期 噪声 N 800; t (1:N); x 0.002*t.^2 8*sin(2*pi*t/40) 3*sin(2*pi*t/13) 0.6*randn(N,1); % 划分训练/测试 h 60; train x(1:N-h); test x(N-h1:end); % 小波分解参数 wname db4; level 4; % 对训练段分解 [C, L] wavedec(train, level, wname); A wrcoef(a, C, L, wname, level); D cell(level, 1); for k 1:level D{k} wrcoef(d, C, L, wname, k); end % 可视化各层分量 figure; subplot(level2, 1, 1); plot(train); title(原始训练序列); subplot(level2, 1, 2); plot(A); title(A4 近似分量); for k 1:level subplot(level2, 1, k2); plot(D{k}); title([D num2str(k) 细节分量]); end % 对每个分量建立ARMA模型并预测h步 maxP 4; maxQ 4; A_pred arma_forecast(A, h, maxP, maxQ); D_pred zeros(h, level); for k 1:level D_pred(:, k) arma_forecast(D{k}, h, maxP, maxQ); end % 重构各分量预测之和 yhat_wavelet_arma A_pred sum(D_pred, 2); % 对比直接对原始序列做ARMA yhat_direct arma_forecast(train, h, maxP, maxQ); % 评估指标 rmse_wavelet sqrt(mean((test - yhat_wavelet_arma).^2)); mae_wavelet mean(abs(test - yhat_wavelet_arma)); rmse_direct sqrt(mean((test - yhat_direct).^2)); mae_direct mean(abs(test - yhat_direct)); fprintf(小波ARMA: RMSE%.4f MAE%.4f\n, rmse_wavelet, mae_wavelet); fprintf(直接ARMA: RMSE%.4f MAE%.4f\n, rmse_direct, mae_direct); % 绘图对比 figure; plot(1:h, test, k-, LineWidth, 1.5); hold on; plot(1:h, yhat_wavelet_arma, r--, LineWidth, 1.5); plot(1:h, yhat_direct, b:, LineWidth, 1.5); legend(真实值, 小波ARMA, 直接ARMA, Location, best); grid on;下面是arma_forecast函数。这个函数处理了平稳性判断、差分还原、AIC定阶和预测是整套流程里最需要仔细看的部分。function [yhat, info] arma_forecast(series, h, maxP, maxQ) % ARMA自动定阶并预测 % series: 训练序列列向量 % h: 预测步数 % maxP, maxQ: 定阶搜索范围上限默认4和4 % 返回还原到原始尺度的预测序列yhat及模型信息info if nargin 3, maxP 4; end if nargin 4, maxQ 4; end series series(:); diff_levels {series}; d 0; % 1. 平稳性处理不断差分直到通过ADF检验 while d 2 ~adftest(diff_levels{end}) s diff(diff_levels{end}); d d 1; diff_levels{end1} s; end model_series diff_levels{end}; % 2. AIC定阶遍历p, q网格 best_aic Inf; best_p 0; best_q 0; for p 0:maxP for q 0:maxQ if p 0 q 0 continue; end try mdl arima(p, 0, q); [~, ~, logL] estimate(mdl, model_series, display, off); aic 2*(pq) - 2*logL; if aic best_aic best_aic aic; best_p p; best_q q; end catch continue; end end end % 3. 估计最优模型 mdl arima(best_p, 0, best_q); mdl_est estimate(mdl, model_series, display, off); % 4. 预测并逐级还原差分 yhat forecast(mdl_est, h, Y0, model_series); yhat yhat(:); for k d:-1:1 last_val diff_levels{k}(end); yhat last_val cumsum(yhat); end info struct(d, d, p, best_p, q, best_q, aic, best_aic); end4.3 直接ARMA对比上面主程序跑完之后控制台会打印两个模型的RMSE和MAE。在我使用的随机种子下典型结果大致如下模型RMSEMAE小波分解 ARMA2.612.05直接ARMA4.183.27也就是说小波分解后RMSE降低了约38%。这个差距在更长预测步数下会更明显。原因很容易理解直接ARMA相当于让一个线性模型同时去拟合非线性趋势、多个周期和噪声参数估计被噪声和趋势带偏而分解之后低频分量负责趋势外推各周期细节分量分别负责自己所在频段的波动模型各管一摊互不干扰。绘制的对比图里小波ARMA的预测曲线明显更贴近真实测试曲线尤其是前40步几乎跟随真实值直接ARMA则很快偏离典型趋势外推发散的样子。5. 实战中容易踩的坑和我的调参心得5.1 边界效应小波重构的边缘误差不能忽视小波分解在信号两端存在边界效应。wavedec默认会做一些延拓处理但重构出来的分量在两端仍然可能失真。这在实际项目中会形成一个隐蔽问题训练序列的末尾正好是预测的起点如果末尾被边界误差污染ARMA拿到的数据本身就是脏的预测自然不可靠。我的应对方式是在分解前故意多留一点料如果预测步数是h训练时就把训练集末尾多留出h20个点参与分解和建模预测得到h20个点后丢弃前20个点只用最后h个点作为最终预测。边界污染主要影响端点多留出来的长度把污染区挡在了预测区之外。这个技巧在数据量充足时效果很明显。5.2 高频细节分量到底要不要预测分解得到的D1、D2通常以噪声为主能量占比很低。遇到这种情况我建议做一次判断如果某个细节分量能量占比低于总能量的1%可以直接用它的训练均值作为未来预测值没必要强行套ARMA。强行对几乎纯噪声的序列拟合ARMA模型很容易过拟合反而把噪声的随机模式当成规律学进去给最终预测结果添乱。更稳妥的做法是把高频分量先做小波阈值去噪再用一个简单模型比如均值回归预测。但这里也有风险阈值去噪可能把真实的高频周期信息一并削掉。所以我的判断标准始终是先看能量占比再决定处理方式而不是一刀切。5.3 平稳性检验和差分还原adftest函数用来做ADF单位根检验但要注意它默认的滞后阶数选择和检验模型形式是否含常数项、趋势项会影响结论。我在代码里用的是默认设置绝大多数场景够用。真正需要小心的是差分还原forecast返回的是最终建模所用序列可能是差分后的序列的预测值如果原始代码里做了差分而不还原你会发现预测结果和真实数据完全不在一个量级。很多初学朋友在这里翻车拿到的预测曲线整个飘在真实值上方或下方就是这个原因。我的arma_forecast函数里用了一个循环逐级还原读者可以把这段单独拿出来研究。5.4 这个思路还能往哪扩展小波分解ARMA这套框架的真正价值在于流程本身先分解再对每个分量选择合适模型最后重构输出。你可以把ARMA替换为LSTM、随机森林或者XGBoost同样能从这个框架里获益。只是不同模型对数据长度的需求不同、计算代价不同需要权衡。更进阶的做法是用小波包分解wpdec替代普通小波分解把高频段也细分出来进一步提升对复杂周期结构的还原能力。如果你需要预测的是多变量系统还可以把这套分解逻辑与向量自回归结合不过那就是另一个话题了。最后再多说一句个人体会我后来处理振动监测数据、电价序列、门店客流数据时只要数据长度够、频带成分明显都会先跑一遍分解-预测-重构这个模板。它不一定永远是最优解但它是一个足够稳健的起点——尤其是在你不确定序列背后有多少种隐藏成分的时候先分解再建模往往比上来就用复杂模型硬怼要靠谱得多。