楊招偉，盧文波，高啟棟，陳 明，胡浩然，嚴 鵬，王高輝

(1.武漢大學水資源與水電工程科學國家重點實驗室,湖北 武漢 430072;2.武漢大學水工巖石力學教育部重點實驗室,湖北 武漢 430072)

近年來，我國正在興建的一批大型水電工程都處在西南、西北地區(qū)的崇山峻嶺之中，均需進行高陡邊坡和地下洞室的爆破開挖，由此導致一系列的巖石工程安全問題，特別是在爆炸荷載、地震荷載等動荷載作用下的變形和穩(wěn)定問題[1-4]。因此，對巖石各種動力學參數的研究就顯得更加重要，如何準確地確定現場巖石動力學參數是工程施工過程中需著重考慮的一個方面?，F階段工程中常用的巖石參數獲取方法有室內實驗法和原位實驗法，但由于室內實驗針對小尺度巖樣進行研究，因而得到的也是小尺度巖樣的動力學參數，由于巖體具有較為明顯的尺度效應，為了將室內實驗測試的小尺度巖樣動力學參數的研究過渡到現場大尺度巖體，因此，現場原位實驗具有不可替代的作用，然而原位實驗也有周期長，費用高等缺點[5]。而通過測定爆破地震波(P、S波)在巖石中的傳播速度反演巖石動力學參數的方法較為簡單快捷。如何精確的判別P、S波震相初至時刻對巖石動力學參數反演工作的效率和精度有著不可忽視的影響，因此根據實測振動波形進行S波的初至時刻識別值得開展深入研究。

迄今為止，水利工程中的爆破地震波的識別研究的方法和技術直接來自于對天然地震活動的研究。對于天然地震，針對P、S波初至時刻判別方面進行了大量的研究。Baer等[6]指出P波初至時刻的自動拾取算法主要是基于監(jiān)測振動信號和噪聲的不同特征作為P波震相識別的判據。Stevenson等[7]首次提出利用長短時窗平均比來自動拾取P波震相初至時刻。如今發(fā)展出的P波到時自動識別方法主要有：自回歸模型、AIC法、波形變化增長法、偏振分析法及長短時間平均方法等[8-10]。但在S波自動識別方面，由于P波尾波以及各種反射波相的相互干擾，以上的方法用于S波震相自動識別的效果較差。Diehl等[11]提出了用于地震層析成像的S波自動識別方法，該方法需要對振動信號做STA/LTA判定、偏振分析等處理，計算量較大，比較適合做后續(xù)分析；Baillard等[12]利用峰態(tài)特征函數和特征向量的分解拾取三分像地震信號的P、S波震相。劉希強等[13-14]提出利用振動信號在不同尺度的小波變換中偏振信息的不同來確定S波震相的小波變換方法；何先龍等[15]提出一種二次自回歸模型拾取P、S波的初至時間；張楚旋等[16]采用短時平均過零率、短時能量等4個指標作為S波自動拾取的指標。

雖然學者們在天然地震領域上S波的拾取上取得了豐富的研究成果，但是與天然地震信號相比，水利工程爆破開挖產生的微震信號持續(xù)時間較短、背景噪音復雜、分布范圍較小，受結構面等地質的不均勻性影響很大，導致P、S波分離不明顯，致使適用于地震領域的S波識別算法在水利工程尺度下識別效果不明顯。因此，對水利工程爆破開挖產生的微震信號的S波自動識別進行更深入的研究具有較為重大的意義。

本文中旨在研究根據實測爆破開挖產生的振動波形進行S波震相的快速識別，在滿足安全生產精度要求的前提下，秉承信號處理上“濾波越少，算法越好”的原則[17]，對原始振動數據不進行任何的濾波處理?；谄渌椒ㄒ呀浀玫降腜波震相到時，在此基礎上，對前人識別方法進行改進，根據P波、S波在頻率、偏振方向和水平能量與總能量上的差異對S波進行精確識別。對于爆破地震波中P波震相的初至時刻識別，推薦采用Baer等[6]提出的自動識別算法。

1 識別方法及特征函數

水利水電工程開挖爆破過程中，炸藥能量的很大一部分以地震波的形式向四周傳播，引起地面振動。爆破地震波在傳播過程中可被分為體波和面波，其中根據傳播過程中質點運動和傳播方向關系，可以將體波進一步分為P波和S波，相較于P波，S波具有傳播速度僅次于P波，并且其周期較長、振幅較大的特點。Diehl等[11]指出，振動信號監(jiān)測儀器會最先監(jiān)測到P波到達引起的振動信號，而相對于較為容易的識別出P波初至震相時刻，在S波震相識別方面，由于P波尾波信號以及各種轉換波信號，如PS波的強干擾導致S波淹沒在P波中難以識別。可根據P、S波周期及偏振特性的差異，對與周期差異相關的短時平均過零率及偏振特性差異相關的偏轉角、偏振度及橫向能量與總能量差異等參數進行分析，現分別討論。

1.1 確定計算窗長

在對上述識別參數進行計算時，當計算窗長選取過長，上述特征參數的值變化不明顯；計算窗長選取過短，則導致識別參數值的不穩(wěn)定。Cichowicz等[18]提出根據P波的卓越頻率確定時窗長度的方法。計算時窗長度與P波卓越頻率成反比關系，具體表達式如下：

(1)

式中：l為計算窗長，Δt為儀器采樣間隔時間，fp為P波的卓越頻率。

Boore等[19]提出P波的卓越頻率與振動波形的速度功率譜密度和相應的位移功率譜密度存在如下式所列關系：

(2)

1.2 短時平均過零率

由于P波的周期小于S波周期，在監(jiān)測信號上就表現出如圖1所示的波形前后疏密不同的變化，為了更好地識別S波震相初至時刻，需對這種疏密變化進行量化。對于平穩(wěn)信號來說，可以直接用周期或頻率直接衡量。而對于地震信號、爆破產生的振動信號、噪聲信號等非平穩(wěn)信號而言，由于其周期的不穩(wěn)定，則用短時平均過零率來進行度量。

振動信號單位時間內過零的次數可以定義為短時平均過零率K1，其值可以作為S波初至震相時刻拾取的指標。其定義式為：

(3)

由式(3)可知，短時平均過零率K1的值與選取的時窗長短有關。關于S波的初至時刻震相，因此主要研究S波到達前后振動信號周期變化的差異。將P波初至時刻到S波初至時刻這一區(qū)間稱為S-P區(qū)間，將S波震相初至時刻后(ns-np)Δt時窗內的區(qū)間，稱為S′-S區(qū)間。定義公式為：

(4)

(5)

式中：K2為S-P區(qū)間內的短時平均過零率，K3為S′-S區(qū)間內的短時平均過零率，ns′=ns+(ns-np)，np、ns分別為P波初至時間和S波初至時間對應的信號點序號。

根據張楚旋等[16]對建立的微震數據庫的信號研究可知，短時平均過零率與振動信號的S波震相的初至時刻存在著對應關系，定義S波初至時刻震相前后的短時平均過零率之比C1=K2/K3，計算發(fā)現，C1的取值范圍大于1.2。

1.3 偏轉角

確定P波初至到時后，結合Kanasewich[20]提出的偏振特性概念，選取P波到時后l長度的計算其協(xié)方差矩陣最大特征值所對應的特征向量，即P波的偏振方向。計算式如下：

(6)

P波初至時刻點為計算起點，取時窗長度為l個采樣點往后開始進行逐點進行滑動，求得各個計算窗口最大特征值所對應的特征向量，并計算其與P波偏振方向的夾角即偏轉角α。對偏轉角α進行歸一化處理得：

(7)

1.4 偏振度

P、S波具有較高的線性偏振度，而P波尾波則更多的表現為橢圓偏振的特性。Cichowicz等[18]提出直接根據實測地震振動數據求取偏振度的方法。計算公式為：

(8)

式中：λ1、λ2、λ3為各個數據窗對應的協(xié)方差矩陣的特征值。

1.5 橫向能量與總能量比值

確定P波偏振方向L后，結合地震波的偏振特性和傳播方向間的關系把原有X(水平徑向)、Y(水平切向)、Z(豎直向)坐標系旋轉為L(P波偏振方向)及Q、T(S波偏正方向)坐標系，轉換式為：

(9)

式中：uij(j=1,2,3)為第i個主方向與坐標軸X、Y、Z夾角的余弦值。

水平能量與總能量的比值，即計算窗長內2個橫向方向能量與3個方向總能量的比值：

(10)

通過對實測波形的偏振轉換，參數F4受震動波形中噪聲影響較小，因此參數F4可以用來提高S波震相識別的精度。P波初至時刻、P波尾波及S波初至時刻所對應的識別參數F2、F3、F4的理論統(tǒng)計平均值如表1所示。

對于實測振動數據而言，由于受P波尾波及P-S

表1 統(tǒng)計平均值Table 1 Statistical averages

轉換波的影響，特征參數F2、F3、F4的統(tǒng)計平均值難以精確計算。因此對于S波初至時刻判定，選取特征參數偏轉角F2、偏振度F3和橫向能量與總能量比值F4三者的平方積作為其判定的特征函數，即

(11)

特征函數C2對識別參數F2、F3、F4進行平方處理，擴大了S波與P波尾波及轉換波對應特征值的數值差異。由表1可知，識別特征函數C2值為1時所對應的時刻為S波震相初至到時。但一般來說，由于受P波尾波以及P波識別精度等因素影響，S波初至時刻對應的特征函數C2值往往不等于1。為了能夠定量識別S波初至時刻，將短時橫向能量值作為權重系數，S波初至時刻加權識別特征函數為:

(12)

S波初至時刻會使特征函數C2與橫向能量Q2(t)+T2(t)急劇增大，因此加權識別特征函數值最大對應時刻即為S波震相初至時刻。對加權識別特征函數進行歸一化處理得：

(13)

1.6 S波初至時刻識別流程

按式(4)、(5)分別計算P-S區(qū)間(1～l)、S-S′區(qū)間(l～l+(l-1))內的短時平均過零率，然后再按照C1=K2/K3計算特征函數C1值。對S波初至時刻進行初步判定，縮小加權識別特征函數Cw的計算空間，提升此方法計算效率。根據式(7)、(8)、(10)分別計算出識別參數F2、F3、F4的值并按式(12)計算加權特征函數Cw的值，最終確定S波震相初至時間。

2 數值模擬及對比分析

為了進一步驗證上述S波識別方法的正確性,基于LS-DYNA動力有限元軟件，模擬半無限巖體中柱狀藥包爆破誘發(fā)的振動傳播規(guī)律。計算中巖石物理力學參數及拉壓損傷模型參數取值見表2[21]。因模型具有對稱性，采取1/4模型模擬巖體中單孔爆破特性。模型尺寸為10 m×9 m×20 m，單元240 304個，節(jié)點260 012個，由于炮孔直徑較小，網格劃分采用從炮孔到邊界漸變方式剖分，炮孔附近單元最小尺寸為0.003 6 m，邊界最大尺寸0.8 m，巖體采用solid 164單元模擬。炮孔直徑90 mm，裝藥直徑70 mm，孔深3 m，堵塞0.9 m，在模型頂面布置振動監(jiān)測點，有限元計算模型如圖2所示。

圖3分別給出了爆心距為10、15 m處三向爆破振動速度時程變化曲線。由表2所給出的巖石動力學參數，根據理論公式:

(14)

式中：E為巖體動彈性模量，kPa；ν為巖體動泊松比，ρ為巖體密度，g/cm3。

可以計算出P波的傳播速度為4 248 m/s，S波的傳播速度為2 545 m/s。由此可以估算出不同爆心距S波的準確初至時間，如表3所示。應用上述S波到時識別方法對以上測點進行S波初至時刻識別，將識別結果與理論到時進行對比分析。根據表3結果，采用該方法拾取與理論初至時間相比，二者拾取S初至時刻相差不大，拾取誤差小于3%，可見該方法具有較高的識別率，可以滿足工程尺度的實際要求。

表3 實驗數據與理論模型結果對比Table 3 Comparison between experimentaland theoretical data

3 工程應用研究

3.1 工程概況

豐寧抽水蓄能電站位于中國河北省豐寧滿族自治縣境內，電站工程分兩期開發(fā)建設,一期工程和二期工程裝機容量分別為1 800 MW，裝機6臺，單機容量300 MW，在電網系統(tǒng)中承擔調峰、調頻、調相和事故備用任務。利用二期工程地下廠房地質勘探洞開挖的時機和條件，進行豎直鉆孔爆破實驗。實測資料表明該地質勘探洞以花崗巖為主，總體質量較好，因此可以認為在本次實驗范圍內巖體為均質、各向同性的彈性體，實驗過程中炮孔及振動監(jiān)測點布置如圖4所示，采用半秒微差延期爆破，逐孔起爆。爆破設計參數見表4，表4中D為孔徑，h為孔深，d為藥卷直徑，L裝藥長度，Lb堵塞長度，m單響藥量。利用振動監(jiān)測儀器得到了各個測點的振動時程曲線，利用上述方法分析實測振動波形圖對S波的初至時刻進行精確識別。

基于《爆破安全規(guī)程》、《水電水利工程爆破安全監(jiān)測規(guī)程》等規(guī)范標準開展爆破振動測試，爆破振動測試系統(tǒng)由三向速度檢波器、信號采集與記錄設備、數據處理系統(tǒng)3個部分組成。主要測試設備為TC-4850爆破振動智能監(jiān)測儀。在實驗過程中，TC-4850爆破振動智能監(jiān)測儀用于監(jiān)測爆破振動速度，其采樣頻率為8 000 Hz。結合爆破實驗進行爆破振動監(jiān)測，在探洞中布置爆破振動測試線，如圖4，每個測點相對于爆源具有水平徑向、水平切向及垂直向3個方向，以獲得水平鉆孔爆破不同工況的爆破振動的傳播規(guī)律。實驗中，TC-4850儀器布置在圖4所示的1#、2#、4#、6#、8#及10#測點位置。

炮孔編號起爆位置D/mmh/cmd/mmL/cmLb/cmm/kgⅠ-1上、底部768005060020012.0Ⅰ-2底部768005060020012.0Ⅱ-1中部76600504201808.4Ⅱ-2底部76600504201808.4Ⅲ-1中部76450502701805.4Ⅲ-2底部76450502701805.4

注：?50 mm藥卷由2節(jié)?32 mm炸藥捆綁而成。

3.2 S波初至時刻識別效果

應用上述S波震相加權識別特征函數結合前文所述的識別流程對豐寧二期工程地下廠房探洞開挖實測振動數據進行S波初至時刻識別，識別結果見表5。限于篇幅，只給出典型測點(10#測點)實測振動波形S波初至時刻的識別效果，如圖5所示。

表5 實驗S波震相初至時刻識別結果Table 5 S-wave arrival identification in experiment

從以上數值模擬分析及豐寧水電站實測波形的應用效果可知，本文中所述方法對S波震相初至時刻的識別率較高，且具有簡單、快捷等優(yōu)點，可以滿足實際水利工程的要求。由于該方法是建立在P波已經被識別的基礎上的，下面對其識別參數做進一步論述。

(1)從操作簡單快速來看，根據周期差異作為識別函數C1對S波進行初步識別速率最快，但由于P波尾波及轉換波對頻率的影響導致該特征函數的識別準確率較低，可用于縮小加權識別特征函數Cw計算空間，提升計算效率。

(2)根據橫向能量與總能量比值識別參數對S波震相的識別效果最優(yōu)。造成該現象的主要原因是因為經過坐標系變換，把P波旋轉到垂直分量(L軸)，S波分別旋轉到橫向坐標軸Q、T方向，能夠有效壓制P波尾波及其他干擾信，并且特征函數C2計算取各參數平方的乘積結果作為識別判據，擴大了P波、P波尾波與S波的差異，使之更加易于識別。

4 結 論

提出一種針對實測爆破地震波S波震相初至時刻識別方法。通過對豐寧抽水蓄能電站地質勘探洞開挖現場實驗實測爆破振動數據分析，可以得到以下結論：

(1)該方法不對爆破振動信號進行任何濾波處理秉承“濾波越少，算法越好”的思想，直接利用P、S波在頻率、偏轉角、偏振度和橫向能量與總能量比值之間的差異對實測爆破地震波S波震相初至到時進行判斷，最終確定S波初至時刻。分析理論數值模擬結果表明，該方法在工程尺度下S波初至時刻識別與理論到時相差不大，誤差小于3%，驗證了該方法工程應用的合理性。

(2)應用的4個判別參數能夠正確的反映波形振動的特點，同時也能正確反映S波初至時刻前后振動波形的頻率，偏振特性及運動方向和能量的變化。

(3)對S波震相初至時刻識別是建立在P波到時已經被拾取的基礎上，其拾取精度受P波拾別精度的影響；同時在研究過程中忽略了R波的傳播及結構面等因素對S波初至時刻識別精度的影響。后期研究中需要更多的結合實例工程實測數據建立綜合考慮R波等影響因素的S波初至時刻識別模型。

[1] 金李,盧文波,陳明,等.節(jié)理巖體的爆破松動機理[J].爆炸與沖擊,2009,29(5):474-480.

JIN Li, LU Wenbo, CHEN Ming, et al. Mechanism of jointed rock loosing under blasting load[J]. Explosion and Shock Waves, 2009,29(5):474-480.

[2] 鐘冬望,吳亮,余剛.鄰近隧道掘進爆破對既有隧道的影響[J].爆炸與沖擊,2010,30(5):456-462.

ZHONG Dongwang, WU Liang, YU Gang. Effect of tunneling blasting on an existing adjacent tunnel[J]. Explosion and Shock Waves, 2010,30(5):456-462.

[3] 羅憶,盧文波,周創(chuàng)兵,等.高地應力條件下地下廠房開挖動態(tài)卸荷引起的變形突變機制研究[J].巖土力學,2011,32(5):1553-1560.

LUO Yi, LU Wenbo, ZHOU Chuangbing, et al. Mechanism study of abrupt deformation of underground powerhouse induced by excavation unloading under high in-situ stress[J]. Rock and Soil Mechanics, 2011,32(5):1553-1560.

[4] 朱俊,楊建華,盧文波,等.地應力影響下隧洞邊墻的爆破振動安全[J].爆炸與沖擊,2014,34(2):153-160.

ZHU Jun, YANG Jianhua, LU Wenbo, et al. Influences of blasting vibration on the sidewall of underground tunnel[J]. Explosion and Shock Waves, 2014,34(2):153-160.

[5] 宋彥輝,巨廣宏,孫苗.巖體波速與壩基巖體變形模量關系[J].巖土力學,2011,32(5):1507-1513.

SONG Yanhui, JU Guanghong, SUN Miao. Relationship between wave velocity and deformation modulus of rock masses[J]. Rock and Soil Mechanics, 2011,32(5),1507-1513.

[6] BAER M, KRADOLFER U. An automatic phase picker for local and tele-seismic events[J]. Bulletin of the Seismological Society of America, 1987,77(4):1437-1445.

[7] STEVENSON P R. Micro-earthquakes at flathead lake, montana: A study using automatic earthquake processing[J]. Bulletin of the Seismological Society of America, 1976,66(1):61-80.

[8] PARK J, VERNON F V, LINDBERG C R. Frequency dependent polarization analysis of high-frequency seismograms[J]. Journal of Geophysical Research, 1987,92(B12):12664-12674.

[9] ALLEN R. Automatic earthquake recognition and timing from single trace[J]. Bulletin of the Seismological Society of America, 1978,68(5):1521-1532.

[10] 田優(yōu)平.近震P波震相自動識別方法研究[D].北京:中國地震局地球物理研究所,2015.

TIAN Youping. Study on the automatic identification of local seismic P-phases[D]. Beijing: Institute of Geophysics China Earthquake Administration, 2015.

[11] DIEHL T, DEICHMANN N, KISSLING E, et al. Automatic S-wave picker for local earthquake to mography[J]. Bulletin of the Seismological Society of America, 2009,99(3):1906-1920.

[12] BAILLARD C, CRAWFORD W C, BALLU V, et al. An automatic kurtosis-based P-and S-phase picker designed for local seismic networks[J]. Bulletin of the Seismological Society of America, 2014,104(1):394-409.

[13] 周彥江,潘一山.基于小波變換的礦震波的P波和S波的識別[J].煤礦開采,2007,12(6):1-4.

ZHOU Yanjiang, PAN Yishan. P and S wave identification of rock-burst wave based on wavelet transform[J]. Coal Mining Technology, 2007,12(6):1-4.

[14] 劉希強,周惠蘭,沈萍,等.用于三分向記錄震相識別的小波變換方法[J].地震學報,2000,22(2):125-131.

LIU Xiqiang, ZHOU Huilan, SHEN Ping, et al. The method of wavelet transformation for 3-component record wave phase identification[J]. Acta Seismologica Sinica, 2000,22(2):125-131.

[15] 何先龍,佘天莉,高峰.一種地震P波和S波初至時間自動拾取的新方法[J].地球物理學報,2016,59(7):2519-2527.

HE Xianlong, SHE Tianli, GAO Feng. A new method for picking up arrival times of seismic P and S waves automatically[J]. Chinese Journal of Geophysics, 2016,59(7):2519-2527.

[16] 張楚旋,李夕兵,董隴軍,等.三函數四指標礦震信號S波到時拾取方法及應用[J].巖石力學與工程學報,2015,34(8):1650-1659.

ZHANG Chuxuan, LI Xibing, DONG Longjun, et al. A S-wave phase picking method with four indicators of three functions for micro-seismic signal in mines[J]. Chinese Journal of Rock Mechanics and Engineering, 2015,34(8):1650-1659.

[17] 曲保安,劉希強,蔡寅,等.近震S波震相實時自動識別方法研究[J].地震學報,2014,36(2):200-208.

QU Baoan, LIU Xiqiang, CAI Yin, et al. Method for real-time automatic identification of S-phase: Application to local seismicity[J]. ACTA Seismologica Sinica, 2014,36(2):200-208.

[18] CICHOWICZ A, GREEN R W E, BRINK A V Z. Coda polarization properties of high-frequency micro-seismic events[J]. Bulletin of the Seismological Society of America, 1988,78(3):1297-1318.

[19] BOORE D M. Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra[J]. Bulletin of the Seismological Society of America, 1983,73(6):1865-1894.

[20] KANASEWICH E R. Time sequence analysis in geophysics[M]. 3rd ed. Alberta: The university of Alberta press, 1981:80-112.

[21] 楊建華,盧文波,嚴鵬,等.全斷面開挖爆破產生的自由面對振動頻率的影響研究[J].振動與沖擊,2016,35(7):192-197.

YANG Jianhua, LU Wenbo, YAN Peng, et al. Influences of blast-created free surfaces on blasting vibration frequencies during full-face excavation[J]. Journal of Vibration and Shock, 2016,35(7):192-197.

[22] AMOROSO O, MAERCKLIN N, ZOLLO A. S-wave identification by polarization filtering and waveform coherence analyses[J]. Bulletin of the Seismological Society of America, 2012,102(2):854-861.DOI:10.1785/0120110140 .