CN112859170B - 一种基于时域高阶有限差分法的高精度波场数值模拟方法 - Google Patents

一种基于时域高阶有限差分法的高精度波场数值模拟方法 Download PDF

Info

Publication number
CN112859170B
CN112859170B CN202110311151.9A CN202110311151A CN112859170B CN 112859170 B CN112859170 B CN 112859170B CN 202110311151 A CN202110311151 A CN 202110311151A CN 112859170 B CN112859170 B CN 112859170B
Authority
CN
China
Prior art keywords
order
time
finite difference
wave field
odd
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.)
Expired - Fee Related
Application number
CN202110311151.9A
Other languages
English (en)
Other versions
CN112859170A (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.)
Southwest Petroleum University
Original Assignee
Southwest Petroleum University
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 Southwest Petroleum University filed Critical Southwest Petroleum University
Priority to CN202110311151.9A priority Critical patent/CN112859170B/zh
Publication of CN112859170A publication Critical patent/CN112859170A/zh
Application granted granted Critical
Publication of CN112859170B publication Critical patent/CN112859170B/zh
Expired - Fee Related 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • 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

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)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于时域高阶有限差分法的高精度波场数值模拟方法,包括以下步骤：S1.在波场方程的基础上，对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法；S2.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式,并将奇数时刻时间高阶有限差分表示；S3.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式,并将偶数时刻时间高阶有限差分表示；S4.对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理。本发明在保证高稳定性的同时有效的压制了频散，可提供精确稳定的波场模拟结果。

Description

一种基于时域高阶有限差分法的高精度波场数值模拟方法
技术领域
本发明涉及地震勘探技术领域，特别是涉及一种基于时域高阶有限差分法的高精度波场数值模拟方法。
背景技术
波波动方程数值模拟已广泛应用于地震正演、逆时偏移、全波形反演等各个方面，当前，波动方程数值模拟的方法主要包括基于傅里叶变换的伪谱法，基于非规则网格剖分的有限元类方法以及基于差分近似的有限差分方法。伪谱法采用傅里叶变换计算空间偏导数，能够压制由于空间离散造成的数值频散；有限元类方法能准确地模拟各种不规则边界，从而避免矩形网格剖分造成的阶梯绕射。但是这两种方法均需要进行大量的计算。相比之下，有限差分方法因为计算效率高、所需内存小、实现简单而广泛应用于地震正演研究。
由于网格的离散，便会从时间，空间两个方面产生频散问题，长期以来解决因网格离散造成空间频散的方法较多，现目前空间已达到2M阶精度，但在时间精度的问题上，任存在改进空间。
发明内容
本发明的目的在于克服现有技术的不足，提供一种基于时域高阶有限差分法的高精度波场数值模拟方法，针对由于时间精度不足造成的频散，给出了基于Taylor公式和奇偶加权算法求解差分系数，在保证高稳定性的同时有效的压制了频散，可提供精确稳定的波场模拟结果。
本发明的目的是通过以下技术方案来实现的：一种基于时域高阶有限差分法的高精度波场数值模拟方法，包括以下步骤：
S1.在波场方程的基础上，对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法；
S2.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式,并将奇数时刻时间高阶有限差分表示；
S3.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式,并将偶数时刻时间高阶有限差分表示；
S4.对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理。
所述步骤S1包括：
S101.给出波动方程
Figure BDA0002989642410000021
其中v是传播速度，p是压力波场，t是传播时间，x和z分别表示空间x轴和z轴；
S102.通过对空间2M阶展开，传统时间二阶空间高阶的有限差分方法可以表示为：
Figure BDA0002989642410000022
其中γm是空间差分系数，
Figure BDA0002989642410000023
是在t时刻，空间位置(x，z)处的离散化波场，Δx、Δz是离散空间间距，Δt是离散时间间距，
Figure BDA0002989642410000024
是离散点(x，z)处的速度。
所述步骤S2包括：
S201.将奇数时刻时间高阶有限差分表示为：
Figure BDA0002989642410000025
其中，
Figure BDA0002989642410000026
λn表示奇数时刻时间高阶的差分系数，下标n＝1,2,3,...,N；
S202.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式：
Figure BDA0002989642410000027
则t+Δt时刻波场为：
Figure BDA0002989642410000028
所述步骤S3包括：
S301.将偶数时刻时间高阶有限差分表示为：
Figure BDA0002989642410000031
μn表示偶数时刻时间高阶的差分系数，下标n＝1,2,3,...,N；
S302.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式：
Figure BDA0002989642410000032
则得到t+Δt时刻波场为：
Figure BDA0002989642410000033
所述步骤S4包括：
对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理，得到：
Figure BDA0002989642410000034
其中，
βn＝ηλn+(1-η)μn,n＝1,2,...,N，βn表示奇偶权重。
本发明的有益效果是：本发明的波场模拟均匀网格有限差分方法适用于任何介质的波动方程的波场模拟，针对其因网格离散造成的频散现象，给出了较好的压制频散的方法。
附图说明
图1为本发明的方法流程图；
图2为实施例中本发明与其他方法的波场快照对比示意图；
图3为实施例中本发明的地震记录示意图。
具体实施方式
下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。
如图1所示，一种基于时域高阶有限差分法的高精度波场数值模拟方法，包括以下步骤：
S1.在波场方程的基础上，对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法；
S2.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式,并将奇数时刻时间高阶有限差分表示；
S3.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式,并将偶数时刻时间高阶有限差分表示；
S4.对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理。
其中，所述步骤S1包括：
S101.给出波动方程
Figure BDA0002989642410000041
其中v是传播速度，p是压力波场，t是传播时间，x和z分别表示空间x轴和z轴；
S102.通过对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法：
Figure BDA0002989642410000042
其中γm是空间差分系数，
Figure BDA0002989642410000043
是在t时刻，空间位置(x，z)处的离散化波场，Δx、Δz是离散空间间距，Δt是离散时间间距；
Figure BDA0002989642410000044
是离散点(x，z)处的速度；
对偶数时刻，我们对时间
Figure BDA0002989642410000045
给出表达式：
Figure BDA0002989642410000046
对其各项进行Taylor展开，并求解差分系数：
Figure BDA0002989642410000051
则可得t+Δt时刻波场：
Figure BDA0002989642410000052
同理，对奇数时刻：
则可得t+Δt时刻波场
Figure BDA0002989642410000053
对其进行奇偶加权可得：
Figure BDA0002989642410000054
如图2所示，在本申请的实施例中，采用的模型网格大小为8m，共501个样点，时间采样间隔为1ms，本发明方法采用时间4阶精度，传统方法采用时间2阶精度，主频为25Hz，速度为2000m/s，同时空间阶数2M＝16的均匀介质模型，图2中，(a)为均匀介质常规方法2ms波场快照(b)为均匀介质常规方法2ms波场快照局部(c)为均匀介质本发明方法2ms波场快照(d)为均匀介质本发明方法2ms波场快照局部，通过图2中的传统时间二阶方法与本发明提出的时间高阶方法的对比可得，在相同的模型参数下，本发明方法得到的波场具有更高的精度
如图3所示，在本申请的实施例中，采用nx*nz＝1001*701，网格大小为10m，时间采样间隔为1ms，时间4阶精度，主频为16Hz，同时空间阶数2M＝16的Marmousi模型，图3中(a)为Marmousi模型的介质参数示意图，(b)为Marmousi模型地面地震记录，(c)为Marmousi模型vsp记录。通过图3可以看出，本发明方法针对复杂的marmousi模型，取得了较为理想的应用效果。
上述说明示出并描述了本发明的一个优选实施例，但如前所述，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (5)

1.一种基于时域高阶有限差分法的高精度波场数值模拟方法，其特征在于：包括以下步骤：
S1.在波场方程的基础上，对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法；
S2.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式,并将奇数时刻时间高阶进行有限差分表示；
S3.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式,并将偶数时刻时间高阶进行有限差分表示；
S4.对波场模拟奇数时刻和偶数时刻的有限差分表示结果进行奇偶加权处理。
2.根据权利要求1所述的一种基于时域高阶有限差分法的高精度波场数值模拟方法，其特征在于：所述步骤S1包括：
S101.给出波动方程
Figure FDA0003683264820000011
其中v是传播速度，p是压力波场，t是传播时间，x和z分别表示空间x轴和z轴；
S102.通过对空间2M阶展开，传统时间二阶空间高阶的有限差分方法表示为：
Figure FDA0003683264820000012
其中γm是空间差分系数，
Figure FDA0003683264820000013
是在t时刻，空间位置(x，z)处的离散化波场，Δx是x轴方向离散空间间距；Δz是z轴方向离散空间间距，Δt是离散时间间距，vx,z是离散点(x，z)处的速度。
3.根据权利要求2所述的一种基于时域高阶有限差分法的高精度波场数值模拟方法，其特征在于：所述步骤S2包括：
S201.将奇数时刻时间高阶有限差分表示为：
Figure FDA0003683264820000014
其中，
Figure FDA0003683264820000015
λn表示奇数时刻时间高阶的差分系数，下标n＝1,2,3,...,N；
S202.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式：
Figure FDA0003683264820000021
则t+Δt时刻波场为：
Figure FDA0003683264820000022
4.根据权利要求3所述的一种基于时域高阶有限差分法的高精度波场数值模拟方法，其特征在于：所述步骤S3包括：
S301.将偶数时刻时间高阶有限差分表示为：
Figure FDA0003683264820000023
μn表示偶数时刻时间高阶的差分系数，下标n＝1,2,3,...,N；
S302.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式：
Figure FDA0003683264820000024
则得到t+Δt时刻波场为：
Figure FDA0003683264820000031
5.根据权利要求4所述的一种基于时域高阶有限差分法的高精度波场数值模拟方法，其特征在于：所述步骤S4包括：
对波场模拟奇数时刻和偶数时刻的有限差分表示结果进行奇偶加权处理，得到：
Figure FDA0003683264820000032
其中，
βn＝ηλn+(1-η)μn,n＝2,3,...,N，βn表示奇偶权重。
CN202110311151.9A 2021-03-24 2021-03-24 一种基于时域高阶有限差分法的高精度波场数值模拟方法 Expired - Fee Related CN112859170B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110311151.9A CN112859170B (zh) 2021-03-24 2021-03-24 一种基于时域高阶有限差分法的高精度波场数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110311151.9A CN112859170B (zh) 2021-03-24 2021-03-24 一种基于时域高阶有限差分法的高精度波场数值模拟方法

Publications (2)

Publication Number Publication Date
CN112859170A CN112859170A (zh) 2021-05-28
CN112859170B true CN112859170B (zh) 2022-08-02

Family

ID=75992458

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110311151.9A Expired - Fee Related CN112859170B (zh) 2021-03-24 2021-03-24 一种基于时域高阶有限差分法的高精度波场数值模拟方法

Country Status (1)

Country Link
CN (1) CN112859170B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103853930A (zh) * 2014-03-19 2014-06-11 中国科学院地质与地球物理研究所 地震矢量波场数值模拟方法和装置
CN107942375A (zh) * 2017-11-17 2018-04-20 河海大学 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法
CN110542923A (zh) * 2019-09-02 2019-12-06 成都理工大学 一种快速高精度叠后地震阻抗反演方法
CN111273345A (zh) * 2020-03-05 2020-06-12 西南石油大学 一种基于高精度时频瞬时相位的地震数据时频谱处理方法
CN111639429A (zh) * 2020-05-29 2020-09-08 中国人民解放军国防科技大学 基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011141440A1 (en) * 2010-05-12 2011-11-17 Shell Internationale Research Maatschappij B.V. Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
US9348044B2 (en) * 2012-04-19 2016-05-24 Cgg Services Sa Vectorization of fast fourier transform for elastic wave propogation for use in seismic underwater exploration of geographical areas of interest
CN108802819B (zh) * 2018-06-26 2019-11-08 西安交通大学 一种深度均匀采样梯形网格有限差分地震波场模拟方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103853930A (zh) * 2014-03-19 2014-06-11 中国科学院地质与地球物理研究所 地震矢量波场数值模拟方法和装置
CN107942375A (zh) * 2017-11-17 2018-04-20 河海大学 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法
CN110542923A (zh) * 2019-09-02 2019-12-06 成都理工大学 一种快速高精度叠后地震阻抗反演方法
CN111273345A (zh) * 2020-03-05 2020-06-12 西南石油大学 一种基于高精度时频瞬时相位的地震数据时频谱处理方法
CN111639429A (zh) * 2020-05-29 2020-09-08 中国人民解放军国防科技大学 基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
A Second-Order Uniformly Stable Explicit Asymmetric Discretization Method for One-Dimensional Fractional Diffusion Equations;Zhu, Lin;《COMPLEXITY》;20190621 *
基于新的差分结构的时-空域高阶有限差分波动方程数值模拟方法;张保庆等;《地球物理学报》;20160515;第59卷(第05期);1804-1814 *
基于高精度时频瞬时相位谱的多尺度曲率及其应用;李勇 等;《石油物探》;20190331;第58卷(第2期);第253-261页 *
基于高阶交错网格地震波数值频散压制方法;范飞;《辽宁化工》;第45卷(第04期);458-460 *
有限差分法正演频散衰减及叠前逆时偏移的研究;刘洋;《中国优秀硕士学位论文全文数据库 基础科学辑》;20131215;A011-553 *
非均匀介质地震波传播交错网格高阶有限差分法模拟;裴正林等;《石油大学学报(自然科学版)》;20040130;第27卷(第06期);17-21 *

Also Published As

Publication number Publication date
CN112859170A (zh) 2021-05-28

Similar Documents

Publication Publication Date Title
CN104122585B (zh) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
Komatitsch et al. The spectral-element method in seismology
CN108763683B (zh) 一种三角函数框架下新weno格式构造方法
Sun et al. Waveform inversion via nonlinear differential semblance optimization
CN104965222A (zh) 三维纵波阻抗全波形反演方法及装置
Xie et al. Singularity cancellation method for time-domain boundary element formulation of elastodynamics: A direct approach
CN105652319A (zh) 一种复杂介质近地表地层q值的估计方法
CN110988988B (zh) 基于垂直裂缝介质的地震波波场模拟方法及装置
CN109946742A (zh) 一种TTI介质中纯qP波地震数据模拟方法
Muratov et al. Grid-characteristic method on unstructured tetrahedral meshes
CN112805598B (zh) 拓展有限差分稳定性条件的波场模拟方法、设备及介质
CN112859170B (zh) 一种基于时域高阶有限差分法的高精度波场数值模拟方法
CN114139335A (zh) 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法
CN111257930B (zh) 一种黏弹各向异性双相介质区域变网格求解算子
CN113311484B (zh) 利用全波形反演获取粘弹性介质弹性参数的方法及装置
CN109725346B (zh) 一种时间-空间高阶vti矩形网格有限差分方法及装置
CN113221392A (zh) 一种非均匀粘声波在无限域内传播模型的构建方法
CN116822297B (zh) 一种应用于弹性波正演的三阶Higdon阻尼吸收边界方法
Assonitis et al. A shock-fitting technique for 2D/3D flows with interactions using structured grids
AlMuhaidib et al. Finite difference elastic wave modeling including surface topography
Capon et al. On the adaptive finite element solution of partial differential equations using hr-refinement
CN110646839A (zh) 地震波模拟方法及系统
CN109307888A (zh) 波场模拟的自由面处理方法及系统
CN112578431B (zh) 一种限存状态全波形反演波场最优化存储方法及系统
CN114924313B (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: 20220802