精准三维激光探测结合可视化建模的采空区稳定性分析

孙刚友， 张卫中， 康钦容*， 夏缘帝， 袁威， 胡清珍

(1. 中核第四研究设计工程有限公司， 石家庄 050021； 2. 武汉工程大学资源与安全工程学院， 武汉 430073)

地下开采过程中，由于采空区的形成，原有的应力平衡被打破，上覆岩层因失去支撑而产生移动变形，导致采空区塌落，地表下沉等[1]。随着矿山的继续扩大，向其更深部位置进行开采后，岩层稳固性逐渐下降,致使该矿山围岩及开采后作业场地条件逐渐恶化, 在强大且复杂多样的围岩地压作用环境下,极易造成采空区的矿山顶板坍塌、冒落等事故。因此，对采空区进行稳定性分析已成为当前矿山生产中厄待解决的问题，其研究对保证矿山持续发展及安全生产具有十分重要的意义[2-5]。

针对采空区稳定性问题，国内外学者对其进行了大量的研究工作，前期采用的方法主要有模糊数学法[6-7]、Bayes判别法[8]、矢量分析法[9]、有限元法[10]等数学方法。由于数值模拟技术在工程等应用领域中的迅速发展和应用,以数值分析技术作为主要分析方法的研究层出不穷。如贺小庆等[11]为解决双王金矿KT9采空区稳定性问题，以FLAC3D为技术平台对采空区周围岩体应力场和位移场的变化情况进行了系统分析,明确了采空区稳定性情况,为后期矿柱回采和采空区处理提供了理论依据。尚振华等[12]针对规模较大采空区的稳定性问题，主张从量化的角度对其进行分析，在数值模拟的基础上结合统计学中的概率基本原理，以FLAC3D计算处理后的数据为基础，建立指标体系。从多方面深入研究采空区的破坏概率，极大地扩展了数值模拟的应用空间。随着研究的深入，单一的分析方法无法精准把握空区稳定性情况。为此，杨勇等[13]在解决金属矿复杂采空区稳定性问题时，以红岭金属矿为背景，选择将现场监测、数值模拟以及经验公式相结合的分析方法，对采空区的稳定性进行综合判定。但传统的空区稳定性数值分析总是基于某一特定范围确定的，且将采空区的边界视为规则的[14-18]，但规则的空区形态往往因外部因素影响呈不规则状，因此，采用传统的分析方法将直接导致数值模拟结果的精确度受到影响。

现以中核赣州某铀矿的818-2#采场采空区为研究对象，根据铀矿采空区的特点采用BLSS-PE三维激光采空区测量系统对复杂采空区进行多点探测，运用DIMINE软件构建三维实体模型，将已建立的采空区三维模型与FLAC3D进行耦合得到基于实测的采空区数值分析模型，经数值模拟计算后，得到应力场、位移场的变化情况，进而分析采空区对周边围岩的稳定性影响，解决实测空间与设计空间不一致的问题，为复杂矿井采空区稳定性分析提供新的研究途径。

1 工程概况

Kt-818-2#矿体位于710#剖面，走向近东西向，长约40 m，矿体平均厚度5 m，赋存高度166～200 m，倾向南，向北倾，倾角约72°。矿体围岩形态与岩矿石类型特征大体一致,主要成分以中粒二云母质花岗岩和中粒小斑状黑云母质花岗岩为主,矿体围岩结构较稳定,岩矿石成分的普氏系数为5～11,岩石成分的平均自然安息角约为37°,矿石容重约大于2.6万t/m3,矿石松散系数为1.65。矿体与围岩无较明显的界线,呈渐变关系，按照正常回采进度，该采场9个月应该回采结束，但后期在回采过程中，发现200 m中段矿化较好，若留顶柱回采，后期损失量较大，为了尽可能回收金属，避免浪费，后期铀矿回采至194 m标高后，又继续上采了16 m，回采至200 m中段至210 m标高，金属回收率很高。以818-2#采场为例，通过数值模型的建立，对矿区围岩的应力、位移以及塑性区分布等方面进行了分析，并估算了采空区的稳定性。技术路线如图1所示。

图1 技术路线图Fig.1 Technology road map

2 三维采空区模型的构建

2.1 原始探测数据的预处理

运用BLSS-PE三维激光扫描测量系统将818-2#采场采空区进行了2次扫描后，获取采空区空间形态位置的点云数据，将点云文件导入BLSS-PE测量系统中，生成“.blss”的点云文件。在三维可视化环境下呈现采空区及附近巷道的空间形态。利用三维激光扫描系统完成采空区空间形态数据的探测处理后,在空区模型的构建之前,必须对已得到的点云数据进行其他必要的数据分析处理，其主要分为以下步骤。

步骤1数据的过滤。通常情况下，经空区探测后的数据由于探测过程中受到外界因素以及空区自身环境影响导致所测数据中存在无效点，因此，在探测后为保证数据的精准度，通常先对数据进行过滤处理。

步骤2数据的拼接。由于空区存留的矿石或矿柱的阻挡作用，导致得到的点云数据在某种程度上可能无法还原真实的空间形态。针对空区的特点，对其采用多点探测法进行探测，再通过建立坐标系的方法将三维点云数据进行拼接。

步骤3数据的拟合。将经过多次探测所获得的点云数据进行模型拟合，以便于使最终获得的点云文件更真实地反映出这个空区的空间形态,而只有当空区点云文件完成了数据的前处理之后，才能进行地下采空区三维模型的构建。

2.2 三维采空区模型构建

将点云数据转换成可导入DIMINE软件的格式后，采用DIMINE软件中的实体建模功能构建采空区的三维图形轮廓，实现采空区的三维可视化。818-2#空区实测三维模型和实测编录模型的拟合图如图2所示，可以看出，空区存在一定的垮塌情况，但经三维激光探测建立的模型基本与其编录模型一致。

图2 818-2#采场采空区实体模型拟合图Fig.2 818-2# stope goaf solid model fitting diagram

3 数值模型的构建

3.1 本构模型的选择

根据研究材料的特点，在进行数值分析时，均以Mohr-Coulomb准则为依据，即选用弹塑性模型为其本构模型。

3.2 围岩物理力学参数及边界条件

根据该矿床的特点，选择矿体和围岩为主要力学介质进行分析。其强度参数的计算根据常用的岩体强度参数折减标准进行折算，内摩擦角φ按0.85折减，内聚力按1/10～1/7折减，弹性模量按2/3折减，得到某铀矿岩体强度力学参数见表1。

表1 某铀矿岩体力学参数Table 1 Mechanical parameters of rock mass of a uranium mine

(1)初始应力。因缺少矿山未做初始应力测试，初始应力采用岩石自重应力进行计算。

(2)计算边界条件。采用位移边界条件，即模型的表面施加位移约束条件，上边界施加根据上部岩体重量换算出的应力，根据地表三维模型和空区三维模型可知，地表距离数值模型上边界水平距离约30 m，施加在上边界的应力为0.9 MPa。

3.3 数值分析模型

利用“DIMINE to Flac3D”数据转换程序将某铀矿空区块体模型质心文件转换成FLAC3D网格数据文件，最终建立某铀矿空区FLAC3D数值计算模型，共生成网络节点531 441个，单元512 000个，如图3所示。

①～⑩为顶板位移监测点图3 818-2#采空区FLAC3D数值计算模型Fig.3 818-2# goaf FLAC3D numerical calculation model

4 数值模拟结果与分析

4.1 应力分析

4.1.1X剖面(X=8 695)分析

图4(a)和图4 (b)为818-2#采场X剖面上采场开挖前后围岩最大主应力分布图，图5 (a)和图5(b)为818-2#采场X剖面上采场开挖前后围岩最小主应力分布图。

图4 采场开挖前后X剖面最大主应力分布图Fig.4 Maximum principal stress distribution before and after quarry excavation of X profile

由图5(b)可知，采场开挖后，最小主应力发生变化，且有增大的趋势，根据开挖前后主应力分布图可知，最小压应力由-4.5 MPa增加到开挖后的-4.9 MPa，在空区顶底板处出现了最大值。最大主应力基本保持不变，最大值为-1.4 MPa，随着采场的开挖，在上下盘位置出现拉应力，应力值为0.5 MPa，低于围岩抗拉强度，空区围岩稳定。

图5 采场开挖前后X剖面最小主应力分布图Fig.5 Minimum principal stress distribution before and after quarry excavation of X profile

4.1.2Y剖面(Y=4 845)分析

图6(a)和图6(b)为818-2#采场Y剖面上采场开挖前后围岩最大主应力分布图，图7(a)和图7(b)为818-2#采场Y剖面上采场开挖前后围岩最小主应力分布图。

图6 采场开挖前后Y剖面最大主应力分布图Fig.6 Maximum principal stress distribution before and after quarry excavation of Y profile

图7 采场开挖前后Y剖面最小主应力分布Fig.7 Minimum principal stress distribution before and after quarry excavation of Y profile

由图6可知，采场开挖后，最大主应力发生较为明显的变化，其值呈逐步增大的趋势。根据图2可知，最小主应力由-4.56 MPa增加到开挖后的-4.67 MPa，且在空区顶底板处出现了最大值。最大主应力基本保持不变，最大值为-1.4 MPa，随着采场的开挖，在上下盘位置出现拉应力，应力值为0.17 MPa，低于围岩抗拉强度，空区围岩稳定。

4.2 位移分析

4.2.1X剖面(X=8 695)分析

由于空区的开挖，围岩的应力和位移发生变化，在围岩稳定性分析中，主要围绕应力和位移的分布变化特征展开研究，相较于应力变化，位移的变化能够更加直接的反映围岩稳定性状态。图 8给出了X剖面上空区形成后的水平和垂直向位移分布变化图，从图8中可见围岩位移分布变化呈现如下特点。

图8 采场开挖后X剖面位移分布图Fig.8 Displacement distribution map after quarry excavation of X profile

(1)818-2#采场回采后，最大水平位移出现在采场上下盘靠近顶底板位置，位移最大值为2 mm。

(2)818-2#采场回采后，最大垂直位移出现在采场顶底板位置，位移最大值为10 mm。

(3)垂直最大位移比水平最大位移大，这与矿山最大主应力方向为垂直方向相符。

(4)空区形成后，空区顶底板和两边位移最大为10 mm左右，空区基本处于稳定状态，但也应当注意顶板安全。

4.2.2Y剖面(Y=4 845)分析

图9分别给出了Y剖面上空区形成后的水平和垂直向位移分布变化图，可见围岩位移分布变化呈现如下特点。

图9 采场开挖后Y剖面位移分布图Fig.9 Displacement distribution map after quarry excavation of Y profile

(1)818-2#采场回采后，最大水平位移出现在采场上下盘靠近顶底板位置，位移最大值为3 mm。

(2)818-2#采场回采后，最大垂直位移出现在采场顶底板位置，位移最大值为9 mm。

(3)垂直最大位移比水平最大位移大，这与矿山最大主应力方向为垂直方向相符。

(4)空区形成后，空区顶底板和两边位移最大为9 mm左右，空区基本处于稳定状态，但也应当注意顶板安全。

4.2.3 采空区顶板位移分析

通过在采空区顶板设置监测点见图3(b)，监测在整个计算过程中顶板的位移随计算步骤的变化过程。图 10和图 11所示为顶板监测点最终位移图和特定开挖步骤顶板各监测点垂直位移曲线图。可以看出，最终顶板最大位移-6.32 mm，为6号监测点，在空区形成后顶板迅速下降，之后趋于稳定，说明形成采空区后顶板处于稳定状态。

图10 顶板监测点位移图Fig.10 Displacement map of roof monitoring points

图11 顶板监测点垂直位移图Fig.11 Vertical displacement curve of roof monitoring points

4.3 塑性区分析

图12(a)为空区形成后X剖面围岩塑性区分布情况，可以看出，随着采场的逐步回采，采空区围岩部分区域出现了塑性区，主要体现在采场下盘靠上部位置，塑性区面积很小，对整个空区稳定性影响不大，空区整体处于稳定状态。

图12(b)为空区形成后Y剖面围岩塑性区分布情况，可以看出，随着采场的逐步回采，采空区围岩部分区域出现了塑性区，主要体现在采场顶板和采场东部位置，塑性区面积很小，对整个空区稳定性影响不大，空区整体处于稳定状态。

图12 采场开挖后空区围岩塑性区分布图Fig.12 Distribution of plastic zone of surrounding rocks in the cavity after quarry excavation

5 结论

(1)运用BLSS-PE三维激光扫描和测量系统探测采空区，经矿用软件DIMINE构建采空区三维模型，实现了818-2#采场采空区三维实体模型的构建，为精确探测空区三维形态提供了新的有效途径。

(2)根据已构建的采空区三维模型运用FLAC3D进行数值计算，实现了采空区实测技术与数值模拟技术与数值模拟技术的耦合。

(3)根据模拟结果分析了采空区对周边围岩的稳定性影响。从应力变化特征、位移变化特征和塑性区变化特征等方面分析，818-2#采空区整体处于稳定状态，但为矿山安全考虑，应及时对采空区进行处理。