CN112305615B - 一种地震资料角度域共成像点道集提取方法及系统 - Google Patents

Abstract

本发明涉及一种地震资料角度域共成像点道集提取方法及系统。该方法包括：基于数值方法求解地震纯纵波波动方程，构建第k炮的背景地震波场；确定第i次迭代的角度域共成像点道集梯度；构建第i次迭代的角度域共成像点道集下降方向；根据下降方向和背景地震波场，得到第i次迭代的地震记录增量；估算最优化步长；根据最优化步长、第i次迭代的角度域共成像点道集下降方向和第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录；判定是否满足收敛标准；若是，则输出更新后的角度域共成像点道集；若否，则确定第i+1次迭代的角度域共成像点道集梯度。本发明能够获得高精度、高分辨率、高信噪比以及振幅保真的角度域共成像点道集。

Description

一种地震资料角度域共成像点道集提取方法及系统

技术领域

本发明涉及地震资料共成像点道集成像处理领域，特别是涉及一种地震资料角度域共成像点道集提取方法及系统。

背景技术

共成像点道集是一种由叠前偏移成像方法所获得的中间数据。共成像点道集所包含的反映地下介质速度和岩性变化的信息是进行偏移速度分析和地震振幅分析的有力工具。在偏移速度分析中，共成像点道集的剩余时差有效地反映了偏移速度场的误差，当速度场准确时，共成像点道集同相轴是水平的，当速度场不准确时，共成像点道集同相轴是弯曲的，而且弯曲的程度与速度场误差量是有关系的，通过这种关系并借助共成像点道集拉平原则进行速度场的更新，最终实现偏移速度建模。在地震振幅分析中，共成像点道集是叠前地震数据和地震属性解释之间的桥梁。作为岩性分析的基础资料，共成像点道集的质量直接决定了地震振幅分析技术能否准确的分析和描述岩性油气藏，也决定了弹性参数及储层参数反演的成败。面对当前复杂的油气勘探环境，传统的偏移距域共成像点道集因多路径问题会产生严重假象，在偏移速度分析和地震振幅分析方面有较大的局限性，为了解决这一问题，学者们提出了角度域共成像点道集。

自提出至今，角度域共成像点道集的研究受到了广泛的关注，发展出了众多提取角度域共成像点道集的方法。基于逆时偏移的角度域共成像点道集因其采用精确的地震波动理论，能够适用于任意复杂构造，且算法精度高、稳健性强，成为了当前油气勘探界应用的主流。然而，由于受速度不准、采集孔径有限、照明不均以及偏移算子不精确等因素影响，基于当前已有方法所提取的角度域共成像点道集分辨率低、振幅不均衡、低波数假象严重，这对后续的数据处理和解释造成很大影响，使其很难直接用于实际生产，也使得角度域共成像点道集的优势难以发挥。为此，必须建立一套新的、可提取高质量的角度域共成像点道集的成像方法及系统。

发明内容

本发明的目的是提供一种地震资料角度域共成像点道集提取方法及系统，能够获得高精度、高分辨率、高信噪比以及振幅保真的角度域共成像点道集，为地震叠前岩性反演及流体识别与预测提供高质量基础数据，指导勘探部署。

为实现上述目的，本发明提供了如下方案：

一种地震资料角度域共成像点道集提取方法，包括：

确定进行成像的多炮观测地震记录、纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型，构造方位角模型及观测系统参数；

针对第k炮，根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建第k炮的背景地震波场；

根据所述多炮观测地震记录和所述第k炮的背景地震波场，确定第i次迭代的角度域共成像点道集梯度；

基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向；

根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量；

估算最优化步长；

根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录；

判定是否满足收敛标准；

若是，则输出更新后的角度域共成像点道集，所述更新后的角度域共成像点道集为最终角度域共成像点道集；

若否，则根据所述多炮观测地震记录和所述背景地震波场，确定第i+1次迭代的角度域共成像点道集梯度。

可选地，所述第k炮，根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建背景地震波场，具体包括：

针对第k炮，基于所述观测系统参数得到第k炮的炮点坐标，在所述第k炮对应炮点位置设置震源子波，利用所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数，基于数值方法求解地震纯纵波波动方程，得到第k炮的背景地震波场。

可选地，所述根据所述多炮观测地震记录和所述第k炮的背景地震波场，确定第i次迭代的角度域共成像点道集梯度，具体包括：

设置当前的迭代次数i，针对第k炮，利用当前炮的观测地震记录及第i-1次迭代更新所得的预测地震记录，得到地震记录残差；

以所述地震记录残差为边值条件，基于数值方法逆时求解地震纯纵波波动方程，得到第k炮的反传地震波场；

在相同的时刻，对所述第k炮的背景地震波场和所述第k炮的反传地震波场应用局部偏移距域共成像点道集梯度计算方程，得到第k炮的单炮局部偏移距域共成像点道集梯度；

将所有所述单炮局部偏移距域共成像点道集依据观测系统位置信息进行叠加，得到第i次迭代的局部偏移距域共成像点道集梯度；

根据所述第i次迭代的局部偏移距域共成像点道集梯度采用局部偏移距域共成像点道集到角度域共成像点道集的转换关系，得到第i次迭代的角度域共成像点道集梯度。

可选地，所述基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向，具体包括：

基于所述第i次迭代的角度域共成像点道集梯度，利用最优化反演算法，得到第i次迭代的角度域共成像点道集下降方向。

可选地，所述根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量，具体包括：

针对第k炮，利用所述第k炮的观测系统位置信息，从所述第i次迭代的角度域共成像点道集下降方向中读取相同位置处的第k炮的下降方向；

根据所述第k炮的下降方向、所述第k炮的背景地震波场、所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型以及所述偏移成像参数，基于数值方法，应用散射条件，求解地震纯纵波散射波动方程，得到第k炮的预测地震波场；

根据所述第k炮的预测地震波场，确定第k炮在第i次迭代的地震记录增量。

可选地，所述根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录，具体包括：

根据所述最优化步长和所述第i次迭代的角度域共成像点道集下降方向，更新第i次迭代的角度域共成像点道集；

所述根据所述最优化步长和所述第i次迭代的地震记录增量，更新第i次迭代的预测地震记录。

可选地，所述收敛标准具体为：

其中，Relerr为迭代停止的阈值标准，Relerr选取1.0^e-3，misfit_i为第i次迭代的目标函数值，

misfit_i-1为第i-1次迭代的目标函数值，d_i为第k炮第i次迭代的预测地震记录，D_w为第k炮的观测纯波地震记录。

一种地震资料角度域共成像点道集提取系统，包括：

参数及模型确定模块，用于确定进行成像的多炮观测地震记录、纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型，构造方位角模型及观测系统参数；

背景地震波场构建模块，针对第k炮，用于根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建第k炮的背景地震波场；

角度域共成像点道集梯度确定模块，用于根据所述多炮观测地震记录和所述第k炮的背景地震波场，确定第i次迭代的角度域共成像点道集梯度；

角度域共成像点道集下降方向确定模块，用于基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向；

地震记录增量确定模块，用于根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量；

最优化步长估算模块，用于估算最优化步长；

更新模块，用于根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录；

判断模块，用于判定是否满足收敛标准；

最优化角度域共成像点道集确定模块，用于当满足收敛标准时，输出更新后的角度域共成像点道集，所述更新后的角度域共成像点道集为最终角度域共成像点道集；

当不满足收敛标准时，根据所述多炮观测地震记录和所述背景地震波场，确定第i+1次迭代的角度域共成像点道集梯度。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

1)本发明与常规地震资料角度域共成像点道集提取方法相比，本发明可以获得高精度、高分辨率、高信噪比、振幅保真的角度域共成像点道集；2)本发明以角度域共成像点道集为反演目标，通过地震波逆时偏移和地震波逆时反偏移，可以获得高质量的角度域叠加偏移剖面，该叠加偏移剖面直接反映了地下介质的反射系数信息，可以直接用于后续的资料解释，极大提高了储层预测及流体识别的精度。

附图说明

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

图1为本发明地震资料角度域共成像点道集提取方法流程图；

图2是本发明提供的二维层状介质模型，(a)纵波偏移速度模型v_p、(b)第一各向异性参数模型ε、(c)第二各向异性参数模型δ、(d)构造倾角模型θ；

图3是图2所示二维层状介质模型的局部偏移距域共成像点道集：其中，图3(a)是利用传统方法所得的局部偏移距域共成像点道集，图3(b)是利用本发明所得的局部偏移距域共成像点道集；

图4是图2所示二维层状介质模型的角度域共成像点道集：其中，图4(a)是利用传统方法所得的角度域共成像点道集，图4(b)是利用本发明所得的角度域共成像点道集；

图5是图2所示二维层状介质模型的角度域共成像点道集叠加偏移剖面：其中，图5(a)是利用图4(a)所示角度域共成像点道集叠加所得的偏移剖面，图5(b)是利用图4(b)所示角度域共成像点道集叠加所得的偏移剖面；

图6是本发明提供的各向异性Marmousi-2模型，(a)纵波偏移速度模型v_p、(b)第一各向异性参数模型ε、(c)第二各向异性参数模型δ、(d)构造倾角模型θ；

图7是图6所示各向异性Marmousi-2模型的地震记录剖面：其中，图7(a)第25炮的观测记录，图7(b)第55炮的观测记录，图7(c)第75炮的观测记录，图7(d)利用传统方法所得的第25炮的预测记录，图7(e)利用传统方法所得的第55炮的预测记录，图7(f)利用传统方法所得的第75炮的预测记录，图7(g)利用本发明所得的第25炮的预测记录，图7(h)利用本发明所得的第55炮的预测记录，图7(i)利用本发明所得的第75炮的预测记录；

图8是图6所示各向异性Marmousi-2模型的局部偏移距域共成像点道集：其中，图8(a)是利用传统方法所得的局部偏移距域共成像点道集，图8(b)是利用本发明所得的局部偏移距域共成像点道集；

图9是图6所示各向异性Marmousi-2模型的角度域共成像点道集：其中，图9(a)是利用传统方法所得的角度域共成像点道集，图9(b)是利用本发明所得的角度域共成像点道集；

图10是图6所示各向异性Marmousi-2模型的角度域共成像点道集叠加偏移剖面：其中，图10(a)是基于各向同性介质模型利用传统方法所得的叠加偏移剖面，图10(b)是基于各向同性介质模型利用本发明所得的叠加偏移剖面，图10(c)是基于各向异性介质模型利用传统方法所得的叠加偏移剖面，图10(d)是基于各向异性介质模型利用本发明所得的叠加偏移剖面。

具体实施方式

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

本发明的目的是提供一种地震资料角度域共成像点道集提取方法及系统，能够获得高精度、高分辨率、高信噪比以及振幅保真的角度域共成像点道集，为地震叠前岩性反演及流体识别与预测提供高质量基础数据，指导勘探部署。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明地震资料角度域共成像点道集提取方法流程图。如图1所示，一种地震资料角度域共成像点道集提取方法包括：

步骤101：确定进行成像的多炮观测地震记录、纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型，构造方位角模型及观测系统参数。

依据油气勘探目标及靶区的地质地球物理条件，确定观测系统参数，建立靶区在笛卡尔坐标系下的观测系统；确定进行角度域共成像点道集提取的观测的多炮地震记录D(x_r,t；x_s)、纵波偏移速度模型v_p、第一各向异性参数模型ε、第二各向异性参数模型δ、构造倾角模型θ、构造方位角模型

以及偏移成像参数；其中，x_s＝(x_s,y_s,z_s)表示震源空间位置矢量，x_r＝(x_r,y_r,z_r)表示检波点空间位置矢量，t表示波传播时间。

步骤102：针对第k炮，根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建第k炮的背景地震波场，具体包括：

针对第k炮，基于所述观测系统参数得到第k炮的炮点坐标，在所述第k炮对应炮点位置设置震源子波，利用所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数，基于数值方法求解地震纯纵波波动方程，得到第k炮的背景地震波场。第k炮为多炮中的任意一炮。

基于数值方法求解地震纯纵波波动方程，实现对该炮点地震波场的数值模拟，获得该炮(第k炮)的每一个时刻的背景地震波场u(x,t；x_s)，从而获得每一炮对应的背景地震波场。

其中所述的基于数值方法求解地震纯纵波波动方程，实现对该炮点地震波场的数值模拟，获得该炮的每一个时刻的背景地震波场u(x,t；x_s)，具体为利用偏微分方程数值方法求解：

所述方程(1)中，x＝(x,y,z)表示介质地下网格位置坐标矢量；f表示震源；x、y和z分别表示笛卡尔坐标系下的x方向、y方向和z方向；等式右端项Q具体为：

所述方程(1)和(2)中的偏微分算子H₁和H₂具体为

步骤103：根据所述多炮观测地震记录和所述第k炮的背景地震波场，确定第i次迭代的角度域共成像点道集梯度，具体包括：

步骤1031：设置当前的迭代次数i，针对第k炮，利用当前炮的观测地震记录及第i-1次迭代更新所得的预测地震记录，得到地震记录残差；

步骤1032：以所述地震记录残差为边值条件，基于数值方法逆时求解地震纯纵波波动方程，得到第k炮的反传地震波场；

步骤1033：在相同的时刻，对所述第k炮的背景地震波场和所述第k炮的反传地震波场应用局部偏移距域共成像点道集梯度计算方程，得到第k炮的单炮局部偏移距域共成像点道集梯度；

步骤1034：将所有所述单炮局部偏移距域共成像点道集依据观测系统位置信息进行叠加，得到第i次迭代的局部偏移距域共成像点道集梯度；

步骤1035：根据所述第i次迭代的局部偏移距域共成像点道集梯度采用局部偏移距域共成像点道集到角度域共成像点道集的转换关系，得到第i次迭代的角度域共成像点道集梯度。

设置当前的迭代次数i，针对每一炮(等同于第k炮)，利用该炮的观测地震记录D及第i-1次迭代更新所得的预测地震记录d_i(x_r,t；x_s)，计算地震记录残差Δd_i(x_r,t；x_s)，其计算方法为Δd_i(x_r,t；x_s)＝d_i-D；以该炮的地震记录残差Δd_i(x_r,t；x_s)为边值条件，基于数值方法逆时求解地震纯纵波波动方程，获得该炮的每一个时刻的反传地震波场u_r(x,t；x_s)；读取步骤102中计算的该炮的每一个时刻的背景地震波场u(x,t；x_s)，在相同的时刻，对该炮的背景地震波场和反传地震波场应用局部偏移距域共成像点道集梯度计算方程，获得该炮的单炮局部偏移距域共成像点道集梯度g(x,h；x_s)；进而获得每一炮对应的单炮局部偏移距域共成像点道集梯度；所有单炮局部偏移距域共成像点道集依据观测系统位置信息进行叠加构成了本次迭代的局部偏移距域共成像点道集梯度g_o(x,h)；利用局部偏移距域共成像点道集到角度域共成像点道集的转换关系，获得本次迭代的角度域共成像点道集梯度g_a(x,ξ,ψ)，也即第i次迭代的角度域共成像点道集梯度g_i(x,ξ,ψ)；其中h＝(h_x,h_y,h_z)表示局部偏移距矢量，h_x、h_y和h_z分别表示局部偏移距矢量h在笛卡尔坐标系下的x方向、y方向和z方向分量，ξ表示地震波传播反射角，ψ表示地震波传播方位角。

所述以该炮的地震记录残差Δd_i(x_r,t；x_s)为边值条件，基于数值方法逆时求解地震纯纵波波动方程，获得该炮的每一个时刻的反传地震波场u_r(x,t；x_s)，具体为利用偏微分方程数值方法求解：

所述读取步骤102中计算的该炮的每一个时刻的背景地震波场u(x,t；x_s)，在相同的时刻，对该炮的背景地震波场和反传地震波场应用局部偏移距域共成像点道集梯度计算方程，获得该炮的单炮局部偏移距域共成像点道集梯度g(x,h)，具体为：

g(x,h)＝∫∫u(x-h,t；x_s)·u_r(x+h,t；x_s)dtdx_s； (5)

所述方程(5)即为局部偏移距域共成像点道集梯度计算方程。

所述的利用局部偏移距域共成像点道集到角度域共成像点道集的转换关系，获得本次迭代的角度域共成像点道集梯度g_a(x,ξ,ψ)，具体为：

首先，对局部偏移距域共成像点道集g_o(x,h)的每一个水平空间位置进行共中心点道集抽取，获得g_o(x′,h)，其中x′＝(m_x,m_y,z)，m_x和m_y分别表示笛卡尔坐标系下x方向和y方向的共中心点坐标。

其次，对所抽取的共成像点道集g_o(x′,h)进行六维傅里叶变换，获得相应的波数域共成像点道集G_o(k′_x,k_h)，其中

表示空间波数，

表示笛卡尔坐标系下x方向共中心点波数，

表示笛卡尔坐标系下y方向共中心点波数，k_z表示笛卡尔坐标系z方向波数；

表示局部偏移距波数，

表示笛卡尔坐标系下x方向局部偏移距波数，

表示笛卡尔坐标系下y方向局部偏移距波数，

表示笛卡尔坐标系下z方向局部偏移距波数。

然后，应用地震波传播反射角ξ和方位角ψ与波数k′_x和k_h之间的关系，对波数域的共成像点道集G_o(k′_x,k_h)进行逐一投影，最终对投影后的共成像点道集进行六维傅里叶逆变换，获得角度域共成像点道集梯度g_a(x,ξ,ψ)，其中，所述的地震波传播反射角ξ和方位角ψ与波数k′_x和k_h之间的关系具体为：

和

步骤104：基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向，具体包括：

基于所述第i次迭代的角度域共成像点道集梯度，利用最优化反演算法，得到第i次迭代的角度域共成像点道集下降方向。

这里设第i次迭代的梯度为g_i(x,ξ,ψ)，第i次迭代的角度域共成像点道集下降方向为dg_i(x,ξ,ψ)。

步骤105：根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量，具体包括：

步骤1051：针对第k炮，利用所述第k炮的观测系统位置信息，从所述第i次迭代的角度域共成像点道集下降方向中读取相同位置处的第k炮的下降方向。即基于第i次迭代的下降方向dg_i(x,ξ,ψ)中读取相同位置处的该炮的下降方向dg(x,ξ,ψ；x_s)。

步骤1052：根据所述第k炮的下降方向、所述第k炮的背景地震波场、所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型以及所述偏移成像参数，基于数值方法，应用散射条件，求解地震纯纵波散射波动方程，得到第k炮的预测地震波场。

即基于该炮(第k炮)的下降方向dg(x,ξ,ψ；x_s)和该炮(第k炮)的背景地震波场u(x,t；x_s)，利用上述确定的纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型、构造方位角模型以及偏移成像参数，基于数值方法，应用散射条件，求解地震纯纵波散射波动方程，实现对该炮点预测地震波场的数值模拟，获得该炮的每一个时刻的预测地震波场u_c(x,t；x_s)。

步骤1053：根据所述第k炮的预测地震波场，确定第k炮在第i次迭代的地震记录增量。即基于该炮的观测系统信息，对该炮的预测地震波场进行采样，获得该炮的在第i次迭代的地震记录增量δd_i(x_r,t；x_s)；进而获得每一炮在第i次迭代的地震记录增量。

所述的读取步骤102中计算的该炮的背景地震波场u(x,t；x_s)，基于该炮的下降方向dg(x,ξ,ψ；x_s)和该炮的背景地震波场u(x,t；x_s)，利用上述确定的纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型、构造方位角模型以及偏移成像参数，基于数值方法，应用散射条件，求解地震纯纵波散射波动方程，实现对该炮点预测地震波场的数值模拟，获得该炮的每一个时刻的预测地震波场u_c(x,t；x_s)，其中，所述散射条件具体为：

所述方程(8)中，F_s(x,t；x_s)表示散射源；基于所述方程(8)表示的散射源方程，所述的地震纯纵波散射波动方程具体为：

所述的基于该炮的观测系统信息，对该炮的预测地震波场进行采样，获得该炮的在第i次迭代的地震记录增量δd_i(x_r,t；x_s)，具体为：

δd_i(x_r,t；x_s)＝u_c(x_r,t；x_s)； (16)

步骤106：估算最优化步长，所述最优化步长用α_i表示。

步骤107：根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录，具体包括：

步骤1071：根据所述最优化步长和所述第i次迭代的角度域共成像点道集下降方向，更新第i次迭代的角度域共成像点道集。

利用步骤106得到的最优化步长α_i及步骤104得到的下降方向dg_i(x,ξ,ψ)，更新第i次迭代的角度域共成像点道集R_i(x,ξ,ψ)＝R_i-1(x,ξ,ψ)+α_idg_i(x,ξ,ψ)，其中，R₀＝0。

步骤1072：所述根据所述最优化步长和所述第i次迭代的地震记录增量，更新第i次迭代的预测地震记录。

利用步骤106得到的最优化步长α_i及步骤105得到的地震记录增量δd_i，更新第i次迭代的预测地震记录d_i＝d_i-1+α_iδd_i，其中，d₀＝0。

步骤108：判定是否满足收敛标准。具体的，基于Huber范数，计算第i次迭代的目标函数值misfit_i，判断当前迭代是否满足收敛标准。

所述收敛标准具体为：

其中，Relerr为迭代停止的阈值标准，Relerr选取1.0^e-3，misfit_i为第i次迭代的目标函数值，

misfit_i-1为第i-1次迭代的目标函数值，d_i为第k炮第i次迭代的预测地震记录，D_w为第k炮的观测纯波地震记录。

步骤109：若满足收敛标准，则输出更新后的角度域共成像点道集R(x,ξ,ψ)，所述更新后的角度域共成像点道集为最终角度域共成像点道集；

步骤110：若不满足收敛标准，则根据所述多炮观测地震记录和所述背景地震波场，确定第i+1次迭代的角度域共成像点道集梯度，直至获得最终角度域共成像点道集。

对应于本发明的地震资料角度域共成像点道集提取系统，本发明还提供一种地震资料角度域共成像点道集提取系统，该系统包括：

参数及模型确定模块，用于确定进行成像的多炮观测地震记录、纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型，构造方位角模型及观测系统参数。

背景地震波场构建模块，针对第k炮，用于根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建第k炮的背景地震波场。

角度域共成像点道集梯度确定模块，用于根据所述多炮观测地震记录和所述第k炮的背景地震波场，确定第i次迭代的角度域共成像点道集梯度。

角度域共成像点道集下降方向确定模块，用于基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向。

地震记录增量确定模块，用于根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量。

最优化步长估算模块，用于估算最优化步长。

更新模块，用于根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录。

判断模块，用于判定是否满足收敛标准。

最优化角度域共成像点道集确定模块，用于当满足收敛标准时，输出更新后的角度域共成像点道集，所述更新后的角度域共成像点道集为最终角度域共成像点道集。

当不满足收敛标准时，根据所述多炮观测地震记录和所述背景地震波场，确定第i+1次迭代的角度域共成像点道集梯度。

实例1：

图2是二维层状介质模型，(a)纵波偏移速度模型v_p、(b)第一各向异性参数模型ε、(c)第二各向异性参数模型δ、(d)构造倾角模型θ。该模型深度为4km，横向宽度为6km。偏移所用空间网格大小为10m，共计41炮，起始炮点位于模型左端，炮点置于地表，炮间隔为150m，每炮最多641道接收，最小偏移距0km，最大偏移距3.2km，道间距10m，记录时间长度4.0s，时间步长1ms，采用主频为15Hz的雷克子波作为震源时间函数。图3是图2所示二维层状介质模型的局部偏移距域共成像点道集：其中，图3(a)是利用传统方法所得的局部偏移距域共成像点道集，图3(b)是利用本发明所得的局部偏移距域共成像点道集。从图3(a)中可以看出，传统方法所得局部偏移距域共成像点道集存在明显的子波旁瓣，而且远离零偏移距处有较强的能量残留。剖面存在比较明显的噪声，剖面的分辨率较低，振幅不均衡。从图3(b)中可以看出，本发明所得的局部偏移距域共成像点道集聚焦性更好，子波旁瓣得到压制，假象更少，也证明了本发明的可行性及有效性。图4是图2所示二维层状介质模型的角度域共成像点道集：其中，图4(a)是利用传统方法所得的角度域共成像点道集，图4(b)是利用本发明所得的角度域共成像点道集。由图4可以看出，利用本发明所得角度域共成像点道集具有更高的分辨率、信噪比及振幅均衡性，这间接证明了本发明的有效性。图5是图2所示二维层状介质模型的角度域共成像点道集叠加偏移剖面：其中，图5(a)是利用图4(a)所示角度域共成像点道集叠加所得的偏移剖面，图5(b)是利用图4(b)所示角度域共成像点道集叠加所得的偏移剖面。由图5可知，基于本发明所得的叠加偏移剖面分辨率高、信噪比高、振幅均衡。上述结果表明了本发明的正确性。

实例2：

图6是本发明提供的各向异性Marmousi-2模型，(a)纵波偏移速度模型v_p、(b)第一各向异性参数模型ε、(c)第二各向异性参数模型δ、(d)构造倾角模型θ。该模型是验证各种成像方法效果的国际标准模型之一。模型深度为5.4km，横向宽度为27.2km。偏移所用空间网格大小为10m，共计109炮，起始炮点位于模型2.8km处，炮点置于地表，炮间隔为200m，中间放炮，两边接收，每炮均为501道接收，最小偏移距0m，最大偏移距2500m，道间距10m，记录时间长度5.5s，时间步长1ms，采用主频为30Hz的雷克子波作为震源时间函数。图7是图6所示各向异性Marmousi-2模型的地震记录剖面：其中，图7(a)第25炮的观测记录，图7(b)第55炮的观测记录，图7(c)第75炮的观测记录，图7(d)利用传统方法所得的第25炮的预测记录，图7(e)利用传统方法所得的第55炮的预测记录，图7(f)利用传统方法所得的第75炮的预测记录，图7(g)利用本发明所得的第25炮的预测记录，图7(h)利用本发明所得的第55炮的预测记录，图7(i)利用本发明所得的第75炮的预测记录。由图7所示地震记录对比可知，利用传统角度域共成像点道集所得的预测地震记录与观测地震记录在同相轴的振幅和相位上有较大的差异，而利用本发明的角度域共成像点道集所得的预测地震记录与观测地震记录一致性较好，间接表明了本发明所得的角度域共成像点道集精度更好。图8是图6所示各向异性Marmousi-2模型的局部偏移距域共成像点道集：其中，图8(a)是利用传统方法所得的局部偏移距域共成像点道集，图8(b)是利用本发明所得的局部偏移距域共成像点道集。图9是图6所示各向异性Marmousi-2模型的角度域共成像点道集：其中，图9(a)是利用传统方法所得的角度域共成像点道集，图9(b)是利用本发明所得的角度域共成像点道集。由图8和图9可以看出，本发明提高了角度域共成像点道集的分辨率、信噪比及振幅随着反射角变化时的能量均衡性。图10是图6所示各向异性Marmousi-2模型的角度域共成像点道集叠加偏移剖面：其中，图10(a)是基于各向同性介质模型利用传统方法所得的叠加偏移剖面，图10(b)是基于各向同性介质模型利用本发明所得的叠加偏移剖面，图10(c)是基于各向异性介质模型利用传统方法所得的叠加偏移剖面，图10(d)是基于各向异性介质模型利用本发明所得的叠加偏移剖面。由图10对比可知，本发明所得图10(d)质量最高，其反射界面的分辨率高，振幅均衡，断层的断面清晰，断点位置准确，绕射能量收敛较好，假象较少。上述结果验证了本发明在复杂模型中的有效性。综上所述，本发明在复杂的地质地球物理模型中具有很好的可行性及实用性。

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

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

Claims (7)

1.一种地震资料角度域共成像点道集提取方法，其特征在于，包括：

确定进行成像的多炮观测地震记录、纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型，构造方位角模型及观测系统参数；

针对第k炮，根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建第k炮的背景地震波场；

根据所述多炮观测地震记录和所述第k炮的背景地震波场，确定第i次迭代的角度域共成像点道集梯度，具体包括：

设置当前的迭代次数i，针对第k炮，利用当前炮的观测地震记录及第i-1次迭代更新所得的预测地震记录，得到地震记录残差；

以所述地震记录残差为边值条件，基于数值方法逆时求解地震纯纵波波动方程，得到第k炮的反传地震波场；

在相同的时刻，对所述第k炮的背景地震波场和所述第k炮的反传地震波场应用局部偏移距域共成像点道集梯度计算方程，得到第k炮的单炮局部偏移距域共成像点道集梯度；

将所有所述单炮局部偏移距域共成像点道集依据观测系统位置信息进行叠加，得到第i次迭代的局部偏移距域共成像点道集梯度；

根据所述第i次迭代的局部偏移距域共成像点道集梯度采用局部偏移距域共成像点道集到角度域共成像点道集的转换关系，得到第i次迭代的角度域共成像点道集梯度；

基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向；

根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量；

估算最优化步长；

根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录；

判定是否满足收敛标准；

若是，则输出更新后的角度域共成像点道集，所述更新后的角度域共成像点道集为最终角度域共成像点道集；

若否，则根据所述多炮观测地震记录和所述背景地震波场，确定第i+1次迭代的角度域共成像点道集梯度。

2.根据权利要求1所述的地震资料角度域共成像点道集提取方法，其特征在于，所述第k炮，根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建背景地震波场，具体包括：

针对第k炮，基于所述观测系统参数得到第k炮的炮点坐标，在所述第k炮对应炮点位置设置震源子波，利用所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数，基于数值方法求解地震纯纵波波动方程，得到第k炮的背景地震波场。

3.根据权利要求1所述的地震资料角度域共成像点道集提取方法，其特征在于，所述基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向，具体包括：

基于所述第i次迭代的角度域共成像点道集梯度，利用最优化反演算法，得到第i次迭代的角度域共成像点道集下降方向。

4.根据权利要求1所述的地震资料角度域共成像点道集提取方法，其特征在于，所述根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量，具体包括：

针对第k炮，利用所述第k炮的观测系统位置信息，从所述第i次迭代的角度域共成像点道集下降方向中读取相同位置处的第k炮的下降方向；

根据所述第k炮的下降方向、所述第k炮的背景地震波场、所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型以及偏移成像参数，基于数值方法，应用散射条件，求解地震纯纵波散射波动方程，得到第k炮的预测地震波场；

根据所述第k炮的预测地震波场，确定第k炮在第i次迭代的地震记录增量。

5.根据权利要求1所述的地震资料角度域共成像点道集提取方法，其特征在于，所述根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录，具体包括：

根据所述最优化步长和所述第i次迭代的角度域共成像点道集下降方向，更新第i次迭代的角度域共成像点道集；

所述根据所述最优化步长和所述第i次迭代的地震记录增量，更新第i次迭代的预测地震记录。

6.根据权利要求1所述的地震资料角度域共成像点道集提取方法，其特征在于，所述收敛标准具体为：

其中，Relerr为迭代停止的阈值标准，Relerr选取1.0e^-3，misfit_i为第i次迭代的目标函数值，

misfit_i-1为第i-1次迭代的目标函数值，d_i为第k炮第i次迭代的预测地震记录，D_w为第k炮的观测纯波地震记录。

7.一种地震资料角度域共成像点道集提取系统，其特征在于，包括：

参数及模型确定模块，用于确定进行成像的多炮观测地震记录、纵波偏移速度模型、第一各向异性参数模型、第二各向异性参数模型、构造倾角模型，构造方位角模型及观测系统参数；

背景地震波场构建模块，针对第k炮，用于根据所述纵波偏移速度模型、所述第一各向异性参数模型、所述第二各向异性参数模型、所述构造倾角模型、所述构造方位角模型及所述观测系统参数基于数值方法求解地震纯纵波波动方程，构建第k炮的背景地震波场；

角度域共成像点道集梯度确定模块，用于根据所述多炮观测地震记录和所述第k炮的背景地震波场，确定第i次迭代的角度域共成像点道集梯度，具体包括：

设置当前的迭代次数i，针对第k炮，利用当前炮的观测地震记录及第i-1次迭代更新所得的预测地震记录，得到地震记录残差；

以所述地震记录残差为边值条件，基于数值方法逆时求解地震纯纵波波动方程，得到第k炮的反传地震波场；

在相同的时刻，对所述第k炮的背景地震波场和所述第k炮的反传地震波场应用局部偏移距域共成像点道集梯度计算方程，得到第k炮的单炮局部偏移距域共成像点道集梯度；

将所有所述单炮局部偏移距域共成像点道集依据观测系统位置信息进行叠加，得到第i次迭代的局部偏移距域共成像点道集梯度；

根据所述第i次迭代的局部偏移距域共成像点道集梯度采用局部偏移距域共成像点道集到角度域共成像点道集的转换关系，得到第i次迭代的角度域共成像点道集梯度；

角度域共成像点道集下降方向确定模块，用于基于所述第i次迭代的角度域共成像点道集梯度，构建第i次迭代的角度域共成像点道集下降方向；

地震记录增量确定模块，用于根据所述第i次迭代的角度域共成像点道集下降方向和所述背景地震波场，得到第i次迭代的地震记录增量；

最优化步长估算模块，用于估算最优化步长；

更新模块，用于根据所述最优化步长、所述第i次迭代的角度域共成像点道集下降方向和所述第i次迭代的地震记录增量更新角度域共成像点道集和预测的地震记录；

判断模块，用于判定是否满足收敛标准；

最优化角度域共成像点道集确定模块，用于当满足收敛标准时，输出更新后的角度域共成像点道集，所述更新后的角度域共成像点道集为最终角度域共成像点道集；

当不满足收敛标准时，根据所述多炮观测地震记录和所述背景地震波场，确定第i+1次迭代的角度域共成像点道集梯度。

Images