使用DInSAR、SBASInSAR、GNSS 和地下测量的方法结合高分辨率LiDAR DEM进行滑坡监测

概述

中国台湾地区70%为山坡。此外，地势波动剧烈，地震和板块运动活跃。山体滑坡是台湾地区普遍的灾害。然而，山体滑坡并不是灾难，除非有人进入山区进行开发。滑坡位移监测是本研究的首要任务。本研究选择了以板岩地质条件为主的潜在滑坡区作为候选地点。该区板岩层理向东南方向倾斜约30～75度，可能因高角度岩层撞击引起的重力变形而发生蠕变。此外，由于研究地点位于植被茂密的地区，数据噪声非常高，不容易获得好的结果。本研究选择Sentinel-1数据、1米LiDAR DEM作为参考高程。高精度的1米LiDAR DEM 可以帮助检测来自 DInSAR 的更复杂的变形。此外，研究区域的道路、桥梁和建筑物的农场区域提供了足够的散射体以保证相干性。分析了2017 年 3 月至 2021 年 6 月的 Sentinel-1 数据，获得了边坡变形并将其转换为垂直方向。将 SAR 得出的形变与其他测量值进行比较，包括 GNSS 和地下边坡测斜仪。SBAS 求解过程提供了更多的差分干涉像对来克服巨大噪声的问题并提高了精度。而且，SBAS方法做了参数修改在植被区域中得到了更多的候选点。GNSS安装位置与升轨SBAS得到的垂直变形对比结果是一致的。此外，还讨论了斜坡朝向SAR卫星的可靠朝向。由于GNSS台站的局限性，本研究提出了一种将边坡测斜仪观测到的形变转换为垂直形变的方法。边坡测斜仪的位移原本是水平位移。假设它固定在地下稳定的地方，自下而上的移动与深度有关。结果表明，本文所提出的将水平位移转换为垂直位移的方程非常适合这种情况。

一、简介

台湾地处两个板块的交界处，地震多发。板块碰撞也导致造山运动，增加了山脉的高度。由于该岛被不同的海水温度所包围，极端天气变化非常快，并引发很多台风和降雨。例如，2009年一次台风事件的最高累积降雨量为3000毫米。降雨量约为年平均降雨量的四分之三。这些降雨给台湾地区带来很多潜在的滑坡。山体滑坡减灾是一个重要课题。地质调查部门在 2009 年台风莫拉克之后的十年间，开始利用 LiDAR DEM 绘制大规模滑坡（深部滑坡）区域。https://data.gov.tw/en（2021年 6 月 1 日访问）现有的区域提供了这些滑坡事件发生的年份，但不能代表当前的状况。本研究下载并详细绘制的滑坡区域如图1所示。区域按面积应大于10公顷的标准绘制。因此，本研究进行了详细的调查和绘图。在实地调查和 LiDAR DEM 测绘过程中发现了几处轻微的滑坡和裂缝。研究区地质构造为板岩，随砂岩含量的不同而变化。板岩地层朝向东南方向，角度为20-70度，可能会因重力而导致缓慢滑动。这些滑坡的变形速度是安装监测系统或工程处理的首要关键问题。

图 1. 研究区 LiDAR 绘制的滑坡痕迹及其在台湾的位置

一般来说，PSInSAR 是一种推导大面积变形的方法，例如地面沉降或断层变形，很多学者使用InSAR方法结合其他地学方法对滑坡活动进行监测，所使用的InSAR方法有PSInSAR、SBASInSAR、DInSAR、A-DInSAR、TCP-InSAR、SqueeSAR等 [参考文献1-28 ]。

台湾国立暨南大学王国龙教授团队使用SBAS方法进行研究。在确认地表运动后进行不同测量值的比较。变形值和速度证明了活动性，可以提供有关滑坡机制的信息。但 DInSAR 应注意滑坡位移大面积分析的正确性。鉴于此，将卫星相对变形（LOS）转换为垂直方向是唯一可行的验证方法。使用的方法是SARScape软件提供的程序。

二、SAR数据处理方法

DInSAR、PS InSAR 和 SBAS InSAR 近年来发展非常迅速。这些技术对于大规模滑坡监测变得更加有效。已有学者提出了各种方法来实现更快的处理和更好的结果 [ 参考文献29-38 ]。

本研究使用Sentinel-1同极化数据，1 m LiDAR DEM作为参考DEM，使用SARscape进行DInSAR处理，处理流程如下图。

图 2 DInSAR处理流程图

DInSAR方法适用于观测基线短、周期短的两个场景的滑坡伤痕。然而，从两个图像中获得的相干性、幅度和条纹很难连续监测位移。因此，在研究区引入了时序InSAR方法。该区域采用 PSInSAR 和 SBASInSAR，但PSInSAR结果的候选点少得多。从PSInSAR中提取的候选点数仅为 SBAS 的 1/10。因此，本研究仅使用了SBASInSAR（以下简称SBAS）方法。

SBAS通过观测时间并跟踪他们在每次观测中的位置来搜索具有相似雷达信号幅度的候选点。处理过程中的大气效应使用时间天数和过滤器来消除。选择2017年3月- 2021年6月的Sentinel-1数据，升轨126时相，降轨121时相，设置时间基线阈值60天，空间基线阈值150 m，生成的SBAS连接图如下图所示：

图3-1 升轨数据集连接图

图3-2 降轨数据集连接图

直接来自 SARScape软件处理的RMS垂直位移误差显示为 8.6 mm，升轨数据集的置信度为 95%，降轨数据集的RMS垂直位移误差显示为9.8 mm，置信度为 95%。升轨数据集的入射角方位角为80.1度，倾角为39.1度。降轨数据集的入射角方位角为281.5度，倾角为36.9度。研究区朝向以东向为主。坡向（方位角）和坡度（陡度）会影响候选点的数量及其准确性。升轨数据集得到的候选点较多，如图4(1). 这是因为研究区的斜坡正对着升轨的 SAR 卫星信号。此外，SARScape 的一些参数针对研究区域进行了优化，以产生最多的候选点。接下来就是验证该研究区域中每个轨道结果的准确性。

图4-1 升轨数据集SBAS处理得到的SBAS候选点

图4-2 降轨数据集SBAS处理得到的SBAS候选点

三、结果

3.1全球导航卫星系统（GNSS）测量

安装了两个双频GNSS站进行分析，以验证SBAS分析结果的准确性，如图1所示，自2017年4月开始。采用位于距离该点1.1公里的稳定位置的GNSS参考站作为节点，进行节点位移计算。假设参考 GNSS 站稳定无位移。GNSS 站址的选择基于初步 DInSAR 分析。差分 GNSS 计算基于开源代码 RTKLIB。由于大雨，道路可能会被阻塞。因此，记录的数据通过互联网发送到实验室，并每小时在另一台服务器上进行计算。Sentinel-1 数据的重访周期为 12 天。每小时的GNSS采样对于 DInSAR 验证来说太密集了。因此，采用了 24 小时的采样频率。

如图5和图6所示，从 2017 年 4 月到 2021 年 6 月有两次强降雨。第一次是2017年6月，即安装 GNSS 站仅两个月后， 6月1日至6月4日降雨量为902.5毫米，6月11 日至6月18日降雨量为640毫米。该事件在GNSS1处向东方向产生了约120毫米，向南方向产生了38毫米，在垂直方向产生了78毫米的位移。第二次强降雨是2019年，每小时降雨量达到了52毫米的峰值，但累积降雨量比2017年那次少。在这次事件中测量了位移。GNSS2 向南的位移大于向东，因为它的坡向是向南的；GNSS2的水平位移与垂直位移几乎相同；但GNSS1的水平位移远大于垂直位移。位移结果意味着GNSS2的位置在单个滑动面，GNSS1在多个滑动面，也就是说下面有更深的滑动面。现场勘查和数值模拟以验证这一观察结果。

图 5. GNSS1 处测量的 GNSS 位移和降雨量。

图 6. GNSS2 处测量的 GNSS 位移和降雨量。

采用垂直位移进行验证。提取GNSS1和GNSS2处的垂直位移与SBAS的结果比较。如图7和图8所示，GNSS1 垂直位移在cm精度内与升轨数据集的SBAS结果非常一致。由于地形的原因，降轨数据集的SBAS结果与实际的位移有很大不同。这一结果表明：在分析滑坡活动时，山体本身对SAR数据轨道的选择起到控制性作用。GNSS2安装在庐山小学内，斜坡朝南。对比结果表明，降轨数据集的SBAS结果比升轨数据集的SBAS结果更加一致。但是，升轨数据集经过四年的计算，得到了更加准确的形变结果。因此，以下结果讨论将集中在本研究区的升轨数据集上。此外，有个现象是：在2017年降雨事件之后，GNSS和SBAS都测量到了隆起现象。该机制将在下一节中讨论。

图 7. 在 GNSS1 处测量的 GNSS 位移和SBAS测量得到的位移。

图 8. GNSS2 处测量的 GNSS 位移和 SBAS测量得到的位移。

GNSS台站的坡度（陡度）和坡向（方位角）如表1所示。对于DInSAR结果的准确性，采用 GNSS 测量作为精确参考。使用DInSAR结果和GNSS当天的数据。如表1所示在 GNSS1 数据中，升轨结果显示 13 mm 的误差和 10 mm 的标准偏差。此外，相关系数为0.95，表示高度相关。对于GNSS1与降轨结果对比，误差为 84mm。相关系数为-0.69，这意味着该位置不能使用降轨数据。对于GNSS2站点，升轨数据的误差和标准差比GNSS1更大，相关系数仅为0.51。结果表明，DInSAR 结果在该位置仍然适用，但精度较低。降轨数据也显示出较低的准确性，不适合该位置。

注：相关系数是 GNSS 和 DInSAR 测量值之间的关系。如果系数接近 1，则表明这两个数据相同。另一方面，如果相关系数小于零，这意味着完全不同的结果，甚至表示较小的平均误差。

表 1 GNSS 台站斜率（陡度）和方位角（方位角）及相关 DInSAR 比较误差

	GNSS1 Location
	坡度 (degree)	坡向 (degree)
	15.6	88.9
	升轨	降轨
误差均值 (mm)	13.21	84.74
误差标准差 (mm)	10.14	37.82
相关系数	0.95	−0.69
	GNSS2 Location
	坡度 (degree)	坡向 (degree)
	19.8	169.8
	升轨	降轨
误差均值 (mm)	24.43	25.76
误差标准差 (mm)	15.17	17.52
相关系数	0.51	−0.20

 

通过以上分析，可以得出结论：在该研究区，SAR升轨数据的结果优于SAR降轨数据的结果。结果表明，升轨的SAR数据适用于东向斜坡。

3.2. 与边坡测斜仪的滑动测量相比较

研究人员十多年来对该研究区一直在进行观测。在此之前，连续的 GNSS 测量非常昂贵。钻孔是了解地下地质和变形测量的传统方法。本研究采用BH1和BH2两个钻孔，如图1所示。BH1 和 BH2 的深度分别为100 m和80 m。形变测量从钻孔的最底部开始。该算法假设最底部没有移动，并测量两个垂直方向的倾斜角。然后可以通过倾角乘以深度得出每个深度的水平位移。与SBAS方法相比，GNSS受限时，水平位移可作为参考。水平位移如图9所示，随着每日降雨。测量的水平方向有两个，分别是A方向和B方向。这两个方向相互垂直。水平位移的方向很难确定，因为地下管道可能会扭曲。因此，我们将两个方向合并为一个位移。结果表明，在这种情况下，A方向是受控方向。边坡测斜仪从 2011 年开始进行测量。2012 年和 2017 年分别发生了两次强降雨事件。两次都是超过1000毫米的强降雨。由于 SBAS 只能检测地表形变，因此本研究以地表随时间的位移作为参考。

图 9. 坡度测斜仪测得的地表水平位移。

本研究提出了将水平位移转换为垂直位移的概念，示意图如图 10。边坡测斜仪测量参考地表深处的水平位移。钻孔处的平均坡度，同一点的垂直位移可根据下列公式估算。计算坡度的推荐距离（半径）是从孔的中心到大于测斜仪深度的半径。

tan θ =→V= H× tan θ

其中, θ : 边坡测斜仪钻孔处的平均坡度；H : 边坡测斜仪观测到的地表位移；V : 估计的垂直位移。

图 10. 坡度测斜仪测量横向地面位移和垂直位移转换的图示

转换后的垂直位移如图 11 所示，将SBAS处理得到的垂直位移在同一时期绘制以进行比较。结果表明，该方法可以正确估计BH1处的垂向位移趋势。但 BH2 处的位移相比只是变形趋势一致，值并不准确。所提出的方法是基于滑坡监测，其中观测点总体上是向下移动的。如果观察点向上移动（隆起），则应仔细检查该过程。这种方法允许比较 DInSAR 和边坡测斜仪数据，尤其是没有 GNSS 或水准数据的传统监测站点。比较这些位置时，有两个有趣的观察结果：一个是在 2017 年强降雨之后，形变趋势出现了反方向的变化；二是这些地点没有受到2019年强降雨的影响。

图 11. 坡度指标与升轨SBAS形变结果相比的转换垂直位移。

3.3. SBAS方法进行潜在滑坡制图

由于 SBAS 方法的准确性在该地区已通过GNSS和边坡测斜仪得到验证，因此在研究区基于该方法进行滑坡制图是可行的。为了建立一个标准的操作程序来绘制滑坡区域，在几次测试后执行以下处理步骤：

（1）SBAS处理之后从 LOS 位移计算垂直位移；

（2）将垂直位移内插为栅格数据；

（3）将栅格垂直位移与从DEM生成的山体阴影图叠加；

（4）将得到的位移与Lidar识别的滑坡区进行比较。

处理后绘制出潜在的滑坡区，如图12所示，同时还叠加了LiDAR识别的滑坡区。图中显示大部分 LiDAR 基础滑坡区都在移动，这说明用LiDAR数据绘制的滑坡是正确的。但是，滑动位移也意味着优先做进一步的工程监测。SBAS 方法可作为选择减灾地点的可能的技术手段。基于SBAS方法绘制的潜在滑坡区比LiDAR绘制的滑坡区域多。这些潜在的滑坡区正在移动，但由于连续的规则滑动，滑动阈值难以定义。这种基于LiDAR DEM的SBAS形变结果表明，大部分沉降发生在滑坡的顶部，此外，在斜坡的低海拔处可以观测到一些抬升。

图 12. 2017年3月到2021年6月Sentinel-1数据SBAS处理得到的垂直位移速率

3.4. 抬升现象的数值模拟与现场调查

在上图5中的垂直位移显示了2017年强降雨事件后的抬升行为。GNSS1的观测表明水平位移远大于垂直位移。检测到的小的形变意味着存在多个滑动平面。对该站点进行现场调查以找到滑动机制。现场调查发现，墙上的裂缝，如图13，在屋顶顶部安装了 GNSS。挡土墙到底部的位移大于顶部，说明了圆形破坏情况，这可能导致滑动底部具有较高高度的旋转。为了模拟这种情况，采用了一个简单的斜率模型：Rocscience 的 RS2 进行有限元分析。坡顶至坡脚剖面进行切割模拟滑坡位移。模拟结果如图14所示，图中标出了GNSS 位置。模拟了位移的水平分量和垂直分量。结果表明，在GNSS位置，水平位移比垂直位移更显著。此外，模拟结果说明该站点至少有两个滑动平面。从数值模拟中识别出两个滑动。主要滑动发生在降雨量更大的情况下，这导致整体滑动和 GNSS 垂直高程下降。小滑坡发生在大滑坡稍停时，小雨袭来，造成圆周旋转，垂直高程增加。

图 13. 挡土墙裂缝显示底部位移大于顶部（GNSS1 安装在建筑物屋顶）

图 14 有限元数值位移模拟结果

四、结论

本研究提出了基于DInSAR和SBAS的滑坡区处理和验证方法。来自LiDAR DEM的滑坡制图显示了过去的滑坡事件，对当前的情况应识别活动滑坡，如果变化速度较快，应进行进一步监测和处理。

DInSAR可以检测滑坡区域并进行位移图绘制。一旦在两个场景中发生了事件，DInSAR 可以使用适当的阈值从边缘和位移图绘制剖面区域。然而，DInSAR 在像对的选择上、不可避免会产生噪声信号，而且很难确定斜坡是否正在滑动等，这些方面都存在问题。DInSAR 的结果在接下来的分析中并不那么精确且易于使用。

GNSS 数据可以检测三维位移并有助于解释滑坡行为。如斜坡中可能的滑动面和滑动位置。SBAS 在假定模型的情况下具有更精确的位移推导，因此可以定义滑坡边界。本研究中使用的 SBAS 方法采用了参数优化，从而在植被区域中产生了更多的候选点。选择了几个深部滑坡进行检测，与DInSAR相比，SBAS结果表明：该方法更适合于蠕滑滑坡测绘与监测。GNSS 的验证表明，SBAS可以得到高精度的滑坡监测结果。根据SBAS与GNSS结果的相关性，升轨数据适合研究区。判断是否得到高精度的位移并不是SBAS或PS方法的唯一问题。还应考虑 SBAS和GNSS之间的相关性。

本研究还提出了一种垂直位移估计方法。许多滑坡现场都有足够的坡度测斜仪数据，这是一种典型的地下地表形变监测方法。坡度测斜仪获得水平地表位移，假设孔底保持在同一位置。尽管有 GNSS 数据，估计的垂直位移仍有助于提供更多的垂直变形。然而，斜坡倾角仪每次测量几个月。结果表明滑动行为和启动时间与 SBAS 垂直位移相匹配。坡度测斜仪数据非常珍贵，可能有几十年的数据。这些数据可以得到更多有效的应用。

数值模拟验证了斜坡的滑动行为，引起垂直位移的不同组分。因此，这种方法可以进行大面积的滑坡监测，并且可将监测仪器选择安装在快速滑动的位置。此外，SBAS方法可以检测不稳定滑坡，为工程处理或监测工作提供预警。

本研究提出了基于SBAS方法的潜在滑坡监测，与野外调查相一致。该方法可进行滑坡制图，为优先做好稳定滑坡的工程治理提供技术支撑。

原文链接：https://www.mdpi.com/2076-3417/11/23/11389/htm （内附参考文献）