李康楠,吳雅琴,杜 鋒,張 翔,王乙橋
(1.中國礦業大學(北京) 共伴生能源精準開采北京市重點實驗室,北京 100083;2.中國礦業大學(北京)應急管理與安全工程學院,北京 100083;3.中國礦業大學(北京) 機電與信息工程學院,北京 100083)
深部脆性巖體受外界擾動影響,發生動力失穩現象致使其內部彈性能極速釋放,從而導致巖體碎片彈射、拋擲與剝落,同時伴隨不同程度的爆炸聲響等現象,稱為巖爆地質災害[1-3]。巖爆災害發生的影響因素多,不確定性大,嚴重干擾現場作業的有序進行;巖爆災害危險性大,一旦發生中等、強烈巖爆,往往會造成一定的人員傷亡和巨大的經濟損失。隨著地下巖土工程逐步加深,巖爆災害事故日益增多,高效準確的巖爆烈度等級預測方法研究迫在眉睫[4-5]。
國內外學者對巖爆烈度等級預測方法的研究大致可分為判據預測方法和綜合分析預測方法2 種[6]。判據預測方法首先會確立巖爆烈度等級的單因素分類標準,然后將預測區域該因素的數值大小與預定標準進行對比,從而預測巖爆的烈度等級,如:Russenes[7]與Hoek[8]與二郎山[9]應力系數判據、Turchaninov[10]判據、Kidybinski[11]與Singh[12]巖爆傾向性指數判據、N-Jhelum[13]判據、Aubertin[14]改進脆性指數判據、Wang[15]最大儲存彈性應變能判據、陶振宇判據[16]等。巖爆誘發因素眾多,單因素很難全面揭示災害發生規律,且各因素之間可能存在相互耦合的現象,故僅以單因素判據預測巖爆災害缺乏可信度。因此,越來越多的學者采用多因素分析預測方法進行綜合研究,以期通過數學方法、機器學習與智能算法達到準確預測巖爆烈度等級的目的。傳統數學方法如貝葉斯模型[17]、功效系數法[18]、支持向量機[19]、決策樹[20]等雖應用廣泛,在一定程度上提升了巖爆烈度等級預測的準確率,但存在一些方面的缺陷:在確定指標權重與定性因素方面受人為主觀影響嚴重,很難客觀預測巖爆烈度等級;影響因素之間相互聯系、共同作用,預測過程往往忽略非線性因素。為彌補上述缺陷,國內外學者展開機器學習與智能算法領域的研究,如SVM (Support Vector Machines)模型[21]、IPP-PNN(Improved Projection Pursuit-Probabilistic Neural Networks)模 型[22]、PCA-PNN(Principle Component Analysis-PNN)模型[23]、PCA-RBF(PCA-Radial Basis Function)[24]模型等,進一步提高了巖爆烈度等級的預測精度,但存在兩方面的局限性:所需訓練樣本量與計算權值量大,導致模型復雜程度提高與運算時間延長;需要配合數據降維方法進行模型訓練,易丟失重要的數據信息從而影響模型預測的準確度,故需要繼續深入研究并探索新的預測方法。
深度學習算法卷積神經網絡(Convolutional Neural Network,CNN)具有局部感知與權重共享的優點,能夠極大減少連接權值的數量并降低模型復雜度,挖掘線性與非線性數據深層次的特征規律與數據之間的內在聯系;無需配合額外的分析法即可降低多重共線性并提高模型的泛化能力,在短期電力負荷預測[25-26]、短時交通流預測[27-28]、電廠存煤量預測[29]與煤層底板突水預測[30]等方面都具有強大的功能。考慮到數據缺失使訓練樣本減少從而導致模型準確率降低的問題,筆者選用鏈式方程多重插補(Multiple Imputation by Chained Equations,MICE)法對空缺數據進行填補[31-33]。該方法能夠通過表達數據的不確定性并進行綜合分析,得到更加可靠與精確的插補結果。基于MICE-CNN,選取120 組典型巖爆案例現場數據,通過參數優選、數據處理與網絡訓練建立模型,并與RBF、SVM、PNN 模型的預測結果進行對比,期望得到更加可靠有效并利于工程實際的巖爆烈度等級預測模型。
基于CNN 的巖爆烈度等級預測模型建立流程如圖1 所示,模型的建立主要分成4 個步驟。第一步,列舉并選取巖爆災害的影響因素,通過綜合分析建立契合本預測模型的指標體系。第二步,搜集原始數據,對數據進行基于拉伊達準則的異常值剔除與基于MICE的缺失值插補,得到完整數據后進行數據集分割與數據轉換,搭建CNN 的初始框架并進行訓練,從而優選模型超參數,建立基于CNN 的巖爆烈度等級預測模型。第三步,輸入現場數據進行驗證模型準確度預測。第四步,建立基于RBF、SVM 與PNN 的巖爆烈度等級預測對比模型,輸入相同數據進行預測,并與CNN模型比較預測結果,并觀察預測結果的混淆矩陣,對比誤判結果的傾向性,最終得到最優模型。

圖1 巖爆烈度等級預測模型建立流程Fig.1 Flow chart for establishing prediction model of rockburst intensity grade
本研究通過對國內外大量文獻及礦井現場資料的查閱收集[7-9,14-16],得到了主流的巖爆分級標準與對應計算公式,詳見表1。

表1 巖爆分級標準Table 1 Rockburst classification standard
巖爆災害的發生與圍巖特征、地質構造、地應力與外界擾動等因素密切相關[34-36],這導致了影響巖爆發生的因素眾多,使指標選取困難。目前國內外尚無統一的巖爆烈度等級預測指標體系,但多數包含3~6 種巖爆判據或巖石力學參數。過多影響因素組成的預測指標體系既不易于獲取,增加計算的冗余程度與預測時間,而且相關度小的因素數據離群點還會在一定程度上影響模型的學習與判斷;過少的指標缺乏代表性,難以支撐災害預測,從而導致預測結果出現較大偏差。由表1 的計算公式可以得到計算巖爆分類標準所常用的具有較好代表性的影響因素,在此基礎上,考慮現場獲取便捷,且能夠全面反映巖爆特征信息等條件,最終選取硐室最大切向應力σθ、巖石單軸抗壓強度σc、巖石單軸抗拉強度σt、巖體應力系數σθ/σc、巖石強度脆性系數σc/σt和彈性變形能系數Wet這6 種指標組成巖爆預測指標體系。
基于本研究所建立的巖爆預測指標體系,通過現場實測、查閱資料等方式搜集巖爆烈度等級工程實例,得到120 組實際發生烈度的原始數據用來建立巖爆預測模型,部分數據見表2,其中,無、弱、中等、強烈巖爆4 種烈度等級分別用數字1、2、3、4 表示。

表2 部分工程實例原始數據[7-24]Table 2 Partial engineeri ng examples raw data[7-24]
CNN 卷積神經網絡的結構如圖2 所示。

圖2 CNN 網絡結構Fig.2 Structure of CNN
1) 卷積層
卷積層線性抽取局部范圍內的神經元信息與特征,然后運用非線性激活函數對神經元進行激活。卷積運算主要為2 個數值矩陣的變換運算,過濾器(卷積核)首先會選擇每次移動的距離(步長),然后與輸入矩陣進行對應的映射點積運算,運算結果組成輸出矩陣(特征圖)。假定在l層的卷積層有n[l]個濾波器,權重為K[l]∈大小為3×3,其中n[l-1]是前一層的過濾器數量。這些過濾器以[1,1]的步伐遍歷整個輸入特征圖,卷積操作可用下式表達,輸出特征Y[l]是:
式中:f(·)為激活函數;b為偏執項。
2) 池化層
池化層可降低模型過擬合并提高模型的泛化能力。常用最大池化法進行降采樣,其表達式為:
式中:X為輸出;α為乘性偏置;S(x)為降采樣函數。
3) 全連接層 (Fully Connected Layers,FC)
全連接層能夠將最后一個池化層所得到的特征圖平鋪為一條一維向量,從而將數據進行由高到低的維度變換,同時不丟棄任何有用的信息,其中每個輸入通過一個可學習的權重連接到每個輸出。
4) 損失函數層
損失函數可以表達模型預測值與真實值之間的不一致程度,損失函數與模型的魯棒性是負相關,在模型的學習過程中起到一定的指導作用。本研究選用Softmax 損失函數[37]。
1) 基于拉伊達準則的異常值剔除
在數據測量與記錄的過程中存在人工操作不當、采動影響等產生的錯誤數據,剔除錯誤數據能夠消除其對模型預測準確性的影響。對120 組原始數據的前4 列數據進行處理,設數據集qi(i=1,2,···,m),計算剩余誤差vi,并計算標準偏差β,若某一值的剩余誤差va(1≤a≤m)滿足下式,則認為該值為異常值并予以剔除。
式中:qa為數據集中第a個數據;為數據集均值。
2) 基于MICE 的缺失值填補
對前4 列中經過異常值剔除的數據在原本空缺的地方進行數據填補。通過在每個空缺位置使用插補模型進行多次填補產生若干完整數據集,然后對若干完整數據集分別進行分析,得到相應的分析結果,對結果做出綜合推斷與比較分析,最終得到最優估計值。本研究使用R 語言中的mice 程序包實現MICE。
本研究選用隨機森林(Random Forest,RF),貝葉斯線性回歸(Bayesian Linear Regression,BLR),極 限樹(Extra Tree,ET),K 臨近(K-Nearest Neighbor,KNN)4 種插補模型作為MICE 的估計器,并結合中值法(Median)與均值法(Mean)2 種傳統插補方法分別進行數據填補,為驗證MICE 的優越性并優選最佳插補模型,對比6 種填補結果的均方根誤差(ERMS),如圖3 所示。

圖3 不同插補模型/方法的ERMSFig.3 ERMS for different imputation models or methods
觀察圖3 可以得出MICE 比傳統插補方法擁有更高的精度,其中ET 模型的精度最高,故本研究選用ET 模型作為估計器對前4 列進行MICE 插補,插補后計算并補全后2 列數據,完整數據如圖4 所示。

圖4 完整數據Fig.4 Complete data
3) 數據集分割
為使模型以有限的現場實測數據進行充分學習,將樣本數據10 等分,其中訓練集占8 份、驗證集與測試集分別為1 份[38],即120 組樣本數據中隨機抽取96 組作為CNN 模型的訓練集樣本,12 組驗證集樣本,共108 組數據為預測模型建立過程中的訓練樣本,剩余12 組測試集樣本作為工程實例以進一步驗證模型準確率;3 部分數據相互獨立且均具有代表性。
4) 數據轉換
CNN 在圖像分類預測方面具有良好的預測精度,故本研究需要將數值數據轉換為包含RGB 三通道的圖像數據[39]。由于本研究的訓練樣本僅包含6 個指標,將數據平鋪為6×1×1 的一維圖像數據,既可以充分提取樣本數據內涵特征信息,又可降低特征維度使神經網絡的訓練與計算更加便捷。
1) CNN 神經網絡構造
第一層為圖像輸入層(Image input)。第二層為卷積層1(Conv_1),由于數據為一維,所以卷積核大小選取為z×1,而z為1 時會喪失感受野的作用,故需z>1;當卷積核為偶數時,輸出特征圖比輸入特征圖尺寸變小,可學習的特征信息減少,故需z為奇數。考慮到卷積核過大可能忽略數據特征信息,故本研究所選用的卷積核大小為3×1;同理,本研究以1 為步長進行卷積操作能夠盡量減少信息的丟失;采用的硬件資源為單GPU,結合硬件配置選擇卷積核的數目為16;卷積層與池化層的特征邊緣處理方法有不填充與前后補零2 種,不填充即忽略最后未卷積的區域并降低特征圖尺寸,前后補零方式可以使最后輸出的特征圖尺寸和原圖尺寸一致。本研究選用前后補零的特征邊緣處理方法,可保證輸入與輸出特征圖尺寸Sc一致。輸出特征圖的尺寸計算公式如下:
式中:Co為輸出通道數;Ho、Io分別為輸出特征圖的高和寬;Hi、Ii分別為輸入特征圖的高和寬;F1與F2分別為卷積核的高與寬;M為步長;P1、P2分別為補零的行數和列數。經過計算,特征邊緣處理后的特征圖尺寸始終保持為6×1×1。
第三層為批量歸一化層(Batch norm)[40],該層對訓練樣本進行局部歸一化,加快模型的收斂與運算速度,并且使模型更加穩定。本研究對是否添加批量歸一化層進行對比訓練,結果顯示,添加該層后預測準確率略微提高,運算速度顯著提高。
第四層為激活函數層1(ReLU_1),激活函數通過添加非線性因素的方法捕捉非線性信息,進一步提高線性模型的表達能力,常用的激活函數有[41-42]:(1) Sigmoid 函數,該函數在其飽和時的梯度值相對低下,梯度的耗散問題將會隨著模型中神經網絡層數的增多而嚴重;(2) Tanh 函數,該函數的缺點是仍然具有飽和的問題;(3) ReLU 函數,可有效解決Sigmoid 函數所存在的梯度耗散問題,網絡稀疏性較大從而降低Tanh 函數存在的飽和問題,且運算效率高。3 個函數公式分別如下:
以訓練樣本的預測準確率為校驗目標對各激活函數進行模型訓練,結果如圖5 所示。觀察圖5,以ReLU 函數進行模型訓練的準確率均較高,特別地,訓練集準確率達到100%,故本研究選取ReLU 函數作為激活層函數。

圖5 激活函數的選取Fig.5 Activation function selection
第五層為池化層1(Maxpool_1),常用的池化方法有平均池化法與最大池化法。平均池化法能夠最大限度地保留數據背景信息,最大池化法能夠盡可能地提取數據的特征紋理,本研究期望能夠挖掘到更多的數據特征,故選取最大池化法進行池化操作。池化核與步長的確定過程類似于卷積層卷積核與步長的確定過程,確定本研究池化核大小為3×1,步長為1,特征邊緣處理方法選擇前后補零。
第六層為卷積層2(Conv_2),由大小為3×1 的卷積核以1 的步長進行卷積操作,生成32 幅特征圖,特征邊緣處理方法選擇前后補零。第七層為批量歸一化層。第八層為激活函數層2。第九層為池化層2(Maxpool_2),特征邊緣處理方法選擇前后補零。
第十層為全連接層(FC),將上一層數據展開,經過激活函數Softmax 激活后挖掘特征信息,得到4 類輸出,即巖爆烈度的4 個等級。本研究所建立的CNN神經網絡共有12 層,各層參數見表3。

表3 CNN 網絡結構參數Table 3 Structural parameters of CNN
2) 超參數優選
優化器的選擇方面,目前常用且性能優秀的優化器有3 種:SGDM、RMSProp 和Adam[43]。SGDM 在梯度下降的過程中引入一階動量解決了可能僅找到局部最優點而不是全局最優的缺陷;RMSProp 加入了迭代衰減,不會在迭代過程中梯度下降過大使自適應梯度出現變化異常等現象;Adam 結合了前2 種優化器的優點,在提高自適應學習率、尋找全局最優點等方面皆較為突出。以訓練樣本的預測準確率為校驗目標,對RMSProp、SGDM 與Adam 進行模型訓練,結果如圖6 所示,經對比,Adam 能夠在較少的訓練次數下達到較高的預測準確率,故本研究選用Adam 優化器進行模型優化。

圖6 優化器的選取Fig.6 Optimizer selection
由圖6 可知,訓練次數達到200 時,預測準確率便不再升高,繼續升高訓練次數只會延長模型訓練時間,故本研究的最大訓練次數選為200。預測模型的超參數詳見表4。

表4 預測模型超參數Table 4 Prediction model hyperparameter
將本研究構造的CNN 神經網絡與選取的超參數相結合建立巖爆烈度等級預測模型,輸入訓練樣本進行訓練,訓練集預測結果如圖7a 所示,驗證集預測結果如圖7b 所示。由兩圖可得,訓練集預測準確率達到100.00%;驗證集的預測準確率達到91.67%,其中第七組真實值為弱巖爆,而預測值為中等巖爆,其他組預測結果均正確,錯誤率在合理范圍內。由訓練樣本預測結果可得,本研究所建立的CNN 巖爆烈度等級預測模型是合理可行的。

圖7 訓練集與驗證集預測結果Fig.7 Prediction results of the training set and verification set
為進一步說明模型的可行性,將包含12 組數據的測試集分別輸入本次所建立的CNN 預測模型和3 個對比模型:RBF、SVM 與PNN 模型中進行預測,預測結果見表5。觀察各模型預測結果準確率可知,CNN模型的預測結果準確率最高,僅第9 組數據預測錯誤,其他組皆能夠正確預測巖爆烈度等級。

表5 各個模型的預測結果與準確率Table 5 Prediction results and accuracy of each model
CNN 模型測試集混淆矩陣[44]如圖8a 所示。圖中可以看到有一組現場案例的實際結果為中等巖爆,而模型預測結果誤判為強烈巖爆,預測結果略高一級,雖在應急防治方面會消耗更多的人力物力,但在誤判的情況下依然能夠保證現場工作人員的生命安全。

圖8 CNN 與PNN 測試集混淆矩陣Fig.8 Confusion matrix of CNN and PNN test set
PNN 模型測試集混淆矩陣如圖8b 所示。PNN 預測模型的準確率與CNN 模型相近,但觀察混淆矩陣,有2 組現場案例的實際結果為中等巖爆,而模型預測結果分別誤判為一組無巖爆與一組弱巖爆。當誤判為弱巖爆時,現場人員的準備與預防措施不夠充分,難以保障人身與財產的安全;而當誤判為無巖爆時,巖爆災害將對施工現場造成極大的破壞并嚴重危害工作人員的生命,造成巨大的財產損失。
通過以上分析可以得到,本研究所建立的基于CNN 的預測模型是可行的,且具有準確率高、安全系數高等優點。由于現場數據收集較為困難,缺乏地下水、微震、聲發射等影響巖爆的關鍵因素與動態指標;需要收集與應用的現場數據量需進一步擴大,未來需要結合更多案例與更多影響指標進行綜合研究。
a.建立巖爆烈度等級預測指標體系,對原始數據進行基于拉依達準則的異常值剔除與MICE 的缺失值填補,建立RF、BLR、ET、KNN 4 種插補模型并結合2 種傳統插補方法,依據ERMS評價指標優選MICE 的估計器,其中,ET 模型得到的觀測值與真值偏差最小,能夠保證模型預測的準確度。
b.搭建CNN 神經網絡初始框架,確定卷積核和池化核大小、特征邊緣處理方法與激活函數等;優化模型超參數,確定優化器函數等,提高了模型預測準確率,降低了運算時間,增強了模型的適用性。運用訓練樣本對模型進行訓練,訓練集準確率為100.00%,驗證集準確率為91.67%。
c.建立基于RBF、SVM 與PNN 對比模型,并對實例應用結果進行綜合分析顯示,基于CNN 的巖爆烈度等級預測模型準確率最高,為91.67%。對比CNN與PNN 模型的混淆矩陣得到,在誤判的情況下CNN預測模型安全性更高。
d.未來需進一步擴大用來建立指標體系與預測模型的現場數據量,結合現場案例的更多影響因素進行綜合分析。有必要運用先進算法進行超參數優化,建立成熟的深度學習模型對巖爆災害進行精確預測。