基于Landsat-8时间序列影像分析西昆仑山地区冰川滑移特征

郑　茜　孙建宝　张　永

1　中国地震局地质研究所地震动力学国家重点实验室,北京市华严里甲1号,100029

基于Landsat-8时间序列影像分析西昆仑山地区冰川滑移特征

郑茜1孙建宝1张永1

1中国地震局地质研究所地震动力学国家重点实验室,北京市华严里甲1号,100029

摘要：利用图像亚像元互相关分析方法处理了Landsat-8卫星获取的时间序列影像数据，得到中国青藏高原西北部地区西昆仑峰区冰川匀速滑移的时空演化过程。利用亚像元影像互相关技术对Landsat-8光学影像精确配准，配准精度达到0.01像元，即该光学影像的水平形变监测精度达到0.15 m。通过对2013-07～2014-08的15景Landsat-8影像进行互相关和形变时间序列反演分析，获得了西昆仑峰区两条冰川的滑动位移场和速度场。研究表明，该区域的冰川基本处于匀速滑移状态(无明显加速和减速现象)；同时也验证了Landsat-8光学影像在监测较大地表位移和地壳形变事件(如沙丘移动、地震、滑坡、火山等)上的应用潜力。

关键词：西昆仑峰区；Landsat-8；配准；冰川滑移；亚像元图像互相关

依据中国冰川编目，西昆仑山共有冰川5 485条，冰川面积8 817.78 km2，约占整个昆仑山冰川总面积的3/4以上[1]。来自昆仑山冰川的冰川水是塔里木盆地的主要水源，为生活在这一流域的近630万人口提供饮用水[2-3]。分析研究这些高山冰川的位移场和它们的动态时空演化过程，对了解青藏高原地区的气候变化、冰川活动对人类活动的影响等有重要意义。

目前，常用的形变监测技术包括差分全球定位系统技术(DGPS)、干涉合成孔径雷达技术 (InSAR)和亚像元光学影像互相关技术等。前两种方法虽然可以获得较高的形变观测精度，但在监测高山冰川方面具有较大的局限性[4]。在过去的几十年里，美国的Landsat-1～7系列卫星所获得的数据广泛用于地球科学研究和人类日常生活的各个领域。特别是在冰川研究领域，通过分析数10 a尺度的Landsat数据，可以直观地测量冰川的演进速度。然而，由于轨道控制方面不够精准，应用Landsat-1～7系列卫星数据需要手工选择地面控制点，以实现精确的地面定位，形变的量算也只能通过肉眼对图像进行比较获得。Landsat-8解决了上述问题，为冰川科学研究提供了非常有价值的数据源。

在本研究中，我们通过亚像元级的图像匹配技术对轨道和帧号为145/35的15景Landsat-8影像数据进行互相关计算，并对其组成的时间序列进行SVD反演，获得了西昆仑地区目标冰川的二维位移场及速度场。

1数据源

研究区为西昆仑峰区,位于青藏高原西北部的新疆自治区境内，是我国冰川分布最密集的区域之一(图1(a))。

图1　研究区域Fig.1　The studied area

Landsat-8图像覆盖范围为185 km×185 km的区域，因此一景图像可以覆盖多条冰川。在一次图像互相关形变计算中，可以得到图像内所有冰川的水平位移场，因而具有很高的空间监测效率。同时，16 d的固定重访周期和不间断的数据获取，使得该数据的时间分辨率比其他卫星数据源具有更大优势。图1(b)为2013-09-11获取的Landsat-8影像；图1(c)显示了本研究中的目标冰川分布，即图1(b)中矩形框所示的范围。本研究将分析该景图像范围内的形变信息。

我们从USGS网站下载了2013-07～2014-08共15景覆盖西昆仑峰区的Landsat-8影像，构成形变研究的时间序列，数据信息如表1所示。为了提高估算冰川速度场的精度，我们只提取空间分辨率为15 m的全色波段图像进行图像互相关计算。同时，选择云量在30%以下的影像数据，以确保图像可以有效地精确配准。

表1　研究区 Landsat-8影像数据

2数据分析方法

从USGS网站下载的Landsat-8影像为使用SRTM DEM数据正射纠正过的图像。只有经过正射校正的图像，排除成像几何畸变，才能互相进行精确配准并分析形变。

图像亚像元互相关分析技术，最先用于SAR图像幅值信息的分析[5]，一方面用于计算SAR斜距向和方位向的位移场，另一方面用于SAR图像的精确配准。其基本原理是在两幅幅值图像上，按照一定间隔选取一定的配准网格点，在每个网格点附近按照一定的窗口大小搜索另一幅图像的相应位置，直到找到窗口内互相关系数最大的点，认为这两个图像互相关系数最大的点是地面同一个散射体目标。通过计算同一个目标在不同时间获取的SAR图像上的位置变化，就能知道该点地面发生的位移大小。为了实现高精度测量，SAR图像的计算都是在单视复图像上、雷达坐标系下进行的。而在光学图像的分析中，图像经过了地理编码和正射校正，所有形变计算直接在地理坐标系下进行。采用同样的方法，通过比较两幅图像的地面反射信息，可以实现亚像元级的图像精确匹配和位移场计算。

使用SAR图像配准方法，将两幅光学图像配准。一般配准精度(或称配准误差、偏差)在0.001～0.01像元范围内，可以达到理想的形变监测信噪比。在图像精确配准过程中，有必要控制获取影像的云量。 在配准之后，把一幅图像(称为“辅图像”或者“从图像”)重采样到另一幅图像(称为“主图像”)的几何坐标下[6]。图像重采样采用SINC插值方法。将任意两幅配准和重采样后的图像在COSI-CORR软件[7]中作互相关计算，可得到初步的位移结果。

2.1图像配准

配准步骤是整个形变计算的关键。使用图像亚像元互相关的方法，寻找两幅图像中共同的地面目标，用最小二乘的方法建立主从图像的多项式变换关系。配准和重采样过程消除了两幅图像的投影偏差，因而能够保证获得精确的地表形变信息[4]。另外一种配准方法是在主从图像之间建立映射关系，即查找表。查找表记录了由辅图像几何坐标到主图像的映射信息。建立查找表的同时，考虑地形的影响，可以进一步优化图像的配准过程。在建立了主从图像的查找表之后，继续用前面的配准多项式来校正查找表，并对其进行更新。基于更新后的查找表，重采样辅图像到参考图像坐标系下。经过该配准步骤，主从图像的配准精度可以进一步提升，配准偏差更小，因而更有利于形变的监测。在本研究中，我们选用获取时间较早的图像作为主图像来配准获取时间较晚的另一幅，当这幅图像精确配准之后，也可将它作为一幅新的主图像来配准其他从图像。以此类推，配准所有获得的Landsat-8数据。总的来说，较短时间间隔内的图像对，其配准结果相对较好，信噪比较高。

2.2图像互相关计算

地表水平位移场的获取是通过多时相正射影像图的像元互相关计算得到的。本研究采用的图像互相关方法是在傅里叶域里对图像进行反复迭代、无偏估计的过程[8]。这一过程已经在COSI-CORR 软件中实现，它将生成两幅互相关形变图像，分别是地表水平位移场的东西向和南北向分量；同时给出一幅评估形变监测结果精度的信噪比图像(图2)。图2给出的形变监测结果中，EW向形变不明显，而NS向位移场中可以看到几条明显的白色高值区域，显示几条最活跃冰川形成的地面位移场。由于这些冰川主要为近南北走向，与地形梯度方向一致，因而冰川位移也以南北向为主，其NS向的形变比较显著，这与本次观测的结果完全一致。整个图像互相关的信噪比SNR比较高，图像显示为高亮色。

图2　主辅图像互相关结果Fig.2　Cross-correlation products from a pair of master and slave images

2.3后处理

2.3.1条带去除

目前，光学图像一般通过多传感器机械扫描得到，这是为了避免单一光电传感器在扫描速率上的局限性[9]。例如，Landsat所使用的传感器系统，每一个频谱波段数据都由6台传感器同时扫描得到。因此，一次单相扫描能够产生图像上的6行，当下一次扫描时，卫星轨道适量前进，之前的6台传感器扫描生成下一个6行图像(LANDSAT 用户手册)。然而，这些传感器并没有统一的增益函数，因此导致互相关结果图上有明显的平行轨道的条带，并在垂直轨道方向起伏波动。例如，在ASTER和SPOT影像上出现过的类似振荡偏差或CCD阵列的错误振幅[10-11]。这种因传感器偏差导致的条带误差，可以通过COSI-CORR软件的后处理工具Destrip移除。

2.3.2数据滤波和异常值消除

互相关计算得到的形变监测结果必须通过滤波来降低噪声，但某些合理可信的结果常会被一些低相关系数的斑块所掩盖，必须先从形变场中移除一些过大或过小的异常值。首先计算位移值的概率密度分布，剔除低信噪比异常值，如图3所示(蓝色标识的部分认为是异常值，予以剔除)。

图3　位移值概率分布Fig.3　The probability distribution of a displacement field

我们选用基于非局部均值算法的滤波方法[12]对结果图像进行滤波处理，其在保证局部细节的同时降低了加性高斯白噪声。非局部均值算法是通过逐步调试来实现数据的滤波，具体方法是：对每个像元点，将其周围的小区域作为一个中心斑块，斑块大小一般为5×5或7×7像元。依据中心斑块与其他斑块的相似度来确定中心斑块周围更大范围的搜索区域，最终以该搜索区域内所有像元的平均值作为中心像元的滤波值。

3观测结果分析

3.1西昆仑峰区冰川位移场

在西昆仑山地区，大部分高山冰川呈南北走向，因此在互相关结果中，东西分量的水平位移信号相对较弱(图2)。所以，我们以南北向水平位移结果为例，研究西昆仑峰区冰川滑移及其随时间的演化。图4(a)给出了西昆仑峰区的地表南北向水平位移场，这是由2013-09-27和2013-12-16获取的两幅Landsat-8正射影像图作互相关得到的，其中一幅覆盖目标冰川的正射影像如图4(b)所示。该地区的位移最大值出现在图4(a)中的灰色矩形框内。在80 d内，两条目标冰川的滑移速率达到了1.5 m/d。此外，该地区内其他冰川上稍小的滑移速率也比较明显。图上红色表示冰川向南移动，蓝色表示向北移动。我们将滑移速率最大的两条目标冰川作为研究对象，其位移场在图4(c)中显示。

图4　西昆仑峰区的地表南北向水平位移场Fig.4　Horizontal displacement field of the glaciers in the west Kunlun mountains

图5　目标冰川位移场的时间序列演化Fig.5　The time-series plot of the displacement fields of the target glaciers

3.2 时间序列分析

我们处理分析了2013-07～2014-08的15景Landsat-8影像(表1)，并且利用奇异值分解(SVD)算法计算了目标冰川的年地表形变[13]。图5显示西昆仑山地区的冰川在观测时间段内，滑移量逐渐累积(南北向)。在384 d内，冰川位移随时间呈明显递增趋势。

为了确定目标冰川精确的滑移速率，我们在两条冰川上选取6个参考点(图4(c))，分析它们的位移(图6)和速率(图7)差异。由于图像处理过程中存在一些不可避免的随机噪声，导致利用奇异值分解算法反演位移场时间序列的结果存在部分“负增长”值(图6)，但依然可以明显看到，在384 d内位移场随时间近似呈线性递增，同时曲线的斜率在不同的参考点上是不同的，说明各点的滑动速率并不完全一致。图7给出了6个参考点的速率变化。可以看出，2013-11～2014-08冰川滑移速率保持在0.2～0.6 m/d，而对于每个参考点，冰川速度场呈匀速运动趋势。因此，研究区的这两条冰川始终保持了极为高速的滑动，其能够被Landsat-8这种图像分辨率的数据检测到。与其相比，InSAR形变监测方法，由于时间去相干效应的影响，监测如此高速的地面冰川滑动存在很大的局限性，这也反映了Landsat-8数据的观测优势。

图6　位移时间序列曲线Fig.6　The curves of the displacement time-series

图7　速度时间序列曲线Fig.7　The curves of the velocity time-series

4结语

本研究首次利用Landsat-8光学图像形变分析方法，获得西昆仑峰区冰川的形变场时间序列演化特征。基于光学影像资料提取地表形变，应用时间序列分析的方法、原理，对Landsat-8光学影像进行形变分析，获得西昆仑峰区冰川滑移的地表水平形变及其时间序列演化过程。结论如下：

1)光学影像互相关技术监测地表形变，特别是在监测较大形变和较快移动的地面目标时具有一定的优势。

2)基于幅值信息的影像配准方法、亚像元影像互相关形变计算、SVD反演时间序列——这些方法能够有效获得研究区时空形变场，提高测量精度。

3)在384 d观测期内，西昆仑峰区冰川位移场随时间近似呈线性变化，滑移速率在0.2～0.6 m/d。

4)本研究的成果不仅对冰川动力学和气候变化方面的研究有一定意义，而且验证了Landsat-8光学影像在形变分析中的应用潜力。该方法也适用于分析大地震和滑坡等形变场。

致谢：USGS为本研究提供了所有的Landsat-8影像数据,在此表示感谢!

参考文献

[1]焦克勤,姚檀栋,李世杰.西昆仑山32 Ka来的冰川与环境演变[J].冰川冻土,2000，22(3),250-256(Jiao Keqin,Yao Tandong,Li Shijie.Evolution of Glaciers and Environment in the West Kunlun Mountains During the Past 32 Ka[J].Journal of Glaciology and Geocryology,2000,22(3),250-256)

[2]李成秀,杨太保,田洪阵.1990～2011 年西昆仑峰区冰川变化的遥感监测[J].地理科学进展,2013,32(4)：548-559 (Li Chengxiu,Yang Taibao,Tian Hongzhen.Variation of West Kunlun Mountains Glacier during 1990-2011[J].Progress in Geography,2013,32(4):548-559)

[3]李成秀. 昆仑山冰川和积雪变化的遥感监测[D]. 兰州:兰州大学,2014(Li Chengxiu. Remote Sensing Monitoring of Glacier and Snow Cover Changes in the Kunlun Mountain[D]. Lanzhou :Lanzhou University,2014)

[4]Berthier E,Vadon H,Baratoux D,et al. Surface Motion of Mountain Glaciers Derived from Satellite Optical Imagery[J].Remote Sensing of Environment,2005,95(1):14-28

[5]刘方坚. 星载雷达影像配准方法的研究与实现[D]. 北京:中国地质大学(北京),2006(Liu Fangjian. The Registration Algorihtm Research of Spaceborne InSAR Images[D]. Beijing :China University of Geosciences(Beijing),2006)

[6]Massonnet D,Feigl K. Radar Interferometry and Its Application to Changes in the Earth’s Surface[J].Reviews of Geophysics,1998,36(4): 441-500

[7]Leprince S,Berthier E,Ayoub F,et al. Monitoring Earth Surface Dynamics with Optical Imagery[J].EOS Transactions,2008,89(1):1-12

[8]Scherler D,Leprince S,Strecker M R. Glacier-Surface Velocities in Alpine Terrain from Optical Satellite Imagery -Accuracy Improvement and Quality Assessment [J].Remote Sensing of Environment,2008,112(10):3 806-3 819

[9]Horn B K P, Woodham R J. Dest-Riping Landsat MSS Images by Histogram Modification[J].Computer Graphics and Image Processing,1979(10):69-83

[10] Leprince S,Muse P,Avouac J P. In-Flight CCD Distortion Calibration for Pushbroom Satellites Based on Subpixel Correlation[J]. IEEE Transactions on Geoscience and Remote Sensing,2008,46(9):2 675-2 683

[11] Teshima Y,Iwasaki A. Correction of Attitude Fluctuation of Terra Correction of Attitude Fluctuation of Terra with Parallax Observation[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(1):222-227

[12] Ayoub F,Leprince S,Keene L.User’s Guide to COSI-CORR[Z].Co-Registration of Optically Sensed Images and Correlation,2009

[13] Berardino P,Fornaro G,Lanari R,et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2 375-2 383

Foundation support:National Natural Science Foundation of China,No. 41374040.

About the first author:ZHENG Qian,postgraduate,majors in time-series analysis of SAR interferometry and crustal deformation application,E-mail: 359300096@qq.com.

Fast and Uniformly Slipping Western-Kunlun Glaciers from Time-series Deformation Analysis Using Periodically Captured Landsat-8 Imagery

ZHENGQian1SUNJianbao1ZHANGYong1

1State Key Laboratory of Earthquake Dynamics,Institute of Geology,CEA,A1 Huayanli,Beijing 100029,China

Abstract:We processed time-series data,using a sub-pixel cross-correlation analysis method, to obtain the temporal evolution of the surface displacements of the west Kunlun mountain glaciers in the north-western Tibetan plateau. This work co-registers periodically acquired Landsat-8 images using the sub-pixel cross-correlation technique at an accuracy of 0.01 pixels, implying a horizontal displacements detecting precision of the optical image 0.15 m. By cross-correlation analysis of the Landsat-8 images acquired from July 2013 to August 2014 and time-series inversion, this study constructs the displacement and velocity fields of two glaciers in the west Kunlun mountains. It demonstrates that the glaciers in the study area are slipping fast and uniformly (no obvious acceleration or deceleration). Based on Landsat-8 images, the method of this study can potentially be used for measuring ground deformation caused by dunes shifting, earthquakes, landslides, or volcanoes etc.

Key words:west Kunlun mountains; Landsat-8; co-registration; glacial motion; sub-pixel cross-correlation

收稿日期：2016-03-16

第一作者简介：郑茜，硕士生，研究方向为InSAR时间序列分析与地壳形变,E-mail: 359300096@qq.com;zhengqianabc@sina.com。

DOI：10.14075/j.jgg.2016.07.010

文章编号：1671-5942(2016)07-0604-05

中图分类号：P228

文献标识码：A

项目来源：国家自然科学基金(41374040)。