CN105549068A - 一种三维各向异性微地震干涉逆时定位方法及系统 - Google Patents

Abstract

本发明涉及一种三维各向异性微地震干涉逆时定位方法及系统，方法包括：步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；步骤4：对得到的所有波场值进行干涉成像，得到成像结果；步骤5：对成像结果进行分析，选取波场值最大的点确定为震源位置。本发明采用波动方程类定位方法进行微地震定位，具有较高的定位精度和可靠性等优点；采用GPU加速三维各向异性微地震逆时定位，大大提高计算速度，提高了计算效率。

Description

一种三维各向异性微地震干涉逆时定位方法及系统

技术领域

本发明涉及一种微地震定位方法，尤其涉及一种三维各向异性微地震干涉逆时定位方法及系统，属于地球物理勘探领域。

背景技术

传统的地震定位工作最早可以追溯到1910年盖戈的工作(Geiger,1912；BrattandBache,1988)。后来进一步完善形成了现在地震定位中最为普遍的方法---盖戈方法。盖戈方法的核心包括两部分，第一，把地震的震源位置与地震波走时的非线性关系线性化；第二，利用最小二乘的方法求解此线性系统；由于其完整的求解及其评价体系，这类方法一直以来被广泛采用。然而这种方法也有其局限性的，它需要从地震记录中获得清晰的纵横波走时信息(Flinn,1965)，要求地震记录具有较高的信噪比.但在实际微地震监控资料中，低信噪比是观测数据的普遍特征，且观测数据量巨大，因此传统的盖戈类方法受到了很大的限制。

同时近年来，人们借鉴地震勘探中的偏移成像原理也发展了不需要拾取震相走时信息、适用于低信噪比数据的类偏移定位方法.这类方法把微地震源类似成偏移成像中的绕射点，利用反射地震学中处理绕射点的成像方法进行震源的定位。相应的定位过程可以分为两步，即首先延拓观测数据“重构”地下波场，之后施加合适的“成像条件”得到震源位置和激发时间。干涉逆时成像定位技术有着适应低信噪比数据、不需要拾取震相走时信息,较高的定位精度和可靠性等优点。

近十年中，利用图形处理单元(GPU)进行计算密集型应用的加速实现已经得到了突飞猛进式的发展。图形处理单元(GPU)由于其具有高速的内存带宽,相较于CPU至少高出两个数量级的计算处理核心，更适合并行计算的单指令多数据(SIMD)计算模式，以及更低的能耗成本，正广泛的应用到计算科学的相关领域。对于勘探地球物理领域，对于使用图形处理单元(GPU)的兴趣也在显著增强，越来越多的研究已经将GPU用于加速地震处理中的核心算法，例如地震数值模拟、地震成像、地震高精度反演等。

发明内容

本发明所要解决的技术问题是提供一种基于GPU加速的三维各向异性微地震干涉逆时定位方法及系统，将GPU加速技术与三维各向异性微地震干涉逆时定位算法结合得到一种高效、快速的微地震定位方法，进而为大规模微地震监测提供更有力的保证。

本发明解决上述技术问题的技术方案如下：一种三维各向异性微地震干涉逆时定位方法，包括以下步骤：

步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；

步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；

步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

步骤4：对得到的所有波场值进行干涉成像，得到成像结果；

步骤5：对成像结果进行分析，选取波场值最大的点确定为震源位置。

本发明的有益效果是：采用波动方程类定位方法进行微地震定位，不需要人工提取初值，相较于传统的走时类方法具有适应低信噪比数据、不需要拾取震相走时信息，具有较高的定位精度和可靠性等优点；采用GPU加速三维各向异性微地震逆时定位，大大提高计算速度，相较于常规CPU实现，可以达到20-30倍的速度提升，提高了计算效率，使得逆时成像定位技术真正可以达到现场实时监控的生产性要求。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述步骤1具体包括以下步骤：

步骤1.1：根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

步骤1.2：采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

采用上述进一步方案的有益效果是，将介质模型离散成网格点，网格点作为波场值的载体，由三维结构的多个网格点构成波场，作为波场传播的初始条件。

进一步，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震源函数的具体公式为：

s(x,y,z,t)＝g(x,y,z)·f(t)公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)公式(2)

$g(x,y,z)=\frac{\exp \{\[-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\]\}}{{\mathrm{\β}}^2}$ 公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀,y₀,z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

进一步，所述步骤3具体包括以下步骤：

步骤3.1：将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长必须满足该数值格式的稳定性条件；

步骤3.2：采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。

进一步，所述步骤4还包括在成像空间域对成像结果数据做干涉处理，得到精确的成像结果。

进一步，所述步骤4中采用类WDF(WignerDistributionFunction)数据变化进行干涉处理。

本发明解决上述技术问题的技术方案如下：一种三维各向异性微地震干涉逆时定位系统，包括建模离散模块、震源函数模块、传播方程模块、干涉成像模块和震源确定模块；

所述建模离散模块用于建立介质模型，对介质模型进行网格离散得到多个网格点；

所述震源函数模块用于计算震源函数，根据震源函数计算每个网格点上的压力值；

所述传播方程模块用于将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

所述干涉成像模块用于对得到的所有波场值进行干涉成像，得到成像结果；

所述震源确定模块用于对成像结果进行分析，选取波场值最大的点确定为震源位置。

本发明的有益效果是：采用波动方程类定位方法进行微地震定位，不需要人工提取初值，相较于传统的走时类方法具有适应低信噪比数据、不需要拾取震相走时信息，具有较高的定位精度和可靠性等优点；采用GPU加速三维各向异性微地震逆时定位，大大提高计算速度，相较于常规CPU实现，可以达到20-30倍的速度提升，提高了计算效率，使得逆时成像定位技术真正可以达到现场实时监控的生产性要求。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述建模离散模块包括建模模块和离散模块；

所述建模模块用于根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

所述离散模块用于采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

采用上述进一步方案的有益效果是，将介质模型离散成网格点，网格点作为波场值的载体，由三维结构的多个网格点构成波场，作为波场传播的初始条件。

进一步，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震源函数的具体公式为：

s(x,y,z,t)＝g(x,y,z)·f(t)公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)公式(2)

$g(x,y,z)=\frac{\exp \{\[-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\]\}}{{\mathrm{\β}}^2}$ 公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀,y₀,z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

进一步，所述传播方程模块包括替代模块和计算模块；

所述替代模块用于将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长必须满足该数值格式的稳定性条件；

所述计算模块采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。

附图说明

图1为本发明实施例1所述的一种三维各向异性微地震干涉逆时定位方法流程图；

图2为本发明实施例1所述的一种三维各向异性微地震干涉逆时定位系统结构框图；

图3为本发明所述方法在GPU设备上的实施流程图；

图4为本发明三维水平层状模型P波速度模型及震源位置图；

图5a为利用三分量数据进行逆时成像定位得到的水平层状模型微地震事件U_Z分量的定位结果；

图5b为利用三分量数据进行逆时成像定位得到的水平层状模型微地震事件U_X分量的定位结果；

图5c为利用三分量数据进行逆时成像定位得到的水平层状模型微地震事件U_Y分量的定位结果。

附图中，各标号所代表的部件列表如下：

1、建模离散模块，2、震源函数模块，3、传播方程模块，4、干涉成像模块，5、震源确定模块。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

如图1所示，为本发明实施例1所述的一种三维各向异性微地震干涉逆时定位方法，包括以下步骤：

步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；

步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；

步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

步骤4：对得到的所有波场值进行干涉成像，得到成像结果；

步骤5：对成像结果进行分析，选取波场值最大的点确定为震源位置。

实施例2中，在实施例1的基础上，所述步骤1具体包括以下步骤：

步骤1.1：根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

步骤1.2：采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

实施例3中，在实施例1或2的基础上，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震源函数的具体公式为：

s(x,y,z,t)＝g(x,y,z)·f(t)公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)公式(2)

$g(x,y,z)=\frac{\exp \{\[-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\]\}}{{\mathrm{\β}}^2}$ 公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀,y₀,z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

实施例4中，在实施例1-3任一实施例的基础上，所述步骤3具体包括以下步骤：

步骤3.1：将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长必须满足该数值格式的稳定性条件；

步骤3.2：采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。

此实施例中，进一步，所述步骤3.2中采用区域分解的方式利用GPU在各个计算节点对三维各向异性弹性波方程进有限差分计算，得到每一刻的波场值的具体方法如下：

首先，给出三维各向异性介质下的弹性波波动方程的数学描述，三维各向异性弹性波动方程其形式为：

$\mathrm{\ρ}{\∂}_{tt}^2{u}_i={\∂}_j{\mathrm{\σ}}_{ij}+F_i,$ 公式(4)

式中，u为波场函数；σ为柯西应力张量；F_i是单位体积的体力项，其可以被等效为应力源，ρ是介质物质密度，表示二阶时间偏导。对于应力张量利用线性弹性理论可以得到：

σ_ij＝c_ijklε_kl，公式(5)

式中ε_kl为线性应变张量，c_ijkl为四阶张量的弹性刚度矩阵，该参数需要根据介质的实际情况给出；而对于ε_kl我们有

$\begin{array}{cc}{\mathrm{\ϵ}}_{kl}=\frac{1}{2}\mathrm{\[}{\∂}_k{u}_l+{\∂}_l{u}_k\mathrm{\]}& k,l=1,2,3\end{array}$ 公式(6)

其中ε_kl代表线性应变张量，代表对kth维求空间导数，u_l表示lth波场位移。根据此式我们就可以得行弹性波动方程的数值模拟，利用GPU进行的弹性波数值模拟的具体计算步骤如下：

3.2.1我们定义一个N_x×N_y×N_z的计算网格，并且指定N_t个时间步长，这样在空间和时间的每一个网格点就可以利用一个四元组进行表示：即，[x,y,z|t]＝[pΔx,qΔy,rΔz|nΔt]。因此在指定点的连续的波场位移记录就可以利用离散的网格点进行表示：

u_i|_x,y,z|t≈u_i ^p,r,q|n公式(7)

我们近似一阶偏导通过一个如下的M阶精度的的中心差商算子D_x[·]，可以得到：

${\∂}_x{u}_j\≈D_x\[u_j^{p,q,r|n}\]=\frac{1}{\mathrm{\Δ}x}\underset{\mathrm{\α}=1}{\overset{M}{\Σ}}{W}_{\mathrm{\α}}(u_j^{p+\mathrm{\α},q,r|n}-u_j^{p-\mathrm{\α},q,r|n})$ 公式(8)

式中W_α是一个权重系数多项式，即有限差分系数，M为有限差分计算阶数；

空间差分算子D_y[·]和D_z[·]与D_x[·]相似分别表示y和z方向的偏导数。将所有这三者差分算子代入方程(6),用以求解本构方程中的偏导数。

3.2.2对于时间离散，我们使用一个标准的空间二阶精度近似来计算方程(6)中的

${\∂}_{tt}^2{u}_j\≈D_{tt}{u}_j^{p,q,r|n}=\frac{1}{{\mathrm{\Δt}}^2}\[u_j^{p,q,r|n+1}-2u_j^{p,q,r|n}+u_j^{p,q,r|n-1}\]$ 公式(9)

将上述的差分算子代入方程(4)并且通过重新排列这些离散项，我们就可以通过时间步增的方式进行未知波场求解。例如对于我们需要给定其当前和之前的波场值以及根据有限差分计算空间计算模版所需要x-,y-,z-方向相邻点。图3表述了空间有限差分计算模版的在当前时间步长关于点所需要的计算网格点。

实施例5中，在实施例1-4任一实施例的基础上，所述步骤4还包括在成像空间域对成像结果数据做干涉处理，得到精确的成像结果。通过在成像空间域对数据做干涉处理，得到精确的定位结果。

实施例6中，在实施例5的基础上，所述步骤4中采用类WDF(WignerDistributionFunction)数据变化进行干涉处理，公式如下：

$I(x,t)={\&Integral;}_{\left|t_h\right|<T}{\&Integral;}_{\left|x_h\right|<X}W(x-\frac{x_h}{2},t-\frac{t_h}{2})\&times;W(x+\frac{x_h}{2},t+\frac{t_h}{2}){\mathrm{dt}}_h{\mathrm{dx}}_h$ 公式(10)

式中x表示成像空间中的每个成像点坐标，W(x,t)表示成像空间内x点在时刻t的波场值，I(x,t)表示干涉处理后的空间域结果，X和T表示选取的类Wigner分布窗口大小。

所述步骤5中随着时间逆推进行每一步的干涉结果累加，最终得到成像结果其中t_a为记录终止时间，t_b为记录开始时间。

如图2所示，为本发明实施例1所述的一种三维各向异性微地震干涉逆时定位系统，包括建模离散模块1、震源函数模块2、传播方程模块2、干涉成像模块4和震源确定模块5；

所述建模离散模块1用于建立介质模型，对介质模型进行网格离散得到多个网格点；

所述震源函数模块2用于计算震源函数，根据震源函数计算每个网格点上的压力值；

所述传播方程模块3用于将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

所述干涉成像模块4用于对得到的所有波场值进行干涉成像，得到成像结果；

所述震源确定模块5用于对成像结果进行分析，选取波场值最大的点确定为震源位置。

本发明的基础是微地震逆时定位以及其高性能GPU计算方法，本发明基于GPU计算加速技术实现了基于GPU的三维弹性波各向异性微地震逆时定位算法，在具体示例中实施步骤分别为：

1.介质模型建立：

我们选择三维水平层状模型作为测试模型，其弹性参数如表1，其中横波速度选取为并采用均匀介质密度取为2000kg/m³。这些参数都通过相应的的变换公式(Thomsen,1986)转换进弹性刚度矩阵中，从而作为模型参数应用的算法中。将震源放置在(x,y,z)＝(760,760,760)m处，即模型的中心处，如图4所示。

	Vp(m/s)	Vs(m/s)	ε	δ	γ
						第一层	2000	1155	0.2	-0.1	0.2
第二层	2500	1443	0.3	-0.15	0.25
						第三层	3000	1733	0.4	-0.2	0.3

表1：层状介质弹性参数表

2.模型离散化：

对我们设计的介质模型进行网格离散，测试地震数据维度大小为N_x×N_y×N_z＝380³，网格间距为Δx＝Δy＝Δz＝0.004km。

3.震源函数：

震源函数在空间上采用高斯函数，时间上采用Ricker子波，形式如下：

s(x,y,z,t)＝g(x,y,z)·f(t)

式中：

f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)

$g(x,y,z)=\frac{\exp \{\&lsqb;-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\&rsqb;\}}{{\mathrm{\&beta;}}^2}$

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀,y₀,z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置；

4.弹性波动方程的数值格式：

我们采用8阶精度进行有限差分数值模拟，则对应的我们利用公式(4)-(6)进行每一步的波场计算。

5.微地震事件成像条件：

我们对每一步获得的三分量波场，施加如下的干涉成像条件：

$I(x,t)={\&Integral;}_{\left|t_h\right|<T}{\&Integral;}_{\left|x_h\right|<X}W(x-\frac{x_h}{2},t-\frac{t_h}{2})\&times;W(x+\frac{x_h}{2},t+\frac{t_h}{2}){\mathrm{dt}}_h{\mathrm{dx}}_h$

并随着时间逆推将每一步的成像结果进行累加，最后得到微地震事件定位结果。

6.数值模拟结果分析：

三分量的定位结果如图5(a)-(c)所示。这些图表明通过逆时成像方法获得的震源定位结果与真实的震源位置相吻合(图4)，验证了我们方法的有效性。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种三维各向异性微地震干涉逆时定位方法，其特征在于，包括以下步骤：

步骤1：建立介质模型，对介质模型进行网格离散得到多个网格点；

步骤2：计算震源函数，根据震源函数计算每个网格点上的压力值；

步骤3：将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

步骤4：对得到的所有波场值进行干涉成像，得到成像结果；

步骤5：对成像结果进行分析，选取波场值最大的点确定为震源位置。

2.根据权利要求1所述的一种三维各向异性微地震干涉逆时定位方法，其特征在于，所述步骤1具体包括以下步骤：

步骤1.1：根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

步骤1.2：采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

3.根据权利要求1所述的一种三维各向异性微地震干涉逆时定位方法，其特征在于，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震源函数的具体公式为：

s(x,y,z,t)＝g(x,y,z)·f(t)公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)公式(2)

$g(x,y,z)=\frac{\exp \{\&lsqb;-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\&rsqb;\}}{{\mathrm{\&beta;}}^2}$ 公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀,y₀,z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

4.根据权利要求1-3任一项所述的一种三维各向异性微地震干涉逆时定位方法，其特征在于，所述步骤3具体包括以下步骤：

步骤3.1：将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长必须满足该数值格式的稳定性条件；

步骤3.2：采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。

5.根据权利要求1所述的一种三维各向异性微地震干涉逆时定位方法，其特征在于，所述步骤4还包括在成像空间域对成像结果数据做干涉处理，得到精确的成像结果。

6.根据权利要求5所述的一种三维各向异性微地震干涉逆时定位方法，其特征在于，所述步骤4中采用类WDF数据变化进行干涉处理。

7.一种三维各向异性微地震干涉逆时定位系统，其特征在于，包括建模离散模块、震源函数模块、传播方程模块、干涉成像模块和震源确定模块；

所述建模离散模块用于建立介质模型，对介质模型进行网格离散得到多个网格点；

所述震源函数模块用于计算震源函数，根据震源函数计算每个网格点上的压力值；

所述传播方程模块用于将三维各向异性弹性波方程转换为传播方程，将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值；

所述干涉成像模块用于对得到的所有波场值进行干涉成像，得到成像结果；

所述震源确定模块用于对成像结果进行分析，选取波场值最大的点确定为震源位置。

8.根据权利要求7所述的一种三维各向异性微地震干涉逆时定位系统，其特征在于，所述建模离散模块包括建模模块和离散模块；

所述建模模块用于根据地质背景条件、实际测得的岩石物理测试数据和测井资料数据建立介质模型；

所述离散模块用于采用形状规则的三维网格对介质模型进行网格离散，得到网格点。

9.根据权利要求7所述的一种三维各向异性微地震干涉逆时定位系统，其特征在于，所述震源函数在空间上采用高斯函数，在时间上采用Ricker子波，所述震源函数的具体公式为：

s(x,y,z,t)＝g(x,y,z)·f(t)公式(1)

其中，f(t)＝(1-2(πf₀t)²)exp(-(πf₀t)²)公式(2)

$g(x,y,z)=\frac{\exp \{\&lsqb;-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\&rsqb;\}}{{\mathrm{\&beta;}}^2}$ 公式(3)

式中：t表示时间，f₀表示Ricker子波的中心频率，模型计算中f₀＝15Hz，β为常数；(x₀,y₀,z₀)为震源的空间位置，x、y和z分别为x轴、y轴和z轴方向上的位置。

10.根据权利要求7-9任一项所述的一种三维各向异性微地震干涉逆时定位系统，其特征在于，所述传播方程模块包括替代模块和计算模块；

所述替代模块用于将三维各向异性弹性波方程中的微分用差分近似替代，得到相应的有限差分格式的传播方程，所述传播方程中空间采样步长和时间采样步长必须满足该数值格式的稳定性条件；

所述计算模块采用区域分解的方式将每个网格点上的压力值带入传播方程进行计算，得到每一刻的波场值。