楊 佳， 薛莎莎， 蘇永恒， 任慶福

（1.航天宏圖信息技術(shù)股份有限公司，北京 100195；2.上海宏圖空間網(wǎng)絡(luò)科技有限公司，上海 201702）

0 引言

全球氣候變化是當(dāng)今世界各國關(guān)注的焦點(diǎn)問題，世界氣象組織（World Meteorological Organization，WMO）在《2020年全球氣候狀況報告》中指出，2020年全球氣候變化進(jìn)一步惡化，隨著全球氣溫的升高、冰川融化退縮、海平面上升，對社會經(jīng)濟(jì)發(fā)展及生態(tài)系統(tǒng)造成嚴(yán)重影響［1］。冰川作為其重要指標(biāo)器，是全球氣候變化的研究熱點(diǎn)［2］。冰川消融不僅會造成海平面上升，而且會在內(nèi)陸引起冰崩、冰湖潰決等災(zāi)害，這對區(qū)域生態(tài)環(huán)境以及人類社會生活均會造成巨大影響［3］，因而實(shí)時監(jiān)測冰川范圍變化具有重要意義［4］。冰川范圍監(jiān)測的傳統(tǒng)研究方法需要進(jìn)行野外考察，外業(yè)調(diào)查的方法對冰川局部范圍變化的精確分析具有一定優(yōu)勢，但存在監(jiān)測范圍較小、耗時費(fèi)力等問題。遙感數(shù)據(jù)具有時間連續(xù)，空間覆蓋廣，數(shù)據(jù)獲取方便等特點(diǎn)，基于遙感數(shù)據(jù)監(jiān)測冰川一定程度上彌補(bǔ)了傳統(tǒng)方法的不足［5］。

基于光學(xué)影像的常規(guī)遙感監(jiān)測方法主要包括冰雪指數(shù)閾值法［6］以及特征分類法［7］。冰雪指數(shù)閾值法基于冰雪在可見光波段的高反射率和近紅外波段的低反射率，進(jìn)行波段組合以及閾值設(shè)定獲取冰川范圍信息，具有簡易、快速、相對準(zhǔn)確的優(yōu)勢，在冰川監(jiān)測中得到廣泛應(yīng)用［8］。Huang等［9］使用“多時態(tài)最小NDSI復(fù)合法”實(shí)現(xiàn)了清潔冰川和非季節(jié)性積雪的提取，對于NDSI指數(shù)閾值法無法直接區(qū)分冰川與冰湖的問題，使用湖泊邊界數(shù)據(jù)集對研究區(qū)影像進(jìn)行湖泊掩膜處理，從而消除了冰湖對提取結(jié)果的影響；季漩等［10］提出增強(qiáng)型雪被指數(shù)（ENDSI）提取了崗加曲巴冰川山體陰影區(qū)的冰川范圍；劉凱等［11］通過比值閾值法和目視解譯法提取冰川邊界，分析了1993—2015年喀喇昆侖山努布拉流域的冰川變化特征。以上冰雪指數(shù)閾值法的應(yīng)用，從結(jié)果而言均實(shí)現(xiàn)了冰川范圍提取，但未能從根本上解決冰雪指數(shù)閾值法無法直接區(qū)分冰川與冰湖的問題。特征分類法主要包括監(jiān)督分類、非監(jiān)督分類、面向?qū)ο蠓诸?、機(jī)器學(xué)習(xí)分類。Aniya等［12］通過非監(jiān)督分類方法提取Patagonia冰原的冰川信息；彥立利等［13］基于面向?qū)ο蟮姆诸惙椒▽Μ敿{斯河上游區(qū)域的冰川及常年積雪進(jìn)行邊界提?。煌踔椅涞龋?4］提出一種上下文感知深度學(xué)習(xí)語義分割網(wǎng)絡(luò)法消除了冰湖及山體陰影的影響完成了冰川范圍提取。以上特征分類法的應(yīng)用，同樣實(shí)現(xiàn)了冰川范圍提取，但仍普遍存在同譜異物現(xiàn)象。為解決以上常規(guī)遙感監(jiān)測方法存在的問題，本文以Landsat-8數(shù)據(jù)為基礎(chǔ)，根據(jù)冰雪的光敏波段以及OLI（operational land imager）傳感器的波段特征構(gòu)建可區(qū)分冰川、冰湖的改進(jìn)冰雪指數(shù)，然后在此基礎(chǔ)上將改進(jìn)冰雪指數(shù)融入面向?qū)ο蠓诸惙ㄖ刑岢雒嫦驅(qū)ο?改進(jìn)冰雪指數(shù)法。以各拉丹冬冰川作為試驗(yàn)區(qū)，運(yùn)用此方法識別冰川邊界，并將此區(qū)域的青藏高原冰川數(shù)據(jù)產(chǎn)品及常規(guī)遙感監(jiān)測方法提取結(jié)果作為參考數(shù)據(jù)，驗(yàn)證方法的有效性和穩(wěn)健性。

1 研究區(qū)與數(shù)據(jù)

各拉丹冬冰川位于青海省格爾木市唐古拉山鎮(zhèn)，屬于內(nèi)陸高寒地區(qū)，地理位置為33°00′～33°40′N、90°45′～91°20′E［15］，東西寬約30 km，南北長約50 km。該地區(qū)常年氣溫在0℃以下，年降水量僅200 mm，且90%集中在5—9月［16］。唐古拉山境內(nèi)有現(xiàn)代冰川1 530條，總面積約2 313 km［17］，各拉丹冬冰川群是其中最大的冰川群。各拉丹冬冰川群中有三個較大冰川群［18］（圖1），其冰川融水是長江源頭重要的補(bǔ)給來源，本研究以這三個冰川群為研究對象。

圖1 研究區(qū)位置及冰川分布Fig.1 Location and glacier distribution of the study area

本研究的基礎(chǔ)數(shù)據(jù)為Landsat-8 OLI數(shù)據(jù)、青藏高原2017年冰川數(shù)據(jù)產(chǎn)品——TPG2017。它來自國家青藏高原科學(xué)數(shù)據(jù)中心（http：//data.tpdc.ac.cn/zh-hans/），此數(shù)據(jù)集是基于2017年的Landsat-8 OLI數(shù)據(jù)，通過目視解譯法獲得的冰川范圍數(shù)據(jù)產(chǎn)品。

Landsat-8數(shù)據(jù)來自USGS（https：//earthexplorer.usgs.gov/），Landsat-8 OLI傳感器與Landsat-7 ETM+傳感器的波段相比新增了海藍(lán)波段（433～453 nm）以及短波紅外波段（2 100～2 300 nm）。Landsat-8影像基于以下幾個方面進(jìn)行篩選：①影像時間各拉丹冬冰川為夏季積累型冰川，影像獲取時間應(yīng)處于或接近于冰川的消融末期；②影像質(zhì)量研究區(qū)夏秋兩季降水較多，云遮擋無法避免，篩選云量較少的影像；③山體陰影影響太陽高度角大較時，山體陰影較小，因而盡量選擇太陽高度角大的影像。綜合以上，選取2017年10月7日的一景Landsat-8數(shù)據(jù)（條帶號為138/37）作為試驗(yàn)數(shù)據(jù)，該影像云量僅為0.25%，太陽高度角達(dá)47.5°影像質(zhì)量較好，滿足研究需求。

2 原理與方法

2.1 光譜分析

遙感影像中的冰川信息主要是通過冰川表層積雪的光譜特性及其差異來反映，不同光譜通道獲取到的冰川信息與冰川表層積雪的特征狀態(tài)（融雪狀態(tài)、雪粒徑）密切相關(guān)。不同融雪狀態(tài)下冰雪的光譜曲線［圖2（a）］存在明顯差異，隨著融雪程度的加深，冰雪中含水量增多，冰雪反射率不斷降低，光譜曲線走勢也趨于平緩；不同粒徑冰雪的光譜曲線［圖2（b）］隨著雪粒徑增大，反射率降低。

圖2 冰雪光譜曲線Fig.2 Spectral curve of ice and snow：different snow melting conditions（a）and different particle sizes（b）

冰雪在不同的特征狀態(tài)下，光譜曲線存在不同，但其“峰—谷”總體變化趨勢表現(xiàn)一致，同時具有明顯的共性特征：①隨著波長的增加，冰雪反射率整體呈下降趨勢，在藍(lán)光波段出現(xiàn)峰值，近紅外長波波段出現(xiàn)谷底；②在可見光波段（390～780 nm），冰雪反射率整體較高且相對穩(wěn)定，在藍(lán)光波段出現(xiàn)峰值，且隨波長的增長逐漸降低；③在近紅外短波波段（780～1 100 nm），冰雪反射率逐漸降低；④在近紅外長波波段（1 100～2 526 nm），冰雪反射率急劇下降，尤其在水的吸收帶（1 400 nm、1 940 nm）反射率降至谷底。

由冰雪光譜曲線的共性特征可知冰川在藍(lán)光波段具有強(qiáng)反射特性，近紅外長波波段具有強(qiáng)吸收特性。冰雪和水體在近紅外長波波段反射率差別較大，特別是在1 100～2 500 nm波段，I類水體反射率很低幾乎趨近于0；由此可基于藍(lán)光波段和近紅外長波波段區(qū)分冰川、冰湖。

2.2 Landsat-8冰川識別指標(biāo)

冰川遙感識別指標(biāo)多為綠波段和短波紅外波段構(gòu)建的冰雪指數(shù)，此冰雪指數(shù)一定程度上可以提取冰川范圍，但是仍存在地物錯分問題，尤其是難以區(qū)分冰川和冰湖。由2.1節(jié)光譜曲線特征分析可知，冰川識別的遙感敏感波段為藍(lán)光波段和近紅外長波波段，這兩個波段可明顯區(qū)分冰川、冰湖，本研究依此構(gòu)建改進(jìn)冰雪指數(shù)。

以Landsat-8為數(shù)據(jù)源，根據(jù)OLI傳感器各波段的光譜范圍，將ρCoastal、ρSWIR1分別作為藍(lán)波段、近紅外長波波段。對于冰川范圍的提取，僅用個別波段或多個波段數(shù)據(jù)的對比分析來提取冰川邊界具有一定的局限性，因此選用波段組合的方式（減、除、線性組合），獲得對冰川狀態(tài)信息具有表征意義的數(shù)值，由此提出如下三種改進(jìn)的冰雪指數(shù)（表1）。

表1 改進(jìn)的冰雪指數(shù)Table 1 The improved ice and snow indexes

2.3 特征分類

基于規(guī)則的面向?qū)ο蠓诸惙ㄊ峭ㄟ^影像分割將影像分為若干同質(zhì)性像元組成的影像對象，然后將影像對象作為分析單元，綜合分析其光譜特征、紋理特征形成對象特征知識庫，最后在此基礎(chǔ)上進(jìn)行目標(biāo)信息的提取。此方法在高海拔區(qū)域的冰川識別方面較傳統(tǒng)分類方法精度更高［19］。采用基于規(guī)則的面向?qū)ο蠓诸惙▽Ρ┲笖?shù)閾值法進(jìn)行優(yōu)化和綜合即面向?qū)ο?改進(jìn)冰雪指數(shù)法，進(jìn)一步明確冰湖和冰川的邊界，精確提取冰川范圍。

面向?qū)ο?改進(jìn)冰雪指數(shù)法技術(shù)流程（圖3）如下：對預(yù)處理后的Landsat-8數(shù)據(jù)使用三種改進(jìn)冰雪指數(shù)（RSI、NDSI*、DSI）分別提取冰川范圍，篩選出最優(yōu)的冰雪指數(shù)影像，并將其與多光譜影像進(jìn)行波段合成構(gòu)建冰川特征影像，然后結(jié)合紋理和光譜等規(guī)則，提取冰川范圍信息。

圖3 面向?qū)ο?改進(jìn)冰雪指數(shù)法技術(shù)流程Fig.3 Technical flowchart of the object-oriented-improved ice and snow index method

3 數(shù)據(jù)處理與結(jié)果分析

3.1 冰川特征影像

本研究使用的Landsat-8影像數(shù)據(jù)級別為L1T，經(jīng)過了系統(tǒng)輻射校正和地面控制點(diǎn)幾何校正，并且完成了數(shù)字高程模型（DEM）校正處理。因此在數(shù)據(jù)預(yù)處理過程中，僅需進(jìn)行輻射定標(biāo)、大氣校正等處理。本研究使用PIE-Basic6.0軟件對Landsat-8 OLI影像進(jìn)行預(yù)處理，消除大氣和光照等因素對地物反射的影響，提高影像質(zhì)量。預(yù)處理后，為了方便展開對各拉丹冬冰川范圍的識別，對影像進(jìn)行了裁剪，保留研究區(qū)內(nèi)較大的三條冰川及其周邊區(qū)域，然后進(jìn)行冰雪指數(shù)閾值提取，并將提取結(jié)果與原始影像疊加對比，通過人機(jī)交互的方式確定三種冰雪指數(shù)的閾值取值范圍。根據(jù)以往研究經(jīng)驗(yàn)，在各指數(shù)閾值范圍內(nèi)，以0.05為間隔，目視解譯原始影像冰川范圍與冰雪指數(shù)提取結(jié)果的重疊度，獲得各冰雪指數(shù)最優(yōu)提取結(jié)果（圖4）。

圖4 冰雪指數(shù)提取結(jié)果Fig.4 Extraction results based on ice and snow indexes：original image（a），extraction results of RSI（b），extraction results of NDSI*（c）and extraction results of DSI（d）

對比分析提取結(jié)果發(fā)現(xiàn)，三種指數(shù)均可很好地區(qū)分云和冰川，但是在冰川與冰湖的識別上，DSI提取結(jié)果明顯優(yōu)于其他。為了進(jìn)一步評價提取結(jié)果精度，將三種冰雪指數(shù)提取結(jié)果分別與TPG2017進(jìn)行對比，以漏分誤差、總體精度為指標(biāo)，篩選最優(yōu)冰雪指數(shù)影像。各冰雪指數(shù)提取結(jié)果的精度如表2所示。

表2 冰雪指數(shù)提取結(jié)果的精度Table 2 Accuracy of extraction results based on ice and snow indexes

綜合以上，DSI冰雪指數(shù)不但可區(qū)分冰川和云以及冰湖，而且提取結(jié)果總體精度最高，因而選擇其作為最優(yōu)冰雪指數(shù)提取結(jié)果。對經(jīng)大氣校正后的Landsat-8數(shù)據(jù)使用Band4、Band3、Band2進(jìn)行波段合成處理形成多光譜影像，然后將其與DSI冰雪指數(shù)影像再次進(jìn)行波段合成，獲得冰川特征影像。

3.2 冰川提取

采用面向?qū)ο?改進(jìn)冰雪指數(shù)法，將DSI冰雪指數(shù)影像和多光譜影像進(jìn)行波段合成構(gòu)建特征影像，根據(jù)特征影像的波段特征，劃分到不同的規(guī)則下并建立規(guī)則屬性，然后各自設(shè)置合適的閾值范圍，獲得最優(yōu)的冰川提取結(jié)果。面向?qū)ο?改進(jìn)冰雪指數(shù)冰川提取的過程如下。

3.2.1 創(chuàng)建影像對象

創(chuàng)建影像對象的過程包括影像分割和合并。冰川在遙感影像中與其背景環(huán)境差異顯著，因而采用邊緣檢測分割算法，基于影像中相鄰區(qū)域像元的光譜、紋理特征，調(diào)整分割參數(shù)實(shí)現(xiàn)不同尺度影像對象邊界的分割，從而產(chǎn)生由細(xì)到粗的多尺度分割結(jié)果。影像分割時，一個特征對象可能會被分成很多個部分，合并算法可解決這個問題，冰川像元具有連續(xù)性、紋理性較強(qiáng)的特征，所以采用Full Lambda Schedule算法［20］進(jìn)行合并。分割合并參數(shù)閾值的設(shè)定，對冰川對象內(nèi)部的差異性以及冰川對象與其他地物類別之間的異質(zhì)性均存在較大影響［21］，如表3所示進(jìn)行多次試驗(yàn)，結(jié)果如圖5所示。當(dāng)冰川對象邊界清晰，影像對象不破碎且不包括其他類別地物時，確定最優(yōu)組合分割合并閾值組合（70，98）。

圖5 影像分割結(jié)果Fig.5 Image segmentation results：Combination 1（a），Combination 2（b），Combination 3（c）and Combination 4（d）

表3 分割合并算法的閾值組合Table 3 Threshold combination of segmentation and merge algorithm

3.2.2 特征提取

根據(jù)特征影像的波段特征，在不同的規(guī)則下建立規(guī)則屬性［22］，多光譜影像數(shù)據(jù)集中的藍(lán)波段在Texture規(guī)則下構(gòu)建Texture Range屬性，DSI冰雪指數(shù)波段在Spectral規(guī)則下構(gòu)建Spectral Mean規(guī)則屬性，然后在此基礎(chǔ)上各自設(shè)置合適的閾值范圍（表4），從而獲得最優(yōu)的冰川提取結(jié)果。

表4 屬性閾值范圍Table 4 Attribute threshold range

3.2.3 結(jié)果對比分析

將TPG2017冰川數(shù)據(jù)產(chǎn)品作為參考，對比分析面向?qū)ο?改進(jìn)冰雪指數(shù)法的冰川提取結(jié)果。結(jié)果表明，面向?qū)ο?改進(jìn)冰雪指數(shù)法的總體精度高達(dá)97.26%，漏分誤差為2.74%（表5），將二者疊加在原始影像上［圖6（a）］，可直觀看到在冰川邊界細(xì)節(jié)的識別上，面向?qū)ο?改進(jìn)冰雪指數(shù)法提取結(jié)果更為精準(zhǔn)。

表5 冰川提取結(jié)果的精度Table 5 Accuracy of glacier extraction results

圖6 冰川提取結(jié)果對比Fig.6 Comparison of glacier extraction results between the object-oriented-improved ice and snow index method with TPG2017（a）and conventional remote sensing methods（b）

將面向?qū)ο?改進(jìn)冰雪指數(shù)法冰川提取結(jié)果與NDSI閾值法以及面向?qū)ο蠓诸惙ǖ谋ㄌ崛〗Y(jié)果相對比［圖6（b）］，與NDSI閾值法相比，面向?qū)ο?改進(jìn)冰雪指數(shù)法的提取結(jié)果，很好地區(qū)分了冰川與冰湖，冰川邊界較為平滑，邊界細(xì)節(jié)部分的識別也更為準(zhǔn)確，零星分布的碎片化區(qū)域更少，冰川范圍提取結(jié)果的總體精度提高了0.12%；與面向?qū)ο蠓诸惐ㄌ崛〗Y(jié)果相比，面向?qū)ο?改進(jìn)冰雪指數(shù)法可以更好地區(qū)分云與冰川，總體精度提高了0.38%。

面向?qū)ο?改進(jìn)冰雪指數(shù)法綜合了改進(jìn)冰雪指數(shù)和面向?qū)ο蠓诸惙ǖ膬?yōu)點(diǎn)，同時也彌補(bǔ)了傳統(tǒng)遙感監(jiān)測方法的不足，與傳統(tǒng)冰雪指數(shù)閾值法相比，此方法精準(zhǔn)地識別了冰川與冰湖的邊界，并解決冰雪指數(shù)閾值法難以識別的“胡椒鹽效應(yīng)”，實(shí)現(xiàn)冰川邊界的精確提??；與面向?qū)ο蠓诸惖姆椒ㄏ啾?，此方法一定程度上解決了面向?qū)ο蠓诸惙ㄍV異物的問題。

4 結(jié)論

針對傳統(tǒng)遙感監(jiān)測方法存在的地物錯分、同譜異物問題，本研究分析冰雪光譜特征提出改進(jìn)冰雪指數(shù)，并以Landsat-8影像為基礎(chǔ)數(shù)據(jù)，采用面向?qū)ο?改進(jìn)冰雪指數(shù)法提取各拉丹冬冰川范圍，然后將提取結(jié)果分別與青藏高原2017年冰川數(shù)據(jù)產(chǎn)品、NDSI閾值法提取結(jié)果、面向?qū)ο蠓诸惙ㄌ崛〗Y(jié)果進(jìn)行對比分析，結(jié)論如下：

（1）冰雪在藍(lán)光波段具有強(qiáng)反射特性，在近紅外長波波段具有強(qiáng)吸收特性，這兩個波段為冰川識別的敏感波段，其不僅對冰川具有較高的辨識度，而且在區(qū)分冰川與冰湖方面也有明顯效果。分析Landsat-8 OLI傳感器的波段特征，其冰川敏感波段為Coastal波段、SWIR1波段。

（2）針對Landsat-8數(shù)據(jù)的冰川敏感波段提出了三種改進(jìn)冰雪指數(shù)RSI、NDSI*、DSI。以各拉丹冬冰川為試驗(yàn)區(qū)，基于三種改進(jìn)冰雪指數(shù)提取冰川范圍，并將提取結(jié)果分別與青藏高原2017年冰川數(shù)據(jù)產(chǎn)品進(jìn)行對比分析，結(jié)果表明三種冰雪指數(shù)均很好地提取了冰川邊界，提取精度均達(dá)到95%以上，在易混淆地物識別方面，三種指數(shù)均可區(qū)分云和冰川，但是在冰川與冰湖的區(qū)別上，DSI冰雪指數(shù)明顯優(yōu)于其他。

（3）采用面向?qū)ο?改進(jìn)冰雪指數(shù)法提取各拉丹冬冰川范圍，提取結(jié)果的總體精度高達(dá)97.26%，與傳統(tǒng)NDSI閾值法、面向?qū)ο蠓诸惙ㄌ崛〗Y(jié)果相比精度均獲得一定的提高。此方法綜合了冰雪指數(shù)閾值提取法和面向?qū)ο蠓诸惙ǖ膬?yōu)點(diǎn)，彌補(bǔ)了二者的缺陷，更為精確地識別了冰川邊界。

（4）面向?qū)ο?改進(jìn)冰雪指數(shù)法一定程度上解決了冰川提取過程中同譜異物的問題，特別是對于冰湖與冰川的區(qū)分。但仍存在一些不足，本研究原始數(shù)據(jù)篩選的是冰川消融末期且太陽高度角較高的影像，對冰川區(qū)和非冰川區(qū)積雪難以區(qū)分的現(xiàn)象以及陰影區(qū)冰川難識別的問題未作詳細(xì)考慮。后續(xù)，針對以上不足，擬采用多源遙感數(shù)據(jù)（高分辨率影像、DEM等），融合多元冰川特征指標(biāo)，進(jìn)一步完善面向?qū)ο?改進(jìn)冰雪指數(shù)法，力求從根本上解決潔凈冰川分類的難點(diǎn)問題。