地球物理学报  2021, Vol. 64 Issue (7): 2394-2404   PDF    
卷积神经网络快速挑选接收函数
甘露1,2, 吴庆举1, 黄清华2, 唐荣江1     
1. 中国地震局地球物理研究所, 北京 100081;
2. 北京大学地球与空间科学学院, 北京 100081
摘要:地震三分量波形数据中提取的接收函数受震源复杂性及随机噪声等因素的影响,往往出现一些波形异常现象,需要在资料解释前予以剔除.当接收函数数量较多时,人为挑选质量合格的接收函数将耗费大量的时间.为了高效的挑选高质量的接收函数,本文利用深度学习卷积神经网络(Convolutional Neural Network,简称CNN)方法来对接收函数的质量进行判断.我们使用华北地区和阿巴嘎地区地震观测台站接收的9833个不同的地震事件建立训练集,并用1521个新的地震事件作为检测数据集,得到的训练集的准确率和召回率均达到99%以上,测试集的准确率和召回率分别达到95.3%和92.4%.我们还使用了训练集和测试集以外的数据进行验证,并对比了不同CNN评估结果所对应的波形图,实验证明评估结果与实际的接收函数波形对应良好.此外,对于某些台站的接收函数,可能存在如下问题:波形虽然具有很好的一致性,但由于不符合常规意义下质量好的标准导致CNN无法识别.为解决该问题,本文首先对一定方位角和震中距范围内的接收函数相互求取二范数,再对二范数较低的结果进行统计,并与CNN的挑选结果进行对比,挑选出合格的数据.
关键词: 接收函数      卷积神经网络      快速挑选     
Quick selection of receiver function based on convolutional neural network
GAN Lu1,2, WU QingJu1, HUANG QingHua2, TANG RongJiang1     
1. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
2. School of Earth and Space Sciences, Peking University, Beijing 100871, China
Abstract: The receiver function extracted from seismic three-component waveform data is often affected by source complexity and random noise, and some waveform anomalies often appear, which need to be eliminated before data interpretation. When there are many receiver functions, it will take a lot of time to select the normal receiver function artificially. In order to select high-quality receiver function efficiently, this paper used convolutional neural network (CNN) to judge the quality of receiver function. We used 9833 different earthquake events received in North China and Abaga area for training, and 1521 earthquake events were used for validation. The accuracy and recall rate of the training datasets are above 99%, and the accuracy and the recall rate of the test datasets are 95.3% and 92.4%, respectively. We also used the CNN to evaluate other data outside the training and the test datasets to get different types of waveforms, which proves that the evaluation results correspond well with the expected receiver function waveform. In addition, however, CNN cannot efficiently evaluate some receiver function which have good consistency but do not meet conventional standard of good quality. In order to solve this problem, this paper calculates the 2-norm of the receiver function within a certain azimuth and epicentral distance, and makes statistics on the results with lower 2-norm, and compares the results with the CNN output to select final qualified data.
Keywords: Receiver function    Convolution neural network    Quick selection    
0 引言

接收函数是利用远震P波波形的垂直分量对径向分量和切向分量作反褶积,直接从R分量中提取出地下间断面处由P波转换而来的S波信号,从而获得台站正下方的S波速度结构(Langston, 1979; Owens et al., 1984; 刘启元等,1996).该方法由于工作原理简单,经济实用,探测深度大,以及较高的垂向分辨率被广泛应用于地球深部速度间断面的研究.

在利用接收函数研究地下结构时,高质量的接收函数是得到可靠结果的重要前提.目前,挑选高信噪比接收函数主要以人为选取为主, 但是这种方法耗时巨大,且容易受到主观因素的影响.此外, 通过设计一套程序来完成这件事是很困难的,因为人脑中对理想接收函数的认识是一个“模糊”的模型,这个模型难以通过数学公式来进行精确衡量.尽管我们可以根据一些准则在程序中对数据质量进行判别, 例如,标准接收函数的初至直达P波的振幅最大;同一个台站的接收函数有较好的一致性;正确的接收函数的主要频段有既定的特征等(朱智昊和王毅,2019).然而实际得到的波形千差万别,且通过上述方案需要设定较多人为参数,最终可能导致数据挑选错误或数据漏选.

深度学习是机器学习研究的一个重要领域.与传统方法相比,深度学习完全是数据驱动的,其通过学习可以自动找到提取图像特征的最佳方法.卷积神经网络是深度学习的代表算法之一,也是首个被真正成功训练的深层神经网络(LeCun et al., 2015),该网络在目前得到了广泛的应用.它由用于特征提取的卷积层、池化层和用于特征处理的全连接层组成.卷积层由一系列特征图组成,每个特征图上的所有神经元共享同一个卷积核的参数,它由卷积核对前一层输入图像做卷积运算得到.池化层又称采样层,其实质是在局部区域进行最大值采样或平均值采样,目的是在减少数据量的同时保留有效信息.

接收函数的质量评估本质是波形图像的识别问题,而CNN在图像处理方面表现出非常优异的性能(Krizhevsky et al., 2012; LeCun et al., 1998; Simonyan and Zisserman, 2014; Szegedy et al., 2015; Zeiler and Fergus, 2014),这为本文使用CNN快速挑选接收函数奠定了基础.与图像、视频、文本和语音处理等其他科学领域相比,深度神经网络在天然地震中的应用仍处于起步阶段.在迄今发表的相关研究中,大多使用了中等复杂但功能强大的网络结构.例如,在20世纪90年代,Wang和Teng(1995)使用具有一个隐藏层的多层感知器进行事件检测和相位拾取;Musil和Plešinger(1996)利用神经网络进行了微震检测,并从宽频带波形记录中提取出微震记录;Del Pezzo等(2003)发表了一项关于识别地震和水下爆炸的研究,仅使用了三层神经网络;Böse等(2008)使用全连接神经网络提取了地震信号中的P波到达时间和传播速度.在2010年前后,科学家们在深度神经网络中的梯度消失和梯度爆炸方面的研究有了重大突破(Glorot and Bengio, 2010),这极大地推动了卷积神经网络和循环神经网络的发展:Zhu和Beroza(2019)利用U型深度卷积神经网络进行了P波和S波到时的提取;Ross等(2018)使用三个卷积层的CNN用于P波到时拾取和初动方向判断,其研究显示经过大规模数据训练的CNN模型对地震事件的拾取非常敏感,即使存在较强的噪声信号.Perol等(2018)基于CNN,使用简单的波形信号实现了地震的粗定位,该方法将地震定位看作分类问题处理,所得的定位结果只能分布在少数区域.此后,Kriegerowski等(2019)利用CNN实现了地震精定位,在这个过程中地震定位被看作回归问题进行处理.

然而,利用神经网络方法对接收函数进行快速准确挑选到目前为止仍没有相关报道.本文首先讨论了接收函数的质量评价标准,然后构建CNN模型,用9833个接收函数数据对模型进行训练,1521个数据进行测试, 并对比了不同CNN评估结果所对应的波形图,实验证明评估结果与实际的接收函数波形对应良好.该方法极大的提高了接收函数的挑选效率,减少了人为挑选接收函数的主观性.

1 接收函数的基本原理和质量评价标准 1.1 基本原理

获取接收函数的先决条件是远震事件,原因在于当远震P波到达台站下方时,可以被近似地看作平面波陡角入射到地表, 此时P波在垂直方向占优,而由P波在速度分界面上形成的一系列转换横波(例如:Ps,PpPms, PpSms,PsPms, PsSms等)在水平方向占优,这为区分P波和S波创造了有利条件.

地震三分量记录是由震源时间函数,传播路径,台站下方介质结构,仪器响应等多种因素共同作用的结果,为了得到台站下方介质的响应,考虑将地表三分量位移近似表示为

(1)

其中D(t)为远震P波波形数据,S(t)表示震源时间函数,I(t)表示仪器响应,Ev(t)、ER(t)、ET(t) 分别代表地下介质结构脉冲响应的垂直、径向和切向分量,符号*表示褶积运算.

许多波形简单的地震事件和理论计算表明,远震事件地表位移的垂直分量表现为尖脉冲的时间函数与仪器响应的褶积,紧随其后的续至震相非常小(Burdick and Langston, 1977; Langston, 1977),因此可以将介质结构响应的垂直分量近似看作Dirac函数.即可以得到:

(2)

因此,利用远震P波波形的垂直分量Dv(t) 对径向分量DR(t) 和切向分量DT(t) 作反褶积运算可以提取出转换S波信号ER(t)和ET(t),即P波接收函数,它消除了震源和射线路径、以及仪器的影响,可用来进一步获得台站正下方的S波速度结构.

1.2 接收函数的质量评价标准与影响因素

本文主要考虑P波径向接收函数.由于远震P波并不完全垂直入射,在径向还存在直达P波分量,且该震相明显强于转换S波震相,因此,结合1.1节的认识可知:P波径向接收函数主要包含较强的直达P波分量,以及台站下方速度间断面产生的转换及多次反射S波分量.因此,对于P波径向接收函数波形信号的质量判断,有两个重要准则:1存在振幅较强的直达P波信号以及振幅明显减弱的多次波和转换波信号,且在很多情况下,能够观察到来自莫霍面的转换Ps波震相,该震相大致位于5 s附近,通常来说Ps波震相振幅明显强于除直达P波以外的其余震相.2同一方位角和震中距范围内的波形有很好的一致性.

理论上来说,如果将介质的横向非均匀性看成一个静态的过程,除非地震波沿射线路径传播时恰好发生如流体活动等介质的瞬间变化,相同路径传播的远震应该有相同的接收函数,并且得到的接收函数都满足以上两个准则.然而,在实际情况中,并不是同一个台站的每个远震事件都能得到高质量的接收函数,部分接收函数不满足准则1和2,其原因大致可以归结为:①仪器工作的不稳定性.②震源的复杂性.远震事件受到台站附近小地震的影响,也就是远震事件叠加了小地震,这使得三分量记录的垂直分量无法近似为Dirac函数,因此直接使用远震P波的垂直分量对径向分量进行反褶积计算不能完全消除震源和仪器的影响.③反褶积方法受到数据噪声的影响较大.由于实际的远震数据是有限带宽的,且含有各种噪声,这使得最终得到的反褶积结果存在较大的干扰信号.

除此之外,还存在部分特殊接收函数,这些接收函数通常来自盆地内部的台站,由于盆地内部巨厚低速沉积层产生的多次波严重干扰了接收函数的有效信号,导致这些接收函数波形满足准则2,却不满足准则1,但它们仍然是地下结构的真实反应,是需要被挑选出来的质量合格的接收函数.然而,在训练集中该类数据被定义为“质量差”,这是因为该类数据仅对特定的台站才被认为质量好,而对于其他绝大部分台站而言属于质量差;且CNN模型判断数据质量时,输入参数仅为接收函数波形信号,无法添加有效的台站信息.

因此,本文在对接收函数进行质量评估时,首先使用训练好的CNN模型挑选出满足准则1的数据,实际上,满足准则1的数据通常也满足准则2.对于盆地内的数据,由于CNN模型的评估结果往往绝大部分为质量差,本文使用统计相关性算法对这部分数据进行处理,挑选出有效数据,以防止有效数据的漏选.

2 CNN挑选接收函数方法 2.1 数据准备

本文使用的接收函数数据大部分来源于中国地震局地球物理研究所在2006年6月—2009年10月在华北地区布设的190套宽频带地震观测仪器的记录数据,为了提高模型的泛化能力,同时选取了部分来自于中国地震局地球研究所在2017年6月—2018年10月在内蒙阿巴嘎地区布设的7套宽频带地震仪(如图 1).

图 1 数据来源的台站分布 1 宽频带台站;2 图 6图 7中被提及的台站信息;3 主要城市;4 华北研究区域;5 阿巴嘎研究区域. Fig. 1 Station distribution in the study region 1 Broadband stations; 2 The stations mentioned in figure 6 and figure 7; 3 Main cities; 4 Research area in North China; 5 Research area in Abaga.

选取震级在5.5级以上, 震中距位于30°~90°区域的远震事件11354个(如图 2).并对三分量波形数据进行预处理:去均值,去趋势,去倾斜,滤波处理,然后将两个水平分量旋转为径向分量和切向分量,最后利用远震P波的垂直分量对径向分量和切向分量作时间域最大熵反褶积(吴庆举等, 2003a, 2003b)来提取接收函数,本文中我们只使用径向分量的接收函数.

图 2 地震事件震中分布 Fig. 2 Epicenter distribution map of earthquake events

按照一定的比例将数据分为训练集和测试集,其中训练集9833个,测试集1521个.训练集中质量好的接收函数有2687个,质量差的7146个;测试集中质量好的数据有582个,质量差的939个.模型的输入为接收函数波形记录,输出为人为标定的质量好坏的标签(0为质量差,1为质量好).

选取直达P波到达后的40 s作为接收函数的有效时窗,并将直达波到达前2 s归零,接收函数的采样频率为10 Hz,采样点419个.由于接收函数的值通常分布在-0.2~0.6之间,分布范围较小,因此可以免去输入数据批量归一化的步骤,因为批量归一化的目的是为了防止输入数据方差太大引起梯度消失或梯度爆炸现象(Ioffe and Szegedy, 2015).在训练中,为了提高模型的泛化能力和降低过拟合的风险,使用随机批量梯度下降方法(Keskar et al., 2016),即每次随机选择100个数据进行训练,每次迭代需要训练98次.

2.2 CNN的构建

本文使用tensorflow机器学习开源库(https://www.tensorflow.org)构建经典的卷积神经网络,该模型由输入层、输出层,两个卷积层,一个池化层,和全连接层组成(图 3), Kriegerowski等(2019)曾利用相同的模型成功实现了地震的精定位.模型输入层为419个采样点的接收函数信号,为了便于计算和理解,将其转化为419×1×1的三维数组,输入层、卷积层和池化层均为三维数组,全连接层为一维数组.两个卷积层中滤波器的内核尺寸为3×1,移动步长分别为1和2,层厚分别为32层和64层;池化层采用最大池化方法,内核尺寸为3×1,移动步长为2,层厚为64层.CNN模型可以在多层结构中形成比较高阶的图像特征,同一个滤波器在每一层的卷积计算中权值共享.因此,相比于全连接神经网络,CNN模型的卷积层和池化层进一步提取了图像的多个特征,有效减小了计算机负载,提高了图像的识别能力.

图 3 用于挑选接收函数的卷积神经网络的系统结构 Fig. 3 The system structure of convolutional neural network for selecting receiver function

池化后的数据被重新转化为一维数组,进入全连接层,最终输出两个神经元,其最大值所在的序列(0或1)为输出标签.此外,我们对每层的输出(池化层除外)施加了Relu激活函数来提高模型的非线性能力,Relu激活函数的非饱和性(值域为[0, +∞])被证明能有效减少梯度消失和梯度爆炸的问题,同时计算速度也很快(Glorot et al., 2011).在训练过程中,所有的卷积层中的滤波器内核和全连接层中的神经元权重利用误差的反向传播进行更新,下降方法采用ADAM算法(Kingma and Ba, 2014),学习率使用0.001,迭代次数为100次.一旦模型训练完成,1521个测试集的预测仅需要0.1 s.可以参考Buduma和Locascio(2017)进一步了解CNN算法细节.

目标函数(或误差函数)采用Softmax-Cross-Entropy函数,具体表达为(Buduma and Locascio, 2017)

(3)

其中L为目标函数,实际计算中使用的目标函数是多个L的平均值,y′i为给定标签(0或1,分别代表质量差和好),yi为神经元的输出归一化后的结果,其各个分量的和为1,记作:

(4)

其中xj为第j个神经元的输出(ij取整,范围0~1), yi∈[0, 1].CNN进行数据评估时,通常以0.5为阈值,y1>0.5时数据质量为好,反之质量为差.

3 数据验证 3.1 质量评估结果

图 4a展示了训练集和测试集的平均误差函数随迭代次数的下降曲线,可以看出,在60次迭代之后,两者的误差函数不再下降,且训练集比测试集明显能达到更低的误差值.

图 4 误差函数,准确率和召回率演化 (a) 训练集和测试集的平均误差函数变化曲线; (b) 训练集和测试集的准确率的变化曲线; (c) 训练集和测试集的召回率的变化曲线; (d) 不同dropout的测试集平均误差函数变化曲线. Fig. 4 Evolution of loss function, accuracy and recall rate (a) Evolution of average loss function for training and test datasets; (b) Evolution of accuracy for training and test datasets; (c) Evolution of recall rate for training and test datasets; (d) Evolution of average loss function with different dropout rate for test datasets.

准确率和召回率是衡量CNN模型好坏的重要指标.准确率是指挑选出来的接收函数中,质量好的数据所占比重;召回率是指所有的质量好的接收函数中,被挑选出来的数据所占比重.图 4bc分别为测试集和训练集的准确率和召回率随迭代次数的变化曲线.图 4b显示了测试集的准确率随迭代次数的增加逐渐增加,最终趋于稳定,且测试集的准确率最终略低于训练集,但是仍达到95%以上;图 4c显示了两者的召回率都随着迭代次数的增加逐渐增大,与测试集准确率曲线相比,测试集召回率的曲线在上升的过程中更加不稳定,但最终仍能达到92%以上.通常情况下,准确率和召回率无法同时达到最佳,较高的准确率往往以较低的召回率为代价,反之亦然,因此,我们在实际计算中通常取两者的折中.图 5更加直观的展示了CNN模型测试集的预测结果.

图 5 CNN模型测试集预测结果与实际情况对比 Fig. 5 Comparison between predicted results of CNN model test dataset and true results

为了降低CNN模型的过拟合风险,在每层神经元中使用正则化方法——dropout算法(Srivastava et al., 2014),该算法的实质是在训练中通过忽略部分神经元,使得任意局部数据不会对整个模型过于敏感,从而达到正则化的目的.Dropout算法被证明具有较高的成功率,即使是最先进的神经网络在运用了dropout算法之后也能提高1%~2%的准确率.对于本文的训练数据,选取忽略不同比例的神经元并绘制出平均误差函数的下降情况于图 4d中,可以看出:黑色加粗曲线(忽略20%神经元)下降程度较多且更加稳定, 而橙色曲线(忽略0%神经元)和黄色曲线(忽略10%神经元)的误差曲线下降过程不稳定,出现明显的跳跃,这说明dropout算法一定程度上提高了模型的收敛性和稳定性.因此,本文选取忽略20%神经元的方案.

为了直观地展示CNN对接收函数的挑选结果,对测试集的数据评估结果y1(公式(4))进行统计,y1越大(越接近于1),代表数据质量越好,反之质量越差,统计结果如图 6a.其横坐标为由y1等间隔分割的10个区间;纵坐标为接收函数的统计数量(由于不同y1区间的接收函数数量差异较大,因此纵坐标采用对数显示);黄色柱状图表示实际质量好的数据,蓝色柱状图表示实际质量差的数据,为了方便展示,我们将数量较少的直方图叠加在数量较多的直方图上面.由图 6a可以看出:绝大部分数据的评估结果都位于y1区间的两端(数量级102),即:接近于0或1;少部分位于y1中间区域(数量级10以内).且大部分数据的实际质量和CNN评估结果对应良好,即:CNN评估结果y1>0.5的数据大部分呈黄色,反之y1<0.5的区间对应的数据大部分呈蓝色.

图 6 CNN评估结果与误差统计直方图 (a) CNN评估结果统计直方图; (b) 实际结果与CNN选取结果的误差统计直方图. Fig. 6 Evaluation results of CNN and error statistical histogram (a) Statistical histogram of CNN evaluation results; (b) Error statistical histogram of true results and CNN outputs.

图 6b展示了实际数据质量与CNN评估结果的误差,横坐标为误差, 其计算公式为

(5)

其中xi为第i(i=0, 1)个输出层神经元的值,纵坐标为接收函数的数量(对数间隔).从图 6b中可以看出,大部分接收函数的误差在0.1以内,误差在0.5以上的数据仅70个,占比4.6%,说明了该方法的准确率较高.

3.2 波形对照

为了直观展示评估结果与实际接收函数波形的对应关系,选取阿巴嘎地区AB24台站(如图 1)一定震中距和方位角范围内的数据进行质量评估,数据总数455个,地震事件均位于经度100°~150°,纬度-10°~10°范围内,且这部分数据未被纳入训练集和测试集.根据评估结果y1,我们将数据分为三类:① y1∈[0.9, 1.0], ②y1∈[0.1, 0.9], ③y1∈[0, 0.1], 并将这三类数据分别叠放在一起,如图 7所示.

图 7 AB24台站的不同CNN评估结果对应的455个接收函数叠放图(图中曲线颜色无实际意义) (a) 所有数据的波形整体显示; (b) 类型①数据的叠放图, 即CNN评估结果y1∈[0.9, 1.0]; (c) 类型②数据的叠放图,即评估结果y1∈(0.1, 0.9); (d) 类型③数据的叠放图,即评估结果y1∈[0, 0.1]. Fig. 7 Superposition diagram of receiver functions in AB24 station with different CNN evaluation results.(The color of the curves in the figure has no practical significance) (a) The overall display of the waveform of all data; (b) The waveform superposition graph of y1 ∈[0.9, 1.0]; (c) The waveform superposition graph of y1 ∈(0.1, 0.9); (d) The waveform superposition graph of y1∈[0, 0.1].

图 7a为选取的所有数据的波形结果,可以看出,不同质量的数据混合在一起,无法有效区分出数据质量的好坏,且绝大部分波形的振幅位于-0.4~0.8之间.图 7b展示了y1∈[0.9, 1.0] 的数据,一共55个,从图中可以看出该部分接收函数的直达P波的振幅最大,来自莫霍面的Ps震相明显,且具有较好的一致性,是比较理想的接收函数.

值得注意的是图 7c, 这部分数据的y1∈[0.1, 0.9],一共16个.该部分接收函数波形的振幅范围总体上比较正常,直达波振幅明显,部分数据存在明显的Ps震相,但是一致性不如图 7b, 且存在少部分数据直达P波振幅过大,接近1,因此该部分数据总体上质量属于中等.存在这部分数据的原因是人为挑选的作为训练集的标签存在一定的主观性,这种主观性来源于部分难以定性的“质量不好也不差”的接收函数,这些数据被随意安置了标签,这类数据对应的评估结果y1通常位于[0.1, 0.9].

图 7d展示了被评估为质量差即y1∈[0, 0.1]的384个数据,可以看出这部分接收函数的波形十分杂乱,无明显的Ps震相,且一致性较差,部分甚至完全偏离了正常的接收函数形态.

综上可以看出,CNN模型的大部分评估结果位于1或0附近,分别对应着质量较好和较差的数据, 少部分数据的评估结果位于0.1到0.9之间,且CNN评估结果与人为预期对应良好.数据从①到③,质量整体上逐渐变差,这说明y1在一定程度上能够定量评估数据质量,尽管也存在少数波形例外(例如图 7c中明显存在质量较好的数据),但是质量好的数据在①到③类的数据中所占比例逐渐下降.在数据的挑选过程中,我们重点关注①类数据,当需要的接收函数数量不够时,我们可以进一步对②类数据进行筛选.为了进一步提升模型,还可以人为挑选出错误分类的数据,重新给定标签,加入训练集进一步做增量训练.

4 部分特殊的接收函数处理方法

正如前文所提到的,CNN无法挑选出一致性很好却不符合常规标准的接收函数,例如盆地内的数据.对于这部分数据,我们利用波形的一致性对数据进行挑选,具体步骤如下:

(1) 从原始波形记录中提取接收函数之后,将每一个台站的数据按照远震事件的经纬度分组,将经纬度相近的数据分为一个组,以保证每组数据具有近似的方位角和震中距,以及相同传播路径.

(2) 在每组数据中,两两之间相互求取二范数,计算公式为

(6)

其中misfit越小,表示两个波形越接近,相似程度越高.将每组数据中所有misfit的平均值乘以0.55作为阈值,统计出低于该阈值的数据,若某个数据与多个其他数据计算的misfit都低于阈值,则它被多次统计.这里的阈值0.55为多次试验后所得.最后得到该组中所有数据被统计次数的分布结果,统计次数较多的数据之间必然存在较好的一致性.

(3) 挑选出统计次数较高的数据,同时使用训练好的CNN模型对数据进行评估,如果这些数据的CNN评估结果为差,即是需要被特殊处理的接收函数,它们具有较好的一致性却不符合常规标准.对于其他情况,则按普通接收函数进行处理.

为了对上述方法进行验证, 选择两个位于华北地区盆地内部的台站A007和A010(位置如图 1)的数据进行验证,所有地震事件均位于经度100°~150°,纬度-10°~10°范围内,且该数据不在训练集或测试集以内.按照上述步骤对两个台站的数据进行处理,并根据一致性算法统计结果进行排序,展示排序结果的前80个数据于图 8ab.图 8中右端的红色柱状图为归一化之后的一致性统计结果,值越大表示一致性越好,柱状图按照统计结果从上到下逐渐增大.蓝色菱形散点图为CNN模型评估结果,即y1(公式(2)),值越大表示越接近标准的接收函数.

图 8 华北地区两个台站(a A007 b A010)的接收函数波形图,以及对应的一致性统计和CNN模型评估结果, 所有地震事件来源于经度100°~150°,纬度-10°~10°范围内. Fig. 8 Receiver function waveforms of two different stations in North China (a A007 b A010) as well as corresponding consistency statistics and CNN evaluation results. All events are in the range of longitude 100°~150° and latitude -10°~10°.

可以看出,图 8a图 8b中的前30道波形数据,对应着较大的一致性统计结果,一致性较好.此外,该部分数据大部分在直达波之前出现了明显的振幅较小的波形, 两者的CNN评估结果(蓝色菱形)绝大部分都趋近于0,其中图 8b中存在三个数据的CNN评估结果大于0.5,它们的波形图明显不符合常规标准,说明这三个数据被CNN错误估计,误差在4% 以内.综上都说明图 8中的接收函数是明显的需要特殊处理的接收函数,并且一致性统计算法能有效地挑选出需要的数据.

尽管在一致性统计结果较大的数据中(前30道),也存在小部分数据与其余数据一致性差异较大,例如图 8a中的15道,图 8b中的13和21道,这些数据需要人为进一步剔除;但是一致性统计算法仍能获取大部分一致性较好的数据,这些数据无法通过CNN得到,因此,一致性统计算法可以很好地弥补CNN模型评估这种数据的不足.

总的来说,结合CNN方法和一致性分析能确保快速挑选出几乎所有的合格的接收函数.值得注意的是,一致性分析中需要人为设定misfit的统计阈值(本文经验得出的平均值的0.55),对于不同的组,该阈值大小需要调整,且会极大的影响统计次数,进而影响挑选结果,因此一致性算法仅可以作为CNN评估结果的辅助分析.

5 结论

深度学习卷积神经网络方法在越来越多的地震数据处理上展现出极大的发展潜力.本文展示了将CNN应用于接收函数质量好坏的评定,仅仅用了较少的训练集数据,就在测试集和训练集的准确率和召回率之间取得很好平衡,训练集的准确率和召回率达到了95%以上,测试集也达到92%以上.我们也对比了CNN评估结果与接收函数波形图,实验证明与预期比较吻合.对于沉积盆地地区某些接收函数一致性较好但是不符合常规标准的数据,CNN无法识别的问题,本文还引入了一致性统计算法进行补充,从而避免了部分数据的漏选.结合两种方法能确保快速挑选出几乎所有的合格的接收函数,这极大提高了人为筛选接收函数的效率,因此该方法具有较大的应用前景.

致谢  非常感谢审稿专家提出的宝贵意见.
References
Böse M, Wenzel F, Erdik M. 2008. PreSEIS: A neural networkbased approach to earthquake early warning for finite faults. Bulletin of the Seismological Society of America, 98(1): 366-382. DOI:10.1785/0120070002
Buduma N, Locascio N. 2017. Fundamentals of Deep Learning: Designing Next-Generation Machine Intelligence Algorithms. Sebastopol: O′Reilly Media.
Burdick L J, Langston C A. 1977. Modeling crustal structure through the use of converted phases in teleseismic body-wave forms. Bulletin of the Seismological Society of America, 67(3): 677-691.
Del Pezzo E, Esposito A, Giudicepietro F, et al. 2003. Discrimination of earthquakes and underwater explosions using neural networks. Bulletin of the Seismological Society of America, 93(1): 215-223. DOI:10.1785/0120020005
Glorot X, Bengio Y. 2010. Understanding the difficulty of training deep feedforward neural networks.//Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. 249-256.
Glorot X, Bordes A, Bengio Y. 2011. Deep sparse rectifier neural networks.//Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics (AISTATS). Fort Lauderdale, USA: JMLR, 315-323.
Ioffe S, Szegedy C. 2015. Batch normalization: Accelerating deep network training by reducing internal covariate shift.//Proceedings of the 32nd International Conference on Machine Learning. Lille, France: JMLR. org.
Keskar N S, Mudigere D, Nocedal J, et al. 2017. On large-batch training for deep learning: Generalization gap and sharp minima.//Proceedings of the 5th International Conference on Learning Representations. Toulon, France: OpenReview.
Kingma D P, Ba J. 2014. Adam: A method for stochastic optimization. arXiv preprint arXiv: 1412.6980.
Kriegerowski M, Petersen G M, Vasyura-Bathke H, et al. 2019. A deep convolutional neural network for localization of clustered earthquakes based on multistation full waveforms. Seismological Research Letters, 90(2A): 510-516. DOI:10.1785/0220180320
Krizhevsky A, Sutskever I, Hinton G E. 2012. ImageNet classification with deep convolutional neural networks.//Proceedings of the 25th International Conference on Neural Information Processing Systems. Red Hook, NY, USA: Curran Associates Inc., 1097-1105.
Langston C A. 1977. The effect of planar dipping structure on source and receiver responses for constant ray parameter. Bulletin of the Seismological Society of America, 67(4): 1029-1050.
Langston C A. 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. Journal of Geophysical Research: Solid Earth, 84(B9): 4749-4762. DOI:10.1029/JB084iB09p04749
LeCun Y, Bottou L, Bengio Y, et al. 1998. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11): 2278-2324. DOI:10.1109/5.726791
LeCun Y, Bengio Y, Hinton G. 2015. Deep learning. Nature, 521(7553): 436-444. DOI:10.1038/nature14539
Liu Q Y, Kind R, Li S C. 1996. Maximal likelihood estimation and nonlinear inversion of the complex receiver function spectrum ratio. Acta Geophysica Sinica (in Chinese), 39(4): 442-443.
Musil M, Plešinger A. 1996. Discrimination between local microearthquakes and quarry blasts by multi-layer perceptrons and Kohonen maps. Bulletin of the Seismological Society of America, 86(4): 1077-1090.
Owens T J, Zandt G, Taylor S R. 1984. Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee: A detailed analysis of broadband teleseismic P waveforms. Journal of Geophysical Research: Solid Earth, 89(B9): 7783-7795. DOI:10.1029/JB089iB09p07783
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, 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
Simonyan K, Zisserman A. 2014. Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv: 1409.1556v4.
Srivastava N, Hinton G, Krizhevsky A, et al. 2014. Dropout: a simple way to prevent neural networks from overfitting. The Journal of Machine Learning Research, 15(1): 1929-1958.
Szegedy C, Liu W, Jia Y Q, et al. 2015. Going deeper with convolutions.//Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. Boston, MA, USA: IEEE, 1-9.
Wang J, Teng T L. 1995. Artificial neural network-based seismic detector. Bulletin of the Seismological Society of America, 85(1): 308-319.
Wu Q J, Tian X B, Zhang N L, et al. 2003a. Receiver function estimated by Wiener filtering. Earthquake Research in China (in Chinese), 19(1): 41-47.
Wu Q J, Tian X B, Zhang N L, et al. 2003b. Receiver function estimated by maximum entropy deconvolution. Acta Seismologica Sinica (in Chinese), 25(4): 382-389.
Zeiler M D, Fergus R. 2014. Visualizing and understanding convolutional networks.//13th European Conference on Computer Vision. Zurich, Switzerland: Springer, 818-833.
Zhu W, Beroza G C. 2019. PhaseNet: a deep-neural-network-based seismic arrival-time picking method. Geophysical Journal International, 216(1): 261-273.
Zhu Z H, Wang Y. 2019. Discussion and application of criterions of picking receiver functions. Computing Techniques for Geophysical and Geochemical (in Chinese), 41(5): 580-585.
刘启元, Kind R, 李顺成. 1996. 接收函数复谱比的最大或然性估计及非线性反演. 地球物理学报, 39(4): 442-443.
吴庆举, 田小波, 张乃铃, 等. 2003a. 用Wiener滤波方法提取台站接收函数. 中国地震, 19(1): 41-47.
吴庆举, 田小波, 张乃铃, 等. 2003. 计算台站接收函数的最大熵谱反褶积方法. 地震学报, 25(4): 382-389. DOI:10.3321/j.issn:0253-3782.2003.04.005
朱智昊, 王毅. 2019. 接收函数挑选准则的探讨及成像效果分析. 物探化探计算技术, 41(5): 580-585. DOI:10.3969/j.issn.1001-1749.2019.05.03