基于集成樹類算法估算農(nóng)田蒸散量

顧信欽，吳立峰

（江西省南昌市南昌工程學(xué)院水利與生態(tài)工程學(xué)院，南昌 330099）

0 引 言

蒸散量（ET）是由地表水分蒸發(fā)和植物蒸騰2 部分組成，是地面水循環(huán)重要的水文過程之一。ET的精確估算有助于合理制定作物的灌溉制度，達(dá)到作物高產(chǎn)節(jié)約用水的目標(biāo)，且對(duì)干旱監(jiān)測(cè)預(yù)警有著至關(guān)重要的作用[1,2]。獲取ET的方法包括雙作物系數(shù)法、渦動(dòng)相關(guān)法和機(jī)器學(xué)習(xí)方法等。

雙作物系數(shù)法得到ET 的前提是正確估算參考作物蒸散量（ET0），計(jì)算ET0的方法有Penman-Monteith 模型和Hargreaves模型等。其中Penman-Monteith 公式考慮因素較多，可作為ET0計(jì)算的標(biāo)準(zhǔn)方法，但對(duì)不能獲取完整氣象數(shù)據(jù)的地區(qū)并不適用[3]。Hargreaves 模型只考慮溫度與輻射參數(shù)，形式簡(jiǎn)潔，但在不同區(qū)域使用時(shí)差異較大，使用前必須予以校正[4]。渦動(dòng)相關(guān)（Eddy correlation, EC）是直接測(cè)量ET的主要方法之一。LIU S M 等報(bào)道了在中國(guó)海河流域典型下墊面用渦動(dòng)相關(guān)儀測(cè)量3年（2008-2010年）的ET，年ET介于430～560 mm[5]。LI S等利用渦動(dòng)相關(guān)儀對(duì)2007年中國(guó)西北地區(qū)覆膜條件下春玉米的ET進(jìn)行了測(cè)定，結(jié)果表明日平均ET為2.96 mm/d[6]。雖然EC是目前公認(rèn)的ET測(cè)量標(biāo)準(zhǔn)方法，能夠在大范圍內(nèi)準(zhǔn)確捕獲短期的ET信息[7]，但是觀測(cè)儀器容易受氣候影響，數(shù)據(jù)通常存在大量缺測(cè)，且設(shè)備昂貴，維護(hù)費(fèi)用高。針對(duì)這些問題，使用機(jī)器學(xué)習(xí)方法對(duì)ET的研究越來越流行，它們?cè)趶?fù)雜非線性結(jié)構(gòu)和特性的方面具有強(qiáng)大的能力。DOU X M 等發(fā)現(xiàn)4種機(jī)器學(xué)習(xí)模型預(yù)測(cè)日尺度ET時(shí)均表現(xiàn)較好，如在森林生態(tài)系統(tǒng)中決定系數(shù)（R2） 值介于0.939 8～0.959 3，納什效率系數(shù)（NSE）值介于0.887 7～0.914 7[8]。YAMA? S S 等使用了3 種模型來預(yù)測(cè)意大利南方馬鈴薯ET，發(fā)現(xiàn)當(dāng)氣象數(shù)據(jù)有限時(shí)推薦K-近鄰（KNN）模型進(jìn)行預(yù)測(cè)，反之推薦使用人工神經(jīng)網(wǎng)絡(luò)（ANN） 模型[9]。SHAN X Q 等報(bào)道了多元自適應(yīng)回歸樣條（MARS）模型預(yù)測(cè)夏玉米全生育期ET均優(yōu)于其他模型[10]。CARTER C等在不同生態(tài)類型站點(diǎn)比較了10種機(jī)器學(xué)習(xí)（ML）方法預(yù)測(cè)ET的精度和運(yùn)行時(shí)間，結(jié)果表明常綠、灌叢和草地站點(diǎn)的結(jié)果最好，濕地站點(diǎn)的結(jié)果最差[11]。ABYANEH H Z 等比較了ANN 和自適應(yīng)神經(jīng)模糊推理系統(tǒng)（ANFIS）模型的估算值與蒸滲儀測(cè)得的大蒜作物蒸散量和作物系數(shù)法的估計(jì)值，發(fā)現(xiàn)ANN和ANFIS適用于預(yù)測(cè)ET[12]。

基于決策樹的機(jī)器學(xué)習(xí)方法，在蒸散量方面得到了廣泛的應(yīng)用[13-16]。目前，關(guān)于該類方法在ET 估算方面表現(xiàn)及因子的貢獻(xiàn)方面還缺乏系統(tǒng)研究。綜上，本文采用了2種基于決策樹的機(jī)器學(xué)習(xí)方法XGBoost 和RF 對(duì)IT-CA2 和US-CRT 站點(diǎn)不同時(shí)間尺度的ET 進(jìn)行預(yù)測(cè)，輸入數(shù)據(jù)為2 種數(shù)據(jù)集，分別是地面通量站觀測(cè)數(shù)據(jù)集和ERA 再分析數(shù)據(jù)集，將氣溫（TA）、凈輻射（Rn）、飽和水汽壓差（VPD）、風(fēng)速（WS）、土壤含水量（SWC）和葉面積指數(shù)（LAI）數(shù)據(jù)按3 種不同組合作為模型輸入。研究2 種算法在3 種不同輸入組合及2 種不同時(shí)間尺度下的預(yù)測(cè)ET的能力，從而選出適宜的ET預(yù)測(cè)模型。

1 材料與方法

1.1 隨機(jī)森林模型

隨機(jī)森林（RF）是BREIMAN L[17]提出的一種基于Bagging和分類回歸樹的算法，它使用許多決策樹來創(chuàng)建森林。該方法的主要原理是：對(duì)于分類器集合h1(X)，h2(X)，…，hk(X)，從隨機(jī)向量X，Y的分布中隨機(jī)抽取訓(xùn)練數(shù)據(jù)集，邊際函數(shù)表示如下：

式中：I(.)為指示函數(shù)，這個(gè)邊際函數(shù)描述了在X，Y上分類正確的平均投票數(shù)超過分類為其他任何類別的平均票數(shù)的程度。

1.2 極端梯度提升模型

極端梯度提升（XGBoost）是由CHEN T Q 和GUESTRIN C[18]提出的一種新的算法。該模型基于“提升”的思想，其目的是通過將一組弱學(xué)習(xí)者的所有預(yù)測(cè)結(jié)合起來，通過特殊訓(xùn)練培養(yǎng)出一個(gè)強(qiáng)的學(xué)習(xí)者。該模型具有并行計(jì)算的特點(diǎn)，以提高計(jì)算速度。其表示方式如下：

式中：ft(xi)為步驟t的學(xué)習(xí)值；fit和fi(t-1)是步驟為t和t-1 的預(yù)測(cè)值；xi是輸入變量。

XGBoost 算法每次分裂構(gòu)建CART 樹時(shí)，為了在一定程度上降低模型的過擬合，選取使得Gain值最大的節(jié)點(diǎn)進(jìn)行分裂。其中，Gain表示目標(biāo)函數(shù)在每次分裂后信息增益，百分比越高意味著預(yù)測(cè)參數(shù)越重要，計(jì)算公式如下式。Cover 表示與此參數(shù)相關(guān)的觀測(cè)次數(shù)的覆蓋度量，由統(tǒng)計(jì)獲得。

式中：γ和λ是調(diào)整參數(shù)，用于防止模型過擬合；GL，GR分別為分裂時(shí)左右葉子一階梯度統(tǒng)計(jì)和；HL，HR為左右葉子節(jié)點(diǎn)二階梯度統(tǒng)計(jì)和的信息增益。

1.3 數(shù)據(jù)來源

本文選取了2 個(gè)農(nóng)田站點(diǎn)IT-CA2 和US-CRT 進(jìn)行研究，數(shù)據(jù)來源于全球通量網(wǎng)https://fluxnet.org/data/fluxnet2015-dataset/，站點(diǎn)的基本信息如表1所示。其中，IT-CA2 站點(diǎn)種植的作物為牧草[19]。US-CRT 站點(diǎn)在2年研究期內(nèi)，分別于2011年（DOY162-296） 和2012年（DOY141-275） 種植大豆，關(guān)于該數(shù)據(jù)集的詳細(xì)信息詳見CHU H S[20]、XIE J[21]和NOORMETS A[22]。

表1 站點(diǎn)信息Tab.1 Site information

模型的輸入數(shù)據(jù)使用了通量站點(diǎn)的氣象觀測(cè)數(shù)據(jù)（氣溫（TA）、凈輻射（Rn）、飽和水汽壓差（VPD）、風(fēng)速（WS）、土壤含水量（SWC））和MODIS衛(wèi)星的葉面積指數(shù)（LAI）產(chǎn)品數(shù)據(jù)(MOD15A2)以及ERA-Interim 再分析數(shù)據(jù)。為了評(píng)估不同情景下模型的表現(xiàn)，設(shè)置了3 個(gè)不同的輸入組合，組合1 包括地面通量站觀測(cè)氣象數(shù)據(jù)和MODIS 衛(wèi)星的LAI 產(chǎn)品數(shù)據(jù)，組合2 相比組合1 少了LAI 數(shù)據(jù)。組合3 與組合1 相比，其中TA、VPD和WS為ERA再分析數(shù)據(jù)（表2）。

表2 不同機(jī)器學(xué)習(xí)模型的輸入?yún)?shù)組合Tab.2 Input parameter combinations of different machine learning models

1.4 統(tǒng)計(jì)指標(biāo)

本研究使用了5 個(gè)統(tǒng)計(jì)指標(biāo)，分別是均方根誤差（RMSE）、決定系數(shù)（R2）、平均絕對(duì)誤差（MAE）、偏差百分比（PBIAS）和最大10%均方根誤差（RMSEmax10%）[23]其公式分別如下：

式中：Xi,實(shí)測(cè)，Xi,預(yù)測(cè)分別為i處實(shí)測(cè)值和預(yù)測(cè)值；，分別為i處實(shí)測(cè)平均值和預(yù)測(cè)平均值；n為測(cè)試數(shù)據(jù)的總數(shù)；Emax10%為測(cè)試數(shù)據(jù)最大10%絕對(duì)誤差；z為測(cè)試數(shù)據(jù)總數(shù)的10%。

RMSE表現(xiàn)了模型預(yù)測(cè)值與實(shí)測(cè)值的不同，模型RMSE值越低越好。R2值越高越好，即越接近1意味這表現(xiàn)越好。MAE和PBIAS可以更好的反應(yīng)預(yù)測(cè)誤差的實(shí)際情況，越低越好。RMSEmax10%用于評(píng)價(jià)模型在極端條件下的表現(xiàn)。

2 結(jié)果與分析

2.1 半小時(shí)尺度ET預(yù)報(bào)精度比較

2種方法在IT-CA2和US-CRT站點(diǎn)的半小時(shí)尺度預(yù)報(bào)精度表現(xiàn)如表3所示。在2 個(gè)站點(diǎn)，不同模型的PBIAS均在2%以內(nèi)，整體上不存在明顯高估或低估問題。在IT-CA2 站點(diǎn)，在組合1 和組合3 條件下，2 個(gè)機(jī)器學(xué)習(xí)模型的精度比較接近，且均好于組合2。相比組合2，XGB 模型的組合1 和組合3 的RMSE分別下降10%和11%，MAE分別下降了12%和11%；RF模型的組合1和組合3的RMSE分別下降9%和8%，MAE分別下降了9%和6%。從RMSEmax10%指標(biāo)可以看出，盡管誤差最大的10%點(diǎn)的誤差相比平均值有所放大，但其值仍在可接受范圍內(nèi)。從散點(diǎn)圖可見，XGB 模型與RF 模型的R2值相差很?。▓D1）。XGB3 的R2值最大為0.85，其次是XGB1（R2值為0.84）和XGB2（R2值為0.81），3 種組合均在實(shí)測(cè)值大于0.2 mm/h 時(shí)，出現(xiàn)一定低估。在RF 模型中，RF1 和RF3 比RF2 的散點(diǎn)更接近1∶1 線，同樣實(shí)測(cè)值大于0.2 mm/d 時(shí)，出現(xiàn)一點(diǎn)低估現(xiàn)象。

表3 半小時(shí)尺度下模型預(yù)測(cè)的統(tǒng)計(jì)指標(biāo)結(jié)果Tab.3 Statistical index results predicted by the model at the halfhour scale at the two stations

在US-CRT 站點(diǎn)，在相同輸入組合下，組合3的XGB 模型與RF 模型精度接近；其他組合的XGB 模型相比于RF 模型RMSE和MAE值分別提高了4%～6%和1%～3%。輸入組合2時(shí)，2 個(gè)機(jī)器學(xué)習(xí)模型表現(xiàn)最差。相比組合2，XGB 模型的組合1和組合3 的RMSE都下降2%，MAE分別下降了2%和1%；RF模型的組合1 和組合3 的RMSE和MAE都下降了2%。US-CRT站點(diǎn)的RMSEmax10%稍高于IT-CA2 站點(diǎn)，但在可接受范圍內(nèi)。從散點(diǎn)圖可見， XGB 模型與RF 模型的R2值也差別不大（圖1），在實(shí)測(cè)值大于0.6 mm/h 時(shí)3 種組合均出現(xiàn)部分低估。XGB1 與XGB3 的R2值相同為0.92，優(yōu)于XGB2（R2值為0.87）。而RF 模型中，RF3（R2值為0.92）略優(yōu)于RF1（R2值為0.91），RF2的R2值最小為0.86。

圖1 半小時(shí)尺度下不同模型ET預(yù)報(bào)值與實(shí)測(cè)值散點(diǎn)圖Fig.1 Scatter diagram of ET forecast value vs.measured value under different models in half an hour scale at two stations

2.2 日尺度ET預(yù)報(bào)精度比較

2種方法在IT-CA2和US-CRT站點(diǎn)的日尺度預(yù)報(bào)精度表現(xiàn)如表4所示。在2 個(gè)站點(diǎn)，不同模型的PBIAS均在5%以內(nèi)，在整體上與半小時(shí)尺度一樣都不存在明顯高估或低估問題。在IT-CA2 站點(diǎn)，輸入相同組合時(shí)，XGB 模型較RF 模型RMSE和MAE分別下降了5%～7%和1%～7%。2 種模型的組合3 優(yōu)于組合1，組合1 優(yōu)于組合2。相比組合2，XGB 模型的組合1 和組合3 的RMSE分別下降8%和14%，MAE分別下降了7%和14%；RF 模型的組合1 和組合3 的RMSE分別下降7%和13%，MAE分別下降了6%和9%。從RMSEmax10%指標(biāo)可以看出，2 個(gè)模型誤差RMSE在1.0 mm/d 以內(nèi)，其值不會(huì)對(duì)短期ET模擬造成影響。從圖2可見，輸入相同組合時(shí)，XGB 模型比RF 模型的散點(diǎn)更接近1∶1 線。XGB3 的R2值最大為0.75，其次是XGB1（R2值為0.71）和XGB2（R2值為0.66），3種組合均在實(shí)測(cè)值大于2 mm/d 時(shí)，出現(xiàn)低估現(xiàn)象。在RF 模型中，RF1 和RF3比RF2的散點(diǎn)更接近1∶1線，同樣實(shí)測(cè)值大于2 mm/d時(shí)，出現(xiàn)一點(diǎn)低估現(xiàn)象。

表4 日尺度下模型預(yù)測(cè)的統(tǒng)計(jì)指標(biāo)結(jié)果Tab.4 Statistical index results predicted by the model at the daily scale at the two stations

在US-CRT 站點(diǎn)，在組合1 和組合2 條件下，XGB 模型的RMSE略低于RF 模型。XGB 模型的組合1 和組合3 相比組合2的RMSE分別下降26%和22%，MAE分別下降了27%和28%；RF 模型的組合1 和組合3 相比組合2 的RMSE分別下降33%和32%，MAE都下降了26%。US-CRT 站點(diǎn)的RMSEmax10%值偏大，其誤差達(dá)1.4 mm/d 以上，可能使短期ET模擬誤差升高。XGB模型與RF 模型的R2值差別不大（圖2），其中除組合1 和組合2外，RF3的R2值略大于XGB3。3種組合在實(shí)測(cè)值大于4 mm/d時(shí)均出現(xiàn)一定低估。在XGB 模型，XGB1 的散點(diǎn)最接近1∶1線，其R2值為0.84，XGB2和XGB3相比于XGB1，R2值下降了14%和3%。而在RF模型中，RF2和RF3相比于RF1，R2值下降了23%和1%。

圖2 日尺度下不同模型ET預(yù)報(bào)值與實(shí)測(cè)值散點(diǎn)圖Fig.2 Scatter plots of ET forecast values vs.measured values of different models at the daily scale at two stations

2.3 氣象因子對(duì)ET預(yù)報(bào)貢獻(xiàn)度比較

為更好地認(rèn)識(shí)不同站點(diǎn)在不同時(shí)間尺度下輸入氣象因子進(jìn)行ET預(yù)報(bào)時(shí)的貢獻(xiàn)度情況，通過算法對(duì)氣象因子預(yù)測(cè)ET的重要性評(píng)分，繪制IT-CA2站點(diǎn)半小時(shí)尺度下的氣象因子與ET之間關(guān)系圖（如圖3所示）。Rn 的貢獻(xiàn)度最大，Gain值為0.758。除Rn外，氣象因子貢獻(xiàn)度從高到低依次是SWC、LAI、TA、VPD 和WS （Gain值為0.056、0.055、0.050、0.042 和0.040）。

圖3 IT-CA2站點(diǎn)半小時(shí)尺度下的輸入因子與ET之間關(guān)系Fig.3 Relationship between input factors and ET at the half-hour scale of IT-CA2 stations

IT-CA2 站點(diǎn)日尺度下的氣象因子與ET之間關(guān)系如圖4所示。Rn 的貢獻(xiàn)度最大，Gain值為0.610。除Rn 外，氣象因子貢獻(xiàn)度從高到低依次是TA、WS、VPD、SWC 和LAI（Gain值為0.099、0.096、0.077、0.062 和0.055）。與半小時(shí)尺度的貢獻(xiàn)度排序相比，日尺度因子貢獻(xiàn)有所不同，但貢獻(xiàn)最高的仍是Rn。

圖4 IT-CA2站點(diǎn)日尺度下的輸入因子與ET之間關(guān)系Fig.4 Relationship between input factors and ET at daily scale of ITCA2 stations

US-CRT 站點(diǎn)半小時(shí)尺度下的氣象因子與ET之間關(guān)系如圖5所示。氣象因子貢獻(xiàn)度從高到低依次是Rn、TA、WS、VPD、SWC 和LAI （Gain值為0.688、0.167、0.054、0.043、0.024和0.023）。

圖5 US-CRT站點(diǎn)半小時(shí)尺度下的輸入因子與ET之間關(guān)系Fig.5 Relationship between input factors and ET at the half-hour scale of US-CRT stations

US-CRT站點(diǎn)日尺度下的氣象因子與ET之間關(guān)系如圖6所示。氣象因子貢獻(xiàn)度排序與日尺度稍有不同，從高到低依次是TA、Rn、LAI、SWC、VPD 和WS（Gain值為0.455、0.278、0.147、0.055、0.038和0.026）。

圖6 US-CRT站點(diǎn)日尺度下的輸入因子與ET之間關(guān)系Fig.6 Relationship between input factors and ET at daily scale of US-CRT stations

3 討 論

機(jī)器學(xué)習(xí)在ET預(yù)測(cè)方面具有較高精度，但不同機(jī)器學(xué)習(xí)方法的精度有所不同。DOU X M 等報(bào)道了不同模型有其最優(yōu)的適用環(huán)境，如在草地和農(nóng)田站點(diǎn)上人工神經(jīng)網(wǎng)絡(luò)模型的表現(xiàn)最好，而在森林和濕地站點(diǎn)上混合極端學(xué)習(xí)機(jī)模型預(yù)測(cè)ET精度最優(yōu)[8]。CARTER C 等發(fā)現(xiàn)核嶺回歸模型可以在較小的訓(xùn)練數(shù)據(jù)集上獲得最佳性能，而在數(shù)據(jù)集較大的情況下，提升樹方法最優(yōu)[11]。本研究結(jié)果顯示，在半小時(shí)尺度下XGB 模型的3種組合均優(yōu)于RF 模型。在日尺度下，IT-CA2站點(diǎn)的XGB模型也均優(yōu)于RF，而在US-CRT 站點(diǎn)XGB3 精度略低于RF3。綜合評(píng)估2 個(gè)時(shí)間尺度預(yù)測(cè)ET的能力，XGB 模型可為準(zhǔn)確估算ET提供參考。

ERA 再分析數(shù)據(jù)的精度在國(guó)內(nèi)外得到了廣泛的研究，PAREDES P 等報(bào)道了ERA 再分析數(shù)據(jù)計(jì)算的ET0在未經(jīng)過偏差矯正時(shí)與Penman-Monteith 公式計(jì)算的ET0有很強(qiáng)的相關(guān)性（R2＞0.80），但存在一定的高估[24]。SRIVASTAVA P K 等發(fā)現(xiàn)ERA 再分析中期數(shù)據(jù)和中分辨率成像光譜儀（MODIS）得出的ET0比美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心（NCEP）估算的所有季節(jié)ET0更準(zhǔn)確[25]。本研究結(jié)果顯示，在輸入組合1 和組合3 時(shí)，2 種模型的精度都十分接近。可見ERA 再分析數(shù)據(jù)與地面通量站觀測(cè)的數(shù)據(jù)的差異對(duì)模型模擬沒有明顯影響，鑒于ERA 數(shù)據(jù)在時(shí)間尺度和空間尺度上的優(yōu)勢(shì)，該數(shù)據(jù)集具有廣泛應(yīng)用于ET估算的潛力。

4 結(jié) 語(yǔ)

本文研究了2種基于決策樹的機(jī)器學(xué)習(xí)方法XGBoost和RF對(duì)2個(gè)農(nóng)田不同時(shí)間尺度ET的模擬情況，得出的結(jié)論如下：

（1）在半小時(shí)尺度下，輸入相同組合時(shí)，XGB 模型整體上優(yōu)于RF 模型。輸入組合1 和組合3 時(shí)，2 個(gè)站點(diǎn)的2 個(gè)機(jī)器學(xué)習(xí)模型的精度比較接近，且均好于組合2。2 個(gè)站點(diǎn)模型的PBIAS均在5%以內(nèi)，整體上不存在高估或低估現(xiàn)象。

（2）在日尺度下，輸入相同組合時(shí)，與半小時(shí)尺度相同，XGB 模型整體上優(yōu)于RF 模型。2 個(gè)站點(diǎn)的2 個(gè)機(jī)器學(xué)習(xí)模型輸入組合1 和組合3 時(shí)精度相差很小。2 個(gè)站點(diǎn)模型的PBIAS均在5%以內(nèi)，整體上不存在高估或低估現(xiàn)象。與IT-CA2 站點(diǎn)不同，在US-CRT 站點(diǎn)的RMSEmax10%達(dá)1.4 mm/d 以上，可能使短期ET模擬誤差升高。

（3）在IT-CA2 站點(diǎn)半小時(shí)尺度和日尺度的因子貢獻(xiàn)度排序不同，但貢獻(xiàn)度最大的因子均為Rn，其Gain值分別是0.758和0.610。在US-CRT 站點(diǎn)上2 個(gè)尺度的因子貢獻(xiàn)度排序也不同，半小時(shí)尺度從高到低為Rn、TA、WS、VPD、SWC 和LAI，日尺度從高到低TA、Rn、LAI、SWC、VPD和WS。