基于隨機(jī)森林偏差校正的農(nóng)業(yè)干旱遙感監(jiān)測模型研究

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

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

0 引言

干旱是一種普遍的自然災(zāi)害，對農(nóng)業(yè)、生態(tài)環(huán)境、社會經(jīng)濟(jì)產(chǎn)生重大影響[1]。在全球氣候變暖的大背景下，未來全球干旱頻率將逐漸上升[2]?；春恿饔蚴俏覈匾霓r(nóng)業(yè)基地，亦是我國的干旱多發(fā)區(qū)[3]。因此，針對淮河流域開展農(nóng)業(yè)干旱監(jiān)測模型研究對抗旱工作的開展具有重要的現(xiàn)實意義。

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

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

本文在綜合比較多種RF偏差校正方法適用性的基礎(chǔ)上，合理選取最優(yōu)偏差校正方法，構(gòu)建基于多源遙感數(shù)據(jù)的淮河流域農(nóng)業(yè)干旱監(jiān)測模型，并分析模型在農(nóng)業(yè)干旱監(jiān)測中的適用性，以期提高空間旱情監(jiān)測的準(zhǔn)確性，為抗旱工作的開展提供科學(xué)支撐。

1 研究區(qū)概況

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

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

2 數(shù)據(jù)與方法

2.1 研究數(shù)據(jù)

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

表1 研究數(shù)據(jù)來源及說明Tab.1 Sources and instructions of data

2.2 研究方法

2.2.1SPEI指數(shù)

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

(1)

Pi——降雨量,mm

PETi——潛在蒸發(fā)量,mm

k——時間尺度

i——分析時段起始時間

j——分析時段終止時間

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

2.2.2遙感指數(shù)

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

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

2.2.3隨機(jī)森林算法

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

2.2.4偏差校正方法

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

(1)簡單線性回歸法

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

yobs=a+bpre

(2)

式中yobs——實際值pre——預(yù)測值

a、b——系數(shù)

(2)偏差估算法

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

r=yobs-pre

(3)

=RFres(X,r)

(4)

cor=pre+

(5)

式中r——實際殘差——預(yù)測殘差

RFres——殘差計算函數(shù)

X——樣本數(shù)據(jù)

(3)殘差旋轉(zhuǎn)法

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

(4)最優(yōu)角度殘差旋轉(zhuǎn)法

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

2.2.5模型精度評估方法

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

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

3 結(jié)果與分析

3.1 流域分區(qū)及RF參數(shù)優(yōu)化

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

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

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

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

3.2 偏差校正方法精度分析

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

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

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

表4 各校正方法在不同區(qū)域?qū)Ω鞲珊档燃壉O(jiān)測準(zhǔn)確率Tab.4 Accuracy of different drought conditions monitored by all bias-correcting methods drought in two regions %

3.3 干旱監(jiān)測能力評估

為評估BRFDC模型在農(nóng)業(yè)干旱監(jiān)測中的適用性，采用農(nóng)業(yè)氣象站相對土壤濕度及流域典型干旱事件對BRFDC模型的監(jiān)測能力進(jìn)行驗證。

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

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

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

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

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

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

3.4 討論

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

4 結(jié)論

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

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

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