有限差分数值模拟方法、系统、介质、计算机设备及应用

摘要

本发明属于地球物理勘探数据处理技术领域，公开了一种有限差分数值模拟方法、系统、介质、计算机设备及应用，通过时间任意偶数阶精度有限差分逼近时间偏导数、空间隐式交错网格有限差分逼近空间偏导数、纵横波波场分离实现弹性P‑SV波高精度数值模拟。本发明的弹性P‑SV波高精度有限差分数值模拟方法，采用“十字+菱形”交错网格有限差分模板逼近弹性波方程中的时间偏导数；采样隐式交错网格有限差分法逼近弹性波方程中的空间偏导数；推导了相应的频散关系；基于Taylor级数展开求取差分系数；基于纵横波分离实现了P‑SV波高精度模拟。本发明可以大幅提高P‑SV波的模拟精度，提高弹性波逆时偏移和全波形反演的精度。

说明书

技术领域

本发明属于地球物理勘探数据处理技术领域，尤其涉及一种有限差分数值模拟方法、系统、介质、计算机设备及应用。

背景技术

目前：波动方程数值模拟是了解地震波在地下的传播规律、帮助解释观测数据的有效手段，也是逆时偏移成像和全波形反演的基本单元。到目前为止，波动方程数值模拟已被广泛应用于地球物理勘探领域的各个环节。常用的偏微分方程数值解法包括有限元素法、有限差分法、伪谱法等。有限差分法容易实现、计算效率高，是最常用波动方程数值模拟算法。地球物理学家对有限差分法进行过大量研究，先后提出显式有限差分和隐式有限差分；规则网格有限差分、交错网格有限差分和旋转交错网格有限差分；基于Taylor级数展开的有限差分和基于最优化的有限差分；空间域有限差分和时空域有限差分等。上述有限差分法大多是基于常规“十字”型差分模版设计的，时间精度始终只有二阶。当模拟时间较长或时间步长较大时，这些差分方法都会产生较强的时间频散，模拟精度不够。

与常规方法相比，时间高阶有限差分法可以更好地压制时间频散，且具有更好的稳定性。Lax-Wendroff方法是最早出现的时间高阶有限差分法，该方法通过空间导数来代替高阶时间导数来提高时间精度。Lax-Wendroff法需要存储多个时刻的波场。随着时间精度的增加，存储量和计算量急剧增大，现有的计算设备很难承受，而且还会出现数值不稳定问题，但Lax-Wendroff法很难推广到任意偶数阶精度、且计算量较大。研究发现有限差分法的时间模拟精度和所采用的差分模版种类密切相关。随后出现了基于“菱形”模版的声波方程规则网格有限差分法，该方法可以同时获得任意偶数阶空间和时间精度。地球物理工作者还将常规“十字型”模板和“菱形”模版相结合来获得空间(2M)阶、时间(2N)阶精度。现有的时间高阶有限差分法多数只适用于声波方程，无法直接用于弹性波方程数值模拟，且空间模拟精度仍然不够，现有的时间高阶有限差分法主要提高时间精度，空间精度仍需进一步提高。

通过上述分析，现有技术存在的问题及缺陷为：

(1)现有时间高阶有限差分法很难推广到任意偶数阶精度、且计算量较大。

(2)现有的时间高阶有限差分法多数只适用于声波方程，无法直接用于弹性波方程数值模拟，且空间模拟精度仍然不够。

解决以上问题及缺陷的难度为：弹性介质中，涉及纵波、横波和转换波等波模式。采用基于纵波或横波频散关系的差分系数进行求解只能分别得到精确的纵波或横波波场，都无法得到高精度的P-SV转换波。提高空间精度会一定程度降低时间精度和稳定性。

解决以上问题及缺陷的意义为：研发弹性P-SV波有限差分数值模拟方法可以得到大幅提高纵波、横波和P-SV转换波的模拟精度，进而改善弹性波逆时偏移和全波形反演的成像和反演效果。

发明内容

针对现有技术存在的问题，本发明提供了一种有限差分数值模拟方法、系统、介质、计算机设备及应用。

本发明是这样实现的，一种有限差分数值模拟方法，所述有限差分数值模拟方法包括：

采用十字和菱形交错网格有限差分模板逼近时间偏导数；

采用隐式交错网格有限差分法逼近空间偏导数；

频散关系推导；

差分系数求取；

得到的时间高阶差分系数和隐式空间差分系数分别代入公式进行弹性波方程数值模拟，实现基于纵横波分离的P-SV波高精度模拟。

进一步，所述十字和菱形交错网格有限差分模板逼近时间偏导数方法包括：

二维弹性波速度应力方程表示为：

其中，(v

其中，

其中d

进一步，所述隐式交错网格有限差分法逼近空间偏导数方法包括：隐式交错网格有限差分逼近空间偏导数的公式为：

其中c

隐式交错网格有限差分中，空间偏导数通过求解三对角线性方程组得到。

进一步，所述频散关系推导方法包括：

平面波理论下，波场表示公式为：

其中，

其中，r

公式

进一步，所述差分系数求取方法包括：

将公式

C+D＝E；

其中，

比较方程C+D＝E两端(k

a

比较方程C+D＝E两端

其中，int(x)为取整函数，通过公式：

a

隐式空间交错网格有限差分系数求取公式为：

进一步，所述基于纵横波分离的P-SV波高精度模拟方法包括：将第四步中得到的时间高阶差分系数和隐式空间差分系数分别代入公式

(1)采用纵波时间高阶差分系数和隐式空间差分系数求解弹性波方程，纵横波分离，只提取纵波场；

(2)采用横波时间高阶差分系数和隐式空间差分系数求解弹性波方程，纵横波分离，只提取横波场；

(3)将(1)中的纵波场和(2)中的横波场求和，得到最终的弹性波场；

(4)重复(1)-(3)，直到最大时刻。

本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行如下步骤：

采用十字和菱形交错网格有限差分模板逼近时间偏导数；

采用隐式交错网格有限差分法逼近空间偏导数；

频散关系推导；

差分系数求取；

得到的时间高阶差分系数和隐式空间差分系数分别代入公式进行弹性波方程数值模拟，实现基于纵横波分离的P-SV波高精度模拟。

本发明的另一目的在于提供一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如下步骤：

采用十字和菱形交错网格有限差分模板逼近时间偏导数；

采用隐式交错网格有限差分法逼近空间偏导数；

频散关系推导；

差分系数求取；

得到的时间高阶差分系数和隐式空间差分系数分别代入公式进行弹性波方程数值模拟，实现基于纵横波分离的P-SV波高精度模拟。

本发明的另一目的在于提供一种实施所述有限差分数值模拟方法的有限差分数值模拟系统，所述有限差分数值模拟系统包括：

时间偏导数获取模块，用于采用十字+菱形交错网格有限差分模板逼近时间偏导数；

空间偏导数获取模块，用于采用隐式交错网格有限差分法逼近空间偏导数；

频散关系获取模块，用于实现频散关系推导；

差分系数求取模块，用于差分系数求取；

数值模拟模块，用于基于纵横波分离实现P-SV波高精度数值模拟。

本发明的另一目的在于提供一种所述有限差分数值模拟方法在弹性介质地震波正演、逆时偏移和全波形反演中的应用。

结合上述的所有技术方案，本发明所具备的优点及积极效果为：本发明通过时间任意偶数阶精度有限差分逼近时间偏导数、空间隐式交错网格有限差分逼近空间偏导数、纵横波波场分离实现弹性P-SV波高精度数值模拟。本发明的弹性P-SV波高精度有限差分数值模拟方法，采用“十字+菱形”交错网格有限差分模板逼近弹性波方程中的时间偏导数；采样隐式交错网格有限差分法逼近弹性波方程中的空间偏导数；推导了相应的频散关系；基于Taylor级数展开求取差分系数；基于纵横波分离实现了P-SV波高精度模拟。

本发明通过采用时间任意偶数阶精度-空间隐式交错网格有限差分法实现弹性P-SV波高精度数值模拟；提高弹性P-SV波的模拟精度，进而提高弹性波逆时偏移和全波形反演的精度。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的有限差分数值模拟方法流程图。

图2是本发明实施例提供的有限差分数值模拟系统的结构示意图；

图2中：1、时间偏导数获取模块；2、空间偏导数获取模块；3、频散关系获取模块；4、差分系数求取模块；5、数值模拟模块。

图3是本发明实施例提供的有限差分数值模拟方法实现流程图。

图4是本发明实施例提供的“十字+菱形”交错网格有限差分模板示意图。

图5是本发明实施例提供的双层弹性介质模型示意图。

图6(a)是本发明实施例提供的不同方法得到的水平分量波场快照；I和II分别表示常规方法和本发明中提出的方法。

图6(b)是本发明实施例提供的不同方法得到的垂直分量波场快照；I和II分别表示常规方法和本发明中提出的方法。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对现有技术存在的问题，本发明提供了一种有限差分数值模拟方法、系统、介质、计算机设备及应用，下面结合附图对本发明作详细的描述。

如图1所示，本发明提供的有限差分数值模拟方法包括以下步骤：

S101：采用十字+菱形”交错网格有限差分模板逼近时间偏导数；

S102：采用隐式交错网格有限差分法逼近空间偏导数；

S103：频散关系推导；

S104：差分系数求取；

S105：得到的时间高阶差分系数和隐式空间差分系数分别代入公式进行弹性波方程数值模拟，实现基于纵横波分离的P-SV波高精度模拟。

本发明提供的有限差分数值模拟方法业内的普通技术人员还可以采用其他的步骤实施，图1的本发明提供的有限差分数值模拟方法仅仅是一个具体实施例而已。

如图2所示，本发明提供的有限差分数值模拟系统包括：

时间偏导数获取模块1，用于采用十字+菱形交错网格有限差分模板逼近时间偏导数；

空间偏导数获取模块2，用于采用隐式交错网格有限差分法逼近空间偏导数；

频散关系获取模块3，用于实现频散关系推导；

差分系数求取模块4，用于差分系数求取；

数值模拟模块5，用于基于纵横波分离实现P-SV波高精度数值模拟。

下面结合附图对本发明的技术方案作进一步的描述。

如图3所示，本发明提供的有限差分数值模拟具体包括以下步骤：

第一步，“十字+菱形”交错网格有限差分模板逼近时间偏导数；

第二步，隐式交错网格有限差分法逼近空间偏导数；

第三步，频散关系推导；

第四步，差分系数求取；

第五步，基于纵横波分离的P-SV波高精度数值模拟。

在本发明中，第一步中“十字+菱形”交错网格有限差分模板逼近时间偏导数方法如下：

二维弹性波速度应力方程可以表示为：

其中，(v

其中，

其中d

在本发明中，第二步中隐式交错网格有限差分法逼近空间偏导数方法如下：隐式交错网格有限差分逼近空间偏导数的公式为(以v

其中c

隐式交错网格有限差分中，空间偏导数通过求解上述三对角线性方程组得到。

在本发明中，第三步中频散关系推导方法如下：

平面波理论假设下，波场可以表示为：

其中，

其中，r

公式(7)和(8)即为纵波和横波的时间频散关系。

在本发明中，第四步中差分系数求取方法如下：

将公式(7)中的三角函数进行Taylor级数展开化简可得：

C+D＝E (10)

其中，

比较方程(10)两端(k

a

比较方程(10)两端

其中，int(x)为取整函数。通过公式(14)和(15)可得到基于Taylor级数展开的纵波时间高阶差分系数，将v

隐式空间交错网格有限差分系数求取公式为：

在本发明中，第五步中基于纵横波分离的P-SV波高精度模拟方法如下：

将第四步中得到的时间高阶差分系数和隐式空间差分系数分别代入公式(2)和(4)可进行弹性波方程数值模拟。但弹性介质中，涉及纵波、横波和转换波等多种波模式。采用纵波时间高阶差分系数进行求解只能得到精确的纵波场；采用横波时间高阶差分系数进行求解只能得到精确的横波场；采用纵波或横波差分系数求解都得不到高精度的P-SV转换波。为实现P-SV波高精度模拟，引入纵横波分离(本发明采用波数域分离方法)。模拟步骤变为：

(1)采用纵波时间高阶差分系数和隐式空间差分系数求解弹性波方程，纵横波分离，只提取纵波场；

(2)采用横波时间高阶差分系数和隐式空间差分系数求解弹性波方程，纵横波分离，只提取横波场；

(3)将(1)中的纵波场和(2)中的横波场求和，得到最终的弹性波场。

(4)重复上面的步骤，直到最大时刻。

下面结合实验对本发明的技术效果作详细的描述。

本实验采用如图5所示的双层弹性介质模型验证本发明中提出的方法。模型计算区域有201*201个网格点，空间间隔为10.0m，时间步长为1.5ms，最大记录时间为1.5s。震源为40Hz的雷克子波，位于(1000m，800m)，201个检波点均匀分布于地表。差分算子长度M＝8，N＝4。图6(a)和图6(b)给出不同有限差分方法得到的波场快照。为了便于比较，将常规方法和本发明的新方法的结果拼接显示(分别位于左右两侧)。由图可知，常规方法结果中出现了明显的时间和空间频散，本发明中提出的方法有效压制了这些数值频散。因此，本发明提出的弹性P-SV波高精度有限差分数值模拟方法方法可以得到精确的纵波、横波和P-SV转换波。

应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、CD或DVD-ROM的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。