一种基于规则网格的声波波动方程正演模拟方法及装置与流程

本发明涉及地震波场数值模拟技术领域，尤其涉及一种基于规则网格的声波波动方程正演模拟方法及装置。

背景技术：

地震波成像和反演需要一种高效率和高精度的算法来模拟波的正向和反向传播，例如当前热门的逆时偏移和全波形反演技术。因此研究高精度和高效率的数值模拟算法很有必要。

有限差分法是一种灵活且简便的数值算法，已被广泛应用在波动方程的数值求解之中。对于有限差分法而言，同时压制波动方程的时间和空间频散是最大的挑战之一。1986年，dablain通过运用lax-wendroff方法来提高声波波动方程的模拟精度，即将波动方程中高阶时间偏导使用空间偏导替代，但这种方法计算量将显著增加；传统的有限差分法的空间差分系数是基于空间域频散关系推导得到的，2007年，finkelstein等提出在时间空间域确定空间差分系数，该方法使得时空域频散方程在指定的若干个频率上严格成立，从而获得若干个方程，进而求解出空间差分系数，虽然获得的差分系数能够减小频散关系式，但因为不同的频率需要不同的差分系数因此难以用于实际。为了在不显著增加计算机内存的条件下提高模拟精度，liu等基于时间空间域频散关系式获得了声波波动方程的有限差分系数，相比于基于空间域频散关系获得的差分系数，在相同的离散条件下，一维模拟无条件稳定且时间空间均为阶精度，二维模拟在波场传播的8个方向时间达到阶精度，三维情况48个方向达到阶精度，而在其他的传播方向上时间精度仍为二阶；为了使波传播的各个方向空间和时间达到阶精度，2013年，liu等发展了一种菱形差分格式，但计算效率大大地降低。为了同时提高模拟精度和效率，2014年，tan等基于新差分结构发展了时间四阶、空间任意偶数阶精度的交错网格有限差分法，其差分系数是经泰勒级数展开获得的，这种新的差分格式保留了传统交错网格有限差分控制空间频散关系式的优势，同时增强了其压制时间频散关系式的能力；为了在保持时间四阶精度的同时进一步的提高空间精度，同年，tan等提出采用非线性优化方法来求得优化的差分系数，然而，优化方法需要重复的迭代计算比较耗时，而chen等则通过在最小二乘意义下使得交错网格有限差分算子与一阶波数-空间算子的波数响应误差最小化来优化差分系数。随后，在2016年，张保庆等发展了规则网格剖分的空间任意阶，时间四阶精度的有限差分法。但模拟精度和稳定性都不够高。

技术实现要素：

本发明的主要目的在于提出一种基于规则网格的声波波动方程正演模拟方法及装置，通过规则网格建立新的差分结构，求取的差分系数能够使声波在更大的波数范围内压制数值频散，进一步提高了声波波动方程的模拟精度。

为实现上述目的，本发明提供的一种基于规则网格的声波波动方程正演模拟方法，包括：

获取地震参数；

建立基于规则网格的声波波动方程；

采用时空域有限差分法计算所述声波波动方程的频散关系式；

根据所述频散关系式获取波场模拟所满足的稳定条件；

采用吸收边界条件对声波波动方程进行波场延拓，获取波场及地震记录。

可选地，所述地震参数包括：正演模拟所需的速度场文件、差分算子阶数、震源函数及其主频、正演所采用的时间与空间步长及地震记录时长，海绵吸收边界条件的参数。

可选地，所述建立基于规则网格的声波波动方程包括：

时间偏导数采用二阶差分离散；其公式为：

其中，x,z为直角坐标系的两个轴，t为时间，v为地震波的传播速度，p代表声波波场；

空间偏导数的公式为：

其中，和分别为x和z方向的差分算子，上角标(2m,4)代表空间2m阶时间4阶精度，代表离散的声波波场h代表x方向或z方向的网格离散间隔，a0,0,a1,1和am,0(m＝1,2,...,m)代表差分系数。

可选地，所述采用时空域有限差分法计算所述声波波动方程的频散关系式包括：

计算声波波动方程的离散形式；

获取时空域有限差分法与lax-wendroff方案的关系式；

采用褶积微分算子法求解所述差分系数；

利用平面波理论，得到时空域频散关系式。

可选地，所述声波波动方程的离散形式计算公式为：

作为本发明的另一方面，提供的一种基于规则网格的声波波动方程正演模拟装置，包括：

获取模块，用于获取地震参数；

建模模块，用于建立基于规则网格的声波波动方程；

差分模块，用于采用时空域有限差分法计算所述声波波动方程的频散关系式；

模拟模块，用于根据所述频散关系式获取波场模拟所满足的稳定条件；

延拓模块，用于采用吸收边界条件对声波波动方程进行波场延拓，获取波场及地震记录。

可选地，所述地震参数包括：正演模拟所需的速度场文件、差分算子阶数、震源函数及其主频、正演所采用的时间与空间步长及地震记录时长，海绵吸收边界条件的参数。

可选地，所述建立基于规则网格的声波波动方程包括：

时间偏导数采用二阶差分离散；其公式为：

其中，x,z为直角坐标系的两个轴，t为时间，v为地震波的传播速度，p代表声波波场；

空间偏导数的公式为：

其中，和分别为x和z方向的差分算子，上角标(2m,4)代表空间2m阶时间4阶精度，代表离散的声波波场h代表x方向或z方向的网格离散间隔，a0,0,a1,1和am,0(m＝1,2,...,m)代表差分系数。

可选地，所述采用时空域有限差分法计算所述声波波动方程的频散关系式包括：

计算声波波动方程的离散形式；

获取时空域有限差分法与lax-wendroff方案的关系式；

采用褶积微分算子法求解所述差分系数；

利用平面波理论，得到时空域频散关系式。

可选地，所述声波波动方程的离散形式计算公式为：

本发明提出的一种基于规则网格的声波波动方程正演模拟方法及装置，该方法包括：获取地震参数；建立基于规则网格的声波波动方程；采用时空域有限差分法计算所述声波波动方程的频散关系式；根据所述频散关系式获取波场模拟所满足的稳定条件；采用吸收边界条件对声波波动方程进行波场延拓，获取波场及地震记录；通过规则网格建立新的差分结构，求取的差分系数能够使声波在更大的波数范围内压制数值频散，进一步提高了声波波动方程的模拟精度。

附图说明

图1为本发明实施例一提供的一种基于规则网格的声波波动方程正演模拟方法的流程图；

图2为图1中步骤s30的方法流程图；

图3为本发明实施例一提供的频散误差对比示意图；

图4为本发明实施例一提供的稳定性对比图；

图5为本发明实施例一提供的均匀速度模型中0.6s时刻的波场快照图；

图6为本发明实施例一提供的复杂的marmousi速度模型；

图7为本发明实施例一提供的marmousi速度模型中的4.0s时刻的波场快照图；

图8为本发明实施例二提供的另一种基于规则网格的声波波动方程正演模拟装置的示范性结构框图。

本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

具体实施方式

应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

在后续的描述中，使用用于表示元件的诸如“模块”、“部件”或“单元”的后缀仅为了有利于本发明的说明，其本身并没有特定的意义。因此，"模块"与"部件"可以混合地使用。

实施例

如图1所示，在本实施例中，一种基于规则网格的声波波动方程正演模拟方法，包括：

s10、获取地震参数；

s20、建立基于规则网格的声波波动方程；

s30、采用时空域有限差分法计算所述声波波动方程的频散关系式；

s40、根据所述频散关系式获取波场模拟所满足的稳定条件；

s50、采用吸收边界条件对声波波动方程进行波场延拓，获取波场及地震记录。

在本实施例中，通过规则网格建立新的差分结构，求取的差分系数能够使声波在更大的波数范围内压制数值频散，进一步提高了声波波动方程的模拟精度。

在本实施例中，所述地震参数包括：正演模拟所需的速度场文件、差分算子阶数、震源函数及其主频、正演所采用的时间与空间步长及地震记录时长，海绵吸收边界条件的参数。

在本实施例中，所述建立基于规则网格的声波波动方程包括：

时间偏导数采用二阶差分离散；其公式为：

其中，x,z为直角坐标系的两个轴，t为时间，v为地震波的传播速度，p代表声波波场；

运用二阶差分对时间偏导数进行离散：

其中，δt表示时间步长，

在本实施例中，空间偏导数的公式为：

其中，和分别为x和z方向的差分算子，上角标(2m,4)代表空间2m阶时间4阶精度，代表离散的声波波场h代表x方向或z方向的网格离散间隔，a0,0,a1,1和am,0(m＝1,2,...,m)代表差分系数。

如图2所示，在本实施例中，所述步骤s30包括：

s31、计算声波波动方程的离散形式；

s32、获取时空域有限差分法与lax-wendroff方案的关系式；

s33、采用褶积微分算子法求解所述差分系数；

s34、利用平面波理论，得到时空域频散关系式。

在本实施例中，所述声波波动方程的离散形式计算公式为：

在本实施例中，新差分结构的时空域有限差分法与lax-wendroff方案的关系式为：

在本实施例中，褶积微分算子的基本原理是：将偏导数算子的傅里叶变换进行傅里叶反变换得到的；本实施例推导得到的差分系数表达式为：

其中，

d2(m)＝amw(m),d4(m)＝bmw(m),

|m|＝0,1,2...,m且0.5≤α≤1.。

对声波波动方程的离散形式计算公式的两边进行傅里叶变换，得到声波波动方程基于新差分结构的时空域有限差分法的频散关系表达式：

由差分系数表达式得到新差分结构的时空域有限差分法的频散误差表达式为：

如图3所示，为频散误差对比示意图，其中，图3(a)和图3(b)为传统有限差分法，图3(c)和图3(d)为本发明方法；图中纵坐标为数值相速度与真实速度比值，如果该比值接近于1，则频散误差小，反之频散误差大，图中横坐标表示归一化波数，算子长度2m＝16，如图3(a)所示，当r＝0.2，传统有限差分法仍然存在明显的频散误差；当r＝0.4(实际上是时间步长δt增大一倍)传统有限差分法的频散显著增加；分别对比图3(a)和图3(c)、图3(b)和图3(d)可以看出本发明方法可以在更大的波数范围内保持频散误差小，特别是在高波数情况下，表明本发明方法的模拟精度优于传统有限差分法。

由频散关系得到本实施例的有限差分法满足：

取那奎思波特波数(kx,kz)＝(π/h,π/h)代入频散误差表达式得到稳定条件为：

采用海绵吸收边界条件对边界反射进行吸收，利用求取的差分系数进行声波波场外推，得到任意时刻的波场及整个地震记录。记录波场快照切片，输出地震记录。

图4为本发明方法与传统有限差分法稳定性对比图，图中纵坐标为稳定性因子，该值越大表明方法更稳定，横坐标为差分算子长度2m，对比发现，当差分算子长度相同时，本方明方法的稳定性条件均大于传统有限差分法，这意味着本发明方法稳定性更好，可以采用更大的时间步长进行波场延拓。

图5为均匀速度模型中0.6s时刻的波场快照图。介质的地震波传播速度为1500m/s，模型网格区域为201×201点，网格的纵横空间点距为15m，震源位置在模型正中央，震源主频14.5hz，差分算子的长度2m＝16。对比图5(a)和图5(b)可知，当r由0.2增至0.4时(实际上是时间步长δt由2ms增大至4ms)，传统有限差分法的数值频散显著增强；再对比图5中(a)和(c)、(b)和(d)可知，本发明方法的模拟精度优于传统有限差分法，这与图3中的频散误差曲线分析一致，证实了本发明方法模拟精度更高。

图6为复杂的marmousi速度模型，模型网格区域为2721×701点，网格的纵横空间点距为20m，震源主频为10hz，位于模型地表的正中央，差分算子长度2m＝16。

图7为marmousi速度模型中的4.0s时刻的波场快照图，(a)为采用传统有限差分法，(b)为采用本发明方法图，从图7中黑色矩形框的放大图可以明显的看出，本发明方法提高了波场模拟精度，说明了本文发展的方法的有效性与稳健性。

在本实施例中，基于规则网格的差分结构除了包含轴向2m+1个离散点外，还包含4个非轴向离散点。与传统的差分结构相比，引入了额外的计算量。然而，新差分结构改善了传统时间2阶方法的精度，且通过稳定性分析可知，基于新差分结构的时空域有限差分法可以采用更大的时间步长。从推导出来的差分系数表达式可以看出，差分系数随着介质速度的变化而改变，本实施例的实现策略是预先计算并存储给定速度模型范围内差分系数。例如，对于给定的模型速度范围为1000m/s至5000m/s，本实施例从1000m/s开始以1m/s的增量计算并存储这些差分系数，时间4阶空间16阶差分格式只需要大约0.1mb存储这些差分系数。因此，本实施例的有限差分法和传统时间二阶方案计算内存几乎相等，本方法能够为计算量巨大的逆时偏移成像和全波形反演提供更为准确的波场信息，具有重要的实际用途。

实施例二

如图8所示，在本实施例中，一种基于规则网格的声波波动方程正演模拟装置，包括：

获取模块10，用于获取地震参数；

建模模块20，用于建立基于规则网格的声波波动方程；

差分模块30，用于采用时空域有限差分法计算所述声波波动方程的频散关系式；

模拟模块40，用于根据所述频散关系式获取波场模拟所满足的稳定条件；

延拓模块50，用于采用吸收边界条件对声波波动方程进行波场延拓，获取波场及地震记录。

在本实施例中，通过规则网格建立新的差分结构，求取的差分系数能够使声波在更大的波数范围内压制数值频散，进一步提高了声波波动方程的模拟精度。

在本实施例中，所述地震参数包括：正演模拟所需的速度场文件、差分算子阶数、震源函数及其主频、正演所采用的时间与空间步长及地震记录时长，海绵吸收边界条件的参数。

在本实施例中，所述建立基于规则网格的声波波动方程包括：

时间偏导数采用二阶差分离散；其公式为：

其中，x,z为直角坐标系的两个轴，t为时间，v为地震波的传播速度，p代表声波波场；

运用二阶差分对时间偏导数进行离散：

其中，δt表示时间步长，

空间偏导数的公式为：

其中，和分别为x和z方向的差分算子，上角标(2m,4)代表空间2m阶时间4阶精度，代表离散的声波波场h代表x方向或z方向的网格离散间隔，a0,0,a1,1和am,0(m＝1,2,...,m)代表差分系数。

在本实施例中，所述采用时空域有限差分法计算所述声波波动方程的频散关系式包括：

计算声波波动方程的离散形式；

获取时空域有限差分法与lax-wendroff方案的关系式；

采用褶积微分算子法求解所述差分系数；

利用平面波理论，得到时空域频散关系式。

在本实施例中，所述声波波动方程的离散形式计算公式为：

在本实施例中，新差分结构的时空域有限差分法与lax-wendroff方案的关系式为：

在本实施例中，褶积微分算子的基本原理是：将偏导数算子的傅里叶变换进行傅里叶反变换得到的；本实施例推导得到的差分系数表达式为：

其中，

d2(m)＝amw(m),d4(m)＝bmw(m),

|m|＝0,1,2...,m且0.5≤α≤1.

对声波波动方程的离散形式计算公式的两边进行傅里叶变换，得到声波波动方程基于新差分结构的时空域有限差分法的频散关系表达式：

由差分系数表达式得到新差分结构的时空域有限差分法的频散误差表达式为：

如图3所示，为频散误差对比示意图，其中，图3(a)和图3(b)为传统有限差分法，图3(c)和图3(d)为本发明方法；图中纵坐标为数值相速度与真实速度比值，如果该比值接近于1，则频散误差小，反之频散误差大，图中横坐标表示归一化波数，算子长度2m＝16，如图3(a)所示，当r＝0.2，传统有限差分法仍然存在明显的频散误差；当r＝0.4(实际上是时间步长δt增大一倍)传统有限差分法的频散显著增加；分别对比图3(a)和图3(c)、图3(b)和图3(d)可以看出本发明方法可以在更大的波数范围内保持频散误差小，特别是在高波数情况下，表明本发明方法的模拟精度优于传统有限差分法。

由频散关系得到本实施例的有限差分法满足：

取那奎思波特波数(kx,kz)＝(π/h,π/h)代入频散误差表达式得到稳定条件为：

采用海绵吸收边界条件对边界反射进行吸收，利用求取的差分系数进行声波波场外推，得到任意时刻的波场及整个地震记录。记录波场快照切片，输出地震记录。

图4为本发明方法与传统有限差分法稳定性对比图，图中纵坐标为稳定性因子，该值越大表明方法更稳定，横坐标为差分算子长度2m，对比发现，当差分算子长度相同时，本方明方法的稳定性条件均大于传统有限差分法，这意味着本发明方法稳定性更好，可以采用更大的时间步长进行波场延拓。

图5为均匀速度模型中0.6s时刻的波场快照图。介质的地震波传播速度为1500m/s，模型网格区域为201×201点，网格的纵横空间点距为15m，震源位置在模型正中央，震源主频14.5hz，差分算子的长度2m＝16。对比图5(a)和图5(b)可知，当r由0.2增至0.4时(实际上是时间步长δt由2ms增大至4ms)，传统有限差分法的数值频散显著增强；再对比图5中(a)和(c)、(b)和(d)可知，本发明方法的模拟精度优于传统有限差分法，这与图3中的频散误差曲线分析一致，证实了本发明方法模拟精度更高。

图6为复杂的marmousi速度模型，模型网格区域为2721×701点，网格的纵横空间点距为20m，震源主频为10hz，位于模型地表的正中央，差分算子长度2m＝16。

图7为marmousi速度模型中的4.0s时刻的波场快照图，(a)为采用传统有限差分法，(b)为采用本发明方法图，从图7中黑色矩形框的放大图可以明显的看出，本发明方法提高了波场模拟精度，说明了本文发展的方法的有效性与稳健性。

在本实施例中，基于规则网格的差分结构除了包含轴向2m+1个离散点外，还包含4个非轴向离散点。与传统的差分结构相比，引入了额外的计算量。然而，新差分结构改善了传统时间2阶方法的精度，且通过稳定性分析可知，基于新差分结构的时空域有限差分法可以采用更大的时间步长。从推导出来的差分系数表达式可以看出，差分系数随着介质速度的变化而改变，本实施例的实现策略是预先计算并存储给定速度模型范围内差分系数。例如，对于给定的模型速度范围为1000m/s至5000m/s，本实施例从1000m/s开始以1m/s的增量计算并存储这些差分系数，时间4阶空间16阶差分格式只需要大约0.1mb存储这些差分系数。因此，本实施例的有限差分法和传统时间二阶方案计算内存几乎相等，本方法能够为计算量巨大的逆时偏移成像和全波形反演提供更为准确的波场信息，具有重要的实际用途。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者装置不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者装置所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者装置中还存在另外的相同要素。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质(如rom/ram、磁碟、光盘)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，空调器，或者网络设备等)执行本发明各个实施例所述的方法。

以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。