某矿山地表移动及岩层监测稳定性研究

王少特， 张耀平， 刘书瑶， 雷大星， 许传金

（1. 江西理工大学资源与环境工程学院，江西 赣州 341000； 2. 赣南科技学院，江西 赣州 341000）

随着矿产资源开采规模不断扩大，采空区所引起的地质灾害随之增加。采空区的形成破坏了原岩应力原有的平衡状态，产生了二次应力场，对地下矿体的回采以及对采空区的稳定性造成了严重的影响[1]，并且威胁到井下工人和设备的安全，因此，矿山迫切需要对矿区内可能存在的问题予以防治，在采空区的勘测工作中应有效掌握开采环境内的地理位置、形态范围及其现状，以便科学判断可能因采空区的存在而导致的工程地质灾害。本文选取安徽某矿山为研究对象，其特点是空区遗留时间过长，围岩应力不稳定，需要及时支护及充填，但空区治理方案没有得到统一。国内外大量学者做了相关的研究，在地表沉降方面重点讨论了InSAR 技术进行地表沉降监测的可能性，刘毅[2]从多个角度概括了中国地表沉降研究取得的新进展，并结合InSAR 技术探讨了引起地表沉降的多个原因，以及我国沉降监测研究方面所面临的挑战。以上两者都取得了丰富成果。但针对特殊空区的研究目前成果较少，何旭等[3]使用SBAS-InSAR 技术利用主成分分析法对北京平原的地面沉降时空特征进行了分析，采用最小二乘法和奇异值分解获取地面形变数据，上述研究表明InSAR 技术和地表监测结合是有效的监测手段，能为矿区变形监测与预计提供全新的手段[4]，但截止目前，针对特殊矿山原岩应力还没有得到深入的分析，因此，本文开展矿山地表移动及岩层监测稳定性研究，对后续矿山开采充填以及进一步的深层开采具有一定的指导意义[5]。

1 基于InSAR地表变形监测

1.1 InSAR形变测量原理

InSAR 技术是对形变前后的2 组干涉对进行差分处理，通过从干涉相位信息中去除因地形和参考面因素造成的相位变化[6]，分离出地表的形变信息。任何形式的干涉测量技术（PS、SBAS等）均是以常规“二轨”法InSAR 技术作为基础，因此以下合成孔径雷达差分干涉测量的基本原理及处理流程均以“二轨”法InSAR为例进行说明[7]。其原理如图1所示。

在重复轨道干涉测量情况下，InSAR 干涉测量相位与所选的参考椭球、地面点所处的位置、地表是否发生形变及形变大小、大气延迟误差和SAR 系统噪声等因素有关[8]。干涉相位可表示如式（1）所示：

式（1）中：φflat为参考椭球面相位；φtopo为地形相位；φdefo为沿视线向（LOS）的形变相位；φorbit为轨道误差引起的相位；φatm为大气延迟相位；φnoise为噪声引起的相位[9]。

其中，轨道误差φorbit可以用一个低阶的多项式进行拟合或利用2次观测影像的精密轨道数据消除；大气延迟相位φatm与噪声引起的相位φnoise可以直接忽略或通过滤波方法在一定程度上消除[10]。

1.2 InSAR监测误差源分析

InSAR 测量技术是利用微波合成孔径雷达图像（SAR）数据，对地表重复观测所形成的微波（1 mm～1 m）相位差计算地表形变，研究InSAR 差分干涉测量的误差来源，可以在数据处理过程中进行有效的控制，以提高数据处理结果的准确性[11]。其误差主要由以下几个方面组成：InSAR 去相干源和相位噪声；基线误差或轨道误差；DEM 本身引入的误差；大气效应。

因此，必须精确地球表面和卫星间的距离、获取影像时的地理坐标以及传感器的位置[12]。进行“二轨法”差分干涉测量时，选用的DEM固定应尽量选取垂直基线较小的影像，一般要求小于100 m[13]。

本文所用的 SAR 影像数据为Sentinel-1 数据，该数据由载有C 波段合成孔径雷达的两颗卫星组成[14]。统计数据如表1 所列，自2020 年1 月12 日至2021 年6 月11 日地表分辨率保持不变，说明在这段时间内数据可靠。

表1 Sentinel-1 数据参数Table 1 Sentinel-1 data parameters

1.3 InSAR 处理及解算结果

采用SBAS-InSAR 技术对开采区域自 2020 年1 月至 2020 年 5 月影像进行解算，并根据生成的干涉图质量及相干性进行挑选，最终获取干涉对用于时间序列分析[15]，具体时间基线如图2所示。为进一步提高相干性，分别在方位向去除不重叠的多普勒频率，在距离向进行频域滤波[16]。然后，利用DEM数据将干涉图中地形相位去除，得到差分干涉，并对其进行地理编码[17]。对时间位置连接图形成影像进行解缠计算，如图3 所示。

图2 时间基线图示Fig.2 Time baseline diagram

图3 时间位置连接图示Fig.3 Time position connection diagram

地表形变如图4 所示，对矿区2020 年1 月至2020 年6 月影像进行时序动态连续解算，由图4 可知，在 2020 年1 月至6 月，地表形变变化不大，约为2 mm；没有形成大量级形变，均为无规则凌乱区域形变，不符合开采沉陷规律，通过现场调查及结合光学影像分析可知，主要原因为植被的生长、砍伐及地表土变化。

图4 2020年1月至2020年6月地表形变Fig.4 Surface deformation map from January 2020 to June 2020

根据2020 年1 月至2021 年6 月雷达影像，矿区内建筑物形变如图5所示：整个监测区域没有发生地表裂缝、塌陷及断裂等现象，说明地表发生的形变不是由井下资源开采所诱发。建筑物出现形变约为2 mm，对建筑物稳定性不会产生影响，通过现场调查建筑物未出现裂缝，截至2021 年6 月建筑物结构稳定。

图5 2020年1月至2021年6月建筑物形变Fig.5 Building deformation map from January 2020 to June 2021

1.4 InSAR 解算分析

从2020 年1 月至2021 年 6 月的解算影像图可知，地表有一定的形变，其主要原因是山坡在进行改造、河道清理、地表土的变化及植被等原因造成。其他区域及建筑物均没有发生较大形变。根据形变图选取研究区域中的5处区域进行分析，剖面分布图如图6所示。

图6 剖面分布Fig.6 Profile distribution diagram

由图4至图7 可知，在 2020 年 1 月至 6 月，地表形变变化约为2 mm，地表植被较密的区域在植被生长时发生隆起变化，当植被枯萎时发生沉陷变化，由剖面图可知隆起最大值约为18 mm，沉陷最大值约为14 mm；在裸露的地表及建筑区域基本没有发生形变，建筑物出现形变约为2 mm，对建筑物稳定性不会产生影响，通过现场调查建筑物未出现裂缝，截至2021年6月建筑物结构稳定。结合光学影像及现场调查说明局部区域发生形变主要原因是地表植被及地表土发生了变化，地表没有发生实际形变。

图7 累计形变剖面图：（a） 1处； （b） 2处； （c） 3处； （d） 4处； （e） 5处Fig.7 Cumulative deformation profile diagram：（a） 1 location； （b） 2 location； （c） 3 location； （d） 4 location； （e） 5 location

2 采空区围岩二次应力变化及内部位移监测

2.1 监测方法及监测点布置

应力计监测点布置如图8所示，为了对围岩可能存在破坏的部位及范围进行分析[18]，本文主要针对-20 m 中段围岩进行围岩应力变化监测和沉降变形监测。通过应力监测获得围岩应力变化量、变化速度、应力分布情况，为对围岩可能存在破坏的部位及其范围进行分析提供依据[19]。结合现场实际情况，本次围岩二次应力监测区域主要布置在-20 m中段[20]。

图8 -20 m中段应力计监测点布置示意Fig.8 Layout of monitoring points of stress gauge at -20 m middle section

多点位移计布置如图9所示，通过多点位移计对各测点的位移变化量和变化速度进行监测，以位移值和位移变化速度的突增作为危险情况发生的判断依据[21]。多点位移计的重点监测布设分布在-20 m。

图9 -20 m 中段多点位移计布置示意Fig.9 Layout of multi-point displacement meter at -20 m middle section

2.2 监测数据结果分析

为了与前文InSAR 地表监测形成对照，并且为后文的数值模拟部分提供可靠的参考对比[22]，分别绘制出YL1 到YL6 应力变化图以及WY1 到WY2 位移变化图，如图10—图11所示。

图10 YL1-YL6 变化曲线Fig.10 Variation curve of YL1-YL6

图11 WY1-WY2 变化曲线Fig.11 WY1-WY2 variation curve

通过3 次的监测数据对比可知，围岩二次应力的变化值最大为 YL6 的1.91 MPa，整体的应力变化约为1.5 MPa。每个监测周期较上一轮相比，围岩二次应力都有所增大，但整体围岩二次应力的绝对值不大。监测结果表明，围岩内部监测位移变化量都在毫米级别以内，最大位移量为2.5 mm。根据图10可得，曲线YL6 上升是因为采场回采以及遗留空区的充填对周边岩体应力集中造成的。整体上围岩内部位移量变化较小，基本可以认定-20 m 中段围岩内部的位移平稳，未出现大的地压活动，YL6 及WY2 变化较大，应当重点监测。对采空区周围围岩及顶板有效实时监测以及对大空区实时预警监测，有利于矿山企业提前采取应对措施，避免灾害性事故发生。

3 采空区稳定性及地表变形数值分析

3.1 采空区三维数值计算模型

由于矿区+20 m 以下采空区全部充填治理完毕，对上部采空区及地表稳定性影响不大，因此在物探勘查结果的基础上[23]，重点分析+20 m 和+60 m 两中段遗留采空区自身稳定性和对地表变形以及对地表构筑物的影响。为了构建三维地质模型，实现矿段地表变形与采空区稳定性评价的数值模拟[24]，对地层进行概化：对实测的+20 m 和+60 m 采空区空间分布及其与矿体的位置关系；对矿段-20 m、-60 m中段矿体平面概化所确定范围内的矿体，将其作为模型中的矿体模拟对象；成矿前断裂构造与成矿时断裂构造时间较为久远，主要对矿床的形成有较大影响，对目前的矿区矿岩完整性影响不大，而成矿后断裂构造对矿段矿岩体影响较小，因此，未考虑断层构造。

3.2 模型介质力学与破坏准则

数值模型初始边界条件为：前后左右4个面及底面都使用位移约束，也就是前后左右4个面限制水平方向上的位移，底面限制垂直方向上的位移，在原始应力场垂直方向施加重力场[25]。假定模型中矿岩体均为理想弹塑性连续介质，使用摩尔-库仑（Mohr-Coulomb）屈服准则如式（2）所示，其物理力学模型以弹塑性力学模型为基准[26]。

式（2）中：σ1为最大主应力，单位为MPa；σ3为最小主应力，MPa；ϕ为内摩擦角，（°）；c为介质黏结力，MPa；fs为破坏判断系数。当fs＜0 时，介质为弹性变形阶段；当fs≥0 时，介质为塑性流动阶段。

3.3 岩矿物理力学参数取值

对矿山原岩进行分析后，进行室内实验及现场实验，可以得出不同矿石、围岩及断层间具体的物理力学参数[27]，如表2 所列。将岩石物理力学参数代入建立的模型当中进行赋值并进行计算，得出计算完的数值模型，并对剖面位移、应力云图进行分析[28]。

表2 矿岩物理力学参数Table 2 Physical and mechanical parameters of ore rock

3.4 地表沉降变形

矿区地表沉降见图12，所得最大值为2.7 mm，移动带范围内沉降基本为0，空区相对较小，说明矿区地表最大沉降为-2.7 mm，最大沉降出现在遗留采空区上方位置。整个矿区的位移变化主要集中在移动带范围内，移动带范围外沉降量基本为0，在遗留采空区上部沉降相对大点，遗留采空区范围外沉降值较小为2 mm；遗留采空区对地表沉降有一定影响，但是相对影响不大，对整个地表的稳定性影响较小。通过云图可知，在临近移动带范围外的建筑物附近地表沉降基本为0，对地表构筑物的稳定性不造成影响。

图12 矿区地表沉降云图Fig.12 Cloud diagram of surface subsidence in the mining area

马路上布设的4 个监测点位移沉降变化曲线如图13 所示，随着开采进行，地表形变也慢慢形成，越靠近空区位置地表的形变越大，3#和4#监测点的采空区上部沉降位移值相对较大，3#监测点位移值最大，沉降值约为2.2 mm。说明遗留采空区对地表的形变的形成有一定的促进作用，而且根据形变监测图可知，后续开采中位移慢慢趋于定值，说明闭坑后不存在较大的工程扰动，地表会趋于稳定状态。

图13 1#、 2#、 3#、 4#监测点沉降变化曲线Fig.13 Settlement change curve of 1#， 2#， 3#， 4# monitoring points

3.5 地表水平位移

矿区地表水平位移如图14 所示，矿区地表沿x方向最大水平位移 0.2 mm，较大水平位移主要集中在移动带范围内左侧，在移动带范围内左侧有与x轴相反方向的位移，最大值为0.24 mm。矿区地表沿y方向最大水平位移为0.9 mm，基本以遗留采空区为界线，分成沿y轴方向上，上下两部分位移，上部位移方向朝下，下部位移方向朝上，位移值均较小，不足1 mm。结合水平方向上分析可得：在水平方向位移值较小，均小于1 mm，遗留采空区对地表在水平方向产生的位移很小，不足以对地表构筑物产生破坏。该部分数值模拟分析结果与2.2 节围岩监测位移变化值基本吻合，说明该部分的数值模拟结果与现场监测结果基本一致，该模拟显示出的规律能反映现场的应力变化值。

图14 x 和y方向水平位移云图Fig.14 Cloud chart of horizontal displacement in x and y directions

3.6 采空区围岩应力分析

高程+20 ～+60 m 范围内采空区中剖面围岩第一主应力分布如图15所示，模拟结果表明：2线剖面采空区围岩第一主应力范围约-25 ～-10 MPa，均为压应力（图15（a））。采空区最大第一主应力出现在采空区右侧矿柱位置处，约-25 MPa；最小第一主应力位于空区周围，约-10 MPa。采空区左右两侧部位应力相差较大，变化约为-15 ～-25 MPa，右侧矿柱应力值大于左侧，顶底板应力值基本相差不大。由图15（b）可得：0～1 线剖面+20 m 采空区，围岩第一主应力范围约-40 ～-5 MPa，均为压应力。采空区最大第一主应力出现在空区左侧角隅位置，约-40 MPa。最小第一主应力位于空区顶部位置，约-5 MPa，顶板位置相对稳定，在矿柱角隅出现应力集中，但是未达到岩体的极限抗压强度，在不受极大的工程扰动情况下，能维持自身的稳定性。0～1 线剖面+60 m 采空区，围岩第一主应力范围约-10～-5 MPa，均为压应力。采空区最大第一主应力出现在空区矿柱位置，约-10 MPa，并且整个矿柱受力均匀；最小第一主应力位于空区顶部位置，约-5 MPa，说明顶板处于相对稳定。

图15 采空区围岩第一主应力云图：（a） 2线剖面；（b） 0～1线剖面Fig.15 Cloud diagram of first principal stress of surrounding rock in goaf：（a） 2 line profile；（b） 0-1 line profile

4 结 论

1）通过地表变形监测结果可知，整个矿山的监测区域地表形变不大，均约为2 mm；移动带构筑物形变约为2 mm，建筑物没有出现裂缝等损害现象。个别区域发生零散性形变，地表未发生裂缝、断裂及塌陷等情况，形变区域不符合沉陷规律，综合分析可知矿区移动带内地表整体稳定，地表未出现实质性变形。

2）通过数值模拟对采空区进行稳定性分析可知，矿区的地表最大沉降为2.7 mm。采空区围岩第一主应力值范围主要在-40 ～-5 MPa之间，最大第一主应力值为-40 MPa，+60 m采空区受深部开采影响相对较好，空区稳定性相对较好；总体分析+20 ～+60 m采空区处于稳定状态，最大第一主应力区域大多集中于采空区矿柱角隅处，存在部分应力集中，但均未达到岩石的极限抗压强度，在不受强扰动情况下，矿柱基本都可以维持自身的稳定性。

3）通过采空区围岩二次应力变化及内部位移监测可知，围岩二次应力变化值最大为1.91 MPa，最大位移量为2.5 mm，围岩二次应力变化值较小，围岩处于相对稳定的状态，整体上围岩内部位移量变化较小，可以认定-20 m中段围岩内部的位移平稳，未出现大的地压活动。

4）矿区整体内部变化较小处于稳定状态，未出现大规模地压活动，说明地表监测、二次围岩数据分析及数值模拟均得出矿山处于稳定状态，对于后续矿山开采充填以及进一步的深层开采有指导意义。