基于Budyko假設(shè)的汾河上游水源區(qū)徑流衰減歸因分析

藺彬彬，張亞瓊，郭維維

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

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

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

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

1 研究區(qū)及數(shù)據(jù)說明

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

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

2 研究方法

2.1 趨勢(shì)性和變異性的檢測(cè)方法

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

2.2 Budyko假設(shè)及Choudhury-Yang公式

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

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

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

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

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

2.3 徑流彈性系數(shù)

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

式中：P為流域多年平均降雨量；R為流域多年平均河川徑流量；E為流域多年平均實(shí)際蒸散發(fā)量；ΔS為時(shí)段內(nèi)流域蓄水量的變化，對(duì)于長(zhǎng)歷時(shí)ΔS可以忽略不計(jì)，近似為0。

結(jié)合式（1）～（3），流域長(zhǎng)歷時(shí)年均徑流量可由下式來計(jì)算：

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

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

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

利用式（4）分別對(duì)參數(shù)P、E0、n求偏導(dǎo)，求得彈性系數(shù)εP，εE0，εn表達(dá)式如下：

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

2.4 徑流變化歸因分析

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

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

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

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

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

3 結(jié) 果

3.1 徑流、降雨及潛在蒸散發(fā)趨勢(shì)分析

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

圖2 汾河水庫站年徑流趨勢(shì)分析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

研究區(qū)年降雨量MK 統(tǒng)計(jì)值為-0.163，說明年降雨量呈現(xiàn)不顯著的下降趨勢(shì)，長(zhǎng)系列趨勢(shì)分析見圖4。研究區(qū)年潛在蒸散發(fā)量MK 統(tǒng)計(jì)值為1.97，說明年潛在蒸散量呈增加的趨勢(shì)，長(zhǎng)系列趨勢(shì)分析見圖5?；鶞?zhǔn)期和變化期研究區(qū)年均徑流深、年均降雨量、年均潛在蒸散發(fā)量的統(tǒng)計(jì)見表1。

表1 基準(zhǔn)期和變化期的參數(shù)統(tǒng)計(jì)表 mmTab.1 The value of R，P，E0 and n in base-period and change-period

圖4 研究區(qū)年降雨量趨勢(shì)分析Fig.4 Long-term trend of precipitation

圖5 研究區(qū)年潛在蒸發(fā)量趨勢(shì)分析Fig.5 Long-term trend of potential evapotranspiration

3.2 徑流彈性系數(shù)計(jì)算

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

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

3.3 徑流變化歸因分析

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

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

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

4 結(jié) 論

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