卫星重力与光学遥感组合的念青唐古拉山脉冰川变化分析

周 苗，常晓涛，朱广彬，瞿庆亮，刘 伟

1. 自然资源部国土卫星遥感应用中心，北京 100048； 2. 山东科技大学测绘与空间信息学院, 青岛 266590

冰川是冰冻圈的主要组成部分，对气候变化极为敏感，被认为是气候变化的最佳天然指示器之一[1]。冰川主要分布在极区和高海拔地区，随着气候的周期性变化处于一种动态平衡当中。随着全球变暖，大多数冰川呈加速消融趋势[2-4]，增加了发生冰湖溃决、泥石流等冰川灾害的可能[5]。念青唐古拉山脉是我国第二大冰川聚集区，冰川消融量巨大，2000—2015年间念青唐古拉山脉冰川总面积减少了463.36 km2,消融量高居亚洲其他山脉之首。开展念青唐古拉山脉冰川变化监测对区域生态保护、灾害预警等具有重要的现实意义。

随着空间对地观测技术的发展，冰川变化监测技术更加多样化，监测效率更加高效。近年来，卫星重力和光学遥感技术在冰川研究中广泛应用。卫星重力可直接获取区域质量变化，在研究冰川变化机制和大、中尺度质量变化方面具有一定优势。文献[6]利用GRACE、GRACE-FO时变重力场模型解算的等效水高序列评估了2002年4月至2019年9月期间世界(不包括南极和格林兰岛)冰川和冰盖的质量平衡。结果表明，世界冰川和冰盖的平均损失速率为281.5±30 Gt/a。文献[7]基于GRACE卫星重力资料指出，2003—2015年间青藏高原冰川质量一直减少。文献[8]利用GRACE卫星重力数据在频域上计算了2002—2014年喜马拉雅山冰川质量变化，发现该区域冰川质量变化整体上呈现加速消融趋势。光学卫星遥感技术能很好地反映冰川变化的细节，是提取冰川面积和冰面储量变化的最快速、最有效手段。文献[9]利用Landsat TM/ETM+/OLI影像数据，结合比值阈值法与目视解译提取冰川边界，对1990—2015年喜马拉雅山冰川的消融情况进行了研究。分析认为，近25年来喜马拉雅山冰川整体呈退缩趋势，冰川面积由23 229.27 km2减少到20 676.17 km2，退缩率为10.99%，并有加速退缩的趋势。文献[10]利用资源三号立体影像和SRTM DEM数据对念青唐古拉山脉2000—2017年的冰川物质平衡变化进行估计。研究发现，念青唐古拉山脉整体冰川物质平衡为负值，其中念青唐古拉山脉西段2013—2017年间冰川储量损失最快。文献[11]基于地形图和Landsat TM/OLI遥感影像，利用目视解译和波段比值法研究了南迦巴瓦峰地区1980—2015年的冰川变化。研究发现南迦巴瓦峰地区冰川面积持续减小并呈加速退缩的趋势。

国内外学者在利用卫星重力和遥感技术实现冰川质量、面积变化方面进行了较为广泛的研究，为冰川质量变化监测做出了重要贡献，但两种技术协同应用方面仍有一些不足，具有一定的研究空间。本文结合GRACE、GRACE-FO卫星重力数据和遥感资料对念青唐古拉山脉冰川变化展开研究，并利用气象数据对冰川的变化机制及影响因素进行分析，为该区域冰川水资源变化监测及生态保护提供参考。

1 研究方法

1.1 地表质量变化反演

利用重力场模型的球谐系数变化计算地球表面质量变化的模型为

(1)

本文采用美国得克萨斯大学空间研究中心提供的GRACE、GRACE-FO RL06版本时变重力场数据，模型的最大阶数为96，时间跨度为2002年8月至2020年4月，共有33个月的数据缺失，其中2017年6月至2018年6月这11月的数据空白是由于卫星任务更替造成的，对除了这11个月的其他缺失数据采用3次样条插值方法补齐。由于GRACE、GRACE-FO卫星对低阶项不敏感，本文使用文献[12]提供的一阶项和卫星激光测距(satellite laser ranging，SLR)获得的C20项[13]代替球谐系数中的相应项。另外，从2016年10月开始，GRACE和GRACE-FO卫星上仅有一颗加速度计数据用于地球时变重力场的计算，这会造成该段时间GRACE、GRACE-FO得到的C30项精度较低，本文使用SLR获得的C30对其进行替换[14]。同时，为了削弱高阶次噪声和奇偶阶相关误差的影响[15]，对位系数进行300 km高斯滤波和滑动窗去相关滤波[16]。最后，结合式(1)得到更精确的地表质量变化反演公式

(ΔClmcos(mλ)+ΔSlmsin(mλ))

(2)

式中，w(l)为高斯平滑函数。

经过以上处理得到的结果中仍包含冰后回弹(glacial isostatic adjustment,GIA)和泄露误差的影响，本文利用ICE-6G-C模型[17]进行GIA改正。同时，采用文献[18]提出的空域法对泄露信号进行恢复。空域滤波法的原理是假定地表质量变化是由一系列位置已知、质量未知的质量块(Mascon)引起的，根据Mascon的球谐展开和综合结果与GRACE结果的线性关系，可推求每个Mascon对应的系数，此系数序列即为各个Mascon对应的区域平均结果。

1.2 冰川质量变化计算

念青唐古拉山脉为冰川聚集区，该地区质量变化主要由冰川消融主导。因此，根据水量平衡原理，冰川质量变化可表示为

ΔICE=ΔTWS-ΔSM-ΔSWE

(3)

式中，ΔICE为冰川储量变化；ΔTWS为陆地水储量变化；ΔSM为土壤水储量变化；ΔSWE为雪水当量变化。其中，ΔTWS为GRACE、GRACE-FO计算结果，主要包括地下水、土壤水、积雪水、地表径流、植被水等。ΔSM和ΔSWE均可在GLDAS水文模型中得到。

GLDAS(global land data assimilation system)水文模型由美国宇航局戈达德太空飞行中心和美国国家环境预报中心联合建立，包括NOAH、VIC、CLM、MOSAIC 4种陆地表面模式。本文采用2003年1月至2019年12月Noah陆面模式的Version1数据[19]，时间分辨率为1个月，空间分辨率为1°×1°。选取该数据的4层土壤水(0～200 mm)、雪水当量和冠层水计算区域地表水储量变化。

1.3 气象数据处理

降水数据采用美国国家宇航局和日本空间发展局提供的TRMM(tropical rainfall measuring mission)全球降雨数据产品[20]，空间分辨率为0.25°×0.25°，时间分辨率为1个月。该产品提供的是月降雨率(mm/h)数据，为了反映降水和水储量变化的关系，将其转换成对应月份的降雨量。为了保持与GRACE数据空间分辨率的一致，采用三次样条插值方法将降雨数据插值成1°×1°，然后扣除各个格网点上的平均值，得到中国区域1°×1°的月降雨异常。

温度数据采用美国国家海洋和大气管理局提供的戈达德太空研究所对全球表面温度分析数据集[21](GISS Surface Temperature Analysis，GISTEMP)，数据的空间分辨率为2°×2°，时间分辨率为1月，覆盖范围是89.0N—80.0S,1.0E—359.5E，该数据集提供了1980年至今的所有数据，具有较好的时间连续性。同降雨数据的处理步骤类似，得到中国区域的月温度异常。

1.4 冰川信息提取

1.4.1 冰川面积提取

冰川面积提取主要借助Landsat光学遥感影像，利用冰雪在可见光波段的高反射特性和在短波红外的强吸收特性，通过波段比值法自动提取冰川边界。波段比值法的关键在于阈值选取，经过大量试验得到阈值为2.3时可以更加有效区分冰川和陆地，但对水体和阴影错分较多。为了获得更加准确的结果，采用目视解译法对错分的水体和阴影进行修正，得到较为精确的冰川边界。

冰川研究中遥感影像筛选应满足少云、无积雪覆盖、消融末期3项原则，基于此原则经过大量筛选，最终选取了2003-12-20(TM)、2017-10-07(OLI)两景Landsat影像进行预处理，主要包括辐射定标、大气校正和影像裁剪等。

1.4.2 冰川储量估计

本文使用国家青藏高原科学数据中心(http:∥data.tpdc.ac.cn/zh-hans/)提供的第三极地区冰川表面高程变化数据产品[22]。该数据基于2000年的SRTM(shuttle radar topography mission)和2015年前后ASTER立体相对数据制成，估计了第三极地区40多个典型冰川表面高程变化，空间分辨率为30 m。通过获取念青唐古拉山脉西段冰川范围内每个像元的高程变化值，结合式(4)可以得到2000—2015年间该地区的冰川体积变化ΔV

(4)

式中，Δhi表示第i个像元的高差；sp为单个像元面积；n为像元的个数。

2 结果与分析

2.1 冰川质量变化

图1为GRACE和GLDAS反演的念青唐古拉山脉陆地水储量和地表水储量变化序列。从图1可看出，陆地水储量变化和地表水储量变化均表现出明显的周期性，分别在冬季和夏季达到最小值和最大值。在时间尺度上，GRACE的结果和GLDAS的结果具有较好的一致性，但GRACE结果的振幅较大，这是因为GRACE的结果是综合各种因素后的水储量变化，而GLDAS水文模型只包含表层土壤水变化，不包括深层土壤水、冰川、地下水、地表河流及湖泊等水储量变化[23]。对GRACE和GLDAS的结果进行最小二乘线性拟合，得到他们的变化速率分别为-2.76 cm/a和-0.48 cm/a，这说明地表水储量变化对念青唐古拉山脉整体质量亏损影响较小。

图1 念青唐古拉山脉陆地水储量变化和地表水储量变化曲线Fig.1 Change curve of land water reserve and surface water reserve in Nyainqentanglha Mountain

念青唐古拉山脉不存在大型河流和湖泊，地下水和深层土壤水在冻土层的保护和补给也相对稳定，因此，计算的冰川质量变化结果受其他因素的影响较小。图2给出了念青唐古拉山脉冰川质量变化序列，其中蓝线表示原始冰川质量变化、红线表示经过泄露误差改正的结果。图2中蓝线的变化速率为-2.48 cm/a，对总质量变化的贡献达到了89%，表明念青唐古拉山脉2003—2019年冰川质量亏损严重，尤其是2013年之后念青唐古拉山脉冰川有加速消融趋势，这一结果与文献[10]的结论一致。此外，在2015年1—3月间念青唐古拉山脉地区存在一个明显的质量变化异常信号，这一异常在地表水储量变化序列中未出现，结合2015年3月前后研究区附近地震的分布情况，认为该质量异常可能与地球内部的物质迁移有关。

图2 冰川质量变化和空域法改正曲线Fig.2 The curve of glacier mass change and spatial domain method correction

对比图2两条曲线能够看出(图中坐标轴不一样)，两条曲线变化趋势基本一致但振幅相差较大，这与空域法建立的Mascon起到的质量约束作用有关。图2红线为泄露改正后的念青唐古拉山脉冰川储量变化曲线，其变化速率为-15.53 cm/a，同时计算了念青唐古拉山脉西段和东段两个Mascon的质量变化速率分别为-20.23 cm/a和-10.82 cm/a，其中念青唐古拉山脉东侧变化速率结果与文献[24]使用空域法的结果基本一致。

为探究冰川质量变化及影响机制，本文提取了夏季温度异常和降水异常的变化曲线(图3、图4)。由图3可以看出，2003—2019年间夏季温度呈缓慢上升趋势，研究时段内平均温度大约升高了1℃，而冰川质量变化明显减少。从整体看，两条曲线具有明显的相关性，其相关系数为-0.63。同时，计算了春季、秋季和冬季温度异常与冰川质量变化曲线的相关系数，相关性均不明显，其中秋季的相关系数较大也仅为-0.36。由图4可以看出，念青唐古拉山脉降水有缓慢减少的趋势，表明温度升高和降水减少是致念青唐古拉山脉冰川质量亏损的主要原因。

图3 念青唐古拉山脉夏季冰川质量变化和温度异常曲线Fig.3 The glacier quality change and temperature abnormal curves in the Nyainqentanglha Mountains in summer

图4 念青唐古拉山脉夏季降水异常Fig.4 Abnormal precipitation in the Nyainqentanglha Mountain in summer

2.2 冰川面积变化

念青唐古拉山脉东段为海洋型冰川，冰川活动剧烈，常年云量较多，很难获取同一尺度、覆盖整个范围的光学影像。因此，在冰川面积及储量变化研究中以念青唐古拉山主峰所在区域的冰川为研究对象，该冰川位于念青唐古拉山脉西段(89.89°E—90.84°E，29.83°N—30.57°N)，分布独立，靠近纳木错湖泊，平均海拔约5500 m，是较典型的亚大陆型冰川[25]。

利用两幅Landsat影像获取了念青唐古拉山脉西段冰川边界，经统计2003年冰川总面积为619.305 km2，2017年时为547.003 km2，这14年间冰川总面积退缩了72.302 km2，减少比例约为11.67%，年均变化速度为-5.16 km2/a。另外，从文献[26]得知，2010年念青唐古拉山脉西段冰川总面积为571.81 km2，就2003—2010年和2010—2017年2个时段而言，存在明显的变化差异，2003—2010年冰川面积的变化速率为-6.785 km2/a，而2010—2017年为-3.544 km2/a。

为探究不同海拔冰川面积变化，利用90 m分辨率的SRTM数字高程模型(digital elevation model，DEM)数据生成200 m等高线，统计相邻两等高线间冰川面积变化，见表1。由表1可以看出，海拔在5600～6200 m区间内分布的冰川数量最多，占比91%以上。从年均变化率上看，低海拔区域冰川年均变化率最大，退缩趋势明显，这与低海拔地区温度较高有关。随着海拔的升高，温度降低，冰川年均变化率逐渐减小。在海拔5600～5800 m区间内，冰川面积变化最大，达到-34.56 km2(图5)，该区间内包含有大量冰川边界，属于易消融区域。在6200 m以上高海拔地区，温度较低有利于冰川的发育，冰川面积不减反而缓慢增加。图5中白色网格表示海拔在5600～5800 m区间范围，底图为2017年冰川区Landsat影像。

表1 念青唐古拉山脉西段不同海拔冰川面积变化统计

2.3 冰川储量变化与不确定性分析

由第三极地区冰川表面高程变化数据产品获取的念青唐古拉山脉西段冰川高程变化如图6所示。由图6可知，冰川的大部分冰舌区域均存在不同程度的表面高程减薄现象，其中北坡减薄现象更加明显。经统计，2000—2015年念青唐古拉山脉西段冰川平均减薄速率为38 cm/a，利用该数据产品计算的冰川储量变化速率为-0.145 km3/a，相当于-0.123 Gt/a。为了与遥感估算的冰川储量变化进行对比，尽可能减小观测时段不一致带来的误差，通过计算得知，2002年8月至2015年12月念青唐古拉山脉西段冰川质量变化速率为-12.78 cm/a，利用式(5)转化为冰川体积变化为-0.088 km3/a，相当于-0.075 Gt/a

(5)

式中，ΔV为冰川体积变化；ΔH为GRACE计算的等效水高的变化；S表示两期冰川面积的平均值(583.154 km2)；ρg为冰川平均密度参数，取850 kg/m3。

卫星重力和遥感估计的冰川储量变化存在较大差异，其原因在于，GRACE重力信号中包含冻土融水、地下水、湖泊水等其他信号，这些信号的不确定性影响了GRACE观测结果的精度，使冰川储量的估计结果存在较大的误差[27]。此外，冰川表面高程变化的不确定性对遥感估计结果影响较大，本文采用的冰川储量估计方法依赖每个像元的高程变化，无论是SAR还是光学立体像对在地形起伏较大的区域获取的DEM精度普遍不高，即便使用了DEM配准方法减少了相位误差和一些系统差的影响，但残余误差对结果的影响依然很大。

图5 5600～5800 m区间冰川分布Fig.5 Distribution of glaciers in the 5600～5800 m interval

图6 念青唐古拉山脉西段2000—2015年间冰川表面高程变化Fig.6 Changes of glacier surface elevation in the western Nyainqentanglha Mountain from 2000 to 2015

3 结 论

利用GRACE、GRACE-FO重力卫星数据，结合GLDAS水文模型，反演了念青唐古拉山脉冰川质量变化，借助气象资料说明念青唐古拉山脉的常年质量亏损和区域温度升高有密切联系。此外，以念青唐古拉山脉西段冰川为重点研究对象，结合Landsat光学遥感数据和表面高程产品对冰川的面积和储量变化进行详细分析。主要结论如下。

(1) 2003—2019年间念青唐古拉山脉冰川质量整体亏损严重，经过空域法泄露误差改正后，念青唐古拉山脉冰川质量变化速率为-15.53 cm/a，其中念青唐古拉山脉西段和东段冰川质量变化速率分别为-20.23 cm/a和-10.82 cm/a。对夏季的温度变化和降水异常分析结果表明，念青唐古拉山脉的冰川退缩和降水减少、温度升高有关。

(2) 念青唐古拉山脉西段冰川在不同海拔处冰川面积分布有所不同。海拔5400 m以下分布大量小型冰川，这部分冰川消融速度最快；海拔大于6200 m的区域温度较低，有利于冰川的发育，冰川面积有所增加。海拔在5600～5800 m区间的冰川面积变化最大，其值为-34.56 km2，占总面积变化的48%，该区间内分布大量冰川边界和冰舌区域，在冰川中最易消融。

(3) 利用卫星重力和遥感技术获得的念青唐古拉山脉西段冰川储量变化速率分别为-0.075 Gt/a和-0.12 Gt/a，二者变化趋势相同，但数值差异较大，其主要原因是GRACE后处理误差和冰川表面高程数据精度不高导致的。

(4) 卫星重力和遥感技术的组合应用可以优势互补，能够充分发挥卫星重力在大、中尺度的快速探测能力和遥感技术精确的冰川信息提取能力，实现对变化剧烈冰川的快速定位和精准监测，但在冰川储量估计方面仍存在一定缺陷，有待进一步研究。