基于Budyko弹性系数法的佳芦河流域径流变化归因识别

陈丽丽，莫淑红，巩 瑶

(1.西安理工大学 省部共建西北旱区生态水利国家重点实验室，陕西 西安 710048；2.甘肃省水利水电勘测设计研究院有限责任公司，甘肃 兰州 730000)

1 研究背景

流域水循环受社会水循环和自然水循环耦合影响[1]，气候变化和人类活动共同作用驱动河川径流变化，其中气候条件主要通过降水、气温、风速、湿度等气象要素的变化影响径流过程[2]；人类活动的影响包括水利工程建设运行、土地利用方式变化、工农业用水、水土保持措施等引起的径流变化[3]。由于区域气候变化和人类活动综合作用于径流形成过程，因此径流动态变化及其归因识别成为水文水资源领域的研究热点问题之一[4-7]。

许多学者开展大量研究以尝试定量区分气候变化和人类活动对径流的影响，目前主要采用的方法有统计分析法、水文模型模拟方法和基于Budyko理论的方法等。如王随继等[8]提出累积量斜率变化率比较方法，并用于黄河一级支流皇甫川流域；陈玫君等[9]采用双累积曲线法和线性回归模型量化北洛河上游降水变化和人类活动对年径流、汛期径流的影响程度；Chiew等[10]采用SIMHYD分析了澳大利亚南部气候变化对径流的影响，结果表明该模型适用于估算气候变化对径流的影响；刘柏君等[11]运用SWAT模型研究发现，人类活动是海流兔河流域基流变化的主要原因；陈宏等[12]采用VIC模型分析栾河流域径流变化的驱动因素，指出气候变化和人类活动对径流影响的贡献率分别为54.9%和45.1%。利用Budyko假设理论可以更为简便有效地探究流域径流变化的驱动机制，杨大文等[13]和孙福宝等[14]验证了Budyko水热耦合平衡方程在黄河流域的适用性；Zheng等[15]利用径流与潜在蒸散发的局部弹性和潜在蒸散发之间的Budyko互补关系，定量划分了气候和集水区变化对黄河流域径流变化的影响；Jiang等[16]研究发现，Budyko方程中参数w与气候变化和人类活动因素均有显著的相关性，气候变化和人类活动是渭河径流减少的重要驱动因子；张丽梅等[17]分析了渭河流域多个水文站近60年来的径流序列，基于Budyko弹性系数法得出人类活动引起的下垫面变化是径流减少的主要因素。

佳芦河流域地处黄土高原，属于典型的干旱半干旱区域，对气候变化敏感，是生态环境脆弱区[18]。20世纪70年代以来，退耕还林还草、土地荒漠化防治、水利工程建设等一系列人类活动不断改变着流域下垫面条件，影响着其水文过程和径流情势，为此分析和研究该流域水文水资源演变特性，厘清气候变化和人类活动对该流域径流变化分项贡献很有必要。本文基于Budyko假设的弹性系数法展开对佳芦河径流变化的归因分析，以期揭示该区域水文过程对环境变化的响应关系，为其水资源的开发利用和管理提供参考。

2 数据来源与研究方法

2.1 研究区概况

佳芦河流域地处黄河流域中游，隶属黄河一级支流，位于东经109°59′10″～110°30′00″，北纬37°59′18″～38°28′25″。河流全长93 km，流域出口断面建有申家湾水文站。佳芦河流域属于典型黄土梁峁丘陵沟壑区，地貌特征为谷寨沟深、梁峁起伏。流域内植被稀疏，整体覆盖率较低[19]。多年平均降水量为410.3 mm，年内分配不均，多集中于6-9月，多年平均潜在蒸散发量为949.2 mm。

2.2 数据来源

本文采用的水文气象数据，主要包括1960-2012年流域降水、径流、蒸发及气温、风速、湿度等要素的时间序列。其中降水、径流资料摘自《黄河流域水文年鉴》，流域邻近榆林气象站的气象资料均来源于国家气象科学数据共享服务平台,并基于Penman-Monteith公式[20]计算获得潜在蒸散发量序列。佳芦河流域地形、水系及测站分布情况如图1所示。

图1 佳芦河流域地形、水系及测站分布图

2.3 研究方法

2.3.1 水文要素演化趋势及突变检验 本文采用Mann-Kendall方法对比分析降水量、径流量、潜在蒸散发等水文气象要素时间序列的趋势变化[21]；运用双累积曲线斜率法[22]和有序聚类法[23]探究径流序列的突变情况。

2.3.2 Budyko弹性系数法 Budyko[24]假设认为流域多年平均蒸散发同时受大气对陆面的水分供给(降水)和能量供应条件(净辐射和潜在蒸散发)两个因素的控制。其边界条件是：在干旱条件下，潜在蒸散发量大于降水量，降水量全部转化为蒸发量；在湿润条件下，潜在蒸散发量小于降水量，可用于蒸散发的能量都将转化为潜热。在前人研究的基础上，Budyko提出了满足边界条件的水热耦合平衡方程。基于Budyko假设，许多学者提出了一系列可以有效评价气候、植被和水文循环相互影响的Budyko修正模型，其中我国气候学家傅抱璞[25]在经过量纲分析和数学推导后得到的Budyko解析表达式被广泛应用。

(1)流域水量平衡方程：

R=P-E-ΔS

(1)

式中：R为径流深，mm；P为降水量，mm；E为实际蒸散发，mm；ΔS为流域蓄水变化量，mm，对于多年时间尺度，通常可以取为零值。

以傅抱璞公式为基础，Yang等[26]把Choundhury公式[27]结合成另一种表达Budyko曲线的方程，即Budyko假设中的Choudhury-Yang流域水热耦合平衡方程，采用该方程计算得到下垫面参数，将其代入公式(1)得到：

(2)

式中：E0为多年平均潜在蒸散发，mm；ω为反映流域下垫面特征的参数。

(2)对径流变化贡献的定量计算:

基于Budyko弹性系数法计算得到气候变量和下垫面条件变化对径流变化的贡献。径流的气候弹性系数定义为单位因素变化导致流域内径流量相应的变化程度，径流的下垫面弹性系数可以定义为由于流域下垫面变化导致的径流变化量，对该研究区而言，流域下垫面变化主要是指人类活动所引起的下垫面变化，严格意义上下垫面变化也受气候因素影响，但目前没有可以定量描述下垫面变化的公式，因此流域径流量影响程度中下垫面变化的研究尚存在不足。

视公式(1)、(2)中R、P、E0为相互独立的变量，由公式(2)可得：

R=f(P，E0，ω)

(3)

径流对特定独立气候变量Xi的弹性系数εXi可表示为：

(4)

式中：Xi为气候变量P、E0或下垫面参数ω。

降水量、潜在蒸散发量、下垫面特征参数的弹性系数计算式分别如下[28]：

(5)

(6)

(7)

各因子对径流变化的影响可以用下面的微分方程来计算：

(8)

式中：dP、dE0、dω分别为突变前后两阶段的多年平均降水量、潜在蒸散发量和下垫面特征参数的差值。

将公式(5)～(7)代入公式(8)可以估算气候和下垫面条件变化引起的径流量变化：

(9)

某单个因素对径流量变化的贡献率可以计算如下：

(10)

式中：CXi为因子Xi对径流量变化的贡献率。

3 结果与分析

3.1 水文气象要素变化趋势分析

采用Mann-Kendall秩次相关检验法对1960-2012年佳芦河流域的降水量、径流量和潜在蒸散发量的变化趋势分析检验结果如表1和图2所示。

表1 1960-2012年佳芦河流域各水文气象要素 M-K趋势检验表(α = 0.05，Uα/2 =1.96)

由表1和图2可知，从5 a滑动平均线来看，年平均气温呈波动式上升，最低值(7.1 ℃)与最高值(10.0 ℃)相差2.9 ℃，说明佳芦河流域年平均气温年际变化较大；潜在蒸散发量在1980年前变化平稳，1980年后历经先减小后增大的变化过程，经检验，增大趋势不显著；最大降水量发生在2012年，最小降水量出现在1965年，降水量总体呈不明显增大趋势；径流深自20世纪70年代以来呈持续显著减小变化趋势，其5 a滑动平均值呈现阶梯状。综合佳芦河流域的水文气象要素分析，年平均气温、潜在蒸散发量和降水量均呈增大趋势，而径流深呈显著减小趋势，这意味着气候变化是该流域径流量变化的重要影响因素，但径流量呈锐减的变化态势，需要深入进行归因分析研究。

图2 1960-2012年佳芦河流域各水文气象要素变化趋势

3.2 径流突变特征分析

采用径流深、降水量双累积曲线法和有序聚类法进行径流深序列的突变分析，结果如图3所示。双累积曲线(图3(a))在1978年出现一个明显的转折点，而且在转折点后的曲线斜率明显变缓，说明径流序列有减小趋势；由有序聚类曲线(图3(b))可看出，离差平方和在1978年达到最小值，出现了明显的突变。因此两种方法均表明径流序列在1978年发生突变。

图3 1960-2012年佳芦河流域径流深序列突变分析

因此，可将研究期径流序列划分为两个阶段，即1960-1977年为基准期，该阶段径流量未受人类活动干扰或干扰较小。1978-2012年为变化期，该阶段多年平均径流深为34.8 mm，较基准期减小49 mm，相对变化率为-58.5%；降水量较基准期减少28.1 mm，相对变化率为-6.6%，潜在蒸散发量较基准期增加5.9 mm，相对变化率为0.6%。

3.3 径流变化的归因分析

3.3.1 弹性系数 佳芦河流域不同时期的降水量、潜在蒸散发量以及下垫面特征参数的弹性系数分析计算结果见表2。由表2可知：(1)从整个研究时段(1960-2012年)来看，径流深对降水量和潜在蒸散发量的弹性系数分别为2.64和-1.64，表示降水量和潜在蒸散发量每增加1%，径流深会随之增加2.64%和减少1.64%；(2)对比基准期和变化期各水文要素的弹性系数可知，下垫面特征参数的弹性系数值从1.73增加至2.40，相对变化率为38.7%，说明在1978-2012年期间人类活动较为强烈，对径流变化产生了深刻影响；(3)降水量的弹性系数从2.05增加到2.65，说明降水量每增加1%，在1978年以前会使径流深增加2.05%，在1978年以后的径流深会随之增加2.65%，表明径流深对降水变化的响应程度发生了较大变化；(4)潜在蒸散发量的弹性系数从基准期的-1.05减少到-1.65，说明潜在蒸散发量每增加1%，在1978年以前会使径流深减小1.05%，在1978年以后的径流深会随之减小1.65%。总之，与基准期相比，佳芦河流域径流对气候变化和下垫面条件变化的敏感程度均有所增加。

表2 1960-2012年佳芦河流域不同时期各特征参数的弹性系数计算结果

气候与下垫面特征参数的弹性系数年际变化情况如图4所示。各因子弹性系数的绝对值呈显著增加趋势(α=0.05)，表明佳芦河流域的径流对气候变化和下垫面条件变化均较为敏感。

图4 1960-2012年佳芦河流域各特征参数弹性系数的年际变化

3.3.2 径流变化的量化归因分析 佳芦河流域气候变化和下垫面条件变化对径流深变化的量化归因分析结果见表3。

表3 佳芦河流域径流深变化的量化归因分析结果

由表3可以看出，与基准期相比较，降水量、潜在蒸散发量、下垫面条件等变化引起的径流深的变化量分别为-6.57、-0.35和-23.91 mm，即变化期气候因素和下垫面条件变化共同作用引起径流深的减少。气候变化对径流深变化的贡献率为22.46%，其中降水量为主要影响因素，贡献率为21.31%；下垫面条件变化对径流深变化的贡献率为77.54%。综上所述，人类活动和气候变化均为引起该流域径流量深显著减小的原因，其中人类活动的影响最大、降水量的影响次之，潜在蒸散发量的影响最小。

4 讨 论

大量研究表明，累积量斜率变化率[8]比较法适用于干旱半干旱区域，采用该方法定量分析佳芦河流域气候变化和人类活动对径流变化的贡献，并与前文基于Budyko假设的弹性系数法计算结果进行比较，以探讨弹性系数法应用结果的可靠性。

根据径流序列的突变点分析结果，可计算并获得基准期(1960-1977年)和变化期(1978-2012年)降水量、潜在蒸散发量和径流深的累积量斜率变化，如表4所示。由表4可知，径流深变化率为-62.35%，呈显著下降趋势，其中降水量对径流深减少的贡献率为9.67%，潜在蒸散发量对径流深减少的贡献率为1.55%，人类活动对径流深减少的贡献率达88.78%。对比弹性系数法，二者计算出的对径流深贡献率略有差异，但均表明，人类活动和气候变化均为引起该流域径流深显著减小的原因，其中人类活动的影响最大。

表4 佳芦河流域不同时段累积降水量、潜在蒸散发量、径流深的斜率变化以及对径流深减少的贡献率

前文基于弹性系数法的径流变化归因分析结果表明，基准期和变化期人类活动的下垫面特征参数的弹性系数分别为-1.73和-2.40。经调研分析，为保护和改善黄土高原生态环境，陕西省已采取多种有效措施，如开展退耕还林试点，其中佳芦河林草地面积小幅增大，林草地的持水性能和截留量均有所增加；从佳芦河四期土地利用数据得出1986-2008年耕地和未利用土地面积持续减小，耕地面积比例从53.12%减小到50.54%，未利用土地面积比例从6.06%减小到3.05%；林地和草地持续增加，林地面积比例从1.56%增加到2.17%；草地面积比例从38.95%增加到43.94%[29]。此外，佳芦河流域在20世纪70年代修建淤地坝速度迅猛(1956-2008年淤地坝建设数量如图5所示)，淤地坝是在水土流失地区各级沟道中，以拦泥淤地为目的而修建的坝工建筑物，前期修建的淤地坝也在20世纪70年代共同发挥作用。以上各种人类活动不断改变着流域下垫面条件，综合作用影响了佳芦河径流的变化态势。

图5 1956-2008年佳芦河流域淤地坝建设数量

5 结 论

本文分析了佳芦河流域水文气象要素演变情势，并定量计算了径流深变化影响要素的贡献程度，得出以下主要结论：

(1)研究区1960-2012年降水量、潜在蒸散发量均无显著变化趋势，但径流深序列存在显著的下降趋势，且在1978年发生突变，变化期较基准期多年平均年径流深减少49 mm，相对变化率为-58.5%。

(2)径流深对降水量、潜在蒸散发量及下垫面变化的弹性系数分别由基准期的2.05、-1.05和-1.73变化为变化期的2.65、-1.65和-2.40，且各因子弹性系数的绝对值呈显著增加趋势，说明佳芦河流域的径流深对气候变化和下垫面条件变化均较为敏感。

(3)基于Budyko假设框架的弹性系数法分析表明，人类活动和气候变化均为引起佳芦河流域径流深显著减小的原因，其中人类活动、降水量和潜在蒸散发量对径流深变化的贡献率分别为77.54%、21.31%和1.15%，累积量斜率变化率法的分析结论与之相似。