CN112083493B - 一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法 - Google Patents

Abstract

本发明公开了一种三维c‑τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法。本发明通过输入速度场、实际观测炮记录、起伏地表高程和观测系统文件；根据起伏地表高程生成曲网格；将三维速度模型映射到三维c‑τ坐标系下，计算正向延拓的震源波场；在三维c‑τ坐标系下，构建基于圆锥波编码T分布的最小二乘逆时偏移的目标泛函，计算逆时延拓的检波点波场和线性近似合成地震记录，确定基于圆锥波编码T分布目标泛函的梯度公式，计算更新步长，求取三维c‑τ坐标系下的圆锥波编码多震源最小二乘逆时偏移成像结果，并变换到三维笛卡尔坐标系下，输出成像结果。本发明提高了最小二乘逆时偏移的运算效率和稳定性，为剧烈起伏地表的地震勘探提供了依据。

Description

一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移 成像方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法。

背景技术

最小二乘逆时偏移作为一种高精度的地下介质成像技术，具有良好的应用前景，但是在日益增加的勘探精度需求下，其在二维向三维发展时不仅遇到了巨额计算量的瓶颈，而且稳定性也受到了挑战，因此常规最小二乘逆时偏移技术无法满足现在高纬度勘探的需求。

发明内容

为了克服三维情况下地震资料高清成像的难题，解决传统最小二乘逆时偏移的不稳定性问题，提高最小二乘逆时偏移的运算效率，本发明提出了一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法。

为了实现上述目的，本发明采用如下技术方案：

一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，具体包括如下步骤：

步骤1：输入速度场、实际观测炮记录、起伏地表高程和观测系统文件；

步骤2：根据起伏地表高程生成曲网格；

步骤3：将三维速度模型映射到三维c-τ坐标系下；

步骤4：确定三维c-τ坐标系下的波动方程，计算三维c-τ坐标系下正向延拓的震源波场；

步骤5：构建三维c-τ坐标系下基于圆锥波编码T分布的最小二乘逆时偏移的目标泛函；

步骤6：根据三维c-τ坐标系下的波动方程，计算三维c-τ坐标系下逆时延拓的检波点波场；

步骤7：根据三维c-τ坐标系下的反偏移算子方程，计算三维c-τ坐标系下的线性近似合成地震记录；

步骤8：在三维c-τ坐标系下，确定基于圆锥波编码T分布目标泛函的梯度公式，求取梯度更新方向；

步骤9：计算更新步长，求取三维c-τ坐标系下的圆锥波编码多震源最小二乘逆时偏移成像结果；

步骤10：将三维c-τ坐标系下的成像结果变换到三维笛卡尔坐标系下，输出三维最小二乘逆时偏移成像结果。

优选地，所述步骤4中，非均质各向同性介质下的三维声波一阶方程计算公式为：

式中，ρ表示密度，t表示时间，v表示速度，p表示压力，(x,y,z)为笛卡尔坐标系中的坐标点，v_p表示纵波速度，(u,v,w)分别对应表示(x,y,z)方向上的波场；

曲线τ变换域(ξ,γ,τ)由式(2)求得：

式中，v_p0表示初始纵波速度，τ_max表示曲坐标系中τ方向的最大采样点数；

计算域(x,y,z)与曲线τ变换域(ξ,γ,τ)间的转换关系为：

其中，

式中，(ξ,γ,τ)表示曲坐标系中的坐标点，v_p0表示初始纵波速度，z₀表示起伏地表的海拔函数；

为了简化表达，令：

确定曲线τ域中的三维一阶声波方程为：

式中，f表示震源。

优选地，所述步骤5中，基于圆锥波编码T分布的最小二乘逆时偏移的目标泛函计算公式为：

式中，m表示反射系数模型，E表示目标泛函，ε表示比例参数，L表示正演算子，d_obs表示观测数据，s表示震源数据；

其中，圆锥波编码的震源数据s由式(8)求得：

式中，s₀表示一炮的震源函数，ξ_i表示ξ方向上的第i炮，γ_i表示γ方向上的第i炮，ξ₀表示参考点在ξ方向上的坐标值，γ₀表示参考点在γ方向上的坐标值，c表示射线参数，N_s表示总炮数；

圆锥波编码观测数据d_obs由式(9)求得：

式中，

表示第i次观测到的数据。

优选地，所述步骤6中，在三维环境下，将基于Born近似的速度扰动导致的波场扰动转换到曲坐标系下，得到的波动方程为：

式中，v_p表示纵波速度，δu表示x方向上的扰动波场，δv表示y方向上的扰动波场、δw表示z方向上的扰动波场；δp表示扰动压力场，p₀表示初始压力场；

基于伴随状态理论，建立三维曲坐标系下的波动方程为：

式中，p^*为p的伴随矩阵，u^*为u的伴随矩阵，v^*为v的伴随矩阵，w^*为w的伴随矩阵，p_cal表示计算压力场，p_obs表示观测压力场。

优选地，所述步骤8中，梯度算子计算公式为：

式中，g表示梯度算子，d_cal表示观测数据，ε表示比例参数；

目标泛函扰动由式(13)求得：

式中，R算子将波场限制在检波点处，R^*算子将检波点处的数据扩展到整个模型空间，U表示波场，δU表示扰动波场；

对梯度算子进行简化，简化后的梯度算子计算公式为：

基于随机优化，确定加权梯度由式(15)求得：

式中，k表示迭代次数，m表示权重系数，a表示衰减系数；

基于共轭梯度算法，确定反射系数模型计算公式为：

其中，

式中，g_c表示共轭梯度方向，α表示补偿参数，C表示前置系数。

优选地，在曲线τ变换域(ξ,γ,τ)计算公式中，令：

将曲线τ变换域(ξ,γ,τ)计算公式由式(2)变为：

化简公式(19)得到：

基于链式法则，对公式(19)进一步化简，得到：

因此，

根据公式(4)，得到：

因此，雅可比矩阵J的计算公式为：

优选地，基于伴随状态理论，建立三维曲坐标系下的波动方程，过程如下所示：

根据伴随状态理论，背景波场U₀由式(25)求得：

式中，p₀表示初始压力场，u₀表示x方向上的初始波场，v₀表示y方向上的初始波场，w₀表示z方向上的初始波场；

波场U由式(26)求得：

式中，δu表示x方向上的扰动波场，δv表示y方向上的扰动波场，δw表示z方向上的扰动波场；

基于泰勒展开得到如式(27)所示的近似关系：

式中，δv_p表示纵波扰动速度；

扰动波场δU计算公式为：

基于伴随状态理论，得到：

<L^*U^*,U>＝<U^*,LU＞ (29)

其中，

式中，L表示正传算子，U表示波场，L^*表示伴随算子，U^*表示伴随波场；

设置边界条件，结合公式(1)，得到：

因此，确定伴随方程为：

将伴随方程从笛卡尔坐标系下转换到曲坐标系下，得到三维曲坐标系下的波动方程，如式(11)所示：

本发明所带来的有益技术效果：

本发明提出了一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，通过引入曲线τ变换技术、圆锥波编码技术和T分布，弥补了传统最小二乘逆时偏移在处理三维地震资料时可实施性较低的不足；

本发明缓解了传统最小二乘逆时偏移存在不稳定性的问题，解决了三维情况下地震资料高清成像的难题，提高了最小二乘逆时偏移的运算效率和稳定性，为地表剧烈起伏区域的地震勘探提供了依据。

附图说明

图1为本发明一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法的流程图。

图2为采用本发明方法时迭代次数与模型残差的关系。

图3为采用本发明方法时迭代次数与数据残差的关系。

图4为两种方法的效果对比情况，图4(a)为两种方法关于网格大小与模型残差关系的效果对比情况，图4(b)为两种方法关于网格大小与计算时间关系的效果对比情况。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

本发明提出了一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，如图1所示，具体包括如下步骤：

步骤1：输入速度场、实际观测炮记录、起伏地表高程和观测系统文件；

步骤2：根据起伏地表高程生成曲网格；

步骤3：将三维速度模型映射到三维c-τ坐标系下，通过使用三维c-τ坐标系映射，克服了剧烈起伏地表的影响以及高速区过采样的难题，使得网格离散均匀采样；

步骤4：确定三维c-τ坐标系下的波动方程，计算三维c-τ坐标系下正向延拓的震源波场；

非均质各向同性介质下的三维声波一阶方程计算公式为：

式中，ρ表示密度，t表示时间，v表示速度，p表示压力，(x,y,z)为笛卡尔坐标系中的坐标点，v_p表示纵波速度，(u,v,w)分别对应表示(x,y,z)方向上的波场；

曲线τ变换域(ξ,γ,τ)由式(2)求得：

式中，v_p0表示初始纵波速度，τ_max表示曲坐标系中τ方向的最大采样点数；

计算域(x,y,z)与曲线τ变换域(ξ,γ,τ)间的转换关系为：

其中，

式中，(ξ,γ,τ)表示曲坐标系中的坐标点，v_p0表示初始纵波速度，z₀表示起伏地表的海拔函数；

为了简化表达，令：

确定曲线τ域中的三维一阶声波方程为：

式中，f表示震源。

其中，在曲线τ变换域(ξ,γ,τ)计算公式中，令：

将曲线τ变换域(ξ,γ,τ)计算公式由式(2)变为式(19)：

化简公式(19)得到：

基于链式法则，对公式(19)进一步化简，得到：

因此，

根据公式(4)，得到：

因此，雅可比矩阵J的计算公式为：

步骤5：构建三维c-τ坐标系下基于圆锥波编码T分布的最小二乘逆时偏移的目标泛函，计算公式为：

式中，m表示反射系数模型，E表示目标泛函，ε表示比例参数，L表示正演算子，d_obs表示观测数据，s表示震源数据；

其中，圆锥波编码的震源数据s由式(8)求得：

式中，s₀表示一炮的震源函数，ξ_i表示ξ方向上的第i炮，γ_i表示γ方向上的第i炮，ξ₀表示参考点在ξ方向上的坐标值，γ₀表示参考点在γ方向上的坐标值，c表示射线参数，N_s表示总炮数；

圆锥波编码观测数据d_obs由式(9)求得：

式中，

表示第i次观测到的数据。

步骤6：根据三维c-τ坐标系下的波动方程，计算三维c-τ坐标系下逆时延拓的检波点波场；

在三维环境下，将基于Born近似的速度扰动导致的波场扰动转换到曲坐标系下，得到波动方程为：

式中，v_p表示纵波速度，δu表示x方向上的扰动波场，δv表示y方向上的扰动波场、δw表示z方向上的扰动波场；δp表示扰动压力场，p₀表示初始压力场；

基于伴随状态理论，建立三维曲坐标系下的波动方程为：

式中，p^*为p的伴随矩阵，u^*为u的伴随矩阵，v^*为v的伴随矩阵，w^*为w的伴随矩阵，p_cal表示计算压力场，p_obs表示观测压力场。

其中，基于伴随状态理论，建立三维曲坐标系下的波动方程，过程如下所示：

根据伴随状态理论，背景波场U₀由式(25)求得：

式中，p₀表示初始压力场，u₀表示x方向上的初始波场，v₀表示y方向上的初始波场，w₀表示z方向上的初始波场；

波场U由式(26)求得：

式中，δu表示x方向上的扰动波场，δv表示y方向上的扰动波场，δw表示z方向上的扰动波场；

基于泰勒展开得到如式(27)所示的近似关系：

式中，δv_p表示纵波扰动速度；

扰动波场δU计算公式为：

基于伴随状态理论，得到：

<L^*U^*,U>＝<U^*,LU> (29)

其中，

式中，L表示正传算子，U表示波场，L^*表示伴随算子，U^*表示伴随波场；

设置边界条件为：

通过将公式(30)与公式(1)联立，计算得到：

因此，确定伴随方程为：

将伴随方程从笛卡尔坐标系下转换到曲坐标系下，得到三维曲坐标系下的波动方程，如下所示：

步骤7：根据三维c-τ坐标系下的反偏移算子方程，计算三维c-τ坐标系下的线性近似合成地震记录。

步骤8：在三维c-τ坐标系下，确定基于圆锥波编码T分布目标泛函的梯度公式，求取梯度更新方向。

梯度算子计算公式为：

式中，g表示梯度算子，d_cal表示观测数据，ε表示比例参数；

目标泛函扰动由式(13)求得：

式中，R算子将波场限制在检波点处，R^*算子将检波点处的数据扩展到整个模型空间，U表示波场，δU表示扰动波场；

对梯度算子进行简化，简化后的梯度算子计算公式如式(14)所示：

基于随机优化，确定加权梯度由式(15)求得：

式中，k表示迭代次数，m表示权重系数，a表示衰减系数，本实施例中，设置权重系数m＝10，衰减系数a＝0.4；

基于共轭梯度算法，确定反射系数模型计算公式为：

式中，共轭梯度方向g_c、补偿参数α计算公式为：

式中，C表示前置系数。

步骤9：计算更新步长，求取三维c-τ坐标系下的圆锥波编码多震源最小二乘逆时偏移成像结果。

步骤10：将三维c-τ坐标系下的成像结果变换到三维笛卡尔坐标系下，输出三维最小二乘逆时偏移成像结果。

应用实验

本发明一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，应用于起伏地表三维层状模型的数据处理，取得了理想的计算效果。

输入速度场、实际观测炮记录、起伏地表高程文件和观测系统文件，确定地表网格形态；根据起伏地表高程生成正交贴体网格，并将速度场转换到曲坐标系，在曲坐标系下进行地震波正演，获取特定时刻的波场快照，再在笛卡尔坐标系下进行地震波正演，得到与曲坐标系相同时刻下的波场快照；在曲坐标系下，利用本发明方法进行有限差分正演模拟获得炮记录，再利用Hestholm和Ruud提出的模拟方法获得炮记录，对比两种方法的差别；在曲坐标系下，输入三维层状模型的偏移速度和反射系数模型；采用本发明方法分别在曲坐标系下和笛卡尔坐标系下进行迭代，获得多次迭代后的成像结果；获取本发明方法迭代次数与模型残差的关系，以及本发明方法迭代次数与数据残差的关系；对比本发明方法与Hestholm和Ruud提出的模拟方法，针对网格大小与模型残差关系和网格大小与计算时间关系进行对比；通过对比得到，本发明在起伏地表3D模型中的成像效果及计算效率得到了明显的改善，克服了传统最小二乘逆时偏移存在不稳定性的问题，解决了三维情况下地震资料高清成像的难题，提高了最小二乘逆时偏移的运算效率和稳定性，为地表剧烈起伏区域的地震勘探提供了依据。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (5)

1.一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，其特征在于，具体包括如下步骤：

步骤1：输入速度场、实际观测炮记录、起伏地表高程和观测系统文件；

步骤2：根据起伏地表高程生成曲网格；

步骤3：将三维速度模型映射到三维c-τ坐标系下；

步骤4：根据非均质各向同性介质下的三维声波一阶方程，确定三维c-τ坐标系下的波动方程，计算三维c-τ坐标系下正向延拓的震源波场；

步骤5：构建三维c-τ坐标系下基于圆锥波编码T分布的最小二乘逆时偏移的目标泛函；

步骤6：根据三维c-τ坐标系下的波动方程，计算三维c-τ坐标系下逆时延拓的检波点波场；

步骤7：根据三维c-τ坐标系下的反偏移算子方程，计算三维c-τ坐标系下的线性近似合成地震记录；

步骤8：在三维c-τ坐标系下，确定基于圆锥波编码T分布目标泛函的梯度算子计算公式，求取梯度更新方向；

步骤9：计算更新步长，求取三维c-τ坐标系下的圆锥波编码多震源最小二乘逆时偏移成像结果；

步骤10：将三维c-τ坐标系下的成像结果变换到三维笛卡尔坐标系下，输出三维最小二乘逆时偏移成像结果；

所述步骤4中，非均质各向同性介质下的三维声波一阶方程计算公式为：

式中，ρ表示密度，t表示时间，v表示速度，p表示压力，(x,y,z)为笛卡尔坐标系中的坐标点，v_p表示纵波速度，(u,v,w)分别对应表示(x,y,z)方向上的波场；

曲线τ变换域(ξ,γ,τ)由式(2)求得：

式中，v_p0表示初始纵波速度，τ_max表示曲坐标系中τ方向的最大采样点数；

计算域(x,y,z)与曲线τ变换域(ξ,γ,τ)间的转换关系为：

其中，

式中，(ξ,γ,τ)表示曲坐标系中的坐标点，v_p0表示初始纵波速度，z₀表示起伏地表的海拔函数；

为了简化表达，令：

确定曲线c-τ域中的三维一阶声波方程为：

式中，f表示震源；

所述步骤5中，基于圆锥波编码T分布的最小二乘逆时偏移的目标泛函计算公式为：

式中，m表示反射系数模型，E表示目标泛函，ε表示比例参数，L表示正演算子，d_obs表示观测数据，s表示震源数据；

其中，圆锥波编码的震源数据s由式(8)求得：

式中，s₀表示一炮的震源函数，ξ_i表示ξ方向上的第i炮，γ_i表示γ方向上的第i炮，ξ₀表示参考点在ξ方向上的坐标值，γ₀表示参考点在γ方向上的坐标值，c表示射线参数，N_s表示总炮数；

圆锥波编码观测数据d_obs由式(9)求得：

式中，

表示第i次观测到的数据。

2.根据权利要求1所述的一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，其特征在于，所述步骤6中，在三维环境下，将基于Born近似的速度扰动导致的波场扰动转换到曲坐标系下，得到的波动方程为：

式中，v_p表示纵波速度，δu表示x方向上的扰动波场，δv表示y方向上的扰动波场、δw表示z方向上的扰动波场；δp表示扰动压力场，p₀表示初始压力场；

基于伴随状态理论，建立三维曲坐标系下的波动方程为：

式中，p^*为p的伴随矩阵，u^*为u的伴随矩阵，v^*为v的伴随矩阵，w^*为w的伴随矩阵，p_cal表示计算压力场，p_obs表示观测压力场。

3.根据权利要求2所述的一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，其特征在于，所述步骤8中，梯度算子计算公式为：

式中，g表示梯度算子，d_cal表示观测数据，ε表示比例参数；

目标泛函扰动由式(13)求得：

式中，R算子将波场限制在检波点处，R^*算子将检波点处的数据扩展到整个模型空间，U表示波场，δU表示扰动波场；

对梯度算子进行简化，简化后的梯度算子计算公式为：

基于随机优化，确定加权梯度由式(15)求得：

式中，k表示迭代次数，m表示权重系数，a表示衰减系数；

基于共轭梯度算法，确定反射系数模型计算公式为：

其中，

式中，g_c表示共轭梯度方向，α表示补偿参数，C表示前置系数。

4.根据权利要求1所述的一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，其特征在于，在曲线τ变换域(ξ,γ,τ)计算公式中，令：

将曲线τ变换域(ξ,γ,τ)计算公式由式(2)变为：

化简公式(19)得到：

基于链式法则，对公式(19)进一步化简，得到：

因此，

根据公式(4)，得到：

因此，雅可比矩阵J的计算公式为：

5.根据权利要求2所述的一种三维c-τ坐标系的圆锥波编码多震源最小二乘逆时偏移成像方法，其特征在于，基于伴随状态理论，建立三维曲坐标系下的波动方程，过程如下所示：

根据伴随状态理论，背景波场U₀由式(25)求得：

式中，p₀表示初始压力场，u₀表示x方向上的初始波场，v₀表示y方向上的初始波场，w₀表示z方向上的初始波场；

波场U由式(26)求得：

式中，δu表示x方向上的扰动波场，δv表示y方向上的扰动波场，δw表示z方向上的扰动波场；

基于泰勒展开得到如式(27)所示的近似关系：

式中，δv_p表示纵波扰动速度；

扰动波场δU计算公式为：

基于伴随状态理论，得到：

<L^*U^*,U>＝<U^*,LU> (29)

其中，

式中，L表示正传算子，U表示波场，L^*表示伴随算子，U^*表示伴随波场；

设置边界条件，结合公式(1)，得到：

因此，确定伴随方程为：

将伴随方程从笛卡尔坐标系下转换到曲坐标系下，得到三维曲坐标系下的波动方程，如式(11)所示：

Images