CN112051609B - 地震波成像方法、系统、存储介质、计算机程序、终端 - Google Patents

Abstract

本发明属于多维成像技术领域，公开了一种地震波成像方法、系统、存储介质、计算机程序、终端，将一维方程拓展到二维，将初始模型通过傅氏变换到频率‑波数域；在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；再按地质模型中多组参考参数与实际参数之间的关系，用适当插值方程求得最终的波场延拓结果；执行频率j＝1…j＝J和深度i＝1…i＝I的运算，完成整个二维剖面的波场延拓运算。本发明适用于渗流介质的有效的多物理场、多参数模型的地震数值模拟和偏移成像新方法；提高了孔、渗介质地震波场特征研究的学术水平和油气储层预测效果，提供理论基础与实用技术的目的。

Description

地震波成像方法、系统、存储介质、计算机程序、终端

技术领域

本发明属于多维成像技术领域，尤其涉及一种地震波成像方法、系统、存储介质、计算机程序、终端。

背景技术

目前，因为地震波的传播过程涉及的物理场来说比较多的，除了岩石自身的矿物组分、地质构造和岩石微构造外，还包括孔隙和裂缝中含的流体成份和性质，甚至目标地质体所处的环境，如压力场，位移场，排水情况下的达西速度场，以及非稳定流动状态下的达西加速度场等多会严重影响地震波的传播，而描述孔、渗介质中地震波传播特征的传播方程也是复杂、多样的。一般分为两大类，一类是含有单个物理量(位移或压力等)和多个参数(波速、黏滞系数、弥散系数、弛豫时间等)的波动方程；另一类是含有多个物理场(位移、达西速度和压力等)和多个参数(速度、黏滞系数、密度和渗透率等)的微分方程。

当前的单物理场多参数的波动方程，可以进行多维成像，但局限于其应用条件，不能完全地描述解释地震波传播过程中的多种物理场，及其受到的多种参数影响的波场响应情况；而多物理场多参数的微分方程，由于方程所含物理场多，参数多，导致其求解方程困难，计算复杂，所以多维的波场成像十分困难，只能对其进行属性研究。

通过上述分析，现有技术存在的问题及缺陷为：含有单个物理量和多个参数的波动方程，不能很好的描述解释地震波传播的多种物理场和受到的多种参数影响的波场情况；多物理场多参数的微分方程，由于方程所含物理场多，参数多，导致其求解方程困难，计算复杂，难以多维的波场成像十分困难，只能对其进行属性研究。

解决以上问题及缺陷的难度为：地震波在含流体孔、渗介质中的传播将引起弹性波和流体流动之间的相互作用，使地震波能量衰减，波形和频谱畸变。地震波在含流体孔、渗介质中的传播是一个十分复杂的过程，涉及的参数和物理场很多，既包括岩石固体骨架和胶结物的性质，又包括流体成分与含量，还包括流体流动与岩石骨架和地震波传播的互相作用。因此描述孔、渗介质中地震波传播特征的传播方程也是复杂、多样的。由于方程所含物理场多，参数多，导致其求解方程困难，计算复杂，多维的波场成像十分困难，目前只能对其进行属性研究。

解决以上问题及缺陷的意义为：通过设立等效参数和等效物理场的方法，把具有多参数和多物理场的传播方程化为等效单参数和单物理场的传播方程，通过改进的波场延拓方法，实现非均匀介质的地震波场延拓，发展出适用于渗流介质的有效的多物理场、多参数模型的地震数值模拟和偏移成像新方法，通过地震波场正演和反演成像来模拟和解释地震波的这些变化特征，有助于人们预测地下油气储层所含流体成分，识别地层岩性，正确部署勘探和开发方案。

发明内容

针对现有技术存在的问题，本发明提供了一种地震波成像方法、系统、存储介质、计算机程序、终端。

本发明是这样实现的，一种地震波成像方法，所述地震波成像方法包括：

第一步，将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域；

第二步，在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；

第三步，再按地质模型中多组参考参数与实际参数之间的关系，用适当插值方程求得最终的波场延拓结果；

第四步，执行频率j＝1…j＝J和深度i＝1…i＝I的运算，完成整个二维剖面的波场延拓运算。

进一步，所述第一步将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域：

式中

ε＝iλω＝ρ_fk_k/η；

孔隙度φ，孔隙度形态因子α_G，渗透率k_k，粘滞系数η，流体密度ρ_f，流体压缩系数β_f，介质速度v_b，介质与颗粒弹性模量K，无量纲的参数K_sg，K_fg。

进一步，所述第二步预先取若干组参数，在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；在ΔZ＝Z_i+1-Z_i的间隔内的第一波场延拓公式为：

式中，K_z＝k±iα，k和α均为实数，α又称为吸收系数；

对于快纵波有：

舍去ε的二次以上高阶项，得：

进一步，所述第二步还包括：预先取若干组参数

l＝1,2…L；若L＝3，则有3组参数C¹，C²，C³，对S_i,j(Z_i,x_j)分别按三组在x方向不变的参数，延拓到Z_i+1，于是在深度Z_i+1，得到

这三个数据。

进一步，所述第三步还包括：ΔZ＝Z_i+1-Z_i的间隔内，采用Gazdag相移加插值PSPI的思路解决参数横向变化问题。按地质模型中该道的具体参数C同C¹，C²，C³的关系，用插值方程求得与C对应的地震振幅S_i+1(Z_i+1,x_j)，下标j为地震道标号。

进一步，所述第四步还包括：采用拉格朗日插值方法进行计算：

本发明的另一目的在于提供一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行权利要求任意一项所述包括下列步骤：

第一步，将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域；

第二步，在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；

第三步，再按地质模型中多组参考参数与实际参数之间的关系，用适当插值方程求得最终的波场延拓结果；

第四步，执行频率j＝1…j＝J和深度i＝1…i＝I的运算，完成整个二维剖面的波场延拓运算。

本发明的另一目的在于提供一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施所述的地震波成像方法。

本发明的另一目的在于提供一种实施所述的地震波成像方法的地震波成像系统，所述地震波成像系统包括：

参数预处理模块，用于对波场延拓时的多组参考参数进行预处理；

波场延拓模块，用于对输入的多组参考参数和初始模型进行计算，得到最终的波场延拓结果。

本发明的另一目的在于提供一种终端，所述终端搭载所述的地震波成像系统。

结合上述的所有技术方案，本发明所具备的优点及积极效果为：本发明通过设立等效参数和等效物理场的方法，研发适用于渗流介质的有效的多物理场、多参数模型的地震数值模拟和偏移成像新方法；提高了孔、渗介质地震波场特征研究的学术水平和油气储层预测效果，提供理论基础与实用技术的目的。

地震波在含流体孔、渗介质中的传播将引起弹性波和流体流动之间的相互作用，使地震波能量衰减，波形和频谱畸变。地震波在含流体孔、渗介质中的传播是一个十分复杂的过程，涉及的参数和物理场很多，既包括岩石固体骨架和胶结物的性质，又包括流体成分与含量，还包括流体流动与岩石骨架和地震波传播的互相作用。因此描述孔、渗介质中地震波传播特征的传播方程也是复杂、多样的。通过设立等效参数和等效物理场的方法，把具有多参数和多物理场的传播方程化为等效单参数和单物理场的传播方程，通过改进的波场延拓方法，实现非均匀介质的地震波场延拓，发展出适用于渗流介质的有效的多物理场、多参数模型的地震数值模拟和偏移成像新方法，通过地震波场正演和反演成像来模拟和解释地震波的这些变化特征，有助于人们预测地下油气储层所含流体成分，识别地层岩性，正确部署勘探和开发方案。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的地震波成像方法流程图。

图2是本发明实施例提供的地震波成像系统的结构示意图；

图中：1、参数预处理模块；2、波场延拓模块。

图3是本发明实施例提供的地震波成像方法实现流程图。

图4是本发明实施例提供的衰减参数与频率振幅相关关系的曲线示意图；

图中：(a)渗透率为定值，黏滞系数改变的频率振幅响应；(b)黏滞系数为定值，渗透率改变的频率振幅响应。

图5是本发明实施例提供的二维模型横向多参数变化的正演与偏移成像示意图；

图中：(a)地质模型；(b)正演结果；(c)偏移成像结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对现有技术存在的问题，本发明提供了一种地震波成像方法、系统、存储介质、计算机程序、终端，下面结合附图对本发明作详细的描述。

如图1所示，本发明提供的地震波成像方法包括以下步骤：

S101：将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域；

S102：在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；再按地质模型中多组参考参数与实际参数之间的关系，用适当插值方程求得最终的波场延拓结果；

S103：执行频率j＝1…j＝J和深度i＝1…i＝I的运算，完成整个二维剖面的波场延拓运算。

如图2所示，本发明提供的地震波成像系统包括：

参数预处理模块1，用于对波场延拓时的多组参考参数进行预处理；

波场延拓模块2，用于对输入的多组参考参数和初始模型进行计算，得到最终的波场延拓结果。

下面结合附图对本发明的技术方案作进一步的描述。

本发明首选Silin方程(Silin and Goloshubin,2010)为地震多物理场多参数传播方程的代表进行波场延拓。其优点在于Silin方程是针对地震波在含流体孔渗介质中的传播推导出来的，推导过程中考虑了达西定律、虎克定律。动量和质量的平衡。且与广泛应用的Biot方程有直接联系。该方程含有Biot快、慢纵波和横波，并且独立传播，有利于快慢纵波和偏振方向互相垂直的两个横波的特征研究。下面的公式(1)中分别为纵波和两个横波的方程。式中的U₀，Q₀和P₀分别为岩石骨架位移，达西速度场和流体压力。参数ρ_b，ρ_f为岩石体密度和流体密度，M为体积模量，μ为剪切模量，k_x为水平波数，k，η和τ分别为渗透率，粘滞系数和弛豫时间，ω为频率，其他为中间参数：

Silin采用了小参数ε的渐近解。他定义

又称

为流体流动性。最后得到快慢纵波和横波的波数和吸收系数。下式上标F代表快纵波，s代表慢纵波：

γ_β，γ_K为复合参数：

而

均为γ_β，γ_K的函数。v_b，v_f为地震波在岩石和流体的传播速度。

吸收系数为：

上述公式均是一维的，因此不用考虑在多维模型中参数的横向变化。但是，难以在实际地震勘探中得到应用。为此必须解决多维地震数据的垂直延拓问题。

本发明对纵波二维地震资料S_i,j(Z_i+1,x_j)来说计算方法和步骤：

第一步将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域：

式中

ε＝iλω＝ρ_fk_k/η；

孔隙度φ，孔隙度形态因子α_G，渗透率k_k，粘滞系数η，流体密度ρ_f，流体压缩系数β_f，介质速度v_b，介质与颗粒弹性模量K，无量纲的参数K_sg，K_fg。

第二步，在ΔZ＝Z_i+1-Z_i的间隔内的波场延拓公式为：

式中，K_z＝k±iα，k和α均为实数。α又称为吸收系数。

对于快纵波有：

舍去ε的二次以上高阶项，得：

将(11)代入(9)式；ΔZ＝Z_i+1-Z_i的间隔内，只考虑波场在z方向的变化。这时可采用Gazdag(GazdagandSguazzero,1984)相移加插值PSPI的思路来解决参数横向变化问题。其具体计算步骤如下：预先取若干组参数

，l＝1,2…L；若L＝3，则有3组参数C¹，C²，C³，对S_i,j(Z_i,x_j)分别按三组在x方向不变的参数，延拓到深度Z_i+1，于是对z方向上的每个地震道j,都可得到

这三个数据。再按地质模型中该道的具体参数C同C¹，C²，C³的关系，用插值方程求得与C对应的地震振幅S_i+1(Z_i+1,x_j)。下标j为地震道标号。

第三步，与第二步方法类似，执行j＝1…j＝J(频率)和i＝1…i＝I(深度方向)的运算，即可完成整个二维剖面的插值运算。可采用拉格朗日插值方法进行计算：

这种插值方法称其为扩展相移加插值法EPSPI。与PSPI有两点主要差别：一是参数多，更复杂。二是许多参数物理意义不同，但相互影响，需要按照具体地质模型参数对地震振幅的影响进行计算分析，并绘制响应曲线，如图4所示。

表1含气砂岩储层模型

表2模型岩石物性表

应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、CD或DVD-ROM的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。

Claims (5)

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

第一步，将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域；

所述第一步将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域：

式中

孔隙度φ，孔隙度形态因子α_G，粘滞系数η，流体密度ρ_f，流体压缩系数β_f，介质与颗粒弹性模量K，无量纲的参数K_sg，K_fg；ω为频率，ρ_b为岩石体密度，U₀为岩石骨架位移，M为体积模量，k_x为水平波数，α为吸收系数，P₀为流体压力，τ为弛豫时间，k为渗透率，

为复合函数γ_β，γ_K的函数，其他为中间参数；第二步，在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；

第三步，再按地质模型中多组参考参数与实际参数之间的关系，用拉格朗日插值方程求得最终的波场延拓结果；所述求得最终的波场延拓结果包括：

ΔZ＝Z_i+1-Z_i的间隔内，采用Gazdag相移加插值PSPI的思路解决参数横向变化问题；按地质模型中每道地震道的具体参数C同C¹，C²，C³的关系，用插值方程求得与C对应的地震振幅S_i+1(Z_i+1,x_j)，下标j为地震道标号；

第四步，执行频率j＝1…j＝J和深度i＝1…i＝I的运算，完成整个二维剖面的波场延拓运算；

所述第二步还包括：预先取若干组参数

v_b为地震波在岩石的传播速度；若L＝3，则有3组参数C¹，C²，C³，对S_i,j(Z_i,x_j)分别按三组在x方向不变的参数，延拓到Z_i+1，于是在深度Z_i+1，得到

这三个数据；

所述第四步还包括：采用拉格朗日插值方法进行计算：

2.如权利要求1所述的地震波成像方法，其特征在于，所述第二步预先取若干组参数，在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；在ΔZ＝Z_i+1-Z_i的间隔内的第一波场延拓公式为：

式中，k_z＝k±iα，k和α均为实数，α又称为吸收系数；

对于快纵波有：

舍去ε的二次以上高阶项，得：

其中，v_b为地震波在岩石的传播速度，η为粘滞系数，

3.一种接收用户输入程序存储介质，所存储的计算机程序使电子设备执行权利要求1-2任意一项所述地震波成像方法包括下列步骤：

第一步，将一维方程拓展到二维，将初始模型通过傅氏变换到频率-波数域；

第二步，在每一个深度间隔内使用多组参考参数进行波场延拓得到多组波场延拓结果；

第三步，再按地质模型中多组参考参数与实际参数之间的关系，用拉格朗日插值方程求得最终的波场延拓结果；

第四步，执行频率j＝1…j＝J和深度i＝1…i＝I的运算，完成整个二维剖面的波场延拓运算。

4.一种实施权利要求1～2任意一项所述的地震波成像方法的地震波成像系统，其特征在于，所述地震波成像系统包括：

参数预处理模块，用于对波场延拓时的多组参考参数进行预处理；

波场延拓模块，用于对输入的多组参考参数和初始模型进行计算，得到最终的波场延拓结果。

5.一种终端，其特征在于，所述终端搭载权利要求4所述的地震波成像系统。

Images