CN103675894A - 一种基于三维高斯束射线追踪和频率域合成地震记录方法 - Google Patents

Abstract

本发明涉及一种基于三维高斯束射线追踪和频率域合成地震记录方法，其包括以下步骤：1）利用普通射线追踪的方法进行运动学追踪，获取射线路径和运动学特征；2）利用高斯射线束进行动力学追踪，计算出动力学特征；3）在频域进行合成地震记录。由于频率域进行合成地震记录可以考虑基于Futterman方程的衰减算法，可以考虑基于频率的格林函数，也可以考虑耦合射线理论的弱各向异性模型，所以可以较常规方法更加精细的描述地震波场的动力学特征，使利用射线理论合成的地震记录接近于波动理论的精度。本发明可以方便地计算射线振幅、相位移等信息，可有效地应用于合成记录制作、保幅叠前深度偏移、波形层析反演等领域。

Description

一种基于三维高斯束射线追踪和频率域合成地震记录方法

技术领域

本发明涉及石油地球物理勘探领域，特别是关于一种基于三维高斯束射线追踪和频率域合成地震记录方法。

背景技术

地震波场数值模拟经过多年的发展，目前已经取得了大量的研究成果，主要包括两大类，即波动方程数值模拟方法和射线理论正演方法。波动方程数值模拟方法能够揭示波在介质中传播的振幅、频率、相位变化，真实地反映出波的动力学特征，是进行波场研究的有效手段，但是占用计算机内存大，耗费机时多，有时对复杂地质构造的情况亦难于研究。射线理论正演方法根据模型结构及速度变化，追踪地震波在地下介质中的运动轨迹及波经过速度界面发生反射、透射变化，并通过计算波场扩散、透射系数、反射系数等波的能量损失和相位变化近似模拟波场特征，射线理论正演方法主要用于研究波传播的轨迹、反射点的分布、波场的分布范围等波的运动学特征，也是进行层析反演的基础。目前射线追踪可以分为运动学射线追踪和动力学射线追踪两类。运动学射线追踪需要计算射线路径、波前面、走时信息，动力学射线追踪在此基础上还需要计算射线振幅、相位移等信息。运动学射线追踪是动力学射线追踪的基础。射线追踪方法在地球物理层析成像中有着极其重要的作用，它是研究介质任意速度分布情况下地震波传播问题的有效方法之一。射线追踪的路径精度和追踪计算的速度直接决定着成像的质量和计算速度。因此，研究快速而精确的射线追踪方法，对于层析成像问题来说有着特别重要的实际意义。

现有技术中的众多求解射线路径和走时的方法中，传统的方法为射线级数解法，此方法是由美国的Karal与Keller(1959)将电磁学的研究成果引入弹性波领域、并导出弹性动力学运动方程的射线级数解的过程中率先提出的，它是渐近射线法的前身。近年来发展了程函方程法、波前重建法、最短路径法、模拟退火法、遗传算法等新方法。现有技术中的射线追踪的具体方法很多，通常使用试射法（打靶法）、弯曲法和波前法，其中，试射法利用Snell定律确定射线路径和接收点，不断的调整射线的入射角，使射线出射位置与检波点之间不断靠近，最终满足给定的精度要求为止，从而实现两点间的射线追踪计算，是解决两点追踪问题的一种重要方法。试射法需要尝试多次初始试射，通过不断迭代修改射线入射角来获得准确的射线路径，比较耗时，在三维情况下尤为突处。弯曲法是固定源点和接收点，基于最小射线旅行时的Fermat原理通过多次迭代得到正确的射线路径，对一般介质，弯曲法其实是求解非线性最优化问题，计算过程中有时会陷入局部收敛，得到非全局的最优解。试射法与弯曲法都可以求解模型多路径走时，通常用于走时层析成像、合成记录制作等方面。波前法，基于惠更斯(Huygens-Fresnel)原理，地震波在传播过程中会形成一个波阵面，波阵面作为新的震源发射射线，经过多个次级源，最终传播到检波点，首先要把介质划分成许多网格点，要求射线必须经过这些网格点，源点和检波点也分别处于网格点上，射线从源点所处的网格点出发，经由各网格点以最小走时到达检波点，射线经过的这些网格点就形成了最小走时射线路径，当网格点数量增加时，计算量也将成比例增加。

20世纪90年代初随着叠前深度偏移逐渐实现工业化，上述几种方法的很多缺陷也逐渐暴露出来，主要有：①难于处理介质中速度变化剧烈的情况；②求取多值走时中的全局最小走时困难；③计算效率较低；④阴影区和复杂构造区射线覆盖密度低。因此，近十几年来射线追踪方法主要围绕上述问题进行工作。Vidale(1988，1990)和Podvin等人(1991)则从程函方程出发，通过先求出时间场分布再计算时间场的最快下降方向的办法，得到每一接收点到震源的射线路径。随后，Qin等人(1992)对Vidale(1988)的方法作了改进，提出了波前扩展方法，该方法虽然减少了原方法的不稳定性，但是却使计算量巨增。

发明内容

针对上述问题，本发明的目的是提供一种基于三维高斯束射线追踪和频率域合成地震记录方法，能够获得精细的动力学波场信息。

为实现上述目的，本发明采取以下技术方案：一种基于三维高斯束射线追踪和频率域合成地震记录方法，其包括以下步骤：1）利用普通射线追踪的方法进行运动学追踪，获取射线路径和运动学特征，即根据地震波的传播规律确定地震波在实际地层中传播的射线路径；2）利用高斯射线束进行动力学追踪，计算出动力学特征；3）在频域进行合成地震记录的计算，具体过程为：3.1）对于某一种类型的波，计算从震源开始，由射线参数

规定的充分密集的射线***，并从初值或两点射线追踪及动态射线追踪的结果中取出与每一射线相应的量A(O_s)、q(O_s)、p(O_s)、v(O_s)并确定

其中，A(O_s)是射线O_s处的复值振幅因子，q(O_s)和p(O_s)是射线在O_s处的方向因子，控制反射和透射波的传播方向，v(O_s)射线在O_s处的速度，

是高斯射线的权函数；3.2）对同一条中心射线分别确定不同检波点S处的射线中心坐标(s,n)，同时确定相因子T(S,O_s)，并取其实部和虚部；3.3）给定频率间隔Δf，在某一频率范围f₁-f₂内计算相应与时间Re[T(S,O_s)]的离散频谱，并储存在数组中；3.4）在3.2）做完后，得到的是某一高斯射线束对所有有关检波点的贡献，分别以离散频谱的形式储存在有关各道中，然后返回到步骤3.1）中计算另一高斯射线束对各道的贡献，把结果叠加起来，直至同一种基波高斯射线束作完；接着再计算进行另一种基波的波场，直到把所有的基波高斯射线束都作完，所得的是各道反射系数序列在频域中的结果，以离散谱的形式储存于各道中；3.5）在所有基波高斯射线束作完后，输入时间信号；3.6）根据输入信号的形式，以等时间间隔Δt在某一时间范围内对其进行取样，得到离散时间信号，进行快速傅氏变换后得到信号离散频谱；3.7）将步骤3.6）中的信号离散频谱分别与步骤3.3）中各道反射系数序列中的离散频谱相乘后，再进行反傅氏变换IFFT，得到时域形式的合成地震记录。

所述步骤2）包括以下步骤：2.1）通过射线追踪计算沿射线的旅行时；2.2）确定射线的偏振矢量；2.3）利用动态射线追踪，计算出动力学特征；2.4）确定复矢量方向上的振幅衰减；2.5）确定焦散量。

所述步骤3.5）中输入信号的形式采用Gabor信号、Berlage信号、Muler信号和Ricker信号中的一种。

本发明由于采取以上技术方案，其具有以下优点：1、本发明首先利用普通射线追踪的方法进行运动学追踪，获取射线路径和运动学特征；然后利用高斯射线束进行动力学追踪，计算出动力学特征；最后在计算出动力学特征的基础上在频率域进行合成地震记录的计算，由于频率域进行合成地震记录可以考虑基于Futterman方程的衰减算法，可以考虑基于频率的格林函数，也可以考虑耦合射线理论的弱各向异性模型，所以与常规方法相比能够更加精细的描述地震波场的动力学特征，使利用射线理论合成的地震记录接近于波动理论的精度。2、本发明利用经典的旁轴射线追踪***或者动力学射线追踪***，可以方便地计算射线振幅、相位移等信息，有效地应用于合成记录制作、保幅叠前深度偏移、波形层析反演等领域，对于模型中存在焦散和阴影的情况，可以用它的改进形式-高斯束射线理论来计算异常波场区的振幅。3、目前，波动问题的数值解法分为两大类，一类是以弹性波动力学方程数值解为基础的直接求解方法，如有限差分法、有限元法等，另一类是以不均匀介质波动方程取高频近似为基础的近似求解的射线方法，两者相比，前者模拟的波场较完整、精度高，适合于尺度小于优势波场的介质模型，但是往往需要大型计算机，并且要耗费往往使人难以忍受的计算时间，且波场复杂，不便于识别和解释。波动方程取高频近似的射线法尽管不如波动方程解精确，但是适合于尺度大于优势波长的模型，且解法速度快，波场易于识别和解释。本发明是在常规射线追踪的进一步发展，研究表明常规的旁轴射线追踪并非完美无缺，在用于复杂模型时，仍然遇到困难，比如焦散区、阴影区等非正则区域，旁轴射线追踪法将失效，其次，旁轴射线追踪常常需要做麻烦的两点射线追踪，计算效率低。本发明是在动态射线追踪的基础上进行，在求解动态射线追踪方程组时，在求解中心射线的几何扩展的同时，还可以求取中心射线邻近点上的旅行时、旅行时二阶导数、几何扩展和射线振幅，同时追踪和确定中心射线邻近的其它射线的位置，在中心射线坐标系下，将弹性波动力学波动方程取高频近似化为抛物线方程，导出动态射线追踪方程和传输方程，对抛物线方程取时间调和解得到高斯射束，将高频地震波场展开或分解成高斯射束***表示的波场。高斯射束近似于旁轴近似射线理论相比的主要优点是：一是可以消除焦散区等失效区的奇异性，二是可以避开两点射线追踪，大大提高计算效率。在高斯束射线追踪计算出动力学特征的基础上，在频率域考虑基于Futterman方程的衰减算法，基于频率的格林函数，实现考虑耦合射线理论的弱各向异性模型，所以与常规方法相比能够更加精细的描述地震波场的动力学特征，使利用射线理论合成的地震记录接近于波动理论的精度。本发明可以广泛地应用于合成记录制作、保幅叠前深度偏移、波形层析反演等领域。

附图说明

图1是本发明的观测***及射线路径示意图；

图2是本发明的Z分量合成地震记录（P波）；

图3是本发明的X2分量合成地震记录（P波）；

图4是本发明的Z分量合成地震记录（P＋SV波）；

图5是本发明的X2分量合成地震记录（P＋SV波）；

图6是本发明的Z分量合成地震记录（P＋SV+多次直达波）；

图7是本发明的X2分量合成地震记录（P＋SV+多次直达波）；

图8是本发明的三维射线路径显示示意图；

图9是本发明的非零偏不同波场类型的射线路径显示示意图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

如图1～9所示，本发明的基于三维高斯束射线追踪和频率域合成地震记录方法，包括以下步骤：

1、利用普通射线追踪的方法进行运动学追踪，获取射线路径和运动学特征，即根据地震波的传播规律确定地震波在实际地层中传播的射线路径。

在地震学中，有两类地震射线追踪问题：一类是一点射线追踪，即已知射线初始点（源点）和初始出射方向，求解地震波的传播路径，又叫初值射线追踪问题；另一类是两点射线追踪，即已知射线初始点（源点）和另一个观察点（接收点）的位置，射线初始出射方向未知，求解两点之间的射线路径。本发明的射线追踪既需要研究初值射线追踪，也需要研究两点射线追踪。

由程函方程可以导出射线追踪方程：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{d\τ}}=v^2{p}_x,\frac{\mathrm{dy}}{\mathrm{d\τ}}{v}^2{p}_y,\frac{\mathrm{dz}}{\mathrm{d\τ}}=v^2{p}_z\\ \frac{d{p}_x}{\mathrm{d\τ}}=-\frac{1}{v}\frac{\∂v}{\∂x},\frac{d{p}_y}{\mathrm{d\τ}}=-\frac{1}{v}\frac{\∂v}{\∂y},\frac{d{p}_z}{\mathrm{d\τ}}=-\frac{1}{v}\frac{\∂v}{\∂z}\end{array}\right.---\left(1.1\right)$

式中， $p_x=\∂\mathrm{\τ}/\∂x=\cos \mathrm{\α}/v,$ $p_y=\∂\mathrm{\τ}/\∂y=\cos \mathrm{\β}/v,$ $p_z=\∂\mathrm{\τ}/\∂z=\cos \mathrm{\γ}/v,$ 其中，p_x、p_y、p_z是射线出射方向，其中，α是出射射线和x轴的夹角、β是出射射线和y轴的夹角，γ是出射射线和z轴的夹角，v代表介质的速度。

初始条件可以表示如下：即当τ＝τ₀时，有：

$\left\{\begin{array}{c}x=x_0,y=y_0,z=z_0\\ p_x=p_{x0},p_y=p_{y0},p_z=p_{z0}\end{array}\right.---\left(1.2\right)$

式中，初始条件{x₀,y₀,z₀}确定震源的位置，{p_x0,p_y0,p_z0}给定射线的初始出射方向。

2、利用高斯射线束进行动力学追踪，计算出动力学特征，具体过程：

2.1）通过射线追踪计算沿射线的旅行时；

假设三维速度模型是定义在一个空间体M内：

$M:x_{\min }^i\≤x^i\≤x_{\max }^i(i=\mathrm{1,2,3})$

式中，

和

表示三维模型的边界坐标，通过某种函数关系，可以给定三维空间介质的参数，例如纵、横波速度V_p，V_s，密度ρ和衰减因子

本发明的射线追踪有以下几步：(a)射线追踪和计算沿射线的旅行时；(b)确定射线的偏振矢量；(c)动态射线追踪（例如计算射线的传播模型）；(d)确定复矢量方向上的振幅衰减。

射线追踪由射线的计算组成，例如沿射线逐点计算坐标。射线的参数由独立变量确定：

$\mathrm{\σ}={\mathrm{\σ}}_0+{\∫}_{{\mathrm{\τ}}_0}^{\mathrm{\τ}}{v}^{\mathrm{NEXPS}}\mathrm{d\τ}={\mathrm{\σ}}_0+{\∫}_{s_0}^s{v}^{\mathrm{NEXPS}-1}\mathrm{ds}---\left(2.1\right)$

式中，τ是旅行时，σ代表的是某条射线，σ₀代表的是某条射线的起始位置。s是沿射线的弧长，v是相应传播类型的波的速度，NEXPS整数指数可以在使用时根据需要进行确定，NEXPS=0,±1,±2,.......。我们认为射线追踪***是六个一阶常微分方程组的形式：

$\left.\begin{array}{c}\frac{d{x}^i}{\mathrm{d\σ}}=v^{2-\mathrm{NEXPS}}{G}^{\mathrm{ij}}{p}_j,\\ \frac{d{p}_i}{\mathrm{d\σ}}=v^{2-\mathrm{NEXPS}}(-v^{-3}\frac{\∂v}{\∂x^i}+{\mathrm{\Γ}}_{\mathrm{ij}}^k{G}^{\mathrm{jl}}{p}_k{p}_l).\end{array}\right\}---\left(2.2\right)$

式中，xⁱ是沿射线上点的坐标，p_i是慢度向量在点xi的协变分量，并且

其中度量张量G^ij和克里斯托弗尔符号

定义如下：对于右旋曲线坐标***（x¹，x²，x³），坐标***的局部特性可以由协变分量G_ij＝G_ij(x^k)和逆变分量G^ij＝G^ij(x^k)来描述，其中i,j,k＝1,2,3，克里斯托弗尔符号

在界面上射线按照斯奈尔定律传播。

方程：

$\mathrm{\τ}={\mathrm{\τ}}_0+{\∫}_{{\mathrm{\σ}}_0}^{\mathrm{\σ}}{v}^{-\mathrm{NEXPS}}\mathrm{d\σ}---\left(2.3\right)$

旅行时的实部（2.3）可以由公式（2.1）推导获得，旅行时的虚部由下述公式计算得出：

${\mathrm{\τ}}^{\mathrm{IM}}=\frac{1}{2}{t}^{*}={\∫}_{{\mathrm{\τ}}_0}^{\mathrm{\τ}}{\left(2Q\right)}^{-1}{v}^{-\mathrm{NEXPS}}\mathrm{d\σ}---\left(2.4\right)$

式中，Q是相关品质因素。

2.2)确定射线的偏振矢量；

任意地选择一条射线用Ω表示，确定沿射线Ω的偏振矢量在多种理论下的完全射线追踪都起到很重要的作用。如果要确定沿Ω传播的波的位移矢量的方向，尤其是在S波（横波）情况下，就必须知道偏振矢量。另外，偏振矢量作为与Ω连接的射线中心坐标系矢量的基点，这个射线中心坐标系在很多情况下很有用，尤其在动态射线追踪时。

Popov和Psencik(1978a,b)将正交的射线中心坐标系引入地震学中。本发明将(q¹,q²,q³＝s)定义这个射线中心坐标，其中，s是沿射线Ω的弧长，q^k是原点在射线Ω上的笛卡尔坐标平面，这个平面正切射线Ω上点s＝q³的波前面。由三个单位矢量e₁,e₂,e₃（偏振矢量）组成射线中心坐标系的基向量，其中，e₃正切射线Ω，e₁和e₂垂直射线Ω，v是地震波传播的速度，这种方式引入是为了使射线中心坐标系正交。

方程：

$\frac{d{H}_{\mathrm{ij}}}{\mathrm{d\σ}}=v^{1-\mathrm{NEXPS}}[\frac{\∂v}{\∂x^k}{G}^{\mathrm{kl}}{H}_{\mathrm{lj}}{p}_i-p_k{G}^{\mathrm{kl}}{H}_{\mathrm{lj}}\frac{\∂v}{\∂x^i}v{\mathrm{\Γ}}_{\mathrm{ik}}^m{G}^{\mathrm{kl}}{p}_l{H}_{\mathrm{mj}}]---\left(2.5\right)$

式中，G是度量张量、Γ是克里斯托弗尔符号，H_ij是射线中心基向量e_j的协变分量，(e_j)_i＝H_ij。

方程（2.5）简化为

$\frac{d{H}_{\mathrm{iJ}}}{\mathrm{d\σ}}=v^{1-\mathrm{NEXPS}}[\frac{\∂v}{\∂x^k}{G}^{\mathrm{kl}}{H}_{\mathrm{iJ}}{p}_i+v{\mathrm{\Γ}}_{\mathrm{ik}}^m{G}^{\mathrm{kl}}{p}_l{H}_{\mathrm{mJ}}]---\left(2.6\right)$

矢量e_J垂直射线Ω。当单位矢量e₃正切于射线Ω(j=3在公式(2.5)中)时，

$H_{i3}=v{p}_i=G_{\mathrm{ij}}\frac{d{x}^j}{\mathrm{ds}}---\left(2.7\right)$

这样公式（2.5）就等价于射线追踪***（2.2）中的第二方程 $d{p}_i/\mathrm{d\σ}=v^{2-\mathrm{NEXPS}}({-v}^{-3}\mathrm{dv}/x^i+{\mathrm{\Γ}}_{\mathrm{ij}}^k{G}^{\mathrm{jl}}{p}_k{p}_l),$ 它足够用来唯一计算矢量e₁的协变分量H_i1(i＝1,2,3)，因为在用公式（2.2）进行射线追踪时，e₃取决于公式（2.7），而e₂总垂直于e₁和e₃：

H_i2＝ε_ijkG^jmH_m3G^knH_n1[det(G^rs)]^-1/2          （2.8）

式中，ε₁₂₃＝ε₂₃₁＝ε₃₁₂＝1，ε₁₃₂＝ε₃₂₁＝ε₂₁₃＝-1，ε_ijk＝0适合其他所有的ijk组合。在界面上，射线中心坐标系绕着垂直于入射平面的矢量旋转，这样旋转后，按照理论，矢量e₃将于产生的正切射线的波重合。

2.3）利用动态射线追踪，计算出动力学特征；

本发明中的动力学特征可以是指射线的传播模型，动态射线追踪由四个沿射线Ω的一阶线性微分方程组的解组成，分量是4×1的矩阵W，

$\frac{\mathrm{dW}}{\mathrm{d\σ}}=v^{2-\mathrm{NEXPS}}\mathrm{SW}---\left(2.9\right)$

其中

$S=\left[\begin{array}{cc}0& I\\ {-v}^{-2}V& 0\end{array}\right]$

0和I是0和1的1×1矩阵，V是一个2×2分量为：

$V_{\mathrm{IJ}}=(\frac{{\∂}^{2}v}{\∂x^k\∂x^l}-{\mathrm{\Γ}}_{\mathrm{kl}}^m\frac{\∂v}{\∂x_m})G^{\mathrm{kr}}{G}^{\mathrm{ls}}{H}_{\mathrm{rI}}{H}_{\mathrm{sJ}}---\left(2.10\right)$

的矩阵。W的前两个分量是旁轴射线的射线中心坐标(q¹,q²)（一阶泰勒展开），W的其他的两个分量对应的是旁轴射线慢矢量的射线中心分量。

动态射线追踪公式（2.9）有四个线性的独立解，采用Π(σ,σ₀)作为在初始条件射线的初始点σ＝σ₀，Π(σ,σ₀)＝I下的线性独立解的基础4×4矩阵。这里I是4×4的单位阵，基础矩阵Π(σ,σ₀)也叫射线传播模型，或者是动态射线追踪***的传播模型。沿射线方向它满足下面关系：

${\mathrm{\Π}}^T\mathrm{\Σ\Π}=\mathrm{\Σ},\mathrm{with\Σ}=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -1& 0& 0& 0\\ 0& -1& 0& 0\end{array}\right]---\left(2.11\right)$

这也叫Π的耦合性。

动态射线追踪***公式（2.9）的任意解W(σ)都可以用以下形式表示：

W(σ)＝Π(σ,σ₀)W(σ₀)             (2.12)

射线追踪***公式（2.9）的解在地震学中有很重要的应用。假定计算的射线符合射线***的两个参数γ¹，γ²。这样，动态射线追踪系可以用来确定沿射线方向的2×2矩阵Q^R和P^R，矩阵单元是：

$\left.\begin{array}{c}{Q}_{\mathrm{IJ}}^R={\left[\frac{\∂q^I}{\∂{\mathrm{\γ}}^J}\right]}_{q^1=q^2=0},\\ P_{\mathrm{IJ}}^R={\left[\frac{{\∂}^2\mathrm{\τ}}{\∂q^I\∂{\mathrm{\γ}}^J}\right]}_{q^1=q^2=0}={\left[\frac{\∂p_I^{\left(q\right)}}{\∂{\mathrm{\γ}}^J}\right]}_{q^1=q^2=0}\end{array}\right\}---\left(1.13\right)$

式中，q^I是射线参数γ¹,γ²确定的近轴射线的射线中心坐标；

是相应近轴慢矢量的射线中心分量，

矩阵Q^R也叫做几何扩散模型，|detQ^R|是几何扩散。它表示射线通道的扩展和收缩。几何扩散依赖于射线的参数。

引入射线坐标(γ¹,γ²,γ³)，γ¹和γ²是射线参数，γ³是沿射线的参数（如弧长s）。这样矩阵Q^R也可以理解为从射线坐标(γ¹,γ²)到沿射线Ω的射线中心坐标(q¹,q²)的转换矩阵。相似的，矩阵P^R可理解为从射线坐标(γ¹,γ²)到相位空间坐标

的转换矩阵，此4×2的矩阵X^R定为：

$X^R=\left[\begin{array}{c}{Q}^R\\ P^R\end{array}\right]---\left(2.14\right)$

满足沿射线Ω方向

$\frac{d{X}^R}{\mathrm{d\σ}}=v^{2-\mathrm{NEXPS}}{\mathrm{SX}}^R---\left(2.15\right)$

的动态射线追踪***。公式(2.15)的解可以表示成以下形式：

X^R(σ)＝Π(σ,σ₀)X^R(σ₀)            (2.16)

如果Q^R和P^R是已知的，也可以确定旅行时场M^R(σ)的二阶偏导数的2×2矩阵，矩阵元素为：

$M_{\mathrm{IJ}}^R={\left[\frac{{\∂}^2\mathrm{\τ}}{\∂q^I\∂q^J}\right]}_{q^1=q^2=0}---\left(2.17\right)$

Q^R,P^R和M^R的关系是

M^R＝P^R(Q^R)^-1            (2.18)

近界面的动态射线追踪***的转换解由近轴射线的斯奈尔定律确定。

2.4)确定复矢量方向上的振幅衰减；

矢量衰减的振幅：

$A_i=A_i^R{\left(\mathrm{\ρv}\right|\det Q^R\left|\right)}^{1/2}---\left(2.19\right)$

是在射线中心坐标系表示的射线振幅（包括由于频散引起的相位偏移）的复数矢量偏移乘以波阻抗（密度ρ和速度ν的积）和几何扩散|detQ^R|乘积的平方根。假定矢量衰减的振幅A_i在射线的起始点是单位一。如果这个起始点位于自由表面或者任意内部界面，反演系数包括在A_i里面。

定义3×3矩阵A_ij由对应于S波在起始点e₁方向上的偏振的矢量衰减的振幅A_i1、对应于S波在起始点e₂方向上的偏振的矢量衰减的振幅A_i2和对应于起始点P波的矢量衰减振幅A_i3组成。矢量衰减振幅的分量A_ij是在射线中心坐标系表示的。在单一的物块内部的射线单元是常数，除了频散，因为这时根据频散的类型它们要乘上-i或-1。在界面上，矢量衰减的振幅可以变换通过降低反射/透射（R/T）系数。

$R^E=R{\left[\frac{\widetilde{\mathrm{\ρ}}\tilde{v}\left|\cos \widetilde{\mathrm{\α}}\right|}{\mathrm{\ρv}\left|\cos \mathrm{\α}\right|}\right]}^{1/2}---\left(2.20\right)$

式中，R是平面波反射/透射系数。这些带有共轭符号的量对应于相关的反射/透射波，那些没有共轭符号的量对应入射波；在入射点的两种情况下，α和是反射/透射波各自的入射角，其中，

表示的是介质的阻抗。

2.5）确定焦散量；

在计算的射线上计算矢量振幅的衰减需要焦散点的信息，这个焦散点不具有单一射线的性质，但是具有射线场的信息。它依赖于被计算射线和偏振射线的相互位置。如果被计算射线在起始点σ＝σ₀的旅行时场的二阶偏导的2×2模型M^INIT＝M^R(σ₀)是已知的，那么他们的相互位置可以确定。这模型M^INIT可以按照矩阵Q^INIT＝Q^R(σ₀)和P^INIT＝P^R(σ₀)确定为M^INIT＝P^INIT(Q^INIT)^-1。

3、在频域进行合成地震记录；

由于合成地震记录是在射线追踪和动态射线追踪完成之后进行的，因而合成中所需要的射线上任意点O_s处的量如τ(O_s)、M(O_s)和A(O_s)等对每条中心射线都是已知的，由这些已知的量利用下式：

即可合成地震记录。式中

由已知的

确定，N是参与叠加的射线条数，权函数为：

且

$\mathrm{Re}{\{-i[M\left(O_s\right)-M^R\left(O_s\right)\left]\right\}}^{\frac{1}{2}}>0---\left(3.3\right)$

q^R(O_s)、M^R(O_s)是q(O_s)和M(O_s)的实部，都可以由动态射线追踪确定。由于射线追踪的结果中给出的是位移矢量的分量，3-D情况下，U(S,ω)＝{U_x,U_y,U_z}，所以，合成必须对U_x(S,ω)、U_y(S,ω)、U_z(S,ω)分别进行，这样就可以得到x、y、z三分量的地震记录，具体的过程为：

3.1）对于某一种类型的波，计算从震源开始，由射线参数

规定的充分密集的射线***，并从初值（或两点）射线追踪及动态射线追踪的结果中取出与每一射线相应的量A(O_s)、q(O_s)、p(O_s)、v(O_s)并确定

其中A(O_s)是射线O_s处的复值振幅因子，q(O_s)和p(O_s)是射线在O_s处的方向因子，控制反射和透射波的传播方向，v(O_s)射线在O_s处的速度，τ(O_s)是射线在O_s处的旅行时，

是高斯射线的权函数；

3.2）对同一条中心射线分别确定不同检波点S处的射线中心坐标(s,n)，同时确定相因子T(S,O_s)并取其实部和虚部，实部是波的到达时间，虚部是垂直于中心射线方向上检波点S处振幅的衰减因子。

上述中心坐标(s,n)是确定附加分量的振幅时必须知道的：

n(S,O_s)＝[x(S)-x(O_s)]t_z(O_s)-[z(s)-z(O_s)]t_x(O_s)             (3.4)

s(S,O_s)＝[x(S)-x(O_s)]t_x(O_s)-[z(s)-z(O_s)]t_z(O_s)              (3.5)

式中，t_x(O_s)，t_z(O_s)是在点O_s处与射线相切的单位矢量t的x分量和z分量。

相因子：

$T(S,O_s)=\mathrm{\τ}\left(O_s\right)+v^{-1}\left(O_s\right)I^T\left(O_s\right)X(S,O_s)+\frac{1}{2}{X}^T(S,O_s)W\left(O_s\right)X(S,O_s)---\left(3.6\right)$

其中： $X(S,O_s)=\left[\begin{array}{c}x\left(S\right)-x\left(O_s\right)\\ z\left(S\right)-z\left(O_s\right)\end{array}\right]---\left(3.7\right)$

$I\left(O_s\right)=\left[\begin{array}{c}{t}_x\left(O_s\right)\\ t_z\left(O_s\right)\end{array}\right]---\left(3.8\right)$

W(O_s)＝Ω^T(O_s)A(O_s)Ω(O_s)            （3.9）

$\mathrm{\Ω}\left(O_s\right)=\left[\begin{array}{cc}{t}_z\left(O_s\right)& -t_x\left(O_s\right)\\ t_x\left(O_s\right)& t_z\left(O_s\right)\end{array}\right]---\left(3.10\right)$

$A\left(O_s\right)=v^{-2}\left(O_s\right)\left[\begin{array}{cc}{v}^2\left(O_s\right)M\left(O_s\right)& -v_n\left(O_s\right)\\ -v_n\left(O_s\right)& -v_s\left(O_s\right)\end{array}\right]---\left(3.11\right)$

$\left\{\begin{array}{c}{v}_n\left(O_s\right)=v_x\left(O_s\right)t_z\left(O_s\right)-v_z\left(O_s\right)t_x\left(O_s\right)\\ v_s\left(O_s\right)=v_x\left(O_s\right)t_x\left(O_s\right)+v_z\left(O_s\right)t_z\left(O_s\right)\end{array}\right.---\left(3.12\right)$

$M\left(O_s\right)=\frac{P\left(O_s\right)}{q\left(O_s\right)}---\left(3.13\right)$

3.3）给定频率间隔Δf，在某一频率范围f₁-f₂内计算相应与时间Re[T(S,O_s)]的离散频谱，并储存在数组中。

3.4）在3.2）做完后，得到的是某一高斯射线束对所有有关检波点的贡献，分别以离散频谱的形式储存在有关各道中，然后返回到步骤3.1）中计算另一高斯射线束对各道的贡献，把结果叠加起来，直至同一种基波高斯射线束作完，接着再计算进行另一种基波的波场，直到把所有的基波高斯射线束都作完，所得的是各道反射系数序列在频域中的结果，以离散谱的形式储存于各道中。

3.5）在所有基波高斯射线束作完后，输入时间信号，其中，除了可以人为提供离散信号形式外，可以供选择的时间信号形式有以下几种：

①Gabor信号

式中，f_m为主频f_m、γ、t₀、

都是输入参数。

②Berlage信号

f(t)＝exp{-β(t-t₀)}sin[2πf_m(t-t₀)](t-t₀)^α

式中，f_m为主频f_m、β、α、t₀都是输入参数。

③Muler信号

$\left\{\begin{array}{c}f\left(t\right)=\sin \left(\frac{\mathrm{n\&pi;t}}{t_p}\right)-\left(\frac{n}{n+2}\right)\sin \left[(n+2)\mathrm{\&pi;t}\right]t_p,(0\&le;t\&le;t_p)\\ f\left(t\right)=0,(t<0,t>t_p)\end{array}\right.$

式中，t_p和n是输入参数。

④Ricker信号

f(t)＝{1-2[β(t-t₀)]²}exp{-[β(t-t₀)]²}

式中，β、t₀是输入参数。

3.6）根据输入信号的形式，以等时间间隔Δt（也就是地震道的采样间隔）在某一时间范围内对其进行取样，得到离散时间信号，进行快速傅氏变换后得到信号离散频谱。

3.7）将步骤3.6）中的信号离散频谱分别与步骤3.3）中各道反射系数序列中的离散频谱相乘后，再进行反傅氏变换IFFT，得到时域形式的合成地震记录。

上述各实施例仅用于说明本发明，其中实施方法的各步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (3)

1.一种基于三维高斯束射线追踪和频率域合成地震记录方法，其包括以下步骤：

1）利用普通射线追踪的方法进行运动学追踪，获取射线路径和运动学特征，即根据地震波的传播规律确定地震波在实际地层中传播的射线路径；

2）利用高斯射线束进行动力学追踪，计算出动力学特征；

3）在频域进行合成地震记录的计算，具体过程为：

3.1）对于某一种类型的波，计算从震源开始，由射线参数

规定的充分密集的射线***，并从初值或两点射线追踪及动态射线追踪的结果中取出与每一射线相应的量A(O_s)、q(O_s)、p(O_s)、v(O_s)并确定

其中，A(O_s)是射线O_s处的复值振幅因子，q(O_s)和p(O_s)是射线在O_s处的方向因子，控制反射和透射波的传播方向，v(O_s)射线在O_s处的速度，

是高斯射线的权函数；

3.2）对同一条中心射线分别确定不同检波点S处的射线中心坐标(s,n)，同时确定相因子T(S,O_s)，并取其实部和虚部；

3.3）给定频率间隔Δf，在某一频率范围f₁-f₂内计算相应与时间Re[T(S,O_s)]的离散频谱，并储存在数组中；

3.4）在3.2）做完后，得到的是某一高斯射线束对所有有关检波点的贡献，分别以离散频谱的形式储存在有关各道中，然后返回到步骤3.1）中计算另一高斯射线束对各道的贡献，把结果叠加起来，直至同一种基波高斯射线束作完；接着再计算进行另一种基波的波场，直到把所有的基波高斯射线束都作完，所得的是各道反射系数序列在频域中的结果，以离散谱的形式储存于各道中；

3.5）在所有基波高斯射线束作完后，输入时间信号；

3.6）根据输入信号的形式，以等时间间隔Δt在某一时间范围内对其进行取样，得到离散时间信号，进行快速傅氏变换后得到信号离散频谱；

3.7）将步骤3.6）中的信号离散频谱分别与步骤3.3）中各道反射系数序列中的离散频谱相乘后，再进行反傅氏变换IFFT，得到时域形式的合成地震记录。

2.如权利要求1所述的一种基于三维高斯束射线追踪和频率域合成地震记录方法，其特征在于：所述步骤2）包括以下步骤：

2.1）通过射线追踪计算沿射线的旅行时；

2.2）确定射线的偏振矢量；

2.3）利用动态射线追踪，计算出动力学特征；

2.4）确定复矢量方向上的振幅衰减；

2.5）确定焦散量。

3.如权利要求1或2所述的一种基于三维高斯束射线追踪和频率域合成地震记录方法，其特征在于：所述步骤3.5）中输入信号的形式采用Gabor信号、Berlage信号、Muler信号和Ricker信号中的一种。

Images