Three-dimensional movement inversion of Shisha Pangma glacier using Sentinel-1 and Landsat8 images
-
摘要: SAR偏移量技术和光学偏移量技术是冰川运动监测重要的技术手段，但目前对于融合不同平台的影像进行三维形变的研究较少。文章选取2019年11月至2021年1月西藏聂拉木县希夏邦马峰地区的大型冰川作为研究对象，基于方差分量估计融合该研究区的Sentinel-1与Landsat8两种数据进行冰川的三维位移解算，选取了同一时期的光学影像对偏移量估计结果进行对比分析，同时选取稳定区域进行精度评估，分析该方法在冰川运动监测中的适用性和精确性。结果表明，该冰川在2019年11月至2021年1月，联合解算的东西向最大流速为21 cm/d，南北向最大流速为68 cm/d，垂直向最大流速为17 cm/d。对比单一影像获取的冰川位移结果，多影像联合解算方法，能够弥补SAR数据的失相干和光学数据的低质量像元值的不足，获得更加完整和详细的冰川信息，监测结果精度更高。可为利用不同平台的数据联合监测山地冰川的多维度和高精度变化提供参考和技术支持。
-
关键词:
- Sentinel-1 /
- Landsat8 /
- 偏移量跟踪技术 /
- 联合解算 /
- 冰川三维监测
Abstract: SAR offset and optical offset are important techniques for glacier movement monitoring, but there are few researches on 3D deformation by integrating images from different platforms. In this paper, large glaciers in The Shisha Pangma area of Nyalam County, Tibet from November 2019 to January 2021 were selected as the research area. Based on variance component estimation, Sentinel-1 and Landsat8 data of the study area were combined for three-dimensional glacier calculation. Optical images of the same period were selected for comparative analysis of the offset estimation results. At the same time, the accuracy of the method was evaluated by selecting the stable area, and the applicability and accuracy of the method in glacier movement monitoring were analyzed. The results show that the maximum velocity of the glacier is 21 cm/day in the east-west direction, 68 cm/day in the north-south direction, and 17 cm/day in the vertical direction. Compared with the glacier displacement results obtained from a single image, the multi-image joint algorithm can compensate for the incoherence of SAR data and the low quality pixel values of optical data, and obtain more complete and detailed glacier information and higher accuracy of monitoring results. This paper can provide reference and technical support for using data from different platforms to jointly monitor multi-dimensional and high-precision changes of mountain glaciers. -
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 实验区数据覆盖图和冰川光学影像图Figure 1. Data coverage map of the experimental area and glacier optical images1.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影像编号 采集日期 空间基线/m 时间基线/d 1 2019-12-08 0.00 0 2 2020-01-25 −39.83 48 3 2020-03-13 −112.08 96 4 2020-04-30 −112.19 144 5 2020-06-17 −51.88 192 6 2020-08-04 −33.98 240 7 2020-09-21 −105.59 288 8 2020-11-08 −101.99 336 9 2020-12-26 125.79 384 表 2 研究所采用的Landsat8影像参数Table 2. The Landsat8 imaging parameters used in the study影像编号 采集日期 太阳高度角/(°) 太阳方位角/(°) 含云量/% 1 2019-11-19 38.52 157.40 1.74 2 2020-01-22 35.61 150.5 4.65 3 2020-02-07 39.10 147.2 3.28 4 2020-11-21 37.78 157.5 2.54 5 2021-01-08 33.93 153.0 4.20 6 2021-01-24 36.10 150.0 1.32 2. 研究方法与数据处理
文中对获取的Sentinel-1影像进行偏移量估计得到其距离向和方位向的形变量，对Landsat8影像进行偏移量估计得到东西向和南北向的形变量。为了融合这四个方向的形变量来求解冰川三维流速，操作步骤如下：(1)对距离向和方位向的形变位移结果进行投影转换和分辨率重采样，并将其分解至东西、南北和垂直向的位移；(2)根据观测误差的来源不同，与光学形变数据组成四个观测方程，基于方差分量估计对计算的SAR偏移量数据和光学偏移量数据进行随机模型的验后估计，从而得到冰川的三维形变速率；(3)对于解算的结果采用同期的光学影像对进行冰川流速及其运动特征的对比分析，并通过计算稳定区域的标准差进行精度估计。技术流程如图2所示，主要分为三个关键的步骤：SAR偏移量估计、光学偏移量估计、联合解算。
2.1 SAR偏移量估计
由于SAR影像易受失相干的影响，文中采用小基线集像素跟踪算法[21]，通过设置时空阈值得到最优的基线组合（图3），并利用对相干性要求较低的强度跟踪算法进行偏移量的估算，其主要的参数设置包括配准窗口、互相关系数和步长。匹配窗口的大小对影像配准的精度非常重要，采用多尺度窗口可以提高计算结果的精度[22]，获得更加精确的SAR影像距离向和方位向偏移结果。为了在偏移量信息和噪声值之间达到较好的平衡，经过多次实验确定最终的搜索窗口为128×128像素（距离向×方位向），搜索步长为5×1像素(距离向×方位向)，将相关系数阈值设置为0.2，以掩膜在失相干严重区域获得的不可靠位移值。
![]()
图 3 SAR影像时空基线组合分布图Figure 3. Space otemporal baseline combination distribution of SAR images2.2 光学偏移量估计
光学影像的偏移量估算结果易受失相关噪声、轨道误差、条带误差以及卫星姿态角误差的影响[23]，其中失相关噪声包括由云雪、地表建筑物以及地形阴影等。为了减少失相关噪声的影响，文中对获取的时间跨度小且含云量少的影像计算每一影像对的太阳高度角和太阳方位角差值，选取其差值较小的影像对（表3），共计7组影像对用于影像对的匹配。
表 3 研究所采用的Landsat8影像对Table 3. The Landsat8 image pairs used in the studyLandsat8影像对的日期 太阳高度角差值/(°) 太阳方位角差值/(°) 2019-11-19—2020-02-07 0.58 10.2 2019-11-19—2020-11-21 0.74 0.1 2020-01-22—2020-11-21 2.17 7.0 2020-01-22—2021-01-08 1.68 2.5 2020-01-22—2021-01-24 0.49 0.5 2020-02-07—2020-11-21 1.32 0.5 2020-11-21—2021-01-24 1.68 7.5 文中以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来说，假设发生了匀速形变，则有：
$$ {Los}^{ri} =[ {t}_{r}\cdot {C}_{e}^{ri}{t}_{r}\cdot {C}_{n}^{ri} {t}_{r}\cdot {C}_{u}^{ri}]\cdot {\left[{v}_{e}^{i}{v}_{n}^{i}{v}_{u}^{i}\right]}^{{{T}}} $$ (1) $$ {A\textit{z}}^{ri} =[ {t}_{r}\cdot {D}_{e}^{ri}{t}_{r}\cdot {D}_{n}^{ri} {t}_{r}\cdot {D}_{u}^{ri}]\cdot {\left[{v}_{e}^{i}{v}_{n}^{i}{v}_{u}^{i}\right]}^{{{T}}} $$ (2) 式中：
$ {Los}^{ri} $ 、${A\textit{z}}^{ri}$ ——第r个干涉对得到的相干点i在Los 向和方位向的形变量；${v}_{e}^{i}、{v}_{n}^{i}\mathrm{、}{v}_{u}^{i}$ ——i在东西、南北和垂直三个方向上的 形变速率；$ {t}_{r} $ ——第r个干涉对的时间间隔；${C}_{e}^{ri}$ 、${C}_{n}^{ri}$ 、${C}_{u}^{ri}$ ——Los向上相干点i在东西、南北和 垂直方向上的投影；$ {D}_{e}^{ri} $ 、$ {D}_{n}^{ri} $ 、$ {D}_{u}^{ri} $ ——方位向上相干点i在东西、南北和 垂直方向上的投影。其值根影的 几何关系（图4）来计算。$$\begin{aligned} &{C}_{e}^{ri} =-\sin{\theta }_{inc}^{i} \sin( {\alpha }_{a\textit{z}i}^{i}-3\pi /2 )\\ &{C}_{n}^{ri} =-\sin {\theta }_{inc}^{i} \cos( {\alpha }_{a\textit{z}i}^{i}-3\pi /2 )\\ &{C}_{u}^{ri} =\cos {\theta }_{inc}^{i} \\ &{D}_{e}^{ri} =-\cos( {\alpha }_{a\textit{z}i}^{i}-3\pi /2 )\\&{D}_{n}^{ri} =-\sin( {\alpha }_{a\textit{z}i}^{i}-3\pi /2 ) \\&{D}_{u}^{ri} =0 \end{aligned}$$ 式中：
$ {\theta }_{inc}^{i} $ 、${\alpha }_{a\textit{z}i}^{i}$ ——为相干点i的入射角和方位角。![]()
图 4 升轨SAR影像的几何关系(箭头的方向为正)Figure 4. Geometric relation of rail SAR image (arrow direction is positive)利用光学偏移量技术获取了f幅影像对，对第k幅影像对的相干点i来说，则有：
$$ {\left[{Q}_{e}^{ki}{Q}_{n}^{ki}\right]}^{{{T}}}={t}_{k}\cdot {\left[{v}_{e}^{i}{v}_{n}^{i}\right]}^{{{T}}} $$ (3) 式中：
$ {Q}_{e}^{ki} $ 、$ {Q}_{n}^{ki} $ ——相干点i在东西、南北方向的形变量。基于最小二乘模型融合SAR干涉对观测量和光学影像对观测量求解三维地表形变：
$$ L = BX + V $$ (4) 式中：X——求解的三维形变速率，
${{X}}$ =${\left[{v}_{e}^{i}\;\;{v}_{n}^{i}\;\;{v}_{u}^{i}\right]}^{{{T}}}$ ；L—2 m个SAR偏移量测量值和2f个光学偏移 量组成的观测量，
$ \underset{2（\text{m}+f）\times 1}{L} $ =[${Los}^{{{1i}}}\cdots$ $ {Los}^{mi} $ ${A{{\textit{z}}}}^{1i} \cdots$ ${A\textit{z}}^{mi}$ ${Q}_{e}^{1i} \cdots$ $ {Q}_{e}^{fi} $ ${Q}_{n}^{1i}\cdots$ $ {Q}_{n}^{fi} $ ]T；V——相对应的观测残差；
B——三个方向上的观测值组成的矩阵。
$$\begin{split} &\underset{2（\text{m}+f）\times 3}{B} =\\& {\left[\begin{array}{ccc}{t}_{1}\cdot {C}_{e}^{1i}\cdots {t}_{m}\cdot {C}_{e}^{mi}& {t}_{1}\cdot {D}_{e}^{1i}\cdots {t}_{m}\cdot {D}_{e}^{mi}& {t}_{1}\cdots {t}_{f}\\ {t}_{1}\cdot {C}_{n}^{1i}\cdots {t}_{m}\cdot {C}_{n}^{mi}& {t}_{1}\cdot {D}_{n}^{1i}\cdots {t}_{m}\cdot {D}_{n}^{mi}& 0\cdots 0\\ {t}_{1}\cdot {C}_{u}^{1i}\cdots {t}_{m}\cdot {C}_{u}^{mi}& {t}_{1}\cdot {D}_{u}^{1i}\cdots {t}_{m}\cdot {D}_{u}^{mi}& 0\cdots 0\end{array}\begin{array}{c}0\cdots 0\\ {t}_{1}\cdots {t}_{f}\\ 0\cdots 0\end{array}\right]}^{{{T}}} \end{split}$$ 假设观测值的方差是已知的，则利用最小二乘平差可以计算得到三维形变速率的最优估值：
$$ \hat X =( {{B}^{T}PB)}^{-1}{B}^{T}PL $$ (5) 其中，
$ P $ 为各个观测量方差组成的权阵：$$ \begin{aligned}\mathop P\limits_{2({m} + f) \times 2(m + f)} =& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(1/{{\sigma }}_{{Los}^{1{i}}}^{2}\cdots 1/{{\sigma }}_{{{L}{o}{s}}^{{m}{i}}}^{2}1/{{\sigma }}_{{{A}{\textit{z}}}^{1{i}}}^{2}\cdots 1/{{\sigma }}_{{{A}{\textit{z}}}^{{m}{i}}}^{2}1/{{\sigma }}_{{{Q}}_{{e}}^{1{i}}}^{2}\cdots\\& 1/{{\sigma }}_{{{Q}}_{{e}}^{{f}{i}}}^{2}1/{{\sigma }}_{{{Q}}_{{n}}^{1{i}}}^{2}\cdots 1/{{\sigma }}_{{{Q}}_{{n}}^{{f}{i}}}^{2}) . \end{aligned}$$ 要得到最优的三维形变速率估值，除函数模型外，还需要观测量的随机模型，也就是先验方差，就可以在平差模型中精确定权。通常观测值的方差往往很难精确获取，所以基于方差分量估计进行权阵的验后估计。由于SAR偏移量的观测误差受失相干、过采样和匹配窗口的影响，光学偏移量的观测误差受各种噪声、匹配窗口和步长等影响，因此根据观测值所受观测误差的不同进行分组，将SAR观测量的距离向和方位向各自分成一组，光学观测量的东西向和南北向各自分为一组，总共四组观测值。设每一组观测值为
$ {L}_{i} $ 、权重为$ {P}_{i} $ (i=1,2,3,4)，基于最小二乘方法得到第一次的估值：$$ \left\{\begin{aligned} &\hat X = {N^{ - 1}}W \\& {V_1} = {B_1}\hat X - {L_1} \\&{V_2} = {B_2}\hat X - {L_2} \\& {V_3} = {B_3}\hat X - {L_3} \\& {V_4} = {B_4}\hat X - {L_4}\end{aligned} \right.$$ (6) 其中：
$ L = {\left[ {\begin{array}{*{20}{c}} {{L_1}}&{{L_2}}&{{L_3}}&{{L_4}} \end{array}} \right]^T} $ ，$ B = {\left[ {\begin{array}{*{20}{c}} {{B_1}}&{{B_2}}&{{B_3}}&{{B_4}} \end{array}} \right]^T} $ ，$P = {\rm{diag}}{\left[ {\begin{array}{*{20}{c}} {{P_1}}&{{P_2}}&{{P_3}}&{{P_4}} \end{array}} \right]^T}$ ，$N = {B^T}PB = \displaystyle\sum\limits_{{{i}} = 1}^4 = {B_i^T} {P_i} {B_i} = $ $ \displaystyle\sum\limits_{{{i}} = 1}^4 {N_i^{}}$ ，$W = {B^T}PL = \displaystyle\sum\limits_{{{i}} = 1}^4 {B_i^T} {P_i}{L_i}$ 假设初始权阵
$ {P}_{i} $ 为单位阵，则方差分量和观测残差的关系：$$ \hat \theta = {S^{{- }1}}{W_\theta } $$ (7) 式中：
$\hat \theta = {\left[ {\begin{array}{*{20}{c}} {\hat \sigma _{{0_1}}^2}& {\hat \sigma _{{0_2}}^2}&{\hat \sigma _{{0_3}}^2}&{\hat \sigma _{{0_4}}^2} \end{array}} \right]^{{T}}}$ 为四组观测值的单位权中误差，${W_\theta } = {\left[ {\begin{array}{*{20}{c}} {V_1^T{P_1}{V_1}}&{V_2^T{P_2}{V_2}}&{V_3^T{P_3}{V_3}}&{V_4^T{P_4}{V_4}} \end{array}} \right]^{{T}}}$ ，并且$$ S = \left[ {\begin{array}{*{20}{c}} {{{m}} - 2tr({N^{ - 1}}{N_1}) + tr{{({N^{ - 1}}{N_1})}^2}}&{tr({N^{ - 1}}{N_1}{N^{ - 1}}{N_2})}&{tr({N^{ - 1}}{N_1}{N^{ - 1}}{N_3})}&{tr({N^{ - 1}}{N_1}{N^{ - 1}}{N_4})} \\ {tr({N^{ - 1}}{N_1}{N^{ - 1}}{N_2})}&{{{m}} - 2tr({N^{ - 1}}{N_2}) + tr{{({N^{ - 1}}{N_2})}^2}}&{tr({N^{ - 1}}{N_2}{N^{ - 1}}{N_3})}&{tr({N^{ - 1}}{N_2}{N^{ - 1}}{N_4})} \\ {tr({N^{ - 1}}{N_1}{N^{ - 1}}{N_3})}&{tr({N^{ - 1}}{N_2}{N^{ - 1}}{N_3})}&{{{f}} - 2tr({N^{ - 1}}{N_3}) + tr{{({N^{ - 1}}{N_3})}^2}}&{tr({N^{ - 1}}{N_3}{N^{ - 1}}{N_4})} \\ {tr({N^{ - 1}}{N_1}{N^{ - 1}}{N_4})}&{tr({N^{ - 1}}{N_2}{N^{ - 1}}{N_4})}&{tr({N^{ - 1}}{N_3}{N^{ - 1}}{N_4})}&{{{f}} - 2tr({N^{ - 1}}{N_4}) + tr{{({N^{ - 1}}{N_4})}^2}} \end{array}} \right] $$ 最后计算新的权阵估计值：
$$ {\hat P_1} = {P_1} \text{，} \hat P{}_2 = \frac{{\sigma _{{0_1}}^2}}{{\sigma _{{0_2}}^2P_2^{ - 1}}} \text{，} \hat P{}_3 = \frac{{\sigma _{{0_1}}^2}}{{\sigma _{{0_3}}^2P_3^{ - 1}}} \text{，} \hat P{}_4 = \frac{{\sigma _{{0_1}}^2}}{{\sigma _{{0_4}}^2P_4^{ - 1}}} $$ (8) 将得到的新权阵估计值代入到式(6)中计算最小二乘估计值
$ \hat X $ ，循环进行迭代式(6)到(8)，直到$$ \hat \sigma _{{0_1}}^2 \approx \hat \sigma _{{0_2}}^2 \approx \hat \sigma _{{0_3}}^2 \approx \hat \sigma _{{0_4}}^2 $$ (9) 当循环迭代结束时，此时的
$ {\hat P_1} $ 、$ {\hat P_2} $ 、$ {\hat P_3} $ 、$ {\hat P_4} $ 就是利用方差分量估计所得到的权阵，将其代入到式(6)中即可得到最优的三维形变速率。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影像对之间的距离向冰川表面位移分布Figure 5. Distance between the SAR image pairs is displaced to the glacial surface![]()
图 6 SAR影像对之间的方位向冰川表面位移分布Figure 6. Azimuthal to glacial surface displacement distribution between SAR image pairs3.1.2 光学偏移量冰川位移结果
为了获取目标冰川的光学偏移量观测值，利用同时期的6景Landsat8影像，计算每一影像对的太阳高度角和太阳方位角差值，选取差值较小的影像对组合成7对影像用于偏移量计算。基于COSI-Corr平台的频率域相关算法获得初始偏移量解算结果，冰川的南北向和东西向表面位移分布分别如图7和图8所示。同样，由于在误差处理过程中，将SNR<0.9的形变值剔除而导致了局部形变值的缺失。因此，与SAR偏移量结果融合可以进行信息互补。
![]()
图 7 Landsat8影像对之间的南北向冰川表面位移分布Figure 7. Distribution of north-south glacial surface displacement between Landsat8 image pairs![]()
图 8 Landsat8影像对之间的东西向冰川表面位移分布Figure 8. Distribution of east-west glacial surface displacement between Landsat8 image pairs为了验证太阳高度角和太阳方位角对结果精度的影响，选取光学偏移量计算结果中东西向和南北向稳定区域（见图1(b)的红色框Roi1和Roi2）的标准差进行精度评定。通过统计发现标准差与太阳高度角差值和太阳方位角差值之间的相互关系如图9所示，当太阳高度角差值和太阳方位角差值较小时，其标准差较低；当太阳高度角差值和太阳方位角差值较大时，其标准差较高。这也说明，在时间跨度较小的条件下，选择太阳高度角差值和太阳方位角差值较小的光学影像对能够提高偏移量估算结果的精度。
![]()
图 9 太阳高度角差值和太阳方位角差值与稳定区域标准差的关系Figure 9. Relationship between the solar height angular difference and the solar azimuth difference and the standard deviation of the stable region3.1.3 SAR和光学偏移量联合解算冰川位移结果
按照文中介绍的方差分量估计联合解算模型，将SAR与Landsat8两种数据源提取结果进行融合，得到了冰川同期的三维形变速率（图10）。结果显示，该冰川在东西向的最大形变速率为21 cm/d，南北向的最大形变速率为68 cm/d，垂直向的最大形变速率为17 cm/d。总体上呈现出该冰川上游的运动速度快，中游的速度比较慢，下游的速度逐渐减小。对于冰川上游三个分支的冰川流速变化较大且不均匀，并且南北向的最大流速达到68 cm/d，主要原因有：地形起伏变化大，下坡运动增大了冰川的流速；上游为冰川消融区域，在监测结果中存在特征点消失导致的局部信息不连续。
![]()
图 10 联合解算的冰川表面运动速率分布Figure 10. Jointly solved distribution of glacier surface motion rates为了更好地分析目标冰川的表面流速，将解算后的结果进行东西向、南北向和垂直向上的流速合成（图11）。总体上该冰川从西南方向往东北方向流动，从上游开始消融到中游汇合，与冰川的实际运动方向具有一致性，其中流速较大的地区包括上游冰川的消融区域，坡度较大的区域和三条上游冰川支流汇合的区域。
![]()
图 11 合成后的冰川表面运动速率分布Figure 11. Distribution of glacial surface motion rates after synthesis3.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月光学影像对间的冰川运动速率分布Figure 12. Glacial motion rate distribution between optical image pairs from November 2019 to January 2021为了进一步分析两种数据监测结果的一致性，按冰川的主流方向在上游和中游提取了这两种监测结果分别在南北向和东西向上冰川剖线MM′、NN′、PP′（图1）的流速剖线图（图13），分别对这三条剖线进行分析：
![]()
图 13 冰川剖线MM′、NN′、PP′在南北向和东西向上的流速和高程图注：(a)(b)分别为同期影像对与融合后MM′剖线的南北向和东西向流速；(c)(d)分别为同期影像对与融合后NN′剖线的南北向和东西向流速；(e)(f)分别为同期影像对与融合后PP′剖线的南北向和东西向流速。Figure 13. Flow velocity and elevation of MM′, NN′, PP′（1）MM′剖线：在南北向上，整体趋势具有一致性，同期光学影像对出现了因低质量点剔除的流速值见图13(a)绿色框，但在数据融合后出现了较大的流速变化可达到60cm/d，根据地形因素分析，流速增大的位置对应于较陡的斜坡，受地形的影响较大；在东西向上，整体的变化趋势较为一致，同样在融合后变化较大的地方见图13(b)绿色框，也是对应的斜坡位置。所以，融合后冰川流速的运动特征与地形起伏的变化相关。由于偏移量估算易受冰川表面特性变化的影响，对于同期影像对低质量的缺失值，通过融合后的数据得到了补偿，从图13中可以看出，融合后其对应的流速也较大。
（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统计值 融合后结果/(cm·d−1) Landsat8影像对结果/(cm·d−1) 南北向 东西向 南北向 东西向 Roi1 0.56 0.91 1.0 0.99 Roi2 0.48 0.60 0.77 0.84 4. 结论
多维、高精度的冰川运动变化监测能够为冰川灾害的预警和气候分析提供理论依据。因此，文中基于强度信息的SAR偏移量和基于亚像素频率域相关性匹配的光学偏移量方法，结合方差分量估计，解算了覆盖研究区冰川的三维流速，分析了冰川在不同方向的形变特征，并与同期光学影像对的计算结果进行对比分析，最后统计了稳定区域的标准差进行精度分析。
（1）不同平台影像偏移量融合可改善冰川位移监测的时间分辨率，提供更多的冰川位移细节信息。
（2）实验区冰川在融合后东西向的最大形变速率为21 cm/d，南北向的最大形变速率为68 cm/d，垂直向的最大形变速率为17 cm/d，其中流速较大的地区包括上游冰川的消融区域，坡度较大的区域和三条上游冰川支流汇合的区域。
（3）融合后的监测结果与同期光学影像对在南北向和东西向上具有较好的一致性，并且山地冰川的运动与坡度有关，陡坡处流速增大。
（4）融合后冰川位移监测结果的不确定性小于同期光学影像对的监测结果，在南北方向和东西方向上监测结果精度分别提升了38%和8%。
总体来说，同期单影像对解算方法虽然计算量小，但时间分辨率低，并且结果为二维形变。而不同平台影像的融合解算具有获取冰川的多维度形变信息，可进一步改善冰川监测的时间分辨率和获得更多形变特征的优势，但对于垂直向的研究及其精度分析需要进一步探索，因为垂直向的位移包含了冰川厚度的变化和沿斜坡的垂直向运动。基于不同平台解算多维度的形变信息研究可能是未来的发展趋势，本文可为利用不同平台的数据联合监测山地冰川的多维和高精度变化提供参考和技术支持。
-
![]()
图 1 实验区数据覆盖图和冰川光学影像图
Figure 1. Data coverage map of the experimental area and glacier optical images
![]()
图 3 SAR影像时空基线组合分布图
Figure 3. Space otemporal baseline combination distribution of SAR images
![]()
图 4 升轨SAR影像的几何关系(箭头的方向为正)
Figure 4. Geometric relation of rail SAR image (arrow direction is positive)
![]()
图 5 SAR影像对之间的距离向冰川表面位移分布
Figure 5. Distance between the SAR image pairs is displaced to the glacial surface
![]()
图 6 SAR影像对之间的方位向冰川表面位移分布
Figure 6. Azimuthal to glacial surface displacement distribution between SAR image pairs
![]()
图 7 Landsat8影像对之间的南北向冰川表面位移分布
Figure 7. Distribution of north-south glacial surface displacement between Landsat8 image pairs
![]()
图 8 Landsat8影像对之间的东西向冰川表面位移分布
Figure 8. Distribution of east-west glacial surface displacement between Landsat8 image pairs
![]()
图 9 太阳高度角差值和太阳方位角差值与稳定区域标准差的关系
Figure 9. Relationship between the solar height angular difference and the solar azimuth difference and the standard deviation of the stable region
![]()
图 10 联合解算的冰川表面运动速率分布
Figure 10. Jointly solved distribution of glacier surface motion rates
![]()
图 11 合成后的冰川表面运动速率分布
Figure 11. Distribution of glacial surface motion rates after synthesis
![]()
图 12 2019年11月—2021年1月光学影像对间的冰川运动速率分布
Figure 12. Glacial motion rate distribution between optical image pairs from November 2019 to January 2021
![]()
图 13 冰川剖线MM′、NN′、PP′在南北向和东西向上的流速和高程图
注：(a)(b)分别为同期影像对与融合后MM′剖线的南北向和东西向流速；(c)(d)分别为同期影像对与融合后NN′剖线的南北向和东西向流速；(e)(f)分别为同期影像对与融合后PP′剖线的南北向和东西向流速。
Figure 13. Flow velocity and elevation of MM′, NN′, PP′
表 1 研究所采用的Sentinel-1数据参数
Table 1 The Sentinel-1 data parameters used in the study
影像编号 采集日期 空间基线/m 时间基线/d 1 2019-12-08 0.00 0 2 2020-01-25 −39.83 48 3 2020-03-13 −112.08 96 4 2020-04-30 −112.19 144 5 2020-06-17 −51.88 192 6 2020-08-04 −33.98 240 7 2020-09-21 −105.59 288 8 2020-11-08 −101.99 336 9 2020-12-26 125.79 384 
下载: 导出CSV
表 2 研究所采用的Landsat8影像参数
Table 2 The Landsat8 imaging parameters used in the study
影像编号 采集日期 太阳高度角/(°) 太阳方位角/(°) 含云量/% 1 2019-11-19 38.52 157.40 1.74 2 2020-01-22 35.61 150.5 4.65 3 2020-02-07 39.10 147.2 3.28 4 2020-11-21 37.78 157.5 2.54 5 2021-01-08 33.93 153.0 4.20 6 2021-01-24 36.10 150.0 1.32
下载: 导出CSV
表 3 研究所采用的Landsat8影像对
Table 3 The Landsat8 image pairs used in the study
Landsat8影像对的日期 太阳高度角差值/(°) 太阳方位角差值/(°) 2019-11-19—2020-02-07 0.58 10.2 2019-11-19—2020-11-21 0.74 0.1 2020-01-22—2020-11-21 2.17 7.0 2020-01-22—2021-01-08 1.68 2.5 2020-01-22—2021-01-24 0.49 0.5 2020-02-07—2020-11-21 1.32 0.5 2020-11-21—2021-01-24 1.68 7.5
下载: 导出CSV
表 4 融合后和同期影像对稳定区域的标准差统计
Table 4 Standard deviation statistics for st able regions after fusion and concurrent images
统计值 融合后结果/(cm·d−1) Landsat8影像对结果/(cm·d−1) 南北向 东西向 南北向 东西向 Roi1 0.56 0.91 1.0 0.99 Roi2 0.48 0.60 0.77 0.84
下载: 导出CSV
-
[1] 姚晓军, 刘时银, 孙美平, 等. 20世纪以来西藏冰湖溃决灾害事件梳理[J]. 自然资源学报,2014,29(8):1377 − 1390. [YAO Xiaojun, LIU Shiyin, SUN Meiping, et al. Study on the glacial lake outburst flood events in Tibet since the 20th century[J]. Journal of Natural Resources,2014,29(8):1377 − 1390. (in Chinese with English abstract) [2] 童龙云, 张继, 孔应德. 西藏定日朋曲流域达仓沟冰湖溃决泥石流特征[J]. 中国地质灾害与防治学报,2019,30(6):34 − 39. [TONG Longyun, ZHANG Ji, KONG Yingde. Characteristics of the dacanggou debris flow induced by breakout of glacier-lake in the Pengqu basin Diri County of Tibet[J]. The Chinese Journal of Geological Hazard and Control,2019,30(6):34 − 39. (in Chinese with English abstract) [3] 王伟财. 藏东南伯舒拉岭地区冰湖变化及危险性与影响分析[D]. 北京: 中国科学院大学, 2013 WANG Weicai. Variation of glacial lakes and assessment of GLOF hazards in the Boshula mountain range, southeast of Tibetan Plateau[D]. Beijing: University of Chinese Academy of Sciences, 2013. (in Chinese with English abstract)
[4] 刘春玲, 童立强, 祁生文, 等. 喜马拉雅山地区冰川湖溃决灾害隐患遥感调查及影响因素分析[J]. 国土资源遥感,2016,28(3):110 − 115. [LIU Chunling, TONG Liqiang, QI Shengwen, et al. Remote sensing investigation and influence factor analysis of glacier lake outburst potential in the Himalayas[J]. Remote Sensing for Land & Resources,2016,28(3):110 − 115. (in Chinese with English abstract) [5] 刘建康, 张佳佳, 高波, 等. 我国西藏地区冰湖溃决灾害综述[J]. 冰川冻土,2019,41(6):1335 − 1347. [LIU Jiankang, ZHANG Jiajia, GAO Bo, et al. An overview of glacial lake outburst flood in Tibet, China[J]. Journal of Glaciology and Geocryology,2019,41(6):1335 − 1347. (in Chinese with English abstract) [6] GOLDSTEIN R M, ENGELHARDT H, KAMB B, et al. Satellite radar interferometry for monitoring ice sheet motion: Application to an Antarctic ice stream[J]. Science,1993,262(5139):1525 − 1530. DOI: 10.1126/science.262.5139.1525
[7] RIGNOT E, MOUGINOT J, SCHEUCHL B. Ice flow of the Antarctic ice sheet[J]. Science,2011,333(6048):1427 − 1430. DOI: 10.1126/science.1208336
[8] KUMAR V, VENKATARAMAN G, HØGDA K A, et al. Estimation and validation of glacier surface motion in the northwestern Himalayas using high-resolution SAR intensity tracking[J]. International Journal of Remote Sensing,2013,34(15):5518 − 5529. DOI: 10.1080/01431161.2013.792965
[9] STROZZI T, LUCKMAN A, MURRAY T, et al. Glacier motion estimation using SAR offset-tracking procedures[J]. IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2384 − 2391. DOI: 10.1109/TGRS.2002.805079
[10] 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):3806 − 3819. DOI: 10.1016/j.rse.2008.05.018
[11] DEHECQ A, GOURMELEN N, TROUVE E. Deriving large-scale glacier velocities from a complete satellite archive: Application to the Pamir-Karakoram-Himalaya[J]. Remote Sensing of Environment,2015,162:55 − 66. DOI: 10.1016/j.rse.2015.01.031
[12] 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. DOI: 10.1016/j.rse.2004.11.005
[13] DEHECQ A, GOURMELEN N, TROUVE E. Deriving large-scale glacier velocities from a complete satellite archive: Application to the Pamir-Karakoram-Himalaya[J]. Remote Sensing of Environment,2015,162:55 − 66. DOI: 10.1016/j.rse.2015.01.031
[14] HU J, LI Z W, DING X L, et al. 3D coseismic Displacement of 2010 Darfield, New Zealand earthquake estimated from multi-aperture InSAR and D-InSAR measurements[J]. Journal of Geodesy,2012,86(11):1029 − 1041. DOI: 10.1007/s00190-012-0563-6
[15] GRAY L. Using multiple RADARSAT InSAR pairs to estimate a full three-dimensional solution for glacial ice movement[J]. Geophysical Research Letters,2011,38(5).
[16] MINCHEW B M, SIMONS M, RIEL B, et al. Tidally induced variations in vertical and horizontal motion on Rutford Ice Stream, West Antarctica, inferred from remotely sensed observations[J]. Journal of Geophysical Research:Earth Surface,2017,122(1):167 − 190. DOI: 10.1002/2016JF003971
[17] HU J, LI Z W, DING X L, et al. Two-dimensional co-seismic surface displacements field of the Chi-Chi earthquake inferred from SAR image matching[J]. Sensors (Basel, Switzerland),2008,8(10):6484 − 6495. DOI: 10.3390/s8106484
[18] LI J, LI Z W, WU L X, et al. Deriving a time series of 3D glacier motion to investigate interactions of a large mountain glacial system with its glacial lake: use of Synthetic Aperture Radar Pixel Offset-Small Baseline Subset technique[J]. Journal of Hydrology,2018,559:596 − 608. DOI: 10.1016/j.jhydrol.2018.02.067
[19] 敖萌, 张路, 廖明生, 等. 基于方差分量估计的多源InSAR数据自适应融合形变测量[J]. 地球物理学报,2020,63(8):2901 − 2911. [AO Meng, ZHANG Lu, LIAO Mingsheng, et al. Deformation monitoring with adaptive integration of multi-source InSAR data based on variance component estimation[J]. Chinese Journal of Geophysics,2020,63(8):2901 − 2911. (in Chinese with English abstract) [20] 李海, 杨成生, 惠文华, 等. 基于遥感技术的高山极高山区冰川冰湖变化动态监测: 以西藏藏南希夏邦玛峰地区为例[J]. 中国地质灾害与防治学报,2021,32(5):10 − 17. [LI Hai, YANG Chengsheng, HUI Wenhua, et al. Changes of glaciers and glacier lakes in alpine and extremely alpine regions using remote sensing technology: a case study in the Shisha Pangma area of southern Tibet[J]. The Chinese Journal of Geological Hazard and Control,2021,32(5):10 − 17. (in Chinese with English abstract) [21] 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):2375 − 2383. DOI: 10.1109/TGRS.2002.803792
[22] 杨丽叶, 赵超英, 杨成生. PO-SBAS技术用于错朗玛冰川三维时序运动特征分析[J]. 地球物理学进展,2020,35(6):2116 − 2123. [YANG Liye, ZHAO Chaoying, YANG Chengsheng. Three-dimensional time series movement analysis of the Cuolangma glacier with PO-SBAS technique[J]. Progress in Geophysics,2020,35(6):2116 − 2123. (in Chinese with English abstract) [23] 贺礼家, 冯光财, 冯志雄, 等. 哨兵-2号光学影像地表形变监测: 以2016年MW7.8新西兰凯库拉地震为例[J]. 测绘学报,2019,48(3):339 − 351. [HE Lijia, FENG Guangcai, FENG Zhixiong, et al. Coseismic displacements of 2016 MW7.8 Kaikoura, New Zealand earthquake, using Sentinel-2 optical images[J]. Acta Geodaetica et Cartographica Sinica,2019,48(3):339 − 351. (in Chinese with English abstract) [24] 柳林, 宋豪峰, 杜亚男, 等. 联合哨兵2号和Landsat 8估计白格滑坡时序偏移量[J]. 武汉大学学报·信息科学版,2021,46(10):1461 − 1470. [LIU Lin, SONG Haofeng, DU Yanan, et al. Time-series offset tracking of the Baige landslide based on Sentinel-2 and Landsat8[J]. Geomatics and Information Science of Wuhan University,2021,46(10):1461 − 1470. (in Chinese with English abstract) [25] 胡俊. 基于现代测量平差的InSAR三维形变估计理论与方法[D]. 长沙: 中南大学, 2013 HU Jun. Theory and method of estimating three-dimensional displacement with InSAR based on the modern surveying adjustment[D]. Changsha: Central South University, 2013. (in Chinese with English abstract)
-
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 323
- HTML全文浏览量: 145
- PDF下载量: 243
- 被引次数: 1
邮件订阅
RSS