长期暴露于外部环境中的边坡浅层土体，受气候、降雨、风力等作用的交替变化，常年处于干湿循环状态，密实度逐年降低，与深层土体存在明显差异。降雨入渗导致土体容重增加，抗剪强度衰减，在渗流耦合作用下容易诱发边坡浅层溜滑失稳破坏^[1]。土质浅层溜滑广泛分布于中国东南和西南部的天然残积土边坡及粉质黏土、粉土和黏土路基边坡中，失稳厚度集中在0.5～1.5 m^[2–3]，具有规模小、频率高的特点，90%以上发生在雨季汛期^[3]。

降雨入渗会在土质边坡内部形成明显的湿润锋，湿润锋下移过程中形成以其为分界面的强度差异土层^[4]，对分析土质边坡浅层失稳模式至关重要。Skempton和Delory^[5]基于无限长斜坡的假定，提出降雨饱和形成坡面径流条件下的“顺坡平面”破坏模式（即滑动体长度L远大于滑动深度z_w的近似无限长斜坡体），以及相应的稳定系数计算公式。马世国等^[6]沿用无限长“顺坡平面”破坏模式，分析了强降雨和地下水对浅层土质边坡稳定性的综合影响。张杰等^[7]基于降雨入渗分层的假定对无限长黄土边坡的稳定性进行分析。韩同春^[8]和Cho^[9]等分析了降雨入渗条件下双层结构无限长土质边坡的稳定性。Dou等^[10]分析了边坡内部渗透特性的空间变化对无限长边坡稳定性的影响。Griffiths^[11]和Milledge^[12]等均基于数值计算，分别讨论了无限长斜坡稳定性分析方法的适用性，指出滑动范围L与入渗深度z_w比值较小时该方法过于保守且受土体性质的影响较大，认为忽略边坡上、下边界效应会严重影响土质边坡浅层失稳判断的准确性。黄尚燕^[13]将土质边坡浅层溜滑简化为“顺坡折线”破坏模式，考虑了滑坡体下端的抗力作用，但假设条件过多，不方便应用，且没有考虑渗流作用。陈善雄等^[14]采用滑坡体中段为顺坡平面、上下为随缘机搜索的3段折线型失稳模式。Liu、Oh等^[15–16]采用了常规的整体圆弧面滑动模式。连继峰等^[4,17]认为边坡浅层溜滑为组合式“顺坡曲面”失稳模式，基于Mohr–Coulomb强度理论分析了斜坡中点的应力状态，推导出反映主滑坡体长度和入渗深度的“顺坡曲面”组合滑动的边坡浅层稳定性分析方法。然而，对于不同的破坏失稳模式，对应不同的滑动面建立相适应的稳定分析方法，才能得到与实际情况相符的结果。

综上所述，合理判断降雨入渗软化下土质边坡的失稳模式是稳定性分析的前提，滑动体范围与土体力学特性、边坡几何参数之间的关联尚不明确。天然土质边坡或路基边坡常处于非饱和状态，因此，作者基于上述土质边坡浅层失稳模式及非饱和土的强度理论^[18–21]，分析非饱和土质边坡任意深度处的土体达到极限应力状态所需要满足的先行条件，并建立相应的稳定性分析方法，进一步分析潜在滑动范围与滑坡深度之间的关系，并考虑土体物理力学性质和边坡几何参数的影响。研究成果将对深入认识土质边坡的失稳模式具有参考价值。

1 非饱和土理论 1.1 土–水特征曲线

非饱和土的孔隙中同时包含水和空气，水–气分界面具有表面张力，孔隙气压力u_a和孔隙水压力u_w并不相等，且有u_a>u_w，该压力差为基质吸力，记为 $\psi$ =u_a–u_w（单位 kPa）。在不排气、不排水的情况下，u_a和u_w均随外部压力的增加而增大。非饱和土中的基质吸力 $\psi$ 与体积含水率θ_w有关，其关系曲线为土–水特征曲线，需通过试验测得。Genuchten^[22]提出一个半经验公式用于描述土–水特征曲线：

$S\!e = \frac{{\theta - {\theta _{\rm{r}}}}}{{{\theta _{\rm{s}}} - {\theta _{\rm{r}}}}} = {\left[ {\frac{1}{{1 + {{\left( {\alpha \psi } \right)}^b}}}} \right]^a}$

(1)

式中：Se为饱和度； $\psi$ 为基质吸力；θ_r为残余体积含水率；θ_s为饱和体积含水率；a、b和 $ \alpha$ 分别为基于试验结果得到的回归系数，a=1–1/b，b>1。

1.2 非饱和土强度理论

由于基质吸力的影响，非饱和土的强度准则与饱和土有所差异。降雨后土体的含水量增加，基质吸力减小、强度降低，基质吸力对强度的影响是非线性的。其中，Fredlund等^[19]基于Bishop有效应力概念提出的非饱和土双线性强度准则被广泛应用：

${\tau _{\rm{f}}} = c' + \left( {{\sigma _{\rm{n}}} - {u_{\rm{a}}}} \right)\tan\; \varphi ' + \left( {{u_{\rm{a}}} - {u_{\rm{w}}}} \right)\tan\; {\varphi _{\rm{b}}}$

(2)

式中， ${\tau _{\rm{f}}}$ 为抗剪强度，c′和 $\varphi '$ 分别为土体的有效黏聚力和内摩擦角， ${\sigma _{\rm{n}}}$ 为土体自重及外部荷载产生的总应力， $\varphi_{\rm b}$ 为抗剪强度随基质吸力变化的摩擦角。

Lu等^[20–21]基于Bishop有效应力概念，提出同时适用于饱和土和非饱和土的有效应力 $ \sigma '$ 的计算公式：

$\sigma ' = \left( {{\sigma _{\rm{n}}} - {u_{\rm{a}}}} \right) - {\sigma _{\rm{s}}}$

(3)

式中， $\sigma_{\rm s}$ 为用归一化体积含水率或饱和度表示的吸应力（Lu和Godt^[20]），可表示为：

${\sigma _{\rm{s}}} = \left\{ {\begin{aligned} & { - S\!e\left( {{u_{\rm{a}}} - {u_{\rm{w}}}} \right) = - S\!e \cdot \psi },\;{{u_{\rm{w}}} < 0}{\text{；}}\\ & {{u_{\rm{w}}}},\;{{u_{\rm{w}}} \ge 0} \end{aligned}} \right.$

(4)

作者针对土质边坡的浅层土体，暂不考虑气压的影响，可认为土体中的孔隙与外界空气相互连通，土体中的孔隙气压力 ${u_{\rm{a}}}{\rm{ = }}0$ 。对于非饱和土体， ${\sigma _{\rm{s}}} = S\!e \cdot {u_{\rm{w}}} < 0$ （ ${u_{\rm{w}}} < 0$ ）；对于饱和土体， ${\sigma _{\rm{s}}} = {u_{\rm{w}}} \ge 0$ （ ${u_{\rm{w}}} \ge 0$ ）。同时适用于饱和土和非饱和土抗剪强度的统一表达式为：

${\tau _{\rm{f}}} = c' + \left[ {\left( {{\sigma _{\rm{n}}} - {u_{\rm{a}}}} \right) - {\sigma _{\rm{s}}}} \right]\tan\; \varphi '$

(5)

2 土质边坡浅层破坏模式 2.1 斜坡土体中点的应力状态

实际土质边坡的失稳破坏形态是由上缘的弧状张拉破坏区、中间的平移主滑动破坏区和下缘的弧状挤压破坏区3部分组成^[17]，如图1所示，将该破坏模式称为“上、下缘顺坡曲面”组合破坏。斜坡土体自重产生的竖向应力为 ${p_{\rm{v}}} = \gamma {{\textit{z}}_{\rm{w}}}\cos\; \beta $ ；竖向有效应力为 ${p'_{\rm{v}}} = \left( {\gamma {{\textit{z}}_{\rm{w}}}\cos\; \beta - {u_{\rm{a}}}} \right) - {\sigma _{\rm{s}}}$ ，为该面上有效正应力 ${\sigma '_{\text{z}}}$ 和剪切应力 ${\tau _{\text{z}}}$ 的合力。中间顺坡滑动平面上B点（图1）的应力状态与Mohr圆的关系如图2所示。

对于无限长斜坡浅层土体沿顺坡平面整体滑动的情况，根据抗滑力与下滑力相等（稳定安全系数FS=1.0），可得到临界深度 ${{\textit{z}}_{{\rm{cr}}}}$ 的计算公式为：

${{\textit{z}}_{{\rm{cr}}}} = \frac{{c' - \left( {{u_{\rm{a}}} + {\sigma _{\rm{s}}}} \right)\tan\; \varphi '}}{{\gamma \cos\; \beta \sin\; \beta - \gamma {{\cos }^2}\;\beta\; \tan\;\varphi '}}$

(6)

当 ${{\textit{z}}_{\rm{w}}} > {{\textit{z}}_{{\rm{cr}}}}$ 时，任意长度为L₂的土条将产生剩余下滑力T，这也是土质边坡发生浅层失稳的必要条件，计算公式为：

$\begin{aligned}[b] T =& \gamma {{\textit{z}}_{\rm{w}}}\cos\; \beta \sin\; \beta {L_2} - c'{L_2} -\\ & \left[ {\left( {\gamma {{\textit{z}}_{\rm{w}}}{{\cos }^2}\;\beta - {u_{\rm{a}}}} \right) - {\sigma _{\rm{s}}}} \right]\tan \;\varphi '{L_2} \end{aligned}$

(7)

式中， $ \gamma $ 表示土体的重度， $\beta$ 为边坡坡度。

如图2所示，A_v、B_v和C_v的横、纵坐标分别代表斜坡内A、B和C点处的有效正应力和剪切应力。根据Mohr–Coulomb强度准则，B_v点在强度包线上，说明斜坡上B点处于临界破坏状态，对应着临界深度z_cr；A_v和C_v点位于强度包线的下方，说明斜坡上A和C点处于安全状态，对应的深度z_w均小于z_cr。以上表明，土体的有效内摩擦角 $\varphi '$ 小于斜坡坡度 $\beta$ 时（图2），斜坡必然存在浅层坡体处于极限状态时对应的临界深度。然而，实际情况是斜坡并非无限长，顺坡面BB′上的土体处于极限状态，滑动体的上、下缘均在BB′面的上方，处于安全状态，对斜坡的稳定性起着关键作用。因此，有必要对滑动体上、下缘土体的极限状态进行讨论分析。

2.2 滑动体上、下缘土体的应力极限状态

如图1所示，中间主滑动体在沿顺坡面向下滑动时，上、下缘土体将分别受到张拉和挤压作用。也即，中间滑动体将分别承受上、下缘土体的主、被动土压力。在临界深度z_cr的顺坡面BB′上方，滑动土体上、下缘任意深度A点（z_a<z_cr）的应力状态有3种可能：过A_v点应力圆与竖向应力线OB_v右侧相割、相切和左侧相割，且均不与强度包线相切，分别如图3中应力圆C₁、C₂和C₃所示。

图3中，相割点D₁、相切点A₁和相割点E₁处的应力 ${p'_{11}}$ 、 ${p'_{12}}$ 和 ${p'_{13}}$ 分别为A点菱形单元顺坡向上的应力。对于右侧应力圆C₁， ${p'_{\rm{v}}} < {p'_{11}}$ ；相切应力圆C₂， ${p'_{\rm{v}}} = {p'_{12}}$ ；左侧应力圆C₃， ${p'_{\rm{v}}} > {p'_{13}}$ 。若斜坡内任意点A的竖向应力 ${p'_{\rm{v}}} \le {p'_1}$ ，将存在某一点C_v（图2）在顺坡平面上B点达到临界状态之前已经处于极限破坏状态，且与顺坡平面成0.5η的夹角。因此，只有 ${p'_{\rm{v}}} < {p'_1}$ 的左侧相割应力圆C₃的情况符合，边坡浅层溜滑基本呈顺坡平滑破坏模式。

滑动体上缘A′点要达到极限平衡状态，必须同时满足式（8）和（9）；滑动体下缘A点要达到极限平衡状态，必须同时满足式（8）和（10）。式（8）为Mohr应力圆上A点的应力状态方程，式（9）和（10）为Mohr–Coulomb强度理论的极限状态方程：

${\left[ {\frac{{{{\sigma _1'}} + {{\sigma_3'}}}}{2} - {{\sigma}_{\text{z}}'}} \right]^2}{\rm{ + }}\tau _{\text{z}}^2 = {\left( {\frac{{{{\sigma_1'}} - {{\sigma_3'}}}}{2}} \right)^2}$

(8)

${\sigma '_3} = {\sigma '_1}{K_{\rm{a}}} - 2c'\sqrt {{K_{\rm{a}}}}\quad$

(9)

${\sigma '_1} = {\sigma '_3}{K_{\rm{p}}} + 2c'\sqrt {{K_{\rm{p}}}}\!\!\!\quad\quad $

(10)

式中： ${\sigma '_1}$ 和 ${\sigma '_3}$ 分别为有效大、小主力， ${\sigma '_1} = \left( {{\sigma _{{\rm{n}}1}} - {u_{\rm{a}}}} \right) -$ $ {\sigma _{\rm{s}}}$ ， ${\sigma '_3} = \left( {{\sigma _{{\rm{n3}}}} - {u_{\rm{a}}}} \right) - {\sigma _{\rm{s}}}$ ， ${\sigma _{\rm n1}}$ 和 ${\sigma _{{\rm{n}}3}}$ 为总的大、小主应力；K_a和K_p分别为主动和被动土压力系数， ${K_{\rm{a}}} = {\tan ^2}( {{45}^ \circ } - $ $0.5\varphi ')$ ， $ {K_{\rm{p}}} = {\tan ^2}\left( {{{45}^ \circ } + 0.5\varphi '} \right)$ ； ${\sigma '_{\text{z}}}$ 和 ${\tau _{\text{z}}}$ 分别为斜坡内A点土体的有效法向应力和切向应力：

${\sigma '_{\text{z}}} = \left( {\gamma {{\textit{z}}_{\rm{w}}}{{\cos }^2}\;\beta - {u_{\rm{a}}}} \right) - {\sigma _{\rm{s}}}$

(11)

${\tau _{\text{z}}} = \gamma {{\textit{z}}_{\rm{w}}}\cos\; \beta \sin\; \beta $

(12)

联立式（8）～（12），可以得到斜坡内A点土体达到极限状态的判别方程如下：

$ \begin{aligned}[b] & {\left[ {\left( {\gamma {{\textit{z}}_{\rm{w}}}{{\cos }^2}\;\beta \!-\! {u_{\rm{a}}}} \right) \!-\! {\sigma _{\rm{s}}}} \right]^2} \!=\! \left[ {\left( {\gamma {{\textit{z}}_{\rm{w}}}{{\cos }^2}\;\beta \!- \!{u_{\rm{a}}}} \right) \!-\! {\sigma _{\rm{s}}}} \right] \times\\ & \quad\left[\!{{{\sigma_1'}}\left(\! {1 \!+\! {K_{\rm{a}}}} \!\right) \!-\! 2c'\!\sqrt {{K_{\rm{a}}}} }\right] - \left[\! {{K_{\rm{a}}}{{{{{\sigma}_1'}}}^2} \!-\! 2c'{{\sigma}_1'}\!\sqrt {{K_{\rm{a}}}} } \right]\!-\\ & \quad\quad\quad\quad\quad\quad {{{\left(\! {\gamma {{\textit{z}}_{\rm{w}}}\sin \;\beta \cos\; \beta } \!\right)}^2}} \end{aligned}$

(13)

$ \begin{aligned}[b] & {\left[{\left( {\gamma {{\textit{z}}_{\rm{w}}}{{\cos }^2}\;\beta \!-\! {u_{\rm{a}}}} \right) \!-\! {\sigma _{\rm{s}}}} \!\right]^2} \!=\!\left[ {\left({\gamma {{\textit{z}}_{\rm{w}}}{{\cos }^2}\;\beta \!-\! {u_{\rm{a}}}} \!\right) \!- \!{\sigma _{\rm{s}}}}\! \right] \times \\ & \quad\left[ \!{{{\sigma_1'}}\left(\! {1 \!+\! {K_{\rm{p}}}} \right) \!+\! 2c'\!\sqrt {{K_{\rm{p}}}} } \right] \!+\! \left[\! {{K_{\rm{p}}}{{{{{\sigma}_1'}}}^2} \!-\! 2c'{{\sigma}_1'}\!\sqrt {{K_{\rm{p}}}} } \right]\! -\\ &\quad\quad\quad\quad\quad\quad {{{\left(\! {\gamma {{\textit{z}}_{\rm{w}}}\sin\; \beta \cos\; \beta }\! \right)}^2}} \end{aligned}$

(14)

根据土体的土–水特征曲线，由式（13）和（14）可求出不同饱和状态下斜坡土体A点的有效大、小主应力 ${\sigma '_1}$ 和 ${\sigma '_3}$ 。故可得到过A点的2个与强度包线相切的Mohr圆，如图4中的应力圆C₄和C₅所示。其中，G₁和H₁分别为线段OB₁与应力圆C₄和C₅的相割点。

当顺坡平面深度z_w>z_cr时（图1），任意长度L₂的土条将产生剩余下滑力T，在该下滑力T的作用下E₁点处的应力p₁将受到张拉，减小至最小G₁点，滑动体上缘发生张拉剪切破坏；在T的作用下E₁点处的应力p₁受到挤压增加到最大H₁点，滑动体下缘发生挤压剪切破坏。图4中，圆心角η₃和η₁的一半分别为滑动体上、下缘A点土体达到极限状态时破裂面与顺坡平面的夹角。

建立坐标系 $\sigma 'O\tau $ ，如图4所示。已知A_v点坐标为( ${\sigma '_{\text{z}}}$ ， ${\tau _{\text{z}}}$ )（式（11）、（12）），设M和N点的坐标分别为( ${\sigma '_{\rm{m}}}$ ， ${\tau _{\rm{m}}}$ )和( ${\sigma '_{\rm{n}}}$ ， ${\tau _{\rm{n}}}$ )，O₁和O₃的坐标分别为( ${\sigma '_{O1}}$ ，0)和( ${\sigma '_{O3}}$ ，0)。由线段A_vO₁和MO₁相等，以及线段A_vO₃和MO₃相等，可得：

$\sqrt {{{\left( {{{\sigma_{\text{z}}'}} - {{\sigma}_{O1,3}'}} \right)}^2} + \tau _{\text{z}}^2} = \left( {{{\sigma}_{O1,3}'} + \frac{{c'}}{{\tan\; \varphi '}}} \right) \cdot \sin\; \varphi '$

(15)

$\begin{array}{l} {{\sigma_{O1,3}'}} = \dfrac{1}{{{{\cos }^2}\varphi '}}\bigg[ {\left( {{{\sigma}_{\text{z}}'} + c'\cos \;\varphi '\sin \;\varphi '} \right) \pm } \bigg.\\ \left. {\sqrt {{{\left( {{{\sigma}_{\text{z}}'} \!+\! c'\cos\; \varphi '\sin\; \varphi '} \right)}^2} \!-\! {{\cos }^2}\varphi '\left[ {{{ {{{\sigma}_{\text{z}}'}} }^2} \!+\! \tau _{\text{z}}^2 \!-\! {{\left( {c'\cos\; \varphi '} \right)}^2}} \right]} } \right] \end{array}$

(16)

${\tau _{{\rm{n}},{\rm{m}}}} = \sqrt {{{\left( {{{\sigma}_{\text{z}}'} - {{\sigma}_{O1,3}'}} \right)}^2} + \tau _{\text{z}}^2} \cdot \cos \;\varphi '$

(17)

${\sigma '_{\rm{n,m}}} = {\sigma_{O1,3}'} - {\tau _{\rm{n,m}}} \cdot \tan\; \varphi '\!\!\!\!\!\quad$

(18)

因而，可得：

$\cos\; {\eta _{1,3}} = 1 - \frac{{{{\left( {{{\sigma}_{\text{z}}'} - {{\sigma}_{\rm{n,m}}'}} \right)}^2} + {{\left( {{\tau _{\rm z}} - {\tau _{\rm{n,m}}}} \right)}^2}}}{{2{{\left[ {{{\left( {{{\sigma}_{\text{z}}'} - {{\sigma}_{O1,3}'}} \right)}^2} + \tau _{\rm{z}}^2} \right]}}}}$

(19)

式中：η₁的取值范围为0到0.5π– $\varphi'$ ；η₃的取值范围为0到0.5π+ $\varphi'$ ，若考虑滑动体上缘张拉裂缝的影响，由主动土压力计算公式可得上缘张拉裂缝的深度 ${{\textit{z}}_0} =$ $ {{2c'} / {\gamma \sqrt {{K_{\rm{a}}}} }}$ ，此时η₃的最大值为0.5π+ $ \varphi'$ – $\beta $ 。如图1所示的“上、下缘顺坡曲面”组合破坏，林鸿洲等^[23]通过室内模型试验发现下缘挤压区位于坡脚处，张拉区裂纹可沿坡面逐渐发展到坡顶，形成全坡面失稳。

2.3 滑动体上、下缘破裂面形态

对滑动体上、下缘区剪切破坏曲面沿斜坡面分别建立直角坐标系(图5)，对应的表征函数分别为y=g(x)和y=f(x)，曲线上任意一点的斜率分别为g′(x)=tan(0.5η_3i)和f ′(x)=tan(0.5η_1i)。采用定积分元素法，任取一点ξ_i(x_i≤ξ_i≤x_i-1)，导数f ′(ξ_i)可取[f ′(x_i)+f ′(x_i-1)]/2代替区间[x_i，x_i-1]的斜率，写成差分形式为：

$\left\{ \begin{aligned} & \frac{{g\left( {{x_i}} \right) - g\left( {{x_{i - 1}}} \right)}}{{\Delta x}} \approx \frac{1}{2}\left( {\frac{{\partial g\left( x \right)}}{{\partial {x_i}}} + \frac{{\partial g\left( x \right)}}{{\partial {x_{i{\rm{ - }}1}}}}} \right){\text{，}}\\ & \frac{{f\left( {{x_i}} \right) - f\left( {{x_{i - 1}}} \right)}}{{\Delta x}} \approx \frac{1}{2}\left( {\frac{{\partial f\left( x \right)}}{{\partial {x_i}}} + \frac{{\partial f\left( x \right)}}{{\partial {x_{i{\rm{ - }}1}}}}} \right) \end{aligned} \right.$

(20)

有初始条件：

${\left. {g\left( x \right)} \right|_{x = {x_n}}} = 0, \;\; {\left. {f\left( x \right)} \right|_{x = {x_n}}} = 0$

(21)

${\left. {\frac{{\partial g\left( x \right)}}{{\partial {x_i}}}} \right|_{x = {x_i}}} = \tan\; {\frac{{{\eta _{3i}}}}{2}}, \;\; {\left. {\frac{{\partial f\left( x \right)}}{{\partial {x_i}}}} \right|_{x = {x_i}}} = \tan\;{\frac{{{\eta _{1i}}}}{2}}$

(22)

联立式（20）～（22），可得上缘张拉区的滑动范围L₃及下缘挤压区滑动范围L₁的表达式，分别为：

$\begin{array}{l} {y_{3n}} = g\left( {{x_n}} \right) = \Delta x\left({\dfrac{1}{2}\tan \;\dfrac{{{\eta _{3n}}}}{2} + \tan \;\dfrac{{{\eta _{3n - 1}}}}{2}} \right.+\\ \quad\quad\tan \;\dfrac{{{\eta _{3n - 2}}}}{2} + \cdots + \tan \;\dfrac{{{\eta _{3i}}}}{2}\left. { + \cdots + \tan \;\dfrac{{{\eta _{31}}}}{2} + \dfrac{1}{2}\tan \;\dfrac{{{\eta _{30}}}}{2}} \right) \end{array}$

(23)

$\begin{array}{l} {y_{1n}} = f\left( {{x_n}} \right) =\Delta x\left({\dfrac{1}{2}\tan\; {\dfrac{{{\eta _{1n}}}}{2}} + \tan\; {\dfrac{{{\eta _{1n - 1}}}}{2}} } \right. +\\\quad\quad\tan \; {\dfrac{{{\eta _{1n - 2}}}}{2}} + \cdots + \tan\; {\dfrac{{{\eta _{1i}}}}{2}}\left. { + \cdots + \tan\; {\dfrac{{{\eta _{11}}}}{2}} + \dfrac{1}{2}\tan \; {\dfrac{{{\eta _{10}}}}{2}} } \right) \end{array}$

(24)

边坡滑坡体坡长 $L = L{}_3 + {L_2} + {L_1}$ ，上、下缘滑坡体的长度 $L{}_3$ 和 $L{}_1$ 分别为：

${L_3} = {y_{3n}} - {{\textit{z}}_{{\rm{cr}}}} \cdot \sin\; \beta $

(25)

${L_1} = {y_{1n}} + {{\textit{z}}_{\rm cr}} \cdot \sin \;\beta $

(26)

为得到滑坡体上、下缘的破裂面形态，相关土体的物理力学指标、斜坡坡比（1∶1.5）、土–水特征曲线相关参数的取值如表1所示。

表1中，土水特征曲线参数是由式（1）对文献[24]试验数据的拟合所得。图6给出了z_w=z_cr处，由上述差分法所得不同饱和度Se、土体强度参数、边坡坡度下滑动体上、下缘的破裂面形态。对某一因素的影响进行分析时，其他参数均保持恒定。

由图6可知：滑动体上、下缘的范围与土体的物理力学指标直接相关，上、下缘滑裂面的形态仅受内摩擦角 $\varphi'$ 的影响，且随 $\varphi'$ 的增加而增大，并逐渐趋于平缓；随着边坡坡度 $\beta $ 的增加，滑动范围逐渐减小，同样逐渐趋于平缓；滑动面的范围随黏聚力c′增加逐渐增大，随饱和度Se的增加呈逐渐减小趋势，但其形态并不受c′ 和Se变化的影响；对于饱和土体，斜坡滑动区域受c′变化的影响远大于 $ \varphi'$ 和Se。以上分析说明，滑动体上、下缘危险滑动面的位置随强度参数的减小、饱和度的增加及坡度的增加逐渐向浅层迁移。

2.4 滑动体上、下缘破裂面的数学表达

虽然根据式（23）和（24）可分别求出滑动体上、下缘破裂面的形态，但计算过程比较复杂，不便于边坡浅层稳定性分析的应用。对数螺旋线形式的滑动面在地基承载力计算、边坡稳定性及支挡结构土压力的分析等方面有着广泛的应用^[25]，作者也将采用对数螺旋线函数对滑动体上、下缘破裂面的形态进行描述：

$R = {R_0}\exp \left( \delta \tau _{\rm{n,m}} \cdot \tan\; \varphi '\right)\!\!\!\!\!\quad$

(27)

式中，R为对数螺旋线上任一点的半径，R₀为初始半径，δ为R与R₀之间的夹角。根据土体强度理论，2组滑移线在任一点处的交角为（90°± $\varphi'$ ）。基于此，图7给出了边坡浅层失稳滑动体上、下缘对数螺旋线破坏形态的几何要素关系，分析如下：

下缘挤压区：在△FBD中，有

${{OD} / {{\rm{sin}}\left( {\omega + \varphi '} \right)}} = {{BD} / {\sin\; {\delta _1}}}$

(28)

$OG = OD \cdot \cos \left( {{\delta _1} + \varphi '} \right) = \left( {OB - BC} \right) \cdot \cos\; \varphi '$

(29)

式中， $\omega = \arctan \;{\dfrac{{{y_{1n}}}}{{{{\textit{z}}_{\rm{w}}}\cos\;\beta }}} {\text{，}}$ $BC = \dfrac{{{{\textit{z}}_{\rm{w}}}\cos\; \beta }}{{\cos\; \varphi '}}{\text{，}}$ $BD =$ $ \sqrt {y_{1n}^2 + {{\left( {{{\textit{z}}_{\rm{w}}} \cdot \cos\; \beta } \right)}^2}} {\text{，}}\!\!$ OB=R₀，OD=R。

联立式（28）和（29）可求得下缘对数螺旋线的参数R₀、R和δ₁。

上缘张拉区：在 $\Delta O'C'D'$ 中，有

${{O'D'} / {{\rm{sin}}\left( {90^\circ - \varphi '} \right)}} = {{C'D'} / {\sin \;{\delta _3}}}$

(30)

$O'B'\cos\; \varphi ' - O'D'\cos \left( {\varphi ' - {\delta _3}} \right) = B'F'$

(31)

式中， $B'D' = \sqrt {y_{3n}^2 + {{\left( {{{\textit{z}}_{\rm{w}}}{\cos}\;\beta } \right)}^2}} {\text{，}}$ $O'D'$ = ${R_0}{\text{，}}$ $\omega ' =$ $ \arctan\; {\dfrac{{{y_{3n}}}}{{{{\textit{z}}_{\rm{w}}}\cos\; \beta }}}{\text{，}}$ $C'D' = {y_{3n}} - \dfrac{{{{\textit{z}}_{\rm{w}}}{\cos}\;\beta }}{{\cos \;\varphi '}}{\text{，}}$ $O'B'$ =R。

联立式（30）和（31）可求得下缘对数螺旋线的参数R₀、R和δ₃。

如图6所示，通过比较不同斜坡坡度及不同物理力学特性的土体，对数螺旋公式得到的上、下缘破裂面形态与定积分元素法所得结果非常接近，曲面形态也保持一致。因此，采用对数螺旋线描述滑动体上、下缘破裂面，能够较好地反映土体物理力学性质及边坡几何特性的影响。

3 边坡浅层稳定性分析

边坡浅层破坏失稳模式为3段组合滑动体，各自的受力情况如图8所示。张天宝^[26]认为挡墙背后填土与墙背之间存在相互作用时，墙背所受土压力合力作用点位于墙高1/3～1/2的范围内。为此，取挤压区p₁的作用点于z_cr/2处，取张拉区p₃的作用点于z_cr/3处。

如图8（a）、（b）所示，下缘挤压区和上缘张拉区的力矩平衡方程分别为：

$\begin{aligned}[b] & {{W_1'}}{l_1} + {p_1}{d_1} + {J_3}{m_3}- \\ & \;\;\;\;\;\;0.5c'R_0^2{\rm{cot}}\;\;\varphi '\left[ {\exp \left( {2{\theta _1}\tan \;\varphi '} \right) - 1} \right] = 0 \end{aligned}$

(32)

$\begin{aligned}[b] &{{W_3'}}{l_3} - {p_3}{d_3} + {J_1}\left( {{m_1}\cos\; \beta + {l_1}\sin \;\beta } \right) -\\ & \;\;\;\;\;\; 0.5c'R_0^2{\rm{cot}}\;\;\varphi '\left[ {\exp (2{\theta _3}\tan\; \varphi ') - 1} \right] = 0 \end{aligned}$

(33)

式中： ${W'_1}$ 和 ${W'_3}$ 分别为上、下缘的重力； ${J_1}$ 和 ${J_3}$ 为渗流力，对于非饱和土及不考虑渗流力时， ${J_1}$ 和 ${J_3}$ 均为0；l₁、d₁、m₁、l₃、d₃和m₃分别为相应的力臂长度。

中间主滑动体分别受到上缘土体的作用推力p₃和下缘土体的作用抗力p₁，由中间主滑动区的受力情况（图8（c））得到稳定安全系数的表达式为：

$FS = \frac{{\left[ {\left( {\gamma {\textit{z}}{{\cos }^2}\;\beta - {u_{\rm{a}}} - {\sigma _{\rm{s}}}} \right) \times \tan \;\varphi ' + c'} \right] \cdot {L_2} + {p_1}}}{{\gamma {L_2}\cos \;\beta \cdot \sin\; \beta + {p_3}}}$

(34)

当降雨强度大于入渗速度时，边坡浅层土体处于饱和状态，则主滑动体FS的表达式为：

$FS = \frac{{\left( {\gamma '{\textit{z}}{{\cos }^2}\;\beta \tan\; \varphi ' + c'} \right) \cdot {L_2} + {p_1}}}{{{\gamma _{{\rm{sat}}}}{\textit{z}}{\cos}\;\beta \sin\; \beta \cdot {L_2} + {p_3}}}$

(35)

$\gamma = \frac{{{G_{\rm{s}}} + S\!e}}{{1 + e}}{\gamma _{\rm{w}}} = {\gamma _{\rm{d}}} + S\!e \cdot n \cdot {\gamma _{\rm{w}}} = {\gamma _{\rm{d}}} + \theta {\gamma _{\rm{w}}}$

(36)

式中， $\gamma$ 、 $\gamma_{\rm sat}$ 和 $\gamma'$ 分别为土体的天然重度、饱和重度和有效重度，G_s为土体的比重，e为土体的孔隙比，n为土体的孔隙率。当L₂增至无穷大时（p₁和p₃可忽略不计），可得到无限长斜坡的稳定安全系数FS；当L₂趋近于0时，可得近似圆弧滑动时的FS，可见无限长斜坡和圆弧滑动的稳定分析方法只是“上、下缘顺坡曲面”组合滑动稳定分析方法的2种特例。对于饱和土质边坡，FS需按式（35）计算；对于非饱和土质边坡，FS需按式（34）计算。

当斜坡内深度z_w大于临界深度z_cr时，土体已发生屈服，难以得到滑动体上、下缘破裂面形态的精确解，此时可按z_w=z_cr时确定的破裂面近似计算。当z_w<z_cr时，中间主滑动区尚处于安全状态，并未产生张拉力，计算浅层边坡稳定时可不考虑上缘张拉区的裂缝。

4 分析与讨论 4.1 结果分析

土体的物理力学参数是边坡稳定分析的前提，采用文献[23]中某边坡实测数据，见表1。在分析某个参数的敏感性时，其他参数均保持不变。本文的“上、下缘顺坡曲面”组合滑动破坏模式与无限长斜坡滑动破坏模式的对比计算结果见图9。

由图9（a）可知，本文破坏模式计算得到的安全系数FS随主滑动区长度L₂的增加逐渐减小，且逐渐逼近于无限长斜坡破坏模式的计算结果。实际工程中斜坡长度有限，在z_w一定且L₂较小时，采用无限长斜坡稳定计算方法明显偏于保守，这与相关有限元分析^[12–13]得到的结论一致。原因在于，L₂/z_w较小时，整个滑坡体的稳定性主要取决于上缘张拉区和下缘挤压区的应力状态；反之，随着L₂/z_w增大，上缘张拉区和下缘挤压的应力状态对整个滑坡体稳定性的影响逐渐减弱。对于非饱和土质边坡，饱和度越小，其稳定性受上、下缘土体应力状态的影响越明显，采用无限长斜坡稳定计算方法得到的结果就更为保守。

中间主滑动区范围L₂=6.0 m时，对于不同饱和度土体，本文组合破坏与无限长斜坡破坏模式的对比结果如图9（b）所示。2种破坏模式得到的安全系数FS均随着滑坡深度z_w的增加呈指数衰减的趋势，且后者的计算结果明显偏小；当滑坡深度一定时，饱和度越小的土质斜坡的稳定性越好，且2种破坏模式的计算结果差距越大；饱和度越大的土质斜坡越容易发生浅层滑坡失稳。因此，对于有限长度的土质斜坡很有必要考虑滑坡体上、下缘承受的张拉和挤推作用，采用本文“上、下缘顺坡曲面”组合破坏模式的稳定方法分析斜坡的稳定性更为合理。

4.2 影响因素分析

在保持其他参数不变的情况下，图10分别给出了土质边坡稳定安全系数FS=1.0时，土体有效黏聚力c′、有效内摩擦角 $ \varphi'$ 、边坡坡度 $\beta$ 和饱和度Se等因素对中间主滑动区范围L₂和滑坡体深度z_w之间变化关系的影响。其中：Se=1.0时，土体为饱和状态；Se<1.0，土体为非饱和状态。由图10可知：L₂随着z_w的增加呈逐渐减小趋势，并趋近于0；当L₂近似为0时，整个滑动体可看作为近似圆弧滑动。

由图10（a）可知：当L₂保持恒定时，滑坡深度z_w随着c′增加逐渐增大；c′相同时，非饱和土质边坡对应的滑坡深度z_w要大得多。由图10（b）、（c）可知：对于饱和土质边坡， $\varphi'$ 和 $\,\beta$ 的变化对L₂与z_w关系的影响很小；相反，对于非饱和土质边坡， $\varphi'$ 和 $\,\beta$ 的变化对L₂与z_w关系的影响非常明显；当L₂为定值时， $\varphi'$ 越大对应的z_w越大， $\beta$ 越大对应的z_w越小， $ \varphi'$ （或 $\,\beta$ ）相同时，非饱和土质边坡对应的滑坡深度远大于饱和土质边坡。原因在于饱和土体会导致滑动面上的抗滑力严重减小，一定程度上弱化了 $\varphi'$ 和 $ \,\beta$ 的影响。由图10（d）可以看出，当L₂相同时，z_w随着土质斜坡土体饱和度的增加逐渐增大。以上分析再次说明，降雨浸润导致饱和度升高的土质边坡更容易发生浅层失稳破坏。

5 结　论

基于非饱和土强度理论，分析了斜坡内单元土体有效应力状态与极限Mohr圆的关系，提出由上缘张拉区、中间主滑动区和下缘挤压区组成的“上、下缘顺坡曲面”组合破坏模式，并开展相应的边坡稳定性分析方法的研究，探讨了土体强度参数、边坡坡度及土体饱和度等因素的影响，得到如下结论：

1）基于非饱和土强度理论，建立了边坡体内单元土体达到临界破坏状态时的判别方程，考虑滑坡体上、下缘的力学作用，采用对数螺旋线方程进行近似描述上、下缘土体剪切破坏曲面的形态，提出满足上、下缘土体力矩平衡和中间主滑动区静力平衡的边坡浅层稳定性分析计算方法，得到了相应的稳定安全系数的表达式。

2）本文组合破坏模式因考虑上、下缘土体的力学作用，使其得到的稳定安全系数FS更为准确合理，无限长斜坡失稳破坏模式得到结果明显偏小，且随着中间主滑动区长度的减小，2种破坏模式所得结果的差异性逐渐增大。

3）对于饱和土质边坡，土体有效黏聚力对其浅层稳定性及滑坡体范围的影响明显大于有效内摩擦角和边坡坡度；对于非饱和土体，内摩擦角和边坡坡度的变化对边坡浅层稳定性及滑坡体范围的影响反而更显著；随着边坡坡度和土体饱和度的增加，潜在危险滑动面逐渐向浅层迁移。

4）对于特定的土质边坡，随着土体饱和度的增加，对应的稳定安全系数逐渐减小，潜在滑坡体的深度随着饱和度的增加逐渐缩小，说明降雨入渗导致饱和度上升的土质边坡更容易发生浅层失稳破坏。