基于HMM的煤礦底板突水短時預測模型

楊 坤，陳伯輝,,施式亮

(1.福州大學環境與資源學院,福建 福州 350108；2.湖南科技大學資源環境與安全工程學院，湖南 湘潭 411201)

煤層底板水害作為常見的煤礦水害類型之一，由于其危險性較大、突發性較強等特性，一直以來是我國煤礦安全研究的重點問題之一。自20世紀60年代至今，我國開展了大量有關煤礦安全方面的研究，并取得了重要成果。近年來，在煤礦底板突水預測方面，常見的方法包括突水系數法、五圖雙系數法、脆弱性指數法、GIS擬合模擬法、突變級數法等[1-8]。 其中，突水系數法簡單易用，但判斷精度不高；五圖雙系數法較為全面，但工作量較大，現場實施較為困難；脆弱性指數法和GIS擬合模擬法較為全面，但存在評測因素選取客觀性不足的問題；突變級數法評測準確度同樣受限于評判指標及其權重的選取。

基于已有研究，本文擬在突水系數法的基礎上，引入隱馬爾科夫模型(HMM)作為煤礦底板突水過程分析的基本原理，并結合Ts-q法[9]，從單位涌水量q、相對安全狀態下最大突水系數Tmax和實際突水系數Ts3個指標構建以條件概率為依托的預測模型，從概率統計的角度出發來規避評測因素選取客觀性不足等問題，使建立的預測模型具有數據獲取難度低、使用簡單、易于推廣等特性。

1 突水系數閾值與單位涌水量

1. 1 突水系數閾值分布

突水系數閾值預測模型在突水系數閾值的設定方面參考了Ts-q法，并在其基礎上做了相關的改進。如圖1所示，常見的Ts-q法通過建立平面直角二維坐標系，確定橫、縱坐標分別為單位涌水量q和突水系數Ts的狀態點分布，當突水系數的實際狀態點處于標定的安全面積內則判定為安全[9]。因此，構成相對安全區域的邊界可視為突水系數閾值即最大突水系數Tmax與單位涌水量q間的函數關系映射。

考慮到Ts-q法中相對安全邊界呈現為折線段形式[10]歸因于突水點分布的統計總結，故在此基礎上探討單位涌水量對應最大突水系數點的概率分布現象。若Ts-q法中相對安全邊界兩側存在突水點隨機分布且可做回歸擬合，則最大突水系數Tmax滿足與單位涌水量q間的線性或非線性關系且表現頻次符合特定概率分布，見圖2。根據突水系數對應的概率分布，進一步構建對應的隱馬爾科夫鏈。

圖2 突水系數概率分布示意圖Fig.2 Schematic diagram of probability distribution of water inrush coefficient

1.2 最大突水系數與單位涌水量關系模型的構建

1.2.1 模型機制

在假定的隱馬爾科夫模型中，單位涌水量作為觀測量O，相對安全狀態下的最大突水系數作為隱藏狀態量X，構成隱馬爾科夫鏈[11-14]，見圖3。其中,X1,X2,…,XT分別指代不同觀測時間段相對安全狀態下的最大突水系數Tmax；O1,O2,…,OT分別指代不同觀測時段的單位涌水量q。

圖3 最大突水系數與單位涌水量的隱馬爾科夫模型機制Fig.3 Mechanism of hidden Markov model between water inrush coefficient and unit water inflow

1.2.2 模型要素構成

突水系數閾值即最大突水系數Tmax預測模型要素主要由單位涌水量觀測區間個數n、狀態轉移矩陣A、表達矩陣B和初始狀態概率向量π構成。

(1) 單位涌水量觀測區間個數n。通過劃分單位涌水量觀測區間，確定觀測集的元素個數及對應隱藏狀態集的元素個數n，若初始模型的單位涌水量觀測區間分為n段，則O1,O2,…,On為觀測事件集合，即單位涌水量觀測值的不同分布區間，X1,X2,…,Xn為隱藏狀態集合，即不同單位涌水量分布區間對應的相對安全狀態下的最大突水系數。

(2) 狀態轉移矩陣A。狀態轉移矩陣A為相對安全狀態下最大突水系數Tmax的變化概率矩陣[11-15]，若模型預先設定的觀測分段數為n，則該矩陣為n×n矩陣：

其中，aij表示第i個涌水量觀測區間對應的相對安全狀態下最大突水系數Tmaxi在下次數據采集時變化為第j個涌水量觀測區間對應最大突水系數Tmaxj的概率，考慮到實際作業過程中隱藏狀態變化的隨機性和涌水量變化的不確定性，可設定對于任意i≤n且j≤n,aij取值為1/n。

(3) 表達矩陣B。表達矩陣B為任意單位涌水量區間qi(1≤i≤n)對應相對安全狀態下最大突水系數Tmaxi的實際單位涌水量觀測值表現為任意單位涌水量區間qj(1≤j≤n)的概率矩陣[11-15]。首先，根據所在礦區歷史突水點的對應單位涌水量q、對應監測點突水系數Ts繪制Tmax為自變量的q值的分布圖,并基于Ts-q中相對安全邊界概念選取可作為臨界參考值的突水點作為參與計算點，引入q值離散點概率分布概念，對q值離散點做回歸擬合，見圖4；然后，根據區間數n劃分單位涌水量觀測區間，基于擬合曲線確定各觀測區間qi對應的突水系數最大值作為Tmaxi(1≤i≤n)，由于礦區突水點分布適用于簡單非線性回歸，故假設在自然條件下單位涌水量的區間分布概率符合以相對安全狀態下最大突水系數Tmax為自變量的高斯分布，見圖5;最后，通過概率分布計算，可得到表達矩陣B：

(1)

其中，Pij為第i個觀測區間對應Tmaxi表現為第j個觀測區間單位涌水量的概率;gi(x)為相對安全狀態下最大突水系數為Tmaxi時單位涌水量q的概率分布函數;max(qj)和min(qj)分別為第j個觀測區間的單位涌水量最大值和最小值，μi為第i個觀測區間對應的單位涌水量均值；σ2表示離散程度，由單位涌水量q與相對安全狀態下最大突水系數Tmax的關聯性強弱決定，可根據斷層構造及性質等條件進行判別，在實際條件不明的狀態下，可根據參與計算突水點最大突水系數-單位涌水量回歸曲線方差確定。

圖4 礦區突水點單位涌水量分布示例圖Fig.4 Example diagram of unit water inflow distribution of water inrush point in mining area

圖5 突水系數對應單位涌水量的概率分布示意圖Fig.5 Probability distribution diagram of unit water inflow corresponding to water inrush coefficient

(4) 初始狀態概率向量π。初始狀態概率向量π=(π1,π2,…,πn)[11-15]，其中πi為初始時刻相對安全狀態下最大突水系數Tmax為單位涌水量觀測區間qi對應最大突水系數Tmaxi的概率，可通過統計礦區參與計算突水點分布，構造單位涌水量q為自變量的回歸方程，計算當前單位涌水量對應最大突水系數的概率分布，進而求得：

(2)

其中，h(x)為q值確定情況下Tmax的概率分布函數;Tmaxi為觀測區間qi對應Tmax回歸曲線的Tmax最小值,當i=n時，Tmax(n+1)取區間qn中觀測點突水系數最大值;k為修正系數,用于放大突水系數計算區間;μ為當前單位涌水量觀測值所屬區間對應回歸曲線段的突水系數均值;σ2可根據參與計算突水點由單位涌水量-最大突水系數回歸曲線方差確定。

1. 3 最大突水系數的預測

通過對突水系數閾值預測模型使用Viterbi算法，根據某時刻至今單位涌水量數據，推算可能性最大的對應相對安全狀態下最大突水系數Tmax的變化過程[16-18]，從而計算得出當前Tmax的預測值。

過程如下：

(1) 輸入變量：模型為λ=(A,B,π)，其中A為狀態轉移矩陣,B為表達矩陣,π為初始狀態概率向量。單位涌水量觀測序列(O1,O2,…,OT),其中O1,O2,…,OT分別指代不同觀測時段的單位涌水量所屬區間。

(2) 輸出變量：相對安全狀態下最大突水系數的變化序列(X1,X2,…,XT),其中X1,X2,…,XT指代與不同觀測時間段相對應的最大突水系數預測值。

(3) 初始化:算法初始值計算公式為

δt(i)=πibi(O1) (i=1,2,…,n)

(3)

φt(i)=0 (i=1,2,…,n)

(4)

其中，δt(i)為在監測時間節點t時刻單位涌水量屬于qi區間的最大突水系數Tmax最優路徑經過結點的聯合概率最大值;πi為初始時刻相對安全狀態下最大突水系數Tmax為Tmaxi的概率；bi(O1)為Tmaxi表現為觀測事件1即初始觀測事件的概率；φt(i)為t時刻Tmax的最大概率變化。

(4) 遞推計算：遞推計算t=2,3,…,T時最優路徑經過結點的聯合概率以及計算該時刻對應隱藏狀態:

(5)

(6)

其中，i=1,2,…,n,j=1,2,…,n為狀態序號;aji為狀態轉移矩陣A中對應元素，表示由上一狀態j轉化為狀態i的概率。

(5) 終止:通過T次計算，得出最終最大路徑概率P*及Tmax最終變化值XT為

(7)

(8)

(6) 最優路徑回溯:對于t=T-1,T-2,…,1有:

Xt=φt+1(Xt+1)

(9)

求得最優路徑為

X*=(X1,X2,…,XT)

(10)

其中，XT即為當前相對安全狀態下最大突水系數的預測值。

2 突水系數與單位涌水量

2.1 突水系數與單位涌水量的關系模型構建

2.1.1 模型機制

對于實際突水系數Ts的預測，同樣采用隱馬爾科夫鏈進行模型構建，實際突水系數預測模型中單位涌水量作為觀測量，而實際突水系數則作為隱藏狀態量，構成隱馬爾科夫鏈。

2.1.2 模型要素構成

實際突水系數預測模型要素主要由單位涌水量觀測區間個數m、狀態轉移矩陣A、表達矩陣B和初始狀態概率向量π構成。但與相對安全狀態下最大突水系數預測模型相比，實際突水系數預測模型的初始參數及使用方法存在差異。由于實際突水系數與單位涌水量之間存在一定關聯但具有不確定性，故實際突水系數預測模型需要適量初始測算數據作為建模依據，包括一定時段各單位涌水量區間qi對應的實際突水系數監測數據即底板水壓和隔水層厚度。

另外，在數據選用方面，該模型無須劃分有效點作為參與計算點；在區間劃分上，考慮到該模型監測數據量較突水點數據量大，可在適當范圍內增大觀測區間數m以提高預測精度；在參數構造方法選擇方面，考慮到該模型實際突水系數測算數據獲取有限，難以形成大量訓練樣本來支撐常見監督學習算法，故優先采用概率分布算法，若實際突水系數監測數據獲取難度較大或數據可靠性低，則采用Baum-Welch算法[18]基于單位涌水量觀測值開展實際突水系數的計算。

(1) 狀態轉移矩陣A:對于狀態轉移矩陣A，根據模型預先設定的觀測區間分段個數m來確定該矩陣結構：

其中，aij表示第i個涌水量觀測區間對應的實際突水系數Tsi在下次數據采集時變化為第j個涌水量觀測區間對應的實際突水系數Tsj的概率,考慮到實際作業中含水層中的靜水壓力與隔水層厚度變化的隨機性，可初步設定對于任意i≤m且j≤m,aij取值為1/m。若缺少實際突水系數Ts觀測數據，則采用Baum-Welch算法[18]基于單位涌水量觀測值開展模型參數的計算。

(2) 表達矩陣B:對于表達矩陣B，其矩陣結構為

若實際突水系數Ts觀測數據充足且可靠，則根據觀測點Ts值計算其回歸曲線，劃分單位涌水量區間qi及對應突水系數區間Tsi(1≤i≤m)，則有：

(11)

其中，Pij為第i觀測區間對應Tsi表現為第j觀測區間單位涌水量的概率；fi(x)為實際突水系數Ts對應Tsi時單位涌水量q的概率分布函數;max(qj)和min(qj)分別為第j觀測區間的單位涌水量最大值和最小值;μi為第i觀測區間對應單位涌水量均值;σ2表示離散程度，可根據觀測點實際突水系數-單位涌水量回歸曲線方差確定。

若缺少Ts觀測數據或Ts對q值分布離散程度較大，可在初始樣本數據積累之上，采用Baum-Welch算法[18]進行表達矩陣計算。故初步設定對于任意i≤m且j≤m,bij取值為1/m。

(3) 初始狀態概率向量π:對于初始狀態概率向量π=(π1,π2,…,πm)，需在初始樣本數據積累之上，構造以q為自變量、Ts為因變量的回歸曲線，計算q值確定情況下Ts的概率分布，進而求得：

(12)

其中,l(x)為q值確定情況下Ts的概率分布函數;max(Tsi)和min(Tsi)分別為突水系數區間Tsi的最大值和最小值，當i=m時，Ts(m+1)為觀測區間qm中觀測點Ts的最大值；k為修正系數，用于放大突水系數計算區間；μ為當前q值所屬區間對應Ts的平均值；σ2可根據觀測點單位涌水量-實際突水系數回歸曲線方差確定。

若觀測工作較為復雜，無法準確測算Ts值，可初步設定πi取值為1/m，采用Baum-Welch算法[18]對模型初始參數進行修正。

2.1.3 模型初始參數修正

基于Baum-Welch算法[18]的模型參數修正計算過程如下：

(1) 輸入變量：單位涌水量觀測序列V(O1,O2,…,OT)。

(2) 輸出變量：隱馬爾可夫模型參數λ=(A,B,π)。

λ(0)=(A(0),B(0),π(0))

(13)

(4) 遞推計算： 遞推，對于s=1,2,3,…,則有：

(14)

(15)

(16)

其中,T為觀測序列長度;ξt(i,j)為給定模型λ和觀測序列V情況下，在時刻t表現為qi且在時刻t+1表現為qj的概率;γt(i)為給定模型λ和觀測序列O情況下，涌水量區間qi與突水系數區間Tsi相匹配的概率;ξt(i,j)和γt(i)由λ(s-1)=(A(s-1),B(s-1),π(s-1))和觀測序列V求出;Ot=qi表示t時刻觀測值Ot屬于涌水量區間qi。

(5) 終止：經過多次迭代，λ(s)=(A(s),B(s),π(s)),即λ的極大似然估計。

2. 2 實際突水系數的預測

通過多次測算不同單位涌水量區間對應的實際突水系數，形成預測模型的初始樣本數據，進而計算得到狀態轉移矩陣A、表達矩陣B和初始狀態概率向量π，采用Viterbi算法[16-18]，根據某時刻至今單位涌水量數據，推算實際突水系數的變化過程，從而計算得出當前實際突水系數的預測值。

3 基于單位涌水量的煤礦底板突水短時預測實例應用與分析

將上述建立的突水系數閾值預測模型和實際突水系數預測模型進行功能組合，則構成基于單位涌水量的煤礦底板突水短時預測模型。通過對比該復合預測模型計算得出的相對安全狀態下最大突水系數預測值與實際突水系數預測值，若當前實際突水系數預測值大于相對安全狀態下最大突水系數預測值，則可初步判定當前狀態下發生煤礦底板突水的可能性較大。實際應用中，該預測模型適用于由于部分客觀因素導致難以開展突水系數定時測算的工作情況。

本文以焦作礦區為案例，通過對韓王煤礦、焦西煤礦、王封煤礦、李封煤礦、演馬莊煤礦、九里山煤礦和中馬村煤礦7個煤礦多個突水點采用基于單位涌水量的煤礦底板突水短時預測模型進行案例計算與分析[19]，以驗證該預測模型的可行性。

根據焦作礦區突水點分布[9,19]，選取用于構建相對安全邊界的參與計算突水點，繪制突水點分布圖，并進行回歸曲線擬合，見圖6。

圖6 焦作礦區突水點最大突水系數安全邊界擬合曲線Fig.6 Fitting safety boundary of maximum water inrush coefficient of water inrush point in Jiaozuo Mining Area

然后，構造單位涌水量觀測區間，代入擬合曲線計算對應單位涌水量區間的最大突水系數，其結果見表1。

表1 觀測區間及對應的最大突水系數

根據擬合曲線繪制各單位涌水量觀測區間所含突水點突水系數均值對應的單位涌水量概率分布,計算得到表達矩陣B如下

設狀態轉移矩陣A為10×10矩陣，且矩陣A所有元素取值為0.1。

若取連續單位涌水量觀測數據向量為：[1.4,1.4,1.4,1.4,1.4,1.4,1.4,1.4,1.4,1.4]，修正系數k設為50，則計算可得初始狀態概率向量π=[0.03,0.12,0.14,0.12,0.09,0.07,0.05,0.19,0.06,0.12]。

使用Viterbi算法[16-18]計算得到相對安全狀態下最大突水系數的最優路徑為：[0.032,0.075,0.075,0.075,0.075,0.075,0.075,0.075,0.075,0.075]，即最大突水系數預測值為0.075 MPa/m。

若取連續單位涌水量觀測數據向量為：[1.4,1.0,0.6,1.2,1.0,1.4,1.2,1.6,1.4,1.6]，則同理使用Viterbi算法[16-18]可計算得到相對安全狀態下最大突水系數的最優路徑為：[0.032,0.112,0.112,0.075,0.112,0.075,0.075,0.048,0.075,0.048]，即最大突水系數預測值為0.048 MPa/m。

本文基于單位涌水量的煤礦底板短時預測方法在相對安全狀態下突水系數的預測最終結果與Ts-q法[9]整體較為相近，說明本文方法具有可行性。部分差異包括：①當前單位涌水量的前期連續觀測值會對最終預測結果產生影響；②預測值隨單位涌水量持續變動，在部分區間相較Ts-q法[9]偏高或偏低。

4 結論與建議

基于隱馬爾科夫模型，結合條件概率理論，探索在監測數據有限條件下通過單位涌水量開展煤礦底板突水短時預測的新方法，得到主要結論如下：

(1) 基于Ts-q法，實現了突水系數閾值的概率分布設定，可用于構建突水系數預測模型。

(2) 基于隱馬爾科夫鏈和概率統計算法，構造以單位涌水量為觀測數據的最大突水系數和實際突水系數預測模型，實現了煤礦底板突水的短時預測。

(3) 在條件允許的情況下，實際突水系數采用現場測算的方法可提升模型的預測精度。此外，結合相關實際數據進行反饋分析，本模型在相關參數的標定算法上仍可做進一步的改善。