CN102103216B - 二维高斯射线束叠前偏移方法 - Google Patents

Abstract

本发明涉及地球物理勘探，是二维高斯射线束叠前偏移技术，对炮集内各道地震数据进行处理后的道数据及其Hilbert变换数据，由积分变换得到离散样板表，对给定的炮点处射线角度，计算炮点射线路径和旅行时，根据旅行时利用高斯射线束动力学方程计算炮点高斯射线束动力学参数，利用运动学方程按时间间隔进行射线追踪，计算检波点射线路径和旅行时及检波点高斯射线束动力学参数，并叠加偏移结果并显示为地层剖面图像。本发明用高斯射线束替代Kirchhoff射线进行偏移，解决了射线的盲区、焦散区和多值走时问题，提高了成像效果，同时保持了射线法灵活高效的优点。

Description

二维高斯射线束叠前偏移方法

技术领域

本发明涉及地球物理勘探技术，属于地震资料叠前偏移处理技术范畴，是一种实现反射界面空间归位，恢复反射波场特征，提高地震水平分辨率和地震信号保真度的一种二维高斯射线束叠前偏移方法。

背景技术

在地球物理勘探地震资料叠前偏移处理技术范畴中，目前业界广泛使用的偏移处理方法可以分为基于波动方程积分解的Kirchhoff射线法和波动方程直接解法两大类。其中Kirchhoff射线法的优点是计算速度快，能适用于任意倾角的反射界面，对复杂观测系统有很强的适应能力；缺点是由于采用了高频近似，Kirchhoff法对复杂构造的处理能力较弱，尤其是在多值走时的处理上存在困难，因而在横、纵向速度变化大以及构造形态复杂等条件下成像质量不够理想。波动方程的直接解法提高了偏移成像精度，克服了Kirchhoff射线法对复杂构造成像能力弱的缺点，但是由于直接求解波动方程工作量巨大，这类方法的计算效率大大低于Kirchhoff方法。

对比Kirchhoff和波动方程方法，Kirchhoff法计算效率上有优势而波动方程直接解法在效果上更胜一筹，但在计算效率和精度上没有一个更加合理的平衡，使实际偏移处理总有缺憾。

发明内容

本发明目的是提供一种计算效率高而便于复杂成像的二维高斯射线束叠前偏移方法。

本发明具体步骤包括：

1)激发地震并记录地震波，取得叠前单炮地震数据；

2)对炮集内各道地震数据f进行积分变换，得到处理后的道数据b及其Hilbert变换数据

其中b由叠前单道记录f的傅立叶频谱F(ω)通过下面的积分变换公式得到：

B(ω，T_i)＝exp(|ω|T_i)F(ω)

$b(T_r,T_i)=\frac{1}{2\mathrm{\π}}\∫\mathrm{d\ωexp}(-\mathrm{i\ω}{T}_r)B(\mathrm{\ω},T_i)---\left(1\right)$

式中T_r和T_i分别代表高斯射线束复值旅行时的实部和虚部，对T_r和T_i进行离散化，取ΔT_r＝地震记录采样间隔/4，

计算得到关于b和

的离散样板表。

步骤2所述的T_r的取值范围为0-原始记录最大采样时间，T_i的取值范围为

其中ω_l为地震数据有效频带的低端，取值为60弧/秒。

3)确定炮点处射线角度范围γ₁至γ_n，和射线角度间隔Δγ；

4)对于给定的炮点处射线角度γ，利用运动学方程按一定时间间隔Δt进行射线追踪，计算炮点射线路径和旅行时；

步骤4所述的运动学射线追踪通过求解如下的射线路径方程进行：

$\frac{\mathrm{dx}}{\mathrm{dt}}=v\sin \mathrm{\σ}$

$\frac{\mathrm{dz}}{\mathrm{dt}}=v\cos \mathrm{\σ}---\left(2\right)$

$\frac{\mathrm{d\σ}}{\mathrm{dt}}=-\frac{\∂v}{\∂x}\cos \mathrm{\σ}+\frac{\∂v}{\∂z}\sin \mathrm{\σ}$

式中v为介质速度，σ为射线角度，(x，z)为射线坐标，给定射线初始位置(x₀，z₀)和角度σ₀，选定时间间隔Δt，利用常微分方程数值解法(如常用的四阶隆格库塔方法)求解上面的运动学方程组，获得中心射线路径以及沿中心射线路径的旅行时。

上述常微分方程数值解法是四阶隆格库塔方法。

5)根据炮点射线路径和旅行时，利用高斯射线束动力学方程，沿射线路径进行动力学射线追踪，计算炮点高斯射线束动力学参数；

步骤5)所述的高斯射线束动力学方程组如下：

$\frac{\mathrm{dq}\left(s\right)}{\mathrm{ds}}=\mathrm{vp}\left(s\right)$

$\frac{\mathrm{dp}\left(s\right)}{\mathrm{ds}}={-v}^2\frac{{\∂}^{2}v}{{\∂n}^2}q\left(s\right)---\left(3\right)$

式中v为介质速度，s为射线弧长，表示速度在垂直于射线方向上的二阶偏导数。

上述方程组利用常微分数值四阶隆格库塔方法对p、q的实部和虚部分别求解得到沿射线变化的动力学参数。

6)利用炮点高斯射线束动力学参数，计算射线束在各成像网格点上的复值振幅和旅行时A_s、T_s；

步骤6所述的高斯射线束式为：

u(ω)＝Aexp(iωT)                            (4)

式中u表示波场振幅分布，ω为角频率，i为虚数符号，为高斯射线束复值振幅，

为高斯射线束复值旅行时，s和n分别表示射线长度和到中心射线的距离，t(s)为s处中心射线旅行时。

7)对于炮集中的地震道，取出与之对应的处理后道数据b和

的样板表，设定检波点处射线角度范围

至

和射线角度间隔

8)对于给定的检波点处射线角度

利用运动学方程按一定时间间隔Δt进行射线追踪，利用步骤4)方法计算检波点射线路径和旅行时；

9)根据检波点射线路径和旅行时，利用高斯射线束动力学方程，沿射线路径进行动力学射线追踪，利用步骤5)方法计算检波点高斯射线束动力学参数；

10)利用检波点高斯射线束动力学参数，利用步骤6)方法计算射线束在各成像网格点上的复值振幅和旅行时A_d、T_d；

11)根据炮点高斯射线束和检波点射线束在各成像网格点上的复值旅行时和振幅A_s、T_s、A_d、T_d，按照高斯射线束偏移成像公式计算偏移结果并叠加；

步骤11)所述的叠前地震记录二维高斯射线束偏移的成像公式为：

式中I_s(x，z)表示地下成像点(x，z)处的偏移成像结果，γ为炮点处射线角度，为检波点处射线角度，b和

由前面步骤2)计算得到；

表示高斯射线束复值振幅和旅行时，下标r和i则表示其实部和虚部，

和

由炮点和检波点高斯射线束的振幅和旅行时共同决定：

式中T_s(γ)为炮点射线复值旅行时，A_s(γ)为炮点射线的复值振幅；T_d(γ)为检波点射线复值旅行时，A_d(γ)为检波点射线的复值振幅。

12)从

到

变化检波点射线角度

重复步骤8)-11)，计算并叠加偏移结果；

13)选取炮集中的下一个地震道，重复步骤7)-11)计算并叠加偏移结果；

14)从γ₁到γ_n变化炮点射线角度γ，重复步骤4)-12)，计算并叠加单炮记录的偏移结果；

15)重复步骤1)-14)，计算并叠加所有单炮的偏移结果，并显示为地层剖面图像。

本发明用高斯射线束替代Kirchhoff射线进行偏移，解决了射线的盲区、焦散区和多值走时问题，提高了成像效果。同时在改进Kirchhoff方法的同时保持了射线法灵活高效的优点。

附图说明

图1是Marmousi速度模型；

图2是本发明高斯射线束偏移与Kirchhoff方法偏移的对比效果图；

(a)Kirchhoff偏移结果，(b)本发明偏移结果，(c)单程波偏移结果。

图3是Kirchhoff、本发明方法对某一实际资料的偏移结果；

(a)Kirchhoff偏移结果，(b)本发明偏移结果，(c)单程波偏移结果。

具体实施方式

本发明的二维高斯射线束叠前偏移技术，其具体实施方式为：

1)激发地震并记录地震波，取得叠前单炮地震数据；

2)对炮集内各道地震数据f进行积分变换，得到处理后的道数据b及其Hilbert变换数据

其中b由叠前单道记录f的傅立叶频谱F(ω)通过下面的积分变换公式得到：

B(ω，T_i)＝exp(|ω|T_i)F(ω)

$b(T_r,T_i)=\frac{1}{2\mathrm{\π}}\∫\mathrm{d\ωexp}(-\mathrm{i\ω}{T}_r)B(\mathrm{\ω},T_i)---\left(1\right)$

式中T_r和T_i分别代表高斯射线束复值旅行时的实部和虚部，对T_r和T_i进行离散化，取ΔT_r＝地震记录采样间隔/4，

计算得到关于b和

的离散样板表。

步骤2所述的T_r的取值范围为0-原始记录最大采样时间，T_i的取值范围为

其中ω_l为地震数据有效频带的低端，取值为60弧/秒。

3)确定炮点处射线角度范围γ₁至γ_n，和射线角度间隔Δγ；

步骤3)是确定炮点处射线角度范围，并用角度间隔对它进行离散化，离散化后我们得到了不同角度的射线，如：确定射线角度范围[-80，80]，角度间隔20，那么每一条具体的射线角度就是{-80，-60，-40，-20，0，20，40，60，80}。

4)对于给定的炮点处射线角度γ，利用运动学方程按一定时间间隔Δt进行射线追踪，计算炮点射线路径和旅行时；

步骤4所述的运动学射线追踪通过求解如下的射线路径方程进行：

$\frac{\mathrm{dx}}{\mathrm{dt}}=v\sin \mathrm{\σ}$

$\frac{\mathrm{dz}}{\mathrm{dt}}=v\cos \mathrm{\σ}---\left(2\right)$

$\frac{\mathrm{d\σ}}{\mathrm{dt}}=-\frac{\∂v}{\∂x}\cos \mathrm{\σ}+\frac{\∂v}{\∂z}\sin \mathrm{\σ}$

式中v为介质速度，σ为射线角度，(x，z)为射线坐标，给定射线初始位置(x₀，z₀)和角度σ₀，选定时间间隔Δt，利用常微分方程数值解法(如常用的四阶隆格库塔方法)求解上面的运动学方程组，获得中心射线路径以及沿中心射线路径的旅行时。

上述常微分方程数值解法是四阶隆格库塔方法。

5)根据炮点射线路径和旅行时，利用高斯射线束动力学方程，沿射线路径进行动力学射线追踪，计算炮点高斯射线束动力学参数；

步骤5)所述的高斯射线束动力学方程组如下：

$\frac{\mathrm{dq}\left(s\right)}{\mathrm{ds}}=\mathrm{vp}\left(s\right)$

$\frac{\mathrm{dp}\left(s\right)}{\mathrm{ds}}={-v}^2\frac{{\∂}^{2}v}{{\∂n}^2}q\left(s\right)---\left(3\right)$

式中v为介质速度，s为射线弧长，

表示速度在垂直于射线方向上的二阶偏导数。

上述方程组利用常微分数值四阶隆格库塔方法对p、q的实部和虚部分别求解得到沿射线变化的动力学参数。

6)利用炮点高斯射线束动力学参数，计算射线束在各成像网格点上的复值振幅和旅行时A_s、T_s；

步骤6所述的高斯射线束式为：

u(ω)＝Aexp(iωT)                            (4)

式中u表示波场振幅分布，ω为角频率，i为虚数符号，

为高斯射线束复值振幅，

为高斯射线束复值旅行时，s和n分别表示射线长度和到中心射线的距离，t(s)为s处中心射线旅行时。

7)对于炮集中的地震道，取出与之对应的处理后道数据b和

的样板表，设定检波点处射线角度范围

至

和射线角度间隔

8)对于给定的检波点处射线角度

利用运动学方程按一定时间间隔Δt进行射线追踪，利用步骤4)方法计算检波点射线路径和旅行时；

9)根据检波点射线路径和旅行时，利用高斯射线束动力学方程，沿射线路径进行动力学射线追踪，利用步骤5)方法计算检波点高斯射线束动力学参数；

10)利用检波点高斯射线束动力学参数，利用步骤6)方法计算射线束在各成像网格点上的复值振幅和旅行时A_d、T_d；

11)根据炮点高斯射线束和检波点射线束在各成像网格点上的复值旅行时和振幅A_s、T_s、A_d、T_d，按照高斯射线束偏移成像公式计算偏移结果并叠加；

步骤11)所述的叠前地震记录二维高斯射线束偏移的成像公式为：

式中I_s(x，z)表示地下成像点(x，z)处的偏移成像结果，γ为炮点处射线角度，

为检波点处射线角度，b和

由前面步骤2)计算得到；表示高斯射线束复值振幅和旅行时，下标r和i则表示其实部和虚部，和

由炮点和检波点高斯射线束的振幅和旅行时共同决定：

式中T_s(γ)为炮点射线复值旅行时，A_s(γ)为炮点射线的复值振幅；T_d(γ)为检波点射线复值旅行时，A_d(γ)为检波点射线的复值振幅。

12)从

到

变化检波点射线角度

重复步骤8)-11)，计算并叠加偏移结果；

13)选取炮集中的下一个地震道，重复步骤7)-11)计算并叠加偏移结果；

14)从γ₁到γ_n变化炮点射线角度γ，重复步骤4)-12)，计算并叠加单炮记录的偏移结果；

15)重复步骤2)-13)，计算并叠加所有单炮的偏移结果。

图1是测试偏移算法常用Mamousi模型的，该模型结构复杂，在纵向和横向上都有较大的速度变化，地层和构造断层的倾角较大，模型上部主要由三大断裂组成，中间下方有一低速目的层。利用本发明对Mamousi模型全部240炮数据进行偏移。具体是：

(1)读取一个单炮地震数据；

(2)对炮集内各道地震数据f按公式(2)进行积分变换并离散化，建立b及其Hilbert变换的样板表，供成像计算时使用；

(3)设定炮点处射线角度范围[-75，75]，射线间隔3度；

(4)对于给定的炮点处射线角度γ，利用运动学方程从炮点位置开始，按原始记录时间间隔进行射线追踪，计算炮点射线路径和旅行时。；

(5)根据炮点射线路径和旅行时，利用高斯射线束动力学方程，沿射线路径进行动力学射线追踪，计算炮点高斯射线束动力学参数；

(6)利用炮点高斯射线束动力学参数，计算射线束在各成像网格点上的复值旅行时和振幅A_s、T_s；

(7)对于炮集中的某个地震道，取出与之对应的处理后道数据b和

样板表，设定检波点处射线角度范围-75至75，角度间隔2度；

(8)对于给定的检波点处射线角度利用运动学方程按原始记录时间间隔进行射线追踪，计算检波点射线路径和旅行时；

(9)根据检波点射线路径和旅行时，利用高斯射线束动力学方程，沿射线路径进行动力学射线追踪，计算检波点高斯射线束动力学参数；

(10)利用检波点高斯射线束动力学参数，计算射线束在各成像网格点上的复值旅行时和振幅A_d、T_d；

(11)根据炮点高斯射线束和检波点射线束在各成像网格点上的复值旅行时和振幅A_s、T_s、A_d、T_d，按公式(6)计算

查表得到对应的b(T_r，T_i)和

按照二维高斯射线束偏移成像公式(5)计算偏移结果并叠加；

(12)从-75度到75度变化检波点射线角度重复步骤(8)-(11)，

计算并叠加偏移结果；

(13)选取炮集中下一个地震道，重复步骤(7)-(11)计算并叠加偏移结果；

(14)从-75到75度变化炮点射线角度γ，重复步骤(4)-(12)，计算并叠加单炮记录的偏移结果；

(15)重复步骤(1)-(13)，计算并叠加所有240个单炮的偏移结果；

(16)将最终的偏移结果显示为地层剖面图像。

Kirchhoff、高斯射线束和单程方法对Marmousi模型进行偏移的结果如图2(a)、(b)、(c)所示，从偏移结果上看Kirchhoff方法在深部目的层上的成像效果较差，而高斯射线束方法能够清晰分辨深部目的层，且偏移噪声较小，成像结果接近于波动方程方法。而在偏移速度方面，高斯射线束方法耗时仅为Kirchhoff方法的两倍左右，远高于单程波偏移方法。

图3是利用Kirchhoff、高斯射线束和单程波偏移方法对某一实际地震资料(共360炮记录，本发明实例)偏移的结果，高斯射线束偏移流程及选用参数与实例(1)一致，高斯射线束偏移效果接近于单程波偏移。

Kirchhoff方法将波场分解到完全理想的高频射线上，而本方法将波场分解到与频率相关的具有一定宽度的射线束上来实现波场的正演和偏移。在高斯射线束偏移中，地下成像点处的波场不是简单地由某条特定的射线所确定，而是将多条相关的高斯射线束叠加而形成的。高斯射线束的性质决定了利用高斯射线束进行正演和偏移，可以在很大程度上克服普通射线法存在的盲区、焦散区和多值走时等问题。

Claims (5)

1.一种二维高斯射线束叠前偏移方法，其特征是具体步骤包括：

1)激发地震并记录地震波，取得叠前单炮地震数据；

2)对炮集内各道地震数据f进行积分变换，得到处理后的道数据b及其Hilbert变换数据

其中b由叠前单道记录f的傅立叶频谱F(ω)通过下面的积分变换公式得到：

B(ω，T_i)＝exp(|ω|T_i)F(ω)

$b(T_r,T_i)=\frac{1}{2\mathrm{\π}}\∫\mathrm{d\ωexp}(-\mathrm{i\ω}{T}_r)B(\mathrm{\ω},T_i)---\left(1\right)$

式中T_r和T_i分别代表高斯射线束复值旅行时的实部和虚部，对T_r和T_i进行离散化，取ΔT_r＝地震记录采样间隔/4，

计算得到关于b和

的离散样板表；

所述的T_r的取值范围为0-原始记录最大采样时间，T_i的取值范围为

其中ω_l为地震数据有效频带的低端，取值为60弧/秒；

3)确定炮点处射线角度范围γ₁至γ_n，和射线角度间隔Δγ；

4)对于给定的炮点处射线角度γ，利用运动学方程按一定时间间隔Δt进行射线追踪，计算炮点射线路径和旅行时；

5)根据炮点射线路径和旅行时，利用高斯射线束动力学方程组，沿射线路径进行动力学射线追踪，计算炮点高斯射线束动力学参数；

6)利用炮点高斯射线束动力学参数，计算射线束在各成像网格点上的复值振幅和旅行时A_s、T_s；

7)对于炮集中的地震道，取出与之对应的处理后道数据b和

的样板表，设定检波点处射线角度范围

至和射线角度间隔

8)对于给定的检波点处射线角度

利用运动学方程按一定时间间隔Δt进行射线追踪，利用步骤4)方法计算检波点射线路径和旅行时；

9)根据检波点射线路径和旅行时，利用高斯射线束动力学方程，沿射线路径进行动力学射线追踪，利用步骤5)方法计算检波点高斯射线束动力学参数；

10)利用检波点高斯射线束动力学参数，利用步骤6)方法计算射线束在各成像网格点上的复值振幅和旅行时A_d、T_d；

11)根据炮点高斯射线束和检波点射线束在各成像网格点上的复值旅行时和振幅A_s、T_s、A_d、T_d，按照高斯射线束偏移成像公式计算偏移结果并叠加；

12)从到

变化检波点射线角度重复步骤8)-11)，计算并叠加偏移结果；

13)选取炮集中的下一个地震道，重复步骤7)-11)计算并叠加偏移结果；

14)从γ₁到γ_n变化炮点射线角度γ，重复步骤4)-12)，计算并叠加单炮记录的偏移结果；

15)重复步骤2)-14)，计算并叠加所有单炮的偏移结果，并显示为地层剖面图像。

2.根据权利要求1的二维高斯射线束叠前偏移方法，其特征是步骤4)所述的射线追踪通过求解如下的射线路径方程进行：

$\frac{\mathrm{dx}}{\mathrm{dt}}=v\sin \mathrm{\σ}$

$\frac{\mathrm{dz}}{\mathrm{dt}}=v\cos \mathrm{\σ}---\left(2\right)$

$\frac{\mathrm{d\σ}}{\mathrm{dt}}=-\frac{\∂v}{\∂x}\cos \mathrm{\σ}+\frac{\∂v}{\∂z}\sin \mathrm{\σ}$

式中v为介质速度，σ为射线角度，(x，z)为射线坐标，给定射线初始位置(x₀，z₀)和角度σ₀，选定时间间隔Δt，利用常微分方程数值解法求解上面的运动学方程组，获得中心射线路径以及沿中心射线路径的旅行时。

3.根据权利要求2的二维高斯射线束叠前偏移方法，其特征是上述常微分方程数值解法是四阶隆格库塔方法。

4.根据权利要求1的二维高斯射线束叠前偏移方法，其特征是步骤5)所述的高斯射线束动力学方程组如下：

$\frac{\mathrm{dq}\left(s\right)}{\mathrm{ds}}=\mathrm{vp}\left(s\right)$

$\frac{\mathrm{dp}\left(s\right)}{\mathrm{ds}}={-v}^2\frac{{\∂}^{2}v}{{\∂n}^2}q\left(s\right)---\left(3\right)$

式中v为介质速度，s为射线弧长，

表示速度在垂直于射线方向上的二阶偏导数；

上述方程组利用常微分数值四阶隆格库塔方法对p、q的实部和虚部分别求解得到沿射线变化的动力学参数；

所述的高斯射线束式为：

u(ω)＝Aexp(iωT)                    (4)

式中u表示波场振幅分布，ω为角频率，i为虚数符号，

为高斯射线束复值振幅，

为高斯射线束复值旅行时，s和n分别表示射线弧长和到中心射线的距离，t(s)为s处中心射线旅行时。

5.根据权利要求1的二维高斯射线束叠前偏移方法，其特征是步骤11)所述的叠前地震记录二维高斯射线束偏移的成像公式为：

式中I_s(x，z)表示地下成像点(x，z)处的偏移成像结果，γ为炮点处射线角度，为检波点处射线角度，b和

由前面步骤2)计算得到；

表示高斯射线束复值振幅和旅行时，下标r和i则表示其实部和虚部，

和

由炮点和检波点高斯射线束的振幅和旅行时共同决定：

式中T_s(γ)为炮点射线复值旅行时，A_s(γ)为炮点射线的复值振幅；T_d(γ)为检波点射线复值旅行时，A_d(γ)为检波点射线的复值振幅。

Images