郝春月，鄭 重

（中國地震局地球物理研究所，北京100081）

1 引 言

2006年10月9日1時35分，在北緯41.29°、東經129.09°、深度0km 處發(fā)生一起爆炸事件［1］，此次爆炸事件距離全球地震臺網（GSN）的牡丹江（MDJ）地震臺約3.35deg（372.49km），方位角為186.4°；2009年5月25日0時54分（UTC），在北緯41.3°、東經129.04°、深度0km 處又發(fā)生一起爆炸事件［2］，該事件的位置距離MDJ約3.34deg（371.37km），方位角為187.1°。這2次爆炸事件距離很近。為研究2次爆炸事件的能量比與相似性，本文中試圖對MDJ地震臺3個分向（Z、N、E）記錄的2次事件進行波形記錄單振幅比、功率譜比、相關值和相干函數(shù)的計算與分析。

2 資 料

2次事件均被MDJ地震臺記錄。MDJ地震臺位于中國黑龍江省東南部，屬于全球地震臺網（GSN），本文中選用的數(shù)據均來自MDJ臺。圖1顯示了MDJ臺3個分向（即垂直向、北南向與東西向）在2006年10月9日與2009年5月25日記錄的爆炸事件波形。

圖1顯示了25s的波形數(shù)據，橫軸表示采樣點，縱軸為振幅。每幅圖上顯示的2個重要震相分別為Pn（首波，即Moho面繞射縱波）與Pg（直達縱波），在每組震相的最大值處標注了最大幅值。2次爆炸記錄的Pn與Pg震相在3個分向（Z、N、E）之間的最大振幅的單振幅比值（A09／A06）分別為4.794、2.716、4.402和3.929、3.192、3.950（見圖2）。

3 功率譜比的計算

為了解2次爆炸事件在頻率域的能量分布情況，對2次爆炸事件進行了功率譜的計算與分析，功率譜分析采用Welch 平均周期圖法。Welch 平均周期圖法是對直接法的改進，即把一長度為N 的數(shù)據xN（n）分成L 段（在分段時可允許每一段的數(shù)據有部分的重疊），每一段的長度為M。分別求每一段的功率譜，然后加以平均。每一段的功率譜為［3］

圖1 MDJ地震臺記錄的2次爆炸事件的波形圖（25s時間段）Fig.1 Waveforms of explosion events recorded by MDJ seismic station

圖2 2次爆炸事件記錄的Pn 與Pg 震相在3個分向之間的最大振幅的單振幅比Fig.2 Ratios of the maximum amplitude of Pnand Pgin three components recorded by the two explosion events

在此分別對2個時間窗進行功率譜估計，一個為P波組的30s時間窗，另一個為S波和面波組的200s時間窗。首先分別對2組時間窗進行1～4Hz和0.1～2Hz的濾波（根據圖3（圖中灰度表示功率譜密度），即P震相在1～4Hz頻段能量最大，S與面波震相在0.1～2Hz頻段能量最大）。進行功率譜計算所用的30sP波組數(shù)據為1 200點（采樣率為0.025Hz），利用的漢寧窗窗長為256點，重疊128點。最后對功率譜進行了5點平滑。另外S波和面波組的時間窗為200s，數(shù)據為8 000點，利用的漢寧窗窗長為1 024點，重疊512點。最后對功率譜進行了5點平滑。2次事件的功率譜計算完畢后，可以得出2次事件的功率譜點對點的比值。

由于計算功率譜是為了對2次事件進行對比，并且MDJ臺的傳遞函數(shù)在2次事件之間沒有變化，所以沒有對波形進行去除儀器傳遞函數(shù)的計算。因為在計算2次事件的譜值比過程中，會抵銷傳遞函數(shù)，所以忽略了這個步驟。忽略這個步驟的結果是，圖4～5顯示的2次事件功率譜的縱軸單位不是（m／s）2，也不是表示能量衰減所用的dB，而是量綱為一。

圖3 2006年與2009年的爆炸事件在MDJ臺垂直向記錄波形的時頻譜圖Fig.3 Time－frequency spectrum analysis of the two explosion events recorded in the vertical component

圖4 MDJ地震臺記錄的2次爆炸之間P波組功率譜與譜值比Fig.4 Power spectrum densities of the P phase group of the two explosion events and the ratios of them

圖5 MDJ地震臺記錄的2次爆炸事件之間的S波和面波組功率譜與譜值比Fig.5 Power spectrum densities of the S－Surface wave group of the two explosioneventsandtheratiosofthem

從圖4中可以看出，MDJ臺記錄的P 波組震相在垂直向的譜比值比較穩(wěn)定，平均值在4.6左右。水平向的譜比值在2Hz左右出現(xiàn)最大值，可達到8.0左右。但在1～4Hz的頻段內，北南向與東西向的譜比平均值分別為3.8和4.8。圖5中，MDJ臺記錄的S波和面波組震相在3個方向的譜比值在4.6左右變化。功率譜代表了能量，由2次爆炸事件的P波組和S波組的功率譜比值，可估計2009年爆炸事件所釋放的能量是2006年爆炸事件釋放能量的4～5倍。

4 2次爆炸事件的相關性分析

為探索2次事件是否為相似事件，對MDJ臺記錄的2次事件的3個方向的數(shù)據進行了相關性分析。首先給出2次事件的波形疊加圖（圖6），從圖中可以看出，2次事件確實具有相似性。

為了對相似性有一個量的描述，對2次事件的互相關值進行計算。首先對波形進行1～4Hz和2～4Hz的帶通濾波（根據圖3，P 波組震相在前幾秒的能量主要集中在1～4 Hz）。然后選用2次爆炸前2s的數(shù)據進行互相關值的計算。選擇相關值計算的時間窗一般為1～2s［4－5］，在此選擇2s，使波形窗能包含更多的震相信息。

2個時間序列的互相關函數(shù)可以由下式表示［6］

圖6 2次事件的3個方向的波形記錄疊加Fig.6 Overlap of the waveforms of the two events recorded

式中：x（i）和y（i）分別是2個測點x和y 記錄的噪聲時間序列，N 表示數(shù)據樣本點的個數(shù)和表示N 個樣本點的平均值。

根據公式（3），計算得出MDJ臺3個方向（Z、N、E）記錄的2次事件在1～4 Hz頻段的互相關值分別為0.889 2、0.814 4、0.930 2（圖7（a）），而在2～4 Hz頻段的互相關值分別為0.931 3、0.884 1、0.935 1（圖7（b））。經過2～4Hz帶通濾波的波形相關值比經過1～4Hz帶通濾波的波形互相關值要高一些，表明了前2s數(shù)據在2～4Hz的能量要高些，并且互相關值平均為0.92。而在1～4Hz的互相關值平均為0.88。

為了能夠在頻率域顯示2次事件波形的相似性，對2次事件進行相干函數(shù)的計算。在頻域中，經常用到可用功率譜表示的相干函數(shù)rxy（f），相干函數(shù)rxy（f）（又稱凝聚函數(shù)）給出了2個隨機信號在頻率域的相似性，其中f 表示2個隨機信號的頻率。一般來說，相干函數(shù)的值在0～1之間變化。若2個信號完全不相關，則相干函數(shù)值為零。對于2個相同的信號，其相干函數(shù)為1。相干函數(shù)定義為［7］

式中：Gxx（f）和Gyy（f）分別為2個隨機信號x（t）和y（t）的自功率譜密度；Gxy（f）為它們的互譜密度。

2次事件的波形經過1～4Hz帶通濾波后，利用10s的P 組震相數(shù)據進行相干函數(shù)計算（圖8）。選擇10s，主要考慮盡可能包含足夠的P 組震相信息，還要考慮不要包含其他震相信息（如S震相）。10s的時間窗主要包含2個P組震相，即Pn與Pg震相，Pg震相在距Pn震相約7.3s處，而Sg（直達橫波）震相在距Pg震相約40s處（圖3），沒有包含在內。所以10s時間窗主要包含P組震相信息。由圖8可以看出，在1～4Hz頻段內，MDJ地震臺3個方向記錄的2次事件的相干函數(shù)均保持在0.8左右。

圖7 2次事件的相關值計算結果Fig.7 The correlation results of the two events

圖8 用來計算2次爆炸事件相干函數(shù)的數(shù)據波形（10s）與2次事件的相干函數(shù)結果（1～4Hz）Fig.8 Waveforms（10s）used to calculate the coherence functions and the coherence function results of the two explosion events（1～4Hz）

5 結 論

（1）2次爆炸事件的能量分布相似，主要分布在0.1～4Hz之間。其中P組震相的能量主要分布在1～4Hz，S波和面波組震相的能量主要分布在0.1～2Hz。（2）2次爆炸事件記錄的Pn與Pg震相在3個分向最大振幅的單振幅比（A09／A06）平均為3.97和3.69（圖2）。MDJ臺3個方向記錄的2次爆炸事件的P波組、S波和面波組功率譜比（P09／P06）平均為4.5。由此估計，2009年爆炸事件所釋放的能量是2006年爆炸事件釋放能量的4～5倍。（3）在1～4Hz和2～4Hz頻段內，2次事件3個方向的平均互相關值分別為0.88和0.92。在1～4Hz頻段內，P波組的相干函數(shù)平均保持在0.8左右。

衷心感謝許忠淮研究員的指導，感謝張爽同志提供的數(shù)據。

［1］ USGS National Earthquake Information Center，US.Geological survey earthquake database［DB／OL］.America：USGS，1940［2010－10－8］.http／／：neic.usgs.gov／cgi－bin／epic／epic.cgi？SEARCHMETHOD＝1＆FILEFORMAT＝4＆SEARCHRANGE＝HH＆SYEAR＝2006＆SMONTH ＝10＆SDAY＝9＆EYEAR＝2006＆EMONTH ＝10＆EDAY＝9＆LMAG＝3＆UMAG＝6＆NDEP1＝0＆NDEP2＝2＆IO1＝＆IO2＝＆SLAT2＝0.0＆SLAT1＝0.0＆SLON2＝0.0＆SLON1＝0.0＆CLAT＝0.0＆CLON＝0.0＆CRAD＝0＆SUBMIT＝Submit＋Search.

［2］ USGS National Earthquake Information Center，US.Geological survey earthquake database［DB／OL］.America：USGS，1940［2010－10－8］.http／／：neic.usgs.gov／cgi－bin／epic／epic.cgi？SEARCHMETHOD＝1＆FILEFORMAT＝4＆SEARCHRANGE＝HH＆SYEAR＝2009＆SMONTH ＝5＆SDAY＝25＆EYEAR＝2009＆EMONTH ＝5＆EDAY＝25＆LMAG＝4＆UMAG＝5＆NDEP1＝0＆NDEP2＝2＆IO1＝＆IO2＝＆SLAT2＝0.0＆SLAT1＝0.0＆SLON2＝0.0＆SLON1＝0.0＆CLAT＝0.0＆CLON＝0.0＆CRAD＝0＆SUBMIT＝Submit＋Search.

［4］ Tormod K.On exploitation of small－aperture noress type arrays for enhanced P－wave detectability［J］.Bulletin of the Seismological Society of America，1989，79（3）：888－900.

［5］ Mykkeltveit S，Asteb?l K，Doornbos D J，et al.Seismic array configuration optimization［J］.Bulletin of the Seismological Society of America，1983，73（1）：173－186.

［6］ Harjes H P.Design and siting of a new regional seismic array in Central Europe［J］.Bulletin of the Seismological Society of America，1990，80（6B）：1801－1817.

［7］ Kulhanek O.Signal and noise coherence determination for the uppsala seismograph array station［J］.Pageoph，1973，109（1）：1653－1671.