CN107229067B - 层状模型各向异性射线追踪方法和系统 - Google Patents

Abstract

本发明提供了层状模型各向异性射线追踪方法和系统，包括：获取模型原始信息，对模型原始信息进行参数化得到参数值；给定震源位置和射线第一相速度方向，根据参数值计算相速度矢量和群速度矢量，从而得到射线层间传播信息；判断射线是否到达模型边界；如果未到达模型边界，则计算射线经过当前界面的第二相速度方向；判断射线是否能够实际继续传播，如果不能，则判断射线是否符合调整条件；如果符合，则对射线进行调整得到第三相速度方向；重复上述计算射线在层间传播的步骤，直至射线到达模型边界或无法继续在模型中传播，得到射线的完整路径。本发明能够对复杂三维各向异性介质进行建模，为动力学射线追踪及振幅计算提供了帮助。

Description

层状模型各向异性射线追踪方法和系统

技术领域

本发明涉及射线追踪技术领域，尤其是涉及层状模型各向异性射线追踪方法。

背景技术

使用网格或者元胞进行介质建模，模型的精确程度与单位网格(元胞)的剖分程度相关，即网格(元胞)剖分越密，获取的模型越精确，模型对存储容量的需求越高，对于各向异性介质的网格建模技术，不仅需要存储速度信息，与弹性介质相关的弹性参数信息(密度、P/S波速比、刚度模量)也是必须的。其中，刚度模量参数数量与各向异性介质复杂程度相关。此外，网格模型通常需要进行平滑，其无法准确表达真实的地球模型特征。网格模型的界面无法通过模型进行区分，在部分应用中需要人为给定。针对TTI介质(倾斜横向各向同性介质)，各向异性与地层走向相关的倾角场相关，网格模型没有与倾角场自然相关的参数，需要增加倾角场模型。综上，对于复杂二维、三维介质来说，各向异性网格模型建模具有一定的局限性。

综上所述，现有的层状、块状模型建模技术一般应用于各向同性介质中，各向异性参数使用层状模型进行建模技术存在空缺。

发明内容

有鉴于此，本发明的目的在于提供层状模型各向异性射线追踪方法，能够对复杂三维各向异性介质进行建模，为动力学射线追踪及振幅计算提供了稳定的导数信息。

第一方面，本发明实施例提供了层状模型各向异性射线追踪方法，其特征在于，包括：

获取模型原始信息，并对所述模型原始信息进行参数化得到参数值；

给定震源位置和射线的第一相速度方向，根据所述参数值计算所述射线的相速度矢量和群速度矢量，以得到射线层间传播信息；

根据所述射线层间传播信息判断所述射线是否到达模型边界；

如果未到达所述模型边界，则计算所述射线经过当前界面的第二相速度方向；

根据所述第二相速度方向与界面的关系，判断所述射线是否能够实际继续传播，如果不能，则判断所述射线是否符合调整条件；

如果符合，则对所述射线进行调整得到第三相速度方向；

重复上述计算射线在层间传播的步骤，直至所述射线到达所述模型边界或无法继续在模型中传播，得到所述射线的完整路径。

结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，所述模型原始信息包括界面数据和弹性参数数据，所述对所述模型原始信息进行参数化得到参数值包括：

对所述界面数据和所述弹性参数数据反算得到非均匀分布的具有三次B样条性质的控制节点；

使用所述控制节点进行拟合得到参数化模型信息，其中，所述参数化模型信息包括基函数、节点信息和背景场；

使用所述基函数、所述节点信息和所述背景场通过三次B样条函数进行插值得到与所述参数化模型对应的所述参数值。

结合第一方面，本发明实施例提供了第一方面的第二种可能的实施方式，其中，所述相速度矢量包括相速度大小和所述第一相速度方向，所述根据所述参数值确定所述射线的相速度矢量和群速度矢量包括：

根据下式计算所述相速度矢量和群速度矢量：

(Γ_ik-δ_ik)U_k＝0,Γ_ik＝a_ijklp_jp_l,a_ijkl＝C_ijkl/ρ

其中，Γ为Christoffel矩阵，δ_ik为所述Kronecker delta函数，U_k为复振幅向量在k方向的分量，p为慢度矢量，a_ijkl为正则化的弹性模量，ρ为密度，C为弹性模量，下标i,j,k,l为爱因斯坦求和的指标；

具体的计算过程为根据下式计算所述慢度矢量：

由此，

根据下式计算所述相速度矢量：

v_n＝n′v_n

其中，n′为所述相速度方向单位矢量，v_n为所述相速度大小，v_n为所述相速度矢量；

另外，

根据下式计算所述群速度矢量：

其中，v_g为所述群速度矢量，g^(m)为所述Christoffel方程的特征值对应的特征向量，p为所述慢度矢量，a_ijkl为所述正则化的弹性模量，下标i,j,k,l为爱因斯坦求和的指标。

结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，所述计算所述射线经过当前界面的第二相速度方向包括：

将所述相速度矢量进行全局坐标系至局部坐标系的转换；

对转化后的所述相速度矢量通过广义斯内尔定理得到所述射线经过所述当前界面的所述第二相速度方向。

结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方式，其中，所述调整条件包括所述射线的参数条件、界面的倾角条件和入射角条件，所述对所述射线进行调整得到第三相速度方向包括：

根据下式计算所述第三相速度方向：

其中，f(ξ_i)为高斯函数，ξ为可调整的最大正弦值，ζ为调整后的最小正弦值，ξ_i为欲调整的入射角正弦值，ζ_i为调整后的所述第三相速度对应的所述入射角正弦值，i为当前界面的层数。

第二方面，本发明实施例提供了层状模型各向异性射线追踪系统，包括:

获取单元，用于获取模型原始信息，并对所述模型原始信息进行参数化得到参数值；

第一计算单元，用于给定震源位置和射线的第一相速度方向，根据所述参数值计算所述射线的相速度矢量和群速度矢量，以得到射线层间传播信息；

第一判断单元，用于根据所述射线层间传播信息判断所述射线是否到达模型边界；

第二计算单元，用于在所述射线未到达所述边界的情况下，计算所述射线经过当前界面的第二相速度方向；

第二判断单元，用于根据所述第二相速度方向与界面的关系，判断所述射线是否能够实际继续传播；

第三判断单元，用于判断所述射线是否符合调整条件；

调整单元，用于对所述射线进行调整得到第三相速度方向；

循环单元，用于重复上述计算射线在层间传播的步骤，直至所述射线到达所述模型边界或无法继续在模型中传播，得到所述射线的完整路径。

结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，所述模型原始信息包括界面数据和弹性参数数据，所述获取单元包括：

对所述界面数据和所述弹性参数数据反算得到非均匀分布的具有三次B样条性质的控制节点；

使用所述控制节点进行拟合得到参数化模型信息，其中，所述参数化模型信息包括基函数、节点信息和背景场；

使用所述基函数、所述节点信息和所述背景场通过三次B样条函数进行插值得到与所述参数化模型对应的所述参数值。

结合第二方面，本发明实施例提供了第二方面的第二种可能的实施方式，其中，所述相速度矢量包括相速度大小和所述第一相速度方向，所述第一计算单元包括：

根据下式计算所述相速度矢量和群速度矢量：

(Γ_ik-δ_ik)U_k＝0,Γ_ik＝a_ijklp_jp_l,a_ijkl＝C_ijkl/ρ

其中，Γ为Christoffel矩阵，δ_ik为所述Kronecker delta函数，U_k为复振幅向量在k方向的分量，p为慢度矢量，a_ijkl为正则化的弹性模量，ρ为密度，C为弹性模量，下标i,j,k,l为爱因斯坦求和的指标；

具体的计算过程为根据下式计算所述慢度矢量：

由此，

根据下式计算所述相速度矢量：

v_n＝n′v_n

其中，n′为所述相速度方向单位矢量，v_n为所述相速度大小，v_n为所述相速度矢量；

另外，

根据下式计算所述群速度矢量：

其中，v_g为所述群速度矢量，g^(m)为所述Christoffel方程的特征值对应的特征向量，p为所述慢度矢量，a_ijkl为所述正则化的弹性模量，下标i,j,k,l为爱因斯坦求和的指标。

结合第二方面，本发明实施例提供了第二方面的第三种可能的实施方式，其中，所述第二计算单元包括：

将所述相速度矢量进行全局坐标系至局部坐标系的转换；

对转化后的所述相速度矢量通过广义斯内尔定理得到所述射线经过所述当前界面的所述第二相速度方向。

结合第二方面，本发明实施例提供了第二方面的第四种可能的实施方式，其中，所述调整条件包括所述射线的参数条件、界面的倾角条件和入射角条件，所述调整单元包括：

根据下式计算所述第三相速度方向：

其中，f(ξ_i)为高斯函数，ξ为可调整的最大正弦值，ζ为调整后的最小正弦值，ξ_i为欲调整的入射角正弦值，ζ_i为调整后的所述第三相速度对应的所述入射角正弦值，i为当前界面的层数。

本发明提供了层状模型各向异性射线追踪方法和系统，包括：获取模型原始信息，对模型原始信息进行参数化得到参数值；给定震源位置和射线第一相速度方向，根据参数值计算相速度矢量和群速度矢量，从而得到射线层间传播信息；判断射线是否到达模型边界；如果未到达模型边界，则计算射线经过当前界面的第二相速度方向；判断射线是否能够实际继续传播，如果不能，则判断射线是否符合调整条件；如果符合，则对射线进行调整得到第三相速度方向；重复上述计算射线在层间传播的步骤，直至射线到达模型边界或无法继续在模型中传播，得到射线的完整路径。本发明能够对复杂三维各向异性介质进行建模，为动力学射线追踪及振幅计算提供了帮助。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的层状模型各向异性射线追踪方法流程图；

图2为本发明实施例提供的步骤S101方法流程图；

图3为本发明实施例提供的步骤S104方法流程图；

图4为本发明实施例提供的层状模型各向异性射线追踪系统结构示意图；

图5(a)为本发明实施例提供的三维层状各向同性模型示意图；

图5(b)为本发明实施例提供的使用各向同性模型计算的射线在3000m深度的走时分布；

图5(c)为本发明实施例提供的使用各向同性模型计算的射线在5000m深度的走时分布；

图5(d)为本发明实施例提供的三维层状各向异性模型示意图；

图5(e)为本发明实施例提供的使用各向异性模型计算的射线在3000m深度的走时分布；

图5(f)为本发明实施例提供的使用各向异性模型计算的射线在3000m深度的走时分布；

图6(a)为本发明实施例提供的射线追踪的结果示意图；

图6(b)为本发明实施例提供的PSPI脉冲响应结果图；

图6(c)为本发明实施例提供的不含经过调整的射线段的脉冲响应结果图；

图6(d)为本发明实施例提供的含有经过调整的射线段的脉冲响应结果图；

图7(a)为本发明实施例提供的使用复杂层状模型建模建立的Sigsbee层状模型图；

图7(b)为本发明实施例提供的在区域1比较层状模型偏移结果、真实速度模型与网格模型偏移结果图；

图7(c)为本发明实施例提供的在区域2比较层状模型偏移结果、真实速度模型与网格模型偏移结果图；

图8(a)为本发明实施例提供的P波速度示意图；

图8(b)为本发明实施例提供的密度示意图；

图8(c)为本发明实施例提供的Thomsen’s各向异性参数ε；

图8(d)为本发明实施例提供的另一Thomsen’s各向异性参数δ；

图9(a)为本发明实施例提供的基于网格建模的各向异性Hess模型高斯束偏移结果示意图；

图9(b)为本发明实施例提供的另一基于网格建模的各向同性Hess模型高斯束偏移结果示意图；

图9(c)为本发明实施例提供的基于层状模型的各向异性Hess模型高斯束偏移结果示意图；

图9(d)为本发明实施例提供的另一基于层状模型的各向同性Hess模型高斯束偏移结果示意图；

图10为本发明实施例提供的另一层状模型各向异性射线追踪方法流程图。

图标：

10-获取单元；20-第一计算单元；30-第一判断单元；40-第二计算单元；50-第二判断单元；60-第三判断单元；70-调整单元；80-循环单元。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

目前，现有的层状、块状模型建模技术一般应用于各向同性介质中，各向异性参数使用层状模型进行建模技术存在空缺，基于此，本发明实施例提供的层状模型各向异性射线追踪方法和系统，能够对复杂三维各向异性介质进行建模，为动力学射线追踪及振幅计算提供了帮助。

为便于对本实施例进行理解，首先对本发明实施例所公开的层状模型各向异性射线追踪方法进行详细介绍。

实施例一：

参照图1，层状模型各向异性射线追踪方法，包括：

步骤S101，获取模型原始信息，并对模型原始信息进行参数化得到参数值；

步骤S102，给定震源位置和射线的第一相速度方向，根据参数值计算射线的相速度矢量和群速度矢量，以得到射线层间传播信息；

步骤S103，根据射线层间传播信息判断射线是否到达模型边界；

具体地，经过参数化的模型是有几何上的边界的，此边界可能是规则的，也可能是不规则的，在三维情况下，边界更为复杂。射线到达模型边界，即经过计算，射线与参数化的目标界面的交点不存在(未到达底层界面时)，或者射线到达底层界面(不再继续向下传播)。

步骤S104，如果未到达模型边界，则计算射线经过当前界面的第二相速度方向；

步骤S105，根据第二相速度方向与界面的关系，判断射线是否能够实际继续传播，如果不能，则判断射线是否符合调整条件；

步骤S106，如果符合，则对射线进行调整得到第三相速度方向；

步骤S107，重复上述计算射线在层间传播的步骤，直至射线到达模型边界或无法继续在模型中传播，得到射线的完整路径。

根据本发明的示例性实施例，模型原始信息包括界面数据和弹性参数数据，对模型原始信息进行参数化得到参数值包括：

步骤S201，对界面数据和弹性参数数据反算得到非均匀分布的具有三次B样条性质的控制节点；

步骤S202，使用控制节点进行拟合得到参数化模型信息，其中，参数化模型信息包括基函数、节点信息和背景场；

步骤S203，使用基函数、节点信息和背景场通过三次B样条函数进行插值得到与参数化模型对应的参数值。

具体地，不均匀各向异性弹性介质由层状的横向不均匀界面分隔。每个界面可以使用不规则节点的二维三次B样条函数进行拟合，如公式(1)所示；

其中，D_0i是第i个界面的背景深度；是x_j＜x＜x_j+1和y_k＜y＜y_k+1范围内的扰动系数,相应的j＝0,1,…,J，k＝0,1,…,K；J、K代表X、Y方向的节点数量，m,n分别是j-3到j,k-3到k的整数；B_m.4和B_n,4是(x,y)点上的四阶三次B样条基函数。界面的数量由介质的不均匀性决定。我们令两个相邻介质的深度增量足够小，进而确保在层间传播的射线路径可以近似为沿直线传播。层间横向变化的弹性参数也可以使用2D三次B样条函数进行拟合，如公式(2)所示；

其中，E_i可以是P波速度、S波速度、密度或者其他任意各向异性参数。E_0i,b_mn,B_m,4(x)和B_n,4(y)分别代表第i层的弹性参数场的背景值、扰动系数、四阶三次B样条基函数。任意各项异性均能够通过上述参数化方法表达。层间的弹性参数可以是不连续的，深度方向的线性或者二次变化的弹性场能够通过若干薄层进行逼近。

根据本发明的示例性实施例，相速度矢量包括相速度大小和第一相速度方向，根据参数值确定射线的相速度矢量和群速度矢量包括：

根据方程组(3)计算群速度矢量和相速度矢量：

(Γ_ik-δ_ik)U_k＝0,Γ_ik＝a_ijklp_jp_l,a_ijkl＝C_ijkl/ρ (3)

其中，Γ为Christoffel矩阵，δ_ik为Kronecker delta函数，U_k为复振幅向量在k方向的分量，p为慢度矢量，a_ijkl为正则化的弹性模量，ρ为密度，C为弹性模量，下标i,j,k,l为爱因斯坦求和的指标；

其中，根据方程组(4)计算慢度矢量：

由此，

根据公式(5)计算相速度矢量：

v_n＝n′v_n (5)

其中，n′为相速度方向单位矢量，v_n为相速度大小，v_n为相速度矢量；

另外，

根据公式(6)计算群速度矢量：

其中，v_g为群速度矢量，g^(m)为Christoffel方程的特征值对应的特征向量，p为慢度矢量，a_ijkl为正则化的弹性模量，下标i,j,k,l为爱因斯坦求和的指标。

具体地，群速度和相速度能够通过Christoffel方程进行求解，由广义胡可定理定义的应力应变关系如公式(7)所示；

σ_ij＝C_ijklε_kl (7)

其中，σ和ε分别是应力和应变张量，C为刚度矩阵。与上述公式相同，使用了Einstein求和。

根据本发明的示例性实施例，计算射线经过当前界面的第二相速度方向包括：

步骤S301，将相速度矢量进行全局坐标系至局部坐标系的转换；

步骤S302，对转化后的相速度矢量通过广义斯内尔定理得到射线经过当前界面的第二相速度方向。

具体地，界面处的射线传播可以通过广义斯内尔定理进行求解，相速度大小v_n和相速度方向α_n的关系写作一般形式为公式(8)所示；

其中，上标u和l分别对应上下层参数,p为射线参数。对于上述公式，其参考系应为与界面相关的局部坐标系。如果全局坐标系与局部坐标系的方位角和倾角的夹角分别是和θ，相关的方向余弦矩阵R_θ如矩阵(9)和矩阵(10)所示；

通过旋转，可以将相速度从全局坐标系转换至局部坐标系中，如公式(11)所示；

同理，相速度通过式(12)也能够从局部坐标系旋转至全局坐标系中，

相速度和群速度的矢量表达式如公式(5)(13)所示，

v_g＝v_gn (13)

其中，n是平行于射线传播方向的单位向量，v_g是群速度的标量值。

根据本发明的示例性实施例，调整条件包括射线的参数条件、界面的倾角条件和入射角条件，对射线进行调整得到第三相速度方向包括：

根据方程组(13)计算第三相速度方向：

其中，f(ξ_i)为高斯函数，ξ为可调整的最大正弦值，ζ为调整后的最小正弦值，ξ_i为欲调整的入射角正弦值，ζ_i为调整后的第三相速度对应的入射角正弦值，i为当前界面的层数。

具体地，当射线遇到临近层间速度和弹性参数剧烈变化且具有较大倾角的界面时，将其交点视作新的震源，利用惠更斯理论对其进行部分调整，令射线能够继续下行。将可调整的超过临界角的正弦函数值指定为[1,ξ],对应的调整后正弦函数值范围设为[ζ,1]。为了获得更好的调整射线分布，使用了高斯函数。可选择ξ小于1.05，ζ大于0.98，具体的情况由模型复杂程度决定。

这里，调整条件包括射线的参数条件、界面的倾角条件和入射角条件，具体来说判断依据包括：

(1)射线经过的界面上下层速度及弹性参数变化剧烈；

(2)当前界面倾角过大(导致很难有射线穿透界面)；

(3)经过计算，折射射线对当前界面的入射角正弦值范围在[1,ξ]之间。

需要说明的是，图10的流程图代表了模型参数化和具体射线追踪过程，即，其为从获取模型数据到生成与某震源和模型相关的完整射线路径的流程。使用此流程可以完整得到我们需要的射线路径及射线走向、速度等相关信息。其中，射线初始化是射线初始值设定的相关内容，包括指定震源点，与相速度相关的初始入射角，是人为设定的，与射线追踪的具体方法无关。流程图中层间传播对应的是计算射线从第i-1层射入第i层界面时，相关的相速度方向、群速度方向与第i层界面交点的步骤，即步骤S102中的层间传播信息。

本发明提供了层状模型各向异性射线追踪方法，包括：获取模型原始信息，对模型原始信息进行参数化得到参数值；给定震源位置和射线第一相速度方向，根据参数值计算相速度矢量和群速度矢量，从而得到射线层间传播信息；判断射线是否到达模型边界；如果未到达模型边界，则计算射线经过当前界面的第二相速度方向；判断射线是否能够实际继续传播，如果不能，则判断射线是否符合调整条件；如果符合，则对射线进行调整得到第三相速度方向；重复上述计算射线在层间传播的步骤，直至射线到达模型边界或无法继续在模型中传播，得到射线的完整路径。本发明能够对复杂三维各向异性介质进行建模，为动力学射线追踪及振幅计算提供了帮助。

实施例二：

参照图4，层状模型各向异性射线追踪系统，包括:

获取单元10，用于获取模型原始信息，并对模型原始信息进行参数化得到参数值；

第一计算单元20，用于给定震源位置和射线的第一相速度方向，根据参数值计算射线的相速度矢量和群速度矢量，以得到射线层间传播信息；

第一判断单元30，用于根据射线层间传播信息判断射线是否到达模型边界；

第二计算单元40，用于在射线未到达所述边界的情况下，计算射线经过当前界面的第二相速度方向；

第二判断单元50，用于根据第二相速度方向与界面的关系，判断射线是否能够实际继续传播；

第三判断单元60，用于判断射线是否符合调整条件；

调整单元70，用于对射线进行调整得到第三相速度方向；

循环单元80，用于重复计算射线在层间传播的步骤，直至射线到达模型边界或无法继续在模型中传播，得到射线的完整路径。

根据本发明的示例性实施例，模型原始信息包括界面数据和弹性参数数据，获取单元10包括：

对界面数据和弹性参数数据反算得到非均匀分布的具有三次B样条性质的控制节点；

使用控制节点进行拟合得到参数化模型信息，其中，参数化模型信息包括基函数、节点信息和背景场；

使用基函数、节点信息和背景场通过三次B样条函数进行插值得到与参数化模型对应的参数值。

根据本发明的示例性实施例，相速度矢量包括相速度大小和第一相速度方向，第一计算单元10包括：

根据方程组(3)计算相速度矢量和群速度矢量，求解相速度矢量和群速度矢量的详细公式包括：

根据方程组(4)计算慢度矢量；

根据公式(5)计算相速度矢量；

根据公式(6)计算群速度矢量；

根据本发明的示例性实施例，第二计算单元40包括：

将相速度矢量进行全局坐标系至局部坐标系的转换；

对转化后的相速度矢量通过广义斯内尔定理得到射线经过当前界面的第二相速度方向。

根据本发明的示例性实施例，调整条件包括射线的参数条件、界面的倾角条件和入射角条件，调整单元60包括：

根据方程组(13)计算第三相速度方向。

具体地，实施例二中所涉及到的方程组与公式均在实施例一中示出，再次不再赘述。

需要说明的是，通过本发明实施例提供的射线追踪方法得到完整的射线路径，可以通过振幅计算方法获取相应的振幅信息，振幅信息服务于后期的与振幅相关的偏移反演、成像等应用的。因此，射线追踪方法是振幅计算的基础，也是地震射线理论的基础工作，以下为振幅计算的具体步骤：

使用基于笛卡尔坐标系的动力学射线追踪系统如公式(14)(15)所示；

dQ_i/dτ＝A_ijQ_j+B_ijP_j, (14)

dP_i/dτ＝-C_ijQ_j-D_ijP_j (15)

其中，(γ₁,γ₂,τ)代表了射线坐标系，A,B,C,D分别表示程函方程的四个二阶导

振幅计算使用公式(16)：

其中，S为炮点，R为接收点，震源子波为s(t)，F(S；γ₁,γ₂)为震源的辐射函数，L(R,S)＝|detQ(R,S)|^1/2为几何扩散系数，R^C为从S到R的完整R/T系数，T(R,S)是S到R的走时，T^c(R,S)是有关焦散的相移。

实施例三：

此层状模型建模方法能够对复杂三维各向异性介质进行建模，与传统的网格模型建模技术相比，减少了存储需求，并且三次B样条能够准确方便地给出界面相关的一阶和二阶导数，为动力学射线追踪及振幅计算提供了稳定的导数信息。

三维复杂各向异性介质的射线追踪如图5(a)(b)(c)(d)(e)(f)。其各项异性参数如表1。图5(a)为三维层状各向同性模型，纵轴对应了深度，色标对应了速度。(b)和(c)分别为使用模型(a)计算的射线在3000m，5000m深度的走时分布。(d)为三维层状各向异性模型。(e)和(f)分别为使用模型(d)计算的射线在3000m，5000m深度的走时分布。(a)(b)(c)对应了各向同性模型射线及其在某深度的走时分布。(d)(e)(f)对应了相同界面位置的TTI各向异性模型及同样深度的射线走时分布，层间倾角场与界面倾角相关。由图可知，各向异性介质的速度沿着射线传播方向变化快，射线相对炮点的偏转更快。

表1 三维层状各向异性模型参数

图6(a)(b)(c)(d)为调整射线路径在层状模型射线追踪的应用及效果。图6(a)为射线追踪的结果，实线为没有经过调整的一般射线路径，虚线为经过调整的射线段，横纵轴分别为距离和深度。(b)为PSPI脉冲响应结果。(c)为不含经过调整的射线段的脉冲响应结果。(d)为含有经过调整的射线段的脉冲响应结果。虚线为经过调整的射线路径，可见，调整以后射线能够继续向下传播，将(c)与(d)比较，可以看出，调整射线对部分深部区域起到了一定的补偿作用。分别将(c)、(d)与(b)的基于波动方程的Gazdag’s相移插值波动方程偏移(PSPI WEM)的脉冲响应结果相比，经过调整的脉冲响应结果与其更加接近。

图7(a)(b)(c)显示了介质层状建模技术对Sigsbee复杂模型的支持情况。图7(a)为使用复杂层状模型建模建立的Sigsbee层状模型，横纵轴分别为模型长度和深度。(b)为在区域1(zone 1)比较层状模型偏移结果、真实速度模型与网格模型偏移结果。(c)为在区域2(zone 2)比较层状模型偏移结果、真实速度模型与网格模型偏移结果。(a)为经过层状建模的Sigsbee模型，其中，散射点、陡倾角岩丘、尖灭、断层均能被较精确确定。(b)、(c)分别为两个区域的高斯束偏移成像(GBM)结果。左边为基于层状模型的成像结果，中间为真实模型，右边为基于网格模型的成像结果。对比(b)可见，由于网格模型经过平滑，岩丘附近速度变化剧烈，无法精确成像出岩丘边界位置，而层状模型由于其对模型的接近程度高，能够更加精确地成像出岩丘边界的位置。由(c)可以看到两种模型对于岩丘下方区域都具有较好的成像效果。

图8(a)(b)(c)(d)显示了复杂各项异性模型层状建模对于Hess模型的应用。图8(a)为P波速度；(b)为密度，(c)(d)为Thomsen’s各向异性参数Epsilon和Delta，横纵轴分别为模型长度和深度。各项异性介质需要对速度、密度、弹性参数等进行建模。通过层状建模技术能够准确保留相应的信息。

图9(a)(b)(c)(d)为基于层状模型和基于网格模型的GBM结果，图9(a)和(b)分别使用了基于网格建模的VTI和各向同性的Hess模型，横纵轴分别为长度和深度。(c)(d)分别使用了基于层状模型的VTI和各向同性Hess模型。其中(a)(b)分别是VTI和各项同性网格模型的高斯束偏移结果，(c)(d)分别是VTI和各向同性层状模型的高斯束偏移结果。由左右两组对比可知，各项异性参数对于各向异性模型成像来说具有重要的影响。通过对比上下两组结果，可知，基于层状模型的GBM结果噪声更小，且层位置更清晰。

本发明实施例提供的层状模型各向异性射线追踪系统，与上述实施例提供的层状模型各向异性射线追踪方法具有相同的技术特征，所以也能解决相同的技术问题，达到相同的技术效果。

本发明实施例所提供的层状模型各向异性射线追踪方法和系统的计算机程序产品，包括存储了程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统和装置的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (8)

1.一种层状模型各向异性射线追踪方法，其特征在于，包括：

获取模型原始信息，并对所述模型原始信息进行参数化得到参数值；所述模型原始信息包括界面数据和弹性参数数据，所述对所述模型原始信息进行参数化得到参数值包括：对所述界面数据和所述弹性参数数据反算得到非均匀分布的具有三次B样条性质的控制节点；使用所述控制节点进行拟合得到参数化模型信息，其中，所述参数化模型信息包括基函数、节点信息和背景场；使用所述基函数、所述节点信息和所述背景场通过三次B样条函数进行插值得到与所述参数化模型对应的所述参数值；

给定震源位置和射线的第一相速度方向，根据所述参数值计算所述射线的相速度矢量和群速度矢量，以得到射线层间传播信息；

根据所述射线层间传播信息判断所述射线是否到达模型边界；

如果未到达所述模型边界，则计算所述射线经过当前界面的第二相速度方向；

根据所述第二相速度方向与界面的关系，判断所述射线是否能够实际继续传播，如果不能，则判断所述射线是否符合调整条件；

如果符合，则对所述射线进行调整得到第三相速度方向；

重复上述计算射线在层间传播的步骤，直至所述射线到达所述模型边界或无法继续在模型中传播，得到所述射线的完整路径。

2.根据权利要求1所述的层状模型各向异性射线追踪方法，其特征在于，所述相速度矢量包括相速度大小和所述第一相速度方向，所述根据所述参数值确定所述射线的相速度矢量和群速度矢量包括：

根据下式计算所述相速度矢量和群速度矢量：

其中，为Christoffel矩阵，为Kronecker delta函数，为复振幅向量在方向的分量，为慢度矢量，为正则化的弹性模量，为密度，为弹性模量，下标为爱因斯坦求和的指标；

具体的计算过程为根据下式计算所述慢度矢量：

由此，

根据下式计算所述相速度矢量：

其中，为所述相速度方向单位矢量，为所述相速度大小，为所述相速度矢量；

另外，

根据下式计算所述群速度矢量：

其中，为所述群速度矢量，为所述Christoffel方程的特征值对应的特征向量，为所述慢度矢量，为所述正则化的弹性模量，下标为爱因斯坦求和的指标。

3.根据权利要求1所述的层状模型各向异性射线追踪方法，其特征在于，所述计算所述射线经过当前界面的第二相速度方向包括：

将所述相速度矢量进行全局坐标系至局部坐标系的转换；

对转化后的所述相速度矢量通过广义斯内尔定理得到所述射线经过所述当前界面的所述第二相速度方向。

4.根据权利要求1所述的层状模型各向异性射线追踪方法，其特征在于，所述调整条件包括所述射线的参数条件、界面的倾角条件和入射角条件，所述对所述射线进行调整得到第三相速度方向包括：

根据下式计算所述第三相速度方向：

其中，为高斯函数，为可调整的最大正弦值，为调整后的最小正弦值，为欲调整的入射角正弦值，为调整后的所述第三相速度对应的所述入射角正弦值，为当前界面的层数。

5.一种层状模型各向异性射线追踪系统，其特征在于，包括:

获取单元，用于获取模型原始信息，并对所述模型原始信息进行参数化得到参数值；所述模型原始信息包括界面数据和弹性参数数据，所述获取单元包括：对所述界面数据和所述弹性参数数据反算得到非均匀分布的具有三次B样条性质的控制节点；使用所述控制节点进行拟合得到参数化模型信息，其中，所述参数化模型信息包括基函数、节点信息和背景场；使用所述基函数、所述节点信息和所述背景场通过三次B样条函数进行插值得到与所述参数化模型对应的所述参数值；

第一计算单元，用于给定震源位置和射线的第一相速度方向，根据所述参数值计算所述射线的相速度矢量和群速度矢量，以得到射线层间传播信息；

第一判断单元，用于根据所述射线层间传播信息判断所述射线是否到达模型边界；

第二计算单元，用于在所述射线未到达所述边界的情况下，计算所述射线经过当前界面的第二相速度方向；

第二判断单元，用于根据所述第二相速度方向与界面的关系，判断所述射线是否能够实际继续传播；

第三判断单元，用于判断所述射线是否符合调整条件；

调整单元，用于对所述射线进行调整得到第三相速度方向；

循环单元，用于重复上述计算射线在层间传播的步骤，直至所述射线到达所述模型边界或无法继续在模型中传播，得到所述射线的完整路径。

6.根据权利要求5所述的层状模型各向异性射线追踪系统，其特征在于，所述相速度矢量包括相速度大小和所述第一相速度方向，所述第一计算单元包括：

根据下式计算所述相速度矢量和群速度矢量：

其中，为Christoffel矩阵，为Kronecker delta函数，为复振幅向量在方向的分量，为慢度矢量，为正则化的弹性模量，为密度，为弹性模量，下标为爱因斯坦求和的指标；

具体的计算过程为根据下式计算所述慢度矢量：

由此，

根据下式计算所述相速度矢量：

其中，为所述相速度方向单位矢量，为所述相速度大小，为所述相速度矢量；

另外，

根据下式计算所述群速度矢量：

其中，为所述群速度矢量，为所述Christoffel方程的特征值对应的特征向量，为所述慢度矢量，为所述正则化的弹性模量，下标为爱因斯坦求和的指标。

7.根据权利要求5所述的层状模型各向异性射线追踪系统，其特征在于，所述第二计算单元包括：

将所述相速度矢量进行全局坐标系至局部坐标系的转换；

对转化后的所述相速度矢量通过广义斯内尔定理得到所述射线经过所述当前界面的所述第二相速度方向。

8.根据权利要求5所述的层状模型各向异性射线追踪系统，其特征在于，所述调整条件包括所述射线的参数条件、界面的倾角条件和入射角条件，所述调整单元包括：

根据下式计算所述第三相速度方向：

其中，为高斯函数，为可调整的最大正弦值，为调整后的最小正弦值，为欲调整的入射角正弦值，为调整后的所述第三相速度对应的所述入射角正弦值，为当前界面的层数。