南芬露天矿滑坡机理及其动态过程

王长军，曹春晖，陶志刚

(1.本钢歪头山铁矿，辽宁 本溪 117005；2.深部岩土力学与地下工程国家重点实验室，北京 100083；3.中国矿业大学(北京)力学与建筑工程学院，北京 100083)

随着露天矿开采深度的增加，矿山边坡面临一系列地质灾害问题[1-2]，比如滑坡、泥石流、崩塌、沉降、台阶失稳等。滑坡是露天矿边坡常见的地质灾害，往往会造成严重的经济损失和人身财产损失。滑坡己经成为影响矿山安全生产和长久健康运营的重要灾害问题，迫切需要对露天矿滑坡的机理和动力学过程做更深入的研究。

国内外学者对不同的露天矿滑坡机理进行了分析研究，开挖扰动、冻融循环、降雨、地震等是露天矿滑坡发生的主要原因[3-6]。多数研究仅考虑其中一种或少数几种因素，滑坡机理缺少系统性综合分析。另外，对露天矿滑坡的动力学全过程分析相对较少。受前人研究的启发，本文以2019年1月4日发生在南芬露天矿下盘394～502 m台阶中部老滑坡体区域的滑坡为例，对其机理进行研究。结合现场勘查、地质资料，并采用光滑粒子流体动力学(SPH)方法结合Mohr-Coulomb弹塑性模型，建立三维数值计算模型，对露天矿滑坡进行动力学全过程分析，以了解其运动过程。SPH方法[7]最初用来研究天体物理的运行规律，并逐渐应用于流体力学、固体力学、冲击爆炸等领域。作为一种无网格纯拉格朗日方法，SPH方法在处理大变形问题如滑坡时具有独特优势，在岩土工程领域得到了广泛应用，并取得了一系列进展[8-9]。相对于有限元等方法，SPH方法由于无网格，可以避免大变形过程中产生的网格畸形等问题。而一些离散元方法，在处理像露天矿滑坡这类工程尺度的大变形问题时，对计算机内存等要求较高，计算效率面临挑战。SPH方法以其高效率和可靠性，可以很好地再现露天矿滑坡的全过程，并根据数值结果为露天矿边坡划定危险区域，制定防灾减灾措施，为露天矿安全生产和运营提供指导性建议。

1 滑坡背景

1.1 滑坡介绍

本溪市南芬露天铁矿是亚洲最大的单体露天矿山。近20年来，在特殊地形和长期矿山开采综合影响下，采场下帮边坡发生过多次滑坡，压矿1 000多万t，严重影响了矿山的正常开采。

2019年1月4日在南芬露天矿下盘394～502 m台阶中部老滑坡体区域发生滑坡，见图1。滑体对矿山安全作业带来了严重影响。由于及时预警，所幸未造成人身和设备损失，但滑落的碎石破坏了部分排水管路，阻挡了采矿车辆的正常出入，影响了正常的开采作业。本次滑坡体积大约为86 400 m3。露天矿边坡整体较陡，滑坡体滑动方向约为264°。滑坡体厚度小于10 m，属于浅层滑坡。滑坡源区上部宽度约为50 m，下部宽度约为140 m，上下垂直落差大约110 m。 滑坡传播区垂直落差大约120 m，底部宽度约为240 m。露天铁矿滑坡源区为老滑坡体堆积区，老滑坡发生后滑坡体整体呈松散状态挂靠于边坡上。老滑坡体堆积层主要为碎石土，整体稳定性较低，在爆破振动，开挖、降雨、冻融等因素下极易发生滑坡。

图1 滑坡特征Fig.1 Feature of the landslide

1.2 地质和气象水文资料

南芬矿区位于华北地台辽东台背斜营口-宽甸隆起的北缘太子河凹陷之中，新华夏系第二隆起与东西天山-阴山构造带的交汇处。地貌为单斜构造，主要由变质岩地层组成，属中高山侵蚀构造地貌。矿区内广泛发育太古界鞍山群，其次为元古界辽河群和震旦系地层及新生界第四系地层。南芬露天铁矿属于大型鞍山式沉积变质铁矿，含铁岩段主要赋存于太古界鞍山群，表现为单斜构造。

该区域水文网络不发育。冬天降雨不多，夏季降雨量增加，河水变得丰富。大气降水是地表水和地下水的主要来源。年平均降雨量、最大降雨量和最小降雨量分别为848～856 mm、1 259.8 mm和518.5 mm。汛期降雨量占全年降雨量的70%。 该区地处北温带季风气候区，四季分明。 年平均气温为7.7～8.2 ℃，最低气温为-24 ℃。 霜期(降雪期)一般为11月下旬至次年3月，最大冻结深度为1.49 m。

1.3 滑坡体土体成分

南芬露天矿下盘滑坡源主要为老滑坡堆积体，这些堆积的碎石土主要成分为砾石、块石、矿石碎屑、细砂和黏土。老滑坡堆积体碎石土的主要岩性为底盘角闪岩、绿泥片岩、绿泥角闪片岩等。碎石土黏度较低，颗粒强度高，水稳定性好。经过现场勘查，结合室内和室外试验，给出滑坡源松散堆积物的各级颗粒含量，见表1。

表1 不同级别颗粒组成情况Table 1 Composition of different grades of particles

2 滑坡机理分析

南芬露天矿年度剥岩量和铁矿石开采量巨大，经过长年的开采，逐渐形成了超高边坡。这些高边坡在自重、降雨、爆破振动、机械开挖等因素作用下，其稳定性逐渐降低，并成为矿区安全生产的重大隐患[10]。露天矿上盘边坡以逆倾向岩体为主，下盘边坡则以顺倾向岩体为主。在特殊地形和长期矿山开采综合影响下，采场下帮边坡形成了多处较大规模的滑坡体；其中，下盘394～502 m台阶中部老滑坡体区域，由于其为散体堆积体，稳定性较差，成为滑坡研究的重点区域。在2019年1月4日滑坡发生前，该区域曾发生明显的变形；其中，在502台阶出现长度约为6 m的裂缝，最宽处20～30 cm(图2(a))。主裂缝附近伴有数条宽度2～20 mm的细小裂缝，长度不等，这些裂缝使得降雨入渗变得容易。在降雨、振动及自重等作用下，这些裂纹逐渐扩展，使得降雨入渗进一步增强，边坡稳定性逐渐降低。另外，在滑坡发生前，502 m台阶滑移量为7～8 m(图2(b))。

“滑坡牛顿力空地监测预警系统”是深部岩土力学与地下工程国家重点实验室何满潮院士团队自主研发的滑坡监测预警系统。该系统集滑坡监测、预警、支护加固等技术为一体。牛顿力[11]的变化可以作为判断灾变发生的条件来对滑坡进行预警，可以得到较为准确的滑坡临滑信息。从本文滑坡发生前的牛顿力监测曲线(图2(c))变化可以看出，在2019年1月3日，牛顿力出现突变，指示滑坡可能即将发生；1月4日滑坡发生，验证了牛顿力监测系统的可靠性。

图2 滑坡前兆与滑坡牛顿力监测曲线Fig.2 Landslide precursor and Newton force monitoring curve of landslide

露天矿边坡的稳定性往往受到开挖岩石结构地质的显著影响，如层理、节理和断层，这些不连续面经过爆破振动、风化等活动逐渐破碎成更小的块体颗粒。这些破碎的碎石土一部分堆积在边坡台阶上，一部分成为老滑坡体的堆积体。老滑坡体堆积在露天矿台阶上，结构松散，稳定性差。

岩质边坡的开挖通常涉及爆破。露天矿爆破的基本要求是产生一个碎渣堆，以适应挖掘作业以及运输设备。在多层理、节理和高裂隙岩石中，爆炸气体喷入、楔入并扩展原有的裂缝。总体破碎程度受构造地质控制。在爆破过程中，受施工技术等多方面制约，台阶完整性和稳定性会随着爆破振动而降低。而老堆积体的松散岩体，受爆破振动影响，内部构造会受到扰动，并处在一个不断变化调整的过程。在爆破振动作用下，老堆积体内的碎石土会逐渐松散，其稳定性亦随之逐渐降低，当土体内部块石之间的互锁作用减弱，土体黏聚力和摩擦力降低到某一值时，滑坡就会发生。

坡脚开挖及台阶碎石土的清理，会降低露天矿边坡稳定性。在老滑坡体下部开挖或清理台阶碎石土时(清理浮石，图3)，上部土体失去一部分阻挡和支持作用，坡脚开挖处的部分临空土体受到扰动而滑落，并牵拉上部土体，当上部土体的自重产生的下滑力大于抵挡下滑作用的摩擦阻力时，滑坡随之发生。

图3 浮石清理Fig.3 Cleaning of float stone

降雨通常被认为是滑坡最关键的诱发因素[16]。老滑坡体内含有细颗粒物，这些细颗粒物可以吸收水分，部分颗粒遇水膨胀，使得土体的持水能力增强。老滑坡堆积体内的水分可以降低颗粒间的摩擦力。在冲刷和渗透作用下，细颗粒会在骨架中迁移并被雨水带走，粗颗粒物之间的间隙增大，从而改变原土体结构，而降雨入渗变得更加容易。当有持续降水时，土体内部的水压力不断增大，导致排土场土体抗剪强度和稳定性下降，并诱发滑坡的发生。此外，冻融也是导致滑坡发生一个原因。冻土层的膨胀会使土体局部隆起；当冻土层融化时，会引起土体收缩。这种膨胀和收缩导致排土场局部损伤、薄弱面和裂缝的增长和扩展。最后，坡面出现不均匀的崩塌和沉降。土体内部的某些冻土层在某些条件下成为软弱面或者潜在滑动面，在自重或外部荷载扰动下，滑坡可能会发生。综上所述，南芬露天矿排土场滑坡是爆破振动、坡脚开挖、降雨、冻融等综合作用下发生的。综合理解滑坡机理，对于制定露天矿防灾减灾措施具有重要意义。

3 滑坡动力学全过程分析

3.1 SPH理论和Mohr-Coulomb模型

SPH法[13]是将计算域离散成一定数量的粒子，每个粒子占据一定质量和体积，并携带密度、速度、加速度和压力、应力等物理参数。对于域内任意的函数f(x)，其核近似形式可以表示为式(1)[17]。

(1)

式中：h为平滑长度；W为平滑函数。

滑坡过程中，粒子的运动和质量守恒的控制方程可以写成式(2)和式(3)。

(2)

(3)

式中：ρ、vα、t、σαβ、xβ和fα分别为密度、速度、时间、应力张量、位置和体力；α和β为笛卡尔坐标方向。

老滑坡体为松散岩土颗粒，本研究采用弹塑性Mohr-Coulomb本构模型来研究其力学特性，其本构模型可以用三个不变量来表示，见式(4)。

f(I1,J2,J3)=Rmcq-ptanφ-c=0

(4)

其流动势函数g可以表示为式(5)[18]。

(5)

式中：I1为应力第一不变量；J2和J3分别为偏应力第二不变量和第三不变量；φ为材料摩擦角；c为黏聚力；p为等效压应力；q为Mises等效应力；为子午线偏心率，它定义了双曲函数接近渐近线的速率；c0为初始黏聚力屈服应力；ψ为在高围压下在p-Rmwq平面上测得的膨胀角；Rmc和Rmw为分别为和φ以及ψ有关的函数。通过引入应力张量和应变率张量的加法分解以及塑性屈服理论，并经过一系列计算(推导过程参考[19])，SPH形式的应力应变关系可以表示为式(6)。

(6)

3.2 滑坡动力学过程数值模拟

露天矿边坡因其规模较大，研究人员常常采用数值方法研究其稳定性。为了清楚地了解露天矿滑坡动力学过程，本文选取2019年1月4日发生在南芬露天矿下盘394～502 m台阶中部老滑坡体区域的滑坡为例，建立数值模型并进行全过程分析。边坡岩体相对比较完整，其黏聚力为20 740 kPa，摩擦角为37°。而附着在边坡表面的老滑坡堆积体为松散的碎石土，其黏聚力大大降低，而摩擦角也适当减小。南芬露天矿排土场滑坡土体形态、粒径和该露天矿边坡老滑坡堆积体亦较为接近。参考南芬露天矿排土场发生的滑坡土体的黏聚力和摩擦角，并根据中国科学院武汉岩土力学研究所进行的室内和现场试验，结合露天矿工程地质资料，并根据试错法获得使数值模拟结果和实际滑坡相符的力学参数(表2)。

表2 SPH滑坡模型材料属性Table 2 Material properties of the SPH model

露天矿滑坡SPH数值模型尺寸如图4所示。滑坡体上部宽50 m，下部宽140 m，上下垂直落差为105 m。整个模型尺寸为：长412 m，宽300 m，高288 m。图4中深色部分为滑坡源，在数值计算开始时自动转换成粒子，粒子总数为15 839。采用ABAQUS/Explicit计算，计算时长为100 s。

图4 SPH滑坡模型尺寸Fig.4 Model size of the landslide by the SPH method

模型边界条件施加如下：限制模型左右两侧Z方向的位移；限制模型前后两侧X方向位移；限制模型下方三个方向的位移。整个模型施加重力荷载。滑坡土体和边坡岩体直接设置通用接触，并设置摩擦系数为0.75，以模拟实际接触情况。

图5为露天矿滑坡的真实形态和SPH数值模拟得到的形态。由图5可知，数值模拟形态和真实形态比较吻合。松散岩土体在露天矿边坡台阶上向下滑动，滑动形态呈倒V字形。由数值结果可知，SPH模型可以很好地模拟露天矿滑坡，并可以根据模拟结果滑动危险区域范围，为露天矿边坡制定防灾减灾策略。

图5 真实滑坡与SPH模拟对比Fig.5 Comparison between real landslide and SPH simulation

为了更好地了解露天矿滑坡动力学全过程，选取t=0 s、10 s、20 s、30 s、50 s和100 s来分析不同时刻的滑坡形态，如图6所示。由图6可知，土体颗粒大概在t=20～30 s时滑动到最低处台阶上；从t=50 s起，滑坡形态不再有很大变化，表明滑坡逐渐趋于停止。另外，从图6可以清晰看出，粒子在滑坡运动过程中逐渐向两侧扩散，这使得滑坡的最终形态呈现近似三角形，这与实际情况相符。

为了了解滑坡过程中土体颗粒的位移信息，露天矿滑坡在t=0 s、10 s、20 s、30 s、50 s和100 s的位移云图如图7所示。图7显示位移在中间线位置处为最大，并呈现出左右对称的分布形式。

图7 露天矿滑坡的总位移云图Fig.7 Total displacement(U-Magnitude) contour of the landslides for the full high waste dump

图8给出了露天矿滑坡从t=1 s到t=100 s的速度云图。由图8可知，粒子经历了先加速后减速至零的过程。从图8(d)和8(e)对比发现，粒子的速度云图形态基本一致，但是图8(e)上显示带有速度的粒子个数明显比图8(d)的数量少。其原因是，粒子在露天矿边坡上向下滑动时，一部分粒子存留在了台阶上，这使得滑动层变薄，进而使得图8(e)中的显示带有速度的粒子数目明显减少。在实际露天矿滑坡过程中，土体一边下滑一边在各级台阶上堆积。数值模拟结果和实际滑坡表现出了很好的一致性。

图8 露天矿滑坡的速度云图Fig.8 Velocity (V-Magnitude) contour of the landslides for the full high waste dump

总的来说，SPH方法凭借其自身的高计算效率和占用较低计算机内存等优势，对于描述露天矿滑坡等大变形问题是很适用的。在本文中，结合SPH方法与Mohr-Coulomb弹塑性理论，将南芬露天矿边坡老滑坡松散体离散成一系列粒子，这些粒子携带了代表散体颗粒的位置、密度和速度等信息，可以完整地描述老滑坡体在自重、降雨、冻融及振动等作用下的动态滑动过程。该方法可以获得任意时刻、任意位置粒子的速度、位移、能量耗散等信息，借此可以全面了解露天矿边坡滑坡的动态全过程。数值结果也可以帮助更好地了解滑坡机理。SPH方法属于无网格法，可以避免出现网格畸形而导致的计算不收敛问题。除了计算效率高之外，其可靠性亦能得到保障，如图6所示，SPH计算得到的滑坡形态和范围非常接近实际滑坡形态。另外，滑动距离和速度等信息同样可以按要求在数值模拟结果中予以显示(图7和图8)。数值计算结果可以为露天矿边坡滑坡灾害提供一些防治减灾策略，并可以根据数值结果中滑坡形态和范围为滑坡灾害划定范围，提前制定防治减灾措施，比如对边坡进行加固支护等，以保障露天矿的安全运营。不过，在SPH方法中，露天矿边坡老滑坡体被视为单一物质，未考虑老滑坡体内的弱面、空隙和水压力等影响。在以后的研究中，如果考虑流固耦合等问题，可能会使得计算精度进一步提高，但计算效率可能会受到一定程度的影响。因次，在SPH方法中更加真实地描述露天矿边坡物质组成，以及如何保障高计算效率，将是下一步的研究重点。

4 结 论

1) 本研究以露天矿下盘394～502 m台阶中部老滑坡体区域的滑坡为例，对其机理进行了研究。老堆积体的自身稳定性较低，开挖或爆破振动、降雨和冻融作用是导致滑坡的最主要原因。

2) 采用SPH方法结合弹塑性Mohr-Coulomb模型对露天矿滑坡进行动力学全过程模拟，数值结果在滑动距离、形态和范围等方面与实际滑坡表现出较好的一致性。SPH方法可以用来描述露天矿滑坡动力学过程。

3) SPH数值结果可以为露天矿滑坡划定危险范围，为矿山安全运营提供防治减灾的指导性建议。