CN112327358A - 一种粘滞性介质中声波地震数据正演模拟方法 - Google Patents

Abstract

本说明书实施例公开了一种粘滞性介质中声波地震数据正演模拟方法。本发明通过将振幅衰减项与相位频散项分离，相比传统的衰减项与频散项耦合在一起的粘声波方程正演模拟，可方便地实现粘声逆时偏移成像。相比传统计算方法，可大大提高计算效率、减少计算内存。

Description

一种粘滞性介质中声波地震数据正演模拟方法

技术领域

本说明书涉及勘探地球物理学领域，尤其涉及一种粘滞性介质中声波地震数据正演模拟方法。

背景技术

由于实际地下介质存在内摩擦力及含流体孔隙，地下介质通常表现为粘滞性。地震波在粘滞性介质中传播时，尤其是强衰减区域，会产生能量衰减和频散现象。在对含衰减介质进行成像时，若不考虑这些黏性衰减及频散效应，则会导致含高衰减特性储层的内部和下部构造成像能量较弱，且成像位置不准。为了实现复杂粘滞性介质中地层的精确成像，考虑黏性振幅衰减及相位频散效应的地震正演模拟极其重要。同时，传统时间域正演模拟方案中，衰减项与频散项耦合在一起不利于逆时偏移成像的实施；频率域正演模拟方法存在计算效率慢、内存消耗大等问题。

因此，需要一种便于成像且高效的正演模拟方法。

发明内容

本发明的目的在于，提供一种准确高效的滞性介质中声波地震数据正演模拟方法。

为解决上述技术问题，本发明采用如下技术方案：

获取初始的参数场，所述参数场包括介质纵波速度场v₀和品质因子Q；

将所述参数场进行矩形的网格划分，确定每个网格的位置，根据所述参数场和所述网格的位置，根据如下方程确定每个网格的模拟系数g和h：

其中，x表示横坐标，z表示纵坐标，ω表示角频率，v₀代表介质纵波速度，l表示虚数单位，k表示空间波数，v表示粘滞性介质的复数值波场传播速度；

确定炮点和检波点位置，根据所述g和h，采用如下方程确定检波点的应力场值：

其中,i表示横向第i个网格，j表示纵向第j个网格，u_p表示应力场值，u_q表示辅助应力场值，k表示离散的时间值，Δx表示横向网格间距，Δz表示纵向网格间距，Δt表示时间采样间隔，a₀和a_nn为预设的差分系数，FFT表示快速傅里叶变换，FFT^-1表示快速傅里叶变换的反变换，k_x和k_z分别表示纵横向的波数值，N表示差分阶数，m和n是计算时的中间变量，无实际物理意义，N表示差分阶数，nn表示差分步长。

本说明书实施例采用的上述至少一个技术方案能够达到以下有益效果：

与现有技术相比，本发明通过将振幅衰减项与相位频散项分离，相比传统的衰减项与频散项耦合在一起的粘声波方程正演模拟，可方便地实现粘声逆时偏移成像。相比传统计算方法，可大大提高计算效率、减少计算内存，实际应用也证明了粘声波模拟很好地应用于地震勘探，为逆时偏移及全波形反演提供一种较好的正演模拟方案。

附图说明

图1为本说明书实施例所提供的一种粘滞性介质中声波地震数据正演模拟方法的流程示意图；

图2为本说明书实施例提供的粘滞性介质中衰减与频散波场模拟示意图；

图3为本说明书实施例所提供的一种模型的粘滞性参数与位置的关系示意图；

图4为本说明书实施例所提供的方法与传统方法所得到的地震炮记录对比的示意图；

图5为本说明书实施例所提供的方法与传统方法所得到的地震炮记录单道对比的示意图。

具体实施方式

为使本申请的目的、技术方案和优点更加清楚，下面将结合本申请具体实施例及相应的附图对本申请技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本申请一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

首先需要对于本申请模拟过程中的计算原理进行说明。

在地震领域中，粘滞性介质的复数值波场传播速度为：

又已知空间频率域粘声波方程表达式为：

其中，x和z分别表示横纵向坐标，ω表示角频率，ω₀表示参考角频率，v₀代表介质纵波速度，γ＝arctan(1/Q)/π是一个无量纲的量，Q为品质因子，l表示虚数单位，U表示频率域波场，s为频率域震源项。

在频率域，利用公式(1)和公式(2)很容易求解粘滞性介质中声波波场传播，但是频率域求解存在计算效率低且内存占用大等缺点。公式(2)又因为其高度复杂性，很难直接转换到时间域。因此，我们对公式(2)中ω²/v²(x,ω)项进行近似转化，得到一种可以方便转换到时间域的近似形式。其近似表达如下所示：

其中，

被用来描述相位频散，

被用来描述振幅衰减。k表示空间波数。系数值g(x,z)，h(x,z)是与位置有关的函数，其值的大小取决于纵波速度和品质因子的大小。

基于方程(3)即可以求解得到g和h的值，例如可用多元线性回归方法进行求解得到g和h的值。

基于公式(1)，(2)，(3)，并将频率域转换到时间域，可得以下时间域粘滞性介质声波方程。

其中，p表示时间域波场，t表示传播时间，

表示拉普拉斯算子，f表示时间域震源项，

可描述地震波在粘滞性介质中传播时的振幅衰减现象，称为衰减项；

可描述地震波在粘滞性介质中传播时的频散现象，称为频散项。因为在波动方程(4)中衰减项和频散项是分离的，因此当舍弃一项为零而保留另外一项时，可分别模拟地震波在粘滞性介质中传播时的衰减和频散特征。例如，当我们忽略相位频散而仅仅考虑振幅衰减时，我们可以得到如下衰减域波动方程：

同样，当我们忽略振幅衰减而仅仅考虑相位频散时，我们可以得到如下频散域波动方程：

为方便数值计算，可将公式(4)做如下变换。

其中，q表示时间域辅助波场，仅作为计算方便，没有具体物理意义。此时，方程(7)可用时间二阶空间高阶差分格式求解，方程(8)可用伪谱法进行求解。因此，基于方程(7)和(8)，得到如下的粘滞性声波方程波场传播算子为：

其中，u_p表示应力场值，u_q表示辅助应力场值，i、j分别表示纵横网格点位置，k表示离散的时间值，Δx表示横向网格间距，Δz表示纵向网格间距，Δt表示时间采样间隔，a₀和a_nn表示差分系数，FFT表示快速傅里叶变换，FFT^-1表示快速傅里叶变换的反变换，k_x表示横向波数值，k_z表示纵向的波数值，N表示差分阶数，当N＝3时就是6阶差分精度，m和n是计算时的中间变量，无实际物理意义，nn表示差分步长，g，h是用粘滞性参数v₀和Q表示的值，其大小在方程(3)中确定。m和n是计算时的中间变量，无实际物理意义。

前述部分对于本说明书实施例中所采用的波场传播算子进行了解释与说明，具体的使用方式如图1所示，图1为本说明书实施例所提供的一种粘滞性介质中声波地震数据正演模拟方法的流程示意图，包括：

S101，获取初始的参数场，所述参数场包括介质纵波速度场v₀和品质因子Q。

Q为品质因子，Q为一种描述介质对地震波衰减强弱的量；具体来说，可以使用野外地质露头、钻井、物探、测井等资料绘制地质地球物理模型，并对其填充数值得到初始粘滞性参数场。

S103，将所述参数场进行矩形的网格划分，确定每个网格的位置，根据所述参数场和所述网格的位置，确定每个网格的模拟系数g和h。

可进行规则网格或不规则网格划分，一般是规则网格划分即可，即将网格划分为矩形。在实际工业应用中，考虑计算成本，因此需要满足更小的计算量，更小内存需求。为实现高效、稳定、高精度的正演模拟结果，网格间距Δx、Δz与时间采样间隔Δt的选取，应在满足稳定性条件下，且数值频散最小的原则。此处需要满足的稳定性条件为：

其中，v₀表示速度场，Δt为时间采样间隔，Δd为网格间距，Δd包括横向网格间距Δx和纵向网格间距Δz，即Δx和Δz都可以设置为满足该条件。

S105，确定炮点和检波点位置，根据所述g和h，采用预设方程确定检波点的应力场值。

该预设方程即为前述的方程(9)所对应的粘滞性声波方程波场传播算子：

其中，u_p表示应力场值，u_q表示辅助应力场值，i、j分别表示纵横网格点位置，k表示离散的时间值，Δx表示横向网格间距，Δz表示纵向网格间距，Δt表示时间采样间隔，a₀和a_nn表示差分系数，FFT表示快速傅里叶变换，FFT^-1表示快速傅里叶变换的反变换k_x表示横向波数值，k_z表示纵向的波数值，N表示差分阶数，m和n是计算时的中间变量，无实际物理意义，N表示差分阶数，可以根据实际精度需要进行设置，当N＝3时就是6阶差分精度,nn表示差分步长，g，h是用粘滞性参数v₀和Q表示的值，其大小在已经被预先确定。

表示横向第i个位置、纵向第j个位置的网格在第k离散时刻的应力场值。

波场传播算子在前文已经进行了详细说明。在实际应用，可以通过功能模块或者算法模块的形式进行预设即可。炮点和检波点位置可设置在网格的任一位置，一般设置在第一层网格上(即地表)。

在一种实施例中，在数值模拟过程中，可以应用时间二阶空间十阶的差分格式，实现基于新推导的粘滞性声波方程的数值模拟，得到一个关于波场传播的应力场数据。

与现有技术相比，本发明通过将振幅衰减项与相位频散项分离，相比传统的衰减项与频散项耦合在一起的粘声波方程正演模拟，可方便地实现粘声逆时偏移成像。相比传统计算方法，可大大提高计算效率、减少计算内存，实际应用也证明了粘声波模拟很好地应用于地震勘探，为逆时偏移及全波形反演提供一种较好的正演模拟方案。

如图2所示，图2为本说明书提供的的波场模拟结果对比示意图。其中，a部分为既不含衰减也不含频散的波场模拟结果(即声波模拟结果)，b部分为只含频散的波场模拟结果，c部分为只含衰减的波场模拟结果，d部分为同时含有衰减和频散的波场模拟结果。从图2中可以看出，本发明的粘滞性波动方程可以分别实现振幅衰减及相位频散模拟，也可以实现振幅衰减及相位频散同时模拟。图中，横坐标为长度x，纵坐标为深度z。

本说明书实施例所提供的方案，通过解耦振幅衰减与相位频散项，可分别实现振幅衰减与相位频散模拟。相比振幅衰减与相位频散耦合在一起的粘滞性介质波场正演模拟，可方便地实现逆时偏移成像，因为在做偏移成像时仅需对振幅衰减进行补偿，而相位频散不需要变化。实际应用也证明了解耦的粘滞性声波正演模拟很好地应用于地震勘探。同时，本实施例提供的时间域正演模拟方案相比频率域正演模拟，可极大提高计算效率，节省计算内存。

为了验证本专利粘滞性声波正演模拟方法对复杂模型的准确性，将本发明提出的方案用于复杂气烟囱模型中。该模型的粘滞性参数如图3所示，a部分为纵波速度参数，b部分为品质因子参数，可从图中看出模型中部存在一个较大的气烟囱，该气烟囱对穿过的地震波具有较强的振幅衰减及相位频散影响。图3为本说明书实施例所提供的一种模型的粘滞性参数与位置的关系示意图。图中，横坐标为长度x，纵坐标为深度z。

此处选择传统振幅衰减与相位频散耦合在一起的粘滞性声波方程作参考。图4为本说明书实施例所提供的方法与传统方法所得到的炮记录对比的示意图，图中，横坐标为长度x，纵坐标为时间t。图5为图4中单道炮记录的波形对比示意图。可以看出本专利提出的方案在波场走时和相位上仍与传统方法非常接近。通过气烟囱模型测试，证明了本方案所提出的粘滞性介质正演模拟方法的准确性。

对应的，本申请实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其中，所述处理器执行所述程序时实现前述的粘滞性介质中声波地震数据正演模拟方法。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置、设备和介质类实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可，这里就不再一一赘述。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置、设备和介质类实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可，这里就不再一一赘述。

上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤或模块可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。

Claims (4)

1.一种粘滞性介质中声波地震数据正演模拟方法，包括：

获取初始的参数场，所述参数场包括介质纵波速度场v₀和品质因子Q；

将所述参数场进行矩形的网格划分，确定每个网格的位置，根据所述参数场和所述网格的位置，根据如下方程确定每个网格的模拟系数g和h：

其中，x表示横坐标，z表示纵坐标，ω表示角频率，v₀代表介质纵波速度，l表示虚数单位，k表示空间波数，v表示粘滞性介质的复数值波场传播速度；

确定炮点和检波点位置，根据所述g和h，采用如下预设方程确定检波点的应力场值：

其中,i表示横向第i个网格，j表示纵向第j个网格，u_p表示应力场值，u_q表示辅助应力场值，k表示离散的时间值，Δx表示横向网格间距，Δz表示纵向网格间距，Δt表示时间采样间隔，a₀和a_nn为预设的差分系数，FFT表示快速傅里叶变换，FFT^-1表示快速傅里叶变换的反变换，k_x表示横向波数值，k_z表示纵向的波数值，N表示差分阶数，m和n是计算时的中间变量，无实际物理意义，N表示差分阶数，nn表示差分步长，

表示横向第i个位置、纵向第j个位置的网格在第k离散时刻的应力场值。

2.如权利要求1所述的方法，根据如下方程确定每个网格的模拟系数g和h，包括：

根据所述方程采用多元线性回归的方式确定每个网格的模拟系数g和h。

3.如权利要求1所述的方法，将所述参数场进行矩形的网格划分，包括：

划分得到的网格之间的间距Δd满足条件：

所述Δd包括Δx和Δz。

4.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其中，所述处理器执行所述程序时实现如权利要求1至3任一所述的方法。

Images