1. 引言
在全球气候变暖的影响下,降水普遍发生变化,深入研究区域降水演变规律,对保障社会经济发展具有重要意义[1]。降雨是有效管理水资源、预防和减少洪水以及城市拥堵问题预警的重要参数之一,因为降雨是调节城市气候的重要因素,有效避免了城市高温天气的发生。雨水可以降低空气中的灰尘和其他污染物,使空气变得清新,有利于居民呼吸健康。合理利用雨水,可以有效缓解城市用水紧张的问题。降雨也会对农业生产造成影响,适时的降雨能够满足农作物的生长需求,增加作物产量,但降雨过多导致地面积水、洪涝灾害以及庄稼受淹或死亡[2]。降雨过多往往会导致道路湿滑、积水,增加交通事故的风险。降雨过多也会导致交通拥堵,对城市居民的出行造成一定的不便,增加了城市居民的出行时间和成本。
有效预测城市降雨不仅可以有效减少降雨、干旱等自然灾害造成的巨大经济损失,而且对于工农业的发展、城市居民的出行等也具有十分重大的实际意义。刘杨,徐淑琴等建立了一个时间序列模型,用于预测1973年至2008年查哈扬农场作物生育期的月降雨量。结果表明,该模型能够显示该地区作物生育期的月降雨量变化规律[3]。程敏等利用arima模型对济南市2010~2015年降雨量进行了预测,预测结果显示5年的平均降雨量接近过去57年的平均水平,这表明济南市在未来五年出现干旱和洪涝灾害的可能性较小[4]。白玉洁提出了一种将小波变换与时间序列预测模型(arima)相结合的降雨预测方法,以提高降雨预测的准确性。结果表明改进后的方法比传统预测方法具有更高的降雨量预测精度,较好地反映了降雨量的变化规律。这为降雨预测提供了一种新的预测方法[5]。陈凯,黄丹等则是考虑arima模型和gm(1,1)模型的组合模型对沥青路面使用性能进行预测,发现组合预测模型表现出了一定的优势,其预测结果可为养护计划决策提供一定的科学依据[6]。应瑶,冯利华等通过对金华市梅汛期降水基本特征进行分析,并运用平稳时间序列的线性外推法建立各个地区降水量的预测模型,结果表明金华地区梅汛期降水量的年际变化呈多峰型,梅雨期降水量的预测较准确[7]。刘永和马飞等建立基于秩次集对分析的年降雨量预测模型,进而对德州市某一年的降雨量进行预测,发现使用该方法对降雨量进行预测结果是可行有效的[8]。
上述研究考虑了对降雨量数据直接建立模型进行预测,或者对数据做一定处理后再建模进行预测。但是降雨量数据作为时间序列数据,季节性趋势并没有考虑过。本文以西宁市2009~2021年降雨量为研究对象,首先分析降雨量的季节性趋势,确定季节性趋势之后建立我们所需要的季节性arima模型。还引入了lstm模型,考虑了降雨量的季节性特点,通过对数据进行适当的预处理,设定模型的参数,建立模型。最后使用所建立的模型做预测分析。通过使用lstm模型与季节性arima模型的预测结果与实测值的比较,进一步说明这两种模型在实际应用中的性能。
2. 基础知识
季节性arima模型是时间序列分析方法之一,考虑了时间序列数据中存在季节性变化周期变化的各种非平稳时间序列数据。这些时间序列数据的值随时间变化,并且经差分后过去的观测值与未来值以及对应季节的数据之间存在较强的相关性。由于季节性arima模型能够很好地捕捉季节性因素对时间序列的影响,因此适用于季节性变化对预测结果有重要影响的情况。在处理具有季节性变化的时间序列数据时,季节性arima模型通过季节性差分来消除季节性因素的影响,使得时间序列数据在差分后呈现平稳性。然后,通过自回归和移动平均来捕捉数据的趋势和随机性。这种数据处理方式使得它能够适应许多不同类型的时间序列数据,并能够捕捉数据中的长期和短期变化。
arima模型即差分整合移动平均自回归模型,时间序列预测分析方法之一。arima模型是是由自回归模型(ar)、差分模型(i)和移动平均模型(ma)三部分组成的[9]。arima模型的基本思想是将非平稳的时间序列数据转化为平稳时间序列数据,然后再对转化后的平稳时间序列进行建模和预测。arima模型的建立需要确定三个参数:p、d、q,分别代表自回归项数、差分次数和移动平均项数。其中,p和q是可以通过自相关图和偏自相关图来确定的,而d则需要通过观察时间序列数据的平稳性来确定。
季节性arima模型即arima(p,d,q) (p,d,q)s模型,该模型共有7个参数,分别是自回归阶数(p)、差分次数(d)、移动平均阶数(q)、季节性自回归阶数(p)、季节性差分次数(d)、季节性移动平均阶数(q)、周期长度(s)。p值从偏自相关图中推断,q值从自相关图中推断,d值一般为0或1。
p阶自回归模型可以表示为:
(1)
式(1)中的p为非负的整数,
为均值为0,方差为
的白噪声序列。
为自相关系数,
为当前值
是过去时刻
到
的值。
q阶移动平均模型可以用下式表示:
(2)
式(2)中
为常数,
为均值为0,方差为
的白噪声序列。
为移动平均模型的系数。
arima
模型可以表示为:
(3)
式(3)中
到
是ar模型的参数,描述的是当前值与过去p个值之间的关系。
到
是ma模型的参数,这些参数用来描述当前值与过去q个时间点的误差之间的关系。
是在t时间点的误差项。
是常数项,
是我们的平稳时间序列数据。
季节性arima模型可以表示为:
(4)
其中
为延迟算子,
表示季节性延迟算子,延迟多项式算子
和
,季节性延迟多项式算子
和
分别表示为:
(5)
(6)
(7)
(8)
简单来说就是对时间序列
做d次季节性差分和d次差分后得到一个新的序列
,然后对
建立季节性arima模型。
自相关函数acf (autocorrelation function),时间序列数据
与
的相关系数记为
。
称作
的间隔为k的自相关系数,计算公式为:
(9)
其中k代表滞后期数,自相关系数
组成的集合称为
的自相关函数,自相关函数acf描述的是时间序列观测值与其过去的观测值之间的线性相关性。
偏自相关函数pacf (partial autocorrelation function),偏自相关函数描述的是在给定中间观测值的条件下,时间序列观测值预期过去的观测值之间的线性相关性,对于平稳时间序列
所谓滞后k阶偏自相关系数是指在给定中间条件下,或者说在删除
个随机变量
的干扰之后,
与
之间影响的相关程度。用公式表示为:
(10)
对于arima模型参数
的确定可以通过对平稳时间序列的acf图和pacf图分析得到,根据差分的次数得到参数d的值。
当样本容量很大时,在aic准则中拟合误差提供的信息就要受到样本容量的放大,而参数个数的惩罚因子却和样本容量没关系(一直是2),因此当样本容量很大时,使用aic准则选择的模型不收敛与真实模型,它通常比真实模型所含的未知参数个数要多。bic (bayesian information criterion)贝叶斯信息准则弥补了aic的不足,计算公式为:
,其中n是样本容量。
均方根误差(root mean squared error),简称rmse,它衡量了预测值与真实值的偏差。rmse的值越大,就说明模型的预测误差越大,模型的拟合效果比较不好。相反rmse的值越小,说明模型的预测误差越小,拟合的也更好。平均绝对百分比误差(mean absolute percentage error),称为ape,是一种衡量预测值准确性的标准。mape的值越小模型的预测精度就越好,相反mape的值越大,说明模型预测性就比较差。
lstm是一种特殊的循环神经网络(rnn),传统的rnn在处理长序列数据时,往往面临梯度消失和梯度爆炸的问题,导致难以有效地捕获长期依赖关系[10]。lstm通过大量的数据训练,能够学习到数据中的普遍规律,并具有一定的泛化能力,能够在一定程度上适应未见过的数据。lstm通过引入门控机制和细胞状态,显著提升了rnn处理长序列数据的能力,使得lstm能够有效地捕捉时间序列数据中的长期依赖关系。lstm的核心在于其独特的网络结构和门控机制。它包含三个门:输入门、遗忘门和输出门,以及一个细胞状态。这些门和细胞状态共同决定了信息在网络中的流动方式。lstm的结构图如图1。
figure 1. lstm structure diagram
图1. lstm结构图
遗忘门:遗忘门负责决定从细胞状态中丢弃哪些信息。将上一时刻的输出
与当前时刻的输入
结合,并通过sigmoid函数计算得到一个阈值为[0, 1]的张量
,该值决定了细胞状态中信息的保留程度,见图2。张量
,其中
是权重,
是偏置,
为激活函数。
figure 2. oblivion gate
图2. 遗忘门
输入门:输入门决定了当前时间步的输入信息有多少应该被保留在细胞状态中。tanh函数会给出一个新的候选向量
,输入门为
中的每一项产生一个在[0, 1]之间的值
,控制新信息被加入的多少,见图3。
,
,其中
是输入门的权重矩阵,
是偏置项,输入的
和
通过sigmoid函数
,将输入值映射到0到1之间。
接近1时,新的输入会被完全保留,当
接近0时,新的输入会被丢掉。
细胞状态:细胞状态是lstm中的关键组件,它承载着序列中的长期依赖信息。在每个时间步,细胞状态都会根据输入门、遗忘门和上一个时间步的细胞状态进行更新,见图4。
,其中将得到
的与输入门的输出值
相乘,并与经过遗忘门
处理的单元状态(
)相加,作为下一个时间步的细胞状态
。
figure 3. input gate
图3. 输入门
figure 4. cell state
图4. 细胞状态
输出门:输出门决定了当前时间步的输出值。其中
,
它基于当前时间步的细胞状态,通过sigmoid函数计算一个输出门值,并使用tanh函数对细胞状态进行缩放,然后将两者相乘得到最终的输出,见图5。
figure 5. output gate
图5. 输出门
3. 研究过程
3.1. 西宁市简介
西宁市,古称青唐城、夏都等,西宁的历史悠久文化底蕴丰厚,还有着丰富的自然资源。西宁是青海省省会,也是国务院批准的中国西北部重要的中心城市。西宁位于青藏高原东北部,地处湟水中游河谷盆地,被誉为青藏高原的东方门户。西宁市下辖5个区和2个县,总面积达到7660平方千米。西宁更是一座旅游胜地,凉爽的环境,舒适的公园为市民提供了休闲放松的去处。每到夏日来临,蓝天白云,绿草如茵,更加让人心旷神怡。西宁市在夏季阳光日照强降雨量较为集中,冬季则相反。
3.2. 数据来源
本文中研究的降雨量数据来源于青海省统计年鉴,收集并汇总了2009年~2021年的月降雨量数据,见表1。
table 1. monthly rainfall of xining city from 2009 to 2021 (unit: mm)
表1. 西宁市2009年到2021年的月降雨量(单位:mm)
时间 |
一月 |
二月 |
三月 |
四月 |
五月 |
六月 |
七月 |
八月 |
九月 |
十月 |
十一月 |
十二月 |
2009 |
4.7 |
0.4 |
15.7 |
18.4 |
66.4 |
66.7 |
68.2 |
106.5 |
87.3 |
16.9 |
7.1 |
0.8 |
2010 |
1.2 |
1.2 |
4.2 |
17.6 |
99 |
50.8 |
58.2 |
65.4 |
72.6 |
33.2 |
0.4 |
1.2 |
2011 |
1 |
0.4 |
2.1 |
7.1 |
49.5 |
63.4 |
54.5 |
85.6 |
93.1 |
21.1 |
11.5 |
1.1 |
2012 |
3 |
2.2 |
8.1 |
6.9 |
76.4 |
58.6 |
147.7 |
61.8 |
57.4 |
20.2 |
3.1 |
0.7 |
2013 |
0 |
1.3 |
0.9 |
18.6 |
67.6 |
68.1 |
74.5 |
107.9 |
54.9 |
10.5 |
6.7 |
2.6 |
2014 |
0 |
2.3 |
1.9 |
48 |
26.1 |
109.3 |
86.4 |
59.3 |
71.6 |
32.2 |
9.4 |
0 |
2015 |
1.1 |
0.6 |
3.8 |
15.4 |
21.7 |
51.8 |
71.9 |
61.1 |
45.7 |
25.9 |
5.9 |
1.3 |
2016 |
0.4 |
2.5 |
28.1 |
16.4 |
59.9 |
41.3 |
84.8 |
67.5 |
81.3 |
59.9 |
0.4 |
1.6 |
2017 |
0 |
2.4 |
9.9 |
29 |
87 |
40.6 |
50.1 |
120.8 |
90.4 |
33 |
0 |
0.8 |
2018 |
2.3 |
0.1 |
2.1 |
30.2 |
37.5 |
55.8 |
113 |
136.6 |
80.5 |
31.4 |
29 |
0.6 |
2019 |
2.4 |
2.4 |
4.1 |
22 |
52.4 |
95.8 |
70.8 |
129.2 |
125.3 |
29.4 |
2.5 |
1.1 |
2020 |
3.8 |
2.2 |
7.1 |
7.5 |
67.8 |
89.9 |
87.7 |
163 |
68.1 |
3.6 |
3.5 |
2.2 |
2021 |
1.6 |
3.5 |
2.4 |
38.1 |
49.6 |
49.6 |
62.1 |
89 |
126.7 |
28.3 |
3.3 |
0 |
根据表1的降雨量数据,使用r软件画出西宁市2009年到2021年的年降雨量时序图,如图6。
由图6可以看到在夏季的时候降雨量明显增加,在春季和冬季降雨量大幅度减少。数据不是很平稳,降雨量有明显的季节性趋势,因此对数据做季节性差分处理减少趋势性的影响。结果如图7所示。在对数据做了一阶差分和一阶季节性差分后数据表现得更加平稳,处理过的数据满足我们建立季节性arima模型的建模要求。
figure 6. monthly rainfall time series of xining city from 2009 to 2021
图6. 西宁市2009年到2021年的月降雨量时序图
figure 7. monthly rainfall time series of xining city from 2009 to 2021 after difference
图7. 西宁市2009年到2021年的月降雨量差分后时序图
3.3. 季节性arima模型参数的估计
首先我们把西宁市2009年~2021年的数据作为我们的原始数据,对原始数据做季节性差分处理,得到平稳化的序列,以及平稳化序列的acf图和pacf图,见图3。
而对于模型中参数的选择则通过图8中自相关函数图和偏自相关函数图来确定。从acf图中发现在滞后阶数为1的时候acf值快速接近0,由此可以确定q取1较为合适。在滞后阶数为12的时候acf值又超出置信区间表现出季节效应,可以考虑季节移动平均的阶数取1。从pacf图中发现在滞后阶数为2后的pacf值大部分都是在置信区间或靠近置信区间的,可以考虑建立arima(2,1,1)(0,1,1)模型。此时建立的模型中有些模型参数并没有确定,且根据acf图和pacf图确定的阶数并不一定能保证模型拟合度较好。因此我们可以根据不同的参数构建模型,通过多次建模并根据r方、rmse、mape和bic的值来选择出拟合度比较好的模型,最后用选出的模型进行预测。多次建模后不同模型的拟合值如表2所示。
figure 8. acf diagram and pacf diagram of the stationary sequence
图8. 平稳化序列的acf图和pacf图
table 2. model fitting values under different parameters
表2. 不同参数下的模型拟合值
季节性arima模型 |
r方 |
mape |
rmse |
bic值 |
arima(1,1,0)(0,1,1) |
0.599 |
229.943 |
24.963 |
6.539 |
arima(1,1,1)(0,1,1) |
0.739 |
152.145 |
20.189 |
6.149 |
arima(1,1,2)(0,1,1) |
0.739 |
144.465 |
20.263 |
6.191 |
arima(2,1,0)(0,1,1) |
0.656 |
230.459 |
23.211 |
6.428 |
arima(2,1,1)(0,1,1) |
0.742 |
157.522 |
20.177 |
6.183 |
arima(2,1,2)(0,1,1) |
0.741 |
144.644 |
20.264 |
6.226 |
通过表2发现r方最大的模型为arima(2,1,1)(0,1,1),r方的值为0.742说明这时候的模型是比较好的。但是也要考虑到过拟合的情况,因此我们再次比较一下其他值。bic值较小的模型拟合优度会好一些但在比较bic值时发现bic值最小的模型为arima(1,1,1)(0,1,1)。rmse值最小的模型是arima(2,1,1)(0,1,1),rmse值为20.177,这也表示该模型的预测能力较好。因此我们选择arima(2,1,1)(0,1,1)作为我们降雨量预测的最终模型。
3.4. lstm模型的设置
lstm模型可以用更深层的神经模型训练长期序列数据,目的是用上一个月的数据来预测下一个月的数据。由于数据是时序数据,在划分训练集和测试集时并不对数据进行打乱,数据集的前80%作为我们的训练集,之后的20%作为测试即。之后对划分后的数据做标准化处理。在lstm模型中选择mse作为损失函数,优化算法选择adam,它结合了momentum和rmsprop算法的思想,以更有效地在训练过程中更新模型的权重。该模型包含一个lstm层有4个单元,一个激活层,以及一个输出层,激活函数默认是线性激活函数。通过选择不同的迭代值来得到损失相对平稳的模型。迭代值选择为200。使用训练数据训练lstm模型,可以通过计算每个验证集上的损失来检测模型的性能。发现在迭代次数为120次的时候损失就相对平稳。在训练完成后,使用测试集评估模型的性能。在此画出了在训练集上的预测结果和实际值的对比情况,见图9。
figure 9. actual and predicted values of the training set
图9. 训练集实际值和预测值
从图9可以观察到预测值代表的红线是比较贴合实际值,但也发现在详细分析预测结果时,我们发现预测值与实际值之间存在一些差异。在对比实际值与lstm模型预测值的结果时,我们发现尽管模型在某些方面能够较好地拟合数据的趋势,但预测值存在明显的滞后性。可能是由于模型在训练过程中未能完全捕捉到数据中的快速变化或突发性事件。但是lstm模型还是能够较好的预测数据。
3.5. 预测结果
经过对数据的处理,模型参数的选择最后我们选择arima(2,1,1)(0,1,1)和lstm模型来预测西宁市2022年的月降雨量,结果如表3所示。由表3我们发现在2022年的前三个月,月降雨量都比较小,在三月的时候甚至没有下雨,这对于我们的预测工作具有很大的挑战性。在用arima(2,1,1)(0,1,1)模型进行预测后发现在前四个月的预测误差会小一些,第五个月之后的预测误差就变大了。而lstm对降雨量预测时前四个月的预测误差相对稳定,并且在未来十二个月的预测中误差都相对稳定。我们观察表3中的预测值时发现在1月的时候值比较小,随着时间的移动在8月降雨量达到最大,之后再逐渐变小。季节性arima在预测步数为4或者小于4的时候预测精度是可以的,误差比较小。同时观测到lstm对降雨量预测时,预测误差略大一些,但结果还是相对符合的。对于一些具有趋势性的数据预测结果是比较准确的,当数据有较大浮动时往往会对预测结果造成影响。2022年降雨量预测值和实际值的对比结果见图10。
table 3. predicted and actual rainfall in 2022 (unit: mm)
表3. 2022年降雨量预测值和实际值(单位:mm)
模型 |
arima(2,1,1)(0,1,1) |
lstm |
时间 |
真实值 |
预测值 |
预测误差 |
预测值 |
预测误差 |
2022年1月 |
1.1 |
3.04 |
1.91 |
11.03 |
9.93 |
2022年2月 |
4.3 |
5.87 |
1.57 |
11.77 |
7.47 |
2022年3月 |
0.0 |
8.51 |
8.51 |
13.92 |
13.92 |
2022年4月 |
35.8 |
29.79 |
−6.01 |
11.04 |
−24.76 |
2022年5月 |
11.4 |
62.71 |
51.31 |
36.09 |
24.69 |
2022年6月 |
41.4 |
72.65 |
31.25 |
18.79 |
−22.61 |
2022年7月 |
94.5 |
84.71 |
9.79 |
40.06 |
−55.44 |
2022年8月 |
250 |
115.26 |
−134.74 |
73.67 |
−176.33 |
2022年9月 |
54.8 |
97.18 |
42.38 |
121.62 |
66.82 |
2022年10月 |
18.9 |
33.08 |
14.18 |
49.34 |
30.44 |
2022年11月 |
2.3 |
13.30 |
11 |
24.06 |
21.76 |
2022年12月 |
0.9 |
8.16 |
7.26 |
12.57 |
11.67 |
从图10中可以发现在前四个月的实际值与预测值非常接近,第五个月预测值与实际值之间的差距逐渐变大,预测效果变差。图10也表明所建立的季节性arima模型适合做西宁市降雨量的短期预测,预测步长小于等于4时效果较好。而在预测步长增加时或者数据没有快速变化或突发性事件使用lstm也是一个好的选择,尤其是在预测长期降雨量时lstm具有显著优势。
figure 10. predicted and actual rainfall values of xining in 2022
图10. 西宁市2022年降雨量预测值和实际值
4. 结语
季节性arima模型是一种时间序列预测模型。本文考虑了研究数据具有季节性趋势的情况,通过对降雨量数据做季节性差分处理,得到平稳后的时间序列。之后通过建立季节性arima模型arima(2,1,1)(0,1,1)和lstm模型对西宁市2022年的月降雨量进行预测,结果表明在短期预测时所建立模型的预测精度比较好,2022年前四个月的预测误差都在10 mm之内,结果可以接受。可以用季节性arima模型做西宁市的短期降雨量预测,能够较好地反应西宁市的降雨规律。
lstm模型也展现出了在降雨量预测中的强大潜力。通过其独特的门控机制,lstm模型能够学习并捕捉到降雨量数据中的长期依赖性和复杂的非线性模式。这使得lstm模型在处理具有复杂变化趋势的降雨量数据时表现优异。无论是选择季节性arima模型还是lstm模型,我们都能够获得较为准确和可靠的降雨量预测结果。通过对降雨量进行预测能够帮助我们更加合理地利用水资源,提高水资源的利用率。
季节arima模型主要依赖于历史数据的线性关系进行预测,当时间序列数据中的趋势或季节性模式发生显著变化时,模型的长期预测精度可能会下降。lstm模型由于其复杂的结构和大量的参数,需要较长的训练时间和较高的计算资源。未来可以考虑采用更高效的优化算法来优化模型的拟合过程,还可以探索季节arima模型和lstm模型的混合使用方式,充分发挥各自的优势,提高时间序列预测的准确性和效率。