基于NDWI-NDSI组合阈值法的布加岗日冰湖提取及其变化分析

闫斌，贾洪果，任文静，吴仁哲，黄心茹

西南交通大学 地球科学与环境工程学院,成都 611756

1 引 言

冰川是重要的水资源，由于其对河川径流起着自然地调节作用，被誉为“高山固体水库”；冰川的运动消融也极易引发自然灾害，如冰川跃动、冰湖溃决、冰川泥石流等，因此冰川研究对人类社会具有重要意义（邬光剑等，2019；张九天等，2012；陈宁生等，2019）。冰川作为一种动态资源，与气候变化密切相关（刘潮海等，2002；赵晓艳等，2022）。相关研究（陈虹举等，2017；安国英等，2019）表明：中国冰川变化趋势主要呈现对夏季平均气温变化表现为正响应，而对年降水量变化表现为负响应，尽管冰川分布区年降水量的增加导致冰川积累量增多，但这并不足以抵消因温度升高而增加的消融量，因此温度的升高现已成为中国西部冰川快速退缩的主导性因素。

冰湖是冰川运动的产物，通常由冰川挖蚀成的洼坑和冰碛物堵塞冰川槽谷积水而形成。冰湖个数增多、面积增大的同时，冰川则往往处于退缩状态，因此冰湖的形成演变与冰川的变化密切相关（姚晓军等，2017）。在全球气候变暖的背景下，由冰湖溃决引发的洪水和泥石流灾害呈现数量增多、危害程度加剧特征（Wang 和Zhou，2017）。由于冰湖多位于高海拔地区，其溃决所产生的地质灾害具有突发性强、破坏性大、持续时间短、波及范围广的特点（Bolch 等，2008；韩金良等，2007；刘建康等，2019；杨瑞敏等，2019）。据《昌都地区志》记载，1974 年7 月6日，丁青县的波戈冰川冰湖溃决，大量湖水泄入河谷，沿河木桥全部冲走，黑（河镇）昌（都）公路部分受损（姚晓军等，2014）。因此对冰湖溃决灾害研究具有重要的社会意义。

本文所选研究区为位于藏东南地区唐古拉山东段的布加岗日冰川，该冰川的冰湖周围地物主要包括冰舌和陆地。目前，在冰湖区基于归一化水指数NDWI（Normalized Difference Water Index）进行冰湖提取已有先例借鉴（李均力等，2011；Li 和Sheng，2012；骆剑承 等，2009；Watson 等，2018；Wendleder 等，2018），此种方法适用于目标冰湖与周围地物信息在NDWI的直方图中呈双峰特性的情况；并实验发现，若冰湖与周围地物信息在NDWI直方图上呈现三峰，如果直接使用阈值分割法极易造成错分现象。针对上述问题，本文提出基于NDWI−NDSI 组合阈值的冰湖面积提取方法。首先利用归一化雪被指数NDSI（Normalized Difference Snow Index）具有分离陆地（目标）与冰水混合地物（背景）的特性制作去陆地掩膜，再对掩膜后的NDWI进行阈值分割，以准确提取冰湖面积信息。为验证方法的可行性，首先基于冰川区域分布的主要地物的NDWI−NDSI 散点图和冰湖区的NDSI 直方图进行陆地与冰水混合地物分离，接着分别基于NDWI 和NDWI−NDSI 两种方法对冰湖提取进行对比分析。

2 研究区概况及数据介绍

如图1所示，布加岗日冰川位于距离川藏北线G317 国道不远的西藏自治区那曲市索县、巴青县及昌都市丁青县的交界处，与唐古拉山脉一致为东西走向。其海拔较高，有16 座海拔超过6000 m的高峰，最高峰海拔为6328 m。受来自印度洋的西南季风作用，年降水量达到600—700 mm，高耸的地势和丰富的降水使这里成为较大的山岳冰川发育集中地（Liu等，2016；王宁练和丁良福，2002；王聪强等，2016）。据光学影像目视解译，该冰川发育面积较为明显的冰湖约为12 个，它们分别对应编号L−1，L−2，…，L−12，其位置如图2所示。

图1 布加岗日冰川地理位置及主要冰川分布Fig.1 Geographical location of Bujiagangri glacier and distribution of major glaciers

图2 布加岗日冰川冰湖分布图（2016−08−19）Fig.2 Distribution map of Bujiagangri glacial and glacial lakes（2016−08−19）

本文所用数据源为美国NASA 的陆地卫星（Landsat）TM/OLI，卫星波段介绍如表1 所示。此次实验所使用的Landsat5/8 卫星波段有：Green 波段、NIR波段、SWIR1波段，所选多光谱波段分辨率均为30 m。基于光学影像成功提取冰湖面积的前提条件是影像获取时间段内研究区上空无云或云雾较少、且冰湖处于非结冰期，而布加岗日所在地区的夏季多数影像数据被云雾覆盖严重无法使用，同时部分冬季影像由于冰湖结冰或被雪覆盖也导致数据不符合要求。考虑以上因素，本文筛选出1988 年—2018 年间共16 景质量较好且覆盖布加岗日冰川的Landsat（TM/OLI）影像。对所选影像首先裁剪处理，并进行辐射定标和大气校正，以消除或减弱由传感器和大气带来的各种误差，避免其对后期处理结果的精确性产生影响。

表1 Landsat 5/8卫星波段介绍Table 1 Introduction of Landsat 5/8 satellite bands

3 研究方法

NDWI（McFeeters，1996）又称为归一化水指数，是基于水体的反射从可见光到中红外波段逐渐减弱的原理，利用水体的高反射率主要集中在可见光中的蓝绿光波段，而在可见光其它波段及近、中红外的反射率很低的特性，从而使水体的亮度信息得到最大程度增强，非水体受到普遍的抑制，进而达到突出水体的目的（Benoudjit 和Guida，2019；徐涵秋，2005），其表达式为

由Hall 等（1995） 提出的归一化雪被指数NDSI 是目前光学遥感提取积雪的通用方法。实验发现，该指数除了能突出积雪信息外，对水体和冰同样具有高敏感性，同时又对陆地信息具有很强的抑制作用，本文正是利用这一特性实现对冰湖区陆地信息的识别。如图3所示，结合冰湖区影像的NDWI−NDSI 散点图中各地物信息表现出的差异性与其NDSI 灰度直方图呈现双峰的条件，可由OTSU 法（Otsu，1979）计算得到分割阈值（红线所示），据此可以分离出陆地信息。计算公式为

图3 基于NDSI分离陆地与冰水的方法验证Fig.3 Method verification of separation of land，ice and water based on NDSI

基于NDWI−NDSI 组合法实施冰湖提取的过程如图4 所示，与通常基于NDWI 提取冰湖不同的是，本文首先利用OTSU 对NDSI 进行阈值分割得到阈值T1，将NDSI 中小于T1 的部分用于制作掩膜以剔除陆地信息，避免了NDWI在灰度直方图呈现三峰的可能，再利用OTSU 对掩膜后的NDWI 所包含的地物信息（积雪、冰舌和冰湖）直方图进行分割得到阈值T2。

图4 冰湖提取流程Fig.4 Flow chart of glacial lake extraction

由于局部影像中积雪量较少，不会在NDWI直方图中产生三峰现象；其次，积雪与冰舌在可见光的蓝绿波段反射率较低，而水体则相反，在NDWI 直方图中能很好地利用阈值分割法进行提取。因此，可继续使用阈值分割法对阈值小于T2 的积雪和冰舌信息进行剔除，并统计大于T2的像素个数，最后结合30 m 的影像分辨率成功实现冰湖面积的提取。为了验证该方法的正确性，本文分别选取冰湖L−10 和L−12 用于分析（图5）。L−10 与冰舌和陆地相连，其NDWI 的直方图呈现三峰，直接基于NDWI 提取的冰湖与RGB 影像对比可知，提取结果出现非冰湖区域；而基于NDWI−NDSI 组合法提取的冰湖则与目视判断的基本一致，表明了该方法的正确性。与L−10 不同的是，L−12 周围地物只有陆地，其NDWI 直方图显示双峰，分别基于NDWI 和NDWI−NDSI 两种方法实施冰湖提取结果基本相同，说明该组合法仍适用于双峰情况；此外，将两种方法提取的冰湖面积与目视解译结果进行对比可知（表2），NDWI−NDSI 组合法在上述两种情况下提取冰湖的精度明显优于NDWI法。

表2 冰湖提取方法精度评价Table 2 Accuracy evaluation of glacial lake extraction method

图5 基于NDWI与NDWI−NDSI的冰湖提取方法对比验证Fig.5 Comparative verification of two glacial lake extraction methods based on NDWI and NDWI−NDSI

4 结果分析

本文利用NDWI−NDSI 组合法成功提取了1988 年—2018 年的布加岗日冰川的冰湖面积，并从单个和整体两方面做了时序研究和统计分析。

4.1 单个冰湖变化分析

由图6 可知，30 年来布加岗日冰川的多数冰湖都有不同程度变化。除L−07、L−09、L−12这3个冰湖面积变化不大外，L−01、L−03、L−11 面积有一定增加，L−02、L−04、L−05、L−06、L−08、L−10面积增加趋势明显，其中L−04和L−06面积从无到有，这说明该冰川冰湖个数在增加、面积在扩大，也表明其融水逐年增多，融化加速。结合图2冰川冰湖分布可知，变化不大的3个冰湖均未与冰舌相连，但面积增加明显的6 个冰湖中有5 个与冰舌相连，而冰舌前端亦即冰川末端，其融水用于补充冰湖，这说明冰川消融是引起冰湖面积变化的主要因素。此外，有必要指出的是冰湖面积出现有→无→有的变化是由于冰湖结冰引起的，如L−01 在1990 年、1992 年、1993 年、2002 年、2018 年均处于结冰状态。

图6 冰湖面积年际变化Fig.6 Interannual change of glacial lakes area

为了进一步分析冰湖变化，对该冰川的12 个冰湖均做面积分级及其变化量分级处理。首先将其逐一按30 年中扩张至最大年份的面积进行分级得到表3，其次以1988 年冰湖面积为基准，用面积最大年份的冰湖面积最减去1988年的冰湖面积，可得到近30 年中每个冰湖面积的最大变化量，再将变化量进行分级得到表4，最后将12 个冰湖的面积及其变化在空间的分布情况进行了分级展示，见图7。由表3 可知，面积超过0.5 km2的冰湖有6个，面积在0.1—0.5 km2之间的冰湖有3个，其余均小于0.1 km2；如表4所示，面积变化最大的冰湖分别为L−04、L−05、L−10，三者变化量均超过0.5 km2；其次为L−02、L−06、L−08，变化量在0.1—0.5 km2之间；L−07、L−09、L−12 这3 个冰湖面积变化均小于0.01 km2。另外对冰湖面积及其变化量取交集可得，面积最大且变化量最大的3个冰湖分别为L−04、L−05、L−10，这3 个冰湖极易发生溃堤且由此引发的地质灾害会更为严重。

图7 冰湖面积及其变化空间分布Fig 7 Spatial distributions of the area of glacial lakes and their changes

表3 冰湖面积分级表Table 3 Classification of glacial lakes area

表4 冰湖面积变化分级表Table 4 Classification of glacial lakes area changes

4.2 整体面积变化分析

在4.1 节中针对单个冰湖面积的时序演变情况做了分析，本节将对12 个冰湖面积累加以分析整个冰川的冰湖变化情况。但由于影像数据受结冰、落雪及云雾遮挡等影响导致部分较大冰湖面积未被正确提取，因此从16 景影像中筛选出1988 年、1994 年、1998 年、2013 年、2016 年 这5 组质 量较好（即冰湖不受上述几类影响因素过多干扰而被正确提取）的冰湖数据进行分析。

如图8（a）与图8（b）所示，布加岗日冰川冰湖总面积在1988年约为2.7666 km²，1994年增至约2.9799 km²，1988年—1994年期间年平均变化速率约为0.0356 km²/a；1998年增至约3.1959 km²，1994年—1998 年期间年平均变化速率约为0.0540 km²/a；2013 年增至约5.0409 km²，1998 年—2013 年期间年平均变化速率约为0.1230 km²/a；2016 年增至约5.2308 km²，2013年—2016年期间年平均变化速率约为0.0633 km²/a。可以看出近30 年布加岗日冰川的冰湖面积变化呈慢→较快→很快→较快趋势，虽然面积变化不能准确反映冰川的消融量及湖水变化量，但其结果仍然可以表明布加岗日冰川退缩明显。结合图9可知，冰湖总面积变化与温度呈现正相关特性（该地区温度数据为中国气象局提供），表明该地区逐年气温的不断升高正是冰川消融过快的主要原因；其次，经统计发现，布加岗日冰川地区近30 年来平均气温在0 ℃以上的年份共7 个，而1998 年—2013 年（红色虚线框内）冰湖面积增长速率之所以最快，与这段时间内年均气温在0 ℃以上的年数较多紧密相关（5 个年份）；2013 年—2016 年（橙色虚线框内）增长速率变化次之，其原因为年均气温在0℃以上的年数相应较少（2个年份）。

图8 冰湖总面积年际变化及其变化速率Fig 8 Interannual change and rate of change of total area of glacial lakes

图9 布加岗日地区年均气温变化Fig 9 Annual average temperature change in Bujiagangri area

5 结 论

本文利用改进的NDWI−NDSI 的组合阈值法，对1988 年—2018 年研究区内的16 景Landsat 影像进行了冰湖提取试验，研究的主要结论如下：

（1）归一化积雪指数（NDSI）除了能够提取积雪外还可以用于分离试验区内陆地信息和冰雪水混合信息，据此可制作陆地掩膜；

（2）基于NDWI的阈值分割法在冰湖区提取含有冰舌、陆地两种地物信息的冰湖时会出现错分现象，改进的NDWI−NDSI 的组合阈值法成功解决了这一问题，增加了阈值分割法提取冰湖的适用范围；

（3）通过对试验区提取的冰湖做时序研究和统计分析可知，30 年来布加岗日冰川的冰湖面积从约2.7666 km²增加至约5.2308 km²，增加最快时的年际变化率约0.1230 km²/a。12 个冰湖中多数冰湖面积都有不同程度增加，其中L−04、L−05 和L−10 这3 个冰湖的面积（0.5—1.0 km²）及其扩张面积（0.5—1.0 km²）最大，这3 个冰湖受夏季受高温天气影响蓄水过快随时有溃堤的可能性，对下游居民的人身财产安全及川藏公路（G317）等形成潜在威胁，希望有关部门给予重视。此外，基于单个和整体的冰湖面积变化分析均表明该冰川30年消融过快，退缩严重。

本文的不足之处为提出的NDWI−NDSI 组合阈值法只能适用于单个冰湖提取，且该方法准确性受冰湖周围地物的复杂性影响，是否能在其它冰川用于冰湖提取有待验证；另外值得说明的是，由于受实验区域在该时间段内所获取的光学影像受云雾遮挡较为严重所限，很难于每年相同或相近月份收集到符合要求的数据。因此，只能通过比较年际带有一定时间跨度的冰湖面积进行后续分析。在今后的研究中若能获取到时间跨度较小的高质量的影像数据，将会考虑如季节性问题等对冰湖提取所带来的影响，从而进行更加精确细致的结果分析。

志 谢此次所用光学数据来自美国NASA的陆地卫星（Landsat 5/8），温度数据为中国气象局提供，在此表示衷心的感谢！