基于生成对抗网络与随机森林组合模型的地震与地脉动区分研究

刘赫奕， 宋晋东， 李山有

(1.中国地震局工程力学研究所 地震工程与工程振动重点实验室，哈尔滨 150080;2.地震灾害防治应急管理部重点实验室，哈尔滨 150080)

随着地震仪在全球范围的大规模部署，人类已经进入地震大数据时代。如何处理地震网络每天收集的大量数据是一个非常棘手的问题，特别是从质量参差不齐的连续波形记录中分离出地震和噪声。由于不同区域下不同类型地震事件的特征不同，地震事件的分类与筛选具有很大的地域性差异。目前这些工作仍主要靠低效的人工完成，因此在地震大数据的时代背景下，发展高效、高精度、普适性强的自动波形检测算法显得越来越重要。

传统的地震检测方法在推进地震数据自动检测方面取得了很大进展，如长短时窗方法(short term averaging/long term averaging,STA/LTA)[1]、自回归方法[2]和高阶统计量法[3]，但由于仅采用波形的部分特征设定阈值，所以到时拾取精度比人工处理结果低，误拾取率较高。根据波形相似性提出的波形自相关方法[4]、模板匹配方法[5]等，虽然可以检测出大量微震，检测灵敏度较高，但需要事先建立模板，且互相关计算量大，处理效率低，难以用于地震数据的实时处理。

随着海量地震台站的建设，人工智能已经引入地震事件识别的研究中。目前，人工智能在地震检测[6-8]中已经取得了很多优秀的成果。Kong等[9]利用人工神经网络方法识别出Myshake系统中10 km内98%的地震记录，从而高效准确地区分出地震事件与人类行为。江汶乡等[10]利用触发后3 s内的波形特征值训练决策树模型，建立了特征值和事件类型之间的映射关系，从而正确地区分强震记录和干扰信号，达到了提高预警系统可靠性的目的。基于不同地区发生的不同类别的地震，Perol等[11]提出了针对地震识别和定位的卷积神经网络，该方法成功在美国俄克拉何马地区检测出了超过原始编目17倍的地震事件。为了检测不同时长的地震事件，Wu等[12]设计了一种名为DeepDetect的级联卷积神经网络，其识别事件的准确率可以达到63.8%，显著超过模板匹配方法的准确率。相比于传统的地震检测手段，这些机器学习方法可以自动提取更多更接近数据本质的特征，从而得到更为准确的地震识别效果。

为了探究机器学习中无监督学习算法对地震识别的效用，本文提出了一种利用生成对抗网络(generative adversarial network，GAN)与随机森林(random forest，RF)组合模型区分地震记录与地脉动噪声的方法。GAN是一种无监督学习算法，与人工标注不同的是，它可以自动提取波形特征，简化训练过程。随机森林是一种以决策树为基学习器的二分类模型，它可以利用GAN提取出的波形特征进行分类，并在做出准确分类的基础上有效解决过拟合的问题。通过该组合模型，本文将复杂的地震识别问题转化为简单的二分类问题，从而避免了识别过程中可能会出现的问题，旨在探索无监督学习方法在地震监测与地震预警中的应用可行性。

1 数 据

对于训练集数据，本文选取了2007年10月—2017年9月共79 762条日本K-NET和KiK-net井上竖直方向的强震记录，其中包括39 097条K-NET记录和40 665条KiK-net井上记录，震级范围为3级≤M≤8级。首先，本文筛选出了信噪比大于5的记录(信噪比定义为P波到时后2 s与前2 s的峰值速度比)，然后利用马强等[13]的方法对记录进行P波到时自动拾取，并对拾取后的到时进行人工校核。接着，本文筛选出了到时后1 s标准差大于前1 s标准差的加速度记录。

在正式训练之前，本文对筛选出的加速度记录进行了如下的预处理：① 基线校正；② 2阶0.075 Hz的巴特沃斯高通滤波；③ 利用插值方法将数据统一为100 Hz的加速度记录；④ 截取P波到时前1 s至后3 s的加速度记录作为地震数据；⑤ 截取KiK-net井上竖直记录中起始点至P波到时点间的记录，并随机选取4 s波形作为地脉动数据；⑥ 将所有截取后的数据归一化(除以最大幅值的绝对值)，图1给出了经过预处理后的地震动和地脉动数据的实例。经过上述处理，本文共得到63 052段地震记录和81 118段地脉动噪声记录(时长均为4 s)，地震记录的台站及震中分布图如图2所示。根据地震记录的数量，本文随机选取了63 052段地脉动噪声，并随机选取1/6的训练集作为验证集，其震级震中距分布图如图3(a)所示。

(a) 地震动记录

(b) 地脉动数据

图1 文中所用4 s数据的示例(横坐标为0的点即为P波到时点，图中数据均已经过归一化处理。)

Fig.1 Examples of 4 s data used in this study (The point where the abscissa is 0 is the P arrival. The data in the figure has been normalized)

对于测试集数据，本文选取了独立于训练集之外的数据集，包含2020年1月—6月共7 774条K-NET竖直方向的强震记录，震级范围为3级≤M≤8级。经过与训练集相同的筛选和预处理方法，我们得到了5 378条地震记录和5 378条地脉动噪声数据，所以测试集共有10 756条数据记录，其震中和台站分布图如图2所示，震级震中距分布图如图3(b)所示。

(a) 本文所用训练集的台站分布图

(c) 本文所用测试集的台站分布图

(a) 训练集震级震中距分布图

(b) 测试集震级震中距分布图

2 方 法

2.1 训练GAN

GAN是Goodfellow等[14]提出的一种无监督机器学习算法，由生成器和判别器组成。当向模型中输入真实波形时，生成器会利用随机样本生成与真实数据相似的模拟波形，判别器负责区分真实波形和生成的模拟波形。图4显示了GAN训练的基本过程，图中粗点线为真实波形，实线代表生成波形，细点线代表判别器做出的决策，下方箭头是用于生成波形的随机样本。图4(a)表示训练刚开始时，生成波形与真实波形差异较大，判别器可以很轻松地将它们区分开；图4(b)和图4(c)表示训练时，生成波形与真实数据越来越接近，判别器渐渐混淆了两种波形；图4(d)表示经过多个步骤的训练，生成器与判别器的能力均已达到饱和点，判别器无法区分两种波形，此时训练结束。

(a)

(c)

本文中GAN的架构如图5(a)所示。生成器的输入为50维服从正态分布的向量，输出为4 s的生成波形。生成器主要包含两个全连接层(神经元数目分别为128和400)和两个Leaky ReLU(rectified linear unit)层。为了得到更好的训练结果，将生成器最后一层的激活函数设置为tanh函数。判别器的输入是4 s的波形，输出是输入波形为真实地震动的概率。如果输入的是真实波形，则输出为1；如果输入的是生成波形，则输出为0。判别器包含一组卷积池化层、两个全连接层(神经元数目分别为128和1)和一个Leaky ReLU函数层。为了保持训练的平衡，本文将判别器的学习率(0.000 2)设为生成器学习率(0.000 1)的2倍，且生成器每迭代1次，判别器迭代5次。本文将批量设置为64，所以生成器共有58 128个参数，判别器共有203 073个参数。本文中的参数设置遵循了其他训练GAN的先例[15]，且可以在验证集上获得相对稳定的结果。

(a)

2.2 区分P波与地脉动

因为GAN的输入只有地震动记录，所以它只能用于区分真实的地震动和生成的地震动，因此本文只利用GAN提取P波特征，并进一步利用RF分类器区分地震动与地脉动噪声，这样地脉动数据也可作为RF的输入参与组合模型的训练。RF是一种基于决策树的分类器[16]，它利用bootstrap sampling[17]方法多次重复抽样构成子集[18]，并融合多个子集决策树的结果确定原始数据集的测试结果。因为它具有良好的二分类性能，所以可以充分利用GAN提取的特征对地震动和地脉动噪声进行区分，从而得到更为准确的分类结果。

因此，当GAN训练结束后，本文去掉了判别器最后一个全连接层，并将其余部分(图5(a)中虚线框中的部分)组成一个P波特征提取器，如图5(b)所示。本文假设训练后的判别器已具有识别地震的能力，所以P波特征提取器可以提取出地震波的关键特征[19]。随后，我们将10 756条测试数据输入到P波特征提取器中提取波形特征，其输出的128维向量即为RF的输入数据。经过RF的训练，本文可以得到最终的分类结果。如果输入样本为地震动数据，那么输出为1；如果输入样本为地脉动数据，则输出为0。

训练时，RF中树的个数为100，树深为45，其余参数设置均遵循sklearn.ensemble.RandomForestClassifier中的初始设置[20]。除此之外，当RF的参数probability设置为True时，模型可以输出一个对分类结果“自信程度”的概率，我们可以通过它了解测试结果的可靠性。

3 结 果

经过训练，GAN与RF组合模型的分类性能如表1所示，其中的标签含义如下：① 真阳性，组合模型将真实的P波判定为P波；② 真阴性，组合模型将真实的地脉动判定为地脉动；③ 假阳性，组合模型将真实的地脉动判定为P波，即为地震预警中的误报；④ 假阴性，组合模型将真实的P波判定为地脉动，即为地震预警中的漏报。

从表1可以看出，组合模型可以识别出测试集中99.80%的P波和99.11%的地脉动，证明该方法可以准确地区分地震事件与地脉动噪声，是一种较为可靠的识别地震事件的方法。此外，因为表1中假阳性样本的比例为0.89%，表明该方法可以有效地降低由地脉动造成的误报数量。

表1 测试集中组合模型的分类性能

本文还启用了RF中的输出概率，从而了解组合模型分类结果的可靠性。输出概率是一个介于0和1的概率值，其中1表示真样本(即P波)，0表示假样本(即地脉动)。输出概率越接近于1，该样本为P波的可能性越大；反之，输出概率越接近于0，该样本为地脉动噪声的可能性越大。本文将测试集中的输出概率全部导出，并统计了输出概率在各间隔之间(假设概率间的间隔为0.1)的样本百分比，如图6所示。横坐标为P波的输出概率，纵坐标为测试数据的百分比，所以在0.9～1.0概率段的P波样本数越多，同时在0～0.1概率段的地脉动样本数也越多时，模型的分类性能越好。本文可以通过图中两端的百分比检验分类结果的可靠性，并确定组合模型阈值与分类结果的关系。由图6可以看出，测试集中输出概率为0～0.1的地脉动数据占97.38%，输出概率为0.9～1.0的地震波数据占94.51%，这表明绝大多数地震数据和地脉动数据的分类结果都具有较高的可信度，证明该组合模型可以准确地区分地震动与地脉动噪声，从而识别地震事件。需要注意的是，本文将分类阈值设定为0.1和0.9，但在预警系统的实际应用过程中也可以通过台站噪声的触发概率来调整此阈值。

图6 测试数据百分比与P波输出概率的关系(需要注意的是，地震动数据集中在0.9～1.0，地脉动数据集中在0～0.1。)

另外，本文分析了分类准确率与震级之间的关系，如图7所示。通过观察可知，所有震级段的准确率均超过了97%，且4～5级地震的准确率达到了100%，6级及6级以上地震的准确率稍差。所以，本文着重分析了6级及6级以上地震的分类效果。测试集共有267次6级≤M<7级的地震，34次M≥7级的地震，图8表示了这301次地震的震中距与震级的关系。由图7可知，6级≤M<7级和M≥7级地震的分类准确率分别为98.13%和97.06%，所以组合模型将其中的6次地震误识别为地脉动，其中包含5次6级≤M<7级地震，1次M≥7级地震。表2显示了这6次地震的地震名称、震级和震中距信息，图8中的叉表示了这6次地震震中距和震级之间的关系。由表2和图8可知，在这6次误识别的地震中，3次地震的震中距超过了200 km，且其它3次地震的震中距也超过了150 km。因为地震预警所研究的范围多为近场地震，所以这6次地震的误分类均在可允许的范围内，所以证明该方法能够较为准确地识别地震事件，减少由地脉动造成的误报，而且在大震的近场台站记录中也同样适用。

4 与传统方法的对比

为了凸显出本方法相比于传统方法的优越性，本文选取了使用范围最为广泛的STA/LTA方法与本方法进行对比，从而突出本方法在识别地震事件方面的优势。STA/LTA方法是Allen提出的一种能量方法，主要用于地震P波、S波初至时间的自动拾取，其主要原理是根据地震波形特征函数的长短时均比值等特征拾取初至[21]，目前广泛应用于地震信号的检测中，特别是对于地震弱信号。本文中采用了马强等研究的计算公式

图7 分类准确率与震级之间的关系(直方柱的颜色对应于该范围内测试数据的数量)

图8 测试集M≥6地震的震中距与震级的关系(图中点为真阳性地震，叉为假阴性地震)

表2 测试集中误识别为地脉动(假阴性样本)的地震信息

(1)

式中：i为当前时刻点；k1和k2为i前某一时刻点；且k2

根据式(1)，在P波到来后，短窗内的平均值STA变化快，长窗内的LTA变化稍缓，所以STA/LTA刻画了记录幅值的瞬时变化，使特征变得更为明显。如果设定合适的触发阈值，当STA/LTA值超过设定阈值，则认为P波触发。本文采取了马强等关于STA/LTA方法的参数设置，相关阈值如表3所示。

表3 本文中STA/LTA窗长及阈值参数选取

因为STA/LTA方法的计算过程需要较长时间(至少30 s)的地震动记录，所以本文并未对输入数据进行截取。换句话说，本文只对输入数据进行了前三步的预处理，随后将整条地震动记录直接输入到STA/LTA方法中。而且，STA/LTA方法无需经过训练过程，所以只需输入测试数据就可以得到最终的结果。需要注意的是，因为STA/LTA方法输出的是P波到时的具体时间，并不是分类的结果，所以本文将STA/LTA方法拾取的P波到时位于实际P波到时前1 s至后3 s内的记录视为地震识别成功(输出为1)的数据，其余位置均视为识别失败(输出为0)。

表4显示了组合方法和STA/LTA方法对地震的识别效果对比。在5 378条测试地震动数据中，STA/LTA方法可以准确地识别出4 100条记录，也就是说，STA/LTA方法可以识别出76.24%的地震动。而组合方法可以识别出99.80%的地震动，与传统方法相比，本文的方法对地震动的识别率提高了23.56个百分点，说明该方法是一种可靠的地震动识别方法，其准确率要优于传统的自动捡拾算法。由于数据量巨大，本文比较了两种方法在计算时间上的优劣，如表5所示。由表5可见，组合模型在训练集中所有记录的训练时间总和较长，但因为训练过程可以预先完成，所以训练时间并不直接影响测试过程中组合模型对单条记录的计算处理时间。换句话说，地震发生后，组合模型在测试集中单条记录的计算处理时间仅为0.000 93 s(可忽略不计 )，而STA/LTA算法的单条记录计算处理时间为0.15 s。所以，从识别效果和单条记录的计算处理时间来看，组合模型相比于传统方法可以有效地提升识别地震动的准确性和时效性，是一种较为可靠的地震动识别方法。

表4 测试集中STA/LTA方法与本文方法的识别准确率对比

表5 测试集中STA/LTA方法与本文方法的计算时间对比

5 结 论

为了准确区分地震动与地脉动噪声，本文利用日本K-NET和KiK-net井上地震训练GAN，并将判别器改进为P波特征提取器提取波形特征，结合RF识别地震事件。测试结果表明，本文的方法可以区分出99.80%的地震波和99.11%的地脉动噪声，相比于传统的长短时窗方法提升了23.56个百分点，证明该方法是一种较为可靠的识别地震事件的方法，且在大震的近场台站中也同样适用。与手动标注波形特征相比，该组合模型可以自动提取出表征地震动的向量，从而简化训练步骤，节约训练成本，提升训练效率。

本文采用GAN对比真实地震动和生成的模拟波形，提取出能够代表P波特征的概率分布，从而获得真实数据中的某种隐含特性。GAN的判别器可以度量这种相似性，并通过与生成器的博弈缩短两种波形的差距。但是本文并没有利用判别器直接区分地震动与地脉动噪声，是因为GAN的输入只有真实地震动，而不包括地脉动数据，所以它并没有真正接触过真实的地脉动噪声，区分结果不具有说服力。但与RF结合后，组合模型的输入包含了输入到RF中的地脉动数据。经过大量训练之后，RF会逐渐学习并区分GAN提取出的两种特征，最终区分出两种波形，从而正确地识别地震事件。

虽然本文的组合模型可以比较准确地区分地震波与地脉动噪声，但是在训练过程中仍存在一些不足之处。一方面，组合模型需要在预先提取波形信号特征用以训练模型，所以与传统方法相比需要较长的训练时间,但是这是在地震发生前预先训练完成的，在地震发生后，测试过程单条记录的计算处理时间可忽略不计，随着计算机性能的提高，训练时间的影响可能会减小，训练效率上可以有所提高；另一方面，生成器生成的模拟波形作为副产品，暂时还未得到具体的应用。未来，本文将会朝着优化现有模型或数据增强的方向发展，从而发掘本方法在地震学和地震预警中更大的应用潜力。

致谢

日本防灾科学技术研究所(NIED)为本研究提供了数据支持，所有数据均从日本防灾科学技术研究所官网下载(网址：http://www.kyoshin.bosai.go.jp/(最后登陆日期：2019年5月20日))，文中图件使用通用制图工具GMT(Genetic Mapping Tools)和Python中的matplotlib模块绘制。