全球大震和中国及邻区中强震地震活动(1998年12月～1999年1月)s

陈培善

引言

除地壳岩石破裂引发的地震之外，地下核试验、工业开采、水库蓄水等人为因素也会引起明显的地面震动。在这些非天然地震中，人工爆炸的波形与地震极为相似，难以有效区分。随着人类工业活动日益频繁，矿爆也随之增加，仅北京及邻近地区每年就会有数百起矿爆发生。如果这些事件未被判识而被纳入台网目录中，无疑会导致虚假的断裂分布以及错误的活动性分析，对我国地震学研究造成极大的干扰。因此，发展精准有效的震源性质判别方法是地震监测任务的重要组成部分。

针对天然地震与人工爆炸的判别研究始于地下核实验监测任务（Argo et al，1995；谢小碧，赵连峰，2018；林鑫等，2019）。从理论上讲：爆炸属于膨胀源，震源较浅，波形表现出持续时间短、横波能量弱等特点；而天然地震属于剪切破裂源，震源较深，横波发育显著。以往研究试图通过波形记录中振幅、周期以及相位的差异来直接区分地震与爆炸（郑治真等，1975；Murphy et al，1997；Bonner，2006；Horasan et al，2009；Selby et al，2012），然而，受噪声干扰及波形叠加的影响，在时间域得到的判据效果并不理想。在实际监测中，傅里叶变换可将地震波形映射到频率域，在分离有效信号的同时，也可提取出更多的震源信息，是识别地震与爆炸的常规手段（Brune，1970；Randall，1973；Fisk，2007）。频域判据主要包括倒谱（Baumgardt，Ziegler，1988）、拐角频率（Dargahi-Noubary，1998；张萍等，2009），以及研究最为深入的纵横波谱值比（Walter et al，1995；Taylor，2011；Zhao et al，2016）。在核试验监测的地震学研究中，谱比占有举足轻重的地位。例如：Richards和Kim （2007）利用区域横纵波振幅谱比值将2006年朝鲜首次核试验判别为爆炸；Chun等（2011）比较了不同台站记录到的朝鲜2006年和2009年两次事件与中朝边界地震的Pg与Lg振幅谱比，结果显示核爆炸谱比值在3—1 Hz频段高于天然地震，区分效果显著；Zhao等（2009）和Zhao （2017）利用我国东北地区丰富的区域地震记录，经过振幅-频率-距离校正后在全台网内进行叠加处理，减少了单台观测所造成的数据起伏，结果显示在高于2 Hz的频带，P波与S波谱比值能够将所有朝鲜地下核试验从地震事件中识别出来。为了有效利用振幅比，需要处理传播路径和震源机制的影响，对于小震级工业爆炸而言，震源机制复杂，区域高频衰减模型对振幅比的影响较大，因此为了提高小当量爆炸的识别率，我们还需要挖掘更多的震源信息。

This page contains the following errors:

error on line 1 at column 1: Start tag expected, '<' not found

Below is a rendering of the page up to the first error.

传统的时频分析方法包括两大类：一类是以短时傅里叶变换和小波变换为代表的线性时频分析方法（Daubechies，1992；高静怀等，1997）；另一类是以维格纳为代表的Cohen类二次时频表征（Hlawatsch et al，1995；Qian，Chen，2002；Sejdi et al，2009）。短时傅里叶变换和维格纳分布发展较早（Mark，1970），也是最为成熟的两种算法，但均存在一定缺陷：前者受窗函数长短影响，无法兼顾时频分辨率；后者虽然分辨率高，但交叉项始终难以消除。1982年，法国科学家Morlet参考地震子波，首次提出了能够自适应调节时窗的小波变换（Morlet et al，1982a，b），此算法以其良好的时频聚焦性曾被赋予“数学显微镜”的美称。小波变换不是完美的，虽然其核函数可以根据频率伸缩，但缺乏相位信息，而且在实际应用中，小波时频谱不能有效地刻画信号的高频展布。之后，Stockwell等（1996）对Morlet小波进行改进，通过引入相移因子得到一种新的时频分析方法即S变换。S变换集短时傅里叶与小波变换优点于一身，其核函数与频率直接相关，同样具备多尺度特征，在刻画细节能力上强于小波变换。由于S变换的基本小波固定不变，为了更好地适应具体信号的分析与处理，国内外研究人员对其进行改进，提出了广义S变换（generalized S-transform，缩写为GST）理论。如今，该方法已广泛应用于多个领域（高静怀等，2003；Pinnegar，Mansinha，2003；Gholami，2013；黄捍东等，2014），也为震源性质判定研究提供了全新的思路。

本文拟应用广义S变换理论，围绕三河地区的地震和矿爆事件展开研究。首先对不同震中距的爆炸和地震信号进行时频分析；然后从视觉角度出发，将所有时频谱统一保存成规格一致的灰度图像，以避免震中距不同所带来的影响；最后使用图像复杂度的计算方法，将地震与爆炸的时频差异性量化，得到新的判据—“谱图二阶矩”，并通过实际资料验证该方法的可靠性。

1. 数据

本文使用的波形数据和事件目录源于北京遥测数字地震台网和国家数字测震台网中心（Zheng et al，2010）。试验数据选自2010—2016年期间发生在三河地区（39.8ºN—40.3ºN，116.8ºE—117.4ºE）已有明确分类的52次地震事件（表1）和52次矿爆事件（表2），震级M_L1.5—3.0，震中距处于5—260 km之间。每个事件均被研究区内的数十个观测台站记录到，对应多组三分量波形记录，其发震位置和台站分布如图1所示。

图 1 本文所用地震和爆炸事件及台站分布

Figure 1. Distribution of earthquakes，explosions and stations used in this study

下载: 全尺寸图片幻灯片

表 1 天然地震事件目录

Table 1. Catalogue of earthquakes

编号	发震时间（UTC+8）	北纬/º	东经/º	M_L	编号	发震时间（UTC+8）	北纬/º	东经/º	M_L
编号	年-月-日　时:分:秒	北纬/º	东经/º	M_L	编号	年-月-日　时:分:秒	北纬/º	东经/º	M_L
1	2010-02-23　23:14:48.0	40.16	116.94	1.9	14	2012-07-12　20:05:46.5	39.94	116.94	1.5
2	2010-03-28　10:20:56.3	40.03	117.02	1.7	15	2012-07-30　18:33:08.1	40.12	117.27	2.1
3	2010-04-23　01:35:50.8	39.99	116.95	1.6	16	2012-08-26　13:17:51.2	40.12	117.27	2.3
4	2010-06-02　15:57:36.0	40.06	117.01	1.6	17	2013-02-20　12:26:25.3	39.96	116.72	1.9
5	2010-06-22　03:15:27.1	39.81	117.32	1.8	18	2013-02-23　12:21:47.8	40.09	116.98	1.8
6	2010-11-23　14:38:06.1	40.08	117.06	1.9	19	2013-06-02　02:59:08.2	39.81	117.36	1.5
7	2011-02-05　17:06:16.4	39.89	116.93	1.7	20	2013-11-07　16:03:49.9	40.10	117.39	1.8
8	2011-02-19　11:42:45.0	39.99	116.96	1.9	21	2013-11-12　01:34:22.4	39.98	117.33	1.8
9	2011-08-18　04:31:50.7	39.86	116.73	1.5	22	2013-12-12　20:23:14.2	39.97	117.30	1.5
10	2012-01-18　23:03:54.6	39.80	117.35	2.0	23	2014-04-17　02:41:38.3	39.83	116.78	2.1
11	2012-02-22　14:40:25.2	40.04	117.02	1.5	24	2014-04-23　17:51:25.3	39.90	117.38	1.5
12	2012-03-22　20:59:07.0	39.97	116.95	1.8	25	2014-05-04　02:29:29.3	39.97	116.94	2.6
13	2012-05-25　13:31:29.2	40.19	116.98	1.9	26	2014-05-20　16:04:38.3	40.20	117.35	2.4
27	2014-08-29　23:17:41.2	40.02	116.96	1.6	40	2016-02-06　22:14:25.0	40.01	116.98	2.2
28	2014-12-21　04:44:17.1	39.92	117.23	2.8	41	2016-03-18　17:36:31.3	40.01	116.99	2.0
29	2014-12-21　16:55:42.0	39.81	116.77	1.6	42	2016-04-02　09:09:11.0	40.05	117.00	2.4
30	2015-01-27　21:58:49.2	40.03	116.96	2.2	43	2016-05-28　15:04:28.8	39.80	117.38	2.5
31	2015-03-10　09:24:51.0	39.92	116.91	1.6	44	2016-08-15　15:32:03.5	39.98	116.96	2.5
32	2015-04-08　17:42:35.9	39.87	116.85	2.3	45	2016-11-02　03:19:52.2	40.06	117.00	1.6
33	2015-04-12　18:54:21.2	40.12	117.24	1.5	46	2016-11-05　10:27:30.2	40.05	116.98	1.5
34	2015-04-14　17:30:50.4	39.82	116.80	2.6	47	2016-11-05　20:51:19.8	40.06	117.00	2.0
35	2015-04-19　08:47:25.6	39.86	116.86	2.4	48	2016-11-06　00:10:11.9	40.00	117.30	1.9
36	2015-04-23　18:39:10.0	39.82	116.81	1.5	49	2016-11-11　12:48:54.8	40.01	116.95	1.7
37	2015-05-25　02:00:48.6	40.14	117.09	1.9	50	2016-12-01　20:10:11.3	40.13	117.06	1.6
38	2015-10-02　11:30:53.3	39.83	117.10	2.3	51	2016-12-25　05:00:16.8	39.91	117.05	1.8
39	2016-01-16　10:52:19.3	40.05	117.04	2.1	52	2016-12-28　18:04:21.7	40.03	117.00	1.7

下载: 导出CSV

| 显示表格

表 2 人工爆炸事件目录

Table 2. Catalogues of artificial explosions

编号	发震时间（UTC+8）	北纬/º	东经/º	M_L	编号	发震时间（UTC+8）	北纬/º	东经/º	M_L
编号	年-月-日　时:分:秒	北纬/º	东经/º	M_L	编号	年-月-日　时:分:秒	北纬/º	东经/º	M_L
1	2010-01-06　14:39:36.2	40.04	117.15	1.9	27	2011-05-14　15:56:05.8	40.02	117.20	2.1
2	2010-01-09　16:25:37.4	40.03	117.13	2.0	28	2011-05-19　16:27:58.5	40.02	117.20	2.3
3	2010-01-13　14:34:53.0	39.78	117.05	2.5	29	2011-06-02　15:30:13.7	40.00	117.15	1.9
4	2010-01-13　14:44:23.1	40.04	117.20	1.9	30	2011-06-11　17:22:46.3	39.96	117.14	2.7
5	2010-01-21　14:36:13.4	40.02	117.12	2.1	31	2011-07-08　20:30:25.8	40.00	117.11	2.3
6	2010-01-22　15:40:26.6	40.05	117.14	2.0	32	2011-08-17　16:40:07.8	39.98	117.15	2.4
7	2010-01-26　13:09:38.6	39.99	117.10	2.5	33	2011-09-08　13:08:33.5	39.96	117.12	2.3
8	2010-02-06　12:59:52.2	40.07	117.15	2.2	34	2011-09-14　10:00:26.0	40.01	117.12	1.9
9	2010-02-06　13:17:33.2	40.03	117.16	1.9	35	2011-12-23　15:06:50.7	40.01	117.14	1.9
10	2010-02-06　16:23:30.3	40.04	117.16	2.4	36	2012-01-11　18:11:04.1	40.02	117.12	2.4
11	2010-03-02　17:40:07.5	40.02	117.18	2.3	37	2012-01-15　13:26:49.7	40.03	117.17	2.1
12	2010-03-03　14:57:27.2	40.05	117.19	1.8	38	2012-02-17　15:02:27.2	40.01	117.15	2.4
13	2010-03-27　18:05:28.5	40.02	117.23	2.4	39	2012-03-24　14:51:41.1	40.01	117.13	2.1
14	2010-03-29　13:31:11.4	40.01	117.08	2.2	40	2012-04-11　15:54:57.0	40.03	117.16	2.0
15	2010-04-01　16:20:22.0	40.05	117.16	2.0	41	2012-05-23　10:30:20.8	40.03	117.15	1.9
16	2010-04-04　14:46:13.5	40.03	117.10	2.0	42	2012-05-30　16:01:22.6	39.98	117.05	2.1
17	2010-04-21　11:12:40.2	40.03	117.19	2.1	43	2012-06-21　16:17:32.1	40.05	117.13	2.3
18	2010-05-06　14:48:13.5	40.04	117.13	2.1	44	2012-07-07　15:56:03.0	40.04	117.17	2.4
19	2010-05-06　15:42:23.1	40.03	117.15	1.9	45	2012-08-10　12:33:46.1	40.02	117.14	2.2
20	2010-05-08　14:57:52.6	39.98	117.09	2.3	46	2012-09-06　11:50:11.7	40.05	117.25	2.1
21	2010-05-09　17:41:13.0	40.02	117.14	2.1	47	2012-09-16　15:56:34.9	40.06	117.13	2.3
22	2010-05-13　11:43:29.8	40.02	117.12	2.2	48	2012-09-24　09:49:40.0	40.25	117.35	2.2
23	2010-05-19　17:44:27.6	40.04	117.15	2.1	49	2012-12-02　10:58:25.2	40.00	117.13	2.0
24	2010-07-24　20:06:58.3	40.07	117.12	2.1	50	2013-05-31　17:42:06.4	40.06	117.12	2.2
25	2011-05-08　14:29:44.7	40.04	117.14	2.4	51	2015-06-24　18:21:27.8	40.06	117.12	2.3
26	2011-05-09　16:46:31.4	40.00	117.15	2.3	52	2016-07-06　02:30:25.2	40.03	117.59	1.8

下载: 导出CSV

| 显示表格

对于小震、微震事件，信噪比是决定试验成败的关键因素。本文沿用短长时窗比对单台数据进行扫描，剔除信噪比很低的波形数据（隗永刚等，2019），且所选事件保留至少5个以上台站的记录，这样最终共得到地震波形记录617条，爆炸波形665条。对数据进行去倾向、去均值等标准化处理之后，利用广义S变换对垂直分量进行时频分析，开展下一步研究。

2. 方法原理

2.1 S变换

对于非平稳信号x（t），时频分析的基本思想是将其与指定窗函数相乘，截取小段信号进行傅里叶变换求取局部谱，窗函数沿着时间轴平移即可获得信号频率随时间的变化函数，即时频谱，其表达式为

$FT(\tau{\text{,}}f) {\text{＝}} \int_{ {\text{－}} \infty }^\infty {x(t)} {{w}}(t {\text{－}} \tau){{\rm{exp}}({ {\text{－}} {\rm{i}}2\pi ft}}){\rm{d}}t{\text{，}}$

(1)

式中，f和τ分别表示频率和时移，w（t）为时窗函数。时频分析的性质取决于时窗函数的选取，当w（t）为常用的汉明窗或高斯窗时，式（1）为短时傅里叶变换（short-time Fourier transform，缩写为STFT）。由于时窗固定，短时傅里叶变换无法调节时间和频率的分辨率：当时窗长度较大时，频率分辨率较高而时间分辨率较低；反之，时窗较短时，则时间分辨率高而频率分辨率低。为了克服短时傅里叶变换这一不足，Stockwell （1996）将时窗函数改为宽度随频率变化的高斯函数，即

${h_f}(t) {\text{＝}} \frac{{\left| f \right|}}{{\sqrt {2{\rm{\pi }}} }}{\rm{exp}}\left({ {\text{－}} \frac{{{f^2}{t^2}}}{2}} \right){\text{．}}$

(2)

将式（2）替换式（1）中的时窗函数w（t）即为S变换。S变换窗函数通过频率伸缩，使得时频分辨率在谱分解过程中同时达到最佳，这一特性与小波变换相似，效果远好于短时傅里叶变换。

2.2 广义S变换

This page contains the following errors:

error on line 1 at column 1: Start tag expected, '<' not found

Below is a rendering of the page up to the first error.

${{h}_{{g}}}(t) {\text{＝}} \frac{{\lambda \left| f \right|^{p}}}{{\sqrt {2\pi } }}{{\rm{exp}}\left({{ {\text{－}} \frac{{{\lambda ^2}{f^{2p}}{{t }^2}}}{2}}} \right){\text{．}}}$

(3)

将式（3）带入式（1）中得到广义S变换表达式（陈学华等，2008），即

$GST(\tau{\text{，}}\!\!\!\!f) {\text{＝}} \frac{\lambda }{{\sqrt {2\pi } }}\int_{ {\text{－}} \infty }^{ {\text{＋}} \infty } {{{\left| f \right|}^p}x(t)} {{\rm{exp}} \left[{{ {\text{－}} \frac{{{\lambda ^2}{f^{2p}}{{(t {\text{－}} \tau)}^2}}}{2}}}\right]}{{\rm{exp}} ({ {\text{－}} {\rm{i}}2{\rm{\pi }}f\tau }}){\rm{d}}t{\text{．}}$

(4)

广义S变换是目前最有效的时频分析方法之一，为了达到最佳的分辨效果，利用参数匹配法来确定最优化窗口参数，具体步骤如下：

This page contains the following errors:

error on line 1 at column 1: Start tag expected, '<' not found

Below is a rendering of the page up to the first error.

2）计算每一对λ和р的时频聚焦度（Radad et al，2015）

$ {{{M}}}(f{\text{，}}\!\!\!\!\lambda{\text{，}}\!\!\!\!p) {\text{＝}} \frac {{\displaystyle\int_{ {\text{－}} \infty }^{ {\text{＋}} \infty } {\left| {GST(t{\text{，}}\!\!\!\!f)} \right|} ^4}{\rm{d}}t}{ {\left[\displaystyle\int_{ {\text{－}} \infty }^{ {\text{＋}} \infty } {{{\left| {GST(t{\text{，}}\!\!\!\!f)} \right|}^2}{\rm{d}}t}\right]^2}}{\text{；}} $

(5)

This page contains the following errors:

error on line 1 at column 1: Start tag expected, '<' not found

Below is a rendering of the page up to the first error.

4）波形信号按照

$GST(t{\text{，}}\!\!\!\!f{\text{，}}\!\!\!\!\lambda{\text{，}}\!\!\!\!p) {\text{＝}} GST\{ t{\text{，}}\!\!\!\!f{\text{，}}\!\!\!\! {\rm{argmax}}[{{{M}}}({{f}}{\text{，}}\!\!\!\!\lambda{\text{，}}\!\!\!\!p)]\} $

(6)

计算，得到最终结果。

3. 数据测试

在现有数据中选取一条波形完整、高频信息丰富的地震信号（垂直分量）进行算法测试，所选事件发生于2014年5月4日，对应于表1中的25号事件，由承德地震台（CDT）记录。首先对波形数据进行带通滤波，去除0.5 Hz以下的长周期噪声；然后利用广义S变换，依据式（4）和式（6）对数据进行时频分析，并与常规的小波变换和S变换进行对比，结果如图2所示。

图 2 原始波形及使用不同方法计算得到的时频谱

（a）原始波形；（b）小波变换时频谱；（c） S变换时频谱；（d）广义S变换时频谱（参数最优化）

Figure 2. Time-frequency spectra for seismic signal calculated by different methods

（a） The original waveform；（b） The wavelet transform；（c） S-transform；（d） Generalized S transform；

下载: 全尺寸图片幻灯片

由图2可见：三种算法均能根据频率自适应调节窗口宽度，得到高清晰度时频谱，与短时傅里叶变换相比，图中看不到任何窗函数截取产生的痕迹；本次事件带宽约为2—25 Hz，能量最大值位于5 Hz附近，高频发育显著。对比图2b与图2c可见：小波变换能够清晰地勾画出5 Hz处的强能量团，但与S变换时频谱相比，10 Hz以上的高频信息能量较弱，不能描绘出完整地震的时频信息；由于S变换对高频具有放大作用，时频谱能够显示出更丰富的信息（图2c），但由于其频变时窗固定，高频部分的频率分辨率较低，可靠性有待验证。图2d是经过参数最优化搜索后的广义S变换时频谱，在保证高频信息捕捉的同时，结果更可靠，频率分辨率高于普通的S变换。需要说明的是，根据海森堡测不准原理，频率分辨率的提升会导致时间分辨率的下降，广义S变换得到的结果只是一个最佳的折中方案，但在本文的研究中却能换来置信度更高的判别效果。

为了验证算法准确性，对以上三种方法进行时间积分，退化为频域振幅谱，与常规的傅里叶变换结果进行比较，结果如图3所示，可见：S变换结果（蓝色曲线）的高频部分出现失真，曲线变化趋势与常规傅里叶谱相异，高频信息被明显抬高；小波变换振幅谱（黑色曲线）的低频部分与傅里叶变换吻合良好，但高频能量较弱，更高频段也出现了失真；广义S变换的时间积分（红色曲线）通过参数优化调整，全频段趋势与傅里叶振幅谱基本保持一致，吻合度最佳。此外，从曲线分布也可以看出小波变换振幅谱的低频能量较强，8—20 Hz的高频能量较弱，导致时频谱图的高频信息被掩盖。

图 3 四种方法所获振幅谱曲线的对比

Figure 3. Comparison of amplitude spectrum curves obtained by four methods

下载: 全尺寸图片幻灯片

采用最佳聚焦准则计算窗口参数，对每一频率分量进行变换时，需在λ和р构成的二维矩阵中进行搜索，找到最佳的参数匹配。这种方法虽然结果准确，但运算量大幅增加。为了提高运行效率，我们在试验过程中浏览了每一频率分量的矩阵M（f，λ，р），发现多数频点的λ和 р 分布特征存在一定规律性。如图4所示，主频能量5 Hz处的M矩阵中，聚焦度高的红色能量团分布类似一条对数曲线，虽然极值点max ［M（f，λ，р）］为（0.2，0.7），但选择任何相加接近1的（λ，р）点也能取得不错的聚焦效果，此规律适合任意频带。因此，在实际应用中可以根据经验直接确定窗口参数，提高计算效率。

图 4 广义S变换主频聚焦度（5 Hz）

Figure 4. The energy concentration of generalized S transform （5 Hz）

下载: 全尺寸图片幻灯片

4. 实际应用

4.1 三河地区地震和爆炸时频特征分析

为了研究地震与爆炸在时频域的差异性，我们选取震源相近、震级相同的一组地震与爆炸数据进行对比。选取的地震事件发生在北京时间2012年8月26日13时，爆炸事件发生时刻为2010年3月2日17时，两事件的震级均为M_L2.3，其参数分别对应于表1中16号和表2中的11号。记录台站四座楼（SZL）、兴隆东（XLD）、迁西台（QIX）、上方山（SFS）代表不同的震中距。首先对数据进行滤波，去除长周期噪声干扰，利用广义S变换对波形数据的垂直分量进行时频分析；然后依据上一节内容，为了兼顾时频分辨率，给定参数λ＝0.8，p＝0.4，相应时频分析结果如图5和图6所示。可见地震与爆炸事件在时频谱中存在显著差异，主要表现为以下几点：

图 5 地震事件波形及相应时频谱

Figure 5. Waveforms of earthquakes and corresponding time-frequency spectra

下载: 全尺寸图片幻灯片

图 6 爆炸事件波形及相应时频谱

Figure 6. Waveforms of explosions and corresponding time-frequency spectra

下载: 全尺寸图片幻灯片

1）地震属于破裂源，震源机制较为复杂，高频能量丰富，在时频谱中表现出更宽广的频带范围以及“多峰”特征；爆炸的时频谱图频带较窄，主能量限于10 Hz以下，随着震中距的增加，图中“亮点”所占比重减小。

2）爆炸震源较为简单，发生在地表，更容易激发面波，时频谱图中主能量集中在低频部分，而地震信号主能量位于更高的频段；受面波影响，爆炸的延续时间相比地震更长，这与徐果明和周慧兰（1982）的结果相反，且震中距越小，持时差异越显著。

3）天然地震在时频谱图中表现出更多的高频背景噪声，图像复杂度更高。造成这种明显差异的原因除了震源因素外，也归功于广义S变换对高频信息的准确捕捉，这也是本文不使用小波变换的最主要原因所在。

4.2 时频特征判据提取与应用

受震源、深度以及衰减机制的影响，近场地震信号的时频图像在能量面积、复杂程度上高于人工爆破。时频谱是一个二维矩阵，元素冗余，无法直接进行震源性质判别，因此需要进行量化，将蕴含的信息提炼出来。传统的时频特征提取方法的可靠性差，将波形复杂度、谱比值等一维判据延伸到时频域时，识别率反而下降。此外，震中距是影响地震波持时的主要因素，不同震中距会导致时频谱在时间轴上长短不一，也为时频特征量化带来困难。

对此，我们将广义S变换理论与图像处理中复杂度计算方法相结合，从视觉角度出发，得到能够有效区分小级别地震与人工爆破的判据−“谱图二阶矩”。首先对不同台站接收到的地震与爆炸信号进行时频分析，得到高分辨率时频谱；然后将所有的时频谱转换为规格一致的灰度图像，以避开震中距的影响；最后，利用

$U {\text{＝}} \mathop \sum \limits_{i {\text{＝}} 1}^m \mathop \sum \limits_{j {\text{＝}} 1}^n {[f(i{\text{，}}\!\!\!\!j) {\text{－}} u]^2}$

(7)

计算谱图二阶矩。式中，m和n分别为时频谱图的行数和列数，f （i，j）是时频谱图在（i，j）点的灰度值，u是以（i，j）为中心的3×3领域像素均值。谱图二阶矩能够反映图像的均一程度，数值越小，图像越简单，相反则图像越复杂。为了测试新判据的识别率，围绕本文第二节方法处理后的全部波形数据开展研究。首先，利用“迪格尔-凯西（Teager-Kaise）”能量算子对每一道波形数据进行分割，截取不同震中距、不同事件类型的有效波段（赵刚，2017）；然后，对波形数据进行时频分析及特征提取，得到上千个计算结果。在计算结果中随机抽取四个事件进行展示。图7和图8为归一化后的地震和爆炸时频谱灰度图，每一行代表相同的地震事件不同台站的结果，同样选取有代表性的四个台站进行展示，一共32张灰度图，台站名称、震中距和计算得到的谱图二阶矩U均标注在右上角。转换成灰度图像后，信号持时的差异被消除，每个图像规格均为580×824个采样点，近震被“放大”，相对较远的地震记录被“缩小”。从图7和图8中可以看到：谱图二阶矩U可以有效地反映时频图的复杂程度，图像中的能量团所占面积越大，结构越“碎”，U值也就越高；相比人工爆破，地震事件的时频灰度图具有更高复杂度，最大U值超过300 （图7a）；爆炸事件除了个别U值高于地震（图8b，8c）以外，其余均在12以下。如果将所有台站的U值求和平均，效果会更显著。求取平均值后，地震事件的U值分别为37.86，104.76，32.52，38.64，爆炸事件的U值分别为9.71，11.16，13.25，9.90，明显低于地震事件的U值。

图 7 地震事件时频谱灰度图像

（a）表1中1号事件；（b）表1中5号事件；（c）表1中6号事件；（d）表1中7号事件

Figure 7. The grayscale images of time-frequency spectra maps for earthquakes

（a） The event No.1 in Table 1；（b） The event No.5 in Table 1；（c） The event No.6 in Table 1；（d） The event No.7 in Table 1

下载: 全尺寸图片幻灯片

图 8 爆炸事件时频谱灰度图像

（a）表2中8号事件；（b）表2中12号事件；（c）表2中24号事件；（d）表2中26号事件

Figure 8. The grayscale images of time-frequency spectra maps for explosions

（a） The event No.8 in Table 2；（b） The event No.12 in Table 2；（c） The event No.25 in Table 2；（d） The event No.26 in Table 2

下载: 全尺寸图片幻灯片

对不同震中距的地震与爆炸垂向波形数据进行分类，不同台站数据混合统计，测试识别效果。利用式（7）分别计算谱图二阶矩，得到该判据随震中距的分布，如图9所示。当U大于阈值时被判定为地震，小于阈值时定义为爆炸，本次计算得到的阈值为25.54。从图9中可以看出：地震记录的谱图二阶矩U值整体高于爆炸信号，并且这种差异随着震中距的增大而增大，误判（地震的U值低于阈值虚线）多集中在震中距100 km以内；爆炸信号的U值样本分布较为集中，多数在阈值之下，随震中距的变化并不明显，错误识别（爆炸样本U值高于阈值）较少，多发生在震中距60—160 km范围内。本次试验中，爆炸、地震的样本总数分别为665和617，正确识别数为638和542，总识别率为91.89%。

图 9 谱图二阶矩随震中距分布图

Figure 9. The distribution of the second-order moment of spectrum with epicentral distance

下载: 全尺寸图片幻灯片

为了提高准确度，将同一事件不同台站求取的U值叠加起来求平均，结果如图10所示。可见：多台平均之后，判据分布趋于稳定，地震的U值分布在40—150之间，而爆炸的数值位于40以下，本次计算得到的判别阈值为37.92，仅有一个地震事件的U均值落在虚线以下，识别准确率达到98.08%，验证了利用谱图二阶矩区分事件性质的可行性。

图 10 多台平均后识别效果图

Figure 10. The second-order moment calculated by averaging the criterion values from all stations

下载: 全尺寸图片幻灯片

5. 讨论与结论

本文首先通过与传统时频分析方法的对比，验证了广义S变换算法的可靠性；然后选择合适的时窗调节参数对三河采石场附近发生的52次矿爆和52次事件波形记录进行时频分析，研究小震级地震与矿爆在时频特征上的差异性；最后，结合图形处理中的复杂度求取方法得到新判据谱图二阶矩，并对比新判据在单台记录以及多台平均后的识别能力，得到的结论如下：

1）广义S变换可以更好地兼顾高低频信息，较小波变换和S变换更加可靠，新判据以此为基础，可以得到最佳的识别效果。基于广义S波变换所提供的爆炸与地震事件的属性判据的识别准确率高达98%；

2）受激发与传播机制影响，爆炸与天然地震信号存在诸多差异，以往判据提取除了需要常规预处理外，还要进行震相标定、震中距分类的复杂程序。基于广义S变换方法，我们总结了地震和爆炸的时频特征，跳过复杂的人为处理流程，将时频图保存为规格一致的灰度图像，直接从视觉角度出发提取判据；

3）对于单一台站记录而言，新判据的识别效果受震中距影响较大，误判多发生在100 km范围内。当震中距增大时，P波、S波到时差增加，地震时频谱中呈现出典型的双宽带能量团，震源与爆炸源本来的频谱分布特点更加凸显，识别效果得到提升。

需要说明的是，本文方法无法兼顾时频分辨率。关于时窗参数优化的相关方法很多，在接下来的研究中将尝试不同的聚焦公式，以求达到最佳的分辨效果。

期刊类型引用(5)

1.	王婷婷，边银菊，任梦依，杨千里，侯晓琳. 地震事件分类识别软件. 地震. 2024(02): 104-119 . 百度学术
2.	刘甜甜，王禄军，张晖. 机器学习在地震事件自动分类中的应用. 华北地震科学. 2024(03): 96-101 . 百度学术
3.	张娜，王霞. 基于S变换的山西地区不同地震事件频谱特征分析. 山西地震. 2022(01): 1-6+11 . 百度学术
4.	杨志鹏，陈秀清，张御阳，颜欢，陈碧洪，阮祥. 基于MATLAB的定点形变观测数据时频分析软件设计及应用研究. 震灾防御技术. 2022(01): 172-180 . 百度学术
5.	孟娟，张家声，李亚南. 基于改进EWT和LogitBoost集成分类器的地震事件分类识别算法. 地震工程学报. 2022(05): 1233-1242 . 百度学术