李 俐 許連香 王鵬新 齊 璇 王 蕾

(1.中國(guó)農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院， 北京 100083； 2.農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)災(zāi)害遙感重點(diǎn)實(shí)驗(yàn)室， 北京 100083；3.中國(guó)農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083)

0 引言

干旱是對(duì)農(nóng)業(yè)生產(chǎn)影響較大的自然災(zāi)害，其波及范圍廣、持續(xù)時(shí)間長(zhǎng)、危害性大[1]。及時(shí)準(zhǔn)確地進(jìn)行干旱預(yù)測(cè)對(duì)提前采取防災(zāi)措施以減小因旱災(zāi)造成的經(jīng)濟(jì)損失、保證糧食安全具有重要意義[2]。

干旱指數(shù)是衡量干旱程度的重要指標(biāo)，研究人員提出了許多干旱指數(shù)來(lái)監(jiān)測(cè)和預(yù)測(cè)干旱程度[3-4]。傳統(tǒng)干旱指數(shù)如帕爾默干旱嚴(yán)重程度指數(shù)(Palmer drought severity index， PDSI)、標(biāo)準(zhǔn)化降水指數(shù)(Standardized precipitation index，SPI)等[5]，通常利用氣象站或水文站等觀測(cè)數(shù)據(jù)，雖然數(shù)據(jù)較為準(zhǔn)確，但其精度受限于區(qū)域內(nèi)站點(diǎn)的位置和分布密度，大范圍的旱情監(jiān)測(cè)評(píng)估應(yīng)用中，其代表性有待提高[6]。遙感技術(shù)的快速發(fā)展使得實(shí)時(shí)、動(dòng)態(tài)的區(qū)域性旱情監(jiān)測(cè)成為可能。常用的遙感干旱監(jiān)測(cè)指數(shù)包括植被指數(shù)(Vegetation index，VI)、地表溫度(Land surface temperature，LST)，以及在此基礎(chǔ)上開發(fā)的條件植被指數(shù)(Vegetation condition index，VCI)和條件溫度指數(shù)(Temperature condition index，TCI)等[7-8]。有研究表明，植被指數(shù)對(duì)降水或土壤水分虧缺具有延遲響應(yīng)[9]，而地表溫度雖然對(duì)初始水分脅迫具有更大的敏感性[10]，但時(shí)空變化受大氣、環(huán)境狀況影響較大，僅利用植被指數(shù)或者地表溫度進(jìn)行干旱監(jiān)測(cè)并不理想。故許多學(xué)者嘗試將植被指數(shù)和地表溫度結(jié)合，利用兩者互補(bǔ)特性提供的作物水分虧缺信息來(lái)監(jiān)測(cè)旱情[11]。王鵬新等[12]在歸一化植被指數(shù)和地表溫度的散點(diǎn)圖呈三角形分布的條件下，提出了條件植被溫度指數(shù)(Vegetation temperature condition index，VTCI)的干旱監(jiān)測(cè)指數(shù)，彌補(bǔ)了單一遙感指數(shù)的不足，并且成功實(shí)現(xiàn)了陜西省關(guān)中地區(qū)冬小麥生長(zhǎng)季的干旱監(jiān)測(cè)和預(yù)測(cè)[13]。

常用的預(yù)測(cè)方法有灰色預(yù)測(cè)模型、神經(jīng)網(wǎng)絡(luò)模型[14]、加權(quán)馬爾可夫模型和求和自回歸移動(dòng)平均(Autoregressive integrated moving average，ARIMA)模型等。其中，ARIMA模型作為時(shí)間序列分析的主流方法，在很多領(lǐng)域中得到廣泛應(yīng)用[15]。艾欣等[16]利用ARIMA模型對(duì)電價(jià)時(shí)間序列進(jìn)行分析，表明ARIMA模型可提高實(shí)時(shí)市場(chǎng)的電價(jià)預(yù)測(cè)準(zhǔn)確性；竇慧麗等[17]應(yīng)用ARIMA模型，以較高精度實(shí)現(xiàn)了小波消噪后交通流的實(shí)時(shí)動(dòng)態(tài)預(yù)測(cè)。在農(nóng)業(yè)干旱預(yù)測(cè)方面，田苗等[18]基于VTCI干旱監(jiān)測(cè)結(jié)果，利用季節(jié)性求和自回歸移動(dòng)平均(Seasonal autoregressive integrated moving average，SARIMA)模型，實(shí)現(xiàn)了冬小麥生長(zhǎng)季的短期干旱預(yù)測(cè)。有研究表明，農(nóng)業(yè)旱情的發(fā)生和發(fā)展除受降水、溫度等影響外，還與農(nóng)作物自身生理機(jī)能有關(guān)[19]。在不同區(qū)域、不同作物生育期內(nèi)，作物對(duì)水分脅迫的抵抗力不同，其反映的干旱特征也存在差異[20]。本文以河北中部平原為研究區(qū)域，基于遙感反演獲取的VTCI時(shí)間序列數(shù)據(jù)，利用ARIMA模型及季節(jié)性ARIMA模型，分別逐像素對(duì)該地區(qū)的VTCI時(shí)間序列進(jìn)行分析建模預(yù)測(cè)，分析對(duì)比兩種模型的預(yù)測(cè)精度，以獲得適合河北中部平原旱情預(yù)測(cè)的模型；同時(shí)，在灌溉類型種植區(qū)，驗(yàn)證基于條件植被溫度指數(shù)(VTCI)的夏玉米生長(zhǎng)季對(duì)農(nóng)業(yè)旱情預(yù)測(cè)的有效性。

1 材料與方法

1.1 研究區(qū)概況

選取位于河北中部平原保定市、石家莊市、廊坊市的部分地區(qū)以及滄州市和衡水市為研究區(qū)(圖1)，其覆蓋范圍介于114°32′～117°36′E，36°57′～39°50′N之間，面積約為5.3×104km2。該地區(qū)屬暖溫帶大陸性季風(fēng)氣候，年平均氣溫10℃以上，主要耕作制度為冬小麥-夏玉米輪作，是我國(guó)重要的玉米生產(chǎn)基地[21]。然而，夏玉米生育期內(nèi)氣溫高且蒸散量大，降水時(shí)空分布不均以及水資源的總體不足使得該地區(qū)夏玉米生育期內(nèi)干旱時(shí)有發(fā)生[22]。灌漿成熟期是夏玉米產(chǎn)量形成的主要階段，需要充足的水分作為溶媒把莖葉中的營(yíng)養(yǎng)物質(zhì)運(yùn)輸?shù)阶蚜Ｖ?，此時(shí)土壤水分狀況較生育前期更具有重要的生理意義[23]。本文重點(diǎn)研究預(yù)測(cè)該期間(9月)的干旱時(shí)空分布情況，以期為當(dāng)?shù)剞r(nóng)業(yè)的防災(zāi)減災(zāi)提供科學(xué)指導(dǎo)和依據(jù)。

圖1 研究區(qū)域位置及氣象站點(diǎn)分布圖Fig.1 Location map of study area and weather stations

1.2 數(shù)據(jù)獲取與處理

條件植被溫度指數(shù)(VTCI)計(jì)算式[24-25]為

(1)

其中

Lmax(Ni)=a+bNi

(2)

(3)

式中Ni——第i個(gè)像素的歸一化植被指數(shù)

L(Ni)——在研究區(qū)域內(nèi)，某一像素的NDVI值為Ni時(shí)的地表溫度

Lmax(Ni)、Lmin(Ni)——當(dāng)NDVI值為某一特定值Ni時(shí)研究區(qū)內(nèi)所有像素的地表溫度的最大值和最小值，稱作熱邊界和冷邊界

a、b、a′、b′——待定系數(shù)，由研究區(qū)域的NDVI和LST的散點(diǎn)圖近似得到

選取河北省中部平原2010—2018年夏玉米生長(zhǎng)季的Aqua-MODIS數(shù)據(jù)，主要包括時(shí)間分辨率均為1 d、空間分辨率均為1 km的地表溫度產(chǎn)品MYD11A1與地表反射率產(chǎn)品MYD09GA。在進(jìn)行鑲嵌、重采樣、投影轉(zhuǎn)換以及裁剪等預(yù)處理的基礎(chǔ)上，將這兩類MODIS數(shù)據(jù)產(chǎn)品批量處理反演得到研究區(qū)域日LST產(chǎn)品和日NDVI產(chǎn)品，進(jìn)一步利用最大值合成技術(shù)，生成以旬為單位的NDVI和LST最大值合成產(chǎn)品。將多年某一旬的NDVI和LST最大值合成產(chǎn)品再一次應(yīng)用最大值合成技術(shù)，生成多年旬尺度NDVI和LST最大值合成產(chǎn)品；同時(shí)，基于最小值合成技術(shù)將多年某一旬的LST最大值合成產(chǎn)品逐像素取最小值，得到多年旬尺度LST最大-最小值合成產(chǎn)品。

利用上述NDVI和LST合成產(chǎn)品，根據(jù)VTCI計(jì)算方法生成2010—2018年每年夏玉米生長(zhǎng)季(7—9月)以旬為單位的VTCI時(shí)間序列數(shù)據(jù)，共81旬。另外，VTCI的取值范圍為[0，1]，VTCI值越小，表明作物受水分脅迫程度越重；VTCI值越大，表明作物受水分脅迫程度較輕或不受水分脅迫?？紤]到云雨因素可能造成的某些時(shí)段或者某些像素?cái)?shù)據(jù)的缺失，采用反距離加權(quán)插值法[26]對(duì)缺失數(shù)據(jù)進(jìn)行插補(bǔ)，以保證VTCI數(shù)據(jù)的完整性。

1.3 研究方法

1.3.1ARIMA預(yù)測(cè)模型

ARIMA(p,d,q)模型是由BOX和JENKINS提出的時(shí)間序列預(yù)測(cè)方法，通過(guò)對(duì)非平穩(wěn)時(shí)間序列進(jìn)行差分處理獲得平穩(wěn)序列，進(jìn)而利用自回歸移動(dòng)平均模型ARMA(p,q)實(shí)現(xiàn)對(duì)差分后平穩(wěn)時(shí)間序列未來(lái)變化的預(yù)測(cè)。p為自回歸階數(shù)，q為移動(dòng)平均階數(shù)，d為差分階數(shù)。利用ARIMA模型進(jìn)行預(yù)測(cè)的基本步驟為[27]：

(1)平穩(wěn)性檢驗(yàn)及平穩(wěn)化處理：首先檢驗(yàn)VTCI時(shí)間序列{Xt；t=1，2，…}的平穩(wěn)性，若為含有趨勢(shì)性的非平穩(wěn)時(shí)間序列，對(duì)其進(jìn)行d階差分處理將其轉(zhuǎn)換為平穩(wěn)序列{xt}，即

(4)

由于xt取值不僅與VTCI時(shí)間序列自身值有關(guān)，而且還與進(jìn)入系統(tǒng)的隨機(jī)擾動(dòng)值有關(guān)，故對(duì)平穩(wěn)時(shí)間序列{xt}擬合ARMA(p,q)模型為

(5)

φ(B)xt=θ(B)εt

(6)

(7)

式中θi、φi——模型參數(shù)εt——白噪聲序列

(2)模型定階：對(duì)于平穩(wěn)化處理后的時(shí)間序列，需要對(duì)自回歸階數(shù)p和移動(dòng)平均階數(shù)q進(jìn)行確定，即實(shí)現(xiàn)模型定階。首先，通過(guò)自相關(guān)函數(shù)(Autocorrelation function，ACF)和偏自相關(guān)函數(shù)(Partial autocorrelation function，PACF)的拖尾或截尾特征建立對(duì)應(yīng)的模型，其結(jié)構(gòu)判定的基本準(zhǔn)則如表1所示，確定p和q的取值范圍。然后，根據(jù)AIC準(zhǔn)則(Akaike information criterion，AIC)對(duì)p和q進(jìn)行優(yōu)選以保障模型的預(yù)測(cè)精度，AIC函數(shù)達(dá)到最小的模型即為最優(yōu)模型。AIC函數(shù)定義為

AIC=-2lnL()+2ω

(8)

式中L()——極大似然函數(shù)值

ω——模型參數(shù)個(gè)數(shù)

表1 ARMA模型定階基本原則Tab.1 ARMA model fixed basic principle

(3)參數(shù)估計(jì)：模型定階后，使用極大似然估計(jì)法對(duì)選定模型中的參數(shù)θi、φi進(jìn)行估計(jì)。

(4)模型檢驗(yàn)：根據(jù)殘差序列是否為白噪聲序列檢驗(yàn)?zāi)Ｐ褪欠癯浞痔崛⌒蛄兄档男畔?。采用正態(tài)分布檢驗(yàn)法，若殘差的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)值均落在95%的置信區(qū)間內(nèi)，則殘差序列為白噪聲序列，表明擬合模型有效提取了序列信息。否則，需要重新擬合模型。

圖2 逐像素建模預(yù)測(cè)流程圖Fig.2 Pixel-by-pixel forecasting flow chart

(5)模型預(yù)測(cè)：經(jīng)過(guò)步驟(1)～(4)，即可確定最優(yōu)預(yù)測(cè)模型，利用VTCI時(shí)間序列運(yùn)行此最優(yōu)模型實(shí)現(xiàn)VTCI預(yù)測(cè)。

為了有效減少金礦開采導(dǎo)致的水資源污染，應(yīng)采用水層隔離方式來(lái)減少水污染現(xiàn)象的發(fā)生，礦井排水也可以充分利用，用來(lái)灌溉農(nóng)田。通過(guò)應(yīng)用阻水技術(shù)和截流技術(shù)封閉進(jìn)水通道，避免礦井水資源泄露，保障礦金周邊能夠正常供水。針對(duì)已污染水源需要及時(shí)進(jìn)行凈化處理和水污染治理。采用物理或化學(xué)方式進(jìn)行污水治理，分段治理污水，逐步改善水環(huán)境，凈化污水。

1.3.2SARIMA預(yù)測(cè)模型

若非平穩(wěn)性時(shí)間序列不僅含有趨勢(shì)性變化，還含有季節(jié)性變化，對(duì)其進(jìn)行平穩(wěn)化處理過(guò)程中，需要在進(jìn)行d階差分的基礎(chǔ)上，再進(jìn)行D階季節(jié)性差分以消除季節(jié)性變化影響得到平穩(wěn)序列[28]，模型表示為ARIMA(p,d,q)(P,D,Q)S，具體公式為

(9)

其中

(10)

式中P——季節(jié)性自回歸階數(shù)

Q——季節(jié)性移動(dòng)平均階數(shù)

Φi、Θi——模型參數(shù)

將2010年7月上旬—2018年8月下旬的VTCI數(shù)據(jù)作為分析建模數(shù)據(jù)，2016—2018年每年9月上旬—9月下旬的VTCI數(shù)據(jù)作為檢驗(yàn)數(shù)據(jù)。逐像素提取多旬VTCI建模數(shù)據(jù)組成一維時(shí)間序列，分別作為兩個(gè)模型的輸入數(shù)據(jù)，預(yù)測(cè)流程如圖2所示，從2016—2018年每年8月下旬分別向前1步、2步和3步得到2016—2018年每年9月上旬1步預(yù)測(cè)結(jié)果、9月中旬2步預(yù)測(cè)結(jié)果和9月下旬3步預(yù)測(cè)結(jié)果。

1.4 預(yù)測(cè)精度評(píng)價(jià)

將VTCI遙感監(jiān)測(cè)結(jié)果作為真值，應(yīng)用絕對(duì)誤差(Absolute error，AE)、平均絕對(duì)誤差(Mean absolute error，MAE)與均方根誤差(Root mean square error，RMSE)評(píng)價(jià)河北中部平原夏玉米生長(zhǎng)季VTCI預(yù)測(cè)結(jié)果的精度，計(jì)算式為

AE=i-Xi

(11)

(12)

(13)

Xi——第i個(gè)像素的VTCI監(jiān)測(cè)值

N——研究區(qū)域內(nèi)所有像素點(diǎn)數(shù)(或氣象站點(diǎn)總數(shù))

2 結(jié)果與分析

2.1 時(shí)間序列長(zhǎng)度的確定

有研究表明，時(shí)間序列長(zhǎng)度是影響模型預(yù)測(cè)準(zhǔn)確性的重要因素[29]。為探究時(shí)間序列長(zhǎng)度對(duì)基于ARIMA模型的VTCI預(yù)測(cè)精度的影響，選取均勻分布在河北中部平原地區(qū)，包括饒陽(yáng)、任丘、河間、獻(xiàn)縣、涿州、雄縣、高碑店、固安、永清、霸州、晉州、辛集、藁城、深州、故城等49個(gè)氣象站點(diǎn)(圖1)，利用49個(gè)氣象站點(diǎn)所在像素的VTCI時(shí)間序列，分別使用9n(n=2，3，…，8)旬不同長(zhǎng)度的VTCI數(shù)據(jù)進(jìn)行建模預(yù)測(cè)，并分析基于ARIMA模型的VTCI預(yù)測(cè)精度隨時(shí)間序列長(zhǎng)度增加的變化特點(diǎn)。

由表2可得，同一時(shí)間序列長(zhǎng)度時(shí)，平均絕對(duì)誤差隨預(yù)測(cè)步長(zhǎng)的增加而增大，表明基于ARIMA模型的VTCI預(yù)測(cè)精度隨預(yù)測(cè)步長(zhǎng)增加而降低。不同時(shí)間序列長(zhǎng)度時(shí)，模型預(yù)測(cè)精度隨時(shí)間序列長(zhǎng)度增加上下波動(dòng)，當(dāng)序列長(zhǎng)度大于36旬時(shí)，平均絕對(duì)誤差波動(dòng)幅度逐漸減小，模型預(yù)測(cè)精度趨于穩(wěn)定。綜上，考慮到模型預(yù)測(cè)精度的穩(wěn)定性，本文分別選取2010年7月上旬至2016年8月下旬、2010年7月上旬至2017年8月下旬、2010年7月上旬至2018年8月下旬的VTCI時(shí)間序列數(shù)據(jù)進(jìn)行建模，以得到2016—2018年每年9月的VTCI預(yù)測(cè)結(jié)果。

表2 ARIMA模型平均絕對(duì)誤差隨時(shí)間序列長(zhǎng)度變化的統(tǒng)計(jì)結(jié)果Tab.2 Statistical results of ARIMA model mean absolute error varied with time series length

2.2 ARIMA模型及SARIMA模型定階結(jié)果分析

根據(jù)ARIMA模型建模方法，首先分析49個(gè)氣象站點(diǎn)所在像素的VTCI時(shí)間序列適合的模型結(jié)構(gòu)，再由點(diǎn)及面，逐像素對(duì)研究區(qū)所有像素的VTCI時(shí)間序列進(jìn)行模型定階。以饒陽(yáng)為例，其平穩(wěn)化處理后VTCI時(shí)間序列的自相關(guān)系數(shù)及偏自相關(guān)系數(shù)(圖3)未隨延遲時(shí)期增加迅速衰減至零值附近作小值波動(dòng)，均表現(xiàn)拖尾特征，表明序列適用ARMA(p,q)模型。自相關(guān)階次p和移動(dòng)平均階次q可由低階向高階逐步試探，p、q的取值范圍分別取1～3和0～2。依據(jù)AIC準(zhǔn)則進(jìn)一步進(jìn)行模型優(yōu)選，AIC值最小的模型即為該序列的最優(yōu)模型。

圖3 饒陽(yáng)氣象站差分序列的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)Fig.3 Autocorrelation and partial autocorrelation function graphs of differenced time series of VTCI in Raoyang weather station

圖4 模型定階結(jié)果Fig.4 Model identification results

逐像素對(duì)研究區(qū)所有像素進(jìn)行模型優(yōu)選，得到ARIMA模型和SARIMA模型面上定階結(jié)果(圖4)?？梢钥闯?，ARIMA模型的定階結(jié)果分布具有區(qū)域性，未出現(xiàn)“椒鹽式”分布現(xiàn)象，表明相鄰像素點(diǎn)干旱變化情況具有良好的相關(guān)性。廊坊市、滄州市、衡水市及石家莊東部等區(qū)域適合ARIMA(1,1,1)模型，模型形式較為一致。然而，保定市的模型形式存在ARIMA(1,1,1)、ARIMA(1,1,2)及ARIMA(2,1,1)等多種情況，表明受客觀環(huán)境及人為因素的影響，同一地區(qū)不同像素VTCI時(shí)間序列反映的旱情特性也存在差異性，適用的模型形式可能不同。綜上，逐像素確定ARIMA模型形式的方法較為合理。SARIMA模型的定階結(jié)果分布雖呈現(xiàn)了類似的區(qū)域性特征，但適用的模型形式更為多樣，大部分地區(qū)適用的模型為ARIMA(1,1,1)(0,1,0)9以及ARIMA(3,1,2)(0,1,0)9。整體來(lái)看，ARIMA模型定階結(jié)果較SARIMA模型區(qū)域分布特征更為明確，相鄰像素點(diǎn)間干旱變化狀況相關(guān)性更強(qiáng)。

2.3 干旱預(yù)測(cè)結(jié)果分析

2.3.1兩種模型VTCI預(yù)測(cè)結(jié)果比較

圖5 2017年9月干旱監(jiān)測(cè)結(jié)果與預(yù)測(cè)結(jié)果Fig.5 Drought monitoring results and forecasting results in September 2017

根據(jù)2.2節(jié)中兩種模型定階結(jié)果，逐像素進(jìn)行參數(shù)估計(jì)和預(yù)測(cè)，分別得到2017年9月上旬1步預(yù)測(cè)結(jié)果、9月中旬2步預(yù)測(cè)結(jié)果、9月下旬3步預(yù)測(cè)結(jié)果(圖5)。對(duì)于ARIMA模型預(yù)測(cè)結(jié)果，從時(shí)間上看，各旬VTCI均存在較大差異，9月上旬和下旬預(yù)測(cè)結(jié)果VTCI值偏高，特別是中部地區(qū)，處于不旱或較為濕潤(rùn)的狀態(tài)；而9月中旬預(yù)測(cè)結(jié)果VTCI值較上旬和下旬預(yù)測(cè)結(jié)果整體偏低，大部分地區(qū)均有旱情發(fā)生，雖然與監(jiān)測(cè)結(jié)果相比，預(yù)測(cè)旱情偏輕，但準(zhǔn)確反映了河北中部平原地區(qū)9月上旬到9月中旬旱情加重，9月中旬到9月下旬旱情有所緩解的變化趨勢(shì)。從空間分布來(lái)看，河北中部平原西北部VTCI值較東南部偏低，特別是保定市和石家莊市東部地區(qū)較其他地區(qū)干旱嚴(yán)重，與實(shí)際監(jiān)測(cè)結(jié)果一致。總體來(lái)說(shuō)，三旬的預(yù)測(cè)結(jié)果基本反映了監(jiān)測(cè)結(jié)果的特征。

比較來(lái)說(shuō)，SARIMA模型預(yù)測(cè)結(jié)果與ARIMA模型預(yù)測(cè)結(jié)果表征的旱情發(fā)展趨勢(shì)較為相似，9月中旬大部分地區(qū)均有旱情發(fā)生，與實(shí)際情況吻合。然而，石家莊及衡水市部分地區(qū)的3步預(yù)測(cè)結(jié)果顯示干旱程度加深，受旱面積也呈擴(kuò)大趨勢(shì)，與監(jiān)測(cè)結(jié)果相差較大，表明SARIMA模型3步預(yù)測(cè)結(jié)果的不確定較大，不適合更長(zhǎng)期的預(yù)測(cè)。整體來(lái)看，SARIAM模型向前1～2步預(yù)測(cè)較為準(zhǔn)確，但向前3步預(yù)測(cè)結(jié)果精度稍差。

綜上所述，ARIMA模型較SARIMA模型對(duì)河北中部平原夏玉米生長(zhǎng)季干旱的預(yù)測(cè)能力更為突出，向前1～3旬預(yù)測(cè)結(jié)果均較為準(zhǔn)確反映旱情的發(fā)展變化趨勢(shì)。

2.3.2兩種模型干旱預(yù)測(cè)結(jié)果的精度評(píng)價(jià)

為定量評(píng)價(jià)兩個(gè)模型的預(yù)測(cè)精度，以2017年9月研究區(qū)VTCI干旱監(jiān)測(cè)結(jié)果為真值，以基于兩模型分別向前預(yù)測(cè)1步、2步和3步得到的2017年9月上旬、中旬和下旬VTCI預(yù)測(cè)結(jié)果為預(yù)測(cè)值，計(jì)算得到兩模型預(yù)測(cè)結(jié)果與監(jiān)測(cè)結(jié)果的絕對(duì)誤差和絕對(duì)誤差頻數(shù)分布圖(圖6)。頻數(shù)大于500時(shí)，ARIMA模型1步預(yù)測(cè)結(jié)果絕對(duì)誤差主要分布在[-0.05，0.14]，而SARIMA模型1步預(yù)測(cè)結(jié)果主要分布在[-0.18，0.27]，較ARIMA模型誤差分布更為分散。兩者峰值雖較為接近，分別為0.04和0.07，但前者峰值頻數(shù)為6 022，而后者僅為4 122。2步和3步預(yù)測(cè)絕對(duì)誤差分布規(guī)律與1步預(yù)測(cè)相似，ARIMA模型預(yù)測(cè)誤差分布較SARIMA模型更為集中，最大頻數(shù)對(duì)應(yīng)的絕對(duì)誤差更接近零值。另外，逐像素計(jì)算并統(tǒng)計(jì)兩模型預(yù)測(cè)結(jié)果與監(jiān)測(cè)結(jié)果的平均絕對(duì)誤差以及均方根誤差(表3)，結(jié)果表明，基于ARIMA模型1步、2步、3步VTCI預(yù)測(cè)結(jié)果的平均絕對(duì)誤差和均方根誤差均低于基于SARIMA模型的誤差，平均絕對(duì)誤差分別降低0.05、0.05、0.08，均方根誤差分別降低0.06、0.07、0.09。綜上，ARIMA模型預(yù)測(cè)結(jié)果整體上精度更高，預(yù)測(cè)結(jié)果反映的旱情與實(shí)際情況更為吻合，可用于研究區(qū)夏玉米生長(zhǎng)季干旱預(yù)測(cè)。

圖6 2017年9月預(yù)測(cè)結(jié)果絕對(duì)誤差頻數(shù)分布Fig.6 Frequency distributions of absolute errors of forecasting results in September 2017

表3 ARIMA模型和SARIMA模型預(yù)測(cè)誤差的統(tǒng)計(jì)分析Tab.3 Statistical results of forecasting errors of ARIMA model and SARIMA model

2.3.3ARIMA模型干旱預(yù)測(cè)結(jié)果分析

在比較ARIMA模型和SARIMA模型預(yù)測(cè)精度的基礎(chǔ)上，利用精度較高的ARIMA模型逐像素建模預(yù)測(cè)，得到研究區(qū)域2016—2018年每年9月的VTCI干旱預(yù)測(cè)結(jié)果，計(jì)算統(tǒng)計(jì)所有像素絕對(duì)誤差絕對(duì)值及誤差在不同區(qū)間的分布情況(表4)，分析不同年份夏玉米生長(zhǎng)季VTCI的預(yù)測(cè)精度。結(jié)果表明2016—2018年各旬VTCI預(yù)測(cè)結(jié)果中，2017年9月上旬1步預(yù)測(cè)結(jié)果精度最高，僅有約0.61%的像素絕對(duì)誤差絕對(duì)值超過(guò)0.20；2016年9月下旬3步預(yù)測(cè)結(jié)果精度最低，有大約17.26%的像素絕對(duì)誤差絕對(duì)值超過(guò)0.20。其中，2016—2018年向前1步的VTCI預(yù)測(cè)結(jié)果中僅有5.84%的像素絕對(duì)誤差絕對(duì)值大于0.20，并且有62.16%的像素絕對(duì)誤差絕對(duì)值小于0.10，表明向前1步的VTCI預(yù)測(cè)精度較高。隨著預(yù)測(cè)步長(zhǎng)的增加預(yù)測(cè)精度略微下降，2步、3步預(yù)測(cè)結(jié)果中像素絕對(duì)誤差絕對(duì)值大于0.20的百分比分別為6.38%、8.72%。整體來(lái)說(shuō)，不同年份夏玉米生長(zhǎng)季ARIMA模型1步、2步、3步的VTCI預(yù)測(cè)精度均較理想。

表4 2016—2018年VTCI預(yù)測(cè)的絕對(duì)誤差區(qū)間的分布Tab.4 Distribution of absolute error interval of VTCI forecasting from 2016 to 2018 %

注：|δ|表示VTCI預(yù)測(cè)的絕對(duì)誤差絕對(duì)值。

3 討論

基于遙感反演的VTCI干旱監(jiān)測(cè)結(jié)果進(jìn)行夏玉米生長(zhǎng)季干旱預(yù)測(cè)，雖然VTCI時(shí)間序列在物理意義上具有周期為9的季節(jié)性，但SARIMA模型預(yù)測(cè)精度整體較ARIMA模型低，主要因?yàn)楹颖敝胁科皆挠衩追N植區(qū)VTCI時(shí)間序列未表現(xiàn)明顯的季節(jié)特性，若對(duì)平穩(wěn)化處理后的序列再進(jìn)行季節(jié)差分，會(huì)因?yàn)椴罘诌^(guò)程中信息的損失出現(xiàn)過(guò)度差分(簡(jiǎn)稱過(guò)差分)的現(xiàn)象，從而影響模型預(yù)測(cè)精度。

另外，前人研究已表明干旱預(yù)測(cè)是屬于帶有強(qiáng)非線性特征的系統(tǒng)問(wèn)題[30]，ARIMA模型作為一種線性預(yù)測(cè)方法，會(huì)忽略VTCI時(shí)間序列數(shù)據(jù)中的非線性部分，在準(zhǔn)確預(yù)測(cè)旱情發(fā)展趨勢(shì)方面具有一定劣勢(shì)。所以在后期的研究中，可利用在非線性系統(tǒng)預(yù)測(cè)方面具有較強(qiáng)優(yōu)勢(shì)的神經(jīng)網(wǎng)絡(luò)等方法對(duì)ARIMA建模過(guò)程中的未擬合的非線性誤差進(jìn)行修正，以提高干旱預(yù)測(cè)的精度。

4 結(jié)論

(1)基于49個(gè)氣象站點(diǎn)所在像素的VTCI時(shí)間序列，使用不同時(shí)間序列長(zhǎng)度的VTCI數(shù)據(jù)進(jìn)行建模預(yù)測(cè)，結(jié)果表明，基于ARIMA模型的VTCI預(yù)測(cè)精度與時(shí)間序列長(zhǎng)度未呈現(xiàn)明顯的相關(guān)關(guān)系，但模型預(yù)測(cè)精度隨序列長(zhǎng)度的增加而逐漸趨于穩(wěn)定。

(2)將ARIMA模型和SARIMA模型分別用于河北中部平原2017年夏玉米生長(zhǎng)季VTCI預(yù)測(cè)，結(jié)果表明，ARIMA模型較SARIMA模型VTCI預(yù)測(cè)精度更高，更適合河北中部平原夏玉米生長(zhǎng)季的干旱預(yù)測(cè)，其1步、2步、3步VTCI預(yù)測(cè)的均方根誤差分別為0.07、0.14、0.10。

(3)應(yīng)用ARIMA模型逐像素對(duì)2016—2018年夏玉米生長(zhǎng)季VTCI進(jìn)行預(yù)測(cè)，結(jié)果表明，不同年份夏玉米生長(zhǎng)季VTCI的預(yù)測(cè)精度穩(wěn)定性較好，其2016—2018年1步、2步和3步VTCI預(yù)測(cè)結(jié)果絕對(duì)誤差絕對(duì)值大于0.20的像素平均占比分別為5.84%、6.38%、8.72%。