基于BP神经网络的土体细观力学参数反演分析

Ｊｏｕｒｎａｌ　ｏｆＥｎｇｉｎｅｅｒｉｎｇ　Ｇｅｏｌｏｇｙ　工程地质学报　１００４－９６６５／２０１５／２３（４）．０６０９．０７　ＤＯＩ：１０．１３５４４／ｊ．ｃｎｋｉ．ｊｅｇ．２０１５．０４．００４　基于ＢＰ神经网络的土体细观力学参数反演分析冰　李澄清①　刘天为②　张海洋①　徐文杰①　（①清华大学水沙科学与水利水电工程国家重点实验室（②中国电建集团成都勘察设计研究院有限公司　成都摘要北京１０００８４）　６１００７２）　利用离散元方法对颗粒材料的细观力学特性研究，目前确定数值计算模型的细观力学参数大多数通过反复调试获　取，效率低、可重复性差。本文采用开源的颗粒离散元程序ＬＭＧＣ开展了土体双轴压缩数值试验，通过２５组土体细观力学参　数计算得到相应的宏观力学参数，建立了ＢＰ人工神经网络反演系统。利用土体物理试验得到的土体宏观力学参数，输入ＢＰ　神经网络，反演得到土体的细观力学参数。将所得细观力学特性参数输入所建立的土体数值计算模型，得到土体破坏过程中　的应力一应变关系曲线，以及土体颗粒的力链图和旋转变形云图。所建立的土体数值试验模型能够较好地模拟土体变形破坏　过程，利用ＢＰ神经网络反演细观力学参数以及数值模型计算得到的土体宏观力学参数与物理试验吻合较好，误差在１０％左　右，土颗粒间力链云图以及旋转变形云图较好地揭示了土体变形破坏的机理。　关键词　ＢＰ神经网络中图分类号：ＴＵ４３　反演离散元细观力学参数颗粒力学　文献标识码：Ａ　ＢＡＣＫ．ＡＮＡＬＹＳＩＳ　０Ｎ　ＭＩＣＲｏＭＥＣＨＡＮＩＣＡＬ　ＰＡＲＡＭＥＴＥＲＳ　ｏＦ　ＳｏＩＬ　ＭＡＳＳ　ＵＳＩＮＧ　ＢＰ　ＮＥＵＲＡＬ　ＮＥＴＷｏＲＫ　ＬＩ　Ｃｈｅｎｇｑｉｎｇ①ＬＩＵ　Ｔｉａｎｗｅｉ②ＺＨＡＮＧ　Ｈａｉｙａｎｇ①ＸＵ　Ｗｅｎｊｉｅ①　（（￣Ｓｔａｔｅ　Ｋｅｙ　Ｌａｂｏｒａｔｏｒｙ　ｏｆＨｙｄｒｏｓｃｉｅｎｃｅ　ａｎｄ　Ｈｙｄｒａｕｌｉｃ　Ｅｎｇｉｎｅｅｒｉｎｇ，Ｔｓｉｎｇｈｕａ　Ｕｎｉｖｅｒｓｉｔｙ，Ｂｅｌｉｔｎｇ　１０００８４）　（Ｐｏｗｅｒ　Ｃｈｉｎａ　Ｃｈｅｎｇｄｕ　Ｅｎｇｉｎｅｅｒｉｎｇ　Ｃｏｒｐｏｒａｔｉｏｎ　Ｌｉｍｉｔｅｄ，Ｃｈｅｎｇｄｕ　６１００７２）　Ａｂｓｔｒａｃｔ　Ｔｈｅ　ｍｉｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　ｏｆ　ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌ　ａｒｅ　ｏｂｔａｉｎｅｄ　ｂｙ　ｒｅｐｅａｔｅｄ　ｔｅｓｔｓ．Ｒｅｓｅａｒｃｈｅｓ　ｏｎ　ｐａｒｔｉｃｌｅ　ｍａｔｅｒｉａｌ　ｂａｓｅｄ　ｏｎ　ＤＥＭ　ａｒｅ　ｄｏｎｅ．Ｔｈｉｓ　ａｒｔｉｃｌｅ　ｂｕｉｌｄｓ　ａ　ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌ　ｏｆ　ｓｏｉｌ　ｕｓｉｎｇ　ＬＭＧＣ，ａｎ　ｏｐｅｎ　ｓｏｕｒｃｅ　ｓｏｆｔｗａｒｅ．Ａｎｄ　ｔｈｅ　２５ｇｒｏｕｐｓ　ｏｆ　ｍａｃｒｏｍｅｃｈａｎｉｃａｌ　ａｎｄ　ｍｉｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　ｃａｌｃｕｌａｔｅｄ　ｂｙ　ｔｈｅ　ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌ　ｂｕｉｌｄ　ａ　ＢＰ　ｎｅｕｒａｌ　ｎｅｔｗｏｒｋ．Ｔｈｅ　ＢＰ　ｎｅｕｒａｌ　ｎｅｔｗｏｒｋ　ｃａｎ　ｂａｃｋ—ｃａｌｃｕｌａｔｅｓ　ｔｈｅ　ｍｉｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　ｉｆ　ｔｈｅ　ｍａｃｒｏｍｅｃｈａｎｉｃａ１　ｐａｒａｍｅｔｅｒｓ　ｏｆ　ｒｅａｌ　ｐｈｙｓｉｃａｌ　ｔｅｓｔｓ　ａｒｅ　ｉｎｐｕｔｔｅｄ　ｉｎｔｏ　ｔｈｅ　ｎｅｕｔｒａｌ　ｎｅｔｗｏｒｋ．Ｔｈｅ　ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌ　ｃａｎ　ｇｅｔ　ｔｈｅ　ｓｔｒｅｓｓ—ｓｔｒａｉｎ　ｌｉｎｅ．ｓｔｒｅｓｓ　ｔｒａｉｎ　ａｎｄ　ｐａｒｔｉｃｌｅ　ｒｏｔａｔｉｏｎ　ｇｒａｐｈ　ｉｆ　ｔｈｅ　ｍｉｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　ａｒｅ　ｓｅｔ．Ｔｈｅ　ｒｅｓｕｌｔｓ　ｓｈｏｗ　ｔｈａｔ　ｔｈｅ　ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌ　ｃａｎ　ｗｅｌｌ　ｓｉｍｕｌａｔｅ　ｔｈｅ　ｓｏｉｌ　ｄｅｓｔｒｕｃｔｉｏｎ　ｐｒｏｃｅｓｓ．Ｔｈｅ　ｓｔｒｅｓｓ—ｓｔｒａｉｎ　ｌｉｎｅ　ｏｆ　ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌ　ｍａｔｃｈｅｓ　ｗｅｌｌ　ｗｉｔｈ　ｔｈｅ　ｒｅｓｕｌｔ　ｏｆ　ｐｈｙｓｉｃａｌ　ｔｅｓｔ．Ａｎｄ　ｅｒｒｏｒ　ｏｆ　ｍａｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　ｃａｌｃｕｌａｔｅｄ　ｂｙ　ＢＰ　ｎｅｕｒａｌ　ｎｅｔｗｏｒｋ　ｉｓ　ａｂｏｕｔ　１０％．Ｂｅｓｉｄｅｓ，ｔｈｅ　ｓｔｒｅｓｓ　ｔｒａｉｎ　ａｎｄ　ｐａｒｔｉｃｌｅ　ｒｏｔａｔｉｏｎ　ｇｒａｐｈ　ｒｅｖｅａｌ　ｔｈｅ　ｍｅｃｈａｎｉｓｍ　ｏｆ　ｓｏｉｌ　ｄｅｓｔｒｕｃｔｉｏｎ　ｐｒｏｃｅｓｓ．　Ｋｅｙ　ｗｏｒｄｓ　ＢＰ　ｎｅｕｔｒａｌ　ｎｅｔｗｏｒｋ，Ｂａｃｋ　ｃａｌｃｕｌａｔｉｏｎ，ＤＥＭ，Ｍｉｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ，Ｐａｒｔｉｃｌｅ　ｍｅｃｈａｎｉｃｓ　收稿日期：２０１５—０３—２２；收到修改稿日期：２０１５—０５—１８．　基金项目：国家自然科学基金青年基金项目（５１１０９１１７），清华大学自主科研计划项目（２０１１１０８１１２５）资助．　第一作者简介：李澄清（１９９２一），男，硕士生，主要从事岩土颗粒力学研究工作．Ｅｍａｉｌ：ｌｉｃｈｅｎｇｑｉｎｇ０２＠１６３．ｃｏｍ　６ｌ０　０引　言　Ｃｕｎｄａｌｌ于１９７１年提出了适用于岩石的离散单　元法（Ｃｕｎｄａｌｌ，１９７１），并于１９７９年提出了适用于土　体的离散单元法（Ｃｕｎｄａｌｌ　ｅｔ　ａ１．，１９７９）。２０世纪８０　年代，王泳嘉引入Ｃｕｎｄａｌｌ的离散元法研究岩石力　学和颗粒系统（王泳嘉，１９８６），开启了国内离散元　研究的先河。离散元法在非连续介质方面有着广阔　的应用前景，较适合岩土体材料力学变形特性的研　究，在隧道工程、边坡工程、地基处理、工程地质、采　矿工程等方面，得到十分广泛的应用。　连续介质方法中采用的力学变形特征参数，如　抗压强度、弹性模量、黏聚力、内摩擦角等，均可通过　力学试验直接获取。而对于非连续介质的离散元　法，其细观力学参数与宏观力学参数存在很大的差　别，两者间的关系较为复杂，细观力学参数无法通过　宏观力学试验或数学表达式直接获取（周喻等，　２０１１）。在离散元法中，建立数值模型需要提前定　义单元间的接触模型，单元问接触模型的细观参数　对数值试验结果影响很大（徐泳等，２００３）。例如：　ＰＦＣ、ＹＡＤＥ（Ｋｏｚｉｃｋｉ　ｅｔ　ａ１．，２００８，２００９）、ＬＭＧＣ　（Ｄｕｂｏｉｓ　ｅｔ　ａ１．，２００３）等离散元程序，需要根据选定　的单元接触模型，假定细观力学参数进行数值试验，　并将计算得到的宏观力学参数与物理试验得到的宏　观力学参数进行对比，通过反复调试，当计算结果与　物理试验结果相近时，选取该组细观力学参数，用于　计算其他数值模型。目前，离散单元间接触模型的　细观力学参数大多是根据物理力学试验反复调试获　取，这一过程耗时长、效率低且带有很大程度的盲目　性（贾学明等，２０１０）。因此，需要引入一种较为科　学的方法来获取离散元数值模型中的细观力学参　数，建立岩土体等颗粒材料宏一细观力学参数之问　的联系。　本文计算采用的ＬＭＧＣ９０由法国Ｕｎｉｖｅｒｓｉｔ６　Ｍｏｎｔｐｅｌｌｉｅｒ开发，为一款致力于数值模拟多种形态、　多种材料、大规模接触问题的开源数值计算平台。　作为开源的离散元计算软件，ＬＭＧＣ９０中的物理模　型（物态法则、接触定理等），数值单元模型（有限　元、离散元），技术特点（接触检测、可视化、并行性　等）和数值策略（时间积分、数值解算器等）在充分　了解掌握程序的基本理论与运行法则的前提下，均　可按使用者的意图进行二次开发（Ｄｕｂｏｉｓ　ｅｔ　ａ１．，　２００３）。　Ｊｏｕｒｎａｌ　ｏｆ　Ｅｎｇｉｎｅｅｒｉｎｇ　Ｇｅｏｌｏｇｙ工程地质学报２０１５　人工神经网络是处理非线性关系的有力工具，　能够建立输入样本和输出样本之间的非线性关系。　岩土工程中存在很多已知量和未知量之间较为复杂　非线性关系，这种非线性关系可以通过神经网络得　到较好的映射（高浪等，２００２；毛健等，２０１１）。本　文基于ＬＭＧＣ离散元软件，采用ＢＰ神经网络方法，　建立了土体宏细观力学参数的联系，输入物理力学　试验得到的宏观力学变形特征参数，即可反演计算　离散元数值模型中的细观力学参数，并针对该种方　法的可靠性进行了分析与探究。最后，根据反演得　到土颗粒细观力学参数，利用离散元方法计算得到　土体应力一应变关系曲线，颗粒问力链图以及颗粒　旋转变形图，从细观角度分析了颗粒的变形破坏特　征　１试验原理和设备　作者从唐家山堰塞体选取了部分典型的土石混　合体，制作直径为１００ｍｍ，高度为２００ｍｍ的圆柱形　试样，对其进行土体三轴试验。试验应力路径为：　（１）采用静水压力条件施加应力ｏｒ　＝ｏｒ　：　１００ｋＰａ（　为轴向应力；　为侧向应力）至预定压　力水平使试样排水固结；（２）排水固结稳定后，试验　采用轴向位移控制方法对试样进行压缩剪切。得到　土体应力一应变关系曲线（图１）。　图１　土体三轴试验的应力一应变关系　Ｆｉｇ．１　Ｓｔｒｅｓｓ—ｓｔｒａｉｎ　ｒｅｌａｔｉｏｎｓｈｉｐ　ｏｆ　ｓｏｉｌ　ｔｒｉａｘｉａｌ　ｔｅｓｔ　为了与室内试验进行对比并且简化计算，建立　了二维平面下的土颗粒双轴压缩离散元数值模型，　宽为１００ｍｍ，高为２００ｍｍ，颗粒直径为２～５ｒａｍ，颗　粒总数量为２７３２个。为避免因刚性约束造成左右　边界上的应力集中，模型的左右边界采用了柔性膜　边界，上下边界则为刚性的加载板（图２）。　６１２　Ｊｏｕｒｎａｌ　ｏｆ　Ｅｎｇｉｎｅｅｒｉｎｇ　Ｇｅｏｌｏｇｙ工程地质学报】Ｉ＼　氇粤　２０１５　图４给出了在土体双轴压缩数值试验中，距离　阈值Ｗ固定为０．０１ｒａｍ、细观黏聚力ｃ固定为１０ｋＰａ、　细观内摩擦角　分别为１０。、１５。、２０。、３０。下的应　力一应变关系。从图中分析可知，细观内摩擦角对　输入层、隐藏层、输出层这３部分，层与层之间通过　权重系数连接。输入层和输出层的神经元数量由具　体要研究问题的参数数量决定，而隐藏层单元数量　由问题的复杂程度、误差收敛程度决定（易宏杰，　数值模型的宏观力学变形特性有明显影响。不同细　观内摩擦角下，土体双轴压缩数值试验，其应力一应　变关系特征基本相似，随着细观内摩擦角的增大，数　值试样的峰值应力有一定增长、弹性模量有较大提　高。细观内摩擦角在细观层面决定了颗粒间发生相　互运动的难易程度与所需作用力大小，于宏观层面　在较大程度上决定了试样的弹性模量与承载能力。　图４细观内摩擦角对应力一应变关系的影响　Ｆｉｇ．４　Ｅｆｆｅｃｔ　ｏｎ　ｓｔｒｅｓｓ—ｓｔｒａｉｎ　ｒｅｌａｔｉｏｎｓｈｉｐ　ｏｆ　ｉｎｔｅｒｎａｌ　ｆｒｉｃｔｉｏｎ　ａｎｇｌｅ　３土体细观力学参数反演分析　３．１　ＢＰ神经网络基本原理　本文所采用的ＢＰ神经网络是一种按误差逆向　传播算法训练得到的多层前馈网络，在岩土工程领　域应用较为广泛。ＢＰ神经网络的优势在于能存储　与学习大量输入一输出映射关系，无需提前揭示该　种映射关系的数学表达式。ＢＰ神经网络的学习法　则为均方误差值最快速下降法，能通过反向传播不　断地调整神经网络的权重与阈值从而使神经网络的　总体均方误差值最小。１９８９年Ｒｏｂｅｒｔ　Ｈｅｃｈｔ—Ｎｉｅｌｓｏｎ　证明了任意复杂的非线性关系，都可以用仅含有一　个隐含层的ＢＰ神经网络来逼近，隐含层过多，会导　致反演计算耗时太长，甚至得不到反演结果，而计算　精度并不会有显著提升。本文根据实际需求，所建　立的神经网络模型中只包含一个隐含层。　图５是本文采用的神经网络模型原理图，包括　２００４）。建立土体宏观力学参数与离散元数值模型　细观力学参数联系，主要分５个工作流程：（１）输入　用于训练人工神经网络的宏观力学参数与细观力学　参数；（２）训练人工神经网络，细调神经网络的相应　参数，使训练样本中的宏一细观力学参数能建立较　好的联系；（３）人工神经网络建立后，输入测试宏观　力学参数，利用人工神经网络反演细观力学参数；　（４）利用反演得到的细观力学参数建立数值计算检　验模型进行数值试验，将数值试验计算得到的宏观　力学参数与输入的测试宏观力学参数进行对比，检　验人工神经网络的合理性；（５）上述工作流程不断　重复，直至效果较好的人工神经网络得以建立。　输入层　隐藏层　输出层　宏　观　力　学　参　数　图５　ＢＰ神经网络模型　３．２　ＢＰ神经网络的建立、训练　利用ＢＰ神经网络分析方法，针对围压为　１００ｋＰａ下的土体双轴压缩数值试验模型，建立人工　神经网络细观参数反演系统。选取细观力学参数特　征值：细观黏聚力Ｃ、细观内摩擦角　…；选取宏观力　学参数特征值：弹性模量Ｅ、弹性阶段的比例极限　、峰后强度　（图６）。其中宏观力学参数作为输　入层参数，细观力学参数作为输出层参数。　选取２５组不同细观参数组合的数值模型作为　土体双轴压缩数值试验的宏细观参数反演系统的训　练样本，训练样本中细观黏聚力Ｃ分别为１０ｋＰａ、　３０ｋＰａ、５０ｋＰａ、７０ｋＰａ、ｌＯＯｋＰａ；细观内摩擦角　分　别为１０。、１５。、２０。、２５。、３０。；细观距离阈值Ｗ因其　２３（４）　李澄清等：基于ＢＰ神经网络的土体细观力学参数反演分析　６１３　对力学性能影响不大，固定为０．０１ｒａｍ。根据所选　取的２５组不同细观力学参数组合，利用建立的离散　元数值计算模型，得到试验的宏观力学参数。宏一　细观力学参数取值范围如表１所示。　表１　细观力学和宏观力学参数取值范围　Ｔａｂｌｅ　１　Ｖａｌｕｅ　ｒａｎｇｅ　ｏｆ　ｍａｃｒｏｍｅｃｈａｎｉｃａｌ　ａｎｄ　ｍｉｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　为了消除量纲不同，便于统计分布规律，将用于　训练神经网络的２５组不同参数组合进行归一化处　理再输入ＢＰ神经网路，处理格式如式（１）所示。　ｐ　＝　—　ｐ　——，　Ｐ　∈Ｌ，ｊ∈［０，１］　（１）　Ｐｍ“一Ｐｍｉｎ　式中，Ｐ　为归一化数据处理值；Ｐ为样本初始值；　Ｐ…为样本初始最大值；Ｐ　ｉ　为样本初始最小值。　图６土体双轴压缩试验宏观力学参数特征值　Ｆｉｇ．６　Ｍａｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　ｃｈａｒａｃｔｅｒｉｓｔｉｃ　ｖａｌｕｅ　ｏｆ　ｓｏｉｌ　ｂｉａｘｉａｌ　ｃｏｍｐｒｅｓｓｉｖｅ　由于输入的细观力学参数和宏观力学参数都为　３个，所以选取的输入层和输出层神经元数目均为　３，经调试神经网络内部设置参数，所选取的人工神　经网络选取隐藏层神经元数量为３，已经可以达到　比较好的精度要求和速度要求。将所得到的２５组　宏一细观力学特性参数作为训练样本，输入所建立　的ＢＰ神经网络，便可以建立宏一细观力学特性参数　之间的非线性关系。　３．３　ＢＰ神经网络可靠性验证　为检测人工神经网络的可靠性，根据神经网络　训练模型中宏观力学特征参数的取值范围，随机生　成了４组宏观力学特征参数的测试组合样本，用以　测试神经网络的反演能力，反演计算得到的相应细　观力学参数经过式（２）的反归一化处理。随机生成　的土体宏观力学特性参数以及神经网络反演计算得　到的细观力学特性参数（表２）。　表２　ＢＰ神经网络测试样本的宏观力学参数　Ｔａｂｌｅ　２　Ｍａｃｒｏｍｅｃｈａｎｉｃａｌ　ｐａｒａｍｅｔｅｒｓ　ｏｆ　ＢＰ　ｎｅｕｔｒａｌ　ｎｅｔｗｏｒｋ　ｔｅｓｔ　ｓａｍｐｌｅ　ｑ＝ｑ　（ｑ　一ｑ　ｉ　）＋ｑ　ｉ　（２）　根据反演得到的土体细观力学特性参数，利用所建　立的土体离散元数值计算模型，计算得到宏观力学　特性参数，并与随机生成的宏观力学特性参数进行　对比（图７）。　计算结果表明，土体双轴数值试验ＢＰ神经网　络的参数反演的总体性能较为良好，具有较高的精　度。数值试验的反演计算参数与测试原始参数在宏　观力学特征参数弹性模量Ｅ、比例极限ｏｒ　、峰后强度　ｒｏ　上吻合较好。　３．４　ＢＰ神经网络参数反演与物理试验对比　已经建立土体双轴压缩试验的离散元数值计算　模型，以及建立土体宏观参数与细观参数的ＢＰ神　经网络反演模型，并且验证了其可靠性。为了进一　步验证ＢＰ神经网络参数反演在土体剪切试验中的　可靠性，对１００ｋＰａ围压下常规三轴压缩的真实物理　力学试验进行细观力学参数反演。根据真实物理试　验得到的宏观力学参数，利用所建立的ＢＰ神经网　络，反演得到的土体细观力学参数为：黏聚力ｃ为　２３．０３２ｋＰａ，内摩擦角　～为１３．１。，距离阈值　为　０．０１ｍｍ。根据反演得到的细观力学参数，利用所建　立的离散元数值计算模型进行计算，得到土体压缩　剪切试验的应力一应变曲线（图８）。表３是根据真　实物理实验得到的宏观力学参数，以及利用反演得　到的细观力学参数进行数值计算得到的宏观力学参　数。　图８和表３表明反演试验与物理试验在力学特　征参数上误差较小，所关心的土体宏观参数弹性模