邵良杉,王 振
(遼寧工程技術(shù)大學(xué) 系統(tǒng)工程研究所,遼寧 葫蘆島 125105)
瓦斯是煤礦安全的主要隱患之一,中國(guó)瓦斯事故占到煤礦事故的80%以上.準(zhǔn)確預(yù)測(cè)瓦斯涌出,不僅是瓦斯防護(hù)和治理、煤礦通風(fēng)系統(tǒng)設(shè)計(jì)的重要基礎(chǔ)[1],也為礦井高質(zhì)量生產(chǎn)提供了安全保障.隨著工業(yè)化、信息化的發(fā)展,機(jī)械程度不斷提高,礦井開采深度不斷增加,工作面的產(chǎn)能加重,隨之帶來工作面瓦斯涌出量的上升,制約著煤礦的安全生產(chǎn)和經(jīng)濟(jì)效益[2].因此研究一種瓦斯涌出量的預(yù)測(cè)模型顯得極為重要.
近年來,許多學(xué)者針對(duì)瓦斯涌出量預(yù)測(cè)進(jìn)行了研究.呂伏等[3]提出將主成分分析法應(yīng)用于采掘工作面瓦斯涌出的多步線性回歸預(yù)測(cè)中.萬仁寶等[4]將灰色系統(tǒng)與神經(jīng)網(wǎng)絡(luò)相結(jié)合,在小信息量的情況下做出準(zhǔn)確預(yù)測(cè).付華等[5]通過優(yōu)化粒子群和Elman神經(jīng)網(wǎng)絡(luò)的權(quán)重和閾值建立絕對(duì)瓦斯涌出量預(yù)測(cè)系統(tǒng)模型.施式亮等[6]通過分解瓦斯實(shí)時(shí)涌出數(shù)據(jù),通過粒子群算法優(yōu)化支持向量機(jī)(SVM),利用優(yōu)化過的SVM預(yù)測(cè)多個(gè)本征模函數(shù),并對(duì)結(jié)果進(jìn)行累加從而得出瓦斯預(yù)測(cè)值.溫廷新等[7]運(yùn)用灰色關(guān)聯(lián)和因子分析對(duì)指標(biāo)進(jìn)行分析提取,利用量子遺傳算法對(duì)最小二乘支持向量機(jī)進(jìn)行了優(yōu)化,建立了煤與瓦斯突出預(yù)測(cè)模型.林海飛等[8]為提高煤層瓦斯含量預(yù)測(cè)的科學(xué)性及準(zhǔn)確性,提出基于粒子群優(yōu)化誤差反向傳播神經(jīng)網(wǎng)絡(luò)的瓦斯含量預(yù)測(cè)模型.
以上研究存在一定局限性,瓦斯涌出情況多變且不穩(wěn)定,傳統(tǒng)預(yù)測(cè)方法的預(yù)測(cè)結(jié)果與實(shí)際情況相差較遠(yuǎn),因此需要更為優(yōu)化的預(yù)測(cè)方法.本文采用因子分析去除數(shù)據(jù)相關(guān)性并濃縮數(shù)據(jù)的優(yōu)點(diǎn),運(yùn)用改進(jìn)遺傳算法,對(duì)BP網(wǎng)絡(luò)進(jìn)行優(yōu)化,建立預(yù)測(cè)回采工作面瓦斯涌出量的優(yōu)化模型,并通過對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行預(yù)測(cè)和誤差分析驗(yàn)證該模型的準(zhǔn)確性.
因子分析(factor analysis,F(xiàn)A)通過研究變量之間的內(nèi)部依賴性,對(duì)具有錯(cuò)綜復(fù)雜關(guān)系的變量進(jìn)行簡(jiǎ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),待觀測(cè)變量通過m(m<k)個(gè)兩兩正交的公共因子F= [f1,f2,… ,fm]T表示,其數(shù)學(xué)表達(dá)式為
式中,每個(gè)變量都是標(biāo)準(zhǔn)化的變量,iε為特殊因子,滿足E(εi) = 0,作用于其對(duì)應(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是對(duì)應(yīng)的征根1λ的任意一個(gè)模為1的特征向量,記;則可證也是所求.同理可求α3,… ,αm.進(jìn)一步解因子值,可得R為X的相關(guān)系數(shù)矩陣,令βj=R-1αj,通過Fi=Xβj計(jì)算,得到公共因子值.
近年來,遺傳算法[10](genetic algorithm,GA)應(yīng)用廣泛,但在應(yīng)對(duì)局部搜索時(shí),有大概率陷入“早熟”的問題,影響其收斂性.針對(duì)該問題,對(duì)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.
對(duì)懲罰系數(shù)γ進(jìn)一步優(yōu)化為
式中,t為迭代次數(shù),當(dāng)t較小時(shí),γ接近于γmax,保證了算法的全局搜索能力;隨著t的增大,γ以非線性減小,從而保障了算法的局部求解能力.
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對(duì)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對(duì)其中的權(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.
經(jīng)過實(shí)際現(xiàn)場(chǎng)調(diào)研及相關(guān)資料的參考,選擇12個(gè)因素作為回采工作面瓦斯涌出量的主要影響因素[11]:原始瓦斯量X1、煤層埋藏深度X2、煤層厚度X3、傾角X4、工作面長(zhǎng)度X5、推進(jìn)速度X6、采出率X7、臨近層瓦斯量X8、臨近層厚度X9、臨近層間距X10、巖性X11、開采深度X12.
參照文獻(xiàn)[12]提供的某礦18組回采工作面瓦斯涌出量數(shù)據(jù),對(duì)以上提及的12個(gè)因素進(jìn)行Pearson相關(guān)性的分析.Pearson系數(shù)及Sig.(2-tailed)檢驗(yàn)結(jié)果表明:煤層傾角、工作面長(zhǎng)度、臨近厚度、臨近層間距與其他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)并對(duì)其進(jìn)行因子分析,從而得到少數(shù)公共因子,實(shí)現(xiàn)對(duì)回采工作面瓦斯涌出量的影響因素降維.為保證指標(biāo)選取的可靠性,通過KMO值和巴特利球體檢驗(yàn)對(duì)8個(gè)高度相關(guān)的指標(biāo)檢測(cè):樣本KMO值為0.8大于0.6,屬于允許域;巴特利球度檢驗(yàn)近似卡方為164.193,顯著性概率p為0.000,小于0.05;可知在各個(gè)指標(biāo)數(shù)據(jù)間具有一定的相關(guān)性并能夠進(jìn)行因子分析.
將實(shí)測(cè)數(shù)據(jù)中原始瓦斯量、煤層埋藏深度、厚度、推速、采出率、臨近層瓦斯量、巖性、開采深度8個(gè)指標(biāo)作為因子分析的變量,采用SPSS軟件對(duì)回采工作面瓦斯涌出量的影響因素進(jìn)行回歸分析,在分析方法上,采用主成分提取法,得出公共因子碎石圖,見圖1.
圖1 公共因子碎石圖Fig.1 common factor scree plot
由圖1可以看出,前2個(gè)成分的特征值有較為明顯坡度變化,而自第3個(gè)成分的特征值開始,其變化幅度快速地趨于穩(wěn)定.由碎石準(zhǔn)則[13]可見,所取的兩個(gè)公共因子能夠更加明顯地實(shí)現(xiàn)對(duì)原變量方差的解釋.
總方差解釋見表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)的綜合性,故需對(duì)所取因子做一定處理從而降低初始因子的綜合性.本文選用正交旋轉(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
由于煤層傾角X4、工作面長(zhǎng)度X5、臨近厚度X9及臨近層間距X10與另外8個(gè)變量之間的相關(guān)性較弱,因此將煤層傾角、工作面長(zhǎng)度、臨近厚度、臨近層間距,以及所選取的2個(gè)公共因子作為IGA-BP神經(jīng)網(wǎng)絡(luò)的特征因子.為便于數(shù)據(jù)處理,加快訓(xùn)練網(wǎng)絡(luò)的收斂速度,對(duì)數(shù)據(jù)進(jìn)行規(guī)范化[14],本文采用Z-score歸一化方法對(duì)其他指標(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
(1)參數(shù)設(shè)置
瓦斯涌出量預(yù)測(cè)結(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ì)延長(zhǎng)訓(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],測(cè)試隱含層節(jié)點(diǎn)從3到13的仿真誤差,對(duì)其擬合后的數(shù)據(jù)特征進(jìn)行比較,結(jié)果見表5.
表5 隱含層節(jié)點(diǎn)數(shù)目測(cè)試結(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ò)的測(cè)試集,并在Matlab R2016b中建立基于FA-IGA-BP算法的回采工作面瓦斯涌出量預(yù)測(cè)模型.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對(duì)BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練是良好的,可以得到一系列適用于回采工作面瓦斯涌出量預(yù)測(cè)的最優(yōu)權(quán)值和閾值,數(shù)據(jù)見表6.
表6 IGA得到的最優(yōu)權(quán)值和閾值Tab.6 optimal weights and thresholds obtained by IGA
(2)預(yù)測(cè)結(jié)果分析
采用回代估計(jì)方法得到訓(xùn)練集瓦斯涌出量的預(yù)測(cè)值,并與實(shí)測(cè)值比較[18-19],見圖2.預(yù)測(cè)相對(duì)誤差見圖3.由圖3可知,模型訓(xùn)練集的相對(duì)誤差均小于5%,證明該預(yù)測(cè)模型具有良好的預(yù)測(cè)準(zhǔn)確度.
圖2 訓(xùn)練樣本瓦斯涌出量預(yù)測(cè)值和實(shí)測(cè)值對(duì)比Fig.2 comparison between predicted value and measured value of gas emission in training samples
圖3 訓(xùn)練樣本預(yù)測(cè)值相對(duì)誤差Fig.3 relative error of training sample prediction values
測(cè)試樣本中預(yù)測(cè)值與實(shí)際值的對(duì)比見圖4.由圖4可知,模型對(duì)瓦斯涌出量的預(yù)測(cè)精度良好.
圖4 測(cè)試樣本瓦斯涌出量預(yù)測(cè)值與實(shí)測(cè)值對(duì)比Fig.4 comparison between the predicted value and the measured value of the gas emission of the test sample
選取5種預(yù)測(cè)模型進(jìn)行對(duì)比分析,結(jié)果見表7,可知,基于FA-IGA-BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測(cè)結(jié)果與實(shí)際的情況更為吻合.樣本16的相對(duì)誤差雖然高于IGA-BP神經(jīng)網(wǎng)絡(luò)模型,但誤差在可接受的范圍內(nèi).FA-IGA-BP模型最大相對(duì)誤差為3.794%,平均誤差為2.945%,低于IGA-BP模型的3.338%、GA-BP模型的5.007%、GA模型的6.271%、BP模型的6.179%.
表7 瓦斯涌出量預(yù)測(cè)結(jié)果及誤差Tab.7 prediction results and error comparison of test samples
平均相對(duì)變化量EARV可以用來表示神經(jīng)網(wǎng)絡(luò)的泛化能力,值越小泛化能力越強(qiáng),計(jì)算公式為
式中,Y(i)、y(i)、分別為實(shí)際值、預(yù)測(cè)值和平均值.
通過計(jì)算得到FA-IGA-BP神經(jīng)網(wǎng)絡(luò)模型的平均相對(duì)變化量為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具有對(duì)預(yù)測(cè)值中的異常誤差反應(yīng)敏銳的特點(diǎn)[20],能良好地反映預(yù)測(cè)的準(zhǔn)確程度,計(jì)算公式為
式中,Y(i)為實(shí)際值;y(i)為預(yù)測(cè)值;n為預(yù)測(cè)值個(gè)數(shù).
利用均方根誤差公式計(jì)算得出,F(xiàn)A-IGA-BP模型的預(yù)測(cè)均方根誤差20.050%,而IGA-BP模型的預(yù)測(cè)均方根誤差為25.800%,GA-BP模型的預(yù)測(cè)均方根誤差為34.442%,GA模型預(yù)測(cè)均方根誤差為40.898%,BP模型均方根誤差為40.604%,說明基于FA-IGA- BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測(cè)更加準(zhǔn)確.
(1)綜合選取影響瓦斯涌出量的12個(gè)因素,應(yīng)用因子分析對(duì)相關(guān)性較強(qiáng)的8個(gè)影響因素進(jìn)行因子分析,提取出兩個(gè)公共因子f1和f2,有效減少判別因子之間的數(shù)據(jù)冗余,使算法的泛化能力得到增強(qiáng).
(2)通過改進(jìn)遺傳算法(IGA)對(duì)BP神經(jīng)網(wǎng)絡(luò)的權(quán)重和閾值進(jìn)行尋優(yōu),優(yōu)化神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)及算法的全局收斂性,建立最佳預(yù)測(cè)模型,選用實(shí)測(cè)數(shù)據(jù)進(jìn)行預(yù)測(cè)分析,同時(shí)與其他預(yù)測(cè)模型對(duì)比.結(jié)果表明,F(xiàn)A-IGA-BP模型的預(yù)測(cè)準(zhǔn)確度較高,最大相對(duì)誤差為3.794%,平均誤差為2.945%,能滿足回采工作面瓦斯預(yù)測(cè)的工程需求.
(3)選取的瓦斯涌出的影響因素還存在一定的局限性,對(duì)人為影響因素的選取及分析較少,下一步工作需要全面分析影響瓦斯涌出的因素,進(jìn)一步提高預(yù)測(cè)的準(zhǔn)確性.