专利名称：叠加前地震数据的向下外推方法
技术领域：
一般地说，本发明涉及地球物理勘探领域。特别是，本发明涉及地震数据处理领域。具体地说，本发明是用于地震偏移的叠加前地震数据的向下外推的方法。
背景技术：
地震偏移是由记录的地震数据构造出能定义感兴趣的地下地层反射体表面的过程。这样，偏移过程提供地球在深度或时间上的图像。这是要解释与地震能量从地震源到地震接收器过程中的透射和反射相关联的位置(运动学的)和振幅(动力学的)两种效应。尽管速度的垂直变化是最常见的，但也会遇到横向变化。这些速度变化源于若干原因，包括地球各层的差异性挤压，盐底辟侧翼上的隆起、以及近端(页岩的)和远端(砂质到碳质的)大陆架位置之间的沉积动力学变化。
多年来在地球物理处理业界已经认识到地震偏移应在叠加前进行并在深度域进行，而不是在叠加后进行，以便在构造复杂的区域得到最佳叠加图像。此处，叠加前深度偏移的一个优点是为其后的AVO(振幅对炮检距)分析提供最佳准备的数据。传统上，叠加前深度偏移是通过应用基尔霍夫方法进行的。然而，因为近来计算硬件的进步和有效的外推器设计上的改进，基于单向波动方程解的方法在变得越来越普及。
文献中已清楚地确认，基于波动方程的方法提供暗含多重路径的运动学优点，从而能更准确地勾画复杂盖层下的结构。然而，过去对于通过基于波动方程的方法可能(或应该)实现的动力学优点却讨论得相当少。这并不令人惊奇，因为基尔霍夫叠加前深度偏移已经经过了比其基于波动方程的同类方法长久得多的演化。这种演化的一部分便是各种因素或策略的研发以考虑在数据采集过程中的几何发散和空间不规则性。
一种叠加前深度偏移方法单个地处理共炮点道集(common shotgather)或炮点记录。因此这一方法被称作“炮点记录”或“炮点剖面”偏移。炮点记录偏移能提供极好的成像准确性和极好的振幅保持能力。这些优点已使炮点记录偏移成为更普遍采用的波动理论偏移方法之一。
炮点记录偏移的一种早期实现是在Moshe Reshef和DanKosloff的论文“共炮点道集的偏移(Migration of common shot gathers)”(Geophysics，第51卷，第2期(1986年2月)，第324-331页)中描述的。
地震偏移包含两个步骤(1)波外推和(2)成像。本申请涉及第一步骤，即在地震偏移中的向下外推。在地震偏移中相应的成像方法在这同一发明者的共同待决美国专利申请“叠加前地震数据成像方法”中被描述，该专利被转让给本专利申请的同一受让人。
向下波外推造成一个波场，该波场是如果震源和接收器都位于深度z时将会记录到的波场的一个近似。各种地震偏移方法的差别之一在于它们实现波外推步骤的方法，这是由于在解波动方程方面的差别。
在Jen Gazdag的论文“利用相移方法的波动方程偏移(Waveequation migration with the phase-shift method)”(Geophysics，第43卷第7期(1978年12月)，第1342-1351页)中描述了Gazdag外推。Gazdag(1978)公开了一种基于使用付立叶变换数值求解波动方程的零点偏移(zero-offset)地震数据的偏移方法。在Gazdag外推方法中，利用相移，即在频率域中地震数据付立叶系数的相角的转动，使源和接收器的位置向下外推。这些计算等效于乘以一个单位模的复数，保证无条件的稳定性。此外，该方法处置大传播(倾)角的情况。然而，Gazdag的相移法只对横向均匀速度有效。
在Jen Gazdag和Piero Sguazzero的论文“利用相移加内插的地震数据偏移(Migration of seismic data by phase shift plusinterpolarion)”(Geophysics，第49卷第2期(1984年2月)，第124-131页)中，Gazdag相移法被扩展到有横向速度变化的介质。Gazdag和Sguazzero(1984)公开说明了一种两步外推法，通常被称作相移加内插。在第一步，由Gazdag相移法使用多个横向均匀速度场对波场进行外推，产生若干参考波场。在第二步，由这些参考波场进行内插，从而计算出实际波场。如Gazdag的相移法一样，Gazdag的相移加内插法是无条件稳定的，且能准确处置很陡的倾斜。随着所用参考速度个数增加，偏移的总体准确度提高，但计算成本相应提高。有很大横向速度变化的地区需要多达5个参考波场，每个需求一对空间付立叶变换。这样，Gazdag相移加内插法会是很费钱的。
分步(split-step)付立叶外推(也称作相位筛选(phase-screen)法)是Gazdag相移法的修改，以在每个偏移间隔接受横向速度变化。在Stoffa、P.L.Fokkema、J.T.、de Luna Freire、R.M.以及Kessinger、W.P.的论文“分步付立叶偏移(Split-stepFourier migration)”(Geophysics，第55卷第4期(1990年4月)，第410-421页)中描述了分步付立叶外推。Stoffa等(1990)公开说明了一种两步相移法以处置横向速度变化首先定义一个参考慢度(速度的倒数)作为该偏移间隔中的平均慢度，然后使用空间变化扰动或校正项来确定源贡献。第一相移与常速度Gazdag相移法完全相同。平均慢度确定一个参考垂直波数，它用于在频率-波数阈内向下续该数据。第二相移是一个校正项，它根据在每个空间位置的实际慢度和参考慢度之差提供一个时移。对于每个偏移间隔重复这两个相移。如果相对速度变化小的话，分步付立叶外推对于一个宽的入射角范围都能工作得很好。它还能接受大的速度变化，但只对于一个有限的角范围有效。
在Lian-Jie Huang和Michael Fehler的论文“分步付立叶传播器的准确度分析对地震建模和偏移的意义(Accuracy analysisof the split-step Fourier propagatorImplications for Seismicmodeling and Migration)”(Bulliten of the SeismoLogicalSociety of America，第88卷第1期(1998年2月)，第18-19页)中描述了对分步付立叶方法的误差分析。分步付立叶法在单向标量波动方程中展开平方根算子和分裂指数算子。Huang和Fehler(1998)表明对于小于约45°传播角和小于约10％的相对速度变化，这些近似给出的误差低于5％。对于小的相对速度变化，分步付立叶方法对于传播角高达约60°仍是准确的。
在横向均匀介质中，通过应用单个相移能实现准确的波外推，这是在频率-波数阈中单向波动方程的解。在频率-波数阈的波动方程有时被称作“平方根方程”。然而，在横向非均匀介质中，是通过利用平方根方程的有理近似实现波外推。这种近似通常是通过有限差分方法实现的。使平方根方程有理化的两种方法是泰勒级数方法和连分式方法。泰勒级数法被称作显式算法，而连分式法被称作隐式算法。隐式有限差分外推法在本技术领域特别有名。
例如，在下面两篇论文中描述了隐式有限差分外推法Myumg W.Lee和Sang Y.Suh的“单向波动方程的最优化(Optimization ofone-way wave equations)”(Geophysics，第50卷，第10期(1985年10月)，第1634-1637页)和Zhiming Li的“补偿3D偏移和建模中的有限差分误差(Compensating finite-difference errorsin 3-D migration and modeling)”(Geophysics，第56卷，第10期(1991年10月)，第1650-1660页)。Lee和Suh(1985)用最小二乘法优化对平方根方程的二阶近似的系数。Li(1991)利用应用于少数几步有限差分外推的相移补偿滤波器来补偿3D偏移中的横向(crossline)和纵向(inline)分裂误差。相移补偿滤波器是基于Gazdag相移或相移加内插方法。
在Dietrich Ristow和Thomas Ruhl的论文“付立叶有限差分偏移(Fourier finite-difference migration)”(Geophysics，第59卷，第12期(1994年12月)，第1882-1893页)中描述了一种混合偏移法。Ristow和Ruhl(1994)将频率-波数域中的Gazdag相移外推与其后在频率-空间域中的有限差分外推组合在一起。这一方法将分步付立叶方法与有限差分方法组合在一起。有限差分步骤校正大横向速度扰动对绕射项的影响。在这方面，有限差分补偿空间域中对薄透镜传播的初始扰动。
如果横向速度变化足够小，则初始的常速度相移造成绕射项中的主要时差(moveout)。这样，有限差分计算可能不需要任何相位校正。在这些环境条件下，对于每个波场，在(ω，κ)域中只有一个Gazdag相移，跟随有在(ω，x)域中的一个薄透镜相移和一个有限差分计算。然而，如果横向速度变化足够大，则在有限差分之后可能需要一个或多个相位校正，其中每一个需要前向和反向空间付立叶变换和一个相移。这些相移的成本增加到初始Gazdag相移的成本中。
在下文中讨论了这种混合偏移方法的新近的一个实现JamesSun、Carl Notfors、Sam Gray和Yu Zhang的“3D叠加前共炮点深度偏移一种适应于结构的实现(3D Pre-stack common shot depthmigrationA structural adaptive implemenation)”(SEG Int′lExp.&Ann.Mtg.，San Antonio，Taxas，2001年9月9-14日，详细摘要，第1017-1020页)。Sun等人(2001)描述了为每一偏移间隔在Gazdag和有限差分方法之间做出决定的过程。他们还研发出他们自己的有限差分系数组。依然，有限差分系数的计算方法看来有某种任意性和困难。此外，尚没有作者提供一种算法能定量识别需要尾随有限差分计算的条件。
这样，存在着适用于以横向速度变化程度表明的均匀、弱不均匀和强不均匀介质的若干外推算法。毫不奇怪，更准确的算法也是更费钱的算法。因此，希望能为特定外推步骤识别出性能价格比最好的方法。这样，需要一种混合偏移方法，它能利用基于地下速度模型复杂性的自动误差分析为每个深度步长选择最有效的外推器。

发明内容
本发明是用于地震偏移的叠加前地震数据的向下外推方法。通过首先确定地震数据集中的偏移间隔，一个叠加前地震数据集被向下外推。一个混淆条件(aliasing condifion)被用于确定偏移间隔的大小。然后计算出相位最大相对误差，作为频率、传播角以及该偏移间隔中速度相对变化的函数。将相位最大相对误差与最大误差判据进行比较。根据相位最大相对误差与最大误差判据的比较结果，确定用于该偏移间隔的外推类型。外推类型是从一组类型中选出的，其中包含Gazdag相移外推、分步付立叶外推和隐式有限差分外推。


通过参考下文的详细描述和附图，可更容易理解本发明及其优点，其中图1是流程图，说明本发明的叠加前地震数据向下外推方法的一个实施例的处理步骤；图2是流程图，说明本发明的叠加前地震数据误差分析方法的一个实施例的处理步骤；图3是流程图，说明本发明的用于进行Gazdag相移外推的方法的一个实施例的处理步骤；以及图4是流程图，说明本发明的用于进行分步付立叶外推的方法的一个实施例的处理步骤。
尽管将结合其优选实施例描述本发明，但应该理解本发明不限于这些实施例。相反，本发明是要覆盖可能被包括在由所附权利要求定义的本发明范围内的所有变体、修改和等效物。
具体实施例方式
本发明是对由炮点记录方法得到的叠加前地震数据进行偏移的一种方法。该偏移包含叠加前地震数据的向下外推和成像。取决于地下速度模型的复杂性，使用三种不同外推器之一，对每个外推步骤模拟(建模)每个点源炮点的波场。对外推器的选择是在Gazdag相移、分步付立叶以及有限差分外推三种算法之间进行。类似地，使用这同一算法将每个炮点的接收器记录集反外推到底面。对这两个波场使用哪个外推器所作的选择由自动误差分析来指导，并保证对每个深度步长应用最有效的外推器。
这样，本发明包含应用一种适应于结构的混合炮点记录偏移法，它在每个偏移步骤，根据误差分析在三种不同类型外推器之间进行选择。这三种不同类型的单向外推器是(1)用于横向均匀介质的Gazdag相移；(2)用于弱横向不均匀介质的分步付立叶方法；以及(3)用于具有中度到强烈横向速度变化的介质的带相位校正的隐式有限差分方法。
图1显示一个流程图，说明本发明的叠加前地震数据向下外推方法的一个实施例的处理步骤。
首先，在步骤101，选择一个叠加前地震数据集。优选地，所选择的地震数据集被配置成要由炮点记录偏移来处理。
在步骤102，为步骤101中选择的叠加前地震数据集选择一个速度模型。所选择的速度模型包含局部速度c(x，y，z)。该速度模型可由本技术领域已知的任何手段确定。
在步骤103，使用步骤102中确定的速度模型进行两种误差分析。误差分析的目的是确定下一个偏移间隔的大小和在那个偏移间隔要采用的外推方法类型。优选地，进行混淆误差分析以确定下一个偏移间隔的大小。然后，优选地，进行相位误差分析以确定在该偏移间隔采用哪类外推。
在步骤104，使用步骤103中进行的混淆误差分析确定下一个偏移间隔。优选地，混淆误差分析是由本专利申请的同一发明者在美国专利申请09/893,793号“地震数据偏移的方法和系统(Method andSystem for Migration of Seismic Data)”中描述的误差分析，该申请于2001年6月28日提交，并转让给本专利申请的相同受让人。
在步骤105，使用步骤103中进行的相位误差分析确定在步骤104中确定的偏移间隔中要采用的外推类型。外推类型从一组外推方法中选择，该组包含Gazdag相移外推、分步付立叶外推以及隐式有限差分外推。优选地，误差分析是下文中参考图2的流程图描述的本发明实施例。在图2所示优选实施例中，误差分析估计相位最大相对误差并将它与相位的最大可允许相对误差进行比较。
如果误差分析确定在下一个偏移间隔要采用的外推类型是Gazdag相移外推，则过程进入步骤106。如果误差分析确定在下一个偏移间隔要采用的外推类型是分步付立叶外推，则过程进入步骤107。如果误差分析确定在下一个偏移间隔要采用的外推类型是隐式有限差分外推，则过程进入步骤108。
在步骤106，Gazdag相移外推被应用于在步骤103中选择的偏移间隔。优选地，Gazdag相移外推由下文中参考图3的流程图描述的实施例给出。然后过程进入步骤109。
在步骤107，分步付立叶外推被应用于在步骤103中选择的偏移间隔。优选地，分步付立叶外推由下文中参考图4的流程图描述的实施例给出。然后过程进入步骤109。
在步骤108，隐式有限差分外推被应用于在步骤103中选择的偏移间隔。优选地，隐式有限差分外推由本技术领域公知的任何适当方法完成。
在步骤109，确定在步骤101中选择的地震数据集中是否仍存在其他感兴趣的偏移间隔。如果回答为“是”，表明仍存在感兴趣的偏移间隔，于是过程返回步骤103，去进行误差分析以选择下一个偏移间隔和要采用的外推技术类型。如果回答为“否”，表明再没有剩下感兴趣的偏移间隔，于是过程继续到步骤110以结束该过程。
在步骤110，过程结束。该叠加前地震数据集已由混合炮点记录方法进行了偏移。
炮点记录偏移可如下文所述。设u代表在均匀、弹性介质中在空间位置(x1，x2，x3)由地震检波器记录的依赖于时间的位移，它是由脉冲式(如爆炸样的)点源引起的。将运动方程与线性弹性体的本构关系组合起来，位移的第i个矢量分量的位移由下式给出ρ∂2ui∂t2=fi+∂∂xi{[λδijδkl+μ(δikδjl+δilδjk)]∂uk∂xl}...(1)]]>其中λ(x)和μ(x)是随空间变化的拉梅(Lame)参数，ρ(x)是随空间变化的密度，fi是与源关联的体力的第i个矢量分量。
方程(1)包括关于位移和弹性性质二者的导数。如果材料性质只相对于深度改变，则建模和偏移二者都可在均匀层叠加的假定下进行。于是只需要保证穿过层边界的应力和位移的连续性。对于均匀区域内的传播，该方程有简化的矢量形式ρ∂2u∂t2=f+(λ+2μ)▿(▿·u)-μ▿×(▿×u)...(2)]]>方程(2)中的梯度项对应于压缩波运动，而最后一项对应于剪切波运动。在地球的声学模型中，我们只考虑压缩运动。在此假定下，位移u和压强P由下式关联p＝-ρc2·u(3)其中，我们设c2=λ+2μρ...(4)]]>取方程(2)的散度并用方程(3)和(4)替换(λ+2μ)和·u，则可导出波动方程。设x≡(x，y.z)，并假定源处在(xs，ys，zs)并有时间依赖关系S(t)。如果该介质没有密度变化，并忽略关于速度的横向导数，则压强记录遵从“双向”标量波动方程
▿2P(x,y,z,t)-1c2∂2∂t2P(x,y,z,t)=-4πc2δ(x-xs)δ(y-ys)δ(z-zs)S(t)...(5)]]>方程(5)是“非齐次”微分方程，因为它对在源位置和离开源位置处都有效。第二个“齐次”方程，它只对离开源区域处有效，由下式给出▿2P(x,y,z,t)-1c2∂2∂t2P(x,y,z,t)=0...(6)]]>如果介质还包括密度的变化ρ(x，y，z)，则方程(5)和(6)必须分别修改为ρ▿·[1ρ▿P(x,y,z,t)]-1c2∂2∂t2P(x,y,z,t)=-4πc2δ(x-xs)δ(y-ys)δ(z-zs)S(t)...(7)]]>和ρ▿·[1ρ▿ρ(x,y,z,t)]-1c2∂2∂t2P(x,y,z,t)=0...(8)]]>其中“·”代表矢量ρ和(1/ρ)之间的点积。由于这些方程被限定于P波运动，如由压强表示的那样，它们往往被称作声波波动方程。
与上文讨论的Ristow和Ruhl(1994)的付立叶有限差分方法不同，本发明有一种自动选择方法，它只进行为达到某一准确度水平所必需的那些计算。对于分步付立叶外推的相位误差，可导出一个解析表达式，它是频率、传播方向以及速度相对变化的函数。这一表达式与一个最大误差判据一起被用于在外推之前确定分步付立叶方法能否提供可接受的结果。设θ代表相对于垂直方向的传播角，并设c(x，y)/c0表示相对于参考值c0的速度横向变化。在一个外推间隔上的相位相对误差由下式给出
ΔΦΦ(c(x,y))=(ωΔzc0)·[(c(x,y)c0)2-sin2(θ)+(1-cos(θ)-c(x,y)c0)]...(9)]]>这是相位相对误差的一个确切方程式。
当相对速度趋于1或传播方向趋于垂直时，相位误差趋于零。最大误差对应于θ和c/c0的极值。为了识别θmax，可在(ω，κ)域中扫描源和接收器波场。值κmax会被识别出来，它包括半径κ＝κmax的圆内例如95％的功率。为了把这κmax转换成传播角，用于下一个外推步骤的速度切片会被扫描以识别一个最小值Cmin。于是，可以(保守地)设θmax=sin-1(kmax·cminω)]]>计算Δφmax并将它与最大可允许的相对相位误差进行比较。后者必须是一个输入参数。如果误差容限被超过，则下一次外推将要由有限差分方法完成，而不是由分步付立叶方法完成。类似类型的分析能用于确定Gazdag外推的最大相位误差。
上述途径不需任何附加的FFT(快速付立叶变换)。如果刚刚完成了Gzadag或分步付立叶外推，则在逆FFT之前会找到κmax。如果刚刚完成了FD外推及相关联的相位校正，则波场在(ω，κ)(相位校正)域，于是会再次在逆FFT之前找到κmax。最后，如果刚刚完成有限差分外推，但尚不需要后续的相位校正，则下一次外推必须默认为另一次有限差分。
不论在给定深度间隔上使用的外推方法如何，确定κmax需要整个(完全被校正的)波场。这意味着薄透镜校正必须先于每次外推或绕射项的相位校正。本发明在外推和绕射项相位校正后进行所有薄透镜校正。这样，相位校正将不得不重新建立公式，以与总波场一致。
最后，因为炮点波场包括较宽的传播角范围，(实际上是靠近表面的所有角度)，上文讨论的自动选择策略只是最适合于接收器波场。考虑到这些限制，简单地让使用者指定一个最大传播角会更有效。于是我们能绕过对κmax的确定，并且方程(9)能用于炮点波场和接收器波场二者。此外，能为相位校正和薄透镜位置校正二者保持当前的编程。
在实现时，必需选择c0使方程(9)的平均根中的量保持非负。对于一个输入θmax，这是通过扫描速度切片以寻找cmin，然后按下式计算c0实现的c0=cminsin(θmax)...(10)]]>于是，使用方程(9)能对最小速度和最大速度二者进行误差估计。由方程(10)计算出的速度可能大于该外推间隔内的最大速度，从而造成对速度相对变化和相位相对误差二者的保守估计。
在下面的流程图中描述图1的步骤103中提到的误差分析。图2显示的流程图说明本发明的叠加前地震数据误差分析方法的一个实施例的处理步骤。
首先，在步骤201，为在图1的步骤101中选择的地震数据集选择最大可允许的相对相位误差。
在步骤202，选择偏移间隔。偏移间隔的大小强烈地影响外推方法的准确度。准确度必须根据下列因素来衡量(1)对于每个深度样本，必须保证内插过程比精确计算快得多，因为对每个外推间隔它只需要确定一个复指数。(2)在同一准确度水平，与薄透镜项相比，能在大得多的外推间隔上准确地内插绕射项。(3)因为格点之间的横向耦合，绕射项的外推比薄透镜项的外推昂贵得多。
因素(2)和(3)表明，外推间隔的大小应由对绕射项的准确度要求来指导。一旦确立了这一间隔的大小，则能对卓越频率和更高频率精确计算出薄透镜项，而对较低频率能调用内插。一个更简单的但只是稍微更昂贵一些的策略需要对所有频率进行对薄透镜项的精确计算。在本发明中，优选地对所有薄透镜计算都遵循这一策略，不管用于求解绕射项所使用的方法如何。为了减小对存储器的需求，薄透镜计算是在应用该成像条件的同一例行程序中进行的。
绕射项能在比薄透镜项大得多的深度间隔上被内插而且保持同样的准确度，特别是对于小入射角的情况。对于大入射角，我们不能期望很准确地恢复振幅，除非外推步长的大小相当小，即使在卓越频率和对于中等大的速度也是如此。在本发明中，下列三个可选项可用于确定外推间隔的大小。
(1)对于每次外推的深度采样间隔(δz)的大小和格点数被指定。在这一选项下，对所有深度，外推间隔的大小保持相同。
(2)深度采样间隔(δz)的大小被指定。通过对该间隔的最大频率和最小速度应用混淆条件来选择外推间隔的大小。然后该同一间隔用于所有频率。这一途径产生大的间隔，但在卓越频率和较低频率从混淆(相位)观点看它们又是保守的。由于非弹性衰减通常减小较深反射的较高频成份，最好是输入一个fmax对深度的剖面。这样，在较大深度，计算出的外推间隔将比较大，这首先是因为最大频率降低，其次是因为速度通常随深度增大而增大。
(3)深度采样间隔(δz)的大小被指定。然后使用外推间隔的最小速度对每个频率计算单个外推间隔。这一构想能以高达二倍因子减小运行时间，尽管由于内插误差导致较大噪声。
在步骤203，在由步骤202选择的偏移间隔中选择穿过波场的频率切片。
在步骤204，识别出最大波数κmax。优选地，识别出κmax以包括在以κ＝κmax为半径的圆圈中功率的95％。
在步骤205，识别出最小速度cmin。优选地，cmin的识别是通过对下一个外推步骤扫描述速度切片以识别最小速度值。
在步骤206，识别出最大传播角θmax。优选地，θmax是通过设θmax=sin-1(kmax·cminω)]]>得到的，其中使用在步骤203中选择的频率ω，在步骤204中识别出的最大波数κmax，以及在步骤205中识别出的最小速度cmin。
在步骤207，选择一个参考速度c0。优选地，通过应用上述方程(10)选择c0。
在步骤208，对该偏移间隔计算出相位最大相对误差Δφmax。优选地，通过使用在上述方程(9)中对θ代入在步骤206中确定的最大传播角θmax以及在步骤207中选择的参考速度c0，从而计算出Δφmax。
在步骤209，在步骤208中计算出的相位最大相对误差Δφmax与在步骤201中选择的最大可允许的相位相对误差进行比较。
在步骤210，根据在步骤209中最大的和最大可允许的相对误差的比较结果，确定在下一个偏移间隔中要采用的外推方法类型。如果相位的最大相对误差大于最大可允许的相位相对误差，则在下一个偏移间隔期间使用下一个(准确度)较高的外推方法。该外推类型是从一组外推类型中选择的，其中包含Gazdag相移外推、分步付立叶外推和隐式有限差分外推，按考虑的递增顺序排列。如果当前使用的是Gazdag相位外推而且必须改变，则下一次使用的是分步付立叶外推。如果当前使用的是分步付立叶外推而且必须改变，则下一次使用的是隐式有限差分外推。
在图2中的处理完成之后，下一个偏移间隔和下一次使用的外推类型已被确定。对相位相对误差的估计确定外推方法类型是否需要改变为更昂贵但更准确的方法。
所描述的第一个外推方法是Gazdag相移法，用于横向均匀介质。在常数速度的假定下，能从方程(2)导出一个单向波动方程。在对向下传播取空间和时间付立叶变换并选择正根之后，造成的频散关系由下式给出kz=ωc1-(cω)2(kx2+ky2)...(11)]]>相应的单向波动方程是∂∂zP(kx,ky,z,ω)=ikzP(kx,ky,z,ω)...(12)]]>
在波数域和空间域分别有如下解P(kx，ky，z+Δz，ω)＝P(kx，ky，z，ω)exp(ikzΔz) (13)P(x，y，z+Δz，ω)＝FT-1{P(kx，ky，z，ω)exp(ikzΔz)} (14)在方程(14)中，FT-1代表在两个水平波数κx和κy上的逆付立叶变换。
垂直波数κz能被分成两个贡献-一项归因于S零炮检距有关的传播，第二项归因于在零炮检距的相位变化。组合方程(11)和(13)，产生P(x,y,z+Δz,ω)=exp(iωΔzc)]]>·FT-1{P(kx,ky,z,ω)exp{(iωΔzc)[1-(cω)2(kx2+ky2)-1]}}...(15)]]>第一个指数项对应于“零炮检距”项，往往被称作薄透镜贡献或薄透镜项。第二个指数项对应于“时差”项，往往被称作绕射贡献或绕射项。
在实践中，首先应用绕射项，使用的外推间隔延伸于若干个输出深度采样上。然后，在成像期间，绕射贡献在每个深度采样被内插。在这一阶段，通过在每个深度采样的准确相移，还应用薄透镜项。对于输出深度采样间隔δz和绕射项偏移间隔Δz＝mδz，由方程(15)代表的解能在深度采样范围m≥n≥1上被分裂成一系列步骤。这一过程示于下面的流程图。这一过程描述图1的步骤104中所指的Gazdag相移外推。
图3显示的流程图说明本发明的进行Gazdag相移外推的方法的一个实施例的处理步骤。首先，在步骤301，为下一个偏移间隔Δz选择输出深度采样间隔δz，使对某一整数m有Δz＝mδz。
在步骤302，从叠加前地震数据集中选出在深度z的压强P(x，y，z，ω)。压强P(x，y，z，ω)要被外推到偏移间隔Δz中的深度采样间隔z+nδz，其中m≥n≥1。
在步骤303，对步骤302中选择的压强P(x，y，z，ω)应用付立叶变换，产生付立叶变换后的压强P(kx，ky，z，ω)，如下式所示P(kx，ky，z，ω)＝FT{P(x，y，z，ω)} (16)在步骤304，在深度采样间隔z+mδz，由来自步骤303的付立叶变换后的压强P(κx，κy，z，ω)计算出绕射项PDF(x，y，z+mδz，ω)。优选地，使用下面的方程计算绕射项PDF(x，y，z+mδz，ω)PDF(x,y,z+mδz,ω)=]]>FT-1{p(kx,ky,z,ω)exp{(iωmδzc)[1-(cω)2(kx2+ky2)-1]}}...(17)]]>在步骤305，使用在深度z的绕射项PDF(x，y，z，ω)和来自步骤304的在深度采样间隔z+mδz的绕射项PDF(x，y，z+mδz，ω)内插得到在深度采样间隔z+nδz的绕射项PDF(x，y，z+nδz，ω)。优选地，绕射项PDF(x，y，z+nδz，ω)使用下式内插得到PDF(x,y,z+nδz,ω)=PDF(x,y,z,ω)(m-nm)+PDF(x,y,z+mδz,ω)(nm)...(18)]]>在步骤306，薄透镜项被应用于在步骤305中计算出的在深度采样间隔z+nδz的绕射项PDF(x，y，z+nδz，ω)，产生在深度采样间隔z+nδz的压强P(x，y，z+nδz，ω)。优选地，压强P(x，y，z+nδz，ω)用下式产生P(x,y,z+nδz,ω)=PDF(x,y,z+nδz,ω)exp(iωnδzc)...(19)]]>
在图3中的处理完成之后，通过使用如方程(16)-(19)给出的Gazdag相移外推，向下外推出穿过下一个偏移间隔的压强。
所描述的第二个外推方法是分步付立叶方法，用于有弱横向不均匀性的介质。Gazdag外推能处理只随深度变化的介质。这是通过迭代地应用指数并在每个外推间隔内保持速度为常数实现的。在横向不均匀的情况中，方程(19)中给出的薄透镜项能被直接实现，无需修改，因为相移是在空间域中。然而，由于绕射解是在波数域中计算的，所以速度必须保持为常数。对于在小入射角时的快横向速度变化，或者可替换地对中等到大入射角时的小横向速度变化，使用外推间隔上的平均慢度是一个合理的近似。分步付立叶方法将初始的类似于Gazdag的外推与随后的在空间域慢度扰动上的外推组合在一起PGAZ(x,y,z+Δz,ω)=]]>FT-1{P(kx,ky,z,ω)exp[(iωΔzc0)1-(c0ω)2(kx2+ky2)]}...(20)]]>P(x，y，z+Δz，ω)＝PGAZ(x，y，z+Δz，ω)exp(iωΔzδ(1/c)) (21)其中δ(1/c)=1c(x,y,z)-1c0(z)...(22)]]>这样，解被分成两步，它们二者都是在频率域。
为了计算的准确性，方程(20)和(21)被重写，以把绕射贡献与薄透镜贡献分开。得到的公式由下式给出PDF(x,y,z+Δz,ω)=]]>FT-1{P(kx,ky,z,ω)exp{(iωΔzc0)[1-(c0ω)2(kx2+ky2)-1]}}...(23)]]>
和P(x,y,z+Δz,ω)=PDF(x,y,z+Δz,ω)exp(iωΔzc)...(24)]]>其中c＝c(x，y，z)，c0＝c0(z)。这些方程实质上能以与方程(16)-(19)中对均匀介质相同的方式来应用。唯一的差别是用于薄透镜计算的相移，只有方程(19)随横向坐标而改变。
在图1的步骤106中提到的用于实现分步付立叶外推方法的本发明的方法是恰如图3中那样应用方程(16)-(19)完成的，一个例外是在步骤305中方程(19)的薄透镜计算中使用空间上可变的速度。
在图4中的处理完成之后，通过使用分步付立叶外推，向下外推出穿过下一个偏移间隔的压强。
所描述的第三个外推方法是带有相位校正的隐式有限差分方法，用于具有中度到大的横向速度变化的介质。如果方程(11)被转换回到空间域，则可得到差分算子的平方根，它未被定义。如果需要空间域中的解，则需要对方程(11)中平方根的近似。连续分式可用于得到适合各种最大倾斜的显式外推器。数据中的倾斜越大，则为了准确代表它所需要的近似阶数就越高。这一途径产生隐式算法，它们费钱但稳定。
最小二乘分析可用于在广大的角范围内对有理近似进行最优化。准确到2n阶的近似由下式给出∂P∂z=(iωc)[1+Σj=1nαjS1+βjS]P...(25)]]>其中S=(c2ω)2(∂2∂x2+∂2∂y2)]]>
i=-1,]]>并且，在上文讨论的Lee和Suh(1985)的表1中可找到系数αi和βi。这样，方程(25)的第零阶、二阶和四阶近似可分别由下列各式给出∂P∂z=(iωc)P...(26)]]>∂P∂z=(iωc)[1+α1S1+β1S]P...(27)]]>其中α1＝0.478242060，β1＝0.376369527以及∂P∂z=(iωc)[1+α1S1+β1S+α2S1+β2S]P...(28)]]>其中α1＝0.040315157，β1＝0.873981642，α2＝0.457289566以及β2＝0.222691983。
方程(26)、(27)和(28)分别指“5”、“65”和“80”度近似。更多的系数集由Curtis Ober、Rob Gjertsen以及DavidWomble在“Salvo用于复杂地质构造的地震偏移代码(SalvoASeismic Migration Code for Complex Gelogy)”(1999，为ACTI组合，Sandia National Lab，Albuquergue，N.M.准备的报告)中给出，它们对应于在角度达15、45、60、75、87和90度的范围上的最优化。能由单个(α，β)对近似的最宽角孔径是75度，这是本发明选择的默认的优选可选项。
如果通过乘以每个分母使在有理表达式中提出公因子的话，则差分算子的作用变得显而易见。这对于二阶近似造成形如2/x2(P/z)和2/y2(P/z)的导数，对四阶近似造成形如4/x4(P/z)，4/x2y2(P/z)以及4/y4(P/z)的导数。表示这些高阶导数是麻烦的，需要长的算子。
能通过相继求解方程(27)或(28)中的每一项计算出更高阶的解。这有时被称作“Marchuk分解”。这样，对于四阶解，要求解如下方程组∂P0∂z=(iωc)P0...(29)]]>∂P1∂z=(iωc)[α1S1+β1S]P1...(30)]]>∂P2∂z=(iωc)[α2S1+β2S]P2...(31)]]>其中P0(x，y，z，ω)＝P(x，y，z，ω)P1(x，y，z，ω)＝P0(x，y，z+Δz，ω)P2(x，y，z，ω)＝P1(x，y，z+Δz，ω)P(x，y，z+Δz，ω)＝P2(x，y，z+Δz，ω)方程(29)是惯用的深度偏移的薄透镜项，而方程(30)和(31)是对时间偏移和深度偏移共同的绕射校正。
方程(30)和(31)仍是产生形如2/x2(P/z)和2/y2(P/z)的导数。这不是一个严重问题。更严重的关切是通过公因子S存在于方程(30)和(31)每一个中的导数2/x2和2/y2。这导致水平面内解的耦合并增加结果差分方程的复杂性和成本。通过在x方向和y方向分别求解有理表达式能消除这一耦合。方程(30)和(31)能被表示为如下通用形式
∂Pj∂z=(iωc)[αjS1+βjS]Pj...(32)]]>忽略形如4/x2y2的因子，方程(32)能通过如下系统近似∂Pj∂z=(iωc)[αjSx1+βjSx]Pj...(33)]]>和∂Pj∂z=(iωc)[αjSy1+βjSy]Pj...(34)]]>其中Sx=(c2ω)2∂2∂x2]]>和Sy=(c2ω)2∂2∂y2]]>现在方程(33)中的差分算子只是x的函数，因此这一方程能对y为常数的每个数据切片求解。类似地，方程(34)中的差分算子只是y的函数，因此这一方程方程能对x为常数的每个数据切片求解。
上文讨论的分步忽略了沿副线方向(subline)和交叉线方向导数间的交叉项。对于传播方向靠近正交坐标轴的情况，这一近似是有效的。
组合方程(25)、(33)和(34)，该分步外推有如Lee和Suh(1985)给出的形式∂P∂z=(iωc)[1+Σj=1nαjSy1+βjSx+Σj=1nαjSy1+βjSy]P...(35)]]>
通过将方程(35)的结果与使用精确频散关系得到的结果进行简单地比较，能定量得到与分步和旁轴近似以及网格频散(griddispersion)关联的总误差。这一构想构成了为调节波场相位可在外推期间间歇应用的校正方法的基础。
在利用等效于算子的波数重写方程(35)并取精确外推器和方程(35)之差以后，沿径向未校正有限差分计算的相位校正由下式(Li，1990)给出Γ(kx,ky,ω)=ωc1-(c2ω)2(kx2+ky2)]]>-ωc[1-Σj=1nαjk^x2(c/2ω)2-βjk^x2-Σj=1nαjk^y2(c/2ω)2-βjk^y2],...(36)]]>其中k^x2=2-2cos(kxΔx)(Δx)2]]>和k^y2=2-2cos(kyΔy)(Δy)2]]>方程(36)是对波数写出的，因为未在平方根中定义差分算子。带“帽子”的波数代表对真实波数κx和κy的有限差分近似，其中Δx和Δy是沿交叉线和副线方向的离散化间隔。为了使方程(36)有意义，在其中不能存在任何横向变化的量。这样，必须假定对于该外推间隔速度c是常数。通过对波场P(x，y，z，ω)进行2D FFT(快速付立叶变换)然后进行相移以应用校正，这要求解如下方程∂∂zP(x,y,z,ω)=iΓ(kx,ky,ω)P(kx,ky,z,ω)...(37)]]>需要最后的2D逆FFT，以使解返回到空间域。对于大的横向速度变化，可在选定的参考速度应用校正。然后所得到的波场能对速度内插。
应该理解，上述只是对本发明具体实施例的详细描述，根据这里公开说明的内容能对所公开的实施例做出大量改变、修改和替换并不离开本发明的范围。所以，上文的描述不意味着对本发明范围的限制。相反，本发明的范围只是由所附权利要求及其等效物确定。
权利要求
1.一种对叠加前的地震数据进行向下外推的方法，包含选择一个叠加前地震数据集；确定该地震数据集的偏移间隔；为该偏移间隔选择最大误差判据；计算相位的最大相对误差，作为频率、传播角和该偏移间隔中速度相对变化的函数；将相位的最大相对误差与最大误差判据进行比较；以及根据相位的最大相对误差与最大误差判据的比较结果确定该偏移间隔中使用的外推类型。
2.权利要求1的方法，其中的外推类型是从一组外推类型中选出的，该组外推类型包含Gazdag相移外推、分步移动付立叶外推以及隐式有限差分外推。
3.权利要求1的方法，其中计算相位的最大相对误差的步骤包含在该偏移间隔中选择频率ω；在该偏移间隔中识别最大波数κmax；在该偏移间隔中识别最小速度cmin；使用频率ω、最大波数κmax和最小速度cmin，确定该偏移间隔中的最大传播角θmax；在该偏移间隔中选择参考速度co；使用最大传播角θmax和参考速度co，计算相位的最大相对误差Δφmax。
4.权利要求3的方法，其中识别出最大波数κmax，以便在半径κ＝κmax的圆圈内包括功率的95％。
5.权利要求3的方法，其中的最小速度cmin是通过对下一个偏移步长扫描速度切片以识别最小速度值而识别出来的。
6.权利要求3的方法，其中最大传播角θmax是由如下方程确定的θmax=sin-1(kmax·cminω)]]>其中κmax是最大波数，cmin是最小速度，ω是频率。
7.权利要求3的方法，其中参考速度co是由如下方程确定的c0=cminsin(θmax)]]>其中cmin是最小速度，θmax是最大传播角。
8.权利要求3的方法，其中用于分步付立叶外推的相位最大相对误差Δφmax，是用如下方程确定的ΔΦΦ(c(x,y))=(ωΔzc0)·[(c(x,y)c0)2-sin2(θ)+(1-cos(θ)-c(x,y)c0)]]]>其中φ是相位，Δφ是在该偏移间隔上的相位误差，θ是相对于垂直方向的传播角，c(x，y)是速度，co是参考速度。
全文摘要
通过首先确定一个叠加前地震数据集中的偏移间隔，该叠加前地震数据集被向下外推。为该偏移间隔选择一个最大误差判据。作为频率、传播角和该偏移间隔中速度相对变化的函数，计算出相位的最大相对误差。该相位最大相对误差与最大误差判据进行比较。根据相位最大相对误差与最大误差判据的比较结果确定该偏移间隔中使用的外推类型。该外推类型是从一组外推类型中选出的，该组外推类型包括Gazdag相移外推、分步移动付立叶外推和隐式有限差分外推。
文档编号G01V1/36GK1782738SQ20041008705
公开日2006年6月7日 申请日期2004年10月22日 优先权日2003年10月23日
发明者S·M·凯利 申请人:Pgs美洲公司