CN117471531B - 一种基于幂律频变q效应的vti黏弹性波方程数值模拟方法 - Google Patents

Abstract

本发明公开了一种基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，包括：S1：推导VTI介质中基于幂律频变Q效应的复数模量表达式；S2：推导二维VTI介质黏弹性波方程；S3：利用数值方法求解VTI黏弹性波方程；S4：震源与接收器设置。本发明能够较好地模拟各向异性介质中的幂律频变Q效应，对提高各向异性介质中地震资料成像与解释具有重要意义。

Description

一种基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法

技术领域

本发明涉及地震波数值模拟与成像技术领域，具体涉及一种基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法。

背景技术

黏弹性是地球介质的基本属性，地震波在黏弹性介质中传播会发生振幅衰减和速度频散，导致能量和旅行时出现变化，影响地震资料的准确成像及解释。另一方面，实验测量和野外观测表明，各向异性普遍存在于地下介质中，地震波沿不同方向的传播速度和衰减特征会出现明显差异。因此，研究地震波在各向异性黏弹性介质中的传播规律具有重要意义。当前，主流方法假设用于表征介质黏弹性强弱的物理量品质因子Q不随频率变化，并利用常Q模型(简称CQM)表征地震波在地球介质中的传播规律。但常Q假设与高温、高压或含流体介质中观测到的Q随频率变化规律不符。因此，进一步发展适用于频变Q效应的各向异性波动方程数值模拟方法，对研究地震波在真实介质中的传播规律及提高地下介质成像分辨率具有重要意义。

在勘探地震学中，地震波衰减模拟方法大致可以分为两类：传统力学模型和分数阶导数模型。前者通过串联或并列多组弹簧(弹性)和阻尼器(黏性)来模拟近似常Q效应。相比之下，分数阶导数模型由常Q模型直接推导而来。基于常Q模型，许多学者研究了地震波在黏弹性介质中的传播规律，并应用于Q补偿逆时偏移成像，在一定程度上提高了地震成像剖面的分辨率。但是上述方法主要基于常Q假设。通常认为，在地球内部高温高压环境下，品质因子Q会随频率发生变化，并符合幂律关系。忽略Q的频变特性，必然会导致模拟误差增大；特别在Q随频率变化剧烈的区域，这一问题更加突出。

另一方面，由于裂缝的定向排列或含流体影响，速度和品质因子均存在一定的各向异性特征。其中，速度各向异性表现在地震波沿不同方向传播的速度不同，直接影响地震波的旅行时间；而品质因子各向异性则影响不同方向品质因子的大小。由于品质因子对地震波传播的影响同时体现在振幅和相位两个方面，忽略品质因子Q的各向异性必然会使地震波能量和旅行时都发生变化，严重影响地震资料的解释。通过引入速度和品质因子各向异性参数，充分考虑地震波沿不同方向的传播特征，发展一种针对幂律频变Q效应的各向异性黏弹性波方程数值模拟方法对提高地震成像与解释的精度具有重要意义。

发明内容

为了克服现有技术方案各向异性黏弹性波方程数值模拟方法主要基于常Q假设、忽略了Q随频率的变化特征的不足，本发明提出了一种基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，能够较好地模拟各向异性介质中的幂律频变Q效应，对提高各向异性介质中地震资料成像与解释具有重要意义。为了实现上述技术目的，本发明采用如下技术方案：

一种基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，包括：

S1：推导VTI介质中基于幂律频变Q效应的复数模量表达式；

S2：推导二维VTI介质黏弹性波方程；

S3：利用数值方法求解VTI黏弹性波方程；

S4：震源与接收器设置。

优选地，步骤S1包括：

S11：在幂律频变Q模型中，品质因子Q会随频率发生变化，表示为：

Q(ω)＝Q₀(ω/ω₀)^χ (1)

其中Q₀为参考品质因子，ω和ω₀分别为角频率和参考角频率，χ为分数阶指数；

S12：由Kramers-Kronig关系，得到幂律频变Q模型中复数模量的近似表达式如下：

其中M₀为参考模量，i为虚数单位，复数模量的实部和虚部分别代表相位频散和振幅衰减；

S13：引入速度和Q各向异性，得到各向异性VTI介质中的复数模量的表达式如下：

其中，为各向异性介质的弹性张量，ρ为密度，c_0,ij为参考速度，Q_0,ij为参考品质因子，χ_ij为各向异性介质的分数阶指数；

S14：结合波数k与角频率ω之间的近似关系ω≈kc₀，其中k＝|k|为波数的模，复数模量又可以改写成如下形式：

其中：

在二维VTI介质中，弹性张量其中c_P0和c_S0分别为纵、横波的参考速度；C₁₁和C₁₃由速度各向异性参数ε和δ计算得到：

C₁₁＝(1+2ε)C₃₃ (6)

Q_0,33＝Q_P0和Q_0,55＝Q_S0分别为纵、横波的参考品质因子；其它参考品质因子张量Q_0,11和Q_0,13由品质因子各向异性参数ε_Q和δ_Q计算得到：

优选地，步骤S2包括：

S21：在二维模型中，一阶动量守恒方程表示为：

其中，v_x,v_z为速度分量，σ_xx,σ_zz,σ_xz为应力分量，ρ为密度，为一阶时间偏导数，和/>为一阶空间偏导数；

S22：结合应力σ_xx,σ_zz,σ_xz与应变ε_xx,ε_zz,ε_xz之间的关系：

σ_xx＝M₁₁ε_xx+M₁₃ε_zz (11)

σ_zz＝M₁₃ε_xx+M₃₃ε_zz (12)

σ_xz＝2M₅₅ε_xz (13)

其中，M₁₁、M₁₃、M₃₃和M₅₅为复数模量在时间域对应的实数模量；

结合近似公式M₁₃≈(1+δ)M₃₃-2M₅₅以及速度-应变关系可得：

优选地，步骤S3包括：

S31：针对方程4中的混合域算子和/>引入一阶泰勒展开近似上述算子对应的波数域算子，即：

S32：结合方程17和18，方程14–16改写为：

其中 和/>分别代表正、反傅里叶变换；

S33：采用有限差分求解方程19–21中的时间偏导数：

其中x＝(x,z)为空间坐标向量，t为当前时刻，Δt为时间采样间隔。

优选地，步骤S4包括：

选择震源子波和子波频率，设置震源激发位置，并设置检波器接收地震信号，完成上述参数设置即可开始数值模拟，以获得地震波场和地震记录。

与现有技术相比，本发明的有益效果：

本发明充分考虑了速度和品质因子的各向异性特征，基于幂律频变Q效应，推导了各向异性VTI介质黏弹性波方程，并求解，较好地实现了VTI介质频变Q效应的数值模拟。发明效果表现在：

(1)Q的频变特性方面：传统黏弹性波方程数值模拟方法主要基于常Q假设，不能模拟Q的频变特性。本发明直接从幂律频变Q模型出发，通过推导新的复数模量表达式(方程4)，建立新的黏弹性波方程(方程14–16)，可以直接在全频段模拟幂律频变Q效应。

(2)各向异性特征方面：通过引入速度和品质因子Q的各向异性参数，对二维VTI介质中的弹性模量和参考品质因子进行表征，可以较好地模拟各向异性黏弹性介质中地震波的传播规律与波场特征。

(3)频散与衰减解耦方面：类似传统常Q模型方法，类似传统常Q模型方法，本发明中基于幂律频变Q效应推导的复数模量表达式(方程3和4)较好地实现了对实部和虚部的分离，其中实部和虚部分别代表相位畸变和振幅衰减。据此，可以独立模拟相位与振幅的变化特征。

附图说明

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

图1(a)为复杂HessVTI模型参数-P波参考速度c_P0参数图。

图1(b)为复杂HessVTI模型参数-参考品质因子Q_P0参数图。

图1(c)为复杂HessVTI模型参数-速度各向异性参数ε参数图。

图1(d)为复杂HessVTI模型参数-速度各向异性参数δ参数图。

图2为复杂HessVTI模型1.1s时刻的波场快照图。

图3为复杂HessVTI模型地震记录图。

图4为复杂HessVTI模型单道记录对比图。

具体实施方式

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

一种基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，包括：

S1：推导VTI介质中基于幂律频变Q效应的复数模量表达式；

S11：在幂律频变Q模型中，品质因子Q会随频率发生变化，表示为：

Q(ω)＝Q₀(ω/ω₀)^χ (1)

其中Q₀为参考品质因子，ω和ω₀分别为角频率和参考角频率，χ为分数阶指数；

S12：由Kramers-Kronig关系，得到幂律频变Q模型中复数模量的近似表达式如下：

其中M₀为参考模量，i为虚数单位，复数模量的实部和虚部分别代表相位频散和振幅衰减；

S13：引入速度和Q各向异性，得到各向异性VTI介质中的复数模量的表达式如下：

其中，为各向异性介质的弹性张量，ρ为密度，c_0,ij为参考速度，Q_0,ij为参考品质因子，χ_ij为各向异性介质的分数阶指数；

S14：结合波数k与角频率ω之间的近似关系ω≈kc₀，其中k＝|k|为波数的模，复数模量又可以改写成如下形式：

其中：

在二维VTI介质中，弹性张量其中c_P0和c_S0分别为纵、横波的参考速度；C₁₁和C₁₃由速度各向异性参数ε和δ计算得到：

C₁₁＝(1+2ε)C₃₃ (6)

Q_0,33＝Q_P0和Q_0,55＝Q_S0分别为纵、横波的参考品质因子；其它参考品质因子张量Q_0,11和Q_0,13由品质因子各向异性参数ε_Q和δ_Q计算得到：

S2：推导二维VTI介质黏弹性波方程；

S21：在二维模型中，一阶动量守恒方程表示为：

其中，v_x,v_z为速度分量，σ_xx,σ_zz,σ_xz为应力分量，ρ为密度，为一阶时间偏导数，和/>为一阶空间偏导数；

S22：结合应力σ_xx,σ_zz,σ_xz与应变ε_xx,ε_zz,ε_xz之间的关系：

σ_xx＝M₁₁ε_xx+M₁₃ε_zz (11)

σ_zz＝M₁₃ε_xx+M₃₃ε_zz (12)

σ_xz＝2M₅₅ε_xz (13)

其中，M₁₁、M₁₃、M₃₃和M₅₅为复数模量在时间域对应的实数模量；

结合近似公式M₁₃≈(1+δ)M₃₃-2M₅₅以及速度-应变关系可得：

S3：利用数值方法求解VTI黏弹性波方程；

S31：针对方程4中的混合域算子和/>引入一阶泰勒展开近似上述算子对应的波数域算子，即：

S32：结合方程17和18，方程14–16改写为：

其中 和/>分别代表正、反傅里叶变换；

S33：采用有限差分求解方程19–21中的时间偏导数：

其中x＝(x,z)为空间坐标向量，t为当前时刻，Δt为时间采样间隔。

S4：震源与接收器设置；选择震源子波和子波频率，设置震源激发位置，并设置检波器接收地震信号，完成上述参数设置即可开始数值模拟，以获得地震波场和地震记录。

以复杂HessVTI模型为实施例，开展基于幂律频变Q效应的VTI黏弹性波方程数值模拟测试，其中模型大小为501×251个计算网格，横向和纵向的网格间距均为10米。具体步骤为：

如图1所示，根据HessVTI模型参数，包括纵、横波参考速度c_P0和c_S0，纵、横波参考品质因子Q_P0和Q_S0，以及速度各向异性参数ε和δ，品质因子各向异性参数ε_Q和δ_Q，c_S0＝c_P0/1.5，Q_S0＝10×(c_S0/1000)^2.2，ε_Q＝4.0×ε，δ_Q＝5.0×δ和设置时间采样间隔为1毫秒，记录长度为2秒，参考角频率为2π弧度/秒；震源采用主频为25Hz的雷克子波，设置在(2700，10)米处激发地震波，检波器设置在地表接收地震信号。依据方程19、20、21分别求解二维VTI介质黏弹性波方程14、15、16，执行数值模拟并存储地震波场和记录。数值模拟结果对比与分析：图2展示了1.1s时刻的波场快照。其中图2(a)和2(b)为常Q方法计算结果，图2(c)和2(d)为幂律频变Q效应数值模拟方法的结果。其中，图2(a)和2(c)为v_x波场分量，图2(b)和2(d)为和v_z波场分量。对比可见，当忽略Q的频变效应时，图2(a)和2(c)之间、图2(b)和2(d)之间存在明显差异，其中图2(a)和2(c)的常Q方法计算波场的振幅衰减更严重。图3展示了对应的地震记录，同样可以观察到上述现象。此外，图4分别从图3地震记录中横向位置x＝1500和3500米处抽取了单道，并对比不同方法计算结果的振幅与相位差异。由图可见，常Q方法(χ＝0)计算的曲线振幅比本发明所述方法计算的曲线振幅更小，相位超前。由此说明常Q方法计算结果遭受的吸收衰减更严重，而本发明所述方法充分考虑了Q的频变效应，其实际模拟的Q值比参考Q值Q₀更大，吸收衰减比常Q方法模拟结果偏弱。本发明所述方法充分考虑了Q的频变特性与地下实际介质的各向异性特征，能够更加准确地模拟地震波在地下介质中的传播规律与波场特征，对实际生产具有较好的指导意义。/>

Claims (4)

1.一种基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，包括：

S1：推导VTI介质中基于幂律频变Q效应的复数模量表达式；

步骤S1包括：

S11：在幂律频变Q模型中，品质因子Q会随频率发生变化，表示为：

Q(ω)＝Q₀(ω/ω₀)^χ (1)

其中Q₀为参考品质因子，ω和ω₀分别为角频率和参考角频率，χ为分数阶指数；

S12：由Kramers-Kronig关系，得到幂律频变Q模型中复数模量的近似表达式如下：

其中M₀为参考模量，i为虚数单位，复数模量的实部和虚部分别代表相位频散和振幅衰减；

S13：引入速度和Q各向异性，得到各向异性VTI介质中的复数模量的表达式如下：

其中，i,j＝1,2,3,4,5,6，为各向异性介质的弹性张量，ρ为密度，c_0,ij为参考速度，Q_0,ij为参考品质因子，χ_ij为各向异性介质的分数阶指数；

S14：结合波数k与角频率ω之间的近似关系ω≈kc₀，其中k＝|k|为波数的模，复数模量又可以改写成如下形式：

其中：

在二维VTI介质中，弹性张量其中c_P0和c_S0分别为纵、横波的参考速度；C₁₁和C₁₃由速度各向异性参数ε和δ计算得到：

C₁₁＝(1+2ε)C₃₃ (6)

Q_0,33＝Q_P0和Q_0,55＝Q_S0分别为纵、横波的参考品质因子；其它参考品质因子张量Q_0,11和Q_0,13由品质因子各向异性参数ε_Q和δ_Q计算得到：

S2：推导二维VTI介质黏弹性波方程；

S3：利用数值方法求解VTI黏弹性波方程；

S4：震源与接收器设置。

2.根据权利要求1所述的基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，其特征在于，步骤S2包括：

S21：在二维模型中，一阶动量守恒方程表示为：

其中，v_x,v_z为速度分量，σ_xx,σ_zz,σ_xz为应力分量，ρ为密度，为一阶时间偏导数，/>和为一阶空间偏导数；

S22：结合应力σ_xx,σ_zz,σ_xz与应变ε_xx,ε_zz,ε_xz之间的关系：

σ_xx＝M₁₁ε_xx+M₁₃ε_zz (11)

σ_zz＝M₁₃ε_xx+M₃₃ε_zz (12)

σ_xz＝2M₅₅ε_xz (13)

其中，M₁₁、M₁₃、M₃₃和M₅₅为复数模量在时间域对应的实数模量；

结合近似公式M₁₃≈(1+δ)M₃₃-2M₅₅以及速度-应变关系可得：

3.根据权利要求1所述的基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，其特征在于，步骤S3包括：

S31：针对方程4中的混合域算子和/>引入一阶泰勒展开近似上述算子对应的波数域算子，即：

S32：结合方程17和18，方程14–16改写为：

其中 和/>分别代表正、反傅里叶变换；

S33：采用有限差分求解方程19–21中的时间偏导数：

其中x＝(x,z)为空间坐标向量，t为当前时刻，Δt为时间采样间隔。

4.根据权利要求1所述的基于幂律频变Q效应的VTI黏弹性波方程数值模拟方法，其特征在于，步骤S4包括：

选择震源子波和子波频率，设置震源激发位置，并设置检波器接收地震信号，完成上述参数设置即可开始数值模拟，以获得地震波场和地震记录。