基于物理信息神经网络的煤体黏滑冲击动力学分析

王民华，闫永敢

（1.山西能源学院 矿业工程系，山西 太原 030006；2.太原理工大学 矿业工程学院，山西 太原 030024）

冲击地压是采掘工作面煤岩体积聚的弹性变形能突然释放，产生强烈震动，造成煤岩体剧烈破坏的动力灾害[1]。煤系地层为层状结构，并存在大量强度相对较低的结构面，在煤岩失稳时，煤层会沿顶底板之间产生水平方向的位移，煤岩结构系统会因摩擦滑动造成失稳而产生冲击地压[2]。煤岩体的主要破坏方式就是结构块沿结构面的剪切滑动[3]，煤岩体的摩擦滑动主要分为稳定滑动和黏滑，其中黏滑是指滑动过程中滑移面上的剪应力不断出现急剧增大或减小的过程[4]。近年来国内相关学者通过煤岩体的黏滑失稳机理研究，来客观认识冲击地压现象，认为冲击地压是煤岩体结构摩擦滑动破坏的一种形式，表现为瞬时的黏滑失稳过程。齐庆新等[5]开展了煤岩摩擦滑动实验研究，并通过煤岩体的黏滑理论对冲击地压发生机理进行了分析；梁冰等[6]认为矿震是煤岩体黏滑摩擦失稳的现象，从而提出了矿震的黏滑失稳理论；潘一山等[7]通过建立黏滑失稳模型，解释了断层冲击地压的间歇性；郭德勇等[8]设计了冲击地压与突出过程摩擦滑动模拟实验系统，进行了摩擦滑动的作用机理模拟研究，得到了黏滑过程中的振动及流变理论解释；尹光志等[9]用双状态变量模型模拟了冲击地压的黏滑特性，对冲击地压系统的黏滑特性和冲击地压系统的动力失稳过程进行了研究；黄滚等[10]通过建立黏滑失稳的双滑块模型，进行了煤岩体黏滑失稳机理的研究，并进行了煤岩体黏滑非线性混沌特性的分析；任晓龙[11]进行了不同静荷载水平和不同扰动荷载幅值作用下的超低摩擦现象的研究；姜耀东等[12]实验研究了煤岩组合体特定条件下的滑动类型；邓志辉等[13]通过实验研究，总结了断层的黏滑位移是台阶式跳跃运动规律；曹彦彦[14]通过系统的实验研究，并结合数值结果的分析，定量总结了黏滑过程中的损伤演化和裂纹扩展规律，进一步加深了对岩石结构黏滑失稳过程的认识。通过对煤岩体黏滑特性和滑动实验的大量研究，目前煤岩体黏滑冲击的物理机理已较为清晰，但理论成果应用到实际冲击地压防治和危险性评价中还较少，主要是基于黏滑冲击分析较为繁琐，且数值分析手段不易快速实现。

目前对于矿山灾害的数据驱动智能预测预警的研究较多，但能在实际生产中进行有效实际运用的数据智能控制模型，鲜有文献报道。目前人工智能驱动的科学研究（AI for science）已经成为工程数值计算的新范式，目前人工智能驱动的科学研究（AI for science）领域中，物理信息神经网络（physics－informed neural networks，PINN）是研究的热点，物理信息神经网络（PINN）是一种科学的机器学习技术，可以进行偏微分方程的求解。物理信息神经网络（PINN）是通过进行损失函数优化问题进行偏微分方程解的逼近，其工作原理是将数学模型集成到网络中，并用控制方程的残差项来强化损失函数，该残差项作为惩罚项来限制可接受解的空间，其可以被认为是一种无监督策略，不需要标记数据[15]。其中RAISSI 等[16]提出用于解决涉及非线性偏微分方程的正反问题的深度学习框架以来，物理信息神经网络（PINN）已经成为建模和科学计算的研究热点。在进行工作面冲击危险性评价过程中，经常遇到相关指数难于合理确定，有关物理机理难以运用到精准科学的评价工作中。基于物理驱动深度学习建模的PDE 求解新方法，具有在不使用任何标记数据（监测数据）的情况下快速解决计算力学问题的优势。为此，以PINN 计算框架为基础，研究使用较小计算代价，在不失物理机理的基础上，分析煤体黏滑冲击动力学特性，为快速动态的冲击地压危险性分析方法提供一种技术手段，以期在生产中为冲击地压危险性评价提供参考。

1 冲击地压的黏滑物理机理

1.1 工作面围岩极限平衡与支承压力分布

煤层采出后在围岩应力重新分布的范围内作用在煤层、岩层和矸石上的垂直压力称为“支承压力”[17]。准确地计算预计采场支承压力分布规律及大小对于冲击地压动力灾害预防具有重要的工程应用价值，尤其确定弹性区长度和支承压力分布大小是进行煤体黏滑冲击动力学分析的基础。

回采工作面前后的支承压力状态，可以划分为应力降低区、应力增高区、应力不变区，并且把工作面前方支承压力峰值到煤壁为极限平衡区，也称为塑性区[18]（X0），由峰值向煤体内侧为弹性区（X1）。计算极限平衡区宽度示意图如图1。

图1 计算极限平衡区宽度示意图Fig.1 Schematic diagram of calculating the width of the limit equilibrium zone

通过建立支承压力分布的平衡微分方程，可进行极限平衡区的计算，极限平衡区的宽度X0为[18]：

式中：λ 为侧压系数；c0、φ 为煤层与顶底板岩石交界面的黏聚力与内摩擦角；pa为支架对煤帮的支护阻力；K 为应力集中系数；ρr为上覆岩层的平均密度；g 为重力加速度；H 为开采深度；M 为开采厚度。

支承压力的影响区域即支承压力长度包括塑性区长度和弹性区长度，塑性区宽度按照式（1）进行计算，支承压力的影响区的长度和支承压力分布大小按照文献［19］进行确定，这样确定支承压力分布函数形式为线性函数，但又不失较为严谨的物理机理，对构建煤体黏滑冲击动力学控制方程，进行物理信息神经网络的无监督训练有较好的收敛性。支承压力线性分布示意图如图2。

图2 支承压力线性分布示意图Fig.2 Linear distribution of abutment pressure

支承压力采用简化的线性分布，支承压力分布方程为[19]：

式中：Sx为采场推进到L0后的支承压力分布范围，m；L0为工作面长度，m；cx为采场推进距离,m；bx为采场推进到L0后的由综合移动角所形成范围，m，bx=H·cot θ；θ 为综合移动角，（°）；p 为弹性区与塑性区分布的比值；Kmax为压力峰值与原始应力场比值。

煤体黏滑冲击动力学分析的关键是进行弹性区煤体动力特性的表征，进行煤体黏滑冲击动力学物理驱动求解的前提需要得到弹性区支承压力的分布规律。通过式（2）和式（3）可得到弹性区长度X1，即：

1.2 煤体黏滑动力学控制方程

通过双剪摩擦实验和三轴摩擦实验，表明了砂岩-煤试件最易产生黏滑，且组合煤岩体的摩擦滑动特性表现为典型的黏滑过程，因此可以认为黏滑失稳破坏是冲击地压的一种特殊形式[5]。

煤体黏滑最可能发生方式为煤体与顶底之间黏滑，尤其是在工作面超前支承压力的作用下，工作面前方煤体出现应力集中，这种顶底板-煤体的组合煤岩体摩擦滑动特性表现为典型的黏滑过程，在一定条件下会导致黏滑失稳导致冲击地压的发生。在超前支承压力的作用下，采煤工作面前方煤体会产生塑性区X0、弹性区X1及原岩应力区。在超前支承压力的作用下，弹性区积聚了大量的弹性能，对于黏滑失稳导致的冲击地压现象，进行弹性区黏滑失稳的动力学特性研究有重要的意义。弹性区黏滑过程的分析示意图如图3。

图3 煤体黏滑冲击动力学分析模型Fig.3 Dynamic analysis model of coal stick-slip impact

按图3 关系可建立如下煤体动力学方程[20]：

式中：f 为动摩擦系数；σx为水平应力；σy为垂直应力；ρc为煤体密度，kg/m3；u（x，t）为煤体动力学方程的解；x 为弹性区水平位置坐标；t 为时间，s。

引入胡克定律:

式中：ε 为轴向应变；E 为弹性模量。

则式（6）变为：

由本文采用的工作面支承压力线性分布，可得在弹性分布区域（X1），垂直应力为：

式（9）代入式（8）可得煤体黏滑动力学偏微分方程：

由式（10）可知，煤体黏滑动力学系统特性与主要与煤弹性模量、煤密度、煤层厚度、动摩擦系数、弹性区域长度、上覆岩层密度、超前应力集中系数、煤层埋深有关。

2 煤体黏滑冲击物理驱动分析原理

2.1 PINN 神经网络架构

基于物理驱动的深度学习（PINN）是近年来开发的一种用于建模PDE 解决方案的新方法，并显示出在不使用任何标记数据（例如测量数据不可用）的情况下解决计算力学问题的前景[21]。

物理信息神经网络（PINN）结合了物理知识和神经网络的方法，可用于解决偏微分方程（组）问题。相比于传统的PDE 求解手段（有限差分、有限元和谱方法等），PINN 无需离散化空间和时间，使用神经网络来逼近未知函数，并使用物理定律作为先验知识来保证解的合理性，具有高效，无网格计算等优势。RAISSI 考虑通过训练全连接神经网络来逼近偏微分方程的解，其具体方法如下所述。

假设非线性偏微分方程：

式中：u（x，t）为偏微分方法的解；N［u］为微分算子。

定义f（t，x）由式（11）的左边给出[15]，通过深度神经网络进行逼近u（t，x），并产生了1 个物理信息的神经网络f（t，x），该网络可以通过应用链式法则来推导，使用自动微分对函数的组成进行微分，并且具有与表示u（t，x）的网络相同的参数。尽管由于微分算子N 的作用而具有不同的激活函数，神经网络u（t，x）和f（t，x）之间的共享参数可以通过最小化均方误差损失来学习。

通过训练神经网络进行逼近，定义损失函数MSE 为：

式中：MSEu为边界和初始条件损失项；MSEf为偏微分方程残差损失项；Nu为边界和初始条件的个数；Nf为内部残差项的个数；ui为目标函数；（tu，xu）为边界上的点；（tf，xf）为内部的配点。

如果得到网络参数能使控制方程的MSE 趋于0，那么就得到了PDE 的逼近值，即可以认为所求结果趋于真实值。这样就把求解偏微分方程的问题转化为了如何优化损失函数的问题。

2.2 煤体黏滑冲击动力学分析的PINN 模型

通过物理神经网络（PINN）对煤体黏滑冲击动力学控制方程进行求解，建立基于PINN 煤体黏滑冲击动力学分析模型，来研究煤体黏滑冲击动力特性。PINN 方法中将偏微分方程边界条件和初始条件的总和作为神经网络的损失函数以达到求解目的。物理信息神经网络（PINN）通过使用自动微分将煤体黏滑冲击动力学控制偏微分方程（式（10））嵌入到神经网络的损失函数中。

煤体动力学控制方程损失函数包括为边界和初始条件损失项和偏微分方程残差损失项，其中偏微分方程残差损失项为：

式中：（xi，ti）、（xj，tj）分别是在初始与边界位置和整个定义域采样的2 组点。

设从x=0 到x=X1的一段煤体(弹性区煤体)发生黏滑，两端固定，发生黏滑前煤体处于静止状态，假设黏滑前煤体位移为0，则边界条件为：

通过基于梯度的优化器的最小化损失来训练网络，直到损失小于阈值，由此构造1 个神经网络u（x，t，θ）。其中：θ 为可训练权重w 和偏差b 的集合；σ 为非线性激活函数；（xi，ti，ui）为u 指定数据；（xj，tj）为PDE 指定残差点。通过将数据和PDE 的加权损失相加，指定损失L，训练NN 通过最小化损失L 来找到最佳参数θ*。

3 算例分析

通过给定初始条件的算例，进行煤体黏滑冲击动力学分析的PINN 模型训练，并验证其求解结果的有效精准性。假设煤弹性模量为5 GPa，煤密度为1.1 g/cm3，煤层厚度2 m，动摩擦系数0.15，弹性区域长度X1为30 m，上覆岩层密度2.5 t/m3，超前应力集中系数为1.5，煤层埋深400 m。边界条件如式（17）。建立的PINN 求解算例如图4。

图4 算例示意图Fig.4 Schematic diagram of calculation example

在传统神经网络中，输入量常为被预测量的相关影响因素，而正向PINN 网络是1 种无监督学习，在研究中，将求解域坐标点（x，t）作为模型输入，Net作为网络中间变量输出，u 作为网络最终输出，从而形成网络的映射关系，具体过程可表述为：

式中：X 为输入变量，为包含x，t 坐标点的张量；σ 为激活函数，采用ELU 激活函数；W、b 分别为网络权重和偏置；D（x，t）为距离函数，满足当x，t 位于定解条件坐标时，函数值为0；BC、IC 分别为边界条件和初始条件。

由上述可知，一旦满足相应坐标，网络输出强制性满足边界条件或初始条件，由于正向PINN 是一种无监督学习求解器，因此需要构造目标函数以满足模型求解需要。根据偏微分方程形式，可以构造1个函数f（x，t）：

如果f（x，t）趋向于0，则原偏微分方程满足，根据f（x，t）以及神经网络损失函数性质，可将f（x，t）作为PINN 的损失函数，f（x，t）中的偏导项由自动微分技术完成。

整个求解网络基于Pytorch 1.10 环境编写，求解GPU 为Tesla V100 32GB，坐标采样点为30×30，网络层输入节点为2，隐含层结构为64×4，输出层节点为1。为了验证PINN 求解的正确性，研究基于相同工况采用显示有限差分法进行了对比求解，物理信息神经网络与显式有限差分求解结果对比如图5，物理信息神经网络与显式有限差分局部求解结果对比曲线（分别为x=5、15、20 m 剖面下2 种模型对比）如图6。为了反映局部求解性能，依次选择x=5、15、20 剖面下的求解结果进行对比，模型损失下降曲线如图7。

图5 物理信息神经网络与显式有限差分求解结果对比Fig.5 Comparison between physical information neural network and explicit finite difference solution

图6 物理信息神经网络与显式有限差分局部求解结果对比曲线（分别为x=5、15、20 m 剖面下2 种模型对比）Fig.6 Comparison curves of local solution results（comparison of two models under x＝5, 15 and 20 m sections）

图7 模型损失下降曲线Fig.7 Model loss decline curve

由图5、图6 可知：PINN 求解结果与显示有限差分求解结果基本相同，PINN 模型能够很好地捕捉形变过程，并且2 个模型求解结果的量级几乎一致，求解均方根误差低至1.198 9，说明PINN 可以很好地应用于冲击地压的动力学方程求解。

图6 反映了PINN 在局部求解性能，通过对比曲线可知，2 个模型的求解变化趋势一致，但在谷值处仍然存在误差，考虑到有限差分法是一种近似数值求解，无法通过局部些许差异判断模型的误差。所以对PINN 网络的损失曲线进行分析。

图7 展示了PINN 模型求解过程的模型损失下降过程，在迭代后期，损失量级达到了0.001，说明此时模型基本处于收敛状态。结合研究中的PINN损失函数理论，模型损失可以视为对偏微分方程及定界条件的拟合程度，图7 中极具收敛的损失可以认为模型已经对冲击地压的动力学方程完成了求解。

4 深度学习物理驱动的煤体黏滑冲击动力特性分析

通过算例分析可知，物理神经网络可以完成对煤体冲击动力学控制方程的求解，并且具有较强逼近控制方程真解的能力。由煤体冲击动力学控制方程所求水平位移u（x，t），即为研究区域的岩体水平变形，通过一定的工程研究背景，研究采煤工作面超前支承弹性区水平变形演化过程，来分析煤体黏滑冲击动力特性。

4.1 矿井概况

晋华宫矿12 号煤层埋深377～415 m，12 号煤层顶板岩性主要为中粒砂岩，局部为砂质泥岩、细粒砂岩，底板岩性为细粒砂岩、砂质泥岩。12 号煤层冲击倾向性综合判定结果为弱冲击倾向性，弹性模量5.9 GPa，煤体密度为1.41 g/cm3，抗压强度13.99 MPa，煤层厚度2.65 m，动摩擦系数为0.15 应力集中系数为1.5。且所研究工作面推进长度大于工作面的长度，X1长度可以按照式（5）进行计算，经过计算并结合现场支承压力实测结果，取值为50 m。

4.2 煤体黏滑动态演化分析

以晋华宫矿12 号煤层回采工作面为研究背景，运用物理神经网络进行煤体冲击动力学控制方程的求解，共训练50 000 步。K=1.5 位移场预测如图8。K=1.5 局部求解结果对比曲线如图9。

图8 K=1.5 位移场预测Fig.8 Displacement field prediction of K=1.5

图9 K＝1.5 局部求解结果对比曲线Fig.9 Comparison curves of local solution results of K＝1.5

通过图8 可知：煤体弹性区黏滑具有明显的波动性，与文献［5，14］实验室进行煤岩组合结构和岩体结构的黏滑试验结果相吻合。并且煤体变形过程出现时间间隔性，变形随着时间可以恢复，进一步说明了煤体处在弹性变形阶段。由于控制方程是基于胡克定律使用位移（变形）的微分替代的水平应力微分，所以求解得到的位移（变形）的变化，也能表明水平应力也存在波动变化性，并且间接表明了煤体黏滑过程是区域内水平应力瞬时的松弛（减小）和增加的过程，煤体黏滑的这种演化过程与冲击地压发生时出现的顶底板瞬时加卸载及工作面煤层出现瞬时破坏的冲击地压孕育过程是一致的。

由图9 可知：不同时刻煤体黏滑区域的位移（变形）大小变化较大，但不同时刻变形最大值所处位置基本相同，即煤体黏滑过程，随着时间的推移水平应力的大小在弹性区内存在一定的分布规律, 不同时刻在距离弹性区起始位置25 m 左右均出现形变最大值。随着时间的推移黏滑区域的位移大小出现波动性。

采场超前支承压力的变化是工作面煤体冲击的重要前兆信息，也是冲击地压监测的主要指标之一，所以研究不同支承压力作用下的煤体黏滑冲击动力特性，对工作面煤体冲击灾害的防治有重要的现实意义。以晋华宫矿12 号煤层回采工作面为研究背景，煤体黏滑冲击动力控制方程其它参数保持不变，只改变应力集中系数，分别进行应力集中系数K=2、K=3 情况下的基于物理神经网络的煤体黏滑冲击动力控制方程求解。

K=2 位移场预测如图10，K=2 局部求解结果对比曲线如图11。

图10 K=2 位移场预测Fig.10 Displacement field prediction of K＝2 at time

图11 K=2 局部求解结果对比曲线Fig.11 Comparison curves of local solution results of K＝2

由图10 可知：当采场支承压力集中系数为2时，煤体黏滑位移（形变）的最大值有明显的增加，尤其初期煤体黏滑位移（形变）最大值比集中系数为1.5 时增加将近1 倍，并且水平应力瞬时的松弛（减小）和增加的过程数值变化幅度较大，初期的2 次位移（形变）的高峰出现明显的大幅度变化间隔。

由图11 可知：在不同时刻沿x 轴22 m 位置截面左右均出现形变最大值，在计算时间域内位移均出现波动性。

K=3 位移场预测如图12，K=3 局部求解结果对比曲线如图13。

图12 K=3 位移场预测Fig.12 Displacement field prediction of K=3

图13 K=3 局部求解结果对比曲线Fig.13 Comparison curves of local solution results of K=3

由图12 可知：当采场支承压力集中系数为3时，采场煤体发生典型的黏滑现象，位移出现准周期性的变化。初期出现位移的急剧增加，然后瞬时松弛（减小），经过短时间的稳定，又进入急剧增加，然后瞬时松弛（减小）的周期过程。每个周期间隔过程都出现位移无法完全恢复的位移量，且随着黏滑周期性的波动，无法恢复的位移量逐渐增大。随着煤岩系统黏滑多周期性的波动演化，系统储存弹性能逐渐增大以及煤岩黏滑过程的强度损伤，在一定的采动影响下，有可能会孕育煤层顶底板瞬时加卸载破坏的冲击地压灾害。

由图13 可知：煤体沿x 轴各截面形变周期性明显，且周期间隔的形变量在逐步增加，在不同时刻沿x 轴24 m 位置截面左右均出现形变最大值。

地质构造异常区、断层带附近易发生冲击地压和其所处煤岩体摩擦系数的变化有较大关系[5]，通过改变动摩擦系数的大小进行煤体黏滑冲击动力控制方程的求解，来研究动摩擦系数对煤体黏滑冲击动力特性的影响。f=0.16 位移场预测如图14，f=0.20 位移场预测如图15。

图14 f=0.16 位移场预测Fig.14 Displacement field prediction of f=0.16

图15 f=0.20 位移场预测Fig.15 Displacement field prediction of f=0.20

由图14 可知：当动摩擦因数f 增加0.01 为0.16时，系统位移场数值变化明显，所以动摩擦因数对系统动力学特性影响较大。

由图15 可知：当动摩擦因数为2 时，位移场相对最大位移减少了20 mm，比原位移场最大值减少了50%，且通过位移场预测结果可知系统黏滑波动性已经不明显，随着时间推移系统较为稳定。

由文献［5］可知砂岩-煤试样的摩擦因数在0.172 8～0.275 3 之间，数值变化范围较大，煤岩结构在地质构造异常区、断层带附近动摩擦因数急剧减小，极易造成黏滑失稳，也为在地质构造异常区、断层带附近易发生冲击地压找到了摩擦滑动失稳的数值模拟解释。

5 结 语

结合煤体黏滑冲击理论和物理神经网络方法，在无任何观测标签数据的情况下建立了基于物理驱动深度学习的煤体黏滑冲击分析方法，并通过与显式有限差分求解方法进行求解结果对比，证明了其求解的有效准确性，并以晋华宫矿12 号煤层工作面为研究背景，进行了煤体黏滑冲击动力学分析。

1）物理信息神经网络（PINN）能够进行煤体黏滑冲击控制方程的数值模拟求解，并且使用较小的计算代价能够得到较为准确有效的数值模拟结果，实现了无网格和无标签数据情况下的煤体黏滑冲击数值模拟。

2）煤体黏滑变形过程出现时间间隔性，水平位移（变形）和水平应力均存在波动变化性；煤体黏滑过程是区域内水平应力瞬时的松弛（减小）和增加的过程，其与冲击地压发生时出现的顶底板瞬时加卸载及工作面煤层出现瞬时破坏一致。

3）当采场支承压力集中系数为2 时，煤体黏滑位移（形变）的最大值有明显的增加，尤其初期煤体黏滑位移（形变）最大值比集中系数为1.5 时增加将近1 倍，并且水平应力瞬时的松弛（减小）和增加的过程数值变化幅度较大；当采场支承压力集中系数为3 时，采场煤体发生典型的黏滑现象，位移出现准周期性的变化，初期出现位移的急剧增加，煤岩系统出现黏滑多周期性的波动演化，系统储存弹性能逐渐增大以及出现煤岩黏滑过程的强度损伤。

4）动摩擦因数对煤体黏滑动力学特性影响较大，当动摩擦因数变化时，煤体黏滑系统位移场数值变化明显。