波场正演模拟方法及装置

摘要

本发明实施例提供了一种波场正演模拟方法及装置，其中，该方法包括：获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；根据最大范数原理和所述窗函数，创建目标函数；求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述最优解代入所述窗函数；利用所述窗函数得到优化后的均匀网格有限差分算子；利用优化后的均匀网格有限差分算子进行波场正演模拟。该方案优化了目标函数中待优化的系数，首次实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散。

说明书

技术领域

本发明涉及数据处理技术领域，特别涉及一种波场正演模拟方法及装置。

背景技术

基于有限差分格式的地震波场正演模拟在逆时偏移中得到了广泛的应用，全波反演、逆时偏移都需要进行多次正演和波场的反传，因此，正演模拟的精度和速度是非常重要的考量因素，其对成像和反演具有重要意义。许多学者对波动方程的有限差分格式进行了广泛的研究。为了减少数值频散，已经开发了多种有限差分格式，例如，可变网格、不规则网格、交错网格和旋转交错网格等。

然而，空间导数上的高阶有限差分算子的经典系数通常是由空间导数项的泰勒级数展开确定的。仅改变有限差分格式不能最小化数值误差。如果传统的有限差分系数用于地震波场正演模拟，强的数值频散是不可避免的，特别是对于波动方程中的高波数范围。在最近的二十年中，学者们已经提出了许多优化方法，如牛顿方法、隐格式、基于时空域频散关系的方法、最小二乘法等。但是这类算法使用的目标函数参数过多，需要对与有限差分算子阶数相同数量的参数进行优化，对优化算法的要求提出了更高的要求。与上述方法相比，基于窗函数的优化方法是非常灵活的，我们认为，截断伪谱方法的窗函数决定了有限差分算子的精度，如今各种窗函数被用来得到有限差分系数。有限差分法的另一缺点是计算成本、可变时间步长和自适应变长空间算子是从不同方面减少计算时间的两种方式。然而，这些方法不能从根本上解决计算速度的问题。

上世纪60年代末期，有人提出将二阶显格式的有限差分技术应用于层状介质的弹性波数值模拟，拉开了有限差分深入研究的序幕。进而发展了适应非均匀介质的有限差分格式以进行弹性波数值模拟，还发展了高阶数的有限差分格式进行声波方程求解。有人将时间域有限差分方法应用于粘声波方程，并进行了基于各向异性介质的逆时偏移(RTM)。高阶有限差分算子在空间导数上的传统系数通常由泰勒级数展开确定。采用常规的有限差分系数进行地震波场正演模拟，会出现较强的数值频散，特别是波动方程中的高波数范围。近二十年来，不同学者提出了牛顿法、隐式格式、模拟退火算法、最小二乘法。但是这些优化方法实施起来非常复杂，并且由于目标函数的影响而不能广泛用于偏移和反演。与上述方法相比，基于窗口函数的优化方法实现起来非常灵活，有限差分法是伪谱法空间卷积序列的空间截断形式。不同学者利用不同窗函数得到有限差分系数。有限差分法的另一个缺点是计算成本，可变时间步长和自适应可变长度空间算子，是从不同方面减少计算时间的两种方法。然而，这些方法不能从根本上解决计算速度问题，但是GPU技术会带来相当大的加速效应，有人提出使用单GPU技术和窗函数优化方法进行弹性波数值模拟。但是对于大型模型或者三维模型，由于内存限制，单GPU不再适用。

综上，目前主流的有限差分算子优化算法有两大类，第一类是构建目标函数，利用最优化算法进行求解，但其缺点是目标函数过于复杂(参数过多)，限制了优化的效果；第二类是窗函数算法，其缺点是不同窗函数的优化效果很难定量控制，使得在进行地震波场正演模拟时存在较强的数值频散的问题。

发明内容

本发明实施例提供了一种波场正演模拟方法，以解决现有技术中有限差分算子优化过程中目标函数复杂、地震波场正演模拟时存在较强的数值频散的技术问题。该方法包括：

获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；

根据最大范数原理和所述窗函数，创建目标函数；

求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；

将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；

利用所述窗函数得到优化后的均匀网格有限差分算子；

利用优化后的均匀网格有限差分算子进行波场正演模拟。

本发明实施例还提供了一种波场正演模拟装置，以解决现有技术中有限差分算子优化过程中目标函数复杂、地震波场正演模拟时存在较强的数值频散的技术问题。该装置包括：

窗函数获取模块，用于获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；

目标函数创建模块，用于根据最大范数原理和所述窗函数，创建目标函数；

求解模块，用于求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；

优化模块，用于利用所述窗函数得到优化后的均匀网格有限差分算子；

模拟模块，用于利用优化后的均匀网格有限差分算子进行波场正演模拟。

在本发明实施例中，提出了一种新的窗函数，该窗函数包括有限差分系数和整调系数，即该窗函数相对于现有技术中传统的窗函数添加了整调系数；进而基于最大范数原理和窗函数，创建目标函数，通过窗函数的加持，使得降低了目标函数中待优化的系数；通过求解目标函数，来得到有限差分系数的解和整调系数的解，将有限差分系数的解和整调系数的解，代入窗函数后，即可基于窗函数得到优化后的均匀网格有限差分算子，由于新的窗函数添加了整调系数，使得实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

图1是本发明实施例提供的一种波场正演模拟方法的流程图；

图2是本发明实施例提供的一种点源相应波场快照的示意图；

图3是本发明实施例提供的一种单道对比波形图；

图4是本发明实施例提供的一种波场正演模拟装置的结构框图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

在本发明实施例中，提供了一种波场正演模拟方法，如图1所示，该方法包括：

步骤101：获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；

步骤102：根据最大范数原理和所述窗函数，创建目标函数；

步骤103：求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；

步骤104：将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；

步骤105：利用所述窗函数得到优化后的均匀网格有限差分算子；

步骤106：利用优化后的均匀网格有限差分算子进行波场正演模拟。

由图1所示的流程可知，在本发明实施例中，提出了一种新的窗函数，该窗函数包括有限差分系数和整调系数，即该窗函数相对于现有技术中传统的窗函数添加了整调系数；进而基于最大范数原理和窗函数，创建目标函数，通过窗函数的加持，使得降低了目标函数中待优化的系数；通过求解目标函数，来得到有限差分系数的解和整调系数的解，将有限差分系数的解和整调系数的解，代入窗函数后，即可基于窗函数得到优化后的均匀网格有限差分算子，由于新的窗函数添加了整调系数，使得实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散，有利于获得更好的模拟精度。

本申请发明人发现，现有技术中通过以下方式获得均匀网格有限差分算子：

均匀网格有限差分算子可以通过sinc函数插值理论推导得到，sinc函数可以重构均匀采样后的带限信号fn(Diniz et al.,2012)，sinc函数可以表示为：

其中，_x表示采样点的位置；Δx为空间采样间隔，代表截止波数，f_n表示采样信号。利用窗函数截断公式(1)的一阶导数和二阶导数，获取了均匀网格有限差分算子，均匀网格有限差分算子求导后可表示为公式(2)和公式(3)，其中，f(n)表示原始信号。

利用窗函数截断公式(2)和(3)后，可以分别表示为公式(4)和公式(5)。

因n＝0为公式(2)和公式(3)的一个奇异点，根据信号采样理论，公式(2)和公式(3)可以表示为公式(6)和公式(7)。

其中，f₀表示中间位置的函数，f_n为中间位置的正方向的函数，f_-n为中间位置的负方向的函数，和为f_n的系数，这里ξ代表黎曼ξ函数，加窗截断后可以表示为公式(8)和公式(9)。

这里，w(n)表示窗函数，c_n表示正方向系数，c_-n表示反方向系数。

将公式(8)和公式(9)变换到频率域后分别为公式(10)和公式(11)

对于一阶导数，可以将误差函数表示成公式(12)。对于二阶导数，可以将误差函数表示成公式(13)。

可见，通过定义不同的窗函数可以得到均匀网格有限差分算子。还可以基于传统二项窗函数定义一个改进二项式窗函数，常规的二项式窗口函数可以表示为方程式(14)，该二项式窗函数改进过程中认为可以用N+M代替N，其中N是阶数，M是偶数。这种方法的最大缺点是不能减少在有限差分阶数较低的情况下的数值频散。

原始的二项式窗函数给出了传统的均匀网格有限差分算子，相当于从泰勒级数展开式得到的均匀网格有限差分算子，它保留了低波数部分的良好精度，但是对于高波数分量却没有改进效果。

本申请发明人提出了一种新的窗函数，该新的窗函数相对于现有技术中的窗函数添加了整调系数，具体的，该整调系数可以是多个参数也可以是1个参数，例如，本申请以整调系数为2个参数为例，该新的窗函数可以表示为等式(15)，从中我们可以发现，在新的窗函数中，我们增加了两个参数，分别为参数m和h，可以称之为整调系数。

其中，表示窗函数；_m和h表示整调系数；_n表示窗函数空间点位；N表示有限差分阶数。

基于新的窗函数可以得到包含新参数的优化后的均匀网格有限差分算子，例如，优化后的均匀网格有限差分算子如公式(16)所示：

其中，bnew_n表示优化后的均匀网格有限差分算子。

本申请发明人发现，由于现有技术中的窗函数中的未知数只包含了有限差分系数，使得现有技术中的窗函数无法与目标函数结合，本申请提出新的窗函数之后，该新的窗函数包含有限差分系数的同时还添加了整调系数，使得新的窗函数可以与目标函数进行结合，本实施例根据最大范数原理(例如，1范数或2范数)和上述新的窗函数创建目标函数，并通过求解目标函数，来得到有限差分系数的解和整调系数的解。该目标函数如公式(17)所示：

其中，k_x表示波数范围；Δx表示空间采样间隔；T表示优化时允许的最大误差。

具体实施时，为了获得更好的求解精度，在本实施例中，采用模拟退火算法求解上述目标函数，求得的有限差分算子如下表1所示，其中，整调系数m＝32.3，h＝15.6。

表1

10阶12阶16阶20阶0.85714285710.88888888890.90909090901.257430489308-0.2678571429-0.311111111-0.3409090909-0.1263778016890.11313131310.139860139860.037148969395-0.03535353534-0.05244755244-0.0139097550480.016783216780.005509987784-0.00437062937-0.0021305265260.000763273991-0.0002411120230.0000626766642-0.000011537920

解得有限差分系数的解和整调系数的解之后，将有限差分系数的解和整调系数的解代入上述公式(15)中的窗函数，再将窗函数代入公式(16)，即可得到优化后的均匀网格有限差分算子。在该优化均匀网格有限差分算子的过程，同时利用了最优化算法和窗函数优化算法对均匀网格有限差分算子进行优化，该算法在业内尚属首次，使得优化均匀网格有限差分算子可以同时获得窗函数算法和最优化算法的特点，在实用阶段可以获得更好的模拟精度。

以下结合示例说明上述波场正演模拟方法的优点。

例如，建立一个600×600的、网格间距为5m的二维均匀模型，但实际计算的模型尺寸为700×700，其中，每边包含50个吸收边界。纵波速度为2000m/s，横波速度为1400m/s，密度为1000。Ricker子波的主频率为50Hz，位于速度模型的中心。二维非均匀介质中的弹性波方程如公式(18)所示：

图2显示了点源相应波场的快照，抽取图2中虚线部位的数据，使用48阶有限差分算子作为参考解。在图3中，虚线为图2抽取的数据使用本申请优化后的均匀网格有限差分算子处理的数据，实线为参考解，图3中的(a)至(d)分别对应了图2的(a)至(d)中虚线波形的处理数据，图3中的(e)至(f)显示了图2的(a)至(d)中抽取的数据采用本申请优化后的均匀网格有限差分算子处理的数据与参考解的差值。显然，本申请中优化后的均匀网格有限差分算子处理的数据比传统方法和改进二项式窗有着更少的数值频散。当我们将这些方法与参考解比较时，可以注意到，优化后的均匀网格有限差分算子处理的数据更接近参考解，这也证明了上述误差分析的结论。

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

图4是本发明实施例的波场正演模拟装置的一种结构框图，如图4所示，该装置包括：

窗函数获取模块401，用于获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；

目标函数创建模块402，用于根据最大范数原理和所述窗函数，创建目标函数；

求解模块403，用于求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；

优化模块404，用于利用所述窗函数得到优化后的均匀网格有限差分算子；

模拟模块405，用于利用优化后的均匀网格有限差分算子进行波场正演模拟。

在一个实施例中，所述窗函数的表达式为：

其中，表示窗函数；m和h表示整调系数；n表示窗函数空间点位；N表示有限差分阶数。

在一个实施例中，优化后的均匀网格有限差分算子的表达式为：

其中，bnew_n表示优化后的均匀网格有限差分算子。

在一个实施例中，所述目标函数的表达式为：

其中，k_x表示波数范围；Δx表示空间采样间隔；T表示优化时允许的最大误差。

在一个实施例中，所述求解模块具体用于采用模拟退火算法求解所述目标函数。

在另外一个实施例中，还提供了一种软件，该软件用于执行上述实施例及优选实施方式中描述的技术方案。

在另外一个实施例中，还提供了一种存储介质，该存储介质中存储有上述软件，该存储介质包括但不限于：光盘、软盘、硬盘、可擦写存储器等。

本发明实施例实现了如下技术效果：提出了一种新的窗函数，该窗函数包括有限差分系数和整调系数，即该窗函数相对于现有技术中传统的窗函数添加了整调系数；进而基于最大范数原理和窗函数，创建目标函数，通过窗函数的加持，使得降低了目标函数中待优化的系数；通过求解目标函数，来得到有限差分系数的解和整调系数的解，将有限差分系数的解和整调系数的解，代入窗函数后，即可基于窗函数得到优化后的均匀网格有限差分算子，由于新的窗函数添加了整调系数，使得实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散。

显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。