CN115598704A - 一种基于最小二乘逆时偏移生成保幅角道集的方法、设备及可读存储介质 - Google Patents

Abstract

本发明公开一种基于最小二乘逆时偏移生成保幅角道集的方法、设备及可读存储介质。其中，一种基于最小二乘逆时偏移生成保幅角道集的方法包括在逆时偏移中提取角道集作为初始成像结果，将初始成像结果作为输入，利用基于Kirchhoff近似的波动方程正演算子得到模拟地震数据，将所述模拟地震数据与观测数据的残差代入逆时偏移成像算子中，得到梯度图像，使用共轭梯度法迭代求解最小二乘反演问题以优化角道集。该方法使得在使用LSRTM进行复杂构造的定量解释时，成像幅度表示与角度相关的反射系数，也就是以成像道集的方式呈现，可直接用于幅度随角度变化(AVA)的地震解释中。采用LSRTM反演成像方法能够使得提取的角道集保幅性更好，可提高AVA分析精确度与可靠性。

Description

一种基于最小二乘逆时偏移生成保幅角道集的方法、设备及 可读存储介质

技术领域

本发明涉及地震勘探技术领域，特别涉及一种基于最小二乘逆时偏移生成保幅角道集的方法、设备及可读存储介质。

背景技术

地震深度偏移旨在通过观测到的地震反射数据生成地下构造的准确图像，在此过程中期望地震图像中的幅度保真，以便图像中的幅度信息可以用于对地下介质的物性参数进行反演解释。然而，实际中受到观测系统和上覆介质速度变化的影响，地震图像中表现出照明不均衡的现象。特别是复杂上覆介质下的照明可能存在剧烈的变化，比如在盐丘下的构造中。此时，就需要采用全波方程偏移方法比如逆时偏移，为了获得具有物理意义的幅度，其中一个解决的方法是将地震偏移建模为一个反问题。假设已知背景速度模型，目标是寻找使得模拟地震数据和观测地震数据之间差异最小的反射系数模型。这种求解线性反问题的成像方法，称为最小二乘逆时偏移(LSRTM)。

近年来，LSRTM的发展使得复杂构造的叠加成像质量大大提升。现有的LSRTM方法通常假设地下反射率独立于反射角。然而，与地下介质弹性性质相关的反射系数是依赖于反射角的，并且有时也依赖于方位角。叠加图像中的幅度不依赖于角度，是不同的观测角度下反射系数幅度的综合。因此，叠加图像中的幅度是不具有明确物理含义的，导致幅度随角度变化分析技术(Amplitude versus angle analysis，AVA)的不保幅问题，从而降低AVA的分析精确度与可靠性。

发明内容

基于此，有必要提供一种基于最小二乘逆时偏移生成保幅角道集的方法、设备及可读存储介质，以解决现有技术中AVA技术的不保幅，而降低AVA的分析精确度和可靠性的问题。

本发明提供的一种基于最小二乘逆时偏移生成保幅角道集的方法，包括如下步骤：

在逆时偏移中提取角道集作为初始成像结果；

将初始成像结果作为输入，利用基于Kirchhoff近似的波动方程正演算子输出得到模拟地震数据；

将所述模拟地震数据与观测数据的残差代入逆时偏移成像算子中，得到梯度图像；

使用共轭梯度法迭代求解最小二乘反演问题以优化角道集。

进一步地，采集原始共炮点地震记录；

对采集的所述原始共炮点地震记录做预处理，称处理后得到共炮地震记录为观测地震数据，记为D(x_r,x_s；t)，其中x_r表示检波点坐标，x_s表示震源坐标，t表示时间；

基于速度分析构建偏移速度v(x)，其中x＝(x,z)为空间坐标，不考虑密度，设密度ρ≡1；

通过对观测地震数据D(x_r,x_s；t)的频谱分析，构建宽频带的偏移子波函数w(t)；

以偏移子波w(t)作为常密度声波方程的震源项进行正向波场延拓，得到震源波场p_F(x；t)；

以观测地震数据D(x_r,x_s；t)作为常密度声波方程的边界条件进行逆时波场延拓，得到接收波场p_B(x；t)；

对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)，根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ，并且利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)的法向量估计地层倾角δ；

对震源波场p_F(x；t)和接收波场p_B(x；t)应用互相关成像条件，并做角度分箱操作将成像结果分配到对应的反射角θ上得到角度域共成像点道集R(x；θ)，所述角度域共成像点道集R(x；θ)为所述初始成像结果；

采用基于Kirchhoff近似的正演算子以

作为边界条件进行正向波场延拓，并在检波点位置记录地震波场得到模拟地震数据d(x_r,x_s；t)；

根据伴随状态法，将所述模拟地震数据与所述观测地震数据之间的残差d_r(x_r,x_s；t)＝d(x_r,x_s；t)-D(x_r,x_s；t)作为边界条件代入到偏移成像生成角道集的步骤中得到对应于最小二乘目标泛函

的梯度图像ΔR(x,θ)，所述梯度图像ΔR(x,θ)＝L^T(LR(x,θ)-D)；

其中，用L表示Kirchhoff正演算子，L^T表示生成角道集的偏移成像算子，D表示观测数据；

利用共轭梯度法对角道集进行迭代更新R_n(x,θ)＝R_n-1(x,θ)+αΔR_g(x,θ)，α是由共轭梯度法计算得到的更新步长，ΔR_g(x,θ)是共轭梯度方向。

进一步地，所述通过对观测地震数据D(x_r,x_s；t)的频谱分析，构建宽频带的偏移子波函数w(t)的步骤包括：

对共炮地震记录逐道沿时间方向做一维傅里叶变换，计算多道平均振幅谱；

确定地震记录的有效频带范围为[ω₁,ω₂]，设计通带为[ω₁,ω₂]的滤波器；

将其通过逆傅里叶变换转换到时间域，得到的时间序列作为偏移子波w(t)。

进一步地，所述以偏移子波w(t)作为常密度声波方程的震源项进行正向波场延拓，得到震源波场p_F(x；t)的步骤还包括：

所述震源波场的正向延拓可表示为

其中，v(x)为速度模型，x为地下位置坐标，p_F为震源波场，w(t)为震源子波，震源子波的积分作为边界条件。

进一步地，所述以观测地震数据D(x_r,x_s；t)作为常密度声波方程的边界条件进行逆时波场延拓，得到接收波场p_B(x；t)的步骤还包括：

接收波场的反向传播为

进一步地，所述对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)，根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ，并且利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)对应的法向量估计地层倾角δ的步骤包括：

对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)；

根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ:

对Poynting矢量进行空间平滑和归一化，得到

其中，Ω为围绕成像点的一个区域，ε为一个正则化参数以避免除零问题；

利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)计算地层倾角δ，设P(x,t)＝(P₁,P₂)，使用与P(x；t)垂直的向量来表示地层倾角，即

进一步地，所述对震源波场p_F(x；t)和接收波场p_B(x；t)应用互相关成像条件，并做角度分箱操作将成像结果分配到对应的反射角θ上得到角度域共成像点道集R(x；θ)的步骤包括：

应用基于波场分离的互相关成像条件和角度分箱操作得到地下结构的像：

其中，m(x,x_s)为叠加成像剖面，T_max为地震记录的总接收时间；

为波场在深度方向的傅里叶变换。

进一步地，所述采用基于Kirchhoff近似的正演算子以

作为边界条件进行正向波场延拓，并在检波点位置记录地震波场得到模拟地震数据d(x_r,x_s；t)的步骤包括：

其中包括背景波场p_F的正向传播，以

为边界条件得到散射波场p_B，记录检波点位置的散射波场，得到仿真地震数据d(x_r；t；x_s)＝p_B(x_r；t；x_s)。

本发明还提供一种计算机设备，包括存储器和处理器，所述存储器存储计算机程序，其特征在于，所述处理器执行所述计算机程序时实现上述的一种基于最小二乘逆时偏移生成保幅角道集的方法的步骤。

本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述的一种基于最小二乘逆时偏移生成保幅角道集的方法的步骤。

本发明提供的一种基于最小二乘逆时偏移生成保幅角道集的方法，在逆时偏移中提取角道集作为初始成像结果，使得在使用LSRTM进行复杂构造的定量解释时，成像幅度作为与角度相关的反射系数，也就是以成像道集的方式呈现，可直接用于幅度随角度变化(AVA)的地震解释中。采用反演成像方法LSRTM使得提取的角道集保幅性更好，提高了AVA分析精确度与可靠性，有助于提高含油气性检测、储层刻画和流体识别的准确度和精度。

附图说明

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

图1为本发明实施例的方法流程图。

图2为本发明实施例的波传播方向示意图。

图3为本发明实施例的Marmousi模型的真实速度模型(a)和偏移速度模型(b)。

图4为本发明实施例的炮点位于1.0km处接收到观测地震记录(a)和模拟地震记录(b)。

图5为本发明实施例的CDP 600处提取的角道集结果比较，(a)为采用逆时偏移提取的角道集，(b)为基于最小二乘逆时偏移迭代6次后生成的角道集。

图6为本发明实施例的叠加成像结果的比较，(a)采用逆时偏移提取的角道集叠加后成像结果，(b)基于最小二乘逆时偏移生成的角道集叠加后成像结果。

图7为本发明实施例的AVA曲线的比较。

本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

具体实施方式

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

需要说明，本发明实施例中所有方向性指示(诸如上、下、左、右、前、后……)仅用于解释在某一特定姿态(如附图所示)下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。

另外，在本发明中涉及“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。另外，全文中的“和/或”包括三个方案，以A和/或B为例，包括A技术方案、B技术方案，以及A和B同时满足的技术方案；另外，各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。

在一些实施例中，一种基于最小二乘逆时偏移生成保幅角道集的方法，包括：在逆时偏移的过程中中应用Poynting矢量法提取角道集作为初始成像结果，采用依赖于角度的Kirchhoff正演算子以角道集作为输入得到模拟地震数据，将模拟地震数据与观测数据的残差代入逆时偏移成像算子中，得到梯度图像，使用共轭梯度法迭代求解最小二乘反演问题以优化角道集。

具体地，Poynting矢量法，坡印亭矢量(Poynting vector)是指电磁场中的能流密度矢量。这里用作方向矢量，指示地震波场的传播方向。“Kirchhoff”正演算子是指通过Kirchhoff积分法推导出的波动方程正演算子。正演算子是指根据模型参数(这里即地下构造反射率模型，也即地震图像)通过波动方程有限差分模拟得到仿真地震数据的过程。

共轭梯度法(Conjugate Gradient)是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。在各种优化算法中，共轭梯度法是非常重要的一种。其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。共轭梯度法最早是由Hestenes和Stiefle提出来的，在这个基础上，Fletcher和Reeves首先提出了解非线性最优化问题的共轭梯度法。由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现在共轭梯度法已经广泛地应用于实际问题中。共轭梯度法是一个典型的共轭方向法，它的每一个搜索方向是互相共轭的，而这些搜索方向d仅仅是负梯度方向与上一次迭代的搜索方向的组合，因此，存储量少，计算方便。

本发明提供的基于最小二乘逆时偏移生成保幅角道集的方法提取的角道集保幅性更好，照明得到均衡(照明更加均衡是指图5中小角度道集和大角度道集之间的幅度变化更加连续光滑，从图7中的曲线对比中也能观察到，图6中图像左右两侧幅度得到提升改善，与中间位置的照明能量平衡，分辨率的提高是指图5和图6中波形的旁瓣得到抑制，模糊效应被部分消除)，分辨率提高，减少了由于地震波在复杂介质中传播时照明不均衡和偏移成像算子不完美所导致的AVA不保幅问题，从而大大提高了AVA分析精度与可靠性。因此，ADCIGs可为AVO/AVA分析以及叠前弹性参数反演提供必要的高质量、高保真度叠前角度域数据体，有助于提高含油气性检测、储层刻画和流体识别的准确度和精度。

如图1所示，基于最小二乘逆时偏移生成保幅角道集的方法的步骤。

具体地，S1、采集原始共炮点地震记录，对采集的原始共炮点做预处理，预处理包括切除直达波，去除折射波、鬼波或地滚波等，称处理后得到共炮地震记录为观测地震数据，记为D(x_r,x_s；t)，其中x_r表示检波点坐标，x_s表示震源坐标，t表示时间。

S2、通过速度分析构建偏移速度v(x)，其中x＝(x,z)为空间坐标，不考虑密度，设密度ρ≡1。

S3、通过对观测地震数据D(x_r,x_s；t)的频谱分析，构建宽频带的偏移子波函数w(t)。

具体步骤包括：对共炮地震记录逐道沿时间方向做一维傅里叶变换，计算多道平均振幅谱，确定地震记录的有效频带范围为[ω₁,ω₂]，设计通带为[ω₁,ω₂]的滤波器，将其通过逆傅里叶变换转换到时间域，得到的时间序列作为偏移子波w(t)。

S4、以偏移子波w(t)作为常密度声波方程的震源项进行正向波场延拓，得到震源波场p_F(x；t)。

震源波场的正向延拓描述为：

其中，v(x)为速度模型，x为地下位置坐标，p_F为震源波场，w(t)为震源子波，将震源子波的积分作为边界条件。

S5、以观测地震数据D(x_r,x_s；t)作为常密度声波方程的边界条件进行逆时波场延拓，得到接收波场p_B(x；t)。

接收波场的反向传播为

S6、对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)，根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ，并且利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)的法向量估计地层倾角δ。

S6的具体步骤包括：

对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)；

根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ:

对Poynting矢量进行空间平滑和归一化，得到

其中，Ω为围绕成像点的一个区域，ε为一个正则化参数以避免除零问题；

利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)计算地层倾角δ，设P(x,t)＝(P₁,P₂)，使用与P(x；t)垂直的向量来表示地层倾角，即

波场传播方向如图2所示。

S7、对震源波场p_F(x；t)和接收波场p_B(x；t)应用互相关成像条件，并做角度分箱操作将成像结果分配到对应的反射角θ上得到角度域共成像点道集R(x；θ)，角度域共成像点道集R(x；θ)为初始成像结果。

S7具体步骤包括：

应用基于波场分离的互相关成像条件和角度分箱操作得到地下结构的像：

其中，m(x,x_s)为叠加成像剖面，T_max为地震记录的总接收时间；

为波场在深度方向的傅里叶变换。

S8、采用基于Kirchhoff近似的正演算子以

作为边界条件进行正向波场延拓，并在检波点位置记录地震波场得到模拟地震数据d(x_r,x_s；t)；由于接收波场未知，接收波场通过观测数据或者数据残差作为波动方程的边界条件进行逆时反向传播获得。在此场景下，基于Kirchhoff近似的正演算子，其输入为地下反射率模型(即地震图像)，输出为散射波场(即仿真模拟得到的接收波场)，通过记录地表接收点位置的接收波场得到仿真地震数据。为了检验生成的地震图像的可靠性，目的是根据反射率模型通过正演模拟得到仿真地震数据，检验其与观测数据的相似性，以评价反射率模型。这是求解反问题的一般范式。正演的过程根据模型参数得到观测数据，其中观测数据未知。反演的过程根据观测数据得到模型参数，其中模型参数是要求解的内容。因此，其中计算反射角的方法更改为利用震源波场p_F(x；t)的方向矢量P_s(x；t)和地层倾角的法向量P(x；t)之间的夹角来表示。

S8具体步骤包括：

其中包括背景波场p_F的正向传播，

为边界条件得到散射波场p_B，记录检波点位置的散射波场，得到仿真地震数据d(x_r；t；x_s)＝p_B(x_r；t；x_s)。

S9、根据伴随状态法(伴随状态法是求解偏微分方程约束下的最优化问题时的一种计算目标函数梯度的方法)，将模拟地震数据与观测地震数据之间的残差d_r(x_r,x_s；t)＝d(x_r,x_s；t)-D(x_r,x_s；t)作为边界条件代入到偏移成像生成角道集的S4-S7中得到对应于最小二乘目标泛函

的梯度图像ΔR(x,θ)，梯度图像ΔR(x,θ)＝L^T(LR(x,θ)-D)；

其中，用L表示Kirchhoff正演算子，L^T表示生成角道集的偏移成像算子，D表示观测数据。

由于散射波场未知，其中步骤S6中计算反射角的方法更改为利用震源波场p_F(x；t)的方向矢量P_s(x；t)和地层倾角δ的法向量P(x；t)的夹角表示，P(x；t)可以用δ表示为P(x,t)＝(-sinδ,cosδ)，则反射角可以表示为：

S10、利用共轭梯度法对角道集进行迭代更新R_n(x,θ)＝R_n-1(x,θ)+αΔR_g(x,θ)，α是由共轭梯度法计算得到的更新步长，αR_g(x,θ)是共轭梯度方向。

S11、重复S8-S10步骤直到数据残差收敛到可以接受的程度或者达到一定的迭代次数。

在本实施例中，利用Marmousi模型进行数值实验，分别采用RTM和LSRTM提取角道集并进行比较。如图3所示，给出了真实速度、和偏移速度模型。模型的横向宽度和深度分别为7.5km和3.75km，网格间距均为10m，网格大小为750*375。震源从模型左侧50.0m处开始，炮间距75.0m，共100炮覆盖。采用双边接收观测系统，横向每个网格均布置检波器，共750个。震源和检波器的深度均为10.0m。震源选用主频为20Hz的Ricker子波。ADCIGs的角度范围为0～60°，角度采样间隔2°。

图4(a)给出了炮点位于1.0km处接收到的共炮地震记录。

首先，在逆时偏移中应用Poynting矢量法提取角道集作为初始成像结果，如图5(a)中所示；然后，不同于常规的LSRTM使用基于Born近似的正演算子以叠加图像作为输入得到模拟地震数据，本实施例采用基于Kirchhoff近似的正演算子以角道集作为输入得到模拟地震数据。

如图4(b)所示，通过共轭梯度法求解最小二乘反演问题以优化角道集，得到的结果如图5(b)中所示。比较图5(a)和5(b)中的结果，可以看到，在圆圈圈出的位置，LSRTM改善了角道集幅度的连续性和照明的均衡性，在箭头指出的位置中，原始RTM角道集中的弱信号在LSRTM的角道集中得到了增强，分辨率提升。将角道集进行叠加，得到叠加成像结果如图6中所示，图6(a)和6(b)中分别是RTM和LSRTM的叠加图像，通过对比可以看到，如箭头所指位置所示，LSRTM有效改善了成像分辨率，部分地缓解了同相轴子波旁瓣效应，照明更加均衡，噪声得到压制，成像幅度更加保真。进一步，通过图7中反射系数随角度变化曲线的比较可以看到，LSRTM得到的AVA曲线与理论曲线吻合地更好。

以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。

Claims (10)

1.一种基于最小二乘逆时偏移生成保幅角道集的方法，其特征在于，包括如下步骤：

在逆时偏移中提取角道集作为初始成像结果；

将初始成像结果作为输入，利用基于Kirchhoff近似的波动方程正演算子得到模拟地震数据；

将所述模拟地震数据与观测数据的残差代入逆时偏移成像算子中，得到梯度图像；

使用共轭梯度法迭代求解最小二乘反演问题以优化角道集。

2.根据权利要求1所述的方法，其特征在于，包括：

采集原始共炮点地震记录；

对采集的所述原始共炮点地震记录做预处理，称处理后得到的共炮地震记录为观测地震数据，记为D(x_r，x_s；t)，其中x_r表示检波点坐标，x_s表示震源坐标，t表示时间；

通过速度分析构建深度域偏移速度v(x)，其中x＝(x，z)为空间坐标，不考虑密度，设密度ρ≡1；

通过对观测地震数据D(x_r，x_s；t)的频谱分析，构建宽频带的偏移子波函数w(t)；

以偏移子波w(t)作为常密度声波方程的震源项进行正向波场延拓，得到震源波场p_F(x；t)；

以观测地震数据D(x_r，x_s；t)作为常密度声波方程的边界条件进行逆时波场延拓，得到接收波场p_B(x；t)；

对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)，根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ，并且利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)的法向量估计地层倾角δ；

对震源波场p_F(x；t)和接收波场p_B(x；t)应用互相关成像条件，并做角度分箱操作将成像结果分配到对应的反射角θ上得到角度域共成像点道集R(x；θ)，将所述角度域共成像点道集R(x；θ)作为初始成像结果；

采用基于Kirchhoff近似的波动方程正演算子以

作为边界条件进行正向波场延拓，并在检波点位置记录地震波场得到模拟地震数据d(x_r，x_s；t)；

根据伴随状态法，将所述模拟地震数据与所述观测地震数据之间的残差d_r(x_r，x_s；t)＝d(x_r，x_s；t)-D(x_r，x_s；t)作为边界条件代入到偏移成像生成角道集的步骤中得到对应于最小二乘目标泛函

的梯度图像ΔR(x，θ)，所述梯度图像ΔR(x，θ)＝L^T(LR(x，θ)-D)；

其中，用L表示Kirchhoff正演算子，L^T表示生成角道集的偏移成像算子，D表示观测数据；

利用共轭梯度法对角道集进行迭代更新R_n(x，θ)＝R_n-1(x，θ)+αΔR_g(x，θ)，α是由共轭梯度法计算得到的更新步长，ΔR_g(x，θ)是共轭梯度方向。

3.根据权利要求2所述的方法，其特征在于，所述通过对观测地震数据D(x_r，x_s；t)的频谱分析，构建宽频带的偏移子波函数w(t)的步骤包括：

对共炮地震记录逐道沿时间方向做一维傅里叶变换，计算多道平均振幅谱；

确定地震记录的有效频带范围为[ω₁，ω₂]，设计通带为[ω₁，ω₂]的滤波器；

将其通过逆傅里叶变换转换到时间域，得到的时间序列作为偏移子波w(t)。

4.根据权利要求3所述的方法，其特征在于，所述以偏移子波w(t)作为常密度声波方程的震源项进行正向波场延拓，得到震源波场p_F(x；t)的步骤还包括：震源波场的正向延拓描述为

其中，v(x)为速度模型，x为地下位置坐标，p_F为震源波场，w(t)为震源子波，将震源子波的积分作为边界条件。

5.根据权利要求4所述的方法，其特征在于，所述以观测地震数据D(x_r，x_s；t)作为常密度声波方程的边界条件进行逆时波场延拓，得到接收波场p_B(x；t)的步骤还包括：

接收波场的反向传播描述为

6.根据权利要求5所述的方法，其特征在于，所述对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)，根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ，并且利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)的法向量估计地层倾角δ的步骤包括：

对震源波场p_F(x；t)和接收波场p_B(x；t)分别应用Poynting矢量法估计波场传播方向的方向矢量P_s(x；t)和P_r(x；t)；

根据P_s(x；t)与P_r(x；t)的夹角计算反射角θ：

对Poynting矢量进行空间平滑和归一化，得到

其中，Ω为围绕成像点的一个区域，ε为一个正则化参数以避免除零问题；

利用P_s(x；t)和P_r(x；t)的和矢量P(x；t)计算地层倾角δ，设P(x，t)＝(P₁，P₂)，使用与P(x；t)垂直的向量来表示地层倾角，即

7.根据权利要求6所述的方法，其特征在于，所述对震源波场p_F(x；t)和接收波场p_B(x；t)应用互相关成像条件，并做角度分箱操作将成像结果分配到对应的反射角θ上得到角度域共成像点道集R(x；θ)的步骤包括：

应用基于波场分离的互相关成像条件和角度分箱操作得到地下结构的像：

其中，m(x，x_s)为叠加成像剖面，T_max为地震记录的总接收时间；

为波场在深度方向的傅里叶变换。

8.根据权利要求7所述的方法，其特征在于，所述采用基于Kirchhoff近似的正演算子以

作为边界条件进行正向波场延拓，并在检波点位置记录地震波场得到模拟地震数据d(x_r，x_s；t)的步骤包括：

其中包括背景波场p_F的正向传播，以

为边界条件得到散射波场p_B，记录检波点位置的散射波场，得到仿真地震数据d(x_r；t；x_s)＝p_B(x_r；t；x_s)。

9.一种计算机设备，包括存储器和处理器，所述存储器存储计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至8任意一项所述的一种基于最小二乘逆时偏移生成保幅角道集的方法的步骤。

10.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现权利要求1至8任意一项所述的一种基于最小二乘逆时偏移生成保幅角道集的方法的步骤。

Images