基于Budyko理论的韩江流域径流变化敏感性分析及归因识别

刘洋, 陈菡, 谭学志*

(1. 中山大学 土木工程学院水资源与环境研究中心，广州510275; 2. 湖北省水利水电规划勘测设计院, 武汉 430064)

0 引言

流域河川径流的归因识别是制定气候变化应对策略和水资源可持续利用的基础[1]。在全球气候变化与人类活动对流域径流的共同作用下，全球诸多流域发生了明显的径流变化[2]。越来越多的学者开始着眼于气候变化与人类活动对流域径流变化影响及归因的研究。其中对径流变化的归因方法有许多。部分研究采用基于Budyko理论的流域水热耦合平衡方程对流域径流变化进行归因的方法，并取得了较好的研究结果[3-8]。该方法与通常采用的分布式流域水文模型相比较，不仅方法简单，且避免了模型参数不确定性对归因结果的影响[8]。

韩江流域河流开发利用较早，流域内修建有多处水利工程，取水引水等人类活动较为频繁。流域内人口稠密，经济发达，经济社会的日益发展对区域水资源的需求日益增加。分析韩江流域径流变化及对其成因的量化对该流域水资源合理开发与可持续利用等方面有着重要意义。在以往的研究中，已有学者对韩江流域及其子流域的河川径流变化进行了分析与研究[9-11]，但关于径流变化归因的量化分析研究相对较少，也缺乏对各影响因子贡献的量化分析。本研究从气候变化和人类活动两方面，基于Budyko理论的水热耦合平衡方程与弹性系数理论，对韩江流域及其子流域即梅江流域、汀江流域的径流变化进行敏感性分析，并对各因子对径流变化造成的贡献进行定量归因分析，旨在为韩江流域在河流开发策略、水灾害防控、水资源合理规划与可持续利用等方面提供科学依据。

1 研究区域与数据来源

1.1 研究区域概况

韩江流域位于粤东和闽西南地区,地理位置处于东经115°13′～117°09′，北纬 23°17′～ 26°05′。韩江流域水系主要由流域上游的梅江子流域与汀江子流域组成，梅江与汀江于大埔三河坝汇合后，形成韩江干流，而后向南延伸，最终于汕头汇入南海。其中，韩江干流全长470 km，总流域面积30 112 km2;梅江子流域面积为13 929 km2;汀江子流域面积为11 802 km2。韩江流域属亚热带季风气候，多年平均气温介于20～21.5 ℃之间，多年平均降雨量约为1 620 mm,气候温和，雨量充沛，流域内大部分地区植被覆盖率高，生态资源丰富[12]。

1.2 数据来源

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

图 1 韩江流域水文站点与气象站点分布 Figure 1 Locations of hydro-climatic stations in Han River Basin

采用韩江流域潮安、溪口、横山3个代表性水文站点的1960—2000年共41年逐日平均流量数据资料。其中，潮安水文站为韩江主要控制站，横山水文站为梅江主要控制站，溪口水文站为汀江主要控制站，站点详细信息见表1。此外，本研究采用韩江流域及其周边区域的7个气象站点1960—2000年共41年降水及气象数据(图1)。其中，用于计算日潜在蒸散发的气象指标包括逐日气压、气温、风速、相对湿度和日照等数据。采用的气象站数据均来自中国气象科学数据共享服务网提供的《中国地面气候资料日值数据集》。

2 研究方法

2.1 径流序列变化突变分析

采用非参数Mann-Kendall突变点检验法与累积距平法来检测年径流系列的突变点，以寻找径流序列变化趋势发生转变的年份。并将两种方法得到的突变点进行对比，互相验证，结合实际情况分析，以确定径流序列变化趋势发生突变的具体年份。其中，(1)非参数Mann-Kendall突变点检验法是一种广泛应用于水文、气象等领域的非参数检验法[13]。(2)累积距平法可以通过观察数据序列差积曲线的变化来判断数据点的离散程度和变化趋势，并找到数据序列的突变位置[13]。在确定突变点后，将径流系列划分为突变点前与突变点后2个子序列，并对突变前后子系列的各水文气象要素进行分析，进而对造成径流突变的各影响因素进行归因分析。

2.2 径流变化敏感性分析与归因识别

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

流域水量平衡方程为：

P=E+R+∆S

(1)

式(1)中：R为径流深，P为降水量，E为实际蒸散发量，ΔS为水量变化量。在较长时间尺度的计算中，ΔS可忽略不计，则式(1)简化为下式：

R=P-E

(2)

(b) Choudhury-Yang公式

Budyko假设[14]认为在一定气候和植被条件下，流域水分和能量之间存在着水热耦合平衡关系，流域多年平均降水量、潜在蒸散发量和实际蒸散发量的关系可以用经验曲线来描述[15]。在此基础上，诸多学者通过数学分析得到了一些衍生公式。其中，Choudhury[16]与Yang[17]等基于经验法或解析法等方法分别得到了形式一致的Choudhury-Yang公式，本研究采用该公式描述流域水量与能量的平衡关系，具体形式如下：

(3)

式(3)中：ET0—长期平均的年潜在蒸散发量，n—反映流域下垫面特征的参数。ET0采用世界粮农组织提出的Penman-Monteith公式[18]计算：

(4)

式(4)中：ET0—可能蒸散量(mm/d)，Δ—饱和水气压曲线斜率(kPa/℃)，Rn—地表净辐射(MJ/(m×d))，G—土壤热通量(MJ/(m2×d))，γ—干湿表常数(kPa/℃)，Tmean—日平均温度(℃)，u2—2 m高处风速(m·s-1)，es—饱和水气压(kPa)，ea—实际水气压(kPa)。

(c) 径流影响因素的弹性系数

(5)

(6)

(7)

在敏感性分析中，为判断各影响因子的弹性系数年际间的变化趋势，本研究采用M-K趋势检验方法来对其趋势进行显著性检验[8]，显著性水平设为α=0.05。

(d) 径流影响因素对径流变化的贡献率

由上述突变点检验方法将径流序列分为突变点前、突变点后两个子序列，记突变点前子序列为R1，突变点后子序列为R2，则△R=R2-R1，降水、潜在蒸散发与下垫面特征指数的突变前后变化同理。将径流的变化量分解为降水变化引起的径流变化量△RCP，潜在蒸散发变化引起的径流变化量△RCET0，下垫面特征指数变化引起的径流变化量△RCn，则

(8)

(9)

(10)

将各影响因子造成的径流变化量相加为△RC=△RCP+△RCET+△RCn，则流域气象要素与下垫面变化对径流变化的贡献率为：

(11)

(12)

(13)

式(11～13)中:CP为降水变化对流域径流变化的贡献率，CET0为潜在蒸散发变化对流域径流变化的贡献率，Cn为下垫面特征指数变化对流域径流变化的贡献率。其中，CP与CET0之和可看作气候变化对流域径流变化的贡献率，Cn可看作人类活动对流域径流变化的贡献率。

3 结果与分析

3.1 径流变化趋势及突变点检验

采用非参数Mann-Kendall突变点检验法对潮安、横山、溪口的流域径流过程发生突变的年份进行识别，得到突变检验结果如图2所示，并以此作为依据分析韩江流域径流突变点发生的年份。综合考虑3个水文站的累积距平趋势，其中3个水文站的累积距平曲线均于1960—1982年呈快速下降趋势，于1982—2000年呈平缓上升趋势。韩江流域3个水文站点均于1961与1973年出现剧烈变化，尤其以溪口站于1973年的起伏最为明显。此外，溪口站还于1990年出现剧烈变化。由累积距平法得到的结果分析，潮安站和横山站径流发生突变的年份为1972、1982年前后，溪口站径流发生突变的年份为1972、1982、1991年前后。由M-K检验法得到的结果分析可知，潮安站径流发生突变的年份为1982、1984、1988年前后，横山站径流发生突变的年份为1975、1979年前后，溪口站径流发生突变的年份为1989、1991年前后。根据3个站点的累积距平曲线与M-K突变检验法结果，以韩江干流的突变时间作为全流域径流发生突变的主要依据，综合实际情况分析，于1960—2000年时间序列中，取1982年为韩江全流域径流发生突变的时间点。

为验证结论的可靠性，与已有的研究进行对比。缪连华[9]等以1951—2010年潮安水文站径流时间序列进行分析，发现其突变点为1980年，与本研究选取1982年为潮安站发生突变年份的结论较为接近，考虑到序列时间长度选取的不同，故结论略有差异，认为可以接受。有关临近流域的突变年份，赖天锃[5]等对珠江流域下东江流域1959—2000年多个水文站径流量进行突变分析，发现其突变时间点为1972年。同时，根据李艳[20]等得出的结论：珠江流域下的北江水系在1980年后的多年平均径流量比1980年以前更大。以上各项研究虽然研究区域与选取的时间序列不尽相同，用于判断突变点发生的方法也不尽相同。但综合来看，珠江流域下各大水系基本于在二十世纪七八十年代年前后发生突变，且突变后流量大致呈上升趋势，故认为本研究取1982年为韩江流域径流突变年份大致合理。

图 2 潮安、横山、溪口水文站年平均径流深累积距平法与M-K法检验结果Figure 2 Cumulative departure curves and M-K test curves of annual runoff depth at three stations

3.2 径流变化影响因素敏感性分析

由公式(5～7)计算韩江流域、梅江流域和汀江流域的径流对降水、潜在蒸散发、下垫面特征指数的弹性系数(1960—2000年)，其结果如图3所示。各影响因子的弹性系数年际波动较大，曲线起伏较为明显。总体而言，径流变化与降水正相关，与潜在蒸散发和下垫面特征指数均呈负相关。降水弹性系数介于1.23～2.22之间，潜在蒸散发弹性系数介于-1.22～-0.23之间，下垫面特征指数弹性系数介于-1.07～-0.18之间。各弹性系数取值表明，降水每增加1%，将导致径流量增加1.23%～2.22%，潜在蒸散发每增加1%，将导致径流量减少0.23%～1.22%，下垫面特征指数每增加1%，将导致径流量减少0.18%～1.07%。

表 2 径流对各影响因子弹性系数的M-K趋势检验结果Table 2 The M-K trend test results of elastic coefficient of runoff on various impact factors

表 3 径流对各影响因子的多年弹性系数Table 3 Multi-year elastic coefficient of runoff on various impact factors

图 3 1960—2000年潮安、横山和溪口水文站实测径流对(a)降水、(b)潜在蒸散发和(c)下垫面特征指数的弹性系数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个站点3种弹性系数可得，梅江子流域的绝对值普遍较高，汀江子流域的绝对值普遍较低，韩江流域则普遍位于二者之间。即梅江流域对各影响因子的敏感性最高，汀江流域对各影响因子的敏感性最低，韩江流域总体的敏感性介于梅江子流域与汀江子流域之间。由表2的M-K趋势检验结果(显著性水平α=0.05)可知，降水弹性系数呈不显著下降趋势，潜在蒸散发与下垫面特征指数的弹性系数呈不显著上升趋势，表明在韩江及其子流域内，降水对径流变化的正相关性有略微减弱的趋势，潜在蒸散发与下垫面特征指数与径流变化的负相关性有略微减弱的趋势，但3种弹性系数的年际波动均较大，总体变化趋势不明显。

由3个流域径流与气候数据的多年平均值，通过公式(5～7)计算得到3个研究区域多年尺度的弹性系数，如表3所示。其结果与年际变化弹性系数的结果分析一致，梅江子流域对3个影响因素的敏感性最高，汀江子流域的敏感性最低，韩江流域处于二者之间。且韩江流域径流变化均与降水呈正相关，与潜在蒸散发与下垫面特征指数均为负相关。

3.3 径流变化归因分析

由公式(11～13)计算出韩江流域、梅江子流域、汀江子流域于1982年前后气候变化与人类活动对径流变化的贡献率(表4)。各影响因子造成的径流变化量之和(△RC=△RCP+△RCET0+△RCn)与实测径流变化量△R差异很小(0.13～0.4 mm)，因此说明本研究采用的归因方法是有效可行的。1960—2000年，气候变化是引起韩江各流域径流变化的主要因素，其中，又以降水变化影响为主，降水变化对韩江流域、梅江子流域和汀江子流域的径流变化贡献率分别为58.24%、52.72%和51.70%，潜在蒸散发变化对韩江流域、梅江子流域和汀江子流域的径流变化贡献率分别为10.80%、5.96%和25.75%，其中汀江子流域的潜在蒸散发变化的贡献较高，其原因是汀江子流域的潜在蒸散发变化较大，达-53.59 mm(表5)。下垫面特征指数变化也对韩江各流域径流变化起到了很大的影响，其中对梅江子流域的影响最大，高达41.32%，对汀江子流域的影响最小，为22.55%，对韩江总流域的径流变化影响为30.96%，其中梅江子流域的下垫面特征指数变化的影响较高，其原因是梅江子流域的下垫面特征指数变化较大，达-0.15(表5)。由此可见，在韩江流域，气候变化是径流变化的主要因素，约占60%～70%，其中大部分贡献由降水变化造成，而人类活动对径流变化的影响也同样重要，约占30%～40%。

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

4 总结与讨论

4.1 气候变化对韩江流域径流变化的影响

韩江流域及其子流域梅江、汀江流域对气候变化均呈现相似响应。降水变化是气候变化对流域径流变化产生影响的关键因素，其变化往往与流域径流量的变化有直接关系。在本研究中，根据敏感性分析结果，降水是韩江流域径流变化的主要影响因子，对径流的贡献呈正相关，且在3个影响因子中其影响程度最高。在本研究中，潜在蒸散发是气候因素中的次要因素，对径流变化的影响相对较小,且与径流呈负相关关系，这与目前许多研究的观点一致[1,4,6-8]，即潜在蒸散发对流域径流变化的弹性系数为负。由表5可知，于1982年前后，韩江流域及其子流域降水量均上涨，潜在蒸散发量则在一定程度上降低，故二者变化皆导致了韩江流域径流量的增加。

4.2 人类活动对韩江流域径流变化的影响

依据表5中n值在两个时间段的变化，韩江流域的n值于1982年前后有所减小，由于敏感性分析的结果，n值与径流变化负相关，因此n值的减小也导致了径流量的增加。根据现有研究，通常认为流域下垫面条件中的地形、地质和土壤条件相对稳定，因此下垫面参数n的变化主要由土地利用变化/植被覆盖变化导致。一般认为n值的增加是由流域内植被覆盖的提高导致的[21]。在本研究中，韩江流域及其子流域的下垫面特征指数均呈现5%～10%左右的下降，可以认为这主要是流域内土地利用变化所致。根据黄苑君的研究[22]，1986—2006年韩江流域未利用地面积减少62 076.0 hm2，建设用地增加48 559.3 hm2，建设用地增加导致地表硬化，流域下垫面不透水面积所占比例增加，致使流域内径流增加。

5 结论

本研究采用累积距平法与M-K突变点检验法分析了韩江流域及其子流域梅江与汀江流域在1960—2000年径流量变化，根据基于Budyko理论的水热耦合平衡方程，结合弹性系数法，对韩江流域干支流对各影响因子的敏感性及其年际变化进行了分析，并从气候变化与人类活动两方面对韩江流域径流变化产生的影响，进行了定量归因分析。主要结论如下：

(1)结合累积距平法与M-K突变点检验法的结果分析，韩江流域干支流于1982年前后发生突变，1982年后韩江流域3个水文站点的多年平均径流量均显著增加。

(2)基于Budyko理论的水热耦合平衡方程，通过弹性系数法对韩江流域干支流进行敏感性分析。在弹性系数的年际变化中，各影响因子的弹性系数年际波动较大，降水弹性系数呈不显著下降趋势，潜在蒸散发与下垫面特征指数的弹性系数呈不显著上升趋势，但总体变化趋势不明显。总体而言，径流变化与降水呈正相关，与潜在蒸散发和下垫面特征指数均呈负相关。

(3)对韩江流域径流变化进行定量归因分析，结果表明：气候变化对韩江流域径流变化的贡献率为59%～77%，其中气候变化的贡献主要体现在降水的变化，对径流变化的贡献率约为50%左右。人类活动对韩江流域径流变化的贡献率为22.6%～41.3%，其中梅江子流域的人类活动贡献率较高，汀江子流域的人类活动贡献率相对较低。