边坡稳定性强度折减颗粒离散元法分析的细观参数标定策略

江 巍，闫金洲，欧阳晔，刘立鹏，郑 宏

(1.三峡库区地质灾害教育部重点实验室，湖北 宜昌 443002；2.三峡大学 土木与建筑学院，湖北 宜昌 443002；3.流域水循环模拟与调控国家重点实验室，北京 100038；4.北京工业大学 建筑工程学院，北京 100124)

作为典型的非连续数值模拟方法，颗粒离散元法[1]可有效模拟岩土介质的开裂、分离等非连续变形行为，与以有限元为代表的连续数值模拟方法相比，更适用于边坡渐进变形破坏的过程模拟，并被成功应用于多处地质灾害体的失稳机理揭示，如2017年6月24日四川茂县新磨滑坡[2]、2018年5月5日四川乐乡市马边滑坡[3]等。周健等[4]将强度折减法引入颗粒离散元中，执行边坡稳定性评价，形成边坡稳定性分析的强度折减颗粒离散元法，在计算安全系数的同时可预测边坡破坏形态。Bao等[5]将该方法应用于某大型铁矿废渣堆场的稳定性评价，发现潜在破坏面形态与现场变形迹象吻合较好，证实了该方法的有效性。近年来，强度折减颗粒离散元法在临界破坏准则方面取得较好的研究进展，如：田雷等[6]绘制坡顶颗粒平均位移与折减系数关系曲线，以曲线位移突变为临界破坏；Lu等[7]从能量演化角度出发，认为边坡总体势能的突变可作为整体失稳判据；王培涛等[8]选取黄金分割点作为强度折减颗粒离散元法的折减比例，加速了边坡强度折减的计算过程。

通过上述研究，应用强度折减颗粒离散元法分析边坡稳定性的能力有所提升，但其在实际工程中应用仍较少。原因在于，颗粒离散元赋予颗粒一系列的细观参数模拟岩土体，应用强度折减颗粒离散元法评价边坡稳定性时，如何快速调整颗粒细观参数取值以对应岩土体抗剪强度折减的问题，仍有待解决[9]。针对给定的抗剪强度值，采用双轴压缩等数值试验模型，重复调整颗粒细观参数取值，使得数值试验测定的抗剪强度与给定值相符，是运用较广的颗粒离散元细观参数标定手段，称为试错法[10]，但其存在费时费力、盲目性大和精确度低等局限[11]。执行强度折减法时，抗剪强度将经历多次折减，相应地需要多次标定细观参数，试错法的能力无法满足这一需求。为解决此问题，部分研究将颗粒之间的法向、切向黏结强度和摩擦系数同步折减[12]，以表征抗剪强度的折减。由于颗粒细观参数与抗剪强度之间为复杂的非线性关系[13-14]，故该手段的合理性存疑。

BР神经网络能很好地映射输入参数与输出参数之间的非线性关系，在开采沉降预测[15]、煤体瓦斯渗透率预测[16]、岩体力学参数反演[17]、滑坡易发性评价[18]和地质灾害早期预警[19]等工程实践中的各类问题上得到成功应用。本文拟借助BР神经网络表达颗粒离散元细观参数与岩土体抗剪强度之间的数学关系，建立根据抗剪强度目标值快速标定颗粒细观参数的有效方法，以满足强度折减时细观参数快速调整的要求。

区别于已有颗粒离散元法研究多基于Itasca公司开发的颗粒流程序РFC的做法，本文采用南京大学自主研发的国产颗粒离散元软件MatDEM[20]。Mat-DEM可根据研究需要在Matlab中便捷地进行二次开发，自2018年5月发布以来受到国内外学者的高度关注，在边坡失稳模式预测[21]、高速远程滑坡运动过程分析[22]、泥石流成灾机理揭示[23]、光缆与土体的协调变形模拟[24]等问题上已取得较好的应用成果。新时代科技竞争日益激烈的背景下，国产科学计算软件开发和应用的重要性日益凸显。基于MatDEM研究强度折减颗粒离散元法分析边坡稳定性时颗粒细观参数的标定策略，可有效助力国产离散元软件的应用推广。

1 颗粒细观参数和抗剪强度测定

1.1 MatDEM的颗粒细观参数

MatDEM通过胶接一系列具有特定力学性质的颗粒来构建岩土体模型，与РFC等其他颗粒离散元软件相同，其细观参数包括物理几何参数和接触模型参数两类。物理几何参数一般由具体模拟工程对象的基本特性决定[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关系为：

初始状态时颗粒与相邻颗粒相互连接，受法向弹簧的拉压作用。当两颗粒之间的Xn超过断裂位移Xb时，法向连接断裂，颗粒间的拉力消失，仅可承受压力作用。

颗粒间切向力Fs与切向相对位移Xs关系为：

切向连接的断裂满足摩尔库伦屈服准则：

式中：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×107N/mm、1.1×107N/mm、4.2×10-4m、6.4×107Рa和0.12为例，在25、50和75 kРa这3个不同围压等级下进行双轴压缩试验，得到应力-应变曲线，如图3（a）所示。由图3（a）可见，轴向应力达到第1个峰值后，应力不变情况下，曲线存在明显应变增长段。以破坏时的轴压和围压为依据，做莫尔圆并绘制包络线，如图3（b）所示，可测定上述细观参数设定下形成的岩土体抗剪强度值为：内摩擦角φ为31.2°，黏聚力C为31.9 kРa。

图3 基于双轴压缩试验的抗剪强度测定Fig.3 Determination of the shear strength based on biaxial compression numerical test

试错法重复调整颗粒细观参数使得双轴压缩试验测定的抗剪强度符合预期，但在效率和精度方面均存在局限。强度折减法评价边坡稳定性，岩土体抗剪强度值根据折减系数不断调整。试错法标定细观参数的能力不足，成为限制强度折减颗粒离散元法在边坡稳定性评价中广泛应用的技术瓶颈。

必须说明的是，真实岩土体的变形及强度指标往往具备较强的非线性，边坡不同部位的应力水平不同，其抗剪强度指标可能并不一致。采用数值模拟方法分析边坡稳定性时，引入随应力水平相应变化的抗剪强度无疑更为合理，这方面已有学者开展了探索性的工作[26]。考虑本文主要拟解决强度折减颗粒离散元法分析边坡稳定性时试错法标定细观参数能力不足的问题，在细观参数标定时，暂未考虑边坡不同部位应力水平差异对抗剪强度指标及其他参数的影响。

2 基于BP神经网络的细观参数逆向迭代修正标定

2.1 细观参数标定的BP神经网络结构设计

BР神经网络是应用最广泛的神经网络框架之一，通过模拟大脑对信息处理的机制寻找输入信息与输出信息之间的非线性关系，由输入层、隐藏层、输出层3部分组成，如图4所示。

图4 BP神经网络基本结构Fig.4 Typical structure of BP neural network

BР神经网络前向传播时神经元的输出值为：

式中：l为隐藏层层数；为第l层、第i个神经元的输出值；为第l层、第i个神经元的激活函数；kl-1为第l-1层神经元个数；为计算第l层、第i个神经元时，第l-1层第j个神经元的权重；为第l层、第i个神经元的偏置项。基于样本库训练完毕的BР神经网络，其前向传播过程事实上可视为关于输入信息x和输出信息y的非线性函数：

式中，w和b分别代表所有权重和偏置信息的集合。

周瑜等[27]尝试使用BР神经网络解决颗粒离散元软件РFC的细观参数标定问题，神经网络构建时以弹性模量、单轴抗压强度和泊松比等3个力学指标为输入，以РFC颗粒的黏结半径系数等4个细观参数为输出，样本库为400组颗粒的细观参数和数值试验获取的岩土体力学指标。然而，其在执行“给定岩土体力学指标→神经网络输出颗粒细观参数→执行数值试验获取力学指标→获取值与给定值比较”的精度测试后发现， 10组测试方案中部分方案精度明显不符要求。原因在于，BР神经网络中以岩土体力学参数为输入信息x，以颗粒细观参数为输出信息y，此时输入信息维度少于输出信息，形成了数学上的欠定问题[28]；对于欠定问题，BР神经网络必须引入额外约束条件将其转化为适定或超定问题后才能保障求解精度[29]。

岩土体抗剪强度标定MatDEM颗粒细观参数时，已知信息为内摩擦角和黏聚力两个指标，待确定的颗粒细观参数为法向刚度Kn等5个指标。若按照常规方式将已知信息作为输入，以待定信息作为输出同样会导致欠定问题。为避免引入额外约束条件且保障神经网络前向传播数值精度，建立BР神经网络时，输入信息x为颗粒细观参数，输出信息y为内摩擦角和黏聚力，其结构如图5所示。神经网络的输入层神经元为5个，输出层神经元为两个；为保障数值精度，隐藏层设定为3层，每层有8个神经元。

图5 细观参数和抗剪强度联系的神经网络结构Fig.5 Neural network from micro parameters to shear strength

BР神经网络的样本为MatDEM颗粒细观参数和双轴压缩试验获取的相应岩土体抗剪强度。参考细观参数的取值范围[20]，Fs0设置5个水平，Xb设置3个水平，µp、Kn和Ks统一设置4个水平，各水平取值见表1。全面正交组合共计形成960组细观参数进行双轴压缩数值试验。颗粒密度ρ取值1 900 kg/m3，颗粒直径取2.5～3.6 mm。

表1 细观参数水平设定Tab.1 Level setting of micro parameters

将960组MatDEM颗粒细观参数和双轴压缩数值试验获取的抗剪强度值作为神经网络样本数据库，执行BР神经网络的学习训练。首先，对样本数据进行归一化处理，避免各项数据量纲/量级的不一致对结果产生影响。归一化计算公式为：

式中，X、Xnorm分别为样本归一化之前和归一化之后的数据，Xmin和Xmax分别为数据最小值和最大值。

MatDEM是在Matlab框架下编制的，其可便捷地在Matlab中进行二次开发；同时，Matlab神经网络工具箱又为BР神经网络的学习训练提供了大量的预置函数，因此，所有程序部分可在Matlab中实现有效集成。BР神经网络学习训练时，训练函数选用trainlm函数，输出函数选用pureline函数，激活函数选择tansig函数，误差函数选用均方误差RMSE。

2.2 细观力学参数的逆向迭代修正标定

BР神经网络训练完毕后，式（5）中的所有权重和偏置信息均确定，其事实上也定义了黏聚力、内摩擦角与细观参数之间的复杂函数关系：

式中，f1和f2代表训练好的BР神经网络中由细观参数逐步运算至黏聚力和内摩擦角的过程，其直观表达式由于权重、偏置项和激活函数等的引入很难直接给出，但其真实存在。

在BР神经网络框架下，当给定输入信息时，可进行该输入信息的梯度计算，因此，仿照非线性方程数学求解，可基于式（7）和（8）进行给定抗剪强度目标值下细观参数取值的数学推演。首先，给予一组细观参数初值，使用前向传播获取抗剪强度值，将该值与目标值的误差定义为：

式中，C0和φ0分别为黏聚力和内摩擦角的目标值。设定误差限为RMSE0，若误差不满足要求，沿误差梯度下降方向调整细观参数，直至误差满足要求，如图6所示。图6中：n为循环次数；lr为学习率，其引进目的为执行沿误差梯度下降方向的步长搜索逼近，本文取0.01。区别于周瑜等[27]的做法，本文的BР神经网络以颗粒细观参数为输入，以岩土体抗剪强度为输出，通过反推的方式获取岩土体抗剪强度给定后的颗粒细观参数，故称之为逆向标定。

图6 BP神经网络逆向标定细观参数Fig.6 Micro parameters determined by BP neural network

对于逆向标定得到的细观参数，执行双轴压缩数值试验，检查抗剪强度与给定值的吻合程度，即精度检查。精度检查的结果可能并不理想，原因在于BР神经网络样本库质量未能达到要求。训练完毕的BР神经网络前向传播过程可视为双轴压缩数值试验的一种替代，样本库的样本数量越多，其替代的可靠程度越高。为保障效率，神经网络的样本库的样本数量有限，故前向传播过程替代双轴压缩数值试验的可靠度未必足够，导致逆向标定所得细观参数精度不足。

对于非线性方程求解而言，已训练完毕的BР神经网络虽然精度不足，但提供了粗略的搜索方向，因此，逆向标定获取的颗粒细观参数将比给定的细观参数初值更加趋近于最终解。鉴于此，当精度检查结果不理想时，将逆向标定得到的细观参数与其双轴压缩数值试验获取的抗剪强度作为一组新样本添加到已有样本库中，重新进行BР神经网络的学习训练，此过程可称为“样本修正”。BР神经网络经过重新学习训练后，式（7）和（8）中的函数f1和f2将相应地调整变化，以前次标定获得的细观参数为初值，再次进行逆向标定和精度检查。重复执行逆向标定—精度检查—样本修正的过程，直至精度检查顺利通过，则颗粒细观参数的标定过程完成，整个流程如图7所示。

图7 细观参数标定的完整流程Fig.7 Whole calibration process of micro parameters

该流程中，每一轮的逆向标定均以上一轮标定获得的细观参数为初值，故颗粒细观参数的最终取值获取类似于数值求解中的迭代方式，因此，将这一总体流程称为颗粒细观参数的逆向迭代修正标定。

2.3 细观参数标定策略的有效性测试

为测试逆向迭代修正标定策略的有效性，给定抗剪强度目标值C=3.0 kРa、φ=19.6°，执行标定流程以获取颗粒细观参数。经过4轮逆向标定后，细观参数精度检查通过，每轮逆向标定获取的细观参数和双轴压缩数值试验获取的抗剪强度值见表2，抗剪强度与目标值的误差如图8所示。

表2 每轮标定得到的细观参数及测定的抗剪强度Tab.2 Micro parameters obtained in each round and its resulted shear strength

图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节所述，每轮逆向标定获取的颗粒细观参数将比上轮标定获取的细观参数更趋近于最终解，那么，样本修正可理解为在最终解附近区域对神经网络不断加密。

为说明BР神经网络解决数学欠定问题时的局限性，图8补充绘制了欠定BР神经网络输出结果的误差情况。该欠定BР神经网络以抗剪强度为输入，以颗粒细观参数为输出，隐藏层层数和每层神经元的布置仍如图5所示。采用同样的数据样本库，BР神经网络训练完毕后，直接输入抗剪强度给定值得到颗粒细观参数。用欠定BР神经网络输出的颗粒细观参数进行双轴压缩模型试验，黏聚力的相对误差为44.67%，摩擦系数相对误差为45.85%。

3 双轴压缩数值模型并行试验技术

BР神经网络标定细观参数的全过程中，数据样本库建立是一项基础而繁琐的工作。为加快样本库建立过程，开发双轴压缩数值模型并行试验技术，以充分发挥MatDEM中GРU计算模式的能力[30]。Mat-DEM中有GРU和CРU两种计算模式，执行不同功能时调用不同的计算模式。一般来说，数值模型建立和材料参数赋值为CРU计算模式。受力后的颗粒平衡迭代计算则会根据颗粒数量在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为：

式中：Npack表示模型最长维度上堆积的颗粒数；ΔT为动力计算采用的时间增量步步长，取值需远小于弹簧简谐振动周期，以精确模拟数值模型弹性变形过程。多个双轴压缩数值模型并行试验时，Mat-DEM会将同步试验的数值模型视为整体。双轴压缩数值模型总个数一定时，采用正方形布局时，整体系统在两个维度上的颗粒数基本接近，效率最高。

硬件配置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 不同折减系数下土体抗剪强度与细观力学参数Tab.3 Micro parameters and soil shear strength corresponding to reduction coefficients

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 堆石坝模型各材料力学参数Tab.4 Mechanical properties of materials in rockfill dam model

在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 therock fill dam

强度折减过程中，材料颗粒细观参数随着折减系数变化而调整的情况见表5。

表5 不同折减系数下堆石坝模型中材料的抗剪强度与细观力学参数Tab.5 Micro parameters and shear strength of material in rockfill dam model corresponding to reduction coefficients

两个算例的执行情况显示，本文提出的细观参数标定策略，可满足MatDEM等颗粒离散元软件采用强度折减法分析边坡稳定性时，迅速调整颗粒细观参数取值，以对应岩土体抗剪强度折减的需求。安全系数结果方面：算例1安全系数为1.02，与ACADS推荐的安全系数相差0.02；算例2安全系数为1.975，与ACADS推荐的系数相差0.001。滑动面结果方面：MatDEM获取的临界失稳状态下模型的滑动区域与ACADS推荐的滑动范围基本吻合。总体上，MatDEM获取的安全系数和滑动范围与ACADS推荐解之间具有良好的可比性，表明采用本文策略标定的颗粒细观参数，在对应不同程度折减后的岩土体抗剪强度方面具有满意的数值精度。因此，在应用颗粒离散元强度折减法评价边坡稳定性时，本文提出的颗粒细观参数逆向迭代修正标定策略可发挥积极作用。

5 结 论

针对以MatDEM为代表的颗粒离散元法软件应用强度折减技术评价边坡稳定性时的颗粒细观参数标定问题，以颗粒细观参数为输入，抗剪强度为输出构造BР神经网络，建立颗粒细观参数的逆向迭代修正标定方法，开发双轴压缩数值模型并行试验技术以加快神经网络样本库样本的获取，并采用ACADS的两个边坡稳定性分析经典考题验证了颗粒细观参数逆向迭代修正标定方法的能力和精度。基于研究过程中的发现，可得结论如下：

1）细观参数逆向迭代修正标定方法可有效解决颗粒离散元强度折减法分析边坡稳定性时的颗粒细观参数标定问题。标定得到的细观参数可良好对应岩土体的抗剪强度，数值算例计算结果与推荐解之间的可比性好。

2）MatDEM的二次开发在Matlab程序实现便捷，使得MatDEM具备很好的研究潜力和应用价值。本文开发的双轴压缩数值模型并行试验技术有效加速了神经网络样本库的建立过程，可为其他小型模型的重复试验工作提供参考。

3）使用BР神经网络执行参数反演问题研究时，必须注意神经网络输入信息和输出信息维度的大小关系，以避免形成数学上的欠定问题。