CN108828659A - 地震波场延拓方法及装置 - Google Patents

Abstract

本申请实施例提供了一种地震波场延拓方法及装置，该方法包括：确定用于弹性波模拟的一阶递归时间弹性波场延拓算子；利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束；利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果；基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。本申请实施例可以实现在满足大时间步长波场延拓且能保证弹性波场延拓计算稳定的前提下，提高弹性波场延拓的精度。

Description

地震波场延拓方法及装置

技术领域

本申请涉及地震勘探技术领域，尤其是涉及一种地震波场延拓方法及装置。

背景技术

地震勘探是寻找油气资源，解决油气勘探、开发问题的有效方法。地震正演、地震成像和地震反演是地震勘探中的重要技术，而这三种技术的计算效率和精度都依赖于所采用的时间域地震波场延拓方法。目前最常用的两类波场延拓方法包括有限差分方法和谱方法。有限差分利用差分代替微分，具有计算量小且容易实现的优点，被广泛应用于地震勘探技术中。有限差分的系数可以通过泰勒(Taylor)级数展开或是优化算法来求取。前者可以看作是对伪谱算子级数展开的截断，而后者可以看作是对伪谱算子在某个或某些特定频率上的最小二乘拟合。有限差分本质上是对伪谱算子的一种近似，这种近似使得有限差分存在计算精度低，频散严重以及计算不稳定等问题。在实际计算过程中，为了提高有限差分的计算精度，往往采用高阶的空间差分算子，但由于计算机存储等限制，在时间上通常还是采用二阶精度的差分，因此波场延拓在时间方向上精度较低，这也导致了频散误差严重。另外，为了保证计算的稳定，需要采用比地震采集数据小得多的时间步长进行波场延拓，增加了计算量。

近年来随着计算机技术的发展，在波数-时空域进行波场延拓成为可能。在波数-时空域构建的波场延拓算子能够补偿时间离散引起的误差，即便采用大时间步长延拓，仍能保持极高的精度和稳定性，对于依赖于波场延拓方法的地震正演、成像和反演技术的发展具有重要意义。但是由于实际地震数据量巨大，直接利用波数-时空域算子处理实际数据依然受到计算速度和计算存储的限制。将空间波数域算子和有限差分方法相结合，发挥空间波数域的波场延拓算子精度高、稳定性好的优点以及有限差分计算速度快的优点，是解决波场延拓问题的一种新思路。交错网格Lowrank有限差分方法就是基于这种思想提出的一种新的波场延拓方法。该方法利用傅里叶有限差分方法处理空间-波数域的算子，有效降低波数-时空域算子矩阵的秩，实现波数-时空域算子的Lowrank分解，构建傅里叶有限差分计算格式，从而在保证精度的前提下，大大节省了计算量。该方法虽然精度较高，但时间步长的选择依然受到稳定性条件的限制。因此，提供一种既能满足大时间步长延拓且能保证计算稳定的高精度波场延拓方法是目前本领域所普遍期望的，这对于实际地震资料的处理具有重要价值。

发明内容

本申请实施例的目的在于提供一种地震波场延拓方法及装置，以实现在满足大时间步长波场延拓且能保证弹性波场延拓计算稳定的前提下，提高弹性波场延拓的精度。

为达到上述目的，一方面，本申请实施例提供了一种地震波场延拓方法，包括：

确定用于弹性波模拟的一阶递归时间弹性波场延拓算子；

利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束；

利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果；

基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。

本申请实施例的地震波场延拓方法，所述用于弹性波模拟的一阶递归时间弹性波场延拓算子，包括：

其中，c_p,s为纵波速度或横波速度，k为波数矢量，k_i为i方向上的波数分量，且满足|k|为取k的绝对值，Δt为时间采样间隔，Δi为笛卡尔坐标系中i方向上的空间采样间隔，x为笛卡尔坐标系空间坐标。

本申请实施例的地震波场延拓方法，所述衰减约束后的一阶递归时间弹性波场延拓算子，包括：

其中，tp(k)为与波数相关的稳定约束算子，且

其中，k为波数，k₀为特征波数，δ为控制衰减的宽度，ε为控制衰减的大小。

本申请实施例的地震波场延拓方法，所述利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，根据如下公式实现：

其中，tp(k)为与波数相关的稳定约束算子，且

k_x为笛卡尔坐标系中x方向上的波数分量，Δx为笛卡尔坐标系中x方向上的空间采样间隔，k为波数矢量，|k|为取k的绝对值，x为笛卡尔坐标系空间坐标，c_0,p/s为参考速度，Δt为时间采样间隔，c_p,s为纵波速度或横波速度，a,b_x,b_z为分解得到的有限差分系数，Δz为笛卡尔坐标系中z方向上的空间采样间隔。

本申请实施例的地震波场延拓方法，所述基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，包括：

令

对进行空间傅里叶反变换，得到q(x,z)；

则获得时空域弹性波场延拓结果为aq(x,z)+2b_x[q(x+1,z)+q(x-1,z)]+

2b_z[q(x,z+1)+q(x,z-1)]；

其中，k_x为笛卡尔坐标系中x方向上的波数分量，k_z为笛卡尔坐标系中z方向上的波数分量，为波数域弹性波场，q为时空域弹性波场，+1表示向前一个网格，-1表示向后一个网格，x为笛卡尔坐标系中x方向上的坐标，z为笛卡尔坐标系中z方向上的坐标，a,b_x,b_z为有限差分系数。

另一方面，本申请实施例还提供了一种地震波场延拓装置，包括：

波场延拓算子确定模块，用于确定用于弹性波模拟的一阶递归时间弹性波场延拓算子；

波场延拓算子约束模块，用于利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束；

波场延拓算子分解模块，用于利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果；

弹性波场延拓模块，用于基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。

本申请实施例的地震波场延拓装置，所述用于弹性波模拟的一阶递归时间弹性波场延拓算子，包括：

其中，c_p,s为纵波速度或横波速度，k为波数矢量，k_i为i方向上的波数分量，且满足|k|为取k的绝对值，Δt为时间采样间隔，Δi为笛卡尔坐标系中i方向上的空间采样间隔，x为笛卡尔坐标系空间坐标。

本申请实施例的地震波场延拓装置，所述衰减约束后的一阶递归时间弹性波场延拓算子，包括：

其中，tp(k)为与波数相关的稳定约束算子，且

其中，k为波数，k₀为特征波数，δ为控制衰减的宽度，ε为控制衰减的大小。

本申请实施例的地震波场延拓装置，所述利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，根据如下公式实现：

其中，tp(k)为与波数相关的稳定约束算子，且

k_x为笛卡尔坐标系中x方向上的波数分量，Δx为笛卡尔坐标系中x方向上的空间采样间隔，k为波数矢量，|k|为取k的绝对值，x为笛卡尔坐标系空间坐标，c_0，p/s为参考速度，Δt为时间采样间隔，c_p,s为纵波速度或横波速度，a,b_x,b_z为分解得到的有限差分系数，Δz为笛卡尔坐标系中z方向上的空间采样间隔。

本申请实施例的地震波场延拓装置，所述基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，包括：

令

对进行空间傅里叶反变换，得到q(x,z)；

则获得时空域弹性波场延拓结果为aq(x,z)+2b_x[q(x+1,z)+q(x-1,z)]+2b_z[q(x,z+1)+q(x,z-1)]；

其中，k_x为笛卡尔坐标系中x方向上的波数分量，k_z为笛卡尔坐标系中z方向上的波数分量，为波数域弹性波场，q为时空域弹性波场，+1表示向前一个网格，-1表示向后一个网格，x为笛卡尔坐标系中x方向上的坐标，z为笛卡尔坐标系中z方向上的坐标，a,b_x,b_z为有限差分系数。

由以上本申请实施例提供的技术方案可见，在本申请实施例中，通过稳定约束算子对用于弹性波模拟的一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束，可以提高弹性波场延拓的稳定性；在此基础上，利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，并基于傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，从而可以实现大时间步长的高精度的弹性波场延拓。

附图说明

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

图1为本申请一实施例中地震波场延拓的流程图；

图2a为本申请一实施例中部分Marmousi2模型的纵波速度场；

图2b为本申请一实施例中部分Marmousi2模型的横波速度场；

图3a为对部分Marmousi2模型采用常规有限差分方法(10m空间采样间隔、时间步长1.2ms)得到的X分量模拟结果；

图3b为对部分Marmousi2模型采用常规有限差分方法(10m空间采样间隔、时间步长1.2ms)得到的Z分量模拟结果；

图4a为对部分Marmousi2模型采用本申请实施例的地震波场延拓方法(10m空间采样间隔、时间步长2.5ms)得到的X分量模拟结果；

图4b为对部分Marmousi2模型采用本申请实施例的地震波场延拓方法(10m空间采样间隔、时间步长2.5ms)得到的Z分量模拟结果；

图5为本申请一实施例中地震波场延拓装置的结构框图；

图6为本申请另一实施例中地震波场延拓装置的结构框图。

具体实施方式

为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。例如在下面描述中，在第一部件上方形成第二部件，可以包括第一部件和第二部件以直接接触方式形成的实施例，还可以包括第一部件和第二部件以非直接接触方式(即第一部件和第二部件之间还可以包括额外的部件)形成的实施例等。

而且，为了便于描述，本申请一些实施例可以使用诸如“在…上方”、“在…之下”、“顶部”、“下方”等空间相对术语，以描述如实施例各附图所示的一个元件或部件与另一个(或另一些)元件或部件之间的关系。应当理解的是，除了附图中描述的方位之外，空间相对术语还旨在包括装置在使用或操作中的不同方位。例如若附图中的装置被翻转，则被描述为“在”其他元件或部件“下方”或“之下”的元件或部件，随后将被定位为“在”其他元件或部件“上方”或“之上”。

参考图1所示，本申请实施例的地震波场延拓方法可以包括以下步骤：

S101、确定用于弹性波模拟的一阶递归时间弹性波场延拓算子。

在本申请一些实施例中，地震波场在地下介质中的传播，可以利用表征弹性波在无限大的、各向同性介质中的二阶弹性位移方程出发，如下所示：

其中，u＝(u_x,u_z)是位移矢量，t为时间变量，c_p和c_s分别为纵、横波速度值，x＝(x,z)为笛卡尔坐标系空间坐标， 和分别代表时空域梯度、散度和旋度运算符。

利用空间傅里叶变换，可得到波数域弹性波动方程：

其中，代表波数域中的位移矢量，k＝(k_x,k_z)是波数矢量，k，k·，k×分别表示波数域梯度、散度和旋度运算符。考虑到波场的矢量特性及线性特性，并应用波场分解理论，即总的分量波场可以描述为对应的P波分量和S波分量之和的形式同时，根据P波为无散场以及S波为无旋场，P波分量可以进一步表示为对波场求散度再求梯度的结果，S波可以进一步表示为对波场求两次旋度的结果，即：

将上述波数域P波和S波的表示形式代入到波数域中的弹性波动方程(2)中，可以得到P波和S波分离传播方程：

分别求解上述微分方程得到纵波及横波的解析解如下所示：

将二阶方程进行时间方向蛙跳格式离散，利用频率ω_p,s与波数的关系ω_p,s＝|k|c_p,s以及傅里叶变换并结合式(5)和式(6)，可以得到准确的二阶时间差分递归形式：

其中，表示傅里叶变换，sin表示正弦函数，Δt为时间采样间隔。

从式(9)及式(10)中可以看出，用于弹性波延拓的二阶递归时间弹性波场延拓算子可定义为：

其中，c_p,s代表纵波速度或横波速度，i＝x,z，k_i表示i方向上的波数分量，并满足|k|表示为取波数矢量k的绝对值。

其次，基于波数域卷积变换性质，并结合一阶算子和二阶算子的关系以及交错网格的定义，可以从二阶算子式(11)中导出一阶递归时间弹性延拓算子：

其中，Δi为笛卡尔坐标系中i方向上的空间采样间隔，x为笛卡尔坐标系空间坐标，指数项使得波场在计算时沿着计算方向的正向或负向交错Δx/2个网格步长，或者沿z轴的正向或负向交错Δz/2个网格步长。

S102、利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束。

在本申请一些实施例中，为了提高大时间步长波场延拓的稳定性，可以给出对传播算子进行衰减约束的方法。交错网格一阶递归时间延拓算子式(12)是一个光滑震荡的函数，根据Von-Neumann稳定性分析可知，计算稳定的充分条件是交错网格一阶递归时间弹性波场延拓算子的取值对于任意情况都位于区间[-1,+1]内。但在实际计算过程中，数值误差使得交错网格一阶递归时间弹性波场延拓算子在时间步长选择较大时，往往不能满足稳定性条件。为此，需对交错网格一阶递归时间弹性波场延拓算子施加衰减约束，施加衰减约束之后的交错网格一阶递归时间弹性波场延拓算子形式可以表示为：

其中，tp(k)是与波数相关的稳定约束算子。所述稳定约束算子为一衰减函数。在本申请一些实施例中，这里衰减函数的选取可满足如下三个条件：

1)、衰减函数是关于波数k的光滑函数；

2)、能够自动识别出一阶递归时间弹性波场延拓算子的稳定临界条件；

3)、能在临界值附近加入可控的衰减。

在本申请一示例性实施例中，稳定约束算子例如可以选用取如下形式：

上式中，tp(k)关于波数k的任意阶导函数存在(即该函数是光滑的)，通常k的取值范围是从负Nyquist波数到正Nyquist波数。k₀，ε和δ是三个控制参数，k₀为特征波数，是使得一阶递归时间延拓算子达到稳定条件所限制的临界值的波数，δ为控制衰减的宽度，ε为控制衰减的大小。

S103、利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果。

虽然可采用式(13)进行波场延拓，但需要在每一个时刻上做多次空间Fourier变换，从而导致计算量巨大。为此，在本申请一些实施例中，可利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，从而可以极大降低初始算子的秩，进而可以大大降低求解计算量。下面给出傅里叶有限差分Lowrank分解的主要原理。

以衰减约束后的一阶递归时间弹性波场延拓算子中x方向的偏导数为例，在空间-波数域求解格式具有如下形式，

上式中的由下列式(16)进行恒等变形

其中，c_0,p/s为参考速度，c_0,p/s的值与模型速度场有关，k_x为笛卡尔坐标系中x方向上的波数分量，c_p,s为纵波速度或横波速度。式(16)中可通过傅里叶变换在波数域进行计算，具体可以是先利用傅里叶变换将时空域波场变换到波数域进行相关计算，计算完后再反变换到时空域中。

对于式(16)中的同模型速度场相关的另一项在k＝0处进行泰勒展开，可得下式：

其中，a,b_x,b_z为分解得到的有限差分系数，Δz为笛卡尔坐标系中z方向上的空间采样间隔，k_z为笛卡尔坐标系中z方向上的波数分量。

对于式(17)，可利用最小二乘拟合计算可得傅里叶有限差分法差分系数a,b_x,b_z，该差分系数随速度场变化。通过式(17)的分解，式(16)中的可在时空域进行计算。具体可以是利用上述波数域计算得到的结果，经过傅里叶反变换得到时空域的量，并通过式(17)计算得到的有限差分系数进行差分计算。

通过如上所示的算子分解，可以将原来秩很大的一阶递归时间弹性波场延拓算子变为秩为1的伪谱算子同有限差分算子结合的形式，进而实现Lowrank分解。相比直接在空波数域求解一阶递归时间弹性波场延拓算子，这种方法可大大降低计算量。

S104、基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。

在本申请一些实施例中，将步骤S103得到的傅里叶有限差分格式用于对目标震源进行弹性波场延拓，可实现大时间步长稳定的地震波场延拓。在本申请一示例性实施例中，下面以波场延拓过程中时空域弹性波场q(x,z)沿着x轴正方向空间偏导数计算为例，介绍傅立叶有限差分计算流程：将时空域弹性波场q(x,z)进行空间傅里叶变换，可以得到波数域波场具体流程可如下所示：

令

对进行空间傅里叶反变换，得到q(x,z)；

则获得时空域弹性波场延拓结果为aq(x,z)+2b_x[q(x+1,z)+q(x-1,z)]+2b_z[q(x,z+1)+q(x,z-1)]。

其中，为波数域弹性波场，q为时空域弹性波场，+1表示向前一个网格，-1表示向后一个网格，x为笛卡尔坐标系中x方向上的坐标，z为笛卡尔坐标系中z方向上的坐标。

地震勘探通常采用较大的时间采用间隔记录地震信号，而常规有限差分方法采用大时间步长进行数值波场延拓时计算结果不稳定。例如图3a所示，为对图2a和图2b所示部分Marmousi2模型采用常规有限差分方法(10m空间采样间隔、时间步长1.2ms)得到的X分量模拟结果。再例如图3b所示，为对图2a和图2b所示部分Marmousi2模型采用常规有限差分方法(10m空间采样间隔、时间步长1.2ms)得到的Z分量模拟结果。

而在本申请的一些实施例中，如图4a所示，为对图2a和图2b所示部分Marmousi2模型采用本申请实施例的地震波场延拓方法(10m空间采样间隔、时间步长1.2ms)得到的X分量模拟结果。再比如图4b所示，为对图2a和图2b所示部分Marmousi2模型采用本申请实施例的地震波场延拓方法(10m空间采样间隔、时间步长1.2ms)得到的Z分量模拟结果。由此可以看出：当使用更大时间步长时，本申请实施例的地震波场延拓方法的计算仍然稳定，而且得到的波场精确，没有数值噪音。

参考图5所示，本申请实施例的一种地震波场延拓装置可以包括：

波场延拓算子确定模块51，可以用于确定用于弹性波模拟的一阶递归时间弹性波场延拓算子；

波场延拓算子约束模块52，可以用于利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束；

波场延拓算子分解模块53，可以用于利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果；

弹性波场延拓模块54，可以用于基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。

参考图6所示，本申请实施例的另一种地震波场延拓装置可以包括存储器、处理器、以及存储在所述存储器上的计算机程序，所述计算机程序被所述处理器运行时执行如下步骤：

确定用于弹性波模拟的一阶递归时间弹性波场延拓算子；

利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束；

利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果；

基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。

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

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

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

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

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

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

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

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

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

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

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

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

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

Claims (10)

1.一种地震波场延拓方法，其特征在于，包括：

确定用于弹性波模拟的一阶递归时间弹性波场延拓算子；

利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束；

利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果；

基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。

2.如权利要求1所述的地震波场延拓方法，其特征在于，所述用于弹性波模拟的一阶递归时间弹性波场延拓算子，包括：

其中，c_p,s为纵波速度或横波速度，k为波数矢量，k_i为i方向上的波数分量，且满足|k|为取k的绝对值，Δt为时间采样间隔，Δi为笛卡尔坐标系中i方向上的空间采样间隔，x为笛卡尔坐标系空间坐标。

3.如权利要求2所述的地震波场延拓方法，其特征在于，所述衰减约束后的一阶递归时间弹性波场延拓算子，包括：

其中，tp(k)为与波数相关的稳定约束算子，且

其中，k为波数，k₀为特征波数，δ为控制衰减的宽度，ε为控制衰减的大小。

4.如权利要求1所述的地震波场延拓方法，其特征在于，所述利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，根据如下公式实现：

其中，tp(k)为与波数相关的稳定约束算子，且

k_x为笛卡尔坐标系中x方向上的波数分量，Δx为笛卡尔坐标系中x方向上的空间采样间隔，k为波数矢量，|k|为取k的绝对值，x为笛卡尔坐标系空间坐标，c_0,p/s为参考速度，Δt为时间采样间隔，c_p,s为纵波速度或横波速度，a，b_x,b_z为分解得到的有限差分系数，Δz为笛卡尔坐标系中z方向上的空间采样间隔。

5.如权利要求1所述的地震波场延拓方法，其特征在于，所述基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，包括：

令

对进行空间傅里叶反变换，得到q(x，z)；

则获得时空域弹性波场延拓结果为aq(x，z)+2b_x[q(x+1,z)+q(x-1,z)]+2b_z[q(x,z+1)+q(x,z-1)]；

其中，k_x为笛卡尔坐标系中x方向上的波数分量，k_z为笛卡尔坐标系中z方向上的波数分量，为波数域弹性波场，q为时空域弹性波场，+1表示向前一个网格，-1表示向后一个网格，x为笛卡尔坐标系中x方向上的坐标，z为笛卡尔坐标系中z方向上的坐标，a,b_x,b_z为有限差分系数。

6.一种地震波场延拓装置，其特征在于，包括：

波场延拓算子确定模块，用于确定用于弹性波模拟的一阶递归时间弹性波场延拓算子；

波场延拓算子约束模块，用于利用预设的稳定约束算子对所述一阶递归时间弹性波场延拓算子的不稳定区域进行衰减约束；

波场延拓算子分解模块，用于利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，获得傅里叶有限差分Lowrank分解结果；

弹性波场延拓模块，用于基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，获得弹性波波场。

7.如权利要求6所述的地震波场延拓装置，其特征在于，所述用于弹性波模拟的一阶递归时间弹性波场延拓算子，包括：

其中，c_p,s为纵波速度或横波速度，k为波数矢量，k_i为i方向上的波数分量，且满足|k|为取k的绝对值，Δt为时间采样间隔，Δi为笛卡尔坐标系中i方向上的空间采样间隔，x为笛卡尔坐标系空间坐标。

8.如权利要求7所述的地震波场延拓装置，其特征在于，所述衰减约束后的一阶递归时间弹性波场延拓算子，包括：

其中，tp(k)为与波数相关的稳定约束算子，且

其中，k为波数，k₀为特征波数，δ为控制衰减的宽度，ε为控制衰减的大小。

9.如权利要求6所述的地震波场延拓装置，其特征在于，所述利用傅里叶有限差分方法对衰减约束后的一阶递归时间弹性波场延拓算子进行Lowrank分解，根据如下公式实现：

其中，tp(k)为与波数相关的稳定约束算子，且

k_x为笛卡尔坐标系中x方向上的波数分量，Δx为笛卡尔坐标系中x方向上的空间采样间隔，k为波数矢量，|k|为取k的绝对值，x为笛卡尔坐标系空间坐标，c_0,p/s为参考速度，Δt为时间采样间隔，c_p,s为纵波速度或横波速度，a,b_x,b_z为分解得到的有限差分系数，Δz为笛卡尔坐标系中z方向上的空间采样间隔。

10.如权利要求6所述的地震波场延拓装置，其特征在于，所述基于所述傅里叶有限差分Lowrank分解结果对目标震源进行弹性波场延拓，包括：

令

对进行空间傅里叶反变换，得到q(x,z)；

则获得时空域弹性波场延拓结果为aq(x,z)+2b_x[q(x+1,z)+q(x-1,z)]+2b_z[q(x,z+1)+q(x,z-1)]；

其中，k_x为笛卡尔坐标系中x方向上的波数分量，k_z为笛卡尔坐标系中z方向上的波数分量，为波数域弹性波场，q为时空域弹性波场，+1表示向前一个网格，-1表示向后一个网格，x为笛卡尔坐标系中x方向上的坐标，z为笛卡尔坐标系中z方向上的坐标，a,b_x,b_z为有限差分系数。