基于地质结构探测的多滑面边坡系统可靠度分析

System Reliability Analysis of Multi-Slip Surface Slopes Based on Geological Structure Detection

摘要

地质结构识别与强度参数不确定性量化是岩质边坡稳定性评估的核心问题.为此，提出了一种基于地质结构探测的多滑面边坡系统可靠度分析方法.该方法首先结合多道勒夫波分析（multichannel analysis of Love waves，MALW）与初至旅行时层析成像（first-arrival travel time tomography，FATT），实现软弱层与断层探测效果的互补.随后，通过弹性波速折减软弱层强度参数，统计获取其概率分布.最后，考虑参数不确定性，计算边坡地表位移及单一滑面与系统失效概率.该方法在两个边坡案例的测试中表明：多阶频散曲线对边坡深部及浅部软弱层的反演精度均优于基阶频散曲线，且勒夫波较瑞雷波受岩层界面起伏的影响更小.边坡断层在初至旅行时记录中表现为特定范围内的波动特征，基于该特征的反演可定位局部断层.对于多滑面边坡案例，结合地表位移，地质结构与各滑面单一失效概率，确定边坡主要受深层滑面控制.且该边坡系统失效概率受内摩擦角变异系数的影响远大于黏聚力.此方法能够有效探测边坡地质结构，定位软弱层与断层.可考虑软弱层内部裂隙结构和风化程度的影响，折减强度参数并量化其不确定性.能够准确识别关键控制滑动面，定量评估各潜在滑动面及系统失效概率，为边坡防治提供了科学依据和参考.

Abstract

The identification of geological structures and the quantification of uncertainties in strength parameters are crucial for assessing the stability of rock slopes. This study proposes a system reliability analysis method for multi-slip surface slopes based on geological structure detection. The method integrates multichannel analysis of Love waves (MALW) and first-arrival travel time tomography (FATT) to achieve complementary detection of weak layers and faults. Elastic wave velocities are used to reduce the strength parameters of weak layers, and their probabilistic distributions are statistically derived. By incorporating parameter uncertainties, the surface displacement of slopes and the failure probabilities of individual slip surfaces and the entire system are calculated. Case studies indicate that multi-mode dispersion curves achieve higher inversion accuracy for both deep and shallow weak layers compared to fundamental-mode dispersion curves, and Love waves are less affected by undulations at rock layer interfaces than Rayleigh waves. Slope faults exhibit characteristic wave fluctuations within specific ranges in first-arrival travel time records. Inversion based on these features enables the localization of local faults.In the case of multi-slip surface slopes, deep slip surfaces are identified as the primary controlling factor, and the system failure probability is more significantly affected by the coefficient of variation of the internal friction angle than by cohesion. This method effectively detects geological structures, locates weak layers and faults, and quantifies uncertainties in strength parameters, providing a scientific basis for slope stability analysis and mitigation measures.

Graphical abstract

关键词

岩质边坡 / 参数不确定性 / 系统可靠度分析 / 地质结构探测 / 工程地质学.

Key words

rock slope / parameter uncertainty / system reliability analysis / geological structure detection / engineering geology

0 引言

岩质滑坡灾害每年都在全国范围内造成大量的人员伤亡和财产损失.例如2014年贵州省黔南州福泉市道坪镇英坪村发生大规模山体滑坡，造成大量房屋被掩埋，23人死亡.2017年四川省茂县叠溪镇新磨村高位岩质滑坡，掩埋了64户农房和1 500 m道路，导致83人死亡或失踪（李为乐等， 2019）.边坡地质结构及强度参数对准确评估边坡稳定性，预防滑坡灾害至关重要.然而，由于地下结构的隐蔽性和复杂性，边坡结构面定位困难，强度参数不确定性难以量化，滑动面亦难以识别.

边坡结构面，如软弱层和断层，可通过侵入式和非侵入式两类方法定位.侵入式方法，是一种直观揭示边坡地下结构的方法.钻孔可以收集不同深度土壤或岩石样本，从而判断分析地下的组成、性质和地质特征.然而，在工程中钻孔的数量往往受到经济因素和现场技术的限制.孔间地层通常通过钻孔处地层线性连接获取，这不可避免地降低了地层划分的准确度，为潜在滑动面识别带来较大不确定性（Wang et al.， 2018）.地震勘探作为一种广泛应用的非侵入式方法，包括多道面波分析（MASW）、多道勒夫波分析（MALW）和初至旅行时层析成像（FATT）等.由于其低成本、高效率，目前已被广泛应用于地下调查和潜在滑动面识别中.这类方法通过地震信号获取边坡速度结构，识别速度显著降低区域，定位软弱层及断层.多道面波分析（MASW）依赖于瑞雷波在垂直异向地层中传播时的频散特性反演速度剖面（Park et al.， 1999）.其频散曲线，包括基阶和高阶曲线，可通过τ-p变换（Dunne and Beresford， 1995）、相移法（Park et al.， 1998）等从瑞雷波信号中提取.多阶频散曲线反演由于其更大的有效频率范围，以及高阶频散曲线对深层参数更为敏感这一特征（Luo et al.， 2007），可提高反演的精度和范围.与MASW相似，多道勒夫波分析方法（MALW）依赖于勒夫波的频散特性（Xia et al.， 2012），与瑞雷波相比，勒夫波是一种纯SH波，能够避免膨胀波的干扰（Xie and Liu， 2015），频散能量更加集中，“模式接吻”现象较少（Xia，2014），更易于提取高阶频散曲线.然而，对于MASW和MALW，两者的理论频散曲线均基于层状速度结构计算（Knopoff， 1964）.因此，当遇到断层等横向非均匀介质时，反演的分辨率会降低，甚至出现反演失效（Pan et al.， 2018）.相比之下，初至旅行时层析成像（FATT）则具有较高的横向分辨率（Wang et al.， 2022）.FATT的原理是搜索一个最佳速度模型，使其计算初至旅行时与观测初至旅行时相匹配.这是一个高度非线性的反演问题，近年来，一些学者提出了几种有效解决方法，如并行竞争粒子群优化（Luu et al.， 2018）、监督下降学习技术（Guo et al.， 2019）和伴随状态法（Taillandier et al.， 2009）.其中，伴随状态法可以避免繁琐的射线追踪技术，适用于处理大规模地震数据.然而，FATT在估算界面深度，即层状结构成像方面分辨率不足（Liu et al.， 2010）.因此，单一使用MALW或FATT对复杂结构边坡难以全面准确定位软弱层和断层位置.

强度参数不确定性表征及滑动面识别是边坡系统可靠度分析中的两个核心问题.强度参数的不确定性通常以概率分布形式表示，其参数特性可通过对现场勘验数据统计估计（杨智勇等， 2022）或基于经验准则确定（刘国锋等， 2023）.两类方法均在不同程度上依赖于钻孔岩样的室内试验强度参数.但受制于勘探成本、场地条件及设备限制，钻孔试验数据通常较为有限，难以全面、准确地估计岩土参数的统计特征（蒋水华等， 2020）.此外，室内试验岩样为完整岩块，不能反映岩体内部的节理、裂隙等对强度参数的削弱作用，导致参数估计结果与实际工程条件存在偏差，这一问题在风化程度较高的岩体中尤为突出.波速与岩体的风化程度密切相关（王亮清等， 2010），因此，基于地震勘探方法获取的岩石波速可为强度参数的折减与不确定性标定提供一种可行的途径.针对具有多个潜在滑动面的边坡，系统可靠度分析可通过比对各滑动面的失效概率及相关系数来识别代表性滑动面（李静萍等， 2016）.在进行系统可靠度分析前，也可基于监测数据对滑动面进行定位.一些高密度监测手段，如钻孔线式位移计，倾斜测井套管及多场一体化观测等，不仅可以推测边坡的滑动特性（Dou et al.， 2025），确定滑动面的深度（Crosta et al.， 2014），还能够进行实现滑坡预测（窦杰等， 2023）.然而，这类密集监测系统需要大量的资金和人力投入.相比之下，地表位移数据易于获取，可通过设定位移阈值来预测滑坡并可有效区分滑体与滑床.因此，地质勘探获取的边坡地质结构与地表位移相结合，可为滑动面定位提供一种可行方案.

为此，本文提出了基于地质结构探测的多滑面边坡系统可靠度分析方法，包括3个模块：多道勒夫波分析（MALW）划分层状结构，定位软弱层；初至旅行时层析成像（FATT）定位局部断层；系统可靠度分析计算系统失效概率，判断控制滑动面.以2个边坡模型为例，对该方法进行了说明与验证.

1 方法原理

本文提出的基于地质结构探测的多滑面边坡系统可靠度分析包括多道勒夫波分析，初至旅行时层析成像，边坡系统失效概率计算3个模块，基本原理如下.

1.1　多道勒夫波分析(MALW)

为进行边坡稳定性分析，首先需要明确边坡结构，包括软弱层断层位置形态，其中软弱层分布可采用勒夫波探测获取，其原理为通过迭代优化，搜索得最优一维横波速度模型，以最小化观测频散曲线与理论频散曲线误差：

{\lambda }_i=\frac{{‖v_i^{\mathrm{*}}-v_i^{}‖}_{\mathrm{2}}}{m}，(1)

式中：{\lambda }_i为目标函数；v_i^{\mathrm{*}}为基于多道勒夫波信号提取的相速度矩阵，即观测频散曲线；v_i 是基于理论计算得到的相速度矩阵，即理论频散曲线；m为频率数量.

观测相速度矩阵v_i^{\mathrm{*}}采用相移法获取（Park et al.， 1998），其核心概念是在频率域内倾斜累积相位信息，具体步骤如下：步骤1：对信号在距离-时间（x-t）域的记录沿时间方向应用傅里叶变换，并将结果表示为相位和振幅谱的乘积（公式（2））；步骤2：在进行与偏移量相关的相移和振幅归一化后，将每条地震道沿偏移方向的能谱进行叠加（公式（3））；步骤3：根据频率、相位和速度之间的关系，通过坐标变换，将频散能量转换到频率-速度域，得到频散能量图，提取能量集中峰值得观测相速度矩阵v_i^{\mathrm{*}}.

U(x,w)=\int u(x,t){\mathrm{e}}^{iwt}\mathrm{d}t=P(x,w)A(x,w)，(2)

式中：w为圆频率；u(x,t)表示勒夫波在距离（x）-时间（t）域的记录；P(x,w)和A(x,w)分别表示相位谱和振幅谱.

\begin{array}{l}V(w,\varphi )=\int {\mathrm{e}}^{i\varphi w}\left[U(x,w)/\left|U(x,w)\right|\right]\mathrm{d}x=\\ \int {\mathrm{e}}^{-i(\mathrm{\Phi }-\varphi )x}\left[A(x,w)/\left|A(x,w)\right|\right],(\mathrm{3})\end{array}

式中：\varphi表示相位；\mathrm{\Phi }=w/c_w，c_w为对应圆频率w的相速度，对于一个给定的圆频率，当\varphi \approx \mathrm{\Phi }时V(w,\varphi )取模后存在极大值.

理论相速度曲线，v_i 采用Knopoff方法（Knopoff， 1964）计算，该方法基于界面处位移和应力的连续性条件、自由边界条件、无穷远处的辐射条件以及拉普拉斯分解定理，可以推导出勒夫波的频散函数，如下所示：

F_{\mathrm{L}}\left(f,V_{\mathrm{L}}\right)=\left[s-i\right]a_{n-\mathrm{1}}\dots a_{\mathrm{2}}\left[\begin{array}{l}{\left(a_{\mathrm{1}}\right)}_{\mathrm{11}}\\ {\left(a_{\mathrm{1}}\right)}_{\mathrm{21}}\end{array}\right]，(4)

a_m=\left[\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}{Q}_m& i\mathrm{s}\mathrm{i}\mathrm{n}{Q}_m/{\mu }_m{r}_{\beta m}\\ i{\mu }_m{r}_{\beta m}\mathrm{s}\mathrm{i}\mathrm{n}{Q}_m& \mathrm{c}\mathrm{o}\mathrm{s}{Q}_m\end{array}\right]，(5)

式中：V_L为理论勒夫波相速度；n为地层数量；Q_m=k{r}_{\beta m}{h}_m， k表示波数，r_{\beta m}=-i\sqrt[]{\mathrm{1}-(V_{\mathrm{L}}-V_{sm}{)}^{\mathrm{2}}}， V_sm为第m层的剪切波速度，h_m 表示该层的厚度，{\mu }_m表示该层的拉梅常数.

目标函数（公式（1））采用人工蜂群算法（ABC）获取其全局最小值（Karaboga and Akay， 2009）.反演解在限定空间中随机生成（公式（6）），每个解均设有对应的适应度（公式（7）），通过贪婪选择，可获取一个迭代周期中适应度最高的解，若在超过限制次数后仍未更新，则在全局重新随机选择.限制次数（Limit）是ABC算法的一个关键参数，影响算法的探索能力.为搜索全局最小值而避免陷入局部极小值，限制参数的设置取决于群体数量（Akay and Karaboga， 2009）.

v_{ij}=x_{ij}+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\left[-\mathrm{1,1}\right](x_{ij}-x_{kj})，(6)

式中：i,k\in \left\{\mathrm{1,2},...,SN\right\}，j\in \left\{\mathrm{1,2},...,D\right\}，SN是解的数量，D是解的维度；k为随机选定且不同于i，v_{ij}为新生成的解；x_{ij}和x_{kj}是前一步中已有解，对于第一步，x_{ij}和x_{kj}在预设范围内随机选择的.

fi{t}_i=\frac{\mathrm{1}}{{\lambda }_i+\mathrm{1}}，(7)

式中：fit_i 是第i个解的适应度值；{\lambda }_i为目标函数.

1.2　初至旅行时层析成像(FATT)

边坡内部断层位置可通过初至旅行时层析成像获取，该方法通过迭代计算找到一个最优的速度模型，从而最小化理论初至旅行时与观测初至旅行时之间的差异.目标函数定义如下（Taillandier et al.， 2009）：

J\left(\mathbf{c}\right)=\frac{\mathrm{1}}{\mathrm{2}}{\int }_{\partial \mathrm{\Omega }}\mathrm{d}r{\left|T\left(\mathbf{c},r\right)-T_{\mathrm{o}\mathrm{b}\mathrm{s}}\left(r\right)\right|}^{\mathrm{2}}，(8)

式中：T_obs（r）是基于真实速度模型计算或从地震数据中提取的观测初至旅行时，采集面为\partialΩ；T（c， r）是基于速度模型c求解Eikonal方程得到的理论初至旅行时：

{\left|\nabla T\left(x\right)\right|}^{\mathrm{2}}=\frac{\mathrm{1}}{c_{}^{\mathrm{2}}\left(x\right)},(9)

T\left(x_{\mathrm{s}}\right)=\mathrm{0}，(10)

式中：x为观测点；x_s为震源位置；符号∇表示梯度算子.Eikonal方程通过快速扫掠法（Fast Sweeping Method）求解（Zhao， 2004）.

采用伴随状态法来解决目标函数优化问题（Taillandier et al.， 2009），结合公式（9）并引入伴随变量，目标函数公式（8）扩展为：

\begin{array}{l}L\left(\mathbf{c},t,\lambda \right)=\frac{\mathrm{1}}{\mathrm{2}}{\int }_{\partial \mathrm{\Omega }}\mathrm{d}r{\left|t\left(r\right)-T_{\mathrm{o}\mathrm{b}\mathrm{s}}\left(r\right)\right|}^{\mathrm{2}}-\\ \frac{\mathrm{1}}{\mathrm{2}}{\int }_{\mathrm{\Omega }}\mathrm{d}x\lambda \left(x\right)\left({\left|\nabla t\left(x\right)\right|}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathbf{c}{\left(x\right)}^{\mathrm{2}}}\right)\end{array}，(11)

式中：\lambda \left(x\right)是拉格朗日乘子，用于确保在最优解处，Eikonal方程的解t（r） 等于T.在搜索最优解时，\lambda \left(x\right)满足以下条件：

\frac{\partial L}{\partial \mathbf{c}}=-{\int }_{\mathrm{\Omega }}\mathrm{d}x\frac{\lambda (x)}{c^{\mathrm{3}}(x)}=\frac{\partial J}{\partial \mathbf{c}}，(12)

\frac{\partial L}{\partial \lambda }=\mathrm{0}，(13)

\frac{\partial L}{\partial \mathrm{t}}=\mathrm{0}.(14)

通过将L的定义（公式（11））引入到公式（14）中，得\lambda的解为：

\nabla \cdot \Lambda (x)\nabla T(x)=\mathrm{0}.(15)

基于公式（12）得到的梯度，采用准牛顿l-BFGS优化算法（Nocedal and Wright， 2006）更新速度模型：

{\mathbf{c}}_{k+\mathrm{1}}={\mathbf{c}}_k-{\mathbf{D}}_k\cdot {\mathbf{g}}_k,(16)

式中：g_k 为梯度，即公式（12）.D_k 为Hessian矩阵的近似逆矩阵.

1.3　边坡系统失效概率计算

对于存在多潜在滑动面边坡，任一失效模式的发生均会导致系统失效（Jiang et al.， 2015）：

P_{F,\mathrm{s}\mathrm{y}\mathrm{s}}=\mathrm{1}-{\displaystyle \prod _{i=\mathrm{1}}^n}\left[\mathrm{1}-P\left(E_i\right)\right]，(17)

式中：P_F，sys为系统失效概率；n为失效模式数量，即潜在滑动面数量；P（E_i ）为模式i的失效概率，可通过蒙特卡洛法估算：

，(18)

式中：N_i 为蒙特卡洛法产生的总样本；I｛•｝为指示函数，当FS<1时，I（FS <1）=1，否则，取0.在蒙特卡洛模拟中，安全系数（FS）通过响应面函数计算（Bucher and Bourgund， 1990）：

G\left(x\right)=a+{\displaystyle \sum _{\mathrm{1}}^n}{b}_i{x}_i+{\displaystyle \sum _{\mathrm{1}}^n}{c}_i{x}_i^{\mathrm{2}}+\epsilon，(19)

式中：x_i 为随机变量，即岩层强度参数；n为随机变量总个数；a， b_i， c_i，为响应面系数；ε为响应面误差，为均值为0的标准正态分布.

由于岩体内部结构复杂，不同性质结构面及风化程度对其强度参数有较大影响，使得岩体强度参数小于室内完整岩样测得的强度参数，并且存在较大不确定性.其中，在未扰动情况下，黏聚力和内摩擦角可基于弹性波波速比进行折减（王亮清等， 2010）：

K_{\mathrm{c}}=\mathrm{0.113}\mathrm{6}{\beta }^{\mathrm{3}}+\mathrm{0.019}\mathrm{6}{\beta }^{\mathrm{2}}+\mathrm{0.390}\mathrm{5}\beta +\mathrm{0.474}\mathrm{7}，(20)

K_{\phi }=\mathrm{0.006}\mathrm{7}{\beta }^{\mathrm{3}}+\mathrm{0.187}\mathrm{9}{\beta }^{\mathrm{2}}+\mathrm{0.381}\mathrm{7}\beta +\mathrm{0.424}\mathrm{1}，(21)

式中：K_c为软岩黏聚力折减系数；K_φ 为软岩内摩擦角折减系数；β为弹性波波速比，本文采用勒夫波反演横波速度与完整岩石横波速度比值.

设各层强度参数服从对数正态分布（Peng et al.， 2014），其均值及标准差通过对各测点折减后强度参数统计获取.用以拟合响应面函数的随机数样本点基于Box-Behnken设计采样（Ferreira et al.， 2007），该方法将变量设置为-1，0，1三个水平，分别表示μ-σ，μ，μ+σ，其中μ为均值，σ为标准差.通过不同组合获取变量间关系，例如组合（+1， +1， 0）， （-1， -1， 0）， （-1， +1， 0）和（+1， -1， 0） 反映前两个因素相互作用.采样点总数量可按照公式（22）计算.样本安全系数采用有限差分软件FLAC3D基于强度折减法计算.综上，提出方法具体流程介绍如图1所示.

N=\mathrm{2}k\left(k-\mathrm{1}\right)+c，(22)

式中：k为因素数量，c为中心点数量，常设为3.

（1）基于多道勒夫分析模块划分边坡层状结构.在无断层散射波记录情况下，采用相移法提取频散曲线，基于人工蜂群算法反演层状速度结构，定位边坡软弱层，并获取软弱层横波速度数据集.

（2）基于初至旅行时层析成像模块定位断层低速带.在存在断层散射波记录情况下，于坡表增设初至旅行时测试.对面波未明显受到断层散射波干扰的部分进行MALW分析，基于泊松比将横波速度转换为纵波速度，获取层状结构作为FATT初始模型，采用伴随状态法成像断层低速带.

（3）基于边坡系统失效概率计算模块，判断边坡控制滑面，计算单一及系统失效概率.首先，基于各测点软弱层反演横波速度与完整岩石横波速度比值，折减软岩强度参数，并通过统计分析计算其均值和方差，构建概率分布模型.随后，采用Box-Behnken设计方法选取随机样本，通过FLAC3D计算地表位移及安全系数用以拟合响应面函数，采用蒙特卡洛模拟得地表位移及安全系数概率分布，识别控制滑动面，评估单一及系统失效概率.

2 算例分析

2.1　边坡模型

算例边坡模型高40 m，宽110 m，坡度设为45°，稳定性控制因素分别为软弱层（图2a）及断层软弱层组合（图2b）.为了更准确地反映真实边坡特征，模型构建时考虑地层起伏及断层区域速度渐变.地层起伏通过在倾斜层面基础上，添加均值为0的随机上下波动（公式（23）），并采用高斯平滑函数对界面向量进行平滑处理，以避免不连续性.断层可看作是由中心断层和周边破损区组成的复杂区域（图2b），其总宽度与断层滑动位移相关（Choi et al.， 2016），如公式（24）所示.在破损区内，次生滑移和裂隙从内部到外部逐渐减少.因此横波速度设置为在中心断层最低，并沿破损区半径逐渐增加，直到达到围岩的速度.模型参数设置如表1所示 （Peng et al.， 2014； 郑志勇等， 2016；Sun et al.， 2019）.

Y_n(x)=a_n\times (x-x_n)+d_n+y_u(x),n=\mathrm{1,2},...,N，(23)

式中：Y_n （x）表示第n层界面的纵坐标向量；x_n， d_n 分别是倾斜界面时选用的第n层界面基准点横坐标和纵坐标；a_n 为第n层界面倾斜斜率；y_u(x)是在波动范围内均值为零的随机数；N为地层界面的数量.

w=\mathrm{1.659}\mathrm{5}\mathrm{\Delta }{z}^{\mathrm{0.687}\mathrm{2}},R^{\mathrm{2}}=\mathrm{0.762}\mathrm{2}，(24)

式中：\mathrm{\Delta }z是滑移位移，w是包括周边破损区的断层总宽度.

勒夫波及初至旅行时的激发和接收都位于坡表.多道勒夫波采用线性阵列，偏移距10 m，道间距2 m（图2c）.线性阵列测点位于接收器中点，设置于坡顶水平面，共12个，间距5 m（图2a）.初至旅行时激发震源设置于坡面，垂直距离为2.5 m，接收器沿边坡顶布设，水平间距2.5 m（图2b）.震源和接收器的数量分别为15个和25个.

2.2　勒夫波模拟

勒夫波在地下介质中的传播过程采用二阶时间导数和十阶空间导数交错有限差分网格模拟（Jean， 1984）.网格大小为0.5 m，时间步长为0.02 ms.震源子波采用60 Hz的Ricker波.边坡倾斜自由表面近似为阶梯状，如图3蓝色折线所示，折线下方为固体介质，上方为空气.自由边界设置在声-弹性方法（AEA）（Xu et al.， 2007）基础上进行了改进以适用与勒夫波的传播.根据AEA方法，自由表面上存在4类网格：（1）水平网格（H），（2）右侧为空气的垂直网格（VR），（3）上方和右侧为空气的内角网格（IR），（4）上方和右侧为空气的外角网格 （OR）.该4类网格在边坡中进行勒夫波模拟时，采用了两步处理方案设置.第1步是参数预设，自由界面上各类网格的密度和拉梅系数如表2所示.{\rho }_{\mathrm{1}}和{\mu }_{\mathrm{1}}分别表示相邻固体网格的密度和拉梅系数.第2步是在有限差分法计算过程中固定剪切应力.水平网格（H）上的应力分量{\tau }_{xy}和垂直网格（VR）上的应力分量{\tau }_{zy}在计算过程中固定为零.在IR和OR点上，{\tau }_{xy}和{\tau }_{zy}分别固定为0.在除自由表面外的其他边缘，例如坡底及右侧界面，应用多轴完美匹配层（MPML）吸收技术来抑制伪反射（Meza-Fajardo and Papageorgiou， 2008）.

2.3　边坡层状结构反演

模型Ⅰ测点P₁消除坡面反射波后的多道勒夫波记录如图4a所示，基于相移法获取的频散能量图存在多阶模式（图4b），频散能量在约70 Hz处从基阶跳跃到高阶.“跳跃”现象是由软弱层内生成的LVL导波（low velocity layer）造成（Mi et al.， 2018），LVL导波在低速层中生成并传播，其中短波长、高频率、低相速度的波无法逸出低速层.因此，各模式在高频段缺乏能量，频散能量图出现“跳跃”现象.由于“跳跃”现象以及高阶模式对深部速度模型的敏感性大于基模式这一特征（Luo et al.， 2007），本文采用多模式反演，以保证分辨率.同时，观察发现，在低频段，频散能量与理论值存在差异，这种现象是由泄漏表面波引起的，该波在低频范围内与基阶能量连续，因此基阶频散能量在低频段不可取（Gao et al. 2014）.

基于多阶频散曲线，采用人工蜂群算法（ABC）在约束范围内进行全局寻优反演得一维横波速度曲线（图5a），可确定软弱层深度，波速及厚度.通过对12个测点一维速度曲线线性插值生成二维横波速度剖面，如图5b所示.边坡分层结构清晰可见，存在两个坡度平缓软弱层，深度和趋势与合成模型较为一致.

2.4　边坡局部断层成像

受断层散射波影响的多道勒夫波信号如图6所示，该散射波的形成原理与LVL导波类似，在断层中生成并传播的导波， 经过多次反射和折射后，具有更长波长和较低频率的波可以从低速层逸出，被地震检波器接收.当存在断层散射波时，可初步判断存在断层等横向速度异常体，增设初至旅行时层析成像分析.

断层导致观测初至波旅行时记录在35~55 m之间出现波动（图7c）.FATT初始模型为不包含断层的层状速度结构，由对面波未明显受到断层散射波干扰的部分进行MALW分析反演，依据趋势延伸，并基于表1中泊松比转换得（图7a），其计算初至旅行时无波动特征（图7c）.经过7次迭代后，误差达到收敛（图8），反演初至旅行时曲线在35~55 m之间的波动区域与观测初至旅行时趋势一致（图7d），这表明FATT有效捕捉了断层结构特征.更新后速度剖面如图7b所示，局部断层位置与倾斜程度与设置模型一致，根据速度剖面，可判断边坡存在两个潜在滑动面，S₁与S₂.

2.5　边坡失效概率分析

边坡模型Ⅱ存在两种破坏模式，即图7b中所示的S₁和S₂滑动面，任一破坏模式的发生将会引起边坡的系统失效.在系统可靠度分析过程中，软弱层2，4和断层的黏聚力及内摩擦角设为符合对数正态分布的随机变量，其中软弱层2，4黏聚力及内摩擦角均值标准差通过统计基于公式（20）和公式（21）折减后强度参数数据集获取，断层黏聚力及内摩擦角均值标准差参考Peng et al.（2014）研究取值，如表3所示.

不同强度参数组合样本点的地表位移及安全系数采用FALC3D计算.FLAC3D模型长140 m，宽1 m，高60 m，由5 158个节点，15 508个四面体单元组成（图9）.考虑地震工况，采用拟静力法，水平系数设为0.075，并在顶部承受300 kPa均匀荷载.模型左侧和底部为固定边界，顶部及右侧斜面为自由面.在边坡顶部和倾斜面上共设置11个测量点，测点x方向的间距为10 m.由于FATT反演结果中的断层厚度被高估，选取低速带的中央条状区域作为断层，其厚度为低速带总宽度的1/3.

基于FALC3D计算结果拟合的测点D₈，D₉及D₁₀处X方向地表位移响应面系数如表4所示，响应面随机变量为2个软弱层及断层的黏聚力与内摩擦角，Box-Behnken设计采样数量为63组.为验证响应面准确性并计算响应面误差，选取12个随机点对比响应面预测位移结果与数值模拟预测结果，如图10a所示.基于响应面函数进行1 000 000次蒙特卡洛模拟得测点D₈处X方向地表位移均值为9.43 mm，标准差为2.40 mm（图11a），测点D₉处均值为8.29 mm，标准差为3.10 mm（图11b），测点D₁₀处为0.61 mm和0.02 mm（图11c）.测点D₈和D₉位移数值接近且远大于D₁₀处位移， 可判断测点D₈和D₉同时发生滑动，边坡沿S₂发生破坏.

为验证地表位移判断结果并计算边坡系统可靠度，对2个滑动面安全系数建立响应面函数，如表4中所示，FS₁对应于滑动面S₁（图12a），变量为第2层及断层的黏聚力内摩擦角，FS₂对应于滑动面S₂（图12b），变量为第4层及断层的黏聚力内摩擦角，两者Box-Behnken设计采样数量均为27组.响应面预测值与数值模拟计算值基本一致，对比如图10b所示.基于拟合的响应面函数，FS₁及FS₂均在10 000次蒙特卡洛模拟内达到收敛（图13a~13c）.基于10 000次蒙特卡洛模拟得边坡沿S₂滑动失效概率为5.14×10^-3 （图13d）.为捕捉小概率失效事件，边坡沿S₂滑动失效概率采用1 000 000次蒙特卡洛模拟计算，得失效概率为1.0×10^-7（图13d）.在本研究中，使用配备18核Intel® Core™ i9-14900HX处理器和32 GB内存的笔记本电脑运行10 000次及1 000 000次蒙特卡罗模拟分别耗时约0.02~890.00 s.最后，根据公式（17），边坡系统失效概率为5.140 099×10^-3，与沿S₂滑动失效概率极为接近，边坡稳定性受潜在滑动面S₂控制，与地表位移判断一致.

3 讨论

3.1　勒夫波与瑞雷波频散能量对比

勒夫波频散能量与瑞雷波频散能量相比，具有能量分布均匀.受边坡岩层界面波动影响小两大优势.构建了3个层面波动范围（R）分别为0 m，5 m和10 m的含软弱层边坡模型（图14a），速度参数列于表5中.震源采用20 Hz的Ricker子波，最小偏移量和地震检波器间距分别为30 m和1 m，震源与检波器均位于坡顶.通过比较相移法获取的瑞雷波频散能量和勒夫波频散能量（图14），可以发现勒夫波频散能量图像有两个优势.第1个优势是其简单性.勒夫波每个模式的能量分布均匀，图像更清晰，便于提取频散曲线.然而，对于瑞雷波，在频散图像中出现了“模式接吻”的现象，这可能导致模式识别错误，导致反演速度不可靠.勒夫波是由SH波的多次反射产生，瑞雷波是由P波和S波相互作用产生，因不涉及P波，故而勒夫波的频散能量更简单清晰（Xia et al.， 2012）；第2个优势是受边坡岩层界面起伏影响较小.随着波动范围的增加，瑞雷波的高阶频散能量变得不连续且混乱，难以提取频散曲线，并且容易导致模式识别错误.然而，对于勒夫波频散能量，前三阶几乎不受影响.这种现象是由两者质点运动轨迹的不同导致.瑞雷波质点运动轨迹是在边坡剖面内的逆行椭圆，勒夫波质点运动轨迹是水平的、垂直于边坡剖面的横向运动.因此，边坡剖面的几何变化对勒夫波传播的影响较小.但在实际边坡勘探应用中，垂直方向的瑞雷波相比水平方向的勒夫波更易于激发和采集，需依据边坡场地条件选择.

3.2　基阶反演与多阶反演对比

多阶反演对边坡浅部及深部软弱层探测准确率均大于基阶反演.图15显示了边坡模型Ⅰ，测点P₁基于基阶和多阶频频散曲线分别反演得到的一维速度剖面.基阶反演浅部软弱层分辨率较差，因为其频散曲线在66 Hz处被截断，导致反映浅层地质信息的高频信息被遗漏.同时高阶频散曲线对深层参数更为敏感（Luo et al.， 2007），多阶反演深部软弱层的探测精度也高于基阶反演.相比仅依赖基阶数据的反演，多阶反演能够获得更优的结果.然而，如上节所示，高阶模式能量存在“模式接吻”及混乱现象，易造成模式提取错误，因此，在边坡地质结构探测中，需避免使用不可靠的高阶模式.

3.3　变异系数对失效概率影响

随着黏聚力和内摩擦角变异系数的增加，系统失效概率逐渐增大（图16）.这是由于变异系数越大，内摩擦角及黏聚力取到最小值的概率也随之增加，导致单一滑动面失效概率增大，进而引起系统失效概率上升.内摩擦角对系统失效概率的影响大于黏聚力，基于摩尔库伦剪切强度公式τ=c+σtanφ，黏聚力对安全系数影响相对直接且线性，而内摩擦角以正切形式影响安全系数.同时，随着深度增加，由自重引起的正应力σ增大，使得内摩擦角影响程度更为显著.

4 结论

本文提出了基于地质结构探测的多滑面边坡系统可靠度分析方法，用于定位边坡结构面，量化强度参数不确定性，分析边坡失效概率及滑动面，主要结论如下：

（1）提出的方法基于多道勒夫波分析（MALW）和初至旅行时层析成像（FATT）探测边坡地质结构，利用两者的互补优势，可定位软弱层与断层；基于弹性波波速比折减软弱层黏聚力和内摩擦角，统计分析可得具有场地特异性的参数概率分布模型；针对多潜在滑动面的边坡模型，结合地表位移概率分布与地质结构信息，可识别控制滑动面，通过可靠度分析方法评估各潜在滑动面及系统失效概率.

（2）多道勒夫波分析（MALW）采用线性整列，震源及检波器均布置于坡顶.多阶频散曲线对边坡浅部及深部软弱层的反演精度均高于基阶.与瑞雷波相比，勒夫波频散能量分布更均匀，受边坡岩层界面波动影响小，更易提取高阶频散曲线.

（3）初至旅行时层析成像（FATT）震源设置于坡面，接收器设置于坡顶.断层对初至旅行时记录的影响表现为特定范围内的波动特征，通过对该波动特征的拟合，可实现边坡局部断层的成像.

（4）基于弹性波波速折减统计所得的强度参数概率分布，充分考虑了软弱层内部裂隙结构和风化程度的影响，体现了场地特异性.通过地表位移概率分布对比与系统可靠度分析相结合，可识别关键控制滑动面，定量评估各潜在滑动面的失效概率，为边坡防治提供了科学参考.在本文边坡案例中，经位移判断，边坡受深层滑动面控制，失效概率为5.14×10^-3，浅层滑动面失效概率仅为1.0×10^-7，故其系统可靠度主要由深层滑动面失效概率所决定，且受内摩擦角变异系数影响远大于黏聚力.

参考文献

[1] Akay, B., Karaboga, D., 2009. Parameter Tuning for the Artificial Bee Colony Algorithm.Springer, Berlin, Heidelberg, 608-619.

[2] Bucher, C. G., Bourgund, U., 1990. A Fast and Efficient Response Surface Approach for Structural Reliability Problems. Structural Safety, 7(1): 57-66. https://doi.org/10.1016/0167-4730(90)90012-E

[3] Choi, J. H., Edwards, P., Ko, K., et al., 2016. Definition and Classification of Fault Damage Zones: A Review and a New Methodological Approach. Earth-Science Reviews, 152: 70-87. https://doi.org/10.1016/j.earscirev.2015.11.006

[4] Crosta, G. B., di Prisco, C., Frattini, P., et al., 2014. Chasing a Complete Understanding of the Triggering Mechanisms of a Large Rapidly Evolving Rockslide. Landslides, 11(5): 747-764. https://doi.org/10.1007/s10346-013-0433-1

[5] Dou, J., Xiang, Z. L., Xu, Q., et al., 2023. Application and Development Trend of Machine Learning in Landslide Intelligent Disaster Prevention and Mitigation. Earth Science, 48(5): 1657-1674 (in Chinese with English abstract).

[6] Dou, J., Xing, K., Wang, L. Z., et al., 2025. Air-Space-Ground Synergistic Observations for Rapid Post-Seismic Disaster Assessment of 2025 Ms6.8 Xigazê Earthquake, Xizang. Journal of Earth Science. https://doi.org/10.1007/s12583-025-0160-2

[7] Dunne, J., Beresford, G., 1995. A Review of the τ-p Transform, Its Implementation and Its Applications in Seismic Processing. Exploration Geophysics, 26(1): 19-36. https://doi.org/10.1071/eg995019

[8] Ferreira, S. L. C., Bruns, R. E., Ferreira, H. S., et al., 2007. Box-Behnken Design: An Alternative for the Optimization of Analytical Methods. Analytica Chimica Acta, 597(2): 179-186. https://doi.org/10.1016/j.aca.2007.07.011

[9] Gao, L. L., Xia, J. H., Pan, Y. D., 2014. Misidentification Caused by Leaky Surface Wave in High-Frequency Surface Wave Method. Geophysical Journal International, 199(3): 1452-1462.https://doi.org/10.1093/gji/ggu337

[10] Guo, R., Li, M. K., Yang, F., et al., 2019. First Arrival Traveltime Tomography Using Supervised Descent Learning Technique. Inverse Problems, 35(10): 105008. https://doi.org/10.1088/1361-6420/ab32f7

[11] Jean, V., 1984. SH-Wave Propagation in Heterogeneous Media; Velocity-Stress Finite-Difference Method. Geophysics, 49(11): 1933-1942. https://doi.org/10.1190/1.1441605

[12] Jiang, S. H., Liu, Y., Zhang, H. L., et al., 2020. Quantitatively Evaluating the Effects of Prior Probability Distribution and Likelihood Function Models on Slope Reliability Assessment. Rock and Soil Mechanics, 41(9): 3087-3097 (in Chinese with English abstract).

[13] Jiang, S.H., Li, D.Q., Cao, Z.J., et al., 2015. Efficient System Reliability Analysis of Slope Stability in Spatially Variable Soils Using Monte Carlo Simulation. Journal of Geotechnical and Geoenvironmental Engineering, 141(2): 04014096. https://doi.org/10.1061/(asce)gt.1943-5606.0001227

[14] Karaboga, D., Akay, B., 2009. A Comparative Study of Artificial Bee Colony Algorithm. Applied Mathematics and Computation, 214(1): 108-132. https://doi.org/10.1016/j.amc.2009.03.090

[15] Knopoff, L., 1964. A Matrix Method for Elastic Wave Problems. Bulletin of the Seismological Society of America, 54(1): 431-438. https://doi.org/10.1785/BSSA0540010431

[16] Li, J. P., Cheng, Y. G., Li, D. Q., et al., 2016. System Reliability Analysis of Spatially Variable Soil Slopes Using the Multiple Response Surfaces Method. Rock and Soil Mechanics, 37(1): 147-155, 165(in Chinese with English abstract).

[17] Li, W. L., Xu, Q., Lu, H. Y., et al., 2019. Tracking the Deformation History of Large-Scale Rocky Landslides and Its Enlightenment. Geomatics and Information Science of Wuhan University, 44(7): 1043-1053 (in Chinese with English abstract).

[18] Liu, G. F., Feng, K., Yan, C. G., et al., 2023. Probabilistic Evaluation of Excavation Unloading Response of Rock Slope Considering the Uncertainty of Mechanical Parameters. Rock and Soil Mechanics, 44(7): 2115-2128 (in Chinese with English abstract).

[19] Liu, H., Zhou, H. W., Liu, W. G., et al., 2010. Tomographic Velocity Model Building of the near Surface with Velocity-Inversion Interfaces: A Test Using the Yilmaz Model. Geophysics, 75(6): U39-U47. https://doi.org/10.1190/1.3502665

[20] Liu,G. F., Feng, K., Yan, C. G., et al., 2023. Probabilistic Evaluation of Excavation Unloading Response of Rock Slope Considering the Uncertainty of Mechanical Parameters.Rock and Soil Mechanics, 44(7): 2115-2128 (in Chinese with English abstract).

[21] Luo, Y. H., Xia, J. H., Liu, J. P., et al., 2007. Joint Inversion of High-Frequency Surface Waves with Fundamental and Higher Modes. Journal of Applied Geophysics, 62(4): 375-384. https://doi.org/10.1016/j.jappgeo.2007.02.004

[22] Luu, K., Noble, M., Gesret, A., et al., 2018. A Parallel Competitive Particle Swarm Optimization for Non-Linear First Arrival Traveltime Tomography and Uncertainty Quantification. Computers & Geosciences, 113: 81-93. https://doi.org/10.1016/j.cageo.2018.01.016

[23] Meza-Fajardo, K. C., Papageorgiou, A. S., 2008. A Nonconvolutional, Split-Field, Perfectly Matched Layer for Wave Propagation in Isotropic and Anisotropic Elastic Media: Stability Analysis. Bulletin of the Seismological Society of America, 98(4): 1811-1836. https://doi.org/10.1785/0120070223

[24] Mi, B. B., Xia, J. H., Shen, C., etal., 2018. Dispersion Energy Analysis of Rayleigh and Love Waves in the Presence of Low-Velocity Layers in Near-Surface Seismic Surveys. Surveys in Geophysics, 39(2): 271-288. https://doi.org/10.1007/s10712-017-9440-4

[25] Nocedal, J., Wright, S. J., 2006. Quasi-Newton Methods. In: Numerical Optimization. Springer Series in Operations Research and Financial Engineering. Springer, New York.

[26] Pan, Y. D., Gao, L. L., Bohlen, T., 2018. Time-Domain Full-Waveform Inversion of Rayleigh and Love Waves in Presence of Free-Surface Topography. Journal of Applied Geophysics, 152: 77-85. https://doi.org/10.1016/j.jappgeo.2018.03.006

[27] Park, C. B., Miller, R. D., Xia, J. H., 1998. Imaging Dispersion Curves of Surface Waves on Multi-Channel Record. SEG Technical Program Expanded Abstracts, 1377-1380. https://doi.org/10.1190/1.1820161

[28] Park, C. B., Miller, R. D., Xia, J. H., 1999. Multichannel Analysis of Surface Waves. Geophysics, 64(3): 800-808. https://doi.org/10.1190/1.1444590

[29] Peng, M., Li, X. Y., Li, D. Q., et al., 2014. Slope Safety Evaluation by Integrating Multi-Source Monitoring Information. Structural Safety, 49: 65-74. https://doi.org/10.1016/j.strusafe.2013.08.007

[30] Sun, Y., Huang, J. S., Jin, W., et al., 2019. Bayesian Updating for Progressive Excavation of High Rock Slopes Using Multi-Type Monitoring Data. Engineering Geology, 252: 1-13. https://doi.org/10.1016/j.enggeo.2019.02.013

[31] Taillandier, C., Noble, M., Chauris, H., et al., 2009. First-Arrival Traveltime Tomography Based on the Adjoint-State Method. Geophysics, 74(6): WCB1-WCB10. https://doi.org/10.1190/1.3250266

[32] Wang, L. Q., Liang, Y., Wu, Q., et al., 2010. Estimation Method of Rockmass Strength Parameters Based on Elastic Wave Velocity Ratio. Coal Geology & Exploration, 38(6): 37-42 (in Chinese with English abstract).

[33] Wang, X. R., Wang, H., Liang, R. Y., 2018. A Method for Slope Stability Analysis Considering Subsurface Stratigraphic Uncertainty. Landslides, 15(5): 925-936. https://doi.org/10.1007/s10346-017-0925-5

[34] Wang, Z. N., Sun, C. Y., Wu, D. S., 2022. Near-Surface Site Characterization Based on Joint Iterative Analysis of First-Arrival and Surface-Wave Data. Surveys in Geophysics, 44(2): 357-386. https://doi.org/10.1007/s10712-022-09747-8

[35] Xia, J. H., 2014. Estimation of Near-Surface Shear-Wave Velocities and Quality Factors Using Multichannel Analysis of Surface-Wave Methods. Journal of Applied Geophysics, 103: 140-151. https://doi.org/10.1016/j.jappgeo.2014.01.016

[36] Xia, J. H., Xu, Y. X., Luo, Y. H., et al., 2012. Advantages of Using Multichannel Analysis of Love Waves (MALW) to Estimate Near-Surface Shear-Wave Velocity. Surveys in Geophysics, 33(5): 841-860. https://doi.org/10.1007/s10712-012-9174-2

[37] Xie, H., Liu, L. B., 2015. Near-Surface Anisotropic Structure Characterization by Love Wave Inversion for Assessing Ground Conditions in Urban Areas. Journal of Earth Science, 26(6): 807-812. https://doi.org/10.1007/s12583-015-0619-7

[38] Xu, Y. X., Xia, J. H., Miller, R. D., 2007. Numerical Investigation of Implementation of Air-Earth Boundary by Acoustic-Elastic Boundary Approach. Geophysics, 72(5): SM147-SM153. https://doi.org/10.1190/1.2753831

[39] Yang, Z. Y., Yin, C. C., Nie, J. Y., et al., 2022. Probability Density Function Estimation of Geotechnical Parameters Considering the Three-Dimensional Spatial Variability Based on Multi-Source Site Investigation Data. Rock and Soil Mechanics, 43(6): 1571-1584 (in Chinese with English abstract).

[40] Zhao, H. K., 2004. A Fast Sweeping Method for Eikonal Equations. Mathematics of Computation, 74(250): 603-627. https://doi.org/10.1090/s0025-5718-04-01678-3

[41] Zheng, Z. Y., Yu, H. B., Xu, H. Q., 2016. Numerical Analysis of Failure Modes and Stability of Soft and Hard Rock Interbedded Slope. Journal of Yangtze River Scientific Research Institute, 33(9): 102-106 (in Chinese with English abstract).

基金资助

国家自然科学基金-联合基金重点项目(U23A2044)

国家自然科学基金-国际合作项目(42061160480)

国家自然科学基金-面上项目(42071010)

国家自然科学基金-面上项目(42477195)