基于现场监测与数值模拟的浅层采空区稳定性分析

贾瀚文 刘洪磊 张忠政 朱万成 王雷鸣 杨 震

（1.东北大学资源与土木工程学院，辽宁沈阳110819；2.东北大学岩石破裂与失稳研究所，辽宁沈阳110819；3.鞍钢集团矿业公司弓长岭有限公司露采分公司，辽宁辽阳111008）

在我国多空区矿山的露天开采过程中，随着台阶开采的不断剥离，露天开采境界内各台阶与地下空区群的隔离层厚度越来越薄，随时有可能发生采空区顶板坍塌事故，从而危及采场工作人员和大型采掘及运输设备的安全。而在多数矿区内，由于历史原因，大多数地下采空区往往资料不详，同时还存在大量的盲釆空区，且由于安全原因，这些采空区已无法进入。因此，露天采空区顶板监测预警以及对其稳定性进行安全评估，已成为当前多空区露天矿山生产中亟待解决的问题。

釆空区稳定性分析研究一直是个热点问题，也是一个极其复杂的问题，众多学者采用数值模拟法和现场监测方法进行了研究。徐恒等［1］采用尖点突变理论和FLAC 数值模拟进行了充填体下采空区稳定性分析。付建新等［2］建立了深部采空区卸荷失稳分析模型，引入局部能量释放率指标，采用数值模拟方法计算了矿房回采过程中局部能量突变判别值。朱卫兵等［3］采用UDEC 数值模拟方法对下煤层工作面采动时上覆房采煤柱群的动态失稳过程进行了研究。李世俊等［4］采用数值模拟与离心模型试验相结合的方法研究了采空条件下边坡的渐进破坏过程及破坏模式。汪文勇等［5］利用FLAC数值模拟方法分析了在回采过程中沿空留巷的应力分布以及变形特征。苏学贵等［6］利用FLAC分析了采空区及其临界巷道的稳定性并提出了相应的支护方案。王朋飞等［7］采用理论分析与数值模拟试验相结合的方法分析了非充分采动采空区与煤岩柱（体）的耦合作用机制。宋许根等［8］依据地质调查分析以及地表水平位移和垂直位移监测成果，结合采空区顶板崩落钻孔监测和崩落区高密度电法勘探等手段分析了程潮铁矿西区顶板崩落特征。赵勇等［9］利用微震监测设备研究了露天转地下矿山坑底采空区稳定性情况及其随时间的劣化过程。宋许根等［10］结合地表GPS、水准测量和三维激光扫描技术，分析了采空区周界变化与移动线和陷落线的对应关系。彭府华等［11］对降雨入渗诱发采场失稳的典型案例以及采场地压活动的影响规律进行了微震监测量化分析。王书文等［12］研究了工作面回采过程中采空区侧向支承压力演化及微震活动全过程的特征。范洪冬等［13］研究了一种基于多轨道SAR 影像的老采空区地表三维形变监测方法。胡静云等［14］利用多通道微震监测、手持式GPS 仪和全站仪等对上覆岩层崩落高度、上覆岩层下沉变形量、地表开裂范围、地表沉降与水平移动等进行了监测研究。

通过上述采空区稳定性研究成果分析可以发现，数值模拟方法和现场监测手段往往是独立存在的，彼此仅仅起到了互相验证的作用。数值模拟计算一般对岩体的现场情况进行简化，模拟结果与现场真实情况存在一定的误差，同时现场监测数据一般通过数据拐点的出现判断岩体失稳，无法进行长期预测。因此需要一种基于现场监测数据的数值模拟计算方法分析采空区稳定性，将现场监测数据作为数值模拟计算的输入端，从而可以充分利用两者的优势。本研究以弓长岭露天铁矿浅层采空区为工程背景，建立采空区地表沉降监测系统，将地表沉降数据反馈至数值模拟计算模型，综合分析该浅层采空区的稳定性。

1 工程地质概况

弓长岭露天铁矿位于为我国东北地区鞍山市东，矿山地理坐标：E 123°30'00″，N 41°07′40″。矿区最低温度可达-28 ℃，昼夜温差最大可达30 ℃。弓长岭露天铁矿分为独木采区、何家采区以及大砬子采区。目前独木采区共有6个生产水平，最高生产水平为+292 m，最低生产水平为+120 m。采区有4个铁矿层，分别为Fe1、Fe2、Fe3 和Fe4，其中Fe2 矿体规模最大，矿体厚度为34 m。最下层矿体为Fe1，最上层矿体为Fe4，矿层之间有岩石夹层，岩层平均厚度为25 m。Fe4矿层距地表150 m，全部为岩石覆盖层，矿体倾角为8°～25°，矿石平均地质品位32.42%，矿石密度3.3 t/m3，岩石密度2.6 t/m3。由于历史原因，弓长岭露天铁矿周边地方小矿点不规范乱采乱挖导致在独木采区形成了大量的浅层采空区，随着露天采场开采范围拓展与延伸，采空区随时可能发生塌陷，从而严重威胁人员及设备安全。为研究浅层采空区塌陷过程中的地表沉降规律，对其塌陷过程进行预测预报，有必要对该区域地表进行沉降监测。由于独木采区空区众多，因此选择该采区+145 m 平台的小北沟区域采空区作为典型空区进行地表沉降监测。小北沟区域采空区共分为两个区域（图1）：第一个区域为北部空区（N1～N6测点），该区域采空区长约65 m，宽约48 m，高度为6 m 左右，空区顶板厚度为22～26 m，顶板浮渣厚度为2～3 m；第二个区域为南部空区（N7～N12测点），该区域空区长约130 m，宽约40 m，高度为5 m 左右，空区顶板厚度为12 ～20 m，顶板浮渣厚度约5 m。

2 地表沉降监测系统建立及监测结果

2.1 静力水准监测系统建立

静力水准监测具有精度高、价格便宜且同时能实现区域性实时、自动、连续的垂直位移监测等优点。因此，静力水准监测技术被广泛应用于桥梁［15］、隧道［16］、地表建筑物［17］及地下厂房［18-19］的沉降位移监测。故本研究选用静力水准仪对典型采空区的地表变形规律进行监测。小北沟区域为浅层空区，其顶板厚度较小，因此在空区塌陷过程中的地表沉降值较小，同时+145 m 平台大部分区域内高差较小（约0.5 m），因此在N1～N10测点区域内选用量程为1.5 m 的传感器，在N11～N12测点区域内选用量程为2.5 m 的传感器。由于该区域最低温度可达到-28 ℃，为防止连通液在冬季结冰，连通液选用冰点为-40 ℃的55%乙二醇水溶液。

静力水准监测系统选用上海盛迪传感技术公司生产的SD-226 型静力水准仪，传感器测量方式采用24 h 自动无间断监测方式，每隔5 min 采集一次地表垂直位移数据。根据矿山静力水准测量原则，以第一个测点为参考点采用“线”形布置方式沿开采沉陷盆地主断面布置了12个水准仪测点。由于现场环境恶劣，故采用PPR 管对各线路进行保护，测点布置如图1所示，现场安装实物如图2所示。

2.2 监测结果

本研究HLS 采样时间间隔为5 min，各测点每日可分别采集288 个累计沉降数据。为了更直接地反映地表累计沉降值随时间的变化规律，对监测数据采用先剔除异常值后求日均值的方式进行处理。图3 为2018 年1 月27 日—2018 年5 月26 日N2～N12各测点数据处理后获得的累计沉降随时间的变化曲线。由该图分析可知，累积沉降的变化规律可以分为3类：第一类为N2～N6测点（图3（a）），在监测期间这5个测点的地表累计沉降值在-2～2 mm范围内上下浮动，可以认为该区域内地表没有发生沉降变化；第2类为N7～N9测点（图3（b）），这3 个测点的累计沉降值随时间变化呈不断增大趋势，N7测点与N9测点的累计沉降值在5 月26 日分别达到-6.6 mm 和-7.2 mm，N8测点最大累计沉降值达-9.8 mm，表明该区域地表发生了微小变形；第3 类为N10～N12测点（图3（c）），这3 个测点在监测期内沉降值波动较大，其中，N10测点的监测值在-5～5 mm范围内波动，N11、N12测点的监测值于2018年3月9日后在沉降值10 mm处上下波动。

3 基于监测数据的FLAC 模拟空区稳定性分析

3.1 模型建立

根据矿山地质资料建立了整个矿区地质模型，计算模型长800 m，宽600 m，模型自上而下厚度不等，均符合实际情况。其中模型长度方向为X 轴方向，宽度方向为Y 轴方向，高度为Z轴方向，高度方向底面为0 m 水平，顶面随地表起伏高度不等，最高312m、最低217 m，共划分749 241 个单元，129 773 个节点。根据岩体的岩性、结构类型、岩体强度、岩块硬度等方面的不同，可将矿床的工程地质条件大致分为5 个区：（I）铁矿、（II）第四系、（III）页岩、（IV）花岗岩、（V）角闪岩，如图4所示。

3.2 参数选取

本研究数值模拟参数采用的是在前期大量室内岩石力学试验得到的岩块基本力学参数基础上，结合现场结构面赋存情况，对岩体稳固性进行分级，而后通过Hoek-Brown 强度进行计算，最终选取的参数如表1所示。

?

3.3 采空区蠕变计算分析

考虑到Norton-Bailey蠕变方程参数少、简单容易获取，故本研究采用该方程进行蠕变计算，根据Wang等［20］研究成果，可得到其张量形式为

该式即为复杂应力状态下的岩石蠕变本构方程，式中，A、m、n为经验常数；σe为有效应力；t为蠕变时间；Si,j为偏应力张量为蠕变应变率。但以上所得到的蠕变模型只能描述蠕变的前两个阶段，无法得到加速蠕变阶段。为了得到蠕变全过程曲线，需要建立一个非线性加速蠕变模型。因此，本研究在蠕变计算的基础上，考虑单元达到破坏准则后会发生损伤，分析方法采用与RFPA 程序类似的弹性损伤本构关系［20-21］。首先假设岩体的某些力学参数，如弹性模量、强度等为非均质分布，以反应天然岩体介质的非均匀性；然后选取摩尔-库伦破坏准则对岩体各节点位置的应力状态加以判断，对于满足破坏准则的位置，基于损伤力学理论，对其力学参数进行相应调整，从而实现蠕变第三阶段——加速蠕变阶段。

该蠕变方程为三参数经验方程，其中A、m、n 等参数取值与蠕变曲线的应力、稳定蠕变速率、减速蠕变速率相关。本研究首先根据现场监测数据与数值模拟结果对比，调整蠕变参数，使得数值模拟的蠕变速率与现场监测结果一致，然后据此进行未来结果预测。本研究参数最终取值为A=1.0×10-12、m=1.75、n=0.35。

3.3.1 空区顶板力学特性时变演化规律

在实际情况中来看，采空区形成后，空区结构力学特性随着时间发生劣化，是导致地表塌陷和空区坍塌的重要因素之一。本研究数值计算过程分为2个阶段，分别为回采过程和流变过程。为了考虑时变力学的影响，在蠕变模拟时，将蠕变过程按年分为若干模拟阶段，在每一阶段中设置时间参数，从而实现时变模拟。

在Y=230 m 处设置剖面进行采空区应力和位移分析，图5（a）为蠕变5 a 后采空区顶板位移的变化情况，位移最大值为20.4 mm，最大位移发生在南部较大空区。图5（b）和图5（c）分别为蠕变5 a 后最小主应力与最大主应力云图，由于空区埋深较浅，故产生的最大与最小应力值均较小，应力最大值不超过25 MPa。

3.3.2 空区上方地表变形时变规律

地表下沉是浅埋空区引起的主要灾害之一，当变形过大时就有可能引起地表塌陷，带来严重的安全隐患［22］。因此，有必要分析浅埋空区引起地表变形的时变规律。本研究面临的实际问题在于，弓长岭露天矿的不少空区是早期形成的，其开始形成时间即空区顶板开始发生明显蠕变变形的时间难以确定。鉴于此，本研究采用蠕变损伤数值计算与现场监测相结合的方式实现。对长达0.5 a 的地表监测获得的沉降变形量求取平均沉降速率，再将其与数值计算中相应监测点的模拟沉降速率对应。

根据矿山实际情况，在采空区上方地表布置了12 个水准仪进行地表位移监测。在一段时间的监测后，发现N8点的实际监测位移值最大，在该点位移的基础上，按年度开展蠕变计算，预测未来5 a的地表位移情况。蠕变过程中，在空区上方地表对应的监测位置均设置了数据提取点，图6为水准仪N8点实际监测的位移值，图7为N8数据提取点蠕变5 a的位移值，图8为空区上方12个数据提取点蠕变5 a的位移变化情况。

由图7 可知：空区上方岩体变形逐渐趋于稳定，未来5 a 年内不会产生加速破坏。由图8 可知：随着蠕变时间的增加，地表下沉情况开始比较明显，而且越靠近采空区中部，位移越大。其中N8数据提取点处在蠕变计算5 a 后，位移为8 mm，而全部12 个数据提取点中，沉降位移值最大达到10 mm。考虑到空区地表没有出现大范围失稳破坏以及加速蠕变情况，认为采空区上方地表在回采5 a内不会发生明显的沉陷变形。

4 结 论

以弓长岭露天铁矿为工程背景，选择小北沟区域采空区为典型浅层空区，采用现场监测和数值模拟相结合的方法对该浅层空区进行了稳定性分析，得出如下结论：

（1）建立了小北沟浅层采空区静力水准沉降监测系统，仅N7、N8、N9测点处地表有微小变形，3 处监测点的变形值分别为-6.6、-9.8、-7.2 mm，现场监测数据表明该监测区域采空区较为稳定，有待于进一步监测。

（2）运用FLAC 数值模拟方法进行了小北沟浅层采空区的蠕变计算分析，计算结果显示，蠕变5 a后采空区顶板位移最大值发生在南部空区为20.4 mm，但应力值最大不超过25 MPa。

（3）根据“基于现场监测数据的数值模拟分析方法”的研究思路，将数值模拟地表监测点与静力水准仪监测点的沉降速率进行对比，获取数值模拟中采空区顶板开始发生明显蠕变变形的时间，从而进行蠕变求解，数值计算结果表明，12 个监测点5 a 蠕变地表发生的最大位移约10 mm，地表不会出现大范围失稳破坏以及加速蠕变情况。