微震偏移叠加定位方法、装置、电子设备及可读存储介质

技术领域

本发明属于微震监测资料处理技术领域，具体涉及一种微震偏移叠加定位方法、装置、电子设备及可读存储介质。

背景技术

近年来，基于偏移成像的定位方法因无需拾取初至到时而得到广泛应用。该基于偏移成像的定位方法的思路是：在整个时间、空间范围内寻找可能发生的微震震源。这是一种新的可对震源分布进行成像的方法，利用监测数据的振幅和到时等波形数据，通过系统的扫描某个假设的震源位置和起始时刻，找出整个地震序列分布。在不用精确拾取到时和计算理论地震图的情况下，达到比较理想的定位结果。随着微震监测技术需求的发展，这类定位方法发展迅速，已经成为了微震监测中主要的定位方法之一。

天然地震及诱发(微)地震事件的震源机制中往往包含双力偶成分，使得布设在地表的检波器所接收到的微地震事件波形出现极性反转现象，即初至极性方向具有正负之分。这导致了常规波形叠加成像函数的失效，使得本该聚集在真实震源处的成像能量发生衰减和分散，造成微地震事件的探测失败和定位误差增大。目前提出的三种解决策略：一是通过反演震源机制来校正极性，但是计算成本高，且震源机制与震源位置耦合性强，不确定性大；二是利用聚焦算子将分散的成像能量归位到真实震源处，但要求规则且密集的台网布设；三是通过特征函数将波形数据转化成正值序列的波形数据，该策略计算成本低且不受台网布设的限制，但却失去了噪音的不相干性特征，导致成像域噪音增多，影响定位精度。

因此，有必要研究一种合适的微震偏移叠加定位方法，使其极性反转的微地震信号波形都一致，且同时保存噪音的不相干性特征。以此达到既增强微震信号叠加的成像能量，将分散的成像能量重新归位到真实震源位置；又能借助噪音的不相干性压制成像域噪音，提高成像能量分辨率，达到提升定位精度的目的。

发明内容

为了解决上述技术问题，本发明提供了一种微震偏移叠加定位方法、装置、电子设备及可读存储介质。

本发明公开一种波形一致性转换的微震偏移叠加定位方法，包括：

步骤1、对实际微震波形数据进行常规带通滤波预处理；

步骤2、利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；

步骤3、将滤波后波形的归一化解析形式对所述滤波后波形进行相乘，得到转换波形的解析形式；

步骤4、将所述转换波形的解析形式的实部减去所述滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；

步骤5、按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算所述三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；

步骤6、设定发震时刻间隔，利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，其中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三维空间的能量成像体；

步骤7、获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；

步骤8、设定成像能量阈值，所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。

作为本发明的进一步改进，步骤2中利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理，包括：

希尔伯特变换公式为以下公式1：

公式1：

其中，H[s(t)]为滤波后波形s(t)的希尔伯特变换，将其作为滤波后波形解析形式的虚部，s(t)作为解析形式的实部，得到各地震道滤波后波形的解析形式的第一种表达方式为以下公式2：

公式2：s

其中，s

各地震道滤波后波形的解析形式的第二种表达方式为以下公式3：

公式3：s

其中，

为波形的瞬时相位；

对各地震道滤波后波形的解析形式s

公式4：

其中，Norm(s

作为本发明的进一步改进，步骤3中将各地震道滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式，包括：

根据以下公式5得到转换波形的解析形式：

公式5：x

＝2A(t)cos(φ(t))*e

＝A(t)[cos(2φ(t))+1]+iA(t)sin(2φ(t))

其中，x

作为本发明的进一步改进，步骤4中将转换波形解析形式的实部减去滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据，包括：

根据以下公式6得到波形一致性转换后的波形数据：

公式6：T(t)＝real(x

＝A(t)[cos(2φ(t))+1]-A(t)

＝A(t)cos(2φ(t))

其中，T(t)为滤波后波形s(t)经过波形一致性转换后的波形数据，即作为后续偏移成像定位的输入波形数据。

作为本发明的进一步改进，步骤6中利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，包括：

根据以下公式7得到所述四维成像能量体；

公式7：

其中，

作为本发明的进一步改进，步骤7中获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线，包括：

根据以下公式78计算微震事件检测曲线；

公式8：

其中，

作为本发明的进一步改进，步骤8中所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，包括：

若设定成像能量阈值THRE，则满足

本发明还公开了一种波形一致性转换的微震偏移叠加定位的装置，包括：

预处理模块，用于对实际微震波形数据进行常规带通滤波预处理；

变换模块，用于利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；

相乘模块，用于将滤波后波形的归一化解析形式对所述滤波后波形进行相乘，得到转换波形的解析形式；

相减模块，用于将所述转换波形的解析形式的实部减去所述滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；

计算模块，用于按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算所述三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；

叠加模块，用于设定发震时刻间隔，利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，其中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三维空间的能量成像体；

获取模块，用于获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；

定位模块，用于设定成像能量阈值，所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。

本发明还公开了一种电子设备，包括：处理器、通信接口、存储器和通信总线，其中，处理器、通信接口、存储器通过通信总线完成相互间的通信；

所述存储器中存储有计算机程序，当所述程序被所述处理器执行时，使得所述处理器执行上述波形一致性转换的微震偏移叠加定位方法的步骤。

本发明还公开了一种存储介质，其存储有可由电子设备执行的计算机程序，当所述程序在所述电子设备上运行时，使得所述电子设备执行上述波形一致性转换的微震偏移叠加定位方法的步骤。

与现有技术相比，本发明的有益效果为：

利用各地震道上微震事件信号的相位差为零或者±π的特性，将原波形的归一化解析形式对原波形进行相乘加权，可以得到波形一致性的转换波形，可以将极性反转引起的分散成像能量归位到震源位置处。此外，还可以使得转换波形中的噪音成分仍然保存着不相干性特征，从而压制成像域中的噪音能量，提高微震事件的探测率和成像能量分辨率，提升定位精度。

附图说明

为了更清楚地说明本申请的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对本申请保护范围的限定。在各个附图中，类似的构成部分采用类似的编号。

图1示出了本申请实施例提供的波形一致性转换的微震偏移成像定位方法的流程示意图；

图2示出了本申请实施例提供的极性相反的微震信号的相位差原理图；

图3示出了本申请实施例提供的原始波形及波形一致性转换的波形数据的对比示意图；

图4示出了本申请实施例提供的实际煤层气水力压裂的地面监测系统及P波速度模型；

图5示出了本申请实施例提供的实际原始微震波形记录及波形一致性转换后的波形记录对比图；

图6示出了本申请实施例提供的微震波形记录进行微地震事件检测结果对比图；

图7示出了本申请实施例提供的利用所探测到的微震事件的发震时刻及发震位置所预测的初至到时曲线在原始波形记录上的投影图；

图8示出了本申请实施例提供的E4微震事件进行定位成像的对比图；

图9示出了本申请实施例提供的90秒长的实际微震记录进行微地震事件检测及定位结果对比图。

具体实施方式

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

本发明提供一种波形一致性转换的微震偏移叠加定位方法、装置、电子设备及可读存储介质，首先对所有台站的微震波形数据进行常规带通滤波预处理；利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；然后，将所述滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式；将所述转换波形解析形式的实部减去滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；将定位成像区域进行三维网格划分，计算每个网格到所有检波器的初至旅行时间，得到走时表数据；设定合理的发震时刻间隔，将输入波形数据根据走时表偏移叠加得到四维成像能量体；取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；最后，设定成像能量阈值，微震事件探测曲线中超过阈值的发震时刻信息，即为成功检测出来的微震事件的发震时刻和发震位置。

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

如图1所示，本发明提供一种波形一致性转换的微震偏移叠加定位方法，包括：

步骤1、对实际微震波形数据进行常规带通滤波预处理；

步骤2、利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；其中，

希尔伯特变换公式为：

其中，H[s(t)]为滤波后波形s(t)的希尔伯特变换，将其作为滤波后波形解析形式的虚部，s(t)作为解析形式的实部，即得到各地震道滤波后波形的解析形式：

s

其中，s

上式也可以表示为：

s

其中，

为波形的瞬时相位。

对各地震道滤波后波形的解析形式s

/>

其中Norm(s

步骤3、将滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式。计算公式为：

x

＝2A(t)cos(φ(t))*e

＝A(t)[cos(2φ(t))+1]+iA(t)sin(2φ(t))

其中，x

步骤4、将转换波形解析形式的实部减去滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；

计算公式如下：

T(t)＝real(x

＝A(t)[cos(2φ(t))+1]-A(t)

＝A(t)cos(2φ(t))

其中，T(t)为滤波后波形s(t)经过波形一致性转换后的波形数据，即作为后续偏移成像定位的输入波形数据。

步骤5、按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；

步骤6、设定合理的发震时刻间隔，利用成像函数将输入波形数据沿着各个网格点到所有检波器的初至到时(发震时刻+初至旅行时间)曲线进行偏移叠加，得到四维(发震时刻+三维空间)成像能量体。成像函数的计算公式为：

其中，

在本实施例中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三位空间的成像能量体。

步骤7、获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线。微震事件检测曲线的计算公式为：

其中，

步骤8、设定成像能量阈值，所述微震事件探测曲线中超过阈值的发震时刻点所包含的信息，即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。

具体的，设定合理的成像能量阈值，若设定成像能量阈值THRE，则满足

本发明提供一种用于实现上述微震偏移叠加定位方法的装置，包括：

预处理模块，用于对实际微震波形数据进行常规带通滤波预处理；

变换模块，用于利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；

相乘模块，用于将滤波后波形的归一化解析形式对所述滤波后波形进行相乘，得到转换波形的解析形式；

相减模块，用于将所述转换波形的解析形式的实部减去所述滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；

计算模块，用于按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算所述三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；

叠加模块，用于设定发震时刻间隔，利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，其中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三维空间的能量成像体；

获取模块，用于获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；

定位模块，用于设定成像能量阈值，所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。

本发明提供一种电子设备，包括：处理器、通信接口、存储器和通信总线，其中，处理器、通信接口、存储器通过通信总线完成相互间的通信；存储器中存储有计算机程序，当程序被处理器执行时，使得处理器执行上述波形一致性转换的微震偏移成像定位方法的步骤。

该电子设备的技术方案与上述波形一致性转换的微震偏移成像定位方法的技术方案属于同一构思，电子设备的技术方案未详细描述的细节内容，均可以参见上述波形一致性转换的微震偏移叠加定位方法的技术方案的描述。

电子设备可以是任何类型的静止或移动电子设备，包括移动计算机或移动电子设备(例如，平板计算机、个人数字助理、膝上型计算机、笔记本计算机、上网本等)、移动电话(例如，智能手机)、可佩戴的电子设备(例如，智能手表、智能眼镜等)或其他类型的移动设备，或者诸如台式计算机或PC的静止电子设备；电子设备还可以是移动式或静止式的服务器。

计算机指令包括计算机程序代码，计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。

本发明提供一种计算机可读存储介质，其存储有可由电子设备执行的计算机程序，当程序在电子设备上运行时，使得电子设备执行上述波形一致性转换的微震偏移叠加定位方法的步骤。

该存储介质的技术方案与上述波形一致性转换的微震偏移叠加定位方法的技术方案属于同一构思，存储介质的技术方案未详细描述的细节内容，均可以参见上述波形一致性转换的微震偏移叠加定位方法的技术方案的描述。

存储介质可以包括：能够携带计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、电载波信号、电信信号以及软件分发介质等。

本发明的优点为：

1、将极性反转的微震事件波形转换成波形一致性的波形数据：现有技术中，初至极性反转现象会导致常规波形叠加成像函数的失效，使得本该聚集在真实震源处的成像能量出现衰减和分散，造成微地震事件的探测失败和定位误差增大。本发明通过利用各地震道上微震事件信号的相位差为零或者±π的特性，利用原波形的归一化解析形式对原波形进行相乘加权，得到波形一致性的转换波形，可以将极性反转引起的分散成像能量归位到震源位置处，提高成像能量值，提升定位精度。

2、转换波形中的噪音成分仍然保存着不相干性特征：波形转换后的噪音成分仍具有不相干性特征，可以压制成像域中的噪音能量。从而增大微震事件的探测率，提高成像能量分辨率，提升定位精度。

实例1：模拟数据

如图2中(a)所示，先合成含噪音的两个极性相反的微地震波形其中黑色竖直虚线内为微震事件波形，外侧为噪音波形。利用希尔伯特变换得到所对应的瞬时相位结果，如图2中(b)所示。将两道微震波形的瞬时相位进行相减，得到两个地震道波形的相位差值，如图2中(c)所示，在微震事件时间段内，相位差具有±π的特性，而噪音波形时段内的相位差杂乱无章。

如图3所示，利用本发明提供的方法后，原始波形变成了波形具有一致性的转换波形，实现了极性一致的校正。

实例2：实际数据

微地震信号因能量低，易被各种复杂噪音所淹没，基于波形叠加的偏移成像定位方法能够自动探测和定位低信噪比微地震事件。但是，由于微地震震源机制中所包含的双力偶成分会使得检波器所接收到的波形出现初至极性反转现象，导致常规偏移成像定位方法失效。图4为某地区的实际煤层气水力压裂的地面监测系统及P波速度模型，其中压裂井口的相对位置位于(0m，0m，0m)，共布设了728个采样率为1ms的检波器来记录压裂过程中所诱发的微地震事件信号。

图5中(a)为5s长的实际微震波形记录，其中有经过本发明方法的波形一致性转换后的波形如图5中(b)所示，从图5中(c)和(d)的放大细节对比可知，图5中(c)箭头所指示的初至极性反转波形经过波形一致性转换后，波形的初至极性变为一致(图5中(d))。然后分别将图5中(a)和(b)的原波形数据和转换波形数据代入到步骤6的成像函数中，并利用步骤7得到常规方法和本发明方法的微震事件归一化探测能量随发震时刻变化的曲线，如图7所示。本发明方法的探测能量明显大于常规方法的检测能量，这有利于增大微震事件检测的准确性。

从图6可知，本发明方法共检测出4个明显微地震事件E1-E4和2个疑似事件S1-S2，根据所检测到的微地震事件的发震时刻和发震位置计算初至到时，并投影在原波形记录上，如图7所示。可以看到预测的初至到时曲线较为完美地追踪了微地震事件的初至。这说明了本发明方法在探测和定位低信噪比且初至极性反转的微震事件的有效性。

再将常规方法(图8中(a))和本发明方法(图8中(b))在E4微震事件的定位成像结果进行对比，可以直观地观察到本发明方法可以将常规方法所导致的分散能量重新聚集于震源位置处，且同时压制了非震源处的成像噪音，得到成像分辨率更高的定位成像结果。

图9为利用常规方法和本发明方法对90s长的实际微震记录中的微地震事件进行检测和定位的结果对比。根据由2倍的前0～7秒背景能量所得到的检测能量阈值，本发明方法可以检测到27个微地震事件(图9中(b))，而常规方法只检测到了19个微地震事件(图9中(a))。并通过对比各自检测到的微地震事件发震位置可得，本发明方法的微地震事件位置更加聚集在井口附近(图9中(c))，且在深度方向更集中在1200～1400m深度的煤层气地层位置，反映了微震事件位置更加合理。以上测试结果表明本发明的定位方法具有更出色的微震事件检测能力和更高的定位精度。

以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。