基于Chebyshev-Galerkin-KL展开的土质边坡稳定可靠度分析

顾鑫; 仉文岗; 欧强; 王林; 覃长兵

doi:10.11779/CJGE20220831

基于Chebyshev-Galerkin-KL展开的土质边坡稳定可靠度分析 English Version

顾鑫^1,,
仉文岗^{1, 2, 3, ,},
欧强^{1, 2, 3},
王林⁴,
覃长兵^{1, 2, 3}

1.
重庆大学土木工程学院, 重庆 400045
2.
重庆大学山地城镇建设与新技术教育部重点实验室, 重庆 400045
3.
重庆大学库区环境地质灾害防治国家地方联合工程研究中心, 重庆 400045
4.
北京师范大学国家安全与应急管理学院, 广东珠海 519087

基金项目:

国家重点研发计划重点专项 2019YFC1509605

重庆市研究生科研创新项目 CYS22054

国家自然科学基金青年基金项目 52008058

国家自然科学基金青年基金项目 52108299

外国文教专家项目 G2022165002

详细信息

作者简介:
顾鑫(1995—)，男，博士研究生，主要从事岩土工程可靠度分析与风险控制方面的研究。E-mail: guxin@cqu.edu.cn

通讯作者:
仉文岗, E-mail: zhangwg@cqu.edu.cn

中图分类号: TU43
计量
- 文章访问数: 468
- HTML全文浏览量: 57
- PDF下载量: 187
出版历程
- 收稿日期: 2022-07-02
- 网络出版日期: 2023-03-05
- 刊出日期: 2023-11-30

Reliability analysis of soil slope stability based on Chebyshev-Galerkin-KL expansion

GU Xin^1,,
ZHANG Wengang^{1, 2, 3, ,},
OU Qiang^{1, 2, 3},
WANG Lin⁴,
QIN Changbing^{1, 2, 3}

1.
School of Civil Engineering, Chongqing University, Chongqing 400045, China
2.
Key Laboratory of New Technology for Construction of Cities in Mountain Area, Chongqing University, Chongqing 400045, China
3.
National Joint Engineering Research Center of Geohazards Prevention in the Reservoir Areas, Chongqing University, Chongqing 400045, China
4.
School of National Safety and Emergency Management, Beijing Normal University, Zhuhai 519087, China

摘要

摘要: 提出了一种基于Chebyshev-Galerkin-KL（Karhunen-Loève）展开的新型随机场离散方法，推导了该方法的理论公式，并研发了基于Python语言的边坡滑体体积计算及其失效模式自动识别的高效程序。采用一个水位上涨时的非饱和土坡算例验证了方法的有效性。结果表明：所提随机场生成方法为求解第二类Fredholm积分方程提供了一种新思路，可准确表征岩土体参数的空间变异性。基于Python提出的边坡风险评估程序，与随机有限元计算过程不耦合，极大地提升了开展边坡风险评估的效率，从而增强预测滑坡致灾风险的时效性。此外，研究中非饱和土坡算例可靠度分析结果表明：水位上涨速度越慢，最大水位高度越高，该边坡的安全程度将越低。岩土体参数的竖向空间变异性对该算例边坡安全系数的影响甚微，但对滑体体积的影响较为显著。在对水位上升条件下的边坡开展可靠度分析时，应关注抗剪强度参数间的负相关性，否则将会低估边坡的安全程度。
- Chebyshev-Galerkin-KL展开法 /
- 空间变异性 /
- 边坡稳定性 /
- 可靠度分析 /
- 失效风险
Abstract: A novel method is put forward for the random field discretization based on the Chebyshev-Galerkin-KL (Karhunen-Loève) expansion, followed by the derivation of equations for the proposed method. By means of Python language, an efficient program is exploited for automatically calculating the slope sliding volume and identifying the slope failure mode. The proposed method is validated through an unsaturated slope example subjected to water rising. The results indicate that the proposed method for the random field generation provides a new way to solve the Fredholm integral equation of the second kind, which can accurately characterize the spatial variability of geotechnical parameters. The Python-based program for risk estimation is decoupled from the random finite element calculations, which ensures the slope risk estimation with sufficient efficiency and promotes the reduction of the required time to predict the landslide risk. In addition, the obtained results from the unsaturated slope example show that a lower water rising velocity and a greater maximum water level will lead to the decrease of slope stability. The vertical spatial variability of geotechnical parameters has marginal effects on the safety factor of slopes. However, the sliding volume may be significantly affected. Attention should be paid to the negative cross-correlation between shear strength parameters when conducting reliability analysis of slope stability under water rising. Otherwise, the slope stability will be underestimated.
- Chebyshev-Galerkin-KL expansion /
- spatial variability /
- slope stability /
- reliability analysis /
- risk of failure

HTML全文

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

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

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

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

下载: 全尺寸图片幻灯片

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

Figure 2. Geometry and finite element mesh of slope

下载: 全尺寸图片幻灯片

图 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 FEM的计算得到的滑面及基于Python提取的滑体信息

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

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

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

Figure 8. Effects of maximum water level on slope stability

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

图 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

下载: 全尺寸图片幻灯片

不同相关系数 ${\rho _{c'\varphi '}}$ 对边坡稳定性的影响

Effects of different cross-correlation coefficients ${\rho _{c'\varphi '}}$ on slope stability

下载: 全尺寸图片幻灯片

不同相关系数 ${\rho _{c'\varphi '}}$ 对滑体体积及边坡失效模式的影响

Effects of different cross-correlation coefficients ${\rho _{c'\varphi '}}$ on sliding volume and failure mode

下载: 全尺寸图片幻灯片

表 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

表 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

参考文献(22)

[1]	李典庆, 蒋水华, 周创兵, 等. 考虑参数空间变异性的边坡可靠度分析非侵入式随机有限元法[J]. 岩土工程学报, 2013, 35(8): 1413-1422. http://www.cgejournal.com/cn/article/id/15248 LI Dingdian, JIANG Shuihua, ZHOU Chuangbing, et al. Reliability analysis of slopes considering spatial variability of soil parameters using non-intrusive stochastic finite element method[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(8): 1413-1422. (in Chinese) http://www.cgejournal.com/cn/article/id/15248
[2]	LIU L L, CHENG Y M, WANG X M, et al. System reliability analysis and risk assessment of a layered slope in spatially variable soils considering stratigraphic boundary uncertainty[J]. Computers and Geotechnics, 2017, 89: 213-225. doi: 10.1016/j.compgeo.2017.05.014
[3]	谭晓慧, 董小乐, 费锁柱, 等. 基于KL展开的可靠度分析方法及其应用[J]. 岩土工程学报, 2020, 42(5): 808-816. doi: 10.11779/CJGE202005002 TAN Xiaohui, DONG Xiaole, FEI Suozhu, et al. Reliability analysis method based on KL expansion and its application[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(5): 808-816. (in Chinese) doi: 10.11779/CJGE202005002
[4]	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
[5]	LIU Y, LI J F, SUN S Y, et al. Advances in Gaussian random field generation: a review[J]. Computational Geosciences, 2019, 23(5): 1011-1047. doi: 10.1007/s10596-019-09867-y
[6]	蒋水华. 水电工程边坡可靠度非侵入式随机分析方法[D]. 武汉: 武汉大学, 2014. JIANG Shuihua. A Non-Intrusive Stochastic Method for Slope Reliability in Hydroelectricity Engineering[D]. Wuhan: Wuhan University, 2014. (in Chinese)
[7]	GHANEM R G, SPANOS P D. Stochastic Finite Element: a Spectral Approach[M]. New York: Dover Publication, 2003.
[8]	PHOON K K, HUANG S P, QUEK S T. Implementation of Karhunen–Loeve expansion for simulation using a wavelet-Galerkin scheme[J]. Probabilistic Engineering Mechanics, 2002, 17(3): 293-303. doi: 10.1016/S0266-8920(02)00013-9
[9]	刘春凤. 应用数值分析[M]. 北京: 冶金工业出版社, 2005. LIU Chunfeng. Applied Numerical Analysis[M]. Beijing: Metallurgical Industry Press, 2005. (in Chinese)
[10]	秦权, 林道锦, 梅刚. 结构可靠度随机有限元: 理论及工程应用[M]. 北京: 清华大学出版社, 2006. QIN Quan, LIN Daojin, MEI Gang. Stochastic Finite Element Method for Structural Reliability: Theory and Engineering Application[M]. Beijing: Tsinghua University Press, 2006. (in Chinese)
[11]	SUDRET B, DER KIUREGHIAN A. Comparison of finite element reliability methods[J]. Probabilistic Engineering Mechanics, 2002, 17(4): 337-348. doi: 10.1016/S0266-8920(02)00031-0
[12]	仉文岗, 王琦, 陈福勇, 等. 考虑岩体空间变异性的边坡可靠度分析及抗滑桩随机响应研究[J]. 岩土力学, 2021, 42(11): 3157-3168. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202111024.htm ZHANG Wengang, WANG Qi, CHEN Fuyong, et al. Reliability analysis of slope and random response of anti-sliding pile considering spatial variability of rock mass properties[J]. Rock and Soil Mechanics, 2021, 42(11): 3157-3168. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202111024.htm
[13]	HUANG J, LYAMIN A V, GRIFFITHS D V, et al. Quantitative risk assessment of landslide by limit analysis and random fields[J]. Computers and Geotechnics, 2013, 53: 60-67. doi: 10.1016/j.compgeo.2013.04.009
[14]	PHOON K K, KULHAWY F H. Evaluation of geotechnical property variability[J]. Canadian Geotechnical Journal, 1999, 36(4): 625-639. doi: 10.1139/t99-039
[15]	FREDLUND D G, XING A Q. Equations for the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994, 31(4): 521-532. doi: 10.1139/t94-061
[16]	KIM J, HWANG W, KIM Y. Effects of hysteresis on hydro-mechanical behavior of unsaturated soil[J]. Engineering Geology, 2018, 245: 1-9. doi: 10.1016/j.enggeo.2018.08.004
[17]	AVERJANOV S F. About permeability of subsurface soils in case of incomplete saturation[J]. English Collection, 1950, 7: 19-21.
[18]	FREDLUND D G, MORGENSTERN N R, WIDGER R A. The shear strength of unsaturated soils[J]. Canadian Geotechnical Journal, 1978, 15(3): 313-321. doi: 10.1139/t78-029
[19]	陈朝晖, 黄凯华. 土质边坡可靠性分析的分层非平稳随机场模型[J]. 岩土工程学报, 2020, 42(7): 1247-1256. doi: 10.11779/CJGE202007008 CHEN Zhaohui, HUANG Kaihua. Non-homogeneous random field model for reliability analysis of slopes[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(7): 1247-1256. (in Chinese) doi: 10.11779/CJGE202007008
[20]	CHEN P, WEI C F. Numerical procedure for simulating the two-phase flow in unsaturated soils with hydraulic hysteresis[J]. International Journal of Geomechanics, 2016, 16(1): 04015030.
[21]	水利水电工程边坡设计规范: SL386—2007[S]. 北京: 中国水利水电出版社, 2007. Design Code for Engineered Slopes in Water Resources and Hydropower Projects: SL386—2007[S]. Beijing: China Water Power Press, 2007. (in Chinese)
[22]	李典庆, 肖特, 曹子君, 等. 基于高效随机有限元法的边坡风险评估[J]. 岩土力学, 2016, 37(7): 1994-2003. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201607021.htm LI Dianqing, XIAO Te, CAO Zijun, et al. Slope risk assessment using efficient random finite element method[J]. Rock and Soil Mechanics, 2016, 37(7): 1994-2003. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201607021.htm

施引文献(5)

期刊类型引用(3)

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