基于改進PSO-BP神經網絡的爆破振動速度峰值預測

范 勇，裴 勇，楊廣棟，冷振東，2，盧文波

(1.三峽大學 湖北省水電工程施工與管理重點實驗室，湖北 宜昌 443002；2.中國葛洲壩集團易普力股份有限公司 重慶市民用爆器材工程技術研究重心，重慶 401121；3.武漢大學 水工巖石力學教育部重點實驗室，武漢 430072)

由于鉆爆法的經濟性和高效性，依然被廣泛地應用于巖體破碎、礦山開采等方面。經研究，在爆破的過程中炸藥所產生的能量只有20%～30%被利用，其余的70%～80%被以其他的方式浪費掉，并產生了多種負面影響，例如爆破飛石、爆破振動、空氣超壓等，其中爆破振動被認為對結構最有害的負面影響。

目前，爆破振動質點速度峰值(peak particle velocity，PPV)被用來評估爆破振動是否給結構帶來破壞的重要指標。準確地預測PPV來優化爆破設計參數，可以有效降低爆破振動帶來的不利影響[1]。以往的研究表明影響PPV的因素有爆破參數(如孔深、孔徑、最大單響藥量等)、巖體性質參數(如巖體的縱波波速[2]等)和地形地貌[3]。為了準確的預測PPV，許多研究學者提出了經驗公式，如我國和俄國使用的是由前蘇聯學者提出的薩道夫斯基公式、美國使用的美國礦務局公式、印度使用的印度標準局公式等。但這些經驗公式主要考慮了2個主要影響參數，即爆心距和最大單響藥量，將其他影響參數歸為了公式中的經驗系數，已有的研究表明經驗公式具有一定的局限性，存在預測的準確度低[4-5]。為了提高預測的準確度也有研究學者提出改正后的預測公式，如駱曉鋒等[6]基于量綱分析理論在傳統薩道夫斯基公式基礎上提出反映高程變化的爆破振動的改進公式，并通過對比改進公式與薩道夫斯基公式的預測誤差進一步論證改進公式的可靠性。

薩道夫斯基公式

(1)

式中:V為爆破振動質點速度峰值；r為爆心距；Q為最大單響藥量；k，α分別為場地系數和衰減系數。

經驗公式使用的影響因素有限，而且PPV與影響因素存在著復雜的非線性關系，預測的結果有時存在較大的誤差，為了解決公式的不足，有必要使用新的預測方法。機器學習在解決非線性關系的問題上體現出明顯的優勢，而且可以將更多的影響因素作為輸入參數，使得預測結果更加接近實測值。例如：Xu等[7]利用主成分分析從垂直距離、炸藥類型、裝藥量等8個輸入因子中提取前4個主成分作為輸入參數，建立了GA-ANN(genetic algorithm-artificial neural network )預測模型，預測結果與實測值吻合更好；史秀志等[8]把最大單響藥量、水平距離、高程差等作為輸入參數，建立了基因表達式編程預測模型對PPV進行預測，結果表明比經驗公式精度更高；邵良杉等[9]建立了LS-SVM(least squares support vector machine)預測模型，預測露天采礦爆破振動對民房的破壞情況，其結果與實際情況吻合較好；趙紅夢等[10]建立了基于BP(back propagation)神經網絡預測模型，預測結果更接近工程實測值。更多的研究表明[11-13]，BP神經網絡適合PPV的預測。但是BP神經網絡存在一定的缺陷，如初始權值和閾值的取值問題會使得預測的結果不穩定、準確度低等。目前，采用智能算法(如遺傳算法(genetic algorithm，GA)[14]、粒子群算法等)來優化BP神經網絡是一種有效的解決方法。

為了提高BP神經網絡預測模型的準確度，采用了粒子群算法優化BP神經網絡，建立改進的PSO-BP(particle swarm optimization-back propagation)神經網絡預測模型，將其應用于爆破振動速度預測。針對BP神經網絡輸入參數的選擇問題，當前，更多的研究是將爆破參數和距離作為主要的輸入參數，缺少了對傳播場地影響的考慮。本文以白鶴灘水電站左岸壩肩槽爆破開挖為依據，考慮了代表場地條件的縱波波速。對比BP神經網絡以及薩道夫斯基公式檢驗結果，驗證了該模型的合理性和適用性，為輸入參數的選擇提供了參考。

1 PSO-BP神經網絡基本原理

1.1 BP神經網絡

BP神經網絡是一種按照誤差逆向傳播算法訓練的多層前饋神經網絡模型，具有很強的非線性映射能力。網絡的結構包含著輸入層、隱含層、輸出層三部分，如圖1所示。

圖1 3層BP神經網絡拓撲結構Fig.1 The topology structure of 3-layer BP neural network

BP神經網絡包括信號的正向傳播和誤差的反向傳播2個過程。輸入信號從輸入層傳入，經隱含層處理后傳向輸出層，若輸出層的實際輸出與期望輸出不相符，則進入誤差的反向傳播。誤差由輸出層向前逐層反傳，分配給各層神經元，連接神經元的權值和閾值通過各層獲得的誤差進行調整，其過程采用梯度下降法。經過反復的傳播，當滿足最小誤差要求或者最大訓練次數時，網絡訓練停止。

定義誤差函數為

(2)

式中：Z為訓練樣本組數；m為輸出層節點。

BP神經網絡訓練的步驟如下：

步驟1初始化網絡權值wnk、vkj，閾值a、b。

步驟2計算隱含層神經元的輸出值qk。

(3)

式中：k=1,2,…,h；f為隱含層激活函數；xn為第n個輸入信號。

步驟3計算輸出層神經元的輸出值oj。

(4)

步驟4權值更新。

(5)

vkj(t+1)=vkj(t)+ηqk(yj-oj)

(6)

式中，η稱為學習率的參數，一般取值(0，1)區間。

步驟5閾值更新。

(7)

bj(t+1)=bj(t)+(yj-oj)

(8)

步驟6判斷算法是否迭代結束，若沒有結束，則返回步驟2。

采用BP神經網絡預測，需要確定隱含層數與其節點數。Hecht-Nielson在理論上已經證明單個隱含層的BP神經網絡能夠逼近任何區間的連續函數，所以本文神經網絡選用單個隱含層結構。隱含層節點數采用經驗式(9)[15],并根據訓練集的均方誤差(歸一化后)最小值來確定隱含層節點數。

(9)

式中：h為隱含層節點數；i為輸入層節點數；m為輸出層節點數；a取1～10的整數。

1.2 粒子群算法原理

粒子群算法是由 Kennedy 和 Eberhart 于1995年提出的一種基于群體智能的優化方法，基本思想源于對鳥群捕食行為的研究。該算法首先在可行解空間初始化一群粒子，每個粒子代表一個優化問題的潛在解，用位置、速度和適應度值表示該粒子特征，粒子在每個維度的速度和位置變化范圍限定在最大邊界。群體中的粒子在每次迭代搜索過程中，通過跟蹤群體的2個最佳位置不斷地調整自己的速度和位置來搜索新解[16]。

假設在一個D維空間中，有n個粒子：

第i個粒子的位置：xi=(xi1,xi2,…,xiD),代表尋優問題在D維空間的一個潛在的解向量,將粒子i代入適應度函數求適應值；

第i個粒子的速度：vi=(vi1,vi2,…,viD)；

第i個粒子經歷過的最佳位置：pbest=(pi1,pi2,…,piD)；

群種所經歷過的最佳位置：gbest=(gi1,gi2,…,giD)；

粒子i的第d維速度和位置更新公式為

(10)

(11)

1.3 PSO-BP神經網絡模型

BP神經網絡雖然有很強的非線性映射能力，但誤差的反向傳播是基于梯度下降法來調整網絡的連接權值和閾值，常常會陷入局部最優解，使得預測的準確度降低。PSO-BP神經網絡模型利用了PSO算法的全局尋優能力，在其解空間內搜索網絡的權值、閾值。粒子在每一維上的位置就是所求的解，要建立粒子的維數和權值、閾值的映射關系。單個隱含層的BP神經網絡維數如下

D=i×h+h+h×m+m

(12)

式中：D為空間維數；i為輸入層節點；h為隱含層節點；m為輸出層節點。

粒子位置的優劣由適應度函數確定，文中選擇訓練集與測試集整體均方誤差的平均值為適應度函數。適應度值越小，表明訓練越準確，且兼顧模型的檢驗結果更好。適應度函數如下公式(單個輸出節點)。

(13)

式中：Z，N分別為訓練集、測試集樣本組數；yp，op分別對應訓練樣本下輸出層第p個實際輸出和期望輸出；Yl，Ol分別對應測試樣本下輸出層第l個實際輸出和期望輸出。

BP神經網絡采用PSO算法優化后的權值和閾值進行訓練和檢驗。具體流程圖如圖2所示。

圖2 PSO-BP神經網絡流程圖Fig.2 PSO-BP neural network flow chart

1.4 模型評估

在這項研究中，需要對模型進行評估。從統計學的角度來看，僅用一個性能指標進行評價是不夠的。因此，采用決定系數(R2)、平均絕對誤差(mean absolute error,MAE)、均方根誤差(rooot mean square error,RMSE)和平均絕對百分比誤差(mean absolute percentage error,MAPE)4個性能指標對模型綜合評估。R2表示實測值與預測值的線性相關性，MAE用來表示結果的偏差，RMSE用來表示結果的離散度，MAPE表示結果的準確度。訓練的評估采用R2、RMSE,檢驗的結果采用RMSE、MAE和MAPE進行評估[17]。計算公式為

(14)

(15)

(16)

(17)

2 爆破振動速度峰值預測

2.1 樣本數據收集

2.1.1 工程概況

數據來自白鶴灘水電站左岸壩肩槽833.7 m，820.5 m，800.0 m，790.0 m和780.0 m爆破開挖所搜集到的63組數據集。白鶴灘水電站位于金沙江下游四川省寧南縣和云南省巧家縣交界的峽谷中，壩型為混凝土雙曲拱壩，壩高289.0 m，設計正常蓄水位825.0 m，壩頂高程834.0 m，是金沙江下游4個梯級電站中的第二級。

金沙江該段總體由南向北流，左岸為大涼山山脈東南坡，整體上呈向金沙江傾斜的斜坡地形，巖體結構主要是在玄武巖巖流層中形成的層間、層內構造錯動帶以及斷層構造裂隙系統，巖體的質量參次不齊，如圖3所示。隨著工程的進展，開挖揭示出新的地質條件。為確保壩肩槽巖體爆破開挖質量，設計了爆前、爆后聲波檢測，并開展了現場爆破振動監測工作，對每次爆破效果進行檢測。目的是對爆破的效果進行評價，為壩肩槽的開挖提供合適的爆破參數。

圖3 白鶴灘水電站工程地質剖面圖Fig.3 Engineering geological profile of Baihetan Hydropower Station

2.1.2 爆破振動監測

白鶴灘水電站爆破振動監測采用成都中科測控有限公司生產的爆破監測儀(TC-4850)。在監測點布置一臺三向速度傳感器(可同時測量豎直向、水平徑向和水平切向PPV，文中采用的是水平徑向PPV)，用石膏將傳感器固定在所需監測的點位，然后將自記儀與其相聯。爆破振動信號傳遞到監測點時，自記儀自動記錄信號。爆后利用爆破振動分析軟件將自記儀采集到的振動信號輸入電腦中進行分析處理。

以高程820.5 m爆破開挖為例，根據現場地形條件，在爆區后方布置監測點如圖4(a)所示，圖4(b)中的數字分別代表測點編號及爆心距。現場爆破統一采用乳化炸藥，起爆順序為預裂孔、爆破孔、緩沖孔，預裂孔為5～7孔一響，爆破孔和緩沖孔為2孔一響，由于預裂孔先起爆產生預裂縫，削減了爆破孔和緩沖孔產生的爆破振動，在裝藥量上緩沖孔小于爆破孔，因此只關注預裂孔和爆破孔的爆破振動速度。具體的爆破參數如表1所示，炮孔布置如圖4(c)所示。記錄到的波形如圖5所示。

圖4 監測點及炮孔布置Fig.4 Layout of monitoring points and blasting parameters

圖5 水平徑向振動波形Fig.5 Horizontal radial vibration waveform

表1 820.5 m爆破參數Tab.1 820.5 m blasting parameters

2.1.3 巖體聲波檢測

聲波檢測儀器為HX-SYB型智能型巖石聲波儀，換能器為40 kHz單孔一發雙收換能器。

具體聲波檢測步驟為：連好設備，將換能器置于檢測孔底部，向檢測孔注水至水流出；操作聲波儀進行檢測、讀數并記錄；移動換能器到下一檢測位置，重復上述步驟。檢測孔布置及現場爆前、爆后聲波檢測如圖6所示(833.7 m聲波檢測)。檢測到的縱波波速如圖7所示。

圖6 聲波孔布置及現場聲波檢測Fig.6 Acoustic wave hole layout and field acoustic wave detection

圖7 爆后單孔縱波波速Fig.7 The longitudinal wave velocity of single hole after explosion

2.2 模型參數的選擇

輸入參數的選擇是BP神經網絡一個重要因素。為了建立一個全面而準確的模型，需要確定影響PPV的主要輸入參數，同時應考慮所選參數必須代表現場條件和爆破設計參數必須是可獲的。

表2、表3是白鶴灘水電站左岸壩肩槽高程833.7 m與高程817.0 m兩次爆破開挖收集到的爆破振動數據，兩次爆破地形不同。采用薩道夫斯基公式擬合。

表2 833.7 m預裂爆破振動數據Tab.2 833.7 m vibration data of presplitting blasting

表3 817 m預裂爆破振動數據Tab.3 817 m vibration data of presplitting blasting

將薩道夫斯基公式兩邊取對數轉化為線性公式

(18)

y=?x+b

(19)

采用最小二乘法對兩次爆破振動數據進行回歸擬合，擬合結果如圖8所示。

圖8 薩道夫斯基公式擬合曲線Fig.8 Curve fitting of Sadovsky formula

高程833.7 m和817.0 m擬合的相關系數分別為R2=0.863，R2=0.104。從擬合的結果可以得出，采用薩道夫斯基公式預測監測點與爆區在同水平面的PPV具有很好的適用性。當爆區與監測點產生高差時，薩道夫斯基公式并不再適用于預測PPV。因此，將高程差(H)作為一個影響參數。

通過擬合得到高程833.7 m的系數k=41.513，α=1.068。回歸分析后得到的薩道夫斯基公式為

(20)

采用高程833.7 m擬合得到的k、α對高程817.0 m的PPV預測，預測結果如表4所示。

表4 817.0 m PPV實測值與預測值對比Tab.4 Comparison between measured value and predicted value of 817.0 m PPV

從表4實測值與預測值的對比結果可以看出，誤差整體在30%以上，準確度低。可以得出振動速度在傳播的時候受到了場地條件的影響。縱波波速可以很好的反應巖體的裂隙、結構面發育情況等場地特征。因此，將縱波波速作為振動速度傳播的影響參數。本文采用的是爆破后所測得的6個檢測孔的平均值作為輸入參數。

經分析選擇了爆心距、最大單響藥量、高程差和縱波波速作為預測模型的輸入參數。

2.3 敏感性分析

為識別各輸入參數對PPV的相對影響，采用Yang等[18]提出的余弦振幅法進行敏感性分析，評估輸入參數對PPV的影響。兩者之間的關系強度(rij)由式(21)表示

(21)

式中：m為數據個數(文中數據為63組)；xik，xjk分別為輸入參數和輸出參數。每個輸入參數的rij值在0～1變化，最高的rij值表示對PPV最有影響的參數。輸入參數與輸出參數的關系強度計算結果如圖9所示。

圖9 輸入和輸出參數之間的關系強度Fig.9 The strength of the relationship between input and output parameters

敏感性分析結果表明，在傳播的過程中代表場地條件的縱波波速對PPV的影響最顯著，其次是最大單響藥量、爆心距，高程差對PPV影響最小。

2.4 預測方法訓練

基于MATLAB語言構建的預測模型，輸入層節點數為4(即爆心距、最大單響藥量、高程差和縱波波速)，輸出層節點為一個(即PPV)，隱含層節點數采用經驗公式(9)，得到[3,12]10個不同的節點數，每次運行通過對比不同節點數得到訓練集的均方誤差來選取節點數。

BP神經網絡參數設置，訓練次數1 000，學習速率0.01，最小誤差1×10-6，隱含層選用Tansig函數，輸出層為purelin函數，訓練函數采用trainlm。

將所收集到的63組數據集隨機排序，選取前53組作為訓練集進行訓練，后10組為測試集進行檢驗。訓練前對數據進行歸一化處理，歸一化可以消除特征指標的量綱和數量級的影響，有助于提高網絡的學習速度，對于檢驗后的數據還要進行反歸一化。利用式(22)歸一化，區間為[0,1]。

(22)

式中：i為歸一化后的輸入數據；x為實測數據；xmin，xmax分別為實測數據中的最大值與最小值。

2.4.1 PSO-BP神經網絡預測模型

PSO-BP神經網絡預測模型的最佳性能與其慣性權重ω、加速度常數c1，c2和種群規模有著直接的影響，需要對參數合理的選取，收斂性分析[19]的討論為參數選擇提供了依據。

對式(10)、式(11)可簡化到一維

v(k+1)=ωv(k)+φ1[pbest-x(k)]+
φ2[gbest-x(k)]

(23)

x(k+1)=x(k)+v(k+1)

(24)

則參數ω，φ需滿足

ω<1,φ>0,2ω-φ+2>0

(25)

文章選取粒子群規模為200，最大迭代次數1 000，以達到搜索全局的目的。慣性權重ω體現了粒子繼承先前速度的能力。針對PSO算法容易早熟以及算法后期易在全局最優解附近產生振蕩現像，研究學者提出自適應調整的PSO算法，讓ω隨著算法的迭代而遞減。常見的有以下4種方法[20-21]。

(26)

(27)

(28)

(29)

式中：ωmax為慣性權重最大值；ωmin為慣性權重最小值；k為當前迭代次數；Tmax為最大迭代次數。研究表明，當慣性權重ωmax=0.9，ωmin=0.4時算法性能最好。

針對4種不同的ω進行算法性能分析。為了減少模型中隨機數的影響，設置了隨機種子，對不同的慣性權重預測模型都具有相同的隨機數，選取加速度常數c1=c2=2，分析結果如表5所示。

表5 4種慣性權重下的算法性能比較Tab.5 Performance comparison of algorithms under 4 inertia weights

從表5和圖10得到的檢驗結果和最小適應度值對模型綜合分析，選取ω3作為慣性權重來改進PSO算法。

圖10 4種慣性權重的收斂曲線Fig.10 Convergence curves of 4 inertia weights

選取ω3作為慣性權重，加速度常數c1，c2采取不同比值[22]情況下的取值。對于不同加速度常數模型的訓練結果如表6所示。

表6 7種不同加速度常數模型訓練評估Tab.6 Training evaluation of 7 different acceleration constant models

2.4.2 BP神經網絡預測模型

由于BP神經網絡誤差的反向傳播采用的是梯度下降法，易陷入局部最優，因此，采用多次預測的方法。文中采用3次預測，選取最優結果。對每次訓練結果評估如表7所示。

表7 模型訓練評估Tab.7 Model training evaluation

2.5 預測方法檢驗

每次訓練后，采用后10組測試集進行檢驗，并對不同模型檢驗的結果進行評估。BP神經網絡和改進的PSO-BP神經網絡預測模型的檢驗結果如表8、表9所示。

由表8可以得出，BP神經網絡預測模型3檢驗結果最好。由表9檢驗結果和不同加速度模型的收斂曲線如圖11綜合評估，改進的PSO-BP神經網絡預測模型6檢驗結果最優。兩種方法隱含層節點得選取如下表10、表11，通過選取的結果可以發現8個隱含層節點是網絡模型的最優節點，搭建了4×8×1的網絡拓撲結構。此外，由檢驗結果可以得出，改進的PSO-BP神經網絡預測模型的準確度均好于BP神經網絡預測模型，說明了改進的PSO算法具有較好的優化能力。

表8 BP神經網絡檢驗結果Tab.8 BP neural network test results

表9 改進的PSO-BP神經網絡檢驗結果Tab.9 The test results of the improved PSO-BP neural network

圖11 7種不同加速度常數模型的收斂曲線Fig.11 Convergence curves of 7 models with different acceleration constants

表10 BP神經網絡預測模型3隱含層節點選取Tab.10 BP neural network model 3 hidden layer node selection

表11 改進的PSO-BP神經網絡預測模型6隱含層節點選取Tab.11 Improved PSO-BP neural network model 6 hidden layer node selection

2.6 薩道夫斯基公式預測

利用薩道夫斯基公式回歸分析。將前53組數據進行回歸擬合，對后10組檢驗。得到的擬合曲線如圖12所示。

圖12 擬合曲線Fig.12 Curve fitting

通過擬合得到k=33.428，α=1.085，相關系數為0.577。回歸分析后得到的薩道夫斯基公式為

(30)

對后10組的檢驗結果如表12所示。

表12 檢驗結果對比分析Tab.12 Comparative analysis of test results

3 結果分析

薩道夫斯基公式、BP神經網絡預測模型3和改進的PSO-BP神經網絡預測模型6對后10組檢驗的結果對比見表12和圖13。

圖13 檢驗結果對比圖Fig.13 Comparison chart of test results

由表12和圖13可得，改進的PSO-BP神經網絡預測模型的最大誤差為20.94%，MAPE為8.27%；BP神經網絡預測模型的最大誤差為44.45%，MAPE為18.71%；薩道夫斯基公式檢驗的最大誤差為78.74%，MAPE為32.56%，準確度最低。對比BP神經網絡預測模型和薩道夫斯基公式的檢驗結果，改進的PSO-BP神經網絡預測模型的準確度更高，與真實值接近，具有較好的泛化能力。

4 結 論

影響PPV的參數較多，而且存在復雜的非線性關系，采用經驗公式預測準確度低。本文建立了改進的PSO-BP神經網絡預測模型，將其應用于白鶴灘水電站左岸壩肩槽爆破開挖PPV的預測，并與BP神經網絡及薩道夫斯基公式檢驗結果比較，初步得出以下結論：

(1) 采用改進的PSO-BP神經網絡預測模型，選取了爆心距、最大單響藥量、高程差和縱波波速作為輸入參數，PPV作為輸出參數，通過分析建立的網絡拓撲結構為4×8×1的最佳性能模型，檢驗結果MAPE為8.27%，遠好于BP神經網絡預測的最佳結果18.71%和薩道夫斯基公式32.56%。檢驗結果表明改進的PSO算法優化BP神經網絡預測模型具有較好的泛化能力，可有效的應用于實際工程預測PPV。

(2) 采用余弦振幅法分析爆心距、最大單響藥量、高程差和縱波波速與PPV的關系強度分別為0.503，0.690，0.310和0.787，結果表明PPV在傳播的過程中，代表場地條件的縱波波速對PPV影響最顯著。薩道夫斯基公式僅將爆心距和最大單響藥量作為主要影響參數，檢驗的結果誤差較大，準確度低，不再適用于PPV的預測。

(3) 與經驗公式相比，BP神經網絡在處理非線性關系時具有較好的性能，而且可以較全面考慮影響PPV的參數，使得預測結果與實際值更加吻合。