基于物理信息神經(jīng)網(wǎng)絡(luò)的煤體黏滑沖擊動力學(xué)分析

王民華，閆永敢

（1.山西能源學(xué)院 礦業(yè)工程系，山西 太原 030006；2.太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原 030024）

沖擊地壓是采掘工作面煤巖體積聚的彈性變形能突然釋放，產(chǎn)生強(qiáng)烈震動，造成煤巖體劇烈破壞的動力災(zāi)害[1]。煤系地層為層狀結(jié)構(gòu)，并存在大量強(qiáng)度相對較低的結(jié)構(gòu)面，在煤巖失穩(wěn)時，煤層會沿頂?shù)装逯g產(chǎn)生水平方向的位移，煤巖結(jié)構(gòu)系統(tǒng)會因摩擦滑動造成失穩(wěn)而產(chǎn)生沖擊地壓[2]。煤巖體的主要破壞方式就是結(jié)構(gòu)塊沿結(jié)構(gòu)面的剪切滑動[3]，煤巖體的摩擦滑動主要分為穩(wěn)定滑動和黏滑，其中黏滑是指滑動過程中滑移面上的剪應(yīng)力不斷出現(xiàn)急劇增大或減小的過程[4]。近年來國內(nèi)相關(guān)學(xué)者通過煤巖體的黏滑失穩(wěn)機(jī)理研究，來客觀認(rèn)識沖擊地壓現(xiàn)象，認(rèn)為沖擊地壓是煤巖體結(jié)構(gòu)摩擦滑動破壞的一種形式，表現(xiàn)為瞬時的黏滑失穩(wěn)過程。齊慶新等[5]開展了煤巖摩擦滑動實驗研究，并通過煤巖體的黏滑理論對沖擊地壓發(fā)生機(jī)理進(jìn)行了分析；梁冰等[6]認(rèn)為礦震是煤巖體黏滑摩擦失穩(wěn)的現(xiàn)象，從而提出了礦震的黏滑失穩(wěn)理論；潘一山等[7]通過建立黏滑失穩(wěn)模型，解釋了斷層沖擊地壓的間歇性；郭德勇等[8]設(shè)計了沖擊地壓與突出過程摩擦滑動模擬實驗系統(tǒng)，進(jìn)行了摩擦滑動的作用機(jī)理模擬研究，得到了黏滑過程中的振動及流變理論解釋；尹光志等[9]用雙狀態(tài)變量模型模擬了沖擊地壓的黏滑特性，對沖擊地壓系統(tǒng)的黏滑特性和沖擊地壓系統(tǒng)的動力失穩(wěn)過程進(jìn)行了研究；黃滾等[10]通過建立黏滑失穩(wěn)的雙滑塊模型，進(jìn)行了煤巖體黏滑失穩(wěn)機(jī)理的研究，并進(jìn)行了煤巖體黏滑非線性混沌特性的分析；任曉龍[11]進(jìn)行了不同靜荷載水平和不同擾動荷載幅值作用下的超低摩擦現(xiàn)象的研究；姜耀東等[12]實驗研究了煤巖組合體特定條件下的滑動類型；鄧志輝等[13]通過實驗研究，總結(jié)了斷層的黏滑位移是臺階式跳躍運動規(guī)律；曹彥彥[14]通過系統(tǒng)的實驗研究，并結(jié)合數(shù)值結(jié)果的分析，定量總結(jié)了黏滑過程中的損傷演化和裂紋擴(kuò)展規(guī)律，進(jìn)一步加深了對巖石結(jié)構(gòu)黏滑失穩(wěn)過程的認(rèn)識。通過對煤巖體黏滑特性和滑動實驗的大量研究，目前煤巖體黏滑沖擊的物理機(jī)理已較為清晰，但理論成果應(yīng)用到實際沖擊地壓防治和危險性評價中還較少，主要是基于黏滑沖擊分析較為繁瑣，且數(shù)值分析手段不易快速實現(xiàn)。

目前對于礦山災(zāi)害的數(shù)據(jù)驅(qū)動智能預(yù)測預(yù)警的研究較多，但能在實際生產(chǎn)中進(jìn)行有效實際運用的數(shù)據(jù)智能控制模型，鮮有文獻(xiàn)報道。目前人工智能驅(qū)動的科學(xué)研究（AI for science）已經(jīng)成為工程數(shù)值計算的新范式，目前人工智能驅(qū)動的科學(xué)研究（AI for science）領(lǐng)域中，物理信息神經(jīng)網(wǎng)絡(luò)（physics－informed neural networks，PINN）是研究的熱點，物理信息神經(jīng)網(wǎng)絡(luò)（PINN）是一種科學(xué)的機(jī)器學(xué)習(xí)技術(shù)，可以進(jìn)行偏微分方程的求解。物理信息神經(jīng)網(wǎng)絡(luò)（PINN）是通過進(jìn)行損失函數(shù)優(yōu)化問題進(jìn)行偏微分方程解的逼近，其工作原理是將數(shù)學(xué)模型集成到網(wǎng)絡(luò)中，并用控制方程的殘差項來強(qiáng)化損失函數(shù)，該殘差項作為懲罰項來限制可接受解的空間，其可以被認(rèn)為是一種無監(jiān)督策略，不需要標(biāo)記數(shù)據(jù)[15]。其中RAISSI 等[16]提出用于解決涉及非線性偏微分方程的正反問題的深度學(xué)習(xí)框架以來，物理信息神經(jīng)網(wǎng)絡(luò)（PINN）已經(jīng)成為建模和科學(xué)計算的研究熱點。在進(jìn)行工作面沖擊危險性評價過程中，經(jīng)常遇到相關(guān)指數(shù)難于合理確定，有關(guān)物理機(jī)理難以運用到精準(zhǔn)科學(xué)的評價工作中。基于物理驅(qū)動深度學(xué)習(xí)建模的PDE 求解新方法，具有在不使用任何標(biāo)記數(shù)據(jù)（監(jiān)測數(shù)據(jù)）的情況下快速解決計算力學(xué)問題的優(yōu)勢。為此，以PINN 計算框架為基礎(chǔ)，研究使用較小計算代價，在不失物理機(jī)理的基礎(chǔ)上，分析煤體黏滑沖擊動力學(xué)特性，為快速動態(tài)的沖擊地壓危險性分析方法提供一種技術(shù)手段，以期在生產(chǎn)中為沖擊地壓危險性評價提供參考。

1 沖擊地壓的黏滑物理機(jī)理

1.1 工作面圍巖極限平衡與支承壓力分布

煤層采出后在圍巖應(yīng)力重新分布的范圍內(nèi)作用在煤層、巖層和矸石上的垂直壓力稱為“支承壓力”[17]。準(zhǔn)確地計算預(yù)計采場支承壓力分布規(guī)律及大小對于沖擊地壓動力災(zāi)害預(yù)防具有重要的工程應(yīng)用價值，尤其確定彈性區(qū)長度和支承壓力分布大小是進(jìn)行煤體黏滑沖擊動力學(xué)分析的基礎(chǔ)。

回采工作面前后的支承壓力狀態(tài)，可以劃分為應(yīng)力降低區(qū)、應(yīng)力增高區(qū)、應(yīng)力不變區(qū)，并且把工作面前方支承壓力峰值到煤壁為極限平衡區(qū)，也稱為塑性區(qū)[18]（X0），由峰值向煤體內(nèi)側(cè)為彈性區(qū)（X1）。計算極限平衡區(qū)寬度示意圖如圖1。

圖1 計算極限平衡區(qū)寬度示意圖Fig.1 Schematic diagram of calculating the width of the limit equilibrium zone

通過建立支承壓力分布的平衡微分方程，可進(jìn)行極限平衡區(qū)的計算，極限平衡區(qū)的寬度X0為[18]：

式中：λ 為側(cè)壓系數(shù)；c0、φ 為煤層與頂?shù)装鍘r石交界面的黏聚力與內(nèi)摩擦角；pa為支架對煤幫的支護(hù)阻力；K 為應(yīng)力集中系數(shù)；ρr為上覆巖層的平均密度；g 為重力加速度；H 為開采深度；M 為開采厚度。

支承壓力的影響區(qū)域即支承壓力長度包括塑性區(qū)長度和彈性區(qū)長度，塑性區(qū)寬度按照式（1）進(jìn)行計算，支承壓力的影響區(qū)的長度和支承壓力分布大小按照文獻(xiàn)［19］進(jìn)行確定，這樣確定支承壓力分布函數(shù)形式為線性函數(shù)，但又不失較為嚴(yán)謹(jǐn)?shù)奈锢頇C(jī)理，對構(gòu)建煤體黏滑沖擊動力學(xué)控制方程，進(jìn)行物理信息神經(jīng)網(wǎng)絡(luò)的無監(jiān)督訓(xùn)練有較好的收斂性。支承壓力線性分布示意圖如圖2。

圖2 支承壓力線性分布示意圖Fig.2 Linear distribution of abutment pressure

支承壓力采用簡化的線性分布，支承壓力分布方程為[19]：

式中：Sx為采場推進(jìn)到L0后的支承壓力分布范圍，m；L0為工作面長度，m；cx為采場推進(jìn)距離,m；bx為采場推進(jìn)到L0后的由綜合移動角所形成范圍，m，bx=H·cot θ；θ 為綜合移動角，（°）；p 為彈性區(qū)與塑性區(qū)分布的比值；Kmax為壓力峰值與原始應(yīng)力場比值。

煤體黏滑沖擊動力學(xué)分析的關(guān)鍵是進(jìn)行彈性區(qū)煤體動力特性的表征，進(jìn)行煤體黏滑沖擊動力學(xué)物理驅(qū)動求解的前提需要得到彈性區(qū)支承壓力的分布規(guī)律。通過式（2）和式（3）可得到彈性區(qū)長度X1，即：

1.2 煤體黏滑動力學(xué)控制方程

通過雙剪摩擦實驗和三軸摩擦實驗，表明了砂巖-煤試件最易產(chǎn)生黏滑，且組合煤巖體的摩擦滑動特性表現(xiàn)為典型的黏滑過程，因此可以認(rèn)為黏滑失穩(wěn)破壞是沖擊地壓的一種特殊形式[5]。

煤體黏滑最可能發(fā)生方式為煤體與頂?shù)字g黏滑，尤其是在工作面超前支承壓力的作用下，工作面前方煤體出現(xiàn)應(yīng)力集中，這種頂?shù)装?煤體的組合煤巖體摩擦滑動特性表現(xiàn)為典型的黏滑過程，在一定條件下會導(dǎo)致黏滑失穩(wěn)導(dǎo)致沖擊地壓的發(fā)生。在超前支承壓力的作用下，采煤工作面前方煤體會產(chǎn)生塑性區(qū)X0、彈性區(qū)X1及原巖應(yīng)力區(qū)。在超前支承壓力的作用下，彈性區(qū)積聚了大量的彈性能，對于黏滑失穩(wěn)導(dǎo)致的沖擊地壓現(xiàn)象，進(jìn)行彈性區(qū)黏滑失穩(wěn)的動力學(xué)特性研究有重要的意義。彈性區(qū)黏滑過程的分析示意圖如圖3。

圖3 煤體黏滑沖擊動力學(xué)分析模型Fig.3 Dynamic analysis model of coal stick-slip impact

按圖3 關(guān)系可建立如下煤體動力學(xué)方程[20]：

式中：f 為動摩擦系數(shù)；σx為水平應(yīng)力；σy為垂直應(yīng)力；ρc為煤體密度，kg/m3；u（x，t）為煤體動力學(xué)方程的解；x 為彈性區(qū)水平位置坐標(biāo)；t 為時間，s。

引入胡克定律:

式中：ε 為軸向應(yīng)變；E 為彈性模量。

則式（6）變?yōu)椋?/p>

由本文采用的工作面支承壓力線性分布，可得在彈性分布區(qū)域（X1），垂直應(yīng)力為：

式（9）代入式（8）可得煤體黏滑動力學(xué)偏微分方程：

由式（10）可知，煤體黏滑動力學(xué)系統(tǒng)特性與主要與煤彈性模量、煤密度、煤層厚度、動摩擦系數(shù)、彈性區(qū)域長度、上覆巖層密度、超前應(yīng)力集中系數(shù)、煤層埋深有關(guān)。

2 煤體黏滑沖擊物理驅(qū)動分析原理

2.1 PINN 神經(jīng)網(wǎng)絡(luò)架構(gòu)

基于物理驅(qū)動的深度學(xué)習(xí)（PINN）是近年來開發(fā)的一種用于建模PDE 解決方案的新方法，并顯示出在不使用任何標(biāo)記數(shù)據(jù)（例如測量數(shù)據(jù)不可用）的情況下解決計算力學(xué)問題的前景[21]。

物理信息神經(jīng)網(wǎng)絡(luò)（PINN）結(jié)合了物理知識和神經(jīng)網(wǎng)絡(luò)的方法，可用于解決偏微分方程（組）問題。相比于傳統(tǒng)的PDE 求解手段（有限差分、有限元和譜方法等），PINN 無需離散化空間和時間，使用神經(jīng)網(wǎng)絡(luò)來逼近未知函數(shù)，并使用物理定律作為先驗知識來保證解的合理性，具有高效，無網(wǎng)格計算等優(yōu)勢。RAISSI 考慮通過訓(xùn)練全連接神經(jīng)網(wǎng)絡(luò)來逼近偏微分方程的解，其具體方法如下所述。

假設(shè)非線性偏微分方程：

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

定義f（t，x）由式（11）的左邊給出[15]，通過深度神經(jīng)網(wǎng)絡(luò)進(jìn)行逼近u（t，x），并產(chǎn)生了1 個物理信息的神經(jīng)網(wǎng)絡(luò)f（t，x），該網(wǎng)絡(luò)可以通過應(yīng)用鏈?zhǔn)椒▌t來推導(dǎo)，使用自動微分對函數(shù)的組成進(jìn)行微分，并且具有與表示u（t，x）的網(wǎng)絡(luò)相同的參數(shù)。盡管由于微分算子N 的作用而具有不同的激活函數(shù)，神經(jīng)網(wǎng)絡(luò)u（t，x）和f（t，x）之間的共享參數(shù)可以通過最小化均方誤差損失來學(xué)習(xí)。

通過訓(xùn)練神經(jīng)網(wǎng)絡(luò)進(jìn)行逼近，定義損失函數(shù)MSE 為：

式中：MSEu為邊界和初始條件損失項；MSEf為偏微分方程殘差損失項；Nu為邊界和初始條件的個數(shù)；Nf為內(nèi)部殘差項的個數(shù)；ui為目標(biāo)函數(shù)；（tu，xu）為邊界上的點；（tf，xf）為內(nèi)部的配點。

如果得到網(wǎng)絡(luò)參數(shù)能使控制方程的MSE 趨于0，那么就得到了PDE 的逼近值，即可以認(rèn)為所求結(jié)果趨于真實值。這樣就把求解偏微分方程的問題轉(zhuǎn)化為了如何優(yōu)化損失函數(shù)的問題。

2.2 煤體黏滑沖擊動力學(xué)分析的PINN 模型

通過物理神經(jīng)網(wǎng)絡(luò)（PINN）對煤體黏滑沖擊動力學(xué)控制方程進(jìn)行求解，建立基于PINN 煤體黏滑沖擊動力學(xué)分析模型，來研究煤體黏滑沖擊動力特性。PINN 方法中將偏微分方程邊界條件和初始條件的總和作為神經(jīng)網(wǎng)絡(luò)的損失函數(shù)以達(dá)到求解目的。物理信息神經(jīng)網(wǎng)絡(luò)（PINN）通過使用自動微分將煤體黏滑沖擊動力學(xué)控制偏微分方程（式（10））嵌入到神經(jīng)網(wǎng)絡(luò)的損失函數(shù)中。

煤體動力學(xué)控制方程損失函數(shù)包括為邊界和初始條件損失項和偏微分方程殘差損失項，其中偏微分方程殘差損失項為：

式中：（xi，ti）、（xj，tj）分別是在初始與邊界位置和整個定義域采樣的2 組點。

設(shè)從x=0 到x=X1的一段煤體(彈性區(qū)煤體)發(fā)生黏滑，兩端固定，發(fā)生黏滑前煤體處于靜止?fàn)顟B(tài)，假設(shè)黏滑前煤體位移為0，則邊界條件為：

通過基于梯度的優(yōu)化器的最小化損失來訓(xùn)練網(wǎng)絡(luò)，直到損失小于閾值，由此構(gòu)造1 個神經(jīng)網(wǎng)絡(luò)u（x，t，θ）。其中：θ 為可訓(xùn)練權(quán)重w 和偏差b 的集合；σ 為非線性激活函數(shù)；（xi，ti，ui）為u 指定數(shù)據(jù)；（xj，tj）為PDE 指定殘差點。通過將數(shù)據(jù)和PDE 的加權(quán)損失相加，指定損失L，訓(xùn)練NN 通過最小化損失L 來找到最佳參數(shù)θ*。

3 算例分析

通過給定初始條件的算例，進(jìn)行煤體黏滑沖擊動力學(xué)分析的PINN 模型訓(xùn)練，并驗證其求解結(jié)果的有效精準(zhǔn)性。假設(shè)煤彈性模量為5 GPa，煤密度為1.1 g/cm3，煤層厚度2 m，動摩擦系數(shù)0.15，彈性區(qū)域長度X1為30 m，上覆巖層密度2.5 t/m3，超前應(yīng)力集中系數(shù)為1.5，煤層埋深400 m。邊界條件如式（17）。建立的PINN 求解算例如圖4。

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

在傳統(tǒng)神經(jīng)網(wǎng)絡(luò)中，輸入量常為被預(yù)測量的相關(guān)影響因素，而正向PINN 網(wǎng)絡(luò)是1 種無監(jiān)督學(xué)習(xí)，在研究中，將求解域坐標(biāo)點（x，t）作為模型輸入，Net作為網(wǎng)絡(luò)中間變量輸出，u 作為網(wǎng)絡(luò)最終輸出，從而形成網(wǎng)絡(luò)的映射關(guān)系，具體過程可表述為：

式中：X 為輸入變量，為包含x，t 坐標(biāo)點的張量；σ 為激活函數(shù)，采用ELU 激活函數(shù)；W、b 分別為網(wǎng)絡(luò)權(quán)重和偏置；D（x，t）為距離函數(shù)，滿足當(dāng)x，t 位于定解條件坐標(biāo)時，函數(shù)值為0；BC、IC 分別為邊界條件和初始條件。

由上述可知，一旦滿足相應(yīng)坐標(biāo)，網(wǎng)絡(luò)輸出強(qiáng)制性滿足邊界條件或初始條件，由于正向PINN 是一種無監(jiān)督學(xué)習(xí)求解器，因此需要構(gòu)造目標(biāo)函數(shù)以滿足模型求解需要。根據(jù)偏微分方程形式，可以構(gòu)造1個函數(shù)f（x，t）：

如果f（x，t）趨向于0，則原偏微分方程滿足，根據(jù)f（x，t）以及神經(jīng)網(wǎng)絡(luò)損失函數(shù)性質(zhì)，可將f（x，t）作為PINN 的損失函數(shù)，f（x，t）中的偏導(dǎo)項由自動微分技術(shù)完成。

整個求解網(wǎng)絡(luò)基于Pytorch 1.10 環(huán)境編寫，求解GPU 為Tesla V100 32GB，坐標(biāo)采樣點為30×30，網(wǎng)絡(luò)層輸入節(jié)點為2，隱含層結(jié)構(gòu)為64×4，輸出層節(jié)點為1。為了驗證PINN 求解的正確性，研究基于相同工況采用顯示有限差分法進(jìn)行了對比求解，物理信息神經(jīng)網(wǎng)絡(luò)與顯式有限差分求解結(jié)果對比如圖5，物理信息神經(jīng)網(wǎng)絡(luò)與顯式有限差分局部求解結(jié)果對比曲線（分別為x=5、15、20 m 剖面下2 種模型對比）如圖6。為了反映局部求解性能，依次選擇x=5、15、20 剖面下的求解結(jié)果進(jìn)行對比，模型損失下降曲線如圖7。

圖5 物理信息神經(jīng)網(wǎng)絡(luò)與顯式有限差分求解結(jié)果對比Fig.5 Comparison between physical information neural network and explicit finite difference solution

圖6 物理信息神經(jīng)網(wǎng)絡(luò)與顯式有限差分局部求解結(jié)果對比曲線（分別為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 求解結(jié)果與顯示有限差分求解結(jié)果基本相同，PINN 模型能夠很好地捕捉形變過程，并且2 個模型求解結(jié)果的量級幾乎一致，求解均方根誤差低至1.198 9，說明PINN 可以很好地應(yīng)用于沖擊地壓的動力學(xué)方程求解。

圖6 反映了PINN 在局部求解性能，通過對比曲線可知，2 個模型的求解變化趨勢一致，但在谷值處仍然存在誤差，考慮到有限差分法是一種近似數(shù)值求解，無法通過局部些許差異判斷模型的誤差。所以對PINN 網(wǎng)絡(luò)的損失曲線進(jìn)行分析。

圖7 展示了PINN 模型求解過程的模型損失下降過程，在迭代后期，損失量級達(dá)到了0.001，說明此時模型基本處于收斂狀態(tài)。結(jié)合研究中的PINN損失函數(shù)理論，模型損失可以視為對偏微分方程及定界條件的擬合程度，圖7 中極具收斂的損失可以認(rèn)為模型已經(jīng)對沖擊地壓的動力學(xué)方程完成了求解。

4 深度學(xué)習(xí)物理驅(qū)動的煤體黏滑沖擊動力特性分析

通過算例分析可知，物理神經(jīng)網(wǎng)絡(luò)可以完成對煤體沖擊動力學(xué)控制方程的求解，并且具有較強(qiáng)逼近控制方程真解的能力。由煤體沖擊動力學(xué)控制方程所求水平位移u（x，t），即為研究區(qū)域的巖體水平變形，通過一定的工程研究背景，研究采煤工作面超前支承彈性區(qū)水平變形演化過程，來分析煤體黏滑沖擊動力特性。

4.1 礦井概況

晉華宮礦12 號煤層埋深377～415 m，12 號煤層頂板巖性主要為中粒砂巖，局部為砂質(zhì)泥巖、細(xì)粒砂巖，底板巖性為細(xì)粒砂巖、砂質(zhì)泥巖。12 號煤層沖擊傾向性綜合判定結(jié)果為弱沖擊傾向性，彈性模量5.9 GPa，煤體密度為1.41 g/cm3，抗壓強(qiáng)度13.99 MPa，煤層厚度2.65 m，動摩擦系數(shù)為0.15 應(yīng)力集中系數(shù)為1.5。且所研究工作面推進(jìn)長度大于工作面的長度，X1長度可以按照式（5）進(jìn)行計算，經(jīng)過計算并結(jié)合現(xiàn)場支承壓力實測結(jié)果，取值為50 m。

4.2 煤體黏滑動態(tài)演化分析

以晉華宮礦12 號煤層回采工作面為研究背景，運用物理神經(jīng)網(wǎng)絡(luò)進(jìn)行煤體沖擊動力學(xué)控制方程的求解，共訓(xùn)練50 000 步。K=1.5 位移場預(yù)測如圖8。K=1.5 局部求解結(jié)果對比曲線如圖9。

圖8 K=1.5 位移場預(yù)測Fig.8 Displacement field prediction of K=1.5

圖9 K＝1.5 局部求解結(jié)果對比曲線Fig.9 Comparison curves of local solution results of K＝1.5

通過圖8 可知：煤體彈性區(qū)黏滑具有明顯的波動性，與文獻(xiàn)［5，14］實驗室進(jìn)行煤巖組合結(jié)構(gòu)和巖體結(jié)構(gòu)的黏滑試驗結(jié)果相吻合。并且煤體變形過程出現(xiàn)時間間隔性，變形隨著時間可以恢復(fù)，進(jìn)一步說明了煤體處在彈性變形階段。由于控制方程是基于胡克定律使用位移（變形）的微分替代的水平應(yīng)力微分，所以求解得到的位移（變形）的變化，也能表明水平應(yīng)力也存在波動變化性，并且間接表明了煤體黏滑過程是區(qū)域內(nèi)水平應(yīng)力瞬時的松弛（減?。┖驮黾拥倪^程，煤體黏滑的這種演化過程與沖擊地壓發(fā)生時出現(xiàn)的頂?shù)装逅矔r加卸載及工作面煤層出現(xiàn)瞬時破壞的沖擊地壓孕育過程是一致的。

由圖9 可知：不同時刻煤體黏滑區(qū)域的位移（變形）大小變化較大，但不同時刻變形最大值所處位置基本相同，即煤體黏滑過程，隨著時間的推移水平應(yīng)力的大小在彈性區(qū)內(nèi)存在一定的分布規(guī)律, 不同時刻在距離彈性區(qū)起始位置25 m 左右均出現(xiàn)形變最大值。隨著時間的推移黏滑區(qū)域的位移大小出現(xiàn)波動性。

采場超前支承壓力的變化是工作面煤體沖擊的重要前兆信息，也是沖擊地壓監(jiān)測的主要指標(biāo)之一，所以研究不同支承壓力作用下的煤體黏滑沖擊動力特性，對工作面煤體沖擊災(zāi)害的防治有重要的現(xiàn)實意義。以晉華宮礦12 號煤層回采工作面為研究背景，煤體黏滑沖擊動力控制方程其它參數(shù)保持不變，只改變應(yīng)力集中系數(shù)，分別進(jìn)行應(yīng)力集中系數(shù)K=2、K=3 情況下的基于物理神經(jīng)網(wǎng)絡(luò)的煤體黏滑沖擊動力控制方程求解。

K=2 位移場預(yù)測如圖10，K=2 局部求解結(jié)果對比曲線如圖11。

圖10 K=2 位移場預(yù)測Fig.10 Displacement field prediction of K＝2 at time

圖11 K=2 局部求解結(jié)果對比曲線Fig.11 Comparison curves of local solution results of K＝2

由圖10 可知：當(dāng)采場支承壓力集中系數(shù)為2時，煤體黏滑位移（形變）的最大值有明顯的增加，尤其初期煤體黏滑位移（形變）最大值比集中系數(shù)為1.5 時增加將近1 倍，并且水平應(yīng)力瞬時的松弛（減?。┖驮黾拥倪^程數(shù)值變化幅度較大，初期的2 次位移（形變）的高峰出現(xiàn)明顯的大幅度變化間隔。

由圖11 可知：在不同時刻沿x 軸22 m 位置截面左右均出現(xiàn)形變最大值，在計算時間域內(nèi)位移均出現(xiàn)波動性。

K=3 位移場預(yù)測如圖12，K=3 局部求解結(jié)果對比曲線如圖13。

圖12 K=3 位移場預(yù)測Fig.12 Displacement field prediction of K=3

圖13 K=3 局部求解結(jié)果對比曲線Fig.13 Comparison curves of local solution results of K=3

由圖12 可知：當(dāng)采場支承壓力集中系數(shù)為3時，采場煤體發(fā)生典型的黏滑現(xiàn)象，位移出現(xiàn)準(zhǔn)周期性的變化。初期出現(xiàn)位移的急劇增加，然后瞬時松弛（減小），經(jīng)過短時間的穩(wěn)定，又進(jìn)入急劇增加，然后瞬時松弛（減?。┑闹芷谶^程。每個周期間隔過程都出現(xiàn)位移無法完全恢復(fù)的位移量，且隨著黏滑周期性的波動，無法恢復(fù)的位移量逐漸增大。隨著煤巖系統(tǒng)黏滑多周期性的波動演化，系統(tǒng)儲存彈性能逐漸增大以及煤巖黏滑過程的強(qiáng)度損傷，在一定的采動影響下，有可能會孕育煤層頂?shù)装逅矔r加卸載破壞的沖擊地壓災(zāi)害。

由圖13 可知：煤體沿x 軸各截面形變周期性明顯，且周期間隔的形變量在逐步增加，在不同時刻沿x 軸24 m 位置截面左右均出現(xiàn)形變最大值。

地質(zhì)構(gòu)造異常區(qū)、斷層帶附近易發(fā)生沖擊地壓和其所處煤巖體摩擦系數(shù)的變化有較大關(guān)系[5]，通過改變動摩擦系數(shù)的大小進(jìn)行煤體黏滑沖擊動力控制方程的求解，來研究動摩擦系數(shù)對煤體黏滑沖擊動力特性的影響。f=0.16 位移場預(yù)測如圖14，f=0.20 位移場預(yù)測如圖15。

圖14 f=0.16 位移場預(yù)測Fig.14 Displacement field prediction of f=0.16

圖15 f=0.20 位移場預(yù)測Fig.15 Displacement field prediction of f=0.20

由圖14 可知：當(dāng)動摩擦因數(shù)f 增加0.01 為0.16時，系統(tǒng)位移場數(shù)值變化明顯，所以動摩擦因數(shù)對系統(tǒng)動力學(xué)特性影響較大。

由圖15 可知：當(dāng)動摩擦因數(shù)為2 時，位移場相對最大位移減少了20 mm，比原位移場最大值減少了50%，且通過位移場預(yù)測結(jié)果可知系統(tǒng)黏滑波動性已經(jīng)不明顯，隨著時間推移系統(tǒng)較為穩(wěn)定。

由文獻(xiàn)［5］可知砂巖-煤試樣的摩擦因數(shù)在0.172 8～0.275 3 之間，數(shù)值變化范圍較大，煤巖結(jié)構(gòu)在地質(zhì)構(gòu)造異常區(qū)、斷層帶附近動摩擦因數(shù)急劇減小，極易造成黏滑失穩(wěn)，也為在地質(zhì)構(gòu)造異常區(qū)、斷層帶附近易發(fā)生沖擊地壓找到了摩擦滑動失穩(wěn)的數(shù)值模擬解釋。

5 結(jié) 語

結(jié)合煤體黏滑沖擊理論和物理神經(jīng)網(wǎng)絡(luò)方法，在無任何觀測標(biāo)簽數(shù)據(jù)的情況下建立了基于物理驅(qū)動深度學(xué)習(xí)的煤體黏滑沖擊分析方法，并通過與顯式有限差分求解方法進(jìn)行求解結(jié)果對比，證明了其求解的有效準(zhǔn)確性，并以晉華宮礦12 號煤層工作面為研究背景，進(jìn)行了煤體黏滑沖擊動力學(xué)分析。

1）物理信息神經(jīng)網(wǎng)絡(luò)（PINN）能夠進(jìn)行煤體黏滑沖擊控制方程的數(shù)值模擬求解，并且使用較小的計算代價能夠得到較為準(zhǔn)確有效的數(shù)值模擬結(jié)果，實現(xiàn)了無網(wǎng)格和無標(biāo)簽數(shù)據(jù)情況下的煤體黏滑沖擊數(shù)值模擬。

2）煤體黏滑變形過程出現(xiàn)時間間隔性，水平位移（變形）和水平應(yīng)力均存在波動變化性；煤體黏滑過程是區(qū)域內(nèi)水平應(yīng)力瞬時的松弛（減?。┖驮黾拥倪^程，其與沖擊地壓發(fā)生時出現(xiàn)的頂?shù)装逅矔r加卸載及工作面煤層出現(xiàn)瞬時破壞一致。

3）當(dāng)采場支承壓力集中系數(shù)為2 時，煤體黏滑位移（形變）的最大值有明顯的增加，尤其初期煤體黏滑位移（形變）最大值比集中系數(shù)為1.5 時增加將近1 倍，并且水平應(yīng)力瞬時的松弛（減?。┖驮黾拥倪^程數(shù)值變化幅度較大；當(dāng)采場支承壓力集中系數(shù)為3 時，采場煤體發(fā)生典型的黏滑現(xiàn)象，位移出現(xiàn)準(zhǔn)周期性的變化，初期出現(xiàn)位移的急劇增加，煤巖系統(tǒng)出現(xiàn)黏滑多周期性的波動演化，系統(tǒng)儲存彈性能逐漸增大以及出現(xiàn)煤巖黏滑過程的強(qiáng)度損傷。

4）動摩擦因數(shù)對煤體黏滑動力學(xué)特性影響較大，當(dāng)動摩擦因數(shù)變化時，煤體黏滑系統(tǒng)位移場數(shù)值變化明顯。