CN116520418A - 一种弹性波角度域共成像点道集高效提取方法 - Google Patents

Abstract

本发明涉及勘探地球物理学领域，尤其涉及一种弹性波角度域共成像点道集高效提取方法。本发明利用交错网格高阶有限差分法求解声波方程，同时记录各网格点最大振幅时刻作为该点成像时刻，并保留该时刻震源波场；利用交错网格高阶有限差分法求解无转换解耦延拓方程，获得各网格点成像时刻的纵、横波检波波场；利用各网格点成像时刻震源波场和检波波场分别计算对应的能流密度矢量，进而根据余弦定理求取各网格点纵、横波反射角；并以反射角信息为索引，对各网格点成像时刻震源波场和检波波场应用激发振幅成像条件，获得纵、横波角度域共成像点道集。通过实施例，本发明能够在不影响精度的情况下，降低技术所需的计算资源和存储资源，提高计算效率。

Description

一种弹性波角度域共成像点道集高效提取方法

技术领域：

本发明属于勘探地球物理学领域，尤其涉及一种弹性波角度域共成像点道集高效提取方法。

背景技术：

角度域共成像点道集是偏移技术重要的输出结果，尤其是基于弹性逆时偏移技术提取的角度域共成像点道集，保留了逆时偏移技术在复杂构造区域内的成像优势，同时输出的纵、横波两种结果相较于纵波结果能够提供更精细的地下构造信息。逆时偏移技术本身存在高计算量、高存储量等计算效率问题；同时模拟弹性波波场相比模拟纵波波场需要占用更多的计算存储资源；此外，计算各个网格点的反射角也会占用较多的计算资源。因此，现在亟需一种弹性波角度域共成像点道集高效提取方法，在不影响角度域共成像点道集精度情况下，解决现有技术存在的高计算量、高存储量等计算效率问题。

发明内容：

为解决现有基于弹性逆时偏移角度域共成像点道集提取技术中高计算量、高存储量等计算效率问题，本发明提供了一种弹性波角度域共成像点道集高效提取方法，对技术流程中存在的高计算量、高存储量问题，在不影响成像精度的情况下，结合多个计算效率改善策略，实现提高基于弹性逆时偏移提取角度域共成像点道集技术的计算效率。

通过利用交错网格高阶有限差分法求解声波方程，同时记录各网格点最大振幅时刻作为该点成像时刻，并保留该时刻震源波场；利用交错网格高阶有限差分法求解无转换解耦延拓方程，获得各网格点成像时刻的纵、横波检波波场；利用各网格点成像时刻震源波场和检波波场分别计算对应的能流密度矢量，进而根据余弦定理求取各网格点纵、横波反射角；并以反射角信息为索引，对各网格点成像时刻震源波场和检波波场应用激发振幅成像条件，获得纵、横波角度域共成像点道集。

为了解决上述技术问题，本文的具体技术方案如下：

一方面，本文实施例提供了一种弹性波角度域共成像点道集高效提取方法，包括，

以介质纵波速度参数和震源子波作为输入，利用交错网格高阶有限差分法求解一阶声波声压速度方程，同时记录各网格点最大振幅时刻作为该点成像时刻，并保留该时刻震源波场；

以介质纵、横波速度参数和多分量地震记录作为输入，利用交错网格高阶有限差分法求解一阶速度-应力形式的无转换解耦延拓方程，获得各网格点成像时刻的纵、横波检波波场；

利用各网格点成像时刻震源波场和检波波场分别计算对应的能流密度矢量，进一步根据余弦定理求取各网格点纵、横波反射角；

对各网格点成像时刻震源波场和检波波场应用激发振幅成像条件，以各网格点反射角信息作为索引，进而提取纵、横波角度域共成像点道集。

进一步地，所述一阶声波声压速度方程为，

其中，是矢量纵波震源波场的质点振动速度分量，P是压力分量，ρ表示密度，V_P表示纵波速度，t表示时间，x和z分别表示二维笛卡尔坐标系x和z方向；

利用交错网格有限差分法进行离散，其离散格式如下：

其中，Δt表示时间采样间隔，(Δx,Δz)表示空间采样间隔，i、j分别表示所述二维笛卡尔坐标系中x、z方向的离散点，t表示时间，表示精度为2N阶的交错网格有限差分系数；

以波场初始值作为各网格点的最大振幅值，在震源波场正向时间延拓过程中，当某个网格点波场绝对值大于该网格点最大振幅值时，记录当前时刻作为该网格点的成像时刻，记录该网格点波场绝对值作为最大振幅，更新所记录网格点成像时刻质点振动速度和声压值。

进一步地，根据预先设定的介质纵、横速度模型，将地面记录的多分量地震资料作为边值条件，利用一阶应力-速度形式的无转换解耦延拓方程，实现地震检波波场的逆时延拓，构建检波矢量波场；其中所述一阶应力-速度形式无转换解耦延拓方程为：

其中，(v_x,v_z)是总波场的速度分量，是纵波波场的质点振动速度分量，是横波波场的质点振动速度分量，τ^P是纵波应力场，/>是横波应力场，x和z分别表示二维笛卡尔坐标系x和z方向，(V_P,V_S)分别代表纵、横波速度；

在检波波场逆时延拓过程中，提取空间上各网格点在成像时刻的质点振动速度和应力值。

进一步地，利用保存的震源波场各网格点成像时刻质点振动速度和检波波场各网格点成像时刻质点振动速度应力，计算各反射点纵、横波的能流密度矢量，具体表达式为：

其中，是各网格点的纵波震源波场能流密度矢量，通过成像时刻压力场P和质点振动速度场/>计算获得；/>是各网格点的纵波检波波场能流密度矢量，通过成像时刻纵波应力场τ^P和质点纵波振动速度场/>计算获得；是各网格点的横波检波波场能流密度矢量，通过成像时刻横波切应力场/>和质点纵波振动速度场/>计算获得；

利用各网格点震源波场和检波波场的能流密度矢量，基于余弦定理计算各网格点的反射角，具体表达式为：

其中，(θ^P,θ^S)分别是纵波反射角和横波反射角，表示纵波震源波场的能流密度矢量，/>分别表示纵波检波波场和横波检波波场的能流密度矢量。纵波反射角为纵波入射波和纵波反射波夹角的一半，而转换横波反射角为纵波入射波和横波反射波夹角减去纵波入射角(纵波入射角等于纵波反射角)。

进一步地，对震源波场和检波波场应用激发振幅成像条件获得纵波入射-纵波反射和纵波入射-横波反射两种成像结果；对两种成像结果分别以纵、横波反射角作为索引，构建纵、横波共成像点道集，具体表达式为：

其中，空间坐标为(x,z)的网格点其单炮纵、横波角度域共成像点道集为(I_PP,I_PS)，为纵、横波角度域共成像点道集上的角度值，选取成像时刻t的纵波震源波场值/>和纵、横检波波场值/>来计算成像值，/>是高斯采样函数，σ函数方差，(θ^P,θ^S)是纵、横波反射角。

利用本说明书实施例，通过基于弹性逆时偏移技术，首先利用交错网格高阶有限差分法求解声波方程，同时记录各网格点最大振幅时刻作为该点成像时刻，并保留该时刻震源波场；再利用交错网格高阶有限差分法求解无转换解耦延拓方程，获得各网格点成像时刻的纵、横波检波波场；然后利用各网格点成像时刻震源波场和检波波场分别计算对应的能流密度矢量，进而根据余弦定理求取各网格点纵、横波反射角；最后应用激发振幅成像条件并以反射角信息为索引，输出纵、横波角度域共成像点道集。通过实施例，本发明能够在不影响精度的情况下，降低技术所需的计算资源和存储储资源，提高计算效率。

附图说明：

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

图1所示为本文实施例一种弹性波角度域共成像点道集高效提取方法的流程示意图；

图2所示为本文实施例交错网格点定义的示意图；

图3所示为本文实施例震源波场正向延拓过程更新成像时刻及该时刻震源波场数据；

图4所示为本文实施例具有水平-倾斜分界面弹性模型的具体参数示意图；

图5-1和图5-2所示为本文实施例对图4所示模型所提取的纵波入射-纵波反射模式和纵波入射-横波反射模式两个不同模式下角度域共成像点道集；；

图6所示为本文实施例盐丘弹性模型的具体参数示意图；

图7-1和图7-2所示为本文实施例对图6所示模型提取的纵波入射-纵波反射模式和纵波入射-横波反射模式两个不同模式下角度域共成像点道集；

具体实施方式：

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

为了解决现有技术中存在的问题，本文实施例提供了一种弹性波角度域共成像点道集高效提取方法，基于弹性逆时偏移技术下，在不降低精度的前提下，能够提高提取弹性波角度域共成像点道集计算效率。图1所示为本文实施例一种弹性波角度域共成像点道集高效提取方法的流程示意图。在本图中描述了基于弹性波逆时偏移高效提取角度域共成像点道集的过程，但基于常规或者无创造性的劳动可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。具体的如图1所示，可以由处理器执行，所述方法可以包括：

以介质纵波速度参数和震源子波作为输入，利用交错网格高阶有限差分法求解声波一阶声压速度方程，同时记录各网格点最大振幅时刻作为该点成像时刻，并保留该时刻震源波场；

以介质纵、横波速度参数和多分量地震记录作为输入，利用交错网格高阶有限差分法求解一阶速度-应力形式的无转换解耦延拓方程，获得各网格点成像时刻的纵、横波检波波场；

利用各网格点成像时刻震源波场和检波波场分别计算对应的能流密度矢量，再根据余弦定理求取各网格点纵、横波反射角；

对各网格点成像时刻震源波场和检波波场应用激发振幅成像条件，以各网格点反射角信息作为索引，进而提取纵、横波角度域共成像点道集。

利用交错网格高阶有限差分法求解声波方程，同时记录各网格点最大振幅时刻作为该点成像时刻，并保留该时刻震源波场；利用交错网格高阶有限差分法求解无转换解耦延拓方程，获得各网格点成像时刻的纵、横波检波波场；利用各网格点成像时刻震源波场和检波波场分别计算对应的能流密度矢量，进而根据余弦定理求取各网格点纵、横波反射角；并以反射角信息为索引，对各网格点成像时刻震源波场和检波波场应用激发振幅成像条件，获得纵、横波角度域共成像点道集。通过实施例，本发明能够在不影响精度的情况下，降低技术所需的计算资源和存储资源，提高计算效率。

本文实施例所述的弹性波角度域共成像点道集高效提取方法包括如下步骤：

步骤1：根据预先给定的介质纵波速度模型，加载震源子波，利用声波一阶声压速度方程，实现地震波场的正向时间延拓，构建纵波震源波场，所述的一阶声波一阶声压速度方程如下：

一阶声压速度方程(1)：

其中，是矢量纵波震源波场的质点振动速度分量，P是压力分量，ρ表示密度，V_P表示纵波速度，t表示时间，x和z分别表示二维笛卡尔坐标系x和z方向；

交错网格网格点分布如图2所示，各网格点定义如公式(2)所示：

基于上述网格定义，利用交错网格高阶有限差分法对声波方程离散构建纵波震源波场，其差分离散格式为公式(3)：

其中，Δt表示时间采样间隔，(Δx,Δz)表示空间采样间隔，i、j分别表示所述二维笛卡尔坐标系中x、z方向的离散点，t表示时间，表示精度为2N阶的交错网格有限差分系数；

以波场初始值作为各网格点的最大振幅值，在震源波场正向时间延拓过程中，当某个网格点波场绝对值大于该网格点最大振幅值时，记录当前时刻作为该网格点的成像时刻，记录该网格点波场绝对值作为最大振幅，更新所记录网格点成像时刻质点振动速度和声压值，其步骤如图3所示。

步骤2：根据预先设定的介质纵、横速度模型，将地面记录的多分量地震资料作为边值条件，利用一阶应力-速度形式的无转换解耦延拓方程，实现地震检波波场的逆时延拓，构建检波矢量波场；其中所述一阶应力-速度形式无转换解耦延拓方程为：

一阶无转换解耦延拓方程为公式(4)：

其中，(v_x,v_z)是总波场的速度分量，是纵波波场的质点振动速度分量，是横波波场的质点振动速度分量，τ^P是纵波应力场，/>是横波应力场，x和z分别表示二维笛卡尔坐标系x和z方向，(V_P,V_S)分别代表纵、横波速度；

在检波波场逆时延拓过程中，提取空间上各网格点在成像时刻的质点振动速度和应力值。

步骤3：利用保存的震源波场各网格点成像时刻质点振动速度和检波波场各网格点成像时刻质点振动速度应力，计算各反射点纵、横波的能流密度矢量，如公式(5)所示：

其中，是各网格点的纵波震源波场能流密度矢量，通过成像时刻压力场P和质点振动速度场/>计算获得；/>是各网格点的纵波检波波场能流密度矢量，通过成像时刻纵波应力场τ^P和质点纵波振动速度场/>计算获得；是各网格点的横波检波波场能流密度矢量，通过成像时刻横波切应力场/>和质点纵波振动速度场/>计算获得；

利用各网格点震源波场和检波波场的能流密度矢量，基于余弦定理计算各网格点的反射角，如公式(6)所示：

其中，(θ^P,θ^S)分别是纵波反射角和横波反射角，表示纵波震源波场的能流密度矢量，/>分别表示纵波检波波场和横波检波波场的能流密度矢量。纵波反射角为纵波入射波和纵波反射波夹角的一半，而转换横波反射角为纵波入射波和横波反射波夹角减去纵波入射角(纵波入射角等于纵波反射角)。

步骤4：对震源波场和检波波场应用激发振幅成像条件获得纵波入射-纵波反射和纵波入射-横波反射两种成像结果；对两种成像结果分别以纵、横波反射角作为索引，构建纵、横波共成像点道集，，具体表达式如公式(7)所示：

其中，空间坐标为(x,z)的网格点其单炮纵、横波角度域共成像点道集为(I_PP,I_PS)，为纵、横波角度域共成像点道集上的角度值，选取成像时刻t的纵波震源波场值/>和纵、横检波波场值/>来计算成像值，/>是高斯采样函数，_σ函数方差，(θ^P,θ^S)是纵、横波反射角。

示例性地，结合本文所述的具有水平-倾斜分界面的弹性模型的应用，该模型及其参数如图4所示，在图4中，横坐标“Distance”表示模型宽度、纵坐标“Depth”表示模型深度、单位均为“km”，以30Hz主频Ricker子波震源函数激发纯纵波为输入，基于弹性逆时偏移提取角度域共成像点道集。

图5为图4模型的弹性波角度域共成像点道集，从第50个CMP开始，每隔50个CMP抽取一个角度域共成像点道集，共选取10个CMP点。从图中可以看出，当速度模型准确时，显示的同相轴都处在正确深度，且同相轴平直。在图5-1纵波角度域共成像点道集的浅层大角度区域，同相轴角度区间大，其原因是震源波场和检波波场波矢量夹角过大形成的低波数噪音。同时在在图5-2横波角度域共成像点道集，虽然在深层的倾斜反射面上显示的同相轴能量相比纵波角度域共成像点道集弱，但是仍有连续性好。在两个角度域共成像点道集的右侧区域，由于照明不足，其显示同相轴角度区间小于左侧区域。通过观察分析，纵、横波角度域共成像点道集都符合预期设想，验证了本发明具有较为理想的成像精度。

为了验证本发明弹性波角度域共成像点道集高效提取方法在复杂模型的适应性，将本文的提取流程用于盐丘模型成像实验，该模型速度如图6所示。

下面以具体弹性波角度域共成像点道集为例说明应用本发明成像精度。

盐丘模型在其中深层处设置了十个绕射点，图7为这十个绕射点处相应的纵、横角度域共成像点道集，可以看到在盐丘左侧区域，这些绕射点同相轴都较为清晰，且主要体现在小角度区间。而在右侧区域，尤其右侧四个点的浅层区域，即对应盐丘区域其收到噪声影响较其他区域明显更为严重。同时在盐丘下方深层的绕射点，横波角度域共成像点道集其同相轴清晰度也不如左侧区域绕射点。但总体而言，在绕射点处成像精度达到预期目标。

为了进一步验证本发明弹性波角度域共成像点道集高效提取方法在具体计算耗时上的改进，对比常规做法：用传统解耦延拓方程构建波场、利用互相关成像条件成像、能流密度矢量计算传播角度。其中正传耗时，其耗时包含构建震源波场耗时和存储震源波场耗时；在反传耗时，该环节包含：构建检波逆向延拓耗时，计算成像时刻的纵、横波反射角耗时，对各网格点进行成像所用耗时，以及输出角度域共成像点道集所用耗时。形成耗时统计表如下所示:

本文中应用了具体实施例对本文的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本文的方法及其核心思想；同时，对于本领域的一般技术人员，依据本文的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本文的限制。

Claims (5)

1.一种弹性波角度域共成像点道集高效提取方法，其特征在于，步骤如下：

步骤1：以介质纵波速度参数和震源子波作为输入，利用交错网格高阶有限差分法求解声波一阶声压速度方程，同时记录各网格点最大振幅时刻作为该点成像时刻，并保留该时刻震源波场；

步骤2：以介质纵、横波速度参数和多分量地震记录作为输入，利用交错网格高阶有限差分法求解一阶速度-应力形式的无转换解耦延拓方程，获得各网格点成像时刻的纵、横波检波波场；

步骤3：利用各网格点成像时刻震源波场和检波波场分别计算对应的能流密度矢量，再根据余弦定理求取各网格点纵、横波反射角；

步骤4：对各网格点成像时刻震源波场和检波波场应用激发振幅成像条件，以各网格点反射角信息作为索引，进而提取纵、横波角度域共成像点道集。

2.根据权利要求1所述的一种弹性波角度域共成像点道集高效提取方法，其特征在于，步骤1的具体实现方法如下：

(1)根据预先给定的介质纵波速度模型，加载震源子波，利用声波一阶声压速度方程，实现地震波场的正向时间延拓，构建纵波震源波场，所述的一阶声波一阶声压速度方程如下：

其中，是矢量纵波震源波场的质点振动速度分量，P是压力分量，ρ表示密度，V_P表示纵波速度，t表示时间，x和z分别表示二维笛卡尔坐标系x和z方向；

利用交错网格有限差分法进行离散，其离散格式如下：

其中，Δt表示时间采样间隔，(Δx,Δz)表示空间采样间隔，i、j分别表示所述二维笛卡尔坐标系中x、z方向的离散点，t表示时间，表示精度为2N阶的交错网格有限差分系数；

(2)以波场初始值作为各网格点的最大振幅值，在震源波场正向时间延拓过程中，当某个网格点波场绝对值大于该网格点最大振幅值时，记录当前时刻作为该网格点的成像时刻，记录该网格点波场绝对值作为最大振幅，更新所记录网格点成像时刻质点振动速度和声压值。

3.根据权利要求1所述的一种弹性波角度域共成像点道集高效提取方法，其特征在于，步骤2的具体实现方法如下：

(1)根据预先设定的介质纵、横速度模型，将地面记录的多分量地震资料作为边值条件，利用一阶应力-速度形式的无转换解耦延拓方程，实现地震检波波场的逆时延拓，构建检波矢量波场；其中所述一阶应力-速度形式无转换解耦延拓方程为：

其中，(v_x,v_z)是总波场的速度分量，是纵波波场的质点振动速度分量，/>是横波波场的质点振动速度分量，τ^P是纵波应力场，/>是横波应力场，x和z分别表示二维笛卡尔坐标系x和z方向，(V_P,V_S)分别代表纵、横波速度；

(2)在检波波场逆时延拓过程中，提取空间上各网格点在成像时刻的质点振动速度和应力值。

4.根据权利要求1所述的一种弹性波角度域共成像点道集高效提取方法，其特征在于，步骤3的具体实现方法如下：

(1)利用步骤1保存的震源波场各网格点成像时刻质点振动速度和步骤2保存的检波波场各网格点成像时刻质点振动速度应力，计算各反射点纵、横波的能流密度矢量，具体表达式为：

其中，是各网格点的纵波震源波场能流密度矢量，通过成像时刻压力场P和质点振动速度场/>计算获得；/>是各网格点的纵波检波波场能流密度矢量，通过成像时刻纵波应力场τ^P和质点纵波振动速度场/>计算获得；是各网格点的横波检波波场能流密度矢量，通过成像时刻横波切应力场/>和质点纵波振动速度场/>计算获得；

(2)利用各网格点震源波场和检波波场的能流密度矢量，基于余弦定理计算各网格点的反射角，具体表达式为：

其中，(θ^P,θ^S)分别是纵波反射角和横波反射角，表示纵波震源波场的能流密度矢量，/>分别表示纵波检波波场和横波检波波场的能流密度矢量。纵波反射角为纵波入射波和纵波反射波夹角的一半，而转换横波反射角为纵波入射波和横波反射波夹角减去纵波入射角(纵波入射角等于纵波反射角)。

5.根据权利要求1所述的一种弹性波角度域共成像点道集高效提取方法，其特征在于，步骤4的具体实现方法如下：

对震源波场和检波波场应用激发振幅成像条件获得纵波入射-纵波反射和纵波入射-横波反射两种成像结果；对两种成像结果分别以纵、横波反射角作为索引，构建纵、横波共成像点道集，具体表达式为：

其中，空间坐标为(x,z)的网格点其单炮纵、横波角度域共成像点道集为(I_PP,I_PS)，为纵、横波角度域共成像点道集上的角度值，选取成像时刻t的纵波震源波场值/>和纵、横检波波场值/>来计算成像值，/>是高斯采样函数，σ函数方差，(θ^P,θ^S)是纵、横波反射角。