許德合，張 棋，黃會平

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

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

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

河南省位于中國中部和黃河中、下游，是中國重要的糧食作物產(chǎn)區(qū)和農(nóng)業(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指數(shù)公式為［23-24］:

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

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

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

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

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)對模型進行定階，差分次數(shù)記為d，ARIMA(p，d，q)模型就是經(jīng)過了d階差分后的ARMA(p，q)模型。如下式所示［25］:

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

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

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

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

1.3 SVR模型

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

式中，w為權重系數(shù)為松弛變量，規(guī)定了模型的誤差要求;C為懲罰參數(shù)，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是樣本數(shù)。RMSE和MAPE越接近0，表示預測值與觀測值越接近。

2 實例應用

2.1 數(shù)據(jù)資料及SPI計算

本文選用1951—2017年河南省國家級氣象站鄭州站逐日降水量數(shù)據(jù)來進行計算，原始數(shù)據(jù)來源于國家氣象信息中心提供的中國地面氣候資料日值數(shù)據(jù)集。運用Matlab數(shù)學建模軟件編寫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是經(jīng)過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)化處理，經(jīng)一階差分后時間序列趨于平穩(wěn)(圖2)，再對一階差分后的時間序列進行單位根檢驗(結果見表3)。表3中P值為0.00128，小于0.05，因此差分之后為平穩(wěn)時間序列。采用純隨機性檢驗(Ljung-Box Test)進行白噪聲檢驗，得P值為0.006025，遠小于0.05，因此經(jīng)一階差分后該序列為平穩(wěn)非白噪聲序列。然后用自相關函數(shù)(Autocorrelation Function，ACF)及偏自相關函數(shù)(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 標準化降水指數(shù)干旱分級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模型參數(shù)估計及適用性檢驗 對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)進行參數(shù)估計(結果見表 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且方差為常數(shù)的正態(tài)分布，QQ殘差圖中散點均落在擬合直線附近(見圖5)，且正態(tài)圖殘差曲線也滿足正態(tài)分布(見圖6)，再通過Ljung-Box檢驗，得到p值為0.88566，遠大于0.1，綜合圖5和圖6，明顯符合白噪聲序列特征，說明該模型適用于擬合與預測SPI12的變化趨勢。選取1951—1995年數(shù)據(jù)作為訓練集，1996—2017年數(shù)據(jù)作為測試集，應用 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 模型參數(shù)檢驗結果Table 4 The test results of the model parameters

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

2.2.4 組合模型的預測與檢驗 在ARIMA模型的基礎上加入SVR模型進行組合預測時通常有兩種方式，一種是并聯(lián)型，一種是串聯(lián)型，并聯(lián)型是通過對兩種模型預測結果分別加上權重后重組得到;串聯(lián)型則是用ARIMA模型進行預測后所得到的殘差值輸入SVR模型再進行預測，將SVR輸出的殘差修正值和ARIMA預測值組合得到最終組合模型。兩種模型在線性和非線性預測中各有優(yōu)勢，并且ARIMA模型隨著時間增加，預測結果越趨于平穩(wěn)，在用并列組合方式給ARIMA模型權重賦值時無法隨著時間長度的增加而改變。因此本文采用串聯(lián)型進行組合預測，選取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。從表中數(shù)據(jù)可以發(fā)現(xiàn)兩種模型在SPI1時間尺度的預測效果最差，ARIMA模型在SPI12尺度預測效果最好，在 SPI3和 SPI6時間尺度預測效果僅次于SPI12。組合模型在各時間尺度預測效果都比單獨ARIMA模型預測效果好，在SPI12預測效果最好，隨著時間尺度減小，預測精度降低。

3 討論與結論

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

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

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