2025, Vol. 42
2. 广东工业大学 土木与交通工程学院, 广东 广州 510006;
3. 广东工业大学 大湾区城市环境安全与绿色发展教育部重点实验室, 广东 广州 510006;
4. 华南农业大学 资源环境学院, 广东 广州 510642
2. School of Civil and Transportation Engineering, Guangdong University of Technology, Guangzhou 510006, China;
3. Key Laboratory for City Cluster Environmental Safety and Green Development of the Ministry of Education, Guangdong University of Technology, Guangzhou 510006, China;
4. College of Natural Resources and Environment, South China Agricultural University, Guangzhou 510642, China
滑坡会破坏周围的环境，给社会带来巨大的经济损失，同时也会严重影响人民的生命安全[1-2]。因此，非常有必要对滑坡易发生区域进行长期监测，用来作为滑坡灾害早期预警的参考。滑坡往往发生在地形起伏比较大，植被覆盖茂密的山区，对其进行大规模的地面监测非常困难。而基于遥感技术的滑坡监测具有范围大、精度高等优势，因而被广泛应用于滑坡的时序监测[3-5]。
合成孔径雷达对地表形变监测的方法一般分为两种：合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar, InSAR) 技术和偏移追踪(Offset Tracking, OT) 技术。其中InSAR在精度方面比较突出，但容易受到时空失相干的影响，适合于缓慢的滑坡[6-7]。偏移追踪技术虽然在精度上比不上InSAR技术，但受相干性的限制较小，可以测量更大范围的梯度形变(米级)，因此偏移追踪技术更适合于中速以上的滑坡监测[8-10]。传统偏移追踪技术(Traditional Offset Tracking, Trad-OT) 是在整个研究区域内以一定的间隔规则地选择同名点，使用基于规则匹配窗口的归一化互相关跟踪方法，找出主影像与从影像中的相同点，以此来计算两个影像的偏移量[11-13]。许多研究表明，在光学遥感中Trad-OT可获取大尺度的地表形变[14-15]，但Trad-OT仍然存在一些局限性。若规则匹配窗口的像素包括两种或多种类型的移动特征：一部分像素具有大位移，另外一部分像素具有小位移或无位移。这些像素是不相关的，将会极大地降低自相关计算结果的可靠性，从而导致偏移量估计出错。为了解决这个问题，Cai等[16]提出了一种自适应偏移追踪算法，通过设计自适应窗口，使得窗口内仅包括与窗口中心像素运动特性相似的像素，并将该方法应用在了滑坡边界区域，获得了更加准确的形变信息。
在光学遥感影像分析中，偏移追踪技术通常会受到水体、阴影和云层等干扰因素的影响，从而导致偏移估计出错。为了解决这个问题，本文提出了一种自适应偏移跟踪算法。在偏移估计前，先进行预处理，识别出研究区域中水体、阴影和云层等干扰因素的位置，并生成相应的水体、阴影和云掩膜。在偏移估计过程中，通过干扰因素的掩膜找出互相关窗口中干扰因素的位置，然后在计算时排除互相关窗口内的干扰因素像素，以此来提高偏移估计的准确性和可靠性。
1 研究区域概述及数据源白格滑坡位于西藏自治区昌都市江达县白格村金沙江右岸，其经纬度为31°4′56.41″N和98°42′17.98″E。该地区位于青藏高原和四川盆地的过渡地带，地势起伏较大，坡度介于35°~55°之间，海拔范围为
|
表 1 卫星影像参数 Table 1 The parameters of radar satellite data |
本文先对包含白格滑坡区域的整景影像进行裁剪，以得到一个大小为3 km×3 km的区域。然后，利用全色波段合成方法合成模拟波段，并将其作为偏移估计的使用波段。接下来，利用陆地卫星掩码(Function of Mask, Fmask)算法找出影像区域内的水体、阴影和云层等干扰因素的位置，并生成对应的掩膜。然后通过水体、阴影和云层的掩膜找出互相关窗口中水体、阴影和云层的位置，并在计算时排除互相关窗口内的水体、阴影和云层像素。最后将自适应偏移估计的结果进行误差处理，获得最终的白格滑坡东西向和南北向地表形变结果图。为了分析白格滑坡的时间演变，本实验对其进行了形变时序监测，使用的Sentinel-2影像对参数如图1所示。
|
图 1 Sentinel-2影像对的时间分布 Figure 1 Distribution of Sentinel-2 image pairs |
由于传感器硬件条件的约束，传感器在某个特定波长范围内的响应无法达到完全响应。因此本文采用波段合成的方法[17-18]，线性合成模拟波段。波段合成计算式为
| $ P=\left(\frac{{\alpha }_{2}}{{\alpha }_{{p}}} \times {B}_{2}\right) +\left(\frac{{\alpha }_{3}}{{\alpha }_{{p}}} \times {B}_{3}\right) +\left(\frac{{\alpha }_{4}}{{\alpha }_{{p}}} \times {B}_{4}\right)$ | (1) |
式中：
Fmask算法[19-20]是一种用于Landsat-8影像的云检测方法，其云检测总体精度达到96.4%。由于Sentinel-2影像的大部分波段与Landsat-8影像相同，所以Fmask也可用于Sentinel-2影像的云层、阴影和水体检测。通过Fmask算法找出影像中水体、阴影和云层的位置，然后对其进行掩膜处理，得到相应的水体、阴影和云层掩膜。
2.4 自适应窗口设计图2为规则窗口的像素分布示意图。绿色方格表示未被干扰因素覆盖的区域，红色方格为匹配窗口中心点(红色边框代表匹配窗口)，黄色和蓝色方格区域表示被不同性质的干扰因素覆盖。本文所提的自适应偏移跟踪算法(Adaptive Offset Tracking，Ad-OT) 在计算过程中只选取那些未被水体、阴影、云层等干扰因素覆盖的像素，这样可以减少水体、阴影、云层对偏移估计的影响，从而提升偏移估计的可靠性。自适应偏移跟踪算法相关系数通过式(2)进行计算。
|
图 2 规则窗口的像素分布示意图 Figure 2 Pixel distribution diagram of regular window |
| $ \rho =\frac{\bigg|{\displaystyle\sum\limits _{(i,j) \in \boldsymbol{W}}}({m}_{i,j}-{\mu }_{\mathrm{m}}) ({s}_{i,j}-{\mu }_{\mathrm{s}}) \bigg|}{\sqrt{{\displaystyle\sum\limits _{(i,j) \in \boldsymbol{W}}}{({m}_{i,j}-{\mu }_{\mathrm{m}}) }^{2}}\sqrt{{\displaystyle\sum\limits _{(i,j) \in \boldsymbol{W}}}{({s}_{i,j}-{\mu }_{\mathrm{s}}) }^{2}}} $ | (2) |
式中：
偏移估计结果包含了轨道误差、条带误差以及失相关噪声，为了得到更准确的结果，需要对误差进行处理。使用一阶多项式曲面模型消除线轨道误差[21]，使用均值相减法消除条带误差[22]，通过偏移量、相关系数、信噪比阈值消除失相关噪声。
3 结果比较白格滑坡主要发生在东西方向上，而在南北方向上的位移较小，因此本文只关注滑坡在东西方向上的位移。由于没有现场测量数据或GPS数据可用，所以使用PlanetScope影像的Trad-OT方法测量结果来验证本文所提的方法。为了方便后续的描述，本文定义了一些缩写。例如，Trad-OT(S2) 表示使用Sentinel-2影像数据的传统偏移追踪计算测量结果。这些缩写将按照以下格式进行表示：偏移估计方法缩写(影像缩写)。文中使用的偏移估计方法缩写主要包括Trad-OT和Ad-OT。而影像缩写则包括S2、L8和PlanetScope。
3.1 滑坡区域的水体干扰分析为了评估Ad-OT抗水体干扰的能力，本文对同一对Sentinel-2影像采用了Trad-OT和Ad-OT两种方法进行滑坡形变监测。图3(a)~(b)分别为两种方法的测量结果，从图中可以明显看出，Trad-OT方法在河流两侧区域出现了错误估计，甚至很多地方的偏移估计都出现了粗差。相反，Ad-OT方法在河流两侧大部分地区的偏移估计值都趋于正常，显示出较好的抗水体干扰能力。此外，图3(c)~(d)都是PlanetScope影像的Trad-OT方法测量结果，区别在于互相关窗口的大小，前者的互相关窗口大小为32×32，后者为102×102。从图中发现河流的干扰范围与互相关窗口大小有关，若窗口尺寸变小，河流的干扰范围则会往河流两侧收缩。
|
图 3 Trad-OT与Ad-OT方法抗水体干扰对比分析 Figure 3 Comparative analysis of water interference resistance between Trad-OT and Ad-OT methods |
为了更直观地对比Trad-OT和Ad-OT两种方法的结果，本实验从主要的滑坡区域和河流两侧区域选取了4条剖面线进行分析，剖面线的位置分布如图4(e)中P1-P2、P3-P4、P5-P6、P7-P8所示。其中，P5-P6、P7-P8处于河流两侧，该区域的测量结果会受到河流的干扰。因此，选择使用互相关窗口大小为32×32的Trad-OT(PlanetScope) 来验证Ad-OT方法在P5-P6、P7-P8区域的测量结果。这种选择的原因在于，互相关窗口大小为32×32时，河流两侧的干扰范围会向内收缩，从而使得Trad-OT方法在P5-P6和P7-P8区域的测量结果不受河流的干扰。
|
图 4 P1-P2、P3-P4、P5-P6、P7-P8位置分布及分析 Figure 4 Position distribution and analysis of P1-P2, P3-P4, P5-P6, P7-P8 |
图4(c)~(d)为P5-P6、P7-P8的剖面线分析，从中可以看出，Trad-OT(S2) 与Trad-OT(PlanetScope) 存在较大差异，而Ad-OT(S2) 整体上与Trad-OT(PlanetScope)接近，局部存在一些差异。这是由于互相关窗口过小导致抗干扰能力减弱，从而使Trad-OT(PlanetScope)结果出现噪声。为进一步量化对比，将Trad-OT(PlanetScope)与Trad-OT(S2) 和Ad-OT(S2) 相差在3 m以内的测量结果定义为正确测量结果。其中，在剖面P5-P6和P7-P8上，Trad-OT(S2) 的正确测量结果占比分别为47%和21%，而Ad-OT(S2) 的正确测量结果占比分别为87%和84%。相较于Trad-OT(S2) ，Ad-OT(S2) 在P5-P6和P7-P8的正确测量结果分别提升了40%和63%。
P1-P2、P3-P4处于主要的滑坡区域，该区域的测量结果不会受到河流的干扰。因此，Trad-OT和Ad-OT两种方法在P5-P6、P7-P8区域的测量结果几乎一样。为了验证Ad-OT方法在P5-P6、P7-P8区域的测量结果，选择使用互相关窗口大小为102×102的Trad-OT(PlanetScope) 进行验证。此外，选择互相关窗口大小为102×102的目的在于确保窗口包含的区域面积与Ad-OT(S2) 窗口保持一致。图4(a)~(b)为P1-P2、P3-P4的剖面线分析，从中可以发现Ad-OT(S2) 与Trad-OT(PlanetScope) 的结果很相似，沿着剖面线的偏移趋势相对一致。
3.2 滑坡区域的阴影干扰分析图5(d)为白格滑坡区域的Sentinel-2光学影像，部分滑坡区域因为云层遮挡太阳光线而出现阴影，阴影区域的位置由图中绿色虚线方框标出。为了评估Ad-OT方法抗阴影干扰的能力，用Trad-OT和Ad-OT两种方法对白格滑坡区域进行形变监测，其结果如图5(a)~(b)所示。从中可以看出，两种方法的测量结果在阴影区域附近存在较大的差异。此外，为了验证Ad-OT在阴影区域附近的测量结果，本实验将互相关窗口大小为102×102的Trad-OT(PlanetScope) 作为验证数据，其结果如图5(c)所示。
|
图 5 Trad-OT与Ad-OT方法抗阴影干扰对比分析 Figure 5 Comparative analysis of shadow interference resistance between Trad-OT and Ad-OT methods |
为了更直观地对比Trad-OT和Ad-OT两种方法的结果，本实验从滑坡区域和阴影区域附近选取了3条剖面线进行分析，剖面线的位置分布如图6(d)中P9-P10、P11-P12、P13-P14所示。图6(a)~(b)为P9-P10、P11-P12的剖面线分析，从中可以看出Ad-OT(S2) 与Trad-OT(S2) 相同，且与Trad-OT(PlanetScope) 大体相近。主要原因是P9-P10、P11-P12远离阴影区域，Trad-OT方法在这个区域的测量结果不会受到阴影的干扰。
|
图 6 P9-P10、P11-P12、P13-P14位置分布及分析 Figure 6 Position distribution and analysis of P9-P10, P11-P12, P13-P14 |
图6(c)为P13-P14的剖面线分析，从中可以看出Trad-OT(S2) 与Trad-OT(PlanetScope) 差异较大，主要原因是P13-P14靠近阴影区域，导致Trad-OT方法在该区域的测量结果受到阴影的干扰。相比之下，Ad-OT(S2) 与Trad-OT(PlanetScope) 更为接近。为进一步量化对比，将Trad-OT(PlanetScope) 与Trad-OT(S2) 和Ad-OT(S2) 相差在3 m以内的测量结果定义为正确测量结果。在剖面P13-P14中，Trad-OT(S2) 的正确测量结果占比为28%，而Ad-OT(S2) 的正确测量结果占比为80%。相较于Trad-OT(S2)，Ad-OT(S2) 在P13-P14的正确测量结果提升了52%。
3.3 滑坡区域的云层干扰分析图7(d)为包含白格滑坡区域的Landsat-8光学影像，其中部分白格滑坡区域被云层所覆盖，这些云层区域的位置已经用绿色虚线方框标出。为了评估Ad-OT方法抗云层干扰的能力，用Trad-OT和Ad-OT两种方法对白格滑坡区域进行形变监测，其结果如图7(a)~(b)所示。从中可以看出，两种方法的测量结果在云层区域附近存在较大的差异。此外，为了验证Ad-OT在云层区域附近的测量结果，本实验将互相关窗口大小为160×160的Trad-OT(PlanetScope) 作为验证数据，其结果如图7(c)所示。
|
图 7 Trad-OT与Ad-OT方法抗云干扰对比分析 Figure 7 Comparative analysis of cloud interference resistance between Trad-OT and Ad-OT methods |
为了更直观地对比Trad-OT和Ad-OT两种方法的结果，本实验从云层区域附近选取了2条剖面线(P15-P16、P17-P18) 进行分析，剖面线的位置分布如图8(c)中黑线所示。图8(a)为P15-P16的剖面线分析，从中可以看出在P15-P16的中间部分，Trad-OT(S2) 与Trad-OT(PlanetScope) 相差较大。主要原因是P15-P16的中间部分靠近云层区域，导致Trad-OT方法在该区域的测量结果受到云层的干扰。相比之下，Ad-OT(S2) 与Trad-OT(PlanetScope) 更为接近。为进一步量化对比，将Trad-OT(PlanetScope) 与Trad-OT(S2) 和Ad-OT(S2) 相差在3 m以内的测量结果定义为正确测量结果。在剖面线P15-P16中，Trad-OT(S2) 的正确测量结果占比为77%，而Ad-OT(S2) 的正确测量结果占比为98%。相较于Trad-OT(S2) ，Ad-OT(S2) 在P13-P14的正确测量结果提升了21%。
|
图 8 P15-P16、P17-P18位置分布及分析 Figure 8 Position distribution and analysis of P15-P16, P17-P18 |
图8(b)为P17-P18的剖面线分析，从中可以发现离P17的越近部分，Trad-OT(S2) 与Trad-OT(PlanetScope) 的差异越大。这是因为剖面线P17-P18的P17端点靠近云层，导致越靠近P17的测量结果受到的干扰越大。相比之下，Ad-OT(S2) 与Trad-OT(PlanetScope) 更为接近。
4 滑坡形变的时间序列分析本实验选取了2016年5月11日至2018年6月10日期间的7景Sentinel-2卫星影像进行了自适应偏移估计，得到18个影像对的偏移估计结果。接着对这些结果进行时间序列形变计算，得到了最终的白格滑坡时序形变结果，具体结果如图9所示。其中，测得的最大累计形变量为30.18 m。
|
图 9 滑坡时间序列形变结果 Figure 9 Time series deformation results of slide area |
为了更好地说明白格滑坡的地表形变趋势，将白格滑坡的时序形变过程分为稳定滑移阶段(形变速率小于23 mm/d)和加速滑移阶段(形变速率大于23 mm/d)。随后，从坡顶到坡底选取了6个特征点(如图9(a)中A～F所示)进行时序分析，其分析结果如图10所示。其中，A点位于滑移区域的顶部，从2017年5月16日才开始进入加速滑移阶段，形变速率从原来的15.68 mm/d加速至25.72 mm/d。B、C两点位于滑坡的中部，从2017年5月16日才开始进入加速滑移阶段，形变速率分别从原来的22.74、20.72 mm/d加速至54.67、51.72 mm/d。D点位于滑移的中下部，从2017年10月18日才开始进入加速滑移阶段，形变速率从原来的15.34 mm/d加速至57.50 mm/d。E点位于滑移区域的底部，从2018年3月2日才开始进入加速滑移阶段，形变速率从原来的10.58 mm/d加速至55.06 mm/d。F点靠近河流，在此期间几乎没有发生形变。从分析结果中可以发现，滑坡的顶部和中部最早进入加速滑移阶段，并且累计形变量最大。随后，滑坡中下部进入加速滑移阶段，最后是滑坡底部。
|
图 10 6个特征点的累计形变 Figure 10 Time-series deformation of six selected points |
本文提出了一种自适应偏移估计方法。该方法通过干扰因素的掩膜找出互相关窗口中的干扰因素像素，并在计算时排除这些像素，从而提高了偏移估计的准确性和可靠性。将该方法应用于白格滑坡，实验结果表明，该方法确实能够减少水体、阴影、云层对周围区域偏移估计的影响。其中，在剖面线P5-P6、P7-P8、P13-P14、P15-P16上，该方法的正确测量结果占比分别为87%、84%、80%、98%。相较于传统偏移估计方法，提升的百分比分别为40%、63%、52%、21%。
通过对白格滑坡形变的时间序列分析，测得在2016年5月11日至2018年6月10日期间该滑坡的最大累计位移形变为30.18 m。滑坡的时序形变过程可分为稳定滑移阶段和加速滑移阶段，在2016年5月11日至2017年5月16日期间，滑坡的主体区域处于稳定滑移阶段。从2017年5月16日开始，滑坡的顶部和中部进入加速滑移阶段。之后滑坡中下部和底部也相继进入加速滑移阶段。
| [1] |
李梦华. 基于星载雷达遥感的大梯度滑坡形变监测与失稳致灾时间预测[D]. 武汉: 武汉大学, 2019.
|
| [2] |
陈炳杰. 基于光学和雷达影像像素偏移跟踪技术的形变监测研究[D]. 广州: 广东工业大学, 2021.
|
| [3] |
STUMPF A, MALET J P, PUISSANT A, et al. Monitoring of earth surface motion and geomorphologic processes by optical image correlation[J]. Land Surface Remote Sensing, 2016: 147-190.
|
| [4] |
ZHAO C, LU Z. Remote sensing of landslides—a review[J].
Remote Sensing, 2018, 10(2): 279.
DOI: 10.3390/rs10020279. |
| [5] |
INTRIERI E, RASPINI F, FUMAGALLI A, et al. The maoxian landslide as seen from space: detecting precursors of failure with Sentinel-1 data[J].
Landslides, 2018, 15(1): 123-133.
DOI: 10.1007/s10346-017-0915-7. |
| [6] |
ZEBKER H A, VILLASENOR J. Decorrelation in interferometric radar echoes[J].
IEEE Transactions on Geoscience and Remote Sensing, 1992, 30(5): 950-959.
DOI: 10.1109/36.175330. |
| [7] |
陈佳炜, 吴希文, 王华, 等. 基于双极化时序InSAR技术的地表形变监测[J].
广东工业大学学报, 2022, 39(3): 77-82.
CHEN J W, NG A H M, WANG H, et al. Surface deformation monitoring based on dual-polarization time series InSAR technology[J]. Journal of Guangdong University of Technology, 2022, 39(3): 77-82. DOI: 10.12052/gdutxb.210084. |
| [8] |
STURKELL E, EINARSSON P, SIGMUNDSSON F, et al. Volcano geodesy and magma dynamics in iceland[J].
Journal of Volcanology and Geothermal Research, 2006, 150(1/3): 14-34.
|
| [9] |
MICHEL R, AVOUAC J P, TABOURY J. Measuring near field coseismic displacements from SAR images: application to the Landers earthquake[J].
Geophysical Research Letters, 1999, 26(19): 3017-3020.
DOI: 10.1029/1999GL900524. |
| [10] |
PATHIER E, FIELDING E, WRINGT T, et al. Displacement field and slip distribution of the 2005 Kashmir earthquake from SAR imagery[J].
Geophysical Research Letters, 2006, 33(20): L20310.
|
| [11] |
YOO J C, HAN T H. Fast normalized cross-correlation[J].
Circuits, Systems and Signal Processing, 2009, 28(6): 819-843.
DOI: 10.1007/s00034-009-9130-7. |
| [12] |
HU X, WANG T, LIAO M. Measuring coseismic displacements with point-like targets offset tracking[J].
IEEE Geoscience and Remote Sensing Letters, 2014, 11(1): 283-287.
DOI: 10.1109/LGRS.2013.2256104. |
| [13] |
BAMLER R, EINEDER M. Accuracy of differential shift estimation by correlation and split-bandwidth interferometry for wideband and delta-k SAR systems[J].
IEEE Geoscience and Remote Sensing Letters, 2005, 2(2): 151-155.
DOI: 10.1109/LGRS.2004.843203. |
| [14] |
STUMPF A, MALET J P, DELACOURT C. Correlation of satellite image time-series for the detection and monitoring of slow-moving landslides[J].
Remote Sensing of Environment, 2017, 189: 40-55.
DOI: 10.1016/j.rse.2016.11.007. |
| [15] |
BONTEMPS N, LACROIX P, DOIN M P. Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in Peru[J].
Remote Sensing of Environment, 2018, 210: 144-158.
DOI: 10.1016/j.rse.2018.02.023. |
| [16] |
CAI J, WANG C, MAO X, et al. An adaptive offset tracking method with SAR images for landslide displacement monitoring[J].
Remote Sensing, 2017, 9(8): 830.
DOI: 10.3390/rs9080830. |
| [17] |
OTAZU X, GONZALEZ-AUDICANA M, FORS O, et al. Introduction of sensor spectral response into image fusion methods. application to wavelet-based methods[J].
IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(10): 2376-2385.
DOI: 10.1109/TGRS.2005.856106. |
| [18] |
BOGGIONE G A. PIRES E G, SANTOS P A, et al. Simulation of a panchromatic band by spectral combination of multispectral ETM+bands[C]//PUTKONEN J. Proceedings of International Symposium on Remote Sensing of Environmental. Honolulu: IEEE, 2003: 321-324.
|
| [19] |
ZHU Z, WANG S, WOODCOCK C E. Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4-7, 8, and Sentinel 2 images[J].
Remote Sensing of Environment, 2015, 159: 269-277.
DOI: 10.1016/j.rse.2014.12.014. |
| [20] |
ZHU Z, WOODCOCK C E. Object-based cloud and cloud shadow detection in Landsat imagery[J].
Remote Sensing of Environment, 2012, 118(6): 83-94.
|
| [21] |
FENG G, DING X, Li Z, et al. Calibration of an InSAR-derived coseimic deformation map associated with the 2011 Mw-9.0 Tohoku-Oki earthquake[J].
IEEE Geoscience & Remote Sensing Letters, 2012, 9(2): 302-306.
|
| [22] |
MICHEL R, AVOUAC J P. Coseismic surface deformation from air photos: the Kickapoo step over in the 1992 Landers rupture[J].
Journal of Geophysical Research-Solid Earth, 2006, 111(B3): 240-247.
|