雅鲁藏布江流域陆地水储量变化特征及归因

刘 浏，王宣宣，牛乾坤，伦玉蕊，程湫雅，程 磊，徐宗学

（1.中国农业大学水利与土木工程学院，北京 100083；2.中国农业大学中国农业水问题研究中心，北京 100083；3.武汉大学水利水电学院，武汉 430072；4.北京师范大学水科学研究院，北京 100875；5.城市水循环与海绵城市技术北京市重点实验室，北京 100875）

0 引 言

陆地水储量（Terrestrial Water Storage，TWS）是指储存在地表或地表以下的各种形式的水，包括地表水、地下水、土壤水、冰川、湖泊、林冠截留水和生物质水等[1]，在陆地和全球水循环中扮演着关键的角色，它的组成成分对水文、气候系统有着重要影响[2]。作为水资源平衡的重要组成部分，陆地水储量变化（Terrestrial Water Storage Change，TWSC）是指陆地水储量（TWS）在时间和空间上的演变过程[3]。陆地水储量及其变化在评估和预测气候变化、农业气象、洪水和干旱等极端事件中起着不可替代作用[4]。

青藏高原位于中国西南部、亚洲中部，是世界上海拔最高、面积最大的高原，面积约 250 万 km2，有地球的“第三极”之称[5]。青藏高原地区是中国境内众多主要河流的发源地，湖泊面积约占中国湖泊总面积的一半左右，因此又被称为“亚洲水塔”[6]。青藏高原作为全球气候变化最为敏感的地区之一，直接关系到中国西南地区的水资源利用和流域的用水安全状况[7]。近几十年来，相对于全球平均温度每10 a升高0.17 ℃，青藏高原的平均温度每10 a升高0.35 ℃左右，其增温幅度远高于全球平均水平[8]。青藏高原地区的温度升高导致了冰川融化速度加快、径流不断增加等一系列的失衡现象，导致其水循环过程随之发生巨大的转变[9]。冰川体积不断缩小导致青藏高原地区的固态水储量减少[10]，与冰川的变化相反，青藏高原地区有超过三分之二的湖泊在扩张，湖泊的水量和面积整体都呈现出明显的增加趋势[11]。同时，青藏高原地区的径流也呈现出不同程度的增加，进而导致了青藏高原及其周边地区的水灾风险加大。开展青藏高原陆地水储量变化规律的研究关乎到中国以及东南亚各国的用水安全，有助于进一步加强对高寒地带缺资料地区的水循环与水资源演变规律的了解，同时为科学合理开发利用青藏高原地区的水资源提供参考[12]。

青藏高原地区缺乏高精度的水文气象实测资料[13]，限制了对其陆地水储量的定量研究。2002年3月发射的GRACE（Gravity Recovery and Climate Experiment）重力卫星能够提供较高精度的水文信号，为水储量的反演开辟了新的途径[14]。然而，GRACE数据反演的水储量的时间序列较短，不足20 a，难以对青藏高原地区过去几十年来在全球气候变暖的背景下水循环加强的进行全面深刻的解析。因此，本研究选取青藏高原东南部的雅鲁藏布江流域为研究区，基于2003—2017年的GRACE重力卫星数据、ERA5（the fifth-generation reanalysis product of the European Centre for Medium Range Weather Forecasts）再分析数据和 GLDAS（Global Land Data Assimilation System）数据，采用水量平衡法和组分相加法两种方法，将反演的4种陆地水储量结果与GRACE重力卫星数据反演的陆地水储量在时间和空间尺度上进行适应性评估，获取陆地水储量反演的较优方法，并基于该方法系统探究雅鲁藏布江流域长序列（1948—2017年）陆地水储量的时空演变规律及主导因素，以期为探究资料匮乏地区的生态水文过程提供可靠参考。

1 数据与方法

1.1 研究区概况

雅鲁藏布江流域位于 27°N～32°N、82°E～97°E 之间，地处青藏高原东南部、喜马拉雅山脉北麓，流域面积约24 万km2，干流长度超过2 000 km[15]，流域地势西高东低，海拔垂直差异大[16]。雅鲁藏布江流域水资源储藏丰富[17]，流域水资源状况关乎到中国西藏以及下游地区的印度、孟加拉国等国家的用水的安全。上游地区平均海拔4 600 m以上，以高山谷地等地貌为主[18]；中游地区主要为属于藏南谷地，主要为温暖半干旱气候[19]，是人类活动集中的地区；下游地区包括高原温带湿润、半湿润气候区和山地亚热带湿润气候区，气候相对温暖湿润[20]。雅鲁藏布江是世界上海拔最高的河流之一，受特殊的地形条件和海拔高度的影响，流域呈现出对气候变化的响应较为敏感和脆弱的特征[21]。尤其是近几十年来，雅鲁藏布江流域冻土不断融化[22]，活动层增厚[23]，冰川退缩幅度明显[24]，这些都对流域生态环境和水资源保护产生了巨大影响。因此，准确反演雅鲁藏布江流域陆地水储量变化，有助于加深在全球气候变暖背景下对流域水文、生态系统变化对人类活动以及气候变化响应机理的理解。研究区域概况如图1所示。

1.2 数据来源

本研究使用的两种数据集分别是欧洲中期天气预报中心第五代再分析资料ERA5和美国的全球陆面数据同化模型GLDAS数据集。利用这两种数据集提供的降水、蒸散发、径流、土壤含水量和雪水当量来估算TWSC，再将得到的TWSC与GRACE数据反演的TWSC进行比较。

1.2.1 GRACE数据

GRACE重力卫星是由美国宇航局（National Aeronautics and Space Administration，NASA）和德国航天局（German Aerospace Center，DLR）合作研发的双星系统，于2002年3月成功发射[25]。目前有德克萨斯大学空间研究中心（Center for Space Research，CSR）、德国地学研究中心（German Research Centre for Geosciences，GFZ）以及美国喷气动力实验室（Jet Propulsion Laboratory，JPL）提供较为成熟的GRACE数据产品[26]。本研究使用的是JPL提供的 2003—2017年的月数据产品（https://grace.jpl.nasa.gov/data/get-data/ jpl_global_mascons/），具体版本为 JPL RL06M Version 2.0，空间分辨率为 0.5°×0.5°。

1.2.2 ERA5再分析数据

本研究使用的 ERA5数据是来自欧洲中期天气预报中心（European Centre for Medium Range Weather Forecasts，ECMWF）最新发布的第五代全球大气再分析资料 ERA5-Land的月平均数据（https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysi s-era5-land-monthly-means?tab=form），其空间分辨率较之前的版本提升至0.1°×0.1°，精度相对ERA-interim亦有所改善[27]。

1.2.3 GLDAS数据

本研究使用的 GLDAS数据是美国国家航空航天局（NASA）戈达德空间飞行中心（Goddard Space Flight Center，GSFC）和美国海洋和大气局（National Oceanic and Atmospheric Administration，NOAA）国家环境预报中心（National Centers for Environmental Prediction，NCEP）联合发布的基于卫星、陆面模式和地面观测数据的同化产品，能够向用户提供多种驱动数据[28]。目前，GLDAS结合MOSAIC、NOAH（NCEP，OSU，Air Force and Office of Hydrology）、CLM（Community Land Model）、VIC（Variable Infiltration Capacity）4个陆面模式提供了大量的陆面数据产品。本研究采用的 GLDAS数据版本为 GLDAS_NOAH025_M（https://disc.gsfc.nasa.gov/datasets? keywords=GLDAS&page=1&temporalResolution=1%20mont）。与其他GLDAS产品相比，GLDAS-NOAH的空间分辨率更高，时间跨度更长，在青藏高原地区的应用较为广泛[29-31]。

1.3 研究方法

本研究采用水量平衡法和组分相加法两种方法，基于ERA5和GLDAS两种数据集，分别计算雅鲁藏布江流域的TWSC，共生成4套TWSC数据集，并基于GRACE数据，定量评估不同数据源、不同方法对雅鲁藏布江流域水储量估算的适用性。为了便于比较，ERA5和GLDAS数据均采用双线性插值，插值为 0.5°×0.5°，与 GRACE数据的空间分辨率保持一致。

1.3.1 水量平衡法

水量平衡法（PER）根据水量平衡方程[32]，用降水（Precipitation）、蒸散发（Evapotranspiration）、径流（Runoff）来计算TWSC，公式如下：

式中S(t)为t时刻的TWS，mm；P(t)为总降水量，mm；E(t)为蒸散发量，mm；R(t)为该流域的径流量，mm。

1.3.2 组分相加法

TWSC也可以用其主要的组成成分的变化来表示[33]，大量研究表明，在全球尺度和流域尺度上，用土壤含水量和雪水当量的变化估算的TWSC与GRACE数据反演的TWSC有很好的一致性[34-36]。下式是组分相加法（SS）用土壤含水量（Soil moisture）和雪水当量（Snow water equivalent）的变化近似表示TWSC：

式中So(t)和Sn(t)分别为某一区域和时间内的土壤含水量和雪水当量，mm。

1.3.3 评估指标

本研究采用线性回归方程来衡量水储量的变化。

式中y代表水储量的变化，x代表时间，a代表斜率，b代表截距[37]。

采用皮尔逊相关系数来衡量不同变量之间的相关性。

式中xi和yi分别为两组序列的数据，n是变量的序列数目，和分别为两组数列的平均值，rxy为两个变量之间的相关系数[38]。

Mann-Kendall检验方法是是世界气象组织推荐的非参数检验方法，能够客观反映时间序列趋势，目前已经被广泛用于气候参数和水文序列的分析中[39]。Mann-Kendall方法通过给定一个时间序列Kt=(k1,k2, …,km)，计算统计量Zc，计算公式如下：

式中kp和kq分别表示第p个和第q个样本，1＜p＜q＜m，m为样本容量；D(s)表示s的方差；β表示变化趋势，β＞0表示增加趋势，β＜0表示减少趋势。|Zc|＞Z(1-α/2)表示在α显著性水平下，变化趋势显著。

2 结果与分析

2.1 陆地水储量反演结果的适应性评估

2.1.1 陆地水储量时间特征对比分析

图2是基于ERA5和GLDAS数据，分别采用水量平衡法和组分相加法估算的2003—2017年雅鲁藏布江流域TWSC逐月变化及其与GRACE数据反演结果的对比图。如图2a所示，GRACE数据反演的雅鲁藏布江流域的TWSC在 2003—2017年间以 0.90 mm/月的速度减少（P＜0.01），相当于 2.46 亿 t/a。不同数据源以及不同方法反演的陆地水储量存在较大差异性，PER（ERA5）、SS（ERA5）和 PER（GLDAS）的结果均呈现周期波动的状态，无显著的增减趋势。4种结果中，SS（GLDAS）以0.397 mm/月的速度减少（P＜0.01），在月尺度变化趋势上与GRACE数据反演的结果最为一致。

如图2b所示，PER（ERA5）和SS（ERA5）变化幅度分别为-44.77～77.15和-50.04～71.81 mm，PER（GLDAS）和SS（GLDAS）范围分别为-62.22～82.97和-69.54～113.36 mm，均与GRACE数据反演的TWSC变幅（-159.23～208.19 mm）存在一定差异，但相较于其他3套结果，SS（GLDAS）反演的水储量变幅与 GRACE结果更为接近。此外，本研究进一步检验了 4套反演结果2003—2017年间180个月的TWSC与GRACE数据反演结果的相关性，PER（ERA5）、SS（ERA5）、PER（GLDAS）和SS（GLDAS）与GRACE结果的相关系数分别为-0.03、0.36、-0.12和0.75。由此看出，组分相加法估算的 TWSC的逐月变化与 GRACE数据反演的TWSC逐月变化的相关性均优于水量平衡法。此外，无论是水量平衡法还是组分相加法，基于GLDAS数据反演的月尺度TWSC均与GRACE的反演结果的相关性更高。综上所述，SS（GLDAS）在月尺度变化趋势和变化幅度上，与GRACE数据反演的结果更为一致。

2.1.2 陆地水储量空间特征对比分析

图3是基于 ERA5和 GLDAS数据的 TWSC与GRACE数据反演结果相关系数的空间分布。基于ERA5数据，采用水量平衡法（PER）和组分相加法（SS）估算的上游和中游地区的相关系数较大，最大值分别为0.37和 0.67，而相关系数较小的地区主要集中在下游的尼洋河流域。基于GLDAS数据，采用水量平衡法估算的大部分地区的相关系数较小，仅流域东南部的相关系数略高，最大值为仅 0.21；采用组分相加法估算的大部分地区的相关系数相对较大，大于0.50的面积超过61.3%，最大值为 0.82，而相关系数较小的地区主要集中在下游的尼洋河流域。由此看出，陆地水储量反演结果较差的地区主要集中在下游的尼洋河流域，即帕隆藏布冰川分布地区[40]，该地区雪线高度的高程变化超过2 000 m[41]，受特殊的海拔和地形的影响，模型存在较大的测量误差，导致 GLDAS数据估算的结果难以准确反映该地区的TWSC。与时间变化特征对比分析结果类似，SS（GLDAS）与GRACE数据反演结果的空间相关性最好。

综合两套数据、两种算法得出的四组TWSC的时空变化的结果表明：无论是基于ERA5数据还是GLDAS数据，采用组分相加法估算的TWSC与GRACE数据反演的TWSC在时间和空间上的一致性均优于水量平衡法。其中，SS（GLDAS）与GRACE重力卫星数据反演的结果最为一致。

2.1.3 组分相加法敏感性分析

为探究不同组分对 TWSC反演结果的影响，本研究基于 GLDAS数据和组分相加法进行了TWSC反演的敏感性分析。分别考虑以下组合：1）土壤含水量、雪水当量、地表水、地下水和植被冠层截留（图4a）；2）土壤含水量、雪水当量和地下水（图4b）；3）土壤含水量和雪水当量（图4c）。三种组合得到的相关系数空间分布特征相似，下游的尼洋河流域相关性较差，统计结果如表1所示。考虑土壤含水量、雪水当量、地表水、地下水和植被冠层截留时，SS（GLDAS）与GRACE数据反演的流域平均 TWSC的相关系数为0.50，48.4%的地区相关系数大于 0.50。考虑土壤含水量、雪水当量和地下水时，SS（GLDAS）与 GRACE数据反演的流域平均TWSC的相关系数为0.51，58.0%的地区相关系数大于 0.50。仅考虑土壤含水量和雪水当量时，SS（GLDAS）与 GRACE数据反演的流域平均TWSC的相关系数最大，为0.53，且该条件下，61.3%的地区相关系数大于 0.50。表明综合考虑地表水、地下水和植物冠层截留时反演的 TWSC结果与 GRACE反演结果的一致性并未有所提升。此外，已有研究表明，在雅鲁藏布江流域，相对于其他组分，土壤含水量和雪水当量的变化更加显著[10,34,37]。因此，基于土壤含水量和雪水当量，采用组分相加法反演雅鲁藏布江流域陆地水储量的变化是合理的。

表1 组分相加法敏感性分析结果对比Table 1 Comparison of sensitivity analysis results of the summation method

因此，基于GLDAS数据的组分相加法更适用于雅鲁藏布江流域陆地水储量的反演，可为准确描述高寒地区水储量变化提供重要参考。

2.2 1948—2017年雅鲁藏布江流域陆地水储量变化特征

鉴于 GLDAS数据的组分相加法反演的 TWSC与GRACE重力卫星数据反演的结果最为一致，本研究基于1948—2017年的GLDAS数据，采用组分相加法进一步探究 1948—2017年雅鲁藏布江流域陆地水储量的变化特征。

2.2.1 1948—2017年TWSC的时间变化特征

1948—2017 年间 TWSC呈现极显著上升趋势（P＜0.01），但TWSC在2002年前后发生突变。如图5所示，2002年之前TWSC以0.024 mm/月的速度呈现极显著上升的趋势（P＜0.01），2002年之后TWSC以0.397 mm/月的速度呈现极显著下降的趋势（P＜0.01）。这与刘浏等[42]在雅鲁藏布江流域干湿转换特征研究中的结果保持一致，即1982—2015年雅鲁藏布江流域总体呈现湿润化趋势，但进入21世纪以后流域表现出极显著的干旱化趋势。

本研究通过小波分析探究1948—2017年雅鲁藏布江流域TWSC的长时间的变化特征，如图6所示。从大尺度周期来看，对应于时频分布图的上部，TWSC存在着40 a左右的大尺度震荡周期，但鉴于时间序列总长度为70 a，因此对40 a左右的大尺度周期不予考虑。从中尺度周期来看，对应于时频分布图的中部，TWSC存在24 a左右的中尺度振荡周期，经历3次高低起伏。从小尺度周期来看，对应于时频分布图的下部，TWSC在1948—1975年期间存在着5 a左右的小尺度振荡周期，在1976—2017年期间存在着8 a左右的小尺度振荡周期，TWSC随着时间尺度的减小表现出更为复杂的周期振荡。小波方差图中显示，在1948—2017年期间，陆地水储量的年变化的主周期为35 a。从中尺度和小尺度来看，TWSC在2002年左右都发生了一次由大到小的突变。

2.2.2 1948—2017年TWSC的空间变化特征

1948—2017 年雅鲁藏布江流域TWSC从上游到下游呈现递增的规律，如图7a所示，陆地水储量最大减少值为-10.26 mm，位于雅鲁藏布江上游西部，最大增加值为38.36 mm，位于下游东南部，中游和地区的陆地水储量基本保持不变或增加。因此，1948—2017年间雅鲁藏布江流域 TWSC呈现显著增加趋势主要与中下游地区的陆地水储量增加有关。如图7b，流域大部分地区的TWSC变化整体呈增加趋势，但在中游的“一江两河”（雅鲁藏布江中游及其支流拉萨河和年楚河）地区和下游冰川分布较多的帕隆藏布地区的多年平均TWSC变化呈现减少的趋势，进一步表明人类活动和冰川积雪集中地带的冰雪消融对水储量变化趋势有着不可忽略的影响。

图8是2002年前后雅鲁藏布江流域多年平均TWSC和变化趋势的空间分布图。在1948—2002年，相对于上游地区，中下游地区的陆地水储量较为丰富，但是中下游大部分地区TWSC呈现减少的趋势，表明虽然中下游地区陆地水储量在1948—2002年间有所增加，但增加的趋势却在变缓。在2003—2017年，流域陆地水储量减少的地区主要分布在除下游东南部以外的地区，占流域总面积的83.7%。陆地水储量增加的地区主要集中在下游东南部，但是增加的趋势也明显变缓。2002年前后，中游地区的陆地水储量转变最为明显，陆地水储量由增加变为减少。这可能是由于近些年来中游地区人类活动强度较大，城市建设用地增加、过度放牧导致植被覆盖面积减少，降低了该地区的持水能力，进而导致了该地区的陆地水储量急剧减少。

2.3 TWSC变化的归因分析

本研究基于1948—2017年的GLDAS数据进行70 a TWSC变化的归因分析，图9是1948—2017年土壤含水量和雪水当量的逐月变化。在2002年前后，土壤水变化和雪水当量的变化趋势与TWSC的变化趋势基本一致。1948—2002年土壤含水量的变化以 0.014 mm/月的速率呈现出极显著上升的趋势（P＜0.01），2003—2017年土壤含水量的变化以 0.394 mm/月的速率呈现出极显著下降的趋势（P＜0.01）。1948—2002年雪水当量的变化以0.009 mm/月的速率呈现出极显著上升的趋势（P＜0.01），2003—2017年雪水当量的变化以0.002 mm/月的速率呈现出极显著下降的趋势（P＜0.01）。根据Meng等[43]计算不同组分对青藏高原陆地水储量变化贡献率的方法，本研究计算2002年前后土壤含水量和雪水当量的变化对TWSC的贡献率。1948—2002年，土壤含水量和雪水当量的变化对TWSC的贡献率分别为61%和39%；2003—2017年，土壤含水量和雪水当量的变化对TWSC的贡献率分别为99%和1%。相较于雪水当量，2002年之后土壤含水量的变化速率更为显著，这可能是由于 2002年以后，青藏高原地区气温增加，帕隆藏布地区冰川和积雪消融转化为土壤水[44]，多年冻土向季节性冻土转化，导致土壤含水量的变化增加[45]，土壤含水量的变化对TWSC的贡献率显著增加。同时，2002年以后流域积雪消融，积雪的面积分布和深度都有所减少，导致雪水当量的变化对TWSC的贡献率急剧减少[46]。

图10是1948—2017年TWSC、土壤含水量变化和雪水当量变化量的空间分布。可以看出，土壤含水量和TWSC的空间分布规律更加一致，2002年以后，TWSC和土壤含水量减少的区域都在从上游向下游逐渐扩张，其中中游地区变化最为明显。中游地区的TWSC和土壤含水量在1948—2002年间基本保持不变，在2003—2017年间有所减少。1948—2002年和2003—2017年土壤含水量减少的地区主要分布在上游，土壤含水量增加的地区主要集中在下游，这可能是由于下游地区由于海拔相对较低且温度较高，植被覆盖密度相对较大[37]，促进了下游地区土壤的持水能力[47]。1948—2002年雪水当量减少的地区主要集中在上游地区，中下游大部分地区主要呈现增加的趋势。2003—2017年雅鲁藏布江流域大部分地区的雪水当量与1948—2002年的雪水当量空间分布整体相似，但下游的帕隆藏布地区的雪水当量却明显减少，这再次表明了气候变暖对冰雪消融的促进作用[10]。

综上所述，在1948—2002年和2003—2017年期间，土壤含水量变化和雪水当量的变化的趋势与TWSC时间和空间上的变化趋势基本一致。在时间上，土壤含水量变化和雪水当量的变化都在2002年前后呈现出相反的变化趋势。在空间分布上，土壤含水量变化和雪水当量的变化呈现出明显的空间异质性，主要体现在人类活动强度较高的“一江两河”地区和冰川分布集中的帕隆藏布地区。但是，相对于雪水当量的变化，土壤含水量的变化对TWSC的变化起主导作用。

3 结 论

本研究利用 2003—2017年 180个月的 ERA5（the fifth-generation reanalysis product of the European Centre for Medium Range Weather Forecasts）和GLDAS（Global Land Data Assimilation System）的月数据，采用水量平衡法和组分相加法计算雅鲁藏布江流域 TWSC（Terrestrial Water Storage Change）变化。基于 GRACE（Gravity Recovery and Climate Experiment）数据反演的陆地水储量，进行不同数据源、不同计算方法在雅鲁藏布江流域的适应性评估。基于GLDAS数据的组分相加法进一步研究1948—2017年的雅鲁藏布江流域TWSC时空演变规律及其主导因素。主要结论如下：

1）综合两套数据、两种算法得出的四组TWSC的时空变化的结果表明：基于GLDAS数据的组分相加法反演的TWSC与GRACE数据反演的结果最为一致，该方法能够有效识别雅鲁藏布江流域陆地水储量变化状况，对进一步探究高寒地区的水-植被-生态复杂互馈机理具有重要意义。

2）相对于水量平衡法，组分相加法估算的雅鲁藏布江流域TWSC效果更好。但是组分相加法未考虑地下水的变化，因此基于 GLDAS数据的组分相加法估算的TWSC与 GRACE数据反演的结果存在一定的差异。例如，基于GLDAS数据的组分相加法估算的TWSC的变化幅度（-69.54～113.36 mm）比 GRACE数据反演的TWSC变化幅度（-159.23～208.19 mm）小，且二者之间存在一定的时空差异。

3）Mann-Kendall检验和小波分析的结果表明，TWSC在2002年发生了突变。且1948—2017年间TWSC、土壤含水量变化和雪水当量变化都在2002年前后呈现出相反的变化趋势，即2002年之前呈现增加的趋势，2002年之后呈现减少的趋势。在空间分布上，TWSC、土壤含水量变化和雪水当量变化都存在明显的空间异质性，主要体现在人类活动强度较高的“一江两河”地区和冰川分布集中的帕隆藏布地区。2002年前后，土壤含水量的变化对TWSC的贡献率分别为61%和99%，表明相对于雪水当量的变化，雅鲁藏布江流域土壤含水量变化对TWSC起主导作用。