CN110018515A - 一种用振幅变化长短时均值比实现地震预警中p波自动拾取的方法 - Google Patents

Abstract

本发明提出一种用振幅变化长短时均值比实现地震预警中P波自动拾取的方法，流程包括：实时输入地震加速度数据记录；通过“振幅变化绝对值”计算特征函数；对特征函数序列进行长短时均值比、用阈值判断是否触发地震预警。本发明改进方法引进“振幅变化”绝对值作为特征函数，具有更强的捕捉地震波到时能力；以“振幅变化绝对值”作为特征函数计算时调用函数少，计算速度快，P波到时计算几乎不占用机时，提高了地震预警的效率；在较大振幅噪音条件下，抗噪音能力强；在较小振幅噪音条件下，拾取的准确度略高。

Description

一种用振幅变化长短时均值比实现地震预警中P波自动拾取 的方法

技术领域

本发明属于地震信号处理领域，具体涉及一种用振幅变化长短时均值比实现地震预警中P波自动拾取的方法。

背景技术

地震预警是针对地震预报能力不足提出的震时预警方法。其一般原理是，把地震探测器安置在活断层附近，当地震发生时，由于P波传播速度比S波快，地震探测器首先接收到P波，在破坏力很强的S波到达前用电磁波把地震消息传达到城市，通知人们在地震S波到达前做好应急准备。最近几年发展迅速的地震预警技术，需要自动、快速、准确确定P波到时。

长短时均值比方法是目前广泛使用的一种地震波初至拾取方法(据参考文献：Allen R V.Automatic earthquake recognition and timing from single traces[J].BSSA，1978，68(5)：1521-1532；Allen R V.Automatic Phase pickers：their Presentuse and future prospects[J].BSSA，1982，72(6B)：S225-S242.)，其原理是取地震波特征函数在短时间内平均值与长时间内平均值之比确定地震波初至，即长短时平均值法，地震波初至时刻就是用长短时平均值法计算得到超过阈值的时刻。长短时均值比方法具有算法简单、计算速度快和特征函数多样等优点，这非常符合地震预警的要求。其不足是特征函数的不同选择使得到时的拾取产生差异，短时间和长时间取得过短、过长都会造成地震触发错误。

本发明为一种用振幅变化长短时均值比实现地震预警中P波自动拾取方法，具体是对目前用长短时均值比实现地震P波自动拾取方法的改进。

发明内容

本发明基于以上技术问题，提出一种用振幅变化长短时均值比实现地震预警中P波自动拾取的方法，解决了原方法中用“振幅变化平方”作为特征函数拾取P波，捕捉地震波到时能力不强的问题；计算相对繁琐，地震预警效率相对较低的问题；在较大振幅噪音条件下，已经不能拾取P波到时的问题；在较小振幅噪音条件下，拾取准确度较低的问题。改进后的方法引进“振幅变化”绝对值作为特征函数，具有更强的捕捉地震波到时能力；计算时调用函数少，计算速度快，P波到时计算几乎不占用机时，提高了地震预警的效率；在较大振幅噪音条件下，抗噪音能力强；在较小振幅噪音条件下，拾取的准确度略高。

一种用振幅变化长短时均值比实现地震预警中P波自动拾取的方法，具体流程如下：

步骤1：实时输入竖直分量地震加速度记录流；

步骤2：计算地震数据的特征函数，具体公式如下：

P_i＝|x_i|+|x_i+1-x_i|，0≤i≤N-2

其中，{x_i}为地震记录的离散序列，N为序列长度，i为每个地震记录的离散序列编号；

步骤3：确定长、短窗长度后，按确定的步长同步移动长、短窗，同时实时对特征函数计算STA/LTA数值，具体公式如下：

其中，{P_i}为对序列Pi进行长短时均值比计算得到序列，N₄为长短时均值比方法的短尺度，N₃为长尺度；

步骤4：判断STA/LTA是否超过阈值P_c，具体公式如下：

P_n＝max{Pi}≥P_c

若是则取其对应的采样时间为地震波到时，同时触发地震预警装置；否则返回到步骤1。

有益技术效果：

(1)本发明提出一种用振幅变化长短时均值比实现地震预警中P波自动拾取的方法，改进方法以“振幅变化绝对值”作为特征函数比原方法以“振幅变化平方”作为特征函数具有更强的捕捉地震波到时能力；

(2)原方法由于以“振幅变化平方”作为特征函数，计算时相对繁琐，地震预警效率相对较低，新方法以“振幅变化绝对值”作为特征函数，计算时调用函数少，计算速度快，P波到时计算几乎不占用机时，提高了地震预警的效率；

(3)在较大振幅噪音条件下，原方法已经不能拾取P波到时，但新方法在P波到时处的异常依然明显，说明改进方法的抗噪音能力强；

(4)在较小振幅噪音条件下，两种方法大都在允许的误差范围内拾取了P波到时，但改进方法拾取的准确度略高。

(5)新方法确定一个触发阈值，地震预警工作中可以排除噪音影响，阈值的确定应考虑仪器和环境因素，还要考虑地震的特点。

附图说明

图1为本发明实施例的地震预警流程图；

图2为本发明实施例的地震信号(其中(1)是噪音信号，(2)是竖直方向地震信号)；

图3为本发明实施例的对图2地震动计算得到的振幅变化(图(1))和长短时均值比曲线(图(2))；

图4为本发明实施例的地震信号(其中(1)是噪音信号，(2)是竖直分量地震信号)；

图5为本发明实施例的噪音振幅为10^-5cms^-2时用原方法和新方法拾取的P波到时(其中(1)和(2)分别是用原方法得到的“振幅变化平方”曲线和对应的长短时均值比曲线；(3)和(4)分别是用改进方法得到的“振幅变化”曲线和对应的长短时均值比曲线)；

图6为本发明实施例的噪音振幅为10^-4cms^-2时用原方法和新方法拾取的P波到时(其中(1)和(2)分别是用原方法得到的“振幅变化平方”曲线和对应的长短时均值比曲线；(3)和(4)分别是用改进方法得到的“振幅变化”曲线和对应的长短时均值比曲线)；

图7为本发明实施例的噪音振幅为10^-3cms^-2时用原方法和新方法拾取的P波到时(其中(1)和(2)分别是用原方法得到的“振幅变化平方”曲线和对应的长短时均值比曲线；(3)和(4)分别是用改进方法得到的“振幅变化”曲线和对应的长短时均值比曲线)；

图8为本发明实施例的噪音振幅为10^-5cms^-2时用原方法和新方法拾取的P波到时偏差(其中(1)和(2)分别是噪音振幅10^-3cms^-2时，用原方法得到的P波自动拾取统计偏差点状图及统计直方图；(3)和(4)分别是噪音波振幅为10^-3cms^-2时，用新方法得到的P波自动拾取统计偏差点状图及统计直方图)；

图9为本发明实施例的噪音振幅为10^-4cms^-2时用原方法和新方法拾取的P波到时偏差(其中(1)和(2)分别是噪音振幅10^-3cms^-2时，用原方法得到的P波自动拾取统计偏差点状图及统计直方图；(3)和(4)分别是噪音波振幅为10^-3cms^-2时，用新方法得到的P波自动拾取统计偏差点状图及统计直方图)；

具体实施方式

下面结合附图和具体实施实例对发明做进一步说明，以某次地震某一台站竖直方向加速度数据记录为例，说明改进方法预警地震的步骤如图1所示，一种用振幅变化长短时均值比实现地震预警中P波自动拾取的方法，具体流程如下：

步骤1：实时输入竖直分量地震加速度记录流，如图2所示：图2(1)是在地震波到达前2s施加的噪音信号，噪音信号幅值小于地震波初至竖直方向振幅一个数量级；图2中(2)是采样周期为0.005s的地震竖直方向加速度记录，P波初至时刻为2s；

步骤2：计算地震数据的特征函数，如图3(1)“振幅变化”曲线所示，图具体公式如下：

P_i＝|x_i|+|x_i+1-x_i|，0≤i≤N-2 (1)

其中，{x_i}为地震记录的离散序列，N为序列长度，i为每个地震记录的离散序列编号；

步骤3：确定长、短窗长度后，按确定的步长同步移动长、短窗，同时实时对特征函数计算STA/LTA数值，如图3(2)对应的长短时均值比曲线所示，具体公式如下：

其中，{P_i}为对序列Pi进行长短时均值比计算得到序列，N₄为长短时均值比方法的短尺度，N₃为长尺度；

步骤4：判断STA/LTA是否超过阈值P_c，若是则取其对应的采样时间为地震波到时，同时触发地震预警装置；否则返回到步骤1，具体公式如下：

P_n＝max{Pi}≥P_c (3)

阈值P_c的确定应考虑环境噪音、仪器噪音和地方震特点。

现有技术中振幅变化长短时均值比方法存在的问题：

拾取P波的特征函数，具体公式如下：

CF_P＝x_ud(t)²+[x_ud(t)-x_ud(t-1)]² (4)

公式中x_ud(t)是t时刻的地震记录。

目前方法存在的问题是：

此公式中用“振幅变化平方”作为特征函数拾取P波，由于P波初至时振幅都比较低，以cms^-2为单位小于1，甚至远远小于1，平方后数值变小，从而使异常降低，即用“振幅变化平方”作为特征函数拾取P波，捕捉地震波到时能力不强的问题；用“振幅变化平方”作为特征函数拾取P波，使计算相对繁琐，地震预警效率相对较低的问题；根据振幅变化拾取P波到时分析可看出，现有方法在较大振幅噪音条件下，已经不能拾取P波到时的问题；在较小振幅噪音条件下，拾取准确度较低的问题。

抗噪音分析：

地震预警的关键技术之一是P波到时自动拾取的抗噪音问题。需要把现有技术方法和改进方法的计算结果进行比较，以说明两者步骤和计算结果的区别。

分别用原有“振幅变化平方长短时均值比“方法和改进”振幅变化长短时均值比方法”拾取P波到时，对图4中某一地震记录在P波到达前2s施加振幅不同的噪音，其中(1)是在地震波到达前2s前施加的噪音信号，其振幅可以从10^-5cms^-2增加到10^-3cms^-2，(2)是采样周期为0.005s的汶川地震竖直分量加速度记录，P波初至振幅为10^-3cms^-2数量级。具体操作内容如下：

图5是分别对图4中的地震P波到达前2s施加振幅为10^-5cms^-2噪音振幅(小于P波初至振幅二个数量级)，用原方法和改进方法得到相关曲线。其中(1)和(2)分别是用原方法得到的“振幅变化平方”曲线和对应的长短时均值比曲线；(3)和(4)分别是用改进方法得到的“振幅变化”曲线和对应的长短时均值比曲线。

图6是分别对图4中的地震施加振幅为10^-4cms^-2噪音振幅(小于P波初至振幅一个数量级)，用原方法和改进方法得到相关曲线。其中(1)和(2)分别是用原方法得到的“振幅变化平方”曲线和对应的长短时均值比曲线；(3)和(4)分别是用改进方法得到的“振幅变化”曲线和对应的长短时均值比曲线。

图7是分别对图4中的地震施加振幅为10^-3cms^-2噪音振幅，此时的噪音振幅与地震P波振幅同数量级，但小于P波初至振幅。用原方法和改进方法得到的相关曲线。其中(1)和(2)分别是用原方法得到的“振幅变化平方”曲线和对应的长短时均值比曲线；(3)和(4)分别是用改进方法得到的“振幅变化”曲线和对应的长短时均值比曲线。

对比分析结果如下：

从图5可以看出在噪音振幅小的条件下，两种方法都准确地拾取了P波到时；从图6可以看出在噪音振幅增大的条件下，虽然两种方法都准确地拾取了P波到时，但新方法在P波到时的异常明显，在规定阈值的条件下，更容易准确地拾取P波；从图7可以看出在噪音振幅继续增大的条件下，原方法由于峰值未出现在P波初至位置已经不能拾取P波到时，但新方法在P波到时处的异常明显。

为进一步验证用“振幅变化长短时均值比方法”拾取P波到时抗噪音能力。选择震级大于Ms4.0的汶川大地震和余震共160条记录(采样周期为0.005s)，首先用手工拾取P波到时，以此为基准，对用原方法和新方法拾取到的P波到时进行比较，偏差小于等于0.1s认为拾取成功。拾取分二种情况进行，地震到达前噪音振幅分别为10^-5cms^-2和10^-4cms^-2。具体操作内容如下：

图8中(1)和(2)分别是噪音振幅10^-5cms^-2时，用原方法得到的P波自动拾取统计偏差点状图及统计直方图；(3)和(4)分别是噪音波振幅为10^-5cms^-2时，用新方法得到的P波自动拾取统计偏差点状图及统计直方图。

图9中(1)和(2)分别是噪音振幅10^-4cms^-2时，用原方法得到的P波自动拾取统计偏差点状图及统计直方图；(3)和(4)分别是噪音波振幅为10^-4cms^-2时，用新方法得到的P波自动拾取统计偏差点状图及统计直方图。

对比分析结果如下：

从图8可以看出两种方法拾取P波到时的偏差都小于等于0.1s，也就是都正确触发了地震预警，但用新方法拾取的偏差更小些；从图9可以看出两种方法拾取P波到时偏差大多数都小于0.1s，但原方法拾取的P波到时偏差大于新方法拾取的偏差。

Claims (1)

1.一种用振幅变化长短时均值比实现地震预警中p波自动拾取的方法，其特征在于，包括如下流程：

步骤1：实时输入竖直分量地震加速度记录流；

步骤2：计算地震数据的特征函数，具体公式如下：

P_i＝|x_i|+|x_i+1-x_i|，0≤i≤N-2

其中，{x_i}为地震记录的离散序列，N为序列长度，i为每个地震记录的离散序列编号；

步骤3：确定长、短窗长度后，按确定的步长同步移动长、短窗，同时实时对特征函数计算STA/LTA数值，具体公式如下：

其中，{P_i}为对序列Pi进行长短时均值比计算得到序列，N₄为长短时均值比方法的短尺度，N₃为长尺度；

步骤4：判断STA/LTA是否超过阈值P_c，具体公式如下：

P_n＝max{Pi}≥P_c

若是，则取其对应的采样时间为地震波到时，同时触发地震预警装置；否则返回到步骤1。