雅鲁藏布江流域径流演变规律与归因分析

徐宗学,班春广,张 瑞

(1. 北京师范大学水科学研究院,北京 100875；2. 北京师范大学城市水循环与海绵城市技术北京市重点实验室,北京 100875)

全球气候和下垫面变化是影响流域水循环的2个重要因素。目前,在下垫面变化和不同气候条件的共同影响下,全球生态水文过程正在逐渐发生变化,成为当前水文学研究的热点领域之一[1]。青藏高原是世界上海拔最高、面积最大的高原,面积为250万km2,有“世界屋脊”之称[2]。青藏高原是亚洲许多河流的发源地及重要水源地,其水循环状况直接关系着周边及下游地区的水资源配置与防洪安全[3]。在气候变化背景下,降水[4]、潜在蒸散发、冰川及下垫面均发生了不同程度的变化[5],最终导致流域径流发生变化。因此,研究气候变化影响下青藏高原高山区流域水循环演变机理,揭示下垫面和径流演变规律,对区域水资源合理利用和规划管理具有十分重要的现实意义。雅鲁藏布江发源于青藏高原东南部,是一条重要的国际性河流[6]。近50 a来,雅鲁藏布江流域下垫面发生了一系列变化,如冰川加速退缩、冻土退化[7],对水循环过程和机制产生了深刻影响。

当前,国内外学者定量识别径流归因的方法主要包括水文模型模拟法、基于Budyko假设的弹性系数法和双累计曲线法[8]。近年来,基于Budyko假设的弹性系数法被广泛应用于国内外各大流域，如MOPEX流域[9]、长江流域[8]和汉江上游流域[4]。目前,针对雅鲁藏布江流域径流变化和径流归因分析等方面已开展大量研究[10- 13],但在全面揭示径流演变规律、定量厘清气候变化和下垫面变化对径流变化的影响方面仍存在不足。

雅鲁藏布江作为中国重要的国际性河流,蕴含丰富的淡水资源,该流域也是中国重要的生态安全屏障和生物多样性重点保护区域,其生态系统结构相对简单,生态环境脆弱,对气候变化尤为敏感。近50 a来,雅鲁藏布江流域水文气象要素发生了一系列变化。有研究表明,雅鲁藏布江流域降水以6.32 mm/(10 a)的速率呈不显著增加,气温以0.32 ℃/(10 a)的速率呈显著增加[14- 15]。流域中游年、季节平均气温显著增高,年降水呈现先减少后增加的趋势[16- 17]。Li等[18]基于MODIS NDVI数据分析了雅鲁藏布江流域的植被覆盖情况,发现该区域在2000年后的显著变绿过程是由气候- 下垫面- 水文过程之间的耦合关系所驱动。本文在分析气象要素和下垫面演变特征的基础上,研究变化环境下雅鲁藏布江流域的径流变化规律及其归因分析,相关研究结论可为高寒山区水资源管理、生态环境保护等工作提供科技支撑。

1 研究区概况

雅鲁藏布江发源于杰玛央宗冰川,在巴昔卡出境,进入印度阿萨姆邦改称布拉马普特拉河(Brahmaputra River),流域范围为28°N—31°N,82°E—97°E(图1)。雅鲁藏布江从杰玛央宗冰川的末端至拉孜为上游段,从拉孜到奴下为中游段,奴下到巴昔卡附近为下游段[12],干流全长2 057 km,流域面积约为24万 km2。雅鲁藏布江水量非常丰沛,年径流总量为1 395.1亿m3,年平均流量为4 425 m3/s[2]。流域平均海拔为4 000 m,地势西高东低。由于复杂的地理条件,各地气温条件差异较大,年降水量存在显著的空间差异,由东南部向西北部从4 000 mm以上降至500 mm左右[19]；降水的垂直分布规律明显,从非冰川区过渡到冰川区年降水量随高程增加的递减率为10～30 mm/(100 m)[20]；同时降水季节分配不均匀,雨季旱季分明,雨季为每年的5—10月,旱季为11月至次年的4月。雅鲁藏布江蕴含丰富的淡水资源,对西藏经济社会和下游国家发展具有十分重要的作用。

图1 雅鲁藏布江流域DEM及气象站点分布Fig.1 Distribution of meteorological stations and DEM in the Yarlung Zangbo River basin

2 研究数据与方法

2.1 研究数据

2.1.1 气象要素数据

选取中国气象局通过拟合数据序列计算并优化薄盘平滑样条函数、最终利用样条函数插值得到的1961—2014年气温和降水数据,该数据时间分辨率为1 d,空间分辨率为0.25°×0.25°。

2.1.2 下垫面数据

植被覆盖数据采用美国航天局开发的GIMMS NDVI第三代数据集(1982—2014年),空间分辨率为8 km×8 km,时间分辨率为15 d,该数据集已被证明是描述植被生长动态变化最好的数据集之一[21]；积雪数据采用MODIS网站提供的2000—2016年时间分辨率为8 d、空间分辨率为500 m的数据；雪深数据采用中国科学院寒区旱区环境与工程研究所提供的1982—2016年中国雪深长时间序列数据集[22],空间分辨率为0.25°,时间分辨率为1 d；冰川数据采用国家青藏高原科学数据中心(http:∥data.tpdc.ac.cn)提供的1976年、2001年、2013年和2017年4期数据[23];冻土数据采用Zou等[24]提供的青藏高原新绘制的冻土分布图(2017)。

2.1.3 径流数据

径流数据采用中国气象局国家气候中心提供的20个国家气象站数据和西藏水文水资源勘测局提供的雅鲁藏布江干支流上的12个水文站(拉孜、日喀则、江孜、奴各沙、唐加、旁多、拉萨、羊村、工布江达、巴河桥、更张、奴下，见图1)1961—2015年实测月流量数据。

2.2 研究方法

2.2.1 统计分析方法

本文采用线性倾向估计方法,通过计算Slope值[21]研究气温、降水、植被、积雪面积、雪深等要素的变化趋势。利用多元线性回归方法,研究冻土区域气温变化特征。

2.2.2 径流变化及归因分析方法

为深入认识雅鲁藏布江流域径流演变规律并进行归因分析,采用集中度和集中期指标分析月径流量的年内分布规律,采用Pettitt检验、滑动t检验和Mann- Kendall检验方法检验年径流量突变点[25],利用去趋势波动分析方法、小波分析方法探讨各站年径流量序列的尺度行为和周期变化。

考虑到雅鲁藏布江流域水文气象数据稀缺,径流产生的物理机制复杂等特点,本文采用Budyko框架[26]对雅鲁藏布江流域径流变化进行归因分析。由于流域径流产生机制独特,冰川径流在流域径流中占有不可忽略的比例,所以冰川径流变化在径流归因中单独考虑。下垫面变化主要受地形、土壤、植被和人类活动变化影响,长期而言地形和土壤可看作固定不变,同时由于该流域人类活动较弱,考虑忽略人类活动。因此,该流域下垫面变化主要由植被覆盖变化引起[25]。

Schaake[27]研究表明,可以通过弹性系数来量化特定变量对径流的敏感性,其定义如下:

(1)

式中:εyi为弹性系数；Rt为水文站观测的年平均总径流深，mm；yi代表特定的影响因素,即年平均降水量(P)、年平均实际蒸散发量(ET0)、冰川径流占总径流的比例(r)以及流域下垫面特征的参数(n)。

假设φ=ET0/P，弹性系数可以分别表示如下:

(2)

(3)

(4)

(5)

一旦确定了总径流量的突变点,就可以通过比较突变点前后2个时期的冰川径流量,并由Budyko框架中的公式(5)来计算冰川径流量对2个时期总径流变化的贡献。

根据弹性系数的定义,可以通过将每个变量的变化引起的径流变化相加得出总径流变化,如下所示:

O(Rt)=C(P)+C(ET0)+C(n)+C(r)+δ

(6)

式中:O(Rt)为观测到的径流变化,mm；C(P)、C(ET0)、C(n)和C(r)分别表示P、ET0、n和r所导致的径流变化,mm；δ为观测和计算的径流变化之间的偏差,mm。

每个变量对径流变化的相对贡献可以表示如下:

(7)

式中:RCyi为yi的相对贡献。

3 结果分析与讨论

3.1 气象要素演变特征

3.1.1 气温演变趋势

在全球变暖的大背景下,青藏高原地区存在明显的温度上升趋势。图2给出了1961—2014年最高气温、平均气温、最低气温在年尺度和季节尺度上的线性趋势变化空间分布。整个流域的气温是一致升高的,升温率在0.20～0.60 ℃/(10 a)之间变化,流域上游和拉萨河流域升温最快。最高气温的升温率最小,平均气温次之,最低气温的升温幅度最大。

图2 1961—2014年最高气温、平均气温、最低气温的年和季节趋势空间分布Fig.2 Distribution of the maximum temperature,mean temperature,minimum temperature and seasonal trends in 1961—2014

3.1.2 降水演变特征

图3为雅鲁藏布江流域1961—2014年降水在年尺度和季节尺度上线性趋势的空间分布。由图3可知,降水整体呈增加趋势,但上游部分区域降水呈减少趋势。春季除上游外,整个流域降水都是增加的,尤其是下游地区；夏季下游地区降水减少,与该流域冰川退缩存在较好的对应关系[28]；秋季的空间分布与年尺度类似；冬季时整个流域一致增加。

图3 1961—2014年降水的年和季节趋势空间分布Fig.3 Distribution of annual and seasonal trends of precipitation in 1961—2014

3.2 下垫面演变规律

3.2.1 植被时空变化特征

图4为1982—2014年雅鲁藏布江流域归一化植被指数(INDV)的变化趋势。由图4可知,1982—2014年INDV变化范围为0.25～0.28,INDV呈不显著波动上升趋势,上升速率为0.002/(10 a),表明近33 a间雅鲁藏布江流域植被整体不断改善,植被覆盖度逐步提高。

图4 1982—2014年流域INDV的变化趋势Fig.4 Variation trend of INDV during the period of 1982—2014

图5(a)为雅鲁藏布江流域归一化植被指数多年平均分布图,由图5可知,下游植被生长状况较中上游好(下游INDV高于中上游),INDV最大值可达0.90,而中上游INDV值整体分布在0～0.46之间。图5(b)反映了流域年尺度INDV的变化情况,统计得到Slope大于0的网格占比为63.7%,表明整体上流域归一化植被指数得到改善,植被覆盖率不断提高。图5(c)统计了通过0.05显著性检验的INDV趋势值,中游流域INDV呈显著上升趋势,上游流域西北部和下游流域INDV呈显著下降趋势,表明中游流域植被得到改善,上游流域西北部和下游流域植被存在退化现象。这与刘浏等[21]对雅鲁藏布江流域归一化植被指数研究的成果较为一致。

图5 1982—2014年雅鲁藏布江流域INDV分布Fig.5 Distribution of INDV in Yarlung Zangbo River basin during the period of 1982—2014

3.2.2 积雪时空变化特征

图6 年内积雪覆盖率Fig.6 Map of annual snow cover fraction

对2000—2016年雅鲁藏布江流域积雪覆盖率年内变化进行统计分析。如图6所示,5—9月流域积雪覆盖率较小,1—3月流域积雪覆盖率较大。2000—2016年流域年内变化呈现出夏季小、冬春大的特点。最大积雪覆盖率、平均积雪覆盖率、最小积雪覆盖率年内分布均值分别为49.12%、28.70%和12.41%。最小积雪覆盖率最小值为3.12%,出现在7月；最大积雪覆盖率最大值为65.29%,出现在3月。

图7(a)为雅鲁藏布江流域积雪覆盖率分布。由图7(a)可知,下游北部嘉黎附近积雪覆盖率较大,均值达到70%左右；下游河谷至流域中游北部山区积雪覆盖率均值逐步减小,均值为10%～20%；中游雅鲁藏布江干流沿线积雪覆盖率最小,在3%以内。这与Ban等[29]研究得到积雪覆盖率在雅鲁藏布江流域的东北部相对高于中部和南部区域的研究成果一致。

图7 雅鲁藏布江流域积雪覆盖率及雪深变化率分布Fig.7 Distribution map of snow cover fraction and variation of snow cover depth in Yarlung Zangbo River basin

图7(b)为雅鲁藏布江流域雪深变化率分布。上游雪深整体以增加为主,部分区域雪深变幅为-1.2～-0.5 m/(10 a)；中下游雪深变化率较上游大,变幅为-1.2～0.8 m/(10 a),中游北部和下游东北部雪深以减少趋势为主,最大可达-1.2 m/(10 a)。

3.2.3 冰川、冻土的演变特征

通过收集雅鲁藏布江流域1976年、2001年、2013年和2017年冰川数据研究流域近40 a冰川变化。结果显示,近40 a随着气温升高,冰川呈现快速退缩的趋势,变幅为-422.9 km2/(10 a),尤其是2013—2017年冰川面积锐减了1 600 km2。随着冰川面积的加速退缩,在未来一段时间内冰川将不断减少乃至消失,这将给流域水资源的开发、利用和配置带来巨大挑战。

为了探究雅鲁藏布江流域等温线与冻土分布的关系,同时由于流域内气象站点稀少,本文参考焦世晖等[30]在青藏高原的研究方法,基于流域及周围20个国家气象站1980—2017年逐日气温数据,采用Excel对流域20个气象站点38 a的年平均气温与海拔、纬度之间的关系进行多元线性回归分析,得到

Y=39.875 59-0.005 2e-0.489 71l

(8)

式中:Y为多年平均气温, ℃；e为海拔,m；l为纬度, °N。

多元线性回归结果的R2=0.7,且该结果通过α为0.05的F值检验,结果较为可信。因此,根据海拔和纬度数值便可利用式(1)计算出某点的平均温度。采用此公式,在流域上游及中游站点稀少的位置进行补点。基于补点和站点气温数据,利用地理分析软件的自然邻域法进行空间插值,计算得到流域年平均气温分布，见图8。由图8可知,季节冻土分布较广,在流域的上中游均有分布,多年冻土主要分布在流域的上游前端和流域的东北部,非冻土主要分布在流域的东南部。流域上游前端的多年冻土年平均气温约-2 ℃,流域东北部的多年冻土年平均气温约为-2～4 ℃,这一研究结果与前人研究结论一致[30]。

图8 雅鲁藏布江流域冻土分布及绘制的多年平均气温等温线Fig.8 Distribution of frozen soil and the drawn multi- year average temperature isotherm in the Yarlung Zangbo River basin

3.3 径流变化及其归因分析

3.3.1 径流变化特征

(1) 年内变化特征。雅鲁藏布江流域奴下、羊村等12个水文站径流的年内分布规律结果表明,月平均径流的集中度值为0.18～0.72(图9(a)),主要集中在7—9月(图9(b)),且在空间上存在一定差异。一般来说,5—9月是雅鲁藏布江流域的雨季,可以看出干流的雨季径流量占年径流总量的73%左右,空间差异不大。黄俊雄等[31]基于雅鲁藏布江奴各沙、羊村、奴下3个水文站1956—1995年逐月流量数据,研究得到雅鲁藏布江流域6—9月径流分配最丰,径流量占年径流总量的65%以上，本文研究结果与其较为一致。流域内3个子流域(尼洋河、拉萨河和年楚河)的雨季径流总量占其年径流总量的比例分别为68%、79%和83%,具有一定的空间差异。

整体上流域径流具有相似的年内分布规律,即主要集中在雨季。然而,从下游到上游,径流高值在流域西部仅集中在8月份,在中部则主要在7—8月,而在东部地区多为6—9月,最大径流出现的月份逐渐延后且持续时间变短。

(2) 年际变化特征。通过对径流进行周期分析,得到干流年径流量存在3～4 a和12～15 a左右的不显著周期变化,部分支流(拉萨、工布江达和日喀则水文站)为20 a左右周期,同时各站均存在一个约30～33 a的长周期且通过了显著性检验。利用去趋势波动分析(DFA)方法分析了12个水文站年径流序列的尺度行为，除工布江达和唐加水文站的α<0.53外,其余10个站点的α最小值为更张站年径流序列的0.745,最大值为江孜站的0.963,表明雅鲁藏布江干流及其主要支流的年径流过程具有明显的长程相关性,并在95%置信水平下表现为正相关性。

图9 各水文站月径流序列Fig.9 Monthly runoff series of each hydrological station

利用小波分析定义的各站年平均径流序列与降水序列的平均周期变化显示:降水序列具有3～4 a、16 a和32 a的周期,年径流序列相应具有3～4 a、12 a、20 a和32 a左右的周期,进一步表明径流的年际变化规律与降水变化具有相似性。

3.3.2 径流变化的量化归因分析

根据雅鲁藏布江干流上设置的3座水文站(奴各沙、羊村、奴下)的地理位置对雅鲁藏布江流域进行子流域划分,由此得到了4个子流域。但由于奴下站以东的区域缺少相关水文数据,无法准确有效地进行径流归因分析,因此本文只选取奴下站以西的3个子流域作为研究区域,这3个子流域由西至东依次被视为雅鲁藏布江流域的上游、中游以及下游(图10)。每个子流域中的径流量均由该子流域上下游两端的控制水文站决定:上游的径流量为奴各沙站实测径流量,中游的径流量为羊村站与奴各沙站径流量的差值,下游的径流量则为奴下站与羊村站径流量的差值。

图10 雅鲁藏布江水文站分布及子流域划分Fig.10 Distribution of the hydrological stations and division of the sub- basins in the Yarlung Zangbo River basin

结合中上游地区径流量的突变检验结果(表1),为了统一雅鲁藏布江流域3个分区的气候变化和下垫面变化对径流变化的影响时期,选择1997年作为突变点,以确定本文的基准期和变化期,即基准期(Ⅰ期)为1966—1997年,变化期(Ⅱ期)为1998—2015年。本文与张建云等[32]研究雅鲁藏布江奴下站实测年径流得到的突变点分别发生在1966年和1998年一致。

表1 通过3种统计方法检测到的变化点

图11表明,气候、冰川径流和下垫面的变化对上游、中游和下游的径流变化表现出不同影响。与Ⅰ期相比,Ⅱ期上游、中游、下游的径流深分别增加了29.08 mm、85.80 mm和20.67 mm。为了评估调整后的Budyko框架在雅鲁藏布江流域的适用性,本文对计算的径流量变化(C(Rt))与观测的径流量变化(O(Rt))之间的相对误差进行了计算,从图11中可以看到,所有3个子流域相对误差均低于0.4%,表明本文所提出的归因分析方法对于量化各因素对径流变化的影响是有效和可信的[25]。在上游,由降水变化(ΔP)、潜在蒸散量变化(ΔET0)、下垫面变化(Δn)和冰川径流变化(Δr)引起的径流增加的贡献分别为-0.81 mm、-1.54 mm、14.34 mm和14.30 mm,对应的贡献率分别为-2.78%、-5.30%、59.61%和49.18%；以上影响因素对中游径流增加的相应贡献分别为33.77 mm、-1.16 mm、31.24 mm和22.34 mm,对应的贡献率分别为39.36%、-1.36%、36.4%和26.04%；而在下游,ΔP、ΔET0、Δn和Δr对径流增加的贡献分别为26.71 mm、-0.04 mm、-12.14 mm和6.16 mm,对应贡献率分别为129.21%、-0.17%、-55.74%和29.8%。

图11 气候变化、冰川径流和下垫面变化对径流变化的贡献Fig.11 Contribution of climate change,glacier runoff and underlying surface change to runoff variation

以上结果表明,以n表示的下垫面变化以及以r为代表的冰川径流变化对上游径流的增加起着主导作用；而降水量的变化对中下游径流的增加影响最大,同时冰川径流也有着超过25%的不可忽略的贡献率。此外,本文还对全流域径流量进行了归因分析,ΔP、ΔET0、Δn和Δr对整个流域径流量增加的贡献率分别为39.62%、-2.74%、32.32%和30.94%,进一步说明了研究下垫面变化和冰川径流变化对径流变化的重要性。

目前,在雅鲁藏布江流域已开展许多径流归因分析工作,本文与以往研究的不同之处在于以往研究多集中在中游和奴下站集水范围内,且作为一个整体进行考虑[33]；本文则尝试把奴下站集水区域划分为上游、中游、下游3个子流域,深入分析了各子流域径流变化的原因,并量化降水变化、潜在蒸散量变化、下垫面变化和冰川变化对径流变化的贡献。本文与Wang等[13]采用Budyko方法进行径流归因得到的结果偏差较大的原因可能是:① 所采取得Budyko方法不同,例如在计算冰川径流时本文采用流域水量平衡计算,Wang等[13]采用度日因子模型；② 采取的数据不同，不同数据由于融合或插值方法的不同,也会引起数据质量具有一定的差异,最终也会使研究结果出现偏差。

4 结 论

本文在分析气象要素和下垫面演变特征的基础上,研究变化环境下雅鲁藏布江流域的径流变化规律及其归因分析,得到以下几点结论:

(1) 1961—2014年,整个流域的气温是一致升高的,升温率在0.20～0.60 ℃/(10 a)之间变化,降水整体呈增加趋势,但上游降水呈减少趋势。过去50 a来,流域气候向暖湿化方向发展。

(2) 流域植被整体不断改善,植被覆盖度逐步提高。中游流域植被得到改善,上游流域西北部和下游流域植被存在退化现象。上游雪深整体以增加为主,中游北部和下游东北部雪深以减少趋势为主。冰川呈快速退缩的趋势,变幅为-422.9 km2/(10 a)。季节冻土在上中游均有分布,多年冻土主要分布在上游前端和流域的东北部,非冻土主要分布在流域的东南部。上游前端的多年冻土年平均气温约-2 ℃,流域东北部的多年冻土年平均气温约为-2～4 ℃。

(3) 整体上流域径流具有相似的年内分布规律,主要集中在雨季。降水序列具有3～4 a、16 a和32 a的周期,年径流序列相应具有3～4 a、12 a、20 a和32 a左右的周期,径流的年际变化规律与降水变化具有相似性。奴下站集水范围内降水变化、潜在蒸散量变化、下垫面变化和冰川变化对径流量增加的贡献率分别为39.62%、-2.74%、32.32%和30.94%。

未来在高寒流域开展相关研究中,应加强以下几个方面的探索:① 加强气温、降水、冰川径流、融雪径流、冻土水文等实测数据的观测工作,创新数据融合算法,提高融合数据的质量；② 加强流域尺度冰冻圈全要素各过程的综合模拟,深入理解各要素在流域径流变化中的影响,深化对冰冻圈变化机理与过程的科学认识。

致谢:本文主要成果来源于项目研究报告,作者对所有项目参加人员表示衷心的感谢。