基于尖点突变理论的采场台阶边坡稳定性研究

Study on the Stability of Mining Bench Slopes Based on Cusp Catastrophe Theory

摘要

离散元法适用于模拟露天开采后的非均质岩质台阶边坡在重力荷载作用下的非线性变化，但缺乏明晰量化的边坡失稳判据，难以精确获取边坡安全系数。针对该问题，基于尖点突变理论结合重力增加法，建立了边坡水平方向平均位移与重力系数的尖点突变模型，将边坡临界失稳判据量化为确定的位移突变特征值。通过对不同节理间距和连通率进行组合模拟，揭示了不同组合下边坡临界位移及安全系数变化规律。分析表明：该模型能够有效确定边坡失稳临界位移及安全系数；当节理连通率为0.70时，边坡的临界失稳位移整体高于节理连通率为0.70以上的组合；随着节理间距的增大，边坡临界失稳时沿节理面形成范围更大的滑体，对应较大的临界失稳位移。该研究结果为节理岩体区域边坡位移监测及位移预警阈值提供了理论依据。

Abstract

The discrete element method is particularly effective for simulating the nonlinear deformation of heterogeneous rock slopes subjected to gravitational loads following open-pit mining. Traditional methodologies often suffer from ambiguous criteria for determining slope instability，thereby complicating the accurate assessment of the slope’s safety factor. To address this issue，a cusp catastrophe model was developed，establishing a relationship between the average horizontal displacement of the slope and the gravitational coefficient. This model quantifies the critical instability criterion as a definitive displacement catastrophe value. Using the discrete element program MatDEM，simulations were conducted with varying combinations of joint spacing and connectivity to examine the changes in critical displacement and safety factors under different conditions. Case analyses indicate that this model effectively determines the critical displacement and corresponding safety factor at the point of slope failure. Notably，when joint connectivity is 0.70，the critical instability displacement is generally higher compared to scenarios with greater connectivity. For joint spacings ranging from 1 to 3 meters，an increase in spacing results in a larger sliding mass along the joint surface，which consequently leads to a rise in the critical instability displacement. The gravity increase method，grounded in cusp catastrophe theory，models the evolution of slopes under realistic conditions by incrementally increasing the gravitational load. This approach is particularly well-suited for jointed rock masses that exhibit complex nonlinear deformation behaviors. Unlike the strength reduction method，it obviates the need for recalibrating material parameters before each calculation，thereby enhancing computational efficiency. The study’s findings indicate that joint spacing and connectivity have a significant impact on slope stability. With uniform joint spacing，increased connectivity results in a lower safety factor，whereas with uniform connectivity，larger joint spacing leads to a higher safety factor. These findings are consistent with practical engineering conditions，underscoring the importance of considering joint spacing and connectivity in slope stability analysis.The developed model and results provide a scientific basis for determining monitoring point placement and setting warning thresholds for slope displacement in jointed rock mass regions.

Graphical abstract

关键词

尖点突变理论 / 重力增加法 / 失稳判据 / 临界位移 / 节理间距 / 节理连通率

Key words

cusp catastrophe theory / gravity increment method / instability criteria / critical displacement / joint spacing / joint connectivity rate

近年来，采矿活动后露天采场台阶边坡失稳事故频发，岩体中的节理和裂隙等不连续面显著影响岩体的力学性能，使得边坡整体稳定性受到严重削弱（盛建龙等，2018；王正堂等，2021；孙杰涛等，2024）。节理面对岩质边坡稳定性的负面作用及其力学机制成为岩土工程领域的热点问题。而露天开采后的非均质岩质台阶边坡因其内部结构复杂，呈现出明显的非线性行为，给边坡安全评价与风险预警带来了巨大挑战（杜时贵，2018；荆永滨等，2018）。

离散元法（Discrete Element Method，DEM）能详细模拟岩体破裂、滑动等局部及整体响应，因此被广泛应用于边坡稳定性工程领域。在以往研究中，研究人员通过引入数值方法和理论模型对边坡稳定性进行了若干探索。一方面，通过将离散元法配合强度折减法和重力增加法等数值技术应用于估计边坡安全系数，揭示了在模拟进程中岩体内部破裂及失稳的演化过程（周健等，2009；Scholtès et al.，2012）。另一方面，针对边坡稳定性的判据问题，提出了基于位移突变的失稳判据，尝试通过监测位移响应来识别临界失稳状态（徐卫亚等，2007；Li et al.，2023）。与此同时，非线性动力学理论尤其是尖点突变模型被引入解释岩体在外界荷载作用下参数的突变特征，为判断斜坡失稳提供了新的理论依据（秦四清，2000；Qin et al.，2001）。尽管上述研究为揭示复杂岩体失稳机制做出了积极探索，但是所采用的方法在量化临界失稳阈值、明确安全系数指标方面仍存在不足，难以直接应用在实际工程中。

针对上述不足，通过对顺层岩质边坡在不同重力系数条件下向临界平衡状态演化的过程进行离散元分析，提出了一种基于尖点突变理论的顺层岩质边坡稳定性评价方法，即重力增加法。将边坡临界失稳判据量化为确定的位移突变特征值，并对不同节理间距和连通率进行组合模拟，揭示了不同组合下边坡临界位移及安全系数变化规律。该研究不仅为复杂地质条件下边坡失稳问题提供了一种明确、量化的判据，也为工程实践中岩质台阶边坡安全性评价提供了新思路。

1 基于突变理论的边坡稳定性重力增加法

1.1 基于尖点突变理论的重力增加法

重力增加法的基本原理与强度折减法正好相反，其基本思想是在控制岩体抗剪强度不变的前提下，通过单调增加重力加速度或土体重度，改变边坡的体荷载（重力），实现边坡的超载直至边坡失稳。以边坡临界失稳破坏时的重力系数与实际值的比值作为边坡的安全系数（Swan et al.，1999），表示为

F=\frac{{\mathrm{g}}^{\mathrm{\text{'}}}}{\mathrm{g}}=\frac{{\gamma }^{\mathrm{\text{'}}}}{\gamma }（1）

式中：F为安全系数；{\mathrm{g}}^{\mathrm{\text{'}}}和{\gamma }^{\mathrm{\text{'}}}分别为临界失稳破坏时的重力加速度和土体重度；\mathrm{g}和\gamma分别为实际的重力加速度和土体重度。

突变理论由数学家Thom于1972年创立（凌复华，1984），利用数学模型来描述系统在不同控制参数下的状态变化，当某些控制变量（如荷载和应力等）达到某一临界值时系统突然失稳。尖点突变模型是滑动失稳分析中常用的一类突变分析模型，系统的状态函数可用包含一个变量的多次项来表示。顺层岩质边坡的失稳过程既是结构面贯通破坏的过程，又是边坡相应破坏部位位移突变的过程。运用重力增加法进行边坡稳定性分析时，通过离散元数值模拟，导出了多级重力系数F_i条件下边坡水平方向平均位移S的变化数据。基于此，建立以S为系统状态变量、以F_i为控制变量的突变模型状态函数S(F_i)。具体地，对每级重力系数下的S值进行拟合，构建S与F_i曲线的泰勒级数形式。为满足精度要求，截取至4次项构建该函数：

S\left(F_i\right)=b_{\mathrm{0}}+b_{\mathrm{1}}{F}_i+b_{\mathrm{2}}{F_i}^{\mathrm{2}}+b_{\mathrm{3}}{F_i}^{\mathrm{3}}+b_{\mathrm{4}}{F_i}^{\mathrm{4}}（2）

式中：b_{\mathrm{0}}、b_{\mathrm{1}}、b_{\mathrm{2}}、b_{\mathrm{3}}和b_{\mathrm{4}}为待定系数。

运用多项式变换平移根，令F_i=p-q,q=\frac{b_{\mathrm{3}}}{\mathrm{4}{b}_{\mathrm{4}}}则可消去三次方项将上述多项式方便转换为尖点突变理论的基本函数形式K\left(p\right)，表示为

K\left(p\right)=d_{\mathrm{4}}{x}^{\mathrm{4}}+d_{\mathrm{2}}{x}^{\mathrm{2}}+d_{\mathrm{1}}x+d_{\mathrm{0}}（3）

其中，d_i和b_i的关系可表示为

\left[\begin{array}{c}{d}_{\mathrm{0}}\\ d_{\mathrm{1}}\\ d_{\mathrm{2}}\\ d_{\mathrm{4}}\end{array}\right]=\left[\begin{array}{ccccc}{q}^{\mathrm{4}}& -q^{\mathrm{3}}& q^{\mathrm{2}}& -q& \mathrm{1}\\ -\mathrm{4}{q}^{\mathrm{3}}& \mathrm{3}{q}^{\mathrm{2}}& -\mathrm{2}q& \mathrm{1}& \mathrm{0}\\ \mathrm{6}{q}^{\mathrm{2}}& -\mathrm{3}q& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}{b}_{\mathrm{4}}\\ b_{\mathrm{3}}\\ b_{\mathrm{2}}\\ b_{\mathrm{1}}\\ b_{\mathrm{0}}\end{array}\right]（4）

式（3）两边同时除以d_{\mathrm{4}}，得到：

V\left(p\right)=p^{\mathrm{4}}+u{p}^{\mathrm{2}}+vp+\mathrm{C}（5）

式中：u=\frac{b_{\mathrm{2}}}{b_{\mathrm{4}}}–\frac{\mathrm{3}{b}_{\mathrm{3}}^{\mathrm{2}}}{\mathrm{8}{b}_{\mathrm{4}}^{\mathrm{2}}},\nu =\frac{b_{\mathrm{1}}}{b_{\mathrm{4}}}-\frac{b_{\mathrm{2}}{b}_{\mathrm{3}}}{\mathrm{2}{b}_{\mathrm{4}}^{\mathrm{2}}}+\frac{b_{\mathrm{3}}^{\mathrm{3}}}{\mathrm{8}{b}_{\mathrm{4}}^{\mathrm{3}}}，为控制变量；p为状态变量；\mathrm{C}为常数项，对突变无意义；V(p)即为尖点突变模型的标准势函数。

对势函数V(p)求导，联立V\mathrm{\text{'}}(p)=\mathrm{0} 、 V″\left(p\right)=\mathrm{0}求解得到尖点突变的分叉集方程：

\Delta =\mathrm{8}{u}^{\mathrm{3}}+\mathrm{27}{v}^{\mathrm{2}}（6）

由突变理论可知（陈忠辉等，1997），当\Delta >\mathrm{0}时，即u、\nu处于分歧点集外，边坡系统位移稳定；当\Delta=0时，即u、\nu处于分歧点边界，边坡系统位移处于临界平衡状态，此时控制变量的微小变化均会引起平衡状态的改变；当\Delta<0时，即u、\nu处于分歧点集内，边坡系统位移突变，边坡会发生失稳滑动。

1.2 基于尖点突变理论的重力增加法的实现过程

岩体的抗剪强度是顺层岩质边坡破坏的主要影响因素。其抗剪强度表达式为

T=c+\sigma \mathrm{t}\mathrm{a}\mathrm{n}\phi（7）

式中：T为抗剪强度；c为黏聚力；\sigma为正应力；\phi为内摩擦角。

使用重力增加法计算边坡稳定性时，随着重力的增加，岩体正应力\sigma会相应增大，导致边坡抗剪强度增加，人为提高了边坡稳定性，相较于强度折减法和极限平衡法稳定性计算结果偏大。为此，孟永东等（2015）通过降低颗粒摩擦系数，弱化模型因重力增加而导致抗剪强度提高的负面影响，得出与强度折减法相似的结果。李雨倩等（2014）在研究中，尝试在增加重力加速度的同时，对内摩擦角进行同等比例的折减计算。结果显示，这种计算方法所得到的结果与传统有限元法的结果十分接近，验证了该方法的可行性。在离散元模型中，黏聚力主要受颗粒间的连接强度影响，而受颗粒单元摩擦系数的影响相对较小。此外，在离散元程序MatDEM中，内摩擦系数与颗粒单元摩擦系数之间存在如下宏微观参数转换公式（周健等，2008；刘春等，2020）：

\mathrm{t}\mathrm{a}\mathrm{n}\phi ={\mu }_{\mathrm{i}}=\frac{P-\mathrm{1}}{\mathrm{2}\sqrt[]{P}},P=\frac{\mathrm{2}+\sqrt[]{\mathrm{2}}{\mu }_{\mathrm{p}}}{\mathrm{1}-\sqrt[]{\mathrm{2}}{\mu }_{\mathrm{p}}}（8）

式中：{\mu }_{\mathrm{i}}为内摩擦系数；{\mu }_{\mathrm{p}}为颗粒单元摩擦系数。

据此，对本次使用基于重力增加法的MatDEM离散元模型进行调整：通过内摩擦系数{\mu }_{\mathrm{i}}与颗粒单元摩擦系数{\mu }_{\mathrm{p}}的宏微观力学参数转换公式计算，使内摩擦系数{\mu }_{\mathrm{i}}的降低倍数取值与重力增加倍数相同；保持边坡的抗剪强度T不变进行计算，得到合理的安全系数。

基于尖点突变理论的顺层岩质边坡稳定性重力增加法的实现过程如下：

（1）首先堆积颗粒、切割生成边坡模型，提取节理范围内的颗粒单元号；

（2）将训练好的初始强度参数分别赋予各单元组，在自然重力条件下进行标准平衡；

（3）采用标准平衡后的模型，逐渐增大重力系数F_i，同时相应降低边坡岩体的内摩擦系数{\mu }_i；

（4）记录不同重力系数F_i下，顺层岩质边坡水平方向的平均位移，对位移S和重力系数F_i的相关数据点进行4次多项式拟合，得到拟合方程S(F_i)；

（5）根据重力增加法下顺层岩质边坡位移尖点突变模型的突变判据，确定顺层岩质边坡的安全系数。

2 顺层岩质边坡安全系数分析

本研究以某磷矿露天采场台阶边坡作为模拟对象，该采场边坡地质结构类型为层状岩体同倾边坡，节理产状为50°~65°，节理间距和贯通程度分布不均。为便于研究，模拟时将其中节理分布简化为倾角与坡角一致且平行均匀分布的节理，采用不同间距和连通率组合的边坡模型来研究顺层岩质边坡的稳定性。

本次顺层岩质边坡的计算模型中，边坡坡角设为60°，高度设为24 m，节理间距为d，节理连通率为k，表示为

k={\displaystyle \sum _{i=\mathrm{1}}^n} l_i\cdot L^{-\mathrm{1}}（9）

其中，各参数及最终形成的顺层岩质边坡模型如图1所示。

模型计算参数取自该磷矿山露天采场台阶边坡，岩石分布主要为中风化白云岩，受采矿开挖振动等因素影响，岩体较破碎，结构面分离，起伏粗糙，岩块、碎屑填充，结合程度一般。除上述内摩擦系数与颗粒单元摩擦系数之外，在MatDEM软件中还存在如下宏观和微观力学参数转换公式，通过MatDEM自动训练材料模块不断调整参数比率进行单轴压缩和剪切等模拟试验，可获得力学性质接近实际值的训练值（刘春等，2020）。

K_{\mathrm{n}}=\frac{\sqrt[]{\mathrm{2}}Ed}{\mathrm{4}\left(\mathrm{1}-\mathrm{2}v\right)}（10）

K_{\mathrm{s}}=\frac{\sqrt[]{\mathrm{2}}\left(\mathrm{1}-\mathrm{5}v\right)Ed}{\mathrm{4}\left(\mathrm{1}+v\right)\left(\mathrm{1}-\mathrm{2}v\right)}（11）

X_{\mathrm{b}}=\frac{\left(\mathrm{3}{K}_{\mathrm{n}}+K_{\mathrm{s}}\right)T_{\mathrm{u}}{d}^{\mathrm{2}}}{\mathrm{6}\sqrt[]{\mathrm{2}}{K}_{\mathrm{n}}\left(K_{\mathrm{n}}+K_{\mathrm{s}}\right)}（12）

F_{\mathrm{s}\mathrm{0}}=\frac{\left(\mathrm{1}-\sqrt[]{\mathrm{2}}{\mu }_{\mathrm{p}}\right)C_{\mathrm{u}}{d}^{\mathrm{2}}}{\mathrm{6}}（13）

式中：宏观参数包括弹性模量（E）、泊松比（\nu）、抗拉强度（T_{\mathrm{u}}）、抗压强度（C_{\mathrm{u}}）和内摩擦系数（{\mu }_{\mathrm{i}}）；细观参数包括正向刚度（K_{\mathrm{n}}）、切向刚度（K_{\mathrm{s}}）、断裂位移（X_{\mathrm{b}}）、初始抗剪力（F_{\mathrm{s}\mathrm{0}}）和颗粒单元直径（d）。岩体具体的物理力学参数见表1，细观参数转换结果见表2，模型颗粒半径取值为0.05 m，各参数训练值与实际值平均相对误差为4.46%。生成节理时使用光滑节理模型，将筛选出的节理范围内颗粒黏结移除（Asadi et al.，2012）。颗粒单元连接情况如图2所示。

2.1 不同节理间距和连通率的顺层岩质边坡安全系数分析

顺层岩质边坡节理间距分别取1.0，1.5，2.0，2.5，3.0 m，节理连通率分别取0.90，0.85，0.80，0.75，0.70，研究不同节理间距和连通率的节理分布对顺层岩质边坡稳定性的影响。以节理间距为2.0 m，连通率为0.70时的分布情况为例，说明求解顺层岩质边坡安全系数的过程。

通过不断增加重力系数，得到节理间距为2.0 m，连通率为0.70时，边坡水平方向平均位移S随重力增加系数F_i的演化规律及泰勒级数拟合曲线（图3）。由图3可知：随着重力系数F_i的增大，边坡水平方向平均位移逐渐增大，当F_i<1.64时，位移增长缓慢，变化平缓；当F_i\ge \mathrm{1.65}时，边坡裂隙迅速发育，节理面出现贯通破坏。将水平方向平均位移S和重力系数F_i的泰勒级数拟合曲线中各参数代入式（5）和突变特征方程［式（6）］，得到不同重力系数F_i下顺层岩质边坡水平方向平均位移突变特征值\Delta。\Delta随重力系数F_i的演化曲线如图4所示，当F_i=\mathrm{1.64}时， u=-\mathrm{0.03475}， v=\mathrm{0.01313}， \Delta =\mathrm{0.00432}>\mathrm{0}， 而当F_i=\mathrm{1.65}时，u=-\mathrm{0.07640}，v=\mathrm{0.00428}，，边坡失稳。因此，认为F_i=\mathrm{1.64}为顺层岩质边坡的临界重力系数，即节理间距为2.0\mathrm{m}，连通率为0.70时，该顺层岩质边坡的安全系数为1.64，拟合曲线S\left(F_i\right)=\mathrm{0.0231}{F_i}^{\mathrm{4}}-\mathrm{0.1143}{b}_{\mathrm{3}}{F_i}^{\mathrm{3}}+\mathrm{0.2112}{F_i}^{\mathrm{2}}-\mathrm{0.1725}{F}_i+\mathrm{0.0525}。

图5所示为F_i=\mathrm{1.65}时边坡失稳的水平方向位移变化过程，随着计算步的增加，边坡的位移分布不均匀，沿节理倾向形成较明显的滑裂面。由图5（a）和图5（b）可知，坡脚处最先出现集中位移，并具有向上发展趋势；图5（c）和图5（d）中共面节理逐渐贯通形成滑裂面，非共面节理间出现裂纹切割岩体，形成位移滑动体，此时边坡形态明显改变，位移区域明显。

由图6可以看出，在坡体的中下部区域，尤其是边坡坡脚位置，单元合力明显大于其他岩体区域。随着计算步数的增加，近坡面节理面及坡脚位置的单元合力逐渐集中，且其方向主要沿着裂纹向坡面外扩展。这一现象与位移分布规律一致，表明坡体变形集中在这些关键区域。

图7所示为边坡新增裂隙数目随重力系数的变化曲线。由图7可知，随着重力系数的增大，顺层岩质边坡的裂隙数目逐渐增多，削弱了岩体的整体强度。表明随着重力作用的增强，坡体内部的应力重新分布，岩体产生更多的裂隙，边坡结构的完整性遭到破坏。当F_i>1.64时，边坡裂隙迅速发育，新增裂隙数目发生突变，如图8所示。当F_i=\mathrm{1.65}时，新增裂隙发育明显，且主要为剪裂隙沿近坡面节理面发育，剪裂隙的形成常与剪切应力有关，是边坡失稳的一个重要标志。此时近坡面共面节理贯通，形成滑裂面，非共面节理间出现裂纹，切割岩体形成潜在滑动体，剪裂隙集中发育的近坡面区域发生剪切破坏，边坡处于失稳状态。

将不同间距和连通率节理面的顺层岩质边坡的安全系数列于表3。由表3可知，在节理间距保持一致的情况下，节理连通率越大，边坡的安全系数就越低；而在节理连通率相同的情形下，节理间距越大，边坡的安全系数则越高。这一规律与实际工程中的情况相契合，表明节理间距和节理连通率是影响边坡安全性的关键因素。

2.2 顺层岩质边坡稳定性和节理间距、连通率的关系

在节理连通率为0.70、0.80和0.90组合下，顺层岩质边坡的临界位移变化曲线如图9所示。当节理连通率为0.70时，边坡失稳临界位移整体较大；随着节理间距的增大，边坡失稳临界位移也逐渐增大。研究结果表明，同一边坡模型在不同节理几何产状的影响下，边坡岩体变形的位置和范围不同，临界失稳状态所对应的临界位移也不同。图10所示为节理间距为1.0 m、节理连通率为0.80和节理间距为2.0 m、节理连通率为0.80时边坡临界失稳状态的位移图，可以看出随着节理间距的增大，边坡失稳时非共面节理间较大的间距将产生较长的裂纹切割岩体，而较高的节理连通率使得岩体破坏形成更大的滑动体。这一结果与王旭一等（2021）的层状岩体试件压缩力学试验研究相吻合，即层面间距的增大，增强了试件抵抗变形的能力，在本文模拟模型中表现为随着节理间距的增大，发生位移的颗粒数量和范围逐渐增加，因此最终边坡失稳所对应的临界位移较大，其试验的破坏变形特性与本文数值模拟基本一致。由二次多项式拟合得到不同节理连通率下顺层岩质边坡安全系数与节理间距的关系式见表4，令F_i为具体值可得到相应安全系数下边坡节理间距。图11所示为节理连通率为0.70、0.75、0.80、0.85和0.90的顺层岩质边坡的安全系数和节理间距的二次多项式拟合关系曲线。

由图11中曲线整体变化趋势可知，随着节理连通率的降低，边坡安全系数增大，说明边坡岩体抗剪强度有所增加，这一结果与胡波等（2011）的节理岩体模型试验研究结果相吻合，即低连通率时抗剪强度下降速度慢，高连通率时抗剪强度下降速度快。值得注意的是，当节理连通率达到0.85时，节理间距的增大对安全系数的影响具有显著的非线性特征，边坡稳定性主要受较高连通率的节理面影响，安全系数表现为上升后趋于平缓；而节理连通率为0.70时，随着节理间距的增大，边坡岩体宏观力学性质逐渐趋向于岩石，安全系数在此次模拟的节理间距变化范围内依然呈上升趋势，数值模拟与试验结果一致，进一步验证了研究结论。

3 结论

（1）基于尖点突变理论的重力增加法，直接模拟了导致边坡失稳的增量荷载条件，得出明确量化的失稳判据，为离散元程序确定边坡临界位移，以及边坡变形稳定性分析提供了高效的安全系数计算方法。

（2）针对不同节理间距和连通率组合对岩质边坡变形稳定性的影响，结合离散元数值模拟得到的边坡临界失稳状态，可应用于实际边坡工程中的安全监测，识别出位移和单元合力集中的区域，如边坡表面、坡脚和节理密集带等，重点关注节理面连通形成的潜在滑移面，尤其是在节理连通率高、节理间距小的区域增加监测点的密度，以捕捉细微的变形和位移变化信息。

（3）设定基于临界失稳位移的预警阈值，可参考数值模拟中边坡失稳前的临界位移值，结合现场监测数据的变化趋势，制定科学合理的预警阈值，分级设定预警水平，当监测位移接近或超过临界值时及时发出预警。

参考文献

[1] Asadi M S， Rasouli V， Barla G，2012.A bonded particle model simulation of shear strength and asperity degradation for rough rock fractures［J］.Rock Mechanics and Rock Engineering，45（5）：649-675.

[2] Chen Zhonghui， Tang Chun’an， Fu Yufang，1997.Study on sudden jump of deformation in unstable failure of rocks［J］.Journal of Engineering Geology，5（2）：143-149.

[3] Du Shigui，2018.Method of equal accuracy assessment for the stability analysis of large open-pit mine slopes［J］.Chinese Journal of Rock Mechanics and Engineering，37（6）：1301-1331.

[4] Hu Bo， Zhang Nan， Wang Sijing，et al，2011.Model test and strength analysis research on intermittent joint rock mass［J］.Chinese Journal of Underground Space and Engineering，7（4）：657-665.

[5] Jing Yongbin， Zhao Xintao， Feng Xinglong，2018.Study on 3D simulation of rock fraghmentation for jointed rock mass［J］.Gold Science and Technology，26（3）：357-364.

[6] Li M， Yue Z F， Ji H G，et al，2023.Numerical analysis of interbedded anti-dip rock slopes based on discrete element modeling：A case study［J］.Applied Sciences，13（23）：12583.

[7] Li Yuqian， Li Guangyuan，2014.Application of FEM gravity increase method in stability analysis of inhomogeneous slopes［J］.Science Technology and Engineering，14（33）：290-293，299.

[8] Ling Fuhua，1984.Catastrophe theory：History，current situation and future［J］.Advances in Mechanics，14（4）：389-404.

[9] Liu Chun， Le Tiancheng， Shi Bin，et al，2020.Discussion on three major problems of engineering application of the particle discrete element method［J］.Chinese Journal of Rock Mechanics and Engineering，39（6）：1142-1152.

[10] Meng Yongdong， Su Qingming， Lu Weiping，et al，2015.Research on improved grain flow gravity increase method［J］.Water Resources and Power，33（2）：149-151.

[11] Qin S， Jiao J J， Wang S，2001.A cusp catastrophe model of instability of slip-buckling slope［J］.Rock Mechanics and Rock Engineering，34（2）：119-134.

[12] Qin Siqing，2000.Nonlinear catastrophy model of slope instability and chaotic dynamics mechanism of slope ecolution process［J］.Chinese Journal of Rock Mechanics and Engineering，19（4）：486-492.

[13] Scholtès L， Donzé F V，2012.Modelling progressive failure in fractured rock masses using a 3D discrete element method［J］.International Journal of Rock Mechanics and Mining Sciences，52：18-30.

[14] Sheng Jianlong， Zhai Mingyang，2018.Reliability analysis and sampling method comparison of Jinjiling rock slop based on stochastic response surface method［J］.Gold Science and Technology，26（3）：297-304.

[15] Sun Jietao， Li Haifeng，2024.Study on stability of rock slope based on safety and stability rate［J］.Chinese Journal of Rock Mechanics and Engineering，43（Supp.2）：3872-3884.

[16] Swan C C， Seo Y K，1999.Limit state analysis of earthen slopes using dual continuum/FEM approaches［J］.International Journal for Numerical and Analytical Methods in Geomechanics，23（12）：1359-1371.

[17] Wang Xuyi， Huang Shuling， Ding Xiuli，et al，2021.Study on the effect of inhomogeneous bedding plane on the mechanical properties of uniaxial compression of layered rock mass［J］.Rock and Soil Mechanics，42（2）：581-592.

[18] Wang Zhengtang， Zhang Qi， Wang Chenlong，et al，2021.Influence of joint geometrical parameters on mechanical properties of rock mass［J］.Chinese Journal of High Pressure Physics，35（6）：145-157.

[19] Xu Weiya， Xiao Wu，2007.Study on slope failure criterion based on strength reduction and gravity increase method［J］.Rock and Soil Mechanics，28（3）：505-511.

[20] Zhou Jian， Jia Mincai，2008.Experimental Investigation and Numerical Simulation of Mesoscopic Models in Geotechnical Engineering［M］.Beijing：Science Press.

[21] Zhou Jian， Wang Jiaquan， Zeng Yuan，et al，2009.Slope safety factor by methods of particle flow code strength reduction and gravity increase［J］.Rock and Soil Mechanics，30（6）：1549-1554.

[22] 陈忠辉，唐春安，傅宇方，1997.岩石失稳破裂的变形突跳研究［J］.工程地质学报，5（2）：143-149.

[23] 杜时贵，2018.大型露天矿山边坡稳定性等精度评价方法［J］.岩石力学与工程学报，37（6）：1301-1331.

[24] 胡波，张楠，王思敬，等，2011.断续节理岩体模型试验及强度分析研究［J］.地下空间与工程学报，7（4）：657-665.

[25] 荆永滨，赵新涛，冯兴隆，2018.节理岩体矿岩块度三维模拟研究［J］.黄金科学技术，26（3）：357-364.

[26] 李雨倩，李光元，2014.改进的重力加载比例法在非均质边坡稳定性分析中的应用［J］.科学技术与工程，14（33）：290-293，299.

[27] 凌复华，1984.突变理论：历史、现状和展望［J］.力学进展，14（4）：389-404.

[28] 刘春，乐天呈，施斌，等，2020.颗粒离散元法工程应用的三大问题探讨［J］.岩石力学与工程学报，39（6）：1142-1152.

[29] 孟永东，苏情明，卢伟平，等，2015.颗粒流重力增加法的改进研究［J］.水电能源科学，33（2）：149-151.

[30] 秦四清，2000.斜坡失稳的突变模型与混沌机制［J］.岩石力学与工程学报，19（4）：486-492.

[31] 盛建龙，翟明洋，2018. 基于随机响应面法的金鸡岭岩质边坡可靠度分析及抽样方法对比［J］.黄金科学技术，26（3）：297-304.

[32] 孙杰涛，李海枫，2024.基于安全稳定率的岩质边坡稳定问题研究［J］.岩石力学与工程学报，43（增2）：3872-3884.

[33] 王旭一，黄书岭，丁秀丽，等，2021.层状岩体单轴压缩力学特性的非均质层面影响效应研究［J］.岩土力学，42（2）：581-592.

[34] 王正堂，张祺，王晨龙，等，2021.节理几何参数对岩体力学特征的影响［J］.高压物理学报，35（6）：145-157.

[35] 徐卫亚，肖武，2007.基于强度折减和重度增加的边坡破坏判据研究［J］.岩土力学，28（3）：505-511.

[36] 周健，贾敏才，2008.土工细观模型试验与数值模拟［M］.北京：科学出版社.

[37] 周健，王家全，曾远，等，2009.颗粒流强度折减法和重力增加法的边坡安全系数研究［J］.岩土力学，30（6）：1549-1554.

基金资助

国家自然科学基金项目“高海拔高寒矿区边坡失稳与地表形变行为模式研究”(53264020)