青藏高原冰雪覆蓋變化監(jiān)測的衛(wèi)星立體影像方法

王聰華，劉 煒，杜鶴娟，林宗堅

(1. 西藏民族大學(xué)西藏光信息處理與可視化技術(shù)重點實驗室，陜西 咸陽 712082； 2. 中國測繪科學(xué)研究院，北京 100830)

青藏高原冰雪覆蓋變化監(jiān)測的衛(wèi)星立體影像方法

王聰華1，劉 煒1，杜鶴娟1，林宗堅2

(1. 西藏民族大學(xué)西藏光信息處理與可視化技術(shù)重點實驗室，陜西 咸陽 712082； 2. 中國測繪科學(xué)研究院，北京 100830)

針對青藏高原冰雪覆蓋變化監(jiān)測問題，以各拉丹冬冰川為例，提出了利用ALOS立體像對提取冰雪DEM，并以同源不同時相的冰雪DEM監(jiān)測冰雪覆蓋變化及體量變化的方法。依據(jù)DEM糾正裁切后的影像，采用ISODATA分類方法得到了冰雪覆蓋面積變化量，由不同時相網(wǎng)格DEM高程點數(shù)據(jù)高差與面積求得體積變化，得出了2009年與2010年冬季冰雪量大小。結(jié)果顯示2010年12月比2009年12月冰雪量減少了19.728 3 km3，雪覆蓋面積減少了349.691 km2。文中還列出了遇到的技術(shù)困難及有待進一步研究的問題。

冰雪覆蓋變化；立體像對；數(shù)字高程模型；冰雪立體信息；ALOS-PRISM

利用立體影像進行青藏高原冰雪覆蓋變化的監(jiān)測一直是攝影測量界十分關(guān)注、想做而有困難的工作。在過去的幾十年間利用遙感影像提取冰雪信息及進行冰雪監(jiān)測的研究很多。文獻[1]從積雪面積制圖、雪密度和雪深估算、雪水當量提取算法等方面分析了合成孔徑雷達(SAR)積雪監(jiān)測研究的國內(nèi)外進展情況。文獻[2]對冰雪監(jiān)測中常用的遙感衛(wèi)星(如Landsat、SPOT1—SPOT5、IKONOS、QuickBird、ALOS等)及其參數(shù)進行了分析和總結(jié)，并研究了冰川信息提取方法的發(fā)展。文獻[3—5]利用MODIS 影像數(shù)據(jù)，通過幾何校正、去云預(yù)處理，應(yīng)用歸一化差分積雪指數(shù)(NDSI)算法和綜合閾值可以獲取積雪覆蓋范圍和面積信息。文獻[6]研究了HJ-1B數(shù)據(jù)的雪蓋提取方法，利用NDVI 或TM 影像反演的林區(qū)輔助判識積雪，采用光譜與紋理信息結(jié)合的SVM 法提取雪蓋，但未涉及冰雪立體信息。目前獲取冰雪的立體信息主要通過地形圖、立體像對、雷達地形測繪及激光高度計，文獻[7]利用MSS/TM影像提取冰川面積變化, 研究利用地形圖DEM 和SRTM DEM 數(shù)據(jù)得到冰雪厚度信息。但目前未見有文獻表明利用青藏高原立體影像進行冰雪覆蓋的立體觀測問題。因此本文擬研究利用衛(wèi)星遙感立體影像觀測青藏高原冰雪覆蓋變化的技術(shù)方法。試驗采用各拉丹東雪山2009年冬季和2010年冬季兩個時相的4個ALOS-PRISM立體像對，研究ALOS立體像對提取冰雪DEM的方法，并以兩個不同時相的DEM研究冰雪覆蓋變化及立體信息變化問題。

1 試驗區(qū)域與數(shù)據(jù)準備

1.1 試驗區(qū)域介紹

試驗區(qū)域是各拉丹東雪山，長江的發(fā)源地，位于西藏北部與青海南部交界處的各拉丹東地區(qū)，其海拔高度為6621 m。區(qū)域地理坐標約為90°49′12″E—91°15′E，33°11′24″N—33°37′48″N。區(qū)域內(nèi)主要有兩座冰雪覆蓋的山峰，如圖1所示。

圖1 試驗區(qū)域各拉丹東

1.2 數(shù)據(jù)準備

從網(wǎng)站(https:∥auig2.jaxa.jp/ips/home?language= en_US#inline1)上搜索研究區(qū)域的ALOS-PRISM遙感影像存檔數(shù)據(jù)，根據(jù)存檔數(shù)據(jù)分析，選擇適合研究區(qū)域的ALOS PRISM遙感影像，全色2.5 m分辨率，帶RPC參數(shù)，云量小于2%。已經(jīng)購置的兩個時段的遙感數(shù)據(jù)PRISM 1B2立體像對見表1和表2。

表1 各拉丹東研究區(qū)域2009冬季影像數(shù)據(jù)

表2 各拉丹東研究區(qū)域2010冬季影像數(shù)據(jù)

2 基于ALOS立體像對的冰雪DEM提取

蔡慶空基于ALOS-PRISM立體影像對，采用ERDAS 9.2 LPS軟件提取了不同地形條件下的DEM[8]，并得到平坦裸露區(qū)DEM的高程中誤差為4.46 m，有植被覆蓋平原為4.56 m，城區(qū)為5.37 m，山區(qū)為9.2 m。研究表明，帶有RPC的ALOS立體像對提取DEM的精度是可靠的。但該文未涉及冰雪覆蓋區(qū)域DEM情況。

本文試驗借助ENVI遙感圖像處理軟件，利用ALOS立體像對實現(xiàn)冰雪覆蓋區(qū)域DEM的提取。

2.1 關(guān)于提取冰雪DEM方法問題

ENVI遙感圖像處理軟件具有立體像對DEM自動提取功能，對于帶RPC參數(shù)的立體像對可以自動完成DEM的提取，具體過程如圖2所示。對于質(zhì)量比較高的立體像對可以按照默認的設(shè)定，自動完成DEM提取或只要少量的人工調(diào)節(jié)。但對冰雪覆蓋面比較大的像對會造成諸多難點，采取了下列解決辦法：

圖2 冰雪DEM提取過程

(1) 同名點匹配困難，冰雪覆蓋的表面自動匹配連接點有很大困難，曾嘗試自動匹配連接點，結(jié)果能夠達到要求的連接點很少，因此需要在連接點上進行人工調(diào)整，使連接點在Y方向的視差盡量小，經(jīng)過不斷試驗，研究區(qū)域的立體像對在選擇連接點時，要求Y方向視差小于2個像素。

(2) 不能完全按默認參數(shù)設(shè)定生成DEM，試驗得出的經(jīng)驗是，移動窗口大小(moving window size)為13×13，地形細部(terrain detail)為Level 5,由于研究區(qū)域為山區(qū)，地形地貌(terrain relief)為High。

(3) 地面控制點問題，由于冰雪覆蓋區(qū)域總有裸露的大石頭，可以將裸露的石頭作為地面控制點。

圖3所示對應(yīng)于表1和表2立體像對提取的2009年和2010年冬季DEM數(shù)據(jù)。

圖3 兩個立體像對提取的冰雪DEM對照

2.2 誤差分析

由于試驗區(qū)域受冰雪覆蓋，獲得的DEM是含冰雪覆蓋的。在研究過程中，從地理空間數(shù)據(jù)云網(wǎng)站(http:∥www.gscloud.cn/)下載了研究區(qū)域的ASTER衛(wèi)星影像的30 m和90 m DEM，基準面為WGS-84。這些DEM數(shù)據(jù)是2009年冬季產(chǎn)品，理論上來說可以把ASTER衛(wèi)星影像的30 m DEM作為參照，用來與2009年ALOS立體像對提取的DEM進行比對。

利用ENVI軟件分別在ASTER影像30 m DEM和ALOS立體像對3 m DEM較平坦區(qū)域，均勻獲取了17個點的坐標值，經(jīng)統(tǒng)計，平面中誤差為2.25 m，高程中誤差為15.33 m。這樣的高程精度并不令人十分滿意，但不影響利用不同時相DEM解決冰雪量變化問題。筆者曾想利用PALSAR雷達影像提取DEM[9]，以此來檢查ALOS立體像對DEM的誤差問題，但經(jīng)查詢ALOS衛(wèi)星存檔數(shù)據(jù)，沒有相應(yīng)覆蓋試驗區(qū)域的PALSAR雷達影像數(shù)據(jù)。

3 冰雪覆蓋變化與立體信息提取

試驗影像分別是2009年12月和2010年12月的兩對立體像對，影像分辨率為2.5 m，生成DEM后分辨率為3 m，分辨率一致，WGS-84基準面和UTM投影，Zone 46 North，研究區(qū)域也是在無地質(zhì)災(zāi)害情況下的討論。

3.1 冰雪覆蓋面積信息提取

冰雪覆蓋及升降變化數(shù)據(jù)處理步驟如圖4所示。

由于研究區(qū)域兩對立體影像的覆蓋區(qū)域有所差別，所生成的DEM區(qū)域不完全相同，又由于DEM靠近邊緣的部分誤差較大，因此在研究冰雪覆蓋變化時需要對兩個時相DEM及影像進行裁剪，使得研究區(qū)域完全重疊。

3.1.1 ALOS影像裁剪與糾正

(1) 首先對2009年DEM裁剪，手動繪制感興趣多邊形區(qū)域，保存為ROI格式文件。利用感興趣多邊形區(qū)域?qū)?009年DEM進行裁剪，得到裁剪后的DEM文件。

圖4 冰雪覆蓋變化數(shù)據(jù)處理過程

(2) 利用前面繪制的感興趣多邊形區(qū)域?qū)?010年DEM進行裁剪，得到2009年和2010年數(shù)據(jù)區(qū)域一致的DEM，如圖5所示。

圖5 經(jīng)裁剪后的DEM對照

(3) 分別利用2009年和2010年的DEM對2009年和2010年影像進行校正。

(4) 將繪制的感興趣多邊形區(qū)域由ROI格式轉(zhuǎn)換成矢量格式EVF，再用矢量數(shù)據(jù)分別對2009年和2010年校正過的ALOS影像進行裁剪?？梢缘玫?009年和2010年數(shù)據(jù)區(qū)域一致影像數(shù)據(jù)，如圖6所示。

圖6 經(jīng)糾正裁剪后的影像對比

3.1.2 采用ISODATA分類方法提取積雪信息

ENVI軟件提供兩種非監(jiān)督分類方法：ISODATA和K-means方法。K-means使用聚類分析方法，隨機查找聚類族的聚類相似度，即中心位置，利用各聚類中對象的均值獲得一個中心對象來進行計算，然后迭代地重新配置，完成分類過程。

ISODATA算法是在K-means算法的基礎(chǔ)上，增加對聚類結(jié)果的合并和分裂兩個操作，并設(shè)定算法運行控制參數(shù)的一種聚類算法。

合并操作：當聚類結(jié)果某一類中樣本數(shù)太少，或兩個類間的距離太近時，進行合并。

分裂操作：當聚類結(jié)果某一類中樣本某個特征類內(nèi)方差太大，將該類進行分裂。

ISODATA算法將誤差平方和作為基本聚類準則，設(shè)定指標參數(shù)來決定是否進行合并或分裂，設(shè)定算法控制參數(shù)來決定算法的運算次數(shù)。

利用ISODATA分類方法提取冰雪覆蓋面積，其方法操作步驟如下：①執(zhí)行ISODATA非監(jiān)督分類；②類別定義與子類合并；③分類后處理與分類統(tǒng)計。

應(yīng)用Classification—Unsupervised—IsoData對試驗區(qū)域裁剪后的不同時相影像進行分類處理，再經(jīng)人工識別手工處理后的結(jié)果如圖7所示。分類統(tǒng)計后得到的冰雪覆蓋變化情況見表3。

圖7 經(jīng)ISODATA方法分類后的對比

統(tǒng)計項目非冰雪比例/(%)冰雪比例/(%)冰雪覆蓋面積/km22009年12月冰雪覆蓋11.4588.55821.3151332010年12月冰雪覆蓋49.1550.85471.624462同期冰雪覆蓋變化37.70-37.70-349.690671

從表3可以看出，2010年12月冰雪覆蓋面積比2009年12月減少了349.691 km2。本文在提取冰雪覆蓋面積時只是考慮了不同時相冰雪覆蓋面積的變化，并沒有著重考慮地形的問題。

3.2 不同時相冰雪變化量計算

基于2009年12月DEM(冰雪覆蓋變化前DEM)和2010年12月DEM(冰雪覆蓋變化后DEM)，利用ENVI的DEM提取功能，可以按像素點間隔規(guī)則劃分DEM，構(gòu)成規(guī)則網(wǎng)格DEM。規(guī)則網(wǎng)格DEM中的每個網(wǎng)格由4個頂點高程差組成近似小立方體，這樣可以遍歷冰雪覆蓋變化前DEM和冰雪覆蓋變化后DEM中所有網(wǎng)格，求對應(yīng)網(wǎng)格高程差值，由高程差值乘以每個格網(wǎng)對應(yīng)的面積值，再將所有網(wǎng)格上的體積變化累加，即可求得兩個時相冰雪變化量[10]。具體計算公式如下

(1)

式中，D是網(wǎng)格間隔， 為對應(yīng)i、j網(wǎng)格上的高程差；n、m對應(yīng)于DEM分割行數(shù)和列數(shù)；ΔV是計算得到的2009年12月與2010年12月間冰雪變化量。

可以利用IDL語言分別讀取2009年12月和2010年12月的DEM文件，以二維矩陣形式存儲，設(shè)定間隔為3個像素，根據(jù)3 m分辨率構(gòu)成9 m×9 m網(wǎng)格(D為9 m)，讀取DEM網(wǎng)格點，獲取網(wǎng)格4個頂點高程，根據(jù)式(1)可以統(tǒng)計出2009年12月—2010年12月冰雪變化量近似為19.728 3 km3。

4 結(jié) 語

本文從ALOS-PRISM立體像對提取冰雪DEM出發(fā)，討論了基于同源不同時相DEM求解冰雪變化量的方法，利用DEM的高程差值求得了2009年冬季與2010年冬季相比冰雪的升降變化量。針對ALOS影像，采用ISODATA分類方法得到了2009年冬季與2010年冬季相比冰雪覆蓋面積變化量。試驗表明，通過光學(xué)立體影像，同源不同時相的冰雪DEM求解冰雪變化量問題的關(guān)鍵在于立體像對獲得的DEM的精度，精度越高求解的冰雪變化量就越精確。盡管本文提取的冰雪DEM還存在一定誤差，但說明通過光學(xué)遙感影像用于立體觀測冰雪覆蓋面積變化量和冰雪升降變化量是可行的。試驗中遇到的若干技術(shù)難題雖然采取了應(yīng)對辦法，但是還有待進一步研究改善，也期望引發(fā)研討。

[1] 孫少波,車濤.基于合成孔徑雷達(SAR)的積雪監(jiān)測研究進展[J].冰川凍土,2013,35(3):636-647.

[2] 彥立利,王建.基于遙感的冰川信息提取方法研究進展[J].冰川凍土,2013,35(1):110-118.

[3] 侯慧姝,楊宏業(yè),王秀梅.基于MODIS影像的內(nèi)蒙古草原積雪監(jiān)測[J].測繪科學(xué),2010,35(4):117-119.

[4] 于秀娟,燕琴,蔡玉林.利用MODIS數(shù)據(jù)進行積雪檢測[J].測繪與空間地理信息,2010,33(4):38-41.

[5] 曾小箕，丁建麗,鄢雪英，等.基于MODIS數(shù)據(jù)的土庫曼斯坦山區(qū)積雪監(jiān)測[J].干旱區(qū)地理,2013,36(4): 717-719,721,722.

[6] 孫志群,劉志輝,邱冬梅.基于HJ-1B數(shù)據(jù)的雪蓋提取方法研究——以軍塘湖流域為例[J].干旱區(qū)地理, 2012,35(1):125-132.

[7] 王祎婷,陳秀萬,柏延臣,等.多源DEM和多時相遙感影像監(jiān)測冰川體積變化——以青藏高原那木納尼峰地區(qū)為例[J]. 冰川凍土,2010,32(1):126-132.

[8] 蔡慶空,蔣金豹,張玲,等. ALOS-PRISM立體像對提取DEM的應(yīng)用研究[J].測繪科學(xué),2014,39(1):70-73.

[9] 倪文儉,過志峰,孫國清.基于PALSAR數(shù)據(jù)的DEM提取方法研究[J].國土資源遙感,2009,81(3):19-23.

[10] 藍秋萍,費立凡,劉一寧.利用多時相DEM數(shù)據(jù)的冰雪覆蓋體量變化計算方法研究[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2010,35(10):1223-1225.

[11] 孫少波,車濤,王樹果,等.C波段SAR山區(qū)積雪面積提取研究[J].遙感技術(shù)與應(yīng)用,2013,28(3):444-452.

[12] HE Liye, LI Dongliang, CHEN Lian. Classification of Snow Cover Days in Western China and Comparison with Satellite Remote Sensing Data[J].Sciences in Cold and Arid Regions, 2012,4(3):0249-0258.

[13] MA Lijuan, QIN Dahe, BIAN Lingen, et al. Assessment of Snow Cover Vulnerability over the Qinghai-Tibetan Plateau[J].Advances In Climate Change Research,2011,2(2): 93-100.

[14] CHEN Jing, LI Rendong, YE Ming, et al. Monitoring Snow-cover Area Change in Antarctic Coastline Region Using MODIS Product Data[J].Chinese Journal of Polar Science, 2009,20(1):64-71.

[15] 張朝忙,劉慶生,劉高煥,等.我國東部地區(qū)ASTER GDEM高程精度評價[J].測繪科學(xué)技術(shù)學(xué)報,2014,31(1):67-72.

ChangeMonitoringofSnow-coverageBasedonALOS-PRISMStereoImages

WANG Conghua1，LIU Wei1，DU Hejuan1，LIN Zongjian2

(1. Xizang Key Laboratory of Optical Information Processing and Visualization Technology，Xizang Minzu University，Xianyang 712082，China; 2. Chinese Academy of Surveying and Mapping,Beijing 100830，China)

For the problem of snow cover land surface change monitoring, a new method was introduced taking Geladandong for instance. Using this method, the extraction of glacier DEM with ALOS-PRISM stereo-pair images was carried out, and snow-covered land surface change as well as volume change monitoring was applied with multi-temporal homology DEM images. Area variation of snow-covered land surface was obtained using ISODATA classing method based on rectified images by DEM images.Volume variation of glacier was obtained by area difference value and elevation difference value, which was derived from elevation-noted points in multi-temporal DEM images. Therefore, estimation of snow volume from 2009 to 2010 was implemented.Results showed that reduced snow volume was 19.728 3 km3and area of snow cover was 349.691 km2from December 2009 to December 2010.A method for extracting glacier DEM was put forward with stereo-pair images, and as well as for calculating volume variation of glacier with multi-temporal homology DEM images.

change of snow cover;stereo images;DEM;information of snow cover volume; ALOS-PRISM

王聰華，劉煒，杜鶴娟，等.青藏高原冰雪覆蓋變化監(jiān)測的衛(wèi)星立體影像方法[J].測繪通報，2017(11)：22-26.

10.13474/j.cnki.11-2246.2017.0341.

P237

A

0494-0911(2017)11-0022-05

2017-02-09;

2017-09-18

國家自然科學(xué)基金(61162025；41361044)

王聰華(1959—)，男，博士，教授，主要研究方向為攝影測量與遙感影像數(shù)據(jù)處理、民族文化數(shù)字化保護。E-mail:wangcshui@126.com