顧及時空特征的參考作物蒸散量集成學習估算

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

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

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

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

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

1 方法與資料

1.1 研究區概況

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

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

1.2 數據說明

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

1.3 時空自相關分析方法

1.3.1空間自相關

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

(1)

(2)

1.3.2時間自相關

時間自相關是衡量時間序列中不同時段的觀測數據間的相關性,能反映研究區域ET0在時間上的變化規律.

(3)

1.4 ET0經驗公式估算方法

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

(4)

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

1.5 精度評定

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

(5)

(6)

(7)

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

2 Stacking模型ET0估算方法

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

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

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

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

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

4) 將數據集按需劃分為訓練集與測試集.在引入空間特征的情況下,應按照地理位置劃分數據集去檢驗模型精度.

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

圖2 Stacking基模型構建過程Fig.2 Stacking base model construction process

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

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

3 試驗結果與分析

首先得到本研究區域ET0空間和時序特征自相關分析結果,考察ET0隨空間和時序的變化情況.之后根據基模型與特征的選取結果構建模型,進而得到顧及時空特征的Stacking模型的試驗結果和精度分析情況.

3.1 研究區ET0時空自相關

3.1.1研究區ET0空間自相關分析結果

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

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

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

3.1.2研究區ET0時間自相關

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

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

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

圖4 ET0時間自相關系數隨時間間隔k變化圖Fig.4 ET0 time autocorrelation coefficient changes with time interval k

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

3.2.1Stacking模型的基模型選擇

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

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

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

3.2.2數據集特征篩選

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

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

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

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

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

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

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

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

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

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

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

圖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%,精度提升程度最小.

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

3.3.3Stacking與其基模型及經驗模型精度對比

將1981—2005年四川省逐日氣象數據作為訓練集,顧及時空特征構建Stacking模型.將2006—2010年數據作為測試集,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模型精度顯著優于彭曼公式,一方面由于彭曼公式作為一種經驗公式,具備一定的地域性,但Stacking模型同時結合時空特征與氣象特征可探索出更通用的ET0估算方法.另一方面由于ET0與氣象特征之間存在著復雜的非線性關系,很難使用數學模型準確地描述他們之間內在聯系.文中則采用多種機器學習模型集成的方式,可適應不同的數據特征,從而更好地處理非線性關系,以此提升ET0估算精度.此外,在精度驗證中,2006—2010年逐年精度最優基模型分別為極限樹,LightGBM,極限樹,LightGBM,XGBoost,但Stacking模型在逐年精度驗證中皆優于其各個基模型,并有著較小的誤差波動范圍,說明模型同時具備較高的精度與穩定性.在逐年驗證中,Stacking模型精度并未隨著年份的遞增而下降,說明其在時間上也具備較高的泛化能力.

4 結 論

文中以四川省52個氣象站點為研究對象,揭示了四川省ET0值的時空分布變異情況,提出了顧及時空特征的Stacking模型ET0估算方法,并將其與經驗模型彭曼公式以及各個基模型包括XGBoost,LightGBM、隨機森林、GBDT、極限樹等進行了全面的精度對比和驗證.得出主要結論如下:

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