基于Stacking-InSAR 和SBAS-InSAR 地面沉降监测方法研究
——以新疆阿勒泰地区为例

孙 鹏

（新疆水利水电勘测设计研究院有限责任公司，新疆乌鲁木齐 830000）

0 引言

合成孔径雷达干涉测量技术（Interferometric Synthetic Aperture Radar，InSAR）因其具有高精度、高分辨率、全天候等优点，是监测地表形变中一种方法。该技术应用在城市地表沉降监测、矿区地表沉降监测、山体滑坡监测等项目中，发挥了极大的作用。InSAR 技术演化出Stacking-InSAR 技术与SBAS-InSAR 技术，并可将二者结合起来作为一种组合方法应用在地表沉降、滑坡监测中［1-2］。

本文利用InSAR 技术对新疆阿勒泰地区布尔津滑坡进行监测与分析。Sentinel-1 系统在周边区域拥有丰富的InSAR 数据积累、高重访周期与大面积的影像覆盖范围，收集基于2014 年10 月20 日至2021 年2 月9 日采集的Sentinel-1 升轨影像，数据处理以增强型SBAS-InSAR 集成时序分析为主要技术手段，获取布尔津滑坡区域形变的趋势、范围和梯度等信息，并从中筛选出地表形变变化量达到预警临界值，和发生形变突变的可疑区域，作为地表形变监测预警的重点工作区域。

1 基本原理与技术优势

1.1 Stacking-InSAR技术

干涉图叠加（Stacking-InSAR）技术的数据处理基础是进行经典二轨差分InSAR 技术处理。Stacking-InSAR 技术是在获取M 幅相位解缠的差分干涉图后，将解缠相位加权叠加，估计出平均相位变化速率。由于大气相位分量在空间域具有低频特性、时间域具有高频特性，故Stacking 技术最大程度地减少了大气误差，提高了形变速率解算精度［3-4］。该方法假设每一幅独立的解缠相位图中，在不考虑其他噪声的分布情况下，区域的形变速率为总形变量与时间的比值，即线性形变速率，大气延迟误差随机且相等。通常情况下，为保证结果的准确性，选择时间和空间基线都比较短的数据进行组合。

1.2 SBAS-InSAR技术

时序分析采用小基线集干涉测量技术（Small Baseline Subset，SBAS）充分考虑大量的干涉图信息，能有效地消除和削弱解缠粗差、大气误差以及DEM误差等因素的影响，从而获取更高精度的地表形变成果［5］。基于SBAS 技术的数据处理主要解算步骤如下：

1）对N 幅SAR 影像数据按一定时空基线条件进行干涉组合处理，形成M 幅干涉图，利用已有DEM（本次采用SRTM）作为外部高程数据，进行差分处理，生成差分干涉图，以去除地形及平地效应影响；

2）对差分干涉图进行自适应滤波处理以去除相位噪声影响，对滤波后的差分干涉图进行相位解缠，生成解缠相位图；

3）根据基线条件和干涉相位信息估算高程误差及线性形变相位，在原始干涉相位中减去估算的线性形变相位得到残余相位，此时的残余相位中主要包括非线性形变及大气相位；

4）解缠此残余相位，并补偿线性形变相位部分得到完整的形变相位，此时相位中还包含有大气相位的影响；

5）对去除线性形变的残余相位进行空域低通和时域高通滤波处理以分离出大气相位；

6）在形变相位中减去大气相位影响，得到形变相位值；

7）基于SVD 的形变求解；

8）对形变进行地理编码，获取WGS-84 坐标系下的形变成果图。

1.3 总体技术路线

本文根据监测区域基本情况，计划出合理可行的技术路线。应用雷达监测数据，采用Stacking-InSAR 技术与SBAS-InSAR 技术相结合方式实现了变形数据提取。根据地质变形情况，及时建立灾害预警与防治措施。本文总体技术路线如图1 所示。

图1 Stacking-InSAR 和SBAS-InSAR 时序分析提取地表形变技术路线图

2 研究区域与数据源

2.1 研究区域

布尔津滑坡位于新疆阿勒泰地区布尔津县境内，滑坡位于水电站库区的上游区域。该区域山体受雨水的冲刷导致地质情况不稳定，很容易出现山体滑坡现象，造成河流堵塞带来的发电站隐患问题。研究区域范围内无居住区，地势起伏较大；最低点海拔+203.4 m，最高点海拔+1125.9 m，平均海拔+841.9 m。该滑坡的研究范围如图2 所示。

图2 研究区域范围

为比较本文提及的Stacking-InSAR 和SBASInSAR 技术在滑坡体监测中的有效性，在滑坡区域布设监测点，并采用传统极坐标法获取监测点位移。为便于体现各监测点位移方向，建立监测坐标系，X方向为平行河道，Y方向为垂直河道，通过坐标变化量直观表示各点位移趋势，监测基准网及监测点布置如图3 所示。

图3 监测基准网及监测点布置图

2.2 数据源

2.2.1 影像数据

为避免SAR 影像几何畸变的影响，所用SAR数据为升轨Sentinel-1 数据，采用TOPS（Terrain Observation With Progressive Scans）扫描方式，共计153 景，覆盖时间段为2014 年10 月20 日至2021年2 月9 日。数据基本参数如表1 所示，影像采集范围如图4 所示。

表1 升轨SAR 数据的基本参数

图4 研究范围数据覆盖图

2.2.2 DEM获取

形变监测过程中，数字高程模型（DEM）是非常重要的外部数据，在生成差分干涉图时需要DEM 来模拟地形相位，此外在雷达坐标系和地理坐标系相互转化（地理编码）的过程中也需要高精度的DEM。

美国太空总署（NASA）和国防部国家测绘局（NIMA）于2000 年联合发射的“奋进号”航天飞机完成了为期11 天的航天飞机雷达制图任务（Shuttle Radar Topography Mission，SRTM）。经过两年的数据处理，最终获得了平面精度±20 m，高程精度±16 m 的全球数字高程模型，分辨率包含1 弧秒和3 弧秒两种，分别对应30 m 和90 m。本项目选用30 m分辨率的DEM 数据用于地形相位的消除。

2.2.3 极坐标差分法获取数据

依据153 景的SAR 影像的成像时间，选取2014 年10 月20 日至2021 年2 月9 日之间的SAR数据，每次监测方法保持一致，采用相同观测人员、仪器使用极坐标差分法观测。结合每次观测成果，分析得出各监测点观测成果相对首期成果水平位移量和垂直位移量。位移量的X方向为正表示沿河向上游方向，反之沿河向下游方向，Y方向为正表示垂直河道向河道方向，反之垂直河道向山顶方向；H方向为0 表示无变化，为负表示下沉。为了验证Stacking-InSAR 与SBAS-InSAR 对沉降边缘与沉降中心的监测能力，选取了图2 中的6 个特征点，如表2 所示。其中，1-1、1-3、2-1 和2-2 属于沉降边缘区域，分布于滑坡区域西北侧；1-3 和1-4属于沉降中心区域，2-3 和2-4 位于滑坡区东南侧。

表2 极坐标法获取的变形点坐标值

3 数据处理与结果分析

对获取的153 景Sentinel-1 数据进行雷达差分干涉测量处理，按照时间序列构建152 个干涉对，其中2020.08.15-2020.08.29 的垂直基线最长，达到了200.89 m，2016.04.15-2016.04.20 的垂直基线最短，数值为43.18 m。按照图1 所示的Stacking-InSAR 和SBAS-InSAR 技术处理流程对152 个干涉对处理，获取到相邻影像之间的形变信息，通过时序累加，得到各时间节点上的累计沉降情况。

依据干涉图的相干性最大原则，选取2020 年8月20 日的影像为主影像，其余影像与其配准，分别设定时间基线阈值为100 d，空间基线阈值为220 m。将所有满足阈值的限定的两景SAR 影像两两组合，形成了152 个时序差分干涉对进行处理。按照本文的技术路线图，采用Stacking-InSAR 与SBAS-InSAR 技术处理152 个干涉对，得到研究区域形变信息。布尔津滑坡153 景SAR 数据计算获得的2014 年10 月至2021 年2 月地表形变累计沉降量如图5 所示。布尔津滑坡2014 年10 月至2021年2 月地表形变速率如图6 所示。从这两幅图可以看出：布尔津滑坡区域这段时期内最大累计形变量达-550 mm，最大形变速率达-76.9 mm/a。图6中，负值表示地表正在远离卫星，正值表示地表正在靠近卫星。从图6 可以看出：新疆阿勒泰布尔津滑坡区域分布有三处明显变形体，变形体1 平均形变-6.0 mm/a，最大形变-25.4 mm/a；变形体2 平均形变-4.8 mm/a，最大形变-17.1 mm/a；变形体3 平均形变-19.0 mm/a，最大形变-76.9 mm/a。

图5 地表累计形变图

图6 2014-2021 平均形变速率图

为了更直观感受研究区域内的时空演变规律，这里展示SBAS-InSAR 的时序结果。2014 年10月、12 月与2021 年1 月、2 月时空演变情况如图7-图10 所示。

图7 2014 年10 月时空演变结果

图8 2014 年12 月时空演变结果

图9 2021 年1 月时空演变结果

图10 2021 年2 月时空演变结果

为进一步研究该滑坡形变的时间演化，在3 个形变体上分别提取3 个特征点，获得滑坡上这些点2014 年10 月至2021 年2 月的时间序列形变，如图11-图16 所示。变形体1 最大的累积形变出现在其P3 点上，这个时间段内累积形变达到了-174mm；变形体2 最大的累积形变出现在其P2 点上，这个时间段内，累积形变达到了-223 mm；变形体3 最大的累积形变出现在其P3 点上，这个时间段内，累积形变达到了-499 mm。3 个变形体在2014 年10 月至2021 年2 月存在明显的非线性形变趋势。可以看到，在2018 年1 月至2021 年2月，3 个变形体均出现快速形变现象，形变量比较大，其中变形体3 从2016 年年初就开始出现加速的情况。

图11 变形体1 特征点示意图

图12 变形体1 地表变形速率及特征点累计变形量

图13 变形体2 特征点示意图

图14 变形体2 地表形变速率及特征点累计变形量

图15 变形体3 特征点示意图

图16 变形体3 地表形变速率及特征点累计变形量

4 精度验证

4.1 传统测量仪器结果比对

为了验证Stacking-InSAR 和SBAS-InSAR 技术在滑坡监测方面的精度，采用传统测绘方法在滑坡体上布置一定的测点。考虑到滑坡体不同位置移动情况不同，在沉降边缘区域选择1-1 号点位，该点位于滑坡体西北侧；沉降中心区域选择1-3、1-4 号点位；2-4 号点位位于滑坡区东南侧。各点位分别比较了水平方向与垂直方向移动量。以传统测绘仪器获取的数值为真值，为了减小位置差错产生的误差，以监测点为中心，将一定范围内的SBAS点求取平均值作为SBAS 监测值。表3、表4 统计了在2018 年1 月至2021 年2 月时间段内二者比较结果。

表3 水平方向移动量比较

表4 垂直方向移动量比较

由表3、表4 可知：不管是用InSAR 技术还是传统测绘仪器技术，变形体位移下降趋势相同，位移量相近、数据相差较小。可以认为采用Stacking-InSAR 和SBAS-InSAR 技术的方式监测地表形变能达到亚厘米级的精度。产生误差的主要原因是两种方法监测点位置存在偏差，不能完全对应。采用InSAR 技术监测地面沉降，精度能够达到厘米级甚至毫米级别，能够满足监测项目需求。

5 结语

本文利用2014 年10 月至2021 年2 月期间覆盖新疆阿勒泰地区的多景Sentinel-1 数据，应用Stacking-InSAR 与SBAS-InSAR 技术相结合的方式获取了该地区该时间段的地表累计形变量图，分析得知该地区最大累计形变量达-550 mm，最大形变速率达-76.9 mm/a，可见最大变形体一直处于持续的形变状态。提出以下建议：

1）实地获取下游水电站相关数据，对该滑坡地质情况进行调查，做好防护工作与监测风险的情况的工作。

2）为保证变形体能够得到实时监测，可结合自动化监测技术，保证变形体全方位、多角度、可视化监测。