基于爆炸地震波分析3次爆炸事件的主要特征*

郝春月，郑 重

(中国地震局地球物理研究所，北京 100081)

2006年10月9日1时35分(世界时,UTC)，全球地震台网(GSN)记录到一次爆炸事件。根据美国地质调查局所属的国家地震信息中心提供的地震目录(PDE)，此次爆炸事件发生的位置为北纬41.294°、东经129.094°,震级为mb4.3， 深度为0 km[1]。该事件距离GSN的牡丹江地震台(MDJ)约3.34°，方位角为186.4°。2009年5月25日0时54分(UTC)，GSN记录到一次爆炸事件，位置为北纬41.303°、东经129.037°，震级为mb4.7，深度为0 km[2]。该爆炸事件发生的位置距离牡丹江地震台约3.34°，方位角为187.1°。 2013年2月12日2时57分(UTC)，GSN记录到一次爆炸事件，位置为北纬41.308°、东经129.076°，震级为mb5.1，深度为0 km[3]。该爆炸事件发生的位置距离牡丹江地震台约3.33°，方位角为186.7°。这3次爆炸事件均位于朝鲜东北部咸镜北道吉州郡附近的丰溪里试验场。W.Y.Kim等[4]曾对2006年爆炸事件和爆炸事件所在区域的地震事件进行过判据方面的研究，范娜等[5]对这3次爆炸事件的震级进行过分析。由于这3次爆炸事件发生的位置较近，震级逐渐递增，本文中主要对这3次爆炸事件的能量比、相关性和相对位置进行研究。本文中，对牡丹江地震台3个分向(垂直向Z/北南向N/东西向E)记录的3次爆炸事件进行了P波振幅比、体波和面波功率谱比、初动震相相关系数的计算与分析。为了更好地研究这3次爆炸事件的相对位置，还引进了GSN白家疃地震台(BJT)记录的波形数据。牡丹江地震台位于中国黑龙江省东南部，白家疃地震台位于北京西北部，同属于GSN。为方便起见，用E06、E09和E13分别表示2006年、2009年和2013年发生的3次爆炸事件。

1 资 料

3次爆炸事件均被牡丹江和白家疃地震台记录到，牡丹江地震台记录的3次爆炸事件都很清晰,信噪比较高,见图1(a)。白家疃地震台由于距离较远(震中距约1 000 km)，只清晰地记录到了E09和E13。对于E06，经过2～4 Hz滤波也能清晰地看到Sg震相，但Pg震相被淹没在噪声里,见图1(b)。

图1(a)显示了120 s的波形数据。每幅图上显示的4个重要震相分别为Pn(首波，即Moho面绕射纵波)、Pg(直达纵波)、Sg(直达横波)和瑞利波。3次爆炸事件的体波清晰，Z向体波最大振幅的比值如表1所示。由于震级的计算也是利用最大振幅，所以最大振幅的比值在一定程度上描述了3次爆炸事件爆发能量的关系。表1给出了最大体波振幅比和最大面波振幅比。根据表1，把最大体波振幅和最大面波振幅比进行平均，可得出E13最大振幅是E09最大振幅的2.3倍，而E13的最大振幅是E06最大振幅的10.1倍。

在白家疃地震台记录的3次爆炸事件中，面波不发育，E06的Pg震相不清晰，所以没有进行能量比和相关系数的计算，引进白家疃地震台数据的主要目的是对这3次爆炸事件的相对位置进行分析。

表1 牡丹江地震台Z向记录的3次爆炸事件的振幅比Table 1 Amplitude ratios between three explosion events recorded by Mudanjiang seismic station in vertical direction (Z)

图1 2个地震台垂直向(Z)记录到的3次爆炸事件的波形Fig.1 Waveforms of three explosion events recorded by two seismic stations in vertical direction (Z)

2 功率谱比的计算

为了了解3次爆炸事件在频率域的能量分布情况，对3次爆炸事件进行了功率谱的计算与分析。功率谱分析采用Welch平均周期图法。Welch平均周期图法是对直接法的改进，即把一长度为N的数据xN(n)分成L段(在分段时可允许每一段的数据有部分的重叠)，每一段的长度为M,分别求每一段的功率谱，然后加以平均。第i段的功率谱可由下式表示[6]：

(1)

(2)

在此，分别对2个时间窗进行了功率谱估计，一个为P波组的15 s时间窗，另一个为面波组的15 s时间窗。根据郝春月等[7]的研究结果可知，前2次事件体波主要能量集中在2～4 Hz，所以首先对3次事件的体波时间窗进行2～4 Hz频段的带通滤波，又由于瑞利波的主要频段集中在1 Hz以下[8]，而后，对面波组时间窗进行了1 Hz以下的低通滤波。计算功率谱所用的15 s P波组和面波组数据均为600个采样点(采样率40 s-1)，利用Hanning窗，窗长为256点，重叠128点。3次事件的功率谱计算完毕后，可以得出3次事件中每对事件的功率谱点对点的比值。图2给出了3次爆炸事件垂直向记录中体波和面波的功率谱和功率谱比值。从图2可以看出E06、E09、E13爆发能量是递增的，谱比值见表2。

图2 牡丹江地震台记录的3次爆炸事件Pg震相和瑞利波之间的功率谱与谱值比Fig.2 Power spectrum densities for the Pg and Rayleigh wave of the three explosion events and the ratios of them

计算功率谱是为了对3次爆炸事件释放的能量进行对比，并且牡丹江地震台的传递函数在3次事件之间没有变化，所以没有对波形进行去除仪器传递函数的计算。因为在计算3次事件的谱值比的过程中，会抵扣掉传递函数，所以省却了这个步骤。省掉这个步骤的结果

表2 E13与E06、E09 Pg和瑞利波的功率谱比Table 2 Power spectrum ratios of E13 to E06 and E09

是图2显示的3次事件功率谱的纵轴单位不是(m/s)2,也不是表示能量衰减所用的dB，而是量纲为一。

表2中分频段给出了3次爆炸事件的功率谱比值。2～4 Hz频段(Pg震相的优势频段)，E13与E06的 Pg震相功率谱比值P13/P06=13.9，而E13与E09的谱比值P13/P09=2.6。瑞利波的主要频段为1 Hz以下，在此频段，P13/P06= 13，P13/P09=2.7。功率谱代表了能量，把Pg与瑞利波震相的功率谱比平均，E13释放的能量是E06的13.5倍，是E09的2.7倍。

3 3次爆炸事件的相关性分析

为了探索3次事件的相似性，对牡丹江地震台记录的3次事件的垂直向、北南向和东西向数据进行了相关性分析。为了对相似性进行定量地描述，对3次事件的相关系数进行了计算。首先对波形进行2～4 Hz的带通滤波，然后，选用了3次爆炸头2 s的数据进行了相关系数的计算。选择相关系数计算的时间窗一般为1～2 s[9-10]，在此选择2 s，使波形窗能包含更多的震相信息。

2个时间序列的互相关函数可以表示为[11]：

(3)

根据郝春月等[7]的计算结果，知道前2次事件的主要能量集中在2～4 Hz，所以对这3次事件进行了该频段的滤波。再根据式(3)，计算得出牡丹江地震台三分向记录的3次爆炸事件波形相关系数的分析结果，由于篇幅限制，图片只给出垂直向结果，见图3。根据计算，得出E13与E06垂直向、北南向、东西向波形的最大相关系数分别为0.92、 0.88和0.92；E13与E09垂直向、北南向、东西向波形的最大相关系数为0.99、 0.99和0.98。可见E13与E09的相关系数较高，接近于1。E13与E06的相关系数较小，这也可能是由于E06的震级较小，波形信噪比较低造成的，因为低信噪比的波形计算会导致精度降低。3次事件的相关性分析表示，3次事件在它们的主要能量分布频段高度相关。这个结果表明3次爆炸事件发生的位置相距很近，根据经验[12]，应不超过5 km。

图3 E13与E06、E09垂直向记录的相关系数计算结果Fig.3 Calculated correlation coefficients between E13 and E06，E09 recorded in vertical direction

4 事件的相对位置

3次事件具有特殊性，对它们进行精确定位具有重要意义。

为了精确定位地震位置，一般在地震活动区建立很多台站，形成一个区域地震台网。台网中台站的监测能力能够覆盖该地区[13]。根据区域地震台网管理规定，一个区域地震台网的定位误差小于5 km时，一般定义为Ⅰ类精度，也就是最好的定位结果。3次爆炸事件发生位置的周围没有区域地震台网，而GSN的地震台站分布稀疏，其中距离3次爆炸事件位置最近的是牡丹江地震台，震中距为300多千米。在台站稀疏并且没有区域台网的情况下很难对3次爆炸事件进行精确定位，也就是利用常规方法很难达到Ⅰ类精度的定位。鉴于确定绝对位置存在困难，本文中准备研究其相对位置。弄清其相对位置，对判断该地区以后爆炸事件的发生位置、侦破该地区的地质构造都具有重要意义。

根据3次爆炸事件的相关性分析, 判断3次爆炸事件组成的区域小于5 km, 而本文中要研究3次爆炸事件的相对位置, 就是要求精度比I类高。 由于E06震级较小，初动不清晰，初动读数误差将使定位精度大大降低，从而不能达到超越I类精度的要求，所以只针对E09和E13，计算其相对位置。

众所周知，假设虚波速度不变，震中距越远，则Sg与Pg的到时差tSg-tPg和Pg与Pn的到时差tPg-tPn就越大[14]。牡丹江地震台记录的波形出现了Pn震相(图1(a))，Pg震相振幅最大。把牡丹江台记录的2次事件的Pn震相两两对齐，最大的Pg震相进行对比,见图4(a)。白家疃地震台距离事件位置1 000多千米,最明显的震相是Pg和Sg(图1(b))，把白家疃地震台记录的2次事件的Pg震相两两对齐，Sg震相的最大振幅进行对比，见图4(b)。

图4 牡丹江与白家疃地震台记录的E13、E09的相对位置分析Fig.4 Analysis of the relative location between E13 and E09 recorded by MDJ and BJT

图5 台站与事件的相对位置Fig.5 Distribution of stations and explosions

图4(a)给出了牡丹江台记录的E13与 E09 Pg震相的波形对比。可以看出，E13的Pg震相比E09的晚到，由于牡丹江地震台位于爆炸事件位置的北偏东方向(图5)，可以得出E13事件在E09事件的南侧；图4(b)给出了白家疃台记录的E13与E09间Sg震相的波形对比。可以看出，E13的Sg震相比E09的晚到，由于白家疃地震台位于爆炸事件位置的西偏南方向(图5)，可以得出E13在E09的东侧。根据牡丹江、白家疃和2次爆炸事件的地理方位，可以判断，E13位于E09的东南方向，见图6。图6只是示意图，并不表示事件的真实位置。

本文结果与美国NEIC给出的PDE结果(图6)不同，尤其是E13事件，牡丹江地震台的波形很清晰，在Pn对齐后，E13的Pg震相比E09的晚到，即E13的震中距比E09的震中距大，又由于牡丹江地震台位于爆炸事件的北偏东方向，所以E13位于E09的偏南方向。

图6 3次爆炸事件的相对位置Fig.6 Relative location of three explosions

5 结 论

2013年爆炸事件(E13)爆发能量是2009年爆炸事件(E09)爆发能量的2倍多，是2006年爆炸事件(E06)爆炸能量的10多倍。

在2～4 Hz频段内，E13与E06和E09三分向波形的最大相关系数平均分别为0.90和0.99。分析3次爆炸事件的相关性可知，在2～4 Hz频段，E13与E06和E09高度相关。这表明这3次事件发生位置相距很近,根据经验，总孔径不超过5 km。

E13发生位置位于E09的东南方向,该结论与2013年2月5日郑州晚报的报道[15]相符，该报道暗示2013年爆炸事件发生在丰溪里试验场西侧坑道(2009年爆炸事件发生地)东南向的南侧坑道。

[1] USGS.Earthquakes: Search EQ archives[EB/OL].[2006-10-09].http:∥earthquake.usgs.gov/earthquakes/eqarchives/epic/ [2] USGS.Earthquakes: Search EQ archives[EB/OL].[2009-5-25].http:∥earthquake.usgs.gov/earthquakes/eqarchives/epic/ [3] USGS.Earthquakes: Search EQ archives[EB/OL].[2013-02-12].http:∥earthquake.usgs.gov/earthquakes/eqarchives/epic/ [4] Kim W Y, Richards P G.North Korean nuclear test: Seismic discrimination at low yield[J].Eos, Transactions American Geophysical Union， 2007,88(13/14):157-161.

[5] 范娜,赵连锋,谢小碧,等.朝鲜核爆的Rayleigh波震级测量[J].地球物理学报,2013,56(3):906-915.Fan Na, Zhao Lian-feng, Xie Xiao-bi, et al.Measurement of Rayleigh-wave magnitudes for North Korean nuclear tests[J].Chinese Journal of Geophysics, 2013,56(3):906-915.

[6] 胡广书.数字信号处理:理论、算法与实现[M].北京:清华大学出版社,2003:324-341.

[7] 郝春月,郑重.2次爆炸事件的相关性与能量比研究[J].爆炸与冲击,2010,30(5):535-540.Hao Chun-yue, Zheng Zhong.Relativity and energy ratio between two explosion events[J].Explosion and Shock Waves, 2010,30(5):535-540.

[8] 时振梁,张少泉,赵荣国,等.地震工作手册[M].北京:地震出版社,1992:41-47.

[10] Mykkeltveit S, Åstebøl K, Doornbos D J, et al.Seismic array configuration optimization[J].Bulletin of the Seismological Society of America, 1983,73(1):173-186.

[11] 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.

[12] Ingate S F, Husebye E S, Christoffersson A.Regional arrays and optimum data processing schemes[J].Bulletin of the Seismological Society of America, 1985,75(4):1155-1177.

[13] 彼得·鲍曼.新地震观测实践手册[M].中国地震局监测预报司,译.北京:地震出版社,2006:321-367.

[14] 傅淑芳,刘宝诚,李文艺.地震学教程[M].北京:地震出版社,1980:137-153.

[15] 郑州晚报数字报:A12版:国际新闻[N/OL].[2013-02-05]http:∥zzwb.zynews.com/html/2013-02/05/content_440861.htm