靳 平 张诚鎏 沈旭峰 王红春潘常周 严 峰 王电源

(中国西安710024西北核技术研究所)

基于信号整体与局部特征的地震数据自动处理方法研究*

(中国西安710024西北核技术研究所)

提出一种基于信号整体与局部特征的地震数据自动处理新方法,该方法不同于以往仅采用包络线互相关来直接检测事件. 新方法依然按照检测、 识别、 关联和定位等4个步骤进行处理,但在进行单个震相信号检测的同时,也检测信号波列并利用其包络线特征来识别和关联震相． 文中详细阐述了数据处理过程中如何定义一个波列及抽取和应用其特征． 相关的数据处理技术目前已成功应用于区域台网的日常数据处理分析中． 作为例子,给出了对新疆区域台网连续16天数据进行测试处理的结果． 实际应用结果表明,这种新方法可以大幅度降低自动处理结果的误检、 漏检率,在实际应用中具有很好的前景．

地震 自动数据处理 地震信号整体特征 地震信号局部特征

引言

可靠的地震数据自动处理方法对提高日常地震监测分析的效率和时效性具有重要意义． 特别是在过去的一二十年里,随着数字化地震仪的普及和通信技术的发展,普通地震台站的建设和运行维护成本明显降低,各地震台网的台站数目大量增加． 而增加的台站数目在提高相关台网的监测能力、 改善台网的区域覆盖性能的同时,也使得数据的分析处理量增加,其中包括需要检测分析的地震信号和地震事件的数目都大量增加． 在缺乏可靠准确的地震数据自动处理技术的情况下,这种急剧增加的数据量可能会使地震监测部门的日常监测分析工作面临严峻挑战．

常见的地震数据自动处理一般采用信号检测、 识别、 关联和定位的方式(singal detection, phase identification, association, and event localization, 简写为DIAL)进行(IDC,1999； 潘科等,2004； 金星等,2007； 侯建民等,2009； Küperkochetal,2010),其基本过程和方法一般可以简单归纳为： ① 信号检测与信号特征参数提取． 在对数据进行适当的预处理(如滤波等)之后,首先根据信号的某种特征,通常是信号幅值的瞬时变化,采用短时平均/长时平均比(STA/LTA)或其它一些方法(Allen,1978； Evans,Allen,1983； Baer,Kradolfer,1987； Joswig,1990； Kedrov,Ovtchinnikov,1990； Earle,Shearer,1994； Cansi,1995； Dai,Macbeth,1995； Withersetal,1998； Sleeman,van Eck,1999； Zhao,Takano,1999； 王海军,刘贵忠,2007； Selby,2008； Arrowsmithetal,2009)来检测信号,然后采用自回归-赤池信息准则(Autoregressive-Akaike’s information criterion,简写为AR-AIC)等一些更加精细的方法来估算信号的到时(Kamigaichi,1992； Inclán,Tiao,1994； Der,Shumway,1999； Saragiotisetal,2002； Tibuleacetal,2003； Zhangetal,2003； Stefanoetal,2006； Ahmedetal,2007； Hildyardetal,2008； Küperkoch,2010),并提取信号的幅值、 周期、 方位角、 慢度或偏振度等特征参数. ② 信号震相类型识别． 根据提取到的信号特征参数,采用推理或神经网络等模糊识别方法来判断检测信号的基本类型(Robertsetal,1989； Suteau-Henson,1991； Tong,1995； Wang,Teng,1995,1997； Anant,Dowla,1997； 张范民等,1998； IDC,1999； Bai,Kennett,2000,2001； 刘希强等,2000； Wang,2002； 李山有等,2006)(如噪声或信号、 P波或S波以及区域震或远震信号等). ③ 震相台站关联． 即在初步震相类型识别基础上,根据信号在时间、 方位角、 幅值等特征的相容性,将同一台站上来自同一事件的信号关联在一起,形成台站震相组(IDC,1999)． 区域震时可同时进行二次震相识别及单台定位. ④ 台网关联和定位． 在信号到时与幅值等特征相容性分析的基础上,采用定位测试的方法将不同台站上的信号关联在一起,并确定事件的发生时间和位置(IDC,1999)．

目前在维也纳国际数据中心(IDC)运行的地震数据自动处理系统是DIAL方法的典型代表． 其处理过程不但非常复杂,而且处理效果也很不理想． 例如,在IDC自动处理结果中,三分向台站上的震相误识率高达30%—50%(Research Required to Support Comprehensive Nuclear Test Ban Treaty Monitoring,1997),其标准事件列表(SEL3)中,通常有50%左右的虚假或分裂事件,并且事件目录(REB)中大约15%的事件是由分析员人工增补的(Pearceetal,2009； Procopioetal,2009； Le Brasetal,2011； Draelosetal,2012)． 虚假事件多数是因震相的识别和关联有误产生．

上述问题已在国际地震监测特别是全面禁止核试验条约核查地震监测技术领域引起广泛关注,并有多位学者尝试提出一些新的处理方法来解决这方面的技术问题(Aroraetal,2012； Draelosetal,2012； Mooreetal,2012)． 但这些新的尝试基本上还是基于IDC已有信号检测结果及其特征参数,其对提高自动处理结果准确可靠性的作用尚有待验证．

实际上,造成传统DIAL方法处理效果不佳的根本原因在于其自动处理技术仅利用了信号的局部特征． 局部特征是指利用很短时间窗口(通常是数秒)内信号提取的特征信息,如信号的幅值、 周期、 方位角、 慢度、 偏振度等． 尽管这些特征能较好地表征信号的瞬时特性,但难以反映信号的全貌及相互之间的联系． 然而,人工分析员判读地震图时却更多地依据整个地震图的包络线形状特征． 通常情况下,分析员首先根据整个地震图的包络线形状对其震中距类型(地方震、 区域震或远震等)作出判读,再从中检测不同的震相. 在地方震或区域震的情况下,根据它们在地震图中的不同位置对震相(如Pn,Pg,Sn和Lg等)进行判断．

受分析员判读地震图过程的启发,本文提出一种新的结合地震图整体特征和信号局部特征的地震信号自动处理技术． 该方法不同于以往有些学者提出的基于地震图包络线互相关来直接进行事件检测的方法(Ryzhikovetal,1996； Mooreetal,2012). 新的处理技术依然保留了检测、 识别、 关联和定位的过程,但在进行信号检测同时,引入了波列检测的概念,并巧妙地将信号检测与波列检测结合在一起. 在利用传统的局部信号特征的同时,依据波列的整体特征来进行震相的识别和关联,从而形成一种全新的基于信号整体特征与局部特征的地震数据处理方法． 我们利用新疆区域地震台网对该方法进行了测试,结果表明,该方法在降低自动处理事件的误检和漏检率、 提高定位精度等方面具有很好的效果．

1 地震图整体特征的定义和测量

图1给出了将地震图的整体特征与传统的DIAL处理方法相结合的原理． 在该方法中,我们采用普通的STA/LTA(STA和LTA分别表示信号的短时和长时平均)方法检测单个信号,即对时刻t,有

(1)

(2)

式中,x(i)为第i个采样点,N和L分别为STA和LTA对应的时间窗口长度. 实际处理过程中,依据Trnkoczy (1999)对STA和LTA的取值讨论,我们分别取为1 s和30 s． 为引入信号整体特征,在检测单个信号时,同步检测和识别各个连续波列,将属于同一个波列的信号自动地联系在一起． 为此,我们在信号检测过程中引入一个反映波列持续状态的变量ΘD,在没有波列持续时ΘD=0,有波列持续时ΘD=1,其具体计算方法如图2所示,即程序开始时将ΘD初始化为0,当有波列到达并触发其第一个检测时ΘD重置为1并使其保持到波列结束,即其信号幅值近似恢复到其到达之前的本底噪声水平时(具体算法见后面关于波列持续时间计算方法的说明),再将ΘD重置为0． 除波列开始时触发的检测外,波列持续过程中通常还有可能触发其它检测． 为在后续处理过程中判断各个信号检测是在一个持续波列的起始位置还是中间位置,我们在检测到一个信号时将该信号触发前的ΘD值作为其特征之一赋予该信号,用符号Ws来表示,并将其称作信号触发时的波列持续状态． 图2给出了对一段实际波形记录的信号检测结果,记录中的4个波列共触发了8个检测． 其中,D1,D3和D5等3个位于相应波列开始处检测的Ws特征值等于0,而其它5个位于有关波列中间的检测的Ws特征值等于1．

图1 地震信号整体特征及其与信号检测的关系

图2 波列持续状态计算方法示意图

我们进一步引入波列持续时间、 包络线峰值大小、 峰值信噪比及峰值延迟时间等特征参数来简单地描述波列特征. 为不增加额外的计算量,包络线简单地由STA构成的时间序列来表示． 仍如图1所示,定义一个波列的持续时间为从其第一个触发检测的到时t0开始到信号幅值重新衰减到本底噪声水平时的时间,即

(3)

式中tε为t0后STA/LTA0降低到某一阈值ε以下时所对应的时刻,而LTA0为t0时波列的第一个检测被触发时的LTA值．

为进一步说明信号包络线峰值大小、 峰值信噪比和峰值延迟时间等参数的定义,假定一个波列在其持续过程中共触发了D0,D1,…,Dn共n+1个检测,并记从t0到tε的整个波列的包络线峰值为P0,相应的到时记作tP0,则P0对应的包络线峰值大小AP0、 峰值信噪比SNRP0及峰值延迟时间TP0可分别定义为

(4)

类似地,记波列中位于检测Di-1与Di到时之间的包络线峰值为Pi,对应的到时为tPi,则相应的峰值特征参数的定义分别为

(5)

式中,ti-1,ti分别为Di-1和Di的到时． 在实际处理过程中,这些特征参数都可以在不同的频带范围内进行测量. 同时,为便于各个特征量的保存和调用,我们将P0的特征参数赋予检测D0,将Pi的特征参数赋予检测Di(i=1,2,…,n)．

利用上述包络线特征参数可以大致勾画出整个波列包络线的基本形状． 例如,对图1所示信号,根据4个触发检测的到时及各检测的信号峰值位置和峰值大小,并假定适当的包络线幅值上升和衰减模式,可以推测出相应的信号波列应具有类似于图3的包络线形状． 这对于确定信号的震中距类型和波列中除初始检测以外的其它检测的震相类型显然是非常重要的信息． 后面我们将对此作进一步的介绍和分析．

图3 根据信号检测的到时及各检测信号的峰值位置和峰值大小推测信号包络线形状

2 信号整体特征在数据自动处理过程中的应用

上述地震图包络线特征可以很容易地集成到通常的DIAL自动处理方法中,并在地震信号震相类型识别和关联过程中起到非常重要的作用． 首先,除极少数有两个或两个以上地震事件的信号波列在时间上相互重叠的情况外,对一个连续波列来说,在其持续时间范围内,被触发的各个检测都应该来自同一个地震事件,因而可以自然地将它们关联在一起． 其次,除多个地震的信号波列重叠在一起的情况外,正常情况下,直达P波,这里包括地方震时的Pg、 区域震时的Pn和远震时的P等,都应该是一个波列的起始检测,即相应的波列持续状态Ws应等于0． 同时,除非是信号信噪比非常低的情况,地方震和区域震的各种后续震相,包括地方震时的Sg和区域震时的Pg、 Sn和Lg等,都应该包含在以直达P为起始检测的同一个波列当中,它们的Ws值都应该等于1． 最后,对地方震和区域震来说,多数情况下整个波列的包络线最大峰值出现应该略晚于S或Lg波到时(图4). 特别是对经较低频带(如0.5—2.0 Hz)带通滤波后的水平分量记录尤其如此,即使是对激发S波能力相对较差的地下爆破来说一般也是这样．

图4 一次区域性地震在不同震中距台站的三分向信号记录 发震时间： 2010-06-01 04:29:29,震中位置： 38.4°N、 74.0°E， 震源深度： 156 km,震级： ML3.1． 为提高在远距离台站上的信噪比,波形全部经过 4—8 Hz带通滤波. 图中的P、 S线分别对应于P波和S波的理论到时

图5给出了实际应用的一个区域台网地震数据处理原型系统的流程原理示意图． 该流程与传统的DIAL处理方法从表面上看起来非常相似,也主要由信号检测、 信号特征参数提取、 信号震相类型识别、 信号台站关联、 台网关联和定位等技术环节组成． 所不同的是在信号检测和特征参数提取过程中增加了波列检测和波列特征参数提取的内容． 这里,波列特征参数在上一节已经进行了阐述,而信号特征参数则包括不同频带内的信号幅值、 信噪比、 信号优势频率、 信号方位角、 偏振特性或慢度大小等． 由于这些传统信号特征参数都是在相对较短的时间窗口内测量的,文中将其统称为信号局部特征．

图5 一个实际应用的同时基于信号和波列检测及其特征参数 的区域台网数据自动处理软件流程示意图

对于区域性地震事件的自动检测和定位来说,可靠地识别出信号中的P波、 S波震相并将其关联在一起非常关键,而本文引入的信号整体特征恰好在此方面可以发挥非常大的作用． 张诚鎏(2010)曾综合利用信号整体和局部特征,应用神经网络方法来对震相类型进行识别,结果表明对普通三分向台站,其准确识别率可达90%． 尽管如此,我们目前实际应用的系统中还是仅采用了信号局部特征来对触发检测进行分类,而信号整体特征则主要应用在信号的台站关联与区域震相二次识别过程中． 事实上,从实际应用情况来看,在引入了信号整体特征之后,台站关联前的震相类型识别已经不再必要,这也是图5中我们将震相类型识别用虚框来表示的原因．

图6是应用信号整体特征进行区域震相台站关联和识别的原理示意图． 首先,从触发的信号检测中搜索可能成为区域或地方性事件初至P波的信号． 正常情况下,这样的信号应该是一个波列的第一个触发检测,所以其触发时的波列持续状态Ws应等于0． 除此之外,实际处理时还可以利用其它一些信号特征,诸如信号水平分量与垂直分量的比值、 信号信噪比、 波列持续时间、 不同频带的信号幅值比等作进一步的筛选．

上述搜索初至P波的过程相当于搜索可能的地方或区域性事件波列． 对于这样的一个波列,正常情况下波列持续范围内的各个检测已经可以自然地关联在一起,而需要进一步完成的是找出其中的S波或Lg波,并与初至P波一起构成P-S(Lg)震相组,以便在此基础上进行单台定位并形成单台事件． 实际处理时,我们综合利用信号局部和整体特征并通过一系列的测试来检验识别波列中的主要S震相,同时也测试当前波列是否真的是一个区域或地方性事件所产生的信号． 测试对波列持续时间范围内的各个触发检测逐个进行,直到找出其中的主要S震相,其中包括地方震的Sg、 区域震的Sn或Lg． 测试的主要内容包括： ① 假设的S波检测与初至P波的幅值相容性,即要求假设的S波检测与初至P波检测的幅值比不大于特定阈值K,否则,则认为假设的两者并非来自同一个事件,而当前波列是由不只一个事件的信号混叠在一起而形成的．K的大小与频带和被测试的信号与初至P波信号的到时差有关,参考IDC数据处理方法(IDC,1999),多数情况下我们取K=40. ② 偏振特性相容性,即假设S波检测的偏振特性参数与其作为S类型震相是否相容. 统计表明S波的水平分量幅值明显比垂直分量幅值大(张诚鎏,2010),根据经验,处理时要求假设的S波水平分量幅值与垂直分量幅值比要大于1.5. ③ 距离相容性,即初至P波及假设S波信号频率成分应当与到时差tS-tP相一致． 由于tS-tP可以等效于震中距,而信号中不同频率成分的信号幅值随距离的衰减速度不同,使得震中距较小亦即tS-tP较小的信号中高频成分较丰富,而震中距或tS-tP较大的低频信号则更为显著． 实际处理时我们用不同频带内的信号幅值比来衡量高、 低频信号的相对强弱． 例如,在tS-tP<15 s的情况下,要求信号高频(8—16 Hz)与低频(0.75—1.5 Hz)频带内的信号幅值比不小于某一阈值η, 而阈值η则由经验确定． 另外,由于S波衰减通常会大于P波的衰减,所以,相容性测试时也要求前者的高、 低频信号幅值比不得超过后者的幅值比． ④ 包络线相容性,即要求至少存在若干分量及频带,在其上面测量的信号包络线峰值位置应该略晚出现于假设的S波到时的地方,即要求

图6 同时基于信号整体与局部特征的区域性P、 S震相对搜索方法示意图

(6)

式中：TP0为前面定义的整个波列包络线的峰值延迟时间；tS-tP为假设的S波与初至P波的到时差；τ为指定的适当阈值,可以与tS-tP的大小有关． 通过这一方法,可以可靠地区分主要S波震相与非初至P波震相(如区域震情况下的Pg),从而避免台站对P-S震相组的误识别,而这种误识别有可能在最后的处理结果中造成分裂事件或大的定位误差．

以上是综合应用信号整体与局部特征进行区域性地震事件台站关联与震相识别的基本思路． 实际应用过程中,还有几种特殊情况需要谨慎处理． 第一种特殊情况是两个事件的波列叠加在一起的情况． 在这种情况下,第二个事件的初至P在触发时的波列持续状态就不等于0而是等于1. 由于相互干扰,两个事件的波列特征参数测量结果都可能出现错误． 对于这类情况,多数情况下可以根据信号在幅值、 频率、 偏振和方位角(如果是台阵型台站)等方面的相容性对信号是否为多事件波列进行识别. 在此基础上将新事件的初至P波的触发状态修正为0,并对相应的包络线特征参数测量结果进行修正． 第二种特殊情况是关于一些震级较小的地震在震中距相对较大的区域台站上的记录． 由于其信号信噪比较低,P波波尾信号的幅值可能在主要S波震相(Sn或Lg)到达之前就已经衰减到本底噪声水平,使得P波与主要S波震相不能连在一起形成一个连续的波列(即主要S波震相的Ws等于0而不是1)． 对于这样的情况,如果是台阵型台站,只要P波与S波之间的方位角一致,且各自的慢度大小与相应的震相类型相符,则也可以将它们关联在一起,以提高对低震级事件的检测能力． 而对于普通三分向台站,由于其信号方位角、 慢度无法准确可靠地进行估算,这样的信号我们目前不允许在台站关联时进行关联,而是等到后面的台网关联时再来测试它们能否关联给已经形成的地震事件． 事实上,在各种区域性台网的正式地震目录中,基本上每一个能够可靠定位的区域性地震事件都会在不止一个台站上呈现出连续的、 包括主要P波和S波震相在内的波列． 因此,在台站关联时要求普通三分向台站上的区域性地震的P波和S波必须形成一个连续波列,在大幅度降低事件误检率的同时,并不会明显影响事件的检测率．

需要说明的是,在区域地震图的情况下,主要的S波震相既有可能是Sn,也有可能是Lg,其具体震相名称通常在台网关联时才能最后解决． 但是,如果波列中同时有Sn和Lg发育并都已经被检测到,在这种情况下,可以根据tSn-tP和tLg-tP的理论关系比较准确地对二者进行区分． 具体方法为,对波列持续时间范围内可能的S波检测,两两之间形成一个S-S震相对． 假定到时较早的S波为Sn,到时较晚的为Lg,计算它们与初至P波(此时即Pn)之间的到时差tSn-tP和tLg-tP． 假定这里的两个S检测的确分别为tSn-tP和tLg-tP，则二维空间中由坐标(tSn-tP,tLg-tP)所定义的观测点应该落在tSn-tP相对于tLg-tP理论关系曲线附近(图7)． 因此,找出离理论关系曲线距离最近的点,如果该距离小于规定阈值,则对应的两个信号可以分别判断为Sn和Lg. 而利用由最先到达台站的Pn、 Sn和Lg组成的台站三震相组合,则可以更加准确地进行单台地震定位．

图7 自动区分Sn和Lg原理示意图Fig.7 Illustration of the method for automatic identification of phases Sn and Lg

尽管本节介绍的实际数据自动处理系统的应用对象为区域性地震台网,但该系统并不局限于检测区域和地方性地震事件,同时也能比较准确可靠地检测和定位远震事件． 对于远震事件,通常情况下P和S等震相并不呈现为连续波列,甚至对占绝大多数的中小震级事件来说基本上只有直达P波信号才能被检测到． 尽管如此,本文引入的波列特征在远震信号的台网关联过程中依然可以起到很好的作用． 限于篇幅,这里不对利用区域性台网进行远震信号关联定位详细讨论,仅指出在正常情况下,要求被关联的各次远震P波的Ws值应等于0,仅这一点就可以大幅度提高关联搜索的效率． 此外,实际处理时除了可以根据信号的频率和偏振特征来测试待关联信号是否符合远震信号的特征外,信噪比与波列持续时间之间的关系也可发挥重要作用． 例如,如果一个波列的信噪比很低而持续时间很大,则一般不会是远震信号．

3 应用实例

为说明相关的数据自动处理技术在应用于较大规模的区域地震台网时也有效,我们用新疆区域地震台网连续半个月的数据来进行测试作为实例,介绍其测试处理结果．

图8为测试台网的台站分布图． 测试所用的数据从2010年6月1—16日共16天． 测试对台网可能记录到的所有事件,包括区域事件、 地方性事件和远震事件,自动进行检测的情况. 在此基础上,对检测到的事件人工逐个进行审核,以确认是否为虚假地震事件并统计其所占的比例． 为分析自动处理结果的漏检率和定位误差,我们对地方震和区域震事件与同样台网条件下完全由人工分析产生的由新疆地震局提供的地震目录(以下简称新疆局地震目录)进行了对比． 根据我们的分析,虽然新疆局地震目录实际上并不完整,许多可以检测和定位的事件并没有在目录中反映出来(其中包括一些采矿爆破被故意剔除)． 尽管如此,由于新疆局地震目录的产生过程与本文自动处理的产生过程完全独立,相当于是对所有事件的一个独立抽样,因此统计分析其中被自动处理程序检测到的事件的比例,足以客观地反映本自动处理系统对地方和区域性事件的检测能力．

图8 测试台网(新疆)台站分布

对上述16天内新疆区域地震台网记录到的波形数据,利用本文所介绍的地震数据自动处理原型系统共检测出3654个地震事件,其中包括地方震、 区域震、 远震和极远震． 根据人工分析结果,其中285个事件为虚假事件,约占全部事件的7.8%．

在这一时间段内,新疆局地震目录中共有1025个地震事件,其中ML≥2.0地震事件136个,ML1.0—2.0的事件325个,ML<1.0的事件48个,另有516个未给出定位结果的微破裂事件． 所谓未给出定位结果的事件是指只有一个台站记录到信号、 信号的tS-tP在数秒以内的事件． 对于这样的事件,新疆局地震目录的处理方法是,如果tS-tP<3 s，直接将记录台站的坐标作为事件的震中位置,而对于tS-tP>3 s的情况,则仅给出其发震时间而不给出具体的震中位置．

图9 本文方法对新疆局地震目录中2010年6月1—16日不同震级范围内地震事件检测结果的统计分析Fig.9 Statistic result of automatically detected rate for events listed in the Xinjiang seismic events satalog in 1—6 June, 2010

图9给出了对新疆局地震目录中不同震级范围内关于事件实际检测结果的统计分析. 图10给出了实际检测到的两个区域性事件的例子,其中第二个事件因为发生在台网以外,新疆局地震目录未给出其相关的发震参数． 图中给出了距离相应事件最近的15个台站记录到的波形,标出了对不同震相的检测、 识别和关联结果. 需要说明的是,按我们的实际处理程序,图10中给出的未关联或未定义检测的震相名称的确定并未用到信号整体特征,而仅仅是基于各种局部信号特征用神经网络方法来初步确定的． 对于ML≥2.0的事件,自动处理系统的检测率为100%； 对于ML1.0—2.0的事件,检测率为98.5%； 对于ML<1.0的事件,检测率为95.8%； 而对新疆局地震目录中未定位事件的检测率为52.1%． 对未定位事件的检测率较低的主要原因除了这些事件仅在一个台站上有信号以外,还因为实际测试的自动处理系统本身并不是针对这样的事件设计的． 例如,系统的信号检测模块在检测信号时要求在触发一个检测时其到时应距上一个检测3 s以上,而在新疆局地震目录中的516个未定位事件中,有174个事件的tS-tP小于这一数值．

需要指出的是,我们的自动处理系统还另外检测到许多新疆局地震目录中没有的区域和地方性事件,其中主要包括位置在新疆区域地震台网外的地震事件(图10b)、 工业爆破以及仅有一个或两个台站记录到的信号事件和一些信噪比很低、 单纯依靠人工分析难以检测的事件．

自动处理系统除具有良好的检测能力外,新的地震数据自动处理方法在自动定位结果的准确性方面也有比较好的表现． 图11给出了区域、 地方性事件的自动定位结果相对于新疆局地震目录中震中位置的误差分布． 对ML≥2.0的136个地震事件,系统的平均定位误差为15.6 km. 其中72个事件(占53%)的相对定位误差小于10 km,110个事件(占81%)的定位误差小于20 km,较大定位误差的地震基本都发生在网缘或网外． 考虑到新疆台网台站相对比较稀疏,上述136个ML≥2.0地震事件到最近新疆台站的距离平均约为70 km,部分事件震中位置到最近台站的距离达200 km左右,自动处理时能达到这一自动定位能力是可以接受的．

除有效检测区域、 地方性事件以外,以新疆台网这样的区域性地震台网的数据为输入,自动处理系统对远震和远区域震事件也表现出了比较好的检测定位能力． 由于对远震和远区域震的关联定位方法不同于对地方和区域震的关联定位方法,相关的方法和处理结果我们将另文阐述．

图10 典型区域性事件自动处理结果实例

(a) 2010-06-09 10:58:59新疆和静ML2.6地震,新疆局地震目录给出震中位置为42.33°N、 84.82°E,自动处理结果误差约11 km； (b) 2010-06-09 16:01:08阿富汗兴都库什山区mb3.4地震,IDC REB中给出的震中位置为35.89°N、 69.57°E,自动处理相对定位误差约31 km. 对每个事件,图中显示的是最近的15个台站经0.5—5 Hz带通滤波后的垂直分量波形,上面标出了对应台站的代码名称、 震中距及自动处理过程中实际触发的各个信号检测的到时位置及震相名称,其中正体字母表示最终关联给相应事件的定义震相,斜体字母代表未关联或未定义检测的震相,其中N表示被识别为噪声的检测．

Fig.10 Examples of automatically processing results for two regional events

(a) Seismograms of the southern XinjiangML2.6 earthquake occurred at 10:58:59 on 9 June 2010. The epicenter is (42.33°N, 84.82°E) in Xinjiang seismic events catalog (XJEC), our automatically determined epicenter is 11 km away from it. (b) Seismograms of anmb3.4 earthquake occurred in the Hindu Kush region, Afghanistan, at 16:01:08 on 9 June 2010. For the epicenter of this earthquake is outside of Xinjiang seismic network, no information is given in the XJEC. However, the IDC REB epicenter of this earthquake is (35.89°N, 69.57°E),and our automatically determined result is 31 km away from it. For each event, 0.5—5 Hz band-pass filtered vertical component seismograms from the 15 closest stations are displayed. For each seismogram, respective recording station and epicentral distance are given at its left side, and all detections triggered within the displayed waveform time window are added, where the labels marked by upright letters are automatically associated defining phases for corresponding event and the labels marked by italics are unassociated ones.Nstands for detections identified as noise by the software

图11 新疆局地震目录中ML≥2.0地震自动定位结果与人工编目结果 的比较(a)及震中相对误差分布(b) 空心圆为新疆局地震目录给出的震中位置,实心圆为本文自动定位结果,实线相连表示是同一个事件

4 结论

本文提出一种基于信号整体和局部特征的地震数据自动处理新方法. 以波列包络线特征为代表的地震信号整体特征对提高地震事件自动处理结果的准确可靠性具有非常重要的意义和作用． 结果表明,利用少数几个定量特征参数可以较准确地反映整个地震图的大致形状. 在此基础上通过引入波列检测的概念,将各个触发检测与相应的波列检测相关联,可以很好地解决地震信号整体特征与局部特征相结合的技术问题． 新的地震信号整体特征可以较好地反映不同地震信号之间的相互联系及它们各自在地震图中的相对位置,再将这些特征与传统的反映地震信号幅值大小、 频率高低、 偏振特性及传播方向的局部特征相结合,可以很好地解决地震信号关联和震相类型识别等方面的问题．

新的处理技术可以很方便地集成到传统的信号检测、 震相识别、 关联、 定位处理模式当中,并明显地降低地震事件的误检、 漏检率和平均定位误差． 实际的运行经验及利用新疆区域地震台网连续16天的数据进行测试的结果表明,基于该方法的自动处理结果能够在保持较高检测率的情况下将事件误检率控制在比较低的水平．

需要说明的是,我们对地震信号整体特征的应用还不够充分,许多判断条件的参数设置都是经验性的． 随着对信号整个特征参数分布特性研究的进一步深入和判断条件的进一步优化,预期可以取得更加令人满意的处理结果．

本研究得到新疆地震局的数据支持,在此表示感谢．

侯建民,黄志斌,余书明,黄静,代光辉,赵永. 2009. 中国国家地震台网中心技术系统[J]. 地震学报,31(6): 684--690.

Hou J M,Huang Z B,Yu S M,Huang J,Dai G H,Zhao Y. 2009. Technical system in National Earthquake Network Center of China[J].ActaSeismologicaSinica,31(6): 684--690 (in Chinese).

金星,廖诗荣,陈绯雯. 2007. 区域数字地震台网实时速报系统研究[J]. 地震地磁观测与研究,28(1): 64--72.

Jin X,Liao S R,Chen F W. 2007. The test running of Real-time Earthquake Information System for province-level seismic network[J].SeismologicalandGeomagneticObservationandResearch,28(1): 64--72 (in Chinese).

李山有,朱海燕,武东坡,宋晋东. 2006. 基于振幅和瞬时频率的震相自动识别方法[J]. 世界地震工程,22(4): 1--4.

Li S Y,Zhu H Y,Wu D P,Song J D. 2006. Automatic recognition of seismic phase based on amplitude and instantaneous frequency[J].WorldEarthquakeEngineering,22(4): 1--4 (in Chinese).

刘希强,周蕙兰,曹文海,李红,李永红,季爱东. 2000. 用于三分向记录震相识别的小波变换方法[J]. 地震学报,22(2): 125--131.

Liu X Q,Zhou H L,Cao W H,Li H,Li Y H,Ji A D. 2000. Identification method of seismic phase in three-component seismograms on the basis of wavelet transform[J].ActaSeismologicaSinica,22(2): 125--131 (in Chinese).

潘科,刘援朝,肖立萍. 2004. 辽宁数字地震台网地震参数自动处理系统的研制[J]. 东北地震研究,20(4): 51--56.

Pan K,Liu Y C,Xiao L P. 2004. Development on the automatic process system of earthquake parameter in Liaoning Digital Seismic Network[J].SeismologicalResearchofNortheastChina,20(4): 51--56 (in Chinese).

王海军,刘贵忠. 2007. 基于支持向量机的信号自动检测算法[J]. 地震学报,29(1): 85--94.

Wang H J,Liu G Z. 2007. Automatic signal detection based on Support Vector Machine[J].ActaSeismologicaSinica,29(1): 85--94 (in Chinese).

2016年冬2017年春，温室辣椒市场销售价格较往年降低不少，且长期在低价位徘徊。农户积极调整种植结构，改种荚豆、茄子、番茄、黄瓜等其他蔬菜，较为有利的避免了因市场价格过低造成的经济损失。顺市而变，变被动为主动，积极开展种植结构调整，努力确保设施蔬菜生产经济效益不滑坡。

张诚鎏. 2010. 地震信号识别和关联技术研究[D]. 西安: 西北核技术研究所: 1--71.

Zhang C L. 2010.StudyonSeismicPhasesIdentificationandAssociation[D]. Xi’an: Northwest Institute of Nuclear Technology: 1--71 (in Chinese).

张范民,李清河,张元生,盛国英,范兵. 1998. 利用人工神经网络理论对地震信号及地震震相进行识别[J]. 西北地震学报,20(4): 43--49.

Zhang F M,Li Q H,Zhang Y S,Sheng G Y,Fan B. 1998. The seismic signal and phase recognition by using artificial neural network theory[J].NorthwesternSeismologicalJournal,20(4): 43--49 (in Chinese).

Ahmed A,Sharma M L,Sharma A. 2007. Wavelet based automatic phase picking algorithm for 3-component broadband seismological data[J].JSEE:SpringandSummer,9(1): 15--24.

Allen R V. 1978. Automatic earthquake recognition and timing from single traces[J].BullSeismolSocAm,68(5): 1521--1532.

Anant K S, Dowla F U. 1997. Waveform transform methods for phase identification in three-component seismograms[J].BullSeismolSocAm,87(6): 1598--1612.

Arora N S,Given J,Tomuta E,Russell S,Spiliopoulos S. 2012. Analyst evaluation of model-based Bayesian seismic monitoring at the CTBTO[C]∥Proceedingsofthe2012MonitoringResearchReview:Ground-BasedNuclearExplosionMonitoringTechnologies. New Mexico: National Nuclear Security Administration, I(LA-UR-12-24325): 91--199.

Arrowsmith S J,Whitaker R,Katz C,Hayward C. 2009. TheF-detector revisited: An improved strategy for signal detection at seismic and infrasound arrays[J].BullSeismolSocAm,99(1): 449--453.

Baer M,Kradolfer U. 1987. An automatic phase picker for local and teleseismic events[J].BullSeismolSocAm,77(4): 1437--1455.

Bai C Y,Kennett B L N. 2000. Automatic phase-detection and identification by full use of a single three-component broadband seismogram[J].BullSeismolSocAm,90(1): 187--198.

Bai C Y,Kennett B L N. 2001. Phase identification and attribute analysis of broadband seismograms at far-regional distances[J].JSeism,5(2): 217--231.

Cansi Y. 1995. An automatic seismic event processing for detection and location: The PMCC method[J].GeophysResLett,22(9): 1021--1024.Dai H,Macbeth C. 1995. Automatic picking of seismic arrivals in local earthquake data using an artificial neural network[J].GeophysJInt,120(3): 758--774.

Der Z A,Shumway R H. 1999. Phase onset time estimation at regional distances using the CUSUM-SA algorithm[J].PhysEarthPlanetInt,113(1): 227--246.

Draelos T J,Ballard S,Young C J,Brogan R A. 2012. Refinement and testing of the probabilistic event detection,association,and location algorithm[C]∥Proceedingsofthe2012MonitoringResearchReview:Ground-BasedNuclearExplosionMonitoringTechnologies. New Mexico: National Nuclear Security Administration, I(LA-UR-12-24325): 221--231.

Earle P,Shearer P. 1994. Characterization of global seismograms using an automatic-picking algorithm[J].BullSeismolSocAm,84(2): 366--376.

Evans J,Allen S. 1983. A teleseismic-specific detection algorithm for single short-period traces[J].BullSeismolSocAm,73(4): 1173--1186.

Hildyard M W,Nippress S E J,Rietbrock A. 2008. Event detection and phase picking using a time-domain estimate of predominante periodTpd[J].BullSeismolSocAm,98(6): 3025--3032.

IDC. 1999. IDCProcessingofSeismic,Hydroacoustic,andInfrasonicData.IDCDocument[M/OL]. [1999-03-12]. http:∥www.rdss.info/librarybox/idcdocs/downloads/521.pdf.

Inclán C,Tiao G C. 1994. Use of cumulative sums of squares for retrospective detection in the changes of variance[J].JAmerStatistAssoc,89(427): 913--923.

Joswig M. 1990. Pattern recognition for earthquake detection[J].BullSeismolSocAm,80(1): 170--186.

Kamigaichi O. 1992. A fully automated method for determining the arrival times of seismic waves and its application to an on-line processing system, GSE/Japan/40[C]∥34thGSESession. July, 1992, Geneva.

Kedrov O K,Ovtchinnikov V M. 1990. An online analysis system for three-component seismic data: Method and preliminary results[J].BullSeismolSocAm,80(1): 2053--2071.

Küperkoch L. 2010.AutomatedRecognition,PhaseArrivalTimeEstimation,andLocationofLocalandRegionalEarthquakes[D]. Bochum: Fakultät für Geowissenschaften, Ruhr-Universität Bochum: 1--133.

Küperkoch L,Meier T,Lee J,Friederich W,EGELADOS Working Group. 2010. Automated determination of P-phase arrival times at regional and local distances using higher order statistics[J].GeophysJInt,181(2): 1159--1170.

Le Bras R,Russell S,Arora N,Miljanovic V. 2011. Testing and evaluation of the false events identification (FEI) and vertically integrated seismic association (VISA) projects[C]∥Proceedingsofthe2011MonitoringResearchReview:Ground-BasedNuclearExplosionMonitoringTechnologies. Arizona: National Nuclear Security Administration, I(LA-UR-11-04823): 313--321.

Moore D A,Mayeda K M,Myers S M,Seo M J,Russell S J. 2012. Progress in singal-based bayesian monitoring[C]∥Proceedingsofthe2012MonitoringResearchReview:Ground-BasedNuclearExplosionMonitoringTechnologies. Arizona: National Nuclear Security Administration, I(LA-UR-11-04823): 263--273.

Pearce R G,staff,Monitoring Data Analysis Section,IDC. 2009. Exploiting the skills of waveform data and analysts in the quest for improved automatic processing[C/OL]∥ISS09Conference. 10--12 June,2009,Vienna. [2009-06-20]. http:∥www.ctbto.org/fileadmin/user_upload/ISS_2009/Poster/DM-12A%20(PTS)%20--%20Robert_Pearce%20etal.pdf.

Procopio M J,Young C J,Lewis J E. 2009. Using machine learning to improve the efficiency and effectiveness of automatic nuclear explosion monitoring system[C]∥Proceedingsofthe2009MonitoringResearchReview:Ground-BasedNuclearExplosionMonitoringTechnologies. Arizona: National Nuclear Security Administration, Ⅰ(LA-UR-09-05276): 788--797.

Research Required to Support Comprehensive Nuclear Test Ban Treaty Monitoring. 1997.PanelonBasicRequirementsinSupportofComprehensiveTestBanMonitoring[M]. Washington DC: National Academy Press: 49--82.

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].GeophysJ,97(3): 471--480.

Ryzhikov G A,Biryulina M S,Husebye E S. 1996. A novel approach to automatic monitoring of regional seismic events[J].IRISNewsletter,XV: 12--14.

Saragiotis C D,Hadjileontiadis L J,Panas S M. 2002. PAI-S/K: A robust automatic seismic P phase arrival identification scheme[J].IEEETransGeosciRemoteSens,40(6): 1395--1404.Selby N D. 2008. Application of a generalizedFdetector at a seismometer array[J].BullSeismolSocAm,98(5): 2469--2481.

Sleeman R, van Eck T. 1999. Robust automatic P-phase picking: An online implementation in the analysis of broadband seismogram recordings[J].PhysEarthPlanetInt,113(1): 265--275.

Stefano R D,Aldersons F,Kissling E,Baccheschi P,Chiarabba C, Giardini D. 2006. Automatic seismic phase picking and consistent observation error assessment: Application to the Italian seismicity[J].GeophysJInt,165(1): 121--134.

Suteau-Henson A. 1991. Three-component analysis of regional phases at NORESS and ARCESS: Polarization and phase identification[J].BullSeismolSocAm,81(6): 2419--2440.

Tibuleac I,Herrin E,Britton J,Shumway R,Rosca C. 2003. Automatic determination of secondary seismic phase arrival times using wavelet transforms[J].SeismResLett, 74(6): 884--892.

Tong C. 1995. Characterization of seismic phases: An automatic analyzer for seismograms[J].GeophysJInt,123(3): 937--947.

Trnkoczy A. 1999. Understanding and parameter setting of STA/LTA trigger algorithm[G]∥IASPEI:NewManualofSeismologicalObservatoryPractice. Postsdam: GeoForschungsZentrum: 1003--1020.

Wang J. 2002. Adaptive training of neural networks for automatic seismic phase identification[J].PureApplGeophys,159(5): 1021--1041.

Wang J, Teng T-L. 1995. Artificial neural network-based seismic detector[J].BullSeismolSocAm,85(1): 308--319.

Wang J, Teng T-L. 1997. Identification and picking of S phase using an artificial neural network[J].BullSeismolSocAm,87(5): 1140--1149.

Withers M,Aster R,Young C,Beirger J,Harris M,Moore S,Trujillo J. 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection[J].BullSeismolSocAm,88(1): 95--106.

Zhang H,Thurber C,Rowe C. 2003. Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings[J].BullSeismolSocAm,93(5): 1904--1912.

Zhao Y,Takano K. 1999. An artificial neural network approach for broadband seismic phase picking[J].BullSeismolSocAm,89(3): 670--680.

A novel technique for automatic seismic data processing using both integral and local features of seismograms

(NorthwestInstituteofNuclearTechnology,Xi’an710024,China)

Reliable automation of data processing is of great significance for modern seismic monitoring. A novel technique for automatic seismic data processing which utilizes both integral and local features of seismograms was presented in this paper. However,unlike some previous efforts which seek to use envelope cross-correlation to detect seismic events directly,our technique keeps to follow the DIAL approach,i.e.,by steps of signal detection,phase identification,association and event localization. However,in addition to detect signals corresponding to individual seismic phases,the new technique also detects continuous wave-trains and explores their envelope features for phase type identification and signal association. More concrete ideas about how to define wave-trains and combine them with various detections as well as how to measure and utilize their features in the seismic data processing were expatiated in the paper. This approach has been applied to the routine data processing of regional seismic networks for several years,and as an application example,here were presented test results for a 16 days’ period using data from the Xinjiang regional seismic network. Practical application results show that the new technique can reduce both false alarm and missed event rate significantly and has good application prospects.

seismic; automatic data processing; integral features of seismic signals; local features of seismic signal

10.3969/j.issn.0253-3782.2014.03.012.

国家自然科学基金(41174033)资助.

2013-05-27收到初稿,2013-10-08决定采用修改稿.

e-mail: jinping668@sohu.com

10.3969/j.issn.0253-3782.2014.03.012

P315.61

A

靳平, 张诚鎏, 沈旭峰, 王红春, 潘常周, 严峰, 王电源. 2014. 基于信号整体与局部特征的地震数据自动处理方法研究. 地震学报, 36(3): 464--479.

Jin P, Zhang C L, Shen X F, Wang H C, Pan C Z, Yan F, Wang D Y. 2014. A novel technique for automatic seismic data processing using both integral and local features of seismograms.ActaSeismologicaSinica, 36(3): 464--479. doi:10.3969/j.issn.0253-3782.2014.03.012.