地球物理学进展  2017, Vol. 32 Issue (3): 1000-1007   PDF    
地震自动识别及震相自动拾取方法研究进展
刘翰林, 吴庆举    
中国地震局地球物理研究所, 北京 100081
摘要:微小地震事件的研究工作应用领域广泛,在诱发地震监测、地颤动低频地震监测、地震预警等研究中有举足轻重的地位.而微震事件的自动识别以及震相的自动拾取的研究则是后期微震数据处理工作的基础,包括地震定位、地震成像、震源机制反演等工作.本文就目前各类地震事件自动识别及震相拾取方法的原理起源与大致分类应用进行了阐述,重点介绍了三大类方法:第一类是基于振幅(能量)判据的长短时窗比、震源扫描叠加等事件识别以及震相拾取系列方法;第二类是基于波形互相关技术的事件识别方法;第三类是基于波形自相关的盲搜索技术.最后还列举了几类较为前沿的综合性方法及其应用分析,并与传统搜索和拾取方法进行了比对,认为各类单一传统方法各有利弊,不具有普适性,应就具体的研究领域进行综合应用.
关键词微震    地震事件自动识别    震相自动拾取    长短时窗比    震源扫描叠加    波形互相关    波形自相关    
Developments of research on earthquake detection and seismic phases picking
LIU Han-lin, WU Qing-ju    
Institute of Geophysics, China Administration of Earthquake, Beijing 100081, China
Abstract: The research on microseism has a wide application area, for instance, the monitoring of induced earthquake and low frequency earthquakes in earth tremor, the earthquake early warning. And the automatic identification of events or variable phases is the fundamental of later data processing, comprising seismic location and imaging, inversion of focal mechanism, etc. Some basic theories, also the typical applications of various approaches to earthquake detection and seismic phases picking are discussed in this paper, especially focusing on the three categories, i.e., the energy detector basing on the amplitude, including Short-Term Average (STA)/Long-Term Average (LTA) and Source Scanning Algorithm (SSA) picker, the template matching basing on waveform cross-correlation and the blind searching method basing on waveform auto-correlation. We also list some of the edge techniques of earthquake detection and seismic phases picking, and compare them with traditional methods. Ultimately, we suggest to apply a multiple algorithm to a specific research instead of a unitary method, which is not universal.
Key words: microseism     earthquake detection     seismic phases picking     STA/LTA     source scanning algorithm     waveform cross-correlation     waveform auto-correlation    
0 引言

最早期的地震学工作中,台网工作人员都是通过肉眼人工来识别地震事件,并结合走时表计算拾取各类震相,但随着日益密集的地震台网分布,地震数据的日益庞大,堪称海量,而地震记录中所包含的内容也越来越丰富复杂,很多微小地震淹没在背景噪声中,人工识别显然已不能满足所需的工作效率及研究精度,因此对于无论是流动台网还是固定台网的地震数据中,关于地震事件的自动识别以及震相的自动拾取研究方法的发展越来越成熟,应用越来越广泛.并且对基于地震速报、地震预警的防震减灾工作,地震事件的自动识别及震相自动拾取工作也是其重要基础(武东坡, 2004).Stevenson(1976)提出了应用长短时窗比值法来判定地震事件初至到时,其后发展了各种基于该方法下的从时间域、频率域、偏振分析、高阶统计量、模式识别等各方面综合分析的震相拾取方法(武东坡, 2004).Kao和Shan(2004)提出震源扫描叠加(SSA)地震定位法后,J.Zahradnik等(2015)在其基础上提出了另一种基于振幅判据的p波震相识别方法.Schaff和Richards(2004)提出了基于波形互相关的重复地震的概念,随后衍生了一种利用台站记录与模板波形波形进行互相关地震事件检测(模板匹配技术)的方法(Gibbons and Ringdal, 2006; Gibbons et al., 2007; Shelly et al., 2007; Peng and Zhao, 2009); Zhang和Wen(2014)接着提出了基于模板匹配技术的M&L(Match&Locate)事件识别定位方法.为了将检测方法推向一般化,摆脱对模板的依赖,Brown等(2008)提出了一种基于波形自相关的彻底的盲搜索技术,随后Clara E. Yoon等(2015)对其算法进行了改进,提出FAST方法,实现了波形自相关盲搜索技术的计算可行性.

1 能量法 1.1 长短时窗比

该方法最初由Stevenson(1976)提出来进行地震事件识别和初至到时拾取.STA/LTA(Short Term Average/Long Term Average)即长短时窗比法的原理中,STA和LTA分别指的是信号在固定长度的短时窗和长时窗中的平均值,然后用前者与后者之比值来描述该段信号的振幅即能量的变化情况.当信号遇到突变时,信号的短时窗平均值变化会比信号的长时窗平均值快,从而比值增加,若设定一个合适的阈值,当比值超过阈值时,则判定地震事件信号的到来.长短时窗及比值定义式为

(1)
(2)
(3)

式中,NM分别为短时窗和长时窗中的采样点数,Y(n)为地震信号,R为长短时窗比值.

在最初的方法中是直接用地震序列值作为输入进行长短时窗平均值计算,该输入值对于地震信号的到达不够灵敏,后来Allen(1978)提出使用特征函数代替Y(n)作为输入来提高算法对信号突变的灵敏度.随之Baer和Kradolfer(1987)又对其进行了改进,提出用一个包络函数代替特征函数并将原始的固定阈值动态化,提高了精度.意识到用以计算长短时窗平均值的序列特征函数是影响算法精度的关键后,Earle和Sheare(1994)高淑芳等(2008)余建华等(2011)刘晗和张建中(2014)等分别提出了不同的特征函数以测试其灵敏度和拾取精度.目前最常用的几种特征函数为

(4)
(5)
(6)
(7)
(8)

式中CF(i)即为特征函数,Y(i)为地震信号.

1.2 基于长短时窗比的改进方法

由于考虑到地震记录中的随机噪声及一些高频脉冲对信号检测的干扰,长短时窗比只能大致确定时间信号到时的位置,并且有时候会产生非地震信号拾取的误判,学者们又纷纷提出基于长短时窗比的各类提高精度的改进方法,在利用长短时窗比大致判定拾取事件到时后,再应用其他方法进行进一步精确拾取到时.

1.2.1 赤池准则法(AIC)

赤池信息量准则(AIC:Akaike Information Criterion)首次由日本统计学家赤池弘次提出,用来描述统计模型与实际数据拟合的好坏的一个标准.该准则建立在熵的基础之上,赤池函数一般可表示为

(9)

式中k为模型独立参数个数,L为似然函数,当AIC值取得最小时,模型最佳.将赤池准则与地震信号判别相结合时,我们认为信号突变的前后是两个不同的趋于稳态的模型,若将n看作两个模型的分界点,然后计算两个模型的AIC函数值,取最小时,模型最佳,也意味着此时的n为两个稳态的分界点(Sleeman and van Eck 1999).表达式为

(10)

式中M为自回归阶数,N为信号长度,k为分界点,C为常数,σ12σ22分别表示前后两个模型的方差.

考虑到自回归模型与赤池准则的结合时,AR(自回归)阶数过大会导致整个算法计算量过大,于是有学者引入数学中的高阶统计量与赤池准则进行结合判定.数学期望(mathematical expectation)即均值为一阶统计量,二、三、四阶统计量分别为方差(variance)、偏斜度(skewness)、峰度(kurtosis),表达式分别为

(11)
(12)
(13)

式(11) 表示随机变量Xk阶统计量,描述X与数学期望的偏离程度;式(12) 表示随机变量的偏斜度,描述随机变量相对于均值的不对称程度;式(13) 表示随机变量的峰度,描述随机变量分布曲线的尖峭程度(赵大鹏等, 2013).Maeda(1985)首次将方差引入与赤池准则进行结合,排除了AR-AIC中的高阶矩阵求逆复杂运算,提高了效率,表达式为

(14)

式中,k为分界点,N为信号长度,分别为分界点前后两个阶段的方差,仍以赤池函数去最小值时的k作为信号到来的判别.王继等(2006)VAR-AIC(方差-赤池)方法与AR-AIC(自回归-赤池)方法进行了对比应用,取得了较好的结果;马强等(2013)牟磊育和吴建平(2015)也分别对AR-AIC方法进行了应用,并对结果进行了Delauney三角时间空间合理性约束,取得了较好的结果;刘希强等(2009)赵大鹏等(2013)分别将VAR-AIC中的方差替换成了偏斜度和峰度进行了对比研究,得出信噪比较高的情况下,方差、偏斜度、峰度三种方法结果大致一致,而信噪比较低的情况下,峰度方法更为精确;刘劲松等(2013)提出了基于长短时窗、AIC及PAI-S/K方法的结合,并应用了大量的微震数据验证与对比,对传统方法进行了有效的改进.

1.2.2 偏振分析法

由于p波尾波的干扰,长短时窗比方法对于s波的识别效果较差,在关于s波震相自动拾取方面,一般采用偏振分析方法(刘希强等, 2014).Roberts等(1989)首次据p、s波在传播过程中的质点振动方向差异,提出用计算信号中质点偏振方向将s波识别出来.该方法的基本思路为(Amoroso et al., 2011):

(1) 选取p波到时后的一定长度信号,计算地震记录三分量的协方差矩阵的最大特征向量,该方向即可估计为p波偏振方向,公式为

(15)

式中dec(t)为三分量的协方差矩阵,cov(n, n)、cov(n, e)、cov(n, z)等为地震记录的三分量两两之间的协方差(Montalbetti and Kanasewich, 1970).

(2) 以p波偏振方向为标准,将地震三分量旋转至LQT分量即垂向(p波振动方向)、径向(sv波振动方向)、切向(sh波振动方向)三个分量.公式为

(16)

式中,φ为射线入射角,β为直达p波反方位角(Plešinger et al., 1986).

(3) 计算地震信号的偏振函数.首先类似第一步选取一定长度时窗计算LQT三分量的协方差矩阵,进而计算得到三个特征值:λ1, λ2, λ3(从大到小),然后信号偏振函数表达式为

(17)

当偏振函数曲线取最大时,为s波信号到时.

Cichowicz(1993)提出联合射线偏斜角度、偏振函数、能量比三个方面的因素进行s波信号判别;Diehl等(2009)对该方法进行了应用并提出了利用拾取到的s波到时进行层析成像的方法;Ortensia Amoroso等(2011)在方法中又加入了加权因子,并对理论地震图与实际数据进行了应用对比,得到了较好的结果;刘希强等(2014)应用此方法对山东数字地震台网数据进行了处理并与人机交互方法进行对比得到了较好的结果.

1.2.3 频率域分析法

当信号信噪比较低时,振幅判据可能对于地震初动信号不太敏感了,但是由于地震信号相比噪声信号的频率特殊性,频谱判据就相对很敏感了,于是有学者提出相关的频率域识别法进行分析.

提到频谱分析,大家自然想到最早的傅里叶变换,但是傅里叶变换只是局限性地把信号分解成无限长的正、余弦波的叠加,显示的只是频率信息,而地震信号是一种频率成分随时间变化具有规律特征性的非稳态信号,因此需要对局部时域进行频谱分析,时频分析法很好地解决了该问题.

在用于震相拾取的研究中,小波变换法是目前国内外较为流行的一种时频分析方法.大致思路都是先应用长短时窗比进行粗略到时拾取,然后再应用小波变换法对信号分解重构以压制噪声,最后再通过赤池准则法或者偏振度方法等进行精确到时拾取.小波变换于1982年由勘探地球物理学家Morlet提出,继承了短时傅里叶变换(STFT)局部化的思想,同时又克服了窗口不随频率变化的缺点,即能提供随频率变化的时频窗口.小波函数的定义如下(武东坡, 2004):设ψ(t)为一平方可积函数,即ψ(t)∈L2(R),若其傅里叶变换满足条件:

(18)

则称为一个基本小波或小波母函数.小波函数族是由同一小波母函数通过伸缩和平移后所得到的,相对于用一系列正余弦函数去表示信号,小波变换就是采用一族小波函数去表示信号.通过小波变换伸缩平移变换对信号逐步进行多尺度细化,高频处时间细分,低频处频率细分,自适应时频分析的要求,对时间频率局部化分析,充分突出信号某些特定特征.该方法在震相识别中的应用一般思路是:利用小波变换系数中包含的信号特征,寻找不同尺度下小波变换系数的显著特征,通过分析小波变换系数的主成分以剔除随机背景噪声成分,得到震相的识别特征,从而进行震相拾取(刘希强等, 2000).Yomogida(1994)Chakraborty和Okaya(1995)Savvaidis等(2002)等人对该方法均有研究应用.

1.2.4 希尔伯特——黄变换法(HHT)

在处理地震记录等非平稳信号的方法中,还有一类叫做希尔伯特黄变换(HHT).1998年中国台湾海洋学家黄锷等人提出经验模态分解方法(EMD),并引入希尔伯特谱分析概念,得到HHT方法.该方法的核心理念是:对于任意一个非平稳信号,可根据信号的特点自适应地将其分解为互不相同的、简单非正弦的一列本征模态函数(IMF)组,并对每一个IMF作希尔伯特变换,继而可求得每一个IMF的瞬时频率.HHT方法在震相自动拾取研究中的应用(贾瑞生等, 2015)如下:将记录到的含噪地震信号进行经验模态分解,并根据背景噪声的特征将其IMF分量剔除,然后对信号进行重构降噪并对其作希尔伯特变换,求得包络信号,归一化处理后设置包络阈值即可开始震相拾取.

1.3 基于震源扫描叠加技术(SSA)的左右时窗比(RPA/ LPA)

震源扫描叠加技术(SSA)由Kao和Shan(2004)提出,是一种试验性搜索扫描地震定位方法,对可能的发震时刻和震源位置进行时间和空间上的网格化,并对每一空间格点计算其与台站的旅时,根据可能的发震时刻将各台站记录到的振幅绝对值进行叠加取平均,并命名其为“明亮度”函数,使得明亮度函数取最大值时的格点位置和格点时刻即为地震震源参数.公式为

(19)

式中br(η, τ)为“明亮度”函数,N为台站个数,n为台站序号,ητ分别为可能的震源位置和发震时刻,tηn为格点位置到第n个台站的旅时.

J.Zahradnik等(2015)基于SSA技术发展了一种新的p波震相自动拾取的方法. (1) 定义震相拾取函数——左右时窗比(RPA/LPA).由两个分别位于某一时间点左右侧的滑动平均时窗构成,左右时窗窗长相等.公式为

(20)

式中τj地震记录采样间隔,J为左右两个时窗窗长,s(t)为地震记录振幅.

(2) 叠加各台站的震相拾取函数曲线求平均.将各台站的震相拾取函数曲线按每一个可能的震源位置试验点与台站之间的旅时进行校正后叠加,并求取平均.公式为

(21)

式中m为可能的震源位置网格点序号,N为台站数量,ttnm为第n个台站与第m个网格点的旅时.

(3) 时间空间同时搜索.将可能的发震时刻亦网格化成oti,分别计算每个发震时刻试验点和震源位置试验点的“明亮度”函数,同时从时间和空间进行搜索,“明亮度”函数取最大值时即为发震时刻和震源位置,从而得到p波震相到时.公式为

(22)

该方法在计算试验点与台站之间的旅时时会比较依赖地下结构信息的精确度.

2 基于波形互相关的模板匹配技术

对于一些发生在同一断层位置的微小重复地震和一些强震过后的余震序列,它们的波形具有很好的相似性,于是便衍生出了一套基于波形互相关的识别微小地震事件的研究方法.

2.1 重复地震与余震序列

Schaff和Richards(2004)对重复地震给出实际操作层面上的定义:两次地震事件至少被同一个台站记录到的p波前5 s到Lg波后40 s的时窗内波形之间互相关系数不小于0.8,则两次事件称为重复地震对.在此之前,关于重复地震是没有一个定性的统一定义的,一般认为发生在同一断层位置上,波形相似性好,震级接近的两次或两次以上的地震事件称为重复地震,它们通常具有十分接近的复发间隔以及相同的震源机制(Poupinet et al., 1984; Rubin, 2002).而余震序列(或地震序列)指的是在一次强震发生的前后,同段时间内(几天、几个月、或几年),同一震源区发生的一系列大大小小的地震,它们有共同的发震构造和相互关联的震源机制,同样的,它们的波形也具有很好的相似性.但是严格地从物理角度来说,重复地震和余震序列是有本质区别的:重复地震发生在同一断层位置,破裂区域重叠,余震序列是在断层的不同位置,破裂区域独立分开,它们在空间位置上是存在差异的(Rubin and Gillard, 2000).

这些重复地震或者余震序列的微小地震事件,对于一些高精度段层面的确定、前震活动性、地震触发、地颤动低频地震检测、小当量核试验检测等工作是十分关键的,因此对于它们的研究是具有重要意义的.

2.2 常规波形互相关叠加法

关于两个信号之间的互相关系数,定义式为

(23)

式中x(t)、y(t)分别为两个信号,分别为两个信号的算数平均值.可以看出该公式的意义为:互相关系数按积差方法计算,两个信号去均值后,通过其分别的离差相乘来反映两个信号的相关程度;亦可理解为通过两个信号的离差组成的n维向量在Rn空间中所成夹角的余弦来表示两个信号的相关程度,具有尺度不变性.当两个信号具有很好的时间一致性时,相关系数会较大,反之则否.在地震记录中两个来自重复地震记录的信号是具有很好时间一致性的,而背景噪声是随机杂乱无章的,因此可通过互相关运算来压制噪声干扰,拾取有效的地震信号.

实际操作互相关运算进行地震信号拾取步骤为:

(1) 确定地震模板.从台网目录中获取地震事件,以各台站记录到的震相清晰、信噪比高的事件信号y(t)为模板.

(2) 滑动互相关运算.获取台站连续波形记录,将其作为待识别信号x(t),与该台站的模板信号y(t)进行滑动互相关运算,公式为

(24)

假设连续信号x(t)和模板信号y(t)分别有NM个采样点,而N»M,从x(t)的起点开始,截取其与y(t)相同长度M的一段,与y(t)进行如上式的逐点滑动互相关运算,得到互相关系数序列C(i).式中i=1, 2,…,N-M+1; n=1, 2, …, M.

(3) 互相关波形叠加判别地震事件.将各台站三分量的互相关系数序列波形进行叠加求平均,设定合理的判别阈值,当互相关系数超过阈值即判定为地震事件信号.

波形匹配技术能有效地压制噪声干扰,并且Schaff(2010)通过对1999年岫岩MS 5.4地震序列的研究,指出波形匹配技术相比于长短时窗比法将完备震级提高了1.3个震级单位.目前该方法应用十分广泛,Gibbons和Ringdal(2010)将其应用于2009年5月25日的核爆检测,Shelly等(2007)将波形匹配技术应用到地颤动中的低频地震检测中;Peng和Zhao(2009)用该方法进行了余震检测;吕鹏等(2011)将其与双差定位结合应用在2008年汶川地震余震序列搜索工作中;郑晨等(2015)在灌县-安县断裂的精细结构研究中应用了该方法;谭毅培等(2015)将该方法应用到2015年河北滦县震群中遗漏地震的检测及精定位,均获得了较好的结果.

3 基于波形自相关的盲搜索技术

模板匹配技术中,先验的模板信息的质量是关键因素,而模板的选取通常是通过台网目录下载和人工拾取,因此该类方法受模板事件约束,不具有很好的普适性.尽管后来的Harris(2006)以及S. A. Barrett和G. C. Beroza(2014)分别提出的子空间检测和经验子空间检测方法,即将匹配检测的模板波形多样化,使其更具有一般性,但并没有从根本上解决检测方法的普适性问题.

波形自相关地震检测技术,是完全独立于检测模板的系统性的“多对多”匹配搜索方法,是不需要先验的地震信号信息,彻底的盲搜索技术.该方法由J.R.Brown等(2008)在解决地颤动中的低频地震搜索问题时提出.方法思路为:

(1) 将各台站各分量的连续时间序列分别按一定的时窗长度和时窗间隔分成若干个时窗.

(2) 将所有可能成对的时窗进行两两互相关,并将对应该两个时窗的所有台站的各分量的互相关波形进行叠加,当互相关系数超过一定阈值时,它们被标记为候选事件,表达式为

(25)

式中n为台站数量,ij分别为两个时窗,u为时窗内的地震位移记录;

(3) 将所有的候选事件对中对应在同一台站的波形进行额外的互相关运算,并按互相关系数高低进行分组并叠加形成高信噪比模板,然后运用模板匹配滤波技术进行补充搜索.

波形自相关方法很好地避免了对模板先验信息的依赖性,实现了方法一般化推广,同时也继承了互相关运算对噪声压制的优点,但是该方法的致命弱点是计算量要求过高,对于大量的连续数据搜索运算来说,该方法是不实用的.

4 几种新的综合性方法

近年来,各学者为了弥补长短时窗、偏振分析、赤池准则及波形相关等各基础类方法的不足,提出了几种新的综合性方法.

4.1 SLPEA拾取算法

由谭玉阳等(谭玉阳等,2016)提出的SLPEA初至拾取新算法,是基于STA/LTA、偏振分析以及引入的边缘检测函数而设计的.该方法同样以滑动时窗搜索的思路,但是将目标函数改为了归一化处理后的长短时窗比值与偏振度的乘积,表达式为

(26)
(27)
(28)

上述三个式子中,前两个分别为长短时窗比函数与偏振度函数的归一化过程,第三个式子即为具有增强的抗噪性的敏感目标函数.结果表明,检测函数Ki的极值点通常滞后初至震相,而震相却更接近离Ki极值点前的斜率最大点.于是研究者又引入了以下边缘检测函数Di,公式为

(29)

式子中N为滑动窗口长度,Ki则是前者的综合检测函数,Di检测的是Ki变化的快慢程度,在Ki变化最快的位置(即之前提到的斜率最大点)存在一个局部极大值,这样Di的极值点比Ki极值点更靠近初至震相,则可选取Di最大的两个局部极值点作为P、S波初至到时初步估计值.然后以两个初步估计点为中心延拓一个短时窗,计算窗内信噪比S,设定信噪比阈值C1,并对短时窗再应用边缘检测函数Di得到对应的检测到时TD,计算短时窗内地震记录的最大似然函数Lk最大值点对应得到检测到时TL,将LkDi作归一化处理后相乘得到检测函数Fk,对应检测到时为TF,如果最初的信噪比S超过阈值,则选取TL作为最后结果,否之再设定误差阈值C2,如果TDTL之间的差值小于C2,则选取TDTL的中点作为最后结果,否之则选取TF作为最后结果.

研究者将该方法分别对模型数据和实际数据进行了应用验证,并与传统方法和人工拾取的结果进行了对比,平均绝对误差小于3个采样点,证明了该方法的可靠性.

4.2 基于能量变化率的二次方自回归模型拾取

何先龙等(2016)提出了一种基于信号能量变化率的长短时窗比与二次方自回归模型综合拾取的新方法.不同于传统方法的直接对信号的振幅序列进行处理,该方法先将三分量信号序列转化成能量变化率序列,公式为

(30)

xi为地震信号序列,dt为采样间隔,Ei即为瞬态震动能量变化率,对Ei应用长短时窗比进行初步拾取,得到一个到时范围区间,然后对区间内序列ei基于二次方自回归模型法精确拾取.设定理想的二次方自回归模型为

(31)

并用以下序列对模型进行逼近,公式为

(32)
(33)
(34)
(35)
(36)

i为序列序数,2s为区间序列长度.对于式(36) 通过最小二乘求解系数ωgf,代入表达式得到序列ci最大值对应的点即为精准到时.

笔者同样通过实际数据对该方法进行了验证,与人工拾取及传统拾取方法结果进行了对比,说明了方法的可靠性及准确性.

4.3 波形互相关叠加——空间扫描定位法(M&L)

在检测到地震信号之后的例如地震定位、地震成像的工作中,一般采用基于绝对走时计算或震相到时精确拾取的常规方法,如双差定位、层析成像等等,这些工作都依赖绝对到时的拾取和地下速度模型的精确度.Zhang和Wen(2015)提出了一种基于波形匹配技术的能同时进行搜索和定位地震的方法.具体思路为:

(1) 与常规波形匹配技术一样,通过模板事件对各台站分量进行滑动互相关搜索,得到互相关波形序列.

(2) 计算走时差并对各台站互相关波形进行走时校正后叠加.围绕着模板事件所在位置,将可能发生地震事件的区域沿着经度、纬度、深度三个方向进行网格化,认为每个格点都是潜在的震源位置;然后运用双差分思想,计算对于每个台站的可能事件位置与模板事件位置在该台站产生的走时差;最后将每个台站的互相关波形序列按每个可能的发震时刻和网格点位置计算得到的走时差进行校正后叠加求平均.

(3) 联合平均互相关系数和信噪比设定阈值进行判别.综合考虑平均平均相关系数和左右信噪比因素设定阈值,当超过阈值时,认为检测到一个地震信号,并把互相关系数最大时所对应的网格点位置认为是地震事件的位置.

该方法避免了依赖于地下速度模型的绝对走时的计算,走时差的计算是用到了双差思想,也避免了绝对震相到时拾取不够精确的问题,采用了互相关约束的搜索定位方法,并且事件识别和定位同时进行,大大提高了工作效率.Zhang and Wen(Zhang and Wen, 2015)该方法应用到2010年5月12日的朝鲜核试验地震学证据检测中,取得了较好的结果.

4.4 波形自相关地震搜索方法的算法改进——FAST (Fingerprint And Similarity Thresholding)

LSH(Locality-sensitive hashing/局部敏感哈希算法)是数据挖掘算法中的一种近似最近邻搜索算法,具有严密坚实理论依据且在高维数据中表现优异,主要作用是从海量的数据中挖掘出相似的数据,最早设计用来在大数据库中识别相似音频片段,大体思路是如果两个信号在原有的数据空间是相似的,那么在分别经过哈希函数转换以后仍能保持很高的相似性,反之则否(Andoni and Indyk, 2006).我们可以通过比较从原信号提取出来的可区别的关键特征——“指纹”,来替代直接的原信号之间比较,并且该算法可以避免比较不相似的“指纹”对,从而大大降低了计算量,提升了搜索效率.Clara E. Yoon等(2015)结合LSH算法提出了一种基于波形自相关的地震搜索的新方法——FAST(Fingerprint And Similarity Thresholding),很好地解决了传统波形自相关地震搜索方法的庞大计算量的瓶颈.首先我们将原始连续地震记录波形进行频谱分析,提取可区别的关键性特征,并对其进行压缩生成二进制“指纹”,然后运用LSH算法对数据库中的二进制“指纹”进行一个快速的相似性搜索成对,最后生成地震事件检测列表.该方法对6个月的连续数据进行搜索只需不到三天的运行时间,而对与传统的自相关搜索方法来说,则需要至少20年的运行时间,可见FAST方法在计算中的强大时效性.但不足的是计算过程中需要较大的内存支撑大量的特征“指纹”的存储,并行运算也是很必要的.

5 结论和讨论 5.1

总的来说,对于基于长短时窗比(sta/lta)的一类方法,它们不依赖于模板事件的选取,可操作性高,计算快速简单,但是又太过于依赖地震记录的信噪比高低,信噪比较高时效果可能较好,但当信噪比很低时,阈值设定过大会导致漏检率增高,很多小震级地震,会检测不到阈值设定过小又会导致误检率增高,很多为无效噪声干扰,如果要对记录信号进行频谱分析然后通过不同频段滤波来压制背景噪声,由于噪声的时间空间随机性,在做时频分析时窗长、窗函数等参数的选择不具有很好的一致性,导致工作量过大.而对于基于波形互相关的模板匹配技术一类方法,它们压制背景噪声的能力强,不太依赖记录的信噪比,检测精度和可信度很高,但是却过于依赖模板事件的选择,对于有些小震多发区可能地震序列和重复地震较少,而有些地区地下断层构造可能过于复杂导致地震事件的多样性,从而模板选取会很困难.基于波形自相关的搜索方法,虽然在一定程度上解决了前两者的短处问题,不仅能通过相关计算很好压制噪声,也能避免模板依赖造成的方法普适性不足,但是计算量和计算所需内存则是最大的瓶颈了.

5.2

鉴于各类单一传统方法各有优弊,我们应可各取所需,综合应用,例如在核爆检测方面的研究工作中,大多数情况下模板先验信息较全,并对事件定位要求较高,则可采用M&L来同时进行精确检测和定位,而在地颤动中的低频地震检测研究工作中,对于被检测事件的多样性要求较高,可采用FAST方法进行更完备普适地检测地震.而对于已检测到的微震事件,我们可以采用类似于SLPEA等综合性算法进行震相精确拾取.

致谢 感谢国家自然科学基金41274088和41474074对本研究的联合资助.
参考文献
[] Allen R V. 1978. Automatic earthquake recognition and timing from single traces[J]. Bulletin of the Seismological Society of America, 68(5): 1521–1532.
[] Amoroso O, Maercklin N, Zollo A. 2011. S-wave identification by polarization filtering and waveform coherence analysis[C].//EGU General Assembly, European Geosciences Union.
[] Andoni A, Indyk P. 2006. Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions[C].//47th Annual IEEE Symposium on Foundations of Computer Science. Berkeley, CA: IEEE.
[] Baer M, Kradolfer U. 1987. An automatic phase picker for local and teleseismic events[J]. Bulletin of the Seismological Society of America, 77(4): 1437–1445.
[] Barrett S A, Beroza G C. 2014. An empirical approach to subspace detection[J]. Seismological Research Letters, 85(3): 594–600. DOI:10.1785/0220130152
[] Brown J R, Beroza G C, Shelly D R. 2008. An autocorrelation method to detect low frequency earthquakes within tremor[J]. Geophysical Research Letters, 35(16): L16305. DOI:10.1029/2008GL034560
[] Chakraborty A, Okaya D. 1995. Frequency-time decomposition of seismic data using wavelet-based methods[J]. Geophysics, 60(6): 1906–1916. DOI:10.1190/1.1443922
[] Cichowicz A. 1993. An automatic S-phase picker[J]. Bulletin of the Seismological Society of America, 83(1): 180–189.
[] Diehl T, Deichmann N, Kissling E, et al. 2009. Automatic S-wave picker for local earthquake tomography[J]. Bulletin of the Seismological Society of America, 99(3): 1906–1920. DOI:10.1785/0120080019
[] Earle P S, Shearer P M. 1994. Characterization of global seismograms using an automatic-picking algorithm[J]. Bulletin of the Seismological Society of America, 84(2): 366–376.
[] Gibbons S J, Ringdal F. 2006. The detection of low magnitude seismic events using array-based waveform correlation[J]. Geophysical Journal International, 165(1): 149–166. DOI:10.1111/gji.2006.165.issue-1
[] Gibbons S J, Ringdal F. 2010. Detecting the DPRK nuclear test explosion on 25 May 2009 using array-based waveform correlation[C].//EGU General Assembly Conference. Vienna, Austria: EGU.
[] Gibbons S J, Sørensen M B, Harris D B, et al. 2007. The detection and location of low magnitude earthquakes in northern Norway using multi-channel waveform correlation at regional distances[J]. Physics of the Earth and Planetary Interiors, 160(3-4): 285–309. DOI:10.1016/j.pepi.2006.11.008
[] Harris D B. 2006. Subspace detectors: Theory[R]. Livermore, CA: Lawrence Livermore National Laboratory.
[] He X L, She T L, Gao F. 2016. A new method for picking up arrival times of seismic P and S waves automatically[J]. Chinese Journal of Geophysics, 59(7): 2519–2527. DOI:10.6038/cjg20160717
[] Kao H, Shan S J. 2004. The source-scanning algorithm: Mapping the distribution of seismic sources in time and space[J]. Geophysical Journal International, 157(2): 589–594. DOI:10.1111/gji.2004.157.issue-2
[] Liu H, Zhang J Z. 2014. STA/LTA algorithm analysis and improvement of Microseismic signal automatic detection[J]. Progress in Geophysics, 29(4): 1708–1714. DOI:10.6038/pg20140429
[] Liu J S, Wang Y, Yao Z X. 2013. On micro-seismic first arrival identification: A case study[J]. Chinese Journal of Geophysics, 56(5): 1660–1666. DOI:10.6038/cjg20130523
[] Ma Q, Jin X, Li S Y, et al. 2013. Automatic P-arrival detection for earthquake early warning[J]. Chinese Journal of Geophysics, 56(7): 2313–2321. DOI:10.6038/cjg20130718
[] Maeda N. 1985. A method for reading and checking phase time in auto-processing system of seismic wave data[J]. Zisin, 38(3): 365–379. DOI:10.4294/zisin1948.38.3_365
[] Montalbetti J F, Kanasewich E R. 1970. Enhancement of teleseismic body phases with a polarization filter[J]. Geophysical Journal International, 21(2): 119–129. DOI:10.1111/j.1365-246X.1970.tb01771.x
[] Mu L Y, Wu J P. 2015. Automatic onset phase picking for China array data set[J]. Progress in Geophysics, 30(6): 2511–2516. DOI:10.6038/pg20150610
[] Peng Z G, Zhao P. 2009. Migration of early aftershocks following the 2004 Parkfield earthquake[J]. Nature Geoscience, 2(12): 877–881. DOI:10.1038/ngeo697
[] Plešinger A, Hellweg M, Seidl D. 1986. Interactive high-resolution polarization analysis of broad-band seismograms[J]. Journal of Geophysics, 59(2): 129–139.
[] Poupinet G, Ellsworth W L, Frechet J. 1984. Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California[J]. Journal of Geophysical Research, 89(B7): 5719–5731. DOI:10.1029/JB089iB07p05719
[] Roberts R G, Christoffersson A, Cassidy F. 1989. Real-time event detection, phase identification and source location estimation using single station three-component seismic data[J]. Geophysical Journal International, 97(3): 471–480. DOI:10.1111/gji.1989.97.issue-3
[] Rubin A M. 2002. Using repeating earthquakes to correct high-precision earthquake catalogs for time-dependent station delays[J]. Bulletin of the Seismological Society of America, 92(5): 1647–1659. DOI:10.1785/0120010180
[] Rubin A M, Gillard D. 2000. Aftershock asymmetry/rupture directivity among central San Andreas fault microearthquakes[J]. Journal of Geophysical Research, 105(B8): 19095–19109. DOI:10.1029/2000JB900129
[] Savvaidis A, Papazachos C, Soupios P, et al. 2002. Implementation of additional seismological software for the determination of earthquake parameters based on MatSeis and an automatic phase-detector algorithm[J]. Seismological Research Letters, 73(1): 57–69. DOI:10.1785/gssrl.73.1.57
[] Schaff D. 2010. Improvements to detection capability by cross-correlating for similar events: A case study of the 1999 Xiuyan, China, sequence and synthetic sensitivity tests[J]. Geophysical Journal International, 180(2): 829–846. DOI:10.1111/gji.2010.180.issue-2
[] Schaff D P, Richards P G. 2004. Repeating seismic events in China[J]. Science, 303(5661): 1176–1178. DOI:10.1126/science.1093422
[] Shelly D R, Beroza G C, Ide S. 2007. Non-volcanic tremor and low-frequency earthquake swarms[J]. Nature, 446(7133): 305–307. DOI:10.1038/nature05666
[] Sleeman R, van Eck T. 1999. Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings[J]. Physics of the Earth and Planetary Interiors, 113(1-4): 265–275. DOI:10.1016/S0031-9201(99)00007-2
[] Stevenson P R. 1976. Microearthquakes at Flathead Lake, Montana: A study using automatic earthquake processing[J]. Bulletin of the Seismological Society of America, 66(1): 61–80.
[] Tan Y P, Deng L, Cao J Q, et al. 2016. Seismological mechanism analysis of 2015 Luanxian swarm, Hebei province[J]. Chinese Journal of Geophysics, 59(11): 4113–4125. DOI:10.6038/cjg20161115
[] Tan Y Y, Yu J, Feng G, et al. 2016. Arrival picking of microseismic events using the SLPEA algorithm[J]. Chinese Journal of Geophysics, 59(1): 185–196. DOI:10.6038/cjg20160116
[] Yomogida K. 1994. Detection of anomalous seismic phases by the wavelet transform[J]. Geophysical Journal International, 116(1): 119–130. DOI:10.1111/gji.1994.116.issue-1
[] Yoon C E, O'Reilly O, Bergen K J, et al. 2015. Earthquake detection through computationally efficient similarity search[J]. Science Advances, 1(11): e1501057. DOI:10.1126/sciadv.1501057
[] Zahradník J, Jansky J, Plicka V. 2015. Analysis of the source scanning algorithm with a new P-wave picker[J]. Journal of Seismology, 19(2): 423–441. DOI:10.1007/s10950-014-9475-7
[] Zhang M, Wen L X. 2014. Seismological evidence for a low-yield nuclear test on 12 may 2010 in north Korea[J]. Seismological Research Letters, 86(1): 138–145.
[] Zhang M, Wen L X. 2015. An effective method for small event detection: Match and locate (M&L)[J]. Geophysical Journal International, 200(3): 1523–1537. DOI:10.1093/gji/ggu466
[] 高淑芳, 李山有, 武东坡, 等. 2008. 一种改进的STA/LTA震相自动识别方法[J].世界地震工程, 24(2): 37–41.
[] 何先龙, 佘天莉, 高峰. 2016. 一种地震P波和S波初至时间自动拾取的新方法[J].地球物理学报, 59(7): 2519–2527. DOI:10.6038/cjg20160717
[] 贾瑞生, 谭云亮, 孙红梅, 等. 2015. 低信噪比微震P波震相初至自动拾取方法[J].煤炭学报, 40(8): 1845–1852.
[] 刘晗, 张建中. 2014. 微震信号自动检测的STA/LTA算法及其改进分析[J].地球物理学进展, 29(4): 1708–1714. DOI:10.6038/pg20140429
[] 刘劲松, 王赟, 姚振兴. 2013. 微地震信号到时自动拾取方法[J].地球物理学报, 56(5): 1660–1666. DOI:10.6038/cjg20130523
[] 刘希强, 周彦文, 曲均浩, 等. 2009. 应用单台垂向记录进行区域地震事件实时检测和直达P波初动自动识别[J].地震学报, 31(3): 260–271.
[] 刘希强, 周蕙兰, 李红. 2000. 基于小波包变换的地震数据时频分析方法[J].西北地震学报, 22(2): 143–146, 176.
[] 刘希强, 蔡寅, 赵瑞, 等. 2014. 区域地震信号自动识别方法及应用[J].应用地球物理:英文版, 11(2): 128–138.
[] 吕鹏, 丁志峰, 朱露培. 2011. 结合波形互相关的双差定位方法在2008年汶川地震余震序列中的应用[J].地震学报, 33(4): 407–419.
[] 马强, 金星, 李山有, 等. 2013. 用于地震预警的P波震相到时自动拾取[J].地球物理学报, 56(7): 2313–2321. DOI:10.6038/cjg20130718
[] 牟磊育, 吴建平. 2015. 利用单台触发检测台阵数据库事件[J].地球物理学进展, 30(6): 2511–2516. DOI:10.6038/pg20150610
[] 谭毅培, 邓莉, 曹井泉, 等. 2016. 2015年河北滦县震群发震机理分析[J].地球物理学报, 59(11): 4113–4125. DOI:10.6038/cjg20161115
[] 谭玉阳, 于静, 冯刚, 等. 2016. 微地震事件初至拾取SLPEA算法[J].地球物理学报, 59(1): 185–196. DOI:10.6038/cjg20160116
[] 王继, 陈九辉, 刘启元, 等. 2006. 流动地震台阵观测初至震相的自动检测[J].地震学报, 28(1): 42–51.
[] 武东坡. 2004. 震相识别的实时方法研究[硕士论文]. 哈尔滨: 中国地震局工程力学研究所.
[] 余建华, 李丹丹, 韩国栋. 2011. 特征函数响应特性分析及STA/LTA方法的改进[J].国外电子测量技术, 30(7): 17–19, 23, 27.
[] 赵大鹏, 刘希强, 刘尧兴, 等. 2013. 高阶统计量及AIC方法在区域地震事件和直达P波初动识别中的应用[J].地震地磁观测与研究, 34(5): 61–69.
[] 郑晨, 丁志峰, 周晓峰, 等. 2015. 利用波形互相关方法识别分析灌县-安县断裂重复地震[J].地震学报, 37(2): 299–311. DOI:10.11939/j.issn:0253-3782.2015.02.010