不同核函數高斯過程回歸算法與不同因子輸入情況下對長江流域蒸散發量應用研究

楊梓涵,崔崢錚,張鵬程

(河海大學 水文水資源學院,南京 210024)

0 引 言

參考作物騰發量(Reference evapotranspiration,ET0)是自然界水分循環與能量平衡的重要組成部分,與降水等因素共同影響著區域的干濕狀況,也是估算生態需水和農業灌溉用水的關鍵因子[1]。長江流域地處亞熱帶季風氣候區,年降水量大,流域內有豐富的自然資源,是我國重要的水利經濟帶。因此,對長江流域內部分城市的參考作物騰發量進行研究,對水資源合理開發利用與水循環有促進和指導意義。

目前,常用的參考作物騰發量估算公式主要有Penman-Monteith、Priestley-Taylor、Hargreaves-Samani和Thornthwaite 公式[2],這4個公式對長江流域都具有一定的適用性。由于每個公式計算因子不同,對不同地區蒸散發計算的結果精度也不同,許多學者針對不同地區對ET0估算公式進行了適用性分析并進行修正改進。如黃強等[2]針對珠江流域42個觀測站點1959-2011年的逐日氣象數據,以實測蒸散發數據為基準,探討了PM公式、PT公式、HS公式和TH公式在珠江流域蒸散發計算中的適用性,結果表明PM公式在不同季節和地區的計算結果都較為穩定,與實測值更接近,最適用于珠江流域的蒸散發計算。朱瀟梟等[3]針對海南省7個氣象站點2000-2014年的逐日氣象數據,以PM公式計算結果為標準,對Priestley-Taylor法,Irmak-Allen法和Hargreaves-Samani法的適用性進行了評價和對比分析,結果表明3種不同的ET0計算方法在海南省的相關性依次為:PT公式>IA公式>HS公式。

ET0的計算可以看作是一個依賴于大量氣象變量的復雜非線性回歸過程。因此,很難建立準確的經驗模型來代表所有復雜的過程[4]。近年來,有學者提出ET0估計的機器學習方法,因為它們可以為非線性和多變量函數提供簡單的解決方案[5]。FAN等[4]基于我國不同氣候區8個氣象站點1960-2010的逐日氣象數據,對SVM、ELM、RF、M5tree、GBDT、XGBoost這6種高斯過程回歸算法模擬ET0的精度進行了對比,結果表明GBDT和XGBoost在模擬ET0中有著顯著的優越性,但SVM和ELM是模擬ET0最穩定的兩個模型。趙文剛等[6]選取BP神經網絡、極限學習法和小波神經網絡3種較為常用的機器學習法,以廣東省典型代表站點1981-2010年的實測逐日氣象數據為研究對象,對比了3種高斯過程回歸算法的可靠性,結果表明ET0計算精度表現為:BP>極限學習機>小波神經網絡。此外,還有部分學者針對有限氣象數據條件下高斯過程回歸算法計算ET0的精度進行了研究[7-8]。

因此,本文以長江流域10個站點1970-2019共50年的逐日氣象數據為研究對象,以PM公式計算的ET0為參考值,對比不同因子輸入組合情況下,采用PT公式、HS公式、TH公式和3種不同核函數的GPR算法為計算方法,計算逐日參考作物騰發量的精度。本研究目的:①對比全氣象因子輸入情況下,3種ET0經驗公式與3種不同核函數的GPR算法計算參考作物騰發量的計算精度。②對比不同因子輸入情況下,3種不同核函數GPR算法的計算精度,并嘗試分析因子輸入改變情況下計算精度改變的原因。

1 材料與方法

1.1 數據資料

本研究所使用的氣象資料由國家氣象科學中心(http://www.nmic.cn/)提供。基于站點空間分布均勻且具有較強代表性為原則,收集了長江流域10個常規氣象站近50年(1970-2019)的逐日氣象數據,主要包括最低氣溫(Tmin,℃)、最高氣溫(Tmax,℃)、平均氣溫(Tmean,℃)、相對濕度(RH,%)、平均本站氣壓(P,kPa)、日照時數(N,h)、風速(U2,m/s)共8項指標。

1.2 參考作物騰發量(ET0)計算方法

1.2.1 FAO-56 Penman-Monteith模型

由于FAO-56 PM公式計算精度較高,能夠較為準確反映參考作物騰發量的真實情況,因此近年來一直被廣泛應用。具體計算公式如下:

(1)

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

1.2.2 Hargreaves-Samani模型

HS公式是一種基于溫度來估算ET0的方法,由于需要的氣象數據較少,常用于氣象設備缺乏地區的ET0預報。其表達式為:

(2)

式中:Ra為大氣輻射,MJ/(m2·d);C、E、T為Hargreaves公式的3個參數,建議值分別為0.0023、0.5和17.8。

1.2.3 Priestley-Taylor模型

Priestley-Taylor(P-T)公式法是在假定周圍濕潤的條件下提出來的,因此該方法一般適用于在濕潤地區計算參考作物騰發量。其公式如下:

(3)

式中:α系數主要考慮空氣動力因素影響,一般情況下取1.26;λ=2.501-0.002361T或直接取2.45。

1.2.4 Irmak-Allen模型

Irmak-Allen(I-A)見式(6):

ET0=0.489+0.289Rn+0.023Tmean

(4)

1.3 高斯過程回歸算法

1.3.1 高斯過程回歸(GPR)

高斯過程回歸(GPR)數學模型是依賴于非參數核的概率模型,它在機器學習編程領域具有重要意義[9]。

假設GPR模型中輸入向量(x)和目標(y)之間的運算關系為:

yi=f(xi)+ε

(5)

式中:f(xi)為任意函數;ε為噪聲(回歸殘差)。

本研究中,認定各選定的氣候指標為輸入向量,參考作物騰發量為輸出向量。

在函數f(x)中,f是因變量,x為自變量,并且任何未觀察到的數對(x*,f*)見式(6)。

(6)

式中:K(X,X) 為訓練數據中所有點之間的n×n協方差矩陣;k(X,x*) 為點x*和訓練數據之間的n×1協方差向量。

p(f*|x*,X,f)=N(k(x*,X)K(X,X)-1f,k(x*,x*)-k(x*,X)K(X,X)-1k(X,x*))

(7)

(6)式通過最大化f*條件下的聯合概率來確定f*。為了利用有噪聲的數據,該模型必須伴隨有測量誤差。因此,式(5)常常被改為如下形式:

(8)

條件似然函數變化為:

方差變化為:

式中:σ2為觀測誤差的方差;I為單位矩陣。

本文研究中運用的3種核函數如下:

1)二次有理核函數。公式為:

(9)

2) 平方指數協方差函數(SE)。公式為:

(10)

3)Matern協方差函數(Matern 5/2)。公式為:

(11)

式中:Kv為修正貝塞爾函數。

1.3.2 評價方法

ET0精度采用平均絕對誤差(MAE)、均方根誤差(RMSE)和決定系數R2這3個指標進行評價[10],各指標的計算公式為:

(12)

(13)

(14)

通常情況下,MAE和RMSE越小,表明模型偏差越小,精度越高;R2越接近1,表明模型的擬合程度越高。

1.4 統計與分析

使用Microsoft Excel 2011對各站點氣象數據進行整理,并采用FAO56-PM公式、HS公式、Priestley-Taylor公式和Irmak-Allen公式對ET0進行計算,采用Matlab 2021b對機器學習進行訓練與預測,并采用Origin 2021對3種經典算法和3種高斯過程回歸算法的預測精度進行圖表繪制。

2 結果與分析

2.1 基于灰色關聯分析法分析不同氣象因子對逐日蒸散發的影響程度

灰色關聯分析法是一種根據因素之間發展趨勢的相似或者相異程度,即“灰色關聯度”,作為用來衡量因素間的關聯程度的一種方法。

對ET0影響程度排序為:日最高氣溫>平均氣溫>日最低氣溫>平均本站氣壓>平均風速>相對濕度>日照時數,將關聯度數據進行可視化,見圖1。

圖1 關聯度可視化

利用灰色關聯度分析可得,日最高溫度對ET0影響較大,日照時數對ET0影響較小,其余子因素對ET0的影響適中。

2.2 基于同一參考公式全因子輸入下的3種經典算法和3種不同核函數高斯過程回歸算法預測精度對比

為了比較Hargreaves-Samani、Priestley-Taylor、Irmak-Allen等3種經典算法和二次有理、平方指數、Matern 5/2等3種不同核函數高斯過程回歸算法在全因子輸入下預測ET0方面的精度情況,本研究首先通過PM公式、HS公式、Priestley-Taylor公式和Irmak-Allen公式,計算出各站點逐日參考作物騰發量ET0,PM、ET0,HS、ET0,PT、ET0,IA,再分別以前40年(1970-2009年)逐日氣象資料(Tmin、Tmax、Tmean、RH、P、N和U2)和ET0值(ET0,PM)為機器學習輸入變量和輸出訓練集,并最終通過輸入后10年(2010-2019)逐日氣象因子,模擬該時間序列上的ET0值,最后與對應的ET0,PM值進行對比分析,其具體結果見圖2。

圖2 PM參考公式下3種經典算法和3種不同核函數高斯過程回歸算法評價指標對比

由圖2可以看出,當以PM公式計算的ET0結果為參考時,Hargreaves-Samani公式計算結果的R2、MAE和RMSE分別維持在0.396～0.871、4.235～5.378和4.900～6.279范圍內;Priestley-Taylor公式計算結果的R2、MAE和RMSE分別維持在0.422～0.951、0.017～0.396和0.348～2.805范圍內;Irmak-Allen公式計算結果的R2、MAE和RMSE分別維持在0.424～0.945、0.284～0.952和0.472～1.419范圍內。二次有理核函數模擬結果的決定系數R2的變化范圍為0.970～0.988,MAE的變化范圍為0.003～0.105,RMSE的變化范圍為0.172～0.232;平方指數核函數模擬結果的決定系數R2的變化范圍為0.957～0.988,MAE的變化范圍為0.002～0.105,RMSE的變化范圍為0.172～0.269;Matern 5/2核函數模擬結果的決定系數R2的變化范圍為0.970～0.988,MAE的變化范圍為0.001～0.103,RMSE的變化范圍為0.171～0.233。

從整體上分析可以得到,在以PM為參考公式情況下,6種算法的模擬ET0表現效果均為:Matern 5/2>二次有理>平方指數>PT>IA>HS,總體上3種不同核函數高斯過程回歸算法的預測精度遠優于3種經典算法。

2.3 基于同一參考公式不同因子輸入下的3種不同核函數高斯過程回歸算法預測精度對比

為了比較二次有理、平方指數、Matern 5/2等3種不同核函數高斯過程回歸算法在不同因子輸入下預測ET0方面的精度情況,本研究首先通過PM公式計算出各站點逐日參考作物騰發量ET0,PM,再根據2.1一節灰色關聯分析結果,確定四因子、五因子、六因子3種因子輸入組合。四因子輸入時,分別以前40年(1970-2009年)逐日氣象資料(Tmin、Tmax、Tmean和P)和ET0值(ET0,PS)作為機器學習輸入變量和輸出訓練集;五因子輸入時,分別以前40年(1970-2009年)逐日氣象資料(Tmin、Tmax、Tmean、P和N)和ET0值(ET0,PM)作為機器學習輸入變量和輸出訓練集;六因子輸入時,分別以前40年(1970-2009年)逐日氣象資料(Tmin、Tmax、Tmean、P、N和U2)和ET0值(ET0,PS)作為機器學習輸入變量和輸出訓練集,并最終通過輸入后10年(2010-2019)逐日氣象因子,模擬該時間序列上的ET0值,最后與對應的ET0,PM值進行對比分析。具體結果見圖3。

由圖3可以看出,當以PM公式計算的ET0結果為參考時,二次有理核函數在四因子輸入下模擬結果的R2、MAE和RMSE分別維持在0.569～0.932、0.025～0.182和0.407～0.869范圍內;二次有理核函數在五因子輸入下模擬結果的R2、MAE和RMSE分別維持在0.804～0.963、0.011～0.117和0.303～0.560范圍內;二次有理核函數在六因子輸入下模擬結果R2、MAE和RMSE分別維持在0.908～0.977、0.027～0.202和0.055～0.454范圍內。

當以PM公式計算的ET0結果為參考時,平方指數核函數在四因子輸入下模擬結果的決定系數R2的變化范圍為0.568～0.932,MAE的變化范圍為0.026～0.181,RMSE的變化范圍為0.407～0.870;平方指數核函數在五因子輸入下模擬結果的決定系數R2的變化范圍為0.801～0.963,MAE的變化范圍為0.011～0.120,RMSE的變化范圍為0.303～0.565;平方指數核函數在六因子輸入下模擬結果的決定系數R2的變化范圍為0.910～0.976,MAE的變化范圍為0.028～0.202,RMSE的變化范圍為0.252～0.453。

當以PM公式計算的ET0結果為參考時,Matern 5/2核函數在四因子輸入下模擬結果的決定系數R2的變化范圍為0.569～0.932,MAE的變化范圍為0.028～0.182,RMSE的變化范圍為0.406～0.868;Matern 5/2核函數在五因子輸入下模擬結果的決定系數R2的變化范圍為0.803～0.963,MAE的變化范圍為0.011～0.118,RMSE的變化范圍為0.303～0.561;Matern 5/2核函數在六因子輸入下模擬結果的決定系數R2的變化范圍為0.909～0.977,MAE的變化范圍為0.026～0.201,RMSE的變化范圍為0.251～0.453。

從整體上分析可以得到,在以PM為參考公式的情況下,3種不同核函數高斯過程回歸算法的模擬ET0表現效果均為:六因子>五因子>四因子。

3 討論與分析

3.1 經典算法和機器學習算法適用性對比

不同經典算法對ET0的計算結果有顯著不同,同樣不同機器學習算法對ET0預測結果有顯著影響。如任傳棟等[10]以山東省為研究區域,選取深度神經網絡(DNN)、時間卷積神經網絡(TCN)和長短期記憶神經網絡(LSTM) 3種深度學習模型,極限學習機模型(ELM)、廣義回歸神經網絡模型(GRNN)和隨機森林模型(RF) 3種傳統機器學習模型,Hargreaves-Samani模型(HS)、Droogres-Allen模型(DA)、Priestley-Tayor模型(PT)、Marrink模型(MK)、WMO模型(WMO)、Trabert模型(TRA) 6種經驗模型,以均方根誤差(RMSE)、決定系數(R2)、平均絕對誤差(MAE)和效率系數(Ens)為精度評價體系,找出了適用于山東省ET0估算的最優模型,結果表明,相同氣象參數輸入條件下,機器學習模型精度普遍優于經驗模型。李可利[11]以陜西省為研究區域,利用陜西省30個氣象站點(1960-2016年)逐日氣象資料,采用FAO P-M公式計算出標準ET0,對ET0進行時空分布特征分析、敏感性分析,并對ET0進行周期分析、預測,對ANFIS、RF、SVM 3種機器學習算法在陜西省ET0計算的適用性與可移植性進行分析,結果表明,在使用相同氣象因子計算ET0時,SVM模型要優于H-S、Iramk、Makkink、P-T等4種傳統ET0計算方法,SVM可以應用于陜西省缺少氣象資料地區的ET0計算,代替其他ET0簡化方法,提高計算精度。上述研究成果與本研究的結果保持一致。高斯過程回歸的本質是通過一個映射把自變量從低維空間映射到高維空間的過程,由于選擇核函數與懲罰參數的不同,該模型的預測精度也不同,所以該模型存在著一定的優化空間,有待進一步深入研究。

3.2 結 論

本文以長江流域10個氣象站點1970-2019年的日氣象數據為基礎,以PM公式計算的ET0結果為參考值,研究了3種經典算法在計算ET0方面和3種不同核函數高斯過程回歸算法在預測ET0方面的表現。結論如下:

1)對于3種不同核函數高斯過程回歸算法和3種經典算法來說,總體上,3種不同核函數高斯過程回歸算法和3種經典算法的預測精度大小關系表現為:Matern 5/2>二次有理>平方指數>PT>IA>HS。

2)利用灰色關聯度分析得到,日最高溫度對參考作物騰發量影響較大;日照時數對參考作物騰發量影響較小;其余氣象因子對參考作物騰發量的影響適中。

3)對于不同因子輸入下同種核函數的高斯過程回歸算法,在以PM為參考公式的情況下,3種不同核函數高斯過程回歸算法的模擬ET0表現效果均為:六因子>五因子>四因子。