CN101545986A - 基于最大能量旅行时计算的三维积分叠前深度偏移方法 - Google Patents

Abstract

本发明涉及石油地震勘探数据处理过程中的叠前深度偏移技术。具体步骤包括：对叠前地震数据建立深度速度模型；设计好观测网格，在球坐标系下求解波动方程模拟地震波的传播过程；将计算出频率域的地震波场变换转换到时间域，然后在时间域拟合地震波场格林函数能量谱，检测出最大能量到达的时间和振幅；将球坐标系下计算的旅行时和振幅场转换到直角坐标系下；根据计算好的旅行时和振幅，完成最大能量积分法叠前深度偏移；对叠前深度偏移成像道集进行速度分析，修改深度速度模型；通过迭代，修改深度速度模型，直到偏移结果满足精度后输出最终结果。本发明提出自适应变化差分网格计算技术，实现差分网格自动随半径增加逐步变细，以保证有限差分计算精度。

Description

基于最大能量旅行时计算的三维积分叠前深度偏移方法

技术领域：

该技术涉及石油地震勘探数据处理过程中的叠前深度偏移技术，具体说是一种采用波场外推技术计算地震波旅行时和振幅以提高三维积分法叠前深度偏移精度的方法。

背景技术：

Kirchhoff积分法叠前深度偏移的核心是复杂介质情况下的地震波旅行时计算，目前国内外常用的旅行时计算方法基本上都是基于程函方程，来计算地震波的初至到达时。虽然地震波初至旅行时方法计算效率高，但是在许多复杂地质构造情况下，使用地震波初至旅行时得不到精确的地震成像。

发明内容：

发明的目的是提供一种基于最大能量旅行时计算的三维积分叠前深度偏移方法，该方法基于波动方程而不是程函方程，在地震波有效范围内计算携带最大能量的地震波旅行时和振幅，在理论和实际应用上都优于基于高频近似地震波初至的常规kirchhoff积分偏移算法。

本发明的技术方案是基于单程波动方程，利用波场外推技术计算携带最大能量传播的地震波的旅行时和振幅，在此基础上完成最大能量积分法叠前深度偏移。其具体实现步骤包括：

1、对叠前地震数据进行去噪、滤波等预处理过程，建立深度速度模型。

2、在地面设计好观测网格。在每一个网格点上使用雷克子波模拟震源，根据输入的三维深度速度模型，在球坐标系下求解波动方程模拟地震波的传播过程。

3、将计算出频率域的地震波场通过反傅立叶变换转换到时间域，然后在时间域拟合地震波场格林函数能量谱，检测出最大能量到达的时间和振幅。

4、通过坐标变换将球坐标系下计算的旅行时和振幅场转换到直角坐标系下，以满足偏移计算的要求。

5、根据计算好的旅行时和振幅，完成最大能量积分法叠前深度偏移。

6、对叠前深度偏移成像道集进行速度分析，修改深度速度模型。

7、通过迭代，修改深度速度模型，重复步骤2-5，直到偏移结果满足精度后输出最终结果。

本发明针对球坐标系的网格呈放射状的特点，提出自适应变化差分网格计算技术，实现差分网格自动随半径增加逐步变细，以保证有限差分计算精度。

为提高计算效率，本发明选择稀疏频率采样间隔下进行波场外推计算，利用球坐标系下描述的地震波传播具有物理因果性的特点，可以合理地避免时间假频对成像结果的影响。本发明采用的具体方法原理如下：

1)在三维球坐标系中波动方程的Helmholtz形式为：

其中U为地震波场，r为外推半径，θ为倾角，φ为方位角，ω为角频率，v地震波传播速度。通过推导，可以得到以下方程组：

$\frac{\∂U}{\∂r}=(-\frac{1}{r}+i\sqrt{\mathrm{\α}})U---\left(2a\right)$

其中a和b是和最大倾角有关的优化系数。α定义为： $\mathrm{\α}=-\frac{{\mathrm{\ω}}^2}{v^2}+\frac{1}{r^2}.$

方程(2b)需要用有限差分法求解，使用因式分解法，将三维问题分解成两个相对简单的二维问题，丢掉高阶交叉项后，得到最终的差分方程。

其中，

${\mathrm{\α}}_{\mathrm{\θ}}=\frac{b}{r^2\mathrm{\α\Δ}{\mathrm{\θ}}^2},$

${\mathrm{\β}}_{\mathrm{\θ}}=\frac{\mathrm{a\Δr}}{2\sqrt{\mathrm{\α}}{r}^2{\mathrm{\Δ\θ}}^2}$ 和

在频率空间域用相移法来求解(2a)，差分法求解方程(3)完成波场外推。所计算的波场包含了地震波传播的能量和旅行时信息。

2)将以上方程计算的波场反变换到时间域，在时间域使用二次多项式拟合地震波格林函数能量谱，并使用牛顿叠代法检测出最大能量到达的时间。公式如下：

$t_{\max }=t-\frac{\mathrm{dE}}{\mathrm{dt}}/\frac{d^{2}E}{{\mathrm{dt}}^2},$ 公式中E是格林函数能量谱。

3)使用下述公式，完成基于数值格林函数计算的最大能量积分法叠前深度偏移：

$R(x,z)=\frac{\∫\∫A(s,x,z)A(x,g,z)e^{i[\mathrm{\φ}(s,x,z)+\mathrm{\φ}(x,g,z)]}P[s,g,t=\mathrm{\τ}(s,x,z)+\mathrm{\τ}(x,g,z)]\mathrm{dsdg}}{\∫\∫{\left[A(s,x,z)A(x,g,z)\right]}^2\mathrm{dsdg}}$

其中，τ(s，x，z)，τ(x，g，z)表示地下网格点到炮点和检波点的旅行时。A(s，x，z)，A(x，g，z)为对应波的振幅值。φ(s，x，z)，φ(x，g，z)为对应波的初相位。公式中的复相位与实波场函数的乘积可用下面方法计算：

e^iφP(t)＝cos(φ)P(t)+sin(φ)H(P(t))其中，H(·)代表Hilbert变换。

发明的效果：本发明能够适应复杂三维介质，不需要对速度场进行任何平滑，基本不存在角度限制和计算盲区。该方法结合了波动方程叠前深度偏移和Kirchhoff叠前深度偏移两种方法的优点，在充分发挥积分法适应观测系统、计算量小、成像灵活等优点的前提下，能够完成和波动方程偏移质量相近的成像效果，因此是一种兼具精度和效率的方法。

附图说明：

图1是本方法简单流程图。

图2分别显示SEG/EAGE盐丘模型INLINE方向和CROSSLINE方向各一条测线对应的速度模型和偏移结果效果图。

图3显示胜利某三维探区本方法和国外同类引进软件叠前深度偏移结果的对比效果图。

具体实施方式：

1、理论模型测试实例：

数据情况：使用SEG/EAGE盐丘模型进行测试，该模型是一个国际上通用的测试复杂构造成像和其它三维处理技术效果的理论地质模型。本次测试选用数据C的窄方位角数据。数据大小共有6个G，1炮8线式拖缆观测，记录长度5秒，8毫秒采样，每道625个样点。

旅行时和振幅计算参数：理论震源在地表激发，根据数据和构造情况选择以下参数：共计算84条炮线，每线84炮，炮间距和线间距都是160米，每炮网格：Nx，Ny，Nz分别是241、241、120，Dx，Dy，Dz都是40米。

叠前深度偏移参数：三维叠前深度偏移一共计算200条线，线间距40米；每线有401个CDP点，CDP间距是20米；偏移深度4200米，深度采样间隔为10米。偏移效果见附图：本方法对盐丘侵入体的成像轮廓清晰，断面准确，精度很高。

2、实际数据测试实例：

数据情况：使用胜利油田某探区的实际三维数据进行测试，处理范围：共有452条测线，580个CDP点；CDP网格：25*25米，面积：164.5平方公里；采样间隔：4毫秒；记录长度：5000毫秒；数据大小：60GB。首先对叠前地震数据进行去噪、滤波等预处理，建立深度速度模型。旅行时和振幅计算参数：理论震源在地表激发，根据构造情况选择以下参数：共计算40条炮线，每线50炮，炮间距和线间距都是300米，每炮网格：Nx，Ny，Nz分别是241、241、120，Dx，Dy，Dz都是100米。

叠前深度偏移参数：三维叠前深度偏移一共计算452条线，线间距25米；每线有580个CDP点，CDP间距是25米；偏移深度8000米，深度采样间隔为10米。

对叠前深度偏移成像道集进行速度分析，修改深度速度模型，迭代三次后输出最终结果。偏移效果见附图：为方便对比，偏移结果已转换到时间域。本算法对深层古潜山内幕形态成像，明显优于国外同类成像软件取得的结果。

Claims (2)

1、基于最大能量旅行时计算的三维积分叠前深度偏移方法，其特征是基于单程波动方程，利用波场外推技术计算携带最大能量传播的地震波的旅行时和振幅，在此基础上完成最大能量积分法叠前深度偏移，其具体实现步骤包括：

1)、对叠前地震数据进行去噪、滤波等预处理过程，建立深度速度模型；

2)、在地面设计好观测网格。在每一个网格点上使用雷克子波模拟震源，根据输入的三维深度速度模型，在球坐标系下求解波动方程模拟地震波的传播过程；

3)、将计算出频率域的地震波场通过反傅立叶变换转换到时间域，然后在时间域拟合地震波场格林函数能量谱，检测出最大能量到达的时间和振幅；

4)、通过坐标变换将球坐标系下计算的旅行时和振幅场转换到直角坐标系下，以满足偏移计算的要求：

5)、根据计算好的旅行时和振幅，完成最大能量积分法叠前深度偏移；

6)、对叠前深度偏移成像道集进行速度分析，修改深度速度模型；

7)、通过迭代，修改深度速度模型，重复步骤2-5，直到偏移结果满足精度后输出最终结果。

2、根据权利要求1所述的基于最大能量旅行时计算的三维积分叠前深度偏移方法，其特征是本发明采用的具体方法如下：

A)在三维球坐标系中波动方程的Helmholtz形式为：

其中U为地震波场，r为外推半径，θ为倾角，φ为方位角，ω为角频率，v地震波传播速度，通过推导，可以得到以下方程组：

$\frac{\∂U}{\∂r}=(-\frac{1}{r}+i\sqrt{\mathrm{\α}})U---\left(2a\right)$

其中a和b是和最大倾角有关的优化系数，α定义为： $\mathrm{\α}=-\frac{{\mathrm{\ω}}^2}{v^2}+\frac{1}{r^2}.$

方程(2b)需要用有限差分法求解，使用因式分解法，将三维问题分解成两个相对简单的二维问题，丢掉高阶交叉项后，得到最终的差分方程；

其中，

${\mathrm{\α}}_{\mathrm{\θ}}=\frac{b}{r^2\mathrm{\α\Δ}{\mathrm{\θ}}^2},$

${\text{\β}}_{\mathrm{\θ}}=\frac{\mathrm{a\Δr}}{2\sqrt{\mathrm{\α}}{r}^2\mathrm{\Δ}{\mathrm{\θ}}^2}$ 和

在频率空间域用相移法来求解(2a)，差分法求解方程(3)完成波场外推；

B)将以上方程计算的波场反变换到时间域，在时间域使用二次多项式拟合地震波格林函数能量谱，并使用牛顿叠代法检测出最大能量到达的时间；公式如下：

$t_{\max }=t-\frac{\mathrm{dE}}{\mathrm{dt}}/\frac{d^{2}E}{{\mathrm{dt}}^2}$ ，公式中E是格林函数能量谱。

C)使用下述公式，完成基于数值格林函数计算的最大能量积分法叠前深度偏移：

$R(x,z)=\frac{\∫\∫A(s,x,z)A(x,g,z)e^{i[\mathrm{\φ}(s,x,z)+\mathrm{\φ}(x,g,z)]}P[s,g,t=\mathrm{\τ}(s,x,z)+\mathrm{\τ}(x,g,z)]\mathrm{dsdg}}{{\∫\∫\left[A(s,x,z)A(x,g,z)\right]}^2\mathrm{dsdg}}$

其中，τ(s，x，z)，τ(x，g，z)表示地下网格点到炮点和检波点的旅行时；A(s，x，z)，A(x，g，z)为对应波的振幅值；φ(s，x，z)，φ(x，g，z)为对应波的初相位；公式中的复相位与实波场函数的乘积可用下面方法计算：

e^iφP(t)＝cos(φ)P(t)+sin(φ)H(P(t))其中，H(·)代表Hilbert变换。

Images