CN112965102A - 一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法 - Google Patents

Abstract

本说明书实施例公开了一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法。基本本申请的方案，获取初始输入数据之后，基于预先确定的包含有阻抗敏感核函数的伴随偏移算子

确定反射系数模型m_k的初始值，进而基于线性正演算子L和所述伴随偏移算子

对反射系数模型m_k进行迭代更新，并最终得到成像结果。与现有地震成像技术相比，该技术方法采用阻抗敏感核函数计算成像条件，主可以很好的避免逆时偏移中出现的背景噪音，改善成像质量；在每次迭代中自动消除梯度中的背景噪音，提高泛函梯度的空间分辨率和信噪比，进而改善总体收敛速率，节省大量计算及人工质控成本。

Description

一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法

技术领域

本说明书涉及勘探地震学领域，尤其涉及一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法。

背景技术

由于地震采集中存在有限频带限制、观测系统不规则、地下照明不均匀等一系列问题，常规的地震成像方法无法满足目前高精度地震勘探需求。通过使用合成数据拟合观测数据，最小二乘偏移技术可以克服上述部分问题，获得高质量地震剖面。但是，当前成像精度最高的逆时偏移技术存在低频噪音干扰，这使得对应的最小二乘逆时偏移收敛较慢，计算量较大，很难在计算成本较小的前提下的获得宽频带、高分辨率偏移结果。

因此，针对常规的最小二乘逆时偏移这一问题，需要进行系统分析，研发成本较低、收敛较快、反演精度较高的稳健最小二乘成像技术。

发明内容

本发明的目的在于，提供一种高效、稳健、准确的最小二乘逆时偏移成像方法。

为解决上述技术问题，本发明采用如下技术方案：

一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法，包括：

获取初始输入数据，所述初始输入数据包括：偏移速度模型v₀(x)、偏移密度模型ρ₀、震源函数子波w(t)、观测数据d^obs、预设最大迭代次数N_iter和冗余误差ε；

根据所述初始输入数据和预先确定的伴随偏移算子

确定反射系数模型的初始值m₀；

对反射系数模型m_k进行如下迭代更新，直至满足预设的迭代终止条件：

对第k次获得的反射系数m_k，进行线性正演模拟，得到合成地震数据d^syn＝Lm_k，其中L为线性正演算子，0≤k≤N_iter；

计算数据残差r＝d^syn-d^obs，根据所述数据残差和所述伴随偏移算子确定第k+1次迭代的更新梯度

根据所述更新梯度计算第k+1次的更新步长α_k+1、共轭步长β_k+1和更新方向h_k+1，对反射系数模型进行第k+1次的更新m_k+1＝m_k+α_k+1h_k+1；

当迭代停止时，输出当前的反射系数模型，作为最终成像结果；

其中，所述线性正演算子L的形式为：

其中，x为成像空间位置，x_s为震源位置，和x_r为检波点位置，t为波场传播时间，T为记录长度，v₀(x)为偏移速度模型，z₀(x)为背景阻抗模型，p(x_s,x,t)为正向传播波场，

为伴随波场，

为空间梯度算子，f(x_s,t)为震源函数，

为虚拟震源；

所述伴随偏移算子

的形式为：

其中，

为伴随波场，r(x_r,t)为数据残差，δ为狄拉克δ函数，K(x_s,x_r,x)为阻抗敏感核函数。

本说明书实施例采用的上述至少一个技术方案能够达到以下有益效果：

基本本申请的方案，获取初始输入数据之后，基于预先确定的包含有阻抗敏感核函数的伴随偏移算子

确定反射系数模型m_k的初始值，进而基于线性正演算子L和所述伴随偏移算子

对反射系数模型m_k进行迭代更新，并最终得到成像结果。与现有地震成像技术相比，该技术方法采用阻抗敏感核函数计算成像条件，主可以很好的避免逆时偏移中出现的背景噪音，改善成像质量；在每次迭代中自动消除梯度中的背景噪音，提高泛函梯度的空间分辨率和信噪比，进而改善总体收敛速率，节省大量计算及人工质控成本。

附图说明

图1为本说明书实施例所提供的流程示意图；

图2为真实速度v模型；

图3为真实密度模型ρ；

图4为偏移速度v₀模型；

图5为偏移密度模型ρ₀；

图6为速度扰动反射率模型；

图7为阻抗扰动反射率模型；

图8为常规最小二乘逆时偏移第一次迭代结果；

图9为常规最小二乘逆时偏移第十五次迭代结果；

图10为本说明书实施例所提出最小二乘偏移第一次迭代结果；

图11为本说明书实施例所提出最小二乘偏移第十五次迭代结果；

图12为四种方法局部放大对比结果；

图13为本说明书实施例与常规方法的迭代收敛速度对比示意图。

具体实施方式

为使本申请的目的、技术方案和优点更加清楚，下面将结合本申请具体实施例及相应的附图对本申请技术方案进行详细、完整地描述。显然，所描述的实施例仅是本申请一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

本说明书的实施例提供的一种基于阻抗敏感核函数最小二乘逆时偏移成像方法，如图1所示，该过程具体包括以下步骤：

S101，获取初始输入数据。

这里的输入数据包括偏移速度模型v₀、偏移密度模型ρ₀、震源函数子波w(t)、观测数据d^obs、预设最大迭代次数N_iter和冗余误差ε。

需要说明的是，在本说明书实施例中，偏移速度模型一般由地震射线类层析或全波形反演获得，观测数据由野外采集得到，震源函数可以从近偏移距数据提取得到，迭代次数一般设置为15次，冗余误差一般设置为5％。

S103，根据所述初始输入数据和预先确定的伴随偏移算子，确定反射系数模型的初始偏移结果。

具体而言，伴随偏移算子

的形式为

其中的K(x_s,x_r,x)即为阻抗敏感核函数，初始偏移结果由所有炮点和检波点的敏感核函数叠加构成，即

在该算子中，x为成像空间位置，x_s和x_r分别为震源与检波点位置，v₀(x)为输入的偏移速度模型，z₀(x)为背景阻抗模型，可通过z₀(x)＝v₀(x)ρ₀(x)计算得到，f(x_s,t)为震源函数，r(x_r,t)为数据残差，第一次迭代设置为观测数据，δ(x)为狄拉克δ函数。在计算成像条件时，采用如前所述的阻抗敏感核函数K(x_s,x_r,x)进行约束，其中，p(x_s,x,t)为正向波场，

为伴随波场。

S105，使用基于敏感核函数计算的模型梯度，对反射系数模型进行迭代更新，直至满足预设的迭代终止条件。

具体而言，进行迭代更新时首先采用使用如下的线性正演模拟算子L进行正演波场模拟：

其中，p₀(x_s,x,t)为背景波场，δp(x,x_r,t)为扰动波场，f(x_s,t)为震源函数。

为虚拟震源，并可以通过如下表达式计算得到

δlnz(x)＝δz(x)/z₀(x)，为高波数阻抗扰动项，也就是希望得到的反射系数模型m。

进而可以计算数据残差为d^syn-d^obs，当数据残差大于预设误差时或小于最大迭代次数时，计算更新梯度方向g_k和步长α_k，并根据最速下降法下降方向和步长计算共轭梯度法下降方向h_k和共轭步长β_k，进而更新反射系数场m_k+1＝m_k-α_kh_k。具体而言，本发明建立在观测地震记录d^obs与线性正演合成数据d^syn＝Lm_k的L₂范数的拟合基础上，所使用的L₂范数目标泛函如下：

在求取反射系数时，可以采用不同的优化方法，不同的优化方法，所对应的求解步骤不同，本专利使用预条件共轭梯度方法。

所对应的基于阻抗敏感核函数的最小二乘偏移成像方法的基本迭代公式如下：

m_k+1＝m_k+α_kh_k (5)

其中，m_k表示在第k次迭代中所求取的反射系数；α_k和h_k分别表示第k次迭代的更新步长和共轭梯度方向(即更新方向)，预条件共轭梯度法是在最速下降法和常规共轭梯度法基础上发展而来，具有更快的收敛速度。对于每次迭代的更新梯度而言，其取值为即在所有炮点和检波点上的敏感核函数的取值的叠加，实际上由于数据残差在不断变化，因此每次迭代时在同一点上所计算得到的敏感核函数的取值也会相应发生变化，而导致更新梯度的变化，即

具体参见公式(7)。预条件共轭梯度法中所需要的步长及更新梯度可通过以下计算得到：

g_k+1＝L^T(Lm_k+1-d^obs) (7)

h_k+1＝H^-1g_k+1+β_k+1h_k+1(9)其中，g_k为第k次迭代时的更新梯度，g_k+1和β_k+1表示第k+1次迭代的梯度和共轭步长。

S107，当迭代停止时，输出当前的反射系数模型，作为最终成像结果。

当数据残差r小于所述冗余误差ε，或者，当前迭代次数超过最大迭代次数N_iter时，迭代即终止。此时即可以输出当前的反射系数模型，作为最终成像结果。

基本本申请的方案，获取初始输入数据之后，基于预先确定的包含有阻抗敏感核函数的伴随偏移算子

确定反射系数模型m_k的初始值，进而基于线性正演算子L和所述伴随偏移算子

对反射系数模型m_k进行迭代更新，并最终得到成像结果。与现有地震成像技术相比，该技术方法采用阻抗敏感核函数计算成像条件，主可以很好的避免逆时偏移中出现的背景噪音，改善成像质量；在每次迭代中自动消除梯度中的背景噪音，提高泛函梯度的空间分辨率和信噪比，进而改善总体收敛速率，节省大量计算及人工质控成本。

以下给出本说明是实施例在模型中的实际效果说明。

将本发明所提供的方法应用于国际标准Marmousi模型成像，取得了较理想的成像效果。真实速度模型(如图2所示)、密度模型(如图3所示)、偏移速度模型(如图4示)、偏移密度模型(如图5所示)、真实速度反射率模型(如图6所示)、阻抗反射率模型(如图7所示)；采用3.2km的观测孔径，使用真实的速度密度模型，进行全波形正演模拟，得到观测地震记录；输入偏移模型和密度模型进行波场正演模拟，并反传观测地震记录，采用常规最小二乘偏移方法，得到常规逆时偏移(如图8所示)及最小二乘逆时偏移成像结果(如图9所示)。

采用本说明书实施例所提供的方案，对常规声波方程使用速度和阻抗双参数进行参数化，并使用高波数阻抗扰动模型作为反射率模型，得到Born近似下的线性正演地震记录，将该地震记录与实际观测地震数据相减，得到L₂范数框架下的数据残差，计算残差的L₂模并判断其是否小于给定冗余误差或者当前迭代是否超过预设的最大迭代次数，如果没有则求取泛函梯度、更新方向、更新步长及共轭步长，然后更新反射率模型，再次进行线性正演并计算合成记录与观测记录的L₂模数据残差，直至数据残差小于给定阈值或迭代次数大于预设最大迭代次数，输出最终的偏移成像剖面(如图10、11及12所示)。

相比常规最小二乘偏移结果(如图8和9所示)，本说明书实施例的方案采用阻抗敏感核函数计算泛函梯度，可以自动消除逆时偏移中的低频噪音，改善成像剖面分辨率和信噪比(如图11和12所示)，大大提高模型反演的收敛速度(如图13所示)。深部油藏位置反射同相轴具有更好的连续性(如图12所示)。相比第一次基于阻抗核偏移成像结果(如图10所示)，浅层偏移噪音和采集脚印得到压制，连续性得到明显提高，中深层振幅更加均衡。在前述示意图中，Distance对应于横坐标x，Depth对应于纵坐标z。

为本说明书实施例所提供的方案，为复杂地质构造及岩性解释工作提供较准确的成像数据基础，为油田高精度勘探开发提供强有力的技术支持。

对应的，本申请实施例还提供一种计算机设备，所述设备包括包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其中，所述处理器执行所述程序时实现如前所述的基于阻抗敏感核函数最小二乘逆时偏移成像方法所述的方法。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置、设备和介质类实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可，这里就不再一一赘述。

上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤或模块可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。

Claims (4)

1.一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法，包括：

获取初始输入数据，所述初始输入数据包括：偏移速度模型v₀(x)、偏移密度模型ρ₀、震源函数子波w(t)、观测数据d^obs、预设最大迭代次数N_iter和冗余误差ε；

根据所述初始输入数据和预先确定的伴随偏移算子

确定反射系数模型m_k的初始值m₀；

对反射系数模型m_k进行如下迭代更新，直至满足预设的迭代终止条件：

对第k次获得的反射系数m_k，进行线性正演模拟，得到合成地震数据d^syn＝Lm_k，其中L为线性正演算子，0≤k≤N_iter；

计算数据残差r＝d^syn-d^obs，根据该数据残差和所述伴随偏移算子确定第k+1次迭代的更新梯度

根据所述更新梯度计算第k+1次的更新步长α_k+1、共轭步长β_k+1和更新方向h_k+1，对反射系数模型进行第k+1次的更新m_k+1＝m_k+α_k+1h_k+1；

当迭代停止时，输出当前的反射系数模型，作为最终成像结果；

其中，所述线性正演算子L的形式为：

其中，x为成像空间位置，x_s为震源位置，和x_r为检波点位置，t为波场传播时间，T为记录长度，v₀(x)为偏移速度模型，z₀(x)为背景阻抗模型，p(x_s,x,t)为正向传播波场，

为伴随波场，

为空间梯度算子，f(x_s,t)为震源函数，

为虚拟震源；

所述伴随偏移算子

的形式为：

其中，

为伴随波场，r(x_r,t)为数据残差，δ为狄拉克δ函数，K(x_s,x_r,x)为阻抗敏感核函数。

2.如权利要求1所述的方法，第k+1次更新步长为

第k+1次的共轭步长为

第k+1次迭代的更新方向为h_k+1＝H^-1g_k+1+β_k+1h_k，其中，H^-1为震源照明预条件算子，<>表示点乘运算，g_k为第k次迭代时的更新梯度。

3.如权利要求1所述的方法，所述预设的迭代终止条件，包括：数据残差L₂模小于所述冗余误差ε，或者，当前迭代次数超过最大迭代次数N_iter。

4.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其中，所述处理器执行所述程序时实现如权利要求1至3任一所述的方法。

Images