考慮“蒸發悖論”的洱海灌區逐日參考作物蒸散發預測

趙 眾，周 密，張劉東，顧世祥,2*，李 靖

（1.云南農業大學 城鄉水安全與節水減排高校重點實驗室，昆明 650201；2.云南省水利水電勘測設計研究院，昆明 650021；3.武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072）

0 引 言

【研究意義】蒸散發（ET）作為表征太陽輻射到地面的水汽和能量轉化的重要參數，是區域水循環中最重要的水文過程之一[1]。地面氣象觀測和衛星遙感等多途徑資料分析顯示，由于植被葉面積指數增加，全球自1978年以來蒸散發平均增長值為0.63 mm/a，其中植被蒸騰增加速度為0.72 mm/a，土壤蒸發則以0.32 mm/a 的速度減少[1-2]。參考作物蒸散發（ET0）又是蒸散發的氣象因素分量，其估算和準確測定具有重要的現實意義[3-5]。【研究進展】逐日ET0的計算模型主要包括Penman-Monteith 公式[5]、Priestley-Taylor模型[6]、Hargreaves 公式[7]、Irmark-Allen 法[8]等方法。Penman-Monteith 公式作為聯合國糧農組織（FAO）推薦的基本公式，其理論完備性、計算精度和適用性已經得到全球的公認，但所需的資料信息過多，故而在數據缺乏的地區難以使用。ET0的變化不是由于單一氣象因素的變化而決定的，Traore 等[9]采用氣象因素基于人工智能網絡對ET0進行研究，發現Tmax是很重要的因素，但其對于預測的準確性還取決于太陽凈輻射的精度和實時的天氣預報信息。

除上述的理論及經驗模型計算法之外，FAO 還建議使用蒸發皿資料來確定ET0（即直接法），通過設計試驗環境，采用蒸發皿等觀測儀器觀測水面蒸發量，再經過折算系數Kp轉換計算得到ET0[10-11]。蒸發皿蒸發作為揭示湖泊流域區氣候變化響應最直接的指標，在我國鄱陽湖、巴丹吉林沙漠湖泊、洱海等不同地理區域的增減變化等研究均有應用[12-15]。在工程實踐中過去單純使用氣溫、日照等單一或多個氣象因子資料，運用人工智能或統計模型的方法得到ET0估計值，近年來已逐漸將蒸發皿蒸發作為新增輸入項以提高ET0估計精度[16-17]。Lennartz 等[18]將A 級蒸發皿資料用于模擬計算逐日、逐小時等不同時間尺度的ET0，發現在逐日的時間尺度下模擬值與標準值的決定系數高達0.90～0.94，但在小時等短時間尺度由于微氣象環境的隨機波動，模擬效果較差。也有采用輻射與溫度等參數法估算ET0[19-20]，效果也很好。張鑫等[21]針對西南地區分析了9 種Kp經驗模型的適用性，得出A98 模型的適用性相對較好。【切入點】折算系數Kp的確定需全面考慮蒸發皿所處的下墊面、局部環境及風場和濕度等條件[5]，數值范圍不適于大尺度的推廣使用，在計算和推廣上帶來一定困難。由于“蒸發悖論”現象的存在，且“蒸發悖論”的發生與時間尺度的長短又有關[22]，因此在采用直接觀測法對其進行分析計算時，若ET0與水面蒸發的變化趨勢相反時，也會導致其結果走向出現誤差。【擬解決的關鍵問題】探究出現“蒸發悖論”時適宜的參考作物蒸散發量估算模型。

1 材料與方法

1.1 數據資料

洱海位于云南省大理白族自治州，屬瀾滄江—湄公河水系，流域面積2 785 km2，年均氣溫15.1 ℃，年降水量1 057 mm，水資源總量10.7 億m3[23]。全球氣候變化及低緯度高原的區域響應，2010—2015年洱海地區遭遇了的連續干旱災害，流域農業用水及農田面源加大、湖水位下降、局部區域藍藻爆發等問題突出，洱海流域水生態保護治理已成為全國關注焦點之一[23]。使用資料包括：①洱海灌區內代表性的大理氣象站1954—2018年逐日平均氣溫（T）、最低氣溫（Tmin）、最高氣溫（Tmax）、日照時間（n）、風速（u）、相對濕度(RH)、降水量、20 cm 蒸發皿蒸發量（Epan）。1954—2001年的20 cm 蒸發皿蒸發量為大理站實測數據，2002—2018年缺乏實測值，為保證不同年代蒸發資料的一致性，采用“云南省地表水資源”書中E601/E20=0.69 等直線圖，按0.69 采用E601的實測資料對20 cm蒸發皿蒸發量進行折算。②大理、洱源、賓川等市縣的自然地理、社會經濟、農業綜合統計年報、水利（或年鑒）等技術成果。

1.2 逐日ET0 模擬預測方法

加入蒸發皿蒸發資料后，逐日參考作物蒸散發ET0預測模擬的精度得到顯著提高，本質上是基于中長期蒸發皿蒸發與ET0標準值的高度相關性。但具體到不同的時期和季節，若出現“蒸發悖論”現象時，二者的變化趨勢完全相反，就不能再繼續使用蒸發皿蒸發來預測模擬逐日ET0了。為此，本文先采用M-K非線性趨勢檢驗法對年、春夏秋冬季節的蒸發皿蒸發和ET0的變化進行趨勢檢驗，識別出“蒸發悖論”時段。其次，判斷模擬預測時段所處時期是否存在“蒸發悖論”現象，以決策采用不同的ET0實時預測模型方法。第三，在“蒸發悖論”時期，由于蒸發皿蒸發與ET0的變化趨勢相反，只能采用隨機方法之Copula聯合分布函數模型預測[24]。第四，未發生“蒸發悖論”時期，蒸發皿蒸發與ET0的變化趨勢一致，可采用多元線性回歸模型和Kp折算系數模型，或進一步加入蒸發皿蒸發項以提高逐日ET0預測模擬的精度，其技術流程圖見圖1。為分析不同方法模型預報ET0的精度，使用4 個常用的統計指數來進行評估，分別為平均偏差（ME）、均方根誤差（RMSE）、符合指數（IA）、Nash-Sutcliffe 效率系數（NSE）[25]。符合指數（IA）在0 和1 之間，IA越大模擬效果越好。Nash-Sutcliffe效率系數（NSE）變化范圍從-∞到1，值越接近1 說明模擬值和實際值越接近。

圖1 技術流程Fig.1 Research technical flow chart

1.2.1 M-K 趨勢檢驗分析

為檢驗ET0、蒸發皿蒸發及氣溫的變化趨勢，采用Mann Kendall（M-K）非線性趨勢檢驗法[26]。具體檢驗原理如下：

設時間序列為{xi}（i=1,2，…，n），{xj}的對偶數S（xi

式中：sgn（）為符號函數；U為M-K 統計量。使用M-K 法進行突變檢驗時，其原理如下：

利用{xi}（i=1,2，…，n）構造一新序列，即：

式中：mi為xi>xj(1≤j≤i)的樣本累積數；dk的均值以及方差定義如下：

在時間序列隨機獨立假設下，定義統計量：

所有的UFk組成1 條UF曲線，把同樣的方法引用到反序列中，得到1 條UB曲線，將統計量UF和UB2 條曲線與±1.96 2 條直圖線繪制到同一坐標表下，當出現UF和UB曲線2 條線的交點時，即為突變點。

1.2.2 多元線性回歸模型

將數據輸入SPSS 統計分析工具中，采用逐步回歸方法進行變量引入，采用F檢驗概率作為判斷標準值，進入概率＜0.05，移除概率＞0.1。

1.2.3 多元Copula 函數預測ET0

構建T、Tmax、Epan3 種氣象因素的邊緣分布函數，選用有代表性的正態分布、gamma 分布、lognormal分布、weibull 分布分別構建T、Tmax、Epan的邊緣分布函數，并從中選擇擬合效果最好的分布函數。T選擇weibull 分布，Tmax選擇weibull 分布，Epan選擇正態分布。之后分別構建T-Tmax-Epan三維Normal Copula模型，T-Tmax、T-Epan、Tmax-Epan的3 種二維Normal Copula 模型。計算出3 種氣象因素之間不同組合對應的聯合分布，并將其記為預測ET0分布概率，在帶回ET0的邊緣分布函數中計算出ET0的值，其中ET0的邊緣分布函數為gamma 函數，即為預測ET0值，優選出預測效果最佳的氣象因子組合模型。詳細的方法模型及應用參見相關文獻[24]。

2 結果與分析

2.1 變化趨勢分析

采用M-K 趨勢檢驗法對大理站1954—2018年共65 a 的蒸散發量、氣溫及20 cm 蒸發皿蒸發量長時間序列變化趨勢，以及各年代際間的變化趨勢進行檢驗，結果如表1。并繪制年值及春夏秋冬季各自65 a 長時間序列的所示變化趨勢線以及5 a 滑動平均曲線（圖2）。

表1 大理站1954—2018年ET0、T 和Epan 的M-K 統計值Table Dali meteorological station’s M-K statistics of ET0 in 1954—2018

由表1 可知，大理站參考作物蒸散發ET0在1954—2018年長時間序列中M-K 統計值為0.07，呈上升趨勢，但除1990年和2000年呈上升趨勢外，其余年代呈下降趨勢，由圖2 可看出，整體呈上升趨勢。氣溫的M-K 統計值為3.13，呈上升趨勢，且其變化通過0.05 的顯著性檢驗，1960、1970年及2010年呈下降趨勢。由圖2 可看出，整體呈先下降后上升趨勢。20 cm 蒸發皿蒸發量的M-K 統計值為-5.46，呈下降趨勢，且其變化通過0.05 的顯著性檢驗，只有1990年呈上升趨勢，并且并由圖1 可看出，整體呈下降趨勢。

春夏秋冬四季按春季2—4月、夏季5—7月、秋季8—10月、冬季11月—翌年1月分段計。對于春季，蒸散發量和氣溫基本上呈上升趨勢，而20 cm 蒸發皿蒸發量呈下降趨勢；參考作物蒸散發及氣溫均只有1960年和2010年呈下降趨勢，20 cm 蒸發皿蒸發量在20 世紀各年代（10 a）變化不大，進入22 世紀后呈下降趨勢。對于夏、秋、冬三季，蒸散發量和20 cm 蒸發皿蒸發量基本呈下降趨勢，而氣溫大多呈上升趨勢。

圖2 1954—2018年ET0、氣溫、20 cm 蒸發皿蒸發量的變化趨勢Fig.2 1954—2018 trend chart of ET0、temperature and evaporation capacity of 20 cm evaporating dish

2.2 突變分析

采用M-K 趨勢檢驗法對大理站蒸散發量、氣溫以及20 cm 蒸發皿蒸發量1954—2018年共65 a 的長時間序列的變化趨勢以及各季節的變化趨勢進行突變分析，其結果見圖3。由圖3 可知，蒸散發量對于全年的長時間序列來說，UF和UB的交點較多，表示其變化較為頻繁，而從M-K 統計值分析可知在1960、1970年及2000年出現“蒸發悖論”，但除1960年出現交點較多外，1970年和2000年均未出現交點。對于春、秋2 季，春季突變在2008年和2016年，秋季突變在1993年和2015年。夏、冬2 季則變化較為頻繁。20 cm 蒸發皿蒸發量對于全年的長時間序列來說，UF和UB的交點出現在1995年，即表示突變出現在1995年，處于1990年，但只有夏季出現“蒸發悖論”；春、冬2 季突變均出現在1996年，夏季出現在1994年，秋季出現在1987年。氣溫對于全年的長時間序列來說，UF和UB的交點出現在2009年，表示其突變出現在2009年，處于2000年。其春、冬2 季出現在2008年，夏季出現在2009年，秋季出現在2007年。

蒸散發量、20 cm 蒸發皿蒸發量和氣溫出現突變的時間不同步，這可能是由于影響蒸散發量和20 cm蒸發皿蒸發量的因素不僅有氣溫，還與日照、輻射、降水等多種因素有關。

圖3 1954—2018年ET0、氣溫、20 cm 蒸發皿蒸發量的變化趨勢Fig.3 1954—2018 mutation analysis of ET0、temperature and evaporation capacity of 20 cm evaporating dish

2.3 逐日ET0 預測方法對比分析

2.3.1 未出現“蒸發悖論”時期

1）單純的氣象因子多元線性回歸模型法（ET0,Metr）

未引入蒸發皿蒸發項時，為方便對比計算，選擇p=5%、25%、50%、75%、97%（第1 組典型年）所對應的1993、1987、1996、1986、2012年共5 a 的數據帶入模型之中進行驗證。構建氣象因素與ET0的線性回歸模型，其氣象因素對ET0影響的顯著性及其模型表達式見表2。表3 為不同水文年進行預測ET0，其結果與實際逐日平均ET0（作為標準值）進行對比的結果，也將Copula 聯合分布函數模型法的結果一并列出（記為ET0,Copl）。由表3 可知，2 種模型的計算結果，其相對誤差小于10%、15%、20%、25%的樣本數比例基本上都是多元線性回歸模型優于Copula 聯合分布函數模型。

表2 各氣象因素在不同年份對ET0 影響顯著性結果（第1 組典型年）Table 2 Significant results of influence of meteorological factors on ET0 in different years（Typical year of the first group）

表3 ET0 預測值與實際值精度統計（第1 組典型年）Table 3 Statistical table of accuracy of ET0 predicted value and ET0 actual value（Typical year of the first group）

2）蒸發皿折算系數Kp模型法（ET0,Kp）

有學者研究蒸發皿系數Kp=ET0/Epan，將實際的Kp作為因變量，以實測的逐日相對濕度和2 m 高度處風速作為自變量。本文根據前人研究[21]，選用A98模型，并根據大理站的氣象因素將其模型修正后為下：

式中：Kp即為蒸發皿系數；U2為2 m 高度處風速；RH為相對濕度；F為生長作物頂風吹程，20 cm 蒸發皿安置與地面距離顯著大于Class-A 型蒸發皿，受F影響較小，F取經驗值20 m。預測結果與標準值見表4。其預測結果總體上不如采用多元線性回歸模型對ET0進行預測的效果好，這可能是由于該模型在預測時間尺度較小時不適用，為進一步探究原因，又用該模型進行了月尺度的ET0預測，其結果顯示，相對誤差＜10%的樣本數增加到45%、＜15%的樣本數為55%、＜20%的樣本數為76.7%、＜25%的樣本數為90%、＜30%的樣本數為95%。由此可知，該模型在月尺度上的ET0預測精度較高，但在日尺度上預測精度較低，可能是由于大理站風速在中短時間尺度下變化各異，高原盆地和湖泊區形成的小氣候影響所。

表4 Kp 蒸發皿折算系數模型ET0 預測值與實際值精度統計（第1 組典型年）Table 4 Statistical table of accuracy of ET0 predicted value and ET0 actual value of conversion coefficient model of Kp evaporating dish（Typical year of the first group）

表5 各氣象因素在不同年份對ET0 影響顯著性結果（引入Epan）Table 5 Significant results of influence of meteorological factors on ET0 in different years（introduce Epan）

表6 ET0 預測值與實際值精度統計（引入Epan）Table 6 Statistical table of accuracy of predicted value and actual value（introduce Epan）

3）加入Epan后的多元回歸模型法（ET0,Epan+Metr）

從大理站氣象因子、蒸發皿蒸發與逐日ET0的線性回歸分析看，將20 cm 蒸發皿蒸發量與氣象因子一并作為輸入，采用多元線性規劃方法推求模型參數、預測逐日ET0，并進行誤差分析，結果如表5、表6所示。對比表3 未引入Epan時的結果，其預測值與實測值的精度均有提升。表7 為單純的氣象因子多元線性回歸模型法（ET0,Metr）、T-Tmax二維Copula 聯合分布函數模型法（ET0,Copl）、蒸發皿折算系數Kp模型法（ET0,Kp）和加入Epan后的多元回歸模型法（ET0,Epan+Metr）等4 種模型方法，進行前述典型年組的逐日ET0預測，相對誤差ERR小于10%、15%、20%和25%的樣本量占比。由表7 可知，在未出現“蒸發悖論”的時期，加入Epan后的多元回歸模型法（ET0,Epan+Metr）所得預測ET0的預測結果與標準值間的誤差最小，其次為單純的氣象因子多元線性回歸模型法（ET0,Metr），最差的為蒸發皿折算系數Kp模型法（ET0,Kp）。另一方面，參與對比的T-Tmax二維Copula聯合分布函數模型（ET0,Copl）的預測效果總體最佳，進一步表明其方法的普適性。

2.3.2 出現“蒸發悖論”時期

為方便計算，選擇p=7%、37%、47%、77%、95%（第2 組典型年，處于“蒸發悖論”期）所對應的2000、2016、2001、2003、1960年共5 a 的數據帶入模型之中進行驗證。構建氣象因素與ET0的線性回歸模型，其氣象因素對ET0影響的顯著性及模型見表8。由表9 可看出，在“蒸發悖論”時期，蒸發皿蒸發與同期ET0的變化趨勢相反，只有采用修正后的T-Tmax二維Copula 模型（記為ET0,Copl），其預測精度也顯著高于多元線性回歸模型。

表7 4 種逐日ET0 預測模型的相對誤差占比Table 7 Relative error ratio of four daily ET0 prediction models

表8 各氣象因素在不同年份對ET0 影響顯著性結果（第2 組典型年）Table 8 Significant results of influence of meteorological factors on ET0 in different years（Typical year of the second group）

表9 逐日ET0 預測值與實際值精度統計(第2 組典型年)Table 9 Statistical table of accuracy of ET0 predicted value and ET0 actual value

3 討 論

“蒸發悖論”的結果與謝平[22]等采用云南省1981—2011年長時間序列共52 個站點的氣象數據對云南省內滇西地區“蒸發悖論”的探討結果趨勢相同，“蒸發悖論”的有無與時間段的長短選取有關，長時段的研究序列可能會掩蓋短時段的“蒸發悖論”現象。

蒸散發量、20 cm 蒸發皿蒸發量和氣溫出現突變的時間不同步，這可能是由于影響蒸散發量和20 cm蒸發皿蒸發量的因素不僅有氣溫，還與日照、輻射、降雨等多種因素有關[27]。孫潔等[28]通過對鄂爾多斯的蒸散發研究，發現風速的下降是潛在蒸散發量減小的主要因素；同時，日照時間的減小和降水量的增加也是鄂爾多斯高原西部潛在蒸散發量減小的關鍵因素。與分析結果相一致。

對于參考作物蒸散發量ET0的多元線性回歸模型預測方法、數值計算方法都有較多研究[5-11]，本文將“蒸發悖論”與幾種計算方法相結合，找出在未出現“蒸發悖論”時采用線性回歸模型方法較優，而在出現“蒸發悖論”時，則采用高維Normal Copula 函數模型進行預測較優。

4 結 論

在1954—2018年的長時間序列上，洱海流域大理站的20 cm 蒸發皿蒸發量均呈下降趨勢，參考作物蒸散發量ET0和氣溫呈上升趨勢，ET0的上升趨勢更平緩，但在短時間序列上結果存在年代各異性。

未出現“蒸發悖論”時期，加入Epan后的多元回歸模型法（ET0,Epan+Metr）所得逐日ET0的預測結果與標準值間的誤差最小，其次為單純的氣象因子多元線性回歸模型法（ET0,Metr），最差為蒸發皿折算系數Kp模型法（ET0,Kp）。因此，可采用加入Epan后的多元回歸模型法（ET0,Epan+Metr）進行逐日ET0預測，氣溫對其影響的顯著性最大。出現在“蒸發悖論”時期，采用T-Tmax二維Normal Copula 模型的精度更高。

采用加入Epan的多氣象因子線性回歸模型進行預測時，其預測結果精度較高。而采用二維Normal Copula 模型預測所需氣象因素卻是最少，可在氣象因素缺失時獲得較好的預測結果。