2次爆炸事件的相关性与能量比研究＊

郝春月，郑 重

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

1 引 言

2006年10月9日1时35分，在北纬41.29°、东经129.09°、深度0km 处发生一起爆炸事件［1］，此次爆炸事件距离全球地震台网（GSN）的牡丹江（MDJ）地震台约3.35deg（372.49km），方位角为186.4°；2009年5月25日0时54分（UTC），在北纬41.3°、东经129.04°、深度0km 处又发生一起爆炸事件［2］，该事件的位置距离MDJ约3.34deg（371.37km），方位角为187.1°。这2次爆炸事件距离很近。为研究2次爆炸事件的能量比与相似性，本文中试图对MDJ地震台3个分向（Z、N、E）记录的2次事件进行波形记录单振幅比、功率谱比、相关值和相干函数的计算与分析。

2 资 料

2次事件均被MDJ地震台记录。MDJ地震台位于中国黑龙江省东南部，属于全球地震台网（GSN），本文中选用的数据均来自MDJ台。图1显示了MDJ台3个分向（即垂直向、北南向与东西向）在2006年10月9日与2009年5月25日记录的爆炸事件波形。

图1显示了25s的波形数据，横轴表示采样点，纵轴为振幅。每幅图上显示的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 的数据xN（n）分成L 段（在分段时可允许每一段的数据有部分的重叠），每一段的长度为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波组数据为1 200点（采样率为0.025Hz），利用的汉宁窗窗长为256点，重叠128点。最后对功率谱进行了5点平滑。另外S波和面波组的时间窗为200s，数据为8 000点，利用的汉宁窗窗长为1 024点，重叠512点。最后对功率谱进行了5点平滑。2次事件的功率谱计算完毕后，可以得出2次事件的功率谱点对点的比值。

由于计算功率谱是为了对2次事件进行对比，并且MDJ台的传递函数在2次事件之间没有变化，所以没有对波形进行去除仪器传递函数的计算。因为在计算2次事件的谱值比过程中，会抵销传递函数，所以忽略了这个步骤。忽略这个步骤的结果是，图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 波组震相在垂直向的谱比值比较稳定，平均值在4.6左右。水平向的谱比值在2Hz左右出现最大值，可达到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个方向的数据进行了相关性分析。首先给出2次事件的波形叠加图（图6），从图中可以看出，2次事件确实具有相似性。

为了对相似性有一个量的描述，对2次事件的互相关值进行计算。首先对波形进行1～4Hz和2～4Hz的带通滤波（根据图3，P 波组震相在前几秒的能量主要集中在1～4 Hz）。然后选用2次爆炸前2s的数据进行互相关值的计算。选择相关值计算的时间窗一般为1～2s［4－5］，在此选择2s，使波形窗能包含更多的震相信息。

2个时间序列的互相关函数可以由下式表示［6］

图6 2次事件的3个方向的波形记录叠加Fig.6 Overlap of the waveforms of the two events recorded

式中：x（i）和y（i）分别是2个测点x和y 记录的噪声时间序列，N 表示数据样本点的个数和表示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数据在2～4Hz的能量要高些，并且互相关值平均为0.92。而在1～4Hz的互相关值平均为0.88。

为了能够在频率域显示2次事件波形的相似性，对2次事件进行相干函数的计算。在频域中，经常用到可用功率谱表示的相干函数rxy（f），相干函数rxy（f）（又称凝聚函数）给出了2个随机信号在频率域的相似性，其中f 表示2个随机信号的频率。一般来说，相干函数的值在0～1之间变化。若2个信号完全不相关，则相干函数值为零。对于2个相同的信号，其相干函数为1。相干函数定义为［7］

式中：Gxx（f）和Gyy（f）分别为2个随机信号x（t）和y（t）的自功率谱密度；Gxy（f）为它们的互谱密度。

2次事件的波形经过1～4Hz带通滤波后，利用10s的P 组震相数据进行相干函数计算（图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次事件的相干函数均保持在0.8左右。

图7 2次事件的相关值计算结果Fig.7 The correlation results of the two events

图8 用来计算2次爆炸事件相干函数的数据波形（10s）与2次事件的相干函数结果（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波组的相干函数平均保持在0.8左右。

衷心感谢许忠淮研究员的指导，感谢张爽同志提供的数据。

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