青藏高原冰雪覆盖变化监测的卫星立体影像方法

王聪华，刘 炜，杜鹤娟，林宗坚

(1. 西藏民族大学西藏光信息处理与可视化技术重点实验室，陕西 咸阳 712082； 2. 中国测绘科学研究院，北京 100830)

青藏高原冰雪覆盖变化监测的卫星立体影像方法

王聪华1，刘 炜1，杜鹤娟1，林宗坚2

(1. 西藏民族大学西藏光信息处理与可视化技术重点实验室，陕西 咸阳 712082； 2. 中国测绘科学研究院，北京 100830)

针对青藏高原冰雪覆盖变化监测问题，以各拉丹冬冰川为例，提出了利用ALOS立体像对提取冰雪DEM，并以同源不同时相的冰雪DEM监测冰雪覆盖变化及体量变化的方法。依据DEM纠正裁切后的影像，采用ISODATA分类方法得到了冰雪覆盖面积变化量，由不同时相网格DEM高程点数据高差与面积求得体积变化，得出了2009年与2010年冬季冰雪量大小。结果显示2010年12月比2009年12月冰雪量减少了19.728 3 km3，雪覆盖面积减少了349.691 km2。文中还列出了遇到的技术困难及有待进一步研究的问题。

冰雪覆盖变化；立体像对；数字高程模型；冰雪立体信息；ALOS-PRISM

利用立体影像进行青藏高原冰雪覆盖变化的监测一直是摄影测量界十分关注、想做而有困难的工作。在过去的几十年间利用遥感影像提取冰雪信息及进行冰雪监测的研究很多。文献[1]从积雪面积制图、雪密度和雪深估算、雪水当量提取算法等方面分析了合成孔径雷达(SAR)积雪监测研究的国内外进展情况。文献[2]对冰雪监测中常用的遥感卫星(如Landsat、SPOT1—SPOT5、IKONOS、QuickBird、ALOS等)及其参数进行了分析和总结，并研究了冰川信息提取方法的发展。文献[3—5]利用MODIS 影像数据，通过几何校正、去云预处理，应用归一化差分积雪指数(NDSI)算法和综合阈值可以获取积雪覆盖范围和面积信息。文献[6]研究了HJ-1B数据的雪盖提取方法，利用NDVI 或TM 影像反演的林区辅助判识积雪，采用光谱与纹理信息结合的SVM 法提取雪盖，但未涉及冰雪立体信息。目前获取冰雪的立体信息主要通过地形图、立体像对、雷达地形测绘及激光高度计，文献[7]利用MSS/TM影像提取冰川面积变化, 研究利用地形图DEM 和SRTM DEM 数据得到冰雪厚度信息。但目前未见有文献表明利用青藏高原立体影像进行冰雪覆盖的立体观测问题。因此本文拟研究利用卫星遥感立体影像观测青藏高原冰雪覆盖变化的技术方法。试验采用各拉丹东雪山2009年冬季和2010年冬季两个时相的4个ALOS-PRISM立体像对，研究ALOS立体像对提取冰雪DEM的方法，并以两个不同时相的DEM研究冰雪覆盖变化及立体信息变化问题。

1 试验区域与数据准备

1.1 试验区域介绍

试验区域是各拉丹东雪山，长江的发源地，位于西藏北部与青海南部交界处的各拉丹东地区，其海拔高度为6621 m。区域地理坐标约为90°49′12″E—91°15′E，33°11′24″N—33°37′48″N。区域内主要有两座冰雪覆盖的山峰，如图1所示。

图1 试验区域各拉丹东

1.2 数据准备

从网站(https:∥auig2.jaxa.jp/ips/home?language= en_US#inline1)上搜索研究区域的ALOS-PRISM遥感影像存档数据，根据存档数据分析，选择适合研究区域的ALOS PRISM遥感影像，全色2.5 m分辨率，带RPC参数，云量小于2%。已经购置的两个时段的遥感数据PRISM 1B2立体像对见表1和表2。

表1 各拉丹东研究区域2009冬季影像数据

表2 各拉丹东研究区域2010冬季影像数据

2 基于ALOS立体像对的冰雪DEM提取

蔡庆空基于ALOS-PRISM立体影像对，采用ERDAS 9.2 LPS软件提取了不同地形条件下的DEM[8]，并得到平坦裸露区DEM的高程中误差为4.46 m，有植被覆盖平原为4.56 m，城区为5.37 m，山区为9.2 m。研究表明，带有RPC的ALOS立体像对提取DEM的精度是可靠的。但该文未涉及冰雪覆盖区域DEM情况。

本文试验借助ENVI遥感图像处理软件，利用ALOS立体像对实现冰雪覆盖区域DEM的提取。

2.1 关于提取冰雪DEM方法问题

ENVI遥感图像处理软件具有立体像对DEM自动提取功能，对于带RPC参数的立体像对可以自动完成DEM的提取，具体过程如图2所示。对于质量比较高的立体像对可以按照默认的设定，自动完成DEM提取或只要少量的人工调节。但对冰雪覆盖面比较大的像对会造成诸多难点，采取了下列解决办法：

图2 冰雪DEM提取过程

(1) 同名点匹配困难，冰雪覆盖的表面自动匹配连接点有很大困难，曾尝试自动匹配连接点，结果能够达到要求的连接点很少，因此需要在连接点上进行人工调整，使连接点在Y方向的视差尽量小，经过不断试验，研究区域的立体像对在选择连接点时，要求Y方向视差小于2个像素。

(2) 不能完全按默认参数设定生成DEM，试验得出的经验是，移动窗口大小(moving window size)为13×13，地形细部(terrain detail)为Level 5,由于研究区域为山区，地形地貌(terrain relief)为High。

(3) 地面控制点问题，由于冰雪覆盖区域总有裸露的大石头，可以将裸露的石头作为地面控制点。

图3所示对应于表1和表2立体像对提取的2009年和2010年冬季DEM数据。

图3 两个立体像对提取的冰雪DEM对照

2.2 误差分析

由于试验区域受冰雪覆盖，获得的DEM是含冰雪覆盖的。在研究过程中，从地理空间数据云网站(http:∥www.gscloud.cn/)下载了研究区域的ASTER卫星影像的30 m和90 m DEM，基准面为WGS-84。这些DEM数据是2009年冬季产品，理论上来说可以把ASTER卫星影像的30 m DEM作为参照，用来与2009年ALOS立体像对提取的DEM进行比对。

利用ENVI软件分别在ASTER影像30 m DEM和ALOS立体像对3 m DEM较平坦区域，均匀获取了17个点的坐标值，经统计，平面中误差为2.25 m，高程中误差为15.33 m。这样的高程精度并不令人十分满意，但不影响利用不同时相DEM解决冰雪量变化问题。笔者曾想利用PALSAR雷达影像提取DEM[9]，以此来检查ALOS立体像对DEM的误差问题，但经查询ALOS卫星存档数据，没有相应覆盖试验区域的PALSAR雷达影像数据。

3 冰雪覆盖变化与立体信息提取

试验影像分别是2009年12月和2010年12月的两对立体像对，影像分辨率为2.5 m，生成DEM后分辨率为3 m，分辨率一致，WGS-84基准面和UTM投影，Zone 46 North，研究区域也是在无地质灾害情况下的讨论。

3.1 冰雪覆盖面积信息提取

冰雪覆盖及升降变化数据处理步骤如图4所示。

由于研究区域两对立体影像的覆盖区域有所差别，所生成的DEM区域不完全相同，又由于DEM靠近边缘的部分误差较大，因此在研究冰雪覆盖变化时需要对两个时相DEM及影像进行裁剪，使得研究区域完全重叠。

3.1.1 ALOS影像裁剪与纠正

(1) 首先对2009年DEM裁剪，手动绘制感兴趣多边形区域，保存为ROI格式文件。利用感兴趣多边形区域对2009年DEM进行裁剪，得到裁剪后的DEM文件。

图4 冰雪覆盖变化数据处理过程

(2) 利用前面绘制的感兴趣多边形区域对2010年DEM进行裁剪，得到2009年和2010年数据区域一致的DEM，如图5所示。

图5 经裁剪后的DEM对照

(3) 分别利用2009年和2010年的DEM对2009年和2010年影像进行校正。

(4) 将绘制的感兴趣多边形区域由ROI格式转换成矢量格式EVF，再用矢量数据分别对2009年和2010年校正过的ALOS影像进行裁剪。可以得到2009年和2010年数据区域一致影像数据，如图6所示。

图6 经纠正裁剪后的影像对比

3.1.2 采用ISODATA分类方法提取积雪信息

ENVI软件提供两种非监督分类方法：ISODATA和K-means方法。K-means使用聚类分析方法，随机查找聚类族的聚类相似度，即中心位置，利用各聚类中对象的均值获得一个中心对象来进行计算，然后迭代地重新配置，完成分类过程。

ISODATA算法是在K-means算法的基础上，增加对聚类结果的合并和分裂两个操作，并设定算法运行控制参数的一种聚类算法。

合并操作：当聚类结果某一类中样本数太少，或两个类间的距离太近时，进行合并。

分裂操作：当聚类结果某一类中样本某个特征类内方差太大，将该类进行分裂。

ISODATA算法将误差平方和作为基本聚类准则，设定指标参数来决定是否进行合并或分裂，设定算法控制参数来决定算法的运算次数。

利用ISODATA分类方法提取冰雪覆盖面积，其方法操作步骤如下：①执行ISODATA非监督分类；②类别定义与子类合并；③分类后处理与分类统计。

应用Classification—Unsupervised—IsoData对试验区域裁剪后的不同时相影像进行分类处理，再经人工识别手工处理后的结果如图7所示。分类统计后得到的冰雪覆盖变化情况见表3。

图7 经ISODATA方法分类后的对比

统计项目非冰雪比例/(%)冰雪比例/(%)冰雪覆盖面积/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提取功能，可以按像素点间隔规则划分DEM，构成规则网格DEM。规则网格DEM中的每个网格由4个顶点高程差组成近似小立方体，这样可以遍历冰雪覆盖变化前DEM和冰雪覆盖变化后DEM中所有网格，求对应网格高程差值，由高程差值乘以每个格网对应的面积值，再将所有网格上的体积变化累加，即可求得两个时相冰雪变化量[10]。具体计算公式如下

(1)

式中，D是网格间隔， 为对应i、j网格上的高程差；n、m对应于DEM分割行数和列数；ΔV是计算得到的2009年12月与2010年12月间冰雪变化量。

可以利用IDL语言分别读取2009年12月和2010年12月的DEM文件，以二维矩阵形式存储，设定间隔为3个像素，根据3 m分辨率构成9 m×9 m网格(D为9 m)，读取DEM网格点，获取网格4个顶点高程，根据式(1)可以统计出2009年12月—2010年12月冰雪变化量近似为19.728 3 km3。

4 结 语

本文从ALOS-PRISM立体像对提取冰雪DEM出发，讨论了基于同源不同时相DEM求解冰雪变化量的方法，利用DEM的高程差值求得了2009年冬季与2010年冬季相比冰雪的升降变化量。针对ALOS影像，采用ISODATA分类方法得到了2009年冬季与2010年冬季相比冰雪覆盖面积变化量。试验表明，通过光学立体影像，同源不同时相的冰雪DEM求解冰雪变化量问题的关键在于立体像对获得的DEM的精度，精度越高求解的冰雪变化量就越精确。尽管本文提取的冰雪DEM还存在一定误差，但说明通过光学遥感影像用于立体观测冰雪覆盖面积变化量和冰雪升降变化量是可行的。试验中遇到的若干技术难题虽然采取了应对办法，但是还有待进一步研究改善，也期望引发研讨。

[1] 孙少波,车涛.基于合成孔径雷达(SAR)的积雪监测研究进展[J].冰川冻土,2013,35(3):636-647.

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

[3] 侯慧姝,杨宏业,王秀梅.基于MODIS影像的内蒙古草原积雪监测[J].测绘科学,2010,35(4):117-119.

[4] 于秀娟,燕琴,蔡玉林.利用MODIS数据进行积雪检测[J].测绘与空间地理信息,2010,33(4):38-41.

[5] 曾小箕，丁建丽,鄢雪英，等.基于MODIS数据的土库曼斯坦山区积雪监测[J].干旱区地理,2013,36(4): 717-719,721,722.

[6] 孙志群,刘志辉,邱冬梅.基于HJ-1B数据的雪盖提取方法研究——以军塘湖流域为例[J].干旱区地理, 2012,35(1):125-132.

[7] 王祎婷,陈秀万,柏延臣,等.多源DEM和多时相遥感影像监测冰川体积变化——以青藏高原那木纳尼峰地区为例[J]. 冰川冻土,2010,32(1):126-132.

[8] 蔡庆空,蒋金豹,张玲,等. ALOS-PRISM立体像对提取DEM的应用研究[J].测绘科学,2014,39(1):70-73.

[9] 倪文俭,过志峰,孙国清.基于PALSAR数据的DEM提取方法研究[J].国土资源遥感,2009,81(3):19-23.

[10] 蓝秋萍,费立凡,刘一宁.利用多时相DEM数据的冰雪覆盖体量变化计算方法研究[J].武汉大学学报(信息科学版),2010,35(10):1223-1225.

[11] 孙少波,车涛,王树果,等.C波段SAR山区积雪面积提取研究[J].遥感技术与应用,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] 张朝忙,刘庆生,刘高焕,等.我国东部地区ASTER GDEM高程精度评价[J].测绘科学技术学报,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

王聪华，刘炜，杜鹤娟，等.青藏高原冰雪覆盖变化监测的卫星立体影像方法[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

国家自然科学基金(61162025；41361044)

王聪华(1959—)，男，博士，教授，主要研究方向为摄影测量与遥感影像数据处理、民族文化数字化保护。E-mail:wangcshui@126.com