基于多源遥感数据的珠峰自然保护区冰川监测研究

2023-11-18张翠兰丁悦凯曹香芹梁文莉

干旱区地理(汉文版) 2023年10期

冀琴，张翠兰，丁悦凯，曹香芹，梁文莉

（重庆师范大学GIS应用研究重庆市高校重点实验室，重庆 401331）

冰川变化与气候变化息息相关，对气候系统有重要影响［1］。研究表明冰川变化对温度响应程度较高，其变化特征能反映出全球气候变化特点［2］。随着全球气候变暖，尤其是21 世纪以来，冰川退缩趋势愈加显著。一定时期内冰川消融有利于水资源补给［3］，而长期的冰川消融不仅会引发冰湖溃决、冰川泥石流等自然灾害，还会导致海平面上升，威胁沿海居民的生命安全和社会经济发展。因此，研究全球变化背景下的冰川变化，对认识区域冰川变化及其影响、评价区域的生态环境、探究冰川与气候的相互作用关系具有重要科学与社会意义，冰川面积的分布与变化特征也一直是科学研究热点［4-5］。尽管冰川面积分布与变化特征研究对冰川区的灾害预测和防治、水资源利用等具有重要意义，但物质平衡作为冰川变化重要指标，可以直观反映区域气候变化特征，也是连接冰川与水资源的重要枢纽［6］。因此，研究冰川面积的同时监测冰川物质平衡，有助于深化对冰川响应气候变化的科学认识，增加对冰冻圈生态系统演变的理解，为全球气候变暖下的环境应对策略提供科学参考［7］。

传统的冰川监测以实地考察为主，精度高，但耗时耗力，难以实现大范围动态监测。遥感技术的发展使冰川研究进入了崭新阶段，光学遥感数据因其直观、清晰、易于判读的优点被广泛应用于冰川面积研究中。而现阶段冰川表面物质平衡的变化监测，大多采用物质平衡模型和多源数字高程模型（DEM）空间匹配的大地测量法。尽管物质平衡模型在冰川物质平衡研究中发挥了重要作用，但模型存在参数不确定性和区域适用性问题。目前，多源DEM 空间匹配的大地测量法发展成熟且被广泛应用，但存在多源数据匹配问题。冰川表面形变作为冰川物质平衡的重要参数，可直接参与反演冰川物质平衡计算，其变化趋势能直观反映冰川物质平衡变化。因此，研究冰川表面形变特征可在一定程度反映冰川物质平衡变化。近年来，雷达遥感尤其是差分干涉测量短基线时序分析（SBAS-InSAR）技术已经在毫米级地表形变应用中取得了多项标志性成果［8-9］，其高精度观测、较强适用范围以及可长时间监测缓慢地表形变的优势［10］，为冰川表面物质平衡研究提供了新的途径。

本文以喜马拉雅山珠峰自然保护区冰川为研究对象，采用Landsat系列光学影像和哨兵1号（Sentinel-1）雷达数据，运用比值阈值法和SBAS-InSAR技术，监测珠峰地区冰川面积分布、变化特征以及冰川表面形变，分析冰川表面形变的时空分异特征，同时结合前人研究成果，探究冰川物质平衡与气候变化的响应关系，以期为珠峰保护区冰川灾害的防治与治理提供科学依据。

1 研究区概况

珠峰自然保护区地处中国和尼泊尔交界处（图1），总面积约3.38×104km2，是世界上海拔最高的自然保护区［11］。据中国第二次冰川编目显示［3］，保护区内共有冰川1049 条，总面积约1450.40 km2，其中绒布冰川分布面积最大。研究区为大陆性高原气候，光照充足，冬寒夏凉。嘎尔气象站资料显示，区域年均气温2.1 ℃，日照百分率75.3%，年均降雨量270.5 mm，雨热同季。珠峰自然保护区因其独特复杂的地理位置和气候条件，成为科学研究关注的热点区域之一。

图1 研究区概况Fig.1 Overview of the study area

2 数据与方法

2.1 数据来源

2.1.1 遥感影像本文使用的Landsat系列遥感影像来源于美国地质调查局（USGS）（https://earthexplorer.usgs.gov/）和地理空间数据云（https://www.gscloud.cn/）。为提高冰川边界提取的精确度，尽量选取少云、少积雪的遥感影像，最终使用的影像如表1所示。Sentinel-1 卫星是欧洲航天局为监测全球尺度陆地、海洋地形和极地等设计的双星观测卫星，搭载能提供连续图像的C波段合成孔径雷达。其升轨影像集中在夜间成像，降轨影像在白天成像。由于冰川热融可能导致影像失相关，故选择更稳定的升轨数据［12］。本研究选取2020 年1 月1 日到12 月26日共72 景升轨数据（表2），轨道文件采用欧空局提供的精确轨道数据（https://e4ftl01.cr.usgs.gov/）以消除大气相位影响和平地效应［13］，轨道数据精度在10 cm以内。

表1 遥感影像数据Tab.1 Remote sensing image data

表2 Sentinel-1A数据Tab.2 Sentinel-1A data

2.1.2 DEM DEM 数据采用ASTER GDEM 和SRTM DEM 2 种，分辨率均为30 m。2 种DEM 均可用于提取冰川边界，但SRTM受镜面反射、回波滞后和雷达阴影影响，数据会出现缺失和异常［14］。而ASTER GDEM 数据缺失较小，可较好反映地貌特征。此外，龙四春等［15］基于不同参考DEM 进行DInSAR 技术处理，结果表明SRTM 的差分干涉图的噪声较少。故本研究将ASTER GDEM 用于冰川边界识别以及海拔和坡度计算，将SRTM DEM作为In-SAR技术参考DEM。

2.2 研究方法

2.2.1 冰川边界提取本研究选用较为精准有效的比值阈值法（band3/band5），提取珠峰自然保护区的冰川边界。对比多次实验结果，确定TM/ETM+影像比值图像的最佳阈值为1.8，OLI影像为1.0。基于决策树分类得到二值图像，最后进行目视解译。同时，基于第二次冰川编目数据对结果进行修正，从而提高冰川边界解译精度（图2）。

图2 1990—2020年部分区域冰川变化Fig.2 Part of glacier variations from 1990 to 2020

2.2.2 冰川总面积计算采用年均变化率评估珠峰自然保护区冰川总面积变化情况，计算公式如下：

式中：APACi为年均变化率；ΔSi为第i时段内冰川面积变化量；Si为第i时段内初始年份冰川面积；ΔTi为第i时段的时间间隔。

2.2.3 冰川形变速率监测Berardino等［16］于2002年提出小基线集技术（SBAS-InSAR技术），该方法无需考虑主影像选择等问题，对时空失相干和大气等因素不敏感，且监测精度达到毫米级别。其原理如下［17-18］：

设研究区在时间段t0～tn内有N+1幅影像，通过基线组合得到M景干涉图像。

设ta和tb时刻（b＞a）得到的图像干涉处理后生成第i幅干涉图，像元（x,y）处的干涉相位表示为：

式中：d(ta,x,y)和d(tb,x,y)为相对于参考时刻t0、ta和tb时间点对应LOS 向地表形变量；φa和φb为对应时间的解缠相位；λ为雷达波长。

解缠后的干涉图相位为：

式中：δφ为干涉相位；Aφ为系数矩阵；(x，y)为像元位置。

本研究利用SBAS-InSAR 技术对Sentinel-1 升轨数据进行处理并获取形变结果，数据处理过程如图3所示。

图3 差分干涉测量短基线时序分析（SBAS-InSAR）技术流程图Fig.3 Workflow diagram of SBAS-InSAR technology processing

2.3 精度评价

以20 m 为距离对1990—2020 年冰川边界做缓冲区处理，用缓冲区面积除以冰川总面积得到误差率（表3）。最终得到的误差率较小，可以满足精度要求［19］。

表3 1990—2020年珠峰自然保护区缓冲区面积及误差率Tab.3 Buffer area and error rate of Qomolangma Nature Reserve from 1990 to 2020

3 结果与分析

3.1 冰川总面积分布与变化特征

1990—2020 年珠峰自然保护区冰川面积整体呈退缩状态，由1306.31 km2减少到1059.15 km2，共退缩247.16 km2，面积变化率为-18.92%。不同时期的冰川退缩程度有所差异（表4）。其中，1990—2000 年冰川面积退缩了11.48 km2，年均变化率为-0.09%·a-1；2000—2010年冰川面积退缩63.51 km2，年均变化率达-0.49%·a-1；2010—2020 年冰川面积退缩了172.17 km2，相较于前2 个时段，该时段冰川退缩速率显著升高。综上，近30 a来，珠峰自然保护区冰川呈加速退缩趋势，且近10 a退缩趋势愈发显著，与喜马拉雅山地区冰川退缩趋势一致［20］。

表4 1990—2020年珠峰自然保护区冰川面积分布和变化Tab.4 Distribution and change of glacier area in Qomolangma Nature Reserve from 1990 to 2020

3.2 不同海拔的冰川分布与变化特征

结合DEM 数据，以200 m 为间隔，统计1990—2020 年珠峰自然保护区冰川面积在随海拔的分布情况（图4a）。随着海拔升高，冰川面积分布均呈现先增后减，且各时期分布面积均在5600～5800 m 处达到最大。以1990 年为例，海拔3600～4800 m 和7000～8200 m处冰川分布面积较少，分别为14.35 km2和12.77 km2；5600～6000 m范围内面积分布较多，达445.51 km2，占冰川总面积的34.10%。海拔4800 m以下冰川面积分布较少，可能是由于温度影响且水汽补给较少，导致低海拔呈干冷气候，不利于冰川的发育与积累；5200～6400 m 处有良好的空间位置和充沛的降水，气温随海拔升高而降低，为冰川的发育提供适宜条件。冰川界认为一座山峰可能存在一个或者更多降水带［21］。高海拔山地冰川的降水量随高度变化有明显变化规律，随着海拔升高，山坡的湿润空气温度下降，饱和度增加，当到达一定高度以降水形式下落，形成第一个降雨带，高海拔地区的冰川和积雪具有较低温度和充足的水汽，可降低过境云团温度，从而形成降水，即第二大降水带。谢自楚等［22］认为受局部环流影响，珠穆朗玛峰地区存在较明显的“第二大降水带”。日本学者Yasunari等［23］对喜马拉雅山研究发现第二大降水带分布在海拔5000 m 以上。因此，5200～6400 m 范围的冰川积累可能是依赖于第二大降水带的补给。

图4 1990—2020年珠峰自然保护区不同海拔冰川面积分布和变化Fig.4 Distribution and variation of glacier area at different altitudes from 1990 to 2020

1990—2020 年珠峰地区冰川面积总退缩量变化表现为先增后减（图4b），在5400～5600 m 处达到最大，为50.06 km2。海拔8000 m 以上几乎无退缩。相比于退缩量，研究区内冰川总退缩率则呈波动下降态势，在5000～5200 m 处出现小高峰，达30.93%。可以发现：低海拔区冰川退缩量较少，但退缩率较高；而高海拔处冰川退缩量和退缩率较低。可能是因为低海拔处温度较高，不利于冰川积累，导致冰川虽然退缩量小，但在较少冰川分布面积影响下，形成高退缩率。而高海拔区温度相对较低，冰川不易消融，但积累空间较小，因此退缩率和退缩量均不明显。

3.3 不同坡度的冰川分布和变化特征

以5°为间隔，将冰川坡度分为12 个等级。4 个时期冰川面积分布特征相似，随坡度增加，冰川分布面积均表现为先增加后减小趋势（图5a）。以1990 年为例，冰川在10°～15°面积分布最大，达196.21 km2，占对应时段冰川总面积的14.94%，55°～60°分布最少，仅为21.11 km2，面积占比为1.62%。可见，珠峰自然保护区冰川主要分布在较平缓地带，极平缓和陡峭区域分布较少。近30 a来研究区冰川总退缩量和总退缩率均呈现单峰形态（图5b），在10°～15°坡度范围内退缩面积最大，达37.95 km2，55°～60°范围内退缩量较小，仅为2.53 km2。而退缩率在20°～25°处最大，55°～60°范围内最小。

图5 1990—2020年珠峰自然保护区不同坡度冰川分布特征和变化特征Fig.5 Distribution and variation characteristics of glacier area at different slopes in Qomolangma Nature Reserve from 1990 to 2020

3.4 冰川形变特征

基于SBAS-InSAR 技术获取2020 年1 月1 日至12月26日珠峰自然保护区冰川表面形变信息，结果表明2020 年珠峰自然保护区冰川平均形变速率介于-129.069～140.252 mm·a-1之间（图6）。

大部分冰川形变呈现下降，仅少部分区域表现为抬升。选取a、b、c 和d 4 个区域，分析可知：a区域形变速率达-129.069～113.944 mm·a-1，b 区域形变速率为-32.811～36.114 mm·a-1，c 区域形变速率为-37.397～36.357 mm·a-1，d 区域形变速率为-21.033～32.013 mm·a-1。可见，大部分冰舌区域呈减薄趋势，而高海拔区呈相反趋势。

由于缺乏实地测量数据，本研究利用非冰川区形变结果对冰川形变结果进行间接验证，在非冰川区域选取4个裸地区（图7），统计形变量（表5）。计算得到均方根误差分别为0.16 mm、0.16 mm、0.17 mm和0.17 mm，证明形变结果具有较高的可靠性。

表5 裸地区域形变统计Tab.5 Deformation statistics of bare land area

图7 裸地区域选取Fig.7 Selection of bare ground areas

3.4.1 不同海拔冰川表面形变特征基于ArcGIS软件，利用分区统计工具获取不同海拔的冰川形变量（图8）。结果表明：不同高程带冰川呈现出不同形变趋势，研究区内绝大部分冰川呈现减薄态势。海拔5800 m 以下冰川表现不同程度沉降趋势，其中4200～4400 m处沉降最为严重，变化量达-17.246 mm。海拔5800～6400 m出现轻微抬升，并在6200～6400 m处抬升较大，形变量为0.347 mm。上文中提及珠峰地区在5200～6400 m 处可能存在“第二大降水带”。根据冰川形变结果分析，冰川在5800～6400 m 呈现抬升，结合前文对冰川面积的探讨，由此推测，珠峰自然保护区的“第二大降水带”大致位于5800～6400 m处，与Yasunari等［23］在喜马拉雅山发现“第二大降水带”的观测事实相符。6400 m以上冰川沉降量随海拔升高波动增加，在8000～8200 m 处形变量达-12.448 mm。低海拔区域受气温较高、湿度较小等因素影响，冰川出现消融，导致冰川减薄。而通常高海拔地区由于空间充足，降水充沛，冰川更容易积累，从而表现抬升。但本研究中高海拔区域冰川沉降严重，推测可能是研究区发生过雪崩导致形变异常。经查证，研究区包含的希夏邦马峰地区曾于2020 年8 月发生过较大雪崩。总体而言，冰川在低海拔地区更容易出现物质消融导致冰川减薄，而高海拔多呈现抬升趋势。

图8 2020年珠峰自然保护区不同海拔冰川形变特征Fig.8 Characteristics of glacier deformation at different elevations in Qomolangma Nature Reserve in 2020

3.4.2 不同坡度的冰川表面形变特征以5°为间隔，将冰川分为13个等级（图9），珠峰自然保护区冰川在各个坡度呈现不同形变程度。坡度0～50°内，研究区冰川呈现沉降态势，形变量随坡度增大呈“减-增-减-增”趋势，其中40°～45°范围内，沉降最为严重，形变量达-1.401 mm。当坡度大于50°后变化趋势减小，在坡度55°左右开始出现抬升趋势，并于大于60°处达到形变量最大值，为1.545 mm。出现上述情况的可能原因是：陡坡一般位于高海拔区，温度较低，冰川消融速率较慢，同时由于积雪的作用，导致冰川积累，抬升趋势显著。

图9 2020年珠峰自然保护区不同坡度冰川形变特征Fig.9 Characteristics of glacier deformation at different slopes in Qomolangma Nature Reserve in 2020

4 讨论

冰川面积变化对气候的响应具有一定滞后性，但气候变化规律可以透过冰川物质平衡直观呈现出来。近几十年来青藏高原气温持续上升，且增温速率大于北半球中纬度区域［24］。但近10 a青藏高原东部增温趋势开始减缓，而中西部地表增温较为强烈［25］。任贾文等［26］研究表明喜马拉雅山中段气温不断波动上升且在夏季升温更为明显。过去几十年珠峰地区的增温速率也显著大于中国和世界的平均水平［27］。青藏高原降水变化受西风、南亚风及高原季风影响，喜马拉雅山中段降水主要源于印度季风，但研究区绝大部分冰川位于喜马拉雅山北坡，由于高大的山体阻挡了季风降水，导致北坡降水减少。此外，冰芯积累量可以直接反映降水量变化。侯书贵等［28］基于远东绒布冰川冰芯记录指出珠峰地区冰川积累量不断减少，与冰川退缩趋势一致。段克勤等［29］根据达索普冰芯积累量研究推测，喜马拉雅山中部地区的季风降水将随气温上升而减少。本研究中，珠峰自然保护区冰川面积在近30 a来不断退缩，且2020 年冰川整体出现减薄，物质平衡呈亏损态势。其呈现的物质平衡亏损状态与前人研究一致［30］，且从时间尺度看，珠峰地区冰川消融在进一步加强。研究区冰川变化趋势与近年气温升高和降水减少具有一定相关性，因此，气温增加、降水减少可能是导致冰川退缩的主要因素。

除气温和降水影响外，冰川所处地理位置、地形差异等因素，会导致冰川接收到的太阳辐射、水汽来源和季风强度不同，进而影响冰川变化。由于研究区降水受印度季风影响，研究区范围较大，而季风降水具有空间差异性［31］，在一定程度上会导致物质平衡差异。此外，冰川退缩与冰湖扩张之间存在密切联系，珠峰自然保护区冰川加速退缩，冰湖迅速扩张［32］。温度升高导致冰川出现一定程度消融，冰川末端表现尤为显著，而末端冰川消融会导致积雪融水增加，促进末端冰湖扩张。同时，温度上升导致冰湖蒸发量增加，当蒸发量远小于融水的补给量，冰湖继续扩张。可见，冰川消融对冰湖扩张有显著影响，同时冰川融水也是冰湖补给的重要来源［33］；而冰湖不断扩张对冰川末端的消融产生正向作用，导致冰川末端加速融化。

珠峰自然保护区内分布大量表碛覆盖型冰川，相比裸冰区域，表碛覆盖型冰川表面反射率较低、吸热快，可向表碛下部冰川和周围冰川传递较多热量，促进冰川消融。厚度较薄的表碛覆盖对冰川消融有一定促进作用，大于一定厚度时有隔热作用［34］。黑碳、粉尘和气溶胶等可通过干/湿沉降富集在冰川表面，当冰川表面污染，表碛覆盖面积进一步扩张，冰川反照率继续下降，吸热能力增强，加速冰川消融。

综上所述，气温升高、降水减少可能是导致珠峰自然保护区大部分冰川物质亏损的主要因素。同时，气候空间差异、地形等因素也可能是导致物质平衡差异的重要因素。