CN115951403A - 微震偏移叠加定位方法、装置、电子设备及可读存储介质 - Google Patents
微震偏移叠加定位方法、装置、电子设备及可读存储介质 Download PDFInfo
- Publication number
- CN115951403A CN115951403A CN202211711621.1A CN202211711621A CN115951403A CN 115951403 A CN115951403 A CN 115951403A CN 202211711621 A CN202211711621 A CN 202211711621A CN 115951403 A CN115951403 A CN 115951403A
- Authority
- CN
- China
- Prior art keywords
- waveform
- microseismic
- imaging
- filtered
- origin
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 52
- 238000003384 imaging method Methods 0.000 claims abstract description 107
- 238000001514 detection method Methods 0.000 claims abstract description 35
- 230000010355 oscillation Effects 0.000 claims abstract description 16
- 238000006243 chemical reaction Methods 0.000 claims description 30
- 238000001914 filtration Methods 0.000 claims description 19
- 238000010606 normalization Methods 0.000 claims description 16
- 238000004458 analytical method Methods 0.000 claims description 15
- 238000004891 communication Methods 0.000 claims description 15
- 230000009466 transformation Effects 0.000 claims description 14
- 230000006870 function Effects 0.000 claims description 13
- 238000012545 processing Methods 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 11
- 238000004590 computer program Methods 0.000 claims description 8
- 238000005070 sampling Methods 0.000 claims description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 6
- 238000007781 pre-processing Methods 0.000 claims description 6
- 230000004807 localization Effects 0.000 claims description 5
- AYFVYJQAPQTCCC-GBXIJSLDSA-N L-threonine Chemical compound C[C@@H](O)[C@H](N)C(O)=O AYFVYJQAPQTCCC-GBXIJSLDSA-N 0.000 claims description 4
- 238000013508 migration Methods 0.000 description 8
- 230000005012 migration Effects 0.000 description 8
- 230000006872 improvement Effects 0.000 description 7
- 238000012544 monitoring process Methods 0.000 description 6
- 238000007796 conventional method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 230000007246 mechanism Effects 0.000 description 4
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 2
- 239000003245 coal Substances 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000004984 smart glass Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开一种微震偏移叠加定位方法、装置、电子设备及可读存储介质，包括：利用希尔伯特变换得到滤波后波形的解析形式，并归一化；将滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式，取其实部再减去滤波后波形的包络，得到输入波形数据；计算各个网格点到所有地震道的初至旅行时间；设定发震时刻间隔，将输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体；取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线，根据微震事件检测曲线及成像能量阈值得到定位结果。本发明既校正波形极性，又保留噪音的不相干性，能提高定位精度。
Description
技术领域
本发明属于微震监测资料处理技术领域，具体涉及一种微震偏移叠加定位方法、装置、电子设备及可读存储介质。
背景技术
近年来，基于偏移成像的定位方法因无需拾取初至到时而得到广泛应用。该基于偏移成像的定位方法的思路是：在整个时间、空间范围内寻找可能发生的微震震源。这是一种新的可对震源分布进行成像的方法，利用监测数据的振幅和到时等波形数据，通过系统的扫描某个假设的震源位置和起始时刻，找出整个地震序列分布。在不用精确拾取到时和计算理论地震图的情况下，达到比较理想的定位结果。随着微震监测技术需求的发展，这类定位方法发展迅速，已经成为了微震监测中主要的定位方法之一。
天然地震及诱发(微)地震事件的震源机制中往往包含双力偶成分，使得布设在地表的检波器所接收到的微地震事件波形出现极性反转现象，即初至极性方向具有正负之分。这导致了常规波形叠加成像函数的失效，使得本该聚集在真实震源处的成像能量发生衰减和分散，造成微地震事件的探测失败和定位误差增大。目前提出的三种解决策略：一是通过反演震源机制来校正极性，但是计算成本高，且震源机制与震源位置耦合性强，不确定性大；二是利用聚焦算子将分散的成像能量归位到真实震源处，但要求规则且密集的台网布设；三是通过特征函数将波形数据转化成正值序列的波形数据，该策略计算成本低且不受台网布设的限制，但却失去了噪音的不相干性特征，导致成像域噪音增多，影响定位精度。
因此，有必要研究一种合适的微震偏移叠加定位方法，使其极性反转的微地震信号波形都一致，且同时保存噪音的不相干性特征。以此达到既增强微震信号叠加的成像能量，将分散的成像能量重新归位到真实震源位置；又能借助噪音的不相干性压制成像域噪音，提高成像能量分辨率，达到提升定位精度的目的。
发明内容
为了解决上述技术问题，本发明提供了一种微震偏移叠加定位方法、装置、电子设备及可读存储介质。
本发明公开一种波形一致性转换的微震偏移叠加定位方法，包括：
步骤1、对实际微震波形数据进行常规带通滤波预处理；
步骤2、利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；
步骤3、将滤波后波形的归一化解析形式对所述滤波后波形进行相乘，得到转换波形的解析形式；
步骤4、将所述转换波形的解析形式的实部减去所述滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；
步骤5、按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算所述三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；
步骤6、设定发震时刻间隔，利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，其中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三维空间的能量成像体；
步骤7、获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；
步骤8、设定成像能量阈值，所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。
作为本发明的进一步改进，步骤2中利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理，包括：
希尔伯特变换公式为以下公式1：
公式1：
其中，H[s(t)]为滤波后波形s(t)的希尔伯特变换，将其作为滤波后波形解析形式的虚部，s(t)作为解析形式的实部，得到各地震道滤波后波形的解析形式的第一种表达方式为以下公式2：
公式2：sa(t)＝s(t)+iH[s(t)]
其中，sa(t)为滤波后波形的解析形式，i为虚数单位；
各地震道滤波后波形的解析形式的第二种表达方式为以下公式3：
公式3：sa(t)＝A(t)eiφ(t)或sa(t)＝A(t)cos(φ(t))+iA(t)sin(φ(t))
其中，
为波形的包络或者瞬时振幅；
对各地震道滤波后波形的解析形式sa(t)进行归一化处理，归一化计算公式为以下公式4：
公式4：
其中，Norm(sa(t))为各地震道滤波后波形的归一化解析形式。
作为本发明的进一步改进，步骤3中将各地震道滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式，包括：
根据以下公式5得到转换波形的解析形式：
公式5：xa(t)＝2s(t)*Norm(sa(t))
＝2A(t)cos(φ(t))*eiφ(t)
＝A(t)[cos(2φ(t))+1]+iA(t)sin(2φ(t))
其中，xa(t)为转换波形的解析形式。
作为本发明的进一步改进，步骤4中将转换波形解析形式的实部减去滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据，包括：
根据以下公式6得到波形一致性转换后的波形数据：
公式6：T(t)＝real(xa(t))-A(t)
＝A(t)[cos(2φ(t))+1]-A(t)
＝A(t)cos(2φ(t))
其中，T(t)为滤波后波形s(t)经过波形一致性转换后的波形数据，即作为后续偏移成像定位的输入波形数据。
作为本发明的进一步改进，步骤6中利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，包括：
根据以下公式7得到所述四维成像能量体；
公式7：
其中，
为网格点rk(x，y，z)在第v个发震时刻
的成像能量，tkn为网格点k到检波器n的初至旅行时，W为设置的包含P或S波形时窗的采样点个数，Δt为采样率，N为检波器数目。
作为本发明的进一步改进，步骤7中获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线，包括：
根据以下公式78计算微震事件检测曲线；
公式8：
其中，
为在第v个发震时刻
的最大成像能量，及其所对应的网格点位置r0。
作为本发明的进一步改进，步骤8中所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，包括：
若设定成像能量阈值THRE，则满足
的发震时刻点所包含的信息为成功检测出来的微震事件的发震时刻和发震位置。
本发明还公开了一种波形一致性转换的微震偏移叠加定位的装置，包括：
预处理模块，用于对实际微震波形数据进行常规带通滤波预处理；
变换模块，用于利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；
相乘模块，用于将滤波后波形的归一化解析形式对所述滤波后波形进行相乘，得到转换波形的解析形式；
相减模块，用于将所述转换波形的解析形式的实部减去所述滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；
计算模块，用于按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算所述三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；
叠加模块，用于设定发震时刻间隔，利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，其中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三维空间的能量成像体；
获取模块，用于获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；
定位模块，用于设定成像能量阈值，所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。
本发明还公开了一种电子设备，包括：处理器、通信接口、存储器和通信总线，其中，处理器、通信接口、存储器通过通信总线完成相互间的通信；
所述存储器中存储有计算机程序，当所述程序被所述处理器执行时，使得所述处理器执行上述波形一致性转换的微震偏移叠加定位方法的步骤。
本发明还公开了一种存储介质，其存储有可由电子设备执行的计算机程序，当所述程序在所述电子设备上运行时，使得所述电子设备执行上述波形一致性转换的微震偏移叠加定位方法的步骤。
与现有技术相比，本发明的有益效果为：
利用各地震道上微震事件信号的相位差为零或者±π的特性，将原波形的归一化解析形式对原波形进行相乘加权，可以得到波形一致性的转换波形，可以将极性反转引起的分散成像能量归位到震源位置处。此外，还可以使得转换波形中的噪音成分仍然保存着不相干性特征，从而压制成像域中的噪音能量，提高微震事件的探测率和成像能量分辨率，提升定位精度。
附图说明
为了更清楚地说明本申请的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对本申请保护范围的限定。在各个附图中，类似的构成部分采用类似的编号。
图1示出了本申请实施例提供的波形一致性转换的微震偏移成像定位方法的流程示意图；
图2示出了本申请实施例提供的极性相反的微震信号的相位差原理图；
图3示出了本申请实施例提供的原始波形及波形一致性转换的波形数据的对比示意图；
图4示出了本申请实施例提供的实际煤层气水力压裂的地面监测系统及P波速度模型；
图5示出了本申请实施例提供的实际原始微震波形记录及波形一致性转换后的波形记录对比图；
图6示出了本申请实施例提供的微震波形记录进行微地震事件检测结果对比图；
图7示出了本申请实施例提供的利用所探测到的微震事件的发震时刻及发震位置所预测的初至到时曲线在原始波形记录上的投影图；
图8示出了本申请实施例提供的E4微震事件进行定位成像的对比图；
图9示出了本申请实施例提供的90秒长的实际微震记录进行微地震事件检测及定位结果对比图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明提供一种波形一致性转换的微震偏移叠加定位方法、装置、电子设备及可读存储介质，首先对所有台站的微震波形数据进行常规带通滤波预处理；利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；然后，将所述滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式；将所述转换波形解析形式的实部减去滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；将定位成像区域进行三维网格划分，计算每个网格到所有检波器的初至旅行时间，得到走时表数据；设定合理的发震时刻间隔，将输入波形数据根据走时表偏移叠加得到四维成像能量体；取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；最后，设定成像能量阈值，微震事件探测曲线中超过阈值的发震时刻信息，即为成功检测出来的微震事件的发震时刻和发震位置。
下面结合附图对本发明做进一步的详细描述：
如图1所示，本发明提供一种波形一致性转换的微震偏移叠加定位方法，包括：
步骤1、对实际微震波形数据进行常规带通滤波预处理；
步骤2、利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；其中，
希尔伯特变换公式为：
其中，H[s(t)]为滤波后波形s(t)的希尔伯特变换，将其作为滤波后波形解析形式的虚部，s(t)作为解析形式的实部，即得到各地震道滤波后波形的解析形式：
sa(t)＝s(t)+iH[s(t)]
其中，sa(t)为滤波后波形的解析形式，i为虚数单位。
上式也可以表示为：
sa(t)＝A(t)eiφ(t)或sa(t)＝A(t)cos(φ(t))+iA(t)sin(φ(t))
其中，
为波形的包络，或者瞬时振幅；
对各地震道滤波后波形的解析形式sa(t)进行归一化处理，归一化计算公示为：
其中Norm(sa(t))为各地震道滤波后波形的归一化解析形式。
步骤3、将滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式。计算公式为：
xa(t)＝2s(t)*Norm(sa(t))
＝2A(t)cos(φ(t))*eiφ(t)
＝A(t)[cos(2φ(t))+1]+iA(t)sin(2φ(t))
其中，xa(t)为转换波形的解析形式。
步骤4、将转换波形解析形式的实部减去滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；
计算公式如下：
T(t)＝real(xa(t))-A(t)
＝A(t)[cos(2φ(t))+1]-A(t)
＝A(t)cos(2φ(t))
其中，T(t)为滤波后波形s(t)经过波形一致性转换后的波形数据，即作为后续偏移成像定位的输入波形数据。
步骤5、按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；
步骤6、设定合理的发震时刻间隔，利用成像函数将输入波形数据沿着各个网格点到所有检波器的初至到时(发震时刻+初至旅行时间)曲线进行偏移叠加，得到四维(发震时刻+三维空间)成像能量体。成像函数的计算公式为：
其中，
为网格点rk(x，y，z)在第v个发震时刻
的成像能量；tkn为网格点k到检波器n的初至旅行时；W为设置的包含P或S波形时窗的采样点个数；Δt为采样率，N为检波器数目。
在本实施例中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三位空间的成像能量体。
步骤7、获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线。微震事件检测曲线的计算公式为：
其中，
为在第v个发震时刻
的最大成像能量，及其所对应的网格点位置r0。
步骤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深度的煤层气地层位置，反映了微震事件位置更加合理。以上测试结果表明本发明的定位方法具有更出色的微震事件检测能力和更高的定位精度。
以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种波形一致性转换的微震偏移成像定位方法，其特征在于，包括：
步骤1、对实际微震波形数据进行常规带通滤波预处理；
步骤2、利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；
步骤3、将滤波后波形的归一化解析形式对所述滤波后波形进行相乘，得到转换波形的解析形式；
步骤4、将所述转换波形的解析形式的实部减去所述滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；
步骤5、按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算所述三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；
步骤6、设定发震时刻间隔，利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，其中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三维空间的能量成像体；
步骤7、获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；
步骤8、设定成像能量阈值，所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。
2.根据权利要求1所述的波形一致性转换的微震偏移叠加定位方法，其特征在于，步骤2中利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理，包括：
希尔伯特变换公式为以下公式1：
公式1：
其中，H[s(t)]为滤波后波形s(t)的希尔伯特变换，将其作为滤波后波形解析形式的虚部，s(t)作为解析形式的实部，得到各地震道滤波后波形的解析形式的第一种表达方式为以下公式2：
公式2：sa(t)＝s(t)+iH[s(t)]
其中，sa(t)为滤波后波形的解析形式，i为虚数单位；
各地震道滤波后波形的解析形式的第二种表达方式为以下公式3：
公式3：sa(t)＝A(t)eiφ(t)或sa(t)＝A(t)cos(φ(t))+iA(t)sin(φ(t))
其中，
为波形的包络或者瞬时振幅；
为波形的瞬时相位；
对各地震道滤波后波形的解析形式sa(t)进行归一化处理，归一化计算公式为以下公式4：
公式4：
其中，Norm(sa(t))为各地震道滤波后波形的归一化解析形式。
3.根据权利要求1所述的波形一致性转换的微震偏移叠加定位方法，其特征在于，步骤3中将各地震道滤波后波形的归一化解析形式对滤波后波形进行相乘，得到转换波形的解析形式，包括：
根据以下公式5得到转换波形的解析形式：
其中，xa(t)为转换波形的解析形式。
4.根据权利要求1所述的波形一致性转换的微震偏移叠加定位方法，其特征在于，步骤4中将转换波形解析形式的实部减去滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据，包括：
根据以下公式6得到波形一致性转换后的波形数据：
公式6：T(t)＝real(xa(t))-A(t)
＝A(t)[cos(2φ(t))+1]-A(t)
＝A(t)cos(2φ(t))
其中，T(t)为滤波后波形s(t)经过波形一致性转换后的波形数据，即作为后续偏移成像定位的输入波形数据。
5.根据权利要求1所述的波形一致性转换的微震偏移叠加定位方法，其特征在于，步骤6中利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，包括：
根据以下公式7得到所述四维成像能量体；
公式7：
其中，
为网格点rk(x，y，z)在第v个发震时刻
的成像能量，tkn为网格点k到检波器n的初至旅行时，W为设置的包含P或S波形时窗的采样点个数，Δt为采样率，N为检波器数目。
6.根据权利要求1所述的波形一致性转换的微震偏移叠加定位方法，其特征在于，步骤7中获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线，包括：
根据以下公式78计算微震事件检测曲线；
公式8：
其中，
为在第v个发震时刻
的最大成像能量，及其所对应的网格点位置r0。
7.根据权利要求6所述的波形一致性转换的微震偏移叠加定位方法，其特征在于，步骤8中所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，包括：
若设定成像能量阈值THRE，则满足
的发震时刻点所包含的信息为成功检测出来的微震事件的发震时刻和发震位置。
8.一种波形一致性转换的微震偏移叠加定位装置，其特征在于，所述装置包括：
预处理模块，用于对实际微震波形数据进行常规带通滤波预处理；
变换模块，用于利用希尔伯特变换得到各地震道滤波后波形的解析形式，并进行归一化处理；
相乘模块，用于将滤波后波形的归一化解析形式对所述滤波后波形进行相乘，得到转换波形的解析形式；
相减模块，用于将所述转换波形的解析形式的实部减去所述滤波后波形的包络，得到波形一致性转换后的波形数据，作为输入波形数据；
计算模块，用于按预设网格间距将定位成像区域划分为三维网格体，并利用射线追踪算法计算所述三维网格体中各个网格点到所有地震道的初至旅行时间，得到走时表数据；
叠加模块，用于设定发震时刻间隔，利用成像函数将所述输入波形数据沿着各个网格点到所有检波器的初至到时曲线进行偏移叠加，得到四维成像能量体，其中，所述初至到时曲线根据发震时刻和所述初至旅行时间确定得到，所述四维能量成像体为表征发震时刻及三维空间的能量成像体；
获取模块，用于获取各个发震时刻的成像能量体中的最大值及所对应的网格点位置，得到各个发震时刻的微震事件检测曲线；
定位模块，用于设定成像能量阈值，所述微震事件探测曲线中超过所述成像能量阈值的发震时刻点所包含的信息即为成功检测出来的微震事件的发震时刻和发震位置，得到定位结果。
9.一种电子设备，其特征在于，包括：处理器、通信接口、存储器和通信总线，其中，处理器、通信接口、存储器通过通信总线完成相互间的通信；
所述存储器存储有计算机程序，当所述程序被所述处理器执行时，使得所述处理器执行如权利要求1至7中任一项所述的波形一致性转换的微震偏移叠加定位方法的步骤。
10.一种计算机可读存储介质，其特征在于，其存储有可由电子设备执行的计算机程序，当所述程序在所述电子设备上运行时，使得所述电子设备执行如权利要求1至7任一项所述的波形一致性转换的微震偏移叠加定位方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211711621.1A CN115951403B (zh) | 2022-12-29 | 2022-12-29 | 微震偏移叠加定位方法、装置、电子设备及可读存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211711621.1A CN115951403B (zh) | 2022-12-29 | 2022-12-29 | 微震偏移叠加定位方法、装置、电子设备及可读存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115951403A true CN115951403A (zh) | 2023-04-11 |
| CN115951403B CN115951403B (zh) | 2023-08-18 |
Family
ID=87287365
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211711621.1A Active CN115951403B (zh) | 2022-12-29 | 2022-12-29 | 微震偏移叠加定位方法、装置、电子设备及可读存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115951403B (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102928873A (zh) * | 2012-10-30 | 2013-02-13 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 基于四维能量聚焦的地面微地震定位方法 |
| US20150316669A1 (en) * | 2012-11-21 | 2015-11-05 | Westerngeco Seismic Holdings Limited | Processing microseismic data |
| US20170097430A1 (en) * | 2014-04-28 | 2017-04-06 | Microseismic, Inc. | Method of using semblance of corrected amplitudes due to source mechanisms for microseismic event detection and location |
| CN107179551A (zh) * | 2017-06-19 | 2017-09-19 | 吉林大学 | 一种利用微震记录对地下构造直接成像的方法 |
| CN110389377A (zh) * | 2018-04-23 | 2019-10-29 | 中国海洋大学 | 基于波形互相关系数相乘的微震偏移成像定位方法 |
| CN114415231A (zh) * | 2021-12-22 | 2022-04-29 | 南方科技大学 | 一种基于台站对edt面概率分布函数的微震定位方法 |
-
2022
- 2022-12-29 CN CN202211711621.1A patent/CN115951403B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102928873A (zh) * | 2012-10-30 | 2013-02-13 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 基于四维能量聚焦的地面微地震定位方法 |
| US20150316669A1 (en) * | 2012-11-21 | 2015-11-05 | Westerngeco Seismic Holdings Limited | Processing microseismic data |
| US20170097430A1 (en) * | 2014-04-28 | 2017-04-06 | Microseismic, Inc. | Method of using semblance of corrected amplitudes due to source mechanisms for microseismic event detection and location |
| CN107179551A (zh) * | 2017-06-19 | 2017-09-19 | 吉林大学 | 一种利用微震记录对地下构造直接成像的方法 |
| CN110389377A (zh) * | 2018-04-23 | 2019-10-29 | 中国海洋大学 | 基于波形互相关系数相乘的微震偏移成像定位方法 |
| CN114415231A (zh) * | 2021-12-22 | 2022-04-29 | 南方科技大学 | 一种基于台站对edt面概率分布函数的微震定位方法 |
Non-Patent Citations (2)
| Title |
|---|
| 曾志毅;张建中;: "利用微地震记录互相关成像的震源定位方法", 石油地球物理勘探, no. 02 * |
| 毕丽飞;曾志毅;张建中;黄忠来;芮拥军;刁瑞;: "一种适于存在极性反转的微震初至到时拾取方法", 石油物探, no. 03 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115951403B (zh) | 2023-08-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Stachnik et al. | Determination of New Zealand ocean bottom seismometer orientation via Rayleigh-wave polarization | |
| Tian et al. | Variable-eccentricity hyperbolic-trace TFPF for seismic random noise attenuation | |
| Rueda et al. | Near-real-time seismic moment-tensor determination in Spain | |
| Baker et al. | Real-time earthquake location using Kirchhoff reconstruction | |
| Hino et al. | Was the 2011 Tohoku-Oki earthquake preceded by aseismic preslip? Examination of seafloor vertical deformation data near the epicenter | |
| US20120041682A1 (en) | Attenuating internal multiples from seismic data | |
| CN107193045B (zh) | 一种地震数据处理方法及装置 | |
| Saad et al. | Automatic arrival time detection for earthquakes based on stacked denoising autoencoder | |
| García et al. | Advances on the automatic estimation of the P-wave onset time. | |
| CN113625337A (zh) | 一种极浅水高精度地震资料快速成像方法 | |
| CN104570116A (zh) | 基于地质标志层的时差分析校正方法 | |
| CN114994754A (zh) | 基于直达波和深度震相初动极性的震源机制联合反演方法 | |
| Yoon et al. | Unsupervised large‐scale search for similar earthquake signals | |
| Ruppert et al. | Enhanced regional earthquake catalog with Alaska amphibious community seismic experiment data | |
| CN111352153B (zh) | 一种基于瞬时相位互相关加权的微地震干涉定位方法 | |
| CN115951403A (zh) | 微震偏移叠加定位方法、装置、电子设备及可读存储介质 | |
| Wang et al. | Automatic onset phase picking for portable seismic array observation | |
| FR2985818A1 (fr) | Dispositif et procede d'estimation de decalages temporels | |
| CN110780346A (zh) | 一种隧道超前探测复杂地震波场的分离方法 | |
| CN109633744B (zh) | 地震子波的提取方法、装置、设备及存储介质 | |
| Yin et al. | Research on interference signal recognition in p wave pickup and magnitude estimation | |
| Yao et al. | Microseismic signal denoising using simple bandpass filtering based on normal time–frequency transform | |
| Kemerait et al. | Quality control metrics for cepstral analysis with homomorphic deconvolution | |
| CN117687077B (zh) | 利用das和检波器混合阵列监测微震的方法和系统 | |
| CN103698812A (zh) | 利用叠前地震道集计算地层品质因数的方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |