基于Budyko假設的湟水徑流變化歸因識別

杜嘉妮，蔡宜晴，劉希勝，王 崗，閔 敏

（青海省水文水資源測報中心，西寧 810001）

0 引 言

地表徑流作為水循環(huán)的重要環(huán)節(jié)，是揭示水資源的重要指標。近年來，在全球暖化背景下，氣侯變化使流域水循環(huán)和水文過程產生改變，水資源的時空分布也隨之發(fā)生改變［1］；同時隨著社會經濟的快速發(fā)展，人類活動如城鎮(zhèn)化發(fā)展、水利工程建設、水土保持工程實施等，使下墊面條件發(fā)生了顯著變化［2］，改變了流域產匯流機制。因此變化環(huán)境下水循環(huán)成為水科學領域研究的熱點問題，研究流域水文過程與演變機制，評價各驅動要素的貢獻，對流域水資源管理、水資源可持續(xù)利用及促進社會經濟發(fā)展具有重要的現實意義。

目前定量識別徑流變化歸因常用的方法主要為水文模型法和基于Budyko 假設的敏感性方法。水文模型法基于物理過程模擬驅動因素的貢獻率，物理概念清晰，但受模型結構、模型參數、初始和邊界條件的影響，存在一定的不確定性，可能會對模擬結果產生較大誤差?；贐udyko 假設的敏感性方法主要基于年尺度進行研究，物理意義明顯，對數據要求低，并且能對各要素的貢獻率單獨進行計算，是一種簡單有效的分析方法，近年來被廣泛應用［3，4］。

湟水位于黃河上游，是青海的母親河，是青海省東部地區(qū)工農業(yè)的主要水源，也是該地區(qū)經濟社會發(fā)展的根本保障。近年來，由于氣侯變化、人口增長和經濟快速發(fā)展，湟水水資源供需矛盾突出、水污染加劇、水能資源無序開發(fā)、水土流失嚴重，水事矛盾尖銳，如何應對環(huán)境變化帶來的影響是水資源管理面臨的關鍵問題之一。目前有關湟水徑流的研究主要著重于徑流模擬與預測［5-7］，研究變化環(huán)境下各要素對徑流變化貢獻率的相對較少，張調風［8］和劉希勝［9］分別采用累積量斜率變化率法和雙累積曲線定量評價氣侯變化和人類活動對湟水徑流變化的貢獻率，上述方法均為統(tǒng)計分析法，缺乏與其他方法的比較。本文進一步采用基于Budyko假設的敏感性方法，定量識別氣侯變化和人類活動對徑流的貢獻（影響），為實現變化環(huán)境下水資源管理、優(yōu)化配置和可持續(xù)利用提供技術支持。

1 研究區(qū)概況及數據來源

1.1 研究區(qū)域概況

湟水是黃河上游一級支流，發(fā)源于青海省海晏縣包呼圖山，流經青海、甘肅兩省，于甘肅省永靖縣注入黃河。全長374 km，流域面積1.77 萬km2，其中青海省境內流域面積1.61 萬km2，多年平均天然徑流量21.56 億m3。氣候屬高原干旱半干旱氣候，多年平均氣溫2.1～8.1 ℃，其中川水地區(qū)多年平均氣溫3.3～8.1 ℃，是流域內最暖地區(qū)；淺山和河谷地區(qū)多年平均氣溫2.1～3.2 ℃，屬冬寒夏涼的半干旱和干旱氣候區(qū)；腦山地區(qū)多年平均氣溫為1 ℃，屬高寒半濕潤山地氣候。降水量分布由河谷向兩側山區(qū)遞增，河谷地區(qū)降水一般在300～400 mm，山區(qū)則達到500～600 mm。水面蒸發(fā)量空間變化與降水量相反，隨海拔的升高而降低，河谷地區(qū)多年平均水面蒸發(fā)量在900 mm 以上，而山區(qū)則在700 mm左右。

1.2 數據來源

民和水文站數據來源于黃河水利委員會水文局，其他水文站和雨量站的數據來源于青海省水文水資源測報中心，氣象數據來源于青海省氣象局。本文選用湟水民和水文站1957-2017年實測徑流數據作為研究徑流變化的基礎數據。降水采用流域內81 個水文氣象站點的降水量，系列長度1957-2017年。潛在蒸散發(fā)采用流域內氣象站的逐日最高氣溫、最低氣溫、日照時數、相對濕度、10 m 處風速及氣壓資料，利用FAO 修訂的Penman-Monteith 方法［10］計算，系列長度1961-2017年。面降水和潛在蒸散發(fā)通過泰森多邊形法計算求得。收集到區(qū)域1986、1990、1995、2000、2005、2010、2015 和2017年月度合成30 m×30 m 地表植被指數數據［11］。該數據采用最大值合成（Max value composition，MVC）方法，利用Landsat5，Landsat8 和sentinel2 紅光和近紅外兩個通道的反射率數據，實現對地表月度NDVI 產品的合成。數據通過Google Earth Engine 云平臺對反射率進行月度合成，來源于國家青藏高原科學數據中心（http：//data.tpdc.ac.cn/）。

圖1 湟水示意圖Fig.1 The overview of the Huangshui River

2 研究方法

2.1 線性傾向估計趨勢分析

yi表示樣本為n的變量，xi為yi所對應的時間，建立兩者之間的線性關系，即：

式中：a為回歸系數；b為回歸常數；a和b利用最小二乘法進行計算。

利用a與相關系數之間的關系，求出yi和xi的相關系數R。

當a＞0（＜0）時，表明變量y隨x變化呈增加（減少）趨勢。相關系數R表示y和x之間相關程度，給定顯著性水平α，若|R|＞Rα，表明變量y隨時間x變化趨勢顯著；否則變化趨勢不顯著。

2.2 突變檢驗

目前突變檢驗的方法很多，不同的方法原理不同，可能得出的結果也不一樣，為了合理確定判別突變點，需對多種方法比較以確定。本文采用累積距平和滑動t檢驗［12］判別徑流突變點。

2.2.1 累積距平

累積距平常用來判斷序列的變化趨勢，有時也用來診斷發(fā)生突變的大致時間。對于序列x，其某一時刻t的累積距平可表示為：

其中：

2.2.2 滑動t檢驗

該方法是通過檢驗兩系列均值是否有顯著性差異來判別突變，如果兩序列的均值差異超過了一定的顯著性水平，可認為均值發(fā)生了質變，有突變發(fā)生。

將一樣本為n的序列x分成兩個樣本為n1和n2的子序列x1和x2，定義統(tǒng)計量：

式中：

給定顯著性水平α，若|t|＞tα，表明存在顯著性差異，|t|最大值對應的點即認為發(fā)生了突變。為避免子序列長度不同引起突變點漂移，反復變動子序列的長度進行比較，以提高計算準確度。文中n1=n2=10，顯著性水平α取0.05。

2.3 氣候和下墊面對徑流的貢獻率

2.3.1 Budyko假設

多年平均情況下，一個流域的水量平衡方程可表示為：

式中：Q、P、E分別代表流域多年平均徑流量、降水量、實際蒸散發(fā)量，mm。

Budyko［13］假設認為流域實際蒸散發(fā)受水分（主要為降水）和能量（主要為潛在蒸散發(fā)）兩個供應條件的影響，并給出了反映水分和能量關系的耦合平衡方程，即：

式中：E0代表流域多年平均潛在蒸散發(fā)量，mm；P、E代表含義同上。

基于Budyko 假設，1981年傅抱璞［14］通過量綱分析和微分理論得到實際蒸散發(fā)的計算公式：

式中：E、E0、P代表含義同上；n是流域下墊面特征參數，主要反映下墊面變化對流域水量平衡的影響。

結合式（9）、（11），流域多年平均水量平衡方式可表示為：

式中：P、E0、Q已知，可求出n。

2.3.2 徑流對氣候的敏感性

由式（12）可知，徑流Q可以表示P、E0、n的函數，結合敏感性系數的定義：ε=表示徑流Q對獨立變量X的敏感性系數，X代表P、E0、n。

徑流對各變量的敏感性系數表示如下：

徑流對各變量的敏感性定義為單位變量的變化導致徑流的變化程度，敏感性系數為正（負）表示徑流隨變量的增加而增加（減小），通過式（13）～（15）可分別求得徑流對各變量的敏感性系數。

2.3.3 徑流變化貢獻率分析

各變量對徑流的影響表示為：

式中：ΔQ＇為計算得到的徑流變化量。

根據徑流突變點，將徑流序列分為基準期和變化期，基準期的徑流為Q1，變化期的徑流為Q2，徑流的變化量ΔQ可以用兩時段的徑流之差來表示：

式中：ΔQ為實測的徑流變化量。

各變量對徑流的貢獻率表示為：

式中：X表示P、E0、n，C表示各變量對徑流的貢獻率。

3 結果分析

3.1 徑流變化

湟水民和站徑流年內分配不均勻，年內最大徑流量2.601億m3，出現在9月；最小徑流量0.609 億m3，出現在3月，最大月徑流量是最小月徑流量的4.27 倍。連續(xù)最大四個月徑流量出現在7-10月，占全年徑流量的57.8%。

湟水民和站年徑流呈減少趨勢，但趨勢不顯著，變化傾向率為-0.60 億m3/10 a，與黃河流域徑流減少的趨勢一致。多年平均徑流量15.91 億m3，其中年徑流量最大值31.12 億m3，出現在1961年，最小值7.089 億m3，出現在1991年，極值比4.39，年際變化較大。20世紀70年代、90年代和21世紀前10 a 徑流較多年平均分別偏少8.10%、13.3%、12.7%；20世紀60年代和80年代徑流較多年平均偏多14.4%和10.8%；2010-2017年徑流與多年平均基本持平，偏少0.27%，徑流變化呈現出豐-枯-豐-枯-枯-平的年代際變化規(guī)律。

3.2 趨勢及突變分析

圖2 湟水徑流年際變化Fig.2 Interannual variation of runoff in the Huangshui River

為了進一步了解湟水徑流變化趨勢，采用累積距平判斷變化趨勢。如圖3（a）所示，湟水1957-2017年徑流經歷了二次明顯的波動，主要表現在1957-1968年和1981-1990年徑流增加，1969-1980年和1991-2004年徑流減少，2004年以后徑流減少趨勢變緩，接近多年平均水平。由此看出，湟水徑流在1968年、1980年、1990年和2004年發(fā)生轉折。

累積距平可以初步判斷出徑流突變點發(fā)生的時段，為了更準確確定突變點，進一步采用滑動t檢驗來識別突變點，如圖3（b）所示，湟水徑流T統(tǒng)計值有三處超出了0.05 顯著性水平，出現在1980年、1990年和2004年，雖然1968年徑流經歷了由增加到減少的轉折，但未超出顯著性水平。綜合確定湟水徑流突變年份為1980年、1990年和2004年。

圖3 湟水徑流累積距平和滑動t檢驗Fig.3 Accumulative anomaly and moving t test of runoff in the Huangshui River

根據湟水徑流突變年份，將徑流系列劃分為基準期和變化期，湟水徑流基準期為1957-1980年，變化期為1981-1990年、1991-2004年和2005-2017年。

3.3 徑流的敏感性分析

表1給出了湟水基準期和變化期各變量的特征值以及徑流對各變量的敏感性系數。從表1中可以看出變化期1981-1990年相比其他時段，降水、徑流、徑流系數較大，而潛在蒸散發(fā)和干旱指數較小，說明變化期1981-1990年氣候較其他時段濕潤。

表1 湟水不同時段水文氣象特征值Tab.1 Hydro-climate characters of different times in the Huangshui River

敏感性系數的符號有正有負，正（負）值表示與徑流呈正（負）相關，說明Q與P正相關，與E0和n負相關。數值的大小則表示對徑流的影響程度，?Q/?P、?Q/?E0和?Q/?n的范圍分別在0.40～0.47、-0.16～-0.13 和-105.6～-77.6，表明P增加1%將導致Q增加0.40%～0.47%，E0增加1%將導致Q減少0.13%～0.16%，n增加1%將導致Q減少77.6%～105.6%。從3個變量的敏感系數絕對值可以看出，徑流對下墊面（n）的變化最為敏感。

3.4 徑流變化歸因識別

各變量對徑流的影響見表2。計算得到的徑流變化值（ΔQ＇）與實測的徑流變化值（ΔQ）相差-1.82～0.21 mm，誤差范圍為4.47%～7.46%，說明本文評估氣侯和人類活動對徑流影響的方法基本可靠。

表2 徑流變化歸因分析Tab.2 The attribution analysis of runoff changes

湟水徑流變化期不同時段降水、潛在蒸散發(fā)和下墊面對徑流變化的影響程度不同。變化期1981-1990年降水對徑流的影響最顯著，其次是潛在蒸散發(fā)，下墊面的影響最小。變化期1991-2004年和2005-2017年下墊面對徑流的影響最顯著，其次是降水，潛在蒸散發(fā)的影響最小。綜上所述，湟水1981-1990年降水是影響徑流增加的主要因素，1990年以后下墊面變化是導致徑流減少的主要因素。這與以往學者的研究基本一致［8，9，15］，如楊大文等人［15］的研究表明，黃河流域各子流域徑流突變發(fā)生在20世紀80年代末和90年代初，對大多數流域，降水對徑流減少的貢獻率為49.3%，潛在蒸散發(fā)對徑流減少的貢獻率為-3.5%，下墊面變化對徑流減少的貢獻率為72.0%，下墊面變化對徑流減少起主導作用。

湟水土地面積僅占青海省土地面積的2.3%，卻承載著全省約60%的人口，人口密集且人類活動頻繁。根據全國第一次水利普查統(tǒng)計的湟水水庫相關數據，截至1980年水庫累積庫容6 730.3 萬m3，1990年水庫累積庫容12 976 萬m3，2004年水庫累積庫容達到35 927 萬m3。水庫的修建使徑流的年內分配發(fā)生變化，攔蓄洪水削減洪峰，另外還會導致蒸發(fā)和滲漏增加。根據青海省第三次水資源調查評價資料，民和水文站以上湟水年取用水量1980年以前為3.534 億m3，1981-1990年為5.266 億m3，1991-2004年為6.244 億m3，2005-2017年為6.718 億m3。由此可以看出，隨著社會經濟的發(fā)展，人口逐年增長和城市化、工業(yè)化進程的不斷加快，以及農業(yè)灌溉面積的增加，生產生活用水增加對湟水徑流的影響越來越大，使徑流量不斷減少。

湟水干流西寧以下的黃土丘陵區(qū)水土流失嚴重，是水土保持的重點區(qū)域，先后實施了一系列生態(tài)保護工程和綠化工程。根據全國第一次水利普查水土保持措施數據，截至2010年湟水累計完成水土保持治理面積2 092.2 km2，其中，基本農田中梯田面積1 246.0 km2，水土保持林（喬木林、灌木林）663.2 km2，經濟林8.1 km2，種草43.2 km2，封禁治理131.7 km2，淤地壩建設661 座，小型蓄水保土工程62 864 個。經過治理，水土流失面積減少，湟水干流輸沙量呈顯著減少趨勢，并且下游站點的減少趨勢強于上游站點［9］，反映了實施水土保持措施的效果。水土保持措施的實施大大改變了徑流原來的產匯流特性與路徑，林草地較荒地截留和持水能力強，梯田較坡耕地入滲能力強，這些使降水不易直接形成徑流，導致地表徑流量減少。

我國從1999年開始推行退耕還林還草政策，青海省按照黨中央、國務院的部署，自2000年開始實施這一重大舉措，湟水過度墾殖和超載放牧導致生態(tài)環(huán)境惡化，被列入重點實施的區(qū)域，受資料獲得限制，未收集到流域退耕還林（草）面積數據，根據收集到的流域NDVI數據，可以看出NDVI年均值和生長季均值均呈增加趨勢（見表3），表明地表植被覆蓋度在增加。張曉明等人［16］認為，黃土丘陵溝壑區(qū)森林植被增加對徑流有顯著的調節(jié)作用，不論是一場暴雨、雨季還是年際尺度，均使流域產流能力下降，導致徑流量減少。

表3 湟水歸一化差值植被指數（NDVI）信息表Tab.3 The NDVI information in the Huangshui River

總的來說，1990年以后人類活動使湟水下墊面發(fā)生變化主要體現在水利工程建設、退耕還林（草）和水土保持措施實施，這些共同作用于河川徑流，使徑流量趨于減少，而這種減流作用有可能隨著人類活動加劇而增強。

4 討 論

氣侯和下墊面是影響流域徑流的重要因素，識別二者對徑流的影響時，關鍵在于將二者的影響進行分離，無論是基于Budyko假設的敏感性方法還是水文模型，在分離氣候變化和下墊面對徑流的影響時，都假設氣候變化和下墊面是兩個相互獨立的變量。而事實上，氣候變化和下墊面變化并非完全獨立，是相互影響的。如氣溫升高，蒸發(fā)加劇，土壤水分及下滲將隨之改變，反之土壤濕度及植被的變化也會影響土壤增發(fā)、植被蒸騰和降水。目前的水量平衡方程中，未考慮蓄水量的變化，忽略了水利工程的攔蓄影響，近幾十年來人類活動加強，蓄水量也會有很大的變化，忽略長期蓄水量的變化有可能會對計算結果產生誤差。還有將水熱耦合平衡方程中的下墊面作為一個整體考慮，下墊面參數（n）的物理意義和數學表達式尚未明確。在未來研究中需進一步考慮水利工程攔蓄影響以及量化下墊面因子，以便更準確地評價人類活動對徑流的影響。

5 結 論

本文以湟水為研究對象，采用線性傾向估計、累積距平、滑動t檢驗分析徑流變化特征，基于Budyko 假設的敏感性方法定量評價氣侯和下墊面對徑流的影響，主要結論如下：

（1）湟水徑流年內分配不均，連續(xù)最大四月徑流量占全年徑流量的57.8%；年徑流總體呈不顯著減少趨勢，變化傾向率為-0.60 億m3/10 a；年代際徑流變化呈現出豐-枯-豐-枯-枯-平的變化規(guī)律。

（2）利用累積距平和滑動t檢驗對徑流進行突變檢測，湟水徑流突變年份為1980年、1990年和2004年。根據突變年份將徑流系列劃分為基準期和變化期，湟水徑流的基準期為1957-1980年，變化期 為1981-1990年、1991-2004年 和2005-2017年。

（3）采用基于Budyko假設的敏感性方法計算徑流對各變量的敏感系數，結果顯示，徑流對下墊面的變化最為敏感；湟水1981-1990年降水是影響徑流增加的主要因素，1990年以后下墊面變化是導致徑流減少的主要因素，人類活動導致下墊面發(fā)生變化主要體現在水利工程建設、退耕還林（草）和水土保持措施實施。