2. 美国哥伦比亚大学应用物理和应用数学系, 纽约 10027;
3. 中国地质大学大数据学院, 武汉 430074
2. University of Minnesota, New York 10027, USA;
3. China University of Geoscience, Wuhan 430074, China
随着地震仪在全球范围的大规模部署, 人类已经进入地震大数据时代.如何处理地震网络每天收集的大量数据, 特别是从质量参差不齐的连续波形记录中分离出地震和噪声是一个非常棘手的问题.例如:大地震发生之后, 抗震救灾与应急救援需要对余震序列进行动态的监测与快速定位, 地震编目, 地球内部速度结构研究如层析成像等都是以大量地震事件的挑选与到时拾取作为基础, 而目前这些工作仍主要靠低效的人工完成(Fang et al., 2015).
为了解决这一问题, 多年来人们发展了各种自动地震识别算法, 比较典型的有:基于振幅的长短窗、基于Akaike信息准则(Akaike information criteria, 简称AIC)的方法(Akaike, 1974), 以及基于波形相似性的波形自相关、互相关、指纹和相似性阈值法.各种方法都有自己的优势和局限性:长短窗法用短时平均(Short-Term Average, STA)和长时平均(Long-Term Average, LTA)之比来反映信号振幅、频率等的变化, 当地震信号到达时, STA/LTA值会有一个突变, 当其比值大于某一个设定阈值时, 则判定有地震事件发生, 这种方法计算快速, 对于脉冲、信噪比高的信号识别准确率较高, 但会漏掉低信噪比的信号(Allen, 1978; Withers et al., 1998; Trnkoczy, 2012); AIC方法, 依据AIC寻求波形曲线的极小点作为震相的到时点, 包括自回归AIC方法(Sleeman and Van Eck, 1999; Akazawa, 2004)、基于峰度的Kurtosis-AIC方法(Küperkoch et al., 2010)以及基于小波变换的wavelet-AIC方法(Zhang et al., 2003)等, 该类方法对震相和频率的变化十分敏感, 常被用于震相到时的精确拾取, 然而, 与STA/LTA类似, AIC的结果也强烈依赖于信噪比以及检测区间; 波形相关性(或称模板匹配)分析可以用来检测源自单个区域具有相同的震源机制的地震, 对于重复发生的地震检测十分有效, 然而, 其检测精度很大程度上取决于使用的模板数量, 此外随着模板数据库的增长, 计算量也呈指数级增加(Skoumal et al., 2014); 指纹和相似性阈值法(简称FAST)使用了结合计算机视觉技术和大数据处理方法的Waveprint音频指纹算法, 对波形进行特征提取, 将每个波形浓缩成一个只保留其主要特征的紧凑“指纹”, 使用局部敏感哈希算法降低相似性搜索量, 该方法在检测灵敏度、适用性、计算效率方面都有不错的表现, 但内存和计算开销较大(Yoon et al., 2015).
神经网络算法是指旨在模拟人类大脑学习过程的一类算法, 它利用彼此互联的非线性“神经元”组成的复杂网状结构, 可以模拟输入特征集合(如地震波形)与预期输出值(如震相类型, 到时等)之间复杂的非线性关系, 并能进一步对新的输入特征给出有效预测(Goodfellow et al., 2016).类似于模板匹配, 神经网络算法需要事先从大量样本学习不同地震波形的特征, 但其独特优势在于训练好的模型具有泛化能力, 可以发现训练集不包含的新特征.早在20世纪90年代, 就有人将人工神经网络引入到地震波形自动识别中来(Wang and Teng, 1995; Zhao and Takano, 1999; Gentili and Michelini, 2006), 由于早期的神经网络均为全连接层, 受限于当时的计算水平, 很难构建深层网络, 其发展受到很大限制.2010年以后, 随着大数据时代的到来, 以及计算机性能的进一步提高, 神经网络算法的发展迎来新的契机.其中以卷积神经网络的广泛应用最具代表性:卷积神经网络通过共享权重, 比人工神经网络参数更少, 可构建的网络层数更深, 提取特征的能力大大提高, 在语音识别、视觉识别、自然语言处理等多个领域, 都取得了前所未有的突破(Krizhevsky et al., 2012; LeCun et al., 2015; He et al., 2016).与这些问题类似, 地震事件的识别本质上也是一个特征提取与参数优化的问题:首先经过大量预设样本的训练, 不断更新优化神经元的权重和偏重参数, 使得损失函数(用于衡量模型输出值与真实值之间的差异)最小, 最终得到保存了不同地震与噪声样本特征的模型, 应用这一最优模型就可以对全新的数据进行特征识别.近几年来, 卷积神经网络已经在地震事件检测方面得到了一系列应用:孔庆开等(Kong et al., 2016)利用人工神经网络区分智能手机网络收集的加速度波形中的地震信号与人类活动噪声; Perol等(2018)应用卷积神经网络对俄克拉荷马地区的连续记录波形进行了识别与定位; Ross等(2018)基于南加州地区的273882多个地震事件建立数据集并训练CNN模型, 其得到的模型具有较强的泛化能力, 即使对与训练数据所属区域构造不同的地区, 如日本Kumamoto地震的事件波形, 也能有效检测.
本研究以汶川地震之后四川及邻区14个台站7—8月期间手动挑选的13839条和8900条地震事件波形分别构建训练数据集和测试集, 搭建深度卷积神经网络进行模型训练和测试.训练好的CNN模型不仅在连续波形识别上精度和召回率优于传统方法, 而且能够有效识别大量人工挑选遗漏掉的小震事件, 其检测灵敏度、适用性、计算效率等综合性能相比传统方法有了很大的提升.
1 CNN方法及性能对比 1.1 数据收集与数据预处理2008年5月12日汶川8.0级地震发生以后, 余震非常频繁, 其震级分布面广(0~6.6级不等), 以0.5~3级的小震居多, 同时前人对该区进行了大量的研究, 积累了充分的余震目录和到时拾取资料, 这些资料具有很高的可信性, 为本研究利用CNN网络地震波形识别提供了必要的训练数据和检验手段.
首先, 深度学习需要收集大量样本数据.数据样本太少会导致模型识别能力的缺陷, 同时易造成过拟合.理想的数据集的每一类别都需要充分的样本数量, 并且能够覆盖地震波形所具有的各种复杂特征.我们的数据集来源于四川及邻区14个中国地震台网固定台站7—8月间的连续波形记录, 分事件和噪声两类.在事件类方面, 从8月连续波形中截取13839个事件波形片段, 作为训练数据, 从7月连续波形记录中截取8900个事件波形片段作为测试数据, 所有事件波形均取P到时前5 s至S到时后10 s; 在噪声类方面, 采用排除掉事件波形之后随机筛选截取.数据集所有样本均统一为30 s数据(长度不足的以零值填充, 长度超过的从尾部截断).此外为尽可能减少事件人工标注和噪声随机筛选带来的误差, 我们手动对数据集进行了清洗, 更正较为明显的标注错误.
为了增加数据量, 我们还对清洗后的数据集进行了增强操作.由于地震三分量波形可以看作一维三通道图像数据, 我们参照图像识别, 采用了如下几种常用的数据增强方式:平移、加噪(模糊处理)、滤波.平移是指我们以P、S到时为中心, 在[-3 s, +3 s]时间窗内沿时间轴滑窗, 对震相位置进行微调处理; 加噪是指对波形的每个采样点进行随机扰动, 本研究采用高斯噪声扰动; 滤波是指我们对波形使用[0, 20 Hz]的带通滤波.通过上述操作, 我们将事件样本数量扩展到65000个, 为下一步训练模型打好了基础.
图 2展示了数据预处理前后不同的训练和损失函数曲线, 可以看出用处理后的数据集算法能更快收敛, 而且我们将预处理之后的数据集时间戳列表放在https://github.com/mingzhaochina/dataset, 供大家测试使用.大家可先向中国地震台网中心提交连续波形数据申请, 并根据时间戳自行截取波形建立数据集.
在用于地震-噪声分类的CNN网络结构上, 参考了在手写数字分类问题上取得成功的经典卷积神经网络“Le-Net5”(LeCun et al., 2015).如图 1所示, 其中卷积层的功能在于提取特征, 池化层通过降采样过滤掉一些特征.与手写数字图像分类的二维卷积不同, 地震波形采用一维卷积, 其卷积核为3, 通道数为32.训练数据按批次输入, 3000×3张量, 对应于30s长三分量波形, 经过逐层特征提取与过滤, 到第八层时减为12×32, 并一维化为384个特征, 经过全连接层之后使用ReLU激活函数计算每一类的概率, 并与预先设置的分类阈值比对, 低于阈值则对应“0”标签(噪声), 高于阈值则为对应于“1”标签(地震).
在所用的算法方面, 我们使用了L2正则化以及随机梯度下降算法来最小化交叉熵损失函数, 并用ADAM优化算法以及可变时间步长, 进一步提高了计算效率.在深度学习训练中, 将所有数据迭代完一遍称为一个epoch, 本研究使用Tensorflow程序(https://www.tensorflow.org), 在一块GTX 1080ti GPU卡训练10个epoch(约40000步迭代)只需要20 min, 其训练达到的最高精度为98.2%, 在验证集上对8900个样本的识别仅需不到20 s, 其验证精度最高达到94.4%(图 2).CNN方法的计算效率介于传统方法(STA/LTA, AR-AIC等)与模板匹配之间, 前者仅需数分钟且无需GPU计算, 后者往往需要几百甚至数千GPU计算机时, 可参见Perol等(2018).
2 连续波形检测 2.1 连续波形检测对比实验我们将在训练集和测试集上精度最高的CNN模型用于连续波形检测, 并与经典STA/LTA算法(Krischer等, 2015)和fbpicker算法(Chen and Holland, 2016)的识别效果对比.所用数据为四川台网MXI台一个月连续波形(2008.8.1 —9.1), 并以该时间段专家手动挑选的1901个事件波形作为参照.我们用精度P和召回率R来衡量不同算法的识别效果(下标e和n分别代表事件和噪声), 其定义为
其中, Tp为真正例, 即算法识别的事件是真实事件, 反之则为假正例Fp, Tn为真反例, 即算法识别的噪声是真实噪音, 反之则为假反例Fn.精度高说明误检率低, 召回率高则证明算法漏检率低, 只有二者同时都高算法或模型才具有实用价值.表 1列出了三种方法各自的精度和召回率.其中CNN网络识别出3763个事件波形, 其中有1803个与参照样本相符.传统算法在精度和召回率上均低于CNN方法.
需要指出的是, 作为参照的专家样本大多经过多台震相关联、定位等后续处理, 其准确度很高, 但由于人眼识别的主观性和局限性, 完整性未必很高(例如很容易漏掉孤立事件和信噪比低的微震事件), 而已有研究表明, CNN模型的一大优势在于能够有效识别微震.因此严格意义上来讲, 这里的精度Pe并不能完全反映CNN算法的分类效果.如图 3所示, 以2008年8月25日MXI台一日的连续波形识别为例, 这一日共有43个专家挑选事件, 但进一步的甑别我们发现该日可能还存在约230个震幅较弱的微震事件, 应用CNN方法可以更好地识别这些微震(图 3b), 如果将这些微震考虑进去, 则经过计算, 算法的精度和召回率大约在73%和62%左右, 同时CNN方法还可以通过将误分类波形代入数据集重新训练, 从而获得越来越好的实际分类效果.
以上我们通过模型训练、验证、连续波形检测三个阶段, 证明了CNN算法和模型具备对复杂地震波形的强大分类能力.为了进一步检验CNN模型对不同地域的连续地震波形的泛化能力, 我们将该模型尝试用于全国台网连续波形的记录.我们一共采用了441个台站的连续波形(仪器工作正常, 数据质量较好), 其分布遍布全国, 时间为2017年12月11日—19日期间.此外, 由于不同的区域具有不同地质构造背景, 其事件波形也各不相同, 因此我们参考历史地震活动性区域划分(图 4a), 按照地震震中之间的欧几里得距离进行了K-Means聚类(图 4b), 将441个台站的事件波形划分为9类(对应9个区域), 检验CNN算法性能与地域的相关性.
我们的识别方法如下:
(1) 应用长短窗方法, 根据各地台站的信噪比条件选择合适的阈值(在大部分情况下阈值>5)做初步筛选, 这样做是为了提高算法效率, 过滤掉大部分环境噪声;
(2) 将筛选出的波形片段用预训练好的CNN网络识别;
(3) 对CNN网络识别为地震的波形片段, 我们进一步采用dbshear震相自动识别程序进行挑取, 该方法的流程是:首先进行P波震相的挑取, 然后使用偏振能量分析将地震波水平分量上的P波能量去除, 再使用STA/LTA进行S波初拾取, 最后应用kurtosis函数法进行S震相的准确拾取(参见Ross et al., 2016);
(4) 我们将拾取到的P、S到时与中国地震台网发布的参考地震目录进行关联:首先利用一维速度模型计算理论参考到时, 最后将理论到时与拾取P、S到时进行比较, < 2 s的被认为是正确识别.
应用以上方法, CNN网络检出7016段波形, 用方法检测到1380对P, S到时, 并与540个地震目录事件成功关联.中国地震台网目录在这一时间段共记录1级以上地震事件1073个, 二级以上246个, 经过比较, 我们的方法对1级以上事件总体识别准确率为54%(图 4b), 二级以上为80%(图 4c).
2.3 地域分析从地域统计来看, ④⑤(大致位于川滇地区)、①③(大致位于新疆地区)、⑦(华北地区)与参考目录匹配较好, 其余地区较差(表 2).我们推测这可能与我们的训练数据大量采用了汶川地震的余震波形有关.以上地区均为我国地震多发区域, 除了大震强震, 3级以下的小震微震活动也非常频繁, 我们的模型对于识别这些地区的微震小震特别有效.
地震学是一门数据驱动的科学, 从海量原始波形数据中分离出地震信号是地震目录产出的第一步, 因此也是地震学研究的基础.与基于单一或多个特征函数的传统方法相比, CNN网络的精度和召回率(误检和漏检率)相比传统算法有非常明显的提升(降低), 并且训练好的模型具有非常稳定的输出能力, 不需要根据不同信噪比水平的数据频繁调整阈值.与模板匹配、FAST等严格基于波形相似性的方法不同, CNN神经网络从训练数据中提取的是抽象特征, 这意味着模型具有更强的泛化能力, 能够与更多样的波形特征相匹配, 所以CNN往往能够检测到新的地震类型, 这也是深度学习方法的主要优势所在.在2.3节中我们基于汶川余震数据训练的CNN网络, 对地质构造背景不同①③⑦的也具有较强的泛化能力, 便是这种能力的证明.
尽管CNN展示了其在地震波形识别领域良好的应用前景, 目前来看, 仍具有以下局限性:首先, 作为一种监督学习方法, 神经网络对于训练样本的依赖性较强, 如果缺乏大量预设样本, 或是预设样本过于单一, 无法涵盖所分类对象所具有的一般特征, 则很容易导致模型过拟合或根本无法训练; 其次, 大量正负样本的标注十分棘手, 并且很难用统一标准评判, 因此, 即使训练好的模型能在测试集上达到很高的精度, 其最终质量好坏往往只能通过实际应用检验, 在实际应用中我们经常遇到这样的情况:有些CNN误识别的波形, 使用长短窗或自回归方法甚至肉眼都可以轻松识别, 因此在实际应用中应该设置某种投票机制, 将基于不同原理的多种方法结合起来, 才能取得最佳的识别效果; 最后, 目前CNN等深度学习算法背后的原理尚不清楚, 尽管已经在地震学领域得到一系列应用, 但不同学者提出的不同神经网络结构均是通过反复试错之后取得, 其依据是模型在训练和测试集上的表现而不是基于对算法本身的理解, 因此很多模型放到实际应用中都会出现“性能瓶颈”, 即精度和召回率在到达一定程度后很难继续提升, 而我们很难分清这到底是由于训练数据样本分布的不足, 还是模型设计本身存在缺陷.
4 结论和展望本研究通过搭建CNN深度神经网络模型, 建立数据集进行模型训练、验证与实际波形检测, 证明了CNN网络的强大复杂特征分类能力可以推广应用于地震波形分类这一应用场景.同时, 在训练样本数量不够大的情况下, 我们对数据集采用了平移、加噪和滤波等增强操作来扩大训练样本, 对防止模型过拟合, 增强模型容错能力起到了一定作用.对不同地质构造背景区域的波形识别, 则证明了我们基于余震数据训练的CNN网络, 对其他构造区域也具有泛化能力.
根据深度学习在其他领域的应用经验, 其优势主要在于处理多分类问题.由于地震的复杂性, 其波形可能呈现多种形态, 很难用地震-噪音这样的简单二分类问题囊括.因此, 如何建立一个包含不同地震类型(例如近震、远震; 天然地震、爆破、塌陷、滑坡、核爆等)、每一类具有充足样本数量的大型数据集是一个非常繁重的任务.在下一步工作中, 我们将增加来自不同地域、不同构造类型的地震和噪声样本, 根据样本的特征不同做更进一步分类, 同时应用更深层的网络进行训练, 获得不但能够识别地震信号, 还能推断地震类型的神经网络模型, 使得地震数据的检测与识别更加智能化.
致谢 感谢地震局球所房立华研究员协助提供人工挑选数据, 感谢天池余震识别AI大赛平台提供连续波形数据.感谢中国地震局第二监测中心提供部分台站实时数据, 以及长安大学孙少波老师、球所张贝博士为实时数据流下载提供帮助.感谢加州理工的Zachary Ross博士提供的dbshear程序作为对比.Fbpicker程序来自https://github.com/austinholland/PhasePApy.
Akaike H. 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6): 716-723. DOI:10.1109/TAC.1974.1100705 |
Akazawa T.2004.A technique for automatic detection of onset time of P-and S-Phases in strong motion records.//13th World Conference on Earthquake Engineering.Vancouver B C, Canada: International Association for Earthquake Engineering.
|
Allen R V. 1978. Automatic earthquake recognition and timing from single traces. Bulletin of the Seismological Society of America, 68(5): 1521-1532. |
Chen C, Holland A A. 2016. PhasePApy:A robust pure python package for automatic identification of seismic phases. Seismological Research Letters, 87(6): 1384-1396. DOI:10.1785/0220160019 |
Fang L H, Wu J P, Wang W L, et al. 2015. Aftershock observation and analysis of the 2013 MS7.0 Lushan earthquake. Seismological Research Letters, 86(4): 1135-1142. DOI:10.1785/0220140186 |
Gentili S, Michelini A. 2006. Automatic picking of P and S phases using a neural tree. Journal of Seismology, 10(1): 39-63. DOI:10.1007/s10950-006-2296-6 |
Goodfellow I, Bengio Y, Courville A, et al. 2016. Deep Learning. The MIT Press.
|
He K M, Zhang X Y, Ren S Q, et al.2016.Deep residual learning for image recognition.//2016 IEEE Conference on Computer Vision and Pattern Recognition.Las Vegas, NV, USA: IEEE.
|
Küperkoch L, Meier T, Lee J, et al. 2010. Automated determination of P-phase arrival times at regional and local distances using higher order statistics. Geophysical Journal International, 181(2): 1159-1170. |
Kong Q K, Allen R M, Schreier L, et al. 2016. Myshake:A smartphone seismic network for earthquake early warning and beyond. Science Advances, 2(2): e1501055. DOI:10.1126/sciadv.1501055 |
Krischer L, Megies T, Barsch R, et al. 2015. ObsPy:A bridge for seismology into the scientific Python ecosystem. Computational Science & Discovery, 8(1): 014003. DOI:10.1088/1749-4699/8/1/014003 |
Krizhevsky A, Sutskever I, Hinton G E.2012.Imagenet classification with deep convolutional neural networks.//Advances in Neural Information Processing Systems, 1097-1105.
|
LeCun Y, Bengio Y, Hinton G. 2015. Deep learning. Nature, 521(7553): 436-444. DOI:10.1038/nature14539 |
Perol T, Gharbi M, Denolle M. 2018. Convolutional neural network for earthquake detection and location. Science Advances, 4(2): e1700578. DOI:10.1126/sciadv.1700578 |
Ross Z E, White M C, Vernon F L, et al. 2016. An improved algorithm for real-time S-wave picking with application to the (augmented) ANZA network in southern California. Bulletin of the Seismological Society of America, 106(5): 2013-2022. DOI:10.1785/0120150230 |
Ross Z E, Meier M A, Hauksson E, et al. 2018. Generalized seismic phase detection with deep learning. Bulletin of the Seismological Society of America, 108(5A): 2894-2901. DOI:10.1785/0120180080 |
Skoumal R J, Brudzinski M R, Currie B S, et al. 2014. Optimizing multi-station earthquake template matching through re-examination of the Youngstown, Ohio, sequence. Earth and Planetary Science Letters, 405: 274-280. DOI:10.1016/j.epsl.2014.08.033 |
Sleeman R, Van Eck T. 1998. Robust automatic P-phase picking:An on-line implementation in the analysis of broadband seismogram recordings. Physics of the Earth and Planetary Interiors, 113(1-4): 265-275. |
Trnkoczy A.2012.Understanding and parameter setting of STA/LTA trigger algorithm.//Bormann P ed.New Manual of Seismological Observatory Practice 2 (NMSOP-2).Potsdam: Deutsches GeoForschungsZentrum GFZ.
|
Wang J, Teng T L. 1995. Artificial neural network-based seismic detector. Bulletin of the Seismological Society of America, 85(1): 308-319. |
Withers M, Aster R, Young C, et al. 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection. Bulletin of the Seismological Society of America, 88(1): 95-106. |
Yoon C E, O'Reilly O, Bergen K J, et al. 2015. Earthquake detection through computationally efficient similarity search. Science Advances, 1(11): e1501057. DOI:10.1126/sciadv.1501057 |
Zhang H, Thurber C, Rowe C. 2003. Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings. Bulletin of the Seismological Society of America, 93(5): 1904-1912. DOI:10.1785/0120020241 |
Zhao Y, Takano K. 1999. An artificial neural network approach for broadband seismic phase picking. Bulletin of the Seismological Society of America, 89: 670-680. |