基于隨機森林偏差校正的農業干旱遙感監測模型研究

劉 冀 張 特 魏 榕 張 茜 劉艷麗 董曉華

(1.三峽大學水利與環境學院， 宜昌 443002； 2.三峽庫區生態環境教育部工程研究中心， 宜昌 443002；3.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室， 南京 210098；4.水利部應對氣候變化研究中心， 南京 210029)

0 引言

干旱是一種普遍的自然災害，對農業、生態環境、社會經濟產生重大影響[1]。在全球氣候變暖的大背景下，未來全球干旱頻率將逐漸上升[2]。淮河流域是我國重要的農業基地，亦是我國的干旱多發區[3]。因此，針對淮河流域開展農業干旱監測模型研究對抗旱工作的開展具有重要的現實意義。

帕爾默指數(Palmer drought severity index，PDSI)、標準化降水指數(Standardized precipitation index，SPI)、標準化降水蒸散指數(Standardized precipitation evapotranspiration index，SPEI)被廣泛應用于干旱監測，并取得了較好的監測效果[4-6]，尤其是SPEI指數與土壤水分之間有較強的相關性[7-8]。干旱是一種區域性現象，這些指數通過站點觀測數據計算得出，僅能準確反映站點附近的旱情。傳統方法通過空間插值實現區域干旱監測，其結果受地理因素與插值方法等影響，存在較大的不確定性[9]。遙感技術能夠提供降水量、地表溫度以及植被長勢等空間分布信息，為實現區域農業干旱監測提供了有力支撐。一些研究嘗試應用多源遙感數據構建綜合干旱監測指數，并證實了多源遙感數據在干旱監測中的應用潛力,但這些研究采用權重方法結合多源遙感干旱因子，難以描述不同干旱因子間的非線性關系[10-12]。近年來，一些學者采用支持向量機、人工神經網絡、分類回歸樹、深度人工神經網絡等機器學習方法，融合多源遙感信息構建模型，預測旱情空間分布，為實現區域干旱監測提供了新思路[13-14]。

隨機森林(Random forest，RF)是基于集成思想的機器學習算法，較單一分類器的誤差更加穩定，具有較強的抗過擬合能力[15]。董婷等[16]采用RF、Cubist和Bagging 3種方法在黃淮地區構建了綜合干旱指數，結果表明，RF算法的擬合能力更強。PARK等[9]研究表明，相比增強回歸樹、Cubist方法，RF算法對旱情預測更加準確。RF算法基于多棵決策樹的平均結果進行預測，其結果相對準確可信，但同時也會導致一定的偏差，尤其是對極值的預測能力較弱[17]。對于干旱災害而言，極端情況往往帶來更大的損失，更應予以充分重視。相關研究表明，對RF算法進行偏差校正能有效提升預測結果的準確性[18]。

本文在綜合比較多種RF偏差校正方法適用性的基礎上，合理選取最優偏差校正方法，構建基于多源遙感數據的淮河流域農業干旱監測模型，并分析模型在農業干旱監測中的適用性，以期提高空間旱情監測的準確性，為抗旱工作的開展提供科學支撐。

1 研究區概況

淮河流域位于我國東部(111°55′～121°25′E， 30°55′～36°36′N)，地跨湖北、河南、安徽、江蘇及山東五省，是我國重要的農業基地。流域面積約270 000 km2，其中約有2/3的面積為平原和1/3的面積為丘陵山地。淮河流域地跨我國南北氣候帶，淮河以北屬暖溫帶區，淮河以南為北亞熱帶區，流域內氣候溫和，年平均降雨量約1 040 mm，年平均氣溫11～16℃，年內降雨集中，汛期降水占全年降水量的60%～80%，其獨特的氣候條件和地理位置，導致流域內的旱澇災害頻發。流域位置、高程及土地利用概況如圖1所示。

圖1 淮河流域地理位置、高程及土地利用情況Fig.1 Geographic location, elevation and land use type of Huaihe River Basin

2 數據與方法

2.1 研究數據

研究數據主要包括1970—2014年的氣象數據以及2001—2014年的遙感數據兩部分。氣象數據主要包括淮河流域內31個氣象站的逐日氣象數據，含14個農業氣象站的逐旬土壤相對濕度數據，統一處理至月尺度。遙感數據包括TRMM(Tropical rainfall measure mission)3B43降雨數據、MODIS(Moderate-resolution imaging spectroradiometer)的MOD11A3地表溫度數據和MOD13A3歸一化植被指數(Normalized vegetation index，NDVI)、STRM(Shuttle radar topography mission)高程數據。此外還包括土壤有效持水量(Available water capacity of soil，AWC)柵格數據，統一處理至1 km空間分辨率。各項數據來源信息如表1所示。

表1 研究數據來源及說明Tab.1 Sources and instructions of data

2.2 研究方法

2.2.1SPEI指數

SPEI指數采用降水與潛在蒸發的虧缺程度來量化干旱程度[19]，計算式為

(1)

Pi——降雨量,mm

PETi——潛在蒸發量,mm

k——時間尺度

i——分析時段起始時間

j——分析時段終止時間

SPEI指數可有效評估不同時間尺度的干旱情況，其中3個月尺度的SPEI指數與土壤濕度有較好的相關性[8]。通常采用Thornthwaite公式計算PET，但該公式受溫度影響較大，因此本文采用Penman-Monteith公式計算的參考蒸發量ET0代替PET，并計算1970—2014年各站點SPEI3指數。當SPEI<-2時，定義為極端干旱(極旱)；當-2.0≤SPEI<-1.5時，定義為嚴重干旱(重旱)；當-1.5≤SPEI<-1.0時，定義為中等干旱(中旱)；當SPEI≥-1.0時，定義為無旱。

2.2.2遙感指數

對各像元1、3個月的平均降水量、地表溫度、NDVI進行歸一化處理，即可計算出不同時間尺度(1、3個月)的降水狀態指數(Precipitation condition index，PCI)、地表溫度狀態指數(Temperature condition index，TCI)以及植被狀態指數(Vegetation condition index，VCI)，計算方法如表2所示。

表2 遙感指數計算方法Tab.2 Computation of different remote sensing drought monitoring indexes

2.2.3隨機森林算法

隨機森林是一種基于決策樹的集成算法，在回歸及分類問題中應用廣泛。該方法通過Bootstrap自助法在原數據集中有放回地隨機抽取n個樣本集，每次抽取2/3樣本容量的數據作為袋內數據，建立n棵決策樹構建隨機森林，并以n棵決策樹回歸結果的平均值進行預測。該方法的優點在于訓練速度相對較快且不需要進交叉驗證，同時抽樣隨機性和特征選擇隨機性使得隨機森林不易陷入過擬合[20]。隨機森林方法中的兩個重要參數分別為決策樹的數量ntree及隨機選取特征變量的數量mtry，算法詳見文獻[17]。本文采用R語言中的Random forest程序包(http:∥cran.at.r-project.org/src/contrib/Archive/randomForest/randomForest_4.6-7.tar.gz)進行計算。

2.2.4偏差校正方法

為提高RF模型精度，采用4種偏差校正方法對RF的回歸結果進行校正。

(1)簡單線性回歸法

簡單線性回歸法(Simple linear regression, SLR)主要通過建立袋內數據預測值與實際值間的線性關系式，并應用該線性關系式對袋外樣本數據的預測值進行修正，以達到偏差校正效果。計算公式為

yobs=a+bpre

(2)

式中yobs——實際值pre——預測值

a、b——系數

(2)偏差估算法

ZHANG等[21]提出的偏差估算法(Bias corrected, BC)以預測值和實際值之間的殘差為因變量，以原樣本數據及預測值為自變量，通過RF構建估計殘差的模型，在原始預測結果上進行偏差校正。計算式為

r=yobs-pre

(3)

=RFres(X,r)

(4)

cor=pre+

(5)

式中r——實際殘差——預測殘差

RFres——殘差計算函數

X——樣本數據

(3)殘差旋轉法

SONG[18]提出基于殘差旋轉(Residual rotation, RR)的偏差校正方法，能有效降低原始RF的均方誤差。該方法通過建立預測值與殘差間的線性關系，將擬合直線旋轉至與X軸重合位置，通過旋轉矩陣估計殘差對原始預測值進行校正。

(4)最優角度殘差旋轉法

最優角度殘差旋轉法(Best-angle residual rotation，BRR)基于RR法，對旋轉角進行尋優，使整體均方誤差最小。例如，將擬合曲線旋轉α可與X軸重合，在(α-β，α+β)區間尋找一個最優角度，采用該角度下的旋轉矩陣對預測值進行校正。

2.2.5模型精度評估方法

以站點SPEI3為因變量，以對應站點的PCI1、TCI1、VCI1、PCI3、TCI3、VCI3(1、3分別代表1、3個月時間尺度)、高程、AWC共8個變量為自變量，建立2001—2014年作物生長季(4—10月)的RF模型，采用上述4種偏差校正方法對RF預測結果進行校正，并選取最優的校正方法構建隨機森林偏差校正干旱監測模型(Bias-corrected random forest drought condition，BRFDC)。采用均方根誤差(Root mean square error，RMSE)最小原則確定RF模型參數mtry，并選取決定系數(Determination coefficient，R2)、RMSE以及干旱等級預測準確度評估BRFDC模型精度，綜合比較后確定最優校正方法，最后通過站點土壤相對濕度數據及流域典型干旱事件記錄對BRFDC模型的監測能力進行評估。此外，考慮到淮河流域特殊的氣候分區特征，采用C均值模糊聚類法將淮河流域劃分為兩個氣候區分別研究，計算流程如圖2所示。

圖2 計算流程圖Fig.2 Flowchart of procedures used

3 結果與分析

3.1 流域分區及RF參數優化

以各氣象站年降雨量與參考蒸發量的差值(P-ET0)為依據，將淮河流域各氣象站劃歸為北(區域1)、南(區域2)兩個分區，站點聚類結果如圖1所示，各分區氣象站特征如表3所示。由表3可知，區域1、區域2的多年平均降雨量分別為755、1 018 mm，其氣候差異主要體現在年降雨量，而年平均氣溫T和ET0的差異較小。

表3 不同區域氣候特征Tab.3 Climatic characteristics of different regions

本文構建的RF模型包含8個變量，從1～8逐一取值以確定最優參數mtry，而ntree參數統一設置為1 000。將不同參數mtry的RF模型獨立運行100次，各參數取值下的平均RMSE如圖3所示。由圖3可知，當RMSE最小時，區域1中4—10月各模型的mtry取值分別為4、7、4、3、5、4、5，區域2的mtry取值分別為5、5、7、6、5、7、5。

圖3 不同區域各參數取值對模型RMSE的影響Fig.3 RMSE of different parameter values on model of region 1 and region 2

3.2 偏差校正方法精度分析

采用SLR、BC、RR及BRR共4種方法對RF模型監測結果進行偏差校正，不同方法校正后監測結果與站點SPEI3的散點圖如圖4所示。由圖4可知，區域1中RF監測值與站點SPEI3的R2和RMSE分別為0.719和0.534，經偏差校正后的R2為0.719～0.897，RMSE為0.335～0.533；區域2中RF監測值與站點SPEI3的R2和RMSE分別為0.751和0.502，經偏差校正后的R2為0.751～0.874，RMSE為0.362～0.501。由此可見，偏差校正后較校正前精度均有一定程度上的提高，其中SLR法的校正效果并不明顯，而BC、RR以及BRR法對模型精度有顯著提升。整體而言，BRR方法的校正效果最優，校正后各區域的R2分別為0.897和0.874，RMSE分別為0.335和0.362。

圖4 不同區域各校正方法的精度評估Fig.4 Scatter diagrams of drought index estimated by different bias-correcting methods and SPEI for three-month time scale

不同校正方法對干旱等級監測的準確率如表4所示。由表4可知，RF對區域1中旱、重旱和極旱的監測準確率分別為45.1%、40.2%和0，經校正后的中旱、重旱和極旱監測準確率分別為44.4%～69.3%、39.2%～74.2%和0～50%；RF對區域2中旱、重旱和極旱的監測準確率分別為36.6%、18.6%和0，經校正后的中旱、重旱和極旱監測準確率分別為39.2%～64.1%、18.6%～43.3%、0～33.3%。整體而言，SLR法的校正效果較弱，而BC、RR及BRR法則有效提升了干旱等級監測的準確率，其中BC方法的校正效果最優，尤其是對于極端干旱監測具有顯著的優勢。綜合考慮R2、RMSE以及干旱等級監測準確率的評估結果，最終采用BC方法構建BRFDC模型。

表4 各校正方法在不同區域對各干旱等級監測準確率Tab.4 Accuracy of different drought conditions monitored by all bias-correcting methods drought in two regions %

3.3 干旱監測能力評估

為評估BRFDC模型在農業干旱監測中的適用性，采用農業氣象站相對土壤濕度及流域典型干旱事件對BRFDC模型的監測能力進行驗證。

采用BRFDC模型計算指數、站點SPEI3分別與農業氣象站10、20 cm的土壤相對濕度進行相關性分析，相關系數如表5所示。由表5可知，多數站點的SPEI3、BRFDC模型計算指數與10、20 cm的土壤濕度顯著相關(P<0.01)，二者都能有效表征土壤墑情的變化。相比SPEI3，共有10個站點10 cm土壤相對濕度與BRFDC模型計算指數相關性更強，此外，9個站點20 cm土壤相對濕度與BRFDC模型計算指數相關性更顯著。由于土壤水分是農業干旱的決定性因素，因此BRFDC模型計算的指數較SPEI3更適宜監測農業干旱。土壤水分變化受諸多因素影響，SPEI3僅采用氣象要素進行計算，在表征土壤水變化過程上略有不足，而BRFDC模型綜合考慮降水、地表溫度、植被長勢及環境因素，其計算指數更貼近于真實土壤水分的變化過程。在構建農業干旱監測模型時，需綜合考慮各種潛在影響因素。

表5 BRFDC模型模擬干旱指數與不同深度土壤相對濕度相關性分析Tab.5 Correlation analysis between drought index simulated by BRFDC and soil relative moisture in different depths

根據文獻[22]記載，2001年淮河流域發生了多年罕見的春、夏、秋連旱。河南省3—5月無有效降水，8月信陽地區創50多年降水最低紀錄；安徽省5—9月降水量較常年偏少50%～60%；江蘇省淮北地區5—6月降水量偏少，7—8月降水與往年基本持平，旱情得到緩解，9—11月降水較同期偏少81%，旱情進一步發展；山東省淮河流域汛期降水偏少，且主要集中在7—8月。

采用反距離平方(Inverse distance weighting，IDW)法對站點SPEI3進行空間插值(簡稱IDW模型)，并應用RF模型和BRFDC模型模擬了2001年5—10月淮河流域旱情空間分布，結果如圖5所示。由圖5可知，3種方法對旱情發展過程的監測結果與真實旱情基本一致，并且在空間上具有較好的一致性。IDW模型通過簡單的數學方式即可實現空間旱情監測，能大致反映干旱的時空特征，但受站點個數及插值方法限制，無站點區域的旱情監測結果存在不確定性。例如9月流域西北部出現中度旱情，由于該區域僅有開封站一個站點，而未能監測到該區域的旱情；10月流域南部旱情減輕，而該區域僅霍山站一個站點，對旱情等級有所誤判，而RF和BRFDC模型的監測結果則有效反映了這些區域的旱情。RF模型以多源遙感數據為輸入，考慮了降水、地表溫度、植被狀態、高程、AWC等因素的空間異質性，其空間監測結果更加可靠。通常極端干旱將帶來更大的損失，需要重視，但RF對6、9月流域南部以及10月流域中部極端干旱的監測效果并不理想，而BRFDC模型則準確識別了流域內的極端旱情。

圖5 淮河流域2001年5—10月IDW、RF、BRFDC模型干旱監測結果Fig.5 Drought maps of Huaihe River Basin from May to October in 2001

整體上看，BRFDC模型與土壤濕度具有較好的相關性，并且能有效反映淮河流域2001年5—10月的干旱事件的時空演變過程，對于極端旱情的監測結果也相對準確，表明該模型可用于淮河流域農業干旱監測。

3.4 討論

農業干旱涉及降水、植被、溫度以及環境因素影響，其過程非常復雜。結合多源遙感數據進行監測，能有效解釋農業干旱的復雜性。相比原始RF而言，本文構建的BRFDC模型具有更準確的干旱監測能力，對極端干旱監測的準確性從0提升到33.3%～50%，模型在淮河流域農業干旱監測中取得了較好的效果。BRFDC模型依據免費的遙感數據即可實現區域農業干旱監測，能有效減少開展干旱監測業務的成本，為缺資料地區提供可靠的旱情空間分布信息。通常降水對干旱程度起決定性作用，而本文僅采用了TRMM衛星的降水產品，其誤差勢必會對模型精度產生一定影響，提高降水產品的精度能顯著改善監測結果的準確性。由于不同的衛星降水產品在不同區域、時間的精度有所差異[23]，融合多種衛星降水產品是提高干旱監測準確性的重要途徑。此外，本文在構建模型時，考慮了降水、地表溫度、植被長勢、高程及AWC等影響因子，但農業干旱的復雜程度遠不止如此。在后續的研究中，仍有一系列其他遙感數據需要考慮，以期提高模型的監測性能，例如MODIS蒸發數據[24]、SMOS(Soil moisture and ocean salinity)土壤水分數據[25]以及土地利用類型等數據。

4 結論

(1)BRR法校正后，監測結果的R2和RMSE分別為0.897、0.874和0.335、0.362，優于其他方法；BC法校正后，預測結果的R2和RMSE分別為0.870、0.854和0.365、0.387，同時顯著提高了模型對極端干旱監測的準確率，達到33.3%～50.0%。

(2)BRFDC模型計算指數與站點實測10、20 cm的土壤相對濕度具有顯著相關性，且BRFDC模型比SPEI3更適于農業干旱監測。

(3)相比IDW空間插值法及原始RF模型，BRFDC模型模擬結果能更加準確地反映淮河流域2001年典型干旱事件的時空演變過程。