都偉冰，李均力，包安明，王寶山，王雙亭

1.河南理工大學測繪與國土信息工程學院，河南 焦作454000；2.中國科學院新疆生態(tài)與地理研究所，新疆 烏魯木齊830011

1 引 言

由于容易受到積雪及山體陰影的影響，在利用遙感技術(shù)提取冰川信息的過程中，遙感數(shù)據(jù)的時相選擇、觀測角度等會對冰川提取的準確度和精度造成較大影響。2002年全球陸地冰空間監(jiān)測 計 劃（global land ice measurements from space，GLIMS）發(fā)布了由遙感技術(shù)提取的全球陸地冰川編目數(shù)據(jù)，然而用于提取新疆山區(qū)冰川的遙感數(shù)據(jù)在時間跨度上長達30多年，冰川邊界的精度和時效性存在較大問題［1］。目前，遙感冰川信息提取方法主要有分類法和閾值分割法［2］。分類方法適用于較小局部區(qū)域精細分類，效果較好，但算法復(fù)雜，難以推廣；分割方法主要采用基于雪蓋指數(shù)或比值圖像［3］的單一閾值提取，然而針對大區(qū)域的冰川制圖，簡單的閾值分割無法適用于高山冰川區(qū)域復(fù)雜多樣的局部環(huán)境。“全域—局部”迭代的閾值分割方法可實現(xiàn)大范圍，復(fù)雜環(huán)境下地物提取閾值的快速選擇［4－5］。

新疆地區(qū)冰川分布在地形陡峭的高山地區(qū)，受太陽入射角度的影響，部分冰川往往被山體陰影遮擋。文獻［6］在提取喜馬拉雅山地區(qū)冰川時采用具有較高太陽高度角的影像來減少山體陰影的影響。另外，由于積雪與冰川很難區(qū)分，而高山區(qū)域終年降雪，很難獲取無積雪覆蓋的遙感影像［7］。因此，采用單一時相的遙感影像很難獲取精確的冰川邊界。采用一定時段內(nèi)的多時相遙感影像對冰川邊界進行觀測，綜合利用多角度信息以減少地形對冰川的影響，可獲取某個時段內(nèi)冰川的最佳邊界信息。為此，本文擬采用一定時段內(nèi)的Landsat多時相遙感影像提取天山托木爾峰西部的冰川，通過不同太陽入射角度下的冰雪信息疊加分析，消除山體陰影和積雪對冰川信息的干擾，以獲得較精確的冰川邊界信息。

2 高山冰川提取的特點

新疆天山地區(qū)的冰川屬于溫帶大陸型冰川，發(fā)育于海拔4000m以上的地形陡峭山區(qū)，其北坡由于受崎嶇地形的遮擋，冰川發(fā)育良好，不少冰川位于巨大山體形成的陰影區(qū)域內(nèi)［8］。在不同的太陽高度條件下，山體形成的陰影分布是不同的，可以通過多角度遙感的方法減少遮擋冰川的陰影面積。

通過對天山冰川變化最為劇烈地區(qū)的冰川觀測記錄發(fā)現(xiàn)，1972—2005年間哈密地區(qū)的廟爾溝冰川平均退縮量為5m／a，廟爾溝冰帽年退縮量為2.3m／a，烏魯木齊東部的博格達冰川長度平均退縮量為4.11m／a，四工河4號冰川末端2006—2009年最大退縮量為13.3m／a［9］。如果采用空間分辨率為30m的Landsat TM數(shù)據(jù)（30m）監(jiān)測冰川的變化，在冰川變化最劇烈的區(qū)域，可觀測到冰川變化的時間尺度約為2.26a。因此，可認為Landsat提取的冰川邊界在2年之內(nèi)是相對穩(wěn)定的。

為此，本文考慮采用某個時段內(nèi)的多期遙感數(shù)據(jù)聯(lián)合確定冰川的邊界，其主要思想是合理利用不同時相遙感影像成像時的不同太陽入射角度來減小山體陰影對冰川遮擋的影響，而多期遙感數(shù)據(jù)聯(lián)合觀測可以獲取冰川上的最小積雪覆蓋面積。該方法基于以下兩個假設(shè)條件：

（1）利用多時相遙感數(shù)據(jù)的多角度信息可以最大程度減小甚至消除陰影對冰川提取的影響。

（2）采用一定時段內(nèi)（2年內(nèi)）的多時相Landsat影像提取冰川邊界時，可認為這個時間段內(nèi)的冰川邊界不變。

這樣就可以用多時相遙感提取的冰川邊界進行疊置分析，從而逼近冰川的真實邊界。本文的多角度遙感冰川提取的技術(shù)路線如圖1所示，包括山體陰影提取、冰雪信息提取、多角度陰影剔除和多時相冰川邊界疊置分析。

圖1 多角度冰川提取流程Fig.1 Multi－angle glacier delineation flow

3 多角度冰川邊界提取方法

3.1 山體陰影的提取

本文基于高程數(shù)據(jù)和遙感成像的太陽入射角度，利用朗伯特余弦定律［10］建立山體陰影提取的物理模型，通過計算數(shù)字高程模型（digital elevation model，DEM）的中的坡面法線，結(jié)合太陽角度得到朗伯特余弦定律的入射角，最終計算出山體陰影圖像。本文在地形建模工具中輸入DEM數(shù)據(jù)和對應(yīng)時相的太陽角度參數(shù)，生成地形暈渲圖像，通過分析得到山體陰影。

圖2以天山西段托木爾峰東北部的幾條冰川［11］為例，通過 DEM（圖2（b））和 Landsat TM 影像（圖2（a）波段組合為7／4／2）對應(yīng)的太陽角度求得該區(qū)域的地形暈渲圖，提取出山體陰影（圖2（d）），所用Landsat TM影像成像時太陽方位角147.22°，高度角47.56°，成像日期為2011－09－13。圖2（c）是單純利用光譜信息提取陰影區(qū)域疊加在Landsat TM影像上顯示。由圖2可以看出，由高程數(shù)據(jù)結(jié)合衛(wèi)星成像時的太陽提取的山體陰影效果較光譜信息提取的陰影好。與單純利用遙感影像的光譜信息提取陰影相比，利用太陽角度結(jié)合DEM提取山體陰影能夠減少因“同物異譜”或“異物同譜”造成的錯分現(xiàn)象。

圖2 山體陰影提取Fig.2 Extraction of mountain shadow

3.2 “全域－局部”迭代冰雪信息提取方法

本文采用“全域—局部”分步迭代的閾值提取方法［12］提取遙感圖像上的冰雪信息。這種方法先用單一閾值對整景遙感圖像進行全域閾值分割，獲取初步的目標與背景信息，然后針對每個目標單元，采用分步迭代的方法對緩沖區(qū)內(nèi)的目標進行精細分割。以圖3為例，“全域－局部”迭代法提取冰雪的技術(shù)流程如下：

（1）歸一化雪覆蓋指數(shù)NDSI的提?。?3］。冰雪在短波紅外波段有強吸收特性，在綠波段有強的反射特性，能很好地區(qū)分冰雪與背景信息［14］。以Landsat TM影像為例，歸一化雪覆蓋指數(shù)NDSI的計算公式如下

式中，ρB2、ρB5分別代表Landsat TM 影像第2、第5波段的反射率。圖3第2行為4個冰川樣例圖中的歸一化雪蓋指數(shù)。

（2）全域閾值分割。采用初始閾值T0對NDSI圖像進行單閾值分割，獲取初步的冰雪單元［16］。已有研究表明，當 NDSI＞0.4時，能夠有效提取絕大部分的冰川和積雪［15］，這里也采用T0＝0.4進行閾值分割。由圖3可知，不同形態(tài)與狀態(tài)的冰川，采用單一閾值提取冰川均取得了一定的效果，明顯的冰川邊界均能精確地提取，部分區(qū)域的背景信息也有較大的NDSI值，被劃分為冰川單元。

（3）局域閾值分割：在全域閾值分割的每個冰雪單元，計算與其面積相等的緩沖區(qū)，作為局域分割單元（圖3第2行），并對局域單元的NDSI進行直方圖統(tǒng)計，可以看到冰雪與背景的雙峰分布圖，通過雙峰分布直方圖閾值選擇方法［12］確定新的分割閾值（圖3第3行）。圖3中，紅色邊界為全域分割冰雪邊界；綠色邊界為緩沖區(qū)邊界，NDSI直方圖中紅色豎線為局部分割閾值；黃色邊界為局部分割的冰雪邊界；圖3（d）中紫色邊界為水體邊界。4個冰川的最終的局部分割閾值分別為：0.376、0.324、0.553、0.475，這就使得不同類型的冰雪單元依據(jù)自身不同的特點實現(xiàn)精細的提取。

因此，針對不同形態(tài)和特征的冰川，采用不同的閾值方法提取可得到較好的效果。然而，圖3（a）、圖3（b）、圖3（d）的冰雪均有陰影的遮擋，且無法獲取陰影區(qū)域的冰雪邊界，部分地區(qū)的山體陰影也被視為冰川［17］

圖3 “全域－局部”冰雪提取以及其NDSI統(tǒng)計直方圖Fig.3 Snow／ice extraction based on“global－local”segmentation and their NDSI histogram

3.3 陰影區(qū)域冰雪邊界提取

為了識別陰影區(qū)域內(nèi)冰雪邊界，這里采用一定時期內(nèi)不同時相的遙感數(shù)據(jù)進行聯(lián)合分析。不同時相數(shù)據(jù)對同一地區(qū)進行觀測時具有不同的太陽角度，所形成的山體陰影也不同。通過對多期陰影信息的空間分析獲取陰影分布最小的冰雪范圍的方法如下：

（1）提取出多期影像受陰影影響的總和，確定多角度遙感影像總的陰影區(qū)域Smax

設(shè)有N個時相的遙感影像用于多角度剔除山體陰影，S（i）表示時相i影像的山體陰影區(qū)域，i∈［1，N］。

（2）為了使總陰影區(qū)域提取的冰川不受積雪干擾，確定總陰影區(qū)內(nèi)的冰川邊界。先求出時相i對應(yīng)的總的陰影區(qū)域的冰雪為GS（i）Gt（i）為3.2節(jié)方法提取的時相i冰雪。

然后，剔除總陰影區(qū)各個時相冰川受積雪的影響，得到總陰影區(qū)的冰川GS

（3）確定i時相剔除陰影后冰雪G（i）

G（i）表示時相i影像確定陰影區(qū)域冰川后的冰雪邊界，即G（i）剔除了山體陰影的影響。判斷多角度方法能完全剔除山體陰影的條件是若無法完全剔除陰影，標記未剔除的陰影。

圖4選取了4幅不同時相的Landsat數(shù)據(jù)進行陰影剔除試驗，成像日期分別為2011－02－17、2011－10－22、2011－09－29、2011－06－25，太陽高度角分別為31.20°、34.56°、42.32°、63.53°，太陽方位角分別為 150.42°、157.50°、152.28°、126.81°。圖4（a）、（b）、（c）、（d）的底圖分別為4個時相Landsat數(shù)據(jù)的波段7、4、2合成影像，圖中無填充和有填充的多邊形區(qū)域分別為利用其余時相依次剔除圖4（a）積雪和山體陰影后的結(jié)果。由圖4可見，試驗區(qū)在春、秋、冬季的影像受山體陰影面積較大，而夏季影像的山體陰影面積最小。陰影區(qū)的面積大小跟太陽高度角成反比。通過對4個陰影區(qū)的空間分析，可以獲取山體陰影的最小范圍（圖4（d））。雖然圖4（d）中仍有未剔除的4段陰影，但陰影對冰雪提取的影響降低到最低程度。在有更理想的夏季時相數(shù)據(jù)的情況下，可以完全消除陰影的影響。

3.4 多時相冰川邊界提取

上節(jié)利用多角度信息確定了山體陰影區(qū)域的冰川邊界，從而得到每個時相剔除山體陰影影響的冰雪信息。由于沒有氣象數(shù)據(jù)的輔助，無法確定單一時相遙感影像中哪些冰川邊界不受積雪影響。本文求出多時相冰雪最小邊界為理想的冰川邊界，在3.3節(jié)得到時相i剔除陰影影響后的冰雪為G（i），則多時相冰川邊界Gmulti的計算公式如下

G（i）表示時相i影像確定陰影區(qū)域冰川邊界后的冰雪邊界。圖4（a）、（b）、（c）、（d）中相近時間段內(nèi)各時相的冰雪分布不同，通過3.3節(jié)中的方法分別剔除四個時相的山體陰影，得到對應(yīng)時相完整的冰雪邊界G（1）、G（2）、G（3）、G（4），利用公式（5）求得冰川邊界Gmulti（見圖4（d）），最大限度地排除了積雪的影響。

圖4 山體陰影剔除與冰川提取示意圖Fig.4 The diagram to eliminate mountain shadow and to extract glacier

4 試驗與討論

為驗證本文提出方法的有效性，本文使用天山托木爾峰西部地區(qū)Landsat TM數(shù)據(jù)和SRTM高程數(shù)據(jù)進行多角度冰川提取試驗，試驗區(qū)的范圍為79.33°E至80.16°E，42.12°N 至42.55°N。遙感數(shù)據(jù)采用了2009年至2010年4個時相Landsat TM數(shù)據(jù)（見表1），軌道號為147／031。

利用第3節(jié)方法提取試驗區(qū)冰川的步驟為：①利用30m分辨率的DEM與表2中各時相Landsat TM數(shù)據(jù)觀測的太陽角度提取該區(qū)域的各個時相的山體陰影，分析得4個時相山體陰影的交集為空集；②通過NDSI“全域—局部”閾值迭代的方法提取單一時相對應(yīng)的冰雪邊界；③多時相的冰雪邊界與陰影邊界疊加運算，計算出每個時相山體陰影區(qū)域內(nèi)的最小的冰雪邊界，確定山體陰影區(qū)域的冰川邊界，得到多角度剔除山體陰影后的冰雪結(jié)果（圖5（a）、（b）、（c）、（d））；④用4個時相的剔除山體陰影后的冰雪邊界求最小值得到冰川邊界，得到的冰川邊界疊加在該區(qū)域2009年的RapidEye高分辨率（分辨率為5m）影像上（圖5（e））。圖5（a）、（b）、（c）、（d）的底圖分別為4個時相遙感數(shù)據(jù)的波段7、4、2合成影像，多邊形區(qū)域為對應(yīng)時相剔除陰影后的冰雪邊界。圖5（e）的底圖為2009年高分辨率影像，多邊形區(qū)域為利用圖5（a）、（b）、（c）、（d）的冰雪邊界提取的冰川邊界。

山體陰影區(qū)域的冰雪的比例表示對應(yīng)時相的冰雪受山體陰影影響的程度［18－19］，計算方法為某一時相山體陰影區(qū)域冰雪面積除以該時相的總冰雪面積。

時相n的冰雪面積Gn減去冰川面積Gmulti，得到該時相的積雪面積，再除以Gn，得到該時相的積雪比例。本文方法提取2009至2010年該區(qū)域的冰川面積為Gmulti＝568.7km2。由表2分析可得試驗區(qū)2009年9月14日山體陰影影響最小，2010年7月15日的積雪范圍最?。ㄒ姳?）。

表1 多角度冰川提取試驗Tab.1 Glacier mapping test with multiangle Landsat imagery

結(jié)合紋理、形態(tài)、鄰近度與DEM數(shù)據(jù)的高度和坡度信息，采用面向?qū)ο蠓椒梢钥焖偬崛「叻直媛视跋裰械谋ǎ?0］。但是對比高分辨率方法和本文方法（見圖6）。

圖5 多角度提取冰川邊界Fig.5 Delineation of glacier boundary using multi－angle glaciers

本文方法對于冰川提取顯示出很好的優(yōu)勢，具體分析如下：

（1）本文方法采用的Landsat數(shù)據(jù)時相較為齊全，而高分辨率影像的時相單一。尤其是在積雪影響較為嚴重的時期，單一時相的高分辨率影像提取的冰川可以通過形狀、光譜、地形等特征將冰川、積雪分離，雖然其描述冰川的輪廓極為精細，但是一些與冰川特征相似的積雪卻無法分離。而連續(xù)的多時相的Landsat影像，可以綜合不同時相的冰川邊界，最大限度地區(qū)分冰川和積雪的邊界。圖6（a）顯示了多角度提取的冰川邊界與受積雪影響較小的高分影像的冰川邊界線非常一致，而且在描述冰川內(nèi)部的裸地時，多角度方法要比面向?qū)ο蟮母叻直媛史椒蚀_。高分辨率影像在描述有很多細小分支的冰川時極為精細。但是，對于較長時間（30年或者更長時間）的冰川監(jiān)測而言，本文方法所用數(shù)據(jù)的時間的長度和密度都顯示出很大優(yōu)勢。

（2）Landsat影像與高分辨率影像的成像條件有一定的差異，尤其是在幾何形變較為嚴重高山地區(qū)表現(xiàn)尤為明顯，因此針對地勢起伏劇烈的高山地區(qū)冰川的提取，選擇近似正射投影（幾何形變較?。┑谋疚谋ㄌ崛》椒ㄝ^好（圖6（b））。

（3）圖6（c）通過面向?qū)ο蠓椒ㄌ崛「叻直媛视跋竦年幱皡^(qū)域的冰川，非冰川信息與冰川信息在光譜和地形特征等方面極為相似，不容易區(qū)分。圖6（d）顯示本文方法在相同山體陰影區(qū)域結(jié)合不同太陽角度下的冰川分布信息，能夠很好地提取出陰影區(qū)域的冰川邊界，并且能夠保證提取的陰影區(qū)域冰川邊界受積雪影響最小。

圖6 結(jié)果對比Fig.6 Comparison results

5 結(jié) 論

針對冰川提取容易受到積雪和山體陰影的影響，本文在“全域－局部”閾值分割提取冰雪的基礎(chǔ)上，用多角度多時相遙感數(shù)據(jù)對短時期內(nèi)的冰川進行提取，并對該方法提取冰川的精度進行了驗證。結(jié)論如下：

（1）本文用DEM和遙感成像的太陽角度信息建立山體陰影的物理模型。這種方法提取的山體陰影準確度高，與單純利用遙感影像的光譜信息提取陰影相比，減少了因“異物同譜”造成的錯分現(xiàn)象。

（2）冰川退縮是天山地區(qū)冰川變化的總體趨勢，冰川提取中把冰雪變化最劇烈的情況作為理想的冰川邊界［6］。從本文多角度山體陰影剔除方法可以看出，山體陰影的分布與太陽的高角度有關(guān)，在選擇多時相數(shù)據(jù)進行冰川提取時，可根據(jù)太陽高度信息來選擇合適的遙感數(shù)據(jù)進行冰川提取。

（3）多期遙感影像可以剔除積雪的影響。由結(jié)果分析可得，多角度遙感提取有效地融合了積雪范圍最小和山體陰影影響最小兩個最佳遙感觀測條件，提取的冰川最接近真實情況。

［1］ LI X，CHENG G D，JIN H J，et al.Cryospheric Change in China［J］.Global and Planetary Change，2008，62（3）：210－218.

［2］ PAN B T，ZHANG G L，WANG J，et al.Glacier Changes from 1966—2009in the Gongga Mountains，on the South－Eastern Margin of the Qinghai－Tibetan Plateau and Their Climatic Forcing ［J］. The Cryosphere，2012，6：1087－1101.

［3］ LI Zhen，SHI Jiancheng.Snow Mapping Algorithm Development and Validation Using Hyperspectral Data［J］.Acta Geodaetica et Cartographica Sinica，2001，30（1）：67－72.（李震，施建成.高光譜遙感積雪制圖算法及驗證［J］.測繪學報，2001，30（1）：67－72.）

［4］ ZHU Changming，LUO Jiancheng，SHEN Zhanfeng，et al.River Liner Water Adaptive Auto－extraction on Remote Sensing Image Aided by DEM［J］.Acta Geodaetica et Cartographica Sinica，2013，42（2）：277－283.（朱長明，駱劍承，沈占鋒，等.DEM輔助下的河道細小線性水體自適應(yīng)迭代提取［J］.測繪學報，2013，42（2）：277－283.）

［5］ HU Xiaodong，LUO Jiancheng，XIA Liegang，et al.Adaptive Water Body Information Extraction Using RS TUPU Computing Model［J］.Acta Geodaetica et Cartographica Sinica，2011，40（5）：544－550.（胡曉東，駱劍承，夏列鋼.圖譜迭代反饋的自適應(yīng)水體信息提取方法［J］.測繪學報，2011，40（5）：544－550.）

［6］ FREY H，PAUL F，STROZZI T.Compilation of a Glacier Inventory for the Western Himalayas from Satellite Data：Methods，Challenges，and Results［J］.Remote Sensing of Environment，2012，124（9）：832－843.

［7］ PAUL F，KAAB A，MAISCH M；et al.The New Remote Sensing Derived Swiss Glacier Inventory：I：Methods［J］.Annals of Glaciology，2002，34（1）：355－361.

［8］ BOLCH T.Climate Change and Glacier Retreat in Northern Tienshan（Kazakhstan／Kyrgyzstan）Using Remote Sensing Data［J］.Global and Planetary Change，2007，56（1）：1－12.

［9］ WANG P Y，Li Z Q，Li H L，et al.Comparison of Glaciological and Geodetic Mass Balance at Urumqi Glacier No.1，Tian shan，Central Asia［J］.Global and Planetary Change，2014，114：14－22.

［10］ KATZIL Y，DOYTSHER Y.A Logarithmic and Sub－pixel Approach to Shaded Relief Representation［J］.Computers＆ Geosciences，2003，29（9）：1137－1142.

［11］ RAUP B，RACOVITEANU A，KHALSA S J S，et al.The GLIMS Geospatial Glacier Database：A New Tool for Studying Glacier Change［J］. Global and Planetary Change，2007，56（1）：101－110.

［12］ LI J L，SHENG Y W.An Automated Scheme for Glacial Lake Dynamic Mapping with Landsat Imagery and Digital Elevation Models：A Case Study in the Himalayas［J］.International Journal of Remote Sensing，2012，33（16）：5194－5213.

［13］ XIAO J Y，JI N，LI X.A New Index for Steel Framed Roof Information Extraction Based on Remote Sensing［J］. Advanced Materials Research， 2013， 726：4682－4685.

［14］ YIN D，CAO X，CHEN X，et al.Comparison of Automatic Thresholding Methods for Snow－cover Mapping Using Landsat TM Imagery［J］.International Journal of Remote Sensing，2013，34（19）：6529－6538.

［15］ BISHOP M P，BONK R，KAMP J U，et al.Terrain Analysis and Data Modeling for Alpine Glacier Mapping［J］.Polar Geography，2001，25（3）：182－201.

［16］ YANG Shuwen，XUE Chongsheng，LIU Tao，et al.A Method of Small Water Information Automatic Extraction from TM Remote Sensing Images［J］.Acta Geodaetica et Cartographica Sinica，2010，39（6）：611－616.（楊樹文，薛重生，劉濤，等.一種利用TM影像自動提取細小水體的方法 ［J］.測繪學報，2010，39（6）：611－616.）

［17］ CHEN W B，F(xiàn)UKUI H，DOKO T，et al.Improvement of Glacial Lakes Detection under Shadow Environment U－sing ASTER Data in Himalayas，Nepal［J］.Chinese Geographical Science，2013，23（2）：217－226.

［18］ LIANG Dong，KONG Jie，HU Gensheng，et al.The Removal of Thick Cloud and Cloud Shadow of Remote Sensing Image Based on Support Vector Machine［J］.Acta Geodaetica et Cartographica Sinica，2012，41（2）：225－231，238.（梁棟，孔頡，胡根生，等.基于支持向量機的遙感影像厚云及云陰影去除［J］.測繪學報，2012，41（2）：225－231，238.）

［19］ CHENG Wei，WANG Liming，TIAN Qingjiu.A Method of Atmospheric Correction Based on Shadow－pixel for Optical Satellite Data［J］.Acta Geodaetica et Cartographica Sinica，2008，37（4）：469－475.（程偉，王黎明，田慶久.一種基于陰影像元的光學遙感大氣校正方法［J］.測繪學報，2008，37（4）：469－475.）

［20］ ZHONG Zhenwei，YE Qinghua.Integrated Method of Glacial Information Extraction around the Mount Naimonacnyi［J］.Journal of Glaciology and Geocryology，2009，13（4）：717－723.（仲振維，葉慶華.納木那尼峰地區(qū)冰川信息的 綜 合 提 取 方 法 ［J］.冰 川 凍 土，2009，13（4）：717－723.）