基于Budyko理論的韓江流域徑流變化敏感性分析及歸因識別

劉洋, 陳菡, 譚學(xué)志*

(1. 中山大學(xué) 土木工程學(xué)院水資源與環(huán)境研究中心，廣州510275; 2. 湖北省水利水電規(guī)劃勘測設(shè)計(jì)院, 武漢 430064)

0 引言

流域河川徑流的歸因識別是制定氣候變化應(yīng)對策略和水資源可持續(xù)利用的基礎(chǔ)[1]。在全球氣候變化與人類活動對流域徑流的共同作用下，全球諸多流域發(fā)生了明顯的徑流變化[2]。越來越多的學(xué)者開始著眼于氣候變化與人類活動對流域徑流變化影響及歸因的研究。其中對徑流變化的歸因方法有許多。部分研究采用基于Budyko理論的流域水熱耦合平衡方程對流域徑流變化進(jìn)行歸因的方法，并取得了較好的研究結(jié)果[3-8]。該方法與通常采用的分布式流域水文模型相比較，不僅方法簡單，且避免了模型參數(shù)不確定性對歸因結(jié)果的影響[8]。

韓江流域河流開發(fā)利用較早，流域內(nèi)修建有多處水利工程，取水引水等人類活動較為頻繁。流域內(nèi)人口稠密，經(jīng)濟(jì)發(fā)達(dá)，經(jīng)濟(jì)社會的日益發(fā)展對區(qū)域水資源的需求日益增加。分析韓江流域徑流變化及對其成因的量化對該流域水資源合理開發(fā)與可持續(xù)利用等方面有著重要意義。在以往的研究中，已有學(xué)者對韓江流域及其子流域的河川徑流變化進(jìn)行了分析與研究[9-11]，但關(guān)于徑流變化歸因的量化分析研究相對較少，也缺乏對各影響因子貢獻(xiàn)的量化分析。本研究從氣候變化和人類活動兩方面，基于Budyko理論的水熱耦合平衡方程與彈性系數(shù)理論，對韓江流域及其子流域即梅江流域、汀江流域的徑流變化進(jìn)行敏感性分析，并對各因子對徑流變化造成的貢獻(xiàn)進(jìn)行定量歸因分析，旨在為韓江流域在河流開發(fā)策略、水災(zāi)害防控、水資源合理規(guī)劃與可持續(xù)利用等方面提供科學(xué)依據(jù)。

1 研究區(qū)域與數(shù)據(jù)來源

1.1 研究區(qū)域概況

韓江流域位于粵東和閩西南地區(qū),地理位置處于東經(jīng)115°13′～117°09′，北緯 23°17′～ 26°05′。韓江流域水系主要由流域上游的梅江子流域與汀江子流域組成，梅江與汀江于大埔三河壩匯合后，形成韓江干流，而后向南延伸，最終于汕頭匯入南海。其中，韓江干流全長470 km，總流域面積30 112 km2;梅江子流域面積為13 929 km2;汀江子流域面積為11 802 km2。韓江流域?qū)賮啛釒Ъ撅L(fēng)氣候，多年平均氣溫介于20～21.5 ℃之間，多年平均降雨量約為1 620 mm,氣候溫和，雨量充沛，流域內(nèi)大部分地區(qū)植被覆蓋率高，生態(tài)資源豐富[12]。

1.2 數(shù)據(jù)來源

表 1 韓江流域水文站Table 1 Hydrological stations in Han River Basin

圖 1 韓江流域水文站點(diǎn)與氣象站點(diǎn)分布 Figure 1 Locations of hydro-climatic stations in Han River Basin

采用韓江流域潮安、溪口、橫山3個代表性水文站點(diǎn)的1960—2000年共41年逐日平均流量數(shù)據(jù)資料。其中，潮安水文站為韓江主要控制站，橫山水文站為梅江主要控制站，溪口水文站為汀江主要控制站，站點(diǎn)詳細(xì)信息見表1。此外，本研究采用韓江流域及其周邊區(qū)域的7個氣象站點(diǎn)1960—2000年共41年降水及氣象數(shù)據(jù)(圖1)。其中，用于計(jì)算日潛在蒸散發(fā)的氣象指標(biāo)包括逐日氣壓、氣溫、風(fēng)速、相對濕度和日照等數(shù)據(jù)。采用的氣象站數(shù)據(jù)均來自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)提供的《中國地面氣候資料日值數(shù)據(jù)集》。

2 研究方法

2.1 徑流序列變化突變分析

采用非參數(shù)Mann-Kendall突變點(diǎn)檢驗(yàn)法與累積距平法來檢測年徑流系列的突變點(diǎn)，以尋找徑流序列變化趨勢發(fā)生轉(zhuǎn)變的年份。并將兩種方法得到的突變點(diǎn)進(jìn)行對比，互相驗(yàn)證，結(jié)合實(shí)際情況分析，以確定徑流序列變化趨勢發(fā)生突變的具體年份。其中，(1)非參數(shù)Mann-Kendall突變點(diǎn)檢驗(yàn)法是一種廣泛應(yīng)用于水文、氣象等領(lǐng)域的非參數(shù)檢驗(yàn)法[13]。(2)累積距平法可以通過觀察數(shù)據(jù)序列差積曲線的變化來判斷數(shù)據(jù)點(diǎn)的離散程度和變化趨勢，并找到數(shù)據(jù)序列的突變位置[13]。在確定突變點(diǎn)后，將徑流系列劃分為突變點(diǎn)前與突變點(diǎn)后2個子序列，并對突變前后子系列的各水文氣象要素進(jìn)行分析，進(jìn)而對造成徑流突變的各影響因素進(jìn)行歸因分析。

2.2 徑流變化敏感性分析與歸因識別

(a) 多年尺度下的流域水量平衡方程

流域水量平衡方程為：

P=E+R+?S

(1)

式(1)中：R為徑流深，P為降水量，E為實(shí)際蒸散發(fā)量，ΔS為水量變化量。在較長時(shí)間尺度的計(jì)算中，ΔS可忽略不計(jì)，則式(1)簡化為下式：

R=P-E

(2)

(b) Choudhury-Yang公式

Budyko假設(shè)[14]認(rèn)為在一定氣候和植被條件下，流域水分和能量之間存在著水熱耦合平衡關(guān)系，流域多年平均降水量、潛在蒸散發(fā)量和實(shí)際蒸散發(fā)量的關(guān)系可以用經(jīng)驗(yàn)曲線來描述[15]。在此基礎(chǔ)上，諸多學(xué)者通過數(shù)學(xué)分析得到了一些衍生公式。其中，Choudhury[16]與Yang[17]等基于經(jīng)驗(yàn)法或解析法等方法分別得到了形式一致的Choudhury-Yang公式，本研究采用該公式描述流域水量與能量的平衡關(guān)系，具體形式如下：

(3)

式(3)中：ET0—長期平均的年潛在蒸散發(fā)量，n—反映流域下墊面特征的參數(shù)。ET0采用世界糧農(nóng)組織提出的Penman-Monteith公式[18]計(jì)算：

(4)

式(4)中：ET0—可能蒸散量(mm/d)，Δ—飽和水氣壓曲線斜率(kPa/℃)，Rn—地表凈輻射(MJ/(m×d))，G—土壤熱通量(MJ/(m2×d))，γ—干濕表常數(shù)(kPa/℃)，Tmean—日平均溫度(℃)，u2—2 m高處風(fēng)速(m·s-1)，es—飽和水氣壓(kPa)，ea—實(shí)際水氣壓(kPa)。

(c) 徑流影響因素的彈性系數(shù)

(5)

(6)

(7)

在敏感性分析中，為判斷各影響因子的彈性系數(shù)年際間的變化趨勢，本研究采用M-K趨勢檢驗(yàn)方法來對其趨勢進(jìn)行顯著性檢驗(yàn)[8]，顯著性水平設(shè)為α=0.05。

(d) 徑流影響因素對徑流變化的貢獻(xiàn)率

由上述突變點(diǎn)檢驗(yàn)方法將徑流序列分為突變點(diǎn)前、突變點(diǎn)后兩個子序列，記突變點(diǎn)前子序列為R1，突變點(diǎn)后子序列為R2，則△R=R2-R1，降水、潛在蒸散發(fā)與下墊面特征指數(shù)的突變前后變化同理。將徑流的變化量分解為降水變化引起的徑流變化量△RCP，潛在蒸散發(fā)變化引起的徑流變化量△RCET0，下墊面特征指數(shù)變化引起的徑流變化量△RCn，則

(8)

(9)

(10)

將各影響因子造成的徑流變化量相加為△RC=△RCP+△RCET+△RCn，則流域氣象要素與下墊面變化對徑流變化的貢獻(xiàn)率為：

(11)

(12)

(13)

式(11～13)中:CP為降水變化對流域徑流變化的貢獻(xiàn)率，CET0為潛在蒸散發(fā)變化對流域徑流變化的貢獻(xiàn)率，Cn為下墊面特征指數(shù)變化對流域徑流變化的貢獻(xiàn)率。其中，CP與CET0之和可看作氣候變化對流域徑流變化的貢獻(xiàn)率，Cn可看作人類活動對流域徑流變化的貢獻(xiàn)率。

3 結(jié)果與分析

3.1 徑流變化趨勢及突變點(diǎn)檢驗(yàn)

采用非參數(shù)Mann-Kendall突變點(diǎn)檢驗(yàn)法對潮安、橫山、溪口的流域徑流過程發(fā)生突變的年份進(jìn)行識別，得到突變檢驗(yàn)結(jié)果如圖2所示，并以此作為依據(jù)分析韓江流域徑流突變點(diǎn)發(fā)生的年份。綜合考慮3個水文站的累積距平趨勢，其中3個水文站的累積距平曲線均于1960—1982年呈快速下降趨勢，于1982—2000年呈平緩上升趨勢。韓江流域3個水文站點(diǎn)均于1961與1973年出現(xiàn)劇烈變化，尤其以溪口站于1973年的起伏最為明顯。此外，溪口站還于1990年出現(xiàn)劇烈變化。由累積距平法得到的結(jié)果分析，潮安站和橫山站徑流發(fā)生突變的年份為1972、1982年前后，溪口站徑流發(fā)生突變的年份為1972、1982、1991年前后。由M-K檢驗(yàn)法得到的結(jié)果分析可知，潮安站徑流發(fā)生突變的年份為1982、1984、1988年前后，橫山站徑流發(fā)生突變的年份為1975、1979年前后，溪口站徑流發(fā)生突變的年份為1989、1991年前后。根據(jù)3個站點(diǎn)的累積距平曲線與M-K突變檢驗(yàn)法結(jié)果，以韓江干流的突變時(shí)間作為全流域徑流發(fā)生突變的主要依據(jù)，綜合實(shí)際情況分析，于1960—2000年時(shí)間序列中，取1982年為韓江全流域徑流發(fā)生突變的時(shí)間點(diǎn)。

為驗(yàn)證結(jié)論的可靠性，與已有的研究進(jìn)行對比?？娺B華[9]等以1951—2010年潮安水文站徑流時(shí)間序列進(jìn)行分析，發(fā)現(xiàn)其突變點(diǎn)為1980年，與本研究選取1982年為潮安站發(fā)生突變年份的結(jié)論較為接近，考慮到序列時(shí)間長度選取的不同，故結(jié)論略有差異，認(rèn)為可以接受。有關(guān)臨近流域的突變年份，賴天锃[5]等對珠江流域下東江流域1959—2000年多個水文站徑流量進(jìn)行突變分析，發(fā)現(xiàn)其突變時(shí)間點(diǎn)為1972年。同時(shí)，根據(jù)李艷[20]等得出的結(jié)論：珠江流域下的北江水系在1980年后的多年平均徑流量比1980年以前更大。以上各項(xiàng)研究雖然研究區(qū)域與選取的時(shí)間序列不盡相同，用于判斷突變點(diǎn)發(fā)生的方法也不盡相同。但綜合來看，珠江流域下各大水系基本于在二十世紀(jì)七八十年代年前后發(fā)生突變，且突變后流量大致呈上升趨勢，故認(rèn)為本研究取1982年為韓江流域徑流突變年份大致合理。

圖 2 潮安、橫山、溪口水文站年平均徑流深累積距平法與M-K法檢驗(yàn)結(jié)果Figure 2 Cumulative departure curves and M-K test curves of annual runoff depth at three stations

3.2 徑流變化影響因素敏感性分析

由公式(5～7)計(jì)算韓江流域、梅江流域和汀江流域的徑流對降水、潛在蒸散發(fā)、下墊面特征指數(shù)的彈性系數(shù)(1960—2000年)，其結(jié)果如圖3所示。各影響因子的彈性系數(shù)年際波動較大，曲線起伏較為明顯?？傮w而言，徑流變化與降水正相關(guān)，與潛在蒸散發(fā)和下墊面特征指數(shù)均呈負(fù)相關(guān)。降水彈性系數(shù)介于1.23～2.22之間，潛在蒸散發(fā)彈性系數(shù)介于-1.22～-0.23之間，下墊面特征指數(shù)彈性系數(shù)介于-1.07～-0.18之間。各彈性系數(shù)取值表明，降水每增加1%，將導(dǎo)致徑流量增加1.23%～2.22%，潛在蒸散發(fā)每增加1%，將導(dǎo)致徑流量減少0.23%～1.22%，下墊面特征指數(shù)每增加1%，將導(dǎo)致徑流量減少0.18%～1.07%。

表 2 徑流對各影響因子彈性系數(shù)的M-K趨勢檢驗(yàn)結(jié)果Table 2 The M-K trend test results of elastic coefficient of runoff on various impact factors

表 3 徑流對各影響因子的多年彈性系數(shù)Table 3 Multi-year elastic coefficient of runoff on various impact factors

圖 3 1960—2000年潮安、橫山和溪口水文站實(shí)測徑流對(a)降水、(b)潛在蒸散發(fā)和(c)下墊面特征指數(shù)的彈性系數(shù)Figure 3 Elastic coefficient of runoff depth on precipitation (a), potential evapotranspiration (b) and underlying surface index (c) at three stations during 1960—2000

表 4 韓江流域徑流變化歸因分析Table 4 Attribution analysis of runoff changes in Han River Basin

對比3個站點(diǎn)3種彈性系數(shù)可得，梅江子流域的絕對值普遍較高，汀江子流域的絕對值普遍較低，韓江流域則普遍位于二者之間。即梅江流域?qū)Ω饔绊懸蜃拥拿舾行宰罡撸〗饔驅(qū)Ω饔绊懸蜃拥拿舾行宰畹?，韓江流域總體的敏感性介于梅江子流域與汀江子流域之間。由表2的M-K趨勢檢驗(yàn)結(jié)果(顯著性水平α=0.05)可知，降水彈性系數(shù)呈不顯著下降趨勢，潛在蒸散發(fā)與下墊面特征指數(shù)的彈性系數(shù)呈不顯著上升趨勢，表明在韓江及其子流域內(nèi)，降水對徑流變化的正相關(guān)性有略微減弱的趨勢，潛在蒸散發(fā)與下墊面特征指數(shù)與徑流變化的負(fù)相關(guān)性有略微減弱的趨勢，但3種彈性系數(shù)的年際波動均較大，總體變化趨勢不明顯。

由3個流域徑流與氣候數(shù)據(jù)的多年平均值，通過公式(5～7)計(jì)算得到3個研究區(qū)域多年尺度的彈性系數(shù)，如表3所示。其結(jié)果與年際變化彈性系數(shù)的結(jié)果分析一致，梅江子流域?qū)?個影響因素的敏感性最高，汀江子流域的敏感性最低，韓江流域處于二者之間。且韓江流域徑流變化均與降水呈正相關(guān)，與潛在蒸散發(fā)與下墊面特征指數(shù)均為負(fù)相關(guān)。

3.3 徑流變化歸因分析

由公式(11～13)計(jì)算出韓江流域、梅江子流域、汀江子流域于1982年前后氣候變化與人類活動對徑流變化的貢獻(xiàn)率(表4)。各影響因子造成的徑流變化量之和(△RC=△RCP+△RCET0+△RCn)與實(shí)測徑流變化量△R差異很小(0.13～0.4 mm)，因此說明本研究采用的歸因方法是有效可行的。1960—2000年，氣候變化是引起韓江各流域徑流變化的主要因素，其中，又以降水變化影響為主，降水變化對韓江流域、梅江子流域和汀江子流域的徑流變化貢獻(xiàn)率分別為58.24%、52.72%和51.70%，潛在蒸散發(fā)變化對韓江流域、梅江子流域和汀江子流域的徑流變化貢獻(xiàn)率分別為10.80%、5.96%和25.75%，其中汀江子流域的潛在蒸散發(fā)變化的貢獻(xiàn)較高，其原因是汀江子流域的潛在蒸散發(fā)變化較大，達(dá)-53.59 mm(表5)。下墊面特征指數(shù)變化也對韓江各流域徑流變化起到了很大的影響，其中對梅江子流域的影響最大，高達(dá)41.32%，對汀江子流域的影響最小，為22.55%，對韓江總流域的徑流變化影響為30.96%，其中梅江子流域的下墊面特征指數(shù)變化的影響較高，其原因是梅江子流域的下墊面特征指數(shù)變化較大，達(dá)-0.15(表5)。由此可見，在韓江流域，氣候變化是徑流變化的主要因素，約占60%～70%，其中大部分貢獻(xiàn)由降水變化造成，而人類活動對徑流變化的影響也同樣重要，約占30%～40%。

表 5 韓江流域各影響要素1982年前后變化Table 5 Changes in the impact factors of the runoff changes in Han River Basin before and after 1982

4 總結(jié)與討論

4.1 氣候變化對韓江流域徑流變化的影響

韓江流域及其子流域梅江、汀江流域?qū)夂蜃兓尸F(xiàn)相似響應(yīng)。降水變化是氣候變化對流域徑流變化產(chǎn)生影響的關(guān)鍵因素，其變化往往與流域徑流量的變化有直接關(guān)系。在本研究中，根據(jù)敏感性分析結(jié)果，降水是韓江流域徑流變化的主要影響因子，對徑流的貢獻(xiàn)呈正相關(guān)，且在3個影響因子中其影響程度最高。在本研究中，潛在蒸散發(fā)是氣候因素中的次要因素，對徑流變化的影響相對較小,且與徑流呈負(fù)相關(guān)關(guān)系，這與目前許多研究的觀點(diǎn)一致[1,4,6-8]，即潛在蒸散發(fā)對流域徑流變化的彈性系數(shù)為負(fù)。由表5可知，于1982年前后，韓江流域及其子流域降水量均上漲，潛在蒸散發(fā)量則在一定程度上降低，故二者變化皆導(dǎo)致了韓江流域徑流量的增加。

4.2 人類活動對韓江流域徑流變化的影響

依據(jù)表5中n值在兩個時(shí)間段的變化，韓江流域的n值于1982年前后有所減小，由于敏感性分析的結(jié)果，n值與徑流變化負(fù)相關(guān)，因此n值的減小也導(dǎo)致了徑流量的增加。根據(jù)現(xiàn)有研究，通常認(rèn)為流域下墊面條件中的地形、地質(zhì)和土壤條件相對穩(wěn)定，因此下墊面參數(shù)n的變化主要由土地利用變化/植被覆蓋變化導(dǎo)致。一般認(rèn)為n值的增加是由流域內(nèi)植被覆蓋的提高導(dǎo)致的[21]。在本研究中，韓江流域及其子流域的下墊面特征指數(shù)均呈現(xiàn)5%～10%左右的下降，可以認(rèn)為這主要是流域內(nèi)土地利用變化所致。根據(jù)黃苑君的研究[22]，1986—2006年韓江流域未利用地面積減少62 076.0 hm2，建設(shè)用地增加48 559.3 hm2，建設(shè)用地增加導(dǎo)致地表硬化，流域下墊面不透水面積所占比例增加，致使流域內(nèi)徑流增加。

5 結(jié)論

本研究采用累積距平法與M-K突變點(diǎn)檢驗(yàn)法分析了韓江流域及其子流域梅江與汀江流域在1960—2000年徑流量變化，根據(jù)基于Budyko理論的水熱耦合平衡方程，結(jié)合彈性系數(shù)法，對韓江流域干支流對各影響因子的敏感性及其年際變化進(jìn)行了分析，并從氣候變化與人類活動兩方面對韓江流域徑流變化產(chǎn)生的影響，進(jìn)行了定量歸因分析。主要結(jié)論如下：

(1)結(jié)合累積距平法與M-K突變點(diǎn)檢驗(yàn)法的結(jié)果分析，韓江流域干支流于1982年前后發(fā)生突變，1982年后韓江流域3個水文站點(diǎn)的多年平均徑流量均顯著增加。

(2)基于Budyko理論的水熱耦合平衡方程，通過彈性系數(shù)法對韓江流域干支流進(jìn)行敏感性分析。在彈性系數(shù)的年際變化中，各影響因子的彈性系數(shù)年際波動較大，降水彈性系數(shù)呈不顯著下降趨勢，潛在蒸散發(fā)與下墊面特征指數(shù)的彈性系數(shù)呈不顯著上升趨勢，但總體變化趨勢不明顯?？傮w而言，徑流變化與降水呈正相關(guān)，與潛在蒸散發(fā)和下墊面特征指數(shù)均呈負(fù)相關(guān)。

(3)對韓江流域徑流變化進(jìn)行定量歸因分析，結(jié)果表明：氣候變化對韓江流域徑流變化的貢獻(xiàn)率為59%～77%，其中氣候變化的貢獻(xiàn)主要體現(xiàn)在降水的變化，對徑流變化的貢獻(xiàn)率約為50%左右。人類活動對韓江流域徑流變化的貢獻(xiàn)率為22.6%～41.3%，其中梅江子流域的人類活動貢獻(xiàn)率較高，汀江子流域的人類活動貢獻(xiàn)率相對較低。