CN105988135A - 基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法 - Google Patents

Abstract

本发明提供一种基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，该基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法包括：步骤1，建立起伏地表的三维地质模型；步骤2，进行运动学射线追踪；步骤3，进行动力学射线追踪；步骤4，在三维介质中用投影菲涅尔带椭圆来约束高斯束在射线终点处能量的分布；以及步骤5，采用高斯波包法合成地震记录。本方法中的基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法在中深部地层具有更好的保幅性，对起伏地表和复杂模型具有更好的适应性。

Description

基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法

技术领域

本发明涉及地震资料处理，特别是涉及到一种基于投影菲涅尔带的三维起伏地表高斯束演模拟方法。

背景技术

近年来，国内外许多学者对起伏地表下的正演模拟和高斯束正演模拟进行了很多研究。首先，在不同介质中应用方面：Cerveny(1982)在引入高斯束理论后，将其应用于单层二维非均匀介质、侧向变速的层状介质、三维弹性非均匀介质以及三维弹性侧向变速层状介质的正演模拟中，Muler(1988)将其应用于二维非均匀吸收介质的正演模拟中，徐盛岩等(1988)将高斯束理论用于二维粘弹性介质地震波场的模拟。其次，在高斯束参数选择方面：Cerveny(1982)给出了最初的动力学追踪初始参数；Muler(1984)给出了二维介质中几种初始参数的选择，并讨论了其对高斯束性质的影响；George(1987)对束参数的选择给出了一些改进，确保了合成记录的稳定性。Cruz(2012)首先提出在二维介质中利用投影菲涅尔带半径来约束高斯束的有效半宽度，并给出了相应的初始参数，但没有给出投影菲涅尔带半径的计算方法。近年来，基于起伏地表的三维高斯束正演方法还处于起步阶段。为此我们发明了一种新的基于投影菲涅尔带的三维起伏地表高斯束演模拟方法，解决了以上技术问题。

发明内容

本发明的目的是提供一种针对起伏地表条件下的三维高斯束正演模拟的基于投影菲涅尔带的三维起伏地表高斯束演模拟方法。

本发明的目的可通过如下技术措施来实现：基于投影菲涅尔带的三维起伏地表高斯束演模拟方法，该基于投影菲涅尔带的三维起伏地表高斯束演模拟方法包括：步骤1，建立起伏地表的三维地质模型；步骤2，进行 运动学射线追踪；步骤3，进行动力学射线追踪；步骤4，在三维介质中用投影菲涅尔带椭圆来约束高斯束在射线终点处能量的分布；以及步骤5，采用高斯波包法合成地震记录。

本发明的目的还可通过如下技术措施来实现：

在步骤2中，运动学射线追踪在三维坐标系下满足如下方程：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{d\τ}}=v\sin \mathrm{\α}*\cos \mathrm{\β}\\ \frac{\mathrm{dy}}{\mathrm{d\τ}}=v\sin \mathrm{\α}*\sin \mathrm{\β}\\ \frac{\mathrm{dz}}{\mathrm{d\τ}}=v\cos \mathrm{\α}\\ \frac{\mathrm{d\α}}{\mathrm{d\τ}}=-\cos \mathrm{\α}\left(v_x\cos \mathrm{\β}{+v}_y\sin \mathrm{\β}\right)+v_z\sin \mathrm{\α}\\ \frac{\mathrm{d\β}}{\mathrm{d\τ}}=\frac{1}{\sin \mathrm{\α}}(v_x\sin \mathrm{\β}-v_y\cos \mathrm{\β})\end{array}\right.---\left(1\right)$

其中：τ为旅行时，v为中心射线的P波速度，v_i(i＝x,y,z)为中心射线处速度函数的偏导数，α和β分别为其坐标系下的倾角与方位角；在层内采用四阶龙格-库塔法进行求解各时间步的坐标，在界面上利用二分法求出交点的近似坐标，并利用矢量Snell定律计算反射或透射方向。

在步骤3中，将动力学追踪方程：

$\frac{\mathrm{dM}\left(s\right)}{\mathrm{ds}}+v\left(s\right)M^2\left(s\right)+v^{-2}\left(s\right)V=0---\left(2\right)$

其中： $V=\left[\begin{array}{cc}\frac{{\∂}^{2}v\left(s\right)}{{\∂n}^2}& \frac{{\∂}^{2}v\left(s\right)}{\∂n\∂m}\\ \frac{{\∂}^{2}v\left(s\right)}{\∂m\∂n}& \frac{{\∂}^{2}v\left(s\right)}{{\∂m}^2}\end{array}\right],M={\mathrm{PQ}}^{-1}$

化为线性动力学追踪方程：

dX_i/ds＝HX_i (i＝1,2) (3)

其中： $X=\left(\begin{array}{c}Q\\ P\end{array}\right),$ X_i为X的列向量， $H=\left(\begin{array}{cc}0& \mathrm{vI}\\ {-v}^2\left(s\right)V& 0\end{array}\right),$ 0为2×2零矩阵，I为2×2单位矩阵；(s,m,n)为三维中心射线坐标，V为速度二阶偏导数矩 阵，Q,P为动力射线追踪方程的参数；

在给定四阶单位矩阵的初始条件下，得到X的通解：

X＝Π(s)·C (4)

其中：Π(s)为4×4传播矩阵，其列向量为方程(3)四组线性独立的解，C为一个2×4复值初始参数矩阵；

把Π(s)和C分块如下：

$\left\{\begin{array}{c}\mathrm{\Π}\left(s\right)=\left(\begin{array}{cc}{\mathrm{\Π}}_{11}& {\mathrm{\Π}}_{12}\\ {\mathrm{\Π}}_{21}& {\mathrm{\Π}}_{22}\end{array}\right)\\ C=\left(\begin{array}{c}{C}_1\\ C_2\end{array}\right)\end{array}\right.---\left(5\right)$

其中：Π_ij和C_i(i,j＝1,2)为2×2子矩阵；

这样，根据 $X=\left(\begin{array}{c}Q\\ P\end{array}\right),$ 得到高斯束的动力学参数矩阵：

$\left\{\begin{array}{c}Q={\mathrm{\Π}}_{11}{C}_1+{\mathrm{\Π}}_{12}{C}_2\\ P={\mathrm{\Π}}_{21}{C}_1+{\mathrm{\Π}}_{22}{C}_2\\ M={\mathrm{PQ}}^{-1}=({\mathrm{\Π}}_{21}{C}_1+{\mathrm{\Π}}_{22}{C}_2){({\mathrm{\Π}}_{11}{C}_1+{\mathrm{\Π}}_{12}{C}_2)}^{-1}\end{array}\right.---\left(6\right)$

由于射线经过界面时传播矩阵Π(s)发生突变，需要在界面处重新计算Π(s)；

在射线与界面交点O点处边界条件为：

$\mathrm{\Π}(\tilde{O},S_0)=\mathrm{\Π}(\tilde{O},O)\mathrm{\Π}(O,S_0)---\left(7\right)$

$\mathrm{\Π}(\tilde{O},O)=\left[\begin{array}{cc}{G}^T\left(\tilde{O}\right)G^{-T}\left(O\right)& 0\\ G^{-1}\left(\tilde{O}\right)[E\left(O\right)-E\left(\tilde{O}\right)-\mathrm{\μD}]G^{-T}\left(O\right)& G^{-1}\left(\tilde{O}\right)G\left(O\right)\end{array}\right]---\left(8\right)$

$G\left(O\right)=\left[\begin{array}{cc}\mathrm{\ϵ}{\cos i}_s\cos \mathrm{\κ}& -\mathrm{\ϵ}{\cos i}_s\sin \mathrm{\κ}\\ \sin \mathrm{\κ}& \cos \mathrm{\κ}\end{array}\right]G\left(\tilde{O}\right)=\left[\begin{array}{cc}\±\mathrm{\ϵ}{\cos i}_R\cos \mathrm{\κ}& \stackrel{-}{+}\mathrm{\ϵ}{\cos i}_s\sin \mathrm{\κ}\\ \sin \mathrm{\κ}& \cos \mathrm{\κ}\end{array}\right]$

$E=\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& 0\end{array}\right]$

E₁₁(O)＝-sini_sv^-2(O)[(1+cos²i_s)v_,z1-εcosi_ssini_sv_,z3]

E₁₂(O)＝E₂₁(O)＝-sini_sv^-2(O)v_,z2 (9)

$E_{11}\left(\tilde{O}\right)=-\sin i_R{v}^{-2}\left(\tilde{O}\right)\left[(1+{\cos }^2{i}_R){\tilde{v}}_{,z1}\stackrel{-}{+}\mathrm{\ϵ}\cos i_R\sin i_R{\tilde{v}}_{,z3}\right]$

$E_{12}\left(\tilde{O}\right)=E_{21}\left(\tilde{O}\right)=-\sin i_R{v}^{-2}\left(\tilde{O}\right){\tilde{v}}_{,z2}$

$\mathrm{\μ}=\mathrm{\ϵ}\left(v^{-1}\left(O\right)\cos i_s\stackrel{-}{+}{v}^{-1}\left(\tilde{O}\right)\cos i_R\right)$

其中：S₀为炮点，O为射线与界面交点，为下一个界面交点，G^T为G的转置，为方向指数；i_s为入射角，i_R为反射角或透射角；D为界面的曲率矩阵；z1、z2和z3为局部笛卡尔坐标系中的三个分量；κ为射线中心坐标系中e_n和局部笛卡尔坐标系中的夹角；当发生透射时取等式(9)上面的符号，发生反射时取等式(9)下面的符号。

征在于，在步骤4中，在三维介质中用投影菲涅尔带椭圆来约束高斯束在射线终点处能量的分布；

当中心射线终点处高斯束的有效半宽度椭圆与射线的投影菲涅尔带椭圆一致时，可得：

$\mathrm{eig}\left(\mathrm{Im}\left(M\left(G_r\right)\right)\right)=\mathrm{eig}\left(\frac{1}{\mathrm{\π}}\mathrm{Hp}\left(G_r\right)\right)---\left(10\right)$

其中：eig为特征值，G_r为中心射线终点，Hp为投影菲涅尔带矩阵，其值可以通过经典射线传播矩阵Π(s)、面面传播矩阵Τ(s)以及两者转换关系求出；

根据传播矩阵Π(s)的辛属性及数学推导可得到：

$\left\{\begin{array}{c}M\left(s_0\right)=\left[\begin{array}{cc}{\mathrm{i\ϵ}}_1& 0\\ 0& {\mathrm{i\ϵ}}_2\end{array}\right]\\ {\mathrm{\ϵ}}_1=\frac{\mathrm{\π}/{\mathrm{\ζ}}_1+\sqrt{{(\mathrm{\π}/{\mathrm{\ζ}}_1)}^2-4{{\mathrm{\λ}}_1}^2{{\mathrm{\η}}_1}^2}}{{{2\mathrm{\η}}_1}^2}\\ {\mathrm{\ϵ}}_2=\frac{\mathrm{\π}/{\mathrm{\ζ}}_2+\sqrt{{(\mathrm{\π}/{\mathrm{\ζ}}_2)}^2-{{4\mathrm{\λ}}_2}^2{{\mathrm{\η}}_2}^2}}{{{2\mathrm{\η}}_2}^2}\end{array}\right.---\left(11\right)$

其中：λ₁、λ₂(λ₁≥λ₂)为Π₁₁的特征值，η₁、η₂(η₁≥η₂)为Π₁₂的特征值，ζ₁、ζ₂(ζ₁≥ξ₂)为投影菲涅尔带矩阵Hp的特征值。

在步骤5中，记以初始角(α_i,β_j)出发在检波点R处t时刻高斯波包为g(R,t,α_i,β_j)，其近似表达式为：

$\begin{array}{c}g(R,t,{\mathrm{\α}}_i,{\mathrm{\β}}_j)={2\mathrm{\πf}}_m\left|\mathrm{\ΦA}\right|\exp \{-{\left[{2\mathrm{\πf}}_m(t-\mathrm{\θ})\mathrm{\γ}\right]}^2+{({2\mathrm{\πf}}_{m}G/\mathrm{\γ})}^2-{2\mathrm{\πf}}_{m}G\}\\ \×\cos ({2\mathrm{\πf}}^{*}(t-\mathrm{\θ})+v-\arg \left(\mathrm{\ΦA}\right)+\mathrm{\π}/2)\end{array}---\left(12\right)$

其中：f_m、γ、ν为高斯子波参数，R'为R在中心射线上的投影， q^T＝(q₁,q₂),q₁,q₂为射线在R点的坐标，A为振幅，在层状介质中 $A\left(R\right)=A_0{\left[\frac{\mathrm{\ρ}\left(s_0\right)v\left(s_0\right)\·\det \left(Q\left(s_0\right)\right)}{\mathrm{\ρ}\left(R\right)v\left(R\right)\·\det \left(Q\left(R\right)\right)}\right]}^{1/2}\·\underset{i=1}{\overset{N}{\mathrm{\Π}}}{R}_i\·\underset{i=1}{\overset{N}{\mathrm{\Π}}}\left[\frac{\widetilde{\mathrm{\ρ}}\tilde{v}}{\mathrm{\ρv}}\right]\·\underset{i=1}{\overset{N}{\mathrm{\Π}}}\left(\frac{\det \left(\tilde{Q}\right)}{\det \left(Q\right)}\right)$ ；ρ为密度，v为速度，R_i为经过第i个界面反射系数，N表示射线经过界面的个数，“～”表示生成射线一侧的量值；Φ＝ω/2π·|det(Q(R'))|^1/2·{-det[M(R')-Re(M(R'))]}，为能量叠加的权系数；

那么，检波点R处离散的能量叠加表达式为：

$u(R,t)=\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}g(R,t,{\mathrm{\α}}_i,{\mathrm{\β}}_j)\mathrm{\Δ\α\Δ\β}---\left(13\right)$

其中：I，J分别表示射线在倾角和方位角上离散的个数，Δα，Δβ为射线间隔。

在步骤1中，结合地震和测井资料及其解释结果包括层位、断层等， 建立模型框架，之后加上各层的速度参数，建立起伏地表的三维地质模型。

本发明中的基于投影菲涅尔带的三维起伏地表高斯束演模拟方法，引入了三维投影菲涅尔带的思想，发展了基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，在中深部地层具有更好的保幅性，对起伏地表和复杂模型具有更好的适应性。相对于常规方法，在投影法菲涅尔带约束下高斯束的传播更稳定，能量聚集性更好；本方法使高斯束能量分布在投影菲涅尔带内，赋予了能量分布范围明确的物理意义，使得基于渐进射线理论的高斯束更符合波动理论；与常规初始参数正演结果相比，本方法在中深部地层具有更好的保幅性，对起伏地表和复杂模型具有更好的适应性。

附图说明

图1为本发明的基于投影菲涅尔带的三维起伏地表高斯束演模拟方法的一具体实施例的流程图；

图2为本发明的一具体实施例中用于测试的简单三维水平地表层状模型；

图3为本发明的一具体实施例中利用两种初始参数对三维水平地表层状模型进行三维高斯束正演模拟得到的结果的示意图；

图4为本发明的一具体实施例中用于测试的经典起伏地表台地模型；

图5为本发明的一具体实施例中起伏地表台地模型盲区示意图；

图6为本发明的一具体实施例中利用两种初始参数对起伏地表台地模型进行三维高斯束正演模拟得到的结果的示意图；

图7为本发明的一具体实施例中用于测试的起伏地表复杂模型；

图8为本发明的一具体实施例中利用两种初始参数对起伏地表复杂模型进行三维高斯束正演模拟得到的结果的示意图。

具体实施方式

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

如图1所示，图1为本发明的基于投影菲涅尔带的三维起伏地表高斯 束演模拟方法的流程图。

在步骤101，建立起伏地表的三维地质模型。结合地震和测井资料及其解释结果包括层位、断层等，建立模型框架，之后加上各层的速度参数，建立起伏地表的三维地质模型。流程进入到步骤102。

在步骤102，进行运动学射线追踪。

运动学射线追踪在三维坐标系下满足如下方程：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{d\τ}}=v\sin \mathrm{\α}*\cos \mathrm{\β}\\ \frac{\mathrm{dy}}{\mathrm{d\τ}}=v\sin \mathrm{\α}*\sin \mathrm{\β}\\ \frac{\mathrm{dz}}{\mathrm{d\τ}}=v\cos \mathrm{\α}\\ \frac{\mathrm{d\α}}{\mathrm{d\τ}}=-\cos \mathrm{\α}(v_x\cos \mathrm{\β}+v_y\sin \mathrm{\β})+v_z\sin \mathrm{\α}\\ \frac{\mathrm{d\β}}{\mathrm{d\τ}}=\frac{1}{\sin \mathrm{\α}}(v_x\sin \mathrm{\β}-v_y\cos \mathrm{\β})\end{array}\right.---\left(1\right)$

其中：τ为旅行时，v为中心射线的P波速度，v_i(i＝x,y,z)为中心射线处速度函数的偏导数，α和β分别为其坐标系下的倾角与方位角；在层内采用四阶龙格-库塔法进行求解各时间步的坐标，在界面上利用二分法求出交点的近似坐标，并利用矢量Snell定律计算反射或透射方向。流程进入到步骤103。

在步骤103，进行动力学射线追踪。

将动力学追踪方程：

$\frac{\mathrm{dM}\left(s\right)}{\mathrm{ds}}+v\left(s\right)M^2\left(s\right)+v^{-2}\left(s\right)V=0---\left(2\right)$

其中： $V=\left[\begin{array}{cc}\frac{{\∂}^{2}v\left(s\right)}{{\∂n}^2}& \frac{{\∂}^{2}v\left(s\right)}{\∂n\∂m}\\ \frac{{\∂}^{2}v\left(s\right)}{\∂m\∂n}& \frac{{\∂}^{2}v\left(s\right)}{{\∂m}^2}\end{array}\right],M={\mathrm{PQ}}^{-1}$

化为线性动力学追踪方程：

dX_i/ds＝HX_i (i＝1,2) (3)

其中： $X=\left(\begin{array}{c}Q\\ P\end{array}\right),$ X_i为X的列向量， $H=\left(\begin{array}{cc}0& \mathrm{vI}\\ {-v}^2\left(s\right)V& 0\end{array}\right),$ 0为2×2零矩阵，I为2×2单位矩阵；(s,m,n)为三维中心射线坐标，V为速度二阶偏导数矩阵，Q,P为动力射线追踪方程的参数；

在给定四阶单位矩阵的初始条件下，得到X的通解：

X＝Π(s)·C (4)

其中：Π(s)为4×4传播矩阵，其列向量为方程(3)四组线性独立的解，C为一个2×4复值初始参数矩阵；

把Π(s)和C分块如下：

$\left\{\begin{array}{c}\mathrm{\Π}\left(s\right)=\left(\begin{array}{cc}{\mathrm{\Π}}_{11}& {\mathrm{\Π}}_{12}\\ {\mathrm{\Π}}_{21}& {\mathrm{\Π}}_{22}\end{array}\right)\\ C=\left(\begin{array}{c}{C}_1\\ C_2\end{array}\right)\end{array}\right.---\left(5\right)$

其中：Π_ij和C_i(i,j＝1,2)为2×2子矩阵；

这样，根据 $X=\left(\begin{array}{c}Q\\ P\end{array}\right),$ 得到高斯束的动力学参数矩阵：

$\left\{\begin{array}{c}Q={\mathrm{\Π}}_{11}{C}_1+{\mathrm{\Π}}_{12}{C}_2\\ P={\mathrm{\Π}}_{21}{C}_1+{\mathrm{\Π}}_{22}{C}_2\\ M={\mathrm{PQ}}^{-1}=({\mathrm{\Π}}_{21}{C}_1+{\mathrm{\Π}}_{22}{C}_2){({\mathrm{\Π}}_{11}{C}_1+{\mathrm{\Π}}_{12}{C}_2)}^{-1}\end{array}\right.---\left(6\right)$

由于射线经过界面时传播矩阵Π(s)发生突变，需要在界面处重新计算Π(s)；

在射线与界面交点O点处边界条件为：

$\mathrm{\Π}(\tilde{O},S_0)=\mathrm{\Π}(\tilde{O},O)\mathrm{\Π}(O,S_0)---\left(7\right)$

$\mathrm{\Π}(\tilde{O},O)=\left[\begin{array}{cc}{G}^T\left(\tilde{O}\right)G^{-T}\left(O\right)& 0\\ G^{-1}\left(\tilde{O}\right)[E\left(O\right)-E\left(\tilde{O}\right)-\mathrm{\μD}]G^{-T}\left(O\right)& G^{-1}\left(\tilde{O}\right)G\left(O\right)\end{array}\right]---\left(8\right)$

$G\left(O\right)=\left[\begin{array}{cc}\mathrm{\ϵ}{\cos i}_s\cos \mathrm{\κ}& -\mathrm{\ϵ}{\cos i}_s\sin \mathrm{\κ}\\ \sin \mathrm{\κ}& \cos \mathrm{\κ}\end{array}\right]G\left(\tilde{O}\right)=\left[\begin{array}{cc}\±\mathrm{\ϵ}{\cos i}_R\cos \mathrm{\κ}& \stackrel{-}{+}\mathrm{\ϵ}{\cos i}_s\sin \mathrm{\κ}\\ \sin \mathrm{\κ}& \cos \mathrm{\κ}\end{array}\right]$

$E=\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& 0\end{array}\right]$

E₁₁(O)＝-sini_sv^-2(O)[(1+cos²i_s)v_,z1-εcosi_ssini_sv_,z3]

E₁₂(O)＝E₂₁(O)＝-sini_sv^-2(O)v_,z2 (9)

$E_{11}\left(\tilde{O}\right)=-\sin i_R{v}^{-2}\left(\tilde{O}\right)[(1+{\cos }^2{i}_R)\tilde{v},z1\stackrel{-}{+}\mathrm{\ϵ}\cos i_R\sin i_R\tilde{v},z3]$

$E_{12}\left(\tilde{O}\right)=E_{21}\left(\tilde{O}\right)=-\sin i_R{v}^{-2}\left(\tilde{O}\right){\tilde{v}}_{,z2}$

$\mathrm{\μ}=\mathrm{\ϵ}\left(v^{-1}\left(O\right)\cos i_s\stackrel{-}{+}{v}^{-1}\left(\tilde{O}\right)\cos i_R\right)$

其中：S₀为炮点，O为射线与界面交点，为下一个界面交点，G^T为G的转置，为方向指数；i_s为入射角，i_R为反射角或透射角；D为界面的曲率矩阵；z1、z2和z3为局部笛卡尔坐标系中的三个分量；κ为射线中心坐标系中e_n和局部笛卡尔坐标系中的夹角；

当发生透射时取等式(9)上面的符号，发生反射时取等式(9)下面的符号。流程进入到步骤104。

在步骤104，在三维介质中用投影菲涅尔带椭圆来约束高斯束在射线终点处能量的分布。

当中心射线终点处高斯束的有效半宽度椭圆与射线的投影菲涅尔带椭圆一致时，可得：

$\mathrm{eig}\left(\mathrm{Im}\left(M\left(G_r\right)\right)\right)=\mathrm{eig}\left(\frac{1}{\mathrm{\π}}\mathrm{Hp}\left(G_r\right)\right)---\left(10\right)$

其中：eig为特征值，G为中心射线终点，Hp为投影菲涅尔带矩阵，其值可以通过经典射线传播矩阵Π(s)、面面传播矩阵Τ(s)以及两者转换关系求出。

根据传播矩阵Π(s)的辛属性及数学推导可得到：

$\left\{\begin{array}{c}M\left(s_0\right)=\left[\begin{array}{cc}{\mathrm{i\ϵ}}_1& 0\\ 0& {\mathrm{i\ϵ}}_2\end{array}\right]\\ {\mathrm{\ϵ}}_1=\frac{\mathrm{\π}/{\mathrm{\ζ}}_1+\sqrt{{(\mathrm{\π}/{\mathrm{\ζ}}_1)}^2-4{{\mathrm{\λ}}_1}^2{{\mathrm{\η}}_1}^2}}{{{2\mathrm{\η}}_1}^2}\\ {\mathrm{\ϵ}}_2=\frac{\mathrm{\π}/{\mathrm{\ζ}}_2+\sqrt{{(\mathrm{\π}/{\mathrm{\ζ}}_2)}^2-{{4\mathrm{\λ}}_2}^2{{\mathrm{\η}}_2}^2}}{{{2\mathrm{\η}}_2}^2}\end{array}\right.---\left(11\right)$

其中：λ₁、λ₂(λ₁≥λ₂)为Π₁₁的特征值，η₁、η₂(η₁≥η₂)为Π₁₂的特征值，ζ₁、ζ₂(ζ₁≥ξ₂)为投影菲涅尔带矩阵Hp的特征值。流程进入到步骤105。

在步骤105，高斯波包法合成地震记录。

记以初始角(α_i,β_j)出发在检波点R处t时刻高斯波包为g(R,t,α_i,β_j)，其近似表达式为：

$\begin{array}{c}g(R,t,{\mathrm{\α}}_i,{\mathrm{\β}}_j)={2\mathrm{\πf}}_m\left|\mathrm{\ΦA}\right|\exp \{-{\left[{2\mathrm{\πf}}_m(t-\mathrm{\θ})\mathrm{\γ}\right]}^2+{({2\mathrm{\πf}}_{m}G/\mathrm{\γ})}^2-{2\mathrm{\πf}}_{m}G\}\\ \×\cos ({2\mathrm{\πf}}^{*}(t-\mathrm{\θ})+v-\arg \left(\mathrm{\ΦA}\right)+\mathrm{\π}/2)\end{array}---\left(12\right)$

其中：f_m、γ、ν为高斯子波参数，R'为R在中心射线上的投影， q^T＝(q₁,q₂),q₁,q₂为射线在R点的坐标，A为振幅，在层状介质中 $A\left(R\right)=A_0{\left[\frac{\mathrm{\ρ}\left(s_0\right)v\left(s_0\right)\·\det \left(Q\left(s_0\right)\right)}{\mathrm{\ρ}\left(R\right)v\left(R\right)\·\det \left(Q\left(R\right)\right)}\right]}^{1/2}\·\underset{i=1}{\overset{N}{\mathrm{\Π}}}{R}_i\·\underset{i=1}{\overset{N}{\mathrm{\Π}}}\left[\frac{\widetilde{\mathrm{\ρ}}\tilde{v}}{\mathrm{\ρv}}\right]\·\underset{i=1}{\overset{N}{\mathrm{\Π}}}\left(\frac{\det \left(\tilde{Q}\right)}{\det \left(Q\right)}\right)$ ；ρ为密度，v为速度，R_i为经过第i个界面反射系数，N表示射线经过界面的个数，“～”表示生成射线一侧的量值；Φ＝ω/2π·|det(Q(R'))|^1/2·{-det[M(R')-Re(M(R'))]}，为能量叠加的权系数；

那么，检波点R处离散的能量叠加表达式为：

$u(R,t)=\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}g(R,t,{\mathrm{\α}}_i,{\mathrm{\β}}_j)\mathrm{\Δ\α\Δ\β}---\left(13\right)$

其中：_I，_J分别表示射线在倾角和方位角上离散的个数，Δα，Δβ为射线间隔。

在应用本发明的一具体实施例中，包括：

1)选取如图2所示的简单三维水平地表层状模型，模型尺寸为4000m×4000m×4000m，各层速度如图2所示，Δα＝2°，Δβ＝10°；

2)观测系统为：三炮四线制(如图2，实线为炮线，虚线为检波点线)，道间距25m，线间距500m，炮间距50m，炮线距500m，采样间隔4ms，记录时间4s，全排列接收。经试算得到的射线路径如图1中红色线条所示；

3)利用两种初始参数进行三维高斯束正演模拟得到的单炮记录如图3所示；

在应用本发明的另一具体实施例中，包括：

1)选取如图4所示的经典起伏地表台地模型，模型尺寸为4000m×4000m×2500m，上层速度为2000m/s，下层速度为3000m/s，Δα＝2°，Δβ＝2°；

2)观测系统与图2所示模型相同，试算得到的射线路径如图4所示，得到射线追踪的盲区如图5；

3)利用两种初始参数进行三维高斯束正演模拟得到的单炮记录如图6所示；

在应用本发明的一具体实施例中，包括：

1)选取如图6所示的经典起伏地表台地模型，模型尺寸为4000m×4000m×4500m，接收时间5s，速度分布如图7；

2)观测系统与图2所示模型相同，试算得到的射线路径如图7所示；

3)利用两种初始参数进行三维高斯束正演模拟得到的单炮记录如图8所示；

对比图3中(a)、(b)两者炮记录，(a)为常规初始参数模拟的炮记录，(b)为投影菲涅尔带约束下模拟的炮记录，可知，(1)两者同相轴的形态和位置是一致的，说明两者的运动学特征一致；(2)在纵向上，a和b同相轴主要能量分布范围(图3中矩形框)从浅层到深层都逐渐增大，但b比a增加的更慢些，这是因为投影菲涅尔带椭圆和常规有效半宽度椭圆都随射线路径增大而增大，但前者增大的速度没有后者速度快；(3)在横向上， 同一深度b中同相轴主要能量分布范围比a小，并且深度越大，两者差异越大；经过此简单模型的试算和分析，验证了基于投影菲涅尔带的高斯束正演方法的正确性，尤其是本方法还具有一定的保幅性。

图6(a)为常规初始参数模拟的炮记录，(b)为投影菲涅尔带约束下模拟的炮记录，由图5和图6可知，(1)对比两种初始参数合成的地震记录，可知本方法对起伏地表角点模型是适应的；(2)在普通射线方法盲区内(如图5中的空白区域)，两种初始参数试算的炮记录都有一定的绕射能量(图6中箭头)，但b中绕射能量范围(约25道)比a中(约40道)的小，并且没有a中的能量干涉现象(图6a矩形框)，究其原因可能是在焦点处投影菲涅尔带范围小于常规初始参数计算的有效范围，这样使得绕射能量分布在合理的范围内。经过起伏地表台地模型的计算结果证实了本方法的对起伏地表和焦点绕射一定的适应性和保幅性。

图8(a)为常规初始参数模拟的炮记录，(b)为投影菲涅尔带约束下模拟的炮记录，对比图8中(a)、(b)两者炮记录可知，本方法对复杂构造也具有较好的适应性，对绕射能量具有很好的控制作用，因而模拟结果较常规三维高斯束正演模拟方法具有更好的保幅性。

Claims (6)

1.基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，其特征在于，该基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法包括：

步骤1，建立起伏地表的三维地质模型；

步骤2，进行运动学射线追踪；

步骤3，进行动力学射线追踪；

步骤4，在三维介质中用投影菲涅尔带椭圆来约束高斯束在射线终点处能量的分布；以及

步骤5，采用高斯波包法合成地震记录。

2.根据权利要求1所述的基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，其特征在于，在步骤2中，运动学射线追踪在三维坐标系下满足如下方程：

其中：τ为旅行时，v为中心射线的P波速度，v_i(i＝x,y,z)为中心射线处速度函数的偏导数，α和β分别为其坐标系下的倾角与方位角；在层内采用四阶龙格-库塔法进行求解各时间步的坐标，在界面上利用二分法求出交点的近似坐标，并利用矢量Snell定律计算反射或透射方向。

3.根据权利要求2所述的基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，其特征在于，在步骤3中，将动力学追踪方程：

其中：M＝PQ^-1

化为线性动力学追踪方程：

dX_i/ds＝HX_i(i＝1,2) (3)

其中：X_i为X的列向量，0为2×2零矩阵，I为2×2单位矩阵；(s,m,n)为三维中心射线坐标，V为速度二阶偏导数矩阵，Q,P为动力射线追踪方程的参数；

在给定四阶单位矩阵的初始条件下，得到X的通解：

X＝Π(s)·C (4)

其中：Π(s)为4×4传播矩阵，其列向量为方程(3)四组线性独立的解，C为一个2×4复值初始参数矩阵；

把Π(s)和C分块如下：

其中：Π_ij和C_i(i,j＝1,2)为2×2子矩阵；

这样，根据得到高斯束的动力学参数矩阵：

由于射线经过界面时传播矩阵Π(s)发生突变，需要在界面处重新计算Π(s)；

在射线与界面交点O点处边界条件为：

E₁₁(O)＝-sini_sv^-2(O)[(1+cos²i_s)v_,z1-εcosi_ssini_sv_,z3]

E₁₂(O)＝E₂₁(O)＝-sini_sv^-2(O)v_,z2

其中：S₀为炮点，O为射线与界面交点，为下一个界面交点，G^T为G的转置，为方向指数；i_s为入射角，i_R为反射角或透射角；D为界面的曲率矩阵；z1、z2和z3为局部笛卡尔坐标系中的三个分量；κ为射线中心坐标系中e_n和局部笛卡尔坐标系中的夹角；当发生透射时取等式(9)上面的符号，发生反射时取等式(9)下面的符号。

4.根据权利要求3所述的基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，其特征在于，在步骤4中，在三维介质中用投影菲涅尔带椭圆来约束高斯束在射线终点处能量的分布；

当中心射线终点处高斯束的有效半宽度椭圆与射线的投影菲涅尔带椭圆一致时，可得：

其中：eig为特征值，G_r为中心射线终点，Hp为投影菲涅尔带矩阵，其值可以通过经典射线传播矩阵Π(s)、面面传播矩阵Τ(s)以及两者转换关系求出；

根据传播矩阵Π(s)的辛属性及数学推导可得到：

其中：λ₁、λ₂(λ₁≥λ₂)为Π₁₁的特征值，η₁、η₂(η₁≥η₂)为Π₁₂的特征值，ζ₁、ζ₂(ζ₁≥ξ₂)为投影菲涅尔带矩阵Hp的特征值。

5.根据权利要求4所述的基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，其特征在于，在步骤5中，记以初始角(α_i,β_j)出发在检波点R处t时刻高斯波包为g(R,t,α_i,β_j)，其近似表达式为：

g(R,t,α_i,β_j)＝2πf_m|ΦA|exp{-[2πf_m(t-θ)γ]²+(2πf_mG/γ)²-2πf_mG}

(12)

×cos(2πf^*(t-θ)+ν-arg(ΦA)+π/2)

其中：f_m、γ、ν为高斯子波参数，R'为R在中心射线上的投影， q^T＝(q₁,q₂),q₁,q₂为射线在R点的坐标，A为振幅，在层状介质中

ρ为密度，v为速度，R_i为经过第i个界面反射系数，N表示射线经过界面的个数，“～”表示生成射线一侧的量值；Φ＝ω/2π·|det(Q(R'))|^1/2·{-det[M(R')-Re(M(R'))]}，为能量叠加的权系数；

那么，检波点R处离散的能量叠加表达式为：

其中：I，J分别表示射线在倾角和方位角上离散的个数，Δα，Δβ为射线间隔。

6.根据权利要求1所述的基于投影菲涅尔带的三维起伏地表高斯束正演模拟方法，其特征在于，在步骤1中，结合地震和测井资料及其解释结果包括层位、断层等，建立模型框架，之后加上各层的速度参数，建立起伏地表的三维地质模型。