回采工作面瓦斯涌出量預(yù)測模型

邵良杉，王 振

（遼寧工程技術(shù)大學(xué) 系統(tǒng)工程研究所，遼寧 葫蘆島 125105）

0 引言

瓦斯是煤礦安全的主要隱患之一，中國瓦斯事故占到煤礦事故的80%以上.準(zhǔn)確預(yù)測瓦斯涌出，不僅是瓦斯防護(hù)和治理、煤礦通風(fēng)系統(tǒng)設(shè)計(jì)的重要基礎(chǔ)[1]，也為礦井高質(zhì)量生產(chǎn)提供了安全保障.隨著工業(yè)化、信息化的發(fā)展，機(jī)械程度不斷提高，礦井開采深度不斷增加，工作面的產(chǎn)能加重，隨之帶來工作面瓦斯涌出量的上升，制約著煤礦的安全生產(chǎn)和經(jīng)濟(jì)效益[2].因此研究一種瓦斯涌出量的預(yù)測模型顯得極為重要.

近年來，許多學(xué)者針對瓦斯涌出量預(yù)測進(jìn)行了研究.呂伏等[3]提出將主成分分析法應(yīng)用于采掘工作面瓦斯涌出的多步線性回歸預(yù)測中.萬仁寶等[4]將灰色系統(tǒng)與神經(jīng)網(wǎng)絡(luò)相結(jié)合，在小信息量的情況下做出準(zhǔn)確預(yù)測.付華等[5]通過優(yōu)化粒子群和Elman神經(jīng)網(wǎng)絡(luò)的權(quán)重和閾值建立絕對瓦斯涌出量預(yù)測系統(tǒng)模型.施式亮等[6]通過分解瓦斯實(shí)時(shí)涌出數(shù)據(jù)，通過粒子群算法優(yōu)化支持向量機(jī)（SVM），利用優(yōu)化過的SVM預(yù)測多個(gè)本征模函數(shù)，并對結(jié)果進(jìn)行累加從而得出瓦斯預(yù)測值.溫廷新等[7]運(yùn)用灰色關(guān)聯(lián)和因子分析對指標(biāo)進(jìn)行分析提取，利用量子遺傳算法對最小二乘支持向量機(jī)進(jìn)行了優(yōu)化，建立了煤與瓦斯突出預(yù)測模型.林海飛等[8]為提高煤層瓦斯含量預(yù)測的科學(xué)性及準(zhǔn)確性，提出基于粒子群優(yōu)化誤差反向傳播神經(jīng)網(wǎng)絡(luò)的瓦斯含量預(yù)測模型.

以上研究存在一定局限性，瓦斯涌出情況多變且不穩(wěn)定，傳統(tǒng)預(yù)測方法的預(yù)測結(jié)果與實(shí)際情況相差較遠(yuǎn)，因此需要更為優(yōu)化的預(yù)測方法.本文采用因子分析去除數(shù)據(jù)相關(guān)性并濃縮數(shù)據(jù)的優(yōu)點(diǎn)，運(yùn)用改進(jìn)遺傳算法，對BP網(wǎng)絡(luò)進(jìn)行優(yōu)化，建立預(yù)測回采工作面瓦斯涌出量的優(yōu)化模型，并通過對實(shí)測數(shù)據(jù)進(jìn)行預(yù)測和誤差分析驗(yàn)證該模型的準(zhǔn)確性.

1 因子分析及IGA-BP理論

1.1 因子分析原理

因子分析（factor analysis，F(xiàn)A）通過研究變量之間的內(nèi)部依賴性，對具有錯(cuò)綜復(fù)雜關(guān)系的變量進(jìn)行簡化.因子分析依據(jù)其運(yùn)算出發(fā)點(diǎn)的不同可以分為兩類，即R型與Q型，前者以相關(guān)系數(shù)矩陣作為基礎(chǔ)，而后者則是通過相似系數(shù)矩陣實(shí)現(xiàn)[9]，本文采用R型因子分析.

設(shè)p個(gè)變量Xi(i= 1,2,… ,p)，待觀測變量通過m(m＜k)個(gè)兩兩正交的公共因子F= [f1,f2,… ,fm]T表示，其數(shù)學(xué)表達(dá)式為

式中，每個(gè)變量都是標(biāo)準(zhǔn)化的變量，iε為特殊因子，滿足E(εi) = 0，作用于其對應(yīng)的Xi(i=1,2,… ,p)變量；αij為Xi在fj上的負(fù)擔(dān)承載，與fj和Xi的關(guān)系的強(qiáng)弱成正比.模型的矩陣形式表示為

式中，Λ為公共因子F的負(fù)載矩陣；F為通過因子載荷矩陣找出，其運(yùn)算過程可以轉(zhuǎn)化為

式中，R(X)*= cov(X)-cov(ε)，若ε未知，可用R(X)代替R(X)*.可證明，就是一個(gè)解，其中α1=α11+α21+ … +αk1，λ1是ΛΛT的一個(gè)特征根，1b是對應(yīng)的征根1λ的任意一個(gè)模為1的特征向量，記；則可證也是所求.同理可求α3,… ,αm.進(jìn)一步解因子值，可得R為X的相關(guān)系數(shù)矩陣，令βj=R-1αj，通過Fi=Xβj計(jì)算，得到公共因子值.

1.2 改進(jìn)遺傳算法

近年來，遺傳算法[10]（genetic algorithm，GA）應(yīng)用廣泛，但在應(yīng)對局部搜索時(shí)，有大概率陷入“早熟”的問題，影響其收斂性.針對該問題，對GA算法改進(jìn)（IGA算法），方法如下.

（1）選擇每個(gè)交叉和突變產(chǎn)生的個(gè)體，并將交叉概率Pc和突變概率Pm設(shè)置為固定值，不隨環(huán)境變化而變化.

式中，f ′為交叉染色體中較大的適應(yīng)值；為種群中平均適應(yīng)值；fmax為最大適應(yīng)函數(shù)值；δ為最優(yōu)的個(gè)體的適應(yīng)函數(shù)值同種群的平均適應(yīng)函數(shù)值之間的差值；f為個(gè)體的適應(yīng)度函數(shù)值.

在交叉和變異動(dòng)態(tài)調(diào)整的過程中，為了避免算法較早地收斂到局部最優(yōu)值和個(gè)體的過度繁殖，染色體的適應(yīng)度函數(shù)值的極大值要無限趨近或等于最大適應(yīng)度函數(shù)值的中間值.

（2）在遺傳算法初期確定適應(yīng)度函數(shù)時(shí)，通過探究其優(yōu)化懲罰系數(shù)γ，懲罰函數(shù)設(shè)計(jì)達(dá)到最佳效果.其極值問題為

式中，f(x)為適應(yīng)度函數(shù)；x為待求解函數(shù)g(x)的一個(gè)可行解；γ為懲罰系數(shù)；hi(x)為松弛變量，i=1,2,… ,m.

對懲罰系數(shù)γ進(jìn)一步優(yōu)化為

式中，t為迭代次數(shù)，當(dāng)t較小時(shí)，γ接近于γmax，保證了算法的全局搜索能力；隨著t的增大，γ以非線性減小，從而保障了算法的局部求解能力.

2 模型設(shè)計(jì)及數(shù)據(jù)分析

2.1 基于IGA的BP神經(jīng)網(wǎng)絡(luò)流程

IGA-BP算法步驟如下.

步驟1BP神經(jīng)網(wǎng)絡(luò)賦初值，設(shè)定網(wǎng)絡(luò)各層神經(jīng)元個(gè)數(shù)，設(shè)置網(wǎng)絡(luò)參數(shù)，完成網(wǎng)絡(luò)連接權(quán)值和閾值的初始化.

步驟2對GA進(jìn)行優(yōu)化，初始化IGA的種群，設(shè)置IGA參數(shù).

步驟3提供相應(yīng)的BP神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)模式.

步驟4計(jì)算適應(yīng)度函數(shù).選擇均方差來作適應(yīng) 度函數(shù)，為

式中，M為訓(xùn)練樣本總數(shù)；kjy為第k個(gè)樣本第j個(gè)網(wǎng)絡(luò)輸出節(jié)點(diǎn)的理論輸出；為第k個(gè)樣本第j個(gè)網(wǎng)絡(luò)輸出節(jié)點(diǎn)的實(shí)際輸出；C為網(wǎng)絡(luò)樣本的數(shù)目.

步驟5染色體解碼.

步驟6群體評(píng)價(jià).

步驟7染色體選擇.將評(píng)價(jià)結(jié)果進(jìn)行比較，選擇性能優(yōu)秀的染色體進(jìn)行下一步操作.

步驟8染色體交叉.

步驟9染色體保留.

步驟10在BP神經(jīng)網(wǎng)絡(luò)中，利用IGA對其中的權(quán)值、閾值進(jìn)行優(yōu)化，同時(shí)導(dǎo)入到BP算法中以初始值進(jìn)行訓(xùn)練，并根據(jù)訓(xùn)練結(jié)果進(jìn)行權(quán)值和閾值的調(diào)整.如果達(dá)到預(yù)定的誤差要求，則停止迭代，輸出最終值，否則進(jìn)入步驟3.

2.2 參數(shù)選取及相關(guān)性分析

經(jīng)過實(shí)際現(xiàn)場調(diào)研及相關(guān)資料的參考，選擇12個(gè)因素作為回采工作面瓦斯涌出量的主要影響因素[11]：原始瓦斯量X1、煤層埋藏深度X2、煤層厚度X3、傾角X4、工作面長度X5、推進(jìn)速度X6、采出率X7、臨近層瓦斯量X8、臨近層厚度X9、臨近層間距X10、巖性X11、開采深度X12.

參照文獻(xiàn)[12]提供的某礦18組回采工作面瓦斯涌出量數(shù)據(jù)，對以上提及的12個(gè)因素進(jìn)行Pearson相關(guān)性的分析.Pearson系數(shù)及Sig.(2-tailed)檢驗(yàn)結(jié)果表明：煤層傾角、工作面長度、臨近厚度、臨近層間距與其他8個(gè)指標(biāo)之間的t值的顯著性概率p大于0.05，相關(guān)系數(shù)不異于0，相關(guān)性較弱.原始瓦斯量、煤層埋藏深度、煤厚、推進(jìn)速度、采出率、臨近層瓦斯量、巖性、開采深度這8個(gè)指標(biāo)t值的顯著性概率p小于0.05，可見相關(guān)系數(shù)顯著異于0.故而變量的關(guān)聯(lián)性較強(qiáng)，可以降維.

在12個(gè)主要因素中提取8個(gè)高度相關(guān)的指標(biāo)并對其進(jìn)行因子分析，從而得到少數(shù)公共因子，實(shí)現(xiàn)對回采工作面瓦斯涌出量的影響因素降維.為保證指標(biāo)選取的可靠性，通過KMO值和巴特利球體檢驗(yàn)對8個(gè)高度相關(guān)的指標(biāo)檢測：樣本KMO值為0.8大于0.6，屬于允許域；巴特利球度檢驗(yàn)近似卡方為164.193，顯著性概率p為0.000，小于0.05；可知在各個(gè)指標(biāo)數(shù)據(jù)間具有一定的相關(guān)性并能夠進(jìn)行因子分析.

3 瓦斯涌出量預(yù)測模型

3.1 實(shí)驗(yàn)數(shù)據(jù)選取

將實(shí)測數(shù)據(jù)中原始瓦斯量、煤層埋藏深度、厚度、推速、采出率、臨近層瓦斯量、巖性、開采深度8個(gè)指標(biāo)作為因子分析的變量，采用SPSS軟件對回采工作面瓦斯涌出量的影響因素進(jìn)行回歸分析，在分析方法上，采用主成分提取法，得出公共因子碎石圖，見圖1.

圖1 公共因子碎石圖Fig.1 common factor scree plot

由圖1可以看出，前2個(gè)成分的特征值有較為明顯坡度變化，而自第3個(gè)成分的特征值開始，其變化幅度快速地趨于穩(wěn)定.由碎石準(zhǔn)則[13]可見，所取的兩個(gè)公共因子能夠更加明顯地實(shí)現(xiàn)對原變量方差的解釋.

總方差解釋見表1.由表1可知，第一和第二個(gè)公共因子的方差累積貢獻(xiàn)率達(dá)到91.137%，包含了總數(shù)據(jù)的大部分信息.提取2個(gè)公共因子f1和f2.參數(shù)空間的維數(shù)由原來的8維降至2維后，消除了冗余，算法的泛化能力得以增強(qiáng).

表1 總方差解釋Tab.1 total variance explained

所選取的初始因子具有過強(qiáng)的綜合性，故需對所取因子做一定處理從而降低初始因子的綜合性.本文選用正交旋轉(zhuǎn)的方法，由旋轉(zhuǎn)成分矩陣（表2）可得，公共因子f1由原始瓦斯量、煤層埋藏深度、厚度、推速、采出率、開采深度構(gòu)成，公共因子f2由臨近層瓦斯量、巖性構(gòu)成.以回歸法實(shí)現(xiàn)因子值計(jì)算過程獲得因子得分矩陣表見表3.

表2 旋轉(zhuǎn)成分矩陣Tab.2 rotated component matrix table

表3 成分得分系數(shù)矩陣Tab.3 componentscore coefficient matrix table

3.2 實(shí)驗(yàn)數(shù)據(jù)處理

由于煤層傾角X4、工作面長度X5、臨近厚度X9及臨近層間距X10與另外8個(gè)變量之間的相關(guān)性較弱，因此將煤層傾角、工作面長度、臨近厚度、臨近層間距，以及所選取的2個(gè)公共因子作為IGA-BP神經(jīng)網(wǎng)絡(luò)的特征因子.為便于數(shù)據(jù)處理，加快訓(xùn)練網(wǎng)絡(luò)的收斂速度，對數(shù)據(jù)進(jìn)行規(guī)范化[14]，本文采用Z-score歸一化方法對其他指標(biāo)數(shù)據(jù)進(jìn)行處理，部分?jǐn)?shù)據(jù)見表4.

表4 回采工作面瓦斯涌出量相關(guān)數(shù)據(jù)Tab.4 data about gas emission quantity of working face

3.3 FA-IGA-BP神經(jīng)網(wǎng)絡(luò)預(yù)測

（1）參數(shù)設(shè)置

瓦斯涌出量預(yù)測結(jié)果準(zhǔn)確性與BP神經(jīng)網(wǎng)絡(luò)隱含層節(jié)點(diǎn)數(shù)目息息相關(guān)[15-16]，節(jié)點(diǎn)數(shù)過少將引起網(wǎng)絡(luò)訓(xùn)練次數(shù)增多，訓(xùn)練結(jié)果的精度也會(huì)隨之下降；網(wǎng)絡(luò)節(jié)點(diǎn)數(shù)目過多會(huì)延長訓(xùn)練時(shí)間，使網(wǎng)絡(luò)過擬合化.為解決這類問題，需要設(shè)計(jì)者進(jìn)行反復(fù)實(shí)驗(yàn)和依據(jù)自身的經(jīng)驗(yàn)來確定[17].本文借鑒“試錯(cuò)法”的思想，由式（9）、式（10）和式（11）求BP網(wǎng)中隱含層神經(jīng)元的個(gè)數(shù)，試驗(yàn)系統(tǒng)的訓(xùn)練結(jié)果，進(jìn)而確定較優(yōu)的BP網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu).

式中，m為隱含層的節(jié)點(diǎn)數(shù)；n為輸入層的節(jié)點(diǎn)數(shù)；l為輸出層的節(jié)點(diǎn)數(shù)；α為1～10之間的調(diào)節(jié)參數(shù).結(jié)合實(shí)驗(yàn)數(shù)據(jù)情況，選取BP神經(jīng)網(wǎng)絡(luò)隱含層節(jié)點(diǎn)數(shù)區(qū)間為[3，13]，測試隱含層節(jié)點(diǎn)從3到13的仿真誤差，對其擬合后的數(shù)據(jù)特征進(jìn)行比較，結(jié)果見表5.

表5 隱含層節(jié)點(diǎn)數(shù)目測試結(jié)果Tab.5 test results of hidden layer nodes

由表5可知，隱含層節(jié)點(diǎn)數(shù)為13時(shí)，3次重復(fù)實(shí)驗(yàn)均具有明顯優(yōu)越性，故選定隱含層節(jié)點(diǎn)數(shù)為13.

選擇1～15組數(shù)據(jù)作為網(wǎng)絡(luò)的訓(xùn)練集，16～18組數(shù)據(jù)作為網(wǎng)絡(luò)的測試集，并在Matlab R2016b中建立基于FA-IGA-BP算法的回采工作面瓦斯涌出量預(yù)測模型.BP神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)為“6-13-1”，即網(wǎng)絡(luò)包含6個(gè)輸入節(jié)點(diǎn)、13個(gè)隱含層節(jié)點(diǎn)和1個(gè)輸出層節(jié)點(diǎn)，傳遞函數(shù)為{'logsig'，'purelin'}.訓(xùn)練算法采用標(biāo)準(zhǔn)梯度算法trainlm，訓(xùn)練目標(biāo)為0.00001，學(xué)習(xí)速率為0.3，訓(xùn)練次數(shù)為6000.IGA每個(gè)種群的維數(shù)為105×(6×13+13×13+1)，其余參數(shù)設(shè)置：群體大小為20；最大允許迭代次數(shù)為1000；最大權(quán)值ωmax為0.9；最小權(quán)值ωmin為0.4；迭代次數(shù)為 1000；交叉概率Pc1為0.4、Pc2為0.9；變異概率Pm1為0.01、Pm2為0.1；δ為0.001.

訓(xùn)練集經(jīng)過訓(xùn)練后，IGA算法的適應(yīng)度值可以達(dá)到3.38×10-6，收斂性能和精度都較高.從進(jìn)化的過程中可以看出，IGA對BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練是良好的，可以得到一系列適用于回采工作面瓦斯涌出量預(yù)測的最優(yōu)權(quán)值和閾值，數(shù)據(jù)見表6.

表6 IGA得到的最優(yōu)權(quán)值和閾值Tab.6 optimal weights and thresholds obtained by IGA

（2）預(yù)測結(jié)果分析

采用回代估計(jì)方法得到訓(xùn)練集瓦斯涌出量的預(yù)測值，并與實(shí)測值比較[18-19]，見圖2.預(yù)測相對誤差見圖3.由圖3可知，模型訓(xùn)練集的相對誤差均小于5%，證明該預(yù)測模型具有良好的預(yù)測準(zhǔn)確度.

圖2 訓(xùn)練樣本瓦斯涌出量預(yù)測值和實(shí)測值對比Fig.2 comparison between predicted value and measured value of gas emission in training samples

圖3 訓(xùn)練樣本預(yù)測值相對誤差Fig.3 relative error of training sample prediction values

測試樣本中預(yù)測值與實(shí)際值的對比見圖4.由圖4可知，模型對瓦斯涌出量的預(yù)測精度良好.

圖4 測試樣本瓦斯涌出量預(yù)測值與實(shí)測值對比Fig.4 comparison between the predicted value and the measured value of the gas emission of the test sample

選取5種預(yù)測模型進(jìn)行對比分析，結(jié)果見表7，可知，基于FA-IGA-BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測結(jié)果與實(shí)際的情況更為吻合.樣本16的相對誤差雖然高于IGA-BP神經(jīng)網(wǎng)絡(luò)模型，但誤差在可接受的范圍內(nèi).FA-IGA-BP模型最大相對誤差為3.794%，平均誤差為2.945%，低于IGA-BP模型的3.338%、GA-BP模型的5.007%、GA模型的6.271%、BP模型的6.179%.

表7 瓦斯涌出量預(yù)測結(jié)果及誤差Tab.7 prediction results and error comparison of test samples

平均相對變化量EARV可以用來表示神經(jīng)網(wǎng)絡(luò)的泛化能力，值越小泛化能力越強(qiáng)，計(jì)算公式為

式中，Y(i)、y(i)、分別為實(shí)際值、預(yù)測值和平均值.

通過計(jì)算得到FA-IGA-BP神經(jīng)網(wǎng)絡(luò)模型的平均相對變化量為0.0137，低于IGA-BP神經(jīng)網(wǎng)絡(luò)的0.0228、GA-BP神經(jīng)網(wǎng)絡(luò)的0.0405、GA模型的0.0823、BP模型的0.0554，表明該模型的泛化能力較強(qiáng).

此外，均方根誤差ERMSE具有對預(yù)測值中的異常誤差反應(yīng)敏銳的特點(diǎn)[20]，能良好地反映預(yù)測的準(zhǔn)確程度，計(jì)算公式為

式中，Y(i)為實(shí)際值；y(i)為預(yù)測值；n為預(yù)測值個(gè)數(shù).

利用均方根誤差公式計(jì)算得出，F(xiàn)A-IGA-BP模型的預(yù)測均方根誤差20.050%，而IGA-BP模型的預(yù)測均方根誤差為25.800%，GA-BP模型的預(yù)測均方根誤差為34.442%，GA模型預(yù)測均方根誤差為40.898%，BP模型均方根誤差為40.604%，說明基于FA-IGA- BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測更加準(zhǔn)確.

4 結(jié)論

（1）綜合選取影響瓦斯涌出量的12個(gè)因素，應(yīng)用因子分析對相關(guān)性較強(qiáng)的8個(gè)影響因素進(jìn)行因子分析，提取出兩個(gè)公共因子f1和f2，有效減少判別因子之間的數(shù)據(jù)冗余，使算法的泛化能力得到增強(qiáng).

（2）通過改進(jìn)遺傳算法（IGA）對BP神經(jīng)網(wǎng)絡(luò)的權(quán)重和閾值進(jìn)行尋優(yōu)，優(yōu)化神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)及算法的全局收斂性，建立最佳預(yù)測模型，選用實(shí)測數(shù)據(jù)進(jìn)行預(yù)測分析，同時(shí)與其他預(yù)測模型對比.結(jié)果表明，F(xiàn)A-IGA-BP模型的預(yù)測準(zhǔn)確度較高，最大相對誤差為3.794%，平均誤差為2.945%，能滿足回采工作面瓦斯預(yù)測的工程需求.

（3）選取的瓦斯涌出的影響因素還存在一定的局限性，對人為影響因素的選取及分析較少，下一步工作需要全面分析影響瓦斯涌出的因素，進(jìn)一步提高預(yù)測的準(zhǔn)確性.