馬 磊，陸衛(wèi)東，魏國營,,3

(1.河南理工大學(xué) 安全科學(xué)與工程學(xué)院，河南 焦作 454000；2.新疆工程學(xué)院 安全科學(xué)與工程學(xué)院，新疆 烏魯木齊830023；3.煤炭安全生產(chǎn)與清潔高效利用省部共建協(xié)同創(chuàng)新中心，河南 焦作 454000 )

0 引言

煤層瓦斯含量是描述單位質(zhì)量煤容納瓦斯體積的量，是預(yù)測瓦斯涌出量和煤層突出危險性的重要參數(shù)，瓦斯含量的精準(zhǔn)預(yù)測對瓦斯災(zāi)害治理具有重要意義[1-2]。

近年來，隨著機器學(xué)習(xí)和智能算法的快速發(fā)展，人工神經(jīng)網(wǎng)絡(luò)、支持向量機和極限學(xué)習(xí)機等方法在煤層瓦斯預(yù)測中得到了廣泛應(yīng)用[3-5]。魏國營等[6]利用主成分分析法(Principle Component Analysis,PCA)提取瓦斯含量影響因素的主成分，并用自適應(yīng)混合粒子群算法優(yōu)化支持向量回歸機(Vector Regression Machine,SVM)的初始參數(shù)構(gòu)建PCA-AHPSO-SVR瓦斯含量預(yù)測模型；林海飛等[7]利用粒子群優(yōu)化算法(Particle Swarm Optimization,PSO)優(yōu)化BP神經(jīng)網(wǎng)絡(luò)(Back Propagation Neural Network,BPNN)的初始權(quán)值和閾值，構(gòu)建PSO-BPNN瓦斯含量預(yù)測模型；汪吉林等[8]將構(gòu)造指數(shù)、煤層埋深、煤層厚度和原始瓦斯含量作為模型輸入，構(gòu)建BPNN瓦斯含量損失量預(yù)測模型；曹博等[9]采用PCA法降低數(shù)據(jù)維度，利用遺傳算法(Genetic Algorithm,GA)搜索BPNN的初始參數(shù)，構(gòu)建GA-BPNN瓦斯含量預(yù)測模型；溫廷新等[10]利用改進(jìn)的果蠅優(yōu)化算法優(yōu)選極限學(xué)習(xí)機(Extreme Learning Machine，ELM)的參數(shù)建立IFOA-ELM瓦斯預(yù)測模型；付華等[11]利用蟻群聚類(Ant Colony Clustering，ACC)算法優(yōu)化Elman神經(jīng)網(wǎng)絡(luò)建立絕對瓦斯涌出量預(yù)測模型。對BPNN易陷入局部極小的陷阱且收斂速度慢的問題，學(xué)者都利用改進(jìn)了的單一算法優(yōu)化模型初始參數(shù)，但瓦斯含量影響因素多且各影響因素之間非線性關(guān)系復(fù)雜，再加上單一智能算法的改進(jìn)程度的局限性導(dǎo)致預(yù)測模型穩(wěn)定性和準(zhǔn)確率低。為解決這類問題，將GA和具有時變且最終趨于零概率突跳性的模擬退火(Simulated Annealing，SA)算法整合成GASA混合優(yōu)化算法，該算法吸取2類算法的優(yōu)點，克服GA易早熟的缺點，提升BPNN的尋優(yōu)訓(xùn)練速度和最優(yōu)解質(zhì)量。

本文擬利用灰色關(guān)聯(lián)分析法篩選瓦斯含量的地質(zhì)影響因素，建立預(yù)測參數(shù)體系，將GA和SA整合為GASA混合算法協(xié)同初始化BPNN參數(shù)建立GASA-BPNN瓦斯含量預(yù)測模型，以期得到更高效和精準(zhǔn)的預(yù)測模型。

1 建立模型

1.1 相關(guān)理論

GA是模擬自然界遺傳機制和達(dá)爾文生物進(jìn)化論的并行隨機搜索算法[12]。GA將待優(yōu)化的參數(shù)編碼為若干個個體形成群體，通過遺傳操作來效仿生物的進(jìn)化過程，利用適應(yīng)度值來評價個體的質(zhì)量，通過循環(huán)迭代和適應(yīng)度函數(shù)引導(dǎo)種群共同進(jìn)化，使最優(yōu)個體逼近問題全局最優(yōu)解。

SA[13]是基于Monte Carlo迭代求解策略的1種隨機搜索算法，出發(fā)點是基于固體物質(zhì)的退火過程與組合優(yōu)化問題之間的相似性。SA由初始高溫開始，利用具有概率突跳性的Metropolis抽樣準(zhǔn)則以一定概率接受鄰域中次優(yōu)解隨機搜索，隨著溫度不斷下降，最終找到全局最優(yōu)解。

依據(jù)SA在搜索過程中逐漸趨于0的概率突跳性和GA優(yōu)勝劣汰的優(yōu)勢，將GA和SA整合成GASA混合算法，避免算法落入局部極小的陷阱；在結(jié)構(gòu)方面，GA是依據(jù)鄰域函數(shù)對群體并行搜索，而SA每次搜索時只針對1個個體，GASA混合算法相當(dāng)于在GA中群體的每個個體下均串行1個SA，使個體的鄰域搜索結(jié)構(gòu)更加豐富，增強算法搜索能力和效率。GASA算法優(yōu)化流程包括以下10個步驟[14]：

1)初始化GA。設(shè)定種群規(guī)模N，初始化種群PK，設(shè)定最大遺傳代數(shù)M，遺傳迭代次數(shù)K=1。

2)定義適應(yīng)度函數(shù)。

3)評價PK適應(yīng)度值。判斷是否滿足遺傳停機條件？若滿足，選擇適應(yīng)度最大的個體作為最優(yōu)解輸出，算法結(jié)束；若不滿足，執(zhí)行步驟4)～9)。

4)對種群PK執(zhí)行遺傳交叉和變異生成種群PK0，并評價PK0適應(yīng)度值。

5)初始化SA參數(shù)。設(shè)定種群PK0為SA的初始解，令i=0，設(shè)定初溫T=Ti(足夠高)，確定每個T狀態(tài)下Metropolis鏈長L，令鏈中迭代次數(shù)Q=1。

6)對當(dāng)前溫度T和Q=1,2,3,…L狀態(tài)下，重復(fù)執(zhí)行步驟7)～9)。

7)對當(dāng)前群體的每個個體隨機擾動產(chǎn)生新的群體，評價新群體適應(yīng)度。

8)計算新群體和當(dāng)前群體各個體適應(yīng)度值的差值ΔE。若-ΔE<0，則接受新個體；若-ΔE≥0，計算接受概率，即判斷P=exp(-ΔE/T)≥rand是否成立，若成立，也接受新解；否則拒絕新解。

9)令Q=Q+1，判斷條件Q>L是否成立，若成立，則i=i+1，按溫度衰減函數(shù)進(jìn)行降溫，即T=Ti+1，其中TL不成立，則轉(zhuǎn)至步驟7)。SA停機條件通常為設(shè)置停機溫度。

10)PK種群經(jīng)遺傳操作產(chǎn)生種群PK0，種群PK0作為SA的初始解經(jīng)退火操作產(chǎn)生PK+1，若新種群PK+1沒有滿足遺傳停機條件，則由SA產(chǎn)生的種群PK+1作為GA的初始種群再參加迭代尋優(yōu)。

1.2 建立GASA-BPNN預(yù)測模型

BPNN是由誤差逆向傳播算法訓(xùn)練的多層前饋神經(jīng)網(wǎng)絡(luò)，在結(jié)構(gòu)上由輸入層、若干隱含層和輸出層組成，各層采用權(quán)重參數(shù)全連接方式，其核心是信息從輸入層經(jīng)隱含層處理后傳遞到輸出層計算誤差，若誤差不滿足要求，將誤差信號反向傳遞給各節(jié)點，節(jié)點更新參數(shù)，反復(fù)迭代直至滿足要求，最終構(gòu)建能逼近復(fù)雜非線性映射的模型[15]。但BPNN初始參數(shù)對模型性能的影響十分敏感，隨機初始化參數(shù)是導(dǎo)致BPNN陷入局部極值、收斂速度慢和模型精度低的主要原因。因此，利用GASA算法初始化BPNN的參數(shù)來代替?zhèn)鹘y(tǒng)方法，有效地提高BPNN訓(xùn)練速度和預(yù)測精度，GASA算法初始化BPNN參數(shù)具體流程如圖1所示。

圖1 GASA初始化BPNN參數(shù)流程Fig.1 Procedure of GASA to initialize BPNN parameters

2 現(xiàn)場應(yīng)用

2.1 瓦斯含量影響因素分析及建立數(shù)據(jù)集

根據(jù)焦作礦區(qū)某礦瓦斯地質(zhì)條件和瓦斯分布規(guī)律研究情況[16]，瓦斯含量(X0,m3/t)影響因素選取可量化處理的以下8個因素：煤層埋深(X1,m)、煤層厚度(X2,m)、傾角系數(shù)(X3)、上覆基巖厚度(X4,m)、圍巖等效系數(shù)(X5)、斷層復(fù)雜系數(shù)(X6)、褶皺復(fù)雜系數(shù)(X7)和底板標(biāo)高(X8,m)[6]。在該礦選取300組瓦斯含量及影響因素的數(shù)據(jù)作為實驗數(shù)據(jù)集，訓(xùn)練集和測試集利用10折交叉驗證的方式劃分，部分?jǐn)?shù)據(jù)見表1。

表1 瓦斯含量及影響因素數(shù)據(jù)Table 1 Data of gas contents and influencing factors

2.2 確定BPNN結(jié)構(gòu)

利用灰色關(guān)聯(lián)分析法篩選主控因素為模型的輸入?；疑P(guān)聯(lián)計算步驟如下[17]：

1)確定參考序列和比較序列：參考序列為瓦斯含量(X0)，比較序列為瓦斯含量的8個影響因素。

2)數(shù)據(jù)按公式(1)進(jìn)行歸一化處理。

(1)

式中：xl(h)為樣本編號為h的第l評價指標(biāo)的數(shù)值；h為樣本編號，h=1,2,…,n；l為評價指標(biāo)，即瓦斯含量的各影響因素，l=1,2,…,m；xmax和xmin分別是評價指標(biāo)的最大和最小值。

3)按公式(2)計算關(guān)聯(lián)度系數(shù)。

(2)

式中：ξe(k)為比較序列xe對參考序列x0在第k指標(biāo)上的關(guān)聯(lián)系數(shù)，k=1,2,…,n；ρ∈[0,1]為分辨系數(shù)。

4)按公式(3)計算關(guān)聯(lián)度。

(3)

式中：n為樣本數(shù)量；wk為指標(biāo)權(quán)重。

關(guān)聯(lián)度計算結(jié)果見表2。由表2可知：影響因素X4和X5的關(guān)聯(lián)度相近，說明這2個因素對含量影響效果一樣，只需保留1個即可說明效果，因此，剔除X5。因素X3和X7的關(guān)聯(lián)度不足0.5，屬低相關(guān)度，因此剔除，從而確定輸入層節(jié)點數(shù)為5。

表2 灰色關(guān)聯(lián)分析結(jié)果Table 2 Results of grey correlation analysis

BPNN隱含層節(jié)點數(shù)過少會造成參數(shù)無法得到充分訓(xùn)練和預(yù)測精度低的問題，節(jié)點過多會造成過擬合現(xiàn)象。依據(jù)經(jīng)驗公式(4)確定隱含層節(jié)點數(shù)范圍為[4,12][18]：

(4)

式中：W為隱含層節(jié)點數(shù)；c是輸入層節(jié)點數(shù)；q為輸出層節(jié)點數(shù)；B是1～10范圍內(nèi)的整數(shù)。

通過比較不同節(jié)點數(shù)模型的誤差試湊出最佳隱含層節(jié)點數(shù)，結(jié)果見表3。由表3可知：當(dāng)W=12時誤差最小，因此，網(wǎng)絡(luò)隱含層節(jié)點確定為12。研究表明三層BPNN能任意逼近任何復(fù)雜非線性系統(tǒng)，因此，GASA-BPNN預(yù)測模型的結(jié)構(gòu)被確定為5-12-1型BPNN。

表3 隱含層節(jié)點數(shù)與平均相對誤差Table 3 Numbers of hidden layer nodes and average relative errors

2.3 模型應(yīng)用

模型應(yīng)用具體流程如下[14]：

1)劃分測試集與訓(xùn)練集。首先將數(shù)據(jù)按公式(1)歸一化到[0,1]范圍內(nèi)。然后將數(shù)據(jù)集分層抽樣劃分10個互斥子集，每次利用某個子集作為測試集，其余作為訓(xùn)練集進(jìn)行10組訓(xùn)練和測試，重復(fù)10次，10折交叉驗證示意如圖2所示[19]。

圖2 10折交叉驗證示意Fig.2 10-fold cross-validation diagram

2)GA種群初始化。采用浮點數(shù)的編碼規(guī)則對85個參數(shù)進(jìn)行編碼，每個個體代表1套權(quán)值和閾值，并隨機生成50個個體的種群PK。

3)GA參數(shù)設(shè)定。設(shè)定遺傳最大迭代次數(shù)為300，迭代次數(shù)K=1，當(dāng)滿足最大遺傳迭代次數(shù)時遺傳停機，交叉概率pc=0.7，變異概率pm=0.005。

4)定義適應(yīng)度函數(shù)。適應(yīng)度函數(shù)根據(jù)預(yù)測均方誤差(Mean Square Error，MSE)進(jìn)行設(shè)定，MSE越小，個體適應(yīng)度值越大。

5)評價初始種群PK每個個體的適應(yīng)度，依據(jù)適應(yīng)度值進(jìn)行排序。

6)GASA停機判定。判斷是否滿足最大遺傳迭代次數(shù)，若滿足，則優(yōu)化結(jié)束，轉(zhuǎn)至步驟7)～11)；否則，K=K+1，轉(zhuǎn)步驟12)。

7)初始化BPNN。將步驟6)優(yōu)化后的解賦值給BPNN，設(shè)定最大訓(xùn)練次數(shù)為H=750，記錄迭代次數(shù)為S=0，隱含層和輸出層的激活函數(shù)分別選擇logsig函數(shù)和purelin函數(shù)。

8)輸入數(shù)據(jù)集。數(shù)據(jù)按照BPNN的正向向前傳輸，第f個隱含層神經(jīng)元的輸入如式(5)所示：

(5)

式中：γf為第f個隱含層神經(jīng)元的輸入；vvf是輸入層第v個神經(jīng)元和隱含層第f個神經(jīng)元之間連接的權(quán)值。

第j個輸出節(jié)點的輸入如式(6)所示：

(6)

式中：βj為第j個輸出神經(jīng)元的輸入；wfj為隱含層第f個神經(jīng)元與輸出層第j個神經(jīng)元的連接權(quán)值；bf為隱含層第f個神經(jīng)元的輸出。

9)更新BPNN權(quán)值和閾值。預(yù)測均方誤差反向傳播給各個節(jié)點，節(jié)點誤差自適應(yīng)學(xué)習(xí)率梯度下降法更新參數(shù)如式(7)所示：

(7)

式中：μ(S+1)為第S+1代的學(xué)習(xí)率；μ(S)為第S代的學(xué)習(xí)率；Sin為增量系數(shù)；Sde為減量系數(shù)；E(S+1)為第S+1代的誤差；E(S)為第S代的誤差。

當(dāng)誤差以較少的趨勢逼近目標(biāo)值時，表明模型訓(xùn)練方向正確，可增大學(xué)習(xí)速度；當(dāng)誤差的增加超過允許范圍時，表明模型訓(xùn)練方向錯誤，應(yīng)降低學(xué)習(xí)率。

10)BPNN停機判定。判斷是否滿足停機條件(設(shè)定最小誤差為停機條件)，若滿足，網(wǎng)絡(luò)結(jié)束學(xué)習(xí)完成模型建立，轉(zhuǎn)步驟11)；否則判斷是否H=S，若是則返回7)；若不是，則S=S+1返步驟8)。

11)模型仿真及評價。將測試集輸入GASA-BPNN模型測試，利用MSE、迭代次數(shù)和預(yù)測相對誤差進(jìn)行模型性能評價。

12)GA選擇。采用輪盤賭選擇法篩選個體，個體xu的適應(yīng)度值為f(xu)，則個體xu被選擇的概率為式(8)：

(8)

式中：N為群體中個體數(shù)目總和。

(9)

式中：λ1+λ2=1，λ1>0，λ2>0。

14)GA變異。隨機選擇個體基因型X=x1x2…xb…xS，將xb作為突變點，按公式(10)進(jìn)行遺傳運算。

(10)

式中：UB,LB分別為xb的上、下邊界值；r為隨機數(shù)；K為進(jìn)化代數(shù)；Δ(t,y)的定義為式(11)：

(11)

式中：q為控制非一致性參數(shù)，取0.8。

15)評價PK0適應(yīng)度。評價PK0的適應(yīng)度值，設(shè)定PK0為退火算法的初始種群。

16)SA參數(shù)初始化。將所有經(jīng)過遺傳運算的個體進(jìn)行退火。設(shè)定初始溫度T=100 ℃，冷卻系數(shù)0.98，停機準(zhǔn)則，令Q=1，馬爾可夫鏈長L=60。

17)對群體中的每個個體進(jìn)行擾動產(chǎn)生1個新解，計算新舊群體每個個體的適度值的差值，根據(jù)Metropolis抽樣準(zhǔn)則判斷是否接受新解。

18)令Q=Q+1，判斷是否在馬爾可夫鏈中，若在則轉(zhuǎn)步驟17)，若不在鏈中，則退溫，再判斷是否達(dá)到SA停機條件，如滿足停機，則SA算法結(jié)束，轉(zhuǎn)步驟4)；否則，Q=1，轉(zhuǎn)步驟17)。

2.4 模型性能對比與結(jié)果分析

為了驗證GASA-BPNN瓦斯含量的預(yù)測效果，在同等條件下分別建立BPNN和GA-BPNN模型對比分析。首先利用GA算法和GASA算法分別初始化BPNN的參數(shù)30次，最優(yōu)個體的平均適應(yīng)度變化曲線如圖3所示，由圖3可得GASA算法和GA優(yōu)化過程平均適應(yīng)度分別在第92次和第177次迭代時完成尋優(yōu)，GASA算法尋優(yōu)速度明顯快于GA算法并且最優(yōu)適應(yīng)度大于GA算法的最優(yōu)適應(yīng)度，說明GASA算法在參數(shù)尋優(yōu)過程中效率和效果均優(yōu)于GA算法。

圖3 平均最優(yōu)個體適應(yīng)度值曲線Fig.3 Curves of average optimal individual fitness value

其次，將經(jīng)過10次10折交叉驗證劃分的100組訓(xùn)練集分別輸入BPNN模型、GA-BPNN模型和GASA-BPNN模型中按照誤差自適應(yīng)學(xué)習(xí)率梯度下降法進(jìn)行參數(shù)學(xué)習(xí)，3種模型分別訓(xùn)練100次，100次訓(xùn)練平均均方誤差曲線如圖4所示，圖4顯示GASA-BPNN和GA-BPNN分別在第332次和第459次訓(xùn)練時達(dá)到目標(biāo)誤差線，而BPNN在第736次訓(xùn)練時仍未達(dá)到目標(biāo)誤差，說明經(jīng)GASA優(yōu)化后的模型對樣本學(xué)習(xí)所需的訓(xùn)練次數(shù)明顯減少并且權(quán)值和閾值更接近最優(yōu)參數(shù)。

圖4 平均均方誤差收斂曲線Fig.4 Curves of average mean square error convergence

最后，將100個測試樣本集分別輸入3種模型中預(yù)測計算，各組預(yù)測平均相對誤差見表4，從表4可知，各組GASA-BPNN模型預(yù)測的平均相對誤差和10組總平均相對誤差明顯低于其他模型。

表4 各組測試集平均相對誤差Table 4 Average relative error of each group of test set %

選取樣本編號為1到30的30個樣本集作為研究對象，分析30個樣本在10組仿真測試中均是測試集時的預(yù)測性能，得到實測值與平均預(yù)測值對比如圖5所示，平均預(yù)測相對誤差對比如圖6所示。由圖5可得GASA-BPNN的預(yù)測值較GA-BPNN和BPNN最接近實測值，由圖6可得GASA-BPNN模型預(yù)測的平均相對誤差最小且相對誤差波動最小，說明GASA-BPNN模型預(yù)測穩(wěn)定性強并且預(yù)測準(zhǔn)確率更高，能滿足預(yù)測目標(biāo)和工程需求。

圖5 實測值與平均預(yù)測值對比Fig.5 Comparison of measured values and average predicted values

圖6 平均預(yù)測相對誤差對比Fig.6 Comparison of average prediction relative error

3 結(jié)論

1)將GA和具有時變概率突跳性的SA整合為GASA算法，該策略突破了單一智能算法改進(jìn)程度的局限性，相同環(huán)境下該算法在BPNN初始參數(shù)尋優(yōu)過程中效率和效果均優(yōu)于GA算法。

2)利用灰色關(guān)聯(lián)分析法消除瓦斯含量影響因素相關(guān)性，建立瓦斯含量預(yù)測參數(shù)體系從而確定BPNN輸入層節(jié)點個數(shù)，減少數(shù)據(jù)維度，降低模型訓(xùn)練的計算量和預(yù)測誤差。

3)利用GASA算法優(yōu)化BPNN的權(quán)值和閾值，有效地解決BPNN收斂速度慢和易陷入極小的問題，與其他模型相比，經(jīng)過GASA優(yōu)化的BPNN對數(shù)據(jù)的學(xué)習(xí)迭代次數(shù)明顯減少并且GASA-BPNN模型能更加準(zhǔn)確地預(yù)測瓦斯含量，對煤礦瓦斯治理具有重要意義。