不同核函數(shù)高斯過程回歸算法與不同因子輸入情況下對長江流域蒸散發(fā)量應(yīng)用研究

楊梓涵,崔崢錚,張鵬程

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

0 引 言

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

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

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

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

1 材料與方法

1.1 數(shù)據(jù)資料

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

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

1.2.1 FAO-56 Penman-Monteith模型

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

(1)

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

1.2.2 Hargreaves-Samani模型

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

(2)

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

1.2.3 Priestley-Taylor模型

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

(3)

式中:α系數(shù)主要考慮空氣動力因素影響,一般情況下取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)數(shù)學(xué)模型是依賴于非參數(shù)核的概率模型,它在機器學(xué)習(xí)編程領(lǐng)域具有重要意義[9]。

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

yi=f(xi)+ε

(5)

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

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

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

(6)

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

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*條件下的聯(lián)合概率來確定f*。為了利用有噪聲的數(shù)據(jù),該模型必須伴隨有測量誤差。因此,式(5)常常被改為如下形式:

(8)

條件似然函數(shù)變化為:

方差變化為:

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

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

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

(9)

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

(10)

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

(11)

式中:Kv為修正貝塞爾函數(shù)。

1.3.2 評價方法

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

(12)

(13)

(14)

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

1.4 統(tǒng)計與分析

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

2 結(jié)果與分析

2.1 基于灰色關(guān)聯(lián)分析法分析不同氣象因子對逐日蒸散發(fā)的影響程度

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

對ET0影響程度排序為:日最高氣溫>平均氣溫>日最低氣溫>平均本站氣壓>平均風(fēng)速>相對濕度>日照時數(shù),將關(guān)聯(lián)度數(shù)據(jù)進行可視化,見圖1。

圖1 關(guān)聯(lián)度可視化

利用灰色關(guān)聯(lián)度分析可得,日最高溫度對ET0影響較大,日照時數(shù)對ET0影響較小,其余子因素對ET0的影響適中。

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

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

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

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

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

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

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

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

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

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

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

3 討論與分析

3.1 經(jīng)典算法和機器學(xué)習(xí)算法適用性對比

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

3.2 結(jié) 論

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

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

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

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