基于逐步回歸算法的邊坡爆破振動(dòng)控制研究

(1.蘭州理工大學(xué) 石油化工學(xué)院，蘭州 730050；2.浙江大學(xué) 海洋學(xué)院，浙江 舟山 316021； 3. 福建省特種設(shè)備檢驗(yàn)研究院，福州 350008)

1 研究背景

爆破振動(dòng)控制是工程爆破研究的主要課題之一，合理的振動(dòng)幅值能改善爆破效果、降低成本、確保周圍建筑物安全、提高工程質(zhì)量[1-2]。國(guó)內(nèi)外學(xué)者對(duì)于爆破振動(dòng)幅值研究大致有4種方法：薩道夫斯基法、數(shù)值法、波動(dòng)法和數(shù)據(jù)擬合法[3]。盧文波等[4]、Tripathy等[5]、李啟月等[6]、張立國(guó)等[7]采用多元回歸法對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行回歸擬合，確定不同影響因素下薩道夫斯基參數(shù)取值；胡建華等[8]利用數(shù)值法的相似理論，推導(dǎo)出修正后的爆破振動(dòng)速率公式，分析不同影響因子和振動(dòng)速率之間的變化關(guān)系以及因子之間的相關(guān)性；Dhakal等[9]亦采用文獻(xiàn)[8]的方法研究了爆破振動(dòng)引起的結(jié)構(gòu)響應(yīng)特性；李夕兵等[10]、楊建華等[11]、謝全民等[12]采用波動(dòng)法對(duì)爆破振動(dòng)速率進(jìn)行研究，將垂直方向爆破速率峰值作為振動(dòng)大小的判據(jù)；徐全軍等[13]、文建華等[14]采用遺傳算法或神經(jīng)網(wǎng)絡(luò)算法對(duì)現(xiàn)場(chǎng)實(shí)測(cè)爆破振動(dòng)數(shù)據(jù)進(jìn)行擬合分析，使得實(shí)測(cè)值與預(yù)測(cè)值之間存在較高擬合度，該方法可較準(zhǔn)確地給出適合具體施工現(xiàn)場(chǎng)的振動(dòng)峰值速度模型，但未考慮變量因子之間的相互影響。

本文采用薩道夫斯基和數(shù)據(jù)擬合相結(jié)合的方法，實(shí)施修正后的多元回歸擬合。振動(dòng)速率回歸分析前將振動(dòng)速率影響因子函數(shù)化，對(duì)因子函數(shù)逐步回歸分析，確保最終入選因子的優(yōu)越性；再將實(shí)測(cè)振動(dòng)數(shù)據(jù)進(jìn)行逐步回歸，剔除不顯著因子，入選顯著因子，確定道夫斯基公式參數(shù)。該算法優(yōu)點(diǎn)在于不僅可以分析自變量和因變量之間的影響關(guān)系，還可以分析自變量之間的相互影響程度，保證了因子的全面性、完整性和初始性，并且代入運(yùn)算的因子都經(jīng)過了循環(huán)篩選，顯著性較高；所以修正后的多元回歸擬合用于爆破振動(dòng)研究更具有優(yōu)越性。呂濤等[3]在對(duì)薩道夫斯基公式研究過程中忽略高程和主頻的影響，而通過相關(guān)研究[15-17]發(fā)現(xiàn)，高程和主頻大范圍的變化會(huì)引起振動(dòng)速率的發(fā)散。故對(duì)于爆破振動(dòng)幅值研究，高程差和主頻可視為主要的研究物理量。

2 逐步回歸分析

逐步回歸分析法主旨為自變量依次進(jìn)入模型，每一個(gè)自變量因子進(jìn)入模型都進(jìn)行F檢驗(yàn)，同時(shí)將已進(jìn)入的自變量因子進(jìn)行t檢驗(yàn)；后引入自變量使前一引入自變量變得不再顯著時(shí)，則將前一自變量剔除，每次循環(huán)后保留的自變量都對(duì)因變量有顯著影響[18]。

2.1 基本原理

逐步回歸類同于多元回歸，假設(shè)有n組自變量和對(duì)應(yīng)的應(yīng)變量為

(1)

式中：xij表示第i組第j個(gè)自變量；yi表示第i組自變量對(duì)應(yīng)的因變量。

多元回歸方程求解過程中，可將多元回歸方程的法方程改寫為

(2)

式中：aij為法方程自變量；bj為未知量參數(shù)；Sj為法方程因變量；cij為變換參數(shù)。

(3)

直到m步(1≤m≤k)，等價(jià)于方程兩端各乘以一轉(zhuǎn)換矩陣，即Dk[aij]{bj}=Dk[cij]{sjy}；故此，向量[aij]的第2列被單位矩陣[δij]相應(yīng)置換，其他元素也進(jìn)行了相應(yīng)變換。同時(shí)，向量[Cij]相應(yīng)單位列向量被另一個(gè)列向量(式(4))所置換，其余各元素向量也進(jìn)行了等價(jià)變換。

(4)

隨著循環(huán)計(jì)算的進(jìn)行，部分預(yù)報(bào)因子xm被選入回歸方程，同時(shí)單位向量已經(jīng)置換成相應(yīng)[aij]中被選入因子xm對(duì)應(yīng)的各列，新向量也引入了[cij]對(duì)應(yīng)的各列中。

2.2 基本步驟

2.2.1 計(jì)算相關(guān)矩陣

為保證計(jì)算的精度與準(zhǔn)確性，將一次均值算法換為二次均值算法，用標(biāo)準(zhǔn)化的相關(guān)矩陣[rij]代替[sij]，擴(kuò)展為(k+1)階矩陣，y用n表示，即

(5)

其中：

式中：xi為研究對(duì)象的自變量；yi為研究對(duì)象的因變量，故有

(7)

2.2.2 因子的篩選和消元變換

由以上循環(huán)計(jì)算后，第m步回歸循環(huán)時(shí)回歸方程各引入和剔除變量所用特征值可表示如下。

(9)

其中，

(2)剩余平方和Q(m)，其表達(dá)式為

(11)

(12)

(4)引入因子xkm+1的檢驗(yàn)。在第m步循環(huán)時(shí)，偏回歸平方和最大的因子對(duì)應(yīng)變量y影響最為顯著，故選入回歸方程，即

(13)

2.2.3 回歸方程

進(jìn)行第m步循環(huán)時(shí)，回歸方程為

(14)

圖1 試驗(yàn)現(xiàn)場(chǎng)及測(cè)振儀安裝Fig.1 Blasting circumstance and vibration instrument

3 工程實(shí)例

某施工平臺(tái)臨近高陡邊坡，緊接最終邊坡的交接位置，巖層較為破碎，因前期爆破振動(dòng)影響，臺(tái)階穩(wěn)定性較差。為了保證安全生產(chǎn)、確保最終邊坡的穩(wěn)固以及提高生產(chǎn)效益和經(jīng)濟(jì)利益，現(xiàn)對(duì)爆破振動(dòng)信號(hào)幅值傳播規(guī)律進(jìn)行研究。起爆藥包為海西東諾產(chǎn)Φ32 mm乳化炸藥，干孔炸藥采用海西東諾產(chǎn)巖石膨化硝銨炸藥，水孔采用Φ90 mm水膠炸藥；5次不同時(shí)期爆破施工過程對(duì)爆區(qū)周圍地表質(zhì)點(diǎn)采用爆破測(cè)振儀進(jìn)行測(cè)振試驗(yàn)，每次試驗(yàn)安裝4個(gè)臨近監(jiān)測(cè)點(diǎn)，共計(jì)20組試驗(yàn)數(shù)據(jù)；但因?yàn)楦鞣N操作和接收原因，部分測(cè)振結(jié)果數(shù)據(jù)雜散，無(wú)法使用，故只有18組數(shù)據(jù)可用于研究分析。試驗(yàn)現(xiàn)場(chǎng)及測(cè)振儀安裝見圖1。

4 基于逐步回歸分析的振動(dòng)速率監(jiān)控

4.1 炸藥量逐步回歸分析

1978年蘭格福爾斯在《巖石爆破現(xiàn)代技術(shù)》中給出了有關(guān)爆破藥量計(jì)算的半經(jīng)驗(yàn)公式[19]，即

Q=k2W2+k3W3+k4W4。

(15)

式中：Q為爆破藥量；W表示最小抵抗線；W2表示與介質(zhì)破裂面積相關(guān)的參量；W3表示與介質(zhì)體積相關(guān)的參量；W4表示與最小抵抗線相關(guān)的參量，又稱之為“膨脹”或“拋擲”分量；k2，k3，k4均為系數(shù)。

此次爆破測(cè)振試驗(yàn)炸藥量以及相關(guān)參數(shù)取值如表1所示。

表1 炸藥量相關(guān)參數(shù)Table 1 Parameters of explosive charge

注:0510-1表示5月10日測(cè)量時(shí)編號(hào)為1的傳感器，以此類推

單孔負(fù)擔(dān)面積、單孔爆破方量、最小抵抗線等諸多因素直接影響炸藥量的取值，并且這些因素之間也存在著不同程度的影響，要在一個(gè)函數(shù)模型里考慮到所有影響因素不太現(xiàn)實(shí)。因此，采用灰色關(guān)聯(lián)理論建立相關(guān)函數(shù)模型分析各因子之間的關(guān)聯(lián)程度。本文采用逐步回歸的理念，入選顯著影響因子，剔除不顯著影響因子，通過逐步回歸循環(huán)達(dá)到最優(yōu)取值的目的。

建立逐步回歸模型，程序運(yùn)行后得到：k2=0.057，k3=0.359，k4=0.084,常值為-0.017；依據(jù)逐步回歸入選與剔除原理，得到選入因子為單孔爆破方量，剔除因子為單孔負(fù)擔(dān)面積和最小抵抗線。對(duì)回歸報(bào)告進(jìn)一步分析，可知回歸方程中R2=0.940，F(xiàn)=248.816；方程確定性系數(shù)R2表示自變量對(duì)因變量的相關(guān)性，其取值在[0，1]之間，越靠近1，表明相關(guān)性越強(qiáng)；回歸方程的線性關(guān)系顯著性由F的顯著性來(lái)檢驗(yàn)，查F(3,14)分布表得λ=6.68，F(xiàn)=248.816>λ，由于Sig=P(F>λ)=0<0.05(置信水平)，表明該回歸方程相關(guān)性特別顯著。

對(duì)炸藥量函數(shù)進(jìn)行逐步回歸分析時(shí)，因變量(炸藥量)與自變量(單孔負(fù)擔(dān)面積、單孔爆破方量)之間的數(shù)據(jù)分布如圖2所示，該圖可以直觀顯示除部分野值點(diǎn)外的相對(duì)合理準(zhǔn)確數(shù)據(jù)的變化趨勢(shì)。

圖2 炸藥量與單孔負(fù)擔(dān)面積、單孔爆破方量之間的關(guān)系曲線Fig.2 Curves of explosive charge vs. burden area ofsingle-hole and unit volume of single-hole blasting

經(jīng)逐步回歸分析后炸藥量函數(shù)關(guān)系可表示為

(16)

工程爆破實(shí)例中，一般最小抵抗線W取值在0.5～10 m范圍內(nèi)。故在W=0.5～10 m范圍，按照0.5 m的步長(zhǎng)取20個(gè)值代入上式進(jìn)行運(yùn)算。最終結(jié)果表明：當(dāng)W>2 m時(shí)，介質(zhì)破裂面積相關(guān)參量(公式第1項(xiàng))與總藥量的比值<10%；當(dāng)W<10 m時(shí)，最小抵抗線相關(guān)分量(公式第3項(xiàng))與總藥量的比值<10%；按照《爆破安全規(guī)程》規(guī)定，爆破允許誤差為±10%以內(nèi)，則當(dāng)2 m≤W≤10 m時(shí)爆破炸藥量主要由爆破方量所決定，此時(shí)炸藥量公式可改寫為

Q=k3W3。

(17)

4.2 振動(dòng)速率逐步回歸分析

在該爆區(qū)進(jìn)行5次爆破測(cè)振試驗(yàn)，共采集18個(gè)點(diǎn)次的振動(dòng)速率相關(guān)數(shù)據(jù)，具體如表2所示。采集0510-1所測(cè)的垂向振動(dòng)速率傳播規(guī)律，如圖3所示。

國(guó)內(nèi)外對(duì)于爆破振動(dòng)峰值速率模型運(yùn)用最為廣泛的是薩道夫斯基公式[20-21]，具體數(shù)學(xué)模型如式(18)所示。

(18)

式中：V為質(zhì)點(diǎn)垂直方向的最大振動(dòng)速率；R為水平距離；k，α均為與爆區(qū)環(huán)境中水文地質(zhì)、地形條件以及介質(zhì)參數(shù)等相關(guān)系數(shù)。

表2 試驗(yàn)爆區(qū)振動(dòng)數(shù)據(jù)Table 2 Vibration data of blasting area

注：ρ為比藥量

圖3 垂向振動(dòng)速率-時(shí)間分布規(guī)律Fig.3 Distribution rule of vibration velocityin Y axis against time

實(shí)踐表明，高程對(duì)于振動(dòng)速率具有重要影響，當(dāng)高程超過一定值時(shí)爆破振動(dòng)信號(hào)會(huì)出現(xiàn)高程放大效應(yīng)[22-23]。為研究高程對(duì)爆破振動(dòng)信號(hào)影響規(guī)律，本文對(duì)含高程參量的薩道夫斯基公式進(jìn)行逐步回歸分析，具體函數(shù)方程為

(19)

式中β為與爆區(qū)環(huán)境中水文地質(zhì)、地形條件以及介質(zhì)參數(shù)等相關(guān)的系數(shù)。

lgV=lgk+αlgρ+βlgH。

(20)

建立逐步回歸模型，程序運(yùn)行后得到：k=200.8，α=1.8，β=-0.563；對(duì)回歸報(bào)告進(jìn)一步分析，可知回歸方程中方程確定性系數(shù)R2=0.973，基本接近1，表明自變量對(duì)因變量有較強(qiáng)的相關(guān)性，該回歸方程相關(guān)性特別顯著。經(jīng)逐步回歸后可得振動(dòng)速率函數(shù)關(guān)系為

(21)

4.3 高程放大影響

邊坡爆破振動(dòng)速度的高程放大效應(yīng)是在一定條件下產(chǎn)生的，受爆破振動(dòng)荷載特性及邊坡坡形和坡高等因素的影響。爆破振動(dòng)荷載作用下，邊坡坡面不同高程臺(tái)階巖體結(jié)構(gòu)的自振主頻處于爆破振動(dòng)荷載主頻帶范圍內(nèi)，臺(tái)階部位巖體結(jié)構(gòu)的振動(dòng)響應(yīng)會(huì)產(chǎn)生“鞭梢效應(yīng)”。大量工程實(shí)例驗(yàn)證表明，運(yùn)用薩道夫斯基公式預(yù)測(cè)地表質(zhì)點(diǎn)峰值振動(dòng)速度時(shí)，只有在同一高程的平整地形下具有較高精度，如存在較大高程差以及地表起伏明顯的情況下，實(shí)測(cè)數(shù)據(jù)與預(yù)測(cè)數(shù)據(jù)之間存在較大誤差。故國(guó)內(nèi)外部分學(xué)者[22-23]對(duì)高程影響振動(dòng)峰值速度進(jìn)行研究，并指出超出一定范圍時(shí)隨高差正向增大高程放大率也出現(xiàn)增大效應(yīng)，反之降低。

圖4 垂向振動(dòng)速率幅值-高程差變化曲線Fig.4 Variation of vibration amplitude in Y axisagainst elevation difference

垂向振動(dòng)速率幅值-高程差變化曲線如圖4所示。由爆區(qū)垂向振動(dòng)速率逐步回歸分析報(bào)告和圖4可知：高程差作為入選因子進(jìn)入回歸方程，而比藥量在循環(huán)過程中逐步被剔除。除去部分野值點(diǎn)外，越靠近炮孔，振動(dòng)速度越大；隨監(jiān)測(cè)點(diǎn)離炮孔距離增大，地表質(zhì)點(diǎn)振動(dòng)速率呈遞減趨勢(shì)。

根據(jù)此次爆破測(cè)振試驗(yàn)數(shù)據(jù)可確定與該爆區(qū)環(huán)境中水文地質(zhì)、地形條件以及介質(zhì)參數(shù)等相關(guān)系數(shù)k，α，β的取值，且已知每次起爆的最大段藥量Q；將這些參數(shù)代入式(17)可得到不考慮高程影響的爆破振動(dòng)峰值速度V。

假設(shè)高程放大率為ε，則高程放大率公式可表示為

ε=(Hv-Pv)/Pv。

(22)

式中：Hv為考慮高程影響時(shí)的最大垂向速度；Pv為不考慮高程影響時(shí)的最大垂向速度。

通過對(duì)高程放大率和高程差回歸分析可確定邊坡爆破高程放大效應(yīng)的臨界值。當(dāng)ε>0時(shí)，會(huì)出現(xiàn)爆破振動(dòng)信號(hào)放大效應(yīng)；針對(duì)本文所研究的施工環(huán)境、地質(zhì)條件和爆破參數(shù)，對(duì)放大率和高程差回歸分析可得ε=1.693H-63.7；故依據(jù)理論分析給出的高程放大率公式可確定該施工條件下當(dāng)高程差超過38 m時(shí)，出現(xiàn)明顯的振動(dòng)信號(hào)放大效應(yīng)。

圖5為高程差-高程放大率變化曲線，通過分析圖5可以看出：拋去部分野值點(diǎn)外，在本文監(jiān)測(cè)范圍內(nèi)，當(dāng)高程差超過38 m時(shí)，隨高程差值增大，高程放大率也呈現(xiàn)出顯著增強(qiáng)趨勢(shì)。

圖5 高程差-高程放大率變化曲線Fig.5 Curve of elevation difference vs. elevationamplification

4.4 主頻預(yù)測(cè)公式的回歸分析

此次測(cè)振試驗(yàn)中采集0510-1時(shí)所測(cè)得的垂向振動(dòng)速率幅值-主頻變化規(guī)律，如圖6所示。

圖6 垂向振動(dòng)速率幅值-主頻變化規(guī)律Fig.6 Distribution rule of dominant vibrationfrequency vs. vibration velocity in Y axis

在爆破振動(dòng)頻率采集過程中，通常選取最大振幅所對(duì)應(yīng)的頻率(主頻f)為研究對(duì)象，原因是爆破振動(dòng)頻率特性很大程度是由主頻所決定的。由爆破振動(dòng)頻率相似準(zhǔn)數(shù)方程可預(yù)測(cè)工程爆破中任意一介質(zhì)點(diǎn)的振動(dòng)主頻函數(shù)模型為[4]

(23)

取測(cè)振試驗(yàn)采樣的18個(gè)測(cè)點(diǎn)y軸方向(垂直)的主頻-速度相關(guān)數(shù)據(jù)，對(duì)主頻函數(shù)模型進(jìn)行逐步回歸計(jì)算，求得k=231.4，α=-0.76，R2=0.940，F(xiàn)=456.19，Sig=0.03；方程的確定性系數(shù)R2基本接近1，線性關(guān)系顯著性F=456.19>λ，Sig=P(F>λ)=0.03<0.05(置信水平)，說(shuō)明此次逐步回歸運(yùn)算中，主頻回歸方程特別顯著，故主頻預(yù)測(cè)函數(shù)關(guān)系式可表示為

(24)

主頻與比藥量的線性回歸結(jié)果和主頻與Q，R的關(guān)系曲線分別如圖7、圖8所示。分析圖7、圖8可知：比藥量ρ為定值時(shí)，隨水平距離R的增大，主頻呈遞減趨勢(shì)；炸藥量Q不變時(shí)，隨水平距離R的增大，主頻呈遞減趨勢(shì)，并且遞減速率逐漸變緩；水平距離R不變時(shí)，隨炸藥量Q的增大，主頻呈遞減趨勢(shì)，并且遞減速率逐漸變緩。

圖7 主頻與比藥量的線性回歸結(jié)果Fig.7 Linear regression relationship between dominantfrequency and specific charge

圖8 主頻與Q，R的關(guān)系曲線Fig.8 Curves of dominant frequency against Q and R

5 結(jié) 論

(1)對(duì)于爆破振動(dòng)幅值這種多預(yù)報(bào)因子影響的函數(shù)研究，可采用改進(jìn)的多元回歸分析模型，即轉(zhuǎn)換為逐步回歸分析，這樣就可以剔除不顯著預(yù)報(bào)因子，避免了全部因子進(jìn)入回歸方程，導(dǎo)致無(wú)法計(jì)算或者降低了計(jì)算精度。

(2)通過對(duì)炸藥量函數(shù)逐步回歸分析發(fā)現(xiàn)：當(dāng)最小抵抗線W>2.0 m時(shí)，炸藥量主要由單孔爆破方量決定；隨單孔藥量的增加，單孔負(fù)擔(dān)面積、單孔爆破方量都呈線性遞增趨勢(shì)。

(3)對(duì)高程差與高程放大率函數(shù)回歸分析發(fā)現(xiàn)：當(dāng)高程差超過一定范圍時(shí)，爆破振動(dòng)信號(hào)才會(huì)出現(xiàn)高程放大效應(yīng)；由高程放大率函數(shù)可求得針對(duì)此次爆破施工條件出現(xiàn)高程放大效應(yīng)的臨界差值為38 m左右；隨高程差值繼續(xù)增大，振動(dòng)信號(hào)放大率也逐步增強(qiáng)，在本文監(jiān)測(cè)范圍內(nèi)高程差取值越大，振動(dòng)信號(hào)放大作用越強(qiáng)。

(4)主頻函數(shù)逐步回歸分析結(jié)果表明：比藥量ρ為定值時(shí)，隨水平距離R增大，主頻呈遞減趨勢(shì)；炸藥量Q為定值時(shí)，隨水平距離R增大，主頻呈遞減趨勢(shì)，并且遞減速率逐漸變緩；水平距離R為定值時(shí)，隨炸藥量Q的增大，主頻呈遞減趨勢(shì)，并且遞減速率逐漸變緩。