李小涛，张孝荣，黄景偲，汤伏全

（1.陕西一八五煤田地质有限公司，陕西 榆林 719000； 2.西安科技大学测绘科学与技术学院，陕西 西安 710000）

D-InSAR作为一种新型的地面观测技术，已逐渐成为雷达遥感技术的热点研究方向[1]，通过InSAR技术可以获得cm甚至mm的高精度表面变形信息[2-3]，因此该技术在地震、冰川漂移、火山喷发、城区地面塌陷、山体滑坡、泥石流等灾害监测中发挥着重要作 用[4]，但在矿区大范围开采沉陷监测中应用相对较少。本文选取2018-01-06～2019-01-13的31期Sentinel-1A雷达影像，基于小基线集（SBAS-InSAR）技术对榆神矿区地表沉降进行动态监测，基于变形序列结果分析地下开采工作面与地表变形的量化关系，实现了榆神矿区大范围的遥感沉降监测。

1 关键技术

1.1 SBAS-InSAR原理

小基线集（SBAS-InSAR）技术由Berardino于2002年最先提出[5]，该方法的基本思想是为了有效抑制时空失相干及大气延迟带来的误差，将获取的 SAR 影像按照时间基线及空间基线阈值分成若干小基线集组合，对小基线组合内的干涉图进行差分干涉处理，并利用奇异值分解法（SVD）计算出研究区的时序形变量。

设有按时间序列t0,t1,…,tn获取N+1幅单视复数影像，首先将他每年以任意影像为主影像进行配准，然后设定垂直基线阈值，将垂直基线小于该阈值的SAR影像归为一组，共分为L组。对每组内的影像进行差分干涉处理，最终L组影像共得到M幅差分干涉图，假设N为奇数，则差分干涉图的个数M可以表示为：

以t0时刻为初始时刻，对于任意时刻ti(i=1,…,N)相对于初始时刻的差分相位φ(ti)为未知参数，观测量为数据处理获取的差分干涉相位δφ(tk)(k=1,…,m)如果所有的差分干涉图都正确解缠，而且差分干涉相位被校正到某个稳定或形变信息已知的高相干像元(x,r)上，该像元则为参考像元。参考像元解缠后的差分干涉相位为[6]：

式中，φ(tb,x,r)和φ(ta,x,r)分别为d(tb,x,r)和d(ta,x,r)对应的解缠相位值；φn为随机噪声；λ为波长；d(tb,x,r)和d(ta,x,r)分别为tb和ta时间内得累计LOS向形变量；v(x)为像元x处的形变速率；Ti为主辅图像的时间间隔。

对任意干涉组合均能表示上式的线性方程组，矩阵表示如下式：

式中，B为系数矩阵，每行对应一幅干涉图，每列对应一个时间上的SAR图像，影像分别为1和-1，其余列为0。如果m≥n，且B的秩是N，为满秩矩阵，则最小二乘解得形变速率v的估计值：

为了保证干涉图的高相干性，通常选取时间基线和空间基线都比较小的干涉组合，使得原有的SAR影像被分成L个独立的干涉组合，由于系数矩阵B的秩 L＜N，导致上式秩亏，方程组的解不唯一，因此需要对B进行奇异值分解。

式中，U是m×m阶酉矩阵；∑是半正定m×n阶对角矩阵，其对角线上的元素∑i为B的奇异值；而VT，即V的共轭转置，是n×n阶酉矩阵。因此形变速率的估值v为：

将形变速率积分即可获得每个SAR影像的累计地表形变相位。

1.2 下沉量提取

由于InSAR技术只能监测雷达卫星视线向的一维位移量[7-8]，如要进一步分析井下工作面开采与地面下沉范围之间的关系，需要获取地面真实的三维位移量。以往许多学者针对确定矿区下沉移动盆地边界的问题时，一般采用忽略水平位移的叠加影响，直接将LOS形变量投影到下沉方向作为下沉量的方法。

本文基于开采沉陷的对称原理（如图1所示），相对于各工作面对应地表下沉盆地中心对称的两点P1和P3之间下沉量相等，水平位移大小相等方向相反的特征，对P1和P3点进行相加就能抵消掉水平位移，从而分别提取出该点的真实下沉量，如公式（8）[9]。

图1 工作面地表对称点的水平位移和下沉示意图

本文在获取地表下沉量时，均根据公式（8）求出采矿区任意地表点的下沉量。该下沉量剔除了水平位移的影响，能够获取真实的三维形变量。

2 研究区监测实验

2.1 研究区概况

榆神矿区跨陕西省榆林市榆阳区和神木市，是国家14个大型煤炭基地中陕北基地的主力矿区之一，总面积约5 265 km2，已探明储量733.5亿吨。矿区位于鄂尔多斯高原东南部，毛乌素沙漠的东南缘与陕北黄土高原北部的接壤地带，总地势西北高、东南低。榆神矿区西北部以陕蒙边界为界，西南部以榆神矿区北边界为界，东北部以神府矿区西边界为界。矿区地理坐标为109 08′24″E～110 27′59″E，38 19′33″N～39 11′23″N。

2.2 实验数据及数据处理

本文的实验数据选取2018-01-06～2019-01-01的31期Sentinel-1A雷达影像，其工作模式为条带模式，幅宽为250 km，距离向分辨率为5 m，方位向分辨率为20 m，时间分辨率约为12 d[10-11]。较高的时间、空间分辨率保证了影像对的相干性和变形监测精度。

DEM数据选用美国地质调查局网站公布的SRTM DEM数据，并且选择覆盖范围为（38 N～39 N，109 E～110 E）的DEM。

本实验基于小基线原理，选取时间基线不超过48 d，垂直基线不超过200 m，30幅影像得到了78个干涉良好的干涉像对。对选取的78个干涉像对分别进行差分干涉，去除地形相位，滤波处理；再利用干涉处理得到的干涉图，依据相干系数法，对高相干点进行相位解缠；然后利用奇异值分解（SVD）算法， 计算出高相干点的线性形变速率及高程残差，通过时空滤波及形变速率的积分累积运算，得到时序累积形变结果。

3 金鸡滩矿地表动态沉陷与地下开采关系分析

本论文研究基于GAMMA软件进行SBAS技术处理，绘制金鸡滩106工作面从2018-01-06～2019-01-01期间的地表时序形变图2，该图中的黑色虚线框为井下工作面开采面。

从图2可以看出2018-03-31之前地表并未出现地表形变，到2018-05-06，地面开始出现了形变，且往后形变范围持续向东北方向延长，呈长条状。根据动态变形情况可推断106工作面的开采影响时间是2018-04，经实地调查矿区的开采资料，该工作面确实于2018年四月中旬开始开采，与监测到的地表形变情况相符。

图2 金鸡滩106工作面地表的时序形变图

为了揭示监测到的地表动态沉降特征与煤矿地下开采工作面推进情况的关系，根据实际采矿资料，该变形区不同时段对应的井下工作面推进距离增量与变形区域前沿位置距工作面开切眼的水平距离之间的相对关系如表1所示。

表1 变形区前沿扩展与工作面实际推进距离的关系

从表1可见，地表变形区域的前沿位置扩展的距离与相应时间段地下工作面的推进距离之比值始终在1∶1左右，其动态特征符合开采沉陷的基本规律，也与该观测站的实测数据一致。说明随着地下工作面的不断推进，地表变形区向前扩展的速度与地下工作面推进的速度基本相同，变形区域随着地下工作面的推进而同步发展。

4 榆神矿区地表动态沉陷监测结果分析

基于以上对金鸡滩106工作面地表形变与地下工作面开采关系的研究基础，可以将SBAS-InSAR技术和本文的下沉量提取方法应用到大范围变形区域的时序形变分析。首先可利用干涉图堆叠技术搜索出整个榆神矿区的形变区域分布，如图3 为整个榆神矿区在2018年3月到2019年1月间的形变累积速率图。然后根据就搜索出的形变区分块分别进行时序InSAR数据处理，并提取出地表下沉量。

图3 榆神矿区2018年InSAR监测变形速率分布

从累积形变结果图中总共监测到20个变形速率大于0.1m/year的变形区域，分别位于杭来湾、榆树湾、曹家滩、小保当、柳巷、千树塔等煤矿开采区，其中 9号形变区为前文的金鸡滩矿区的106工作面。

为了统计出现的变形区范围的动态变化，利用SBAS-InSAR技术继续处理31期SAR影像，将榆神矿区地表变形区域分成9个区块，分别统计9个变形区在2018年度不同时段累积的变形区面积统计如表2所示。

表2 各变形区累计变形面积统计/km2

从表2可以看出，位于金鸡滩矿区的2号、9号形变区，杭来湾的3号形变区、榆树湾的5号6号形变区的形变范围均发生了明显的变化，结合图3中各个形变区的形态及表2中形变区面积变化，可进一步推断地下工作面的开采动态情况。

其中，1号和2号形变区位于金鸡滩煤矿区域，根据其变形区的西南端窄，东北端宽的形态可推断该地区开采前进方向是东偏北方向，并且从2018年1月至2019年1月的形变范围基本无变化，推断该地下工作面于2018年初停止开采。2号变形区呈长条状，2018年4月到10月间形变范围逐渐增大，根据其西南端宽而东北端窄，且10月9日到11月2日间的变形区域基本保持不变的表现特征，可以推断该地下工作面沿西偏南方向进行，并且于10月份的地下开采工作出现了暂停。结合采矿资料与上述对地下工作面开采情况的推断完全符合。9号在该监测时间范围内其形变面积逐渐增大。

3号形变区位于杭来湾煤矿，该形变区呈相邻两个长条状的形态，分成左上角和右下角两块，从表2中的数据显示3号变形区面积也在逐渐变大，并且从形变图3上可以看出该形变区中心出现了明显的失相干现象，说明在该时间段内地下正在处于开采工作状态，并且该地下工作面的开采造成了其对应地表出现了大梯度变形。第4变形区域呈圆状，且一年内形变范围没有变化，可推断该区域地下处于停采状态且存在残余变形。

5、6号变形区域位于榆树湾煤矿，其中5号变形区域面积基本没有变化，推断在2018年初地下开采工作逐渐停止，地表处于残余变形。6号变形区在2018年5月前形变面积逐渐增大，但在5月以后变形区范围不再变化，推断该区地下工作面开采在5月后停止，地表属于残余变形。经实地调查该矿6工作面在5月停采，与监测到的形变情况相符。

7、8号变形区域位于柳巷矿区，其中7号变形区面积和位移量基本没有变化，可推断该时间段内地下开采工作已停止，地表处于残余变形阶段。8号区域于2018年5月期间地表出现变形，呈长条形状，7月以后变形区范围没有随时间扩展，推断该区域于5月期间开始地下开采工作，然而在七月后地下停止开采，属于小面积的地下开采。

在沉陷盆地中央下沉梯度较大的情况下，由于SAR影像的波长有限，导致时序累积D-InSAR监测结果不能真实地反映沉陷区中央的绝对变形量。但是利用其微小形变结果能够较为准确地反演地下的开采时间以及工作面的开采位置及其范围情况。

5 结 语

1）基于时序D-InSAR的煤矿区地表沉降监测相比于传统测量而言，具有全天时、全天候、覆盖广、监测结果更直观等诸多优势。可以从时间和空间上反映出井下工作面开采后地表的变形范围和动态发展趋势，能够实现矿区大范围、长时序的地表变形区域探测，该技术在西部煤矿大规模开采地表沉陷遥感监测中具有普适性。

2）D-InSAR技术对于煤矿地表沉陷区边缘的小梯度变形监测非常有效，通过地表变形信息可揭示动态变形边界与井下工作面推进边界之间的量化关系。实验表明，2018年1月至2019年1月榆神矿区地表动态变形区域的范围和分布特征与地下工作面开采推进范围和方向呈显著的线性相关。