随着地下金属矿山的开采不断增加，在采用空场法和嗣后充填法进行开采的矿山中，受开采技术条件等因素的影响，出现了大量采空区，这些采空区往往存在着隐蔽性强、不规整、空间分布规律性差和空间形态变化大等特点（马巍，2021；郭甲一等，2021；何荣兴等，2020），对矿山安全生产造成极大的威胁。因此，开展采空区的稳定性分析及隐患治理研究，对矿山安全生产具有重大意义。

目前，针对采空区的稳定性分析方法主要有4类，分别是理论分析、实验室相似模拟、数值模拟和现场监测及工程类比（刘海林等，2018）。对于大型复杂采空区群，理论分析方法往往很难解释各采空区之间的扰动影响；由于采空区的复杂性，又很难建立相似模拟试验，现场监测则需要较长的时间成本；而数值模拟可利用计算机在短时间内对大型复杂采空区进行计算分析，以应力、应变及塑性区等多种表现形式来解析采空区的稳定性问题。因此，采用数值模拟方法对采空区群进行分析成为采空区稳定性研究的最佳选择（王伟，2021；王悦青等，2019；李建红等，2019；陈光飞等，2019；李作栋等，2021；孔学伟等，2021；宋佳等，2020）。李杨等（2021）运用UDEC对钱家营煤层工作面采空区进行数值模拟，验证了“正—倒置三角形”垮落顶板形态存在的合理性；叶光祥等（2017）利用FLAC3D对矿柱资源回采进行数值模拟，通过对应力、应变及地表监测点进行分析，提出了矿柱最大安全回采的技术方案；韦文蓬等（2020）利用FLAC3D对大新锰矿西北采空区群进行稳定性分析，并根据数值计算结果对采空区安全处理提出相关建议；李鸿飞等（2019）利用FLAC3D对复杂采空区进行模拟分析，为采空区群治理提供了理论依据。Liu et al.（2015）采用FLAC3D模拟了无渗漏、普通渗漏和严重渗漏3种工况下巷道和采空区的变形和应力，结果表明垂直位移和影响深度随渗漏的增大而增大，但附加值减小；注浆使应力场重新分布，采空区两端存在应力集中现象；在渗漏工况下，水和渠道荷载对采空区无影响。Pan et al.（2020）利用FLAC3D模拟采空区和非采空区条件下工作面的回采过程，得到工作面前承载压力和塑性区分布规律，结果表明采空区下工作面采动时顶板来压步距和超前支承压力均小于非采空区下工作面，为采空区下开采提供了理论依据。Chen et al.（2016）采用FLAC3D建立了采空区封闭墙的三维数值模型，对不同条件下封闭墙的应力和变形进行了数值模拟，得到了不同条件下密封墙的应力场、位移场和塑性区的分布。

数值模拟计算结果的准确性很大程度上依赖于计算模型的精度，而模型的建立主要是依据地质资料和实测数据。近年来，三维激光扫描技术在矿山测量中得到了广泛应用（何赞碧，2021）。郭庆等（2021）对大红山铁矿采空区进行无人机三维激光扫描，构建出精确可靠的采空区模型；刘明江等（2020）对汤丹铜矿采空区进行三维激光扫描，构建出可直观获取采空区边界、体积、内部形态及分布结构信息的采空区三维地质模型，并将扫描计算结果与矿山实测统计数据进行了对比，数据测量误差为0.9%~9.4%，验证了三维激光扫描系统在复杂采空区探测方面的适用性。然而，在数值分析过程中，对于复杂大型采空区群，直接将三维激光点云数据转变为可计算的网格模型存在着网格数量巨大、计算时间长、无法收敛甚至无法计算等困难。

本文以某矿山大型复杂采空区群为例，使用无人机进行三维激光扫描探测，首先利用Geomagic软件对模型进行处理，保证模型精度，然后利用FLAC3D进行数值模拟分析并对计算结果进行验证，最后采用无人机三维激光扫描技术对采空区群的稳定性进行进一步分析，为矿山安全生产提供指导建议。

某矿山经过多年的开采，形成了大量的采空区。由于采空区存在时间较长，部分采空区的上下中段及相邻盘区发生冒落贯通，形成了大型复杂采空区群。因此，准确掌握采空区的位置和空间形态尤为重要。

综合采用井下手持式三维激光设备、架站式三维激光扫描系统和无人机载三维激光系统等多种方法进行了扫描探测，共探测出25个采空区，总体积超过150×10⁴ m³。由于采空区数量较多，为了能够从众多采空区群中快速找出大中型采空区，作为研究对象进行重点分析，根据采空区体积差异，将采空区划分为大型采空区（体积>40×10⁴ m³）、中型采空区（10×10⁴ m³<体积<40×10⁴ m³）和小型采空区（体积<10×10⁴ m³），并根据体积从大到小依次编号，如表1所示。

经实测，得到采空区群的分布范围：X=500 m，Y=320 m，Z=180 m。采空区群的空间关系及赋存情况如图1所示。

对于大型复杂采空区群的数值模拟，主要问题集中在采空区群模型的建立，尤其是网格模型的建立。由于采空区的形状极不规则，若采用常规手段建模，模型的形状与实际采空区群相差较大，进而影响计算结果的准确性。本文采用三维激光扫描获取采空区群模型，为保证模型的准确性，在测量时通过增加扫描点的密度，使得模型形状能够真实反映采空区群的空间形态。由于无人机载三维激光扫描系统获得的点云数据量增大，构建的采空区实体模型存在着三角片数量多、相交面多和模型粗糙等缺陷，导致无法建立网格模型进行计算，因此需要对模型进行优化以及局部简化。在优化及简化过程中，须保持模型的基本形态不失真，从而保证数值计算结果的准确性。本文利用Geomagic Control模型优化软件，对采空区模型进行降噪、删除钉状物等处理。

以3号中型采空区为例，初始模型如图2所示，可以看出该模型整体较为粗糙，三角片数量为1 274 534。

在保留初始边界的条件下，按最小尺寸重画网格，此时三角片数量降低为71 137。进一步通过去除钉状物、减噪等处理，使模型更加圆滑，最后进行二阶段重画网格，网格尺寸由最小逐渐增大，保证模型整体形态基本不失真。最终优化结果的三角片数量为19 659，如图3所示。

重复上述操作，依次对所有采空区进行处理，得到采空区群模型，如图4所示。

将优化后的采空区模型导入MIDAS软件，根据圣维南原理确定地下坑硐影响区为其开挖范围的3~5倍，从而建立了一个长3 528 m、宽2 268 m、高1 089 m的数值计算模型，该模型中网格单元数量为4 986 842，节点数量为3 085 923。

本文研究的矿山地质条件较为简单，矿体为磁铁矿，围岩主要为变钠质熔岩。根据矿岩物理力学参数与岩体质量评价结果，按照Hoek-Brown强度准则对矿岩物理参数进行折减处理，获得了相应的岩体物理参数，如表2所示。

将网格模型导入FLAC3D中，采用分阶段弹塑性求解法生成初始应力场；随后进行采空区开挖，求解结束后，获得采空区周边应力云图、位移云图及塑性区分布图，如图5所示。

由图5（a）~5（b）可知，各采空区顶板主要受压应力影响，且压应力较小，对顶板影响较小；而少数空区特别是1号大型采空区顶板出现较大拉应力，拉应力接近顶板围岩抗拉强度（0.230 MPa），极易产生拉伸破坏，发生冒落。由图5（c）可知，塑性区数量少且分布零散，仅少量塑性区发生局部贯通，存在一定的冒落风险。由图5（d）可知，1号大型采空区和3号中型采空区顶板垂直位移较大，极易发生冒落。

对于大型复杂采空区群的数值模拟，目前还没有形成统一的岩体破坏判据。在以往的研究成果中，主要结合矿山的实际情况（如岩体地质条件和采场跨度等因素），制定出矿山顶板失稳冒落评判标准。本文在参考相关文献的基础上，结合该矿山的工程地质条件和地压监测经验，发现该矿山采空区稳定性与顶板位移的相关度较高。综合考虑后，选择了顶板下沉量制定该判据，建立了符合该矿山情况的大型复杂采空区群顶板位移失稳判据（表3）。

针对3号中型采空区，由于该空区下方存在生产活动，顶板易受爆破震动等影响而发生冒落。根据此情况，将3号中型采空区的位移失稳判据下调，认为当该空区的顶板位移大于30 mm时，顶板极易破坏。

根据位移云图并结合位移判据，圈定出1号大型采空区和3号中型采空区顶板的隐患区域，如图6所示。由图6可知，1号和3号采空区存在着较大面积的隐患区域，顶板易发生冒落。

根据上述数值模拟方法，分析确定出1号和3号采空区存在顶板冒落风险的隐患位置。由于采空区群范围大，实际中1号大型采空区和3号中型采空区发生了局部区域冒落。

为了掌握冒落后的采空区群形态，对冒落区域及周边采空区情况进行了现场调查，并完成了三维激光扫描探测，获得1号大型采空区和3号中型采空区的冒落范围，如图7所示。

将图6与图7进行对比，发现数值计算结果圈定的1号大型采空区隐患区域位于边缘处（沿X轴、Z轴方向）且比实测结果略大，但中心冒落范围基本吻合；而数值计算结果圈定的3号中型采空区隐患区域位于下端（沿X轴、Z轴方向）且比实测结果略小，但上部主要冒落范围基本吻合。

数值模拟与实测结果出现差异化的原因主要是：①建立计算模型的过程中，对模型的形态进行了边缘合并等简化操作，导致局部范围出现差异性，3号中型采空区尤为明显，但主要冒落区并未出现大偏差；②顶板冒落过程是一个灾害孕育过程，具有明显的时效特点，在此过程中，顶板未冒落区域发生了应力重分布并达到稳定，导致计算结果与实测范围在边缘处出现偏差。尤其对于 1号大型采空区，根据现场监测，该采空区时常发生小范围冒落情况，此次冒落后扫描结果范围虽小于数值计算预测范围，但核心区域发生大面积冒落，整体预测结果较为准确；③3号中型采空区下方正在进行生产活动，爆破震动对下部顶板影响较大，导致在数值模拟结果中，3号采空区冒落范围偏小，致使该部位实际冒落范围大于预测结果。此外，采空区冒落时，由于岩体内部裂隙等结构较为复杂，很容易产生连续垮落，这也是3号采空区实际冒落情况大于数值模拟结果的一个重要因素。

综上分析，认为顶板冒落的实际位置与根据数值模拟计算结果圈定出的隐患位置基本一致，数值模拟结果与实际情况基本吻合，因此，可以利用数值计算结果来指导采空区安全处理方案的制定。

由于采空区群体积大、数量多且处理时间长，因此需对不同阶段采空区群的稳定性进行数值分析，从而提前预测隐患区域并制定出科学的防治措施。根据第2节所采用的方法，重新构建冒落后的采空区群模型并进行数值计算，得到冒落后的采空区群模型，如图8所示。

通过数值计算，获得采空区冒落后应力云图、位移云图及塑性区分布图，如图9所示。

由图9可知，随着1号大型采空区的冒落，大部分采空区顶板拉应力得到释放，而3号中型采空区发生冒落后，顶板跨度增加，导致顶板拉应力在重新分布后仍存在较大的拉应力，且拉应力影响范围也增大；同时，3号中型采空区顶板产生了较大面积的拉伸塑性区，极易发生冒落。与发生顶板冒落前相比，1号大型采空区和3号中型采空区的顶板位移增加，其中最大位移达到73 mm。根据表3中的位移破坏判据，1号和3号采空区极易发生冒落，其余采空区位移变化较小，整体上处于相对稳定状态。根据位移判据，并结合实际情况，圈定出1号大型采空区、2号中型采空区和3号中型采空区的隐患区域，如图10所示。

由图10可知，1号大型采空区、2号中型采空区和3号中型采空区顶板均存在较大范围的隐患区域。考虑到1号大型采空区和3号中型采空区周边赋存有较多其他中小型采空区，若发生冒落，极易引发相邻采空区的垮落连通，从而影响周边中小型采空区的稳定性，因此首先对1号大型采空区和3号中型采空区进行充填处理。待充填后，采用数值模拟方法对采空区群的稳定性进行数值分析，实现隐患区域的预测，从而为采空区群的隐患治理提供科学依据。

（1）通过无人机等各类三维激光扫描仪获得大型复杂采空区的点云数据模型，利用Geomagic-Midas对模型进行处理，建立了大型复杂采空区群数值计算模型。

（2）利用FLAC3D对采空区群进行数值计算，根据位移云图和位移失稳判据，确定出1号大型采空区和3号中型采空区为隐患采空区，并将模拟计算结果与现场冒落探测结果进行对比，验证了数值模拟方法的准确性。

（3）对冒落后的采空区群进行二次数值模拟，实现了隐患采空区的预测，为采空区群隐患治理方案的制定提供了指导。