2. 美国哥伦比亚大学应用物理和应用数学系, 纽约 10027;
3. 中国地质大学大数据学院, 武汉 430074
2. Department of Applied Physics and Applied Mathematics, Columbia University, New York 10027, USA;
3. Department of Big Data, China University of Geoscience, Wuhan 430074, China
基于地震台网的观测数据,编制地震观测报告(简称“地震编目”),是中国地震局地震监测工作的核心任务之一.这些数据是地震速报、地震预警等防震减灾工作的重要环节,也是地震层析成像、震源机制反演和地球内部结构等研究的重要基础.地震编目的主要内容包括震相识别、到时拾取、振幅和周期测量、初动极性标注和地震定位等工作,其中最重要的环节是震相识别和到时拾取.在区域台网中,最常见的震相是Pg和Sg.Pg、Sg分别指来自上地壳内震源的上行P、S波,或射线底部到达上地壳的P、S波,其震中距一般不超过10°(中国地震局,2017).
传统的震相识别和到时拾取主要有人工处理和自动处理两种.人工识别主要依靠地震台网分析人员通过浏览连续波形,识别地震.这种方法效率低、主观性强、一致性差,尤其是遇到余震序列时,经常会漏掉地震.自动处理算法是根据波形在振幅、频率等方面的不同特征建立特征方程,进行检测.常用方法有:利用短时平均和长时平均比值的STA/LTA方法(Allen, 1978; Baer and Kradolfer, 1987; Withers et al., 1998),利用赤池信息准则寻找全局最小值作为震相到时的AIC方法(Akaike, 1974; Sleeman and Van Eck, 1999; Akazawa, 2004; 刘劲松等, 2013),以及利用偏度(Skewness)和峰度(Kurtosis)等反映数据分布形态的高阶统计量的偏度/峰度方法(Saragiotis et al., 2002; Küperkoch et al., 2010, 2012),等等.这些方法前人做过较多综述(刘翰林和吴庆举, 2017; 李力等, 2018).由于地震波形受仪器类型、环境噪声以及震源复杂性的影响,这些自动处理方法局限性也很明显,主要表现在如下方面:稳定性较差,不同的参数和阈值选择对结果影响很大;灵敏度较低,对信噪比较低的微震信号,检测效果往往难以达到预期;经常混淆P、S震相;S震相到时拾取精度较低.
近年来,由于机器学习和深度学习技术的快速发展,国内外地震学家越来越多利用机器学习和深度学习对地震和噪声进行分类,如随机森林(Maggi et al., 2017)、支持向量机(Kortström et al., 2016; 蒋一然和宁杰远, 2019)、卷积神经网络(CNN)(Perol et al., 2018)等方法,参见Kong等(2018)的综述.其中,CNN方法近两年来得到了较多应用,其强大之处在于多层结构的设计能自动学习不同层次的特征:较浅的卷积层能学习到一些局部特征,而较深的卷积层能够学习到更加抽象的特征,从而有助于同时提升模型分类能力和泛化能力.CNN在地震与噪声分类、震相分类等问题上得到了较好应用,即使对一些低信噪比的数据,甚至肉眼难以识别的微震信号,也能有效甄别(如Ross et al., 2018;赵明等,2019).然而,在震相到时拾取方面,机器学习和深度学习方法的拾取精度有待进一步提高.
2014年,Long等提出全卷积网络(FCN),使得CNN无需全连接层即可进行密集的像素预测,CNN从而得到进一步普及(Long et al., 2015).Ronneberger等在FCN基础上提出U网络并应用在医学图像分割上,取得了媲美专业医生的精准度(Ronneberger et al., 2015).U网络的主要特点是不仅实现了输出层与输入层在分辨率上的一致性,还通过对称式的结构设计融合了FCN网络中的低维和高维特征.图像分割是对图像的每个像素进行分类,根据像素点对应的类别概率进行定位和边缘检测,与此十分类似,震相分类则是确定P、S波形的位置,到时拾取则是寻找P、S对应的起震点.我们可以利用标注好到时及误差范围的地震P、S波形数据作为样本训练U网络模型,使其输出波形每个采样点所对应的P、S起震点概率,概率最大的采样点即为P、S到时.基于以上基本原理,本研究利用汶川地震余震和首都圈地震台网记录的Pg、Sg到时数据建立数据集,应用U网络进行Pg、Sg震相分类与到时拾取,并与STA/LTA、Kurtosis方法比较了震相识别正确率、拾取误差、漏检率,以及在不同信噪比水平下的表现进行了对比分析.
1 方法与数据 1.1 U形神经网络本研究在Perol等(2018)和Zhu and Beroza(2019)的研究基础上,编写了基于U网络的震相识别程序(Unet_cea,https://github.com/mingzhaochina/unet_cea),相比于经典U网络,主要区别如下:
(1) 网络降维:U网络默认为二维或三维图像数据所设计,而地震波形数据为一维时间序列,因此所有卷积层、池化层、drop out层,以及输入输出张量均需降为一维.
(2) 损失函数:深度学习算法进行分类的数学原理实际是对目标函数进行最优化,通过反复比较当前网络的预测值与我们期望的目标值的差异,不断调整每一层的权重和偏差等超参数,使得差值最小化.衡量预测值与期望值差异的方程被称为损失函数(或目标函数).本研究采用的损失函数为
其中Y′为采用二值化编码的标签,i=1, 2, 3分别表示噪音、Pg、Sg,n为波形采样点数,j=1, …, n为采样点序号,其表达式为
Y为最后一层softmax函数计算得到概率值,其表达式如下:
其中z为最后一层的输出张量([m, n, 3]),m为输入数据个数.
如图 1所示,U网络的基本组件为增采样层和降采样层.降采样层由两个一维卷积层、一个池化层组成,中间会随机加入一些drop out层来防止过拟合;增采样层由一个转置卷积层、一个裁切层、一个卷积层组成,同样根据情况适当加入drop out层.网络输入为三分量地震波形.根据前人在图像分割上的应用经验,左半部分4个降采样层执行卷积和池化操作,提取震相的抽象特征,有助于解决震相定位的问题,右半部分4个增采样层执行转置卷积和左右对称层之间互连通等操作,可以逐步恢复震相的细节特征,从而有助于解决震相分类的问题.最后通过激活函数计算Pg、Sg或噪音概率值,并与预设阈值比对确定该采样点的类别.
与CNN网络一样,U网络也需要大量预设样本进行训练.本研究使用的样本数据来源于两部分:一部分来自于2017年中国地震局地球物理研究所与阿里云联合举办的“余震捕捉—AI大赛”(Fang et al., 2017), 包括2008年6—8月间汶川地震后震源区附近15个宽频带地震台站记录的数据,数据信噪比水平集中于2~5之间(图 2a);另一部分数据由北京数字遥测地震台网提供,包含了首都圈地区79个地震台站2014年1月至2015年8月的连续波形数据和震相观测报告,样本事件的震级主要分布在M0.0~2.0,信噪比集中分布在1.6~5区间(图 2b).
如图 3的数据标注样本所示,我们取Pg到时前3 s至Sg到时后10 s的波形,并统一裁切为30 s长度(长度不足的以零值填充),取30 s的原因是由于上述两个区域地震事件的震中距小于150 km,Pg和Sg到时差最大不会超过25 s,可以保证波形同时包含Pg和Sg震相.然后对所有数据进行了去均值、归一化等预处理.为了保证数据集质量,我们手动对数据集进行了清洗,更正较为明显的标注错误.最终,对汶川余震我们截取和标注的波形为41722段,首都圈为47622段.在训练集与验证集分配比例上,我们使用80%的数据做训练,剩余20%作为测试.汶川余震数据质量和到时精度很高,震源较集中,震级相对高,适用于检测U网络模型的精度;而首都圈地区台站分布广,震源较分散,微震事件多,适合用于检测模型泛化能力.
此外,我们采用了与赵明等(2019)相同的数据增强操作,包括平移、加噪(模糊处理)、滤波,并将预处理之后的数据集时间戳列表放在https://github.com/mingzhaochina/dataset,方便科研人员下载使用与测试.
1.3 U网络模型训练与测试我们使用Keras(https://keras.io)构建U网络,并在tensorflow(https://www.tensorflow.org)上进行训练和测试.在算法选择上使用了L2正则化以及随机梯度下降(SGD)算法来计算交叉熵损失函数,以及ADAM优化算法、阶梯型学习率进一步优化计算过程.U网络的训练在中国科学院网络中心(http://www.cnic.cn)进行,使用1个计算节点(12核CPU,64GB内存,1块Telsa P100 GPU卡)上运行8个周期(epoch,表示将所有样本遍历一遍)只需15 min.训练好的神经网络模型可以实时用于测试集,对9524个测试样本的验证仅需20 s.
训练及测试结果如图 4所示.其中在训练集上训练完第一个周期(约5000时间步)后精度即达到98%,之后的7个周期均稳定在这一水平(图 4a).在测试集上Pg、Sg震相的识别正确率均随着训练时间的增加稳步上升:Pg由78%(20000步)至86%(40000步),Sg由74%(20000步)到83.9%(40000步),这说明模型并没有发生过拟合,并且训练时间越长,精度越高.
我们将在训练集和测试集上表现最优的U网络模型,应用于首都圈台网2015年9—12月间事件连续波形的震相挑取.首先根据河北省地震局提供的这一期间的人工编目报告(包括发震时刻,各台站记录的Pg、Sg到时等),以地震事件的发震时间为开始,以最远台站接收到Sg震相后30 s所对应的时间为止,对各个台站统一截取三分量连续事件波形数据.然后每隔2 s依次取30 s长波形,作为U网络的输入.由于输入的截断连续波形很大一部分是噪声,或者只包含一部分震相,因此其识别难度更大.为了客观评估U网络模型在连续波形应用中的性能,我们选取了两种常用的震相拾取算法进行对比:STA/LTA算法(Krischer et al., 2015)和Kurtosis算法(Ross et al., 2016).对于三种方法的检测结果,参照人工编目到时作为标准答案评价.为保证结果可靠性,我们仅选取人工编目结果中≥3台以上同时监测到的事件,最终确定了2538段标注波形.计算各种方法的均方根误差(RMSE)、命中率(H)和查全率(R),三个指标的计算公式如下:
(1) RMSE为算法挑取到时Tij与参考到时tij之间均方根误差,用于反映算法到时挑取的精度,n为台站数,m是每个台的震相数,i, j为序号;
(2) H为正确率,表示算法拾取的震相在误差范围内且与参考震相类型一致.Te为自动算法到时检测正确的次数,Fe为检测错误的次数.当自动拾取的到时与参考答案的差异小于一定阈值(Pg为±0.25 s,Sg为±0.5 s),且震相类型一致时,Te增加一次,否则Fe增加一次;
(3) R为查全率,用于检查不同算法拾取的到时与参考到时相比的完整程度,其中Fn为漏检的震相个数.
计算结果如表 1所示.总体而言,在震相识别正确率上,U网络对Pg、Sg震相的得分分别为81%和79.1%,均高于STA/LTA和Kurtosis方法的情况.在反映震相到时拾取精度的RMSE值上,U网络拾取平均误差也比STA/LTA和Kurtosis方法的低.在反映震相识别完整性的查全率上,三种方法都不是很高,只有U网络接近一半,这说明仍有相当数量的微震震相未被检测出来.图 5给出了不同情形下三种方法各自的震相到时拾取结果,以及与人工处理结果的对比.
实时地震学和地震早期预警都建立在对震相的准确快速识别上,区域地震台网的日常编目工作也需要对大量波形进行分类与震相识别,这些工作都迫切需要实现对地震数据的快速、高效和自动化处理.由于地震的复杂多样性,实现震相的自动化识别难度很大.目前的震相自动识别技术很难同时在检测精度、稳定性和计算效率上都达到令人满意的水平.
本研究显示,U网络具有良好的泛化能力.对9524个具有不同信噪比、不同仪器类型、不同震级和震中分布的地震波形测试取得了Pg识别正确率最高86%和Sg最高83.9%的结果.在识别难度更高的首都圈地震台网实际波形检测上,U网络对Pg、Sg震相的正确识别率分别达到了81%和79.1%,与参考到时的均方根误差分别为0.41 s和0.54 s,这些指标均优于STA/LTA、Kurtosis等传统的自动检测方法.
尽管目前仍没有充分的理论阐明深度学习算法的具体工作过程,但在测试数据集和连续波形数据上较高的检测正确率和较小的均方根误差表明,U网络在获取震相特征上,远比传统方法仅仅利用振幅或频率等的一两种特征更加有效.例如对于Pg、Sg到时相近的波形,Sg常常受Pg尾波干扰甚至淹没,传统方法基本很难给出有效的Sg到时,而U网络不但能学习到Pg、Sg的特征,还能区别二者到时的细微差异.在计算量和耗时方面,使用GPU卡加速的U网络算法明显优于模板匹配和波形相关等方法.其主要计算开销来自模型的训练,其余无论是测试集测试还是实际波形检测,计算量和耗时都很低,而且还可以迁移到其他地区进行应用.
目前,将U网络用于实际波形的震相到时检测还存在两个主要问题:
(1) 漏检的事件偏多.虽然U网络在实际波形检测中展示了与测试集一致的精度,但这一定程度上是由于我们的概率阈值较高(Yi>95%),以及比对的人工标注波形质量较高(≥3台同时检测到).如图 6所示首都圈实际波形检测结果,一部分波形未能预测或预测完全错误,另一部分只给出了Pg或Sg预测.我们对总计1100余段未能完全正确识别的波形进行了分析,发现绝大多数(>70%)属于噪声干扰较大的复杂波形,即使是人工检测结果也因人而异.对这部分波形,可以通过震相关联和定位来辅助确定其到时,从而降低漏检率.同时还有一小部分(<30%)属于信噪比较高、震相也比较清晰.我们采用波形互相关方法计算其与训练样本之间的相关系数,发现相关度均较低,这说明这部分波形是我们训练样本中所缺失的特征,应该进一步在类型上扩充训练样本数据,重新加以训练.此外由于在本研究中我们只使用了北京数字遥测地震台网一年的数据,在后续研究中使用5—10年的数据有望显著降低漏检率.
(2) 震相到时挑取仍不够精确.如图 5中不同方法的特征曲线所示,我们训练的U网络目前能够预测的Pg、Sg到时均有一定误差区间,分别为±0.25 s和±0.5 s,而STA/LTA和Kurtosis等传统方法对波形的突变特征比较敏感,可达到±0.2 s甚至更低的拾取误差.当然,传统算法的劣势主要在于不够稳定,当数据信噪比低时,其精度大大降低,而此时U网络仍能保证相对稳定的输出.因此,U网络于传统方法具有互补性,在实际应用中应当结合使用.例如,可以先设定一个较小的阈值,采用STA/LTA检测可能的震相到时,然后使用U网络确定震相类型和到时范围,然后再用AIC、Kurtosis等方法精确拾取震相到时.
4 结论本研究使用四川地震台网和首都圈地震台网共89344个包含Pg、Sg标注震相的事件波形样本,训练和测试了一个可以进行震相识别和到时拾取的U网络.与STA/LTA、Kurtosis等传统震相识别方法的对比结果表明,U网络的均方根误差、命中率和查全率均优于STA/LTA、Kurtosis等传统的自动检测方法,与人工挑选结果的比对表明其准确识别出了81%和79.1%的Pg、Sg震相.
U网络在检测精度和稳定性上均明显优于传统方法,同时可以极小的计算开销应用于首都圈台网的实时地震数据快速Pg、Sg震相检测.当然,U网络目前仍存在漏检和拾取震相不够精确的问题,后续工作包括扩充数据样本量(包括收集漏检样本与使用5—10年的波形数据)重新训练模型、使用传统的震相关联和定位等手段排除震相误检率和漏检率,并结合传统方法进行到时精确拾取.深度学习和人工智能技术未来将在实时地震学、地震早期预警和地震台网自动编目等业务上发挥越来越重要的作用.
致谢 感谢“余震捕捉—AI大赛”提供汶川地震余震的波形和震相到时数据,感谢北京数字遥测地震台网和河北省地震局提供连续波形数据和震相观测报告.感谢加州理工学院Zachary Ross博士提供震相到时自动拾取程序.
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. |
Baer M, Kradolfer U. 1987. An automatic phase picker for local and teleseismic events. Bulletin of the Seismological Society of America, 77(4): 1437-1445. |
Fang L H, Wu Z L, Song K. 2017. Seism Olympics. Seismological Research Letters, 88(6): 1429-1430. DOI:10.1785/0220170134 |
Jiang Y R, Ning J Y. 2019. Automatic detection of seismic body-wave phases and determination of their arrival times based on support vector machine. Chinese Journal of Geophysics (in Chinese), 62(1): 361-373. DOI:10.6038/cjg2019M0442 |
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. |
Küperkoch L, Meier T, Brüstle A, et al. 2012. Automated determination of S-phase arrival times using auto regressive prediction:Application to local and regional distances. Geophysical Journal International, 188(2): 687-702. DOI:10.1111/j.1365-246X.2011.05292.x |
Kong Q K, Trugman D T, Ross Z E, et al. 2018. Machine learning in seismology:turning data into insights. Seismological Research Letters, 90(1): 3-14. DOI:10.1785/0220180259 |
Kortström J, Uski M, Tiira T. 2016. Automatic classification of seismic events within a regional seismograph network. Computers & Geosciences, 87: 22-30. |
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 |
Liu H L, Wu Q J. 2017. Developments of research on earthquake detection and seismic phases picking. Progress in Geophysics (in Chinese), 32(3): 1000-1007. DOI:10.6038/pg20170308 |
Liu J S, Wan Gym Yao Z X. 2013. On micro-seismic first arrival identification:A case study. Chinese Journal of Geophysics (in Chinese), 56(5): 1660-1666. DOI:10.6038/cjg20130523 |
Long J, Shelhamer E, Darrell T. 2015. Fully convolutional networks for semantic segmentation. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. Boston, MA, USA: IEEE, 3431-3440.
|
Maggi A, Ferrazzini V, Hibert C, et al. 2017. Implementation of a multistation approach for automated event classification at Piton de la Fournaise volcano. Seismological Research Letters, 88(3): 878-891. DOI:10.1785/0220160189 |
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 |
Ronneberger O, Fischer P, Brox T. 2015. Invited Talk: U-net Convolutional networks for biomedical image segmentation. International Conference on Medical Image Computing and Computer-Assisted Intervention, Cham: Springer, 234-241.
|
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. 2018. P wave arrival picking and first-motion polarity determination with deep learning. Journal of Geophysical Research:Solid Earth, 123(6): 5120-5129. DOI:10.1029/2017JB015251 |
Saragiotis C D, Hadjileontiadis L J, Panas S M. 2002. PAI-S/K:A robust automatic seismic P phase arrival identification scheme. IEEE Transactions on Geoscience and Remote Sensing, 40(6): 1395-1404. DOI:10.1109/TGRS.2002.800438 |
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. |
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. |
Zhao M, Chen S, Yuen D. 2019. Waveform classification and seismic recognition by convolution neural network. Chinese Journal of Geophysics (in Chinese), 62(1): 374-382. DOI:10.6038/cjg2019M0151 |
Zhu W, Beroza G C. 2019. PhaseNet:a deep-neural-network-based seismic arrival-time picking method. Geophysical Journal International, 216(1): 261-273. |
蒋一然, 宁杰远. 2019. 基于支持向量机的地震体波震相自动识别及到时自动拾取. 地球物理学报, 62(1): 361-373. DOI:10.6038/cjg2019M0442 |
李力, 张建国, 李盛乐. 2018. 余震自动识别技术研究进展. 地震研究, 41(1): 1-13. DOI:10.3969/j.issn.1000-0666.2018.01.001 |
刘翰林, 吴庆举. 2017. 地震自动识别及震相自动拾取方法研究进展. 地球物理学进展, 32(3): 1000-1007. DOI:10.6038/pg20170308 |
刘劲松, 王赟, 姚振兴. 2013. 微地震信号到时自动拾取方法. 地球物理学报, 56(5): 1660-1666. DOI:10.6038/cjg20130523 |
赵明, 陈石, Yuen D. 2019. 基于深度学习卷积神经网络的地震波形自动分类与识别. 地球物理学报, 62(1): 374-382. DOI:10.6038/cjg2019M0151 |
中国地震局. 2017. DB/T66-2016地震编目规范. 北京: 地震出版社.
|