CN112068195B - 微地震p＆s波匹配事件初至自动拾取方法和计算机存储介质 - Google Patents

Abstract

本发明公开了一种微地震P&S波匹配事件初至自动拾取方法和计算机存储介质。该方法通过两两分量偏振分析与旋转处理，从原始微地震三分量矢量波场中分离出P波数据和S波数据以及S波剩余分量T；然后，利用长短时窗能量比、改进能量比，构建P波初至拾取非线性特征曲线，设定门槛值，自动搜索出P波初至时间；最后，利用S波分量与S波剩余分量T，构建S波协方差矩阵，计算最大特征值，形成以S波特征值为中心的各种能量比曲线，构建S波初至拾取非线性特征曲线，设置P波和相匹配的S波时间范围，自动搜索出区域峰值对应的时间，即为与P波相匹配的S波初至时间，从而实现微地震P&S波初至时间自动拾取。

Description

微地震P＆S波匹配事件初至自动拾取方法和计算机存储介质

技术领域

本发明涉及微地震监测数据处理技术领域，尤其涉及一种微地震P&S波匹配事件初至自动拾取方法和计算机存储介质。

背景技术

微地震监测数据处理重点在于事件的准确定位，影响因素主要包括速度模型建立、反演算法的适用性、正演算法的精度、初至拾取等方面。其中，纵、横波初至准确拾取是震源精确定位首要条件之一。井中微地震初至拾取方法，一般采用常规地震直达波初至拾取方法，其原理主要是根据有效波与噪声在能量、偏振特性以及其他一些统计特性上存在区别，获得稳定可靠的初至时间，如能量分析法、自回归AR模型法、偏振分析法等等。

每种类型方法都有各自算法特点，能量分析法是基于长短时窗能量比，当信号到达时,能量比变化快,相应的值会有一个明显的突跳，将该点时间定义为有效事件初至时间，但是信噪比较低情况容易出现误拾、漏拾；自回归AR模型法是基于信号与背景噪音属于不同AR模型，信号到达时AIC值极小值，但是不能直接判断是否是有效信号；偏振分析法是基于有效信号偏振度高、随机信号偏振度低，对应不同线性偏振系数曲线，但是该方法无法单独进行有效信号的检测。

由于在水力压裂裂缝过程中，岩石破碎产生的微地震信号，既有弹性波又有剪切波，导致监测到的微地震事件震相类型复杂，可能是P波、S波组合，也可能只产生单一的P波或S波，并且信号能量不等，而且受到地层因素、井周边环境影响及各类噪音的干扰，微地震初至相位容易受到影响波动，使得事件初至时间更难拾取，这直接影响微地震事件定位准确性。另外，从微地震定位精度上，对于同一微地震事件，P&S波联合定位精度往往比单一P波或S波定位精度高，而P&S波联合定位需要准确拾取P&S初至时间。因此，如何从复杂的微地震监测数据记录中识别出更加准确、稳定的相匹配的微地震事件P&S波初至时间，是亟须解决的技术问题。

发明内容

针对上述现有技术的不足，本发明提供了一种微地震P&S波匹配事件初至自动拾取方法和计算机存储介质。

本发明的微地震P&S波匹配事件初至自动拾取方法，主要包括以下步骤：

S100，针对微地震P波信号，从原始微地震三分量矢量波场中分解出沿波传播方向P波分量；

S200，针对微地震S波信号，从原始微地震三分量矢量波场中分解出垂直传播方向S波分量和切向S波剩余分量T；

S300，根据P波分量确定P波长短时窗能量比和改进能量比，进而确定P波能量比特征系数，根据P波长短时窗能量比和P波能量比特征系数构建P波初至拾取非线性特征曲线；

S400，从P波初至拾取非线性特征曲线上搜索出大于给定的P波初至识别门槛值的峰值，拾取所述峰值对应的时间点作为P波初至时间；

S500，根据S波分量与切向S波剩余分量T构建S波协方差矩阵，计算其最大特征值，基于最大特征值计算S波特征值前后时窗能量比和改进能量比，进而确定S波能量比特征系数，根据S波特征值前后时窗能量比和S波能量比特征系数构建S波初至拾取非线性特征曲线；

S600，确定与P波初至时间相匹配的S波初至时间的可能的时间范围，在此时间范围内，自动搜索出S波初至拾取非线性特征曲线上的极值，拾取所述极值对应的时间点作为与P波初至时间相匹配的S波初至时间；

S700，输出P波初至时间和S波初至时间，从而完成微地震P&S波匹配事件初至自动拾取。

根据本发明的实施例，上述步骤S100包括以下步骤：

S110，针对微地震P波信号，通过矢端曲线分析，对原始微地震三分量中的X分量和Y分量进行偏振分析与旋转处理，获得波传播径向分量Rp和切向分量Tp；

S120，再通过矢端曲线分析，对波传播径向分量Rp和原始微地震三分量中的Z分量进行偏振分析与旋转处理，获得沿波传播方向P波分量和垂直传播方向N分量。

其中，首先，针对P波信号，对原始微地震三分量X、Y分量进行偏振分析，获得P波水平方位角α，并作旋转处理，获得P波新的水平两分量，即旋转后的P波水平径向分量、旋转后的P波水平切向分量：

R_Pi＝X_i cos(α)+Y_i sin(α)

T_Pi＝-X_i sin(α)+Y_i cos(α)

其中i为时间样点，X_i、Y_i分别为原始X、Y分量瞬时振幅，R_Pi、T_Pi分别为旋转后P波水平径向分量、P波水平切向分量瞬时振幅。

然后，再次针对P波信号，对原始微地震三分量Z分量和上式中的新的R分量进行偏振分析，获得垂向偏振角β，并进一步作旋转处理，得到沿传播方向P波分量与垂直传播方向N分量：

P_i＝Z_i cos(β)+R_Pi sin(β)

N_Pi＝-Z_i sin(β)+R_Pi cos(β)

其中，Z_i为原始Z分量瞬时振幅，P_i、N_Pi分别为旋转后沿传播方向P波分量、垂直传播方向N分量瞬时振幅。

最终目的是通过以上两两垂直分量偏振分析与旋转，从原始微地震三分量X、Y、Z分量矢量求和，获得P波信号最强的、沿传播方向的P波分量。

根据本发明的实施例，上述步骤S200包括以下步骤：

S210，针对微地震S波信号，通过矢端曲线分析，对原始微地震三分量中的X分量和Y分量进行偏振分析与旋转处理，获得波传播径向分量Rs和切向分量Ts；

S220，再通过矢端曲线分析，对波传播的S波水平径向分量R_S和原始微地震三分量中的Z分量进行偏振分析与旋转处理，获得垂直传播方向S波分量和切向S波剩余分量T。

其中，首先，针对S波信号，对原始微地震三分量X、Y分量进行偏振分析，获得S波水平方位角α^*，并作旋转处理，获得S波相应的新的水平两分量，即旋转后的S波径向分量、旋转后的S波切向分量：

R_Si＝X_i cos(α^*)+Y_i sin(α^*)

T_Si＝-X_i sin(α^*)+Y_i cos(α^*)

其中i为时间样点，X_i、Y_i分别为原始X、Y分量瞬时振幅，R_Si、T_Si分别为旋转后S波径向分量和旋转后的S波切向分量瞬时振幅。

然后，再次针对S波信号，对原始微地震三分量Z分量和上式中的S波新的R分量进行偏振分析，获得S波垂向偏振角β^*，并进一步作旋转处理，得到垂直于沿传播方向S波分量与剩余垂向N分量：

S_i＝Z_i cos(β^*)+R_Si sin(β^*)

N_Si＝-Z_i sin(β^*)+R_Si cos(β^*)

其中，Z_i为原始Z分量瞬时振幅，S_i、N_Si分别为旋转后沿垂直传播方向的S波分量瞬时振幅、剩余垂向N分量瞬时振幅。

虽然S波偏振特性不如P波，但是通过类似于P波两两垂直分量偏振分析与旋转，最终目的是可以从原始微地震三分量X、Y、Z分量中获得信号较强的、垂直于传播方向的S波分量以及S波剩余水平分量T。

根据本发明的实施例，上述步骤S300包括以下步骤：

S310，根据P波分量按照下式计算P波长短时窗能量比：

其中，ER_Pi为时间样点i对应的P波长短时窗能量比，P_i为时间样点i对应的P波分量，L₁和L₂分别为给定的长时窗和短时窗；

S320，根据P波长短时窗能量比ER_Pi，通过引入P波分量的振幅作为权重系数，按照下式计算一级改进能量比MER1_Pi和二级改进能量比MER2_Pi：

MER1_Pi＝ER_Pi·|P_i|

MER2_Pi＝ER_Pi ²·|P_i|²。

其中，|P_i|为时间样点i对应的P波分量的振幅；

S330，根据P波长短时窗能量比ER_Pi、一级改进能量比MER1_Pi和二级改进能量比MER2_Pi，按照下式计算P波能量比特征系数AA_Pi、BB_Pi、CC_Pi：

S340，根据所述P波长短时窗能量比ER_Pi和P波能量比特征系数AA_Pi、BB_Pi、CC_Pi，构建以下P波初至拾取特征曲线TER_Pi：

TER_Pi＝ER_Pi*(AA_Pi+BB² _Pi+CC³ _Pi)。

根据本发明的实施例，上述步骤400包括以下步骤：

S410，设置P波初至识别门槛值K_P，搜索所述P波初至拾取特征曲线TER_Pi上大于P波初至识别门槛值K_P时对应的时间样点的位置T_KP，其中：

当i<T_KP时TER_Pi<K_P，

当i＝T_KP时TER_Pi≥K_P；

S420，设置初至搜索时窗的大小W_P，搜索在初至搜索时窗[T_KP,T_KP+W_P]内P波初至拾取特征曲线TER_Pi的峰值F_P：

F_P＝max{TER_Pi}且i∈[T_KP,T_KP+W_P]

S430，将搜索出的峰值F_P所对应的时间样点t_P作为P波初至时间。

根据本发明的实施例，上述步骤500包括以下步骤：

S510，根据S波分量和切向S波剩余分量T，构建以下S波瞬时协方差矩阵CC_Si：

其中，S_i、T_Si分别为样点i处的S波分量瞬时振幅、切向S波剩余分量T瞬时振幅，∑表示对样点时窗内数值求和；

S520，对S波协方差矩阵CC_Si进行奇异值SVD分解，计算出其最大特征值λ_Si：

CC_Si＝{λ₁,λ₂}U_SU_S ^-1

λ_Si＝max{λ₁,λ₂}

其中，{λ₁,λ₂}、U_S分别为S波特征值、特征向量。

S530，按照下式计算S波特征值前后时窗能量比,：

其中，EER_Si为时间样点i的S波特征值前后时窗能量比，L为给定的时窗。

S540，将S波最大特征值作为权重系数引入，按照下式计算一级改进能量比MEER1_Si和二级改进能量比MEER2_Si：

MEER1_Si＝EER_Si·λ_Si

MEER2_Si＝EER² _Si·λ_Si ²

S550，通过S波特征值前后时窗能量比EER_Si、一级改进能量比MEER1_Si和二级改进能量比MEER2_Si，按照下式计算S波能量比特征系数AA_Si、BB_Si、CC_Si：

S560，基于S波特征值前后时窗能量比EER_Si、S波能量比特征系数AA_Si、BB_Si、CC_Si，构建以下S波初至自动拾取非线性特征曲线TER_Si：

TER_Si＝EER_Si*(AA_Si+BB² _Si+CC³ _Si)。

根据本发明的实施例，上述步骤S600包括以下步骤：

S610，根据已有微地震观测系统和测井数据，定义与P波初至时间t_P相匹配的S波初至时间的可能的时间范围[t_P+Δt₁,t_P+Δt₂]；

S620，在时间范围[t_P+Δt₁,t_P+Δt₂]内，从S波特征值前后时窗能量比曲线上搜索出极值F_S：

F_S＝max{EER_Si}且i∈[t_P+Δt₁,t_P+Δt₂]

S630，将搜索出的极值F_S对应的时间样点t_S定义为与P波初至时间t_P相匹配的S波初至时间。

根据本发明的实施例，在上述步骤S410中，所述P波初至识别门槛值K_P是通过基于已知射孔信号识别进行试算来确定：根据本发明P波初至拾取特征曲线TER_Pi，对已知射孔数据进行微地震信号识别，将门槛值K_P(大于0)从较小数值逐步增大进行信号识别试算，一开始会出现很多个识别结果，但是由于射孔信号是已知的，也就是说，射孔只有一个P波信号且初至是已知的，即当门槛值K_P达到某个值，刚好识别出一个P波初至且时间位置刚好或接近真实位置，此时门槛值K_P数值则为最优门槛值，被应用于其他事件识别。

根据本发明的实施例，上述时窗L₁和L₂以及L根据选择已知射孔或者人为选择一个事件进行试算来确定。不同地区数据有不同经验值，长时窗L₁与短时窗L₂比值常常选择在3到6之间(如L₁＝120，L₂＝20)，时窗L常常取值20到100之间，要求计算出的初至拾取特征曲线尽可能突出有效信号、同时压制周围干扰，目的是便于微地震事件信号识别。

此外，本发明还提供一种计算机存储介质，其中存储有用于实现上述方法的计算机程序。

与现有技术相比，本发明的一个或多个实施例可以具有如下优点：

本发明提供了一种同一微地震事件P波、S波初至联合自动拾取方法，首先，通过两两分量偏振分析与旋转处理，从原始微地震三分量矢量波场中分离出P波数据和S波分量以及切向S波剩余分量T；然后，利用长短时窗能量比、改进能量比，构建P波初至拾取非线性特征曲线，设定门槛值，自动搜索出P波初至时间；最后，利用S波分量以及切向S波剩余分量T，构建S波协方差矩阵，计算最大特征值，形成以S波特征值为中心的各种能量比曲线，构建S波初至拾取非线性特征曲线，设置P波和相匹配的S波时间范围，自动搜索出区域峰值对应的时间，即为与P波相匹配的S波初至时间，至此实现微地震P&S波初至时间自动拾取。本发明能够实现对微地震事件快速、准确地定位，整个定位过程简单便捷。

本发明的其它特征和优点将在随后的说明书中阐述，并且部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1是本发明P&S波初至自动拾取操作流程图；

图2是微地震事件三分量数据：(a)Z分量；(b)X分量；(c)Y分量；

图3是微地震事件三分量数据矢量分解出沿传播方向P波分量；

图4是微地震事件三分量数据矢量分解出垂直传播方向S波分量与其剩余分量；

图5是微地震事件P波分量长短时窗能量比曲线(a)、一级改进能量比(b)、二级改进能量比(c)；

图6是微地震事件P波分量特征系数：(a)图5(a)特征系数；(b)图5(b)特征系数；(c)图5(c)特征系数；

图7是微地震事件P波初至自动拾取非线性特征曲线；

图8是微地震事件S波最大特征值对应的时窗能量比曲线(a)、一级改进能量比(b)、二级改进能量比(c)；；

图9是微地震事件S波最大特征值对应的特征系数：(a)图8(a)特征系数；(b)图8(b)特征系数；(c)图8(c)特征系数；

图10是微地震事件S波初至自动拾取非线性特征曲线；

图11是识别出的P&S波初至时间在微地震事件P波、S波分量上箭头标识；

具体实施方式

微地震事件P&S波初至时间准确拾取是微地震事件P&S联合定位前提条件之一，也是微地震三分量处理关键技术之一。本发明通过两两分量偏振分析与旋转处理，从原始微地震三分量矢量波场中分离出P波数据和S波数据以及切向S波剩余分量T；然后，利用长短时窗能量比、改进能量比，构建P波初至拾取非线性特征曲线，设定门槛值，自动搜索出P波初至时间；最后，利用S波分量与切向S波剩余分量T，构建S波协方差矩阵，计算最大特征值，形成以S波特征值为中心的各种能量比曲线，构建S波初至拾取非线性特征曲线，设置P波和相匹配的S波时间范围，自动搜索出区域峰值对应的时间，即为与P波相匹配的S波初至时间，以此来实现微地震P&S波初至时间自动拾取。

下面结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

实施例一

本发明提供了一种基于矢量波场分离的微地震P&S波匹配事件初至自动拾取方法。如图1所示，在本实施例中，该方法包含以下五个步骤：

第一步，针对P波信号，通过矢端曲线分析，对原始微地震三分量中的X、Y分量进行偏振分析与旋转处理，分别获得波传播径向分量Rs和切向分量Ts与P波切向分量Tp，再对S波径向分量Rs与原始Z分量数据进行偏振分析与旋转处理，分别获得沿波传播方向P波分量、垂直传播方向N分量，同样地，针对S波信号，通过两两分量偏振分析与旋转，可以获得垂直传播方向S波分量和切向S波剩余分量T；

第二步，输入第一步矢量波场分离出的P波分量数据，计算长短时窗能量比ER_Pi，同时将P波振幅作为权系数引入，计算一级改进能量比MER1_Pi＝ER_Pi·|P_i|、二级改进能量比MER2_Pi＝ER² _Pi·P_i ²，再进一步获得对应的P波能量比特征系数

第三步，输入第二步长短时窗能量比ER_Pi与三个P波能量比特征系数，构建P波初至自动拾取非线性特征曲线TER_Pi＝ER_Pi*(AA_Pi+BB² _Pi+CC³ _Pi)，根据门槛值，自动搜索特征曲线TER_Pi大于门槛值极值，其对应的时间，即为所求事件P波初至时间t_P；

第四步，将第一步矢量波场分离出的S波分量与切向S波剩余分量T，构建S波瞬时协方差矩阵ST_Si，计算其最大特征值λ_Si，在此基础上进一步计算特征值相应的瞬时能量比EER_Si，计算S波特征值一级改进能量比MEER1_Si＝EER_Si*λ_Si、二级改进能量比MEER2_Si＝EER² _Si*λ² _Si，以及对应的S波能量比特征系数

第五步，输入第四步S波特征值能量比EER_Si与三个S波特征值能量比特征系数，构建S波初至自动拾取非线性特征曲线TER_Si＝EER_Si*(AA_Si+BB² _Si+CC³ _Si)，再根据第二步识别出的微地震事件P波初至时间t_P，设计同一微地震事件与P波相匹配的S波初至可能时间范围[t_P+Δt₁,t_P+Δt₂]，搜索该范围内S波非线性特征曲线TER_Si上区域峰值对应时间，即为与P波相匹配S波初至时间t_S，结合识别出的P波初至时间t_P，实现本发明微地震P&S波匹配事件初至自动拾取。

下面具体介绍各步骤。

首先，本实施例优选采用已知算法--矢端曲线分析法，针对P波信号，对原始微地震三分量X、Y分量进行偏振分析，获得水平方位角α。然后，根据旋转式，对原始X、Y分量数据进行旋转处理，获得新的两分量S波径向分量Rs和切向分量Ts：

其中，i为时间样点，X_i、Y_i分别为原始X、Y分量瞬时振幅，R_Pi、T_Pi分别为旋转后S波径向分量和旋转后S波切向分量瞬时振幅。

以同样的操作，针对P波信号，对原始微地震三分量Z分量和新的R分量进行偏振分析，获得垂向偏振角β，并进一步作旋转处理，得到沿传播方向P波分量与垂直传播方向N分量：

其中，Z_i为原始Z分量瞬时振幅，P_i、N_Pi分别为旋转后沿传播方向P波分量、垂直传播方向N分量瞬时振幅。

同样地，针对S波信号，通过两两分量偏振分析与旋转，可以获得垂直传播方向S波分量和切向S波剩余分量T，对应的瞬时振幅分别为S_i、T_Si。

然后，本发明通过长短时窗能量比、改进能量比构建初至识别非线性特征曲线来识别出P波初至时间。定义长时窗L₁、短时窗L₂，对P波分量数据进行长短时窗能量比计算：

其中，ER_Pi为时间样点i对应的P波长短时窗能量比。将P波分量振幅作为权系数引入，计算一级改进能量比MER1_Pi、二级改进能量比MER2_Pi：

MER1_Pi＝ER_Pi·|P_i| (4)

MER2_Pi＝ER_Pi ²·P_i ² (5)

通过上述三个能量比ER_Pi、MER1_Pi、MER2_Pi，分别计算出P波能量比特征系数AA_Pi、BB_Pi、CC_Pi：

构建本发明P波初至自动拾取非线性特征曲线TER_Pi：

TER_Pi＝ER_Pi*(AA_Pi+BB² _Pi+CC³ _Pi) (9)

设定门槛值K_P，搜索P波特征曲线TER_Pi大于门槛值时对应的时间样点位置T_KP，即出现事件判断标准：当i<T_KP且TER_Pi<K_P；当i＝T_KP且TER_Pi≥K_P。

同时，定义初至搜索时窗大小W_P，找出时窗内最大峰值F_P：F_P＝max{TER_Pi}且i∈[T_KP,T_KP+W_P]。将搜索出的最大峰值F_P对应的时间样点t_P输出，即实现微地震事件P波初至时间自动拾取。

其次，本发明通过构建协方差矩阵、奇异值分解、时窗能量比来构建S波特征值能量比特征系数AA_Si、BB_Si、CC_Si。定义样点时窗大小，由矢量分解出的S波分量和切向S波剩余分量T，构建S波瞬时协方差矩阵：

其中，S_i、T_Si分别为样点i处的S波分量瞬时振幅、切向S波剩余分量T瞬时振幅，∑表示对样点时窗内数值求和。

对S波协方差矩阵CC_Si奇异值SVD分解，计算出相应最大特征值λ_Si：CC_Si＝{λ₁,λ₂}U_SU_S ^-1且λ_Si＝max{λ₁,λ₂}。其中，{λ₁,λ₂}、U_S分别为S波特征值、特征向量。

定义时窗L，对所有S波最大特征值λ_Si数据进行能量比计算：

其中，EER_Si为样点i的S波最大特征值前后时窗能量比。同时，将S波最大特征值作为权系数引入，计算与其相关的一级改进能量比MEER1_Si、二级改进能量比MEER2_Si：

MEER1_Si＝EER_Si·λ (12)

MEER2_Si＝EER² _Si·λ² (13)

通过S波特征值瞬时能量比EER_Si、一级改进能量比MEER1_Si、二级改进能量比MEER2_Si，构建S波能量比特征系数AA_Si、BB_Si、CC_Si：

最后，构建S波初至识别非线性特征曲线，自动识别出与其P波初至相匹配的波初至时间。通过S波特征值瞬时能量比EER_Si及其特征系数AA_Si、BB_Si、CC_Si，构建本发明S波初至自动拾取非线性特征曲线TER_Si：

TER_Si＝EER_Si*(AA_Si+BB² _Si+CC³ _Si) (17)

根据已有微地震观测系统和测井数据，定义同一微地震事件与识别出的P波初至t_P相匹配的S波初至可能出现的时间范围[t_P+Δt₁,t_P+Δt₂]，搜索出该时间段内TER_Si数据最大峰值F_S：

F_S＝max{TER_Si}且i∈[t_P+Δt₁,t_P+Δt₂] (18)

其中，最大峰值F_S对应的时间样点t_S，定义为与P波相匹配S波初至时间t_S，结合识别出的P波初至时间t_P，即最终实现本发明微地震P&S波匹配事件初至自动拾取。

实施例二

下面从微地震实际事件资料应用中说明本发明P&S波初至自动拾取过程。

图2为本实施例的实际微地震事件三分量单道数据，它包含多个P波、S波微地震信号。本实施例的主要目的是根据图1所示的工作流程，准确地拾取P波初至以及相匹配的S博初至：

首先，根据图1所示的工作流程，对如图2所示的原始微地震事件三分量数据开展矢量分解，即通过式(1)、(2)，进行两两分量偏振分析和旋转处理，获得事件沿传播方向、达到三分量检波器之前的P波分量，如图3所示，在能量上，它是原始三分量P波矢量求和，并且主要包含P波信号、基本无S波。同样操作，针对S波，开展两两分量偏振分析与旋转处理，矢量分解出垂直传播方向S波和切向S波剩余分量T，如图4所示。其中，S分量主要包含S波信号，剩余分量T则主要是噪音以及极少量由于S波不具有严格偏振关系而残余微地震信号。

然后，针对P波分量，根据式(3)～(5)，分别计算长短时窗能量比、一级改进能量比、二级改进能量比(如图5和图6所示)，并根据式(6)～(8)，计算出P波分量三个特征系数(如图2所示)，进而根据式(9)构建出P波初至拾取非线性特征曲线(如图7所示)，给定门6槛值(如门槛值为“25”)，可以识别出5个事件P波初至。

其次，针对S波分量，根据式(10)，构建S波协方差矩阵，计算S波最大特征值，目的是为了提高S波初至拾取准确性，用最大特征值代替原始微地震数据，计算出S波基于特征值的时窗能量比、一级改进能量比、二级改进能量比(式(11)～(13)、如图8所示)，再根据式(14)～(16)，计算出相应的S波基于特征值的能量比特征系数(如图9所示)，进而根据式(17)，构建出高信噪比的S波初至拾取的非线性特征曲线(如图10所示)。

最后，针对S波初至拾取的非线性特征曲线，不是通过门槛值拾取S波初至，而是通过与P波匹配关系自动搜索出S波初至。该过程需根据已有微地震观测系统和测井数据，定义同一微地震事件与识别出的P波初至t_P相匹配的S波初至可能出现的时间范围[t_P+Δt₁,t_P+Δt₂]，搜索出该时间段内S波初至拾取特征曲线数据最大峰值F_S(式(18))。其最大峰值F_S对应的时间样点t_S，即为与P波相匹配S波初至时间t_S，结合识别出的P波初至时间t_P，最终实现本发明微地震P&S波匹配事件初至自动拾取。将这5对P&S波初至时间标识在P波、S波分量上，如图11所示箭头，都是真实的P波、S波信号显示。

上述实施例较好地验证了本发明的P&S波初至自动拾取方法的准确性，具有较好的实际应用价值。

实施例三

此外，本发明还提供一种计算机存储介质，其特征在于，其中存储有用于实现上述方法的计算机程序,在此不再赘述。

应该理解的是，本发明所公开的实施例不限于这里所公开的特定处理步骤或材料，而应当延伸到相关领域的普通技术人员所理解的这些特征的等同替代。还应当理解的是，在此使用的术语仅用于描述特定实施例的目的，而并不意味着限制。

说明书中提到的“实施例”意指结合实施例描述的特定特征、或特性包括在本发明的至少一个实施例中。因此，说明书通篇各个地方出现的短语“实施例”并不一定均指同一个实施例。

本领域的技术人员应该明白，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。本领域的技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储介质(RAM)、内存、只读存储介质(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM或技术领域内所公知的任意其它形式的存储介质中。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (8)

1.一种微地震P&S波匹配事件初至自动拾取方法，其特征在于，包括以下步骤：

S100，针对微地震P波信号，从原始微地震三分量矢量波场中分解出沿波传播方向P波分量；

S200，针对微地震S波信号，从原始微地震三分量矢量波场中分解出垂直传播方向S波分量和切向S波剩余分量T；

S300，根据P波分量确定P波长短时窗能量比和改进能量比，进而确定P波能量比特征系数，根据P波长短时窗能量比和P波能量比特征系数构建P波初至拾取非线性特征曲线；

S400，从P波初至拾取非线性特征曲线上搜索出大于给定的P波初至识别门槛值的峰值，拾取所述峰值对应的时间点作为P波初至时间；

S500，根据S波分量与切向S波剩余分量T构建S波协方差矩阵，计算其最大特征值，基于最大特征值计算S波特征值前后时窗能量比和改进能量比，进而确定S波能量比特征系数，根据S波特征值前后时窗能量比和S波能量比特征系数构建S波初至拾取非线性特征曲线；

S600，确定与P波初至时间相匹配的S波初至时间的可能的时间范围，在此时间范围内，自动搜索出S波初至拾取非线性特征曲线上的极值，拾取所述极值对应的时间点作为与P波初至时间相匹配的S波初至时间；

S700，输出P波初至时间和S波初至时间，从而完成微地震P&S波匹配事件初至自动拾取；

所述步骤S300包括以下步骤：

S310，根据P波分量按照下式计算P波长短时窗能量比：

其中，ER_Pi为时间样点i对应的P波长短时窗能量比，P_i为时间样点i对应的P波分量，L₁和L₂分别为给定的长时窗和短时窗；

S320，根据P波长短时窗能量比ER_Pi，通过引入P波分量的振幅作为权重系数，按照下式计算一级改进能量比MER1_Pi和二级改进能量比MER2_Pi：

MER1_Pi＝ER_Pi·|P_i|

MER2_Pi＝ER_Pi ²·|P_i|²

其中，|P_i|为时间样点i对应的P波分量的振幅；

S330，根据P波长短时窗能量比ER_Pi、一级改进能量比MER1_Pi和二级改进能量比MER2_Pi，按照下式计算P波能量比特征系数AA_Pi、BB_Pi、CC_Pi：

S340，根据所述P波长短时窗能量比ER_Pi和P波能量比特征系数AA_Pi、BB_Pi、CC_Pi，构建以下P波初至拾取特征曲线TER_Pi：

TER_Pi＝ER_Pi*(AA_Pi+BB² _Pi+CC³ _Pi)；

所述步骤S500包括以下步骤：

S510，根据S波分量和切向S波剩余分量T，构建以下S波瞬时协方差矩阵CC_Si：

其中，S_i、T_Si分别为样点i处的S波分量瞬时振幅、切向S波剩余分量T瞬时振幅，∑表示对样点时窗内数值求和；

S520，对S波协方差矩阵CC_Si进行奇异值SVD分解，计算出其最大特征值λ_Si：

CC_Si＝{λ₁,λ₂}U_SU_S ^-1

λ_Si＝max{λ₁,λ₂}

其中，{λ₁,λ₂}、U_S分别为S波特征值、特征向量，

S530，按照下式计算S波特征值前后时窗能量比：

其中，EER_Si为时间样点i的S波特征值前后时窗能量比，L为给定的时窗，

S540，将S波最大特征值作为权重系数引入，按照下式计算一级改进能量比MEER1_Si和二级改进能量比MEER2_Si：

MEER1_Si＝EER_Si·λ_Si

MEER2_Si＝EER² _Si·λ_Si ²

S550，通过S波特征值前后时窗能量比EER_Si、一级改进能量比MEER1_Si和二级改进能量比MEER2_Si，按照下式计算S波能量比特征系数AA_Si、BB_Si、CC_Si：

S560，基于S波特征值前后时窗能量比EER_Si、S波能量比特征系数AA_Si、BB_Si、CC_Si，构建以下S波初至自动拾取非线性特征曲线TER_Si：

TER_Si＝EER_Si*(AA_Si+BB² _Si+CC³ _Si)。

2.如权利要求1所述的微地震P&S波匹配事件初至自动拾取方法，其特征在于，所述步骤S100包括以下步骤：

S110，针对微地震P波信号，通过矢端曲线分析，对原始微地震三分量中的X分量和Y分量进行偏振分析与旋转处理，获得波传播径向分量Rp和切向分量Tp；

S120，再通过矢端曲线分析，对波传播径向分量Rp和原始微地震三分量中的Z分量进行偏振分析与旋转处理，获得沿波传播方向P波分量和垂直传播方向N分量。

3.如权利要求1所述的微地震P&S波匹配事件初至自动拾取方法，其特征在于，所述步骤S200包括以下步骤：

S210，针对微地震S波信号，通过矢端曲线分析，对原始微地震三分量中的X分量和Y分量进行偏振分析与旋转处理，获得波传播径向分量Rs和切向分量Ts；

S220，再通过矢端曲线分析，对波传播的S波水平径向分量R_S和原始微地震三分量中的Z分量进行偏振分析与旋转处理，获得垂直传播方向S波分量和切向S波剩余分量T。

4.如权利要求1所述的微地震P&S波匹配事件初至自动拾取方法，其特征在于，所述步骤S400包括以下步骤：

S410，设置P波初至识别门槛值K_P，搜索所述P波初至拾取特征曲线TER_Pi上大于P波初至识别门槛值K_P时对应的时间样点的位置T_KP，其中：

当i<T_KP时TER_Pi<K_P，

当i＝T_KP时TER_Pi≥K_P；

S420，设置初至搜索时窗的大小W_P，搜索在初至搜索时窗[T_KP,T_KP+W_P]内P波初至拾取特征曲线TER_Pi的峰值F_P：

F_P＝max{TER_Pi}且i∈[T_KP,T_KP+W_P]

S430，将搜索出的峰值F_P所对应的时间样点t_P作为P波初至时间。

5.如权利要求1所述的微地震P&S波匹配事件初至自动拾取方法，其特征在于，所述步骤S600包括以下步骤：

S610，根据已有微地震观测系统和测井数据，定义与P波初至时间t_P相匹配的S波初至时间的可能的时间范围[t_P+Δt₁,t_P+Δt₂]；

S620，在时间范围[t_P+Δt₁,t_P+Δt₂]内，从S波特征值前后时窗能量比曲线上搜索出极值F_S：

F_S＝max{EER_Si}且i∈[t_P+Δt₁,t_P+Δt₂]

S630，将搜索出的极值F_S对应的时间样点t_S定义为与P波初至时间t_P相匹配的S波初至时间。

6.如权利要求4所述的微地震P&S波匹配事件初至自动拾取方法，其特征在于：

在所述步骤S410中，所述P波初至识别门槛值K_P通过基于已知射孔信号识别进行试算来确定。

7.如权利要求1所述的微地震P&S波匹配事件初至自动拾取方法，其特征在于：

所述时窗L₁和L₂以及L通过基于已知射孔信号识别进行试算来确定。

8.一种计算机存储介质，其特征在于，其中存储有用于实现如权利要求1至7中任意一项所述方法的程序。

Images