0.   引言

我国进入新发展阶段后，对煤炭资源仍有较大需求，采空地面塌陷问题随着煤矿的开采逐渐严重^[1]，地面裂缝等地质灾害也随之发生，严重影响着矿区的安全生产^[2]。若在采空地面塌陷前期实施有效及时的措施进行监测和治理，不仅能预防大型事故，也能保护矿区生态环境的可持续发展^[3]。

水准测量、全站仪测量和GPS测量等作为采空地面塌陷传统的监测方式，由于受天气变化影响较大、对人力物力消耗较多、对监测环境要求较高及监测点不易维护等缺陷^{[4 − 8]}，难以进行长时间、大跨度的监测^{[9 − 11]}。另外，上述监测方式的结果都是基于离散点的稀疏数据^[12]，难以精确地反映区域整体的采空地面塌陷情况^[13]。

近年来，合成孔径雷达干涉测量技术（ interferometric synthetic aperture radar，InSAR ）迅速发展^{[14 − 15]}，在采空地面塌陷监测等方面具备低成本、高精度、大范围和全天候等优势^{[16 − 19]}，相较于传统的监测方式，该技术明显提升了监测工作效率^{[20 − 22]}。差分干涉测量技术（differential InSAR, D-InSAR）是传统InSAR技术的延伸^[23]，该技术可基于2景SAR影像获取地表形变的空间分布特征，但是监测结果会受到各种大气等误差的影响，精度有限。为了突破D-InSAR技术的限制，采用连续多幅影像进行时间序列分析处理的永久散射体（persistent scatterer InSAR, PS-InSAR）技术和短基线集（SBAS-InSAR）技术^[24]，二者均有监测精度高、抗相位误差干扰和揭示监测目标时序形变规律等优点。但是，PS-InSAR技术需要选取散射特性稳定、回波信号较强的PS点^[25]，且不适合大范围地区高分辨率的沉降监测，因而在城市重点设施的形变监测方面有着较多应用^[26]。SBAS-InSAR技术则更适合在缺少PS点的矿区进行监测^{[27 − 28]}，近几年SBAS-InSAR被广泛应用在采空地面塌陷信息提取和动态特征分析等方面^[29]，其监测结果的可靠性也逐渐被验证，何清等^[30]将SBAS-InSAR结果与同时期水准点测量值对比，结果显示，两者误差值在1 mm以内，其中最精确的监测值与实际值的误差仅为5.5%，这表明SBAS-InSAR技术监测采空地面塌陷的结果有较强可信度。

高家梁矿区虽仅有十几年开采历史，但是采空地面塌陷问题逐渐严重，地表草地、林地、耕地及房屋建筑受到影响。同时，高家梁矿区地表多为山地丘陵，地面高低起伏、沟谷纵横，在丘陵山地中难以使用传统的采空地面塌陷监测方式，因此矿区采空地面塌陷分布位置以及形变量等具体情况的调查较少，这对矿区安全生产和环境保护存在潜在隐患。因此，本文基于哨兵一号（Sentinel-1）卫星雷达影像数据，通过SBAS-InSAR技术，对高家梁矿区开展采空地面塌陷监测，并对20314、20313和20312这3个工作面的采空地面塌陷动态特征进行研究，为高家梁矿区的地质灾害防治工作提供数据、方法支持和措施建议。

1.   研究区概况

高家梁矿区处于内蒙古自治区鄂尔多斯市东北部，矿区范围为一近似规则的长方形（图1），东西宽约4.56 km，南北长约10.55 km，面积47.01 km²，开采标高1330～1040 m，开采方式为井工开采，采用长壁后退式采煤方法，综合机械化采煤工艺，全部垮落法顶板管理。工作面采用后退式回采，开采面积约为47.01 km²。

图  1  高家梁矿区位置

Figure  1.  Location of Gaojialiang Mining Area

下载: 全尺寸图片 幻灯片

矿区气候特征属于冬寒夏凉、日照雨少的干旱半干旱的温带高原大陆性气候。土地利用类型以耕地、林地和草地为主，分别占矿区总面积的4.28%、14.52%和74.27%，其余土地利用类型为房屋建筑。井田处于东胜煤田的东缘，区域地层属于华北地层区鄂尔多斯分区，地层倾角小于5°，区域构造总体为一向南西倾斜的近水平产状的单斜构造，由于新生代地质应力作用强烈，因此研究区地貌属构造剥蚀丘陵区，地形起伏较大，沟谷发育，呈树枝状分布。

矿区煤层顶底板岩性主要为砂质泥岩、粉砂岩和细粒砂岩等软质岩，地质构造为泥质结构和碎裂结构。经过十几年的开采，采空区引发了一定程度的地面塌陷伴生地裂缝，无明显塌陷坑，引发地面塌陷及沉陷地质灾害面积6.73 km²。裂缝宽度10～30 cm，局部可达50 cm，裂缝间距约15～40 m，裂缝可见深度0.50～2.00 m，长度100～200 m，局部塌陷裂缝较为严重，呈离层错动台阶状，台阶落差一般为0.5～1.5 m。矿区目前形成36处采空区，面积共计12.41 km²，其中就包含了位于矿区东部的20314、20313和20312工作面。

本文选高家梁矿区203盘区的3个工作面为研究对象（图1），开采顺序依次为20314、20313、20312，走向为南北走向。各工作面的开、停采时间、开采深度、开采厚度、采深采厚比和工作面长、宽等详细情况见表1。由于获取的研究区SAR数据最早时间为2018年4月，为了较为完整监测3个工作面的沉降情况，本文选取2018年4月到2020年12月作为监测时间段。

表  1  工作面详细情况

Table  1.  Details of working faces

工作面	长/m	宽/m	开采时间	停采时间	开采深度/m	煤层厚度/m	采深采厚比
20314	1600	280	2017年	2018年初	171.30～186.76	3.98～3.99	42.9～46.9
20313	2000	290	2018年初	2019年初	178.69～189.98	4.20～4.24	42.5～44.8
20312	2600	290	2019年初	2019年末	140.45～205.07	3.78～4.24	33.1～54.3

下载: 导出CSV 

| 显示表格

2.   数据和方法

2.1   数据

2.1.1   SAR数据

选取2018年4月—2020年12月期间12景同轨多时相的C波段的Sentinel-1升轨宽幅单视复数（single look complex, SLC）雷达影像为数据源，极化方式为VV极化，所选SAR影像区域位置图见图2（a）。本次研究所用的SAR影像的编号、轨道号、日期、成像模式、极化方式和飞行方式等详细信息见表2。

图  2  研究区SAR影像和DEM

Figure  2.  SAR image and DEM of the study area

下载: 全尺寸图片 幻灯片

表  2  SAR数据参数表

Table  2.  Parameters of SAR data

编号	轨道号	日期	成像模式	极化方式	飞行方式	入射角/（°）
1	026158	2018-04-12	IW	VV	升轨	42.02
2	026683	2018-07-05	IW	VV	升轨	42.02
3	027208	2018-10-09	IW	VV	升轨	42.02
4	027558	2018-12-08	IW	VV	升轨	42.02
5	028083	2019-03-02	IW	VV	升轨	42.02
6	028433	2019-06-06	IW	VV	升轨	42.02
7	028958	2019-09-10	IW	VV	升轨	42.02
8	029308	2019-11-21	IW	VV	升轨	42.02
9	029833	2020-02-13	IW	VV	升轨	42.02
10	030008	2020-04-25	IW	VV	升轨	42.02
11	030183	2020-09-04	IW	VV	升轨	42.02
12	030708	2020-11-15	IW	VV	升轨	42.02

下载: 导出CSV 

| 显示表格

2.1.2   DEM数据

所用外部数字高程模型（digital elevation model, DEM）是美国国家航空航天局（NASA）向公众发放的30 m空间分辨率的SRTM DEM数据，见图2（b），在SBAS-InSAR处理中用于去除平地效应和地形相位。

2.2   研究方法

2.2.1   SBAS-InSAR技术原理

SBAS-InSAR技术原理可简述为选取等时间间隔的N+1景SAR影像，依据设定的时间与空间基线参数，自动选取最符合条件的1景影像作为公共主影像，剩下的N景作为辅影像与公共主影像配准。得到若干个短基线组合，然后差分干涉处理每个组合内的像对，假设$ {t_{\text c}} $为公共主影像对应时刻，$ {t_{\text a}} $和$ {t_{\text b}} $时刻相对于$ {t_{\text c}} $时刻所获取的第$ i $幅差分干涉图的任意像元$ (x,y) $的差分干涉相位可表示为：

式中：$ \varphi ({t}_{\text b},x,y) $——$ {t_{\text b}} $时刻对应的相位值；

$ \varphi ({t}_{\text a},x,y) $——$ {t_{\text a}} $时刻对应的相位值；

$ \lambda $——雷达波长；

$ \mathrm{d}(t_{\text b},x,y)\overline{\overline{ }} $——$ {t_{\text b}} $时刻相对于参考影像在雷达视线 方向形变量；

$ \mathrm{d}(t_{\text a},x,y) $——$ {t_{\text a}} $时刻相对于参考影像在雷达视线 方向的形变量；

$ \Delta \varphi _i^{\text{topo}} $、$ \Delta \varphi _i^{\text{atm}} $、$ \Delta \varphi _i^{\text{noi}} $——地形残差相位、大气延迟 相位和噪声相位。

后续通过建立多项式去除地形残差，再经过相位解缠估算形变速度，进而对速度进行时间域上积分得到形变量。

2.2.2   SBAS-InSAR关键步骤

（1）生成连接图：设置空间基线阈值和时间基线阈值分别为临界基线长度的45%和72 d，共生成若干组干涉像对，其中配准时的公共主影像成像时间为2019年9月10日。图3为干涉像对的时空基线具体连接方式，黄色节点和绿色节点分别表示公共主影像和其他时间的影像，连接节点的线段为一个干涉像对。

图  3  时间-位置和时间-基线连接图

Figure  3.  Time-position and time-baseline plot

下载: 全尺寸图片 幻灯片

（2）干涉工作流：差分干涉处理上述（1）步骤中生成的干涉对，包括干涉图的去地形、去平、滤波和解缠等具体操作。经过不断调试，选取金斯坦滤波去除时空失相干噪声，解缠算法选择最小费用流法。

（3）轨道精炼和重去平：选取形变区远处的形变量较小的稳定的地面控制点（ground control points, GCP）用于重去平，利用卫星精密轨道（precise orbit ephemerides, POD）数据来进行轨道精炼，进行去除残余地形相位等操作。

（4）反演：用多项式模型估算残余地形和地面形变速度，进而去除大气相位，对相干性在0.2以上的像元执行相位解缠，反演出平均形变速度和累计形变量。

（5）地理编码：将SBAS-InSAR结果由SAR坐标系统转换为地理坐标系统。

2.2.3   技术路线

通过SBAS-InSAR技术处理Sentinel-1雷达影像，获取研究区采空地面塌陷时间序列以及平均位移速率数据，并对采空地面塌陷动态特征分析，具体的技术路线见图4。

图  4  技术路线图

Figure  4.  Technology roadmap

下载: 全尺寸图片 幻灯片

3.   结果

3.1   采空地面塌陷平均位移速度

经SBAS-InSAR处理，获取3个工作面的平均位移速度（图5），其中由暖色到冷色分别代表了沉降速度由慢到快，白色区域是由于地形因素导致的失相干地区。SBAS-InSAR技术所获取的地面形变是雷达视线（LOS）方向，需据SAR成像的几何参数将视线方向转换成垂直方向：设LOS入射角为θ，则视线方向的沉降值是$ \Delta \mathrm{L}\mathrm{O}\mathrm{S} $，垂直向的沉降值$ \Delta\mathit{\mathrm{\mathit{h}}}=\Delta\mathrm{L}\mathrm{O}\mathrm{S}/\mathrm{cos}\mathrm{\theta} $。

图  5  研究区平均位移速度及分布直方图

Figure  5.  Average displacement rate and distribution histogram of the study area

下载: 全尺寸图片 幻灯片

整体来看，203盘区平均位移速度较快的区域位于盘区的中北部，最大速度约为−17.2 mm/a。每个工作面的平均位移速度特征各不相同，20314工作面北部平均位移速度较快，大多为−10～−15 mm/a，而南半部分沉降速度较慢，为−1～−15 mm/a；20313工作面是整体沉降速度最大的工作面，位移速度的最大值区域也处于该工作面；20312工作面也呈现出北部快，南部慢的特征，位移速度较大。

为了进一步分析203盘区的平均位移速度特征，统计得到速度分布直方图（图5），由图5可知，速度在(−10，0] mm/a的区域占塌陷总面积的71%，速度大于−10 mm/a的区域占比为29%。

3.2   采空地面塌陷时序形变量

通过SBAS-InSAR技术，获取了监测时间段内各时期的采空地面塌陷相对累计形变量（图6），研究区采空地面塌陷范围和形变量随时间不断扩张和增大，塌陷情况整体呈现出“由南向北，由西向东”的特征，与工作面实际开采方向和顺序一致。

图  6  研究区采空塌陷时序图

Figure  6.  The collapse sequence in the study area

下载: 全尺寸图片 幻灯片

截至2020年末，203盘区的3个工作面相对的累计形变量如图7所示，采空地面塌陷都比较严重，主要塌陷区分布在20314和20313工作面的北部，最大形变量达到了−106 mm左右。

图  7  研究区累计形变量及分布直方图

Figure  7.  Cumulative subsidence and distribution histogram of the study area

下载: 全尺寸图片 幻灯片

统计得到形变量分布直方图（图7），呈正态分布，研究区绝大部分区域形变量处于−20 ~ −90 mm，占研究区沉降区域的90.73%；形变量处于−90 ~ −106 mm的沉降较大区域，占比5.23%；形变量处于0 ~ −20 mm的沉降较小区域最少，仅占比4.04%。

为进一步分析3个工作面的采空地面塌陷特征，分别沿着平行于3个工作面的开采方向和垂直于开采方向分别选取4条剖面观测线，其中南北向的观测线分别命名为Z1、Z2和Z3，东西向的观测线命名为Q1（图8），分别绘制各测线累计形变量曲线。

图  8  测线布置图

Figure  8.  Survey line layout map

下载: 全尺寸图片 幻灯片

由南北向剖面线累计沉降曲线（图9）可知，随着时间的推移，20312和20314工作面剖面线的形变量由南向北不断增大，20313工作面剖面线的形变量由北向南不断增大，与实际开采状况一致。20314工作面已于2018年4月之前停采，因此该工作面于观测时间段的沉降情况主要受到临近的工作面开采活动的影响。2018年4月12日至2018年12月8日是20313工作面的开采时间段，从沉降曲线来看，曲线间隔明显增大，反映了工作面沉降剧烈，沉降中心由北向南移动，形成了明显的漏斗形状的沉降盆地，这说明此时20313工作面的开采活动已经影响到了上覆地面，该工作面的开采沉降也影响到了临近工作面，在20314和20312工作面相应时段的沉降曲线间隔也都明显变大。2019年3月2日至2019年11月21日是20312工作面的开采时间段，沉降中心由南向北移动，最终在工作面北半部分形成了明显的漏斗形状的沉降盆地，20314和20313工作面沉降曲线上的间隔增大也正是因为受到20312工作面开采沉降的影响。在2019年11月21日之后，附近再无开采活动的影响，203盘区处于沉降缓慢阶段，直至2020年11月27日，工作面沉降趋于稳定，20314、20313和20312工作面剖面线最大沉降累计量分别为−32 ，−70 ，−76 mm，分别位于工作面的北边界和由北向南三分之一处。

图  9  20314、20313和20312工作面南北向测线（Z1、Z2和Z3）累计沉降曲线图

Figure  9.  Cumulative Subsidence Curve Chart of North-South Survey Lines (Z1, Z2, Z3) for Working Faces 20314, 20313, and 20312

下载: 全尺寸图片 幻灯片

由东西向测线（Q1）累计沉降曲线（图10）可知，东西方向上，整体时序曲线成“V”字型，呈现出沉降盆地的明显特征，图中两次沉降曲线间隔较大的剧烈沉降期发生在2018年4月12日至2019年3月2日和2019年6月6日至2019年11月21日分别与20313和20312工作面开采停采时间吻合。其中，20314工作面下沉量很小，下沉不明显；20313工作面下沉量最大，明显大于两侧，处于沉降中心，最大值达到了约−60 mm；20312工作面下沉量较小，处于塌陷盆地的边缘位置，下沉曲线陡峭，下沉量变化快。

图  10  东西向测线（Q1）累计沉降曲线图

Figure  10.  Cumulative subsidence curve chart of East-west survey line (Q1)

下载: 全尺寸图片 幻灯片

3.3   采空区地面验证

为了对SBAS-InSAR结果进行验证，前往高家梁矿区203盘区的沉降情况进行实地探查。从上述累计沉降曲线图可知，采空区的沉降呈现为漏斗状，实地探查过后发现，大部分区域地表沉降肉眼无法分辨，但是在沉降漏斗边缘的形变量变化较大的区域发现了明显的沉降裂缝。本次探查到的典型地表裂缝位于20314工作面北边界附近，见图11（a），与上述剖面沉降分析中20314工作面的最大沉降位置吻合。从现场拍摄的沉降裂缝照片来看，见图11（b），裂缝规模较大，沉降区地面标高下降明显，对地表的草地造成了严重破坏，对附近房屋及地面人员存在着安全隐患。

图  11  （a）研究区典型裂缝位置；（b）典型裂缝照片

Figure  11.  (a) Location of typical cracks in the study area; (b) photo of typical cracks

下载: 全尺寸图片 幻灯片

4.   讨论

总体来说，监测期内高家梁矿区203盘区整体地质情况相对稳定，采空区地面塌陷以缓慢沉降为主，局部地区有较快的沉降，煤矿开采是造成采空地面塌陷的主导因素，井工开采形成的采空区，使得煤层上覆地层岩体受重力作用发生弯曲、断裂和垮落，进而影响到地表，形成采空地面塌陷。然而，203盘区呈现出采空地面塌陷不均匀，这是因为煤层上覆地层岩性是强度低、易破碎的砂质泥岩和粉砂岩等软质岩，随着顶底板的掉落，围岩发生指向采空区的变形，最终在地表形成裂缝和塌陷坑等不连续变形。

每个工作面的采空地面塌陷特征都不尽相同。20314工作面作为最早开采的工作面，开采时间为2017年，而监测时间是2018年4月开始，这使得20314工作面从开采到2018年4月这段时间的采空地面塌陷情况未能监测，由于20314工作面是由南向北开采，因此这段时间的地面塌陷主要发生在工作面南部，这也解释了20314工作面的平均位移速度为何呈现“北部快，南部慢”的特征，以及剖面形变量曲线呈现出沉降漏斗不明显且北部形变量大的特征。20313工作面位于其余两工作面之间，会受到其余两个工作面采空地面塌陷的叠加影响，因此20313工作面既是平均位移速度最快的工作面，又是累计形变量最大的工作面。20312工作面的采空地面塌陷也被完整监测，因此平均位移速度较快且形变量较大，与203盘区整体采空地面塌陷形变量特征一致，呈现出“北部大，南部小”的特征。

203盘区的塌陷中心位于距离北部边界约900 m、距离东部边界约300 m处，沉降漏斗形成的盆地呈现近似椭圆形，范围是距203盘区北边界200~1500 m、距东边界100 ~600 m，实地验证发现的典型地裂缝就是在203盘区北部边界附近，可以预计未来地裂缝等地质灾害的发生也很有可能在这范围内，矿区应对此范围内的地质灾害进行细致排查和治理，同时加强监测，谨防地质灾害的发生。

5.   结论

（1）203盘区平均位移速度呈现出“北部快，南部慢”的特征，速率在(0，10]（mm/a）的区域占盘区总面积的71%，速率在10 mm/a以上的区域占盘区总面积的29%，研究区在监测时段的采空地面塌陷以缓慢沉降为主，局部地区有较快的沉降，最大沉降速率约为17.2 mm/a，位于20313工作面北部三分之一区域。

（2）研究区沉降范围和形变量随时间不断扩张和增大，沉降情况整体呈现出“由南向北，由西向东”的特征，与实际的工作面开采方向和顺序一致，主要沉降区分布在20314.20313盘区北部，最大形变量达到了106 mm左右，应当作为矿区后续灾害治理的重点关注区域，谨防地裂缝或者地面塌陷等地质灾害的产生。

（3）在对实验结果进行实地验证后，发现203盘区产生裂缝的实际位置与实验结果基本吻合，表明了SBAS-InSAR技术也可以良好地适用于地表高低起伏、沟谷纵横的丘陵地貌的矿区，且具备较强的技术优势，可以给同类型矿区采空地面塌陷监测提供技术支持。