基于历史边界的喜马拉雅山脉冰湖提取方法对比研究

陈 方， 王金晓， 张美美， 于 博

（1. 可持续发展大数据国际研究中心，北京 100094； 2. 中国科学院 空天信息创新研究院 中国科学院数字地球重点实验室，北京 100094； 3. 中国科学院大学 资源与环境学院，北京 100049）

0 引言

受气候变化和冰川消退的影响，从1990 年到2018 年，全球范围内的冰湖经历了极速的扩张和变化［1］。冰湖作为高山区域重要的水资源，也引发了许多冰川和地质灾害［2］。当地震、滑坡、雪崩等因素扰动导致湖坝崩溃时，冰湖会突然释放大量湖水，即暴发由于冰湖溃决而产生的洪水灾害［3］，给下游基础设施和人民的生命财产带来严重危害［4］。冰湖扩张一方面提高了高山区域水资源的可用性［5］，但另一方面则增大了冰湖溃决风险。在全球变暖的背景下，调查不同区域冰湖的储水量、研究冰湖的变化规律、评估冰湖溃决灾害风险对区域水资源管理、气候变化响应研究和防灾减灾等具有重要意义，而冰湖提取则是以上所有研究的重要基础。

水体在近红外波段的地表反射率显著低于蓝光和绿光波段，由蓝/绿波段和近红外波段构成的水体指数（Normalized Difference Water Index， NDWI）可以凸显遥感图像上的水体区域。因此，基于NDWI 的阈值分割是常见的冰湖提取方法［6-7］。阈值法简单快速，适用于较小的区域，当冰湖制图区域扩大时，由于不同类型、不同环境下的冰湖光谱并不相同，仅使用单一阈值分割往往造成部分区域冰湖的过高或过低估计。全域-局部迭代分割法［8］解决了这一问题，该方法根据初始分割得到的湖泊单元及其缓冲区生成自适应于该湖泊的最佳水体指数分割阈值，能更好地将冰湖从其背景中分离出来，得到了广泛的应用［9-10］。在局部分割阶段，除基于阈值的分割方法外，需要更复杂计算的活动轮廓模型也能很好地提取冰湖［11-12］。全域-局部迭代法的初始分割结果一般由简单的阈值法得到，由于冰湖和山体阴影的光谱较为相似［13］，仅对水体指数或波段比值设置阈值无法较好地区分冰湖和山体阴影。因此，该方法需要利用高程数据（Digital Elevation Model，DEM）来消除被错分的山体阴影，受高程数据的精度和相应阈值的设定影响较大。随着计算机硬件和深度学习算法的发展，近些年涌现出许多利用深度学习方法（如U-NET 网络［14］）提取冰湖的研究［15-18］，深度学习方法能自适应地从数据中学习到适用于目标任务的特征表征，但该类方法需要大量样本驱动，与常规方法相比，所需计算量极大且运行速度慢。

在实际研究中，冰湖的自动化提取过程可分为两个步骤，首先解决冰湖在哪的问题，即定位冰湖；然后解决冰湖具体形状如何的问题，即生成精准的冰湖边界。如在上述方法中，山体阴影的去除是为了让冰湖的定位更准确，而局部迭代分割是为了让冰湖的形状更精确，深度学习方法则同时解决这两个问题。近些年来，在诸多学者的努力下，越来越多的冰湖编目数据被生产和公开，如高亚洲地区1990年和2018年的冰湖编目［19］、高亚洲地区2008—2017年间每年的冰湖编目［20］。由于冰湖每年的新增数量有限，可以说这些高质量的冰湖编目在很大程度上解决了地物是不是冰湖的问题，在此基础上提取冰湖只需进行高质量的图像分割。因此，这些冰湖编目可以作为生产未来冰湖数据的基础，尤其当学者感兴趣的冰湖是某种特定类型的冰湖时，如要对曾发生过溃决灾害的冰湖进行动态监测，这些冰湖的位置通常是已知的，并且存在一定年份的相应历史边界数据。

图1展示了冰湖历史边界（位置）已知情况下人工阈值法提取冰湖的流程，尽管由于设置的NDWI阈值偏低导致图像上大面积的积雪、阴影也随冰湖一起被提出，但该图像瓦片是以待提取冰湖历史边界为中心的（详见2.1 节），只需保留提取结果中心位置的连通区，即可剔除与真实冰湖不相连的积雪、阴影等地物。这表明，冰湖的历史边界（位置）信息，极大地方便了冰湖当前边界的更新，然而，目前仍然缺乏在已知冰湖历史边界情况下的冰湖提取方法对比分析研究，而且即使已知冰湖历史边界，在冰冻、积雪、云、山体阴影等因素影响下各方法提取冰湖的适用性也尚不明确。因此，本研究对比评估了冰湖历史边界已知情况下人工阈值法、OTSU阈值法、U-NET模型和三种全域-局部迭代分割法的冰湖提取效果，为已知历史边界冰湖的长期监测提供技术支持。

图1 位置已知情况下人工阈值法的冰湖提取流程Fig. 1 Glacial lake extraction process of artificial threshold method with the known location

1 研究区概况和数据来源

本研究区域位于喜马拉雅山地区，地理位置为27°23′～28°28′ N、85°55′～87°58′ E，如图2 所示。该区域地形起伏较大，平均海拔约为4 500 m，最高海拔可达到8 771 m。区域内共有1 509 条冰川，总面积为2 078.17 km2。王欣团队2018 年的冰湖编目［19］显示该区域发育种类多样、大小不一的冰湖：与冰川连接的湖泊109个，与冰川不连接的湖泊540个，非冰川供给的湖泊296 个，冰面湖118 个；这些冰湖大小各异，最小的冰湖仅有0.0054 km2，最大的冰湖为3.8831 km2；近30年间，该地区冰湖面积迅速扩张，1990年、2018年冰湖面积分别为67.5057 km2、83.5022 km2，相对面积增长了23.7%。

图2 研究区及冰湖分布Fig. 2 Study area and the distribution of glacial lakes

Landsat 系列卫星数据具有长达几十年的地球观测记录，广泛用于冰湖监测研究。为了综合评估各方法的性能，研究采用了5景Landsat-8 OLI影像，其元数据详见表1，采用该区域1990 年的冰湖编目［19］（https：//doi. org/10.12072/casnw. 064.2019.db）作为初始冰湖边界，受气候变化影响，该编目记录的冰湖形状、面积与2015 年的冰湖有较大差异，有利于评估历史边界已知情况下各方法的冰湖提取性能。序号为1 的影像，包含488 个已知1990 年边界的冰湖，用于各方法性能的整体评估，序号为2到5 的影像，分别包含20 个已知1990 年边界的冰湖，用于评估各方法对冰冻状态、积雪覆盖、与山体阴影连接、云覆盖的冰湖的提取性能。这些冰湖（共568个）的对应矢量均由有经验的操作员在Arc-GIS 10.6软件中绘制。

表1 使用影像详细信息Table 1 Details of the images used

2 研究方法

2.1 技术路线

本研究的整体技术路线如图3 所示。首先对Landsat-8 OLI 影像进行辐射定标，即将影像DN 值转换为具有物理意义的大气顶层（Top of Atmosphere，TOA）反射率，用于生成NDWI 图像。TOA反射率影像已经可以用于提取冰湖边界，与文献［22］一致，本研究不对Landsat-8 OLI 影像进行大气校正。DN值到TOA值的转化公式如下：

图3 技术路线图Fig. 3 The technical flow chart

式中：增益Gain 和偏置Offset 可以从影像的元数据中获取。本研究采用McFeeters提出的NDWI［21］：

式中：ρGreen和ρNIR分别为绿光波段和近红外波段的TOA反射率。

由于具有相对稳定的阈值以及较好的分类性能，该水体指数已被广泛应用于各类冰湖的提取［22］。为便于算法处理，根据1990年的冰湖边界生成冰湖10倍面积的缓冲区（可完全包含冰湖的所有变化），然后以该缓冲区矢量的最小外包矩形对由Landsat-8 影像生成的NDWI 图像进行裁剪，得到待处理的568个图像瓦片。接下来采用六种不同方法提取这些图像瓦片上的冰湖，由于已知待提取冰湖的位置（历史边界），只需保留所提冰湖二值图中位于中心位置的连通区作为最终提取的冰湖。2.2 节详细介绍了本研究对比的六种方法。

2.2 对比方法

（1）人工阈值法。一般而言，冰湖在近红外波段的反射率总是低于绿光波段。因此，本研究在利用人工阈值法提取冰湖时将NDWI 的阈值设置为0，即判定所有绿光波段的TOA 反射率大于近红外波段TOA反射率的像元为冰湖像元。

（2）OTSU 阈值法。有研究利用OTSU［23］阈值法结合图像的近红外波段进行冰湖提取［24］。为公平对比各方法，本研究将OTSU 阈值法应用于NDWI 图像。OTSU 方法根据最大化类间距离且最小化类内方差的原则确定最优的灰度图像二分类阈值，可以避免阈值选择的主观性。

（3）全域-局部分割法。全域-局部的方法在全域分割阶段总是设定较低的NDWI 阈值以尽可能地分割出图像包含的所有冰湖。在得到初始冰湖单元后，可应用不同的方法对每个冰湖单元进行局部分割，下面介绍三种不同的局部分割方法。

原创的全域-局部方法［22］中，局部分割分为三步：（i）对当前冰湖单元建立等面积的缓冲区；（ii）对缓冲区内像元的NDWI 值进行直方图统计，按照双峰分布准则确定分割阈值并进行分割得到新的冰湖单元；（iii）重复（i）到（ii）步骤，直至冰湖像元数稳定。双峰分布准则的阈值确定公式为：

式中：μWater和μLand分别为当前冰湖单元和其缓冲区内像元的NDWI 均值；σWater和σLand分别为其对应的标准差。

部分研究在第（ii）步中使用OTSU 阈值分割法［24］或C-V 模型［12］进行冰湖局部优化提取。其中，OTSU 方法的原理见2.2（2）小节，C-V 模型通过变分法和梯度下降法不断演化初始边界C使得能量函数E最小化，此时的C即可将前景和背景分割。该能量函数E定义为：

式中：第一项是对边界C长度的正则化，第二项和第三项是灰度图像强度和区域均值差异的积分。f(x，y)代表图像在x，y处的灰度值；c1代表边界C内部的灰度均值；c2代表边界C外部的灰度均值。参照文献［11］，本研究将μ设置为0.1，λ1、λ2均为1。

在应用全域-局部法提取冰湖时，上述两项研究在利用OTSU 方法和C-V 模型分割出冰湖后，并未进行第（iii）步中的迭代过程。当初始冰湖单元的等面积缓冲区完全包含真实冰湖的所有像元时，迭代过程不是必需的。而本研究中冰湖的初始边界（1990 年）和冰湖真实边界（2014 年以后）可能相差很大，若不进行第（iii）步中的迭代分割，则可能最终只提取出真实冰湖的小部分区域，因此，本研究增加了第（iii）步的迭代过程。在使用全域-局部法时，本文利用1990年冰湖编目数据作为初始冰湖单元，这使得全域分割不再必要，且局部分割都包含迭代分割的过程，因此将本节介绍的三种方法命名为：双峰迭代法、OTSU 迭代法和C-V 迭代法。将迭代分割的终止条件设置为：相邻两次分割所得冰湖的相对面积差小于等于1%。

（4）深度学习方法U-NET。U-NET 是一种对称的编码解码语义分割网络，其核心力量来源于跳连接，结构详见图4。U-NET 在编码阶段通过堆叠卷积层和池化层提取输入图像的多层级特征，池化层以保留特征图中各子区域内最强特征的方式减小特征图尺寸，同时使特征图丢失了很多细节信息，U-NET 在解码阶段通过不断堆叠上采样层和卷积层逐步得到和输入图像长宽一致的预测图，跳连接使得上采样后的粗分辨率高级语义特征和编码阶段的高分辨率浅层细节特征相融合，从而恢复因池化损失的目标细节信息，实现地物的精细分割。UNET 这类基于学习的模型总是需要一定规模的样本驱动，本研究将序号为1 的影像分为A、B 两个区域（图2），分别包含237、251 个冰湖，使用A 区域的237 个图像瓦片训练U-NET 并对B 区域的251 个冰湖进行预测，再使用B 区域的251 个图像瓦片训练U-NET 并对A 区域的237 个冰湖进行预测，由此得到序号为1 的影像上488 个冰湖的预测结果。为了提取冰冻状态、积雪覆盖、云覆盖、与山体阴影连接的共80 个冰湖，使用A 区域和B 区域的488 个图像瓦片训练U-NET。由于本研究中用于训练的图像瓦片数量较少，且为NDWI 单波段图像，研究将UNET 编码阶段滤波器（卷积核）的数量设置为32、64、128、256、512（图4）。在训练U-NET 时采用SGD（Stochastic Gradient Descent）作为优化器，动量设置为0.9，由于冰湖和背景地物的面积占比极不均衡，使用Dice 损失作为网络的损失函数，Batch 大小设置为4，初始学习率设置为0.0025，采用Poly策略逐渐降低学习率，共训练50轮。考虑到样本数量较少，在训练过程中采用随机水平翻转和竖直翻转进行样本增强。图5 展示了利用A 区域和B 区域图像瓦片训练U-Net的损失曲线。

图4 U-NET结构图Fig. 4 The Structure of U-NET

图5 U-Net的训练损失曲线Fig. 5 Training loss curve for U-Net： the training process of image tiles in region A （a），the training process of image tiles in region B （b）

3 结果分析

3.1 成像条件较好情况下各方法的冰湖提取结果

冰湖提取结果的精度评价采用精确率（Precision）、召回率（Recall）和F1 分数三个指标，其定义参考文献［25］。表2 呈现了2.2 节中所述六种方法对序号为1 的Landsat 影像上488 个冰湖的提取精度。人工阈值法提取的冰湖常包含与真实冰湖相连的冰、雪［图6（c）第3 行］、阴影和部分陆地［图6（c）第2行］等NDWI值大于0的地物，因此该方法的精确率只有66.02%，且F1 分数最低。OTSU 阈值法利用了图像瓦片全局的NDWI 统计信息，当真实冰湖和图像其他区域的NDWI 值有明显差异时，该方法可以很好地分割出冰湖［图6（d）第2 行和第4行］。然而，当图像可按NDWI 值从高到低分为三类（冰湖、阴影/云/雪、陆地）时，该方法提取的冰湖也会包含与真实冰湖连接的阴影、雪［图6（d）第1行和第3 行］。因此，OTSU 阈值法的精确率也比较低，只有75.76%。值得注意的是，这两种方法的召回率都很高，分别为96.70%和94.25%，即能够尽可能地提取出所有冰湖。

表2 六种方法的冰湖提取精度对比Table 2 Accuracies of the six glacial lake extraction methods

图6 六种方法的冰湖提取结果示例Fig. 6 Examples of glacial lake extraction results by six methods： true color images （a）， NDWI （b）， results of artificial threshold method （c）， results of OTSU threshold method （d）， results of bimodal iterative method （e）， results of OTSU iterative method （f）， results of C-V iterative method （g）， predictions of U-NET （h） and ground truths （i）

以1990 年的冰湖边界作为初始边界的OTSU迭代法和C-V 迭代法取得了非常高的整体精度，F1分数分别达到88.89%和89.30%，这是因为这两种方法在当前冰湖及其等面积的缓冲区内迭代地提取冰湖。冰湖和其等面积缓冲区内像元的NDWI值往往呈现双峰分布，与基于全局NDWI 的OTSU阈值法相比，OTSU 迭代法能给出更恰当的分割阈值。对于C-V 迭代法而言，等面积缓冲区使得式（4）中第二、三项在数量上更均等，更利于边界的演化。同样以1990 年冰湖作为初始冰湖单元的双峰迭代法却无法完整地提取出冰湖［图6（e）］，该方法的F1 分数只有77.67%，原因是双峰分布的阈值确定准则式（3）不适用于初始冰湖边界（1990 年）和真实冰湖边界（2015年）相差很大的情况，详见4.1节。U-NET 方法的冰湖提取精度最高，其F1 分数为89.8%，但对于面积>0.1 km2的冰湖，其最终提取的冰湖与真实冰湖相比总是存在缺失［图6（h）第1 行和第2行］，这可能是由小冰湖在训练样本中占比较高（图7）导致。

图7 冰湖面积和提取精度间的关系Fig. 7 The relationship between glacial lake areas and extraction accuracies of OTSU iterative method，C-V iterative method and U-NET

为了进一步评估各方法的冰湖提取性能，按照面积将上述488个冰湖分为三类：小型冰湖（≤0.01 km2），中型冰湖（0.01～0.1 km2），大型冰湖（>0.1 km2）。由表3可知，受影像空间分辨率（30 m）的影响，各方法对大型冰湖的提取精度（F1 分数）往往大于对中型冰湖的提取精度，对小型冰湖的提取精度则最低。对三类不同大小的冰湖，OTSU 迭代法、C-V 迭代法和U-NET 取得的F1 分数明显高于人工阈值法、OTSU 阈值法和双峰迭代法，这与表2 呈现的特点一致。具体地，OTSU 迭代法、C-V 迭代法和UNET 在小型冰湖、中型冰湖和大型冰湖上取得的F1分数分别集中在82%、88%和94%附近，对小型冰湖而言，OTSU 迭代法取得的F1 分数最高（82.79%），对中型冰湖而言，U-NET 取得的F1 分数最高（89.69%），对大型冰湖而言，C-V 迭代法取得的F1 分数最高（94.91%），整体上，这三种方法都适用于对各类不同大小冰湖的提取。图7更详细地展示了OTSU 迭代法、C-V迭代法、U-NET对不同大小冰湖的提取精度，尽管这三种方法取得的F1分数都在89%左右，仍然有一些小面积冰湖的F1分数低于0.6，原因是小冰湖无法提供足够的统计信息。

表3 六种方法对不同大小冰湖的提取精度Table 3 Extraction accuracy of six methods for different sizes of glacial lakes

3.2 复杂成像条件下各方法的冰湖提取结果

表4～7 展示了本文所述六种方法对冰冻状态、积雪覆盖、云覆盖、与山体阴影连接的冰湖的提取精度。尽管冰冻状态冰湖的NDWI 值有所降低，与周围背景地物的NDWI 值仍有较明显差异［图8（b）］，OTSU 迭代法、C-V 迭代法和U-NET 都能很好地提取冰冻状态的冰湖，且与表2所示结果一致，这三种方法的冰湖提取精度明显优于人工阈值法、OTSU 阈值法和双峰迭代法，其中C-V 迭代法取得了最高的F1 分数（90.61%）。与双峰迭代法和OT-SU 迭代法这类基于阈值分割迭代提取冰湖的方法相比，C-V 迭代法不受冰湖内部NDWI 值较低的部分像素影响，能完整地提出冰湖（图8 第三行）。不同于冰冻状态的冰湖，积雪覆盖下冰湖的NDWI 值极低，与周围背景的NDWI 值十分接近［图9（b）］，因此，人工阈值法、OTSU 阈值法、双峰迭代法、OTSU 迭代法和U-NET 的冰湖提取结果都较差，其中U-NET 取得的F1 分数仅有13.57%，只能将冰湖中未被积雪覆盖的一小部分提出［图9（i）］。虽然积雪覆盖下冰湖的NDWI 值极低，但冰湖边界仍然较为清晰，故C-V 迭代法取得了比其他方法显著更高的F1分数（81.40%）。

表4 六种方法对冰冻状态冰湖的提取精度Table 4 Extraction accuracy of six methods for frozen glacial lakes

表5 六种方法对积雪覆盖冰湖的提取精度Table 5 Extraction accuracy of six methods for glacial lakes covered with snow

表6 六种方法对云覆盖冰湖的提取精度Table 6 Extraction accuracy of six methods for glacial lakes covered by clouds

表7 六种方法对与山体阴影连接冰湖的提取精度Table 7 Extraction accuracy of six methods for glacial lakes connected to the mountain shadows

图8 冰冻状态冰湖提取结果示例Fig. 8 Extraction results of frozen glacial lakes： true color images （a）， NDWI （b）， lake contours in 1990 （c）， results of artificial threshold method （d）， results of OTSU threshold method （e）， results of bimodal iterative method （f）， results of OTSU iterative method （g）， results of C-V iterative method （h）， predictions of U-NET （i） and ground truths （j）

图9 积雪覆盖冰湖提取结果示例Fig. 9 Extraction results of glacial lakes covered by snow： true color images （a）， NDWI （b）， lake contours in 1990 （c）， results of artificial threshold method （d）， results of OTSU threshold method （e）， results of bimodal iterative method （f）， results of OTSU iterative method （g）， results of C-V iterative method （h）， predictions of U-NET （i） and ground truths （j）

云覆盖下冰湖的NDWI值与周围背景的NDWI值仍有较明显的差异［图10（b）］，因此，OTSU 迭代法和C-V 迭代法都能较好地提出冰湖，其中C-V 迭代法取得了最高的F1 分数（89.51%），U-NET 则几乎无法提出云覆盖下的冰湖［图10（i）］，取得的F1分数仅有2.26%。当冰湖与山体阴影相连接时，由于两者的NDWI 值十分接近，人工阈值法和OTSU阈值法都会误提大量山体阴影［图11（d）和（e）］，即使OTSU 迭代法和C-V迭代法也无法较好地分离二者［图11（g）和（h）］，仅U-NET这种基于学习的方法能较完整地分割出冰湖，取得的F1 分数高达85.19%，显著优于其他方法。

图10 云覆盖冰湖提取结果示例Fig. 10 Extraction results of glacial lakes covered by clouds： true color images （a）， NDWI （b）， lake contours in 1990 （c），results of artificial threshold method （d）， results of OTSU threshold method （e）， results of bimodal iterative method （f），results of OTSU iterative method （g）， results of C-V iterative method （h）， predictions of U-NET （i） and ground truths （j）

4 讨论

4.1 双峰迭代法与OTSU迭代法的对比

双峰迭代法和OTSU 迭代法都以1990 年冰湖作为初始冰湖单元，且都基于阈值对缓冲区内的像元进行分割，二者唯一的区别在于分割阈值的计算公式不同。本节以图12展示的冰湖为例，探究两种方法表现差异的原因。

图12 两种迭代法的冰湖提取结果示例Fig. 12 Examples of glacial lake extraction results by two iterative methods： true color image （a）， NDWI image （b）， the glacial lake of 1990s （c）， result of bimodal iterative method （d）， result of OTSU iterative method （e） and ground truth （f）

图13 和图14 展示了双峰迭代法和OTSU 迭代法的冰湖迭代提取过程，在以1990年冰湖为初始冰湖单元进行第一次迭代分割时，两种方法要处理的冰湖等面积缓冲区的NDWI 值完全一致，但基于双峰分布的阈值确定准则式（3）给出的NDWI 分割阈值（0.551）明显高于OTSU 的阈值（0.357）。这是因为1990 年冰湖只是2015 年冰湖的一小部分，由此得到的等面积缓冲区除包含陆地像元外还有大量的水体像元，导致式（3）计算的阈值更偏向直方图中的水体（图13中迭代次数i=0的直方图）。较高的阈值会导致冰湖漏分。因此，本次分割得到的新冰湖及其等面积缓冲区会包含更高比例的水体像元，并再次由式（3）得到更高的分割阈值（0.647），导致冰湖进一步漏分。

图13 双峰迭代法的冰湖迭代提取过程，黑白二值图代表第i次迭代分割后的冰湖，与其对应的直方图显示该冰湖等面积缓冲区内NDWI的分布，红色竖线代表分割阈值Fig. 13 The iterative process for glacial lake extraction of bimodal iterative method， binary images represent segmentation results after the i-th iteration， histograms show the corresponding distributions of NDWI in equal-area buffers of glacial lakes and segmentation thresholds displayed in red lines

图14 OTSU迭代法的冰湖迭代提取过程，黑白二值图代表第i次迭代分割后的冰湖，与其对应的直方图显示该冰湖等面积缓冲区内NDWI的分布，红色竖线代表分割阈值Fig. 14 The iterative processes for glacial lake extraction of OTSU iterative method， binary images represent segmentation results after the i-th iteration， histograms show the corresponding distributions of NDWI in equal-area buffers of glacial lakes and segmentation thresholds displayed in red lines

表8 展示了双峰迭代法对图12 所示冰湖进行提取时，每次迭代前冰湖和缓冲区的像元数、冰湖和缓冲区的NDWI 均值和标准差以及迭代后冰湖的像元数。如图13 所示，双峰迭代法从第1 次到第5 次迭代给出的分割阈值从0.551 增加至0.686，这期间分割所得的冰湖范围不断增大；从第6 次到第10 次迭代，分割阈值稳定在0.69 左右，并且分割所得冰湖面积不断减小，直至达到稳定的漏分状态（图13 迭代次数i=11）。在此过程中，冰湖缓冲区内的NDWI 直方图从不太显著的双峰分布逐渐演化为只包含水体像元的单峰分布。与双峰分布的阈值确定准则相比，OTSU 在缓冲区内水体像元较多的情况下仍能给出恰当的分割阈值，这是由于OTSU 准则在计算二分类阈值过程中不涉及由当前冰湖或缓冲区计算的参量。恰当的分割阈值使得分割后冰湖的等面积缓冲区包含更多的陆地像元，因此OTSU 迭代法给出的分割阈值从第1 次迭代的0.357降至第5次迭代的0.3（图14），且直方图呈现的双峰分布态势在迭代过程中逐步明显。

表8 双峰迭代法的过程参数Table 8 The process parameters of the bimodal iterative method

4.2 C-V模型的参数设置

在上述三种迭代法中，C-V 迭代法取得的F1 分数最高（89.30%），与OTSU 迭代法相比仅提升了0.40%，但在每次迭代过程中C-V 模型都会对初始边界进行多次演化，因此C-V 迭代法比OTSU 迭代法需要的时间成本更高。OTSU 方法确定分割阈值的准则是无参数的，而C-V 模型的演化则需要设定三个参数，即式（4）中的μ、λ1、λ2。通常λ1、λ2都设置为1，μ则根据待分割对象的几何特征设定。待分割对象边界越圆滑则μ越大，待分割对象边界越不规则（周长相对较长）或面积较小，则μ越小。本研究中待分割的488 个冰湖大小各异，μ值的设置直接影响C-V 迭代法的分割结果（如表9 所示）。μ为0.05 时，C-V 迭代法取得的F1 分数最高（89.62%），随着μ值的增大，C-V 迭代法的整体精度不断下降。

表9 μ值不同时C-V迭代法的冰湖提取精度Table 9 Extraction accuracies of glacial lakes by C-V iterative method with different μ values

图15 展示了μ值不同时C-V 迭代法的冰湖提取结果。对于小冰湖（第一行），仅当μ=0.05 时C-V迭代法才能完整地提取冰湖；对于大冰湖（第二行），不论μ取何值，C-V 迭代法都能将冰湖完整地提出，但随着μ值的增大，提取的冰湖边界越来越光滑。如表10所示，对于小型冰湖，μ=0.05时的F1分数比μ=0.20 时高6.86%；对于中型冰湖，μ=0.05 时的F1分数比μ=0.20 时高2.67%；当冰湖面积>0.1 km2时，μ=0.05 时的F1 分数比μ=0.20 时低0.21%。由此可知，μ的取值对面积越小的冰湖影响越大，小冰湖应设置较低的μ值，大冰湖则对μ的取值不敏感。因此，建立冰湖面积与μ值之间的关系，在对不同大小的冰湖进行提取时设置相应的最佳μ值，将进一步提升C-V迭代法的冰湖提取精度。

表10 μ=0.05和μ=0.20时C-V迭代法提取不同大小冰湖的精度统计Table 10 Extraction accuracies of glacial lakes in different sizes by C-V iterative method with different μ values

图15 不同μ值的C-V迭代法的冰湖提取结果Fig. 15 Extraction results of glacial lakes by C-V iterative method with different μ values： true color image （a），NDWI image （b）， μ=0.05 （c）， μ=0.10 （d）， μ=0.15 （e）， μ=0.20 （f）

4.3 U-NET的适用性

在成像条件较好的488个图像瓦片上训练的UNET能较好地提取冰冻状态的冰湖和与山体阴影连接的冰湖，其取得的F1 分数分别为88.20% 和85.19%，但对于积雪覆盖和云覆盖冰湖的提取，UNET表现十分糟糕，与其他方法相比其取得的F1分数最低。这可能是由用于训练和测试的冰湖的NDWI值分布差异较大导致，经统计，冰冻状态和与山体阴影连接冰湖的NDWI 均值分别为0.4912 和0.7908，与序号为1的图像中的488个冰湖的NDWI均值0.6175较为接近，而积雪覆盖和云覆盖下冰湖的NDWI 均值分别为0.0841 和-0.0003，与0.6175差距极大。因此，若想使得U-NET 这类深度学习方法在各种场景下都取得优异的冰湖提取效果，构建包含各种类型、不同成像条件的冰湖样本集至关重要。

5 结论

本研究利用喜马拉雅地区2014年后的Landsat-8 OLI 影像，基于1990 年冰湖编目数据提供的历史边界（位置）信息，对成像条件较好的488 个冰湖和受积雪、冰冻、云、山体阴影影响的80个冰湖开展了冰湖提取方法对比实验，结论如下：

（1）U-NET、OTSU 迭代法、C-V 迭代法都能精确地提取成像条件较好的冰湖，F1 分数均达到88.80%以上，显著优于人工阈值法和OTSU 阈值法。由于1990 年与2015 年冰湖范围差异较大，利用冰湖及其缓冲区内统计信息的双峰迭代法易造成过分割（漏分）现象。

（2）在提取冰冻状态、积雪覆盖和云覆盖下的冰湖时，应优先采用C-V 迭代法。而U-NET 受样本分布影响较大，尽管无法有效提取积雪覆盖和云覆盖下的冰湖，但相比其他方法，能够更有效地分离冰湖和与之连接的山体阴影。

（3）在提取大型、中型和小型冰湖时，优先采用的方法分别为C-V 迭代法、U-NET 和OTSU 迭代法。由于本研究的实验数据为中分辨率Landsat 8影像，六种方法对大型冰湖的提取精度显著高于中、小型冰湖，可采用更高空间分辨率的遥感影像，如Sentinel-2（10 m）或高分2 号数据（4 m）精细刻画中、小型冰湖的边界。

（4）冰湖在不同分辨率影像上的位置是相对不变的，本研究所述基于历史边界的冰湖提取方法也适用于更高分辨率的影像。现存的冰湖历史边界多来自Landsat 数据，尽管这些边界信息不够精细，但可作为OTSU 迭代法和C-V 迭代法的迭代起点。今后将开展基于高分辨率遥感影像的冰湖提取方法对比和优化研究，以实现中、小型冰湖边界的精细提取。