粒子群優化混合核極限學習機的構造煤厚度預測方法

范 君，王 新，徐 慧

(中國礦業大學計算機科學與技術學院，江蘇徐州221116)

(*通信作者電子郵箱344298099@qq.com)

0 引言

構造煤是指原生結構煤在構造應力作用下，煤層物理結構甚至化學成分發生明顯變化的煤體。早在19世紀20年代就有學者研究構造煤的產生和發育，至今已有研究表明煤礦瓦斯突出區域與構造煤存在著必然的聯系［1－2］。變形強的構造煤儲層的煤層結構及其厚度變化是影響煤與瓦斯突出的主要因素，在受到強烈流變的構造煤煤層中，煤層厚度值越大，煤層含氣量越大。理論上，煤層越破碎，煤粒表面能越增加，隨著煤層受到構造應力的增強，煤中大分子結構和空隙結構發生變化，從而改變了甲烷等氣體的吸附能力，就會導致瓦斯壓力增高，增大突出條件。事實也證明，煤體結構受到破壞的地區是瓦斯含量增高的地方，在高瓦斯礦區，順煤層滑動構造發育的區域即是瓦斯突出區域。我國幅員遼闊，地質活動頻繁，含煤地區大都經歷了復雜的地質構造演化，這為我國煤層構造煤提供了發育的環境，瓦斯爆炸事故時有發生。具不完全統計，自從1950年吉林省遼源礦務局富國二礦發生我國首次有記載的瓦斯突出事故以來，我國煤礦瓦斯突出事故時有發生，人員傷亡情況嚴重，經濟損失慘重。通過定量預測構造煤厚度來劃分瓦斯突出區域，將對煤礦安全管理和煤層氣的開發與利用起到至關重要的作用。傳統的構造煤厚度的預測方法有煤壁觀測和鉆孔取芯等，其中煤壁觀測法只適用于對已揭露區的構造煤進行編錄，而鉆孔取芯法由于其缺點明顯，例如取芯率低等，也無法保證煤礦的安全生產需求［3］。

當煤層含有構造煤時，其彈性性質明顯區別于正常煤層，含構造煤煤層的地震屬性也明顯區別于正常煤層的地震屬性，因此運用地震勘探技術解決煤礦構造煤的預測成為目前解決煤礦地質問題的新方法。地震屬性與煤層構造煤厚度的關系具有明顯的非線性特征，而機器學習算法很適合處理這種非線性特征的相關數據［4］。極限學習機(Extreme Learning Machine，ELM)是一種單隱含層前饋神經網絡，與其他神經網絡算法相比，ELM有明顯的優點:具有更快的學習速度;輸入參數更少;可以避免局部極值問題;具有較好的泛化性能。但是由于傳統的ELM算法本身存在一些問題:例如隱含層節點數量無法確定、易產生奇點問題等。針對這些問題，將核函數與傳統的極限學習機相結合，提出一種核極限學習機模型，直接采用核函數代替隱含層節點的顯式映射，無需給定隱含層節點數，計算速度快，泛化性能好，得到廣泛應用。

本文在以上研究的基礎上，將多項式核和高斯核相結合，構建混合核極限學習機(Hybrid Kernel Extreme Learning Machine，HKELM)預測模型，利用粒子群算法(Particle Swarm Optimization，PSO)較好的尋優能力，優化混合核極限學習機的核參數，并將預測模型應用至實際采區。在傳統的粒子群算法［5－6］中加入模擬退火的思想、隨進化次數逐漸減小的慣性權重以及基于反向學習變異操作，從而使得粒子群算法更加容易跳出局部最優，將地震屬性數據運用主成分分析技術(Principal Component Analysis，PCA)降維消除相關性，然后輸入構建的預測模型中，以提高采區構造煤厚度預測的精度和可靠性。

1 基本原理

1.1 粒子群算法基本原理

由于混合核參數較多，構造混合核極限學習機模型預測構造煤厚度關鍵的一步為優選核函數的參數。粒子群算法是一種群體智能的優化算法。

PSO算法首先在可解空間中初始化一群粒子，算法中每個粒子表示最優化問題的一個潛在解，利用位置、速度和適應度三種指標來表示一個已有粒子，其中適應度表示該粒子的優劣，其值由適應度函數計算得到［7］。粒子在解空間中運動，通過個體極值和群體極值更新個體的位置，其中個體極值是指個體所經歷位置中計算所得適應度值最優位置，群體極值指的是所有粒子適應度最優的位置。粒子每更新一次位置就計算一次適應度值，并選取最優適應度的粒子更新個體極值和群體極值。

假設在一個D維的搜索空間中有個粒子組成的種群X=(X1，X2，…，Xn)，其中第 i個粒子表示一個 D維向量 Xi=［Xi1，Xi2，…，XiD］T，代表尋優問題的一個潛在解。根據適應度函數初始化每一個粒子Xi對應的適應度值。設第i個粒子的速度為Vi= ［Vi1，Vi2，…，ViD］T，其個體極值為Pi= ［Pi1，Pi2，…，PiD］T，種群的全局極值為 Pg= ［Pg1，Pg2，…，PgD］T。粒子在每一次迭代過程中更新自身的速度和位置，更新公式如式(1)～(2):

其中:ω為慣性權重;d=1，2，…，D;i=1，2，…，n;k為當前迭代次數;Vid為當前粒子速度;c1和c2為非負的常數，稱為加速度因子;r1和r2為分布于［0，1］的隨機數。為了防止粒子的盲目搜索，通常將粒子的位置和速度限制在［－Xmax，Xmax］、［－ Vmax，Vmax］。

1.2 核函數

核函數方法的原理是通過一個特征映射可以將低維輸入空間中的線性不可分數據映射到高維特征空間中的線性可分數據，其本質是對應于高維空間中的內積，從而與生成的高維空間的特征映射一一對應。這種對應關系不僅使待解決的問題變得線性可分，同時又避免了高維空間可能帶來的維數災難。

設x，z∈X，X屬于R(n)空間，非線性函數Φ實現從輸入空間X到特征空間F的映射，其中F屬于R(m)，n m。根據核函數理論，有:

其中:〈Φ(x)，Φ(z)〉為內積;K(x，z)為核函數。從式(3)可以看出，核函數將m維高維空間的內積運算轉化為n維低維輸入空間的核函數計算，從而解決了在高維特征空間中維數災難等問題。傳統構造核函數是以Mercer定理為基礎，滿足Mercer條件的函數均可以作為核函數。目前使用較多的核函數主要有四類:

1)高斯核函數:

2)多項式核函數:

3)感知器核函數:

4)線性核函數:

此外，將簡單核函數構造成復雜混合核函數仍然滿足Mercer定理對核函數的要求。

1.3 HKELM的基本原理

傳統的ELM是一個三層前饋神經網絡，在這個網絡中輸入權值和隱含層閾值是隨機給的，影響模型性能的輸入參數是隱含層節點數，而整個模型需要訓練確定的是隱含層權值，其計算過程是通過計算一個矩陣的偽逆來實現的［8－9］。

若采用(a，b)表示輸入權值和隱含層閾值，樣本的訓練集用(x，t)表示，隱含層映射函數為g(x)，輸出權重為β，隱含層節點數為L，則極限學習機模型可表示為式(8):

令隱含層輸出矩陣為H，則式(8)可轉化為O=‖Hβ－T‖。極限學習機不斷學習的過程中，誤差O不斷減小，當無誤差學習時Hβ=T。隱含層節點數和樣本數相等時，H為非奇異矩陣，β可由H的逆矩陣求出，但一般隱含層節點數要遠小于樣本數，H不存在逆矩陣。此時廣義逆矩陣可以用于求解奇異矩陣的逆:

其中H+表示隱含層輸出矩陣的廣義逆矩陣［10］。

將核學習理論加入ELM模型中，隱含層將之前傳統激活函數顯式映射的任務交給核函數處理［11］。核函數的類型可以分為全局型核函數和局部型核函數，例如常用的高斯徑向基核函數為局部型核函數，在測試點附近數據點對核函數有影響，學習能力較強但泛化能力較弱;感知器核和多項式核為全局型核函數，遠離測試點的數據對核函數也有影響，泛化能力較強但學習能力較弱。

因此，本文將從兩類核函數中分別選取一種核函數組合成混合核，該混合核函數將作為混合核極限學習機模型的核函數［12－13］。

2 模型優化

2.1 粒子群算法的改進

針對HKELM模型參數較多的問題，本文提出用粒子群算法進行參數尋優［14］。為了提高粒子群算法尋優能力，針對粒子群算法容易陷入局部最小的問題，對傳統的粒子群算法主要作出三點改進:首先，利用模擬退火算法具有較強全局收斂性的特點，將模擬退火算法的思想引入粒子群算法，在加溫、等溫、冷卻的過程中增強粒子的多樣性，使得粒子群算法更有機會跳出局部最優的陷阱。其次，由于傳統的PSO算法中，慣性權重ω表示繼承先前粒子速度的能力，當ω較大時有利于算法的全局搜索能力，反之則有利于算法的局部搜索，ω值設置過大或過小容易造成PSO算法陷入局部最優或是無法真正地找出最優解。因此，加入隨迭代次數減小的慣性權重，慣性權重如式(10)所示:

其中:ωs為初始慣性權重;ωe為迭代至最大次數時的慣性權重;k為當前迭代代數;Tmax為最大迭代代數。ω隨迭代次數的增加而減小，既保證了迭代前期較強的全局搜索能力，又保證了在迭代后期粒子的局部搜索能力。最后，在粒子群算法中引入基于隨機的反向學習變異操作。變異概率如式(11):

其中:P 為變異概率，且當P ＞ 0.3時，P=0.5;Ps和Pe為常數，實驗中取 Ps=0.35，Pe=0.1。變異概率隨迭代次數先增大后減小，從而在一定程度上保留優秀粒子的同時增強粒子多樣性。變異公式為基于隨機的反向學習公式，如式(12)所示:

其中z為全局最佳粒子。

2.2 混合核極限學習機

極限學習機的核函數按照對數據的影響可以分為兩大類，局部型核函數和全局型核函數。局部型核函數學習能力強，但泛化性能相對弱;全局型核函數學習能力一般，但泛化性能強［15］。將局部核和全局核混合成為的混合核函數具備較強的學習能力，同時又有良好的泛化能力。

通常將局部的高斯徑向基核函數與全局的多項式核函數混合生成混合核構建的HKELM模型既具有較強的學習能力，又具有不錯的泛化性能，因此本文在此基礎上構建預測構造煤厚度的預測模型。高斯徑向基核函數具有較強的學習能力，能夠很容易地將訓練樣本在特征空間中線性分開，并能夠對相距一定范圍內的樣本準確預測，但對超過一定范圍的樣本則無法準確預測。

如圖1 所示，以X=0.4為測試點，分別取σ為0.1、0.2、0.3、0.4。徑向基核函數具有局部特性，體現在它只對測試點附近的值產生影響，且離測試點越近的點影響越大，由此可知，局部性核函數學習能力較強，泛化能力較弱。

圖1 徑向基核函數曲線Fig.1 Radial basis kernel function curve

多項式核函數具有較強的泛化能力，對即使相距很遠的樣本點也能產生影響，但對局部的樣本能力較弱［16］。

如圖2所示，以X=0.4為測試點，m=1，n=1，分別取d為1、2、3、4。多項式核函數的全局性特性體現在與測試點相距較遠的樣本點也能對函數值產生一定的影響，具有較強的泛化能力，學習能力較弱［16］。

圖2 多項式核函數曲線Fig.2 Polynomial kernel function curve

將多項式核函數與徑向基核函數結合，得到混合核函數:

如圖3 所示，以 X=0.4 為測試點，分別取 λ 為0.2、0.3、0.4、0.5，σ取值為0.1，d為2。由圖3可以看出，混合核函數不僅對測試點附近的樣本具有一定的影響，并且當樣本點遠離測試點時，仍能產生影響。因此，混合核函數結合了多項式核函數和高斯徑向基核函數的特點，有效地彌補了單核核函數的缺點，從而提高了擬合能力。

此外，為了防止模型復雜度過高引起過擬合現象，在核函數的基礎上加入L2正則項，有效地避免了噪聲和異常點對模型泛化性能的影響。

圖3 混合核函數曲線Fig.3 Hybrid kernel function curve

2.3 算法描述

改進PSO優化混合核極限學習機參數的算法具體過程如下:1)數據標準化，進行主成分分析處理，保存結果集備用;2)將σ、m、n、λ設為粒子(d=2)，隨機初始化粒子的位置和速度;

3)根據式(13)計算訓練集隱含層節點的輸出，并加入L2正則項;

4)計算隱含層輸出權值矩陣;

5)計算驗證集隱含層節點的輸出，并根據隱含層節點輸出計算得到測試集的預測值;

6)將均方誤差mse作為粒子適應度，計算個體極值和種群極值;

7)根據更新公式更新粒子位置和速度;

8)根據變異概率進行粒子變異;

9)計算更新后粒子的適應度，更新個體極值和種群極值，若未滿足最大迭代次數，則返回7)，否則繼續10);

10)保存種群最優適應度對應的粒子，即最佳混合核函數參數，將參數代入HKELM模型，獲得預測構造煤厚度的預測模型，并用測試集進行測試。

由于混合核極限學習機參數較多，而用傳統的粒子群算法優化混合核極限學習機無法找出最優參數下的預測模型，加入隨迭代次數逐漸減小的慣性權重、模擬退火思想以及基于反向學習的變異公式使得PSO在尋優能力和收斂速度上都有了改善。相比傳統ELM激活函數的顯示映射，混合核函數在較優的參數下，兼顧了數據的局部性和全局性，既具有較強的學習能力，又具有較好的泛化性能。因此，通過改進PSO尋找到最優的核函數參數，獲得最優核參數下的HKELM模型，理論上可以作為預測構造煤厚度的預測模型。

2.4 預測模型建立

2.4.1 測試數據

為了建立適合構造煤厚度預測的HKELM模型，首先建立了構造煤模型如圖4所示，煤層的頂底板參數為最薄處為0 m，最厚處為10 m;原生煤最薄處為0 m，最厚處也為10 m;直接頂/底是厚度為2 m的泥巖，老頂/底為砂巖。根據褶積模型，利用主頻為50 Hz的Ricker子波進行正演模擬，獲得其對應的正演地震剖面。利用圖4所示的構造煤模型正演地震剖面如圖5所示，分別提取出剖面中的正負相位的地震屬性。

利用圖4、圖5建立的模型，以所提取出的負相位地震屬性為例，分別進行有噪聲與無噪聲數據進行測試。地震屬性維數較大，為去除屬性中的相關屬性，實驗使用PCA對地震屬性進行降維去噪，取累計貢獻率大于95%的主成分并求其主成分得分作為實驗數據進行后續的實驗分析。

圖4 構造煤模型Fig.4 Tectonic coal model

2.4.2 參數設置

在運用改進的粒子群算法優化混合核極限學習機構建的構造煤厚度預測模型中，學習因子c1、c2均設置為2，最大迭代數MAXGEN=300，種群粒子數SIZEPOP=20，模擬退火算法初始溫度T=100，退火常數L=0.25，初始慣性權重ωs為1.2，迭代至最大次數時的慣性權重ωe為0.4，核極限學習機中核函數的參數由改進粒子群算法尋優得到。此外，為了探究本文提出改進的模型預測性能，將改進后的模型與傳統的BP神經網絡和用K折交叉驗證方法優化參數的支持向量機(Support Vector Machines，SVM)進行預測結果對比。其中BP神經網絡采用Matlab自帶函數newff()建立模型，輸入層神經元為9，輸出層神經元設為1，隱含層神經元數設為10，輸入層和隱含層、隱含層和輸出層之間的傳遞函數分別選擇tansig和logsig函數，最大訓練次數設為500，學習速率為0.1，訓練算法設置為 LM算法，動量因子為0.9，期望誤差設置為10－3，性能函數采用均方誤差函數(Mean Squared Error，MSE)和決定系數(Coefficient of Correlation)記為R2。利用Libsvm工具箱建立SVM預測模型，交叉驗證數的參數v設置為5，即參數c和參數g由5折交叉驗證得出，其中搜索空間為［－30，30］，步長設為0.5，核函數類型參數t設置為2，即采用徑向基函數作為核函數，支持向量機的類型參數s設置為3，即使用εSVR類型，即ε-支持向量機回歸算法，εSVR中損失值 ε 設置為0.1。

2.4.3 測試結果

無噪聲數據提取出26個地震屬性，通過PCA屬性降維獲得7個線性不相關屬性。含噪聲數據獲得11個地震屬性，進行屬性降維獲得9個線性不相關屬性。兩種數據集分別取98個樣本進行實驗并隨機抽取其70%作為訓練集，剩余30%作為測試集測試改進PSO-HKELM模型。預測結果如圖6所示。

圖6預測結果顯示，在無噪聲數據集下，改進的粒子群算法優化混合核極限學習機構建的構造煤厚度預測模型有很好的擬合效果，擬合度可以達到0.989，均方誤差為0.0041971;而含噪聲數據集的預測結果比無噪聲數據集的預測結果要差，擬合度為0.9293，均方誤差為0.0195，但在噪聲數據的影響下，同樣體現了較好的擬合效果。

圖6 模擬不同數據預測結果Fig.6 Prediction results of simulating different data

目前國內外已有利用智能算法預測構造煤厚度的研究［17］，為進一步驗證模型的優劣，在當前學者研究的基礎上，將改進的粒子群算法優化混合核極限學習機構建的構造煤厚度預測模型分別與用K折交叉驗證方法優化參數的SVM模型以及BP神經網絡進行實驗比較［18］。取98個含噪聲數據集樣本進行實驗構建模型，每個模型分別運行10次，且每次均隨機抽取其70%作為訓練集，剩余30%作為測試集，預測結果如圖7所示。

圖7 不同算法預測結果比較Fig.7 Comparison of prediction results of different algorithms

由圖7可以看出，采用的改進PSO-HKELM模型的預測準確率明顯高于另兩個模型預測準確率。在10次實驗中，改進PSO-HKELM模型的平均預測準確度達到90%以上。此外，將運行時間作為判斷預測模型時間復雜度的標準，模擬含噪聲數據在相同環境和參數下運行5次，運行時間如表1所示(保留兩位小數)。

表1 模擬含噪聲數據的運行時間 sTab.1 Running time of simulated data with noise s

由表1可以發現，改進模型運行時間較長，算法復雜度較高，但考慮到其預測精度提高較為明顯，預測結果較為穩定，且在實際開采中，預測構造煤厚度對實時要求不高，因此，本文提出的一種改進PSO-HKELM預測模型在預測構造煤厚度上的運用是可行的，且預測效果較好。

3 實例預測

3.1 研究區概況

本文選取陽煤集團新景煤礦蘆南二采區中部的15#煤層作為研究區域，結合實際地震屬性，對構造煤厚度進行預測。陽煤集團新景煤礦底層層滑構造特征明顯，瓦斯含量較高，極易發生煤與瓦斯突出，因此陽煤集團新景煤礦成為研究煤與瓦斯突出課題的理想區域。新景煤礦二采區面積為3.42 km2，其中15#煤層埋深 575.9 m ～640.93 m，平均深埋623.8 m，煤層厚度為6.25 m ～10.4 m，平均厚度 7.925 m;構造煤厚度為0 m～4.3 m，平均厚度為1.57 m。鉆孔共有10口，詳細情況見表2。

表2 15#煤層鉆孔揭露煤層信息 mTab.2 Exposing coal information of drills in 15#coal seam m

根據鉆孔揭露的煤層厚度和構造煤厚度，插值生成煤層厚度分布圖，如圖8所示。

3.2 構造煤厚度的預測

利用第2章建立的預測模型，對研究目標區域的構造煤厚度進行預測。如果只采取已知的10口鉆孔數據作為訓練樣本，則訓練樣本太少，預測模型的可靠性難以保證。因此，將鉆孔附近較小范圍內(15 m×15 m)的地震屬性提取出作為訓練集(共2232個樣本)，以提高訓練集的數量。15#煤層收集地震屬性共14維，首先對14維數據進行主成分分析，降維后的數據變為6維。然后采用2232×6的數據訓練構造煤厚度預測模型。將BP、SVM和改進PSO-HKELM模型分別運用到整個15#煤層的構造煤厚度預測中，將三種方法預測出的10口鉆孔處的構造煤厚度提取出來，與實際的鉆孔數據作對比。10口鉆孔的值與預測誤差如表3所示。

圖8 15#煤層厚度分布Fig.8 15#coal seam thickness distribution

表3 不同方法鉆孔處預測值對比 mTab.3 Prediction value comparison of drilling holes of different methods m

由表3可以得出，其中改進PSO-HKELM模型預測出的相對誤差最小，誤差均值為0.05(其中最大誤差為0.06，最小誤差為0.00)，而且預測結果均比較穩定;K折交叉驗證方法優化參數的SVM的預測精度其次，誤差均值為0.21(其中最大誤差為0.97，最小誤差為0.00);BP效果最弱，誤差均值為0.27(其中最大誤差為 1.02，最小誤差為 0.01)。因此，利用本文提出的改進PSO-HKELM模型在構造煤厚度的預測上效果更好。

數據經過PCA處理，運用改進PSO-HKELM模型預測出的構造煤厚度分布圖如圖9所示。

圖9 預測目標區構造煤厚度分布Fig.9 Distribution of tectonic coal thickness in prediction target area

根據圖9總體上可以看出，預測分布圖與煤層整體的構造煤分布具有很高的一致性。構造煤主要分布在井3-1381、井3-146和井3-137為中心的三個獨立區域。

4 結語

本文利用地震屬性預測構造煤厚度，首先運用主成分分析的方法對數據進行預處理，然后針對傳統預測方法的不足，建立了改進的PSO-HKELM的構造煤厚度預測的模型。采用較為常用的高斯核和多項式核組合成為混合核極限學習機的核函數，并用改進的粒子群算法對混合核參數進行優化。在改進的粒子群算法中加入模擬退火的思想、隨迭代次數減小的慣性權重以及基于隨機反向學習的變異操作，能有效地彌補粒子群算法較易陷入局部最優的缺點。最后，將建立預測模型應用于陽煤集團新景煤礦#15采區的實際煤層。通過模擬數據和新景煤礦#15區煤層數據的實驗分析，采用改進的PSO-HKELM算法對構造煤厚度建立預測模型，其性能將優于傳統的BP神經網絡和用K折交叉驗證優化參數SVM預測模型，且其精度較高，為預測煤礦采區煤與瓦斯突出形勢提供了一種可行的新途徑。構建的混合核極限學習機的構造煤厚度預測模型可以較為準確地預測實際采區構造煤厚度，但預測結果受地震屬性信噪比的影響較為嚴重，當訓練數據中含有噪聲的樣本較多時，預測結果不理想，因此尋找一種降低數據信噪比的處理方法需要進一步的研究。