Calibration of Micro Parameters of Particles in Granular Discrete Element Method to Assess Slope Stability by Strength Reduction Method
-
摘要: 根据岩土体力学指标标定颗粒细观参数是应用颗粒离散元解决岩土工程问题的一项基础工作。利用颗粒离散元法执行边坡稳定性强度折减法分析时，岩土体抗剪强度随折减系数变化不断调整，使用试算法标定颗粒细观参数效率严重不足。为解决此问题，采用国产颗粒离散元软件MatDEM，以颗粒细观参数为输入，以岩土体抗剪强度指标为输出，构建超定BP神经网络，开发双轴压缩数值模型并行试验技术，加速神经网络样本数据的获取过程，重复执行“逆向标定—精度检查—样本修正”以实现颗粒细观参数的逆向迭代修正标定。结果表明：该策略的精度明显优于直接以抗剪强度指标为输入、以颗粒细观参数为输出的欠定BP神经网络，标定的颗粒细观参数与抗剪强度目标高度匹配，数值试验结果与目标值相比误差可控制在不超过1%的范围内。采用澳大利亚计算机应用协会（ACADS）的均质土质边坡和堆石坝等两个边坡稳定性分析经典考题，以模型平均位移突变为极限状态判据执行强度折减法分析，检验该策略的应用效果。结果显示：该策略的标定能力可满足抗剪强度不断调整时颗粒细观参数重新设定的需要，安全系数计算结果与ACADS推荐解具有良好可比性。该策略可为MatDEM及其他颗粒离散元程序执行边坡稳定性强度折减法分析时颗粒细观参数标定提供可靠途径。Abstract: Calibration of micro parameters of particles based on the mechanical properties of rock and soil mass is essential to solve various geotechnical problems by using the granular discrete element method. If the strength reduction method is employed in the granular discrete element method to assess slope stability, the shear strength of rock and soil mass will be continuously adjusted according to the reduction factor, and then, micro parameters of particles should be constantly calibrated. In this case, the trial-and-error method is cumbersome and time-consuming. To solve the issue, take MatDEM as an example, an overdetermined BP neural network was built by setting micro parameters of particles as input and the shear strength of rock and soil mass as output. A technique was developed to carry out biaxial compression numerical tests in parallel. Then, a reverse-iterative-correct strategy was proposed to calibrate micro parameters of particles for a prescribed value of shear strength, by executing the “reversely determination—error check—sample correction” process repeatedly. Results of numerical tests verified that the new strategy apparently had a higher precision than the underdetermined BP neural network by taking the shear strength as input and micro parameters as output. Micro parameters obtained by the new strategy had a remarkable level of compliance with the specified shear strength, and the relative error between the prescribed values and the numerical test results was less than one percent. Two exam problems, a homogeneous soil slope and a rockfill dam, suggested by ACADS are adopted to verify the ability of the proposed strategy in slope stability analysis based on the strength reduction method. Results showed that the new strategy satisfied the calibration requirement on micro parameters of particles when executing the strength reduction method in MatDEM and the resulting safe factors agree with the safe factors recommended by ACADS. The proposed strategy provides a reliable approach for MatDEM and other granular discrete element software to calibrate micro parameters when assessing slope stability by the strength reduction method.
-
作为典型的非连续数值模拟方法，颗粒离散元法[1]可有效模拟岩土介质的开裂、分离等非连续变形行为，与以有限元为代表的连续数值模拟方法相比，更适用于边坡渐进变形破坏的过程模拟，并被成功应用于多处地质灾害体的失稳机理揭示，如2017年6月24日四川茂县新磨滑坡[2]、2018年5月5日四川乐乡市马边滑坡[3]等。周健等[4]将强度折减法引入颗粒离散元中，执行边坡稳定性评价，形成边坡稳定性分析的强度折减颗粒离散元法，在计算安全系数的同时可预测边坡破坏形态。Bao等[5]将该方法应用于某大型铁矿废渣堆场的稳定性评价，发现潜在破坏面形态与现场变形迹象吻合较好，证实了该方法的有效性。近年来，强度折减颗粒离散元法在临界破坏准则方面取得较好的研究进展，如：田雷等[6]绘制坡顶颗粒平均位移与折减系数关系曲线，以曲线位移突变为临界破坏；Lu等[7]从能量演化角度出发，认为边坡总体势能的突变可作为整体失稳判据；王培涛等[8]选取黄金分割点作为强度折减颗粒离散元法的折减比例，加速了边坡强度折减的计算过程。
通过上述研究，应用强度折减颗粒离散元法分析边坡稳定性的能力有所提升，但其在实际工程中应用仍较少。原因在于，颗粒离散元赋予颗粒一系列的细观参数模拟岩土体，应用强度折减颗粒离散元法评价边坡稳定性时，如何快速调整颗粒细观参数取值以对应岩土体抗剪强度折减的问题，仍有待解决[9]。针对给定的抗剪强度值，采用双轴压缩等数值试验模型，重复调整颗粒细观参数取值，使得数值试验测定的抗剪强度与给定值相符，是运用较广的颗粒离散元细观参数标定手段，称为试错法[10]，但其存在费时费力、盲目性大和精确度低等局限[11]。执行强度折减法时，抗剪强度将经历多次折减，相应地需要多次标定细观参数，试错法的能力无法满足这一需求。为解决此问题，部分研究将颗粒之间的法向、切向黏结强度和摩擦系数同步折减[12]，以表征抗剪强度的折减。由于颗粒细观参数与抗剪强度之间为复杂的非线性关系[13–14]，故该手段的合理性存疑。
BP神经网络能很好地映射输入参数与输出参数之间的非线性关系，在开采沉降预测[15]、煤体瓦斯渗透率预测[16]、岩体力学参数反演[17]、滑坡易发性评价[18]和地质灾害早期预警[19]等工程实践中的各类问题上得到成功应用。本文拟借助BP神经网络表达颗粒离散元细观参数与岩土体抗剪强度之间的数学关系，建立根据抗剪强度目标值快速标定颗粒细观参数的有效方法，以满足强度折减时细观参数快速调整的要求。
区别于已有颗粒离散元法研究多基于Itasca公司开发的颗粒流程序PFC的做法，本文采用南京大学自主研发的国产颗粒离散元软件MatDEM[20]。MatDEM可根据研究需要在Matlab中便捷地进行二次开发，自2018年5月发布以来受到国内外学者的高度关注，在边坡失稳模式预测[21]、高速远程滑坡运动过程分析[22]、泥石流成灾机理揭示[23]、光缆与土体的协调变形模拟[24]等问题上已取得较好的应用成果。新时代科技竞争日益激烈的背景下，国产科学计算软件开发和应用的重要性日益凸显。基于MatDEM研究强度折减颗粒离散元法分析边坡稳定性时颗粒细观参数的标定策略，可有效助力国产离散元软件的应用推广。
1. 颗粒细观参数和抗剪强度测定
1.1 MatDEM的颗粒细观参数
MatDEM通过胶接一系列具有特定力学性质的颗粒来构建岩土体模型，与PFC等其他颗粒离散元软件相同，其细观参数包括物理几何参数和接触模型参数两类。物理几何参数一般由具体模拟工程对象的基本特性决定[11]，包括数值模型尺寸、颗粒直径、颗粒集合孔隙率、颗粒密度等。其中，颗粒直径的取值一方面影响模型的宏观力学行为，另一方面关系到模型的计算规模和效率，因此必须慎重选取。宿辉等[25]统计研究后指出，当特征长度比L/R（模型长边尺寸L和颗粒平均半径R的比值）达到200以上时，模型的宏观力学性质与破坏模式将不受颗粒尺寸效应的影响。
MatDEM颗粒之间依靠弹簧相互接触产生力的作用，如图1所示。接触模型参数包括法向刚度Kn、切向刚度Ks、断裂位移Xb、切向初始抗剪力Fs0和摩擦系数μp。
图 1 MatDEM接触模型示意图Fig. 1 Schematic diagram of contact model in MatDEM
下载:
全尺寸图片
颗粒间法向力Fn与法向相对位移Xn关系为：
$$ {F}_{\text{n}}=\left\{ \begin{array}{l}{K}_{\text{n}}{X}_{\text{n}},{X}_{\text{n}} < 0; 受压状态;\\ {K}_{\text{n}}{X}_{\text{n}},{X}_{\text{n}} < {X}_{\text{b}}; 受拉状态\text{, }连接未断裂;\\ 0,{X}_{\text{n}} > 0; 连接断裂\text{, }无抗拉作用 \end{array}\right. $$ (1) 初始状态时颗粒与相邻颗粒相互连接，受法向弹簧的拉压作用。当两颗粒之间的Xn超过断裂位移Xb时，法向连接断裂，颗粒间的拉力消失，仅可承受压力作用。
颗粒间切向力Fs与切向相对位移Xs关系为：
$$ {F_{\text{s}}} = {K_{\text{s}}}{X_{\text{s}}} $$ (2) 切向连接的断裂满足摩尔库伦屈服准则：
$$ {F_{{\text{s,max}}}} = {F_{{\text{s0}}}} - {\mu _{\text{p}}}{F_{\text{n}}} $$ (3) 式中：Fs,max为抗剪力；切向初始抗剪力Fs0为没有施加法向压力时颗粒间的抗剪力，类似于岩土体的黏聚力。法向压力越大，抗剪力越大。当切向力Fs超过抗剪力时，切向连接断裂，颗粒间仅剩滑动摩擦力–μpFn；对应真实颗粒之间的胶结断裂，法向连接断裂时，切向连接也相应断开，反之亦然。
边坡颗粒离散元数值模型建立完毕后，模型物理几何参数将保持不变。应用强度折减法评价稳定性时，必须对接触模型参数进行调整，以对应宏观的岩土体抗剪强度折减。因此，本文细观参数标定主要针对接触模型参数进行。
1.2 基于双轴压缩模型的抗剪强度测定
颗粒离散元方法中，双轴压缩模型常用于测定颗粒胶结形成的岩土体抗剪强度。图2（a）所示的MatDEM双轴压缩模型数值试样的尺寸为150 mm×300 mm，由5 000个颗粒组成，直径为2.5～3.6 mm。模型的顶部和底部为刚性簇单元组成的刚性压力板，侧壁为柔性簇单元组成的柔性压力板，保持围压不变的情况下可轻微移动变形。先施加两侧压力板上的围压，再给顶部压力板逐步施加垂直向下位移，单步施加量0.05 mm，平衡后记录顶部压力板受力情况，直至试样最终破坏，如图2（b）所示。
图 2 双轴压缩试验示意图Fig. 2 Schematic diagram of biaxial compression numerical test下载:
全尺寸图片
以细观参数Kn、Ks、Xb、Fs0和μp分别取值2.9×107 N/mm、1.1×107 N/mm、4.2×10–4 m、6.4×107 Pa和0.12为例，在25、50和75 kPa这3个不同围压等级下进行双轴压缩试验，得到应力–应变曲线，如图3（a）所示。由图3（a）可见，轴向应力达到第1个峰值后，应力不变情况下，曲线存在明显应变增长段。以破坏时的轴压和围压为依据，做莫尔圆并绘制包络线，如图3（b）所示，可测定上述细观参数设定下形成的岩土体抗剪强度值为：内摩擦角ϕ为31.2°，黏聚力C为31.9 kPa。
图 3 基于双轴压缩试验的抗剪强度测定Fig. 3 Determination of the shear strength based on biaxial compression numerical test下载:
全尺寸图片
试错法重复调整颗粒细观参数使得双轴压缩试验测定的抗剪强度符合预期，但在效率和精度方面均存在局限。强度折减法评价边坡稳定性，岩土体抗剪强度值根据折减系数不断调整。试错法标定细观参数的能力不足，成为限制强度折减颗粒离散元法在边坡稳定性评价中广泛应用的技术瓶颈。
必须说明的是，真实岩土体的变形及强度指标往往具备较强的非线性，边坡不同部位的应力水平不同，其抗剪强度指标可能并不一致。采用数值模拟方法分析边坡稳定性时，引入随应力水平相应变化的抗剪强度无疑更为合理，这方面已有学者开展了探索性的工作[26]。考虑本文主要拟解决强度折减颗粒离散元法分析边坡稳定性时试错法标定细观参数能力不足的问题，在细观参数标定时，暂未考虑边坡不同部位应力水平差异对抗剪强度指标及其他参数的影响。
2. 基于BP神经网络的细观参数逆向迭代修正标定
2.1 细观参数标定的BP神经网络结构设计
BP神经网络是应用最广泛的神经网络框架之一，通过模拟大脑对信息处理的机制寻找输入信息与输出信息之间的非线性关系，由输入层、隐藏层、输出层3部分组成，如图4所示。
图 4 BP神经网络基本结构Fig. 4 Typical structure of BP neural network下载:
全尺寸图片
BP神经网络前向传播时神经元的输出值为：
$$ u_i^l = \sigma _i^l\left( {\sum\limits_{j = 1}^{{k_{l - 1}}} {w_{ji}^lu_j^{l - 1} + } b_i^l} \right) $$ (4) 式中：l为隐藏层层数；
$u_i^l$ 为第l层、第i个神经元的输出值；$\sigma _i^l$ 为第l层、第i个神经元的激活函数；${k_{l - 1}}$ 为第l–1层神经元个数；$w_{ji}^l$ 为计算第l层、第i个神经元时，第l–1层第j个神经元的权重；$b_i^l$ 为第l层、第i个神经元的偏置项。基于样本库训练完毕的BP神经网络，其前向传播过程事实上可视为关于输入信息x和输出信息y的非线性函数：$$ {\boldsymbol{y}} = f\left( {{\boldsymbol{x}},{\boldsymbol{w}},{\boldsymbol{b}}} \right) $$ (5) 式中，w和b分别代表所有权重和偏置信息的集合。
周瑜等[27]尝试使用BP神经网络解决颗粒离散元软件PFC的细观参数标定问题，神经网络构建时以弹性模量、单轴抗压强度和泊松比等3个力学指标为输入，以PFC颗粒的黏结半径系数等4个细观参数为输出，样本库为400组颗粒的细观参数和数值试验获取的岩土体力学指标。然而，其在执行“给定岩土体力学指标→神经网络输出颗粒细观参数→执行数值试验获取力学指标→获取值与给定值比较”的精度测试后发现， 10组测试方案中部分方案精度明显不符要求。原因在于，BP神经网络中以岩土体力学参数为输入信息x，以颗粒细观参数为输出信息y，此时输入信息维度少于输出信息，形成了数学上的欠定问题[28]；对于欠定问题，BP神经网络必须引入额外约束条件将其转化为适定或超定问题后才能保障求解精度[29]。
岩土体抗剪强度标定MatDEM颗粒细观参数时，已知信息为内摩擦角和黏聚力两个指标，待确定的颗粒细观参数为法向刚度Kn等5个指标。若按照常规方式将已知信息作为输入，以待定信息作为输出同样会导致欠定问题。为避免引入额外约束条件且保障神经网络前向传播数值精度，建立BP神经网络时，输入信息x为颗粒细观参数，输出信息y为内摩擦角和黏聚力，其结构如图5所示。神经网络的输入层神经元为5个，输出层神经元为两个；为保障数值精度，隐藏层设定为3层，每层有8个神经元。
图 5 细观参数和抗剪强度联系的神经网络结构Fig. 5 Neural network from micro parameters to shear strength下载:
全尺寸图片
BP神经网络的样本为MatDEM颗粒细观参数和双轴压缩试验获取的相应岩土体抗剪强度。参考细观参数的取值范围[20]，Fs0设置5个水平，Xb设置3个水平，μp、Kn和Ks统一设置4个水平，各水平取值见表1。全面正交组合共计形成960组细观参数进行双轴压缩数值试验。颗粒密度ρ取值1 900 kg/m3，颗粒直径取2.5～3.6 mm。
表 1 细观参数水平设定Table 1 Level setting of micro parameters变量 取值 Ks/(N·m–1) 1.0×103、1.0×104、1.0×105、1.0×106 Kn/(N·m–1) 1.0×103、5.0×104、1.0×106、1.0×107 Fs0/Pa 1.0×10–3、2.0×10–1、6.0×101、2.0×103、8.0×103 up 0、0.2、0.4、0.6 Xb/m 1.0×10–8、1.0×10–5、1.0×10–3 将960组MatDEM颗粒细观参数和双轴压缩数值试验获取的抗剪强度值作为神经网络样本数据库，执行BP神经网络的学习训练。首先，对样本数据进行归一化处理，避免各项数据量纲/量级的不一致对结果产生影响。归一化计算公式为：
$$ {X_{{\text{norm}}}} = 2 \times \frac{{X - {X_{\min }}}}{{{X_{\max }} - {X_{\min }}}} - 1 $$ (6) 式中，X、Xnorm分别为样本归一化之前和归一化之后的数据，Xmin和Xmax分别为数据最小值和最大值。
MatDEM是在Matlab框架下编制的，其可便捷地在Matlab中进行二次开发；同时，Matlab神经网络工具箱又为BP神经网络的学习训练提供了大量的预置函数，因此，所有程序部分可在Matlab中实现有效集成。BP神经网络学习训练时，训练函数选用trainlm函数，输出函数选用pureline函数，激活函数选择tansig函数，误差函数选用均方误差RMSE。
2.2 细观力学参数的逆向迭代修正标定
BP神经网络训练完毕后，式（5）中的所有权重和偏置信息均确定，其事实上也定义了黏聚力、内摩擦角与细观参数之间的复杂函数关系：
$$ C{\text{ = }}{f_1}\left( {{K_{\text{n}}},{K_{\text{s}}},{X_{\text{b}}},{F_{{\text{s0}}}},{\mu _{\text{p}}}} \right) $$ (7) $$ \varphi {\text{ = }}{f_2}\left( {{K_{\text{n}}},{K_{\text{s}}},{X_{\text{b}}},{F_{{\text{s0}}}},{\mu _{\text{p}}}} \right) $$ (8) 式中，f1和f2代表训练好的BP神经网络中由细观参数逐步运算至黏聚力和内摩擦角的过程，其直观表达式由于权重、偏置项和激活函数等的引入很难直接给出，但其真实存在。
在BP神经网络框架下，当给定输入信息时，可进行该输入信息的梯度计算，因此，仿照非线性方程数学求解，可基于式（7）和（8）进行给定抗剪强度目标值下细观参数取值的数学推演。首先，给予一组细观参数初值，使用前向传播获取抗剪强度值，将该值与目标值的误差定义为：
$$ R_{\text{MSE}} = {\left( {C - {C_0}} \right)^2} + {\left( {\varphi - {\varphi _0}} \right)^2} $$ (9) 式中，C0和ϕ0分别为黏聚力和内摩擦角的目标值。设定误差限为RMSE0，若误差不满足要求，沿误差梯度下降方向调整细观参数，直至误差满足要求，如图6所示。图6中：n为循环次数；lr为学习率，其引进目的为执行沿误差梯度下降方向的步长搜索逼近，本文取0.01。区别于周瑜等[27]的做法，本文的BP神经网络以颗粒细观参数为输入，以岩土体抗剪强度为输出，通过反推的方式获取岩土体抗剪强度给定后的颗粒细观参数，故称之为逆向标定。
图 6 BP神经网络逆向标定细观参数Fig. 6 Micro parameters determined by BP neural network下载:
全尺寸图片
对于逆向标定得到的细观参数，执行双轴压缩数值试验，检查抗剪强度与给定值的吻合程度，即精度检查。精度检查的结果可能并不理想，原因在于BP神经网络样本库质量未能达到要求。训练完毕的BP神经网络前向传播过程可视为双轴压缩数值试验的一种替代，样本库的样本数量越多，其替代的可靠程度越高。为保障效率，神经网络的样本库的样本数量有限，故前向传播过程替代双轴压缩数值试验的可靠度未必足够，导致逆向标定所得细观参数精度不足。
对于非线性方程求解而言，已训练完毕的BP神经网络虽然精度不足，但提供了粗略的搜索方向，因此，逆向标定获取的颗粒细观参数将比给定的细观参数初值更加趋近于最终解。鉴于此，当精度检查结果不理想时，将逆向标定得到的细观参数与其双轴压缩数值试验获取的抗剪强度作为一组新样本添加到已有样本库中，重新进行BP神经网络的学习训练，此过程可称为“样本修正”。BP神经网络经过重新学习训练后，式（7）和（8）中的函数f1和f2将相应地调整变化，以前次标定获得的细观参数为初值，再次进行逆向标定和精度检查。重复执行逆向标定—精度检查—样本修正的过程，直至精度检查顺利通过，则颗粒细观参数的标定过程完成，整个流程如图7所示。
图 7 细观参数标定的完整流程Fig. 7 Whole calibration process of micro parameters下载:
全尺寸图片
该流程中，每一轮的逆向标定均以上一轮标定获得的细观参数为初值，故颗粒细观参数的最终取值获取类似于数值求解中的迭代方式，因此，将这一总体流程称为颗粒细观参数的逆向迭代修正标定。
2.3 细观参数标定策略的有效性测试
为测试逆向迭代修正标定策略的有效性，给定抗剪强度目标值C=3.0 kPa、ϕ=19.6°，执行标定流程以获取颗粒细观参数。经过4轮逆向标定后，细观参数精度检查通过，每轮逆向标定获取的细观参数和双轴压缩数值试验获取的抗剪强度值见表2，抗剪强度与目标值的误差如图8所示。
表 2 每轮标定得到的细观参数及测定的抗剪强度Table 2 Micro parameters obtained in each round and its resulted shear strength逆向标
定轮数Kn/(107·
N·m–1)Ks/(106·
N·m–1)Xb/
(10–5 m)Fs0/
(105 Pa)μp C/kPa ϕ/(°) 1 0.27 5.27 1.95 1.25 0.220 2.17 30.2 2 0.47 1.84 2.72 3.67 0.120 2.84 22.8 3 2.03 2.32 1.84 6.42 0.043 2.95 18.9 4 1.27 2.09 2.04 5.95 0.041 3.02 19.6 
图 8 两种神经网络误差对比Fig. 8 Comparison on results of the two neural networks下载:
全尺寸图片
以抗剪强度给定值为标准，细观参数逆向迭代修正标定执行过程中，第1轮标定的黏聚力相对误差为27.52%，摩擦系数相对误差为59.70%；迭代1次后黏聚力相对误差降至4.19%，摩擦系数相对误差降至15.21%；4次迭代后，黏聚力相对误差降至0.70%，摩擦系数相对误差降至0.36%。因此，该方法在颗粒细观参数的标定上具有可靠的数值精度。同时，必须注意的是，前面2轮标定得到的细观参数精度不足。但是，逐轮的样本修正后，标定结果的数值精度将迅速提升。其原因在于，如第2.2节所述，每轮逆向标定获取的颗粒细观参数将比上轮标定获取的细观参数更趋近于最终解，那么，样本修正可理解为在最终解附近区域对神经网络不断加密。
为说明BP神经网络解决数学欠定问题时的局限性，图8补充绘制了欠定BP神经网络输出结果的误差情况。该欠定BP神经网络以抗剪强度为输入，以颗粒细观参数为输出，隐藏层层数和每层神经元的布置仍如图5所示。采用同样的数据样本库，BP神经网络训练完毕后，直接输入抗剪强度给定值得到颗粒细观参数。用欠定BP神经网络输出的颗粒细观参数进行双轴压缩模型试验，黏聚力的相对误差为44.67%，摩擦系数相对误差为45.85%。
3. 双轴压缩数值模型并行试验技术
BP神经网络标定细观参数的全过程中，数据样本库建立是一项基础而繁琐的工作。为加快样本库建立过程，开发双轴压缩数值模型并行试验技术，以充分发挥MatDEM中GPU计算模式的能力[30]。MatDEM中有GPU和CPU两种计算模式，执行不同功能时调用不同的计算模式。一般来说，数值模型建立和材料参数赋值为CPU计算模式。受力后的颗粒平衡迭代计算则会根据颗粒数量在GPU和CPU两种模式中自动切换，单元数量较少时选择CPU模式，单元数量增加到一定程度则选择GPU模式。单个双轴压缩数值模型试验盒中，颗粒数量一般在5 000左右，无法发挥GPU模式计算优势。
本文对MatDEM进行二次开发，实现双轴压缩数值模型的并行试验；同时，构建多个双轴压缩模型，每个模型的细观参数单独赋值，压力板位移同步施加，颗粒受力后的平衡迭代同步进行，压力板数据同步记录。此时，参与同步计算的颗粒数量将达数十万，利用专业工作站的高性能GPU，可有效提高计算效率。图9为54个双轴压缩数值模型的排列分布情况和并行试验获取的试样破坏时的位移云图。
图 9 54个双轴压缩模型的排列分布和位移云图Fig. 9 Distribution and displacement diagrams of 54 biaxial compression models下载:
全尺寸图片
试验研究发现，双轴压缩数值模型并行试验所需要的时间一方面取决于GPU的计算能力，另一方面取决于双轴压缩数值模型的排列分布。当双轴压缩数值模型的排列分布呈现图9所示的正方形布局时，并行试验的效率最高。其原因在于，以颗粒离散元为理论基础的MatDEM，执行逐时间步的增量计算，必然涉及动力系统的能量平衡和耗散问题。默认阻尼下，MatDEM数值模型完成一次标准平衡，即应力波从模型一端传至另一端所需的时间Ttotal为：
$$ {T_{{\text{total}}}} = \Delta T \times {N_{{\text{pack}}}} $$ (10) 式中：Npack表示模型最长维度上堆积的颗粒数；
$ \Delta T $ 为动力计算采用的时间增量步步长，取值需远小于弹簧简谐振动周期，以精确模拟数值模型弹性变形过程。多个双轴压缩数值模型并行试验时，MatDEM会将同步试验的数值模型视为整体。双轴压缩数值模型总个数一定时，采用正方形布局时，整体系统在两个维度上的颗粒数基本接近，效率最高。硬件配置Intel（R） Xeon（R） E-2224@3.40 GHz CPU，NIVIDIA Quadro RTX 4000显卡的情况下，单个双轴压缩模型的试验耗时约50 min。同样配置条件下，54个双轴压缩模型并行试验总耗时约480 min，节约82%的计算时间。因此，采用并行试验技术进行双轴压缩数值模型试验，可充分发挥GPU计算模式的能力。
4. 边坡稳定性算例分析
以澳大利亚计算机应用协会（简称ACADS）的两个边坡经典考题为例，采用MatDEM执行边坡稳定性的强度折减法分析，考察本文提出的颗粒细观参数标定策略的实用性。稳定性分析时，临界破坏准则参照田雷等[6]的做法，记边坡模型中所有颗粒的位移平均值为模型平均位移S，以其发生突变为边坡达到极限状态的判据。
4.1 算例1：均质土质边坡
图10为均质土质边坡，土体力学参数为C=3.0 kPa，ϕ=19.6°，重度γ=20.0 kN/m3。ACADS给出的安全系数推荐值为1.00，对应的滑动范围如图10橙色区域所示。
图 10 算例1模型及ACADS推荐解的滑动范围Fig. 10 Model of Example 1 and the failure region by ACADS下载:
全尺寸图片
在MatDEM中建立该边坡的颗粒离散元模型如图11所示，模型含82 698个直径67～96 mm的颗粒。此时，根据特征长度比L/R要求[25]将用于测定抗剪强度的双轴压缩数值模型尺寸调整为4 000 mm×8 000 mm，以避免颗粒尺寸效应影响细观参数标定。
图 11 边坡的MatDEM模型和临界失稳状态位移Fig. 11 MatDEM model and its critical displacements下载:
全尺寸图片
分析稳定性时，初始折减系数设为0.70；若边坡处于稳定状态，则将折减系数增加0.10继续计算，折减系数达到1.10时边坡失稳破坏。在折减系数0.95到1.10范围内进一步细化分析，最终确定边坡稳定系数为1.02。边坡模型平均位移S随着强度折减系数Fsr增大的变化曲线如图12所示。由图12可知：Fsr达到1.02前，S基本呈水平直线；数据检查发现，S随着Fsr增加仍略有增加，表明抗剪强度下降会影响颗粒平衡后的位置，宏观上体现为变形轻微增长；Fsr达到1.02后，S值随Fsr增加而急剧变大，表明边坡进入临界失稳状态。边坡达到临界失稳状态时的位移，如图11所示。此时边坡自坡脚到坡顶已出现明显贯通的滑动面，边坡沿该滑动面将开始发生滑动。
图 12 算例1的S–Fsr曲线Fig. 12 S–Fsr curves for Example 1下载:
全尺寸图片
强度折减过程中，抗剪强度不同折减系数折减后对应的细观参数见表3。由表3可知，由于抗剪强度与细观参数之间的非线性映射关系非常复杂，抗剪强度的折减很难通过直接对细观参数进行折减来实现。
表 3 不同折减系数下土体抗剪强度与细观力学参数Table 3 Micro parameters and soil shear strength corresponding to reduction coefficients折减
系数C/kPa ϕ/(°) Kn/(107·
N·m–1)Ks/(106·
N·m–1)Xb/
(10–5 m)Fs0/
(105 Pa)μp 0.70 4.28 28.0 4.49 7.74 4.46 22.40 0.072 0.80 3.75 24.5 2.56 5.76 3.47 17.20 0.063 0.90 3.33 21.7 1.42 4.81 2.82 12.40 0.054 0.95 3.15 20.6 3.36 4.22 2.46 9.15 0.048 0.97 3.09 20.2 1.21 2.09 2.37 7.27 0.043 0.99 3.03 19.8 1.25 2.09 2.16 6.68 0.042 1.00 3.00 19.6 1.27 2.09 2.04 5.95 0.040 1.02 2.94 19.2 3.68 7.40 1.91 5.07 0.038 1.03 2.91 19.2 4.02 7.20 1.85 4.74 0.036 1.04 2.88 18.8 3.96 7.52 1.72 3.92 0.033 1.05 2.85 18.6 4.57 6.89 1.61 3.47 0.031 4.2 算例2：堆石坝模型
算例2为一由堆石体、过渡层、心墙和反滤层组成的堆石坝，坝高146 m，坝底宽648 m，其形态如图13所示。
图 13 堆石坝模型及ACADS推荐解的滑动区域Fig. 13 A rock fill dam and the failure region by ACADS下载:
全尺寸图片
坝体模型各种材料的力学参数见表4，由于算例2是用于稳定性分析，其参数与实际工程中各种材料的参数取值存在一定差异。ACADS对该堆石坝模型给出的安全系数推荐值为1.976，对应的滑动面如图13所示。
表 4 堆石坝模型各材料力学参数Table 4 Mechanical properties of materials in rockfill dam model材料 C/kPa ϕ/(°) r/(kN·m–3) 堆石 0 45 20.4 过渡层 0 45 20.4 反滤层 0 45 20.4 心墙 8.5 23 20.4 在MatDEM中建立该模型的颗粒离散元模型如图14所示，模型含89 768个直径为0.33～0.48 m的颗粒。颗粒细观参数标定时，按照特征长度比L/R的要求[25]，将双轴压缩数值模型尺寸调整为20 m×40 m。
图 14 堆石坝的MatDEM模型Fig. 14 MatDEM model for the rock fill dam下载:
全尺寸图片
采用强度折减法分析该堆石坝模型的稳定性，模型的平均位移S和强度折减系数Fsr关系曲线如图15（a）所示。以S突变为模型临界失稳判据，确定该堆石坝安全系数为1.975。模型达到临界失稳状态时的位移如图15（b）所示。图15（a）显示，即便模型处于稳定状态，模型平均位移S已超过2 m。其原因在于，颗粒离散元边坡模型建立完毕之后，由于存在重力势能，颗粒自平衡过程中会导致模型产生一定的变形，类似于重力作用下的沉降变形。事实上，算例1的边坡模型处于稳定状态时，也存在0.002 m的模型平均位移S。算例2模型平均位移S偏大，其原因主要为，施加于模型上的边界约束条件，堆石坝模型仅在底部边界设置位移约束，而算例1边坡模型在底部和左右两侧边界均设置位移约束，因此堆石坝模型重力效应可充分发挥，变形相对较大；此外，堆石坝模型高度达146 m，大尺寸规模下，颗粒自平衡过程中的模型变形也相对较大。
图 15 堆石坝的安全系数结果和临界失稳状态位移Fig. 15 Safety factors and the critical displacements for the rock fill dam下载:
全尺寸图片
强度折减过程中，材料颗粒细观参数随着折减系数变化而调整的情况见表5。
表 5 不同折减系数下堆石坝模型中材料的抗剪强度与细观力学参数Table 5 Micro parameters and shear strength of material in rockfill dam model corresponding to reduction coefficients折减
系数堆石、过渡层、反滤层 心墙 C/kPa ϕ/(°) Kn/(108·
N·m–1)Ks/(108·
N·m–1)Xb/(10–4 m) Fs0/Pa μp C/kPa ϕ/(°) Kn/(108·
N·m–1)Ks/(108·
N·m–1)Xb/(10–2 m) Fs0/(106 Pa) μp 0.70 0 64.3 4.76 2.94 1.47 152.1 0.68 12.14 32.86 11.63 4.32 2.74 9.86 0.336 0.90 0 50.0 2.52 1.94 1.36 147.0 0.66 9.44 25.56 8.74 4.65 2.13 7.67 0.261 1.00 0 45.0 1.41 2.33 1.27 128.0 0.64 8.50 23.00 4.47 3.49 1.92 6.90 0.235 1.10 0 40.9 1.74 3.07 1.30 132.0 0.65 7.73 20.91 4.47 3.49 1.75 6.27 0.214 1.30 0 34.6 1.41 2.33 1.26 126.0 0.63 6.54 17.69 4.95 1.79 1.48 5.31 0.181 1.50 0 30.0 1.57 2.33 1.23 117.0 0.62 5.67 15.33 9.29 1.98 1.28 4.60 0.157 1.70 0 26.5 1.53 1.86 1.21 104.0 0.61 5.00 13.53 6.05 3.71 1.13 4.06 0.138 1.90 0 23.7 1.53 1.86 1.17 92.0 0.60 4.47 12.11 6.05 3.71 1.01 3.63 0.124 1.92 0 23.4 1.96 0.99 1.02 86.0 0.58 4.43 11.98 3.51 2.42 1.00 3.59 0.123 1.94 0 23.2 1.20 0.68 0.99 80.2 0.57 4.38 11.86 5.99 1.40 0.99 3.56 0.121 1.95 0 23.1 1.76 0.55 0.98 70.1 0.58 4.36 11.79 11.21 2.39 0.99 3.54 0.121 1.96 0 23.0 2.31 0.61 0.97 67.2 0.59 4.34 11.73 6.88 2.75 0.98 3.52 0.120 1.97 0 22.8 1.27 0.60 0.97 61.5 0.57 4.31 11.68 11.08 4.24 0.98 3.50 0.119 1.98 0 22.7 1.66 0.77 0.96 55.2 0.55 4.29 11.62 4.61 4.32 0.97 3.48 0.119 两个算例的执行情况显示，本文提出的细观参数标定策略，可满足MatDEM等颗粒离散元软件采用强度折减法分析边坡稳定性时，迅速调整颗粒细观参数取值，以对应岩土体抗剪强度折减的需求。安全系数结果方面：算例1安全系数为1.02，与ACADS推荐的安全系数相差0.02；算例2安全系数为1.975，与ACADS推荐的系数相差0.001。滑动面结果方面：MatDEM获取的临界失稳状态下模型的滑动区域与ACADS推荐的滑动范围基本吻合。总体上，MatDEM获取的安全系数和滑动范围与ACADS推荐解之间具有良好的可比性，表明采用本文策略标定的颗粒细观参数，在对应不同程度折减后的岩土体抗剪强度方面具有满意的数值精度。因此，在应用颗粒离散元强度折减法评价边坡稳定性时，本文提出的颗粒细观参数逆向迭代修正标定策略可发挥积极作用。
5. 结 论
针对以MatDEM为代表的颗粒离散元法软件应用强度折减技术评价边坡稳定性时的颗粒细观参数标定问题，以颗粒细观参数为输入，抗剪强度为输出构造BP神经网络，建立颗粒细观参数的逆向迭代修正标定方法，开发双轴压缩数值模型并行试验技术以加快神经网络样本库样本的获取，并采用ACADS的两个边坡稳定性分析经典考题验证了颗粒细观参数逆向迭代修正标定方法的能力和精度。基于研究过程中的发现，可得结论如下：
1）细观参数逆向迭代修正标定方法可有效解决颗粒离散元强度折减法分析边坡稳定性时的颗粒细观参数标定问题。标定得到的细观参数可良好对应岩土体的抗剪强度，数值算例计算结果与推荐解之间的可比性好。
2）MatDEM的二次开发在Matlab程序实现便捷，使得MatDEM具备很好的研究潜力和应用价值。本文开发的双轴压缩数值模型并行试验技术有效加速了神经网络样本库的建立过程，可为其他小型模型的重复试验工作提供参考。
3）使用BP神经网络执行参数反演问题研究时，必须注意神经网络输入信息和输出信息维度的大小关系，以避免形成数学上的欠定问题。
-
图 1 MatDEM接触模型示意图
Fig. 1 Schematic diagram of contact model in MatDEM
下载:
全尺寸图片
图 2 双轴压缩试验示意图
Fig. 2 Schematic diagram of biaxial compression numerical test
下载:
全尺寸图片
图 3 基于双轴压缩试验的抗剪强度测定
Fig. 3 Determination of the shear strength based on biaxial compression numerical test
下载:
全尺寸图片
图 4 BP神经网络基本结构
Fig. 4 Typical structure of BP neural network
下载:
全尺寸图片
图 5 细观参数和抗剪强度联系的神经网络结构
Fig. 5 Neural network from micro parameters to shear strength
下载:
全尺寸图片
图 6 BP神经网络逆向标定细观参数
Fig. 6 Micro parameters determined by BP neural network
下载:
全尺寸图片
图 7 细观参数标定的完整流程
Fig. 7 Whole calibration process of micro parameters
下载:
全尺寸图片
图 8 两种神经网络误差对比
Fig. 8 Comparison on results of the two neural networks
下载:
全尺寸图片
图 9 54个双轴压缩模型的排列分布和位移云图
Fig. 9 Distribution and displacement diagrams of 54 biaxial compression models
下载:
全尺寸图片
图 10 算例1模型及ACADS推荐解的滑动范围
Fig. 10 Model of Example 1 and the failure region by ACADS
下载:
全尺寸图片
图 11 边坡的MatDEM模型和临界失稳状态位移
Fig. 11 MatDEM model and its critical displacements
下载:
全尺寸图片
图 12 算例1的S–Fsr曲线
Fig. 12 S–Fsr curves for Example 1
下载:
全尺寸图片
图 13 堆石坝模型及ACADS推荐解的滑动区域
Fig. 13 A rock fill dam and the failure region by ACADS
下载:
全尺寸图片
图 14 堆石坝的MatDEM模型
Fig. 14 MatDEM model for the rock fill dam
下载:
全尺寸图片
图 15 堆石坝的安全系数结果和临界失稳状态位移
Fig. 15 Safety factors and the critical displacements for the rock fill dam
下载:
全尺寸图片
表 1 细观参数水平设定
Table 1 Level setting of micro parameters
变量 取值 Ks/(N·m–1) 1.0×103、1.0×104、1.0×105、1.0×106 Kn/(N·m–1) 1.0×103、5.0×104、1.0×106、1.0×107 Fs0/Pa 1.0×10–3、2.0×10–1、6.0×101、2.0×103、8.0×103 up 0、0.2、0.4、0.6 Xb/m 1.0×10–8、1.0×10–5、1.0×10–3 表 2 每轮标定得到的细观参数及测定的抗剪强度
Table 2 Micro parameters obtained in each round and its resulted shear strength
逆向标
定轮数Kn/(107·
N·m–1)Ks/(106·
N·m–1)Xb/
(10–5 m)Fs0/
(105 Pa)μp C/kPa ϕ/(°) 1 0.27 5.27 1.95 1.25 0.220 2.17 30.2 2 0.47 1.84 2.72 3.67 0.120 2.84 22.8 3 2.03 2.32 1.84 6.42 0.043 2.95 18.9 4 1.27 2.09 2.04 5.95 0.041 3.02 19.6 表 3 不同折减系数下土体抗剪强度与细观力学参数
Table 3 Micro parameters and soil shear strength corresponding to reduction coefficients
折减
系数C/kPa ϕ/(°) Kn/(107·
N·m–1)Ks/(106·
N·m–1)Xb/
(10–5 m)Fs0/
(105 Pa)μp 0.70 4.28 28.0 4.49 7.74 4.46 22.40 0.072 0.80 3.75 24.5 2.56 5.76 3.47 17.20 0.063 0.90 3.33 21.7 1.42 4.81 2.82 12.40 0.054 0.95 3.15 20.6 3.36 4.22 2.46 9.15 0.048 0.97 3.09 20.2 1.21 2.09 2.37 7.27 0.043 0.99 3.03 19.8 1.25 2.09 2.16 6.68 0.042 1.00 3.00 19.6 1.27 2.09 2.04 5.95 0.040 1.02 2.94 19.2 3.68 7.40 1.91 5.07 0.038 1.03 2.91 19.2 4.02 7.20 1.85 4.74 0.036 1.04 2.88 18.8 3.96 7.52 1.72 3.92 0.033 1.05 2.85 18.6 4.57 6.89 1.61 3.47 0.031 表 4 堆石坝模型各材料力学参数
Table 4 Mechanical properties of materials in rockfill dam model
材料 C/kPa ϕ/(°) r/(kN·m–3) 堆石 0 45 20.4 过渡层 0 45 20.4 反滤层 0 45 20.4 心墙 8.5 23 20.4 表 5 不同折减系数下堆石坝模型中材料的抗剪强度与细观力学参数
Table 5 Micro parameters and shear strength of material in rockfill dam model corresponding to reduction coefficients
折减
系数堆石、过渡层、反滤层 心墙 C/kPa ϕ/(°) Kn/(108·
N·m–1)Ks/(108·
N·m–1)Xb/(10–4 m) Fs0/Pa μp C/kPa ϕ/(°) Kn/(108·
N·m–1)Ks/(108·
N·m–1)Xb/(10–2 m) Fs0/(106 Pa) μp 0.70 0 64.3 4.76 2.94 1.47 152.1 0.68 12.14 32.86 11.63 4.32 2.74 9.86 0.336 0.90 0 50.0 2.52 1.94 1.36 147.0 0.66 9.44 25.56 8.74 4.65 2.13 7.67 0.261 1.00 0 45.0 1.41 2.33 1.27 128.0 0.64 8.50 23.00 4.47 3.49 1.92 6.90 0.235 1.10 0 40.9 1.74 3.07 1.30 132.0 0.65 7.73 20.91 4.47 3.49 1.75 6.27 0.214 1.30 0 34.6 1.41 2.33 1.26 126.0 0.63 6.54 17.69 4.95 1.79 1.48 5.31 0.181 1.50 0 30.0 1.57 2.33 1.23 117.0 0.62 5.67 15.33 9.29 1.98 1.28 4.60 0.157 1.70 0 26.5 1.53 1.86 1.21 104.0 0.61 5.00 13.53 6.05 3.71 1.13 4.06 0.138 1.90 0 23.7 1.53 1.86 1.17 92.0 0.60 4.47 12.11 6.05 3.71 1.01 3.63 0.124 1.92 0 23.4 1.96 0.99 1.02 86.0 0.58 4.43 11.98 3.51 2.42 1.00 3.59 0.123 1.94 0 23.2 1.20 0.68 0.99 80.2 0.57 4.38 11.86 5.99 1.40 0.99 3.56 0.121 1.95 0 23.1 1.76 0.55 0.98 70.1 0.58 4.36 11.79 11.21 2.39 0.99 3.54 0.121 1.96 0 23.0 2.31 0.61 0.97 67.2 0.59 4.34 11.73 6.88 2.75 0.98 3.52 0.120 1.97 0 22.8 1.27 0.60 0.97 61.5 0.57 4.31 11.68 11.08 4.24 0.98 3.50 0.119 1.98 0 22.7 1.66 0.77 0.96 55.2 0.55 4.29 11.62 4.61 4.32 0.97 3.48 0.119 -
[1] Cundall P A,Strack O D L.A discrete numerical model for granular assemblies[J].Géotechnique,1979,29(1):47–56. doi: 10.1680/geot.1980.30.3.331 [2] Fan Xuanmei,Xu Qiang,Scaringi G,et al.Failure mechanism and kinematics of the deadly June 24th 2017 Xinmo landslide,Maoxian,Sichuan,China[J].Landslides,2017,14(6):2129–2146. doi: 10.1007/s10346-017-0907-7 [3] Wei Jiangbo,Zhao Zhou,Xu Chong,et al.Numerical investigation of landslide kinetics for the recent Mabian landslide(Sichuan,China)[J].Landslides,2019,16(11):2287–2298. doi: 10.1007/s10346-019-01237-0 [4] 周健,王家全,曾远,等.颗粒流强度折减法和重力增加法的边坡安全系数研究[J].岩土力学,2009,30(6):1549–1554. doi: 10.3969/j.issn.1000-7598.2009.06.003 Zhou Jian,Wang Jiaquan,Zeng Yuan,et al.Slope safety factor by methods of particle flow code strength reduction and gravity increase[J].Rock and Mechanics,2009,30(6):1549–1554 doi: 10.3969/j.issn.1000-7598.2009.06.003 [5] Bao Yiding,Sun Xiaohui,Chen Jianping,et al.Stability assessment and dynamic analysis of a large iron mine waste dump in Panzhihua,Sichuan,China[J].Environmental Earth Sciences,2019,78(2):1–17. doi: 10.1007/s12665-019-8043-4 [6] 田雷,杨永森,李敬军,等.基于离散元法植物根系对黏性土坡稳定性的影响[J].中国农业大学学报,2018,23(7):133–140. doi: 10.11841/j.issn.1007-4333.2018.07.15 Tian Lei,Yang Yongsen,Li Jingjun,et al.Influence of plant roots on cohesive soil slope stability based on discrete element method[J].Journal of China Agricultural University,2018,23(7):133–140 doi: 10.11841/j.issn.1007-4333.2018.07.15 [7] Lu Ye,Tan Yong,Li Xiang.Stability analyses on slopes of clay-rock mixtures using discrete element method[J].Engineering Geology,2018,244:116–124. doi: 10.1016/j.enggeo.2018.07.021 [8] 王培涛,杨天鸿,朱立凯,等.基于PFC2D岩质边坡稳定性分析的强度折减法[J].东北大学学报(自然科学版),2013,34(1):127–130. doi: 10.3969/j.issn.1005-3026.2013.01.030 Wang Peitao,Yang Tianhong,Zhu Likai,et al.Strength reduction method for rock slope stability analysis based on PFC2D[J].Journal of Northeastern University(Natural Science),2013,34(1):127–130 doi: 10.3969/j.issn.1005-3026.2013.01.030 [9] Cheng Yi,Wong L N Y.A study on mechanical properties and fracturing behavior of Carrara marble with the flat-jointed model[J].International Journal for Numerical and Analytical Methods in Geomechanics,2020,44(6):803–822. doi: 10.1002/nag.3040 [10] 张社荣,孙博,王超,等.双轴压缩试验下岩石裂纹扩展的离散元分析[J].岩石力学与工程学报,2013,32(2):3083–3091. Zhang Sherong,Sun Bo,Wang Chao,et al.Discrete element analysis of crack propagation in rocks under biaxial compression[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(2):3083–3091 [11] 李坤蒙,李元辉,徐帅,等.PFC2D数值计算模型微观参数确定方法[J].东北大学学报(自然科学版),2016,37(4):563–567. doi: 10.3969/j.issn.1005-3026.2016.04.023 Li Kunmeng,Li Yuanhui,Xu Shuai,et al.Method to determine microscopic parameters of PFC2D numerical model[J].Journal of Northeastern University(Natural Science),2016,37(4):563–567 doi: 10.3969/j.issn.1005-3026.2016.04.023 [12] 戴健健,朱振学,叶楠,等.基于位移突变准则的颗粒流强度折减法[J].科学技术与工程,2021,21(28):12052–12057. doi: 10.3969/j.issn.1671-1815.2021.28.018 Dai Jianjian,Zhu Zhenxue,Ye Nan,et al.Particle flow strength reduction method based on displacement mutation criterion[J].Science Technology and Engineering,2021,21(28):12052–12057 doi: 10.3969/j.issn.1671-1815.2021.28.018 [13] Ye Fangping,Wheeler C,Chen Bin,et al.Calibration and verification of DEM parameters for dynamic particle flow conditions using a backpropagation neural network[J].Advanced Powder Technology,2019,30(2):292–301. doi: 10.1016/j.apt.2018.11.005 [14] Sun Wei,Wu Shunchuan,Cheng Ziqiao,et al.Interaction effects and an optimization study of the microparameters of the flat-joint model using the Plackett-Burman design and response surface methodology[J].Arabian Journal of Geosciences,2020,13:1–24. doi: 10.1007/s12517-019-5018-4 [15] Freitag S,Cao B T,Ninić J,et al.Recurrent neural networks and proper orthogonal decomposition with interval data for real-time predictions of mechanised tunnelling processes[J].Computers & Structures,2018,207:258–273. doi: 10.1016/j.compstruc.2017.03.020 [16] 尹光志,李铭辉,李文璞,等.基于改进BP神经网络的煤体瓦斯渗透率预测模型[J].煤炭学报,2013,38(7):1179–1184. doi: 10.13225/j.cnki.jccs.2013.07.026 Yin Guangzhi,Li Minghui,Li Wenpu,et al.Model of coal gas permeability prediction based on improved BP neural network[J].Journal of China Coal Society,2013,38(7):1179–1184 doi: 10.13225/j.cnki.jccs.2013.07.026 [17] 金长宇,马震岳,张运良,等.神经网络在岩体力学参数和地应力场反演中的应用[J].岩土力学,2006,27(8):1263–1266. doi: 10.3969/j.issn.1000-7598.2006.08.007 Jin Changyu,Ma Zhenyue,Zhang Yunliang,et al.Application of neural network to back analysis of mechanical parameters and initial stress field of rock masses[J].Rock and Soil Mechanics,2006,27(8):1263–1266 doi: 10.3969/j.issn.1000-7598.2006.08.007 [18] Van Dao D,Jaafari A,Bayat M,et al.A spatially explicit deep learning neural network model for the prediction of landslide susceptibility[J].CATENA,2020,188:104451. doi: 10.1016/j.catena.2019.104451 [19] Niu Hongtao.Smart safety early warning model of landslide geological hazard based on BP neural network[J].Safety Science,2020,123:104572. doi: 10.1016/j.ssci.2019.104572 [20] 刘春.地质与岩土工程矩阵离散元分析[M].北京:科学出版社,2019. [21] Zhan Qinghua,Wang Shimei,Wang Li,et al.Analysis of failure models and deformation evolution process of geological hazards in Ganzhou city,China[J].Frontiers in Earth Science,2021,9:731447. doi: 10.3389/feart.2021.731447 [22] Xia Guoqing,Liu Chun,Xu Chong,et al.Dynamic analysis of the high-speed and long-runout landslide movement process based on the discrete element method:A case study of the Shuicheng landslide in Guizhou,China[J].Advances in Civil Engineering,2021,2021:1–16. doi: 10.1155/2021/8854194 [23] 高相波,李丽慧.热水河流域典型泥石流灾害成因机制与协同防治研究[J].工程地质学报,2020,28(5):1039–1048. doi: 10.13544/j.cnki.jeg.2020-288 Gao Xiangbo,Li Lihui.Initiation mechanism of typical debris flow hazard in hot water river basin and synergetic prevention study[J].Journal of Engineering Geology,2020,28(5):1039–1048 doi: 10.13544/j.cnki.jeg.2020-288 [24] 向伏林,杨天亮,顾凯,等.钻孔全断面分布式光纤监测中光缆–土体变形协调性的离散元数值模拟[J].岩土力学,2021,42(6):1743–1754. doi: 10.16285/j.rsm.2020.1420 Xiang Fulin,Yang Tianliang,Gu Kai,et al.Discrete element numerical simulation of cable-soil deformation compatibility in borehole distributed optical fiber monitoring[J].Rock and Soil Mechanics,2021,42(6):1743–1754 doi: 10.16285/j.rsm.2020.1420 [25] 宿辉,杨家琦,胡宝文,等.颗粒离散元岩石模型的颗粒尺寸效应研究[J].岩土力学,2018,39(12):4642–4650. doi: 10.16285/j.rsm.2017.0945 Su Hui,Yang Jiaqi,Hu Baowen,et al.Study of particle size effect of rock model based on particle discrete element method[J].Rock and Soil Mechanics,2018,39(12):4642–4650 doi: 10.16285/j.rsm.2017.0945 [26] Kong Xiaoang,Cai Guoqing,Cheng Yongfeng,et al.Numerical implementation of three-dimensional nonlinear strength model of soil and its application in slope stability analysis[J].Sustainability,2022,14(9):5127. doi: 10.3390/SU14095127 [27] 周喻,吴顺川,焦建津,等.基于BP神经网络的岩土体细观力学参数研究[J].岩土力学,2011,32(12):3821–3826. doi: 10.3969/j.issn.1000-7598.2011.12.045 Zhou Yu,Wu Shunchuan,Jiao Jianjin,et al.Research on mesomechanical parameters of rock and soil mass based on BP neural network[J].Rock and Soil Mechanics,2011,32(12):3821–3826 doi: 10.3969/j.issn.1000-7598.2011.12.045 [28] 向文,张强勇,张建国.坝区岩体蠕变参数解析–智能反演方法及其工程应用[J].岩土力学,2015,36(5):1505–1512. doi: 10.16285/j.rsm.2015.05.035 Xiang Wen,Zhang Qiangyong,Zhang Jianguo.Analytical intelligence inversion method of rheological parameters for dam zone rock mass and its application to engineering[J].Rock and Soil Mechanics,2015,36(5):1505–1512 doi: 10.16285/j.rsm.2015.05.035 [29] 季策,姜雨田.基于方向幅值比的欠定盲源分离算法[J].东北大学学报(自然科学版),2019,40(7):920–924. doi: 10.12068/j.issn.1005-3026.2019.07.002 Ji Ce,Jiang Yutian.Underdetermined blind source separation algorithm based on directional amplitude ratio[J].Journal of Northeastern University(Natural Science),2019,40(7):920–924 doi: 10.12068/j.issn.1005-3026.2019.07.002 [30] Liu Chun,Liu Hui,Zhang Hongyong.MatDEM-fast matrix computing of the discrete element method[J].Earthquake Research Advances,2021,1(3):100010. doi: 10.1016/j.eqrea.2021.100010
