氣象要素缺失條件下不同機(jī)器學(xué)習(xí)模型計算參考作物蒸散量比較

郝林如，郭向紅，2，雷 濤，鄭利劍，馬娟娟，孫西歡，蘇媛媛，胡飛鵬

（1.太原理工大學(xué)水利科學(xué)與工程學(xué)院，太原 030024；2.中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038）

0 引 言

參考作物蒸散量（ET0）是指導(dǎo)農(nóng)業(yè)灌溉的重要參數(shù)，它被定義為具有特定特征的假設(shè)作物參考表面的蒸散量，將ET0與不同作物系數(shù)相乘可得不同作物需水量[1]。因此ET0的準(zhǔn)確計算對灌溉規(guī)劃和農(nóng)業(yè)用水管理具有重大意義[2]。

Penman-Monteith（PM）公式是聯(lián)合國糧農(nóng)組織推薦計算ET0的標(biāo)準(zhǔn)方法[3]，所需氣象要素包括最高氣溫、最低氣溫、相對濕度、2 m 高風(fēng)速和日照時數(shù)[4]。但是在很多地區(qū)，由于降雨、儀器故障、人為誤操作等原因會導(dǎo)致氣象要素缺失，無法使用PM 公式計算[5]。因此，很多學(xué)者提出了基于較少氣象要素的簡化模型，例如Hargreaves-Samani 模型、Priestly-Taylor 模型和Irmak-Allen 模型等，但是這些模型的模擬精度不高[6]。隨著計算機(jī)技術(shù)和機(jī)器學(xué)習(xí)算法的發(fā)展，越來越多的機(jī)器學(xué)習(xí)算法應(yīng)用于氣象要素缺失條件下ET0的預(yù)測，該方法簡便高效，易于推廣。WEN等[7]在氣象要素缺失條件下，使用支持向量機(jī)、人工神經(jīng)網(wǎng)絡(luò)和經(jīng)驗法Priestley-Taylor、Hargreaves-Samani 對蒸散量進(jìn)行預(yù)測，得出支持向量機(jī)在蒸散量的估計方面比人工神經(jīng)網(wǎng)絡(luò)和經(jīng)驗法更準(zhǔn)確的結(jié)果。程才華[8]在26個氣象要素缺失條件下，使用支持向量機(jī)對蒸散量進(jìn)行模擬，得出模擬精度最高的氣象要素組合為氣溫、風(fēng)速、日照時數(shù)、相對濕度。張薇等[9]在4 個氣象要素缺失組合下，使用3種基于樹型算法，包括梯度提升決策樹(GBDT)、隨機(jī)森林(RF)和回歸樹(CART)模型預(yù)測蒸散量，結(jié)果表明GBDT 和RF有較高模擬精度。

然而盡管國內(nèi)外在氣象要素缺失條件下使用支持向量機(jī)、隨機(jī)森林、人工神經(jīng)網(wǎng)絡(luò)等機(jī)器學(xué)習(xí)模型預(yù)測ET0已經(jīng)有了大量的研究，但是在不同氣象要素缺失組合下系統(tǒng)比較不同機(jī)器學(xué)習(xí)模型的研究還較少，因此本文以山西省果樹研究所Adcon-Ws 無線自動氣象站為例，把PM 公式的計算值作為標(biāo)準(zhǔn)值，用決策樹（CART）、隨機(jī)森林（RF）、梯度提升決策樹(GBDT)、極端梯度提升樹（XGBoost）、支持向量機(jī)（SVR）、BP 神經(jīng)網(wǎng)絡(luò)（BPNN）和深度學(xué)習(xí)（DL）7 種機(jī)器學(xué)習(xí)模型模擬不同氣象要素缺失條件下的ET0，為氣象要素缺乏的地區(qū)預(yù)測蒸散量提供理論依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

研究站點位于山西省晉中市太谷縣的農(nóng)科院果樹研究所（東經(jīng)112°32’,北緯37°23’,海拔約800 m），所在區(qū)域為典型的大陸性半干旱氣候，多年平均氣溫為9.8 ℃，年平均降雨量約460 mm，無霜期175 d。受季風(fēng)季影響，年內(nèi)降雨嚴(yán)重不均[10]。

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

本文所用基礎(chǔ)氣象數(shù)據(jù)來自山西省果樹研究所Adcon-Ws無線自動氣象監(jiān)測站，包括最高氣溫（Tmax）、最低氣溫（Tmin）、2 m高風(fēng)速（u2）、相對濕度（RH）、日照時數(shù)（n）。

1.3 研究方法

（1）Penman Monteith 模型。PM 公式綜合考慮各種氣象要素，具有較強(qiáng)的普適性，其計算結(jié)果可作為模擬ET0的標(biāo)準(zhǔn)值[3]。公式如下：

式中：ET0為參考作物蒸散量，mm/d；Rn為凈輻射，MJ/（m2·d）；G為土壤熱通量，MJ/（m2·d）；Tmean為2 m 高處日平均氣溫，℃；u2為2 m高處的風(fēng)速，m/s；es為飽和水汽壓，kPa；ea為實際水汽壓，kPa；Δ為飽和水汽壓曲線的斜率；γ為濕度計常數(shù)，kPa/℃。

（2）決策樹（CART）。CART 是一種廣泛應(yīng)用的決策樹學(xué)習(xí)方法，用于回歸時稱為CART回歸樹。它包含特征選擇、樹的生成與剪枝。該模型計算復(fù)雜度不高，解釋性強(qiáng)，對中間缺失值不敏感，但方差高，容易過擬合[9]。

（3）隨機(jī)森林算法（RF）。隨機(jī)森林是一種基于Bagging的集成學(xué)習(xí)算法。首先用bootstrap方法有放回的從數(shù)據(jù)集中重復(fù)取樣來選取訓(xùn)練樣本，在得到訓(xùn)練集后分別獨立訓(xùn)練決策樹，然后將其組合，由決策樹預(yù)測結(jié)果的均值作為最終結(jié)果[11]。RF 算法數(shù)據(jù)采樣以及特征采樣隨機(jī),并且它訓(xùn)練速度快，容易做成并行化方法，實現(xiàn)比較簡單，能通過減小模型方差提高性能，被廣泛用來解決分類與回歸問題[11]。

（4）梯度提升決策樹(GBDT)算法。GBDT 是一種基于Boosting 的集成學(xué)習(xí)算法，它的基本思想是用當(dāng)前模型損失函數(shù)的負(fù)梯度信息來訓(xùn)練新加入的弱分類器，然后將訓(xùn)練好的弱分類器以累加的形式結(jié)合到現(xiàn)有的模型中，弱分類器一般選擇CART 模型[12]。GBDT 的性能在RF 的基礎(chǔ)上又有一步提升，它能靈活的處理各種類型的數(shù)據(jù)，在相對較少的調(diào)參時間下，預(yù)測的準(zhǔn)確度較高，并且對數(shù)據(jù)缺失不敏感，能通過減小模型偏差提高性能。但由于它是Boosting，基學(xué)習(xí)器存在串行關(guān)系，難以并行訓(xùn)練數(shù)據(jù)[9]。

（5）極端梯度提升（XGBoost）算法。XGBoost 是一種高效的GBDT 算法。在算法的損失函數(shù)上，XGBoost 加上了正則化部分。在算法的優(yōu)化方式上，GBDT 的損失函數(shù)只對誤差部分做負(fù)梯度（一階泰勒）展開，而XGBoost的損失函數(shù)對誤差部分做二階泰勒展開，使得結(jié)果更加準(zhǔn)確[13]。同時XGBoost 具有模塊結(jié)構(gòu)可以支持并行計算[14]。

（6）支持向量機(jī)（SVR）。支持向量機(jī)算法用于回歸時稱為支持向量回歸（SVR），它比較好的實現(xiàn)了結(jié)構(gòu)風(fēng)險最小化思想，得到的是全局最優(yōu)解[5]。為解決非線性回歸問題，模型通常用轉(zhuǎn)換函數(shù)φ將變量x映射到高維特征空間，再對其進(jìn)行計算，并引入核函數(shù)可以將原來線性算法非線性化[15]。SVR能較好解決小樣本以及非線性問題。本研究核函數(shù)類型選擇為徑向基函數(shù)（RBF）。

（7）BP 神經(jīng)網(wǎng)絡(luò)（BPNN）。BP 神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)過程由信號的正向傳播和誤差的反向傳播組成，要素信息由輸入層輸入，經(jīng)隱含層處理后，傳向輸出層，若輸出層的實際輸出和期望輸出不符，就進(jìn)入誤差的反向傳播階段，此過程不斷重復(fù)直至達(dá)到要求[16]。Kolmogorov 定理表明，具有一個隱含層的3 層BPNN 已經(jīng)有足夠的逼近能力[17]。因此本文采用3 層網(wǎng)絡(luò)結(jié)構(gòu)，分別由輸入層、隱含層和輸出層組成，每層的激活函數(shù)均設(shè)為Relu函數(shù)，學(xué)習(xí)率為0.01。

（8）深度學(xué)習(xí)（DL）。深度學(xué)習(xí)是在傳統(tǒng)人工神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上添加更多的隱含層，以建立更加復(fù)雜的網(wǎng)絡(luò)結(jié)構(gòu)[5]。本文采用增加BP 神經(jīng)網(wǎng)絡(luò)隱含層數(shù)量的深度學(xué)習(xí)算法，把隱含層設(shè)為3 層，它的具體訓(xùn)練過程與BP 神經(jīng)網(wǎng)絡(luò)相同。每層激活函數(shù)設(shè)為Relu 函數(shù)，優(yōu)化方法采用Adam 算法，學(xué)習(xí)率為0.01。為了防止過擬合，正則化系數(shù)L1、L2設(shè)為0.01。

（9）經(jīng)驗法。選取4種經(jīng)驗法與機(jī)器學(xué)習(xí)模型進(jìn)行對比分析，具體計算公式見表1[18]。

表1 經(jīng)驗法及計算公式Tab.1 Empirical method and calculation formula

1.4 氣象要素缺失條件下不同機(jī)器學(xué)習(xí)預(yù)測模型構(gòu)建

首先，選取與PM 公式計算所需的氣象要素最高氣溫（Tmax）、最低氣溫（Tmin）、相對濕度（RH）、2 m 高風(fēng)速（u2）和日照時數(shù)（n）一致的數(shù)據(jù)構(gòu)成組合1。為了構(gòu)建氣象要素缺失條件下機(jī)器學(xué)習(xí)模型，在組合1的基礎(chǔ)上分別缺失相對濕度（RH）、日照時數(shù)（n）、2 m 高風(fēng)速（u2）、最高最低氣溫（Tmax、Tmin）構(gòu)成組合2～5。將氣象要素缺失最多，僅包含Tmax、Tmin數(shù)據(jù)組合構(gòu)成組合9。再在組合9 的基礎(chǔ)上分別增加了RH、n、u2構(gòu)成組合6～8。氣象要素缺失組合見表2。將9組氣象要素組合分別作為7 種機(jī)器學(xué)習(xí)模型的輸入，ET0作為模型的輸出，構(gòu)建CART1～CART9，RF1～RF9，GBDT1～GBDT9，XGBoost1～XGBoost9，SVR1～SVR9，BPNN1～BPNN9，DL1～DL9共63個機(jī)器學(xué)習(xí)模型。

表2 不同氣象要素缺失組合Tab.2 Combinations of different input meteorological factors

CART、RF、GBDT、XGBoost、SVR 算法是基于Python 語言的Skicit-Learn 模塊建立的模型，BPNN 和DL 是基于Python語言Keras 模塊建立的模型。在選取的687 組樣本數(shù)據(jù)中，隨機(jī)選取其中515組（約占總樣本的3/4）作為訓(xùn)練樣本，172組（約占總樣本的1/4）作為測試樣本。使用StandardScaler 對數(shù)據(jù)集進(jìn)行標(biāo)準(zhǔn)化和歸一化處理。CART、RF、GBDT、XGBoost、SVR 5 種模型使用網(wǎng)格搜索和5 折交叉驗證相結(jié)合的方法選擇最優(yōu)參數(shù)組合。網(wǎng)格搜索其本質(zhì)就是將參數(shù)空間劃分成若干網(wǎng)格，通過遍歷網(wǎng)格交叉點處所有的參數(shù)組合來優(yōu)化需要訓(xùn)練的模型[19]；5 折交叉驗證是將初始數(shù)據(jù)集分成5份，每次將其中一份保留，其他4 份用于訓(xùn)練，以最終5 折交叉驗證所得的平均誤差指標(biāo)作為模型最終結(jié)果[20]。BPNN和DL模型的參數(shù)有隱含層節(jié)點數(shù)和最大迭代次數(shù)，本研究在根據(jù)經(jīng)驗確定出隱含層節(jié)點數(shù)的范圍的基礎(chǔ)上，將其放大，通過不斷調(diào)整節(jié)點數(shù)和迭代次數(shù)，使得訓(xùn)練目標(biāo)方差達(dá)到要求，最終確定最優(yōu)參數(shù)組合。

1.5 模型評價指標(biāo)

選用平均絕對誤差（MAE）、均方根誤差(RMSE)、平均絕對百分比誤差（MAPE）和決定系數(shù)（R2）作為模型評價指標(biāo)，公式如下[21]：

若計算精度越高，MAE、RMSE和MAPE越接近0，R2越接近于1。

2 結(jié)果與分析

2.1 氣象要素缺失條件下不同機(jī)器學(xué)習(xí)模型模擬精度對比

為了對比氣象要素缺失條件下不同機(jī)器學(xué)習(xí)模型模擬精度，本文選用訓(xùn)練誤差和預(yù)測誤差對比來說明模型的泛化能力，同時選用預(yù)測誤差結(jié)果來說明模型的模擬精度。各模型在不同氣象要素缺失組合下的ET0計算的精度對比見表3。

表3 7種模型在不同氣象要素缺失組合下的訓(xùn)練誤差和測試誤差Tab.3 Training and testing errors of seven models under different combinations of meteorological elements

續(xù)表

當(dāng)輸入全部5個氣象要素時，所有模型的測試集R2達(dá)到最高，MAE、RMSE和MAPE達(dá)到最低。7 種機(jī)器學(xué)習(xí)模型中DL的模擬精度最高，R2達(dá)到0.995，MAE、RMSE和MAPE分別為0.015 mm/d、0.018 mm/d 和5.124%，其他6 種模型按模擬精度由高到低依次為BPNN、SVR、GBDT、XGBoost、RF 和CART。同時DL 的訓(xùn)練誤差與預(yù)測誤差相差較小，R2差值為0.002，MAE和RMSE的相差范圍均不超過0.005 mm/d，MAPE的差異值在1.5%以內(nèi)，泛化能力較好，其他6 種模型按泛化能力由高到低依次為BPNN、SVR、XGBoost、GBDT、CART和RF。

當(dāng)缺乏一個氣象要素時（組合2～4），機(jī)器學(xué)習(xí)模型模擬精度相比組合1有所降低，在缺失相同氣象要素條件下，不同機(jī)器學(xué)習(xí)模型模擬精度也有所差異。組合2～4中，DL和BPNN模擬精度較高，排名為第1 或第2，其他5 種模型模擬精度在組合2、4 中排名相同，由高到低依次為SVR、GBDT、XGBoost、RF、CART,組合3 中排名依次為XGBoost、RF、GBDT、SVR、CART。在這3 個組合中，BPNN 訓(xùn)練誤差與預(yù)測誤差相差均較小，泛化能力較好，其他6種模型泛化能力排名在組合2 中由高到低依次為DL、SVR、XGBoost、GBDT、RF、CART，組合3 中依次為DL、SVR、GBDT、XGBoost、CART、RF，組合4 中依次為DL、SVR、XGBoost、GBDT、CART、RF。

當(dāng)缺乏2 個氣象要素（組合5～8）時，機(jī)器學(xué)習(xí)模型模擬精度都進(jìn)一步降低。在組合5～8 中，DL 的預(yù)測誤差均最小，模擬精度排名第1。其他6種模型模擬精度在組合5～6，8中排名相同，由高到低依次為BPNN、SVR、XGBoost、GBDT、RF、CART，組合7 中排名依次為BPNN、SVR、GBDT、XGBoost、RF、CART。在這4 個組合中，BPNN 訓(xùn)練誤差與預(yù)測誤差相差均較小，泛化能力排名第1。其他6 種模型泛化能力在組合5、7 中排名相同，由高到低依次為DL、SVR、XGBoost、GBDT、CART、RF，在組合6、8 中排名相同，由高到低依次為DL、SVR、GBDT、XGBoost、CART、RF。

當(dāng)輸入2個氣象要素（組合9），所有機(jī)器學(xué)習(xí)模型模擬精度達(dá)到最低，即模型測試集的R2均小于0.8，MAE、RMSE除BPNN 和DL 外均大于0.5 mm/d，MAPE均大于25%。7 種機(jī)器學(xué)習(xí)模型中DL 的模擬精度最高，其他6 種模型按模擬精度由高到低依次為BPNN、SVR、XGBoost、GBDT、CART、RF。同時BPNN 的訓(xùn)練誤差與預(yù)測誤差相差較小，泛化能力較好，其他6 種模型按泛化能力由高到低依次為DL、SVR、CART、GBDT、XGBoost、RF。

綜上，在所有組合中DL 模擬精度普遍較高，BPNN 泛化能力普遍較好，表明DL、BPNN 相比其他模型可以有效避免參數(shù)的隨機(jī)性選擇，具有高度自學(xué)習(xí)自適應(yīng)能力，能夠快速有效地反映多因素之間高度非線性的復(fù)雜關(guān)系[22]，提高了模擬精度和泛化能力。同時DL 相比BPNN 建立了更復(fù)雜的網(wǎng)絡(luò)結(jié)構(gòu)使得模擬結(jié)果更加準(zhǔn)確，但可能因為訓(xùn)練樣本數(shù)量或噪聲的影響導(dǎo)致訓(xùn)練過程出現(xiàn)過擬合現(xiàn)象，從而泛化能力稍遜于BPNN。劉堃等[5]在干旱半干旱區(qū)生態(tài)系統(tǒng)缺失蒸散量插補(bǔ)研究中得出DL和BPNN 相比SVR、CART等模型預(yù)測效果更好的結(jié)論，與本研究結(jié)論一致。XGBoost、GBDT、RF 是以決策樹為基學(xué)習(xí)器的集成模型，模擬精度和泛化能力相比CART都有了一定的提高。但這3種模型與SVR相比訓(xùn)練誤差較小或相差不多，而預(yù)測誤差卻大于SVR，這可能是因為訓(xùn)練的樣本數(shù)量不足或數(shù)據(jù)有噪聲，導(dǎo)致XGBoost、GBDT、RF 在訓(xùn)練過程中容易出現(xiàn)過擬合現(xiàn)象，從而降低了預(yù)測精度。本研究中除組合3 以外其他組合XGBoost 和GBDT 模擬精度和泛化能力均大于RF，說明XGBoost、GBDT 在預(yù)測蒸散量時相比RF 更有優(yōu)勢，這與牛曼麗等[23]、張薇等[9]的研究結(jié)果基本一致。本文不同組合中XGBoost、GBDT 模型的模擬精度和泛化能力排名不同，但整體上XGBoost 效果更好，這與理論上XGBoost 是一種高效的GBDT算法一致。

2.2 不同氣象要素對機(jī)器學(xué)習(xí)模型模擬的影響

為了分析不同氣象要素對機(jī)器學(xué)習(xí)模型模擬的影響，以隨機(jī)森林（RF）為例將氣象要素缺失條件下的模型模擬結(jié)果進(jìn)行對比，其他6種模型情況與之相同。

組合2～5 分別缺失氣象要素相對濕度（RH）、日照時數(shù)（n）、2 m 高風(fēng)速（u2）和最高最低氣溫（Tmax、Tmin），這4 種組合下隨機(jī)森林模型（RF2～RF5）的模擬精度由圖1可得，RF3、RF5 的模擬精度稍遜于RF2、RF4，說明最高最低氣溫（Tmax、Tmin）和日照時數(shù)（n）對于模型模擬蒸散量的影響比其他氣象要素更顯著，同時RF2的模擬精度低于RF4，RF3的模擬精度低于RF5，說明日照時數(shù)（n）對模型的影響最大，然后是最高最低氣溫（Tmax、Tmin）、相對濕度（RH），最后是2 m高風(fēng)速（u2）。

圖1 RF2～RF5與PM模型計算結(jié)果對比Fig.1 Comparison of calculation results between RF2～RF5 and PM model

組合6～8 是在組合9（Tmax、Tmin）的基礎(chǔ)上分別增加了RH、n、u2，這3 種組合下隨機(jī)森林模型（RF6～RF8）的模擬精度見圖2，RF7 的MAE和RMSE較小，模擬精度較高，其次是RF6，最后是RF8，因此RH、n、u2對模型預(yù)測ET0的影響程度為n＞RH＞u2。

圖2 RF6～RF9與PM模型計算結(jié)果對比Fig.2 Comparison of calculation results between RF6～RF9 and PM model

綜上，氣象要素對機(jī)器學(xué)習(xí)模型模擬參考作物蒸散量（ET0）的影響程度由大到小依次為n、Tmax、Tmin、RH、u2。

2.3 機(jī)器學(xué)習(xí)模型與經(jīng)驗法比較

為了比較機(jī)器學(xué)習(xí)模型和經(jīng)驗法的模擬精度，選取機(jī)器學(xué)習(xí)模型中模擬精度較差的決策樹模型與H-S、I-A、Makkink、P-T 經(jīng)驗法進(jìn)行對比。組合9 與H-S 的輸入?yún)?shù)相同，組合7 與I-A、Makkink、P-T 的輸入?yún)?shù)相同，故將這2種組合下決策樹算法（CART7、CART9）與對應(yīng)經(jīng)驗法的MAE、RMSE進(jìn)行比較，見圖3。

圖3 決策樹模型與經(jīng)驗法精度比較Fig.3 Comparison of accuracy between decision tree model and empirical method

當(dāng)Tmax、Tmin（組合9）作為輸入?yún)?shù)，H-S 經(jīng)驗法的MAE和RMSE大于CART9；當(dāng)Tmax、Tmin、n（組合7）作為輸入?yún)?shù)，P-T 經(jīng)驗法的精度高于I-A、Makkink 但MAE和RMSE仍大于CART7。因此在輸入氣象要素相同的條件下，決策樹模型的模擬精度均大于經(jīng)驗法，由此說明機(jī)器學(xué)習(xí)模型模擬精度更高，故建議在氣象要素缺失條件下優(yōu)先考慮使用機(jī)器學(xué)習(xí)模型。

3 結(jié) 論

本文構(gòu)建了9 種氣象要素缺失組合下的7 種機(jī)器學(xué)習(xí)ET0計算模型進(jìn)行模擬精度和泛化能力的對比分析，并與經(jīng)驗法比較，結(jié)果表明：

（1）在9 種氣象要素缺失組合下，深度學(xué)習(xí)（DL）和BP神經(jīng)網(wǎng)絡(luò)（BPNN）的模擬精度均較高且有較好的泛化能力，其他模型在不同氣象要素缺失組合中模擬精度和泛化能力有不同的排名，但整體效果較好的是支持向量機(jī)（SVR）。

（2）在氣象要素缺失條件下，各氣象要素對模型模擬參考作物蒸散量（ET0）的影響程度不同，由大到小排序依次為日照時數(shù)、最高最低氣溫、相對濕度、2 m 高風(fēng)速（n、Tmax、Tmin、RH、u2）。

（3）與H-S、I-A、Makkink、P-T 經(jīng)驗法相比，機(jī)器學(xué)習(xí)模型的模擬精度均大于輸入相同組合的經(jīng)驗法，在氣象要素缺失條件下可優(yōu)先考慮使用機(jī)器學(xué)習(xí)模型。