CN106154332A - 一种井中微地震纵横波事件初至识别方法 - Google Patents

Abstract

本发明提供了一种井中微地震纵横波事件初至识别方法，属于井中微地震三分量信号处理领域。所述方法利用两次矢段曲线偏振分析与旋转处理，完成井中微地震三分量矢量波场分离，然后构造协方差矩阵，再根据长短时窗能量比以及引入权系数，获得事件特征曲线，最后根据门槛值搜索识别出P波、S波事件初至到达时间位置。

Description

一种井中微地震纵横波事件初至识别方法

技术领域

本发明属于井中微地震三分量信号处理领域，具体涉及一种井中微地震纵横波事件初至识别方法。

背景技术

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

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

由于在水力压裂裂缝过程中，产生的井中微地震事件震相类型是随机的，可能是P波、S波组合，也可能只产生单一的P波或S波，而且受到地层因素、井周边环境影响及各类噪音的干扰，微地震初至相位容易受到影响波动，使得事件初至时间更难拾取。

井中微地震事件初至时间准确识别是井中P/S波震源精确定位前提。目前常用的长短时窗能量比法容易受背景噪音干扰，使得相位初至拾取不稳定，出现了同一事件相邻检波器初至拾取时间值波动较大，特别是井中S波微地震事件。因此，如何从复杂的井中微地震三分量记录中识别分离出更加准确、稳定的有效微地震事件初至，是本发明要实现的主要目标。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种井中微地震纵横波事件初至识别方法，利用矢量分解、结合最大特征值与改进能量比法，尽可能地突出井中微地震P波、S波事件初至起跳位置。

本发明是通过以下技术方案实现的

一种井中微地震纵横波事件初至识别方法，利用两次矢段曲线偏振分析与旋转处理，完成井中微地震三分量矢量波场分离，然后构造协方差矩阵，再根据长短时窗能量比以及引入权系数，获得事件特征曲线，最后根据门槛值搜索识别出P波、S波事件初至到达时间位置。

所述方法包括：

第一步，通过矢端曲线分析方法，对井中微地震原始三分量中的X、Y分量数据进行偏振分析与旋转处理，分别获得波传播水平方位角α、径向分量R与切向分量T；

第二步，通过矢端曲线分析方法，对径向分量R与原始Z分量数据进行偏振分析与旋转处理，分别获得波传播垂向方位角β、沿波传播方向P分量、垂直波传播方向S分量；

第三步，利用P分量、T分量构建包含P波信息的瞬时协方差矩阵CC_Pi，利用S分量、T分量构建包含S波信息的瞬时协方差矩阵CC_Si，最后识别出不同震相P波、S波达到检波点接收器初至时间。

所述第一步是这样实现的：

建立坐标系，横坐标为X分量振幅值，垂直坐标为Y分量振幅值，将X、Y分量所有样点值(x_i,y_i)投影到坐标系上，并按时间顺序将所有样点坐标连线，形成一个闭合曲线，该曲线称为矢端曲线；线性拟合矢端曲线得到近似线性关系式y_i＝tan(α)·x_i，通过如下公式估算出与横坐标斜率tan(α)

$\tan \left(\mathrm{\α}\right)=\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_{M-1}\\ y_M\end{array}\right][x_1,x_2,...,x_{M-1},x_M]\·{\left[\left[\begin{array}{c}{x}_1\\ x_2\\ .\\ .\\ .\\ x_{M-1}\\ x_M\end{array}\right]\right[x_1,x_2,...,x_{M-1},x_M\left]\right]}^{-1},$

其中M为所有样点个数，最后反正切即可求出X、Y分量水平方位角α；

根据旋转公式，对原始X、Y分量数据进行旋转处理：

R_i＝x_i cos(α)+y_i sin(α)

(1)

T_i＝-x_i sin(α)+y_i cos(α)

获得径向分量R与切向分量T。

所述第二步是这样实现的：

建立坐标系，横坐标为原始Z分量振幅，垂直坐标为R振幅，将Z、R分量所有样点值(z_i,R_i)投影到坐标系上，并按时间顺序将所有样点坐标连线，形成一个闭合曲线，即矢端曲线；线性拟合矢端曲线R_i＝tan(β)·z_i，估算出波传播垂直偏振角β。计算过程类似于上面水平方向线性拟合y_i＝tan(α)·x_i，计算线性斜率公式为：

$\tan \left(\mathrm{\β}\right)=\left[\begin{array}{c}{R}_1\\ R_2\\ .\\ .\\ .\\ R_{M-1}\\ R_M\end{array}\right][z_1,z_2,...,z_{M-1},z_M]\·{\left[\left[\begin{array}{c}{z}_1\\ z_2\\ .\\ .\\ .\\ z_{M-1}\\ z_M\end{array}\right]\right[z_1,z_2,...,z_{M-1},z_M\left]\right]}^{-1},$

其中M为所有样点个数；

根据旋转公式，对原始Z、R分量数据进行旋转处理：

P_i＝z_i cos(β)+R_i sin(β)

S_i＝-z_i sin(β)+R_i cos(β) (2)

获得沿波传播方向P分量、垂直波传播方向S分量。

所述第三步包括：

(1)，定义样点时窗大小，分别构建P波瞬时协方差矩阵CC_Pi、S波瞬时协方差矩阵CC_Si：

${\mathrm{CC}}_{\mathrm{Pi}}=\left[\begin{array}{c}\mathrm{\Σ}{P}_i\·P_i,\mathrm{\Σ}{P}_i\·T_i\\ \mathrm{\Σ}{T}_i\·P_i,\mathrm{\Σ}{T}_i\·T_i\end{array}\right]---\left(3\right)$

${\mathrm{CC}}_{\mathrm{Si}}=\left[\begin{array}{c}\mathrm{\Σ}{S}_i\·S_i,\mathrm{\Σ}{S}_i\·T_i\\ \mathrm{\Σ}{T}_i\·S_i,\mathrm{\Σ}{T}_i\·T_i\end{array}\right]---\left(4\right)$

其中，∑表示对样点时窗内数值求和；

(2)分别对矩阵CC_Pi、CC_Pi进行奇异值SVD分解，计算出相应最大特征值λ_Pi、λ_Si：

CC_Pi＝{λ_P1,λ_P2}U_PU_P ^-1

(5)

CC_Si＝{λ_S1,λ_S2}U_SU_S ^-1

λ_Pi＝max{λ_P1,λ_P2} (6)

λ_Si＝max{λ_S1,λ_S2}

其中，{λ_P1,λ_P2}、{λ_S1,λ_S2}分别表示P波、S波特征值矩阵，U_P、U_S则表示对应的特征向量矩阵；

(3)，再定义长时窗L₁、短时窗L₂，分别对所有λ_Pi、λ_Si数据进行能量比计算：

${\mathrm{ER}}_{\mathrm{Pi}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Pi}}}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Pi}}}^2/L_2}---\left(7\right)$

${\mathrm{ER}}_{\mathrm{Si}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Si}}}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Si}}}^2/L_2}---\left(8\right)$

其中，ER_Pi、ER_Si为样点i最大特征值曲线P波、S波长短时窗能量比。

(4)，为了获得更加稳定的识别函数，将最大特征值λ_Pi、λ_Si作为权系数，对ER_Pi、ER_Si进行修正，获得改进的能量比公式：

MER_Pi＝ER_Pi ³·λ_Pi ³

(11)

MER_Si＝ER_Si ³·λ_Si ³

(5)，识别出不同震相P波、S波达到检波点接收器初至时间。

所述步骤(5)是通过设定门槛值，搜索大于门槛值区域内最大峰值对应的样点时间来实现的，具体包括：

(A)，定义设定门槛值K_P、K_S，搜索P波、S波改进的能量比MER_Pi、MER_Si曲线第一次大于门槛值时对应的时间样点位置T_KP、T_KS，即出现事件判断标准：

(B)，分别定义初至搜索时窗大小W_P、W_S，找出时窗内最大峰值F_P、F_S：

F_P＝max{MER_Pi}且Pi∈[T_KP,T_KP+W_P] (14)

F_S＝max{MER_Si}且Pi∈[T_KS,T_KS+W_S] (15)

(C)，将搜索出的最大峰值F_P、F_S对应的时间样点t_P、t_S输出，即为最终所求的井中微地震P波、S波事件初至时间。

与现有技术相比，本发明的有益效果是：在井中微地震三分量处理中，传统长短时窗能量比法与本发明在P波初至拾取上，都能够获得较好的识别曲线，但是，本发明识别曲线波动性要小得多，设置简单门槛值就能实现初至拾取；在S波初至拾取上，传统方法识别曲线S波初至时间峰值不突出、较难拾取，而本发明却仍然能够获得较好的“尖脉冲”识别曲线，优势明显，这一点体现了本发明有较强的实用性。

附图说明

图1是本发明实现井中微地震纵横波初至识别操作流程图；

图2是建立井中微地震模型示意图；

图3a是模拟井中微地震三分量数据Z分量；

图3b是模拟井中微地震三分量数据X分量；

图3c是模拟井中微地震三分量数据Y分量；

图4a是模型数据偏振分析与旋转处理后获得新的三分量数据P分量；

图4b是模型数据偏振分析与旋转处理后获得新的三分量数据S分量；

图4c是模型数据偏振分析与旋转处理后获得新的三分量数据T分量；

图5a是初至非常明显的P波数据；

图5b是初至比较明显的P波数据；

图5c是初至不明显的P波数据；

图6a是初至非常明显的S波数据；

图6b是初至比较明显的S波数据；

图6c是初至不明显的S波数据；

图7a是P波图5a长短时窗能量比识别曲线；

图7b是P波图5b长短时窗能量比识别曲线；

图7c是P波图5c长短时窗能量比识别曲线；

图8a是S波图6a长短时窗能量比识别曲线；

图8b是S波图6b长短时窗能量比识别曲线；

图8c是S波图6c长短时窗能量比识别曲线；

图9a是P波图5a最大特征值曲线；

图9b是P波图5b最大特征值曲线；

图9c是P波图5c最大特征值曲线；

图10a是P波图9a改进能量比法识别曲线；

图10b是P波图9b改进能量比法识别曲线；

图10c是P波图9c改进能量比法识别曲线；

图11a是S波图6a最大特征值曲线；

图11b是S波图6b最大特征值曲线；

图11c是S波图6c最大特征值曲线；

图12a是S波图11a改进能量比法识别曲线；

图12b是S波图11b改进能量比法识别曲线；

图12c是S波图11c改进能量比法识别曲线；

图13a是实际资料井中微地震三分量数据：Z分量；

图13b是实际资料井中微地震三分量数据：X分量；

图13c是实际资料井中微地震三分量数据：Y分量；

图14a是实际资料偏振分析与旋转处理后获得新的三分量数据：P分量；

图14b是实际资料偏振分析与旋转处理后获得新的三分量数据：S分量；

图14c是实际资料偏振分析与旋转处理后获得新的三分量数据：T分量；

图15a是P波实际资料初至较明显拾取方法对比中P波单道微地震波形显示；

图15b是P波实际资料初至较明显拾取方法对比中P波传统长短时窗能量比识别曲线；

图15c是P波实际资料初至较明显拾取方法对比中P波最大特征值曲线；

图15d是P波实际资料初至较明显拾取方法对比中P波本发明改进能量比识别曲线；

图16a是S波实际资料初至较明显拾取方法对比中S波单道微地震波形显示；

图16b是S波实际资料初至较明显拾取方法对比中S波传统长短时窗能量比识别曲线；

图16c是S波实际资料初至较明显拾取方法对比中S波最大特征值曲线；

图16d是S波实际资料初至较明显拾取方法对比中S波本发明改进能量比识别曲线；

图17a是P波实际资料初至不明显拾取方法对比中P波单道微地震波形显示；

图17b是P波实际资料初至不明显拾取方法对比中P波传统长短时窗能量比识别曲线；

图17c是P波实际资料初至不明显拾取方法对比中P波最大特征值曲线；

图17d是P波实际资料初至不明显拾取方法对比中P波本发明改进能量比识别曲线；

图18a是S波实际资料初至不明显拾取方法对比中S波单道微地震波形显示；

图18b是S波实际资料初至不明显拾取方法对比中S波传统长短时窗能量比识别曲线；

图18c是S波实际资料初至不明显拾取方法对比中S波最大特征值曲线；

图18d是S波实际资料初至不明显拾取方法对比中S波本发明改进能量比识别曲线。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明方法包括：

第一步，通过矢端曲线分析方法，对井中微地震原始三分量中的X、Y分量数据进行偏振分析与旋转处理，分别获得波传播水平方位角α、径向分量R与切向分量T；

第二步，通过矢端曲线分析方法，对径向分量R与原始Z分量数据进行偏振分析与旋转处理，分别获得波传播垂向方位角β、沿波传播方向P分量、垂直波传播方向S分量；根据P波、S波沿传播方向与质点振动方向物理特征，可以判断，P分量主要包含了微地震P波信息，S分量主要包含了微地震S波信息，而切向分量T理主要是噪音。

第三步，构建P分量、切向分量T二维瞬时协方差矩阵CC_Pi，计算其瞬时最大特征值λ_Pi；对所有最大特征值样点值计算长短时窗能量比曲线(又称为瞬时能量比数据)ER_Pi，并将最大特征值λ_Pi作为权系数，获得P分量改进的能量比MER_Pi曲线，最后根据门槛值在特征曲线MER_Pi数据上搜索极大值出现的时间位置，即为识别出的P波事件初至时间t_P：

第四步，操作同第三步，构建S分量、T分量瞬时协方差矩阵CC_Si，相应计算出包含S波信息的瞬时最大特征值λ_Si、长短时窗能量比曲线(又称能量比数据)ER_Si，同样，将λ_Si作为权系数，获得S分量改进的能量比MER_Si曲线，即特征曲线MER_Si，最后根据门槛值搜索识别出S波事件初至时间t_S：设定P波、S波初至识别门槛值，分别在特征曲线MER_Pi、MER_Si上搜索大于门槛值的区域最大峰值，而其峰值对应的时间样点，即为所求的P波、S波事件初至时间。

其中，对井中微地震三分量数据进行偏振分析与旋转处理，目的是实现矢量分解，从原始三分量X、Y、Z中获得新的三分量P、S、T，本发明采用已有的算法--矢端曲线分析法：

首先对原始三分量X、Y、Z数据进行两次偏振分析与旋转处理，获得新的三分量P、S、T，本发明采用已有算法--矢端曲线分析法：

建立坐标系，横坐标为X分量振幅值，垂直坐标为Y分量振幅值，将X、Y分量所有样点值(x_i,y_i)投影到坐标系上，并按时间顺序将所有样点坐标连线，形成一个闭合曲线，称为矢端曲线。线性拟合矢端曲线(y_i＝α·x_i)，估算与横坐标交角，即为所求X、Y分量水平方位角α。

根据旋转公式，对原始X、Y分量数据进行旋转处理：

R_i＝x_i cos(α)+y_i sin(α)

(1)

T_i＝-x_i sin(α)+y_i cos(α)

获得新的两分量R、T。

以同样操作，建立坐标系，横坐标为原始Z分量振幅，垂直坐标为R振幅，画出(z_i,R_i)矢端曲线，线性拟合(R_i＝β·z_i)，估算出波传播垂直偏振角β，并进一步作旋转处理：

P_i＝z_i cos(α)+R_i sin(α)

(2)

S_i＝-z_i sin(α)+R_i cos(α)

再一次获得新的两分量组合P、S。

从物理意义上，经过两次偏振分析与旋转后，沿波传播方向的P分量主要包含了震相P波微地震数据，垂直波传播方向的S分量则主要包含了震相S波类型微地震数据，而第一次获得的水平面切向分量T则主要包含随机噪音。因此，按照本次操作，可以实现将原始三分量X、Y、Z矢量分解后获得新的三分量P、S、T。

然后，构建不同震相P波、S波初至识别的特征曲线，本发明通过建立协方差、最大特征值、长短时窗能量比、权系数能量比等计算过程来实现，具体如下：

首先，定义样点时窗大小，分别构建P波瞬时协方差矩阵、S波瞬时协方差矩阵：

${\mathrm{CC}}_{\mathrm{Pi}}=\left[\begin{array}{c}\mathrm{\Σ}{P}_i\·P_i,\mathrm{\Σ}{P}_i\·T_i\\ \mathrm{\Σ}{T}_i\·P_i,\mathrm{\Σ}{T}_i\·T_i\end{array}\right]---\left(3\right)$

${\mathrm{CC}}_{\mathrm{Si}}=\left[\begin{array}{c}\mathrm{\Σ}{S}_i\·S_i,\mathrm{\Σ}{S}_i\·T_i\\ \mathrm{\Σ}{T}_i\·S_i,\mathrm{\Σ}{T}_i\·T_i\end{array}\right]---\left(4\right)$

其中，∑表示对样点时窗内数值求和。

然后分别对矩阵CC_Pi、CC_Pi奇异值SVD分解，计算出相应最大特征值λ_Pi、λ_Si：

CC_Pi＝{λ_P1,λ_P2}U_PU_P ^-1

(5)

CC_Si＝{λ_S1,λ_S2}U_SU_S ^-1

λ_Pi＝max{λ_P1,λ_P2} (6)

λ_Si＝max{λ_S1,λ_S2}

其中，{λ_P1,λ_P2}、{λ_S1,λ_S2}分别表示P波、S波特征值矩阵，U_P、U_S则表示对应的特征向量矩阵。

其次，再定义长时窗L₁、短时窗L₂，分别对所有λ_Pi、λ_Si数据进行能量比计算：

${\mathrm{ER}}_{\mathrm{Pi}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Pi}}}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Pi}}}^2/L_2}---\left(7\right)$

${\mathrm{ER}}_{\mathrm{Si}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Si}}}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Si}}}^2/L_2}---\left(8\right)$

其中，ER_Pi、ER_Si为样点i最大特征值曲线P波、S波长短时窗能量比。

这一点，不同于传统长短时窗能量比计算方法：

${\mathrm{ER}}_{\mathrm{PPi}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{P_i}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{P_i}^2/L_2}---\left(9\right)$

${\mathrm{ER}}_{\mathrm{SSi}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{S_i}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{S_i}^2/L_2}---\left(10\right)$

从公式(9)、(10)中可以看出，传统方法直接利用P波、S波振幅信息计算识别曲线ER。

最后，为了获得更加稳定的识别函数，将最大特征值λ_Pi、λ_Si作为权系数，对ER_Pi、ER_Si进行修正，获得改进的能量比公式：

MER_Pi＝ER_Pi ³·λ_Pi ³

(11)

MER_Si＝ER_Si ³·λ_Si ³

最后，识别出不同震相P波、S波达到检波点接收器初至时间，本发明通过设定门槛值，搜索大于门槛值区域内最大峰值对应的样点时间来实现：

首先，定义设定门槛值K_P、K_S，搜索P波、S波改进的能量比MER_Pi、MER_Si曲线第一次大于门槛值时对应的时间样点位置T_KP、T_KS，即出现事件判断标准：

然后，再分别定义初至搜索时窗大小W_P、W_S，找出时窗内最大峰值F_P、F_S：

F_P＝max{MER_Pi}且Pi∈[T_KP,T_KP+W_P] (14)

F_S＝max{MER_Si}且Pi∈[T_KS,T_KS+W_S] (15)

最后，将索搜出的最大峰值F_P、F_S对应的时间样点t_P、t_S输出，即为最终所求的井中微地震P波、S波事件初至时间。

根据上述步骤可知，与传统长短时窗不同在于，本发明不是直接利用长短时窗，而是在偏振分析、旋转处理之后，获得P波、S波信号占优的基础上，利用协方差，奇异值分解计算最大特征值，也就是说，在此基础上再利用长短时窗能量比，同时引入权系数，为了获得更加稳定的抗噪能力较强的初至识别特征曲线MER。简言之，本发明有两个创新点：一是在偏振分析矢量分解后纵横波最大特征曲线上利用长短时窗；而是引入权系数，提出改进能量比思想。

本发明提供的井中微地震震相P波、S波事件初至拾取方法，具体实施需要分四个步骤，如图1所示操作流程：从原始三分量数据，利用偏振分析与旋转处理，获得新的三分量P、S、T；构建P与T、S与T协方差，通过奇异值分解，获得瞬时最大特征值；利用长短时窗，同时引入权系数，获得改进的能量比特征曲线；给定门槛值，搜索纵横波改进的能量比特征曲线上，符合“大于门槛值且区域内最大峰值”，其对应的样点时间即为所求的初至时间。

下面分别从井中微地震理论模型数据和实际资料应用中说明本发明事件初至拾取效果。

开展理论模型初至拾取测试，为了验证本发明的可行性。设计了一个在一级3C检波器上接收到的井中微地震有效事件记录，如图2所示。在这里，模拟压裂微地震为双震源，即同时激发P波与S波震相信号。如图3a-图3c所示，模拟出井中单道微地震三分量Z、X、Y数据，且P波初至时间为100ms、S波初至时间为250ms。根据图1操作流程，利用已知矢端曲线算法，偏振分析出水平方位角α、垂直方位角β，对原始三分量进行旋转处理，获得沿波传播方向P波、垂直波传播方向S波以及切向T分量，如图4a-图4c所示。为了对比与传统方法，突出本发明抗噪能力，在P波、S波数据上增加三种不同程度随机噪音，使得波的初至相位非常明显、较明显、不明显等三种不同信噪比微震数据(如图5a至图6c所示)。

传统初至拾取方法是直接在微地震信号上计算长短时窗能量比，如图7a-图7c为P波三种不同信噪比微震数据长短时窗能量比曲线，图8a-图8c为S波三种不同信噪比微震数据长短时窗能量比曲线。从图7a-图7c中可以看出，设置一个简单的门槛值，就能够准确拾取三种不同信噪比数据P波初至，而相对于S波，图8a、图8b也能够容易地拾取S波初至，但是对于图8c，识别曲线波动较大，必须设置较大门槛值，才能准确拾取S波初至。传统方法能够较好地拾取P波初至，相比较S波，初至相位受噪音影响较大，识别曲线波动大，特别是低信噪比数据，可能因为门槛值设置不准确，导致错误的结果。

按照流程图1，用本发明来验证模型数据初至拾取情况，如图9a至图12c所示。其中，图9a-图9c是P波三种不同信噪比最大特征值曲线，图10a-图10c是图9a-图9c对应的三种改进能量比法识别曲线，图11a-图11c是S波三种不同信噪比最大特征值曲线，图12a-图12c是图11a-图11c对应的三种改进能量比法识别曲线。从图10a-图10c、图12a-图12c中可以看出，无论是P波还是S波，其识别初至特征曲线在“初至”时间表现出较好的“尖脉冲”且曲线波动非常小，也就是说，设置简单门槛值，就能够准确拾取各种信噪比数据P波、S波微地震初至时间。相比较传统方法，验证了本发明在低信噪比、特别是S波初至拾取方面稳定性较好，保证了初至时间拾取的准确性较可靠。

而井中微地震实际资料由于受到周边环境影响，使得微地震数据较为复杂，可能会出现大部分微地震事件处于中、低信噪比数据中，如图13a-图13c为某煤层气区井中微地震实际资料三分量单道数据波，存在4组振幅不同的P波、S波微地震对，即代表了4种不同信噪比微震数据。在水平方向、垂直方向上分别进行矢端曲线偏振分析法，然后对原始井中微地震三分量进行旋转处理，获得沿波传播方向的P波分量、垂直波传播方向S波分量以及切向T分量，如图14a-图14c所示，基本实现了原始三分量矢量分解处理。

首先，当实际资料P波、S波初至比较明显时，比较本发明与传统方法识别曲线，如图15a-图16d所示。图15b为传统长短时窗能量比法P波初至识别曲线，图15d为本发明改进能量比法P波初至识别曲线，两者相比较，虽然传统方法在初至最大峰值较为明显，但是本发明识别曲线初至“尖脉冲”却更加突出且曲线波动非常少。对于S波，传统方法识别曲线(如图16b)波动较大，已经不能简单地设置门槛值来识别S波初至了，而本发明识别曲线(如图16b)S波初至“尖脉冲”非常突出，很容易搜索到S波初至时间。

如果实际资料P波、S波初至不明显时，本发明是否还能够获得较好识别曲线呢？如图17a-图18d所示，比较图17b、图17d，虽然传统方法仍然能够识别出P波初至时间位置，但是本发明识别曲线波动非常小，更容易识别初至；同样，比较图18b、图18d，传统方法已经无法识别出S波初至峰值了，而本发明识别曲线波动性较小，比较容易通过设置简单门槛值就能找到S波初至时间“最大峰值”。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (6)

1.一种井中微地震纵横波事件初至识别方法，其特征在于：所述方法利用两次矢段曲线偏振分析与旋转处理，完成井中微地震三分量矢量波场分离，然后构造协方差矩阵，再根据长短时窗能量比以及引入权系数，获得事件特征曲线，最后根据门槛值搜索识别出P波、S波事件初至到达时间位置。

2.根据权利要求1所述的井中微地震纵横波事件初至识别方法，其特征在于：所述方法包括：

第一步，通过矢端曲线分析方法，对井中微地震原始三分量中的X、Y分量数据进行偏振分析与旋转处理，分别获得波传播水平方位角α、径向分量R与切向分量T；

第二步，通过矢端曲线分析方法，对径向分量R与原始Z分量数据进行偏振分析与旋转处理，分别获得波传播垂向方位角β、沿波传播方向P分量、垂直波传播方向S分量；

第三步，利用P分量、T分量构建包含P波信息的瞬时协方差矩阵CC_Pi，利用S分量、T分量构建包含S波信息的瞬时协方差矩阵CC_Si，最后识别出不同震相P波、S波达到检波点接收器初至时间。

3.根据权利要求2所述的井中微地震纵横波事件初至识别方法，其特征在于：所述第一步是这样实现的：

建立坐标系，横坐标为X分量振幅值，垂直坐标为Y分量振幅值，将X、Y分量所有样点值(x_i,y_i)投影到坐标系上，并按时间顺序将所有样点坐标连线，形成一个闭合曲线，即矢端曲线；线性拟合矢端曲线y_i＝α·x_i，估算出与横坐标交角，即为所求X、Y分量水平方位角α；

根据旋转公式，对原始X、Y分量数据进行旋转处理：

R_i＝x_i cos(α)+y_i sin(α)

(1)

T_i＝-x_i sin(α)+y_i cos(α)

获得径向分量R与切向分量T。

4.根据权利要求2所述的井中微地震纵横波事件初至识别方法，其特征在于：所述第二步是这样实现的：

建立坐标系，横坐标为原始Z分量振幅，垂直坐标为R振幅，将Z、R分量所有样点值(z_i,R_i)投影到坐标系上，并按时间顺序将所有样点坐标连线，形成一个闭合曲线，即矢端曲线；线性拟合矢端曲线R_i＝β·z_i，估算出波传播垂直偏振角β；

根据旋转公式，对原始Z、R分量数据进行旋转处理：

P_i＝z_i cos(α)+R_i sin(α)

(2)

S_i＝-z_i sin(α)+R_i cos(α)

获得沿波传播方向P分量、垂直波传播方向S分量。

5.根据权利要求4所述的井中微地震纵横波事件初至识别方法，其特征在于：所述第三步包括：

(1)，定义样点时窗大小，分别构建P波瞬时协方差矩阵CC_Pi、S波瞬时协方差矩阵CC_Si：

${\mathrm{CC}}_{\mathrm{Pi}}=\left[\begin{array}{c}\mathrm{\Σ}{P}_i\·P_i,\mathrm{\Σ}{P}_i\·T_i\\ \mathrm{\Σ}{T}_i\·P_i,\mathrm{\Σ}{T}_i\·T_i\end{array}\right]---\left(3\right)$

${\mathrm{CC}}_{\mathrm{Si}}=\left[\begin{array}{c}\mathrm{\Σ}{S}_i\·S_i,\mathrm{\Σ}{S}_i\·T_i\\ \mathrm{\Σ}{T}_i\·S_i,\mathrm{\Σ}{T}_i\·T_i\end{array}\right]---\left(4\right)$

其中，∑表示对样点时窗内数值求和；

(2)分别对矩阵CC_Pi、CC_Pi进行奇异值SVD分解，计算出相应最大特征值λ_Pi、λ_Si：

CC_Pi＝{λ_P1,λ_P2}U_PU_P ^-1

(5)

CC_Si＝{λ_S1,λ_S2}U_SU_S ^-1

λ_Pi＝max{λ_P1,λ_P2} (6)

λ_Si＝max{λ_S1,λ_S2}

其中，{λ_P1,λ_P2}、{λ_S1,λ_S2}分别表示P波、S波特征值矩阵，U_P、U_S则表示对应的特征向量矩阵；

(3)，再定义长时窗L₁、短时窗L₂，分别对所有λ_Pi、λ_Si数据进行能量比计算：

${\mathrm{ER}}_{\mathrm{Pi}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Pi}}}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Pi}}}^2/L_2}---\left(7\right)$

${\mathrm{ER}}_{\mathrm{Si}}=\frac{\underset{i-L_1}{\overset{i}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Si}}}^2/L_1}{\underset{i}{\overset{i+L_2}{\mathrm{\Σ}}}{{\mathrm{\λ}}_{\mathrm{Si}}}^2/L_2}---\left(8\right)$

其中，ER_Pi、ER_Si为样点i最大特征值曲线P波、S波长短时窗能量比；

(4)，为了获得更加稳定的识别函数，将最大特征值λ_Pi、λ_Si作为权系数，对ER_Pi、ER_Si进行修正，获得改进的能量比公式：

MER_Pi＝ER_Pi ³·λ_Pi ³

(11)

MER_Si＝ER_Si ³·λ_Si ³

(5)，识别出不同震相P波、S波达到检波点接收器初至时间。

6.根据权利要求5所述的井中微地震纵横波事件初至识别方法，其特征在于：所述步骤(5)是通过设定门槛值，搜索大于门槛值区域内最大峰值对应的样点时间来实现的，具体包括：

(A)，定义设定门槛值K_P、K_S，搜索P波、S波改进的能量比MER_Pi、MER_Si曲线第一次大于门槛值时对应的时间样点位置T_KP、T_KS，即出现事件判断标准：

(B)，分别定义初至搜索时窗大小W_P、W_S，找出时窗内最大峰值F_P、F_S：

F_P＝max{MER_Pi}且Pi∈[T_KP,T_KP+W_P] (14)

F_S＝max{MER_Si}且Pi∈[T_KS,T_KS+W_S] (15)

(C)，将搜索出的最大峰值F_P、F_S对应的时间样点t_P、t_S输出，即为最终所求的井中微地震P波、S波事件初至时间。