0 引言

数字高程模型(digital elevation model，DEM)，是用一组有序数值阵列形式表示地面高程的一种实体地面模型[1]。由于DEM数据变化检测分析可以为地表监测、灾情分析、气候分析等方面提供数据支持，国内外众多专家学者都从不同方面对这一问题进行了研究探讨。李德仁等在2006年提出利用新时期的航空影像和旧时期的DEM数据，基于正射影像匹配进行地形变化检测和数据更新的自动化算法[2]。夏松等在2007年提出利用多源空间数据进行地形变化检测的方法，利用空间数据库中已有的DEM、DOM、DLG和新影像进行三维变化检测[3]。冯甜甜等在2008年提出基于物方影像匹配的DEM聚类变化检测方法，采用了基于密度的聚类分析方法对匹配所得的初始变化区域进行分析，实现了地表高程变化区域的检测[4]。曹建君等在2011年将面向对象的变化检测方法运用于舟曲泥石流受灾区域DEM的变化检测之中，探讨改进了利用多源空间数据进行DEM变化检测的方法[5]。王祎婷等在2012年提出利用多源DEM和多时相遥感影像数据，监测青藏高原那木纳尼峰地区冰川体积变化方法[6]。Di Wang与Andreas Kääb 于2015年通过建立DEM数据的时间序列，计算出福克斯冰川和弗兰兹-约瑟夫冰川在2000年至2014年间的海拔和体积变化情况[7]。

为保证DEM变化检测分析的精度，必须在此之前对数据进行配准处理，将来自不同时间、不同传感器成像条件的数据统一到相同的地理坐标系中。针对多期DEM数据配准对匹配算法精度与效率要求较高的特点，本文采用了尺度不变特征变换(scale-invariant feature transform，SIFT)匹配算法。SIFT是David Lowe于1999年提出的局部特征描述子[8]，并于2004年进行了更深入的发展和完善[9]。SIFT算法具有尺度不变性，可在图像中检测出大量关键特征点，处理2幅图像之间发生平移、旋转、仿射变换情况下的匹配问题，匹配能力较强，具有稳定性、独特性、多量性、高速性等优势。在Mikolajczyk于2005年对包括SIFT算子在内的多种局部描述子所做出的不变性对比实验中，证实SIFT及其扩展算法在同类描述子中具有最强的健壮性[10]。

本文采用SIFT算法在多期DEM数据中提取多组坐标匹配点并进行配准，将多期DEM数据统一到一致的地理参考系下。在成像时间最早的一期数据中提取出冰川边界，作为变化检测的基准。在多时相数据间进行数学运算构造差异，生成差分图像，根据统计学原理减少变化图像中的系统误差，设定合理阈值后进行变化检测，生成专题图、统计表格等变化检测分析成果，并对变化检测结果的精度进行评定。

1 研究区域与数据

1.1 研究区域

研究区域位置地处我国西藏自治区山南地区浪卡子县附近，影像数据覆盖面积约为3 600 km2，区域内主要山峰为宁金岗桑峰与卡鲁雄峰，冰川数目多，分布广。由于研究区域属于人口聚集区，区域内冰川具有典型中低纬度山地冰川的特征，对其进行变化检测对当地气候条件变化、当地居民的生产生活、下游河流及湖泊水量平衡等研究方面都具有重要意义。研究区域位置如图1所示。

1.2 研究数据

本研究所使用的数据集为高级星载热发射和反射辐射仪(advanced spaceborne thermal emission and reflection radiometer，ASTER)数据。由于ASTER数据传感器获取立体像对的方向为沿轨道方向，能够增强光照条件的一致性，使得地形数据质量更好[11]。

选择了同一区域内的4组多时相ASTER数据，其成像时间分别为2004年、2007年、2011年与2014年，数据集时间分布均匀，冰川范围变化较为直观，可以满足冰川变化检测的需求。

2 研究方法

本文提出的冰川变化检测研究方法通过采用SIFT特征提取与特征匹配算法，从多期DEM数据中提取同名点坐标并配准，在此基础上对具有相同地理坐标参考的多期DEM数据进行变化检测分析，其总体流程如图2所示。

2.1 数据预处理

为得到可用于SIFT匹配算法的数据，需要通过线性内插方法建立高程数据与图像灰度域(0～255)间的映射关系，将包含高程信息的DEM数据渲染为灰度影像。

2.2 SIFT算法配准多期地形数据

SIFT算法首先在尺度空间进行特征检测，寻找极值点，并确定极值点的位置及其所处的尺度，过滤出稳定的特征点，生成形成与尺度和方向的无关性局部描述算子[12]，并将其应用于匹配计算中。

本文使用opencv计算机视觉类库中的自带SIFT特征点检测与特征点匹配函数，提取出DEM影像中的特征点，经过筛选计算后生成匹配结果。以2004年与2007年DEM影像配准为例，实验中生成18组匹配点，配准的总体均方根误差RMS为0.458个像元，符合变化检测要求。

通过在4组地形数据间进行3次DEM匹配运算，生成3组同名匹配点坐标值。提取同名匹配点后，对4组DEM地形数据分别进行3次地理配准操作，将地形数据统一到同一个地理坐标参考系之下。

2.3 提取冰川范围

以成像时间最早(2004年)的地形数据作为基准提取冰川边界范围，在此基础上研究冰川变化的情况。为提高提取算法的自动化程度，可以通过计算ASTER影像的归一化雪被指数(normalized difference snow index，NDSI)并对生成的结果设定合理阈值，进行二值化并处理后得到冰川边界数据。

NDSI计算公式如(1)式所示：

(1)

式中：DNASTER1、DNASTER4分别为ASTER1波段与ASTER 4波段图像的灰度值。通过对NDSI计算结果进行密度分割，将其阈值选择为0.4，即NDSI值大于0.4的地物可认为是冰雪。

由于使用NDSI法提取出的冰川范围时会产生较多细小图斑，还需要结合目视解译方法进行修改，生成符合要求的冰川边界多边形矢量数据。

获得冰川边界多边形矢量数据后，通过对比多组DEM数据冰川分布规律，设置缓冲距离为1 500 m建立缓冲区。缓冲区内的地表高程变化归类为冰川变化，缓冲区外的则认为是其他地物种类发生的变化。下文进行的变化检测即在缓冲区范围内进行。

2.4 冰川变化检测

按照成像时间顺序，在4组DEM数据间进行差分运算，获得3张差值图像序列，如图3所示。

图3(a)为2004年与2007年数据生成的差值图像，图3(b)为2007年与2011年数据生成的差值图像，图3(c)为2011年与2014年数据生成的差值图像。

为由于消除光照强度、传感器性能、气候条件等因素差异产生的系统误差，本文研究中利用统计学原理，通过建立差值图像灰度值均值和标准差的数学关系来确定变化阈值[13]，具体方法如式(2)：

(2)

式中：μ为差值图像灰度值均值；σ为差值图像灰度值标准差；N为门限值，默认值可取为1，通过多次比对并参考图像灰度直方图分布情况调整其取值；T1、T2为计算得到的阈值。

变化阈值确定后，即可对差值图像进行分类，分类具体思路如式(3)所示：

(3)

以第一组差值图像数据为例进行分析，首先生成该差值图像灰度分布直方图，如图4所示。

由统计信息可以得到该差值图像灰度值均值为34，标准差为86。从灰度直方图分布情况发现该图像灰度分布集中程度较高，需要通过选取一个略小的门限值，增强变化的检测灵敏度。通过对比，将门限值N值确定为0.9，求得变化阈值T1=121，T2=-33。根据变化阈值对差值图像灰度值进行划分区间，获得如图5所示的冰川变化检测结果图，其中冰川地区黑色部分为冰川消减区域，灰色部分为冰川增长区域。

由变化检测结果统计得到2004至2007年间，研究区域发生减退的冰川面积为215.83 km2，增长的冰川面积为28.04 km2，冰川总体减退面积为187.79 km2，年均减退速率为62.60(km2·a-1)，冰川退化情况较为严重。

对另外2组差值图像进行同样的变化检测流程，得到以下2组冰川变化检测结果，如图6所示。图6(a)为2007年至2011年间冰川变化结果，图6(b)为2011年至2014年冰川变化结果，其中冰川地区黑色部分为冰川消减区域，灰色部分为冰川增长区域。

3 研究结果与分析

3.1 变化检测结果专题图

在变化检测结果影像的基础上，添加图名、指北针、图例、比例尺等制图要素，生成多张变化检测结果专题图，如图7所示。其中图7(a)为2004—2007年冰川变化检测结果，图7(b)为2007—2011年冰川变化检测结果，图7(c)为2011—2014年冰川变化检测结果。

3.2 冰川变化面积与类型统计

2004年至2014年间冰川变化面积与类型统计结果如表1所示。

从表1可以得出，研究区域内的冰川在2004年至2014年发生了较为显著的面积变化，变化幅度上下波动较大。其中2004年至2007年间冰川面积减少52.23%，年均变化率为17.41%；2007年至2011年间冰川面积增加148.49%，年均变化率为37.12%；2011年至2014年间冰川面积减少32.42%，年均变化率为10.81%。

冰川的变化情况可分为增长与消减2种类型，通过统计2004年至2014年间冰川变化不同类型的面积所占的百分比，发现2004年至2007年与2011年至2014年2个时间段内，冰川变化类型以消减为主；2007年至2011年时间段内，冰川范围的增长在变化类型中占据优势。冰川总体的变化速率在2004年至2007年间为81.29(km2·a-1)，2007年至2011年间为87.07(km2·a-1)，2007年至2011年间为118.51(km2·a-1)，随着时间发展，冰川面积变化的速率不断加快，这也体现了近年来极端气候频发等因素对地形地貌产生了重大影响[14]。

结合研究区域2004年至2014年间的年均气温数据与冰川面积，生成统计图，如图8所示。

由统计图10可以发现冰川面积的变化与气温变化因素有着密切联系[15]。由年均气温变化折线可以看出，2004年至2007年间，年均气温呈上升趋势，冰川发生了程度较大的消减。2007年至2011年平均气温下降程度较大，2011年的平均气温在10年间处于最低水平，故该年的冰川面积与2007年相比大幅度增长。2011年至2014年间年均气温重新升高，冰川范围再次进入消减状态。

3.3 变化检测精度统计

变化检测精度统计采用抽样调查方法原理，每次变化检测后在冰川缓冲区范围内生成50个随机点，结合目视解译方法分析各次变化检测结果的正确性，最后统计变化检测精度，统计结果如表2所示。

统计出变化检测总体精度为85.33%。

4 结束语

DEM数据具有能够精确反映地表高程的特性，通过对其变化状态进行检测和分析，可以得到许多具有重要应用价值的信息。本论文完成的主要研究工作包括以下几个方面：

1)使用SIFT匹配算法，提取多组地形数据中的大量特征点，生成多组匹配，并按照距离比例准则去除错误的匹配。根据匹配点坐标，对多组地形数据进行配准操作，进行几何变换与灰度重采样，统一到一致的地理坐标参考系之下。

2)计算基准地形数据的NDSI值，提取冰川边界多边形数据，建立缓冲区确定变化检测范围。在多期地形数据间进行差分运算，生成多组差值图像，通过建立图像灰度值与标准差间数学关系确定变化阈值，生成冰川变化检测结果。

3)根据多期冰川变化检测结果生成不同类型的专题图与专题表，结合当地气候条件分析冰川变化的原因。通过变化检测精度评价，证明本文所采用的研究方法将变化检测误差控制在可接受范围内。

本论文在变化检测方法与结果分析等方面还存在着一些问题需要改进。基于图像像素的变化检测方法由于受成像条件、配准精度等因素影响，难以完全消除系统误差的干扰。若在变化检测前预先对地形数据进行聚类计算等处理，能够在一定程度上提高变化检测的精度。变化结果分析中若能采用时间间隔更短的地形数据和精度更高的气象数据，则可以建立气象数据与冰川变化数据的数学函数，进行更为精确的定量分析。

参考文献

[1] 刘锦中，马辉.基于等高线的DEM生成算法研究和实现[J].现代测绘，2004，27(3)：34-35.

[2] 李德仁，夏松，江万寿.基于正射影像匹配的地形变化检测与更新算法[J].地理与地理信息科学，2006，22(6)：9-11.

[3] 夏松，李德仁，巫兆聪.利用多源空间数据进行地形的三维变化检测[J].测绘科学，2007，32(1)：49-50.

[4] 冯甜甜，王密.基于物方影像匹配的DEM聚类变化检测[J].测绘信息与工程，2008，33(1)：35-36.

[5] 曹建君，骆莹，李德仁，等.利用多源空间数据进行的DEM变化检测方法研究[J].矿山测量，2011(5)：17-20.

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

[7] WANG D，KAAB A.Modeling glacier elevation change from DEM time series[J].Remote Sens,2015，7，10117-10142.

[8] LOWE D G.Object recognition from local scale-invariant features[C]∥International Conference on Computer Vision，Corfu，Greece，1999(9)：1150-1157.

[9] LOWE D G.Distinctive image features from scale-invariant keypoints[J].International Journal of Computer Vision，2004，60(2)：91-110.

[10] MIKOLAJCZYK K，TUYTELAARS T，SCHMID C，et al.A comparison of affine region detectors[J].International Journal of Computer Vision，2005，65(1/2)：43-72.

[11] 李海涛，田庆久.ASTER数据产品的特性及其计划介绍[J].遥感信息，2004,19(3)：53-55.

[12] 李铁军，陈哲，王任享.尺度不变特征变换算法在图像配准中的应用[J].弹箭与制导学报，2008，28(2)：183-185.

[13] 杨希，刘国祥，秦军，等.基于多时相遥感图像灰度差值法的地表变化检测[J].四川测绘，2008，31(3)：99-103.

[14] KHALSA S，DYURGEROV M，KHROMOVA T，et al.Space-based mapping of glacier changes using ASTER and GIS tools[J].IEEE Transactions on Geoscience & Remote Sensing，2004，42(10)：2177-2183.

[15] 拉巴卓玛，邱玉宝，除多，等.1972～2010年西藏卡若拉冰川面积变化遥感研究[J].遥感技术与应用，2015，30(4)：784-792.