联合Sentinel-1与Landsat8影像的希夏邦马峰冰川三维运动反演

魏春蕊，杨成生，魏云杰，熊国华，李晓阳

（1. 长安大学地质工程与测绘学院，陕西 西安 710054；2. 中国地质环境监测院（自然资源部地质灾害技术指导中心），北京 100081）

0 引言

全球气候变暖导致冰川和冻土消融，进而引发滑坡、泥石流等自然灾害。西藏是我国山地冰川分布最广泛的区域，据资料统计，20世纪以来西藏地区发生过多次冰湖溃决事件[1−2]。1991年由冰川活动诱发西藏八宿县吉达乡江查村各布马错冰湖溃决引发泥石流，冲毁10余处川藏公路路基、55 ha农田，冲走11户民房、144头牲畜、18座木桥[3]；2000年位于康马县萨马达乡雅拉山口北坡的冲巴吓错湖因后缘冰川滑坡引发涌浪溃坝洪水和泥石流，使冰湖下游的萨马达、康马、少岗和南尼等4个乡受灾[4]；2014年由于冰川泥石流导致波密县的那隆藏布支沟冰湖溃决造成易贡乡2.7ha农田和550 m乡村道路被洪水冲毁等[5]。冰川运动变化的精确监测能够为冰川灾害的预警和气候分析提供理论依据，具有重大意义。

近年来，偏移量跟踪技术在冰川位移的研究中得到了广泛应用[6−7]。目前，偏移量跟踪技术根据影像来源不同可分为两大类：SAR偏移量技术、光学偏移量技术。由于SAR影像的获取不受天气条件的限制，且测量周期短，监测范围广，许多国内外学者将SAR偏移量跟踪技术应用到测量冰川流速的研究中[8−9]。然而SAR影像侧视成像的特点使结果受叠掩和阴影等几何畸变的影响较大，且计算结果为距离向和方位向形变，因此一些学者也集中研究基于光学影像的像素偏移量跟踪技术[10−11]。例如，BERTHIER等[12]利用SPOT5光学影像获取了阿尔卑斯地区山地冰川流速场[12]，DEHECQ等[13]利用Landsat-5和7影像，获取了帕米尔-喀喇昆仑-喜马拉雅(PKH)三年期间的冰川年速度场。

随着多源数据和多平台影像的发展，针对大量级的位移监测，也有学者相继提出多孔径差分干涉测量技术(MAI)[14]、融合InSAR升降轨技术的冰川二维或三维运动监测[15]，融合多个轨道的Offset-tracking或MAI形变测量技术[16−17]、小基线时序处理的PO-SBAS技术[18]和融合多源InSAR数据的三维形变测量[19]，但目前对于融合SAR影像与光学影像监测多维度形变的研究较少。因此，文中选取西藏聂拉木县希夏邦马峰地区的大型冰川作为实验区，基于方差分量估计，利用Sentinel-1与Landsat8两种影像进行了冰川位移的三维分解研究，分析该方法在冰川运动监测中的适用性和精确性。

1 研究区与数据

1.1 研究区概况

希夏邦马峰为喜马拉雅山脉中段的高峰之一，空间位置如图1(a)所示，是一座完全在中国西藏自治区聂拉木县境内的8 000 m级高峰。希夏邦玛峰地区山地冰川持续退缩，并且冰湖加速扩张，该地区冰湖溃决等自然灾害[20]增多。文中的研究对象属于希夏邦马峰的大型冰川，如图1(b)所示，该冰川有3条上游支流并且在中游汇集，为冰川三维运动特征分析提供了典型实验区域。

1.2 研究区数据

文中采用了Sentinel-1和Landsat8两种数据源。考虑到计算的高效率以及冰川在较长时间间隔里其位移量较明显的因素，采用了从欧空局获取的间隔48 d的2019年12月至2020年12月9景Sentinel-1升轨影像，详细参数信息如表1所示，其距离向分辨率为5 m、方位向分辨率为20 m。对于光学数据，采用了2019年11月至2021年1月之间的6景美国Landsat8影像，参数信息如表2所示，Landsat8的全色影像空间分辨率为15 m。这两种数据的影像覆盖范围如图1所示。

表1 研究所采用的Sentinel-1数据参数Table 1 The Sentinel-1 data parameters used in the study

表2 研究所采用的Landsat8影像参数Table 2 The Landsat8 imaging parameters used in the study

图1 实验区数据覆盖图和冰川光学影像图Fig.1 Data coverage map of the experimental area and glacier optical images

2 研究方法与数据处理

文中对获取的Sentinel-1影像进行偏移量估计得到其距离向和方位向的形变量，对Landsat8影像进行偏移量估计得到东西向和南北向的形变量。为了融合这四个方向的形变量来求解冰川三维流速，操作步骤如下：(1)对距离向和方位向的形变位移结果进行投影转换和分辨率重采样，并将其分解至东西、南北和垂直向的位移；(2)根据观测误差的来源不同，与光学形变数据组成四个观测方程，基于方差分量估计对计算的SAR偏移量数据和光学偏移量数据进行随机模型的验后估计，从而得到冰川的三维形变速率；(3)对于解算的结果采用同期的光学影像对进行冰川流速及其运动特征的对比分析，并通过计算稳定区域的标准差进行精度估计。技术流程如图2所示，主要分为三个关键的步骤：SAR偏移量估计、光学偏移量估计、联合解算。

图2 联合解算冰川三维形变速率技术流程图Fig.2 Technical flow chart of 3D deformation rate

2.1 SAR偏移量估计

由于SAR影像易受失相干的影响，文中采用小基线集像素跟踪算法[21]，通过设置时空阈值得到最优的基线组合（图3），并利用对相干性要求较低的强度跟踪算法进行偏移量的估算，其主要的参数设置包括配准窗口、互相关系数和步长。匹配窗口的大小对影像配准的精度非常重要，采用多尺度窗口可以提高计算结果的精度[22]，获得更加精确的SAR影像距离向和方位向偏移结果。为了在偏移量信息和噪声值之间达到较好的平衡，经过多次实验确定最终的搜索窗口为128×128像素（距离向×方位向），搜索步长为5×1像素(距离向×方位向)，将相关系数阈值设置为0.2，以掩膜在失相干严重区域获得的不可靠位移值。

图3 SAR影像时空基线组合分布图Fig.3 Space otemporal baseline combination distribution ofSAR images

2.2 光学偏移量估计

光学影像的偏移量估算结果易受失相关噪声、轨道误差、条带误差以及卫星姿态角误差的影响[23]，其中失相关噪声包括由云雪、地表建筑物以及地形阴影等。为了减少失相关噪声的影响，文中对获取的时间跨度小且含云量少的影像计算每一影像对的太阳高度角和太阳方位角差值，选取其差值较小的影像对（表3），共计7组影像对用于影像对的匹配。

表3 研究所采用的Landsat8影像对Table 3 The Landsat8 image pairs used in the study

文中以COSI-Corr软件为偏移量估算平台，基于亚像素相关性匹配算法获取影像对的东西向和南北向形变量。为了得到较好的数据处理结果，进行多次实验选择合适的参数：初始搜索窗口设置为64×64，最终相关搜索窗口设置为32×32，步长设置为4，迭代次数为2，掩膜阈值为0.9。处理过程中，采用了以下措施来削弱或消除误差的影响：设置信噪比阈值>0.9去除失相关噪声引起的误差；利用一次多项式曲面拟合模型去除轨道误差；利用“均值相减法”去除条带误差；利用改进后的“均值相减法”去除卫星姿态角误差[24]。

2.3 联合解算

在多类观测数据的联合平差中，Helmert方差分量估计是对不同类型的观测值进行定权的常用方法之一。它的基本思想：将观测值按不同的观测来源分类，确定各类观测值的初始权，进行预平差；然后根据一定的原则计算各类观测值的验后单位权方差，进行观测量方差的迭代运算，直至各类观测值的单位权中误差相等或各类单位权方差之比等于1。之前有学者已经用该方法探索过InSAR观测量与GPS数据的融合[25]。同样，文中基于方差分量估计对计算的SAR偏移量数据和光学偏移量数据进行随机模型的验后估计，具体步骤和原理：

首先，以Landsat8数据为参考对象，将SAR数据的坐标系统转换到Landsat8数据的坐标系统下，并对SAR干涉对重采样进行分辨率的统一，使得两种图像的像元能够一一对应，便于后续的计算。

假设有N幅研究区的SAR影像组成了m个干涉对，对第r（r≤m）个影像对的相干点i来说，假设发生了匀速形变，则有：

式中：Losri、Azri——第r个干涉对得到的相干点i在Los向和方位向的形变量；

在东西、南北和垂直三个方向上的形变速率；

tr——第r个干涉对的时间间隔；

——Los向上相干点i在东西、南北和垂直方向上的投影；

——方位向上相干点i在东西、南北和垂直方向上的投影。其值根影的几何关系（图4）来计算。

图4 升轨SAR影像的几何关系(箭头的方向为正)Fig.4 Geometric relation of rail SAR image (arrow direction is positive)

利用光学偏移量技术获取了f幅影像对，对第k幅影像对的相干点i来说，则有：

基于最小二乘模型融合SAR干涉对观测量和光学影像对观测量求解三维地表形变：

式中：X——求解的三维形变速率，

L—2m个SAR偏移量测量值和2f个光学偏移量组成的观测量，L=[Los1i···LosmiAz1i···

V——相对应的观测残差；

B——三个方向上的观测值组成的矩阵。

假设观测值的方差是已知的，则利用最小二乘平差可以计算得到三维形变速率的最优估值：

其中，P为各个观测量方差组成的权阵：

要得到最优的三维形变速率估值，除函数模型外，还需要观测量的随机模型，也就是先验方差，就可以在平差模型中精确定权。通常观测值的方差往往很难精确获取，所以基于方差分量估计进行权阵的验后估计。由于SAR偏移量的观测误差受失相干、过采样和匹配窗口的影响，光学偏移量的观测误差受各种噪声、匹配窗口和步长等影响，因此根据观测值所受观测误差的不同进行分组，将SAR观测量的距离向和方位向各自分成一组，光学观测量的东西向和南北向各自分为一组，总共四组观测值。设每一组观测值为Li、权重为Pi(i=1,2,3,4)，基于最小二乘方法得到第一次的估值：

假设初始权阵Pi为 单位阵，则方差分量和观测残差的关系：

最后计算新的权阵估计值：

3 结果与分析

3.1 冰川解算结果

3.1.1 SAR偏移量冰川位移结果

对于覆盖研究实验区冰川的9景Sentinel-1影像，基于小基线集的思想得到了14对SAR影像的优化组合（图3），利用SAR偏移量强度跟踪技术获取了目标冰川每对影像的距离向位移量（图5）和方位向位移量（图6）。从目标冰川的SAR偏移量组合对结果中可以得出，在2020年，1−3月整体形变量较小，4−8月的形变量增加，9月份形变量减缓，11−12月形变量较小，这也符合冰川季节性的变化趋势。对于图中形变值缺失的部分，由于SAR影像的失相干特性，在计算过程中掩膜了低相干区域的值。因此，形变值缺失的区域部分可以通过与光学偏移量的融合得到形变信息的补偿。

图5 SAR影像对之间的距离向冰川表面位移分布Fig.5 Distance between the SAR image pairs is displaced to the glacial surface

图6 SAR影像对之间的方位向冰川表面位移分布Fig.6 Azimuthal to glacial surface displacement distribution between SAR image pairs

3.1.2 光学偏移量冰川位移结果

为了获取目标冰川的光学偏移量观测值，利用同时期的6景Landsat8影像，计算每一影像对的太阳高度角和太阳方位角差值，选取差值较小的影像对组合成7对影像用于偏移量计算。基于COSI-Corr平台的频率域相关算法获得初始偏移量解算结果，冰川的南北向和东西向表面位移分布分别如图7和图8所示。同样，由于在误差处理过程中，将SNR<0.9的形变值剔除而导致了局部形变值的缺失。因此，与SAR偏移量结果融合可以进行信息互补。

图7 Landsat8影像对之间的南北向冰川表面位移分布Fig.7 Distribution of north-south glacial surface displacement between Landsat8 image pairs

图8 Landsat8影像对之间的东西向冰川表面位移分布Fig.8 Distribution of east-west glacial surface displacement between Landsat8 image pairs

为了验证太阳高度角和太阳方位角对结果精度的影响，选取光学偏移量计算结果中东西向和南北向稳定区域（见图1(b)的红色框Roi1和Roi2）的标准差进行精度评定。通过统计发现标准差与太阳高度角差值和太阳方位角差值之间的相互关系如图9所示，当太阳高度角差值和太阳方位角差值较小时，其标准差较低；当太阳高度角差值和太阳方位角差值较大时，其标准差较高。这也说明，在时间跨度较小的条件下，选择太阳高度角差值和太阳方位角差值较小的光学影像对能够提高偏移量估算结果的精度。

图9 太阳高度角差值和太阳方位角差值与稳定区域标准差的关系Fig.9 Relationship between the solar height angular difference and the solar azimuth difference and the standard deviation of the stable region

3.1.3 SAR和光学偏移量联合解算冰川位移结果

按照文中介绍的方差分量估计联合解算模型，将SAR与Landsat8两种数据源提取结果进行融合，得到了冰川同期的三维形变速率（图10）。结果显示，该冰川在东西向的最大形变速率为21 cm/d，南北向的最大形变速率为68 cm/d，垂直向的最大形变速率为17 cm/d。总体上呈现出该冰川上游的运动速度快，中游的速度比较慢，下游的速度逐渐减小。对于冰川上游三个分支的冰川流速变化较大且不均匀，并且南北向的最大流速达到68 cm/d，主要原因有：地形起伏变化大，下坡运动增大了冰川的流速；上游为冰川消融区域，在监测结果中存在特征点消失导致的局部信息不连续。

图10 联合解算的冰川表面运动速率分布Fig.10 Jointly solved distribution of glacier surface motion rates

为了更好地分析目标冰川的表面流速，将解算后的结果进行东西向、南北向和垂直向上的流速合成（图11）。总体上该冰川从西南方向往东北方向流动，从上游开始消融到中游汇合，与冰川的实际运动方向具有一致性，其中流速较大的地区包括上游冰川的消融区域，坡度较大的区域和三条上游冰川支流汇合的区域。

图11 合成后的冰川表面运动速率分布Fig.11 Distribution of glacial surface motion rates after synthesis

3.2 结果对比

为了检验联合模型解算的冰川位移及其运动特征的可靠性，对联合解算的南北向和东西向的冰川流速与同一时期光学影像对的冰川流速进行对比分析。这里选择光学影像对的原因是，光学偏移量的结果与联合解算后的南北和东西向为一致的方向，具有对比性，而SAR偏移量的结果为距离向和方位向，无可对比性。因此，选择 2019年11月19日与2021年1月24日的Landsat8影像对，从表2可知该影像对的太阳高度角差值和太阳方位角差值较小，其差值分别为2.42°和7.4°，结合图9分析得到的结果具有较高的精度。在进行误差消除处理和去除低质量点后，对影像对进行了偏移量估计，得到该冰川南北向和东西向的表面运动速率。如图12所示，从整体上看，冰川在南北向和东西向的运动速率达到了60 cm/d，并且运动速率较大的区域为冰川上游消融区和中游支流汇集区，与联合结算后的结果与图10(a)、图10(b)较为符合。由于光学影像偏移量技术不能获取垂直向位移，且没有其他监测数据可以被利用，因此我们无法对冰川的垂直向位移进行对比。

图12 2019年11月—2021年1月光学影像对间的冰川运动速率分布Fig.12 Glacial motion rate distribution between optical image pairs from November 2019 to January 2021

为了进一步分析两种数据监测结果的一致性，按冰川的主流方向在上游和中游提取了这两种监测结果分别在南北向和东西向上冰川剖线MM′、NN′、PP′（图1）的流速剖线图（图13），分别对这三条剖线进行分析：

（1）MM′剖线：在南北向上，整体趋势具有一致性，同期光学影像对出现了因低质量点剔除的流速值见图13(a)绿色框，但在数据融合后出现了较大的流速变化可达到60cm/d，根据地形因素分析，流速增大的位置对应于较陡的斜坡，受地形的影响较大；在东西向上，整体的变化趋势较为一致，同样在融合后变化较大的地方见图13(b)绿色框，也是对应的斜坡位置。所以，融合后冰川流速的运动特征与地形起伏的变化相关。由于偏移量估算易受冰川表面特性变化的影响，对于同期影像对低质量的缺失值，通过融合后的数据得到了补偿，从图13中可以看出，融合后其对应的流速也较大。

图13 冰川剖线MM′、NN′、PP′在南北向和东西向上的流速和高程图Fig.13 Flow velocity and elevation of MM′, NN′, PP′

（2）NN′剖线：从图13(c)和图13(d)可以看出，融合后与同期影像对的结果整体趋势上的变化是一致的，在南北向上，流速变化较为平缓；在东西向上，在坡度较陡的位置冰川流速较大。根据地形因素分析，冰川的运动与坡度的变化有关。

（3）PP′剖线：在冰川的中游区域，从图13(e)和图13(f)可以看出，总体趋势具有一致性，并且地势起伏变化小，同期影像对与融合后的结果没有出现较大的流速变化差异，微小的差异是因为研究方法的不同。对于图13(f)中东西向凸起的流速变化，是因为冰川的三条支流汇集增大了流速的变化。

3.3 精度分析

由于缺乏实测数据，采用先验信息辅助进行精度分析是常用的方法，比如利用稳定区域的标准差。由于光学偏移量的结果与联合解算后的南北和东西方向一致，具有可对比性。因此，对于两种不同方法的冰川流速监测结果，我们统计了融合后与同期影像对在稳定区域见图1(b)的红色框Roi1和Roi2的标准差（表4）。在南北方向上，Landsat8影像对在稳定区域的统计结果分别为1.0 cm/d和0.77 cm/d，而融合解算结果在稳定区域的标准差分别为0.56 cm/d和0.48 cm/d，相比于同期光学影像对结果提升了44%和38%；在东西方向上，Landsat8影像对在稳定区域的统计结果分别为0.99 cm/d和0.84 cm/d，而融合解算结果在稳定区域的标准差分别为0.91 cm/d和0.6 cm/d，相比于同期影像对结果提升了8%和29%。总体而言，融合后冰川位移监测结果的不确定性小于同期光学影像对的监测结果。

表4 融合后和同期影像对稳定区域的标准差统计Table 4 Standard deviation statistics for st able regions after fusion and concurrent images

4 结论

多维、高精度的冰川运动变化监测能够为冰川灾害的预警和气候分析提供理论依据。因此，文中基于强度信息的SAR偏移量和基于亚像素频率域相关性匹配的光学偏移量方法，结合方差分量估计，解算了覆盖研究区冰川的三维流速，分析了冰川在不同方向的形变特征，并与同期光学影像对的计算结果进行对比分析，最后统计了稳定区域的标准差进行精度分析。

（1）不同平台影像偏移量融合可改善冰川位移监测的时间分辨率，提供更多的冰川位移细节信息。

（2）实验区冰川在融合后东西向的最大形变速率为21 cm/d，南北向的最大形变速率为68 cm/d，垂直向的最大形变速率为17 cm/d，其中流速较大的地区包括上游冰川的消融区域，坡度较大的区域和三条上游冰川支流汇合的区域。

（3）融合后的监测结果与同期光学影像对在南北向和东西向上具有较好的一致性，并且山地冰川的运动与坡度有关，陡坡处流速增大。

（4）融合后冰川位移监测结果的不确定性小于同期光学影像对的监测结果，在南北方向和东西方向上监测结果精度分别提升了38%和8%。

总体来说，同期单影像对解算方法虽然计算量小，但时间分辨率低，并且结果为二维形变。而不同平台影像的融合解算具有获取冰川的多维度形变信息，可进一步改善冰川监测的时间分辨率和获得更多形变特征的优势，但对于垂直向的研究及其精度分析需要进一步探索，因为垂直向的位移包含了冰川厚度的变化和沿斜坡的垂直向运动。基于不同平台解算多维度的形变信息研究可能是未来的发展趋势，本文可为利用不同平台的数据联合监测山地冰川的多维和高精度变化提供参考和技术支持。