基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法

摘要

本发明提供了一种基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法，属于石油勘探领域，包括如下步骤：将速度和品质因子进行不规则网格剖分，并映射到曲网格坐标系下，将传统的一阶黏声速度‑应力方程映射到曲网格坐标系下，并应用曲网格坐标系下的一阶黏声速度‑应力方程推导了一种同一坐标系下不含记忆变量的二阶黏声拟微分方程；通过新的黏声方程可以更好地控制黏声介质造成的振幅损失与相位频散，以更准确地模拟地震波在存在起伏地表的黏声介质中传播的精度，在方程求解时，提出了一种混合空间偏导数差分方法进行求解，并提出了一种新的适用于起伏地表曲网格坐标系下的边界条件吸收人工边界反射。

说明书

技术领域

本发明属于石油勘探领域，具体涉及一种基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法。

背景技术

地球介质中广泛存在着黏滞性，在进行地震波正演模拟时，必须考虑地层黏滞性对地震波传播的振幅衰减影响和相位频散影响。随着油气勘探开发的深入，所面临的地质构造越来越复杂，特别是复杂的起伏地表构造为地震勘探带来了巨大挑战。针对我国探区地质构造的复杂性，需要针对复杂起伏地表构造的地震波传播规律进行研究，从而更好地为油气藏的识别、描述、挖潜等提供更加重要的地球物理依据。

传统的黏声介质正演模拟基于GSLS模型的黏声方程进行波场计算的，但该方法无法准确处理剧烈起伏地表的影响，同时传统黏声方程存在记忆变量而且无法将黏滞性介质对地震波的衰减和相位频散的影响分离开，因此无法准确地补偿振幅衰减、校正相位频散的影响。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法，设计合理，克服了现有技术的不足，具有良好的效果。

为了实现上述目的，本发明采用如下技术方案：

基于黏声拟微分方程的黏声起伏地表正演模拟系统，包括输入模块、坐标变换模块、曲坐标下黏声方程计算模块、黏声波场延拓模块、边界吸收模块和输出模块；

输入模块，被配置为用于输入含有起伏地表速度模型和品质因子Q模型；

坐标变换模块，被配置为用于在笛卡尔坐标系下进行曲网格剖分，并转换到曲网格坐标系下的水平地表速度模型和品质因子模型；

曲坐标下黏声方程计算模块，被配置为用于将一阶黏声速度-应力方程坐标变换到曲网格坐标系下的一阶黏声速度-应力方程，得到曲网格坐标系下的不含记忆变量的二阶黏声拟微分方程；

黏声波场延拓模块，被配置为用于利用混合空间偏导数差分方法求解波场；

边界吸收模块，被配置为用于利用起伏地表曲网格坐标系下的边界条件吸收人工边界反射；

输出模块，被配置为用于输出笛卡尔坐标系下的炮记录和波场快照。

此外，本发明还提到一种基于黏声拟微分方程的黏声起伏地表正演模拟方法，该方法采用如上所述的基于黏声拟微分方程的黏声起伏地表正演模拟系统，包括如下步骤：

步骤1：通过输入模块，输入含有起伏地表速度模型和品质因子Q模型，并建立观测系统；

步骤2：通过坐标变换模块，在笛卡尔坐标系下进行曲网格剖分，并采用如下所示的坐标变换方程转换到曲网格坐标系下的水平地表速度模型和品质因子模型；

其中，x和z分别表示笛卡尔坐标系下的水平坐标和垂直坐标；ξ和η分别表示曲网格坐标系下的水平坐标和垂直坐标；

步骤3：通过曲坐标下黏声方程计算模块，将笛卡尔坐标系下的一阶黏声速度-应力方程坐标变换到曲网格坐标系下的一阶黏声速度-应力方程，如公式(5)所示；

其中，v_x和v_z分别表示x分量和z分量的速度场；p是声波的压力场；t为时间；ρ表示密度；K^R和K^μ分别表示松弛模量和非松弛模量；τ_σ和τ_ε分别表示应力松弛时间和应变松弛时间，τ_σ和τ_ε由下式求得

其中，Q为品质因子，ω为角频率；

步骤4：根据曲网格坐标系下的一阶黏声速度-应力方程，推导得到曲网格坐标系下不含记忆变量的二阶黏声拟微分方程，如公式(17)所示；

其中，v₀表示零频速度，τ由下式求得

步骤5：通过黏声波场延拓模块，利用混合空间偏导数差分方法，对曲网格坐标系下的不含记忆变量的二阶黏声拟微分方程进行求解；

步骤6：通过边界吸收模块，利用起伏地表曲网格坐标系下的边界条件公式(23)，吸收人工边界反射；

其中，由如式(24)所示的变换求得

其中，ε_ξ和ε_η由式(25)求得

其中，α为衰减系数；κ和σ为辅助衰减系数；i表示虚部；

步骤7：将波场快照变换到笛卡尔坐标系下；

步骤8：通过输出模块输出笛卡尔坐标系下的炮记录和波场快照。

优选地，在步骤3中，笛卡尔坐标系下传统的基于GSLS模型的一阶黏声速度-应力方程为：

采用链式法则映射后，方程(2)可以转换成曲网格坐标系下的一阶黏声速度-应力方程，如方程(5)所示：

优选地，在步骤4中，将曲网格坐标系下的一阶黏声速度-应力方程转换到频率域：

其中，上标～表示频率域中的变量，k_ξ和k_η分别表示ξ和η方向的波数，v₀表示零频速度，根据方程(6)中的前两个方程得到：

将方程(7)和(8)代入方程(6)中的后两个方程可得：

消除方程(8)和(10)中的变量得到曲网格坐标系下的频散关系式：

整理方程(11)可得

其中，

当Q不是特别小的时候，τ＜＜1，得到如下近似式：

因此，方程(12)可近似为

方程(15)可等价地改写为

对方程(16)作傅里叶反变换，得到曲网格坐标系下不含记忆变量的二阶黏声拟微分方程：

优选地，在步骤5中，使用基于混合空间偏导数的O(Δt²,Δξ^2M)阶公式对曲网格坐标系下的不含记忆变量的二阶黏声拟微分方程进行求解，混合空间偏导数差分公式为：

其中，和分别表示傅里叶变换和傅里叶反变换；Δt为时间步长；Δξ和Δη分别为ξ和η方向的空间网格间距；c¹和c²分别表示一阶和二阶的差分系数；对于变量p，上标表示时间坐标，下标表示空间坐标；k_ξ和k_η分别表示ξ和η方向的波数；M为差分精度，递推公式为

本发明所带来的有益技术效果：

本发明能够同时考虑起伏地表影响和地下介质黏滞性影响，通过实现一种曲网格坐标系下不含记忆变量的黏声拟微分方程，既能克服山区地形剧烈起伏地表对地震波模拟的影响，又能同时准确地模拟黏声介质对地震波的衰减和相位频散的影响，有助于得到更准确的起伏地表黏声介质的地震波波场特征，克服了传统黏声正演模拟方法无法准确处理剧烈起伏地表且无法准确地补偿振幅衰减、校正相位频散的缺点，开发基于一种基于新的黏声拟微分方程的黏声起伏地表正演模拟方法，为山前带存在剧烈起伏地表及强黏弹性探区提供准确的正演模拟波场，有助于分析地震波在该探区的传播规律。

附图说明

图1为本发明的基于新的黏声拟微分方程的黏声起伏地表正演模拟方法的流程图；

图2为黏声起伏地表Marmousi模型。(a)笛卡尔坐标系下的速度模型；(b)笛卡尔坐标系下的Q模型；(c)曲网格坐标系下的速度模型；(d)曲网格坐标系下的Q模型。

图3为网格剖分图。(a)笛卡尔坐标系下的曲网格；(b)曲网格坐标系下的矩形网格。

图4为不同时刻曲网格坐标系下的波场快照。(a-c)黏声介质；(d-f)声波介质。

图5为不同时刻笛卡尔坐标系下的波场快照。(a-c)黏声介质；(d-f)声波介质。

图6为炮记录。(a)黏声介质；(b)声波介质。

图7为从黏声介质和声波炮记录中抽取的波形图。(a)直达波；(b)反射波。

图8为从黏声介质和声波炮记录中抽取的频谱图。(a)直达波；(b)反射波。

图9为本发明中基于新的黏声拟微分方程的黏声起伏地表正演模拟系统的结构示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

实施例1：

一种基于黏声拟微分方程的黏声起伏地表正演模拟系统，其结构如图9所示，包括输入模块、坐标变换模块、曲坐标下黏声方程计算模块、黏声波场延拓模块、边界吸收模块和输出模块；

输入模块，被配置为用于输入含有起伏地表速度模型和品质因子Q模型；

坐标变换模块，被配置为用于在笛卡尔坐标系下进行曲网格剖分，并转换到曲网格坐标系下的水平地表速度模型和品质因子模型；

曲坐标下黏声方程计算模块，被配置为用于将一阶黏声速度-应力方程坐标变换到曲网格坐标系下的一阶黏声速度-应力方程，得到曲网格坐标系下的不含记忆变量的二阶黏声拟微分方程；

黏声波场延拓模块，被配置为用于利用混合空间偏导数差分方法求解波场；

边界吸收模块，被配置为用于利用起伏地表曲网格坐标系下的边界条件吸收人工边界反射；

输出模块，被配置为用于输出笛卡尔坐标系下的炮记录和波场快照。

实施例2：

在上述实施例的基础上，本发明还提到一种基于黏声拟微分方程的黏声起伏地表正演模拟方法，其流程如图1所示，具体包括如下步骤：

步骤1：输入含有起伏地表速度模型和品质因子Q模型，并建立观测系统；

步骤2：笛卡尔坐标系下进行曲网格剖分，并采用如下所示的坐标变换方程转换到曲网格坐标系下的水平地表速度模型和品质因子模型；

其中，x和z分别表示笛卡尔坐标系下的水平坐标和垂直坐标；ξ和η分别表示曲网格坐标系下的水平坐标和垂直坐标；

步骤3：将笛卡尔坐标系下的一阶黏声速度-应力方程坐标变换到曲网格坐标系下的一阶黏声速度-应力方程；

笛卡尔坐标系下传统的基于GSLS模型的一阶黏声速度-应力方程为：

其中，v_x和v_z分别表示x分量和z分量的速度场；p是声波的压力场；t为时间；ρ表示密度；K^R和K^μ分别表示松弛模量和非松弛模量；τ_σ和τ_ε分别表示应力松弛时间和应变松弛时间，τ_σ和τ_ε由下式求得

其中，Q为品质因子，ω为角频率，采用链式法则映射后，方程(2)可以转换成曲网格坐标系下的一阶黏声速度-应力方程：

步骤4：得到曲网格坐标系下的不含记忆变量的二阶黏声拟微分方程；

将方程(5)转换到频率域：

其中，上标～表示频率域中的变量，k_ξ和k_η分别表示ξ和η方向的波数，v₀表示零频速度，根据方程(6)中的前两个方程得到：

将方程(7)和(8)代入方程(6)中的后两个方程可得：

消除方程(8)和(10)中的变量得到曲网格坐标系下的频散关系式：

整理方程(11)可得

其中，

当Q不是特别小的时候，τ＜＜1，得到如下近似式：

因此，方程(12)可近似为

方程(15)可等价地改写为

对方程(16)作傅里叶反变换，得到曲网格坐标系下不含记忆变量的二阶黏声拟微分方程：

步骤5：利用混合空间偏导数差分方法求解波场；

使用基于混合空间偏导数的O(Δt²,Δξ^2M)阶公式对曲网格坐标系下的不含记忆变量的二阶黏声拟微分方程进行求解，混合空间偏导数差分公式为：

其中，和分别表示傅里叶变换和傅里叶反变换；Δt为时间步长；Δξ和Δη分别为ξ和η方向的空间网格间距；c¹和c²分别表示一阶和二阶的差分系数；对于变量p，上标表示时间坐标，下标表示空间坐标；k_ξ和k_η分别表示ξ和η方向的波数；M为差分精度，递推公式为

步骤6：利用起伏地表曲网格坐标系下的边界条件吸收人工边界反射；

其中，由如下所示的变换求得

其中，ε_ξ和ε_η由下式求得

其中，α为衰减系数；κ和σ为辅助衰减系数；i表示虚部；

步骤7：将波场快照变换到笛卡尔坐标系下；

步骤8：输出笛卡尔坐标系下的炮记录和波场快照。

本发明基于新的黏声拟微分方程的黏声起伏地表正演模拟能够同时考虑起伏地表影响和地下介质黏滞性影响，通过实现一种曲网格坐标系下不含记忆变量的黏声拟微分方程，既能克服山区地形剧烈起伏地表对地震波模拟的影响，又能同时准确地模拟黏声介质对地震波的衰减和相位频散的影响，有助于得到更准确的起伏地表黏声介质的地震波波场特征。

应用实验

本发明基于新的黏声拟微分方程的黏声起伏地表正演模拟方法，应用于起伏地表黏声Marmousi模型数据，取得了理想的计算效果。输入含有起伏地表速度模型(如图2a所示)和品质因子Q模型(如图2b所示)，并建立观测系统；笛卡尔坐标系下进行曲网格剖分(如图3所示)，并转换到曲网格坐标系下的水平地表速度模型(如图2c所示)和品质因子模型(如图2d所示)；将一阶黏声速度-应力方程坐标变换到曲网格坐标系下的一阶黏声速度-应力方程；得到曲网格坐标系下的不含记忆变量的二阶黏声拟微分方程；利用混合空间偏导数差分方法求解波场；利用起伏地表曲网格坐标系下的边界条件吸收人工边界反射；将波场快照(如图4a-4c所示)变换到笛卡尔坐标系下；输出笛卡尔坐标系下的波场快照(如图5a-5c所示)和炮记录(如图6a所示)。从图中可以看出，采用本发明提出的起伏地表黏声正演模拟方法可以很好地对地震波的传播特征进行模拟，包括直达波和反射波。为了比较本发明的效果，给出了地震波在声波介质中传播的波场快照(如图4d-4f和图5d-5f所示)和炮记录(如图6b所示)。从两种介质中抽取的单道记录的波形和频谱如图7和图8所示。图中虚线所示的为本发明的黏声正演模拟的结果，实现所示的为声波正演模拟的结果。从这些结果中可以看出，本发明可以准确地模拟出地震波在存在起伏地表的黏声介质中的传播特征，振幅衰减与相位频散都可以准确地模拟出来，证明了本论文提出的基于曲网格坐标系下不含记忆变量的二阶黏声拟微分方程的黏声起伏地表正演模拟方法能够准确且稳定地模拟地震波在存在剧烈起伏地表的复杂黏声介质中的传播特征。

地球介质中广泛存在着黏滞性，在进行地震波正演模拟时，必须考虑地层黏滞性对地震波传播的振幅衰减影响和相位频散影响。随着油气勘探开发的深入，所面临的地质构造越来越复杂，特别是复杂的起伏地表构造为地震勘探带来了巨大挑战。针对我国探区地质构造的复杂性，需要针对复杂起伏地表构造的地震波传播规律进行研究，从而更好地为油气藏的识别、描述、挖潜等提供更加重要的地球物理依据。提高存在起伏地表的黏滞性介质的地震波正演模拟精度对适应国内地震勘探的发展趋势是非常必要的。

为此本发明提出了一种基于新的黏声拟微分方程的黏声起伏地表正演模拟方法，开发存在起伏地表的黏滞性介质的地震波正演技术，为后续的复杂地质构造成像、反演工作提供地震波场计算基础。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。