基于不同LSTM模型和Hargreaves模型估算魯中地區(qū)參考作物蒸散量

馬 釗， 任傳棟， 劉 靜， 王志真

(1.山東省水利勘測設(shè)計院有限公司，山東濟南250013；2.山東省農(nóng)業(yè)交流合作中心，山東濟南250013)

隨著全球變暖現(xiàn)象的日益加劇，區(qū)域干旱發(fā)生頻率逐漸提高，這嚴重影響居民生活及人身安全[1-2]。同時，水資源供需矛盾的逐漸加劇，給區(qū)域農(nóng)業(yè)發(fā)展帶來了嚴峻挑戰(zhàn)，嚴重影響了區(qū)域糧食產(chǎn)量及農(nóng)業(yè)產(chǎn)值[3-4]。參考作物蒸散量(ET0)是影響區(qū)域水資源平衡和能量傳遞的重要參數(shù)之一，其值的準確估算對區(qū)域經(jīng)濟發(fā)展有著十分重要的意義[5]。Penman-Monteith(PM)模型是計算ET0的標準模型，但利用該模型計算ET0時需用到溫度、日照時數(shù)、相對濕度、風(fēng)速等多個氣象數(shù)據(jù)，使得該模型在氣象資料缺乏地區(qū)的應(yīng)用受限[6-7]。因此，在保證ET0估算精度的前提下，找尋區(qū)域ET0簡化計算模型成為國內(nèi)外的研究熱點。

Hargreaves等[8]基于溫度資料及大氣頂層輻射構(gòu)建了Hargreaves(HS)模型，該模型在僅需要溫度資料的條件下就能以較高精度估算區(qū)域ET0，但該模型中的經(jīng)驗參數(shù)需針對不同區(qū)域進行屬地化處理，因此許多學(xué)者對HS模型進行了參數(shù)校正。胡慶芳等[9]基于全局優(yōu)化的洗牌算法對中國不同區(qū)域HS模型參數(shù)進行了校正，得到了全國不同分區(qū)的HS模型參數(shù)參考值。李晨等[10]基于貝葉斯原理在川中丘陵區(qū)對HS模型進行了參數(shù)校正，指出改進后的HS模型精度明顯提高。夏興生等[11]基于回歸修正原理在中國農(nóng)業(yè)區(qū)對HS模型進行了修正，同樣取得了較高的精度。

近年來，機器學(xué)習(xí)算法憑借計算的簡便性及高精度被廣泛應(yīng)用于區(qū)域ET0估算中。Tabari等[12]基于支持向量機模型估算了伊朗半干旱區(qū)的ET0，指出模型精度高于HS、Priestley-Taylor和Blaney-Criddle等經(jīng)驗?zāi)Ｐ?。武劍飛等[13]在氣象資料受限條件下基于NARX模型對攀枝花站點ET0進行了估算，指出了該模型精度普遍優(yōu)于HS模型。張薇等[14]基于梯度提升決策樹、隨機森林、回歸樹3種機器學(xué)習(xí)模型估算江蘇省ET0，并將計算結(jié)果與3種經(jīng)驗?zāi)Ｐ蛯Ρ龋赋鎏荻忍嵘龥Q策樹精度最高，可作為江蘇省簡化ET0計算模型使用。不同研究結(jié)果均表明，相同氣象參數(shù)輸入條件下，機器學(xué)習(xí)模型的精度普遍優(yōu)于經(jīng)驗?zāi)Ｐ汀?/p>

長短期記憶神經(jīng)網(wǎng)絡(luò)(LSTM)模型將時間概念引入模型結(jié)構(gòu)中，已被廣泛應(yīng)用于數(shù)據(jù)預(yù)測和模擬中[15-16]。但傳統(tǒng)的LSTM模型精度易受到參數(shù)取值的影響，因此，需采用智能算法對其進行優(yōu)化。本研究擬基于貝葉斯理論(BA)、粒子群算法(PSO)和遺傳算法(GA)3種算法優(yōu)化LSTM模型，選取魯中地區(qū)6個氣象站點，在僅輸入溫度數(shù)據(jù)的條件下，估算區(qū)域ET0，將3種優(yōu)化模型的計算結(jié)果與傳統(tǒng)LSTM模型、傳統(tǒng)HS模型和5種HS優(yōu)化模型進行對比，在氣象數(shù)據(jù)受限條件下，得出區(qū)域ET0簡化估算最優(yōu)模型，以期為氣象數(shù)據(jù)缺乏地區(qū)灌溉制度的制訂提供科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)域概況及數(shù)據(jù)來源

魯中地區(qū)地處山東省中部，屬于暖溫帶大陸性季風(fēng)氣候區(qū)，區(qū)域四季分明，主要地形包括山區(qū)及平原，是山東省重要的糧食生產(chǎn)基地[17]。為得出魯中地區(qū)ET0估算的最優(yōu)簡化模型，本研究選擇魯中地區(qū)濟南、泰山、濰坊、淄博、淄川、沂源6個氣象站點1961-2019年的逐日氣象數(shù)據(jù)，數(shù)據(jù)包括日最高氣溫(Tmax)、日最低氣溫(Tmin)、日照時數(shù)(n)、相對濕度(RH)和離地面10 m高處風(fēng)速(U10)，氣象數(shù)據(jù)均來自國家氣象中心，數(shù)據(jù)質(zhì)量控制良好。研究區(qū)域地形及站點分布情況見圖1。

1.2 研究方法

1.2.1 Penman-Monteith模型 1998年FAO-56分冊中規(guī)定Penman-Monteith(PM)模型為計算ET0的標準模型[18]，具體公式為：

(1)

式中，Rn為作物表面的凈輻射[MJ/(m2·d)]；T為離地面2 m高處的日均氣溫(℃)；G為土壤熱量通量密度[MJ/(m2·d)]；U2為離地面2 m高處風(fēng)速(m/s)；es為飽和水汽壓差(kPa)；ea為實際水汽壓差(kPa)；△為蒸汽壓曲線的斜率(kPa/℃)；γ為干濕計常數(shù)(kPa/℃)。

圖1 研究區(qū)域概況及站點分布Fig.1 Research area overview and site distribution

1.2.2 長短期記憶神經(jīng)網(wǎng)絡(luò)模型 長短期記憶神經(jīng)網(wǎng)絡(luò)(LSTM)模型引入時間序列概念，使模型識別訓(xùn)練數(shù)據(jù)集的時間規(guī)律，保證了該模型的精度。模型由輸入門、遺忘門、更新門和輸出門4個部分組成，具體步驟可見文獻[19]。

1.2.3 貝葉斯優(yōu)化LSTM模型 貝葉斯優(yōu)化基于高斯回歸原理，對模型進行參數(shù)尋優(yōu)。假設(shè)待優(yōu)化的參數(shù)組合為X={x1,x2,…,xn}，模型允許的計算輸出結(jié)果誤差函數(shù)為f(x)，在預(yù)測過程中判定輸出參數(shù)能否符合條件，其誤差函數(shù)符合二項分布：

f(x)=GP[μ,k(x,x)]

(2)

式中，GP為高斯分布；μ為序列均值；k(x,x)為序列協(xié)方差。

通過貝葉斯優(yōu)化的第t-1組數(shù)據(jù)為[xt-1,f(xt-1)]，則第t組數(shù)據(jù)的預(yù)測值f(xt)為符合高斯分布的一個樣本，公式為：

(3)

(4)

k=[(xt,x1),(xt,x2),…,(xt,xt-1)]

(5)

式中，x1、xt-1、xt為樣本參數(shù)；μ為序列均值；ft、f1:t-1為樣本計算值。

聯(lián)立式(3)至式(5)即可解出第t個觀測值。在LSTM模型訓(xùn)練過程中，通過貝葉斯理論進行迭代計算，得出滿足要求的預(yù)測值。以原始數(shù)據(jù)集為起點，通過本研究提出的貝葉斯目標函數(shù)及LSTM模型，即可得出優(yōu)化后的LSTM模型(BA-LSTM)。

1.2.4 粒子群算法優(yōu)化LSTM模型 粒子群算法(PSO)基于群鳥喂養(yǎng)的原理，將每個訓(xùn)練數(shù)據(jù)視為一個粒子，全部數(shù)據(jù)集即為粒子群，應(yīng)用每個粒子本身的運動屬性完成算法參數(shù)尋優(yōu)過程[20]。在參數(shù)迭代過程中，算法通過計算數(shù)據(jù)個體的極值從而更新粒子運動屬性，直到參數(shù)滿足誤差要求為止。將PSO算法用于優(yōu)化LSTM模型，構(gòu)建PSO-LSTM模型，可快速選取最優(yōu)參數(shù)，提高模型的精度。具體步驟可參照文獻[21]。

1.2.5 遺傳算法優(yōu)化LSTM模型 遺傳算法(GA)通過交叉、變異處理計算訓(xùn)練數(shù)據(jù)的個體適應(yīng)度，通過比較個體適應(yīng)度的值能否滿足要求，判定是否輸出最終結(jié)果。該算法基于生物界適者生存的原理，保留最優(yōu)個體，從而提高模型參數(shù)的迭代速率。構(gòu)建的GA-LSTM模型可在傳統(tǒng)LSTM模型的基礎(chǔ)上提高模型的運行速率和精度，模型具體步驟可見參考文獻[22]。

1.2.6 Hargreaves模型 Hargreaves(HS)模型僅基于溫度資料便能以較高精度估算區(qū)域ET0[23]，具體公式為：

(6)

式中，ET0表示參考作物蒸散量(mm/d)；Tmax、Tmin、Tmean分別表示日最高溫度(℃)、日最低溫度(℃)、日平均溫度(℃)；Ra表示大氣頂層輻射[MJ/(m2·d)]；λ表示汽化潛熱常數(shù)。

1.2.7 Droogres改進的Hargreaves模型 Droogres等[24]基于氣象數(shù)據(jù)觀測的隨機性，提出了2種Hargreaves改進模型，分別為HSM1、HSM2模型，具體公式分別為：

ET0=0.408 0×0.003 0×[(Tmax-Tmin)0.400 0](Tmean+20.000 0)Ra

(7)

ET0=0.408 0×0.002 5×[(Tmax-Tmin)0.500 0](Tmean+16.800 0)Ra

(8)

式中，ET0表示參考作物蒸散量(mm/d)；Tmax、Tmin、Tmean分別表示日最高溫度(℃)、日最低溫度(℃)、日平均溫度(℃)；Ra表示大氣頂層輻射[MJ/(m2·d)]。

1.2.8 胡慶芳等改進的Hargreaves模型 胡慶芳等[25]基于洗牌算法對全國不同分區(qū)HS模型參數(shù)進行了校正，改進的模型為HSM3模型，基本公式為：

ET0=0.001 0×Ra[(Tmax-Tmin)0.660 0](Tmean+34.500 0)

(9)

式中，ET0表示參考作物蒸散量(mm/d)；Tmax、Tmin、Tmean分別表示日最高溫度(℃)、日最低溫度(℃)、日平均溫度(℃)；Ra表示大氣頂層輻射[MJ/(m2·d)]。

1.2.9 夏興生等改進的Hargreaves模型 夏興生等[11]對中國農(nóng)業(yè)區(qū)不同分區(qū)HS模型進行了回歸修正，改進的模型為HSM4模型，基本公式為：

ET0=a+bET0-HS

(10)

式中，ET0表示參考作物蒸散量(mm/d)；ET0-HS表示HS模型的參考作物蒸散量計算結(jié)果(mm/d)；a、b表示回歸經(jīng)驗參數(shù)。

1.2.10 基于貝葉斯原理改進的Hargreaves模型 李晨等[10]基于貝葉斯原理對川中丘陵區(qū)HS模型參數(shù)進行了修正，改進的模型為HSM5模型，基本公式為：

ET0=C×Ra[(Tmax-Tmin)m](Tmean+a)

(11)

式中，ET0表示參考作物蒸散量(mm/d)；Tmax、Tmin、Tmean分別表示日最高溫度(℃)、日最低溫度(℃)、日平均溫度(℃)；Ra表示大氣頂層輻射[MJ/(m2·d)]；C、m、a為修正后的經(jīng)驗參數(shù)。本研究同樣基于貝葉斯原理修正了魯中地區(qū)6個站點HS模型的經(jīng)驗參數(shù)，魯中地區(qū)不同站點修正后的參數(shù)取值見表1。

表1 不同站點HSM5模型參數(shù)取值

1.2.11 模型訓(xùn)練與驗證 僅使用溫度數(shù)據(jù)，以最高溫度(Tmax)和最低溫度(Tmin)為模型輸入數(shù)據(jù)，以1961-2000年的數(shù)據(jù)訓(xùn)練模型，2001-2019年的數(shù)據(jù)驗證模型，不同模型及算法的參數(shù)取值見表2。以均方根差(RMSE)、決定系數(shù)(R2)、平均絕對誤差(MAE)、效率系數(shù)(Ens)和綜合性指標指數(shù)(GPI)評價不同模型精度，計算公式分別為:

(12)

(13)

(14)

(15)

(16)

2 結(jié)果與分析

2.1 不同模型估算的ET0日值擬合效果分析

將魯中地區(qū)不同站點不同模型估算的ET0日值與PM模型ET0日值進行對比，表3顯示，不同模型精度不同，在相同氣象參數(shù)輸入的條件下，LSTM模型精度普遍優(yōu)于HS模型及HS改進模型。BA-LSTM模型估算的ET0日值在不同站點與PM模型的ET0日值擬合效果最好，其在不同站點的擬合方程斜率最接近標準值1，決定系數(shù)R2均在0.850以上。PSO-LSTM模型和GA-LSTM模型估算的ET0日值與PM模型的ET0日值擬合效果次之，在不同站點的擬合方程斜率分別為0.775～0.979和0.746～0.978，R2分別為0.795～0.919和0.754～0.915。在經(jīng)驗?zāi)Ｐ椭校琀SM5模型估算的ET0日值擬合效果最優(yōu)，其在不同站點的擬合方程斜率為0.710～0.901，決定系數(shù)R2在6種經(jīng)驗?zāi)Ｐ椭休^高。

表2 不同算法參數(shù)取值情況

表3 不同模型估算的不同站點ET0日值與PM模型ET0日值擬合效果對比

2.2 不同模型估算的ET0月值擬合效果分析

不同模型不同站點ET0月值與PM模型ET0月值的擬合結(jié)果(圖2)表明，不同模型估算的ET0月值有所差異，但在不同月份的變化趨勢基本一致，均呈拋物線趨勢變化。4種機器學(xué)習(xí)模型估算的ET0月值與PM模型的ET0月值擬合效果優(yōu)于其余模型，其中BA-LSTM模型估算的ET0月值與PM模型的擬合效果最優(yōu)。在經(jīng)驗?zāi)Ｐ椭?，HSM5模型估算的ET0月值與PM模型的擬合效果最優(yōu)，HS模型計算結(jié)果精度較低。HS模型和5種HS改進模型基本均高估了ET0月值，HS模型、HSM3模型和HSM4模型估算的ET0月值與PM模型的ET0月值擬合效果較差。

2.3 不同模型估算的ET0日值精度指標對比

為進一步比較不同模型的計算精度，本研究比較了不同模型估算的ET0日值，結(jié)果(圖3)顯示，4種LSTM模型估算的ET0日值比6種HS模型估算結(jié)果精度更高，綜合考慮，在所有模型中BA-LSTM模型估算的ET0日值精度最高，在6種經(jīng)驗?zāi)Ｐ椭蠬SM5模型估算的ET0日值精度最高。在ET0日值計算結(jié)果對比中，PM模型估算的ET0日值中位線為2.43 mm/d，BA-LSTM模型估算的ET0日值中位線為2.45 mm/d，其余LSTM模型估算的ET0日值中位線為2.46～2.48 mm/d；6種經(jīng)驗?zāi)Ｐ椭?，?yōu)化后的HS模型估算的ET0日值精度高于傳統(tǒng)HS模型，HSM5模型估算的ET0日值中位線與PM模型較接近，為2.65 mm/d，其余HS模型估算的ET0日值中位線為3.15～6.56 mm/d。此外，從5項精度指標綜合考慮，BA-LSTM模型估算的ET0日值精度最高，其均方根差(RMSE)、平均絕對誤差(MAE)、決定系數(shù)(R2)、效率系數(shù)(Ens)的中位線分別為0.378 mm/d、0.276 mm/d、0.904和0.902，其GPI的中位線為1.837，在所有模型中排名第一。PSO-LSTM模型和GA-LSTM模型估算ET0日值的精度次之，均優(yōu)于LSTM模型，GPI中位線分別為1.724和1.566。4種LSTM模型估算ET0日值的精度高于6種HS模型，在經(jīng)驗?zāi)Ｐ椭?，HSM5模型估算的ET0日值精度最高，其RMSE、MAE、R2、Ens的中位線分別為1.122 mm/d、0.838 mm/d、0.716和0.609，GPI中位線為0.166，在所有模型中排名第5位，其次為HSM1、HSM2、HSM3、HSM4模型，傳統(tǒng)HS模型估算ET0日值的精度最低。因此，在僅輸入溫度這一氣象參數(shù)條件下，BA-LSTM模型可用于估算魯中地區(qū)ET0。

ET0：參考作物蒸散量；PM：Penman-Monteith模型。LSTM、PSO-LSTM、BA-LSTM、GA-LSTM、HS、HSM1、HSM2、HSM3、HSM4、HSM5見表3注。圖2 不同模型估算的ET0月值與PM模型ET0月值的擬合對比Fig.2 Comparison of fitting effects between monthly ET0 estimated by different models and monthly ET0 estimated by PM model

a:PM;b:PSO-LSTM;c:BA-LSTM;d:GA-LSTM;e:LSTM;f:HS;g:HSM1;h:HSM2;i:HSM3;j:HSM4;k:HSM5。PM：Penman-Monteith模型。LSTM、PSO-LSTM、BA-LSTM、GA-LSTM、HS、HSM1、HSM2、HSM3、HSM4、HSM5見表3注。ET0：參考作物蒸散量；RMSE：均方根差；MAE：平均絕對誤差；R2：決定系數(shù)；Ens：效率系數(shù)；GPI：綜合性指標指數(shù)。圖3 不同模型ET0日值精度箱線圖Fig.3 Boxplots of daily ET0 accuracy of different models

2.4 不同模型估算的ET0日值相對誤差對比

圖4顯示，不同模型計算魯中地區(qū)ET0日值的相對誤差(RE)不同，但分布趨勢基本一致，均呈現(xiàn)西部相對誤差較低，中部相對誤差較高的趨勢。其中，BA-LSTM模型在整個魯中地區(qū)的RE均較低，RE最高值出現(xiàn)在沂源站附近，僅為1.75%。PSO-LSTM模型和GA-LSTM模型精度次之，RE最高值分別為2.01%和2.42%。在經(jīng)驗?zāi)Ｐ椭校琀SM5模型的RE較低，為2.08%～16.96%，其余改進的HS模型的精度較低，但均優(yōu)于傳統(tǒng)HS模型。HSM5模型可在保證計算簡便的前提下，明顯提高模型計算精度，機器學(xué)習(xí)模型的精度普遍優(yōu)于HS模型，BA-LSTM模型在整個區(qū)域均可保證較高的精度。

RE：相對誤差。LSTM、PSO-LSTM、BA-LSTM、GA-LSTM、HS、HSM1、HSM2、HSM3、HSM4、HSM5見表3注。圖4 不同模型估算的ET0日值相對誤差空間分布對比Fig.4 Comparison of spatial distribution of relative errors of daily ET0 estimated by different models

3 討論

溫度資料是氣象數(shù)據(jù)中最容易獲得的，并且其觀測精度較高，HS模型僅需溫度資料即可估算區(qū)域ET0。本研究發(fā)現(xiàn)，未優(yōu)化的HS模型在魯中地區(qū)的計算精度較低，這說明HS模型在應(yīng)用時需要對參數(shù)進行屬地化處理。HSM5模型基于貝葉斯原理，可將新計算得出的數(shù)據(jù)作為后驗分布，將后驗分布作為新一輪計算的先驗分布，以此保證模型計算的精確性。馮禹等[26]在四川盆地比較了不同HS改進模型的計算精度，指出基于貝葉斯原理優(yōu)化的HS模型計算精度最高；趙凡萱等[27]在寧夏地區(qū)比較了貝葉斯優(yōu)化的HS模型與其他修正模型的計算精度，得出的結(jié)論與本研究結(jié)果基本一致。

LSTM模型作為深度學(xué)習(xí)模型的一種[28]，基于傳統(tǒng)機器學(xué)習(xí)模型，在模型訓(xùn)練時引入時間序列，可以更好地反映訓(xùn)練數(shù)據(jù)在時間尺度上的變化規(guī)律，擺脫傳統(tǒng)機器學(xué)習(xí)模型參數(shù)隨機性選取的缺點，降低了模型訓(xùn)練的隨機性。從本研究的箱線圖中可以看出，4種LSTM模型精度指標的標準差較小，表明LSTM模型無需進行參數(shù)屬地化處理即可應(yīng)用于區(qū)域ET0的估算。李莉等[29]在研究番茄蒸散量估算時比較了LSTM模型與Elman模型、RNN模型的精度，指出LSTM模型具有較高的預(yù)測精度和泛化性能。謝家興等[30]比較了LSTM模型、ELM模型、GRNN模型的精度，也得出了相同的結(jié)論。

本研究結(jié)果表明，經(jīng)算法優(yōu)化的LSTM模型可進一步提高精度，其中，BA-LSTM模型精度最優(yōu)。PSO算法和GA算法已被證明可用于區(qū)域ET0估算的機器學(xué)習(xí)算法優(yōu)化中，算法的應(yīng)用可普遍提高機器學(xué)習(xí)模型精度和迭代速率。貝葉斯理論可在傳統(tǒng)模型參數(shù)基礎(chǔ)上，自動調(diào)節(jié)模型參數(shù)，優(yōu)化模型結(jié)構(gòu)，提高模型的泛化能力。綜上所述，在僅有溫度這一氣象參數(shù)條件下，可推薦使用BA-LSTM模型估算魯中地區(qū)ET0。

4 結(jié)論

本研究在僅輸入溫度數(shù)據(jù)的條件下，比較了BA-LSTM、PSO-LSTM、GA-LSTM、LSTM共4種機器學(xué)習(xí)模型和6種HS模型的精度，找出適用于魯中地區(qū)ET0的簡化估算模型，得出以下結(jié)論：

(1)在ET0日值估算中，4種LSTM模型的估算效果明顯優(yōu)于6種HS模型，其中，BA-LSTM模型對ET0日值的估算效果最優(yōu)，其在不同站點的擬合方程斜率為0.812～0.982，決定系數(shù)均在0.850以上。

(2)在ET0月值估算中，不同模型估算結(jié)果的變化趨勢基本一致，以BA-LSTM模型估算的結(jié)果與PM模型的擬合效果最優(yōu)。

(3)從誤差指標中可以看出，BA-LSTM模型估算精度最高(RMSE、MAE、R2、Ens、GPI的中位線分別為0.378 mm/d、0.276 mm/d、0.904、0.902和1.837)，經(jīng)驗?zāi)Ｐ椭蠬SM5模型估算精度最高(RMSE、MAE、R2、Ens、GPI的中位線分別為1.122 mm/d、0.838 mm/d、0.716、0.609、0.166)。

(4)從全區(qū)來看，不同模型估算結(jié)果的RE空間分布趨勢基本一致，BA-LSTM模型在全區(qū)的估算精度均最高，RE僅為0.01%～1.75%。

綜上所述，BA-LSTM模型可作為魯中地區(qū)ET0簡化估算的推薦模型。