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

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

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

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

∫_R│F(ω)│²dω=1/(C_ψ)_R²[W_f(a,b)]²(da)/(a²)db.(1)

其中,C_ψ=∫_R(│ψ^{^}(ω)│²)/(│ω│)dω,

W_f(a,b)=1/((│a│)^1/2)∫^∞∫_-∞f(t)ψ((t-b)/a)^-dt, a为伸缩系数, b为平移分量, ψ(t)表示小波母函数, ψ^{^}(ω)表示ψ(t)的傅立叶变换。

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

∫_R│F(ω)│²dω=∑_m2^m∑ⁿΔb_m[W_f(a_m,b_m,n)]².(2)

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

∫_ωj│F^j_k(ω)│²dω=[W_f(a_j,b_j,k]².(3)

取规范正交小波包基{2^j/2W_n(2^jt-k),n∈Z₊,k∈Z},则记录在t=t_k时刻n频带内的能量和为

∫_ωn│Fⁿ_k(ω)│²dω=[C^n,j_k]².(4)

其中C^n,j_k=∫_Rf(t)2^j/2W_n(2^jt-k)^-dt, j=1,2,3,…; n=0,1,2,…,2^j-1; k=1,2,…。该式展现了样本函数能量在频率和时间上的分布情况,因此得到了信号的时频局部特性——局部谱密度(曹晖等,2004)。

本研究范围在内蒙古阿拉善左旗与宁夏石嘴山交界地带(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的记录信号与典型的地震信号较为接近,可以较为可靠地判定为地震事件。

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

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

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

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

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

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

笔者仅选一对实例进行具体分析,分别选取表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

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

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

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

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

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

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

图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

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

地震的发生是岩石破裂过程的结果(陈运泰,顾浩鼎,2004),其能量释放相对爆破要缓慢得多,所以相对而言地震信号能量多集中于低频段,

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

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

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