王明陽(yáng)，董前進(jìn)，張艷敏，黃 馗

（1. 武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072； 2. 中工武大設(shè)計(jì)集團(tuán)有限公司，湖北 武漢 430070；3. 廣西電網(wǎng)電力調(diào)度控制中心，廣西 南寧 530013）

0 引 言

近年來(lái)，全球變暖趨勢(shì)加劇，打破了持續(xù)穩(wěn)定的水循環(huán)系統(tǒng)，加劇了水資源分布不平衡現(xiàn)狀［1］。同時(shí)，人類活動(dòng)對(duì)流域產(chǎn)匯流等水資源循環(huán)過(guò)程產(chǎn)生了重大影響［2］。河川徑流是水循環(huán)過(guò)程的重要環(huán)節(jié)，它反映了集水區(qū)在時(shí)間和空間上的綜合變化［3］。同時(shí)，作為生態(tài)、經(jīng)濟(jì)和社會(huì)發(fā)展的基本需求，河川徑流的變化給人類社會(huì)的生產(chǎn)、生活、生態(tài)帶來(lái)了極大影響。因此定量分析河川徑流的變化及驅(qū)動(dòng)機(jī)制是水資源管理中的重要基礎(chǔ)工作，可以從理論上幫助我們適應(yīng)未來(lái)氣候條件的演化，拓展非一致性研究成果［4］。

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

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

由此，基于Budyko假設(shè)對(duì)六硐河流域進(jìn)行年徑流變化歸因分析。根據(jù)六硐河流域數(shù)據(jù)獲取情況及水文情勢(shì)變化產(chǎn)生的物理成因，分別建立由一定物理機(jī)制推求的和由時(shí)間推求的兩類時(shí)變Budyko 公式，以Nash 效率系數(shù)為評(píng)價(jià)指標(biāo)，驗(yàn)證最適應(yīng)研究區(qū)的水熱耦合公式，該研究將有助于認(rèn)識(shí)貴州省黔南地區(qū)流域徑流變化機(jī)理。

1 研究區(qū)概況

六硐河為珠江水系西江干流紅水河段支流，位于東經(jīng)106°54'～107°38'，北緯25°11'～26°21'。流域自北向南流經(jīng)貴州黔南獨(dú)山縣、平塘縣、都勻等市，其干流發(fā)源于獨(dú)山縣拉林鄉(xiāng)八壩，流出貴州省于三堡匯合曹渡河，最終匯入紅水河。六硐河全長(zhǎng)164 km，研究區(qū)總面積1 441 km2。流域?qū)儆趤啛釒Ъ撅L(fēng)性濕潤(rùn)氣候，夏季高溫多雨，冬季溫和少雨。年均氣溫15.0～17.0 ℃，流域年降水量1 100～1 300 mm，河口多年平均流量56 m3/s［12］。

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

2 研究方法

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

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

2.1 Budyko假設(shè)

在研究全球水量和能量平衡分析時(shí)，Budyko 認(rèn)為地表年尺度上蒸散發(fā)量由當(dāng)?shù)亟邓a(bǔ)給和當(dāng)?shù)剌椛淠芰抗┙o的平衡關(guān)系決定，進(jìn)而提出Budyko理論［15］，公式表示為：

式中：φ定義為某地區(qū)的干燥指數(shù)，其值為；Bk（）為Budyko理論下的函數(shù)形式，其中廣泛應(yīng)用的Choudhury-Yang公式［17］表達(dá)如下：

式中相關(guān)符號(hào)含義同公式（2），ω為流域特征參數(shù)，表示流域下墊面特征的綜合參數(shù)。

2.2 時(shí)變Budyko公式

原始Budyko公式是不含參數(shù)的［17］，現(xiàn)有研究表明引入特定參數(shù)得到改進(jìn)的Budyko公式更能反映流域?qū)嶋H情況，且模擬性能更佳。時(shí)變Budyko 公式是對(duì)公式（4）的一種改進(jìn)。在時(shí)變Budyko 公式中，流域特征參數(shù)ω不再為定值，而是隨時(shí)間或者其他解釋變量變化的量。

構(gòu)建的時(shí)變Budyko公式可用如下公式表示：

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

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

解釋變量x表示流域降水、輻射能量、溫度、人均取用水量等有一定物理含義的變量，使得在該公式下的流域特征參數(shù)ω與地區(qū)人類活動(dòng)有一定的聯(lián)系，β為相關(guān)系數(shù)矩陣。結(jié)合1981-2012年國(guó)家統(tǒng)計(jì)年報(bào)，選擇年平均降水、年平均溫度數(shù)據(jù)、地區(qū)人均生產(chǎn)總值這3個(gè)變量作為有一定物理機(jī)制時(shí)變參數(shù)的解釋變量。

隨時(shí)間變化的Budyko 公式一般以線性和二次多項(xiàng)式函數(shù)為主如下式所示：

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

2.3 徑流變化歸因分析

在長(zhǎng)時(shí)間尺度（如年）徑流變化歸因分析上通?？珊雎粤饔蛐钏康淖兓?，由此結(jié)合水量平衡方程和公式（4），可得到徑流變化的表達(dá)式如下［18，19］：

其中各項(xiàng)的偏微分形式如下：

為統(tǒng)一彈性系數(shù)公式，以干燥指數(shù)φ對(duì)式（10）進(jìn)行簡(jiǎn)化，有：

在已有流域潛在蒸散發(fā)量和降水的前提下，由時(shí)變Budyko公式模擬徑流的表達(dá)式為：

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

由式（16）并通過(guò)離差平方和最小可得到最優(yōu)的流域特征參數(shù)值ω。

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

其中ΔR代表兩個(gè)時(shí)期徑流變化的差值；ΔRp、ΔREo和ΔRω分別代表P、Eo和ω變化導(dǎo)致的徑流變化，為簡(jiǎn)便表達(dá)下文均以xi替代。

則流域各影響因素對(duì)流域徑流變化的貢獻(xiàn)率η為［21］：

3 結(jié)果與討論

3.1 六硐河流域水文氣象要素趨勢(shì)變化與突變檢驗(yàn)

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

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

3.2 多種Budyko公式擬合參數(shù)結(jié)果與性能對(duì)比

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

由擬合結(jié)果數(shù)據(jù)來(lái)看，有一定物理機(jī)制變量推求的時(shí)變Budyko 公式估算流域徑流深的Nash 系數(shù)最高，為0.923。該模型將年降水、年均溫度、年地區(qū)人均產(chǎn)值數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理后并擬合。相比于不考慮參數(shù)ω變化的基本Budyko 公式提高了0.24，比考慮ω隨時(shí)間變化的Budyko 公式提高0.22 左右。由此說(shuō)明，引入物理成因變量的ω徑流擬合效果精度更優(yōu)。以下，采用有一定物理機(jī)制變量推求的時(shí)變Budyko公式分析該流域各影響因子對(duì)徑流變化的影響程度。

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

3.3 有一定物理機(jī)制的時(shí)變Budyko徑流模擬結(jié)果

結(jié)合式（10）～（15）計(jì)算各站基準(zhǔn)期與變化期對(duì)應(yīng)的彈性系數(shù)值，結(jié)果見(jiàn)表2。

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

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

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

3.4 徑流歸因分析

按基準(zhǔn)期和變化期分別計(jì)算流域特征參數(shù)的平均值作為各因子彈性系數(shù)的計(jì)算值，并根據(jù)公式（17）～（19）計(jì)算得到較準(zhǔn)確的流域各因子對(duì)流域徑流變化的貢獻(xiàn)率，結(jié)果見(jiàn)表3。

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

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

4 結(jié) 論

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

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

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

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