关守安， 李 皓， 金长宇， 刘 冬

(1.东北大学 资源与土木工程学院， 辽宁 沈阳 110819； 2.辽宁省安全科学研究院， 辽宁 沈阳 110004；3.汕头大学 工学院， 广东 汕头 515063； 4.汕头大学 智能制造技术教育部重点实验室， 广东 汕头 515063)

矿产资源是国民经济的重要物质基础，为社会的发展提供最基本的原材料，是人类社会发展和文明进步的重要组成部分.矿产资源开发分为露天开采和地下开采，地下开采占全部开采的 80%以上，其开采深度一般为数百米，而且随着对矿石需求增长、采矿技术的发展，开采深度会越来越深[1]，与此同时，地下开采带来的岩层移动和地表塌陷问题也日趋严重.因此，矿山资源开采导致的岩层移动成为岩石力学领域的热点研究问题.

为了研究地下矿山开采对覆岩运移规律的影响，国内外诸多学者开展了大量研究.Greco[2]于1996年便提出了基于Knothe时间函数的地表移动变形预测模型.随后，一些学者对不同采矿方法下覆岩运移规律进行了统计分析[3-5].随着计算机技术的快速发展，数值模拟法被广泛运用到岩层移动和地表沉陷问题的研究中[6-8].可见，已有的研究成果侧重理论分析与数值模拟，而数值模拟方法也多采用连续介质力学分析手段，较少利用可以更真实模拟岩层移动的非连续力学以及随机介质的方法进行分析.

本文基于粗榆金矿的急倾斜极薄矿体开采工程背景，利用多种不同数值模拟手段开展采矿工程中上部岩层运移规律研究，所得成果为急倾斜极薄矿体的矿山开采有效控制岩层移动问题提供重要指导.

1 工程概况

粗榆金矿位于吉林省磐石市，属于急倾斜极薄矿体，采用干式充填采矿方法.-325 m标高以上为空区，品位低的采场留设矿石顶柱.形成的空区沿矿体走向及延深方向上不连续，最长空区走向长度约为360 m.区内矿体厚度为极薄至薄矿体，所形成的空区在一定深度及暴露面积下通常稳定性较好.在空区走向方向，无大的构造弱面相交.

由于原明竖井和盲竖井都在圈定的岩石移动带之内，原选厂和职工宿舍等基础设施也需废弃，变更计划后在新明竖井处计划建一选矿厂，按照原本圈定的岩石移动带，选厂离移动带不足20 m，不符合安全规定.

由于我国金属矿山设计规范关于地表移动带的规定一直沿用前苏联矿山设计标准，导致很多金属矿进入深部开采后，地表岩石移动界限范围巨大，致使设计的地表岩石移动范围与我国诸多大型矿山在深部开采中的实际情况相悖.如果依然按照这一设计标准进行岩石移动带的确定会导致多余征地，甚至地下矿体无法开采等问题，造成土地资源与矿产资源的浪费.为此，本文利用岩石力学分析方法对粗榆金矿地下矿体开采过程的岩移问题进行详细分析，通过不同的计算方法划定岩移界限，为粗榆金矿地下开采提供理论指导.

2 现场节理调查

岩体是由结构面和岩块构成的复合体，因此结构面对岩体稳定性的影响至关重要.本文利用CAE立体图像采集仪对1#～4#测点进行了节理裂隙统计，为围岩稳定性分析提供理论依据.

1#测点合成后的岩体三维壳体如图1所示.该测点范围岩体所揭露节理的极点图和玫瑰花图如图2所示.

图1 1#测点岩体三维壳体

图2 1#测点范围内的结构面产状信息

根据数据聚类结果，对测点范围内节理的产状以及平均节理间距进行分析，结果如表1所示.

表1 节理产状调查结果

3 岩体力学参数确定

通过现场实际考察，在X4矿脉的+126 m和+99 m中段、-181 m中段的上下盘以及X3矿脉-473 m开拓系统内开展了矿岩体的现场采样工作.该采样工作尽可能靠近采场，保证后续计算的评价工作能够反映采场围岩的破坏情况.

通过对标准矿岩样进行室内岩石力学测试工作，得到表2所示的矿岩石力学参数.

表2 矿岩石力学参数

结合前期地质资料、节理三维扫描与室内力学试验的结果，确定了Hoek-Brown强度准则中的评价指标，计算得到了矿岩体力学参数，分别如表3、表4所示.

表3 Hoek-Brown强度准则所需指标

表4 矿岩体力学参数

4 数值计算分析

4.1 边界条件

数值计算的边界条件采用位移控制，通过神经网络反演得到，具体方法见文献[9],得到的数值计算的位移边界条件和重力修正系数，见表5.

表5 反演得到的位移边界条件和重力修正系数

4.2 基于连续介质力学的计算

考虑到粗榆金矿以浅部硬岩为主，因此FLAC3D计算中主要采用了摩尔-库伦本构模型.建立大型的三维精细数值模型分析粗榆金矿开采过程中的地表沉陷问题，其中矿体的开采最深水平是-473 m，矿体的上覆岩层厚度达到400 m左右，模型的大小为883 m×700 m×1 290 m(长×宽×高).计算范围的坐标为x方向在665.6～1 548 m，y方向在492.8～1 191 m，z方向在-789.4～500 m，模型的计算原点坐标对应大地坐标为(42 502 229.48，4 773 992.64，0).计算单元为1 116 185个，节点为186 557个，其中模型计算参数参见表4与表6.

表6 碎石充填体计算参数

根据矿体埋深、厚度、地质条件以及矿岩体的变形、应力大小，选取了某个典型断面进行分析，典型断面的位置如图3所示.

图3 典型断面位置

根据数值计算结果，统计下行开采过程中不同水平的最大竖直位移，如图4所示.当开采深度≤-325 m时，随着开采深度的增加，各水平的最大竖直位移逐渐增大；当开采深度>-325 m时，随着开采深度的增加，最大位移量增长缓慢，最终趋于稳定.当开采至-473 m水平时，矿区地表的最大竖直位移为2.58 mm.图5给出了-5 m岩层随开采深度的倾斜量变化曲线，由图可知，地表的最大倾斜量为0.015 mm/m.根据我国《有色金属采矿设计规范》中的规定，I类建筑物(如竖井)的倾斜量应小于±3 mm/m.因此，可认为粗榆金矿的地下开采活动对地表沉降变形影响较小，不会导致明显的沉陷问题.

图4 不同水平最大竖直位移变化曲线(FLAC3D计算)

图5 -5 m岩层倾斜量变化曲线(FLAC3D计算)

4.3 基于非连续介质力学的计算

为了分析粗榆金矿下行回采过程中围岩与矿体的非连续力学行为，本文采用3DEC进行数值模拟分析.

根据表1中节理产状的测试结果，对图3中典型断面进行节理切割，建立准三维模型，模型水平向长度为883 m，模型上下高度为1 248 m，厚度为8 m，如图6所示.基于文献[10]中经验公式，确定了粗榆金矿结构面参数的取值，如表7所示.其中矿岩体参数见表4.

图6 3DEC数值计算模型

表7 粗榆金矿结构面参数的取值

图7给出了下行开采过程中不同水平的最大竖直位移，其变化规律与FLAC3D计算结果基本一致，当开采至-473 m水平时，矿区地表的最大竖直位移为3.54 mm.图8给出了-5 m岩层随开采深度的倾斜量变化曲线，由图可知，地表的最大倾斜量为0.003 5 mm/m，小于±3 mm/m.分析-440～-473 m阶段开挖完毕后的变形云图(图9)可以看出，矿体开采后在其顶部围岩内形成明显的拱形承载结构，受压力拱作用影响，围岩及地表变形均不明显，所产生的变形主要集中在矿体上盘和顶部，呈筒形陷落形态，且向上盘倾斜.

图7 不同水平最大竖直位移变化曲线(3DEC计算)

图8 -5 m岩层倾斜量变化曲线(3DEC计算)

图9 -440～-473 m阶段开采完毕后位移云图

4.4 基于随机介质方法的计算

利用东北大学自主开发的矿体开采地表沉降模拟系统[11-12]对粗榆金矿开采过程中岩层运移规律展开研究.该系统基于随机介质理论，可构建不同地质模型，制作不同形状、大小和倾角(产状)的采空区，并对沉降位移、沉降范围以及时间进行计算分析.

计算方法采用散体移动时空演化模型.将散体分成正方形模块，用模块之间从下向上随机递补运动，来模拟矿岩的运动过程.模块之间的递补通过“空位”向相反方向的随机传递来实现.每沉降1个模块就在其原来位置产生1个“空位”，该“空位”由其上面相邻9块模块按给定的概率和时间因素随机递补，在递补模块下移之后，其原来的位置又变为“空位”.因此，随机介质方法更适合对急倾斜极薄矿体开采导致的岩层运移进行模拟研究.

图10为松散系数1.5时的模拟结果.由图可知，粗榆金矿地表潜在沉陷区主要集中在矿体上盘区域，其上盘岩移角在88°左右，下盘岩移角接近90°.如图11所示，计算结果较原设计资料(中上盘围岩移动角75°，下盘围岩移动角80°)岩移范围大大缩小，其中下盘岩移界限向沉陷中心移动了215 m，上盘岩移界限向沉陷中心移动了195 m，沉陷区面积缩小了75%.综合上述各项评价指标，2#竖井、新建的选矿厂以及职工宿舍浴池均在岩移带之外.

图10 随机介质方法计算结果

图11 移动带的确定

5 结 论

1) 连续、非连续介质力学计算显示，粗榆金矿开采过程中，矿体上盘岩层沉降明显，岩层倾斜量小于允许值.根据承压拱理论，矿体上覆岩层稳定性较好，矿体覆岩沉降变形区域形态呈筒形.

2) 随机介质方法计算得到的岩移范围较原设计资料缩小75%.可判定目前二号竖井、新建的选矿厂以及职工宿舍浴池均在岩移带之外.

3) 岩体破坏过程是由连续变形向非连续断裂破坏转变.本文通过三种计算方法表明：随机介质方法在模拟急倾斜极薄矿体开采导致的岩层运移中更具优势.