CN110007340B - 基于角度域直接包络反演的盐丘速度密度估计方法 - Google Patents

Abstract

本发明公开了基于角度域直接包络反演的盐丘速度密度估计方法，具体为：S1、设置炮点和检波点，采集得到观测数据；S2、构建矩形网格地质模型；S3、给定直接包络反演阶段待优化的目标函数J_e(v(x),ρ(x))；S4、给定用于速度长波长速度分量反演的散射角范围；S5、通过直接包络反演和迭代优化算法对目标函数J_e(v(x),ρ(x))进行优化，得到速度模型的长波长分量v_l(x)；S6、给定下一阶段新的待优化目标函数J(v(x),ρ(x))；S7、通过最速下降法对目标函数J(v(x),ρ(x))进行优化，得到精细的速度结构v_f(x)；S8、给定用于密度反演的散射角范围；S9、通过最速下降法对目标函数进行进一步优化，得到精细的密度结构ρ_f(x)。本方法能够使用直接包络中新的Fréchet导数进行能量传播。

Description

基于角度域直接包络反演的盐丘速度密度估计方法

技术领域

本发明属于地球物理勘探技术领域，具体涉及一种基于角度域直接包络反演的盐丘速度密度估计方法。

背景技术

石油和天然气是事关国家经济发展及国家安全的重要战略资源，随着石油工业的发展和勘探开发的不断深入，从构造勘探阶段逐渐走向岩性勘探阶段，勘探难度逐渐加大。为了顺应这种要求，全波形反演迅速发展起来，并成为当今地球物理界的研究热点。全波形反演是进行地层参数估计的有效方法，不但能够为地震成像提供高分辨率高精度的速度模型，而且具有同时反演各种不同参数的能力。

传统的全波形反演基于弱散射理论假设，导致该方法对初始模型的依赖性极强。如果初始模型与真实模型相差太大，则容易陷入局部极值。为了克服这一问题，一些学者选择采用全局优化的方法。然而，全局优化需要在整个模型空间中进行搜索，最终确定一个最优的模型参数，计算量相当大。虽然近年来已有不少改进算法，但是巨大的计算量仍然是阻碍全局优化发展的主要问题，使得该类方法不能在高维问题中得到广泛应用。全波形反演的重心仍然是局部优化方法。在这一领域，学者们研究和发展了多种方法，如多尺度反演方法，层析全波形反演方法，自适应全波形反演方法等等。这些方法通过使用不同的波场信息来构建目标函数从而降低反演的非线性程度，可以在一定程度上缓解局部极值问题。然而这些方法仍然是基于弱散射假设，当模型参数中包含强散射扰动，如大尺度的盐丘模型时，这些方法将失去其功效。针对大尺度的盐丘模型，工业界通常采用成像，拾取边界，再对盐丘进行填充的策略来构建盐丘体，然而这一技术需要大量的人工干预。近年来，一些新的技术逐渐被提出，来对盐丘内部的速度参数进行估计。如在自适应全波形反演中引入全变差方法进行约束，以及在反演中定义最优传输距离，都取得了不错的效果。

上述方法都是针对速度估计，在所有的地层属性中，密度是一种非常重要的属性，对于地震解释和岩性分析具有重要意义，因此对于密度参数的准确估计是很有必要的。然而，虽然全波形反演具有对各种地层参数进行估计的能力，准确的密度估计却是非常困难的，这是因为速度和密度对于波场的响应相互耦合，密度反演受到速度反演的严重影响，而对于包含大尺度盐丘体的强散射模型，这两者之间的相互影响变得更加严重。

发明内容

本发明的目的是提供一种基于角度域直接包络反演的盐丘速度密度估计方法，该方法能够使用直接包络中新的Fréchet导数进行能量传播，并在角度域对速度和密度进行解耦，实现对盐丘速度和密度的同时反演。

本发明所采用的技术方案是，基于角度域直接包络反演的盐丘速度密度估计方法，具体按照以下步骤实施：

步骤1、设置炮点和检波点，采集得到观测数据u_obs(t,x_r；x_s)，其中t表示时间变量；x_r,x_s分别表示检波点和炮点的位置；

步骤2、构建矩形网格地质模型，根据技术指标及工区要求设定模型大小Nx×Nz，并设定正演模拟的空间采样间隔、时间采样间隔dt、最大采样时间Nt，所述空间采样间隔包括横向采样间隔dx和纵向采样间隔dz；

步骤3、给定速度模型v(x)和密度模型ρ(x)的初始速度模型v₀(x)和初始密度模型ρ₀(x)，给定直接包络反演阶段待优化的目标函数J_e(v(x),ρ(x))；

步骤4、给定用于速度长波长速度分量反演的散射角范围；

步骤5、固定密度模型为初始密度模型ρ₀(x)，通过直接包络反演和迭代优化算法对目标函数J_e(v(x),ρ(x))进行优化，得到速度模型的长波长分量v_l(x)；

步骤6、通过Gardner公式得到密度模型的长波长分量ρ_l(x)，并给定下一阶段新的待优化目标函数J(v(x),ρ(x))；

步骤7、固定密度模型为长波长分量ρ_l(x)，通过最速下降法对目标函数J(v(x),ρ(x))进行优化，得到精细的速度结构v_f(x)；

步骤8、给定用于密度反演的散射角范围；

步骤9、固定速度模型v_f(x)，通过最速下降法对目标函数进行进一步优化，得到精细的密度结构ρ_f(x)。

本发明的特点还在于，

步骤3中，设定目标函数J_e(v(x),ρ(x))为观测数据包络平方与计算数据包络平方之间残差的二范数，具有如下形式：

式(1)中，e_obs(t)和e_cal(t)分别表示观测波场u_obs(t,x；x_s)和计算波场u_cal(t,x；x_s)的包络，其中计算波场通过有限差分正演模拟计算得到；r(t)为包络平方残差；上标*表示取共轭；T为波场最大观测时间；

表示关于炮点和检波点进行求和运算。

步骤3中，所述计算波场u_cal(t,x；x_s)满足如下声波方程：

式(2)中，v(x)和ρ(x)分别为空间位置x处的介质速度和密度值；δ(x-x_s)f(t)表示位于x_s处的震源函数，信号包络可通过希尔伯特变换得到，定义如下：

式(4)中H{u(t)}表示信号u(t)的希尔伯特变换。

步骤4中，散射角通过坡印廷矢量计算得到，形式如下：

式(5)中，Θ表示入射波场与散射波场之间的夹角，即散射角；P_i和P_g分别表示入射波场的坡印廷矢量和散射波场的坡印廷矢量，P_i和P_g通过如下计算得到：

式(6)中，

和

分别表示来自于炮点的正传波场和来自于检波点的反传波场，正传波场的源信号即公式(2)中所示的震源信号；反传波场通过采用新的Fréchet导数进行能量传播得到；

在进行速度长波长分量反演阶段，选取散射角范围依次为0°～30°和0°～60°。

步骤5的具体实施方式如下：

直接包络反演针对强散射问题，直接基于能量散射计算新的Fréchet导数：

式(7)中，

表示关于波场包络的新的Fréchet导数，因此直接包络反演方法利用包络扰动来反演模型参数，这样得到的速度梯度为：

式(8)中，G_e表示能量传播的格林函数，利用G_e对包络残差进行反传，得到公式(6)中的反传信号

引入角度信息，得到基于角度域直接包络反演的速度梯度公式：

式(9)中，F_v(θ_i,θ_g)为对应于速度反演的角度滤波器，控制散射角在选定的范围之内，利用该公式计算速度梯度进而实现速度更新；

利用公式(7)～(9)计算出速度梯度之后，利用如下准则对速度进行迭代更新：

vⁿ⁺¹(x)＝vⁿ(x)-α_nζ_v ⁿ(x) (10)

式(10)中，vⁿ(x),vⁿ⁺¹(x)分别为第n步和第n+1步迭代的速度值；α_n为迭代步长，为一个正常数；

通过上述方法计算速度梯度并进行迭代，在迭代收敛之后可以得到速度模型的长波长分量v_l(x)。

步骤6中，得到长波长分量v_l(x)和长波长分量ρ_l(x)之后，下一阶段目标函数J(v(x),ρ(x))定义为观测波场与计算波场残差的二范数，形式如下：

式(12)中，δu表示观测波场与计算波场之间的残差。

步骤7中，对目标函数J(v(x),ρ(x))进行优化的这一阶段速度模型梯度公式如下：

式(13)中，

和

分别表示来自于炮点的正传波场和来自于检波点的反传波场，其中正传波场源信号为公式(2)中所述的震源信号，反传波场源信号为公式(12)中的残差波场δu；

和

分别为关于

和

的时间导数。

步骤8中，选取的散射角范围为120°～180°。

步骤9中，涉及到的密度参数的梯度表达式为：

式(14)中，F_ρ(Θ)为密度反演的角度滤波器。

本发明一种基于角度域直接包络反演的盐丘速度密度估计方法的有益效果是：首先，从线性初始模型出发，固定密度参数，利用小散射角信息和直接包络反演方法得到盐丘长波长速度模型；其次，基于长波长速度模型通过Gardner经验公式得到长波长密度模型；再次，基于长波长速度模型和长波长密度模型，利用包含全散射角信息的传统全波形反演得到精细速度结构；最后，基于精细速度结构和长波长密度模型，利用大散射角信息和传统全波形反演得到精细密度结构。

与常用的反演方法相比，该方法更适用于强散射介质参数的反演，并且能够同时对速度参数和密度参数进行准确估计。该方法具有上述优势的原因在于：首先，因为使用了直接包络反演方法，采用新的Fréchet导数进行能量传播，避免了传统反演的弱散射近似；其次，因为使用了角度域分解，针对速度和密度的散射特性对这两者在角度域进行解耦，所以能够避免速度和密度的相互影响，实现同时反演。

附图说明

图1是本发明方法的流程图；

图2是目标区真实速度模型；

图3是目标区真实密度模型；

图4是线性初始速度模型；

图5是线性初始密度模型；

图6是使用传统方法反演得到的速度长波长分量；

图7是使用本发明所提出方法反演得到的速度长波长分量；

图8是反演得到的精细速度模型；

图9是反演得到的精细密度模型。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明基于角度域直接包络反演的盐丘速度密度估计方法，如图1所示，具体步骤分别如下：

步骤1、设置炮点和检波点，对原始地震数据进行采集得到观测数据u_obs(t,x_r；x_s)，并对采集到的地震数据进行常规预处理，包括动静校正、噪声消除等处理后得到叠前炮集地震数据d_obs(t,x_r；x_s)，其中t表示时间变量，x_r,x_s分别为检波点和炮点的位置；

步骤2、构建矩形网格地质模型，根据技术指标及工区要求设定模型大小Nx×Nz，并设定正演模拟的空间采样间隔(包括横向采样间隔dx和纵向采样间隔dz)、时间采样间隔dt、最大采样时间Nt。实际中可以根据地震数据有效频带范围、地震记录时间长度、地震数据观测系统确定相关参数；选定网格参数的标准是使得基于该网格进行有限差分正演模拟时，可以不仅满足稳定性条件而且有效压制数值频散。

步骤3、给定速度模型v(x)和密度模型ρ(x)的初始速度模型v₀(x)和初始密度模型ρ₀(x)，给定直接包络反演阶段待优化的目标函数J_e(v(x),ρ(x))；速度模型v(x)和密度模型ρ(x)的初始模型v₀(x)和ρ₀(x)是对真实速度模型和真实密度模型的一个粗略估计，一般可设定为横向均匀纵向线性增加的线性初始模型；目标函数J_e(v(x),ρ(x))定义了观测数据与计算数据之间的匹配程度，对目标函数J_e(v(x),ρ(x))进行优化以实现对模型参数的更新。本发明中，在直接包络反演阶段设定目标函数J_e(v(x),ρ(x))为观测数据包络平方与计算数据包络平方之间残差的二范数，具有如下形式：

式(1)中，e_obs(t)和e_cal(t)分别表示观测波场u_obs(t,x；x_s)和计算波场u_cal(t,x；x_s)的包络，其中计算波场通过有限差分正演模拟计算得到；r(t)为包络平方残差；上标*表示取共轭；T为波场最大观测时间；

表示关于炮点和检波点进行求和运算。式(1)中计算波场u_cal(t,x；x_s)满足如下声波方程：

式(2)中，v(x)和ρ(x)分别为空间位置x处的介质速度和密度值；δ(x-x_s)f(t)表示位于x_s处的震源信号。信号包络可通过希尔伯特变换得到，定义如下：

式(4)中H{u(t)}表示信号u(t)的希尔伯特变换。

上述目标函数J_e(v(x),ρ(x))中使用信号的包络信息，是因为信号包络中含有丰富的低频分量，有助于模型长波长分量的求解。

步骤4、给定用于速度长波长速度分量反演的散射角范围；速度和密度对于波场响应存在耦合，速度扰动同时产生前向散射和后向散射，而密度扰动只产生后向散射。前向散射对应于小散射角信息和介质长波长分量，后向散射对应于大散射角信息和介质短波长分量，即精细结构。可以通过使用不同的角度信息实现对速度和密度的解耦。本发明中的散射角通过坡印廷矢量计算得到，形式如下：

式(5)中，Θ表示入射波场与散射波场之间的夹角，即散射角；P_i和P_g分别表示入射波场的坡印廷矢量和散射波场的坡印廷矢量，P_i和P_g分别通过如下计算得到：

式(6)中，

和

分别表示来自于炮点的正传波场和来自于检波点的反传波场。正传波场的源信号即公式(2)中所示的震源信号；反传波场通过采用新的Fréchet导数进行能量传播得到。

计算出散射角之后，可以为速度反演和密度反演分配不同的角度信息以克服反演中的窜扰和相互影响。本发明中，在进行速度长波长分量反演阶段，选取散射角范围依次为0°～30°和0°～60°。

步骤5、固定密度模型为初始密度模型ρ₀(x)，通过直接包络反演和迭代优化算法对目标函数进行优化，得到速度模型的长波长分量v_l(x)；基于公式(1)所定义的目标函数J_e(v(x),ρ(x))进行直接包络反演，求取速度模型的长波长分量v_l(x)。传统反演方法基于弱散射近似并使用链式法则计算目标函数J(v(x),ρ(x))关于模型参数的导数，本发明中的直接包络反演针对强散射问题，直接基于能量散射计算新的Fréchet导数：

式(7)中，

表示关于波场包络的新的Fréchet导数，因此直接包络反演方法利用包络扰动来反演模型参数，这样得到的速度梯度为：

式(8)中，G_e表示能量传播的格林函数。利用G_e对包络残差进行反传，可以得到公式(6)中的反传信号

引入角度信息，可以得到基于角度域直接包络反演的速度梯度公式：

式(9)中，F_v(θ_i,θ_g)为对应于速度反演的角度滤波器，控制散射角在选定的范围之内，利用该公式计算速度梯度进而实现速度更新。

利用上述公式计算出速度梯度之后，利用如下准则对速度进行迭代更新：

vⁿ⁺¹(x)＝vⁿ(x)-α_nζ_v ⁿ(x) (10)

式(10)中，vⁿ(x),vⁿ⁺¹(x)分别为第n步和第n+1步迭代的速度值；α_n为迭代步长，一般可以取为一个小的正常数。

通过上述方法计算速度梯度并进行迭代，在迭代收敛之后可以得到速度模型长波长分量v_l(x)。

步骤6、通过Gardner公式得到密度模型的长波长分量ρ_l(x)，并给定下一阶段新的待优化目标函数J(v(x),ρ(x))；基于上述步骤5得到的长波长分量v_l(x)，使用Gardner公式计算长波长密度模型。密度模型与速度模型具有如下关系：

得到长波长分量v_l(x)和长波长分量ρ_l(x)之后，下一阶段进行精细结构反演，目标函数定义为观测波场与计算波场残差的二范数，形式如下：

式(12)中，δu表示观测波场与计算波场之间的残差。

步骤7、固定密度模型为长波长分量ρ_l(x)，通过最速下降法对目标函数J(v(x),ρ(x))进行优化，得到精细的速度结构v_f(x)；精细结构反演阶段，根据公式(12)所定义的目标函数，对模型参数进行迭代优化。这一阶段速度模型梯度公式如下：

式(13)中，

和

分别表示来自于炮点的正传波场和来自于检波点的反传波场，其中正传波场源信号为公式(2)中所述震源信号，反传波场源信号为公式(12)中的残差波场δu；

和

分别为关于

和

的时间导数。在这一阶段，首先将密度固定为由公式(11)得到的长波长分量，根据公式(13)计算得到的ζ_v(x)对速度进行迭代更新，更新准则仍然采用公式(10)中所述准则。迭代收敛之后可以得到高精度的速度模型。

步骤8、由于密度扰动只产生后向散射，对应于短波长密度分量即精细密度结构和大散射角信息，为了减小速度对密度的影响，仅选取大散射角信息进行密度反演。本发明关于密度反演选取的散射角范围为120°～180°。

步骤9、固定速度模型v_f(x)，通过最速下降法对目标函数进行进一步优化，得到精细的密度结构ρ_f(x)；引入角度信息之后，密度参数的梯度表达式为：

式(14)中，F_ρ(Θ)为密度反演的角度滤波器，将角度信息限定在选定的范围之内。根据上述公式对密度进行迭代更新，迭代收敛之后可以得到高精度的密度模型ρ_f(x)。

模型算例

本发明的具体实施过程被应用于某强散射盐丘地质模型。该模型横向共有6600米，纵向共有1600米。模型中包含大范围的高速盐丘体，而背景介质呈现低速分布，表现为典型的强散射特性。目标区的速度模型和密度模型分别如图2和3所示。

速度模型和密度模型的初始模型设定为线性初始模型，分别如图4和5所示，其中速度模型从浅到深逐渐从1500m/s线性增长到4000m/s，密度模型从浅到深逐渐从1500kg/m³线性增长到2000kg/m³。

目标区速度模型长波长分量反演对比结果如图6和图7所示。其中图6为采用传统反演方法得到的长波长速度分量；图7是经过本发明所提出的方法的步骤1到步骤5进行反演得到的目标区速度模型的长波长分量。长波长分量即大尺度的背景模型，在反演过程中，长波长分量的准确获取是整个参数反演的关键，因为它保证了迭代能够收敛到正确的全局最优解以及后续精细结构的进一步反演。而在涉及速度密度双参数反演中，速度长波长分量的正确性又关系到密度长波长分量的准确性。由图6和图7所示结果可以看出，传统方法不能很好的反演出盐丘的整体结构，而本发明中所提出的方法很好得给出了盐丘的轮廓以及速度分布。

基于上述所获得的长波长速度模型，使用本发明中所提出的方法反演得到的速度模型和密度模型的精细结构如图8和图9所示。其中图8所示为速度模型，图9所示为密度模型。可以看出无论是速度还是密度，都与真实模型接近，说明了本发明中所提出方法的有效性。

Claims (9)

1.基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，具体按照以下步骤实施：

步骤1、设置炮点和检波点，采集得到观测数据u_obs(t,x_r；x_s)，其中t表示时间变量；x_r,x_s分别表示检波点和炮点的位置；

步骤2、构建矩形网格地质模型，根据技术指标及工区要求设定模型大小Nx×Nz，并设定正演模拟的空间采样间隔、时间采样间隔dt、最大采样时间Nt，所述空间采样间隔包括横向采样间隔dx和纵向采样间隔dz；

步骤3、给定速度模型v(x)和密度模型ρ(x)的初始速度模型v₀(x)和初始密度模型ρ₀(x)，给定直接包络反演阶段待优化的目标函数J_e(v(x),ρ(x))；

步骤4、给定用于速度长波长速度分量反演的散射角范围；

步骤5、固定密度模型为初始密度模型ρ₀(x)，通过直接包络反演和迭代优化算法对目标函数J_e(v(x),ρ(x))进行优化，得到速度模型的长波长分量v_l(x)；

步骤6、通过Gardner公式得到密度模型的长波长分量ρ_l(x)，并给定下一阶段新的待优化目标函数J(v(x),ρ(x))；

步骤7、固定密度模型为长波长分量ρ_l(x)，通过最速下降法对目标函数J(v(x),ρ(x))进行优化，得到精细的速度结构v_f(x)；

步骤8、给定用于密度反演的散射角范围；

步骤9、固定速度模型v_f(x)，通过最速下降法对目标函数进行进一步优化，得到精细的密度结构ρ_f(x)。

2.根据权利要求1所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，步骤3中，设定目标函数J_e(v(x),ρ(x))为观测数据包络平方与计算数据包络平方之间残差的二范数，具有如下形式：

式(1)中，e_obs(t)和e_cal(t)分别表示观测波场u_obs(t,x；x_s)和计算波场u_cal(t,x；x_s)的包络，其中计算波场通过有限差分正演模拟计算得到；r(t)为包络平方残差；上标*表示取共轭；T为波场最大观测时间；

表示关于炮点和检波点进行求和运算。

3.根据权利要求2所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，步骤3中，所述计算波场u_cal(t,x；x_s)满足如下声波方程：

式(2)中，v(x)和ρ(x)分别为空间位置x处的介质速度和密度值；δ(x-x_s)f(t)表示位于x_s处的震源函数，信号包络可通过希尔伯特变换得到，定义如下：

式(4)中H{u(t)}表示信号u(t)的希尔伯特变换。

4.根据权利要求3所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，步骤4中，散射角通过坡印廷矢量计算得到，形式如下：

式(5)中，Θ表示入射波场与散射波场之间的夹角，即散射角；P_i和P_g分别表示入射波场的坡印廷矢量和散射波场的坡印廷矢量，P_i和P_g通过如下计算得到：

式(6)中，

和

分别表示来自于炮点的正传波场和来自于检波点的反传波场，正传波场的源信号即公式(2)中所示的震源信号；反传波场通过采用新的Fréchet导数进行能量传播得到；

在进行速度长波长分量反演阶段，选取散射角范围依次为0°～30°和0°～60°。

5.根据权利要求4所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，所述步骤5的具体实施方式如下：

直接包络反演针对强散射问题，直接基于能量散射计算新的Fréchet导数：

式(7)中，

表示关于波场包络的新的Fréchet导数，因此直接包络反演方法利用包络扰动来反演模型参数，这样得到的速度梯度为：

式(8)中，G_e表示能量传播的格林函数，利用G_e对包络残差进行反传，得到公式(6)中的反传信号

引入角度信息，得到基于角度域直接包络反演的速度梯度公式：

式(9)中，F_v(θ_i,θ_g)为对应于速度反演的角度滤波器，控制散射角在选定的范围之内，利用该公式计算速度梯度进而实现速度更新；

利用公式(7)～(9)计算出速度梯度之后，利用如下准则对速度进行迭代更新：

vⁿ⁺¹(x)＝vⁿ(x)-α_nζ_v ⁿ(x) (10)

式(10)中，vⁿ(x),vⁿ⁺¹(x)分别为第n步和第n+1步迭代的速度值；α_n为迭代步长，为一个正常数；

通过上述方法计算速度梯度并进行迭代，在迭代收敛之后可以得到速度模型的长波长分量v_l(x)。

6.根据权利要求5所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，步骤6中，得到长波长分量v_l(x)和长波长分量ρ_l(x)之后，下一阶段目标函数J(v(x),ρ(x))定义为观测波场与计算波场残差的二范数，形式如下：

式(12)中，δu表示观测波场与计算波场之间的残差。

7.根据权利要求6所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，步骤7中，对目标函数J(v(x),ρ(x))进行优化的这一阶段速度模型梯度公式如下：

式(13)中，

和

分别表示来自于炮点的正传波场和来自于检波点的反传波场，其中正传波场源信号为公式(2)中所述的震源信号，反传波场源信号为公式(12)中的残差波场δu；

和

分别为关于

和

的时间导数。

8.根据权利要求7所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，步骤8中，选取的散射角范围为120°～180°。

9.根据权利要求8所述的基于角度域直接包络反演的盐丘速度密度估计方法，其特征在于，步骤9中，涉及到的密度参数的梯度表达式为：

式(14)中，F_ρ(Θ)为密度反演的角度滤波器。

Images