CN111257930A - 一种黏弹各向异性双相介质区域变网格求解算子 - Google Patents

一种黏弹各向异性双相介质区域变网格求解算子 Download PDF

Info

Publication number
CN111257930A
CN111257930A CN202010097267.2A CN202010097267A CN111257930A CN 111257930 A CN111257930 A CN 111257930A CN 202010097267 A CN202010097267 A CN 202010097267A CN 111257930 A CN111257930 A CN 111257930A
Authority
CN
China
Prior art keywords
anisotropic
viscoelastic
representing
equation
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.)
Granted
Application number
CN202010097267.2A
Other languages
English (en)
Other versions
CN111257930B (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.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
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 East China filed Critical China University of Petroleum East China
Priority to CN202010097267.2A priority Critical patent/CN111257930B/zh
Publication of CN111257930A publication Critical patent/CN111257930A/zh
Application granted granted Critical
Publication of CN111257930B publication Critical patent/CN111257930B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/58—Media-related
    • G01V2210/586—Anisotropic media
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/66—Subsurface modeling

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Vibration Prevention Devices (AREA)

Abstract

本发明公开了一种黏弹各向异性双相介质区域变网格求解算子，包括以下步骤：输入速度场及各种参数场，建立观测系统；选择目标区域进行局部网格加密；在全局区域采用黏声各向异性介质进行波场延拓；在目标区域采用黏弹各向异性双相介质进行波场延拓；在目标区域的边界处进行黏声各向异性和黏弹各向异性双相波场传递；应用人工边界条件对边界反射进行吸收；输出炮记录和波场快照。本发明黏弹各向异性双相介质区域变网格求解算子，在全局区域采用黏声各向异性介质的声波算子进行波场延拓，在目标区域采用黏弹各向异性双相介质的弹性波算子进行延拓，提高了地震波场的模拟准确度并提高了计算效率。

Description

一种黏弹各向异性双相介质区域变网格求解算子
技术领域
本发明属于石油勘探技术领域，具体涉及一种黏弹各向异性双相介质区域变网格求解算子。
背景技术
地震波波场正演模拟越来越受到地球物理学者们的重视，它的地位也显得越来越重要。其中，地球介质具有黏弹性、各向异性和双相性，目前针对地震波波场正演模拟大都基于以上其中一种或两种性质考虑。但是，这种基于一种或两种性质模拟得到的地震波数据与基于三种性质模拟出的多分量的地震波数据在流体识别与描述、裂缝分布估计、岩性估计、各向异性分析、亮点反射设别等方面都有所差距，基于三种性质模拟出的多分量的地震波数据更加有效。
虽然基于三种性质模拟出的多分量的地震波数据更加有效，但是全地下构造均采用黏弹各向异性双相介质弹性波正演模拟方法进行地震波波场模拟时，模拟计算量巨大，现有的计算机硬件难以承受。因此，本发明提出了一种黏弹各向异性双相介质区域变网格求解算子，既同时考虑了地球介质的黏弹性、各向异性和双相性，又只针对目标区域进行弹性纵横波场的模拟，从而降低模拟的计算量。
发明内容
本发明的目的是为克服上述现有技术的不足，提供一种黏弹各向异性双相介质区域变网格求解算子。
为实现上述目的，本发明采用如下技术方案：
一种黏弹各向异性双相介质区域变网格求解算子，包括以下步骤：
步骤1：输入速度场及各种参数场，建立观测系统；
步骤2：选择目标区域进行局部网格加密；
步骤3：在全局区域采用黏声各向异性介质进行波场延拓；
步骤4：在目标区域采用黏弹各向异性双相介质进行波场延拓；
步骤5：在目标区域的边界处进行黏声各向异性和黏弹各向异性双相波场传递；
步骤6：应用人工边界条件对边界反射进行吸收；
步骤7：输出炮记录和波场快照。
优选的，所述步骤2中，在全局区域进行粗网格剖分，在目标区域进行细网格剖分。
优选的，所述步骤2中，全局区域的粗网格尺寸为目标区域细网格尺寸的奇数倍。
优选的，所述步骤3中，在全局区域采用黏声各向异性波动方程(1)进行P波传播；
Figure BDA0002385881370000021
方程(1)中，p表示各向异性声波波场；q表示各向异性声波辅助波场；vpz表示对称轴方向的相速度；ε表示用以表征纵波各向异性强度的Thomson参数，δ表示用以表征动校正速度的Thomson参数；t表示时间；fx表示水平分量的震源；fz表示垂直分量的震源；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；τ表示松弛时间，τ由公式(2)求得：
τ＝τε/τσ-1 (2)
公式(2)中，τσ为应力松弛时间，τε为应变松弛时间；其中，τσ和τε与纵波品质因子Qp的关系如公式(3)、公式(4)所示；
Figure BDA0002385881370000031
Figure BDA0002385881370000032
公式(3)和公式(4)中，ω代表角频率。
优选的，所述步骤3中全局区域采用的黏声各向异性波动方程(1)由公式(5)递推求得；
Figure BDA0002385881370000033
公式(5)中，n表示差分计算次数；F和F-1分别为傅里叶变换与傅里叶反变换；Δt是时间步长；kx为水平分量的波数；kz为垂直分量的波数；
Figure BDA0002385881370000034
为水平分量的空间二阶偏导数；
Figure BDA0002385881370000035
为垂直分量的空间二阶偏导数；
Figure BDA0002385881370000036
和
Figure BDA0002385881370000037
由公式(6)～公式(9)求得：
Figure BDA0002385881370000038
Figure BDA0002385881370000039
Figure BDA00023858813700000310
Figure BDA00023858813700000311
公式(6)～公式(9)中，Δx表示水平分量的空间网格步长；Δz表示垂直分量的空间网格步长；c2,0和c2,m为二阶偏导数的有限差分系数；m为积分变量；M为差分阶数；i为水平分量的空间坐标计数，j为垂直分量的空间坐标计数。
优选的，所述步骤4中，在目标区域采用黏弹各向异性双相介质一阶速度应力方程(10)进行多分量波场延拓；
Figure BDA0002385881370000041
方程(10)中，ux为固相水平分量的质点速度；uz为固相垂直分量的质点速度；Ux为流相水平分量的质点速度；Uz为流相垂直分量的质点速度；τxx为水平分量的正应力；τzz为垂直分量的正应力；τxz表示切应力；S表示流相正应力；Q为固体与流体之间相耦合弹性系数；R为流相弹性系数；b表示耗散系数；t表示时间；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；L1表示体积模量的标准线性固体个数；L2表示剪切模量的标准线性固体个数；
Figure BDA0002385881370000042
表示胀缩模式的记忆变量；
Figure BDA0002385881370000043
表示体积模量的第l个标准线性固体的记忆变量；
Figure BDA0002385881370000044
表示剪切模量的第l个标准线性固体的记忆变量；
Figure BDA0002385881370000045
Figure BDA0002385881370000046
为各向异性粘弹介质的弹性系数；e11、e55为各向异性弹性系数；D1、D2和D3分别为三个中间变量；
方程(10)中，D1、D2和D3由公式(11)求得：
Figure BDA0002385881370000051
公式(11)中，ρ11表示单位体积中固体相对流体运动时固体部分总的等效质量，ρ22表示单位体积中流体相对固体运动时流体部分总的等效质量，ρ12表示单位体积中流体和固体之间的质量耦合系数；
方程(10)中，
Figure BDA0002385881370000052
和
Figure BDA0002385881370000053
由公式(12)求得；
Figure BDA0002385881370000054
公式(12)中，
Figure BDA0002385881370000055
表示体积模量的应力松弛时间常数；
Figure BDA0002385881370000056
表示剪切模量的应力松弛时间常数；
Figure BDA0002385881370000057
表示体积模量的应变松弛时间常数；
Figure BDA0002385881370000058
表示剪切模量的应变松弛时间常数；
方程(10)中，e11、e55由公式(13)求得：
Figure BDA0002385881370000059
公式(13)中，ε表示用以表征纵波各向异性强度的Thomson参数；ρ表示密度；vp表示纵波速度；vs表示横波速度；
方程(10)中，
Figure BDA00023858813700000510
由公式(14)求得：
Figure BDA0002385881370000061
公式(14)中，λ为与纵横波速度和密度有关的拉梅常数；μ为与横波速度和密度有关的拉梅常数；
Figure BDA0002385881370000062
为体积模量应力松弛时间常数与应变松弛时间常数的比值；
Figure BDA0002385881370000063
分别为剪切模量应力松弛时间常数与应变松弛时间常数的比值和；
Figure BDA0002385881370000064
和
Figure BDA0002385881370000065
由公式(15)求得；
Figure BDA0002385881370000066
优选的，所述步骤4中，目标区域采用的黏弹各向异性双相介质一阶速度应力方程(10)由公式(16)递推求得；
Figure BDA0002385881370000071
公式(16)中，Δt1为目标区域的时间采样间隔，Δt1与时间步长Δt之间的关系如公式(17)所示；
Δt＝(2k+1)Δt1 (17)
公式(17)中，k为任一正整数；
公式(16)中，Dx表示水平分量的空间一阶偏导数，Dz表示垂直分量的空间一阶偏导数，Dx和Dz由公式(18)求得：
Figure BDA0002385881370000072
公式(18)中，U为波场值，可表示ux、uz、Ux、Uz、τxx、τzz、τxz、
Figure BDA0002385881370000073
Figure BDA0002385881370000074
和
Figure BDA0002385881370000075
c1,m为一阶偏导数的有限差分系数；Δx1为目标区域水平分量的网格大小；Δz1为目标区域垂直分量的网格大小；Δx1和水平分量的空间网格步长Δx的关系及Δz1和垂直分量的空间网格步长Δz的关系如公式(19)所示；
Figure BDA0002385881370000076
优选的，所述步骤5中，在目标区域的左右边界采用公式(20)所示的边界条件进行声压P与应力τxx的波场传递；
-P＝τxx (20)
在目标区域的上下边界采用公式(21)进行声压P与应力τzz的波场传递；
-P＝τzz (21)。
优选的，所述步骤5中，目标区域边界处全局区域延拓的黏声各向异性介质声波波场与目标区域延拓的黏弹各向异性双相介质弹性波波场互相传递的步骤如下：
步骤51：将目标区域上下边界上处于粗网格点上的声压P传递给对应的细网格点上的τxx，粗网格点上的uz传递给对应细网格点上的uz；将目标区域左右边界上处于粗网格点上的声压P传递给对应的细网格点上的τzz，粗网格点上的ux传递给对应细网格点上的ux；
步骤52：利用插值公式计算目标区域边界上其他细网格点上的波场值；
步骤53：采用降阶方法求解过渡区域的波场值；
步骤54：在目标区域其他细网格点黏弹各向异性双相介质波动方程(16)进行波场延拓；
步骤55：将目标区域的过渡区域上位于细网格点上的波场值传递给对应的粗网格；
步骤56：将目标区域上下边界上位于粗网格位置上的细网格点波场值τxx传递给粗网格点波场值P，将目标区域左右边界上位于粗网格位置上的细网格点波场值τzz传递给粗网格点波场值P，位于粗网格位置上的细网格点波场值ux和uz传递给粗网格点波场值ux和uz。
优选的，所述步骤6中，在全局区域采用PML边界条件对波场的衰减进行吸收，对目标区域的边界采用CFS-PML边界条件进行波场延拓；
使用拉伸函数将速度-应力方程变换到拉伸坐标域；在频率域中，2D拉伸函数的具体表述形式：
Figure BDA0002385881370000091
Figure BDA0002385881370000092
式中，σx代表CFS-PML中的水平分量衰减系数；σz代表CFS-PML中的垂直分量衰减系数；εx为水平分量的拉伸函数；εz为垂直分量的拉伸函数；ω为角频率；I表示复数虚单位；
笛卡尔坐标系(x，z)与拉伸坐标系
Figure BDA0002385881370000093
的坐标变换方程如下：
Figure BDA0002385881370000094
将方程(24)代入方程(10)，得到拉伸坐标系下一阶速度-应力方程的表达式：
Figure BDA0002385881370000101
其中，
Figure BDA0002385881370000102
为拉伸坐标系下的固相水平分量的质点速度；
Figure BDA0002385881370000103
为拉伸坐标系下的固相垂直分量的质点速度；
Figure BDA0002385881370000104
为拉伸坐标系下的流相水平分量的质点速度；
Figure BDA0002385881370000105
为拉伸坐标系下的流相垂直分量的质点速度；
Figure BDA0002385881370000106
表示拉伸坐标系下水平分量的正应力；
Figure BDA0002385881370000107
表示拉伸坐标系下垂直分量的正应力；
Figure BDA0002385881370000108
表示拉伸坐标系下的切应力；
Figure BDA0002385881370000109
表示拉伸坐标系下的流相正应力；
Figure BDA00023858813700001010
表示拉伸坐标系下的胀缩模式的记忆变量；
Figure BDA00023858813700001011
表示拉伸坐标系下的体积模量的第l个标准线性固体的记忆变量；
Figure BDA00023858813700001012
表示拉伸坐标系下的剪切模量的第l个标准线性固体的记忆变量。
本发明的有益效果是：
本发明黏弹各向异性双相介质区域变网格求解算子，同时考虑了地下介质的黏弹性、各向异性和双相性，在全局区域采用黏声各向异性介质的声波算子进行波场延拓，在目标区域采用黏弹各向异性双相介质的弹性波算子进行延拓，在目标区域的边界上进行稳定的声波和弹性波场传递，并在目标区域采用细网格剖分，提高了地震波场的模拟准确度并提高了计算效率。
附图说明
构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。
图1是本发明黏弹各向异性双相介质区域变网格求解算子的结构示意图；
图2(a)为实施例中黏弹各向异性双相介质模型的双相介质参数；
图2(b)为实施例中黏弹各向异性双相介质模型的速度密度参数；
图2(c)为实施例中黏弹各向异性双相介质模型的品质因子参数；
图2(d)为实施例中黏弹各向异性双相介质模型的各向异性参数；
图3(a)为采用本发明模拟得到的地震波正演模拟炮记录的固相水平分量图；
图3(b)为采用本发明模拟得到的地震波正演模拟炮记录的固相垂直分量图；
图3(c)为采用本发明模拟得到的地震波正演模拟炮记录的流相水平分量图；
图3(d)为采用本发明模拟得到的地震波正演模拟炮记录的流相垂直分量图；
图4(a)为采用本发明模拟得到的波场快照的固相水平分量图；
图4(b)为采用本发明模拟得到的波场快照的固相垂直分量图；
图4(c)为采用本发明模拟得到的波场快照的流相水平分量图；
图4(d)为采用本发明模拟得到的波场快照的流相垂直分量图；
图5(a)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的固相水平分量图；
图5(b)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的固相垂直分量图；
图5(c)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的流相水平分量图；
图5(d)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的流相垂直分量图；
图6(a)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的固相水平分量图；
图6(b)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的固相垂直分量图；
图6(c)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的流相水平分量图；
图6(d)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的流相垂直分量图。
具体实施方式
应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
下面结合附图和实施例对本发明进一步说明。
如图1所示，一种黏弹各向异性双相介质区域变网格求解算子，包括以下步骤：
步骤1：输入速度场及各种参数场，建立观测系统；
具体地，输入纵波速度场、横波速度场、密度场、耗散系数场、流相弹性系数场、固体与流体之间相耦合弹性系数场、纵波品质因子场Qp、横波品质因子场Qs、各向异性参数场；
步骤2：在全局区域进行粗网格剖分，在目标区域进行细网格剖分，从而实现整个计算区域的非均匀网格划分；
具体地，全局区域的粗网格尺寸为目标区域细网格尺寸的奇数倍，如3、5、7等倍数；
步骤3：在全局区域采用黏声各向异性介质进行波场延拓；
具体地，所述步骤3中，在全局区域采用黏声各向异性波动方程(1)进行P波传播；
Figure BDA0002385881370000131
方程(1)中，p表示各向异性声波波场；q表示各向异性声波辅助波场；vpz表示对称轴方向的相速度；ε表示用以表征纵波各向异性强度的Thomson参数，δ表示用以表征动校正速度的Thomson参数；t表示时间；fx表示水平分量的震源；fz表示垂直分量的震源；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；τ表示松弛时间，τ由公式(2)求得：
τ＝τε/τσ-1 (2)
公式(2)中，τσ为应力松弛时间，τε为应变松弛时间；其中，τσ和τε与纵波品质因子Qp的关系如公式(3)、公式(4)所示；
Figure BDA0002385881370000141
Figure BDA0002385881370000142
公式(3)和公式(4)中，ω代表角频率。
具体地，所述步骤3中全局区域采用的黏声各向异性波动方程由公式(5)递推求得；
Figure BDA0002385881370000143
公式(5)中，n表示差分计算次数；F和F-1分别为傅里叶变换与傅里叶反变换；Δt是时间步长；kx为水平分量的波数；kz为垂直分量的波数；
Figure BDA0002385881370000144
为水平分量的空间二阶偏导数；
Figure BDA0002385881370000145
为垂直分量的空间二阶偏导数；
Figure BDA0002385881370000146
和
Figure BDA0002385881370000147
由公式(6)～公式(9)求得：
Figure BDA0002385881370000148
Figure BDA0002385881370000149
Figure BDA00023858813700001410
Figure BDA00023858813700001411
公式(6)～公式(9)中，Δx表示水平分量的空间网格步长；Δz表示垂直分量的空间网格步长；c2,0和c2,m为二阶偏导数的有限差分系数；m为积分变量；M为差分阶数；i为水平分量的空间坐标计数，j为垂直分量的空间坐标计数。
步骤4：在目标区域采用黏弹各向异性双相介质进行波场延拓；
具体地，所述步骤4中，在目标区域采用黏弹各向异性双相介质一阶速度应力方程(10)进行多分量波场延拓；
Figure BDA0002385881370000151
方程(10)中，ux为固相水平分量的质点速度；uz为固相垂直分量的质点速度；Ux为流相水平分量的质点速度；Uz为流相垂直分量的质点速度；τxx为水平分量的正应力；τzz为垂直分量的正应力；τxz表示切应力；S表示流相正应力；Q为固体与流体之间相耦合弹性系数；R为流相弹性系数；b表示耗散系数；t表示时间；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；L1表示体积模量的标准线性固体个数；L2表示剪切模量的标准线性固体个数；
Figure BDA0002385881370000152
表示胀缩模式的记忆变量；
Figure BDA0002385881370000153
表示体积模量的第l个标准线性固体的记忆变量；
Figure BDA0002385881370000154
表示剪切模量的第l个标准线性固体的记忆变量；
Figure BDA0002385881370000155
Figure BDA0002385881370000156
为各向异性粘弹介质的弹性系数；e11、e55为各向异性弹性系数；D1、D2和D3分别为三个中间变量；
方程(10)中，D1、D2和D3由公式(11)求得：
Figure BDA0002385881370000161
公式(11)中，ρ11表示单位体积中固体相对流体运动时固体部分总的等效质量，ρ22表示单位体积中流体相对固体运动时流体部分总的等效质量，ρ12表示单位体积中流体和固体之间的质量耦合系数；
方程(10)中，
Figure BDA0002385881370000162
和
Figure BDA0002385881370000163
由公式(12)求得；
Figure BDA0002385881370000164
公式(12)中，
Figure BDA0002385881370000165
表示体积模量的应力松弛时间常数；
Figure BDA0002385881370000166
表示剪切模量的应力松弛时间常数；
Figure BDA0002385881370000167
表示体积模量的应变松弛时间常数；
Figure BDA0002385881370000168
表示剪切模量的应变松弛时间常数；
方程(10)中，e11、e55由公式(13)求得：
Figure BDA0002385881370000169
公式(13)中，ε表示用以表征纵波各向异性强度的Thomson参数；ρ表示密度；vp表示纵波速度；vs表示横波速度；
方程(10)中，
Figure BDA00023858813700001610
由公式(14)求得：
Figure BDA0002385881370000171
公式(14)中，λ为与纵横波速度和密度有关的拉梅常数；μ为与横波速度和密度有关的拉梅常数；
Figure BDA0002385881370000172
为体积模量应力松弛时间常数与应变松弛时间常数的比值；
Figure BDA0002385881370000173
分别为剪切模量应力松弛时间常数与应变松弛时间常数的比值和；
Figure BDA0002385881370000174
和
Figure BDA0002385881370000175
由公式(15)求得；
Figure BDA0002385881370000176
具体地，所述步骤4中，目标区域采用的黏弹各向异性双相介质一阶速度应力方程(10)由公式(16)递推求得；
Figure BDA0002385881370000181
公式(16)中，Δt1为目标区域的时间采样间隔，Δt1与时间步长Δt之间的关系如公式(17)所示；
Δt＝(2k+1)Δt1 (17)
公式(17)中，k为任一正整数；
公式(16)中，Dx表示水平分量的空间一阶偏导数，Dz表示垂直分量的空间一阶偏导数，Dx和Dz由公式(18)求得：
Figure BDA0002385881370000182
公式(18)中，U为波场值，可表示ux、uz、Ux、Uz、τxx、τzz、τxz、
Figure BDA0002385881370000183
Figure BDA0002385881370000184
和
Figure BDA0002385881370000185
c1,m为一阶偏导数的有限差分系数；Δx1为目标区域水平分量的网格大小；Δz1为目标区域垂直分量的网格大小；Δx1和水平分量的空间网格步长Δx的关系及Δz1和垂直分量的空间网格步长Δz的关系如公式(19)所示；
Figure BDA0002385881370000186
步骤5：在目标区域的边界处进行黏声各向异性和黏弹各向异性双相波场传递；
具体地，所述步骤5中，在目标区域的左右边界采用公式(20)所示的边界条件进行声压P与应力τxx的波场传递；
-P＝τxx (20)
在目标区域的上下边界采用公式(21)进行声压P与应力τzz的波场传递；
-P＝τzz (21)
具体地，所述步骤5中，目标区域边界处全局区域延拓的黏声各向异性介质声波波场与目标区域延拓的黏弹各向异性双相介质弹性波波场互相传递的步骤如下：
步骤51：将目标区域上下边界上处于粗网格点上的声压P传递给对应的细网格点上的τxx，粗网格点上的uz传递给对应细网格点上的uz；将目标区域左右边界上处于粗网格点上的声压P传递给对应的细网格点上的τzz，粗网格点上的ux传递给对应细网格点上的ux；
步骤52：利用插值公式计算目标区域边界上其他细网格点上的波场值；
步骤53：采用降阶方法求解过渡区域的波场值；如采用10阶差分精度，在边界区域逐渐采用8阶、6阶、4阶、2阶差分精度；
步骤54：在目标区域其他细网格点黏弹各向异性双相介质波动方程(16)进行波场延拓；
步骤55：将目标区域的过渡区域上位于细网格点上的波场值传递给对应的粗网格；
步骤56：将目标区域上下边界上位于粗网格位置上的细网格点波场值τxx传递给粗网格点波场值P，将目标区域左右边界上位于粗网格位置上的细网格点波场值τzz传递给粗网格点波场值P，位于粗网格位置上的细网格点波场值ux和uz传递给粗网格点波场值ux和uz。
在进行波场传递时，直接一对一传递会引起不稳定，因此本发明采用的是九点加权法，利用粗网格点周围若干个细网格点波场值计算该点相应的粗网格点波场值；变时间步长的基本原理是在全局区域采用大尺度时间步长，在目标区域采用小尺度时间步长，边界区域的小尺度时间步长利用双线性插值计算得到；
步骤6：应用人工边界条件对边界反射进行吸收；
具体地，所述步骤6中，在全局区域采用PML边界条件对波场的衰减进行吸收，对目标区域的边界采用CFS-PML边界条件进行波场延拓；
使用拉伸函数将速度-应力方程变换到拉伸坐标域；在频率域中，2D拉伸函数的具体表述形式：
Figure BDA0002385881370000201
Figure BDA0002385881370000202
式中，σx代表CFS-PML中的水平分量衰减系数；σz代表CFS-PML中的垂直分量衰减系数；εx为水平分量的拉伸函数；εz为垂直分量的拉伸函数；ω为角频率；I表示复数虚单位；
笛卡尔坐标系(x，z)与拉伸坐标系
Figure BDA0002385881370000203
的坐标变换方程如下：
Figure BDA0002385881370000204
将方程(24)代入方程(10)，得到拉伸坐标系下一阶速度-应力方程的表达式：
Figure BDA0002385881370000211
其中，
Figure BDA0002385881370000212
为拉伸坐标系下的固相水平分量的质点速度；
Figure BDA0002385881370000213
为拉伸坐标系下的固相垂直分量的质点速度；
Figure BDA0002385881370000214
为拉伸坐标系下的流相水平分量的质点速度；
Figure BDA0002385881370000215
为拉伸坐标系下的流相垂直分量的质点速度；
Figure BDA0002385881370000216
表示拉伸坐标系下水平分量的正应力；
Figure BDA0002385881370000217
表示拉伸坐标系下垂直分量的正应力；
Figure BDA0002385881370000218
表示拉伸坐标系下的切应力；
Figure BDA0002385881370000219
表示拉伸坐标系下的流相正应力；
Figure BDA00023858813700002110
表示拉伸坐标系下的胀缩模式的记忆变量；
Figure BDA00023858813700002111
表示拉伸坐标系下的体积模量的第l个标准线性固体的记忆变量；
Figure BDA00023858813700002112
表示拉伸坐标系下的剪切模量的第l个标准线性固体的记忆变量。
步骤7：输出炮记录和波场快照。
本发明公开了一种黏弹各向异性双相介质区域变网格求解算子，该方法在全局区域采用黏声各向异性介质的声波算子进行波场延拓，在目标区域采用黏弹各向异性双相介质的弹性波算子进行延拓，在目标区域的边界上进行稳定的声波和弹性波场传递，并在目标区域采用细网格剖分，提高了地震波场的模拟准确性并提高了计算效率。
实施例：
本发明黏弹各向异性双相介质区域变网格求解算子，应用于黏弹各向异性双相介质模型，取得了理想的计算效果。
首先输入速度场及各种参数场，建立观测系统，其中实施例中黏弹各向异性双相介质模型如图2(a)～图2(d)所示，图2(a)为黏弹各向异性双相介质模型的双相介质参数，图2(b)为黏弹各向异性双相介质模型的速度密度参数，图2(c)为黏弹各向异性双相介质模型的品质因子参数，图2(d)为黏弹各向异性双相介质模型的各向异性参数。
然后采用本发明黏弹各向异性双相介质区域变网格求解算子进行地震波正演模拟，模型大小为3km×3km。在实施例中，全局网格间距为10m×10m，全局时间采样间隔为1ms，其中，全局区域的粗网格尺寸为目标区域细网格尺寸的3倍；记录时间为1.5s，炮点位于地表位置激发，水平位置为1500m，检波点均匀分布于地表处，共300个检波点。
采用本发明模拟得到的地震波正演模拟炮记录如图3(a)～图3(d)所示，其中图3(a)为固相水平分量，图3(b)为固相垂直分量，图3(c)为流相水平分量，图3(d)为流相垂直分量；
采用本发明模拟得到的波场快照图如图4(a)～图4(d)所示，从图中可以看出，地震波在目标区域以弹性纵横波形式传播，在其他区域以声学纯纵波形式传播，波场快照证明了本发明方法的正确性。
为了作为对比，将采用本发明模拟得到的地震波正演模拟炮记录与其他地震波正演模拟炮记录进行对比，具体如下：
图5(a)～图5(d)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录，其中图5(a)为固相水平分量，图5(b)为固相垂直分量，图5(c)为流相水平分量，图5(d)为流相垂直分量；
图6(a)～图6(d)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录，其中图6(a)为固相水平分量，图6(b)为固相垂直分量，图6(c)为流相水平分量，图6(d)为流相垂直分量；
通过将3(a)～图3(d)与图5(a)～图5(d)对比可以看出，采用本发明模拟得到的地震波正演模拟炮记录相比于各向异性双相介质地震波正演模拟炮记录，地震波的形态基本一致，但存在明显的地震波场衰减，说明本发明除了能准确模拟地震波的各向异性和双相性之外，还能模拟出地下介质中的黏弹衰减性；通过将3(a)～图3(d)与图6(a)～图6(d)进行对比可以看出，采用本发明模拟得到的地震波正演模拟炮记录相比于黏弹各向同性双相介质地震波正演模拟炮记录，振幅能量一致，不同的是，本发明能够将各向异性引起的波形准确地模拟出来。因此可以说明本发明方法同时考虑了地下介质的黏弹、各向异性和双相介质的性质，得到了准确的正演模拟结果。
本发明提出的一种黏弹各向异性双相介质区域变网格求解算子，既同时考虑了地下介质的黏弹性、各向异性和双相性，又只针对目标区域进行弹性纵横波场的模拟；为后续的复杂地质构造成像、反演工作提供地震波场计算基础。
上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (10)

1.一种黏弹各向异性双相介质区域变网格求解算子，其特征是，包括以下步骤：
步骤1：输入速度场及各种参数场，建立观测系统；
步骤2：选择目标区域进行局部网格加密；
步骤3：在全局区域采用黏声各向异性介质进行波场延拓；
步骤4：在目标区域采用黏弹各向异性双相介质进行波场延拓；
步骤5：在目标区域的边界处进行黏声各向异性和黏弹各向异性双相波场传递；
步骤6：应用人工边界条件对边界反射进行吸收；
步骤7：输出炮记录和波场快照。
2.如权利要求1所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤2中，在全局区域进行粗网格剖分，在目标区域进行细网格剖分。
3.如权利要求2所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤2中，全局区域的粗网格尺寸为目标区域细网格尺寸的奇数倍。
4.如权利要求3所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤3中，在全局区域采用黏声各向异性波动方程(1)进行P波传播；
Figure FDA0002385881360000011
方程(1)中，p表示各向异性声波波场；q表示各向异性声波辅助波场；vpz表示对称轴方向的相速度；ε表示用以表征纵波各向异性强度的Thomson参数，δ表示用以表征动校正速度的Thomson参数；t表示时间；fx表示水平分量的震源；fz表示垂直分量的震源；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；τ表示松弛时间，τ由公式(2)求得：
τ＝τε/τσ-1 (2)
公式(2)中，τσ为应力松弛时间，τε为应变松弛时间；其中，τσ和τε与纵波品质因子Qp的关系如公式(3)、公式(4)所示；
Figure FDA0002385881360000021
Figure FDA0002385881360000022
公式(3)和公式(4)中，ω代表角频率。
5.如权利要求4所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤3中全局区域采用的黏声各向异性波动方程(1)由公式(5)递推求得；
Figure FDA0002385881360000023
公式(5)中，n表示差分计算次数；F和F-1分别为傅里叶变换与傅里叶反变换；Δt是时间步长；kx为水平分量的波数；kz为垂直分量的波数；
Figure FDA0002385881360000024
为水平分量的空间二阶偏导数；
Figure FDA0002385881360000025
为垂直分量的空间二阶偏导数；
Figure FDA0002385881360000026
和
Figure FDA0002385881360000027
由公式(6)～公式(9)求得：
Figure FDA0002385881360000028
Figure FDA0002385881360000029
Figure FDA0002385881360000031
Figure FDA0002385881360000032
公式(6)～公式(9)中，Δx表示水平分量的空间网格步长；Δz表示垂直分量的空间网格步长；c2,0和c2,m为二阶偏导数的有限差分系数；m为积分变量；M为差分阶数；i为水平分量的空间坐标计数，j为垂直分量的空间坐标计数。
6.如权利要求5所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤4中，在目标区域采用黏弹各向异性双相介质一阶速度应力方程(10)进行多分量波场延拓；
Figure FDA0002385881360000033
方程(10)中，ux为固相水平分量的质点速度；uz为固相垂直分量的质点速度；Ux为流相水平分量的质点速度；Uz为流相垂直分量的质点速度；τxx为水平分量的正应力；τzz为垂直分量的正应力；τxz表示切应力；S表示流相正应力；Q为固体与流体之间相耦合弹性系数；R为流相弹性系数；b表示耗散系数；t表示时间；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；L1表示体积模量的标准线性固体个数；L2表示剪切模量的标准线性固体个数；
Figure FDA0002385881360000041
表示胀缩模式的记忆变量；
Figure FDA0002385881360000042
表示体积模量的第l个标准线性固体的记忆变量；
Figure FDA0002385881360000043
表示剪切模量的第l个标准线性固体的记忆变量；
Figure FDA0002385881360000044
Figure FDA0002385881360000045
为各向异性粘弹介质的弹性系数；e11、e55为各向异性弹性系数；D1、D2和D3分别为三个中间变量；
方程(10)中，D1、D2和D3由公式(11)求得：
Figure FDA0002385881360000046
公式(11)中，ρ11表示单位体积中固体相对流体运动时固体部分总的等效质量，ρ22表示单位体积中流体相对固体运动时流体部分总的等效质量，ρ12表示单位体积中流体和固体之间的质量耦合系数；
方程(10)中，
Figure FDA0002385881360000047
和
Figure FDA0002385881360000048
由公式(12)求得；
Figure FDA0002385881360000049
公式(12)中，
Figure FDA00023858813600000410
表示体积模量的应力松弛时间常数；
Figure FDA00023858813600000411
表示剪切模量的应力松弛时间常数；
Figure FDA00023858813600000412
表示体积模量的应变松弛时间常数；
Figure FDA00023858813600000413
表示剪切模量的应变松弛时间常数；
方程(10)中，e11、e55由公式(13)求得：
Figure FDA0002385881360000051
公式(13)中，ε表示用以表征纵波各向异性强度的Thomson参数；ρ表示密度；vp表示纵波速度；vs表示横波速度；
方程(10)中，
Figure FDA0002385881360000052
由公式(14)求得：
Figure FDA0002385881360000053
公式(14)中，λ为与纵横波速度和密度有关的拉梅常数；μ为与横波速度和密度有关的拉梅常数；
Figure FDA0002385881360000054
为体积模量应力松弛时间常数与应变松弛时间常数的比值；
Figure FDA0002385881360000055
分别为剪切模量应力松弛时间常数与应变松弛时间常数的比值和；
Figure FDA0002385881360000056
和
Figure FDA0002385881360000057
由公式(15)求得；
Figure FDA0002385881360000058
7.如权利要求6所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤4中，目标区域采用的黏弹各向异性双相介质一阶速度应力方程(10)由公式(16)递推求得；
Figure FDA0002385881360000061
公式(16)中，Δt1为目标区域的时间采样间隔，Δt1与时间步长Δt之间的关系如公式(17)所示；
Δt＝(2k+1)Δt1 (17)
公式(17)中，k为任一正整数；
公式(16)中，Dx表示水平分量的空间一阶偏导数，Dz表示垂直分量的空间一阶偏导数，Dx和Dz由公式(18)求得：
Figure FDA0002385881360000062
公式(18)中，U为波场值，可表示ux、uz、Ux、Uz、τxx、τzz、τxz、
Figure FDA0002385881360000063
Figure FDA0002385881360000064
和
Figure FDA0002385881360000065
c1,m为一阶偏导数的有限差分系数；Δx1为目标区域水平分量的网格大小；Δz1为目标区域垂直分量的网格大小；Δx1和水平分量的空间网格步长Δx的关系及Δz1和垂直分量的空间网格步长Δz的关系如公式(19)所示；
Figure FDA0002385881360000066
8.如权利要求7所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤5中，在目标区域的左右边界采用公式(20)所示的边界条件进行声压P与应力τxx的波场传递；
-P＝τxx (20)
在目标区域的上下边界采用公式(21)进行声压P与应力τzz的波场传递；
-P＝τzz (21)。
9.如权利要求8所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤5中，目标区域边界处全局区域延拓的黏声各向异性介质声波波场与目标区域延拓的黏弹各向异性双相介质弹性波波场互相传递的步骤如下：
步骤51：将目标区域上下边界上处于粗网格点上的声压P传递给对应的细网格点上的τxx，粗网格点上的uz传递给对应细网格点上的uz；将目标区域左右边界上处于粗网格点上的声压P传递给对应的细网格点上的τzz，粗网格点上的ux传递给对应细网格点上的ux；
步骤52：利用插值公式计算目标区域边界上其他细网格点上的波场值；
步骤53：采用降阶方法求解过渡区域的波场值；
步骤54：在目标区域其他细网格点黏弹各向异性双相介质波动方程(16)进行波场延拓；
步骤55：将目标区域的过渡区域上位于细网格点上的波场值传递给对应的粗网格；
步骤56：将目标区域上下边界上位于粗网格位置上的细网格点波场值τxx传递给粗网格点波场值P，将目标区域左右边界上位于粗网格位置上的细网格点波场值τzz传递给粗网格点波场值P，位于粗网格位置上的细网格点波场值ux和uz传递给粗网格点波场值ux和uz。
10.如权利要求9所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤6中，在全局区域采用PML边界条件对波场的衰减进行吸收，对目标区域的边界采用CFS-PML边界条件进行波场延拓；
使用拉伸函数将速度-应力方程变换到拉伸坐标域；在频率域中，2D拉伸函数的具体表述形式：
Figure FDA0002385881360000081
Figure FDA0002385881360000082
式中，σx代表CFS-PML中的水平分量衰减系数；σz代表CFS-PML中的垂直分量衰减系数；εx为水平分量的拉伸函数；εz为垂直分量的拉伸函数；ω为角频率；I表示复数虚单位；
笛卡尔坐标系(x，z)与拉伸坐标系
Figure FDA0002385881360000083
的坐标变换方程如下：
Figure FDA0002385881360000084
将方程(24)代入方程(10)，得到拉伸坐标系下一阶速度-应力方程的表达式：
Figure FDA0002385881360000091
其中，
Figure FDA0002385881360000092
为拉伸坐标系下的固相水平分量的质点速度；
Figure FDA0002385881360000093
为拉伸坐标系下的固相垂直分量的质点速度；
Figure FDA0002385881360000094
为拉伸坐标系下的流相水平分量的质点速度；
Figure FDA0002385881360000095
为拉伸坐标系下的流相垂直分量的质点速度；
Figure FDA0002385881360000096
表示拉伸坐标系下水平分量的正应力；
Figure FDA0002385881360000097
表示拉伸坐标系下垂直分量的正应力；
Figure FDA0002385881360000098
表示拉伸坐标系下的切应力；
Figure FDA0002385881360000099
表示拉伸坐标系下的流相正应力；
Figure FDA00023858813600000910
表示拉伸坐标系下的胀缩模式的记忆变量；
Figure FDA00023858813600000911
表示拉伸坐标系下的体积模量的第l个标准线性固体的记忆变量；
Figure FDA00023858813600000912
表示拉伸坐标系下的剪切模量的第l个标准线性固体的记忆变量。
CN202010097267.2A 2020-02-18 2020-02-18 一种黏弹各向异性双相介质区域变网格求解算子 Active CN111257930B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010097267.2A CN111257930B (zh) 2020-02-18 2020-02-18 一种黏弹各向异性双相介质区域变网格求解算子

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010097267.2A CN111257930B (zh) 2020-02-18 2020-02-18 一种黏弹各向异性双相介质区域变网格求解算子

Publications (2)

Publication Number Publication Date
CN111257930A true CN111257930A (zh) 2020-06-09
CN111257930B CN111257930B (zh) 2022-03-29

Family

ID=70945865

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010097267.2A Active CN111257930B (zh) 2020-02-18 2020-02-18 一种黏弹各向异性双相介质区域变网格求解算子

Country Status (1)

Country Link
CN (1) CN111257930B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112162038A (zh) * 2020-09-28 2021-01-01 中国石油大学(华东) 从套管井声波测井中获取套管与水泥界面剪切耦合刚度的方法
CN113406698A (zh) * 2021-05-24 2021-09-17 中国石油大学(华东) 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5999488A (en) * 1998-04-27 1999-12-07 Phillips Petroleum Company Method and apparatus for migration by finite differences
CN102183790A (zh) * 2011-02-12 2011-09-14 中国石油大学(华东) 基于时空双变网格的弹性波正演模拟技术
US20130060544A1 (en) * 2010-05-12 2013-03-07 Petrus Maria Bakker Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
CN105676280A (zh) * 2016-01-21 2016-06-15 中国矿业大学(北京) 基于旋转交错网格的双相介质地质数据获取方法和装置
US20170336522A1 (en) * 2014-12-31 2017-11-23 Landmark Graphics Corporation Seismic Elastic Wave Simulation For Tilted Transversely Isotropic Media Using Adaptive Lebedev Staggered Grid
CN108108331A (zh) * 2017-12-13 2018-06-01 国家深海基地管理中心 一种基于拟空间域弹性波方程的有限差分计算方法
CN108646293A (zh) * 2018-05-15 2018-10-12 中国石油大学(华东) 基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法
CN110609325A (zh) * 2018-06-15 2019-12-24 中国石油化工股份有限公司 弹性波场数值模拟方法及系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5999488A (en) * 1998-04-27 1999-12-07 Phillips Petroleum Company Method and apparatus for migration by finite differences
US20130060544A1 (en) * 2010-05-12 2013-03-07 Petrus Maria Bakker Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
CN102183790A (zh) * 2011-02-12 2011-09-14 中国石油大学(华东) 基于时空双变网格的弹性波正演模拟技术
US20170336522A1 (en) * 2014-12-31 2017-11-23 Landmark Graphics Corporation Seismic Elastic Wave Simulation For Tilted Transversely Isotropic Media Using Adaptive Lebedev Staggered Grid
CN105676280A (zh) * 2016-01-21 2016-06-15 中国矿业大学(北京) 基于旋转交错网格的双相介质地质数据获取方法和装置
CN108108331A (zh) * 2017-12-13 2018-06-01 国家深海基地管理中心 一种基于拟空间域弹性波方程的有限差分计算方法
CN108646293A (zh) * 2018-05-15 2018-10-12 中国石油大学(华东) 基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法
CN110609325A (zh) * 2018-06-15 2019-12-24 中国石油化工股份有限公司 弹性波场数值模拟方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
曲英铭等: "基于多尺度双变网格的时间域全波形反演", 《石油物探》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112162038A (zh) * 2020-09-28 2021-01-01 中国石油大学(华东) 从套管井声波测井中获取套管与水泥界面剪切耦合刚度的方法
CN112162038B (zh) * 2020-09-28 2022-07-29 中国石油大学(华东) 从套管井声波测井中获取套管与水泥界面剪切耦合刚度的方法
CN113406698A (zh) * 2021-05-24 2021-09-17 中国石油大学(华东) 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法

Also Published As

Publication number Publication date
CN111257930B (zh) 2022-03-29

Similar Documents

Publication Publication Date Title
CN104122585B (zh) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
Manolis et al. Elastic waves in continuous and discontinuous geological media by boundary integral equation methods: A review
Ladopoulos Elastodynamics for Non-linear Seismic Wave Motion in Real-Time Expert Seismology
CN108646293B (zh) 基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法
CN111257930B (zh) 一种黏弹各向异性双相介质区域变网格求解算子
CN104965223B (zh) 粘声波全波形反演方法及装置
CN105652321A (zh) 一种粘声各向异性最小二乘逆时偏移成像方法
CN112327358B (zh) 一种粘滞性介质中声波地震数据正演模拟方法
CN109946742B (zh) 一种TTI介质中纯qP波地震数据模拟方法
CN109143340A (zh) 一种基于常q模型的粘弹介质地震波模拟方法及系统
CN105652319A (zh) 一种复杂介质近地表地层q值的估计方法
Balasubramanyam et al. A finite-difference simulation of ultrasonic Lamb waves in metal sheets with experimental verification
CN104965222A (zh) 三维纵波阻抗全波形反演方法及装置
CN110888159B (zh) 基于角度分解与波场分离的弹性波全波形反演方法
Zhang et al. Pseudospectral modeling and dispersion analysis of Rayleigh waves in viscoelastic media
Huang et al. A multi-block finite difference method for seismic wave equation in auxiliary coordinate system with irregular fluid–solid interface
Blanch et al. Viscoelastic finite-difference modeling
CN109738944B (zh) 基于广角反射的地震采集参数确定方法及装置
Shivaprasad et al. Numerical modelling methods for ultrasonic wave propagation through polycrystalline materials
CN110658558A (zh) 吸收衰减介质叠前深度逆时偏移成像方法及系统
Wang et al. Improved hybrid iterative optimization method for seismic full waveform inversion
CN110764146B (zh) 一种基于声波算子的空间互相关弹性波反射波形反演方法
CN115327626A (zh) 一种用于三维ati介质的弹性波场矢量分解方法及系统
CN113866823A (zh) 一种粘声各向异性介质中的正演成像方法
CN114139335A (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