基于多源遙感數據和隨機森林的綜合旱情指標構建

董 婷 任 東 邵 攀 孟令奎

(1.三峽大學計算機與信息學院， 宜昌 443002； 2.武漢大學遙感信息工程學院， 武漢 430079)

0 引言

干旱是全球最常見的自然災害之一，持續時間長，波及范圍廣[1-2]。受全球氣候變暖影響，近年來干旱災害呈上升趨勢[3]。旱災不僅影響農業生產，造成直接經濟損失，而且給社會發展和人民生活帶來巨大影響[4-5]。研究干旱時空演變特征、提高預警預報水平和能力具有重要的現實意義。遙感技術具有快速、經濟、宏觀、數據連續等優勢，是當前實現大面積農業旱情監測的重要手段[6-7]。

植被長勢與土壤水分密切相關，通過植被指數比較不同時期植被長勢情況，是遙感監測旱情狀態的重要途徑，常用植被指數包括歸一化植被指數(Normalized difference vegetation index, NDVI)[8]、增強型植被指數(Enhanced vegetation index, EVI)[9]、植被狀態指數(Vegetation condition index, VCI)[10]等。地表溫度(Land surface temperature, LST)也是旱情監測需考慮的一個重要因素，KOGAN[11]基于地表溫度建立溫度狀態指數(Temperature condition index, TCI)，并將VCI和TCI進行線性加權，提出植被健康指數(Vegetation health index, VHI)。農業干旱是因長時間降水異常短缺造成土壤缺水，導致作物產量減少甚至無收的一種氣象災害[12]，氣候降水和土壤濕度是農業干旱的重要關聯因子。DU等[13]在構建綜合旱情指數時考慮降水因子，利用TRMM降水數據建立降水狀態指數(Precipitation condition index, PCI)，效果顯著。ZHANG等[14]對AMSR-E土壤濕度數據進行歸一化，得到土壤濕度狀態指數(Soil moisture condition index, SMCI)；通過對SMCI、TCI和PCI進行線性加權構建微波集成旱情指數(Microwave integrated drought index, MIDI)，實驗表明，由MIDI反演得到的旱情空間分布與實測值有較好的一致性。

隨著長時序遙感數據的積累，除了簡單線性加權方法外，利用機器學習方法構建旱情監測模型已成為實施遙感旱情監測的另一主要途徑[15-17]。然而，當前基于機器學習的旱情監測方法多采用單一分類器，由于農業旱情關聯因子的多樣性和不確定性，單分類器旱情監測方法的監測效果不夠理想。集成學習通過將多個分類器聯合，能夠減少單分類器的泛化誤差[18-19]，為提高遙感旱情監測精度提供一種可能途徑。隨機森林是BREIMAN[20]提出的一種多決策樹集成學習方法，具有較高的分類精度和擬合能力[21-22]，已在土地覆蓋分類[23-25]、生物量估計[26-27]等環境遙感領域得到廣泛應用，但在旱情監測方面的相關研究較少。

本文以實測標準化降水蒸散指數作為旱情參考指標，利用隨機森林方法構建一種綜合旱情監測模型。首先通過相關性分析對不同遙感旱情單因子指數的旱情監測效果進行對比，選取特定時間的建模因子，然后利用隨機森林回歸算法對多源旱情相關信息進行融合，建立一種多因子集成的旱情狀態指數(Integrated drought condition index, IDCI-RF)，最后利用IDCI-RF對中國北部區域大范圍旱情空間分布進行評估，以期為實施旱災狀態監測、制定防災減災決策提供參考和依據。

1 研究區與數據

1.1 研究區概況

本文針對大空間范圍農業旱情監測展開研究，選擇北方地區作為研究區(圖1)，處于34°38′～45°46′N、108°7′～130°44′E，覆蓋范圍包括吉林省、遼寧省、北京市、天津市、河北省、山西省、山東省，以及內蒙古自治區、黑龍江省、陜西省、河南省部分區域。

圖1 研究區概況Fig.1 Maps of study area

依據全球MODIS土地覆蓋分類產品MCD12C1的分類結果，研究區主要土地覆蓋類型為草原、作物區和林地，其中，草原地區主要分布在內蒙古自治區、山西省北部和河北省北部部分地區，山東省以作物覆蓋為主，作物區還分布在河北省南部、吉林省和遼寧省部分區域(圖1b)。根據柯本氣候分類法，研究區的氣候環境主要是干帶草原氣候和冬干冷溫氣候，不同氣候類型區域處的降水與氣溫存在較大差異(圖1c)。受東亞季風氣候影響，研究區降水不均，易受旱澇等自然災害影響。

1.2 氣象數據

從氣象數據共享服務網(http:∥cdc.cma.gov.cn/)獲取研究區內所有氣象站點1961—2012年間月降水量以及平均氣溫，通過對所有站點氣象數據整理和分析，確定132個包含所有數據記錄的氣象站(圖1a)。分別計算各省份所有站點的月降水量和平均氣溫算術平均值，以評估該省份的區域旱情狀況。

VICENTE-SERRANO等[28]提出標準化降水蒸散指數(Standardized precipitation evapotranspiration index, SPEI)。SPEI結合了帕默爾干旱指數PDSI對蒸散的響應以及標準化降水指數SPI多時間尺度等優點，其監測效果已在全球區域的應用中得到了很好的驗證[29-30]。本文利用站點長時序月降水量和平均氣溫計算132個站點處3種時間尺度(1個月、3個月和6個月)的SPEI值，采用THORNTHWAITE方法[31]估算潛在蒸發量PET。

1.3 遙感數據及預處理

遙感數據包括MODIS植被指數產品MOD13C2和地表溫度產品MOD11C3(2003—2012年)，分別提供歸一化植被指數NDVI、增強型植被指數EVI以及地表溫度LST。利用質量控制文件對數據中噪聲點進行掩膜，并對NDVI和LST數據進行歸一化處理得到植被狀態指數VCI和溫度狀態指數TCI。利用空間數據聚合方法對5種指數(NDVI、EVI、VCI、LST、TCI)進行重采樣，使空間分辨率變為0.25°。

遙感降水量數據采用2003—2012年的熱帶降雨測量衛星TRMM 3B43數據，首先將降水速率轉換為單月降水總量數據，同時計算出1個月、3個月和6個月的累積降水量，記為PRCP1、PRCP3和PRCP6，利用3種累積降水量建立3種降水狀態指數，記為PCI1、PCI3和PCI6。

遙感土壤濕度數據采用歐空局氣候變化計劃項目提供的全球土壤濕度數據(2003—2012年)，包括3類數據集：主動微波數據集、被動微波數據集以及主、被動微波集成數據集。計算3類土壤濕度數據月均值作為月土壤濕度信息(記為SM-A、SM-P以及SM-M)，在此基礎上進行研究區裁剪和歸一化處理得到土壤濕度狀態指數，記為SMCI-A、SMCI-P以及SMCI-M。

歸一化方法如表1所示，歸一化后的指數取值范圍是[0,1]，其中0和1分別代表最干旱和最濕潤狀態。

本文還獲取土地覆蓋類型、地表高程和氣候分類數據作為模型輸入參數。其中土地覆蓋類型采用2003—2012年的MODIS土地覆蓋類型產品MCD12C1 IGBP數據；地表高程數據采用西部環境與生態科學數據中心提供的1 km分辨率中國數字高程模型DEM數據集；氣候分類信息主要依據柯本氣候分類法，使用PEEL等[32]提供的Koppen氣候分類圖作為研究區氣候分類參考數據，其空間分辨率為0.1°。利用眾數算法將土地覆蓋類型和氣候分類數據空間分辨率重采樣為0.25°，利用3次卷積法將地表高程數據空間分辨率重采樣為0.25°。

表1 干旱指數計算方法Tab.1 Remote sensing drought index formulas

2 研究方法

2.1 技術流程

利用隨機森林回歸算法集成多源遙感數據進行農業旱情監測研究，技術流程如圖2所示。

首先利用Pearson相關性分析方法對比17種不同遙感旱情指數在不同時期的旱情監測效果，確定最優建模因子；然后利用隨機森林回歸方法對不同旱情因子進行學習并對實測指數SPEI進行擬合得到回歸模型；基于回歸模型對研究區2003—2012年4—10月期間所有遙感影像進行處理，得到多因子集成旱情指數IDCI-RF；最后對IDCI-RF指數進行驗證與評價。驗證與評價方法包括：對IDCI-RF指數與實測指數進行散點圖分析和相關性驗證，利用驗證樣本集評估IDCI-RF的泛化能力；與Bagging和Cubist回歸方法進行對比，檢驗IDCI-RF指數旱情監測精度；討論IDCI-RF指數在不同空間尺度的旱情監測適用性，同時對比IDCI-RF旱情監測圖與站點實測指數分布情況，探討IDCI-RF指數在大空間范圍內的旱情監測能力。

圖2 總體技術流程Fig.2 Flow chart of technical process

2.2 隨機森林回歸方法

隨機森林方法是建立在分類回歸樹基礎上的多決策樹集成學習方法，在構建過程中隨機選擇訓練樣本子集和特征變量，使得該方法具有較高的抗噪能力，能夠有效避免分類回歸樹方法的過擬合問題[20]。

利用R語言實現隨機森林回歸，基于該方法構建實測指數反演模型，將樣本按1∶1的比例隨機劃分，得到訓練樣本集和驗證樣本集，為了驗證模型的可靠性，本文重復進行5次隨機分組，依次記為樣本組1～5。其中，訓練樣本用于建模，驗證樣本用于評估模型精度以及泛化能力等。

2.3 模型評價方法

采用均方誤差(Mean square error, MSE)作為確定隨機森林模型參數取值的統計依據；采用決定系數(Coefficient of determination,R2)、均方根誤差(Root mean square error, RMSE)和平均絕對誤差(Mean absolute error, MAE)作為評價和驗證模型精度的指標。決定系數R2反映數據擬合程度，取值越大，說明擬合值與實測值的相關性越好；均方誤差MSE、均方根誤差RMSE和平均絕對誤差MAE可判斷擬合值與實測值的差異，取值越小，說明擬合值越接近實測值。

3 結果與分析

3.1 旱情指數敏感性分析及建模因子確定

在不同旱情階段植被、降水、氣溫、土壤等不同旱情關聯因子的表現不同。分別提取2003—2012年間研究區內各站點位置處的17種遙感指數，指數包括NDVI、VCI、EVI、LST、TCI、PRCP1、PRCP3、PRCP6、PCI1、PCI3、PCI6、SM-A、SM-P、SM-M、SMCI-A、SMCI-P和SMCI-M，按照不同月份將各遙感指數分別與3種時間尺度(1個月、3個月和6個月)的標準化降水蒸散指數SPEI進行Pearson相關性分析，依據相關性結果確定建模因子。

圖3(圖中*表示P>0.01)為相關性分析結果。由圖3可知，與實測數據的相關性，除了植被指數(NDVI、EVI、VCI)在植被生長季初期與末期以及土壤濕度指數SM-P在5月外，其他的相關性均通過0.01顯著性檢驗。

圖3 遙感指數與不同時間尺度SPEI的相關性分析結果(4—10月)Fig.3 Correlations between remote sensing drought indices and in situ different time scales SPEI from April to October

3種植被指數中，VCI與不同時間尺度SPEI的相關性都強于NDVI和EVI；TCI與實測指數相關性也優于LST。對于降水指數，與歸一化前降水指數相比，歸一化后的降水指數(PCI3和PCI6)與SPEI-3、SPEI-6的相關性更高。此外，由對應月份累積起來的降水指數與其相對應時間尺度SPEI的相關性強于其與其他時間尺度SPEI相關性，如由3個月累積降水量建立起來的歸一化降水指數PCI3與SPEI-3的相關性都高于其與SPEI-1和SPEI-6的相關性，PCI1與SPEI-1的相關性也高于其與另外兩種時間尺度SPEI的相關性。土壤濕度類指數也表現出類似的規律，歸一化后的指數與3種尺度SPEI的相關性均優于歸一化前的指數。

由圖3可知，從4月到7月再到10月，植被類指數與不同時間尺度SPEI的相關性經歷先增大后減小的過程，以SPEI-6為例，在4月，VCI與SPEI-6相關性未通過0.01顯著性檢驗，而在6月，相關性增強(r為0.31)，到7月植被生長旺盛時期，相關系數達到最大(r為0.4)，8月開始，VCI與SPEI-6的相關性又逐漸減弱。這說明以植被指數建立起來的旱情模型在植被生長較為旺盛的階段監測結果更具參考價值，與JI等[33]的研究結果一致。溫度類指數在不同時間也表現出不同的相關性，其在植被生長季初期與實測指數的相關系數比在植被生長季中后期高，如TCI與3種不同時間尺度SPEI的相關性均在6月達到最高。不同類型的土壤濕度指數與實測指數相關性表現也不同，以SPEI-3為例，在6—10月期間，主動土壤濕度指數與SPEI-3的相關性明顯優于被動土壤濕度指數和主被動合成土壤濕度指數，說明主動土壤濕度指數適合在植被覆蓋度較高的區域進行旱情監測，YUAN等[34]的研究工作也得出相類似的結論。

不同時間尺度的實測標準化降水蒸散指數SPEI可以分析不同的干旱類型，小時間尺度數據適合研究土壤含水率和河道流量，大時間尺度數據可分析地下水變化情況[35-36]。本文重點研究農業旱情，IDCI-RF綜合指數的構建以SPEI-3為依據?；诟黝愡b感指數與SPEI-3的相關性分析結果確定不同月份的模型輸入數據，如表2所示。此外，旱情的發生發展與下墊面等環境因素也相關，IDCI-RF除考慮能表征植被、地表溫度、降水和土壤濕度狀態的4類遙感指數外，還考慮地表高程、土地覆蓋類型和氣候分類信息3類輔助參數。

3.2 RF參數優化及回歸模型構建

表2 4—10月隨機森林模型輸入參數Tab.2 Inputting parameters of random forest model from April to October

圖4 不同模型參數與模型袋外誤差的關系Fig.4 Relationship between changing parameters and out-of-bag error

經驗值不一定能獲得最優模型，本文在RF回歸模型構建中，為獲得最優mtry，采取逐一增加變量的方法，即對于所有mtry取值分別建立回歸模型，找到對應模型誤差最小的mtry值作為其最優參數值。當分枝創建中自變量數目mtry最優參數值確定后，通過調節ntree，即決策樹的數量，基于可視化分析觀察模型誤差變化趨勢以確定ntree最優參數值。圖4分別對應6月5組不同訓練模型得到的不同模型參數與模型袋外數據集均方誤差MSE的關系，可以發現，當mtry取值依次為2、3、3、2、3時，對應各自模型的MSE最??；當ntree達到1 000時，模型MSE都趨于穩定，因此，6月由5組訓練集訓練的回歸模型中ntree確定為1 000，mtry分別設置為2、3、3、2、3，同樣方法得到另外30組模型參數。

最終得到7個月的實測指數回歸模型，其中每月基于不同訓練集依次對應5個回歸模型，模型擬合結果稱為多因子集成旱情狀態指數IDCI-RF?；谝陨匣貧w模型，得到驗證樣本IDCI-RF值，與驗證樣本SPEI-3實測值進行對比，評估回歸模型泛化能力。圖5為各月份第3組驗證樣本集的實測指數SPEI-3與其擬合數據IDCI-RF的散點圖，可以看出所有月份的相關系數都不小于0.75，r依次為0.81、0.81、0.76、0.79、0.75、0.79和0.84，所有相關關系都通過0.005顯著性水平檢驗，說明由隨機森林算法得到的擬合模型能很好地估測實測指數，可用于實際旱情監測應用中。

圖5 4—10月驗證樣本實測數據SPEI-3與擬合數據IDCI-RF散點圖分析Fig.5 Scatter plots and correlation coefficient r values between IDCI-RF and SPEI-3 of validation samples from April to October

3.3 不同回歸方法結果對比

為檢驗隨機森林算法對旱情指數的擬合精度，選擇Cubist和Bagging兩種集成算法與RF進行對比，Cubist和Bagging模型輸入數據可以是連續值和離散值，在土地利用分類、環境監測等領域得到廣泛應用[37-38]。分別利用Cubist和Bagging模型對不同月份的5組訓練集進行學習得到回歸模型，模型擬合結果依次記為IDCI-Cubist和IDCI-Bagging。表3為4—10月不同模型估測的實測指數統計結果。

由表3可以看出，由RF構建的35組回歸模型中，IDCI-RF與實測SPEI-3的決定系數R2在0.49～0.71之間，其中，4月、9月和10月的R2都在0.6以上；由Cubist和Bagging回歸得到的IDCI與實測值的R2范圍依次為0.46～0.69和0.48～0.69，低于RF模型，說明RF估測的擬合值與實測值的相關性更好。3種算法各自建立的35組回歸模型中，RMSE和MAE最大值分別為0.69和0.54(RF)、0.72和0.55(Cubist)、0.71和0.56(Bagging)，說明RF估測的擬合值與實測值的差異更小。

從不同月份的不同樣本組統計得到的模型精度來看，除了10月(樣本組1和樣本組4)以及8月(樣本組1)以外，其他所有的回歸模型中，由RF得到的IDCI-RF與實測值決定系數R2在3種算法中都能達到最高，RF模型表現出更強的預測能力。本文也針對不同月份求得5個樣本組的統計指標平均值(表3)，4—10月IDCI-RF與實測值的平均R2依次為0.64、0.63、0.57、0.59、0.54、0.62和0.68，均優于Cubist和Bagging擬合指數，IDCI-RF對應的平均誤差指標也最小。分析表明與Cubist和Bagging方法相比，由RF算法得到的IDCI-RF與SPEI-3的一致性更強，用IDCI-RF指數估測SPEI-3指數合理。

3.4 IDCI-RF在不同空間尺度上的敏感性分析

根據各月份5組不同樣本集得到的統計評價指標結果，將由第3、3、1、3、1、1、3個樣本組構建的回歸模型依次確定為4—10月的最優模型。分別提取研究區內各省份所有站點的由最優模型得到的IDCI-RF以及實測指數SPEI-3值，利用散點圖分析和統計指標驗證方法對IDCI-RF指數在不同省級區域的監測精度進行評估。

圖6為7個主要省份(河北省、山西省、內蒙古自治區、遼寧省、吉林省、山東省和陜西省)2003—2012年植被生長季IDCI-RF指數與實測指數的散點圖以及統計值驗證結果，可以看出IDCI-RF指數在不同省份都能很好地擬合實測指數，所有相關關系都通過0.005顯著性檢驗，決定系數R2都在0.7以上，其中在內蒙古自治區(n=1 557)和遼寧省(n=1 039)(圖6c、6d)，IDCI-RF與SPEI-3的相關性最強(R2=0.82)。除吉林省和陜西省外，其他省的實測值與IDCI-RF值的誤差指標均在0.5(RMSE)和0.35(MAE)以下。

表3 不同模型實測指數回歸統計結果Tab.3 Statistics results for different regression models of in situ reference data

圖6 IDCI-RF指數與實測指數SPEI-3的散點圖Fig.6 Scatter plots between IDCI-RF and SPEI-3

為驗證站點尺度IDCI-RF指數效果，選取研究區部分站點，針對站點處IDCI-RF與實測指數的變化趨勢展開分析，選取的站點分散在研究區各區域，能較好地代表研究區不同土地覆蓋類型。圖7為21個站點在2005—2012年7月時的IDCI-RF與實測SPEI-3變化趨勢。

從圖7可以看出，IDCI-RF指數能較好地擬合實測指數，但也存在少部分站點處的IDCI-RF與SPEI-3變化不一致，如53480站點的實測SPEI-3在2009年7月至2010年7月呈上升趨勢，但其在2010年7月的IDCI-RF值比2009年同期稍低(圖7d)；54311站點處SPEI-3在2007年7月到2008年7月呈上升趨勢，而IDCI-RF指數在相應時間的變化趨勢是略下降的(圖7p)；54826站點處的實測指數在2007年7月到2009年7月呈先下降后上升趨勢，但IDCI-RF指數表現的變化趨勢與實際相反(圖7u)。盡管存在少數不一致現象，但在絕大部分情況下，IDCI-RF指數的變化規律能與實測指數SPEI-3保持一致，還有部分站點處的IDCI-RF值與實測值接近(圖7a、7e、7o、7r、7t)，說明IDCI-RF指數能較精準地擬合實測指數，可用于判斷實際旱情狀態變化趨勢。

圖7 不同站點IDCI-RF與SPEI-3的變化趨勢(2005—2012年7月)Fig.7 Year-to-year changes in IDCI-RF and SPEI-3 in representative meteorological sites for July from 2005 to 2012

3.5 旱情空間分布格局

基于研究區2003—2012年多源遙感數據集和3類輔助參數，利用不同月份最優回歸模型得到研究區所有位置處IDCI-RF值，分析區域內旱情空間分布特征。以干旱年2009年為例，圖8為該年度4—10月的IDCI-RF監測圖，為評估IDCI-RF旱情監測精度，不同站點處的實測指數SPEI-3值也依據點的大小進行區分，站點位置處點越大，說明實測指數值越小，旱情越嚴重，反之亦然。

圖8 2009年4—10月IDCI-RF空間分布圖以及實測指數SPEI-3分布情況Fig.8 Seasonal changes in drought detected by IDCI-RF and SPEI-3 from April to October in 2009

由圖8可知，內蒙古自治區在2009年8—10月旱情有加重趨勢，尤其是研究區內內蒙古自治區北部區域，這與實測SPEI-3空間分布一致，該區域站點位置處點變大，說明旱情越來越嚴重(圖8e～8g)。IDCI-RF旱情圖顯示2009年6月陜西省北部出現嚴重旱情，從站點大小看出，站點處SPEI-3值較小(圖8c)，也說明旱情嚴重。2009年4月遼寧省所有站點處SPEI-3指數都大于-0.5，IDCI-RF監測圖反映的2009年4月遼寧省旱情與實測數據一致(圖8a)；從7月開始，IDCI-RF監測圖反映出遼寧省出現旱情并表現出加重趨勢，SPEI-3也反映出同樣的旱情形勢(圖8d～8f)。

從以上分析可知，由IDCI-RF監測圖反映的研究區旱情狀態與實測指數SPEI-3分布特征有較好的一致性，2009年4—10月的IDCI-RF旱情空間分布圖能較精確地指出研究區干旱受災區域，也能在一定程度上描繪出不同區域旱情程度及變化趨勢，說明IDCI-RF指數可應用于大范圍區域農業旱情狀態監測業務中。

4 討論

農業干旱是一種異常復雜的自然災害，以往的遙感旱情監測指數多側重于考慮植被、溫度、土壤等單一因素，這類指數僅是旱情某方面、某個層次特征的表征，影響監測效果。隨著歷史數據的積累，機器學習技術被應用于旱情監測中，目前研究多數采用單一分類器對旱情特征進行挖掘。

選擇隨機森林方法，對多源遙感旱情信息進行集成構建綜合旱情指數，取得較優的監測效果。原因為：綜合考慮旱情關聯因子的多樣性，選取特征集包括能表征植被狀態、地表溫度、氣候降水、土壤濕度的多種旱情關聯因子，以及3種生態環境因素(地表高程、土地覆蓋類型、氣候類型)；選擇的建模方法隨機森林依據分類回歸樹構建，是一種集成學習方法，與單分類器相比，能夠實現信息優勢互補，獲得更高的擬合精度。下一步嘗試將其他集成學習方法應用到遙感旱情監測中。

采用的遙感數據空間分辨率較低(0.25°×0.25°)，混合像元對旱情監測效果有一定影響，后期工作中將采用更高空間分辨率的遙感影像作為基礎數據源來開展旱災監測研究。

5 結論

(1) 利用RF、Cubist和Bagging 3種方法分別建立回歸模型，利用RF構建的IDCI-RF與實測SPEI-3的平均決定系數R2為0.54～0.68，優于另外兩種方法，RF模型表現出更強的預測能力。

(2) 研究區主要省份的IDCI-RF指數與實測

SPEI-3的決定系數R2都在0.7以上，大部分站點的IDCI-RF變化規律與實測SPEI-3保持一致，IDCI-RF值與實測值接近；由IDCI-RF監測圖反映的旱情狀態與SPEI-3空間分布特征吻合度較高，IDCI-RF監測圖能較準確地指出干旱受災區域，表明IDCI-RF指數在實際大范圍旱情監測中具有應用潛力。