基于歷史邊界的喜馬拉雅山脈冰湖提取方法對比研究

陳 方， 王金曉， 張美美， 于 博

（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反射率。

由于具有相對穩定的閾值以及較好的分類性能，該水體指數已被廣泛應用于各類冰湖的提?。?2］。為便于算法處理，根據1990年的冰湖邊界生成冰湖10倍面積的緩沖區（可完全包含冰湖的所有變化），然后以該緩沖區矢量的最小外包矩形對由Landsat-8 影像生成的NDWI 圖像進行裁剪，得到待處理的568個圖像瓦片。接下來采用六種不同方法提取這些圖像瓦片上的冰湖，由于已知待提取冰湖的位置（歷史邊界），只需保留所提冰湖二值圖中位于中心位置的連通區作為最終提取的冰湖。2.2 節詳細介紹了本研究對比的六種方法。

2.2 對比方法

（1）人工閾值法。一般而言，冰湖在近紅外波段的反射率總是低于綠光波段。因此，本研究在利用人工閾值法提取冰湖時將NDWI 的閾值設置為0，即判定所有綠光波段的TOA 反射率大于近紅外波段TOA反射率的像元為冰湖像元。

（2）OTSU 閾值法。有研究利用OTSU［23］閾值法結合圖像的近紅外波段進行冰湖提?。?4］。為公平對比各方法，本研究將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）冰湖在不同分辨率影像上的位置是相對不變的，本研究所述基于歷史邊界的冰湖提取方法也適用于更高分辨率的影像?，F存的冰湖歷史邊界多來自Landsat 數據，盡管這些邊界信息不夠精細，但可作為OTSU 迭代法和C-V 迭代法的迭代起點。今后將開展基于高分辨率遙感影像的冰湖提取方法對比和優化研究，以實現中、小型冰湖邊界的精細提取。