ARIMA-SVR組合模型在基于標準化降水指數干旱預測中的應用

許德合，張 棋，黃會平

(1.華北水利水電大學地球科學與工程學院，河南鄭州450000;2.華北水利水電大學測繪與地理信息學院，河南鄭州450000)

旱災被認為是世界上最嚴重的自然災害類型之一［1］，極大程度地影響了人們的日常生活以及農業(yè)產量［2］。干旱是指水分收支或供求不平衡形成的水分短缺現象［3］，全球氣候變暖、碳排放量超標等問題將加劇未來農業(yè)干旱情況，嚴重威脅全球糧食生產，因此準確評估、監(jiān)測、分析干旱情況一直是國內外學者的熱門話題［4］。對干旱進行量化研究有助于研究干旱時空變化特征，提升干旱監(jiān)測能力，開展干旱預報工作，尋求干旱治理及應對策略，對未來我國農業(yè)生產以及防旱抗旱等方面具有重要意義［5］。我們通常選用便于計算的干旱指標來監(jiān)測評估干旱發(fā)生的強度、持續(xù)時間和受災范圍［6］。由于干旱指標種類多、運用范圍廣，且不同專業(yè)和學科對干旱理解不同，因此出現了多種干旱指標［7］。標準化降水蒸散指數(Standard Precipitation Evaporation Index，SPEI)、帕默爾干旱指數(Plamer Drought Severity Index，PDSI)、降雨Z指數(ZIndex)、標準化降水指數 (Standard Precipitation Index，SPI)、綜合干旱指數(Colligation Drought Index，CI)等在氣象干旱、農業(yè)干旱、水文干旱等領域已得到廣泛應用［8-12］。其中，SPI是用于表征某時段降水量出現概率多少的指標，計算結果對干旱分級精度相對較高，所需源數據少(僅利用日或月降水量數據就可進行計算)，適用范圍廣，并且不同時間尺度的SPI值可以適用于不同類型的干旱。由于該指標使用靈活，易于計算，已成為實際應用最廣泛且適用于所有氣候狀況的干旱指標［13-16］。

加強干旱預測方面的研究對相關部門預防干旱災害、減少干旱損失具有重要意義［17］。常用來預測干旱的模型有很多，如人工神經網絡(Artificial Neural Network，ANN)、支持向量機(Support Vector Machine，SVM)和差分自回歸移動平均模型(Autoregressive Integrated Moving Average，ARIMA)等。其中ARIMA模型是時間序列中常用的模型，通常用來預測線性數據;SVM模型是一種二分類模型，通常用來處理非線性數據;支持向量回歸機(Support Vector Regression，SVR)是SVM的一種拓展，多用來進行非線性數據的回歸預測，而ARIMA與SVR組合模型分別對線性模型及非線性模型處理具有優(yōu)勢，所以它們之間存在優(yōu)勢互補［18］。

河南省位于中國中部和黃河中、下游，是中國重要的糧食作物產區(qū)和農業(yè)大省，對于保障國家糧食安全發(fā)揮著至關重要的作用［19］。本研究以河南省國家級氣象觀測站鄭州站為例，計算不同時間尺度的SPI值，利用ARIMA模型及ARIMA-SVR組合模型對其進行預測，并采用均方根誤差(Root Mean Squared Error，RMSE)和平均絕對百分比誤差(Mean Absolute Percentage Error，MAPE)對 2種模型預測能力進行分析。

1 研究方法

1.1 SPI的計算

降水量通常是一種偏態(tài)分布，常采用Γ分布概率描述降水量的變化，再將偏態(tài)概率分布進行正態(tài)標準化處理，最后用標準化降水累計頻率分布劃分干旱等級［20-22］。SPI可以定量化研究多時間尺度的降水量不足。

SPI指數公式為［23-24］:

式中，Y(x)2為與Γ函數相關的降水量分布概率;x為樣本值(即降水量);G為正負系數;u0、u1、u2和l1、l2、l3為常數:

當Y(x)＞0.5時，G=1，當Y(x)≤0.5時，G=-1。

Y(x)由Γ函數求得，其中Γ為概率密度積分公式:

式中，γ，β為Γ分布函數的形狀和尺度參數。

1.2 ARIMA模型

ARIMA分為自回歸模型(Autoregressive model，AR)、滑動平均模型(Moving average model，MA)以及自回歸移動平均模型(Autoregressive Moving Average Model，ARMA)，是傳統(tǒng)的時間序列預測模型。其建模流程是首先判斷模型平穩(wěn)程度，其次利用差分法對非平穩(wěn)時間序列進行平穩(wěn)化處理，然后選取AR(p)，MA(q)對模型進行定階，差分次數記為d，ARIMA(p，d，q)模型就是經過了d階差分后的ARMA(p，q)模型。如下式所示［25］:

式中，φ(L)和θ(L)分別為:

ARIMA(p，d，q)模型的一般式:

上式中，d為差分次數，Δ=1-L，(φ1，φ2，…，φp)為自回歸系數，p為自回歸階數，(θ1，θ2，…，θp)為移動平均系數;ut為白噪聲序列(服從0均值、正態(tài)分布且相互獨立的白噪聲序列)。

其中p、q階數采用赤池信息準則(Akaike Information Criterion，AIC)和貝葉斯信息準則(Bayesian Information Criterion，BIC)來確定，當樣本數N固定時，選擇AIC和BIC最小值來確定p，q。公式如下:

1.3 SVR模型

SVR模型是SVM的推廣，SVR的本質屬性不再是原來的二分類方法，而是回歸方法。由于SVR模型在對非線性數據預測方面具有優(yōu)勢，因此采用ARIMA模型對線性數據SPI進行預測，將所得殘差(非線性數據)傳入SVR模型，再利用SVR模型對殘差進行預測。通過引入徑向基核函數(RBF)來把訓練樣本映射到高維空間下來進行回歸預測，再將回歸問題轉為優(yōu)化問題［26-28］:

式中，w為權重系數為松弛變量，規(guī)定了模型的誤差要求;C為懲罰參數，C越大則支持向量的決策邊界越大;b為偏置項。RBF公式中σ和γ關系如下:

1.4 ARIMA-SVR組合模型

由于ARIMA模型和SVR模型在線性和非線性預測中各有優(yōu)點，因此本文分別采用 ARIMA與SVR模型的優(yōu)點建立組合模型ARIMA-SVR，假設時間序列Yt可視為線性自相關部分Lt與非線性殘差Nt兩部分的組合，利用ARIMA模型對SPI值進行預測，將結果與實際值相減得到殘差，將殘差記為非線性部分，帶入SVR模型進行預測，最后把兩者預測結果相加得到組合結果，即:

1.5 評價驗證指標

在回歸模型當中，平均絕對誤差(Mean Absolute Error，MAE)、均方誤差(Mean Square Error，MSE)、RMSE和MAPE是常見的回歸預測評估指標，其中MAE和MSE是最基礎的評估方法，RMSE和MAPE是回歸任務最常用的性能度量，是前兩種指標的擴展，不但在多場景下可以使用，而且相對前兩種更加準確，因此本文采用RMSE和MAEP作為模型評價的指標。

RMSE均方根誤差

MAPE平均絕對百分誤差

式中，xi是觀測值，yi是預測值，N是樣本數。RMSE和MAPE越接近0，表示預測值與觀測值越接近。

2 實例應用

2.1 數據資料及SPI計算

本文選用1951—2017年河南省國家級氣象站鄭州站逐日降水量數據來進行計算，原始數據來源于國家氣象信息中心提供的中國地面氣候資料日值數據集。運用Matlab數學建模軟件編寫SPI計算程序，分別計算了 1951—2017年的 1、3、6、12個月共4 個尺度的SPI值，記為 SPI1，SPI3，SPI6，SPI12，并通過國家標準氣象干旱等級(GB/T20481-2006)規(guī)定的干旱分級標準(表1)來表征干旱情況［29］。

2.2 時間序列生成及建模流程

本文使用Python 3.6平臺對ARIMA進行建模，并利用Python中matplotlib可視化庫對多尺度SPI計算結果進行可視化展示，如圖1所示。

2.2.1 平穩(wěn)化處理及ARIMA模型定階 由于ARIMA是經過d次差分的平穩(wěn)時間序列ARMA模型，且通常針對平穩(wěn)時間序列進行建模，因此在建模前首先應對時間序列的平穩(wěn)性進行判斷，本文采用觀察時間序列的時序圖和單位根檢驗(Augmented Dickey-Fuller Test，ADF)來判斷平穩(wěn)性。由圖1可得，SPI1、SPI3、SPI6序列無明顯的上升和下降趨勢，SPI12序列略有上升趨勢，進一步對SPI1、SPI3、SPI6和 SPI12進行 ADF檢驗，在 ADF檢驗中，原假設為非平穩(wěn)時間序列且存在單位根，給定顯著水平α=0.05，如果檢驗統(tǒng)計量對應的概率值P＜0.05，則拒絕原假設。其中SPI1、SPI3、SPI6 的ADF檢驗P值均小于0.05，SPI12的P值為0.36753且顯著大于0.05，檢驗結果見表2，因此判斷SPI12序列為非平穩(wěn)時間序列(非平穩(wěn)時間序列一定不是白噪聲序列)，SPI1、SPI3、SPI6為平穩(wěn)時間序列。利用差分法對SPI12非平穩(wěn)時間序列進行平穩(wěn)化處理，經一階差分后時間序列趨于平穩(wěn)(圖2)，再對一階差分后的時間序列進行單位根檢驗(結果見表3)。表3中P值為0.00128，小于0.05，因此差分之后為平穩(wěn)時間序列。采用純隨機性檢驗(Ljung-Box Test)進行白噪聲檢驗，得P值為0.006025，遠小于0.05，因此經一階差分后該序列為平穩(wěn)非白噪聲序列。然后用自相關函數(Autocorrelation Function，ACF)及偏自相關函數(Partial Autocorrelation Function，PACF)來為 ARMA模型定階(圖3)。這里采用觀察法為模型ARIMA(p，d，q)定階，由于圖3中ACF在二階之后均落在置信區(qū)間上，所以直接定q=2或1或0;而PACF圖置信區(qū)間較小，看到在二階后落在置信區(qū)間上，但最后幾階則在置信區(qū)間外，因此需要通過時間序列相關性來判斷，判斷結果如圖4所示。由圖4可見，曲線在橫軸0處相關性最強，最接近0.5，在2處小于-0.25，即p可取的值有0和2。由于d取1，且q可取3個值，p可取2個值，所以ARIMA(p，d，q)模型有6個可取值。

表1 標準化降水指數干旱分級Table 1 Drought classification based on SPI

圖1 1951—2017年鄭州站多時間尺度SPI變化趨勢Fig.1 Trend charts of multi-time scale SPI of Zhengzhou Station from 1951 to 2017

表2 原始序列SPI12單位根檢驗Table 2 Unit root test of original sequence SPI12

圖2 一階差分(SPI12)Fig.2 First-order difference of SPI12

表3 一階差分后的SPI12單位根檢驗Table 3 Unit root test after first order difference of SPI12

2.2.2 ARIMA模型參數估計及適用性檢驗 對ARIMA(p，d，q)的6個可取值分別進行計算，將相關檢驗結果匯總(表4)。本文采用AIC、BIC準則和標準誤差來選取最優(yōu)模型，由表4可得，ARIMA(2，1，0)的AIC、BIC和標準誤差值最小，因此判定為最優(yōu)模型，即一階差分后的AR(2)模型為最優(yōu)模型。對 ARIMA(2，1，0)進行參數估計(結果見表 5)，可得ARIMA(2，1，0)模型的具體形式:

圖3 ACF和PACF(SPI12)Fig.3 Autocorrelation function and partial autocorrelation function of SPI12

圖4 PACF折線圖(SPI12)Fig.4 Line chart of PACF(SPI12)

對一階差分后的序列擬合AR(2)模型進行殘差檢驗，以評價所建模型的穩(wěn)定程度，在這里選擇QQ殘差圖及正態(tài)分布圖來檢驗模型殘差是否是平均值為0且方差為常數的正態(tài)分布，QQ殘差圖中散點均落在擬合直線附近(見圖5)，且正態(tài)圖殘差曲線也滿足正態(tài)分布(見圖6)，再通過Ljung-Box檢驗，得到p值為0.88566，遠大于0.1，綜合圖5和圖6，明顯符合白噪聲序列特征，說明該模型適用于擬合與預測SPI12的變化趨勢。選取1951—1995年數據作為訓練集，1996—2017年數據作為測試集，應用 ARIMA(2，1，0)來對 SPI12進行 22 a預測(1996-2017)，對于SPI1、SPI3、SPI6 的ARIMA 建模流程同 SPI12，選定模型分別為:SPI1(0，0，2)、SPI3(2，0，0)和 SPI6(1，0，0)，預測結果見圖 7。

表4 模型參數檢驗結果Table 4 The test results of the model parameters

表 5 ARIMA(2，1，0)模型參數Table 5 Parameters of ARIMA(2，1，0)

圖5 QQ殘差圖Fig.5 Residual chart of QQ model

圖6 殘差正態(tài)分布Fig.6 Normal distribution map of residual

2.2.3 SVR模型殘差定階及參數尋優(yōu) SVR模型對ARIMA模型殘差預測之前，首先要知道過去幾個時期的殘差會對下一個時期的殘差產生的影響，即SVR模型殘差定階。從起始時間開始每次取k個按時間順序排列的殘差數據，并依次排列下來，將排列好的矩陣作為SVR模型的輸入，k+1個數據作為模型的輸出并保留誤差。利用交叉驗證的方法，隨機取80%的數據作為模型的訓練集，其余作為測試集，利用SVR模型訓練后對測試集的結果進行預測。在選定階數為k的情況下，若RMSE第k階小于第k+1階的時候，停止循環(huán)并輸出k;反之，則繼續(xù)增加階數。對于SVR模型中的懲罰參數C和徑向基函數(Radial basis function，RBF)中的參數γ采用網格尋優(yōu)(Grid Search，GS)算法進行率定，并完成SVR模型的建模。以1951—1995年作為訓練集，1996—2017年作為測試集，利用SVR模型對ARIMA模型所預測的SPI12值的殘差進行預測，結果見圖7。

2.2.4 組合模型的預測與檢驗 在ARIMA模型的基礎上加入SVR模型進行組合預測時通常有兩種方式，一種是并聯型，一種是串聯型，并聯型是通過對兩種模型預測結果分別加上權重后重組得到;串聯型則是用ARIMA模型進行預測后所得到的殘差值輸入SVR模型再進行預測，將SVR輸出的殘差修正值和ARIMA預測值組合得到最終組合模型。兩種模型在線性和非線性預測中各有優(yōu)勢，并且ARIMA模型隨著時間增加，預測結果越趨于平穩(wěn)，在用并列組合方式給ARIMA模型權重賦值時無法隨著時間長度的增加而改變。因此本文采用串聯型進行組合預測，選取1951—1995年4個時間尺度SPI值作為訓練集，1996—2017年作為測試集，將ARIMA模型與ARIMA+SVR組合模型預測結果與實際4個時間尺度SPI值進行比對，結果見圖8。

圖7 SVR模型殘差預測(SPI12)Fig.7 Residual prediction of SVR model(SPI12)

圖8 ARIMA模型與ARIMA-SVR組合模型的多時間尺度SPI值的預測(1996-2017)Fig.8 Prediction of multi time-scale SPI based on ARIMA and ARIMA-SVR models

表6 ARIMA模型與ARIMA-SVR組合模型的均方根誤差和平均絕對百分誤差Table 6 RMSE and MAPE values for ARIMA and ARIMA-SVR model

利用RMSE和MAPE兩種評價指標對4個時間尺度的兩種模型預測結果進行評價，結果見表6。從表中數據可以發(fā)現兩種模型在SPI1時間尺度的預測效果最差，ARIMA模型在SPI12尺度預測效果最好，在 SPI3和 SPI6時間尺度預測效果僅次于SPI12。組合模型在各時間尺度預測效果都比單獨ARIMA模型預測效果好，在SPI12預測效果最好，隨著時間尺度減小，預測精度降低。

3 討論與結論

本文引入了反映干旱強度和持續(xù)時間的SPI指數，以鄭州市氣象站點為例，利用ARIMA模型和ARIMA-SVR組合模型對不同時間尺度的SPI序列進行建模預測，并利用RMSE和MAPE兩種回歸模型預測指標進行評價，得到如下結論:

1)從ARIMA模型的預測結果來看，該模型對較長時間尺度的預測精度較高，對較短時間尺度的預測精度較低，隨著時間尺度增加，擬合精度提高。對于SPI12預測效果最好，對SPI1預測效果最差，主要由以下原因導致:ARIMA模型本質上是一種整體線性自回歸模型，該模型預測趨勢會隨著測試集時間增長而逐漸趨于平穩(wěn)。因為SPI1相對于其他3個時間尺度數據量最多，整體趨于嚴平穩(wěn)(嚴平穩(wěn)表示的分布不隨時間的改變而改變)，所以預測精度最低，同理，SPI3、SPI6和SPI12時間尺度逐漸增加，數據量逐漸減少，越來越趨于弱平穩(wěn)(期望與相關系數不變，未來時刻值依賴于過去時刻的值)，因此擬合精度整體逐漸提高。

2)從組合模型的預測結果來看，對多尺度SPI值利用ARIMA模型預測線性部分，用SVR模型預測非線性部分，疊加在一起的組合模型在各個時間尺度的預測精度均比單一ARIMA模型的預測精度高。從多時間尺度來看，該模型預測精度隨著時間尺度的增加而提高，對SPI12預測效果最佳，主要是由于ARIMA模型對SPI12預測精度最高，殘差最小，所以傳入SVR模型中的誤差就小，最后組合預測結果也更加準確。同理，隨著時間尺度的減小，SPI6、SPI3和SPI1的ARIMA模型預測精度也降低，則組合模型的預測精度也隨之降低。