考虑各向异性空间变异性抗滑桩加固边坡可靠度上限分析

Upper-bound Analysis of Slope Reliability Considering Anisotropic Spatial Variability and Reinforcement with Anti-slide Piles

摘要

土在自然沉积过程中其强度经常呈现出各向异性空间变异性特征，其中一般旋转各向异性在自然界土体中分布广泛，该特性体现土体在空间上力学特性的差异。为了考虑该特征对抗滑桩加固边坡稳定性的影响，首先，采用乔列斯基分解法，生成了具有代表性的横观各向异性和一般旋转各向异性的随机场；然后，基于离散机构上限分析理论，在空间变异土体中生成离散的速度间断面，使间断面严格满足塑性力学相关联流动法则，即空间上某处速度间断面的切向速度与真实速度的夹角等于该空间位置土体的内摩擦角，根据弹塑性分析理论求解抗滑桩侧极限土压力，继而建立抗滑桩加固边坡的能量平衡方程，通过强度折减法求解边坡破坏临界状态的安全系数；最后，从波动尺度、变异系数和互相关系数等方面，对比两种各向异性空间变异性随机场对抗滑桩加固边坡可靠度的影响。结果表明：一般旋转各向异性空间变异性桩加固边坡的安全系数均值低于横观各向异性桩加固边坡；一般旋转各向异性空间变异性桩加固边坡最优桩位X_f/L_x（桩安装位置到坡脚的距离与坡顶到坡脚距离的比值）趋近0.8，而横观各向异性桩加固边坡最优桩位X_f/L_x在0.6～0.7之间，这表明了在设计抗滑桩桩位时将土体的各向异性空间变异性纳入考量的必要性。

Abstract

Objective Accurately assessing the probability of slope failure induced by the anisotropic spatial variability of soil is crucial for slope reinforcement and risk evaluation. A current limitation in existing research is the insufficient discussion of the reinforcement effects of anti-slide piles and limit analysis under anisotropic random fields. This study introduces the discrete mechanism of limit analysis into general rotational anisotropic random fields and employs the upper-bound analysis method to evaluate slope reliability. Methods Firstly, the random field was generated using the Cholesky decomposition method. The autocorrelation function was transformed through a coordinate system conversion before the decomposition of the autocorrelation coefficient matrix to ensure that the random field demonstrated general rotational anisotropy. This procedure generated representative random fields for the soil strength parameters c and φ, which exhibited transverse isotropy and general rotational anisotropy, respectively. Then, the discrete mechanism upper-bound analysis theory was employed to generate velocity discontinuity surfaces in spatially variable soils, which rigorously satisfied the associated flow rule in plasticity mechanics. Specifically, at any spatial location, the angle between the tangential velocity along the velocity discontinuity surface and the actual velocity equaled the soil's internal friction angle. Next, to evaluate the reinforcement effect of anti-slide piles, the elastic-plastic analysis theory was applied to determine the ultimate lateral soil pressure acting on the piles. Assuming a linear distribution of lateral soil pressure along the piles, an energy balance equation was established for the slope reinforced with anti-slide piles. Finally, the strength reduction method was employed to calculate the safety factor at the critical state of slope failure. The influences of fluctuation scales, coefficients of variation, and correlation coefficients on the reliability of slopes reinforced with anti-slide piles were compared for the two types of anisotropic spatial variability random fields. Results and Discussions The verification of the calculation for the homogeneous pile-reinforced slope was conducted by setting all unit grid parameters of the random field to equal values, simulating a homogeneous condition. The safety factors for Examples 1 and 2 were calculated as 1.63 and 1.46, respectively, with an average error of 1.8% in the verification results. For the anisotropic spatially variable slope, the error in the failure probability was on the order of 10^-3. In addition, the variation trends of the failure probability for the anisotropic random field based on coordinate transformation remained consistent regardless of the differences between the minimum and maximum fluctuation angles. When analyzing the influence of fluctuation scales, the probability density function (PDF) curves of the safety factor became higher and narrower as the minimum fluctuation scale decreased. The influence of the minimum fluctuation scale on the PDF curve was more significant than that of the maximum fluctuation scale. As the minimum and maximum fluctuation scales increased, the PDF curves for general rotational anisotropy exhibited a greater tendency to become shorter and wider compared to those for transverse anisotropy. Regarding the influence of the coefficient of variation, as V_c and V_φ increased, the corresponding PDF curves became shorter and wider while shifting to the left, resulting in a significant increase in failure probability. As the coefficient of variation of the internal friction angle (V_φ ) increased, the distribution of the critical sliding surface of the slope became more dispersed. This result indicated that when the soil's V_φ value was high, ensuring slope stability could require increasing the burial depth of the anti-slide piles. Regarding the influence of the cross-correlation coefficient, the failure probability (P_f) increased as ρ_c,φ increased. Under identical cross-correlation coefficients, transverse anisotropic slopes exhibited greater reliability than general rotational anisotropic slopes. Regarding the installation position of the anti-slide piles, for transverse anisotropic random fields, P_f initially decreased as the height of the pile arrangement increased, before slightly increasing. For general rotational anisotropic random fields, P_f generally decreased as the height of the pile arrangement increased. The failure probability of the general rotational anisotropic spatially variable pile-reinforced slope was higher than that of the transverse anisotropic pile-reinforced slope. The optimal pile position X_f/L_x (the ratio of the distance from the pile installation position to the toe of the slope to the distance from the top of the slope to the toe of the slope) for the general rotational anisotropic spatially variable pile-reinforced slope approached 0.8. In contrast, for the transverse anisotropic pile-reinforced slope, it ranged between 0.6 and 0.7. When designing spatially variable pile-reinforced slopes, neglecting the influence of anisotropy on design parameters can lead to risks in determining the anti-slide pile installation position. Conclusions The results demonstrate that the optimal pile position is relatively insensitive to spatial variability parameters. However, a significant difference in the optimal pile position is observed between the two typical anisotropic random fields. In particular, under general rotational anisotropy conditions, the use of the X_f/L_x=0.6 pile position in practical engineering can pose potential risks. In addition, the internal friction angle of soil shows more significant fluctuations in spatial distribution. When the slope becomes unstable, the slip surface distribution tends to become more discrete and deeper. Therefore, increasing the embedding depth of anti-slide piles in engineering projects is essential to improve the reliability of slope stability.

Graphical abstract

关键词

土质边坡 / 各向异性空间变异性 / 蒙特卡洛模拟 / 上限分析 / 抗滑桩加固 / 离散机构

Key words

soil slope / anisotropic spatial variability / Monte Carlo simulation / upper-bound analysis / reinforcement with anti-slide piles / discrete mechanisms

本刊网刊

自然界中的土体在形成过程中普遍受到风化、侵蚀、搬运、沉积以及地质运动等多重因素的影响，使得看似均一的土层内部，在矿物组成、属性参数等方面呈现出明显的空间变异性^[1‒3]。这种变异性不仅是自然界复杂作用的直接体现，而且是工程实践中不可忽视的重要因素。在理论研究中，能使用随机场理论更准确地描述和模拟土体的空间变异性^[4]。为了具体实现随机场，常用的方法有局部平均法、乔列斯基分解法和级数展开法等。其中：局部平均法在处理非高斯随机场时，效果较差；级数展开法需要求解Fredholm积分方程，求解难度较高，特别是对于复杂形式的随机场函数；乔列斯基分解（Cholesky）法计算过程简单易于编程实现^[5]，推导过程更加具体、易于理解，尽管该方法在展开非平稳随机场方面尚存在局限性，但是本研究主要针对各向异性随机场的分布形式，同时，该方法在处理非高斯随机场（对数正态分布随机场）时等效自相关系数与原始自相关系数间的差别对计算结果影响较小，可忽略不计^[6]，所以，采用乔列斯基分解法对对数正态分布随机场进行展开。

如图1^[7]所示，自然界土体的空间变异性类型多样，除常见的各向同性、横观各向异性外，旋转各向异性及一般旋转各向异性同样是重要的变异形式，在实际工程中不可忽视。Liu等^[8‒9]提出，最大波动尺度和最小波动尺度的方向并不总是为水平和竖向。Chen等^[10]则基于此结论，采用随机有限差分法（RFDM）将各向异性分为水平情况和非水平情况，深入研究了一般旋转各向异性空间变异性对边坡的可靠度的影响。Qin等^[11]结合大漾云高速公路路基边坡的工程实际，对一般旋转各向异性空间变异性的多级边坡进行了有限元（FEM）计算，然而这一特性对桩加固边坡的影响尚需深入研究。

在滑动面或潜在滑动面以下的稳定地层中插入抗滑桩，桩体嵌入坚固地层，借助岩土体对桩的锚固约束作用，平衡上部滑体的下滑推力，从而稳定滑体或潜在滑体^[12‒13]。抗滑桩在岩土工程领域具有广泛的应用，特别是在山区、河岸等地质条件复杂、易发生滑坡的地区。通过合理设计和施工抗滑桩，可以显著减少滑坡及其次生地质灾害^[14]的发生，保护人民生命财产安全，维护工程设施的稳定运行。Sun等^[15]对横观各向异性随机场的桩加固边坡可靠度进行了参数分析，但对于一般旋转各向异性的研究尚显不足。宁帅等^[16]对各向异性和非均质性的抗滑桩边坡稳定性进行了分析，但其分析仍局限于确定性分析层面，尚未涉及到空间变异性的可靠度分析。Li等^[17]采用极限平衡法对旋转各向异性空间变异性桩加固边坡可靠度展开分析，但极限平衡需要进行条间力假设，无法考虑土体的塑性特征，因此，极限分析的相关研究尚需补充。

目前，较多边坡可靠度的研究都是采用数值分析的方法，计算过程耗时且对计算资源要求较高，并且Michalowski等^[18]指出数值方法更适合用于分析特殊且复杂的工况，而不利于进行广泛的参数分析。理论方法则主要分为极限平衡法和极限分析法，根据塑性力学的基本原理，边坡稳定性分析极限分析法通过构建机动速度场对边坡极限承载力进行求解，该方法成熟，理论基础严格，被广泛应用于非饱和边坡^[19]、裂隙边坡^[20‒21]、非均质边坡^[22‒23]的稳定性研究，为边坡稳定分析提供了坚实理论支撑^[24]。对于非均质边坡，孙志斌^[25]、孙志豪^[26]等提出了离散对数螺旋线的边坡上限分析机构，该机构对边坡滑移面进行了离散，使得速度间断面严格满足相关联流动法则，为空间变异性边坡的上限分析提供了新的解决思路。

上限分析离散机构已被广泛应用于空间变异性横观各向异性边坡的可靠度分析^[15,26]，但是尚未考虑一般旋转各向异性的影响，并且这一影响在数值方法和极限平衡法中进行了系统的讨论^[10,17]，因此，本文将一般旋转各向异性随机场引入极限分析上限法离散机构可靠度评估框架；另外，目前针对空间变异性边坡的加固方式及加固效果的研究较少，对于空间变异性边坡的治理与防护的认识和方法不够成熟，因此，本文针对抗滑桩的加固效果开展研究。首先，介绍了离散机构的生成过程及在考虑抗滑桩抗滑力作用下的功率计算方法。然后，利用乔列斯基分解法，对两种典型的各向异性随机场进行了离散。在蒙特卡洛模拟（MCS）过程中，使用MATLAB编写了上限分析安全系数求解程序。最后，进行参数分析，探讨了空间变异性及抗滑桩安装位置对结果的影响。

1 边坡稳定性的上限分析

1.1 离散机构

如图2所示的边坡上限分析离散机构中，AC为破坏机构的速度间断面，也可以被视为边坡的潜在滑动面，在经典旋转破坏机构中，该间断面为连续且光滑的对数螺旋线。假设速度间断面上下都是刚体，在这个设定中，滑动块体ABC绕点O旋转，其角速度为ω。为了方便进行后续的间断面离散，以坡脚C为原点建立直角坐标系。θ₀和θ_h分别为OC和OA与X轴负方向的夹角，r₀表示点C到旋转中心的半径，θ_i 为OP_i 与X轴负方向的夹角，φ_i 为点P_i 处土的内摩擦角，H为边坡高度， v_i 为点P_i 处的速度，P_i+1为求解点，i=1,2,3,…,n。

使用光滑的对数螺旋线时，由于整个空间的内摩擦角φ在空间上非均匀分布，因此存在局限性。离散机构采用点到点的方式构建速度间断面，将AC离散为一系列的单元 P₀P₁, P₁P₂,…， P_iP_i+1。此时，使每个间断面都满足相关联流动法则的要求，即速度间断面的切向速度与真实速度的夹角的大小等于土体的内摩擦角（即向量 P_i+1P_i 与点P_i 的速度 v_i 夹角为φ_i 的大小）。于是，已知点P_i 的坐标(x_i,y_i)时，可以按式（1）、（2）求解点P_i+1的坐标(x_{i+\mathrm{1}},y_{i+\mathrm{1}})：

\begin{array}{l}{x}_{i+\mathrm{1}}=x_i+\frac{\sqrt[]{(x_i-x_O{)}^{\mathrm{2}}+(y_i-y_O{)}^{\mathrm{2}}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Delta }\theta }{\mathrm{s}\mathrm{i}\mathrm{n}(\frac{\mathrm{\pi }}{\mathrm{2}}+{\phi }_i-\mathrm{\Delta }\theta )}\cdot \\ \mathrm{c}\mathrm{o}\mathrm{s}({\theta }_i-\frac{\mathrm{\pi }}{\mathrm{2}}+{\phi }_i)\end{array}（1）

\begin{array}{l}{y}_{i+\mathrm{1}}=y_i+\frac{\sqrt[]{(x_i-x_O{)}^{\mathrm{2}}+(y_i-y_O{)}^{\mathrm{2}}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\Delta }\theta }{\mathrm{s}\mathrm{i}\mathrm{n}(\frac{\mathrm{\pi }}{\mathrm{2}}+{\phi }_i-\mathrm{\Delta }\theta )}\cdot \\ \mathrm{s}\mathrm{i}\mathrm{n}({\theta }_i-\frac{\mathrm{\pi }}{\mathrm{2}}+{\phi }_i)\end{array}（2）

式（1）、（2）中：x_O 和y_O 分别为旋转中心点O的x坐标和y坐标；Δθ为OP_i 与OP_i+1夹角，离散机构的具体离散过程参见文献[25]，并将Δθ设定为1°。

1.2 能耗计算

桩加固边坡，仅考虑自重时，依据上限分析方法，如图3所示，破坏机构的能耗计算包括重力功率、抗滑桩提供的抗力做功功率，以及滑动面内部的能量耗散率3个方面^[27]。图3中，z代表距离地表的深度，p(z)表示桩身在深度z处所受侧向力的大小，R_i 为旋转中心O到P_i 点的距离，R_Gi 为旋转中心O到离散块体重心P_Gi 点的距离，θ_Gi 为OP_Gi 与X轴负方向的夹角。

其中：重力功率可以通过计算各离散单元的重力功率W之和来求解；滑动面的内部能量耗散率则可通过计算各离散单元的能量耗散率D_si 之和来求解。W和D_s的计算公式分别如下：

W=-\gamma \omega \sum _i{S}_i{R}_{Gi}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{Gi}（3）

D_{\mathrm{s}}=\omega \sum _i{c}_i{L}_i{R}_i\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_i（4）

式（3）、（4）中，γ为土体重度，S_i 为离散区域P_iBP_i+1的面积，L_i 为间断面P_iP_i+1的长度，c_i 为P_i 处的黏聚力。

抗滑桩引起的能量耗散与抗力的分布直接相关，Sun等^[15]在进行空间变异性桩加固边坡的可靠度分析中，根据弹塑性分析理论采用Ito等^[28]给出的公式求解了抗滑桩侧极限土压力大小（式（5））。本文也采用上述相同的方法，计算抗滑桩桩侧抗力p(z)。

\begin{array}{l}p(z)={(\frac{D_{\mathrm{1}}}{D_{\mathrm{2}}})}^{(N_{\phi (z)}^{\mathrm{1}/\mathrm{2}}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)+N_{\phi (z)}-\mathrm{1})}\left[\frac{\mathrm{1}}{N_{\phi (z)}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)}\left\{{\mathrm{e}}^{(\frac{D_{\mathrm{1}}-D_{\mathrm{2}}}{D_{\mathrm{1}}})N_{\phi (Z)}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)\mathrm{t}\mathrm{a}\mathrm{n}(\frac{\mathrm{\pi }}{\mathrm{8}}+\frac{\phi (z)}{\mathrm{4}})}-\mathrm{2}{N}_{\phi (z)}^{\mathrm{1}/\mathrm{2}}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)+N_{\phi (z)}-\mathrm{1}\right\}\right.+\\ \left.\frac{\mathrm{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)+\mathrm{2}{N}_{\phi (z)}^{\mathrm{1}/\mathrm{2}}+N_{\phi (z)}^{-\mathrm{1}/\mathrm{2}}}{N_{\phi (z)}^{\mathrm{1}/\mathrm{2}}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)+N_{\phi (z)}-\mathrm{1}}\right]\cdot c(z)D_{\mathrm{1}}-c(z)\left\{D_{\mathrm{1}}\frac{\mathrm{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)+\mathrm{2}{N}_{\phi (z)}^{\mathrm{1}/\mathrm{2}}+N_{\phi (z)}^{-\mathrm{1}/\mathrm{2}}}{N_{\phi (z)}^{\mathrm{1}/\mathrm{2}}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)+N_{\phi (z)}-\mathrm{1}}-\mathrm{2}{D}_{\mathrm{2}}{N}_{\phi (z)}^{-\mathrm{1}/\mathrm{2}}\right\}+\\ \frac{\gamma z}{N_{\phi (z)}}\left\{D_{\mathrm{1}}{(\frac{D_{\mathrm{1}}}{D_{\mathrm{2}}})}^{(N_{\phi (z)}^{\mathrm{1}/\mathrm{2}}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)-\mathrm{1})}\cdot {\mathrm{e}}^{\frac{D_{\mathrm{1}}-D_{\mathrm{2}}}{D_{\mathrm{1}}}{N}_{\phi (z)}\mathrm{t}\mathrm{a}\mathrm{n}\phi (z)\mathrm{t}\mathrm{a}\mathrm{n}(\frac{\mathrm{\pi }}{\mathrm{8}}+\frac{\phi (z)}{\mathrm{4}})}-D_{\mathrm{2}}\right\}\end{array}（5）

式中：c(z)和φ(z)分别为桩身在深度z处的平均黏聚力和摩擦力；D₁为桩间中心距；D₂为两桩外边缘之间的距离；N_φ(z)为与深度z处内摩擦角φ(z)相关的承载力系数，具体计算为：

N_{\phi (z)}=(\mathrm{t}\mathrm{a}\mathrm{n}(\mathrm{\pi }/\mathrm{4}+\phi {(z)/\mathrm{2}))}^{\mathrm{2}}（6）

因此，抗滑桩引起的能量耗散率可表示为：

D_{\mathrm{p}}=\omega {\int }_{\mathrm{0}}^h\frac{p(z)}{D_{\mathrm{1}}}l(z)\mathrm{d}z（7）

式中，l(z)为桩位于深度为z的点到旋转轴的垂直距离，h为滑移面与桩轴线的交点到地面的深度。

此时，滑动面能量耗散率与抗滑桩抗力做功功率之和D可表示为：

D=D_{\mathrm{p}}+D_{\mathrm{s}}（8）

1.3 边坡稳定性分析

在上限分析中，可通过强度折减法求解边坡的安全系数，即通过将土体强度指标折减F_s倍，使边坡处于临界破坏状态。该折减系数F_s为边坡的安全系数，且折减后的土体强度参数c_R和φ_R表示如下：

\left\{\begin{array}{l}{c}_{\mathrm{R}}=c/F_{\mathrm{s}},\\ {\phi }_{\mathrm{R}}=\mathrm{t}\mathrm{a}{\mathrm{n}}^{-\mathrm{1}}(\mathrm{t}\mathrm{a}\mathrm{n}\phi /F_{\mathrm{s}})\end{array}\right.（9）

对于边坡稳定性分析，本文采用安全系数优化方法^[26]，该方法适用于离散速度间断面的上限分析。在离散机构中，功率的计算是通过对离散单元求和得到的，因此，在计算过程中结合了二分法对安全系数进行优化^[26]。采用该方法，本文的安全系数优化流程图如图4所示。具体流程如下：

第1步：设置θ₀和r₀的搜索范围以及二分法中折减系数的范围[s₁,s₂]，其中，s₁和s₂分别代表二分法的初始下限和上限，本文取s₁=0.1，s₂=3，令折减系数F_s=(s₁+s₂)/2，代入式（9）求得折减后的土体强度参数。

第2步：在上限分析中，找到一组参数θ₀和r₀（使用序列二次规划法或全局法搜索），使f_min=min|W-D|<ε，其中，ε为预先设定的能量平衡收敛阈值。若F_s是安全系数的上限解，则而后令s₂=F_s；否则，说明在当前F_s下无法找到满足能量平衡的参数组合，边坡已处于破坏状态，令s₁=F_s。其中，阈值ε=2即可满足精度要求。

第3步：重复第1步和第2步，直到s₂-s₁<Δs，结束循环。其中，二分法区间收敛阈值Δs=0.01。

2 土体空间变异性和边坡可靠度分析

2.1 各向异性空间变异类型

Chen等^[10]对一般各向异性与一般旋转各向异性的边坡可靠度开展了研究。在一般各向异性中，最大波动尺度l_max固定于水平方向，最小波动尺度l_min与最大波动尺度l_max的夹角即最小波动尺度方向角θ_min在(0°,180°)区间变动；特别地，当θ_min=90°（即l_min与竖直方向重合）时，表现为横观各向异性。

在一般旋转各向异性中，l_max与l_min夹角在(0°,180°)区间变动，且l_max和l_min各自与水平方向的夹角也在变动，这种各向异性更能反映真实土壤的空间变异性，在实际土体中广泛存在。为了便于研究，Chen等^[10]固定l_min的方向与水平方向重合，令l_max与l_min的夹角即最大波动尺度方向角θ_max在(0°,180°)区间变动，本文延续这一思路，开展一般旋转各向异性的研究。

2.2 基于坐标转换的各向异性随机场生成

为了实现各向异性随机场，对笛卡儿直角坐标系进行转换，坐标变换示意图如图5所示，把点P的直角坐标(x,y)以平行投影的方式分解到斜坐标系X^{\text{'}}和Y^{\text{'}}坐标轴上，从而得到其在斜坐标轴上的斜坐标为(x^{\text{'}},y^{\text{'}})。

图5中，α为斜坐标轴X^{\text{'}}与直角坐标轴X的夹角，β为斜坐标轴X^{\text{'}}和Y^{\text{'}}的夹角。此时，直角坐标与斜坐标的转换关系表达式^[10]如下：

\left\{\begin{array}{l}{x}^{\text{'}}=(x\mathrm{s}\mathrm{i}\mathrm{n}(\alpha +\beta )-y\mathrm{c}\mathrm{o}\mathrm{s}(\alpha +\beta ))/\mathrm{s}\mathrm{i}\mathrm{n}\beta ,\\ y^{\text{'}}=(y\mathrm{c}\mathrm{o}\mathrm{s}\alpha -x\mathrm{s}\mathrm{i}\mathrm{n}\alpha )/\mathrm{s}\mathrm{i}\mathrm{n}\beta \end{array}\right.（10）

将式（10）中的斜坐标替换到自相关函数式（11）中的直角坐标位置，即通过相应的各向异性自相关函数\rho ({\tau }_x,{\tau }_y)来实现更多种类的各向异性随机场。

\rho ({\tau }_x,{\tau }_y)=\mathrm{e}\mathrm{x}\mathrm{p}[-\mathrm{2}(\frac{{\tau }_x}{{\delta }_{\mathrm{h}}}+\frac{{\tau }_y}{{\delta }_{\mathrm{v}}})]（11）

式中，τ_x 和τ_y 分别为空间任意两点在X轴方向和Y轴方向上的相对距离，δ_h和δ_v分别为波动尺度与X轴和Y轴方向重合。

根据式（11）中波动尺度δ_h和δ_v分别与X和Y轴方向重合的特性，通过调整式（10）改变X轴与Y轴的夹角来改变两个波动尺度之间的夹角。

采用乔列斯基分解法展开随机场，需要对n阶自相关系数矩阵进行乔列斯基分解，其中，n为随机场单元数量，也即模型划分单元网格的数量。然而，当随机场单元数量过多，例如，对于拥有4×10⁵个单元数量的三维随机场，自相关系数矩阵的分解需要大约1.28 TB的物理内存^[29]，这超出了绝大多数计算机的容量。针对这一挑战，Li等^[30]提出并完善了一种逐步分解的方法，将物理内存的占用从1.28 TB降低到5 GB，并且，在1 s内即可完成一次随机场的展开。在本文中，单元数量为1.2×10³，不进行逐步分解，物理内存占用约为4 GB，2×10^-1 s即可完成一次随机场的展开。综上所述，即使在高维场景下乔列斯基分解法也可以高效展开随机场。另外，该方法实现过程简单且不受具体随机场空间形态的限制，说明该方法具有较高的泛用性，适用于本文的研究。为了引入一般旋转各向异性，只需将建立的坐标系转换关系嵌入乔列斯基分解技术中即可。基于坐标转换的乔列斯基分解技术模拟相关非高斯各向异性随机场的流程概述如下。

将离散随机场空间划分为若干个单元，并通过转换关系将每个单元中心点的直角坐标P_i=(x_i,y_i)转换为斜坐标{P_i}^{\mathrm{\text{'}}}=(x_i^{\text{'}},y_i^{\text{'}})。

与简单随机抽样相比，拉丁超立方抽样（LHS）在相同的样本数量下，能够提供更全面的变量覆盖，从而更准确地反映变量的分布特性。故采用LHS从标准正态空间中抽取n个独立的土体强度参数c和φ的独立随机样本：

\mathit{\xi }=\left[{\mathit{\xi }}_c=({\xi }_{c,\mathrm{1}},{\xi }_{c,\mathrm{2}},\cdots ,{\xi }_{c,n}{)}^{\mathrm{T}},{\mathit{\xi }}_{\phi }=({\xi }_{\phi ,\mathrm{1}},{\xi }_{\phi ,\mathrm{2}},\cdots ,{\xi }_{\phi ,n}{)}^{\mathrm{T}}\right]（12）

互相关系数{\rho }_{c,\phi }反映了土体强度参数c和φ的相关程度，其取值范围为[-1,1]。当{\rho }_{c,\phi }为负时，表示负相关；相反地，当{\rho }_{c,\phi }为正时，表示正相关；当{\rho }_{c,\phi }=\mathrm{0}时，则表示不相关。

乔列斯基分解法展开随机场矩阵运算示意图，如图6所示。具体步骤如下。

第1步，对互相关系数矩阵{\mathit{\rho }}_{c,\phi }^{\mathrm{r}}进行Cholesky分解得到下三角矩阵 L₁，使得{\mathit{L}}_{\mathrm{1}}{\mathit{L}}_{\mathrm{1}}^{\mathbf{T}}={\mathit{\rho }}_{c,\phi }^{\mathrm{r}}。

第2步，使用Cholesky分解对自相关系数矩阵 C 进行分解得到下三角矩阵 L₂，使得{\mathit{L}}_{\mathrm{2}}{\mathit{L}}_{\mathrm{2}}^{\mathbf{T}}=C。

第3步， L₁左乘土体强度的独立随机样本 ξ，得到相关标准正态随机样本矩阵 ξ^D=ξL₁。

第4步， L₂右乘矩阵 ξ^D，得到相关高斯参数随机场{\mathit{\chi }}_{\mathrm{R}}^{\mathrm{G}}：

{\mathit{\chi }}_{\mathrm{R}}^{\mathrm{G}}(x^{\text{'}},y^{\text{'}})={\mathit{\chi }}_{\mathrm{R}}^{\mathrm{G}}(x,y)={\mathit{L}}_{\mathrm{2}}{\mathit{\xi }}^{\mathrm{D}}={\mathit{L}}_{\mathrm{2}}\mathit{\xi }{\mathit{L}}_{\mathrm{1}}（13）

{\mathit{\chi }}_{\mathrm{R}}^{\mathrm{G}}经等概率变换为相关对数正态参数随机场：

{\mathit{\chi }}_{\mathrm{R}}^{\mathrm{L}\mathrm{N}}(x^{\text{'}},y^{\text{'}})=\mathrm{e}\mathrm{x}\mathrm{p}({\mu }_{\mathrm{l}\mathrm{n}\mathrm{R}}+{\sigma }_{\mathrm{l}\mathrm{n}\mathrm{R}}{\mathit{\chi }}_{\mathrm{R}}^{\mathrm{G}}(x^{\text{'}},y^{\text{'}}))（14）

式中，μ_{ln R}和σ_{ln R}分别为高斯参数随机场ln c、ln φ的均值和标准差。

在横观各向异性随机场中，参数设置为α =0°，β =90°，δ_h=l_max，δ_v=l_min；一般旋转各向异性随机场中，参数设置为α=0°，β=150°，δ_h=l_min，δ_v=l_max。为了使随机场的土体参数真实地反映实际土壤参数分布情况，研究者们对现实土体进行了统计测量，并据此得到表1的土体参数统计特性^[5]。基于表1中的统计特性，得到随机场的一次典型实现，如图7所示。

2.3 可靠度分析方法

在对空间变异性边坡进行可靠度分析时，常采用失效概率作为评价指标^[10,26,31]，在选定特定的设计安全系数后，对生成的N个边坡随机场样本进行蒙特卡洛模拟，并通过上限分析方法计算出N个边坡安全系数F_s值。如果安全系数小于设计安全系数，则判定边坡失稳。当N足够大时，失稳边坡的数量占总样本数量N的比例即为边坡的失效概率。在边坡可靠度分析中，失效概率越大边坡可靠度越低，概率上更容易发生失稳。桩加固边坡的破坏属于小概率事件，由于样本数量受到计算效率的限制，当样本数量不够庞大时，采用直接蒙特卡洛的方法评估这种小概率事件会存在较大的统计不确定性。此外，评估小概率事件的相关统计方法在精度和泛用性上也存在一定的局限性。顾鑫等^[32]在土质边坡可靠度分析中，参照相关规范选取了设计安全系数，并在有限的蒙特卡洛模拟次数下有效评估了水位对边坡可靠度的影响。本文根据《建筑边坡工程技术规范》（GB50330—2013）^[33]，取F_s=1.15为设计安全系数，以平衡计算精度与计算效率。

为更准确估计边坡的失效概率并提高计算效率，在不同黏聚力c变异系数（V_c ）下观测随机统计量的收敛情况，如图8所示，确定N=1 500，即认为在设计安全系数为1.15时，1 500次MCS所统计的失效概率P_f是稳定的。

在本文的可靠度分析中，首先，根据边坡的几何形状和土体的空间变异特征生成N个随机场（在本文中N=1 500）。然后，联合使用强度折减法与二分法对每个随机场求解边坡的安全系数。最后，对所得到的结果采用直接蒙特卡洛方法计算边坡的失效概率P_f。

2.4 算例验证

2.4.1 均质桩加固边坡验证

对均质桩加固边坡进行计算验证时，使随机场的每个单元网格参数相等，以模拟均质状态。将本文方法与已有的算例^[34‒36]所得的均质桩加固边坡安全系数进行对比，结果如表2所示。

由表2可见：算例1中，Li等^[34]采用了水平方向抗滑力分布模式，其中，桩的阻力呈三角形分布，并通过极限分析得出安全系数F_s=1.70。谭捍华等^[35]额外考虑了有效土压力弯矩效应，得出算例1和2的安全系数分别为1.63和1.46。Won等^[36]则考虑了桩土耦合效应的影响，并计算出算例2的安全系数为1.49。本文计算结果与已有解答非常接近。

2.4.2 空间变异性边坡验证

选取边坡高度H=10 m，坡角β_s=45°；随机场参数中，摩擦角均值φ_mu=30°，波动尺度l_min=4 m、l_max=40 m，互相关系数ρ_c,φ =0；其他参数与表1保持一致。在一般各向异性和一般旋转各向异性下，将本文方法与Chen等^[10]采用随机有限差分法的边坡失效概率结果进行对比验证，结果如图9所示。由图9可以看出，本文的计算得到的失效概率与Chen等^[10]所得失效概率基本吻合，误差在10^-3数量级以内。此外，基于坐标转换的各向异性随机场在θ_min和θ_max不同时失效概率的变动趋势也相同。

此外，如图9（a）所示，一般各向异性中，当θ_min=90°时，边坡的失效概率达到最高，为0.069。Chen等^[10]也指出，横观各向异性（θ_min=90°）是一般各向异性中最危险的情况，因此，横观各向异性是本文所考虑的第1种随机场。

如图9（b）所示，一般旋转各向异性中，当θ_max=150°时，边坡的失效概率最高，为0.151。Chen等^[10]也指出，当l_max与l_min的夹角为150°时，边坡处于最危险的状态。从边坡最危险的角度出发，本文考虑θ_max=150°这一特定夹角下的一般旋转各向异性（后文中都简称为一般旋转各向异性），因此，一般旋转各向异性是本文所考虑的第2种随机场。

3 算例分析

分析了横观各向异性随机场和一般旋转各向异性随机场的不同参数对桩加固边坡稳定性的影响。这些参数包括随机场参数以及抗滑桩安装位置等。分析过程中，考虑相应参数的变化，同时保持其余参数不变，具体数值参照表1，其中边坡高度H=10 m，坡角β_s=45°。

3.1 波动尺度的影响

波动尺度表示了随机场土体强度参数的相关性。在实际中，最大波动尺度的范围为10～30 m，最小波动尺度的范围为1～3 m^[37]。图10（a）、（b）分别展示了最小、最大波动尺度对F_s分布的影响，其中，实线代表横观各向异性的情况，虚线则代表一般旋转各向异性的情况。

由图10可见：在两种各向异性随机场中，随着最小波动尺度的减小，概率函数密度（PDF）曲线都呈现出更高更窄的趋势，而且，相较于最大波动尺度，最小波动尺度对PDF曲线的影响更为显著。从物理的角度解释，最小波动尺度的改变直观地体现在土体强度参数的分层情况上：当最小波动尺度减小时，在对应方向上，若边坡尺寸保持不变，则土的层数会增加，每一层土的厚度会相应减小，这时边坡的稳定性会趋于增强。而最大波动尺度则主要体现在每一层土的分布情况中：当最大波动尺度增大时，每层土更趋于均质状态；尤其是当最大波动尺度和边坡尺寸相近时，这种趋于均质化的趋势则更加明显。

观察图10（a），随着最小波动尺度的增大，一般旋转各向异性的PDF曲线趋于矮胖的程度要大于横观各向异性的情况；同样地，如图10（b）所示，随着最大波动尺度的增大，一般旋转各向异性的PDF曲线趋于矮胖的程度也大于横观各向异性的情况。特别是在横观各向异性的情况下，最大波动尺度的影响相对较小；而在一般旋转各向异性的情况下，最大波动尺度的影响才更加显著。

3.2 变异系数的影响

图11绘制了不同土体强度指标变异系数下F_s的PDF曲线，实线为横观各向异性，虚线为一般旋转各向异性。由图11（a）可见，随着V_c 的增大，相应的PDF曲线变得矮胖，并向左偏移，失效概率明显增大。这是由于在V_c 增大的情况下，若摩擦角φ变异系数V_φ 保持不变，则土体的总抗剪强度方差也会增大，安全系数的方差增大而均值减小。同样地，如图11（b）所示，随着V_φ 的增大，也出现相同的变化趋势，根本原因与V_c 变化时的原因一致。

按照这一趋势，当V_φ 增大时，边坡的可靠度会降低，即V_φ=0.15时的可靠度低于V_φ =0.10时的可靠度。而在图15（b）中，对于一般旋转各向异性的情况，V_φ =0.10的曲线和在横观各向异性情况下的V_φ =0.15的曲线基本重合。因此，如果只考虑横观各向异性的空间变异性而未考虑一般旋转各向异性的空间变异性，可能会高估桩加固边坡的可靠度水平。

在离散滑移面迭代生成过程中，向量 P_i+1P_i 与点P_i 的速度 v_i 的夹角等于该点处土体的内摩擦角。均质边坡与空间变异边坡滑移面对比如图12所示。

由图12可见，均质边坡内摩擦角为10°，空间变异边坡内摩擦角均值虽同为10°，但由于其随空间位置的变动呈现随机波动特性，滑移面轨迹有所偏移且不光滑。值得注意的是，在不同的随机场中，滑移面的分布也具有随机性，这充分体现了内摩擦角空间变异性对滑移面形状具有显著影响。不同内摩擦角变异系数条件下边坡临界滑移面分布如图13所示。

由图13可见，随着内摩擦角的变异系数V_φ 增大，边坡临界滑移面的分布会更加离散。这表明，当土体V_φ 值较大时，为确保边坡的稳定性，需要再增加抗滑桩的埋置深度。

3.3 互相关系数的影响

在自然界中，土体的黏聚力和摩擦角具有一定相关性。为表征土体强度参数之间的相关性，常采用互相关系数{\rho }_{c,\phi }进行量化。Lumb^[38]提供了一系列测量的互相关系数数据，指出在实际土体中，互相关系数为负数的情况所占比例更多，例如，Yucemen^[39]给出互相关系数的变化范围为-0.49～-0.24。另外，Wolff^[40]在土体固结排水试验中得出的互相关系数为0.25。图14展示了当桩位X_f/L_x=0.6时失效概率P_f随互相关系数变化的关系。其中，其他土体参数参考表1。

从图14中可以发现，P_f随着ρ_c,φ 的增大而增大。在互相关系数相同的情况下，横观各向异性边坡比一般旋转各向异性边坡更可靠。

3.4 抗滑桩安装位置的影响

桩的安装位置是影响加固效果的重要设计参数。不同波动尺度、不同土体强度指标变异系数下桩的位置X_f/L_x与失效概率P_f之间的关系分别如图15、16所示，图15、16中，实线为横观各向异性，虚线为一般旋转各向异性。

对于横观各向异性随机场，随着桩位布置高度的提升，P_f先降低，而后略有增加。对于一般旋转各向异性随机场，随着桩位布置高度的提升，P_f则基本呈现下降趋势，但存在一个例外情况，在图16（a）中，当摩擦角变异系数较小时，P_f也表现出先降低后略有增加的趋势。

从P_f的最小值所对应的桩位可以看出，在一般旋转各向异性情况下，最安全桩位更靠近坡顶。在进行空间变异性桩加固边坡的设计时，如果仅仅考虑随机场参数，而忽视各向异性对设计的影响，那么可能会给抗滑桩安装位置的设计带来风险。

4 结 论

为了在抗滑桩加固边坡稳定性分析中考虑土体的空间变异性，本文采用了前人提出的上限分析离散机构方法，并将其与乔列斯基分解法相结合来展开随机场。这种方法在生成速度间断面时能考虑土体内摩擦角随空间改变的影响，并同时考虑了两种典型的空间各向异性情况。在进行边坡可靠度分析时，得到了边坡稳定性与土体强度、空间变异性和抗滑桩安装位置之间的关系。主要结论如下：

1）土体空间变异性的3个参数（波动范围、互相关系数和变异系数）对桩加固边坡稳定性的影响显著。在两种典型的旋转各向异性随机场中，P_f都随这3个参数的增大而增大。土体强度的波动尺度和变异系数越小，得到的F_s值离散性越小。此外，最小波动尺度对F_s分布的影响大于最大波动尺度，V_c 的影响比V_φ 更显著。

2）桩体设计参数X_f/L_x对桩加固边坡的P_f有较大的影响。将桩安装在适当位置可以有效地提高桩对边坡的承载力。虽然最优的桩体位置受空间变异性参数的影响较小，但是在两种典型各向异性随机场中，最优桩体位置却出现明显不同。特别是在一般旋转各向异性的情况下，实际工程中采取X_f/L_x=0.6的桩位可能会带来风险。

参考文献

[1] Dasaka S M, Zhang L M.Spatial variability of in situ weathered soil[J].Géotechnique,2012,62(5):375-384. doi:10.1680/geot.8.p.151.3786

[2] Mitchell J K, Soga K.Fundamentals of soil behavior[M].Hoboken:John Wiley & Sons,2005.

[3] Liu Dongsheng, Liu Hanlong, Wu Yue,et al.Genetic features of geo-materials and their testing method[J].Journal of Civil and Environmental Engineering,2022,44(4):1-9. doi:10.11835/j.issn.2096-6717.2021.183

[4] 刘东升,刘汉龙,吴越,等.岩土材料的基因特征及其测试方法[J].土木与环境工程学报(中英文),2022,44(4):1-9. doi:10.11835/j.issn.2096-6717.2021.183

[5] Vanmarcke E H.Probabilistic modeling of soil profiles[J].Journal of the Geotechnical Engineering Division,1977,103(11):1227-1246. doi:10.1061/ajgeb6.0000517

[6] Jiang Shuihua, Li Dianqing, Zhou Chuangbing,et al.Slope reliability analysis considering effect of autocorrelation functions[J].Chinese Journal of Geotechnical Engineering,2014,36(3):508-518. doi:10.11779/CJGE201403014

[7] 蒋水华,李典庆,周创兵,等.考虑自相关函数影响的边坡可靠度分析[J].岩土工程学报,2014,36(3):508-518. doi:10.11779/CJGE201403014

[8] Cho S E, Park H C.Effect of spatial variability of cross-correlated soil properties on bearing capacity of strip footing[J].International Journal for Numerical and Analytical Methods in Geomechanics,2010,34(1):1-26. doi:10.1002/nag.791

[9] Zhu H, Zhang L M.Characterizing geotechnical anisotropic spatial variations using random field theory[J].Canadian Geotechnical Journal,2013,50(7):723-734. doi:10.1139/cgj-2012-0345

[10] Liu W F, Leung Y F.Characterising three-dimensional anisotropic spatial correlation of soil properties through in situ test results[J].Géotechnique,2018,68(9):805-819. doi:10.1680/jgeot.16.p.336

[11] Liu W F, Leung Y F.Spatial variability of saprolitic soil properties and relationship with joint set orientation of parent rock:Insights from cases in Hong Kong[J].Engineering Geology,2018,246:36-44. doi:10.1016/j.enggeo.2018.09.015

[12] Chen Longlong, Zhang Wengang, Chen Fuyong,et al.Probabilistic assessment of slope failure considering anisotropic spatial variability of soil properties[J].Geoscience Frontiers,2022,13(3):101371. doi:10.1016/j.gsf.2022.101371

[13] Qin Hui, Yin Xiaotao, Tang Hua,et al.Reliability analysis and geometric optimization method of cut slope in spatially variable soils with rotated anisotropy[J].Engineering Failure Analysis,2024,158:108019. doi:10.1016/j.engfailanal.2024.108019

[14] Peng Yipu, Li Jian, Huang Yi,et al.Three-dimensional deformation monitoring and numerical simulation on slope excavation and support[J].Journal of Shenyang Jianzhu University Natural Science,2023,39(6):961-969.

[15] 彭仪普,李剑,黄毅,等.边坡开挖支护三维变形监测与数值模拟研究[J].沈阳建筑大学学报(自然科学版),2023,39(6):961-969.

[16] Liu Jing.Application of micro-cluster steel pipe piles for foundation pit support of landslide mass engineering[J].Journal of Shenyang Jianzhu University Natural Science,2015,31(6):1031-1040. doi:JournalArticle/5b3b5e61c095d70f006f94e7

[17] 刘敬.微型集群钢管桩在滑坡体基坑支护工程中的应用[J].沈阳建筑大学学报(自然科学版),2015,31(6):1031-1040. doi:JournalArticle/5b3b5e61c095d70f006f94e7

[18] Wu Fengyuan, Qin Pengfei, Li Xinchao,et al.Process simulation and characteristics analysis of long run-out debris flow movement based on material point method[J].Journal of Shenyang Jianzhu University Natural Science,2024,40(1):86-95.

[19] 吴凤元,秦鹏飞,李鑫超,等.基于物质点法的远程滑坡碎屑流运动过程模拟及特征分析[J].沈阳建筑大学学报(自然科学版),2024,40(1):86-95.

[20] Sun Zhibin, Huang Guoxian, Hu Yining,et al.Reliability analysis of pile-reinforced slopes in width-limited failure mode considering three-dimensional spatial variation of soil strength[J].Computers and Geotechnics,2023,161:105528. doi:10.1016/j.compgeo.2023.105528

[21] Ning Shuai, Zhuang Yan, Cui Xiaoyan,et al.Stability analysis of slope reinforced with piles in anisotropic and heterogeneity soils[J].Journal of Harbin Engineering University,2022,43(4):529-535.

[22] 宁帅,庄妍,崔晓艳,等.各向异性和非均质性的抗滑桩边坡稳定性分析[J].哈尔滨工程大学学报,2022,43(4):529-535.

[23] Li Jingze, Zhang Shaohe, Liu Leilei,et al.Probabilistic analysis of pile-reinforced slopes in spatially variable soils with rotated anisotropy[J].Computers and Geotechnics,2022,146:104744. doi:10.1016/j.compgeo.2022.104744

[24] Michalowski R L, Martel T.Stability charts for 3D failures of steep slopes subjected to seismic excitation[J].Journal of Geotechnical and Geoenvironmental Engineering,2011,137(2):183-189. doi:10.1061/(asce)gt.1943-5606.0000412

[25] Huang Wengui, Leong E C, Rahardjo H.Upper-bound limit analysis of unsaturated soil slopes under rainfall[J].Journal of Geotechnical and Geoenvironmental Engineering,2018,144(9):04018066. doi:10.1061/(asce)gt.1943-5606.0001946

[26] Gao Yufeng, Song Wenzhi, Zhang Fei,et al.Limit analysis of slopes with cracks:Comparisons of results[J].Engineering Geology,2015,188:97-100. doi:10.1016/j.enggeo.2015.01.013

[27] Tang Gaopeng, Zhao Lianheng, Li Liang,et al.Stability charts of slopes under typical conditions developed by upper bound limit analysis[J].Computers and Geotechnics,2015,65:233-240. doi:10.1016/j.compgeo.2014.12.008

[28] Wang Zhen, Cao Lanzhu, Wang Dong.Evaluation on upper limit of heterogeneous slope stability[J].Rock and Soil Mechanics,2019,40(2):737-742. doi:10.16285/j.rsm.2017.1276

[29] 王珍,曹兰柱,王东.非均质边坡稳定性上限分析评价研究[J].岩土力学,2019,40(2):737-742. doi:10.16285/j.rsm.2017.1276

[30] Luan Maotian, Tingkai Nian, Yang Qing.Stability analysis of pile-stabilized slopes considering both nonhomogeneity and anisotropy of soil strength using upper bound method of limit analysis[J].Rock and Soil Mechanics,2006,27(4):530-536. doi:10.1016/S1872-1508(06)60035-1

[31] 栾茂田,年廷凯,杨庆.考虑非均质各向异性效应的阻滑桩加固土坡稳定性分析[J].岩土力学,2006,27(4):530-536. doi:10.1016/S1872-1508(06)60035-1

[32] Liu Fengtao, Zhang Shaofa, Dai Beibing,et al.Upper bound limit analysis of soil slopes based on rigid finite element method and second-order cone programming[J].Rock and Soil Mechanics,2019,40(10):4084-4091. doi:10.16285/j.rsm.2019.0226

[33] 刘锋涛,张绍发,戴北冰,等.边坡稳定分析刚体有限元上限法的锥规划模型[J].岩土力学,2019,40(10):4084-4091. doi:10.16285/j.rsm.2019.0226

[34] Sun Zhibin, Pan Qiujing, Yang Xiaoli,et al.Discrete mechanism for upper bound analysis of nonhomogeneous slopes[J].Chinese Journal of Rock Mechanics and Engineering,2017,36(7):1680-1688. doi:10.13722/j.cnki.jrme.2016.1565

[35] 孙志彬,潘秋景,杨小礼,等.非均质边坡上限分析的离散机构及应用[J].岩石力学与工程学报,2017,36(7):1680-1688. doi:10.13722/j.cnki.jrme.2016.1565

[36] Sun Zhihao, Tan Xiaohui, Sun Zhibin,et al.Reliability of spatially variable earth slopes based on the upper bound analysis[J].Rock and Soil Mechanics,2021,42(12):3397-3406.

[37] 孙志豪,谭晓慧,孙志彬,等.基于上限分析的空间变异土质边坡可靠度[J].岩土力学,2021,42(12):3397-3406.

[38] Chen W F.Limit analysis and soil plasticity[M].Amsterdam:Elsevier,1975:99-106.

[39] Ito T, Matsui T.Methods to estimate lateral force acting on stabilizing piles[J].Soils and Foundations,1975,15(4):43-59. doi:10.3208/sandf1972.15.4_43

[40] 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

[41] 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

[42] Zhu Desheng, Xia Lei, Ke Lijun,et al.Reliability analysis of infinite slopes subjected to seismic loadings considering spatially variable shear strength parameters[J].Advanced Engineering Sciences,2025,57(3):106-114.

[43] 朱德胜,夏磊,柯力俊,等.考虑土体参数空间变异性的地震作用下无限长边坡可靠度研究[J].工程科学与技术,2025,57(3):106-114.

[44] Gu Xin, Zhang Wengang, Qiang Ou,et al.Reliability analysis of soil slope stability based on Chebyshev-Galerkin-KL expansion[J].Chinese Journal of Geotechnical Engineering,2023,45(12):2472-2480.

[45] 顾鑫,仉文岗,欧强,等.基于Chebyshev-Galerkin-KL展开的土质边坡稳定可靠度分析[J].岩土工程学报,2023,45(12):2472-2480.

[46] 中华人民共和国住房和城乡建设部. 建筑边坡工程技术规范:GB 50330—2013 [S].北京:中国建筑工业出版社,2013.

[47] Li X P, He S M, Wang C H.Stability analysis of slopes reinforced with piles using limit analysis method[C]//Proceedings of the Advances in Earth Structures.Shanghai:American Society of Civil Engineers,2006:105-112. doi:10.1061/40863(195)8

[48] Tan Hanhua, Zhao Lianheng, Li Liang,et al.Energy analysis method for pre-reinforcing slopes with anti-slide piles[J].Rock and Soil Mechanics,2011,32(Supp2):190-197.

[49] 谭捍华,赵炼恒,李亮,等.抗滑桩预加固边坡的能量分析方法[J].岩土力学,2011,32():190-197.

[50] Won J, You K, Jeong S,et al.Coupled effects in stability analysis of pile-slope systems[J].Computers and Geotechnics,2005,32(4):304-315. doi:10.1016/j.compgeo.2005.02.006

[51] El-Ramly H, Morgenstern N R, Cruden D M.Probabilistic stability analysis of a tailings dyke on presheared clay-shale[J].Canadian Geotechnical Journal,2003,40(1):192-208. doi:10.1139/t02-095

[52] Lumb P.Safety factors and the probability distribution of soil strength[J].Canadian Geotechnical Journal,1970,7(3):225-242. doi:10.1139/t70-032

[53] Yucemen M S.A probabilistic study of safety and design of earth slopes[M].Twin Cities Urbana:University of Illinois at Urban-Champaign,1973:89-117.

[54] Wolff T F.Analysis and design of embankment dam slopes:A probabilistic approach[M].West Lafayette:Purdue University,1985:124-131.

基金资助

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

辽宁省自然科学基金计划项目（2022‒BS‒195）

教育部“春晖计划”合作科研项目(HZKY20220416)

辽宁省教育厅高校基本科研项目(JYTMS20231608)

沈阳市优秀中青年科技人才项目(RC230805)