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

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

（1. 可持續(xù)發(fā)展大數(shù)據(jù)國際研究中心，北京 100094； 2. 中國科學(xué)院 空天信息創(chuàng)新研究院 中國科學(xué)院數(shù)字地球重點實驗室，北京 100094； 3. 中國科學(xué)院大學(xué) 資源與環(huán)境學(xué)院，北京 100049）

0 引言

受氣候變化和冰川消退的影響，從1990 年到2018 年，全球范圍內(nèi)的冰湖經(jīng)歷了極速的擴(kuò)張和變化［1］。冰湖作為高山區(qū)域重要的水資源，也引發(fā)了許多冰川和地質(zhì)災(zāi)害［2］。當(dāng)?shù)卣稹⒒?、雪崩等因素擾動導(dǎo)致湖壩崩潰時，冰湖會突然釋放大量湖水，即暴發(fā)由于冰湖潰決而產(chǎn)生的洪水災(zāi)害［3］，給下游基礎(chǔ)設(shè)施和人民的生命財產(chǎn)帶來嚴(yán)重危害［4］。冰湖擴(kuò)張一方面提高了高山區(qū)域水資源的可用性［5］，但另一方面則增大了冰湖潰決風(fēng)險。在全球變暖的背景下，調(diào)查不同區(qū)域冰湖的儲水量、研究冰湖的變化規(guī)律、評估冰湖潰決災(zāi)害風(fēng)險對區(qū)域水資源管理、氣候變化響應(yīng)研究和防災(zāi)減災(zāi)等具有重要意義，而冰湖提取則是以上所有研究的重要基礎(chǔ)。

水體在近紅外波段的地表反射率顯著低于藍(lán)光和綠光波段，由藍(lán)/綠波段和近紅外波段構(gòu)成的水體指數(shù)（Normalized Difference Water Index， NDWI）可以凸顯遙感圖像上的水體區(qū)域。因此，基于NDWI 的閾值分割是常見的冰湖提取方法［6-7］。閾值法簡單快速，適用于較小的區(qū)域，當(dāng)冰湖制圖區(qū)域擴(kuò)大時，由于不同類型、不同環(huán)境下的冰湖光譜并不相同，僅使用單一閾值分割往往造成部分區(qū)域冰湖的過高或過低估計。全域-局部迭代分割法［8］解決了這一問題，該方法根據(jù)初始分割得到的湖泊單元及其緩沖區(qū)生成自適應(yīng)于該湖泊的最佳水體指數(shù)分割閾值，能更好地將冰湖從其背景中分離出來，得到了廣泛的應(yīng)用［9-10］。在局部分割階段，除基于閾值的分割方法外，需要更復(fù)雜計算的活動輪廓模型也能很好地提取冰湖［11-12］。全域-局部迭代法的初始分割結(jié)果一般由簡單的閾值法得到，由于冰湖和山體陰影的光譜較為相似［13］，僅對水體指數(shù)或波段比值設(shè)置閾值無法較好地區(qū)分冰湖和山體陰影。因此，該方法需要利用高程數(shù)據(jù)（Digital Elevation Model，DEM）來消除被錯分的山體陰影，受高程數(shù)據(jù)的精度和相應(yīng)閾值的設(shè)定影響較大。隨著計算機硬件和深度學(xué)習(xí)算法的發(fā)展，近些年涌現(xiàn)出許多利用深度學(xué)習(xí)方法（如U-NET 網(wǎng)絡(luò)［14］）提取冰湖的研究［15-18］，深度學(xué)習(xí)方法能自適應(yīng)地從數(shù)據(jù)中學(xué)習(xí)到適用于目標(biāo)任務(wù)的特征表征，但該類方法需要大量樣本驅(qū)動，與常規(guī)方法相比，所需計算量極大且運行速度慢。

在實際研究中，冰湖的自動化提取過程可分為兩個步驟，首先解決冰湖在哪的問題，即定位冰湖；然后解決冰湖具體形狀如何的問題，即生成精準(zhǔn)的冰湖邊界。如在上述方法中，山體陰影的去除是為了讓冰湖的定位更準(zhǔn)確，而局部迭代分割是為了讓冰湖的形狀更精確，深度學(xué)習(xí)方法則同時解決這兩個問題。近些年來，在諸多學(xué)者的努力下，越來越多的冰湖編目數(shù)據(jù)被生產(chǎn)和公開，如高亞洲地區(qū)1990年和2018年的冰湖編目［19］、高亞洲地區(qū)2008—2017年間每年的冰湖編目［20］。由于冰湖每年的新增數(shù)量有限，可以說這些高質(zhì)量的冰湖編目在很大程度上解決了地物是不是冰湖的問題，在此基礎(chǔ)上提取冰湖只需進(jìn)行高質(zhì)量的圖像分割。因此，這些冰湖編目可以作為生產(chǎn)未來冰湖數(shù)據(jù)的基礎(chǔ)，尤其當(dāng)學(xué)者感興趣的冰湖是某種特定類型的冰湖時，如要對曾發(fā)生過潰決災(zāi)害的冰湖進(jìn)行動態(tài)監(jiān)測，這些冰湖的位置通常是已知的，并且存在一定年份的相應(yīng)歷史邊界數(shù)據(jù)。

圖1展示了冰湖歷史邊界（位置）已知情況下人工閾值法提取冰湖的流程，盡管由于設(shè)置的NDWI閾值偏低導(dǎo)致圖像上大面積的積雪、陰影也隨冰湖一起被提出，但該圖像瓦片是以待提取冰湖歷史邊界為中心的（詳見2.1 節(jié)），只需保留提取結(jié)果中心位置的連通區(qū)，即可剔除與真實冰湖不相連的積雪、陰影等地物。這表明，冰湖的歷史邊界（位置）信息，極大地方便了冰湖當(dāng)前邊界的更新，然而，目前仍然缺乏在已知冰湖歷史邊界情況下的冰湖提取方法對比分析研究，而且即使已知冰湖歷史邊界，在冰凍、積雪、云、山體陰影等因素影響下各方法提取冰湖的適用性也尚不明確。因此，本研究對比評估了冰湖歷史邊界已知情況下人工閾值法、OTSU閾值法、U-NET模型和三種全域-局部迭代分割法的冰湖提取效果，為已知歷史邊界冰湖的長期監(jiān)測提供技術(shù)支持。

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

1 研究區(qū)概況和數(shù)據(jù)來源

本研究區(qū)域位于喜馬拉雅山地區(qū)，地理位置為27°23′～28°28′ N、85°55′～87°58′ E，如圖2 所示。該區(qū)域地形起伏較大，平均海拔約為4 500 m，最高海拔可達(dá)到8 771 m。區(qū)域內(nèi)共有1 509 條冰川，總面積為2 078.17 km2。王欣團(tuán)隊2018 年的冰湖編目［19］顯示該區(qū)域發(fā)育種類多樣、大小不一的冰湖：與冰川連接的湖泊109個，與冰川不連接的湖泊540個，非冰川供給的湖泊296 個，冰面湖118 個；這些冰湖大小各異，最小的冰湖僅有0.0054 km2，最大的冰湖為3.8831 km2；近30年間，該地區(qū)冰湖面積迅速擴(kuò)張，1990年、2018年冰湖面積分別為67.5057 km2、83.5022 km2，相對面積增長了23.7%。

圖2 研究區(qū)及冰湖分布Fig. 2 Study area and the distribution of glacial lakes

Landsat 系列衛(wèi)星數(shù)據(jù)具有長達(dá)幾十年的地球觀測記錄，廣泛用于冰湖監(jiān)測研究。為了綜合評估各方法的性能，研究采用了5景Landsat-8 OLI影像，其元數(shù)據(jù)詳見表1，采用該區(qū)域1990 年的冰湖編目［19］（https：//doi. org/10.12072/casnw. 064.2019.db）作為初始冰湖邊界，受氣候變化影響，該編目記錄的冰湖形狀、面積與2015 年的冰湖有較大差異，有利于評估歷史邊界已知情況下各方法的冰湖提取性能。序號為1 的影像，包含488 個已知1990 年邊界的冰湖，用于各方法性能的整體評估，序號為2到5 的影像，分別包含20 個已知1990 年邊界的冰湖，用于評估各方法對冰凍狀態(tài)、積雪覆蓋、與山體陰影連接、云覆蓋的冰湖的提取性能。這些冰湖（共568個）的對應(yīng)矢量均由有經(jīng)驗的操作員在Arc-GIS 10.6軟件中繪制。

表1 使用影像詳細(xì)信息Table 1 Details of the images used

2 研究方法

2.1 技術(shù)路線

本研究的整體技術(shù)路線如圖3 所示。首先對Landsat-8 OLI 影像進(jìn)行輻射定標(biāo)，即將影像DN 值轉(zhuǎn)換為具有物理意義的大氣頂層（Top of Atmosphere，TOA）反射率，用于生成NDWI 圖像。TOA反射率影像已經(jīng)可以用于提取冰湖邊界，與文獻(xiàn)［22］一致，本研究不對Landsat-8 OLI 影像進(jìn)行大氣校正。DN值到TOA值的轉(zhuǎn)化公式如下：

圖3 技術(shù)路線圖Fig. 3 The technical flow chart

式中：增益Gain 和偏置Offset 可以從影像的元數(shù)據(jù)中獲取。本研究采用McFeeters提出的NDWI［21］：

式中：ρGreen和ρNIR分別為綠光波段和近紅外波段的TOA反射率。

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

2.2 對比方法

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

（2）OTSU 閾值法。有研究利用OTSU［23］閾值法結(jié)合圖像的近紅外波段進(jìn)行冰湖提取［24］。為公平對比各方法，本研究將OTSU 閾值法應(yīng)用于NDWI 圖像。OTSU 方法根據(jù)最大化類間距離且最小化類內(nèi)方差的原則確定最優(yōu)的灰度圖像二分類閾值，可以避免閾值選擇的主觀性。

（3）全域-局部分割法。全域-局部的方法在全域分割階段總是設(shè)定較低的NDWI 閾值以盡可能地分割出圖像包含的所有冰湖。在得到初始冰湖單元后，可應(yīng)用不同的方法對每個冰湖單元進(jìn)行局部分割，下面介紹三種不同的局部分割方法。

原創(chuàng)的全域-局部方法［22］中，局部分割分為三步：（i）對當(dāng)前冰湖單元建立等面積的緩沖區(qū)；（ii）對緩沖區(qū)內(nèi)像元的NDWI 值進(jìn)行直方圖統(tǒng)計，按照雙峰分布準(zhǔn)則確定分割閾值并進(jìn)行分割得到新的冰湖單元；（iii）重復(fù)（i）到（ii）步驟，直至冰湖像元數(shù)穩(wěn)定。雙峰分布準(zhǔn)則的閾值確定公式為：

式中：μWater和μLand分別為當(dāng)前冰湖單元和其緩沖區(qū)內(nèi)像元的NDWI 均值；σWater和σLand分別為其對應(yīng)的標(biāo)準(zhǔn)差。

部分研究在第（ii）步中使用OTSU 閾值分割法［24］或C-V 模型［12］進(jìn)行冰湖局部優(yōu)化提取。其中，OTSU 方法的原理見2.2（2）小節(jié)，C-V 模型通過變分法和梯度下降法不斷演化初始邊界C使得能量函數(shù)E最小化，此時的C即可將前景和背景分割。該能量函數(shù)E定義為：

式中：第一項是對邊界C長度的正則化，第二項和第三項是灰度圖像強度和區(qū)域均值差異的積分。f(x，y)代表圖像在x，y處的灰度值；c1代表邊界C內(nèi)部的灰度均值；c2代表邊界C外部的灰度均值。參照文獻(xiàn)［11］，本研究將μ設(shè)置為0.1，λ1、λ2均為1。

在應(yīng)用全域-局部法提取冰湖時，上述兩項研究在利用OTSU 方法和C-V 模型分割出冰湖后，并未進(jìn)行第（iii）步中的迭代過程。當(dāng)初始冰湖單元的等面積緩沖區(qū)完全包含真實冰湖的所有像元時，迭代過程不是必需的。而本研究中冰湖的初始邊界（1990 年）和冰湖真實邊界（2014 年以后）可能相差很大，若不進(jìn)行第（iii）步中的迭代分割，則可能最終只提取出真實冰湖的小部分區(qū)域，因此，本研究增加了第（iii）步的迭代過程。在使用全域-局部法時，本文利用1990年冰湖編目數(shù)據(jù)作為初始冰湖單元，這使得全域分割不再必要，且局部分割都包含迭代分割的過程，因此將本節(jié)介紹的三種方法命名為：雙峰迭代法、OTSU 迭代法和C-V 迭代法。將迭代分割的終止條件設(shè)置為：相鄰兩次分割所得冰湖的相對面積差小于等于1%。

（4）深度學(xué)習(xí)方法U-NET。U-NET 是一種對稱的編碼解碼語義分割網(wǎng)絡(luò)，其核心力量來源于跳連接，結(jié)構(gòu)詳見圖4。U-NET 在編碼階段通過堆疊卷積層和池化層提取輸入圖像的多層級特征，池化層以保留特征圖中各子區(qū)域內(nèi)最強特征的方式減小特征圖尺寸，同時使特征圖丟失了很多細(xì)節(jié)信息，U-NET 在解碼階段通過不斷堆疊上采樣層和卷積層逐步得到和輸入圖像長寬一致的預(yù)測圖，跳連接使得上采樣后的粗分辨率高級語義特征和編碼階段的高分辨率淺層細(xì)節(jié)特征相融合，從而恢復(fù)因池化損失的目標(biāo)細(xì)節(jié)信息，實現(xiàn)地物的精細(xì)分割。UNET 這類基于學(xué)習(xí)的模型總是需要一定規(guī)模的樣本驅(qū)動，本研究將序號為1 的影像分為A、B 兩個區(qū)域（圖2），分別包含237、251 個冰湖，使用A 區(qū)域的237 個圖像瓦片訓(xùn)練U-NET 并對B 區(qū)域的251 個冰湖進(jìn)行預(yù)測，再使用B 區(qū)域的251 個圖像瓦片訓(xùn)練U-NET 并對A 區(qū)域的237 個冰湖進(jìn)行預(yù)測，由此得到序號為1 的影像上488 個冰湖的預(yù)測結(jié)果。為了提取冰凍狀態(tài)、積雪覆蓋、云覆蓋、與山體陰影連接的共80 個冰湖，使用A 區(qū)域和B 區(qū)域的488 個圖像瓦片訓(xùn)練U-NET。由于本研究中用于訓(xùn)練的圖像瓦片數(shù)量較少，且為NDWI 單波段圖像，研究將UNET 編碼階段濾波器（卷積核）的數(shù)量設(shè)置為32、64、128、256、512（圖4）。在訓(xùn)練U-NET 時采用SGD（Stochastic Gradient Descent）作為優(yōu)化器，動量設(shè)置為0.9，由于冰湖和背景地物的面積占比極不均衡，使用Dice 損失作為網(wǎng)絡(luò)的損失函數(shù)，Batch 大小設(shè)置為4，初始學(xué)習(xí)率設(shè)置為0.0025，采用Poly策略逐漸降低學(xué)習(xí)率，共訓(xùn)練50輪?？紤]到樣本數(shù)量較少，在訓(xùn)練過程中采用隨機水平翻轉(zhuǎn)和豎直翻轉(zhuǎn)進(jìn)行樣本增強。圖5 展示了利用A 區(qū)域和B 區(qū)域圖像瓦片訓(xùn)練U-Net的損失曲線。

圖4 U-NET結(jié)構(gòu)圖Fig. 4 The Structure of U-NET

圖5 U-Net的訓(xùn)練損失曲線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 結(jié)果分析

3.1 成像條件較好情況下各方法的冰湖提取結(jié)果

冰湖提取結(jié)果的精度評價采用精確率（Precision）、召回率（Recall）和F1 分?jǐn)?shù)三個指標(biāo)，其定義參考文獻(xiàn)［25］。表2 呈現(xiàn)了2.2 節(jié)中所述六種方法對序號為1 的Landsat 影像上488 個冰湖的提取精度。人工閾值法提取的冰湖常包含與真實冰湖相連的冰、雪［圖6（c）第3 行］、陰影和部分陸地［圖6（c）第2行］等NDWI值大于0的地物，因此該方法的精確率只有66.02%，且F1 分?jǐn)?shù)最低。OTSU 閾值法利用了圖像瓦片全局的NDWI 統(tǒng)計信息，當(dāng)真實冰湖和圖像其他區(qū)域的NDWI 值有明顯差異時，該方法可以很好地分割出冰湖［圖6（d）第2 行和第4行］。然而，當(dāng)圖像可按NDWI 值從高到低分為三類（冰湖、陰影/云/雪、陸地）時，該方法提取的冰湖也會包含與真實冰湖連接的陰影、雪［圖6（d）第1行和第3 行］。因此，OTSU 閾值法的精確率也比較低，只有75.76%。值得注意的是，這兩種方法的召回率都很高，分別為96.70%和94.25%，即能夠盡可能地提取出所有冰湖。

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

圖6 六種方法的冰湖提取結(jié)果示例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 迭代法取得了非常高的整體精度，F(xiàn)1分?jǐn)?shù)分別達(dá)到88.89%和89.30%，這是因為這兩種方法在當(dāng)前冰湖及其等面積的緩沖區(qū)內(nèi)迭代地提取冰湖。冰湖和其等面積緩沖區(qū)內(nèi)像元的NDWI值往往呈現(xiàn)雙峰分布，與基于全局NDWI 的OTSU閾值法相比，OTSU 迭代法能給出更恰當(dāng)?shù)姆指铋撝?。對于C-V 迭代法而言，等面積緩沖區(qū)使得式（4）中第二、三項在數(shù)量上更均等，更利于邊界的演化。同樣以1990 年冰湖作為初始冰湖單元的雙峰迭代法卻無法完整地提取出冰湖［圖6（e）］，該方法的F1 分?jǐn)?shù)只有77.67%，原因是雙峰分布的閾值確定準(zhǔn)則式（3）不適用于初始冰湖邊界（1990 年）和真實冰湖邊界（2015年）相差很大的情況，詳見4.1節(jié)。U-NET 方法的冰湖提取精度最高，其F1 分?jǐn)?shù)為89.8%，但對于面積>0.1 km2的冰湖，其最終提取的冰湖與真實冰湖相比總是存在缺失［圖6（h）第1 行和第2行］，這可能是由小冰湖在訓(xùn)練樣本中占比較高（圖7）導(dǎo)致。

圖7 冰湖面積和提取精度間的關(guān)系Fig. 7 The relationship between glacial lake areas and extraction accuracies of OTSU iterative method，C-V iterative method and U-NET

為了進(jìn)一步評估各方法的冰湖提取性能，按照面積將上述488個冰湖分為三類：小型冰湖（≤0.01 km2），中型冰湖（0.01～0.1 km2），大型冰湖（>0.1 km2）。由表3可知，受影像空間分辨率（30 m）的影響，各方法對大型冰湖的提取精度（F1 分?jǐn)?shù)）往往大于對中型冰湖的提取精度，對小型冰湖的提取精度則最低。對三類不同大小的冰湖，OTSU 迭代法、C-V 迭代法和U-NET 取得的F1 分?jǐn)?shù)明顯高于人工閾值法、OTSU 閾值法和雙峰迭代法，這與表2 呈現(xiàn)的特點一致。具體地，OTSU 迭代法、C-V 迭代法和UNET 在小型冰湖、中型冰湖和大型冰湖上取得的F1分?jǐn)?shù)分別集中在82%、88%和94%附近，對小型冰湖而言，OTSU 迭代法取得的F1 分?jǐn)?shù)最高（82.79%），對中型冰湖而言，U-NET 取得的F1 分?jǐn)?shù)最高（89.69%），對大型冰湖而言，C-V 迭代法取得的F1 分?jǐn)?shù)最高（94.91%），整體上，這三種方法都適用于對各類不同大小冰湖的提取。圖7更詳細(xì)地展示了OTSU 迭代法、C-V迭代法、U-NET對不同大小冰湖的提取精度，盡管這三種方法取得的F1分?jǐn)?shù)都在89%左右，仍然有一些小面積冰湖的F1分?jǐn)?shù)低于0.6，原因是小冰湖無法提供足夠的統(tǒng)計信息。

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

3.2 復(fù)雜成像條件下各方法的冰湖提取結(jié)果

表4～7 展示了本文所述六種方法對冰凍狀態(tài)、積雪覆蓋、云覆蓋、與山體陰影連接的冰湖的提取精度。盡管冰凍狀態(tài)冰湖的NDWI 值有所降低，與周圍背景地物的NDWI 值仍有較明顯差異［圖8（b）］，OTSU 迭代法、C-V 迭代法和U-NET 都能很好地提取冰凍狀態(tài)的冰湖，且與表2所示結(jié)果一致，這三種方法的冰湖提取精度明顯優(yōu)于人工閾值法、OTSU 閾值法和雙峰迭代法，其中C-V 迭代法取得了最高的F1 分?jǐn)?shù)（90.61%）。與雙峰迭代法和OT-SU 迭代法這類基于閾值分割迭代提取冰湖的方法相比，C-V 迭代法不受冰湖內(nèi)部NDWI 值較低的部分像素影響，能完整地提出冰湖（圖8 第三行）。不同于冰凍狀態(tài)的冰湖，積雪覆蓋下冰湖的NDWI 值極低，與周圍背景的NDWI 值十分接近［圖9（b）］，因此，人工閾值法、OTSU 閾值法、雙峰迭代法、OTSU 迭代法和U-NET 的冰湖提取結(jié)果都較差，其中U-NET 取得的F1 分?jǐn)?shù)僅有13.57%，只能將冰湖中未被積雪覆蓋的一小部分提出［圖9（i）］。雖然積雪覆蓋下冰湖的NDWI 值極低，但冰湖邊界仍然較為清晰，故C-V 迭代法取得了比其他方法顯著更高的F1分?jǐn)?shù)（81.40%）。

表4 六種方法對冰凍狀態(tài)冰湖的提取精度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 冰凍狀態(tài)冰湖提取結(jié)果示例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 積雪覆蓋冰湖提取結(jié)果示例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 分?jǐn)?shù)（89.51%），U-NET 則幾乎無法提出云覆蓋下的冰湖［圖10（i）］，取得的F1分?jǐn)?shù)僅有2.26%。當(dāng)冰湖與山體陰影相連接時，由于兩者的NDWI 值十分接近，人工閾值法和OTSU閾值法都會誤提大量山體陰影［圖11（d）和（e）］，即使OTSU 迭代法和C-V迭代法也無法較好地分離二者［圖11（g）和（h）］，僅U-NET這種基于學(xué)習(xí)的方法能較完整地分割出冰湖，取得的F1 分?jǐn)?shù)高達(dá)85.19%，顯著優(yōu)于其他方法。

圖10 云覆蓋冰湖提取結(jié)果示例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 年冰湖作為初始冰湖單元，且都基于閾值對緩沖區(qū)內(nèi)的像元進(jìn)行分割，二者唯一的區(qū)別在于分割閾值的計算公式不同。本節(jié)以圖12展示的冰湖為例，探究兩種方法表現(xiàn)差異的原因。

圖12 兩種迭代法的冰湖提取結(jié)果示例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年冰湖為初始冰湖單元進(jìn)行第一次迭代分割時，兩種方法要處理的冰湖等面積緩沖區(qū)的NDWI 值完全一致，但基于雙峰分布的閾值確定準(zhǔn)則式（3）給出的NDWI 分割閾值（0.551）明顯高于OTSU 的閾值（0.357）。這是因為1990 年冰湖只是2015 年冰湖的一小部分，由此得到的等面積緩沖區(qū)除包含陸地像元外還有大量的水體像元，導(dǎo)致式（3）計算的閾值更偏向直方圖中的水體（圖13中迭代次數(shù)i=0的直方圖）。較高的閾值會導(dǎo)致冰湖漏分。因此，本次分割得到的新冰湖及其等面積緩沖區(qū)會包含更高比例的水體像元，并再次由式（3）得到更高的分割閾值（0.647），導(dǎo)致冰湖進(jìn)一步漏分。

圖13 雙峰迭代法的冰湖迭代提取過程，黑白二值圖代表第i次迭代分割后的冰湖，與其對應(yīng)的直方圖顯示該冰湖等面積緩沖區(qū)內(nè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次迭代分割后的冰湖，與其對應(yīng)的直方圖顯示該冰湖等面積緩沖區(qū)內(nè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 所示冰湖進(jìn)行提取時，每次迭代前冰湖和緩沖區(qū)的像元數(shù)、冰湖和緩沖區(qū)的NDWI 均值和標(biāo)準(zhǔn)差以及迭代后冰湖的像元數(shù)。如圖13 所示，雙峰迭代法從第1 次到第5 次迭代給出的分割閾值從0.551 增加至0.686，這期間分割所得的冰湖范圍不斷增大；從第6 次到第10 次迭代，分割閾值穩(wěn)定在0.69 左右，并且分割所得冰湖面積不斷減小，直至達(dá)到穩(wěn)定的漏分狀態(tài)（圖13 迭代次數(shù)i=11）。在此過程中，冰湖緩沖區(qū)內(nèi)的NDWI 直方圖從不太顯著的雙峰分布逐漸演化為只包含水體像元的單峰分布。與雙峰分布的閾值確定準(zhǔn)則相比，OTSU 在緩沖區(qū)內(nèi)水體像元較多的情況下仍能給出恰當(dāng)?shù)姆指铋撝担@是由于OTSU 準(zhǔn)則在計算二分類閾值過程中不涉及由當(dāng)前冰湖或緩沖區(qū)計算的參量。恰當(dāng)?shù)姆指铋撝凳沟梅指詈蟊牡让娣e緩沖區(qū)包含更多的陸地像元，因此OTSU 迭代法給出的分割閾值從第1 次迭代的0.357降至第5次迭代的0.3（圖14），且直方圖呈現(xiàn)的雙峰分布態(tài)勢在迭代過程中逐步明顯。

表8 雙峰迭代法的過程參數(shù)Table 8 The process parameters of the bimodal iterative method

4.2 C-V模型的參數(shù)設(shè)置

在上述三種迭代法中，C-V 迭代法取得的F1 分?jǐn)?shù)最高（89.30%），與OTSU 迭代法相比僅提升了0.40%，但在每次迭代過程中C-V 模型都會對初始邊界進(jìn)行多次演化，因此C-V 迭代法比OTSU 迭代法需要的時間成本更高。OTSU 方法確定分割閾值的準(zhǔn)則是無參數(shù)的，而C-V 模型的演化則需要設(shè)定三個參數(shù)，即式（4）中的μ、λ1、λ2。通常λ1、λ2都設(shè)置為1，μ則根據(jù)待分割對象的幾何特征設(shè)定。待分割對象邊界越圓滑則μ越大，待分割對象邊界越不規(guī)則（周長相對較長）或面積較小，則μ越小。本研究中待分割的488 個冰湖大小各異，μ值的設(shè)置直接影響C-V 迭代法的分割結(jié)果（如表9 所示）。μ為0.05 時，C-V 迭代法取得的F1 分?jǐn)?shù)最高（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 迭代法的冰湖提取結(jié)果。對于小冰湖（第一行），僅當(dāng)μ=0.05 時C-V迭代法才能完整地提取冰湖；對于大冰湖（第二行），不論μ取何值，C-V 迭代法都能將冰湖完整地提出，但隨著μ值的增大，提取的冰湖邊界越來越光滑。如表10所示，對于小型冰湖，μ=0.05時的F1分?jǐn)?shù)比μ=0.20 時高6.86%；對于中型冰湖，μ=0.05 時的F1分?jǐn)?shù)比μ=0.20 時高2.67%；當(dāng)冰湖面積>0.1 km2時，μ=0.05 時的F1 分?jǐn)?shù)比μ=0.20 時低0.21%。由此可知，μ的取值對面積越小的冰湖影響越大，小冰湖應(yīng)設(shè)置較低的μ值，大冰湖則對μ的取值不敏感。因此，建立冰湖面積與μ值之間的關(guān)系，在對不同大小的冰湖進(jìn)行提取時設(shè)置相應(yīng)的最佳μ值，將進(jìn)一步提升C-V迭代法的冰湖提取精度。

表10 μ=0.05和μ=0.20時C-V迭代法提取不同大小冰湖的精度統(tǒng)計Table 10 Extraction accuracies of glacial lakes in different sizes by C-V iterative method with different μ values

圖15 不同μ值的C-V迭代法的冰湖提取結(jié)果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個圖像瓦片上訓(xùn)練的UNET能較好地提取冰凍狀態(tài)的冰湖和與山體陰影連接的冰湖，其取得的F1 分?jǐn)?shù)分別為88.20% 和85.19%，但對于積雪覆蓋和云覆蓋冰湖的提取，UNET表現(xiàn)十分糟糕，與其他方法相比其取得的F1分?jǐn)?shù)最低。這可能是由用于訓(xùn)練和測試的冰湖的NDWI值分布差異較大導(dǎo)致，經(jīng)統(tǒng)計，冰凍狀態(tài)和與山體陰影連接冰湖的NDWI 均值分別為0.4912 和0.7908，與序號為1的圖像中的488個冰湖的NDWI均值0.6175較為接近，而積雪覆蓋和云覆蓋下冰湖的NDWI 均值分別為0.0841 和-0.0003，與0.6175差距極大。因此，若想使得U-NET 這類深度學(xué)習(xí)方法在各種場景下都取得優(yōu)異的冰湖提取效果，構(gòu)建包含各種類型、不同成像條件的冰湖樣本集至關(guān)重要。

5 結(jié)論

本研究利用喜馬拉雅地區(qū)2014年后的Landsat-8 OLI 影像，基于1990 年冰湖編目數(shù)據(jù)提供的歷史邊界（位置）信息，對成像條件較好的488 個冰湖和受積雪、冰凍、云、山體陰影影響的80個冰湖開展了冰湖提取方法對比實驗，結(jié)論如下：

（1）U-NET、OTSU 迭代法、C-V 迭代法都能精確地提取成像條件較好的冰湖，F(xiàn)1 分?jǐn)?shù)均達(dá)到88.80%以上，顯著優(yōu)于人工閾值法和OTSU 閾值法。由于1990 年與2015 年冰湖范圍差異較大，利用冰湖及其緩沖區(qū)內(nèi)統(tǒng)計信息的雙峰迭代法易造成過分割（漏分）現(xiàn)象。

（2）在提取冰凍狀態(tài)、積雪覆蓋和云覆蓋下的冰湖時，應(yīng)優(yōu)先采用C-V 迭代法。而U-NET 受樣本分布影響較大，盡管無法有效提取積雪覆蓋和云覆蓋下的冰湖，但相比其他方法，能夠更有效地分離冰湖和與之連接的山體陰影。

（3）在提取大型、中型和小型冰湖時，優(yōu)先采用的方法分別為C-V 迭代法、U-NET 和OTSU 迭代法。由于本研究的實驗數(shù)據(jù)為中分辨率Landsat 8影像，六種方法對大型冰湖的提取精度顯著高于中、小型冰湖，可采用更高空間分辨率的遙感影像，如Sentinel-2（10 m）或高分2 號數(shù)據(jù)（4 m）精細(xì)刻畫中、小型冰湖的邊界。

（4）冰湖在不同分辨率影像上的位置是相對不變的，本研究所述基于歷史邊界的冰湖提取方法也適用于更高分辨率的影像?，F(xiàn)存的冰湖歷史邊界多來自Landsat 數(shù)據(jù)，盡管這些邊界信息不夠精細(xì)，但可作為OTSU 迭代法和C-V 迭代法的迭代起點。今后將開展基于高分辨率遙感影像的冰湖提取方法對比和優(yōu)化研究，以實現(xiàn)中、小型冰湖邊界的精細(xì)提取。