0.   引言

缓倾层状结构的自然边坡破坏通常以局部塌陷变形为主，不易发生重大地质灾害。然而对于矿区的缓倾层状边坡，在采空区变形与降雨的共同作用下，使采空区产生冒落破坏，上覆岩体变形松动，极大程度降低了边坡的稳定性，进而诱发滑坡等自然灾害，造成严重的人员伤亡与财产损失。1980年宜昌盐池河滑坡摧毁磷矿矿区，导致284人死亡和2500万元财产损失^[1]；1996年5月31日至6月3日，元阳县老金山金矿群采区接连发生两次滑坡，最终导致327人死亡，直接经济损失1.4亿元^[2]；2009年6月5日重庆武隆鸡尾山铁矿开采区发生滑坡-碎屑流特大灾害，造成包含27名矿工在内的74人遇难，8人受伤^[3]；近年的贵州水城尖山营滑坡、抚顺西露天矿南帮滑坡也均属于此类矿区地质灾害，破坏大量基础设施并掩埋道路，直接威胁人民生命财产安全^{[4 − 5]}。

采空区诱发山体滑坡是较为常见的一种地质灾害致灾模式，目前已有一定研究基础^[6]。王玉川等^[7]针对缓倾煤层采空区上覆山体滑坡的形成机制进行了研究，发现采空区的冒落作用是诱发此类灾害的主要原因。赵建军等^{[8 − 9]}通过物理模型试验及基于FLAC^3D的数值模拟方法探究了此类型岩体破坏进而诱发地质灾害的过程，进一步阐释了此类地质条件下滑坡诱发的机理。此外，节理发育会产生性质较差的结构面对岩土体进行切割，极大程度影响了边坡的稳定性，加速了岩土体破坏过程中的解体和破碎。目前已有研究表明，滑体的碰撞破碎解体会对滑坡的运动特性与致灾范围均产生较大影响^{[10 − 12]}。

鉴于节理与下伏采空区冒落变形的作用会对边坡稳定性产生极大影响，研究此类含节理及下伏采空区的失稳破坏过程，揭示滑坡破碎体堆积特征对相应的灾害预防工作将具有重要的指导意义。本文在对贵州省都匀市缓倾层状结构下煤矿采空的马达岭滑坡进行野外调查的基础上，在 PFC 中建立相关的地质模型，研究了马达岭滑坡在不同阶段的破碎特征，分析了破碎体数量的变化规律。研究结果对此类型的滑坡灾害防治提供了理论依据及科学指导。

1.   工程地质概况

马达岭滑坡位于贵州省都匀市江洲镇富溪村，所在地属溶蚀侵蚀低中山地形地貌。研究区域内斜坡为一典型的缓倾斜横向斜坡。滑源区位于山脊的东南侧，滑体前缘高程1400 m，后缘高程约1550 m，高差150 m。滑坡后壁陡立，前缘坡度25°，马达岭滑坡的滑后影像如图1所示。滑坡主滑方向南西8°，平面面积3 × 10⁴ m²，滑面深度40～60 m，体积约1.9 × 10⁶ m³。

图  1  马达岭滑坡碎屑流滑后无人机航摄影像

Figure  1.  Aerial image of Madaling landslide-debris flow

下载: 全尺寸图片 幻灯片

研究区域内岩层据现场调查从古至今主要包括：上泥盆统者王组（D₃z），下石炭统汤粑沟组（C₁t）、祥摆组（C₁x），第四系冲洪积（${\mathrm{Qh}}^{al+pl} $）等物质。岩层为灰、浅灰色薄至厚层状细粒石英砂岩间夹暗灰色薄至中厚层状泥质粉砂岩、黑色炭质泥岩、含炭泥质粉砂岩。岩层近南北走向，其中石英砂岩、白云岩为硬岩，泥质、炭质页岩为软岩，其间有多层煤层呈软硬互层状结构。其滑坡详细地质剖面图如图2所示。

图  2  马达岭滑坡碎屑流地质剖面图

Figure  2.  Geological cross-section profile of Madaling landslide-debris flow

下载: 全尺寸图片 幻灯片

就构造特征而言，根据节理统计结果岩体内主要发育两组近垂直的节理，滑源区主要包括两组产状分别为120°∠82°、330°∠87°的节理（图3）。两组节理较发育，且近乎直立，与坡面呈小角度相交，以30～50 cm不等的间距分布，有泥质填充，且节理面较为平整，将滑源区内岩体切割破碎为块体，极大地影响了边坡的稳定性。斜坡距坡顶200 m深度范围内发育3层可采煤层，且大量开采形成采空区，厚度1～2 m。

图  3  马达岭滑坡源区节理发育特征

Figure  3.  Joints development characteristics of Madaling landslide-debris flow

下载: 全尺寸图片 幻灯片

结合现场水文地质条件勘察及无人机航摄影像发现，采空区冒落作用下，上部岩体应力分布发生改变，应力平衡破坏，对滑源区内岩体形成张拉应力，孕生张拉裂隙，加之节理发育，岩体被裂隙、节理、岩层面切割强烈破碎。随着强降雨的发生，雨水渗入后缘裂缝，孔隙压力使裂缝不断增大，推动坡体向临空面变形，最终使得坡体滑动面抗剪强度无法抵抗下滑力，坡体失稳破坏。

织金县“6•8”白岩崩塌、金沙县“6•20”安洛乡崩塌等此类具有相似地质条件、采动降雨等多因素耦合诱发的地质灾害多发于贵州省六盘水市、毕节市，变形破坏过程与马达岭滑坡具有的相似性。

2.   离散元方法

2.1   细观参数标定

PFC是一种常用的离散单元方法程序，模拟中主要通过颗粒间内力和位移的不断更新来模拟介质体的变形，解决了常规数值方法分析大变形的局限性，能够更好地揭示微观及细观介质材料的破坏机理。计算中颗粒之间的接触可以选用不同作用定律的接触模型来施加。平行黏结模型是 PFC中一种较为常用的接触模型，该黏结可以在不同实体之间传递力和力矩，这种力和力矩作用在两个接触块上，与胶结材料接触周围的最大正应力和剪应力有关。当应力超过黏结强度时，平行键断裂，颗粒之间的接触转化为线性模型。

研究过程中选取研究区域内砂岩，经已有室内岩石力学试验，测定其物理力学参数^[13]，结果如表1所示。而在PFC中，宏观物理力学参数与颗粒的细观参数一般不能直接联系，但两者可以通过相应的数值模拟试验来建立联系。采用 PFC 内置参数标定程序来模拟室内岩石力学试验，对模型材料颗粒的参数进行标定。

表  1  马达岭砂岩物理力学性质室内试验结果

Table  1.  Laboratory experiment results of mechanical properties of sandstone at Madaling

指标	σt/MPa	σc/MPa	E/GPa	μ	c/MPa	φ/（°）	ρ/（g·cm−3）
数值	0.91	80	15.59	0.19	28.02	42.08	2.63

下载: 导出CSV 

| 显示表格

经过不断试算，PFC^3D数值模型中，颗粒半径取0.03～0.06 m，共生成3095个颗粒（图4a）。利用 PFC 程序内置的伺服机制控制试样压缩过程中的围压，进行多次试算，在不同围压作用下获得岩体的强度参数。获取模拟过程中的应力应变数据，绘制数值试样在围压分别为30，50 ，70 MPa 时的应力应变曲线。根据Mohr - Coulomb理论，通过绘制摩尔破坏包络线得到土体的强度参数，模拟得到φ = 41.49°，c = 28.48 MPa，与室内试验所得砂岩试样的强度参数接近，见图4（b）。最后基于相关处理方法对颗粒的粒径及参数进行相应的转换，获得的模型需要的细观参数如表2所示。

图  4  马达岭滑坡颗粒流数值模拟参数标定

Figure  4.  Parameter calibration of granular flow numerical simulation for Madaling landslide

下载: 全尺寸图片 幻灯片

表  2  PFC标定数值模拟微观参数

Table  2.  Micromechanical parameters for numerical simulation calibrated from PFC

材料参数	取值		材料参数	取值
颗粒密度（ρ）/（kg·m−3）	2630		胶结法向切向刚度比（$\kappa_b^* $）	1.0
颗粒有效模量（E*）/（N·m−2）	1.7e10		胶结抗拉强度（pb_ten）/（N·m−2）	7.2e8
法向切向刚度比（κ*）	1.0		胶结黏聚力（pb_coh ）/（N·m−2）	3.5e8
胶结有效模量（$E_b^* $）/（N·m−2）	3.9e7		胶结内摩擦角（pb_fa）/（°）	40

下载: 导出CSV 

| 显示表格

2.2   滑坡模型建立

数值模拟过程中的滑坡模型主要由滑面模型及滑体模型构成，两部分模型通过无人机航摄及现场调查结合的方式建模生成。无人机进行摄影控制航向重叠度达到90%以上及旁向重叠度达到80%以上，生成三维点云并进行三角剖分形成高精度数字表面模型（digital surface model, DSM），作为数值模拟滑面模型。滑体模型通过滑前历史地形地貌数据与滑后地形对比，在建模过程中通过布尔交集运算建立。通过现场调查发现，滑源区岩体内主要发育两组近垂直的节理，通过在离散元软件中对滑源区模型进行颗粒填充模型平衡后采用颗粒删除法来模拟。完成滑体内节理的制备后，使模型在静力作用下再次平衡，完成本研究所需的最终滑体模型，模型共计39459个颗粒。

3.   动力破碎模拟与分析

为了全面分析马达岭滑坡运动全过程中的破碎解体情况，本文从滑体破坏过程以及运动破碎过程两个阶段对其进行分析。颗粒离散元中，对破碎体（fragment）的定义是具有连接属性的一系列颗粒、簇或者墙体的集合体。在下述的分析中，本文采用PFC中fragment这一元素对不同时刻斜坡的破坏情况进行讨论。

3.1   滑体失稳过程

图5为滑体失稳过程中的颗粒破碎分布情况，图中不同颜色的碎片为PFC内置算法判定的破碎体。可以看出，T = 2 s时，斜坡内的破碎体主要分布在滑体前缘两侧以及后缘，表现为局部的小规模开裂破坏，滑动面尚未贯通，斜坡整体仍保持相对稳定。随后滑体后缘拉裂处产生粒径5 m以上的块状破碎体，滑体前缘深处岩体沿节理开裂并向临空面延伸，贯通后岩体沿节理形成层状破碎体，滑面进一步向后缘拉裂处发育扩展。在 T = 10 s时，后缘拉裂与前缘开裂扩展裂隙完全贯通形成滑面，滑体开始运动，滑面出露，滑体前缘下部层状破碎体在冒落作用下进一步受压破碎，产生大量粒径约为 1 m 的破碎体向主滑方向两侧运动散落。

图  5  PFC^3D模拟2～10 s马达岭滑坡碎屑流滑体破坏破碎过程 （时间间隔ΔT = 2 s）

Figure  5.  Simulation of fragmentation with PFC^3D from 2 to 10 s (ΔT = 2 s)

下载: 全尺寸图片 幻灯片

不难看出，节理在滑动面贯通、滑体失稳的过程中提供了不利条件，并将滑体切割为层状破碎体。下层岩体在上覆岩体的重力与冒落作用下受压进一步解体，使得滑体在失稳过程中就产生了较强烈的破碎，为滑坡碎屑流的远程运动提供了基础。

3.2   滑体运动过程

在滑体失稳后续过程中，随着滑动面出露，滑体主体产生较大位移并在重力作用下沿山体沟谷方向运动，伴随大量粒径约为1 m的破碎体呈扇形向垂直主滑方向发散。

图6为该类型滑坡-碎屑流的运动堆积过程。当T = 20 s时，滑体整体解体，后缘拉裂与前缘裂隙贯通加之节理切割形成粒径为10～20 m的大破碎体，前缘下部岩体由于重力与下伏采空区冒落作用被挤压强烈破碎产生粒径为1～4 m的小破碎体，颗粒分布形成下小上大的整体趋势，形成滑坡-碎屑流沿滑动面进行滑动。此时滑体内部裂隙已经完全贯通，节理不再对滑体的开裂破碎产生影响。

图  6  PFC^3D模拟0～180 s马达岭滑坡碎屑流滑体破碎堆积过程（时间间隔ΔT = 20 s）

Figure  6.  Simulation of fragmentation accumulation with PFC^3D from 0 to 180 s (ΔT = 20 s)

下载: 全尺寸图片 幻灯片

随后，滑体下部受挤压区域进一步破碎产生细小破碎体向两侧蔓延，前缘细小颗粒部分沿沟谷滑动（T = 40 s），在T = 60 s时刻滑体大部分颗粒离开源区，体积较大破碎体以较慢速度运动，仍主要分布在源区；较小颗粒以较快速度冲入沟谷中继续向前运动。由图6可知，从T = 80 s开始，已破碎的大量滑体在滑至坡脚停止运动后，其在地表的运动过程中仍有破碎解体现象发生，这一定程度上反映了此类斜坡滑体在滑动至坡脚过程中的冲击刮铲效应。最终滑体冲毁下方道路及农田，主要堆积在沟谷出口位置，其中堆积体破碎体体积大小不一，通过现场调查发现的实际破碎体堆积情况较为相似（图7）。

图  7  堆积区与PFC^3D数值模拟结果对比

注：a为前缘已被清理的堆积区；b为离散元模拟堆积区范围；c为堆积区巨石；d为离散元模拟产生的相似形状破碎体。

Figure  7.  Comparison between deposition area and simulation results with PFC^3D

下载: 全尺寸图片 幻灯片

从整个堆积体内不同大小破碎体的分布来看，沿主滑方向，尺寸在5 m以上的大粒径破碎体在堆积区内均有分布，而尺寸在1～2 m小粒径破碎体主要分布在堆积体前缘及下部。而在垂直主滑方向的剖面上，不难发现岩土体颗粒在沟谷内滑动的过程中在纵剖面上形成小颗粒在下大颗粒在上的反粒序分布，与现场所观察到的情况基本一致（图8）。已有学者研究表明^{[14 − 15]}，对于完全破碎的滑体，底部的颗粒因频繁碰撞而剧烈混合，可作为润滑层，以降低滑坡碎屑流的整体有效摩擦，从而促进其远程运动。本文对滑坡运动过程、最终破碎体信息进行了研究。

图  8  堆积体内存在的反粒序现象

Figure  8.  Reverse grading phenomenon within deposits

下载: 全尺寸图片 幻灯片

3.3   破碎体信息分析

3.3.1   破碎体数量变化情况

模型计算过程中，通过PFC内置History命令监测了滑体内破碎体的数目变化，如图9所示。从图中可知，在滑体失稳过程中（0～10 s），滑体内的破碎体数目增长极为缓慢，主要表现为节理将滑体切割为较完整的层状破碎体，此时模型内破碎体共105个。

图  9  破碎体数量在滑体破坏与堆积过程中的变化曲线

Figure  9.  Variation curve of fragment quantity during landslide failure and deposition process

下载: 全尺寸图片 幻灯片

随着滑体前后缘产生的裂隙逐步扩展贯通，滑体失稳开始运动，滑体运动过程可以分为两个阶段。10～60 s时，岩体剧烈解体，模型内的破碎体数目迅速增加，此时模型内破碎体从105个陡增至约2925个。结合图9所示的滑体运动过程中的破碎特征，可以总结为重力与冒落作用挤压将下层岩体破坏解体形成小粒径破碎体，与此同时运动过程中岩体相互碰撞摩擦加速了破碎的发展，而该阶段（0～60 s）也是滑体运动破碎的过程中主要阶段。在此之后的时间（60～180 s），模型内破碎体数目的缓慢增加至约3800个后基本达到稳定，随后破碎体数目略有上下波动。破碎数目曲线上下浮动的特征说明了在滑体后续的堆积过程中，既存在着5～10 m大粒径破碎体再次解体形成新的破碎体过程，也存在1～ 2 m小粒径破碎体不断破碎分裂成土粒的过程。T = 180 s 时刻，整个模型的破碎体数量为3820个。

3.3.2   破碎体数量变化情况

基于从PFC模拟结果导出的破碎体体积信息，通过对破碎体的等效粒径进行计算，即可得到斜坡不同时刻的所有破碎体的粒径信息。采用目前常用的一种岩块等效粒径定义标准^{[16 − 17]}，对不同时刻的破碎体的粒径进行统计，如式（1）所示：

式中：d——破碎体的等效粒径/m；

V_f——破碎体的体积/m³；

V₀——模型滑体的总体积/m³。

在上述基础上，利用双参数Weibull分布模型来描述不同时刻的破碎体粒径分布^[18]。碎片粒径分布的双参数Weibull模型可用式（2）表示：

式中：P——小于等效粒径d的颗粒累计质量百分比。

d_c——尺寸参数，可用来衡量破碎体中细粒的占比，决定曲线尺寸比例。

β——形状参数，表征破碎体的级配。

基于上述对马达岭滑坡不同阶段破碎情况的总结，本文对T = 10 s（滑体失稳过程）与T = 180 s（滑体运动过程）两个典型时刻的破碎体粒径信息进行分析。图10展示了不同时刻破碎体小于某一粒径百分含量级配分布的Weibull拟合结果。从图中不同时刻的形状参数 也能得出这一结论，即T = 10 s时刻的形状参数β（2.973）较T = 180 s时刻（ β = 2.785）的大。从图9中可得知，在T = 10 s时刻的破碎体粒径范围明显较大，同时存在小粒径破碎体及大量大粒径破碎体。这表明在滑体破坏过程中，滑体内的破碎体粒径分布较广，级配较好，对照了图8及图9中的结果。基于拟合结果得到的尺度参数可以得知，10～180 s 堆积体中的碎片尺度参数产生了小幅的减小，其表明在滑体停止滑动时（T>180 s），经历在坡脚前水平地表的运动堆积，滑体内的细粒径破碎体占比略有增加。双参数 Weibull 模型的拟合结果表明，滑体在水平地表滑动堆积过程中的仍会有一定的破碎现象发生，破碎作用会导致堆积体内的细粒径破碎体增多，破碎体的粒径分布也变得更为集中。

图  10  滑体破坏与堆积过程双参数Weibull破碎体粒径百分含量分布曲线拟合曲线

Figure  10.  Dual-parameter Weibull distribution fitting curve of particle size percentage content distribution in landslide failure and deposition process

下载: 全尺寸图片 幻灯片

4.   结论

（1） 此类下伏采空区含节理的滑坡碎屑流滑体动力致灾过程可归纳为：后缘拉裂及前缘开裂扩展，裂纹沿节理一端向临空面扩展，滑体沿节理被分割成不同的块体，重力及下伏采空区的冒落作用进一步加剧滑体破碎，坡脚剪切及坡肩拉张作用同时发生，节理面前方主滑体形成，滑体沿节理面破碎解体滑下。

（2） 此类滑坡碎屑流在滑体失稳及运动过程中，均存在明显的破碎现象，且整个滑坡碎屑流的主要破碎发生在滑体运动过程中。具体表现为滑体失稳过程中（0～10 s）的切割与在运动堆积过程中（10～180 s）的碰撞、摩擦、拉裂等产生的剧烈破碎。滑坡碎屑流从失稳到停止运动堆积的整个过程中伴随有新破碎体的产生及旧破碎体的粉碎消失。

（3） 采用 Weibull 双参数模型拟合的不同时刻的破碎体粒径分布结果表明滑体在堆积的过程中仍有破碎产生。滑体内的细粒径破碎体不断增多，且破碎体整体以中小粒径为主。