農業分區框架下Hargreaves-Samani公式的逐月回歸修正

夏興生 潘耀忠 朱秀芳 張錦水

(1.北京師范大學遙感科學國家重點實驗室，北京 100875；2.北京師范大學遙感科學與工程研究院，北京 100875；3.青海師范大學地理科學學院，西寧 810016；4.北京師范大學地表過程與資源生態國家重點實驗室，北京 100875)

0 引言

參考作物需水量(Reference crop evapotranspiration，ET0)是計算作物需水量、制定農田水利管理策略的核心參數，對農業水資源的開發利用和農業生產潛力的評估具有重要意義[1]。目前，常用的參考作物需水量計算方法有Penman-Monteith公式法(PM公式)[2]、Makkink方法[3]、Hargreaves-Samani公式(HS公式)[4]、Droogers and Allen方法[5]等10余種。其中，PM公式因符合氣候學和空氣動力學機理而被聯合國糧農組織(FAO)推薦為計算ET0的標準方法[2]。然而，PM公式要求研究或應用區域具有相對完善的氣象觀測數據，但在實際研究和應用中，全球范圍內均存在氣象站點觀測數據的時空差異，可獲得的數據很難滿足PM公式的輸入要求[6]。因此，結果精確、可行性強、適用范圍廣的ET0簡化計算方法得以發展，例如Makkink方法[3]、HS公式[4]、Droogers and Allen方法[5]等。PM公式則作為標準算法為這些簡化算法提供標準值[6]。

國內外的研究顯示[7-11]，基于溫度和天頂輻射的HS公式與PM公式計算的ET0具有較好的吻合度，是一種計算ET0較好的簡化氣候學經驗公式，FAO也將其作為氣候、氣象資料缺失條件下計算ET0的推薦方法[2]。但是，考慮到HS公式是HARGREAVES和SAMINA在干旱半干旱條件下觀測和實驗得出的經驗方程[12]，國內外眾多學者針對不同研究區或應用區基于HS公式的精度評價、參數校正等問題進行了大量研究[13-19]，結果均表明，直接使用HS公式會產生誤差，并建議在新的地區使用HS公式計算ET0時，需對其進行地域化修正。對此，FAO也給出了通過線性回歸分析方法進行HS公式地域化修正的建議方案[2]。

現有研究中主要包括基于溫度的HS校正方法[20]、基于降水的HS校正方法[5]、基于全局優化的HS校正方法[21-22]、基于濕潤條件下的局部校正方法[14]等，這些方法不僅要考慮地形、氣候要素對每一個參數的影響[23-25]，還要考慮參數之間的組合問題[26-27]，且能夠考慮到其他地理因子對ET0的影響。但是，在中國復雜的農業地域類型條件下，仍以小區域研究為主[27-29]，對應研究結果的時空尺度差異較大，因此，并未得到廣泛的推廣應用。對于FAO建議的線性回歸修正方案，有關研究案例卻鮮有報道，僅有李丹陽等[1]基于四川地區逐日氣象資料對HS公式進行了線性回歸修正，結果表明，修正后的HS公式可以反映該地區的ET0實際狀況。

本研究以中國九大陸地綜合農業區38個子區為研究區，采用FAO建議的線性回歸分析方法逐月逐區進行HS公式的修正，驗證在中國農業綜合分區框架下HS公式線性回歸修正方案的適用性，并嘗試確定全國不同農業區逐月的HS公式最佳線性回歸校正系數，為大規模實踐中快速獲得較高精度的ET0提供參考，豐富HS公式修正應用的案例，以期為未來更為精準的ET0空間化研究奠定基礎。

1 參考作物需水量計算方法

1.1 Penman-Monteith公式

PM公式的原型是1948年PENMAN采用能量平衡與質量傳遞相結合的方法以日照、溫度、濕度和風速4個標準氣象數據推導的水面蒸發量方程[30]。后經MONTEITH引入阻力因素改進并將其擴展至計算作物表面的騰發量[31]。1990年開始，經FAO專家咨詢委員會會同國際灌排委員會和世界氣象組織對其驗證和修正后，誕生了FAO推薦的PM公式[2]，并作為新的確定參考作物需水量的標準方法，其計算式為

(1)

(2)

(3)

(4)

(5)

(6)

Rn=Rns-Rnl

(7)

Rns=(1-α)Rs

(8)

(9)

Rso=(0.75+2×10-5Z)Ra

(10)

式中ET0-PM——基于PM公式計算的參考作物需水量，mm/d

ea——平均水汽壓，kPa/d

u2——2 m高處風速，m/(s·d)

T——大氣溫度，℃

Tmean——2 m高處日平均氣溫，℃/d

Tmax——計算時間區間的平均最高氣溫，℃/d

Tmin——計算時間區間的平均最低氣溫，℃/d

es——以日為區間計算的飽和水汽壓，kPa/d

e(T)——空氣溫度T時的飽和水汽壓，kPa

Δ——氣溫T下的飽和水汽壓曲線的傾率，kPa/(℃·d)

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

P——平均大氣壓，kPa/d

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

CP——常壓下的比熱容，取1.013×10-3MJ/(kg·℃)

ε——水蒸氣分子量與干燥空氣分子量的比，取0.662

G——以月為計算區間的土壤熱通量，MJ/(m2·d)

Rn——參考作物表面上的凈輻射，MJ/(m2·d)

Rns——地表凈太陽短波輻射，MJ/(m2·d)

Rs——地表太陽總輻射，MJ/(m2·d)

α——發射率或冠層發射系數，以草為假想的參考作物時，取0.23

Rnl——凈輸出長波輻射，MJ/(m2·d)

σ——Stefan-Boltzmann常數，取4.903×10-9MJ/(K4·m2·d)

Tmax,k——24 h內最高絕對溫度，K

Tmin,k——24 h內最低絕對溫度，K

Rso——晴空條件下到達地表的太陽總輻射，也稱晴空太陽輻射，MJ/(m2·d)

Rs/Rso——相對短波輻射值，其值小于等于1.0

Z——測站的海拔，m，因2×10-5Z實際值較小，在計算中可忽略不計

Ra——天頂輻射，由FAO建議的以日為時段的天頂輻射程序計算獲得，MJ/(m2·d)

以月為計算區間的土壤熱通量計算式為

Gmonth,i=0.07(Tmonth,i+1-Tmonth,i-1)

(11)

式中Gmonth,i——第i月的土壤熱通量，MJ/(m2·d)

Tmonth,i-1——第i-1月的大氣平均溫度，℃/d

Tmonth,i+1——第i+1月的大氣平均溫度，℃/d

1.2 Hargreaves-Samani公式

HS公式是HARGREAVES和SAMANI[4]1985年基于干旱半干旱地區的研究提出的經驗公式，FAO推薦其為氣象觀測數據缺失條件下的ET0簡化計算方法，同時也建議在有條件的地區使用該公式計算ET0時要基于月或年的經驗系數進行修正，相關算式為

ET0-HC=a+bET0-HS

(12)

其中

ET0-HS=α(Tmean+17.8)(Tmax-Tmin)0.5Ra

(13)

式中ET0-HC——校正后的ET0-HS，mm/d

a、b——校正系數

ET0-HS——基于HS公式計算的參考作物需水量，mm/d

α——常數，取0.002 3

校正系數a和b通過回歸分析或目測的方法擬合。本研究就是基于線性回歸分析方法逐月求取a、b系數在中國各農業區的取值，達到校正HS公式的目的。

2 數據及預處理

2.1 數據來源

根據上述PM公式和HS公式的參數算法，本研究使用的數據包括氣象站點數據和農業區劃數據，具體如下：

(1)氣象站點數據：逐月的平均氣壓P、平均最低氣溫Tmin、平均最高氣溫Tmax、平均風速u、平均水汽壓ea，均來自中國地面氣候資料月值數據集，包含756個站點的數據；月總太陽輻射Rs，來自中國輻射月值數據集，包含130個站點的數據；站點經緯度數據，來自中國地面氣候資料月值數據集的元數據文件。兩個數據集均下載自中國氣象數據網(http:∥data.cma.cn/)，時間尺度1957—2016年。

(2)中國農業綜合區劃數據：來自全國農業委員會發布的《中國農業綜合區劃圖》(圖1，圖中A表示東北區：A1興安嶺區，A2松嫩三江平原區，A3長白山山地區，A4 遼寧平原區；B表示內蒙古及長城沿線區：B1內蒙古北部區，B2內蒙古中南部，B3長城沿線區；C表示黃淮海區：C1燕山太行山山麓平原區，C2冀魯豫低洼平原區，C3黃淮平原區，C4山東丘陵區；D表示黃土高原區：D1晉東豫西丘陵區，D2汾渭谷地區，D3晉陜甘黃土丘陵區，D4隴中青東丘陵區；E表示長江中下游區：E1長江下游平原區，E2豫鄂皖低山區，E3長江中游平原區，E4江南丘陵區，E5浙閩丘陵區，E6南嶺丘陵區；F表示西南區：F1秦嶺大巴山區，F2四川盆地，F3川鄂湘黔邊境區，F4黔桂高原山地區，F5川滇高原山地區；G表示華南區：G1閩南粵中區，G2粵西桂南區，G3滇南區，G4瓊雷及南海諸島區，G5臺灣；H表示甘新區：H1蒙寧甘區，H2北疆區，H3南疆區；I表示青藏區：I1藏南區，I2川藏區，I3青甘區，I4青藏高寒區)。

圖1 農業綜合區劃及氣象站點位置示意圖

2.2 數據預處理

(1)將月總太陽輻射數據轉換為逐月日均值，其他氣象數據校正為PM公式要求的2 m高度標準數據，并統一數據單位尺度。

(2)通過站點編號將Rs、點經緯度及其他氣象數據一一對應關聯，共得到124個交集站點(圖1)。

(3)根據元數據說明剔除無效數據，并基于FAO[2]建議的以日為時段的天頂輻射計算程序，求得各站點逐月的日平均Ra，以Ra>Rs的理論關系為規則進行數據質量檢查，剔除Rs異常值。

(4)將數據分為系數a、b的確定和驗證兩部分，即1957—2010年的數據進行系數a、b的回歸確定，2011—2016年的數據則用來進行系數a、b的應用驗證。

(5)以《中國農業綜合區劃圖》上的經緯度和中國行政矢量數據為基準，進行空間校正并數字化得到農業分區矢量數據。

3 技術路線

本研究過程主要包括ET0-PM和ET0-HS的計算、HS公式校正系數的計算、誤差分析3部分，具體技術路線如圖2所示。

(1)ET0-PM和ET0-HS計算

按照PM公式和HS公式分別逐站點逐年逐月計算得到ET0-PM和ET0-HS，其中ET0-PM作為本研究的真值。

(2)HS公式校正系數計算

以各站點多年的逐月ET0-PM為因變量，以與之對應的ET0-HS為自變量，以最小二乘回歸算法計算得到各站點各月份的系數a、b。

(3)誤差分析

以農業區劃矢量數據區域統計得到各農業子區內系數a、b的站點平均值，以式(13)計算并統計得到2011—2016年逐年逐月各農業子區的ET0-HC，通過計算比較ET0-HC、ET0-HS的6年平均相對誤差，說明HS校正系數在中國各農業區的適用性，并篩選出系數a、b的最佳建議取值。

4 結果與分析

4.1 HS公式校正系數的回歸結果

圖3是逐月統計的各區1957—2010年ET0-PM與ET0-HS相關系數均值，其中，D1區因整體沒有有效的站點觀測數據，故D1區缺失，下同。從圖可以看出，大部分農業區在年內有超過一半的月份ET0-PM、ET0-HS相關系數均超過0.6，相關性顯著，只有G4、H3和I1相關系數相對較小的月份較多，可能是因為有效數據較少且數據時間序列不連續導致ET0-PM和ET0-HS存在隨機誤差。整體上，ET0-PM與ET0-HS具有較好的相關性，可以直接基于最小二乘回歸方法計算HS公式的校正系數a、b。

圖3 逐月統計的各區1957—2010年ET0-PM、ET0-HS相關系數均值

圖4是基于最小二乘回歸方法逐月計算的各區系數a、b平均值。從圖可以明顯看出，九大農業區系數a、b在時間上存在明顯的變化規律。其中，系數a統一表現為冬季月份最大，夏季月份最小，春、秋季月份為過渡變化時段的規律；系數b則在年內統一表現比較平穩，集中在1.0～2.0之間。a、b的相對大小卻沒有明顯的對應關系。由此可見，系數a的變化在很大程度上決定了ET0-HS的校正精度，也說明在中國ET0-HS和ET0-PM之間存在明顯的季節性差異，且夏季差異較大，冬季差異較小。這一差異顯然與影響ET0變化的主導因子隨中國“雨熱同期”的季風氣候變化規律有關。HS公式只考慮了氣溫對ET0的影響，PM公式則是綜合考慮了氣溫、風速、濕度等各項氣候條件，當影響ET0的主導因子隨著季節發生變化時，ET0-HS和ET0-PM之間的差異性也隨之發生變化。二者在冬季差異小，說明冬季影響ET0的主導因子為氣溫，在夏季差異大，說明夏季影響ET0的主導因子不僅僅是氣溫，還有其他因子，但主要包含哪些主導因子則需要進一步探討。

在圖4所示的九大農業區內部，各子區之間的系數a、b則存在明顯的差異，且沒有規律可循。例如在A3、C2、H3和I4區，系數a的年內變化幅度較其所在大區的其他子區明顯要??；G4和I1區的系數a則與整體變化規律完全不同，大部分月份a均大于0，且波動變化明顯。這一結果說明九大農業區內部的局部地理因子可能也會引起ET0-HS和ET0-PM之間的差異，還需要在精細化的尺度上再進行深入研究。

理論上，最小二乘回歸結果a越接近于0，b越接近于1，說明ET0-HS就越接近ET0-PM。但綜合圖4的分析結果，在中國九大農業區，ET0-PM和ET0-HS之間整體上存在差異，特別是夏季月份的差異還較大，不僅驗證了引言中說明的已有研究結果表明直接使用HS公式計算的ET0-HS會產生誤差的結論，也驗證了FAO建議進行HS公式地域化修正的必要性[2]，即，在實際應用中，應以盡可能得到的觀測數據進行ET0-HS的本地化校正。

圖4 校正系數a、b逐月的回歸結果

4.2 誤差分析和最優系數的確定

基于上述的系數a、b回歸結果，以式(13)計算并統計得到2011—2016年逐年逐月各農業子區的站點平均ET0-HC，進一步以ET0-PM站點均值為真值，計算了ET0-HC、ET0-HS的6年逐月平均相對誤差，結果如圖5所示。由圖5可以看出，ET0-HS和ET0-HC相對于ET0-PM均存在不同程度的誤差，但是，ET0-HS的相對誤差年內較為分散，特別是在A、B兩區的分散度較大，進一步驗證了ET0-HS和ET0-PM之間存在顯著差異，而ET0-HC的相對誤差則較為集中，說明校正后的ET0-HS可以縮小ET0-HS和ET0-PM的差異。在具體數據上，ET0-HS相對誤差集中在50%以內，最大值高達250%(A1區)，而ET0-HC的相對誤差則集中在20%以內，最大值也未達到40%。因此，以線性回歸方法得到HS公式校正系數a、b能夠代表除溫度外的其他因子對ET0的影響，從而使得ET0-HS的校正結果能夠顯著縮小其對應的年內相對誤差范圍，是實踐中可行的方案。

圖5 ET0-HS校正前后相對誤差各區站點均值年內分布

表1是ET0-HC和ET0-HS在2011—2016年不同區域逐月的相對誤差比較結果，其中小于0的數字表示ET0-HC優于ET0-HS，值越小表示ET0-HC的相對精度越高，大于等于0的數字則相反。從表1可以看出，在具體的驗證應用中，并不是所有農業區在所有月份均是ET0-HC優于ET0-HS，而是大部分農業區均不同程度存在ET0-HS優于ET0-HC的情況。具體來說，A1、B1、B2、F3、H1、I1、I3區全年均表現為ET0-HC優于ET0-HS，而C2、D4、G4區有11個月，A2、C3、C4、E2、F1、I4區有10個月，A3、C1、D3、F2、H2、H3區有9個月，B3、E6、G1區有8個月，A4、F5區有7個月，E4、I2區有6個月，E1、E3、E5、F4區有5個月，D2、G2有4個月，G3僅有1個月ET0-HC優于ET0-HS，這些區域剩余的時間段均ET0-HS優于ET0-HC。對于這一結果的原因，最大的可能是系數a、b是長時間序列數據線性回歸的結果，本質上是一種整體的趨勢擬合結果，在短期內的驗證應用中，存在氣象/氣候條件的隨機性誤差，故而造成了在具體的驗證應用中存在ET0-HS優于ET0-HC的情況，但具體情況有待進一步分析。

表1 2011—2016年ET0-HS校正前后精度比較結果

綜上誤差分析結果，整體上ET0-HC要優于ET0-HS，但在各農業子區的具體應用中，ET0-HC的相對誤差結果則參差不齊。因此，優選使用ET0-HS的修正(HS公式的回歸校正)系數才是當前的最佳方案，即根據本研究的驗證結果，整體上，中國各農業區使用回歸得到的本地化校正系數a、b進行ET0-HS的修正，但在校正結果較差的時段，仍然直接使用ET0-HS作為ET0-HC簡化計算的最終結果(即a=0,b=1)，從而使HS公式計算的ET0-HC在不同農業區的各個時段均達到最優。基于這一方案，本研究最終推薦使用的中國各農業區HS公式逐月的校正系數如表2所示。

表2 中國各農業區HS公式逐月的校正系數建議值

5 討論

ET0是影響作物需水量的重要因素，也是區域灌溉制度和水資源管理策略制定的重要依據。本研究基于可獲得的數據，逐月回歸校正了中國各農業管理區利用HS公式計算的ET0-HS，結果表明基于FAO建議的HS公式線性回歸修正方案可行性強、適用范圍廣，且具有較好的精度，對大規模實踐中快速獲得較高精度的ET0以指導農業水資源的開發利用、農業生產潛力的評估和強化節水農業的發展具有重要的參考意義。但從精度分析和驗證應用的結果來看，還需進一步討論如下問題：

(1)由圖5可知，ET0-HC的相對誤差年內變化范圍已經顯著小于ET0-HS，說明ET0-HS和ET0-PM存在線性關系，基于最小二乘回歸逐月求得的ET0-HS校正系數a、b能夠在應用中使ET0-HS更接近于ET0-PM，也就是說系數a、b已經反映了除溫度外的其他氣象因子對ET0的影響，可在實踐中應用，也驗證了HS公式線性校正方案在中國的可行性。但是，分析ET0-PM的計算過程可知，在Rso的計算中，使用了式(10)，該式中的常數0.75是晴空條件下Ra到達地面的比例?，F實中，在大氣組分和厚度時空變化的影響下，這一比例數值可能存在時空差異，而且在中國并未發現有研究或行業部門發布該式中的本地化參數，因此，基于此計算的ET0-PM本身就存在系統誤差，只能稱作準真值，并不是真正的標準值。所以，后續研究中有必要探討計算Rso的公式參數在中國各農業區的具體取值。

(2)由表1可知，大部分農業區均不同程度存在ET0-HS優于ET0-HC的情況。存在這一結果的原因，除了長時間序列數據的線性回歸增強了氣候變化引起的誤差和短期驗證應用存在的隨機誤差外，如從數據的角度分析，可能數據質量也是這一結果的一個影響因子。本研究數據的獲取時間尺度是1957—2016年。因氣象觀測中存在的不可控因素，站點數據均存在不同程度缺失，即，實際上站點的數據時間序列是不連續的，只是在研究中將其假設為時間序列數據。另外，不同站點的有效觀測時間尺度也不一致，所以，最后各站點回歸求得系數a、b時的樣本量不一致。之后，又以區域的站點平均值來代表整個研究區的取值，但各農業區內站點密度明顯不均(圖1)，系數a、b取站點的平均值沒有考慮其他地理因子的影響，所以，以區域內站點的平均結果代表整個區域的取值也可能進一步致使結果存在偏差。因此，未來不僅要考慮驗證不同氣候尺度下ET0-HS的校正精度，以確定最優的氣候尺度從而削弱氣候變化對ET0的影響，還需要討論如何進一步弱化現有數據不連續的影響，而最優系數在具體驗證年份及未來應用中是否存在變化，也值得進一步的探討。

6 結論

(1)在年內，大部分農業區超過一半的月份ET0-PM與ET0-HS的相關系數超過0.6，因此，可以直接以ET0-PM為真值進行ET0-HS的回歸校正。

(2)九大農業區回歸校正得到的系數a整體上存在顯著的季節變化規律，冬季月份最大，夏季月份最小，系數b則表現較為平穩。從a、b變化可知，ET0-PM和ET0-HS之間存在差異，且季節性明顯，在實際應用中進行ET0-HS校正十分必要。

(3)基于驗證應用年份的結果表明，校正前后的ET0-HS相對于ET0-PM均存在不同程度的誤差，但校正后的ET0-HS的誤差范圍相比校正前的ET0-HS誤差范圍已經顯著縮小，因此，進行HS公式的線性回歸校正在中國各農業區具有一定的適用性，且是實踐中簡便、可行的方案。

(4)在各農業子區的具體驗證應用中，ET0的相對誤差結果表現參差不齊，因此，校正系數a、b的優選使用才是以HS方程快速計算的ET0在不同農業區的各個時段均達到最優的最佳方案。