基于多阈值目标提取的时序InSAR矿区地表沉降监测研究

刘泽洲，卢才武，章 赛，李 萌，和郑翔

(1.西安建筑科技大学资源工程学院，陕西 西安 710055；2.西安建筑科技大学建筑学院，陕西 西安 710055)

矿区资源长期地下开采会导致地表沉降和形变，从而引发地表沉陷、裂缝以及滑坡等地质灾害，不但带来严重的经济损失,而且还会危害居民正常生活及人身安全。为了有效控制矿区灾害的发生，必须通过监测手段持续准确地对矿区地表变化情况进行监测，全面并及时地获取地表沉降情况。传统煤矿形变的监测方法虽然能精确获取测量点的形变信息[1]，但这种仅依赖有限离散观测数据的方法无法真实地反映地表形变的动态进程。 近年来，合成孔径雷达差分干涉测量技术(differential interferometric synthetic aperture radar，D-InSAR)以其突出优势发展迅猛[2-4]，尤其是多种时间序列InSAR方法的提出[5-7]，有效弥补了D-InSAR方法监测精度和可靠性受时间、空间去相干以及大气效应等因素影响的问题[8-10]。

受自身特殊条件影响，矿区地表沉降呈现出地形条件复杂、梯度大、沉降速率大、形变范围小等特点，使得监测难度较大[11-12]。短基线集(small baseline subsets，SBAS)技术能够有效解决矿区地表缓慢形变监测问题[13]。在SBAS技术处理过程中，一般会通过人工选取目标点进行轨道精炼与沉降反演[14]，这种方式虽然可以有效去除残留相位与平地效应，但若没有提前对研究区进行考查，获取地面目标点信息就盲目地选择，反而会引入误差[15]。由于矿区开采极易发生大量级地表沉陷现象，使得地表像元相干性较差，无法轻易准确获取稳定像元点，因此目标点的选取是当下矿区InSAR监测中的重点研究问题。

针对以上问题，本文提出一种基于多阈值目标提取的SBAS矿区地表监测方法。通过设定离差阈值参数、区域窗口阈值参数与相干性阈值参数提取最为稳定的目标点，并将其应用于SBAS-InSAR处理中。以兰州窑街矿区为研究区，探究多阈值目标提取的时序InSAR方法在矿区地表沉降监测中的应用效果。

1 方法原理

1.1 基于多阈值目标点的提取方法

地物目标的干涉相位中一般含有大量噪声相位，通过衡量噪声贡献程度即可选取干涉相位较为纯净的地物目标。SAR影像是由多个像元组成，假设任意像元的高斯噪声标准差为σn，相应的振幅∂信息为Rice分布，可得到式(1)。

(1)

式中：A为某一像元；I为Bessel函数；g为地面目标的正实数反射能量。

若该分布近似接近于高斯分布，则需要g/σn>4。此时，振幅离差指数可以近似为式(2)。

(2)

式中：σ∂为地物目标所对应的振幅标准差；σnI为虚部标准差；m∂为地物目标所对应的时序振幅均值；D∂为时序振幅离差指数。

相干系数数学表达式见式(3)。

(3)

式中：S1、S2为一组干涉对影像；M、N为局部窗口大小；i、j为提取像元的像素坐标。

本文提出的方法首先通过设定合理的离差指数阈值参数对地物目标进行首次筛选，提取干涉相位较为纯净的地物目标。其次，由于相干性能直接反应出数据处理生成的干涉相位质量的高低程度，因此设定相干性阈值参数对地物目标进行第二次筛选，消除大部分相干性差的像素点，最终获得探测结果。

在计算相干系数时，局部移动窗口大小的设定是保证相干系数精度的重要因素之一。一般情况下，窗口越大，可靠性越高，但在一定程度上降低了分辨率，易将附近的不稳定点误判为稳定目标点，且易遗漏离散分布单独存在的个别目标点。因此，为了平衡图像分辨率和相干系数估计精度这两个互为影响的因素，在进行第二次筛选之前对研究区域进行窗口划分。本文采用通过设定局部移动窗口阈值参数将整个区域划分成若干个子区域，再对每个子区域的地物目标进行相干性筛选，提取大于相干性阈值的目标点为参考点，最后对所有子区域进行镶嵌。

1.2 SBAS-InSAR方法原理

假设SAR影像在tA和tB(tA>tB)两个时刻干涉得到第j幅干涉图，在方位向和距离向的坐标为(x，r)，则在该点的相位可以表示为式(4)。

δφj(x)=φ(tB,x)-φ(tA,x)≈

(4)

式中：λ为雷达波长；d(tB，x)和d(tA，x)为雷达视线方向(LOS)的累计形变，d(t0,x)=0。

在(x，r)这一像元点，假设主影像IE、辅影像IS对应生成M期干涉图可表示为式(5)。

(5)

第j(j=1，2，…，M)幅干涉图的主影像获取时间较辅影像晚，即IEj>ISj，则其相位组成的观测方程见式(6)。

δφj=δφ(tIEj)-δφ(tISj)

(6)

式(6)表示含有N个未知干涉相位的M个观测方程，其矩阵形式见式(7)。

Aφ=δφ

(7)

A是一个M×N的近似关联矩阵，数据中生成的干涉图可以决定矩阵A。数据都在一个小基线子集里，当M接近或者小于N时，采用奇异值分解的方法可以求得该方程最小范数意义上的最小二乘解，最终获取时序形变量。

2 研究区概况与数据

2.1 研究区概况

窑街煤矿地处甘肃青海两省交界处，经度范围为102°50′～102°56′，纬度范围为36°20′～36°26′，行政区划属兰州市红古区，该区域地处黄土高原西部，海拔1 800～2 400 m，地形较复杂[15]。

该矿区现有窑街三矿、金河煤矿、海石湾煤矿等三个开采矿区。由于多年来不断的地下开采，矿区内地表沉陷、滑坡等地质灾害频发，严重影响矿区内各类建设设施以及采矿生产安全[16]，研究区信息如图1所示。

图1 研究区影像图Fig.1 Image map of the study area

2.2 实验数据准备

本研究选用具有较短时间及空间基线的哨兵一号卫星(Sentinel-1 A)干涉宽幅模式(IW)下的单视复数数据，该卫星搭载5.6 cm的C波段合成孔径雷达，重访周期为12 d，影像分辨率为5 m×20 m。本文选取的实验影像时间跨度为2018年1月—2020年4月，影像雷达入射角为39.06°，方位向和距离向分辨率为13.91 m与3.70 m。此外，为提高影像配准的精度，加入精密轨道星历数据降低轨道误差。参考DEM数据采用90 m分辨率SRTM-DEM数据。

3 实验数据处理

3.1 多阈值目标点提取

传统SBAS-InSAR方法采用人工选取目标点，但在未掌握地面目标点信息的情况下进行选择很有可能会增大误差。因此，本文提出一种基于多阈值的目标点提取方法。对原始单视复数数据进行多视配准之后，得到研究区强度信息图，采用本文提出的方法对所有影像数据进行分析，提取地面目标点。首先，为了最大程度提取最优点作目标点，设定离差阈值参数为3.2，进行第一次筛选，结果如图2(a)所示，实验共筛选出目标点29 453个。其次，由于研究区窗口较大，计算像元相干性时，在较高相干性区域易产生簇成团的控制点，在较低相干性区域易将附近的不稳定点误判为控制点，且易遗漏离散分布单独存在的个别控制点。因此，为了平衡高图像分辨率与高相干系数精度，设定局部窗口阈值参数将整个区域分解成若干个子区域，再进行相干性计算，经过反复实验，设定区域窗口阈值为16 sqkm，将整个区域划分为132个子区间，在每个子区间内进行筛选。最后，设定相干性阈值参数为0.9，在第一次筛选结果的基础上进行第二次筛选，筛选后的地面控制点如图2(b)所示，实验最终总共筛选地面控制点103个，将最终筛选结果作为后续轨道精炼的目标点用以进行地面沉降反演。

图2 筛选的目标点Fig.2 Target point for screening

3.2 改进SBAS-InSAR处理

对研究区数据进行SBAS-InSAR处理，空间基线阈值设定临界极限的5%，时间基线阈值设定为180 d，得到的时间和空间基线连接图如图3和图4所示。

图3 空间基线连接图Fig.3 Diagram of spatial baseline connection

图4 时间基线连接图Fig.4 Diagram of time baseline connection

在干涉处理结束后，检查所得到的所有干涉像对，去除相干性差的像对，避免形变提取时产生误差。而后将基于多阈值方法提取的目标点代替人工选点方法，用作SBAS-InSAR轨道精炼与重去平处理。最后，利用奇异值分解法估算形变速率与残余地形，并通过空间域低通滤波算法与时间域高通滤波算法去除大气相位与残余地形相位，获取时序形变速率。

4 结果分析

4.1 沉降速率分析

经过上述处理后，得到研究区2018年1月—2020年4月期间的雷达视线方向(LOS)年平均速率图。由式(8)可以得到研究区年平均垂直沉降速率情况，如图5所示。

Δh=Δr/cosθ

(8)

式中：θ为雷达入射角；Δr为雷达视线方向沉降量；Δh为垂直方向沉降量。

图5 研究区年平均垂直沉降速率图Fig.5 Annual average vertical settlement rate map

由图5可知，研究区内大部分地区地表沉降较为稳定，年平均沉降速率主要集中在-24～13 mm/a。A、B、C三处区域存在明显地表沉降，具体表现为由边缘区域向中心逐渐递增的沉降漏斗状，沉降速率为-25～-156 mm/a沉陷和14～87 mm/a抬升，这与文献[16]得出的结果基本一致。结合光学影像资料对比，图5中三处沉降漏斗与实际该煤矿开采工作面基本吻合。图6为沿HD方向剖面侧视图，结合矿区地下开采资料进行对比分析可知，矿区地下开采是导致三个沉降区发生沉陷的主要原因。

图6 HD方向剖面侧视图Fig.6 Side view of section in HD direction

4.2 时序沉降分析

为分析矿区地表沉降时序演变规律，对整个研究区进行时间序列分析计算后，获得研究期内各个时间点相对于第一期影像(2018年1月2日)的沉降累计量，取其中较为重要的11个时间节点沉降累计图(图7)。由图7可知，2018年5月之前研究区没有明显沉降，沉降量集中在-5～-15 mm之间。从2018年5月开始，研究区出现较为明显沉降，并且随着时间推移，沉降值和沉降范围逐渐增大，到2020年4月沉降值达到-376 mm。其中，2018年9月、2019年1月、2020年4月沉降程度相对较大，沉降值较前分别有58 mm、54 mm和45 mm的增长。

图7 研究区时序沉降量Fig.7 Time series subsidence in the study area

为了进一步分析开采沉陷区沉降发育过程，在研究区内选取三处(A区、B区、C区)沉降最为明显的区域，对其做时序沉降趋势分析，分析结果如图8所示。与研究区煤矿地下开采资料对比分析，三处区域分别位于窑街煤矿区域中的窑街三矿、金河煤矿、海石湾煤矿三个开采矿区。三处矿区在整个研究时段内大体呈现线性下沉趋势，其中A区和B区整体下沉速率较为稳定，C区呈阶段性变化趋势，在2018年3月之前下沉趋势较为严重，沉陷值可达-78 mm，2018年4月之后下沉趋势较为稳定。

图8 时序沉降趋势图Fig.8 Time series subsidence trend

4.3 监测结果可靠性分析

为验证多阈值目标提取的SBAS方法的有效性，使用相同的原始数据，利用传统SBAS-InSAR方法获取A区、B区、C区的时序地表沉降值，将两方法监测结果进行对比分析，对比结果见图9。

图9 两种方法监测结果对比Fig.9 Comparison of monitoring results of two methods

由图9可知，A区两种方法监测结果中差值最大为9.30 mm，最小为0.12 mm，差值的平均值为0.76 mm。对24期数据的差值离散程度进行计算，得到标准差范围在0.06～4.66 mm之间，标准差平均值为2.06 mm。B区两种方法监测结果中差值最大为18.80 mm，最小为0.39 mm，差值的平均值为6.80 mm。对24期数据的差值离散程度进行计算，得到标准差范围在0～9.40 mm之间，标准差平均值为3.96 mm。C区两种方法监测结果中差值最大为42.20 mm，最小为0.80 mm，差值的平均值为28.70 mm。对24期数据的差值离散程度进行计算，得到标准差范围在0.40～21.10 mm之间，标准差平均值为12.30 mm。

综上所述，基于多阈值目标提取的SBAS方法与传统SBAS-InSAR得到的监测结果大致相吻合，基于以往SBAS-InSAR方法的监测精度，可认为本文提出方法得出的结果是可靠的。

5 结 论

由于矿区地表环境复杂，很难有效提取稳定的目标点用以SBAS处理。因此，本文提出了一种基于多阈值目标提取的时序InSAR方法，利用此方法反演出研究区的年平均沉降速率和时序累计沉降值等时间序列沉降信息，并与传统SBAS方法监测结果进行实例对比。研究得出结论如下所述。

1) 两种方法监测结果对比分析，三处矿区沉降差值平均值分别为0.76 mm、6.80 mm、28.70 mm，标准差平均值分别为2.06 mm、3.96 mm、12.30 mm，两种方法监测结果较为吻合，证明基于多阈值目标提取的SBAS方法与传统SBAS-InSAR方法在矿区地表沉降提取精度均可达到亚厘米级。

2) 传统SBAS-InSAR方法为保证提取到更多地表相位较为稳定的目标点，在选择研究区方面表现较为苛刻，这使得其在实际应用过程中存在许多局限性。而本文提出方法可在保证精度的同时，克服传统SBAS-InSAR方法选取目标点在实际应用中存在的局限性，尤其在矿区发生大量级地表沉陷现象，导致地表像元较差的情况下，本文提出方法具有更好的适用性，可为后续矿区沉降监测分析提供新的监测手段。

3) 监测结果表明，研究区范围内有三处开采沉陷盆地，分别位于三矿(A区)、金河煤矿(B区)、海石湾煤矿(C区)。最大沉陷值为-376 mm，最大年平均沉降速率为-156 mm/a，且沉降下沉趋势在逐步增加，因此需对该矿区进行长时间序列的监测，用以提供灾害预警。