CN114139335A - 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法 - Google Patents

基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法 Download PDF

Info

Publication number
CN114139335A
CN114139335A CN202111166518.9A CN202111166518A CN114139335A CN 114139335 A CN114139335 A CN 114139335A CN 202111166518 A CN202111166518 A CN 202111166518A CN 114139335 A CN114139335 A CN 114139335A
Authority
CN
China
Prior art keywords
model
discrete
relaxation time
lattice boltzmann
wave field
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111166518.9A
Other languages
English (en)
Inventor
夏木明
周辉
姜春涛
唐瑾璇
王灿云
杨长春
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum Beijing
Institute of Geology and Geophysics of CAS
Original Assignee
China University of Petroleum Beijing
Institute of Geology and Geophysics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum Beijing, Institute of Geology and Geophysics of CAS filed Critical China University of Petroleum Beijing
Priority to CN202111166518.9A priority Critical patent/CN114139335A/zh
Publication of CN114139335A publication Critical patent/CN114139335A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法。所述方法包括：根据待模拟区域的结构特征，建立所述待模拟区域的数学模型；建立基于单松弛时间格子玻尔兹曼方程的离散模型，并根据所述离散模型中的物理参数构建各个声波波场参数的迭代计算格式；对所述离散模型中各个预设方向的波场参数赋初值，定义震源函数，并制定波阻抗界面处理方法、流‑固边界条件以及计算区域外部边界条件；根据所述迭代计算格式，利用单松弛时间格子玻尔兹曼方程，迭代更新计算得到预设采样时刻的声波波场参数的数值。利用本发明中各个实施例，可以有效提高粘滞声波模拟的准确度。

Description

基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法
技术领域
本发明涉及地球物理勘探技术领域，特别涉及一种基于单松弛时间格子玻尔兹曼模型的 粘滞声波模拟方法。
背景技术
声波波数值模拟，是根据已知的传播介质的几何结构和相应的物理参数，计算波场在各 种介质中的传播规律，并计算得到在各观测点所观测到的波场参数数值。声波波场数值模拟， 不仅是地质结构复杂地区的声波测井与地震勘探资料采集、处理和解释的有效辅助手段，也 是地球物理反演和成像的基础。
现有技术中，常用的声波波场数值模拟方法，比如有限差分方法(FDM)、有限元法等， 是基于波动方程描述波动现象在介质中的传播情况，通过求解波动方程，得到声波波场参数 数值。在求解波动方程的过程中，通常需要将波动方程离散化为代数方程进行求解。但是， 利用上述方法计算波场参数数值时，往往会受到波动方程本身的连续性假设条件的限制，而 波动方程的这些假设条件与实际情况往往存在一定偏差。这样就会导致模拟得到的波场与真 实波场存在较大偏差，进而导致求得的波场参数的数值与真实值之间存在较大的误差。
现有技术中至少存在如下问题：计算波场参数数值时，往往会受到波动方程本身一些连 续性假设条件的限制。由于波动方程的假设条件与实际情况存在偏差，会导致模拟得到的波 场与真实波场存在较大偏差，使求得的波场参数数值与真实值之间存在较大的误差，并进一 步导致地球物理反演所得地质模型/结构的精度不够高，这限制了其在实际生产中的应用。
发明内容
为了解决现有技术存在的问题，本发明引入格子玻尔兹曼方法(LBM)，并在此基础上 开发出一套新颖的粘滞声波波场数值模拟方法。格子玻尔兹曼方法是一种以当代统计物理学 为理论基础的数值模拟方法，它是一种求解纳维-斯托克斯方程的简化数值模拟方法。通过追 踪介质微观尺度上大量离散粒子的相互作用，实现模拟宏观层面上的复杂物理现象。该方法 稳定性好、计算精度高、物理意义明确、不用考虑非线性项、算法简单易实现且扩展迁移性 强、具备天然的并行特性以及内部边界条件处理灵活等优势，因此，其在计算流体力学、气 动声学、热力学、电磁学等领域都有了一定的发展及应用。
具体的，本发明提供了一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法，包 括：S1根据待模拟区域的结构特征，建立所述待模拟区域的数学模型；S2建立基于单松弛 时间格子玻尔兹曼方程的离散模型，并根据所述离散模型中的物理参数构建各个声波波场参 数的迭代计算格式；S3对所述离散模型中各个预设方向的波场参数赋初值，定义震源函数， 并制定波阻抗界面处理方法、流-固边界条件以及计算区域外部边界条件；S4根据所述迭代 计算格式，利用单松弛时间格子玻尔兹曼方程，迭代更新计算得到预设采样时刻的声波波场 参数的数值。
优选的，所述数学模型可以包括待模拟区域介质的声波速度模型、密度模型、品质因子 模型等物理参数模型。
优选的，离散模型为DdQq模型，所述DdQq模型表示d维空间q个离散速度的离散格子玻 尔兹曼模型。
优选的，单松弛时间格子玻尔兹曼方程为：
Figure BDA0003291534740000021
式中，fi(x,t)表示x位置、t时刻、i方向的粒子数密度；fi (eq)(x,t)表示x位置、t时刻、 i方向的平衡态粒子数密度；τ表示松弛时间；ci表示i方向的离散速度。
优选的，所述震源函数包括Ricker子波函数。
本发明提供的模拟方法通过求解离散格子玻尔兹曼方程来模拟波场传播过程，避免了波 动方程对波场模拟的限制，可以提高复杂介质中传播的声波波场参数的计算精度；同时，采 用本发明提供的新方法，不需要知道介质的品质因子，只需要已知流体介质的运动粘度，即 可模拟流体介质中粘滞声波的传播过程；此外，利用本发明公布的LBM方法，不仅能够准确 地模拟均匀介质中的粘滞声波波场传播现象，而且针对水平/倾斜层状介质、修改的Marmousi 模型等不规则/复杂介质模型，也可以获得高精度的多分量的波场模拟结果，这极大地扩展了 LBM系列数值模拟方法在地震波正演及相关领域的应用范围。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术 描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明中记 载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根 据这些附图获得其他的附图。
图1是本发明一个实施例提供的一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟 方法的方法流程图。
图2是本发明一个实施例提供的一种D2Q9模型的离散速度示意图。
图3是本发明一个实施例提供的一种倾斜层状介质模型几何结构及物理参数示意图。
图4是本发明一个实施例提供的一种基于倾斜层状介质模型模拟得到的LBM与FDM的波 场快照对比图。
图5是本发明一个实施例提供的一种基于倾斜层状介质模型模拟得到的LBM与FDM的波 剖面对比图。
图6是本发明一个实施例提供的一种修改的Marmousi速度模型。
图7是本发明一个实施例提供的一种基于修改的Marmousi模型模拟得到的LBM与FDM 的波场快照对比图。
具体实施方式
本发明实施例提供一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法。
为了使本技术领域的人员更好地理解本发明中的技术方案，下面将结合本发明实施例中 的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅 是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人 员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
图1是本发明所述一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法一种实施 例的方法流程图。虽然本发明提供了如下述实施例或附图所示的方法操作步骤或装置结构， 但基于常规或者无需创造性的劳动在所述方法或装置中可以包括更多或者更少的操作步骤 或模块单元。在逻辑性上不存在必要因果关系的步骤或结构中，这些步骤的执行顺序或装置 的模块结构不限于本发明实施例或附图所示的执行顺序或模块结构。所述的方法或模块结构 的在实际中的装置或终端产品应用时，可以按照实施例或者附图所示的方法或模块结构进行 顺序执行或者并行执行(例如并行处理器或者多线程处理的环境、甚至包括分布式处理的实 施环境)。
本发明的思路是用介观动态模拟方法在介观层面上追踪粒子集合的运动状态，达到在宏 观层面上模拟复杂物理现象的目的。当理想流体中具有微小扰动时，会产生声波，而地震纵 波可用声波方程来近似，因此，该方法可用于地震波/声波的正演模拟。
格子玻尔兹曼方法(LBM)是本发明发明人选择的一种最适合解决地震波/声波的正演 模拟的介观动态模拟方法，它与常规的基于波动方程的FDM方法不同，LBM是一种基于纳 维-斯托克斯方程的数值模拟方法，它不依赖于经典的波动方程，不仅能够用于计算含有强间 断面、多相流体等复杂情形的地质模型中的波动传播问题，而且具有灵活多样的流-固耦合边 界处理方式，可用于孔洞介质中的波场模拟、岩石物理建模、地震资料解释、流体预测、地 震监测与岩土防护工程等多个工程应用领域。
具体的如图1所述，本发明提供的一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模 拟方法包括：
S1：根据待模拟区域的结构特征，建立所述待模拟区域的数学模型。
本发明一个实施例中，所述数学模型可以包括待模拟区域介质的声波速度模型、密度模 型、品质因子模型、流体粘度模型等物理参数模型。
进一步的，所述数学模型可以是二维的，也可以是三维的。
S2：建立基于单松弛时间格子玻尔兹曼方程的离散模型，并根据所述离散模型中的物理 参数构建各个声波波场参数的迭代计算格式。
采用不同的离散格子玻尔兹曼模型进行数值模拟时，对介质的物理特性的刻画精度是不 一样的，相应的计算量也不同；特别地，在针对波动现象的正演模拟过程中，当模拟过程中 采用的离散格子玻尔兹曼模型不同时，LBM对粒子在内部界面、流-固边界或者计算区域外 部边界处的运动状态的模拟是有区别的，本发明人通过实验发现，离散速度的个数越多，LBM 对波动现象的描述就越精细，相应的计算量也越大。
本申请一个实施例中，针对二维地质模型的地震波场模拟问题，通常采用D2Q9模型来 描述离散粒子的运动状态，具体的示意图如图2所示，其中编号为0的圆圈表示某一时刻离散 粒子所处的网格节点，与其相邻的8个网格节点分别用编号为1、2、3……8的圆圈表示，由 编号为0的中心位置的圆圈指向相邻8个网格节点的箭头代表离散粒子可能的运动速度，分别 用向量e1、e2、e3……e8表示。特别地，当粒子在下一时刻在原地不动时，其运动速度用e0表 示。也就是说，当用D2Q9模型进行波场模拟时，计算区域内每个节点上的离散粒子在某一 时刻，既可以向其8个邻点迁移，也可以在原地静止不动。
本申请其他实施例中，针对三维地质模型的地震波场模拟问题，可以采用D3Q19模型、 D3Q15模型、D3Q27模型等LBM模型来描述离散粒子的运动状态，粒子迁移的原理与二维情 况下类似，这里不再赘述。
在实际应用过程中，具体选用哪种离散模型来做数值模拟，需要综合考虑实际问题所涉 及的求解精度、计算量等因素，选择最合适的离散格子玻尔兹曼模型。
在本发明一个实施例中，通过对纳维-斯托克斯方程进行单松弛时间格子简化处理，得到 所述单松弛时间格子玻尔兹曼方程为：
Figure BDA0003291534740000051
式中，fi(x,t)表示x位置、t时刻、i方向的粒子数密度；
fi (eq)(x,t)表示x位置、t时刻、i方向的平衡态粒子数密度；
τ表示松弛时间；
ci表示沿着i方向的离散速度。
在本发明一个实施例中，所述LBM方程中的松弛时间τ与待模拟区域介质的品质因子Q 之间的映射关系(简称Q-τ映射模型)为：
Figure BDA0003291534740000052
式中，fm表示震源函数的主频；Δt*表示物理系统中的时间采样间隔；τ表示松弛时间。
在本发明一个实施例中，所述单松弛时间格子玻尔兹曼方程中的松弛时间τ与流体的运 动粘度υ之间的定量关系为：
τ＝3υ+0.5
在采用LBM进行粘滞声波波场正演模拟时，τ的取值范围通常在0.5到2.0之间且不包 含0.5，因为介质的运动粘度值不可能为零。考虑到地震波是一种特殊的粘滞声波，故本发 明的说明书或者实施例中提到的地震波或者声波，一般特指具有粘滞特性的声波(即粘滞声 波)。
本发明另一个实施例中，所述平衡态粒子数密度的计算方式可以包括：利用平衡态粒子 数分布函数，计算得到所述平衡态粒子数密度。本发明再一个实施例中，所述平衡态粒子数 分布函数的函数表达式可以包括：
Figure BDA0003291534740000053
优选的，在求解某些完全线性的物理问题时，平衡态粒子数分布函数也可以省略高阶项 而采用如下函数表达式：
Figure BDA0003291534740000054
在平衡态时，所述粒子震动速度和粒子数密度应满足以下关系式：
Figure BDA0003291534740000055
Figure BDA0003291534740000061
式中，ρ表示宏观流体密度；
cs表示声速；
ci表示i方向的离散速度；
u表示粒子震动速度；
wi表示i方向的权系数；
n表示离散粒子在单位时间内跳跃的网格个数，且该数值可以不是整数。
在本发明一个实施例中，所述迭代计算格式，指的是与离散模型对应的迭代计算的模版， 通过建立离散迭代格式，可以确定相关参数的取值，比如离散速度、权系数、松弛时间、流 体粘度等。
比如，D2Q9模型对应的离散迭代格式，采用的离散速度可取值为：
Figure BDA0003291534740000062
D3Q19模型对应的离散迭代格式，采用的离散速度可取值为：
Figure BDA0003291534740000063
D2Q9模型对应的离散迭代格式，权系数取值为：
Figure BDA0003291534740000064
D3Q19模型对应的离散迭代格式，权系数取值为：
Figure BDA0003291534740000065
S3：对所述离散模型中各个预设方向的波场参数赋初值，定义震源子波函数，并制定波 阻抗界面处理方法、流-固边界条件以及计算区域外部边界条件。
所述赋初值，就是设定迭代计算格式的初始值，所述初值的确定可以由实施人员根据实 际情形确定，比如宏观流体密度的初始值通常定义为1、宏观流体速度的初始值定义为0，那 么，与之相对应的各个方向的粒子密度函数的初始值定义为各个方向的权系数值。
在本发明的一个实施例中，所述震源函数添加在各个方向粒子数密度函数的迭代公式 中，进一步地，所述震源函数包括Ricker子波函数，其表达式为：
Figure BDA0003291534740000071
式中，fm为震源函数的主频；
在本发明的其他实施例中，所述震源函数也可以添加在某些特定方向的粒子数密度函数 的迭代公式中，或者直接添加在宏观流体密度或者宏观流体速度的各个分量上。
在本发明的其他实施例中，震源函数也可以定义为高斯函数、单个/多个周期的简谐波函 数以及多个不同幅度值的简谐波函数的叠加等。
在本发明的一个实施例中，所述波阻抗界面处理方法通过定义如下反射(R)与透射系 数(T)实现：
Figure BDA0003291534740000072
Figure BDA0003291534740000073
其中，ρ1和ρ2表示波阻抗界面两侧介质的体积密度；n1和n2表示波阻抗界面两侧离散 粒子在单位时间内跳跃的网格个数，该数值不局限于取整数值；
在本发明的一个实施例中，所述流体粒子遇到模型内部边界时的运动状态，不仅需要满 足反射与透射定律，当界面两侧的粒子速度不再是整数(即粒子每次跳跃的距离不再是整数 网格长度)时，则需要通过插值算法来求得相应的整数网格节点上的粒子数密度，以做进一 步迭代计算。
在本发明的一个实施例中，所述流-固边界条件特指流体粒子遇到流-固边界时的运动状 态：在实际模拟中，粒子可以发生标准反弹，也可以发生镜面反弹，或者部分粒子标准反弹、 部分粒子镜面反弹；至于具体选用何种流-固边界条件，可根据实际工程问题的具体情况选择 合适的边界条件来做模拟计算。
在本发明的一个实施例中，所述计算区域外部边界条件特指流体粒子运动到计算区域外 部边界时的运动状态：为了防止计算区域外部边界的反射波干扰内部区域的波场，可以添加 指数衰减型吸收边界或者完全匹配层(PML)类型的吸收边界来压制人工边界反射波。当然， 在进行某些特殊情况下的地震波场模拟时，也可以采用自由表面边界条件、刚性边界条件等 其他不具备吸收波动能量特性的外部边界条件。
S4：根据所述迭代计算格式，利用单松弛时间格子玻尔兹曼方程，迭代更新计算得到预 设采样时刻的声波波场参数的数值。
所述预设采样时刻，可以由实施人员根据实际情况自行选择，比如，本发明一个实施例 中，所述采样时刻之间的间隔确定为0.5毫秒。本发明其他实施例中，所述采样时刻之间的间 隔也可以选择1.0毫秒、0.75毫秒、0.25毫秒等其他的时间采样间隔。
所述声波波场参数至少可以包括压力、粒子震动速度的各个分量。
其中，压力和粒子数密度直接相关，粒子震动速度与宏观的流体流动速度和介质体积密 度直接相关。
在本发明的一个实施例中，可以采用LBM计算得到地质模型中任意观测点处的压力、粒 子震动速度的各个分量在所计算时间段内的波场幅度值，即单道地震记录；进一步地，将多 个观测点(比如地表一条直线上的各个观测点)处的单道地震记录按照某种特殊的规则排列 在一起，即可获得相应的地震剖面。需要特别指出的是，本发明中介绍的波场传播速度与流 体流动速度是完全不同概念，其中涉及的研究区域内的压力、粒子震动速度等地震波场参数 是在非平衡态或者扰动状态下的实时数值大小，因而波场传播过程也是一个动态演化过程。
本申请一个实例中，对二维倾斜层状介质模型进行地震波场数值模拟，所述二维倾斜层 状介质模型的示意图如图2所示，其中，介质模型的尺寸为801m×801m，倾斜界面与左右 两侧边界的交点坐标分别为(0m，300m)和(801m，500m)。
所述二维倾斜层状介质模型中，倾斜界面上下两侧介质的纵波速度分别为1155m/s和 2310m/s，上下两侧介质的品质因子分别为13和64，对应的LBM松弛时间分别为1.0和0.6。
本实例中，LBM数值模拟采用的时间采样间隔为0.5毫秒，沿着x轴和z轴方向的空间采样 间隔均为1.0米，主频为25Hz的Ricker子波震源函数施加在9个方向的粒子数密度函数上，震 源点坐标为(361m，401m)。
对应的，采用二维情况下的D2Q9模型以及，基于本发明公布的Q-τ模型分别采用LBM与 FDM方法进行地震波场正演计算，得到的波场快照如图4(a)、图4(b)、图4(c)、图4(d)、图4(e)、 图4(f)所示，所述六幅图中，图4(a)中LBM-P表示图4(a)对应的是LBM方法计算的压力的波场 模拟结果，图4(b)中FDM-P表示图4(b)对应的是FDM方法计算的压力的波场模拟结果，图4(c) 中LBM-Vx表示图4(c)对应的是LBM方法计算的x轴方向上的震动速度的波场模拟结果，图4(d) 中FDM-Vx表示图4(d)对应的是FDM方法计算的x轴方向上的震动速度的波场模拟结果，图4(e) 中LBM-Vz表示图4(e)对应的是LBM方法计算的z轴方向上的震动速度的波场模拟结果，图4(f) 中FDM-Vz表示图4(f)对应的是FDM方法计算的z轴方向上的震动速度的波场模拟结果。以上6 幅图中，Depth(m)表示深度，单位是米。Distance(m)表示距离，单位是米。NormalizedAmplitude 代表归一化的波场振幅值。
所示的LBM方法模拟的压力以及x轴与z轴方向上的震动速度分量对应的波场快照与FDM计算的相应波场变量的波场快照几乎一样，表明本申请所述方法在层状介质中模拟获得 了与传统的FDM方法一致的波场分布，证明了新方法(LBM)的有效性和可行性。
进一步地，分别从图4的各个波场快照中抽出了600米深度、500米距离位置的波剖面叠 加在一起放在图5相应的子图中进行比较分析，发现两种数值模拟方法(LBM和FDM)计算 得到的波剖面(虚线和实线)几乎重合在一起，且通过计算图5中各个子图所示的波剖面之 间的相对误差，发现其最大的相对误差为1.14％，这说明在均匀介质模型、两层介质模型等 相对简单的介质中，LBM计算得到的波场与FDM求得的结果是非常吻合的。
本申请另一个实施例中，对修改的Marmousi模型(见图6)进行了地震波场正演模拟测 试。在模拟过程中，采用513×767的离散计算网格，纵横向空间网格边长均为1米，时间离 散采样间隔为0.5毫秒，LBM采用的松弛时间(τ)为0.6，相应地，根据本发明提供的Q–τ映 射模型推算出FDM中的品质因子(Q)为32。在修改的Marmousi模型的正演数值模拟计算中， 分别采用了D2Q9模型的离散迭代计算格式以及经典的二维FDM离散计算方法，以做对比验 证。地震波正演模拟时，震源均采用主频为50Hz的雷克子波震源，激发点坐标为(201m， 384m)，模拟得到的波场快照见图7。
图6显示的是修改的Marmousi速度模型，其中横轴Distance(m)代表距离，纵轴Depth(m) 表示深度，色标Velocity(m/s)代表速度值，图中每个像素点显示的不同颜色代表不同的速度 值。
图7中各个子图显示的是修改的Marmousi速度模型计算得到的LBM与FDM的各个物理 量的波场快照。所述六幅图中，图7(a)中LBM-P表示图7(a)对应的是LBM方法计算的压力的 波场模拟结果，图7(b)中FDM-P表示图7(b)对应的是FDM方法计算的压力的波场模拟结果， 图7(c)中LBM-Vx表示图7(c)对应的是LBM方法计算的x轴方向上的震动速度的波场模拟结 果，图7(d)中FDM-Vx表示图7(d)对应的是FDM方法计算的x轴方向上的震动速度的波场模拟 结果，图7(e)中LBM-Vz表示图7(e)对应的是LBM方法计算的z轴方向上的震动速度的波场模 拟结果，图7(f)中FDM-Vz表示图7(f)对应的是FDM方法计算的z轴方向上的震动速度的波场模 拟结果。以上6幅图中，Depth(m)表示深度，单位是米。Distance(m)表示距离，单位是米。 NormalizedAmplitude代表归一化的波场振幅值。
在本实施例中，LBM方法计算得到的波场快照与FDM的参考解在整体形态上是一致的， 但是细节上却有一些差异。在Marmousi模型中进行模拟测试时，两种数值模拟方法求得的两 组波场快照之间之所以会有一些差异，是因为两种方法的理论基础是不一样的：LBM是由纳 维-斯托克斯方程演化而来的一种新颖的波场正演模拟方法，它不依赖于传统的波动方程也不 受其连续性假设条件等限制，因而能够对地震波在界面处的传播过程做更加细致地刻画；而FDM是一种直接求解波动方程的正演模拟方法，其逃不出波动方程的限制条件。
总之，通过上述两个实施例的测试与分析，充分说明新方法(LBM)能够适应均匀介质 以及复杂的非均匀介质中的波场数值模拟应用需求，且LBM方法对地震波在复杂介质中的传 播过程描述得更加细致。此外，与传统的正演模拟方法不同，本发明提供的LBM方法在进行 粘滞声波数值模拟时，不需要知道介质的品质因子(通常情况下，品质因子是一个不太好提 取的介质物理参数)，LBM采用介质的运动粘度代替品质因子，再配合常规的地震物理模型 参数，即可进行地震波场的数值模拟并得到可靠的地震记录，这也是LBM方法相比较传统的 FDM等方法的一个优势。
利用上述各实施例提供的一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法 的实施方式，可以通过追踪微观粒子的运动及其相互作用，实现模拟宏观的地震波动的演化 过程，并迭代更新计算出每个预设采样时刻所述空间网格中各网格点的地震波场参数数值。 由于所述微观粒子的运动及相互作用情况不受波动方程限制，因此根据离散格子玻尔兹曼方 程模拟出的地震波动更接近实际的地震波动，计算出的地震波场参数数值的精度也更高，因 此，可以有效提高地震波场数值模拟的准确度。
尽管本发明内容中提到不同的地震波场数值模拟处理方式，从建立待模拟区域的几何模 型、在所述几何模型中建立空间网格、建立待模拟区域的离散模型、建立离散迭代格式到迭 代更新计算得到预设采样时刻所述空间网格中各网格点的地震波场参数数值的各种时序方 式、数据获取/处理/输出方式等的描述，但是，本发明并不局限于必须是行业标准或实施例 所描述的情况等，某些行业标准或者使用自定义方式或实施例描述的实施基础上略加修改后 的实施方案也可以实现上述实施例相同、等同或相近、或变形后可预料的实施效果。应用这 些修改或变形后的数据获取、处理、输出、判断方式等的实施例，仍然可以属于本发明的可 选实施方案范围之内。
虽然本发明提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的 手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序 中的一种方式，不代表唯一的执行顺序。在实际中的装置或客户端产品执行时，可以按照实 施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境， 甚至为分布式数据处理环境)。术语“包括”、“包含”或者其任何其他变体意在涵盖非排 他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而 且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有 的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备 中还存在另外的相同或等同要素。

Claims (5)

1.一种基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法，其特征在于，所述方法包括：
S1根据待模拟区域的结构特征，建立所述待模拟区域的数学模型；
S2建立基于单松弛时间格子玻尔兹曼方程的离散模型，并根据所述离散模型中的物理参数构建各个声波波场参数的迭代计算格式；
S3对所述离散模型中各个预设方向的波场参数赋初值，定义震源函数，并制定波阻抗界面处理方法、流-固边界条件以及计算区域外部边界条件；
S4根据所述迭代计算格式，利用单松弛时间格子玻尔兹曼方程，迭代更新计算得到预设采样时刻的声波波场参数的数值。
2.如权利要求1所述的方法，其特征在于，所述数学模型可以包括待模拟区域介质的声波速度模型、密度模型、品质因子模型等物理参数模型。
3.如权利要求1所述的方法，其特征在于，所述离散模型为DdQq模型，所述DdQq模型表示d维空间q个离散速度的离散格子玻尔兹曼模型。
4.如权利要求1所述的方法，其特征在于，所述单松弛时间格子玻尔兹曼方程为：
Figure FDA0003291534730000011
式中，fi(x,t)表示x位置、t时刻、i方向的粒子数密度；
fi (eq)(x,t)表示x位置、t时刻、i方向的平衡态粒子数密度；
τ表示松弛时间；ci表示i方向的离散速度。
5.如权利要求1所述的方法，其特征在于，所述震源函数包括Ricker子波函数。
CN202111166518.9A 2021-09-30 2021-09-30 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法 Pending CN114139335A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111166518.9A CN114139335A (zh) 2021-09-30 2021-09-30 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111166518.9A CN114139335A (zh) 2021-09-30 2021-09-30 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法

Publications (1)

Publication Number Publication Date
CN114139335A true CN114139335A (zh) 2022-03-04

Family

ID=80394052

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111166518.9A Pending CN114139335A (zh) 2021-09-30 2021-09-30 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法

Country Status (1)

Country Link
CN (1) CN114139335A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115345082A (zh) * 2022-07-06 2022-11-15 中山大学 一种面向冲击系统的二维九速离散玻尔兹曼方法及装置
CN116559942A (zh) * 2023-04-07 2023-08-08 中国科学院地质与地球物理研究所 一种流-固耦合的波场模拟方法、装置、设备及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1696356A1 (en) * 2005-02-24 2006-08-30 Siemens Aktiengesellschaft Flow acoustic simulation with the Lattice-Boltzmann method
KR20150091592A (ko) * 2014-02-03 2015-08-12 한국과학기술연구원 래티스 볼츠만 이론을 이용한 유체 유동 시뮬레이션 방법 및 이를 실현하기 위한 기록 매체
CN106021828A (zh) * 2016-07-15 2016-10-12 华中科技大学 一种基于格子-玻尔兹曼模型的流体模拟方法
CN106709191A (zh) * 2016-12-29 2017-05-24 中国石油大学(北京) 一种地震波场数值模拟方法及装置
CN108038304A (zh) * 2017-12-08 2018-05-15 西安交通大学 一种利用时间局部性的格子玻尔兹曼方法并行加速方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1696356A1 (en) * 2005-02-24 2006-08-30 Siemens Aktiengesellschaft Flow acoustic simulation with the Lattice-Boltzmann method
KR20150091592A (ko) * 2014-02-03 2015-08-12 한국과학기술연구원 래티스 볼츠만 이론을 이용한 유체 유동 시뮬레이션 방법 및 이를 실현하기 위한 기록 매체
CN106021828A (zh) * 2016-07-15 2016-10-12 华中科技大学 一种基于格子-玻尔兹曼模型的流体模拟方法
CN106709191A (zh) * 2016-12-29 2017-05-24 中国石油大学(北京) 一种地震波场数值模拟方法及装置
CN108038304A (zh) * 2017-12-08 2018-05-15 西安交通大学 一种利用时间局部性的格子玻尔兹曼方法并行加速方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
MUMING XIA: "Modelling viscoacoustic wave propagation with the lattice Boltzmann method", 《SCIENTIFIC REPORTS》, 31 August 2017 (2017-08-31) *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115345082A (zh) * 2022-07-06 2022-11-15 中山大学 一种面向冲击系统的二维九速离散玻尔兹曼方法及装置
CN115345082B (zh) * 2022-07-06 2023-05-16 中山大学 一种面向冲击系统的二维九速离散玻尔兹曼方法及装置
CN116559942A (zh) * 2023-04-07 2023-08-08 中国科学院地质与地球物理研究所 一种流-固耦合的波场模拟方法、装置、设备及存储介质

Similar Documents

Publication Publication Date Title
De La Puente et al. Discontinuous Galerkin methods for wave propagation in poroelastic media
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
Ladopoulos Elastodynamics for Non-linear Seismic Wave Motion in Real-Time Expert Seismology
US10495768B2 (en) Method of operating a data-processing system for the simulation of the acoustic wave propagation in the transversely isotropic media comprising an hydrocarbon reservoir
CN109725345B (zh) 一种初至波正演模拟方法及装置
CN109188519B (zh) 一种极坐标下的弹性波纵横波速度反演系统及方法
MXPA05006835A (es) Metodos para determinar los parametros de yacimiento y pozo de sondeo utilizando tomografia de volumen de fresnel.
Antonietti et al. Numerical modeling of seismic waves by discontinuous spectral element methods
CN114139335A (zh) 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法
Boehm et al. A semismooth Newton-CG method for constrained parameter identification in seismic tomography
Sjögreen et al. Source estimation by full wave form inversion
CN114114403B (zh) 一种基于分数阶拉氏算子的各向异性衰减介质模拟方法
Muratov et al. Grid-characteristic method on unstructured tetrahedral meshes
CN108680968B (zh) 复杂构造区地震勘探数据采集观测系统评价方法及装置
Gao et al. Parallel 3-D simulation of seismic wave propagation in heterogeneous anisotropic media: a grid method approach
Ma et al. Topography-dependent eikonal traveltime tomography for upper crustal structure beneath an irregular surface
CN113671570B (zh) 一种地震面波走时和重力异常联合反演方法与系统
KR20160012922A (ko) 직접 파형 역산의 반복 적용을 이용한 탄성파 영상화 장치 및 방법
CN111948708A (zh) 一种浸入边界起伏地表地震波场正演模拟方法
CN109738944B (zh) 基于广角反射的地震采集参数确定方法及装置
CN111505714A (zh) 基于岩石物理约束的弹性波直接包络反演方法
KR101614138B1 (ko) 모델링 영역 분해를 이용한 3차원 모멘트 텐서 역산방법
Sompotan et al. Comparing models GRM, refraction tomography and neural network to analyze shallow landslide
Shin et al. Laplace-domain full waveform inversion using irregular finite elements for complex foothill environments
Tadeu et al. 2.5 D elastic wave propagation in non-homogeneous media coupling the BEM and MLPG methods

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination