基于Budyko假设的永定河流域上游不同支流径流变化归因分析

李 秀， 郎 琪， 雷 坤， 程全国， 刘中伟， 商放泽， 孟翠婷

(1.中国环境科学研究院， 北京 100012； 2.沈阳大学, 辽宁 沈阳 110044； 3.水利部河湖保护中心，北京 100012； 4.中电建生态环境集团有限公司， 广东 深圳 518102)

1 研究背景

河川径流作为水循环过程的核心环节和水资源的重要组成部分，从多方面推动着自然生态系统的发育与演化[1]。近年来，气候变化引发的极端天气事件频发、降水量分布不均，人口增长和用水量增多造成我国区域水资源供需矛盾、水质恶化等问题日益突出[2]。印度[3]、中国[4]等多个国家河流径流情势已发生显著变化，其中自20世纪50年代以来我国海河流域的实测径流量减幅超过80%，严重威胁了区域水资源开发利用和优化配置。因此，探究径流演变规律及其驱动因子对制定应对气候变化的策略和水资源规划与管理具有重要支撑作用。

气候变化和人类活动被认为是导致水文循环系统发生变化的两大驱动因素。近年来，学者从定性转变为定量角度开展河川径流变化的归因分析，常用方法有水文模型法、水文统计法和基于Budyko假设理论的水热耦合模型法等[5]。Ricardo等[6]通过SWAT(soil and water assessment tool)模型研究发现，西班牙西北部径流受到气温和降雨变化的影响，大气中二氧化碳浓度的增加会适当减弱气候变化对径流的干扰；陈鑫等[7]基于SWAT模型还原了海河流域天然径流量，并发现人类活动对平水年和丰水年径流的影响更为显著。虽然水文模型可以评估下垫面变化对径流的影响，但模型运行需要大量的数据支撑，若缺失数据通过模型模拟进行插补，则存在较大的误差。基于Budyko假设理论框架的水热平衡关系具有物理意义明确、参数易获取等优点，可准确量化气候因子和人为因子对径流变化的影响[8]，国内外学者已应用该方法在加拿大、中国等多个国家开展了相关研究。如Xu等[9]选取海河流域33个子流域，基于Budyko假设定量评估了气候和下垫面对径流变化的影响，发现土地利用和植被覆盖对径流减少的贡献率为73.1%。Tan等[10]基于Budyko水热耦合方程评估了人类活动和气候变化对加拿大96个子流域年均径流的贡献，结果表明土地利用和覆盖变化(land use/cover change, LUCC)是近几十年径流量减少的普遍原因。在归因分析中，水文序列趋势诊断和变异识别是Budyko框架的关键步骤，主要采用Mann-Kendall突变检验、Pettitt检验和双累积曲线等方法[11]，通过统计特征值的变化来确定突变时间，但不同检验方法所依据的理论基础不同，对序列长度及数值分布要求也相异，目前并未有统一的划分原则，致使同一序列应用不同检测方法得到的结果略有差异[12]。

永定河上游是整个流域的主要水源涵养区与产流区，是中下游北京、天津等重点城市生活、工业用水的重要来源，研究永定河流域上游的河川径流演变对流域水生态环境良性发展与水资源开发利用具有重要意义。20世纪80年代以后，流域面临水资源短缺、下游河道生态来水量不足、水生态功能丧失、水体自净能力薄弱以及无法净化随着经济生活不断增加的污染物等一系列水生态环境问题[13]，因此，亟需深入揭示永定河流域气象水文要素的演变趋势并定量分析径流量锐减的主要原因。部分学者如莫崇勋等[14]、张利平等[15]针对永定河全流域径流演变和影响因素开展了相关研究，而关于流域上游不同支流径流变化的对比研究相对较少。因此，本文选取永定河上游桑干河和洋河子流域为研究区域，采用Mann-Kendall趋势法和双累积曲线法辨析实测径流量的演变规律及突变特征，并基于Budyko假设的互补关系模型量化气候因子和人为因子对径流量变化的贡献率，旨在为永定河流域水资源规划与管理提供理论依据。

2 研究区概况

永定河上游处于38°55′～41°26′N, 111°55′～116°20′E之间，包括山西省、内蒙古自治区、河北省部分地区，是重要水源涵养区和生态廊道。桑干河和洋河是永定河上游的两大支流，其中桑干河河长390 km，流域面积为2.58×104km2，约占研究区域总面积的55%，该地区煤矿资源丰富，是重要的能源示范区；洋河河长101 km，流域面积为1.67×104km2，占研究区域总面积的36%，该地区人口密集，城市化率较高，是未来工业科技发展的聚集地。研究区地理位置、流域水系及水文和气象站点分布见图1。

图1 研究区地理位置、流域水系及水文和气象站点分布

3 数据来源与研究方法

3.1 数据来源

1957-2017年径流数据来自海河流域水文年鉴。响水堡和石匣里水文站分别位于洋河和桑干河下游(图1)，其记录数据可体现两个子流域多年来径流量的变化，因此选取响水堡和石匣里水文站径流数据并根据各流域集水面积计算径流深。气象数据类型包括降水量、气温和日照时数等，来源于中国气象局国家气象信息中心(http://www.data.cma.cn)。利用泰森多边形法根据桑干河、洋河子流域气象站所占面积系数计算子流域降水量，桑干河流域气象站包括右玉、朔州、大同、代县、灵丘、蔚县；洋河流域气象站包括天镇、张北、张家口、怀来。潜在蒸散量根据Penman-Monteith公式计算求得[16]。

3.2 基于Budyko假设的互补关系法

前苏联气候学家Budyko[17]认为陆面多年年均蒸发量的主控因子为陆面水分供给和大气蒸发需求。其中水分供给可用降水量(P)表示，大气蒸发需求可用潜在蒸散量(E0)表示，则有：

(1)

式中：E和E0分别为实际年均蒸散量和潜在年均蒸散量，mm；P为年均降水量，mm；E0/P为干旱指数，用符号φ表示。

多年来，以我国学者傅抱璞[18]为代表的众多学者基于Budyko假设推导出表征下垫面特征的经验公式。本文选用2008年Yang等[19]引入可供水量概念，通过量纲分析和数学推导提出可应用于任意时间尺度的解析表达式：

(2)

式中：n为下垫面特征参数(体现人类活动、相对入渗能力、植被盖度、土壤性质和平均坡度等)，取值范围为(0，+∞)。

将永定河流域官厅水库以上看作一个闭合流域，假设在长时间尺度上(超过10 a)土壤含水量的变化可忽略不计，可采用降水量与实际蒸散量的差值来估算径流量[20]。Zhou等[21]假设降水量(P)和潜在蒸散量(E0)相对独立，提出径流量(R)对降水量和潜在蒸散量的弹性系数具有互补关系：

(3)

2016年，Zhou等[22]又基于互补关系提出了引入权重因子α划分径流量变化的互补法，该方法与全微分法相比，具有适用性广泛、计算结果更为准确并且可定义气候变化和人类活动对径流影响程度的范围等优点，其表达式如下：

(4)

式中：ΔR为径流量的变化量，mm；α为权重因子(0≤α≤1)，表示气候变化和人类活动不同路径对径流量变化的贡献量，本文使用作者推荐的α=0.5。角标1、2代表突变点前后的基准期和突变期。降水量、蒸散发和下垫面的变化量，即ΔRP、ΔRE0和ΔRc的计算公式如下：

(5)

(6)

(7)

(8)

(9)

4 结果与分析

4.1 径流变化特征分析

图2为1957-2017年桑干河和洋河子流域年均径流深变化趋势。分析图2(a)可知，1957-2017年桑干河子流域径流深均值为13.63 mm，减小速率为6.236 mm/10a，M-K趋势分析Z值为-8.41，|Z|>1.96，通过了置信度0.01的显著性检验，表明桑干河流域近60年径流深呈显著减小趋势。其中，1957-1969年径流深均值为35.24 mm，减小速率为13.790 mm/10a；1970-1983年径流深均值为14.91 mm，减小速率为0.250 mm/10a；1984-1998年径流深均值为7.20 mm，增大速率为1.540 mm/10a；1999-2017年径流深均值为2.96 mm，减小速率为0.100 mm/10a。分析图2(b)可知，1957-2017年洋河子流域径流深均值为18.23 mm，减小速率为6.520 mm/10a，M-K趋势分析Z值为-7.42，|Z|>1.96，通过了置信度0.01的显著性检验，表明洋河流域近60年径流深呈显著减小趋势。其中，1957-1983年径流深均值为30.66 mm，减小速率为5.270 mm/10a；1984-1998年径流深均值为14.02 mm，增大速率为4.230 mm/10a；1999-2017年径流深均值为3.88 mm，增大速率为1.700 mm/10a。

图2 1957-2017年桑干河和洋河子流域年均径流深变化趋势

本文基于双累积曲线法判别桑干河和洋河子流域降水量-径流深双累积曲线斜率突变点[24]，图3为1957-2017年桑干河和洋河子流域降水量-径流深双累积曲线。由图3(a)可知，桑干河子流域曲线斜率分别在1969、1983和1998年发生突变，因此将1957-1969年划分为天然基准期，该时期水文气候要素变化相对稳定，径流受人类活动影响程度较轻。1970-2017年为人类活动干扰期，该时期人口逐渐稠密，经济迅速发展，兴建水利、植树造林等人类活动对径流的干扰较为剧烈。为了进一步分析不同时期气候变化和人为活动对径流变化的影响，根据突变年份将桑干河干扰期划分为1970-1983年、1984-1998年、1999-2017年3个阶段。由图3(b)可知，洋河子流域曲线斜率在1983和1998年发生突变，因此将1957-1983年划分为天然基准期，1984-2017年为人类活动干扰期，该时期工业化进程加快，人为取用水量增加，人类活动对径流影响较大，进一步将干扰期分为1984-1998年和1999-2017年两个阶段。

图3 1957-2017年桑干河和洋河子流域降水量-径流深双累积曲线

4.2 气候变化特征分析

图4和5分别为1957-2017年桑干河和洋河子流域年降水量和年平均气温变化趋势。

图4 1957-2017年桑干河和洋河子流域年降水量变化趋势

由图4(a)可知，1957-2017年桑干河子流域年降水量均值为398.4 mm，最大值为592.7 mm(1964年)，最小值为217.5 mm(1965年)，两者相差375.2 mm；桑干河子流域年降水量总体呈微弱减少趋势，变化速率为-1.348 mm/10a，年际间波动较为剧烈。由图4(b)可知，1957-2017年洋河子流域年降水量均值为308.5 mm，最大值为446.1 mm(1959年)，最小值为211.6 mm(1965年)，两者相差234.5 mm；洋河子流域年降水量总体呈微弱减少趋势，变化速率为-2.390 mm/10a，年际间波动较为剧烈。由图5可知，1957-2017年桑干河和洋河子流域年均气温分别为7.08和7.48 ℃，气温均呈显著升高趋势，升高速率分别达0.321、0.296 ℃/10a。

图6为桑干河和洋河子流域1982-2017年年均NDVI值(归一化植被指数)变化趋势。图6表明，该时段内桑干河和洋河子流域年均NDVI值分别在0.222～0.324和0.227～0.327之间，多年均值分别为0.27和0.28；桑干河和洋河子流域年均NDVI值均呈增加趋势，增长速率分别为0.019/10a、0.016/10a。

图6 1982-2017年桑干河和洋河子流域年均NDVI值变化趋势

表1为研究时段内桑干河和洋河子流域降水量、平均气温及NDVI的Mann-Kendall趋势分析Z值。由表1可看出，桑干河和洋河子流域年降水量减少趋势并不明显，未通过0.05显著性检验；年均气温和年均NDVI的Z值绝对值均大于1.96，通过0.01显著性检验，表现为显著上升趋势。

表1 桑干河和洋河子流域年降水量、年平均气温及NDVI值的Mann-Kendall趋势检验Z值

4.3 气候和人类活动对径流变化贡献的量化分离

将基于水量平衡理论计算得出的实际蒸发量代入MCY(Mezentsev-Choudhury-Yang)公式(公式(2))反推出下垫面特征参数n，并将n值代入公式(8)、(9)得出敏感性系数。1957-2017年不同阶段桑干河和洋河子流域计算结果见表2。

表2 1957-2017年不同阶段桑干河和洋河子流域下垫面参数和径流敏感性系数计算结果

由表2可知，桑干河子流域基准期径流对降水量和潜在蒸散量的敏感性系数分别为0.24和-0.07，表明该时期降水量和潜在蒸散量每增加1 mm，将导致径流深分别增加0.24 mm和减少0.07 mm；干扰期径流对气候要素的敏感性持续减弱，整个干扰期降水量每增加1 mm仅引起径流深增大0.08 mm，潜在蒸散量每增加1 mm，仅使径流深减少0.03 mm。洋河子流域基准期径流对降水量和潜在蒸散量的敏感性系数分别为0.23和-0.05，表明该时期降水量和潜在蒸散量每增加1 mm，将导致径流深分别增加0.23 mm和减少0.05 mm；整个干扰期降水量和潜在蒸散量每增加1 mm，将使径流深增加0.09 mm和减少0.02 mm。总体来说，在研究时段内永定河上游气候要素对径流的影响程度持续减弱，其中降水量是径流变化的主导气候因子。

基于公式(5)～(7)计算得到不同时期桑干河、洋河子流域气候变化和下垫面变化对径流的影响量，并根据各因子影响量与径流总变化量(模拟ΔR)的比值计算出各因子对径流量变化的贡献率，计算结果如表3所示。

表3 桑干河和洋河子流域气候变化和下垫面变化对径流贡献的量化分离计算结果

表3表明：在桑干河子流域，1970-1983年(干扰期Ⅰ)与基准期相比径流深减少了20.3 mm，降水量、潜在蒸散量和下垫面变化均使径流减少，三者贡献率分别为32.5%、2.0%和65.5%；1984-1998年(干扰期Ⅱ)下垫面变化对径流的影响程度大幅增加，贡献率达到88.9%，降水量对径流的影响程度减弱，贡献率为12.5%，潜在蒸散量变化导致径流增加，贡献率为-1.4%；1999-2017年(干扰期Ⅲ)下垫面仍为影响径流的主要因素，但与干扰期Ⅱ相比略有下降，为85.4%，潜在蒸散量对径流的影响有所增大，这可能与该区域气候呈暖干化变化有关(参见图4、5)；从整个干扰期来看，与基准期相比，1970-2017年径流深共减少了27.5 mm，其中下垫面变化对径流减少的贡献率最大，为82.2%，其次降水量的贡献率为16.7%，潜在蒸散量的贡献率仅为1.1%。总体来说，下垫面变化对桑干河流域径流的影响自20世纪80年代开始大幅增强，是导致该地区径流减少的主要因素。在洋河子流域，1984-1998年(干扰期Ⅰ)与基准期相比径流深减少了16.6 mm，下垫面变化使径流深减少17.4 mm，贡献率为104.8%，而气候变化的贡献率仅为-4.8%，其中降水量使径流深减少1.1 mm，贡献率为6.6%，潜在蒸散量导致径流深增加1.9 mm，贡献率为-11.4%；1999-2017年(干扰期Ⅱ) 与基准期相比径流深减少了26.8 mm，下垫面变化对径流深的影响量为25.2 mm，贡献率为94.0%，气候变化的贡献率仅为6.0%，其中潜在蒸散量仍导致径流深增加0.4 mm，但贡献率与干扰期Ⅰ相比明显减弱，为-1.5%；1984-2017年(总干扰期)降水量、潜在蒸散量和下垫面对径流变化的贡献率分别为7.2%、-4.5%和97.3%。总体来说，洋河子流域下垫面变化对径流的作用占主导因素。

通过桑干河和洋河子流域的对比可知，下垫面变化对两个子流域径流变化均具有巨大影响，但洋河子流域由于人口、经济社会发展和人为取用水量增加等因素，人类活动对径流的影响更为明显，贡献率均超过90%。总体来看，桑干河和洋河子流域径流减少主要归因于人类活动的强烈干扰。

5 讨 论

本文基于Budyko假设互补关系法对比分析气候变化和人类活动对桑干河和洋河子流域径流的影响量及贡献率。研究发现下垫面变化是永定河上游径流减少的根本原因，降水量、潜在蒸散量和下垫面变化对桑干河和洋河子流域的径流贡献率分别为16.7%、1.1%、82.2%和7.2%、-4.5%、97.3%。

由统计分析可知，桑干河流域1970-1983年、1984-1998年和1999-2017年的降水量与1957-1969年相比分别减少了9%、5%、7%，径流深分别减少了57%、80%、91%，说明永定河上游降水总量与径流量变化趋势存在非一致性。但降水量年内分配和降水强度也与产流有关，汛期降水占比减少及大雨、暴雨等强降水频次降低均会使产汇流发生变化[25]，因此降水量虽然对径流变化具有一定的贡献率，但降水量并不是该区域径流变化的主导因素。此外，潜在蒸散量对径流减少的影响程度最小，多为负贡献，这可能与该区域存在“蒸发悖论”现象有关[26]，在潜在蒸散量远大于降水量的条件下，潜在蒸散量的适度减少对径流的影响较小[27]。

关于下垫面特征参数n的解析，本文主要围绕土地利用、植被覆盖两方面展开。研究发现桑干河和洋河子流域参数n在研究期内分别增加了1.11、0.85(见表2)，这表明随着人类活动范围和强度的加大，下垫面发生了剧烈变化，对径流的影响程度也呈加剧趋势。通过对永定河流域年均NDVI值的分析表明，20世纪80年代以来永定河上游植被覆盖率呈显著增大趋势(见图6)。永定河上游是我国水土流失重点治理区，自新中国成立以来采取了一系列生态治理与修复手段，《京津风沙源治理工程》(2001-2004，2009)开启了大规模退耕还林(草)、封山造林工作，此外《21世纪初首都水资源可持续利用规划》(2001-2005)、《永定河流域综合治理与修复工程》(2016)等项目的实施也使森林植被覆盖状况明显好转[28]。

不同的土地利用类型具有不同的产汇流机制，对地表径流的汇流贡献有所差异[29]。如城镇用地多为硬化地面，强降雨发生时地表径流短期内迅速增加并蒸散流失，下渗水量减少；耕地面积变化影响水资源的利用结构，进而改变其产汇流条件；草地和水域面积变化影响土壤含水量，对径流具有调节能力[30]。有关数据表明，1980-2015年期间永定河流域城镇用地面积持续增加，增幅达102%，农村居民用地增加了11.8%，其他建设用地增幅为50.6%；高覆盖度草地面积增加了4 459.5 km2，增幅为85.0%；水田面积从82.18 km2减少至8.2 km2[31]。总体来说，下垫面的变化是土地利用与植被覆盖等因素综合作用的结果，而目前研究无法定量区分土地利用和植被覆盖对径流变化的贡献率，对参数n的解析仍需该领域学者由定性解析向定量离析作进一步探索。

6 结论与建议

本文采用水文气象基础数据，基于Budyko假设的互补法从气候因素和人为因素两个方面定量揭示了桑干河和洋河子流域径流量锐减的主要原因及贡献率，主要结论如下：

(1)1957-2017年桑干河、洋河子流域径流深呈显著减小趋势，减小速率分别为6.236和6.520 mm/10a，主要在1969、1983和1998年发生突变。桑干河子流域基准期径流深变化速率为-1.370 mm/a，干扰期径流深减少了27.5 mm，变化速率为-0.320 mm/a；洋河子流域基准期径流深变化速率为-0.520 mm/a，干扰期径流深减少了22.2 mm，变化速率为-0.380 mm/a。

(2)1957-2017年桑干河、洋河子流域气温和NDVI值均呈显著上升趋势，两子流域气温上升速率分别为0.321、0.296 ℃/10a，NDVI值增大速率分别为0.019/10a、0.016/10a，降水量变化速率为-1.348、-2.390 mm/10a，减小趋势并不显著，年际间波动较为剧烈。

(3)基于Budyko互补关系法量化了气候变化和下垫面变化对径流变化的贡献率。桑干河流域干扰期1970-2017年下垫面变化对径流深减小的贡献率达82.2%，降水量和潜在蒸散量的贡献率分别为16.7%、1.1%。其中在1984-1998年下垫面影响大幅增强，贡献率达88.9%；洋河子流域干扰期1984-2017年下垫面、降水量和潜在蒸散量对径流深减小的贡献率分别为97.3%、7.2%、-4.5%。人类活动干扰是永定河上游径流减少的根本原因，其中洋河子流域径流受人类活动的影响更为显著。

本文基于Budyko假设，即假设长时间尺度下水库蓄水量等对流域径流的影响作用忽略不计的前提下开展分析，因此并未考虑水库闸坝等水利工程的影响。建议在未来研究径流变化时探寻新的特征参数，将水库调蓄及河道外取用水等人为活动因素考虑其中，细化各项人为影响因素对径流变化的贡献率，为流域水文循环过程研究与水资源开发利用提供技术支撑。