CN110018517A - 一种多尺度地面微地震逆时干涉定位方法 - Google Patents

Abstract

本发明公开了一种多尺度地面微地震逆时干涉定位方法，我们将逆时干涉定位与多尺度变分模态分解相结合，通过多尺度分解技术将原始微地震数据分解到多个具有不同振动特征的尺度域，随着分解尺度的增加，噪声幅度逐渐降低，有效信号的幅度逐渐增加，就能够降低背景噪声的比重，进而提高震源定位的精度；另外，多尺度变分模态分解能够使原始微地震数据分解到多个尺度域，多个尺度提取的有效信号的特征得到增强，因此能够弥补稀疏观测系统引起的采样信息的不足，有效尺度逆时干涉定位的乘积能够提高逆时干涉定位成像椭圆的聚焦性和分辨率，进而压制成像假象。

Description

一种多尺度地面微地震逆时干涉定位方法

技术领域

本发明涉及一种地面微地震定位方法，具体涉及一种多尺度地面微地震逆时干涉定位方法，属于微地震监测技术领域。

背景技术

微地震监测技术近年来被广泛地应用于致密油气藏的水力压裂监测、矿山岩爆监测、二氧化碳地质封存监测和常规油气田注采诱发地震监测。尤其在致密油气藏的开发中，微地震水力压裂监测是唯一能够对诱导裂缝的几何形状进行成像的远场技术。通过实时微地震监测震源的时间和空间分布，可估算诱导裂缝的几何产状、地应力的变化和裂缝网格的连通性，以此可来估测压裂井的最优射孔区间和储层改造体积，从而达到优化压裂工程的目的。

震源定位是微地震监测的核心技术，震源定位的结果直接影响反演诱导裂缝的空间展布、时域分布及震源机制的精度。常规的震源定位方法需要拾取准确的初至，通过射线追踪最小化理论与拾取初至之差来获取震源的位置。然而，实际采集的微地震数据普遍包含强噪声干扰，所以难以拾取准确可靠的初至，从而降低了这类方法的精度和可靠性。

近年来，适用于低信噪比数据的逆时干涉定位方法(TRI)成为了新的研究热点。逆时干涉定位方法基于弹性波动方程，无需拾取震相走时，同时利用微地震波场的运动学和动力学特征，使逆时传播的波场能量在真实震源处干涉叠加至最大值，避免了人工拾取初至所引入的定位误差。波场互易原理和波动方程在非耗散介质中的时不变性是逆时干涉定位方法的前提。逆时干涉定位方法借鉴逆时偏移原理，以实际采集的波形作为震源子波逆时加载在各个接收器(检波器)，通过有限差分或相移等数值方法进行波场的逆时延拓，重构逆时传播的波场，再施加一定的成像条件即可获得震源的位置。常规的成像条件是：比较不同时刻的波场幅度，使波场幅度达到最大的时刻和位置为推定震源的起震时刻和震源位置。

以二维弹性波方程为例，逆时传播的波场可表示为：

其中，d为实际观测到的波形数据，T为记录时间总长，x，z为空间坐标，λ和μ为拉梅系数，ρ为介质密度，u为逆时外推的波场，震源的位置(x₀，z₀)和起震时刻t₀可由下式得出：

(x₀,z₀,t₀)＝max[u(x,z,t)]

为了得到完整的震源信息，逆时加载的波形应包含整个有效事件。接收器布置在地表(z＝0)不同时刻的逆时波场快照如图1(a)、图1(b)和图1(c)所示。当波场传播到震源的起震时刻(t＝0ms)时，波场能量会聚焦到真实的震源位置。受实际采集条件的限制，震源附近完整的采样信息难以得到，导致逆时定位的最终结果为一个椭圆而非一个点。这个椭圆的半轴长给出了逆时干涉定位的精度，提供了评价逆时干涉定位结果的依据。

然而，微地震逆时干涉定位方法目前还未得到大规模的商业化应用，究其原因主要是：较强的背景噪声和稀疏的观测系统这两个客观缺点影响了逆时延拓波场能量的聚焦性，从而限制了该方法的应用。

(1)背景噪声

实际测量的微地震数据存在很强的背景噪声干扰，使得逆时干涉的波场能量难以聚焦在震源位置。

图2(a)是信噪比＝0.1dB的常规逆时干涉定位结果的俯视图。

图2(b)是信噪比＝0.1dB的常规逆时干涉定位结果的侧视图。

由图2(a)和图2(b)可知，强背景噪声导致逆时干涉定位成像的能量分散，难以聚焦在震源位置，所以难以确定震源位置。

(2)观测系统

实际的微地震观测系统往往比较稀疏，无法提供来自震源附近充足的采样信息，使得逆时干涉定位成像包含较强的成像噪声，难以确定真实震源的位置。

图3(a)是水平层状纵波速度模型。

图3(b)是现有逆时干涉定位结果。

由图3(a)和图3(b)可知，稀疏观测系统导致逆时干涉定位成像产生较强成像噪声，在浅层区域有明显的成像假象，该成像假象在三维复杂介质中更为严重，干扰真实震源位置的确定。

发明内容

为解决现有技术的不足，本发明的目的在于提供一种可以提高震源定位精度的多尺度地面微地震逆时干涉定位方法。

为了实现上述目标，本发明采用如下的技术方案：

一种多尺度地面微地震逆时干涉定位方法，其特征在于，包括以下步骤：

Step 1：读入地面微地震波场记录、观测系统文件、速度模型文件及相关参数文件；

Step 2：采用时频分析技术将原始微地震数据分解到多个尺度域，通过求解带约束的变分优化问题，将原始信号分解成为多个尺度的内蕴模式函数的集合：

其中，u_k和w_k分别是第k个尺度的内蕴模式函数和其对应的中心频率，δ(t)是Dirac函数，t是时间采样点，j是复数，e是自然常数，f(t)是读入的地面微地震波场记录；

Step 3：对分解得到的各尺度内蕴模式函数进行独立成分分析以压制背景噪声；

Step 4：将多尺度分解和独立成分分析噪声压制后的各尺度波形分量加载在检波器位置，进行逆时延拓，以波场模的自相关作为成像条件：

其中，I为波场u的自相关逆时干涉定位成像结果；

Step 5：通过峰值信噪比判断某个尺度的逆时干涉定位成像结果是否有效，当某个尺度逆时定位成像的峰值信噪比大于事先给定的阈值时，判定该尺度为有效尺度；

Step 6：拾取有效尺度的逆时干涉定位成像结果的乘积作为最终的定位成像

其中，n为有效尺度的个数，Iⁱ为第i个有效尺度的波场u的自相关成像。

前述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step1中，前述相关参数文件包括：参与计算地震工区的网格点数、网格间距、采样频率、震源子波参考频率和地震记录时长。

前述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step2中，采用变分模态分解、经验模态分解、S变换或小波变换将原始微地震数据分解到多个尺度域。

前述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step2中，为增加求解过程的稳定性，在公式中加入惩罚因子α和Lagrangian乘法算子λ，将公式改写成：

改写后的公式可用交替方向乘子算法进行求解。

前述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step3中，采用基于负熵最大Fast ICA算法对各尺度内蕴模式函数进行独立成分分析，提取有效信号压制背景噪声。

前述的多尺度地面微地震逆时干涉定位方法，其特征在于，对各尺度内蕴模式函数进行独立成分分析，提取有效信号压制背景噪声，具体方法为：

首先，将观测信号X中心化，使其均值为0；

其次，初始化需要估计的分量个数m，迭代求解白化矩阵W，使得观测信号X在新的子空间的投影为白化向量Z(t)＝WX(t)，W_p的赋值方式为：

其中，p为迭代求解的次数，W_p ^T为W_p的转置，E为取均值计算，g(y)＝tanh(1,y)，g’为g的导数，W₀为随机的初始权矢量，通过下式对W_p进行更新：

若W_p收敛，则W＝W_p，Y(t)＝WX(t)即为去噪后的各尺度内蕴模式函数，否则p＝p+1，重新初始化随机向量W₀，继续迭代求解，若p>m，则W＝W₀。

本发明的有益之处在于：

(1)我们将逆时干涉定位与多尺度变分模态分解相结合，通过多尺度分解技术将原始微地震数据分解到多个具有不同振动特征的尺度域，随着分解尺度的增加，噪声幅度逐渐降低，有效信号的幅度逐渐增加，以此理论为基础，对分解得到的各尺度信号分量施加一定的自适应噪声压制，就能够突出有效信号的能量，降低背景噪声的比重，对多尺度分解后的信号分量进行逆时干涉定位，选取有效尺度的逆时干涉定位结果进行重构，能够降低背景噪声对微地震逆时干涉定位成像图的干扰，进而提高震源定位的精度；

(2)多尺度变分模态分解能够使原始微地震数据分解到多个尺度域，多个尺度提取的有效信号的特征得到增强，因此能够弥补稀疏观测系统引起的采样信息的不足，有效尺度逆时干涉定位的乘积能够提高逆时干涉定位成像椭圆的聚焦性和分辨率，进而压制成像假象。

附图说明

图1(a)至图1(c)分别是接收器布置在地表(z＝0)不同时刻的逆时波场快照；

图2(a)和图2(b)分别是信噪比＝0.1dB的常规逆时干涉定位结果的俯视图和侧视图；

图3(a)是水平层状纵波速度模型；

图3(b)是现有逆时干涉定位结果；

图4是本发明提供的多尺度地面微地震逆时干涉定位方法的流程图；

图5(a)是西南油气田某区块的地面微地震观测系统；

图5(b)是定位用的一维层状介质速度模型；

图6是原始地面微地震波场记录；

图7(a)至图7(f)分别是原始地面微地震波场记录经过多尺度分解之后第1尺度分量波形、第2尺度分量波形、第3尺度分量波形、第4尺度分量波形、第5尺度分量波形、第6尺度分量波形；

图8(a)至图8(f)分别是第1尺度分量波形、第2尺度分量波形、第3尺度分量波形、第4尺度分量波形、第5尺度分量波形、第6尺度分量波形的逆时干涉定位成像俯视图；

图9(a)至图9(f)分别是第1尺度分量波形、第2尺度分量波形、第3尺度分量波形、第4尺度分量波形、第5尺度分量波形、第6尺度分量波形的逆时干涉定位成像侧视图；

图10(a)至图10(c)分别是多尺度逆时干涉定位结果的X-Y切片图、X-Z切片图、Y-Z切片图；

图11(a)至图11(c)分别是常规逆时干涉定位结果的X-Y切片图、X-Z切片图、Y-Z切片图。

具体实施方式

以下结合附图和具体实施例对本发明作具体的介绍。

参照图4，本发明提供的多尺度地面微地震逆时干涉定位方法包括以下步骤：

Step 1：读入地面微地震波场记录、观测系统文件、速度模型文件及相关参数文件，该相关参数文件包括：参与计算地震工区的网格点数、网格间距、采样频率、震源子波参考频率和地震记录时长。

Step 2：利用随机噪声能量随多尺度分解后逐渐发散而有效信号能量却逐渐突出的特点，采用变分模态分解(VMD)(或经验模态分解、S变换、小波变换等时频分析技术)将原始微地震数据分解到多个尺度域，通过求解带约束的变分优化问题，将原始信号分解成为多个尺度的内蕴模式函数(IMF)的集合：

其中，u_k和w_k分别是第k个尺度的IMF和其对应的中心频率，δ(t)是Dirac函数，t是时间采样点，j是复数，e是自然常数，f(t)是输入信号，这里f(t)即为读入的地面微地震波场记录。

为增加求解过程的稳定性，在上式中常加入惩罚因子α和Lagrangian乘法算子λ，因此上式可改写成：

该式可用交替方向乘子算法(ADMM)进行求解。

Step 3：对VMD分解得到的各尺度IMF进行独立成分分析(ICA)以压制背景噪声。

采用基于负熵最大Fast ICA算法对各尺度IMF进行独立成分分析，提取有效信号压制背景噪声，具体的：

首先，将观测信号X中心化，使其均值为0；

其次，初始化需要估计的分量个数m，迭代求解白化矩阵W，使得观测信号X在新的子空间的投影为白化向量Z(t)＝WX(t)，W_p的赋值方式为：

其中，p为迭代求解的次数，W_p ^T为W_p的转置，E为取均值计算，g(y)＝tanh(1,y)，g’为g的导数，W₀为随机的初始权矢量，通过下式对W_p进行更新：

若W_p收敛，则W＝W_p，Y(t)＝WX(t)即为去噪后的各尺度IMF，否则p＝p+1，重新初始化随机向量W₀，继续迭代求解，若p>m，则W＝W₀。

Step 4：将多尺度VMD分解和ICA噪声压制后的各尺度波形分量加载在检波器位置，进行逆时延拓，以波场模的自相关作为成像条件：

其中，I为波场u的自相关逆时干涉定位成像结果。

Step 5：通过峰值信噪比(PSNR)判断某个尺度的逆时干涉定位成像结果是否有效，即有效尺度指该尺度波形分量的逆时延拓波场能量足够聚焦形成清晰的定位椭圆：

其中，N是定位成像I中所有点的个数，D_n是I的一个子集，它由数值小于某个阈值的所有点构成，即：

当某个尺度逆时定位成像的PSNR大于事先给定的阈值时，判定该尺度为有效尺度。

Step 6：拾取有效尺度的逆时干涉定位成像结果的乘积作为最终的定位成像以此加强逆时干涉波场在震源位置的聚焦性，压制稀疏观测系统和背景噪声所引起的成像噪声：

其中，n为有效尺度的个数，Iⁱ为第i个有效尺度的波场u的自相关成像。

为了验证本发明提供的多尺度地面微地震逆时干涉定位方法具有很好的震源定位精度，我们在西南油气田某区块进行了地面微地震逆时干涉定位试验。

西南油气田某区块的地面微地震观测系统如图5(a)所示，其中，三角代表检波器，共9个检波器，观测系统较为稀疏。

定位用的速度模型为一维层状介质速度模型，如图5(b)所示。

我们采用本发明提供的方法进行震源定位，其中：

(1)读入的文件：地面微地震原始波场记录(OriginalWave.seed)、观测系统文件(geometry.dat)、平滑后的层状速度模型文件(layer_vel.dat)。

原始地面微地震波场记录如图6所示。

由图6可知，原始微地震波场的背景噪声干扰严重，难以观测到有效信号，信噪比约为0.5db。

(2)设定的相关参数文件：参与计算地震工区的网格点数、网格间距、采样频率、震源子波参考频率和地震记录时长。

(3)利用变分模态分解(VMD)对原始地震记录进行6层的多尺度分解并进行独立成分分析去噪，得到各尺度分量波形。

原始地面微地震波场记录经过多尺度分解之后的各尺度分量波形如图7(a)至图7(f)所示。

由图7(a)至图7(f)可知，第1尺度分量波形和第2尺度分量波形主要为背景噪声，第3尺度分量波形至第6尺度分量波形包含幅度较强的有效信号。

(4)对各尺度分量波形分别进行逆时干涉定位。

各尺度分量波形的逆时干涉定位成像俯视图如图8(a)至图8(f)所示。

各尺度分量波形的逆时干涉定位成像侧视图如图9(a)至图9(f)所示。

由图8(a)至图8(f)以及图9(a)至图9(f)可知，第1尺度分量波形和第2尺度分量波形的逆时干涉定位结果主要为成像噪声，第3尺度分量波形至第6尺度分量波形的逆时干涉定位结果除浅层的成像噪声之外，包含能量很强的定位聚焦椭圆。

(5)利用峰值信噪比(PSNR)判断有效尺度分量波形为第3尺度分量波形至第6尺度分量波形，将第3尺度分量波形至第6尺度分量波形的逆时定位结果进行连乘，得到最终的多尺度逆时干涉定位成像结果。

图10(a)、图10(b)和图10(c)分别是多尺度逆时干涉定位结果的X-Y切片图、X-Z切片图和Y-Z切片图。

图11(a)、图11(b)和图11(c)分别是常规逆时干涉定位结果的X-Y切片图、X-Z切片图和Y-Z切片图。

对比图10(a)和图11(a)、对比图10(b)和图11(b)、对比图10(c)和图11(c)可知，应用本发明提供的多尺度地面微地震逆时干涉定位方法，对稀疏观测系统和强背景噪声条件下的地面微地震记录，能够得到能量强、聚焦性好的定位成像椭圆，为辨别真实震源的位置奠定了基础，也检验了本发明提供的方法的正确性和有效性。

由此可见，本发明提供的地面微地震定位方法，将逆时干涉定位与多尺度变分模态分解相结合，通过将原始微地震数据分解到多个具有不同振动特征的尺度域，压制了强背景噪声和稀疏观测系统引起的定位成像假象，从而提高了震源的定位精度，这将会有力的推动微地震逆时干涉定位方法在商业中的大规模应用。

需要说明的是，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

Claims (6)

1.一种多尺度地面微地震逆时干涉定位方法，其特征在于，包括以下步骤：

Step1：读入地面微地震波场记录、观测系统文件、速度模型文件及相关参数文件；

Step2：采用时频分析技术将原始微地震数据分解到多个尺度域，通过求解带约束的变分优化问题，将原始信号分解成为多个尺度的内蕴模式函数的集合：

其中，u_k和w_k分别是第k个尺度的内蕴模式函数和其对应的中心频率，δ(t)是Dirac函数，t是时间采样点，j是复数，e是自然常数，f(t)是读入的地面微地震波场记录；

Step3：对分解得到的各尺度内蕴模式函数进行独立成分分析以压制背景噪声；

Step4：将多尺度分解和独立成分分析噪声压制后的各尺度波形分量加载在检波器位置，进行逆时延拓，以波场模的自相关作为成像条件：

其中，I为波场u的自相关逆时干涉定位成像结果；

Step5：通过峰值信噪比判断某个尺度的逆时干涉定位成像结果是否有效，当某个尺度逆时定位成像的峰值信噪比大于事先给定的阈值时，判定该尺度为有效尺度；

Step6：拾取有效尺度的逆时干涉定位成像结果的乘积作为最终的定位成像

其中，n为有效尺度的个数，Iⁱ为第i个有效尺度的波场u的自相关成像。

2.根据权利要求1所述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step1中，所述相关参数文件包括：参与计算地震工区的网格点数、网格间距、采样频率、震源子波参考频率和地震记录时长。

3.根据权利要求1所述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step2中，采用变分模态分解、经验模态分解、S变换或小波变换将原始微地震数据分解到多个尺度域。

4.根据权利要求1所述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step2中，为增加求解过程的稳定性，在公式中加入惩罚因子α和Lagrangian乘法算子λ，将公式改写成：

改写后的公式可用交替方向乘子算法进行求解。

5.根据权利要求1所述的多尺度地面微地震逆时干涉定位方法，其特征在于，在Step3中，采用基于负熵最大Fast ICA算法对各尺度内蕴模式函数进行独立成分分析，提取有效信号压制背景噪声。

6.根据权利要求5所述的多尺度地面微地震逆时干涉定位方法，其特征在于，对各尺度内蕴模式函数进行独立成分分析，提取有效信号压制背景噪声，具体方法为：

首先，将观测信号X中心化，使其均值为0；

其次，初始化需要估计的分量个数m，迭代求解白化矩阵W，使得观测信号X在新的子空间的投影为白化向量Z(t)＝WX(t)，W_p的赋值方式为：

其中，p为迭代求解的次数，W_p ^T为W_p的转置，E为取均值计算，g(y)＝tanh(1,y)，g’为g的导数，W₀为随机的初始权矢量，通过下式对W_p进行更新：

若W_p收敛，则W＝W_p，Y(t)＝WX(t)即为去噪后的各尺度内蕴模式函数，否则p＝p+1，重新初始化随机向量W₀，继续迭代求解，若p>m，则W＝W₀。