基于集成樹類算法估算農田蒸散量

顧信欽，吳立峰

（江西省南昌市南昌工程學院水利與生態工程學院，南昌 330099）

0 引 言

蒸散量（ET）是由地表水分蒸發和植物蒸騰2 部分組成，是地面水循環重要的水文過程之一。ET的精確估算有助于合理制定作物的灌溉制度，達到作物高產節約用水的目標，且對干旱監測預警有著至關重要的作用[1,2]。獲取ET的方法包括雙作物系數法、渦動相關法和機器學習方法等。

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

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

1 材料與方法

1.1 隨機森林模型

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

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

1.2 極端梯度提升模型

極端梯度提升（XGBoost）是由CHEN T Q 和GUESTRIN C[18]提出的一種新的算法。該模型基于“提升”的思想，其目的是通過將一組弱學習者的所有預測結合起來，通過特殊訓練培養出一個強的學習者。該模型具有并行計算的特點，以提高計算速度。其表示方式如下：

式中：ft(xi)為步驟t的學習值；fit和fi(t-1)是步驟為t和t-1 的預測值；xi是輸入變量。

XGBoost 算法每次分裂構建CART 樹時，為了在一定程度上降低模型的過擬合，選取使得Gain值最大的節點進行分裂。其中，Gain表示目標函數在每次分裂后信息增益，百分比越高意味著預測參數越重要，計算公式如下式。Cover 表示與此參數相關的觀測次數的覆蓋度量，由統計獲得。

式中：γ和λ是調整參數，用于防止模型過擬合；GL，GR分別為分裂時左右葉子一階梯度統計和；HL，HR為左右葉子節點二階梯度統計和的信息增益。

1.3 數據來源

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

表1 站點信息Tab.1 Site information

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

表2 不同機器學習模型的輸入參數組合Tab.2 Input parameter combinations of different machine learning models

1.4 統計指標

本研究使用了5 個統計指標，分別是均方根誤差（RMSE）、決定系數（R2）、平均絕對誤差（MAE）、偏差百分比（PBIAS）和最大10%均方根誤差（RMSEmax10%）[23]其公式分別如下：

式中：Xi,實測，Xi,預測分別為i處實測值和預測值；，分別為i處實測平均值和預測平均值；n為測試數據的總數；Emax10%為測試數據最大10%絕對誤差；z為測試數據總數的10%。

RMSE表現了模型預測值與實測值的不同，模型RMSE值越低越好。R2值越高越好，即越接近1意味這表現越好。MAE和PBIAS可以更好的反應預測誤差的實際情況，越低越好。RMSEmax10%用于評價模型在極端條件下的表現。

2 結果與分析

2.1 半小時尺度ET預報精度比較

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

表3 半小時尺度下模型預測的統計指標結果Tab.3 Statistical index results predicted by the model at the halfhour scale at the two stations

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

圖1 半小時尺度下不同模型ET預報值與實測值散點圖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預報精度比較

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

表4 日尺度下模型預測的統計指標結果Tab.4 Statistical index results predicted by the model at the daily scale at the two stations

在US-CRT 站點，在組合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 站點的RMSEmax10%值偏大，其誤差達1.4 mm/d 以上，可能使短期ET模擬誤差升高。XGB模型與RF 模型的R2值差別不大（圖2），其中除組合1 和組合2外，RF3的R2值略大于XGB3。3種組合在實測值大于4 mm/d時均出現一定低估。在XGB 模型，XGB1 的散點最接近1∶1線，其R2值為0.84，XGB2和XGB3相比于XGB1，R2值下降了14%和3%。而在RF模型中，RF2和RF3相比于RF1，R2值下降了23%和1%。

圖2 日尺度下不同模型ET預報值與實測值散點圖Fig.2 Scatter plots of ET forecast values vs.measured values of different models at the daily scale at two stations

2.3 氣象因子對ET預報貢獻度比較

為更好地認識不同站點在不同時間尺度下輸入氣象因子進行ET預報時的貢獻度情況，通過算法對氣象因子預測ET的重要性評分，繪制IT-CA2站點半小時尺度下的氣象因子與ET之間關系圖（如圖3所示）。Rn 的貢獻度最大，Gain值為0.758。除Rn外，氣象因子貢獻度從高到低依次是SWC、LAI、TA、VPD 和WS （Gain值為0.056、0.055、0.050、0.042 和0.040）。

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

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

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

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

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

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

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

3 討 論

機器學習在ET預測方面具有較高精度，但不同機器學習方法的精度有所不同。DOU X M 等報道了不同模型有其最優的適用環境，如在草地和農田站點上人工神經網絡模型的表現最好，而在森林和濕地站點上混合極端學習機模型預測ET精度最優[8]。CARTER C 等發現核嶺回歸模型可以在較小的訓練數據集上獲得最佳性能，而在數據集較大的情況下，提升樹方法最優[11]。本研究結果顯示，在半小時尺度下XGB 模型的3種組合均優于RF 模型。在日尺度下，IT-CA2站點的XGB模型也均優于RF，而在US-CRT 站點XGB3 精度略低于RF3。綜合評估2 個時間尺度預測ET的能力，XGB 模型可為準確估算ET提供參考。

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

4 結 語

本文研究了2種基于決策樹的機器學習方法XGBoost和RF對2個農田不同時間尺度ET的模擬情況，得出的結論如下：

（1）在半小時尺度下，輸入相同組合時，XGB 模型整體上優于RF 模型。輸入組合1 和組合3 時，2 個站點的2 個機器學習模型的精度比較接近，且均好于組合2。2 個站點模型的PBIAS均在5%以內，整體上不存在高估或低估現象。

（2）在日尺度下，輸入相同組合時，與半小時尺度相同，XGB 模型整體上優于RF 模型。2 個站點的2 個機器學習模型輸入組合1 和組合3 時精度相差很小。2 個站點模型的PBIAS均在5%以內，整體上不存在高估或低估現象。與IT-CA2 站點不同，在US-CRT 站點的RMSEmax10%達1.4 mm/d 以上，可能使短期ET模擬誤差升高。

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