阅读说明：本技术 一种基于时域高阶有限差分法的高精度波场数值模拟方法 (High-precision wave field numerical simulation method based on time domain high-order finite difference method ) 是由 梁晨曦 张固澜 段景 李勇 罗一梁 何承杰 罗帆 杜皓 曾梦 詹熠宗 王佳 于 2021-03-24 设计创作，主要内容包括：本发明公开了一种基于时域高阶有限差分法的高精度波场数值模拟方法,包括以下步骤：S1.在波场方程的基础上,对空间2M阶展开,发展得到的传统时间二阶空间高阶的有限差分方法；S2.基于均匀网格有限差分方法以及Taylor公式,对二阶时间导数进行N阶精度的展开,得到奇数时刻时间高阶差分求解公式,并将奇数时刻时间高阶有限差分表示；S3.基于均匀网格有限差分方法以及Taylor公式,对二阶时间导数进行N阶精度的展开,得到偶数时刻时间高阶差分求解公式,并将偶数时刻时间高阶有限差分表示；S4.对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理。本发明在保证高稳定性的同时有效的压制了频散,可提供精确稳定的波场模拟结果。(The invention discloses a high-precision wave field numerical simulation method based on a time domain high-order finite difference method, which comprises the following steps of: s1, on the basis of a wave field equation, a space 2M order is expanded, and an obtained traditional time second-order space high-order finite difference method is developed; s2, based on a uniform grid finite difference method and a Taylor formula, carrying out N-order precision expansion on a second-order time derivative to obtain an odd-number moment time high-order difference solving formula, and expressing the odd-number moment time high-order finite difference; s3, based on a uniform grid finite difference method and a Taylor formula, carrying out N-order precision expansion on the second-order time derivative to obtain an even-number moment time high-order difference solving formula, and expressing the even-number moment time high-order finite difference; and S4, carrying out odd-even weighting processing on the finite difference method for simulating odd-numbered time and even-numbered time of the wave field. The invention effectively suppresses frequency dispersion while ensuring high stability and can provide accurate and stable wave field simulation results.)

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

技术领域

本发明涉及地震勘探技术领域，特别是涉及一种基于时域高阶有限差分法的高精度波场数值模拟方法。

背景技术

波波动方程数值模拟已广泛应用于地震正演、逆时偏移、全波形反演等各个方面，当前，波动方程数值模拟的方法主要包括基于傅里叶变换的伪谱法，基于非规则网格剖分的有限元类方法以及基于差分近似的有限差分方法。伪谱法采用傅里叶变换计算空间偏导数，能够压制由于空间离散造成的数值频散；有限元类方法能准确地模拟各种不规则边界，从而避免矩形网格剖分造成的阶梯绕射。但是这两种方法均需要进行大量的计算。相比之下，有限差分方法因为计算效率高、所需内存小、实现简单而广泛应用于地震正演研究。

由于网格的离散，便会从时间，空间两个方面产生频散问题，长期以来解决因网格离散造成空间频散的方法较多，现目前空间已达到2M阶精度，但在时间精度的问题上，任存在改进空间。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于时域高阶有限差分法的高精度波场数值模拟方法，针对由于时间精度不足造成的频散，给出了基于Taylor公式和奇偶加权算法求解差分系数，在保证高稳定性的同时有效的压制了频散，可提供精确稳定的波场模拟结果。

本发明的目的是通过以下技术方案来实现的：一种基于时域高阶有限差分法的高精度波场数值模拟方法，包括以下步骤：

S1.在波场方程的基础上，对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法；

S2.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式,并将奇数时刻时间高阶有限差分表示；

S3.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式,并将偶数时刻时间高阶有限差分表示；

S4.对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理。

所述步骤S1包括：

S101.给出波动方程

其中v是传播速度，p是压力波场，t是传播时间，x和z分别表示空间x轴和z轴；

S102.通过对空间2M阶展开，传统时间二阶空间高阶的有限差分方法可以表示为：

其中γ_m是空间差分系数，是在t时刻，空间位置(x，z)处的离散化波场，Δx、Δz是离散空间间距，Δt是离散时间间距，是离散点(x，z)处的速度。

所述步骤S2包括：

S201.将奇数时刻时间高阶有限差分表示为：

其中，

λ_n表示奇数时刻时间高阶的差分系数，下标n＝1,2,3,...,N；

S202.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式：

则t+Δt时刻波场为：

所述步骤S3包括：

S301.将偶数时刻时间高阶有限差分表示为：

μ_n表示偶数时刻时间高阶的差分系数，下标n＝1,2,3,...,N；

S302.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式：

则得到t+Δt时刻波场为：

所述步骤S4包括：

对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理，得到：

其中，

β_n＝ηλ_n+(1-η)μ_n,n＝1,2,...,N，β_n表示奇偶权重。

本发明的有益效果是：本发明的波场模拟均匀网格有限差分方法适用于任何介质的波动方程的波场模拟，针对其因网格离散造成的频散现象，给出了较好的压制频散的方法。

附图说明

图1为本发明的方法流程图；

图2为实施例中本发明与其他方法的波场快照对比示意图；

图3为实施例中本发明的地震记录示意图。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。

如图1所示，一种基于时域高阶有限差分法的高精度波场数值模拟方法，包括以下步骤：

S1.在波场方程的基础上，对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法；

S2.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到奇数时刻时间高阶差分求解公式,并将奇数时刻时间高阶有限差分表示；

S3.基于均匀网格有限差分方法以及Taylor公式，对二阶时间导数进行N阶精度的展开，得到偶数时刻时间高阶差分求解公式,并将偶数时刻时间高阶有限差分表示；

S4.对波场模拟奇数时刻和偶数时刻有限差分方法进行奇偶加权处理。

其中，所述步骤S1包括：

S101.给出波动方程

其中v是传播速度，p是压力波场，t是传播时间，x和z分别表示空间x轴和z轴；

S102.通过对空间2M阶展开，发展得到的传统时间二阶空间高阶的有限差分方法：

其中γ_m是空间差分系数，是在t时刻，空间位置(x，z)处的离散化波场，Δx、Δz是离散空间间距，Δt是离散时间间距；是离散点(x，z)处的速度；

对偶数时刻，我们对时间给出表达式：

对其各项进行Taylor展开，并求解差分系数：

则可得t+Δt时刻波场：

同理，对奇数时刻：

则可得t+Δt时刻波场

对其进行奇偶加权可得：

如图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模型，取得了较为理想的应用效果。

上述说明示出并描述了本发明的一个优选实施例，但如前所述，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。