郝春月，鄭 重

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

1 引 言

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

2 資 料

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

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

3 功率譜比的計(jì)算

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

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

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

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

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

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

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

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

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

4 2次爆炸事件的相關(guān)性分析

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

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

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

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

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

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

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

式中：Gxx（f）和Gyy（f）分別為2個(gè)隨機(jī)信號(hào)x（t）和y（t）的自功率譜密度；Gxy（f）為它們的互譜密度。

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

圖7 2次事件的相關(guān)值計(jì)算結(jié)果Fig.7 The correlation results of the two events

圖8 用來計(jì)算2次爆炸事件相干函數(shù)的數(shù)據(jù)波形（10s）與2次事件的相干函數(shù)結(jié)果（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 結(jié) 論

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

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

［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.