引言

地震是一种极具破坏性的自然灾害，准确预测地震的发生目前还是一个尚未解决的科学难题（Geller，1997；Peng et al，2017），因此地震预警（earthquake early warning，缩写为EEW）作为一种新型防灾技术近年来愈加受到重视。地震预警的基本任务是利用震源区附近台站监测到的地震波早期信息估算出该地震的震级、位置、深度等震源参数，并快速评估出目标预警区域潜在的地震破坏，及时发布地震警报信息（Allen et al，2009；Minson et al，2018），目前已有多个国家开发出地震预警系统并成功运行（Wu，Teng，2002；Hoshiba et al，2008；Espinosa-Aranda et al，2009），我国也在2008年之后启动了“国家地震烈度速报与预警工程”项目建设，目前该项目已经进入最后的建设阶段（李山有，2018）。

地震预警系统中准确快速地估算震级大小是最核心的关键问题。地震震级的快速估算有多种方法，传统方法是从初始地震动信息中提取关键指标作为特征参数，然后利用特征参数与震级之间的经验公式来估算地震震级的大小，其中特征参数有很多种选择，可以利用初始地震动的峰值、周期（频率）、能量及其组合等信息来建立。基于初至P波的频域信息与震级之间比较强的相关性，提出了一系列特征频率类参数，如最大卓越周期$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{p}} $ （Allen，Kanamori，2003）、平均卓越周期τ_c （Wu，Zhao，2006）等。Wolfe （2006）认为利用$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{p}} $估算震级的稳定性值得怀疑，其研究结果显示该参数实际上是谱幅值和频率的非线性函数，因此在震级估计时$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{p}} $会受到振幅谱和相位谱的共同作用因而造成较大估算误差。Lockman和Allen （2005）等在日本Hi-net （high sensitivity seismograph network）台阵记录上使用$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{p}} $估计震级，其结果表明该参数受区域特性影响较大，在低密度台站地区需要联合多台站记录以提高精度，因此时效性得不到保障。王延伟等（2020）利用日本KiK-net台网井下基岩强震数据研究的结果显示，τ_c受加速度积分运算影响较大，该过程会引入长周期成分致使震级估算偏大。除特征频率类参数外，初至P波的位移（P_d）、速度（P_v）或加速度（P_a）的峰值也被证明与震级有关，也可用于震级估算，此类方法被称为特征幅值类参数法（Wu，Kanamor，2005；Zollo et al，2006；Melgar et al，2015）。这类方法中，P_d参数方法的结果表现最优（马强，2008），该方法受场地条件影响较小，对地区差异性不敏感，但当震级M≥6.5时，利用初始P波信息计算出的P_d参数与震级之间的关系存在震级饱和效应。$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{p}} $，$ {\tau }_{\mathrm{c}} $和$ {P}_{{\rm{d}}} $等基于特征参数的传统方法，需要提前定义特征参数，之后仍然需要依赖特征参数与震级之间的经验关系，所以在震级估算的精度和时效性方面难免存在弊端。因此，能否放弃人为定义的特征参数和震级估算经验公式，利用深度学习方法直接从地震记录数据中学习到震级的相关特征从而直接估算出地震震级成为一个值得期待的新的研究方向。

大数据的深度学习方法很多，其中卷积神经网络算法（convolutional neural network，缩写为CNN）因其能够自主进行标签特征的提取，再加上对复杂非线性行为具有良好的分析能力，被应用于地震预警领域的研究（于子叶等，2018；Dokht et al，2019；Zhu，Beroza，2019；Wang et al，2021）。近些年，CNN方法也被用来解决传统参数在震级估算中存在的问题，但在精度和时效性方面仍然存在较大的改进空间。Lomax等（2019）提出的ConvNetQuake_INGV深度卷积模型，将单台站记录的三分量加速度时程数据作为第一卷积层的输入数据，对三分量波段中的加速度峰值进行归一化处理，并将归一化值作为全连接层的辅助输入参数。利用上述两部分数据对模型进行迭代训练，以达到自动化检测并快速表征地震震级的目的，但其不足之处在于该方法截取波段为P波到时前5 s和到时后45 s的长时窗加速度时程数据，时效性较差，无法在实际地震预警应用中给公众预留足够的行动时间。此外，ConvNetQuake_INGV模型仅能对中等震级范围内的地震事件进行估算，误差较大，且该研究中的CNN模型会出现一定程度的过拟合现象。Saad等（2020）构建了一个CNN深度网络分类器，选择地震到达时间差小于6 s的三个台站，将各个台站P波到时前6 s和到时后2 s的加速度垂直分量数据联合组成3×800的二维矩阵，以此作为训练数据集对CNN模型进行训练，从而对地震事件的震级、深度、位置和发震时间进行分类预测。测试结果表明该分类器的分类准确性较高，但就震级而言因分类标准宽泛，仅能区分震级是否在M_W2.0—3.49，M_W3.5—4.5或M_W≥4.5，这样的分类标准容易导致震级预警误差偏大，减弱地震预警的效果，并且训练数据由三个台站联合构建，时效性上略逊于单台站震级估算法。Mousavi和Beroza （2020）尝试使用将CNN与长短时记忆神经网络（long short-term memory neural networks，缩写为LSTM）拼接的复合模型，利用单台站P波到时后30 s的三分量加速度时程信息进行震级估算，得到了较高的震级回归相关系数，但同样存在截取时窗尺度过长、时效性无法确保的问题。总之，上述几种利用卷积神经网络估算震级的方法在时效性和震级精度方面无法较好地满足地震预警的需求。

鉴于此，本文提出一种基于多层全连接卷积网络架构的深度学习方法，利用P波到达后3 s至9 s内的地震频域数据、震中距、震源深度和v_S30等信息直接进行震级估算，以提高地震预警中震级估算的精度和时效性，为未来利用多台站记录进一步提高震级估算精度提供参考。

1.   数据集设置

数据质量是决定人工智能神经网络模型能否学习到各输入特征之间隐式关系的关键因素，因此有必要预先对地震记录数据进行处理。

1.1   数据来源

本文所用数据来自KiK-net和K-NET台网，从日本防灾科学技术研究所的强震动数据库中获取（National Research Institute for Earth Science and Disaster Resilience，2020）。收集1997年10月8日至2019年3月11日期间的地表地震动记录，每组数据包括各台站对应的三分量地震波信息，此次研究仅使用垂直分量加速度数据。数据预处理前需对地震记录按照以下原则进行筛选：① 日本气象厅震级M_JMA≥4.0；② 震源距介于25 km与200 km之间，这样一方面能保证所筛选的地震记录至少有3 s初至P波信息，另一方面200 km的震源距上限可包含破坏性的日本近海地震（王延伟等，2020）；③ 三分量井上地震波合成后的地震动峰值加速度大于等于2 cm/s²；④ 信噪比（signal noise ration，简写为SNR）大于等于20 dB。筛选后的地震事件震级、所选测站的震源距及场地条件分布情况见图1，场地分类见表1。

图  1  震级、场地剪切波速及记录数随震源距的变化

Figure  1.  Distribution of magnitude，site shear wave velocity and event counts with hypocentral distance

下载: 全尺寸图片 幻灯片

表  1  场地条件分类（引自Zeynalov et al，2013）

Table  1.  Site condition classification （after Zeynalov et al，2013）

场地分类	场地分类	vS30/（m·s−1）
A	硬土	［ 600，＋∞）
B	中土	［ 200，600）
C	软土	（0，200）
注：vS30代表地表下30 m土体的平均剪切波速。

下载: 导出CSV 

| 显示表格

1.2   数据处理流程

经筛选后的井上地震动记录格式不统一，且含有冗杂信息，因此数据在输入神经网络前要经过净化流程处理，具体步骤如下：① 对加速度垂直分量数据进行基线校正，同时为确保全部台站在相同时间尺度下的采样点数一致，将采样频率统一至100 Hz；② 采用长短时平均方法（short term averaging/long term averaging，缩写为STA/LTA）结合赤池信息量准则（Akaike information criterion，缩写为AIC）震相到时估算法对加速度时程记录中的P波进行初步拾取；③ 人工选波并剔除P波拾取异常的地震记录；④ 对加速度记录时程进行快速傅里叶变换得到对应的频谱记录；⑤ 将频谱记录取对数后的结果与震中距、震源深度以及v_S30合并保存，作为模型分析的最终输入数据。

加速度时程数据经以上过程筛选并整理后，最终得到3 065次地震的16万4 547条垂直分量加速度记录的频谱信息。考虑到地震发生时间的序列性，选用1997年至2014年所记录到的12万4 955条加速度频谱信息作为模型训练数据集（不包含2011年3月11日地震记录），将2015年至2019年所记录到的3万9 592条加速度频谱信息作为模型测试数据集（包含2011年3月11日地震记录）。按照时间划分数据集，一方面确保训练集与测试集的地震数目大致符合深度学习领域的八二开原则，另一方面又确保地震事件的完全独立性，即同一地震事件的不同台站记录不会混杂于训练集和测试集中，具体分布见图2。另外，此次研究将2011年3月11日所获的24条特大型地震记录全部放进测试集，以测试本文所提出的模型在面对超灾难性地震事件发生时的表现。

图  2  训练集与测试集分布情况对比图

Figure  2.  Distribution of training and test datasets

下载: 全尺寸图片 幻灯片

2.   模型架构及训练超参数设置

本研究编译环境选择基于数据流图的机器学习框架tensorflow-gpu 1.13.1，编程语言是python 3.7。构建CNN模型时使用keras提供的函数式模型定义法，以便有效地避免传统序列模型缺乏自由度的问题。

2.1   模型架构

模型体系架构由输入区、卷积区、全连接区三部分组成，如图3所示。

图  3  CNN模型架构图

Figure  3.  Architecture of CNN model proposed in this study

下载: 全尺寸图片 幻灯片

1） 输入区。分为主输入区和辅助输入区。加速度垂直分量的频域信息作为主输入区的输入数据设置在卷积池化层前，震中距、震源深度和v_S30作为辅助输入参数附加到后置全连接层的第一层中。这样分区设置输入参数的优势在于频谱数据经过卷积池化层自动提取震级特征图后，可以结合全连接层引入的震中距、震源深度和v_S30等信息，使模型进行震级估算时将地震动在不同路径上的衰减效应考虑在内，从而获得更准确的结果。

2） 卷积区。前置八层为卷积池化层，最后一层为全连接层。前两个卷积层的滤波器设置为64个，后两个卷积层的滤波器设置为128个，四个卷积层的滤波器尺寸统一取2，移动步长统一取1。每个卷积层后都连接同参数池化（下采样）层，选用最大池化（max pooling）操作进行层间数据处理，并规定池化层滤波器尺寸和移动步长均取2。另外，选择相同填充（same padding）作为全部卷积池化层中的填充属性（padding），使数据经过某一隐含层前后的维度保持一致。全连接结构实际上是对最后一个池化层提取出的局部特征图进行展开（flatten）处理，实现数据一维化。这样将卷积池化层学习到的分布式特征整合并输出为特征值，空间结构特性则可忽略，从而减少特征位置对输出结果的影响。

3） 全连接区。设置4个全连接层，相应的神经元节点分别设置为259，64，16，1个。实际上在构建全连接区第一层时仅设置了256个节点，但为了将震中距、震源深度和v_S30这三个参数引入，额外设置了仅包含三个节点的全连接层与之拼接。整体采用λ＝0.001的L2正则化方法在损失函数中添加惩罚项，防止模型过拟合。卷积网络后拼接多层全连接网络的目的是为了更好地拟合特征值、辅助参数与输出目标之间的非线性关系。

在整体的网络设置中，除最后输出节点外，Relu非线性激活函数应用于每一个卷积和全连接层的输出。

2.2   模型超参数设置及训练过程说明

本文在配置训练方法时选择Adam随机梯度下降优化器（η＝0.001）并采用均方误差公式（1）作为模型的损失函数，平均绝对误差公式（2）作为模型性能评估指标。

式中，n为样本数，$\stackrel{\sim}{{Y}_{l}}$为预测震级，Y_i为真实震级，文中所用震级均采用M_JMA。

模型训练采取五折交叉验证的方式来实现模型内部参数的权重更新，具体过程为：将全部训练数据集随机划分成五个数量相等但不重复的子集，每次训练把其中一个子集作为另外四个子集训练的测试集，通过迭代训练完成神经网络中的权重值更新，流程如图4所示。每批输入2 048组数据，设置100轮训练并配备早停机制，当连续10轮验证损失不再下降时停止训练。整个过程所使用的算力来自于Nvidia GTX 2080Ti，每轮训练平均耗时3分钟。

图  4  模型训练五折交叉验证示意图

Figure  4.  Schematic diagram of five-fold cross-validation in model training

下载: 全尺寸图片 幻灯片

3.   结果

3.1   3 s固定时窗下的模型预测结果

确定时窗长度为3 s，根据已选时窗长对训练测试数据进行截取，获得300个样本数据点，截取后的训练数据用于卷积模型内部超参数的确定，测试集用于模型性能测定。鉴于测试集和训练集中的地震事件完全独立，并且测试数据集不参与任何训练过程，因此基于测试集的模型预测结果评估比较合理，具有代表性。

震级预测值与真实值之间关系如图5所示。假定容错范围为±0.5 （陈运泰，2009），那么容错范围内包含的样本为3万5 602个，占总测试集样本的比例为89.92%，误差分布如图6所示，可见：3 s时窗下估算震级的绝对平均误差为0.255 0，标准差为0.244 3。基于特征参数预测地震震级的相关研究（王延伟等，2020）表明，P_d参数的预警震级估计结果表现出的准确率总体优于其它参数，且在P_d与震级之间的经验公式的拟合过程中使用了震源距参数，与本文模型的特征参数选取具有共同因子，故本研究与利用基于P_d预警参数计算结果进行对比研究，以评估本文提出的模型在震级估算的精度上相较于传统特征参数算法是否有提升。如图6所示，利用$ {P}_{\mathrm{d}} $参数估算震级，有2万4 168个样本在误差允许区间内，所占比例为62.18%，对应误差分布如图6中黑边条带所示，3 s时窗下P_d估算震级的绝对平均误差为0.456 0，标准差为0.361 9。两者对比表明，本文提出的估算震级方法相较于利用参数P_d估算震级的方法，其震级估算误差小、离散度低，具有更高的震级预测精度。需要指出的是，两种方法在估算震级大于M_JMA6.5的地震事件时均存在震级低估现象，但本文模型的误差均值略小于P_d参数方法，对预测结果具有一定的改善，详细的原因待下文探讨。

图  5  3 s时窗下CNN估算震级（a）、P_d估算震级（b）与实际震级的比较分析

Figure  5.  Comparison of the estimated magnitude by CNN （a） and by P_d （b） with actual magnitude for the time window of 3 s

下载: 全尺寸图片 幻灯片

图  6  3 s时窗下震级估算的误差分布（a）和误差均值图（b）

Figure  6.  Error distribution （a） and error mean （b） of CNN and P d magnitude estimation with 3 s time window

下载: 全尺寸图片 幻灯片

图  7  不同时间窗长t下CNN （上）与P_d （下）估算震级与实际震级的比较分析

Figure  7.  Comparison of the estimated magnitude by CNN （upper） and P_d （lower） with actual ones for different time windows

下载: 全尺寸图片 幻灯片

3.2   不同时窗下的模型能力评估

在实际的地震预警应用中，台站采集样本点的数据量随时窗增长而增加。为了深入探究本文模型是否能够随可用数据的增加而提高地震震级的预测精度，本文将P波触发后4—9 s各个时窗长度下的震级预估结果与同条件下的P_d参数方法结果进行对比分析，具体数据列于表2，其中5 s，7 s和9 s的样本预测结果与实际震级的对比结果如图7所示。

表  2  神经网络（CNN）模型和P_d参数方法在不同时窗下的预测结果统计

Table  2.  Prediction result statistics of neural network model （CNN） and P_d with different time windows

时窗长度/s	误差标准差			误差平均值			准确率
	CNN	Pd		CNN	Pd		CNN	Pd
4	0.234 0	0.329 1		0.224 4	0.419 8		92.46%	66.09%
5	0.218 0	0.305 9		0.222 8	0.396 7		93.32%	68.83%
6	0.197 8	0.282 9		0.199 3	0.376 6		94.91%	71.17%
7	0.185 2	0.267 5		0.194 7	0.361 7		95.41%	73.02%
8	0.173 9	0.259 1		0.189 7	0.351 5		95.78%	74.43%
9	0.164 6	0.254 4		0.191 3	0.343 1		96.08%	75.52%

下载: 导出CSV 

| 显示表格

从表2可以看出，当时窗长度从4 s增加至9 s时，模型估算震级的准确率从92.46%提高至96.08%，增幅为3.62%，估算误差标准差从0.234 0降低至0.164 6，均值从0.224 4降低至0.191 3，这表明本文所提出的模型存在随着地震记录数据的增加震级估算结果愈加接近真实震级的趋势。同样，P_d参数方法也存在相同趋势，但其震级估算的精度明显低于本文提出的CNN模型，例如，当t＝9 s时P_d震级估计准确率仅为75.52%，甚至小于3 s时窗下本文模型估计的准确率89.92%。

从图7可以观察到，随时窗长度增加，M_JMA＞6.5地震事件的模型预测误差散点逐步接近1 ∶ 1回归线，即估算误差逐步减小，说明震级低估现象有一定程度的改善。特别是针对2011年3月11日的M_JMA9地震事件，随时窗从3 s增大到9 s时，模型对24条记录的震级预估结果的误差均值从最初3.417减小至2.271，虽然预测误差减小了很多，但是模型对巨大地震的预测精度仍未达到要求。

4.   讨论与结论

本研究提出的基于深度学习估算震级的方法在震级估算的准确性方面显著优于传统的基于特征参数的P_d参数方法。究其原因，其一是因为卷积全连接网络处理地震频域数据的过程完全由模型内部算法运行（Hinton，Salakhutdinov，2006），有效规避了人为定义的特征参数和震级估算经验公式，直接由模型内置卷积池化层从地震记录数据中自动提取特征图后通过多层全连接层能更好地拟合上层特征图、辅助输入参数与输出标签震级之间的非线性关系；其二是在使用数据方面，本文模型较P_d参数额外考虑了v_S30场地信息，一定程度上也有助于模型拟合不同场地条件下的地震记录特征，使其在震级估算中的表现更优。

虽然模型在各震级区间的误差均小于P_d参数方法，但当震级超过M_JMA6.5时，模型估算震级明显低于实际的震级，其原因有三：一是大地震的断层破裂时间长，需要更长时段的波形数据；二是大地震的样本较少，没有足够的训练和检测数据；三是大震中的断层破裂情况极其复杂，需引入更多的地震特征参数进行信息补充。例如：对于2011年3月11日的东日本M_JMA9大地震而言，其地震破裂是在多个断层多次发生，Asano和Iwata （2012）认为是由四次破裂叠加而成，虽然模型可随地震持续时间的增长来提取更多有效信息进行内部超参数的权重值更新，以此获取更高精度的震级估算结果，但是对于大地震的预测结果仍有较大的提升空间，这方面的改进尚需更多大地震数据的积累，同时需要融合地震破裂的一些物理规律，或是引入额外的诸如GPS数据等。

时效性方面，同样是采用深度学习方法，相较于Lomax等（2019）以及Mousavi和Beroza （2020）等需要几十秒的地震波形数据，而本研究的模型只需前几秒（3—9 s）的数据，所以本文模型的时效性好，更能满足地震预警快速及时的要求。

数据方面，按年份制作的KiK-net和K-NET训练测试数据集虽然在时序顺序上保证了事件的独立性和合理性，但是本研究只选取KiK-net和K-NET的井上地表记录。相比于井下台站记录来说，井上的数据保有量多，但数据噪声大，质量相对较差。后期计划加入井下的数据进行相关对比分析。另外，本研究仅基于单台的垂直分量加速度的频域信息开展相关研究，未来将加入其余两分量信息并同时考虑多台站联合分析，期望能进一步提高本模型的震级预测精度。

本文采用深度学习方法建立了一种基于地震动频域信息的多层全连接的卷积神经网络，通过由KiK-net和K-NET台网16万4 547条地震动记录制作的数据集开展训练和测试，并与传统P_d预警参数法的计算结果进行了对比，得出如下结论：

1） 将CNN深度学习回归模型应用于地震频域数据集，并辅助以震源信息和场地v_S30信息，能取得较好的震级估算结果。

2） 本文基于单台站竖向地震动频域信息所提出的深度学习方法，不依赖人为定义的特征参数和经验公式，能够实现端到端的震级估算，较基于特征参数P_d的传统地震震级估算方法具有更高的精度。

3） 本文所提出的模型只需要3 s的地震动数据就可以达到89.92%的估算精度，并且随着P波数据的增加，能稳步提升准确率。

4） 对于M_JMA＞6.5地震事件，模型的震级估算结果虽然优于P_d方法的结果，精度也随着时窗尺度的增大而不断提高，但是明显存在着低估的趋势，这是未来需要研究解决的重点。

日本防灾科学技术研究所提供强震数据下载服务，在此表示感谢。