CN105629299A - 角度域叠前深度偏移的走时、角度表获取方法及成像方法 - Google Patents

Abstract

本发明涉及一种各向异性介质角度域基尔霍夫叠前深度偏移的走时表和角度表获取方法及其成像方法。所述获取方法包括：根据各向异性介质声波方程推导出各向异性介质中的程函方程，通过程函方程推导得到各向异性介质声学近似的射线方程组并求解，从而获得各向异性介质声学近似意义下的射线路径、走时及传播方向信息，进行角度域射线追踪，并建立走时表与出射角度表。本发明的优点包括：能够基于声学规律近似推导出各向异性介质的射线方程；能够为基尔霍夫叠前深度偏移提供可靠的走时表和角度表，从而有利于准确实现角度域成像。

Description

角度域叠前深度偏移的走时、角度表获取方法及成像方法

技术领域

本发明涉及一种各向异性(VTI)介质角度域基尔霍夫(Kirchhoff)叠前深度偏移的走时表和角度表获取方法以及一种各向异性介质角度域基尔霍夫叠前深度偏移的成像方法。

背景技术

目前，叠前深度偏移是强横向非均匀介质复杂构造成像与速度模型建立依赖的关键技术。其算法实现要么基于射线理论，如基尔霍夫(Kirchhoff)偏移和高斯束偏移，要么基于波动理论，如单程波方程深度延拓偏移和双程波方程逆时延拓偏移。近十多年来，各向异性介质深度偏移方法也得到了极大的发展，先后出现了横向各项同性(TI)介质基尔霍夫偏移、高斯束偏移、单程波方程偏移与逆时偏移等深度域成像方法。尽管波动方程偏移存在精度上的优势，射线理论基础上的偏移方法因其在灵活性、面向局部目标的成像能力以及计算成本等优势，在复杂构造成像尤其是速度模型建立过程中得到广泛应用。目前主要地震数据处理软件中的深度域偏移速度模型构建都仍以基尔霍夫偏移作为引擎。

在复杂介质条件下，即使偏移速度是合理的，传统的偏移距域和炮域共成像点道集都可能存在假象干扰。为此，近十几年来人们一直在致力于研究射线理论或波动理论基础上的角度域成像方法。基于射线理论和广义拉冬变换(GRT)，deHoop提出了共散射角偏移/反演理论。

就目前而言，如何高保真地提取地下方位角和反射角信息，虽然经过了10余年的研究，发表了不少的文献，但是做到真正地下反射点的方位角信息和反射角信息，也只是最近一些年才成为现实。针对叠前时间偏移的基尔霍夫的高保真方位角和反射角道集信息的获取，在2008年和2011年才相继完成(Chengetal.，2008，2011)。最近Koren和Ravve(2011)实现和应用了方位保真局部角度域基尔霍夫深度成像方法，他们利用了广义拉冬变换抽取全方位、高分辨率角度有关的反射率信息，并产生了角度域道集。该偏移技术已集成到Paradigm公司的商业软件中，取得了良好的应用效果。

局部角度域基尔霍夫叠前深度偏移算法的核心在于稳健、快速地计算地震射线的走时与方向信息。在传统基尔霍夫叠前深度偏移过程中，程函方程有限差分解法与波前重建算法被广泛用于走时表的计算。然而，对局部角度域成像与反射走时层析基础上的偏移速度分析而言，射线追踪算法显得更有吸引力，因为它除了计算走时，还可以显式地得到射线路径及其方向信息。不过，传统的各向异性介质射线追踪方程是以刚度系数而不是Thomsen参数表示的，不方便数值计算，效率也较低。

发明内容

针对现有技术存在的诸多不足中的至少一项，本发明提出一种适用于各向异性介质角度域叠前深度偏移的走时表与角度表计算方法，以及基于该走时表与角度表计算方法的各向异性介质角度域叠前深度偏移的成像方法。

本发明一方面提供了一种各向异性介质角度域叠前深度偏移的走时表和角度表获取方法。所述方法包括以下步骤：

A、根据各向异性介质声波方程推导出各向异性介质中的程函方程：

${V_{nmo}}^2(1+2\mathrm{\η})({\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂y}\right)}^2)+V_{p0}^2{\left(\frac{\∂\mathrm{\τ}}{\∂z}\right)}^2\×(1-2{V_{nmo}}^2\mathrm{\η}({\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂y}\right)}^2\left)\right)=1---\left(1\right)$

其中，V_nmo为NMO速度，V_p0为qP波垂直速度，η为反椭圆系数，τ表示沿射线的走时，x，y和z表示各向异性介质空间上的三个彼此垂直的位移分量(即，沿着笛卡尔坐标系方向的三个位移分量)；

B、通过程函方程推导得到各向异性介质声学近似的射线方程组：

$\frac{{\mathrm{dx}}_i}{d\mathrm{\τ}}=\frac{1}{2}\frac{\∂F}{\∂p_i}$ 和 $\frac{{\mathrm{dp}}_i}{d\mathrm{\τ}}=-\frac{1}{2}\frac{\∂F}{\∂x_i},$

其中，设上面的式(1)即为F，p_i为慢度矢量，i对应x，y和z分量；

C、求解步骤B中的射线方程组，获得各向异性介质声学近似意义下的射线路径、走时及传播方向信息；

D、进行角度域射线追踪，并建立走时表与出射角度表。

在本发明的各向异性介质角度域叠前深度偏移的走时表和角度表获取方法的一个示例性实施例中，优选地，所述步骤D中进行角度与射线追踪通过从地下多个成像点中的每个成像点以起飞角β_s或β_r与方位角α_s或α_r等间隔向上发射一族射线到达地表各观测点，将所述一族射线中不同方向起飞射线的走时与角度信息保存在数值表中，从而形成走时表与出射角度表。

在本发明的各向异性介质角度域叠前深度偏移的走时表和角度表获取方法的一个示例性实施例中，优选地，所述步骤C中的传播方向信息通过入射角γ、散射方位角φ、照明矢量的倾角与方位角来表征，并且通过下列方程得到：

p_s＝(p_x,p_y,p_z)＝(sinβ_scosα_s,sinβ_ssinα_s,cosβ)，

p_r＝(p_x,p_y,p_z)＝(sinβ_rcosα_r,sinβ_rsinα_r,cosβ_r)，

$cos\mathrm{\θ}=cos\left(2\mathrm{\γ}\right)=\frac{p_s\·p_r}{\left|p_s\right|\left|p_r\right|},$

$cos\mathrm{\φ}=\frac{(p_m\×y)\·(p_r\×p_s)}{|p_m\×y||p_r\×p_s|},$

其中，入射慢度矢量p_s和散射慢度矢量p_r共同描述了散射点m处波的传播方向特征，入射慢度矢量与散射慢度矢量之和p_m，x，y和z表示各向异性介质空间上的三个彼此垂直的位移分量(即，沿着笛卡尔坐标系方向的三个位移分量)，为照明矢量的垂向分量。

本发明的另一方面提供了一种各向异性介质角度域叠前深度偏移的成像方法。所述方法在局部角度域成像时，根据炮点、成像点、接收点关系，从上述获取方法中所获取的走时表和角度表中读取走时和角度数据，以完成成像。

与现有技术相比，本发明的有益效果包括：能够基于声学规律近似推导出各向异性介质的射线方程；能够为基尔霍夫叠前深度偏移提供可靠的走时表和角度表，从而有利于准确实现角度域成像。

附图说明

图1示出了本发明的示例性实施例中的某一成像点处地震波局部角度特征示意图。

具体实施方式

在下文中，将结合示例性实施例来详细说明本发明的各向异性介质角度域基尔霍夫叠前深度偏移的走时表和角度表获取方法以及各向异性介质角度域基尔霍夫叠前深度偏移的成像方法。

在本发明的一个示例性实施例中，各向异性介质角度域叠前深度偏移的走时表和角度表获取方法可通过以下步骤来实现：

1、根据各向异性介质声波方程推导出各向异性介质中的程函方程

所谓声学近似，就是假设沿对称轴方向qSV波(这里，qSV波全称为拟横波，大概意思就是类似于横波的，但不是真正意义上各向同性介质中的横波)的传播速度为零，即V_S0＝0，这样就可将原始的各向异性介质弹性波动方程及其频散关系简化。假设地下介质为声学介质，由各向异性介质弹性波动方程及其频散关系可推导出近似的qP波(这里，qP波全称为拟纵波)标量波动方程，进而得到相应的程函方程和射线方程。根据各向异性介质qP波的频散关系，声学近似qP波波动方程满足：

$\frac{{\∂}^{4}F}{\∂t^2}-(1+2\mathrm{\η}){V_{nmo}}^2(\frac{{\∂}^{4}F}{\∂x^2\∂t^2}+\frac{{\∂}^{4}F}{\∂y^2\∂t^2})=V_{p0}^2\frac{{\∂}^{2}F}{\∂z^2\∂t^2}-2{{\mathrm{\ηV}}_{nmo}}^2{V}_{p0}^2(\frac{{\∂}^{4}F}{\∂x^2\∂z^2}+\frac{{\∂}^{4}F}{\∂y^2\∂z^2})---\left(1\right)$

其中，V_p0为qP波垂直速度，V_nmo为NMO速度(NMO速度全称是动矫正速度)，η为反椭圆系数，且与Thomsen参数ε与δ存在如下关系：

$V_{nmo}=V_{p0}\sqrt{1+2\mathrm{\δ}}---\left(2a\right)$

$\mathrm{\η}=\frac{\mathrm{\ϵ}-\mathrm{\δ}}{1+2\mathrm{\δ}}---\left(2b\right)$

将平面波解带入方程(1)可推导出各向异性介质的程函方程：

${V_{nmo}}^2(1+2\mathrm{\η})({\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂y}\right)}^2)+V_{p0}^2{\left(\frac{\∂\mathrm{\τ}}{\∂z}\right)}^2\×(1-2{V_{nmo}}^2\mathrm{\η}({\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂y}\right)}^2\left)\right)=1---\left(30\right)$

2、通过程函方程推导得到各向异性介质声学近似的射线方程

通过特征值方法可进一步推导出描述射线路径的常微分方程组。为此，将式(3)改写为如下形式：

F(x,y,z,p_x,p_y,p_z)＝0(4)

其中，和为慢度矢量的三个分量。通过特征值法可得到射线方程组：

$\frac{{\mathrm{dx}}_i}{d\mathrm{\τ}}=\frac{1}{2}\frac{\∂F}{\∂p_i}---\left(5a\right)$

$\frac{{\mathrm{dp}}_i}{d\mathrm{\τ}}=-\frac{1}{2}\frac{\∂F}{\∂x_i}---\left(5b\right)$

其中τ代表沿着射线的走时，i对应x，y和z分量。该方程组描述了各向异性介质声学近似意义下的射线路径、走时及传播方向信息。

3、求解射线方程

例如，通过4阶龙格库塔法求解射线方程组(5a)和(5b)。这里，若用低阶的龙格库塔法求解射线方程组，则精度略有降低，因此，优选用4阶龙格库塔法求解射线方程组(5a)和(5b)。

4、角度域射线追踪

在图1中，点S表示某炮点的位置，点R表示某检波点的位置，点S和点R均位于地表面；曲线Sisochrone表示某炮点等时线，曲线Risochrone表示某检波点等时线。

如图1所示，三维情况下，入射慢度矢量p_s和散射慢度矢量p_r共同描述了散射点m处波的传播方向特征。入射慢度矢量与散射慢度矢量之和p_m称为照明矢量。根据地震勘探的需要，可用两类、四个角度共同定义局部传播方向。第一类是描述入射与散射(包括绕射和反射)方向特征的两个角度，即入射角γ(散射张角θ的一半)和散射方位角(即局部入射慢度与散射慢度所在平面的方位角)φ。第二类是描述局部照明方向的两个角度，即照明矢量的倾角与方位角基于射线理论，这四个角度参数可由走时的空间梯度计算得到。

设入射射线的起飞角为β_s，方位角为α_s，散射射线的起飞角为β_r，方位角为α_r，可得入射射线与散射射线的单位慢度矢量：

p_s＝(p_x,p_y,p_z)＝(sinβ_scosα_s,sinβ_ssinα_s,cosβ)(6a)

p_r＝(p_x,p_y,p_z)＝(sinβ_rcosα_r,sinβ_rsinα_r,cosβ_r)(6b)

式(6a)中β即为β_s。根据矢量运算法则，前文所述四个局部角度参数分别满足：

$cos\mathrm{\θ}=cos\left(2\mathrm{\γ}\right)=\frac{p_s\·p_r}{\left|p_s\right|\left|p_r\right|}---\left(7a\right)$

$cos\mathrm{\φ}=\frac{(p_m\×y)\·(p_r\×p_s)}{|p_m\×y||p_r\×p_s|}---\left(7b\right)$

式中x、y与z分别代表沿坐标轴的单位矢量，其中y指向正北方向并作为定义方位角的参考方向，为照明矢量的垂向分量。可见，只要根据起飞角及其方位角计算得到入射与散射慢度矢量，就可根据上述方程求取四个局部角度参数γ、φ、和

在局部角度域进行射线追踪时，从地下成像点m(如图1)以起飞角β_s(或β_r)与方位角α_s(或α_r)等间隔向上发射一族射线到达地表各观测点，将这些不同方向起飞射线的走时与角度信息保存在数值表中，形成走时表和出射角度表。

在局部角度域成像时，便可根据炮点-成像点-接收点关系，在数值表中读取计算好的样本通过插值获得实际射线路径的走时，起飞角及其方位角，进而按公式(6a)至(7d)从射线路径的局部角度参数(β_s,α_s；β_r,α_r)转换成局部角度域成像需要的角度参数(γ；φ,)。

总体来讲，本发明能够给出适用于各向异性介质角度域叠前深度偏移的走时计算方法。本发明能够基于声学近似推导出各向异性介质的射线方程；而且与传统射线追踪计算走时方法不同，该方法从地下每个成像点处激发一簇射线到达地表，同时，记录该成像点m达到地表检波点与炮点的走时及出射角度并建表。对实际地震资料的偏移结果表明，该方法能够为基尔霍夫叠前深度偏移提供一个可靠的走时表和角度表，从而实现角度域成像。

综上所述，本发明能够提供一种适用于各向异性介质基尔霍夫叠前深度偏移的走时计算方法，其除了像传统基尔霍夫叠前深度偏移那样输出成像剖面和炮检距域的共成像点道集，还遵循地震波在成像点处的局部方向特征、基于扩展的脉冲响应叠加原理获得入射角度域和照明角度域的成像结果。本发明的方法既适用于复杂构造成像，又可为各向异性介质深度域偏移速度分析与模型建立提供高效的偏移引擎。

尽管上面已经结合附图和示例性实施例描述了本发明，但是本领域普通技术人员应该清楚，在不脱离权利要求的精神和范围的情况下，可以对上述实施例进行各种修改。

Claims (4)

1.一种各向异性介质角度域叠前深度偏移的走时表和角度表获取方法，其特征在于，所述方法包括步骤：

A、根据各向异性介质声波方程推导出各向异性介质中的程函方程：

${V_{nmo}}^2(1+2\mathrm{\η})({\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂y}\right)}^2)+V_{p0}^2{\left(\frac{\∂\mathrm{\τ}}{\∂z}\right)}^2\×(1-2{V_{nmo}}^2\mathrm{\η}({\left(\frac{\∂\mathrm{\τ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\τ}}{\∂y}\right)}^2\left)\right)=1,$

其中，V_nmo为NMO速度，V_p0为qP波垂直速度，η为反椭圆系数，τ表示沿射线的走时，x，y和z表示各向异性介质空间上的三个彼此垂直的位移分量；

B、通过程函方程推导得到各向异性介质声学近似的射线方程组：

$\frac{{\mathrm{dx}}_i}{d\mathrm{\τ}}=\frac{1}{2}\frac{\∂F}{\∂p_i}$ 和 $\frac{{\mathrm{dp}}_i}{d\mathrm{\τ}}=-\frac{1}{2}\frac{\∂F}{\∂x_i},$

其中，F为……，p_i为……，i对应x，y和z分量；

C、求解步骤B中的射线方程组，获得各向异性介质声学近似意义下的射线路径、走时及传播方向信息；

D、进行角度域射线追踪，并建立走时表与出射角度表。

2.根据权利要求1所述的各向异性介质角度域叠前深度偏移的走时表和角度表获取方法，其特征在于，所述步骤D中进行角度与射线追踪通过从地下多个成像点中的每个成像点以起飞角β_s或β_r与方位角α_s或α_r等间隔向上发射一族射线到达地表各观测点，将所述一族射线中不同方向起飞射线的走时与角度信息保存在数值表中，形成走时表与出射角度表。

3.根据权利要求1所述的各向异性介质角度域叠前深度偏移的走时表和角度表获取方法，其特征在于，所述步骤C中的传播方向信息通过入射角γ、散射方位角φ、照明矢量的倾角与方位角来表征，并且通过下列方程得到：

p_s＝(p_x,p_y,p_z)＝(sinβ_scosα_s,sinβ_ssinα_s,cosβ)，

p_r＝(p_x,p_y,p_z)＝(sinβ_rcosα_r,sinβ_rsinα_r,cosβ_r)，

$cos\mathrm{\θ}=cos\left(2\mathrm{\γ}\right)=\frac{p_s\·p_r}{\left|p_s\right|\left|p_r\right|},$

$cos\mathrm{\φ}=\frac{(p_m\×y)\·(p_r\×p_s)}{|p_m\×y||p_r\×p_s|},$

其中，入射慢度矢量p_s和散射慢度矢量p_r共同描述了散射点m处波的传播方向特征，入射慢度矢量与散射慢度矢量之和p_m，x、y与z分别代表沿坐标轴的单位矢量，为照明矢量的垂向分量。

4.一种各向异性介质角度域叠前深度偏移的成像方法，其特征在于，所述方法在局部角度域成像时，根据炮点、成像点、接收点关系，从权利要求1至3中任意一项的方法中所获取的走时表和角度表中读取走时和角度数据，进行成像。