基于遙感影像的山岳冰川信息計算機解譯方法探討
——以梅里雪山為例

鐘鼎杰，楊存建**

(1. 四川師范大學 西南土地資源評價與監(jiān)測教育部重點實驗室，四川 成都 610068；2. 四川師范大學 地理與資源科學學院，四川 成都 610068)

冰川不僅存貯了大量的水資源，同時也是全球氣候變化的指示器[1]. 在全球氣候變暖的背景下對于冰川變化的研究迫在眉睫，在信息化時代如何運用高效、便捷、精準的方法對冰川變化做出監(jiān)測與預(yù)測是目前的研究熱點. 相比大陸型冰川，海洋型冰川作用能大、運動速度快，對于氣候變化的響應(yīng)也更為迅速；海洋型冰川在大量降水與劇烈消融的作用下會產(chǎn)生破壞性的泥石流與雪崩，所以對于海洋型冰川變化的研究具有重要的意義[2]. 在冰川信息提取的方法上，近年來遙感技術(shù)、地理信息系統(tǒng)技術(shù)的不斷發(fā)展和集成應(yīng)用，可利用遙感影像和DEM 等數(shù)據(jù)快速高效地提取冰川信息并監(jiān)測冰川面積、高程和冰川末端變化. 1999 年國際全球陸地冰 空 間 監(jiān) 測 計 劃（Global Land Ice Measurements from Space Project，GLIMS）運 用 了 衛(wèi) 星 傳 感 器ASTER 數(shù)據(jù)和Landsat 數(shù)據(jù)對世界范圍內(nèi)的冰川進行動態(tài)監(jiān)測[3]. Sidjak 等[4]利用了監(jiān)督分類法、TM 影像的4、5 波段比值以及歸一化積雪指數(shù)（Normalized Difference Snow Index，NDSI）等方法對加拿大冰川國家公園進行了冰川信息提取.Kumar 等[5]采用干涉合成孔徑雷達（Interfermetric Synthetic Aperture Radar，InSAR）技術(shù)對喜馬拉雅地區(qū)冰川表面速度進行了估算. Robson 等[6]利用基于面向?qū)ο蟮膱D像分析（Object-based Image Analysis，OBIA）技術(shù)，結(jié)合Landsat 8 影像、DEM數(shù)據(jù)和ALOS-PALSAR 影像對喜馬拉雅山脈馬納斯魯峰進行了冰川信息的自動提取. 張世強等[7]利用TM 圖像并采用閾值法、譜間關(guān)系法、監(jiān)督分類和非監(jiān)督分類提取青藏高原喀喇昆侖山地區(qū)冰川信息，證明利用比值圖像取閾值是提取冰川信息最有效的手段. 懷保娟等[8]利用Landsat TM/ETM+遙感影像和DEM 數(shù)據(jù)，并運用了面向?qū)ο蟮膱D像信息自動提取方法，對黑河流域近50 年來的冰川變化進行了分析. 聶勇等[9]利用多光譜遙感數(shù)據(jù)，采用面向?qū)ο蟮慕庾g方法，結(jié)合專家知識分類規(guī)則自動提取1976、1988 年和2006 年3 個時期珠穆朗瑪峰地區(qū)的冰川信息. 李宗省等[10]對1900—2007 年橫斷山區(qū)部分海洋型冰川變化進行了研究，并提出各冰川的變化幅度因冰川規(guī)模、緯度位置等因素而存在明顯差異. 李霞等[11]對貢嘎山地區(qū)近40 年來的冰川變化進行了遙感監(jiān)測研究，并從不同的角度研究分析了冰川變化的特征.

雖然冰川信息提取方法有多種，但因冰川類型的不同和地域性差異，并沒有通適性較高的方法，冰川監(jiān)測仍以目視解譯為主. 因此，研究分類精度較高、通適性強且高效的冰川信息的計算機解譯方法意義重大. 本文在已有的計算機自動解譯冰川信息方法研究上，納入了神經(jīng)網(wǎng)絡(luò)分類方法，通過對比分析監(jiān)督分類、分監(jiān)督分類、比值閾值法、雪蓋指數(shù)法（NDSI）、基于多尺度分割的面向?qū)ο蠓ǎ∣BIA）和基于神經(jīng)網(wǎng)絡(luò)的冰川識別方法的信息提取結(jié)果，探討每種方法的優(yōu)缺點及適用情況，并將梅里雪山地區(qū)冰川信息的目視解譯作為地面真實數(shù)據(jù)，對6 種方法的提取結(jié)果進行精度驗證，發(fā)現(xiàn)了神經(jīng)網(wǎng)絡(luò)分類法提取的冰川信息精度最高. 在此研究基礎(chǔ)上，基于ENVI 深度學習模塊，利用神經(jīng)網(wǎng)分類法解譯了1989、1998、2009 年和2019 年的梅里雪山地區(qū)冰川信息. 最終，結(jié)合Google Earth數(shù)據(jù)及DEM 數(shù)據(jù)對梅里雪山地區(qū)冰川信息進行目視修正，并通過DEM 數(shù)據(jù)分析獲取梅里雪山地區(qū)山脊線，將提取的結(jié)果分割為55 條獨立的冰川，得到了1989—2019 年梅里雪山地區(qū)較為精確的4期冰川信息的提取成果，并在此基礎(chǔ)上分析統(tǒng)計其冰川面積的變化量與變化速率. 結(jié)果表明，基于神經(jīng)網(wǎng)絡(luò)的冰川識別方法分類精度、解譯效果明顯優(yōu)于傳統(tǒng)冰川信息計算機解譯方法，且相較于目視解譯該方法效率更高. 在全球氣候變暖背景下，本文基于神經(jīng)網(wǎng)絡(luò)的冰川識別方法可為全球冰川時空演化的長時序監(jiān)測提供有效的技術(shù)支撐.

1 研究區(qū)概況

橫斷山地處青藏高原東緣，是中國地勢一、二階梯的過渡地帶，區(qū)域大致包括了藏東南、云南西北部和四川的阿壩、涼山及甘孜地區(qū)，總面積約為5×105km2. 橫斷山區(qū)南北走向的高山和深谷相間排列，是本區(qū)域最為顯著的地貌特征. 梅里雪山位于橫斷山東南緣，也是橫斷山地區(qū)集中低緯度海洋性/暖性冰川最多的山段[12]（圖1）. 梅里雪山高度大于6 000 m 的山峰有13 座，最高峰卡瓦格博海拔6 740 m，是云南省第一高峰，山腳瀾滄江海拔2 020 m，相對高差4 720 m.

圖1 研究區(qū)位置與冰川積雪分布圖Fig. 1 Location of study area and distribution map of glacier and snow cover

據(jù)中國第二次冰川編目數(shù)據(jù)集（V1.0）[13-16]，梅里雪山地區(qū)共有冰川55 條. 梅里雪山地區(qū)的冰川受來自印度洋西南季風帶來的暖濕氣流影響，隨著暖濕氣流的不斷深入，給梅里雪山提供了充沛的降水，為冰川提供了豐富的物質(zhì)來源.

2 數(shù)據(jù)獲取與處理

2.1 數(shù)據(jù)來源本文選用的遙感數(shù)據(jù)均為通過美國地質(zhì)勘探局（https://earthexplorer.usgs.gov）獲取的Landsat Level 2 產(chǎn) 品. 數(shù) 據(jù) 為Landsat 5 搭 載 的TM 傳感器、Landsat 8 搭載的陸地成像儀（Operational Land Imager，OLI）傳感器所獲取的影像（表1）.在實際的數(shù)據(jù)選擇過程中，需要特別注意研究區(qū)的積雪覆蓋是否會影響冰川解譯精度以及云層是否覆蓋了研究區(qū).

表1 梅里雪山遙感影像數(shù)據(jù)信息Tab. 1 Remote sensing image data information in Meili Snow Mountain

DEM 數(shù)據(jù)來源于美國地質(zhì)勘探局的ASTER GDEM V2 版數(shù)據(jù)，其空間分辨率為30 m. 高分辨率可見光影像主要來自于Google Earth，用于輔助目視解譯. 地圖的底圖來源于國家科技資源共享服務(wù)平臺的國家青藏高原科學數(shù)據(jù)中心（http://data.tpdc.ac.cn/zh-hans/）.

2.2 傳統(tǒng)分類方法提取冰川信息

2.2.1 非監(jiān)督分類 是指在沒有先驗類別樣本的情況下，根據(jù)像元間地物光譜相似度統(tǒng)計分布規(guī)律將其自動分類. 常用的非監(jiān)督分類方法有ISODATA 和K-Means 兩種，K-means 算法 初始隨機選取數(shù)據(jù)集中K個點作為聚類中心，以所有樣本到聚類中心距離平方和最小的為分類結(jié)果. 在Kmeans 中，K值需要預(yù)先人為確定，并且在整個算法過程中無法更改；ISODATA 算法針對這個問題進行了改進，并增加了“合并”和“分裂”操作. 本文選用ISODATA 算法，對梅里雪山地區(qū)冰川信息進行提取. ISODATA 算法思想為輸入N個樣本，預(yù)選Nc個初始聚類中心和各參數(shù)指標K、θN、θS、θc、L和I（表2）. 若 θc小于閥值則歸為一類，否則為不同類；若 θN小于閥值，則將其取消[17].

表2 初始聚類中心和參數(shù)指標Tab. 2 Initial cluster center and parameter index

2.2.2 監(jiān)督分類 是在已知類別的訓練場地上提取不同種類的訓練樣本，通過選擇特征變量、確定判別函數(shù)或判別規(guī)則，從而把圖像中的各個像元點劃歸到各個給定類的分類方法[18]. 本文選用最大似然法對梅里雪山地區(qū)進行監(jiān)督分類，在訓練樣本時將訓練區(qū)內(nèi)的地物分為裸地、植被和冰川. 最大似然法算法思想是通過遙感影像的統(tǒng)計特征，假設(shè)各類的分布函數(shù)為正態(tài)分布，在多變量空間中形成橢圓（球）分布，按正態(tài)分布規(guī)律用最大似然判別規(guī)則進行判決，得到較準確的分類結(jié)果.

2.2.3 比值閾值法 是將兩個具有一定差異的光譜波段反射率相除，比值增強波譜特征的微小差異，該算法對不同波段反射率差異大的地物效果明顯. 對于冰川提取，通常使用的是對于冰川反射率較高的可見光波段以及對于冰川反射率較低的波段，一般是近紅外波段[19]. 其公式表示為

2.2.4 雪蓋指數(shù)法 受到植被指數(shù)法的啟發(fā)，Dozier[20]提出了雪蓋指數(shù)法（NDSI）作為冰川信息提取的方式，主要原理基于冰川表面在可見光波段具有強反射的特性，以及在近紅外波段具有的強吸收特性. 將兩個波段反射率進行歸一化處理，以突出冰、雪的信息特性. 其公式表示為

常用的Landsat TM 影像波段有TM2 和TM5，Landsat OLI 影像波段選用Band 3 和Band 6. NDSI的閾值一般根據(jù)具體的情況結(jié)合經(jīng)驗值確定，本文根據(jù)前人經(jīng)驗及多次試驗結(jié)果，選取的二值化分類閾值為0.6[21].

2.3 基于多尺度分割的面向?qū)ο蠓ǎ∣BIA）提取冰川信息基于多尺度分割的面向?qū)ο蠓ǎ∣BIA）提取冰川信息是選用多尺度分割算法并結(jié)合遙感影像的光譜特征進行. 當全部待分割對象的綜合加權(quán)值大于某一特定閾值時分割完成；若小于閾值，則重復進行迭代運算，直至條件成立. 本文使用Edge 分割算法和Fast Lambda 合并算法對遙感影像進行處理. 圖像分割完成后，對各波段的光譜、紋理特征等參照直方圖探尋相關(guān)信息的閾值，綜合梅里雪山冰川與其它地物差異性較明顯的冰川信息指數(shù)分別進行閾值設(shè)置并建立提取知識規(guī)則[22]，即①比值閾值TM3/TM5，以閾值2.4 提取冰川信息；②冰雪指數(shù)NDSI，設(shè)置閾值為0.6.

2.4 基于神經(jīng)網(wǎng)絡(luò)的冰川識別方法近些年來，深度學習發(fā)展迅速，在圖像分類方面的應(yīng)用表現(xiàn)優(yōu)異[23-27]. 它通過學習大量樣本，提取樣本中的特征值，形成分類模型，進而對遙感影像進行分類. 在冰川信息提取方面，可以提高冰川識別的精度，并對有冰磧覆蓋物的冰舌部分提取效果較好. 本文通過ENVI 深度學習模塊對梅里雪山地區(qū)冰川信息進行提取研究，該模塊基于TensorFlow 框架，利用ENVI Net5 網(wǎng)絡(luò)構(gòu)建深度學習模型. ENVI Net5 是一種U-Net 架構(gòu)，它是基于掩碼（mask-based）、編碼器-解碼器（encoder-decoder）的體系結(jié)構(gòu)，用于對圖像中每個像素進行分類. 具體的分類步驟如圖2所示.

圖2 基于ENVI 深度學習模塊的圖像信息提取步驟Fig. 2 Image information extraction steps based on ENVI deep learning module

影像分類完成后會得到一個CAM 類激活灰度圖，灰度圖中每個像元大致表示屬于目標類別的概率，值閾范圍為0～1. 通過密度分割，將類激活灰度圖（圖3）按不同顏色顯示，紅色部分表示深度學習識別為冰川的概率較高.

圖3 密度分割后梅里雪山地區(qū)冰川神經(jīng)網(wǎng)絡(luò)分類結(jié)果Fig. 3 Neural network classification results of glaciers in Meili Snow Mountain after density segmentation

3 比較與討論

3.1 分類結(jié)果對比為驗證本文所用方法的有效性，基于目視解釋法梅里雪山地區(qū)冰川信息，以目視解譯結(jié)果作為地面真實數(shù)據(jù)，對計算機解譯結(jié)果進行精度驗證（表3）. 精度評價指標為冰川總面積、面積誤差、重疊部分面積、冰川精度[28-29].

表3 梅里雪山地區(qū)冰川分類精度Tab. 3 Accuracy of glacier classification in Meili Snow Mountain

圖4 至圖9 分別為不同分類方法提取的冰川信息結(jié)果. 其中，黃色為基于目視判讀方法手動提取梅里雪山地區(qū)冰川信息，紅色為對應(yīng)的每種分類方法所提取的冰川信息. 如圖4 所示，非監(jiān)督分類提取冰川信息結(jié)果中，裸地被錯分為冰川的情況較多（如圖4（d）、4（f）），有冰磧覆蓋物的冰川被漏分（如圖4（a）、4（c）），山體陰影區(qū)中的冰川未被識別出（如圖4（b）、4（e））. 如圖5 所示，監(jiān)督分類提取冰川信息結(jié)果中，水體被錯分為冰川（如圖5（d）），山體陰影區(qū)中的冰川未被全部識別出（如圖5（b）、5（e）），有冰磧物覆蓋物的冰川被漏分（如圖5（a）、5（c））. 由圖6 可見，比值閾值法消除了山體陰影的影響，但提取結(jié)果中“椒鹽”現(xiàn)象也較嚴重（如圖6（b）、6（e）），同時也造成了大量陰影被錯分為冰川（如圖6（a）、6（c）），水體被錯分（如圖6（d））、冰磧覆蓋區(qū)的冰川被漏分（如圖6（a）、6（c））的情況仍未得到有效解決. 如圖7 所示，利用NDSI 法提取的冰川信息結(jié)果與比值閾值法差異不大，陰影和水體被錯分為冰川（如圖7（a）、7（c）、7（d）），冰磧覆蓋區(qū)的冰川被漏分（如圖7（a）、7（c））. 由圖8 可見，在尺度選取合適時，面向?qū)ο蠓诸惙椒ǎ∣BIA）提取冰川信息的效果較好，可以區(qū)分水體與冰川（如圖8（d））. 面向?qū)ο蠓诸惙椒ㄔ诮⑿畔⒁?guī)則時，有效地結(jié)合了比值閾值法及NDSI，部分有冰磧覆蓋物的冰川被識別，但未完全正確劃分（如圖8（d））.如圖9 所示，神經(jīng)網(wǎng)絡(luò)分類的整體提取冰川信息的效果好，有效地消除了山體陰影的影響，可以區(qū)分陰影、水體與冰川，且有冰磧覆蓋物的冰川被正確劃分，但整體分類精度仍待提高.

圖4 非監(jiān)督分類提取梅里雪山地區(qū)冰川信息結(jié)果Fig. 4 Glacier information extraction of Meili Snow Mountain based on unsupervised classification

圖5 監(jiān)督分類提取梅里雪山地區(qū)冰川信息結(jié)果Fig. 5 Glacier information extraction of Meili Snow Mountain based on supervised classification

圖6 比值閾值法提取梅里雪山地區(qū)冰川信息結(jié)果Fig. 6 Glacier information extraction of Meili Snow Mountain based on ratio threshold method

圖7 NDSI 法提取梅里雪山地區(qū)冰川信息結(jié)果Fig. 7 Glacier information extraction of Meili Snow Mountain based on NDSI

圖8 面向?qū)ο蠓诸愄崛∶防镅┥降貐^(qū)冰川信息結(jié)果Fig. 8 Glacier information extraction of Meili Snow Mountain based on object oriented classification

圖9 神經(jīng)網(wǎng)絡(luò)分類提取梅里雪山地區(qū)冰川信息結(jié)果Fig. 9 Glacier information extraction of Meili Snow Mountain based on neural network classification

基于遙感影像的山岳冰川信息計算機解譯方法從不同的原理角度出發(fā)提取冰川信息，不同方法提取的結(jié)果各異，但每種方法都有各自的優(yōu)缺點.表4 結(jié)合已有研究成果與本文試驗結(jié)果，綜合比較了6 種分類方法優(yōu)缺點. 非監(jiān)督分類法不需要人為干涉，提取冰川信息的效率最高，但是整體提取效果較差、精度較低，裸地錯分嚴重，山體陰影區(qū)及冰磧覆蓋區(qū)中的冰川未被有效識別. 監(jiān)督分類法非冰川錯分的概率較低，并解決了裸地錯分的問題，但冰川信息提取效果受樣本質(zhì)量影響，主要問題是水體被錯分為冰川、山體陰影區(qū)中的冰川未被全部識別出、有冰磧覆蓋物的冰川被漏分. 比值閾值法、雪蓋指數(shù)法（NDSI）均有效地消除了山體陰影的影響，但無法區(qū)分水體與冰川，且有冰磧覆蓋物的冰川被漏分的問題仍未得到有效地解決. 面向?qū)ο蠓诸惙ǎ∣BIA）可以有效地區(qū)分水體與冰川，精度相對較高，但操作相對復雜，精度依賴于知識規(guī)則的建立. 本文選取的6 種提取方法中，神經(jīng)網(wǎng)絡(luò)分類法是精度最高的一種方法，但是由于樣本數(shù)量較少、訓練深度學習模型的迭代次數(shù)較少，整體的冰川精度未達到理想水平. 如何有效地結(jié)合傳統(tǒng)的計算機自動解譯方法，提高深度學習冰川信息提取精度，利用有效的標簽對大量數(shù)據(jù)進行深度挖掘分析且減輕神經(jīng)網(wǎng)絡(luò)對于數(shù)據(jù)的依賴程度，需要進一步的研究.

表4 梅里雪山地區(qū)冰川分類方法比較Tab. 4 Comparison of glacier classification methods in Meili Snow Mountain

3.2 目視修正為了得到更精準的梅里雪山地區(qū)冰川信息，本文結(jié)合Google Earth 數(shù)據(jù)和DEM 數(shù)據(jù)，在神經(jīng)網(wǎng)絡(luò)分類法提取冰川信息的基礎(chǔ)上，對其進行了二次目視解譯，最終得到梅里雪山地區(qū)1989—2019 年較為精確的4 期冰川信息提取成果.利用DEM 數(shù)據(jù)，通過水文分析得到梅里雪山地區(qū)的山脊線對冰川進行分割，得到了55 條獨立的冰川，并通過冰川制圖得到梅里雪山1989 年—2019年冰川邊界變化圖（圖10）. 通過分析得出，1989 年梅里雪山地區(qū)冰川面積達139.58 km2，2019 年冰川面積退縮至115.81 km2. 近30 年來，梅里雪山地區(qū)冰川面積共退縮23.77 km2，年均退縮0.79 km2，面積相對退縮率為17.03%，年均相對退縮率為0.57%（表5）. 計算公式如下：

圖10 1989—2019 年梅里雪山地區(qū)冰川邊界變化Fig. 10 Glacier boundary changes in Meili Snow Mountain from 1989 to 2019

表5 梅里雪山冰川面積變化統(tǒng)計Tab. 5 Statistics of glacier area change in Meili Snow Mountain

其中，RS代表面積絕對變化速率，S代表面積變化量，t代表時段年長；RS′代表面積相對變化率，A代表時段起始年份冰川面積.

4 結(jié)論

本文利用TM、OLI 影像，通過非監(jiān)督分類法、監(jiān)督分類法、比值閾值法、雪蓋指數(shù)法（NDSI）及基于多尺度分割的面向?qū)ο蠓ǎ∣BIA）及基于神經(jīng)網(wǎng)絡(luò)的冰川識別方法對梅里雪山地區(qū)冰川信息進行提取，綜合比較6 種方法的分類效果. 在此基礎(chǔ)上，對梅里雪山地區(qū)1989、1998、2009 年和2019年的冰川信息進行了提取，得出了以下結(jié)論：

（1）計算機自動解譯方法中，基于神經(jīng)網(wǎng)絡(luò)的冰川識別方法不僅消除了山體陰影的影響，而且可以正確區(qū)分陰影、水體及冰川，有冰磧覆蓋物的冰川信息也能被有效地提取，冰川信息提取精度最高. 但由于樣本數(shù)量較少、訓練深度學習模型的迭代次數(shù)較少，整體的冰川信息提取精度仍待提高.

（2）非監(jiān)督分類法提取冰川信息效果最差、精度較低，裸地錯分嚴重，山體陰影區(qū)、冰磧覆蓋區(qū)中的冰川未被有效識別. 監(jiān)督分類法解決了裸地錯分的問題，但未提取出山體陰影區(qū)、冰磧覆蓋區(qū)中的冰川，冰川信息提取效果受樣本質(zhì)量影響較大.比值閾值法和雪蓋指數(shù)法（NDSI）均有效消除了山體陰影的影響，但無法區(qū)分水體與冰川，且有冰磧覆蓋物的冰川被漏分的問題仍未得到有效地解決.面向?qū)ο蠓诸惙梢詤^(qū)分水體與冰川，精度相對較高，但操作相對復雜，精度依賴于知識規(guī)則的建立.

（3）目前基于遙感的冰川信息提取的過程中，目視解譯仍然是精度最高的方法，并且是冰磧覆蓋部分冰川信息解譯效果最好的方法，在計算機自動解譯后進行二次目視修正，仍是冰川信息解譯中最重要的部分.

（4）1989—2019 年梅里雪山地區(qū)冰川面積退縮了23.77 km2，年均退縮0.79 km2，面積相對退縮率為17.03%，年均相對退縮率為0.57%.

本文在解譯冰川信息的過程中發(fā)現(xiàn)，目前傳統(tǒng)的計算機自動解譯方法提取的冰川信息效果較差.而目視解譯的效率較低，并且受解譯人員的經(jīng)驗和主觀因素的影響，不同人員的解譯會有一定的偏差. 在深度學習迅速發(fā)展的背景下，如何有效結(jié)合傳統(tǒng)的計算機自動解譯方法，提高深度學習冰川信息提取精度，利用有效的標簽對大量數(shù)據(jù)進行深度挖掘分析且減輕神經(jīng)網(wǎng)絡(luò)對于數(shù)據(jù)的依賴程度是未來研究的方向.