引言

地震观测台网记录到的事件中包含天然地震事件和非天然地震事件（如爆破、塌陷、矿震、核试验等），二者均是地震学的研究对象。天然地震是地壳应力积累和释放的结果，携带了地球深部构造、介质和应力的信息。天然地震的活动性是地震危险性判定的重要依据，但很多情况下天然地震不容易辨识，特别是天然地震与非天然地震事件的区分，一直是监测工作中的重要问题（刘莎等，2012）。天然地震与爆破的波形具有相似性，而在震源性质（曾融生等，2000）、频谱特征等方面有一定的差异，监测实践中，事件类型的判别主要依赖监测工作人员的经验（张帆等，2016），因而受到主观性的影响，无统一标准，也会出现误判的情形。随着地震台网规模的增加，地震台网监测能力的增强，地震记录数据量也大幅增加，人工处理方式在处理效率上面临挑战，因此产生了自动化处理和快速识别的需求。

通过不同事件类型在时域、频域等方面的差异（杨选辉等，2005；刘莎等，2012；王婷婷等，2013），可以设计基于特征和判别规则的自动分类方法，这些方法的判据有较明确的物理意义，且规则相对简单。随着计算机技术的发展，机器学习算法也被应用在地震事件分类中，以数据驱动的方式训练模型参数，获得事件分类模型，如支持向量机、神经网络等（边银菊，2002；黄汉明等，2010；毕明霞等，2011；张博等，2014；范晓易等，2019；蔡杏辉等，2020）。传统机器学习应用的主要困难在于特征工程中有大量需要人工完成的工作，且特征提取过程是基于专家知识的降维过程，可能丢失原始数据蕴含的有价值信息。

深度学习是机器学习和人工智能研究的最新发展方向之一，在图像和语音识别、自然语言处理等方面已取得较大进步，也被推广应用到各个领域并取得成效。深度学习通过与人脑相似的自下而上的层级结构，将输入端的数据逐层次提取（余凯等，2013），建立了底层数据到高层语义结构的映射。相较传统机器学习，深度学习具有以下主要特点：在输入数据方面，深度学习对原始数据人工加工较少，保留了大部分原始信息；在模型结构方面，主要使用深层神经网络，具有较复杂的模型结构和较多的模型参数；在流程方面，深度学习往往采用端到端的解决方案，而传统机器学习将任务分割为几个模块，分别训练。深度学习在地震事件分类问题中获得了较好的效果（陈润航等，2018；隗永刚等，2019）。

内蒙古地区有丰富的人工爆破事件记录，本文拟基于内蒙古地区天然地震和爆破的波形资料构建试验数据集，通过短时傅里叶变换将事件波形转换为时频域的振幅谱图，作为卷积神经网络的输入，使用计算机视觉领域的卷积神经网络作为分类器实现地震事件分类，并对结果进行测试，以评估该算法的有效性和稳定性。

1.   数据

1.1   数据选取

本文采用的原始数据为2016年3月至2020年5月内蒙古地震台网记录的事件波形，由于其中天然地震事件较爆破事件多，为了使数据集平衡，仅选取M_L2.7以上的天然地震事件。最终采用的数据包含417个爆破事件的900条台站记录和519个地震事件的900条记录，共计1 800条单台事件记录，选取垂直向记录作为原始数据。图1给出了事件空间分布图，天然地震事件在空间上分布较均匀，爆破的空间分布相对集中。图2为爆破事件的震级频次统计，爆破事件震级分布在M_L1.5—3.5之间；天然地震的震级分布在M_L2.7—4.5之间。

图  1  地震和爆破事件的空间分布图

Figure  1.  Spatial distribution of natural earthquakes and explosions

下载: 全尺寸图片 幻灯片

图  2  天然地震和爆破的震级频次分布

Figure  2.  Magnitude-frequency distribution of natural earthquakes and explosions

下载: 全尺寸图片 幻灯片

1.2   数据预处理

本文所使用的原始数据为单台垂直向单通道地震记录，事件波形存储在seed格式的文件中，使用MSDP软件解压缩为sac格式，通过obspy读取后在python环境下分析。obspy是一个针对地震领域开发的python库，支持常用的地震数据文件格式的读写，可以实现文件转换功能，集成了许多地震学专有工具包，并提供了简单易用的通用接口（Beyreuther et al，2010）。

选用数据的震中距在200 km以内，对每条记录截取100 s记录，截取起始点为初至震相前1 s，包含事件完整信息。原始数据的采样率为100 Hz。在数据预处理中，对信号进行高通滤波，保留1 Hz以上的信号并过滤掉低频噪声，之后进行去倾斜和归一化，最后通过短时傅里叶变换，将其转换为时频域的振幅谱。为了使数据分布均匀，再进行取对数操作和归一化操作，取对数操作在保持数据相对关系的同时，可缩小数据尺度，使数据的变化相对平稳，同时也削弱了模型的共线性和方差性等问题。转换后最终的数据为50×100的矩阵，对应频域1—50 Hz，时域0—100 s的对数振幅谱。图3为数据预处理的过程。

图  3  数据预处理过程

Figure  3.  The process of data preprocessing

下载: 全尺寸图片 幻灯片

短时傅里叶变换（short-time Fourier transform）是一种常用的时频分析方法，其优点在于对信号具有一定程度的自适应性，可以还原输入信号在时频域的能量分布，对信号的低频和高频两端都有较为一致的局部描述，得到的非平稳信号的时频分布较为平滑（郝国成等，2019）。STFT是关于频率f和事件t的二元函数，即

式中：x（$ \tau $）为信号；$\tau $表示窗口$ {g}^{*}$中心坐标对应的时刻；t表示待分析时刻；$ {g}^{*} $为窗函数，本文选用汉明（Hamming）窗，以减少频谱泄露问题。窗长决定频谱图的时间分辨率和频率分辨率，窗长越长，截取的信号越长，变换后频率分辨率越高，时间分辨率越差。为了克服数据的关联性，相邻的窗部分重叠，本文选择的窗长为2 s，重叠1 s。

STFT在时域将信号截取为多段，对每一段作傅里叶变换，求出每一段的频率特性，组合后得到时间-频率域的二维图像信号（时频域的振幅谱），该图像反映了信号能量在时频域的分布特征。

1.3   爆破与天然地震时频特征对比

本文分析了天然地震与爆破在频域和时域能量分布的差异。图4为地震与爆破的原始波形和时频谱的对比，由图可见：爆破的尾波有较多低频信号，含有面波成分（图4a）；天然地震的频率较丰富，而爆破的能量集中在低频成分（图4b）。

图  4  天然地震与爆破的波形（a）及时频域对数振幅谱（b）对比

Figure  4.  Comparison of waveform （a） and time-frequency log amplitude spectrum （b） between earthquake and explosion

下载: 全尺寸图片 幻灯片

天然地震和人工爆破可以看作随机信号，频率分布服从一定的统计规律。功率谱密度作为识别判据可以放大地震和爆炸的频域差异（隗永刚等，2019）。本文分别选取天然地震和爆破各300条记录，对归一化后的信号使用周期图法获得功率谱曲线，周期图法将随机序列观测数据视为能量有限的序列，直接计算离散傅里叶变换，取其幅值的平方并除以观测点数，作为真实功率谱的估计（王凤瑛，张丽丽，2006）。图5为天然地震与爆破的功率谱曲线（经平滑）对比，爆破在低频部分能量比例较高，随着频率增加逐渐衰减，地震在低频部分能量比例比爆破低，逐渐升高，在5 Hz左右达到峰值，随后衰减，天然地震能量随频率升高而衰减的速度比爆破缓慢。

图  5  天然地震与爆破的功率谱对比

Figure  5.  Power spectrum comparison between natural earthquakes and explosions

下载: 全尺寸图片 幻灯片

由于天然地震震源深度较爆破更深，天然地震P波频率成分更复杂，爆破的震源较浅，一般在地表附近，其产生的振动的高频成分被松散地层波吸收，因此记录到的爆破P波的周期比较大（张帆，2006）。人工爆破的震源为膨胀源，其破裂机制较简单，频率成分也比较单一，能量快速释放，P波极为发育，当距震中距50 km以上时，会产生瑞雷面波（赵永等，1995）。

2.   模型

2.1   卷积神经网络基本原理

人工神经网络受物神经网络概念的启发，由神经元层组成，这些神经元层通过权重值相互连接。每个神经元从上一层接收输入值，计算这些值的加权总和，加上一个偏置项，经过一个非线性的激活函数的作用，并将该函数值输出到下一层，即

式中，x_l－1为前一层输出特征图输入，w^l为权重系数，b^l为偏置.

该过程模拟了生物神经元的行为。训练人工神经网络的过程包括调整权重和偏置项，以使输入信号与所需输出值一致。基于对输出误差进行量化的损失函数对权重和偏差的更新称为反向传播的过程。在此过程中，使用梯度下降算法来迭代更新权重和偏差以使量化损失函数最小化，每个完整数据传递和参数更新过程称为一个迭代周期。

卷积神经网络与标准神经网络的主要区别在于卷积层的使用。卷积神经网络的前半部分一般为堆叠的卷积层，其间插入池化层构成特征提取器，后半部分将全连接层作为分类器，构成一个端到端的网络模型。卷积层中的神经元从前一层的较小的连续区域接收输入，乘以在整个输入数据上共享权重的滤波器（卷积核），一组共享的滤波器对输入图像进行卷积操作后得到一组特征图（张号奎等，2018）。卷积层输出可以表示为

式中，${\boldsymbol{x}}_j^l$表示第l层卷积层第j个卷积核对应的特征向量，M_j表示当前神经元的接受域，${\boldsymbol{w}}_{ij}^l $表示第l层第j个卷积核的第i个加权系数，${\boldsymbol{b}}_j^l $表示第l层第j个卷积核对应的偏置系数，f是一个非线性激活函数。

相比全连接层，卷积层具有局部连接、权值共享和输入/输出数据结构化等特点，减少了参数数量和计算复杂度。池化层的操作类似于卷积层，只是不计算加权和，而是对每个区域的输出进行简单的处理，如取最大值（最大池化）或平均值（平均池化），这样减小了特征空间尺度，也保持特征对旋转的不变性，提高了网络的计算效率。

卷积神经网络的训练过程分为数据前向传播和误差反向传播两个阶段。在前向传播阶段，数据由低层向高层传播。在反向传播阶段，误差从高层向低层进行传播，并根据误差对各层进行权值更新。反向传播算法主要基于梯度下降方法，网络参数首先被初始化为随机值，然后通过梯度下降法向训练误差减小的方向调整，优化卷积核参数k、下采样层网络权重β、全连接层网络权重w和各层的偏置参数b等。

2.2   模型构建

本文主要利用地震与爆破在频域和时域的差异，将地震事件分类问题转化为图像识别问题，将事件信号经过短时傅里叶变换得到的时频域对数振幅谱图作为输入，使用卷积神经网络实现天然地震事件和爆破事件的分类。

经过多次试验，对卷积层数量、卷积核数量和全连接层单元数等超参数进行调整。对比不同模型的效果，我们最终选用一个20层的卷积神经网络，包括3个卷积层和2个全连接层，图6为神经网络的结构图示，表1给出了神经网络各层的参数。每个卷积层之后加入一个批量正则化层（batch normalization，缩写为BN）和一个池化层，卷积层使用修正线性单元（rectified linear unit，缩写为RELU）（式（4））。堆叠的卷积层用于特征提取，输入数据为50×100矩阵，经过第一个卷积层的16个卷积核作卷积，经原始输入转换为16个50×100的特征图（50×100×16），卷积层的输出经正则化后作为池化层的输入，池化层对16个特征图降采样，将相邻的四个元素合并为一个元素，经池化层处理后输出为16个25×50的特征图（25×50×16）。第一个池化层的输出作为第二个卷积层的输入，经过32个3×3×16的卷积核处理后转换为32个25×50的特征图（25×50×32），经过第二个池化层降采样处理后输出为12×25×32。第二个池化层的输出作为第三个卷积层的输入，经过卷积层和池化层后输出为6×12×64，合并为长度4 608的特征向量，这个向量提取了时频图的特征，作为全连接层的输入。BN层对中间层的输入（激活函数前）作标准化处理，使得输出服从均值为0，方差为1的正态分布，减少网络训练过程中的中间协变量迁移，使网络训练过程对参数初始值不再敏感（刘建伟等，2020），从而改善反向传播时训练速度变慢和低层神经网络的梯度消失等问题。两个全连接层实现了分类：第一个全连接层包含128个单元，之后设置一个批量正则化层，也采用RELU激活函数；第二个全连接层实现分类功能，采用Softmax激活函数 （柔性最大值传输函数），Softmax函数将第二个全连接层的输出映射成为（0,1）之间的值，两个输出节点之和为1，每个节点的输出可理解为输入信号为该类型的概率，概率最大结点作为分类结果。第一个全连接层前添加了一个丢弃层。丢弃层通过忽略一定比例的特征检测器（将一定比例的神经元的输出置零）减少过拟合现象（Krizhevsky et al，2017）。

图  6  卷积神经网络结构

Figure  6.  CNN network structure

下载: 全尺寸图片 幻灯片

表  1  模型的超参数设置

Table  1.  Super parameter setting of the model

序号	类型	参数	输出
1	输入层		50×100×1
2	卷积层	16个 3×3卷积核	50×100×16
3	批量正则化层		50×100×16
4	RELU激活层		50×100×16
5	最大池化层	2×2	25×50×16
6	卷积层	32个3×3×16卷积核	25×50×32
7	批量正则化层		25×50×32
8	RELU激活层		25×50×32
9	最大池化层	2×2	12×25×32
10	卷积层	64个3×3×32卷积核	12×25×64
11	批量正则化层		12×25×64
12	RELU激活层		12×25×64
13	最大池化层	2×2	6×12×64
14	Dropout层	比例50%	6×12×64
15	全连接层	单元数128	128
16	批量正则化层		128
17	Relu激活层		128
18	全连接层	单元数2	2
19	Softmax激活层		2
20	输出层

下载: 导出CSV 

| 显示表格

3.   训练和检验

本文实验的硬件环境：英特尔i7-9850H处理器，64 G内存，英伟达RTX3000显卡，显存6 G。操作系统为Windows10，使用基于python的anaconda科学计算环境，模型的搭建和训练在tensowflow2.0开源框架进行。从1 800个样本中随机选取1 500个样本作为训练集，300个样本作为测试集。训练集是用于模型拟合的数据样本，测试集用来评估模型的泛化能力。

本文选用交叉熵作为损失函数，交叉熵刻画的是实际输出（概率）与期望输出（概率）的距离，也就是交叉熵的值越小，两个概率分布就越接近，交叉熵H的算法为

式中，P为期望输出，Q为实际输出。

采用的训练算法为Adam算法（Kingma，Ba，2015），Adam优化算法是随机梯度下降算法（stochastic gradient descent）的扩展，可以更加有效地更新网络权重，使用动量和自适应学习率来加快收敛速度。

为了评估学习率对训练效果的影响，分别将学习率依次设置为10⁻³，10⁻⁴，···，10⁻⁷进行训练，设置迭代次数为880次，图7a为不同学习率的损失曲线（经平滑），图7b为不同学习率的准确率曲线，准确率随着学习率的降低而下降，学习率较高时训练曲线收敛较快，最终准确率也较高，因此我们选取10⁻³作为初始学习率，并设置分段降低学习率。

图  7  学习率测试

（a） 损失函数曲线； （b） 准确率曲线； （c） 准确率随学习率变化

Figure  7.  Learning rate test

（a） Loss curves ；（b） Accuracy curves ；（c） Accuracy rate changes with learning rate

下载: 全尺寸图片 幻灯片

交叉验证方法将训练数据集本身细分成不同的验证数据集去训练模型，验证集可以用于超参数的优化和对模型效果的初步测试。为了评估算法的性能和稳定性，用5折交叉验证检验，将1 500个样本均分为5组，每组300个样本，每轮选取一组作为验证集，将留下的4组作为训练集（图8）。5轮训练的验证集准确率可达95%—99.44%，验证集平均准确率为97.33%。

图  8  5折交叉验证结果

Figure  8.  Five fold cross validation results

下载: 全尺寸图片 幻灯片

将全部训练集样本用于模型参数的训练，使用测试集的300个样本评估网络的泛化能力，最终测试的准确率达到98.02%，模型有较好的泛化能力和准确率。

4.   讨论与结论

选取内蒙古地区2016年3月至2020年5月记录的天然地震和爆破事件中天然地震和爆破的单台垂直向记录各900条，经过预处理后转换成时频域的对数振幅谱分布图作为实验数据集，共包含1 800条样本。使用一个含有3个卷积层的卷积神经网络作为分类器，将实验数据集划分为训练集和测试集进行训练和检验。在训练集上使用5折交叉验证，结果显示，验证集的平均准确率达到97.3%。使用全部训练集数据对网络训练后，测试集的准确率达到98.02%。

相比于传统机器学习，本文采用的模型在特征提取方面仅进行了简单处理和时频分析，保留了信号在时域和频域分布特征较完整的信息，使用卷积层自动提取特征，避免了使用先验的降维规则，以数据驱动的方式得到分类模型。实验结果表明本文使用的基于短时傅里叶变换和卷积神经网络的分类算法有效提取了地震信号的能量分布的特征，在地震事件分类任务上有较好的效果，与已有研究对比，本文模型达到了较高的准确率。

在地震类型识别研究中，有多种数据处理方式和分类模型可供选择，在进一步研究中可以对比不同模型的效果。非天然事件包括爆破、塌陷和滑坡等多种类型，在现有数据条件下，本文仅研究了天然地震和爆破的分类问题，在积累了充足的多种类型数据样本后，可以开展更加细致的分类研究。

评审专家的建议对稿件的提升帮助很大，作者在此表示感谢。