西藏高原灌區(qū)參考作物蒸散量模型的適用性研究

蒙 強，劉靜霞，李玉慶，張文賢

(西藏農(nóng)牧學(xué)院水利土木工程學(xué)院，西藏 林芝 860000)

參考作物蒸散量(ET0)是作物系數(shù)法計算作物需水量、制定農(nóng)業(yè)灌溉制度及評價區(qū)域水資源的重要參數(shù)。近年來，隨著蒸散理論研究的深入，一系列基于溫度、輻射、風(fēng)速及綜合氣象參數(shù)的ET0計算模型被廣泛應(yīng)用于蒸散量的估算中。其中，F(xiàn)AO 56 Penman-Monteith(FAO 56 PM)模型綜合考慮了所有氣象因子對ET0的影響且具有較完整的理論，被國際公認為是ET0計算的標準方法，且常被用于修正及評價其他ET0模型的精度和適用性[1]。此外，國內(nèi)學(xué)者關(guān)于ET0計算方法及其適用性的研究結(jié)果已充分證明了FAO 56 PM法在我國具有較高的計算精度[2-4]，但其計算過程較為復(fù)雜且需要完整的氣象資料，故無法應(yīng)用于氣象資料缺失條件下的ET0計算。對此，學(xué)者發(fā)展了多種ET0的簡化模型，例如FAO 24Radiation、Hargreaves-Samani、Makkink等，但上述模型大多是針對特定的區(qū)域及氣候條件而提出的，已有的研究結(jié)果表明其計算精度及適用性存在一定的地域差異。例如，樊軍等[5]通過對比多種ET0模型在陜西省的適用性，指出Priestley-Taylor模型的模擬值最接近于FAO 56 PM模型的計算結(jié)果，F(xiàn)AO-24 Blaney-Criddle、FAO-24 Radiation、Hargreaves-Samani及Makkink 4種模型與FAO 56 PM模型的模擬值間存在明顯的差異。李志[6]研究認為FAO-24 Blaney-Criddle、Hargreaves-Samani模型在黃土高原區(qū)有較好的適用性，但FAO-24 Radiation、Priestley-Taylor、Makkink模型對ET0的模擬效果較差。余婷等[7]的研究結(jié)果表明，在我國西北地區(qū)Kimberly Penman、FAO 1979 Penman、FAO 24 Penman、FAO 1948 Penman、FAO 79 Penman、Hargreaves-Samani、Mc Cloud、Hargreaves-Samani、Priestley-Taylor-1、Priestley-Taylor-2、FAO-24 Radiation、Makkink、Irmark-Allen、Irmark 這14種模型的模擬精度存在明顯差異，其中FAO 1948 Penman、Hargreaves-Samani、Priestley-Taylor-1的模擬精度最高。徐冰[8]等關(guān)于西藏高寒牧區(qū)ET0模型適用性的研究結(jié)果表明，F(xiàn)AO 1979 Penman、Irmark-Allen模型模擬值較FAO 56 PM總體偏大，Priestley-Taylor模型的模擬值誤差最小，并指出Hargreaves-Samani模型可作為高海拔區(qū)域的ET0在氣象資料短缺條件下的簡化算法。

然而，西藏高原特殊的地形地貌及氣候環(huán)境給作物需水量的準確估算帶來了一定的困難，為降低作物需水量估算的不確定性，提高農(nóng)業(yè)灌溉用水效率，開展ET0計算模型的適用性評價具有重要的現(xiàn)實意義。但西藏地區(qū)有關(guān)ET0的理論研究起步較晚，已有的研究大多是關(guān)于ET0的時空分布及變化特性的分析，而關(guān)于ET0計算方法的分析對比及適用性探討的研究則明顯不足。因此，本研究選取西藏高原3座典型灌區(qū)為研究區(qū)域，基于長系列的灌區(qū)氣象數(shù)據(jù)，利用目前常用的FAO 17 Penman、Hargreaves-Samani、Makkink、Priestley-Taylor、Irmark-Allen 5種模型計算ET0值，并以FAO 56 Penman-Monteith公式的計算值為標準，對5種ET0計算模型的適用性進行評價，旨在為西藏高原灌區(qū)開展高效節(jié)水灌溉及優(yōu)化配置農(nóng)業(yè)水資源提供技術(shù)參考。

1 數(shù)據(jù)與方法

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

西藏的灌區(qū)主要分布在日喀則(2 416 座)、山南地區(qū)(1 957 座)。其中，滿拉、墨達、江北灌區(qū)分別位于日喀則、拉薩、山南地區(qū)，是西藏高原上的大型灌區(qū)，灌溉面積均超過2 萬hm2。根據(jù)《西藏農(nóng)業(yè)氣候資源區(qū)劃》結(jié)果，滿拉、墨達、江北灌區(qū)位于半干旱農(nóng)業(yè)區(qū)。本研究在綜合考慮了灌區(qū)面積、行政區(qū)劃、農(nóng)業(yè)種植等情況下，選取滿拉、墨達、江北3個灌區(qū)為本文的典型研究區(qū)域，各灌區(qū)基本資料見表1。滿拉、墨達、江北灌區(qū)附近的氣象站點分別為江孜、墨竹工卡和貢嘎，各站點的逐日氣象資料來源于中國氣象數(shù)據(jù)服務(wù)網(wǎng)。其中，江孜站點為1956-2016年61年逐日氣象資料；貢嘎及墨竹工卡站點為1991-2016年26年逐日氣象資料。

1.2 計算參考作物蒸散量的6種模型

(1)FAO 56 Penman-Monteith(FAO 56 PM)。FAO 56 Penman-Monteith模型是基于水汽擴散理論和能量平衡為基礎(chǔ)建立的并被FAO推薦為目前計算ET0的標準方法。該模型綜合考慮了所有氣象因子對ET0的影響且具有較完整的理論依據(jù)，常被用來評價其他ET0計算模型的精度，屬綜合法范疇[1]。國內(nèi)學(xué)者的研究成果已證明了該方法在我國具有較高的普適性[4]。因此，本研究以FAO 56 PM模型計算的ET0值為標準，對本文中其他5種模型進行適用性評價。

(1)

式中：ET0為參考作物蒸散量，mm/d；Δ為飽和水汽壓-溫度曲線斜率，kPa/℃；Rn為太陽凈輻射，MJ/(m2·d)；G為土壤熱通量，MJ/(m2·d)；γ為濕度計常數(shù)，kPa/℃；T為2 m高度處日平均氣溫，℃；es為溫度為T時的實際水汽壓，kPa；ea為溫度為T時的飽和水汽壓，kPa；U2為2 m高處的日平均風(fēng)速,m/s。以下與此公式中的參數(shù)、符號相同者意義均相同。

(2)FAO 17 Penman(FAO 17 PM)。FAO 17 Penman模型是基于能量平衡法并引入干燥力的概念后得到的僅需普通氣象資料可計算ET0的方法，屬綜合法范疇[12]。

(2)

式中：p0為海平面氣壓，hPa；p為本站地面氣壓，MPa；a1=0.26[13]；b1=0.14[13]。

(3)Hargreaves-Samani(HS)。Hargreaves-Samani模型是基于溫差來反映輻射項的ET0計算法，該模型被證明是缺失輻射資料情況下的估算ET0的有效方法之一，屬溫度法范疇[14]。

ET0=?(T+17.8)(Tmax-Tmin)0.5Ra

(3)

式中：?為參數(shù)，取0.002 3；Tmax為2 m高度處日最高氣溫，℃；Tmin為2 m高度處日最低氣溫，℃；Ra為大氣頂層輻射，MJ/(m2·d)，計算公式為[13]：

3.3.3 人類活動影響分析 人類通過實踐活動，不斷影響并改造著大氣圈、水圈、巖石圈及生物圈.人類活動于遙感影像上最直接的體現(xiàn)即為地表覆被的變化.通過對土壤侵蝕與土地利用類型的研究，進一步探討人類活動對土壤侵蝕的影響.基于Landsat影像，通過監(jiān)督分類的方法，獲取椒江流域各年份地表覆被情況椒江流域土壤侵蝕與地表覆被統(tǒng)計見表4(其他主要為居民地所產(chǎn)生土壤侵蝕的統(tǒng)計，微度侵蝕中包含水體統(tǒng)計).

Ra=37.6dr(wssinφsinδ+cosφcosδsinws)

(4)

dr=1+0.033 cos(0.017 2J)

(5)

ws=arccos (-tanφtanδ)

(6)

式中：dr為日地相對距離，無量綱；ws為日照時數(shù)角，rad；φ為地理緯度，rad；δ為太陽磁偏角，rad；J為日序數(shù)，例如1月1日為1，逐日累加。

(4)Makkink[15](Mak)模型，屬輻射法范疇。

(7)

式中：Rs為太陽輻射，MJ/(m2·d)，計算公式為[11]：

(8)

(5)Priestley-Taylor(PT)。Priestley-Taylor模型是假定在濕潤環(huán)境條件下，忽略了空氣動力學(xué)項并引進參數(shù)而得出的簡化算法，是Penman方程的簡化式，屬輻射法范疇[16]。

(9)

(6)Irmark-Allen。Imark-Allen(IA)模型是基于美國濕潤地區(qū)資料得到的經(jīng)驗式，計算時需輸入地理位置、氣溫、日照時數(shù)等參數(shù)。

ET0=0.489+0.289Rn+0.023T

(10)

1.3 數(shù)據(jù)分析方法

本研究依據(jù)典型性及全面性原則在各灌區(qū)范圍內(nèi)或附近選取一氣象站點為代表站[17]，以FAO 56Penman-Monteith公式計算的ET0值為標準，對其他5種模型進行誤差分析，其中誤差統(tǒng)計變量包括平均絕對誤差(MAE)、均方根誤差(RMSE)和納什效率系數(shù)(NSE)，并用Microsoft Office Excel 2010 計算統(tǒng)計變量并作圖；對氣象數(shù)據(jù)缺測值及異常值，本研究進行了合理插值，以期延長數(shù)據(jù)系列的長度。各統(tǒng)計變量計算公式為：

(11)

(12)

(13)

2 結(jié)果與分析

2.1 不同模型模擬的日ET0值比較

圖1～3為HS、FAO 17 PM、PT、IA和Mak模型模擬的墨達、江北及滿拉灌區(qū)的ET0日值變化。由圖1～3可知，不同模型模擬的ET0日值在年際間呈先增后減的變化規(guī)律，且ET0最大值、最小值分別出現(xiàn)于6-7月、12-次年1月。HS、FAO 17 PM模型模擬的ET0最大值、最小值均大于FAO 56 PM模型的計算結(jié)果；較HS、FAO 17 PM模型，PT模型對ET0最小值的模擬結(jié)果吻合度較高；Mak、IA模型模擬的ET0最大值、最小值較其他模型的模擬值更接近于FAO 56 PM模型計算結(jié)果，但IA模型對ET0最小值的模擬結(jié)果存在高估問題。

圖1～3直觀地反映了各模型模擬值的變化趨勢，但不足以定量評價其適用性。因此，本研究以FAO 56 PM計算值為標準，對其他5種模型進行相關(guān)性分析，結(jié)果見表2。由表2可知，F(xiàn)AO 56 PM計算值與其他5種模型模擬結(jié)果間線性關(guān)系極顯著。各灌區(qū)中FAO 17 PM、HS、PT模型的斜率a介于1.58～3.06之間，決定系數(shù)R2介于0.83～0.91之間，說明這3種模型的模擬值偏大，存在高估現(xiàn)象，其高估程度依次排序為：PT模型>HS模型>FAO 17 PM模型；Mak、IA模型的斜率a介于0.82～0.98之間，決定系數(shù)R2介于0.84～0.89之間，說明這2種模型的模擬值較其他模型更準確。

表3為FAO 56 PM模型與5種ET0模型模擬結(jié)果間的誤差統(tǒng)計結(jié)果。由表3可知，F(xiàn)AO 56 PM模型模擬的墨達、江北和滿拉灌區(qū)多年平均ET0日值介于2.72～2.84 mm/d之間；HS、FAO 17 PM、PT、IA和Mak模型的模擬值介于2.91～7.75 mm/d之間，明顯高于FAO 56 PM模擬結(jié)果。在各灌區(qū)中，F(xiàn)AO 17 PM模型的RMSE值介于1.83～2.03 mm/d之間，MAE值介于1.62～1.79 mm/d之間；HS模型的RMSE值介于4.80～5.30 mm/d之間，MAE值介于4.32～4.85 mm/d之間；Mak模型的RMSE值介于0.42～0.48 mm/d之間，MAE值介于0.35～0.39 mm/d之間；PT模型的RMSE值介于5.35～5.54 mm/d之間，MAE值介于4.69～4.91 mm/d之間；IA模型的RMSE值介于0.68～0.75 mm/d之間，MAE值介于0.62～0.68 mm/d之間。由各模型的MAE、RMSE值排序(PT模型>HS模型>FAO 17 PM模型>IA模型>Mak模型)可以看出，PT模型模擬的ET0日值距FAO 56 PM計算結(jié)果的離散度最大，計算精度最低；Mak模型的離散度最小，精度最高；其他模型的離散度及計算精度介于PT與Mak模型之間。分析NSE值發(fā)現(xiàn)，F(xiàn)AO 17 PM、HS、PT模型在3個灌區(qū)的NSE值均小于0，表明模型的計算可信度有待提高；Mak、IA模型的NSE值均大于0，表明模型的計算結(jié)果具有一定的可信度。其中，Mak模型的NSE均值(0.84)較IA(0.60)更接近于1，說明Mak模型的模擬質(zhì)量、計算可信度更高，IA模型次之。

圖1 江北灌區(qū)ET0日值變化

圖2 墨達灌區(qū)ET0日值變化

圖3 滿拉灌區(qū)ET0日值變化

綜上認為，各模型在模擬墨達、江北和滿拉灌區(qū)ET0日值方面，其適用性存在明顯差異。其中，Mak模型的模擬效果最優(yōu)，適用性最好；IA模型次之；FAO 17 PM、HS和PT模型較差，且較FAO 56 PM模擬結(jié)果，ET0日值存在不同程度的高估問題。

表2 FAO 56 PM與其他5種模型模擬的ET0日值間的回歸系數(shù)

注：表中“a”為相關(guān)分析擬合直線斜率；R2為確定性系數(shù)；“**”表示極顯著相關(guān)(p<0.01)。

表3 FAO 56 PM與其他5種模型模擬的ET0日值間的誤差分析

2.2 不同模型模擬的ET0年值比較

各ET0模型模擬的灌區(qū)ET0年值變化見圖4～6。由圖4～6可知，F(xiàn)AO 17 PM、HS、Mak、PT及IA模型模擬的江北、墨達和滿拉灌區(qū)ET0年值普遍高于FAO 56 PM模擬結(jié)果。對比分析發(fā)現(xiàn)，Mak模型模擬值最接近于FAO 56 PM的模擬結(jié)果，IA模型次之，這與日尺度的模擬結(jié)果一致。從ET0年值變化趨勢可以看出，F(xiàn)AO 17 PM、Mak、PT、IA模型模擬值的變化趨勢與FAO 56 PM模型的基本一致，但HS模型在江北、墨達灌區(qū)的模擬情況卻與之相反。此外，HS、PT模型模擬的ET0年值在年際間的浮動變化較其他模型大；FAO 17 PM和IA模型的變化較平穩(wěn)，但Mak模型的浮動變化與FAO 56 PM模型的相近。

誤差分析結(jié)果(表4)表明，江北、墨達和滿拉灌區(qū)中PT模型的MAE值介于1 692.45～1 786.01 mm/a之間，RMSE值介于1 693.49～1 786.65 mm/a之間；HS模型的MAE值介于1 535.71～1 768.79 mm/a之間，RMSE值介于1 536.88～1 733.80 mm/a之間；FAO 17 PM模型的MAE值介于568.28～647.08 mm/a之間，RMSE值介于575.05～647.51 mm/a之間；IA模型的MAE值介于194.07～214.69 mm/a之間，RMSE值介于195.32～218.00 mm/a之間；Mak模型的MAE值介于76.19～80.90 mm/a之間，RMSE值介于81.08～83.53 mm/a之間。對比分析發(fā)現(xiàn)，PT模型的MAE、RMSE值最大；Mak模型的MAE、RMSE值最低，這表明PT模型的模擬精度最低，Mak模型的最高，其他模型模擬精度介于二者之間。由各灌區(qū)不同模型的NSE值可知，F(xiàn)AO 17 PM、HS、PT模型的NSE值介于-3571.76～-118.00之間，其結(jié)果遠遠小于0，表明該模型的模擬結(jié)果不可信；Mak模型的NSE值較其他模型更接近于0，表明該模型模擬結(jié)果接近于FAO 56 PM的模擬值，即總體模擬結(jié)果可信，但模擬過程的誤差較大。

圖4 江北灌區(qū)ET0年值變化

圖5 墨達灌區(qū)ET0年值變化

圖6 滿拉灌區(qū)ET0年值變化

表4 FAO 56 PM與其他5種模型模擬的年ET0值間的誤差分析

3 討 論

國內(nèi)外學(xué)者研究并提出了多種基于溫度、輻射、風(fēng)速及綜合氣象參數(shù)的ET0計算模型[18,19]，但起初因模型建立的需要，模型的結(jié)構(gòu)中引進了不同的估算系數(shù)、空氣動力學(xué)項及輻射項參數(shù)，致使多數(shù)模型在實際應(yīng)用中受區(qū)域下墊面及氣候變化的影響，導(dǎo)致其普適性變差。因此，在實際應(yīng)用中需對ET0計算模型進行適用性評價，避免盲從引起作物需水量計算的精度降低。

本研究以FAO 56 PM模型模擬的西藏江北、墨達和滿拉灌區(qū)的ET0值為標準，分別采用5種參考作物蒸散量模型模擬ET0變化，并從不同時間尺度進行了分析。日尺度上，Mak、IA、FAO 17 PM、HS和PT模型的模擬值變化趨勢與FAO 56 PM模型的基本一致，在年際間呈先增后減的變化規(guī)律，且模擬值分別在6-7月、12月-次年1月出現(xiàn)峰值和谷值。但各模型的模擬效果存在明顯差異，其中Mak模型的模擬質(zhì)量及適用性最好，RMSE、MAE、NSE值分別為0.45 mm/d、0.37 mm/d、0.84；IA模型次之，RMSE、MAE、NSE分別為0.71 mm/d、0.65 mm/d、0.62；PT模型的最差，MAE值最大達4.91 mm/d，NSE值為-21.65，模擬結(jié)果的可信度低。年尺度上，各模型的模擬值較FAO 56 PM計算結(jié)果均存在高估現(xiàn)象，以PT模型最為突出。誤差分析結(jié)果(表4)表明，PT模型的MAE、RMSE值最大，Mak模型的MAE、RMSE值最低，這說明PT模型模擬的ET0年值的精度最低，Mak模型的模擬精度最高。此外，年NSE值進一步也證實了PT模型的模擬結(jié)果不可信，而Mak模型結(jié)果總體可信，但在年尺度的模擬過程中Mak模型的誤差較大。由徐冰等[8]在西藏當雄、那曲、改則等地研究認為的PT模型在日、年尺度上的模擬效果均最優(yōu)的結(jié)論可知，同一種ET0計算模型在西藏高原的適用性存在區(qū)域差異，這與海拔、氣象條件及區(qū)域下墊面的差異有關(guān)。

西藏高原平均海拔4 000 m以上，特殊的地域環(huán)境決定了西藏高原區(qū)氣候具有低氣溫、強輻射的區(qū)域特征。因此，完全基于溫度參數(shù)的HS模型難以反映出藏區(qū)高原實際強輻射條件下的ET0變化。此外，HS模型也未考慮空氣動力學(xué)項影響[20]，這是引起該模型計算結(jié)果無論是在日尺度還是在年尺度中均較FAO 56 PM模型計算結(jié)果存在較大偏差的主要原因。而Mak、IA模型則是以輻射參數(shù)為基礎(chǔ)的估算方法，在本研究中較其他3種模型在日尺度和年尺度上均表現(xiàn)出較好的適用性，且MAE、RMSE及NSE分析也證實了Mak、IA模型模擬結(jié)果的準確性及可靠性。然而，IA模型結(jié)構(gòu)中由于僅考慮了平均溫度的影響，未考慮水汽壓差、凈輻射、總輻射等關(guān)鍵氣象因子的影響[21,22]，致使其模擬精度及可靠性較Mak模型低。PT模型的模擬結(jié)果較FAO 56 PM存在較大偏差，其主要原因與PT模型計算時運用了不同的修正參數(shù)計算的空氣動力學(xué)項和未考慮相對濕度變化的影響有關(guān)[16]，即PT模型是假定在濕潤條件下并忽略了空氣動力學(xué)項影響后僅引進常數(shù)α=1.26進行代替而得出的簡化方法，并未考慮不同區(qū)域特征對α值的影響[23]；已有的研究結(jié)果表明ET0與相對濕度間呈負相[24]，即當未考慮相對濕度對ET0的影響時，將引起PT模型計算結(jié)果偏大，這在本研究中已得到證實。雖然，F(xiàn)AO 17 PM模型與FAO 56 PM同屬于綜合法，也考慮了空氣動力學(xué)項及輻射項的影響，但2種模型因采用了不同的空氣動力學(xué)和輻射項計算公式，致使其結(jié)果間存在較大偏差。而造成FAO 17 PM模型日、年尺度模擬值均大于FAO 56 PM模型計算結(jié)果的主要原因與FAO 17 PM模型的結(jié)構(gòu)中未考慮土壤熱通量的影響[25]有關(guān)。然而，綜合法因其模型結(jié)構(gòu)中需估算的中間參數(shù)較多且計算過程繁雜，同時細節(jié)誤差也可能會影響計算結(jié)果，使得目前關(guān)于綜合法的評價研究仍存在不確定性[18]。

4 結(jié) 論

本文根據(jù)西藏高原江北、墨達和滿拉灌區(qū)的長系列氣象資料，以FAO 56 PM模型計算值為參照標準，對Mak、IA、FAO 17 PM、HS和PT模型進行分析驗證后認為：

(1)Mak、IA、FAO 17 PM、HS和PT模型模擬的西藏高原江北、墨達和滿拉灌區(qū)的日尺度ET0變化趨勢與FAO 56 PM模型的一致，均在年際間呈先增后減的變化規(guī)律，且各模型模擬值分別在6-7月、12月-次年1月出現(xiàn)峰值和谷值。

(2)線性回歸及誤差分析表明，Mak、IA、FAO 17 PM、HS和PT模型在江北、墨達和滿拉灌區(qū)的適用性存在明顯差異。其中，Mak模型的日尺度模擬精度及可信度最高，MAE、RMSE、NSE值分別為0.37 mm/d、0.45 mm/d和0.84；IA模型次之，RMSE、MAE、NSE分別為0.71 mm/d、0.65 mm/d、0.62；PT模型最差，MAE值最大達4.91 mm/d且NSE值小于0，其模擬結(jié)果可信度低。年尺度中，各模型較FAO 56 PM均存在高估現(xiàn)象，以PT模型最為突出。但相對而言，PT模型模擬精度最低，模擬質(zhì)量不可信；Mak模型精度最高，模擬質(zhì)量總體可信，但模擬過程誤差較大。

(3)Mak模型結(jié)構(gòu)簡單，所需氣象資料少，且模擬結(jié)果無論在日尺度還是年尺度中均最接近于FAO 56 PM模型的模擬結(jié)果。綜合評定，可以考慮將Mak模型作為代替FAO 56 PM模型來計算西藏高原灌區(qū)氣象數(shù)據(jù)缺測下ET0值。