频率域粘弹性波正演模拟方法、装置、设备及存储介质

技术领域

本说明书涉及粘弹性波正演模拟技术领域，尤其是涉及一种频率域粘弹性波正演模拟方法、装置、设备及存储介质。

背景技术

针对垂向变化粘弹性介质的波传播模拟在勘探地震领域有着十分广泛的应用。比如，研究波在层状介质中的反射和透射规律、在积分方程法中获取背景场等均涉及到一维模型的全波场计算。在早期，半解析法被广泛应用于研究地震波在一维层状介质中的传播特征。根据Muller(1985)，用于研究水平层状介质波传播模拟和合成人工地震记录的半解析法主要可以分为射线理论、全波形理论、反射率法以及波数累加法。在这些方法中，反射率法应用最为广泛。该方法由Thomseon(1950)最早提出，之后诸多学者对该方法进行了改良和完善(Kennett,1980；Schimide和Tango,1986；Mallick和Frazer)。反射率法的最大优势在于它能给出三维点源在均匀层状模型中激发的全波场解析表达式。然而，尽管反射律法发展十分完善，但由于该方法涉及复杂的汉克尔变换，因此导致不同算法的计算精度和计算效率差距显著。此外，半解析法仅适用于层状均匀介质，对于具有垂向变化特征的粘弹性各向异性介质，通常难以推导其解析表达式。在这种情况下，数值方法显得尤为重要。

目前，大多数的数值模拟方法均针对二维和三维情况。在处理具有复杂分布的实际介质或地形条件时，这些二维或三维算法具有显著优势。然而，采用二维和三维正演模拟方法对垂向变化粘弹性介质进行波传播模拟，则计算代价过高。为此，JafarGandomi和Takenaka(2007)提出了一种模拟一维衰减介质平面波场的有限差分方法，该方法利用广义Zener模型描述介质衰减，并通过在τ-p域求解波动方程获得各个波场分量。之后，JafarGandomi和Takenaka(2013)将该方法进行了拓展，使其能够处理速度变化剧烈的一维模型。然而，这种方法仅能获得平面波解，而无法得到完整的波场解；此外，该方法只局限于各向同性介质。针对垂向变化各向异性粘弹性介质的数值方法对研究地震波在复杂一维介质中的传播规律具有重要价值，而目前仍被普遍忽略。因此，如何提高垂向变化粘弹性介质的波传播模拟的精度，已成为目前亟待解决的技术问题。

发明内容

本说明书实施例的目的在于提供一种频率域粘弹性波正演模拟方法、装置及设备，以提高频率域粘弹性波正演模拟的精度。

为达到上述目的，一方面，本说明书实施例提供了一种频率域粘弹性波正演模拟方法，包括：

获取频率波数域一维粘弹性波动方程；

将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；

将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；

加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；

对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；

根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；

对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。

本说明书一实施例中，所述频率波数域一维粘弹性波动方程，包括：

其中，

本说明书一实施例中，所述频率波数域一维粘弹性波动方程的等效积分弱形式为：

其中，Γ为计算区域；

本说明书一实施例中，所述有限元控制方程包括：

其中，K为全局刚度矩阵；M为全局质量矩阵；F为全局震源矢量；K,M,F的单元矩阵分别为

本说明书一实施例中，所述加载有所述三维点力源和所述吸收边界条件的有限元控制方程，包括：

其中，F的单元矩阵为

本说明书一实施例中，所述自适应波数采样策略包括：

将波数区间

最大波数设置为：k

在波数区间

三角单元的最大尺寸设置为：(Δk)

其中，Δk为三角单元的平均尺寸；N

本说明书一实施例中，所述频率域三维弹性波波场包括：

其中，u(ξ,ζ)为(ξ,ζ)下的频率域三维弹性波波场；n为二维波数域中第n个离散节点的编号，n＝1,2,…,N

另一方面，本说明书实施例还提供了一种频率域粘弹性波正演模拟装置，包括：

获取模块，用于获取频率波数域一维粘弹性波动方程；

转换模块，用于将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；

离散模块，用于将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；

加载模块，用于加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；

求解模块，用于对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；

确定模块，用于根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；

变换模块，用于对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。

另一方面，本说明书实施例还提供了一种计算机设备，包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时，执行上述方法的指令。

另一方面，本说明书实施例还提供了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被计算机设备的处理器运行时，执行上述方法的指令。

由以上本说明书实施例提供的技术方案可见，在本说明书实施例中，可以根据自适应波数采样策略确定频率波数域波场对应的波数集合，并对波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场；如此，不仅考虑了频率波数域粘弹性波场的振荡效应，并且辅之以自适应波数采样策略，从而实现三维点力源在任意一维各向异性粘弹性介质中的波场模拟；因此，与现有技术相比，本说明书实施例能够大大提升震源附近粘弹性波场的模拟精度。

附图说明

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

图1示出了本说明书一些实施例中频率域粘弹性波正演模拟方法的流程图；

图2示出了本说明书一实施例中刚性弱化(Strength Reduction Method，SRM) 吸收边界示意图；

图3a～图3d示出了本说明书一实施例中二维波数域波场的变化特征图；

图4示出了说明书一实施例中二维波数空间线单元剖分示意图；

图5示出了本说明书一实施例中各向异性全空间模型示意图；

图6示出了本说明书一实施例中各向异性全空间模型的解析解与数值解对比图；

图7示出了本说明书一实施例中层状各向异性粘弹性模型示意图；

图8a～图8d分别示出了本说明书一实施例中层状各向异性粘弹性模型的频率域位移场分量示意图；

图9示出了本说明书一些实施例中频率域粘弹性波正演模拟装置的结构框图；

图10示出了本说明书一些实施例计算机设备的结构框图。

【附图标记说明】

91、获取模块；

92、转换模块；

93、离散模块；

94、加载模块；

95、求解模块；

96、确定模块；

97、变换模块；

1102、计算机设备；

1104、处理器；

1106、存储器；

1108、驱动机构；

1110、输入/输出模块；

1112、输入设备；

1114、输出设备；

1116、呈现设备；

1118、图形用户接口；

1120、网络接口；

1122、通信链路；

1124、通信总线。

具体实施方式

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

参考图1所示，本说明书实施例的频率域粘弹性波正演模拟方法可以包括以下步骤：

S101、获取频率波数域一维粘弹性波动方程。

S102、将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式。

S103、将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程。

S104、加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程。

S105、对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场。

S106、根据自适应波数采样策略确定所述频率波数域波场对应的波数集合。

S107、对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。

在本说明书实施例中，可以根据自适应波数采样策略确定频率波数域波场对应的波数集合，并对波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场；如此，不仅考虑了频率波数域粘弹性波场的振荡效应，并且辅之以自适应波数采样策略，从而实现三维点力源在任意一维各向异性粘弹性介质中的波场模拟；因此，与现有技术相比，本说明书实施例能够大大提升震源附近粘弹性波场的模拟精度。

在本说明一些实施例中，在频率域中，根据弹性波运动方程、本构方程以及应变 -位移关系可以表示为：

其中：σ＝(σ

e＝(e

L为偏微分算子矩阵；σ为应力向量；f为震源矢量；u为频率域位移矢量；e为应变向量；C为粘弹性介质的复刚度矩阵；x，y，z分别为空间三维方向；ρ为介质密度；ω为角频率；ρ为介质密度；ω为角频率，u为频率域位移矢量；e为应变向量； T为矩阵转置；在应力向量σ中，每个元素的第一个下标表示曲面的法向量方向，每个元素的第二个下标表示应力方向；例如，以σ

对于沿深度方向的一维介质，C和ρ仅与z有关。因此，将频率波数域一维粘弹性波动方程关于x方向和y方向进行一维空间傅里叶变换，可以得到如下所示的频率波数域一维粘弹性波动方程：

其中，

在本说明书实施例中，将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式可以方便后续进行离散化处理。在本说明一些实施例中，所述将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式，可以包括：

1)、确定

2)、根据Galerkin加权余量法，取

3)、将

其中，Γ为计算区域；

在本说明一些实施例中，所述将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程，可以包括：

1)、将深度方向的计算区域Γ剖分为NE个线单元，(这里以线单元分析为例)，则

在每个线单元内，

其中，

2)、将上述

其中，K为全局刚度矩阵；M为全局质量矩阵；F为全局震源矢量；u为频率域位移矢量；K,M,F的单元矩阵分别为

在本说明书一些实施例中，三维点力源的加载可以选择雷克子波作为震源函数。频率域雷克子波函数可以表示为

f＝AR(f)δ(x)δ(y)δ(z-z

其中，

将f＝AR(f)δ(x)δ(y)δ(z-z

将

即在加载三维点力源后，全局震源矢量F的单元矩阵变为：

在本说明书一些实施例中，为了消除人工边界引起的虚假反射，可以在计算区域边界处施加SRM吸收边界条件(如图2所示)。在此吸收边界条件中，阻尼矩阵D可以由质量矩阵M和刚度矩阵K的线性组合表示为D＝D

式中，D

此外，为了增强SRM边界的吸收效果，还需要减小吸收层介质的刚度，即粘弹性介质的复刚度矩阵需满足：

其中，C(z′)为SRM区域的计算点到内边界的距离为x′时对应的粘弹性介质的复刚度矩阵，α(z′)＝α

在加入阻尼后，有限元控制方程可以表示为：

如此，加载有所述三维点力源和所述吸收边界条件的有限元控制方程可以表示为：

其中，全局震源矢量F的单元矩阵变为：

在本说明书一些实施例中，可以采用直接解法等求解

在获得频率波数域波场

根据频率波数域波场

基于此，在本说明书一些实施例中，利用三角单元对二维波数空间进行离散化(如图4所示)可以采用如下波数采样策略确定区间

式中，Δk为三角单元的平均尺寸；N

另外，为提高震源附近的计算精度，可以设置最大波数为k

在本说明书一些实施例中，基于形函数构建非等间隔二维傅里叶变换算法，可以获得频率域三维弹性波波场。为了方便处理，

其中，C＝1/(2π)且ξ＝x,ζ＝y,η＝k

在常规技术中，通常采用二维快速傅里叶变换计算。尽管这种做法计算效率较高，但由于二维快速傅里叶变换技术仅适用于矩形网格，因此无法准确描述频率-波数域 位移场的变化特征，从而导致数值精度较低。为了改善二维傅里叶变换的精度，这里 采用三角网格对波场进行精确刻画，具体步骤如下：

将二维波数区域

或

式中，I

为了更好地描述频率波数域位移场的振荡特性，这里仅对核函数

式中，m＝1,2,…,6分别为第e个三角单元中第m个离散节点编号；

其中，

将

其中，

考虑到三角单元形状的任意性，直接由

首先，根据关系式

其中：

χ

χ

χ

在此基础上，将

其中，u(ξ,ζ)为(ξ,ζ)下的频率域三维弹性波波场；n为二维波数域中第n个离散节点的编号，n＝1,2,…,N

获得的频率域三维弹性波波场，可以用于研究地震波在目标地球介质中的传播规律，从而为后续对目标地球介质的油气资源进行勘探开发提供参考依据。

图5示出了本说明书一示例性实施例中的各向异性全空间介质模型。基于此模型，利用上述的频率域粘弹性波正演模拟方法，可以获得频率波数域波场；然后在自适应波数采样策略下，对不同波数下的频率波数域波场进行非等间隔二维空间反傅里叶变换，即可以获得对应的频率域三维弹性波场。图6中示出了各向异性全空间介质模型的数值解(图6中的虚线)和解析解(图6中的实线)。在图6的(a)中，Re(u

虽然上文描述的过程流程包括以特定顺序出现的多个操作，但是，应当清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行(例如使用并行处理器或多线程环境)。

参考图9所示，与上述的频率域粘弹性波正演模拟方法对应，本说明书一些实施例中还提供了频率域粘弹性波正演模拟装置，其可以包括：

获取模块91，可以用于获取频率波数域一维粘弹性波动方程；

转换模块92，可以用于将所述频率波数域一维粘弹性波动方程转换为等效积分弱形式；

离散模块93，可以用于将频率波数域一维粘弹性波动方程的等效积分弱形式离散化为有限元控制方程；

加载模块94，可以用于加载目标工区的三维点力源和吸收边界条件至所述有限元控制方程；

求解模块95，可以用于对加载有所述三维点力源和所述吸收边界条件的有限元控制方程进行求解，获得所述目标工区的频率波数域波场；

确定模块96，可以用于根据自适应波数采样策略确定所述频率波数域波场对应的波数集合；

变换模块97，可以用于对所述波数集合中每个波数作为输入时对应的频率波数域波场，进行非等间隔二维傅里叶反变换，获得对应的频率域三维弹性波波场。

为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本说明书时可以把各单元的功能在同一个或多个软件和/或硬件中实现。

本说明书的实施例还提供一种计算机设备。如图10所示，在本说明书一些实施例中，所述计算机设备1102可以包括一个或多个处理器1104，诸如一个或多个中央处理单元(CPU)或图形处理器(GPU)，每个处理单元可以实现一个或多个硬件线程。计算机设备1102还可以包括任何存储器1106，其用于存储诸如代码、设置、数据等之类的任何种类的信息，一具体实施方式中，存储器1106上并可在处理器1104上运行的计算机程序，所述计算机程序被所述处理器1104运行时，可以执行根据上述方法的指令。非限制性的，比如，存储器1106可以包括以下任一项或多种组合：任何类型的RAM，任何类型的ROM，闪存设备，硬盘，光盘等。更一般地，任何存储器都可以使用任何技术来存储信息。进一步地，任何存储器可以提供信息的易失性或非易失性保留。进一步地，任何存储器可以为计算机设备1102的固定或可移除部件。在一种情况下，当处理器1104执行被存储在任何存储器或存储器的组合中的相关联的指令时，计算机设备1102可以执行相关联指令的任一操作。计算机设备1102还包括用于与任何存储器交互的一个或多个驱动机构1108，诸如硬盘驱动机构、光盘驱动机构等。

计算机设备1102还可以包括输入/输出模块1110(I/O)，其用于接收各种输入(经由输入设备1112)和用于提供各种输出(经由输出设备1114)。一个具体输出机构可以包括呈现设备1116和相关联的图形用户接口1118(GUI)。在其他实施例中，还可以不包括输入/输出模块1110(I/O)、输入设备1112以及输出设备1114，仅作为网络中的一台计算机设备。计算机设备1102还可以包括一个或多个网络接口1120，其用于经由一个或多个通信链路1122与其他设备交换数据。一个或多个通信总线1124将上文所描述的部件耦合在一起。

通信链路1122可以以任何方式实现，例如，通过局域网、广域网(例如，因特网)、点对点连接等、或其任何组合。通信链路1122可以包括由任何协议或协议组合支配的硬连线链路、无线链路、路由器、网关功能、名称服务器等的任何组合。

本申请是参照本说明书一些实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理器的处理器以产生一个机器，使得通过计算机或其他可编程数据处理器的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理器以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理器上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

在一个典型的配置中，计算机设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。

内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。

计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算机设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体 (transitory media)，如调制的数据信号和载波。

本领域技术人员应明白，本说明书的实施例可提供为方法、系统或计算机程序产品。因此，本说明书实施例可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本说明书实施例可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书实施例，在这些分布式计算环境中，由通过通信网络而被连接的远程处理器来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本说明书实施例的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

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