一种改进偏斜度的自动拾取S波震相方法与流程

本发明属于地震信号处理技术领域，具体涉及一种改进偏斜度的自动拾取s波震相方法。

背景技术：

天然地震s波震相到时的拾取是地震波研究的一项基础工作，是精确确定震源位置的必要参数。大量的地震记录，如果用人工拾取，则费时费力。由于受到p波震相的干扰，s波震相到时自动拾取准确度相对较低，研究的也相对较少。有论文作者通过判定偏斜度曲线极值点，以极值点前曲线斜率最大的位置作为微地震动p波震相到时。存在的问题是，偏斜度曲线总体来说是不断起伏下降的，期间会出现好多极值点，用手工分析、拾取、可以确定我们需要的极值点，但如果用计算机实现自动拾取，很难准确确定哪个是我们需要的极值点，这是拾取出现误差的重要原因；他们认为极值点前曲线斜率最大位置是地震p波震相到时，但根据我们的研究，p波震相和s波震相到时往往都在曲线斜率最大位置之后；他们在拾取前对记录进行了5-30hz带通滤波，然后对滤波后的记录进行偏斜度计算，但做了滤波，地震波震相到时会偏离原始位置，拾取的结果也会出现较大误差。

技术实现要素：

基于目前用偏斜度拾取地震波震相到时存在的问题，提出一种改进偏斜度的自动拾取s波震相方法，并应用于拾取震级大于ms4.0天然地震s波震相到时。一种改进偏斜度的自动拾取s波震相方法，具体流程如下：

步骤1：分别输入地震东西向、南北向地震记录；

步骤2：根据东西、南北向地震记录，用公式(1)计算特征函数xi(t)，即地震动水平向合振幅，其中，xew(t)为东西向地震记录，xns(t)为南北向地震记录：

步骤3：在整个地震记录时段内，用公式(2)、(3)和(4)计算特征函数改进偏斜度ki(t)：

其中，m是改进偏斜度计算区间内的记录点数，我们用的地震加速度记录采样时间间隔为0.005s，所以m也是该区间内记录时长；改进偏斜度是计算区间m内的偏斜度，规定为第一个记录点对应时刻的改进偏斜度；为改进偏斜度计算区间m内记录数据平均值，σ(t)为该区间m内的记录数据xi(t)的标准差，n为大于等于5的正整数；

步骤4：在整个地震记录时段内，用公式(5)计算改进偏斜度相对变化率绝对值，得到序列{ri(t)}：

其中，ki-1(t-1)为前一时刻改进偏斜度，ki(t)为当前时刻改进偏斜度；

步骤5：用公式(6)拾取改进偏斜度相对变化率绝对值序列{ri(t)}的最大值rmax，其对应的时刻就是s波初至时刻。

rmax＝max{ri(t)}(6)

有关改进偏斜度改进意义的描述：

目前偏斜度计算公式如下：

我们改进的偏斜度计算公式如下：

公式(7)，k(t)是偏斜度，公式(2)称ki(t)为改进的偏斜度；xi(t)是接收到的波动信号；m是记录点数，也就是记录时长，我们用的地震记录采样时间间隔为0.005s；改进偏斜度是记录点1到记录点m的偏斜度，实际是一个时段的偏斜度，我们规定，该时段的改进偏斜度定义为记录点1时刻的改进偏斜度，根据实践，我们m取200；为在采样时段m内记录数据平均值；σ(t)为时段m内的记录数据xi(t)标准差，n为正整数，改进偏斜度n≥5，目前的偏斜度n＝4，见公式(7)。在本发明方法中，n的取值可以变化，我们研究的地震震级≥ms4.0，n＝5，对于大地震，n越大，拾取到时的精度越高。

在没发生地震动时，接收到的是地脉动等信号，为干扰信号；发生地震时接收到的地震信号为有效信号。我们研究采用加速度记录仪单位为国际通用的cm/s²，地脉动信号的数量级在10^-6-10^-2cm/s²之间变化，平均处于10^-4cm/s²数量级；地震p波初至信号数量级一般在10^-2cm/s²左右，而ms4.0以上地震的s波初至合振幅数量级往往会达到10⁰cm/s²。公式(2)中σ(t)是标准差，其物理意义是表示信号的涨落程度，当接收到地脉动信号时，其涨落水平低，因此其改进偏斜度ki(t)较大；当接收到地震信号时，其涨落水平增大，改进偏斜度变小，甚至会变成负值；改进的偏斜度公式(2)没有目前偏斜度的分子项，突出了标准差的影响，其物理意义十分明确，即：当无地震信号输入时，涨落程度低，标准差小，其改进偏斜度大，而当有地震信号输入时，涨落加剧，标准差大，其改进偏斜度小，震相到时对应着偏斜度的突变时刻。在改进的公式(2)中用标准差的5次方，而不是用原来4次方的目的是增强标准差对改进偏斜度数值影响，当改进偏斜度的计算区间1～m一旦完全进入地震信号时段就会涨落加剧，标准差变大，改进偏斜度立即变小。需要强调的是由于改进偏斜度公式(2)中有个“-3”项，无地震信号输入时，标准差非常小，改进偏斜度数值是远大于“-3”这一数值，“-3”对改进偏斜度或为正、或为负起的作用很小；当改进偏斜度计算区间包含地震信号时，标准差变大，“-3”对改进偏斜度或为正、或为负起的作用增大。在改进偏斜度计算区间逐渐进入地震信号区间时，其标准差逐渐变大，改进偏斜度逐渐变小。当改进偏斜度计算区间完全进入、且刚好完全进入地震信号区间时，改进偏斜度有一个大幅度下降，此时相邻点间改进偏斜度相对变化率绝对值最大，可以用公式(5)计算地震记录时段内的改进偏斜度相对变化率绝对值，得到其相对变化率绝对值序列。在此基础上用公式(6)拾取改进偏斜度相对变化率最大的绝对值，其对应的t时刻就是s波震相到时。但如果地震s波初至比较强，改进偏斜度计算区没有完全被地震信号占据就满足公式(6)，此时拾取的s波会提前；如果地震s波初至比较弱，改进偏斜度计算区进入地震信号区一段时间才满足公式(6)，此时拾取的s波震相到时会落后。

rmax＝max{ri(t)}(6)

公式(5)是t时刻的改进偏斜度的相对变化率，在公式(2)中选择适当大的n可以使改进偏斜度计算区间一旦完全进入地震记录区，立即可以满足公式(5)和(6)。

关于公式(2)中n的作用需要进一步说明，当改进偏斜度计算区间处于非地震信号区间时，包括p波的干扰,，振幅一般小于10^-2cm/s²，其标准差小于1，大于等于5的n次方后使得标准差更小，其改进偏斜度非常大；当计算区间进入天然地震s波信号区间时，由于天然地震以剪切破裂为主，s波初至振幅很强，特别是ms≥4.0的，其水平向初至合振幅会达到10⁰cm/s²数量级，也就是1cm/s²数量级，其标准差往往会大于1，n次方后使得标准差变得更大，其改进偏斜度非常小。这样，改进偏斜度在进入地震信号区后，其值会一直减小。

由于地震p波震相受到的干扰相对小，目前拾取p波的诸多方法准确度较高；s波震相受到p波先到影响，其干扰相对较大，使得目前拾取s波准确度不高。我们提出的改进方法，由于天然地震s波初至绝对振幅大，特别是对于ms≥4.0的地震，往往会大于10⁰cm/s²数量级，即使与干扰波振幅相对差异不大，但由于干扰波振幅小于10⁰cm/s²数量级，标准差经过数次方后，干扰波的标准差会很小，改进偏斜度会很大；由于s波初至水平方向合振幅数量级大于10⁰cm/s²，其标准差多次方后迅速增大，改进的偏斜度迅速减小，与干扰波的差距变大，所以改进的方法很适合于拾取s波到时。

即使s波初至振幅达不到10⁰cm/s²数量级，而是10^-1cm/s²数量级，如果我们假设干扰波的振幅为10^-2cm/s²数量级，公式(2)中的n＝5，通过粗略计算可以得到干扰波的改进偏斜度大约为10¹⁰cm/s²数量级，而地震信号的改进偏斜度大约为10⁵cm/s²数量级，两者的相对差异也十分明显。

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

s波振动方向与传播方向垂直，其水平方向振动强，所以拾取s波震相采用两水平方向，即：东西、南北向地震记录作为地震动信号输入，用公式(1)计算水平向地震记录幅值，得到特征函数xi(t)，其中，xew(t)为东西向地震记录，xns(t)为南北向地震记录。

地震波震相到时定义为地震波的第一个波峰或波谷到时，所以我们用东西、南北向合振幅拾取地震波s震相到时，其物理意义是清楚的。有些论文作者用能量变化率做为特征函数拾取震相是不可靠的，因为在波峰或波谷处的能量变化率往往是不大的。

有益技术效果：

本发明提出一种改进偏斜度的自动拾取s波震相方法，改进偏斜度的主要部分与标准差的高次方成反比，突出了标准差对描述涨落的作用；改进方法标准差高次方的使用，使得改进偏斜度计算区间进入地震信号区时单调变化，也使得地震信号与非地震信号的改进偏斜度的差异明显，这种明显的差异使得我们可以通过改进偏斜度相对变化率确定震相到时，从而使得拾取震相到时的误差小；定义水平方向地震记录合振幅为特征函数，突出了s波到时的异常表现；改进方法拾取震相到时不进行任何滤波处理，不会使波形失真；此发明可以推广到p波震相到时的拾取。

附图说明

图1为本发明实施例的一种改进偏斜度的自动拾取s波震相方法流程图；

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

图3为本发明实施例的地震加速度东西向记录；

图4为本发明实施例的水平向合振幅；

图5为本发明实施例的改进偏斜度；

图6为本发明实施例的改进偏斜度的相对变化率曲线；

图7为本发明实施例的拾取s波到时误差；

图8为本发明实施例的拾取到时误差分析；

图9为本发明实施例的未改的偏斜度；

图10为本发明实施例的未改进偏斜度的相对变化比曲线；

图11为本发明实施例的改进偏斜度；

图12为本发明实施例的改进偏斜度的相对变化率曲线。

具体实施方式

下面结合附图和具体实施实例对发明做进一步说明，一种改进偏斜度的自动拾取s波震相方法，如图1所示，具体流程如下：

步骤1：分别输入地震东西向、南北向地震记录，其中，用xew(t)表示东西向地震记录，如图3所示，用xns(t)表示南北向地震记录，如图2所示；

步骤2：根据东西、南北向地震记录，用公式(1)计算特征函数xi(t)，即地震动水平向合振幅，，其中，xew(t)为东西向地震记录，xns(t)为南北向地震记录，如图4所示：

特征函数选择与allen不同，参考allenrv.automaticearthquakerecognitionandtimingfromsingletraces[j].bssa,1978,68(5):1521-1532，和allenrv.automaticphasepickers:theirpresentuseandfutureprospects[j].bssa,1982,72(6b):s225-s242。allen的特征函数选择了振幅幅值与振幅变化之和。由于地震波到时的定义是波峰或波谷处，波峰或波谷处振幅极大，但振幅的变化一般不是极大，allen的做法可能会是特征函数的极值点不是地震震相到时；我们特征函数选择与有关论文不同，他们用能量变化率做为特征函数拾取震相是不可靠的，因为在波峰或波谷处的能量变化率往往不是极大的。我们选择s波，亦称为剪切波或横波，以水平东西、南北向，这两分向合振幅作为特征函数，而不是有关论文中的东西、南北和竖直的三分量记录，突出了s波水平向振动强的特点，把东西、南北向的振幅合成后，也突出了震相到时的异常表现。

步骤3：在整个地震记录时段内，用公式(2)、(3)和(4)计算特征函数改进的偏斜度ki(t)，如图5所示：

改进偏斜度突出了没有地震信号输入时标准差小，偏斜度大、有地震信号输入时标准差大，偏斜度小，即当有无地震信号时，改进偏斜度差异明显，物理意义也十分明确，改进偏斜度表示信号涨落水平，此时选择n＝5。

步骤4：在整个地震记录时段内，用公式(5)计算改进偏斜度相对变化率绝对值序列{ri(t)}，如图6所示：

与有关论文作者不同，本发明s波震相到时不在偏斜度曲线斜率最大处，而在相邻改进偏斜度相对变化率绝对值最大处。

步骤5：用公式(6)拾取改进偏斜度相对变化率绝对值序列中的最大值，最大值对应的时刻就是s波初至时刻，如图6所示，s波初至时刻为18.975s。

rmax＝max{ri(t)}(6)

误差分析：

选取100个震级大于ms4.0地震，用手工拾取其s波到时做为标准，用我们改进偏斜度方法拾取s波到时，其绝对误差曲线如图7所示。当拾取的s波震相提前标准值时误差为负，否则为正。

从图7可以看出，拾取的误差都在0.1s内。我们用图8分析误差原因。我们定义震相到时为地震波波峰或波谷对应的时刻，即当改进偏斜度计算区间完全进入、且刚好完全进入地震信号区间时，改进偏斜度计算区的起始端所处的波峰位置，即图8中1时刻的位置，此时由于标准差变得很大，改进的偏斜度突然变得很小，相邻改进偏斜度相对变化率绝对值很大，最大值处对应着s波震相到时，此时准确地拾取了到时；如果地震波初至振幅相对较大，其改进偏斜度减小得很快，使得改进偏斜度计算区还没有完全进入地震信号区，改进偏斜度相对变化率绝对值已经达到最大值，如图8中2时刻对应的位置，此时拾取的s波震相到时提前；如果地震波的初至振幅相对较小，其改进偏斜度减小得慢，使得改进偏斜度计算区完全进入地震信号区、并超过图8中1时刻，改进偏斜度相对变化率绝对值才达到最大值，如图8中3时刻对应的位置，此时拾取的s波震相到时落后。

前已述及，我们定义震相到时为地震波波峰或波谷对应的时刻，这是我们统一的标准，但实际上在波峰或波谷之前地震波震相已经到达了，图8中，定义的s波震相是在1时刻到达，但2时刻震相已经到达了，所以没必要要求过于准确的到时，只要误差较小，我们改进的方法就是有意义的。

有论文作者用偏斜度拾取了p波到时，没研究s波到时拾取。我们把该方法推广，应用于震级大于ms4.0的破坏性地震s波震相拾取，提出了改进偏斜度拾取s波震相方法，而以往s波震相由于受到p波先到干扰使得拾取准确度不高，我们提出的方法丰富了震相拾取理论。

标准差对改进偏斜度的影响：

图9是未改进的偏斜度函数，由于函数用的幂次是n＝4，计算区间进入到地震信号时段，涨落增加相对不十分明显，在s波震相到时附近偏斜度不是一直减小，而是有一个偏斜度的增加。图10是未改进的偏斜度相对变化率绝对值曲线，该曲线说明为改进偏斜度变化起伏严重，期间出现很多极大值，且比值绝对值相差不大，给确定震相到时带来很多困难。图11是改进偏斜度曲线，n＝8，由于函数用的幂次是n＝8，计算区间进入到地震信号时段，涨落增加相对十分明显，在s波震相到时附近改进偏斜度的一直在减小。图12是改进偏斜度的相对变化率绝对值曲线，n＝8，在全程时段内，整个曲线只有一个十分明显的极大值。此处改进偏斜度总体不大，但由于进入地震信号区后，出现一个十分接近于零的改进偏斜度，使得其相对变化率绝对值最大。此例说明适当大的n值对准确拾取震相到时有益。