基于條件植被溫度指數的夏玉米生長季干旱預測研究

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

(1.中國農業大學土地科學與技術學院， 北京 100083； 2.農業農村部農業災害遙感重點實驗室， 北京 100083；3.中國農業大學信息與電氣工程學院， 北京 100083)

0 引言

干旱是對農業生產影響較大的自然災害，其波及范圍廣、持續時間長、危害性大[1]。及時準確地進行干旱預測對提前采取防災措施以減小因旱災造成的經濟損失、保證糧食安全具有重要意義[2]。

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

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

1 材料與方法

1.1 研究區概況

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

圖1 研究區域位置及氣象站點分布圖Fig.1 Location map of study area and weather stations

1.2 數據獲取與處理

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

(1)

其中

Lmax(Ni)=a+bNi

(2)

(3)

式中Ni——第i個像素的歸一化植被指數

L(Ni)——在研究區域內，某一像素的NDVI值為Ni時的地表溫度

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

a、b、a′、b′——待定系數，由研究區域的NDVI和LST的散點圖近似得到

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

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

1.3 研究方法

1.3.1ARIMA預測模型

ARIMA(p,d,q)模型是由BOX和JENKINS提出的時間序列預測方法，通過對非平穩時間序列進行差分處理獲得平穩序列，進而利用自回歸移動平均模型ARMA(p,q)實現對差分后平穩時間序列未來變化的預測。p為自回歸階數，q為移動平均階數，d為差分階數。利用ARIMA模型進行預測的基本步驟為[27]：

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

(4)

由于xt取值不僅與VTCI時間序列自身值有關，而且還與進入系統的隨機擾動值有關，故對平穩時間序列{xt}擬合ARMA(p,q)模型為

(5)

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

(6)

(7)

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

(2)模型定階：對于平穩化處理后的時間序列，需要對自回歸階數p和移動平均階數q進行確定，即實現模型定階。首先，通過自相關函數(Autocorrelation function，ACF)和偏自相關函數(Partial autocorrelation function，PACF)的拖尾或截尾特征建立對應的模型，其結構判定的基本準則如表1所示，確定p和q的取值范圍。然后，根據AIC準則(Akaike information criterion，AIC)對p和q進行優選以保障模型的預測精度，AIC函數達到最小的模型即為最優模型。AIC函數定義為

AIC=-2lnL()+2ω

(8)

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

ω——模型參數個數

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

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

(4)模型檢驗：根據殘差序列是否為白噪聲序列檢驗模型是否充分提取序列值的信息。采用正態分布檢驗法，若殘差的自相關函數和偏自相關函數值均落在95%的置信區間內，則殘差序列為白噪聲序列，表明擬合模型有效提取了序列信息。否則，需要重新擬合模型。

圖2 逐像素建模預測流程圖Fig.2 Pixel-by-pixel forecasting flow chart

(5)模型預測：經過步驟(1)～(4)，即可確定最優預測模型，利用VTCI時間序列運行此最優模型實現VTCI預測。

為了有效減少金礦開采導致的水資源污染，應采用水層隔離方式來減少水污染現象的發生，礦井排水也可以充分利用，用來灌溉農田。通過應用阻水技術和截流技術封閉進水通道，避免礦井水資源泄露，保障礦金周邊能夠正常供水。針對已污染水源需要及時進行凈化處理和水污染治理。采用物理或化學方式進行污水治理，分段治理污水，逐步改善水環境，凈化污水。

1.3.2SARIMA預測模型

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

(9)

其中

(10)

式中P——季節性自回歸階數

Q——季節性移動平均階數

Φi、Θi——模型參數

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

1.4 預測精度評價

將VTCI遙感監測結果作為真值，應用絕對誤差(Absolute error，AE)、平均絕對誤差(Mean absolute error，MAE)與均方根誤差(Root mean square error，RMSE)評價河北中部平原夏玉米生長季VTCI預測結果的精度，計算式為

AE=i-Xi

(11)

(12)

(13)

Xi——第i個像素的VTCI監測值

N——研究區域內所有像素點數(或氣象站點總數)

2 結果與分析

2.1 時間序列長度的確定

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

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

表2 ARIMA模型平均絕對誤差隨時間序列長度變化的統計結果Tab.2 Statistical results of ARIMA model mean absolute error varied with time series length

2.2 ARIMA模型及SARIMA模型定階結果分析

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

圖3 饒陽氣象站差分序列的自相關函數和偏自相關函數Fig.3 Autocorrelation and partial autocorrelation function graphs of differenced time series of VTCI in Raoyang weather station

圖4 模型定階結果Fig.4 Model identification results

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

2.3 干旱預測結果分析

2.3.1兩種模型VTCI預測結果比較

圖5 2017年9月干旱監測結果與預測結果Fig.5 Drought monitoring results and forecasting results in September 2017

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

比較來說，SARIMA模型預測結果與ARIMA模型預測結果表征的旱情發展趨勢較為相似，9月中旬大部分地區均有旱情發生，與實際情況吻合。然而，石家莊及衡水市部分地區的3步預測結果顯示干旱程度加深，受旱面積也呈擴大趨勢，與監測結果相差較大，表明SARIMA模型3步預測結果的不確定較大，不適合更長期的預測。整體來看，SARIAM模型向前1～2步預測較為準確，但向前3步預測結果精度稍差。

綜上所述，ARIMA模型較SARIMA模型對河北中部平原夏玉米生長季干旱的預測能力更為突出，向前1～3旬預測結果均較為準確反映旱情的發展變化趨勢。

2.3.2兩種模型干旱預測結果的精度評價

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

圖6 2017年9月預測結果絕對誤差頻數分布Fig.6 Frequency distributions of absolute errors of forecasting results in September 2017

表3 ARIMA模型和SARIMA模型預測誤差的統計分析Tab.3 Statistical results of forecasting errors of ARIMA model and SARIMA model

2.3.3ARIMA模型干旱預測結果分析

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

表4 2016—2018年VTCI預測的絕對誤差區間的分布Tab.4 Distribution of absolute error interval of VTCI forecasting from 2016 to 2018 %

注：|δ|表示VTCI預測的絕對誤差絕對值。

3 討論

基于遙感反演的VTCI干旱監測結果進行夏玉米生長季干旱預測，雖然VTCI時間序列在物理意義上具有周期為9的季節性，但SARIMA模型預測精度整體較ARIMA模型低，主要因為河北中部平原夏玉米種植區VTCI時間序列未表現明顯的季節特性，若對平穩化處理后的序列再進行季節差分，會因為差分過程中信息的損失出現過度差分(簡稱過差分)的現象，從而影響模型預測精度。

另外，前人研究已表明干旱預測是屬于帶有強非線性特征的系統問題[30]，ARIMA模型作為一種線性預測方法，會忽略VTCI時間序列數據中的非線性部分，在準確預測旱情發展趨勢方面具有一定劣勢。所以在后期的研究中，可利用在非線性系統預測方面具有較強優勢的神經網絡等方法對ARIMA建模過程中的未擬合的非線性誤差進行修正，以提高干旱預測的精度。

4 結論

(1)基于49個氣象站點所在像素的VTCI時間序列，使用不同時間序列長度的VTCI數據進行建模預測，結果表明，基于ARIMA模型的VTCI預測精度與時間序列長度未呈現明顯的相關關系，但模型預測精度隨序列長度的增加而逐漸趨于穩定。

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

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