CN104297789A - 一种三维倾角域稳相叠前时间偏移方法及系统 - Google Patents

Abstract

本发明公开了一种三维倾角域稳相叠前时间偏移方法及系统，包括：读取叠前地震资料，确定偏移速度场和拉伸切除参数；依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集；基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带；利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果；收集三维稳相叠前时间偏移的偏移结果，形成三维成像数据体；将三维成像数据体转换为地下反射构造的剖面图像。本发明有效地解决了存在陡倾角构造时现有的偏移成像技术得到的偏移图像信噪比低的问题。

Description

一种三维倾角域稳相叠前时间偏移方法及系统

技术领域

本发明涉及地震勘探技术领域，特别是涉及一种三维倾角域稳相叠前时间偏移方法及系统。

背景技术

地震勘探反射地震资料处理过程中，叠前偏移成像是关键的环节，而叠前时间偏移是叠前偏移成像中的一种重要方法。叠前时间偏移方法可对一类断层较为复杂但速度横向变化不是很剧烈的地质构造较好成像，与叠前深度偏移方法相比，除具有较高的计算效率外，其主要的优点是只需使用叠加(均方根)速度；这样可简单地通过速度扫描等方式得到恰当的速度模型，规避了使用叠前深度偏移方法面临的一个主要困难：速度建模。因此，叠前时间偏移方法已成为地震勘探领域广泛应用的重要技术。

影响叠前时间偏移成像效果的因素包括：偏移速度、地震波走时计算、偏移孔径、计算偏移幅值的权系数和偏移算法实现流程。走时计算与偏移速度共同决定了反射波能否正确归位，偏移孔径及其应用方式决定了偏移噪音和偏移算法的计算量，权系数决定了成像幅值能否正确反应地下界面的反射特征，偏移算法实现流程对偏移的计算效率和存储需求有重要影响。对偏移方法而言，成像效果、计算效率和存储需求是评价偏移方法的重要指标。现有偏移方法的上述指标与理想状态还有一定的差距。

偏移孔径对叠前时间偏移是重要的。较好的偏移孔径对压制偏移噪音和减少偏移计算量有重要作用。较小的偏移孔径可减少偏移计算量，但存在着不能对陡倾角构造正确成像的风险；过大的孔径又带来了偏移噪音和较大的偏移计算量。由于偏移算法实现方式的限制，在一个偏移作业中，总是选取一个统一的偏移孔径；这样，由于构造倾角是空变的，对某一成像点合适的孔径，对其他成像点就可能是过大或过小。此外，由于偏移前人们很难对拟成像构造的倾角有一个准确的估计，偏移孔径的选取只能采用保守的方法：选取较大的孔径。这就导致了较大的偏移噪音，特别是在存在陡倾角构造时。

最大程度地减少偏移噪音且又可对陡倾角构造正确成像，一直是叠前偏移技术的努力目标；而实现基于菲涅耳带叠加的偏移成像就是实现这一目标的努力方向。由于实际地质构造中速度和反射构造的复杂性，直接从地质模型出发估计准确的菲涅耳带几乎是不可能实现的。通过在偏移过程中构建倾角域偏移道集，即，将偏移结果根据拟成像反射界面的倾角大小进行分选、叠加，可直接将菲涅耳带形象地展示在这一偏移道集中，从而可容易地从这一道集中确定菲涅耳带，避免了直接估计准确菲涅耳带的巨大困难。得到准确的菲涅耳带，这可以在叠前时间偏移中计算中，通过计算拟成像反射界面的倾角，来判断该偏移结果是否参与叠加计算，从而实现基于菲涅耳带叠加的偏移成像。

发明内容

本发明要解决的技术问题是提供一种三维倾角域稳相叠前时间偏移方法及系统，用以提高偏移成像的信噪比和解决存在陡倾角构造时现有的偏移成像技术得到的偏移图像信噪比低的问题。

为解决上述技术问题，一方面，本发明提供一种三维倾角域稳相叠前时间偏移方法，包括：

读取叠前地震资料，确定偏移速度场和拉伸切除参数；

依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集；

基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带；

利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果；

收集三维稳相叠前时间偏移的偏移结果，形成三维成像数据体；

将三维成像数据体转换为地下反射构造的剖面图像。

进一步，读取叠前地震资料，确定偏移速度场和拉伸切除参数的步骤包括：

读取叠前地震资料，对叠前地震资料做线性干扰衰减；对完成了线性干扰衰减的叠前地震资料，选取部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到的动校正速度做横向平均，作为初始偏移速度；

将完成了线性干扰衰减的叠前地震资料按偏移距大小排序，基于给定的共反射点道集中成像道的偏移距间距分组，将不同组叠前地震资料存放到集群计算机的不同计算节点上，利用初始偏移速度，对已存放到集群计算机各个计算节点上的叠前地震资料，应用三维叠前时间偏移方法，进行并行的偏移计算，收集各计算节点的偏移结果，形成常规叠前时间偏移的共反射点道集；

对常规叠前时间偏移的共反射点道集，利用初始偏移速度做常规的反动校，再做常规的动校正得到新的速度，对新速度的倒数做空间平滑处理，空间平滑处理后得到的平滑速度场可作为偏移速度场；根据动校正后共反射点道集中大偏移距地震道的动校拉伸情况，确定拉伸切除参数。

进一步，在拟成像的区块上选定成像线，定义成像线是平行于测线的直线，CDP位置指示成像线上等间距的离散点；在选定的成像线上的每个CDP处，定义两个二维数组存放偏移结果，一个二维数组的一个维对应时间深度，另一个维是反射界面沿平行测线方向的倾角，该数组存储的偏移结果即是该CDP的平行测线方向的倾角道集；第二个二维数组的一个维对应时间深度，另一个维是反射界面沿垂直测线方向的倾角，该数组存储的偏移结果即是该CDP的垂直测线方向的倾角道集；通过偏移计算，在选定的成像线上每一个CDP处得到上述的两个数组，即是生成分别对应平行测线和垂直测线方向倾角的两个倾角道集。

进一步，通过偏移计算，在选定的成像线上每一个CDP处得到上述的两个数组，即是生成分别对应平行测线和垂直测线方向倾角的两个倾角道集的步骤包括：在偏移计算中，求取反射界面的倾角；对每个地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为从偏移速度场中拾取的成像点处的偏移速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量b＝x₁x₂+y₁y₂+2(V_rmsT₀)²，然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处的偏移结果所对应的反射界面的两个倾角θ_x和θ_y分别满足：

tanθ_x＝(x₁+g cosφ)/(V_rmsT₀)；

tanθ_y＝(y₁+g sinφ)/(V_rmsT₀)；

式中θ_x是反射界面沿平行测线方向的倾角，θ_y是反射界面沿垂直测线方向的倾角；

使用叠前地震资料，进行并行的偏移计算；对每一个选定的成像线，依据事先确定的偏移孔径，确定对该成像线有贡献的叠前地震资料；对这些叠前地震资料中的每一个地震道，根据偏移孔径和拉伸切除参数，确定该成像线在不同CDP处的起始成像的时间深度；利用偏移速度场，求得走时，再通过插值计算得到该地震道在成像点(x₀,y₀,T₀)处的偏移结果v；就这一偏移结果，x₀和y₀指示了CDP位置，T₀指示了在时间深度维上的位置，由计算得到的tanθ_x，可进一步确定在平行测线方向倾角维上的位置；这样，可将偏移结果v累加到第一个二维数组，即平行测线方向的倾角道集上；由计算得到的tanθ_y确定在垂直测线方向倾角维上的位置，再将偏移结果v累加到第二个二维数组，即垂直测线方向的倾角道集上；完成对全部有贡献的叠前地震资料中全部地震道的循环，收集各计算节点上的倾角道集并累加，即得到分别对应平行测线和垂直测线方向倾角的两个倾角道集。

进一步，基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带的步骤包括：对每一成像线的两个倾角道集，分别沿倾角方向进行边缘平滑衰减，在每个CDP处，分别沿两个倾角方向进行叠加，最后将两个倾角道集的叠加结果累加；这样可得到该成像线的偏移叠加剖面；

根据偏移叠加剖面上的断层和同相轴的变化情况，确定一组倾角拾取CDP位置；

对每个倾角拾取CDP位置，在偏移叠加剖面上拾取该CDP的成像道，做增益补偿，然后计算其幅值绝对值的平均值A₀；令a_i-1,a_i,a_i+1是成像道的连续3个时间样点上的幅值，若|a_i|≥0.2A₀，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i；对全部时间样点循环，可得到该CDP的波峰时间样点集；

利用波峰时间样点集，在该CDP的两个倾角道集中搜寻稳相点和菲涅耳带；分别对每个道集进行如下步骤，确定初步菲涅耳带，步骤如下：a)对倾角道集做带通滤波和增益补偿，在时间上取4个时窗，在每个时窗内计算倾角道集的幅值绝对值的平均值；b)对每个倾角的采样点，对波峰时间样点集循环，令i是时间样点的索引值，令a_i-1,a_i,a_i+1是该角度处在i左右连续3个时间样点上的幅值，令A是i所在的时窗内的幅值绝对值的平均值，若|a_i|≥0.2A，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i，得到该角度处的角度时间样点集；若角度时间样点集内样点间的时间采样点数小于2，剔除一个样点，保留幅值绝对值大的样点；c)令i和j为角度时间样点集内样点的时间和角度索引值，a_i,j为道集上该点的幅值；根据叠前地震资料的主频f₀，单位1/秒，用i和j对应的时间深度T_i和角度θ_j，角度θ_j的单位是弧度，可近似计算出倾角域菲涅耳带的两个方向的角度宽度m₁和m₂，即

$m_1=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j-\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

$m_2=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j+\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

式中Δθ为角度的采样间距，单位弧度；计算连续3个时间样点上的最大幅值b_ij，若a_i,j＞0，则b_ij＝max(a_i-1,j,a_i,j,a_i+1,j)，否则b_ij＝min(a_i-1,j,a_i,j,a_i+1,j)；引入两个无量纲参数α和β，一般取α＝0.3，β＝0.5；先计算菲涅耳带的平均幅值

$d=\frac{\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}\left|b_{\mathrm{ik}}\right|}{m_2-m_1+1}$

式中k代表近似菲涅耳带内的角度索引值，b_ik是角度k对应的最大幅值；若d≥αA，则继续计算菲涅耳带内的相关系数

$c=\frac{\left|\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}{b}_{\mathrm{ik}}\right|}{(m_2-m_1+1)d}$

c是无量纲数，若c≥β，则判断深度T_i和角度θ_j对应一个稳相点；对全部角度的角度时间样点集循环，记录得到的稳相点；d)对每个稳相点，重复m₁＝m₁-1和m₂＝m₂+1并计算d和c，当不能满足d≥αA和c≥β时，记录下对应的m₁和m₂；e)在整个时间深度上定义等间距的连续时窗，取50个样点为一个时窗；在每个时窗内，对全部的稳相点，取最小的m₁和其对应的稳相点的角度θ₁和时间深度T₁，取最大的m₂和其对应的稳相点的角度θ₂和时间深度T₂；角度θ₁-m₁·Δθ和T₁，角度θ₂+m₂·Δθ和T₂就分别是不同时间深度上菲涅耳带的下、上界；时窗内没有稳相点就忽略该时窗；f)对得到的各个菲涅耳带的下、上界分别沿时间深度线性插值，即得到该CDP点的初步菲涅耳带的界限曲线；

对倾角拾取CDP位置处的两个倾角道集，分别将拾取的初步菲涅耳带的上、下界曲线投影到倾角道集上，观察曲线是否正确保留了反射同相轴的菲涅耳带，对曲线做必要的修改，即得到该CDP的菲涅耳带；

分别对选定成像线的全部倾角拾取CDP位置的两个菲涅耳带的四组曲线进行空间插值，就可获得拟成像区块上全部成像点的菲涅耳带；将存储上、下界倾角对应的正切值。

进一步，利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果，包括：对叠前地震资料的全部地震道循环；对每一地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为成像点处的偏移速度，即均方根速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处对应的反射界面的两个倾角θ_x和θ_y的正切值为

tanθ_x＝(x₁+g cosφ)/(V_rmsT₀)

tanθ_y＝(y₁+g sinφ)/(V_rmsT₀)

读取成像点的菲涅耳带在两个倾角方向上的上、下界角度的正切值α₊，α_-，β₊和β_-，计算正切值的平均值a₀＝(α₊+α_-)/2和b₀＝(β₊+β_-)/2以及其变化的宽度a₁＝(α₊-α_-)/2和b₁＝(β₊-β_-)/2；计算菲涅耳带切除时的衰减带宽度和其中Δθ为倾角道集的角度采样间距，单位弧度；计算无量纲的衰减系数数组

$d_i=0.3+0.7\sin \left(\frac{(i-1)\mathrm{\π}}{100}\right),i=\mathrm{1,51}$

式中d_i是数组的元素，i是数组的索引数；

计算正切值的差x_A＝|tanθ_x-a₀|，若x_A＞a₁，不需近一步计算该成像点的偏移幅值；否则，进一步计算正切值的差y_A＝|tanθ_y-b₀|，判断是否满足y_A≤b₁；若满足，计算炮点和检波点到成像点的走时和成像权系数，由走时和成像权系数以及该地震道的地震记录确定偏移幅值；计算平行测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为c₁＝a₂/50和c₂＝a₁+0.1；计算正整数M＝2-int(x_A/(a₁-a₂))，进一步计算用于确定衰减系数的正整数

N₁＝int((a₁-x_A)/c_M)+1+50(M-1)

计算垂直测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为d₁＝b₂/50和d₂＝b₁+0.1；计算正整数M＝2-int(y_A/(b₁-b₂))，进一步计算用于确定衰减系数的正整数

N₂＝int((b₁-y_A)/d_M)+1+50(M-1)

依据N₁和N₂从衰减系数数组中提取对应的元素和求得无量纲的衰减系数将衰减系数η乘上前面得到的偏移幅值，将该结果累加到存放该成像点偏移结果的数组中对应的偏移距上，就得到仅对倾角域偏移道集中的菲涅耳带部分叠加成像的高信噪比偏移结果。

进一步，收集三维稳相叠前时间偏移在各计算节点上的偏移结果，得到三维稳相叠前时间偏移的共反射点道集；对得到的三维稳相叠前时间偏移的全部共反射点道集做剩余动校正；偏移和动校正过程将使得偏移距成像道出现拉伸和畸变，将拉伸超过50％的部分和畸变部分的幅值置为零，将不同偏移距的偏移结果叠加，形成三维成像数据体。

另一方面，本发明还提供一种三维倾角域稳相叠前时间偏移系统，其特征在于，包括：

参数准备单元，用于读取叠前地震资料，确定偏移速度场和拉伸切除参数；

倾角道集生成单元，用于依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集；

菲涅耳带生成单元，用于基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带；

偏移成像单元，用于利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果；

数据体形成单元，用于收集三维稳相叠前时间偏移的偏移结果，形成三维成像数据体；

数据体转换单元，用于将三维成像数据体转换为地下反射构造的剖面图像。

本发明有益效果如下：本发明通过生成倾角道集、在倾角道集上拾取菲涅耳带和通过插值获得全部成像点的菲涅耳带等手段，有效压制了偏移噪音和叠前资料中残存的噪音，可生成高信噪比的共反射点偏移道集和偏移叠加剖面；所生成的高信噪比共反射点偏移道集可更好地服务于叠前反演等油气和流体检测技术，高信噪比的偏移叠加剖面可更好地指示地下构造的形态、断裂以及地层沉积样式；该方法为提高地震成像的信噪比提供了一个有效途径。

附图说明

图1是本发明实施例涉及的一种三维倾角域稳相叠前时间偏移方法流程示意图；

图2是本发明实施例涉及的方法中偏移速度场在1085号成像线上的等值线图，图中数字是偏移速度值；

图3a是本发明实施例涉及的方法中1085号成像线在CDP为14.5km处沿测线方向倾角的倾角道集；图中横坐标是倾角的角度，单位度，从负角度到正角度；1085号成像线就是图2偏移速度场对应的成像线；

图3b是本发明实施例涉及的方法中1085号成像线在CDP为14.5km处垂直测线方向倾角的倾角道集；图中横坐标是倾角的角度，单位度，从负角度到正角度；1085号成像线就是图2偏移速度场对应的成像线；

图4a是本发明实施例涉及的方法中图3a的倾角道集对应的菲涅耳带，图中的一对实线折线指示了自动确定的菲涅耳带的角度上、下限；

图4b是本发明实施例涉及的方法中图3b的倾角道集对应的菲涅耳带，图中的一对实线折线指示了自动确定的菲涅耳带的角度上、下限；

图5是本发明实施例涉及的方法中1085号成像线的沿测线方向倾角的最终菲涅耳带上界的等值线图，图5中数字是倾角上界对应的正切值；

图6是本发明实施例涉及的方法中三维倾角域稳相叠前时间偏移方法得到的三维偏移数据体在1085号成像线上的剖面图；

图7a是本发明实施例涉及的方法中常规叠前时间偏移方法得到的1085号成像线上的局部成像结果；

图7b是本发明实施例涉及的方法中图6的三维倾角域稳相叠前时间偏移方法的偏移结果在图7a相同局部区域的放大图；

图8是本发明实施例涉及的一种三维倾角域稳相叠前时间偏移系统结构示意图。

具体实施方式

为了提高偏移成像的信噪比和解决存在陡倾角构造时现有的偏移成像技术得到的偏移图像信噪比低的问题，本发明提供了一种三维倾角域稳相叠前时间偏移方法及系统，以下结合附图以及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不限定本发明。

如图1所示，本发明实施例中的一种三维倾角域稳相叠前时间偏移方法，包括：

S101，读取叠前地震资料，确定偏移速度场和拉伸切除参数；

S102，依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集；

S103，基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带；

S104，利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果；

S105，收集三维稳相叠前时间偏移的偏移结果，形成三维成像数据体；

S106，将三维成像数据体转换为地下反射构造的剖面图像

下面以以东部陆上油田为例，对本发明实施例涉及的三维倾角域稳相叠前时间偏移方法加以详细描述，包括：

(1)用多条测线记录人工震源激发产生的地震波信号，记录到磁带上，得到叠前地震资料。具体采集参数是，12条测线同时记录地震信号，每条测线有144个检波器组，检波器组间的间距为50m；沿垂直测线方向布炮线，炮间距50m，炮线间的线间距200m，共采集39240炮，记录时长5s，时间采样1ms。

(2)从磁带上读取叠前地震资料，对叠前地震资料做线性干扰衰减；对完成了线性干扰衰减的叠前地震资料，选取部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到的动校正速度做横向平均，作为初始偏移速度。

(3)将完成了线性干扰衰减的叠前地震资料按偏移距大小排序，基于给定的共反射点道集中成像道的偏移距间距分组，将不同组叠前地震资料存放到集群计算机的不同计算节点上。具体是，将叠前地震资料按偏移距大小分为76组，放置到集群计算机的19个节点上，每个计算节点放置4组叠前地震资料。

(4)利用初始偏移速度，对已存放到集群计算机各个计算节点上的叠前地震资料，应用三维叠前时间偏移方法，进行并行的偏移计算，收集各计算节点的偏移结果，形成常规叠前时间偏移的共反射点道集。

(5)对常规叠前时间偏移的共反射点道集，利用初始偏移速度做常规的反动校，再做常规的动校正得到新的速度，对新速度的倒数做空间平滑处理，空间平滑处理后得到的平滑速度场可作为偏移速度场；根据动校正后共反射点道集中大偏移距地震道的动校拉伸情况，确定拉伸切除参数。图2是偏移速度场在1085号成像线上的等值线图，1085号成像线距成像区域边界5.35km；具体的拉伸切除参数为15ms处188m,437ms处826m,845ms处1608m,1932ms处3633m,2000ms以下全保留。

(6)在拟成像的区块上选定成像线，依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集。具体是，在拟成像的区块上，选定425、475、525、565、625、675、725、779、825、875、925、975、1025、1085、1130号共15条成像线，最大和最小号成像线距成像区域边界的距离分别是3.1km和4.225km。图3a和图3b给出了1085号成像线在CDP为14.5km处的两个倾角道集。

在拟成像区块的水平表面定义二维坐标，令平行测线方向为x轴，垂直测线方向为y轴，则成像线为y坐标是常量的直线，即是平行于测线的直线，成像线按y坐标大小顺序给予编号；CDP(Common Depth Point，共深度点)位置指示成像线上等间距的离散点，其大小代表x坐标的变化；在选定的成像线上的每个CDP处，定义两个二维数组存放偏移结果，一个二维数组的一个维对应时间深度，即单程旅行时，另一个维是反射界面沿平行测线方向的倾角，该数组存储的偏移结果即是该CDP的平行测线方向的倾角道集；第二个二维数组的一个维对应时间深度，另一个维是反射界面沿垂直测线方向的倾角，该数组存储的偏移结果即是该CDP的垂直测线方向的倾角道集；通过偏移计算，在选定的成像线上每一个CDP处得到上述的两个数组，即是生成了对应平行测线和垂直测线方向倾角的两个倾角道集；

获得倾角道集的关键是在偏移计算中，利用炮点、检波点和成像点坐标和偏移速度场，求得成像点处反射界面的两个倾角，其算法如下：对每个地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为从偏移速度场中拾取的成像点处的偏移速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量b＝x₁x₂+y₁y₂+2(V_rmsT₀)²，然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处的偏移结果所对应的反射界面的两个倾角θ_x和θ_y分别满足：

tanθ_x＝(x₁+g cosφ)/(V_rmsT₀)；

tanθ_y＝(y₁+g sinφ)/(V_rmsT₀)；

式中θ_x是反射界面沿平行测线方向的倾角，θ_y是反射界面沿垂直测线方向的倾角；

使用已存放到集群计算机各个计算节点上的叠前地震资料，进行并行的偏移计算；对每一个选定的成像线，依据事先确定的偏移孔径，确定对该成像线有贡献的叠前地震资料；对这些叠前地震资料中的每一个地震道，根据偏移孔径和拉伸切除参数，确定在不同CDP处的起始成像的时间深度；利用偏移速度场，求得走时，再通过插值计算得到该地震道在成像点(x₀,y₀,T₀)处的偏移结果v；就这一偏移结果，x₀和y₀指示了CDP位置，T₀指示了在时间深度维上的位置，由计算得到的tanθ_x，可进一步确定在平行测线方向倾角维上的位置；这样，可将偏移结果v累加到第一个二维数组，即平行测线方向的倾角道集上；由计算得到的tanθ_y确定在垂直测线方向倾角维上的位置，可再将偏移结果v累加到第二个二维数组，即垂直测线方向的倾角道集上；完成对全部有贡献的叠前地震资料中全部地震道的循环，收集各计算节点上的倾角道集并累加，即可得到选定成像线的平行测线和垂直测线方向倾角的两个倾角道集。

(7)根据步骤6得到的选定成像线的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带。具体是，取Δθ＝4度，α＝0.3，β＝0.5；图4a和图4b中的道集就是图3a和图3b中的倾角道集，图中的一对实线折线代表了自动确定的菲涅耳带的角度上、下限；图5给出了最终确定的1085号成像线平行测线方向倾角的菲涅耳带上界的等值线图。

采用以下6个步骤：1)针对选定的成像线中的每个成像线，对得到的两个倾角道集进行如下处理：分别沿倾角方向进行边缘平滑衰减，在每个CDP处，分别沿两个倾角方向进行叠加，最后将两个倾角道集的叠加结果累加；这样可得到该成像线的偏移叠加剖面；2)根据偏移叠加剖面上的断层和同相轴的变化情况，从可利用线性插值确定全部成像点的菲涅耳带的角度出发，确定一组倾角拾取CDP位置；3)对每个倾角拾取CDP位置，在偏移叠加剖面上拾取该CDP的成像道，做增益补偿，然后计算其幅值绝对值的平均值A₀；令a_i-1,a_i,a_i+1是成像道的连续3个时间样点上的幅值，若|a_i|≥0.2A₀，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i；对全部时间样点循环，可得到该CDP的波峰时间样点集；

4)利用波峰时间样点集，在该CDP的两个倾角道集中搜寻稳相点和菲涅耳带；分别对每个道集进行如下步骤，自动确定初步菲涅耳带，步骤如下：a)对倾角道集做带通滤波和增益补偿，在时间上取4个时窗，在每个时窗内计算倾角道集的幅值绝对值的平均值；b)对每个倾角的采样点，对波峰时间样点集循环，令i是时间样点的索引值，令a_i-1,a_i,a_i+1是该角度处在i左右连续3个时间样点上的幅值，令A是i所在的时窗内的幅值绝对值的平均值，若|a_i|≥0.2A，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i，得到该角度处的角度时间样点集；若角度时间样点集内样点间的时间采样点数小于2，剔除一个样点，保留幅值绝对值大的样点；c)令i和j为角度时间样点集内样点的时间和角度索引值，a_i,j为道集上该点的幅值；根据地震资料的主频f₀，单位1/秒,可用i和j对应的时间深度T_i和角度θ_j，角度θ_j的单位是弧度，近似计算出倾角域菲涅耳带的两个方向的角度宽度m₁和m₂，即

$m_1=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j-\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

$m_2=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j+\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

式中Δθ为角度的采样间距，单位弧度；计算连续3个时间样点上的最大幅值b_ij，若a_i,j＞0，则b_ij＝max(a_i-1,j,a_i,j,a_i+1,j)，否则b_ij＝min(a_i-1,j,a_i,j,a_i+1,j)；引入两个无量纲参数α和β，一般取α＝0.3，β＝0.5；先计算菲涅耳带的平均幅值

$d=\frac{\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}\left|b_{\mathrm{ik}}\right|}{m_2-m_1+1}$

式中k代表近似菲涅耳带内的角度索引值，b_ik是角度k对应的最大幅值；若d≥αA，则继续计算菲涅耳带内的相关系数

$c=\frac{\left|\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}{b}_{\mathrm{ik}}\right|}{(m_2-m_1+1)d}$

c是无量纲数，若c≥β，则可判断深度T_i和角度θ_j对应一个稳相点；对全部角度的角度时间样点集循环，记录得到的稳相点；d)对每个稳相点，重复m₁＝m₁-1和m₂＝m₂+1并计算d和c，当不能满足d≥αA和c≥β时，记录下对应的m₁和m₂；e)在整个时间深度上定义等间距的连续时窗，一般取50个样点为一个时窗；在每个时窗内，对全部的稳相点，取最小的m₁和其对应的稳相点的角度θ₁和时间深度T₁，取最大的m₂和其对应的稳相点的角度θ₂和时间深度T₂；角度θ₁-m₁·Δθ和T₁，角度θ₂+m₂·Δθ和T₂就分别是不同时间深度上菲涅耳带的下、上界；时窗内没有稳相点就忽略该时窗；f)对得到的各个菲涅耳带的下、上界分别沿时间深度线性插值，即得到该CDP点的初步菲涅耳带的界限曲线；

5)对倾角拾取CDP位置处的两个倾角道集，分别将自动拾取的初步菲涅耳带的上、下界曲线投影到倾角道集上，观察曲线是否正确保留了反射同相轴的菲涅耳带，对曲线做必要的修改，即得到该CDP的菲涅耳带；6)分别对选定成像线的全部倾角拾取CDP位置的两个菲涅耳带的四组曲线进行空间插值，就可得到拟成像区块上全部成像点的菲涅耳带；将存储上、下界倾角对应的正切值。

(8)用步骤7得到的菲涅耳带，进行三维稳相叠前时间偏移，在各计算节点上得到对应不同偏移距的偏移结果。

在每个计算节点，对叠前地震资料的全部地震道循环；对每一地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为成像点处的偏移速度，即均方根速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处对应的反射界面的两个倾角θ_x和θ_y的正切值为

tanθ_x＝(x₁+g cosφ)/(V_rmsT₀)

tanθ_y＝(y₁+g sinφ)/(V_rmsT₀)

读取成像点的菲涅耳带在两个倾角方向上的上、下界角度的正切值α₊，α_-，β₊和β_-，计算正切值的平均值a₀＝(α₊+α_-)/2和b₀＝(β₊+β_-)/2以及其变化的宽度a₁＝(α₊-α_-)/2和b₁＝(β₊-β_-)/2；计算菲涅耳带切除时的衰减带宽度和其中Δθ为倾角道集的角度采样间距，单位弧度；计算无量纲的衰减系数数组

$d_i=0.3+0.7\sin \left(\frac{(i-1)\mathrm{\π}}{100}\right),i=\mathrm{1,51}$

式中d_i是数组的元素，i是数组的索引数；

计算正切值的差x_A＝|tanθ_x-a₀|，若x_A＞a₁，不需近一步计算该成像点的偏移幅值；否则，进一步计算正切值的差y_A＝|tanθ_y-b₀|，判断是否满足y_A≤b₁；若满足，计算炮点和检波点到成像点的走时和成像权系数，由走时和成像权系数以及该地震道的地震记录确定偏移幅值；计算平行测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为c₁＝a₂/50和c₂＝a₁+0.1；计算正整数M＝2-int(x_A/(a₁-a₂))，进一步计算用于确定平行测线方向衰减系数的正整数

N₁＝int((a₁-x_A)/c_M)+1+50(M-1)

下面计算垂直测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为d₁＝b₂/50和d₂＝b₁+0.1；计算正整数M＝2-int(y_A/(b₁-b₂))，进一步计算用于确定垂直测线方向衰减系数的正整数

N₂＝int((b₁-y_A)/d_M)+1+50(M-1)

依据N₁和N₂从衰减系数数组中提取对应的元素和可求得无量纲的衰减系数将衰减系数η乘上前面得到的偏移幅值，将该结果累加到存放该成像点偏移结果的数组中对应的偏移距上，就可得到仅对倾角域偏移道集中的菲涅耳带部分叠加成像的高信噪比偏移结果。

(9)收集三维稳相叠前时间偏移在各计算节点上的偏移结果，得到三维稳相叠前时间偏移的共反射点道集；对得到的三维稳相叠前时间偏移的全部共反射点道集做剩余动校正；偏移和动校正过程将使得大偏移距成像道出现拉伸和畸变，将拉伸超过50％的部分和畸变部分的幅值置为零，将不同偏移距的偏移结果叠加，形成三维成像数据体。

(10)通过显示软件将三维成像数据体转换为地下反射构造的剖面图像，剖面图像将指示地下构造的形态、断裂部位、断距大小和地层沉积样式及地层的波阻抗特征，用于确定地下生、储油构造和识别油气储。如图6、图7a和图7b所示，应用本发明实施例涉及的方法后成像的信噪比明显提高，不清楚的同相轴和断层更加清晰。

为了更好的理解本发明实施例涉及的方法，简述本发明实施例的具体实现原理：

1.在叠前时间偏移的过程中直接生成倾角域偏移道集

将三维非均匀介质近似为层状介质。在波数－频率域，单个检波点记录的地震波场的深度延拓可用旅行时(时间深度)表示为：

$\tilde{P}(p_x,p_y,\mathrm{\ω},T=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{T}_i)=F\left(\mathrm{\ω}\right)\exp (-\mathrm{j\ω}(p_x{x}_g+p_y{y}_g))\exp \left(\mathrm{j\ω}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i\sqrt{1-v_i^2(p_x^2+p_y^2)}\right)---\left(1\right)$

式中ΔT_i是各层介质用单程旅行时表达的厚度，n是某一深度包含的层数，是该深度对应的时间深度(单程旅行时)，v_i是各层介质的速度，p_x和p_y是射线参数，j是单位虚数，x_g和y_g是检波点的坐标，ω是频率，F(ω)是检波点记录的时域信号的傅立叶变换。

式(1)中第二个指数函数的指数项中的累加部分可近似表达为：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left({\mathrm{\ΔT}}_i\sqrt{1-v_i^2(p_x^2+p_y^2)}\right)\≈\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}\·\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i\right)---\left(2\right)$

对(2)式两边做Taylor展开，并截取到第二项，可得：

$V_{\mathrm{rms}}^2=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_i^2{\mathrm{\ΔT}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔT}}_i}---\left(3\right)$

V_rms即是叠前时间偏移所需要的偏移速度，即均方根速度。

将(2)式代入(1)式并做空间傅立叶反变换，得空间－频率域的波场为：

$P(x,y,\mathrm{\ω},T)=\frac{{\mathrm{\ω}}^2}{4{\mathrm{\π}}^2}\∫\∫F\left(\mathrm{\ω}\right)\exp \left[\mathrm{j\ω}(T\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}+p_x(x-x_g)+p_y(y-y_g))\right]{\mathrm{dp}}_x{\mathrm{dp}}_y---\left(4\right)$

根据稳相点原理，式(4)中积分的主要贡献来自射线参数的稳相点而这一稳相点就是检波点反传波场的传播方向。若定义

$\mathrm{\φ}(p_x,p_y)=T\sqrt{1-V_{\mathrm{rms}}^2(p_x^2+p_y^2)}+p_x(x-x_g)+p_y(y-y_g)---\left(5\right)$

可由

$\frac{\∂\mathrm{\φ}(p_x,p_y)}{{\∂p}_x}=0,\frac{\∂\mathrm{\φ}(p_x,p_y)}{{\∂p}_y}=0$

求得稳相点为：

$p_y^0=\frac{y-y_g}{V_{\mathrm{rms}}\sqrt{V_{\mathrm{rms}}^2{T}^2+{(x-x_g)}^2+{(y-y_g)}^2}}---\left(7\right)$

式中(x,y,T)即是成像点的坐标。若在成像点处忽略介质层速度与偏移速度V_rms的不同，可得成像点处检波点反传波场，即反射波场的方向余弦为

$l_g={\mathrm{vp}}_x^0=\frac{(x-x_g)}{\sqrt{{(x-x_g)}^2+{(y-y_g)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(8\right)$

$m_g={\mathrm{vp}}_y^0=\frac{(y-y_g)}{\sqrt{{(x-x_g)}^2+{(y-y_g)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(9\right)$

$n_g=\sqrt{1-l_g^2-m_g^2}=\frac{V_{\mathrm{rms}}T}{\sqrt{{(x-x_g)}^2+{(y-y_g)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(10\right)$

同理，可得成像点处炮点入射波场的方向余弦为

$l_s=\frac{(x-x_s)}{\sqrt{{(x-x_s)}^2+{(y-y_s)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(11\right)$

$m_s=\frac{(y-y_s)}{\sqrt{{(x-x_s)}^2+{(y-y_s)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(12\right)$

$n_s=\frac{V_{\mathrm{rms}}T}{\sqrt{{(x-x_s)}^2+{(y-y_s)}^2+V_{\mathrm{rms}}^2{T}^2}}---\left(13\right)$

定义成像点处反射界面法线的方向余弦为(l₀,m₀,n₀)，由于入射波场和反射波场是在这一界面发生反射，因此三组方向余弦对应的矢量满足：1)三个矢量共面，2)两者间的夹角相等。由第一个条件可得

由第二个条件可得

l_sl₀+m_sm₀+n_sn₀＝l_gl₀+m_gm₀+n_gn₀     (15)

再考虑到方向余弦的模为1，可由式(14)和(15)解得界面法线的方向余弦(l₀,m₀,n₀)。令反射界面在沿测线方向的倾角为θ_x，则矢量(cosθ_x,0,sinθ_x)恰落在该反射界面中，将该矢量与界面法线的方向余弦点乘，有

l₀cosθ_x+n₀sinθ_x＝0

可解得

tanθ_x＝-l₀/n₀     (16)

同理可解得垂直测线方向的倾角θ_y为

tanθ_x＝-m₀/n₀     (17)

为避免炮点、检波点、成像点间水平距离太近时式(16)和(17)中除法带来的不稳定，应在公式推导中约去这些距离量。定义如下距离变量：x₁＝x_s-x，y₁＝y_s-y，x₂＝x_s+x_r-2x，y₂＝y_s+y_r-2y；定义φ是方位角；计算a＝x₁cosφ+y₁sinφ，b＝x₁x₂+y₁y₂+2(V_rmsT)²，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}T\right)}^2}+\sqrt{{(x_r-x)}^2+{(y_r-y)}^2+{\left(V_{\mathrm{rms}}T\right)}^2},$ 然后计算g＝(a·d²-b·c)/(c²-d²)；则式(16)和(17)可表达为

tanθ_x＝(x₁+g cosφ)/(V_rmsT)     (18)

tanθ_y＝(y₁+g sinφ)/(V_rmsT)     (19)

若在每一成像点，不考虑偏移距，仅将偏移结果按θ_x和θ_y的大小进行分选和叠加，就可在每个水平位置(即CDP点)形成两个分别针对反射界面沿测线方向倾角和垂直测线方向倾角的倾角域偏移道集。

2.由倾角道集确定菲涅耳带

在倾角域偏移道集中，反射同相轴呈现双曲线形态，而这双曲线的顶点就对应反射界面的真实倾角。由于常规偏移结果是倾角域偏移道集沿倾角方向的叠加，从叠加的角度来看，这个顶点就是稳相点，而顶点的邻域就是菲涅耳带。在倾角域偏移道集中，菲涅耳带被形象地展示，因此可以容易从这一道集中确定菲涅耳带，这避免了直接估计准确菲涅耳带的巨大困难。

菲涅耳带实际上对应双曲反射同相轴的近似水平的部分。令a_i,j为时间样点i和角度样点j处道集的幅值，b_ij＝min(a_i-1,j,a_i,j,a_i+1,j)，可用如下指标c的大小来判断i和j是否是稳相点，即过该点的同相轴是否是水平的，定义

$c=\frac{\left|\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}{b}_{\mathrm{ik}}\right|}{\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}\left|b_{\mathrm{ik}}\right|}---\left(20\right)$

式中m₁和m₂代表近似的菲涅耳带的角度下、上界，c等于1表明同相轴是绝对水平的；m₁和m₂的选择很重要，若两个值太接近，在许多样点处都将有接近1的c，若两个值相差太大，真正的稳相点也不能满足条件。

本发明通过用零偏移距资料近似估计倾角域的菲涅耳带，来确定m₁和m₂。其核心思路是，在反射界面上定义稳相点的边缘成像点，该成像点处法向出射的地震射线，其到接收点的走时应与该接收点到稳相点上方时间深度差四分之一周期的成像点处相同；利用稳相点处已知的反射界面的两个倾角和走时条件，可在反射界面上求得菲涅耳带的边缘曲线。在反射界面的平面上，该曲线是一半径为的圆，式中f₀是地震资料的主频，T是稳相点的时间深度。

计算菲涅耳带的边缘曲线上的反射波在稳相点上方成像时对应的两个倾角，就可得到用两个倾角表达的菲涅耳带边缘曲线。由于本发明用的两个倾角道集实际是有关两个倾角的二维倾角道集分别沿对应倾角方向叠加的结果，对沿测线方向倾角道集而言，其角度下、上界是菲涅耳带边缘曲线对应的沿测线方向倾角的最小、最大值；而垂直测线方向倾角道集的角度下、上界则是菲涅耳带边缘曲线对应的垂直测线方向倾角的最小、最大值。对两个倾角道集，有统一的m₁和m₂的表达式：

$\begin{array}{c}{m}_1=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j-\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}\\ m_2=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j+\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}\end{array}---\left(21\right)$

式中T_i是样点的时间深度，Δθ为角度样点的采样间距，对沿测线方向倾角道集，θ_j是沿测线方向倾角，而对垂直测线方向倾角道集，θ_j则是垂直测线方向倾角。

式(20)的判断指标对孤立的大值点等特殊情况也会得到较大的值，因此应首先排除那些明显不是稳相点的样点；为此本发明采用平均值判断标准，先计算平均值

$d=\frac{\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}\left|b_{\mathrm{ik}}\right|}{m_2-m_1+1}---\left(22\right)$

当平均值大于某一给定值时才进行式(20)的计算。

逐样点进行式(20)和(22)的计算是很不经济的，实际上也没必要逐点确定稳相点。为此，我们进行如下的前期处理，以筛选确定一批可能的样点，具体步骤如下：1)对每个倾角拾取CDP点，在偏移叠加剖面上拾取该CDP点处的成像道，做增益补偿，然后计算其幅值绝对值的平均值A₀；令a_i-1,a_i,a_i+1是成像道的连续3个时间样点上的幅值，若|a_i|≥0.2A₀，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i；对全部时间样点循环，可得到一个该CDP点的峰值时间样点集；2)对倾角道集做带通滤波和增益补偿，在时间上取4个时窗，在每个时窗内计算道集的幅值绝对值的平均值；3)在每个角度采样点，对峰值时间样点集循环，令i是时间样点，令a_i-1,a_i,a_i+1是道集在i左右连续3个时间样点上的幅值，令A是i所在的时窗内的幅值绝对值的平均值，若|a_i|≥0.2A，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i，得到该角度的时间样点集；若时间样点集内样点间的时间采样点数小于2，剔除一个样点，保留幅值绝对值大的样点。

式(20)确定的仅是稳相点，还需近一步准确决定菲涅耳带的下、上界。本发明仍采用式(20)和(22)来确定菲涅耳带的下、上界，方法如下：对每个稳相点，重复m₁＝m₁-1和m₂＝m₂+1并由式(20)和(22)计算d和c，当不能满足判断标准时，记录下对应的m₁和m₂；在整个时间深度上定义等间距的连续时窗，在每个时窗内，对全部稳相点，取最小的m₁和其对应的稳相点的角度θ₁和时间深度T₁，取最大的m₂和其对应的稳相点的角度θ₂和时间深度T₂，而角度θ₁-m₁·Δθ和T₁，角度θ₂+m₂·Δθ和T₂就分别是不同时间深度上菲涅耳带的下、上界。对得到的各个菲涅耳带的下、上界分别沿时间深度线性插值，即得到该CDP点的菲涅耳带曲线。由于倾角道集上已直观地展示了菲涅耳带，可进一步将插值得到的菲涅耳带下、上界的曲线投影到倾角道集上，以此修正自动拾取产生的错误。

3.三维稳相叠前时间偏移

利用得到的菲涅耳带，这可在叠前时间偏移计算中，通过计算拟成像反射界面的倾角，来判断该偏移结果是否应累加到相应的道集中。这样，就可保证仅是菲涅耳带内的结果参与偏移成像。如此，即可避免常规偏移中采用较小的偏移孔径导致的陡倾角构造不能正确成像，又可避免采用过大的偏移孔径带来的偏移噪音问题。

与常规叠前时间偏移相比，稳相叠前时间偏移在对每一地震道的偏移计算中，先利用式(18)和(19)计算在每一成像点处的正切值tanθ_x和tanθ_y，同时读取该成像点的菲涅耳带在两个倾角方向上的上、下界角度的正切值α₊，α_-，β₊和β_-；计算a₀＝(α₊+α_-)/2，a₁＝(α₊-α_-)/2，b₀＝(β₊+β_-)/2和b₁＝(β₊-β_-)/2。利用近似关系

tan(θ+Δθ)≈tanθ+Δθ(1+tan²θ)     (23)

可计算菲涅耳带切除时对应角度宽度Δθ的衰减带宽度；因此，有衰减带宽度和其中Δθ为倾角道集的角度采样间距，单位弧度。事先计算衰减系数数组

$d_i=0.3+0.7\sin \left(\frac{(i-1)\mathrm{\π}}{100}\right),i=\mathrm{1,51}---\left(24\right)$

计算x_A＝|tanθ_x-a₀|，若x_A＞a₁，不对该成像点的偏移幅值做任何操作；否则，进一步计算y_A＝|tanθ_y-b₀|，判断是否满足y_A≤b₁；若满足，计算炮点和检波点到成像点的走时和成像权系数，由走时和成像权系数以及该地震道的幅值确定偏移幅值。然后，计算沿测线方向角度的幅值衰减系数；定义一个有两个元素的数组，c₁＝a₂/50，c₂＝a₁+0.1；计算正整数M＝2-int(x_A/(a₁-a₂))，进一步计算衰减系数索引数

N₁＝int((a₁-x_A)/c_M)+1+50(M-1)     (25)

当x_A≤a₁-a₂，即正切值没落在衰减带内，由(25)式可得N₁＝51，此时由式(24)可知衰减系数为1，即不衰减；当x_A＞a₁-a₂，即正切值落在衰减带内，由(25)式可得N₁≤50，此时可由式(24)取得小于1的正数。

同理，可得到垂直测线方向角度的幅值衰减系数，将两个衰减系数相乘，就可得到总的衰减系数。将总的衰减系数乘上前面得到的偏移幅值，将该结果累加到存放该成像点偏移结果的数组中对应的偏移距上，就可得到仅对倾角域偏移道集中的菲涅耳带部分叠加成像的高信噪比偏移结果。

本发明实施例涉及的方法，能大幅压制偏移噪音，得到更高信噪比的共反射点偏移道集和偏移图像；通过计算拟成像反射界面的倾角，来判断该偏移结果是否参与叠加计算，从而实现了基于菲涅耳带叠加的偏移成像。

本发明实施例涉及的方法得到的倾角域偏移道集是将偏移结果根据拟成像反射界面的倾角大小进行分选、叠加；在倾角域偏移道集中，反射界面的真实倾角对应着弯曲同相轴的顶点，由于常规偏移实际是倾角域偏移道集沿倾角方向的叠加，从叠加的角度来看，这个顶点就是稳相点，而顶点的邻域就是菲涅耳带。由于菲涅耳带已被形象地展示在倾角域偏移道集中，就可以容易地从这一道集中确定菲涅耳带，这避免了直接估计准确菲涅耳带的巨大困难。

得到准确的菲涅耳带，这可以在叠前时间偏移中计算中，通过计算拟成像反射界面的倾角，来判断该偏移结果是否参与叠加计算，从而实现基于菲涅耳带的偏移成像。这样，就可以避免常规偏移孔径出现的问题，获得更高信噪比的偏移结果，特别是当地质构造存在陡倾角构造时。该方法也可理解为采用了空变的最优偏移孔径的叠前时间偏移方法。

如图8所示，本发明实施例中的一种三维倾角域稳相叠前时间偏移系统，包括：参数准备单元801，用于读取叠前地震资料，确定偏移速度场和拉伸切除参数；

倾角道集生成单元802，用于依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集；

菲涅耳带生成单元803，用于基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带；

偏移成像单元804，用于利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果；

数据体形成单元805，用于收集三维稳相叠前时间偏移的偏移结果，形成三维成像数据体；

数据体转换单元806，用于将三维成像数据体转换为地下反射构造的剖面图像。

参数准备单元，还用于读取叠前地震资料，对叠前地震资料做线性干扰衰减；对完成了线性干扰衰减的叠前地震资料，选取部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到的动校正速度做横向平均，作为初始偏移速度；

将完成了线性干扰衰减的叠前地震资料按偏移距大小排序，基于给定的共反射点道集中成像道的偏移距间距分组，将不同组叠前地震资料存放到集群计算机的不同计算节点上，利用初始偏移速度，对已存放到集群计算机各个计算节点上的叠前地震资料，应用三维叠前时间偏移方法，进行并行的偏移计算，收集各计算节点的偏移结果，形成常规叠前时间偏移的共反射点道集；

对常规叠前时间偏移的共反射点道集，利用初始偏移速度做常规的反动校，再做常规的动校正得到新的速度，对新速度的倒数做空间平滑处理，空间平滑处理后得到的平滑速度场可作为偏移速度场；根据动校正后共反射点道集中大偏移距地震道的动校拉伸情况，确定拉伸切除参数；

倾角道集生成单元，还用于在拟成像的区块上选定成像线，定义成像线是平行于测线的直线，CDP位置指示成像线上等间距的离散点；在选定的成像线上的每个CDP处，定义两个二维数组存放偏移结果，一个二维数组的一个维对应时间深度，另一个维是反射界面沿平行测线方向的倾角，该数组存储的偏移结果即是该CDP的平行测线方向的倾角道集；第二个二维数组的一个维对应时间深度，另一个维是反射界面沿垂直测线方向的倾角，该数组存储的偏移结果即是该CDP的垂直测线方向的倾角道集；通过偏移计算，在选定的成像线上每一个CDP处得到上述的两个数组；

在偏移计算中，求取反射界面的倾角；对每个地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为从偏移速度场中拾取的成像点处的偏移速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处的偏移结果所对应的反射界面的两个倾角θ_x和θ_y分别满足：

tanθ_x＝(x₁+g cosφ)/(V_rmsT₀)；

tanθ_y＝(y₁+g sinφ)/(V_rmsT₀)；

式中θ_x是反射界面沿平行测线方向的倾角，θ_y是反射界面沿垂直测线方向的倾角；

使用叠前地震资料，进行并行的偏移计算；对每一个选定的成像线，依据事先确定的偏移孔径，确定对该成像线有贡献的叠前地震资料；获对这些叠前地震资料中的每一个地震道，根据偏移孔径和拉伸切除参数，确定在不同CDP处的起始成像的时间深度；利用偏移速度场，求得走时，再通过插值计算得到该地震道在成像点(x₀,y₀,T₀)处的偏移结果v；就这一偏移结果，x₀和y₀指示了CDP位置，T₀指示了在时间深度维上的位置，由计算得到的tanθ_x，进一步确定在平行测线方向倾角维上的位置；这样，可将偏移结果v累加到第一个二维数组，即平行测线方向的倾角道集上；由计算得到的tanθ_y确定在垂直测线方向倾角维上的位置，再将偏移结果v累加到第二个二维数组，即垂直测线方向的倾角道集上；完成对全部有贡献的叠前地震资料中全部地震道的循环，收集各计算节点上的倾角道集并累加，即得到分别对应平行测线和垂直测线方向倾角的两个倾角道集；

菲涅耳带生成单元，还用于对每一成像线的两个倾角道集，分别沿倾角方向进行边缘平滑衰减，在每个CDP处，分别沿两个倾角方向进行叠加，最后将两个倾角道集的叠加结果累加；这样可得到该成像线的偏移叠加剖面；

根据偏移叠加剖面上的断层和同相轴的变化情况，确定一组倾角拾取CDP位置；

对每个倾角拾取CDP位置，在偏移叠加剖面上拾取该CDP的成像道，做增益补偿，然后计算其幅值绝对值的平均值A₀；令a_i-1,a_i,a_i+1是成像道的连续3个时间样点上的幅值，若|a_i|≥0.2A₀，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i；对全部时间样点循环，可得到该CDP的波峰时间样点集；

利用波峰时间样点集，在该CDP的两个倾角道集中搜寻稳相点和菲涅耳带；分别对每个道集进行如下步骤，确定初步菲涅耳带，步骤如下：a)对倾角道集做带通滤波和增益补偿，在时间上取4个时窗，在每个时窗内计算倾角道集的幅值绝对值的平均值；b)对每个倾角的采样点，对波峰时间样点集循环，令i是时间样点的索引值，令a_i-1,a_i,a_i+1是该角度处在i左右连续3个时间样点上的幅值，令A是i所在的时窗内的幅值绝对值的平均值，若|a_i|≥0.2A，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i，得到该角度处的角度时间样点集；若角度时间样点集内样点间的时间采样点数小于2，剔除一个样点，保留幅值绝对值大的样点；c)令i和j为角度时间样点集内样点的时间和角度索引值，a_i,j为道集上该点的幅值，根据叠前地震资料的主频f₀，单位1/秒，用i和j对应的时间深度T_i和角度θ_j，角度θ_j的单位是弧度，近似计算出倾角域菲涅耳带的两个方向的角度宽度m₁和m₂，即

$m_1=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j-\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

$m_2=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j+\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{4T_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

式中Δθ为角度的采样间距，单位弧度；计算连续3个时间样点上的最大幅值b_ij，若a_i,j＞0，则b_ij＝max(a_i-1,j,a_i,j,a_i+1,j)，否则b_ij＝min(a_i-1,j,a_i,j,a_i+1,j)；引入两个无量纲参数α和β，一般取α＝0.3，β＝0.5；先计算菲涅耳带的平均幅值

$d=\frac{\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}\left|b_{\mathrm{ik}}\right|}{m_2-m_1+1}$

式中k代表近似菲涅耳带内的角度索引值，b_ik是角度k对应的最大幅值；若d≥αA，则继续计算菲涅耳带内的相关系数

$c=\frac{\left|\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}{b}_{\mathrm{ik}}\right|}{(m_2-m_1+1)d}$

c是无量纲数，若c≥β，则判断深度T_i和角度θ_j对应一个稳相点；对全部角度的角度时间样点集循环，记录得到的稳相点；d)对每个稳相点，重复m₁＝m₁-1和m₂＝m₂+1并计算d和c，当不能满足d≥αA和c≥β时，记录下对应的m₁和m₂；e)在整个时间深度上定义等间距的连续时窗，取50个样点为一个时窗；在每个时窗内，对全部的稳相点，取最小的m₁和其对应的稳相点的角度θ₁和时间深度T₁，取最大的m₂和其对应的稳相点的角度θ₂和时间深度T₂；角度θ₁-m₁·Δθ和T₁，角度θ₂+m₂·Δθ和T₂就分别是不同时间深度上菲涅耳带的下、上界；时窗内没有稳相点就忽略该时窗；f)对得到的各个菲涅耳带的下、上界分别沿时间深度线性插值，即得到该CDP点的初步菲涅耳带的界限曲线；

对倾角拾取CDP位置处的两个倾角道集，分别将拾取的初步菲涅耳带的上、下界曲线投影到倾角道集上，观察曲线是否正确保留了反射同相轴的菲涅耳带，对曲线做必要的修改，即得到该CDP的菲涅耳带；

分别对选定成像线的全部倾角拾取CDP位置的两个菲涅耳带的四组曲线进行空间插值，就可获得拟成像区块上全部成像点的菲涅耳带；将存储上、下界倾角对应的正切值；

偏移成像单元，还用于对叠前地震资料的全部地震道循环；对每一地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为成像点处的偏移速度，即均方根速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处对应的反射界面的两个倾角θ_x和θ_y的正切值为

tanθ_x＝(x₁+g cosφ)/(V_rmsT₀)

tanθ_y＝(y₁+g sinφ)/(V_rmsT₀)

读取成像点的菲涅耳带在两个倾角方向上的上、下界角度的正切值α₊，α_-，β₊和β_-，计算正切值的平均值a₀＝(α₊+α_-)/2和b₀＝(β₊+β_-)/2以及其变化的宽度a₁＝(α₊-α_-)/2和b₁＝(β₊-β_-)/2；计算菲涅耳带切除时的衰减带宽度和其中Δθ为倾角道集的角度采样间距，单位弧度；计算无量纲的衰减系数数组

$d_i=0.3+0.7\sin \left(\frac{(i-1)\mathrm{\π}}{100}\right),i=\mathrm{1,51}$

式中d_i是数组的元素，i是数组的索引数；

计算正切值的差x_A＝|tanθ_x-a₀|，若x_A＞a₁，不需近一步计算该成像点的偏移幅值；否则，进一步计算正切值的差y_A＝|tanθ_y-b₀|，判断是否满足y_A≤b₁；若满足，计算炮点和检波点到成像点的走时和成像权系数，由走时和成像权系数以及该地震道的地震记录确定偏移幅值；计算平行测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为c₁＝a₂/50和c₂＝a₁+0.1；计算正整数M＝2-int(x_A/(a₁-a₂))，进一步计算用于确定衰减系数的正整数

N₁＝int((a₁-x_A)/c_M)+1+50(M-1)

计算垂直测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为d₁＝b₂/50和d₂＝b₁+0.1；计算正整数M＝2-int(y_A/(b₁-b₂))，进一步计算用于确定衰减系数的正整数

N₂＝int((b₁-y_A)/d_M)+1+50(M-1)

依据N₁和N₂从衰减系数数组中提取对应的元素和求得无量纲的衰减系数将衰减系数η乘上前面得到的偏移幅值，将该结果累加到存放该成像点偏移结果的数组中对应的偏移距上，就得到仅对倾角域偏移道集中的菲涅耳带部分叠加成像的高信噪比偏移结果；

数据体形成单元，还用于收集三维稳相叠前时间偏移在各计算节点上的偏移结果，得到三维稳相叠前时间偏移的共反射点道集；对得到的三维稳相叠前时间偏移的全部共反射点道集做剩余动校正；偏移和动校正过程将使得偏移距成像道出现拉伸和畸变，将拉伸超过50％的部分和畸变部分的幅值置为零，将不同偏移距的偏移结果叠加，形成三维成像数据体。

本发明实施例涉及的系统通过生成倾角道集、在倾角道集上拾取菲涅耳带和通过插值获得全部成像点的菲涅耳带等手段，有效压制了偏移噪音和叠前资料中残存的噪音，可生成高信噪比的共反射点偏移道集和偏移叠加剖面；所生成的高信噪比共反射点偏移道集可更好地服务于叠前反演等油气和流体检测技术，高信噪比的偏移叠加剖面可更好地指示地下构造的形态、断裂以及地层沉积样式。

尽管为示例目的，已经公开了本发明的优选实施例，本领域的技术人员将意识到各种改进、增加和取代也是可能的，因此，本发明的范围应当不限于上述实施例。

Claims (9)

1.一种三维倾角域稳相叠前时间偏移方法，其特征在于，包括：

读取叠前地震资料，确定偏移速度场和拉伸切除参数；

依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集；

基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带；

利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果；

收集三维稳相叠前时间偏移的偏移结果，形成三维成像数据体；

将三维成像数据体转换为地下反射构造的剖面图像。

2.如权利要求1所述的方法，其特征在于，读取叠前地震资料，确定偏移速度场和拉伸切除参数的步骤包括：

读取叠前地震资料，对叠前地震资料做线性干扰衰减；对完成了线性干扰衰减的叠前地震资料，选取部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到的动校正速度做横向平均，作为初始偏移速度；

将完成了线性干扰衰减的叠前地震资料按偏移距大小排序，基于给定的共反射点道集中成像道的偏移距间距分组，将不同组叠前地震资料存放到集群计算机的不同计算节点上，利用初始偏移速度，对已存放到集群计算机各个计算节点上的叠前地震资料，应用三维叠前时间偏移方法，进行并行的偏移计算，收集各计算节点的偏移结果，形成常规叠前时间偏移的共反射点道集；

对常规叠前时间偏移的共反射点道集，利用初始偏移速度做常规的反动校，再做常规的动校正得到新的速度，对新速度的倒数做空间平滑处理，空间平滑处理后得到的平滑速度场可作为偏移速度场；根据动校正后共反射点道集中大偏移距地震道的动校拉伸情况，确定拉伸切除参数。

3.如权利要求1所述的方法，其特征在于，依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集的步骤包括：

在拟成像的区块上选定成像线，定义成像线是平行于测线的直线，CDP位置指示成像线上等间距的离散点；在选定的成像线上的每个CDP处，定义两个二维数组存放偏移结果，一个二维数组的一个维对应时间深度，另一个维是反射界面沿平行测线方向的倾角，该数组存储的偏移结果即是该CDP的平行测线方向的倾角道集；第二个二维数组的一个维对应时间深度，另一个维是反射界面沿垂直测线方向的倾角，该数组存储的偏移结果即是该CDP的垂直测线方向的倾角道集；通过偏移计算，在选定的成像线上每一个CDP处得到上述的两个数组，即是生成分别对应平行测线和垂直测线方向倾角的两个倾角道集。

4.如权利要求3所述的方法，其特征在于，通过偏移计算，在选定的成像线上每一个CDP处得到上述的两个数组，即是生成分别对应平行测线和垂直测线方向倾角的两个倾角道集的步骤包括：在偏移计算中，求取反射界面的倾角；对每个地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为从偏移速度场中拾取的成像点处的偏移速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处的偏移结果所对应的反射界面的两个倾角θ_x和θ_y分别满足：

tanθ_x＝(x₁+gcosφ)/(V_rmsT₀)；

tanθ_y＝(y₁+gsinφ)/(V_rmsT₀)；

式中θ_x是反射界面沿平行测线方向的倾角，θ_y是反射界面沿垂直测线方向的倾角；

使用叠前地震资料，进行并行的偏移计算；对每一个选定的成像线，依据事先确定的偏移孔径，确定对该成像线有贡献的叠前地震资料；对这些叠前地震资料中的每一个地震道，根据偏移孔径和拉伸切除参数，确定在不同CDP处的起始成像的时间深度；利用偏移速度场，求得走时，再通过插值计算得到该地震道在成像点(x₀,y₀,T₀)处的偏移结果v；就这一偏移结果，x₀和y₀指示了CDP位置，T₀指示了在时间深度维上的位置，由计算得到的tanθ_x，可进一步确定在平行测线方向倾角维上的位置；这样，可将偏移结果v累加到第一个二维数组，即平行测线方向的倾角道集上；由计算得到的tanθ_y确定在垂直测线方向倾角维上的位置，再将偏移结果v累加到第二个二维数组，即垂直测线方向的倾角道集上；完成对全部有贡献的叠前地震资料中全部地震道的循环，收集各计算节点上的倾角道集并累加，即得到分别对应平行测线和垂直测线方向倾角的两个倾角道集。

5.如权利要求1所述的方法，其特征在于，基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带的步骤包括：对每一成像线的两个倾角道集，分别沿倾角方向进行边缘平滑衰减，在每个CDP处，分别沿两个倾角方向进行叠加，最后将两个倾角道集的叠加结果累加；这样可得到该成像线的偏移叠加剖面；

根据偏移叠加剖面上的断层和同相轴的变化情况，确定一组倾角拾取CDP位置；

对每个倾角拾取CDP位置，在偏移叠加剖面上拾取该CDP的成像道，做增益补偿，然后计算其幅值绝对值的平均值A₀；令a_i-1,a_i,a_i+1是成像道的连续3个时间样点上的幅值，若|a_i|≥0.2A₀，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i；对全部时间样点循环，可得到该CDP的波峰时间样点集；

利用波峰时间样点集，在该CDP的两个倾角道集中搜寻稳相点和菲涅耳带；分别对每个道集进行如下步骤，确定初步菲涅耳带，步骤如下：a)对倾角道集做带通滤波和增益补偿，在时间上取4个时窗，在每个时窗内计算倾角道集的幅值绝对值的平均值；b)对每个倾角的采样点，对波峰时间样点集循环，令i是时间样点的索引值，令a_i-1,a_i,a_i+1是该角度处在i左右连续3个时间样点上的幅值，令A是i所在的时窗内的幅值绝对值的平均值，若|a_i|≥0.2A，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i，得到该角度处的角度时间样点集；若角度时间样点集内样点间的时间采样点数小于2，剔除一个样点，保留幅值绝对值大的样点；c)令i和j为角度时间样点集内样点的时间和角度索引值，a_i,j为道集上该点的幅值；根据叠前地震资料的主频f₀，单位1/秒，用i和j对应的时间深度T_i和角度θ_j，角度θ_j的单位是弧度，可计算出倾角域菲涅耳带的两个方向的角度宽度m₁和m₂，即

$m_1=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j-\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{{4T}_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

$m_2=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j+\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{{4T}_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

式中Δθ为角度的采样间距，单位弧度；计算连续3个时间样点上的最大幅值b_ij，若a_i,j＞0，则否则b_ij＝min(a_i-1,j,a_i,j,a_i+1,j)；引入两个无量纲参数α和β，一般取α＝0.3，β＝0.5；先计算菲涅耳带的平均幅值

$d=\frac{\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}\left|b_{\mathrm{ik}}\right|}{m_2-m_1+1}$

式中k代表菲涅耳带内的角度索引值，b_ik是角度k对应的最大幅值；若d≥αA，则继续计算菲涅耳带内的相关系数

$c=\frac{\left|\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}{b}_{\mathrm{ik}}\right|}{(m_2-m_1+1)d}$

c是无量纲数，若c≥β，则判断深度T_i和角度θ_j对应一个稳相点；对全部角度的角度时间样点集循环，记录得到的稳相点；d)对每个稳相点，重复m₁＝m₁-1和m₂＝m₂+1并计算d和c，当不能满足d≥αA和c≥β时，记录下对应的m₁和m₂；e)在整个时间深度上定义等间距的连续时窗，取50个样点为一个时窗；在每个时窗内，对全部的稳相点，取最小的m₁和其对应的稳相点的角度θ₁和时间深度T₁，取最大的m₂和其对应的稳相点的角度θ₂和时间深度T₂；角度θ₁-m₁·Δθ和T₁，角度θ₂+m₂·Δθ和T₂就分别是不同时间深度上菲涅耳带的下、上界；时窗内没有稳相点就忽略该时窗；f)对得到的各个菲涅耳带的下、上界分别沿时间深度线性插值，即得到该CDP点的初步菲涅耳带的界限曲线；

对倾角拾取CDP位置处的两个倾角道集，分别将拾取的初步菲涅耳带的上、下界曲线投影到倾角道集上，观察曲线是否正确保留了反射同相轴的菲涅耳带，对曲线做必要的修改，即得到该CDP的菲涅耳带；

分别对选定成像线的全部倾角拾取CDP位置的两个菲涅耳带的四组曲线进行空间插值，就可获得拟成像区块上全部成像点的菲涅耳带；将存储上、下界倾角对应的正切值。

6.如权利要求1所述的方法，其特征在于，利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果，包括：对叠前地震资料的全部地震道循环；对每一地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为成像点处的偏移速度，即均方根速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量b＝x₁x₂+y₁y₂+2(V_rmsT₀)²，然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处对应的反射界面的两个倾角θ_x和θ_y的正切值为

tanθ_x＝(x₁+gcosφ)/(V_rmsT₀)

tanθ_y＝(y₁+gsinφ)/(V_rmsT₀)

读取成像点的菲涅耳带在两个倾角方向上的上、下界角度的正切值α₊，α_-，β₊和β_-，计算正切值的平均值a₀＝(α₊+α_-)/2和b₀＝(β₊+β_-)/2以及其变化的宽度a₁＝(α₊-α_-)/2和b₁＝(β₊-β_-)/2；计算菲涅耳带切除时的衰减带宽度和其中Δθ为倾角道集的角度采样间距，单位弧度；计算无量纲的衰减系数数组

$d_i=0.3+0.7\sin \left(\frac{(i-1)\mathrm{\π}}{100}\right),$ i＝1,51

式中d_i是数组的元素，i是数组的索引数；

计算正切值的差x_A＝|tanθ_x-a₀|，若x_A＞a₁，不需近一步计算该成像点的偏移幅值；否则，进一步计算正切值的差y_A＝|tanθ_y-b₀|，判断是否满足y_A≤b₁；若满足，计算炮点和检波点到成像点的走时和成像权系数，由走时和成像权系数以及该地震道的地震记录确定偏移幅值；计算平行测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为c₁＝a₂/50和c₂＝a₁+0.1；计算正整数M＝2-int(x_A/(a₁-a₂))，进一步计算用于确定衰减系数的正整数

N₁＝int((a₁-x_A)/c_M)+1+50(M-1)

计算垂直测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为d₁＝b₂/50和d₂＝b₁+0.1；计算正整数M＝2-int(y_A/(b₁-b₂))，进一步计算用于确定衰减系数的正整数

N₂＝int((b₁-y_A)/d_M)+1+50(M-1)

依据N₁和N₂从衰减系数数组中提取对应的元素和求得无量纲的衰减系数将衰减系数η乘上前面得到的偏移幅值，将该结果累加到存放该成像点偏移结果的数组中对应的偏移距上，就得到仅对倾角域偏移道集中的菲涅耳带部分叠加成像的高信噪比偏移结果。

7.如权利要求1所述的方法，其特征在于，收集三维稳相叠前时间偏移在各计算节点上的偏移结果，得到三维稳相叠前时间偏移的共反射点道集；对得到的三维稳相叠前时间偏移的全部共反射点道集做剩余动校正；偏移和动校正过程将使得偏移距成像道出现拉伸和畸变，将拉伸超过50％的部分和畸变部分的幅值置为零，将不同偏移距的偏移结果叠加，形成三维成像数据体。

8.一种三维倾角域稳相叠前时间偏移系统，其特征在于，包括：

参数准备单元，用于读取叠前地震资料，确定偏移速度场和拉伸切除参数；

倾角道集生成单元，用于依据偏移速度场和拉伸切除参数，对选定的成像线，生成分别对应平行测线和垂直测线方向倾角的两个倾角道集；

菲涅耳带生成单元，用于基于选定的成像线处的两个倾角道集，拾取两个倾角方向的菲涅耳带，通过插值获得拟成像区块上全部成像点的菲涅耳带；

偏移成像单元，用于利用获得的全部成像点的菲涅耳带，进行三维稳相叠前时间偏移，得到对应不同偏移距的偏移结果；

数据体形成单元，用于收集三维稳相叠前时间偏移的偏移结果，形成三维成像数据体；

数据体转换单元，用于将三维成像数据体转换为地下反射构造的剖面图像。

9.如权利要求8所述的系统，其特征在于，

参数准备单元，具体用于读取叠前地震资料，对叠前地震资料做线性干扰衰减；对完成了线性干扰衰减的叠前地震资料，选取部分共中心点，抽取共中心点道集，对抽取的道集作常规的动校正速度拾取，对所得到的动校正速度做横向平均，作为初始偏移速度；

将完成了线性干扰衰减的叠前地震资料按偏移距大小排序，基于给定的共反射点道集中成像道的偏移距间距分组，将不同组叠前地震资料存放到集群计算机的不同计算节点上，利用初始偏移速度，对已存放到集群计算机各个计算节点上的叠前地震资料，应用三维叠前时间偏移方法，进行并行的偏移计算，收集各计算节点的偏移结果，形成常规叠前时间偏移的共反射点道集；

对常规叠前时间偏移的共反射点道集，利用初始偏移速度做常规的反动校，再做常规的动校正得到新的速度，对新速度的倒数做空间平滑处理，空间平滑处理后得到的平滑速度场可作为偏移速度场；根据动校正后共反射点道集中大偏移距地震道的动校拉伸情况，确定拉伸切除参数；

倾角道集生成单元，还用于在拟成像的区块上选定成像线，定义成像线是平行于测线的直线，CDP位置指示成像线上等间距的离散点；在选定的成像线上的每个CDP处，定义两个二维数组存放偏移结果，一个二维数组的一个维对应时间深度，另一个维是反射界面沿平行测线方向的倾角，该数组存储的偏移结果即是该CDP的平行测线方向的倾角道集；第二个二维数组的一个维对应时间深度，另一个维是反射界面沿垂直测线方向的倾角，该数组存储的偏移结果即是该CDP的垂直测线方向的倾角道集；通过偏移计算，在选定的成像线上每一个CDP处得到上述的两个数组；

在偏移计算中，求取反射界面的倾角；对每个地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为从偏移速度场中拾取的成像点处的偏移速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处的偏移结果所对应的反射界面的两个倾角θ_x和θ_y分别满足：

tanθ_x＝(x₁+gcosφ)/(V_rmsT₀)；

tanθ_y＝(y₁+gsinφ)/(V_rmsT₀)；

式中θ_x是反射界面沿平行测线方向的倾角，θ_y是反射界面沿垂直测线方向的倾角；

使用叠前地震资料，进行并行的偏移计算；对每一个选定的成像线，依据事先确定的偏移孔径，确定对该成像线有贡献的叠前地震资料；对这些叠前地震资料中的每一个地震道，根据偏移孔径和拉伸切除参数，确定在不同CDP处的起始成像的时间深度；利用偏移速度场，求得走时，再通过插值计算得到该地震道在成像点(x₀,y₀,T₀)处的偏移结果v；就这一偏移结果，x₀和y₀指示了CDP位置，T₀指示了在时间深度维上的位置，由计算得到的tanθ_x，进一步确定在平行测线方向倾角维上的位置；这样，可将偏移结果v累加第一个二维数组，即到平行测线方向的倾角道集上；由计算得到的tanθ_y确定在垂直测线方向倾角维上的位置，再将偏移结果v累加到第二个二维数组，即垂直测线方向的倾角道集上；完成对全部有贡献的叠前地震资料中全部地震道的循环，收集各计算节点上的倾角道集并累加，即得到分别对应平行测线和垂直测线方向倾角的两个倾角道集；

菲涅耳带生成单元，还用于对每一成像线的两个倾角道集，分别沿倾角方向进行边缘平滑衰减，在每个CDP处，分别沿两个倾角方向进行叠加，最后将两个倾角道集的叠加结果累加；这样可得到该成像线的偏移叠加剖面；

根据偏移叠加剖面上的断层和同相轴的变化情况，确定一组倾角拾取CDP位置；

对每个倾角拾取CDP位置，在偏移叠加剖面上拾取该CDP的成像道，做增益补偿，然后计算其幅值绝对值的平均值A₀；令a_i-1,a_i,a_i+1是成像道的连续3个时间样点上的幅值，若|a_i|≥0.2A₀，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i；对全部时间样点循环，可得到该CDP的波峰时间样点集；

利用波峰时间样点集，在该CDP的两个倾角道集中搜寻稳相点和菲涅耳带；分别对每个道集进行如下步骤，确定初步菲涅耳带，步骤如下：a)对倾角道集做带通滤波和增益补偿，在时间上取4个时窗，在每个时窗内计算倾角道集的幅值绝对值的平均值；b)对每个倾角的采样点，对波峰时间样点集循环，令i是时间样点的索引值，令a_i-1,a_i,a_i+1是该角度处在i左右连续3个时间样点上的幅值，令A是i所在的时窗内的幅值绝对值的平均值，若|a_i|≥0.2A，则再判断是否满足|a_i+1-a_i-1|＜0.02|a_i|，若满足，进一步判断是否满足(a_i+1-a_i)(a_i-a_i-1)＜0，若还满足，则记录时间样点i，得到该角度处的角度时间样点集；若角度时间样点集内样点间的时间采样点数小于2，剔除一个样点，保留幅值绝对值大的样点；c)令i和j为角度时间样点集内样点的时间和角度索引值，a_i,j为道集上该点的幅值；根据叠前地震资料的主频f₀，单位1/秒，用i和j对应的时间深度T_i和角度θ_j，角度θ_j的单位是弧度，计算出倾角域菲涅耳带的两个方向的角度宽度m₁和m₂，即

$m_1=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j-\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{{4T}_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

$m_2=\frac{1}{\mathrm{\Δ\θ}}\{{\mathrm{\θ}}_j-\arctan [\tan {\mathrm{\θ}}_j+\sqrt{\frac{1}{T_i{f}_0}(1-\frac{1}{{4T}_i{f}_0})(1+{\tan }^2{\mathrm{\θ}}_j)}\left]\right\}$

式中Δθ为角度的采样间距，单位弧度；计算连续3个时间样点上的最大幅值b_ij，若a_i,j＞0，则否则b_ij＝min(a_i-1,j,a_i,j,a_i+1,j)；引入两个无量纲参数α和β，一般取α＝0.3，β＝0.5；先计算菲涅耳带的平均幅值

$d=\frac{\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}\left|b_{\mathrm{ik}}\right|}{m_2-m_1+1}$

式中k代表菲涅耳带内的角度索引值，b_ik是角度k对应的最大幅值；若d≥αA，则继续计算菲涅耳带内的相关系数

$c=\frac{\left|\underset{k=j-m_1}{\overset{j+m_2}{\mathrm{\Σ}}}{b}_{\mathrm{ik}}\right|}{(m_2-m_1+1)d}$

c是无量纲数，若c≥β，则判断深度T_i和角度θ_j对应一个稳相点；对全部角度的角度时间样点集循环，记录得到的稳相点；d)对每个稳相点，重复m₁＝m₁-1和m₂＝m₂+1并计算d和c，当不能满足d≥αA和c≥β时，记录下对应的m₁和m₂；e)在整个时间深度上定义等间距的连续时窗，取50个样点为一个时窗；在每个时窗内，对全部的稳相点，取最小的m1和其对应的稳相点的角度θ₁和时间深度T₁，取最大的m₂和其对应的稳相点的角度θ₂和时间深度T₂；角度θ₁-m₁·Δθ和T₁，角度θ₂+m₂·Δθ和T₂就分别是不同时间深度上菲涅耳带的下、上界；时窗内没有稳相点就忽略该时窗；f)对得到的各个菲涅耳带的下、上界分别沿时间深度线性插值，即得到该CDP点的初步菲涅耳带的界限曲线；

对倾角拾取CDP位置处的两个倾角道集，分别将拾取的初步菲涅耳带的上、下界曲线投影到倾角道集上，观察曲线是否正确保留了反射同相轴的菲涅耳带，对曲线做必要的修改，即得到该CDP的菲涅耳带；

分别对选定成像线的全部倾角拾取CDP位置的两个菲涅耳带的四组曲线进行空间插值，就可获得拟成像区块上全部成像点的菲涅耳带；将存储上、下界倾角对应的正切值；

偏移成像单元，还用于对叠前地震资料的全部地震道循环；对每一地震道，令x_s和y_s是炮点的水平坐标，x_r和y_r是检波点的水平坐标，φ是方位角；令V_rms为成像点处的偏移速度，即均方根速度，x₀、y₀和T₀分别是成像点的水平坐标和时间深度；定义如下距离变量：x₁＝x_s-x₀，y₁＝y_s-y₀，x₂＝x_s+x_r-2x₀，y₂＝y_s+y_r-2y₀；用这些已知量计算距离变量a＝x₁cosφ+y₁sinφ，c＝x₂cosφ+y₂sinφ， $d=\sqrt{x_1^2+y_1^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}+\sqrt{{(x_r-x_0)}^2+{(y_r-y_0)}^2+{\left(V_{\mathrm{rms}}{T}_0\right)}^2}$ 和距离的平方量b＝x₁x₂+y₁y₂+2(V_rmsT₀)²，然后再计算距离g＝(a·d²-b·c)/(c²-d²)；则该地震道在成像点处对应的反射界面的两个倾角θ_x和θ_y的正切值为

tanθ_x＝(x₁+gcosφ)/(V_rmsT₀)

tanθ_y＝(y₁+gsinφ)/(V_rmsT₀)

读取成像点的菲涅耳带在两个倾角方向上的上、下界角度的正切值α₊，α_-，β₊和β_-，计算正切值的平均值a₀＝(α₊+α_-)/2和b₀＝(β₊+β_-)/2以及其变化的宽度a₁＝(α₊-α_-)/2和b₁＝(β₊-β_-)/2；计算菲涅耳带切除时的衰减带宽度和其中Δθ为倾角道集的角度采样间距，单位弧度；计算无量纲的衰减系数数组

$d_i=0.3+0.7\sin \left(\frac{(i-1)\mathrm{\π}}{100}\right),$ i＝1,51

式中d_i是数组的元素，i是数组的索引数；

计算正切值的差x_A＝|tanθ_x-a₀|，若x_A＞a₁，不需近一步计算该成像点的偏移幅值；否则，进一步计算正切值的差y_A＝|tanθ_y-b₀|，判断是否满足y_A≤b₁；若满足，计算炮点和检波点到成像点的走时和成像权系数，由走时和成像权系数以及该地震道的地震记录确定偏移幅值；计算平行测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为c₁＝a₂/50和c₂＝a₁+0.1；计算正整数M＝2-int(x_A/(a₁-a₂))，进一步计算用于确定衰减系数的正整数

N₁＝int((a₁-x_A)/c_M)+1+50(M-1)

计算垂直测线方向倾角有关的幅值衰减系数，为此，定义一个有两个元素的数组，其两个元素为d₁＝b₂/50和d₂＝b₁+0.1；计算正整数M＝2-int(y_A/(b₁-b₂))，进一步计算用于确定衰减系数的正整数

N₂＝int((b₁-y_A)/d_M)+1+50(M-1)

依据N₁和N₂从衰减系数数组中提取对应的元素和求得无量纲的衰减系数将衰减系数η乘上前面得到的偏移幅值，将该结果累加到存放该成像点偏移结果的数组中对应的偏移距上，就得到仅对倾角域偏移道集中的菲涅耳带部分叠加成像的高信噪比偏移结果；

数据体形成单元，还用于收集三维稳相叠前时间偏移在各计算节点上的偏移结果，得到三维稳相叠前时间偏移的共反射点道集；对得到的三维稳相叠前时间偏移的全部共反射点道集做剩余动校正；偏移和动校正过程将使得偏移距成像道出现拉伸和畸变，将拉伸超过50％的部分和畸变部分的幅值置为零，将不同偏移距的偏移结果叠加，形成三维成像数据体。