吳順川，張晨曦，成子橋

(1.北京科技大學 土木與資源工程學院，北京 100083; 2.昆明理工大學 國土資源工程學院，云南 昆明 650093; 3.中電建路橋集團有限公司,北京 100048)

巖爆是高地應力狀態(tài)下地下硐室開挖過程中，硬脆性圍巖因開挖卸荷導致洞壁切向應力超出圍巖強度極限，引發(fā)巖體彈性能急劇猛烈釋放而產(chǎn)生碎片剝落、巖塊動力彈射拋擲現(xiàn)象，并伴有不同程度的爆炸、撕裂聲，直接威脅現(xiàn)場作業(yè)人員和設(shè)備安全，影響工程進度[1-4]。

巖爆預測是巖爆災害的有力防控手段[5]。巖爆預測方法主要分為兩大類，即單因素預測方法以及多因素綜合預測方法。前者利用已建立的巖爆判據(jù)實現(xiàn)巖爆烈度預測[2,6-13]，判據(jù)不同，巖爆烈度判別結(jié)果稍有差別;后者利用非線性理論綜合多種巖爆判據(jù)、巖石力學參數(shù)等實現(xiàn)巖爆預測，應用的非線性理論主要分為兩類:數(shù)學方法和智能算法。前者主要結(jié)合模糊數(shù)學綜合評判法[6]、可拓學理論[14]、逼近理想解排序法(TOPSIS)[15-16]、灰色系統(tǒng)理論[17]、云模型理論[18-19]、距離判別法[2]、理想點法[5,20]、未確知測度理論[21]等;后者主要應用BP神經(jīng)網(wǎng)絡[22]、支持向量機[23]、人工神經(jīng)網(wǎng)絡[24]、廣義回歸神經(jīng)網(wǎng)絡[25]等。與單因素預測方法相比，多因素綜合預測方法考慮巖爆影響因素更為全面、判別結(jié)果較為準確，具有一定優(yōu)勢。此外，基于微震監(jiān)測技術(shù)的巖爆預測也正逐步推廣應用，該技術(shù)通過獲取巖體內(nèi)部產(chǎn)生微破裂的前兆信息，可實時反映出圍巖應力狀態(tài)及強度等物理力學參數(shù)的變化特征，實現(xiàn)巖爆預測[3,26-27]。

以上理論方法、技術(shù)從不同角度揭示了巖爆現(xiàn)象發(fā)生規(guī)律及指標特征，豐富了巖爆預測理論。由于巖爆影響因素較多，其孕育演化機制尚不明朗，需要充分利用指標信息、挖掘指標特征與巖爆烈度等級的內(nèi)在聯(lián)系，并結(jié)合多種預測理論方法、技術(shù)才有可能實現(xiàn)巖爆的準確預測，故仍有必要引入新的預測模型開展巖爆烈度分級預測方法研究。

概率神經(jīng)網(wǎng)絡(PNN)是由Specht于1988年提出的，具有結(jié)構(gòu)簡單、訓練容易、收斂速度快，容錯能力強等優(yōu)點，并可用線性學習算法實現(xiàn)非線性學習算法的功能[28-29]，主要應用于故障診斷等模式分類問題[30-31]。PNN網(wǎng)絡的激活函數(shù)采用高斯函數(shù)時，指標變量需滿足互不相關(guān)且同分布[32]。事實上，大多數(shù)巖爆預測指標間均存在一定的相關(guān)性，故應用PNN網(wǎng)絡前需先消除巖爆預測指標間相關(guān)性。

消除指標間相關(guān)性的常用方法有限制指標數(shù)量、分離重迭元、修正指標權(quán)重、主成分分析法、因子分析法等[33]?？紤]到本文選用預測指標較多，選用主成分分析法(PCA)對巖爆預測指標數(shù)據(jù)預處理，除可消除指標間的相關(guān)性外，還能實現(xiàn)指標數(shù)據(jù)降維，提高PNN網(wǎng)絡訓練和分類速度。由此筆者基于主成分分析-概率神經(jīng)網(wǎng)絡模型開展巖爆烈度分級預測方法研究。

1 方法原理

1.1 主成分分析法

主成分分析法的核心是將原指標通過線性組合，使之成為一組新的、互相獨立的、并包含大部分原始信息的綜合指標，以此達到降維的目的[34]。在保證原始數(shù)據(jù)足夠信息量的前提下，對原多維指標變量數(shù)據(jù)進行降維，可壓縮指標個數(shù)、簡化數(shù)據(jù)，縮減巖爆烈度分級預測所消耗的時間。該法模型如下:

假定m個樣本和n個評價指標中第i個樣本的第j個指標取值為xij，構(gòu)造原始評價指標數(shù)據(jù)矩陣X=(xij)m×n：

(1)

則原變量指標為x1,x2,…,xn，設(shè)降維后綜合指標即主成分為y1,y2,…,yn，用X的n個指標向量作線性組合，即

(2)

式(1),(2)滿足:① 各等式系數(shù)平方和為1，如式(3)所示;② 任何兩個主成分相互獨立，互不相關(guān);③y1是x1,x2,…,xn的全部線性組合中方差最大者，y2是與y1不相關(guān)的x1,x2,…,xn的全部線性組合中方差最大者，yn是與y1,y2,…,yn-1均不相關(guān)的x1,x2,…,xn的全部線性組合中方差最大者[34]。

(3)

從相關(guān)矩陣出發(fā)求解主成分的步驟[35-36]。

步驟1:將原始數(shù)據(jù)標準化。

(4)

(5)

步驟2:計算指標間的Pearson相關(guān)系數(shù)矩陣，即

R=(rkl)n×n(k,l=1,2,…,n)

(6)

其中，rkl為第k個指標和第l個指標間的相關(guān)系數(shù)，且rkl=rlk(即對稱矩陣)，計算公式為

(7)

步驟3:計算相關(guān)矩陣R的特征值和特征向量。特征值記為λ1,λ2,…,λn，滿足λi≥0(i=1,2,…,n)，特征值對應的單位化特征向量記為p1,p2,…,pn。

步驟4:確定主成分個數(shù)。計算主成分的累計貢獻率，一般取特征值大于1、累計貢獻率達85%～95%所對應的前k個主成分。

(8)

式中，vs為第s個主成分的方差貢獻率。

(9)

式中，vsumk為前k個主成分的累計貢獻率。

步驟5:計算提取主成分的對應得分。主成分系數(shù)矩陣為:U=(p1,p2,…,pn)，若從原指標中提取了前k個主成分，則有

(10)

1.2 概率神經(jīng)網(wǎng)絡

概率神經(jīng)網(wǎng)絡以徑向基神經(jīng)網(wǎng)絡為基礎(chǔ)，采用統(tǒng)計學方法推導的激活函數(shù)替換S型激活函數(shù)，在某些易滿足的條件下由PNN實現(xiàn)的決策邊界漸近趨向貝葉斯最優(yōu)決策面[28]，是徑向基網(wǎng)絡與貝葉斯決策理論、概率密度函數(shù)的非參數(shù)估計的有機融合。PNN網(wǎng)絡由輸入層、模式層、累加層和輸出層組成，其網(wǎng)絡拓撲結(jié)構(gòu)如圖1[37]所示。輸入層接收訓練樣本數(shù)據(jù)，神經(jīng)元個數(shù)與輸入向量長度相等，輸入向量經(jīng)加權(quán)處理后傳遞給模式層見式(11):

Zi=xwi

(11)

式中，Zi為模式層第i個神經(jīng)元的輸入；x為輸入向量;wi為權(quán)值向量。

圖1 概率神經(jīng)網(wǎng)絡拓撲結(jié)構(gòu)Fig.1 Probabilistic neural network topology

模式層用來計算輸入向量與各模式的匹配關(guān)系，并返回一個標量值。該層神經(jīng)元個數(shù)與輸入的訓練樣本個數(shù)相同。向量x輸入到模式層，則模式層中第i類模式的第j個神經(jīng)元輸出關(guān)系為

φij=1/[(2π)0.5σd]e[-(x-xij)(x-xij)T/σ2]

(12)

式中，d為度量空間的維數(shù);σ為平滑因子;xij為第i類樣本的第j個中心;φij為模式層中第i類模式的第j個神經(jīng)元輸出值。

累加層的主要作用是線性求和、加權(quán)平均。將模式層中同一模式的神經(jīng)元的輸出做加權(quán)平均，該層神經(jīng)元數(shù)目與樣本模式總數(shù)相同，且與模式層神經(jīng)元建立連接關(guān)系的前提條件是兩者屬于同一模式分類。

(13)

式中，vi為第i類類別的輸出，(i=1，2，…，M)，M為訓練樣本模式總數(shù);L表示第i類神經(jīng)元的個數(shù)。

輸出層由競爭神經(jīng)元構(gòu)成，神經(jīng)元個數(shù)與訓練樣本模式總數(shù)相同，每個神經(jīng)元分別對應一種模式。用于為累加層的輸出做臨界值判別，將輸出層內(nèi)具有最大后驗概率密度的神經(jīng)元輸出為1，其余輸出為0。

2 巖爆烈度預測模型

2.1 PNN適用性分析

將巖爆烈度分級預測視為模式分類問題的依據(jù)(即PNN適用性)如下:① 巖爆烈度有分級。關(guān)于巖爆烈度分級國內(nèi)外目前雖尚無統(tǒng)一標準，目前常用的巖爆判據(jù)見表1。由表1可知，巖爆烈度常被分為4級，即巖爆烈度預測結(jié)果共有4類;② 巖爆預測問題有相應的指標量化數(shù)據(jù)，且一組樣本有惟一確定的巖爆烈度;③ 巖爆預測問題的基本思路主要是通過對巖爆指標變量數(shù)據(jù)的非線性處理，挖掘預測指標與巖爆烈度的內(nèi)在聯(lián)系，使模型判定結(jié)果能夠與實際巖爆等級一致，達到預測的目的。該過程與PNN模式分類問題相似，即先通過巖爆典型樣本訓練學習預測指標與巖爆等級間的對應關(guān)系，進而將訓練后的PNN網(wǎng)絡判定測試樣本巖爆等級。

表1 巖爆判據(jù)及烈度分級
Table 1 Rockburst criterion and intensity classification

類型巖爆判據(jù)公式巖爆等級無輕微中等強烈Russenes判據(jù)[6]σθ/σc<0.20.2～0.30.30～0.55>0.55Hoek判據(jù)[6]σθ/σc0.340.420.560.70強度應力比法Turchaninov判據(jù)[6]c=(σθ+σl)/σθ≤0.30.3～0.50.5～0.8>0.8陶振宇判據(jù)[7]σc/σ′1>14.514.5～5.55.5～2.5<2.5二郎山判據(jù)[8]σθ/σc<0.30.3～0.50.5～0.7>0.7N-Jhelum判據(jù)[9]σ′rm/σmax>74～72～41～2能量比[10]η=(φ1/φ0)×100%<3.53.5～4.24.2～4.74.7彈性能量指數(shù)[2]Wet=φsp/φst<2.0—2.0～5.0≥5.0能量法沖擊能量指數(shù)[10]Wcf=e1/e2<1.5—1.5～5.0≥5.0彈性應變能指標[10]Es=R2c/2E<0.20.2～0.50.50～0.75≥0.75峰值能量沖擊性指數(shù)[11]A′CF=Ue/Ua<22～5>5剩余彈性能指數(shù)[11]AEF=Ue-Ua<5050～150150～200>200脆性系數(shù)[12]B=σc/σt>4040.0～26.726.7～14.5<14.5巖性法能量儲耗指數(shù)[13]k=σcεf/σtεb<20—20～130>130變形脆性系數(shù)[12]Ku=u/ul<2.02.0～6.06.0～9.0≥9.0

注:σθ為硐室最大切向應力，MPa;σc為巖石單軸抗壓強度，MPa;σl為硐室軸向應力，MPa;σ′1為最大主應力，MPa;σ′rm為巖體強度，MPa;σmax為垂直于隧道軸向的最大水平主應力，MPa;η為能量比;φ1為破壞后巖石碎片彈射的動能，kJ/m3;φ0為巖石受載后儲存的最大彈性應變能，kJ/m3;Wet為彈性能量指數(shù);φsp為巖石試件在單軸壓縮過程中，達到σc以前(0.8σc～0.9σc)所存儲的彈性應變能，kJ/m3;φst為一個加載循環(huán)中產(chǎn)生微破裂和巖石塑性變形所耗損能量，kJ/m3;Wcf為沖擊能量指數(shù);e1為峰值前貯存的變形能;e2為破壞過程損耗的變形能;Es為彈性應變能指標;E為巖石的彈性模量，GPa;A′CF為峰值能量沖擊性指數(shù);AEF為剩余彈性能指數(shù);Ue,Ua分別為峰前彈性能密度和峰后破壞能密度，kJ/m3;B為脆性系數(shù);k為能量儲耗指數(shù);Ku為變形脆性系數(shù);σt為巖石單軸抗拉強度,MPa;εf,εb為峰值前和峰值后的總應變量;u,ul為巖石峰值荷載前的總變形和永久變形。

2.2 巖爆預測指標選取

巖爆的發(fā)生與巖體特征、地應力、地質(zhì)構(gòu)造和開挖擾動等因素密切相關(guān)[9]，具有顯著的隨機性、突發(fā)性和復雜性等特點[38-39]，給巖爆預測帶來了較大挑戰(zhàn)。巖爆預測的數(shù)據(jù)信息源于選取的指標，目前巖爆預測指標體系多數(shù)基于3～6種巖爆判據(jù)或巖石力學參數(shù)?？偟膩碚f，預測指標的選取應滿足:① 在室內(nèi)或現(xiàn)場易于獲取，具有較好的可操作性;② 具有較好的代表性;③ 能夠從不同角度反映巖爆特征信息。利用上述已有研究成果，根據(jù)巖爆的影響因素、特點及成因，選取圍巖最大切應力σθ、單軸抗壓強度σc、單軸抗拉強度σt、應力系數(shù)σθ/σc、脆性系數(shù)σc/σt和彈性能量指數(shù)Wet六種指標組成巖爆預測指標體系。前3種為圍巖應力參數(shù)，后3種為巖爆判據(jù)(巖爆烈度分級標準見表1)。

2.3 PCA-PNN巖爆烈度預測模型

主要過程為:① 根據(jù)選定指標搜集巖爆案例數(shù)據(jù);② 為消除指標間量綱、數(shù)量級差異，將原始數(shù)據(jù)標準化，再作相關(guān)性分析;③ 采用主成分分析法(PCA)對數(shù)據(jù)預處理，消除指標間相關(guān)性并降維，得到主成分數(shù)據(jù);④ 創(chuàng)建PNN網(wǎng)絡，設(shè)置不同的平滑因子;⑤ 將訓練樣本數(shù)據(jù)輸入PNN進行訓練，獲取各類樣本的特征信息;⑥ 訓練完成后，分別輸入訓練樣本和測試樣本對網(wǎng)絡性能測試，計算其分類正確率。

由表1可知，巖爆烈度分級有無、輕微、中等、強烈4類，筆者將該分級分別設(shè)定為1，2，3，4級，則PNN輸出向量為1×4行向量，類別編碼采用第i類在行向量第i個元素位置輸出為1、其余為0，如輸出向量為(0，1，0，0)表示該樣本屬于2級，即輕微巖爆。

3 預測模型實現(xiàn)過程

3.1 PCA預處理

依據(jù)所選指標，筆者搜集國內(nèi)外46組典型工程巖爆案例數(shù)據(jù)[12,40-42](表2)，前36組為訓練樣本，后10組作為測試樣本。

各指標間相關(guān)系數(shù)見表3，σc與σt,σt與σc/σt,σθ與σc,σt,σθ/σc的相關(guān)系數(shù)均超過0.5，預測指標間有較為明顯的相關(guān)性，巖爆烈度分級預測結(jié)果必定會受其影響，故有必要采用主成分分析法對原始數(shù)據(jù)作相關(guān)性消除處理。

表2 原始巖爆案例數(shù)據(jù)與主成分數(shù)據(jù)
Table 2 Data of original rockburst cases and principal components

注:37～46組為測試樣本,表7同。

表3 各指標間相關(guān)系數(shù)
Table 3 Correlation coefficients between different indexes

采用標準化的46組巖爆案例數(shù)據(jù)作主成分分析，其計算步驟如前文所述。各主成分方差貢獻率及累計貢獻率見表4，前3個主成分的累計貢獻率達86.718%，表明前3個主成分共包含86.718%的原始信息，且對應的特征值均大于1，滿足前述條件，因此選取這3個主成分。表5為主成分系數(shù)矩陣，由此得出主成分y1，y2，y3與6個指標變量之間的關(guān)系為

(14)

(15)

(16)

表4 PCA處理后的結(jié)果
Table 4 Results by PCA

成分特征值主成分貢獻率/%累計貢獻率/%12.76346.05246.05221.36722.79068.84231.07317.87686.71840.5689.45896.17750.2033.38199.55760.0270.443100

表5 主成分系數(shù)矩陣
Table 5 Coefficient matrix of principal components

指標主成分123σθ0.5000.3570.259σc0.425-0.4340.271σt0.521-0.146-0.298σθ/σc0.2510.7650.093σc/σt-0.3000.0080.804Wet0.381-0.2810.339

用式(1)，(2)將指標數(shù)據(jù)標準化后代入式(14)，(15)，(16)，可得主成分數(shù)據(jù)(表2)。3個主成分包含了原指標中的絕大部分信息，故可將3個主成分視為巖爆綜合預測指標RCI1,RCI2,RCI3。

3.2 PNN網(wǎng)絡參數(shù)確定

各組巖爆綜合預測指標RCI1,RCI2,RCI3構(gòu)成PNN網(wǎng)絡輸入向量。訓練樣本數(shù)據(jù)有36組，巖爆烈度有4種分類模式，因此PNN輸入層、模式層、累加層、輸出層的神經(jīng)元個數(shù)分別為3，36，4，4，其拓撲網(wǎng)絡結(jié)構(gòu)為3×36×4×4。先輸入36組訓練樣本對PNN進行訓練，再分別輸入訓練樣本、10組測試樣本對網(wǎng)絡性能進行測試。

平滑因子的選取對網(wǎng)絡性能起著關(guān)鍵作用。平滑因子取值過小易造成網(wǎng)絡過度擬合[43];平滑因子取值過大則不能完全區(qū)分細節(jié)，PNN網(wǎng)絡接近于線形分類器[44]。本文考慮平滑因子Spread值取值區(qū)間為[0.02，1.00]，滿足區(qū)間內(nèi)均勻分布，共50個Spread值。經(jīng)測試，巖爆訓練樣本、測試樣本的預測正確率隨Spread值的變化曲線如圖2所示。

圖2 模型預測正確率隨Spread值的變化Fig.2 Accuracy rate of the model prediction changes with Spread values

由圖2可知，在0.02～1.00，訓練樣本、測試樣本預測正確率總體上隨平滑因子的減小而趨增高;在某些小范圍內(nèi)，訓練樣本正確率、測試樣本正確率不隨平滑因子變化，說明PNN網(wǎng)絡對這些區(qū)間內(nèi)的spread值變化不敏感。訓練樣本、測試樣本正確率最大值及首次達到最大時的spread值見表6。

表6 模型最優(yōu)結(jié)果及對應參數(shù)
Table 6 Optimal results and corresponding parameters

項別最高正確率/%Spread值訓練樣本1000.36測試樣本900.52

3.3 PCA-PNN模型預測結(jié)果

取Spread值為0.36，46組巖爆案例判定結(jié)果如圖3,4所示。由圖3，4可知，訓練樣本測試誤判率為0，測試樣本中第8組巖爆案例誤判，預測結(jié)果比實際發(fā)生的巖爆等級略大一級，誤判率為10%。

為進一步說明模型的可行性，將預測結(jié)果與隨機森林(RF)模型、支持向量機(SVM)模型及人工神經(jīng)網(wǎng)絡(ANN)模型的兩組預測結(jié)果[40]對比，具體見表7。由表7可知，4種模型的預測結(jié)果相差不大，與實際巖爆烈度基本吻合。RF模型、SVM模型及ANN模型的訓練樣本誤判率、測試樣本誤判率見表8。由表8可知，PCA-PNN模型比SVM模型及ANN模型預測效果稍好，誤判率與RF模型訓練樣本平均誤判率、測試樣本平均誤判率一致。此外PNN收斂速度快，通常在數(shù)秒內(nèi)即可完成，且分類效果較為理想，表明PCA-PNN巖爆烈度分級預測模型是合理可行的。

圖3 訓練樣本預測結(jié)果Fig.3 Training sample prediction results

圖4 測試樣本預測結(jié)果Fig.4 Test sample prediction results

表7 與其它模型預測結(jié)果對比
Table 7 Contrast with other model predictions

序號模型輸出向量ⅠⅡⅢⅣ本文結(jié)果實際等級RF模型a組b組SVM模型a組b組ANN模型a組b組1001033333333200103333334430010333343444001033333333500103333333360010333333337001033333333800103333333390010333333331000103333333311001033333333120010333333331300014444444414100011111111151000111111111600103333333317001033333333180010333333331900103333333320010022222222210001444434442210001111111123001033333333240001444443442501002222332326010022222222270001444444442801002222332329000144444444

續(xù) 表

表8 3種模型誤判率匯總[40]
Table 8 Summary of misjudgment rates of three other models[40]

誤判率RF模型a組b組SVM模型a組b組ANN模型a組b組訓練樣本誤判率/%0011.18.35.611.1測試樣本誤判率/%02010102030訓練樣本平均誤判率/%09.708.35測試樣本平均誤判率/%101025

4 結(jié) 論

(1)選取圍巖最大切應力σθ、單軸抗壓強度σc、單軸抗拉強度σt、應力系數(shù)σθ/σc、脆性系數(shù)σc/σt和彈性能量指數(shù)Wet共6個指標構(gòu)成巖爆預測指標體系?？紤]到指標間的相關(guān)性及PNN網(wǎng)絡對指標變量間相互獨立的要求，采用主成分分析法對原始數(shù)據(jù)預處理，消除指標相關(guān)性并降維，得出3個線性無關(guān)的主成分即巖爆綜合預測指標RCI1,RCI2,RCI3。

(2)平滑因子σ(Spread值)對概率神經(jīng)網(wǎng)絡性能影響較大。在滿足均勻分布的前提下，選取[0.02，1.00]內(nèi)的50個Spread值，經(jīng)測試訓練樣本、測試樣本正確率最大值分別為100%，90%，兩者首次達到最大時對應的Spread值分別為0.36，0.52。

(3)為使訓練樣本、測試樣本正確率同時最優(yōu)，取Spread值為0.36。并將判別結(jié)果與RF模型、SVM模型及ANN模型的兩組預測結(jié)果比較，發(fā)現(xiàn)4種模型之間預測結(jié)果總體相差不大，但PCA-PNN模型比SVM模型、ANN模型的預測結(jié)果稍好，誤判率與RF模型的訓練樣本平均誤判率、測試樣本平均誤判率一致。此外PNN收斂速度快，通常在數(shù)秒內(nèi)即可完成，效果較為理想。以上表明PCA-PNN巖爆烈度分級預測模型是合理可行的。