ARIMA-SVR組合模型在基于標(biāo)準(zhǔn)化降水指數(shù)干旱預(yù)測(cè)中的應(yīng)用

許德合，張 棋，黃會(huì)平

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

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

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

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

1 研究方法

1.1 SPI的計(jì)算

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

SPI指數(shù)公式為［23-24］:

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

當(dāng)Y(x)＞0.5時(shí)，G=1，當(dāng)Y(x)≤0.5時(shí)，G=-1。

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

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

1.2 ARIMA模型

ARIMA分為自回歸模型(Autoregressive model，AR)、滑動(dòng)平均模型(Moving average model，MA)以及自回歸移動(dòng)平均模型(Autoregressive Moving Average Model，ARMA)，是傳統(tǒng)的時(shí)間序列預(yù)測(cè)模型。其建模流程是首先判斷模型平穩(wěn)程度，其次利用差分法對(duì)非平穩(wěn)時(shí)間序列進(jìn)行平穩(wěn)化處理，然后選取AR(p)，MA(q)對(duì)模型進(jìn)行定階，差分次數(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)為移動(dòng)平均系數(shù);ut為白噪聲序列(服從0均值、正態(tài)分布且相互獨(dú)立的白噪聲序列)。

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

1.3 SVR模型

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

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

1.4 ARIMA-SVR組合模型

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

1.5 評(píng)價(jià)驗(yàn)證指標(biāo)

在回歸模型當(dāng)中，平均絕對(duì)誤差(Mean Absolute Error，MAE)、均方誤差(Mean Square Error，MSE)、RMSE和MAPE是常見的回歸預(yù)測(cè)評(píng)估指標(biāo)，其中MAE和MSE是最基礎(chǔ)的評(píng)估方法，RMSE和MAPE是回歸任務(wù)最常用的性能度量，是前兩種指標(biāo)的擴(kuò)展，不但在多場(chǎng)景下可以使用，而且相對(duì)前兩種更加準(zhǔn)確，因此本文采用RMSE和MAEP作為模型評(píng)價(jià)的指標(biāo)。

RMSE均方根誤差

MAPE平均絕對(duì)百分誤差

式中，xi是觀測(cè)值，yi是預(yù)測(cè)值，N是樣本數(shù)。RMSE和MAPE越接近0，表示預(yù)測(cè)值與觀測(cè)值越接近。

2 實(shí)例應(yīng)用

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

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

2.2 時(shí)間序列生成及建模流程

本文使用Python 3.6平臺(tái)對(duì)ARIMA進(jìn)行建模，并利用Python中matplotlib可視化庫對(duì)多尺度SPI計(jì)算結(jié)果進(jìn)行可視化展示，如圖1所示。

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

表1 標(biāo)準(zhǔn)化降水指數(shù)干旱分級(jí)Table 1 Drought classification based on SPI

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

表2 原始序列SPI12單位根檢驗(yàn)Table 2 Unit root test of original sequence SPI12

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

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

2.2.2 ARIMA模型參數(shù)估計(jì)及適用性檢驗(yàn) 對(duì)ARIMA(p，d，q)的6個(gè)可取值分別進(jìn)行計(jì)算，將相關(guān)檢驗(yàn)結(jié)果匯總(表4)。本文采用AIC、BIC準(zhǔn)則和標(biāo)準(zhǔn)誤差來選取最優(yōu)模型，由表4可得，ARIMA(2，1，0)的AIC、BIC和標(biāo)準(zhǔn)誤差值最小，因此判定為最優(yōu)模型，即一階差分后的AR(2)模型為最優(yōu)模型。對(duì) ARIMA(2，1，0)進(jìn)行參數(shù)估計(jì)(結(jié)果見表 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)

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

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

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

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

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

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

利用RMSE和MAPE兩種評(píng)價(jià)指標(biāo)對(duì)4個(gè)時(shí)間尺度的兩種模型預(yù)測(cè)結(jié)果進(jìn)行評(píng)價(jià)，結(jié)果見表6。從表中數(shù)據(jù)可以發(fā)現(xiàn)兩種模型在SPI1時(shí)間尺度的預(yù)測(cè)效果最差，ARIMA模型在SPI12尺度預(yù)測(cè)效果最好，在 SPI3和 SPI6時(shí)間尺度預(yù)測(cè)效果僅次于SPI12。組合模型在各時(shí)間尺度預(yù)測(cè)效果都比單獨(dú)ARIMA模型預(yù)測(cè)效果好，在SPI12預(yù)測(cè)效果最好，隨著時(shí)間尺度減小，預(yù)測(cè)精度降低。

3 討論與結(jié)論

本文引入了反映干旱強(qiáng)度和持續(xù)時(shí)間的SPI指數(shù)，以鄭州市氣象站點(diǎn)為例，利用ARIMA模型和ARIMA-SVR組合模型對(duì)不同時(shí)間尺度的SPI序列進(jìn)行建模預(yù)測(cè)，并利用RMSE和MAPE兩種回歸模型預(yù)測(cè)指標(biāo)進(jìn)行評(píng)價(jià)，得到如下結(jié)論:

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

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