一种以振幅为特征函数的自动拾取P波震相方法与流程

本发明属于地震信号处理技术领域，具体涉及一种以振幅为特征函数的自动拾取p波震相方法。

背景技术：

目前常用的自动拾取地震动p波震相到时的主要有长短时均值比法，即sta/lta方法；aic准则；高阶统计计量方法，即pai-s/k方法；自回归方法定义进行地震震相判别的方法，除此之外还有分形维法和神经网络法等。

aic准则适合于地震波初至位置大致已知，并与其它方法结合使用才能较准确地拾取地震波震相。pai-s/k方法引进地震波形峰度和偏斜度函数，计算简单、速度快，但抗噪音能力不强。自回归方法将地震信号分为两个统计时段，取自相关最小值作为地震波震相到时，该方法拾取精度高，但计算速度慢。分形维数方法认为地震波与噪音信号叠加时分形维数将发生变化，并以此为依据拾取地震波到时，该方法的显著缺点是计算速度慢。神经网络法用峰值振幅、时窗内均方根振幅比、峰值与其前后峰值包络斜率及噪声与信号比值等作为神经网络输入条件、综合拾取p波到时。其优点是适用性广泛，缺点是计算用时长。

长短时均值比法是目前广泛使用的一种地震波初至拾取方法，其原理是取地震波特征函数在短时间内平均值与长时间内平均值之比确定地震波初至，即长短时均值比法。长短时均值比法具有算法简单、计算速度快和特征函数多样等优点，不足之处，一是目前特征函数的选择与地震波震相初至定义不吻合，使得到时的拾取在理论上就不够准确；二是短时间和长时间采用的时长不明确，而不同的时长，特别是短时间时长对到时的拾取结果有重要影响；三是做为分子的短时间项包含在分母中，而短时间项作用是反映异常的，出现在分母中减弱了异常表现；四是地震震相到时阈值确定原则不够明确。

技术实现要素：

针对以上技术问题，本发明提出一种以振幅为特征函数的自动拾取p波震相方法，具体流程如下：

步骤1：输入地震竖直方向记录xi(t)，i＝1...n，其中，xi(t)为第i个竖直方向加速度记录，其单位是cm/s²，n为地震记录长度；

步骤2：用公式(1)计算特征函数pi；

pi＝|xi(t)ⁿ|(1)

其中，pi为第i时刻xi(t)对应的特征函数，n是大于或等于4的整数；

步骤3：根据特征函数pi，用公式(2)计算长短时均值比序列{pi}；

其中，pi是与地震记录相关的特征函数、{pi}是经过长短时均值比变换得到的序列、i是地震记录点序号，n4为长短时均值比方法的短尺度，n3为长短时均值比方法的长尺度；j是计算过程中叠加运算需要的地震记录序号。

步骤4：在长短时均值比序列{pi}中，判断长短时均值比计算结果，即判断序列{pi}中第i时刻长短时均值比是否大于或等于阈值，如果大于或等于阈值，继续步骤5，如果小于阈值，则i＝i+1，判断下一时刻长短时均值比结果是否大于或等于阈值；

步骤5：在长短时均值比序列{pi}中,从长短时均值比结果大于或等于阈值的对应时刻开始，按时间顺序，寻找长短时均值比结果的极大值，该极大值对应的时刻就是p波到时；具体方法是：在长短时均值比结果达到阈值后，标记该时刻，把该时刻的长短时均值比结果与下一时刻的长短时均值比结果做比较，如果下一时刻的长短时均值比结果小于或等于该时刻的长短时均值比结果，则该时刻的长短时均值比结果对应的时刻就是p波到时；如果下一时刻长短时均值比结果大于该时刻的长短时均值比结果，则i＝i+1,重复步骤5，直到找到极大的长短时均值比结果，该极大的长短时均值比结果对应的时刻就是p波震相到时。

所述特征函数仅用竖直方向地震记录参与计算，东西、南北向地震记录不参与计算。

所述步骤4中的阈值为阈值下限的2倍，阈值上下限定义为：定义用长短时均值比拾取有效波即p波的第一个峰值平均值为拾取到时阈值的上限，拾取干扰波最大峰值平均值为阈值下限，明确了阈值的确定方法。

所述步骤3中公式(2)中分子和分母的特征函数是不同时段的、是不互相包含的。目前地震波震相自动拾取常常用到长短时变换，其原理用公式(3)表示，参考论文：allenr.v.automaticearthquakerecognitionandtimingfromsingletraces[j].bssa,1978,68(5):1521-1532.及allenrv.automaticphasepickers:theirpresentuseandfutureprospects[j].bssa,1982,72(6b):s225-s242。

公式(3)中pi是与地震记录相关的特征函数、{pi}是经过长短时变换得到的序列、i是地震记录序号、n4为长短时均值比方法的短尺度，n3为长尺度，公式(3)中，n3与n4的数值差异较大。地震记录可以是加速度、速度和位移，由于地震加速度变化相对迅速，有利于震相的自动拾取，本发明重点研究在加速度记录中拾取p波震相初至。

本发明提出用公式(2)进行长短时变换。

注意到公式(2)和(3)中的n4和n3的不同，在公式(2)中，两者数值差异很小，在公示(3)中两者数值差异大。

公式(3)和(2)中，分子是短时项(求和的数据点少)；分母是长时项(求和的数据点多)，公式(3)中短时项包含(n3-n4+i)～(i+n3)之间的数据，长时项包含i～(i+n3)之间的数据；公式(2)中短时项包含(n4+1)～n3之间的数据，长时项包含i～(i+n4)之间数据。公式(3)中的短时项中的数据包含在长时项数据中，而公式(2)中短项数据不包含在长时项数据中。公式(2)中短时项是作为异常项出现的，而长时项是作为正常情况下的正常项出现的，其物理意义是异常项与正常项的均值比。公式(3)分母长时项包含了异常项，其物理意义不够明确。理论上公式(3)和(2)在没有异常输入情况下比值在1左右波动，当有异常波输入时，公式(3)中分母也包含了异常项，而公式(2)中分母不包含异常项，公式(2)计算得到的数值一定大于公式(3)的计算结果，所以公式(2)有更强表现异常的能力。

公式(2)和(3)分子和分母都涉及除以数据点的数目，以公式(3)为例，其分母中除以了n3，但该公式分母求和的数据点是n3+1。类似情况还出现在标准差的定义中，有的学者用的是n，即数据数目，有的用的是n-1。在我们的应用中怎么定义都不影响计算结果，定义的不同仅仅使结果差一个常数。

规定公式(2)中的分母是作为未发生地震时的微振动信号记录出现的，可以在地震波未到达前计算完毕，在实际应用时可以不参与运算，从这个意义上讲，该公式的分母是一个常数。但该发明如果应用于地震预警必须参与实时运算，即输入实时记录。在本发明的工作中，我们通过预留在100个地震记录中的微振动信号，分别取10s时长计算平均值得到。需要说明的是，由于地震记录预留信号长度会远大于n3的长度，短时项取在n3的后段，不会错过p波震相到时。

在实际工作中由于各种随机干扰的存在，拾取震相到时需要确定一个阈值。由于阈值的确定与特征函数有密切关系，所以我们首先选择进行长短时变换的特征函数。

震相到时定义在波峰或波谷处，p波也称纵波到时是地震记录p波第一个波峰(波谷)对应的时刻。根据定义我们取如下与振幅有关的函数p为特征函数：

pi＝|xi(t)ⁿ|(1)

其中，xi(t)为第i个竖直方向加速度记录，n是大于或等于4的整数。与allen提出的特征函数相比，我们的仅有振幅，而不是allen的既有振幅，又有振幅的变化，即振幅平方与振幅变化平方和。这是因为我们在研究过程中发现，振幅达到峰值时，振幅的变化并不一定大，用振幅和振幅变化同时作为特征函数表现异常会使异常分散；我们的特征函数也有别其他论文的能量变化率。这是因为在p波或s波到达的波峰或波谷的能量变化率未必是最大的。公式(1)使得地震波震相到时就是特征函数对应的最大值，所以确定地震波到时的问题变成寻找与特征函数相关的最大值问题。公式(2)相当于一个有异常记录点特征函数与正常情况下特征函数比值，当有地震波输入时，该比值会明显大于1；没有地震波输入时，该比值在1附近波动。在我们的工作中公式(1)的特征函数采用n＝4的形式，小于1的数四次方后会更小，大于1的数四次方后数值会更大，这样地震波到达时，用四次方特征函数会突出异常表现。关于这一点，随后在具体实施方式中还会进一步阐述。

在实际工作中，公式(2)中的n3＝300，n4＝295，作为分子的短时间n4+1～n3取的数据点较少，作为分母的长时间取的数据点较多。短时间取得短的原因是考虑到长短时均值比对异常的敏感性，短时间取得太长使比值失去了敏感性。我们定义长短时均值比是短时间第一个记录对应时刻的长短时均值比。

地震波在均匀、各向同性介质中传播时是沿直线传播的。但地壳介质是层状分布的，离地面越近，介质密度越小，波速也变小。按折射定律地震波入射角正弦与波速成正比，当地震波从深部传向地表时，波速的减小会使地震波的入射角减小，使得在地表接收到的地震波都是几乎垂直地面入射的。p波的传播方向与振动方向平行，故拾取p波只采用竖直方向记录。

拾取p波到时，需要确定一个触发阈值，如果长短时均值比得到的数值达到阈值，即大于或等于阈值，立即在长短时均值比序列中按时间顺序寻找长短时均值比的极大值，其对应的时刻就是p波震相到时。具体方法是在长短时均值比达到阈值后，标记该时刻，把该时刻的长短时均值比与下一时刻的长短时均值比做比较，如果下一时刻的数值小于或等于该时刻的，则该时刻的长短时均值比数值对应的时刻就是p波到时；如果下一时刻数值大于该时刻的，重复上述步骤，直到找到极大的长短时均值比，其对应的时刻就是p波震相到时。

有益技术效果：

本发明特征函数定义与地震波p波震相到时的定义一致，使得寻找震相到时问题转化成寻找与特征函数极大值相关的问题；无地震信号输入时，长短时均值比在1附近波动，有地震信号输入时，比值会明显增大，明显大于1，在我们的工作中，公式(1)中选定n＝4，突出了异常表现；重新定义的长短时均值比的长、短时项不产生交叉，增强了长短时均值比表现异常能力；根据p波的在竖直方向振动强的物理事实，只用竖直方向地震记录形成特征函数；明确了p波到时触发阈值的确定方法，减少了错误触发或遗漏的机会。本发明拾取p波到时误差小，可以减轻人工拾取的工作量；如果输入实时地震记录，并进行相应的实时计算，去掉步骤5，就可以应用于地震预警工作，即长短时均值比到达阈值，就可以触发地震预警系统；本发明在确定s波到时阈值后，可用于拾取s波震相到时自动拾取。

附图说明

图1为本发明实施例的一种以振幅为特征函数的自动拾取p波震相方法流程图；

图2为本发明实施例的地震加速度竖直向记录；

图3为本发明实施例的p波一次特征函数；

图4为本发明实施例的p波四次特征函数；

图5为本发明实施例的用一次特征函数拾取的p波到时；

图6为本发明实施例的用四次次特征函数拾取的p波到时；

图7为本发明实施例的用一次特征函数拾取p波得到的阈值；

图8为本发明实施例的用四次特征函数拾取p波得到的阈值；

图9为本发明实施例的拾取p波震相到时误差。

具体实施方式

下面结合附图和具体实施实例对发明做进一步说明，本发明提出一种以振幅为特征函数的自动拾取p波震相方法，如图1所示，具体流程如下：

步骤1：输入地震竖直方向记录xi(t)，i＝1...n，其中，xi(t)为第i个竖直方向加速度记录，n为地震数据的记录长度；

步骤2：用公式(1)计算特征函数pi；

pi＝|xi(t)ⁿ|(1)

其中，pi为第i时刻xi(t)对应的特征函数，n是大于或等于4的整数；

步骤3：用公式(2)计算长短时均值比序列{pi}；

其中，pi是与地震记录相关的特征函数、{pi}是经过长短时变换得到的序列、i是地震记录序号，n4＝295为长短时均值比方法的短尺度，n3＝300为长短时均值比方法的长尺度；j是计算过程中迭代用的地震记录序号。

步骤4：在长短时均值比序列{pi}中，判断长短时均值比结果，即判断序列{pi}中第i时刻长短时均值比是否大于或等于阈值，如果大于或等于阈值，继续步骤5，如果小于阈值，则i＝i+1，判断下一时刻长短时均值比结果是否达到阈值；

步骤5：在长短时均值比序列{pi}中,从长短时均值比结果大于或等于阈值的对应时刻开始，按时间顺序，寻找长短时均值比结果的极大值，该极大值对应的时刻就是p波到时；具体方法是：在长短时均值比结果大于或等于阈值后，标记该时刻，把该时刻的长短时均值比结果与下一时刻的长短时均值比结果做比较，如果下一时刻的长短时均值比结果小于或等于该时刻的长短时均值比结果，则该时刻的长短时均值比结果对应的时刻就是p波到时；如果下一时刻长短时均值比结果大于该时刻的长短时均值比结果，则i＝i+1,重复步骤5，直到找到极大的长短时均值比结果，该极大的长短时均值比结果对应的时刻就是p波震相到时。

如图2是某次地震某台站的竖向加速度记录，采样时间间隔为0.005s。通过人工拾取得到其p波到时为16.99s，位于3398个采样点。图3是n＝1时的特征函数；图4是n＝4时的特征函数。从上述两图可以看出，无论是一次特征函数，还是四次特征函数都清晰地反映了p波的到时，都处于振幅峰值对应的时刻，但相对而言，一次特征函数的p波到时峰值相对干扰差异小，而四次特征函数的p波到时峰值相对干扰差异大。这就是选择四次函数作为拾取p波到时特征函数的原因之一。

图5是以一次函数，即n＝1，作为特征函数，用长短时均值比拾取的p波到时；图6是以四次函数，即n＝4，作为特征函数，用长短时均值比拾取的p波到时。两种方法都准确地拾取了p波到时，都位于峰值对应的时刻。但是，用一次函数拾取的p波到时对应的峰值(约10cm/s²)与干扰波的峰值(约4cm/s²)差异小，两者峰值在同一个数量级上，不利于确定p波到达的阈值，如图5。而用四次函数拾取的p波到时对应的峰值(约1500(μm/s²)⁴)与干扰波的峰值(约47(μm/s²)⁴)差异十分明显，容易确定p波到达的阈值，如图6。这是选择四次函数为特征函数的另外一个原因。

以四次函数作为特征函数，确定了p波到时的阈值后，按时间次序搜寻长短时均值比的极大值，其对应的时刻就是p波到时，如图6。

所述阈值的确定：

拾取p波到时到时还有一个关键问题，就是阈值的选择问题，阈值的选择在一定程度上决定了拾取的成败。我们把要拾取的p波称为有效波，其大量地震记录对应的长短时均值比第一个峰值的平均值称为阈值上限；在有效波之前到达的波称为干扰波，其大量记录对应的长短时均值比最大值的平均值称为阈值下限。阈值确定的难易程度取决于有效波和干扰波峰值的差异，差异越大，越容易确定阈值。

采用用手工可以拾取到p波到时的地震记录100条，分别用一次函数和四次函数作为特征函数，用长短时均值比法拾取p波到时，并得到阈值上、下限。图7是用一次特征函数拾取p波到时得到的阈值上、下限。可以看出，p波到时阈值上、下限有一定差异，100个地震阈值上限平均值为10.2991，下限平均值为4.2407，但两曲线有交叉部分，从上、下阈值曲线中确定一个p波触发阈值比较困难。

图8是用四次特征函数拾取p波到时得到的阈值上、下限。可以看出，p波到时阈值的上、下限差异明显，100个地震阈值上限平均值为1563.3529，下限平均值为48.7805，两者差异明显，接近2个数量级，从上、下阈值曲线中确定一个p波的触发阈值比较容易。在本发明工作中，暂定拾取p波到时阈值为下限阈值的2倍，约等于97.5。

一种以振幅为特征函数的自动拾取p波震相方法的误差：

选取100个地震，用手工拾取其p波到时做为标准，用我们改进的方法拾取p波到时，其绝对误差曲线如图8所示。当拾取的p波震相提前标准值时，误差为负，否则为正。

从图9可以看出，拾取的误差全部为负或等于0，说明拾取的到时只有准确的和提前的，没有落后的。这是因为定义p波震相到时为波峰或波谷对应的时刻，长短时均值比达到或超过阈值后，在搜取p波到时对应的极大值时，遇到了小的波峰或波谷，程序会判定其对应时刻为p波到时。

人们人为定义波峰或波谷对应的时刻为震相到时，但实际上在波峰或波谷到达前，地震波已经到达了。所以没必要过分强调拾取到时的准确性，只要误差可以接受，工作就是有意义的。

本发明拾取震相到时误差小的另外一个原因是分子取得数据点少n3-n4＝5,取第一个数据点对应时刻为长短时均值比时刻。以往分子数据点取得多，无论取哪个数据点对应时刻为长短时均值比时刻，都会带来较大误差。