2019, Vol. 62
2. 北京大学地球与空间科学学院, 北京 100871
2. School of Earth and Space Sciences, Peking University, Beijing 100871, China
很多地震学问题都依赖于地震震相的准确识别.随着地震资料积累速度的快速增长, 人工震相识别已难以高效地处理这些资料, 迫切需要进行自动识别.
在早期的探索中, 人们提出了不同的震相自动识别方法, 例如长短窗(Allen, 1978)、自回归(Sleeman and Van Eck, 1999; Leonard and Kennett, 1999)和kurtosis(Saragiotis et al., 2004)等.之后也有人联合使用几种方法进行震相自动识别(如尹得余, 2012).这些震相自动识别方法都是利用一个或几个统计特征进行的.尽管人们十分谨慎地选取了参数, 但却很难同时降低误识别率和漏识别率.还有人利用分形技术确定初至(如李庆忠, 1996; 常旭和刘伊克, 1998; 曾富英等, 2002), 但只适用于简单情况下的震相拾取.
为了在尽量减少误识别的前提下增加对地震的检测能力, 很多学者都尝试使用模板识别来寻找地震目录中未被记录的地震(如Peng and Zhao, 2009; Wang et al., 2015; 吴梦羽等, 2018).这种方法大幅度地增加了检测到的地震数目, 但极大地依赖于已有地震模板.对于因种种原因而没有模板的某些位置的地震, 尽管信号很强, 也会被该方法漏识别.
近几十年来, 计算机技术迅猛发展, 使得很多机器学习算法从纸上谈兵走向了实际应用.日益成熟的机器学习算法和强大的计算能力, 提供了将过去积累的大量数据利用起来的可能.很早就有人尝试使用机器学习的思路来研究震相的自动拾取, 其中大多数研究通过选取多个特征, 并利用神经网络等方法依据这些特征进行判断(Wang and Teng, 1997; Bai and Kennett, 2000; 刘希强等, 2009; 王彩霞, 2014).近一两年, 有学者将深度学习方法用于直接从原始波形中提取特征以识别地震信号, 如Perol等(2018)基于卷积深度神经网络进行地震自动识别和分区.这些方法一方面具有比传统方法更高的准确率, 另一方面, 由于机器学习具有一定的泛化能力, 在一个地方的数据中训练得到的模型, 在另一个地区也可以继续使用, 是非常有潜力的方法.
实际上, 利用机器学习的方法在检测地震方面已经取得了初步成功, 但如何在复杂的实际情况下进行震相的准确识别及到时的精确提取依然是需要进一步研究的问题.支持向量机作为一种机器学习方法(见Cristianini and Shawe-Taylor, 2005), 能够很好地处理高维度、非线性的分类问题, 已被广泛应用于模式识别问题中(Bahlmann et al., 2002; 傅军栋等, 2016; 张舒雅等, 2017; Batta et al., 2018).本文设计了两个基于支持向量机技术的分类器SSD和SPS, 并利用这两个分类器对地震体波震相进行了自动识别和到时的自动拾取.为进一步提高准确率, 我们还利用台阵资料进行了约束.模拟数据和实际数据的测试结果表明了方法的正确性和实用性.
1 支持向量机方法简介支持向量机是一种基于统计学习理论的机器学习算法.一般地, 可以用n维向量x来表示对象的n个被量化了的特征(feature), 并用x(i)i=1, 2, …, m}来表示m个对象的特征向量(feature vector); 用y∈-1, 1表示两个不同的类别, 并用y(i)i=1, 2, …, m}来表示m个对象所属的类别.支持向量机考虑特征向量是处在n维空间中的m个点, 如果能够寻找到一个存在于此空间的超平面, 将分属于两类对象的点分开, 则可以根据特征向量x相对于超平面的位置建立到其对应的类别y的映射:
|
(1) |
其中, 当z> 0时, g(z)=1, 否则g(z)=-1; w为X空间内一个矢量, 是超平面的一个法向量, 和x的维数相同; T表示对w进行转置; b是一个标量, 确定了超平面的具体位置; wT·x+b的正负反映特征向量相对于超平面的位置.超平面确定得越好, hw, b(x(i))和相应的y(i)一致的比例越高.
一般地, 将寻找最好的超平面转化为求解以下条件极值问题:
|
(2) |
为了使用到更多数据的信息, 将限制放宽, 得到下面的条件极值问题:
|
(3) |
在这个极值问题中, y(i)(wTx(i)+b)的值可以小于1, 这个小于1的程度ξi也在加了权重C后被引入了求解极值的过程.可以采用SMO算法(Platt, 1999)快速求解上述问题, 得到最优化的分隔超平面.
如果面对的是一个非线性的分类问题, 可以将x映射到更高维的空间, 使其变为线性可分的.由于支持向量机训练和判断的时候只关心向量间的内积, 所以仅需要计算出投影后向量对应的内积, 从而减少计算代价.用于计算投影后向量内积的函数被称为核函数; 而在保证计算代价可接受的前提下, 将原本不可分的向量投影到更高维度以使其可分的方法称为支持向量机的核技巧(Cristianini and Shawe-Taylor, 2005).
本文主要研究初至P波和初至S波的拾取问题, 其关键在于两个二分类问题:在连续波形记录中, 对震相和噪音进行分类; 在拾取到的震相中对P波和S波进行分类.首先, 支持向量机在处理高维度、非线性的二分类问题上表现十分好, 如果能够构造出合适的特征向量, 很适合这两个二分类问题; 再者, 支持向量机对于样本数量的要求不高, 仅需要少量的人工标注数据, 就可以达到很好的分类效果, 这使得我们可以节省标注所需的人力成本(丁世飞等, 2011).其中, 最关键的是特征向量的构造, 需要基于对地震波的深入理解进行设计.
2 利用支持向量机识别震相和拾取震相到时本文把震相拾取分解为两个二分类问题:首先从连续的波形记录中对震相和噪音进行分类; 然后, 从震相中对P波和S波进行分类.我们设计了两个支持向量机来分别进行这两个分类:区分震相和噪音的支持向量机, 作为地震信号探测器SSD(Seismic Signal Detector); 区分P波和S波的支持向量机, 作为震相分离器SPS(Seismic Phases Separator).
在我们的训练过程中, 我们将9/10的输入数据用于训练, 剩下1/10的数据用于测试训练出来的机器:利用训练得到的支持向量机对未参与训练的数据进行分类, 并将分类结果和标注进行比较, 得到相应正确率; 当该正确率不再上升时, 终止训练, 并返回训练好的机器和相应的正确率.根据这个返回的正确率, 可以评估训练好的机器的识别效果.本文尝试不同的特征向量、权重因子C和核函数形式, 选取具有较高正确率的方案.经过比较, 最终选定C值和核函数形式如下:
|
(4) |
|
(5) |
其中, p(i)和p(j)表示x(i)和x(j)投影到高维空间后对应的向量.
2.1 地震信号探测器(SSD)如前所述, 支持向量机效果好坏极大地取决于特征向量的构造是否合理.由于地震波波到达前、到达时和通过后地震台站波形记录有明显的差别, 本文中设计的特征向量同时包括了三个相应时间段的波形信息.因为地震波到达时的波形是本文关注的重点, 且水平方向和垂直方向的波形特征有所不同, 所以在特征向量中保留了该时段水平和垂直方向高采样率的波形记录.考虑当样本中地震的方位分布不均匀时, 水平的两个方向上由于方位角变化造成的能量分配差异可能会降低泛化能力, 所以将水平两个方向的波形按均方根进行合并.同时, 为了降低模型的复杂度以提高泛化能力并节约计算量, 将地震波到达前和通过后的地震波形三分量合并, 并进行了降采处理.
根据以上考虑并经过大量试验, 本文提出利用单台的全部三道信号构造特征向量.对于一个采样率为50 Hz的连续记录的每一个时间点, 本文中将它的特征用四个部分表示(见图 1):第一部分为这个时间点前后共3 s的一个时间窗内每一个采样点上水平运动的矢量和的模, 共151个; 第二部分为与第一部分相同的窗内, 垂直分量的波形记录, 共151个; 第三部分为考察的时间点之前40个窗内的振幅均方根, 涵盖时间点前160 s的平均振幅信息; 第四部分为考察的时间点之后40个窗内振幅均方根, 涵盖时间点后160 s的平均振幅信息.这四部分, 共382个点共同描述了研究的时间点的波形特征.为了减少振幅的影响, 将特征向量归一化, 使得前三部分302个值的平方和为1.
|
图 1 SSD的特征向量示例 特征向量由A、B、C、D四个部分构成, 横坐标N表示特征向量所在的高维空间各个维度的序号, 纵坐标Vf为对应特征归一化后的值. Fig. 1 Example showing SSD′s feature vector The feature vector is composed of four different parts, i.e., A, B, C, and D.The abscissa N shows the serial number of each dimension of the high-dimensional space where the feature vector is located, and the ordinate Vf is the normalized value of the corresponding feature. |
设定对于特征值x(i), 如果是震相, y(i)=1;如果是噪音, y(i)=-1.训练集由余震捕捉AI大赛提供的人工拾取的1763条震相和随机抽取的噪声6237条共8000条构成.用由震相和噪音共20715个例子构成的测试集对训练好的支持向量机进行测试.测试结果表明, 分类的准确率为96.89%.
由于SSD的输出值(wTx+b)分布在一个很广的范围内, 为了在不改变判断结果的前提下, 本文使用了如下形式(sigmoid函数)将输出值映射到一个更小的区间里(Gibbs and Mackay, 2000)以方便表达:
|
(6) |
其中, p(y|x)表示对象x属于y分类的概率; 当wTx+b>0, y=1的概率大于0.5且会随着wTx+b的增大而逐渐趋于1, y=-1的概率则小于0.5并趋于0;反之亦然.对于SSD, 本文利用p(1|x)来代替支持向量机的输出值.
其取值频率(输出值在范围内的对象占总数的比值)分布如图 2所示.可以看出能够很好地区分震相和噪音.
|
图 2 SSD对震相和噪音输出值的频率分布图 横坐标是SSD的输出值, 纵坐标表示相应取值的出现频率. Fig. 2 Frequency distribution of phases′ and noises′ output from SSD The abscissa is SSD′s output, and the ordinate indicates the output′s frequency. |
使用训练好的区分震相和噪音的支持向量机, 可以对连续信号的每一点进行分类并返回一个估值.可以选取估值中大于
考虑到在同一记录中, P波和S波射线在多数情况下入射角度接近, 但偏振方向有所不同(Cichowicz, 1993).同时, 在天然地震中, 同频率的S波比P波波长短、Q值低(Dziewonski and Anderson, 1981), 能量衰减更快, 而这种现象在高频部分体现得更加明显, 所以S波的主频通常会比P波低一些.当初至P波已确定, 通过比较后续震相的偏振方向及主频与初至P波的差异, 则有可能将后续震相区分为P波和S波.所以, 本文的SSD特征向量要同时包含初至P波和需作分类的震相的偏振和主频信息.
震相的偏振信息可以通过分析如下三分量波形记录的互相关矩阵的本征向量和本征值来提取(Jurkevics, 1988):
|
(7) |
其中si表示地震信号的第i个分量.当然, 除了偏振信息, 互相关矩阵还包含了极化程度等震相的其他信息, 也能够帮助我们对P波、S波进行分类.为了尽可能地保留互相关矩阵所包含的震相信息, 本文直接将互相关矩阵作为特征向量的一部分, 而不具体地计算出其偏振方向.
下面介绍引入震相主频信息的方法.首先将三分量记录中的第i个分量si和其导数s′i, 在[N0, N1]的部分做傅里叶展开, 可以得到
|
(8) |
|
(9) |
通过下式考察震相的主频特征:
|
(10) |
其中, ωk=2πk/N.(10)式中, 分母对信号的平方求和对应于信号互相关矩阵的对角元素; 分子对信号导数的平方求和对应于信号导数的互相关矩阵的对角元素.利用支持向量机的抽象能力, 将(10)式的分子和分母分别作为特征向量的一部分, 就可以将震相的主频信息包含进去了.信号的互相关矩阵已经包含在特征向量中, 还需要在特征向量中加入其导数的互相关矩阵.
有些地震的能量主要分布在低频, 这个时候对震相的主频估计有可能会被高频噪声等因素干扰, 对信号进行平滑会改善这种情况.本文先对s用长度为5的窗平滑后得到信号sw, 分别定义原始记录s和经过平滑的记录sw中, P波和其相应的S波的主频比值R和Rw如下:
|
(11) |
|
(12) |
其中, Tp0、Tp1表示对初至P震相计算主频的时窗起始和结束的时刻, Ts0、Ts1表示对S震相计算主频的时窗起始和结束的时刻.将不同记录计算得到的2666个R绘制到图 3a中, 其中有80.12%的R值大于1, 即大多数P波的主频率比S波高; 对同一P波和S波震相的R、Rw值, 将它绘制出来在图 3b中, 仅有346个记录两者同时小于1, 占总数的12.98%, 说明同时考虑原始信号和平滑后信号的主频信息后, 对P、S的分类效果更好.所以本文将平滑后的信号sw和其导数s′w的互相关矩阵补充到特征向量里面.
|
图 3 R值的频率分布直方图(a)和R-Rw散点图(b) (a)横坐标表示R值, 纵坐标为其对应的频率; (b)横坐标表示R值, 纵坐标为其对应的Rw. Fig. 3 Frequency distribution histogram of R (a) and R-Rw scatter plot (b) (a) The abscissa is the R value, and the ordinate indicates the corresponding frequency; (b) The abscissa is the R value, and the ordinate is the corresponding Rw value. |
定义D(s; T0, T1)表示三分量的地震信号s的互相关矩阵对应的向量:
|
(13) |
其中, [T0T1]表示计算互相关的时窗起始和结束的时刻; 互相关矩阵只有六个独立分量, 该向量的长度也为6.
我们用如上形式的向量构建SPS的特征向量, 包含八个部分, 分别表示了初至P波和该震相的信息.初至P波:{D(s; Tp0, Tp1); D(s′; Tp0, Tp1); D(sw; Tp0, Tp1); D(s′w; Tp0, Tp1)}; 该震相:{D(s; Tk0, Tk1); D(s′; Tk0, Tk1); D(sw; Tk0, Tk1); D(s′w; Tk0, Tk1)}.其中, Tp0, Tp1表示对初至P震相时窗起始和结束的时刻, Tk0, Tk1表示对需作分类的震相时窗起始和结束的时刻; s′, s′w表示s, sw对时间的导数.它们共同构成特征向量时, 需要分别地, 对初至P波的部分归一化使得D(s; Tp0, Tp1)的前三个元素和为1, 对该震相的部分归一化使得D(s; Tk0, Tk1)的前三个元素和为1.图 4展示了SPS的特征向量的两个例子.
|
图 4 P波震相和S波震相在SPS里的特征向量 横坐标N为特征向量所在的高维空间各个维度的序号, 纵坐标为对应特征归一化后的值; 第一道表示后续P震相的特征向量, 第二道表示S震相的特征向量; P0表示特征向量关于初至P波的部分. Fig. 4 P and S phases′ feature vectors in SPS The abscissa N indicates the serial number of each dimension of the high-dimensional space where the feature vector is located, and the ordinate is the normalized value of the corresponding feature; the first trace represents the feature vector of the following P-seismic phase and the second trace represents the feature vector of the S-seismic phase; P0 represents the part of the first arrival P phase of feature vectors. |
在SPS中, 设定对于特征值x(i), 如果是P震相, y(i)=1;如果是S震相, y(i)=-1;使用和SSD相同的核, 并用已知的2000个S震相和3000个后续P波震相进行训练, 得到训练好的SPS.用训练好的机器, 对P和S震相做同样的统计, 结果如图 5a.大部分P波和S波分列在0值的两边.由于存在震相相距较近的情况, 少部分震相有越界的情况.在由3999个后续P波震相和2666个S波震相构成的训练集中, 分类准确率达到了96.14%.Cichowicz(1993)曾利用互相关矩阵提取震相偏振信息和极化程度的方式对P波和S波进行区分.他的方法定义了F值表示震相的偏振方向和极化程度与初至P波的差异, F值越大其为S波的可能性越大.我们对于3999个后续P波震相和2666个S波震相的F值进行了统计, 其频率分布如图 5b.可以看出P波和S波计算得到的F值分布有一定差异, 但是很难找到一个明确的分界点, 分类效果明显不如本文方法.
|
图 5 P震相和S震相输出值的频率分布 (a) SPS对P波和S波震相输出值的频率分布.横坐标是SPS的输出值, 纵坐标表示相应的频率; (b)根据Cichowicz(1993)方法计算得到的F值频率分布.横坐标表示F值, 纵坐标为其对应的频率. Fig. 5 Frequency distributions of P and S phases′ output (a) The frequency distribution of P and S phases′ output from SPS.The abscissa is the output value of SPS, and the ordinate is the corresponding frequency; (b) P and S phases′ F values′ frequency distribution according to Cichowicz′s method (1993).The abscissa is the F value, and the ordinate is the corresponding frequency. |
利用两个支持向量机可以在单台上完成对震相的拾取:首先, 使用SSD, 拾取到一系列震相, 取一定时间窗内到时最早的震相作为初至P; 利用初至P, 结合后面震相的信息, 使用SPS, 将震相分为P和S两个部分.S波里面, 到时最早的震相作为S波的初至.拾取结果见图 6.
|
图 6 实际数据拾取结果示意图 上面3道波形, 分别是台站的三分量; 第4道为SSD的输出; 第5道为根据SSD的输出选出的存在震相的一系列时间点; 第6道为使用SPS从第5道的震相点中选出的S震相; 在第4道中, 用圆圈表示挑选出的初至P, 十字表示挑选出的初至S. Fig. 6 Schematic diagram of real data picking-up results The top three traces of waveform are the three components of a seismic station.The fourth trace is the probability whether there being a phase given by SSD.The fifth trace shows the series of phases′ arrival times selected based on the output given by SSD.The sixth trace shows the S phases selected from phases in the fifth trace using SPS.In the fourth trace, the first-arrival P is indicated by a circle and the first-arrival S is indicated by a cross. |
因为对实际数据的人工标注存在误差, 难以准确地测试本文自动拾取程序的性能, 所以本文采用FK程序(Zhu and Rivera, 2002)生成的模拟数据进行测试, 以理论到时为标准, 比较了本文的算法和几种经典的自动拾取方法拾取结果的准确度与精确度.
长短窗(Allen, 1978)、AR-AIC(Sleeman and Van Eck, 1999)和kurtosis(Saragiotis et al., 2004)是比较经典的可以进行P波拾取的算法.本文的SSD也可以单独地进行P波拾取.本文对500个添加了噪音的模拟地震记录(采样率50 Hz, 震源深度:10~60 km, 震中距:30~200 km), 利用上述几个算法对P波进行了拾取.
为了提高长短窗的识别率, 本文将三个分量的信号都进行了初至P波的拾取, 并取时间最早的作为初至P波的到时.长窗长度为100个采样点, 2 s; 短窗长度为5个采样点点, 0.1 s; γ阈值设定为5;在初至P波理论到时前后10 s, 共20 s内进行拾取最早超过阈值的时间点作为初至P波.
对于AR-AIC算法, 如果时间窗过长包含到其他后续震相, 会导致震相拾取不准确, 所以本文仅选取目标时间前后3 s, 共6 s的窗进行研究.选取似然性log(L)最好的分隔点作为P波初至.
对于kurtosis方法, 本文使用了50种组合的(M, q)(M=50, 55, …, 70, q=0.90, 0.91, …, 0.99)进行震相拾取, 并取其中出现频率最高的到时作为最终确定的到时.取初至P波理论到时前后10 s, 共20 s的窗进行拾取.
本文使用SSD, 对初至P波理论到时附近前后10 s, 共20 s的波形进行判断和拾取.拾取最早的一个震相作为初至P波的到时.在上述设定下进行拾取, 并对结果进行了统计, 对比结果见表 1.
|
表 1 本文中的方法与三种经典方法拾取P波震相的结果比较 Table 1 Comparison between the method proposed by this paper and three classic automatic P phase picking methods |
可以看出, 本文给出的基于支持向量机的算法能够比较准确地确定震相初至, 对模拟数据的拾取和理论到时误差在1 s以内的比例高达95.2%, 明显优于其他算法, 也就是说在拾取震相的正确率上本文的算法具有很大的优势.
比较平均误差和误差方均根, 可以发现, 本文的方法在准确确定起振点上效果不够好, 而AR-AIC算法在确定起振点上的效果最好.但是AR-AIC算法不管时间窗内有无震相, 都会给出一个到时, 实际上无法判断时间窗内是否含有震相, 必须依赖于其他方法来给定时间窗.于是本文将两个算法结合起来, 首先使用本文提出的支持向量机方法进行初步拾取, 然后再根据这个结果设定时间窗, 使用AR-AIC算法进行到时的精确拾取.
长短窗使用能量作为判据, 有时候会被一些非地震信号触发, 单独作为拾取方法并不好.但在大多数情况下, 长短窗可以很好地略过地震记录中的平静时间.同时, 长短窗具有递推的表达形式, 计算速度非常快.考虑到这一点, 本文使用长短窗算法作为触发条件.当长短窗被触发后, 使用支持向量机算法在其附近搜索存在的震相到时点.
本文还比较了本文方法和Cichowicz(1993)在拾取S波时的差异.运用本文的方法和Cichowicz的方法在S波理论到时前后10 s, 共20 s的时间窗内进行搜索.其具体结果如表 2.在S波的拾取上, 本文的方法明显地比Cichowicz的方法拾取的正确率更高、系统偏差更小、误差的方均根也更小.
|
|
表 2 本文中的方法与Cichowicz拾取S波震相方法的结果比较 Table 2 Comparison between the method proposed in this paper and Cichowicz′s S phase picking method |
综合上面的结果, 本文基于支持向量机的方法, 在拾取准确的P波和S波上面, 比传统的方法具有较大的优势, 特别是S波的正确拾取进一步提高了实用的可能性.当然, 本文的方法在准确确定起振点上面不够好, 但是可以结合AR-AIC方法进行改进.
4 基于台阵的震相识别和震相到时拾取策略在上面展示的单台拾取方法中, 本文使用了一个假设:在一定时间窗内, 到时最早的震相就是P波初至.在地震活动性比较弱的地方, 采用这种假设可以取得比较好的效果.但是在一些地震活动性比较强的地区, 比如余震捕捉AI大赛里面涉及的汶川地震后龙门山断裂带附近的地震, 这种假设就不能保证正确.所以, 单台方法不能处理不同地震的时空分布很接近的情况.另外, 在单台结果中也存在少数噪音被误识别为地震的情况.为此, 本文部分地借鉴了Chen和Holland(2016)的工作, 结合本文的两个支持向量机, 设计了一套基于台阵的拾取策略:首先将不同台站上对应于同一地震的拾取结果联系起来, 共同地确定相应台站上的应被拾取的P波、S波震相; 然后对地震进行定位; 接着根据定位结果计算每个台站的理论到时; 最后在理论到时附近重新拾取震相.
本文按照同一台站对同一地震记录的初至P波、S波的走时做了拟合, 它们近似地满足:
|
(14) |
其中, To表示发震时刻, Tp表示P波的绝对到时, Ts表示S波的绝对到时, 单位都是秒.
由(14)式, 可以利用P波和S波的绝对到时来推算发震时刻:
|
(15) |
同时, 根据发震时刻, 可以进一步计算出P波和S波的走时, 从而估计台站对应的震中距.利用震中距将地震可能的位置约束在以该台站为中心、震中距附近的一个环带内.不同台站得到的发震时刻越接近、环带区域共同重叠越多, 相应震相来自同一地震的概率越大.
因此, 我们将研究的区域划分为多个次区域, 根据单台结果中地震可能的震中位置, 将单台的震相对应到这些可能位置所在几个次区域.如果多个台站上, 震相对应的次区域相互重叠且发震时刻接近, 可以认为相应震相来源于同一地震.根据这种时空关系, 本文将不同台站对应于同一地震的拾取结果联系起来, 具体的做法如下:
(1) 先使用SSD, 在每个台站上拾取可能存在的震相到时点, 得到震相到时列表T.对T中的每一个震相到时Ti, 我们都假定其为可能的P波初至, 然后再运用SPS, 将Ti到时之后一定时间内的震相分为P波和S波, 选取出其中的S震相作为可能的S震相.
(2) 对于每一个可能的S震相到时Tj, 我们将其与Ti组成一个候选的P-S震相到时对(Ti, Tj).对每一组P-S震相, 都可以通过(15)式算出相应的发震时刻和P波走时, 并估计震中距.
(3) 本文将所研究的区域划分为若干个次区域并编号.每个台站都用一个二维数组Sec记录对发震时刻和震中位置的估算结果, 数组的一个元素对应于该台站在某一次区域内某一秒的估算结果:如果有震相所对应的发震时刻在第m秒、且根据震中距其震中可能位于第n个次区域, 则该台站的Sec数组的第m行第n列的元素被赋为1, 否则为0.考虑到本文的拟合结果存在残差, 也为了增加方法在其他地方的适用性, 发震时刻附近一定范围内Sec数组元素都会被赋值为1, P-S到时差越大, 这个范围也越大, 比如, 这里是按速度模型误差10%设计这个范围.
(4) 本文将多个台站的Sec数组中对应同一次区域、同一秒的元素叠加起来, 并取相应时间段、所有次区域内最高值所对应的时刻和次区域作为最可能的发震时刻和区域.根据最可能的发震时刻和区域, 在每个台站的拾取结果里选取最可能的P-S组合, 然后拾取P波、S波到时.这样就可以初步地选取出不同台站上对应于同一地震的震相.图 7是多台拾取结果叠加的示意图.
|
图 7 多台拾取结果叠加示意图 (a)网格展示每个台站的Sec数组对应同一次区域在该时间段的部分; 横坐标是时间, 纵坐标是台站序号; 当值为1时, 被填涂为蓝色, 否则为白色.(b)则表示的是将每一秒内所有台站对于该次区域的Sec值叠加的结果; 横坐标的时间与(a)的网格对应.(c)则展示了每一个次区域的叠加结果; 纵坐标是次区域的编号, 横坐标的时间与(a)(b)网格对应, 颜色表示叠加的结果. Fig. 7 Diagram showing multiple stacking results (a) The grid shows each station′s Sec array′s part of this time period in the same area; the abscissa is time, the ordinate is the station number; when the value is 1, it is filled with blue, otherwise it is white.(b) shows the result of stacking all the stations′ Sec arrays′ value in this sub-area; the time of the abscissa corresponds to the grid of (a).(c) shows the stacking result in each sub-area; the ordinate is the number of the sub-areas, the time of the abscissa corresponds to the (a) (b) grid, and the color indicates the value of stacking results. |
初步的结果中可能存在以下问题:其它地震的震相由于发震时刻的巧合而被错误地联系到该地震; 只拾取到P波而无法计算发震时刻等原因造成的有该地震的震相记录却没有被联系到该地震.对地震进行定位和再次拾取可以利用空间信息剔除错误震相、补全正确震相.具体做法如下:
使用程序velest(Kissling et al., 1994), 对找到的地震进行初步的定位; 在获得了地震的位置和发震时间的信息后, 使用seizmo工具包中基于taup(Crotwell et al., 1999)开发的MATLAB下的mattaup软件, 根据当地的速度结构(如Pasyanos et al., 2014), 计算地震在每个台站上的P波、S波的理论到时; 在理论到时附近的一个短的时窗内, 利用前述的两个支持向量机, 首先对初至P波进行拾取, 再对S波进行拾取.图 8展示了一个多台联合拾取的例子.
|
图 8 一个地震在台阵上的理论到时和自动拾取结果 按照震中距的远近依次绘出每个台站记录的垂直分量.蓝色和红色的圆圈分别表示P波和S波的理论到时, 蓝色和红色的竖线表示自动拾取程序拾取得到的P波、S波到时. Fig. 8 One quake′s theoretical arrival time and automatic pick result in the array The vertical component of each station record is drawn in turn according to the distance from the hypocenter.The blue and red circles respectively indicate the theoretical arrival time of P-wave and S-wave, and the blue and red vertical lines respectively indicate the P-wave and S-wave arrival times picked up by the automatic pick-up program. |
由于需要被至少三个台记录到, 在多台联合拾取策略下, 我们能够避免大多数的误识别情况; 结合了其他台站的信息计算得到理论到时, 当记录被其他时间相近的其他地震的震相干扰时, 也能够做出正确的拾取.
在拾取的过程中, 如果有两个或多个时空位置比较接近的地震, 本文的方法会倾向于识别被记录到的台站多的地震, 而忽略被记录到较少的地震.为了解决这个问题, 我们会再次进行前面的拾取过程.在重新拾取时, 前面所有的已经被拾取并联系到具体地震的震相不会再参与到拾取中来.
5 在实际数据中的应用余震捕捉AI大赛的主办方中国地震局地球物理研究所和阿里云计算有限公司于2017年10月21日公布了复赛阶段使用的数据的人工拾取结果2834个地震、20647条震相记录.本文将机器拾取的8135个地震、104090条震相记录与其比较.对每条人工拾取的结果都匹配上到时最接近的机器拾取结果, 并统计机器结果与人工拾取的误差大小.设立了不同的最大误差, 并统计了在该要求下的拾取率、平均值和误差均方根, 结果如表 3.本文将与人工拾取结果差异在2 s以内的记录作为拾取到了相应震相的记录结果.据此, 对满足这个最大误差要求的拾取结果统计了P波、S波的误差分布(见图 9).
|
|
表 3 自动拾取结果与人工拾取结果的比较 Table 3 Comparison between the automatic results and manual results |
|
图 9 自动拾取和人工拾取结果中P波(a)和S波(b)到时差分布直方图 横坐标表示自动拾取结果与手动拾取结果的差, 纵坐标为其对应的频率. Fig. 9 Picking difference distribution histograms for P-wave (a) and S-wave (b) arrival time between automatic and manual picking results The abscissa indicates the difference between the automatic picking result and the manual picking result, and the ordinate indicates the corresponding frequency. |
本文的方法中, 一个地震要被识别出来, 必须至少被三个台站记录到.人工拾取的结果中, 还有一些地震只被1~2个台记录到, 这部分是不能触发本文的方法进行拾取的.为了更加客观地反映本文方法在其设计范围内的识别能力, 剔除掉了被少于三台记录到的地震, 并做了同样的统计(见表 4).
|
|
表 4 自动拾取结果与经过筛选的人工结果比较 Table 4 Comparison between the automatic results and selected manual results |
本文还具体给出由程序自动拾取到的震相(图 10), 同时也给出相应的人工拾取结果.在实际的资料中, 地震记录可能被各种噪声干扰, 本文的方法能够很好地处理这种情况(图 10a).在有些情况下, 本文的自动拾取程序的结果比人工结果更加准确.如图 10b展示的例子中, 本文的自动拾取算法拾取的到时更加准确, 而人工拾取的P波到时前已经有一段明显的震相了.在汶川地震的余震活动非常活跃, 在一个小的时窗内可能同时记录到多个地震, 由于本文采用了多台联合拾取的策略, 能够处理这种情况, 图 10c就是其中的一个例子:在展示的时窗内有两个地震, 本文正确地将后一个P波匹配到该地震上, 和人工拾取的结果一致; 前一个P波也被正确拾取并匹配到另一个地震上.
|
图 10 自动拾取结果的四个例子 三道波形分别是地震信号的三分量.E分量上绘出了SSD经过sigmoid函数映射后的结果; N分量上绘出了区分SPS的结果; Z分量上绘出SSD结果的转换值; 转换值具有和SPS结果相同的正负号.自动拾取的结果(pTime和sTime)用竖直实线给出; 人工拾取结果(pTime0和sTime0)用带端点的竖虚线给出. Fig. 10 Four examples of automatic picking results The three waveform traces are the three components of the seismic signal respectively.We plot SSD′s sigmoid function mapping results at E component, the SPS′s results at N component and SSD′s converted value at Z component; the converted value has the same sign with SPS′s result.The auto-picked results (pTime and sTime) are given by the vertical solid line; the manual picking results (pTime0 and sTime0) are given by the vertical dotted line with end points. |
本文的自动拾取程序识别到了人工拾取结果约4倍的地震和约5倍的震相记录.由于评价一个自动拾取算法的能力不仅仅在于正确地拾取到多少的真正的震相, 同时也在于误拾取到了多少假的震相.我们检查了比人工拾取结果多拾取出来的部分(如图 10d所示), 发现基本上都是地震信号.因此, 本文的方法可以在保证识别正确性的前提下, 提供比人工拾取更完整的地震目录和到时数据.
6 结论与讨论本文利用支持向量机, 结合地震波特点及人工震相识别经验, 设计了两个震相分类器SSD和SPS.结果对比表明, 这两个分类器在震相识别准确度方面明显优于传统方法.同时, 本文基于支持向量机设计了一套利用单台资料并结合传统方法精确拾取体波震相到时的方法, 兼顾了震相识别的准确度和震相到时提取的精确度.
进一步地, 本文开发了一套可以用于台阵数据的P波、S波震相识别及震相到时拾取方法.利用这套方法, 可以在没有地震目录的情况下, 从多台连续记录中, 识别地震并对其P、S震相到时进行精确拾取, 可以为地震学研究提供可靠的地震目录和到时数据.
当然, 深度学习在使用了好的网络结构和大的训练样本的情况下可以不太依赖专业知识.即使没有专业经验, 也有潜能取得不错的识别效果.而支持向量机方法需要在一定程度上依赖于使用者的专业知识和经验.不过, 对于地震专业人员来说, 没有地震专业知识方面的障碍.一但选取了好的特征向量, 支持向量机就能够在较小的样本下获得很高的准确率和较强的泛化能力, 并降低人工标注的成本, 在数据量不是绝对大的情况下有明显优势.同时, 处理海量的地震资料的计算成本不可忽视, 支持向量机简单的数学模型和小的计算成本, 可以以较小的计算代价处理大量的数据.和传统方法结合使用, 在计算效率方面有优势.
对于地震活动性较强的区域, 特别是大震后很近的时空范围内, 大量地震频繁发生.这时, 台阵策略尤其必要.当然必须指出的是, 我们的台阵策略还是基于传统的定位加走时的时空约束.如果在台阵策略上也引入机器学习的方法, 研究同一地震不同台站间波形的关联与差异, 可以更好地检测地震和拾取震相到时.
我们将开发的程序代码分享到了我们GitHub的主页下, 网址是:https://github.com/baogegeJiang/jPicker, 欢迎大家下载使用.
致谢 中国地震局地球物理研究所提供了地震资料, 中国地震台网中心赵永高级工程师和苏金蓉研究员提供了人工震相识别经验, 感谢匿名审稿人的建设性意见, 感谢所有为本文做出贡献的人们！
Allen R V. 1978. Automatic earthquake recognition and timing from single traces. Bulletin of the Seismological Society of America, 68(5): 1521-1532. |
Bahlmann C, Haasdonk B, Burkhardt H. 2002. Online handwriting recognition with support vector machines-a kernel approach.//Proceedings Eighth International Workshop on Frontiers in Handwriting Recognition. Niagara on the Lake, Ontario, Canada: IEEE, 49-54.
|
Bai C Y, Kennett B L. 2000. Automatic phase-detection and identification by full use of a single three-component broadband seismogram. Bulletin of the Seismological Society of America, 90(1): 187-198. DOI:10.1785/0119990070 |
Batta P, Singh M, Li Z D, et al. 2018. Evaluation of support vector machine kernels for detecting network anomalies.//2018 IEEE International Symposium on Circuits and Systems. Florence, Italy: IEEE, 1-4.
|
Chang X, Liu Y K. 1998. Distinguishing seismic first break by means of hausdorff fractal dimension. Chinese Journal of Geophysics (in Chinese), 41(6): 826-832. |
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 |
Cichowicz A. 1993. An automatic S-phase picker. Bulletin of the Seismological Society of America, 83(1): 180-189. |
Cristianini N, Shawe-Taylor J. 2005. An Introduction to Support Vector Machines and other Kernel-Based Learning Methods. Beijing: China Machine Press: 93-124.
|
Crotwell H P, Owens T J, Ritsema J. 1999. The TauP toolkit:Flexible seismic travel-time and ray-path utilities. Seismological Research Letters, 70(2): 154-160. DOI:10.1785/gssrl.70.2.154 |
Ding S F, Qi B J, Tang H Y. 2011. An overview on theory and algorithm of support vector machines. Journal of University of Electronic Science and Technology of China (in Chinese), 40(1): 2-10. |
Dziewonski A M, Anderson D L. 1981. Preliminary reference Earth model. Physics of the Earth and Planetary Interiors, 25(4): 297-356. DOI:10.1016/0031-9201(81)90046-7 |
Fu J D, Chen L, Kang S H, et al. 2016. Transformer fault diagnosis based on dragonfly optimization algorithm and support vector machine. Journal of East China Jiaotong University (in Chinese), 33(4): 103-112. |
Gibbs M N, Mackay D J C. 2000. Variational Gaussian process classifiers. IEEE Transactions on Neural Networks, 11(6): 1458-1464. DOI:10.1109/72.883477 |
Jurkevics A. 1988. Polarization analysis of three-component array data. Bulletin of the Seismological Society of America, 78(5): 1725-1743. |
Kissling E, Ellsworth W L, Eberhart-Phillips D, et al. 1994. Initial reference models in local earthquake tomography. Journal of Geophysical Research:Solid Earth, 99(B10): 19635-19646. DOI:10.1029/93JB03138 |
Leonard M, Kennett B L N. 1999. Multi-component autoregressive techniques for the analysis of seismograms. Physics of the Earth and Planetary Interiors, 113(1-4): 247-263. DOI:10.1016/S0031-9201(99)00054-0 |
Li Q Z. 1996. How to correctly use fractal and fractal-dimension techniques?. Oil Geophysical Prospecting (in Chinese), 31(1): 136-163. |
Liu X Q, Zhou Y W, Qu J H, et al. 2009. Real-time detection of regional events and automatic P-phase identification from the vertical component of a single station record. Acta Seismologica Sinica (in Chinese), 31(3): 260-271. |
Pasyanos M E, Masters T G, Laske G, et al. 2014. LITHO1.0:An updated crust and lithospheric model of the Earth. Journal of Geophysical Research:Solid Earth, 119(3): 2153-2173. DOI:10.1002/2013JB010626 |
Peng Z G, Zhao P. 2009. Migration of early aftershocks following the 2004 Parkfield earthquake. Nature Geoscience, 2(12): 877-881. DOI:10.1038/ngeo697 |
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 |
Platt J C. 1999. Sequential minimal optimization: A fast algorithm for training support vector machines.//Advances in Kernel Methods-Support Vector Learning. Cambridge: MIT Press, 185-208.
|
Saragiotis C D, Hadjileontiadis L J, Rekanos I T, et al. 2004. Automatic P phase picking using maximum kurtosis and/spl kappa/-statistics criteria. IEEE Geoscience and Remote Sensing Letters, 1(3): 147-151. DOI:10.1109/LGRS.2004.828915 |
Sleeman R, Van Eck T. 1999. 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. DOI:10.1016/S0031-9201(99)00007-2 |
Wang C X. 2014. Automatic onset time picking algorithms for multiple seismic arrivals[Master's thesis] (in Chinese). Xi'an: Chang'an University.
|
Wang J, Teng T L. 1997. Identification and picking of S phase using an artificial neural network. Bulletin of the Seismological Society of America, 87(5): 1140-1149. |
Wang W J, Meng X F, Peng Z G, et al. 2015. Increasing background seismicity and dynamic triggering behaviors with nearby mining activities around Fangshan Pluton in Beijing, China. Journal of Geophysical Research:Solid Earth, 120(8): 5624-5638. DOI:10.1002/2015JB012235 |
Wu M Y, Mao S J, Ning J Y, et al. 2018. Fault structure detection by matched filter technique. Acta Scientiarum Naturalium Universitatis Pekinensis (in Chinese), 54(4): 730-738. |
Yin D Y. 2012. Automatic P-wave arrival time detection and real-time magnitude estimation for earthquake early warning[Master's thesis] (in Chinese). Harbin: Institute of Engineering Mechanics.
|
Zeng F Y, Li M F, Shen W. 2002. The fractal study on detecting arrival time of the first break. Geoscience (in Chinese), 16(2): 209-213. |
Zhang S Y, Wu K Y, Huang Y Z, et al. 2017. Fall detection algorithm based on SVM_KNN. Computer and Modernization (in Chinese), (12): 49-55. |
Zhu L P, Rivera L. 2002. A note on the dynamic and static displacements from a point source in multilayered media. Geophysical Journal International, 148(3): 619-627. DOI:10.1046/j.1365-246X.2002.01610.x |
常旭, 刘伊克. 1998. Hausdorff分数维识别地震道初至走时. 地球物理学报, 41(6): 826-832. DOI:10.3321/j.issn:0001-5733.1998.06.013 |
丁世飞, 齐丙娟, 谭红艳. 2011. 支持向量机理论与算法研究综述. 电子科技大学学报, 40(1): 2-10. |
傅军栋, 陈俐, 康水华, 等. 2016. 基于蜻蜓算法和支持向量机的变压器故障诊断. 华东交通大学学报, 33(4): 103-112. DOI:10.3969/j.issn.1005-0523.2016.04.017 |
李庆忠. 1996. 怎样正确对待分形、分维技术?. 石油地球物理勘探, 31(1): 136-163. |
刘希强, 周彦文, 曲均浩, 等. 2009. 应用单台垂向记录进行区域地震事件实时检测和直达P波初动自动识别. 地震学报, 31(3): 260-271. DOI:10.3321/j.issn:0253-3782.2009.03.003 |
王彩霞. 2014.多震相初至自动检测识别方法技术[硕士论文].西安: 长安大学.
|
吴梦羽, 毛淑娟, 宁杰远, 等. 2018. 基于模板识别方法探测断层结构. 北京大学学报(自然科学版), 54(4): 730-738. |
尹得余. 2012. P波震相自动捡拾与地震预警震级实时测定技术研究[硕士论文].哈尔滨: 中国地震局工程力学研究所.
|
曾富英, 李敏锋, 申维. 2002. 地震波初至拾取的分形研究. 现代地质, 16(2): 209-213. DOI:10.3969/j.issn.1000-8527.2002.02.016 |
张舒雅, 吴科艳, 黄炎子, 等. 2017. 基于SVM_KNN的老人跌倒检测算法. 计算机与现代化, (12): 49-55. DOI:10.3969/j.issn.1006-2475.2017.12.010 |