应用StaMPS-PS 监测惠州地表沉降时空演化

罗 莉 ， 王 斌

（1. 广东省测绘工程公司， 广州 510670； 2. 广东省国土资源测绘院， 广州 510500）

0 引言

地面沉降是一种典型的地质灾害[1]。 当自然或人为因素诱发地面垂直形变超过一定量时， 会造成滑坡、 崩塌、 地貌塌陷、 泥石流等多种灾害破坏[2]。 传统地面沉降监测技术如水准测量、 传感器网络监测[3]等技术手段， 点位测量精度优于毫米级， GPS 高精度似大地水准测量精度达到或者优于厘米级[4-5]， 在建筑变形测量、 边坡监测、 区域地表沉降[6]等领域发挥重要作用。 对于大范围地面沉降监测， 需耗费大量人力、 物力资源， 劳动强度大， 精度不均匀。 InSAR 是一种新兴的空间大地测量技术[7]， 目前已被广泛用于地表沉降监测， 自然灾害 （如地震、 火山、 滑坡、 水土流失等）监测及环境监测（如海上溢油、 土地利用变化等）等民用、 军用领域。 相较传统监测手段， 具有监测范围广、 监测效率高、 监测网络空间密度高、 提供整体时序形变等特点[8]。 刘一霖等利用短基线集INSAR 技术（SBAS）分析黄河三角洲地表形变特征场和时间序列特征， 对地下水抽取等人为因素进行了相关分析[9]， 但难以获取有效地表INSAR 信息， 对沉积物固结压实产生的地面沉降因素未纳入相关分析。 开展大范围沉降监测和预警， 陈炳乾等基于D-INSAR 和永久散射体雷达干涉技术（PS-INSAR）获取和分析开采区沉陷影响和发展趋势[10-12]， 不需要先验知识获取沉陷区影响范围和发展趋势， 但也存在在有植被的地区提取的PS 点目标形变结果可靠不足， 在未联合GPS 连续运行参考站数据解算情况下， 还难以全面掌握精细的地表三维形变场。

本文采用StaMPS-PS（Stanford Method for Persistent Scatters PS-INSAR）技术[13-14]， 分析研究区地表沉降变化。 通过振幅离差法和相干系数法筛选出相位稳定的PS， 然后将PS 点的相位重新采样到格网上， 利用形变相位、 大气延迟相位、 轨道误差相位在时间域和空间域上频谱特性的不同， 通过一定的滤波手段分离出形变相位和其他误差相位，获得精度较高的地表形变场。

1 研究区与数据预处理

1.1 研究区域

惠州市地势总体北高南低， 地貌类型复杂多样， 有山地、 丘陵、 台地、 平原和湖泊， 北部多为山地和高丘陵， 南部则为平原和台地。 构成惠州市各类地貌的基岩岩石以花岗岩最为普遍， 砂岩和变质岩也较多， 此外局部还有红色岩系地貌，沿海沿河地区多优质沙滩及复杂形态的珊瑚礁， 地质构造多为第四纪沉积层。 这些不同的地貌类型及地质构造导致惠州市地表沉降变化复杂， 尤其是连续性地表沉降可以发展为一种严重的地质灾害， 能持续对道路、 房屋建筑、 地下管线等造成破坏。 尤其位于珠江三角洲平原区的惠州市是广东省经济发达地区之一， 人口稠密， 经济繁荣， 但近几年出现的地面沉降给当地公路、 桥梁、 水利设施、 地下管网设施等造成了不同程度的破坏。 如图1 所示。

表1 Sentinal-1A 主要参数Table 1 Main parameters of Sentinal-1A

1.2 数据预处理

1.2.1 SAR 影像数据

欧空局Sentinel-1A 作为欧洲空间局哥白尼计划的首颗环境监测卫星， 于2014 年4 月3 日发射升空。 在 2015 年 4 月至 5 月期间， 卫星系统稳定运行后开始采用12 d 的重访周期进行不间断全球成像任务。 成像模式主要有4 种： 条带模式（Strip Map， SM）、 干涉宽幅模式（Interferometric Wide，IW）、 极宽幅模式（Extra-Wide Swath， EW）和波谱模式（Wave）。 本项目采用 IW 模式下的 VH 数据， 具体参数如表1。

选取合适的SLC 影像作为原始SAR 影像。 IW模式下 Sentinel-1A 采用 TOPS 技术， 一景 SAR 影像可以成像三个子条带及多个Burst， 根据惠州市范围对原始影像进行裁剪及拼接， 形成最终符合要求的整幅单视复数影像。 如图2 所示。

图1 惠州市辖区分布图Fig.1 Distribution map of Huizhou City

图2 惠州市单视复数影像Fig.2 The single-view complex image of Huizhou City

1.2.2 DEM 数据

采用 ALOS Global Digital Surface Model， 30 m分辨 率 DEM 高程 数 据。 AW3D30（ALOS World 3D－30 m）是日本宇宙航空研究开发机构（JAXA）于2016 年5 月发布的全球数字表面模型数据集， 水平分辨率约30 m。 由陆地观测卫星ALOS 上搭载的PRISM 立体相机对全球陆地±80°纬度地区进行覆盖观测。

图3 惠州市地理编码后DEMFig.3 The geocoding DEM of Huizhou City

AW3D30 存在一定范围内的无效数据值， 在处理过程中， 对原始的DEM 进行缺失数据的探测， 沿海范围有少量的数据缺失， 将其高程置为0值。 然后对预处理后的DEM 数据进行拼接及地理编码， 使其与SAR 影像统一坐标系。 如图3 所示。

2 时序StaMPS 模型及处理方法

2.1 时序StaMPS

StaMPS 技术是基于统计方法筛选PS 点， 利用三维解缠算法， 首先在时间上计算每个PS 点的相位差异， 然后设置参考点用最小二乘法在空间上进行解缠。 对于空间相关误差去除， 利用相位解缠使PS 点的相位值恢复到绝对值， 但同时除了包含形变相位， 也包含了大气相位、 卫星轨道误差、DEM 残差相位、 随机噪声相位， 而这些相位在时间上是不相关的， 采用时间域上的高通滤波分离出形变相位， 并用空间域上的高通滤波分离出辅影像的大气相位和轨道误差。 其具体算法流程如下：

（1）在N 幅影像中选取1 幅作为主影像， 其余影像均与主影像配准， 并进行辐射定标， 将副影像的几何信息和幅度信息与主影像匹配。

（2）作干涉及差分处理， 差分主要去除平地相位， 并利用外部DEM 去除地形相位， 得到N-1 幅干涉图和差分干涉图。

（3）利用 M 幅影像的幅度信息， 设置阈值， 进行PS 点的选取。

（4）进行相位解缠， PSInSAR 较为常用的空间域解缠方法为稀疏格网解缠。 估计相邻PS 点的缠绕相位梯度， 对梯度进行积分， 空间搜索， 进行解缠。

（5）计算相应的形变量、 DEM 误差， 并从残差相位中分离大气扰动和非线性形变量。

（6）得到各个分量之后， 进行后处理工作， 如克里金插值， 最终得到整个实验区的形变场。

2.2 三维解缠及最小二乘解算

StaMPS 中三维解缠分别体现在空间（二维）和时间（一维）上， 基于最小化L-P 范框架， 将相位解缠看做最优化问题， 找到目标函数的最小值，从而获取解缠后的相位值。 目标函数为：

图4 StaMPS 处理流程Fig.4 The processing flow of StaMPS

式（1）中， ΔΦ 为解缠相位梯度， Δψ 为缠绕相位梯度， x 和y 分别表示二维的两个方向， z 表示第三维方向， i 和 j 表示点坐标， w 表示权重。 和传统L-P 范分布不同之处在于， Hooper 在StaMPS加入了第三维度的相位梯度信息， 通常为时间序列信息。 这里， 参数P 决定如何处理ΔΦ 和Δψ 的差值关系， 当P=0 时， 若仅考虑二维信息， 上式为类似枝切解缠算法的L-0 范数目标函数； 当P=2 时，即为最小二乘解缠算法。 在解缠完成并去除空间不相关的误差后， 构建观测方程， 加入最小二乘准则， 求取线性形、 主影像的大气延迟和轨道误差以及空间相关的 DEM 误差

参数估值及协方差阵可以解得：

3 试验结果及分析

3.1 平均形变速率

惠州市覆盖范围跨度较大， 若在不做多视的情况下整体处理， 对内存的要求很高且时间效率低。 本文对原始影像裁剪分为龙门县、 博罗县、惠城、 惠阳区和惠东县分别进行处理。 每个区域的控制点选取稳定的水准点作为参考， 惠州市整体形变速率图根据各区域的形变结果进行拼接，选取了龙门 277 092、 博罗 974 064、 惠城 794 967、 惠阳 1 001 252、 惠东 487 997 个 PS 点， 获取其形变结果。 试验发现： 2015 年 12 月至2018年2 月期间， 惠州市地表平均形变速率整体形变较小， 大部分区域在-3～3 mm/y 之间， 龙门县南部地区和博罗县城区有部分下沉（红圈范围），最大形变速率达到-12 mm/y （如图 5 所示）。 惠城、 惠阳区及惠东县近海域有部分地区出现轻微的抬升，但不超过3 mm/y （蓝圈范围）。 由于C 波段数据在多植被覆盖区相干性较差， 惠东县东部选点较少，城区部分地形较稳定， 基本无形变趋势， 山区有部分地区存在沉降。 博罗县和惠阳区稳定性较差，相对沉降速率较高； 惠东县近海域地区以及惠城区有少量抬升情况， 如图6 所示。

图5 惠州区域平均形变速率Fig.5 The average deformation rate of Huizhou

3.2 形变序列

2015 年 12 月至 2018 年 2 月期间， 惠州市地表整体形变趋势随时间变化总体较稳定。 累积沉降量 （年平均形变速率与监测期乘积）最大为-24 mm， 抬升量最大为 8 mm， 如图 7 所示。

3.3 结果讨论与分析

图6 惠州市形变速率Fig.6 The deformation rate of Huizhou

（1）利用StaMPS 技术对惠州市地表沉降进行分区监测， 获得了沉降速率、 时序累计沉降值。根据监测结果， 龙门县南部地区和博罗县城区有部分下沉情况， 最大形变速率达到-12 mm/y， 应该予以重视。 博罗县最大累积形变量量达-24 mm、惠阳区北部万里工业区附件形变速率达到-6 mm/y，稳定性较差， 相对沉降速率较高， 对于这些地区的地表状况应予以关注； 惠东县近海域地区以及惠城区速率最大达到4 mm/y 有少量抬升情况。

（2）重点监测惠州市整体地表时序沉降， 并分析整个目标区内的沉降时空特征及变化， 对惠州市地表沉降做出预判， 反映地表形变的动态变化信息， 为城市地表沉降防治提供依据。 为惠州市经济社会发展提供科学可靠的地理国情信息服务。

（3）本项目使用的SAR 卫星影像较为单一， 主要是Sentinel-1A 的C 波段数据， 没有使用不同卫星对具有不同地表特征的区域分别进行形变监测分析比较。 Sentinel-1A 数据介于X 波段和L 波段之间， 虽然均衡性较好， 但是在城区监测方面不如X 波段数据； 另一方面， 惠州市龙门县和惠东县地貌多变， 植被覆盖范围较广， C 波段数据失相干严重， 造成了选点过少因此无法有效监测的问题。 建议对城区采用X 波段数据， 对山区采用L波段数据能更加有效精确地反演地表形变特征。

（4）本项目主要使用的StaMPS 的技术， 该技术在选点和解缠方面优势较大， 但是StaMPS 技术是以线性模型求解形变相位和沉降速率， 若实际形变并不是以线性模型沉降， 这会导致StaMPS 技术得出的沉降值偏小。 建议在选点和解缠方面可以利用该技术， 后续的解算可以考虑与其他时序技术联合处理。

图7 惠州市时序形变速率图（部分）Fig.7 The time series deformation rate diagram of Huizhou City （part）