前言

0. 引言

天然土体在形成过程中由于矿物组成、沉积条件、应力历史的不同，土体参数呈现出明显的空间变异性，即土体物理力学参数的固有变异性，因而将边坡土体参数假定为空间上均质的传统确定性分析方法已不再适用。随着岩土工程可靠度理论的发展，边坡概率稳定性分析方法逐渐成为一种趋势，并且研究发现边坡的失效概率及破坏模式受到土体参数空间变异性的显著影响^[1-3]。为了科学合理地评估土体参数的空间变异性对岩土工程结构的影响，Vanmarcke^[4]率先提出利用随机场理论表征土体参数的空间变异性。此后，许多学者针对随机场的离散技术进行了有益的探索。

常见的随机场离散方法主要包括谱分解法、傅里叶变换法、中点法、局部平均法、KL级数展开法等^[5]。其中，较常使用的是中点法和KL级数展开法。中点法的计算效率较高且易于编程实现，但针对少数自相关矩阵非正定的情形，该方法无法对其进行乔列斯基分解，因而中点法存在一定的局限性。相比之下，KL（Karhunen-Loève）展开法不会受到上述制约，且该法不要求随机场网格与有限元网格完全对应^[6]。KL级数展开法的关键在于通过第二类Fredholm积分方程求解自相关函数的特征值与特征函数，该方法仅当自相关函数为三角型或指数型时第二类Fredholm积分方程才有解析解^[7]，其余情形均需通过数值方法求解该方程的特征值及特征函数，如利用伽辽金技术（Galerkin）近似求解该方程^[8]。基于Galerkin技术求解第二类Fredholm积分方程的近似解时，当前研究中多数学者通常将特征函数表示为一组基函数（如三角函数、小波函数等）的线性组合。目前利用KL展开法开展可靠度分析的研究均是直接通过该法针对特定岩土工程问题开展可靠度评估，鲜有学者对求解第二类Fredholm积分方程的方法进行探索。鉴于此，围绕建立第二类Fredholm积分方程的新型求解方法，简化KL展开法的计算过程是当前岩土工程可靠度评估领域多数学者所关注的热点。此外，潜在滑动面是边坡失稳时最可能的滑动位置，决定了滑体的分布范围及滑体体积的大小，对于评估边坡的失效后果、开展后续支护设计具有重要意义。然而通过有限元分析仅能提供一个剪切带，尚不能给出具体的滑面位置信息。通过等效塑性应变图人为主观判断滑面位置，会引入大量的不确定性，且存在效率低的不足。为此，亟需发展一种高效的滑面定位及滑体体积计算的一体化程序。

鉴于高斯型求积公式对数值积分具有较高代数精度的优点^[9]，本文提出一种基于Chebyshev-Galerkin- KL展开的新型随机场离散方法。首先，简单介绍了本文所提方法的基本原理及关键实现过程，然后基于Python提出了边坡滑面识别及滑体体积计算的高效方法。将本文所提方法生成的随机场与传统KL展开法的结果进行对比，验证本文所提方法的准确性。在此基础上通过一个边坡算例，探讨了所提方法在空间变异土质边坡稳定可靠度分析中的应用。

1. 基于Chebyshev-Galerkin-KL展开的随机场离散方法

1.1 KL级数展开法

随机场的生成，即随机场的离散过程，主要基于空间变异参数的统计特征，如概率分布类型、相关函数、波动范围、不同岩土体参数之间的互相关性等生成参数的随机场。本文首先介绍用于生成随机场的KL级数展开法原理及其计算过程。

若研究问题的随机参数为黏聚力c和内摩擦角φ，假定二者均服从正态分布，则基于KL展开的随机场为^{[1, 3]}

。

${{{H}}_{{i}}}({{x}}) \approx {\mu _{{i}}} + \sum\limits_{{{j}} = 1}^{{M}} {{\sigma _{{i}}}\sqrt {{\lambda _{{j}}}} } {{{f}}_{{j}}}({{x}}){\chi _{{{ji}}}} 。$

(1)

式中： ${{{H}}_{{i}}}({{x}})$ 为随机场离散结果；i=c或φ，x为研究区域的空间位置坐标； ${\mu _{{i}}}$ ， ${\sigma _{{i}}}$ 分别为随机参数的均值和标准差；M为随机场级数展开的项数； ${\lambda _{{j}}}$ ， ${{{f}}_{{j}}}({{x}})$ 分别为自相关函数的特征值与特征函数，通常由求解自相关函数的第二类Fredholm积分方程所得； ${\chi _{{{ji}}}}$ 为考虑不同随机场参数之间互相关性的随机向量矩阵。若随机参数服从对数正态分布，则相应的随机场仅需对式（1）求指数。

1.2 Galerkin技术

由式（1）可知，特征值 ${\lambda _{{j}}}$ 和特征函数 ${{{f}}_{{j}}}(x)$ 是离散随机场的关键，为此，引入Galerkin法求解第二类Fredholm积分方程：

$\int_{{\mathit{\Omega} }} {{\rm{Cov}}({{x}}, {{y}})} {{{f}}_{{j}}}(y){\rm{d}}{{y}} = {\lambda _{{j}}}{{{f}}_{{j}}}(x) 。$

(2)

式中：Cov(x, y)为自相关函数；Ω为积分区域。

常见的做法是将特征函数 ${{{f}}_{{j}}}(x)$ 近似表示为一组正交多项式的线性组合，

${{{f}}_{{j}}}(x) \approx {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{f} }}_{{j}}}(x) = \sum\limits_{{{k}} = 1}^{{N}} {\mathop {{d}}\nolimits_{{k}}^{\left( {{j}} \right)} } {\psi _{{{k}} - 1}}(x) = {\left[ {{d^{\left( {{j}} \right)}}} \right]^{\rm{T}}}\psi (x) 。$

(3)

式中：ψ_k-1(x)， $\mathop {{d}}\nolimits_{{k}}^{\left( {{j}} \right)}$ 分别为正交多项式及其对应的系数。

对比式（2），（3）可得式（2）的残差为

${{{r}}^{\left( {{j}} \right)}}(x) = \sum\limits_{{{k}} = 1}^{{N}} {{{d}}_{{k}}^{\left( {{j}} \right)}} \left[ {\int_{{\mathit{\Omega} }} {{\rm{Cov}}} (x, y){\psi _{{{k}} - 1}}(y){\rm{d}}y - {\lambda _{{j}}}{\psi _{{{k}} - 1}}(x)} \right] 。$

(4)

假定残差与基函数正交^[3]，则

$\left\langle {{r^{\left( {{j}} \right)}}(x), {\psi _{{{i}} - 1}}(x)} \right\rangle = 0\quad ({{i}} = 1, \cdots , {{N}}) 。$

(5)

根据式（4），（5）可进一步表示为

$\begin{array}{l} \sum\limits_{{{k}} = 1}^{{N}} {{{d}}_{{k}}^{\left( {{j}} \right)}} \int_{{\mathit{\Omega} }} {\int_{{\mathit{\Omega} }} {{\rm{Cov}}(x, y){\psi _{{{k}} - 1}}(y){\psi _{{{i}} - 1}}(x){\rm{d}}x{\rm{d}}y} } -\\ \ \ \ \ \ \ \ \ \ \ \ {\lambda _{{j}}}\sum\limits_{{{k}} = 1}^{{N}} {{{d}}_{{k}}^{\left( {{j}} \right)}} \int_{{\mathit{\Omega} }} {{\psi _{{{k}} - 1}}(x){\psi _{{{i}} - 1}}(x){\rm{d}}x} = 0 。 \end{array}$

(6)

记

${A_{{{ik}}}} = \int_{{\mathit{\Omega} }} {\int_{{\mathit{\Omega} }} {{\rm{Cov}}(x, y){\psi _{{{k}} - 1}}(y){\psi _{{{i}} - 1}}(x){\rm{d}}x{\rm{d}}y} } {, }$

(7)

$\left. {\begin{array}{*{20}{l}} {{D_{kj}} = d_k^{\left( j \right)},}\\ {\Lambda = {\rm{diag}}\left[ {{\lambda _1},{\lambda _2}, \cdots ,{\lambda _N}} \right],} \end{array}} \right\}$

(8)

${B_{{{ik}}}} = \int_{{\mathit{\Omega} }} {{\psi _{{{k}} - 1}}(x){\psi _{{{i}} - 1}}(x){\rm{d}}x} 。$

(9)

根据式（7）~（9），式（6）可简化为如下矩阵形式：

${\boldsymbol{AD}} = {\boldsymbol{BD\varLambda}} 。$

(10)

本文采用第二类Chebyshev多项式作为正交基函数

${U}_{{n}}(x)=\frac{\sin\left[({n}+1){\rm{arccos}}x\right]}{\sqrt{1-{x}^{2}}}\ \ \ \ ({n=}1{, }2{, }\cdots ) 。$

(11)

将式（7）积分区间变化至[-1, 1]，则式（7）可表示为

$\begin{array}{l} {A_{{{ik}}}} \approx {S_x}{S_y}\sum\limits_{{{{l}}_x} = 1}^{{{{n}}_x}} {\sum\limits_{{{{l}}_y} = 1}^{{{{n}}_y}} {{w_{{{{l}}_x}}}{w_{{{{l}}_{_y}}}}{\rm{Cov}}({S_x}{t_{{{{l}}_x}}} + {T_{{x}}}, {S_y}{t_{{{{l}}_y}}} + {T_y})} } \cdot\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {U_{{{k}} - 1}}({S_y}{t_{{{{l}}_y}}} + {T_y}){U_{{{i}} - 1}}({S_x}{t_{{{{l}}_x}}} + {T_x}) 。 \end{array}$

(12)

式中： ${w_{{{{l}}_x}}}$ ， ${w_{{{{l}}_{_y}}}}$ 与 ${t_{{{{l}}_x}}}$ ， ${t_{{{{l}}_y}}}$ 分别为Gauss-Legendre积分的权重和节点。在编程计算该积分时，取n_x=n_y=min{16，4M}^[10]。S_x，S_y为缩放系数，T_x，T_y为平移系数，分别可由下式计算：

$\left. \begin{array}{l}{S}_{x}=({x}_{max}-{x}_{min})/2{ }{, }\\ {T}_{x}=({x}_{max}+{x}_{min})/2{ }。\end{array} \right\}$

(13)

类似可得S_y，T_y。原理同式（12），矩阵B可如下简化计算：

$\begin{array}{l} {B_{{{ik}}}}{{ = }}\int_{{\mathit{\Omega} }} {{U_{{{k}} - 1}}(x){U_{{{i}} - 1}}(x)} {\rm{d}}x \\ \approx \sum\limits_{{{{l}}_{{x}}} = 1}^{{{{n}}_{{x}}}} {{S_x}{U_{{{k}} - 1}}({S_x}{t_{{{{l}}_{{x}}}}} + {T_x}){U_{{{i}} - 1}}({S_x}{t_{{{{l}}_{{x}}}}} + {T_x})} 。\end{array}$

(14)

对一具体问题，一旦矩阵A和B确定，特征函数矩阵D和特征值矩阵Λ便可相应确定（即通过MATLAB中的函数eigs求解）。

由于Galerkin法近似得到的特征函数往往并非正交的，为此引入一个归一化系数 ${{{c}}_{{j}}}$ ，

$\begin{array}{l} {{{c}}_{{j}}} = \int_{{\mathit{\Omega} }} {{{\widehat {{f}}}_{{j}}}^2(x)} {\rm{d}}x = \sum\limits_{{{m}} = 1}^{{N}} {\sum\limits_{{{n}} = 1}^{{N}} {\mathop {{d}}\nolimits_{{m}}^{\left( {{j}} \right)} } } \mathop {{d}}\nolimits_{{n}}^{\left( {{j}} \right)} \int_{{\mathit{\Omega} }} {{U_{{{m}} - 1}}(x)} {U_{{{n}} - 1}}(x){\rm{d}}x \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{ = }}{\left[ {{d^{\left( {{j}} \right)}}} \right]^{\rm{T}}}{\boldsymbol{B}}{d^{\left( {{j}} \right)}} 。 \end{array}$

(15)

基于式（15），可得归一化后的特征函数为

$\widehat {{f}}_{{j}}^*(x) = \frac{1}{{\sqrt {{{{c}}_{{j}}}} }} \cdot \widehat {{f}}(x) = \frac{{{{\left[ {{d^{\left( {{j}} \right)}}} \right]}^{\bf{T}}}U(x)}}{{\sqrt {{{\left[ {{d^{\left( {{j}} \right)}}} \right]}^{\bf{T}}}{\boldsymbol{B}}{d^{\left( {{j}} \right)}}} }} 。$

(16)

由式（10）~（16），可得基于Chebyshev-Galerkin- KL展开的一维相关高斯型随机场：

${{{H}}_{{i}}}\left( {{x}} \right) \approx {\mu _{{i}}} + {\sigma _{{i}}}\sum\limits_{{{j}} = 1}^{{M}} {\frac{{\sqrt {{\lambda _{{j}}}} {{\left[ {{d^{\left( {{j}} \right)}}} \right]}^{\rm{T}}}U(x)}}{{\sqrt {{{\left[ {{d^{\left( {{j}} \right)}}} \right]}^{\rm{T}}}{\boldsymbol{B}}{d^{\left( {{j}} \right)}}} }}} {\chi _{{{ji}}}} 。$

(17)

对于二维可分离型自相关函数，其特征函数、特征值仅需对一维x，y方向的结果相乘，由此二维相关高斯随机场便可相应求得。

1.3 Chebyshev-Galerkin-KL展开项数的确定

很显然，M的值越大越有利于准确离散随机场，实际应用时考虑到计算效率，通常仅取前M项较大的特征值及其对应的特征函数对随机场进行离散。为保证随机场离散的精度，参考文献[3，11]采用逐点均方误差绝对值的平均值来衡量随机场离散结果的准确性。

$\varepsilon = \frac{1}{{{{{n}}_{\rm{e}}}}}\sum\limits_{{{i}} = 1}^{{{{n}}_{\rm{e}}}} {\left| {\varepsilon ({x_{{i}}})} \right|} 。$

(18)

式中：n_e为研究问题在空间的离散点个数；ɛ(x_i)为空间离散点x_i处的均方误差，

$\varepsilon ({x_{{i}}}) = 1 - \sum\limits_{{{p}} = 1}^{{M}} {{\lambda _{{p}}}} {\widehat {{f}}_{{j}}}^2({x_{{p}}}) 。$

(19)

2. 基于Python的滑体体积计算及失效模式识别方法

边坡一旦发生垮塌，必会对国民经济和人民生命财产安全产生重要影响^[12]，因此开展边坡失效后果评估，对于应急管理部门确定风险等级，制定抢险救灾方案具有重要意义。此外，基于本文所提的随机场生成方法开展有限元计算，仅能提供边坡的安全系数、等效塑性应变云图、位移云图等信息，无法直接输出临界滑动面的位置信息以及相应的滑体体积，且对任一目标问题，通常需要开展几千甚至更多次随机有限元分析，直接对每个结果文件进行后处理也是不现实的，且不利于工程推广。考虑到Python语言运行效率高，本文研发了基于Python的滑面识别及滑体体积计算的程序，其主要原理如下：

以ABAQUS开展随机有限元分析为例。首先，编写Python脚本，根据ABAQUS计算的结果文件，提取边坡全局的节点信息及相应的位移场；其次，由于边坡处于破坏或临界状态时，滑体节点的位移往往显著大于下部滑床的节点位移，为此，本文沿用Huang等^[13]提出利用k-means完成滑面识别的方法，在Python环境下实现了基于该算法对滑体区域的自动识别，并利用Python自带的“ConvexHull”函数提取了滑体轮廓点坐标；然后，根据向量叉乘原理实现对滑体凸多边形面积的计算（即二维边坡单位长度内的滑体体积）。此外，本文基于Python语言，定义了边坡深层破坏和浅层破坏两类失效模式。具体来讲，当滑面滑出点位于坡脚以外或者滑面最低点位于坡脚水位线以下时，视为深层破坏；反之，则视为浅层破坏。

3. 基于Chebyshev-Galerkin-KL展开法的可靠度分析流程

图 1给出了基于Chebyshev-Galerkin-KL展开的可靠度分析流程，该流程主要包括6步。

图 1 基于Chebyshev-Galerkin-KL展开的边坡稳定可靠度分析流程图

Figure 1. Flow chart of reliability analysis of slope stability using Chebyshev-Galerkin-KL expansion

下载: 全尺寸图片幻灯片

（1）确定边坡的几何尺寸，空间变异岩土体参数的统计特征（如均值、变异系数、概率分布类型等）、相关函数以及波动范围。

（2）建立边坡有限元模型，提取网格中心点坐标，用于随机场的离散，并确定S_x，S_y，T_x，T_y。

（3）对M赋初值，令M=1，根据式（12），（14）计算矩阵A和B，由式（10）求解矩阵D和Λ；在每个随机场网格的中心点处，由式（15），（16）计算特征函数 ${{\hat f}}_{{j}}^*(x)$ 。

（4）由式(18)，(19)，计算随机场的离散误差ɛ，若ɛ小于规定允差ɛ₀（如ɛ₀=5%），表明已满足随机场离散的精度要求，否则令M=M+1，重复步骤(3)，(4)，直至ɛ < ɛ₀。

（5）由式（17），生成服从正态分布的随机场，若随机参数服从对数正态分布，则对式（17）取指数即可得到目标随机场。

（6）将生成的随机场映射到有限元软件ABAQUS中，由此开展边坡稳定性分析，评估边坡的安全系数并基于自主研发的Python程序完成边坡滑体体积的计算与失效模式的识别。

4. 算例分析

4.1 算例介绍

本文以黏性/摩擦土质边坡为例开展可靠度分析，验证提出方法的有效性。图 2为该边坡的有限元模型，其坡高为20 m，坡角为40°，初始水位线的高度为15 m，在水位上涨期间，假定最大水位高度h_max为25 m。采用ABAQUS开展边坡渗流及稳定性分析，边坡网格类型为四节点平面应变四边形孔隙流体/应力单元（CPE4P），共计划分网格3244个，节点3375个。在进行有限元分析时，边坡左右边界设为水平约束，底部边界设为固定约束。边坡的岩土体参数如，本研究仅考虑有效黏聚力 $c'$ 、有效内摩擦角 $\varphi '$ 及饱和渗透系数k_s的空间变异性，并假定具有相同的相关结构（本文采用高斯型自相关函数），水平波动范围δ_h及竖直波动范围δ_v分别为50，10 m，且 $c'$ 和 $\varphi '$ 的互相关系数 ${\rho _{{\rm{c}}'\varphi '}}$ 为-0.5^[14]。

图 2 边坡几何尺寸及其有限元网格

Figure 2. Geometry and finite element mesh of slope

下载: 全尺寸图片幻灯片

表 1 土体参数的统计特征

Table 1. Statistics of soil properties

参数	均值	变异系数	分布类型
E	100 MPa	—	—
$\nu$	0.3	—	—
$\gamma$	13 kN/m³	—	—
$c'$	50 kPa	0.3	对数正态分布
$\varphi '$	10^°	0.2	对数正态分布
k_s	5×10^-6 m/s	0.4	对数正态分布
注：E， $\nu$ ， $\gamma$ 分别为弹性模量、泊松比和土体重度。

下载: 导出CSV

| 显示表格

需说明的是，在开展渗流分析时采用Fredlund- Xing模型^[15]来刻画土体的有效饱和度S_e随基质吸力ψ的变化关系。

${{{S}}_{\rm{e}}} = \frac{{\theta - {\theta _{\rm{r}}}}}{{{\theta _{\rm{s}}} - {\theta _{\rm{r}}}}} = \frac{1}{{{{\left\{ {\ln \left[ {{\rm{e}} + {{(\psi \cdot {{a}})}^{{n}}}} \right]} \right\}}^{{m}}}}} 。$

(20)

式中：θ为体积含水量； ${\theta _{\rm{s}}}$ ， ${\theta _{\rm{r}}}$ 分别为土体的饱和体积含水量和残余体积含水量；e为自然底数；a，m，n该模型的拟合参数，各参数的取值见表 2^[16]。

表 2 Fredlund-Xing模型参数^[17]

Table 2. Parameters for Fredlund-Xing model^[17]

参数	取值	参数	取值
θ_s	0.387	α/kPa^-1	0.643
θ_r	0.032	n	1.353
S_r	0.083	m	0.261
注：S_r为残余饱和度。

下载: 导出CSV

| 显示表格

渗透系数k的方程为^[17]

${{k}} = {{{k}}_{\rm{s}}}{({{{S}}_{\rm{e}}})^\alpha } 。$

(21)

式中：α为常数，本研究取ABAQUS中默认值3。此外，在考虑渗流条件下开展边坡稳定性分析时，本构模型采用扩展莫尔-库仑准则 ^[18]：

${\tau _{\rm{f}}}{{ = c'}} + (\sigma - {u_{\rm{a}}})\tan \varphi ' + \chi ({u_{\rm{a}}} - {u_{\rm{w}}})\tan \varphi ' 。$

(22)

式中： $\sigma$ 为总应力；χ为有效应力参数，本研究中χ=S_e；u_a，u_w为孔隙气压力和孔隙水压力，因与大气连通，取u_a=0。

4.2 结果验证

为了验证本研究所提随机场离散方法在表征岩土体参数空间变异性方面的准确性，本节以图 2边坡为例，基于Chebyshev-Galerkin-KL法生成随机场的离散误差随展开项数M的变化关系如图 3。可见，对于高斯型相关结构，波动范围分别为δ_h=50 m，δ_v=10 m时，级数展开项数M为6可保证本算例中随机场的离散误差（ɛ=2.54%）小于规定允差（ɛ₀=5%）。同时文献[3]采用传统KL展开法离散的随机场误差如图 3中虚线所示，总体而言，两种随机场表征方法的离散误差较为一致。上述结果说明本文所提Chebyshev- Galerkin-KL展开法在表征岩土体参数空间变异性方面是有效的。为进一步说明本文所提方法在离散随机场方面的可靠性，图 4以有效内摩擦为例，对比了上述两种方法离散的随机场结果。根据图 4（a），（b）可知，两种方法生成的随机场结果较接近，且图 4（c）揭示了这两种方法生成随机场的误差绝对值较小。因此，本文所提Chebyshev-Galerkin-KL展开法在离散随机场时具有较高的准确性。

图 3 Chebyshev-Galerkin-KL展开法的离散误差与M的关系

Figure 3. Relationship between discretization error (ɛ) of Chebyshev- Galerkin-KL expansion method and M

下载: 全尺寸图片幻灯片

图 4 不同随机场离散方法生成的有效内摩擦角对比图

Figure 4. Comparison between random fields of effective friction angle by different discretization methods

下载: 全尺寸图片幻灯片

图 5为水位以0.4 m/d的速度上涨至最高水位25 m，所有岩土体参数均取其均值时，ABAQUS基于强度折减法得到的等效塑性应变云图，其能够定性地反映边坡临界滑动面的位置且相应的边坡安全系数F_S为1.250。此外，图 5中红色曲线为基于本文研发的Python脚本及有限元计算的位移场识别的潜在滑动面示意图。可以看出，本研究基于k-means算法研发的Python脚本可准确识别滑体区域及其滑面位置，并且还可以自动计算该边坡的单位长度内滑动体积约为511.48 m³/m（后文以m²计）。另一方面，该Python脚本识别的边坡破坏模式为深层破坏，显然与有限元的计算结果一致，在配置为Intel(R) Core(TM) i5-8500 CPU@3.00GHz的电脑上基于该脚本对1000次随机有限元计算文件提取滑体体积与识别失效模式耗时约15 s。可见，本文研发的Python脚本在边坡风险评估方面是准确高效的。

图 5 FEM的计算得到的滑面及基于Python提取的滑体信息

Figure 5. Slip surface via FEM and sliding area obtained by Python script

下载: 全尺寸图片幻灯片

4.3 边坡可靠度分析

（1）水位上升速度的影响

本节将研究不同水位上升速度对边坡可靠度分析结果的影响。为此，计算了4种水位上升速度（如0.1，0.2，0.3，0.4 m/d）至最高水位对应的边坡安全系数、失效模式和滑体体积，结果见图 6。

图 6 水位上升速度对边坡稳定性的影响

Figure 6. Effects of different water rising velocities on slope stability

下载: 全尺寸图片幻灯片

图 6为基于1000次蒙特卡洛模拟（MCS）得到的安全系数累积分布图（CDF）^[19]，很显然不同水位上升速度对应的边坡可靠度分析结果存在一定差异。水位上升速度越大，边坡安全系数的CDF曲线朝FS增大的方向移动，这是因为库水渗入坡体存在滞后效应^[20]，上升速度越大，滞后效应越显著，使得坡体中非饱和区域更大，导致非饱和边坡土体的抗剪强度增大，从而有利于边坡稳定。此外，由图 6还可以发现，4种工况下基于可靠度评估得到的安全系数最小值均大于1.10，满足《水利水电工程边设计规范》^[21]对5级边坡稳定性的规定（即要求F_S达到1.10~1.05），由此判断该边坡在25 m高的水位下基本是稳定的。

图 7比较了不同水位上升速度时基于随机有限元计算得到的滑体体积最大值和最小值以及失效模式。可以看出，不同水位上升速度几乎对滑体体积和失效模式的影响较为微弱，该边坡可能的滑动体积基本稳定在270~750 m²，且均为深层滑坡。

图 7 水位上升速度对滑体体积及边坡失效模式的影响

Figure 7. Effects of different water rising velocities on sliding volume and failure mode

下载: 全尺寸图片幻灯片

（2）最大水位高度的影响

汛期，库岸边坡通常会面临水位上涨的现象，库水对岸坡岩土体的浸泡是造成其强度发生劣化或衰减的重要原因之一，有必要研究最大水位高度h_max对边坡稳定可靠度分析结果的影响。图 8为3种不同最大水位高度（h_max为20，25，30 m）时边坡安全系数的CDF图。需要指出的是，由上节的分析可知，当水位上升速度v=0.1 m/d时，边坡的安全系数整体上更倾向于较小值，为此，本节探讨的3种工况均设定水位上升速度为0.1 m/d。

图 8 最大水位高度对边坡稳定性的影响

Figure 8. Effects of maximum water level on slope stability

下载: 全尺寸图片幻灯片

由图 8可知，最大水位高度对边坡的稳定性具有显著影响。可靠度评估结果表明：库水位最大高度增大将导致边坡安全系数朝F_S减小的方向移动。当h_max由20 m上升至30 m时，边坡最小安全系数由1.3左右骤降至0.9以下。因此，在汛期须为库岸边坡设计相应的泄洪防涝措施。

图 9给出了h_max取不同值时边坡的滑体体积及失效模式，同样可以得出该边坡在水位上涨期间会发生深层滑坡的结论。类似图 7，1000次MCS得到的边坡滑体体积最大值与最小值之间存在显著差异，这意味着忽略土体参数的空间变异性将会对滑体体积产生不准确的估计。此外，观察图 9，还可以发现：当h_max由20 m上升至30 m时，边坡滑体体积整体上呈减小的趋势。这是由于相较于水位较低时（如20，25 m），30 m的水位会产生更大的静水压力作用在边坡坡面上，在一定程度上抑制了滑动面向深部发展。

图 9 不同最大水位高度对滑体体积及边坡失效模式的影响

Figure 9. Effects of different maximum water level on sliding volume and failure mode

下载: 全尺寸图片幻灯片

（3）波动范围的影响

众多学者研究发现，岩土体参数在水平方向的空间变异性并不会对边坡可靠度评估结果产生显著影响^[22]。为此，笔者在本节重点探讨其竖向空间变异性对边坡可靠度分析结果的影响，共计考虑了3种工况，各工况下的水平波动范围δ_h均为50 m，竖向波动范围δ_v依次为5，10，20 m，其余参数同4.1节。

图 10给出了F_S的CDF随δ_v的变化曲线。可以看出，当δ_v为5 m时，F_S的分布范围较为集中，随着δ_v增大，F_S的分布范围变宽，但总体而言，δ_v对该边坡在水位上升条件下F_S的分布影响并不显著。同理，图 11比较了不同竖向波动范围时1000次MCS计算的滑体体积最小值与最大值。显然，随着δ_v的增大，滑体体积的最大值或最小值均呈现出逐渐增大的趋势^[22]，且在每个δ_v下最大滑体体积约为最小滑体体积的3倍。以上结果表明土体参数的竖向空间变异性对该边坡的失效风险具有显著影响，因而在评估库岸边坡的失效风险时不可忽略土体参数的竖向空间变异性。另外，观察图 11还可以发现，土体参数的竖向空间变异性并不会影响该边坡的失效模式，其恒为深层失稳。

图 10 不同竖向波动范围δ_v对边坡稳定性的影响

Figure 10. Effects of different vertical scales of fluctuation δ_v on slope stability

下载: 全尺寸图片幻灯片

图 11 不同竖向波动范围对滑体体积及边坡失效模式的影响

Figure 11. Effects of different vertical scales of fluctuation on sliding volume and failure mode

下载: 全尺寸图片幻灯片

（4）相关系数的影响

下面探讨土体参数的互相关性对水位上升条件下边坡可靠度的影响规律，本节以 $c'$ ， $\varphi '$ 间的互相关性 ${\rho _{c'\varphi '}}$ 为例进行阐述。需要指出的是，本节考虑 ${\rho _{c'\varphi '}}$ 取-0.5，0，0.5等3种情况进行计算，其余参数均为4.1节的参考值。给出了 ${\rho _{c'\varphi '}}$ 取上述3个值时边坡安全系数的CDF图。可以看出，随着 ${\rho _{c'\varphi '}}$ 的增大，该边坡安全系数的分布区间逐渐变窄，且最大安全系数呈现出减小的趋势，可见忽略岩土体强度参数间的负相关性，将会低估边坡的稳定性，从而造成偏保守的设计。比较了不同 ${\rho _{c'\varphi '}}$ 所得到的边坡滑体体积最大、最小值以及相应的失效模式。不难发现，边坡的最小滑体体积亦随 ${\rho _{c'\varphi '}}$ 呈现出逐渐增大的趋势，这进一步说明了忽略岩土体强度参数间的负相关性将会高估边坡的失效风险，但最大滑体体积受 ${\rho _{c'\varphi '}}$ 的影响较为微弱，基本稳定在740 m²左右。此外，由图 13还可以发现，该算例边坡的失效模式不会受岩土体强度参数间相关系数的影响，均为深层滑移失效。

图 12 不同相关系数

${\rho _{c'\varphi '}}$ 对边坡稳定性的影响

Figure 12. Effects of different cross-correlation coefficients

${\rho _{c'\varphi '}}$ on slope stability

下载: 全尺寸图片幻灯片

图 13 不同相关系数

${\rho _{c'\varphi '}}$ 对滑体体积及边坡失效模式的影响

Figure 13. Effects of different cross-correlation coefficients

${\rho _{c'\varphi '}}$ on sliding volume and failure mode

下载: 全尺寸图片幻灯片

5. 结论

针对岩土体参数具有空间变异性的特点，及利用KL展开法离散随机场时难以直接求解第二类Fredholm积分方程、中点法存在相关系数矩阵非正定时无法对其进行乔列斯基分解等难点，基于Galerkin技术和第二类Chebyshev多项式，提出了新的随机场离散方法，并建立了基于Python语言自动进行边坡滑体体积计算及失效模式识别的高效程序。通过水位上涨下的边坡算例探讨了水位上升速度、最大水位高度、岩土体参数的竖向波动范围以及相关系数对边坡稳定可靠度分析结果的影响，得到4点主要结论。

（1）提出基于Chebyshev-Galerkin-KL展开的随机场离散方法，为求解第二类Fredholm积分方程提供了一种新的思路，且能准确地表征岩土体参数的空间变异性。

（2）基于Python研发的可靠度分析结果后处理程序，具有高效求解滑体体积以及准确识别边坡失效模式的优点。此外，该程序与随机有限元分析并不耦合，显著提升该程序的实用性。

（3）水位的上升速度、最大水位高度对算例边坡的可靠度分析结果具有显著影响。涨速越小，最大水位高度越高，该边坡的安全程度越低。水位上升速度对滑体体积的影响较为微弱，相较之下，最大水位高度的影响更为突出。岩土体参数的竖向空间变异性几乎对该边坡安全系数的分布特点无影响，但滑体体积受其影响较为明显。此外，抗剪强度参数间的相关性亦会影响该边坡的稳定性，忽略 $c'$ ， $\varphi '$ 间的负相关性将会低估边坡的稳定性。

（4）通过参数敏感性分析，本研究研发的边坡失效模式自动识别程序表明，算例边坡在汛期一旦发生滑坡，将会造成深层失稳。

开展了水位上升条件下考虑岩土体参数空间变异性的边坡稳定可靠度分析，对于水位下降以及水位周期性涨落时的边坡稳定可靠度分析有待于进一步研究。

期刊类型引用(3)

1.	徐小兵. 改进粒子群算法+BP神经网络在边坡可靠度分析中的应用. 黑龙江交通科技. 2024(08): 41-45 . 百度学术
2.	汪卫，廖志浩，廖杰林. 基于随机场模型的土石坝失效概率分析研究进展. 人民珠江. 2024(10): 114-123 . 百度学术
3.	丁盛环. 基于单根锚杆失效的边坡稳定性可靠度分析. 四川建材. 2024(12): 105-106+112 . 百度学术

前言 English Version

计量

0. 引言

1. 基于Chebyshev-Galerkin-KL展开的随机场离散方法

1.1 KL级数展开法

1.2 Galerkin技术

1.3 Chebyshev-Galerkin-KL展开项数的确定

2. 基于Python的滑体体积计算及失效模式识别方法

3. 基于Chebyshev-Galerkin-KL展开法的可靠度分析流程

4. 算例分析

4.1 算例介绍

4.2 结果验证

4.3 边坡可靠度分析

5. 结论

期刊类型引用(3)

其他类型引用(2)

计量

目录

0. 引言

1. 基于Chebyshev-Galerkin-KL展开的随机场离散方法

1.1 KL级数展开法

1.2 Galerkin技术

1.3 Chebyshev-Galerkin-KL展开项数的确定

2. 基于Python的滑体体积计算及失效模式识别方法

3. 基于Chebyshev-Galerkin-KL展开法的可靠度分析流程

4. 算例分析

4.1 算例介绍

4.2 结果验证

4.3 边坡可靠度分析

5. 结论

前言 English Version

计量

出版历程

0. 引言

1. 基于Chebyshev-Galerkin-KL展开的随机场离散方法

1.1 KL级数展开法

1.2 Galerkin技术

1.3 Chebyshev-Galerkin-KL展开项数的确定

2. 基于Python的滑体体积计算及失效模式识别方法

3. 基于Chebyshev-Galerkin-KL展开法的可靠度分析流程

4. 算例分析

4.1 算例介绍

4.2 结果验证

4.3 边坡可靠度分析

5. 结论

期刊类型引用(3)

其他类型引用(2)

计量

出版历程

目录

0. 引言

1. 基于Chebyshev-Galerkin-KL展开的随机场离散方法

1.1 KL级数展开法

1.2 Galerkin技术

1.3 Chebyshev-Galerkin-KL展开项数的确定

2. 基于Python的滑体体积计算及失效模式识别方法

3. 基于Chebyshev-Galerkin-KL展开法的可靠度分析流程

4. 算例分析

4.1 算例介绍

4.2 结果验证

4.3 边坡可靠度分析

5. 结论