基于季节性自回归积分滑动平均与深度学习长短期记忆神经网络的降水量预测

张丽婷，李鹏飞，庞文静，惠雯，秦孟晟

(1.江苏省扬州市气象局，扬州 225600；2.哈尔滨工业大学电气工程及自动化学院，哈尔滨 150030；3.中国气象局气象探测中心，北京 100081；4.国家气象卫星中心，北京 100081)

中国幅员辽阔，跨越多个气候带，人口活动分布范围较广，每年都会因为气象灾害造成人员伤亡和财产损失[1-2]。相关统计结果表明：中国主要的气象灾害分别是干旱、暴雨、洪涝等。而暴雨作为主要气象灾害之一，其主要的评价指标是降水量。特别是在近年来受全球气候变暖影响，极端灾害事件发生频率增加的背景下[3]，如何准确地预测季节性降水量并进行数据评估，对于汛期的防汛准备工作具有重要指导意义。

影响区域降水量的变化因素繁多，比如地形、气温、日照以及河流等。而降水的内部机理和形成机制也比较复杂，从机理入手的长期预测难度比较大。因此，对于不同时间尺度的降水量的预测问题，国内外学者主要采取以时序序列的统计学分析为主的方法[4-21]。余霖等[4]通过一种十字交叉选择算法，构建一种基于降水量序列的平稳性与周期性的概率预测模型。周明圆等[5]采用Morlet小波分析与Hurst指数法对长江源区近48年来降水量进行了时空分析与预测。张先起等[6]以郑州市年降水量为例，通过算法改进构建了一种CEEMD-Elman(complementary ensemble empirical mode decomposition-Elman)模型，并将其应用于该地区的年降水量预测。龙云等[7]采用粒子群优化的方法构建小波神经网络模型，对洞庭湖流域的月降水量进行了应用评估。王充等[8]以固原市原州区为例，采用Mann-Kendall检验、滑动t检验等方法对该地区近60年的降水量进行了时空与突变分析。郭家力等[9]通过结合快速傅里叶变换(fast Fourier transform，FFT)与集成经验模态分解(ensemble empirical mode decomposition，EEMD)方法对鄱阳湖流域的降水量周期进行了预测研究。胡虎等[10]通过将反向传播神经网络(back propagation neural network，BPNN)、广义回归神经网络(generalized regression neural network，GRNN)与EEMD-GRNN模型进行对比，提出EEMD-GRNN模型在降水量预测方面具有更高的准确度。苗正伟等[11]通过将加权马尔可夫链和改进型的模糊聚类方法结合的方式，对承德市年降水量进行了综合预测研究。熊俊楠等[12]通过将长序列趋势拟合与灰色预测模型方法结合的方式，对西藏地区的暴雨洪涝预测与时空分布进行了研究分析。陈沪生等[13]利用自回归积分移动平均(autoregressive integrated moving average，ARIMA)与小波分析组合的方式，对黄山市的年降水量进行了预测并验证了组合模型的预测效果要优于单一ARIMA模型。

事实上，由于降水量数据周期性与时空性强，关联因素众多，不确定度高等自身原因。传统的统计方法对降水量数据分析具有一定的局限性，特别是这些模型很难对输入序列的保持长期记忆。随着近年来机器学习与深度学习的快速推进，针对降水量等气象要素非线性预测问题有了大幅提升。Hochreiter等[22]提出了长短期记忆神经网络(long short-term memory，LSTM)，解决了循环神经网络(recurrent neural networks，RNN)模型容易出现的梯度消失的问题。Gers等[23]在此的基础上引入了遗忘门机制，使LSTM能够将自身状态重置，减少了网络崩溃的概率。Gers等[24]则在LSTM内部增加了窥视孔(peephole)连接，极大地增强网络对输入序列之间细微特征的区分能力。Graves等[25]在前人的基础上提出了一种双向长短期记忆神经网络(bidirectional long short-term memory，BLSTM)，即当前应用最广泛的一种LSTM模型。现有的中文文献中，通过LSTM进行降水量预测的中文文献较少[26-27]，针对LSTM模型的多种模式的分析应用更是屈指可数，因此如何充分展示LSTM在气象数据预测分析中的多种应用形态是着重考虑的问题。

现首先对原始数据进行分解与周期值提取，其次对季节性自回归积分滑动平均(seasonal autoregressive integrated moving average，SARIMA)模型的构建以及多类别输入输出模式变换的长短期记忆神经网络时间预测模型构建进行系统研究，最后以扬州市区1960—2019年8种基本气象要素为例进行应用与比较分析，对于理解与快速应用深度学习LSTM模型具有借鉴参考价值。

1 资料与方法

1.1 资料来源

数据主要来源为扬州市区气象观测站(编号：58245)1960—2019年8种基本气象要素的21 915行日数据，如表1所示。其中将日数据按照时间序列重采样处理成月数据、季数据和年数据。以月降水量均值数据为例，是指对整个月份中日降水量求和后再除以总天数的均值，将月降水量均值数据简称为月降水量数据。

表1 扬州市区1960—2019年气象要素日数据的样例表

1.2 方法简介

1.2.1 季节自回归积分滑动平均模型(SARIMA)

ARIMA模型全称为自回归积分滑动平均模型，是一种将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型，通俗化可理解为一个未来值与过去的若干个观测值以及高斯白噪声值呈现某种相关性。ARIMA模型建立过程中需要对原数据进行误差趋势季节性法(error trend seasonality，ETS)分解，确定周期与趋势等信息。ARIMA(p,d,q)模型中p为自回归项，q为移动平均项数，d为将非平稳的时间序列转化成平稳时间序列时所做的差分次数，数学表达式为

φ(B)(1-B)dYt=c+θ(B)ε(t)

(1)

式(1)中：Yt为时间序列数据；φ(B)为p阶自回归系数多项式；θ(B)为q阶滑动平均系数多项式；ε(t)为白噪声，服从均值为0且方差为常数的正态分布；d为差分次数；c为常数。

对于具有季节周期性的序列，需要基于季节性将ARIMA模型进行改进，称为SARIMA模型，标记为SARIMA(p,d,q)(P,D,Q)S，表达式为

(2)

(3)

1.2.2 长短期记忆神经网络模型(LSTM)

从网络结构层理解，标准神经网络结构是包含输入层、隐含层、输出层，层与层之间通过权值连接，通过激活函数控制输出。标准的RNN网络结构中，是在隐含层内神经元之间也建立权连接。一般的RNN结构中只有激活函数tanh或者ReLU，如图1(a)所示。该结构的每一个隐含层cell的状态仅有该时刻输入xt和上一个隐含层cell输出决定，这种单线程的模式会导致在长链接模式中出现梯度消失或者梯度爆炸等问题，很难将长时期的依赖关系联系起来。Bengio等[28]以及Hochreiter[29]曾经对这个问题进行过深入的研究，发现RNN的确很难解决这个问题。然而，LSTM神经网络[22]是一种基于RNN网络基础上发展而来的算法，并加入了一种特别的自循环结构，解决了RNN固有的“长期依赖”问题[30-31]。LSTM的链式结构中重复存储块具有4个相互作用的层，如图1(b)所示，包括1个遗忘门(ft)、1个输出门(ot)和1个输入门(it)以及一个记忆单元(C)。LSTM网络通过门控制将加法运算带入网络中，一定程度上解决了梯度消失的问题。sigmoid函数即σ函数。作为逻辑函数，sigmoid函数和tanh函数均为饱和激活函数，在输入值出现极限情况下，输出分别保持在0和1、-1和1，而保证门的开关作用。

符号⊗表示向量的点乘运算；符号⨁表示向量的加法运算；σ表示sigmoid激活函数；tanh表示生成输出时的双曲正切激活函数

输入门：

it=σ(xtUi+ht-1Wi)

(4)

遗忘门：

ft=σ(xtUf+ht-1Wf)

(5)

输出门：

ot=σ(xtUo+ht-1Wo)

(6)

候选记忆细胞：

(7)

单元细胞记忆输出：

(8)

单元隐藏状态输出：

ht=tanh(Ct)ot

(9)

从时间跨度层理解，LSTM增加了一种携带信息跨越多个时间步的方法。假设有一条传送带，其运行方向平行于你所处理的序列，LSTM序列中的信息可以在任意位置跳上传送带，然后被传送到更晚的时间步，并在需要时原封不动地跳回来。LSTM保存了信息以便后面随时使用，同时防止较早期的信号在处理过程中逐渐消失。而RNN仅在短链接中有较好表现，因为随着层数的增加，网络最终变得无法训练。如图2所示，RNN仅在短链接模块可以有效地进行数据的学习。

图2 RNN序列信息传输过程示意图

1.2.3 评价指标

在LSTM神经网络模型给出预测值之后，需要对模型的预测精度以及有效性进行评价，主要采用均方根误差(root mean squared error，RMSE)来衡量这预测值与真实值的偏差程度，其计算方法[32]为

(10)

1.2.4 预测模式

时间序列的预测模型按照预测步数主要分为两种类别，分别是如图3(a)所示的单步(静态)预测(one-step ahead)模型和如图3(b)所示的多步(动态)预测(multi-step ahead)模型。对于单步(静态)预测模式分为两种类型，第一种相邻预测，可以理解为向前外推1个时间步长，比如输入为t时刻状态，输出为t+1时刻状态。这种预测模型能较好地避免误差累积问题，预测效果较好，但是对实际数据依赖性较强。第二种间接预测，可以理解为以现在或者之前若干个时间状态直接外推n个时间步长，比如输入为t、t+1、t+2时刻状态，输出为t+n、t+n+1、t+n+2时刻状态。采取的是多输入多输出的预测模型，从而避免误差累积问题，但是这种情况则需要进行较复杂的数据准备以及大量的数据支持。对于多步(动态)预测模型，事实上它是以单步(静态)预测为基础，将初始的单步(静态)预测的结果作为历史数据，并作为输入值进行下个时刻的预测形成的一种动态模式，这种情况下容易造成误差累积，但是多步(动态)预测模型对于实际数据的依赖性弱，时效性长。

字母上方的“^”符号表示回归值，表示对实际数据的预测；φt为t时刻真实值；为t时刻预测值；Zi、Zj为组合方式；为已训练的预测模型

2 结果与分析

2.1 SARIMA模型预测研究

2.1.1 月降水量数据分解与周期性特征

图4所示为扬州市区近60年间月降水量数据分解图。原始月降水量[图4(a)]均值为2.878 mm，并无明显的周期性规律。原始月降水量分解后的总体趋势[图4(b)]可以看出整体斜率为2.4×10-5，趋于极其微弱的上升稳定状态。周期值[图4(c)]稳定为12个月，与实际的季节变化周期吻合，为后期消除因周期性带来的非稳定问题提供依据。周期[图4(c)]内的月降水量峰值为3.74 mm，由于算法自分解中无数据条件限制，因此会出现月降水量谷值出现负值。针对剩余项[图4(d)]的时间序列值，采用迪基-福勒检验(augmented Dickey-Fuller test，ADF)法检测其稳定性，测试统计值小于临界阈值1%时对应的P值，即检验统计量为-11.595﹤-3.440，因此是平稳序列。采用LB(Ljung-Box)检验白噪声检验，其P值大于置信区间阈值，即13.96%> 5%，因此是无自相关性的白噪声序列，由此，可以判定对原始数据分解无需再进一步深入。

图4 1960—2019年扬州市区月降水量分解图

2.1.2 月降水量数据差分变换

根据月降水量序列差分图(图5)可知，月降水量数据具有周期性(T=12)，并对数据进行12阶差分消除周期性。通过迪基-福勒检验，即检验统计量为-14.619<-3.440，说明12阶差分处理降水数据具有平稳性，且没有单位根。根据12阶差分数据绘制其滚动均值曲线和滚动标准差曲线，从曲线的趋势可以看出数据的整体变化没有明显趋势，处于稳定状态。

图5 1960—2019年扬州市区月降水量序列差分图

2.1.3 降水量数据SARIMA模型预测

在处理后的稳定态数据中截取的2000年1月—2019年12月作为时间数据样例。采用单步静态预测模式与多步动态预测模式进行对比分析，结果如图6所示。对比图6(a)、图6(b)两种模式的结果可知，SARIMA模型动态模式的预测结果基本处于稳定的周期性变化，其月降水量峰值约为6 mm，可以将6(a)中纵坐标蓝色虚线作为参考线。而静态模式的预测结果则出现明显的波动变化，2009年之后多次出现超出蓝色虚线(>6 mm)的峰值数据，对实际月降水量变化趋势的表征效果显示更好。实际上，SARIMA模型中的动态模式下的训练算法是在初始数据输入阶段采用小部分原始历史数据，在得到初始预测数据后，则将预测数据加入历史数据当作实际数据来对下个时间步进行预测，容易造成误差累积，整体趋于周期变化趋势。而SARIMA模型静态模式下，使用实际值作为历史数据，属于直接相邻预测，预测效果更好。

在图6对比分析的基础上，图7左侧显示为通过训练过程(1960—1999年数据)和测试过程(2000—2019年数据)之后建立的SARIMA模型单步静态模式月降水量预测图，图7右侧显示为2020年的月降水量预测的放大图。图7左侧图中绿色区域是单步静态模式下基于测试数据的95%置信区间预测结果，与实际的月降水量变化趋势一致性较好，且变化范围的差距也较小。为了提高对未来降水量预测的准确率，将单步静态模式和多步动态模式二者结合。先使用单步静态模式进行时序模型训练，并使用该模型作为基础，采用多步动态模式进行2020年降水量数据预测。预测结果显示2020年的预测降水量峰值出现在7月为10.88 mm(95%置信上限)，而实际月降水量为双峰值，分别出现在6月为12.39 mm和8月为7.95 mm。预测值与实际值对比结果有一定偏差。

图6 2000—2019年扬州市区月降水量预测模式对比图

图7 基于SARIMA模型的月降水量预测图

2.2 LSTM模型预测研究

2.2.1 数据准备

LSTM神经网络模型的建立中需要对数据进行取样准备。其主要方式是将时间序列中过去的时刻值与现在的时刻值进行一一映射，并依据大量映射样本建立映射函数，从而实现数据预测。因此预测前需要原数据将分解成多个相同长度的数据段，以下列举主要使用的两种取样。假设给定一条需要预测的部分时间序列以及可能该序列相关的另外一条序列如下。

预测序列：1,2,3,4,5,6,7,8,9,10，…。

相关序列：a,b,c,d,e,f,g,h,i,j，…。

假设输入取样中选择的是3个连续的过去时刻的数据作为输入，而当前的1个或者是连续的3个时刻作为输出值，通过神经网络训练学习，可以得出未来时刻的数据预测，如表2所示。

表2 LSTM数据取样示例分类表

2.2.2 LSTM多输入单输入多步动态预测

对于LSTM模型，如果时间步长选择太长，模型训练曲线会出现振荡，算法结构会近似认为前后距离很远的数据还具有联系。如果时间步长选择太短，则会无法充分学习周期性规律。最后都会导致预测准确率降低。因此对于原始数据首先需要进行周期确定。图8(a)～图8(c)分为预测回顾(look_back)为1步、3步和12步时预测结果图。红色曲线与绿色曲线是以训练数据和测试数据为基础，在预测回顾(look_back)为1步、3步和12步时所训练出的模型后，反过来在静态模式下分别对训练数据和测试数据进行预测的结果。图8(a)与图8(b)中红色曲线与绿色曲线的上下限波动幅度显示较弱，图8(c)中红色与绿色曲线的上下限波动幅度显示较强，这是因为图8(c)中预测回顾(look_back)选取的时间步数为12，这也恰好是原始降水量数据的周期值。由此，充分说明了完整的周期数据能更好地表达时序序列内部的信息。图8(d)是不同预测回顾步数情况下的预测曲线合集。

图8 基于LSTM多输入单输入多步动态模式的月降水量预测图

当回顾步数为1步、3步或者低于1个周期(12步)，模型并不能充分学习并提取数据中的信息。同样，当回顾步数超出1个周期(12步)时，模型也可能会造成部分数据信息的过度提取，导致后期模型应用中出现过拟合。

如图9所示，为不同回顾步数和cell单元维数下训练集的损失值(loss)与测试集的损失值(val_loss)的变化曲线。通过图9(a)、图9(b)对比显示，预测回顾步数为3步或者12步的预测显示曲线变化趋势一致，不断稳定下降，说明神经网络状态良好仍有学习空间。通过图9(c)、图9(d)显示，当预测回顾步数为24步时，调整不同神经元个数(Num_units)(即神经网络维度)后，训练集仍然下降，测试集出现缓慢上升的趋势，说明可能模型发生过拟合。虽然算法中使用了Dropout等丢弃函数，对神经网络训练时每次输入数据的时候以一定的概率删除部分神经元，但是效果不明显，并不能很好地处理过拟合问题。综上，可以通过未来预测曲线和损失函数综合对比，对数据的预测回顾(周期性)进行较好的确定。

图9 LSTM多输入单输入多步动态预测模式下不同回顾步数损失曲线图

2.2.3 LSTM 多输入多输出单步静态预测

LSTM模型中除了图8所示的多输入单输入多步动态预测，也可以通过调整输入数据结构为多输入多输出，进行未来12个月数据的静态预测。结果如图10所示，月降水量预测峰值在9月份为3.58 mm。损失函数训练集的损失值(loss)曲线为减少的趋势，测试集的损失值(val_loss)曲线减少的速率相对较缓慢。

图10 基于LSTM多输入多输出单步静态模式的月降水量预测图

2.2.4 M-LSTM多输入单输出单步静态预测

大多数的多源序列数据之间由于不同介质的存在往往都呈现网状分布，正是这种分布使得多源序列数据之间会有错综复杂的关系。研究表明多种气象数据之间也是一种网状关系，即要素之间是有着一定的联系。使用不同类别的时间序列数据对一种数据进行预测，数据源涵盖的基本气象要素主要分为8种，将8种气要素作为输入量，降水量作为输出量进行数据样例准备，如图11所示。通过M-LSTM(即一种含有不同类型的多输入LSTM)神经网络数据训练与建模，对未来12个月降水量进行预测。

图11 M-LSTM多输入单输出建模过程示意图

M-LSTM无法进行多输入单输出动态预测循环。因为模型的输入数据除了月降水量，还需要其他的变量，但是每次预测输出的数据只有降水量，无法构成前后闭合的输出输入循环。因此，本节使用M-LSTM多输入多输出静态模式对未来12个月的降水量进行预测。如图12(a)所示，训练数据预测曲线(红色)与测试数据预测(绿色)曲线波动均不明显，整体变化趋近于周期性稳定变化。对与未来12个月降水量的预测，峰值出现在8月，数值为3.83 mm。如图12(b)所示，训练损失函数和测试损失函数的变化较为一致，M-LSTM神经网络状态良好，且仍具有一定的学习空间。

图12 基于M-LSTM多输入单输出单步静态模式预测图

2.3 算法评估

将以上多种预测模型的不同模式下的评价情况，按照均方根误差(RMSE)的计算结果进行展示，如表3所示，在LSTM多输入单输出动态预测模式下，预测回顾(look_back)选取12步时，其RMSE测试数据取值最优。

表3 RMSE结果统计表

3 结论

以扬州市区1960—2019年气象要素为例，从原始数据的分解与周期性提取，到季节性的自回归积分滑动平均模型(SARIMA)构建以及多类别输入输出模式变换的长短期记忆神经网络时间预测模型(LSTM)构建进行了系统研究与比较分析，主要结论如下。

(1)传统的SARIMA模型中静态模式预测结果较动态模式能更好地反映出降水量数据变化趋势，且与实际值差距较小。动态模式容易造成误差累积或整体易呈现周期性稳态变化，实时性欠缺。

(2)深度学习LSTM 多输入单输出动态预测模式下，完整周期的数据输入可以让神经网络更好地学习数据的变化规律，体现了确定周期值的重要性。然而将多个周期数据作为一个输入单位，易造成模型过拟合。该LSTM模型(look_back=12)对扬州市区月降水量预测准确度优于传统的SARIMA模型，RMSE训练值低0.02，且仍然有继续调参提升的空间。

(3)LSTM多输入单输出动态模式(look_back=12)较LSTM多输入多输出静态模式，RMSE测试值可低0.33，体现出该模式对扬州市区月降水量预测准确度更高。

(4)M-LSTM多输入多输出静态模式预测准确度优于LSTM多输入多输出静态模式，RMSE训练值低0.62，RMSE测试值低0.19，反映出M-LSTM多输入多输出静态模式的优点，但该模式对数据输入的参数结构要求较高。