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

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

0 引言

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

韓江流域河流開發(fā)利用較早，流域內(nèi)修建有多處水利工程，取水引水等人類活動(dòng)較為頻繁。流域內(nèi)人口稠密，經(jīng)濟(jì)發(fā)達(dá)，經(jīng)濟(jì)社會(huì)的日益發(fā)展對(duì)區(qū)域水資源的需求日益增加。分析韓江流域徑流變化及對(duì)其成因的量化對(duì)該流域水資源合理開發(fā)與可持續(xù)利用等方面有著重要意義。在以往的研究中，已有學(xué)者對(duì)韓江流域及其子流域的河川徑流變化進(jìn)行了分析與研究[9-11]，但關(guān)于徑流變化歸因的量化分析研究相對(duì)較少，也缺乏對(duì)各影響因子貢獻(xiàn)的量化分析。本研究從氣候變化和人類活動(dòng)兩方面，基于Budyko理論的水熱耦合平衡方程與彈性系數(shù)理論，對(duì)韓江流域及其子流域即梅江流域、汀江流域的徑流變化進(jìn)行敏感性分析，并對(duì)各因子對(duì)徑流變化造成的貢獻(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′。韓江流域水系主要由流域上游的梅江子流域與汀江子流域組成，梅江與汀江于大埔三河壩匯合后，形成韓江干流，而后向南延伸，最終于汕頭匯入南海。其中，韓江干流全長(zhǎng)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個(gè)代表性水文站點(diǎn)的1960—2000年共41年逐日平均流量數(shù)據(jù)資料。其中，潮安水文站為韓江主要控制站，橫山水文站為梅江主要控制站，溪口水文站為汀江主要控制站，站點(diǎn)詳細(xì)信息見表1。此外，本研究采用韓江流域及其周邊區(qū)域的7個(gè)氣象站點(diǎn)1960—2000年共41年降水及氣象數(shù)據(jù)(圖1)。其中，用于計(jì)算日潛在蒸散發(fā)的氣象指標(biāo)包括逐日氣壓、氣溫、風(fēng)速、相對(duì)濕度和日照等數(shù)據(jù)。采用的氣象站數(shù)據(jù)均來自中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)提供的《中國(guó)地面氣候資料日值數(shù)據(jù)集》。

2 研究方法

2.1 徑流序列變化突變分析

采用非參數(shù)Mann-Kendall突變點(diǎn)檢驗(yàn)法與累積距平法來檢測(cè)年徑流系列的突變點(diǎn)，以尋找徑流序列變化趨勢(shì)發(fā)生轉(zhuǎn)變的年份。并將兩種方法得到的突變點(diǎn)進(jìn)行對(duì)比，互相驗(yàn)證，結(jié)合實(shí)際情況分析，以確定徑流序列變化趨勢(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ì)，并找到數(shù)據(jù)序列的突變位置[13]。在確定突變點(diǎn)后，將徑流系列劃分為突變點(diǎn)前與突變點(diǎn)后2個(gè)子序列，并對(duì)突變前后子系列的各水文氣象要素進(jìn)行分析，進(jìn)而對(duì)造成徑流突變的各影響因素進(jìn)行歸因分析。

2.2 徑流變化敏感性分析與歸因識(shí)別

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

流域水量平衡方程為：

P=E+R+?S

(1)

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

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—長(zhǎng)期平均的年潛在蒸散發(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ù)年際間的變化趨勢(shì)，本研究采用M-K趨勢(shì)檢驗(yàn)方法來對(duì)其趨勢(shì)進(jìn)行顯著性檢驗(yàn)[8]，顯著性水平設(shè)為α=0.05。

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

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

(8)

(9)

(10)

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

(11)

(12)

(13)

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

3 結(jié)果與分析

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

采用非參數(shù)Mann-Kendall突變點(diǎn)檢驗(yàn)法對(duì)潮安、橫山、溪口的流域徑流過程發(fā)生突變的年份進(jìn)行識(shí)別，得到突變檢驗(yàn)結(jié)果如圖2所示，并以此作為依據(jù)分析韓江流域徑流突變點(diǎn)發(fā)生的年份。綜合考慮3個(gè)水文站的累積距平趨勢(shì)，其中3個(gè)水文站的累積距平曲線均于1960—1982年呈快速下降趨勢(shì)，于1982—2000年呈平緩上升趨勢(shì)。韓江流域3個(gè)水文站點(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個(gè)站點(diǎn)的累積距平曲線與M-K突變檢驗(yàn)法結(jié)果，以韓江干流的突變時(shí)間作為全流域徑流發(fā)生突變的主要依據(jù)，綜合實(shí)際情況分析，于1960—2000年時(shí)間序列中，取1982年為韓江全流域徑流發(fā)生突變的時(shí)間點(diǎn)。

為驗(yàn)證結(jié)論的可靠性，與已有的研究進(jìn)行對(duì)比。繆連華[9]等以1951—2010年潮安水文站徑流時(shí)間序列進(jìn)行分析，發(fā)現(xiàn)其突變點(diǎn)為1980年，與本研究選取1982年為潮安站發(fā)生突變年份的結(jié)論較為接近，考慮到序列時(shí)間長(zhǎng)度選取的不同，故結(jié)論略有差異，認(rèn)為可以接受。有關(guān)臨近流域的突變年份，賴天锃[5]等對(duì)珠江流域下東江流域1959—2000年多個(gè)水文站徑流量進(jìn)行突變分析，發(fā)現(xiàn)其突變時(shí)間點(diǎn)為1972年。同時(shí)，根據(jù)李艷[20]等得出的結(jié)論：珠江流域下的北江水系在1980年后的多年平均徑流量比1980年以前更大。以上各項(xiàng)研究雖然研究區(qū)域與選取的時(shí)間序列不盡相同，用于判斷突變點(diǎn)發(fā)生的方法也不盡相同。但綜合來看，珠江流域下各大水系基本于在二十世紀(jì)七八十年代年前后發(fā)生突變，且突變后流量大致呈上升趨勢(shì)，故認(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ì)算韓江流域、梅江流域和汀江流域的徑流對(duì)降水、潛在蒸散發(fā)、下墊面特征指數(shù)的彈性系數(shù)(1960—2000年)，其結(jié)果如圖3所示。各影響因子的彈性系數(shù)年際波動(dòng)較大，曲線起伏較為明顯?？傮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 徑流對(duì)各影響因子彈性系數(shù)的M-K趨勢(shì)檢驗(yàn)結(jié)果Table 2 The M-K trend test results of elastic coefficient of runoff on various impact factors

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

圖 3 1960—2000年潮安、橫山和溪口水文站實(shí)測(cè)徑流對(duì)(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

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

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

3.3 徑流變化歸因分析

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

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

4.2 人類活動(dòng)對(duì)韓江流域徑流變化的影響

依據(jù)表5中n值在兩個(gè)時(shí)間段的變化，韓江流域的n值于1982年前后有所減小，由于敏感性分析的結(jié)果，n值與徑流變化負(fù)相關(guān)，因此n值的減小也導(dǎo)致了徑流量的增加。根據(jù)現(xiàn)有研究，通常認(rèn)為流域下墊面條件中的地形、地質(zhì)和土壤條件相對(duì)穩(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ù)法，對(duì)韓江流域干支流對(duì)各影響因子的敏感性及其年際變化進(jìn)行了分析，并從氣候變化與人類活動(dòng)兩方面對(duì)韓江流域徑流變化產(chǎn)生的影響，進(jìn)行了定量歸因分析。主要結(jié)論如下：

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

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

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