考虑土体参数空间变异性的3维土质边坡可靠度分析

Reliability Analysis of Three-dimensional Soil Slopes Considering Spatial Variability of Soil Parameters

摘要

3维土质边坡可靠度分析通常面临临界滑动面搜索效率低下的挑战。为此，本文采用逐步协方差矩阵分解方法生成边坡土体参数3维对数正态随机场，改进粒子群算法的迭代终止条件并与3维Bishop法结合搜索最小安全系数，最后通过蒙特卡洛模拟得到3维土质边坡的失效概率。通过算例的对比分析验证了本文方法和程序的正确性，针对不同滑动面形式（“圆柱”和“圆柱+椭球”组合）、多种相关长度和变异系数对失效概率的影响进行系统研究。结果表明：改进粒子群算法在保证计算精度的同时可以显著提高其计算效率；土体参数空间变异性是3维土质边坡稳定性高于2维边坡的主要因素之一，空间变异性和3维效应之间的耦合作用显著影响失效概率；土体参数变异系数、相关系数及相关长度的改变，均能导致失效概率的显著波动。研究成果为大规模随机场条件下的3维土质边坡可靠度分析提供了一种高效的方法。

Abstract

Objective Three-dimensional (3D) slope reliability analysis encounters two primary challenges: 1) the stochastic modeling of soil properties involves complex spatial variability, requiring extensive data processing and imposing substantial computational demands; and 2) brute-force trial-and-error methods for locating critical failure surfaces are computationally inefficient, hindering real-world implementation despite the need for precise stability evaluations in civil infrastructure. Therefore, this study employs the covariance matrix decomposition method to generate 3D lognormal random fields for soil parameters, enabling efficient modeling of spatial variability. The particle swarm optimization (PSO) algorithm is refined with enhanced termination criteria and integrated with the 3D Bishop method to search for the minimum factor of safety (F). This approach significantly improves the accuracy and computational efficiency of slope reliability analysis. Methods The core of this framework lies in employing a stepwise covariance matrix decomposition method to generate 3D lognormal random fields for slope soil parameters. The stepwise covariance matrix decomposition efficiently partitions the problem into tractable components, drastically reducing computational demands while enabling high-resolution random field generation with limited resources by systematically decomposing the covariance matrix. This decomposition facilitated rapid simulation of spatially variable soil parameters without excessive overhead. An enhanced PSO algorithm was proposed for 3D slope stability analysis to complement the random field modeling. The algorithm integrated the 3D Bishop method, a well-established limit equilibrium technique, into PSO's global search mechanism, combining PSO's exploratory capacity with the Bishop method's precise stability calculations. This interaction enhanced the identification of critical slip surfaces. Refined iteration termination criteria were incorporated to expedite convergence toward the minimum factor of safety (F) and the corresponding critical surface, improving computational efficiency for reliability analysis. The failure probability of 3D slopes was assessed through Monte Carlo simulation, which accounted for inherent soil property uncertainties. The method's accuracy and effectiveness were validated through numerical examples that systematically analyzed diverse slip surface geometries (for example, cylindrical and cylindrical-ellipsoidal combinations). Results and Discussions The proposed method was applied to analyze Example 1. Excluding spatial variability of soil parameters and assuming a cylindrical sliding surface, a deterministic analysis yielded an F of 1.352 4, which closely matched previous result of 1.352 5. For a sliding surface combining cylindrical and ellipsoidal shapes, F decreased as the failure surface width (B) increased, gradually approaching the two-dimensional analysis outcome. When spatial variability of soil parameters was incorporated, the calculated failure probability (P_f) range of 0.090 5‒0.091 6 (average 9.11%) for Example 1 closely aligned with previous finding of 9.25%, with a relative error of approximately 1.5%. The comparative analysis of safety factors and failure probabilities confirmed the accuracy of the proposed method and the computational procedure. Systematic studies were conducted for various slip surface forms that encompassed combinations of cylindrical and cylindrical + ellipsoidal shapes, as well as different correlation lengths and coefficients of variation. When the sliding surface was cylindrical, and the characteristic length scale in the out-of-plane direction (l_y ) was set to 20 m, P_f gradually decreased as B increased. This observation indicated that the spatial variability of soil parameters was one of the pivotal factors contributing to the enhanced stability assessment accuracy of 3D slope analysis compared to its two-dimensional counterparts. When l_y approaches infinity (l_y →+∞), as B increases, P_f gradually rises while the average safety factor (F_av ) decreases. These findings indicated that B significantly influences the stability of 3D slopes. When the coefficient of variation for cohesion (V_c ) increased from 0.1 to 0.6, P_f rose from 0.0 to 0.25, reflecting an approximate increment of 0.25. When the coefficient of variation for the internal friction angle (V_φ ) increased from 0.1 to 0.4, P_f increased from 0.04 to 0.14, corresponding to an approximate increment of 0.10. When the correlation coefficient between cohesion and the internal friction angle (ρ_c,φ ) increased from -0.5 to 0.5, P_f rose from 0.055 to 0.070, indicating an approximate increment of 0.015. These findings demonstrated that V_c, V_φ, and ρ_c,φ all exert significant influences on the failure probability. When the characteristic length scale in the x-direction (l_x ) increased from 10 m to 40 m, P_f increased by approximately 0.03. In contrast, when l_y increased from 10 to 40 m, P_f increased by approximately 0.075. This result indicated that l_y had a more pronounced impact on slope failure probability compared to l_x . In addition, when the characteristic length scale in the z-direction (l_z ) increased from 1 to 4 m, P_f increased by approximately 0.05, indicating an influence on failure probability. Conclusions The results demonstrate that the enhanced PSO algorithm, when integrated with the 3D Bishop method, markedly improves the computational efficiency of 3D slope reliability analysis without compromising accuracy. When only spatial variability is considered, the failure probability decreases as the width of the failure surface increases. In contrast, when only the 3D effect is considered, the failure probability increases with the width of the failure surface. When both factors are considered simultaneously, their effects on the failure probability tend to offset each other, resulting in a significantly lower failure probability for 3D slopes compared to 2D slopes. When accounting for the variability of soil parameters and 3D effects, the failure probability of slopes initially increases and then decreases with the widening of the failure surface. This pattern indicates that 3D effects dominate when the failure surface is narrow, because geometric constraints and stress redistribution play a critical role. In contrast, spatial variability of soil properties becomes the primary factor when the failure surface exceeds a critical width, because material heterogeneity governs the failure process. The analysis reveals a critical failure surface width that marks the transition from geometry-controlled to material-controlled failure mechanisms in slope stability assessment.

Graphical abstract

关键词

边坡 / 失效概率 / 协方差矩阵分解方法 / 3维效应 / 空间变异性

Key words

soil slope / failure probability / CMD / three-dimensional effect / spatial variability

本刊网刊

边坡工程通常采用安全系数评估其稳定性，但土体因地质历史和矿物成分差异呈现空间变异性^[1‒3]，相同安全系数边坡的失效概率可能显著不同^[4]，单一安全系数难以全面反映边坡稳定状态。为此，引入可靠度理论，通过失效概率指标综合评估土体空间变异性对边坡稳定性的影响^[5‒6]。

工程实践表明，边坡破坏面具有显著的3维空间效应^[7]，而2维分析方法忽略了破坏面的实际宽度和边界条件，通常得出偏于保守的分析结果^[8‒9]。为此，学者结合随机场理论提出了多种3维可靠度分析方法。Griffiths等^[10]将随机场理论和3维有限单元法相结合，使用3维随机有限单元法（RFEM）揭示了2维分析方法低估边坡失效概率的现象。Hicks等^[11]采用RFEM分析了长边坡的稳定性，发现了2维分析无法体现的3种3维失稳模式^[12]。Kasama等^[13]将RFEM和随机响应面法相结合，分析了地震诱发的3维滑坡的可靠度。与此同时，Varkey等^[14]提出了一种改进的半解析方法，用于精确计算3维边坡的可靠度；Shu等^[15]基于塑性极限分析上限解法分析了边坡不同失效模式的发生概率；Jiang等^[16]采用随机物质点法攻克了3维边坡大变形概率特征描述难题。随机极限平衡法（RLEM）也拓展到3维土坡^[17]、3维高聚物加固边坡^[18]、3维加筋边坡^[19]的失效概率分析。在既有方法中，3维RLEM与3维RFEM应用最为广泛，且两种方法在安全系数、滑动长度及滑动体积方面具有一致性^[20]。RLEM计算效率高，这一优势使得结合蒙特卡洛模拟（MCS）评估边坡失效概率成为可能，拓宽了其应用范围^[17]。

较2维方法，3维方法虽能提供更精确的结果，但计算效率的显著降低成为制约3维边坡失效概率分析的主要因素。Li等^[21]提出了逐步协方差矩阵分解方法（CMD），可使用较小的计算机内存实现大规模3维随机场的快速模拟^[22]。然而，除随机场生成效率外，3维RLEM中临界滑动面的搜索效率亦是影响计算效率的关键^[23]。传统试算法需遍历大量潜在滑动面以确定最小安全系数，导致效率与精度难以兼顾。为此，学者引入智能优化算法^[23‒25]，以高效搜索边坡最小安全系数，其中，粒子群优化（PSO）算法因其操作简便、高效稳定^[26]，被广泛应用于均质土坡^[27]、加筋土坡^[28]、地震力作用土坡^[29]、岩质边坡^[30]、3维不对称边坡^[31]等各类边坡分析。Wan等^[23]改进了PSO算法，将其应用于2维边坡失效概率计算中的最小安全系数搜索，显著提升了计算精度与效率。但该方法目前主要适用于2维圆弧滑动面的搜索，其在3维滑动面下的适用性与有效性仍有待进一步验证。

因此，本研究提出一种改进的CMD方法，生成大规模3维相关对数正态随机场，以期在提升方法的计算效率；结合改进PSO算法与简化Bishop法求解最小安全系数，并采用MCS评估边坡失效概率，揭示3维效应及土体空间变异性对边坡稳定性的影响机制。

1 方 法

1.1 随机场生成方法

对Li等^[21]提出的CMD方法进行改进，用于生成相关3维对数正态随机场。采用学者广泛使用的指数型相关函数描述土体的空间相关性\rho ({\tau }_x,{\tau }_y,{\tau }_z)^[2,22]。

\rho ({\tau }_x,{\tau }_y,{\tau }_z)=\mathrm{e}\mathrm{x}\mathrm{p}(-\left|{\tau }_x\right|/l_x-\left|{\tau }_y\right|/l_y-\left|{\tau }_z\right|/l_z)（1）

式中，τ_x 、τ_y 、τ_z 分别为任意两点在x、y、z方向的距离，l_x 、l_y 、l_z 分别为x、y、z方向的相关长度。

CMD方法通过Cholesky分解将协方差矩阵 C 转化为下三角矩阵 L，并结合独立随机变量ξ与 L 构造高斯随机场 H^[22]：

\mathit{H}=\mathit{L}\xi =({\mathit{L}}_{xz}\otimes {\mathit{L}}_y)\xi（2）

式中， L_xz 为xz平面内自相关矩阵 C_xz 的下三角矩阵， L_y 为y方向1维自相关矩阵 C_y 的下三角矩阵，⊗表示矩阵的克罗内克积运算。

土体黏聚力c和内摩擦角φ之间存在着一定的内在关联，忽视c和φ之间的关联会导致分析结果的不准确。因此，本研究进行了两项改进：其一，在生成标准正态随机场前通过互相关变换嵌入c和φ之间的相关系数；其二，结合拉丁超立方抽样以减少抽样方差并提升蒙特卡洛模拟的收敛效率。改进CMD方法模拟3维相关对数正态随机场的流程如下。

1）将随机场离散成若干单元，结合式（1）生成 C_xz 和 C_y，通过Cholesky分解生成 L_xz 和 L_y。

2）采用拉丁超立方抽样生成独立随机样本，并构成一个n×2维（n为随机场单元数目）的矩阵 ξ =[ ξ_c， ξ_φ ]， ξ_c =[ξ_c1,ξ_c2,…,ξ_cn ]^T， ξ_φ=[ξ_φ1,ξ_φ2,…,ξ_φn ]^T，ξ_ci 、ξ_φi 为服从标准正态分布的随机数。

3）对c和φ的互相关系数矩阵 ρ_c,φ 进行Cholesky分解得其上三角矩阵 L₁，相关标准正态随机样本矩阵 ξ^G=ξL₁=[ ξ_c^G， ξ_φ^G]其中， ξ_c^G=ξ_c， ξ_φ^G=ρ_c,φ × ξ_c +(1－ ρ{}_{c,\phi }^{\mathrm{0.5}}× ξ_φ ），代入式（3）即可得到相关高斯随机场{\mathit{H}}_R^{\mathrm{G}}。

{\mathit{H}}_R^{\mathrm{G}}=({\mathit{L}}_{xz}\otimes {\mathit{L}}_y){\mathit{\xi }}_R^{\mathrm{G}},R\in \{c,\phi \}（3）

式中，R包括黏聚力c和内摩擦角\phi两个随机场。

4）将相关高斯随机场等概率变换为相关对数正态随机场 H_R。

{\mathit{H}}_R=\mathrm{e}\mathrm{x}\mathrm{p}({\mu }_{\mathrm{l}\mathrm{n}R}+{\sigma }_{\mathrm{l}\mathrm{n}R}{H}_R^{\mathrm{G}}),R\in \{c,\phi \}（4）

式中，μ_{ln R} 、σ_{ln R} 分别为高斯随机场ln R的均值和方差。

以上步骤重复N次，可以得到N个3维相关对数正态随机场样本。

1.2 安全系数的计算方法

3维Bishop法在边坡稳定性分析中应用广泛，且与严格极限平衡法的计算结果十分接近^[32]。图1为3维简化Bishop法计算简图。

Bishop法将滑体划分成若干条柱，不考虑条间作用力（图1），通过整体力矩平衡计算安全系数F，计算式如下：

\begin{array}{l}F={\displaystyle \sum _{k=\mathrm{1}}^K}\left(c_k{A}_k+N_k\mathrm{t}\mathrm{a}\mathrm{n}{\phi }_k-u_k{A}_k\mathrm{t}\mathrm{a}\mathrm{n}{\phi }_k\right)\times (m_{zk}{x}_k-m_{xk}{z}_k)·\\ {\displaystyle \sum _{k=\mathrm{1}}^K}\left(W_k{x}_k-N_k\left(n_{zk}{x}_k-n_{xk}{z}_k\right)\right)\end{array}（5）

N_k=\frac{W_k-\left(c_k-u_k\mathrm{t}\mathrm{a}\mathrm{n}{\phi }_k\right)A_k{m}_{zk}/F}{n_{zk}+\mathrm{t}\mathrm{a}\mathrm{n}{\phi }_k{m}_{zk}/F}（6）

式（5）、（6）中，K为条分数目，c_k 、φ_k 分别为第k个随机场单元的土体的黏聚力和内摩擦角，A_k 为条柱底部面积，W_k 为条柱重力，u_k 为孔隙水压力，x_k 、z_k 分别为条柱底部中点到圆心的水平和竖直距离，(n_xk,n_yk,n_zk )、(m_xk,m_yk,m_zk )分别为条柱底部的法向作用力N_k 和切向作用力T_k 的方向余弦。

1.3 最小安全系数搜索方法

1.3.1 滑入滑出法

图2、3分别为滑动面形式示意图和3维潜在滑动面生成示意图。滑动面形状选取以下两种假定：1）滑动面为圆柱形状的一部分，且不考虑两端约束，简称“圆柱滑动面”（图2（a））；2）滑动面为椭球的一部分（两端）+圆柱形状的一部分（中间），简称“组合滑动面”（图2（b））。假定滑动面的宽度为B，高度为H，坡角为β，则组合滑动面通过滑入点D、滑出点E、水平向切线l₁和圆柱体宽度b确定（图3），图3（b）为图3（a）中的截面P‒P。将滑入点范围L₁等分成n₁份、滑出点范围L₂等分成n₂份、切线点范围L₃等分成n₃份、圆柱体宽度范围L₄等分成n₄份（图3），通过排列组合形成n₁×n₂×n₃×n₄个潜在滑动面。圆柱滑动面是2维滑动面在3维空间上的直接拓展，通常被认为是“假3维”，因此，无需对L₄进行等分，即n₄=1。

1.3.2 改进的粒子群优化算法

粒子群算法通过迭代计算搜索边坡的最小安全系数，每一次迭代需要根据个体和群体的经验更新粒子的位置和速度。粒子位置和速度的更新式如下：

v_{sj}(t+\mathrm{1})=c_{\mathrm{1}}{r}_{\mathrm{1}}(b_{sj}^{\mathrm{P}}(t)-x_{sj}(t))+c_{\mathrm{2}}{r}_{\mathrm{2}}(b_j^{\mathrm{G}}(t)-x_{sj}(t))+\omega v_{sj}(t)（7）

x_{sj}(t+\mathrm{1})=x_{sj}(t)+v_{sj}(t+\mathrm{1})（8）

式（7）、（8）中：v_sj (t+1)、x_sj (t+1)分别为粒子s第j个维度在第t+1次迭代时的速度和位置；t为迭代次数；c₁、c₂为学习因子，取值为c₁=c₂=2.0；ω为惯性权重，取值为0.7；r₁、r₂为[0,1]之间的随机数；b_s^{\mathrm{P}}(t)为第t次迭代时粒子s的历史最优位置；b_j^{\mathrm{G}}(t)为第t次迭代时第j个维度的全局最优位置。

使用安全系数F评价边坡的稳定性，当F<1.0时，边坡处于不稳定状态。边坡失效概率P_f计算式为N_f/N，其中，N_f和N分别为蒙特卡洛模拟中F<1.0发生的次数和模拟总次数。因此，计算失效概率只需判定F和1.0的关系，并不需要F的精确值。基于此，对PSO算法迭代终止条件进行改进，以期提高计算效率。改进后的收敛条件为：1）若b_j^{\mathrm{G}}(t)<1.0，结束迭代；2）若b_j^{\mathrm{G}}(t)>1.05且t\ge5，结束迭代；3）若b_j^{\mathrm{G}}(t)>1.02且t\ge10，结束迭代；4）若t\geT_max（T_max为最大迭代次数），结束迭代。

2 算例分析

2.1 验证

Cho^[2]提出一个纯黏性土边坡用于分析考虑土体参数空间变异性边坡的可靠度，后被许多研究者^{[5, 23, 33‒34]}用于验证所提边坡可靠度分析方法的合理性。本文将此算例拓展到3维，用于验证本文方法和程序的正确性。算例1的几何形状如图4所示，算例1、2的土体重度γ、黏聚力c的均值u_c 及其变异系数V_c 、内摩擦角φ的均值u_φ 及其变异系数V_φ 、c和φ的相关系数ρ_c，φ 等如表1所示。

不考虑土体参数的空间变异性，对此算例进行确定性分析，使用3维Bishop法结合PSO算法计算其安全系数。图5为安全系数F随破坏面宽度B的变化。

由图5可知，当滑动面为圆柱时，则该3维分析方法退化为2维，本文所得安全系数F=1.352 4和Wan等^[23]的结果F=1.352 5非常接近。当滑动面为圆柱和椭球组合时，F随破坏面宽度B的增大而减小，并逐渐趋近于2维分析结果。

考虑土体参数的空间变异性，使用改进的CMD方法生成随机场，结合3维Bishop法和改进粒子群算法搜索边坡的最小安全系数，通过10 000次蒙特卡洛模拟计算边坡的失效概率P_f。图6对比了不同l_y 下圆柱滑动面边坡失效概率P_f和平均安全系数F_av。当l_y 趋近于无穷大且滑动面为圆柱形时，则该3维分析退化为2维，不同宽高比（B/H=0.25、0.50、1.00、2.50、5.00、10.00、15.00、20.00、25.00、50.00）下的P_f位于[0.090 5，0.091 6]（图6），平均P_f为9.11%，与Wan等^[23]计算结果平均P_f（9.25%）非常接近，相对误差约为1.5%。安全系数和失效概率的对比分析结果证明了本文方法和程序的正确性。

2.2 土体空间变异性和3维效应对失效概率的影响

本文在算例1的基础上进行拓展，算例2与算例1几何形状相同，考虑不同的土体参数及其变异性，以及c和φ的相关性（表1）。将算例1由2维拓展到3维，水平方向增加一个参数l_y 以表征横向扩展尺度。参考3维边坡可靠度研究成果^{[5, 23, 33‒34]}，l_y 取值为20 m。特别地，l_y 趋近于无穷大时模型退化为2维平面应变状态，由此构建3维与2维分析的关联性对比，揭示土体空间变异性和3维效应对边坡失效概率的影响规律。使用本文所述方法计算算例1和算例2的失效概率。图7为不同l_y 下组合滑动面边坡失效概率和平均安全系数。

从图6可以看出，当滑动面为圆柱形且l_y 趋近于无穷大时，随破坏面宽度的增大，边坡的失效概率P_f和平均安全系数F_av基本保持不变。这是因为此种条件下，边坡的破坏面退化为2维且土体参数沿y轴分布均匀且恒定，使得边坡的稳定性分析降至2维范畴，安全系数的计算不再受破坏面宽度B的影响，进而B的变化对失效概率亦无影响。然而，当滑动面位圆柱形且l_y 为20 m时，随B的增大，P_f逐渐减小，F_av逐渐增大。这是因为当l_y 为20 m时，土体参数沿y轴存在变异性，随B的增大，任意一点的土体参数逐渐趋于平均值，相当于土体参数的空间变异性减小，因而P_f逐渐减小，F_av逐渐增大。图8为截面土体黏聚力平均值随破坏面宽度的变化，随着B的增大，截面不同位置土体的黏聚力均值c_u更加趋近于10 kPa。上述结果表明，土体参数的空间变异性是3维边坡稳定性高于2维分析的关键因素之一。众多研究明确指出，3维边坡之所以展现出更高的安全系数，主要归因于3维效应的作用。本文研究成果进一步补充并完善了3维边坡稳定性优于2维分析结果的机理，为理解和评估边坡稳定性提供了更为全面的视角。

从图7中可以看出：当l_y 趋近于无穷大时，随着B的增大，P_f逐渐增大，F_av逐渐减小，这说明B是影响3维边坡稳定性的重要因素；l_y 趋近于无穷大时的失效概率P_f均大于l_y 为20 m时的P_f，这说明了土体空间变异性是影响边坡失效概率另一重要因素；此外，l_y 为20 m时，失效概率随着B的增大呈现先增大后减小的趋势。这说明B较小时，较空间变异性的影响，3维效应的影响更加显著，随着B的增大，空间变异性对失效概率的影响逐渐增大。

2.3 变异系数对失效概率的影响

为研究参数变异性对边坡失效概率的影响，本小节分析V_c 、V_φ 及ρ_c,φ 对失效概率的影响。参考已有研究成果^{[3, 6]}，V_c 的取值范围为[0.1,0.6]，V_φ 的取值范围为[0.1,0.4]，ρ_c,φ 的取值范围为[－0.5,0.5]，算例2具体参数组合见表2。破坏面形式为组合滑动面，宽度设置为25 m。

图9～11分别为黏聚力变异系数、内摩擦角变异系数、黏聚力和内摩擦角的相关系数对失效概率和平均安全系数的影响。从图9～11可以看出，失效概率P_f随V_c 、V_φ 和ρ_c,φ 的增大而增大。V_c 从0.1增大到0.6时，其P_f从0增大到了0.250，增大了约0.250。V_φ 从0.1增大到0.4时，其P_f从0.040增大到了0.140，增大了约0.10。ρ_c,φ 从-0.5增大到0.5时，其P_f从0.055增大到了0.070，增大了约0.015。上述研究结果表明，V_c 、V_φ 和ρ_c,φ 对失效概率均有显著影响。

2.4 相关长度对失效概率的影响

参考已有研究成果^[1,6]，相关长度l_x 、l_y 、l_z 设置范围分别为[10，40]、[10，40]、[1，4] m，算例2具体参数组合见表3。破坏面形式为组合滑动面，宽度设置为25 m。

图12～14分别为x、y、z方向相关长度分别对失效概率的影响。

由图12～14可知：失效概率P_f随着相关长度的增大而增大。l_x 从10 m增大到40 m时，P_f增大了约0.030；l_y 从10 m增大到40 m时，P_f增大了约0.075；l_z 从1 m增大到4 m时，其P_f增大了约0.050。这一计算结果说明与l_x 相比，l_y 对P_f的影响较大。

因本文使用改进PSO算法计算安全系数，精度相对较差，尤其是安全系数小于1.0的情况，这会造成平均安全系数F_av计算结果的精度较差，其随相关长度的变化规律也不明显。

3 讨 论

3.1 粒子群的计算精度和计算效率

选择算例1中l_y 为20 m的3个典型随机场，采用组合滑动面假定，使用滑入滑出法和PSO算法搜索其最小安全系数，结果如表4所示。滑入滑出法设置n₁×n₂×n₃×n₄=41×41×21×21=741 321个滑动面。粒子群算法设置粒子个数M和迭代次数T_max的组合（M、T_max）为（20、60）、（30、40）、（40、30）、（20、24）、（60、20），潜在滑动面M×T_max=1 200个。

由表4可以看出，PSO算法进行1 200个潜在滑动面分析所得结果和滑入滑出法进行741 321个滑动面分析所得结果十分接近，PSO方法搜索滑动面的个数是滑入滑出法的0.16%，显著减少了计算量。结合Wan等^[23]的研究成果，本文选择（M、T_max）=（40、30）可以兼顾计算精度和计算效率。在此条件下，改进PSO算法平均搜索约250个滑动面，是粒子群算法的20%左右，显著提高了3维边坡失效概率的计算效率。

3.2 改进粒子群算法的终止条件

图15为PSO算法搜索最小安全系数时的迭代过程。从图15中可以看出，经过5次迭代与100次迭代安全系数的误差小于0.05；经过10次迭代和100次迭代的安全系数的误差在0.02以内。这一结果证明了本文改进迭代收敛条件的合理性。

4 结 论

本文使用改进的CMD方法生成大规模3维随机场，结合极限平衡法和改进粒子群算法搜索边坡的安全系数，采用蒙特卡洛模拟计算边坡的失效概率，分析了3维边坡的失效概率，得出以下结论：

1）改进粒子群算法可显著提高计算效率，为考虑土体参数空间变形边坡失效概率的计算提供了一个可靠高效的方法。

2）3维效应和空间变异性是影响3维边坡稳定性的两个主要因素。只考虑土体空间变异性时，失效概率随破坏面宽度的增大而减小；只考虑3维效应时，失效概率随破坏面宽度的增大而增大；两者同时考虑时，它们对失效概率的影响相互抵消，因而3维边坡的失效概率（最大为0.20%）远远低于2维边坡的失效概率（9.11%）。

3）考虑土体参数变异性和3维效应时，边坡的失效概率随破坏宽度的增大呈现先增大后减小的趋势，说明破坏面宽度较小时，3维效应起主要作用；破坏面宽度较大时，空间变异性起主要作用；3维效应和空间变异性对边坡稳定性的影响存在一个临界破坏面宽度。

参考文献

[1] El‒Ramly H, Morgenstern N R, Cruden D M.Probabilistic slope stability analysis for practice[J].Canadian Geotechnical Journal,2002,39(3):665‒683. doi:10.1139/t02-034

[2] Cho S E.Probabilistic assessment of slope stability that considers the spatial variability of soil properties[J].Journal of Geotechnical and Geoenvironmental Engineering,2010,136(7):975‒984. doi:10.1061/(asce)gt.1943-5606.0000309

[3] Wu Xingzheng, Jiang Liangwei, Luo Qiang,et al.Analysis of model uncertainty for stability reliability of embankment slope[J].Rock and Soil Mechanics,2015,36(Supp2):665-672.

[4] 吴兴正,蒋良潍,罗强,等.路堤边坡稳定可靠度计算中的模型不确定性分析[J].岩土力学,2015,36():665‒672.

[5] Bai Tao, Huang Xiaoming, Li Chang.Slope stability analysis considering spatial variability of soil properties[J].Journal of Zhejiang University (Engineering Science),2013,47(12):2221‒2226.

[6] 白桃,黄晓明,李昶.考虑土体参数空间变异性的边坡稳定性研究[J].浙江大学学报(工学版),2013,47(12):2221‒2226.

[7] Liu Leilei, Deng Zhiping, Zhang Shaohe,et al.Simplified framework for system reliability analysis of slopes in spatially variable soils[J].Engineering Geology,2018,239:330‒343. doi:10.1016/j.enggeo.2018.04.009

[8] Wang Yuke, Shao Linlan, Wan Yukuai,et al.Reliability analysis of earth rock dam slope considering spatial variability of soil parameters based on K‒L expansion method[J].Journal of Hydraulic Engineering,2024,55(12):1417‒1427.

[9] 王钰轲,邵琳岚,万愉快,等.基于K‒L展开法考虑土体参数空间变异性的土石坝坝坡可靠度分析[J].水利学报,2024,55(12):1417‒1427.

[10] Chen Zuyu, Mi Hongliang, Wang Xiaogang,et al.A three-dimensional limit equilibrium method for slope stability analysis[J].Chinese Journal of Geotechnical Engineering,2001,23(5):525‒529.

[11] 陈祖煜,弥宏亮,汪小刚,等.边坡稳定三维分析的极限平衡方法[J].岩土工程学报,2001,23(5):525‒529.

[12] Zhang Liang, Gong Wenping, Li Xinxin,et al.A comparison study between 2D and 3D slope stability analyses considering spatial soil variability[J].Journal of Zhejiang University:Science A,2022,23(3):208‒224. doi:10.1631/jzus.a2100139

[13] Wu Di, Wang Yuke, Chen Xin,et al.Limit analysis of 3d soil slopes considering the earthquake and nonlinear strength[J].Chinese Journal of Theoretical and Applied Mechanics,2024,56(5):1426‒1438.

[14] 吴迪,王钰轲,陈欣,等.考虑地震作用和非线性强度的三维土坡极限分析[J].力学学报,2024,56(5):1426‒1438.

[15] Griffiths D V, Huang Jinsong, Fenton G A.On the reliability of earth slopes in three dimensions[J].Proceedings of the Royal Society A:Mathematical,Physical and Engineering Sciences,2009,465(2110):3145‒3164. doi:10.1098/rspa.2009.0165

[16] Hicks M A, Spencer W A.Influence of heterogeneity on the reliability and failure of a long 3D slope[J].Computers and Geotechnics,2010,37(7/8):948‒955. doi:10.1016/j.compgeo.2010.08.001

[17] Hicks M A, Nuttall J D, Chen J.Influence of heterogeneity on 3D slope reliability and failure consequence[J].Computers and Geotechnics,2014,61:198‒208. doi:10.1016/j.compgeo.2014.05.004

[18] Kasama K, Furukawa Z, Hu Lihang.Practical reliability analysis for earthquake-induced 3D landslide using stochastic response surface method[J].Computers and Geotechnics,2021,137:104303. doi:10.1016/j.compgeo.2021.104303

[19] Varkey D, Hicks M A, Vardon P J.An improved semi-analytical method for 3D slope reliability assessments[J].Computers and Geotechnics,2019,111:181‒190. doi:10.1016/j.compgeo.2018.12.020

[20] Shu Shuang, Ge Bin, Wu Yongxin,et al.Probabilistic assessment on 3D stability and failure mechanism of undrained slopes based on the kinematic approach of limit analysis[J].International Journal of Geomechanics,2023,23:06022037. doi:10.1061/(asce)gm.1943-5622.0002635

[21] Jiang Shuihua, Li Jianping, Ma Guotao,et al.Probabilistic assessment of 3D slope failures in spatially variable soils by cooperative stochastic material point method[J].Computers and Geotechnics,2024,172:106413. doi:10.1016/j.compgeo.2024.106413

[22] Hu Lihang, Takahashi A, Kasama K.Effect of spatial variability on stability and failure mechanisms of 3D slope using random limit equilibrium method[J].Soils and Foundations,2022,62(6):101225. doi:10.1016/j.sandf.2022.101225

[23] Wang Yuke, Shao Linlan, Wan Yukuai,et al.Three-dimensional reliability stability analysis of earth-rock dam slopes reinforced with permeable polymer[J].Probabilistic Engineering Mechanics,2023,74:103537. doi:10.1016/j.probengmech.2023.103537

[24] Wang Yuke, Shao Linlan, Wan Yukuai,et al.Reliability analysis of three-dimensional reinforced slope considering the spatial variability in soil parameters[J].Stochastic Environmental Research and Risk Assessment,2024,38(4):1583‒1596. doi:10.1007/s00477-023-02636-5

[25] Lee S W, Ching J.Simplified risk assessment for a spatially variable undrained long slope[J].Computers and Geotechnics,2020,117:103228. doi:10.1016/j.compgeo.2019.103228

[26] Li Dianqing, Xiao Te, Zhang Limin,et al.Stepwise covariance matrix decomposition for efficient simulation of multivariate large-scale three-dimensional random fields[J].Applied Mathematical Modelling,2019,68:169‒181. doi:10.1016/j.apm.2018.11.011

[27] Ng C W W, Qu Chuanxiang, Ni Junjun,et al.Three-dimensional reliability analysis of unsaturated soil slope considering permeability rotated anisotropy random fields[J].Computers and Geotechnics,2022,151:104944. doi:10.1016/j.compgeo.2022.104944

[28] Wan Yukuai, Xu Renhao, Yang Rong,et al.Application of the improved particle swarm optimization method in slope probability analysis[J].Marine Georesources & Geotechnology,2024,42(10):1531‒1541. doi:10.1080/1064119x.2023.2291173

[29] Zhou X P, Huang X C, Zhao X F.Optimization of the critical slip surface of three-dimensional slope by using an improved genetic algorithm[J].International Journal of Geomechanics,2020,20(8):04020120. doi:10.1061/(asce)gm.1943-5622.0001747

[30] Sun Cong, Zheng Hong, Li Chunguang,et al.Critical slip surface search by using rigorous maximum residual thrust method based on ant colony algorithm[J].Rock and Soil Mechanics,2014,35(10):3021‒3026.

[31] 孙聪,郑宏,李春光,等.基于蚁群算法的严格最大剩余推力临界滑动面搜索[J].岩土力学,2014,35(10):3021‒3026.

[32] Cheng Y M, Li L, Chi S C.Performance studies on six heuristic global optimization methods in the location of critical slip surface[J].Computers and Geotechnics,2007,34(6):462‒484. doi:10.1016/j.compgeo.2007.01.004

[33] Kalatehjari R, Ali N, Hajihassani M,et al.The application of particle swarm optimization in slope stability analysis of homogeneous soil slopes[J].International Review on Modelling and Simulations,2012,5(1):458‒465.

[34] Shinoda M, Miyata Y.PSO‒based stability analysis of unreinforced and reinforced soil slopes using non-circular slip surface[J].Acta Geotechnica,2019,14(3):907‒919. doi:10.1007/s11440-018-0678-x

[35] Himanshu N, Burman A.Determination of critical failure surface of slopes using particle swarm optimization technique considering seepage and seismic loading[J].Geotechnical and Geological Engineering,2019,37(3):1261‒1281. doi:10.1007/s10706-018-0683-8

[36] Kumar V, Burman A, Himanshu N,et al.Rock slope stability charts based on limit equilibrium method incorporating Generalized Hoek‒Brown strength criterion for static and seismic conditions[J].Environmental Earth Sciences,2021,80(6):212. doi:10.1007/s12665-021-09498-6

[37] Kalatehjari R, Arefnia A, Rashid A S A,et al.Determination of three-dimensional shape of failure in soil slopes[J].Canadian Geotechnical Journal,2015,52(9):1283‒1301. doi:10.1139/cgj-2014-0326

[38] Hungr O, Salgado F M, Byrne P M.Evaluation of a three-dimensional method of slope stability analysis[J].Canadian Geotechnical Journal,1989,26(4):679‒686. doi:10.1139/t89-079

[39] Zhu Lei, Gao Xinyue, Wan Yukuai,et al.Effect of seismic force on slope reliability[J].Science Technology and Engineering,2023,23(10):4324‒4330.

[40] 朱磊,高欣悦,万愉快,等.地震力对边坡可靠度的影响[J].科学技术与工程,2023,23(10):4324‒4330.

[41] Wang Bin, Liu Leilei, Li Yuehua,et al.Reliability analysis of slopes considering spatial variability of soil properties based on efficiently identified representative slip surfaces[J].Journal of Rock Mechanics and Geotechnical Engineering,2020,12(3):642‒655. doi:10.1016/j.jrmge.2019.12.003

基金资助

国家自然科学基金项目(52408373)

宁夏自然科学基金优秀青年项目(2024AAC05036)

宁夏高等学校一流学科建设（水利工程）项目(NXYLXK2021A03)