基于Budyko假設的汾河上游水源區徑流衰減歸因分析

藺彬彬，張亞瓊，郭維維

（1.太原理工大學水利科學與工程學院，太原030000；2.山西漳河水務有限公司，太原030000）

氣候變化和人類活動的雙重影響是導致徑流規律發生變化的兩大主要因素，氣候變化尤其是降雨導致徑流在時間和數量上都發生了變異，破壞了徑流序列的一致性；人類活動通過改變下墊面條件，使流域產匯流過程發生變化［1］。全球平均氣溫在20世紀約升高了0.6 ℃，IPCC的研究表明全球氣溫在21世紀末可能增高1.1～6.4 ℃。在1957-2003年期間山西省降水量總體呈減少趨勢，減少速率為-17.3 mm/（10 a），顯著高于全國水平；氣溫總體呈上升趨勢，增長率為0.15 ℃/（10 a）［2］。汾河是山西的母親河，黃河的第二大支流，進入20世紀80年代，隨著氣候變化及經濟社會的不斷發展，汾河入黃徑流量衰減明顯［3］，天然徑流量的減少，嚴重影響了流域內經濟社會的發展，加劇了流域內生態環境的惡化。汾河上游作為汾河的重要水源區，也是山西省會太原市重要的地表水水源地，同時也處在巖溶地下水水源涵養區和保護區［4］，近幾十年來汾河上游流域徑流衰減明顯，針對這一事實進行定量分析，對深入理解汾河流域水文演變規律，對未來氣候變化和人類活動加劇背景下水資源適應性管理都具有重要的意義。

針對流域徑流變化的原因，不同學者采用不同的方法進行了分析。例如，劉昌明等［5］應用SWAT 分布式水文模型研究了氣候變化對黃河河源區徑流及蒸散發的影響；馬歡等［6］基于GBHM 模型分析了氣候變化和人類活動對密云水庫入庫徑流量急劇減少的貢獻率；張樹磊等［7］應用Budyko 假設水熱平衡理論對1960-2010年期間我國主要河流上游山區小流域的徑流衰減進行了定量歸因分析，研究發現降雨量的減少和人類活動引起的下墊面的變化是徑流減少的主要原因；張連鵬等［8］以渭河的北洛河流域為研究對象，應用Budyko 假設和TOPMODEL 水文模擬方法定量分析徑流衰減的原因，并兩種方法進行了對比，發現分析的結果具有較好的一致性。

總的來看，目前針對徑流衰減歸因分析的方法主要有兩種：一是基于水文模擬法，二是基于Budyko假設的彈性系數法。而基于水熱平衡理論的Budyko假設方法，由于其方法簡單且輸入的參數少，在定量解析徑流衰減原因方面已得到了廣泛的應用。在流域尺度下，實際蒸散發除受能量供給條件和水分供應條件的影響外，植被、土地利用等下墊面條件也是影響蒸散發變化的重要因素［9］。因此，涉及流域下墊面條件的Budyko 修正模型逐漸發展起來，實現了Budyko 假設的參數化。2006年，楊大文教授［10，11］在已有的蒸散發互補理論研究的基礎上，基于Budyko 假設提出了流域水熱耦合平衡方程，即Choudhury-Yang公式。該公式引入了反映流域下墊面特征的參數n，且表達式相對簡單，已得到廣泛的應用。本文將該公式應用于汾河流域上游水源區，來定量解析汾河上游水源區徑流衰減的原因，為山西省正大力開展的汾河流域清水復流及水源區保護工程提供技術支撐。

1 研究區及數據說明

以汾河上游汾河水庫水文站控制流域為研究區，地理位置如圖1，控制流域面積5 268 km2，屬亞熱帶大陸性季風氣候，為半干旱、半濕潤型氣候過渡區，四季分明，春季多風干燥，夏季多雨炎熱，秋季少晴早涼，冬季少雪寒冷。雨熱同期，光熱資源較為豐富，有利于農業發展。多年平均溫度7.19 ℃，多年平均降水量465 mm，降水年際變化較大，無霜期大于130 d。

采用汾河水庫水文站1961-2016年期間的徑流數據；流域內29 個雨量站的雨量數據；流域內3 個氣象站及周邊8 個氣站的數據，包括降雨、氣溫、日照時長、風速和相對濕度等。根據聯合國糧農組織推薦的Penman-Monteith 公式來，利用11 個氣象站點的氣象數據，計算氣象站點的潛在蒸散發量，利用反距離加權法（IDW）插值生成網格數據求平均，得到流域平均潛在蒸散發量?；?9個雨量站的雨量數據，利用泰森多邊形來計算流域的面雨量。

2 研究方法

2.1 趨勢性和變異性的檢測方法

Mann-Kendall 趨勢檢驗法是世界氣象組織（WMO）推薦并已廣泛應用的一種非參數統計檢測方法，非參數不要求樣本遵循一定的分布，也不受少數異常值的干擾，且計算簡便，常被用來檢測水文氣象長時間序列參數的顯著性趨勢，因此選取此方法來判定徑流、降雨、潛在蒸散發的變化趨勢及徑流變化的突變點。由MK 檢驗法得到統計值Z，當Z＞0 時，說明參數系列呈增加趨勢；當Z＜0時，說明參數系列呈減少趨勢。

2.2 Budyko假設及Choudhury-Yang公式

前蘇聯著名氣候學家Budyko通過研究發現，陸面長期實際蒸散量是由陸面的水分條件和能量條件之間的平衡決定的［10］，并認為可用潛在蒸散發量（PET，簡稱E0）表征流域水循環的能量條件，降水量P表征流域水循環的水分條件?；诖薆udyko提出了陸面實際蒸散量的兩個邊界條件，一個是像沙漠地區的極端干旱情況（E0/P→∞），所有的大氣降水都被用于蒸散發（E/P→1）；另一個是在極端濕潤情況下（E0/P→0），水分供給充分所有可用于蒸散的能量都被用于蒸散發，全部轉化為潛熱（E/E0→1）。

在假定邊界條件的基礎上，Budyko 提出了滿足上述邊界條件的水熱耦合平衡方程的一般形式：

式中：E、P分別為流域多年平均的年實際蒸散發量和降雨量；φ為干旱指數（φ=E0/P），是氣候帶和植被帶劃分的基礎；E0為流域多年平均的年潛在蒸散發量。

理論上，Budyko 框架的水量-能量耦合平衡方程具有普適性，這點得到很多研究的證實，然而仍有很多流域的觀測資料與Budyko 理論曲線存在一定偏差。因此，很多研究者對Budyko 理論曲線模型不斷發展與豐富，提出了不同的估算公式，但大多數公式是基于特定流域推算出來的，具有一定的局限性，至今尚未獲得全球普適的估算方法。對此楊大文等［11］引進了一個參數n來調整因下墊面差異引起的偏差，經過推導得到新的公式：

式中：n為下墊面特征參數，表征了流域植被、土地利用、地形地貌的情況，并認為P，E0，n是相互獨立的變量。

2.3 徑流彈性系數

在一個閉合流域，流域的水量平衡可用下式來表示：

式中：P為流域多年平均降雨量；R為流域多年平均河川徑流量；E為流域多年平均實際蒸散發量；ΔS為時段內流域蓄水量的變化，對于長歷時ΔS可以忽略不計，近似為0。

結合式（1）～（3），流域長歷時年均徑流量可由下式來計算：

年徑流量R的變化可以表示為如下全微分形式：

Schaake［12］將徑流的降雨彈性系數（εP）、徑流的潛在蒸散發彈性系數（εE0）、徑流的下墊面彈性系數（εn）；分別定義為εP=

將式（5）除以多年平均徑流深R，可以得到：

利用式（4）分別對參數P、E0、n求偏導，求得彈性系數εP，εE0，εn表達式如下：

這3個彈性系數反映了流域多年平均的水文氣候和下墊面特征，如果設定εP，εE0，εn的值分別為a，b，c，那么εP表示：如果P增加1%，將驅動徑流量R增加a%（或減少b%）；εE0表示：如果E0增加1%，將驅動徑流量R減少b%；εn表示：如果n增加1%，將驅動徑流量R減少c%。

2.4 徑流變化歸因分析

在人類活動和氣候變化的影響下，徑流總的變化可以表示為：

式中：ΔRtot為徑流總的變化量；為基準期多年平均徑流量；為變化期多年平均徑流量。

總的徑流變化量可以表示為：

根據彈性系數εP，εE0，εn，可分別按下式求得ΔRP，ΔRE0，ΔRn。

式中：ΔP，ΔE0，Δn分別表示流域基準期和變化期年均降雨量、年均潛在蒸發量、下墊面參數的變化量；nbas和nvar分別代表基準期和變化期下墊面參數，可由式（2）反推得到。

3 結 果

3.1 徑流、降雨及潛在蒸散發趨勢分析

對1961-2016年汾河水庫站56年的年徑流系列進行統計分析，通過趨勢線及5年滑動平均分析（圖2），通過M-K 檢驗方法，計算年徑流系列的MK 統計值為-3.79，都說明年徑流呈顯著下降趨勢，且通過了0.05顯著性水平檢測，年徑流深以8 mm/（10 a）的速率在遞減。通過M-K 突變分析（圖3）發現1961-2016年期間汾河水庫站的年徑流量在1980年發生了突變，因此將1980年設置為突變點。據此，在進行徑流衰減歸因分析時，將1961-1980年劃分為基準期，將1981-2016年劃分為變化期，基準期和變化期年徑流深的平均值分別為77.8 mm 和50.2 mm，徑流深衰減了27.6 mm，將近35.5%。

圖2 汾河水庫站年徑流趨勢分析Fig.2 Long-term trend of annual-runoff at Fenhe reservoir station

圖3 汾河水庫站年徑流M-K突變分析Fig.3 M-K mutation analysis of annual-runoff at Fenhe reservoir station

研究區年降雨量MK 統計值為-0.163，說明年降雨量呈現不顯著的下降趨勢，長系列趨勢分析見圖4。研究區年潛在蒸散發量MK 統計值為1.97，說明年潛在蒸散量呈增加的趨勢，長系列趨勢分析見圖5。基準期和變化期研究區年均徑流深、年均降雨量、年均潛在蒸散發量的統計見表1。

表1 基準期和變化期的參數統計表 mmTab.1 The value of R，P，E0 and n in base-period and change-period

圖4 研究區年降雨量趨勢分析Fig.4 Long-term trend of precipitation

圖5 研究區年潛在蒸發量趨勢分析Fig.5 Long-term trend of potential evapotranspiration

3.2 徑流彈性系數計算

根據研究區1961-2016年的年降雨量、潛在蒸散發、徑流深，推求長系列年均值，計算干旱指數，根據公式（4）求解研究區下墊面參數n，根據公式（7）～（9）分別求得徑流的3 個彈性系數，計算結果見表2，說明當流域年降雨量增加（減少）1%時，將導致徑流量增加（減少）2.64%；年潛在蒸散發量增加（減少）1%時，將導致徑流量減少（增加）1.64%；下墊面參數n增加（減少）1%時，徑流流量減少（增加）1.89%。

表2 研究區特征及徑流彈性系數Tab.2 The characteristic and runoff elasticity coefficient of research area

3.3 徑流變化歸因分析

在1961-2016年期間，徑流發生突變的1980年前后，基準期和變化期徑流深（R）、降雨（P）、潛在蒸散發（E0）及下墊面參數（n）的統計值見表3。變化期和基準期的下墊面參數n，分別根據公式（4）求解得到?；谘芯繀^基準期和變化期年均降雨、年均潛在蒸散發及下墊面彈性系數的值，利用求得的彈性系數，根據公式（14）分別計算由于三者驅動引起的徑流變化量，并分別計算貢獻率。從表3可以看出汾河水庫站控制流域，變化期相對于基準期多年平均降雨量減少了31.6 mm，驅動徑流量減少了10.6 mm，貢獻率為39.7%；潛在蒸散量增加了5.7 mm，驅動徑流量減少了0.7 mm，貢獻率為2.5%；兩者之和就是氣候變化對徑流衰減的貢獻率為42.2%。而人類活動引起的下墊面改變，導致變化期下墊面特征參數n增加了0.269，驅動徑流量減少了15.4 mm，對徑流衰減的貢獻率為57.8%。由此可說明，1961-2016年汾河流域上游水源區徑流衰減的主要原因是人類活動引起的下墊面變化，其次是降雨的減少。

表3 徑流變化歸因分析 mmTab.3 Attribution analysis of runoff attenuation

變化期與基準期相比下墊面參數n增加了0.269，說明在人類活動的影響下汾河上游下墊面情況發生了較大的變化。從20世紀80年代，在政府主導下的汾河上游進行了大規模的水土保持措施，主要包括大規模的植樹造林及退耕還林還草，建設淤地壩，修建基本農田。邸富宏［13］基于MODIS數據對汾河上游植被動態進行了監測，研究發現2000-2010年來，汾河上游區域NDVI最大值呈上升趨勢，NDVI隨年份的增長率為7.8%/（10 a），植被覆蓋明顯改善；黨晉華等［14］對汾河上游區域土地類型變化進行了分析，發現2000-2013年森林為最活躍的土地利用類型，濕地和森林在空間上呈現擴張的發展趨勢。植被覆蓋的增加在改善區域生態環境、水源涵養及治理水土流失方面具有重要的意義，但由于流域實際蒸散發的增加，導致產流量的減少。另外汾河上游部分區域地處山西六大煤田之一的寧武煤田，當地的煤礦開采，形成大面積的采空區，導致地表變形、塌陷，地層中的裂隙增多增大，地表形成大量的裂縫，降雨入滲量增大產流系數減少。

4 結 論

汾河上游作為汾河流域重要的水源區，1961-2016年期間年徑流量呈現顯著的下降趨勢，變化期相比基準期減少了35.5%，應用基于Budyko假設的水熱耦合平衡理論，對徑流衰減的原因進行了解析。研究發現，汾河上游由于人類活動引起的土地利用、植被覆蓋、地形等下墊面特征參數的改變是導致上游徑流衰減的主要原因，相應的貢獻率接近60%，其次是降雨減少導致的。汾河上游大規模開展的水土保持措施及長期的煤礦開采，導致下墊面條件發生了較大的變化，從而增加了流域的實際蒸散發，導致降雨入滲量增多。研究表明制定科學合理水土保持措施，低影響的煤炭開采方式及采空區的修復治理，對現在山西省正在大力開展的汾河流域水生態修復及清水復流工程具有重要的意義。