CN107798156A - 一种频率域2.5维粘弹性波数值模拟方法及装置 - Google Patents

一种频率域2.5维粘弹性波数值模拟方法及装置 Download PDF

Info

Publication number
CN107798156A
CN107798156A CN201610838492.0A CN201610838492A CN107798156A CN 107798156 A CN107798156 A CN 107798156A CN 201610838492 A CN201610838492 A CN 201610838492A CN 107798156 A CN107798156 A CN 107798156A
Authority
CN
China
Prior art keywords
wave
frequency domain
viscoelastic
equation
ties
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.)
Granted
Application number
CN201610838492.0A
Other languages
English (en)
Other versions
CN107798156B (zh
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.)
Individual
Original Assignee
Individual
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 Individual filed Critical Individual
Priority to CN201610838492.0A priority Critical patent/CN107798156B/zh
Publication of CN107798156A publication Critical patent/CN107798156A/zh
Application granted granted Critical
Publication of CN107798156B publication Critical patent/CN107798156B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Complex Calculations (AREA)

Abstract

本发明实施例提供一种频率域2.5维粘弹性波数值模拟方法及装置，所述方法包括：建立频率域2.5维粘弹性波波动方程；确定所述频率域2.5维粘弹性波动方程的等效积分弱形式；用双线性三角形单元离散得到有限元控制方程；加载频率域3维点源；加载刚度弱化吸收边界条件；确定波数采样策略；所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率‑波数域波场；采用复化辛普森数值积分法近似空间反傅里叶变换，得到频率域3维粘弹性波波场。本发明为2.5维粘弹性波多尺度全波形反演提供了一种准确的正演模拟方法，该方法能适应任意起伏地表复杂介质，其吸收边界条件的吸收效果好，计算效率高。

Description

一种频率域2.5维粘弹性波数值模拟方法及装置
技术领域
本发明涉及勘探地球物理领域，具体的说涉及一种频率域2.5维粘弹性波数值模拟方法及装置。
背景技术
全波形反演被认为是解决地震波成像问题的终极方法，其核心之一就是地震波正演模拟。通常由于2维地震波正演模拟假设震源为线源，其模拟记录与实际现场数据在动力学信息上有较大差异，所以2维正演模拟不能直接用于对现场数据的全波形反演中。3维正演模拟虽然更接近于实际现场数据，但目前计算量巨大。2.5维方法就是在2维介质中加载点源，从而获得3维地震波场，同时计算效率相比3维数值模拟大大提高。2.5正演模拟能直接应用于对一个剖面观测数据的全波形反演。
近几十年来，地震波数值模拟方法飞速发展。常用的数值解法有：有限差分法(FDM)，有限单元法(FEM)，谱元法(SEM)和间断有限元法(DG-FEM)。传统的有限差分法具有计算效率高，容易编程实现的特点，但在自由表面条件下需要特殊处理，不易准确模拟面波。此外，常规有限差分法在起伏地表条件下，会出现阶梯效应，在地表附近需要极密网格才能消除此效应，但这样会大大降低计算效率。有限元法能适应任意起伏地表复杂构造情形，且天然满足自由表面条件，能精确模拟面波。然而常规有限元法在低阶近似时数值频散较为严重。谱元法则是一种改进的高阶有限元法，在相同空间离散条件下，一般比有限元法更精确，但谱元法通常采用四边形网格(2维情形)或六面体网格(3维情形)，难以准确拟合地表尖灭等特殊复杂地形。间断有限元法在近十年来得到广泛发展，其主要优势是易于实现单元插值函数阶数的变化(P型自适应)，但由于单元边界上节点的重复，间断有限元法往往使得波场矢量的维数剧增，增加了较大的计算量。
实际观测数据与2维正演模拟数据在振幅和相位上有较大差异，而3维正演模拟虽更接近观测数据，但计算量巨大。为了得到接近实际的动力学信息，同时避免3维正演的计算量，2.5维正演模拟技术应运而生。2.5维正演模拟技术主要经历了两个发展阶段。第一阶段是基于声波介质的2.5维数值模拟。这期间，Bleistein(1986)定义了2.5维问题的三要素：i)物性参数仅在2维平面内变化，而在平面外的坐标系保持不变；ii)震源具有球对称性；iii)震源和检波点位于垂直于走向的平面内。Song和Williamson(1995)证明了傅里叶变换法解决2.5维问题的有效性。Zhou和Greenhalgh(1998)用有限元法计算了频率域2.5维的格林函数。由于声波模拟没能考虑界面上的转换波，2.5维问题逐渐过渡到第二个阶段-基于弹性波的2.5维数值模拟。Zhou(2012)基于一般各向异性弹性介质，采用谱元法，计算了频率域的格林函数。Liu(2011)实现了基于双孔隙孔弹性介质的2.5维数值模拟。
由于弹性波动方程形式简洁，故通常近似地球为弹性固体。但这往往与油藏工程和测井观测情况不一致，说明实际地下介质更应该是非完全弹性体。因此，在数值模拟时有必要考虑波场的物理衰减，即粘弹性波正演模拟。目前，基于粘弹性介质的正演模拟已经在2维或3维得到发展，但是针对勘探地震学的2.5维粘弹性波模拟仍未有文献记载。
发明内容
本发明实施提供了一种频率域2.5维粘弹性波数值模拟方法及装置，所提出的2.5维有限元解法可以高效准确计算2维任意起伏地表复杂构造介质中的3维波场，使用的粘弹性体模型能够近似地震频带(10-100Hz)的波场衰减特性，加载的刚度弱化(SRM)吸收边界条件可以有效压制边界虚假反射。
一方面，本发明实施例提供了一种频率域2.5维粘弹性波数值模拟方法，所述频率域2.5维粘弹性波数值模拟方法，包括：
建立频率域2.5维粘弹性波波动方程；
确定所述频率域2.5维粘弹性波动方程的等效积分弱形式；
用双线性三角形单元离散得到有限元控制方程；
加载频率域3维点源；
在频率-波数域加载刚度弱化(SRM)吸收边界条件；
确定波数采样策略；
所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率-波数域波场；
采用复化辛普森数值积分法近似空间反傅里叶变换，得到频率域3维粘弹性波波场；
可选的，在本发明一实施例中所述建立频率域2.5维粘弹性波波动方程，包括：采用Kelvin-Voigt流变模型表征粘弹性介质，将对应本构方程变换到频率域，然后代入频率域运动平衡方程中，得到频率域3维粘弹性波动方程。再对频率域3维粘弹性波动方程进行y方向的空间傅里叶变换，从而得到频率域2.5维粘弹性波波动方程。
可选的，在本发明一实施例中所述确定频率域2.5维粘弹性波动方程的等效积分弱形式，包括：基于伽辽金的加权余量法，对频率域2.5维粘弹性波动方程的余量关于形函数作加权积分，再根据分部积分原理，得到频率域2.5维粘弹性波动方程的等效积分弱形式。
另一方面，本发明实施例提供了一种频率域2.5维粘弹性波数值模拟的装置，所述频率域2.5维粘弹性波数值模拟的装置，包括：
有限元控制方程单元，用于得到大型稀疏线性方程组；
点源单元，用于在2维介质中加载点源；
吸收边界单元，用于加载SRM吸收边界条件；
波数采样单元，用于确定离散的波数；
稀疏方程组求解单元，用于采用直接解法求解器PARDISO求解稀疏线性方程组得到频率-波数域波场；
空间反傅里叶变换单元，用于得到频域域3维粘弹性波波场。
可选的，在本发明一实施例中所述点源单元，包括：用频率域雷克子波和3维空间震源的方向矢量表示3维方向力源。
可选的，在本发明一实施例中所述波数采样单元，包括：在有效波数区域内，根据最大波长，最大偏移距和波数区间的有效分段数确定波数采样的个数，然后等间隔采样。
可选的，在本发明一实施例中所述稀疏方程组求解单元，包括：利用PARDISO直接解法求解器，多线程并行技术，快速求解大型稀疏线性方程组。
可选的，在本发明一实施例中所述空间反傅里叶变换单元，包括：在求得空间某个检波点所有波数情形下的波场后，对这些波场沿y方向进行空间反傅里叶变换，用复化辛普森数值积分法近似该变换，最终得到该检波点处频率域3维波场值。
上述技术方案具有如下有益效果：因为建立频率域2.5维粘弹性波波动方程；确定所述频率域2.5维粘弹性波动方程的等效积分弱形式；用双线性三角形单元离散得到有限元控制方程；加载频率域3维点源；在频率-波数域加载SRM吸收边界条件；确定波数采样策略；所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率-波数域波场；采用复化辛普森数值积分法近似空间反傅里叶变换，得到频率域3维粘弹性波波场，所以达到了如下的技术效果：(1)能准确求算2维任意起伏地表复杂构造粘弹性介质中3维粘弹性波波场；(2)刚度弱化吸收边界条件加载简单，吸收效果好；(3)使用直接解法求解器PARDISO，多线程并行求解稀疏方程组，计算效率高，且不受系数矩阵的正定性影响。
附图说明
为了更清楚的说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单的介绍。显而易见，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供了一种频率域2.5维粘弹性波数值模拟的方法流程图；
图2为本发明实施例提供了一种频率域2.5维粘弹性波数值模拟的装置结构示意图；
图3为点震源矢量激发方向示意图；
图4(a)为本发明实施例加载SRM吸收边界条件后x方向激发，x分量虚部的频率域波场结果；
图4(b)为本发明实施例加载卷积完美匹配层(CPML)边界条件后x方向激发，x分量虚部的频率域波场结果；
图4(c)为本发明实施例SRM和CPML在相同厚度条件下，接收排列上x分量虚部波场值的对比图；
图5为本发明实施例用复化辛普森数值积分求得的2.5维数值解与3维解析解的对比图；
图6为本发明实施例起伏地表复杂粘弹性介质模型示意图；
图7为本发明实施例应用起伏地表复杂粘弹性介质模型的2.5维数值模拟结果图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅仅是发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，为本发明实施例提供了一种频率域2.5维粘弹性波数值模拟的方法流程图，所述频率域2.5维粘弹性波数值模拟的方法，包括：
101、建立频率域2.5维粘弹性波波动方程；
102、确定所述频率域2.5维粘弹性波动方程的等效积分弱形式；
103、用双线性三角形单元离散得到有限元控制方程；
104、加载频率域3维点源；
105、在频率-波数域加载刚度弱化(Stiffness Reduction Method)吸收边界条件；
106、确定波数采样策略；
107、所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率-波数域波场；
108、采用复化辛普森数值积分法近似空间反傅里叶变换，得到频率域3维粘弹性波波场；
可选的，在本发明一实施例中所述建立频率域2.5维粘弹性波波动方程，包括：采用Kelvin-Voigt流变模型表征粘弹性介质，将对应本构方程变换到频率域，然后代入频率域运动平衡方程中，得到频率域3维粘弹性波动方程。再对频率域3维粘弹性波动方程进行y方向的空间傅里叶变换，从而得到频率域2.5维粘弹性波波动方程。
可选的，在本发明一实施例中所述确定频率域2.5维粘弹性波动方程的等效积分弱形式，包括：基于伽辽金的加权余量法，对频率域2.5维粘弹性波动方程的余量关于形函数作加权积分，再根据分部积分原理，得到频率域2.5维粘弹性波动方程的等效积分弱形式。
如图2所示，本发明实施例提供了一种频率域2.5维粘弹性波数值模拟的装置，所述频率域2.5维粘弹性波数值模拟的装置，包括：
获取有限元控制方程单元201，用于得到大型稀疏线性方程组；
点源单元202，用于在2维介质中加载点源；
吸收边界单元203，用于加载SRM吸收边界条件；
波数采样单元204，用于确定离散的波数；
稀疏方程组求解单元205，用于采用直接解法求解器PARDISO求解稀疏线性方程组得到频率-波数域波场；
空间反傅里叶变换单元206，用于得到频域域3维粘弹性波波场。
可选的，在本发明一实施例中所述点源单元，包括：用频率域雷克子波和3维空间震源的方向矢量表示3维方向力源。
可选的，在本发明一实施例中所述波数采样单元，包括：在有效波数区域内，根据最大波长，最大偏移距和波数区间的有效分段数确定波数采样的个数，然后等间隔采样。
可选的，在本发明一实施例中所述稀疏方程组求解单元，包括：利用PARDISO直接解法求解器，多线程并行技术，快速求解大型稀疏线性方程组。
可选的，在本发明一实施例中所述空间反傅里叶变换单元，包括：在求得空间某个检波点所有波数情形下的波场后，对这些波场沿y方向进行空间反傅里叶变换，用复化辛普森数值积分法近似该变换，最终得到该检波点处频率域3维波场值。
本发明实施例上述技术方案具有如下有益效果：因为建立频率域2.5维粘弹性波波动方程；确定所述频率域2.5维粘弹性波动方程的等效积分弱形式；用双线性三角形单元离散得到有限元控制方程；加载频率域3维点源；在频率-波数域加载SRM吸收边界条件；确定波数采样策略；所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率-波数域波场；采用复化辛普森数值积分法近似空间反傅里叶变换，得到频率域3维粘弹性波波场，所以达到了如下的技术效果：(1)能准确求算2维任意起伏地表复杂构造粘弹性介质中3维粘弹性波波场；(2)刚度弱化吸收边界条件加载简单，吸收效果好；(3)使用直接解法求解器PARDISO，多线程并行求解稀疏方程组，计算效率高，且不受系数矩阵的正定性影响。
以下举应用实例进行详细说明：
本发明应用实例的目的在于提供一种准确计算2维任意起伏地表复杂构造粘弹性介质中粘弹性波波场的方法。加载的SRM吸收边界条件，实施简单，吸收效果好。借用MKL(Math Kernel Library)的PARDISO直接解法求解器模块，多线程快速求解稀疏方程组。
第一步，建立频率域2.5维粘弹性波波动方程。
采用Kelvin-Voigt流变模型表征粘弹性介质，其频率域本构方程为：
其中：
σ为应力张量，ε为应变张量，D为弹性系数张量，D′为粘弹性系数张量，ω为圆频率，i为虚单位。字母的上横线表示关于时间的傅里叶变换(下文同理)。那么频率域3维粘弹性波动方程可写为：
其中ρ表示密度，u表示频率域的位移矢量，表示应力张量的散度，f为震源矢量。
关于y方向对方程(3)进行空间傅里叶变换，并用字母上波浪线表示(下文同理)。根据变换算子(表示关于y方向的偏导，ky表示y方向的波数)得：
其中表示2.5维偏导算子。式(4)即为频率域2.5维粘弹性波波动方程。
第二步，确定所述频率域2.5维粘弹性波动方程的等效积分弱形式。
方程(4)的余量形式为：
基于伽辽金方法，令位移矢量的变分为权函数，则解满足下面的积分形式：
将式(5)代入式(6)得，并使用格林公式得：
上式(7)即为频率域2.5维粘弹性波动方程的等效积分弱形式。
第三步，用双线性三角形单元离散得到有限元控制方程。
首先将建立好的速度模型用合适的双线性三角形单元离散，则某一单元上的位移矢量可近似为：
其中为第e个单元第i个节点上的位移矢量。Ni(i＝1，2，3)为形函数，I表示3×3的单位矩阵。
在各向同性介质中，根据Voigt法则，由弹性张量和粘弹性张量组合的张量可以写成如下矩阵形式：
其中λ和μ为弹性拉梅常数，λ′和μ′为粘性拉梅常数。在粘弹性介质中主要用P波和S波的品质因子(QP和QS)来表征波的衰减特性。设P波波速为VP，S波波速为VS，那么在Kelvin-Voigt介质中λ，μ，λ′和μ′可以用下式表示：
λ＝ρVP 2-2ρVS 2，μ＝ρVS 2 (10)
同样依据Voigt法则应力张量和应变张量可以写成：
其中L为空间微分算子矩阵，矩阵B＝LN。于是频率域2.5维的本构方程可以表示为：
再将式(8)，(13)和(14)代入(7)中即可得到有限元控制方程：
其中M为全局质量矩阵，K为全局刚度矩阵，对应的单元矩阵分别为：
F为震源矢量，对应的单元矢量为：
第四步，加载频率域3维点源。
频率域雷克子波可表示为：
其中fp为峰值频率。于是震源矢量的模可定义为：
其中x0表示震源的位置矢量。如图3所示，为在XOZ平面上的投影；α为和之间的夹角；θ为与x轴之间的夹角，则震源矢量的分量可以表示为：
所以式(18)可以写成：
其中：
第五步，加载SRM吸收边界条件。
加载SRM吸收边界条件分两个方面。一方面令阻尼矩阵为：
C＝CMM (23)
则有限元控制方程(15)可以写为：
(23)式中的CM为阻尼系数，由下式估算：
CM(x)＝CMmacX(x)3 (25)
其中X(x)在SRM吸收层内逐渐从0增大到1，CM max＝ωλmax，λmax为入射波最大波长。
另一方面，需要使介质刚度在SRM吸收层内减小，即杨氏模量满足：
其中：
α(x)＝αmaxX(x)3
如图4所示，图4(a)为一均匀介质中加载SRM吸收边界条件后x方向激发，x分量虚部的频率域波场结果；图4(b)为相同均匀介质中加载卷积完美匹配层(CPML)边界条件后x方向激发，x分量虚部的频率域波场结果。对比图4(a)和(b)，定性说明了SRM吸收效果类似CPML。再看图4(c)，为SRM和CPML在相同厚度条件下，接收排列上x分量虚部波场值的对比图，从图中可以定量看出本实施例加载的SRM边界吸收效果良好。
第六步，确定波数采样策略。
频率域2.5维数值解是通过逐个波数求解2维数值解，然后通过空间域反傅里叶变换将这些波数域得2维数值解变换到空间域，可以表示为：
由于波数域存在一个截断波数当或波场会迅速衰减为0。这里的截断波数为：
于是(27)式变为：
由于震源和检波点的关系，波场具有对称性，因此上式可简化为：
波数采样策略的选取就是关于波数域的离散问题，对2.5维数值计算的效率有重要的影响。本实施例采用的波数采样方案为等间隔均匀离散，表示为：
其中Nky表示波数采样个数，rmax表示最大偏移距，λmax为最大波长。为波数采样的有效间隔数，一般大于等于10。
第七步，所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率-波数域波场。
在波数域离散后，固定一个频率，根据有限元控制方程(24)，借用Intel MKL(MathKernel Library)的PARDISO直接解法求解器模块，逐个波数多线程快速求解稀疏方程组。PARDISO的使用方法说明参见Intel MKL用户手册。
第八步，采用复化辛普森数值积分法近似空间反傅里叶变换，得到频率域3维粘弹性波波场。
对于单个频率，在得到一系列波数的2维数值解后，根据(29)式，利用复化辛普森数值积分法近似空间反傅里叶变换，即可得到频率域3维粘弹性波波场。如图5所示，为一均匀粘弹性介质的频率域2.5维数值解(红色点画线)与3维解析解(黑色实线)的对比。从图中可以看出，本实施例计算的3维粘弹性波波场是准确的。如图6所示，为起伏地表复杂粘弹性介质模型示意图，可以看出三角单元网格能很好的拟合起伏地表；如图7所示，为该粘弹性介质模型的2.5维数值模拟结果图。
本发明应用实例与现有技术对比，具有如下的有益效果：(1)能准确求算2维任意起伏地表复杂构造粘弹性介质中3维粘弹性波波场；(2)刚度弱化吸收边界条件加载简单，吸收效果好；(3)使用直接解法求解器PARDISO，多线程并行求解稀疏方程组，计算效率高，且不受系数矩阵的正定性影响。
本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块(illustrative logical block)，单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性(interchangeability)，上述的各种说明性部件(illustrative components)，单元和步骤已经通用的描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。
本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动硬盘、CD-ROM或本领域中其它任意形式的存储媒介中。示例性的，存储媒介可以与处理器连接，以使得处理器可以从存储媒介中读取信息，并可以向存储媒介存写信息。可选的，存储媒介还可以集成到处理器中。处理器和存储媒介可以设置与ASIC中，ASIC可以设置于用户终端中。可选的，处理器和存储媒介也可以设置于用户终端中的不同的部件中。
在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储于电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当的定义为电脑可读媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电缆、双绞线、数字用户线(DSL)或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片(disk)和磁盘(disc)包括压缩磁盘、镭射盘、光盘、DVD、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光学复制数据。上述的组合也可以包含在电脑可读媒介中。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种频率域2.5维粘弹性波数值模拟的方法，其特征在于，所述的频率域2.5维粘弹性波数值模拟的方法包括：
建立频率域2.5维粘弹性波波动方程；
确定所述频率域2.5维粘弹性波动方程的等效积分弱形式；
用双线性三角形单元离散得到有限元控制方程；
加载频率域3维点源；
在频率-波数域加载刚度弱化(SRM)吸收边界条件；
确定波数采样策略；
所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率-波数域波场；
采用复化辛普森数值积分法近似空间反傅里叶变换，得到频率域3维粘弹性波波场。
2.如权利要求1所述频率域2.5维粘弹性波数值模拟方法，其特征在于，所述建立频率域2.5维粘弹性波波动方程，包括：
采用Kelvin-Voigt流变模型表征粘弹性介质，将对应本构方程变换到频率域，然后代入频率域运动平衡方程中，得到频域3维粘弹性波动方程。再对频率域3维粘弹性波动方程进行y方向的空间傅里叶变换，从而得到频率域2.5维粘弹性波波动方程。
3.如权利要求1所述频率域2.5维粘弹性波数值模拟方法，其特征在于，所述确定频率域2.5维粘弹性波动方程的等效积分弱形式，包括：
基于伽辽金的加权余量法，对频率域2.5维粘弹性波动方程的余量关于形函数作加权积分，再根据分部积分原理，得到频率域2.5维粘弹性波动方程的等效积分弱形式。
4.如权利要求1所述频率域2.5维粘弹性波数值模拟方法，其特征在于，所述用双线性三角形单元离散得到有限元控制方程，包括：
非结构化和结构化的混合网格剖分，将位移场插值函数表示为形函数的线性组合，然后代入到等效积分弱形式，从而得到有限元控制方程。
5.如权利要求1所述频率域2.5维粘弹性波数值模拟方法，其特征在于，所述加载频率域3维点源，包括：
用频率域雷克子波和3维空间震源的方向矢量表示3维方向力源。
6.如权利要求1所述频率域2.5维粘弹性波数值模拟方法，其特征在于，所述确定波数采样策略，包括：
在有效波数区域内，根据最大波长，最大偏移距和波数区间的有效分段数确定波数采样的个数，然后等间隔采样。
7.如权利要求1所述频率域2.5维粘弹性波数值模拟方法，其特征在于，所述有限元控制方程为大型稀疏线性方程组，用直接解法求解所述方程组，得到频率-波数域波场，包括：
利用PARDISO直接解法求解器，多线程并行技术，求解大型稀疏线性方程组。
8.如权利要求1所述频率域2.5维粘弹性波数值模拟方法，其特征在于，所述采用复化辛普森数值积分法近似空间反傅里叶变换得到频率域3维粘弹性波波场，包括：
在求得空间某个检波点所有波数情形下的波场后，对这些波场沿y方向进行空间反傅里叶变换，用复化辛普森数值积分法近似该变换，最终得到该检波点处频率域3维波场值。
CN201610838492.0A 2016-09-02 2016-09-02 一种频率域2.5维粘弹性波数值模拟方法及装置 Expired - Fee Related CN107798156B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610838492.0A CN107798156B (zh) 2016-09-02 2016-09-02 一种频率域2.5维粘弹性波数值模拟方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610838492.0A CN107798156B (zh) 2016-09-02 2016-09-02 一种频率域2.5维粘弹性波数值模拟方法及装置

Publications (2)

Publication Number Publication Date
CN107798156A true CN107798156A (zh) 2018-03-13
CN107798156B CN107798156B (zh) 2020-12-11

Family

ID=61530114

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610838492.0A Expired - Fee Related CN107798156B (zh) 2016-09-02 2016-09-02 一种频率域2.5维粘弹性波数值模拟方法及装置

Country Status (1)

Country Link
CN (1) CN107798156B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109460587A (zh) * 2018-10-22 2019-03-12 中国地震局地壳应力研究所 火山和地震粘弹性形变自动建模有限元计算方法
CN110263434A (zh) * 2019-06-20 2019-09-20 中国石油大学(华东) 一种基于多尺度混合有限元的流动单元数值模拟方法
CN110296327A (zh) * 2019-06-19 2019-10-01 常州大学 一种基于瞬变流频率响应分析的管道泄漏检测方法
CN110988993A (zh) * 2019-11-27 2020-04-10 清华大学 一种偏移成像方法、装置及电子设备
CN112906267A (zh) * 2021-02-07 2021-06-04 中国石油大学(北京) 频率域粘弹性波正演模拟方法、装置、设备及存储介质
CN112966412A (zh) * 2021-02-07 2021-06-15 中国石油大学(北京) 频率域2.5维地震波正演模拟方法、装置及设备
CN113221392A (zh) * 2021-01-26 2021-08-06 中国地震局工程力学研究所 一种非均匀粘声波在无限域内传播模型的构建方法
CN115828692A (zh) * 2022-12-05 2023-03-21 东莞理工学院 一种计算任意横截面结构超声导波频散特性的方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011048535A1 (en) * 2009-10-19 2011-04-28 Padia Bhadresh K Sustainable chemical process for reduction of nitro compounds (r-no2) or nitroso compounds (r-no) containing sulphonic or carboxylic group into corresponding amino compounds (r-nh2) with inherent recycle of all acidic streams generated in synthesis
EP2543313A1 (en) * 2010-10-13 2013-01-09 Kabushiki Kaisha Toshiba Magnetic resonance imaging apparatus and magnetic resonance imaging method
CN102914799A (zh) * 2012-10-12 2013-02-06 中国石油天然气股份有限公司 非等效体波场正演模拟方法及装置
CN103226197A (zh) * 2013-04-16 2013-07-31 哈尔滨工程大学 一种基于音色参数模型的水下目标回波分类方法
CN104598672A (zh) * 2014-12-30 2015-05-06 中国科学院地质与地球物理研究所 一种计算发射源姿态改变引起的电磁响应误差的方法
CN104749628A (zh) * 2015-03-30 2015-07-01 西安交通大学 一种基于弥散黏滞性波动方程的吸收边界反射方法
CN105093265A (zh) * 2014-05-09 2015-11-25 中国石油化工股份有限公司 一种模拟地震波在ti介质中传播规律的方法
CN105676280A (zh) * 2016-01-21 2016-06-15 中国矿业大学(北京) 基于旋转交错网格的双相介质地质数据获取方法和装置

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011048535A1 (en) * 2009-10-19 2011-04-28 Padia Bhadresh K Sustainable chemical process for reduction of nitro compounds (r-no2) or nitroso compounds (r-no) containing sulphonic or carboxylic group into corresponding amino compounds (r-nh2) with inherent recycle of all acidic streams generated in synthesis
EP2543313A1 (en) * 2010-10-13 2013-01-09 Kabushiki Kaisha Toshiba Magnetic resonance imaging apparatus and magnetic resonance imaging method
CN102914799A (zh) * 2012-10-12 2013-02-06 中国石油天然气股份有限公司 非等效体波场正演模拟方法及装置
CN103226197A (zh) * 2013-04-16 2013-07-31 哈尔滨工程大学 一种基于音色参数模型的水下目标回波分类方法
CN105093265A (zh) * 2014-05-09 2015-11-25 中国石油化工股份有限公司 一种模拟地震波在ti介质中传播规律的方法
CN104598672A (zh) * 2014-12-30 2015-05-06 中国科学院地质与地球物理研究所 一种计算发射源姿态改变引起的电磁响应误差的方法
CN104749628A (zh) * 2015-03-30 2015-07-01 西安交通大学 一种基于弥散黏滞性波动方程的吸收边界反射方法
CN105676280A (zh) * 2016-01-21 2016-06-15 中国矿业大学(北京) 基于旋转交错网格的双相介质地质数据获取方法和装置

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
J.Y.SHIH 等: "The effect of boundary conditions, model size and damping models in the finite element modelling of a moving load on a track/ground system", 《SOIL DYNAMICS AND EARTHQUAKE ENGINEERING》 *
NIE JIAN-XIN 等: "BISQ model based on a Kelvin-Voigt viscoelastic frame in a partially saturated porous medium", 《APPLIED GEOPHYSICS》 *
李芳 等: "2.5D傅氏变换法声波方程数值模拟及精度分析", 《石油地球物理勘探》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109460587A (zh) * 2018-10-22 2019-03-12 中国地震局地壳应力研究所 火山和地震粘弹性形变自动建模有限元计算方法
CN109460587B (zh) * 2018-10-22 2020-02-28 中国地震局地壳应力研究所 火山和地震粘弹性形变自动建模有限元计算方法
CN110296327A (zh) * 2019-06-19 2019-10-01 常州大学 一种基于瞬变流频率响应分析的管道泄漏检测方法
CN110263434A (zh) * 2019-06-20 2019-09-20 中国石油大学(华东) 一种基于多尺度混合有限元的流动单元数值模拟方法
CN110988993A (zh) * 2019-11-27 2020-04-10 清华大学 一种偏移成像方法、装置及电子设备
CN113221392A (zh) * 2021-01-26 2021-08-06 中国地震局工程力学研究所 一种非均匀粘声波在无限域内传播模型的构建方法
CN113221392B (zh) * 2021-01-26 2023-12-19 中国地震局工程力学研究所 一种非均匀粘声波在无限域内传播模型的构建方法
CN112906267A (zh) * 2021-02-07 2021-06-04 中国石油大学(北京) 频率域粘弹性波正演模拟方法、装置、设备及存储介质
CN112966412A (zh) * 2021-02-07 2021-06-15 中国石油大学(北京) 频率域2.5维地震波正演模拟方法、装置及设备
CN115828692A (zh) * 2022-12-05 2023-03-21 东莞理工学院 一种计算任意横截面结构超声导波频散特性的方法
CN115828692B (zh) * 2022-12-05 2023-11-17 东莞理工学院 一种计算任意横截面结构超声导波频散特性的方法

Also Published As

Publication number Publication date
CN107798156B (zh) 2020-12-11

Similar Documents

Publication Publication Date Title
CN107798156A (zh) 一种频率域2.5维粘弹性波数值模拟方法及装置
CN106842306B (zh) 一种全局优化的交错网格有限差分正演模拟方法和装置
Bonev et al. Discontinuous Galerkin scheme for the spherical shallow water equations with applications to tsunami modeling and prediction
CN111123354B (zh) 基于频变反射振幅衰减预测致密气层的方法及设备
Wang Sharp nonlinear stability criterion of viscous non-resistive MHD internal waves in 3D
Manteuffel et al. A root-node--based algebraic multigrid method
Yao et al. An effective absorbing layer for the boundary condition in acoustic seismic wave simulation
CN110879412A (zh) 地下横波速度反演方法、装置、计算设备及存储介质
Sousa et al. A finite difference method with meshless interpolation for incompressible flows in non-graded tree-based grids
CN112578450A (zh) 基于频散介质标量波方程有限差分的正演模拟方法及装置
Fang et al. Analysis of GPR Wave Propagation Using CUDA‐Implemented Conformal Symplectic Partitioned Runge‐Kutta Method
CN108897036A (zh) 一种地震数据处理方法及装置
CN116611309A (zh) 物理场仿真方法、系统、介质和电子设备
Ginting et al. Artificial viscosity technique: A Riemann-solver-free method for 2D urban flood modelling on complex topography
CN104122584A (zh) 根据地震数据确定方向性的方法及装置
Li et al. Rapid Inference of Reservoir Permeability From Inversion of Travel Time Data Under a Fast Marching Method Based Deep Learning Framework
Wu et al. Development of a Cell‐Centered Godunov‐Type Finite Volume Model for Shallow Water Flow Based on Unstructured Mesh
Fu A high-order velocity-based discontinuous Galerkin scheme for the shallow water equations: Local conservation, entropy stability, well-balanced property, and positivity preservation
CN104750954B (zh) 一种在复杂各向异性介质中模拟地震波的方法及装置
Markowskei et al. Measuring the effect of rounding the corners of scattering structures
CN106570201A (zh) 粘弹介质数值模拟方法及系统
CN109800508A (zh) 嵌岩桩桩端的空洞顶板厚度的计算方法及终端设备
Ortega et al. An adaptive finite point method for the shallow water equations
Wu et al. Trapezoid-grid finite difference frequency domain method for seismic wave simulation
CN107102353B (zh) 基于高阶差分方法的弹性波方程逆时偏移成像方法

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
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20201211

Termination date: 20210902