CN104765064B - 一种微地震干涉成像的方法 - Google Patents

Abstract

本发明涉及一种微地震干涉成像的方法，所述方法包括：采用一定大小的网格，将地下模型进行离散，并在地下模型中设置n个检波器；建立地下速度模型，根据速度模型，分别计算每个网格点到每个检波器的初至波的走时；分别获取微地震的主事件在n个检波器中的第一记录信息以及目标事件在n个检波器中的第二记录信息；根据第一记录信息和第二记录信息，对微地震的主事件和目标事件进行互相关计算，并求包络，获取互相关包络道集；在每一个网格点，按照特定的时间步长，设定不同的激发时刻值，对互相关包络道集进行加时窗干涉成像处理；对处理后的图像进行剖面处理，获取不同剖面图像中的最大值；对不同剖面图像的最大值进行比对，获取目标最大值。

Description

一种微地震干涉成像的方法

技术领域

本发明涉及地球物理微地震监测定位技术领域，尤其涉及一种微地震干涉成像的方法。

背景技术

水力压裂是目前实现非常规油气藏经济开发的关键技术之一。微地震监测技术是实时水力压裂效果监测最有效的方法之一。通过对水力压裂等工程作业诱发的微弱地震进行监测，可以揭示压裂裂缝的空间尺寸及延展形态等，有助于优化压裂的施工设计和效果。

传统的井中微地震定位方法，需要进行有效事件的初至拾取，而这种方法可能存在拾取误差，从而产生定位误差，并且效率相对较低。地面微地震监测方法则由于信噪比极低，一般无法直接对微地震记录进行准确的初至拾取。近年来，人们借鉴地震勘探中的偏移方法，发展了一些无需拾取初至波的走时的成像类微地震定位方法。该类基于能量聚集或扫描的方法由于没有考虑震源机制对波形幅度及相位的影响，可能造成定位精度较低的情况。如逆时成像和绕射叠加分别对方向力源和双力偶源成像时，成像精度都明显受到震源机制的影响。

发明内容

针对常规的能量聚集或扫描方法在微地震资料处理中可能受震源机制影响的问题，本发明提供了一种进行震源定位的微地震干涉成像方法，通过引入主事件，提高了定位分辨率，并对包含走时差信息的互相关道集求包络来解决微地震记录中的极性变化问题，以提高微地震定位精度，为水力压裂裂缝监测提供可靠的依据和指导。

本发明实施例提供了一种微地震干涉成像的方法，其特征在于，所述方法包括：

采用一定大小的网格，将地下模型进行离散，并在地下模型中设置n个检波器，n为自然数；

建立地下速度模型，根据所述速度模型，分别计算每个网格点到所述n个检波器的每个检波器的初至波的走时，获取走时表；

分别获取微地震的主事件在n个检波器中的第一记录信息，其中，第一记录信息包括微地震的主事件的初至波的走时以及波形；

分别获取微地震的目标事件在n个检波器中的第二记录信息，其中，第二记录信息包括目标事件的波形；

根据第一记录信息以及第二记录信息，对微地震的主事件和目标事件进行互相关计算，获取互相关值；计算与互相关值相对应的包络值，获取互相关包络道集，互相关包络道集包括主事件和目标事件之间的初至波的走时差值和目标事件的激发时刻信息；

在每一个网格点，按照特定的时间步长，设定不同的激发时刻值，对互相关包络道集进行加时窗干涉成像处理；

对加时窗干涉成像处理后的图像进行剖面处理，获取不同剖面图像中的最大值；

对不同剖面图像的最大值进行比对，获取目标最大值。

优选的，目标最大值对应的时间和位置为微地震的目标事件的激发时刻和位置。

优选的，分别获取微地震的主事件在n个检波器中的第一记录信息具体包括：获取微地震的主事件在n个检波器的每个检波器中的初至波的走时以及波形，初至波的走时包括直达纵波的走时和直达横波的走时。

优选的，根据第一记录信息以及第二记录信息，对微地震的主事件和目标事件进行互相关计算，获取互相关值具体包括：

计算互相关值如下式所示：

其中，C(tc)为互相关系数，x(t)为第一记录信息中微地震的主事件波形的时间离散序列，y(t)为第二记录信息中微地震目标事件波形的时间离散序列，tc为微地震目标事件波形的时间离散序列的时间轴移位值，y(t-tc)为微地震目标事件波形的时间离散序列y(t)经过移位长度tc后的离散序列，t∈[1,N]，tc∈[-N,N]，N为每一个检波器中的地震波波形的时间长度。

优选的，计算与互相关值相对应的包络值具体包括：

计算包络值如下式所示：

其中，E(t)为微地震的主事件和目标事件进行互相关计算后，获取的与互相关值C(t)相对应的包络值，H{C(t)}为互相关值C(t)的希尔伯特变换。

优选的，对互相关包络道集进行加时窗干涉成像处理具体包括：

进行加时窗干涉成像处理具体如下式所示：

其中，M(x,t₀)为不同激发时刻t₀对应的干涉偏移成像值，x为网格点的位置向量，i为目标事件中的第i个检波器，为主事件与目标事件的互相关包络道集，t_w为[-W,W]大小的时窗，i^m为主事件中的第i个检波器，m代表主事件的位置向量，τ_im和τ_ix分别为主事件的第i个检波器中所检测到的初至波的走时和目标事件到第i个检波器的理论初至波的走时，t₀为设定的不同的激发时刻值。

本发明提供的微地震干涉成像的方法，无需对微地震的目标事件的初至波进行拾取，因此，对震源的定位效率明显提高；根据主事件和目标事件在检波器中的记录信息进行互相关运算和求包络运算，而后进行加时窗干涉成像处理。在互相关运算和加时窗干涉成像处理过程中均应用了叠加处理，从而更加方便适应低信噪比数据；通过引入主事件，大大增加了定位分辨率；同时进行求包络运算，有效的克服不同震源机制导致波形极性变化对定位结果的影响，进一步提高了定位精度，对低信噪比的实时微地震监测有广阔的应用前景。

附图说明

图1为本发明实施例提供的微地震干涉成像方法的流程示意图；

图2为速度模型及震源-检波器布置示意图；

图3为模拟的主事件记录波形图；

图4为模拟的目标事件记录波形图；

图5为单脉冲互相关示意图；

图6为互相关包络道集波形图；

图7为部分不同激发时刻对应的成像结果对比图；

图8为不同激发时刻对应的目标最大值的曲线示意图；

图9为微地震震源位置示意图。

具体实施方式

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

图1为本发明实施例提供的一种微地震干涉成像方法流程示意图。如图1所示，本发明实施例具体包括以下步骤：

步骤101，采用一定大小的网格，将地下模型进行离散，并在地下模型中设置n个检波器。

具体的，根据地下目标监测区域的实际情况，将水力压裂周围的待监测区域建立一个地下模型，利用一定大小的网格，将地下模型进行离散。并且，在地下模型中，按照一定的规律，在不同的位置设置n个检波器。其中，网格的尺寸根据实际需要进行设定，网格越小，定位精度越高，然而计算时间越长。

步骤102，建立地下速度模型，根据速度模型，分别计算每个网格点到n个检波器的每个检波器中的初至波的走时，获取走时表。

具体的，建立一个地下速度模型，并且根据地下速度模型，计算地下模型中所有网格点到n个检波器的每一个检波器中的理论初至波的走时，获取一个走时表。其中，计算初至波的走时可以采用射线追踪的方法或者解程函方程的方法。

步骤103，分别获取微地震的主事件在n个检波器中的第一记录信息。

具体的，选取或者模拟一个微地震主事件，将其在n个检波器的每个检波器中所检测到的波形以及初至波的走时提取出来，即获取第一记录信息，其中，初至波的走时包括直达纵波(P波)的走时和直达横波(S波)的走时。

在一个具体的例子中，当有已知的射孔事件或者微地震事件时，则将该事件作为主事件，直接提取该事件在检波器中所检测到的波形和初至波的走时；

否则，则将正演模拟监测区域中特定位置在特定震源机制条件下产生的微地震记录作为主事件，并提取该事件在检波器中所检测到的波形和初至波的走时。

在本实施例中，参见图2所示，图2为速度模型及震源-检波器布置示意图，采用的是二维的水平三层速度模型，模拟井中监测的情况，在模型左侧布置51个间距10m的检波器，模型大小为500m×500m，网格尺寸为2.5m。主事件为正演模拟监测区域中特定位置在特定震源机制条件下产生的微地震记录，并且将其设置在模型中心，即图中星形中标注m符号的位置(位置坐标为：250m，250m)处。图3为模拟的主事件记录波形图，图3所显示的是主事件在各个检波器中所提取的波形图，震源采用垂直方向激发的集中力源，左图为地震波的振动速度在x方向的分量V_x，由图可以看出，该模型下井中记录的特点为：分量V_x的主要极性为负，其中，时间轴200ms前后为起点所对应的波为P波，300ms前后为起点所对应的波为S波，S波能量明显强于P波。右图为地震波的振动速度在z方向的分量V_z，由图可以看出，该分量V_z的主要极性为正，其中，时间轴200ms前后为起点所对应的波同样为P波，300ms前后为起点所对应的波为S波，S波能量明显强于P波。

步骤104，分别获取微地震的目标事件在n个检波器中的第二记录信息。

具体的，提取目标事件在n个检波器中所记录的n个波形信息。

在本实例中，参见图2和图4所示，在图2中，未标注m符号的星形的位置则代表目标事件的位置(位置坐标为：350m，325m)，预设的激发时刻为15ms。图4为模拟的目标事件记录波形图，考虑到实际微地震不同震源机制对定位结果的影响，所以目标事件的震源并不是采用与主事件相同的垂直方向激发的集中力源，而是采用水平方向激发的集中力源。由图中可以看出，该模型下井中记录的特点为：左图为地震波振动速度在x方向的分量V_x，V_x的主要极性为正，其中，起点在时间轴200ms与300ms之间的波为P波，P波极性为正；起点在时间轴300ms与400ms之间，并且更接近400ms所对应的波为S波，S波极性为负，P波能量明显强于S波。右图为地震波振动速度在z方向的分量V_z，由图可以看出，起点在时间轴200ms与300ms之间的波为P波，其主要极性为正，起点在时间轴300ms与400ms之间，并且更接近400ms所对应的波为S波，其极性为负，P波与S波的能量几乎相等。

步骤105，根据第一记录信息以及第二记录信息，对微地震的主事件和目标事件进行互相关计算，获取互相关值；计算与互相关值相对应的包络值，获取互相关包络道集。

具体的，首先假设目标微地震事件中的记录只有单一事件，即所记录的每个检波器中的波形有两个峰值(若只含P波或S波则只有一个峰值)。然后，对主事件在每个检波器中的第一记录信息和相对应的目标事件在每个检波器中的第二记录信息进行互相关，获取互相关值。计算与互相关值相对应的包络值，获取互相关包络道集。其中，主事件和目标事件在每个检波器中的互相关运算可以由公式(1-1)表示：

C(tc)＝Σ_tx(t)y(t-tc) (1-1)

其中，C(tc)为互相关系数，x(t)为第一记录信息中微地震的主事件波形的时间离散序列，y(t)为第二记录信息中微地震目标事件波形的时间离散序列，tc为微地震目标事件波形的时间离散序列的时间轴移位值，y(t-tc)为微地震目标事件波形的时间离散序列y(t)经过移位长度tc后的离散序列，t∈[1,N]，tc∈[-N,N]，N为每一个检波器中的地震波波形的时间长度。

包络值的计算可由公式(1-2)表示：

其中，E(t)为微地震的主事件和目标事件进行互相关计算后，获取的与互相关值C(t)相对应的包络值，H{C(t)}为互相关值C(t)的希尔伯特变换。

根据对微地震的主事件和目标事件进行互相关计算以及求包络运算后，获取互相关包络道集，其中，互相关包络道集中包含了主事件和目标事件之间的初至波的走时差值和目标事件的激发时刻信息，但是激发时刻信息为未知的。

通过互相关运算，获取主事件中在每个检波器的初至波的走时与目标事件在每个检波器中的初至波的走时之间的实际走时差值，参见图5所示，图5为单脉冲互相关波形示意图。互相关系数在互相关波形图中所对应为波形的峰值，该峰值所对应的时间为T＝ΔT+N，其中ΔT＝t_a-t_b，t_a和t_b分别为不同波形对应的走时，ΔT为走时差值，N为一个检波器中地震波形的时间长度。

而不同位置的检波器在检测同一个地震波中，所检测到的波形幅值可能为正值或者负值，如果直接进行叠加则有可能正负值叠加后相互抵消，从而丢失掉波形中幅值最大值，使得定位精度受到影响。所以通过计算包络值，将波形的幅值全部转换成正值后，再进行叠加，可以消除波形正负极性变化对定位结果的影响。参见图6所示，图6为对主事件和目标事件进行互相关计算和求包络之后所获取的互相关包络道集波形图。

步骤106，在每一个网格点，按照特定的时间步长，设定不同的激发时刻值，对互相关包络道集进行加时窗干涉成像处理。

具体的，根据实际情况，预设一个激发时刻时窗，在每一个网格点，按照特定的时间步长，设定不同的激发时刻值，在每一个激发时刻，对互相关包络道集进行加时窗干涉成像处理，如式(1-3)所示：

其中，M(x,t₀)为不同激发时刻t₀对应的干涉偏移成像值，x为网格点的位置向量，i为目标事件中的第i个检波器，为主事件与目标事件的互相关包络道集，t_w为[-W,W]大小的时窗，i^m为主事件中的第i个检波器，m代表主事件的位置向量，τ_im和τ_ix分别为主事件的第i个检波器中所检测到的初至波的走时和目标事件到第i个检波器的理论初至波的走时，t₀为设定的不同的激发时刻值

假设，地震波P波的传播速度为3km/s，地震波的震源位置到检波器之间的距离为300m-600m，实际地震波的走时为100ms-200ms，但在接收的波形中拾取到P波的起跳时刻约为250ms，则可以预估可能的激发时刻时窗为[250ms-200ms,250ms-100ms]，即[50ms,150ms]。考虑到实际的速度和走时拾取存在一定误差，可以将激发时刻时窗稍加扩大为[30ms,180ms]。激发时刻的步长可以按照实际情况来设定，通常来说，步长越短，精度越高，但是计算量相对的也会越多。

而在本实施例中，参见图7所示，图7为部分不同激发时刻对应的成像结果对比图，在图7中，根据实际情况，预设的激发时刻时窗设定为[0ms,80ms],时间步长为1ms，图中仅仅是选取了激发时刻分别为0ms，15ms，30ms，45ms，60ms以及75ms时的图像。

步骤107，对加时窗干涉成像处理后的图像进行剖面处理，获取不同剖面图像中的最大值。

步骤108，对不同剖面图像的最大值进行比对，获取目标最大值。

具体的，对不同的剖面图像的最大值进行比对，确定不同剖面图像的最大值中的最终最大值作为目标最大值。并且该目标最大值所对应的时间和位置为微地震的目标事件的激发时刻和位置。

参见图7所示，将成像结果进行二维切片显示。获取每一个切片图像中的最大值，并且将每一个切片图像中的最大值进行比对，获取目标最大值。而目标最大值所对应的时间和位置就是微地震的激发时刻和位置。参见图8所示，图8为不同激发时刻对应的目标最大值的曲线示意图，其中最大值对应的时刻为14ms，即求得的激发时刻为14ms，而实际预设的激发时刻为15ms，绝对误差为-1ms。图9为微地震震源位置示意图，在图9中显示了最终微地震干涉成像结果，即激发时刻14ms所对应的成像结果，如图中亮点所示的位置。最大值对应的位置坐标为(350m，322.5m)，即得到的微地震目标事件位置为(350m，322.5m)，绝对误差为(0m，-2.5m)(实际测定的目标事件的位置与预设的位置之间的差值：(350-350)m，(322.5-325)m)。

需要说明的是，在本实施例加时窗干涉成像结果是通过对V_x分量和V_z分量各自成像结果的叠加，而在实际运用中也可以为单独对不同的分量的成像结果进行比较分析，从而减少定位的误差和不确定性。

并且，在实际震源成像的过程中，需要考虑主事件微地震记录的特点，并假设目标事件同时含有纵波和横波，根据互相关记录的特点进行干涉成像。

本发明实施例提供的微地震干涉成像的方法，无需对微地震的目标事件的初至波进行拾取，因此，对震源的定位效率明显提高；根据主事件和目标事件在检波器中的记录信息进行互相关运算和求包络运算，而后进行加时窗干涉成像处理。在互相关运算和加时窗干涉成像处理过程中均应用了叠加处理，从而更加方便适应低信噪比数据；通过引入主事件，大大增加了定位分辨率；同时进行求包括运算，有效的克服不同震源机制导致波形极性变化对定位结果的影响，进一步提高了定位精度，对低信噪比的实时微地震监测有广阔的应用前景。

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

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

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种微地震干涉成像的方法，其特征在于，所述方法包括：

采用一定大小的网格，将地下模型进行离散，并在所述地下模型中设置n个检波器，n为自然数；

建立地下速度模型，根据所述速度模型，分别计算每个网格点到所述n个检波器的每个检波器中的初至波的走时，获取走时表；

分别获取微地震的主事件在所述n个检波器中的第一记录信息，所述第一记录信息包括所述微地震的主事件的初至波的走时以及波形；

分别获取所述微地震的目标事件在所述n个检波器中的第二记录信息，所述第二记录信息包括所述目标事件的波形；

根据所述第一记录信息以及所述第二记录信息，对所述微地震的主事件和目标事件进行互相关计算，获取互相关值；计算与所述互相关值相对应的包络值，获取互相关包络道集，所述互相关包络道集包括所述主事件和目标事件之间的初至波的走时差值和所述目标事件的激发时刻信息；

在每一个所述网格点，按照特定的时间步长，设定不同的激发时刻值，对所述互相关包络道集进行加时窗干涉成像处理；

对所述加时窗干涉成像处理后的图像进行剖面处理，获取不同剖面图像中的最大值；

对所述不同剖面图像的最大值进行比对，获取目标最大值。

2.根据权利要求1所述的方法，其特征在于，所述目标最大值对应的时间和位置为所述微地震的目标事件的激发时刻和位置。

3.根据权利要求1所述的方法，其特征在于，所述分别获取微地震的主事件在所述n个检波器中的第一记录信息具体包括：获取所述微地震的主事件在所述n个检波器的每个检波器中的初至波的走时以及波形，所述初至波的走时包括直达纵波的走时和直达横波的走时。

4.根据权利要求1所述的方法，其特征在于，根据所述第一记录信息以及所述第二记录信息，对所述微地震的主事件和目标事件进行互相关计算，获取互相关值具体包括：

计算所述互相关值如下式所示：

$C\left(tc\right)=\underset{t}{\mathrm{\Σ}}x\left(t\right)y(t-tc)$

其中，C(tc)为互相关系数，x(t)为所述第一记录信息中所述微地震的主事件波形的时间离散序列，y(t)为所述第二记录信息中所述微地震目标事件波形的时间离散序列，tc为所述微地震目标事件波形的时间离散序列的时间轴移位值，y(t-tc)为所述微地震目标事件波形的时间离散序列y(t)经过移位长度tc后的时间离散序列，t∈[1,N]，tc∈[-N,N]，N为每一个所述检波器中的地震波波形的时间长度。

5.根据权利要求1所述的方法，其特征在于，所述计算与所述互相关值相对应的包络值具体包括：

计算所述包络值如下式所示：

$E\left(t\right)=\sqrt{C^2\left(t\right)+H^2\left\{C\left(t\right)\right\}}$

其中，E(t)为所述微地震的主事件和目标事件进行互相关计算后，获取的与互相关值C(t)相对应的包络值，H{C(t)}为所述互相关值C(t)的希尔伯特变换。

6.根据权利要求1所述的方法，其特征在于，所述对所述互相关包络道集进行加时窗干涉成像处理具体包括：

进行加时窗干涉成像处理具体如下式所示：

其中，M(x,t₀)为不同激发时刻t₀对应的干涉偏移成像值，x为所述网格点的位置向量，i为所述目标事件中的第i个检波器，为所述主事件与目标事件的互相关包络道集，t_w为[-W,W]大小的时窗，i^m为所述主事件中的第i个检波器，m代表所述主事件的位置向量，τ_im和τ_ix分别为所述主事件的第i个检波器中所检测到的初至波的走时和所述目标事件到第i个检波器的理论初至波的走时，t₀为设定的不同的激发时刻值。