基于數值天氣預報后處理的參考作物蒸散量預報改進

姚付啟 董建華 范軍亮 曾文治 吳立峰

(1.魯東大學水利工程學院， 煙臺 264010； 2.南昌工程學院水利與生態工程學院， 南昌 330099；3.西北農林科技大學水利與建筑工程學院， 陜西楊凌 712100；4.武漢大學水資源與水電工程科學國家重點實驗室， 武漢 430072)

0 引言

準確的蒸散量預報對作物需水量、驅動水文和作物生長模型評估以及農業和水資源管理決策具有重要意義。蒸散量可以通過蒸滲儀、渦動相關和波文比等系統觀測，但這些設備通常價格昂貴，且維護成本很高[1]。作物的實際蒸散量更多采用參考作物蒸散量(ET0)函數進行推導[2]。

Penman-Monteith方程(以下簡稱FAO56 PM)[3]是國際公認的計算ET0的標準方法。FAO56 PM方程是一種基于物理的方法，結合了生理和空氣動力學參數[2]，并不需要參數本地化校正。然而，FAO56 PM方程要求太陽輻射、氣溫、風速和相對濕度等氣象資料齊全。

近年來，學者們不斷嘗試使用公共天氣預報信息進行未來ET0的預測。CAI等[4]將公共天氣預報的多云情況、每日最高和最低溫度以及風速等級轉換為FAO56 PM公式所需的信息，獲得了理想的ET0預報精度。LUO等[5]使用1～7 d公共天氣預報和Hargreaves-Samani模型進行ET0預報，其均方根誤差(RMSE)介于0.87～1.36 mm/d，相關系數介于0.64～0.86。CHANG等[6]使用改用Thornthwaite模型基于公共天氣預報的氣溫數據對ET0進行了預報。YANG等[7]研究發現，天氣預報風速對ET0估算精度低于多年的日平均風速。

進入21世紀以來，數值天氣預報(Numerical weather prediction，NWP)模型在高空探測影響天氣過程、氣候變化和全球尺度天氣異常方面取得了巨大的進展[8]。如同基于衛星和雷達觀測的數據同化技術[9-10]，NWP的進步是計算機能力增強、模型的時空高分辨率以及預報技術和工具改進的直接結果。NWP基于流體動力學方程通過確定性或概率的方法來描述大氣運動，需要明確的邊界條件。初始條件的偏差、方程的簡化或某些因素的低預測性都可能造成天氣預報失靈。因此，預報的結果需要對逐個變量和逐個站點進行檢查。通常使用地面氣象站點數據對NWP的網格值進行偏差校正或后處理[11]。MEDINA等[12]采用歐洲中尺度天氣預報系統(EC)、英國氣象局預報系統(BMO)和美國全球集合預報系統(GEFS)集合預報的方法來改進單個NWP模型預報能力的不足。

機器學習算法由于具有極強的非線性擬合能力，已被應用于氣象數據的空間校正、ET0預測預報等領域。LUO等[13]將ANN和天氣預報相結合，預測我國不同地區的ET0。YIN等[14]評估了4種神經網絡方法在預報ET0中的表現，發現多層感知器的效果最好。TRAORE等[15]使用雙向長短記憶神經網絡對ET0進行了7 d預報。LU等[16]使用蝙蝠算法耦合XGBoost對未來1～3個月的ET0進行了預測。

目前，基于NWP對ET0預報的預報期多為7 d以內，從灌溉決策的角度，迫切需要預報期為15 d以上的可靠方法。以往研究多為基于單氣象要素的偏差校正，但NWP輸出氣象因子的預報結果與相應觀測值是否匹配尚不明確，其他預報因子是否蘊藏著提高ET0預報性能的重要信息等相關研究還未見報道。因此，本研究以第二代全球集合預報系統在我國西北地區的輸出產品為例，以1～16 d為預報期，基于LightGBM機器學習模型，檢驗NWP氣象預報因子與相應觀測值的匹配程度，探索使用多輸出結果提高單個氣象要素預報精度的可能性，以提高ET0的預報能力。

1 材料與方法

1.1 研究區域概況

西北干旱區位于中國內陸西北部，屬大陸性內陸氣候區。氣候干燥，年降水量主要分布在200～400 mm，水資源比較緊缺。天然植被的分布以典型干草原、荒漠草原、荒漠為主。年平均氣溫5～10℃。四季分明，晝夜溫差大，全年日照達3 000 h，無霜期170 d左右，是全國太陽輻射資源最充足的地區之一，特別適宜農作物及瓜果生長。

1.2 數據來源

本研究數據來源有2個，其中氣象站點數據來自中國氣象數據共享網(http:∥www.cma.cn)，選取了中國西北地區9個站點(圖1)作為研究對象，氣象數據包括：每日2 m高處最高和最低氣溫(Tmax和Tmin)、地表總輻射(Rs)、相對濕度(RH)和2 m高處風速(U2)(表1)。數值天氣預報數據來自NOAA的第二代全球集合預報系統(GEFSv2)[17]，1～8 d和9～16 d的空間分辨率分別為0.5°×0.5°和1°×1°。數據包括未來1～16 d的2 m最高和最低氣溫、向下短波輻射、比濕度、U和V方向風速，時間分辨率為4 h和6 h，其中比濕度轉換為相對濕度，U和V方向風速經加權平均獲得平均風速。2個數據集數據分為2部分， 2005—2014年用于模型建模(訓練)，2015—2019年用于檢驗模型。

表1 本研究所選氣象站點概況

1.3 研究方法

1.3.1反距離權重法

反距離權重(IDW)法基于兩個假設:未知值附近的點影響比遠的控制點影響更大；影響程度(權重)與彼此之間的距離的倒數成指數關系[18]，其計算式為

(1)

式中Z——被插值點的值

zi——控制點的值

本研究使用氣象站點周圍NWP的4個控制點信息并使用IDW法進行插值。

1.3.2等距累積分布函數法

累積分布函數匹配(CDFm)方法[19]是一種廣泛應用的偏差校正方法，它需要假設模型和觀測數據的累積分布一致。為了克服這一問題，LI等[20]提出了等距累積分布函數匹配(EDCDFm)法，該方法對模型和觀測的累積分布曲線進行了偏差的校正。EDCDFm定義為

(2)

式中F——累積頻率函數

F-1——累積頻率函數的反函數

Xm′——校正后的未來模型預測值

Xm——數值天氣預報的未來模型預測值

下角o表示歷史訓練期觀測值，s表示歷史訓練期模擬值，s′表示未來模型模擬值。

本研究中，首先使用IDW方法將NWP數據插值到氣象站點，然后使用3次函數對氣象因子的累積頻率曲線進行擬合。

1.3.3LightGBM法

LightGBM是傳統GBDT的一種高效改進算法，擅長處理高維數據，能提高計算效率，同時保證高模型精度[21]。傳統的決策樹算法采用層次策略對樹進行生長，對同一層的葉子不加區分地進行處理，造成了不必要的計算成本(圖2)。為了降低訓練數據的維數，該算法的決策樹采用葉向生長策略。每一次從所有的葉子中，找到一個最大的分裂增益，然后分裂完成一個循環[22]。為了避免在樣本量不足時過擬合，有必要增加樹的最大深度限制。LightGBM還可以使用Gini得分對輸入特征的重要性進行排序，在建模時剔除干擾因子。本研究中，LightGBM分別基于兩種方式進行建模(表2)，第1種為使用單因子對單因子的方式對氣象因子進行偏差校正，記為M2方法，該方法每個站點訓練樣本數為4×3 652=14 608個，檢驗樣本數為4×1 826=7 304個；第2種同時使用5個氣象因子依次對其中一個氣象因子進行偏差校正，記為M3方法。使用預報期1～16 d分別作為輸入，觀測值作為輸出來建立模型，M2方法的輸入特征數量是4站點×6個時刻(第4天開始4個時刻)，共24(16)個，M3方法的輸入特征數量是4站點×6個時刻(第4天開始4個時刻)×5要素，共120(80)個。本研究中所需要調參的參數有：最大樹深度、學習率和樹的個數，調參使用網格搜索法[23]。

表2 本研究所使用偏差校正方法

1.4 FAO-56 Penman-Monteith 模型

采用FAO-56 PM法[23]計算西北地區9個站點的ET0，計算公式為

(3)

式中ET0——參考作物蒸散量，mm/d

Rn——地表凈輻射，MJ/(m2·d)

G——土壤熱通量密度，MJ/(m2·d)

Tmean——2 m高處的平均氣溫，℃

es、ea——飽和與實際水汽壓，kPa

Δ——蒸汽壓曲線的斜率，kPa/℃

γ——溫度計常數，kPa/℃

1.5 評估指標

采用均方根誤差(RMSE)、決定系數(R2)和平均絕對誤差(MAE)來評估氣象因子誤差。

2 結果與分析

2.1 氣象因子預報精度比較

2.1.1氣溫

氣溫是影響飽和水汽壓和地面長波輻射收支的重要因子，從而間接影響了ET0的變化。3種方法在未來16 d的Tmax預報精度如圖3a～3c所示。隨著預報天數的增加，Tmax的預報精度呈下降趨勢。M1、M2和M3方法在9站點平均RMSE介于2.57～5.03℃、2.31～4.82℃和2.01～4.36℃，MAE介于1.90～3.89℃、1.71～3.75℃和1.48～3.40℃;R2介于0.87～0.97、0.87～0.97和0.89～0.98。由此可以看出，3種方法的精度由高到低依次是M3、M2和M1。

Tmin的預報精度也隨預報期的延長而下降(圖3d～3f)，但總體上誤差小于Tmax，其M1、M2和M3方法在9站點平均RMSE介于 2.11～4.03℃、2.65～4.10℃和2.26～3.76℃，MAE介于1.90～3.41℃、2.28～3.39℃和1.94～3.06℃，R2介于0.88～0.97、0.87～0.95和0.90～0.96。盡管M1方法在站點6、9等表現優于M2方法，但是在站點1、2、3精度明顯低于M2方法，特別是在預報期9～16 d精度較差。總體上，M3方法優于M1和M2方法，其在預報期前8 d精度稍高于M1方法，但是9～16 d明顯優于M1方法。

2.1.2地表總輻射

輻射是蒸散的最主要能量來源，其直接主導地表氣溫變化。3種方法在不同站點預報Rs的表現如圖3g～3i所示。3種方法在預報Rs的前8 d精度逐漸下降，在第9天開始由于網格尺度變大造成了Rs明顯下降，其后精度又逐漸下降。其中，M1、M2和M3方法在9站點平均RMSE介于3.44～5.23 MJ/(m2·d)、3.53～5.22 MJ/(m2·d)和3.22～4.91 MJ/(m2·d)，MAE介于2.52～3.76 MJ/(m2·d)、2.57～3.83 MJ/(m2·d)和2.37～3.64 MJ/(m2·d)，R2介于0.82～0.59、0.81～0.57和0.85～0.61(從前期到后期)?？傮w上，M3方法表現優于M1和M2方法，而M1和M2方法誤差比較接近。此外，3種方法在站點1、2、9的誤差較大，特別是9～16 d，RMSE超過了5 MJ/(m2·d)，R2也在0.5以下。

2.1.3風速

風速是影響空氣阻力的重要因子。3種方法在不同站點預報U2表現如圖3j～3l所示。3種方法在預報U2不同站點呈現明顯差異。在站點2、3、6、7，3種方法在預報期1～16 d都呈現較好的表現。而在站點4總體表現較差。其他站點前8 d比后8 d的精度高50%以上。其中，M1、M2和M3方法在9個站點平均RMSE介于0.67～0.87 m/s、0.53～0.69 m/s和0.49～0.64 m/s，MAE介于0.48～0.62 m/s、0.39～0.51 m/s和0.36～0.47 m/s，R2介于0.03～0.24、0.04～0.37和0.10～0.43。

2.1.4相對濕度

相對濕度(RH)是估算飽和水汽壓虧缺的重要因子。3種方法在不同站點預報RH表現如圖3m～3o所示。3種方法在預報RH的前8 d精度快速下降，在第9天開始由于網格尺度變大造成了精度急劇下降，部分站點的精度已經接近隨機分布，可以看出GEFSv2總體上預報RH的能力比較一般。其中，M1、M2和M3方法在9站點平均RMSE介于15.0%～17.2%、11.7%～14.4%和10.0%～13.1%，MAE介于12.0%～13.5%、9.0%～11.4%和7.6%～10.2%，R2介于0.19～0.42、0.25～0.53和0.36～0.66?？梢钥闯觯琈3方法優于M2和M1方法，而M2方法優于M1方法。

2.2 ET0預報精度比較

3種方法在不同站點預報ET0表現如圖3p～3r所示。3種方法在預報ET0誤差隨預報期延長呈上升趨勢，其中，M1、M2和M3方法在9站點平均RMSE介于0.66～0.93 mm/d、0.57～0.83 mm/d和0.53～0.79 mm/d，MAE介于0.44～0.61 mm/d、0.38～0.56 mm/d和0.35～0.53 mm/d，R2介于0.82～0.91、0.84～0.93和0.86～0.94?？梢钥闯?，3種方法的精度從高到低依次是M3、M2和M1。從站點看，站點3、7的精度最高，1～16 d的誤差都沒有超過0.8 mm/d，而站點4、9的誤差較大，特別是9～16 d，其RMSE已經超過了1 mm/d。

3種方法在站點5不同預報時段FAO56 PM估算值與預報值的散點圖如圖4所示。M1方法在預報期1～16 d的R2從0.91下降到0.82，從第4天起存在著嚴重的高估問題，特別是PM 估算值大于3 mm/d時，有大量點預報值達到了8 mm/d以上。M2方法在預報期1～16 d的R2從0.92下降到0.84，與M1方法相比，散點更接近于1∶1線附近，在PM估算值大于8 mm/d時，存在一定的低估問題。M3方法在預報期1～16 d的R2從0.93下降到0.85，與M1和M2方法相比，其散點更接近于1∶1線附近，且在預報期1 d提升最明顯，此后提升效果逐漸減弱。

從農業灌溉角度出發，一個灌溉周期連續多日的平均誤差比單日誤差更適宜為灌溉決策提供參考，因此，計算了1～8 d和1～16 d 3種方法的統計指標，并列出了未進行偏差校正(僅使用IDW方法)的統計指標作為對照，其結果如表3所示。從表3可以看出，IDW方法下只有站點2、3的RMSE在1 mm/d以內，說明該數據集整體誤差很大難以直接使用，需要偏差校正。對比以上4種方法統計指標可以看出，M1、M2、M3 3種方法均可以大幅度降低模型的預報誤差。與IDW方法相比，M1、M2和M3方法的RMSE在9個站點1～8 d分別下降26%～65%、24%～70%和31%～73%，在9～16 d分別下降20%～63%、20%～68%和24～70%。盡管M3方法相比M1方法在1～8 d和1～16 d分別有13%和15%的提升，但是兩者RMSE的絕對誤差在0.15 mm/d以內。從全年角度來看可能并沒有明顯的優勢。因此，有必要評估不同方法在不同季節，特別是灌溉季節的表現。

表3 不同方法下ET0在1～8 d和1～16 d累計誤差

2.3 ET0預報季節差異

3種方法在不同季節的RMSE如圖5所示。在春季，M1、M2和M3方法的1～16 d平均RMSE為0.89、0.90、0.80 mm/d。3種方法的RMSE隨預報期的延長呈增加趨勢，其中M1方法的RMSE從0.74 mm/d增加到1.03 mm/d，M2方法從0.86 mm/d緩慢增加到0.94 mm/d，M3方法的RMSE在前8 d從0.61 mm/d快速增加到0.83 mm/d，之后增速逐漸放緩，在16 d達到0.9 mm/d。對比M1和M2方法可以看出，在前8 d，M1方法表現優于M2方法，而在后8 d，M2方法優于M1方法。

在夏季，3種方法的表現與春季相似，1～16 d平均RMSE為1.21、1.18、1.04 mm/d。不同之處在于由于夏季ET0比春季更高，因此夏季的誤差更大。M1方法的RMSE在0.98～1.4 mm/d，而M2和M3方法可以把最大誤差控制在1.2 mm/d附近。此外，在前6 d，M1方法表現優于M2方法，而在6 d以后M2方法優于M1方法。秋季和冬季呈現與春季和夏季相似的趨勢。M3方法在預報初期優于其余2種方法，但隨著預報期增長，其優勢逐漸變小，在8 d以后，3種方法的RMSE相差約在0.1 mm/d。

2.4 氣象因子對誤差的影響

為了發現NWP系統的不足之處，進一步評估了單個氣象因子的預報誤差對ET0預報精度的影響。使用M3方法并分別使用1種預報數據與其余4種觀測數據組合代入PM公式，來評估氣象因子預報誤差對ET0誤差的影響，結果見圖6。從圖6a、6b可以看出，最高和最低氣溫的誤差對ET0誤差的影響在不同季節有所差異，由大到小為春季、夏季、秋季、冬季。此外，Tmax預報誤差對ET0造成的誤差比Tmin大。在夏季，Rs的誤差對ET0影響更大，ET0的RMSE在1～16 d介于0.55～0.83 mm/d(圖6c)，與前文相比，其誤差已經達到了總誤差的70%以上。而春季，由于Rs預報誤差導致ET0的最大RMSE僅為0.42 mm/d。U2預報誤差對ET0影響在春季和夏季相差不大，其ET0的最大RMSE約為0.4 mm/d，秋季則高于冬季(圖6d)。RH預報誤差對ET0影響較小，在夏季與春季相差不大，其ET0的最大RMSE約為0.2 mm/d(圖6e)。對比不同氣象因子造成的影響可以看出，Rs對ET0誤差貢獻最大，其后依次是U2、Tmax、RH和Tmin。從2.1節可以看出，各氣象因子中預報誤差最大的是U2，其后依次是RH、Rs、Tmax和Tmin。由于不同因子對ET0的貢獻度不同，造成了其自身預報誤差與對ET0的誤差貢獻有所不同。盡管U2的誤差最大，但由于其對ET0影響較小，因此對ET0誤差影響僅排在第2位，而Rs盡管預報精度排第3，但由于其對ET0貢獻最大，因此其對ET0誤差的貢獻也最大。這表明提高Rs的預報精度是更準確預報ET0的關鍵。

2.5 預報氣象因子與觀測氣象因子之間關系

為更好地認識數值天氣預報系統氣象因子內在的聯系，基于LightGBM模型對氣象因子進行偏差校正基礎上，通過算法對因子貢獻度評分(Gini值)繪制預報氣象因子之間的相關關系圖，其各輸入因子對輸出氣象因子的得分結果如圖7所示?？梢钥闯觯珿EFS_Tmax信息對于提升Tmax、Tmin、RH和U2的預報精度都有較大的幫助。GEFS_Tmax的Gini值為2.1，其中僅有約0.9的得分用于預報Tmax，約0.6的得分用于預報Tmin，甚至高于GEFS_Tmin對Tmin的預報得分(0.25)，這說明與Tmin變化規律最接近的是GEFS_Tmax，而非GEFS_Tmin。

GEFS_Rs除主要用于預報Rs外，對RH和U2的預報稍有幫助。此外，GEFS_RH對RH預報得分只有0.2，而GEFS_Tmax對RH預報得分為0.45，高于前者。這說明使用GEFS_Tmax來預報RH比GEFS的GEFS_RH預報更有效。綜合來看，Tmax對其他因子預報的貢獻最大，其后依次是Rs、Tmin、U2和RH。另一方面，多因子預報對太陽輻射預報的幫助最小。以上信息說明，預報因子與觀測氣象因子之間存在不匹配的問題，不匹配程度與因子有關，其中Rs的匹配度最高，RH匹配度最低。這也證明了本研究使用多因子偏差校正的合理性。

3 討論

目前，ET0預報主要基于公共天氣預報和數值天氣預報。公共天氣預報容易獲得、只需要簡單轉換就可以使用，跟經驗模型結合更能簡化計算流程[24]，但由于數據的完整度存在問題，公共天氣預報的預報期通常較短，多為7 d以內[5-7,25-26]。數值天氣預報具有更完整的數據集，且預報期更長，通常可以預報10 d以上，甚至1個月。但數值天氣預報通常前處理比較復雜，并且需要偏差校正，這給其廣泛使用帶來了一定的困難，在NWP基礎上進行后處理，構建預報數據集有助于簡化其復雜度。YANG等[7]比較了不同4種風速輸入方式在中國不同氣候區ET0預報的表現，其中在干旱區RMSE介于0.95～0.99 mm/d，MAE介于0.67～0.70 mm/d。本研究中，M3方法在9站點1～16 d平均RMSE介于0.53～0.79 mm/d，MAE介于0.35～0.53 mm/d，說明NWP結合后處理方法比公共天氣預報的精度更高。此外，本研究結果還證明了機器學習方法比統計方法精度更高，在相同輸入條件下，M2方法比M1方法RMSE低12.0%～15.7%，MAE低8.9%～15.8%。

對NWP數據進行后處理可以提高模型精度，MEDINA等[27]比較了3種基于統計的后處理方法(非其次高斯回歸、仿射核修整和貝葉斯模型平均)對ET0預報的提升效果，發現3種方法都能有效提高ET0預報精度。統計方法的不足之處在于通常只能單因素偏差校正，如使用輻射預報值跟輻射實測值建立關系。本研究發現，預報因子和對應的觀測因子之間存在不匹配的問題，如與Tmin、RH最匹配的預報因子是Tmax，同樣RH最相關的預報因子也是Tmax。也就是說名義上預報的某氣象因子由于預測誤差原因與實測的另一個氣象因子的關系更加緊密。這首先由于氣象因子之間本身就存在著互相影響的關系，如輻射是造成氣溫日變化的主要原因，同樣氣溫會影響空氣中容納的水汽含量，進而影響到相對濕度；其次，可能是由于NWP產品基于物理模型，需要對氣象要素之間復雜的非線性關系進行簡化；最后，NWP的邊界條件往往比較難以控制，這就造成了一定程度預報精度的下降[28]，而部分有價值的信息被隱藏在了其他預報因子之中，使用LightGBM這類基于提升樹的集合算法，可以將多維輸入共同作為輸入特征，發掘其中的非線性關系并剔除干擾特征。綜上，推薦在對NWP數據進行偏差校正時，應考慮使用更多的相關因子建模，來改善NWP模型的不足?？紤]到本研究中太陽輻射的預報精度還有待進一步提高，今后應使用數值天氣預報中更多預報因子，如降水量、云量、反照率、土壤熱通量等來提高太陽輻射的預報精度，進而進一步提升ET0預報能力。

4 結論

(1)預報因子與實測氣象因子之間存在不匹配問題，其中太陽輻射的匹配度較高，最低氣溫、相對濕度的匹配度較低，需要借助最高氣溫等信息進行校正，使用多因子共同預報氣象因子有助于氣象因子校正。

(2)校正后統計指標顯示，M3方法在太陽輻射、最高和最低氣溫、相對濕度和風速預報上均優于M1和M2方法。M1、M2和M3方法在9站點預報ET0平均RMSE分別介于0.66～0.93 mm/d、0.57～0.83 mm/d和0.53～0.79 mm/d，MAE分別介于0.44～0.61 mm/d、0.38～0.56 mm/d和0.35～0.53 mm/d，R2分別介于0.82～0.91、0.84～0.93和0.86～0.94。3種方法均在夏季誤差最大，1～16 d平均RMSE分別為1.21、1.18、1.04 mm/d。

(3)各預報因子中太陽輻射對ET0預報誤差影響最大，其后依次是風速、最高氣溫、相對濕度和最低氣溫。最高氣溫預報對其他因子預報精度的貢獻最大、對相對濕度預報精度的貢獻最小。