CN112163328B - 一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法 - Google Patents

Abstract

本发明提供了基于深度学习和数据驱动建立岩土颗粒材料本构模型的方法，包括：采用离散元数值试验代替室内试验进行岩土颗粒材料的宏细观力学特性研究；对具有不同结构特征的颗粒材料数值试样进行不同加载路径下的离散元数值试验，得到大量的应力应变关系数据；采用累积绝对应变作为岩土颗粒材料的状态变量，用于描述颗粒材料的当前状态；采用改进的长短期记忆神经网络单元，搭建深度学习网络，将数值试样的结构特征和初始状态等非时序数据输入到网络的首个单元，将加载过程中的应变增量依次输入到网络中，输出颗粒材料当前的应力和其他状态量；采用基于奇异值分解的降噪方法对训练数据进行降噪处理，防止网络训练不收敛或误差过大。

Description

一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法

技术领域

本发明属于岩土颗粒材料本构模型研究领域，涉及到岩土颗粒材料的本构理论、离散元数值试验和深度学习，采用改进的长短时记忆神经网络挖掘基于离散元数值试验得到的颗粒材料应力应变关系数据，建立基于数据驱动的岩土颗粒材料本构关系，为岩土颗粒材料宏观和细观的多尺度耦合提供了新思路。

背景技术

岩土工程、水利工程、道路桥梁工程等的存在环境或对象相当大一部分是由离散形式的颗粒材料组成，如黏土、砂土、碎石、粗粒土等，可统称为岩土颗粒材料。岩土颗粒材料的力学特性直接影响地(路)基、堤防、堆石坝等工程的变形和稳定。岩土颗粒材料看似简单，但由于其离散性、摩擦性和多相性，导致其物理力学特性复杂，存在非线性、压硬性、剪胀性、各向异性、应力路径相关性等特性。这些复杂的力学特性与其离散性、多尺度结构和能量耗散机制有关。

为了描述岩土颗粒材料的复杂力学特性，学者们在连续介质力学框架下建立了不同的本构模型。建立本构模型与相关的试验研究，成为现代土力学研究的核心话题之一。关于本构模型的研究最早可以追溯到1773年库伦提出库伦屈服准则，直到1963年Roscoe等人建立了第一个土的弹塑性模型—剑桥模型，标志着人们在土体力学特性认识上的一次飞跃。自此以后，学者提出了数百个土的本构模型，如修正剑桥模型、邓肯EB模型、清华模型、南水双屈服面模型、多屈服面模型等，用以反映岩土颗粒材料的剪胀性、压硬性、非线性、各向异性、流变和率相关等特性。

近年来，岩土颗粒材料的本构模型研究还在不断地发展和创新，取得了丰富的研究成果。由于颗粒破碎会引起级配连续变化并导致压缩性增大和强度降低，因此其对岩土颗粒材料的力学特性影响较大。学者们分别基于损伤力学、颗粒破碎耗散能、广义弹塑性力学框架、亚塑性本构理论等建立了考虑颗粒破碎的本构模型。除了颗粒破碎，流变变形、湿化变形、缩尺效应等也是岩土颗粒材料力学特性的重要方面，需要反映到本构模型中。陈生水等提出了一个可以统一考虑由加载产生的静力变形、随时间的流变和应力松弛等特性的粗粒土弹塑性本构模型。姚仰平等提出引入统一硬化参数的系列UH模型，将其扩展到考虑时间效应的粒状材料，并成功应用到实际工程的应力变形模拟。Zhou等在热力学框架下，提出了一个基于率相关破碎力学和破碎耗散能尺寸效应的本构模型，能反映缩尺效应、颗粒破碎、流变等特性，在水布垭堆石坝的模拟结果与实测变形吻合的较好。本构模型的发展使其能更全面地考虑岩土颗粒材料复杂的力学特性，但也引入了更多的模型参数，有些参数能够通过试验率定，而有些参数的确定较为困难。

除了模型参数取值问题外，基于连续介质力学的宏观唯象本构模型还存在以下几个问题：宏观唯象的建模思路只考虑了宏观的应力应变关系，不能描述颗粒层面的力学行为；常规试验装置和试验技术不能探究岩土颗粒材料在复杂加载路径下的力学特性；宏观唯象本构模型多是基于传统本构理论建立，这些理论大多根据金属材料试验提出的，其中包含一些物理假设，诸如Drucker公设、关联或非关联流动准则等。因此，基于弹塑性理论的岩土颗粒材料本构模型在实际应用中还存在众多亟待解决的问题。

发明内容

相较于宏观唯象本构模型基于室内物理实验，本发明基于离散元数值试验得到具有不同结构特征的颗粒材料在不同加载路径下的应力应变数据，形成大规模的应力应变关系数据集，采用深度学习网络挖掘应力应变数据中潜藏的本构关系，得到基于深度学习和数据驱动的岩土颗粒材料本构模型。不同于宏观唯象本构模型的建模思路，该方法不需要引用过多的假设，纯粹通过数据驱动的方式，建立具有较好泛化能力和较强预测能力的岩土颗粒材料本构模型。

本发明提供了一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其包括以下步骤：

步骤1，首先进行离散元数值模拟的敏感性分析，确定离散元数值试验中加载速度、时间步、颗粒集合体中颗粒数目；

步骤2，制备具有不同结构特性的颗粒材料数值试样，进行不同加载路径下的离散元数值试验，获得不同结构特征的数值试样在不同加载路径下的应力应变特性数据；

步骤3，整理不同数值试样在不同加载路径得到的应力应变数据，提取所需数据，包括三个方向应变、应力、颗粒材料孔隙比、颗粒体系结构描述参数、颗粒材料杨氏模量和摩擦系数，以及累积绝对应变，对提取的所需数据进行降噪处理，建立原始训练数据库；

步骤4，按间隔从每一试样原始数据中抽取指定帧数，组合成时序数据序列构造深度学习训练集，非时序输入、时序输入和时序输出组合起来表示一个训练样本，对训练集进行对数化和归一化处理；

步骤5，搭建深度学习网络模型，深度学习网络模型由改进的LSTM网络单元和全连接层组合而成；

深度网络模型的输入包含时序输入和非时序输入，其中时序输入包含岩土颗粒材料在加载过程中x，y，x三个方向的轴向应变和累积绝对应变；输入时序数据前，将初始细胞状态C₀和初始隐含状态h₀定义为当前岩土颗粒材料的孔隙比、应力，并在输入初始化信息前加入全连接层，对初始化信息进行预处理，初始化信息包括时序数据的第一个时间步的孔隙比、应力和应变，用于初始化细胞状态和隐含状态；

非时序输入包括杨氏模量、摩擦系数和颗粒体系结构描述参数β，所述非时序输入不随加载进行发生变化，且对加载过程持续产生影响，因此，非时序输入每个时间步均输入给LSTM单元；

时序输出为材料的各向应力

和颗粒材料孔隙比时序数据，上标i表示荷载步；

步骤6，采用早停策略防止深度学习网络过拟合，将训练数据分批次输入深度网络模型进行训练，网络收敛后保存网络结构；

步骤7，调用训练好的深度学习网络模型，输入颗粒材料的结构信息和时序应变数据，预测岩土颗粒材料的时序应力数据和结构参数。

进一步的，步骤1中确定模拟加载速度为0.1m/s，时间步大小为5e^-7s，颗粒数目为10000。

进一步的，步骤2中制备具有不同结构特性的颗粒材料数值试样，即制备具有不同级配的颗粒材料数值试样，级配曲线由具有最小粒径截断的分形公式(3)确定，

上式中，F(d)表示累计质量分数，d表示颗粒粒径，d_max和d_min分别表示最大和最小颗粒粒径，β为颗粒体系结构描述参数。

进一步的，步骤2中进行不同加载路径下的离散元数值试验的具体实现方式如下；

加载路径设置：等p等b加载、等σ₃等b、排水常规三轴压缩、排水常规三轴拉伸、不排水常规三轴压缩、不排水常规三轴拉伸，其中

为平均静水压力，σ₃为小主应力，b＝(σ₂-σ₃)/(σ₁-σ₃)为中主应力系数；

加卸载路径设置：设置了在常规三轴状态下的加卸载循环路径；

围压、孔隙比设置：设定10级加载围压，分别为0.01、0.05、0.1、0.2、0.5、1.0、2.0、4.0、8.0和16.0MPa；通过在制样阶段设置不同颗粒间摩擦来制备具有不同初始孔隙比的数值试样，相同围压下摩擦系数越小得到的数值试样孔隙比越小，设定制样摩擦系数分别为0.001、0.05、0.1、0.2、0.3。

进一步的，步骤3中累积绝对应变，用于表示岩土颗粒材料当前状态的唯一性，计算方式如下，

应变增量由塑性应变和弹性应变组成，如下公式(4)所示，

dε＝dε^e+dε^p (4)

其中，dε表示应变增量，dε^e表示弹性应变增量，dε^p表示塑性应变增量；

累积绝对应变的计算公式如式(5)所示，

上式中，ε^m表示第m荷载步时，累积应变增量。

进一步的，步骤3中基于奇异值分解对应力应变数据进行降噪，具体实现方法如下，

一个m×n阶矩阵，其中所有元素属于实全域，则存在一个分解使得下式成立，

M＝UΛV(6)

其中，Λ为对角矩阵，U和V分别为m×m和n×n阶的复数方块矩阵，其具有以下性质

U^TU＝UU^T＝I (7)

U^-1＝U^T (8)

上式中，I表示m阶单位阵；

奇异值分解得到对角阵Λ，将对角阵从大到小排序，通过选取前k个对角元素，其他元素置为0，得到简化后的对角阵∑，采用以下公式恢复数据，

N＝UΣV(9)

矩阵N即为降噪后的数据。

进一步的，步骤5中改进的LSTM网络单元由初始状态、输入门、遗忘门和输出门组成，数据从进入单元到输出经历遗忘阶段、更新阶段和输出阶段，具体过程如下，

1.初始状态

时序训练数据输入单元前，首先定义t₀时刻的LSTM单元的细胞状态和隐含层结构，由于初始化信息：孔隙比e、应力σ_x,σ_y,σ_z、应变ε_x,ε_y,ε_z维度为7，无法与LSTM单元隐含层的hidden_unit_num一致，因此在输入初始化信息前加入全连接层对信息进行预处理，如下公式所示，

C₀＝W_c·C₀+b_c (10)

h₀＝W_h·h₀+b_h (11)

上述公式中，W_c和W_h分别为细胞状态和隐含状态对应的权值矩阵，b_c和b_h分别为细胞状态和隐含状态对应的偏置矩阵；

2.遗忘阶段

这一阶段“遗忘”多余信息，保留重要信息，LSTM单元首先通过遗忘门,对上一步输入h_t-1和该步输入x_t进行筛选，然后输出一个0～1间的向量f_t，如公式(12)所示，

f_t＝σ(W_f·[h_t-1,x_t]+b_f) (12)

上式中，σ表示Sigmoid函数，W_f表示权值矩阵，b_f表示偏置值；

3.更新阶段

更新阶段综合历史信息和当前步输入信息更新当前步的细胞状态，让上一个节点输入的h_t-1和x_t通过输入门以此决定更新后的信息i_t；然后再让h_t-1和x_t信息通过一个tanh层得到新的候选细胞信息

旧的细胞信息C_t-1通过遗忘门选择忘记旧细胞信息的一部分，通过输入门选择添加候选细胞信息，变为新的细胞信息C_t，如公式(13)所示：

4.输出阶段

输出阶段综合更新后的细胞状态信息和该步输入信息计算输出，更新完细胞状态后需要根据输入的h_t-1和x_t来判断输出细胞的状态特征，因此需要将信息通过输出门得到判断条件O_t，然后将更新后的细胞状态C_t经过tanh层得到一个-1～1之间值的向量，该向量与输出门得到的判断条件相乘就得到了最终该LSTM单元的输出h_t，如公式(14)所示：

与现有技术相比，本发明的优点和有益效果：本发明基于离散元数值试验得到具有不同结构特征的颗粒材料在不同加载路径下的应力应变响应，采用深度学习和数据驱动的方法建立岩土颗粒材料的本构关系。

受物理试验装置和试验技术的限制，难以进行复杂加载路径、高围压下的颗粒材料试验，本方法采用的离散元数值试验能够进行多种加载路径下的三轴试验，试验条件更加灵活，便于调整，获取应力应变数据的成本较低，为深度学习方法应用于岩土颗粒材料本构建模提供了高质量数据基础。

相较于宏观唯象本构模型，本方法从数据本身入手，利用深度学习和数据驱动方法进行数据挖掘，不需做过多物理假定，能够更准确地描述岩土颗粒材料在复杂加载路径下的应力应变特性。

附图说明

图1为本发明的技术路线图。

图2为不同颗粒级配试样的级配曲线。

图3为不同降噪阶数奇异值分解降噪效果。

图4为深度学习网络结构。

图5为过拟合现象。

图6为基于深度学习的岩土颗粒材料本构模型预测效果。

图7为时间序列数据构造方法。

图8为时间序列训练数据。

图9为本实例采用深度网络结构。

具体实施方式

相较于宏观唯象本构模型基于室内物理实验，本发明基于离散元数值试验得到具有不同结构特征的颗粒材料在不同加载路径下的应力应变数据，形成大规模的应力应变关系数据集，采用深度学习网络挖掘应力应变数据中潜藏的本构关系，得到基于深度学习和数据驱动的岩土颗粒材料本构模型。

本发明提供了一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其包括以下步骤：步骤1，首先进行离散元数值模拟的加载速度、时间步、颗粒集合体规模的敏感性分析，确定合适的模拟参数能再现岩土颗粒材料RVE尺度的力学特性；步骤2，离散元数值试验，通过改变级配的方式，得到具有不同结构特征的颗粒材料，根据给定的级配曲线生成若干个颗粒材料数值试样，对每个数值试样进行不同加载路径下的离散元数值模拟；步骤3，离散元数据处理，计算生成的离散元数据由文本格式导出，通过Python脚本读取不同路径下的数据文件，提取训练所需数据，计算岩土材料内变量，建立训练数据库；步骤4，按间隔从每一试样原始数据中抽取指定帧数构造深度学习数据集，对数据集进行对数化、归一化处理；步骤5，搭建深度学习网络，采用改进的长短期记忆神经网络(Long Short TermMemory,LSTM)为深度学习网络的主要组成部分，采用网格搜索方法进行深度学习网络超参数优化，确定网络结构；步骤6，采用早停策略防止深度学习网络过拟合，数据分批次输入网络进行训练，网络收敛后保存网络结构；步骤7，调用训练好的深度学习网络，输入颗粒材料的结构信息和时序应变数据，预测岩土颗粒材料的时序应力数据和结构参数。本方法的技术路线如图1所示。

以下结合附图对本发明涉及的基于深度学习和数据驱动的岩土颗粒材料本构建模方法的具体实施方案进行详细说明，技术路线如图1所示。

步骤1，首先进行离散元数值模拟的敏感性分析，确定离散元数值试验中加载速度、时间步、颗粒集合体中颗粒数目。

步骤1采用离散元数值模拟替代室内物理试验，具有以下特征：为保证深度学习网络具有较强的泛化能力和较高的预测精度，需要进行具有不同结构特征的颗粒材料数值试样在不同加载路径下的离散元数值试验，得到大量的应力应变关系数据，用于训练深度学习网络。建立宏观唯象本构模型主要依赖室内物理试验数据，而常规试验装置难以进行复杂应力路径、高围压、大尺寸的物理试验，且试验成本较高耗时费力。因此，使用离散元数值试验代替物理试验具有重要意义，既节省了人力物力，大幅度提高试验效率，又增强试验的可重复性和调试灵活性。

步骤1-1，进行离散元数值试验时，需要保证颗粒试样在准静态条件下加载，因此需要将加载速度控制在一定范围内。加载速度过低，会大幅增加计算成本，因此，需要找到合适的加载速度兼顾计算效率和计算精度。本方法采用离散元模拟加载速度为0.1m/s；

步骤1-2，对离散元数值试验的时间步大小进行控制，时间步过大将导致离散元数值试验中能量不守恒、模拟失真且得到数据波动(噪声)大，不利于深度学习网络训练，时间步太小则计算成本成倍提升，通过分析本实例采用时间步大小为5e^-7s；

步骤1-3，颗粒材料数值试样的规模越大，即颗粒个数越多，其在加载过程中的应力应变数据越稳定，同时计算成本大幅提升，通过颗粒个数为5000、10000、20000的三种数值试样的颗粒数目和数据波动的敏感性分析，最终确定离散元数值试验采用颗粒数目为10000。

步骤2，离散元数值试验，制备具有不同结构特性的颗粒材料数值试样，进行不同加载路径下的离散元数值试验。

步骤2进行具有不同结构特性的颗粒材料数值试样在不同加载路径下的离散元数值试验，具有以下特征：

步骤2-1，加载停止条件：由于轴向应变加载到0.3～0.35附近时候，试样即达到临界状态，临界状态如公式(1)所示，试样的静水压力、剪应力和体积不随剪应变的增加而发生变化，因此将数值试样的加载到轴向应变0.4。

上式中，ε_s为剪应变，p为静水压力，q为等效偏应力，ε_v为体积，计算公式如下所示。

为反映不同结构特征的岩土颗粒材料在多种加载路径下的应力应变响应，挖掘岩土颗粒材料的应力变化本构关系，训练数据应包含不同结构特征试样在不同加载路径下的数据。

步骤2-2，制备了具有不同级配的试样：本方法为研究具有不同结构特征的颗粒材料数值试样在加载过程中的应力应变关系，制备了具有不同级配的颗粒材料数值试样，不同颗粒材料的级配曲线如图2所示。

为了让深度学习网络能挖掘颗粒材料结构特征与宏观应力变形响应之间的关系，训练数据应该包括具有不同结构特征的颗粒材料数值试验的离散元数值试验结果。制备具有不同级配的颗粒材料数值试样：本方法为研究具有不同结构特征的数值试样在加载过程中的应力应变特性，采用了不同的颗粒级配，级配曲线由具有最小粒径截断的分形公式(3)所示。

上式中，F(d)表示累计质量分数，d表示颗粒粒径，d_max和d_min分别表示最大和最小颗粒粒径，β为颗粒体系结构描述参数。

步骤2-3，为了反映岩土颗粒材料在复杂加载路径下的应力变形特性，需要进行多种加载路径的离散元数值试验，以便深度学习网络能挖掘到岩土颗粒材料的应力应变本构关系。

加载路径设置：本方法准备了6种加载路径下的应力应变数据，等p等b加载、等σ₃等b、排水常规三轴压缩、排水常规三轴拉伸、不排水常规三轴压缩、不排水常规三轴拉伸。其中

为平均静水压力。σ₃为小主应力，b＝(σ₂-σ₃)/(σ₁-σ₃)为中主应力系数。

步骤2-4，加卸载路径设置：本方法准备了在常规三轴状态下的加卸载循环路径。

步骤2-5.加载初始围压、孔隙比设置：本方法设定10级加载围压，分别为0.01、0.05、0.1、0.2、0.5、1.0、2.0、4.0、8.0和16.0MPa；通过在制样阶段设置不同颗粒间摩擦系数来制备不同初始孔隙比的试样，相同围压下摩擦越小得到的试样孔隙比越小，设置制样摩擦系数分别为0.001、0.05、0.1、0.2、0.3。

步骤2-6，离散元数值试验，由于试验需要调用不同的试样在各种路径下进行加载试验，计算量庞大且计算复杂，通过编写Python脚本控制加载试验在超算集群上进行计算。

步骤3，离散元数据处理，计算生成的离散元数据由文本格式导出，通过Python脚本读取不同路径下的数据文件，提取训练所需数据，计算岩土材料内变量，建立训练数据库。

步骤3-1，读取的数据文本，转化成浮点数类型；

步骤3-2，提取DEM数据中的颗粒集合体宏观力学信息，包括三个方向应变、应力、集合体孔隙比、集合体颗粒体系结构信息、颗粒材料杨氏模量和摩擦系数，还包括累积绝对应变，该数据作为深度学习训练数据，；

根据应变数据，计算累积绝对应变，采用累积绝对应变作为岩土材料内变量输入，将岩土颗粒材料在变形过程中累积的绝对应变作为状态变量，用于表示材料当前状态的唯一性，具有以下特征：

岩土颗粒材料具有摩擦性、压硬性和剪胀性三个基本特征。由于岩土颗粒材料是典型的多体相互作用的耗散体系，其在荷载作用下的变形主要是由颗粒尺度重排列、接触滑移(摩擦耗散)和颗粒破碎引起，这三种颗粒尺度的力学行为均是不可逆的，因此岩土颗粒材料在加载过程中必然会产生不可逆的变形和结构变化。在弹塑性理论框架内，岩土颗粒材料在未经历过任何塑性变形情况下，首次进入屈服时应该满足屈服条件，此时应力状态在π平面上的投影落在初始屈服面上。材料的压硬性表现为，进入屈服后，随着塑性变形的增长屈服极限会提高，即产生硬化，通俗的说就是随着塑性变形的产生，材料的屈服面发生变化。在宏观唯象本构模型中，将材料的累积塑性应变或塑性功作为内变量，建立屈服面与该内变量的关系，然而在数值试验和物理试验中难以得到塑性应变和塑性功。应变增量由塑性应变和弹性应变组成，如下公式(4)所示，

dε＝dε^e+dε^p (4)

其中，dε表示应变增量，dε^e表示弹性应变增量，dε^p表示塑性应变增量。

利用深度学习网络的数据挖掘能力，将岩土颗粒材料在加载过程中累积的应变增量绝对值作为状态量，用来表示材料的当前状态，内变量如下公式(5)所示，

上式中，ε^m表示第m荷载步时，累积应变增量。

步骤3-3，采用奇异值分解方法完成数据降噪，降噪后的数据能够更好的应用于深度学习网络的训练，降低模型训练收敛难度；

采用基于奇异值分解的方法进行应力应变数据降噪，具有以下特征：离散元数值试验得到的数据中具有一定程度的“噪声”，这些“噪声”部分来源于加载过程中的“雪崩”效应，本方法建立的本构模型不考虑这些“噪声”，且该“噪声”会增大深度学习网络的训练难度，因此需要对数据进行降噪。

奇异值分解方法如下所述，例如一个m×n阶矩阵，其中所有元素属于实全域(可能存在复数)，则存在一个分解使得下式成立，

M＝UΛV (6)

其中，Λ为对角矩阵，U和V分别为m×m和n×n阶的复数方块矩阵，其具有以下性质

U^TU＝UU^T＝I (7)

U^-1＝U^T (8)

上式中，I表示m阶单位阵。

奇异值分解得到对角阵Λ，将对角阵从大到小排序，通过选取前k个对角元素，其他元素置为0，得到简化后的对角阵∑，采用以下公式恢复数据，

N＝UΣV (9)

矩阵N即为降噪后的数据，降噪效果实例展示如图3所示。

步骤4，采用时间窗口方法构建训练数据集，对训练数据集进行对数化和归一化处理。

步骤4-1，按间隔从每一试样原始数据中抽取指定帧数(如图8所示)，组合成时序数据序列构造深度学习数据集，如图8所示，非时序输入、时序输入和时序输出组合起来表示一个训练样本，网格中行数代表不同时间步，网格中列数代表不同输入输出；

步骤4-2，训练数据归一化，训练数据归一化有利于深度神经网络训练收敛。由于输入输出数据差异较大，其中应力数据甚至相差几个数量级，在进行归一化之前，需要对应力数据进行对数化处理，将对数化之后的应力数据进行归一化处理。

步骤5，深度网络模型搭建，深度学习模型主要由改进的LSTM网络单元和全连接层组合而成，网络搭建过程中涉及到多个超参数的优化，例如：网络层数、单层所含改进的LSTM单元个数、单层网络的激活函数选择(relu、sigmoid、tanh、linear等)、优化算法的比选等。其网络结构如图4所示。

本步骤设计初始化信息，时序输入、非时序输入和时序输出。初始化信息包括时序数据的第一个时间步的孔隙比、应力和应变，用于初始化细胞状态和隐含状态；时序输入包括应变、内变量；非时序输入包括杨氏模量、摩擦系数和颗粒体系结构描述参数β；时序输出为孔隙比和应力。初始信息和输入输出的具体情况，如附图7-9所示。

改进具体措施：tensorflow和keras两个库中未提供有效API来根据不同训练样本为经典LSTM单元设置初始细胞状态C₀和初始隐含状态h₀，因此本方法调用tensorflow的tf.Varible类完成LSTM单元结构的编写。时序训练数据输入单元前，首先定义t₀时刻的LSTM单元的细胞状态和隐含层结构。由于本方法中初始化信息(孔隙比e、应力σ_x,σ_y,σ_z、应变ε_x,ε_y,ε_z)维度为7，无法与LSTM单元隐含层的hidden_unit_num一致，因此在输入初始化信息前加入全连接层(Full connected layer)对信息进行预处理，如下公式所示。

C₀＝W_c·C₀+b_c (10)

h₀＝W_h·h₀+b_h (11)

上述公式中，W_c和W_h分别为细胞状态和隐含状态对应的权值矩阵，b_c和b_h分别为细胞状态和隐含状态对应的偏置矩阵。

在采用深度学习建立岩土颗粒材料本构模型的过程中，不仅涉及到时间序列数据的输入(即应力应变关系数据、时序相关的颗粒材料结构特征)，还包括了非时序数据的输入(初始孔隙率、围压、加载路径等)。岩土颗粒材料的应力变形特性不仅受加载历史的影响，还与材料的初始状态和结构特征有关。经典LSTM单元未对初始细胞状态C₀和初始隐含状态h₀进行定义，无法与岩土材料初始特征相关性结合，本方法对通过在LSTM单元中加入全连接层，并在输入时序数据前，将初始细胞状态C₀和初始隐含状态h₀定义为当前岩土材料的孔隙比、应力。改进后的LSTM单元更能够体现岩土材料的初始特性影响，预测精度显著提升。

改进的LSTM网络单元由初始状态、输入门、遗忘门和输出门组成，数据从进入单元到输出经历遗忘阶段、更新阶段和输出阶段，具体过程如下，

1.初始状态

输入初始状态信息，分别赋值给LSTM单元初始细胞状态C₀和初始隐含状态h₀，再分别经过全连接层的预处理，如公式(10-11)所示。

2.遗忘阶段

这一阶段“遗忘”多余信息，保留重要信息。LSTM单元首先通过遗忘门(一个Sigmoid单元)对上一步输入h_t-1和该步输入x_t进行筛选，然后输出一个0～1间的向量f_t，如公式(12)所示。

f_t＝σ(W_f·[h_t-1,x_t]+b_f) (12)

上式中，σ表示Sigmoid函数，W_f表示权值矩阵，b_f表示偏置值。

3.更新阶段

更新阶段综合历史信息和当前步输入信息更新当前步的细胞状态，让上一个节点输入的h_t-1和x_t通过输入门以此决定更新后的信息i_t。然后再让h_t-1和x_t信息通过一个tanh层得到新的候选细胞信息

这些信息可能会被更新到细胞信息中。然后旧的细胞信息C_t-1通过遗忘门选择忘记旧细胞信息的一部分，通过输入门选择添加候选细胞信息，变为新的细胞信息C_t，如公式(13)所示：

4.输出阶段

输出阶段综合更新后的细胞状态信息和该步输入信息计算输出。更新完细胞状态后需要根据输入的h_t-1和x_t来判断输出细胞的状态特征，因此需要将信息通过输出门得到判断条件O_t，然后将更新后的细胞状态C_t经过tanh层得到一个-1～1之间值的向量，该向量与输出门得到的判断条件相乘就得到了最终该LSTM单元的输出h_t，如公式(14)所示：

本方法采用的网络结构如图9。网络结构包含时序输入和非时序输入，其中时序输入包含岩土颗粒材料在加载过程中x，y，x三个方向的轴向应变

(上标i表示荷载步)和累积应变

(上标m表示荷载步)，累积应变充当岩土颗粒材料内变量的角色，其中，加载过程中，轴向变化的长度除以原始长度就是轴向应变，累积应变通过公式(5)获得。

非时序输入包括杨氏模量、摩擦系数和颗粒体系结构描述参数β。非时序输入不随加载进行发生变化，且对加载过程持续产生影响。因此，该信息每个时间步均输入给LSTM单元。

时序输出为材料的各向应力

和颗粒材料孔隙比时序数据。

本方法从计算得到的DEM数据中按时间步抽取出应力应变信息，如图7中第一张图所示，抽取三个方向应力、应变一级孔隙比等时序数据，颗粒集合体的级配信息、摩擦系数、杨氏模量作为非时序数据，如图7所示为时间序列数据构造方法。全部DEM数据集构造得到的训练数据集如图8所示，输入数据中包括时序数据和非时序数据，输出数据为应力和孔隙比两个时序数据。

本方法以网络收敛后的收敛误差为目标函数，通过网格搜索法寻找较优的参数组合。【注：此处使用“较优”是因为单次训练成本较大，采用的网格搜索法仅在能够接受的计算成本上进行搜索。本方法表明，“较优”的网络结构已经能够达到能够接受的验证精度】

步骤6，深度网络训练和保存，深度网络训练防过拟合策略。

步骤6采用早停策略训练深度学习网络，具有以下特征：

过拟合是指对于有限的数据而言，网络参数过多或者结构过于复杂，训练过程中网络已经开始不自觉从数据噪声中获取信息并表现在网络中，导致训练偏差不断减小而在验证集上表现出的误差逐渐增大，模型泛化能力变差。过拟合现象体现在数据和误差上如图5所示，(a)表示的为一个二分类，图中的曲线将两种不同的点完全分开，曲线已经捕捉噪声，而黑色线条将两种点主要部分分开，准确率达到95％甚至更高，更能代表样本的一般性，得到的模型泛化能力更强；(b)中训练的段半段时间，训练误差和验证误差均稳定下降，训练到一定程度之后，训练误差持续下降，但验证误差开始上升，此时网络开始从噪声中获取信息。

防止网络训练过拟合的措施有：模型选择、正则化、drop out、交叉验证和提前停止等。其中模型选择一般表示选择复杂程度与所需训练数据的复杂程度相当的模型，网络在训练得到数据信息之后，没有剩余容量来学***衡关系，实际操作难度大；正则化主要有L1和L2正则化方法，通过控制权值和偏置来控制模型的局部切线斜率，将该值控制在一个什么范围内依然存在问题；drop out方法与前两种方法类似，难以准确控制drop out的程度；因此本方法采用简单实用的提前停止策略，该策略理论简单，实施起来容易控制，当训练偏差持续下降，而验证误差已经开始趋于稳定或者逐渐增大时，备份上一步验证误差最小的模型，继续进行一定的训练尝试，尝试次数可格局需要设定。本方法设定尝试次数为10000次，若在训练10000次的过程中，得到更小的验证误差，则更新已尝试次数，备份最好的模型，继续训练；否则，保存上一次训练备份的模型，停止训练。

步骤7，深度网络恢复与调用，网络训练完成后将其保存为四个文件，分别为checkpoint文件、modelname_stop.data文件、modelname_.ckpt.index文件、modelname_.ckpt.meta文件，四个文件包含了网络结构、网络结构索引、网络权值偏置值等信息。

再次调用文件时首先恢复网络结构，再用过导入checkpoint文件，恢复网络权值偏置值，通过网络结构索引调用输入输出接口，进行不同结构特征的岩土颗粒材料多种复杂加载路径下的应力和孔隙比预测。

本发明的最终的预测效果如图6所示。

本发明能够较好地预测不同级配的颗粒集合体在各种加载路径下的应力以及孔隙比发展情况，该模型通过数据驱动的方法能够将颗粒材料的杨氏模量、摩擦角、级配特性和颗粒集合体的初始特性与颗粒材料的应力应变、孔隙比变化联系起来，较好地预测不同加载路径下应力和孔隙比的发展过程。

该方法训练得到的深度网络能够从数据中挖掘岩土材料的本质特性，再进行预测时将改特性体现在预测数据中。

图6可以看出，小应变的情况下，材料应力随着轴向应变线性增长，网络能够预测材料在弹性阶段的应力和孔隙比变化；材料加载进入塑性阶段后，预测曲线斜率减小，网络能够挖掘得到塑性屈服信息；达到峰值状态前，材料孔隙比轻微下降过程，网络能够模拟到材料的减缩特性；超过峰值状态后，材料孔隙比急剧增大，网络能够预测材料的剪胀特性；颗粒材料加载到轴向应变接近0.3时发生临界行为，与临界状态理论相吻合。图6中(d,e)可以看出，网络不仅能够预测材料加载过程的应力应变，同样能够预测材料卸载过程，卸载完再次加载时，峰值应力不能达到首次的峰值状态，证明该方法训练得到的深度学习网络具有岩土材料的记忆效应，能够通过材料之前的输入信息综合归纳集合体的加载历史，结合当前加载，预测材料应力和孔隙比。

如图6所示，该深度学习网络能够抓住颗粒材料的压硬性、减缩剪胀性、历史依赖性等本质特性，较好预测材料应力和孔隙比，基于深度学习的本构模型不同于传统本构模型引入众多唯象假设，通过DEM模拟数据即可训练得到。该模型具有替代传统本构模型能力，将该模型与有限元计算软件结合，能够发展得到基于深度学习的有限元计算方法。该方法将宏观尺度的有限元模拟与细观尺度的颗粒材料结构紧密连接，有效提高有限元计算精度，

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其特征在于，包括如下步骤：

步骤1，首先进行离散元数值模拟的敏感性分析，确定离散元数值试验中加载速度、时间步、颗粒集合体中颗粒数目；

步骤2，制备具有不同结构特性的颗粒材料数值试样，进行不同加载路径下的离散元数值试验，获得不同结构特征的数值试样在不同加载路径下的应力应变特性数据；

步骤3，整理不同数值试样在不同加载路径得到的应力应变数据，提取所需数据，包括三个方向应变、应力、颗粒材料孔隙比、颗粒体系结构描述参数、颗粒材料杨氏模量和摩擦系数，以及累积绝对应变，对提取的所需数据进行降噪处理，建立原始训练数据库；

步骤4，按间隔从每一试样原始数据中抽取指定帧数，组合成时序数据序列构造深度学习训练集，非时序输入、时序输入和时序输出组合起来表示一个训练样本，对训练集进行对数化和归一化处理；

步骤5，搭建深度学习网络模型，深度学习网络模型由改进的LSTM网络单元和全连接层组合而成；

深度网络模型的输入包含时序输入和非时序输入，其中时序输入包含岩土颗粒材料在加载过程中x，y，x三个方向的轴向应变和累积绝对应变；输入时序数据前，将初始细胞状态C₀和初始隐含状态h₀定义为当前岩土颗粒材料的孔隙比、应力，并在输入初始化信息前加入全连接层，对初始化信息进行预处理，初始化信息包括时序数据的第一个时间步的孔隙比、应力和应变，用于初始化细胞状态和隐含状态；

非时序输入包括杨氏模量、摩擦系数和颗粒体系结构描述参数β，所述非时序输入不随加载进行发生变化，且对加载过程持续产生影响，因此，非时序输入每个时间步均输入给LSTM单元；

时序输出为材料的各向应力

和颗粒材料孔隙比时序数据，上标i表示荷载步；

步骤6，采用早停策略防止深度学习网络过拟合，将训练数据分批次输入深度网络模型进行训练，网络收敛后保存网络结构；

步骤7，调用训练好的深度学习网络模型，输入颗粒材料的结构信息和时序应变数据，预测岩土颗粒材料的时序应力数据和结构参数。

2.如权利要求1所述的一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其特征在于：步骤1中确定模拟加载速度为0.1m/s，时间步大小为5e^-7s，颗粒数目为10000。

3.如权利要求1所述的一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其特征在于：步骤2中制备具有不同结构特性的颗粒材料数值试样，即制备具有不同级配的颗粒材料数值试样，级配曲线由具有最小粒径截断的分形公式(3)确定，

上式中，F(d)表示累计质量分数，d表示颗粒粒径，d_max和d_min分别表示最大和最小颗粒粒径，β为颗粒体系结构描述参数。

4.如权利要求1所述的一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其特征在于：步骤2中进行不同加载路径下的离散元数值试验的具体实现方式如下；

加载路径设置：等p等b加载、等σ₃等b、排水常规三轴压缩、排水常规三轴拉伸、不排水常规三轴压缩、不排水常规三轴拉伸，其中

为平均静水压力，σ₃为小主应力，b＝(σ₂-σ₃)/(σ₁-σ₃)为中主应力系数；

加卸载路径设置：设置了在常规三轴状态下的加卸载循环路径；

围压、孔隙比设置：设定10级加载围压，分别为0.01、0.05、0.1、0.2、0.5、1.0、2.0、4.0、8.0和16.0MPa；通过在制样阶段设置不同颗粒间摩擦来制备具有不同初始孔隙比的数值试样，相同围压下摩擦系数越小得到的数值试样孔隙比越小，设定制样摩擦系数分别为0.001、0.05、0.1、0.2、0.3。

5.如权利要求1所述的一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其特征在于：步骤3中累积绝对应变，用于表示岩土颗粒材料当前状态的唯一性，计算方式如下，

应变增量由塑性应变和弹性应变组成，如下公式(4)所示，

dε＝dε^e+dε^p (4)

其中，dε表示应变增量，dε^e表示弹性应变增量，dε^p表示塑性应变增量；

累积绝对应变的计算公式如式(5)所示，

上式中，ε^m表示第m荷载步时，累积应变增量。

6.如权利要求1所述的一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其特征在于：步骤3中基于奇异值分解对应力应变数据进行降噪，具体实现方法如下，

一个m×n阶矩阵，其中所有元素属于实全域，则存在一个分解使得下式成立，

M＝UΛV (6)

其中，Λ为对角矩阵，U和V分别为m×m和n×n阶的复数方块矩阵，其具有以下性质

U^TU＝UU^T＝I (7)

U^-1＝U^T (8)

上式中，I表示m阶单位阵；

奇异值分解得到对角阵Λ，将对角阵从大到小排序，通过选取前k个对角元素，其他元素置为0，得到简化后的对角阵∑，采用以下公式恢复数据，

N＝UΣV (9)

矩阵N即为降噪后的数据。

7.如权利要求1所述的一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法，其特征在于：步骤5中改进的LSTM网络单元由初始状态、输入门、遗忘门和输出门组成，数据从进入单元到输出经历遗忘阶段、更新阶段和输出阶段，具体过程如下，

1.初始状态

时序训练数据输入单元前，首先定义t₀时刻的LSTM单元的细胞状态和隐含层结构，由于初始化信息：孔隙比e、应力σ_x,σ_y,σ_z、应变ε_x,ε_y,ε_z维度为7，无法与LSTM单元隐含层的hidden_unit_num一致，因此在输入初始化信息前加入全连接层对信息进行预处理，如下公式所示，

C₀＝W_c·C₀+b_c (10)

h₀＝W_h·h₀+b_h (11)

上述公式中，W_c和W_h分别为细胞状态和隐含状态对应的权值矩阵，b_c和b_h分别为细胞状态和隐含状态对应的偏置矩阵；

2.遗忘阶段

这一阶段“遗忘”多余信息，保留重要信息，LSTM单元首先通过遗忘门,对上一步输入h_t-1和该步输入x_t进行筛选，然后输出一个0～1间的向量f_t，如公式(12)所示，

f_t＝σ(W_f·[h_t-1,x_t]+b_f) (12)

上式中，σ表示Sigmoid函数，W_f表示权值矩阵，b_f表示偏置值；

3.更新阶段

更新阶段综合历史信息和当前步输入信息更新当前步的细胞状态，让上一个节点输入的h_t-1和x_t通过输入门以此决定更新后的信息i_t；然后再让h_t-1和x_t信息通过一个tanh层得到新的候选细胞信息

旧的细胞信息C_t-1通过遗忘门选择忘记旧细胞信息的一部分，通过输入门选择添加候选细胞信息，变为新的细胞信息C_t，如公式(13)所示：

4.输出阶段

输出阶段综合更新后的细胞状态信息和该步输入信息计算输出，更新完细胞状态后需要根据输入的h_t-1和x_t来判断输出细胞的状态特征，因此需要将信息通过输出门得到判断条件O_t，然后将更新后的细胞状态C_t经过tanh层得到一个-1～1之间值的向量，该向量与输出门得到的判断条件相乘就得到了最终该LSTM单元的输出h_t，如公式(14)所示：

Images