CN113866823B - 一种粘声各向异性介质中的正演成像方法 - Google Patents

Abstract

本说明书实施例公开了一种粘声各向异性介质中的正演成像方法。发明通过波动方程系数描述了传播过程中的相位频散或者振幅衰减影响的性质，可很方便的实现振幅补偿逆时偏移成像。相比传统的波动方程，具有相同的数值模拟精度，且在极其复杂的介质中仍能稳定模拟波场传播。此外，新方程具有解耦的振幅衰减项与相位频散项，相比传统波动方程，在实现逆时偏移成像时更加稳定，可很好地应用于地震勘探。

Description

一种粘声各向异性介质中的正演成像方法

技术领域

本说明书涉及勘探地球物理学领域，尤其涉及一种粘声各向异性介质中的正演成像方法。

背景技术

由于实际地球介质岩石间存在内摩擦力及岩石内充填的流体使得地下介质通常表现出粘滞性。当地震波在含粘滞性介质中传播时会表现出振幅衰减及相位频散特性。若在地震成像时，不考虑这些粘滞性影响，则得到的成像结果振幅较弱，地震波同相轴连续性差。此外，连续沉积的地层又通常表现出各向异性特性，地震波在各向异性介质中传播时会表现出传播速度沿某一方向更快的特性。若在偏移成像时不对这些各向异性特性进行校正，则偏移成像结果会出现成像位置不准，及能量发散的情况，从而产生成像噪声。在很多地层中这两种性质总是同时存在的，传统的粘声各向异性波动方程在模拟复杂介质中波传播时表现出不稳定。

因此，需要一种可同时考虑粘滞性及各向异性特性的，且能稳定模拟复杂介质中地震波传播的成像方法。

发明内容

本发明的目的在于，提供一种能稳定模拟复杂介质中地震波传播的成像方法。

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

一种粘声各向异性介质中的正演成像方法，包括：

获取初始参数场，其中，所述初始参数场中的参数包括粘滞性参数和各向异性参数，所述粘滞性参数包括v_p、Q、ε_Q和δ_Q，其中，v_p表示介质中纵波的传播速度，Q表示品质因子，ε_Q和δ_Q表示衰减各向异性参数对，所述各向异性参数包括ε，δ，φ，其中，ε，δ表征速度各向异性强度，φ为各向异性倾角参数

对所述参数场进行网格剖分；

针对每个网格点，确定与粘滞性参数和各向异性参数相关的波动方程系数a₁、a₂、a₃、b₁、b₂和b₃，其中，

其中，γ_ij＝arctan(1/Q_ij)/π；

采用预定的波场传播算子，根据所述波动方程系数和所述初始参数场进行计算，生成在每个网格点的纵波场值，其中，所述波场传播算子的形式为：

其中，p表示纵波场值，t表示波场传播时间，表示拉普拉斯算子，f表示震源项，x为横坐标，z为纵坐标。

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

与现有技术相比，本发明通过波动方程系数描述了传播过程中的相位频散或者振幅衰减影响的性质，可很方便的实现振幅补偿逆时偏移成像。相比传统的波动方程，具有相同的数值模拟精度，且在极其复杂的介质中仍能稳定模拟波场传播。此外，新方程具有解耦的振幅衰减项与相位频散项，相比传统波动方程，在实现逆时偏移成像时更加稳定，可很好地应用于地震勘探。

附图说明

图1为本说明书实施例所提供的一种粘声各向异性介质中的正演成像方法的流程示意图；

图2为本说明书实施例提供的粘声各向异性介质中解耦的振幅衰减与相位频散波场模拟示意图；

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

图4为本说明书实施例所提供的方法与传统方法所得到的波场快照对比的示意图；

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

图6为本说明书实施例所提供的方法与传统方法所得到的波场快照对比的示意图。

具体实施方式

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

首先需要对于本申请模拟过程中使用的波场传播算子进行说明。

已知粘声VTI介质(VTI表示不含倾角各向异性影响的特殊的各向异性介质)中复数值精确相速度公式可表示为：

其中，V_p表示qP波复数值相速度，θ表示相位角，M_i,j表示复数值刚度系数，其可表示为：

其中，ω表示角频率，ω₀表示参考角频率，i为虚部单位，γ_ij＝arctan(1/Q_ij)/π是无量纲的量。对于任何的正值Q_ij，存在0＜γ_ij＜0.5，质量因子Q_ij可表示为：

Q₃₃＝Q_p, (5)

其中，ε_Q和δ_Q表示衰减相关的Thomsen参数对，C_ij为弹性刚度系数，其可表示为：

其中，v_p表示沿着对称轴方向的纵波速度，ε和δ表示Thomsen参数。因为公式(1)中存在根号项，无法直接使用方程(1)进行推导，继而得到显式的波动方程表达式。对公式(1)应用一阶泰勒近似可去掉根号项，Z表示一个复值函数，因此可得：

又因为存在关系sinθ＝V_p(θ)k_x/ω和cosθ＝V_p(θ)k_z/ω，其中k_x和k_z为沿着水平方向和垂直方向的空间波数，经过一些数学运算可得空间波数域粘声各向异性波动方程为：

其中，

方程(10)为粘声VTI介质空间波数域波动方程，通过在公式(10)中引入一个各向异性倾角φ，可得粘声TTI介质(TTI介质表示含有倾角各向异性影响的各向异性介质)空间波数域波动方程为：

将公式(11)变换到时间空间域可得以下时间域粘声TTI介质波动方程为：

其中，p表示纵波场值，t表示波场传播时间，表示拉普拉斯算子，f表示震源项。注意，与系数a₁，a₂，a₃相关的波动方程项可以描述地震波在粘滞性介质中传播过程中产生的相位频散影响，与系数b₁，b₂，b₃相关的波动方程项可以描述地震波在粘滞性介质中传播过程中产生的振幅衰减影响。当我们假设b₁，b₂，b₃的值均为0时，方式(12)描述的是波传播过程中的相位频散特征。当我们假设γ₁₁，γ₁₃，γ₃₃的值均为0时，方式(12)描述的是波传播过程中的振幅衰减特征。正是因为波动方程(12)存在的这种可单独描述相位频散或者振幅衰减影响的性质，该方程可很方便的实现振幅补偿逆时偏移成像。

在利用方程(12)模拟波场传播时，为方便利用计算机进行数值离散，我们将方程(12)等价为如下方程：

其中，

q₁＝FFT^-1(ln(k)FFT(p)),q₂＝FFT^-1((ln(k))²FFT(p)),

需要说明的是这些q变量仅为辅助方程计算引入的量，其没有实际的物理意义。

此时对方程(13)使用时间二阶，空间高阶差分格式进行离散，从而得到粘声各向异性介质中的地震波传播算子为：

其中，

其中，p为地震波应力场值，q为辅助地震波应力场值，i和j分别表示横向和纵向网格点的位置，k表示时间离散，Δx表示离散网格的横向间距，Δz表示离散网格的纵向间距，Δt表示差分离散的时间采样间距，a₀和a_nn表示有限差分离散的差分系数，FFT表示快速傅里叶变换算符，FFT^-1表示快速傅里叶变换的反变换算符，k_x和k_z分别表示沿着横向和纵向的波数，N表示有限差分离散的差分阶数。m，n，u，mx，nz，mz，ux，uz表示定义的有限差分法数值计算时的中间变量，其没有实际的物理意义，当N值取5时，代表了10阶空间差分精度。

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

S101，获取初始参数场，其中，所述初始参数场中的参数包括粘滞性参数和各向异性参数。粘滞性参数包括v_p和Q，ε_Q，δ_Q，其中，v_p表示介质中纵波的传播速度，Q表示品质因子，Q的大小决定了地震波能量衰减强弱，ε_Q，δ_Q表示衰减各向异性参数对，用于表征衰减各向异性的强度，各向异性参数包括ε，δ，φ，其中，ε，δ为Thomsen各向异性参数对，用于表征速度各向异性强度，φ为各向异性倾角参数，描述了速度各向异性的方向特性；

具体来说，就是根据野外地质考察，机械钻井分析，地球物理反演及测井探测等方法获取的资料绘制地质地球物理模型，然后对该模型填充粘滞性及各向异性参数值可得所需的用于计算机数值模拟的参数场。

S103，对所述初始参数场进行网格剖分；

根据实际需要可对模型进行规则网格或不规则网格划分，一般对需要处理起伏地表的情况进行非规则网格剖分，对水平地表的情况，简单的规则网格划分即可，即将网格划分为矩形。在实际工业生产中，为了减少计算成本，那么计算量越小，内存需求越少越好，也就意味着剖分的网格数越少越好。同时，为了保证有限差分算法的稳定模拟，时间采样间隔与网格间距需要满足一定的条件。因此，为保证高效，精确，且稳定的有限差分数值模拟，在满足稳定性模拟，且数值频散最小的条件下，速度、时间采样间隔、网格间距应满足以下关系：其中，v_pmax表示最大的纵波速度，Δt表示时间采样间隔，Δd表示划分的网格间距。

S105，针对每个网格点，确定与粘滞性参数和各向异性参数相关的波动方程系数，这些系数包括：a₁，a₂，a₃，b₁，b₂，b₃，可根据方程(10)确定，此处不再赘述。

S107，采用预定的波场传播算子，根据所述波动方程系数和所述初始参数场进行计算，生成在每个网格点的纵波场值。

即根据公式(12)中所给定的粘声各向异性介质地震波场传播算子，采用已计算得到的与空间网格位置有关的波场传播方程的系数和所述粘滞性及各向异性参数场，在确定布设的炮点及检波点位置后，使用计算机进行有限差分数值模拟，从而记录检波器位置的地震波应力场值。

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

如前所述，在数值模拟过程中，使用时间二阶、空间十阶差分格式求解方程(14)，并使用伪谱法计算方程(23)-方程(33)，实现基于新推导的粘声各向异性波动方程的数值模拟，得到一个关于地震波场传播的应力场数据(即纵波场值)。

图2为本说明书提供的模拟不含振幅衰减与相位频散影响的、只含相位频散影响的、只含振幅衰减影响的、同时含有振幅衰减与相位频散影响的波场模拟对比示意图。其中，图2中的a部分为模拟的既不含振幅衰减也不含相位频散的波场，即各向异性声波模拟结果，如图中所示，波场传播速度含有方向性，沿某一方向传播速度更快。图2中的b部分为只含频散影响的地震波场模拟结果，相比图2中的a部分振幅一致，但波前面存在滞后效应。图3中的部分为只含振幅衰减影响的地震波场模拟结果，与图2中的a部分相比其波前面一致，但存在振幅衰减。图2中的d部分为同时含有相位频散和振幅衰减的地震波场模拟结果，与图2中的a部分相比既有波前面滞后，又存在振幅能量降低。根据上述分析可知，本发明提出的粘声各向异性波动方程可以模拟地下波场传播过程中的粘滞性及各向异性影响，同时可单独进行粘滞性引起的振幅衰减影响及相位频散影响模拟，这种特点便于实现衰减补偿的逆时偏移成像。图中，横坐标为长度x，纵坐标为深度z。

为了验证本专利粘声各向异性正演模拟方法在复杂模型中地震波场模拟的准确性，首先构建一个复杂的气烟囱模型，并将其应用于本发明所提出的方案。该模型含有较强的粘滞性及各向异性特征，其参数模型如图3所示，当地震波通过含粘滞性区域时，地震波会产生振幅衰减及相位频散特点，当地震波穿过含各向异性区域时，会由于各向异性的影响产生波场传播速度沿某一方向传播更快的特点。图3为本说明书实施例所提供的一种复杂模型的粘滞性及各向异性参数与位置的关系示意图。图中，横坐标为长度x，纵坐标为深度z。

此处选择传统的粘声各向异性波动方程作参考，其具有较高的地震波场模拟精度，但是其振幅衰减影响与相位频散影响耦合在一起的。图4为本说明书实施例所提供的方法与传统方法所得到的某一时刻波场快照对比的示意图，图中，横坐标为长度x，纵坐标为深度z。其中，图4中的a部分为传统方法模拟的波场快照，图4中的b部分为本说明书实施例模拟的波场快照，图4中的c部分为用红色曲线表示的图4中的b部分叠加在用黑色曲线表示的图4中的a部分上得到的波场叠加结果。黑色的曲线和红色曲线匹配很好，从而证明了本专利所提出的粘声各向异性介质正演模拟方法的准确性。

为了验证本专利粘声各向异性介质正演模拟方法在极其复杂模型中地震波场模拟的稳定性，首先根据BP2007模型新建一个含粘滞性及各向异性的参数模型，并将其应用于本发明所提出的方案。该模型含有剧烈倾角变化的各向异性倾角模型。当使用传统粘声各向异性地震波方程模拟地震波在该含剧烈倾角变化的模型中传播时，会产生随时间以指数形式增加的高频噪声，最终产生不稳定。图5为本说明书实施例所提供的第二种模型的粘滞性及各向异性参数与位置的关系示意图。图中，横坐标为长度x，纵坐标为深度z。

此处仍然选择传统的粘声各向异性波动方程作参考，使用其对含复杂各向异性倾角的模型中的波场进行模拟时易产生不稳定。图6为本说明书实施例所提供的方法与传统方法所得到的某一时刻波场快照对比的示意图，图中，横坐标为长度x，纵坐标为深度z。其中，图6中的a部分为传统方法模拟的波场快照，图6中的b部分为本说明书实施例模拟的波场快照。图6中的a部分中含有高频不稳定，如箭头所示，而图6中的b部分中没有高频不稳定产生，从而证明了本专利所提出的粘声各向异性介质正演模拟方法在极其复杂模型中进行数值模拟的稳定性。

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

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

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

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

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

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

Claims (3)

1.一种粘声各向异性介质中的正演成像方法，包括：

获取初始参数场，其中，所述初始参数场中的参数包括粘滞性参数和各向异性参数，所述粘滞性参数包括v_p、Q、ε_Q和δ_Q，其中，v_p表示介质中纵波的传播速度，Q表示品质因子，ε_Q和δ_Q表示衰减各向异性参数对，所述各向异性参数包括ε，δ，φ，其中，ε，δ表征速度各向异性强度，φ为各向异性倾角参数

对所述初始参数场进行网格剖分；

针对每个网格点，确定与粘滞性参数和各向异性参数相关的波动方程系数a₁、a₂、a₃、b₁、b₂和b₃，其中，

其中，

γ_ij＝arctan(1/Q_ij)/π；

采用预定的波场传播算子，根据所述波动方程系数和所述初始参数场进行计算，生成在每个网格点的纵波场值，其中，所述波场传播算子的形式为：

其中，p表示纵波场值，t表示波场传播时间，表示拉普拉斯算子，f表示震源项，x为横坐标，z为纵坐标。

2.如权利要求1所述的方法，对所述参数场进行网格剖分，包括：

划分得到的网格之间的间距Δd满足条件：其中，v_pmax表示最大的纵波速度，Δt表示时间采样间隔。

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