CN105911584B - 一种隐式交错网格有限差分弹性波数值模拟方法及装置 - Google Patents

Abstract

本发明公开了一种隐式交错网格有限差分弹性波数值模拟方法及装置。其中，该方法包括：基于隐式交错网格有限差分格式和平面波理论公式，得到频散关系表达式；在频散关系表达式中β的取值范围内均匀地取(M+1)个离散点，并设置在离散点上频散关系表达式左右两边相等，得到线性代数方程组；求解线性代数方程组，得到隐式交错网格有限差分系数a_m和常量b；基于隐式交错网格有限差分系数a_m、常量b和隐式交错网格有限差分格式，求解弹性波动方程；根据求解结果生成波场图，并根据波场图进行波场分析。通过本发明能减小弹性波动方程的求解误差，提高其数值解精度，即可得到更准确的波场模拟记录。当约束相同的计算误差时，本发明能减少计算时长，提高计算效率。

Description

一种隐式交错网格有限差分弹性波数值模拟方法及装置

技术领域

本发明涉及地震波勘探技术领域，尤其涉及一种隐式交错网格有限差分弹性波数值模拟方法及装置。

背景技术

地震波场的数值模拟是连结以反射波为主要研究课题的反射地震学和以岩层为主要研究课题的钻井、测井、地质和油藏工程等学科的纽带，其技术的应用几乎贯穿于地震资料采集、处理、解释和油气藏开发工程的各个环节。数值模拟方法就是对特定的地质、地球物理问题作适当的简化，从而形成数学模型，通过数值方法求解该模型的数学物理方程，进而据此模拟研究地震波传播。地震波数值模拟研究各种复杂地质条件下地震波的传播，对认识地震波的规律和特征，解释实际地震资料等均具有重要的理论和实际意义，已在地震勘探和天然地震领域中得到广泛的应用。

地震波数值模拟的核心问题就是如何求解波动方程。常用的地震波波动方程求解方法有有限差分法、有限元法、伪谱法等。其中，由于有限差分法具有很强的灵活性，并且容易实施，其已经被广泛地应用到地震波数值模拟和偏移成像中。有限差分法中最常用的方法包括规则网格有限差分法和交错网格有限差分法等，由于交错网格有限差分法具有更高的精度和稳定性，其在数值模拟中更受欢迎。

有限差分法的计算精度和效率很大程度上取决于计算空间导数的格式，其格式主要包括显格式和隐格式。其中，显格式具有较小的计算量，而隐格式具有更高的精度和稳定性。因此，针对隐格式有限差分，学者们提出了一系列实现方法，例如：针对空间二阶导数在各个方向上的数值解，提出了一种分裂的隐式差分格式；通过采用Fourier分析和最小二乘计算差分系数，为空间二阶导数的求解提出了一系列的优化的隐格式；同时考虑波动方程的时间导数项和空间导数项，提出了一种新的隐式差分格式，并证明了隐格式能够极大地提高计算精度。

然而，上述这些方法存在的问题是：1)上述这些方法仅仅用于声波模拟，而真实的地下介质更接近于弹性介质，在岩石中产生的机械振动可以看成是弹性介质中的弹性振动，地震波可以看成是在岩层中传播的弹性波，因而发展弹性波数值模拟更具有发展意义；2)上述这些方法难以运用于交错网格空间一阶导数的数值求解。

下面对现有技术中的地震勘探中的地震波数值模拟方法进行介绍。

地震波数值模拟方法的核心问题是通过数值方法求解模型的地震波波动方程，模拟研究地震波传播。

在均匀各向同性介质中，二维弹性波动方程(公式1)可表示为：

其中，t表示时间，x和z表示空间坐标，(v_x,v_z)是速度矢量，(τ_xx,τ_zz,τ_xz)是应力张量，ρ是密度，λ和μ是Lame常数。弹性波动方程(公式1)中的时间导数项(即公式1中等号左边各项)通过二阶差分格式求解，而为了提高求解精度，空间导数项(即公式1中等号右边各项)通过下述隐式交错网格有限差分格式(公式2)求解。

函数p(x)的空间一阶导数展开为2M阶的隐式交错网格有限差分格式，如下：

其中，代表一阶差分算子，代表二阶差分算子，b是任意一个常量，h是网格大小，a_m(m＝1,2,...,M)是隐式交错网格的差分系数，M是阶数，m＝1,2,……,M。

令

将公式3代入公式2中，得：

由平面波理论公式(公式5)：

p(x)＝p₀e^ikx， (5)

其中，p₀是一个常数，k是波数，e是自然对数的底数。

将公式5代入公式4，得：

令

β＝kh/2， (7)

则公式6可以化简为：

其中，β∈[0,π/2]，公式8是频散关系表达式。

现有技术差分格式中的差分系数通常采用泰勒级数展开法获得，即对于频散关系表达式(公式8)，将其左右两端作泰勒级数展开，得：

对比β两端的系数，得：

上式可以被写为如下的矩阵形式：

由泰勒级数展开法获得的隐式交错网格有限差分系数可以通过求解上述线性方程组获得。将此系数代入隐式交错网格有限差分格式(公式2)，并用其求解弹性波动方程(公式1)，便可以进行数值模拟。

泰勒级数展开法对频散关系表达式(公式8)两端的表达式作泰勒展开时，只取展开项的前几项，由于截断将会带来误差。同时，在采用泰勒级数展开法求解差分格式的系数时，并没有考虑频率的范围，因而其只能在局部获得较好的精度。

定义如下衡量频散大小的表达式：

该表达式表明，当δ_M(β)＝1，不存在数值频散；当δ_M(β)与1存在较大的差值，将会出现较大的数值频散。

图1是根据相关技术的基于泰勒级数展开法获得的差分系数的数值频散曲线示意图，如图1所示，尽管泰勒级数展开法获得的隐式交错网格差分系数在中低频存在较好的精度，但其在高频区域急速下降，存在极大的误差。

图2和图3分别是根据相关技术的基于泰勒级数展开法的8阶和16阶隐式交错网格有限差分方法获得的1.5s时刻的波场快照示意图，其地质模型和数值模拟参数如下：模型为双层各向同性介质，大小为8000m×8000m，网格间距为20m×20m，分界面深度在5000m，界面上下的纵波速度分别为2100m/s和2600m/s，横波速度1200m/s和1500m/s，密度分别为1800kg/m3和1900kg/m3，震源是主频为20Hz的雷克子波，置于模型的中间，时间步长为0.001s。在数值模拟过程中，对空间导数分别采用8阶和16阶格式进行计算，并展示了x分量和z分量。如图2所示，8阶的模拟结果存在着明显的频散现象，其预示着较大的误差。如图3所示，16阶的差分格式提高了模拟精度，但仍具有少量的频散，并且由于需要更多的网格点参与运算，消耗了更多的计算量。

如上所述，传统的基于泰勒级数展开法的隐式交错网格有限差分方法因在中高频的精度较差，导致地震波数值模拟方法存在严重的误差。

针对相关技术中地震波数值模拟方法的精度较差的问题，目前尚未提出有效的解决方案。

发明内容

本发明提供了一种隐式交错网格有限差分弹性波数值模拟方法及装置，以至少解决相关技术中地震波数值模拟方法的精度较差的问题。

根据本发明的一个方面，提供了一种隐式交错网格有限差分弹性波数值模拟方法，其中，该方法包括：基于隐式交错网格有限差分格式和平面波理论公式，得到频散关系表达式：其中，β＝kh/2，k是波数，h是网格大小，b是常量，M是阶数，a_m是隐式交错网格有限差分系数，m＝1,2,……,M；在所述β的取值范围内均匀地取(M+1)个离散点：β(1),β(2),...β(M+1)，并设置在所述离散点上所述频散关系表达式左右两边相等，得到线性代数方程组；求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b；基于所述隐式交错网格有限差分系数a_m、所述常量b和所述隐式交错网格有限差分格式，求解弹性波动方程，生成求解结果；根据所述求解结果生成波场图，并根据所述波场图进行波场分析。

优选地，所述β的取值范围是：[0,π/2]，或者[0,π/2]中的任一区间。

优选地，所述线性代数方程组是：

优选地，求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b，包括：通过高斯Gauss消元法或者LU分解法求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b。

优选地，基于所述隐式交错网格有限差分系数a_m、所述常量b和所述隐式交错网格有限差分格式，求解弹性波动方程，生成求解结果，包括：将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，求解所述弹性波动方程中的空间偏导数；将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，经过时间迭代，生成所述求解结果。

优选地，将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，求解所述弹性波动方程中的空间偏导数，包括：将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，将隐式交错网格有限差分格式运用到每一个网格点，联立所有网格点上的差分方程，得到三对角矩阵方程；求解所述三对角矩阵方程，得到所述弹性波动方程中的空间偏导数；其中，所述隐式交错网格有限差分方程是经所述隐式交错网格有限差分格式转换得到。

优选地，将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，经过时间迭代，生成所述求解结果，包括：将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，通过对时间偏导数采用二阶有限差分，对时间进行迭代，得到任一时刻空间任意网格点的数值，以生成所述求解结果。

根据本发明的另一个方面，提供了一种隐式交错网格有限差分弹性波数值模拟装置，其中，该装置包括：频散关系计算模块，用于基于隐式交错网格有限差分格式和平面波理论公式，得到频散关系表达式：其中，β＝kh/2，k是波数，h是网格大小，b是常量，M是阶数，a_m是隐式交错网格有限差分系数，m＝1,2,……,M；方程组构造模块，用于在所述β的取值范围内均匀地取(M+1)个离散点：β(1),β(2),...β(M+1)，并设置在所述离散点上所述频散关系表达式左右两边相等，得到线性代数方程组；系数计算模块，用于求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b；方程求解模块，用于基于所述隐式交错网格有限差分系数a_m、所述常量b和所述隐式交错网格有限差分格式，求解弹性波动方程，生成求解结果；模拟分析模块，用于根据所述求解结果生成波场图，并根据所述波场图进行波场分析。

优选地，所述β的取值范围是：[0,π/2]，或者[0,π/2]中的任一区间。

优选地，所述线性代数方程组是：

优选地，所述系数计算模块，还用于通过高斯Gauss消元法或者LU分解法求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b。

优选地，所述方程求解模块，包括：空间偏导数计算单元，用于将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，求解所述弹性波动方程中的空间偏导数；时间迭代单元，用于将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，经过时间迭代，生成所述求解结果。

优选地，所述空间偏导数计算单元，具体用于将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，将隐式交错网格有限差分格式运用到每一个网格点，联立所有网格点上的隐式交错网格有限差分方程，得到三对角矩阵方程；求解所述三对角矩阵方程，得到所述弹性波动方程中的空间偏导数；其中，隐式交错网格有限差分方程是经隐式交错网格有限差分格式转换得到。

优选地，所述时间迭代单元，具体用于将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，通过对时间偏导数采用二阶有限差分，对时间进行迭代，得到任一时刻空间任意网格点的数值，以生成所述求解结果。

本发明提出了一种具有更高模拟精度的基于采样逼近法的隐式交错网格有限差分方案。与传统的基于泰勒级数展开法的隐式交错网格有限差分方法相比，当离散参数和算子长度相同时，采用本发明的方案能够减小弹性波动方程的求解误差，提高其数值解精度。当约束相同的计算误差时，采用本发明的方案可采用较少的差分算子，从而减少计算时长，提高计算效率。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的限定。在附图中：

图1是根据相关技术的基于泰勒级数展开法获得的差分系数的数值频散曲线示意图；

图2是根据相关技术的基于泰勒级数展开法的8阶隐式交错网格有限差分方法获得的1.5s时刻的波场快照示意图；

图3是根据相关技术的基于泰勒级数展开法的16阶隐式交错网格有限差分方法获得的1.5s时刻的波场快照示意图；

图4是本发明一实施例的隐式交错网格有限差分弹性波数值模拟方法的流程图；

图5是本发明一实施例的采用采样逼近法的隐式交错网格有限差分弹性波数值模拟方法的流程图；

图6是本发明一实施例的基于采样逼近法获得的差分系数的数值频散曲线示意图；

图7是本发明一实施例的基于采样逼近法的8阶隐式交错网格有限差分方法获得的1.5s时刻的波场快照示意图；

图8是本发明一实施例的基于采样逼近法的16阶隐式交错网格有限差分方法获得的1.5s时刻的波场快照示意图；

图9是本发明一实施例的基于泰勒级数展开法的隐式交错网格有限差分方法和基于采样逼近法的隐式交错网格有限差分方法的计算时间对比图；

图10是本发明一实施例的隐式交错网格有限差分弹性波数值模拟装置的结构示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

隐式交错网格有限差分弹性波数值模拟的精度和效率一直是地震勘探领域不懈追求的目标。为了提高隐式交错网格有限差分弹性波数值模拟方法的数值模拟精度，减小数值误差，本发明提出了一种隐式交错网格有限差分弹性波数值模拟方法及装置。下面通过实施例进行介绍。

实施例一

本实施例提供了一种隐式交错网格有限差分弹性波数值模拟方法，图4是本发明一实施例的隐式交错网格有限差分弹性波数值模拟方法的流程图，如图4所示，该方法包括以下步骤(步骤S401-步骤S405)：

步骤S401，基于隐式交错网格有限差分格式(即上述公式2)和平面波理论公式(即上述公式5)，得到频散关系表达式(即上述公式8)：

其中，β＝kh/2，k是波数，h是网格大小，b是常量，M是阶数，a_m是隐式交错网格有限差分系数，m＝1,2,……,M。

步骤S402，对频散关系表达式实行采样逼近法，得到线性代数方程组。

具体地：在频散关系表达式(即上述公式8)中的β的取值范围内均匀地取(M+1)个离散点：β(1),β(2),...β(M+1)，并设置在上述(M+1)个离散点上，强制设定频散关系表达式左右两边相等，得到线性代数方程组。

步骤S403，求解线性代数方程组，得到隐式交错网格有限差分系数a_m和常量b。

基于采样逼近法得到的线性代数方程组，其求解方法至少包括以下两种：通过高斯Gauss消元法或者LU分解法求解线性代数方程组，得到隐式交错网格有限差分系数a_m和常量b。本实施例并不限定具体采样何种方式求解线性代数方程组，只要能够准确生成常量b，以及与现有技术中采用泰勒级数展开法求得的隐式交错网格有限差分系数相比更为优化的隐式交错网格有限差分系数a_m即可。

步骤S404，基于隐式交错网格有限差分系数a_m、常量b和隐式交错网格有限差分格式(即上述公式2)，求解弹性波动方程(即上述公式1)，生成求解结果。

步骤S405，根据求解结果生成波场图，并根据波场图进行波场分析。

传统的隐式交错网格有限差分系数由泰勒级数展开法获得，本实施例提出了一种新的隐式交错网格有限差分弹性波数值模拟方法，将采样逼近法运用到隐式交错网格有限差分格式的系数求取中，得到优化的隐式交错网格有限差分系数a_m和常量b，以减小空间导数项的计算误差，进而提高弹性波动方程的求解结果的精度。在不增加计算时长的同时，提高数值模拟精度，减小数值误差。

在上述步骤S402中，均匀取样的取值范围可以是β的整个取值范围(β∈[0,π/2])，出于增强某一个特定频段内的数值模拟精度的目的，也可以是上述整个取值范围[0,π/2]内的任一区间，例如[0,π/4]，本发明不以此为限，具体区间的选择根据实际需求确定。由于β＝kh/2，k是波数，在k＝Nyquist波数时，β的值为上述取值范围的最大值：π/2。

在上述步骤S404中，对于基于隐式交错网格有限差分系数a_m、常量b和隐式交错网格有限差分格式求解弹性波动方程的具体操作方法，本实施例提供了一种优选实施方式：

第一步，将隐式交错网格有限差分系数a_m和常量b代入隐式交错网格有限差分格式，求解弹性波动方程中的空间偏导数；

该步骤具体包括：将隐式交错网格有限差分系数a_m和常量b代入隐式交错网格有限差分格式，将隐式交错网格有限差分格式运用到每一个网格点，联立所有网格点上的差分方程(该隐式交错网格有限差分方程是经上述隐式交错网格有限差分格式转换得到)，得到三对角矩阵方程；求解三对角矩阵方程，得到弹性波动方程中的空间偏导数。

第二步，将利用有限差分法计算得到的空间偏导数代入弹性波动方程，经过时间迭代，生成求解结果；

该步骤具体包括：将利用有限差分法计算得到的空间偏导数代入弹性波动方程，通过对时间偏导数采用二阶有限差分，对时间进行迭代，得到任一时刻空间任意网格点的数值，以生成求解结果。

通过上述优选实施方式，将采样逼近法运用到隐式交错网格有限差分格式的系数求取中，由于将精度更高的隐式交错网格有限差分系数a_m和常量b代入隐式交错网格有限差分格式，因此得到相较于现有技术更为优化的隐式交错网格有限差分格式，运用优化的隐式交错网格有限差分格式进行弹性波动方程的求解和数值模拟，可以达到在不增加计算时长的同时，提高数值模拟精度，减小数值误差的效果。

需要说明的是，在地震勘探领域，数值模拟是指：对特定的地质、地球物理问题作适当的简化，从而形成数学模型，通过数值方法求解该模型的数学物理方程，模拟研究地震波传播。对于步骤S405中的根据求解结果所生成的波场图以及波场分析操作，与现有技术中的具体操作手段相同，在此不做赘述。只是该求解结果相对于现有技术中的求解结果而言，误差更小，精度更高，因此能够更逼真的模拟存图，更有利于提高地质构造模型模拟分析的准确性。

实施例二

本实施例对隐式交错网格有限差分弹性波数值模拟方法的完整流程进行介绍。

在均匀各向同性介质中，二维弹性波动方程(公式1)表示为：

其中，t表示时间，x和z表示空间坐标，(v_x,v_z)是速度矢量，(τ_xx,τ_zz,τ_xz)是应力张量，ρ是密度，λ和μ是Lame常数。弹性波动方程(公式1)中的时间导数项(即公式1中等号左边各项)通过二阶差分格式求解，而为了提高求解精度，空间导数项(即公式(1)中等号右边各项)通过下述隐式交错网格有限差分格式(公式2)求解。

函数p(x)的空间一阶导数展开为2M阶的隐式交错网格有限差分格式，如下：

其中，代表一阶差分算子，代表二阶差分算子，b是任意一个常量，h是网格大小，a_m(m＝1,2,...,M)是隐式交错网格的差分系数，M是阶数，m＝1,2,……,M。

令

将公式3代入公式2中，得隐式交错网格有限差分方程：

由平面波理论公式(公式5)：

p(x)＝p₀e^ikx， (5)

其中，p₀是一个常数，k是波数，e是自然对数的底数。

将公式5代入公式4，得：

令

β＝kh/2， (7)

则公式6可以化简为：

其中，β∈[0,π/2]，公式8是频散关系表达式。

公式8表明，在β的取值范围内，如果频散关系表达式的左右两边完全相等，那么隐式交错网格有限差分不会引入误差。有限差分法引入的误差是不可避免的，假设β∈[0,u]，u是一个常数，且u∈(0,π/2]；为了在全局范围内减小求解空间一阶偏导数的误差，在β的取值范围内均匀地取(M+1)个离散的点β(1),β(2),...β(M+1)，强制在这些离散点上频散关系表达式的左右两边相等，得到线性方程组，如下：

采用高斯Gauss消元法或者LU分解法求解公式13，即得到常量b，以及相比于现有技术中采用泰勒级数展开法获得的隐式交错网格差分系数更为优化的隐式交错网格差分系数a_m。

上述在β的取值范围内均匀采样，并强制在这些采样离散点上使频散关系表达式的左右两边相等的方法称作采样逼近法。根据采样逼近法得到的隐式交错网格有限差分系数a_m和常量b列于表1。

表1

将表1中所示的隐式交错网格差分系数a_m和常量b代入隐式交错网格有限差分格式(公式2)，将隐式交错网格有限差分格式运用到每一个网格点，联立所有网格点上的隐式交错网格有限差分方程，便可以得到一个三对角矩阵方程。求解此三对角矩阵方程，便可以求得在某一时刻的波场空间偏导数。

弹性波动方程(公式1)中的时间导数项(即时间偏导数)通过二阶差分格式求解，空间导数项(即空间偏导数)就通过本实施例提出的基于采样逼近法的隐式交错网格有限差分格式求解。利用弹性波动方程(公式1)进行时间迭代，便可以得到传播时间内任意时刻的波场值，进而完成弹性波数值模拟，进行弹性波场的研究。

弹性波数值模拟的核心工作就是求解弹性波动方程。传统的弹性波动方程的空间导数项在采用隐式交错网格有限差分格式进行求解时，其差分系数通常采用泰勒级数展开法得到，这样得到的差分系数仅在一个较小的频带范围内取得较高的精度，当频带范围较宽、模型复杂程度较高时，就会出现较大的数值误差。考虑到此，本实施例采用采样逼近法计算隐式交错网格差分系数，从而提出一种具有更高精度的隐式交错网格有限差分弹性波数值模拟方法。

图5是本发明一实施例的采用采样逼近法的隐式交错网格有限差分弹性波数值模拟方法的流程图，如图5所示，该流程包括以下步骤(步骤S501-步骤S505)：

步骤S501，基于隐式交错网格有限差分格式和平面波理论公式，得到频散关系表达式。

步骤S502，对频散关系表达式实行采样逼近法。

具体包括：在频散关系表达式中β的取值范围内对其进行离散，均匀采样；在离散点上，强制频散关系表达式左右两边相等，构造方程，得到关于隐式交错网格差分系数的线性代数方程组。

步骤S503，采用Gauss消元法或者LU分解法求解线性代数方程组，得到隐式交错网格有限差分系数a_m和常量b。

步骤S504，将隐式交错网格有限差分系数a_m和常量b代入隐式交错网格有限差分格式，求解弹性波动方程中的空间偏导数，进行时间迭代，进一步得到弹性波动方程的数值解。

步骤S505，根据求解结果进行弹性波场正演模拟(例如：生成波场图)，对模拟结果进行波场分析。

图6是本发明一实施例的基于采样逼近法获得的差分系数的数值频散曲线示意图，如图6所示，基于采样逼近法获得的隐式交错网格差分系数不仅在低频处获得较高的精度，与传统的泰勒级数展开法比较，其显著地提高了高频处的精度，拓宽了有效频带，减小了数值误差。

采用基于采样逼近法的隐式交错网格有限差分方法进行数值模拟，假设同样采用8阶和16阶格式进行计算，绘制1.5s时刻x分量和z分量的波场快照。图7是本发明一实施例的基于采样逼近法的8阶隐式交错网格有限差分方法获得的1.5s时刻的波场快照示意图，图8是本发明一实施例的基于采样逼近法的16阶隐式交错网格有限差分方法获得的1.5s时刻的波场快照示意图，如图7和图8所示，8阶隐式交错网格有限差分方法的模拟结果的频散现象并不明显，与基于泰勒级数展开法的隐式交错网格有限差分方法相比，其极大的提高了求解精度，减小了数值解误差。同时，基于采样逼近法的8阶隐式交错网格有限差分方法的模拟效果可与基于泰勒级数展开法的16阶隐式交错网格有限差分方法获得相当的模拟效果。而在基于采样逼近法的16阶隐式交错网格有限差分方法的模拟效果中，已基本看不出频散现象。

图9是本发明一实施例的基于泰勒级数展开法的隐式交错网格有限差分方法和基于采样逼近法的隐式交错网格有限差分方法的计算时间对比图，如图9所示，本实施例提供的隐式交错网格有限差分方法并不增加计算时长，具有相同阶数的两种方法具有接近的计算时间。

综上，本实施例提出了一种具有更高模拟精度的基于采样逼近法的隐式交错网格有限差分方法。较之于传统的基于泰勒级数展开法的隐式交错网格有限差分方法，当离散参数和算子长度相同时，降低了计算误差，提高了计算精度；当约束相同的计算误差，本实施例的方法可采用较少的差分算子，就可减少计算时长，提高计算效率。

实施例三

基于同一发明构思，本发明实施例中还提供了一种隐式交错网格有限差分弹性波数值模拟装置，可以用于实现上述实施例一所描述的方法，如下面的实施例所述。由于隐式交错网格有限差分弹性波数值模拟装置解决问题的原理与隐式交错网格有限差分弹性波数值模拟方法相似，因此隐式交错网格有限差分弹性波数值模拟装置的实施可以参见隐式交错网格有限差分弹性波数值模拟方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

图10是本发明一实施例的隐式交错网格有限差分弹性波数值模拟装置的结构示意图，如图10所示，该系统包括：频散关系计算模块10、方程组构造模块20、系数计算模块30、方程求解模块40、模拟分析模块50，下面对该结构进行具体说明。

频散关系计算模块10，用于基于隐式交错网格有限差分格式和平面波理论公式，得到频散关系表达式：其中，β＝kh/2，k是波数，h是网格大小，b是常量，M是阶数，a_m是隐式交错网格有限差分系数，m＝1,2,……,M；

方程组构造模块20，连接至频散关系计算模块10，用于在上述β的取值范围内均匀地取(M+1)个离散点：β(1),β(2),...β(M+1)，并设置在离散点上上述频散关系表达式左右两边相等，得到线性代数方程组；

系数计算模块30，连接至方程组计算模块20，用于求解上述线性代数方程组，得到隐式交错网格有限差分系数a_m和常量b；

方程求解模块40，连接至系数计算模块30，用于基于隐式交错网格有限差分系数a_m、常量b和隐式交错网格有限差分格式，求解弹性波动方程，生成求解结果；

模拟分析模块50，连接至方程求解模块40，用于根据上述求解结果生成波场图，并根据波场图进行波场分析。

本实施例提出了一种隐式交错网格有限差分弹性波数值模拟装置，将采样逼近法运用到隐式交错网格有限差分格式的系数求取中，得到优化的隐式交错网格有限差分系数a_m和常量b，以减小空间导数项的计算误差，进而提高弹性波动方程的求解结果的精度。在不增加计算时长的同时，提高数值模拟精度，减小数值误差。

上述方程组计算模块20在确定β的取值范围时，可以确定为β的整个取值范围(β∈[0,π/2])，出于增强某一个特定频段内的数值模拟精度的目的，也可以确定为上述整个取值范围[0,π/2]内的任一区间，具体区间的选择根据实际需求确定。由于β＝kh/2，k是波数，在k＝Nyquist波数时，β的值为上述取值范围的最大值：π/2。

上述系数计算模块30求解线性代数方程组时，至少可以采用下述两种方法：高斯Gauss消元法或者LU分解法。基于此，本实施例提供了一种优选实施方式，即：上述系数计算模块30，还用于通过高斯Gauss消元法或者LU分解法求解线性代数方程组，得到隐式交错网格有限差分系数a_m和常量b。

本实施例并不限定具体采样何种方式求解线性代数方程组，只要能够准确生成常量b，以及与现有技术中采用泰勒级数展开法求得的隐式交错网格有限差分系数相比更为优化的隐式交错网格有限差分系数a_m即可。

对于方程求解模块40基于隐式交错网格有限差分系数a_m、常量b和隐式交错网格有限差分格式，求解弹性波动方程的具体操作，本实施例提供了一种优选实施方式，即上述方程求解模块40包括：空间偏导数计算单元，用于将隐式交错网格有限差分系数a_m和常量b代入隐式交错网格有限差分格式，求解弹性波动方程中的空间偏导数；时间迭代单元，用于将利用有限差分法计算得到的空间偏导数代入弹性波动方程，经过时间迭代，生成求解结果。

具体地，上述空间偏导数计算单元，还用于将隐式交错网格有限差分系数a_m和常量b代入隐式交错网格有限差分格式，将隐式交错网格有限差分格式运用到每一个网格点，联立所有网格点上的隐式交错网格有限差分方程，得到三对角矩阵方程；求解上述三对角矩阵方程，得到弹性波动方程中的空间偏导数。

具体地，上述时间迭代单元，还用于将利用有限差分法计算得到的空间偏导数代入弹性波动方程，通过对时间偏导数采用二阶有限差分，对时间进行迭代，得到任一时刻空间任意网格点的数值，以生成求解结果。

通过上述优选实施方式，将采样逼近法运用到隐式交错网格有限差分格式的系数求取中，由于将精度更高的隐式交错网格有限差分系数a_m和常量b代入隐式交错网格有限差分格式，因此得到相较于现有技术更为优化的隐式交错网格有限差分格式，运用优化的隐式交错网格有限差分格式进行弹性波动方程的求解和数值模拟，可以达到在不增加计算时长的同时，提高数值模拟精度，减小数值误差的效果。

对于模拟分析模块50中根据求解结果所生成的波场图以及波场分析操作，与现有技术的具体操作手段相同，在此不做赘述。只是该求解结果相对于现有技术中的求解结果而言，误差更小，精度更高，因此能够更逼真的模拟成像，更有利于提高地质构造模拟分析的准确性。

当然，上述模块划分只是一种示意划分，本发明并不局限于此。该装置还可以仅包括：计算模块和成像分析模块，计算模块执行与计算相关的功能，成像分析模块执行与模拟成像、模拟分析相关的功能，只要能实现本发明的目的的模块划分，均应属于本发明的保护范围。

从以上的描述中可知，现有技术中为了达到数值模拟的精度和效率要求，采用隐式交错网格有限差分方法求解弹性波动方程。首先，有限差分法，尤其是交错网格有限差分，具有灵活性，易实施，拥有较高的精度和稳定性；同时，交错网格有限差分的隐格式在几乎不增加计算量的同时，大幅度的提高求解精度，提高计算效率。

在采用隐式交错网格有限差分方法求解弹性波动方程时，传统的隐式交错网格差分系数通过泰勒级数展开法获得。这样得到的隐式交错网格有限差分系数，当用其求解弹性波动方程，进而进行数值模拟时，仅在中低频具有较高的模拟精度，而高频领域精度较差，使数值模拟存在严重的误差。针对这一问题，本发明将采样逼近法运用到隐式交错网格有限差分系数的求取中，达到了在不增加计算时长的同时，提高数值模拟精度，减小数值误差的效果。

在本说明书的描述中，参考术语“实施例一”、“本实施例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (14)

1.一种隐式交错网格有限差分弹性波数值模拟方法，其特征在于，所述方法包括：

基于隐式交错网格有限差分格式和平面波理论公式，得到频散关系表达式：其中，β＝kh/2，k是波数，h是网格大小，b是常量，M是阶数，a_m是隐式交错网格有限差分系数，m＝1,2,……,M；

在所述β的取值范围内均匀地取M+1个离散点：β(1),β(2),…β(M+1)，并设置在所述离散点上所述频散关系表达式左右两边相等，得到线性代数方程组；

求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b；

基于所述隐式交错网格有限差分系数a_m、所述常量b和所述隐式交错网格有限差分格式，求解弹性波动方程，生成求解结果；

根据所述求解结果生成波场图，并根据所述波场图进行波场分析。

2.根据权利要求1所述的隐式交错网格有限差分弹性波数值模拟方法，其特征在于，所述β的取值范围是：[0,π/2]，或者[0,π/2]中的任一区间。

3.根据权利要求1所述的隐式交错网格有限差分弹性波数值模拟方法，其特征在于，所述线性代数方程组是：

$\begin{array}{c}\left[\begin{array}{ccccc}\sin \[(2\×1-1)\mathrm{\β}\left(1\right)\]& \sin \[(2\×2-1)\mathrm{\β}\left(1\right)\]& \mathrm{...}& \sin \[(2\×M-1)\mathrm{\β}\left(1\right)\]& 2\mathrm{\β}\left(1\right)-2\mathrm{\β}\left(1\right)\cos \left(2\mathrm{\β}\right(1\left)\right)\\ \sin \[(2\×1-1)\mathrm{\β}\left(2\right)\]& \sin \[(2\×2-1)\mathrm{\β}\left(2\right)\]& \mathrm{...}& \sin \[(2\×M-1)\mathrm{\β}\left(2\right)\]& 2\mathrm{\β}\left(2\right)-2\mathrm{\β}\left(2\right)\cos \left(2\mathrm{\β}\right(2\left)\right)\\ .& & & & .\\ .& & & & .\\ .& & & & .\\ \sin \[(2\×1-1)\mathrm{\β}(M+1)\]& \sin \[(2\×2-1)\mathrm{\β}(M+1)\]& \mathrm{...}& \sin \[(2\×M-1)\mathrm{\β}(M+1)\]& 2\mathrm{\β}(M+1)-2\mathrm{\β}(M+1)\cos \left(2\mathrm{\β}\right(M+1\left)\right)\end{array}\right]\\ \×\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ b\end{array}\right]=\left[\begin{array}{c}\mathrm{\β}\left(1\right)\\ \mathrm{\β}\left(2\right)\\ .\\ .\\ .\\ \mathrm{\β}(M+1)\end{array}\right].\end{array}$

4.根据权利要求1所述的隐式交错网格有限差分弹性波数值模拟方法，其特征在于，求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b，包括：

通过高斯消元法或者LU分解法求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b。

5.根据权利要求1所述的隐式交错网格有限差分弹性波数值模拟方法，其特征在于，基于所述隐式交错网格有限差分系数a_m、所述常量b和所述隐式交错网格有限差分格式，求解弹性波动方程，生成求解结果，包括：

将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，求解所述弹性波动方程中的空间偏导数；

将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，经过时间迭代，生成所述求解结果。

6.根据权利要求5所述的隐式交错网格有限差分弹性波数值模拟方法，其特征在于，将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，求解所述弹性波动方程中的空间偏导数，包括：

将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，将所述隐式交错网格有限差分格式运用到每一个网格点，联立所有网格点上的隐式交错网格有限差分方程，得到三对角矩阵方程；

求解所述三对角矩阵方程，得到所述弹性波动方程中的空间偏导数；其中，所述隐式交错网格有限差分方程是经所述隐式交错网格有限差分格式转换得到。

7.根据权利要求5所述的隐式交错网格有限差分弹性波数值模拟方法，其特征在于，将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，经过时间迭代，生成所述求解结果，包括：

将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，通过对时间偏导数采用二阶有限差分，对时间进行迭代，得到任一时刻空间任意网格点的数值，以生成所述求解结果。

8.一种隐式交错网格有限差分弹性波数值模拟装置，其特征在于，所述装置包括：

频散关系计算模块，用于基于隐式交错网格有限差分格式和平面波理论公式，得到频散关系表达式：其中，β＝kh/2，k是波数，h是网格大小，b是常量，M是阶数，a_m是隐式交错网格有限差分系数，m＝1,2,……,M；

方程组构造模块，用于在所述β的取值范围内均匀地取M+1个离散点：β(1),β(2),…β(M+1)，并设置在所述离散点上所述频散关系表达式左右两边相等，得到线性代数方程组；

系数计算模块，用于求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b；

方程求解模块，用于基于所述隐式交错网格有限差分系数a_m、所述常量b和所述隐式交错网格有限差分格式，求解弹性波动方程，生成求解结果；

模拟分析模块，用于根据所述求解结果生成波场图，并根据所述波场图进行波场分析。

9.根据权利要求8所述的隐式交错网格有限差分弹性波数值模拟装置，其特征在于，所述β的取值范围是：[0,π/2]，或者[0,π/2]中的任一区间。

10.根据权利要求8所述的隐式交错网格有限差分弹性波数值模拟装置，其特征在于，所述线性代数方程组是：

$\begin{array}{c}\left[\begin{array}{ccccc}\sin \[(2\×1-1)\mathrm{\β}\left(1\right)\]& \sin \[(2\×2-1)\mathrm{\β}\left(1\right)\]& \mathrm{...}& \sin \[(2\×M-1)\mathrm{\β}\left(1\right)\]& 2\mathrm{\β}\left(1\right)-2\mathrm{\β}\left(1\right)\cos \left(2\mathrm{\β}\right(1\left)\right)\\ \sin \[(2\×1-1)\mathrm{\β}\left(2\right)\]& \sin \[(2\×2-1)\mathrm{\β}\left(2\right)\]& \mathrm{...}& \sin \[(2\×M-1)\mathrm{\β}\left(2\right)\]& 2\mathrm{\β}\left(2\right)-2\mathrm{\β}\left(2\right)\cos \left(2\mathrm{\β}\right(2\left)\right)\\ .& & & & .\\ .& & & & .\\ .& & & & .\\ \sin \[(2\×1-1)\mathrm{\β}(M+1)\]& \sin \[(2\×2-1)\mathrm{\β}(M+1)\]& \mathrm{...}& \sin \[(2\×M-1)\mathrm{\β}(M+1)\]& 2\mathrm{\β}(M+1)-2\mathrm{\β}(M+1)\cos \left(2\mathrm{\β}\right(M+1\left)\right)\end{array}\right]\\ \×\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ b\end{array}\right]=\left[\begin{array}{c}\mathrm{\β}\left(1\right)\\ \mathrm{\β}\left(2\right)\\ .\\ .\\ .\\ \mathrm{\β}(M+1)\end{array}\right].\end{array}$

11.根据权利要求8所述的隐式交错网格有限差分弹性波数值模拟装置，其特征在于，所述系数计算模块，还用于通过高斯消元法或者LU分解法求解所述线性代数方程组，得到所述隐式交错网格有限差分系数a_m和所述常量b。

12.根据权利要求8所述的隐式交错网格有限差分弹性波数值模拟装置，其特征在于，所述方程求解模块，包括：

空间偏导数计算单元，用于将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，求解所述弹性波动方程中的空间偏导数；

时间迭代单元，用于将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，经过时间迭代，生成所述求解结果。

13.根据权利要求12所述的隐式交错网格有限差分弹性波数值模拟装置，其特征在于，所述空间偏导数计算单元，具体用于将所述隐式交错网格有限差分系数a_m和所述常量b代入所述隐式交错网格有限差分格式，将所述隐式交错网格有限差分格式运用到每一个网格点，联立所有网格点上的隐式交错网格有限差分方程，得到三对角矩阵方程；求解所述三对角矩阵方程，得到所述弹性波动方程中的空间偏导数；其中，所述隐式交错网格有限差分方程是经所述隐式交错网格有限差分格式转换得到。

14.根据权利要求12所述的隐式交错网格有限差分弹性波数值模拟装置，其特征在于，所述时间迭代单元具体用于将利用有限差分法计算得到的所述空间偏导数代入所述弹性波动方程，通过对时间偏导数采用二阶有限差分，对时间进行迭代，得到任一时刻空间任意网格点的数值，以生成所述求解结果。