基于CEEMD的LSTM和ARIMA模型干旱预测适用性研究——以新疆为例

丁 严， 许德合， 曹连海， 管相荣

（1.华北水利水电大学测绘与地理信息学院，河南郑州 450046；2.河南省自然资源电子政务中心，河南郑州 450046）

干旱对农业生产、经济运行、现代生活造成的危害与日俱增，也使得在气候变化过程中确保用水安全、能源安全、粮食安全变得更加困难。近百年来，中国陆地区域平均增温0.9～1.5°C，且气温将在未来持续上升，年均降雨量虽未见显著变化，但不同区域的降雨量差异日趋明显，由此可预见大范围干旱的发生频次将会增加、强度将会增强［1-3］。随着极端天气对人类社会影响的日渐显著，如何针对极端天气的发生进行准确评估、监测和分析，成为了国内外学者关注的重点问题。

现阶段，相关研究常使用干旱指数对干旱发生的程度、持续时间和影响范围进行定量评价［4-5］。目前，学界多使用的评价指标有标准化降水指数（Standardized Precipitation Index，SPI）、帕默尔干旱指数（Plamer Drought Severity Index，PDSI）和综合干旱指数（Composite Index，CI）［6-8］。其中SPI可用于多种时间尺度下的干旱分析，干旱分级精度高且仅用降水数据即可计算，因而广泛应用于干旱研究［9-10］。降水量数据和由此计算得到的SPI 具有非平稳、非线性的特征。应用这一数据进行预测，难以达到精准的预测效果。信号分解能够提取序列的局部特征并使序列平稳，国内外学者通过经验模态分解（Empirical Mode Decomposition，EMD）、集合经验模态分解（Ensemble Empirical Mode Decomposition，EEMD）、互补集合经验模态分解（Complementary Ensemble Empirical Mode Decomposition，CEEMD）对时间序列进行分解，得到了一组较为平稳的分量和一个趋势项，降低了原始时间序列的复杂度，提高了可预测性［11-13］。在干旱预测的过程中，用于预测的模型有很多，如差分自回归移动平均模型（Autoregressive Integrated Moving Average，ARIMA）、人工神经网络（Artificial Neural Network，ANN）、支持向量机（Support Vector Machine，SVM）等，其中ARIMA 模型是最常见的用于时序预测的模型［14-15］。随着机器学习的发展，长短期记忆（Long Short-Term Memory，LSTM）网络在时间序列预测中得到了应用，LSTM 在处理具有很长间隔和延迟的序列上具有优势［16-17］。单一模型在时间序列的预测中容易出现局部最优的情况，预测效果不理想，因此，许多学者将信号分解与预测模型组合用于时序数据的预测，例如EMD-LSTM［18］、EEMD-ARIMA［19］、EEMDLSTM［20］、CEEMD-LSTM［21］均得到了较好的预测结果。目前，对于组合模型预测结果适用性的评价和对比大多是组合模型与传统ARIMA模型的对比［20,22］，缺乏组合模型之间的对比、组合模型与LSTM 的对比。新技术新方法是否优于传统方法仍待考证。CEEMD 解决了EMD模态的混叠问题以及EEMD模态的残留白噪声问题，因此，本文基于CEEMD 构建CEEMD-ARIMA 组合模型和CEEMD-LSTM 组合模型。分别通过ARIMA、LSTM、CEEMD-ARIMA 和CEEMD-LSTM模型进行预测，对其结果进行分析对比，研究其在干旱预测中的适用性。

本文选取新疆32 个站点的1960—2019 年逐日降水量数据，计算1、3、6、9、12 个月及24 个月时间尺度SPI。利用ARIMA、LSTM、CEEMD-ARIMA 和CEEMD-LSTM组合模型对各SPI序列进行预测。通过对4 种模型预测结果和实际计算值的对比，结合决定系数（Coefficient of Determination，R2）、均方根误差（Root Mean Square Error，RMSE）、平均绝对误差（Mean Absolute Error，MAE）3 种评价指标，分析4种模型的干旱预测精度。结合ArcGIS 的经验贝叶斯克里金插值法，展示4 种模型预测的干旱空间分布情况。从模型预测结果的精度和空间分布情况探索模型在干旱预测中的适用性，以期能为气象防灾减灾工作提供决策依据，减少旱灾损失。

1 数据与方法

1.1 研究区概况及数据来源

新疆地处欧亚大陆腹地，地理坐标位为73°40′～96°18′E、34°25′～48°10′N，自北向南有阿尔泰山、天山和昆仑山系，呈“三山夹两盆”的地貌格局。该区远离海洋，降水稀少，干旱频发，是典型的干旱半干旱地区。研究区域的地理位置及气象站点分布如图1所示。本文所用的逐日降水量数据来源于国家气象科学数据中心（http://data.cma.cn/）中新疆气象站观测数据。所用新疆地理高程数据来源于地理空间数据云（http://www.gscloud.cn/search）。

1.2 研究方法

1.2.1 标准化降水指数 降水量是影响干旱的重要因素。标准化降水指数考虑了降水量分布为偏态分布的情况，假定降水量分布服从Γ分布，计算出降水量的分布概率，之后进行正态标准化处理，将处理得到的结果依据气象干旱等级（GB/T20481-2017）中的干旱分级标准，进行干旱等级划分（表1）。SPI 能够计算出不同时间尺度的值，满足多种水资源状况监测的需要，其中1、3、6、9、12、24 个月时间尺度下的SPI 可用于描述区域的气象干旱、农业干旱、水文干旱情况［23-25］。SPI 易于计算，具体计算过程参见气象干旱等级（GB/T20481-2017）。

1.2.2 CEEMD 分解 1998 年，Huang 等［26］提出了EMD，EMD 在处理非线性、非平稳信号上具有优势。原始序列输入EMD 进行分解能够得到有限个固有模态函数（Intrinsic Mode Function，IMF）和趋势项，各分量包含了原始序列在不同尺度上的局部特征。经过EMD分解后的结果具有相当高的信噪比，但这种分解方法存在模态混叠的问题。EEMD作为EMD的进一步改进，通过向原始信号添加高斯白噪声，有效减少了模态混叠的发生，但白噪声的添加，使各分量含有残留白噪声［27］。Yeh 等［28］提出了CEEMD，通过向原始信号中添加n组符号相反的白噪声，减少分量数据中噪声的残余量，达到残余白噪声可以忽略不计的目的，其算法步骤如下［28-29］：

（1）向原始序列B(t)中加入n组包括正噪声和负噪声的辅助白噪声，从而得到正噪声序列H1和负噪声序列H2，此时得到的序列总数为2n。

（4）将得到的IMF 值作为最终分解结果，即原始序列分解为：

式中：r(t)为残留趋势项。

1.2.3 LSTM网络 LSTM网络是一种特殊的循环神经网络（Recurrent Neural Network，RNN），能够学习数据传递中长期依赖的信息，并有效解决梯度问题。LSTM 网络有着比RNN 更复杂的重复模块（图2），其中σ、tanh 分别为sigmoid 函数和双曲正切函数。细胞状态是这个重复的神经网络模块链的关键，即穿过每个模块的水平线，它类似于传送带，贯穿了整个链条，保证了信息传输的不变性。通过“门”，LSTM 向细胞状态添加或移除信息。遗忘门决定了要从细胞状态中移除哪些信息，这是由1 个sigmoid层决定的。输入门用来更新状态信息，由两部分组成，通过sigmoid 层决定哪些信息需要更新，并在tanh 创建1 个包含新的待添加信息的向量，由此对细胞状态进行更新。输出门用sigmoid 层决定了要输出的细胞状态的部分［16］。通过运算（图2 的圆圈部分），将结果继续传递给下1个单元结构。

图2 LSTM结构图Fig.2 Structure diagram of LSTM

1.2.4 ARIMA 模型 Box 等［30］提出了能够进行非平稳非白噪声序列预测的ARIMA模型，通过d次差分使序列平稳，然后利用自回归滑动平均（Autoregressive Moving Average，ARMA）模型预测。ARMA模型假定原始序列为一组随机序列，通过改变模型的参数对该序列近似描述，选出最符合该序列的模型参数，之后依据原始数据对未来情况进行预测［31］。ARIMA（p,d,q）模型的一般式为［15］：

式中：Yt为时间序列值；ωi（i=1,2,…,p）和θj（j=1,2,…,q）分别为自回归系数和滑动平均系数；ut为白噪声序列，且ut～N(0,σ2)。

ARIMA模型的建模流程为：

（1）平稳性检验。本文通过单位根检验（Augmented Dickey-Fuller Test，ADF）判断时间序列的平稳性［32］。若为非平稳时间序列则需对原始序列d次差分。

（2）确定模型阶数的取值范围。根据数据的自相关函数（Autocorrelation Function，ACF）和偏自相关函数（Partial Autocorrelation Function，PACF）确定p、q的取值范围。

（3）模型定阶。利用赤池信息准则（Akaike Information Criterion，AIC）、贝叶斯信息准则（Bayesian Information Criterion，BIC）对模型定阶，AIC、BIC 公式如下：

式中：N为参数个数。选择AIC、BIC 值最小时对应的p、q值。

1.2.5 基于CEEMD 的组合模型 波动性强的原始序列经过CEEMD 分解，能够得到一组波动较低的IMF分量，这提高了序列的可预测性。通过Python，将CEEMD 分别与LSTM 和ARIMA 模型结合组成CEEMD-LSTM 组合模型和CEEMD-ARIMA 组合模型。通过组合模型进行预测的步骤如下：

（1）CEEMD 分解。通过CEEMD 对原始SPI 序列进行分解，得到从高频到低频的IMF1、IMF2、···、IMFn以及Res。

（2）LSTM 或ARIMA 模型预测。将IMF1、IMF2、···、IMFn以及Res 分别导入LSTM 或ARIMA模型进行预测，预测结果分别记为P1、P2、···、Pn+1。

（3）对预测结果相加求和。

基于CEEMD的组合模型建模流程如图3所示。

图3 组合模型建立流程Fig.3 Workflow of combined model

1.2.6 评价指标 本文选取RMSE、MAE、R2作为4种模型的评价指标。RMSE 和MAE 的取值范围为［0，+∞］，值越小，模型效果越好。R2越大，表示拟合效果越好，最大值为1。

2 结果与分析

2.1 LSTM网络模型训练及预测

本文以库尔勒站点为例，利用LSTM 网络模型对1、3、6、9、12 个月及24 个月时间尺度SPI 序列进行建模，步骤如下：

（1）数据归一化处理

对输入的SPI 数据进行归一化处理，以提高模型的训练速度。

（2）网络模型训练

LSTM 网络的激活函数通常有sigmoid、tanh 和ReLU。sigmoid 存在着随神经网络层数加深，梯度后向传播到浅层网络时易出现梯度消失的缺点；tanh也存在梯度消失的情况，且sigmoid和tanh的随机梯度下降收敛速度较慢，因此激活函数选用了ReLU。1次训练选取的样本数为1，即每训练1个样本，更新1 次权重。损失函数则采用均方误差（Mean Squared Error，MSE），优化算法采用了Adam。通过“早停法”防止训练过拟合，即随着迭代次数增加，MSE逐渐下降，模型精度逐渐提高；当MSE值上升时，停止训练。为确保模型精度达到最高，迭代次数设置为300。采用了黄金分割法选择隐藏神经元数量，隐藏层神经元数为25［17］。

（3）输出预测数据

由于之前对数据进行了归一化处理，因此，此处需要采取反归一化处理，以得到模型的实际预测数据（图4）。

2.2 ARIMA模型建模及预测

依据32 个气象站点1960—2019 年的逐日降水量数据进行SPI 计算。不同时间尺度的SPI 适用于干旱研究的不同方面，因此本文计算了1、3、6、9、12、24 个月共6 个时间尺度的SPI。将计算得到的SPI 中1960—2007 年数据作为训练集，2008—2019年数据作为测试集。本文以库尔勒站点为例对ARIMA 建模，在预测前，需要对测试集数据的平稳性进行判断。若数据平稳，则可通过ARMA 模型进行预测；若不平稳，则需进行差分，ADF检验结果见表2。表2中6个时间尺度SPI的P值均小于0.05，即时间序列均为平稳时间序列，因此，可进行下一步。

表2 原始序列单位根检验Tab.2 ADF test of the original sequence

通过ACF、PACF确定各时间序列p、q的可能取值。利用AIC、BIC准则选取最优模型。各序列的模型定阶结果见表3。分别通过6 个时间尺度SPI 的最优模型进行预测，预测结果见图4。

表3 6个尺度SPI的ARIMA模型定阶Tab.3 ARIMA model order based on SPI values of six time scales

图4 LSTM、ARIMA、CEEMD-LSTM与CEEMD-ARIMA模型多时间尺度SPI预测（2008—2019年）Fig.4 Forecast of multi-time scale SPI of LSTM,ARIMA,CEEMD-LSTM and CEEMD-ARIMA model(2008-2019)

2.3 利用组合模型对SPI序列进行预测

经过参数的多次修改和对比，最终选定将Nstd设置为0.2，NE 设置为100，TNM 设置为8。利用CEEMD分解多尺度SPI，得到8个IMF分量和1个趋势项。以SPI3分解为例，原始序列和分解得到的子序列见图5。由图5 可知，原始序列波动范围较大，而分解得到的IMF 分量波动范围较小，随着分解的逐步进行，分量的波动趋于平缓，说明通过CEEMD分解能够降低原始序列的非平稳性。

图5 CEEMD分解SPI3序列Fig.5 CEEMD decomposition results of SPI3 sequence

选取1960—2007 年数据作为训练集，2008—2019年数据作为测试集。利用组合模型进行预测，预测结果见图4。由图4 可知，在1 个月时间尺度下，LSTM 和ARIMA 模型的预测值与实际观测计算值相差较大。CEEMD-LSTM 和CEEMD-ARIMA 组合模型的预测值与实际值则较接近，其中CEEMDARIMA 能准确预测到2011 年的干旱发生强度。在3个月尺度下，2个单一模型的预测值与实际值差距缩小，预测的SPI变化趋势与实际趋势相符。此时，CEEMD-ARIMA 已能准确预测2011 年和2017 年的干旱情况，整体预测结果与实际情况最为一致。在6 个月尺度下，与LSTM 相比，ARIMA 模型对干旱发生时间和强度的预测更为准确。4 个模型中，ARIMA 和CEEMD-ARIMA 模型对干旱的预测较精准。在9个月尺度和12个月尺度下，除LSTM外的其他3种模型预测情况接近实际情况，较1、3、6 个月尺度下，对干旱事件的发生及强度和持续时间的预测更为准确。在24 个月尺度下，4 种模型的预测结果与实际情况近乎一致，从干旱发生强度的预测情况来看，ARIMA和CEEMD-ARIMA模型的预测结果分别优于LSTM和CEEMD-LSTM模型。对模型在6个时间尺度SPI的预测结果进行对比，在1个月时间尺度下，4 种模型的预测结果均为6 个时间尺度中最差的，与实际结果相差最大。随着时间尺度的增大，4种模型预测的准确性有所提升。

通过R2、RMSE、MAE 共3 种评价指标对预测结果进行评价，进一步分析4种模型的预测精度。表4中LSTM 在SPI1 的RMSE、MAE 值分别为0.8681 和0.6478，随着时间尺度的增加RMSE、MAE 值逐渐减小。在24 个月时间尺度下达到最小，SPI24 的RMSE、MAE 值分别为0.4266 和0.2700。R2值则呈现相反趋势，表明随着时间尺度增大，模型的预测精度逐渐提高。ARIMA、CEEMD-ARIMA、CEEMDLSTM 模型预测精度随时间尺度的变化趋势与LSTM一致。对各时间尺度SPI进行预测，ARIMA模型预测结果的R2值均略高于LSTM，RMSE、MAE 的值则均略低于LSTM，说明ARIMA 模型的预测精度优于LSTM。CEEMD-LSTM 和CEEMD-ARIMA 模型的R2值在各时间尺度均高于单一模型，LSTM、ARIMA、CEEMD-LSTM 和CEEMD-ARIMA 模型在SPI24的R2值分别为0.8882、0.9103、0.9403 和0.9846。其中，CEEMD-ARIMA模型除对SPI1的预测结果外，R2值均在0.8 以上，具有较高的预测精度。在各个时间尺度下，预测精度从低到高为：LSTM、ARIMA、CEEMD-LSTM、CEEMD-ARIMA 模型，说明ARIMA的预测精度高于LSTM，CEEMD 能够有效提高模型的预测精度。

表4 4种模型预测结果的R2、RMSE、MAE值Tab.4 R2,RMSE and MAE values of the predicted results of four models

使用ArcGIS 对32 个站点在2019 年SPI 的实际观测计算值和预测值进行可视化展示（图6）。由于新疆的干旱在一年四季皆有发生，此处选择能够进行降雨量季节变化分析的SPI3 对区域干旱情况进行展示。从图6 中可以看出，CEEMD-ARIMA 组合模型对干旱空间分布情况的预测与实际情况最为接近。2019 年2 月的北疆降水量偏多，全疆其余大部分偏少。4种模型在冬季的预测情况与实际情况都存在着偏差，其中CEEMD-ARIMA 组合模型的预测结果与实际计算结果较为一致。

图6 使用克里金插值对实际值和4种模型的预测结果可视化展示Fig.6 Kriging interpolation results of the actual calculated values and the predicted values of four models

3 讨论

SPI 时间序列是非平稳序列，而单一模型预测结果的精度受原始数据平稳性影响较大。Liu 等［33］利用ARMA 对山东省5 个站点的SPI9 序列进行预测，预测结果的平均相对误差最低为20.39%，最高为43.69%，预测精度较低且不同站点间存在很大差异。单独通过LSTM预测SPI，同样有着较差的预测结果［34］。CEEMD分解能够为模型预测提供平稳性，从而提高序列的可预测性［13］。通过CEEMD 分解，原始序列在不同尺度的局部特征被提取出来，非平稳时间序列转化为平稳的分量。因此，本研究利用CEEMD 降低SPI 序列的非平稳性，确保LSTM 和ARIMA模型能够有效预测SPI序列。

在4 种模型的预测结果中，SPI1 的预测精度相较于其他5个时间尺度最差。数据的平稳性与预测结果有密切关系，1个月时间尺度的数据量是6个时间尺度中最大的，并且数据序列趋于严平稳（序列的分布结构不随时间改变），随着时间尺度的增大，数据量减少，并且数据序列趋于宽平稳（未来值与过去值相关），模型的预测情况变好。LSTM在高频序列的预测中具有较高的预测精度，ARIMA在低频序列中有较好的预测效果，因此，LSTM和ARIMA模型分别适用于高频序列和低频序列的预测，同时也造成了LSTM 在SPI 序列预测中预测效果略差于ARIMA模型［16］。CEEMD分解得到的子序列可预测性要高于原始序列，因此，在1个月时间尺度下，2个组合模型的预测情况明显优于单一模型的预测情况。在3 个月和6 个月时间尺度下，组合模型优于单一模型。随着时间尺度的增大，优势逐渐缩小，长时间尺度的SPI序列集合了原始数据中更多的信息，整个序列趋于平稳，单一模型的预测精度随之提高。

SPI易于计算，且能够描述地区的气象干旱、农业干旱、水文干旱情况，但对于新疆这一研究区而言，SPI 具有一定的局限性。新疆农业所耗水分不仅来源于降水，也来源于当地的灌溉用水。地下水位的变化与山区河流径流及新疆农业耗水有着很大的关系。干旱的发生是多种因素的共同作用，除降水外，需要考虑的因素还有很多。在年降水量未有显著变化的情况下，随着全球温度的逐渐上升，干旱发生的频次势必会增加。因此，若只考虑降水因素的影响，干旱发生的预测将会变得越来越困难，还需在研究中考虑多种因素的干旱指数在干旱预测中的适用性。

4 结论

本文分别利用LSTM、ARIMA、CEEMD-LSTM和CEEMD-ARIMA 模型对1、3、6、9、12 个月及24 个月时间尺度的SPI 进行预测，通过对预测结果的对比分析，主要得到以下结论：

（1）4 种模型预测精度随时间尺度的增大而提高，即在1个月尺度下最低，在24个月尺度下最高，此时R2值均在0.85 以上，表明4 种模型在干旱预测中的适用性随着时间尺度的增大逐渐提高。

（2）CEEMD-LSTM和CEEMD-ARIMA组合模型在1、3、6、9、12 个月及24 个月时间尺度下，均有着比单一模型更高的精度。说明CEEMD 在处理非平稳、非线性数据上具有优势，通过CEEMD 分解，原始数据序列变得平稳，序列的可预测性提高。

（3）CEEMD-ARIMA 模型的预测精度最高，除SPI1外，其余5个时间尺度的R2值均在0.80以上，且在SPI24 时达到了0.98。CEEMD-ARIMA 模型预测的干旱空间分布情况与实际情况较为吻合，说明CEEMD-ARIMA 模型能够很好地拟合不同尺度的SPI序列，适用于干旱预测。