基于SBAS-InSAR 技术的杨伙盘矿区地表形变监测

张 帆，常 乐，荀张媛

（1.西安煤航遥感信息有限公司，陕西 西安 710100；2.沈阳城市建设学院，辽宁 沈阳 110167）

0 引言

我国是矿产资源大国，煤炭等资源一直在能源开发利用中占据重要地位。煤炭资源的长期开采不仅引起地面塌陷等一系列问题，还会造成建筑物受损、农田被破坏[1]，以及诱发崩塌或滑坡等次生灾害[2]，对人民的生产生活造成严重的经济损失，甚至会危及生命安全[3-4]。因此，重视矿区地表形变，实现矿区生态环境可持续发展是科学开采煤炭等矿产资源的重要任务。

矿区地表形变是一个基于面区域的大范围的地表运动，但是在实际生产过程中，更多的还是依赖于传统的水准测量、GPS 测量等点监测技术，在生产工作面上每隔一段距离布设一个水准点或GPS 点。一方面，布设较多的点，人工消耗成本高、周期长；另一方面，水准测量和GPS 测量都是基于点观测，对于矿区大范围大区域监测，需要大量加密点，通过拟合获取最终结果，无法保证精度[5-7]。

合成孔径雷达干涉测量（InSAR）技术是一种高效对地观测技术，具有全天时、全天候对地成像，监测结果精确度高，测量范围广的优势，大大降低人力、物力、财力的消耗[8-10]。特别是时间序列InSAR 技术能够更快速获取地表大范围高精度时间序列形变信息，为矿区开采沉陷监测提供了新的方向[11-13]。

越来越多的学者已经将InSAR 技术应用于煤矿地表监测研究。尹宏杰等[14]利用时间序列InSAR 技术对湖南冷水江锡煤矿地面形变进行监测，证实了煤矿地表时间序列形变漏斗的发展过程；赵伟颖等[15]利用SBAS-InSAR 技术获取了徐州矿区地表的时间序列形变及累积形变量，并与水准测量结果进行了外部精度符合验证，且两者的相关性较高；张香凝等[16]采用SBAS-InSAR 技术获取了煤矿地表形变的时空分布特征，并利用数值模拟技术对观测形变结果进行模拟分析，获取了地面沉降在时间和空间上的变形规律和机制；姚佳明等[17]选用升轨、降轨L波段PALSAR-2 影像数据，利用InSAR 技术对煤矿采空区开展了短期动态地表沉降监测，结合研究区开采信息对煤矿采空范围及开采时间进行反演，验证了InSAR 技术对煤矿采区反演的合理性与可靠性；何荣兴等[18]介绍了采空区灾害类型及不同类型的相应案例，分析了采空区灾害发生的一般规律和特征，提出了尽量选择不产生采空区的采矿方法。

本文以杨伙盘煤矿为例，采用SBAS-InSAR 技术，利用中分辨率Sentinel-1A 数据对煤矿地表进行时间序列形变监测分析，提取矿区开采沉陷形变场，并结合光学影像数据，研究其沉陷变化规律。

1 研究区概况

陕北煤炭生产基地是全国煤炭生产基地之一，尤其是神木、府谷等区域，煤炭资源丰富。但是长期高强度的资源开采，加之陕北地区大多为黄土地貌，最终导致矿区的土地大面积沉陷，并引发了滑坡、崩塌等多种次生灾害。因此，为保护矿区生态环境，减少人民的财产损失，需要对矿区的采空区形变进行长时间的有效监测；结合多种技术，加大对采煤塌陷区域的监测监管力度，形成一套行之有效的陕北榆神府矿区采空区沉陷监测体系。

杨伙盘煤矿位于陕西省神木市城北约30 km 的府店一级公路北侧，行政区划属神木市店塔镇，靠近府谷县（图1），工作区面积约26.92 km2。矿区地表广泛覆盖着现代风积沙及第四系黄土层、新近系红土层等，仅部分地区有基岩出露，基岩以砂岩、泥岩为主。根据某一钻孔资料，3-1煤层主要位于侏罗系中统延安组第三段，煤层距地表约160 m。根据现有巷道布置情况，结合各煤层开采范围及煤层赋存特点，一水平3-1煤层划分4 个盘区。图2 展示了3-1煤层开采工作面分布情况。

图1 杨伙盘煤矿地理位置及高程Fig.1 Geographical location and elevation of Yanghuopan Coal Mine

图2 工作面分布图Fig.2 Distribution map of working face

2 InSAR 数据及技术方法

2.1 InSAR 技术方法

BERARDINO 等[19]提出的SBAS-InSAR 技术，是目前较为常用的一种时间序列InSAR 技术，技术流程如图3 所示。首先，根据SAR 数据量，通过设置时间基线阈值和空间基线阈值组合干涉对，一般根据相干系数图和干涉图挑选质量较高的干涉对；其次，通过多视和滤波等操作提高所选干涉图的信噪比，基于具有稳定后向散射特性的高相干点去除地形相位和平地相位，并进行解缠；之后，根据干涉相位及时空基线去除DEM 误差，并分离线性形变相位，解缠残余相位，通过时间域高通滤波和空间域低通滤波分离出残余相位中的大气相位，剩下即为非线性相位，补偿线性形变相位即可得到完整的形变相位；最后，采用奇异值分解，连接多个子集，求得影像序列地表形变速率的最小范数、最小二乘解[20]。该方法能有效去除或减弱时空失相干等误差的影响，并在一定程度消除大气相位和限制地形误差的影响。

图3 SBAS 技术流程图Fig.3 Flow chart of SBAS technology

获取研究 区t0······tn时间段内N+1幅SAR 影像，将影像配准到同一个坐标系，设定时间基线阈值及空间基线阈值，得到M个差分干涉对。假设其中第i幅干涉图由tA、tB（tA＜tB）两时期影像生成，此干涉图上雷达坐标系下（a,r）处的干涉相位可表示为式（1）。

式中：δφi（disp）为形变相位；δφi（topo）为残余地形相位；δφi（aps）为大气延迟相位；δφi（noise）为噪声相位。忽略残余地形相位、大气延迟相位及噪声相位，将形变相位以LOS向形变表示，式（1）可简化为式（2）。

式中：λ为雷达波长；dtA（a,r）、dtB（a,r）为（a,r）处在tA时期、tB时期相对于初始t0时期的形变。针对（a,r）这一像元处，假设IE为主影像，IS为从影像，对应全部的M个干涉对，见式（3）。

式（3）满足IEi＞ISi，其中，i=1,···,M。则对于任意干涉图，对应的有式（4）。

式（4）为n个未知数M个方程所组成的方程组，以矩阵形式表示 δφ=Aφ，其中，A为M*n矩阵。对于小基线子集内部，M≥n，方程可由最小二乘法求解。各子集之间M接近或者小于n时方程会出现多解，采用奇异值分解以及最小范数进行解算获得唯一解[19]。

2.2 SAR 数据

本次研究使用的数据为欧洲航天局Sentinel-1A卫星在TOPS 成像模式下的宽幅干涉SLC 数据，为目前常用的免费中分辨率SAR 数据，数据产品详细参数见表1。根据监测周期，共收集了2020 年8 月—2022 年11 月共58 景升轨数据，轨道号为113/126。

表1 Sentinel-1A IW 产品基本参数Table 1 Basic parameters of Sentinel-1A IW product

本次采用的DEM 数据为ALOS Global Digital Surface Model“ALOS World 3D-30 m”（AW3D30），是由日本宇宙航空研究开发机构（JAXA）2015 年5 月免费提供的高精度全球数字地表模型数据，水平分辨率为30 m（1"），高程精度5 m，是目前世界上最精确的3D 地图，覆盖全球所有的土地尺度，工作区的DEM 数据如图1 所示。精密轨道数据采用成像21 d之后发布的POD 精密轨道数据。

3 结果与分析

3.1 InSAR 形变分析

工作区年平均形变速率图如图4 所示。形变部分主要分布在①、②、③三个区域，沿着工作面均表现为长条状分布，主要覆盖的工作面有30112 工作面、30114 工作面、30116 工作面、30115 工作面、30117工作面、30301 工作面和30302 工作面。三个形变中心可探测的最大形变速率分别约为-52 mm/a、-39 mm/a 和-39 mm/a，形变最大的区域位于区域①，分别对三个形变区进行时间序列和剖线形变分析。

图4 工作区年平均形变速率图Fig.4 Annual average deformation rate of the working area

图5 为区域①的时间序列形变图。由图5 可知，从2020 年8 月开始，30112 工作面处于开采状态，图5（a）为基准面，在图5（b）中地表开始出现形变。随着开采面的推进，地表形变范围逐渐扩大，直至监测期结束，最终形变主要覆盖的工作面为30112 工作面、30114 工作面和30116 工作面。其中，30114 工作面地表形变量级及范围都较大，最大累积形变量约为-130 mm，为整个工作区的可探测到的最大形变量，在30114 工作面上（图4）画一条剖线AA′以及AA′剖线上的形变最值特征点a 进行时间序列分析。

图5 区域①时间序列形变图Fig.5 Time series deformation of region ①

图6 为剖线AA′的时间序列剖线图，其中主要选择了6 个时间段显示，共形成了4 个形变漏斗。由图6 可知，靠近A 点的第一个形变漏斗的变化相对较为平稳；第二个形变漏斗在2021 年1 月后发生突变；第三个漏斗在2021 年5 月后发生突变；第四个形变漏斗在2021 年8 月发生突变，这三个突变情况与工作面的开采时间都较为吻合。其中，第四个形变漏斗的形变量远远超过其他区域的形变量。图7 为剖线AA′上的形变最值点a 的时间序列形变折线图。由图7 可知，2021 年8 月，该点发生了突变，而后速率逐渐减小。2022 年5 月前后，再次出现较小幅度的突变，可能是受到30116 工作面开采的影响。

图6 AA'时间序列剖线图Fig.6 Time series deformation of AA' profile

图7 a 点时间序列图Fig.7 Time series deformation of point a

图8 为区域②的时间序列形变图。由图8 可知，从2020 年8 月开始，30115 工作面处于开采状态。由于初期的形变较小，绝对量值在20 mm 之内，在图8（d）中才显示出来。最终形变主要覆盖30115 工作面和30117 工作面。在30117 工作面上（图4）画一条剖线BB′以及BB′剖线上的形变最值特征点b 进行时间序列分析。

图8 区域②时间序列形变图Fig.8 Time series deformation in region ②

图9 为剖线BB′的时间序列剖线图，其中，选择了6 个时间段显示，共形成多个形变漏斗。由图9 可知，靠近B 点的形变突变基本都发生在2021 年1 月之后，后面时间段的形变变化相对都较为平稳。图10 为剖线BB′上的形变最值点b 的时间序列形变折线图。由图10 可知，在2021 年3 月前后，b 点发生了突变，与工作面开采的时间基本吻合，4 月之后速率逐渐减小。

图9 BB'时间序列剖线图Fig.9 Time series deformation of BB' profile

图10 b 点时间序列图Fig.10 Time series deformation of point b

图11 为区域③的时间序列形变图。由图11 可知，2021 年此区域的30301 工作面才开始开采；2021年9 月，工作面的形变才开始显示。最终，形变主要覆盖的工作面为30301 工作面和30302 工作面。在30302 工作面上（图4）画一条剖线CC′以及CC′剖线上的形变最值特征点c 进行时间序列分析。

图11 区域③时间序列形变图Fig.11 Time series deformation in region ③

图12 为剖线CC′的时间序列剖线图，其中，选择了6 个时间段显示，只形成1 个形变漏斗，漏斗中心基本靠近剖线CC′的中间位置。由于剖线的位置靠近30301 工作面，在2021 年8 月和2022 年1 月，剖线上的形变先后发生两次突变，第一次突变是受到30301 工作面开采的影响，第二次突变是受到30302工作面开采的影响。图13 为剖线CC′上的形变最值点c 的时间序列形变折线图。由图13 可知，在2021年9 月前后及2022 年3 月前后，c 点发生了两次突变，与剖线整体的形变特征基本吻合。

图12 CC'时间序列剖线图Fig.12 Time series deformation of CC' profile

图13 c 点时间序列图Fig.13 Time series deformation of point c

随着时间的推移，矿区内主要形变区域的形变愈加明显，时序InSAR 监测结果数据对比反映了矿区内局部地区的地表变形情况。从获取的结果来看，矿区的形变区域是时刻发生变化的，这是因为矿区开采造成的形变随着时间的推移而逐渐演化，可以清楚地看到矿区形变的缓慢变化过程，矿区其他区域地表大部分处于稳定状态。

3.2 精度评定

精度评定是验证获取地表形变量精度的一种方法，一般分为内符合精度评定和外符合精度评定。内符合精度评定一般是获取的形变的中误差（标准差）、形变量分布直方图或者不同数据源之间形变值的相互比较。外符合精度评定一般是借助外部的形变数据，如GPS 数据或水准数据，与InSAR 获取的形变数据进行对比验证。

由于本次试验没有收集到外部数据，且SAR 数据源只有Sentinel-1A。因此，通过形变量值的标准差的计算，进行内符合精度评定。图14 为形变速率像元分布图。由图14 可知，85.26%的像元的形变速率的绝对值小于10 mm。图15 为工作区形变速率标准差分布图。由图15 可知，大部分区域的像元的标准差均在6 mm 范围内，在形变量级大的区域的像元的标准差相对较大，在6～15 mm 范围之间。

图14 形变速率像元分布图Fig.14 Pixel distribution of deformation rate

图15 形变速率标准差分布图Fig.15 Distribution map of standard deviation of deformation rate

3.3 光学影像解译

通过无人机正射飞行获取到区域①的高分辨率光学影像，将InSAR 形变与光学影像叠加到一起进行光学解译，如图16 所示。矿区地面形变的表现特征一般为地面塌陷造成的地裂缝和塌陷坑等。

图16 区域①光学影像与形变叠加图Fig.16 Optical image and deformation superposition in region ①

由图16 可知，区域①的形变主要为（a）（b）（c）三部分，地表形变最主要的特征表现为地裂缝展布，部分区域分布有多个塌陷坑，而形变量级较小的区域地表的光学特征不甚明显。结果显示，光学影像特征与InSAR 监测的形变结果的分布较为吻合，但InSAR 技术更适合小量级的形变监测，尤其是光学影像上无明显特征的区域。

4 结论

本文基于SBAS-InSAR 方法对杨伙盘煤矿在2020 年8 月至2022 年10 月间的地表进行了形变监测，并利用高分辨率光学影像进行了特征解译。

1）矿区范围内共分布了三处形变区域，均呈长条状分布，与工作面开采范围高度吻合。整个区域可探测的最大形变值达到-130 mm，最大年平均沉陷速率约为-52 mm/a。

2）通过时间序列形变分析，每处形变区在对应或相邻的工作面开采时发生突变，持续时间在1～2个月之间，随即开始缓慢的持续性形变。

3）通过高分辨率光学影像和InSAR 形变结果的综合分析，在形变量级较大的区域，在光学影像上均展布多条地裂缝，部分区域还分布有塌陷坑，表明InSAR 监测结果的可靠性。而形变量级较小的区域，光学影像上无明显的特征，表明InSAR 技术在地表小量级形变监测方面的优势。