路堤边坡的稳定性是高速公路建设中质量控制的关键环节，目前工程设计中应用较为广泛的方法是极限平衡分析法。用该方法进行稳定性分析分两步：1）对某一形状和位置已知的潜在滑动面，计算其稳定安全系数；2）对许多潜在的滑动面，确定相应最小安全系数的临界滑动面^[1]。现有边坡稳定性分析的方法往往只是上述某一个步骤的解决方法^[1–2]，造成了两个步骤脱离开来的状况：如瑞典法等各种条分法是求解步骤1）的成熟且较精确的方法，但因其需预先确定滑动面、划分土条及公式为多项式求和的特点，使其无法真正融入到步骤2）临界滑动面的搜索过程中，目前往往仅作为安全系数的计算手段；又如遗传算法等各种搜索方法可用于确定临界滑动面，但其本身却又无法获得安全系数，且各搜索方法的搜索性都有显著的差异，都有其自身不可逾越的缺陷^[2]，如最突出的易陷入局部极小值的问题等。反过来，这种脱离也导致了 $ M$ 种条分法和 $ N$ 种搜索方法之间有 $M \times N$ 种组合方式之多，不同的组合所得临界滑动面和最小安全系数也不尽相同，同时这种组合也累加了条分法和搜索方法自身的缺陷。

有学者将两步骤结合起来，以安全系数为目标函数，运用数学理论求其极小值，以期不需搜索直接解得临界滑动面和最小安全系数，但均因各自的局限性无法在实际工程上应用。如Baker^[3]、Revilla^[4]、Ramamurthy^[5]等运用变分法求其极小值，但求解出的结果偏差较大，并不理想；之后De Josselin^[6]，Luceno^[7]等对变分法理论应用的正确性和分析问题的局限性也都提出了质疑。为避免条分法多项式求和的特点造成变分和求偏导不易的情况，罗文强^[8]和杨庚宇^[9–10]等均对条分法做了改进，将安全系数表达成积分的形式，分别采用变分法和多元函数的极值条件求极小值，但其仍是条分法的一种形式，且极值条件需要对至少3个自变量求偏导解超越方程组，计算异常复杂且往往得不出结果；蒋斌松等^[11]利用此法也仅仅给出了纯黏性土边坡（ $\varphi {\rm{ = }}0^\circ $ ）的结果。

为此，针对均质地基路堤，以不采用条分和搜索，将两步骤有机结合起来并给出唯一最小安全系数和临界滑动面为目的，通过对滑动面形状和位置的描述、无条分安全系数公式的推导及运用数学理论直接确定临界滑动面3方面的研究，提出一种可用于工程实际的均质地基路堤稳定性分析的新方法。

1 滑动面形状和位置的描述 1.1 自变量的选择和滑动面函数的确定

路堤填土一般是均质的，当地基也为均质土时，滑动面往往为圆弧形，选取3个独立维度的自变量确定圆弧滑动面函数，方法如下：

如图1所示，路堤边坡坡率为 $k$ ，路堤高度为 $H$ ，路堤底宽为 $L$ ；圆弧 $\overset\frown{ABCD}$ 为一假定滑动面，与坡脚内地平线的交点为 $ C$ ，点 $O'$ 为圆心，滑出点为 $ A$ ，滑入点为 $ D$ ，点 $O$ 为圆心 $O'$ 在地平线上的投影；以点 $O$ 为原点，水平向、竖直向为 $ X$ 、 $ Y$ 轴建立直角坐标系。

选取滑弧与坡脚内地平线的交点 $ C$ 与坡脚 $ F$ 的距离 $l$ 、该交点 $ C$ 与圆心在地平线投影点 $ O$ 的距离 $a$ 、该交点 $ C$ 和投影点 $ O$ 所在半径的夹角 $\alpha $ ，共3个维度的自变量表征滑动面的几何特性，表示为 ${{x}} = {\left( {l,a,\alpha } \right)^{\rm T}}$ 。当其确定时，可唯一地确定一个圆弧滑动面，坡面函数 ${f_1}$ 和滑动面函数 ${f_2}$ 分别为：

${f_1}:y = kx + k\left( {l - a} \right)$

(1)

${f_2}:{x^2} + {\left( {y - a\cot \;\alpha } \right)^2} = \frac{{{a^2}}}{{{{\sin }^2}\alpha }}$

(2)

1.2 自变量各维度取值范围的确定 1.2.1 $l$ 的取值范围

如图1所示， $l$ 为线段 $ FC$ 的长度，坡脚点 $ F$ 为固定的，可见 $l$ 仅与点 $ C$ 的位置有关。根据对称性，滑弧与坡脚内地平线的交点 $ C$ 只可能在坡脚点 $ F$ 和路堤底部中心点 $ G$ 之间变化，故 $l$ 的取值范围为：

$0 < l < \frac{L}{2}$

(3)

1.2.2 $a$ 的取值范围

如图1所示， $a$ 为线段 $ OC$ 的长度，也即线段 $ AC$ 长度的1/2，当 $l$ （点 $ C$ 的位置）确定时， $a$ 仅与滑出点 $ A$ 的位置有关。由于滑出点 $ A$ 不可能在坡脚内侧，故当其与坡脚点 $ F$ 重合，即滑弧过坡脚时，线段 $ AC$ 取最小值 $l$ ， $a$ 也取最小值 ${l / 2}$ ；同时，滑出点 $ A$ 也不可能无限远，即线段 $ AC$ 不可能无限长，取线段 $ AC$ 长度为 $2l$ 时点 $ A$ 的位置为最远的滑出点，也即圆心位于坡脚内侧，这也是潘家铮^[12]、殷宗泽^[13]等大多数学者研究和工程实践的结果，故 $a$ 的取值范围为：

$\frac{l}{2} < a < l$

(4)

1.2.3 $\alpha $ 的取值范围

如图1所示， $\alpha $ 为 $a$ 所对应的圆心角，也即线段 $ AC$ 对应圆心角的1/2；当 $l$ （点 $ C$ 的位置）和 $a$ （点 $ A$ 的位置）确定时， $\alpha $ 仅与圆弧的半径有关，且成反比，也即与滑入点 $ D$ 的位置有关。

1） $\alpha $ 的最小值。

取路堤左幅做示意图2，根据对称性，点 $ D$ 只可能在路肩点 $ E$ 和路顶中心点 $ I$ 之间变化；当滑入点 $ D$ 与路顶中心点 $ I$ 重合时，圆弧半径最大， $\alpha $ 取最小值。

此时 $ D$ 点的坐标为 ${D_{\alpha \min }}\left( L/2 - l + a,H\right)$ ，设此时圆心 $O'$ 的纵坐标为 ${y_{O'{\alpha _{\min }}}}$ ，则有：

${a^2} +y_{O'\alpha _{\min}}^2 = x_{D\alpha _{\min }}^2 + {\left( {{y_{O'{\alpha _{\min }}}} - H} \right)^2},$

解得： $ {y_{O'{\alpha _{\min }}}} = \displaystyle\frac{1}{{2H}}\left( {{{\left( L/2 - l + a \right)}^2} + {H^2} - {a^2}} \right)$ 。

则 $\alpha $ 的最小值为：

${\alpha _{\min }} = \arctan \frac{a}{{{y_{O'{\alpha _{\min }}}}}} = \arctan \frac{{2aH}}{{{{\left( L/2 - l + a \right)}^2} + {H^2} - {a^2}}}$

(5)

2） $\alpha $ 的最大值应根据 $ C$ 点的位置分两段考虑。

如图3所示，当 $ C$ 点位于路肩 $ E$ 点内侧，即 $l > H/k$ 时，随着 $ D$ 点由路堤中心点 $ I$ 逐渐向外侧路肩 $ E$ 点移动，圆弧半径愈小， $\alpha $ 愈大；当路堤顶面 $ EI$ 刚好为圆弧在 $ D$ 点的法线，也即圆心纵坐标为 $H$ 时， $\alpha $ 取得最大值。 $ D$ 点越过该点继续向 $ E$ 点移动，滑动面不再为单一凹函数，而是呈上凸下凹（图3中虚线圆弧），这在现实中也是不可能的^[14–15]，此时 $\alpha $ 最大值为：

${\alpha _{\max }} = \arctan \;\frac{a}{H},l > \frac{H}{k}$

(6)

如图4所示，当 $ C$ 点位于路肩 $ E$ 点外侧，随着 $ D$ 点由路堤中点 $ I$ 逐渐向路肩 $ E$ 点移动， $\alpha $ 单调增加，因此，当 $ D$ 点与路肩 $ E$ 点重合时， $\alpha $ 取得最大值。

此时 $ D$ 点的坐标为 ${D_{\alpha _{\max }}}\left(H/k - l + a,H \right)$ ，设此时圆心 $O'$ 的纵坐标为 ${y_{O'{\alpha _{\max }}}}$ ，则有：

${a^2} + y_{O'\alpha _{\max}}^2 = x_{D\alpha _{\max }}^2 + {\left( {{y_{O'{\alpha _{\max }}}} - H} \right)^2},$

解得： $ {y_{O'{\alpha _{\max }}}} = \displaystyle\frac{1}{{2H}}\left( {{{\left(H/k - l + a \right)}^2} + {H^2} - {a^2}} \right)$ 。

此时 $\alpha $ 的最大值为：

${\alpha _{\max }} \!=\! \arctan \frac{a}{{{y_{O'{\alpha _{\max }}}}}} \!=\! \arctan \frac{{2aH}}{{{{\left( {{H/k} \!-\! l \!+\! a} \right)}^2} \!+\! {H^2} \!-\! {a^2}}},l < \frac{H}{k}$

(7)

2 无条分稳定安全系数公式

安全系数公式的构造即已知滑动面函数 $f\left( {x} \right)$ 求其对应安全系数 $F$ 的过程。以滑动面上部的潜在滑动体为一整体作为研究对象，需要将滑面或分段滑面上的平均正应力和平均切应力设为未知数，并由极限平衡原理建立积分形式的平衡方程（此时积分域为面积分），最终求得整个滑动面上的平均抗剪强度（能发挥的最大平均抗剪应力）与平均剪应力（实际发挥的平均抗剪应力），两者之比即为安全系数。

如图5所示，潜在滑动体受到自身重力、下部土体对其的抗剪力和法向力的作用。

设滑动面上的平均剪应力为 $\tau $ ，其在数值上等于滑动面上实际发挥的平均抗剪应力 $\tau '$ ；一般而言，路堤填土和地基土的材料不同，且滑动面 $\overset\frown{ABC}$ 上的平均正应力比滑动面 $\overset\frown{CD}$ 上大，故将整个滑动面上的平均正应力分两段考虑，即设滑动面 $\overset\frown{CD}$ 和 $\overset\frown{ABC}$ 上的平均正应力分别为 ${\sigma _1}$ 和 ${\sigma _2}$ ，共3个未知数，可由3个平衡方程完全求解。

2.1 力矩平衡

根据力矩平衡条件，潜在滑动体所受的力对圆心 $O'$ 的力矩之和为0，即

$\sum {M = } {M_\sigma } + {M_G}{\rm{ + }}{M_{\tau '}} = 0$

(8)

式中： ${M_\sigma }$ 为法向力对圆心 $O'$ 的力矩，滑动面上每个点的法向力均过圆心，对圆心 $O'$ 的力矩 ${M_\sigma }$ 为0； ${M_G}$ 为潜在滑动体自重对圆心 $O'$ 的力矩，潜在滑动体自重对圆心 $O'$ 的力矩 ${M_G}$ 为 $ AOB$ 、 $ BOC$ 和 $ CDEF$ 这3部分土体的重力对圆心 $O'$ 的力矩之和。由对称性可知，均质地基中 $ AOB$ 和 $ BOC$ 两部分土体的重力对圆心 $O'$ 的力矩大小相等、方向相反，则 ${M_G}$ 在数值上等于土体 $ CDEF$ 的重力对圆心 $O'$ 的力矩。

在土体 $ CDEF$ 中取一微元体，如图5所示。该微元体重力的大小为 ${\gamma _1}{\rm d}x{\rm d}y$ ，方向为竖直向下，对圆心 $O'$ 取矩的力臂长度为微元体的横坐标 $ x$ ，在土体 $ CDEF$ 范围内进行二重积分可得潜在滑动体自重对圆心 $O'$ 的力矩 ${M_G}$ 为：

${M_G} = - \iint\limits_{{S_{CDEF}}} {x{\gamma _1}{\rm d}x{\rm d}y} = - {\gamma _1}\int\limits_0^H {{\rm d}y\int\limits_{{f_1}\left( y \right)}^{{f_2}\left( y \right)} {x{\rm d}x} } $

(9)

式中， ${f_1}\left( y \right) = \displaystyle\frac{y}{k} + a - l,{f_2}\left( y \right) = \sqrt {\displaystyle\frac{{{a^2}}}{{{{\sin }^2}\alpha }} - {{\left( {y - a\cot \;\alpha } \right)}^2}}$ 。

抗剪力对圆心 $O'$ 的力矩 ${M_{\tau '}}$ 为：

${M_{\tau '}} = \left( {2\alpha + \theta } \right){R^2}\tau $

(10)

式中， $R$ 为滑动面的半径， $R = \displaystyle\frac{a}{{\sin\;\alpha }}$ 。

$\theta $ 为半径 $ O'C$ 与 $ O'D$ 的夹角，将 $y = H$ 代入式（2），可得 $ D$ 点的横坐标为：

${X_D} = \sqrt {\frac{{{a^2}}}{{{{\sin }^2}\alpha }} - {{\left( {H - a\cot \;\alpha } \right)}^2}} {\rm{ = }}\sqrt {{a^2} - {H^2} + 2Ha\cot \;\alpha } $

(11)

则

$\theta = \arcsin \frac{{\sqrt {\left( {{a^2} - {H^2}} \right){{\sin }^2}\alpha + Ha\sin\left( 2\alpha \right)} }}{a} - \alpha $

(12)

将式（9）～（12）代入式（8），解得滑动面上的平均剪应力 $\tau $ 为：

$\tau {\rm{ = }}\tau '{\rm{ = }}\frac{{{\gamma _1}{{\sin }^2}\alpha }}{{2{a^2}}}\frac{{ - \displaystyle\frac{{1 + {k^2}}}{{3{k^2}}}{H^3} + \left( {a\cot \;\alpha - \frac{{a - l}}{k}} \right){H^2} + \left( {2a - l} \right)lH}}{{\left( {\alpha + \arcsin \displaystyle\frac{{\sqrt {\left( {{a^2} - {H^2}} \right){{\sin }^2}\alpha + Ha\sin \left(2\alpha \right)} }}{a}} \right)}}$

(13)

2.2 水平方向的静力平衡

根据水平方向的静力平衡条件，潜在滑动体所受的力在 $ X$ 方向的分力之和为0，考虑到重力 $ G$ 和滑动面 $\overset\frown{ABC}$ 上的法向力在 $ X$ 方向的分力均为0，则建立水平方向的静力平衡方程为：

$\int_{ - \alpha }^{\alpha + \theta } {\tau 'R\cos\; \eta {\rm d}\eta } - \int_\alpha ^{\alpha + \theta } {{\sigma _1}R\sin \;\eta {\rm d}\eta } {\rm{ = }}0$

(14)

解得滑动面 $\overset\frown{CD}$ 上的平均正应力 ${\sigma _1}$ 为：

${\sigma _1}{\rm{ = }}\frac{{\sin\; \alpha + \sin \left( {\alpha + \theta } \right)}}{{\cos\; \alpha - \cos \left( {\alpha + \theta } \right)}}\tau '$

(15)

由几何关系对式（15）进行简化得：

${\sigma _1}{\rm{ = }}\frac{{\displaystyle\frac{a}{R} + \displaystyle\frac{{{X_D}}}{R}}}{{\cos \;\alpha - \displaystyle\frac{{R\cos\; \alpha - H}}{R}}}\tau ' = \frac{{a + \sqrt {{a^2} - {H^2} + 2Ha\cot \;\alpha } }}{H}\tau '$

(16)

2.3 竖直方向的静力平衡

根据竖直方向的静力平衡条件，潜在滑动体所受的力在 $ Y$ 方向的分力之和为0，建立竖直方向的静力平衡方程：

$\sum {{F_Y}} = {F_{\tau 'Y}} + {F_{{\sigma _1}Y}} + {F_{{\sigma _2}Y}} + {F_{GY}} = 0$

(17)

式中， ${F_{\tau 'Y}}$ 为抗剪力在 $ Y$ 方向的分力，在滑动面上任取一角度为 ${\rm d}\eta $ 的微圆弧，在整个滑动面范围内进行积分，得：

${F_{\tau 'Y}} = \int_{ - \alpha }^{\alpha + \theta } {\tau 'R\sin\; \eta {\rm d}\eta } = \tau 'R\left( {\cos \;\alpha - \cos \left( {\alpha + \theta } \right)} \right)$

(18)

${F_{{\sigma _1}Y}}$ 为滑动面 $\overset\frown{CD}$ 上法向力在 $ Y$ 方向的分力，在滑动面 $\overset\frown{CD}$ 上任取一角度为 ${\rm d}\eta $ 的微圆弧，在滑动面 $\overset\frown{CD}$ 范围内进行积分，得：

${F_{{\sigma _1}Y}} = \int_\alpha ^{\alpha + \theta } {{\sigma _1}R\cos\; \eta {\rm d}\eta } = {\sigma _1}R\left( {\sin \left( {\alpha + \theta } \right) - \sin \;\alpha } \right)$

(19)

${F_{{\sigma _2}Y}}$ 为滑动面 $\overset\frown{ABC}$ 上法向力在 $ Y$ 方向的分力，在滑动面 $\overset\frown{ABC}$ 上任取一角度为 ${\rm d}\eta $ 的微圆弧，在滑动面 $\overset\frown{ABC}$ 范围内进行积分，得：

${F_{{\sigma _2}Y}} = \int_{ - \alpha }^\alpha {{\sigma _2}R\cos \;\eta {\rm d}\eta } = 2{\sigma _2}R\sin \;\alpha $

(20)

${F_{GY}}$ 为潜在滑动体重力在 $ Y$ 方向的分力，如图5所示，潜在滑动体的重力由 $ ABC$ 和 $ CDEF$ 两部分土体的重力组成，即

${F_{GY}}{\rm{ = }}{\gamma _2}{S_2} + {\gamma _1}{S_1}$

(21)

其中：

$ ABC$ 部分土体的面积 ${S_2}$ 可由扇形 $O'ABC$ 与三角形 $O'AC$ 的面积之差得到：

$ \begin{aligned}[b]& {S_2} = {S_{O'ABC}} - {S_{O'AC}} = \frac{1}{2} \times 2\alpha \times {R^2} - \\& \quad\quad \frac{1}{2} \times 2a \times a\cot \;\alpha = {a^2}\left( {\frac{\alpha }{{{{\sin }^2}\alpha }} - \cot \;\alpha } \right)\end{aligned} $

(22)

在土体 $ CDEF$ 中取一微元体 $ {\rm d}x{\rm d}y$ ，在面积 $ CDEF$ 范围内二重积分，得 $ CDEF$ 部分土体面积 ${S_1}$ 为：

$\begin{aligned} {S\!_1}\!\! =\!\!\!\! & \iint\limits_{{S_{CDEF}}} \!\!\!{{\rm d}x{\rm d}y}\!\! =\!\!\!\!\int\limits_0^H \!\!\!{{\rm d}y\!\!\!\int\limits_{{f_{1\left( y \right)}}}^{{f_{2\left( y \right)}}}\!\!\!{{\rm d}x} } \!\!=\!\!\!\!\int\limits_0^H \!\!{\left(\!\!\! {\sqrt {\frac{{{a^2}}}{{{{\sin }^2}\alpha }}\!\!-\!\! V{{\left( {y\!\! -\!\! a\cot \;\alpha } \right)}^2}}\!\!-\! \frac{y}{k}}\right.} \!- \\ & \Bigg.{a + l} \Bigg)\cdot {\rm d}y=- \frac{1}{{2k}}{H^2}\! +\! \left( {l \!-\! a \!+\! \frac{1}{2}\sqrt {{a^2} \!-\! {H^2}\! +\! 2aH\cot \;\alpha } } \right)H + \\ & \frac{{{a^2}}}{{2{{\sin }^2}\alpha }}\left( {\arctan \left( {\cot \;\alpha } \right) - \arctan \frac{{a\cot \;\alpha - H}}{{\sqrt {{a^2} \!-\! {H^2} \!+\! 2aH\cot \;\alpha } }}} \right)+\\ & \frac{1}{2}a\cot \;\alpha \left( {a - \sqrt {{a^2} - {H^2} + 2aH\cot \;\alpha } } \right)\end{aligned} $

(23)

将式（18）～（23）代入式（17），解得滑动面 $\overset\frown{ABC}$ 上的平均正应力 ${\sigma _2}$ 为：

${\sigma _2}{\rm{ = }}\frac{{{\gamma _1}{S_1} + {\gamma _2}{S_2}}}{{2a}} - \tau '\cot\; \alpha $

(24)

2.4 安全系数公式

设 ${c_1}$ 、 ${\varphi _1}$ 和 ${c_2}$ 、 ${\varphi _2}$ 分别为路堤填土和地基土的黏聚力和内摩擦角，由式（16）和（24），得滑动面 $\overset\frown{CD}$ 和滑动面 $\overset\frown{ABC}$ 上的平均抗剪强度 ${\tau _{\rm f1}}$ 、 ${\tau _{\rm f2}}$ 为：

$\left\{ \begin{aligned}& {\tau _{\rm f1}} = {c_1} + \frac{{a + \sqrt {{a^2} - {H^2} + 2Ha\cot \;\alpha } }}{H}\tau '\tan \;{\varphi _1},\\& {\tau _{\rm f2}}{\rm{ = }}{c_2} + \left( {\frac{{{\gamma _1}{S_1} + {\gamma _2}{S_2}}}{{2a}} - \tau '\cot \;\alpha } \right)\tan\; {\varphi _2}\end{aligned} \right.$

(25)

则整个滑动面上的平均抗剪强度 ${\tau _{\rm f}}$ 为：

${\tau _{\rm f}} = \frac{{\theta R{\tau _{\rm f}}_1 + 2\alpha R{\tau _{\rm f}}_2}}{{\left( {2\alpha + \theta } \right)R}} = \frac{{\theta {\tau _{\rm f}}_1 + 2\alpha {\tau _{\rm f}}_2}}{{2\alpha + \theta }}$

(26)

则安全系数公式 $F$ 为：

$\begin{aligned} F\!\! =\!\! & \frac{{{\tau _{\rm f}}}}{\tau }\!\! =\!\! \frac{\theta }{{\left( {2\alpha \!\!+\!\! \theta } \right)\tau }}\!\left(\!{{c_1} \!\!+\!\! \frac{{\tau '\tan\; {\varphi _1}}}{H} \!\left(\! {a \!\!+\!\! \sqrt {{a^2} - {H^2} + 2Ha\cot \;\alpha } } \right)}\!\!\right)\!+ \\ & \quad\quad \frac{{2\alpha }}{{\left( {2\alpha + \theta } \right)\tau }}\left( {{c_2} + \left( {\frac{{{\gamma _1}{S_1} + {\gamma _2}{S_2}}}{{2a}} - \tau '\cot \;\alpha } \right)\tan \;{\varphi _2}} \right)\end{aligned}$

(27)

式中， ${S_1}$ 、 ${S_2}$ 、 $\theta $ 和 $\tau (\tau ')$ 分别见式（23）、（22）、（12）和（13）。

3 临界滑动面的确定 3.1 基本思路

对目标函数中自变量各维度求偏导，利用多元函数极值条件求得所有的极小值点，再对所有极小值点进行比较确定最小值点，可确保临界滑动面确定过程的全局性。由于对3个维度求偏导解3元超越方程组相当复杂，采取降维度求极值（只对2个维度求偏导）与迭代计算相结合的办法简化确定过程。

基本思路为（图1）：先固定点 $ C$ 的位置不变（即 $l$ 为一定值），根据二元（ $a$ 和 $\alpha $ ）函数的极值条件求得过 $ C$ 点的所有滑动面中安全系数的极小值 ${F_{SC}}$ ， ${F_{SC}}$ 对应的滑动面为过 $ C$ 点的临界滑动面；在坡脚点 $ F$ 和路堤底部中点 $ G$ 的范围内变换点 $ C$ 的位置（在式（3）的范围内变换 $l$ 的大小）进行迭代计算，则每一个 $ C$ 点（l）所对应的 ${F_{SC}}$ 中，最小值即为该边坡的安全系数 ${F_S}$ ，相应的滑动面即为该边坡的临界滑动面。

3.2 降维度求极值

当 $l$ 取 $0 < l < L/2$ 中某一值时，式（27）中仅有 $a$ 和 $\alpha $ 两个未知量，因式（27）为比值的形式，对其求偏导较复杂，在此引入目标函数 $G$ ：

$G\left( {a,\alpha } \right) = {\tau _{\rm f}} - \tau = \frac{{\theta {\tau _{\rm f}}_1 + 2\alpha {\tau _{\rm f}}_2}}{{2\alpha + \theta }} - \tau $

(28)

求 $F$ 的极小值转化为求 $G$ 的极小值，由二元函数的极值条件， $G$ 取得极小值时， $a$ 和 $\alpha $ 必须满足式方程组（29）。

将式（30）～（34）中各变量的偏导代入式（29）即可求解，但式（29）为超越方程组目前尚无法得出解析解，故本文用Newton法求 $a$ 和 $\alpha $ 的数值解，精度取小数点后10位。据经验，均质土边坡一般仅有一个数值解（ $G$ 仅有一个极小值），均质地基和路堤的材料性质不同时有两个数值解，随着地基土层数的增加，数值解的个数也相应增加。但即使为数值解，笔算求解也相当复杂，为此编制了所有变量的代码，并编程了各函数和方程的求解程序，限于篇幅，相关代码和编程部分不做介绍。

在 $\left( {a,\alpha } \right)$ 的各组数值解中，令函数 $G$ 取得最小值的解，也即令函数 $F$ 取得最小值的解，该解所确定的滑动面即为该 $l$ 值下的临界滑动面（过点 $ C$ 的临界滑动面），对应安全系数即为该 $l$ 值下的安全系数 ${F_{SC}}$ 。

$\left\{ \begin{aligned} & \frac{{\partial G}}{{\partial a}} = \frac{{\left( {{\tau _{\rm f}}_1\displaystyle\frac{{\partial \theta }}{{\partial a}} + \theta \displaystyle\frac{{\partial {\tau _{\rm f}}_1}}{{\partial a}} + 2\alpha \displaystyle\frac{{\partial {\tau _{\rm f}}_2}}{{\partial a}}} \right)\left( {2\alpha + \theta } \right) - \left( {\theta {\tau _{\rm f}}_1 + 2\alpha {\tau _{\rm f}}_2} \right)\displaystyle\frac{{\partial \theta }}{{\partial a}}}}{{{{\left( {2\alpha + \theta } \right)}^2}}} - \displaystyle\frac{{\partial \tau }}{{\partial a}} = 0 , \\ & \frac{{\partial G}}{{\partial \alpha }} = \displaystyle\frac{{\left( {{\tau _{\rm f}}_1\displaystyle\frac{{\partial \theta }}{{\partial \alpha }} + \theta \frac{{\partial {\tau _{\rm f}}_1}}{{\partial \alpha }} + 2{\tau _{\rm f}}_2 + 2\alpha \displaystyle\frac{{\partial {\tau _{\rm f}}_2}}{{\partial \alpha }}} \right)\left( {2\alpha + \theta } \right) - \left( {\theta {\tau _{\rm f}}_1 + 2\alpha {\tau _{\rm f}}_2} \right)\left( {2 + \displaystyle\frac{{\partial \theta }}{{\partial \alpha }}} \right)}}{{{{\left( {2\alpha + \theta } \right)}^2}}} - \displaystyle\frac{{\partial \tau }}{{\partial \alpha }} = 0 \end{aligned} \right.$

(29)

$\left\{ \begin{aligned}& \frac{{\partial \tau }}{{\partial a}} = \frac{{\gamma H\sin \;\alpha }}{{6{a^3}{k^2}{{(2\alpha + \theta )}^2}}}\left( \Bigg({\left( {3Hk(a - l) + 3{k^2}l(l - 2a) + {H^2}\left( {{k^2} + 1} \right)} \right)\sin\; \alpha - 3aH{k^2}\cos \;\alpha } \right)a\frac{{\partial \theta }}{{\partial a}}+\Bigg. \\& \quad\quad \Bigg. (2\alpha + \theta )\left( {\left( {3Hk(a - 2l) + 6{k^2}l(l - a) + 2{H^2}\left( {{k^2} + 1} \right)} \right)\sin \;\alpha - 3aH{k^2}\cos \;\alpha } \right)\Bigg) , \\& \frac{{\partial \tau }}{{\partial \alpha }} = \frac{{\gamma H}}{{6{a^2}{k^2}{{(2\alpha + \theta )}^2}}}\Bigg( \left( {\left( {3Hk(a - l) + 3{k^2}l(l - 2a) + {H^2}\left( {{k^2} + 1} \right)} \right)\sin \;\alpha - 3aH{k^2}\cos \;\alpha } \right)\sin \;\alpha \left(\frac{{\partial \theta }}{{\partial \alpha }} + 2\right)-\Bigg.\\& \quad\quad \Bigg. (2\alpha + \theta )\left( {\left( {3Hk(a - l) + 3{k^2}l(l - 2a) + {H^2}\left( {{k^2} + 1} \right)} \right)\sin \left(2\alpha\right) - 3aH{k^2}\cos \left(2\alpha \right)} \right)\Bigg) \end{aligned} \right.$

(30)

$\left\{ \begin{aligned}& \frac{{\partial {\tau _{\rm f1}}}}{{\partial a}} = \frac{{\tan \;{\phi _1}}}{H}\left( {\left( {\sqrt {{a^2} - {H^2} + 2aH\cot \;\alpha } + a} \right)\frac{{\partial \tau }}{{\partial a}} + \tau \left( {\frac{{a + H\cot \;\alpha }}{{\sqrt {{a^2} - {H^2} + 2aH\cot \;\alpha } }} + 1} \right)} \right) , \\& \frac{{\partial {\tau _{\rm f1}}}}{{\partial \alpha }} = \tan\; {\phi _1}\left( {\frac{{\left( {\sqrt {{a^2} - {H^2} + 2aH\cot\; \alpha } + a} \right)}}{H}\frac{{\partial \tau }}{{\partial \alpha }} - \frac{{\tau a{{\csc }^2}\alpha }}{{\sqrt {{a^2} - {H^2} + 2aH\cot\; \alpha } }}} \right) \end{aligned} \right.$

(31)

$\left\{ \begin{aligned}& \frac{{\partial {\tau _{\rm f2}}}}{{\partial a}} = \frac{{\tan\; {\phi _2}}}{{2{a^2}}}\left( {a\gamma \frac{{\partial {S_1}}}{{\partial a}} - 2{a^2}\cot \;\alpha \frac{{\partial \tau }}{{\partial a}} - \gamma \left( {{a^2}\left( {\cot\; \alpha - \alpha {{\csc }^2}\alpha } \right) + {S_1}} \right)} \right) , \\& \frac{{\partial {\tau _{\rm f2}}}}{{\partial \alpha }} = \frac{{\tan \;{\phi _2}}}{{2a}}\left( {\gamma \frac{{\partial {S_1}}}{{\partial \alpha }} - 2a\cot\; \alpha \frac{{\partial \tau }}{{\partial \alpha }} + 2a{{\csc }^2}\alpha (a\gamma - a\gamma \alpha \cot\; \alpha + \tau )} \right) \end{aligned} \right.$

(32)

$\left\{ \begin{aligned} & \frac{{\partial {S_1}}}{{\partial a}} = \left( {a - \sqrt {{a^2} - {H^2} + 2aH\cot \;\alpha } } \right)\cot \;\alpha + \left( {\arctan \frac{{H - a\cot\; \alpha }}{{\sqrt {{a^2} - {H^2} + 2aH\cot\; \alpha } }} + \arctan \left( {\cot \;\alpha } \right)} \right)a{\csc ^2}\alpha - H , \\& \frac{{\partial {S_1}}}{{\partial \alpha }} = - a{\csc ^2}\alpha \left( {\left( {\arctan \frac{{H - a\cot\; \alpha }}{{\sqrt {{a^2} - {H^2} + 2aH\cot\; \alpha } }} + \arctan \left( {\cot \;\alpha } \right)} \right)a\cot \;\alpha + a - \sqrt {{a^2} - {H^2} + 2aH\cot \;\alpha } } \right) \end{aligned} \right.$

(33)

$\left\{ \begin{aligned}& \frac{{\partial \theta }}{{\partial a}} = \frac{{H\sin \;\alpha \left( {H\sin \;\alpha - a\cos\; \alpha } \right)}}{{a\sqrt {{{\left( {a\cos \alpha - H\sin \alpha } \right)}^2}} \sqrt {\left( {{a^2} - {H^2}} \right){{\sin }^2}\alpha + aH\sin \left(2\alpha \right)} }} , \\& \frac{{\partial \theta }}{{\partial \alpha }}\!\! =\!\! - 1 \!\!+\!\! \frac{{aH\cos \left(2\alpha \right)\!\!+\!\! \left( {{a^2}\!\! - \!\!{H^2}} \right)\cos\; \alpha \sin\; \alpha }}{{\sqrt {{{\left( {a\cos \;\alpha \!\!-\!\! H\sin\; \alpha } \right)}^2}} \sqrt {\left( {{a^2} \!\!-\!\! {H^2}} \right){{\sin }^2}\alpha \!\!+\!\! aH\sin \left(2\alpha \right)} }} \end{aligned} \right.$

(34)

3.3 迭代计算

在式（3）的范围内变换 $l$ 的大小，重复第3.2节进行迭代计算，并在迭代计算的过程中逐步增大 $l$ 值的精度。据试算，对于相对简单的路堤边坡，收敛区域一般较为明显，且在进入收敛区域后，当 $l$ 的精度取0.1 m时，安全系数的区分精度已可达到10^–6，效果较理想。故在编制计算程序时，进入收敛区域前后， $l$ 的精度分别取1.0和0.1 m。

由程序求得每一个 $l$ 值（每一个 $ C$ 点）下的安全系数 ${F_{SC}}$ 之后进行对比，其最小值即为该边坡的安全系数 ${F_S}$ ，相应的滑动面即为该边坡的临界滑动面。

4 算例验证 4.1 安全系数公式 $F$ 的正确性验证

算例选取《土力学》^[12]中的例题7-2。已知条件为：土坡高度 $H$ =15 m，坡率 $k$ =0.5，土容重 $\gamma $ =19.5 kN/m³，土黏聚力 $c$ =40 kPa，内摩擦角 $\varphi $ =8°；已知潜在滑动面圆心坐标为（20，–10.08），半径为27 m。例题采用瑞典条分法（土条划分为10个），计算得到对应该滑动面的安全系数为1.43。现采用本文方法求解该算例，将已知潜在滑动面的位置换算到图1中的直角坐标系得 $a$ =10 m， $l$ =20 m， $\alpha = $ 0.379 4，计算过程和结果见表1。

作为对比，用简化Bishop法、Spencer法和Morgenstern-Price法3种条分法分别求解该算例，为统一精度，土条个数均取10个（与例题瑞典法相同），具体计算过程此处略，将不同方法的结果列于表2。

由表2可知，不同方法所得的安全系数很接近，进一步验证了满足平衡条件下，所得的安全系数对滑面正应力分布不甚敏感，可见设平均正应力和平均切应力的方法可取，本文方法和严格条分法的Spencer法、Morgenstern-Price法的精度很接近，偏差0.68%。

4.2 临界滑动面确定的准确性验证

设该算例路堤顶宽34 m（双向6车道），则路堤底宽 $L$ =86 m，由式（3）得 $l$ 的取值范围为 $0 < l < 47\;{\rm{m}}$ 。现将临界滑动面确定程序的部分极小值求解和迭代计算的过程及结果举例列于表3，对每一个 $l$ 的值，由式（4）～（7）得到 $a$ 和 $\alpha $ 的取值范围见表3第2～5列。

由表3可知，对 $l$ 的每个取值，即滑动面与坡脚内地表的每个交点，均有对应该值（该点）下的临界滑动面和安全系数（第6～8列），其中安全系数的收敛区域大约为 $l=24\sim 26 $ m；当 $l$ 、 $a$ 和 $\alpha $ 分别为24.9 m、13.232 6 m和0.452 9时，安全系数为最小值1.427 151，其所确定的滑动面即该路堤边坡的临界滑动面。

作为对比，利用Geo-Slope商业软件，分别基于瑞典法、简化Bishop法、Spencer法和M-P法与传统搜索方法（本文采用穷举法）的组合方式，计算各潜在滑动面（实际运算数量为2 160个）的安全系数，并取最小值为该算例路堤边坡的最小安全系数，对应滑动面为临界滑动面。建模过程中的几点关键处理有：考虑对称性取路堤的一半建立模型，滑入点范围为路堤顶面，滑出点范围为坡脚至坡脚外23 m，各条分法计算安全系数时土条个数均取30个以提高精度，计算模型见图6。

以模型左下角为坐标原点，将各方法的路堤边坡最小安全系数、临界滑动面的圆心坐标和半径列于表4，并将各法的临界滑动面位置绘于图7。

由表4和图7可知：瑞典法的潜在滑动体范围更大，滑弧入土更深，滑入点和滑出点的水平距离更大，其安全系数也最小，且较其他各法而言小的幅度较大。其根本原因是瑞典法完全忽略了条间力对滑面底部正应力的贡献，正应力最小，抗剪强度也最小，工程上更偏于保守；本文方法与简化Bishop法、Spencer法、M-P法解得的临界滑动面位置（潜在滑动体范围）相差不大，安全系数也很接近，如比M-P法大0.87%。

采用平均正应力这种形式，其实质是为了在整个潜在滑动体积分域内实现整体极限平衡分析；且在严格满足极限平衡条件的情况下，安全系数对正应力的分布并不敏感，只要精度满足工程需要，对正应力进行合理的处理或假设是可行的。由算例可见本文方法精度较高，且其最大的优点在于经自变量的选择和正应力的处理，可直接运用数学理论确定临界滑动面，即确定安全系数的最小值。

5 结　论

提出一种均质地基路堤无条分稳定安全系数计算及临界滑动面确定的新方法，并编制了计算程序，该法不需搜索即可直接确定临界滑动面的位置，并可给出唯一的最小安全系数；同时该法也适用于其他类似的填挖方工程，如均质的土坡和基坑边坡、均质地基上修建的土坝、防洪/防波堤等。

1）以滑弧与坡脚内地平线的交点与坡脚的距离、该交点与圆心在地平线投影点的距离、该交点和圆心投影点所在半径的夹角，共3个维度的自变量可唯一确定圆弧形滑动面的位置，其取值范围不需经验假定，可实现临界滑动面确定范围的全局性。

2）严格满足极限平衡条件的情况下，安全系数对滑面底部正应力的分布并不敏感，对正应力进行合理的处理或假设，以潜在滑动体为研究对象进行整体极限平衡分析，可得出合理的安全系数值。

3）基于二元函数的极值条件，对积分形式的平衡方程解出的安全系数函数表达式，通过降维度求极小值与迭代计算结合，可直接确定临界滑动面位置，克服了传统搜索方法易陷入局部极小值的缺点。

4）算例结果表明：同一滑动面的安全系数，本文方法比M-P法小0.68%；本文方法确定的临界滑动面位置与严格满足极限平衡条件的条分法+传统搜索方法所确定的临界滑动面位置偏差较小；临界滑动面对应的路堤边坡最小安全系数比M-P条分法+传统搜索方法所得结果大0.87%。