基于時變Budyko 公式的年徑流變化機理分析
——以貴州省六硐河為例

王明陽，董前進，張艷敏，黃 馗

（1. 武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072； 2. 中工武大設計集團有限公司，湖北 武漢 430070；3. 廣西電網電力調度控制中心，廣西 南寧 530013）

0 引 言

近年來，全球變暖趨勢加劇，打破了持續穩定的水循環系統，加劇了水資源分布不平衡現狀［1］。同時，人類活動對流域產匯流等水資源循環過程產生了重大影響［2］。河川徑流是水循環過程的重要環節，它反映了集水區在時間和空間上的綜合變化［3］。同時，作為生態、經濟和社會發展的基本需求，河川徑流的變化給人類社會的生產、生活、生態帶來了極大影響。因此定量分析河川徑流的變化及驅動機制是水資源管理中的重要基礎工作，可以從理論上幫助我們適應未來氣候條件的演化，拓展非一致性研究成果［4］。

當前已有較多關于徑流變化歸因分析的研究，主要分為兩類：一類是利用傳統數學分析方法，如邱玲花等［5］采用雙累積曲線法和敏感性系數法對黑河中游流域的徑流進行氣候變化和人類活動影響的定量識別與評估；一類則利用已有的分布式水文模型，如孫興國等［6］基于SWAT 模型從下墊面變化和水利工程影響對莊浪河流域徑流變化進行分析；李靜［7］等基于MIKE SHE 分布式模型對喀斯特流域進行徑流模擬。這些方法對輸入數據一般有較高的要求，需要對多參數進行擬合，但在長時間尺度上往往不能獲得較好的效果。

基于水熱耦合平衡原理（Budyko 假設）的徑流變化歸因分析方法因其一方面考慮了一定的物理成因，另一方面參數較為簡單，由此在徑流變化歸因長時間尺度分析上得到了廣泛應用。Budyko 假設自1974 年提出后，許多學者基于該假設提出新的理論，如傅抱璞結合水量平衡并利用量綱分析推導出Fu equation［8］，在實踐方面被驗證具有很好的模擬效果；LI［9］等人基于Budyko 假設并使用遙感數據定量估算了無定河流域氣候變化和人類活動對徑流量的影響；ZHANG［10］等人基于Budyko假設利用徑流受氣候的敏感性，估算了黃河源區氣候變化和土地覆蓋植被的變化；XU［11］等利用Choudhury-Yang 公式，分析了徑流的敏感性。

由此，基于Budyko假設對六硐河流域進行年徑流變化歸因分析。根據六硐河流域數據獲取情況及水文情勢變化產生的物理成因，分別建立由一定物理機制推求的和由時間推求的兩類時變Budyko 公式，以Nash 效率系數為評價指標，驗證最適應研究區的水熱耦合公式，該研究將有助于認識貴州省黔南地區流域徑流變化機理。

1 研究區概況

六硐河為珠江水系西江干流紅水河段支流，位于東經106°54'～107°38'，北緯25°11'～26°21'。流域自北向南流經貴州黔南獨山縣、平塘縣、都勻等市，其干流發源于獨山縣拉林鄉八壩，流出貴州省于三堡匯合曹渡河，最終匯入紅水河。六硐河全長164 km，研究區總面積1 441 km2。流域屬于亞熱帶季風性濕潤氣候，夏季高溫多雨，冬季溫和少雨。年均氣溫15.0～17.0 ℃，流域年降水量1 100～1 300 mm，河口多年平均流量56 m3/s［12］。

六硐河流域地勢北高南低，上游地區多為丘陵區，植被覆蓋程度高，且河谷開放。流域周邊是小盆地區，多為耕地，是主要農業種植區，中下游多為高山，流域普遍被大山包圍，且河谷狹窄，河床險峻。主要巖類為石炭系、碳酸鹽巖，呈明顯的喀斯特地貌。

2 研究方法

由于六硐河流域面積較小，可用的水文站和氣象站點只有平湖站。因此，本文徑流數據采用六硐河流域平湖站1960-2012 年逐日平均流量，氣象數據包括平塘站1960-2012 年逐日降水量、日最高、最低氣溫、日平均溫度、日照時數等。采用Hamon 經驗公式［13，14］計算獲得流域1980-2012 年潛在蒸散發量，該公式如下：

式中：N為日照時間數，h；Tα為日平均溫度，℃；E0為潛在蒸散發量，mm。

2.1 Budyko假設

在研究全球水量和能量平衡分析時，Budyko 認為地表年尺度上蒸散發量由當地降水補給和當地輻射能量供給的平衡關系決定，進而提出Budyko理論［15］，公式表示為：

式中：φ定義為某地區的干燥指數，其值為；Bk（）為Budyko理論下的函數形式，其中廣泛應用的Choudhury-Yang公式［17］表達如下：

式中相關符號含義同公式（2），ω為流域特征參數，表示流域下墊面特征的綜合參數。

2.2 時變Budyko公式

原始Budyko公式是不含參數的［17］，現有研究表明引入特定參數得到改進的Budyko公式更能反映流域實際情況，且模擬性能更佳。時變Budyko 公式是對公式（4）的一種改進。在時變Budyko 公式中，流域特征參數ω不再為定值，而是隨時間或者其他解釋變量變化的量。

構建的時變Budyko公式可用如下公式表示：

式中：ωt為時變參數，x=(x1，x2，…，xn)表示時變參數的解釋變量；α和β是解釋變量的系數；f(α，β，x)表示時變參數與解釋變量的函數關系。

根據解釋變量與函數類型的不同，將時變參數ωt與解釋變量的函數關系劃分為有一定物理機制的Budyko 公式和隨時間變化的Budyko公式。前者采用較為簡單的線性關系，其時變參數ωt可按下式表示：

解釋變量x表示流域降水、輻射能量、溫度、人均取用水量等有一定物理含義的變量，使得在該公式下的流域特征參數ω與地區人類活動有一定的聯系，β為相關系數矩陣。結合1981-2012年國家統計年報，選擇年平均降水、年平均溫度數據、地區人均生產總值這3個變量作為有一定物理機制時變參數的解釋變量。

隨時間變化的Budyko 公式一般以線性和二次多項式函數為主如下式所示：

采用最小二乘法對以上兩類表達式進行擬合，并對比不同公式在模擬六硐河流域徑流的精度，確定最適合研究區的徑流變化計算方法，從而進行最終的徑流變化分析。

2.3 徑流變化歸因分析

在長時間尺度（如年）徑流變化歸因分析上通常可忽略流域蓄水量的變化，由此結合水量平衡方程和公式（4），可得到徑流變化的表達式如下［18，19］：

其中各項的偏微分形式如下：

為統一彈性系數公式，以干燥指數φ對式（10）進行簡化，有：

在已有流域潛在蒸散發量和降水的前提下，由時變Budyko公式模擬徑流的表達式為：

式中：Rt表示估算的各年徑流深；Pt指各年降水量；φt即各年干燥指數。

由式（16）并通過離差平方和最小可得到最優的流域特征參數值ω。

通過M-K 趨勢檢驗與突變分析將已有時間序列區劃為兩個時期，即基準期和變化期［20］。同時令基準期的年平均降雨量為P1，變化期年平均降雨量為P2，則流域徑流量變化ΔR受地區氣候變化以及流域特征參數變化的影響如下式表示［18］：

其中ΔR代表兩個時期徑流變化的差值；ΔRp、ΔREo和ΔRω分別代表P、Eo和ω變化導致的徑流變化，為簡便表達下文均以xi替代。

則流域各影響因素對流域徑流變化的貢獻率η為［21］：

3 結果與討論

3.1 六硐河流域水文氣象要素趨勢變化與突變檢驗

基于M-K趨勢分析與突變檢驗得到六硐河流域降水、徑流深及潛在蒸散發量等水文氣象要素的變化趨勢與突變點檢驗成果如圖2 所示。由圖2 可知：六硐河流域降水呈下降趨勢，其突變時間在2005 年左右［圖2（a）］；流域徑流深呈下降趨勢，其突變時間為1997 年［圖2（b）］；而潛在蒸散發總體呈顯著（顯著性水平為0.05）上升趨勢。考慮到降水量對流域徑流的影響，以1981-1996 年為基準期，1997-2012 年為變化期進行本文的研究。

圖2 各水文氣象要素趨勢分析與突變檢驗Fig.2 Trend analysis and mutation test of hydrological elements from 1980 to 2012

3.2 多種Budyko公式擬合參數結果與性能對比

根據站點基準期與變化期多年平均降水量、多年平均徑流深以及計算得到的多年平均潛在蒸散發量，利用公式（16）計算各時期對應的流域特征參數ω，并得到ω基于各類不同影響因素考慮下的擬合方程。為了客觀對比各類方程下參數ω的擬合效果，以Nash-Sutcliffe 模型效率系數為評估模型性能的方法。若Nash系數越接近1，則表明擬合效果越好［22］。

由擬合結果數據來看，有一定物理機制變量推求的時變Budyko 公式估算流域徑流深的Nash 系數最高，為0.923。該模型將年降水、年均溫度、年地區人均產值數據進行標準化處理后并擬合。相比于不考慮參數ω變化的基本Budyko 公式提高了0.24，比考慮ω隨時間變化的Budyko 公式提高0.22 左右。由此說明，引入物理成因變量的ω徑流擬合效果精度更優。以下，采用有一定物理機制變量推求的時變Budyko公式分析該流域各影響因子對徑流變化的影響程度。

表1 多種模型擬合ω的公式比較Tab.1 Results of multiple model and Nash coefficients

3.3 有一定物理機制的時變Budyko徑流模擬結果

結合式（10）～（15）計算各站基準期與變化期對應的彈性系數值，結果見表2。

表2 有一定物理機制的時變參數各期特征值Tab.2 Time-varying parameters derived from the physical mechanism calculate the characteristic values of each period

由表2 可知：流域變化期的降水（P）、徑流值（R）較基準期有所下降，下降速率分別為1.71%和5.13%，流域特征參數（ω）較基準期呈現上升趨勢，增長速率為3.14%，而干燥指數較基準期有所提高，徑流系數則略有減小。其中變化期各影響因子的彈性系數εP、εEo和εω分別為2.06、-0.57 和-0.58，表明降水（P）、潛在蒸散發量（Eo）和流域特征參數（ω）每增加1%時分別引起徑流增加2.06%、減少0.57%和減少0.58%。同時由于彈性系數絕對值的大小反映了徑流對各影響因子的敏感程度，因此可以發現流域徑流變化對降水量最敏感，對流域特征參數次之，對潛在蒸散發量最不敏感。相應的時變參數計算值如圖3所示。

圖3 有一定物理機制的歷年時變參數ωt變化Fig.3 Time-varying parameter calculation value derived from the physical mechanism

3.4 徑流歸因分析

按基準期和變化期分別計算流域特征參數的平均值作為各因子彈性系數的計算值，并根據公式（17）～（19）計算得到較準確的流域各因子對流域徑流變化的貢獻率，結果見表3。

表3 有一定物理機制的時變參數公式得到的各因素對徑流變化的貢獻率Tab.3 The contribution rate of each factor to the change in runoff derived from the time-varying parameter formula including physical causes

結果顯示，相對于基準期（1981-1996），變化期（1997-2012）實測數據徑流深減少了近23 mm。采用模型計算的因降水、潛在蒸散發量和流域特征參數變化引起的徑流深變化量分別為-21.052、-10.747、-11.081 mm，與之相對應的降水、潛在蒸散發量和下墊面變化對徑流影響的貢獻率分別為49%、25%、26%。降水的減少是導致流域徑流變化的主導因素，且其貢獻率接近50%；其次是潛在蒸散發量，但潛在蒸散發量的貢獻率與流域特征參數（人類活動）相當。其原因主要在于該流域面積整體較小，加上隨著地區人口的增加，人均用水以及耕種面積的增加，未來人類活動的影響程度可能會進一步加深。另一方面，就喀斯特氣候變化影響而言，喀斯特地區植被可利用水量較其他地區少，整個流域系統更依賴降水［23］，降水帶來的影響更大，本文的研究與上述結論一致。

4 結 論

（1）六硐河流域1960-2012年年徑流深呈現顯著減小趨勢，其突變點在1997年。年降水量與年徑流深年際變化相一致，表明降水在流域徑流變化中有重要影響。

（2）有一定物理機制的時變Budyko公式模擬徑流的效果要顯著高于其他Budyko公式，表明徑流變化的自然地理因素有一定的非平穩變化特征。

（3）通過對比彈性系數絕對值大小，發現降水量εP最大，表明流域徑流對降水的變化最敏感，其次是下墊面變化εω情況。同時相較于基準期，變化期降水量的彈性系數絕對值逐漸增大，而下墊面變化和潛在蒸散發量的彈性系數絕對值減小。說明隨著時間推移，降水的減少將進一步影響該流域的徑流。

（4）降水量的變化對六硐河流域徑流變化的貢獻率為49%，其次是下墊面的變化，貢獻率為26%，而潛在蒸發量對徑流變化貢獻率較小，因此，可認為徑流量在1997 年后顯著減少的主要原因是降水量的減少。