鄧安妮，康銀紅，王君勤，武劍飛

(1.四川農(nóng)業(yè)大學(xué)水利水電學(xué)院，四川 雅安 625014； 2.四川省水利科學(xué)研究院，四川 成都 610072)

0 引言

參考作物騰發(fā)量(ET0)是計算作物需水量的關(guān)鍵參數(shù)，ET0的準(zhǔn)確計算對灌溉系統(tǒng)設(shè)計運行、灌溉制度制訂及農(nóng)業(yè)用水管理具有重要意義[1]。ET0是指一種假想的參考作物冠層蒸發(fā)蒸騰速率，假設(shè)作物高度0.12 m，固定的葉面阻力70 s/m，反射率0.23，非常類似于表面開闊、高度一致、生長旺盛、完全覆蓋地面且不缺水的綠色草地的蒸發(fā)蒸騰速率[2-3]。

聯(lián)合國糧農(nóng)組織(FAO)推薦的Penman-Monteith(PM)方法是作為唯一定義和計算ET0的標(biāo)準(zhǔn)方法，但該方法要求氣象資料過多且計算冗雜，對于氣象資料不完整的地區(qū)用PM公式計算ET0會受到限制[4]。因此，尋找更為簡便、數(shù)據(jù)要求較低的方法非常有必要，國內(nèi)外學(xué)者提出了許多適用于不同地區(qū)的簡便模型，如溫度法、輻射法、綜合法和蒸發(fā)皿法等，同時對這些公式進(jìn)行了修正[5]。朱瀟梟等[6]采用Priestley-Taylor(PT)法、Irmak-Allen(IA)法和Hargreaves(Har)法對海南省ET0進(jìn)行計算，發(fā)現(xiàn)PT公式的精度最優(yōu)。李志[7]研究了Turc、Makkink(Mak)等7種方法在黃土高原的適用性，結(jié)果表明FAO24 Blaney-Criddle(FAO24 BC)法和Har法適用性最好。XYSTRAKIS F等[8]對Mak等13種方法在希臘南部進(jìn)行適用性分析，最終認(rèn)為輻射法更適用于計算ET0。胡慶芳等[9]針對Har法進(jìn)行了適用性評價和修正。李晨等[10]在四川省對具有代表性的6種簡易算法48 Penman(48PM)法、Har法、PT法等進(jìn)行計算精度對比，得出推薦在東部盆地區(qū)使用Har法，在盆周山地區(qū)、川西南地區(qū)與川西高原區(qū)使用PT法計算ET0的結(jié)論。張戴煒等[11]基于貝葉斯原理改進(jìn)Har模型，ET0計算精度顯著提高。李丹陽等[12]在四川地區(qū)利用線性擬合獲得Har法修正參數(shù)，有效減小了計算結(jié)果的誤差。TEMESGEN B等[13]比較了風(fēng)速對Har法的影響，發(fā)現(xiàn)低風(fēng)速情況下ET0計算誤差低于強風(fēng)速情況。XU C Y等[14]在瑞士地區(qū)評價了Abtew、Har、Turc等5種輻射法，并對各方法進(jìn)行參數(shù)校正，研究表明校正后的方法精度有所提高。以上研究表明，由于各地區(qū)地理環(huán)境不同，ET0計算方法的適用性具有差異。因此，在實際應(yīng)用時需要進(jìn)行適用性評價和參數(shù)校正。

四川省瀘州市降雨量時空分布不均，限制了當(dāng)?shù)剞r(nóng)業(yè)生產(chǎn)發(fā)展，提高該地區(qū)的農(nóng)業(yè)用水效率、確定合理的灌溉規(guī)模、制定科學(xué)的灌溉制度尤為重要[15]。本文針對瀘州市納溪區(qū)、敘永縣ET0的計算方法進(jìn)行適用性分析與參數(shù)修正，使用5種方法包括PM法、溫度法中的Har法、FAO24 BC法、輻射法中的Mak法和PT法計算瀘州市的ET0，且以PM法得到的ET0作為標(biāo)準(zhǔn)，評價了其余4種簡便方法的適用性，并基于線性修正和貝葉斯公式對簡便方法修正，以期進(jìn)一步提高各方法的計算精度，為該地區(qū)水資源管理提供理論依據(jù)和數(shù)據(jù)支持。

1 材料與方法

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

瀘州市位于27°39′～29°20′N、105°08′41″～106°28′E，屬亞熱帶濕潤氣候區(qū)，氣候溫和，雨量充沛但四季分明，受四川盆地地形影響，該地全年少有大風(fēng)，多為0～2 m/s的微風(fēng)，年平均氣溫17.5～18.0 ℃，年際之間的變化為16.8～18.6 ℃[16]。瀘州市的納溪區(qū)地處四川盆地南緣丘陵與低山區(qū)之間，海拔230～960 m，年降雨量1 200 mm；敘永縣位于四川盆地南緣，海拔247～1 902 m，年降雨量1 172.6 mm，都屬于濕潤區(qū)[17-18]。

本文選取瀘州市納溪和敘永兩個氣象站1967—2017年的逐日氣象資料來計算ET0，氣象資料包括平均氣溫、平均相對濕度、日照時間、日最低氣溫、日最高氣溫和風(fēng)速。數(shù)據(jù)來源于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http：//cdc.cma.gov.cn/home.do)。

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

(1)PM綜合方程式[19]。

(1)

式中ET0(PM)——參考作物騰發(fā)量，mm/d

Δ——飽和水汽壓與溫度關(guān)系曲線的斜率，kPa/℃

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

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

γ——濕度計常數(shù)，kPa/℃

T——空氣平均溫度， ℃

u2——地面2 m高處的風(fēng)速，m/s

es——溫度為T時的飽和水汽壓，kPa

ea——溫度為T時的實際水汽壓，kPa

es-ea——水汽壓差，kPa

(2)Har法在缺少輻射資料的地區(qū)得到廣泛的應(yīng)用，公式如下[20]。

(2)

式中ET0(Har)——應(yīng)用Har法計算的參考作物騰發(fā)量，mm/d

Ra——天頂輻射，MJ/(m2·d)

λ——汽化潛熱，MJ/kg，取值2.45

Tmax——日最高氣溫， ℃

Tmin——日最低氣溫， ℃

(3)FAO24 BC法可用于只有氣溫觀測資料的地區(qū)，且只能用于計算月或更長時段的ET0，公式如下[21]。

ET0(FAO24 BC)=a+bp(0.46T+8.18)

(3)

a=0.004 3RHmin-n/N-1.41

b=a0+a1RHmin+a2n/N+a3ud
+a4RHminn/N+a5RHminu2

a0=0.819，a1=-0.004 09，a2=1.071

a3=0.065 6，a4=-0.005 97，a5=-0.000 597

式中ET0(FAO24 BC)——月平均參考作物蒸發(fā)蒸騰量，mm/d

p——月內(nèi)日平均可能日照時間占全年可能日照時間的百分比

RHmin——最低相對濕度， %

n——實際日照時間，h

N——最大可能的日照持續(xù)時間，h

ud——地表以上dm處的風(fēng)速，m/s

(4)Mak法屬輻射法范疇，是在太陽輻射量基礎(chǔ)上提出的，該方法在寒冷地區(qū)的適用性得到廣泛驗證，公式如下[22]。

(4)

式中ET0(Mak)——應(yīng)用Mak法計算的參考作物騰發(fā)量，mm/d

Rs——太陽輻射或短波輻射，MJ/(m2·d)

(5)PT法是假定濕潤環(huán)境條件下，忽略空氣動力學(xué)項并引進(jìn)參數(shù)而得出的簡化算法，公式如下[23]。

(5)

式中ET0(PT)——PT法計算的參考作物騰發(fā)量，mm/d

1.3 修正方法

1.3.1 線性修正

PM法和Har、PT法計算的ET0結(jié)果有一定偏差，實際中常建立以下線性關(guān)系式進(jìn)行修正，并確定修正參數(shù)a、b[12，24]。

(6)

a、b——線性修正參數(shù)

1.3.2 貝葉斯公式修正

貝葉斯理論方法包括先驗分布和后驗分布，可以用已有數(shù)據(jù)來估計某事物將來發(fā)生的可能性。應(yīng)用到修正參數(shù)中時，可以將新的數(shù)據(jù)代入公式得到結(jié)果，并將此作為下一次估算的依據(jù)，從而獲得新的后驗分布，以此得到更加準(zhǔn)確的估計參數(shù)，提高模型準(zhǔn)確度[25]。

利用貝葉斯公式對Har法和PT法進(jìn)行修正，由文獻(xiàn)[25]知，改進(jìn)后的方法可以表示為

(7)

(8)

C*、α*——修正后的經(jīng)驗系數(shù)，可由貝葉斯公式進(jìn)行修正

1.4 評價指標(biāo)

以PM法計算結(jié)果作為依據(jù)，采用均方根誤差RMSE、誤差百分比PE、平均偏差MBE、確定系數(shù)R2共4種評價指標(biāo)來評價上述計算ET0的不同方法，確定其在不同地區(qū)的適用性[26]。各評價指標(biāo)的計算公式如下。

(9)

(10)

(11)

(12)

式中ETi——不同方法計算的ET0，mm/d

ET0(PM)——PM綜合方程式計算的ET0，mm/d

n——樣本總數(shù)

其中，RMSE、PE、MBE越接近于0，R2越高，則該方法與PM法一致性越好，計算精度越高、誤差越小。

2 結(jié)果與分析

2.1 月ET0結(jié)果分析

以PM法計算的ET0作為參考標(biāo)準(zhǔn)值，對各簡便方法進(jìn)行適用性評價，PM法和4種ET0計算模型的逐月ET0結(jié)果如圖1所示。Har法和PT法在納溪站的1—4月、10—12月時計算結(jié)果基本重合，其中5—8月的ET0值與標(biāo)準(zhǔn)值相差較大，誤差范圍為0.7～0.8 mm/d；而在敘永站，1—3月、11—12月結(jié)果基本接近，5—6月的ET0值相差較大，誤差為0.55～0.89 mm/d；在兩站點使用Mak法和FAO24 BC法均得到較大誤差。各方法計算得到的誤差都在7月和8月達(dá)到最大，1月和12月最小，此結(jié)論與李晨等[10]計算四川省不同地區(qū)的ET0得到的變化趨勢保持一致。

圖1 Har法、FAO24 BC法、Mak法、PT法與PM法的月ET0對比Fig.1 Comparison of monthly ET0 calculated by Har，F(xiàn)AO24 BC，Mak，PT methods and PM method

2.2 各方法適用性評價

根據(jù)月ET0計算結(jié)果，利用評價指標(biāo)進(jìn)一步對簡化算法進(jìn)行誤差分析，結(jié)果如表1所示。在納溪站，溫度法中Har法表現(xiàn)出較低的誤差值，其RMSE、PE、MBE、R2分別為1.09 mm/d、10.02%、0.26 mm/d及0.64，與PM標(biāo)準(zhǔn)值最接近，而FAO24 BC法計算出的結(jié)果有較大誤差，誤差百分比超過70%，計算效果最不理想。輻射法中PT法比Mak法更好，其RMSE、PE、MBE、R2分別為0.53 mm/d、15.27%、0.40 mm/d及0.97。對于敘永站而言，Har法和PT法同樣表現(xiàn)出較低的誤差，RMSE、PE、MBE、R2分別為0.94 mm/d、12.94%、0.35 mm/d及0.79和0.55 mm/d、11.18%、0.31 mm/d及0.94，但Har法較PT法的RMSE偏大，R2偏小。綜上所述，Har法和PT法在氣象資料缺少的情況下，更接近PM法，結(jié)果表現(xiàn)更理想。但Har法、PT法計算的ET0普遍都高于標(biāo)準(zhǔn)值，這與劉靜霞等[27]、YODER R E等[28]計算在濕潤地區(qū)都得出高估ET0的研究結(jié)果相同。

表1 Har法、FAO24 BC法、Mak法、PT法與PM法的評價指標(biāo)結(jié)果Tab.1 Evaluation index results of Har，F(xiàn)AO24 BC，Mak，PT methods and PM method

2.3 方法修正分析

為進(jìn)一步提高ET0簡便公式在瀘州市計算的準(zhǔn)確性，選用精度較高的Har、PT法，利用1967—2006年ET0值對各方法進(jìn)行修正，利用2007—2017年對改進(jìn)后的方法進(jìn)行驗證計算，分析Har、PT公式改進(jìn)前后計算ET0的精度。

2.3.1 線性修正

使用線性關(guān)系對Har、PT法進(jìn)行修正，結(jié)果如表2所示[29-30]。

從表2可以看出，兩個站點的相關(guān)系數(shù)都在0.98以上，Har法修正參數(shù)a、b分別是0.878 1、0.088 6和0.866 4、0.059 0，PT法修正參數(shù)a、b分別是0.846 1、0.064 5和0.897 3、0.005 5。因此，在氣象資料缺少的情況下，可以使用表2的修正式對ET0進(jìn)行修正。

表2 Har法和PT法修正參數(shù)值及修正式Tab.2 Corrected parameter values and formulas of Har and PT methods

2007—2017年修正前、后計算得到結(jié)果對比值如圖2所示。從圖2對比中可以觀察到，各方法修正后的計算值明顯比修正前更接近PM標(biāo)準(zhǔn)值，在兩站點使用修正后的PT法最終計算的ET0幾乎與標(biāo)準(zhǔn)值完全擬合，相差最大的值分別出現(xiàn)在3月和6月，僅相差0.15和0.11 mm/d，修正后的Har法計算結(jié)果除5月和6月相差最大，分別為0.33、0.38 mm/d之外，其余月份得到的ET0基本接近標(biāo)準(zhǔn)值，可見修正前后PT法的誤差下降了80%，而Har法下降了50%。

圖2 2007—2016年線性修正前后月ET0比較Fig.2 Comparison of monthly ET0 before and after linear correction from 2007 to 2016

對Har、PT法修正后，大大減小了各方法的誤差，且發(fā)現(xiàn)修正后的PT法比Har法更適合計算瀘州市的ET0。

2.3.2 貝葉斯公式修正

修正數(shù)據(jù)同樣使用1967—2006年兩站點的ET0值，使用貝葉斯公式修正后Har法、PT法在納溪和敘永的經(jīng)驗系數(shù)C*、α*分別為0.002 276、1.206 350及0.002 254、1.176 490。修正后ET0結(jié)果值對比如圖3所示，利用貝葉斯公式修正參數(shù)后沒有明顯改善Har方法和PT方法計算出的ET0，較PM標(biāo)準(zhǔn)值仍是偏大，相差最大的ET0值僅下降0.04～0.20 mm/d。

圖3 2007—2017年貝葉斯公式修正前后月ET0比較Fig.3 Comparison of monthly ET0 before and after Bayesian formula correction from 2007 to 2016

通過對比兩種方法修正前后的評價指標(biāo)(表3)可以看出，兩種修正方法的R2無變化，基于貝葉斯公式在兩個站點修正后Har法的RMSE、PE、MBE分別下降了1.4%～3.2%、11.7%～17.5%、11.5%～17.2%，修正后PT法的RMSE、PE、MBE分別下降了14.3%～22.6%、32.2%～65.4%、32.3%～64.5%，使用線性關(guān)系修正后的Har法和PT法的RMSE、PE、MBE在兩個站點都降低50%以上，貝葉斯公式的修正幅度比利用線性關(guān)系得到的小得多。

表3 Har法、PT法修正前后評價指標(biāo)對比Tab.3 Comparative analysis of evaluation indexes before and after correction of Har and PT methods

線性關(guān)系比貝葉斯公式更適合修正Har法和PT法，且線性修正PT法優(yōu)于修正Har法。

3 結(jié)論

本文選取納溪和敘永站點1967—2017年氣象資料，以PM公式計算的ET0作為標(biāo)準(zhǔn)值，分析評價了4種計算模型的適用性，并對各方法進(jìn)行修正，得出如下結(jié)論。

(1)根據(jù)逐月ET0計算結(jié)果可知，4種方法計算的ET0變化趨勢與PM法相同，都呈先上升后下降趨勢，為開口向下的拋物線。

(2)由誤差分析可知，Mak法和FAO24 BC法計算結(jié)果與PM法相比誤差較大，RMSE在0.7～2.7 mm/d、PE在22%～78%。Har法和PT法的計算結(jié)果與PM法的計算結(jié)果最接近，RMSE在0.5～1.1 mm/d、PE在10%～15%，誤差相對較小，可用于計算當(dāng)?shù)谽T0。

(3)為進(jìn)一步提高ET0簡便公式在瀘州市的計算精度，對Har法和PT法使用線性關(guān)系和貝葉斯公式進(jìn)行修正，線性修正后的Har法和PT法在兩個站點誤差分別下降50%、80%左右，而貝葉斯公式的修正效果不明顯，在對比中發(fā)現(xiàn)使用線性關(guān)系比貝葉斯公式修正效果更好，且修正后PT法比Har法在瀘州市的適用性更好，可以提高ET0計算準(zhǔn)確度。