0 引言

自 2018 年以来，受全球生态环境改变、极端天气、城市建设等因素的影响，中国地质灾害越发频繁，严重制约了中国西部地区的工程建设和运营。仅 2020年，中国就发生 7840处地质灾害，其中危岩崩塌 1797 处，占比高达 37.4%。因此，开展对危岩崩塌前的稳定性研究和评价对中国地质灾害的防治具有重要意义。

目前，国内外学者对危岩体稳定性的研究方法主要分为：半定量分析方法、定量分析方法和数值模拟计算。谢全敏和夏元友（2002）在综合考虑了岩体完整性、滑面倾角及风化系数等因素下，利用灰色系统理论对清江危岩体的稳定性进行了评价，并取得了很好成果；陈洪凯等（2004）利用极限平衡理论推导出了危岩体在不同工况、不同破坏模式下的稳定性系数；王根龙等（2012）通过材料力学的弯曲梁应力计算方法，得到了悬臂-拉裂式崩塌的变形破坏模式和稳定性计算方法；郑安兴等（2013）通过对危岩体主控结构面变形、扩展、破坏的有限元模拟分析，得到了岩石抗拉强度、主控结构面倾角与危岩稳定性的关系；董金玉等（2010）通过 UDEC 离散元数值模拟，得到了砂泥岩软弱互层形成的危岩崩塌的破坏模式，据此提出了相应的防治措施； 朱利辉等（2021）从地形地貌、物质组成、崩塌形成机理等方面对独山子崩塌的稳定性进行了评价，并得出了临泽县地质灾害具有不均匀性、突发性、周期性的特点；贾超等（2021）通过对八仙山自然保护区崩塌等地质灾害的形成条件、发育特征及诱发因素的研究，提出了地质灾害的防治措施，为该区域防灾减灾、生态修复等提供了依据。

从上述文献的研究可知，学者们对危岩体的稳定性、破坏模式以及防治措施等进行了系统的研究分析，但对危岩体的稳定性评价和形成、破坏机理的分析多采用单一方式：公式计算或数值模拟计算，使用联合评价的方法较少，不能做到在简单高效获得危岩体稳定性的同时，通过岩土体自身的应力-应变关系、变形情况以及其渐进破坏过程获得岩土体的形成、破坏机理。基于此，本文利用改进后的水压分布假设（舒继森等，2004），以龙寺矿区 CIV-1 边坡危岩为算例，推导了危岩体在多因素影响下的滑移破坏、倾倒破坏的无量纲表达式，并联合使用离散元强度折减法对危岩体的渐进破坏过程进行分析，为后续防护提供设计依据。

1 研究区概况

1.1 自然地理情况

龙寺矿山研究区位于焦作市北部解放区上白作乡龙寺村，地理坐标为：东经 113°11'44″~113° 13'48″，北纬：35° 15'39″~35° 17'27″，面积约 4.61 km（2 图1）。矿区公路既有焦辉主干线，又有通往附近各县的支线通过，交通便利。

1.2 地形地貌

龙寺矿山研究区主要为山地地貌，全区地形西北高，东南低，最高标高 480 m，最低标高 203 m，最大相对高差 277 m。区内沟谷发育，沟谷呈“V”字型，其切割深度约 10~100 m，纵坡降 9%~15%，山顶呈浑圆状。研究区地形地貌如图1所示。

1 —中山；2—低山；3—丘陵；4—冲洪积倾斜平原；5—冲积平原；6—黄土台原；7—洼地；8—盐碱地；9—研究区位置

1.3 地质构造

焦作市位于新华夏太行山隆起的南端与晋东南山字型构造东翼及射弧的前缘和东秦岭纬向构造带之北北缘相交联合弧地带，区内广泛发育以燕山运动以来所形成的多种构造形迹，多以正断层为主，倾角多在60º以上，走向多为北东向。

1.4 水文地质

根据矿区的地层岩性、含水介质及地下水动力条件，地下水类型可分为碳酸盐岩岩溶水和第四系松散岩类孔隙水两类，流向总体上为东西向。

第四系松散岩类孔隙水主要蕴藏在 2~10 m 不等的沟谷或山前的残坡积物层内，主要依靠灌溉回渗、降水入渗和山区洪水补给；碳酸盐岩岩溶水主要赋存于中奥陶统灰岩中，可通过降水直接入渗补给，也可将水汇入河流沟谷后，沿河床线状渗漏补给，若沿裂隙、溶洞、断层破碎带等向山前运移汇集，则形成强径流带；由于富水性好，主要以泉、矿坑排水等方式排泄。

1.5 人类活动

龙寺矿山由于计划经济时期胡乱无序地露天开采，加之当地村民的乱采乱挖，矿区周边引发了诸多环境地质问题。其中，未按台阶开挖的高陡边坡在极端降雨或地震作用下随时可能发生崩塌灾害，急需治理。

1.6 龙寺矿区CIV-1边坡

为分析矿区危岩体的稳定性及破坏模式，以 CIV-1边坡为例展开分析，具体位置如图2所示。

CIV-1 边坡经过人工开挖后，形成了高约 105 m，坡度 75°~85°的圈椅状高陡边坡。坡体岩性主要为浅灰色中厚层状奥陶系灰岩，由于受断层等构造影响，岩体裂隙发育，且在坡向 190°的坡面形成了较多柱状、块状危岩体（图3）。

根据三维激光扫描对结构面数据的采集、分析，CIV-1 边坡内的 IV 结构面较为发育，长数米到数十米不等，裂隙面受岩溶作用强烈，主要钙质填充，局部泥质填充，其类型主要为卸荷裂隙、层面、构造裂隙。经统计分析得到了坡体内的优势结构面和岩层产状（表1）。

J1、J3是由于开挖卸荷产生的卸荷裂隙，近乎与坡面平行，产状201°∠81°、281°∠81°，裂隙面较为平直，多数张开度为0.5 cm，最大处可达到3 cm。

J2 和 J4 为构造节理面，它们倾向相反，倾角在 80°~85°范围内，但走向基本一致，近似为南北向，它们在平面上呈一对共轭的“X”型，裂隙面整体较平直，局部倾角变化，微粗糙，在边坡上延伸可达75～85 m。

层面，较平直，局部有泥质充填，产状210°∠6°。该组结构面可形成边坡块体倾覆的底部边界和块体坠落的控制性结构面。

在节理岩质边坡中，危岩体的破坏模式往往受坡面、层面、构造裂隙以及卸荷裂隙组合关系的控制，不同的组合方式常常形成不同形式的块体破坏模式。基于此，本文利用赤平投影对优势结构面与临空面的组合关系进行了分析（图4），由此可知， CIV-1 边坡岩体被层面、构造裂隙和卸荷裂隙切割成长方体、正方体等形状的块体（图4）。在此基础上，对该边坡危岩体最可能出现的破坏模式进行初判（表2）。

2 基于极限平衡法的危岩稳定性分析

2.1 危岩体的简化地质模型

为推导出危岩体在多因素影响下的滑移破坏、倾倒破坏的无量纲表达式，本节根据上文对危岩体野外地质调查资料的分析，将 CIV-1 边坡存在的 3 个危岩体进行了地质模型概化，简化的地质模型如图5所示。

a—1号危岩体；b—2号危岩体；c—3号危岩体

2.2 基本假设及计算模型

危岩体由于受卸荷裂隙和结构面的共同作用，多形成块状、长柱状的危岩体，但这些危岩体的形状不规则，为便于计算，对模型做出以下假设：①破坏面为一平面，其走向与坡面平行；②后缘陡倾裂隙面和陡倾坡面均直立，裂隙面走向与坡面平行； ③将危岩主控结构面视为贯通面（重庆市国土资源和房屋管理局，2004）；④反倾结构面 EF 在纵向的投影应位于 ED范围。根据大多危岩真实的坡面形态和坡体假设得到了边坡柱状危岩体的失稳概化理论模型，如图6所示。

危岩体的几何要素包括：块体高度 H（m），坡前直立面高度 H_A（m），裂隙深度 Z（m），裂隙积水深度 Z_w （m），裂隙积水深度到地面高度 H_w （m），张裂隙距坡肩的距离 S（m），倾角 β（°），层面 θ（°）。当坡前直立面 AE与边坡坡面的倾向相反时，上图可简化为 2 号危岩体的失稳模型，即块体 BCDEF，此时反倾节理面长度为 L（m），与竖向节理面的夹角为 α（°）；当 H_A（m）为零时，上图又可简化为1号危岩体柱状体的失稳模型。

2.3 力学模型及受力分析

根据危岩体的受力分析知：其稳定性受冻胀力 N、静水压力 V、U、坡脚出溜缝是否阻塞和地震力 K 等因素的影响；不同的受力分布假设、取值对危岩体的稳定性影响较大（刘富有等，2016；师华鹏等， 2016）；而对于动水压力，由于柱状危岩的高宽比大，受到的动水压力较小，忽略其对危岩稳定性的影响（师华鹏等，2015）。根据危岩体的实际情况，利用改进后的水压分布假设，对危岩体在不同水力状态下的受力情况进行分析，如图7所示。

据图5、图6所示的几何关系及相关力学知识分别对各外力及其绕坡址 A（F）发生倾覆破坏时的力臂进行了计算。

a—出溜缝未堵塞且0 ≤ Z_w ≤ 0.5H_w；b—出溜缝未堵塞且0.5H_w ≤ Z_w ≤ H_w；c—出溜缝堵塞且0 ≤ Z_w ≤ H_w

2.3.1 重力计算

对于具备反倾结构面类似2号危岩体的危岩体重力计算时，EFG 部分的重力 W^*（kN/m）、距坡址的水平距 x^* _G （m）、纵向距 y ^* _G （m）应按上式确定；对于坡前直立面垂直的危岩体重力计算时，W^*（kN/m）应取 0，如1号危岩体、3号危岩体；μ为冻胀系数，取值详见文献（刘富有等，2016；师华鹏等，2016）。

2.3.2 冻胀力及地震力计算

式（2）中，f（ kN/m）、K_h（kN/m）、K_v（kN/m）、k_h、k_v 分别代表冻胀力均布荷载、水平地震力、竖向地震力、水平地震及竖向地震系数。对于k_h、k_v的取值及方向规定已做了详细说明（刘富有等，2016；师华鹏等，2016）。

2.3.3 静水压力计算

当取图7（a）的力学模型时，当取图7（b）的力学模型时，当取图7（c）的力学模型时：

2.4 稳定性计算

2.4.1 滑移式危岩体

根据危岩体滑移稳定性的定义，危岩体滑移稳定性系数为坡体抗滑力和下滑力的比值，即危岩体的滑移稳定性系数F_s：

将式（1）～（5）带入式（6），并将分子分母同时除以 γH²，可得到危岩体滑移稳定性系数 F_s 的无量纲表达式，如式（7）所示：

式（7）中，M、Q在不同状态下取值不同，应满足以下关系：

当取图7（a）的力学模型时，

当取图7（b）的力学模型时，

当取图7（c）的力学模型时：

$\left\{\begin{array}{c}M=\frac{{\gamma }_w}{4\gamma \mathrm{s}\mathrm{i}\mathrm{n}\theta }\left[{\left(\frac{H_w}{H}\right)}^2-2{\left(\frac{Z_w}{H}\right)}^2\right]\\ Q=\frac{{\gamma }_w}{2\gamma }{\left(\frac{Z_w}{H}\right)}^2\end{array}\right.$

$\left\{\begin{array}{c}M=\frac{{\gamma }_w}{2\gamma \mathrm{s}\mathrm{i}\mathrm{n}\theta }{\left(\frac{H_w}{H}-\frac{Z_w}{H}\right)}^2\\ Q=\frac{{\gamma }_w}{4\gamma }\left(\frac{4H_w{Z}_w}{H^2}-\frac{H_w^2}{H^2}-2\frac{Z_w^2}{H^2}\right)\end{array}\right.$

$\left\{\begin{array}{c}M=\frac{{\gamma }_w}{2\gamma \mathrm{s}\mathrm{i}\mathrm{n}\theta }\left[{\left(\frac{H_w}{H}\right)}^2-{\left(\frac{Z_w}{H}\right)}^2\right]\\ Q=\frac{{\gamma }_w}{2\gamma }{\left(\frac{Z_w}{H}\right)}^2\end{array}\right.$

2.4.2 倾覆式危岩体

根据倾覆稳定性的定义，块体的倾覆稳定性系数为绕A（F）的抗倾覆力矩与促倾覆力矩之比，即：

将式（1）～（5）带入式（7），并将分子分母同时除以 γH² /6，可得到危岩体倾覆稳定性系数 F_W 的无量纲表达式，如式（9）所示：

式中，M_W、Q_W在不同状态下取值不同，应满足以下关系：

当为图7（a）的力学模型时，

当为图7（b）的力学模型时，

当为图7（c）的力学模型时：

$\left\{\begin{array}{c}{M}_W=\frac{{\gamma }_w}{4\gamma {\mathrm{s}\mathrm{i}\mathrm{n}}^2\theta }\left\{\begin{array}{c}3\left(1-\frac{Z}{H}\right)\left[{\left(1-\frac{Z}{H}\right)}^2-{\left(\frac{Z_w}{H}\right)}^2\right]-\\ {\left(\frac{Z_w}{H}\right)}^3+9{\left(1-\frac{Z}{H}\right)}^2\frac{Z_w}{H}\end{array}\right\}\\ Q_W=\frac{{\gamma }_w}{\gamma }\left[3\left(1-\frac{Z}{H}\right){\left(\frac{Z_w}{H}\right)}^2+{\left(\frac{Z_w}{H}\right)}^3\right]\end{array}\right.$

$\left\{\begin{array}{c}{M}_W=\frac{2{\gamma }_w}{\gamma {\mathrm{s}\mathrm{i}\mathrm{n}}^2{\theta }_2}{\left(1-\frac{Z}{H}\right)}^3\\ Q_W=\frac{{\gamma }_w}{4\gamma }\left[\begin{array}{c}\left.9\frac{Z_w}{H}\left(1-\frac{Z}{H}\right)\left(1-\frac{Z}{H}+\frac{Z_w}{H}\right)+\right]\\ 3{\left(\frac{Z_w}{H}\right)}^3-5{\left(1-\frac{Z}{H}\right)}^3\end{array}\right]\end{array}\right.$

$\left\{\begin{array}{c}{M}_W=\frac{{\gamma }_w}{\gamma {\mathrm{s}\mathrm{i}\mathrm{n}}^2{\theta }_2}{\left(1-\frac{Z}{H}\right)}^2\left(1-\frac{Z}{H}+3\frac{Z_w}{H}\right)\\ Q_W=\frac{{\gamma }_w}{\gamma }\left[3\left(1-\frac{Z}{H}\right){\left(\frac{Z_w}{H}\right)}^2+{\left(\frac{Z_w}{H}\right)}^3\right]\end{array}\right.$

2.4.3 工况及参数选取

根据现场调查，危岩体可能受到暴雨、地震等外荷载影响，为评价其稳定性，本文从自然工况、暴雨工况、地震工况以及暴雨+地震组合工况分别对其稳定性进行了计算。在计算静水压力时，自然状态下不考虑水的影响，暴雨状态时，根据矿区强降雨时现场观测数据，危岩后缘裂隙内水柱高取三分之二后缘裂隙高度；对于地震力，水平及竖向地震系数分别取 0.15、0.1；对于冻胀作用，μ取 0.5进行危岩稳定性的计算；根据现场岩体结构面直剪试验结果，层面、结构面的饱和内聚力和内摩擦角分别取 0.21 MPa、22.8°，0. 09 MPa、17°；各个危岩体的几何尺寸如图5 所示。经计算，各危岩体的稳定性系数如表3所示。

3 基于强度折减法的危岩稳定性分析

根据上文分析计算，得到了危岩体在各工况下的稳定性系数，但由于其理论的局限性，并不能反应岩体的破坏过程和变形情况。基于此，为查清危岩体的形成、破坏机理，本文利用可以处理大量结构面，并允许块体沿结构面发生大位移和转动的 UDEC 软件对岩体的稳定性进行分析。UDEC 利用离散元强度折减法计算边坡的稳定性时，其基本思路是保持边坡的原始状态不变，主要对材料参数如：黏聚力和内摩擦角正切值进行折减，折减后然后进行稳定性计算，不断重复上述过程，直至系统达到平衡状态。当边坡处于平衡状态时，折减系数的倒数即为边坡的安全系数。

3.1 数值模型建立

根据实测的岩层和节理产状、层面和节理的平均间距以及卸荷裂隙的发育情况建立危岩体的离散元计算模型，由于篇幅有限，以表2初判的可能发生滑移破坏的 1号危岩体和可能发生倾倒破坏的 2 号危岩体为例，如图8所示。

a—危岩体1；b—危岩体2

3.2 模型参数选取

本文根据室内试验获得的岩体和结构面参数对模型进行计算，如表4所示。

3.3 模拟成果

通过对参数折减后的模拟计算，得到了危岩体的位移云图、位移矢量图、最大位移值等数据，并据此得到了危岩体的稳定性系数F^*（表3）。由于模拟计算得到的图幅较多，本文仅列出了危岩体在自然工况下的位移矢量图和位移等值线图进行分析 （图9、10）。

a—位移矢量图；b—位移等值线图

4 数据分析

4.1 模拟成果分析

从图9中可以看出 1号危岩体在边坡开挖后的自然工况下具有沿层面和陡倾节理面向临空面方向滑移的趋势，但最终位移和不平衡力收敛，表明了危岩体在自然状态时是处于稳定状态的；而从图9b 上可以看出危岩体由中部向两边位移量逐渐降低，加之结构面并不是完全光滑，就导致了滑移面上存在多个锁固点，继而在锁固点附近产生拉应力高度集中，形成多个近乎垂直滑移面的喇叭状张拉裂隙（图3），这也正解释了边坡开挖后应力重分布导致的卸荷拉裂隙这一现象。通过对1号危岩体在不同荷载效应下的稳定性模拟计算，结果表明危岩体在不同工况下的变化规律与自然工况下的变化规律相同，总体表现为卸荷-拉裂（剪切）-滑移式的破坏模式，利用 PFC 得到的 CIV-1 号边坡 1 号危岩体的破坏模式是一致的（陈鹏宇，2016）。

a—位移矢量图；b—位移等值线图

从图10中可以看出，在自然工况下危岩体 2具有向临空面方向外倾的趋势，但最终位移和不平衡力收敛，表明了危岩体在自然工况时是处于稳定状态的；从图10b上可以看出，危岩体上部的位移要大于其下部位移，这是因为上部裂隙处于完全贯通而后缘裂隙下部存在锁固段；但在地震、暴雨等外因影响下，2号危岩体的位移量逐渐增大，且后缘裂隙面底部的位移量也在增加，表现为后缘裂隙面的扩展与锁固段的张拉破坏，总体表现为卸荷-拉裂-倾倒式的破坏模式，这与利用 PFC 得到的 CIV-1 号边坡2号危岩体的破坏模式是相同的（陈鹏宇，2016）。

4.2 综合成果分析

由表3和（中华人民共和国国土资源部.2006） 可知：在自然状态下，各危岩体处于基本稳定-稳定状态，在外力下，稳定性系数发生了显著的降低，可能发生破坏，为评价外力对危岩体稳定性的影响程度，规定各工况的滑移稳定性为 F_si，倾倒稳定性为F_wi，数值计算的滑坡稳定性为F^* _i，i=1，2，3，4依次代表天然、地震、暴雨以及组合工况。据此，可得到不同工况对危岩稳定性系数的影响程度表以及不同工况情况下危岩稳定性系数变化曲线，如表5、图11 所示。

由图11、表5可知：

（1）各工况对危岩体稳定性的影响程度为：组合工况>暴雨工况>地震工况，且组合工况对其稳定性的影响并不等同于暴雨工况和地震工况的简单叠加。

（2）在各工况影响下，危岩体的滑移稳定性系数下降幅度为：3 号危岩体>1 号危岩体>2 号危岩体；而倾倒稳定性系数下降幅度为：2号危岩体>1号危岩体>3号危岩体；结合1号、3号危岩体模拟计算结果与理论计算结果对比值小于1，而2号危岩体的对比值大于 1 的结果可知，3 号、1 号危岩体在外动力影响下最有可能发生滑移破坏，而 2 号危岩体最有可能发生倾倒破坏，这与表2 利用结构面与临空面相互组合时对危岩体的初判是一致的。

（3）2 号危岩体的 F^* 比 F_s、F_w都大，这是由于模拟计算时后缘裂隙存在锁固段，而理论计算时为简化公式将危岩主控结构面视为贯通面导致的。

（4）由上述分析可知，通过UDEC强度折减法与极限平衡法得到的危岩体破坏模式，稳定性系数及其在不同工况下的变化规律等结果都较为相近，这表明利用 UDEC 在分析危岩体的形成、破坏机理方面是可行的。

a—1号危岩体；b—2号危岩体；c—3号危岩体

5 结论

（1）在暴雨、地震等外力影响下，危岩体的稳定性发生显著减低，极易失稳，发生崩塌等地质灾害。

（2）UDEC在分析危岩体的形成、破坏机理方面是可行的；基于此，得到了CIV-1边坡1号、3号危岩体的形成与破坏模式为卸荷-拉裂（剪切）-滑移式，稳定性受层面与后缘陡倾节理的组合关系控制；2 号危岩体的形成与破坏模式为卸荷-拉裂-倾倒式，并且后缘节理的贯通程度对于危岩体的稳定性具有控制作用。

（3）通过上文极限平衡法和强度折减法的计算结果可知：危岩体受外力影响时，无论采用何种计算方式得到的稳定性系数的变化规律均相同，这也证明了基于上文假设推导的危岩体稳定性无量纲表达式是正确的，具有实用性。

（4）通过极限平衡法和 UDEC 强度折减法相结合的手段对危岩体的稳定性进行评价，不仅可以弥补极限平衡法无法获取危岩渐进破坏过程的缺陷，还可以以公式的形式快速高效地获取危岩体的稳定性，对危岩体的防治具有很好的指导意义。

参考文献