趙紅夢，姜志俠

(長春理工大學(xué)理學(xué)院，長春130000)

隨著科技的發(fā)展與進(jìn)步，礦場開采等活動所引起的地面爆破振動對生產(chǎn)和生活環(huán)境的影響引起大眾的廣泛關(guān)注；目前評估地面振動以及損害最重要的兩個(gè)參數(shù)為爆破質(zhì)點(diǎn)峰值振動速度(PPV)和頻率[1]。傳統(tǒng)的地面振動研究方法是通過對爆破振動相關(guān)有效數(shù)據(jù)進(jìn)行回歸分析，然后通過回歸分析得到的薩道夫斯基經(jīng)驗(yàn)公式進(jìn)行預(yù)測[2]，現(xiàn)實(shí)中質(zhì)點(diǎn)峰值振動速度本身與很多因素有關(guān)，而常用的經(jīng)驗(yàn)公式考慮的影響因素較少，導(dǎo)致其預(yù)測結(jié)果的誤差較大。

基于常用公式的局限，眾多學(xué)者提出應(yīng)用BP神經(jīng)網(wǎng)絡(luò)算法解決預(yù)測爆破質(zhì)點(diǎn)峰值振動速度問題[3]。其中郭連軍等[4]利用優(yōu)化的神經(jīng)網(wǎng)絡(luò)建立爆破的預(yù)測模型；施建俊等[5]將MATLAB與VB技術(shù)相結(jié)合，開發(fā)得到爆破質(zhì)點(diǎn)峰值振動速度預(yù)測模型；孟陸波等[6]基于人工神經(jīng)網(wǎng)絡(luò)算法建立了爆破質(zhì)點(diǎn)峰值振動速度預(yù)測模型；張藝峰等[7]建立了針對爆破質(zhì)點(diǎn)峰值振動速度以及主頻的預(yù)測模型。

針對目前研究現(xiàn)狀，筆者將主成分分析(PCA)與BP神經(jīng)網(wǎng)絡(luò)算法結(jié)合在一起，提出了一種新的預(yù)測爆破質(zhì)點(diǎn)峰值振動速度的算法,記為PCA-BP算法；其中主成分分析為一種利用降維思想提取數(shù)據(jù)中方差貢獻(xiàn)率較大特征的方法[8]。該算法整體步驟為首先采用主成分分析方法對影響爆破質(zhì)點(diǎn)峰值振動速度的眾多因素進(jìn)行分析和簡化，再結(jié)合BP神經(jīng)網(wǎng)絡(luò)算法進(jìn)行預(yù)測，并且將預(yù)測結(jié)果與傳統(tǒng)經(jīng)驗(yàn)公式預(yù)測結(jié)果進(jìn)行對比分析。

1 PCA-BP算法提出

1.1 PCA算法

PCA算法最核心的思想就是降維，是通過正交變換將指標(biāo)減少的多元統(tǒng)計(jì)分析方法[9]。

現(xiàn)有p個(gè)變量x1,x2,…xp，n個(gè)樣品的數(shù)據(jù)陣為

(1)

(2)

將p個(gè)變量綜合為p個(gè)新的綜合變量，即

(3)

1.2 BP神經(jīng)網(wǎng)絡(luò)算法

BP算法由輸入層、輸出層、隱藏層3部分構(gòu)成，BP算法的學(xué)習(xí)過程由正向、反向傳播組成[11]。

設(shè)輸入模式為x=(x1,x2,…,xn)T，算法設(shè)置隱藏層有h個(gè)單元，隱藏層的輸出為y=(y1,y2,…,yh)T，輸出層有m個(gè)單元，網(wǎng)絡(luò)輸出為z=(z1,z2,…,zm)T，目標(biāo)輸出為t=(t1,t2,…,tm)T，隱藏層到輸出層的傳遞函數(shù)為f,輸出層的傳遞函數(shù)為g，于是可得：

(4)

式中：yj為隱藏層第j個(gè)神經(jīng)元的輸出;w0j=θ,x0=-1。

(5)

式中：zk為輸出層第k個(gè)神經(jīng)元的輸出;此時(shí)網(wǎng)絡(luò)輸出與目標(biāo)輸出的誤差為

(6)

1.3 PCA-BP算法

首先采用PCA方法對影響爆破質(zhì)點(diǎn)峰值振動速度的參數(shù)進(jìn)行降維處理[12]，再以PCA算法計(jì)算的綜合評價(jià)得分作為BP算法輸入的初始數(shù)據(jù)，最后通過BP神經(jīng)網(wǎng)絡(luò)算法對爆破質(zhì)點(diǎn)峰值振動速度進(jìn)行預(yù)測。

PCA-BP算法主要步驟如下：

1)對地面爆破振動的輸入數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理；

2)計(jì)算數(shù)據(jù)相關(guān)系數(shù)矩陣的特征向量；

3)選擇p個(gè)特征向量構(gòu)成主成分；

4)由主成分綜合評價(jià)得分和BP參數(shù)隨機(jī)選取第k個(gè)輸入數(shù)據(jù)作為輸入及對應(yīng)期望輸出；

5)利用隱藏層各神經(jīng)元的輸入和輸出，計(jì)算誤差函數(shù)ε對輸出層的各神經(jīng)元的偏導(dǎo)數(shù)和誤差函數(shù)對隱含層各神經(jīng)元的偏導(dǎo)數(shù)，分別為

6)計(jì)算從隱藏層到輸出層的權(quán)值

wjk(t+1)=wjk(t)+ηδkyj，

再計(jì)算從輸入層到隱藏層的權(quán)值

wij(t+1)=wij(t)+ηδjxi；

7)判斷目標(biāo)誤差是否滿足要求，當(dāng)誤差達(dá)到預(yù)設(shè)精度則停止計(jì)算；否則，返回到第5步，進(jìn)入下一輪學(xué)習(xí)。

2 案例分析

現(xiàn)引用文獻(xiàn)[5]中的20組實(shí)測地面爆破振動數(shù)據(jù)(見表1)進(jìn)行研究，影響爆破振動的參數(shù)比較多，將爆心距x1、高程差x2、總藥量x3、炮孔深度x4、單段最大量x5作為研究的自變量，以爆破質(zhì)點(diǎn)峰值振動速度作為因變量，通過PCA方法提取主成分，利用BP神經(jīng)網(wǎng)絡(luò)算法進(jìn)行預(yù)測。

表1 爆破振動數(shù)據(jù)

表1中，選擇前15組數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，后5組作為測試數(shù)據(jù)，通過訓(xùn)練數(shù)據(jù)得出的結(jié)論預(yù)測后5組數(shù)據(jù)的爆破質(zhì)點(diǎn)峰值振動速度；并將預(yù)測結(jié)果與傳統(tǒng)薩道夫斯基經(jīng)驗(yàn)公式所得結(jié)果進(jìn)行對比分析。

2.1 相關(guān)性分析

首先進(jìn)行原始數(shù)據(jù)的標(biāo)準(zhǔn)化處理，標(biāo)準(zhǔn)化處理公式為

(7)

式中：x′為標(biāo)準(zhǔn)化后的數(shù)據(jù);x為原始數(shù)據(jù);u為樣本均值;σ為樣本標(biāo)準(zhǔn)差.

根據(jù)正態(tài)分布的基本性質(zhì)，有x～N(u,σ2)，則有

(8)

原始數(shù)據(jù)集標(biāo)準(zhǔn)化為接近于正態(tài)分布(均值為0，方差為1)的數(shù)據(jù)集，然后對其數(shù)據(jù)進(jìn)行相關(guān)性分析，得到相關(guān)系數(shù)如表2所示。

表2 相關(guān)系數(shù)

由表2知單段最大藥量與總藥量、炮孔深度相關(guān)系數(shù)較高，所以適合做主成分分析。

2.2 主成分分析

1)主成分選取。首先利用碎石圖(見圖1)選取主成分。

圖1 碎石圖Fig.1 Break store chart

由圖1知，組件序號在4到5時(shí)圖形區(qū)域平緩，初步判定可以選擇3～4個(gè)主成分,其次進(jìn)行累計(jì)方差貢獻(xiàn)率研究(見表3)。

表3 累計(jì)方差貢獻(xiàn)率

可以看出在選擇3個(gè)主成分時(shí)候，保留了原始變量92.641%的信息。利用新的3個(gè)主成分代替原來的5個(gè)變量，從而起到降維的作用[13]。

2)主成分得分。第一步進(jìn)行因子載荷陣的求取，利用軟件得到因子載荷矩陣(見表4)。

表4 因子載荷陣

根據(jù)因子載荷陣計(jì)算主成分得分公式：

(9)

從而根據(jù)主成分得分公式以及方差貢獻(xiàn)率分別計(jì)算20組數(shù)據(jù)的綜合評價(jià)得分為：-7.04、3.21、-11.25、-5.56、-18.33、-11.53、2.38、10.6、18.76、2.38、10.66、18.82、14.15、22.38、-29.47、-21.27、-13.17、1.47、7.29、5.51。

2.3 PCA-BP算法

將主成分分析提取的3個(gè)主成分，經(jīng)計(jì)算得綜合評價(jià)得分作為BP算法初始數(shù)據(jù)的輸入，通過訓(xùn)練組數(shù)據(jù)(前15組)得出的結(jié)論預(yù)測測試組(后5組)數(shù)據(jù)。構(gòu)建一個(gè)3層BP神經(jīng)網(wǎng)絡(luò)對地面振動后5組數(shù)據(jù)(見表1)的爆破質(zhì)點(diǎn)峰值振動速度進(jìn)行預(yù)測，其中輸入層結(jié)點(diǎn)數(shù)為3個(gè)，隱藏層結(jié)點(diǎn)數(shù)為8個(gè)，激活函數(shù)選擇‘tansig’;輸出層結(jié)點(diǎn)數(shù)為1個(gè)，激活函數(shù)選擇‘logsing’。

采用梯度下降動量和自適應(yīng)lr算法‘tansig’訓(xùn)練BP網(wǎng)絡(luò)。學(xué)習(xí)率為0.05，目標(biāo)誤差為0.625×10-3，最大迭代次數(shù)為1 000，使用MATLAB軟件進(jìn)行計(jì)算。前15個(gè)網(wǎng)絡(luò)輸出與實(shí)際輸出對比如圖2所示。

圖2 前15個(gè)網(wǎng)絡(luò)輸出與實(shí)際輸出對比Fig.2 Comparison of the first 15 network outputs and the actual output

地面爆破質(zhì)點(diǎn)峰值振動速度預(yù)測結(jié)果與期望目標(biāo)的關(guān)系如圖3所示。

注：實(shí)直線為預(yù)測結(jié)果擬合的線性曲線，散點(diǎn)圖構(gòu)成的虛線為實(shí)際值擬合而成。圖3 振速預(yù)測與期望目標(biāo)關(guān)系Fig.3 Relationship between vibration speed prediction and desired target

利用PCA-BP算法經(jīng)MATLAB計(jì)算可得測試數(shù)據(jù)(后5組)預(yù)測結(jié)果值分別為2.349 2、4.929 6、0.771 2、1.150 3、1.421 4。將預(yù)測結(jié)果與期望目標(biāo)擬合成線性函數(shù)，計(jì)算得二者的相關(guān)系數(shù)達(dá)到 0.896 75(理想狀態(tài)是相關(guān)系數(shù)等于1時(shí))[13]，說明預(yù)測結(jié)果的可信度較高。

3 對比分析

傳統(tǒng)薩道夫斯基經(jīng)驗(yàn)公式[14]為

(10)

式中：k,α為經(jīng)驗(yàn)系數(shù)；Q為單段藥量；R為測點(diǎn)到爆心的距離。

從而利用對數(shù)變換

(11)

令

公式變?yōu)榛貧w方程

y=b+αx

從而得到預(yù)測公式為

y=2.551+1.747x

(12)

利用傳統(tǒng)經(jīng)驗(yàn)公式計(jì)算后5組數(shù)據(jù)預(yù)測結(jié)果為1.60、1.34、0.91、0.97、0.83，具體回歸分析結(jié)果如表5所示。

表5 回歸分析結(jié)果

將傳統(tǒng)薩道夫斯基經(jīng)驗(yàn)公式[15]以及PCA-BP算法預(yù)測結(jié)果進(jìn)行對比分析(見表6)，這里以測試數(shù)據(jù)(后5組)預(yù)測結(jié)果的相對誤差、平均誤差作為參照，可看出利用PCA-BP算法預(yù)測的相對誤差分別為2.52%、2.83%、7.94%、14.80%、10.65%，平均相對誤差為7.748%；傳統(tǒng)經(jīng)驗(yàn)公式的相對誤差較高，分別為33.61%、72.03%、21.97%、1.02%、34.64%，平均相對誤差為32.654%；除了第19組的相對誤差略高于用傳統(tǒng)經(jīng)驗(yàn)公式預(yù)測的相對誤差外，其他組利用PCA-BP算法預(yù)測的相對誤差均較?。贿@說明PCA-BP算法預(yù)測模型有效性較高。

表6 預(yù)測結(jié)果對比分析

4 結(jié)語

PCA-BP算法解決了在預(yù)測爆破質(zhì)點(diǎn)峰值振動速度中輸入數(shù)據(jù)較多，各個(gè)輸入變量之間存在多重共線性和相關(guān)性的問題。在案例分析中，通過對輸入數(shù)據(jù)進(jìn)行主成分分析提取主成分特征向量，計(jì)算綜合評價(jià)得分作為BP算法的輸入，并與傳統(tǒng)薩道夫斯基經(jīng)驗(yàn)公式預(yù)測相比，PCA-BP 算法綜合考慮影響地面爆破振動的因素，使得預(yù)測結(jié)果更加合理，結(jié)合主成分分析和BP神經(jīng)網(wǎng)絡(luò)算法的預(yù)測結(jié)果平均相對誤差為7.748%，小于用薩道夫斯基經(jīng)驗(yàn)公式進(jìn)行預(yù)測的平均相對誤差32.654%，故說明，將PCA-BP算法應(yīng)用到爆破振動工作中是比較可行的。