一种高效率多分量地震资料真振幅偏移成像方法，应用于多分量地震数据处理技术领域，该方法是在弹性波逆时偏移方法的基础上进行了改进，特征是直接以多分量地震资料为输入，通过编码算法、分段成像等措施，再在反演的框架下，实现基于反演的弹性波逆时偏移。本发明将编码、时间分段及反演的思想引入弹性波逆时偏移中，与常规弹性波逆时偏移相比，本发明可以获得高精度且振幅保真的叠前深度偏移剖面，同时算法的效率很高。本发明完整地保持纵横波的矢量特性、振幅以及相位特性，偏移剖面具有保真性；有效地消除了纵横波之间串扰造成的偏移假象，提高了成像的精度；有效降低了弹性波逆时偏移的存储量，提高了计算效率。
1.一种高效率多分量地震资料真振幅偏移成像方法，其特征在于包括以下步骤：
(a)：通过野外观测获得的多炮多分量地震记录D(xr,t；xs)＝(Dx,Dy,Dz)，给定背景模型、预设参数，确定进行偏移成像的多分量地震记录、观测系统参数及震源函数S(x,t；xs)，其中，x表示空间位置矢量，xs表示震源位置矢量，xr表示检波点位置矢量，t表示时间变量；
(b)：基于读取的观测系统参数，设计平面波编码函数C(xs,p)，具体为：
其中，p为射线参数变量，fδ表示脉冲函数；此处平面波编码函数C(xs,p)等价于时移函数，时移量p·xs随震源点位置xs不同而改变，射线参数p具体为：其中，α为地表纵波速度，θ为入射角；通过不同的入射角θ，可以获得不同的射线参数p；
通过所述编码函数，并利用每一个射线参数p对观测的多炮多分量地震记录D(xr,t；xs)进行编码，获得关于每一个射线参数p的编码多分量炮集DC(xr,t；p)，具体为：DC(xr,t；p)＝∫D(xr,t；xs)·C(xs,p)dxs； (3)同样，采用相同的编码函数，并利用每一个射线参数p对步骤(a)中所述观测系统中的震源函数S(x,t；xs)进行编码，获得每个射线参数p的编码震源函数SC(x,t；p)，具体为：SC(x,t；p)＝∫S(x,t；xs)·C(xs,p)dxs； (4)(c)：利用背景模型及观测的多炮多分量地震记录的最大记录时间tmax，设计记录时间t的分段区间，具体为：将记录时间t分为三段，用下式表示：其中，tR1，tR2和tR3分别表示记录时间的三个区间范围，t1(x)和t2(x)的具体计算式为：其中，αmax为背景模型的最大速度，hz代表深度，hx代表水平距离；
(d)：针对每一个射线参数p，均执行以下操作：利用每一个射线参数p，通过步骤(b)中的编码方式获得与该射线参数对应的编码震源函数，基于常规的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的每个时刻的震源波场，在该震源波场延拓的每一个时间步中，利用步骤(c)中的记录时间分段区间作为时间约束条件，利用常规的弹性波场能量计算方程计算震源波场每个时刻、每个位置的波场能量值，通过能量数值大小判断，获得每个时间区间内、每个位置波场能量最大值所对应的时刻，将此时刻作为该时间区间内、该位置所对应的成像时间T(x)，在该时刻、该位置处，利用成像条件公式，保存成像条件中与震源波场相关的项，同时保存该时刻、该位置的波场最大能量值；
进而获得每一个射线参数对应的成像时间、成像条件中与震源波场相关的项及波场最大能量值；最终将每一个射线参数对应的每一个时间段内的波场最大能量值进行叠加，获得照明度补偿因子；
(e)：针对每一个射线参数p，均执行以下操作：对每一个射线参数p，同样基于所述的各向同性介质的弹性波方程，对该射线参数p对应的步骤(b)所述编码多分量炮集DC(xr,t；p)进行波场逆时延拓，获得该射线参数对应的每一个时刻的检波点波场，在该检波点波场延拓的每一个时间步中，利用步骤(d)中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤(d)中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的偏移剖面；
进而获得每一个射线参数对应的偏移剖面；将所有的射线参数对应的偏移剖面进行叠加，同时利用步骤(d)中照明度补偿因子对叠加的偏移剖面进行能量补偿，获得初始偏移剖面，也即第0次迭代的偏移剖面(f)：针对每一个射线参数p，执行以下操作：对每一个射线参数p，通过步骤(b)中的编码方式获得与该射线参数对应的编码震源函数，同样基于所述的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的背景波场；利用步骤(e)中所述的初始偏移剖面I0、所述的该射线参数的背景波场以及步骤(a)中给定的背景模型，构建该射线参数对应的虚拟震源；基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该射线参数对应的编码震源函数进行波场正向延拓，获得该射线参数对应的反偏移波场，对该射线参数对应的反偏移波场进行记录采样，获得该射线参数对应的编码多分量反偏移炮集；
进而获得每一个射线参数对应的编码多分量反偏移炮集，所有炮的编码多分量炮集构成了初始编码多分量反偏移炮集，也即第0次迭代的编码多分量反偏移炮集此时，设置当前的迭代次数为i＝0；
(g)：设置当前的迭代次数i＝i+1，基于第i-1次迭代所得的所述偏移剖面
和 第 i - 1 次 迭 代 所 得 的 所 述 编 码 多 分 量 反 偏 移 炮 集进行反演迭代；
在第i次迭代中，针对每一个射线参数p，将第i-1次迭代所得的该射线参数对应的编码多分量反偏移炮集di-1与该射线参数对应的观测的编码多分量炮集DC作差，计算该射线参数对应的编码多分量炮集残差；进而获得所有射线参数对应的编码多分量炮集残差所述的将第i-1次迭代所得的该射线参数对应的编码多分量反偏移炮集di-1与该射线参数对应的观测的编码多分量炮集DC作差，计算该射线参数对应的编码多分量炮集残差，具体为：Δdi-1＝di-1-DC； (7)
(h)：针对每一个射线参数p，执行以下操作：对每一个射线参数p，同样基于所述的各向同 性 介质的 弹性波 方程 ，对该 射线 参数对 应的 编码多 分量炮 集残 差进行波场逆时延拓，获得该射线参数对应的每一个时刻的
检波点波场，在该检波点波场延拓的每一个时间步中，利用步骤(d)中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤(d)中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的梯度剖面；
进而获得每一个射线参数对应的梯度剖面；将所有的射线参数对应的梯度剖面进行叠加，同时利用步骤(d)中照明度补偿因子对叠加的梯度剖面进行能量补偿，获得本次迭代的梯度剖面，也即第i次迭代的梯度剖面(i)：利用最优化反演算法，构建第i次迭代的下降方向剖面
(j)：针对每一个射线参数p，通过步骤(b)中的编码方式获得与该射线参数对应的编码震源函数，同样基于所述的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的背景波场；利用步骤(i)中所获得的第i次迭代的下降方向剖面ri、所述的该射线参数的背景波场以及步骤(a)中给定的背景模型，构建该射线参数对应的虚拟震源；同样基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该射线参数对应的编码震源函数进行波场正向延拓，获得该射线参数对应的扰动反偏移波场，对该射线参数对应的扰动反偏移波场进行记录采样，获得该射线参数对应的编码多分量扰动反偏移炮集；
进而获得每一个射线参数对应的编码多分量扰动反偏移炮集，所有炮的编码多分量扰动炮集构成了第i次迭代的编码多分量扰动反偏移炮集 再利用步长公式计算本次迭代的优化步长αi；
所述方程(8)中，求和变量k包含笛卡尔坐标的x，y和z三个方向， 表示第i次迭代的编码多分量扰动反偏移炮集δdi的k方向分量， 表示第i次迭代的编码多分量炮集残差Δdi的k方向分量，k＝x，y，z；
(k)：利用步骤(j)得到的优化步长αi及步骤(i)得到的下降方向剖面ri，更新第i次迭代的偏移剖面Ii＝Ii-1+αiri；利用述优化步长αi及步骤(j)得到的编码多分量扰动反偏移炮集δdi，更新第i次迭代的所述编码多分量反偏移炮集di＝di-1+αiδdi；
所述更新第i次迭代的偏移剖面Ii＝Ii-1+αiri具体为：
所述更新第i次迭代的多分量反偏移记录di＝di-1+αiδdi具体为：
(l)：第i次迭代完后，计算第i次迭代的目标泛函值fi，判断当前迭代是否满足收敛标准，如果满足则输出最新的偏移剖面为最终的偏移剖面 否则重复步骤(g)至(k)，直至获得最终偏移剖面；
所述方程(11)中，求和变量k包含笛卡尔坐标的x，y和z三个方向；
所述方程(12)中，theshold表示迭代停止的阈值标准，通常选取0.001。
2.根据权利要求1所述的一种高效率多分量地震资料真振幅偏移成像方法，其特征在于步骤(d)，步骤(e)，步骤(f)，步骤(h)和步骤(j)中所述的各向同性介质的弹性波方程为：所述方程(13)中，V＝(Vx,Vy,Vz)T表示矢量速度波场；σ＝(σxx,σyy,σzz,σxy,σyz,σxz)T表示矢量应力波场；λ和μ表示介质拉梅常数；ρ表示介质密度；x，y和z分别表示笛卡尔坐标的x，y和z三个方向。
3.根据权利要求2所述的一种高效率多分量地震资料真振幅偏移成像方法，其特征在于步骤(d)，步骤(e)和步骤(h)中所述的成像条件为：其中，α为介质纵波速度，β为介质横波速度，介质的纵横波速度与所述介质拉梅常数、介质密度之间满足一下关系：和
所述方程(14)中，在步骤(e)中， 在步骤(h)中，
所述方程(14)中， 表示震源波场中的矢量速度波场，
4.根据权利要求3所述的一种高效率多分量地震资料真振幅偏移成像方法，其特征在于步骤(d)所述的在该震源波场延拓的每一个时间步中，利用步骤(c)中的记录时间分段区间作为时间约束条件，利用弹性波场能量计算方程计算每个时刻、每个位置的波场能量值，具体为：弹性波场能量计算方程为：E＝K+W， (15)
其中，E为弹性波场的能量，K为弹性波场的动能，W为弹性波场的势能；弹性波场的动能和势能可以由其速度场和应力场计算获得，具体为：其中，U为位移场，该位移场可以通过对所述矢量速度波场V做时间积分而获得；
步骤(d)中所述的通过能量数值大小判断，获得每个时间区间内、每个位置波场能量最大值所对应的时刻，将此时刻作为该时间区间内、该位置所对应的成像时间T(x)，具体为：以步骤(c)中的第一个时间区间tR1为例，在每一个位置，如果当前时刻t在tR1之内，那么当前时刻的波场能量Et(x)若大于该时间区间内、之前所有时刻的波场能量的最大值Emax(x)，那么，Emax(x)＝Et(x)，该时间区间的成像时间T(x)＝t；否则，该时间区间内，该位置的波场能量的最大值和成像时间保持不变；对于另外两个时间区间，则同理；
步骤(d)中所述的在该时刻、该位置处，利用成像条件公式，保存成像条件中与震源波场相关的项，同时保存该时刻、该位置的波场最大能量值，具体为：成像条件中与震源波场相关的项包含： 和 共九项，计算该时刻、该位置处这九项的值并将计算的值进行保存；同时，保存该时刻、该位置的波场最大能量值Emax(x)。
5.根据权利要求3所述的一种高效率多分量地震资料真振幅偏移成像方法，其特征在于步骤(e)和步骤(h)中所述的利用步骤(d)中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤(d)中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的偏移剖面，具体为：以当前时刻t满足步骤(c)中第一个时间区间tR1为例，在每一个位置，如果当前时刻t满足：|T(x)-t|＜5，则按照所述成像条件，对保存的震源波场相关项和当前时刻所得的检波点波场的相关项应用所述成像条件，获得m1和m2，并将该时刻的成像结果叠加到该位置的成像剖面上；否则，该时刻、该位置不进行成像；对于另外两个时间区间，则同理。
6.根据权利要求5所述的一种高效率多分量地震资料真振幅偏移成像方法，其特征在于步骤(f)和步骤(j)中所述基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该射线参数对应的编码震源函数进行波场正向延拓，此处的带有所述虚拟震源的各向同性介质弹性波方程表示为：所述方程(17)中，δV＝(δVx,δVy,δVz)T表示反偏移矢量速度波场；δσ＝(δσxx,δσyy,δσzz,δσxy,δσyz,δσzx)T表示反偏移矢量应力波场；
所述方程(17)中，虚拟震源矢量F＝(Fxx,Fyy,Fzz,Fxy,Fyz,Fxz)T表示为：所述方程(18)中，中间变量m′1和m′2与所述的m1和m2满足如下关系：
一种高效率多分量地震资料真振幅偏移成像方法
技术领域
[0001] 本发明属于地球物理勘探领域，涉及多分量地震资料叠前偏移成像处理，特别涉及一种高效率、低存储多分量地震资料真振幅偏移成像方法。
背景技术
[0002] 传统油气勘探中绝大多数以纵波勘探方法为主，该方法认为地下介质中只存在纵波，然而，实际地下介质中传播的波场不仅有纵波，而且还有横波以及转换波等波型，即地震波场是弹性波场。多分量地震资料自然地包含了纵波，横波以及转换波等多种波型，采用传统单分量地震数据处理方法来处理多分量地震资料必然存在较大的误差，甚至导致处理结果错误。叠前深度偏移是地震数据处理流程中至关重要的一个环节。目前已有的多分量地震资料叠前深度偏移成像方法可以分为三类：射线类偏移方法、单程波类偏移方法和弹性波逆时偏移方法。其中，弹性波逆时偏移是精度最高，算法最稳健的一种。尽管如此，弹性波逆时偏移依然存在诸多的劣势，如：纵横波串扰造成严重的偏移假象、横波极性反转造成偏移剖面存在相消干涉、偏移剖面振幅不保真、偏移剖面分辨率低、偏移算法的计算效率低等，此外，若应用纵横波分离技术，则会引起波场的振幅及相位信息发生变化，改变波场的动力学特性。所述这些问题使得弹性波逆时偏移较难应用于实际多分量地震资料的叠前深度偏移中。为此，须建立一套新的、以“高效率、高精度、保幅”为特点的高效率多分量地震资料真振幅叠前深度偏移成像方法。
发明内容
[0003] 本发明的目的在于针对现有技术存在的上述缺陷，提供一种基于弹性波方程的、以“高精度、高效率、低存储、保幅”为特点的多分量地震资料叠前深度偏移成像方法。
[0004] 为了实现以上目的，本发明提供一种高效率多分量地震资料真振幅偏移成像方法，包括以下步骤：
[0005] (a)：通过野外观测获得的多炮多分量地震记录D(xr，t；xs)＝(Dx，Dy，Dz)，给定背景模型、预设参数，确定进行偏移成像的多分量地震记录、观测系统参数及震源函数S(x，t；xs)，其中，x表示空间位置矢量，xs表示震源位置矢量，xr表示检波点位置矢量，t表示时间变量；
[0006] (b)：基于读取的观测系统参数，设计平面波编码函数C(xs，p)，具体为：
[0007] C(xs，p)＝fδ(t-p·xs)， (1)
[0008] 其中，p为射线参数变量，fδ表示脉冲函数；此处平面波编码函数C(xs，p)等价于时移函数，时移量p·xs随震源点位置xs不同而改变，射线参数p具体为：
[0009]
[0010] 其中，α为地表纵波速度，θ为入射角；通过不同的入射角θ，可以获得不同的射线参数p；
[0011] 通过所述编码函数，并利用每一个射线参数p对观测的多炮多分量地震记录D(xr,t；xs)进行编码，获得关于每一个射线参数p的编码多分量炮集DC(xr,t；p)，具体为：
[0012] DC(xr,t；p)＝∫D(xr,t；xs)·C(xs,p)dxs； (3)
[0013] 同样，采用相同的编码函数，并利用每一个射线参数p对步骤(a)中所述观测系统中的震源函数S(x,t；xs)进行编码，获得每个射线参数p的编码震源函数SC(x,t；p)，具体为：
[0014] SC(x,t；p)＝∫S(x,t；xs)·C(xs,p)dxs； (4)
[0015] (c)：利用背景模型及观测的多炮多分量地震记录的最大记录时间tmax，设计记录时间t的分段区间，具体为：将记录时间t分为三段，用下式表示：
[0016]
[0017] 其中，tR1，tR2和tR3分别表示记录时间的三个区间范围，t1(x)和t2(x)的具体计算式为：
[0018]
[0019] 其中，αmax为背景模型的最大速度，hz代表深度，hx代表水平距离；
[0020] (d)：针对每一个射线参数p，均执行以下操作：利用每一个射线参数p，通过步骤(b)中的编码方式获得与该射线参数对应的编码震源函数，基于常规的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的每个时刻的震源波场，在该震源波场延拓的每一个时间步中，利用步骤(c)中的记录时间分段区间作为时间约束条件，利用常规的弹性波场能量计算方程计算震源波场每个时刻、每个位置的波场能量值，通过能量数值大小判断，获得每个时间区间内、每个位置波场能量最大值所对应的时刻，将此时刻作为该时间区间内、该位置所对应的成像时间T(x)，在该时刻、该位置处，利用成像条件公式，保存成像条件中与震源波场相关的项，同时保存该时刻、该位置的波场最大能量值；
[0021] 进而获得每一个射线参数对应的成像时间、成像条件中与震源波场相关的项及波场最大能量值；最终将每一个射线参数对应的每一个时间段内的波场最大能量值进行叠加，获得照明度补偿因子；
[0022] (e)：针对每一个射线参数p，均执行以下操作：对每一个射线参数p，同样基于所述的各项同性介质的弹性波方程，对该射线参数p对应的步骤(b)所述编码多分量炮集DC(xr,t；p)进行波场逆时延拓，获得该射线参数对应的每一个时刻的检波点波场，在该检波点波场延拓的每一个时间步中，利用步骤(d)中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤(d)中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的偏移剖面；
[0023] 进而获得每一个射线参数对应的偏移剖面；将所有的射线参数对应的偏移剖面进行叠加，同时利用步骤(d)中照明度补偿因子对叠加的偏移剖面进行能量补偿，获得初始偏移剖面，也即第0次迭代的偏移剖面
[0024] (f)：针对每一个射线参数p，执行以下操作：对每一个射线参数p，通过步骤(b)中的编码方式获得与该射线参数对应的编码震源函数，同样基于所述的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的背景波场；利用步骤(e)中所述的初始偏移剖面I0、所述的该射线参数的背景波场以及步骤(a)中给定的背景模型，构建该射线参数对应的虚拟震源；基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该炮点进行波场正向延拓，获得该射线参数对应的反偏移波场，对该射线参数对应的反偏移波场进行记录采样，获得该射线参数对应的编码多分量反偏移炮集；
[0025] 进而获得每一个射线参数对应的编码多分量反偏移炮集，所有炮的编码多分量炮集构成了初始编码多分量反偏移炮集，也即第0次迭代的编码多分量反偏移炮集此时，设置当前的迭代次数为i＝0；
[0026] (g)：设置当前的迭代次数i＝i+1，基于第i-1次迭代所得的所述偏移剖面和 第 i - 1 次 迭 代 所 得 的 所 述 编 码 多 分 量 反 偏 移 炮 集进行反演迭代；
[0027] 在第i次迭代中，针对每一个射线参数p，将第i-1次迭代所得的该射线参数对应的编码多分量反偏移炮集di-1与该射线参数对应的观测的编码多分量炮集DC作差，计算该射线参数对应的编码多分量炮集残差；进而获得所有射线参数对应的编码多分量炮集残差[0028] 所述的将第i-1次迭代所得的该射线参数对应的编码多分量反偏移炮集di-1与该射线参数对应的观测的编码多分量炮集DC作差，计算该射线参数对应的编码多分量炮集残差，具体为：
[0029] Δdi-1＝di-1-DC； (7)
[0030] (h)：针对每一个射线参数p，执行以下操作：对每一个射线参数p，同样基于所述的各向同性介质的弹性波方程，对该射线参数对应的编码多分量炮集残差进行波场逆时延拓，获得该射线参数对应的每一个时刻的
检波点波场，在该检波点波场延拓的每一个时间步中，利用步骤(d)中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤(d)中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的梯度剖面；
[0031] 进而获得每一个射线参数对应的梯度剖面；将所有的射线参数对应的梯度剖面进行叠加，同时利用步骤(d)中照明度补偿因子对叠加的梯度剖面进行能量补偿，获得本次迭代的梯度剖面，也即第i次迭代的梯度剖面
[0032] (i)：利用最优化反演算法，构建第i次迭代的下降方向剖面
[0033] (j)：针对每一个射线参数p，通过步骤(b)中的编码方式获得与该射线参数对应的编码震源函数，同样基于所述的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的背景波场；利用步骤(i)中所获得的第i次迭代的下降方向剖面ri、所述的该射线参数的背景波场以及步骤(a)中给定的背景模型，构建该射线参数对应的虚拟震源；同样基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该炮点进行波场正向延拓，获得该射线参数对应的扰动反偏移波场，对该射线参数对应的扰动反偏移波场进行记录采样，获得该射线参数对应的编码多分量扰动反偏移炮集；
[0034] 进而获得每一个射线参数对应的编码多分量扰动反偏移炮集，所有炮的编码多分量扰动炮集构成了第i次迭代的编码多分量扰动反偏移炮集再利用步长公式计算本次迭代的优化步长αi；
[0035] 所述的步长公式为
[0036]
[0037] 所述方程(8)中，求和变量k包含笛卡尔坐标的x，y和z三个方向；
[0038] (k)：利用步骤(j)得到的优化步长αi及步骤(i)得到的下降方向剖面ri，更新第i次迭代的偏移剖面Ii＝Ii-1+αiri；利用述优化步长αi及步骤(j)得到的编码多分量扰动反偏移炮集δdi，更新第i次迭代的所述编码多分量反偏移炮集di＝di-1+αiδdi；
[0039] 所述更新第i次迭代的偏移剖面Ii＝Ii-1+αiri具体为：
[0040]
[0041] 所述更新第i次迭代的多分量反偏移记录di＝di-1+αiδdi具体为：
[0042]
[0043] (l)：第i次迭代完后，计算第i次迭代的目标泛函值fi，判断当前迭代是否满足收敛标准，如果满足则输出最新的偏移剖面为最终的偏移剖面 否则重复步骤(g)至(k)，直至获得最终偏移剖面；
[0044] 所述计算第i次迭代的目标泛函值fi具体为：
[0045]
[0046] 所述方程(11)中，求和变量k包含笛卡尔坐标的x，y和z三个方向；
[0047] 所述的收敛标准具体为：
[0048]
[0049] 所述方程(12)中，theshold表示迭代停止的阈值标准，通常选取0.001。
[0050] 本发明实施采用的技术方案还包括：在步骤(d)，步骤(e)，步骤(f)，步骤(h)和步骤(j)中所述的各向同性介质的弹性波方程为：
[0051]
[0052] 所述方程(13)中，V＝(Vx,Vy,Vz)T表示矢量速度波场；σ＝(σxx,σyy,σzz,σxy,σyz,σxz)T表示矢量应力波场；λ和μ表示介质拉梅常数；ρ表示介质密度；x，y和z分别表示笛卡尔坐标的x，y和z三个方向。
[0053] 本发明实施采用的技术方案还包括：在步骤(d)，步骤(e)和步骤(h)中所述的成像条件为：
[0054]
[0055] 其中，α为介质纵波速度，β为介质横波速度，介质的纵横波速度与所述介质拉梅常数、介质密度之间满足一下关系：
[0056] 和
[0057] 所述方程(14)中，在步骤(e)中， 在步骤(h)中，
[0058] 所述方程(11)中， 表示震源波场中的矢量速度波场，表示检波点波场中的矢量应力波场。
[0059] 本发明实施采用的技术方案还包括：在步骤(d)所述的在该震源波场延拓的每一个时间步中，利用步骤(c)中的记录时间分段区间作为时间约束条件，利用弹性波场能量计算方程计算每个时刻、每个位置的波场能量值，具体为：弹性波场能量计算方程为：
[0060] E＝K+W， (15)
[0061] 其中，E为弹性波场的能量，K为弹性波场的动能，W为弹性波场的势能；弹性波场的动能和势能可以由其速度场和应力场计算获得，具体为：
[0062]
[0063] 其中，U为位移场，该位移场可以通过对所述矢量速度波场V做时间积分而获得；
[0064] 步骤(d)中所述的通过能量数值大小判断，获得每个时间区间内、每个位置波场能量最大值所对应的时刻，将此时刻作为该时间区间内、该位置所对应的成像时间T(x)，具体为：以步骤(c)中的第一个时间区间tR1为例，在每一个位置，如果当前时刻t在tR1之内，那么当前时刻的波场能量Et(x)若大于该时间区间内、之前所有时刻的波场能量的最大值Emax(x)，那么，Emax(x)＝Et(x)，该时间区间的成像时间T(x)＝t；否则，该时间区间内，该位置的波场能量的最大值和成像时间保持不变；对于另外两个时间区间，则同理；
[0065] 步骤(d)中所述的在该时刻、该位置处，利用成像条件公式，保存成像条件中与震源波场相关的项，同时保存该时刻、该位置的波场最大能量值，具体为：成像条件中与震源波场相关的项包含： 和 共九项，计算该时刻、该位置处这九项的值并将计算的值进行保存；同时，保存该时刻、该位置的波场最大能量值Emax(x)。
[0066] 本发明实施采用的技术方案还包括：在步骤(e)和步骤(h)中所述的利用步骤(d)中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤(d)中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的偏移剖面，具体为：以当前时刻t满足步骤(c)中第一个时间区间tR1为例，在每一个位置，如果当前时刻t满足：|T(x)-t|＜5，则按照所述成像条件，对保存的震源波场相关项和当前时刻所得的检波点波场的相关项应用所述成像条件，获得m1和m2，并将该时刻的成像结果叠加到该位置的成像剖面上；否则，该时刻、该位置不进行成像；对于另外两个时间区间，则同理。
[0067] 本发明实施采用的技术方案还包括：在步骤(f)和步骤(j)中所述基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该炮点进行波场正向延拓，此处的带有所述虚拟震源的各向同性介质弹性波方程表示为：
[0068]
[0069] 所述方程(17)中，δV＝(δVx,δVy,δVz)T表示反偏移矢量速度波场；δσ＝(δσxx,δσyy,δTσzz,δσxy,δσyz,δσzx) 表示反偏移矢量应力波场；
[0070] 所述方程(17)中，虚拟震源矢量F＝(Fxx,Fyy,Fzz,Fxy,Fyz,Fxz)T表示为：
[0071]
[0072] 所述方程(18)中，中间变量m′1和m′2与所述的m1和m2满足如下关系：
[0073] 和
[0074] 在步骤(f)中， 在步骤(j)中，
[0075] 本发明由于采取以上技术方案，其具有以下优点：1)本发明将反演的思想引入弹性波逆时偏移中，与常规弹性波逆时偏移相比，本发明可以获得高精度、高分辨率、高信噪比、振幅保真的偏移剖面；2)本发明所用成像条件自然地校正了横波极性，无需极性反转校正运算便可有效克服横波极性反转造成的偏移剖面相消干涉；3)本发明采用的成像方式，无需额外的纵横波分离运算，而且完整地保持了纵横波的矢量特性、振幅以及相位特性，有效地消除了纵横波之间串扰造成的偏移假象，极大地提高了成像的精度；4)本发明所得成像剖面物理意义非常明确，便于后期地质解译；5)本发明可以广泛用于油气勘探领域中，特别是对于复杂构造深部的成像效果更加明显。
附图说明
[0076] 图1是本发明提供的高效率多分量地震资料真振幅偏移成像方法流程示意图；
[0077] 图2是本发明提供的二维洼陷模型图；其中，图2(a)是纵波速度模型；图2(b)是横波速度模型；图2(c)是背景纵波速度模型；图2(d)是背景横波速度模型；
[0078] 图3是本发明提供的切除直达波后的编码炮集，对应的射线参数p＝0μs/m：其中，图3(a)是编码的观测记录水平分量；图3(b)是编码的观测记录垂直分量；图3(c)是利用本发明方法所得的编码的初始反偏移记录水平分量；图3(d)是利用本发明方法所得的编码的初始反偏移记录垂直分量；图3(e)是利用本发明方法所得的编码的最终反偏移记录水平分量；图3(f)是利用本发明方法所得的编码的最终反偏移记录垂直分量；
[0079] 图4是图2所示二维洼陷模型的多炮叠加偏移剖面：其中，图4(a)是利用传统方法所得的Iα剖面，图4(b)是利用传统方法所得的Iβ剖面，图4(c)是利用本发明方法所得的Iα剖面，图4(d)是利用本发明方法所得的Iβ剖面；
[0080] 图5是本发明提供的Marmousi-ii模型：其中，图5(a)是纵波速度模型；图5(b)是横波速度模型；图5(c)是背景纵波速度模型；图5(d)是背景横波速度模型；
[0081] 图6是图5所示Marmousi-ii模型的多炮叠加偏移剖面：其中，图6(a)是利用传统方法所得的Iα剖面；图6(b)是利用传统方法所得的Iβ剖面；图6(c)是利用本发明方法所得的Iα剖面；图6(d)是利用本发明方法所得的Iβ剖面；
[0082] 图7是本发明提供的SEG/EAGE Salt模型：其中，图7(a)是纵波速度模型；图7(b)是横波速度模型；图7(c)是背景纵波速度模型；图7(d)是背景横波速度模型；
[0083] 图8是图7所示SEG/EAGE Salt模型的多炮叠加偏移剖面：其中，图8(a)是利用传统方法所得的Iα剖面；图8(b)是利用传统方法所得的Iβ剖面；图8(c)是利用本发明方法所得的Iα剖面；图8(d)是利用本发明方法所得的Iβ剖面；
具体实施方式
[0084] 为了是本发明的目的，技术方案及优点更加清楚明白，一下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用于解释本发明，并不用于限定本发明。
[0085] 请参阅图1，是本发明提供的高效率多分量地震资料真振幅偏移成像方法流程示意图。本发明实施例的高效率多分量地震资料真振幅偏移成像方法包括以下步骤：
[0086] 步骤S100：通过野外观测获得的多炮多分量地震记录D(xr,t；xs)＝(Dx,Dy,Dz)，给定背景模型、预设参数，确定进行偏移成像的多分量地震记录、观测系统参数及震源函数S(x,t；xs)，其中，x表示空间位置矢量，xs表示震源位置矢量，xr表示检波点位置矢量，t表示时间变量；
[0087] 步骤S200：基于读取的观测系统参数，设计平面波编码函数C(xs,p)，具体为：
[0088] C(xs,p)＝fδ(t-p·xs)， (1)
[0089] 其中，p为射线参数变量，fδ表示脉冲函数；此处平面波编码函数C(xs,p)等价于时移函数，时移量p·xs随震源点位置xs不同而改变，射线参数p具体为：
[0090]
[0091] 其中，α为地表纵波速度，θ为入射角；通过不同的入射角θ，可以获得不同的射线参数p；
[0092] 通过所述编码函数，并利用每一个射线参数p对观测的多炮多分量地震记录D(xr,t；xs)进行编码，获得关于每一个射线参数p的编码多分量炮集DC(xr,t；p)，具体为：
[0093] DC(xr,t；p)＝∫D(xr,t；xs)·C(xs,p)dxs； (3)
[0094] 同样，采用相同的编码函数，并利用每一个射线参数p对步骤(a)中所述观测系统中的震源函数S(x,t；xs)进行编码，获得每个射线参数p的编码震源函数SC(x,t；p)，具体为：
[0095] SC(x,t；p)＝∫S(x,t；xs)·C(xs,p)dxs； (4)
[0096] 步骤S300：利用背景模型及观测的多炮多分量地震记录的最大记录时间tmax，设计记录时间t的分段区间，具体为：将记录时间t分为三段，用下式表示：
[0097]
[0098] 其中，tR1，tR2和tR3分别表示记录时间的三个区间范围，t1(x)和t2(x)的具体计算式为：
[0099]
[0100] 其中，αmax为背景模型的最大速度，hz代表深度，hx代表水平距离；
[0101] 步骤S400：针对每一个射线参数p，均执行以下操作：利用每一个射线参数p，通过步骤S200中的编码方式获得与该射线参数对应的编码震源函数，基于常规的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的每个时刻的震源波场，在该震源波场延拓的每一个时间步中，利用步骤S300中的记录时间分段区间作为时间约束条件，利用常规的弹性波场能量计算方程计算震源波场每个时刻、每个位置的波场能量值，通过能量数值大小判断，获得每个时间区间内、每个位置波场能量最大值所对应的时刻，将此时刻作为该时间区间内、该位置所对应的成像时间T(x)，在该时刻、该位置处，利用成像条件公式，保存成像条件中与震源波场相关的项，同时保存该时刻、该位置的波场最大能量值；进而获得每一个射线参数对应的成像时间、成像条件中与震源波场相关的项及波场最大能量值；最终将每一个射线参数对应的每一个时间段内的波场最大能量值进行叠加，获得照明度补偿因子；
[0102] 在步骤S400中所述的常规的各向同性介质的弹性波方程为：
[0103]
[0104] 所述方程(9)中，V＝(Vx,Vy,Vz)T表示矢量速度波场；σ＝(σxx,σyy,σzz,σxy,σyz,σxz)T表示矢量应力波场；λ和μ表示介质拉梅常数；ρ表示介质密度；x，y和z分别表示笛卡尔坐标的x，y和z三个方向；
[0105] 在步骤S400中所述的成像条件为：
[0106]
[0107] 其中，α为介质纵波速度，β为介质横波速度，介质的纵横波速度与所述介质拉梅常数、介质密度之间满足一下关系：
[0108] 和
[0109] 在步骤S400中所述的在该震源波场延拓的每一个时间步中，利用步骤S300中的记录时间分段区间作为时间约束条件，利用常规的弹性波场能量计算方程计算震源波场每个时刻、每个位置的波场能量值，具体为：弹性波场能量计算方程为：
[0110] E＝K+W， (15)
[0111] 其中，E为弹性波场的能量，K为弹性波场的动能，W为弹性波场的势能；弹性波场的动能和势能可以由其速度场和应力场计算获得，具体为：
[0112]
[0113] 其中，U为位移场，该位移场可以通过对所述矢量速度波场V做时间积分而获得；
[0114] 所述的通过能量数值大小判断，获得每个时间区间内、每个位置波场能量最大值所对应的时刻，将此时刻作为该时间区间内、该位置所对应的成像时间T(x)，具体为：以步骤S300中的第一个时间区间tR1为例，在每一个位置，如果当前时刻t在tR1之内，那么当前时刻的波场能量Et(x)若大于该时间区间内、之前所有时刻的波场能量的最大值Emax(x)，那么，Emax(x)＝Et(x)，该时间区间的成像时间T(x)＝t；否则，该时间区间内，该位置的波场能量的最大值和成像时间保持不变；对于另外两个时间区间，则同理；
[0115] 所述的在该时刻、该位置处，利用成像条件公式，保存成像条件中与震源波场相关的项，同时保存该时刻、该位置的波场最大能量值，具体为：成像条件中与震源波场相关的项包含： 和 共九项，计算该时刻、该位置处这九项的值并将计算的值进行保存；同时，保存该时刻、该位置的波场最大能量值Emax(x)；
[0116] 步骤S500：针对每一个射线参数p，均执行以下操作：对每一个射线参数p，同样基于所述的各项同性介质的弹性波方程，对该射线参数p对应的所述编码多分量炮集DC(xr,t；p)进行波场逆时延拓，获得该射线参数对应的每一个时刻的检波点波场，在该检波点波场延拓的每一个时间步中，利用步骤S400中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤S400中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的偏移剖面；进而获得每一个射线参数对应的偏移剖面；将所有的射线参数对应的偏移剖面进行叠加，同时利用步骤S400中照明度补偿因子对叠加的偏移剖面进行能量补偿，获得初始偏移剖面，也即第0次迭代的偏移剖面
[0117] 在步骤S500中所述的利用步骤S400中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤S400中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的偏移剖面，具体为：以当前时刻t满足步骤S300中第一个时间区间tR1为例，在每一个位置，如果当前时刻t满足：|T(x)-t|＜5，则按照所述成像条件，对保存的震源波场相关项和当前时刻所得的检波点波场的相关项应用所述成像条件，获得m1和m2，并将该时刻的成像结果叠加到该位置的成像剖面上；否则，该时刻、该位置不进行成像；对于另外两个时间区间，则同理；
[0118] 在步骤S500中所述应用成像条件，具体为基于方程(14)，偏移剖面具体为：
[0119] 步骤S600：针对每一个射线参数p，执行以下操作：对每一个射线参数p，通过步骤S200中的编码方式获得与该射线参数对应的编码震源函数，同样基于所述的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的背景波场；利用步骤S500中所述的初始偏移剖面I0、所述的该射线参数的背景波场以及步骤S100中给定的背景模型，构建该射线参数对应的虚拟震源；基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该炮点进行波场正向延拓，获得该射线参数对应的反偏移波场，对该射线参数对应的反偏移波场进行记录采样，获得该射线参数对应的编码多分量反偏移炮集；进而获得每一个射线参数对应的编码多分量反偏移炮集，所有炮的编码多分量炮集构成了初始编码多分量反偏移炮集，也即第0次迭代的编码多分量反偏移炮集 此时，设置当前的迭代次数为i＝0；
[0120] 在步骤S600中所述的带有虚拟震源的各向同性介质弹性波方程表示为：
[0121]
[0122] 所述方程(17)中，δV＝(δVx,δVy,δVz)T表示反偏移矢量速度波场；δσ＝(δσxx,δσyy,δσzz,δσxy,δσyz,δσzx)T表示反偏移矢量应力波场；所述方程(17)中，虚拟震源矢量F＝(Fxx,Fyy,Fzz,Fxy,Fyz,Fxz)T表示为：
[0123]
[0124] 所述方程(18)中，中间变量m′1和m′2与所述的m1和m2满足如下关系：
[0125] 和
[0126] 在步骤S600中，所述方程(18)中，
[0127] 步骤S700：设置当前的迭代次数i＝i+1，基于第i-1次迭代所得的所述偏移剖面和 第 i - 1 次 迭 代 所 得 的 所 述 编 码 多 分 量 反 偏 移 炮 集进行反演迭代；
[0128] 在第i次迭代中，针对每一个射线参数，将第i-1次迭代所得的该射线参数对应的编码多分量反偏移炮集di-1与该射线参数对应的观测的编码多分量炮集DC作差，计算该射线参数对应的编码多分量炮集残差；进而获得所有射线参数对应的编码多分量炮集残差[0129] 所述的将第i-1次迭代所得的该射线参数对应的编码多分量反偏移炮集di-1与该射线参数对应的观测的编码多分量炮集DC作差，计算该射线参数对应的编码多分量炮集残差，具体为：
[0130] Δdi-1＝di-1-DC； (7)
[0131] 步骤S800：针对每一个射线参数p，执行以下操作：对每一个射线参数p，同样基于所述的各向同性介质的弹性波方程，对该射线参数对应的编码多分量炮集残差进行波场逆时延拓，获得该射线参数对应的每一个时刻的
检波点波场，在该检波点波场延拓的每一个时间步中，利用步骤S400中保存的每个时间段、每个位置的成像时间T(x)，判断当前时刻是否满足成像时间，如果不满足，则不进行成像，如果满足，则按照成像条件，利用步骤S400中保存的成像条件中与震源波场相关的项以及检波点波场进行成像，获得该射线参数对应的梯度剖面；进而获得每一个射线参数对应的梯度剖面；将所有的射线参数对应的梯度剖面进行叠加，同时利用步骤S400中照明度补偿因子对叠加的梯度剖面进行能量补偿，获得本次迭代的梯度剖面，也即第i次迭代的梯度剖面
[0132] 在步骤S800中，所述应用成像条件，具体为基于所述方程(14)，其中，
[0133] 步骤S900：利用最优化反演算法，构建第i次迭代的下降方向剖面
[0134] 步骤S1000：针对每一个射线参数p，通过步骤S200中的编码方式获得与该射线参数对应的编码震源函数，同样基于所述的各向同性介质的弹性波方程，利用上述给定的背景模型对该编码震源进行波场正向延拓，获得该射线参数对应的背景波场；利用步骤S900中所获得的第i次迭代的下降方向剖面ri、所述的该射线参数的背景波场以及步骤S100中给定的背景模型，构建该射线参数对应的虚拟震源；同样基于带有所述虚拟震源的各向同性介质弹性波方程，利用给定的背景模型对该炮点进行波场正向延拓，获得该射线参数对应的扰动反偏移波场，对该射线参数对应的扰动反偏移波场进行记录采样，获得该射线参数对应的编码多分量扰动反偏移炮集；进而获得每一个射线参数对应的编码多分量扰动反偏移炮集，所有炮的编码多分量扰动炮集构成了第i次迭代的编码多分量扰动反偏移炮集再利用步长公式计算本次迭代的优化步长αi；
[0135] 在步骤S1000中所述的利用步骤S900中所获得的第i次迭代的下降方向剖面ri、所述的该射线参数的背景波场以及步骤S100中给定的背景模型，构建该射线参数对应的虚拟震源，具体为基于所述方程(17)，其中，
[0136] 在步骤S1000中所述的步长公式为
[0137]
[0138] 所述方程(8)中，求和变量k包含笛卡尔坐标的x，y和z三个方向；
[0139] 步骤S1100：利用步骤S1000得到的优化步长αi及步骤S900得到的下降方向剖面ri，更新第i次迭代的偏移剖面Ii＝Ii-1+αiri；利用述优化步长αi及步骤S1000得到的编码多分量扰动反偏移炮集δdi，更新第i次迭代的所述编码多分量反偏移炮集di＝di-1+αiδdi；
[0140] 所述更新第i次迭代的偏移剖面Ii＝Ii-1+αiri具体为：
[0141]
[0142] 所述更新第i次迭代的多分量反偏移记录di＝di-1+αiδdi具体为：
[0143]
[0144] 步骤S1200：第i次迭代完后，计算第i次迭代的目标泛函值fi，判断当前迭代是否满足收敛标准，如果满足则输出最新的偏移剖面为最终的偏移剖面 否则重复步骤S700至S1100，直至获得最终偏移剖面；
[0145] 所述计算第i次迭代的目标泛函值fi具体为：
[0146]
[0147] 所述方程(11)中，求和变量k包含笛卡尔坐标的x，y和z三个方向；
[0148] 所述的收敛标准具体为：
[0149]
[0150] 所述方程(12)中，theshold表示迭代停止的阈值标准，通常选取0.001。
[0151] 本发明涉及一种高效率多分量地震资料真振幅偏移成像方法，本发明是在弹性波逆时偏移方法的基础上进行了改进，特征是直接以多分量地震资料为输入，通过编码算法、分段成像等措施，再在反演的框架下，实现基于反演的弹性波逆时偏移。本发明将编码、时间分段及反演的思想引入弹性波逆时偏移中，与常规弹性波逆时偏移相比，本发明可以获得高精度且振幅保真的叠前深度偏移剖面，同时算法的效率很高。本发明完整地保持纵横波的矢量特性、振幅以及相位特性，偏移剖面具有保真性；有效地消除了纵横波之间串扰造成的偏移假象，提高了成像的精度；有效降低了弹性波逆时偏移的存储量，提高了计算效率。
[0152] 为进一步说明本发明的可行性和有效性，下面举三个实例：
[0153] 实例1：
[0154] 图2是二维洼陷模型图；其中，图2(a)是纵波速度模型；图2(b)是横波速度模型；图2(c)是背景纵波速度模型；图2(d)是背景横波速度模型。在此模型上设置21个射线参数，起始射线参数p＝-117μs/m，射线参数间隔为p＝15μs/m，所用震源为爆炸震源，震源子波设定为雷克子波，主频为20赫兹。图3是切除直达波后的编码炮集，对应的射线参数p＝0μs/m：其中，图3(a)是编码的观测记录水平分量；图3(b)是编码的观测记录垂直分量；图3(c)是利用本发明方法所得的编码的初始反偏移记录水平分量；图3(d)是利用本发明方法所得的编码的初始反偏移记录垂直分量；图3(e)是利用本发明方法所得的编码的最终反偏移记录水平分量；图3(f)是利用本发明方法所得的编码的最终反偏移记录垂直分量。从图3中可以看出，本发明方法所得反偏移记录与观测记录匹配非常好，这验证了本发明方法的有效性。图α
4是图2所示二维洼陷模型的多炮叠加偏移剖面：其中，图4(a)是利用传统方法所得的I 剖面，图4(b)是利用传统方法所得的Iβ剖面，图4(c)是利用本发明方法所得的Iα剖面，图4(d)是利用本发明方法所得的Iβ剖面。从图4(a)和图4(b)中可以看出，两个剖面均存在比较明显的串扰噪声及假象，剖面的分辨率较低，振幅不均衡，从图4(c)和图4(d)中可以看出，本发明方法所得的剖面很好地压制了串扰噪声及假象，偏移剖面具有很高的精度、分辨率和信噪比，而且振幅是保真的。由图3和图4可以证明，本发明方法能够获得高质量的偏移剖面，同时也证明了本发明方法的可行性及有效性。
[0155] 实例2：
[0156] 图5是Marmousi-ii模型：其中，图5(a)是纵波速度模型；图5(b)是横波速度模型；图5(c)是背景纵波速度模型；图5(d)是背景横波速度模型。该模型是验证各种偏移方法成像效果的国际标准模型之一。在此模型上设置31个射线参数，起始射线参数p＝-333μs/m，射线参数间隔为p＝20μs/m，所用震源为爆炸震源，震源子波设定为雷克子波，主频为15赫兹。图6是图5所示Marmousi-ii模型的多炮叠加偏移剖面：其中，图6(a)是利用传统方法所得的Iα剖面；图6(b)是利用传统方法所得的Iβ剖面；图6(c)是利用本发明方法所得的Iα剖面；图6(d)是利用本发明方法所得的Iβ剖面。从图6(a)和图6(b)中可以看出，传统方法所得的Iα剖面和Iβ剖面在压制纵横波串扰方面能力不足，而且成像剖面振幅严重失衡，深部剖面相对振幅并不保真。与其相比，本发明方法所得成像剖面纵横波串扰几乎被完全消除，偏移剖面的分辨率和精度更高，信噪比也更好，振幅均衡性更佳，因此可以得出结论，本发明方法成像效果更好。
[0157] 实例3：
[0158] 图7是SEG/EAGE Salt模型：其中，图7(a)是纵波速度模型；图7(b)是横波速度模型；图7(c)是背景纵波速度模型；图7(d)是背景横波速度模型。该模型是验证各种偏移方法成像效果的国际标准盐丘模型之一。在此模型上设置21个射线参数，起始射线参数p＝-333μs/m，射线参数间隔为p＝30μs/m，所用震源为爆炸震源，震源子波设定为雷克子波，主频为15赫兹。图8是图7所示SEG/EAGESalt模型的多炮叠加偏移剖面：其中，图8(a)是利用传统方α β α
法所得的I剖面；图8(b)是利用传统方法所得的I 剖面；图8(c)是利用本发明方法所得的I剖面；图8(d)是利用本发明方法所得的Iβ剖面。从图8中可以看出，本发明方法所得偏移剖面具有更高的分辨率、更好的信噪比、精度也更高，此外我们也可以清晰地看出本发明方法所得剖面具有很高的振幅保真性。
[0159] 以上描述和实施例仅用于说明本发明，凡是在本发明技术方案的基础上进行的等同变化和改进，均应包含在本发明的保护范围之内。
