基于Box-Jenkins 随机模型的滑坡稳定性预测模型

Nonlinear Prediction of Landslide Stability Based on Machine Learning

摘要

滑坡稳定性含有非线性特征，机器学习算法相对于传统算法在滑坡稳定性预测中准确性更高.因此，为了更加准确地分析顺层岩质边坡在循环地震荷载作用下的稳定性，结合室内物理模型试验和离散元数值模拟软件PFC3D相互对比的研究手段，得到了滑带土的应变软化过程；并利用滑坡变形的非线性特点，和数值模拟得到的滑坡稳定性系数数据，经过模式识别、拟合检验等提出了基于机器学习算法（Box-Jenkins 随机模型）的滑坡稳定性预测模型.结果表明：（1）剪切应力的逐渐减小促进了滑带土应变的软化过程，滑带土的围压虽然能抑制滑带土裂缝的增加，但对应变软化的抑制作用有限；（2）本研究所建立的标准BIC值为8.160的ARIMA（1，1，0）（0，1，1）模型，可以对边坡稳定系数时间序列数据进行精准预测.基于边坡稳定系数和应力场的现场观测，进一步描述了两种可能的滑坡触发机制，同时时间序列的机器学习能准确预测循环荷载作用下边坡稳定系数的变化规律.

关键词

机器学习 / 滑坡 / 稳定性计算 / 时间序列 / 工程地质

Key words

machine learning / landslide / stability calculation / time series / engineering geology

滑坡及其相关的边坡破坏是世界上最主要的地质灾害之一，每年都会造成巨大的经济损失和人员伤亡.大量学者对滑坡灾害的预测与防治进行了深入的研究，在滑坡区间预测取得了较大进展.基于机器学习的非线性预测方法因其强大的非线性映射能力在滑坡变形预测领域越来越受欢迎（俞布等， 2012； Guo et al.， 2019， 2022； Vm et al.，2021； 陈志荣等， 2022；孙畅等，2022；王闯等，2022）.

滑带土失稳破坏的动态变形是一个复杂的、随时间渐进的过程，与滑带的岩土及力学性质密切相关（Chigira and Yagi， 2006； Kirschbaum et al.， 2015； Xiao et al.， 2022），随着科学研究尺度的日益细化，滑带土的细观尺度特征逐渐受到广大学者的关注.关于滑带土的细观劣化特征研究，主要有两方面的研究成果：（1）一些学者利用现代高科技仪器对滑带土的矿物组成进行了分析测试，观察了滑带土的微观结构特征.如Wang et al. （2019）研究发现滑带土的微观结构在剪应力作用下形成一种定向排列的扁平粘土颗粒组；Yang et al. （2020） 采用CT电子扫描仪等设备研究得到了滑带土在滑坡发生的不同阶段的定向排列情况.（2） 数值模拟研究在某种意义上是宏观现象的精细补充研究.如Shan and Lai （2019）采用数字图像处理技术和离散元软件结合进一步研究得到了土石混合体的自动生成方法.Zhang et al. （2016）利用PFC3D模拟了无黏性粗颗粒与充填料浆介质的混合滑带土，得到了滑带土的力学变形机制，包括填充介质引起的粗颗粒之间的碰撞和颗粒粘附.胡峰等（2018）研究了土石混合体中剪切带的力学行为，并在随机多边形岩石块体的基础上发展了细观结构.然而现有研究对滑带土细观变形特性的描述非常有限，而且对滑带土的细观结构特征、细观剪切机制及其对宏观力学特性的影响研究也较少.

而对于边坡的安全防护，首先要做的是明确边坡的安全状态，但现有的研究工具大多是通过传统的定量求解方法获得边坡的稳定状态，具有一定的局限性.边坡安全包含很多参数和随机变化规律，计算过程相对复杂，而随着机器学习技术的迅速兴起，为边坡安全稳定研究提出了一种新的思路.决策树（Thai Pham et al.， 2018）、随机森林、支持向量机（Samui， 2008； Peethambaran et al.， 2020； Balogun et al.， 2021）和普通贝叶斯（蒋水华等，2020；黄发明等， 2018， 2021）等算法在边坡研究中得到了广泛的应用.其中，人工神经网络和多元回归分析作为机器学习的两个经典模型，被广泛应用于解决各种预测问题，但这两种模型的影响因素权重分析不够全面，缺乏基于相关检验的权重分析与基于智能模型的权重分析的综合比较，同时没有考虑到时间序列数据的随机性和周期性，导致预测误差较大.而Box-Jenkins 随机模型有效的避免了上述提到的问题，Pinzón-Hassan et al. （2021）和 Ratshefola-Lerefolo and Mpeta （2023）研究发现Box-Jenkins 随机模型能很好地研究对象和时间之间的随机关系，同时能结合时间序列数据，较容易得到一种精度较高且切实有效的方法模型，如陈卫雄（2021）基于Box-Jenkins建模方法，确定了时间序列模型阶数，准确地预测了青藏铁路路基高程的变化.

机械学习模型种类较多，如果机器模型选择不完美的模型结构，在参数选择和数据噪声误差等方面均会导致预测结果较差，而本文采用Box-Jenkins 建模方法对边坡稳定性时间序列数据进行预测，不仅对非平稳的时间序列数据进行了平稳处理，同时增加了预测结果的可靠程度.综上所述，本文结合数值模拟和模型试验，对比分析滑带土在循环荷载作用下的宏细观劣化特征；并采用机器学习的方法建立数学预测模型，研究边坡循环荷载作用下的稳定性特征.

1 数值方法

本文为了研究循环荷载作用下含有预定剪切面的滑坡的演化机制，结合数值模拟和室内物理模型试验，深入分析了循环荷载作用下剪切应力与土体细尺度损伤变量之间的关系；同时，将离散元数值模拟的结果与边坡稳定系数的变化进行比较，分析模拟结果对滑坡形成的影响，最后，基于数值模拟和模型试验得到的研究结果，结合Box-Jenkins建模方法，利用机器学习从边坡稳定系数的时间序列数据出发，全面系统地预测不同循环荷载作用下的边坡稳定性，研究路线如图1所示.

1.1　滑动带土循环剪切破坏行为的概念模型

基于循环剪切破坏行为的假设，概念模型必须能够（i）允许颗粒间链接破裂，同时可以（ii）表征土壤的微观结构特征.对于第一个方面，考虑到离散元法在表征微裂缝方面的独特能力，采用线性平行键模型来模拟土体的应变软化特征，该方法避免了从传统经验数据中获得固有模型来研究材料的力学性能，而从细尺度角度研究块状介质的力学行为，克服了传统连续介质力学模型的宏观连续性假设（Alaei and Mahboubi， 2012）.因此，为了考虑循环荷载作用下滑带土的剪切破坏力学特性，提出了滑带土数值模拟的详细方案，在实验室尺度观测土体应变软化行为试验和滑带土体数值模拟的概念如图2所示，校准的微观参数如表1所示.

PFC模拟的一个优点是能够获得从基于连续的技术或物理实验中无法获得的信息.接触数的变化描述了断层带的破坏损伤过程，假设在数值模型中的测量区域或边界上有N个具有质心的粒子，n_c 为每个粒子的接触数 （Itasca， 2014；张志飞等， 2022）.C_n 定义为每个体的平均活动接触数为配位数（Zhang et al.， 2021），一旦两个颗粒之间的接触（颗粒内或颗粒间）由于过度的接触力而破裂，即而引发微断裂.接触数的监测范围如图3所示，接触数的减少意味着颗粒接触破裂，即压裂，因此，接触数的减少速率和减少幅度分别反映了压裂速率和压裂程度，如式（1）所示：

C_n = \raisebox{1ex}{$\sum _N{n}_c$}\!\left/ \!\raisebox{-1ex}{$N$}\right.,(1)

为了比较不同条件下配位数的变化，取配位数变化百分比P_C，如式（2）所示.

P_C = ( C_n - C_{ni} ) / C_{ni} \times \mathrm{100},(2)

式中C_n 、C_ni 分别为仿真模型中的当前配位数和初始配位数.

1.2　模型校准

滑带土数值模拟试验过程包括三部分：（1）数值模型的应力初始化；（2）数值模型在预定围压下的伺服控制；（3）对数值模型施加循环荷载.室内实验室试验中循环加载主要是通过输入振动波来控制试验装置，而在PFC3D中主要是通过壁速的定义来实现循环剪切加载.

在PFC3D模拟中，通过监测受荷墙与颗粒之间的应力值，获得受荷过程中滑带土剪切应力的大小.根据数值模型中的振动参数，主要对下部剪切盒进行直剪试验，并输入正弦波实时调整墙体运动方向，从而实现对滑带土的循环加载.其中图2中的1#、2#、3#、4#、5#、9#和10#墙指定了相同的剪切速率.循环加载数值模拟中考虑的参数中循环加载次数包括12个水平，即750次、1 500次、 2 250次、3 000次、3 750次、4 500次、5 250次、 6 000次、7 500次、9 000次、11 500次、12 000次，循环加载幅值包括4因素，即0.5 cm/s、0.9 cm/s、1.5 cm/s、2.2 cm/s，并且按照正交试验的原理开展4因素12水平的正交试验.在数值模拟过程中，PFC3D程序主要包括数值模型剪应力的测量和循环加载过程中墙体速度的应用.剪力盒如图3所示，下部剪力盒受力为墙体8上的接触力与墙体6上的接触力之和，上部剪力盒受力为墙体9上的接触力与墙体4上的接触力之和.

PFC模型的宏观行为是粒子及其接触的微观性质的表现.为模拟材料选择一组合适的微观参数，即模型校准，通常是通过将数值响应与实验室观测匹配来实现（Potyondy and Cundall， 2004）.PFC3D所涉及的力学参数是表征颗粒性质的精细参数，具有随机性，与宏观力学性能的关系复杂，当通过基础力学试验得到的宏观力学性能与实际试验结果一致时，一般认为标定参数是合理的（Itasca， 2014）.因此，在滑带土宏观力学参数的基础上，模拟滑带土的剪切应力‒应变曲线，找出并确定相应曲线对应的数值模型的抗剪强度的细尺度参数.大多数校准程序通常遵循Hofmann et al. （2015）提出的方案，即通过正交试验原理，采用试错法原理，将不同参数得到的结果进行正交分析，继而逐步逼近最终计算得到的参数，其中数值应力应变响应与实验室实验结果的良好匹配如图4所示.校准的微观参数如表1所示，通过图4可以看出，数值模拟结果和试验结果相差不大，数值模拟得到的参数可以进行循环荷载作用下表征滑带土软化过程模拟.

1.3　基于Box-Jenkins建模的时间序列数学模型

时间序列的随机分析通常采用Box-Jenkins建模方法，使用Box-Jenkins方法建模的步骤如下（Pinzón-Hassan et al.， 2021； Ratshefola- Lerefolo and Mpeta， 2023）：

（1） 计算观测序列样本的相关系数和偏倚相关系数：当我们有一个样本序列x_t 时，我们可以计算样本的协方差如式（3）到式（7）所示.

\overline{\mu } = \frac{{\sum }_{t=\mathrm{1}}^n{x}_t}{n},(3)

s^{\mathrm{2}} = \frac{{\sum }_{t=\mathrm{1}}^n(x_t-\overline{\mu })(x_{(t)}-\overline{\mu })}{n-\mathrm{1}},(4)

\begin{array}{l}{\gamma }_k=\mathrm{1}/n{\sum }_{t=k+\mathrm{1}}^n(x_t-\overline{\mu })(x_{(t-k)}-\overline{\mu }),\\ {\rho }_{\mathrm{k}}=\frac{{\gamma }_k}{s^{\mathrm{2}}},\end{array} (5)

式中n是样本序列的个数， \overline{\mu }是平均值， s^{\mathrm{2}}是方差， {\gamma }_k是协方差， {\rho }_k是样本自相关函数.

X_t = {\mathrm{\Phi }}_{k\mathrm{1}} x_{t-\mathrm{1}} + {\mathrm{\Phi }}_{k\mathrm{2}} x_{t-\mathrm{2}} + \dots + {\mathrm{\Phi }}_{kk} x_{t-k} + u_t,(6)

自相关函数中的每个回归系数表示的是剔除中间变量x_{t- 1}、x_{t- 2}、…的影响后，x_t 与x_t-k 之间的自相关系数.x_{t-k+ 1}，即偏自相关函数.

（2） 模式识别：如果被检查的序列为白噪声序列，则无需调整，建模已完成.若序列为非平稳时间序列，自回归模型（ARIMA模型）被用来建模，反之则采用滑动平均模型（ARMA模型）来建模.

（3） 初始阶数和参数估计：在模式识别后，确定其对应数学模型的最高阶数；然后对模型进行从低阶到高阶的拟合和检验.有许多标准可用于建模｛εt｝.标准主要基于如下函数：

（4）

\delta \left(p^{\text{'}}\right) = n \mathrm{l} \mathrm{o} \mathrm{g} \left(\widehat{{\sigma }_{p^{\text{'}}}^{\mathrm{2}}}\right) + p^{\text{'}} g ( n )，（7）

式中，｛εt｝是数据x_t 的残差，是常阶函数， p '是常阶，是取p '时残差方差的估计，p '是g（n）= logn时的贝叶斯信息判据（BIC），该值越小，则预测模型越好.

（5） 拟合优度检验：对不同模型进行数学比较，并依据定阶法确定最适合的模型.

（6） 拟合检验：采用数学拟合和参数检验的方法，优秀的拟合模型应该能提取观测值序列中大多数样本的相关信息，即残差序列应该为白噪声序列.因此，模型的有效性检验即为残差序列的白噪声检验，通过检验进一步确定最合适的模型并采用最终预测模型对数据进行预测.

2 滑带土细观劣化分析

2.1　应变软化特性的配位数分析

统计不同循环荷载作用下，数值模型的P_C 值随循环加载次数的变化，如图5所示.在本研究中，我们认为5%的颗粒接触破碎率（即P_C =-5%）是颗粒运动所能达到的最大极限（Itasca， 2014）.当P_C 折减量大于5%时，假定有明显破损.仔细观察图5，配位数的变化出现波动主要原因为监测区域为固定不变的，而数值模型在进行循环荷载加载时，其数值模型中的颗粒之间的链接会出现一定的重分布，故由于荷载加载形式导致数据出现波动性变化.另一方面，当循环加载应力幅值为0.5 cm/s时，P_C 的初始变化不超过10%，这是由于这是一个很小的加载速率，只使数值土模型中的部分颗粒移动，即使随着循环加载的增加，P_C 的变化也不显著.同时，当循环加载应力幅值从0.9 cm/s增加到2.2 cm/s时，P_C 的降低率逐渐增加.为了突出循环加载下滑带土的应变软化程度，当循环加载次数为12 000次时，图6中的100%叠加柱状图显示了不同工况下滑带土的最终P_C 值，其值分别在对应分段中表示.例如，当滑带土在预定围压为50 kPa的情况下，循环加载应力幅值分别为1.5 cm/s和 2.2 cm/s时，可以看到滑带土的应变软化对应于图6中的P_C 值34.8和-53.0.另一方面，从图6可以明显看出，循环加载应力幅值是控制裂缝产生数量的主要因素.观察不同加载应力幅值对裂缝P_C 的影响，即图5a~5c中不同加载应力幅值下的P_C 值，不难发现滑带土裂缝在加载应力幅值为2.2 cm/s时，相对于其他循环加载应力幅值变化最大，说明循环加载应力幅值越大，循环加载次数越多，土体应变软化效应越明显.这说明循环加载应力幅值对断裂的发生（应变软化程度）具有决定性作用.

另一方面，随着滑带土体的预先围压的增大，即土体的围压从50 kPa增加到200 kPa（如图5a和图5c所示），不难发现，当加载应力幅值为0.9 cm/s时，围压分别为50 kPa和200 kPa时，P_C 值分别为-16.6和4.12，由上述讨论可知，当围压降低4倍时，裂缝增加了4倍.但观察图5a~5c时，加载应力幅值增加4倍时P_C 值的增加至少高出5倍.因此，滑带土的围压虽然可以抑制裂缝的增加，但对应变软化的抑制作用有限，循环加载应力幅值是土体应变软化最关键的因素.

2.2　滑带土应变软化特性的应力分析

C_n 的变化是循环荷载和预定围压作用下滑带土应变软化特性的指示参数，而滑带土应变软化特性的根本原因是剪切应力的变化.利用程序伺服控制在数值模拟中得到的滑动带土剪应力的变化趋势，讨论了不同预定围压和不同循环荷载次数下滑动带土剪应力的变化规律，如图7和图8所示，其中红色裂缝为拉裂缝，绿色裂缝为剪切裂缝.

从图7中可以看出，滑带土的峰值抗剪强度随着围压的增大而增大，而滑带土的应变软化程度随着循环荷载应力的增大和循环荷载次数的增加而逐渐增大.当循环加载应力幅值较低（循环加载强度小于或等于1.5 cm/s）时，滑带土的抗剪强度衰减以累积弱化效应为主.即随着循环加载的运行，滑带土的应变软化逐渐发展，剪切应力逐渐减小，而在较高的循环加载幅值（循环加载强度达到 2.2 cm/s）下，滑带土的力学参数弱化表现为附加荷载作用下的惯性破坏特征.在这种形式下，滑带土力学参数的弱化只有两个阶段.在第1阶段，随着循环加载的进行，滑带土有快速衰减的趋势，这意味着该阶段滑带土立即受到加载的破坏.在第二阶段，滑带土的强度参数已达到破坏极限范围，随着循环加载次数的增加基本没有明显变化.通过图8 a所示的 A1不难发现，上下剪切箱的交点处首先出现了一道裂缝，这是由于循环荷载作用造成的，上下剪切箱的左右两端同时以一定的速度运动.这是由于左上挡土墙和右上挡土墙的联合运动，使墙体带动附着在两侧挡土墙上的颗粒进行水平运动，从而间接将其转移到中心.因此，通过图8a所示的A1可以很容易地发现，裂缝的第一个位置是在上下剪盒的交界处.当循环加载次数增加到6 000次时，剪切盒结处出现的裂缝数量增加，同时，裂缝存在的位置从结处向数值模型中部扩展，这是由于数值模型中部颗粒的重排和移动造成的.

在循环荷载的作用下，数值模型粒子之间的相互作用使越来越多的粒子受到循环荷载的作用，它们相互运动引起的连锁反应呈几何级数增加.因此，当循环荷载次数为12 000时，滑带土数值模型中裂缝数量突然增加，模型内部出现透缝区，如图8a所示的C1.另一方面，循环加载应力幅值的增大也促进了滑带土裂缝的渗透.从图8a到图8d不难看出，循环加载应力幅值越大，滑带土的破坏行为越明显，裂缝增加越多，因此循环加载下滑带土的宏观破裂区范围越大.当循环加载强度的变化从0.5 cm/s到2.2 cm/s时，数值模型中裂纹的增长速率增长巨大，但最终均达到了平衡收敛的状态，这是由于地震波加载的方式导致了下剪切盒的墙体有一定的运动位移范围，从而导致振动波的加载对于断层带的影响范围有限，即无论振动强度的大小，数值模型中的断层带的颗粒小球，均随着振动波的加载，表现出一定的残余剪应力变化特征，将200 kPa的分析和扩展应用到 100 kPa和50 kPa.毫无疑问，预定围压越小，对裂缝生长速率的抑制作用越小，滑带土中裂缝的增加越明显，即滑带土应变软化程度越大.

3 滑坡稳定性的非线性预测

3.1　边坡稳定性系数的计算

选取有滑带土的边坡截面作为研究截面，如图9所示，其中假设边坡为层状结构，滑带土的厚度相对均匀，滑带土在不同位置的抗剪强度参数相同.为了合理分析滑带土抗剪强度与边坡稳定性之间的关系，选取极限平衡分析法中的Sarma法来计算含断层的边坡的稳定性，并进一步讨论和分析了不同循环加载次数和应力幅值下边坡稳定性的变化.滑动剖面包括滑动带土、卸载边界线和基岩.潜在滑坡体被划分为条状，潜在滑坡体被划分为8个垂直条状.为简化计算过程，1~8条强度参数采用均匀性值，其中各条带粘聚力为 466 kPa，摩擦角为29°，滑带土强度参数取上述2.2节数据.根据边坡稳定性分析结果，统计分析了不同加载次数和不同加载应力幅值下边坡稳定性系数的变化情况，其变化情况如图10所示.

随着应力幅值和循环加载次数的增加，边坡的稳定系数减小，如图9所示.由2.2节计算可知，在循环荷载作用下，滑带土体的应变逐渐变软.剪切应力的变化不仅是由土体围压引起的，而且循环加载应力对裂缝的扩展起关键作用.然而，有时剪应力的逐渐减小并不一定都是循环荷载的作用，如地下水的持续侵蚀和孔隙水压力的逐渐增大，会导致滑带土壤剪应力的逐渐减小.同时，仔细分析图8可以发现，裂缝的不断发育和延伸导致滑坡锚固截面不断缩小，边坡稳定系数不断减小，是滑坡发生的关键驱动因素之一.

因此，本文描述了两种可能的滑坡触发机制.对于模式1，通过分析图8a~ 8c中的A1和A2可知，随着循环加载次数的增加，剪切应力逐渐减小，裂缝也逐渐扩展为贯通断裂区，即随着循环加载，应变软化行为发生缓慢.图10中循环加载应力幅值小于等于1.5 cm/s时，边坡稳定系数变化也符合渐进破坏规律.另一方面，通过分析图8d中的A3和图10中的循环加载应力幅值等于2.2 cm/s时曲线的变化规律，可以很容易地得到模型2，可以看出，当剪应力不再随着循环加载次数的增加而增加时，裂缝的数量仍然在继续增加，表明滑坡已经触发.模式2表明，当循环加载应力较大时，随着循环加载次数的增加，滑带土剪切应力呈快速衰减趋势；因此，这种现象推断如果循环荷载的振幅和次数足够大，滑带土的应变速率会逐渐加速，并伴随一个大的裂隙贯通带，其次导致滑坡发生.

3.2　基于机器学习的时间序列分析在边坡稳定性预测中的应用

基于时间序列随机分析中的Box-Jenkins建模方法，结合上文中滑坡稳定性系数的计算结果，对该时间序列数据进行模式识别，确定最终建模数学模型，并进一步确定初始阶数，进而对模型进行从低阶到高阶的拟合，检验得到了最合适的数学模型，并采用最终预测模型对数据进行预测.

先将模型识别过程中的阶数设为1，然后建立ARIMA（1，1，1）（1，1，1）模型，通过式（3）~式（8）选择显示模型拟合测度，最终得到参数估价值，其中ARIMA（1，1，1）（1，1，1）模型中的自回归模型（AR模型）的参数估算值为0.019，数据统计概率为0.762，而ARIMA（1，1，1）（1，1，1）模型中的移动平均模型（MA模型）的参数估算值为0.714，数据统计概率为0.另一方面，ARIMA（1，1，1）（1，1，1）模型中的自回归模型的标准BIC值为10.16，对应的数据统计概率为0.748，而移动平均模型的标准BIC值为12.16，对应的数据统计概率为0.648.

AR（1）的参数估计为0.019，数据统计的概率为0.762，即接受AR（1）为零的原假设.MA（1）的参数估计为0.714，数据统计量的概率为0，否定了MA（1）为零的原始假设.因此，模型不是最优的，对数据的分析也不是很合适.因此，我们选择AR模型进行另一次调整，以得到最合适的模型.经过以上分析，可以看出数据是一个非平稳时间序列，因此建立了ARIMA（1，1，0）（0，1，1）模型的拟合优度，其中R平方调整后的拟合优度为0.747，标准BIC值为8.16，数据统计概率为0.848.

上述分析可以看出，标准BIC值为8.160，略小于ARIMA（1，1，1）（1，1，1）模型的标准BIC值.同时，数据统计的概率为0.848，因此，该模型是一个合适的模型，根据这个验证过的模型对图10中的数据进行建模和预测，如图11所示.由图11可以看出，对不同荷载情况下的边坡稳定系数进行了很好的预测，但两者的误差很小.因此该模型可用于边坡稳定系数的预测，为循环荷载作用下边坡的稳定性分析提供了一种准确、高效的学习方法.

4 结论

对于含有潜在滑带的岩质边坡，滑带土的不断劣化是滑坡发生的关键因素.定量研究滑带土中的裂隙分布与剪应力之间的关系，使得对于滑坡的触发机理的研究更有帮助.从这一角度出发，我们建立了一个三维滑带土体数值模型，旨在考虑滑带土应变软化导致滑坡发生的变化机制.同时提出了一种机器学习方法，可用于预测具有预定剪切面的边坡受循环载荷时的稳定性系数变化规律.

（1）滑带土剪切应力变化随着循环加载强度的增加，主要表现为累积渐进破坏和惯性破坏两种形式，不同的循环加载应力幅值导致边坡出现不同的破坏形式，累积渐进破坏含有时间延续性的特点，而惯性破坏主要发生在受到振动幅值较大的边坡.

（2）根据边坡稳定系数的变化与滑带土中裂缝的发育情况，推断出两种可能的滑坡触发机制，较大的循环加载强度和循环次数能加快土体软化的程度，边坡的稳定性变化和滑带土的软化存在非线性的对应关系，在实际边坡防护工程中，应加强滑带土强度变化的监测工作.

（3）通过机器学习，得到边坡稳定系数的数据是一个非平稳时间序列，通过建立Box-Jenkins模型的拟合优度，得到了验证过的预测模型，进而对边坡稳定系数数据进行建模和预测.解决了实际工程边坡稳定性难预测的关键科学问题，为循环荷载作用下的边坡稳定性分析提供了一种准确、高效的机器学习方法.

参考文献

[1] Alaei, E., Mahboubi, A., 2012. A Discrete Model for Simulating Shear Strength and Deformation Behaviour of Rockfill Material, Considering the Particle Breakage Phenomenon. Granular Matter, 14(6): 707-717.

[2] Balogun, A., Rezaie, F., Pham, Q. B., et al., 2021. Spatial Prediction of Landslide Susceptibility in Western Serbia Using Hybrid Support Vector Regression (SVR) with GWO, BAT and COA Algorithms. Geoscience Frontiers, 12(3): 384-398.

[3] Chen, W.X., 2021. Prediction of Roadbed Deformation of Qinghai-Tibet Railway Based on Time Series Model. Science Technology and Engineering, 21(35): 15203-15208 (in Chinese with English abstract).

[4] Chen, Z. R., Song, D. Q., Liu, X. L., et al., 2022. Seismic Dynamic Response Characteristics of a Layered Slope at Tunnel Entrance Using Shaking Table Test. Earth Science, 47(6): 2069-2080 (in Chinese with English abstract).

[5] Chigira, M., Yagi, H., 2006. Geological and Geomorphological Characteristics of Landslides Triggered by the 2004 Mid Niigta Prefecture Earthquake in Japan. Engineering Geology, 82: 202-221.

[6] Guo, Z., Chen, L., Gui, L., et al., 2019. Landslide Displacement Prediction Based on Variational Mode Decomposition and WA-GWO-BP Model. Landslides, 17(2):567-583.

[7] Guo, Z.Z., Tian, B.X.,Li, G.G., 2022. Using and Comparing Three Data-Driven Techniques to Generate Effective Regional Landslide Susceptibility Maps in the Loess Plateau of Northwest, China. Frontiers in Earth Science, 1-23.

[8] Hofmann, H., Babadagli, T., Zimmermann, G., 2015. A Grain Based Modeling Study of Fracture Branching during Compression Tests in Granites. Intertnational Journal of Rock Mechanics and Mining Sciences, 77: 152-162.

[9] Hu, F., Li, Z.Q., Hu, R. L., et al., 2018. Research on the Deformation Characteristics of Shear Band of Soil-Rock Mixture Based on Large Scale Direct Shear Test. Chinese Journal of Rock Mechanics and Engineering, 37(3): 766-778 (in Chinese with English abstract).

[10] Huang, F. M., Chen, J. W., Tang, Z. P., et al., 2021. Uncertainties of Landslide Susceptibility Prediction due to Different Spatial Resolutions and Different Proportions of Training and Testing Datasets. Chinese Journal of Rock Mechanics and Engineering, 40(6):1155-1169 (in Chinese with English abstract).

[11] Huang, F. M., Yin, K. L., Jiang, S. H., et al., 2018. Landslide Susceptibility Assessment Based on Clustering Analysis and Support Vector Machine. Chinese Journal of Rock Mechanics and Engineering, 37(1):156-167 (in Chinese with English abstract).

[12] Itasca, 2014. Version 5.0. Itasca Consulting Group Inc., Minneapolis.

[13] 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).

[14] Kirschbaum, D., Stanley, T., Zhou, Y., 2015. Spatial and Temporal Analysis of a Global Landslide Catalog. Geomorphology, 249: 4-15.

[15] Peethambaran, B., Anbalagan, R., Kanungo, D.P., et al., 2020. A Comparative Evaluation of Supervised Machine Learning Algorithms for Township Level Landslide Susceptibility Zonation in Parts of Indian Himalayas. Catena, 195:104751.

[16] Pinzón-Hassan, A.D., Tique-Ortiz, V., Zafra-Mejía, C.A., 2021. Box-Jenkins Stochastic Models for Studying Air Pollutants in a Latin American Megacity. Journal of Physics: Conference Series, 2139(1): 012003.

[17] Potyondy, D.O., Cundall, P.A., 2004. A Bonded-Particle Model for Rock. International Journal of Rock Mechanics and Mining Sciences, 41(8): 1329-1364.

[18] Ratshefola-Lerefolo, L., Mpeta, K.N., 2023. Analysis and Forecasting of Final Household Expenditure in South Africa Using Box-Jenkins-ARIMA Model. In: International Conference on Intelligent Computing & Optimization. Springer, Cham.

[19] Samui, P., 2008. Slope Stability Analysis: A Support Vector Machine Approach. Environmental Geology, 56: 255-261.

[20] Shan, P., Lai, X., 2019. Mesoscopic Structure PFC2D Model of Soil Rock Mixture Based on Digital Image. Journal of Visual Communication and Image Representation, 58:407-415.

[21] Sun, C., Tang, C. S., Cheng, Q., et al., 2022. Stability of Soil Slope under Soil-Atmosphere Interaction. Earth Science, 47(10): 3701-3722 (in Chinese with English abstract).

[22] Thai Pham, B., Tien Bui, D., Prakash, I., 2018. Landslide Susceptibility Modelling Using Different Advanced Decision Trees Methods. Civil Engineering and Environmental Systems, 35(1/2/3/4): 139-157. https://doi.org/10.1080/10286608.2019.1568418

[23] Vm, A., Mh, A., Zga, B., et al., 2021. Fast Physically-Based Model for Rainfall-Induced Landslide Susceptibility Assessment at Regional Scale-Science Direct. CATENA, 201:105213.

[24] Wang, C.，Dong,J.Y.，Liu, H.D.，et al.，2022. Shaking Table Model Test on Dynamic Response Characteristics and Failure Mechanism of Three Sections Locked Rock Slope.Earth Science，47(12): 4428-4441 (in Chinese with English abstract).

[25] Wang, F., Fan, X., Yunus, A.P., et al., 2019. Coseismic Landslides Triggered by the 2018 Hokkaido, Japan (Mw 6.6), Earthquake: Spatial Distribution, Controlling Factors, and Possible Failure Mechanism. Landslides, 16: 1551-1566.

[26] Xiao, T., Segoni, S., Liang, X., et al., 2022. Generating Soil Thickness Maps by Means of Geomorphological-Empirical Approach and Random Forest Algorithm in Wanzhou County, Three Gorges Reservoir. Geoscience Frontiers, 1-23.

[27] Yang, J., Yin, Z., Laouafa, F., et al., 2020. Hydromechanical Modeling of Granular Soils Considering Internal Erosion. Canadian Geotechnical Journal, 57(2): 157-172. https://doi.org/10.1139/cgj-2018-0653

[28] Yu, B., Pan, W. Z., Song, J., et al., 2012. Risk Zonation and Evaluation of Landslide Geohazards in Hangzhou. Rock and Soil Mechanics, 33(S1): 193-199, 216 (in Chinese with English abstract).

[29] Zhang, Y., Wong, L. N. Y., Meng, F., 2021. Brittle Fracturing in Low-Porosity Rock and Implications to Fault Nucleation. Engineering Geology, 285: 106025 .

[30] Zhang, Z., Zhang, X., Qiu, H., et al., 2016. Micro Parameters Sensitivity of Coarse Grained Material in Mechanics and Deformation Based on PFC3D. Materials Science Forum, 873:115-119.

[31] Zhang, Z.F., Huang, M., Tang, Z.C., 2022. Numerical Experimental Study on Shear Mechanical Properties of Particles in Rock Discontinuities. Chinese Journal of Geotechnical Engineering, 1-10 (in Chinese with English abstract).

[32] 陈卫雄, 2021. 基于时间序列模型的青藏铁路路基变形预测. 科学技术与工程, 21(35)：15203-15208.

[33] 陈志荣, 宋丹青, 刘晓丽, 等, 2022. 隧道口顺层斜坡地震动力响应特征振动台试验. 地球科学, 47(6): 2069-2080.

[34] 胡峰, 李志清, 胡瑞林, 等, 2018. 基于大型直剪试验的土石混合体剪切带变形特征试验研究. 岩石力学与工程学报, 37(3): 766-778.

[35] 黄发明, 陈佳武, 唐志鹏, 等, 2021. 不同空间分辨率和训练测试集比例下的滑坡易发性预测不确定性. 岩石力学与工程学报, 40(6): 1155-1169.

[36] 黄发明, 殷坤龙, 蒋水华, 等, 2018. 基于聚类分析和支持向量机的滑坡易发性评价. 岩石力学与工程学报, 37(1): 156-167.

[37] 蒋水华, 刘源, 章浩龙, 等, 2020. 先验概率分布及似然函数模型的选择对边坡可靠度评价影响的定量评估. 岩土力学, 41(9): 3087-3097.

[38] 孙畅,唐朝生,程青,等,2022. 土体‒大气相互作用下土质边坡稳定性研究.地球科学,47(10): 3701-3722.

[39] 王闯，董金玉，刘汉东，等，2022. 三段式锁固型岩质边坡动力响应特性及破坏机制振动台模型试验研究. 地球科学，47(12): 4428-4441.

[40] 俞布, 潘文卓, 宋健, 等, 2012. 杭州市滑坡地质灾害危险性区划与评价. 岩土力学, 33(S1): 193-199, 216.

[41] 张志飞,黄曼,唐志成, 2022. 岩石不连续面内颗粒剪切力学性质的数值试验研究. 岩土工程学报, 1-10.

基金资助

重庆市自然科学基金项目(cstc2021jcyjbsh0047)

自然资源部丘陵山地地质灾害防治重点实验室（福建省地质灾害重点实验室）开放基金资助项目(FJKLGHK2022K004)

浙江丽水地区灾害地质调查项目(DD20190648)

浙江飞云江流域地质灾害调查项目(DD20160282)

安徽理工大学引进人才基金项目(2022yjrc54)