雷昌奎，江莉娟，鄧存寶，鄧 軍，馬 礪，王偉峰，張勇敢

（1.太原理工大學 安全與應急管理工程學院，山西 太原 030024；2.陜西山利科技發(fā)展有限責任公司，陜西 西安 710075；3.匈牙利佩奇大學 商業(yè)與經(jīng)濟學院，匈牙利 佩奇 7626；4.西安科技大學 安全科學與工程學院，陜西 西安 710054）

煤自燃是一個復雜的動態(tài)氧化過程，是煤氧復合作用的結(jié)果，嚴重威脅我國礦井的安全生產(chǎn)[1-3]。研究表明，只有當煤體氧化產(chǎn)生的熱量大于其向環(huán)境散失的熱量時，熱量的積聚才能導致煤體自燃[4]。煤氧化作用產(chǎn)熱與環(huán)境散熱之間達到平衡狀態(tài)時的極限外界條件稱為煤自燃極限參數(shù)，包括最小浮煤厚度、下限氧體積分數(shù)和上限漏風強度[5]。由于煤礦開采環(huán)境的復雜多變性、煤層地質(zhì)賦存條件多樣性、采空區(qū)隨著工作面推進動態(tài)發(fā)展，致使煤自燃極限參數(shù)呈現(xiàn)出時變性、非線性、模糊性等特點，因此，采用科學準確的煤自燃極限參數(shù)測算方法對采空區(qū)煤自燃危險區(qū)域判定、高溫點范圍確定和預防煤自燃災害的發(fā)生有重要意義。

很多專家學者針對煤自燃極限參數(shù)測算進行了研究。徐精彩等[6]應用熱平衡法推導出煤自燃極限參數(shù)的計算方法，并應用BP 神經(jīng)網(wǎng)絡進行了煤自燃極限參數(shù)的預測[7]；孟倩等[8]提出了煤自燃極限參數(shù)的支持向量機預測模型; 周西華等[9]、郭亞軍等[10]以程序升溫實驗為基礎，研究了不同供風量條件下煤自燃極限參數(shù)變化；馬礪等[11-12]研究了氯鹽阻化劑、粒度對煤自燃極限參數(shù)的影響；易欣等[13]采用實驗方法研究了咪唑類離子液體處理對煤自燃極限參數(shù)的影響規(guī)律；Ma L 等[14]實驗分析了淮南礦區(qū)7 個煤礦煤樣的煤自燃極限參數(shù)；Wang C 等[15]、譚波等[16]研究了不同變質(zhì)程度煤的自燃極限參數(shù)變化規(guī)律及特征；Zhang Y 等[17]探討了煤自然發(fā)火升溫和降溫過程極限參數(shù)變化規(guī)律；章飛[18]分析了上、下分層采空區(qū)遺煤自燃極限參數(shù)并劃分了煤自燃危險區(qū)域；張辛亥等[19]采用分析解法計算煤自燃極限參數(shù)，并以凱達煤礦為例進行了驗證分析。

綜上所述，針對煤自燃極限參數(shù)的研究主要集中于實驗測試分析，由于影響煤自燃極限參數(shù)的因素十分復雜，對這些因素的影響程度以及影響規(guī)律等目前了解得還很少，同時，煤自燃極限參數(shù)影響因素在時空上呈現(xiàn)動態(tài)變化，數(shù)學上表現(xiàn)為復雜非線性，在實際計算煤自燃極限參數(shù)時通常簡化了條件，因而計算結(jié)果與現(xiàn)場實際還有一定的差距。為此，在對煤自燃極限參數(shù)影響因素分析的基礎上，引入灰色綜合關聯(lián)分析方法對煤自燃極限參數(shù)及其影響因素之間的關系進行定量分析，找出影響煤自燃極限參數(shù)的主要因素；提出一種改進的粒子群優(yōu)化算法（IPSO）對支持向量回歸（SVR）參數(shù)進行優(yōu)化選擇，建立參數(shù)優(yōu)化的支持向量機回歸模型來預測煤自燃極限參數(shù)，從而得到更好的逼近效果，提高預測精度。

1 采空區(qū)煤自燃及其影響因素

煤自燃是煤氧化產(chǎn)熱和周圍環(huán)境散熱共同作用的結(jié)果，是一個復雜的動態(tài)熱量變化過程。煤自燃的發(fā)生受各種復雜因素的影響，是多因素綜合作用的結(jié)果，除了與自身的氧化放熱性相關外，還與所處環(huán)境條件有關，采空區(qū)煤自燃影響因素如圖1。

由圖1 可以看出，影響采空區(qū)煤自燃的因素十分復雜，從而也決定了影響煤自燃極限參數(shù)因素的復雜性，因此，對煤自燃極限參數(shù)及其影響因素進行灰色關聯(lián)分析，確定影響煤自燃極限參數(shù)的主要因素和條件，這對煤自然發(fā)火的防治十分必要。

圖1 采空區(qū)煤自燃影響因素Fig.1 Factors of coal spontaneous combustion in goaf

2 采空區(qū)煤自燃極限參數(shù)灰色綜合關聯(lián)分析

2.1 煤自燃極限參數(shù)模型

煤體氧化升溫蓄熱量大于散熱量是保障煤自燃發(fā)生的前提，采空區(qū)遺煤自燃必須有較高的氧氣體積分數(shù)，確保遺煤能夠發(fā)生氧化放熱反應；遺煤須達到一定厚度，使產(chǎn)生的熱量能夠積聚；且漏風強度不能太大，避免產(chǎn)生的熱量被風流帶走。這些外界條件的極限值就是煤自燃極限參數(shù)，其數(shù)學模型如下[5]：

上限漏風強度：

最小浮煤厚度：

下限氧氣體積分數(shù)：

由煤自燃影響因素分析及煤自燃極限參數(shù)數(shù)學模型可知，影響煤自燃極限參數(shù)的因素十分復雜，并且在工作面開采過程中，不同時空條件下的影響因素實時動態(tài)變化，為了避免在實際應用時的簡化計算帶來的誤差，十分必要采用機器學習方法挖掘各影響因素與極限參數(shù)之間的非線性關系，建立煤自燃極限參數(shù)預測模型，以更好指導現(xiàn)場煤自燃防控。

2.2 灰色綜合關聯(lián)分析

在客觀世界中，有許多因素之間的關系是灰色的，難易辨別各因素之間關系的密切程度，這樣就難以找到主要矛盾和主要特性[20]。灰色因素關聯(lián)分析，目的是定量地表征諸因素之間的關聯(lián)程度，從而揭示灰色系統(tǒng)的主要特性。

以文獻［7］中表1 和表2 中的25 組數(shù)據(jù)為基礎，以上限漏風強度為系統(tǒng)特征行為序列X0=（x0（1），x0（2），…，x0（n）），以采空區(qū)距工作面距離（X1）、氧氣體積分數(shù)（X2）、煤溫（X3）、放熱強度（X4）、浮煤厚度（X5）作為相關因素序列Xi=（xi（1），xi（2），…，xi（n））。X0與Xi之間的相關關系即上限漏風強度與影響因素關系如圖2。

圖2 上限漏風強度與影響因素關系Fig.2 Relationship between the superior intensity limit of air leakage and factors

由圖2 可知，上限漏風強度X0與其影響因素Xi之間的線性關系非常弱，無法通過線性分析確定X0與Xi之間的緊密程度。灰色綜合關聯(lián)度分析注重事物發(fā)展過程的相似性與相近性，更加客觀地定量表達事物之間的關聯(lián)程度。

灰色綜合關聯(lián)度計算公式為：

式中：ρ0i為X0與Xi的灰色綜合關聯(lián)度，ρ0i∈［0，1］；θ 為灰色綜合關聯(lián)度系數(shù)，θ∈［0，1］，一般地，可取θ=0.5；ε0i為X0與Xi的灰色絕對關聯(lián)度，ε0i∈［0，1］；r0i為X0與Xi的 灰 色 相 對 關 聯(lián) 度，r0i∈［0，1］。

灰色絕對關聯(lián)度計算公式為：

其中：

灰色相對關聯(lián)度計算公式為：

其中，

灰色綜合關聯(lián)度既體現(xiàn)了折線X0與Xi的相似程度，又反映出X0與Xi相對于始點的變化率的接近程度，是較為全面地表征序列之間聯(lián)系是否精密的1 個定量指標[21]。根據(jù)以上理論計算公式，取θ=0.5，得出的各因素序列與上限漏風強度的灰色綜合關聯(lián)度見表1。

由表1 可知，浮煤厚度X5與上限漏風強度的灰色綜合關聯(lián)度最大，氧氣體積分數(shù)X2與煤溫X3次之，放熱強度X4與采空區(qū)距工作面距離X1最小。說明浮煤厚度對采空區(qū)煤自燃的上限漏風強度影響最大，這與預防采空區(qū)煤自然發(fā)火技術手段之一的提高工作面回采率，減少采空區(qū)遺煤與丟煤相吻合。浮煤越薄，越不利于熱量的蓄積，煤自熱環(huán)境更容易被破壞，無法達到煤自然發(fā)火的極限條件，從而煤自燃就不會發(fā)生。

表1 各因素與上限漏風強度的灰色綜合關聯(lián)度Table 1 The grey synthetic degree between the superior intensity limit of air leakage and factors

3 采空區(qū)煤自燃極限參數(shù)預測與驗證

3.1 支持向量回歸機

支持向量機[22]（Support Vector Machine，SVM）是由Vapnik 領導的研究小組于1995 年首先提出的一種基于統(tǒng)計學習理論的VC 維理論和結(jié)構(gòu)風險化原理的機器學習方法，可用于模式分類和非線性回歸。SVM 方法是把樣本數(shù)據(jù)“升維”，即將樣本空間線性不可分數(shù)據(jù)映射到線性可分的高維空間中，解決樣本空間中的高維度非線性問題[23-24]。低維空間到高維空間映射關系如圖3。

圖3 低維空間到高維空間映射關系Fig.3 Mapping relationship from low-dimensional space to high-dimensional space

支持向量機方法應用于回歸問題，得到支持向量回歸機（Support Vector Regression，SVR），在支持向量機回歸中，輸入樣本X 通過非線性映射Φ（X）映射到1 個高維的特征空間，然后在這個特征空間中構(gòu)造優(yōu)化的線性回歸函數(shù)。由于無論樣本大小、維數(shù)高低，徑向基核函數(shù)均適用，且收斂域?qū)?，只?個待優(yōu)化參數(shù)γ，易于實現(xiàn)，為此選用徑向基核函數(shù)。對于徑向基核函數(shù)的SVR，懲罰因子C 和核參數(shù)γ 是影響SVR 性能的主要參數(shù)，SVR 預測結(jié)果的準確性和可靠性依賴于（C，γ）的最佳選擇[25-27]。為此利用改進的PSO 對SVR 的懲罰因子C 和γ 核參數(shù)進行優(yōu)化選擇，來尋求最優(yōu)的平衡參數(shù)對（C，γ）。

3.2 粒子群優(yōu)化算法改進策略

粒子群優(yōu)化源于對自然界中鳥類捕食行為的研究，并應用于求解優(yōu)化問題。在每一次搜索迭代過程中，粒子通過種群個體極值與全局極值更新速度與位置：

式中：Vi、Li分別為第i 個粒子在D 維搜索空間的速度與位置；k 為當前迭代次數(shù)；ω 為慣性權重；c1、c2為非負常數(shù)，稱為加速度因子；r1、r2為分布于（0，1）之間的隨機數(shù)；Pi、Pg分別為個體極值與全局極值。

粒子群優(yōu)化算法雖然收斂速度快，但在算法的收斂后期，隨著種群多樣性的喪失，局部搜索能力不強，易陷入局部最優(yōu)，為此分別從以下3 個方面進行改進：

1）自適應慣性權重ω。對于適應度函數(shù)值優(yōu)于平均適應度函數(shù)值的粒子，其對應的慣性權重較小，使得該粒子局部搜索能力增強，反之對于適應度函數(shù)值差于平均適應度函數(shù)值的粒子，其對應的慣性權重較大，使得該粒子向全局較好的搜索區(qū)域靠攏。因此，根據(jù)個體適應度值的不同，慣性權重自適應調(diào)整，保證了粒子慣性權重的多樣性。

式中：ωmax、ωmin分別為慣性權重最大值0.9、最小值0.4；f 為粒子當前適應度函數(shù)值；fmin、favg分別為粒子適應度函數(shù)值的最小值、平均值。

2）替代機制。每當全局極值更新之后，將所有粒子按適應適度值大小排序，將最差的一半粒子的位置和速度用最好的一半粒子替代，以增強種群的多樣性和粒子的局部搜索能力。

3）網(wǎng)格搜索。為了進一步增強PSO 算法的局部搜索能力，提高算法的收斂精度，在PSO 算法中引入網(wǎng)格搜索。由于網(wǎng)格搜索是一種窮舉法，如果將網(wǎng)格搜索置于PSO 內(nèi)循環(huán)，勢必造成尋優(yōu)時間過長[28]。為此，在利用PSO 對SVR 的參數(shù)（C，γ）進行優(yōu)化后，再以PSO 尋優(yōu)的參數(shù)（Cp，γp）為中心，在C 維和γ 維上分別以步長cstep 和rstep 布下網(wǎng)格，然后計算二維網(wǎng)格［0.5Cp，1.5Cp］×［0.5γp，1.5γp］區(qū)域內(nèi)各交點的適應度函數(shù)值并與全局極值比較，若優(yōu)于全局極值，則用該點替換全局最優(yōu)粒子，否則保持全局最優(yōu)粒子不變。

步長取值方法：

式中：Cmax、γmax分別為參數(shù)C 與γ 取值范圍中的最大值。

網(wǎng)格搜索示意圖如圖4。

圖4 網(wǎng)格搜索示意圖Fig.4 Schematic diagram of grid search

3.3 基于IPSO-SVR 的煤自燃極限參數(shù)預測模型

3.3.1 樣本數(shù)據(jù)

以煤自燃極限參數(shù)中的上限漏風強度預測為例，選取影響上限漏風強度的主要因素即煤的放熱強度、煤體溫度、實測氧氣體積分數(shù)、采空區(qū)距工作面的距離和浮煤厚度作為輸入，上限漏風強度作為輸出，采用文獻［7］表1 中的20 組數(shù)據(jù)作為訓練樣本，表2 中的5 組數(shù)據(jù)作為測試樣本。為了取消各維數(shù)據(jù)間數(shù)據(jù)級差別，避免因為輸入輸出數(shù)據(jù)數(shù)量級差別較大而造成預測結(jié)果誤差較大，對輸入輸出數(shù)據(jù)進行歸一化處理，歸一化公式為：

3.3.2 IPSO-SVR 算法流程

根據(jù)SVR 原理和PSO 算法理論，IPSO 優(yōu)化SVR 參數(shù)建立煤自燃極限參數(shù)預測模型，IPSO-SVR建模流程如圖5。

圖5 IPSO-SVR 建模流程Fig.5 Flowchart of IPSO-SVR modelling

IPSO-SVR 算法具體計算步驟如下：

步驟1：訓練樣本和測試樣本的載入及其歸一化。

步驟2：粒子群優(yōu)化算法參數(shù)初始化。其中，c1＝1.5，c2＝1.7，C∈［0.01，50］，γ∈［0.01，50］，迭代次數(shù)為200，種群規(guī)模為20，采用5 折交叉驗證。

步驟3：產(chǎn)生初始粒子和速度，以交叉驗證均方誤差作為適應度函數(shù)值。依據(jù)適應度值，初始化粒子全局極值與個體極值。

步驟4：依據(jù)式（8）計算慣性權重，根據(jù)式（7）更新粒子的速度和位置。計算適應度值，更新個體極值與全局極值。

步驟5：執(zhí)行替代機制。

步驟6：判斷是否滿足終止條件。若滿足，則執(zhí)行步驟7，否則，返回步驟4。

步驟7：執(zhí)行網(wǎng)格搜索，輸出最優(yōu)參數(shù)對（C，γ）。

步驟8：將最優(yōu)參數(shù)賦給SVR，構(gòu)建參數(shù)優(yōu)化的SVR 預測模型。

3.4 預測結(jié)果

3.4.1 性能指標

為了進一步研究模型的預測性能，在保持訓練樣本和測試樣本不變的條件下，同時利用標準PSOSVR 模型（SPSO-SVR）、多元線性回歸模型（MLR）對煤自燃極限參數(shù)中的上限漏風強度進行訓練建模和測試分析，并與文獻［7］、文獻［8］的預測結(jié)果進行對比。上述模型預測性能評估指標如下：

平均絕對百分誤差（MAPE）：

絕對百分誤差在3%以內(nèi)的比例（MP）：

決定系數(shù)（R2）：

式中：fi為預測值；yi為真實值；n0為絕對百分誤差在3%以內(nèi)的數(shù)量；yˉ為真實值的平均值。

3.4.2 結(jié)果對比驗證

根據(jù)算法設計，經(jīng)過改進的PSO 對SVR 參數(shù)優(yōu)化選擇，尋優(yōu)結(jié)果為C=7.188 4，γ=0.617 8。同時，利用標準PSO 優(yōu)化SVR，參數(shù)尋優(yōu)結(jié)果為C=35.763 2，γ=0.510 1。將尋優(yōu)參數(shù)賦于SVR，建立參數(shù)優(yōu)化的SVR 預測模型，對測試樣本進行預測，不同模型預測結(jié)果對比結(jié)果見表2。

表2 不同模型預測結(jié)果對比Table 2 Comparison of prediction results for different models

為了更形象地分析數(shù)據(jù)，給出幾種模型的預測結(jié)果，預測結(jié)果散點圖如圖6。

圖6 預測結(jié)果散點圖Fig.6 Scatter plot of the predicted results

由圖6 可以看出，IPSO-SVR 預測結(jié)果基本與零誤差線重合，說明其預測結(jié)果與真實值相吻合，具備很好的擬合效果，而其他幾種模型相對于離散分布于零誤差線周圍，尤其是MLR 模型，預測結(jié)果距離零誤差線較遠，說明其預測結(jié)果產(chǎn)生了較大誤差。不同模型的預測結(jié)果相對誤差對比如圖7。

圖7 不同模型的預測結(jié)果相對誤差對比Fig.7 Comparison of relative error for different models

由圖7 可以看出，IPSO-SVR 模型相對誤差更小，在2.6%以內(nèi)，具有更高的預測精度，MLR 模型與文獻［7］神經(jīng)網(wǎng)絡模型預測結(jié)果明顯波動較大，預測結(jié)果與真實值之間產(chǎn)生了較大偏差。

不同模型性能指標匯總見表3，不同模型性能對比如圖8。

由表3 和圖8 可以看出，相比于其他模型，MLR預測結(jié)果最差，最大相對誤差達到19.788%，MAPE達到9.064 7%，說明線性回歸模型用于煤自燃極限參數(shù)的預測并不理想；從而也凸顯出煤自燃極限參數(shù)與其影響因素之間存在更明顯的非線性關系，而線性MLR 對非線性關系解釋能力有限，不宜用于煤自燃極限參數(shù)預測。

表3 不同模型性能指標匯總Table 3 Summary of performance indicators for different models

圖8 不同模型性能對比Fig.8 Performance comparison of different models

MAPE＜3%往往作為評估預測模型性能的標準，而工程應用中也以此作為衡量指標[29-30]?？傮w而言，支持向量機模型的MAPE 均小于3%，而用神經(jīng)網(wǎng)絡模型預測結(jié)果的MAPE 為6.15%，同時其MP 和R2均小于支持向量機模型。預測結(jié)果表明，與神經(jīng)網(wǎng)絡方法相比，支持向量回歸方法預測煤自燃極限參數(shù)精度更高、效果更好，更適合于解決小樣本、多維非線性映射問題。相比于采用變步長網(wǎng)格搜索方法優(yōu)化SVR 參數(shù)[8]，采用PSO 優(yōu)化SVR 具有更好的預測效果，而且，IPSO-SVR 預測結(jié)果的MAPE 減小達50%以上，說明智能優(yōu)化算法用于SVR 參數(shù)選擇具有提高其預測精度的潛力。同時，對標準PSO 進行改進，能夠顯著提高標準PSO-SVR 模型的非線性處理能力和預測精度。由此可以看出，提出的IPSOSVR 方法能夠更加精準地預測煤自燃極限參數(shù)。

4 結(jié) 語

1）煤自燃極限參數(shù)與其影響因素之間的線性關系較弱，為此引入灰色綜合關聯(lián)分析方法對煤自燃極限參數(shù)及其影響因素之間的關聯(lián)程度進行分析，浮煤厚度與上限漏風強度的灰色綜合關聯(lián)度最大，氧氣體積分數(shù)與煤溫次之，放熱強度與采空區(qū)距工作面距離最小。

2）基于PSO 原理，從迭代初期的自適應慣性權重調(diào)整、迭代中期的粒子替代機制和迭代結(jié)束后網(wǎng)格再搜索3 個方面對標準PSO 進行了改進，加快PSO 收斂速度，提高其局部搜索能力和穩(wěn)定性。提出了基于改進PSO 優(yōu)化SVR 參數(shù)的煤自燃極限參數(shù)預測模型，并對上限漏風強度進行了預測分析。

3）通過IPSO-SVR、SPSO-SVR、變步長網(wǎng)格搜索優(yōu)化SVR 參數(shù)方法、BP 神經(jīng)網(wǎng)絡和MLR 模型對比分析，IPSO-SVR 可獲得較高的預測精度，預測結(jié)果的平均相對誤差僅為1.140 1%，明顯優(yōu)于其他4種模型。IPSO-SVR 模型具有更精準的預測能力，更好的預測效果，更適合于影響因素復雜的煤自燃極限參數(shù)預測。