基于分步傅里叶算法的叠前时间偏移方法

摘要

本发明提供一种基于分步傅里叶算法的叠前时间偏移方法，该方法包括：输入预处理后的地震数据；通过交互方式拾取均方根速度；计算成像点到震源点的射线传播时间；将时间-空间域的地震数据变换到频率-波数；得到背景波场的延拓值；将延拓后的波场经傅里叶变换到频率-空间域，通过时间域的分步傅里叶时移算子进行时移运算，实现对扰动波场的延拓；以及利用激发时间成像条件，对频率-空间域地震波场作相位移动并进行有效频段的积分，得到炮记录的地震成像结果。该基于分步傅里叶算法的叠前时间偏移方法可以应用均方根速度进行偏移，在不增加速度建模难度的情况下，提高了叠前时间偏移对速度纵横向变化的适应能力以及对偏移成像的保幅能力。

说明书

技术领域

本发明涉及地震资料处理中的偏移成像领域，特别是涉及到一种基于分步傅里叶算法的叠前时间偏移方法。

背景技术

地震偏移是一种将地震信息进行重排的反演运算，以便使地震波能量归位到其空间的真实位置，获取地下真实构造形态。波动方程叠前偏移成像主要分为叠前时间偏移成像和叠前深度偏移成像两类。叠前深度偏移对地下速度横向变化适应性强，可提供深度域的优质成像结果，但它对速度模型的精度要求高且计算量大，因此，是一个高投入高回报的方案。叠前时间偏移成像计算效率高，同时其对速度模型应用条件的要求低，在构造复杂、速度横向变化不大的情况下有较好的成像效果。因此，叠前时间偏移是现阶段的实用技术。

叠前时间偏移成像按照算法实现的原理可分为两大类：基于射线理论的偏移方法和基于波场延拓的偏移方法。常规的射线法叠前时间偏移采用均方根速度，它能够实现纵横向变速，计算量较小，而且计算均方根速度并提高其精度比较容易，但是其数值响应特性不如波场法。尽管如此，目前主流的地震数据处理软件里叠前时间偏移(PSTM)仍然只有射线法而没有波场法。这是因为波场法PSTM在递推过程中引用层速度，计算时间-层速度模型具有很大难度。因此，在实践上，如果需要花费同等代价的话，人们宁愿建立深度-层速度模型并选择精度更高的叠前深度偏移，所以，引用层速度模型成为波场法PSTM的最大障碍。

随着油气勘探难度的增加，地震资料保幅性处理越来越重要。近年来出现了平面波、相移加内插和常速扫描等波场类叠前时间偏移技术，但是其诸多指标无法满足实用要求。

在频率-波数域，Dubrulle(1983)、Ekren和Ursin(1999)实现了对共炮检距的平面波波场延拓。但是，理论上共偏移距道集平面波叠前时间偏移只能处理具有垂向速度变化的介质，而且对于地震波动力学特征的保持还需进一步提高。

从单道格林函数的积分公式(Ehinger A.,Lailly P.,and Marfurt K.J.1996)出发，Jianhua P.和Dan N.(2008)借用相移法得到具有波场法响应特性的Kirchhoff型叠前时间偏移结果。这种方式既能保留射线法的主要优点，还能改善振幅和波形特性。但是，以道为单位进行波场延拓会导致非常可观的计算量。

王有新(2008)通过相移加内插的方法实现对单炮记录的叠前时间偏移。该方法成像角度高，但是计算量大，通过线性插值计算波场，不能适应横向速度的较大变化。

Sergey fomel(2013)提出将波动方程从笛卡尔坐标系下转到成像射线坐标系下进行求解的方法，这种方法避免了常规时间偏移中对格林函数计算的近似求解，提高了时间偏移计算精度，但是其对速度的横向变化适应性差。

由于波场法具有良好的数值特性，所以，波场法叠前时间偏移具有重要的实用价值，许多研究工作试图利用波场法的计算特性，但是受制于引用均方根速度与适应速度变化两个问题，难以满足实用要求。为此我们发明了一种新的基于分步傅里叶算法的叠前时间偏移方法，解决了这两个技术问题。

发明内容

本发明的目的是提供一种可以产生更加真实可靠的成像结果的基于分步傅里叶算法的叠前时间偏移方法。

本发明的目的可通过如下技术措施来实现：基于分步傅里叶算法的叠前时间偏移方法，该基于分步傅里叶算法的叠前时间偏移方法包括：步骤1，输入预处理后的地震数据p_u(t,x,y)；步骤2，利用预处理后的地震数据p_u(t,x,y)计算得到均方根速度谱，通过交互方式拾取均方根速度，并获得工区范围内均方根速度场v_rms(t,x,y)；步骤3，基于均方根速度场v_rms(t,x,y)，计算成像点到震源点的射线传播时间τ；步骤4，利用傅里叶变换，将时间-空间域的地震数据p_u(t,x,y)变换到频率-波数域p^*_u(w,k_x,k_y)，其中，p^*_u(w,k_x,k_y)＝FFT(p_u(t,x,y))；步骤5，利用分步傅里叶相移算子将p^*_u(w,k_x,k_y)直接延拓到地下某成像点的时间τ上，得到背景波场的延拓值；步骤6，将延拓后的波场经傅里叶变换到频率-空间域，通过时间域的分步傅里叶时移算子 进行时移运算，实现对扰动波场的延拓；以及步骤7，根据成像点到震源点的射线传播时间τ，利用激发时间成像条件，对频率-空间域地震波场作相位移动并进行有效频段的积分，得到炮记录的地震成像结果。

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

在步骤1中，对深度域层速度模型正演得到的共炮集p_u(t,x,y)进行去噪、静校正这些常规预处理。

在步骤2中，通过交互方式拾取种子点位置的均方根速度，经过三维插值平滑获得目标区的均方根速度场v_rms(t,x,y)，所述均方根速度场用于对地震数据的分步傅里叶叠前时间偏移。

在步骤2中，利用预处理后的地震数据p_u(t,x,y)计算均方根速度谱，通过交互方式，按照100米等间隔拾取种子点位置的均方根速度，然后通过插值平滑获得间隔为10米的模型的时间域均方根速度场v_rms(t,x,y)，即叠前时间偏移所需的速度场。

在步骤3中，基于均方根速度场v_rms(t,x,y)，采用射线法计算出成像点到震源点的射线传播时间其中，τ为成像点的时间，x为成像点的位置，x_s为震源点的位置，v_rms为均方根速度。

在步骤5中，通过应用时间域的均方根速度场，通过单步延拓得到，在延拓时间τ上针对背景速度的相移项，所述相移算子表达式如下：

${p^{*}}_u(w,k_x,k_y,\tau )={p^{*}}_u(w,k_x,k_y,0)·e^{i\tau \sqrt{w^2{s^2}_{rms0}\left(\tau \right)-k_x^2-k_y^2}}$

其中，v_rms0(τ)为τ时刻的背景均方根场速度值，s_rms0(τ)为背景场的慢度，且s_rms0(τ)＝1/v_rms0(τ)，w为频率，k_x为x方向的波数，k_y为y方向的波数。

在步骤6中，将步骤5计算得到的波场经傅里叶反变换转换到频率-空间域，在此基础上通过时间域的分步傅里叶时移算子计算得到针对扰动场的时移项，应用时移算子，实现对扰动场的延拓，其表达式为：

$p{{\_}^{*}}_u(w,x,y,\tau )={p^{*}}_u(w,x,y,\tau )·e^{{\mathrm{iw\Delta s}}_{rms}(\tau ,x,y)}$

其中，Δs_rms(τ,x,y)为扰动场的慢度，至此，完成对检波点一侧地震记录在成像点位置处τ时刻的波场延拓。

在步骤7中，对步骤6得到的频率-空间域地震波场延拓值，按照激发时间成像条件实施相位移动并对有效频段进行积分，得到最终的时间-空间域地震波场值，其成像条件表达式为：

$I(x,y,t_s)={\left(\frac{1}{2\pi }\right)}^2{\int }_{w1}^{w2}I(x,y,\tau ,w)e^{iw(\tau -t_s)}dw$

其中，I为相位移动后的时间-空间域地震波场，w1至w2为地震波的有效频带范围。

步骤7还包括：将炮记录偏移成像转换成共成像点道集，合并相同炮检距；将所有单炮成像结果按照空间位置叠加，或者对共成像点道集直接叠加，获得目标区的叠前时间偏移成像结果。

本发明提出一种基于分步傅里叶算法的叠前时间偏移方法，其中的关键是根据均方根速度实现波场外推。通常，射线法引用均方根速度，波场法引用层速度，而建立层速度模型有很大难度，所以，叠前时间偏移主要采用射线法。然而，射线法叠前时间偏移的数值响应特性不如波场法且保幅性较差，为此，能够简单快捷的应用均方根速度实现波场法偏移是叠前时间偏移需要解决的一个重要问题。在波场类算法中，相移法、平面波法、分步傅里叶及屏算法能把逐步递推变成一步外推，从而引用均方根速度，其中分步傅里叶算法简洁并具有一阶横向变速能力，恰好符合叠前时间偏移的前提条件。本发明基于分步傅里叶算法，应用均方根速度实现了波场法叠前时间偏移，其主要特点是：①引用均方根速度，通过单步波场外推实现速度纵向变化，避免了层速度在逐层向下递推时引入的误差；②通过相对时移补偿实现横向速度一阶变化；③通过炮点射线法和检波点波场法的混合，实现射线法高效性与波场法高保真的优势互补。无论是射线法还是波场法，在理论上既可以用于叠前时间偏移，也可以用于叠前深度偏移。实际中，如果已知层速度，叠前深度偏移比叠前时间偏移的成像效果好，但是，由于计算准确的层速度比计算均方根速度困难，因此，采用均方根速度场实现波场法叠前时间偏移不 失为一种高效快捷、保真性高的成像方法。

本发明的基于分步傅里叶算法的叠前时间偏移方法，波场法叠前时间偏移技术应用前，首先对地震数据进行常规去噪、静校正等预处理，通过分析地震数据，得到准确的均方根速度场。然后，利用射线法，计算震源点到成像点的射线时间。此时，将地震数据转换到频率-波数域，利用背景均方根速度对频率-波数域的地震数据进行一次相移，实现对背景场的波场延拓；将延拓后的地震数据变换到频率-空间域，利用扰动均方根速度场进行一次时移，实现对扰动场的波场延拓；根据震源点到成像点的传播时间，对延拓后的频率-空间域波场实施相位移动；经波场相位移动后，对地震波场数据进行有效频段的积分，得到最终的地震波场延拓值；最后，将所有单炮成像结果按照空间位置叠加，或者对共成像点道集直接叠加，获得目标区的叠前时间偏移成像结果。

该基于分步傅里叶算法的叠前时间偏移方法采用低成本的均方根速度场实现了基于分步傅里叶算法的叠前时间偏移，通过单步波场延拓避免了层速度在向下递推时重复不断的正反变换，在不增加计算成本的同时提高偏移精度并改善振幅特性。与此同时，通过射线法和波场法的混合运算，既体现了波场法的良好振幅响应特性也采纳了射线法适应速度时空变化的优点，取得优势互补的效果，大幅提高了计算效率。与常规叠前时间偏移算法相比，本发明可以应用均方根速度进行偏移，在不增加速度建模难度的情况下，提高叠前时间偏移对速度纵横向变化的适应能力以及对偏移成像的保幅能力。

附图说明

图1为本发明的基于分步傅里叶算法的叠前时间偏移方法的具体实施例的流程图；

图2为本发明的具体实施例中marmousi模型对应的深度-层速度示意图；

图3为本发明的具体实施例中marmousi模型对应的时间-均方根速度示意图；

图4为本发明的具体实施例中marmousi模型的射线法叠前时间偏移结果示意图；

图5为本发明的具体实施例中marmousi模型的基于分步傅里叶算法的波场法叠前时间偏移结果示意图；

具体实施方式

为使本发明实施例的目的、技术方案和优点能更明显易懂，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述。

如图1所示，图1为本发明的基于分步傅里叶算法的叠前时间偏移方法的流程图。

在步骤101，输入地震数据，该地震数据为经过常规去噪、静校正等预处理的共炮点道集。在一实施例中，对由图2所示的marmousi深度域层速度模型正演得到的共炮集p_u(t,x,y)进行去噪、静校正等常规预处理。流程进入到步骤102。

在步骤102，拾取均方根速度，利用预处理后的地震数据计算得到均方根速度谱，通过交互方式拾取种子点位置的均方根速度，经过三维插值平滑获得目标区的均方根速度场，所述均方根速度场用于对地震数据的分步傅里叶叠前时间偏移。在一实施例中，利用预处理后的地震数据p_u(t,x,y)计算均方根速度谱，通过交互方式，按照100米等间隔拾取种子点位置的均方根速度，然后通过插值平滑获得如图3所示的间隔为10米的marmousi模型的时间域均方根速度场v_rms(t,x,y)，即叠前时间偏移所需的速度场。流程进入到步骤103。

在步骤103，利用射线旅行时计算方法并采用所述均方根速度场，计算成像点到震源点的射线传播时间。在一实施例中，基于均方根速度场v_rms(t,x,y)，采用射线法计算出成像点到震源点的射线传播时间其中，τ为成像点的时间，x为成像点的位置，x_s为震源点的位置，v_rms为均方根速度，流程进入到步骤104。

在步骤104，利用傅里叶变换，将时间-空间域的共炮点道集变换到频率-波数域。在一实施例中，通过二维傅里叶变化,将经过预处理的检波点一侧地震数据转变为频率-波数域的p^*_u(w,k_x,k_y)，其中p^*_u(w,k_x,k_y)＝FFT(p_u(t,x,y))。流程进入到步骤105。

在步骤105，利用分步傅里叶相移算子并采用所述均方根速度场对所述共炮点道集地震数据进行波场外推，得到外推后的频率-波数域地震波场数据。

对地震数据应用一个由背景参考慢度s_rms0(τ)获得的相移算子。在叠前深度偏移中，相移延拓的步长为一固定的深度值Δz，需要通过逐层递推获得成像点位置的波场延拓值。在本发明实施例中，通过应用时间域的均方根速度场，通过单步延拓得到，在延拓时间τ上针对背景速度的相移项，所述相移算子表达式如下：

$p{*}_u(w,k_x,k_y,\tau )=p{*}_u(w,k_x,k_y,0).e^{i\tau \sqrt{w^2{s^2}_{\mathrm{rms}0}\left(\tau \right)-k_x^2-k_y^2}}$

其中，v_rms0(τ)为τ时刻的背景均方根场速度值，s_rms0(τ)为背景场的慢度，且s_rms0(τ)＝1/v_rms0(τ)，w为频率，k_x为x方向的波数，k_y为y方向的波数。

应用上述公式，就可以将上行波p^*_u(w,k_x,k_y)直接延拓到地下某成像点的时间τ上，得到背景波场的延拓值。流程进入到步骤106。

在步骤106，利用傅里叶反变换，将外推后的频率-波数域地震波场数据变换到频率-空间域，得到变换后的频率-空间域地震波场数据。利用分步傅里叶时移算子并采用所述均方根速度场对变换后的频率-空间域地震波场数据进行延拓处理，得到二次延拓的频率-空间域地震波场数据。

将步骤105计算得到的波场经傅里叶反变换转换到频率-空间域，在此基础上通过时间域的分步傅里叶时移算子计算得到针对扰动场的时移项，应用时移算子，实现对扰动场的延拓，其表达式为：

$p{*}_{-u}(w,x,y,\tau )=p{*}_u(w,x,y,\tau )·e^{{\mathrm{iw\Delta s}}_{rms}(\tau ,x,y)}$

其中，Δs_rms(τ,x,y)为扰动场的慢度。

至此，完成对检波点一侧地震记录在成像点位置处τ时刻的波场延拓。流程进入到步骤107。

在步骤107，根据所述震源点到成像点的射线传播时间，利用激发时间成像条件，对频率-空间域地震波场作相位移动并进行有效频段的积分，得 到炮记录的地震成像结果；将炮记录偏移成像转换成共成像点道集，合并相同炮检距，该成像点道集可以用于后期AVO分析和岩性反演；将所有单炮成像结果按照空间位置叠加，或者对共成像点道集直接叠加，获得目标区的叠前时间偏移成像结果。

在一实施例中，对步骤106得到的频率-空间域地震波场延拓值，按照激发时间成像条件实施相位移动并对有效频段进行积分，得到最终的时间-空间域地震波场值，其成像条件表达式为：

$I(x,y,t_s)={\left(\frac{1}{2\pi }\right)}^2{\int }_{w1}^{w2}I(x,y,\tau ,w)e^{iw(\tau -t_s)}dw$

其中，I为相位移动后的时间-空间域地震波场，w1至w2为地震波的有效频带范围。

根据所述最终的时间-空间域地震波场成像值，将其所有单炮成像结果按照空间位置叠加，或者将炮记录偏移成像转换成共成像点道集，对共成像点道集直接叠加，，就能得到如图5所示的分步傅里叶叠前时间偏移成像结果，与图4所示采用射线法叠前时间偏移得到的marmousi偏移结果相比，采用本发明方法得到的偏移结果中波场信息丰富，地层内部弱反射清楚，体现了波场法叠前偏移的保幅性，同时在方框所示位置处，采用本发明方法得到的断层成像结果较射线法偏移有较大提高，断点干脆，且分辨率也较高。