基于Landsat—8 OLI影像的南伊内里切克冰川运动速度提取与分析

田毅++李建学++何秦++杨森++李毅

摘 要：当前全球气候变化问题被日益关注，冰川作为气候变化的敏感指示器之一，对分析气候变化过程有重要科学意义。冰川表面运动是冰川的一个重要基本体征，它能够及时快速地反映冰川动态变化信息，为研究冰川对气候变化的响应提供了有效途径。本次研究以天山南伊内里切克冰川为研究对象，利用Landsat-8 OLI卫星影像，采用归一化互相关算法对2013-2016年的四期遥感影像进行分析处理，得到冰川表面的速度时空分布。通过分析得出，该冰川具有以下运动特征：冰川的轴部是冰川运动的主流线，在冰川轴线上，速度呈现先增大后减小的趋势，最大速度出现在平衡线附近；流速在轴部最大，在两侧呈现减小的趋势；2013-2016年间，该冰川主体流速约为34cm/d，呈现出较稳定的态势。

关键词：光学遥感；相关性分析；COSI-Corr；冰川流速

中图分类号：P343.6 文献标志码：A 文章编号：2095-2945（2017）20-0008-03

1 概述

全球气候变化是当今科学研究的热点问题。冰川作为冰冻圈系统的主要组成部分之一，也是天然的气候变化指示器，为了解水资源现状提供了研究途径[1]。迄今为止，前人的研究区域大多集中在南极和格陵兰两大冰盖等地区[2-4]，但只有少数研究是关于内陆山岳冰川的。谢自楚等在珠穆朗玛峰地区科学考察报告中指出，我国大多数的冰川数据资料都是于1960-1980年间获取的[5]；王欣等通过进一步研究，发现目前我国对于冰川的连续监测，仅局限于天山乌鲁木齐河源1号冰川等少数几条冰川[6]。而对于其他的冰川覆盖地区，则缺乏对它们运动的深入研究。

目前，冰川表面运动速度的提取主要有两种方式。一种是实地花杆数据测量[7]，这种方法的精度高，但是受时间和空间的局限性比较大；另一种方法则基于遥感影像进行提取[8]，相比传统野外实测方法，遥感获取冰川表面流速具有快速和成本低的优势，对于大面积的以及人类难以踏足的冰川地区，利用遥感影像的方式获取冰川表面速度的优势更加明显。

2 研究区介绍

托木尔峰地区南伊内里切克冰川（South Inilchek Glacier，中国冰川编目编号：5Y673K1；42°12.5′N，80°12.2′E）是我国典型的树枝状山谷冰川，冰川长约60.5千米，面积约567.2平方千米，是中低纬度区域长度超过50千米的八大冰川之一，图1是研究区域示意图。

3 研究方法

3.1 基本原理

对不同时期天山冰川的遥感影像进行相关性分析，进而提取同名像点的位移。经过图像预处理之后，可以去除太阳辐射、轨道偏差、地形影响等误差。再利用归一化互相关算法得到的同名像点的东西向和南北向位移，经过合成，就可以得到冰川表面的运动总位移，可以认为就是冰川运动的结果[9]，结合时间可以得到该阶段内冰川的日平均速度。

3.2 相关性分析及算法介绍

归一化互相关算法（NCC，Normalized Cross Correlation）[10]，也被称为泊松系数法或者归一化相关系数法。对于归一化互相关算法的表达式定义如下：

其中，f（x，y）搜索影像的搜索窗口，g（x，y）是原始影像的模板窗口，u，v是坐标的偏移量，u，v）是两个窗口各自的灰度平均值。如图2所示，m，n为模板窗口的长和宽，M，N是搜索窗口的长和宽，搜索窗口一般大于模板窗口。归一化处理能够增强该算法的鲁棒性，减少其他因素对相关系数带来的影响，同时降低相关系数范围对图像特征的依赖[11]。

在矩阵NCC（u，v）中，值为1的点表示搜索窗口中该点及其周围像素与模板窗口完全相同；值为-1的点表示搜索窗口中该点及周围像素与模板窗口完全相反。所以，通过选取值最大的点即可初步获得特征点的同名点，之后可以根据图像的坐标信息，计算出关键点的位移[12]。

4 数据选择与处理

与其他遥感影像数据相比，Landsat系列数据价格较低，获取途径多，分辨率较高时间序列较长的Landsat影像适合大区域长周期的冰川变化研究[13]。本次研究选用了4组Landsat-8 OLI影像，用来反演冰川的表面速度。

孔繁司[14]等在对四款常用的光学影像运动软件进行比较后发现，COSI-Corr軟件作为一款基于图像配准和相关系数计算的软件，该软件的配准精度可达到0.1个像元，因而获取到的冰川运动速度精度较高，而且花费的时间也较短。在计算最大相关系数时，采用频率域算法，参考窗口设为128，搜索窗口设为32。将位移数据分成了3个图层：东西向位移、南北向位移以及信噪比（Signal Noise Ratio，SNR）。为了提高研究结果的可信度，对窗口内的像素进行重采样，本研究在南伊内里切克冰川的冰舌区选取大量连续的点，选择SNR≥0.80的部分作为可信区域。

由于受到云雾和阴影的干扰，部分位移值不能真实反映表面运动，需要进行后处理。剔除云、阴影覆盖部分，即可得到冰川表面的位移量信息如图2所示，其中，图（a）（b）（c）分别对应2013-2014年总位移变化量、东西向位移变化量、南北向位移变化量；图（d）（e）（f）分别对应2014-2015年间总位移变化量、东西向位移变化量、南北向位移变化量；图（g）（h）（i）分别对应2015-2016年间总位移变化量、东西向位移变化量、南北向位移变化量。

由于获取的每组影像时间间隔不同，上述流程获取的只是对应时间间隔内的位移量，不方便对这3组冰川速度进行对比。因此，需要计算出年平均运动速度。首先计算夏季日平均运动速度，可由下式得出：

式中：D为两期遥感影像计算的位移数据；Days为两期影像时间段内夏季天数；Dayw为两期影像时间段内冬季天数。由此换算出年平均运动速度Vyr：

式中：A为冬季日平均运动速度与夏季日平均运动速度的比率，本研究采用邻区台兰冰川1978年测得的平均夏日改算系数0.699，由此换算A为0.482[15]。下表为三组影像间的夏季天数，冬季天数情况

通过公式对数据进行处理，得到图3，其中图（a）（b）（c）分别对应2013-2014年间、2014-2015年间、2015-2016年间冰川每天的位移量。

5 分析与讨论

5.1 冰川表面运动速度空间分布特征

冰川表面运动是底部变形和冰川底部滑动等共同作用的综合表现，受冰川厚度、地表坡度、物质平衡、冰温、冰内冰下水压等因素的影响，是一种十分复杂的机制。对比分析3期托木尔峰地区南伊内里切克冰川表面运动流速结果，我们可以发现该冰川具有如下运动特征：

（1）该冰川表面运动速度整体呈现平稳趋势。通过对冰川区域的6356个像素点进行统计分析，得到三年的日平均速度分别为34.92cm/d，33.45cm/d，34.92cm/d。

（2）对冰川在A线所示位置进行剖面分析，如图4。发现冰川主干区域在年间速度变化差异不大，但是2015-2016年间速度相对较高，可能与该地区的气候变化有关；冰川主干所在地区地形起伏变化较小，坡度在2°到17°之间，所以冰川的运动速度较平缓；随着地势降低，三个时间段的冰川速度出现分化趋势，在剖面线所在位置的中部分化最为明显，原因是在低海拔地区，冰温变化也会相对剧烈，导致冰川中部地区的年际变化差异。

（3）冰川运动的主流线在冰川的轴部，流速由轴部向两侧递减，由冰川源头向下至雪线处运动速度逐渐增加，然后再向冰川末端逐渐递减，这符合山地冰川运动的一般规律[16]。

（4）对比三期运动结果，由于受到上部冰川的重力和挤压力的作用，加之末端大量冰碛物的阻碍作用，导致冰川末端向北部移动，冰川速度和冰川的覆盖范围也呈现递减趋势。

5.2 精度评定

由于冰川区位于高寒高海拔的偏远山区，再加上恶劣自然条件的限制，难以长时间开展野外实地验证工作。为此，本文根据遥感冰川运动估算的基本原理，利用非冰川区残余位移对冰川表面运动监测结果开展了精度分析。理想状态下，非冰川地带是稳定、不存在位移的，因而我们可以选定非冰川地区的采样点作为精度评定的依据，将非冰川地带采样点的流速视作该方法的位移偏量，反映图像在匹配过程中产生的误差大小，从而用于评定结果精度。

采用圈定感兴趣区域的手段，选取了非冰川区采样点，以数理统计的方法，对其日运动速度进行了统计，发现非冰川区域点平均速度为1.37cm/d，约83%的取样点日平均速度在3cm/d以下。可见速度观测误差远比速度值小得多，说明了结果的可靠性。

6 结束语

本文以Landsat8-OLI卫星影像为数据源，借助归一化互相关算法，对天山托木尔峰地区南伊内里切克冰川表面运动速度进行了提取，监测结果表明：该地区冰川运动符合一般运动规律，2013-2016年间，该冰川主体流速约为34cm/d，呈现出较稳定的态势。另外，要想全面掌握该冰川的运动特征需要进一步扩展监测周期，获取更多的遥感影像，从而更好地研究该冰川在气候变暖条件下的物质平衡和运动时空演变特征，为深入了解该冰川动力学特征和预防冰川运动导致的地质灾害提供更为丰富的监测资料。

参考文献：

[1]K b A. Combination of SRTM3 and repeat ASTER data for deriving alpine glacier flow velocities in the Bhutan Himalaya[J]. Remote Sensing of Environment， 2005，94（4）：463-474.

[2]任贾文. 全球冰冻圈现状和未来变化的最新评估：IPCC WGI AR5 SPM发布[J].冰川冻土，2013，35（5）：1065-1067.

[3]沈永平，王国亚.IPCC第一工作组第五次评估报告对全球气候变化认知的最新科学要点[J].冰川冻土，2013，35（5）：1068-1076.

[4]Messerli A， Karlsson N B， Grinsted A. No slowing down of Jakobshavn Isbr in 2014.

[5]中国科学院西藏科学考察队.珠穆朗玛峰地区科学考察报告：1966-

1968现代冰川与地貌[M].北京：科学出版社，1975.

[6]王欣，刘琼欢，蒋亮虹，等.基于SAR影像的喜马拉雅山珠穆朗玛峰地区冰川运动速度特征及其影响因素分析[J].冰川冻土，2015，37（3）：570-579.

[7]井哲帆，叶柏生，焦克勤，等.天山奎屯河哈希勒根51号冰川表面运动特征分析[C].全国冰川冻土学大会暨冻土工程国际学术研讨会，20

02：563-566.

[8]蒋宗立，刘时银，韩海东，等.基于SAR数据的山地冰川表面运动速度分析[J].遥感技术与应用，2011，26（5）：640-646.

[9]黄磊，李震.光学遥感影像的山地冰川运动速度分析方法[J].冰川冻土，2009，31（5）：935-940.

[10]粱珊珊，闫世勇.基于合成孔径雷达影像估计慕士塔格峰地区冰川速度场[J].科学技术与工程，2014，14（23）：140-145.

[11]闫世勇.山地冰川表面运动雷达遥感监测方法研究[D].中国科学院大学，2013：15-18.

[12]牛牧野，周春霞，刘婷婷，等.基于改进NCC算法的东南极极记录冰川流速提取研究[J].极地研究，2016，28（2）：243-249.

[13]张国梁，王杰，潘保田，等.冰川变化遥感监测的研究进展[J].兰州大學学报（自然科学版），2010，46（6）：1-10.

[14]孔繁司，乔刚，王卫安.基于光学影像的冰流速测量软件比较与分析[J].中国科技论文在线精品论文，2016，9（12）：1240-1252.

[15]许君利，张世强，韩海东，等.天山托木尔峰科其喀尔巴西冰川表面运动速度特征分析[J].冰川冻土，2011，33（2）：268-275.

[16]曹泊，王杰，潘保田，等.祁连山东段宁缠河1号冰川和水管河4号冰川表面运动速度研究[J].冰川冻土，2013（06）：1428-1435.