郭 璇， 馬思遠(yuǎn)， 郭一帆， 張曉新

(1. 北京交通大學(xué) 城市地下工程教育部重點(diǎn)實(shí)驗(yàn)室，北京100044; 2. 北京交通大學(xué) 土木建筑工程學(xué)院， 北京 100044； 3. 華誠(chéng)博遠(yuǎn)工程技術(shù)集團(tuán)有限公司, 北京 100052)

鑒于爆炸事件(近年如2017年4月俄羅斯地鐵爆炸等)對(duì)地下公共服務(wù)設(shè)施產(chǎn)生的巨大危害和嚴(yán)重后果，各國(guó)紛紛出臺(tái)相應(yīng)的防護(hù)設(shè)計(jì)規(guī)范：如印度國(guó)家標(biāo)準(zhǔn)局1973年制定地下爆炸結(jié)構(gòu)安全性設(shè)計(jì)規(guī)范，美國(guó)陸軍1983年制定《常規(guī)武器防護(hù)設(shè)計(jì)規(guī)范》[1]，美國(guó)國(guó)防部2008年公布結(jié)構(gòu)抵抗偶然爆炸荷載設(shè)計(jì)手冊(cè)等，地下公共設(shè)施的安全防范等級(jí)不斷升級(jí)。

地下結(jié)構(gòu)防護(hù)設(shè)計(jì)中爆炸荷載的確定是關(guān)鍵性難題[2-5]，考慮爆炸波在非線性地層中傳播的復(fù)雜性如：巖土圍巖特性、地下特殊空間中固-液-氣對(duì)爆炸及測(cè)點(diǎn)參數(shù)產(chǎn)生的重要隨機(jī)影響、復(fù)雜的時(shí)空演化效應(yīng)等多因素交叉影響，圍巖介質(zhì)各點(diǎn)的爆炸應(yīng)力即難確定又難復(fù)核。各經(jīng)驗(yàn)公式存在適用范圍離散，地層參數(shù)隨機(jī)或誤差修正及量化困難等問(wèn)題。

由反射系數(shù)法簡(jiǎn)化求取圍巖介質(zhì)爆炸荷載峰值壓力的既有研究如：Ляхов等[6]對(duì)巖石和土中不同形式?jīng)_擊波傳播規(guī)律及計(jì)算方法的研究及實(shí)驗(yàn)；周豐峻等[7]對(duì)飽和土應(yīng)力波傳播黏塑性模型的計(jì)算方法及應(yīng)變硬化特性實(shí)驗(yàn)；Henrych等[8]區(qū)分巖石、黏性土和非黏性土介質(zhì)，以各相體變?yōu)榛A(chǔ)的變形機(jī)理推導(dǎo)絕熱狀態(tài)方程確定其壓力值；葉亞齊等[9]通過(guò)對(duì)砂黏土不同深度的爆炸試驗(yàn)，給出炸藥埋深對(duì)圍巖介質(zhì)壓力的影響及不同埋深峰值壓力經(jīng)驗(yàn)式；王占江等[10]通過(guò)花崗巖化爆實(shí)驗(yàn)，給出按巖性分區(qū)圍巖介質(zhì)的應(yīng)力衰減規(guī)律；Al-Qasimi等[11]對(duì)不同土質(zhì)系列爆炸試驗(yàn)測(cè)定土中一點(diǎn)峰值振動(dòng)速度和峰值孔壓，提出峰值振速的預(yù)測(cè)公式；William等[12]開(kāi)展地下爆炸試驗(yàn)，通過(guò)60組不同土質(zhì)試驗(yàn)測(cè)得不同濕度凝灰?guī)r、花崗巖和淤泥圍巖介質(zhì)的峰值壓力、振速、加速度及相應(yīng)經(jīng)驗(yàn)公式；Leong等[13]對(duì)新加坡殘積土的多組小比例爆炸試驗(yàn)，測(cè)定飽和/非飽和土3種爆炸比例距離的峰值壓力，對(duì)比試驗(yàn)結(jié)果與TM-5-855-1計(jì)算公式的誤差，修正衰減系數(shù)，采用爆炸比例距離立方根的數(shù)據(jù)處理方法替代平方根擬合結(jié)果，提出無(wú)量綱經(jīng)驗(yàn)公式；市野宏嘉等[14]通過(guò)小比例爆炸試驗(yàn)，研究不同飽和度的中砂、山砂、黃土爆炸波的傳播規(guī)律，給出含參數(shù)密度和飽和度的經(jīng)驗(yàn)公式，并與Leong結(jié)果對(duì)比吻合較好；Gefken等[15]通過(guò)爆炸試驗(yàn)箱不同測(cè)點(diǎn)壓力、振動(dòng)速度和加速度等參數(shù)研究歸納四種飽和度砂性土的爆炸波傳播規(guī)律等。

基于以上研究基礎(chǔ)和成果可以看出，科學(xué)分析爆炸圍巖介質(zhì)波的傳播規(guī)律，快速得到大量峰值壓力的關(guān)鍵特征值，提取隨機(jī)特征，迫切需要提升一種穩(wěn)定可靠的快速實(shí)用算法，以便開(kāi)展大量實(shí)測(cè)試驗(yàn)數(shù)據(jù)和理論歸納式之間的關(guān)系歸納和規(guī)律總結(jié)。

BP網(wǎng)絡(luò)[16-17]是一種按誤差逆?zhèn)鞑ニ惴ㄓ?xùn)練的多層前饋網(wǎng)絡(luò)，是應(yīng)用較廣泛的神經(jīng)網(wǎng)絡(luò)模型之一，方便提供簡(jiǎn)易算法和結(jié)果。BP網(wǎng)絡(luò)或算法能學(xué)習(xí)和存貯大量輸入-輸出模式映射關(guān)系，無(wú)需事先揭示描述這種映射關(guān)系的數(shù)學(xué)方程，使用最速下降法作為學(xué)習(xí)規(guī)則，通過(guò)反向傳播不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，使網(wǎng)絡(luò)誤差平方和最小，迅速收斂，提高計(jì)算效率。作為一種基本的智能信息處理系統(tǒng)，BP神經(jīng)網(wǎng)絡(luò)法的核心功能及算法思想由基本拓?fù)浣Y(jié)構(gòu)包括輸入層(input)、隱層(hide layer)和輸出層(output layer)實(shí)現(xiàn)，本文嘗試?yán)迷摲椒ㄌ幚肀▽?shí)驗(yàn)數(shù)據(jù)，預(yù)測(cè)爆炸圍巖地層中任一點(diǎn)圍巖介質(zhì)峰值壓力特征區(qū)間，并考察該簡(jiǎn)易方法的實(shí)效性。

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

人工神經(jīng)網(wǎng)絡(luò)模擬人腦結(jié)構(gòu)和功能的信息處理系統(tǒng)，其近似信息處理的有效性取決于層與層間高度的互連結(jié)構(gòu)，每層可由多神經(jīng)元組成，各神經(jīng)元通過(guò)權(quán)重分配連接至下一層，形成的基本模型通過(guò)數(shù)學(xué)計(jì)算工具實(shí)現(xiàn)結(jié)構(gòu)傳遞和層間有機(jī)聯(lián)系。

圖1 神經(jīng)網(wǎng)絡(luò)模型示意

神經(jīng)網(wǎng)絡(luò)模型的基本流程和示意如圖1所示，x1～xn由上一層神經(jīng)元傳來(lái)輸入信號(hào)，wi表示神經(jīng)元之間的連接權(quán)值，θ表示初值為1的閾值，f為傳輸函數(shù)，第j

個(gè)單元的輸出值可由公式1給出：

(1)

利用輸出層的各輸出值誤差，限定給出直接前導(dǎo)層的誤差，并對(duì)權(quán)值和閾值更新，層層反向調(diào)整，直至誤差減小到預(yù)設(shè)值。主要算法結(jié)構(gòu)由輸入層、隱層和輸出層組成，流程主要分為兩部分。

(1) 數(shù)據(jù)正向輸入：權(quán)值和閾值初始化之后，網(wǎng)絡(luò)隨機(jī)選取一組輸入向量(x1，x2,…,xn)和一組輸出向量(y1，y2,…,yn),隱含層的輸出值如公式1所示，則輸出層第k個(gè)單元的輸出值為：

(2)

式中：標(biāo)號(hào)(1)為輸入層到隱含層的傳輸函數(shù)、權(quán)值和閾值，標(biāo)號(hào)(2)為隱含層到輸出層的傳輸函數(shù)、權(quán)值和閾值，m為隱含層單元數(shù)。

(2) 誤差反向傳遞：用輸出層誤差對(duì)前向權(quán)值和閾值進(jìn)行網(wǎng)絡(luò)訓(xùn)練修正。本文選取Levenberg-Marquardt[18]算法作為訓(xùn)練函數(shù)，結(jié)合高斯-牛頓算法和梯度下降法的優(yōu)點(diǎn)，具有較快的收斂速度。e(X)表示每個(gè)輸入樣本X在輸出層的殘差，X為權(quán)值和閾值組成的向量，XN表示第N次迭代，則基本算法可表示為：

XN+1=XN-(JTJ+λE)-1JTe(X)

(3)

式中：E為單位矩陣；λ為正數(shù)，在迭代中不斷調(diào)整，使JTJ+λE保持正定矩陣；J為雅各比矩陣：

(4)

式中:K為X向量維度，l為輸出層單元數(shù)。

完成初始修正后，隨機(jī)選取下一組數(shù)據(jù)繼續(xù)上述過(guò)程，直至輸出誤差達(dá)到預(yù)設(shè)值或趨于穩(wěn)定，直至整個(gè)網(wǎng)絡(luò)收斂。

2 BP神經(jīng)網(wǎng)絡(luò)構(gòu)建

2.1 數(shù)據(jù)集

本文收集130個(gè)地下封閉爆炸實(shí)驗(yàn)數(shù)據(jù)(William&Robert、E.C.Leong、Hiroyoshi Ichino、Paul R.Gefken等)，根據(jù)實(shí)驗(yàn)土質(zhì)，將圍巖數(shù)據(jù)初分為：黏土、砂性土和巖石三組，各組自動(dòng)隨機(jī)5個(gè)數(shù)據(jù)作為測(cè)試樣本，取值在表1中列出。輸入?yún)?shù)包括介質(zhì)參數(shù)：縱波波速c、密度ρ和飽和度S，使用爆炸比例距離D來(lái)描述位置參數(shù)和炸藥參數(shù)，(此處定義D=R/W1/3，R為距離爆炸中心的距離，W為炸藥的TNT當(dāng)量)?？紤]爆炸比例距離和峰值壓力的數(shù)值量級(jí)差別較大，將其對(duì)數(shù)化處理。

表1 測(cè)試樣本

輸入和輸出參數(shù)及其范圍如表2所示。

表2 輸入輸出參數(shù)取值范圍

2.2 網(wǎng)絡(luò)構(gòu)建

輸入和輸出層分別有4和1個(gè)神經(jīng)元，網(wǎng)絡(luò)性能受隱含層神經(jīng)元數(shù)量影響較大，通過(guò)經(jīng)驗(yàn)公式(5)限定范圍，逐步試驗(yàn)得出最佳節(jié)點(diǎn)數(shù)。

(5)

式中：α為1～10之間的常數(shù)。

為試驗(yàn)得到最佳網(wǎng)絡(luò)性能，共建立了30個(gè)BP神經(jīng)網(wǎng)絡(luò)，隱含層神經(jīng)元個(gè)數(shù)分別為3～12。網(wǎng)絡(luò)選用雙曲正切S型函數(shù)tansig作為輸入層到隱含層的傳輸函數(shù)，隱含層到輸出層選用線性函數(shù)purelin。通過(guò)軟件Matlab?函數(shù)newff命令建立網(wǎng)絡(luò)。將測(cè)試樣本輸入網(wǎng)絡(luò)，得到結(jié)果與實(shí)測(cè)值均值的絕對(duì)誤差MAE為：

(6)

得到不同網(wǎng)絡(luò)的MAE在圖2中列出，可以看出，黏土和巖石隱含層節(jié)點(diǎn)數(shù)為6時(shí)，網(wǎng)絡(luò)均值的絕對(duì)誤差最小；砂性土的隱含層節(jié)點(diǎn)數(shù)為7時(shí)，預(yù)測(cè)性能最好。選用各組最佳網(wǎng)絡(luò)預(yù)測(cè)峰值壓力，通過(guò)net.w和net.b命令可查看訓(xùn)練完成的權(quán)值和閾值，得到預(yù)測(cè)表達(dá)式：

ln(P)=[tan sig(W1X+θ1)]W2+θ2

(7)

(8)

式中：P為峰值壓力；X為輸入矩陣[ln(D)，c，ρ，S]T；Wi和θi(i=1,2)為權(quán)值和閾值矩陣，取值在表3中列出。

圖2 不同網(wǎng)絡(luò)的MAE

3 結(jié)果及分析

3.1 現(xiàn)有預(yù)測(cè)方法

(1) 經(jīng)驗(yàn)公式：現(xiàn)有峰值壓力公式的能力表現(xiàn)可通過(guò)對(duì)比美軍常規(guī)武器防護(hù)規(guī)范TM5-855-1的建議方法進(jìn)行。Robert等[12]根據(jù)不同介質(zhì)的實(shí)驗(yàn)數(shù)據(jù)給出擬合公式，Leong等[13]通過(guò)爆炸試驗(yàn)提出了無(wú)量綱化預(yù)測(cè)模型，市野宏嘉等[14]通過(guò)爆炸試驗(yàn)和參數(shù)分析得出含飽和度S和密度ρ的預(yù)測(cè)公式，各具體公式的形式對(duì)比在表4中列出。

(2) 選擇多元回歸分析(Multiple Value Regression Analysis，MVRA)作為建議方法效果的檢驗(yàn)方法。MVRA方法建立多自變量和一因變量之間的數(shù)學(xué)關(guān)系，通過(guò)判斷幾個(gè)特定變量間是否存在相關(guān)關(guān)系，找出函數(shù)間合適的數(shù)學(xué)表達(dá)式。將網(wǎng)絡(luò)輸入?yún)?shù)作為自變量，輸出參數(shù)作為因變量，通過(guò)軟件Origin Pro?得到本算例的多元線性回歸表達(dá)式：

表3 網(wǎng)絡(luò)參數(shù)表

表4 經(jīng)驗(yàn)公式

ln(P)={-2.804 91×D[kg/m1/3]+0.144 49×

c[km/s]+0.693 87×ρ[g/cm3]+

2.354 39×S[104%]-2.659 25}

(9)

3.2 結(jié)果檢驗(yàn)

為驗(yàn)證預(yù)測(cè)方法的可靠性，使用15個(gè)未參與訓(xùn)練過(guò)程的數(shù)據(jù)對(duì)BP神經(jīng)網(wǎng)絡(luò)和多元回歸分析(Multiple Value Regression Analysis，MVRA)預(yù)測(cè)方法的預(yù)測(cè)效果進(jìn)行相互測(cè)試和對(duì)比；考察各方法計(jì)算結(jié)果的參數(shù)效應(yīng)。

圖3列出基于各預(yù)測(cè)方法的預(yù)測(cè)值和實(shí)測(cè)值對(duì)比分布情況，供效果判斷。

對(duì)比各圖分布效果可以看出，BP神經(jīng)網(wǎng)絡(luò)具有較好的擬合效果，對(duì)不同參數(shù)的預(yù)測(cè)偏離度較??；各經(jīng)驗(yàn)公式在某些區(qū)間或節(jié)點(diǎn)上的準(zhǔn)確度較好，但整體誤差和離散程度較大，在變換不同參數(shù)時(shí)的準(zhǔn)確度變化大，預(yù)測(cè)效果出現(xiàn)不穩(wěn)定。

此外在黏土、砂土、巖石不同介質(zhì)各方法預(yù)測(cè)效果的對(duì)比結(jié)果可知，MVRA的擬合性整體優(yōu)于經(jīng)驗(yàn)公式，圖4中6種預(yù)測(cè)模型預(yù)測(cè)值的MAE值對(duì)比可以看到，BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)方法平均誤差最小，在0.5以下，預(yù)測(cè)效果最好，對(duì)砂性土的預(yù)測(cè)最為準(zhǔn)確；各經(jīng)驗(yàn)公式由于在不同參數(shù)區(qū)間的準(zhǔn)確度不一，個(gè)別區(qū)間偏離嚴(yán)重，整體平均誤差比BP神經(jīng)網(wǎng)絡(luò)法及MVRA均大，三種方法對(duì)比而言，在此例巖石類數(shù)據(jù)的預(yù)測(cè)結(jié)果中誤差出現(xiàn)了偏大現(xiàn)象。MVRA整體預(yù)測(cè)誤差能力介于中游，對(duì)于黏土的預(yù)測(cè)誤差最小。

3.3 地下結(jié)構(gòu)防護(hù)設(shè)計(jì)的實(shí)例分析

為進(jìn)一步考察方法的實(shí)用性，開(kāi)展方法的實(shí)例應(yīng)用分析。將以上各方法在具體巖層條件下地下結(jié)構(gòu)爆沖荷載的防護(hù)設(shè)計(jì)中進(jìn)行對(duì)比，進(jìn)一步考察模型的預(yù)測(cè)效果。

假定爆炸圍巖介質(zhì)中地下結(jié)構(gòu)防護(hù)設(shè)計(jì)的爆沖峰值荷載Pr的主要影響因素和條件為：

(1) 圍巖介質(zhì)峰值壓力Pso：主要影響因素為位置參數(shù)和介質(zhì)參數(shù)；

(2) 炸彈等震源耦合系數(shù)f：與爆炸埋深和炸藥量有關(guān)，封閉爆炸耦合系數(shù)取1；

(3) 考慮爆炸波在土-結(jié)構(gòu)交界面多次反射、相互疊加，峰值增加效應(yīng)。

Pr通過(guò)下式簡(jiǎn)化計(jì)算：

Pr=f×N×Pso

(10)

式中：N為放大倍數(shù)或反射系數(shù)。

根據(jù)Hailong等對(duì)埋入式鋼筋混凝土拱結(jié)構(gòu)在爆炸荷載作用下的響應(yīng)問(wèn)題進(jìn)行的試驗(yàn)研究[19]，在不同TNT炸藥量下測(cè)得結(jié)構(gòu)和圍巖介質(zhì)的壓力時(shí)程變化，選取P2和P3位置的測(cè)點(diǎn)數(shù)據(jù)為研究對(duì)象，P2測(cè)得圍巖介質(zhì)壓力，P3測(cè)得拱頂反射荷載值，反射系數(shù)在6 kg和17.4 kg藥量下分別取1.55和2.1，耦合系數(shù)為1。測(cè)點(diǎn)和炸藥位置示意，如圖5所示，試驗(yàn)場(chǎng)地土質(zhì)為風(fēng)化砂質(zhì)黏土，密度ρ為1 810 kg/m3，縱波波速為1 025 m/s，含水率ω為25.2%，取顆粒密度ρs=2 600 kg/m3，水密度ρw=1 000 kg/m3，則，飽和度：

(a) BP神經(jīng)網(wǎng)絡(luò)

(b) TM5-855-1

(c) Leong

(d) Robert

(e) Hiroyoshi

(f) MVRA

圖4 預(yù)測(cè)結(jié)果的MAE

(11)

將以上參數(shù)輸入到BP神經(jīng)網(wǎng)絡(luò)模型，訓(xùn)練完成后和MVRA模型分別得出P2，P3的圍巖介質(zhì)峰值壓力，進(jìn)而得到P3點(diǎn)的爆炸荷載，計(jì)算值和相對(duì)誤差在表5

圖5 測(cè)點(diǎn)和炸藥位置

中列出，對(duì)比得到BP神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)結(jié)果與實(shí)測(cè)值較接近。

3.4 敏感度分析

為考察各參數(shù)對(duì)峰值壓力結(jié)果的影響程度，進(jìn)而開(kāi)展多自變量參數(shù)變化的敏感度層次分析。使用余弦振幅法評(píng)估輸入特征參數(shù)：縱波波速、密度和飽和度等參數(shù)對(duì)結(jié)果的影響程度[20]。每參數(shù)的所有數(shù)據(jù)用空間V表示。

V={v1,v2,v3,…,vn}

(12)

vi={vi1,vi2,vi3,…,vim}

(13)

式中:v為空間V的長(zhǎng)度;為m向量;n為參數(shù)個(gè)數(shù);m為數(shù)據(jù)個(gè)數(shù)。

表5 P2和P3誤差預(yù)測(cè)

每個(gè)數(shù)據(jù)集可被認(rèn)為是m維空間中的點(diǎn)，任意兩個(gè)數(shù)據(jù)集vi和vj間的影響程度可由式(14)得出。

(14)

黏土、砂土、巖石介質(zhì)的特征參數(shù)敏感度對(duì)比分析如圖6。對(duì)比可以看出，三種地層介質(zhì)中縱波波速對(duì)結(jié)果的影響都最為明顯和顯著；黏土和巖石地層中飽和度對(duì)結(jié)果影響最??；砂性土中敏感度最低的特征參數(shù)是密度。

圖6 特征參數(shù)的敏感度分析

4 結(jié) 論

以上基于BP(Back Propagation)神經(jīng)網(wǎng)絡(luò)法，通過(guò)分析130個(gè)地下封閉爆炸實(shí)測(cè)數(shù)據(jù)，區(qū)分黏土、砂性土和巖石三組地層介質(zhì)，給出一種爆炸圍巖介質(zhì)峰值壓力的簡(jiǎn)易預(yù)測(cè)方法和評(píng)價(jià)；通過(guò)分析特征參數(shù)：爆炸比例距離、土體縱波波速、密度和飽和度等產(chǎn)生的影響，將輸出BP神經(jīng)網(wǎng)絡(luò)模型的爆炸圍巖介質(zhì)峰值壓力特征和實(shí)例結(jié)果對(duì)比，認(rèn)為：

(1) BP神經(jīng)網(wǎng)絡(luò)模型能夠快速處理復(fù)雜非線性映射，區(qū)分黏土、砂性土和巖石三組介質(zhì)數(shù)據(jù)，根據(jù)爆炸比例距離、縱波波速、密度和飽和度等特征參數(shù)，可分別建立相應(yīng)網(wǎng)絡(luò)，當(dāng)隱含層單元數(shù)分別為6、7、6時(shí)，輸出的峰值壓力值獲得最佳網(wǎng)絡(luò)性能。

(2) 利用測(cè)試樣本對(duì)BP模型、美軍規(guī)范及類似經(jīng)驗(yàn)公式的效果對(duì)比中，BP神經(jīng)網(wǎng)絡(luò)得到預(yù)測(cè)結(jié)果的平均誤差最小，其三組測(cè)試數(shù)據(jù)的誤差都小于0.5，相比各經(jīng)驗(yàn)公式優(yōu)勢(shì)明顯，預(yù)測(cè)性能也好于MVRA。

(3) 在工程實(shí)例的應(yīng)用中，考慮介質(zhì)反射系數(shù)進(jìn)行風(fēng)化砂質(zhì)黏土地層爆炸圍巖介質(zhì)峰值壓力的預(yù)測(cè)，輸入工程實(shí)例參數(shù)的BP神經(jīng)網(wǎng)絡(luò)和MVRA模型的結(jié)果對(duì)比，BP神經(jīng)網(wǎng)絡(luò)的爆炸荷載預(yù)測(cè)值相對(duì)誤差在20%以內(nèi)，優(yōu)于MVRA。進(jìn)行了三種介質(zhì)地層各特征參數(shù)的敏感度分析，各圍巖介質(zhì)的縱波波速對(duì)其峰值壓力的影響均最為明顯，黏土和巖石地層飽和度的敏感度較小，對(duì)砂性土影響最小的介質(zhì)參數(shù)是密度。

以上基于BP(Back Propagation)神經(jīng)網(wǎng)絡(luò)法,提出一種土中任一點(diǎn)圍巖介質(zhì)峰值壓力的簡(jiǎn)易評(píng)價(jià)和預(yù)測(cè)方法，并通過(guò)與經(jīng)驗(yàn)法對(duì)比，部分驗(yàn)證了方法的實(shí)效性。建議方法可為預(yù)測(cè)地下結(jié)構(gòu)防護(hù)設(shè)計(jì)中爆炸荷載的響應(yīng)峰值提供一種算例參考。