基于波段運(yùn)算最優(yōu)方法的冰川遙感信息提取及其變化監(jiān)測

屈 冉，申文明，王昌佐，劉曉龍，柏延臣，高彥華

(1.環(huán)境保護(hù)部衛(wèi)星環(huán)境應(yīng)用中心，北京100094;2.北京師范大學(xué)地理學(xué)與遙感科學(xué)學(xué)院，北京100875)

0 引言

冰川是氣候的產(chǎn)物，又是全球氣候變化靈敏的指示器，對氣候變暖高度敏感［1-3］。對冰川變化的監(jiān)測，是了解全球氣候變化的重要手段之一。隨著全球氣候變暖，冰川普遍處于退縮狀態(tài)［4-6］。冰川退縮將直接影響區(qū)域水資源平衡與經(jīng)濟(jì)可持續(xù)發(fā)展［7-8］。因此，冰川變化已經(jīng)受到國內(nèi)外學(xué)者的眾多關(guān)注。隨著遙感技術(shù)的發(fā)展和遙感數(shù)據(jù)的豐富，遙感技術(shù)在冰川監(jiān)測中發(fā)揮著越來越重要的角色。能夠精確提取冰川的輪廓對于計(jì)算冰川的面積和變化至關(guān)重要。然而，受到冰川所處環(huán)境和冰川地區(qū)下墊面的復(fù)雜性影響，基于遙感的冰川面積提取在不同的方法下，提取的結(jié)果具有很大的差異。

當(dāng)前基于遙感影像的冰川面積提取方法主要分為人工目視解譯方法［9-11］、分類器分類方法和波段運(yùn)算方法［12］3類。冰川目視解譯方法主要在遙感應(yīng)用的早期，由于受到遙感影像數(shù)據(jù)的限制(全色航片或者低分辨率影像)，其他方法很難應(yīng)用于冰川提取。這種方法優(yōu)點(diǎn)是解譯結(jié)果可以實(shí)時(shí)受到人的限制，隨時(shí)進(jìn)行修正;缺點(diǎn)是效率低，勞動(dòng)輸入較大，解譯結(jié)果受到人的主觀判斷的限制。但是隨著遙感傳感器的發(fā)展，影像的空間分辨率不斷提高，目視解譯方法可以應(yīng)用于高空間分辨率影像，或者結(jié)合分類器分類法、波段運(yùn)算法［12］，以提高精度。分類器分類方法在冰川分類中也有較多的應(yīng)用。該分類方法受到影像內(nèi)地物復(fù)雜度的限制，影像內(nèi)部越復(fù)雜，可分性越差。例如，中低緯度地區(qū)在影像獲取時(shí)段有陰影，使用分類器分類方法很難將陰影內(nèi)冰川提取出來。但是分類器可以將處于不同狀態(tài)的冰川分為不同的類別，比如冰舌位置冰川位于雪線以下，冰川內(nèi)部含有融水或以冰的狀態(tài)存在，而雪線以上的部分，冰川往往有大量的純凈雪存在，二者在全色波段上存在顏色差別。分類器分類可以界定純凈雪與冰川冰的界限，這種方法已經(jīng)被應(yīng)用于雪線高程的提取。

波段運(yùn)算方法包括比值法［13-14］、差值法和 NDSI等。其中比值法和NDSI方法是這些方法中應(yīng)用最多的。比值法的基本原理就是利用冰雪在近紅外(NIR)和短波紅外(SWIR)波段具有較低的反射率而在可見光波段具有較高的反射率的特點(diǎn)，可以利用冰雪與其他地物差別較大的反射特征來識(shí)別冰雪。比如，在可見光波段，云具有和冰雪基本相同的反射率，但是在近紅外和短波紅外波段，由于云是液態(tài)水蒸汽，其溫度高于冰雪(固態(tài)水)，所以在這兩個(gè)波段內(nèi)具有比冰雪高的反射率，在該波段的影像上表現(xiàn)為比冰雪亮。根據(jù)這一特征，可以將冰川與云區(qū)別開來。正由于波段運(yùn)算方法具有區(qū)分云與冰川的能力，所以，該方法更廣泛地被采納。

本研究將基于波段運(yùn)算的冰川面積提取方法與分類方法進(jìn)行對比，以評價(jià)波段運(yùn)算方法中比較適合應(yīng)用于Landsat TM/ETM+數(shù)據(jù)的冰川面積提取方法，并以希夏邦馬峰研究區(qū)為例，對提取結(jié)果進(jìn)行對比與評價(jià)。

1 研究區(qū)概況

希夏邦馬峰冰川是佩枯錯(cuò)湖的主要水源，南部是尼泊爾首都加德滿都。研究區(qū)南北地形具有較大的差異:北邊是青藏高原，南邊是恒河盆地的邊緣，南北方向地形起伏較大;氣候受印度季風(fēng)氣候影響，夏季降水多，冬季降水少，所以，研究區(qū)內(nèi)的冰川物質(zhì)積累在南坡主要集中在夏季，北坡集中在冬季。

研究區(qū)冰川提取所使用TM/ETM+數(shù)據(jù)(表1)源自美國地質(zhì)調(diào)查局(http://glovis.usgs.gov/)和美國馬里蘭大學(xué)(http://glcfapp.glcf.umd.edu:8080/esdi/index.jsp)，數(shù)據(jù)空間分辨率為30 m。

表1 遙感影像的特點(diǎn)Tab.1 Property of remote sensing images

應(yīng)用NDSI方法，基于Landsat TM/ETM+數(shù)據(jù)，對希夏邦馬峰地區(qū)(圖1)1990—2010年間的冰川面積變化進(jìn)行了分析。

圖1 研究區(qū)示意圖Fig.1 Schematic diagram of study area

2 波段運(yùn)算最優(yōu)方法選擇

2.1 最大似然分類

由于冰川在可見光各波段均具有較高的反射率，所以在Landsat TM/ETM+假彩色合成影像(5，4，3波段合成)中，冰川為較亮的藍(lán)色，與周圍地物具有較大的差別。因此，分別選取假彩色合成影像中冰川、湖泊、植被陰影和裸地作為訓(xùn)練樣本進(jìn)行分類。對比原始影像(圖2a)，分類結(jié)果(圖2b)表明，最大似然分類法提取冰川時(shí)，可以較好地將冰川與湖泊、裸地和陰影等分開。但是，當(dāng)冰川位于陰影內(nèi)部時(shí)，分類器分類方法很難將冰川提取出來。

2.2 比值運(yùn)算

比值運(yùn)算法計(jì)算多光譜影像可見光近紅外波段與短波紅外波段地物反射率的比值(r)，也稱像元值，計(jì)算公式為:r=Nd4/Nd5。式中:Nd4，Nd5分別表示TM/ETM+影像第4波段與第5波段的DN值，計(jì)算結(jié)果為ratio灰度值圖像，像元值r＞0。在利用ratio灰度圖進(jìn)行冰川分類時(shí)，設(shè)置合適的閾值，將大于該閾值的部分定義為冰川，小于該值的部分定義為非冰川，最后得到冰川分布二值圖2c。閾值的設(shè)置根據(jù)影像的不同而有所差別，比如，源自不同傳感器的影像，在計(jì)算時(shí)閾值有差別;同一地區(qū)相同傳感器在不同的時(shí)相和不同的天氣條件下獲得影像也會(huì)影響該閾值的大小。但是總的來說，差別不大，一般在2.5左右。本研究用閾值2.5對冰川與非冰川進(jìn)行了分割，其結(jié)果見圖2c。

圖2 冰川提取各算法比較Fig.2 Comparison of various glacier extraction algorithms

2.3 歸一化雪覆蓋指數(shù)

NDSI(normalized difference snow index)，即歸一化雪覆蓋指數(shù)。與歸一化植被指數(shù)相似，利用冰雪在可見光波段具有強(qiáng)反射率而吸收短波紅外輻射這一差異區(qū)分冰雪與周圍地物類別。計(jì)算公式:Inds=(Nd2-Nd5)/(Nd2+Nd5)。式中:Nd2，Nd5分別表示TM/ETM+影像第2波段與第5波段的DN值，經(jīng)該式計(jì)算可以獲得NDSI灰度圖，取值范圍為［-1，1］?；叶葓D中，高值部分為冰川(圖3)。利用NDSI進(jìn)行冰川分類時(shí)，設(shè)置合適的閾值可以取得較為準(zhǔn)確的結(jié)果。閾值的選取同比值法一樣，也要考慮影像獲取時(shí)的成像條件和數(shù)據(jù)源。本研究取0.4作為NDSI冰川分類的閾值，NDSI灰度圖中，大于該閾值的區(qū)域被劃分為冰川，小于該值的部分視為非冰川。分類結(jié)果見圖2d。

圖3 NDSI灰度圖Fig.3 Grayscale image of NDSI

2.4 不同方法精度比較和最優(yōu)方法選擇

以目視解譯為真實(shí)值，評估分析其他方法的精度。通過對比3種冰川提取方法得到的二值圖冰川面積，可以比較3種方法的差異。首先，將二值圖進(jìn)行中值濾波。波段運(yùn)算的分類結(jié)果中含有噪聲，這些噪聲源自冰磧物和非常小的冰川［14］。將二值圖進(jìn)行窗口大小為5的中值濾波，以消除噪聲的影響(圖4);然后，將濾波后的二值圖轉(zhuǎn)化為矢量圖，轉(zhuǎn)化為矢量圖的目的是便于編輯和計(jì)算冰川面積。上述分類方法的分類結(jié)果中，結(jié)冰的冰川湖和冰磧物常常被誤分類為冰川，所以，需要結(jié)合假彩色影像(TM5，4，3波段合成)，將矢量圖內(nèi)誤分部分進(jìn)行目視辨別和手動(dòng)修改。最后，利用修改后的矢量圖進(jìn)行冰川面積計(jì)算。

對比3種分類方法下得到的二值圖冰川面積(表2)可以看出，最大似然法由于難以提取陰影內(nèi)部的冰川，導(dǎo)致提取的冰川面積偏低。波段運(yùn)算法的比值法、NDSI法與目視解譯法所提取的冰川面積接近。

表2 不同方法提取的冰川面積 km2Tab.2 Areas of glacier extracted by means of various methods

圖4 濾波前后對比圖Fig.4 Comparison of non-filtered images with filtered images

3 基于NDSI的冰川變化監(jiān)測

對冰川提取結(jié)果(圖5)及冰川變化(圖6)進(jìn)行統(tǒng)計(jì)分析，結(jié)果表明，1990，2000，2010年研究區(qū)冰川總面積分別是 985.1，860.4，740.1 km2，其中北坡分別是252.7，235.3，205.1 km2，南坡分別是 732.4，625.1，535.0 km2(圖7)，南坡冰川面積大于北坡;1990—2000，2000—2010 年間冰川退縮面積分別是124.8，120.3 km2，其中南坡分別退縮107.3，90.1 km2，北坡分別退縮17.5，30.2 km2，南坡退縮面積多于北坡，南北坡冰川面積均呈退縮趨勢，這與之前的研究結(jié)果［14］一致;兩個(gè)時(shí)期內(nèi)，南坡冰川的退縮率分別是14.7%，14.4%，北坡冰川的退縮率分別是6.9%，12.8%(圖8)，但北坡冰川退縮率有增加的趨勢，南坡相反。

圖5 冰川提取結(jié)果Fig.5 Result of glacier extraction

圖6 研究區(qū)1990—2000，2000—2010年冰川變化Fig.6 Figures of glacier change in 1990—2000 and 2000—2010

圖7 研究區(qū)冰川面積變化Fig.7 Variation of glacier area in the study area

圖8 研究區(qū)冰川變化率Fig.8 Variation rate of glacier in the study area

4 結(jié)論與討論

對比分析表明，分類方法最大似然分類法提取冰川時(shí)，可以較好地將冰川與湖泊、裸地和陰影等分開。但是，當(dāng)冰川位于陰影內(nèi)部時(shí)，分類器分類方法很難將冰川提取出來。比值運(yùn)算法計(jì)算地物多光譜影像近紅外波段與短波紅外波段反射值的比值，在利用灰度圖進(jìn)行冰川提取時(shí)，設(shè)置合適的閾值，將大于該閾值的部分定義為冰川，小于該值的部分定義為非冰川，最后得到冰川分布二值圖。閾值的設(shè)置根據(jù)影像的不同而有所差別，同一地區(qū)相同傳感器在不同的時(shí)相和不同的天氣條件下獲得影像也會(huì)影響閾值的大小?；贜DSI波段運(yùn)算方法具有運(yùn)算速度快、結(jié)果準(zhǔn)確和算法簡單等優(yōu)勢。

基于NDSI方法提取的希夏邦馬峰地區(qū)的冰川變化結(jié)果表明，1990—2010年間該地區(qū)冰川面積在不斷退縮，并且有加劇的趨勢，這與之前關(guān)于該地區(qū)冰川研究的結(jié)論［15-16］具有一致性。

本研究中基于NDSI方法的冰川面積提取雖然取得了較好的結(jié)果，但是，卻不適合直接用于提取研究區(qū)內(nèi)冰磧物。冰磧物具有更為復(fù)雜的波譜特征，在可見光波段，具有與周圍地標(biāo)相似的波普反射特征。但是由于冰磧物下面有冰覆蓋，所以，在近紅外和短波紅外波段波譜反射強(qiáng)度低于周圍地表。利用這一特征將可以實(shí)現(xiàn)冰磧物的自動(dòng)提取。

［1］ IPCC.Climate Change 2007:The Scientific Basis［M］//Solomon S D，Qin D，Manning M，et al.Contribution of Working GroupⅠto the Fourth Assessment Report of the Intergovernmental Panel on Climate Change.Cambridge:Cambridge University Press，2007:1-996.

［2］ Bolch T，Buchroithner M，Pieczonka T，et al.Planimetric and Volumetric Glacier Changes in the Khumbu Himal，Nepal，since 1962 Using Corona，Landsat TM and AS-TER Data［J］.Journal of Glaciology，2008，54:592-600.

［3］ Haeberli W，Hoelzle M，Paul F，et al.Integrated Monitoring of Mountain Glaciers as Key Indicators of Global Climate Change:The European Alps［J］.Annals of Glaciology，2007，46(1):150-160.

［4］ Dyurgerov M B，Meier M F.Twentieth Century Climate Change:Evidence from Small Glaciers［J］.Proceedings of the National Academy of Sciences of the United States of America，2000，97(4):1406-1411.

［5］ 陳璇，覃建雄，李嫻.海螺溝冰川地質(zhì)特征與冰川退縮防治［J］.地域研究與開發(fā)，2010，29(6):85-88.

［6］ 葉慶華，陳鋒，姚檀棟，等.近30年來喜馬拉雅山脈西段納木那尼峰地區(qū)冰川變化的遙感監(jiān)測研究［J］.遙感學(xué)報(bào)，2007，11(4):511-520.

［7］ Huggel C，K??b A，Haeberli W，et al.Remote Sensing Based Assessment of Hazards from Glacier Lake Outbursts:A Case Study in the Swiss Alps［J］.Canadian Geotechnical Journal，2002，39:316-330.

［8］ 李開明，李忠勤，高聞?dòng)?，?近期新疆東天山冰川退縮及其對水資源影響［J］.科學(xué)通報(bào)，2011，56(32):2708-2716.

［9］ 許君利，劉時(shí)銀，張世強(qiáng)，等.塔里木盆地南緣喀拉米蘭河克里雅河流內(nèi)流區(qū)近30 a來的冰川變化研究［J］.冰川凍土，2006，28(3):312-319.

［10］ 王葉堂，侯書貴，魯安新，等.近40年來天山東段冰川變化及其對氣候的響應(yīng)［J］.干旱區(qū)地理，2008，31(6):813-820.

［11］ 張明杰，秦翔，杜文濤，等.1957—2009年祁連山老虎溝流域冰川變化遙感研究［J］.干旱區(qū)資源與環(huán)境，2013，27(4):70-75.

［12］ 田洪陣，楊太保，劉沁萍.1976—2010年祁連山中段崗格爾肖合力雪山冰川退縮和氣候變化的關(guān)系研究［J］.干旱區(qū)資源與環(huán)境，2012，26(7):41-46.

［13］ Paul F，Huggel C，K??b A，et al.Comparison of TM-derived Glacier Areas with Higher Resolution Data Sets［J］.EARSeL eProceedings，2002(2):15-21.

［14］ Pellikka P，Rees W G.Remote Sensing of Glaciers:Techniques for Topographic，Spatial and Thematic Mapping of Glaciers［M］.London:Taylor＆ Francis Group，2010:147-148.

［15］ 車濤，李新，Mool P K，等.希夏邦馬峰東坡冰川與冰川湖泊變化遙感監(jiān)測［J］.冰川凍土，2005，27(6):801-805.

［16］ 馬凌龍，田立德，蒲健辰，等.喜馬拉雅山中段抗物熱冰川的面積和冰儲(chǔ)量變化［J］.科學(xué)通報(bào)，2010，55(18):1766-1774.