滨海矿山深部开采水力调查及稳定性分析

陈倩男 郭奇峰 万思达 翟 济 杜伟嘉 李 钊

(1.北京科技大学土木与环境工程学院，北京100083;2.金属矿山高效开采与安全教育部重点实验室，北京100083)

海洋蕴藏着丰富的海底矿产资源，包括海底表层沉积物和海底岩层中贮藏的矿产。目前开发利用最为广泛的是海洋石油和天然气，在整个海底矿产资源已开发产值中的占比超过90%。陆地矿产资源经过多年开采，浅部资源枯竭，国内外很多矿山进入深部开采阶段［1-2］，并且诸多近海矿山开始进行海底矿产资源的开采［3-4］，其中日本的海底煤矿开采量达到总产量的30%，海底开采将成为未来矿产资源开发的发展方向。

海底基岩岩体存在大量节理裂隙，水和应力通过节理裂隙进行耦合作用，对海底资源开采过程的安全具有重要影响。因此，在海底资源开采必须深入调查水、节理裂隙和地应力的分布特征，研究水和应力的耦合作用机理。诸多学者开展了大量针对裂隙岩体中水和应力等因素耦合作用的研究，主要工作集中在理论模型研究，试验分析和数值模拟应用3 个方面。王启平［5］、马海涛等［6］基于理论研究和试验分析，采用数值模拟软件，对强夯地基和露天采矿边坡稳定进行了水力耦合分析，验证了水力耦合在工程应用中的必要性和适用性。本课题以滨海矿山三山岛金矿深井开采为背景开展水力研究，为水力耦合作用下的海底资源安全开采提供指导和保障。

1 水和地应力调查

1.1 地下水分布

矿体及上覆顶板中的基岩裂隙含水是三山岛金矿深部开采区域直接充水层。根据含水层与隔水控矿构造F1的相对位置关系，直接影响矿坑充水的基岩裂隙含水带可划分为F3断裂含水带、F1上盘裂隙含水带和F1下盘构造裂隙含水带。

上述含水带中，F1下盘构造裂隙含水带分布广泛，直接影响深部开采和井巷工程施工。根据构造发育特征和补给条件的不同，划分为F1下盘构造裂隙富水区和F1下盘构造裂隙含水区。F1下盘构造裂隙富水区内，北西向导水裂隙发育，且裂隙规模和延伸较大，含水性较强;北西端有海水补给，东南方向与F3断裂构造含水带沟通，有区域卤水补给。F1下盘构造裂隙含水区逐渐远离F3，富水性较弱，北西向含水裂隙密度和规模逐渐变小，北西部有海水越流补给，并有少量的第四系孔隙水渗透补给。

1.2 地应力场分布

矿山深部开采面临“三高一扰动”的问题，其中地应力是一个重要影响因素，其关系到井巷工程和开采工程的布置，直接影响深部开采的安全性。三山岛金矿深部开采中，高应力显现明显，局部地区甚至出现岩爆现象。为保证现有开采工程的安全，为更深层次的开拓和开采设计提供安全指导和保障，用套孔应力解除法在深部－500 m 以下6 个水平开展9 个测点的地应力测量工作，结合2002 年在浅部3 水平4个测点的地应力测量结果［7］，共获取三山岛金矿13个测点的地应力测量数据，通过回归分析，绘制三山岛金矿深部主应力随深度变化曲线如图1 所示。

图1 主应力随深度变化曲线Fig.1 Variation curve of principal stresses with depth

矿区最大水平主应力、最小水平主应力和垂直主应力的回归方程为

2 岩体渗透性

2.1 节理裂隙分布

对中段运输平巷、采场、联络巷等进行现场节理裂隙调查，记录节理倾向、倾角、走向、节理间距、充填物及渗水特征，共获取854 组节理记录，其倾向及倾角分布统计如表1 所示，绘制节理倾向玫瑰花图和走向玫瑰花图如图2 和图3 所示。

表1 节理分布特征Table 1 Joints distribution

图2 节理倾向玫瑰花图Fig.2 Rose diagram of joints＇strike

图3 节理走向玫瑰花图Fig.3 Rose diagram of joints＇dip

统计结果显示，深部节理发育具有5 组优势倾向。第1 组优势倾向为20° ～30°，平均倾角62°;第2组优势倾向为110° ～150°，平均倾角57°; 第3 组优势倾向为190° ～200°，平均倾角66°; 第4 组优势倾向为270° ～280°，平均倾角54°; 第5 组优势倾向为290° ～320°，平均倾角55°。走向具有5°、20° ～40°、270° ～280°3 组优势方位，倾角多分布在40° ～90°。深部节理裂隙的倾角均偏大，平均倾角超过55°。

2.2 岩体渗透特性

2.2.1 几何测量计算

根据Snow［8］提出的渗透张量几何计算方法，假定每个裂隙内没有填充物，均无限延伸且相互连通，但各裂隙的渗流互不影响。以地理N 向为x 轴，地理E 向为y 轴，则各裂隙渗透系数张量Ki的表达式为

式中，ρ 为地下水密度，kg/m3; g 为重力加速度，m/s2; bi为第i 组裂隙宽度，m; μ 为地下水动力黏滞系数，kg/(m2·s); Kxxi= 1 － sin2αicos2βi; Kxyi= Kyxi=－sin2αicosβisinβi; Kzxi= Kxzi=－sinαicosαicosβi; Kyyi=1 －sin2αisin2βi; Kzzi= sin2αi;αi为第i 组裂隙倾向，0≤αi≤360°; βi为第i 组裂隙倾角，0≤βi≤90°。

作为二维张量，渗透系数张量的特征方程可表达为

渗透张量的3 个主值是上述特征方程( 5) 的3个特征根A1、A2和A3，则综合渗透系数K［9-11］为

选取多个有代表性的上盘、矿体、下盘岩体典型表面裂隙特征调查窗口，详细统计窗口内的节理裂隙产状和分布特征。依据上述渗透系数几何计算原理，采用MATLAB 软件计算测量区域内上盘、矿体、下盘的综合渗透系数，结果如表3 所示。

表2 渗透系数Table 3 Permeability coefficients m/s

2.2.2 渗透系数修正

Snow 几何量测法在进行渗透系数计算时，假定所有节理裂隙在裂隙网络交叉空间相互连通，但现实岩体中节理的延伸范围有限，各组节理并未呈现相互交切的形态。为了更准确地掌握深部矿岩渗透特性，依据节理裂隙产状、间距等分布特征，采用Monte －Carlo 随机模拟法［12-13］对上盘、矿体、下盘3 类区域节理裂隙的渗透系数进行修正。模拟结果显示，上盘岩体的连通率为84.76%，矿体连通率为88.49%，下盘岩体连通率为85.13%。张有天［14］提出采用连通系数即连通率对理想条件下计算出的渗透系数进行修正:

式中，k 为修正后的渗透系数; k0为计算值; η 为连通率。将几何计算得出的渗透系数乘以节理裂隙连通率即获得修正后的上盘渗透系数为6.236 1 ×10－7m/s，矿体渗透系数为8.349 2 ×10－7m/s，下盘渗透系数为6.743 3 ×10－7m/s。

3 深部开采水力耦合分析

3.1 计算模型

采用FLAC3D数值模拟软件建立深部开采水力耦合数值模型，模型标高自－430 m 至－630 m。坐标系以矿体走向为x 轴，矿体厚度方向为y 轴，铅垂方向为z 轴。计算模型在x 方向上的长度为400 m，y方向上的长度为300 m，z 方向上的高度为200 m，共划分58 310 个单元，63 900 个节点。对中段F3断层北侧矿体进行房柱交替式盘区上向分层充填模拟开采，盘区长度选取80 m，盘区沿矿体走向共布置8 个采场进行矿房矿柱交替上升式开采，每个采场宽10 m，采场长度为矿体厚度，中段开采共分6 步进行，如图4 所示的1 ～6。

图4 开采步骤Fig.4 Mining steps

3.2 边界条件

根据三山岛金矿深部详查报告中地下水位长期观察和水位降低值的研究结果［15］，在模型上部和底部施加水头边界。依据深部地应力场分布规律计算深部最大主应力、最小主应力和垂直主应力值，对模型施加梯度应力边界。

3.3 力学参数

模型划分上盘、矿体、下盘以及F3断层4 大组，F3断层倾角近90°，垂直矿体走向，矿体倾角45°，开采后直接充填。依据物理力学试验数据和渗透性计算数据，选取水力耦合数值模拟采用的力学参数和渗透系数如表3 所示。

表3 物理力学参数Table 3 Physico-mechanical parameters

3.4 模拟结果分析

3.4.1 采场位移变化

盘区内最大竖向位移均发生在新开采区域的中心位置，竖向位移朝开采区域集中，由图5 所示开采至第二步的采场矿体和充填体竖向位移等值线图可见，各开采步骤竖向位移出现明显分区，即1 到4 号采场为1 组，5 ～8 号为1 组，在4 和5 采场交界带竖向位移值相对较小，呈现出对称性。开采过程中，采场顶板沉降逐渐增大，变化较为稳定，无明显突变发生，在各采场开采接近顶柱时，顶板沉降达到最大值。

图5 开采第二步的竖向位移等值线Fig.5 Vertical displacement isoline in 2nd mining step

为研究开采接近顶柱时的采场安全性，在中段每个采场上方对应的顶柱设置监测点共计8 个。各监测点在中段开采过程中的竖向位移因房柱交替开采呈现阶段跳跃式增长，当开采高度超过中段高度一半之后，其沉降增幅明显，开采结束后2 号监测点的最大沉降达到9.276 cm，存在安全隐患。塑性区云图也显示各采场回采接近顶柱时出现塑性区贯通，应加大支护密度和管理力度，关注采场顶板的安全稳定性。

3.4.2 采场应力变化

在6 个开采步骤过程中最大主应力值分别为50.073、60.024、52.398、67.959，52.603，43.048 MPa，具有代表性的开采第2 步的最大主应力等值线如图6 所示。模拟结果显示在整个盘区开采过程中，二步采引起的应力集中要大于一步采，中间和两翼的最大主应力值较高，二步采阶段采场稳定性相对较差。最大主应力的极值主要出现在1 号采场围岩部位，5 号采场顶板应力集中显现也较为明显，开采过程中需要对上述2 个区域进行重点支护和加固，确保生产安全。

图6 开采第2 步的最大主应力等值线Fig.6 Maximum principal stress isoline in 2nd mining step

各个开采步骤过程中均存在拉应力区，极值分别为2.656、1.532、1.484、1.515、1.362、1.225 MPa，拉应力主要集中在一步采和二步采的交替部位，其分布区域具有一定的时空变化规律，即一步采阶段主要存在于2 和6 号采场，二步采阶段主要存在于3 和7 号采场。针对上述拉应力区的分布，回采过程中应加大支护力度，保证充填体的强度和充填接顶效果。

4 结 语

(1) 针对海底资源深井开采面临的富水和高应力问题，对三山岛金矿地下水和地应力进行调查和测量，掌握了地下水的分布，建立了矿区地应力场模型，为矿山开采和优化提供了基础资料。

(2) 开展深部采场和巷道围岩节理裂隙调查，统计深部岩体节理裂隙分布特征，采用几何测量计算和结构面随机模拟修正的方法研究了深部矿岩体的渗透特性。

(3) 建立水力耦合数值模型，对深部矿体房柱交替式上向分层充填开采进行动态模拟。分析了分步开采过程中采场位移、应力随开采推进的时空变化规律，对应力集中、拉应力破坏等危险区域进行了时空分析，确定了开采过程中需要重点支护的区域和时间，为海底资源水力耦合作用下的安全开采提供了保障。

［1］ 蔡美峰. 金属矿山采矿设计优化与地压控制——理论与实践［M］北京:科学出版社，2001.

Cai Meifeng.Optimization of Mining Design and Control of Ground Pressure in Metal Mines: Theory and Practice［M］. Beijing: Science Press，2001.

［2］ 王春来，吴爱祥，刘晓辉，等. 深井开采微震活动容量维Df 变化特征［J］.北京科技大学学报，2010，32(11) :1379-1382.

Wang Chunlai，Wu Aixiang，Liu Xiaohui，et al.Variation characteristics of capacity dimension Df with microseismicity in deep mining［J］.Journal of University of Science and Technology Beijing，2010，32(11) :1379-1382.

［3］ 乔卫国，吕言新，魏烈昌，等. 三山岛金矿海底高效开采工业试验［J］.金属矿山，2011(7) :43-46.

Qiao Weiguo，Lu Yanxin，Wei Liechang，et al.Industrial test of highefficient submarine exploitation in Sanshandao Gold Mine［J］.Metal Mine，2011(7) :43-46.

［4］ 彭 康，李夕兵，赵国彦，等. 岩层微扰上向分层充填法在海底开采中的应用［J］.金属矿山，2010(4) :50-52.

Peng Kang，Li Xibing，Zhao Guoyan，et al.Application of the upward slicing and filling with perturbation in undersea exploration［J］.Metal Mine，2010(4) :50-52.

［5］ 王启平，谢能刚，史小路.基于强夯大变形的地基流固动力耦合分析［J］.北京科技大学学报，2004，26(4) :345-348.

Wang Qiping，Xie Nenggang，Shi Xiaolu. Fluid-solid dynamic coupling analysis of foundation soil during large deformation by dynamic consolidation［J］. Journal of University of Science and Technology Beijing，2004，26(4) :345-348.

［6］ 马海涛，李迎佳，王云海.港里铁矿水下开采流固耦合数值模拟与风险分析［J］.金属矿山，2014(2) :36-39.

Ma Haitao，Li Yingjia，Wang Yunhai.Fluid-solid coupling numerical simulation and risk analysis of underwater mining of Gangli Iron Mine［J］.Metal Mine，2014(2) :36-39.

［7］ 乔 兰，欧阳振华，来兴平，等. 三山岛金矿采空区地应力测量及其结果分析［J］.北京科技大学学报，2004，26(6) :569-571.

Qiao Lan，Ouyang Zhenhua，Lai Xingping，et al.In-situ stress measurement and its result analysis in Sanshandao Gold Mine of China［J］.Journal of University of Science and Technology Beijing，2004，26(6) :569-571.

［8］ Snow D T. Anisotropic permeability of fractured media［J］. Water Resources Research，1969，5(6) :1273-1289.

［9］ 仵彦卿.岩体水力学基础( 七) ——岩体水力学参数的确定方法［J］.水文地质工程地质，1998，25(2) :42-48.

Chu Yanqing. Basic of rock mass hydraulics ( 7) : Determination methods of hydraulic parameters of rock mass［J］.Hydrogeology and Engineering Geology，1998，25(2) :42-48.

［10］ 黄克智，薛明德，陆明万. 张量分析［M］. 北京: 清华大学出版社，2005.

Huang Kezhi，Xue Mingde，Lu Mingwan.Tensor Analysis［M］.Beijing:Tsinghua University Press，2005.

［11］ 周志芳，王锦国.裂隙介质水动力学［M］.北京: 中国水利水电出版社，2004.

Zhou Zhifang，Wang Jinguo.Theory on Dynamics of Fluids in Fractured Medium［M］.Beijing:China Water ＆ Power Press，2004.

［12］ 方 涛，柴军瑞，胡海浪，等.Monte Carlo 方法在岩体裂隙结构面模拟中的应用［J］.露天采矿技术，2007(1) :7-9.

Fang Tao，Chai Junrui，Hu Hailang，et al. Application of Monte Carlo method to simulating the fracture interconnectivity of rock mass［J］.Opencast Mining Technology，2007(1) :7-9.

［13］ 杜景灿，陈祖煜，弥宏亮，等. 三维条件下应用遗传算法与Monte－Carlo 法确定节理岩体的综合抗剪强度［J］. 岩石力学与工程学报2004，23(13) :2157-2163.

Du Jingcan，Chen Zuyi，Mi Hongliang，et al.Determination of comprehensive shear strength for jointed rock massed in 3D condition using genetic algorithm and Monte-Carlo method［J］. Chinese Journal of Rock Mechanics and Engineering，2004，23(13) :2157-2163.

［14］ 张有天.岩石水力学与工程［M］.北京: 中国水利水电出版社，2005.

Zhang Youtian. Rock Hydraulics and Engineering［M］. Beijing:China Water ＆ Power Press，2005.

［15］ 三山岛金矿. 山东省莱州市三山岛矿区深部金矿详查报告［R］.济南:山东黄金集团有限公司，2007.

Sanshandao Gold Mine.Detailed Deep Survey Report of Sanshandao Gold Mine in Laizhou，Shandong Province［R］. Ji＇nan: Shandong Gold Group Co.，Ltd.，2007.