基于Budyko假设的黄河源区径流演变与量化归因分析

马明卫, 王召航, 王文川, 罗党

(1.华北水利水电大学 水资源学院,河南 郑州 450046; 2.华北水利水电大学 水利学院,河南 郑州 450046;3.华北水利水电大学 数学与统计学院,河南 郑州 450046)

在全球气候变化和人类活动影响加剧的背景下,黄河源区径流量日趋减少,威胁着流域的生态安全,可能损害区域经济的高质量发展,这一问题引起了政府决策部门和社会公众的广泛关注[1]。影响径流形成和变化的因素主要有三大类,即气候因素(降水、气温和蒸散发等)、流域下垫面条件(地形、土壤和植被状况等)以及人类活动(淤地坝建设、过度农业灌溉和直接取用水等)。其中,植被作为陆地生态系统的重要组成部分,在过去几十年经历了巨大的变化,北极地区、欧亚北部地区、北美和发展中国家(特别是中国和印度),出现了植被绿化的现象[2]。植被变化不仅影响气候,而且可以通过改变截留蒸发、植被蒸腾、降水入渗、地下水补给等水文过程,进而影响径流的形成与演变[3]。因此,在径流变化的成因分析中,准确把握植被变化对径流的影响十分重要。

目前,关于径流演变与量化归因分析的研究大多采用以下两类方法:第一类是基于过程的水文模型法,虽然水文模型法分析精度高,物理机制强,可以对降水、蒸散发、径流等要素进行水文过程的模拟,但模型参数率定困难且受人的主观因素影响较大,增加了结果的不确定性[4];第二类方法是基于Budyko假设的弹性系数法,与基于过程的水文模型法相比,结构简单、使用方便,而且对参数的选择和模型构建的要求相对较低,它在降低不确定性方面具有一定优势[5]。以往人们通常采用基于Budyko假设的弹性系数法来定量区分气候变化和人类活动对径流的影响,其中人类活动对径流变化的贡献是通过径流变化总量与气候对径流变化影响分量的差值来估计的,即流域下垫面变化(包括植被状况变化)也包含在人类活动影响中,这样就不能精细地定量分析植被变化是如何对径流产生影响的[6]。在Budyko方程中,下垫面参数n反映的是流域地形、土壤、土地利用和植被覆盖等综合效应[7],有学者通过计算下垫面参数对径流变化的贡献率,来粗略地分析植被变化对径流演变的影响[8-10]。但是,关于植被变化对径流影响的相关研究还比较少且不够深入,不能精细地计算出植被变化对径流变化的贡献量。

归一化植被指数(Normalized Difference Vegetation Index,NDVI)是国际上公认的反映大尺度地表植被覆盖和生长状况的有效指标,通常用来表征植被活动的强弱[11]。本文选取近年来植被恢复较为明显的黄河源区,根据弹性系数法的定义,在Budyko框架中引入植被指数NDVI,将影响径流的因素进一步细化,扩展Budyko框架下的具体应用公式,量化分析气候因素、流域下垫面条件(包括植被状况)和人类活动对径流变化的贡献,以期为黄河源区水资源保护和利用提供参考。

1 研究区概况

黄河源区位于青藏高原东北部,地理位置位于32°10′N～36°05′N和95°50′E～103°30′E范围内,如图1所示。其中唐乃亥水文站以上流域面积约12.2万km2,占黄河流域总面积的16.2%;源区年径流量为227.5亿m3(2018年),占黄河流域年径流量的35%。源区内平均海拔在3 000 m以上,以山地地貌为主,山脉绵延,地势高耸,地形复杂。区域气候属于高原大陆性气候,表现为冷热两季交替,干湿两季分明,年温差小,日温差大,日照时间长,辐射强烈,无四季区分的气候特征,多年平均气温为-0.35 ℃,年均降水量为361.60 mm[12]。

图1 黄河源区概况图

2 数据与研究方法

2.1 数据来源与处理

黄河源区径流资料为唐乃亥水文站1956—2018年的逐日实测流量数据,气象数据来自中国气象局国家气象中心(http://www.data.cma.cn/)。综合考虑黄河源区所有气象站的气象数据序列长度、数据完整性和空间位置等特性,最终选用了河南、久治、达日、玛曲、若尔盖、玛多和清水河7个气象站1956—2018年的降水和气温观测数据。

归一化植被指数NDVI能够表征陆面植被覆盖和植被生长情况,是植被宏观检测的重要指标,故采用NDVI来表征植被变化。本研究使用的GIMMS NDVI数据来自美国国家航空航天局(http://ecocast.arc.nasa.gov/data/pub/gimms/),时间范围为1982—2015年,空间分辨率为8 km,时间分辨率为15 d。根据研究区的空间范围,提取相应网格单元逐月NDVI数据,再对每年12个月的月尺度NDVI数据取平均值得到黄河源区1982—2015年的年均NDVI序列。

潜在蒸散发数据采用Hamon公式计算求得。Hamon方法认为气温是驱动蒸散发的主要动力,同时也考虑了日照时长与饱和蒸气压对蒸散发的影响,该方法的优点在于只需要气温数据便可以估算潜在蒸散发量[13]。

2.2 研究方法

2.2.1 趋势及突变检验

针对长时间序列数据,采用线性趋势法和Mann-Kendall(M-K)非参数统计检验法综合分析黄河源区气温、降水、潜在蒸散发、径流和NDVI的变化趋势。利用Mann-Kendall突变检验法识别相应气象、水文、植被要素演变中的突变现象,结合滑动T突变检验法[14]来去除M-K突变检验法中的干扰点,从而进一步确定合理的突变点。根据年径流序列的突变点,划分基准期和变化期。

2.2.2 植被变化对径流影响的研究方法

流域水量和能量平衡是Budyko理论的基础,其中水量平衡表达式为:

R=P-Ea-ΔS。

(1)

式中:R为年径流深,mm;P为年降水量,mm;Ea为年实际蒸散发量,mm;ΔS为年蓄水变化量,mm。对于闭合流域,长时间尺度的蓄水变化量很小(近似为0),为了减少黄河源区流域ΔS的影响,取各变量时间尺度为10年的滑动平均值。

Budyko理论认为实际蒸散发量主要取决于降水量和潜在蒸散发量之间的平衡,其两个边界条件为:在极端干燥的环境下,所有降水量都将蒸发;在极端湿润的环境下,蒸散发消耗的能量都将转化为潜热[15]。

YANG H B等[16]在Budyko理论的基础上,推导出了实际蒸散发量的经验方程:

(2)

式中:E0为年潜在蒸散发量,mm;n是综合反映流域下垫面特征的无量纲参数(简称下垫面参数),与地形、土壤、植被覆盖和土地利用等一系列因素有关。

公式(2)中n反映的是流域下垫面变化情况,所以n值的大小在年或者多年尺度上应该是变化的。对于一个特定流域,土壤和地形等状况较少会发生变化,故下垫面参数n主要受流域植被状况变化的影响,根据LI D等[17]的研究结果,下垫面参数n与植被指数NDVI有如下线性关系:

n=a·NDVI+b。

(3)

其中,对于特定流域,a和b是常数,a的大小表征了n对NDVI变化的敏感性程度,b则反映除了NDVI以外其他因素对n的影响。

2.2.3 弹性系数法与贡献率计算

根据弹性系数法的定义[18]推导出径流量R对降水量P、潜在蒸散发量E0和下垫面参数n的弹性系数计算公式:

(4)

(5)

(6)

式中:εP、εE0、εn分别为径流量R对降水量P、对潜在蒸散发量E0、对流域下垫面参数n的弹性系数。

弹性系数为正值,表示径流量会随着该影响因子的增加而增加;为负值,则说明径流量会随着该影响因子的增加而减小。弹性系数的绝对值越大,该影响因子对径流量的影响程度越大。

根据全微分的概念,可以用微分方程表示各要素对径流量R变化的贡献量,计算公式如下:

(7)

将公式(7)两边同时除以R,结合公式(4)(5)(6)得到每个因子对径流量R相对变化量的表达式:

(8)

式中:等号右边3项分别是P、E0、n变化引起径流量R的相对变化量。

本文根据弹性系数的定义,并结合公式(3)将植被指数NDVI引入到基于Budyko假设的弹性系数法中,以Budyko方程中的下垫面参数n作为中间变量,来定量分析植被变化对径流变化的影响程度,其中径流量R对NDVI的弹性系数为:

(9)

(10)

各要素对径流变化的贡献率ηxi的计算公式如下:

(11)

式中:ΔRxi为要素xi变化引起的径流量变化;ΔR为所有要素变化引起径流量变化的总和。

3 结果分析

3.1 水文气象要素演变趋势与径流突变点检测

黄河源区1956—2018年年平均气温和年平均潜在蒸散发量的年际变化趋势如图2所示。由图2可知:气温变化倾向率为0.032 ℃/年,潜在蒸散发量变化倾向率为1.026 mm/年;通过气温和潜在蒸散发量的线性拟合结果可以看出,黄河源区1956—2018年气温和潜在蒸散发量上升趋势明显。M-K趋势检验结果见表1,该表显示,黄河源区1956—2018年气温和潜在蒸散发量呈显著上升趋势(显著水平p>0.05),与线性趋势法检验得出的结果相一致。

表1 黄河源区各要素M-K趋势检验结果

黄河源区1956—2018年降水量和径流深的年际变化趋势如图3所示。黄河源区在过去60年间降水相对丰沛,多年平均降水量达到361.65 mm,且有逐年上升的趋势,但是黄河源区1956—2018年的径流量却并没有随着降水量的增加而增加,其变化倾向率为-0.189 mm/年。这说明降水量的变化不是导致黄河源区径流量减少的因素,潜在蒸散发、其他环境因素和人类活动可能是造成径流量减少的原因。

图3 黄河源区1956—2018年降水量和径流深变化趋势

为了进一步研究黄河源区径流的变化特性,采用M-K和滑动T两种突变检验方法综合识别径流变化的突变点。M-K突变检验结果如图4所示,UF和UB两条曲线相交于1991年和1993年两点,均在0.05显著水平范围内。结合滑动T突变检验结果(图5,子序列长度大于10年),最终确定1991年为黄河源区径流变化的突变点。从图5中可以看出,径流量在1991年前后发生突变后开始持续下降。为了进一步厘清其中的机理,下文将研究时段划分为1956—1991年和1992—2018年两个时段进行详细讨论。

图4 黄河源区1956—2018年径流M-K突变检验结果

图5 黄河源区1956—2018年径流滑动T突变检验结果

3.2 下垫面参数n与植被变化的关系

为了定量分析黄河源区植被变化对下垫面参数n的影响,首先对黄河源区植被变化趋势进行分析。图6给出了黄河源区1982—2015年年均NDVI变化情况。由图6可知:年均NDVI呈波动增加趋势;在研究时段内,黄河源区年均NDVI的变化范围是0.337 2～0.383 5,最小值和最大值分别出现在1983年和2010年。总体上,黄河源区1982—2015年年均NDVI线性变化斜率为0.000 4/年(p>0.05),植被恢复较为明显,黄河源区整体上呈变绿趋势。

图6 1982—2015年黄河源区年均NDVI年际变化

通过公式(1)和公式(2)联立,并由已知的径流量R、降水量P和潜在蒸散发量E0,即可推求出黄河源区在不同年份的下垫面参数n。然后根据公式(3),利用植被指数NDVI和下垫面参数n拟合出方程中的系数a和b,结果如图7所示。经拟合,得到a=25.292,b=-6.989 4,即下垫面参数n与植被指数NDVI的关系为n=25.292NDVI-6.989 4(p>0.05)。p=0.017 4>0.05,达到95%的置信度,判定系数R2=0.836 6,说明模型有效且拟合的效果较好。由参数n与NDVI的关系式可以看出,植被变化是流域下垫面变化的主要控制因子。根据系数a(a=25.292)可知,NDVI每增加1%,下垫面参数n就增加25%以上,这意味着黄河源区流域产流性能对植被变化有着极强的敏感性。

图7 黄河源区年均NDVI与下垫面参数n的关系

3.3 径流演变量化归因

根据上文研究结果,将研究时段划分为基准期(1956—1991年)和变化期(1992—2018年)。表2详细列出了黄河源区基准期和变化期的多年平均径流深R、降水量P、潜在蒸散发量E0和下垫面参数n的具体数值以及相应的变化量。其中,变化期的多年平均径流深和降水量较基准期的减少,多年平均潜在蒸散发量和下垫面参数较基准期的增加。

表2 径流突变点前后各要素的变化量

为了进一步对径流变化进行归因分析,根据公式(4)(5)(6)和公式(9)计算径流量对降水量、潜在蒸散发量、下垫面参数和植被指数的弹性系数,结果见表3。由表3可以看出,径流量变化与降水量变化呈正相关,与潜在蒸散发量、下垫面参数和植被指数变化呈负相关。具体表现为,当降水量、潜在蒸散发量、下垫面参数和植被指数NDVI均增加10%时,在变化期,径流量将分别增加15.9%、-6.2%、-13.0%和-9.8%。

表3 径流突变点前后各要素的弹性系数

根据公式(11)求得气候因素、下垫面条件和人类活动对径流变化的贡献率见表4。结果表明,降水量对径流量的减少为负贡献,潜在蒸散发量、植被指数和人类活动对径流量的减少为正贡献。其中,1992—2018年降水量对径流量减少的贡献率为-35.45%,潜在蒸散发量对径流量减少的贡献率为9.24%,下垫面条件对径流量减少的贡献率为55.61%(植被指数NDVI占25.30%),人类活动对径流量减少的贡献率为70.60%。总体上看,人类活动是黄河源区径流量减少的首要因素。气候因素对径流量减少的贡献率变化不大,对径流量减少为负贡献,抑制了径流量减少的趋势。流域下垫面变化是造成径流量减少的另一个重要因素。由表4可知,植被恢复导致径流量的减少抵消了大部分降水量增加引起的径流量增加,并且由于潜在蒸散发量的增加,导致过去几十年径流量的总体小幅下降。这一结果强调了植被变化对控制黄河源区长期径流演变的重要性,需要引起一定的重视。

表4 各影响因子对径流量变化的贡献率 %

4 讨论

近几十年来,黄河源区径流量总体上呈下降趋势[19]。张成凤等[20]采用基于Budyko假设的弹性系数法定量识别黄河源区径流变化的原因,其研究结果表明,气候因素使径流量先减少后增加,其他环境驱动因素是黄河源区径流量减少的主要因素,但其研究只能将冰川冻土和人类活动归类为其他环境驱动因素,不能更精细地量化分离各要素是如何对径流产生影响的。商滢等[21]认为,降水不是引起黄河源区径流变化的原因,降水对径流量减少的贡献率为6.1%,人类活动和其他气候因子对黄河源区径流产生了巨大影响。叶培龙等[22]的研究显示,从年际波动来看,降水是黄河源区年径流量变化的主要影响因子,而生态植被、冻土退化、水储量变化和人类活动对径流量的影响也不容忽视。由此可见,以往的研究在径流演变驱动机制方面存在一定的分歧。

自1998年我国实施退耕还林还草工程以来,黄河源区整体上呈变绿趋势,目前还较少有学者研究植被变化是如何对径流演变产生影响的。本文根据弹性系数法的定义,在Budyko框架中引入植被指数NDVI,将影响径流的因素进一步细化。其中,气候因素、下垫面条件和人类活动是影响径流变化的3个主要因素,将气候因素划分为降水和潜在蒸散发两个因素,将下垫面条件划分为植被指数NDVI和其他会对下垫面变化造成影响的因素。通过量化分离造成径流变化的影响因素,发现植被指数NDVI对黄河源区径流的演变存在不可忽视的作用。具体表现为,植被恢复导致径流的减少量可以抵消大部分降水增加引起径流的增加量。

需要强调的是,本研究仍存在一些不足之处。在本研究中,流域下垫面条件对黄河源区径流减少的贡献率为55.61%,其中植被指数NDVI占25.30%,其他占30.31%。这里的“其他”可能包括流域内的土壤地形变化、湖泊水库变化和土地利用变化等。一般情况下对于特定的流域,其土壤、地形和湖泊水库较难发生改变,但近几十年来,黄河源区的土地利用却发生了较为明显的变化(图8)。

图8 1980年、2015年黄河源区土地利用情况

从图8中可以看出:黄河源区主要土地利用类型是耕地、林地、草地、水域、建设用地和未利用土地,其中又以草地、林地和未利用土地为主;2015年相较于1980年,耕地、林地和建设用地面积分别增加了46、461、23 km2,相对增幅分别为10.27%、5.59%、62.16%;草地和水域面积面积分别减少了669、162 km2,相对减幅分别为0.72%、5.85%。本研究中的植被指数NDVI对径流变化的影响可以反映出耕地、林地增加和草地减少对径流变化的影响,却还不能得出建设用地增加和水域减少是如何对黄河源区径流变化产生影响的。在未来的研究中,应进一步细化“其他”因素对径流变化产生的影响,使研究结果更加准确。

5 结论

本文利用黄河源区水文气象观测数据和卫星遥感NDVI数据,采用基于Budyko假设的弹性系数法,定量分析了黄河源区气候变化、流域下垫面条件变化和人类活动对径流演变的相对贡献,主要结论如下:

1)1956—2018年黄河源区降水量显著增加(p>0.05),但径流量呈现不显著下降趋势,下降速率为-0.189 mm/年,径流量在1991年前后发生突变。

2)在多年时间尺度上优化了Budyko假设下的具体应用公式,拟合得到黄河源区下垫面参数n与植被指数NDVI的函数关系为n=25.292NDVI-6.989 4(p>0.05)。根据系数a为25.292可知NDVI每增加1%,下垫面参数n则增加25%以上,这意味着黄河源区产流性能对植被变化有着极强的敏感性。

3)黄河源区的径流变化主要是由人类活动引起的,在径流量的总变化中占70.60%;气候变化对径流量减少的贡献率为-26.21%,抑制了径流量减少的趋势;流域下垫面条件变化是造成径流量减少的另一个重要因素,植被恢复可能会使黄河源区产流量和径流量进一步减少。