李 海，杨成生，惠文华，朱赛楠，张 勤

（1.长安大学地质工程与测绘学院，陕西 西安 710054；2.中国地质环境监测院（自然资源部地质灾害技术指导中心），北京 100081）

0 引言

山地冰川作为冰冻圈的重要组成部分，其变化与气候变化、区域水循环有着密切的关系[1−2]。对于属于干旱与半干旱区的冰川发育区，冰川融水作为其径流的主要补给源，是区域内可靠的水源保证和生态保障[3−4]。同时冰川的持续消融大大提高了冰湖溃决的风险，危及下游百姓的生命财产安全[5−7]。中国冰川集中分布在青藏高原及其周围的喜马拉雅山脉、天山、喀喇昆仑等山脉，各地冰川分布见表1[8]。研究表明，除喀喇昆仑山地区少量冰川以外，中国西部冰川普遍处于退缩、变薄状态，但是受气候和地形差异的影响，不同区域存在差别[9−12]。其中青藏高原地区冰川退缩速率总体呈现由边缘向内陆腹地、东部向西部变慢的规律[13]。冰湖是由冰川挖蚀、融化形成洼坑或者冰碛物堵塞致冰川槽谷积水而成的一类湖泊[14−15]。王欣等[16]分析了近30年来我国喜马拉雅山区冰湖变化的特征，剔除在气候变化的作用下，冰川变化是冰湖变化的主要影响因素。近年来喜马拉雅地区冰川退缩严重，一些以冰川融水为主要补给的湖泊，面积扩张较快，冰湖溃决危险性急剧升高，引起了研究者的广泛关注。

表1 中国各地冰川面积分布Table 1 Glacier area distribution in China

喜马拉雅地区是中国重要的冰川作用区，也是冰湖溃决等地质灾害的多发区[17−18]。在已有的15次冰湖溃决事件记录中，12次发生在本区。其中藏南的希夏邦玛峰地区就曾多次发生冰湖溃决灾害[19−20]，如1968年的阿亚错、1981年波曲的章臧布错和2002年聂拉木县北侧的嘉龙错，并且上述冰湖溃决引发的洪水和泥石流均造成多处人工设施损坏和人员伤亡[21−22]。鉴于希夏邦玛峰地区冰川广布、作用强烈、冰湖溃决灾害多的特点，文章利用陆地资源卫星（Landsat系列）影像，通过面向对象的阈值分类方法提取了希夏邦玛峰地区1994—2018年共9期冰川冰湖的面积，并探究了冰川冰湖的变化规律。对指导该地区的水资源调控和冰川冰湖灾害的防治有重要的参考意义。

1 研究区与数据

1.1 研究区概况

研究区位于西藏南部喜马拉雅山脉中段，内有希夏邦玛峰。该地区主要受到大陆西北季风和南亚西南季风环流影响，由于地势阻隔，山脉南北两侧气候相差较大。南侧夏季高温多雨，北侧相对寒冷干旱[23−24]。研究区山脉南北两侧皆以沉积构造带为主，地层疏松、板块活跃，在频繁的地震和冰川作用的影响下，滑坡、冰湖溃决等地质灾害频繁发生[25]。此外，研究区内有中国和尼泊尔重要边境通商口岸城市吉隆镇和聂拉木县，政治经济地位独特。

1.2 数据介绍

本次冰川冰湖的监测数据包括遥感数据、DEM数据、气象数据、冰川编目数据4种。其中遥感影像主要来源于美国陆地资源卫星（Landsat系列）TM5、ETM、ETM+、OLI8传感器影像，皆由美国地质调查局免费提供（https://earthexplorer.usgs.gov/）。首先，经过影像筛选发现该地区每年10月份前后的影像具有云量少、积雪少的特点。其次，以往研究表明1992年后该地区气温上升速率明显加快[26]，综合影像资源、云量、积雪等多方面因素，选取1994—2018年跨度24年9期，每期大致10月份的影像开展研究（表2）。冰川参考数据选用2013年伦道夫冰川编目数据(RGI, Version 3.2，2013)，基本反映2010年前的冰川分布情况。地形数据选用SRTM DEM，高程精度为16 m（95%置信区间）[27]。

表2 光学影像数据详情Table 2 Details of optical image data

对于气象数据的选取，由于该地区气象站点稀疏，仅聂拉木县设有气象站，但是聂拉木县位于山脉南侧峡谷之间，受南亚暖湿季风影响显著，不能准确反映山脉北侧冰川区域的气候[28−29]。全球网格气象再分析数据集，采用特殊的插值方法，充分考虑了地形、距离、海拔等因素对气候影响。因此采用来源于美国国家海洋和大气管理局（https://psl.noaa.gov/data/gridded/）提供的NOAANCEP陆地表面0.5°×0.5°网格气温再分析数据和CMP 0.5°×0.5°网格降水再分析数据。其中气温数据均方根误差为0.18 ℃，日降水量偏差小于1 mm，精度较高[30−31]。

2 数据处理

文章在面向对象的基础上采用阈值分类方法提取冰川冰湖面积。首先在ENVI5.3软件中对Landsat系列影像进行大气校正、影像融合、拼接等影像预处理工作。为了更准确的提取冰川边界，运用面向对象的阈值分类方法对影像进行了分割与分类工作，全部处理流程如图1所示。为了方便后期的交叉验证和人工修正，采用伦道夫冰川编目数据或前一期冰川提取矢量结果参与的多尺度分割方法。通过多次试验，最终确定分割因子为60、形状因子为0.4、紧致度因子为0.5。

图1 数据处理流程Fig.1 Data processing flow

对山体阴影区的分类主要采用目视解译进行。为了避免山体阴影对冰川分类结果的影响，首先对山体阴影进行提取，并在结合Google-Earth多期历史影像交叉验证阴影区冰川范围。利用山体阴影与冰湖具有相似的光谱响应特征，以及其在近红外波段的超低反射，并经多次试验，阴影提取采用归一化水指数NDWI>0.18、近红外光谱NIR<1 000。随后，采用NDWI>0.18对冰湖进行了提取。考虑到提取精度和实际意义，冰湖提取排除了面积小于3个像素冰湖。

采用波段比值法提取净冰川。根据冰川在红波段的高反射和短波红外高吸收的光谱特性，实验采用R/SWIR1>2.5进行冰川提取。然而研究区地面常存在的少量残留积雪，影响冰川提取的准确度。而以往的研究往往只采用2～3期影像对冰川进行长时间序列监测，时间间隔较长，无法发现积雪的变化以降低积雪的影响或者需要大量的人工目视判读（表2）。而本文收集了研究区1994—2018年共9期影像，时间间隔短，可以通过前后对比发现积雪的变化。在冰川持续退缩的大前提下，自动去除前时相冰川边界以外的积雪。其中1994年的冰川提取则参考RGI数据。由于RGI数据反映的是2000—2010年间的冰川状况，故仅以此作为参考并辅以人工目视判读，剔除积雪。2000年的冰川提取参考1994年冰川提取结果。以此类推提取了9期的冰川冰湖面积，结果如表3所示。冰川冰湖整体提取结果如图2所示。

图2 1994年冰川冰湖轮廓提取图（TM影像5、4、3合成）和单个冰川、冰湖9期变化结果Fig.2 Overall extraction of glacial and glaciers lake in 1994 (TM images 5, 4, 3 are synthesized),the results of 9-stage change of single glacial and glaciers lake

由于本次实验没有地面测量资料，因此冰川冰湖误差评估参照第二次冰川编目误差评估办法，即将空间分辨率的一半作为缓冲区域确定面积不确定度[11]。因此冰川冰湖面积变化误差е的公式（1）为：

式中：ε1、ε2——前后时相冰川、冰湖面积的误差。

由误差传播定律可知，面积变化率误差ε的计算公式为：

式中：x、y——分别指冰川、冰湖前后两期的面积。

表碛覆盖型冰川采用较为可靠的目视解译。研究区南侧受到南亚暖湿季风影响，属于温带海洋性冰川区，表碛覆盖型冰川广布。然而由于表碛冰川在光谱特性上和裸地相似，难以实现自动化提取。近年来，不少学者针对表碛冰川的自动化提取提出了许多方法，但是都存在准确度不高且需要专家目视修正的问题，不适用于变化分析[32−33]。鉴于该地区表碛覆盖型冰川较集中、面积大且数量少、退缩缓慢等特点，参考伦道夫冰川编目数据，以目视解译方式提取1994、2018年首尾两期的表碛冰川轮廓。1994年表碛冰川面积为（168.41±17.6) km2。至2018年，在部分净冰川退缩裸露出表碛冰川的形势下，表碛冰川面积增加至(223.97±19.48 ) km2。而表碛冰川的退缩主要表现为表碛冰川向末端冰湖演变，面积为11.03 km2。

3 分析与讨论

3.1 冰川冰湖总体变化分析

文章共提取了1994—2018年间共9期冰川冰湖的轮廓，面积变化如图3所示。结果显示24年间冰川总体上持续退缩，面积从(1 012.66±68.42) km2退缩至(699.89±61.50 ) km2、退缩比例为(30.83±7.66) %、年平均退缩速率为(1.28±0.32) %，冰川的退缩速率经历了加速减速的轮换。该地区的冰湖在24年间面积从(36.15±4.32 ) km2扩张至(50.65±7.09 ) km2、扩张比率为(47.23±25.79) %、年平均扩张速率为(1.889±1.07)%，冰湖扩张速率高于冰川退缩速率，并且其扩张的快慢变化与冰川退缩的快慢变化基本吻合，相关系数为−0.98。较高的相关性说明该地区冰湖扩张水源基本来源于冰川融水，间接证明了冰川、冰湖的提取方法可靠。

图3 24年冰川冰湖面积变化Fig.3 Changes of glaciers and glacial lakes area in 24 years

3.2 冰川变化与气候变化的响应

文章探究了该地区冰川冰湖变化与气候变化的响应关系。研究表明，冰川冰湖变化和当地的气候变化存在显著的关联，其中气温和降水影响最大[34]。

用1992—2018年近27年的夏季气温和全年气温变化分别与冰川面积变化作比较。为了防止因部分异常气温影响拟合结果，采用鲁棒性更好的稳健估计法拟合气温变化。从图4（a）可以看出全年平均气温和夏季平均气温均呈上升趋势，与冰川面积变化呈负相关，这可能是该地区冰川退缩速率加快的重要原因。夏季平均气温以每年0.035 ℃的速率上升，并从 2010年以来上升了1.8 ℃左右，为冰川在夏季的消融创造了良好的条件，这也导致2010年后冰川退缩速率加快。

图4 研究区气温和降水变化Fig.4 Variation of temperature and precipitation in the study area

同时，分析了该地区1990—2019年的近30年的降水变化，将其细分为夏季、冬季、全年降水总量三组。如图4(b)所示，该地区降水主要集中在夏季，且全年降水和夏季降水均呈现不断上升趋势。与此同时，夏季降水以15.8 mm/a的速率上升，趋势明显，但冬季降水较少且无明显上升。30年间降水量增加了约500 mm，年增量为17.2 mm，且大部分来自夏季降水，对应冰川在消融季的退缩。而春冬冰川物质积累季，降水并无明显增多。综上所述，夏季气温和降水的升高是本地区冰川的退缩加快的重要原因。与此同时降水增加和冰川的加速退缩，共同促进冰湖的快速扩张，提高了该地区冰湖溃决的风险。

3.3 影响冰川变化的其他因素

为了进一步研究该地区冰川变化的规律，本文进一步讨论了冰川面积的大小、冰川经纬度位置、流向多个因素对冰川变化的影响。首先利用SRTM DEM提取研究区域内的山脊线作为分冰岭，将冰川与冰川之间按山脊线分开，划分成独立的冰川。山脊线主要采用水文分析法和表面分析法相结合的方法自动化提取[11]。由于DEM精度问题，部分山脊线会偏离山脊位置，需要少量人工修定。整个冰川区经过分冰岭分割之后，大约有2 000个冰川，其中细碎冰川较多，分别在1994、2018年的结果中选择其中100个冰川作为统计样本，从面积大小、冰川纬度、经度位置、流向四个角度进行统计，统计结果如表3所示。

表3 9期冰川冰湖面积表Table 3 Area of glaciers and glacial lakes in Phase 9

通过统计不同大小冰川的变化率，探究冰川大小对冰川变化的影响。根据冰川的数量和面积特征，将100个冰川按面积大小分为6组，分别是0～1 km2、1～2 km2、2～3 km2、3～4 km2、4～7 km2、7～20 km2。统计结果如表4所示，随着面积的增大，冰川的退缩速率变小，说明：①面积较大的冰川应对外界气候变化更具有稳定性，同时面积较大的冰川对气候具有较大的调节作用；②小面积的冰川冰量损失极为严重，因此随着冰川的持续退缩，相同的气候的变化将带来更严重的退缩。

表4 不同角度下冰川的退缩率Table 4 Retreat rates of glaciers with multipectives

通过统计不同经纬度位置的冰川的变化率，探究位置对冰川变化的影响。从经纬度位置来看，研究区北侧冰川和西侧冰川退缩速率较慢，这与其他研究者获得的两大规律基本吻合：①喜马拉雅山脉冰川退缩速率整体呈现由东侧向西侧逐步降低；②喜马拉雅山脉冰川退缩速率外侧高于内侧[35−36]。该结果与地形差异而导致的气候差异有着密切的关系。研究区受南亚印度洋温湿季风和大陆西北干冷季风双重影响[37]，由于山体的阻隔印度洋季风由南向北迅速减弱，而由西北而来的大陆性季风因希夏邦玛山脉的阻隔向东减弱。致使冰川退缩速率出现了南快北慢、东快西慢的现象。

通过统计不同流向冰川的变化率，探究流向对冰川变化的影响。将100个冰川大致按北、东北、东南等6个方向分类。冰川变化结果如表3所示。结合图1，北向冰川多位于研究区南侧，冰川面积小，南亚季风影响严重，退缩较大。东北流向冰川主要分布于希夏邦玛峰东北侧，该地区冰川退缩最慢，其次是西北流向冰川，主要位于右侧的错朗玛冰川群北侧，主要原因是希夏邦玛冰川群相对纬度高，气温相对较低。东南流向冰川主要分布在错朗玛冰川群的南侧地区，西南流向冰川为希夏邦玛峰西南侧冰川，两者皆为南亚季风的迎风坡，但是西南流向冰川较慢可能是左侧存在山脉阻隔，南亚季风较弱所致。综上所述，经纬度位置和流向对冰川变化的影响，主要原因在于局部气候的差异。

因此根据本文研究结果，提出以下建议：①持续关注气温和降水对冰川的影响；②对面积扩张严重的高危冰湖进行更细致的监测；③走健康绿色发展道路，减缓全球变暖的速率。

4 结论

本文利用Landsat系列影像对冰川变化显著的希夏邦玛峰区域进行了长时间序列的冰川冰湖变化规律研究。结果表明：①1994—2018年该地区冰川面积退缩比例为30.83%，冰湖面积扩张比率为47.23%，冰湖扩张严重；②该地区冰川退缩速率加快的主要原因在于夏季气温的升高；③随着冰川面积的逐渐减小，相同的气温变化将带来更加严重的冰川退缩；④南亚季风、西北大陆季风、地形等多重因素致使该地区冰川退缩速率整体呈现南快北慢、东快西慢的现象。