基于EEMD-ARIMA的年降水预测拟合模型研究

李智强 邹红霞 齐 斌 郭江昆

(航天工程大学航天信息学院 北京 101400)

0 引 言

年降水量变化一直都是气象科学界对气候变化研究的热点问题。年降水量的变化对河流径流量有着直接影响，同时还紧密地关系到农业发展和粮食安全，对年降水量序列变化的研究可以为洪涝、干旱灾害的预测预防提供参考，为农业生产提供有利条件。年降水量变化受到地理、季风、海洋和太阳等多方面的影响[1]，影响降水的特征指标众多，比如高原积雪、海表温度、ENSO、季风等多种特征因子[2]，其内部机理异常复杂，因此目前大部分的年降水量预测方法主要为气象学方法。虽然现有的气象学预测方法已经有了广泛的应用[3-5]，但是由于气象学方法的自变量难以确定，模型的计算量也异常巨大。时间序列的ARIMA分析方法不用考虑复杂的影响因素，计算量小，因此在降水研究中也有应用[6-7]，但是ARIMA方法对于非线性序列的拟合效果不好，而年际降水序列大多都是非线性的[8]。

为解决上述问题，有学者提出了基于集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)的ARIMA时间序列模型[9]，主要得益于EEMD能够自适应地处理任意非线性非平稳时间序列的特点。文献[10]的研究表明，EEMD-ARIMA模型在短期径流上的预测精度高于传统ARIMA模型。

1 相关模型

1.1 EMD和EEMD

由于EEMD是在经验模态分解EMD基础上的改进，因此先介绍EMD的方法原理。EMD比小波和傅里叶分解变换具有更好的时频分辨率和自适应的特点，且其在非平稳、非线性的信号时频序列处理上表现更加优秀，该方法在海洋、大气和天体观测等领域有应用。

经验模态分解的关键就是将复杂的信号分解为有限个本征模函数(Intrinsic Mode Function，IMF)，其计算方法如下：

首先取原始序列信号X(t)的所有极值点，通过三次样条插值函数[11]分别连接所有极大值点和所有极小值点，形成上包络线和下包络线，使得整个序列信号都在上下包络线中，将原序列减去上下包络线的均值m1，得到新的序列h1(t)，即：

X(t)-m1=h1(t)

(1)

一般来说，h1(t)还不是IMF分量，通过一次的筛选不一定能得到平稳序列，因此需要重复以上过程，第二次筛选将h1(t)作为被分解序列。假设h1(t)的包络均值为m11，则得到的新的序列h11(t)，即：

h1(t)-m11=h11(t)

(2)

筛选的目的是消除噪声波，使波形轮廓对称，直到第k次的筛选，得到的波形第一个满足IMF定义的分量h1k(t)，将第一个IMF分量h1k(t)用c1表示。

c1(t)=h1k

(3)

一个原始信号一般由多个IMF分量组成，由式(3)可得剩余分量r1(t)，然后将剩余分量r1(t)作为新的被分解信号，此过程一直重复，直到满足以下条件之一：(1) 最后一个剩余量rn(t)或IMF分量cn(t)小于某个给定的常量；(2) 剩余量rn(t)为单调函数。

X(t)-c1(t)=r1(t)
r1(t)-c2(t)=r2(t)

(4)

⋮
rn-1(t)-cn(t)=rn(t)

最终原始信号被分解为n个不同周期和不同频率的IMF分量和一个残差量：

(5)

虽然EMD基于数据本身变化，使得IMF与采样频率相关，克服了傅里叶变换的局限性，但是EMD存在的一个重要问题就是IMF分量存在模态混叠[12]，为了更好地解决这一问题，Wu等[13]提出了集合经验模态分解(EEMD)。

EEMD是基于噪声辅助分析的方法，原理较为简单，通过在年降水序列中加入白噪声，白噪声频谱均匀分布，使得信号自动分布到合适的尺度上。由EMD的零均值特性，噪声在多次平均计算后会相互抵消，从而得到的集成均值计算结果即可视为最终结果，最终结果与原始序列的误差随着集成平均次数的增加而减小。因此从理论上讲，EEMD对抑制EMD的模态重叠有较好的效果。EEMD的算法流程如图1所示。

图1 EEMD算法流程示意图

1.2 ARIMA模型

通过EEMD分解得到不同周期的IMF分量后，对分量建立拟合模型。差分自回归移动平均模型ARIMA是在ARMA模型针对非平稳序列模型拟合效果差的基础上作出的改进，通过将非平稳序列进行一次或多次差分转换成平稳序列，再进行ARMA拟合，其构成如下：

ARIMA(p,d,q)=AR(p)+Difference(d)+MA(q)

(6)

式中：AR(p)、Difference(d)、MA(q)分别为自回归模型、差分模型和移动平均模型；p、d、q分别为自回归项、差分阶数和移动平均项数，这三个参数决定了模型的好坏。ARIMA模型预测方程如下：

xt=φ0+φ1xt-1+φ2xt-2+…+φpxt-p+

εt+γ1εt-1+γ2εt-2+…+γqεt-q

(7)

式中：xt为样本值；εt为当期随机误差干扰；φi和γj为模型参数；p、d、q是模型阶数。xt为前p阶xt-1,xt-2,…,xt-p和前q阶εt-1,εt-2,…,εt-q的多元线性函数。ARIMA的建模流程如图2所示。

图2 ARIMA建模流程

2 EEMD-ARIMA模型

2.1 模型框架

EEMD是在对EMD做白噪声统计实验时提出的，由于EEMD分解克服了EMD分解过程中产生模态混叠的问题，因此在序列拟合预测中更具优势。本文针对降水时间序列的特点，首先对降水序列EEMD分解，然后对分解后的序列做白噪声检验，一般EEMD会将高频的白噪声分量直接删除，但这样会直接影响模型拟合的精度。本文的做法是保留分离出来的白噪声分量，对其也建立ARIMA模型以降低拟合误差，对其他的分量进行平稳性判别，然后对平稳序列建立ARMA模型，对非平稳序列建立ARIMA模型，最后对剩余模型拟合分量重构和评价。

本文基于EEMD-ARIMA的年降水序列模型，其整体框架如图3所示。首先通过EEMD分解得到n-1个分量和一个剩余分量，通过白噪声检验后单独对白噪声分量建立ARIMA模型，然后对各IMF分量进行平稳性判别，平稳序列采取ARMA模型拟合，非平稳序列采用ARIMA模型拟合，最后对各个分量进行叠加重构并进行模型评价。其模型框架如图3所示。

图3 EEMD-ARIMA模型框架

2.2 评价指标

为了验证EEMD-ARIMA模型拟合效果是否优于EMD-ARIMA模型以及单一的EMD、ARIMA、EEMD模型，引入均方根误差(RMSE)和平均相对误差(MRE)两个指标来检验模型的预测效果。均方根误差(RMSE)能够从总体上反映预测值的离散程度，其公式如下：

(8)

平均相对误差(MRE)能够反映模型误差的大小以及模型预测的准确程度，其计算公式如下：

(9)

在实际的年际降水预测中，通常采用趋势分析的方式，离精准预测的水平还相差甚远。由于中长期预报的复杂性，当预测值与实际值的相对误差低于20%时即视为正确[14]，正确率公式如下：

(10)

3 模型实证分析

首先通过数据预处理得到原始年降水序列，然后根据EEMD-ARIMA模型框架进行模型的建立和训练，为了检验该模型的应用效果，还建立了EMD、EEMD分解-重构模型、ARIMA模型和EMD-ARIMA模型，以便进行对比分析。

3.1 数据来源及实验环境

本文的实验数据来源于《中国地面气候月值数据集3.0》，选取[52984]临夏地面气象观测站1953年—2012年共60年的逐日降水资料，通过数据预处理得到60年的年降水序列。由于时间序列预测精度随时间的增加而降低的特性，为了保证模型的稳定性和可靠性，选取前57年数据进行建模训练，后3年数据用于预测。实证分析使用Python开发语言，使用Jupyter Notebook交互式编程平台，ARIMA模型构建选取Statsmodels库。

3.2 数据预处理

《中国地面气候资料月值数据集3.0》中的文件命名为：SURF_CLI_CHN_MUL_MON-YYYY.txt，其中“SURF_CLI_CHN_MUL_MON”为数据集代码，文件按年份拆分，文件名中的“YYYY”表示年份。

数据集包括了中国752个基本、基准地面气象观测站及自动站自1951年以来的气候资料月值数据集。本文从数据集中选取了降水数据集，其格式如图4(a)所示；降水数据集中文件内部数据形式如图4(b)所示；截取出[52984]临夏气象站共60年的逐日降水量数据，其数据格式如图4(c)所示，可见数据中存在“32766”、“32700”等异常值，分别表示缺测、微量数据，因此均选取零值填充；通过缺失填充、数据集成，将日降水量以年为单位分别累加求和，最终得到年降水序列，如图4(d)所示。

图4 数据预处理过程

3.3 模型建立

1) EEMD的分解。根据EEMD的算法流程，将57年的年降水序列按照添加噪声、EMD分解、计算集成的顺序，通过Python语言实现EEMD的分解，部分代码如下：

eemd=EEMD()

eemd.trials=50

eemd.noise_seed(12345)

E_IMFs=eemd.eemd(S, T, max_imf)

imfNo=E_IMFs.shape[0]

最终得到4个IMF分量序列(IMF1-IMF4)和1个剩余分量IMF5，如图5所示，其中：横坐标为时间，纵坐标为降水量。

图5 EEMD年降水序列分解图

可以看出：IMF1频率最高，IMF2、IMF3、IMF4频率依次降低，IMF5为残余分量，为递减的单调函数。最高频率的IMF1代表时间周期较短的随机波动，一般作为噪声处理。但是噪声的直接删除会导致拟合和预测的精度大大降低，因此同样将随机波动也进行拟合，能够降低模型误差。

2) 白噪声检测。白噪声检测又称为纯随机性检测，通过白噪声检测，检测各IMF分量中是否含有没有意义的纯随机序列，一般对于纯随机序列，都采取舍弃的原则。常用的白噪声检测方法包括Ljung Box假设检验和Box Pierce假设检验[15]。除此之外，通过自相关和偏自相关函数图，也能从经验上判断出序列是否是白噪声序列。首先通过各个IMF分量的自相关和偏自相关图，可以确定是否存在白噪声序列，如图6所示。由于IMF1分量的自相关和偏自相关阶数几乎都落在置信区间内，其自相关性较小，因此判断IMF1分量为可能的随机性序列。

图6 IMF分量自相关和偏自相关图

3) 平稳性判别。判断序列的平稳性是为了建立更好的模型，平稳性判断除了根据图5进行主观判断外，常用扩展迪基-福勒检验(Augmented Dickey-Fuller，ADF)[16]进一步确定平稳性。对各个IMF的ADF检验结果如表1所示。

表1 ADF 平稳性检验结果

可以看出，IMF1、IMF2、IMF3的T统计值均小于各自1%显著性水平，同时P值接近0，IMF1、IMF2、IMF3分量为平稳序列，IMF4、IMF5的T统计值均大于10%显著性水平，IMF4、IMF5为非平稳序列，与图5的结果一致。因此，对前3个IMF分量用ARMA模型拟合，后两个IMF分量用ARIMA模型进行拟合。

4) 模型拟合。在年际降水预测中，大多学者选择使用单一ARIMA时间序列模型。本文在考虑时间序列数据含有噪声和多重频率特性的基础上，将序列进行EMD分解，将复杂序列变为多个简单序列和一个残差序列。由于EMD分解会产生模态重叠，因此使用EMD的改进算法EEMD，对得到的简单序列和残差序进行白噪声和平稳性检验后，得到平稳分量IMF1、IMF2、IMF3和非平稳序列IMF4、IMF5，并分别采用ARMA模型和ARIMA模型进行拟合，通过参数优化。建模及优化部分代码如下：

order0=st.arma_order_select_ic(imf0,max_ar=5,max_ma=5,ic=[′aic′, ′bic′, ′hqic′])

arma_mod10=ARMA(imf0,(2,0)).fit()

pimf0=arma_mod10.fittedvalues

ax1.plot(imf0)

笔者做完自我介绍和简单的开场白，考虑和新员工年龄差距较小，设计的PPT可以自动播放和用激光笔远程控制，对于所讲内容成竹在胸，于是就尝试走进新员工中间去说，拉近了彼此之间的物理距离。这个尝试，有比较好的收效，全程半个小时的授课过程中，基本没有看到新员工们有低头看手机的现象。

ax1.plot(pimf0,′--′,color=′red′)

对IMF1、IMF2、IMF3、IMF4、IMF5分别选取ARMA(2,0)、ARMA(4,2)、ARMA(2,4)、ARIMA(5,0,1)、ARIMA(1,1,1)5个模型进行拟合，拟合结果如图7所示。

图7 分量拟合结果

可以看出，ARMA对于IMF1分量的拟合效果还不够理想，相比于直接舍弃高频分量，该IMF1分量将大大降低对拟合结果的误差，IMF2、IMF3、IMF4、IMF5分量的拟合效果较好。EEMD和EMD都具有完备性，即将全部分解的分量直接相加就能重构出原信号。因此通过对各个分量的重构，就能得到非常接近原始序列的新序列，误差接近0。从模型训练的结果来看，低频的分量拟合效果比高频分量的拟合效果好。

3.4 模型对比分析

3.4.1模型结果

为了得到直观的对比，同时对同一年降水序列建立了EMD模型、ARIMA模型、EEMD模型和EMD-ARIMA模型，其中EMD模型和EEMD模型均采用了分解-删除高频分量-重构的方法，ARIMA模型是直接对年降水序列建立的，经过参数调优，设定为ARIMA(1,0,1)。得到的各个模型最终结果对比如图8所示。

图8 各模型结果对比

通过图8的模型结果对比分析，有以下4点结论：(1) 5个模型都能反映出真实年降水序列的波动趋势；(2) 从EEMD模型和EMD模型的拟合结果看，二者的拟合效果都比较平滑，且EEMD的拟合效果要略好于EMD模型；(3) EMD-ARIMA模型和EEMD-ARIMA模型拟合效果均好于单一模型，单一模型的拟合值波动幅度范围偏小；(4) EEMD-ARIMA模型和EMD-ARIMA模型结果的走势大致相同，但EMD-ARIMA在极值点处的拟合效果更加激进，单从拟合效果图分析，还不能确定两类模型的优劣。

3.4.2模型评价

5个模型的拟合结果的误差分析如表2所示。

表2 模型评价结果

可以看出EEMD-ARIMA模型的拟合效果最好，其RMSE和MRE均为最小值。本文提出的模型效果较EMD-ARIMA均方根误差小14.53，平均相对误差小0.006。在5个模型中，EEMD-ARIMA的预测准确率最高，达到了82.46%。综合三种模型评价指标来看，EEMD-ARIMA模型的预测结果是最好的。

4 结 语

本文运用结合了EEMD分解的算法和ARIMA时间序列分析模型，针对年降水时间序列进行了分析建模，初步探索了EEMD-ARIMA模型在年降水时间序列的分析应用方法。该模型对年降水量时间序列的预测指导性较高，在年际降水预测中具有实际应用的意义，可以应用于中长期天气业务工作中。

本文对于EEMD分解后的高频分量模型选取简单，未来可选取神经网络模型，可能对高频分量拟合得更好。同时，本文使用的EEMD算法解决了模态重叠问题，但还存在分解出现的端点效应，未来可以进一步探索有效的分解算法，提高准确率。