基于季節性自回歸積分滑動平均與深度學習長短期記憶神經網絡的降水量預測

張麗婷，李鵬飛，龐文靜，惠雯，秦孟晟

(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多輸入多輸出靜態模式的優點，但該模式對數據輸入的參數結構要求較高。