利用小波包变换时频谱识别宁夏及邻区的地震和爆破＊

曾宪伟,赵卫明,李鸿庭,师海阔,姚 琳

(宁夏回族自治区地震局,银川 750001)

利用小波包变换时频谱识别宁夏及邻区的地震和爆破＊

曾宪伟,赵卫明,李鸿庭,师海阔,姚 琳

(宁夏回族自治区地震局,银川 750001)

采用 dmey小波基函数分别对地震和爆破事件的垂直向记录信号进行小波包变换,计算各事件信号的归一化时频谱值以及 P波和 S波时频谱值达到最大时的频率fmp和fms,比较地震信号和爆破信号 P波段(0～6.25 Hz)和 S波段(0～6.25 Hz)在各相同分解频带内的瞬时谱最大值差异,寻找合适的单项定量识别指标,并综合各单项识别指标形成综合识别判据。运用综合识别判据对银川台记录到的宁夏及邻区 14个地震事件和 19个爆破事件进行判别,结果表明,各单项定量识别指标的识别率均在 80%以上,综合判别结果均与事件的真实类型一致。

小波包变换;时频谱;地震;爆破;识别;宁夏及邻区

0 引言

地震记录易受外界因素干扰,这在一定程度上影响了人们对地震事件的识别和分析的有效性和准确性。在干扰因素中,像人工爆破这样的震动源干扰比较难以排除或识别。有些爆破的波形记录特征类似于天然地震,为测震分析人员的震相识别工作带来了一定的困难。为了解决爆破干扰造成的地震漏报或误报的问题,找到一种切实可行的识别天然地震和人工爆破的定量方法就显得尤为重要。

自 20世纪 70年代以来,国内外地震工作者在天然地震与非天然地震识别方面进行了广泛而深入的研究(杨选辉等,2005)。经过近 40年的探索,人们基于地震波的辐射图形、地震波谱分析和震相特征提出了很多识别地下非天然地震(特别是核爆破)和天然地震的判据,其中包括 P波初动振幅与 P波最大振幅比、勒夫波和瑞利波振幅比、P波与 S波谱振幅比(Plafcanet al,1997)、P波谱振幅与勒夫波谱振幅比(Walteret al,1995;Taylor,1996;Kimet al,1997)、瞬态谱等分析方法(刘希强等,2003)。另外,沈萍和郑治真(1999)采用非稳态的离散WD理论,对哈萨克斯坦的 5个地下核爆破事件及其邻近地区的 7个地震事件进行瞬态谱和瞬时频率特征分析,取得了区域地震和区域核爆破识别的半定量化结果。杨选辉等(2005)利用小波包对 7个地震事件和 5个核爆破事件进行多尺度分解,得到了信号在不同频带上的能量分布。刘希强等(2000)、林大超等(2003)还将小波包变换应用于地震信号或爆破信号的时频谱分析。

以上多数方法将不同震相之间的谱比值作为寻找识别指标的切入点,然后比较地震与爆破在该比值上的差异,而本文将小波包变换时频谱分析应用于天然地震与小当量人工爆破的识别,研究地震信号与爆破信号同一震相之间各分解频带内的时频谱值差异。

1 基于小波包变换的时频谱分析方法

连续小波变换频域能量守恒公式(飞思科技产品研发中心,2005)

Wf(a,b)为伸缩系数,b为平移分量,ψ(t)表示小波母函数,表示ψ(t)的傅立叶变换。

对(1)式进行离散化,参数a和b均取离散值,am=2-m,bm,n=nΔbm,m,n∈Z,Z表示整数集,则m和n分别代表尺度参数和时间参数,可得

将总能量依尺度m分段离散,再将频带上的能量和依时间离散,取m=-j,n=k,j∈Z+,k∈Z,Z+表示非负整数集,则在t=bm,n时刻,记录在尺度m频带的能量和为

2 资料选取

本研究范围在内蒙古阿拉善左旗与宁夏石嘴山 交 界 地 带(39°～ 39.25°N,105.75°～106.5°E),资料为银川地震台记录到的地震和爆破事件,事件记录选取信噪比较高的台站记录,记录采样率为 50 Hz。

笔者选取了研究区域内银川地震台记录到的波形较为清晰的 14个地震事件,地震目录取自宁夏地震观测报告(表 1)①宁夏回族自治区地震局 .2007.宁夏地震观测报告(2003～2007).,其中事件 1、6、7和10均具有 P、S波波列清晰、震相较为单一、高频成分较为丰富的波形特征,说明事件 1、6、7和10的记录信号与典型的地震信号较为接近,可以较为可靠地判定为地震事件。

表 1 研究范围内的地震事件目录Tab.1 Earthquake catalogue in the study area

同样,笔者选取了研究区域内银川地震台记录到的波形较为清晰的 19个爆破事件(表 2),其中序号为 1～7的事件为经过落实后确认的爆破事件,序号为 8～19的事件为分析人员通过对波形特征作直观分析后识别出的爆破事件。

表 2 研究范围内的爆破事件目录Tab.2 Explosion catalogue in the study area

另外,对所选的全部爆破和地震事件的记录波形进行分析,发现它们的面波均不发育,因为这些事件的震级都较小(均小于ML3.0),震中距也较小(均小于 80 km)。

3 资料处理

在某个给定的分解尺度上,小波包变换的每组变换系数分别对应于一个频率区间上信号的变化特征。原始信号在给定频率范围内关于时间的变化情况直接通过其对应的变换系数得到刻画,因此这些变换系数就能给出信号的时间频率特征。因此,我们也可以把用变换系数表示的信号特征作重构。显然,这与直接应用变换系数进行分析没有本质区别,只是从表面上更易于被人们接受而已。不过,已建立的小波基函数只能在一定程度上满足这个要求,因此二者的局部差异是不可避免的,但不影响关于信号特征的分析。

事实上,一个时间信号可以表示为不同分解层的重构信号的叠加(多分辨分析),这些叠加信号中每个频率区间的信号与其它频率区间上的信号彼此正交。于是,信号的总能量可以通过相关频率区间上信号的能量进行完全描述。而能量(密度)是相关点分解后导出的速度值的平方,因此,通过计算速度值的平方就可以给出不同事件的时频谱值。

笔者仅选一对实例进行具体分析,分别选取表 1和表 2中的第一个事件,这两次事件的银川台垂直向原始记录波形如图 1a和图 2a所示。首先对原始记录信号进行归零、去倾以及小波包降噪处理(图 1b,图 2b),然后对降噪后的地震信号和爆破信号进行尺度j=5的小波包分解,此时整个信号频带 0～25 Hz平均分成了 32个频带。计算每个频带内的分解信号速度的平方值,即时频谱值,并计算每个事件信号的时频谱值与各自信号的总能量之比,即进行谱值归一化处理,以便于对不同事件信号间进行谱值大小比较。最后通过比较地震信号与爆破信号的 P波段时频谱值达到最大时的频率差异,以及 S波段时频谱值达到最大时的频率差异,寻找定量识别的阈值大小。通过比较地震与爆破 P波信号(0～6.25 Hz)在各分解频带内的时频谱最大值差异以及地震与爆破 S波信号(0～6.25 Hz)在各分解频带内的时频谱最大值差异,寻找不同震相信号在各分解频带内的识别阈值大小。

图 1 地震原始波形(a)和滤波后的波形(b)Fig.1 Original signal(a)and filtered signal(b)of the earthquake

图 2 爆破原始波形(a)和滤波后的波形(b)Fig.2 Original signal(a)and filtered signal(b)of the explosion

4 资料分析

采用 dmey小波基函数对所选地震和爆破事件的垂直向记录进行小波包变换,计算各事件信号的归一化时频谱值、P波和 S波时频谱值达到最大时的频率fmp和fms、地震信号(0～6.25 Hz)P波和 S波在各分解频带内的时频谱最大值和爆破信号(0～6.25 Hz)P波和 S波在各分解频带内的时频谱最大值。

以事件序号为横坐标,分别以频率fmp和fms为纵坐标绘图(图 3)。由图 3可以看出,fmp和fms分别以 4 Hz和 3 Hz为界,识别效果较为明显,可作为两项识别指标。

图 3 地震与爆破银川台垂向记录信号 P波(a)和S波(b)时频谱值达到最大时的频率差异Fig.3 Differences of frequency when the value of ti me-frequency spectrum of P(a)and S(b)wave of the earthquake and explosion from Yinchuan station's UD records is maximum

分析 0～6.25 Hz频带内的地震信号和爆破信号,以每个小波包分解频带为单元,事件序号为横坐标,时频谱最大值为纵坐标绘图(图 4,图5)。地震与爆破的 P波在各分解频带内的时频谱最大值差异见图 4,S波在各分解频带内的时频谱最大值差异见图 5。可以看出,地震与爆破 P波时频谱最大值仅在 0.78125～1.5625 Hz和 1.5625～2.34375 Hz两个频带内存在较显著差异,分界阈值分别为 0.2和 0.35;而 S波时频谱最大值在 0～0.78125 Hz、0.78125～1.5625 Hz、1.5625～2.34375 Hz、3.125 ～3.90625 Hz、3.90625 ～4.6875 Hz和 4.6875～5.46875 Hz 6个频带内均存在显著差异,分界阈值分别为 0.5、1、3.2、4、0.5和 0.7。因此,P波两个频带内的时频谱最大值差异和 S波 6个频带内的时频谱最大值差异均可作为地震和爆破的识别指标。

图 4 银川台垂向记录的地震与爆破信号 P波(0～6.25 Hz)各分解频带内的时频谱最大值差异Fig.4 Differences of maximum value of time-frequency spectrum in decomposing bands of P signal(0～6.25 Hz)of the earthquake and explosion from Yinchuan station's UD records

图 5 银川台垂向记录的地震与爆破信号 S波(0～6.25 Hz)各分解频带内的时频谱最大值差异Fig.5 Differences of maximum value of time-frequency spectrum in decomposing bands of S signal(0～6.25 Hz)of the earthquake and explosion from Yinchuan station's UD records

表 3给出了各单项定量识别指标的识别阈值以及识别率,各单项识别指标的识别率均在 80%以上。为提高识别指标的有效性,笔者将 10个单项识别指标结合在一起,原则为超过半数的识别指标给出的识别结果为事件类型的判别结果。依据该原则,对本文中的 14个地震事件和 19个爆破事件重新进行判别,其判别结果均与事件的原类型一致(表 4和表 5),尤其是 6个落实为爆破的疑似爆破事件均判别为爆破。事件 3的速报定性结果为地震,经落实是爆破,利用该综合识别判据可以较好地将事件 3判断为爆破事件。另外,利用该综合识别判据,较为典型的地震事件 1、6、7和10也判定为地震。

表 3 银川台地震事件和爆破事件的单项识别指标及阈值Tab.3 Identification index and thresholds of Yinchuan station

表 4 石嘴山大峰煤矿爆破的银川地震台垂直向记录的识别结果Tab.4 Identification result of explosion of Dafengmine in Shizuishan

5 判据例证

2007年 12月 20日 11时 30分,宁夏石嘴山大峰煤矿羊齿采区进行硐室爆破,爆破炸药量为5499 t,爆破点距离银川地震台 48 km,宁夏测震台网测定的震级为MS3.8。笔者采用上文所述处理方法,对该爆破的银川地震台垂直向记录进行小波包变换分解,并计算相应的识别指标参数取值(表 4),得到 10项单项识别指标中有 8项的判别结果为爆破,因此认定该事件为爆破事件,判定类型与实际类型一致,说明本文给出的识别判据是有效的。

6 结论与讨论

地震的发生是岩石破裂过程的结果(陈运泰,顾浩鼎,2004),其能量释放相对爆破要缓慢得多,所以相对而言地震信号能量多集中于低频段,而爆破信号能量多集中于高频段。然而,爆破信号的一个特点是距离爆源较近时地震波的高频成分较为丰富,且持续时间较短(阳生权,2002)。随着爆破地震波在岩土介质体内从爆源向四周的传播,介质体的内阻尼使地震波幅值衰减,且对于高频振动阻尼作用较大。由于爆破信号在传播过程中有较长的路程在较为松散的浅层,远距离范围的质点振动高频成分衰减得非常快,低频成分则相对增大。

由于本文所选资料仅为银川地震台所记录的某一区域范围内的地震和爆破,得到的结果也仅是区域地震和区域爆破的识别判据,如果对不同台站同一区域以及同一台站不同区域的事件进行研究,有可能得到更多有意义的结果。另外,随着分解尺度增加,信号在频率域分解得更细,同一事件不同频带内的时频谱差异性也会更大,信号的细节特征会得到更详细的描述,这对于地震和爆破识别或许会得到更多有价值的结果。

曹晖,赖明,白绍良 .2004.地震地面运动局部谱密度的小波变换估计[J].工程力学,21(5):109-115.

陈顺云,杨润海,王赟赟,等 .2006.小波分析在声发射资料处理中的初步应用[J].地震研究,25(4):328-334.

陈运泰,顾浩鼎 .2004.震源理论基础(上册)[M].北京:地震出版社 .

飞思科技产品研发中心 .2005.小波分析理论与 MATLAB7实现[M].北京:电子工业出版社 .

李弼程,罗建书 .2003.小波分析及其应用[M].北京:电子工业出版社.

林大超,施惠基,白春华,等 .2003.爆炸地震效应的时频分析[J].爆炸与冲击,23(1):31-36.

刘希强,沈萍,张玲,等 .2003.用小波变换能量线性度方法识别天然地震与爆破或塌方[J].西北地震学报,25(3):204-209.

刘希强,周惠兰,李红 .2000.基于小波包变换的地震数据时频分析方法[J].西北地震学报,22(2):143-176.

沈萍,郑治真 .1999.瞬态谱在地震与核爆识别中的应用[J].地球物理学报,42(2):233-240.

阳生权 .2002.爆破地震累积效应理论和应用初步研究[D].长沙:中南大学 .

杨选辉,沈萍,刘希强,等 .2005.地震与核爆破识别的小波包分量比方法[J].地球物理学报,48(1):148-156.

曾宪伟,赵卫明,盛菊琴,等 .2008.应用小波包识别宁夏及邻区的地震和爆破[J].地震研究,31(2):142-148.

Plafcan D,Sandvol E,Seber D,et al.1997.Regional discrimination of chemical explosions and earthquakes:A case study in Morocco[J].BSSA,87(5):1126-1139.

Taylor S R.1996.Analysis of high-frequency Pg/Lg ratios from NTS explosions and western U.S.earthquakes[J].BSSA,86(4):1042-1053.

WalterW R,Mayeda KM,Patton H J.1995.Phase and spectral ratio discrimination between NTS earthquakes and explosions.Part I:Empirical observations[J].BSSA,85(4):1050-1067.

Kim W Y,Aharonian V,Lerner-Lam A L,et al.1997.Discrimination of earthquakes and explosions in southern Russia using regional highfrequency three-component data from the IRIS/JSP Caucasus network.[J].BSSA,87(3):569-588.

Discrimination between Earthquakes and Explosions in Ningxia and Its Neighboring Region Using Time-frequency Spectrum Wavelet Packet Transform

ZENG Xian-Wei,ZHAO Wei-ming,L I Hong-ting,SHI Hai-kuo,YAO L in
(Earthquake Administration of Ningxia Hui Autonomous Region,Yinchuan750001,Ningxia,China)

Usingdmey wavelet-base function,We conduct a wavelet-packet transform of the signals of vertical component recordings of the earthquakes and the explosions to calculate the values of the normalized time-frequency spectrum of the signal of each earthquake and explosion,and to calculate the frequency fmpon condition that the value of time-frequency spectrum of P wave is maxim um.We also calculate the frequency fmson condition that the value of time-frequency spectrum of S wave is maximum.Then We compare the maximum value of instantaneous spectrum of the seismic signals with that of the blasting signals in the same decomposition frequency band of P wave within frequency0-6.25Hz.And We compare the maximum value of instantaneous spectrum of the seismic signals with that of the blasting signals in the same decomposition frequency band of S wave within frequency0-6.25Hz.Thus We get some quantitative indexes to identify the earthquake and the explosion.Combining all these indexes We get a synthetic criterion.Then We apply the single index and the synthetic criterion to the identification test of14quakes and19explosions recorded by Yinchuan seismic Station.We find that using the single index,80%of the quakes and the explosions can be identified,and100%of the quakes and the explosions can be identified using the synthetic criterion.

wavelet packet transform; ti m e-frequency spectrum; earthquakes; explosions; identification;Ningxia and the neighboring region

P315.63

A

1000-0666(2010)03-0300-08

2009-09-29.

中国地震台网中心震情跟踪项目(2009A51)资助 .