土质边坡稳定性分析是岩土结构设计重要内容之一。传统上采用稳定安全系数F_s表达的确定性分析方法进行分析，F_s依据极限平衡法或有限元法计算。确定性分析方法简便直观，但不能考虑边坡中诸多不确定性因素的影响，因而，某一特定的F_s的大小并不能表征边坡的安全程度。因此，将以随机变量为研究对象，考虑各种工程设计因素中的不确定性，以可靠指标β或失效概率P_f为表达的概率分析方法被引入边坡工程。

概率分析常采用1次2阶矩法^[1-2]、1阶可靠度法^[3-5]、基于Monte Carlo模拟的随机有限元法^[6-7]等方法，存在计算复杂和耗时较长的不足。近年来，为提高边坡可靠度计算效率，学者们提出子集模拟法^[8-11]、非侵入式随机有限元法^[12]、代理模型方法^[13-14]、稳定性分析图法等改进方法。子集模拟法采用Bayes公式，将小概率失效事件的P_f表达为统一失效模式下具有较大P_f的多个中间失效事件的条件概率乘积，提高了计算效率，但在求解多失效模式可靠度时需对每个模式单独计算，即重复执行随机模拟^[15]；非侵入式随机有限元法采用非侵入式理念将有限元法与随机响应面法结合，无需修改有限元软件源代码，但在运用随机响应面法对F_s随机多项式展开的过程中，需根据有限元法得到F_s线性方程组，并求解获得随机多项式展开系数，耗时较长的问题依然存在；代理模型方法需要在变量空间中随机采样，计算量较大，且该方法在基于极限平衡方法进行边坡可靠度计算时，需要搜索临界滑裂面来确定边坡稳定性系数，这亦会增加额外计算量；稳定性分析图法建立了F_s 与P_f关系曲线，采用图表法获取P_f，如：Javankhoshdel和Bathurst^[16]考虑土体重度γ变异性，构建了饱和黏性土（ $\varphi$ =0）与一般黏性土边坡的P_f ～ F_s关系曲线，简化了边坡稳定性的概率分析步骤，但需采用极限平衡法与Monte Carlo模拟（MCS）重复计算，获得不同边坡参数组合下的F_s及P_f，最后建立无量纲参数条件下的P_f ～ F_s曲线，在工程中推广应用较复杂。

由于稳定可靠指标β能考虑边坡土体强度参数不确定性，相比F_s更为全面，但所需数据量和计算量较大，并缺乏F_s所具有的丰富的工程实际经验，制约了实际应用。因此，在边坡几何、土体强度参数均值及变异系数等数据已知的条件下，研究β的快速获得方法，对指导工程设计具有重要意义。

为此，针对工程条件多样性特点，基于常遇坡高下1∶1.5典型坡率的简单土质边坡模型与土体强度参数变异水平，制定以坡高和土体强度参数为主要因素的稳定分析模型及计算方案；运用Fellenius极限平衡法及Monte Carlo模拟，开展不同方案下可靠指标计算，讨论安全系数与可靠指标的基本关系及变化规律；基于模型参数与计算数据，运用无量纲化与多元非线性回归技术，构建基于稳定安全系数的边坡可靠指标估算方程；分析回归系数与坡高、土体强度参数变异系数关系特征，讨论估算偏差随安全系数和土体强度的变化趋势。研究成果对完善路基工程设计的边坡稳定可靠性评价具有参考价值。

1.   稳定分析模型及计算方案

1.1   安全系数与可靠指标

边坡稳定安全系数F_s通常定义为抗滑力矩与滑动力矩之比。对于简单土质边坡，一般基于圆弧滑动面假设，忽略土条间竖向相互作用，采用式（1）的Fellenius极限平衡法计算F_s：

式中，i为土条序号，n为土条数，c_i、 $\varphi _i $ 分别为土条i的黏聚力和内摩擦角，l_i、θ_i分别为土条i底边弧长与中点处切线的水平夹角，W_i为土条i的重力。

在边坡可靠性分析中，以F_s < 1表示边坡处于失稳状态，对应事件发生概率 P[F_s < 1]为边坡失效概率 P_f。采用Monte–Carlo法获取N个F_s，统计F_s< 1的数量为M；当N足够大时，由Bernoulli大数定理可知频率M/N收敛于真实P_f，如式（2）所示：

F_s服从正态分布时，按式（3）计算可靠指标β：

式中， $\; {\mu _{{F_{\rm{s}}}}}$ 、 ${\sigma _{{F_{\rm{s}}}}}$ 分别为正态分布变量F_s的均值与标准差， $ {\phi }^{-1} (\cdot) $ 为累计标准正态分布函数的逆函数。

1.2   边坡参数及计算方案

边坡几何参数包括坡高h、坡率1∶m。路基工程在一般条件下的土质单级边坡坡高h≤8 m，方案中取坡高范围4～8 m。坡率变化主要受土质和坡高影响，单级土质边坡条件下一般按1∶1.5设计。

边坡土体强度参数包括黏聚力c、内摩擦角 $\varphi $ 、重度γ等，按各参数相互独立且服从正态分布考虑。路基工程一般为土质边坡，黏聚力均值 $\overline c $ 多为5～25 kPa；内摩擦角均值 $\overline \varphi $ 多取20°～40°；重度均值 $\overline \gamma $ 变化不大，一般取 $\overline \gamma $ =20 kN/m³。土体强度参数的变异性水平，采用变异系数δ表示，按式（4）计算：

式中，μ为土性参数均值，σ为土性参数标准差。

国内外大量工程建设对土性参数变异性做了系统勘查，表1为c、 $\varphi $ 、γ的变异系数δ_c、 $\delta_\varphi $ 、δ_γ的变化范围。

由于δ_γ波动范围较小，分析中取δ_γ=0.05，不考虑δ_γ变化影响。δ_c和 $\delta_\varphi $ 的变化相对较大，并多呈现“同大同小”变化规律^[23]，一般有δ_c/ $\delta_\varphi $ ≈2的关系特征，为此将δ_c和 $\delta_\varphi $ 按一致性变化原则进行简化处理，划分为同变化的小、中、大3个等级，结果见表2。

为提高边坡稳定可靠性分析效率，避免β计算耗时过长，选取坡高4、6、8 m和强度小、中、大变异水平进行组合，制定如表3所示的9种计算方案。

2.   计算数据及分析

2.1   可靠指标变化特性

针对表3计算方案，采用Fellenius法分别求得各方案下441个F_s值；同样地，通过Monte Carlo模拟进行可靠度计算，各组合下的试验次数N=50 000，获得相应的441个β值。图1为β与F_s的关系。

图  1  不同计算方案下 β与F_s的关系

Fig.  1  Relationship between β and F_sunder different calculation scheme

下载: 全尺寸图片

针对土体强度参数在不同变异水平下β ～ F_s的映射特征，采用式（5）幂函数描述可靠度指标的上、下界限值 $ \;{\beta ^{\text{u}}} $ 、 $\; {\beta _{\text{d}}} $ 随F_s变化关系：

式中，β^u为β的上界限值，β_d为β的下界限值，a^u与b^u为对应于β^u的回归系数，a_d与b_d为对应于β_d的回归系数。各回归系数取值见表4。

由表4可知，回归系数与土体强度参数变异水平呈负相关。分析变异系数对回归系数影响时，为避免样本数量过少导致较大的回归偏差，增加了δ_c=0.15、 $\delta_\varphi $ =0.075及δ_c=0.25、 $\delta_\varphi $ =0.125两组变异水平，并获得相应水平下的β界限回归方程系数，如图2和3所示。根据图2和3，可分别采用式（6）的负指数幂函数和式（7）的线性函数描述a^u及a_d、b^u 及b_d与δ_c的关系：

图  2  a^u、a_d与δ_c关系

Fig.  2  Relationship between a^u, a_d and δ_c

下载: 全尺寸图片

图  3  b^u、b_d与δ_c关系

Fig.  3  Relationship between b^u,b_d and δ_c

下载: 全尺寸图片

由式（5）、（6）和（7），根据安全系数、土体强度参数变异系数，可获得β波动范围的界限值，对土体强度参数不同变异水平下边坡稳定性给予一个区间评估，但不能确定边坡稳定的可靠指标β，故需构建一个多因素影响的β函数估算方程，对应为可靠指标估算值 $\; \hat \beta$ ，区别于可靠指标计算值β。

由于β既受土体强度参数均值及变异系数的影响，也与坡高有关，简单分析黏聚力均值 $\overline c $ 、内摩擦角均值 $\overline \varphi $ 对β影响，不能全面反映多因素影响下β的变化规律，且难以消去参数量纲，与无量纲指标β不匹配。为此，运用无量纲化技术处理 $\overline c $ 、 $\overline \varphi $ ，将 $\overline \varphi $ 正切化得摩擦系数并与F_s作商，得到内摩擦角无量纲参数tan $\overline \varphi $ /F_s，各计算方案下可得441个数据；黏聚力无量纲参数的获得，可在tan $\overline \varphi $ /F_s不变的条件下，寻得能较好地描述β变化规律且以 $\overline c $ 为变量的函数表达。分别取坡高h为4、6、8 m，得到对应h下的tan $\overline \varphi $ /F_s∈[0.125,0.406]、[0.163,0.452]、[0.194,0.481]，依据tan $\overline \varphi $ /F_s的共同变化范围为[0.194,0.406]，取tan $\overline \varphi $ /F_s分别为0.22、0.28、0.36时对应上下波动0.02范围内的数据，绘制tan $\overline \varphi $ /F_s变化较小时的β或 $\; \hat \beta $ ～ $\overline c $ 关系曲线，如图4所示；表5为相应的土体强度参数变异水平、坡高h及tan $\overline \varphi $ /F_s这3个因素变化组合情况。

图  4  β或 $\;\hat {\boldsymbol{\beta }}$ ～ $\overline {\boldsymbol{c}}$ 关系曲线

Fig.  4  Curves of β or $\;\hat{\boldsymbol{ \beta}}$ versus $\overline {\boldsymbol{c }}$

下载: 全尺寸图片

由图4和表5可以看出，在tan $\overline \varphi $ /F_s基本相同的条件下，β与 $\overline c $ 呈正相关，可采用对数函数曲线进行描述。为避免 $\overline c $ 接近0时β趋于负无穷，采用ln( $\overline c $ +1)为基函数表达，表5中组合下的拟合曲线R²>0.994，回归效果较好。此外，β与 $\overline c $ 也具有一定的线性关系， $\overline c $ 接近于0时，ln( $\overline c $ +1)/ $\overline c $ 趋于1。因此，运用 $\overline c $ 的对数与线性函数的基函数组合对 $\overline c $ 无量纲化，可得黏聚力无量纲参数ln( $\overline c $ +1)/ $\overline c $ 。

综上，采用无量纲参数x( $\overline c $ )=ln( $\overline c $ +1)/ $\overline c $ 、y( $\overline \varphi $ )=tan $\overline \varphi $ /F_s描述β(x,y)，并依据控制变量法进行影响可靠指标变化的主要因素分析。首先，在y基本不变的条件下，对应表5组合下的x与β近似呈线性关系，可整体表达为拟合值 $\; \hat \beta $ (x)=–Ax+C₁，拟合曲线R²>0.994，如图5所示。另外，在x不变的条件下，分析y=tan $\overline \varphi $ /F_s与β关系，可选取土体强度参数变异水平与坡高的一种典型组合进行讨论，图6为中变异水平下，h=6 m， $\overline c $ 取6、12、18、24 kPa下的β ～ y关系。由图6可以看出， $\overline c $ 取值相同时，β与y呈正相关，采用幂函数得拟合值 $\; \hat \beta $ (y)=By²+C₂，拟合曲线R²值为0.981～0.992。

图  5  β或 $\; \hat {\boldsymbol{\beta}} $ 与x关系曲线

Fig.  5  Curves of β or $\; \hat{\boldsymbol{ \beta}} $ versus x

下载: 全尺寸图片

图  6  β或 $\; \hat {\boldsymbol{\beta }} $ 与 y关系曲线

Fig.  6  Curves of β or $\; \hat{\boldsymbol{ \beta}} $ versus y

下载: 全尺寸图片

根据上述分析，无量纲参数x=ln( $\overline c $ +1)/ $\overline c $ 、y=tan $\overline \varphi $ /F_s在单变量控制条件下，分别具有 $ \;\hat \beta $ (x)=–Ax+C₁线形函数和 $\;\hat \beta$ (y)=By²+C₂幂函数拟合关系。因此，在x、y的综合影响下，可构建组合函数 $\; \hat \beta$ (x,y)=By²–Ax+C对β进行描述；通过数学等式变换，得 $\; \hat \beta$ =C{1+(A/C)[(B/A)y²–x]}，令ζ=C、η₁=A/C，η₂=B/A，得到式（8）可靠指标估算值 $\; \hat \beta$ 与 $\overline c $ 、 $\overline \varphi $ 、F_s 的关系式：

式中，F_s为变量 $\overline c $ 、 $\overline \varphi $ 的函数，系数ζ、η₁、η₂受坡高与土体强度参数变异性影响。

2.2   影响可靠指标的主要因素

在坡高h相同条件下，可靠指标β主要受土体强度参数变异性影响。图7为h=6 m，土体强度参数处于小、中、大变异水平下，β计算值散点及运用式（8）对散点进行拟合的 $\; \hat \beta$ 曲面，拟合优度R²分别为0.975、0.976、0.975，回归方程构建良好。由图7可知，β拟合曲面随强度参数变异水平增加呈顺时针方向旋转，即可靠指标与土体强度参数变异水平具有负相关性。

图  7  强度参数变异性对可靠指标影响

Fig.  7  Influence of the variability of strength parameters on reliability index

下载: 全尺寸图片

在土体强度参数均值及变异系数相同的条件下，可靠指标β主要由坡高h决定。图8为土体强度参数为中变异水平，h取4、6、8 m，由β计算值采用式（8）进行回归的 $ \;\hat \beta $ 曲面，对应的R²值为0.969、0.976、0.965，拟合效果较好。结果表明，β回归曲面随坡高h增大整体向下平移，即可靠指标与坡高呈负相关。

图  8  坡高对可靠指标影响

Fig.  8  Influence of slope height on reliability index

下载: 全尺寸图片

3.   可靠指标估算及偏差分析

3.1   估算方程回归参数

由式（8）组合函数可知， $\; \hat \beta $ 与 $\overline c $ 、tan $\overline \varphi $ 、F_s具有显著的非线性，主要受坡高与土体强度参数变异性综合影响。根据表3中9种边坡计算方案，可得可靠指标估算值 $\; \hat \beta $ 与 $\overline c $ 、tan $\overline \varphi $ 、F_s函数关系式的回归系数ζ、η₁、η₂，见表6。

回归系数ζ主要随土体强度参数变异性增大而减小，定义其为材料效应因数。h取4、6、8 m时的(δ_c,ζ)散点如图9所示；基于Levenberg–Marquardt算法，按式（9）双曲函数描述δ_c与ζ的关系：

图  9  ζ～δ_c关系曲线

Fig.  9  Curve of ζ versus δ_c

下载: 全尺寸图片

η₁、η₂仅受坡高h影响，定义为其几何效应因数。取坡高h为4、5、6、7、8 m，得到图10所示的散点图及拟合曲线。

图  10  η₁、η₂随h变化曲线

Fig.  10  Curves of η₁ and η₂ versus h

下载: 全尺寸图片

由图10可知：η₁与h呈良好的线性关系，回归方程为式（10）；η₂随h增大而非线性减小，采用式（11）负幂函数回归。

3.2   估算偏差分析及区域划分

根据坡高、安全系数、土体强度参数均值及变异系数，由式（8）获得边坡稳定可靠指标估算值 $\; \hat \beta $ 。同时，采用Monte Carlo模拟可获得边坡参数相同条件下的可靠指标计算值β。按表3计算方案，获得的 $\; \hat \beta $ 与β散点位于1∶1均等线附近，如图11所示。这表明 $\; \hat \beta $ 与β具有良好正比关系。

图  11  β与 $ \;\hat {\boldsymbol{\beta}} $ 关系

Fig.  11  Relationship between β and $\; \hat{\boldsymbol{ \beta}} $

下载: 全尺寸图片

采用式（12）表达的相对分析误差RPD（relative percent deviation）进行估算模型预测能力评价，得RPD =11.419 $\gg $ 2，方程的估算效果较好。

式中，N_s为( $\; \hat \beta$ _j, β_j)样本数，j为样本序号， $\overline \beta $ 为β_j平均值。

土体强度参数均值 $\overline c $ 、 $\overline \varphi $ 取值为工程常遇范围时，可靠指标的估算偏差主要受土体强度参数变异系数及坡高影响，见表7。

由表7可知：土体强度参数变异水平位于小～大区间，坡高分别为4、6、8 m，对应的RPD为14.115、12.519、8.983，呈单调递减趋势；坡高在4～8 m范围，土体强度参数分别为小、中、大变异水平，RPD的相应值为6.183、6.295、6.187，波动幅值较小。可靠指标的估算偏差受土体强度参数变异水平影响大于坡高，且随坡高提升呈增大趋势。

图12为坡高h取4～8 m，土体强度参数为小～大变异水平下(F_s, Δβ)散点图，其中，估算偏差Δβ= $\; \hat \beta $ –β值分布在以Δβ=0为中心的–1.25～1.37范围。为此，将Δβ值分布范围划分为4个区域，即：区域I₁，Δβ∈[–0.5,0]；区域I₂，Δβ∈[0,0.5]；区域I₃，Δβ<–0.5；区域I₄，Δβ>0.5。可进一步将I₄划分为I_4s和I_4f，对于I_4s区域，F_s∈[2.051,3.064]，明显大于1.15～1.25的安全系数工程设计值^[24]，相应的β∈[7.096,10.504]，大幅超过了一般边坡工程设计的2.7～3.7目标可靠指标^[25]，即区域I_4s内的边坡稳定可靠指标估算值虽然高于计算值0.5，但边坡仍然具有较高安全性；对于F_s≤1.377的I_4f区域，F_s∈[0.926,1.377]，相应的β∈[–1.627,5.148]，计算工况中的边坡稳定性评价指标F_s与β均有低于工程设计控制值的数据，边坡存在一定失稳风险。进一步统计分析可知，土体强度参数分别为小、中、大变异水平时，Δβ分布于风险区域I_4f范围的百分比分别仅有2.94%、0.91%、0.15%，总体仅占1.33%。

图  12  Δβ与F_s的关系

Fig.  12  Relationship between Δβ and F_s

下载: 全尺寸图片

在土体强度参数处于小～大变异区间，坡高为4～8 m条件下，土体强度参数均值 $\overline c $ 、 $\overline \varphi $ 与估算偏差Δβ关系等值线及区域I₁、I₂、I₃、I_4s、I_4f分布范围如图13所示。由图13可知，区域I₁面积最大，多分布于 $\overline c $ ∈[5 kPa,25 kPa]、 $\overline \varphi $ ∈[25°,35°]范围；I₂面积次之，相应强度参数 $\overline c $ ∈[5 kPa,25 kPa]， $\overline \varphi $ ∈[20°,25°]∪[35°,40°]；估算值显著偏低的I₃面积微小，分布在 $\overline c $ ≈5 kPa， $\overline \varphi $ ∈[27°,40°]狭长范围；估算值显著偏高但稳定安全系数F_s ≥2.051的区域I_4s分布在 $\overline \varphi $ ≈20°、 $\overline c $ ∈[21 kPa,25 kPa]和 $\overline \varphi $ ≈40°、 $\overline c $ ∈[10 kPa,22 kPa]两个窄条形区域；对于估算值显著偏高且F_s ≤1.377的I_4f存在风险区域，强度参数均值分布于 $\overline c $ ∈[5.00 kPa,7.69 kPa]、 $\overline \varphi $ ∈[20.00°,21.77°]较小范围。风险区域I_4f与I₂的分布范围，Δβ=0.50或Δβ=0.25界线近似满足式（13）或（14）的2次抛物线函数：

图  13  Δβ随 $\overline{\boldsymbol{ c}} $ 、 $\overline {\text{φ}} $ 变化等值线图

Fig.  13  Contour map of Δβ versus $\overline {\boldsymbol{c }}$ and $\overline {\text{φ}} $

下载: 全尺寸图片

4.   讨　论

为进一步评价边坡稳定可靠指标估算模型的适用性，选用文献[26-28]中坡率为1∶m=1∶1.5的简单土质边坡算例进行比较分析，边坡参数及稳定可靠指标见表8，参数变异系数取δ_c=0.2、δ_φ=0.1、δ_γ=0.05。依据本文模型获得的估算值 $\;\hat \beta$ 与文献算例可靠指标β之间的误差较小，估算偏差Δβ∈[–0.3,0.1]，相对偏差ε_β∈[–4.89%,5.08%]。

事实上，铁路和公路工程中的路堤还会承受车辆及上部结构荷载，对边坡稳定性有一定影响。进行路堤边坡稳定安全系数计算时，一般将路基面以上的车辆等荷载简化为等效静荷载，其中，公路工程换算成分布于路基面的满布荷载^[29]，铁路工程多按条形荷载布设于路基面^[30-31]。对于铁路路堤边坡，列车和轨道的变异系数δ_q=0.07^[32]，与边坡土体重度的变异系数δ_γ=0.05差异不大，因此，计算考虑路基面荷载作用影响的路堤边坡稳定可靠指标β_q，直接按等效条形荷载处理；而采用本文模型估算可靠指标 $\; \hat \beta $ _q时，可将q换算成满布于路基面的等土体重度及等变异性的土层，通过增加土层换算厚度h₀至原路堤高度h反映。

铁路路堤的单线等效条形荷载q平均值和分布宽度为μ_q=60.625 kN/m²、b=3.6 m，荷载边缘距离坡顶的宽度b'=2.35 m。依据表3计算方案，基于Fellenius法与Monte Carlo模拟，获得考虑荷载作用的路堤边坡可靠指标计算值β_q，同时采用本文模型获得相应估算值 $ \;\hat \beta $ _q。图14为β_q与 $ \;\hat \beta $ _q关系散点，绝大多数都位于1∶1均等线附近，相对分析误差RPD=8.29 $\gg $ 2，表明建立的模型对路堤承受车辆荷载作用的边坡稳定可靠指标估算仍具有良好适用性。

图  14  β_q ～ $\; \hat {\boldsymbol{\beta }}$ _q关系

Fig.  14  Relationship between β_q and $ \;\hat{\boldsymbol{ \beta}} $ _q

下载: 全尺寸图片

5.   结　论

针对常用坡率下简单土质边坡稳定分析模型，依据Fellenius极限平衡法与Monte Carlo模拟，开展了基于稳定安全系数的路基边坡可靠指标快速估算方法研究，主要有以下结论：

1) 边坡稳定可靠指标β随稳定安全系数F_s增加总体呈负幂函数增大趋势，其波动范围亦相应扩大；β的上、下界限值 $\; {\beta ^{\text{u}}}$ 、 $\; {\beta _{\text{d}}}$ 与F_s间存在 $\; {\beta ^{\text{u}}}$ = $ {a^{\text{u}}} $ (1– $ F_{\text{s}}^{{{ - }}{b^{\text{u}}}} $ )、 $\; {\beta _{\text{d}}}$ = $ {a_{\text{d}}} $ (1– ${F_{\text{s}}^{ - {b_{\text{d}}}}}$ )的拟合关系式。回归系数主要由土体强度参数变异性决定， $ {a^{\text{u}}} $ 、 $ {a_{\text{d}}} $ 呈负指数幂函数衰减， $ {b^{\text{u}}} $ 、 $ {b_{\text{d}}} $ 近似线性降低。

2）运用无量纲化和多元非线性回归技术，建立了可靠指标估算值 $\; \hat \beta$ 与安全系数F_s、土体强度参数均值 $\overline c $ 、 $\overline \varphi $ 的组合函数表达式，主要受强度参数变异性影响的材料效应因数ζ随变异系数增大而呈双曲线形减小，而边坡几何效应因数η₁、η₂随坡高增加分别表现出线性增大、负幂函数减小趋势。

3）基于边坡土体和几何参数建立的可靠指标估算模型具有良好预测能力，估算偏差受强度参数变异性影响大于坡高。常遇边坡参数下，可靠指标估算偏差Δβ <0.5的区域占比超过99%；对于Δ β≥0.5且F_s ≤1.377的风险区域，强度参数均值满足不等式 $\overline \varphi $ ≤18.52+1.5 $\overline c $ –0.17 ${\overline c ^2}$ 。