基于Budyko假设的三江源径流变化特性及量化分离

商放泽，王可昳，黄跃飞，魏加华，3

（1.清华大学水沙科学与水利水电工程国家重点实验室，北京100084；2.中电建生态环境集团有限公司，广东深圳518102；3.青海大学水利电力学院，青海西宁810016）

近130年（1880—2012年）以来，地球表面温度上升了0.85℃，预计到21世纪末全球地表平均温度将增加1.1～6.4℃［1-2］。近50年来中国地表气温升高了1.1℃，增温速率为每10年上升0.22℃，明显高于全球或北半球同期的增温速率［3］。研究发现全球气候变化极大影响了水循环和水文过程［4-5］，气候变化还会引起降水的巨大变化，并导致区域干湿极端现象频发［6］。三江源区由于青藏高原独有的热力和动力作用，形成了东亚特殊的季风气候系统，对气候变化更为敏感。

从20世纪90年代开始，在全球气候变化、人类活动以及青藏高原地壳不断隆升等自然和人类活动的共同影响下，三江源区出现了源区气温升高、永久冻土层厚度减少、降雨和蒸发强度改变、地表植被退化、冰川退缩、高原湖泊和湿地水源补给减少而造成湖泊湿地萎缩、土地荒漠化、生物多样性降低等一系列生态环境问题，对长江、黄河和澜沧江中下游地区的生态环境也带来了诸多不利影响［7-9］。有研究指出，气候变化将影响依赖青藏高原供水的长江、黄河和澜沧江的径流［10］。1958—2007年期间三江源区温度和降水量明显增加，长江源区和澜沧江源区的年径流量呈上升趋势［11］，也有研究表明最近14年长江源区的直门达水文站年径流共减少了96亿m3［12］；黄河源区径流从20世纪90年代至2007年出现下降趋势，枯季、雨季降水量与流量分别呈负、正反馈机制［13-14］。

气候变化和下垫面变化是径流变化的2个主要驱动因素。人类活动和植被覆盖变化造成下垫面的变化，在地形、地质、土壤和植被等因素中，植被覆盖是影响径流最为重要的下垫面因素［15-16］，而人类活动中的开垦和放牧对植被的影响最为直接。归一化植被指数（normal difference vegetation index，NDVI）是表征植被覆盖的主要参数之一。NDVI是遥感影像的近红外波段反射率与红光波段反射率比值的归一化参数，对环境变化反应十分敏感，且能够较好反映植物生长变化和植被覆盖状况，是植被生长状态和植被覆盖度的最佳指示因子［17-18］。有研究利用GIMMS和MODIS 2个卫星源的遥感监测值，以NDVI为基础，发现三江源地区植被生产力趋于好转，气候变化是影响植被生产力的决定性因素，而人类活动在一定程度上加快了其变化速率［19］。

在流域径流影响因素的贡献离析中，如何准确评估不同因素对径流变化的贡献是研究的难点。定量分析气候变化对径流变化的贡献（ΔRclimate）以及流域下垫面变化对径流变化的贡献（ΔRcatchment）的方法主要有水文模型法和基于Budyko假设水热耦合平衡理论的水量平衡法2种［20-22］。有学者使用水文模型定量分离气候变化和人类活动对汾河径流的影响，发现人类活动对径流减少的贡献占了64.1%［23］。虽然水文模型的方法可以分离出ΔRclimate和ΔRcatchment，但是存在水文循环的部分物理过程缺失、部分模型结构缺陷、系统误差等问题，因此这类模型的准确性有待提高［24］。利用Budyko假设方法方面，有研究基于Budyko假设，选择受人类直接取用水影响较少的山区小流域分析了松花江、辽河、海河、黄河和汉江等中国主要流域径流量减小的原因［25］。但是Budyko假设忽略了下垫面的影响。由于三江源区受到地形及其高原加热场的作用，气候和下垫面条件独特且复杂，径流及相关因子的趋势分析、影响因素贡献评估及模型构建方法和手段都需要有针对性地修正和改进，构建的模型要能有效反映三江源区气候因子和下垫面之间的相关关系和定量关系，目前尚缺乏有效的方法以定量离析各驱动因子对径流变化的贡献。

三江源区由于特殊的水文、地理和气候条件，是我国淡水资源的重要来源区和西部地区的生态屏障。但是三江源区的自然生态系统非常敏感和脆弱，三江源区生态环境的变化状况直接关系到我国西部地区的可持续发展，乃至中、东部江河下游地区的水和生态安全。因此，本文围绕三江源区径流对气候变化和下垫面变化的响应特征及定量表达这一科学问题开展研究，根据三江源地区水文气象站点实测资料，分析三江源、黄河源、澜沧江源气候、下垫面和径流的变化特性，并利用基于Budyko互补关系的新的径流量化离析方法分析气候变化和人类活动对三江源区径流的影响及贡献。以期为三江源区水源涵养和生态环境保护措施的制定提供依据。

1 研究区概况

三江源区位于青海省南部、青藏高原腹地，是长江、黄河和澜沧江的发源地。汇水区流域面积共31.26×104km2，其中长江流域面积13.77×104km2（直门达水文站以上），占源区总面积的44.0%；黄河流域面积12.20×104km2（唐乃亥水文站以上），占源区总面积的39.0%；澜沧江流域面积5.29×104km2（昌都水文站以上），占源区总面积的16.9%。为了保护三江源区生态环境，经国务院批准，2011年实施《青海省三江源国家生态保护综合试验区总体方案》，试验区范围为89°41′～102°40′E、31°53′～37°10′N。

2 数据来源与研究方法

2.1 数据来源与预处理

降雨、气温、蒸发和径流数据来源于青海省水文水资源勘测局和中国气象数据网（http：//data.cma.cn/）。部分站点在1956—2012年间的降雨、气温、蒸发和径流数据缺失，缺失的数据根据青海省水文水资源勘测局提供的相邻站点的数据进行插补。三江源降雨、气温、蒸发监测的气象和水文站点见图1，共33个气象站点和12个水文站点。降雨量利用泰森多边形进行计算，因此包括了三江源区外的气象站点。径流量为根据实测径流量进行还原的天然径流量。三江源监测径流的水文站点分布见图2，从图中可知，长江源区由直门达水文站和班玛水文站以上流域2个部分组成，黄河源区包括唐乃亥水文站和同仁水文站以上流域，澜沧江源区包括香达水文站和下拉秀水文站以上流域。因此，长江源径流量为直门达水文站和班玛水文站的叠加，黄河源径流量为唐乃亥水文站和同仁水文站的叠加，澜沧江源径流量为香达水文站和下拉秀水文站的叠加。

归一化植被指数利用卫星数据计算。卫星数据来自于相关组织的数据服务平台所公布的公开数据。1982—2012年月尺度GIMMS（global inventory modelling and mapping studies）AVHRR（advanced very high resolution radiometer）NDVI数据集来源于NASA Earth Exchange （http：//nex.nasa.gov/nex/），其空间分辨率8km。为了避免大气溶胶以及太阳高度产生的信号噪声，引入最大值合成法（maximum value composite）对NDVI数据进行预处理。最大值合成法假设云层是不断变化的，因此栅格点的最大NDVI值被选为该栅格点的NDVI值［26］。

图1 三江源气象和水文站点Fig.1 Meteorological and hydrological stations in the Three-River Headwaters Region

图2 三江源径流计算使用的水文站点Fig.2 Hydrological stations used for runoff calculation in the Three-River Headwaters Region

2.2 研究方法

2.2.1 线性趋势法

气候要素的趋势变化一般采用一元线性回归模型描述，即

式中：y为气候、下垫面和径流要素序列；x为时间序列；m为线性趋势项；c表示截距。

2.2.2 非参数Mann-Kendall法

利用非参数Mann-Kendall进行趋势分析，分析基于秩的非参数检验方法，辨别三江源区气候、下垫面和径流变化是自然特性的波动还是存在确定性的变化趋势。检验统计量S计算公式为［27］

式中：xi、xj分别为数据序列中的第i和第j个数据点；n为数据序列长度，也即样本数。符号函数表达式为

Mann-Kendall标准检验统计量Z的计算公式如下：

其中ti是第i组的数据点数目。

标准检验统计量Z服从标准正态分布。为了检验数据序列是否存在趋势，如果标准检验量Z的绝对值大于标准正态累计分布表中的，那么数据序列在显著水平σ上存在趋势性变化。90%、95%和99%置信度水平对应的Z的绝对值判别值分别为1.64、1.96和2.58。

衡量趋势大小的指标为

式中：Mz为中位数；1＜j＜i＜n。β为正值表示上升趋势，β为负值表示下降趋势。

2.2.3 Spearman法

Spearman秩相关系数为非参数性质的统计参数，与样本分布无关。Spearman秩相关系数［28］为

Spearman秩相关系数的符号表示X和Y之间的关系。如果Y随着X增大而增大，那么秩相关系数则为正；如果Y随着X的增加而减少，那么秩相关系数则为负。X和Y越接近严格单调的关系，Spearman秩相关系数的数值则越大。

Spearman秩相关的变化趋势检验统计变量，即检验统计量为

其服从自由度为v=N—2的student分布。当时，认为无明显变化趋势。

2.2.4 基于Budyko假设的互补关系权重因子法

采用基于Budyko假设的互补关系权重因子法［29］分析气候变化对径流变化的贡献（ΔRclimate）以及流域下垫面变化（或人类活动）对径流变化的贡献（ΔRcatchment）。

常用的Budyko方程有傅抱璞公式［30］以及MCY方程［31］，如下：

式中：E代表蒸散发量；P代表降水量；E0为潜在蒸发量，表示可利用的能量；E/P表示蒸散发系数；E0/P表示干旱指数Φ；参数ω和n代表下垫面特征，例如人类活动、相对入渗能力、植被盖度、土壤性质和平均坡度等。

假设P和E0是相互独立的，R对P和E0的偏弹性系数存在以下关系［29］：

根据式（11）可得出R的表达式为

基于R与P和E0的互补关系，即式（12），可得出计算径流变化的互补关系，如式（13）：

假设土壤含水量的变化量可以忽略不计，通过计算降水和径流的差值来估算蒸散发量［32］为

时段末径流量为

图3 状态空间（P，E0，E）示意Fig.3 State space(P,E0,E)schematic diagram

总径流变化（ΔR）计算公式为

代入径流变化的互补关系式（13），即可得到

或者

其中0≤α≤1。得到ΔRclimate和ΔRcatchment的计算公式为

随着权重因子α的变化，ΔRclimate和ΔRcatchment也随之线性变化。

在过去的30多年中［29-32］，由于气候条件和下垫面的不断变化，α并非一个定值，因此分析了α为0、0.5、1.0这3种变化路径中气候变化和下垫面变化对径流变化的影响，并分别阐释不同变化路径下气候变化和下垫面变化对径流变化的影响，其中α=0或α=1.0代表径流变化完全由下垫面变化引起或完全由气候变化引起。

潜在蒸散发量根据在中国干旱半干旱地区有着良好的应用的哈格里夫斯公式［33］计算，其计算公式为

式中：Ra表示大气层顶辐射，mm·d-1；Tmax和Tmin表示最高温度和最低温度；T代表平均温度。

3 结果与分析

3.1 三江源气候和下垫面变化特性分析

三江源降雨量变化趋势见图4，降雨量年际间波动变化。长江源、黄河源和澜沧江源降雨量呈增加趋势，年均降雨量分别为365.26、429.28和537.50 mm。线性趋势分析表明，10年的降雨量倾向率分别为10.83、5.51和11.14 mm。

三江源气温变化趋势见图5，气温年际间波动变化。长江源、黄河源和澜沧江源气温呈增加趋势，多年平均气温分别为—1.59、1.36和2.45℃。线性趋势分析表明，10年的气温倾向率分别为0.31、0.29和0.33℃。

利用归一化植被指数NDVI表征植被覆盖度，下垫面的变化情况用生态植被覆盖度表示。三江源NDVI变化见图6，生态植被覆盖情况年际间波动变化。长江源、黄河源和澜沧江源NDVI呈增加趋势，长江源、黄河源和澜沧江源年均NDVI值分别为0.230 5、0.314 1和0.320 4。线性趋势分析表明，NDVI值10年增加幅度分别为0.002、0.004和0.002。年均降雨量、气温和NDVI值从大到小排序为：澜沧江源、黄河源、长江源。

图4 三江源年降雨量变化Fig.4 The change of annual precipitation in the Three-River Headwaters Region

图5 三江源年平均气温变化Fig.5 The change of annual average air temperature in the Three-River Headwaters Regions

图6 三江源年NDVI变化Fig.6 The change of annual NDVI in the Three-River Headwaters Regions

三江源降雨、气温和NDVI的显著性检验及变化趋势见表1。Mann-Kendall、Spearman检验结果表明长江源和澜沧江源降雨量呈显著增加趋势，而黄河源降雨量增加趋势不显著；三江源气温呈显著增加趋势；黄河源NDVI呈显著增加趋势，长江源和澜沧江源NDVI增加趋势不显著。

3.2 三江源径流变化特性分析

三江源年径流量变化见图7，径流深变化见图8。长江源、澜沧江源年径流量呈增加趋势，年均径流量分别为140.45亿m3和58.22亿m3，线性趋势分析表明，10年径流增加幅度分别为6.54亿m3和1.54亿m3。黄河源年径流量呈减少趋势，年均径流量为206.26亿m3，线性趋势分析表明，10年终径流减少幅度为2.05亿m3。长江源、黄河源和澜沧江源区年均径流深分别为98.89、165.27、264.23 mm，年径流深从大到小排序为：澜沧江源、黄河源、长江源，这与长江源、黄河源和澜沧江源年降雨量、气温、NDVI的大小顺序相同。

表1 三江源降雨、气温和NDVI的Mann-Kendall和Spearman秩相关检验Tab.1 Mann-Kendall and Spearman rank correlation tests of precipation,air temperature and NDVI in the Three-River Headwaters Regions

图7 三江源年径流量变化Fig.7 The change of annual runoff amount in the Three-River Headwaters Regions

图8 三江源年径流深变化Fig.8 The change of annual runoff depth in the Three-River Headwaters Regions

三江源径流量的显著性检验见表2。Mann-Kendall、Spearman检验结果表明三江源区径流变化趋势都不显著。

3.3 三江源气候和下垫面对径流变化贡献的量化分离

为了显示径流的动态变化过程，同时为了满足“土壤含水量的变化量可以忽略不计”的前提假设［32］，选择5年作为时间步长。以5年的均值作为每个子时间段的取值，并认为相邻2个子时间段内土壤含水量的变化量忽略不计，当相邻2个子时间段土壤含水量基本保持不变时，降水和径流的差值可以作为蒸发量。因此，将总时间序列（1956—2012年）分为12个子时间段，前11个子时间段以5年为时间步长，最后1个子时段以2年为时间步长。

表2 三江源径流量的Mann-Kendall和Spearman秩相关检验Tab.2 Mann-Kendall and Spearman rank correlation tests of runoff amount in the Three-River Headwaters Regions

长江源、黄河源和澜沧江源每一个子时间段的气候条件和下垫面特征变化分别见表3、表4和表5。其中E为年均蒸散发量，E0为潜在蒸散发量，R为径流量，n为式（10）的参数，R1，R2，… ，R10为各子时间段的年径流均值。定义ΔR为相邻2个子时间段之间径流均值的差值，例如ΔR=R2—R1。通过互补关系法进行计算，得到相邻2个子时间段的径流差值以及气候变化和下垫面变化对径流变化的贡献值，长江源、黄河源和澜沧江源气候和下垫面变化对径流变化的贡献分别见表6、表7和表8。结果表明，三江源ΔR的计算值与观测值几乎相同，表明互补关系法可以准确地计算出径流年均值的变化。虽然ΔRclimate和ΔRcatchment的总和等于ΔR，但是从表6至表8可以看出，长江源、黄河源和澜沧江源ΔRclimate和ΔRcatchment的值随着权重因子α的改变而改变。

现大部分研究假定气候和下垫面变化对径流变化具有相同的贡献特征，即气候和下垫面同时促进径流量或者同时减少径流量。本文未假定气候和下垫面变化对径流变化具有相同的贡献特征，实际上，在某些时段或年份，气候和下垫面变化对径流变化具有相反的贡献特征。如长江源ΔR3中，ΔRclimate，1=18.61和ΔRcatchment，1=—7.01，气候变化对径流量贡献为正，而下垫面变化对径流量贡献为负。黄河源ΔR3中，ΔRclimate，1=12.94 及ΔRcatchment，1=—18.72，气候变化对径流量的贡献为正，而下垫面变化对径流量的贡献为负。澜沧江源ΔR9中，ΔRclimate，1=—9.57和ΔRcatchment，1=42.31，气候变化对径流量的贡献为负，而下垫面变化对径流量的贡献为正。

表3 长江源气候和下垫面参数Tab.3 The climate and underlying surface parameters in the source region of the Yangtze River

表4 黄河源气候和下垫面参数Tab.4 The climate and underlying surface parameters in the source region of the Yellow River

当α分别等于1.0、0.5和0时，长江源ΔRclimate的变化范围分别为（—25.19，29.84）、（—24.89，32.15）和（—24.59，34.46）mm，ΔRcatchment的变化范围分别为（—14.25，23.99）、（—15.22，21.75）和（—16.18，19.51）mm；黄河源ΔRclimate的变化范围分别 为（ —31.18，48.86）、（ —29.77，52.00）和（—28.36，55.15）mm，ΔRcatchment的变化范围分别为（—21.27，32.80）、（—21.94，32.53）和（—23.35，32.27）mm；澜沧江源的ΔRclimate的变化范围分别为（—57.74，62.96）、（—55.30，67.03）和（—52.86，71.11）mm，ΔRcatchment的变化范围分别为（—28.76，42.31）、（—30.23，41.15）和（—31.70，39.98）mm。澜沧江源径流年际变化较大，气候变化和下垫面变化对径流变化影响的绝对值最大；长江源径流年际变化较小，气候变化和下垫面变化对径流变化影响的绝对值最小。

表5 澜沧江源气候和下垫面参数Tab.5 The climate and underlying surface parameters in the source region of the Lancang River

3.4 三江源气候和下垫面对径流的影响

由于计算得到的ΔR为相邻2个子时间段之间的差值，整体上没有可比性，因此将计算得到的ΔR依次进行累加。即以子时间段T1为时间基准，分析ΔRclimate和 ΔRcatchment的变化规律。ΔR’climate，i和ΔR’catchment，i表示子时间段Ti+1与子时间段T1的差值。长江源、黄河源和澜沧江源区ΔR’climate，i和ΔR’catchment，i随时间的变化趋势分别见图9、图10和图11。

长江源ΔR’climate和降水具有显著的正相关关系，α为 1.0、0.5、0时，Pearson相关系数在 0.971～0.982，显著水平小于0.01。黄河ΔR’climate和降水具有显著的正相关关系，α为1.0、0.5、0时，Pearson相关系数在0.954～0.980，显著水平小于0.01。澜沧江源ΔR’climate和降水具有显著的正相关关系，α为1.0、0.5、0时，Pearson相关系数在0.906～0.937，显著水平小于0.01。由此可知，降雨是导致径流变化的主要来源与控制因子，降雨变化越大，则径流变化越大。

长江源ΔR’catchment，i和参数n具有显著的负相关关系，α分别为1.0、0.5和0时，Pearson相关系数在—0.994至—0.937之间，显著水平小于0.01。黄河源ΔR’catchment，i和参数n具有显著的负相关关系，α同样时，Pearson相关系数在—0.997至—0.994之间，显著水平小于0.01。澜沧江源ΔR’catchment，i和参数n具有显著的负相关关系，α同样时，Pearson相关系数在—0.988至—0.997之间，显著水平小于0.01。三江源自然保护区的建立，人类活动受到限制，若忽略人类活动的影响，则参数n的变化可以归因于植被覆盖度的变化。一般来说，植被覆盖越大，降雨入渗量就越多，产生的径流量就越少，因此ΔR’catchment，i与参数n的变化呈相反的趋势。植被覆盖是下垫面影响径流变化的主要控制因子。1956—2012年总体来看，气候变化和下垫面变化增加了长江源径流量，当α分别等于1.0、0.5和0时，长江源气候变化对径流的贡献分别为55.21%、74.03%和92.85%。气候变化和下垫面变化增加了黄河源径流量，当α等于1.0、0.5和0时，黄河源区气候变化对径流的贡献分别为52.47%、60.49%和68.50%。当α等于1.0时，气候变化和下垫面变化增加了澜沧江源径流量，澜沧江源区气候变化对径流的贡献为89.91%；当α分别等于0.5和0时，气候变化增加了澜沧江源径流量，但下垫面变化减少了澜沧江源径流量，因此澜沧江源区气候变化对径流的影响分别为107.25%和124.58%。长江源、黄河源气候和下垫面变化对径流的影响为正贡献，澜沧江源α=1.0时气候和下垫面变化对径流的影响为正贡献，但是澜沧江源α等于0.5和0时气候对径流的影响为正贡献，而下垫面对径流的影响为负贡献。

表6 长江源ΔRclimate和ΔRcatchment量化分离结果Tab.6 The quantitative separation results of ΔRcliamteand ΔRcatchmentin the source region of the Yangtze River

表7 黄河源ΔRclimate和ΔRcatchment量化分离结果Tab.7 The quantitative separation results of ΔRcliamteand ΔRcatchmentin the source region of the Yellow River

表8 澜沧江源ΔRclimate和ΔRcatchment量化分离结果Tab.8 The quantitative separation results of ΔRcliamteand ΔRcatchmentin the source region of the Lancang River

图9 长江源ΔR’climate和P的变化特性及ΔR’catchment和n的变化特性Fig.9 The variation characteristics of ΔR’climateand P,and ΔR’catchmentand n in the source region of the Yangtze River

4 结论

（1）长江源和澜沧江源降雨量呈显著增加趋势，黄河源降雨量增加趋势不显著，长江源、黄河源和澜沧江源10年降雨量倾向率分别为10.83、5.51和11.14 mm。三江源气温呈显著增加趋势，长江源、黄河源和澜沧江源10年气温倾向率分别为0.31、0.29和0.33℃。黄河源NDVI呈显著增加趋势，长江源和澜沧江源NDVI增加趋势不显著，长江源、黄河源和澜沧江源10年NDVI值增加幅度分别为0.002、0.004和0.002。

图10 黄河源ΔR’climate和P的变化特性及ΔR’catchment和n的变化特性Fig.10 The variation characteristics of ΔR’climateand P,and ΔR’catchmentand n in the source region of the Yellow River

（2）长江源、澜沧江源径流量增加趋势不显著，黄河源径流量减少趋势不显著，长江源、黄河源和澜沧江源10年径流量变化为6.54亿m3、—2.05亿m3和1.54亿m3。降雨量、气温、NDVI和径流深由大到小排序为：澜沧江源、黄河源、长江源。

（3）基于Budyko假设的互补关系权重因子法不仅可以准确地计算出径流年均值的变化，并且量化分离时无需假定气候和下垫面变化对径流变化具有相同的贡献特征。总体来说，长江源、黄河源气候和下垫面变化对径流的影响为正贡献，澜沧江源α（权重因子）等于1.0时下垫面变化对径流的影响为正贡献，但是α等于0.5和0时下垫面对径流的影响为负贡献。长江源、黄河源和澜沧江源气候变化对径流的贡献分别为55.21%～92.85%、52.47%～68.50%、89.91%～124.58%。

图11 澜沧江源ΔR’climate和P的变化特性及ΔR’catchment和n的变化特性Fig.11 The variation characteristics of ΔRclimateand P,and ΔR’catchmentand n in the source region of the Lancang River

（4）三江源气候变化对径流的累积贡献值与降雨量具有显著正相关关系，降雨量是气候变化影响径流变化的主要控制因子；三江源下垫面变化对径流的累积贡献值与参数n具有显著负相关关系，植被覆盖是下垫面变化影响径流变化的主要控制因子。

致谢：感谢清华大学水沙科学与水利水电工程国家重点实验室的郑裕彤博士、刘喆硕士在部分数据处理中给予的帮助。