郭延華，趙 帥

(河北工程大學(xué) 土木工程學(xué)院，河北 邯鄲 056038)

巖爆是隧道開挖時，因巖體受到擾動，其中聚積的彈性能突然釋放，導(dǎo)致巖體脫落、爆裂、彈射的動力失穩(wěn)現(xiàn)象，嚴(yán)重時可造成重大財產(chǎn)損失和人員傷亡[1]。隨著我國經(jīng)濟的快速發(fā)展，對地下深部資源開采和地下工程建設(shè)的需求愈發(fā)強烈，對巖爆預(yù)測與防治的研究也愈發(fā)重要。

在巖爆烈度預(yù)測方面，國內(nèi)外學(xué)者進行了大量的研究[2-3]，神經(jīng)網(wǎng)絡(luò)[4]、功效系數(shù)法[5]、模糊數(shù)學(xué)[6]、灰色理論[7]、支持向量機[8]、數(shù)據(jù)挖掘[9]、PCA-PNN法[10]、PCA-OPF法[11]等數(shù)學(xué)模型與智能算法被提出。但上述預(yù)測方法一般基于少量巖爆案例的小樣本巖爆數(shù)據(jù)，而巖爆數(shù)據(jù)的數(shù)量與質(zhì)量往往影響巖爆預(yù)測的準(zhǔn)確性和可靠性[12]，所以由此建立的預(yù)測模型的可靠性與泛化性較差。此外上述方法的預(yù)設(shè)參數(shù)往往由人工經(jīng)驗設(shè)置，預(yù)測結(jié)果受人主觀因素影響較大，或僅對巖爆數(shù)據(jù)作線性降維，未考慮數(shù)據(jù)間的非線性關(guān)系，從而未充分保留巖爆數(shù)據(jù)的特征信息。因此在巖爆預(yù)測方面需要進一步研究，并提出經(jīng)濟可靠的巖爆預(yù)測方法。

核極限學(xué)習(xí)機(Kernel Based Extreme Learning Machine，KELM)是一種分類預(yù)測算法，該算法將核函數(shù)與極限學(xué)習(xí)機(Extreme Learning Machine，ELM)相結(jié)合，以提高ELM算法的預(yù)測性能。故本文巖爆烈度預(yù)測將基于KELM實現(xiàn)，同時為解決因人工設(shè)置參數(shù)而影響KELM預(yù)測效果的問題，使用鯨魚優(yōu)化算法(Whale Optimization Algorithm，WOA)尋優(yōu)KELM的懲罰系數(shù)C和核函數(shù)參數(shù)s以提高其預(yù)測性能。此外，考慮到巖爆評價指標(biāo)之間的多維交叉冗余性，使用核主成分分析(Kernel Principal Component Analysis，KPCA)對巖爆數(shù)據(jù)集做特征壓縮，相比于傳統(tǒng)的主成分分析(Principal Component Analysis，PCA)，KPCA能對數(shù)據(jù)集做非線性降維，從而能更加充分地保留特征信息，以提高KELM的訓(xùn)練速度和預(yù)測速度，故本文建立KPCA-WOA-KELM的巖爆烈度預(yù)測模型。

1 模型原理與建立

1.1 KPCA理論基礎(chǔ)

PCA是對數(shù)據(jù)做線性降維的傳統(tǒng)降維算法，而對具有非線性關(guān)系的數(shù)據(jù)做特征提取時，卻無法充分保留數(shù)據(jù)的特征信息。KPCA在其基礎(chǔ)上結(jié)合核函數(shù)，將數(shù)據(jù)映射至高維空間中做數(shù)據(jù)壓縮。KPCA能在充分保留數(shù)據(jù)特征信息的基礎(chǔ)上做數(shù)據(jù)壓縮[13]。

1.1.1 KPCA基礎(chǔ)理論

使用函數(shù)φ將樣本映射至高維特征空間RF，映射值為φ(x1)，φ(x2)，…，φ(xn)并使用PCA方法得到協(xié)方差矩陣為：

(1)

其特征方程為：

Cν=λν

(2)

式中，λ為協(xié)方差矩陣的特征值，ν為特征向量，由式(1)、式(2)可得ν：

(3)

k(xi,xj)=φ(xi)Tφ(xj)

(4)

對于式(2)，任意的k=1，2，…，n，有：

φ(xk)Cν=λφ(xk)ν

(5)

分別將式(1)、式(3)與式(4)帶入式(5)可得：

Kα=λnα

(6)

式中，K為k對應(yīng)的核矩陣，K=k(xi,xj)，α=(α1,α2,...,αn)。

通過式(6)，求得特征值λ1≥λ2≥…λn及其對應(yīng)特征向量α1,α2,…,αn。選取p(p≤n)個特征值，滿足累計貢獻(xiàn)率≥85%。新樣本φ(xj)投影后的第j(j=1,2,…,p)維坐標(biāo)為：

(7)

(8)

1.1.2 核函數(shù)的選擇

核函數(shù)的類別一般分為局部核函數(shù)和全局核函數(shù)，而核函數(shù)的選擇往往影響KPCA數(shù)據(jù)降維的效果。為更加充分地保留數(shù)據(jù)特征信息，本文將兩類核函數(shù)加權(quán)組合從而建立新的核函數(shù)。其中核函數(shù)參數(shù)通過多重實驗的方法得出，組合核函數(shù)比例系數(shù)采用網(wǎng)格搜索法選取最優(yōu)組合系數(shù)，使其貢獻(xiàn)率達(dá)到最大值。核函數(shù)選用高斯核函數(shù)(局部核函數(shù))和多項式核函數(shù)(全局核函數(shù))：

(9)

k(x,x′)=(γ(x,x′)+c)d

(10)

式中，σ為高斯核函數(shù)的待定參數(shù)，γ、c、d為多項式核函數(shù)的待定參數(shù)。

將式(9)與式(10)加權(quán)組合得：

η2×(γ(x,x′)+c)d

(11)

式中，η1、η2為組合核函數(shù)的比例系數(shù)。

1.2 WOA理論基礎(chǔ)

WOA[14-15]是受座頭鯨捕食行為的啟發(fā)而提出的一種群智能算法。座頭鯨在發(fā)現(xiàn)獵物后，會下潛并旋繞，對獵物進行泡泡網(wǎng)攻擊，如圖1所示。該過程主要經(jīng)歷三個階段：隨機捕獵、環(huán)繞獵物以及捕食獵物。

圖1 座頭鯨泡泡網(wǎng)攻擊捕獵行為示意Fig.1 Humpback whale bubble net attack and hunting behavior

1.2.1 隨機捕獵

模擬鯨魚群隨機捕獵的方式求可行解，數(shù)學(xué)模型如下：

D=|C·X*(t)-X(t)|

(12)

X(t+1)=X*(t)-A·D

(13)

式中，X*為當(dāng)前最優(yōu)位置向量，t為當(dāng)前迭代數(shù)，D為位置衡量參數(shù)，X為位置向量，A和C為控制參數(shù)向量，可通過下式計算：

A=2a·r-a

(14)

C=2·r

(15)

(16)

式中：r為任意向量，tMaxlter為最大迭代數(shù)，a隨迭代次數(shù)的增加由2減少至0，從而實現(xiàn)收縮環(huán)繞機制。

1.2.2 環(huán)繞獵物

環(huán)繞更新位置通過螺旋方程來表示：

D′=|X*(t)-X(t)|

(17)

X(t+1)=D′·ebl·cos(2πl(wèi))+X*(t)

(18)

式中：D′為最佳位置距當(dāng)前位置的距離，l為[-1,-2]中任意數(shù)；b為對數(shù)螺旋形狀的常數(shù)，通過公式(19)計算。

(19)

求解時收縮環(huán)繞和螺旋更新同時進行，各賦予兩種方式1/2的概率執(zhí)行：

(20)

式中：p為[0,1]的隨機數(shù)。

1.2.3 捕食獵物

探索過程中，通過全局搜索更新最佳位置以達(dá)到最優(yōu)，數(shù)學(xué)描述為：

D=|C·Xrand-X|

(21)

X(t+1)=Xrand-A·D

(22)

式中：Xrand是當(dāng)前一代中的隨機位置向量。其中|A|<1時，使用收縮環(huán)繞局部尋優(yōu)；|A|≥1時，使用探索全局尋優(yōu)。

1.3 KELM理論基礎(chǔ)

ELM[16]是一種單隱含層前饋神經(jīng)網(wǎng)絡(luò)，其學(xué)習(xí)目標(biāo)函數(shù)F(x)可用矩陣表示為：

F(x)=h(x)×β=H×β=L

(23)

β=H*·L

(24)

式中：x為輸入向量，h(x)、H為隱層節(jié)點輸出，β為輸出權(quán)重，L為期望輸出H*為H的廣義逆矩陣。

引入懲罰系數(shù)C和單位矩陣I增強神經(jīng)網(wǎng)絡(luò)的穩(wěn)定性，則輸出權(quán)值的最小二乘解為：

(25)

引入核函數(shù)到ELM中，核矩陣為：

ΩELM=HHT=h(xi)h(xj)=K(xi,xj)

(26)

式中：xi，xj為試驗輸入向量，則可將式(23)表達(dá)為：

(27)

式中：(x1,x2…,xn)為給定訓(xùn)練樣本，n為樣本數(shù)量。核函數(shù)選用高斯核函數(shù)即：

(28)

懲罰系數(shù)C和核函數(shù)參數(shù)s的設(shè)置在一定程度上影響KELM[17]的預(yù)測性能。

1.4 KPCA-WOA-KELM模型構(gòu)建

本文巖爆案例基于Zhou等收集的國內(nèi)外地下工程和礦山中245組巖爆案例[18]。去除部分埋深數(shù)據(jù)缺失案例后，選取巖爆烈度判別指標(biāo)并隨機將數(shù)據(jù)集分為訓(xùn)練集174例與測試集30例。

首先利用KPCA對巖爆數(shù)據(jù)做特征提取，然后將訓(xùn)練集輸入KELM進行訓(xùn)練，并使用WOA對其參數(shù)優(yōu)化，最后利用測試集測試該預(yù)測模型的準(zhǔn)確性?；贙PCA-WOA-KELM的巖爆烈度預(yù)測模型的流程如圖2所示：

圖2 巖爆烈度預(yù)測模型流程Fig.2 The process of rockburst intensity prediction model

2 巖爆烈度仿真實驗與分析

2.1 巖爆烈度預(yù)測指標(biāo)選取

影響巖爆的因素具有顯著的突發(fā)性、隨機性與復(fù)雜性的特點，其主要包括巖性、地質(zhì)條件、開挖方式、圍巖初始應(yīng)力狀態(tài)、隧道斷面型式等。

目前國內(nèi)外學(xué)者對巖爆預(yù)測指標(biāo)的選取通常基于巖石力學(xué)參數(shù)或巖爆判據(jù)。一般來說，選取巖爆預(yù)測指標(biāo)應(yīng)符合易于獲取、代表性與可操作性強的特點，能從多個方面反映巖爆的特征信息。

根據(jù)巖爆的影響因素、特點及成因，本文選取巖爆埋深D，圍巖最大切應(yīng)力σθ，單軸抗壓強度σc，單軸抗拉強度σt，應(yīng)力集中系數(shù)SCF(最大切向應(yīng)力比單軸抗壓強度)，脆性系數(shù)(單軸抗壓強度比單軸抗拉強度B1和單軸抗壓、抗拉強度之差與兩者之和的比值B2)，彈性能量指數(shù)Wet作為巖爆烈度的主要評判指標(biāo)。將巖爆分為4級：無巖爆、輕微巖爆、中等巖爆、強烈?guī)r爆，其級別分別對應(yīng)Ⅰ級、Ⅱ級、Ⅲ級、Ⅳ級，表1為不同理論判據(jù)下巖爆烈度分類標(biāo)準(zhǔn)表[19]。

表1 不同理論判據(jù)的巖爆烈度分類標(biāo)準(zhǔn)

2.2 KPCA特征提取

KPCA特征提取步驟為：

(1)選取巖爆數(shù)據(jù)。

(2)歸一化初始數(shù)據(jù)矩陣。

(3)使用核技巧將樣本空間映射至高維空間構(gòu)建核矩陣Kn×n。

(4)核矩陣中心化，即：

(29)

其中，In為元素均為1/n的矩陣。

(5)求解核矩陣的特征值λ，排序特征向量α。選取出p(p≤n)個特征向量，并根據(jù)式(8)規(guī)范化。

(6)通過特征值和特征向量得出主成分，樣本xj在特征向量αk上的投影為：

k=1,2,…,p

(30)

構(gòu)造特征提取后的矩陣Rn×p。

分別通過組合核函數(shù)、高斯核函數(shù)和多項式核函數(shù)實現(xiàn)KPCA特征提取，其中高斯核函數(shù)參數(shù)σ=1.81，多項式核函數(shù)中d=3，γ=c=1，組合系數(shù)η1=0.7，η2=0.3。

表2為通過不同核函數(shù)分別得到的特征值，貢獻(xiàn)率和累計貢獻(xiàn)率?？梢钥闯?，組合核KPCA得到70.43%的第一主成分貢獻(xiàn)率，遠(yuǎn)大于高斯核、多項式核KPCA，說明相比于單一核KPCA，組合核KPCA具有更好的特征提取能力。

表2 不同核函數(shù)特征提取結(jié)果對比

2.3 WOA優(yōu)化KELM模型

由于KELM對其核函數(shù)參數(shù)的敏感性，直接使用KELM建模，往往導(dǎo)致預(yù)測結(jié)果產(chǎn)生較大誤差。因此，為提高KELM的預(yù)測性能，使用WOA算法對KELM的懲罰系數(shù)C和核函數(shù)參數(shù)s進行全局尋優(yōu)。具體步驟如下：

(1)歸一化KPCA降維后的數(shù)據(jù)集。

(2)參數(shù)設(shè)置：變量數(shù)dim=2、鯨魚數(shù)量SearchAgents_no=30、最大迭代次數(shù)tMaxlter=150，變量下限lb=[30,0.2]，變量上限ub=[500,0.7]。

(3)種群初始化：令n=SearchAgents_no，鯨魚位置為：

(31)

通過下式計算鯨魚初始隨機種群位置X0，迭代數(shù)t=1：

X0(i,j)=(ub(i)-lb(i))×rand(i,j)+lb(i)

(32)

式中：X0(i,j)為式(31)第i行第j列的值，lb(i)和ub(i)為第i個座頭鯨的下限和上限，rand(i,j)為[0,1]內(nèi)的隨機數(shù)。

(5)利用WOA尋優(yōu)得到最優(yōu)Cbest和sbest如表3所示。

表3 不同核KPCA下KELM最優(yōu)參數(shù)

(6)將Cbest和sbest代入KELM得到KPCA-WOA-KELM的巖爆烈度等級預(yù)測模型。

2.4 KPCA-WOA-KELM模型實測

將30例測試數(shù)據(jù)輸入模型中，結(jié)果如圖3所示。

圖3 不同核KPCA-WOA-KELM巖爆烈度預(yù)測結(jié)果Fig.3 KPCA-WOA-KELM rockburst intensity prediction results for different core

該模型相關(guān)評估指標(biāo)如表4所示。分別采用F值、準(zhǔn)確率(Accuracy)、精確率(Precision)、召回率(Rcall)量化KPCA-WOA-KELM的預(yù)測性能。其中Accuracy可量化模型整體的預(yù)測性能，其余指標(biāo)則反映模型對4類巖爆等級中單獨的某一級的預(yù)測能力。

表4 不同核KPCA-WOA-KELM性能評價

組合核、高斯核和多項式核的KPCA-WOA-KELM的準(zhǔn)確率分別為90%、86.66%、83.33%，精確率分別為87.05%、83.92%、84.37%，召回率分別為95%、93.18%、87.62%，F(xiàn)值分別為90.06%、86.69%、83.3%，說明該模型有較高的預(yù)測性能。其中，組合核KPCA-WOA-KELM預(yù)測能力最優(yōu)，對4類巖爆等級分類精確率分別為85.71%～100%，85.71%～100%，87.5%～100%，召回率分別為80%～100%，72.73%～100%，72.73%～100%，F(xiàn)值分別為85.71%～93.33%，76.92%～93.33%，66.67%～100%。其中Ⅰ級巖爆與Ⅱ級巖爆的F值均大于Ⅲ級巖爆與Ⅳ級巖爆的F值，說明該模型對Ⅰ級巖爆與Ⅱ級巖爆的預(yù)測性能較強，而對于Ⅲ級巖爆與Ⅳ級巖爆的預(yù)測性能較弱。

綜上所述，綜合考慮該模型的準(zhǔn)確率及相關(guān)評價指標(biāo)，組合核KPCA-WOA-KELM具有較好的預(yù)測性能。

3 工程實踐

3.1 工程概況

秦嶺終南山公路隧道全程18.02 km，最大縱坡11%。該工程采用雙洞雙線設(shè)計，凈寬10.5 m、限高5 m，安全等級一級，結(jié)構(gòu)設(shè)計基準(zhǔn)期100 a。該隧道深埋地段的最大水平主應(yīng)力達(dá)到21.04 MPa，地應(yīng)力水平較高，巖爆發(fā)生概率較大[20]。

3.2 模型應(yīng)用

將KPCA-WOA-KELM模型應(yīng)用于秦嶺終南山隧道工程巖爆烈度的預(yù)測。該工程巖爆相關(guān)指標(biāo)如表5所示，預(yù)測結(jié)果與實際情況基本一致。表明KPCA-WOA-KELM在實際工程中有一定的可靠性，為實際工程提供一定的參考價值。

表5 終南山隧道巖爆實測數(shù)據(jù)與預(yù)測結(jié)果

4 結(jié)論

1)在確定巖爆烈度主要評判指標(biāo)后，使用KPCA對巖爆數(shù)據(jù)做非線性降維，相比于PCA,KPCA能更充分保留特征信息，簡化KELM模型的輸入?yún)?shù)，同時借助WOA優(yōu)化KELM的懲罰系數(shù)C和核函數(shù)參數(shù)s使其達(dá)到最優(yōu)，從而提高模型的訓(xùn)練速度與預(yù)測精度。

2)對局部核函數(shù)和全局核函數(shù)采用加權(quán)方式構(gòu)造組合核KPCA，并與WOA優(yōu)化的KELM相結(jié)合，建立組合核KPCA-WOA-KELM模型。使用204例巖爆案例對模型進行訓(xùn)練與測試，并用準(zhǔn)確率、精確率、召回率、F值等指標(biāo)綜合評價KPCA-WOA-KELM的預(yù)測性能。結(jié)果表明，組合核KPCA-WOA-KELM相較于單一核KPCA-WOA-KELM有更高的預(yù)測性能。

3)為檢驗該模型的預(yù)測性能，將其應(yīng)用于終南山公路隧道中對巖爆等級進行預(yù)測，結(jié)果表明與實際巖爆等級基本一致，說明該模型具有一定的工程應(yīng)用價值。