利用小波包變換時頻譜識別寧夏及鄰區的地震和爆破＊

曾憲偉,趙衛明,李鴻庭,師海闊,姚 琳

(寧夏回族自治區地震局,銀川 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)資助 .