顧及時空特征的參考作物蒸散量集成學(xué)習(xí)估算

劉傲,趙東保*,魏義長,肖煉

(1. 華北水利水電大學(xué)測繪與地理信息學(xué)院,河南 鄭州 450046; 2. 自然資源部四川基礎(chǔ)地理信息中心,四川 成都 610041)

參考作物蒸散量(reference crop evapotranspiration,ET0)是田間水量平衡法中計算作物需水量的關(guān)鍵參數(shù),也是反映土壤水分供應(yīng)能力的重要指標.隨著道爾頓蒸發(fā)定律的提出,學(xué)者們構(gòu)建了多種ET0經(jīng)驗估算公式,如Penman-Monteith公式、Blaney-Criddle公式、Thornthwaite公式、Hargreaves-Samani公式與Priestley-Taylor公式[1-3].這些公式皆是基于ET0與氣象要素在特定的地理條件上的關(guān)系發(fā)展而來[4],在不同的氣候區(qū)這種關(guān)系可能會發(fā)生改變.例如,Penman-Monteith公式在寒冷半干旱與熱帶草原氣候區(qū)域精度較差[5],Hargreaves-Samani公式在氣候濕潤的區(qū)域精度欠佳[6].并且,丘陵等復(fù)雜地形上的ET0也會受到坡度、坡向等地形因子的影響,在這類地理條件復(fù)雜的區(qū)域,將經(jīng)驗公式與地形因子相結(jié)合,可獲得更精準的估算結(jié)果[7].但空間要素、氣象要素與ET0之間關(guān)系十分復(fù)雜,使用數(shù)學(xué)模型去描述這種非線性關(guān)系仍具有局限性.

計算機硬件的提升與機器學(xué)習(xí)技術(shù)的發(fā)展為描述復(fù)雜的非線性關(guān)系提供了一種新的解決思路.許多研究[8-9]均表明,相比經(jīng)驗公式,機器學(xué)習(xí)方法得出的ET0估算結(jié)果與實測值更相近.現(xiàn)有研究大多是通過對比尋求ET0估算精度更好的單一模型,但單一模型只能發(fā)現(xiàn)氣候要素與ET0之間的部分關(guān)系,在地形復(fù)雜、氣候多變的區(qū)域泛化能力欠缺.近些年集成學(xué)習(xí)技術(shù)的飛速發(fā)展,為ET0的估算帶來了新的解決方案.多個機器學(xué)習(xí)模型進行組合構(gòu)建新模型即為集成學(xué)習(xí)技術(shù),大致可分為Bagging,Boosting與Stacking模型.其中,Boosting模型與隨機森林等Bagging模型已被證實能有效提升ET0的估算精度[10-11].與這兩者相比,盡管Stacking模型在ET0估算研究中起步較晚,但也取得了初步成效,WU等[12]基于Stacking模型估算中國不同氣候區(qū)的ET0,指出Stacking模型可以顯著提升ET0估算精度;FAN等[13]發(fā)現(xiàn)在改變ET0估算區(qū)域時,Stacking模型的性能變異性最低.但上述研究在使用Stacking模型估算ET0時僅考慮了氣象要素,這會導(dǎo)致模型在地理條件復(fù)雜的地區(qū)性能下降.

四川省是中國地形最為復(fù)雜的省份,最高峰與最低點落差近7 000 m,氣候差異明顯[14],結(jié)合空間特征去分析ET0與氣候要素之間的關(guān)系會更加可靠.另外,氣候要素在時序上也有一定的規(guī)律[15],可結(jié)合時序特征對ET0估算結(jié)果進行修正.因此,文中在顧及時空特征的情況下構(gòu)建Stacking模型,對四川省氣候要素與ET0之間的關(guān)系進行模擬,并綜合評價在考慮不同特征時,Stacking模型與其他機器學(xué)習(xí)模型的估算精度,最終確定最佳輸入特征及模型,為四川省區(qū)域灌溉決策、農(nóng)業(yè)水資源管理提供一定的理論依據(jù).

1 方法與資料

1.1 研究區(qū)概況

四川省是中國重要的農(nóng)業(yè)區(qū),雖然四川省水資源總量豐富,但時空分布不均,形成季節(jié)性和區(qū)域性缺水.四川省各氣象站點分布如圖1所示,圖中h為海拔.

圖1 四川省數(shù)字高程模型及氣象站位置Fig.1 Digital elevation model and weather station location in Sichuan Province

1.2 數(shù)據(jù)說明

研究所用氣象數(shù)據(jù)均來自國家氣象數(shù)據(jù)中心,使用四川省52個國家基本站1981—2010年的逐日氣象數(shù)據(jù).研究所用(digital elevation model, DEM)數(shù)據(jù)分辨率為30 m,來源于中國科學(xué)院計算機網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺.坡度和坡向數(shù)據(jù)皆基于DEM數(shù)據(jù)提取而來.

1.3 時空自相關(guān)分析方法

1.3.1空間自相關(guān)

空間自相關(guān)是衡量觀測數(shù)據(jù)在空間中是否聚集的標準,能反映研究區(qū)域上相鄰站點ET0的平均關(guān)聯(lián)程度.文中采用全局莫蘭指數(shù)I表示ET0空間自相關(guān)系數(shù).標準差倍數(shù)Z用來檢驗空間自相關(guān)是否存在.置信水平為95%時,如果Z>1.96,可認為I值通過檢驗.

(1)

(2)

1.3.2時間自相關(guān)

時間自相關(guān)是衡量時間序列中不同時段的觀測數(shù)據(jù)間的相關(guān)性,能反映研究區(qū)域ET0在時間上的變化規(guī)律.

(3)

1.4 ET0經(jīng)驗公式估算方法

彭曼公式(FAO 56 Penman-Monteith)常被用作ET0估算方法精度檢驗[17].其公式為

(4)

式中:Rn為作物表面凈輻射,MJ/(m2·d);G為土壤熱通量,MJ/(m2·d);γ為干濕球濕度計常數(shù),kPa/℃;Tmean為距地面高2 m處日平均氣溫;u2為距地面高2 m處風(fēng)速,m/s;es為飽和水汽壓,kPa;ea為實際水汽壓,kPa;Δ為飽和水汽壓與溫度曲線的斜率,kPa/℃.

1.5 精度評定

采用決定系數(shù)(R2)、平均絕對值誤差(MAE)和均方誤差(MSE)來評價所有方法的ET0估算精度.R2為模型估算值與實測值的擬合度,MSE為估算值中誤差較大值的偏離程度.MAE和MSE越趨于0,R2越趨于1,說明該模型估算精度越高.

(5)

(6)

(7)

式中:ypredi為ET0估算值;yi為ET0實測值;l為驗證數(shù)據(jù)樣本數(shù)量;i為樣本序數(shù),i=1,2,3,…,l;yavg為所有驗證數(shù)據(jù)ET0實測值的平均值.

2 Stacking模型ET0估算方法

Stacking模型一般為2層結(jié)構(gòu).第1層為基模型,在給定的訓(xùn)練集上進行訓(xùn)練,將輸出結(jié)果作為新的特征供第2層的元模型學(xué)習(xí),得到最終的估算結(jié)果.相比Bagging與Boosting模型,Stacking模型的精度與泛化能力更高[18].

XGBoost,LightGBM,GBDT、極限樹與隨機森林均為較流行的集成學(xué)習(xí)模型,圍繞這些學(xué)習(xí)模型構(gòu)建Stacking模型,以獲取更高的精度與泛化能力.顧及時空特征的Stacking模型ET0估算方法構(gòu)建過程如下:

1) 對研究對象四川省ET0分別進行空間與時間上的自相關(guān)分析,當研究結(jié)果呈正相關(guān)關(guān)系時,可引入時空特征,與1981—2010年四川省氣象數(shù)據(jù)組成原始數(shù)據(jù)集.

2) 確定Stacking模型的基模型與元模型.

3) 對原始數(shù)據(jù)集構(gòu)建特征工程.由于氣象觀測儀器會出現(xiàn)損壞與維護等情況,氣象數(shù)據(jù)中會包含異常值,需對氣象數(shù)據(jù)進行清洗.此外,初始選用的時空與氣象特征數(shù)據(jù)并非都有益于提升ET0估算精度,需對這些特征進行評估,從而剔除無用特征.將處理過的數(shù)據(jù)集作為構(gòu)建Stacking模型所用數(shù)據(jù)集.

4) 將數(shù)據(jù)集按需劃分為訓(xùn)練集與測試集.在引入空間特征的情況下,應(yīng)按照地理位置劃分數(shù)據(jù)集去檢驗?zāi)Ｐ途?

5) 將訓(xùn)練集劃分為訓(xùn)練集1與測試集1,利用五折交叉驗證訓(xùn)練Stacking基模型.具體過程如圖2所示.

圖2 Stacking基模型構(gòu)建過程Fig.2 Stacking base model construction process

6) 對剩余所有基模型重復(fù)步驟5的操作,并將每個基模型對ET0的估算值都視作元模型的一個特征進行建模.

7) 使用Stacking模型對最初劃分的測試集進行ET0估算,并對其進行精度評價.

3 試驗結(jié)果與分析

首先得到本研究區(qū)域ET0空間和時序特征自相關(guān)分析結(jié)果,考察ET0隨空間和時序的變化情況.之后根據(jù)基模型與特征的選取結(jié)果構(gòu)建模型,進而得到顧及時空特征的Stacking模型的試驗結(jié)果和精度分析情況.

3.1 研究區(qū)ET0時空自相關(guān)

3.1.1研究區(qū)ET0空間自相關(guān)分析結(jié)果

基于2006年四川省各氣象站點月均ET0進行空間自相關(guān)分析如表1所示.

表1 2006年四川省月均ET0全局莫蘭指數(shù)ITab.1 Monthly average ET0 global Moran index I in Sichuan Province in 2006

由表1可知,1—12月I值均大于0,且Z值皆大于1.96.說明在置信水平為95%的情況下,四川省ET0呈現(xiàn)空間正相關(guān)關(guān)系.因此,將空間特征用于模型構(gòu)建,具體選用經(jīng)度、緯度、海拔、坡度、坡向這5項空間特征.

3.1.2研究區(qū)ET0時間自相關(guān)

基于四川省臺站號為56038氣象站1981—2010年逐日ET0進行時間自相關(guān)分析.圖3為ET0隨時間的變化圖,由圖可知,ET0始終在0～6 mm/d以年為單位上下波動,說明該時序數(shù)據(jù)為明顯的平穩(wěn)序列.這意味著四川省歷史ET0具備代表性與可延續(xù)性,在未來一段時期內(nèi)仍遵循歷史變化趨勢.并且,平穩(wěn)序列有助于模型準確捕捉數(shù)據(jù)的潛在趨勢,提升模型的穩(wěn)定性與準確性.

圖3 ET0隨時間變化圖Fig.3 ET0 changes with time

圖4為ET0時間自相關(guān)系數(shù)隨時間間隔k變化圖,由圖可知,當k值在365,730 d等年的倍數(shù)周圍時,ACF趨于1.0,說明ET0具有以年為單位的周期性.根據(jù)分析結(jié)果,將時序特征用于構(gòu)建模型,具體選用年份、日序及歷史的3天共5項時序特征.

圖4 ET0時間自相關(guān)系數(shù)隨時間間隔k變化圖Fig.4 ET0 time autocorrelation coefficient changes with time interval k

3.2 Stacking模型的基模型選擇與特征篩選

3.2.1Stacking模型的基模型選擇

將臺站號為56038的氣象站1981—2010年逐日ET0作為數(shù)據(jù)集,基于隨機森林、極限樹、GBDT,XGBoost,LightGBM創(chuàng)建不同的基模型組合來構(gòu)建Stacking模型.

表2為不同基模型組合的精度評定與時間成本,表中t1為訓(xùn)練時間.當基模型為隨機森林與極限樹2個Bagging模型時,訓(xùn)練時間成本為25 916 ms.當基模型為GBDT, XGBoost, LightGBM這3個Boosting模型時,訓(xùn)練時間成本為6 018 ms,僅為前者的1/4.可見,不同基模型的訓(xùn)練時長主要取決于基模型類型,并非基模型數(shù)量越多,時間成本就越高.基模型數(shù)量增多雖有助于提升精度,但過多的基模型意味著過高的訓(xùn)練時間成本,且精度的提升也有限.當基模型組合為5個時(隨機森林、極限樹、GBDT,XGBoost,LightGBM),比基模型組合為2個(隨機森林,極限樹)的時間成本約增加了22%,但MSE和MAE減少了約12%和5%.當在5個基模型上又添加強模型Adaboost時,其各項精度評定指標僅略微增加,但訓(xùn)練時間成本更高,故文中采用隨機森林、極限樹、GBDT,XGBoost,LightGBM等5個模型的組合作為最終的基模型.

表2 不同基模型組合的精度評定與時間成本Tab.2 Accuracy estimation and time cost of different base model combinations

3.2.2數(shù)據(jù)集特征篩選

數(shù)據(jù)集中用來估算ET0的初始特征為最大溫度Tmax、最小溫度Tmin、日序Rd、日照時長N、平均溫度Tmean、相對濕度RH、氣壓p、緯度φ1、海拔h、經(jīng)度φ2、坡向A、平均風(fēng)速ua、最大風(fēng)速umax、坡度C、年份ta以及前3天的ET0,該值分別用ET(t-1),ET(t-2),ET(t-3)表示.但是這些特征并非都有助于提升ET0估算精度,故需要對這些特征進行篩選.文中基于XGBoost評估這些特征的重要度排序如圖5所示.圖中D為特征的重要度.

XGBoost的基模型為決策樹,能夠根據(jù)結(jié)構(gòu)分數(shù)的增益情況計算出選擇哪個特征作為分割點,每個特征的重要度就是其在所有樹中出現(xiàn)的次數(shù)之和.重要度相對較高的特征會被更多地用于模型中構(gòu)建決策樹.之后通過遞歸特征消除算法(recursive feature elimination, RFE)消除噪聲特征,即按照特征重要度從大到小排名,每次增加一個特征,反復(fù)構(gòu)建隨機森林模型.選用決定系數(shù)R2作為模型的評價分數(shù),當評價分數(shù)開始下降時,即可將新增特征視為噪聲.

圖6為RFE特征篩選.圖中R2為模型評價分數(shù).由圖可知,模型評價分數(shù)R2從添加平均風(fēng)速ua后開始下降,并在之后一直保持下降趨勢.從模型角度考慮,在四川省地區(qū),平均風(fēng)速ua、最大風(fēng)速umax、坡度C、年份ta這些特征并不能幫助提升模型性能,故剔除這些特征.

圖6 RFE特征篩選Fig.6 RFE characteristics screening

3.3 Stacking模型ET0估算精度驗證和分析

分別驗證空間和時間特征對Stacking模型的ET0估算精度改善情況.首先,在氣象特征與時序特征的基礎(chǔ)上引入空間特征,分別訓(xùn)練顧及空間特征模型和沒有顧及空間特征模型,并對四川省52個站點進行五折交叉精度驗證;其次,在氣象特征與空間特征的基礎(chǔ)上引入時間特征,將四川省1981—2010年的逐日ET0數(shù)據(jù)集按照8∶2分成訓(xùn)練集與測試集,分別訓(xùn)練顧及時序特征模型和沒有顧及時序特征模型,并探討時序特征對模型性能的改善效果;最后,在顧及時空特征的情況下,將Stacking模型與經(jīng)驗?zāi)Ｐ鸵约案鞣N基模型(隨機森林、極限樹、LightGBM,XGBoost)進行全面的精度對比,以驗證文中模型的各方面性能.

3.3.1空間特征對ET0估算的影響

圖7為顧及空間特征模型與沒有顧及空間特征模型精度對比.如圖所示,顧及空間特征的模型在訓(xùn)練集上的R2與測試集幾乎一致,平均MSE與MAE分別降低了68%和46%,而沒有顧及空間特征的模型在訓(xùn)練集上的平均R2相較于測試集提升了4%,平均MSE與MAE分別下降了92%和73%.沒有顧及空間特征的模型與顧及空間特征的模型在訓(xùn)練集上精度無明顯變化,相較于沒有顧及空間特征的模型,顧及空間特征的模型在測試集上平均R2提升了3%,平均MSE降低了76%,平均MAE降低了51%.基模型中,XGBoost與LightGBM在加入了空間特征后精度提升最為明顯,其中,XGBoost的R2提升了5%,MSE和MAE分別降低了85%和62%;LightGBM的R2提升了4%,MSE和MAE分別降低了84%和59%.

四川省各氣象要素存在極強的空間異質(zhì)性.在缺少空間特征的情況下,模型無法準確地描述四川省ET0與各特征之間內(nèi)在關(guān)系,使得模型泛化能力受限,在測試集的精度大幅降低.相比于Bagging模型,XGBoost與LightGBM這類Boosting模型往往偏差更低,但在訓(xùn)練過程中也更容易過擬合,導(dǎo)致模型在測試集的精度明顯低于訓(xùn)練集.而在加入空間特征后,XGBoost與LightGBM在測試集的精度與訓(xùn)練集相比無明顯下降,說明顧及空間特征可改善XGBoost與LightGBM的過擬合問題,從而提升了Stacking模型的性能.

3.3.2時序特征對ET0估算的影響

圖8為顧及時序特征模型與沒有顧及時序特征模型精度對比.由圖可知,所有模型在測試集的精度和訓(xùn)練集相比皆無明顯降低.

圖8 顧及時序特征模型與沒有顧及時序特征模型精度對比Fig.8 Comparison of accuracy of models that take into account temporal characteristics and models that do not take into account temporal characteristics

在測試集中,相較于沒有顧及時序特征的模型,顧及時序特征的模型的平均R2提升了4%,平均MSE與MAE降低了92%和72%.其中,在引入時序特征后,GBDT的R2提升了5%,MSE與MAE分別降低了94%和75%;XGBoost的R2提升了4%,MSE與MAE分別降低了93%和74%;LightGBM的R2提升了4%,MSE與MAE分別降低了93%和75%,精度提升程度高;隨機森林在引入時序特征后的R2提升了3%,MSE與MAE分別降低了85%和62%,精度提升程度最小.

在引入時序特征后,所有模型在訓(xùn)練集和測試集上的精度都得到了明顯提升.說明時序特征可協(xié)助模型挖掘四川省地區(qū)ET0周期性變化規(guī)律,從而提升四川省地區(qū)ET0的估算精度.其中,LightGBM,XGBoost與GBDT精度提升程度顯著高于隨機森林與極限樹,這得益于Boosting這類模型在訓(xùn)練過程中會不斷調(diào)整未能充分利用時序特征的基模型權(quán)重.而Bagging模型的每個基模型都是獨立的,訓(xùn)練過程中無法自主調(diào)整每個基模型的權(quán)重,使得利用時序特征程度不同的基模型對整個Bagging模型的貢獻都相同,降低了時序特征帶來的收益.

3.3.3Stacking與其基模型及經(jīng)驗?zāi)Ｐ途葘Ρ?/p>

將1981—2005年四川省逐日氣象數(shù)據(jù)作為訓(xùn)練集,顧及時空特征構(gòu)建Stacking模型.將2006—2010年數(shù)據(jù)作為測試集,Stacking模型與其各個基模型與彭曼公式ET0估算精度對比如表3所示.

表3 Stacking及其基模型與彭曼公式ET0估算精度對比Tab.3 Comparison of estimation accuracy between Stacking and its base model with the Penman-Monteith formula ET0

由表3可知,Stacking模型的平均R2相較于彭曼公式提升了39%,平均MSE和MAE分別降低了95%和77%.Stacking模型精度顯著優(yōu)于彭曼公式,一方面由于彭曼公式作為一種經(jīng)驗公式,具備一定的地域性,但Stacking模型同時結(jié)合時空特征與氣象特征可探索出更通用的ET0估算方法.另一方面由于ET0與氣象特征之間存在著復(fù)雜的非線性關(guān)系,很難使用數(shù)學(xué)模型準確地描述他們之間內(nèi)在聯(lián)系.文中則采用多種機器學(xué)習(xí)模型集成的方式,可適應(yīng)不同的數(shù)據(jù)特征,從而更好地處理非線性關(guān)系,以此提升ET0估算精度.此外,在精度驗證中,2006—2010年逐年精度最優(yōu)基模型分別為極限樹,LightGBM,極限樹,LightGBM,XGBoost,但Stacking模型在逐年精度驗證中皆優(yōu)于其各個基模型,并有著較小的誤差波動范圍,說明模型同時具備較高的精度與穩(wěn)定性.在逐年驗證中,Stacking模型精度并未隨著年份的遞增而下降,說明其在時間上也具備較高的泛化能力.

4 結(jié) 論

文中以四川省52個氣象站點為研究對象,揭示了四川省ET0值的時空分布變異情況,提出了顧及時空特征的Stacking模型ET0估算方法,并將其與經(jīng)驗?zāi)Ｐ团砺揭约案鱾€基模型包括XGBoost,LightGBM、隨機森林、GBDT、極限樹等進行了全面的精度對比和驗證.得出主要結(jié)論如下:

四川省ET0在空間與時間上的分布都具有明顯的自相關(guān)性.時空特征的引入,可顯著提升模型在整個四川省地區(qū)的泛化能力與估算精度.在顧及時空特征的情況下,Stacking模型精度顯著優(yōu)于彭曼公式.此外,雖然文中所用數(shù)據(jù)較老,但是模型精度在2006—2010年的逐年驗證中并未出現(xiàn)隨時間推移而下降的情況,說明其具備較好的泛化能力,對當前時段也有一定的參考意義.綜上所述,顧及時空特征的Stacking模型可有效提升ET0,并給其他地理條件復(fù)雜的地區(qū)ET0估算方法提供一定的參考價值.