王民華，閆永敢

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

引入胡克定律:

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

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

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

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

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

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

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

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

物理信息神經(jīng)網(wǎng)絡(luò)（PINN）結(jié)合了物理知識(shí)和神經(jīng)網(wǎng)絡(luò)的方法，可用于解決偏微分方程（組）問題。相比于傳統(tǒng)的PDE 求解手段（有限差分、有限元和譜方法等），PINN 無需離散化空間和時(shí)間，使用神經(jīng)網(wǎng)絡(luò)來逼近未知函數(shù)，并使用物理定律作為先驗(yàn)知識(shí)來保證解的合理性，具有高效，無網(wǎng)格計(jì)算等優(yōu)勢(shì)。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 個(gè)物理信息的神經(jīng)網(wǎng)絡(luò)f（t，x），該網(wǎng)絡(luò)可以通過應(yīng)用鏈?zhǔn)椒▌t來推導(dǎo)，使用自動(dòng)微分對(duì)函數(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為邊界和初始條件損失項(xiàng)；MSEf為偏微分方程殘差損失項(xiàng)；Nu為邊界和初始條件的個(gè)數(shù)；Nf為內(nèi)部殘差項(xiàng)的個(gè)數(shù)；ui為目標(biāo)函數(shù)；（tu，xu）為邊界上的點(diǎn)；（tf，xf）為內(nèi)部的配點(diǎn)。

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

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

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

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

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

設(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 個(gè)神經(jīng)網(wǎng)絡(luò)u（x，t，θ）。其中：θ 為可訓(xùn)練權(quán)重w 和偏差b 的集合；σ 為非線性激活函數(shù)；（xi，ti，ui）為u 指定數(shù)據(jù)；（xj，tj）為PDE 指定殘差點(diǎn)。通過將數(shù)據(jù)和PDE 的加權(quán)損失相加，指定損失L，訓(xùn)練NN 通過最小化損失L 來找到最佳參數(shù)θ*。

3 算例分析

通過給定初始條件的算例，進(jìn)行煤體黏滑沖擊動(dòng)力學(xué)分析的PINN 模型訓(xùn)練，并驗(yàn)證其求解結(jié)果的有效精準(zhǔn)性。假設(shè)煤彈性模量為5 GPa，煤密度為1.1 g/cm3，煤層厚度2 m，動(dòng)摩擦系數(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ù)測(cè)量的相關(guān)影響因素，而正向PINN 網(wǎng)絡(luò)是1 種無監(jiān)督學(xué)習(xí)，在研究中，將求解域坐標(biāo)點(diǎn)（x，t）作為模型輸入，Net作為網(wǎng)絡(luò)中間變量輸出，u 作為網(wǎng)絡(luò)最終輸出，從而形成網(wǎng)絡(luò)的映射關(guān)系，具體過程可表述為：

式中：X 為輸入變量，為包含x，t 坐標(biāo)點(diǎn)的張量；σ 為激活函數(shù)，采用ELU 激活函數(shù)；W、b 分別為網(wǎng)絡(luò)權(quán)重和偏置；D（x，t）為距離函數(shù)，滿足當(dāng)x，t 位于定解條件坐標(biāo)時(shí)，函數(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個(gè)函數(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)項(xiàng)由自動(dòng)微分技術(shù)完成。

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

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

圖6 物理信息神經(jīng)網(wǎng)絡(luò)與顯式有限差分局部求解結(jié)果對(duì)比曲線（分別為x=5、15、20 m 剖面下2 種模型對(duì)比）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 個(gè)模型求解結(jié)果的量級(jí)幾乎一致，求解均方根誤差低至1.198 9，說明PINN 可以很好地應(yīng)用于沖擊地壓的動(dòng)力學(xué)方程求解。

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

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

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

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

4.1 礦井概況

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 語

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

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

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

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

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