SBAS-InSAR在矿区地表形变监测中的应用与分析

王道顺，李伟*，孙建，郑国栋

(1.山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队)，山东 兖州 272100；2.高分辨率对地观测系统济宁数据与应用中心，山东 兖州 272100；3.自然资源部采煤沉陷区综合治理工程技术创新中心，山东 兖州 272100)

0 引言

煤炭资源作为我国的重要能源之一，极大地推动了国家的经济和社会发展。随着开采强度的日益增加，由此引起的生态环境问题更加凸显[1-3]。煤炭开采造成的地表沉陷属于典型地质灾害之一，往往导致矿区及一定范围内的地表产生不同程度的沉降，对地面房屋设施、道路、农田水利基础设施等建(构)筑物造成不同程度的损毁，为经济发展和社会稳定带来重大隐患，因此对煤炭开采引起的地表沉降进行监测具有重要意义。

以往的矿区地表沉降监测多采用布设岩移观测站的方式，利用水准测量、GNSS测量等传统方法进行监测，此类方法监测精度可观，但作业成本较高，获取的监测数据量有限，不能完全反映矿区地表沉降特征，难以实现大范围的地表沉降监测。合成孔径雷达差分干涉测量(D-InSAR)通过微波遥感主动测量的方式，可以克服传统监测方法的不足，且具有全天时、全天候作业的优点，可获取大范围、高精度的空间延续性地表形变信息。已有学者将其应用于矿区的地表形变监测研究[4-9]，但该方法使用数据量较少，在长时间序列的监测过程中容易受到时空失相干、大气延迟等因素的影响，降低其监测精度。在此基础上发展而来的时序InSAR技术弥补了上述不足，极大地提高了形变监测精度，已广泛应用于矿区形变监测、城市沉降监测等方面，其中以永久散射体技术(Permanent Scatterer InSAR，PS-InSAR)和短基线集方法(Small Baseline Subset InSAR，SBAS-InSAR)应用最为普遍[10-15]。

PS-InSAR技术已成功应用于矿区的地表沉降监测，为了获得在时间序列上散射特性稳定的PS点，PS-InSAR往往需要20景以上的数据量参与计算。该方法基于大量PS点的迭代回归计算，导致运算效率偏低[16]。矿区多位于郊区、农田等低相干地区，若所用数据量偏少、时间间隔较长，利用PS-InSAR技术难以获取足够的PS点，导致无法完整揭示矿区沉降规律。本文在分析SBAS-InSAR技术方法的基础上，以济宁北部矿区为研究区域(图1)，采用C波段Sentinel-1A雷达数据和SBAS-InSAR时序方法开展矿区地表沉降监测和时序分析，获取监测范围内的地面沉降时空分布特征，可为采煤沉陷区的调查和治理提供技术支撑。

图1 研究区域影像图

1 SBAS-InSAR原理

SBAS-InSAR由Berardino等人[17]于2002年提出，不同于PS-InSAR的单主影像，该方法是一种基于多主影像的InSAR时间序列方法。它通过短基线原则，将大量SAR数据组合为具有多个主影像的干涉子集，每个子集内的干涉对基线长度均低于临界基线值，时间基线也尽可能短，集合间的SAR影像基线距大，通过这种方式克服了时间和空间上的失相关，因此SBAS-InSAR可以通过较少的数据量来获取较可靠的监测结果，且其监测对象为分布式散射体，该方法更适合于矿区的地表形变监测。

基于短基线集原则，SBAS-InSAR具有多个主影像，但仍需选取其中一景影像作为公共主影像进行配准；组成各干涉子集后，利用外部参考数字高程模型(Digital Elevation Model，DEM)数据模拟并去除每个干涉相对的地形相位，然后生成时间序列差分干涉图集；相位解缠后得到每个相干目标的相位信息，包括形变相位、大气延迟相位、轨道误差相位等信息，各误差相位可在时间序列上采用滤波方法或多项式模型予以去除；由于SBAS-InSAR具有多个主影像，各干涉子集在联合求解时容易出现方程秩亏现象，因此引入奇异值分解方法，利用最小二乘原理得到地表时间序列形变信息。

2 研究区域与数据源

2.1 研究区域

济宁市位于山东省鲁西南地区，矿产资源丰富，是鲁西煤田基地的重要组成部分，是全国重要的煤炭基地之一[18-20]。本次实验研究区域位于济宁市北部地区，覆盖范围大致为E116°28′01″～116°44′34″、N35°29′40″～35°42′35″，研究区域内矿区集中分布，主要有义桥煤矿、唐阳煤矿、义能煤矿、鲁西煤矿、新驿煤矿等矿区(图1)。

2.2 数据源

本文使用覆盖研究区域2019年9月19日—2020年3月29日间的17景Sentinel-1A C波段单视复数(SLC)数据，该数据为干涉宽幅(IW)成像模式，极化方式VV，空间分辨率5m×20m(方位向×距离向)，具体信息见表1。Sentinel-1A卫星在轨运行高度693km，重访周期12d，若与Sentinel-1B卫星组网观测，可将重访周期缩短至6d，可有效提高监测区域的时间相干性，此外由于Sentinel-1A卫星采用先进的轨道控制技术，其干涉对空间基线最大为162m，多数都在100m范围内[21]，远远小于其空间临界基线值，因此可以忽略其干涉对空间基线的差异。外部参考DEM选择由日本太空发展署(JAXA)免费分发的ALOS World 3D数据，空间分辨率30m，用于模拟去除地形相位和地理编码。

表1 研究区Sentinel-1A影像参数

3 矿区形变监测与分析

3.1 短基线集干涉组合

实验选择2019年11月30日的影像作为公共主影像，其余影像作为从影像与主影像进行配准并采样到主影像像素空间。设置空间基线阈值100m，时间基线阈值60d，共得到55组干涉相对，其中最大空间基线87.47m，最小空间基线4.11m，时空基线连接情况见图2。

图2 时空基线图

3.2 干涉图生成与相位解缠

将得到的55组像对进行差分干涉处理，得到时间序列干涉图集，同时生成对应的相干系数图，由于设置的阈值相对严格，多数差分干涉图干涉效果较好，可以看到明显的形变干涉条纹。图3是2019年11月30日—2019年12月24日内生成的差分干涉结果(SAR坐标系)，图3a为相干系数图，颜色越亮代表相干性越高；图3b为差分干涉图，可以看出，监测区域内有多处明显的形变干涉条纹。相位解缠选择最小费用流(MCF)解缠算法，相干性阈值设为0.2，低于相干阈值的像元进行掩膜处理，不参与计算，然后逐个干涉图进行相位解缠，图3c为对应的相位解缠图，黑色代表相干性过低而没有解缠结果的地方。

3.3 轨道精炼和重去平

对得到的55组干涉像对结果进行目视检查，剔除相干性过低或解缠错误的像对。选取一组干涉效果良好的像对作为参考影像，选择远离形变区域且相干性较高的若干稳定地面点作为控制点，使用多项式模型进行轨道精炼，去除残余的恒定相位和相位坡道。

3.4 获取监测结果

上述步骤完成后反演区域形变速率，同时对干涉图进行二次解缠优化。在第一次形变速率反演的基础上，利用大气相位在时间上的高通滤波和空间上的低通滤波进行去除，获取时间序列上的位移结果(图4)。

a—相干系数图；b—差分干涉图；c—解缠相位图图3 干涉对相干系数图、差分干涉图、解缠相位图

图4 监测区域年均沉降速率图

3.5 结果分析

从监测结果可以看出，监测区域整体基本处于稳定状态，济宁北部矿区内明显分布有6处不同沉降程度的区域(A—F)，分别位于义桥煤矿、唐阳煤矿、义能煤矿、新驿煤矿、鲁西煤矿、何岗煤矿，其中C区域沉降最大，年平均沉降速率达到-242mm/y，呈现出双沉降漏斗特征；A、B两处沉降区域相距仅640m，已呈现接连成片沉降趋势，沉降范围最大；D区域最大沉降速率达到-85.7mm/y；E，F两处区域沉降范围最小，E区域最大沉降速率为-70.1mm/y，F区域最大沉降速率达到-91.6mm/y。部分沉降区域中心不同程度的出现了“监测空值”，主要是矿区沉降过大，超出了InSAR可探测的最大形变梯度，致使沉降中心出现失相干，但其沉降边缘范围探测不受影响。

4 结论

(1)基于2019年9月19日—2020年3月29日的17景Sentinel-1A数据，运用SBAS-InSAR技术获取了济宁北部矿区的年均沉降速率和沉降漏斗分布。

(2)济宁北部矿区内分布有6处明显沉降区域，分布于义桥煤矿、唐阳煤矿、义能煤矿、新驿煤矿、鲁西煤矿和何岗煤矿，形变范围与矿区实际空间分布位置一致，其中义能煤矿沉降特征最为明显，年均沉降速率达到-242mm/y。

(3)受形变梯度影响，部分沉降区域中心出现“监测空值”，但沉降范围仍可以完整探测。相比较于C波段数据，L波段数据波长更长，在低相干区域更容易穿透地表植被覆盖，得到更高质量的干涉结果和更大中心沉降量。下一步工作将收集该区域L波段数据，开展更高质量的InSAR矿区地表形变监测，为矿区的地面沉降防治提供决策依据。