基于Budyko理論的黑河流域中游地區(qū)實際蒸散發(fā)估算及其變化歸因分析

余加男，李占玲，馮雅茹

（1.中國地質(zhì)大學(xué)（北京）水資源與環(huán)境學(xué)院，北京100083；2.水利部牧區(qū)水利科學(xué)研究所，呼和浩特010020）

0 引言

近年來，隨著全球氣候變化和人類活動的影響，水文循環(huán)發(fā)生顯著變化。實際蒸散發(fā)作為水文循環(huán)中的一環(huán)，是水量平衡中非常重要的組分，也是直接受氣候變化和下墊面等因素綜合作用的水文要素，其大小受降水條件、蒸發(fā)能力等多種因素的共同制約。實際蒸散發(fā)可以反映不同流域水量和能量的變化特征，一直以來都是水文界研究的熱點問題。

當(dāng)前，常用的實際蒸散發(fā)估算方法有直接觀測法、遙感反演法、水文模擬法等[1-4]。直接觀測法是根據(jù)蒸發(fā)皿、蒸散儀等儀器測定實際蒸散發(fā)量，這種方法比較簡單，但由于該方法無法直接測定大尺度上的實際蒸散發(fā)，因此在應(yīng)用中具有一定的局限性。遙感反演法是利用遙感資料結(jié)合地面資料計算流域?qū)嶋H蒸散發(fā)，其優(yōu)點在于時空分辨率高、快速且適用于大面積長期觀測，近年來該方法得到了廣泛應(yīng)用，但遙感影像容易受到氣象因素的影響，在一定程度上會影響到估算結(jié)果的精度[4-7]。水文模擬法是利用水文氣象資料，建立概念性或基于物理機制的水文模型等，以此求得實際蒸散發(fā)量；然而，水文模型結(jié)構(gòu)往往比較復(fù)雜且需要較多的驅(qū)動數(shù)據(jù)，并需要率定多個模型參數(shù)才能達(dá)到較好的擬合效果[8-9]。相對而言，基于Budyko 理論的水熱耦合模型，由于其模型結(jié)構(gòu)簡單，模型參數(shù)少，目前已成為流域或區(qū)域?qū)嶋H蒸散發(fā)估算的有效方法之一[10-12]。

很多學(xué)者將基于Budyko 理論的一系列水熱耦合模型應(yīng)用于特定流域或區(qū)域的實際蒸散發(fā)估算。李鴻雁等從Budyko 水熱耦合平衡理論入手，采用傅抱璞模型估算了嫩江流域及典型匯水區(qū)的實際蒸散發(fā)，結(jié)果表明，傅抱璞模型估算的逐年蒸散發(fā)具有較高精度[11]；張靜等考慮土壤水蓄變量因子，改進(jìn)了傅抱璞模型，驗證了該模型在松花江流域?qū)嶋H蒸散發(fā)估算中的適應(yīng)性[12]；趙勇等基于Choudhury-Yang 模型計算了黃土高原涇河流域2000-2014年蒸散發(fā)變化，并定量區(qū)分了氣候因子和下墊面因子對流域蒸散發(fā)變化的貢獻(xiàn)[13]。這些研究極大豐富了水熱耦合模型在實際蒸散發(fā)估算中的應(yīng)用。然而，這些研究都是將降水量作為蒸散發(fā)的水分供給條件，這對于農(nóng)業(yè)綠洲而言，顯然是不合理的?；诖耍n松俊等在傅抱璞模型的基礎(chǔ)上，將降水量和上游來水量之和作為蒸散發(fā)的水分供給，提出了一種適用于綠洲的水熱耦合模型，并以塔里木盆地為例，驗證了其水熱耦合平衡關(guān)系，效果良好[14]。

黑河流域是我國西北地區(qū)第二大內(nèi)陸河流域。中游地區(qū)擁有大面積農(nóng)業(yè)綠洲，其中張掖綠洲擁有國家現(xiàn)代農(nóng)業(yè)試驗示范區(qū)，是我國西北地區(qū)重要的糧食和瓜果蔬菜生產(chǎn)基地。該區(qū)域農(nóng)業(yè)綠洲主要以黑河地表水灌溉為主，蒸散發(fā)是綠洲水分消耗的主要途徑。以往針對黑河流域開展的蒸散發(fā)研究多集中于整個流域[3,15,16]；然而，閆宇會等[7]研究表明，黑河流域多年平均實際蒸散發(fā)量具有明顯的空間差異性。因此，將蒸散發(fā)過程強烈的中游綠洲地區(qū)作為研究區(qū)，探討其實際蒸散發(fā)的變化規(guī)律并對不同影響因素的貢獻(xiàn)度進(jìn)行定量分析，對于科學(xué)指導(dǎo)農(nóng)業(yè)灌溉、優(yōu)化區(qū)域水資源配置與管理、緩解水資源短缺等方面均具有重要的指導(dǎo)意義。

鑒于此，本文采用適用于綠洲的改進(jìn)傅抱璞經(jīng)驗?zāi)Ｐ蛯诤恿饔蛑杏蔚貐^(qū)的實際蒸散發(fā)量進(jìn)行估算，采用Mann-Kendall 方法對實際蒸散發(fā)序列隨時間變化的趨勢特征和變點特征進(jìn)行分析，并通過彈性系數(shù)法對影響實際蒸散發(fā)變化的不同因素的貢獻(xiàn)度進(jìn)行評價，以期對研究區(qū)灌溉用水管理等提供依據(jù)。

1 研究區(qū)概況及數(shù)據(jù)資料

黑河是我國西北地區(qū)第二大內(nèi)陸河，位于河西走廊中部。黑河流域中游地區(qū)位于鶯落峽水文站到正義峽水文站之間（99°19′~100°35′E，38°58′~39°52′N），跨張掖市轄區(qū)內(nèi)的甘州區(qū)、臨澤縣和高臺縣三區(qū)縣，面積約2.56 萬km2，河道全長204 km。中游地區(qū)綠洲面積約占總面積50%以上，其中張掖綠洲東依武威、金昌，西接酒泉、嘉峪關(guān)，南與青海省毗鄰，北與內(nèi)蒙古接壤，位于巴丹吉林沙漠和騰格里沙漠邊緣，集中了全流域90%以上的耕地和人口，多年平均氣溫7 ℃左右，年內(nèi)降水分布不均，夏多冬少，東多西少，蒸散發(fā)與之相反[17]。

文中采用的氣象數(shù)據(jù)包括黑河流域中游高臺和張掖氣象站2006-2015年逐日降水和氣溫數(shù)據(jù)，鶯落峽和正義峽水文站同期逐日流量數(shù)據(jù)。研究區(qū)及站點位置示意圖如圖1所示。

圖1 黑河流域中游及站點位置示意圖Fig.1 Locations of meteorological and hydrological stations and the Heihe River Basin

2 研究方法

2.1 基于Budyko理論的改進(jìn)傅抱璞經(jīng)驗?zāi)Ｐ?/h3>

著名氣候?qū)W家Budyko 在年或多年時間尺度上，以降水量表示蒸散發(fā)的水分供給條件，潛在蒸散發(fā)表示能量供給條件，對陸面蒸散發(fā)限定了極端干旱、極端濕潤兩個邊界條件，提出了水熱耦合平衡方程式[18]。對于黑河流域中游地區(qū)來說，除降水量以外，還有上游來水量，因此以上游來水量與降水量之和作為總供水量[14]，則研究區(qū)水熱耦合平衡方程式可以寫為：

改進(jìn)的傅抱璞經(jīng)驗?zāi)Ｐ涂杀硎緸椋?/p>

式中：P為研究區(qū)指定時段內(nèi)降水量，mm；R1為研究區(qū)同時段來水量，mm；Ea表示實際蒸散發(fā)量，mm；E0表示潛在蒸散發(fā)量，mm；ω為模型參數(shù)，反映下墊面條件對蒸散發(fā)的作用，取值范圍為（1，+∞）；本文采用Mc Cloud 方法對研究區(qū)潛在蒸散發(fā)量進(jìn)行估算[19]。

在多年平均時間尺度上，流域水量均衡方程可以表示為：

式中：R2為研究區(qū)多年平均流出水量，mm。

結(jié)合研究區(qū)多年平均降水?dāng)?shù)據(jù)和徑流數(shù)據(jù)，計算出研究區(qū)多年平均實際蒸散發(fā)量，以該結(jié)果作為基準(zhǔn)，代入改進(jìn)的傅抱璞模型，根據(jù)試算法求出模型中的下墊面參數(shù)ω。

2.2 趨勢及變點分析方法

Mann-Kendall（M-K）檢驗法由Mann 和Kendall 提出，常用于水文、氣候等變量的趨勢特征分析[20-22]；其優(yōu)點在于樣本不需要遵從一定的分布，也不受少數(shù)異常值的干擾，計算比較簡便，是一種有效的非參數(shù)統(tǒng)計檢驗方法。此方法不僅可以檢驗時間序列的變化趨勢，也可以檢驗時間序列是否發(fā)生突變，方法介紹詳見參考文獻(xiàn)[23]。

2.3 歸因分析方法

彈性系數(shù)法也稱為全微分法，是對徑流變化進(jìn)行歸因分析時常采用的方法。對實際蒸散發(fā)變化的歸因分析也可以借助該方法實現(xiàn)[23]。本文中，蒸散發(fā)的水分供給條件用來水量和降水量共同表示，參數(shù)ω為定值，故假設(shè)E0、P和R1是自變量，實際蒸散發(fā)Ea的全微分可用于評價不同因素的貢獻(xiàn)，具體計算公式可表示為：

忽略式中泰勒展開的高階方程，上式可改寫為：

其中3個偏導(dǎo)數(shù)的計算公式為：

定義實際蒸散發(fā)對不同因素變化的彈性系數(shù)為：

彈性系數(shù)越大，表示該因素變化對實際蒸散發(fā)變化的影響越大，反之亦然。

根據(jù)變點分析結(jié)果，將實際蒸散發(fā)序列分為2個時段，分別為時段I 和時段II，定義不同因素對實際蒸散發(fā)變化的貢獻(xiàn)率為：

式中：Δx為各因素兩個時段的變化量；EaI為時段I 的實際蒸散發(fā)量。

3 結(jié)果與分析

3.1 模型參數(shù)的確定

研究區(qū)水分輸入包括降水和上游來水量，水分輸出包括蒸散發(fā)、流域出口流出水量，在多年平均時間尺度上，研究區(qū)水量儲存量的變化可忽略不計。通過計算，研究區(qū)多年平均年降水量為129 mm，多年平均來水量197 mm，多年平均流出水量33 mm，由水量均衡方程計算得到研究區(qū)多年平均實際蒸散發(fā)為293 mm；由Mc Cloud 方法計算得到多年平均潛在蒸散發(fā)量為552 mm。當(dāng)給定模型參數(shù)ω某一具體數(shù)值時，即可得到由改進(jìn)的傅抱璞模型估算的研究區(qū)實際蒸散發(fā)量。通過不斷調(diào)整ω取值，使得由水量均衡方程和改進(jìn)傅抱璞模型得到的研究區(qū)實際蒸散發(fā)估計值的相對誤差不斷減小，當(dāng)該誤差達(dá)到最小且在允許范圍之內(nèi)時，即可認(rèn)為此時ω值為模型參數(shù)的最優(yōu)值。圖2給出了由兩種方法計算的研究區(qū)實際蒸散發(fā)相對誤差隨ω取值的變化曲線?？梢钥闯?，當(dāng)模型參數(shù)ω=3.1時，相對誤差達(dá)到最小，不足0.1%，因此模型參數(shù)ω的最優(yōu)值為3.1。將確定后的ω值代入改進(jìn)的傅抱璞模型，計算出研究區(qū)逐年實際蒸散發(fā)量。

圖2 研究區(qū)多年平均實際蒸散發(fā)估算的相對誤差隨ω取值的變化曲線Fig.2 Variations of relative error of multi-year average actual evaporation with the values of parameter ω

3.2 實際蒸散發(fā)估算及特征分析

根據(jù)公式（2）計算得到2006-2015年研究區(qū)逐年實際蒸散發(fā)量，如圖3所示，其中，2007年實際蒸散發(fā)量最大，達(dá)342 mm，2011年實際蒸散發(fā)量最小，為255 mm，多年平均實際蒸散發(fā)量293 mm。此前，馬心依研究發(fā)現(xiàn)，黑河流域中游林地蒸散量最大，達(dá)到342.4 mm，低覆被草地蒸散發(fā)量最小，為211.2 mm[24]；閆宇會等根據(jù)MOD16 產(chǎn)品估算得到黑河流域多年實際蒸散發(fā)的波動范圍為222.1~352.1 mm，多年平均實際蒸散發(fā)量值為292.0 mm[7]；均與本文計算結(jié)果接近。

圖3 2006-2015年研究區(qū)逐年實際/潛在蒸散發(fā)量、降水量及流域來水量Fig.3 Variations of the actual/potential evapotranspiration，precipitation and the incoming water of the study area from 2006 to 2015

對年尺度和季節(jié)尺度實際蒸散發(fā)序列進(jìn)行Mann-Kendall檢驗，結(jié)果如表1所示?？梢钥闯觯昶骄鶎嶋H蒸散發(fā)序列以1.19 mm/10a 的速率呈下降趨勢，由于Zc統(tǒng)計量的絕對值小于1.28（置信度為90%時對應(yīng)的臨界值），這表明該時間序列的變化趨勢沒有通過置信度為90%的顯著性檢驗，即該下降趨勢變化不顯著；季節(jié)尺度上，夏季變化幅度最大，其次是春季，冬季變化幅度最小；閆宇會等在其研究中也發(fā)現(xiàn)，黑河流域夏季的實際蒸散發(fā)相對變化率最大，而冬季最小[7]；馮雅茹研究表明，黑河流域上游夏季實際蒸散發(fā)變化最為明顯，冬季最不明顯[10]。根據(jù)表1 得出，研究區(qū)夏、冬季實際蒸散發(fā)以2.68 mm/10a 和0.01 mm/10a 的速率呈不顯著上升趨勢，而春、秋季則分別以1.15 mm/10a和0.94 mm/10a的速率呈不顯著下降趨勢。

表1 基于Mann-Kendall方法的研究區(qū)年和季尺度實際蒸散發(fā)序列趨勢檢驗結(jié)果Tab.1 Trend results of the annual and seasonal actual evapotranspiration based on Mann-Kendall test

利用Mann-Kendall 突變檢驗法對實際蒸散發(fā)序列的變點特征進(jìn)行分析，結(jié)果如圖4所示?？梢钥闯?，UF和UB兩條曲線在臨界范圍內(nèi)，且存在交點，因此可判斷變點在2010-2011年之間，結(jié)合實際蒸散發(fā)序列隨時間的變化曲線，將2010年確定為變點年份，由此將整個研究期劃分為兩個時段：2006-2010年為時段I，2011-2016年為時段II。時段I實際蒸散發(fā)量為302.66 mm，時段II 實際蒸散發(fā)為281.63 mm，根據(jù)公式（5）得出實際蒸散發(fā)變化量為-22.08 mm，即時段II 較時段I實際蒸散發(fā)量減少7.3%。

圖4 研究區(qū)實際蒸散發(fā)序列的變點特征統(tǒng)計曲線Fig.4 Change point analysis of the actual evapotranspiration in the study area

3.3 實際蒸散發(fā)序列變化的歸因分析

實際蒸散發(fā)的變化受多種影響因素的影響和制約，是能量條件、供水條件、下墊面條件共同作用的結(jié)果。根據(jù)Budyko 理論，研究區(qū)的實際蒸散發(fā)過程與水分和能量密切相關(guān)，水分主要由降水和流域來水量描述，能量主要由潛在蒸散發(fā)描述；降水和流域來水量決定了水分限制的水平，潛在蒸散發(fā)決定了能量限制的水平。

通過計算，實際蒸散發(fā)量與降水量具有良好的正相關(guān)關(guān)系，實際蒸散發(fā)量會隨著降水量的增加而增加，隨著降水量的減少而減少，二者的相關(guān)系數(shù)達(dá)到0.94；其與流域來水量也呈正相關(guān)關(guān)系，但相關(guān)性較低，相關(guān)系數(shù)僅為0.29；降水量和流域來水量是研究區(qū)主要的水分輸入變量，其多寡決定了研究區(qū)實際蒸散發(fā)量的大小以及可利用水分的多少。根據(jù)圖3和表2 可知，時段II 較時段I 降水量減少38.16 mm，減少約25.7%，降水量的減少會導(dǎo)致流域?qū)嶋H蒸散發(fā)量的減少；雖然來水量時段II 較時段I 有所增加，增加7.19 mm，約3.7%，但來水量的增幅遠(yuǎn)小于降水量的降幅，因此由來水量增加導(dǎo)致的實際蒸散發(fā)增加量遠(yuǎn)小于由降水量減少導(dǎo)致的實際蒸散發(fā)減少量。

表3給出了不同因素的彈性系數(shù)以及對實際蒸散發(fā)變化的貢獻(xiàn)率。降水量和來水量的彈性系數(shù)分別為0.31和0.48，即降水量/來水量每增加或減少1%，實際蒸散發(fā)分別對應(yīng)增加或減少0.31%和0.48%；雖然來水量的彈性系數(shù)比降水量的大，但由于來水量的變化較小，因此其對實際蒸發(fā)量變化的貢獻(xiàn)率不及降水量的貢獻(xiàn)率。根據(jù)彈性系數(shù)的結(jié)果，降水量變化對實際蒸散發(fā)變化的貢獻(xiàn)率是-8.94%（負(fù)號表示該因素導(dǎo)致實際蒸散發(fā)減少），來水量變化對實際蒸散發(fā)變化的貢獻(xiàn)率是+1.68%（正號表示該因素導(dǎo)致實際蒸散發(fā)增加），即時段II較時段I 降水量的減少（25.7%）導(dǎo)致研究區(qū)實際蒸散發(fā)量減少了8.94%；來水量的增加（3.7%）導(dǎo)致研究區(qū)實際蒸散發(fā)量增加1.68%。

潛在蒸散發(fā)在兩個時段基本保持不變，減少僅0.2%（圖3和表2）；實際蒸散發(fā)對潛在蒸散發(fā)的彈性系數(shù)為0.21（表3），即潛在蒸散發(fā)每增加或減少1%，則實際蒸散發(fā)增加或減少0.21%。雖然該彈性系數(shù)較大，但由于兩個時段潛在蒸散發(fā)變化很小，因此潛在蒸散發(fā)變化對實際蒸散發(fā)變化的貢獻(xiàn)率很小，僅為-0.04%。

表2 研究區(qū)兩個時段的Ea、P、R1、E0值及變化量Tab.2 The values and the variations of Ea，P，R1 and E0 in the two periods in the study area

表3 研究區(qū)不同因素的彈性系數(shù)及其對實際蒸散發(fā)變化的貢獻(xiàn)率Tab.3 Elasticity coefficients of different factors and their contributions to the changes of the actual evapotranspiration in the study area

除能量和水分這兩個主要控制因子以外，研究區(qū)實際蒸散發(fā)的變化還受到下墊面條件這一次要控制因子的影響；在改進(jìn)的傅抱璞模型中，下墊面條件采用參數(shù)ω來描述，需要說明的是，這里參數(shù)ω被假設(shè)是固定不變的，下墊面條件處于相對穩(wěn)定的狀態(tài)，即這里未考慮下墊面條件對實際蒸散發(fā)變化的影響。

根據(jù)以上分析可知，研究區(qū)供水條件是影響研究區(qū)實際蒸散發(fā)變化的主導(dǎo)因素。由于時段II 較時段I 研究區(qū)降水量有所減少，而來水量有所增加，因此二者對實際蒸散發(fā)的影響恰好相反；降水量減少導(dǎo)致研究區(qū)實際蒸散發(fā)量有所減少，而來水量增加則導(dǎo)致研究區(qū)實際蒸散發(fā)量有所增加；雖然降水量的彈性系數(shù)不及來水量，但由于其減少的幅度更大，因此對實際蒸散發(fā)減少的貢獻(xiàn)最大；來水量的增加促使實際蒸散發(fā)量有所增加，但因來水量增加幅度較小，因此其促使實際蒸散發(fā)增加的比重相對較小。能量條件對研究區(qū)實際蒸散發(fā)變化的影響較供水條件偏弱，且兩個時段相比，潛在蒸散發(fā)的變化不大，因此其對研究區(qū)實際蒸散發(fā)變化的影響非常微弱。

綜上所述，研究區(qū)實際蒸散發(fā)量的減少主要由降水量減少所導(dǎo)致，降水量的大幅度減少抵消了來水量增加對實際蒸散發(fā)的影響，潛在蒸散發(fā)變化對實際蒸散發(fā)減少的影響最為微弱。

4 結(jié)論

基于Budyko 理論以及改進(jìn)的傅抱璞經(jīng)驗?zāi)Ｐ?，結(jié)合流域水文氣象數(shù)據(jù)，估算黑河流域中游地區(qū)2006-2015年的實際蒸散發(fā)量，并對實際蒸散發(fā)序列的趨勢特征、變點特征及變化歸因進(jìn)行了分析。主要結(jié)論如下。

（1）由改進(jìn)的傅抱璞模經(jīng)驗型得到的研究區(qū)多年平均實際蒸散發(fā)量為293 mm，與水量均衡方程得到的結(jié)果基本一致。

（2）年尺度和春、秋季研究區(qū)實際蒸散發(fā)序列呈不顯著下降趨勢，夏、冬季實際蒸散發(fā)序列呈不顯著上升趨勢；根據(jù)Mann-Kendall 突變檢驗結(jié)果，確定2010年為研究區(qū)實際蒸散發(fā)序列的變點年份，2011-2016年較2006-2010年實際蒸散發(fā)量減少7.3%。

（3）對研究區(qū)實際蒸散發(fā)量減少的貢獻(xiàn)最大的是降水量，降水量減少導(dǎo)致實際蒸散發(fā)量減少約8.94%；來水量的增加導(dǎo)致實際蒸散發(fā)量有所增加，增加比例約1.68%，潛在蒸散發(fā)減少導(dǎo)致實際蒸散發(fā)量減少約0.04%；供水條件是影響研究區(qū)實際蒸散發(fā)變化的主導(dǎo)因素，能量條件次之。