冰川湖泊(简称冰湖)，既是全球生态重要的水资源，同时也是诸多冰川次生灾害发生的诱因^[1-7]。在中国，冰湖主要分布在喜马拉雅山脉、昆仑山脉、天山山脉以及藏东南等地区^[8]。自19世纪中叶以来，青藏高原气温增长趋势已超过北半球同纬度地区，造成了青藏湖区的数量和面积发生了较大改变^{[4, 7]}。在全球气温变暖的背景下，近些年来，因冰川湖泊溃决所导致的洪水、涌浪、泥石流等灾害时有发生。据不完全统计，仅西藏地区近80 a来，因冰湖溃决导致的重大灾情事件就高达27次^[9]；进入本世纪后，因冰湖溃决引发的灾害呈现频次增加和规模增大的潜在趋势^{[8, 10]}，给下游人员和基础设施安全造成巨大的威胁。

因冰湖多分布于高海拔偏远山区，现场勘察工作存在难度大、成本高的问题。为此，国内外学者基于光学遥感影像解译和计算机辅助判识的思路开展了大量系统深入的研究^{[7, 11-13]}。常用的技术途径一方面为人工干预的冰湖监督分类^[14-15]，涉及单波段影像阈值分割、归一化差异水体指数、波谱关系法等方法；另一方面为基于机器学习的非监督自动分类，包括聚类分析、神经网络、支持向量机和决策树等提取方法^[15-18]。随着国内外高分辨率卫星系统的快速发展，为冰湖的编目等相关研究提供了更多可选择的数据源，在危险性冰湖的识别与动态监测方面发挥了重要作用^[16]。然而，基于光学影像的冰湖专题分析受气候环境的制约，尤其是中国藏东南和横断山脉等海洋型冰川分布区域受印度洋西南季风、太平洋东南季风和西风带影响，常年云雾较多、潮湿多雨，可用数据极为匮乏，数据源瓶颈的问题尤为凸显。

近年来，卫星合成孔径雷达(synthetic aperture radar，SAR)传感器因其不受云、雾、雨等恶劣天气条件的限制，获得各行业部门的广泛关注和高度重视。世界各国正在竞相研发和部署新型SAR卫星平台，如德国宇航局的TerraSAR-X/TanDEM-X卫星编队、意大利宇航局的COSMO-SkyMed四星星座、日本航天局的ALOS/PALSAR-1/2卫星系统、加拿大太空署的RadarSat-1/2卫星系统、欧空局的Sentinel-1A/B卫星星座以及中国的GF-3卫星系统等，能够为区域性灾害防治和预警科研等提供稳定可靠的数据支撑。随着国内外相关研究的不断深入，基于SAR遥感的大范围地表位移、形变与地表覆盖变化监测已展现出巨大的应用潜力^[19-20]。

目前，国内外基于SAR遥感的冰川运动、消融监测等研究已取得了长足进展^[20-21]，而SAR数据应用于冰湖发育、溃决等方面的研究仍处于起步阶段，相关科研成果尚不多见。在洪灾和水域变化监测领域，基于SAR极化分类或强度阈值分割实施水域提取和动态评估的相关应用已有先例借鉴^[22-24]。但极化数据的获取极为有限，长期、全面的系统研究目前仍难以开展；基于强度信息的分类亦由于高斯白噪声的存在而无法回避阈值难以确定、分类精度低等问题^[24]。

针对上述问题，本文提出了基于时序SAR影像的强度归一化比值实施冰湖提取与动态监测的方法。首先选取SAR影像时间域强度均值作为参考基准，然后通过SAR影像强度图的标准化比值开展冰湖分类提取和时序变化监测分析。为验证方法的可行性，选取位于贡嘎山西坡的贡巴冰川作为典型实验区，基于日本航天局的ALOS/PALSAR-1与欧空局的Sentinel-1A影像序列获取了贡巴冰川末端冰湖2007-2018年间的时空变化，并结合实地调研和无人机影像对实验对象近年来的变化规律进行了分析。

1 研究区域和实验数据

贡巴冰川位于青藏高原东南缘横断山大雪山主峰贡嘎山的西坡，属季风海洋性冰川，是环绕贡嘎山西坡最大的山谷冰川^[25]。17世纪初至19世纪中期，贡巴冰川分离成大、小贡巴冰川，如图 1(a)所示，其中大贡巴冰川全长10.6 km，面积19.1 km²；小贡巴冰川全长6.46 km，面积5.73 km^2[26]。根据1982-06至1983-06的中国科学院青藏高原综合科考实测数据^[27]，贡巴冰川雪线高度为4 900~5 100 m，是典型的表碛覆盖型冰川，消融区内表碛发育明显，延伸至雪线以下的冰瀑底部，厚度可达1 m^[28]。在贡巴冰川冰瀑底部至冰川末端发育了较多的冰湖，通过无人机影像调查，仅小贡巴冰川海拔4 300 m至冰舌末端0.38 km²内就发育了21个大小冰湖。参考气象观测资料，贡巴冰川末端多年平均气温为2.4℃，低于东坡的海螺沟冰川4.4℃，年降雨量约1 057 mm，降雨多集中于6月至9月间，夏季降雨集中而且多为暴雨，冰川的融水和降雨补充经冰川侧碛垄外沟谷和地下暗河流入田湾河后，最终汇入大渡河^[26-27]。

大、小贡巴冰川在消融退缩过程中，在冰碛物所构成的终碛垄和大小贡巴冰川末端之间储水发育了两个较大的冰碛阻塞湖，总面积超过120 000 m²，如图 1(a)、1(c)、1(d)中的L₁和L₂所示。经现场调研判断，冰湖L₁的湖水多为小贡巴冰川雪冰融水补给，而冰湖L₂的湖水则为大贡巴冰川融水径流因冰碛垄阻塞形成。

为获取贡嘎山贡巴冰川末端冰湖L₁和L₂的蓄水和溃决时序变化，选取两个SAR卫星系统自2007-01至2018-11期间获取的103景SAR卫星影像开展实验。其中，包括2007-01至2011-03期间由ALOS/PALSAR-1卫星获取的20景SAR影像(方位向分辨率为3.18 m，距离向分辨率为4.69 m)，以及2014-11至2018-11期间由Sentinel-1A卫星获取的83景SAR影像(方位向分辨率为13.99 m，距离向分辨率为2.33 m)，各卫星影像覆盖范围如图 1(b)所示。

使用SRTM(shuttle radar topography mission)V3.0 30 m分辨率的数字高程模型(digital elevation model, DEM)作为参考；此外，分别于2018-06-09和2018-10-16前往贡巴冰川现场调研，使用无人机进行航空影像采集，获取了贡巴冰川末端6 cm分辨率的正射影像和DEM。

2 SAR影像强度标准化比值提取冰湖技术方法

为了克服冰湖时序变化监测中SAR影像因斑点噪声的影响而导致阈值难以确定、分类难度高、分类精度低的缺点，提出了基于时序SAR影像的强度标准化比值开展冰湖提取与动态监测的方法，其技术流程如图 2所示。

2.1 基于多时相滤波的公共标准影像获取

在解译SAR影像前，需要完成影像的预处理，主要包含以下步骤：PALSAR影像聚焦和Sentinel-1A影像的镶嵌处理、SAR影像辐射校正、SAR影像集子像元级配准、SAR影像裁剪和SAR影像地理编码等。

卫星SAR成像过程是在雷达斜距-多普勒(range-Doppler，RD)坐标系下，完成回波信号的记录与处理，而RD坐标系是一个由卫星飞行方向(即方位向)和雷达波侧视方向(即距离向)构成的二维平面坐标系。具体到某一地面目标，每个雷达分辨单元由许多散射体组成，而这些散射体对于入射的雷达波都会产生后向散射，各散射体后向散射波相干叠加的结果将会在SAR影像上产生斑点噪声。由于斑点噪声的存在，单个像素的后向散射强度和实际地物分布特性表现为失相关；对于时序SAR影像集，由于斑点噪声的系统误差特征，单个像元在时间分布上呈现随机性，无论是干涉处理还是地物分类，斑点噪声的存在往往会导致边缘可见度降低或出现局部噪声。本文基于SAR影像的强度标准化比值对冰湖变化进行检测，为降低斑点噪声的影响，关键是要确立一个排除了斑点噪声影响的公共标准影像。

现有研究主要采取滤波的方法降低斑点噪声的影响，从而获取相对可信的强度图，常用的技术途径包括空间域(如均值滤波、中值滤波、Lee滤波等)、频率域(如Goldstein滤波、小波滤波等)和时间域多时相滤波^[29]。尽管通过滤波可高效消除局部斑点噪声，但也会导致不同程度的空间分辨率降低或边缘失真。为此，本文首先在目视解译的基础上筛选出冰湖未蓄水的影像子集，然后通过多时相滤波方法逐像元过滤失相关信号，最后计算各像元相关信号的均值，并作为公共的标准影像。

2.2 基于散射强度比值的冰湖提取

雷达卫星成像时，由于入射的雷达波在水体表面发生镜面反射而呈现较低的散射强度，现有研究多采取阈值分割标记水体(如Otsu算法、熵阈值法等)^[30]，从而实现水体区域的分离。然而在SAR影像中，并非低散射强度的区域都是水体，且对于同一地物单元，SAR影像中的后向散射信号不仅受卫星入射角的影响，同时也受地形的影响，使用雷达信号的强度分离水体容易造成误判；再者，斑点噪声的存在也给水体分割阈值的精确确定产生了较大的影响。

针对上述问题，本文在同一强度公共标准S_ij下，通过不同时间获取的雷达影像强度信号T_ij与其作商，如式(1)所示，获取标准化强度比值R_ij，进而检测各成像时刻间的时序变化。

$ R_{i j}=S_{i j} . / T_{i j}(i=1, 2, 3 \cdots ; j=1, 2, 3 \cdots) $

(1)

式中，./表示影像矩阵对应像元相除。

除了可以有效消除噪声的影响外，标准化强度比值的另一个突出优势是可以抵消因地形变化引起的信号分量，即消除因地形不同导致的散射强度变化(包括叠掩、阴影和透视收缩)。对于周期性蓄水和溃决的冰湖，使用溃决后裸露湖床的影像集经多时相滤波处理后，作为公共标准影像，能够保证蓄水后信号变化效应的放大。而对于长期蓄水的冰湖区域，可采用冰湖周边散射信号对湖区内插来作为公共标准影像。

本实验针对贡巴冰川末端冰湖的周期性蓄水和溃决的特点，基于未蓄水影像集经多时相滤波结果后的均值作为公共标准影像，以式(1)对时序SAR影像各像元的强度值进行比值运算，计算各像元的标准化强度比值。

3 冰湖水体的时序分类提取

通过对时序SAR影像配准、裁剪、滤波和目视解译后，确认了冬季L₁和L₂均未蓄水时的数据子集，经多时相滤波后的公共标准影像如图 3(a)所示。经计算，公共标准影像和时序SAR影像散射强度的比值如图 4所示。由于雷达对于水体比较敏感，蓄水前后的冰湖区域比值变化非常显著。

以2017-07-23获取的Sentinel-1A影像为例，该期比值影像如图 3(b)所示，贡巴冰川末端冰湖L₁和L₂蓄水后，在水体和非水体区域对比明显，P-P′剖面显示，非水体区域在(0, 2)值域内波动，而水体区域集中于(5, 33)值域内。分别选取2015-2018年期间冰湖L₁和L₂同时蓄水的比值影像提取P-P′剖面信息，如图 5(a)所示。在不同的成像时刻，非水体区域比值在(0, 2)区间波动，这是由于在非水体区域，冰碛物的散射特性相对稳定；而对于水体区域，尽管湖水的存在造成了湖区内散射信号的差异，但水域的分布范围在各期的比值影像内较为一致。

为了精确确定比值影像在水体和非水体的分割阈值，如图 5(b)标定的位置，在冰湖L₂边界处选取9个像元，提取2014-11至2018-11的时序变化，如图 5(c)所示。像元T₃、T₆和T₉随时间变化梯度明显，而T₁、T₄和T₇几乎与时间完全失相关，冰湖边缘处T₂、T₅和T₈在2015年消融期未发现明显相关性，而在2016-2018年消融期内相对敏感，这是由于冰湖面积的不断增加，引起像元的散射面中湖水的比例逐渐增加。基于此，选择T₁、T₄、T₇的全部时序点及T₂、T₅和T₈在2016年消融期前的比值最大值，同时考虑T₃、T₆和T₉在消融期内的最小值来确定分割阈值，经计算，阈值介于(2.32, 2.59)值域内。实验最终选取2.59作为水体和非水体标准化强度比值的分割阈值，研究区域内冰湖水体和非水体的分类结果如图 6所示，黄色区域为冰湖分布范围。

4 实验结果分析 4.1 贡巴冰川末端冰湖L₁和L₂的时空变化

根据标准化强度比值的阈值分割结果获取冰湖的动态分布，再经过地理编码后，提取冰湖L₁和L₂的时序变化，获得如图 7(a)所示的冰湖蓄水面积变化曲线。由图 7(a)可知，2007-2018年期间，冰湖L₁和L₂的蓄水面积存在显著的膨胀趋势。纵观11 a间的蓄水面积变化，2018年冰湖L₁的蓄水面积较2007年增加约5倍，冰湖L₂的蓄水面积增加约2倍。考虑到2007-2011年PALSAR数据采样间隔较稀疏，仅以2015-2018年间的变化进行考量，冰湖L₁和L₂的最大蓄水面积仍呈现出明显的增长趋势；尤其是2015-2016年间，冰湖L₁和L₂的增长率分别为16%和26%；2016年消融期后至今，冰湖L₁和L₂的最大蓄水面积增长率逐渐变缓。

借助无人机影像进行分析，冰湖L₁和L₂分别承载了小贡巴冰川和大贡巴冰川较大程度的融水汇集，冰湖L₁和L₂蓄水面积的增加，不仅是贡巴冰川冰雪融水对终碛堆积体不断冲刷导致湖盆面积增加的表现，同时也是在全球气温变暖的背景下，贡巴冰川加速消融的直接反映。图 7(b)和图 7(c)分别横向对比了不同年份下冰湖L₁和L₂的蓄水和溃决变化情况。从冰湖L₁和L₂的时间轴线分析，冰湖L₁每年开始蓄水的时间发生在4月下旬至5月初，而冰湖L₂的蓄水时间则开始于每年的6月初至6月下旬，蓄水时间的差异性是由于两者水源的不同所导致。除去降雨的影响，小贡巴冰川的雪冰融水一部分从侧碛外沟谷流出，一部分经冰下暗河在4 200 m处融出后，再经冰面河道汇集至L₁；而在冰湖L₂的上游，未发现明显的河道冲刷痕迹，其水源多为冰川上游冰面融水被侧碛垄阻挡后，受重力影响向下流动，冰湖L₂的周围存在多处低洼，可以推测冰湖L₂是大贡巴冰川冰面融水的最后汇集处之一。

对比冰湖L₁和L₂的蓄水过程，冰湖L₂的蓄水速度明显高于L₁，且冰湖L₁蓄水前期的速度低于后期，这是由于在消融期初始阶段，雪冰消融较慢，进入6月后，雨热同季的气象特征给冰湖L₁和L₂的蓄水提供了充足的水源。对比L₁和L₂的生存周期，湖水水源补充的不同也导致了两者不同的生存周期。冰湖L₂于每年10月份初开始溃决，到11月中下旬完全干涸；而冰湖L₁的溃决相对缓慢，水位在每年9月中下旬开始下降，至次年的2月下旬湖水完全消失。

4.2 贡巴冰川末端冰湖L₁和L₂与气温和降雨量的关系

原中国科学院兰州冰川冻土研究所于1982-1983年在贡巴冰川进行了科学考察，详尽记录了贡嘎寺站的气温、降雨和水文状况，结合1980年以来甘孜州气象站记录的气温和降雨数据，近30多年来，该区域除气温呈上升趋势外，降雨基本保持稳定^[26-27]。贡巴冰川在12月至次年3月间，气温在零度以下，冰川几无消融，且降雨分布较少，冰湖L₂保持着湖床的干涸，冰湖L₁水位逐步降低，至2月后湖水全部消失；进入4月后，降雨量开始增加，气温回暖至零度以上，雪冰开始消融，5月小贡巴冰川融水汇集至冰湖L₁，6月高山带积雪消融，径流出现春汛，加速冰湖L₁和L₂的蓄水，同时，冰湖底部的暗河通道受温度和水压力的影响开始疏通，湖水通过湖下暗河排出，经现场实地踏勘，在冰湖L₂的湖盆底部有多处如图 8所示的地下暗河入口；在8月和9月间，冰湖L₁和L₂的补充和流出达到平衡，湖水保持动态稳定；贡巴冰川的降雨在10月以后开始逐步减少，气温也开始降低，冰湖L₁和L₂的蓄水量开始逐渐减少，到12月后，冰湖L₂完全干涸，冰湖暗河也开始冻结。从图 7(b)和图 7(c)冰湖L₁、L₂湖水面积和温度、降雨量的曲线分布可知，冰湖的蓄水和溃决与温度、降雨量的变化具有较高的相关性，以此可推测贡巴冰川的消融期为4月至10月，比青藏高原冰川的消融期(5月至9月)长两个月。

4.3 贡巴冰川末端冰湖2018年溃决异常

2018-10-17，在贡巴冰川末端实地调研时发现，冰湖L₂已完全干涸，湖床也已龟裂。在冰湖L₂西侧，一个较小湖盆被冲刷形成了一个宽47 m、深20 m、底部长度75 m、断面面积580 m²的倒三角形缺口，经量测，坝体缺口体积达16 000 m³，如图 9中的A区域所示。对比2018-10-16和2018-06-09的无人机航拍影像，此次溃决产生的涌浪引起的下游河道走移与变化显著(图 9)。考虑到冰湖L₂比往年同期蓄水更快，推测在2018-06-24至2018-07-18间有较强的降雨，强降雨作用下，水体汇集挤压松散的冰碛坝体，导致坝体溃口。

从图 7(c)2018年冰湖L₂的变化可知，冰湖L₂在2018-06-24至07-18间蓄水后，于2018-07-18至07-30间溃决，此后直到2018-11，冰湖L₂再未蓄水。通过无人机影像对比，冰湖L₂坝体在溃决前后未发生较大改变，而在湖盆底部新增了一个周长67 m、面积188 m²的地下暗河入口，如图 9中的B区域所示。与湖盆底部发现的其他暗河入口(图 8)相比较，新增暗河入口的断面面积较原入口扩大了约400倍，当湖盆底部的死冰通道渗透潜蚀贯通后，湖水在水压力的作用下，短时间内迅猛排出。

图 9所示的溃决和冰湖L₂的管涌溃决改变了大贡巴冰川末端径流的平衡，雪冰融水和降雨在贡巴冰川末端难以再汇集，尤其是冰湖L₂新增的暗河涌道在新的消融期到来时，会造成冰湖再次溃决，甚至导致洪水和泥石流发生，给下游人员和各类设施带来安全隐患。

5 结语

本文提出了基于时序SAR影像强度标准化比值的冰湖提取与动态监测方法，并选取位于贡嘎山西坡的贡巴冰川作为实验对象，基于日本航天局的ALOS/PALSAR-1与欧空局的Sentinel-1A先后获取的103景SAR影像，提取了贡巴冰川末端两个较大冰湖2007-2018年间的时空分布，并对实验对象近年来的变化过程和趋势进行了归纳和分析。

首先，本文对贡巴冰川的实验研究凸显了SAR卫星不受云雾恶劣气候条件制约的突出优势。尤其是Sentinel-1A/B卫星星座组网后，能够直接获取最快6 d重访周期的影像数据，为冰湖溃决的动态监测以及高原山地的相关遥感应用提供了稳定可靠的高现势性数据支撑。

本文提出以多时相滤波获取参考基准影像，使用时序SAR影像的强度标准化比值提取冰湖面积变化，不仅可以有效消除热噪声的影响，还能够抵消因复杂地形造成的干扰，放大了冰湖水体和非水体信号的差异，从而保证了水体分割的精度。实验结果表明，贡巴冰川末端冰湖L₁近10余年来的蓄水面积从10 000 m²增长到58 000 m²，而冰湖L₂也从32 000 m²增加了近1倍。这一现象充分反映了全球气温持续增长背景下贡巴冰川正在加速消融的现状。

通过长时序的年际变化分析，发现近10 a间贡巴冰川末端冰湖L₁和L₂具有较为稳定的生存周期。尽管冰湖一定程度上受冰川区域的温度变化及降雨量波动的影响，但冰湖L₂通常于每年10月份初开始溃决，至11月中下旬完全干涸；而冰湖L₁的溃决过程较为缓慢延续，每年9月中下旬水位开始下降，直至次年的2月下旬才干涸。然而监测数据表明，2018年消融期内，冰湖L₂底部新形成的188 m²暗河涌道已打破了这一相对稳定的生存周期。新的消融期到来后，湖水的溃决和释放可能会更加迅猛，诱发洪水和泥石流的风险激增，给下游人员和各类设施造成安全隐患，有必要重视并加强监测与防范。

以上结论可供冰湖领域研究参考，也可为该区域相关次生地质灾害的研究提供基础数据。