爆破地震波中S波识别方法及其应用*

杨招伟，卢文波，高启栋，陈 明，胡浩然，严 鹏，王高辉

(1.武汉大学水资源与水电工程科学国家重点实验室,湖北 武汉 430072;2.武汉大学水工岩石力学教育部重点实验室,湖北 武汉 430072)

近年来，我国正在兴建的一批大型水电工程都处在西南、西北地区的崇山峻岭之中，均需进行高陡边坡和地下洞室的爆破开挖，由此导致一系列的岩石工程安全问题，特别是在爆炸荷载、地震荷载等动荷载作用下的变形和稳定问题[1-4]。因此，对岩石各种动力学参数的研究就显得更加重要，如何准确地确定现场岩石动力学参数是工程施工过程中需着重考虑的一个方面。现阶段工程中常用的岩石参数获取方法有室内实验法和原位实验法，但由于室内实验针对小尺度岩样进行研究，因而得到的也是小尺度岩样的动力学参数，由于岩体具有较为明显的尺度效应，为了将室内实验测试的小尺度岩样动力学参数的研究过渡到现场大尺度岩体，因此，现场原位实验具有不可替代的作用，然而原位实验也有周期长，费用高等缺点[5]。而通过测定爆破地震波(P、S波)在岩石中的传播速度反演岩石动力学参数的方法较为简单快捷。如何精确的判别P、S波震相初至时刻对岩石动力学参数反演工作的效率和精度有着不可忽视的影响，因此根据实测振动波形进行S波的初至时刻识别值得开展深入研究。

迄今为止，水利工程中的爆破地震波的识别研究的方法和技术直接来自于对天然地震活动的研究。对于天然地震，针对P、S波初至时刻判别方面进行了大量的研究。Baer等[6]指出P波初至时刻的自动拾取算法主要是基于监测振动信号和噪声的不同特征作为P波震相识别的判据。Stevenson等[7]首次提出利用长短时窗平均比来自动拾取P波震相初至时刻。如今发展出的P波到时自动识别方法主要有：自回归模型、AIC法、波形变化增长法、偏振分析法及长短时间平均方法等[8-10]。但在S波自动识别方面，由于P波尾波以及各种反射波相的相互干扰，以上的方法用于S波震相自动识别的效果较差。Diehl等[11]提出了用于地震层析成像的S波自动识别方法，该方法需要对振动信号做STA/LTA判定、偏振分析等处理，计算量较大，比较适合做后续分析；Baillard等[12]利用峰态特征函数和特征向量的分解拾取三分像地震信号的P、S波震相。刘希强等[13-14]提出利用振动信号在不同尺度的小波变换中偏振信息的不同来确定S波震相的小波变换方法；何先龙等[15]提出一种二次自回归模型拾取P、S波的初至时间；张楚旋等[16]采用短时平均过零率、短时能量等4个指标作为S波自动拾取的指标。

虽然学者们在天然地震领域上S波的拾取上取得了丰富的研究成果，但是与天然地震信号相比，水利工程爆破开挖产生的微震信号持续时间较短、背景噪音复杂、分布范围较小，受结构面等地质的不均匀性影响很大，导致P、S波分离不明显，致使适用于地震领域的S波识别算法在水利工程尺度下识别效果不明显。因此，对水利工程爆破开挖产生的微震信号的S波自动识别进行更深入的研究具有较为重大的意义。

本文中旨在研究根据实测爆破开挖产生的振动波形进行S波震相的快速识别，在满足安全生产精度要求的前提下，秉承信号处理上“滤波越少，算法越好”的原则[17]，对原始振动数据不进行任何的滤波处理。基于其他方法已经得到的P波震相到时，在此基础上，对前人识别方法进行改进，根据P波、S波在频率、偏振方向和水平能量与总能量上的差异对S波进行精确识别。对于爆破地震波中P波震相的初至时刻识别，推荐采用Baer等[6]提出的自动识别算法。

1 识别方法及特征函数

水利水电工程开挖爆破过程中，炸药能量的很大一部分以地震波的形式向四周传播，引起地面振动。爆破地震波在传播过程中可被分为体波和面波，其中根据传播过程中质点运动和传播方向关系，可以将体波进一步分为P波和S波，相较于P波，S波具有传播速度仅次于P波，并且其周期较长、振幅较大的特点。Diehl等[11]指出，振动信号监测仪器会最先监测到P波到达引起的振动信号，而相对于较为容易的识别出P波初至震相时刻，在S波震相识别方面，由于P波尾波信号以及各种转换波信号，如PS波的强干扰导致S波淹没在P波中难以识别。可根据P、S波周期及偏振特性的差异，对与周期差异相关的短时平均过零率及偏振特性差异相关的偏转角、偏振度及横向能量与总能量差异等参数进行分析，现分别讨论。

1.1 确定计算窗长

在对上述识别参数进行计算时，当计算窗长选取过长，上述特征参数的值变化不明显；计算窗长选取过短，则导致识别参数值的不稳定。Cichowicz等[18]提出根据P波的卓越频率确定时窗长度的方法。计算时窗长度与P波卓越频率成反比关系，具体表达式如下：

(1)

式中：l为计算窗长，Δt为仪器采样间隔时间，fp为P波的卓越频率。

Boore等[19]提出P波的卓越频率与振动波形的速度功率谱密度和相应的位移功率谱密度存在如下式所列关系：

(2)

1.2 短时平均过零率

由于P波的周期小于S波周期，在监测信号上就表现出如图1所示的波形前后疏密不同的变化，为了更好地识别S波震相初至时刻，需对这种疏密变化进行量化。对于平稳信号来说，可以直接用周期或频率直接衡量。而对于地震信号、爆破产生的振动信号、噪声信号等非平稳信号而言，由于其周期的不稳定，则用短时平均过零率来进行度量。

振动信号单位时间内过零的次数可以定义为短时平均过零率K1，其值可以作为S波初至震相时刻拾取的指标。其定义式为：

(3)

由式(3)可知，短时平均过零率K1的值与选取的时窗长短有关。关于S波的初至时刻震相，因此主要研究S波到达前后振动信号周期变化的差异。将P波初至时刻到S波初至时刻这一区间称为S-P区间，将S波震相初至时刻后(ns-np)Δt时窗内的区间，称为S′-S区间。定义公式为：

(4)

(5)

式中：K2为S-P区间内的短时平均过零率，K3为S′-S区间内的短时平均过零率，ns′=ns+(ns-np)，np、ns分别为P波初至时间和S波初至时间对应的信号点序号。

根据张楚旋等[16]对建立的微震数据库的信号研究可知，短时平均过零率与振动信号的S波震相的初至时刻存在着对应关系，定义S波初至时刻震相前后的短时平均过零率之比C1=K2/K3，计算发现，C1的取值范围大于1.2。

1.3 偏转角

确定P波初至到时后，结合Kanasewich[20]提出的偏振特性概念，选取P波到时后l长度的计算其协方差矩阵最大特征值所对应的特征向量，即P波的偏振方向。计算式如下：

(6)

P波初至时刻点为计算起点，取时窗长度为l个采样点往后开始进行逐点进行滑动，求得各个计算窗口最大特征值所对应的特征向量，并计算其与P波偏振方向的夹角即偏转角α。对偏转角α进行归一化处理得：

(7)

1.4 偏振度

P、S波具有较高的线性偏振度，而P波尾波则更多的表现为椭圆偏振的特性。Cichowicz等[18]提出直接根据实测地震振动数据求取偏振度的方法。计算公式为：

(8)

式中：λ1、λ2、λ3为各个数据窗对应的协方差矩阵的特征值。

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的理论统计平均值如表1所示。

对于实测振动数据而言，由于受P波尾波及P-S

表1 统计平均值Table 1 Statistical averages

转换波的影响，特征参数F2、F3、F4的统计平均值难以精确计算。因此对于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区间(1～l)、S-S′区间(l～l+(l-1))内的短时平均过零率，然后再按照C1=K2/K3计算特征函数C1值。对S波初至时刻进行初步判定，缩小加权识别特征函数Cw的计算空间，提升此方法计算效率。根据式(7)、(8)、(10)分别计算出识别参数F2、F3、F4的值并按式(12)计算加权特征函数Cw的值，最终确定S波震相初至时间。

2 数值模拟及对比分析

为了进一步验证上述S波识别方法的正确性,基于LS-DYNA动力有限元软件，模拟半无限岩体中柱状药包爆破诱发的振动传播规律。计算中岩石物理力学参数及拉压损伤模型参数取值见表2[21]。因模型具有对称性，采取1/4模型模拟岩体中单孔爆破特性。模型尺寸为10 m×9 m×20 m，单元240 304个，节点260 012个，由于炮孔直径较小，网格划分采用从炮孔到边界渐变方式剖分，炮孔附近单元最小尺寸为0.003 6 m，边界最大尺寸0.8 m，岩体采用solid 164单元模拟。炮孔直径90 mm，装药直径70 mm，孔深3 m，堵塞0.9 m，在模型顶面布置振动监测点，有限元计算模型如图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 工程概况

丰宁抽水蓄能电站位于中国河北省丰宁满族自治县境内，电站工程分两期开发建设,一期工程和二期工程装机容量分别为1 800 MW，装机6台，单机容量300 MW，在电网系统中承担调峰、调频、调相和事故备用任务。利用二期工程地下厂房地质勘探洞开挖的时机和条件，进行竖直钻孔爆破实验。实测资料表明该地质勘探洞以花岗岩为主，总体质量较好，因此可以认为在本次实验范围内岩体为均质、各向同性的弹性体，实验过程中炮孔及振动监测点布置如图4所示，采用半秒微差延期爆破，逐孔起爆。爆破设计参数见表4，表4中D为孔径，h为孔深，d为药卷直径，L装药长度，Lb堵塞长度，m单响药量。利用振动监测仪器得到了各个测点的振动时程曲线，利用上述方法分析实测振动波形图对S波的初至时刻进行精确识别。

基于《爆破安全规程》、《水电水利工程爆破安全监测规程》等规范标准开展爆破振动测试，爆破振动测试系统由三向速度检波器、信号采集与记录设备、数据处理系统3个部分组成。主要测试设备为TC-4850爆破振动智能监测仪。在实验过程中，TC-4850爆破振动智能监测仪用于监测爆破振动速度，其采样频率为8 000 Hz。结合爆破实验进行爆破振动监测，在探洞中布置爆破振动测试线，如图4，每个测点相对于爆源具有水平径向、水平切向及垂直向3个方向，以获得水平钻孔爆破不同工况的爆破振动的传播规律。实验中，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节∅32 mm炸药捆绑而成。

3.2 S波初至时刻识别效果

应用上述S波震相加权识别特征函数结合前文所述的识别流程对丰宁二期工程地下厂房探洞开挖实测振动数据进行S波初至时刻识别，识别结果见表5。限于篇幅，只给出典型测点(10#测点)实测振动波形S波初至时刻的识别效果，如图5所示。

表5 实验S波震相初至时刻识别结果Table 5 S-wave arrival identification in experiment

从以上数值模拟分析及丰宁水电站实测波形的应用效果可知，本文中所述方法对S波震相初至时刻的识别率较高，且具有简单、快捷等优点，可以满足实际水利工程的要求。由于该方法是建立在P波已经被识别的基础上的，下面对其识别参数做进一步论述。

(1)从操作简单快速来看，根据周期差异作为识别函数C1对S波进行初步识别速率最快，但由于P波尾波及转换波对频率的影响导致该特征函数的识别准确率较低，可用于缩小加权识别特征函数Cw计算空间，提升计算效率。

(2)根据横向能量与总能量比值识别参数对S波震相的识别效果最优。造成该现象的主要原因是因为经过坐标系变换，把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] 金李,卢文波,陈明,等.节理岩体的爆破松动机理[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] 罗忆,卢文波,周创兵,等.高地应力条件下地下厂房开挖动态卸荷引起的变形突变机制研究[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] 田优平.近震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 .