利用立体影像进行青藏高原冰雪覆盖变化的监测一直是摄影测量界十分关注、想做而有困难的工作。在过去的几十年间利用遥感影像提取冰雪信息及进行冰雪监测的研究很多。文献[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.2 数据准备

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

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提取或只要少量的人工调节。但对冰雪覆盖面比较大的像对会造成诸多难点，采取了下列解决办法：

(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数据。

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文件。

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

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

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

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

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

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

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

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

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

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

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

从表 3可以看出，2010年12月冰雪覆盖面积比2009年12月减少了349.691 km²。本文在提取冰雪覆盖面积时只是考虑了不同时相冰雪覆盖面积的变化，并没有着重考虑地形的问题。

3.2 不同时相冰雪变化量计算

基于2009年12月DEM(冰雪覆盖变化前DEM)和2010年12月DEM(冰雪覆盖变化后DEM)，利用ENVI的DEM提取功能，可以按像素点间隔规则划分DEM，构成规则网格DEM。规则网格DEM中的每个网格由4个顶点高程差组成近似小立方体，这样可以遍历冰雪覆盖变化前DEM和冰雪覆盖变化后DEM中所有网格，求对应网格高程差值，由高程差值乘以每个格网对应的面积值，再将所有网格上的体积变化累加，即可求得两个时相冰雪变化量^[10]。具体计算公式如下

式中，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 km³。

4 结语

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