代碧波 孙丽军 赵兴东

（1.东北大学采矿地压与控制研究中心，辽宁沈阳，110004；2.金属矿山安全与健康国家重点实验室，安徽马鞍山，243000；3.中钢集团马鞍山矿山研究总院股份有限公司，安徽马鞍山，243000；4.武汉科技大学资源与环境工程学院，湖北武汉，430081）

缓倾斜薄矿体的非充填法回采会形成大量的采空区，这些空区群的稳定性直接关系到井下生产的安全［1-2］。故而对空区群实际三维形态的测定和稳定性评价至关重要［3］。

采空区的实际三维形态和准确的空间位置是影响空区稳定性数值模拟结果准确性的关键因素。过去，在对空区稳定性进行数值模拟时，往往将采空区的形态视为空区形成之前（即开采前）进行采场设计时确定的规则形态，部分则是依据传统的测量方法测定的简化空区形态进行模拟计算［3-5］。然而，采空区的实际形态往往在实际开采过程中，因测量、爆破等因素的影响而呈不规则状，加之采空区形成后受爆破震动的影响，以及空区围岩的冒落、片帮，致使空区的实际形态和开采设计中确定的相差甚远、这就势必严重影响空区稳定性数值模拟结果的精度和可靠性［6-8］。因此，如何精确地确定空区的实际形态和空间位置，并以此为基础建立空区群稳定性数值模拟模型进行数值模拟计算分析，是空区群稳定性数值模拟研究中应加以解决的一个重要问题［4］。

针对上述问题，本项目开展以精密探测获得的空区实际三维形态和空间位置为基础的空区群稳定性数值模拟研究，形成一种新型的基于空区实际三维形态精密探测的空区群稳定性数值模拟技术，以提高数值模拟结果的精度和可靠性，并结合矿山井下开采现状进行采空区的精密三维数字化及数值模拟技术研究，将所形成的技术应用于采空区群稳定性数值模拟分析中，并在此基础上对其进行稳定性分析，为空区群的综合治理提供技术支撑。

1 采空区探测技术

对于地下不明的采空区探测，我国目前以钻探为主，物探为辅。常用的物探方法主要有以下几类：①重力测量；②电法测量；③地震测量；④MT法、CSAMT法和EMAP法；⑤放射性测量。由于物探手段对采空区探测不能明确空区体积与形态，这些方法只能是大概探测地下有无空区以及空区的大致位置［9］。近年来发展起来的激光3D扫描法是利用激光测距原理对采空区进行精密探测的一种技术手段，只要能将激光探测设备伸入到空区内部，便可对空区进行360°扫描，经过后处理建立空区三维实体模型便可获取空区的具体参数，包括空区边界坐标，容积等［8，10］。

2 空区群三维数字化

2.1 工程概况

本研究依托的矿山为中钢南非铬业有限公司Dilokong铬铁矿，该矿山特点：①矿石品质优良、开发利用价值巨大；②矿体倾角缓、厚度薄，仅开采单一矿层，开采环境差、效率低；③矿山空区分布广泛，仅采用木支护，井下开采安全条件差。Dilokong铬铁矿目前开采的资源为LG6铬矿层。LG6矿层厚度约1.0～1.2 m，上覆有厚度约0.5 m厚的LG6a矿层，中间辉岩厚度大约0.5～0.6 m。区内矿体变化不大，LG6与LG6a矿层连续性较好，成矿后断裂影响不大，似层状产出。矿体整体走向南北，倾向西，倾角14°左右，浅部矿体较陡，约18°，深部倾角较缓，约7～11°。矿体遍布矿权范围，长约7 800 m，宽约500～3 700 m，埋深约在300～400 m，矿体规模巨大。矿石Cr2O3品位为41.07%～49.3%，平均品位44.62%，品位变化均匀。矿体密度为4.25 t/m3。

2.2 矿区地表三维数字化

根据矿区地形地质图，对矿区地形等高线进行三维坐标赋值，然后根据三维空间关系生成面域，并根据高程范围进行着色渲染和水流方向的绘制。经三维数字化的地表模型如图1所示。

2.3 采空区群三维数字化

根据矿山开采实际，将已有采空区群进行精密三维数字化，实现空区群的数值三维体态查看和分析。空区群位于名称为Hanekom的老矿山范围内。根据开采现状，矿山共进行了4个中段的回采，形成了大量的采空区。根据采空区位置和开采现状，利用激光3D扫描进行空区数字化，建立空区群的三维形态模型，见图2。整个矿区空区群总体积为135.708万m³。

3 空区群的数值模拟分析

应用三维拉格朗日差分分析程序模拟分析空区群的稳定性、顶板覆岩的变化规律以及留设矿柱参数对采空区稳定性所产生的影响［11，12］。

3.1 三维计算模型的建立

本次计算结合矿山工程地质条件、矿体赋存条件及所用采矿法的特点，以三维化空区群为依据建立模型，从采空区顶板的应力变化和位移，矿柱承载能力变化、采空区顶板应力及位移状态和地表变形等方面进行空区群稳定性的数值模拟。最终确定的模型如下：模型长度3 100 m、宽度3 000 m、高度1 100 m，模型共有287 928个单元，434 768个节点。生成主体计算模型之后，其中空区群模型的构建采用实测三维数字化后的立体模型与本主体模型进行耦合，如图3所示。

3.2 岩体力学参数

根据矿山的工程地质特征，矿体和围岩特性，本次数值模拟过程中选取了对模型计算有影响的3种矿岩体，即辉岩、斜长岩和矿体。数值模拟过程中所用岩体力学参数见表1。原岩应力按常规进行假设计算，忽略构造应力，按照自重应力计算。

3.3 数值模拟结果

对建立的空区群模型进行计算，空区群内留设10 m×6 m间柱，按照下行式和前进式的回采顺序逐步形成空区群，并在每个回采步骤记录矿岩体应力和位移变化特征，经过模拟分析，得到了次生应力平衡之后的最大主应力云图、最小主应力云图、位移云图。以采空区中央为中心，创建剖面，得到剖面位置的最大主应力云图、最小主应力云图、位移云图。其数值分析模型计算第1步骤回采和第12步骤回采结果如图4和图5所示。

3.4 模拟结果分析

由以上数值模拟结果可知，原始地应力呈现层状分布，在进行回采过程中，明显扰动了原岩应力场，采空区部分的上覆岩层中最大主应力和最小主应力均存有共性，呈现似层状分布，在临近空区的矿岩交界处应力变化明显，应力变化沿矿体倾向逐渐过渡。

进行第1步骤回采时，最大主应力为8.44 MPa，出现在上盘与底板的交界处，并在角隅处伴随有应力集中现象，最大位移出现在下盘中央部位，为12.81 mm，表现为底鼓现象。进行第2步骤回采时，次生应力场的分布规律与第1步骤的相似，最大主应力为15.44 MPa，出现在下盘与顶柱的交界处，最大位移为16.99 mm，出现在本步骤内的上盘中央位置。

后续回采过程中，最大主应力值和最大位移值随着回采步骤逐渐增大。在回采包括第10步骤在内的矿体时，最大主应力均出现在本步骤上盘和底板的交界处，表现为应力集中现象，但进行第11和第12步骤回采时，最大主应力始终在第10步骤中的上盘和底板交界处，至第12步骤回采结束并达到次生应力平衡时，最大主应力为66.14 MPa。最大主应力值随回采步骤的变化如图6所示。

在第5回采步骤之前，最大位移出现在上个步骤内的上盘中央位置，最大值为33.14 mm，第6至第8回采步骤中，最大位置始终在第4步骤的上盘中央位置，最大值为40.85 mm，第9回采步骤至第12回采步骤中，最大位移始终在第9步骤上盘中央位置，至回采结束时，最大位移为116.06 mm。最大位移随回采步骤的变化如图7所示。

从图7可以看出，最大位移在第10回采步骤后趋于稳定。

另一方面，从塑性区计算可以看出，回采至第9步骤时，在上盘中央位置出现拉伸屈服区域，在计算结束后，大部分退出屈服状态。至回采结束，未见有塑性区的贯通现象。表明现有空区群表现为稳定状态。

由上述模拟结果可知，在形成空区群过程中明显扰动原岩应力场，在相同条件下，采空区最大压应力和最大拉应力随矿柱宽度增大而减小，随采空区顶板暴露面积的增大而增大。在数值模拟过程中，应力分布状态以及塑性区发生了比较明显的变化。当矿柱达不到与采空区顶板暴露面积最佳耦合的尺寸时，模型在计算过程中出现了屈服区域，计算结束后仍有较大范围的屈服区域，此时采空区处于不稳定状态，有发生顶板冒落的危险，当矿柱达到与采空区暴露面积最佳耦合状态的尺寸时，有部分采空区围岩在开采过程中曾经进入到拉伸和剪切屈服状态，但在次生应力平衡之后又退出了屈服状态。

另外，当增大采空区顶板暴露面积时，采空区上覆岩层的压力传递到矿柱上，使得顶板中的压力拱跨度增大，在采空区顶板产生拉应力区域，当此压力拱跨度大于合理拱轴线对应拱形跨度，拉应力超过岩体抗拉强度时，便会在顶板位置产生拉伸屈服区域，造成采空区的破坏。在矿柱尺寸一定的情况下增大顶板暴露面积时，采空区顶板中央的下沉量也随之增大，矿柱中的位移也明显增大，与理论分析结果相符。

4 结论

（1）按照现有空区群内留设10 m×6 m间柱的情况下，经过数值模拟分析，空区群表现为稳定状态。

（2）数值计算模型中，矿柱靠近顶板的位置为应力集中部位，并且由于矿体缓倾斜，使得应力集中部位偏向于矿柱靠近上部空区一侧。

（3）给空区群模型留设合理尺寸的矿柱时，采空区顶板和矿柱在模型计算结束后均没有拉伸屈服和剪切屈服区域的贯通，并且采空区顶板中央的拉应力均没有超过矿体抗拉强度，顶板和矿柱内的位移趋于定值，并处于工程可控范围，因此处于稳定状态。

（4）采空区顶板中央的下沉量及矿柱内的最大主应力和位移对矿柱尺寸表现敏感，并且顶板暴露面积越大，此敏感度越高，表明在空场条件下，矿柱尺寸对采空区群的稳定性影响很大。