考虑参数旋转各向异性空间变异性的边坡大变形概率分析

Probabilistic Analysis of Large Slope Deformation Considering Soil Spatial Variability with Rotated Anisotropy

摘要

为揭示土体参数旋转各向异性空间变异性对边坡大变形特征的影响规律，将土体不排水抗剪强度参数模拟为旋转各向异性随机场，提出采用多重响应面法对随机场样本进行边坡稳定性安全系数高效求解和升序排列，进而使用随机物质点法按序模拟随机场样本的边坡大变形过程.以一饱和不排水粘土边坡为例，研究了旋转各向异性随机场的旋转角度β和弱主方向自相关长度θ ₂对边坡大变形特征和破坏模式的影响.结果表明：提出的基于多重响应面的随机物质点方法可以高效开展边坡大变形概率分析；β和θ ₂对边坡大变形特征和破坏模式均有显著影响；边坡大变形破坏的影响距离、滑动距离和滑动体积的平均值与标准差均随θ ₂的增加而增加；边坡大变形过程中可能产生4种破坏模式，其中深层滑动和渐进式滑动为该边坡模型的主要概率失稳模式.因此，提出的方法为边坡大变形概率分析提供了一条有效途径，考虑参数旋转各向异性空间变异性的边坡大变形概率分析对准确评估边坡失稳风险具有一定的理论参考意义.

关键词

边坡大变形 / 物质点法 / 随机场 / 旋转各向异性 / 空间变异性 / 多重响应面法 / 工程地质

Key words

large slope deformation / material point method / random field / rotated anisotropy / spatial variability / multiple response surface / engineering geology

0 引言

受构造运动、沉积过程、应力历史、物质组成和风化程度等因素的影响，天然土体物理力学性质（如抗剪强度）存在一定的空间变异性（张抒等， 2018； 何成等， 2021； 肖景红等， 2021）.在土体强度参数空间变异性影响下，边坡往往沿着坡体中土性分界面或阻力最小的路径发生失稳破坏（Hicks and Samy，2002； 吴益平等， 2020； Cheng et al.， 2021）.因此，科学定量表征土体参数的空间变异性，对有效揭示边坡稳定性和概率失稳机制具有重要意义.根据Zhu and Zhang （2013）研究总结，土体参数的空间变异性特征可被划分为6种主要形式：各向同性、横观各向异性、旋转各向异性、一般各向异性、一般旋转各向异性和复合各向异性.早期，考虑土体参数空间变异性的边坡可靠度分析以参数各向同性或横观各向异性研究为主（Wang et al.， 2011； Li et al.， 2016a， 2016b； Jiang et al.， 2017； Liu et al.， 2019）.近年来，其他4种各向异性结构对边坡稳定性的影响开始逐渐被学者们重视.例如，Cheng et al. （2018）提出了考虑参数旋转各向异性的随机有限差分方法，并将其运用于边坡稳定性概率分析和风险评价研究；Huang et al. （2019）使用随机极限平衡法研究了旋转各向异性对考虑条件随机场的边坡可靠性评价的影响；Zhu et al. （2019）基于随机极限平衡法分析了旋转各向异性特征对边坡稳定性和概率破坏模式的影响；Huang et al. （2021）将旋转各向异性与非平稳随机场相结合，使用随机极限平衡法研究了不排水条件下边坡的可靠度和破坏机理；Ng et al. （2021）使用随机极限平衡法研究了Copula函数的选择对考虑土体参数旋转各向异性的土坡稳定性的影响.然而，上述研究采用的均是传统的边坡稳定性分析方法，主要集中于边坡的破坏概率分析，而对边坡失稳后滑坡的过程、影响范围以及滑体对周边环境的影响程度等研究较少.究其原因，主要是因为极限平衡法只能计算得到边坡处于极限平衡状态时的潜在滑动面和相应的安全系数，而无法模拟滑坡的运动形态.同时，有限单元法/有限差分法也因为存在网格畸变等问题，难以模拟滑坡运动全过程.

然而，科学合理地量化边坡滑动后的滑体体积、滑动距离和影响范围等大变形参数，对开展边坡滑坡风险定量评价具有重要的理论意义（郭子正等， 2020）.目前，模拟滑坡大变形过程的方法主要包括光滑粒子流体动力学法（Wang et al.， 2019）、离散元法（Mao et al.， 2021； 李坤等， 2022）、物质点法（material point method， MPM）（Qu et al.， 2021； Ying et al.， 2021）等.其中，光滑粒子流体动力学法对流体力学模拟效果较好，而模拟固体力学问题时往往需要处理复杂的边界条件和克服拉伸不稳定性等问题，一定程度上限制了其在岩土工程中的应用（Xu and Yu， 2018； Wang et al.， 2021）.离散元法是一种非连续方法，在处理边界值问题时计算量较大，因而无法有效解决实际滑坡大变形问题（Augarde et al.， 2021）.而MPM作为一个应用范围广、与有限元有着类似的计算流程的大变形模拟方法（Sulsky et al.， 1994； 廉艳平等， 2013），在分析滑坡问题时具有一定优势.该方法同时使用欧拉网格和拉格朗日质点双重描述，在每个时间步重新划分计算网格，有效避免了网格畸变问题，可实现滑坡全过程模拟.因此，使用MPM对边坡滑坡大变形特征进行量化分析是滑坡灾害风险评估的一个有效途径.但是，相比目前较为成熟的边坡稳定性分析方法，MPM也存在计算效率低的缺点.特别是考虑土体参数空间变异性时，需要开展大量（如10⁴量级）的边坡大变形蒙特卡洛模拟，导致该问题变得更加突出.为此，Liu et al. （2019）提出了随机极限平衡法和随机物质点法（random material point method， RMPM）相结合的方法，以提高RMPM的计算效率.但当计算工况较多时，基于极限平衡法的蒙特卡洛模拟仍然十分耗时.另外，目前基于RMPM的边坡概率稳定性分析主要侧重于土体参数横观各向异性空间变异性对边坡概率失稳过程和破坏后果的影响，而鲜少考虑更符合实际的旋转各向异性空间变异性的影响.

鉴于上述分析，本文基于随机场理论、多重响应面法（multiple response surface method， MRSM） （Li et al.， 2016c）和MPM，提出一种考虑参数旋转各向异性空间变异性的高效边坡大变形概率分析方法——随机多重响应面-物质点法.以一个饱和不排水粘土边坡为例，采用随机场理论表征不排水抗剪强度的旋转各向异性空间变异性，在蒙特卡洛模拟框架内，利用MRSM高效判别边坡稳定性状态，进而采用RMPM模拟边坡滑坡运动全过程，以此定量分析土体参数旋转各向异性对边坡大变形的概率响应规律.

1 随机多重响应面-物质点法

1.1　旋转各向异性空间变异性随机场模拟

自然界中土体参数空间变异性类型十分复杂，图1给出了前述常见土体参数空间变异性类型中的4种（Zhu and Zhang， 2013）.其中，θ ₁和θ ₂所示方向为土体参数在空间上变化的大小主方向，其值越大，表示土体参数在该方向变化越平稳，反之则变化越剧烈.两个主方向的夹角为α，通常为90°.如图1a所示，各向同性表示在空间任意方向上土体参数的相关性都相同；横观各向异性则表示α为90°时水平方向的土体具有更强的物理参数自相关性，而垂直方向的土体参数自相关性则相对较弱，如图1b所示；图1c所示为旋转各向异性，这种形式是指在横观各向异性的基础上逆时针旋转一定的角度β，并保持自相关性最强和最弱的两个主方向正交的空间变异性形式；图1d所示为一般各向异性，指当一个主方向保持平行于坐标轴，另一个主方向与该主方向夹角为α时出现的空间变异性形式.目前，考虑土体参数空间变异性的边坡可靠度研究主要集中于图1a和图1b所示的各向同性和横观各向异性，而对于旋转各向异性和一般各向异性涉及较少，且图1d所示的自相关形式表征较为复杂.因此，本文主要考虑图1c所示的土体参数的旋转各向异性空间变异性.

在岩土工程可靠度分析中，土体参数空间变异性通常可以利用随机场理论（Cami et al.， 2020）来表征和模拟.根据随机场理论，空间变异性可由自相关长度和自相关函数表征：自相关长度用于表征土体参数在空间中的变化程度，自相关长度越大，表示参数的空间自相关程度越强，反之则越弱；自相关函数则用于量化空间中任意两点土体参数的自相关性.由于岩土工程中，可利用的参数数据有限，通常无法准确获取参数的空间自相关结构（即自相关函数和自相关长度），因此一般采用理论模型作为替代（Cami et al.， 2020）.Li et al.（2015）给出了岩土工程可靠度分析中常用的自相关函数类型和自相关长度范围，包括单指数型、高斯型、二阶马尔科夫型等.其中，高斯自相关函数在原点处连续可导，模拟的随机场更光滑，因此本文采用该函数表征土体参数的自相关性，其旋转各向异性形式如下：

式中：ρ表示土体参数在任意两点之间的自相关系数； \mathrm{\Delta } x^{\text{'}} = \mathrm{\Delta } x \mathrm{c} \mathrm{o} \mathrm{s} \beta + \mathrm{\Delta } y \mathrm{s} \mathrm{i} \mathrm{n} \beta， \mathrm{\Delta } y^{\text{'}} = - \mathrm{\Delta } x \mathrm{s} \mathrm{i} \mathrm{n} \beta + \mathrm{\Delta } y \mathrm{c} \mathrm{o} \mathrm{s} \beta， \mathrm{\Delta } x和 \mathrm{\Delta } y分别代表空间中任意两点在水平和垂直方向上的欧氏距离；θ ₁和θ ₂分别代表自相关性最强和最弱两主方向上的自相关长度； \beta为旋转角度，以逆时针方向为正，如图1c所示.图2a所示为基于上述高斯自相关函数生成的各向同性随机场自相关系数三维曲面.图2b和图2c所示则分别为基于高斯自相关函数生成的横观各向异性随机场旋转前后的自相关系数三维曲面，其中自相关性最强和最弱的两个主方向夹角为90°，旋转角度β为45°.据图2可观察到：各向同性随机场变量间的自相关系数在两个变量的水平和垂直空间距离均为 0 m的时候最大，且空间距离越大，自相关系数越小.横观各向异性随机场变量间的自相关系数在垂直空间距离为 0 m的时候最大，而旋转各向异性随机场的变量间的自相关系数在水平距离和垂直距离相等的时候最大.根据上述结果可知，受旋转角度控制，旋转各向异性随机场空间中两点的土体强度自相关性与横向各向异性随机场差别较大.

\rho = \mathrm{e} \mathrm{x} \mathrm{p} \left(-\frac{{\left|\mathrm{\Delta }{x}^{\text{'}}\right|}^{\mathrm{2}}}{{\theta }_{\mathrm{1}}^{\mathrm{2}}}-\frac{{\left|\mathrm{\Delta }{y}^{\text{'}}\right|}^{\mathrm{2}}}{{\theta }_{\mathrm{2}}^{\mathrm{2}}}\right) , (1)

目前，已提出多种随机场模拟方法，如局部平均细分法（Jiang et al.， 2022）、Karhunen-Loève级数展开法（Li et al.， 2019）和Cholesky分解法（Huang et al.， 2019； Liu et al.， 2019）.其中，Cholesky分解法由于概念简单、原理清晰，应用广泛，因此本文采用该方法进行土体参数旋转各向异性空间变异性随机场模拟.由于文献（Huang et al.， 2019； Liu et al.， 2019）对该方法的介绍已较多，本文对此不再赘述.

1.2　RMPM

MPM是Sulsky et al. （1994）提出的一种数值计算方法，该方法同时采用拉格朗日质点和欧拉网格对连续体进行描述，并将需要计算的连续体离散为一系列质点.离散后的每个质点代表连续体中的一部分区域，并携带了该区域物体的大部分信息，如质量、速度、应力和应变等.因此，所有质点的集合便可代表整个连续体.欧拉网格则仅用于质点信息映射到网格后动量方程的求解与空间导数的计算，不记录任何信息.在单个时间步内，MPM计算流程可归纳为以下3个步骤，如图3所示：

（1）映射阶段：将物质点的质量和动量映射到背景网格节点上，形成描述问题所需的外荷载、质量矩阵等.

（2）拉格朗日阶段：采用更新的拉格朗日有限元方法求解背景网格节点动量方程获得节点的速度与加速度，以此更新物质点的位置和速度.

（3）对流阶段：抛弃变形后的网格，使用初始阶段的规则的背景网格开始新时间步的计算.

MPM的应力更新方式包括USF （update stress first）、USL （update stress later）和MUSL （modified update stress later）三种.由于USF具有较好的能量守恒性（Nairn， 2013），本文采用USF进行应力更新.在每一个时间步内，使用USF进行应力更新的具体计算流程包括下列7个步骤：

（1） 将t时刻物质点p的质量 m_p^t和更新前t-1时刻的节点动量映射至背景网格节点i，获得t时刻的节点速度 {\mathit{v}}_i^t:

{\mathit{v}}_i^t = {\displaystyle \sum _{p=\mathrm{1}}^N} N_{ip} m_p^t {\mathit{v}}_p^{t-\mathrm{1}} / {\displaystyle \sum _{p=\mathrm{1}}^N} N_{ip} m_p^t , (2)

式中，N为物质点总数，N_ip 为节点i的形函数在物质点p处的值， {\mathit{v}}_p^{t-\mathrm{1}}为t-1时刻的物质点p的速度.

（2） 计算时刻为t+1物质点应变增量 \mathrm{\Delta } {\mathit{\epsilon }}_p和旋量增量 \mathrm{\Delta } {\mathit{\Omega }}_p，更新物质点的密度和应力：

\mathrm{\Delta } {\mathit{\epsilon }}_p = \left(\mathrm{\Delta }t/\mathrm{2}\right) {\displaystyle \sum _{i=\mathrm{1}}^M} \left[{\mathit{G}}_{ip}^t{\mathit{v}}_i^t+{\left({\mathit{G}}_{ip}^t{\mathit{v}}_i^t\right)}^{\mathrm{T}}\right] , (3)

\mathrm{\Delta } {\mathit{\Omega }}_p = \left(\mathrm{\Delta }t/\mathrm{2}\right) {\displaystyle \sum _{i=\mathrm{1}}^M} \left[{\mathit{G}}_{ip}^t{\mathit{v}}_i^t-{\left({\mathit{G}}_{ip}^t{\mathit{v}}_i^t\right)}^{\mathrm{T}}\right] , (4)

{\mathit{G}}_{ip}^t = {\left.\nabla N_{ip}\left(\mathit{x}\right)\right|}_{x={x^t}_p} , (5)

式中， \mathrm{\Delta } t \mathrm{为时} \mathrm{间步} \mathrm{长}，M为物质点ｐ影响域内的节点个数， {\mathit{G}}_{ip}^t \mathrm{为}形函数梯度.

（3） 基于本构模型更新物质点应力：

{\dot{\mathit{\sigma }}}_p = \mathit{T} : {\dot{\mathit{\epsilon }}}_p , (6)

式中， \mathit{T}为四阶张量， {\dot{\mathit{\sigma }}}_p \mathrm{与} {\dot{\mathit{\epsilon }}}_p分别为物质点p的应力率和应变率.

描述土体应力应变关系的本构模型包括Mohr-Coulomb模型（MC模型）、Drucker-Prager模型（DP模型）、剑桥模型等.为了避免出现计算不收敛的问题，本文采用非相关流动准则的DP模型模拟土体弹塑性行为.同时，土体大变形过程中通常会发生应变软化.为了模拟这一现象，本文参考了文献（Wang et al.， 2016）中的应变软化模型，得到不排水抗剪强度 S_u随土体有效塑性应变的变化规律如下：

(7)

S_u \left({\overline{\mathit{\epsilon }}}_p\right) = S_{ur} ; {\overline{\mathit{\epsilon }}}_p \ge {\overline{\mathit{\epsilon }}}_{pr} , (8)

式中， {\overline{\mathit{\epsilon }}}_p为有效塑性应变， S_{up}为土体峰值不排水抗剪强度， H为软化模量， {\overline{\mathit{\epsilon }}}_{pr}为 S_u折减到残余值 S_{ur}时的有效塑性应变阈值.

（4） 计算背景网格节点t+1时刻的内力 {\mathit{f}}_i^{int}和外力 {\mathit{f}}_i^{ext}：

{\mathit{f}}_i^{int} = - {\displaystyle \sum _{p=\mathrm{1}}^N} {\left.m_p{\mathit{\sigma }}^s\left({\mathit{x}}_p\right):\left(\nabla N_{ip}\right)\right|}_{x=x_p} , (9)

{\mathit{f}}_i^{ext} = {\displaystyle \sum _{p=\mathrm{1}}^N} m_p{\overline{\mathit{t}}}^s\left({\mathit{x}}_p\right)N_{ip}\left({\mathit{x}}_p\right)/h+ {\displaystyle \sum _{p=\mathrm{1}}^N} m_p\mathit{b}\left({\mathit{x}}_p\right)N_{ip}\left({\mathit{x}}_p\right) , (10)

式中， m_p为物质点p的质量， {\mathit{\sigma }}^{\mathit{s}}为Cauchy应力张量， {\mathit{x}}_p为物质点的位置， {\overline{\mathit{t}}}^{\mathit{s}}为边界域上的外力，h则是为了将左端最后一项边界积分转化为体积分而引入的假想边界层厚度， b 是单位体力.

（5） 求得t+1时刻节点加速度，则t+1时刻节点速度为：

{\mathit{v}}_i^{t+\mathrm{1}} = {\mathit{v}}_i^t + {\dot{\mathit{v}}}_i^{t+\mathrm{1}} \mathrm{\Delta } t . (11)

（6） 将节点信息映射回物质点，获得物质点的位置、速度以及加速度：

{\mathit{\Phi }}_p = {\displaystyle \sum _{i=\mathrm{1}}^M} N_{ip}{\mathit{\Phi }}_i . (12)

式中，M为物质点p影响域内的节点个数.

（7） 重置背景网格，下一个时间步从步骤（1）开始重新计算.

RMPM是在MPM的基础上，结合随机场理论描述岩土体参数的空间变异性，并在蒙特卡洛模拟方法框架内，用MPM模拟岩土体结构对岩土体参数空间变异性随机实现的响应的一种理论方法（Wang et al.， 2016）.在考虑土体参数空间变异性时，进行MPM应力更新计算需要将1.1节生成的土体强度参数随机场样本映射到岩土体计算模型中的每一个物质点上.如图4所示，本文参考文献（Liu et al.， 2019）中的做法，在进行RMPM计算时，使用与随机场单元尺寸一致的背景网格单元，因此背景网格单元内的物质点具有相同的土体强度参数.同时，由于本文使用的随机场单元为平面正方形单元，而MPM为了模拟平面应变，使用的背景网格单元为正方体单元，所以与边坡剖面垂直方向的物质点被赋予相同的土体抗剪强度参数.

1.3　本文方法实现步骤

基于蒙特卡洛模拟法的边坡大变形概率分析需要多次反复调用MPM计算程序，因而计算效率十分低下.另外，由于并非所有随机场样本均会导致边坡失稳产生大变形，因此减少不使边坡发生破坏的样本的物质点大变形模拟可以极大提高计算效率.MRSM作为一种简单高效的二次函数代理模型，已被广泛应用于边坡稳定性计算中.该方法通常采用中心复合设计法构建训练样本，将中心复合设计样本作为函数输入值，安全系数作为函数输出值，构建随机场变量和安全系数之间的显示函数关系（Li et al.， 2016c）.为此，本文基于随机场理论、MRSM和MPM，提出了考虑参数旋转各向异性空间变异性的随机多重响应面-物质点法.该方法首先基于随机场理论模拟生成一定数量的强度参数旋转各向异性随机场样本，随后采用MRSM高效预测边坡在所有随机场样本下的安全系数，并将样本按安全系数大小进行升序排列，最后使用RMPM计算程序对升序排列后的样本进行大变形计算分析，且在破坏概率P _f收敛后，终止计算.该方法计算流程如图5所示，可分为以下8个步骤：

（1） 确定旋转各向异性随机场与MPM的计算模型与计算参数.

（2） 利用GeoStudio和MATLAB进行MPM边坡建模并生成MPM计算文件.

（3） 建立土体强度参数旋转各向异性随机场并生成 N_{\mathrm{m}\mathrm{c}\mathrm{s}}组样本数据.

（4） 构建MRSM中心复合样本点矩阵，基于极限平衡法计算样本点对应的安全系数，并根据上述样本点和对应的安全系数计算MRSM系数矩阵，建立该边坡稳定性的MRSM替代模型.

（5） 将步骤3生成的 N_{\mathrm{m}\mathrm{c}\mathrm{s}}组样本数据代入MRSM中计算样本安全系数 FS，并将样本按 FS 大小升序排列，获得排序后的样本，分别标记为 s ¹， s ²， …， {\mathit{s}}^{N_{\mathrm{m}\mathrm{c}\mathrm{s}}}，用于下一步迭代计算.

（6） 对于第一次迭代，设置k=1，N _f=0，并对样本 {\mathit{s}}^k进行RMPM计算，如果失稳（即最大相对位移大于等于1 m（Liu et al.， 2019）），则令N _f=N _f+1，P _f=N _f /N；如果稳定（即最大相对位移小于1 m），则令N _f=N _f，P _f=N _f /N.

（7） 重复步骤（6），判断破坏概率P _f是否收敛（即500个样本内P _f不变化（Liu et al.， 2019））.如果P _f收敛，则终止计算；否则继续循环，且k=k+1，直至k= N_{\mathrm{m}\mathrm{c}\mathrm{s}}.

（8） 统计RMPM计算结果，并进行边坡破坏样本大变形特征和破坏模式概率分析.

为实现上述过程，本文同时使用MATLAB和基于Visual Studio平台的Intel Fortran软件作为开发工具.其中，旋转各向异性随机场模型和多重响应面模型使用MATLAB实现，RMPM则是在Zhang et al. （2016）编写的开源MPM3D代码（https://github.com/xzhang66/MPM3D-F90）的基础上使用Intel Fortran软件二次开发实现.数值模型计算使用的是硬件配置为Intel（R） Xeon（R） Gold 6226 CPU @ 2.70 GHz 2.69 GHz （2处理器）和128 GB RAM的计算机.

2 饱和不排水粘土边坡算例

为验证所提出的随机多重响应面-物质点法的有效性和利用该方法揭示不同旋转各向异性空间变异性特征对边坡大变形特征和破坏模式的影响规律的可行性，以文献中的一个均质饱和不排水粘土边坡为例进行说明（Liu et al.， 2019）.该边坡几何模型如图6所示，土体物理力学参数如表1所示.在进行随机场离散时，土体单元大小设置为 1 m×1 m，共得到1 816个随机场单元，一次典型随机场样本实现如图9所示.MPM背景网格六面体单元大小为1 m×1 m×1 m，物质点间距为 0.5 m，共得到14 432个物质点.MPM计算模型底部设置为固定边界条件，两侧边设置为对称边界条件.本文MPM主要计算参数如表1所示，计算时间步长为7.51×10^-4 s.采用分步加载重力（Andersen and Andersen， 2010）的方式生成初始应力场，加载时间为2 s.由于弹性模量E、泊松比 \upsilon和土体重度γ等参数对边坡大变形特征影响较小（Ma et al.， 2022），本文在进行概率分析时仅考虑了不排水抗剪强度S _u的空间变异性对边坡大变形特征的影响，并将表1中的弹性模量E、泊松比 \upsilon和土体重度γ等土体参数作为确定值赋予每一个物质点.

为定量评价边坡的破坏后果，本文采用影响距离、最大相对位移、滑动距离、滑动深度、滑动体积5个指标来描述边坡的大变形特征（Liu et al.， 2019），如图7所示.其中影响距离为从边坡破坏前的坡肩到破坏后的滑坡后缘的距离；最大相对位移为所有物质点滑动前后的相对位移的最大值；滑动距离为从边坡破坏前的坡脚到破坏后的滑坡前缘的距离；滑动深度为滑坡坡顶到滑面最低点的垂直距离；滑动体积为单位厚度所有滑动物质点（即相对位移大于1 m）的体积总和.另外，由于滑动距离为滑坡灾害风险评估中的一个重要指标，所以采用滑动距离超越概率（Sun et al.， 2021）对边坡的破坏后果进行进一步评价，该评价指标被定义为滑坡滑动距离超越某阈值的概率.

3 方法验证

3.1　边坡确定性分析

为验证MPM的准确性，本节对算例边坡进行确定性分析.在GeoStudio软件中利用Bishop法进行基于极限平衡法（LEM）的安全系数计算，得到的安全系数为1.35.基于物质点进行边坡确定性分析时使用的方法为物质点强度折减法（王双等， 2016）.该方法通过不断对土体强度参数进行折减，使边坡达到由稳定向失稳转变的临界状态，此时的折减系数即为安全系数.根据GeoStudio计算结果，本文采取折减系数初始值为1.30，折减系数增量为0.001进行折减.基于物质点强度折减法判断边坡处于临界状态的判据较多（蒋先平等， 2022），包括计算不收敛、特征点位移突变和塑性应变贯通等.本文以塑性应变贯通作为边坡处于临界状态的判据.如图8b所示，在折减系数为1.384时，边坡产生大变形，由此得到的安全系数为1.384.由上述结果可知，极限平衡法得到的安全系数与MPM接近，且得到的滑动面也基本一致，表明利用基于极限平衡法构建的MRSM进行随机物质点法破坏样本的高效预测具有一定的可行性.由图8b可知，用MPM计算得到的滑动面为圆弧形，这与粘性土滑坡通常为旋转滑动这一事实相符.需要说明的是，此时的滑动面深度接近边坡模型底边，边坡的基础底部深度可能会影响滑动面的扩展，但是由于本文的研究对象为边坡的影响距离、滑动距离、滑动体积以及破坏模式，其中只有破坏模式涉及到滑动深度，并且在4.4节中将滑动深度大于24 m且不为渐进式滑动的样本划分为深层滑动，所以算例边坡的基础底部深度对本文研究结论影响较小.

3.2　MRSM精度验证

本节首先以常见的S _u横观各向异性空间变异性情况为例，验证MRSM进行边坡稳定性分析的有效性.根据表3，将S _u模拟为均值100 kPa、变异系数为0.3的对数正态随机场，强主方向自相关长度θ ₁取24 m，弱主方向自相关长度θ₂ 取6 m.采用1.1节叙述的随机场模拟方法构建S _u随机场模型并随机生成500个随机场样本，图9为其中一次S _u随机场模型典型实现，图中颜色较深部分表示S _u值较大，颜色较浅部分表示S _u值较小.并且，水平方向变化较竖直方向更平滑，这是由于水平方向的自相关长度较垂直方向的自相关长度更大，说明了本文随机场模拟的有效性.本文MRSM的构建采用二次响应面函数，函数的系数矩阵通过中心复合样本点训练求得（Li et al.， 2016c）.在MRSM构建完成后，分别用随机极限平衡法和MRSM对上述500个样本进行安全系数计算，结果如图10所示.由图10可知：随机极限平衡法和MRSM得到的结果基本一致，且使用y=x的直线拟合上述结果得到决定系数R ²=1，表明MRSM对该工况下边坡安全系数的预测效果较好，能够代替随机极限平衡法进行随机场样本安全系数的高效计算.

3.3　随机多重响应面-物质点法准确性与计算效率

为验证随机多重响应面-物质点法的有效性，在3.2节随机场模型的基础上，生成40 000个样本进行MRSM计算，得到的破坏概率P _{f， MRSM}为6.33×10^-2.根据MRSM预测结果，将随机场样本进行升序排列，并用RMPM进行计算，其中随机极限平衡法计算单个样本需要约20 s，RMPM计算单个样本需要约0.8 h，结果如图11所示.由图11可知：P _{f， RMPM}在第4 071个样本收敛于5.08×10^-2，得到2 033个破坏样本，与MRSM在样本数为40 000时的计算结果（P _{f， MRSM}=6.33×10^-2）接近.由此可知，经过样本排序后的RMPM计算仅需计算少量样本 （4 071+500=4 571个）即能获得高精度的结果，并且在RMPM计算过程中节省了约（40 000- 4 571）×0.8-（2×1 816+1） ×20/3 600= 28 323 h，其中后者为构建MRSM需要计算中心复合设计样本安全系数的时间，因此极大地提高了计算效率.同时，图12绘制了破坏样本的滑动距离均值和标准差随排序后的计算样本数的变化曲线，由该图可知滑动距离的均值和标准差在第4 071个样本时也均收敛于定值（分别为 18.46 m和6.00 m），说明随机多重响应面-物质点法适用于边坡大变形特征的统计.需要说明的是，本文所提方法在一定程度上减少了样本计算量，为边坡风险定量评估提供了一条有效途径.

综合3.1、3.2和3.3节的验证结果可知，本文提出的随机多重响应面-物质点法能够在计算结果精度较高的情况下，极大减少直接利用RMPM计算的时间，为研究旋转各向异性空间变异性对边坡大变形特征和破坏模式的影响提供一条高效准确的途径.

4 边坡大变形概率分析结果

4.1　工况设置

为研究旋转各向异性空间变异性对边坡大变形特征和破坏模式的影响，仍然以前述饱和不排水粘土边坡为例，设置17个工况进行边坡大变形概率分析，各个工况随机场参数设置参考文献（Liu et al.， 2019），如表2所示.由于强主方向自相关长度θ ₁相比弱主方向自相关长度θ ₂对边坡的破坏概率影响更小，因此将θ ₂作为边坡大变形概率分析的主要影响参数.另外，为研究θ ₂在不同土体强度带状分布倾向下对边坡大变形特征和破坏模式的影响，分别在θ ₂=3 m、9 m和12 m时选取了旋转角度β=0°、45°和135°的情况作为典型工况进行概率分析.需要说明的是，由于本文研究的是边坡破坏后的大变形行为，所以进行分析统计的样本均为破坏样本.由于β=180°工况与β=0°的工况一致，对应的均是横观各向同性的情况，因此为了节约计算成本，下文中β=180°的工况的计算结果取β=0°的工况的计算结果.

4.2　破坏概率

图13为在不同工况下，基于随机多重响应面-物质点法计算得到的边坡破坏概率.由图13可知：P_f 随旋转角度的变化不呈现单调趋势，在0°<β<135°时，P _f较β=0°时小；135°≤β<180°时，P _f较β=0°时大，这与文献（程红战等， 2017）中的结论一致.β的大小反映了边坡土体强度分布倾向.当0°<β<90°时，边坡为逆层边坡，滑面剪切方向不易出现软弱面，因此破坏概率较β=0°的工况小.当90°<β<180°时，边坡为顺层边坡，且β接近135°时，滑面剪切方向较易出现软弱面，因此破坏概率相比较β=0°的工况大；而β远离135°时，滑面剪切方向不易出现软弱面，故破坏概率较β=0°的工况小.并且，随着θ ₂增大，P _f呈现整体增大的趋势.θ ₁与θ ₂的数值大小反映土体参数的空间自相关程度，值越大，空间自相关性越大，越容易出现贯通的软弱带，从而造成较大的破坏概率.当θ ₂=3 m，β=45°时，P _f接近于0，故不进行该工况下的滑坡失稳后果定量统计（图16）、滑动距离超越概率分析（图17）和破坏模式占比统计（图19）.

4.3　大变形特征

4.3.1　旋转角度的影响

图14为强主方向自相关长度θ ₁=24 m，弱主方向自相关长度θ ₂=6 m时，各工况下边坡大变形特征统计值随旋转角度β的变化曲线.由图14可知：影响距离均值和滑动体积均值随β的增加先减小后增大，并且在β=90°时达到最小值；而滑动距离均值随β的增加先减小后增大再减小，并且分别在β=45°和β=135°时达到最小值和最大值.滑动距离标准差随β增加先减小后增大，并且在β=90°时达到最小值；而影响距离标准差和滑动体积标准差随β的增加先减小后增大再减小，并且分别在β=45°和β=135°时达到最小值和最大值.β的变化反映了边坡土体强度带状分布倾向的变化，在β=0°时，由于土体强度带状分布倾向为水平方向，滑面剪切方向较β=45°和90°更易出现贯通的软弱带，因此该工况下的边坡大变形特征（即影响距离、滑动距离和滑动体积）均较大，并且较大的标准差表明了该工况下大变形特征的不确定性也相对其他工况较高；在0°<β<90°时，土体强度带状分布倾向与滑面剪切方向相反，此时边坡为逆层边坡，在滑面剪切方向不易穿透多个高强度带，因此难以形成整体滑动，导致该工况下样本整体的滑动距离较小；当β=90°时，土体强度带状分布倾向为垂直方向，边坡滑面剪切方向较难出现贯通的软弱带，滑体难以进行滑动且不易发生渐进式滑动造成二次破坏，因此该工况下的边坡大变形特征（即影响距离、滑动距离和滑动体积）的统计值均相对其他工况较小；当90°<β<180°时，土体强度带状分布倾向与滑面剪切方向相同，此时边坡为顺层边坡，导致边坡体内更容易出现贯通的软弱带，造成了滑动距离的增大.并且在该工况下，边坡不易发生渐进式滑动造成二次破坏，所以影响距离和滑动体积取决于滑面位置（即导致边坡滑动的土体强度软弱带位置）.上述结果表明：β的变化对边坡大变形特征有显著影响.β=0°导致边坡相对其他工况具有较大的大变形特征（即影响距离、滑动距离和滑动体积）和较高的不确定性，而在β=90°的工况下，边坡的大变形特征均值和不确定性均较小；β=135°导致边坡具有最大的滑动距离，而在β=45°的工况下，边坡的滑动距离最小.因此，实际工程应用中不考虑β的影响，将无法准确判断边坡的破坏后果.

图15为θ ₁=24 m，θ ₂=6 m时，不同β下边坡滑动距离超越概率曲线.在绘制该曲线时，将未发生破坏的样本的滑动距离取为0 m，故超越概率在0 m时的值近似等于各工况的破坏概率.由图15可知：90°<β≤180°时边坡的滑动距离超越概率整体大于0°<β≤90°的滑动距离超越概率；并且随着滑动距离的增大，各个工况的超越概率均逐渐降低.该结果表明：对边坡进行风险评估时，需考虑边坡土体旋转各向异性随机场的旋转角度β对边坡滑动距离超越概率的影响.取某一超越概率为设计值（即超越概率为确定值）时，顺层边坡的滑动距离更大，所以其影响范围也更大；而对边坡滑动范围内的某个建筑物（即进行超越概率计算的滑动距离阈值为确定值）进行风险性评估时，顺层边坡的滑动距离超越概率更大，从而对该建筑物造成影响的概率也更大.

4.3.2　弱主方向自相关长度的影响

图16为强主方向自相关长度θ ₁=24 m，旋转角度β=0°、45°和135°时，各工况下边坡大变形特征统计值随弱主方向自相关长度θ₂ 的变化情况.由图16可知：β=0°、45°和135°时，影响距离、滑动距离及滑动体积的均值和标准差均随θ₂ 的增大而增大.这是因为θ₂ 的增加会导致软弱带或者高强度带的产生，而这种情况下容易产生极端情况，即边坡大变形程度较大或较小的情况产生，从而增加土体大变形特征（影响距离、滑动距离和滑动体积）的不确定性，这与文献（Wang et al.， 2016； Ma et al.， 2022）的结果一致.上述结果表明，θ ₂越大（即土体参数空间自相关程度越大）时，边坡大变形程度越大，不确定性更强，因此在实际工程设计时应关注土体参数空间自相关程度对边坡风险的影响.

图17为θ ₁=24 m，β=0°、45°和135°时，不同θ ₂下滑动距离超越概率曲线.由图17可知：β=0°、45°和135°时，随着θ ₂的增大，边坡在不同旋转角度β下的滑动距离超越概率均呈现整体增大的趋势.上述结果表明，对边坡进行风险评估时，也需考虑土体参数空间自相关程度对边坡滑动距离超越概率的影响.在取某一超越概率作为设计值时，θ ₂越大（即土体空间自相关程度越强），边坡滑动距离更大，所以其影响范围也更大；而对边坡滑动范围内的某个建筑物（即滑动距离阈值为确定值）进行风险性评估时，土体参数空间自相关程度较强的边坡的滑动距离超越概率更大，从而对该建筑物造成影响的概率也更大.

4.4　破坏模式

土体强度空间分布的差异会对边坡造成不同的破坏模式和大变形特征.本文参考文献（Liu et al.， 2019），根据滑动深度和滑动面数划分了4种边坡破坏模式，分别为浅层滑动（即滑动深度在（0 m，18 m］内，且滑面数为1）、中层滑动（即滑动深度在（18 m，24 m］内，且滑面数为1）、深层滑动（即滑动深度在（24 m，28 m］内，且滑面数为1）和渐进式滑动（即滑面数大于1）.其中，每个破坏样本的滑动面数，可根据每个破坏样本的边坡整体动能随时间变化的波峰数进行判断.

4.4.1　旋转角度的影响

图18为强主方向自相关长度θ ₁=24 m，弱主方向自相关长度θ ₂=6 m时，不同旋转角度β下边坡4种破坏模式占各工况总破坏样本数的比例.由图18可知：深层滑动和渐进式滑动为该边坡的主要破坏模式，且深层滑动占比最大；β=45°时，深层滑动占比较其他工况小，而渐进式滑动占比较大，表明β=45°时边坡发生的破坏为渐进式滑动的可能性较大，这是因为该工况下破坏样本的土体强度软弱带与坡面近似垂直，滑坡发生后易发生二次滑动，从而促进了渐进式滑动的发生；β=135°时，渐进式滑动占比较其他工况小，而深层滑动占比较大，这是因为该工况下破坏样本的土体强度软弱带与坡面近似平行，虽有利于滑坡的发生，但不易发生二次滑动.

4.4.2　弱主方向自相关长度的影响

图19为强主方向自相关长度θ ₁=24 m，旋转角度β=0°、45°和135°时，边坡4种破坏模式占各工况总破坏样本数的比例随弱主方向自相关长度θ ₂的变化规律.由图19可知：β=0°时，渐进式滑动的占比随θ ₂的增大而增大，浅层滑动和中层滑动占比较小且中层滑动占比随θ ₂的增大而减小；β=45°时，θ ₂的增大使深层滑动占比显著增加，使渐进式滑动占比显著下降，且浅层滑动和中层滑动占比接近于0；β=135°时，主要破坏模式为深层滑动，其他破坏模式占比较小，这可能与该工况下的边坡为顺层边坡，土体强度参数沿坡面分布有关.上述变化情况表明，该边坡可能发生的破坏模式不仅与旋转角度有关，而且还受土体强度参数空间自相关程度的影响.

4.4.3　不同破坏模式下的大变形特征概率分析

图20为强主方向自相关长度θ ₁=24 m，弱主方向自相关长度θ ₂=6 m和旋转角度β=0°时，不同破坏模式下边坡大变形特征分布直方图.由图20可知：深层滑动和渐进式滑动为该工况的主要破坏形式，浅层滑动模式占比为0；中层滑动的影响距离主要分布在［13 m， 17 m］，深层滑动的影响距离主要分布在［17 m， 25 m］和40 m；渐进式滑动的影响距离在［13 m， 39 m］分布相对均匀，而在40 m处分布集中.这是因为渐进式滑动的滑面有向远离坡面方向发展的趋势且可能产生多条滑面，从而使影响距离增大.且边坡模型边界条件也会阻碍滑坡的发展，从而导致滑动模式为深层滑动和渐进式滑动的样本的影响距离在40 m处的频数较大；中层滑动的滑动距离主要分布在［17 m， 22 m］，深层滑动的滑动距离主要分布在［11 m， 20 m］，而渐进式滑动的滑动距离分布相对均匀，主要分布在［20 m，30 m］，这是因为渐进式滑动为多次破坏，导致滑动距离增大；中层滑动的滑动体积主要分布在600 m²左右，深层滑动的滑动体积主要分布在［800 m²， 1 000 m²］，渐进式滑动的滑动体积主要分布在［600 m²， 1 400 m²］，分布较均匀.上述结果表明不同破坏模式具有不同边坡大变形特征分布，因此在实际工程设计时应考虑不同破坏模式发生的概率和产生的风险，对滑坡可能造成的危害做出准确评估，以确保设计更安全.

表3为θ ₁=24 m、θ ₂=6 m和β=0°时，不同破坏模式下边坡大变形特征统计值.其中，浅层滑动占比为0，故不进行边坡大变形特征统计.由表3可知，不同破坏模式的边坡大变形特征统计值具有一定差异：深层滑动模式相比中层滑动模式具有更大的影响距离均值与标准差、滑动体积均值与标准差和滑动距离标准差以及更小的滑动距离均值；渐进式滑动的影响距离、滑动距离和滑动体积均值和标准差相比其他破坏模式均最大.上述结果表明在进行边坡风险评估时，应当特别关注渐进式滑动占比较大的情形并采取有效的防治措施，以防造成较大的经济损失和人员伤亡.

5 结论

（1） 将MRSM和RMPM相结合，提出了一种考虑土体强度参数旋转各向异性空间变异性的边坡大变形概率分析方法——随机多重响应面-物质点法.该方法通过使用MRSM提高样本安全系数计算效率和减小RMPM的计算量，达到提高计算效率的目的，为边坡大变形概率分析提供了一条高效准确的途径.

（2） 旋转角度β和弱主方向自相关长度θ ₂对边坡大变形特征有一定影响.当β=0°时，边坡大变形特征（即影响距离、滑动距离和滑动体积）均较大，不确定性也相对较高，而β=90°时，边坡大变形特征和不确定性均较小；当β=45°和135°时，边坡的滑动距离均值分别具有最小值和最大值，并且滑动距离超越概率在β=0°和β=135°时较β=45°和β=90°时大；边坡大变形特征统计值均随θ ₂的增加而增加，并且随着θ ₂的增加，不同β下的滑动距离超越概率也均增大.

（3） β和θ ₂对边坡破坏模式占比有较大影响，且深层滑动和渐进式滑动为该边坡的主要概率失稳模式.β=45°时，深层滑动占比较其他旋转角度工况下显著减小，渐进式滑动占比显著提高，且在该旋转角度下，深层滑动和渐进式滑动占比随θ ₂的增加有明显的增加和减少；β=135°时，渐进式滑动相比其他旋转角度工况占比最小，深层滑动占比最大，为该工况下的主要破坏模式.

（4） 不同破坏模式反映了边坡不同的大变形特征.深层滑动模式相比中层滑动模式具有更大的影响距离均值与标准差、滑动体积均值与标准差和滑动距离标准差以及更小的滑动距离均值，而渐进式滑动的影响距离、滑动距离和滑动体积均值和标准差相比其他破坏模式均最大.

参考文献

[1] Andersen, S., Andersen, L., 2010. Modelling of Landslides with the Material-Point Method. Computational Geosciences, 14: 137-147. https://doi.org/10.1007/s10596-009-9137-y

[2] Augarde, C. E., Lee, S. J., Loukidis, D., 2021. Numerical Modelling of Large Deformation Problems in Geotechnical Engineering: A State-of-the-Art Review. Soils and Foundations, 61(6): 1718-1735. https://doi.org/10.1016/j.sandf.2021.08.007

[3] Cami, B., Javankhoshdel, S., Phoon, K. K., et al., 2020. Scale of Fluctuation for Spatially Varying Soils: Estimation Methods and Values. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering, 6(4): 03120002. https://doi.org/10.1061/AJRUA6.0001083

[4] Cheng, H.Z., Chen, J., Chen, R. P., et al., 2018. Risk Assessment of Slope Failure Considering the Variability in Soil Properties. Computers and Geotechnics, 103: 61-72. https://doi.org/10.1016/j.compgeo.2018.07.006

[5] Cheng, H. Z., Chen, J., Wang, Z. S., et al., 2017. Stability Analysis of a Clay Slope Accounting for the Rotated Anisotropy Correlation Structure. Chinese Journal of Rock Mechanics and Engineering, 36(S2): 3965-3973 (in Chinese with English abstract).

[6] Cheng, H.L., Zhou, J.M., Chen, Z.Y., et al., 2021. A Comparative Study of the Seismic Performances and Failure Mechanisms of Slopes Using Dynamic Centrifuge Modeling. Journal of Earth Science, 32(5): 1166–1173. https://doi.org/10.1007/s12583-021-1481-4

[7] Guo, Z.Z., Yin, K.L., Liu, Q.L., et al., 2020. Rainfall Warning of Creeping Landslide in Yunyang County of Three Gorges Reservoir Region Based on Displacement Ratio Model. Earth Science, 45(2): 672-684 (in Chinese with English abstract).

[8] He, C., Tang, H. M., Shen, P. W., et al., 2021. Progressive Failure Mode and Stability Reliability of Strain- Softening Slope. Earth Science, 46(2): 697-707 (in Chinese with English abstract).

[9] Hicks, M. A., Samy, K., 2002. Influence of Heterogeneity on Undrained Clay Slope Stability. Quarterly Journal of Engineering Geology and Hydrogeology, 35(1): 41-49.

[10] Huang, L., Cheng, Y. M., Leung, Y.F., et al., 2019. Influence of Rotated Anisotropy on Slope Reliability Evaluation Using Conditional Random Field. Computers and Geotechnics, 115: 103133. https://doi.org/10.1016/j.compgeo.2019.103133

[11] Huang, L., Cheng, Y. M., Li, L., et al., 2021. Reliability and Failure Mechanism of a Slope with Non-Stationarity and Rotated Transverse Anisotropy in Undrained Soil Strength. Computers and Geotechnics, 132: 103970. https://doi.org/10.1016/j.compgeo.2020.103970

[12] Jiang, S. H., Huang, J., Griffiths, D. V., et al., 2022. Advances in Reliability and Risk Analyses of Slopes in Spatially Variable Soils: A State-of-the-Art Review. Computers and Geotechnics, 141: 104498. https://doi.org/10.1016/j.compgeo.2021.104498

[13] Jiang, S. H., Huang, J., Yao, C., et al., 2017. Quantitative Risk Assessment of Slope Failure in 2-D Spatially Variable Soils by Limit Equilibrium Method. Applied Mathematical Modelling, 47: 710-725. https://doi.org/10.1016/j.apm.2017.03.048

[14] Jiang, X. P., Zhang, P., Lu, Y. W., et al., 2022. Slope Failure Criterion for the Strength Reduction Material Point Method. Bulletin of Geological Science and Technology, 41(2): 113-122 (in Chinese with English abstract).

[15] Li, D. Q., Jiang, S. H., Cao, Z. J., et al., 2015. A Multiple Response-Surface Method for Slope Reliability Analysis Considering Spatial Variability of Soil Properties. Engineering Geology, 187: 60-72. https://doi.org/10.1016/j.enggeo.2014.12.003

[16] Li, D. Q., Xiao, T., Cao, Z. J., et al., 2016a. Enhancement of Random Finite Element Method in Reliability Analysis and Risk Assessment of Soil Slopes Using Subset Simulation. Landslides, 13: 293-303. https://doi.org/10.1007/s10346-015-0569-2

[17] Li, D. Q., Xiao, T., Cao, Z. J., et al., 2016b. Efficient and Consistent Reliability Analysis of Soil Slope Stability Using Both Limit Equilibrium Analysis and Finite Element Analysis. Applied Mathematical Modelling, 40(9-10): 5216-5229. https://doi.org/10.1016/j.apm.2015.11.044

[18] Li, D. Q., Zheng, D., Cao, Z. J., et al., 2016c. Response Surface Methods for Slope Reliability Analysis: Review and Comparison. Engineering Geology, 203: 3-14. https://doi.org/10.1016/j.enggeo.2015.09.003

[19] Li, D. Q., Zheng, D., Cao, Z. J., et al., 2019. Two-Stage Dimension Reduction Method for Meta-Model Based Slope Reliability Analysis in Spatially Variable Soils. Structural Safety, 81: 101872. https://doi.org/10.1016/j.strusafe.2019.101872

[20] Li, K., Cheng, Q. G., Lin, Q. W., et al., 2022. State of the Art on Rock Avalanche Dynamics from Granular Flow Mechanics. Earth Science, 47(3): 893-912 (in Chinese with English abstract).

[21] Lian, Y. P., Zhang, F., Liu, Y., et al., 2013. Material Point Method and Its Applications. Advances in Mechanics, 43(2): 237-264 (in Chinese with English abstract).

[22] Liu, X., Wang, Y., Li, D., 2019. Investigation of Slope Failure Mode Evolution during Large Deformation in Spatially Variable Soils by Random Limit Equilibrium and Material Point Methods. Computers and Geotechnics, 111: 301-312. https://doi.org/10.1016/j.compgeo.2019.03.022

[23] Ma, G. T., Rezania, M., Nezhad, M. M., 2022. Uncertainty Quantification of Landslide Runout Motion Considering Soil Interdependent Anisotropy and Fabric Orientation. Landslides, 19 (2022): 1231-1247. https://doi.org/10.1007/s10346-021-01795-2

[24] Mao, J., Liu, X., Zhang, C., et al., 2021. Runout Prediction and Deposit Characteristics Investigation by the Distance Potential-Based Siscrete Element Method: The 2018 Baige Landslides, Jinsha River, China. Landslides, 18: 235-249. https://doi.org/10.1007/s10346-020-01501-8

[25] Nairn, J. A., 2013. Modeling Imperfect Interfaces in the Material Point Method Using Multimaterial Methods. CMES: Computer Modeling in Engineering & Sciences, 92(3): 271-299.

[26] Ng, C. W., Qu, C., Cheung, R. W., et al., 2021. Risk Assessment of Soil Slope Failure Considering Copula-Based Rotated Anisotropy Random Fields. Computers and Geotechnics, 136: 104252. https://doi.org/10.1016/j.compgeo.2021.104252

[27] Qu, C. X., Wang, G., Feng, K., et al., 2021. Large Deformation Analysis of Slope Failure Using Material Point Method with Cross-Correlated Random Fields (Cover Paper). Journal of Zhejiang University - Science A: Applied Physics & Engineering, 22: 856-869. https://doi.org/10.1631/jzus.A2100196

[28] Sulsky, D., Chen, Z., Schreyer, H. L., 1994. A Particle Method for History-Dependent Materials. Computer Methods in Applied Mechanics and Engineering, 118(1-2): 179-196.

[29] Sun, X., Zeng, P., Li, T., et al., 2021. From Probabilistic Back Analyses to Probabilistic Run-out Predictions of Landslides: A Case Study of Heifangtai Terrace, Gansu Province, China. Engineering Geology, 280: 105950. https://doi.org/10.1016/j.enggeo.2020.105950

[30] Wang, B., Hicks, M. A., Vardon, P. J., 2016. Slope Failure Analysis Using the Random Material Point Method. Geotechnique Letters, 6(2): 113-118. https://doi.org/10.1680/jgele.16.00019

[31] Wang, L., Xu, F., Yang, Y., 2021. An Improved Total Lagrangian SPH Method for Modeling Solid Deformation and Damage. Engineering Analysis with Boundary Elements, 133: 286-302. https://doi.org/10.1016/j.enganabound.2021.09.010

[32] Wang, S., Li, X.C., Shi, L., et al., 2016.Material Point Strength Reduction Method and Its Application to Slope Engineering. Rock and Soil Mechanics, 37(9):2672-2678 (in Chinese with English abstract).

[33] Wang, Y., Cao, Z., Au, S. K., 2011. Practical Reliability Analysis of Slope Stability by Advanced Monte Carlo Simulations in a Spreadsheet. Canadian Geotechnical Journal, 48(1): 162-172. https://doi.org/10.1007/978-3-662-52914-0_7

[34] Wang, Y., Qin, Z. W., Liu, X., et al., 2019. Probabilistic Analysis of Post-Failure Behavior of Soil Slopes Using Random Smoothed Particle Hydrodynamics. Engineering Geology, 261: 105266. https://doi.org/10.1016/j.enggeo.2019.105266

[35] Wu, Y. P., Lu, L. E., Xue, Y., 2020. Application of Landslide Progressive Failure Mechanical Model Based on the Critical Stress State. Bulletin of Geological Science and Technology, 39(5): 1-7 (in Chinese with English abstract).

[36] Xiao, J H., Wang, M., Wang, C., et al., 2021. Reliability Analysis of Slope with Dominant Seepage Interlayer under Rainfall Infiltration. Bulletin of Geological Science and Technology, 40(6): 193-204 (in Chinese with English abstract).

[37] Xu, X., Yu, P., 2018. A Technique to Remove the Tensile Instability in Weakly Compressible SPH. Computational Mechanics, 62: 963-990. https://doi.org/10.1007/s00466-018-1542-4

[38] Ying, C., Zhang, K., Wang, Z. N., et al., 2021. Analysis of the Run-out Processes of the Xinlu Village Landslide Using the Generalized Interpolation Material Point Method. Landslides, 18 (4): 1519-1529. https://doi.org/10.1007/s10346-020-01581-6

[39] Zhang, S., Tang, H. M., Liu, X., et al., 2018. Seepage and Instability Characteristics of Slope Based on Spatial Variation Structure of Saturated Hydraulic Conductivity. Earth Science, 43(2): 622-634 (in Chinese with English abstract).

[40] Zhang, X., Chen, Z., Liu, Y., 2016. The Material Point Method: A Continuum-Based Particle Method for Extreme Loading Cases. Elsevier Inc., London, 165-190.

[41] Zhu, H., Zhang, L. M., 2013. Characterizing Geotechnical Anisotropic Spatial Variations Using Random Field Theory. Canadian Geotechnical Journal, 50 (7): 723-734. https://doi.org/10.1139/cgj-2012-0345

[42] Zhu, H., Zhang, L. M., Xiao, T., 2019. Evaluating Stability of Anisotropically Deposited Soil Slopes. Canadian Geotechnical Journal, 56 (5): 753-760. https://doi.org/10.1139/cgj-2018-0210

[43] 程红战， 陈健， 王占盛， 等， 2017. 考虑旋转各向异性相关结构的黏土边坡稳定性分析. 岩石力学与工程学报， 36（S2）： 3965-3973.

[44] 郭子正， 殷坤龙， 刘庆丽， 等， 2020. 基于位移比模型的三峡库区云阳县域内蠕变型滑坡降雨预警. 地球科学， 45（2）： 672-684.

[45] 何成， 唐辉明， 申培武， 等， 2021. 应变软化边坡渐进破坏模式及稳定性可靠度. 地球科学， 46（2）： 697-707.

[46] 蒋先平， 张鹏， 卢艺伟， 等， 2022. 物质点强度折减法边坡失稳判据选择方法. 地质科技通报， 41（2）： 113-122.

[47] 李坤， 程谦恭， 林棋文， 等， 2022. 高速远程滑坡颗粒流研究进展. 地球科学， 47（3）： 893-912.

[48] 廉艳平， 张帆， 刘岩， 等， 2013. 物质点法的理论和应用. 力学进展， 43（2）： 237-264.

[49] 王双， 李小春， 石露， 等， 2016.物质点强度折减法及其在边坡中的应用. 岩土力学， 37（9）： 2672-2678.

[50] 吴益平， 卢里尔， 薛阳， 2020. 基于临界状态的边坡渐进破坏力学模型分析及应用. 地质科技通报， 39（5）： 1-7.

[51] 肖景红， 王敏， 王川， 等， 2021. 含优势渗流层边坡降雨入渗下的可靠度分析. 地质科技通报， 40（6）： 193-204.

[52] 张抒， 唐辉明， 刘晓， 等， 2018. 基于饱和渗透系数空间变异结构的斜坡渗流及失稳特征. 地球科学， 43（2）： 622-634.

基金资助

国家自然科学基金项目(41902291)

湖南省自然科学基金项目(2022JJ20058;2020JJ5704;2020JJ5015)

中南大学中央高校基本科研业务费专项资金项目(2022ZZTS0501)