基于InSAR-FLAC3D的采空区稳定性分析*

兰 洋，袁利伟，李素敏，冯家宁，刘大荣，王 华

(1.昆明理工大学 公共安全与应急管理学院，云南 昆明 650093；2.昆明理工大学 国土资源工程学院，云南 昆明 650093；3.昭通永善金沙铅锌矿，云南 昭通 657000)

0 引言

因采空区导致的地质灾害事故频发，采空区稳定性研究越来越受到专家学者的关注。AGIOUTANTIS等[1]运用有限元法和边界元法研究了采动影响下覆岩垮落的开采条件和垮落高度、覆岩产生离层裂缝的力学条件及离层裂缝的位置和高度；李夕兵等[2]对金属矿地下采空区的探测、处理与安全评判进行了研究；JAISWAL等[3]为解决采空区的顶板控制问题，在开采前对其结构进行了模拟，并对地层参数进行了校准。万文[4]综合运用遗传算法、强度折减法和FLAC3D软件研究了采空区对边坡稳定性的影响；刘晓明等[5]使用CMS激光对采空区进行了探测，利用获取的采空区三维边界数据，通过自编的程序接口导入FLAC3D软件进行采空区围岩稳定性分析。

德国、波兰等国家率先将InSAR技术应用于矿区开采引起的地表沉陷监测。我国在InSAR监测地表变形方面的研究起步较晚，且前期主要将其用于城市地表沉陷监测。张学东等[6]利用SBAS-InSAR技术生成差分干涉图，以监测地面沉降；独知行等[7]综合比较了传统测量与InSAR测量的优势和缺陷，提出了将传统测量手段与InSAR技术结合共同监测矿区地表沉陷的构想；陈继伟等[8]利用Sentinel-1A TOPS 模式影像对黄河三角洲地区进行了变形监测，证明了TOPS 模式 SBAS 时序分析方法应用于变形监测的可行性和有效性。

目前多采用单一方式对采空区稳定性进行评价，缺乏系统的综合监测。基于此，本文拟采用InSAR与FLAC3D协同分析采空区的稳定性。

1 工程概况

永善县金沙厂207铅锌矿体(简称“207矿体”)处于金沙铅锌矿官房矿段西北部，形态呈似层状，总体倾向305°～46°，倾角11°～38°；矿体分布标高1 005～1 062 m，矿体厚3.04～6.80 m，平均厚5.94 m。

207矿体采矿方法为房柱法，沿矿体走向布置连续矿块，每个中段高20 m，矿块宽27 m。中段矿柱、矿块间柱和采场顶(底)柱尺寸均为4 m，矿房宽23 m，矿房高为矿体垂直厚度。随着开采的逐步推进，PD207、PD207-1坑道之间形成了采空区。经过多年回采，采空区的形态不断变化，空区规模逐渐扩大，其中一个采空区的最大跨度达67 m，采空区顶板暴露面积超过2 000 m2。采空区统计见表1，采空区位置见图1。

表1 采空区统计

由于采空区和上覆岩层在地应力作用下发生了移动[9]，导致207矿体上方出现了一定程度的地表变形，可能会引发地面塌陷、滑坡等灾害。因此，对207矿体地表进行了InSAR变形监测。以开采境界为依据，对矿区范围内受地下开采影响的地表区域进行变形监测。

图1 采空区位置

2 InSAR地表变形监测

2.1 InSAR技术原理

星载合成孔径雷达(InSAR)可实现对矿山采空区沉陷、边坡滑坡等重大灾害的远程监测、监控、排查、预警等[10]。SBAS-InSAR监测技术原理是将SAR影像在一定时间和空间上用很多小的集合来表示，然后采用最小二乘法将每一个小集合的地表相位进行解算，采用SVD联合求解的方式得到最小二乘解[11-12]。

本文通过小基线集(SBAS-InSAR)技术对矿区地表进行变形监测，具有能克服时空基线失相干问题、显著提高影像利用率等优点[13]。其主要通过对卫星参数及研究区地形参数进行雷达可视区划分，并结合SBAS-InSAR技术反演得到研究区的垂直变形场，进一步获取有效的观测区域，再结合矿区与地表的空间关系，综合分析识别潜在的风险区域。

2.2 数据选择与处理

本文采用的影像数据是Sentinel-1(哨兵1号)卫星获取的时间跨度为2018年3月22日-2020年9月19日的升轨Sentinel-1A数据。为了缩短数据处理时间，需对原始数据进行裁剪，裁去非研究区域；同时在数据裁剪时，为了提高数据处理的可靠性，应适当扩大研究区范围。为了更好地显示矿区每年的地表形态变化，3年的数据分3次处理。

2.3 数据分析

通过对3年矿区变形量图的分析，获得207矿体上方地表沉降量分布云图(见图2-图4)。由图2-图4可以看出：207矿体上方地表沉降量呈逐年增大的趋势，且采空区上方沉降量明显更大；随着时间的推移，2#和3#采空区上方地表出现了较大变形，2020年其沉降量达30 mm。

图2 2018年3月-2019年1月207矿体地表沉降量分布云图

图3 2019年1月-2020年1月207矿体地表沉降量分布云图

图4 2020年1月-2020年9月207矿体地表沉降量分布云图

通过对InSAR监测数据的分析，得到了整个207矿体上方地表的变形情况，识别其变形较严重的区域为2#、3#采空区上方地表。因此，在地表选取1-1′、2-2′两条剖面线(见图1)，在剖面线上每隔10 m选取一个特征点，通过对每个特征点的沉降时间序列图计算累计沉降量，拟合各特征点累计沉降量便可得到地表剖面线的沉降量曲线(见图5、图6)。

图5 1-1′剖面沉降量曲线

图6 2-2′剖面沉降量曲线

由图5、图6可知，地表沉降量最大值为38 mm，在剖面线上呈现出不均匀的凹形沉降，局部下沉非常明显。地表沉降区域位于2#、3#采空区上方，推测是地下开挖导致的地表沉降所致，故对2#、3#采空区开展采空区稳定性分析。

3 采空区稳定性分析

采用Midas/GTX建立地表和地下采空区三维可视化模型，然后将模型的单元、节点和组信息导入FLAC3D定义本构模型和材料参数并施加载荷及边界条件[14-15]。

3.1 采空区三维模型的建立

首先用Midas/GTX将AutoCAD地表等高线生成地表模型，然后确定2#、3#采空区与地表的空间关系，生成对应的地下采空区，建立实体模型。对实体模型进行分组，划分网格。

为提高计算效率，需对实际条件作适当简化：视缓倾斜的207矿体为平行于XY的平面，2#、3#采空区存在于矿体层。采空区高度取平均值5.5 m，顶板暴露面积按实际面积确定。最终构建了X轴长300 m，Y轴长260 m的实体模型。采空区单元尺寸为3 m×3 m×3 m，采空区矿柱单元尺寸为2 m×2 m×2 m，围岩单元尺寸为5 m×5 m×5 m；模型共划分为91 174个单元、87 538个节点(见图7)。

图7 三维实体模型

3.2 数值分析模型

将实体模型导入FLAC3D，综合分析采空区的应力、位移、塑性区的分布规律[16-17]。数值模拟模型如图8所示。

图8 数值模拟模型

1)模型基本假设

视所有岩体为各向同性的连续均匀介质，且仅受重力作用，忽略节理、裂隙、断裂带等因素对岩体稳定性的影响。

2)选取力学模型

摩尔-库伦模型适用于在剪应力下屈服但剪应力仅取决于最大、最小主应力，而第二主应力对屈服不起作用的材料[18]。本文选用摩尔-库伦模型模拟采空区开挖。

3)岩体力学参数

根据岩石力学试验结果确定各项物理力学参数，并采用Hoek-Brown强度折减法对参数进行折减[19-20]。最终计算模型所采用的岩体力学参数见表2。

表2 岩体力学参数

4)边界条件和初始条件

计算模型四周采用水平约束，限制模型侧面横向位移；底部采用固定约束，限制模型竖向位移；顶部为自由面。设置重力加速度为-9.82 m/s2，指向Z轴负方向。

3.3 数值模拟结果分析

分别对2#、3#采空区顶板位移、最大主应力、压应力和剪应力分布情况进行分析，以判断其是否存在失稳风险。将数值模拟得到的地表沉降过程曲线与InSAR监测得到的沉降量曲线进行对比，以判断采空区是否会继续引发地表沉降。

3.3.1 监测与数值模拟协同分析

根据1-1′、2-2′剖面线上的特征点在计算模型中布置相应的位移监测点。通过拟合数值模拟计算结果得到1-1′、2-2′剖面沉降过程曲线(见图9)，并与InSAR监测的剖面沉降过程曲线(见图10)作对比分析。

图9 1-1′剖面沉降过程曲线

图10 2-2′剖面沉降过程曲线

由图9、图10可知：2#、3#采空区在开挖初期就对地表产生了影响，导致地表沉降；在数值模拟计算结果逐渐稳定的过程中，采空区上方岩体变形和地表变形也逐渐趋于稳定；在采空区开挖的初期和中期，采空区对地表产生的沉降趋势与InSAR监测的沉降趋势具有一致性，但随着数值模拟计算结果稳定以后，2#、3#采空区导致地表变形的沉降量明显大于InSAR监测结果，表明采空区沉降尚未稳定，可能会导致地表沉降量继续加大。

3.3.2 位移分析

相关研究表明，采空区顶板位移会随着采空区跨度的增大而增大。因此，在完成数值模拟计算后，对2个采空区最大跨度方向取剖面，根据位移分布云图判断采空区的稳定性。

根据采场开挖后在无支护条件或临时支护条件下的现场经验[23]，得出顶底板围岩位移和采场稳定性的对应关系(见表3)。

表3 围岩位移和采场稳定性的对应关系

数值模拟计算后得到的采空区顶板位移分布云图如图11所示。由图11可以看出，2#、3#采空区直接顶板的位移分别为30～60 mm、30～70 mm，最大位移集中在采空区跨度中心区域，与工程经验相符。根据最大容许位移判断，当顶板位移大于50 mm时，采空区由稳定状态变为不稳定状态，存在失稳风险。2#采空区直接顶板最大位移达到了60 mm，3#采空区直接顶板最大位移达到了70 mm，均超过了50 mm。因此，判断2个采空区均存在因顶板位移过大而失稳的风险。

(a)2#采空区纵剖面位移云图

(b)3#采空区纵剖面位移云图

3.3.3 剪应力及压应力分析

Mohr-Coulomb强度准则认为岩体材料会因受剪应力作用而破坏，且其破坏所需剪应力和作用在破坏面上的正应力之间存在函数关系。该强度准则表达式为

(1)

式中：σ1、σ3分别表示岩体发生破坏时的最大主应力、最小主应力，MPa；C表示黏聚力，MPa；φ表示内摩擦角。

当作用在矿柱上的最大剪应力和最大压应力大于其抗剪、抗压强度时，矿柱存在失稳风险。最大剪应力和最大压应力模拟结果分别见图12、图13。由图12、图13可知，剪应力和压应力主要分布在矿柱上，且随着矿柱尺寸的减小而增大。2#、3#采空区矿柱上的剪应力为5.50～8.13 MPa，压应力为12.00～23.63 MPa。根据表2和式(1)可得矿柱的抗剪强度和抗压强度分别为7.29、18.65 MPa，采空区部分尺寸较小的矿柱上的剪应力和压应力超过了其抗剪、抗压强度，存在失稳风险。当某个矿柱失稳时，原有载荷将由其相邻矿柱承担，使得相邻矿柱上的剪应力和压应力极易超过其自身强度，继而诱发大面积的矿柱破坏，因此2#、3#采空区存在矿柱失稳风险。

图12 剪应力分布云图

图13 压应力分布云图

4 结论

a.利用InSAR技术监测金沙铅锌矿207矿体上方2018-2020年的地表沉降发现，2#、3#采空区上方地表的沉降量较大，且出现了不均匀沉降，存在采空区失稳风险；1#、4#采空区上方地表沉降量较小且均匀，暂无失稳风险。

b.数值模拟与InSAR监测协同分析结果显示，数值模拟中的2#、3#采空区导致的地表沉降量明显大于InSAR监测结果，说明地表沉降量还可能会继续增大。

c.数值模拟结果显示，2#、3#采空区顶板最大位移分别达到了60、70 mm，且采空区内部尺寸较小的矿柱上的剪应力和压应力超过了其抗剪、抗压强度，故2#、3#采空区存在失稳风险。

d.建议在后期采用废石充填的方式对2#、3#采空区进行治理，控制周围爆破作业的一次起爆药量，加强地表沉降监测，一旦沉降量出现突变，应立即停止采空区周围采场的回采作业。