CN103995290B - 一种微震p波震相初至自动拾取方法 - Google Patents

Abstract

本发明提出了一种微震P波震相初至自动拾取方法，包括如下步骤：微震信号的Hilbert变换；对Hilbert包络信号的包络分析；微震P波信号震相初至的预拾取；微震P波信号震相初至的精确计算；微震P波信号震相初至自动拾取结束。本发明具有抗噪性能好、计算精度高和算法实时性较强等优点，具有很好的技术价值和应用前景。

Description

一种微震P波震相初至自动拾取方法

技术领域

本发明涉及信息处理技术领域，特别是指一种微震P波震相初至自动拾取方法。

背景技术

微震监测系统是冲击地压、煤与瓦斯突出、矿井突水等煤矿灾害预警的主要技术手段之一。微震发生时产生纵波和横波，即P波和S波，由于P波比S波传播速度快，在长距离传输时，P波与S波在图形上容易区分，而对于矿区较小范围内，S波在波形上叠加于P波尾波中不宜区分，工程上一般选择易于辨识的P波进行初至到时的标记。微震P波的初至到时拾取是震源定位及震相识别技术的关键环节，其精度直接影响震源定位精度并对震源机制解释具有重要意义。

鉴于微震监测系统拾取的震动信号具有突发性、多样性及不确定性等特点，目前的长短时均值比法(STA/LTA)、赤池信息准则(AIC)、PAI-S/K等震相初至拾取方法精度低、抗噪性能差、算法实时性不强。譬如单纯应用AIC函数求解震相初至需要在震相初至附近选择合适的时窗来计算AIC值，不同的时窗将产生不同的AIC值，时窗选择不合理时会出现错误的震相初至拾取结果。

为实现微震P波震相初至的准确拾取，本发明针对微震信号的非平稳性和非线性特点，引入Hilbert变换求解信号包络，通过设置包络信号阈值搜索震相初至的大致位置，以该位置为基础为AIC函数选择合适的计算时窗，并在选取的时窗内计算AIC函数，AIC取最小值的位置即是震相初至对应的时刻，形成一种高精度自动拾取P波震相初至的算法。

发明内容

本发明提出一种微震P波震相初至自动拾取方法，解决了现有技术中微震P波信号震相初至自动拾取精度低及算法实时性不强的问题。

本发明的技术方案是这样实现的：一种微震P波震相初至自动拾取方法，包括如下步骤：

Step 1：微震信号的Hilbert变换，

微震信号时序序列设为x(t)，t＝1，2，…，n，其中n为微震信号的采样点个数，对微震信号x(t)进行Hilbert变换后获得微震信号x(t)的Hilbert包络信号a(t)；

Step 2：对Hilbert包络信号a(t)进行包络分析；

Step 3：微震信号x(t)包括微震P波信号和微震S波信号，对微震P波信号震相初至预拾取；

Step 4：微震P波信号震相初至的精确计算；

Step 5：微震P波信号震相初至自动拾取结束。

进一步地，Hilbert变换定义为：

$H\[x\left(t\right)\]=\hat{x}\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}d\mathrm{\τ}---\left(1\right)$

微震信号x(t)的解析信号z(t)定义为：

其中，a(t)为解析信号z(t)的幅值，为解析信号z(t)的相位。

进一步地，幅值a(t)和相位分别定义为：

$a\left(t\right)=\sqrt{x^2\left(t\right)+{\hat{x}}^2\left(t\right)}---\left(3\right)$

其中，a(t)为微震信号x(t)的Hilbert包络信号。

进一步地，令y(t)＝2a(t)，得到时序序列y(t)，将y(t)在[0，1]范围内进行归一化处理，归一化处理采用公式(5)，计算得到时序序列

$\hat{y}\left(t\right)=\frac{y\left(t\right)-ymin}{y\max -ymin}---\left(5\right)$

其中，ymax，ymin分别表示时序序列y(t)中的最大值和最小值。

进一步地，设置阈值，中第一个大于阈值的数值对应的时刻t₀定义为微震P波信号震相初至的预拾取值。

优选地，阈值范围为[0.15，0.25]。

优选地，阈值为0.20。

进一步地，微震信号x(t)以t₀为基准向前及向后分别取1000点和1500点数据作为待分析数据，应用公式(6)给出的AIC函数计算AIC函数值，在AIC函数值中搜索最小值，最小值对应的时刻即为微震P波信号震相初至的精确拾取值，

AIC(k)＝klg(Var(x[1,k]))+(n-k-1)lg(Var(x[k+1,n])) (6)

其中，x(i)，i＝1，2，…，n为微震信号采样点数据，k的取值范围为[1，n]，n为采样点个数，Var(x[1，k])是指x(1)到x(k)之间k个采样点数据的方差。

本发明引入Hilbert变换求解信号包络，首先通过设置包络阈值粗略求取震相初至，再合理选择时窗应用AIC函数精确求取震相初至，从而达到微震P波震相初至高精度拾取的目的。

本发明的有益效果为：

本发明具有抗噪性能好、计算精度高和算法实时性较强等优点，具有很好的技术价值和应用前景。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明一种微震P波震相初至自动拾取方法一个实施例的流程示意图；

图2为本发明一种微震P波震相初至自动拾取方法一个实施例的微震信号波形图；

图3为本发明一种微震P波震相初至自动拾取方法一个实施例的震相初至预拾取示意图；

图4为本发明一种微震P波震相初至自动拾取方法一个实施例的震相初至AIC精确拾取示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明算法思想的描述如下：微震信号x(t)的Hilbert变换，Hilbert变换具体为微震信号x(t)和时间倒数的卷积，通过包络分析和预拾取值的参照，结合AIC函数进行数据分析，进而对微震P波震相初至精确自动拾取。

如图1所示，本发明一种微震P波震相初至自动拾取方法中算法的步骤为：

Step 1：微震信号的Hilbert变换，

微震信号时序序列设为x(t)，t＝1，2，…，n，其中n为微震信号的采样点个数，对微震信号x(t)进行Hilbert变换后获得微震信号x(t)的Hilbert包络信号a(t)；

其中Hilbert变换定义为：

$H\[x\left(t\right)\]=\hat{x}\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}d\mathrm{\τ}---\left(1\right)$

微震信号x(t)的解析信号z(t)定义为：

公式(2)中，a(t)为解析信号z(t)的幅值，为解析信号z(t)的相位。

幅值a(t)和相位分别定义为：

$a\left(t\right)=\sqrt{x^2\left(t\right)+{\hat{x}}^2\left(t\right)}---\left(3\right)$

公式(3)中，a(t)为微震信号x(t)的Hilbert包络信号。

Step 2：对Hilbert包络信号a(t)进行包络分析；

令y(t)＝2a(t)，得到时序序列y(t)，将y(t)在[0，1]范围内进行归一化处理，归一化处理采用公式(5)，计算得到时序序列

$\hat{y}\left(t\right)=\frac{y\left(t\right)-ymin}{y\max -ymin}---\left(5\right)$

公式(5)中，ymax，ymin分别表示时序序列y(t)中的最大值和最小值。

Step 3：微震信号x(t)包括微震P波信号和微震S波信号，对微震P波信号震相初至预拾取；

设置阈值，中第一个大于阈值的数值对应的时刻t₀定义为微震P波信号震相初至的预拾取值。阈值范围为[0.15，0.25]，本实施例阈值取0.20。

Step 4：微震P波信号震相初至的精确计算；

微震信号x(t)以t₀为基准向前及向后分别取1000点和1500点数据作为待分析数据，应用公式(6)给出的AIC函数计算AIC函数值，在AIC函数值中搜索最小值，最小值对应的时刻即为微震P波信号震相初至的精确拾取值，

AIC(k)＝klg(Var(x[1,k]))+(n-k-1)lg(Var(x[k+1,n])) (6)

其中，x(i)，i＝1，2，…，n为微震信号采样点数据，k的取值范围为[1，n]，n为采样点个数，Var(x[1，k])是指x(1)到x(k)之间k个采样点数据的方差。

Step 5：微震P波信号震相初至自动拾取结束。

如图2所示，Step 1以时刻t为横轴，振幅Amplitude为纵轴，微震信号时序序列表示为x(t)，t＝1，2，…，3000，微震信号采样点数据见表1。

表1微震信号采样点数据(可以存储于Excel中)

按照Step 1提供的算法编程对微震信号时序序列x(t)进行Hilbert变换，本实施例通过调用MATLAB 7.0版本中Hilbert函数求解出微震信号x(t)的Hilbert包络信号a(t)，t＝1，2，…，3000。

如图3所示，Step2按照Step1提供的方法对包络a(t)进行处理，计算得到时序序列t＝1，2，…，3000，以t为横坐标，为纵坐标。

Step 3按照Step 2将时序序列中的数据按时序与所设阈值0.20顺序进行比较，当大于0.20时，t的取值t₀即为震相初至的预拾取值，如图3所示，求得t₀＝1239ms。

如图4所示，Step 4在微震波形时序序列x(t)上，以t₀＝1239ms为基准向前及向后分别取1000点和1500点数据作为待分析数据，应用公式(6)给出的AIC函数计算AIC函数值，在AIC函数值中搜索最小值，最小值对应的时刻即为微震P波信号震相初至的精确拾取值，也就是微震P波信号震相初至时刻，即1241ms。

为实现微震P波震相初至的准确拾取，本发明引入Hilbert变换对微震信号进行包络分析，先通过设置包络阈值粗略判断震相初至所在的时间点，再以此时间点为基准向前向后各截取一段数据形成时间窗口，在时窗内应用赤池信息准则(Akaike InformationCriteria，AIC)精确求解震相初至，达到微震P波震相初至高精度自动拾取的目的。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种微震P波震相初至自动拾取方法，其特征在于，包括如下步骤：

Step1：微震信号的Hilbert变换，

所述微震信号时序序列设为x(t)，t＝1，2，…，n，其中n为所述微震信号的采样点个数，对所述微震信号x(t)进行所述Hilbert变换后获得所述微震信号x(t)的Hilbert包络信号a(t)，

所述Hilbert变换定义为：

$H\[x\left(t\right)\]=\hat{x}\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}d\mathrm{\τ}---\left(1\right)$

所述微震信号x(t)的解析信号z(t)定义为：

其中，a(t)为所述解析信号z(t)的幅值，为所述解析信号z(t)的相位，

所述幅值a(t)和所述相位分别定义为：

$a\left(t\right)=\sqrt{x^2\left(t\right)+{\hat{x}}^2\left(t\right)}---\left(3\right)$

其中，a(t)为所述微震信号x(t)的Hilbert包络信号；

Step2：对所述Hilbert包络信号a(t)进行包络分析；

令y(t)＝2a(t)，得到时序序列y(t)，将所述y(t)在[0，1]范围内进行归一化处理，所述归一化处理采用公式(5)，计算得到时序序列

$\hat{y}\left(t\right)=\frac{y\left(t\right)-ymin}{y\max -ymin}---\left(5\right)$

其中，ymax，ymin分别表示时序序列y(t)中的最大值和最小值；

Step3：所述微震信号x(t)包括微震P波信号和微震S波信号，对所述微震P波信号震相初至预拾取，

设置阈值，所述中第一个大于所述阈值的数值对应的时刻t₀定义为所述微震P波信号震相初至的预拾取值；

Step4：所述微震P波信号震相初至的精确计算，

所述微震信号x(t)以所述t₀为基准向前及向后分别取1000点和1500点数据作为待分析数据，应用公式(6)给出的AIC函数计算AIC函数值，在所述AIC函数值中搜索最小值，所述最小值对应的时刻即为所述微震P波信号震相初至的精确拾取值，

AIC(k)＝klg(Var(x[1,k]))+(n-k-1)lg(Var(x[k+1,n])) (6)

其中，x(i)，i＝1，2，…，n为所述微震信号采样点数据，k的取值范围为[1，n]，n为采样点个数，Var(x[1，k])是指x(1)到x(k)之间k个采样点数据的方差；

Step5：所述微震P波信号震相初至自动拾取结束。

2.根据权利要求1所述的一种微震P波震相初至自动拾取方法，其特征在于，所述阈值范围为[0.15，0.25]。

3.根据权利要求2所述的一种微震P波震相初至自动拾取方法，其特征在于，所述阈值为0.20。