龙 笛，李雪莹，李兴东，韩鹏飞，赵凡玉，洪仲坤，王一鸣，田富强

(清华大学水沙科学与水利水电工程国家重点实验室，北京 100084)

青藏高原被誉为“亚洲水塔”、“地球第三极”，平均海拔约4 000 m，是亚洲主要大河的发源地，包括亚洲重要的跨境河流印度河、恒河、雅鲁藏布江-布拉马普特拉河、怒江-萨尔温江、澜沧江-湄公河以及中国长江和黄河等。同时，青藏高原在亚洲气候系统的稳定和高原生态系统(山地森林、高寒草甸、高寒草原、农田等)的碳收支动态变化方面具有重要的生态安全屏障作用，是亚洲乃至北半球环境变化的调控器[1]。以南亚季风及西风带为主导的大气环流[2]和青藏高原独特的高山地形，共同形成了该地区丰富的水资源储量，区域内湖泊、冰川众多，湖泊群总面积约5万km2，占中国湖泊面积的一半以上，其中大于1 km2的湖泊数量接近1 400个[3-4]，是世界上湖泊分布最密集的地区之一；
共有冰川约7.7万条，总面积达8.3万km2，是除两极和格陵兰岛之外冰川数目最多的地区[5]。其中，藏东南地区是中国最主要的海洋性冰川集中区，7 700余条冰川的面积超过7 200 km2。此外，青藏高原冬、春两季的季节性积雪覆盖范围广、深度较深、持续时间长，是中国三大稳定积雪区之一[6-7]。青藏高原各类水储量(包括湖泊、冰川、积雪、冻土等)变化不仅与区域能量、大气循环密切相关，而且对亚洲主要大河的冰雪融水补给有重要影响，进一步对泛第三极地区的20 多个国家和30 多亿人口的生存发展具有重要意义[8]。

气候变化对青藏高原具有显著而深刻的影响。研究表明，青藏高原是全球范围内对气候变化最敏感的地区之一[9-10]。总体而言，近年来青藏高原正经历显著的变暖变湿趋势[11-13]，持续性的气候变化使“亚洲水塔”水储量发生失衡[8]，主要表现为藏东南地区冰川快速消融[14-15]、喀喇昆仑地区冰川正平衡(即喀喇昆仑异常)[16-17]、内流区湖泊显著扩张[18-19]等。这些现象表明，“亚洲水塔”水储量及相关水文过程正发生深刻改变，冰湖溃决、冰崩、泥石流等自然灾害频率和强度不断增加，对“亚洲水塔”及泛第三极地区的水资源利用、水灾害防治和可持续发展带来一系列重大挑战[20]。

区域水储量变化的地面观测手段极其有限。近年来在美国对地观测系统(Earth Observing System)、欧洲哥白尼计划(Copernicus Programme)哨兵(Sentinel)卫星系列、中国高分等重大对地观测计划开展下，遥感的时间、空间、光谱和辐射分辨率不断提高，水循环要素的反演算法不断发展，多卫星组网观测和数据融合技术不断成熟，为水储量变化的大范围时空动态监测提供了可能[21-22]。本研究在作者前期工作及数据集的基础上[18，23-26]，集成卫星测高、ICESat-2(Ice，Cloud and land Elevation Satellite)陆冰产品、青藏高原500 m分辨率雪深产品、重力恢复与气候实验卫星(Gravity Recovery and Climate Experiment，GRACE)及其后续卫星GRACE-Follow On(GRACE-FO)等多源遥感信息，对部分水储量变化数据进行延长，并结合再分析数据，解析2000—2020年(以下简称近20 a)青藏高原湖泊水储量、冰川水储量、雪水当量和总水储量的变化及其气候驱动机制。本研究将对认识“亚洲水塔”各类水储量对气候变化的响应和驱动机制提供关键科学支撑，并为“亚洲水塔”失稳的适应性对策制定提供科学依据。

本研究以青藏高原及其典型冰川区藏东南地区(图1)作为研究区域。青藏高原是世界上平均海拔最高的高原，是北半球冰冻圈的重要组成部分。青藏高原位于亚洲中部，地理范围为68°E—104°E，26°N—40°N，总面积约300万km2，其中中国境内面积约250万km2[27]。该地区是中国重要的水能、水资源战略储备区，主要发源于青藏高原的中国西南诸河(包括藏西和藏南诸河、雅鲁藏布江、滇西诸河、怒江、澜沧江和元江)年出境水量约5 800亿m3，约占中国年用水总量6 100亿方(2020年)的95%。高原内降水、气温空间差异显著，藏东南地区年平均气温达20 ℃，年降水量超过1 000 mm；
而西北部地区年平均气温低于0 ℃，年降水量不足100 mm[28]。

图1 青藏高原及藏东南地理位置Fig.1 Location information of the Tibetan Plateau and Southeastern Tibetan Plateau

藏东南地区主要包括帕隆藏布(雅鲁藏布江一级支流)流域及其周边地区，受丰沛的降水和高海拔带来的低温影响，该区域是中国海洋性冰川分布最为集中的地区，冰川水资源储量十分丰富，同时也是青藏高原冰川消融最为剧烈的地区之一。藏东南地区的冰川多发育于山谷，冰舌较长，通常延伸至海拔非常低的区域，主要分布在念青唐古拉山中东段、横断山区的伯舒拉岭、喜马拉雅山东端的南迦巴瓦峰区域。该地区山高谷深，平均海拔超过4 000 m，最高峰为喜马拉雅山东端的南迦巴瓦峰，海拔7 782 m，而最低的雅鲁藏布江大峡谷海拔不到1 000 m。受地形影响，藏东南地区形成了罕见的水汽通道，降水丰沛，是整个青藏高原最为湿润的地区，喜马拉雅山脉及岗日嘎布山以南地区年降水量在1 000 mm以上；
念青唐古拉山南麓的察隅、波密和易贡年降水量也在790 mm以上[29]。

2.1 湖泊水储量

本研究在前期工作基础上，使用了包括ICESat(2003—2009年)、Envisat(2002—2011年)、CryoSat-2(2010—2020年)和Jason-1/2/3(2002—2020年)系列等多源测高数据，补充延长了52个大、中型湖泊测高水位。利用Landsat 5/7/8影像(2002—2017年)提取上述湖泊岸线变化和面积变化，并分别用于生成光学水位和建立湖泊面积—水位关系[18]。不同测高卫星覆盖的湖泊数量不同，ICESat覆盖42个湖泊，Envisat覆盖35个湖泊，CryoSat-2覆盖全部52个湖泊，Jason-1/2/3覆盖12个湖泊，光学水位覆盖全部湖泊。在各种测高卫星中，Envisat和Jason-1/2/3是传统的脉冲有限雷达测高卫星，地面足迹通常在2～4 km，其中Jason-1/2/3系列卫星拥有现有测高数据中最短的重访周期(10 d)；
CryoSat-2是采用了InSAR技术的新一代雷达测高卫星，具有较小的地面足迹(沿轨道方向约300 m，垂直轨道方向约1.65 km)；
ICESat是激光雷达测高卫星，地面足迹约170 m，具有较高的精度。Envisat和CryoSat-2数据均可从欧洲航天局(ESA)官网公开下载(https:∥earth.esa.int/eogateway)，Jason-1/2/3数据可通过AVISO+(Archiving，Validation and Interpretation of Satellite Oceanographic data)官网下载(https:∥www.aviso.altimetry.fr/en/home.html)，ICESat数据可通过美国国家冰雪数据中心(NSIDC)官网获取(https:∥nsidc.org/)。

2.2 冰川表面高程

本研究利用测高卫星ICESat-2计算冰川表面高程变化[23]。ICESat-2(观测时段为2018年10月至今)是ICESat的后续卫星，于2018年9月15日成功发射。其有效载荷是先进地形激光高度计系统(Advanced Topographic Laser Altimeter System，ATLAS)，是目前高程精度最高的星载激光雷达。ICESat-2的轨道高度为496 km，重访周期为91 d。ICESat-2 发射3对(6束)激光束，每对光束包括一强一弱2束光。每个波束对之间的距离为3.3 km，同一波束对中强弱波束之间的距离为跨轨90 m，沿轨2.5 km。同一波束的相邻足迹点之间的距离为 0.7 m。与ICESat相比，ICESat-2在数据覆盖密度、数据精度和空间分辨率方面都有显著提高。ICESat-2 共有 1 387 条轨道覆盖整个地球表面，具有0.1 m的高程精度，可获取精细的冰面地形。本研究使用L3A 陆冰高程数据(ATL06)，空间分辨率约为 20 m，数据可通过美国宇航局(NASA)数据平台EARTHDATA(https:∥search.earthdata.nasa.gov/search)公开获取。

2.3 积雪面积和雪水当量

本研究利用Dai等[24]开发的青藏高原500 m分辨率雪深产品对青藏高原积雪时空变化规律进行研究。青藏高原500 m雪深产品综合利用了MODIS积雪面积、地表温度和被动微波亮度温度来反演雪深，利用地表温度与被动微波K和Ka波段水平极化亮度温度的关系，消除土壤对亮度温度差(TBD)的贡献，然后将被动微波格网上的TBD分配到有雪像元，最后利用TBD和雪深的关系提取雪深，得到的雪深产品分辨率和准确度都较高。此外，本研究根据由欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts，ECMWF)官网(https:∥www.ecmwf.int/en/era5-land)公开获取的ERA5-Land再分析资料，利用其月均积雪变量相关数据(雪深、雪密度、雪水当量)和上述500 m分辨率雪深数据进行对比，分析青藏高原积雪变化规律。ERA5-Land数据的空间分辨率为0.1°×0.1°，时间分辨率为1 h。

2.4 总水储量

重力卫星GRACE(观测时段为2002年4月至2017年6月)及其后续卫星GRACE-FO(观测时段为2018年5月至今)的研制和发射，为陆地总水储量变化(以下简称总水储量变化，包括地表水储量(如冰川、湖泊、积雪等)、土壤水储量和地下水储量的总变化)的反演开辟了新途径。重力卫星由2颗完全一样的卫星构成，轨道距离地球约450 km，2颗卫星的距离约220 km，将双星之间的距离变化作为反演地球时变重力场的基本信息，可进一步获得由地球质量重分布引起的总水储量变化。重力卫星可提供月尺度，空间分辨率约10万km2(约300 km×300 km)的总水储量变化信息[30]。本研究对比4种不同的GRACE三级产品在青藏高原的表现，包括2种基于球谐系数的网格产品(JPL-SH和CSR-SH)和2种基于质量浓度斑块(mascon)的网格产品(JPL-M和CSR-M)。其中，2种球谐系数产品为RL05版本，其名义分辨率为1°×1°；
2种mascon产品为RL06版本，其名义分辨率分别为0.5°×0.5°(JPL-M)和0.25°×0.5°(CSR-M)。重建GRACE和GRACE-FO之间的缺失数据时，采用最新的GRACE JPL RL06版本。所有重力卫星数据均为总水储量相对于2004—2009年时段的距平值，可由NASA官网公开获取(http:∥grace.jpl.nasa.gov)。

3.1 湖泊水储量变化

采用湖泊面积变化和湖泊水位变化观测计算湖泊水储量变化。湖泊面积的获取主要通过遥感影像识别湖泊边界，湖泊水位的获取可通过测高卫星直接观测。此外，可以通过对湖泊岸坡局部地形的分析，建立特定区域岸线位置与水位之间的关系来推求湖泊水位。本研究采用了基于卫星测高(测高水位)和基于光学影像的岸线变化(光学水位)2种方法来推求湖泊水位[18]。

基于云平台GEE(Google Earth Engine)和水体指数法提取湖泊面积和岸线变化包括4个基本步骤：① 构建包含湖泊边界/岸线的缓冲区；
② 构建水体指数；
③ 在缓冲区内进行水体指数计算，并选取合适的阈值将水体指数灰度图二值化；
④ 根据二值化图像计算水体面积后导出。常用的水体指数有归一化差异水体指数(NDWI)和改进的归一化差异水体指数(MNDWI)。其中，MNDWI可以更有效识别浑浊水体，更适用于夏季；
NDWI可以一定程度区分湖岸的积雪和湖水，更适用于冬季。自动阈值的选取通常基于最大类间方差法，即大津法(Otsu′s Method)。

卫星测高最初应用于海平面和冰盖监测，近年来在内陆水体以及河、湖、冰研究中得到广泛应用[18，31-33]。雷达测高卫星数据通常还需要进行波形重定，目的是消除复杂地形引起的雷达回波信号异常，从而提高水位的反演精度。Huang等[31]提出的50% Threshold and Ice-1 Combined(TIC)波形重定算法，在改进阈值法的基础上，优化了测高足迹的选取，并针对不同季节使用不同的阈值，该方法可以显著改善河湖水位的反演精度。在湖泊水位提取中，由于整体地形平坦、水位变动幅度较小，仅使用TIC算法的50%阈值模式即可满足应用需求。另一方面，不同的测高卫星轨道密度、重访周期以及精度各有不同，使用多源测高数据能够延长和加密湖泊测高水位序列，但往往需要消除不同传感器之间的系统误差。在获取多源雷达测高水位和光学水位后，可以通过比较数据重叠期的平均水位值，移除不同源测高水位之间的系统偏差[18]。经与实测水位检验以及理论推导，修正系统误差后的光学和测高水位精度(均方根误差)均可达到0.2 m以内[18，34]。

通过建立湖泊岸线变化与同期测高水位之间的统计回归关系，可将岸线变化转化为水位信息，即光学水位。由于光学水位覆盖的时段最长，以光学水位为参照，通过比较数据重叠期的平均水位值，可将多源测高卫星数据之间的系统误差全部移除，形成时空一致性良好的高时间分辨率水位序列。获取湖泊水位和面积变化时间序列后，可以通过积分得到湖泊水量变化(常见的方法有台体公式，面积—水位关系积分得到体积—水位关系等)。

3.2 冰川质量平衡

本研究利用ICESat-2 L3A陆冰高程(ATL06 V3)产品，计算藏东南地区的冰川表面高程变化[23]。使用的参考DEM是最新发布的NASADEM(空间分辨率约30 m)，通过计算ICESat-2与NASADEM之间的高程差，即可得到2000—2019年间冰川表面高程变化值，进而求得冰川厚度的变化速率和冰川物质平衡。经与野外考察获取的GPS测点高程以及无人机获取的帕隆四号冰川超高分辨率DSM(数字表面模型)相比，ICESat-2数据在相同位置的高程误差在3 m以内，具有极高的精度。

Δh=hICESat-2-hNASADEM+hp

(1)

式中：Δh为冰川表面高程变化值；
hICESat-2和hNASADEM分别对应ICESat-2足迹点的高程值和对应NASADEM的高程值；
hp为NASADEM对冰川的穿透深度， 在本研究中hp=1.5 m。

(2)

(3)

式中：R为冰川厚度变化速率；
M为冰川物质平衡；
ρg和ρw分别为冰川和水的密度，ρg=850 kg/m3。

3.3 总水储量变化

3.3.1 时间序列分解方法

本研究采用STL(Seasonal and Trend decomposition using Loess)时间序列分解方法，获取总水储量变化的长期趋势(式(4) )，以消除季节性变化和残差对评估总水储量变化的影响。

Stot=Slt+Ssea+R

(4)

式(4) 将总水储量距平(Stot)分解为长期趋势项(Slt)、 季节性变化(Ssea)和残差项(R)。

对分解后的长期趋势项序列作线性拟合， 可获取总水储量变化的趋势。

STL以鲁棒局部加权回归作为平滑方法， 被广泛应用于总水储量时间序列分析。

3.3.2 重力卫星缺失数据重建

GRACE和GRACE-FO之间的数据缺失(2017年7月至2018年5月)，造成总水储量变化监测连续性不足。近年来，机器学习模型在水文研究中的应用日趋成熟，通过学习输入层及输出层变量之间的非线性相关关系，机器学习算法可以基于相关输入变量对重力卫星缺失数据进行重建[35]。应用于重建重力卫星数据的机器学习算法经历了由相对简单的自回归移动平均模型(ARIMA)、多元线性及非线性回归模型(multiple linear and nonlinear regressions)向较复杂的人工神经网络(ANN)、深度神经网络(DNN)及卷积神经网络(CNN)发展。相关研究中，Mo等[36]在全球尺度利用贝叶斯卷积神经网络对GRACE缺失数据进行重建，该研究参照Long等[37]的方法将总水储量变化信号分解为趋势项和非趋势项，可改进前人研究中气候因子较难捕捉总水储量长期变化趋势的这一缺陷，提升重建结果的可信度。本研究补充Mo等[36]的重建结果，用以分析2002—2020年完整的青藏高原水储量变化，具体实施步骤如下：

(1) 确定输入层变量。所有变量均由再分析资料ERA5-Land数据获取，包括降水、气温、模拟的总水储量距平(土壤水、雪水及冠层水储量的加和)，以及累积总水储量变化量(一段时间内总降水减去总蒸散及总径流)；

(2) 对输入层所有变量去除线性趋势，仅保留年际和季节性变化；

(3) 确定输出层变量。针对不同的GRACE三级产品，将对应的GRACE及GRACE-FO时间序列去除线性趋势，保留年际和季节性变化信号作为神经网络的目标变量；

(4) 基于贝叶斯卷积神经网络，依据输入的气候变量重建缺失时段总水储量年际和季节性变化信号；

(5) 缺失时段最终总水储量变化信号为第(3)步确定的线性趋势与第(4)步确定的年际和季节性变化信号的和。

4.1 2000—2020年青藏高原各类水储量变化

4.1.1 湖泊水储量显著扩张

本研究研制了2000—2020年青藏高原52个大、中型湖泊的高时间分辨率水位、水量数据集。52个湖泊水储量总体的不确定性约为±6%[18]，整体呈扩张趋势，2000—2020年间湖泊总水储量增加约130 Gt，超过3个三峡水库的总库容。湖泊的水量变化可划分为3个阶段：① 2000—2012年水量增加速率较快，为6.35 Gt/a，其中内流区湖泊水量增长速度约为5.95 Gt/a；
② 2012—2017年整体增速放缓至1.42 Gt/a，其中内流区湖泊为1.20 Gt/a；
③ 2017年后湖泊再次出现快速扩张，整体增速达到10.59 Gt/a，其中内流区湖泊达到8.18 Gt/a(图 2(a))。从内流区的年降水量时间序列上可以看出，年降水量显著高于多年平均降水量时，内流区湖泊整体水量会快速增加，二者具有较高的一致性(图 2(b))。

图2 青藏高原52个大、中型湖泊累积水量与青藏高原内流区年降水量变化时间序列Fig.2 Accumulated lake water storage change for 52 large and middle lakes and annual precipitation in the endorheic basin on the Tibetan Plateau

从湖泊水量变化趋势的空间分布看(图 3)，内流区湖泊扩张趋势比外流区更为显著，但在部分流域呈现出“大湖扩张，小湖收缩”的情形[18]。以色林措流域为例，色林措水量快速扩张(0.93 Gt/a)，而其上游的吴如措(-0.007 Gt/a)、格仁措(-0.007 Gt/a)和错愕湖(-0.004 Gt/a)水量则缓慢下降，这些湖泊与色林措之间有河道相连。可能原因是随着水流对上游湖泊下泄口的侵蚀，下泄口高程缓慢降低，湖泊蓄水能力下降，而时间分辨率较低的数据集难以反映这种现象。在外流区，黄河上游流域湖泊水量增加较为显著，例如青海湖(0.54 Gt/a)、鄂陵湖(0.08 Gt/a)、扎陵湖(0.03 Gt/a)等。雅鲁藏布江流域的主要湖泊出现水量下降情况，如羊卓雍措(-0.14 Gt/a)。

图3 青藏高原52个大、中型湖泊2000—2020年水量变化趋势Fig.3 Spatial pattern of water storage trends of 52 large and middle lakes for 2000—2020 in the Tibetan Platean

湖泊持续扩张可能引发湖泊溃决或溢流洪水灾害，如青海可可西里地区的卓乃湖-库赛湖-海丁诺尔湖-盐湖系统。位于上游的卓乃湖于2011年发生溃决，湖水下泄到库赛湖并引起库赛湖溢流，进而注满了海丁诺尔湖，最后注入盐湖。据估算，溃决后从卓乃湖下泄进入库赛湖的水量达到2.47 Gt，造成库赛湖水位上升近8 m并引发库赛湖溢流。盐湖的面积近年来扩大了数倍，也存在溢流的风险，其下泄点距离青藏铁路沿线仅10 km左右，对铁路运输安全可能构成威胁[18]。

4.1.2 藏东南地区冰川后退

通过2020年的ICESat-2数据与2000年的NASADEM，可得到2000—2020年间20 a整个藏东南地区的平均冰川表面高程变化速率，该数据集的总体误差约为±25%。20 a间该地区平均冰川表面高程变化速率为(-0.71±0.18)m/a，相当于每年损失44.7亿t的冰川水资源。分析2018年10月至今近3 a的ICESat-2数据，得到近3 a藏东南地区季节平均的冰川表面高程变化时间序列如图4所示，其中逐月季平均指取该月前后3个月数据作平均；
季平均指取每个季度的数据作平均；
滑动平均是对季平均时间序列作周期为4的滑动平均；
线性回归是对季平均时间序列的线性拟合。可知近3 a藏东南地区冰川表面高程变化速率为(-0.97±0.59)m/a。Zhao等[23]得出藏东南地区10 a(2011—2020年)的冰川表面高程变化速率为(-0.87±0.13)m/a，故藏东南地区的冰川消融在不断加速。通过高空间覆盖范围的ICESat-2数据，分析藏东南地区冰川物质平衡的空间分布特征(图5)。在0.5°×0.5°的网格尺度，藏东南地区冰川物质损失最严重的区域位于横断山区域(藏东南地区的东北部)，2000—2020年间物质平衡速率达到了(-1.43±0.36)m w.e./a(w.e.表示水当量)。藏东南地区东部的冰川物质损失比西部地区要快很多，东部区域的表面高程变化速率约为-1.0 m/a，而西部区域的表面高程变化速率在-0.5 m/a左右。此外，藏东南地区成片分布的大型冰川的表面高程变化速率，要低于周边零星分布的小型冰川的速率[23]。

图 4 近3 a藏东南地区冰川表面高程季节变化Fig.4 Seasonal time series of cumulative glacier elevation changes in the Southeastern Tibetan Plateau in recent three years

图5 2000—2020年藏东南地区冰川表面高程变化空间分布(0.5°×0.5°网格)Fig.5 Spatial distribution of glacier elevation change during 2000—2020 in the Southeastern Tibetan Plateau

4.1.3 积雪深度和雪水当量

2001—2015年青藏高原逐日500 m分辨率雪深产品季节性分布显示(图 6)，总体上积雪主要分布在山区，冬春季雪深较深，其次是秋季，夏季积雪零星分布在山区。雪深年内变化一般从10月开始增加，2月达到峰值，而后逐渐消融减少，至6月完全融化。对比500 m分辨率雪深产品和ERA5-Land再分析产品(图6)可知，2种产品的雪深空间分布一致性较高，然而ERA5-Land产品雪深数值明显高于青藏高原500 m分辨率雪深产品。这与再分析数据存在对雪深或积雪面积的系统性高估有关[38]，而过大的降雪量可能是ERA5-Land再分析产品中雪深高估的主要因素。

图6 青藏高原500 m分辨率雪深产品2001—2015年四季平均雪深空间分布Fig.6 Spatial distributions of multi-year (2001—2015) averaged snow depth derived from a 500 m snow depth product for four seasons

图7 ERA5-Land产品2001—2015年四季平均雪深空间分布Fig.7 Spatial distributions of multi-year (2001—2015) averaged snow depth derived from the ERA5-Land product for four seasons

根据ERA5-Land产品的积雪密度数值，进而通过雪深、雪密度与雪水当量之间的关系，计算得到青藏高原500 m分辨率雪深产品的雪水当量时间序列，进一步对比其和ERA5-Land生成的雪水当量时间序列可知(图 8)，ERA5-Land产品数值明显高于500 m雪深产品，但二者相关性很高。其中，雪深时间序列相关系数达0.83，雪水当量时间序列相关系数达0.86。因此，ERA5-Land产品可一定程度反映青藏高原积雪时空变化规律，但其数值存在显著高估。

图8 2001—2015年青藏高原500 m分辨率雪深产品和ERA5-Land产品平均雪深和雪水当量Fig.8 Time series of regionally averaged snow depth and snow water equivalent over the Tibetan Plateau during 2001—2015 for a 500 m snow depth product and the ERA5 product

对比2001—2015和2016—2020年间ERA5-Land产品反映的青藏高原多年月平均雪水当量(图9)可以发现，相比2001—2015年，2016—2020年间除7—10月积雪较少的月份，其余月份多年月平均雪水当量呈增加趋势。从年平均雪水当量看，2001—2015年间青藏高原多年平均雪水当量为13.6 cm，而2016—2020年间多年平均雪水当量为14.1 cm，增加了3.7%。

图9 2001—2015年和2016—2020年间ERA5-Land产品反映的青藏高原多年月平均雪水当量Fig.9 Multi-year monthly averaged snow water equivalent derived from the ERA5-Land product for 2001—2015 and 2016—2020 over the Tibetan Plateau

4.1.4 总水储量

本研究在获取2002—2020年完整的总水储量变化时间序列后，利用STL时间序列分解获得水储量多年变化趋势(图 10)。4种GRACE反演结果均显示2002—2020年青藏高原外流区(印度河、恒河、雅鲁藏布江、怒江、澜沧江等)总水储量呈显著下降趋势；
内流区(塔里木、柴达木和羌塘盆地)的总水储量呈显著增加趋势。但GRACE不同产品在估算总水储量变化趋势的量级方面差异较大，2种球谐系数产品的结果比较一致，2种mascon产品的结果在外流区比较一致，但JPL-M比CSR-M能反映更丰富的水储量空间变化信息。

图10 2002—2020年GRACE不同三级产品反演的青藏高原水储量变化趋势Fig.10 Trends in total water storage over the Tibetan Plateau during 2002—2020 estimated from different GRACE solutions

基于GRACE JPL-M估算的总水储量变化趋势(图 11)，2002—2020年间青藏高原总水储量以-4.50 Gt/a的速率下降。其中，外流区(印度河-恒河-雅鲁藏布江-怒江-澜沧江-长江-黄河流域)总水储量的下降速率达到-10.90 Gt/a，内流区(阿姆河-塔里木-羌塘盆地-柴达木-河西区域)总水储量则以6.40 Gt/a的速率上升。在冰川广泛分布的印度河-恒河-雅鲁藏布江流域(青藏高原海洋性冰川集中分布的地区)，总水储量以-9.50 Gt/a的速度下降，冰川快速剧烈消融是该区总水储量下降的主要原因。例如，藏东南地区(念青唐古拉山中东段,横断山区的伯舒拉岭,喜马拉雅山东端的南迦巴瓦峰区域)的冰川消融速率达到-4.50 Gt/a，占该地区总水储量下降速率(-5.90 Gt/a)的76%。印度河-恒河-雅鲁藏布江流域的冰川是极为重要的固态水资源，冰川与积雪融水是下游印度河-恒河平原的主要灌溉水源，影响约1.3亿人的生产生活用水[39]，冰川快速消融已成为该地区最严重的威胁之一。对于湖泊广泛分布的羌塘盆地，总水储量增加速率为5.20 Gt/a，主要由湖泊扩张导致(湖泊水储量增加速率为4.90 Gt/a，占总水储量变化的94%)。尽管内流区湖泊水储量与下游地区的用水量联系较弱，但显著的湖泊扩张可能导致湖泊溃决，进而威胁周边基础设施和居民生命财产安全，如2011年9月青藏高原内流区卓乃湖和库塞湖的溃决事件[40-41]。

图11 基于GRACE JPL-M获取的2002—2020年总水储量变化时间序列Fig.11 Time series of total water storage changes over the Tibetan Plateau，Indus-Ganges-Brahmaputra，and Inner Tibetan Plateau basins during 2002—2020 estimated from GRACE JPL-M

4.2 水储量变化驱动机制

4.2.1 湖泊

降水是驱动湖泊水量变化最主要的因素。在定量解析气候变量对湖泊水量平衡贡献的研究中(包括色林措、纳木措、当若雍措、玛旁雍措、佩古措、可可西里湖、勒斜武担措等湖泊)[42-43]，降水对不同湖泊贡献率为59%～95%。青藏高原内流区是湖泊最集中的区域，近30 a降水量呈增加趋势[44]。色林措、纳木措和当若雍措均位于青藏高原内流区南部，这一区域的湖泊在2000—2010年前后经历了快速扩张期，之后进入缓慢扩张/相对稳定期。可可西里湖和勒斜武担措位于内流区的北部，呈现持续的水量上涨趋势。整体上这些湖泊的水量变化与降水量的变化保持了较高的一致性。

冰川融水的贡献可以通过模型模拟或遥感观测获得。冰川融水对湖泊水量的贡献大多为10%～20%，对佩古措水量的贡献较大，达到了30%[4]。从冰川分布看，内流区北部的冰川面积大于南部地区，尤其喀喇昆仑地区冰川较为集中，但由于该区冰川质量整体保持稳定，冰川融水对其流域内湖泊扩张的贡献相对有限。由于观测资料匮乏，分析冻土对湖泊水量的贡献通常基于模型模拟。针对色林措流域的模拟结果显示，冻土对色林措水量变化的贡献达到了28%[45]。

蒸散也是影响湖泊水量平衡的重要变量，是水离开内流流域的主要途径。在暖湿化背景下，青藏高原蒸散整体呈增加趋势。温度升高、冻结期缩短都将导致湖泊总蒸发量增大。同时，青藏高原湖泊面积扩张也可能导致水面蒸发总量增加，但青藏高原的湖泊水量不可能无限增长，将达到某种平衡状态[44]。这种平衡状态所对应的湖泊面积、湖泊水量、气候条件以及对气候变化的响应规律有待进一步探索。

4.2.2 冰川

气温和降水是影响冰川物质平衡的重要因素。降水的多少决定了冰川物质积累量的大小，气温的高低决定了冰川消融(积累)的快慢，降水和气温的变化是冰川变化的重要驱动因素[46]。基于观测数据的统计结果表明，1993—2012年全球平均升温速率约为0.14 ℃/(10 a)[47]；
藏东南地区的气象站点数据资料显示，2003—2018年间，藏东南地区年均气温以0.23 ℃/(10 a)的速率上升，显著高于全球平均水平。年降水量则基本保持稳定，但降水的年际变化幅度明显高于以往。气温升高使冰川消融速度加快，同时，降水相态也随之发生改变[48-50]。当越来越多的降雪转变为降雨时，冰川的物质积累量减少，且降雨对冰川的冲刷作用增强。近些年来，气温高和降水少的年份多次出现(2009年、2014年和2018年)，这些年份的冰川消融尤为剧烈。气温升高和降水波动是近年来藏东南地区冰川后退的主要原因。

4.2.3 总水储量

在人类活动影响相对较小的青藏高原，气候变化是引起该地区近20 a来总水储量显著变化的主要原因。总体而言，降水和气温是影响总水储量变化最重要的气候因子。降水作为总水储量的输入，是引起冰川物质正平衡、湖泊扩张、土壤和地下水储量增加的重要变量；
气温与冰川融化及蒸散发密切相关，气温上升是引起总水储量下降的主要原因。研究表明近年来南亚季风逐渐减弱[51]，导致青藏高原南部(喜马拉雅地区)由南亚季风带来的降水减少，进而引起该区总水储量下降。此外，近年来西风带增强[52]导致帕米尔高原降水有所增加，是青藏高原西部阿姆河流域冰川物质基本保持平衡的主要原因之一[2]。季风变化深刻影响了青藏高原降水变化的时空格局，进一步影响青藏高原冰川等固态水储量和总水储量的变化。

基于观测数据统计结果，青藏高原1955—1996年升温速率约为0.16 ℃/(10 a)，其中冬季升温速率高达0.23 ℃/(10 a)[13]。显著的升温导致青藏高原冰川快速消融，进一步引起总水储量变化。例如，在海洋性冰川集中分布的藏东南地区，气温升高是冰川融化和总水储量减少的主要气候原因。

本研究主要基于多源遥感反演(包括卫星重力、卫星测高和光学遥感影像等)对2000—2020年青藏高原各类水储量变化规律及其气候驱动机制进行了研究，主要结论如下：

(1) 气候暖湿化下降水量和冰川融水增加，引发了青藏高原湖泊水量整体上涨。52个大、中型湖泊水量变化过程分为3个阶段：2000—2012年为平稳增长阶段(6.35 Gt/a)，2012—2017年为相对稳定阶段(1.42 Gt/a)，2017年后进入快速增长阶段(10.59 Gt/a)。

(2) 受气温上升及其带来的降水相态改变等多种因素影响，藏东南地区的冰川在快速消融，达到-4.50 Gt/a，且有加快趋势。

(3) 2001—2015年间青藏高原多年平均雪水当量为13.6 cm，而2016—2020年间多年平均雪水当量为14.1 cm，增加了3.7%，积雪变化主要受累积期平均气温和降水影响。

(4) 在气温上升、南亚季风减弱和西风带增强综合影响下，青藏高原外流区总水储量下降速率达到-10.90 Gt/a，内流区总水储量则以6.40 Gt/a的速率上升。

在2000—2020年间，气候变化显著改变了青藏高原水储量的状态，由此带来水资源量改变及其上下游分配等问题将产生安全、生态、经济、政治连锁效应，制定全球和区域气候变化的减缓和适应性对策并付诸实践迫在眉睫。

猜你喜欢 储量青藏高原湖泊 青藏高原上的“含羞花”今日农业(2022年3期)2022-11-1613.22亿吨资源导刊(2021年10期)2021-11-05耳钉小天使·五年级语数英综合(2020年8期)2020-12-23摸清甘肃迭部县不同生态系统碳储量家底：丰富且价值重大各级政府部门应加快建立与完善碳汇补贴机制时代人物(2020年24期)2020-11-28“害羞”的湖泊阅读(低年级)(2018年4期)2018-05-14为什么南极会有温水湖泊？文理导航·科普童话(2017年4期)2018-02-10奇异的湖泊小学阅读指南·低年级版(2017年2期)2017-03-232013年我国主要矿产资源勘查情况能源研究与信息(2014年1期)2014-04-29青藏高原筑“天路”小哥白尼·趣味科学画报(2006年1期)2006-02-15