CN101937100A - 一种叠前深度偏移方法 - Google Patents

Abstract

本发明涉及地震勘探领域，具体的是指一种叠前深度偏移方法。该方法包括：对地震数据进行偏移速度分析；通过射线追踪计算走时、弧长、出射角、入射角；进行偏移孔径计算；利用Kirchhoff积分法矢量偏移公式进行叠前深度偏移。本方法将偏移速度分析、偏移孔径选取、射线追踪、Kirchhoff积分公式结合于一体，不需要进行波场分离，并实现多分量同时偏移，实现转换波的准确归位。从实际勘探数据来看，采用本发明的叠前深度偏移方法，不仅地表成像效果好，成像分辨率高，而且转换波成像剖面上的深层反射同相轴的连续性提高，构造更加清晰。

Description

一种叠前深度偏移方法

技术领域

本发明涉及地震勘探领域，具体的是指一种叠前深度偏移方法。

背景技术

偏移是地震数据处理的最后一个步骤，也是最重要、关键的一个过程，对于复杂的地下介质而言更是如此。

由于地下介质层位的非水平性，造成了地震叠加剖面并不能反映真实的层位形态，因此需要对地震数据在叠加的基础上进一步进行偏移处理，校正由于地层倾角存在所产生的层位反射的假象时移，使地下层位的反射归位，反映地层的真实起伏形态和产状。

偏移处理方法分为两类，一类在叠加的基础上进行，称为叠后偏移；一类在叠加前进行偏移处理，称为叠前偏移。叠前偏移精度要高于叠后偏移，对复杂绕射波的归位能力较好；但叠前偏移处理需要的计算工作量很大，随着近十多年来计算机技术的提高而得到了广泛的发展和应用。其中叠前偏移处理按照偏移方法的不同，主要是方法所在的时空域的差异，分为叠前时间偏移和叠前深度偏移。叠前深度偏移由于输出剖面直接反映真实空间深度域的地下介质形态而在应用领域得到更大的关注。

目前业界关于多分量地震数据的叠前深度偏移研究非常少，相对成熟应用的还是叠后时间偏移。多分量地震数据的叠前偏移处理是目前地震勘探领域的攻关热点和难点，而且当前研究与尝试应用还大多局限于叠前时间偏移，也远没有走入生产实践应用中。而且大部分研究还局限于在纵波偏移的基础上，沿用与纯纵波偏移相同的思路进行C波或等效转换波的偏移处理，将多分量地震信号间的空间矢量关系完全割裂开来，进行单纯波的偏移处理，既损失了波场中蕴含的丰富地下介质信息，成像效果也较差。

英国爱丁堡大学的李向阳教授所领导的研究小组(Hengchang Dai，and Xiang-Yang Li.Effect of errors in the migration velocity model of PS-converted waves on traveltime accuracy in prestack Kirchhoff time migration in weak anisotropic media.Geophysics，2008，73(5)：S195-S205)(Xiang-Yang Li，Jianxin Yuan.Converted-wave moveout and conversion-point equations in layered VTI media：Theory and applications.Journal of Applied Geophysics，2003，34：297-318)初步试验成功了多分量转换波的叠前时间偏移技术，在一些地区的多分量地震勘探中初步试验和应用，但尚需进一步攻关完善。此外，戴和李针对弱各向异性介质转换横波叠前时间偏移中速度对偏移结果的影响进行数值模拟和误差的讨论。他们分别讨论了转换横波速度、垂向速度比、有效速度比、各向异性参数的精度对旅行时计算误差的贡献，对于时间域转换波的成像提高具有较好的参考价值；且一些思路可以为转换波的叠前深度偏移所借鉴。

上世纪80年代王妙月在Kuo(Kuo，J.T.，and Dai，T.，1984，Kirchhoff elastic wave migration for the case of non-coincident source and receiver：Geophysics，49(8)，1223-1238)的基础上最早提出了矢量偏移的思想，并推出了弹性波Kirchhoff积分公式(秦福浩、郭亚曦、王妙月，弹性波克希霍夫积分偏移法，地球物理学报，1988，31(5)：577-587)。这是一种均匀各向同性介质假设下的弹性波的多分量同时偏移方法，通过PP、PS波的走时将二分量的数据逆时延拓到地下成像点，得到每个成像点刚刚被激发时的三个分量上的振动形式。通过数值模拟的合成地震记录，他们对该方法进行了测试，获得相对较好的结果，在这一方向进行了初步的探索，但截至目前没有进一步的发展。

该方法有以下两个不足：1)得到的是总的二个分量振动能量，没有实现波场分离；2)该方法得到的偏移结果是某一炮点的能量传播到对于地下反射界面时激发的震动。不同炮点对于地下同一反射点激发时的振动矢量不相同，因此不能用于多炮记录的叠加。而且他们的研究是一种近似的叠前深度偏移，即在偏移过程中仍然在时间域进行，只是后来进行时-深转换把输出的偏移剖面转换成深度剖面。而且对于走时的计算他们只停留在单界面直射线追踪这一层面，对于不同的深度成像也没有考虑到偏移孔径的差异，偏移速度的获得相对简单，只是简单借用叠加速度然后进行层速度的换算，因而效果并不理想。

公开号为US 2009/0257308A1的发明专利保护偏移速度分析领域的一种技术，该技术基于动校正误差来不断的更新速度模型，从而使得最终的速度模型逼近于地下介质的真实速度。该专利是基于剩余曲率分析方法的速度分析技术，而且该技术基于小炮检距假设和水平层位介质假设，即要求剩余旅行时可以近似为双曲线或抛物线。而当速度横向变化剧烈或大炮检距情况下，该技术所基于的算法近似不合适的；而且该专利技术只针对常规纵波数据和纵波速度进行处理分析。

亓雪冬等(亓雪冬、贾光华、仝兆歧，Kirchhoff叠前时间偏移处理技术及应用，物探化探计算技术，2009，31(2)：126-130)研究的是针对常规纵波资料处理的叠前时间偏移技术及其应用，这一技术在油气地震勘探领域已是成熟应用的、商业化的生产技术。

张猛等(张猛、孟祥宾、匡斌等，积分法叠前深度偏移技术在BS6地区的应用，勘探地球物理进展，2009，32(1)：48-55)介绍的是成熟的、商业化应用的常规纵波地震数据的叠前深度偏移技术及其软件对胜利油田某区块的地震数据进行处理试验和应用，并与常规纵波的叠后偏移(时间)进行成像效果的对比，从而说明叠前深度偏移的优势和应用效果。

叶月明等(叶月明、李振春、仝兆歧，基于稳定成像条件的保幅叠前深度偏移，石油地球物理勘探，2009，44(1)：28-32)主要针对常规纵波叠前深度偏移处理中阻尼型反褶积处理存在的问题，将分母趋于零时所引起的计算不稳定现象通过高斯平滑处理，提高偏移处理的计算稳定性，并同时达到部分压噪的目的。

上述有关技术所针对的地震数据都是单分量、常规纵波数据，无一涉及多分量地震数据的处理，所涉及的技术与理论方法也不适用于三分量地震数据的处理。

公开号为CN 101598805A的发明专利只是针对纵波与转换横波传播时间的差异，通过速度比扫描和层位追踪的方法，将同一层位在纵波与转换横波剖面的不同反射相位进行标定、拾取，从而将各对应相位的纵横波速度比确定下来，进而实现将转换横波剖面按照纵波传播时间压缩，使得在相同的纵波传播时间标度上，方便纵波与转换横波的层位对比。其保护的内容并不涉及三分量地震数据的成像技术。

公开号为CN 1797031A的发明专利针对常规纵波数据的叠前深度偏移处理中由于海量数据处理耗时、计算效率低的问题，通过相同相位炮的匹配和组合，实现多炮合并一炮，从而减少偏移处理的炮集数，在不损伤偏移效果的基础上实现偏移处理效率的提高。其保护的内容不涉及三分量地震数据的成像技术。

公开号为CN 101419292A的发明专利保护的内容是针对三分量采集的地震数据，通过叠前预处理、静校正、速度分析、动校正、叠加这样一个常规的处理流程对转换横波进行处理成像。该专利保护的内容在三分量地震数据处理领域一般又称为常规叠后处理方法，仅涉及时间偏移技术。

公开号为US 2004/0117123A1的发明专利保护的内容是针对常规纯纵波数据的叠前时间偏移处理技术，而且该技术是在频率-波数域实现的。该专利是将地震数据变换到频率-波数域，然后根据旅行走时图对给定空间任一点的射线参数进行相移的线性拟合使得叠前时间偏移算法适应速度的横向变化剧烈的复杂介质成像问题。

公开号为US 2010/0054082A1的发明专利保护的内容不涉及具体的偏移算法，更与多分量叠前深度偏移无关。在任何偏移过程中均涉及波场的模拟，其并不是通过射线追踪走时计算完成的，而是通过声波方程的有限差分模拟实现。该专利主要针对叠前逆时偏移处理中需要存储空间任一点的波场而造成偏移处理中内存与硬盘空间的大量消耗问题，且使得偏移效率低下，为此提出一种节省计算机内存与硬盘空间的技术改进措施，即在正演波场的模拟中只保存阻抗界面点上的波场值，其他空间点不参与计算，从而大大改进了偏移处理的速度，降低了计算资源的消耗。

叠前深度偏移是对大地层倾角和速度横向变化剧烈的复杂构造成像的有效工具。但是叠前深度偏移依赖更精确的速度模型。常规的速度建模方法无法达到足够的精度。同时，因为叠前深度偏移对速度模型十分敏感，所以它作为一种有力的速度分析工具，发展出了两种偏移速度分析方法：深度聚焦分析(DFA，depth-focusing analysis)和剩余曲率分析(RCA，residual-curvature analysis)。要进行偏移速度分析就必须解决两个问题：1)建立一个评价速度是否准确的标准；2)如何进行速度更新。

深度聚焦分析的原理是根据叠加能量测量速度误差。当偏移深度和聚焦深度相同时，速度达到要求。剩余曲率分析是根据剩余旅行时差测量速度误差。当不同炮检距的成像深度相同时，速度达到要求。

因为速度分析通常是基于一个粗糙的初始迭代速度，因此一个健壮的速度分析方法在减少迭代次数的同时还要保证速度的收敛性。传统的速度分析方法，通常是基于以下三个假设：①层状均匀介质；②小炮检距；③水平地层。

MacKay和Abma(MacKay S.and Abma R.Imaging and velocity estimation with depth-focusing analysis.Geophysics.1992，57：1608-1622)的DFA方法用到了全部以上三个近似。

在RCA中，速度分析时所采用的偏移方法不同，基于的假设也不同：1.共炮点或共检波点偏移：Al-Yahay(Al-Yahya K.1989.Velocity analysis by iterative profile migration，Geophysics，1989，54：718-729)提出的方法基于以上三个假设；Lee和Zhang(Lee W.and Zhang L.Residual shot profile migration，Geophysics，1992，57：815-822)的方法将第三个假设拓展到小倾角。2.共炮检距偏移：Deregowski(Deregowsiki S.M.Common-offset migration and velocity analysis.First Break，1990，8(6)：225-234)的方法使用了前两个假设。

以上的方法都基于小炮检距假设，在小炮检距情况下，剩余旅行时可以近似为双曲线或抛物线。当速度横向变化剧烈时，使用这个近似是不合适的。因此，Liu(Liu Z.and Bleistein N.Migration velocity analysisi：Theory and an iterative algorithm，Geophysics，1995，60：142-153)在Lafond 和Levander(Lafond C.F.and Levander A.R.Migration moveout analysis and depth focusing，Geophysics，1993，58：91-100)的基础上提出了一种新的RCA方法——基于深度聚焦的共成像点速度分析方法，这种方法适用于任何炮检距、倾角和速度变化。

发明内容

本发明所要解决的技术问题是针对多分量地震数据，提供一种基于矢量波场的三分量叠前深度偏移方法，不需要进行波场分离，并实现多分量偏移，实现转换波的准确归位。

本发明解决上述技术问题的技术方案如下：一种叠前深度偏移方法，包括以下步骤：

步骤1、对地震数据进行偏移速度分析；

步骤2、通过射线追踪计算走时、弧长、出射角、入射角；

步骤3、进行偏移孔径计算；

步骤4、利用Kirchhoff积分法矢量偏移公式进行叠前深度偏移。

本发明的有益效果是：本方法将偏移速度分析、偏移孔径选取、射线追踪、Kirchhoff积分公式结合于一体，不需要进行波场分离，并实现多分量偏移，实现转换波的准确归位。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，步骤1中所述偏移速度分析采用基于深度聚焦的共成像点速度分析方法。

采用上述进一步方案的有益效果是，基于深度聚焦的共成像点速度分析方法是在剩余曲率分析方法上修改而成，其收敛速度快，适用于各种速度模型和炮检距。

进一步，所述步骤1包括如下步骤：

步骤a1、初始地质模型，并进入步骤a2；

步骤a2、建立初始速度模型，并进入步骤a3；

步骤a3、进行基于目标层位底界面的PSDM(pre-stack depth migration，叠前深度偏移)，并进入步骤a4；

步骤a4、进行偏移速度分析，并进入步骤a5；

步骤a5、判断是否建立了完整的速度模型，如果为“是”则进入步骤2，如果为“否”则进入步骤a6；

步骤a6、进行基于目标层位底界面的PSDM，并进入步骤a7；

步骤a7、拾取目标层底界面，并进入步骤a8；

步骤a8、更新速度模型，将目标层移至下一层位，并进入步骤a3。

其中，步骤a5中所述完整的速度模型为：从浅到深各层位齐全的以速度、密度、层位埋深作为参数的地震地质模型。

进一步，步骤2中，采用直射线追踪法或者采用基于费马原理的最小走时射线追踪法。

采用上述进一步方案的有益效果是，采用直射线追踪法计算速度快，采用基于费马原理的最小走时射线追踪法可以避免射线盲区。

进一步，步骤4中采用适用于非均匀各向异性介质的弹性波Kirchhoff积分法矢量偏移公式进行叠前深度偏移。

采用上述进一步方案的有益效果是，适用于非均匀各向异性介质的弹性波Kirchhoff积分法矢量偏移公式，具有波场分离的特性，可以做到偏移的同时分离PP波和PS波；偏移结果可分别输出PP波剖面和PS波剖面，利于分别分析PP波和PS波的AVO(Amplitude Versus Offset，振幅随偏移距的变化)特征；使用了弹性波动方程的格林函数，提高算法的保幅性；适用于非均匀各向异性介质。

本发明在王妙月等人(秦福浩、郭亚曦、王妙月，弹性波克希霍夫积分偏移法，地球物理学报，1988，31(5)：577-587)研究的基础上推出新的偏移公式，并针对其所存在的两个不足进行改进。在偏移过程中直接输出PP波和PS波，在进行多炮叠加的时候实现了相干加强的效果，得到PP和PS波的剖面。

附图说明

图1为本发明叠前深度偏移方法的原理图一；

图2为本发明叠前深度偏移方法的原理图二；

图3为本发明叠前深度偏移方法的原理图三；

图4为本发明叠前深度偏移方法的流程图；

图5为本发明中偏移速度分析的实现流程；

图6A为采用基于共CCP(Common Converted Point，共转换点)道集的叠后时间偏移方法对大庆一条2D3C地震数据进行偏移之后的PP波效果图；

图6B为采用基于共CCP道集的叠后时间偏移方法对大庆一条2D3C地震数据进行偏移之后的PS波效果图；

图7A为采用本发明方法对大庆的一条2D3C地震数据进行的偏移处理之后的PP波效果图；

图7B为采用本发明方法对大庆的一条2D3C地震数据进行的偏移处理之后的PS波效果图；

图8为本发明叠前深度偏移方法的模拟数据模型；

图9A为图8的数据模型的地震数据Z分量记录图；

图9B为图8的数据模型的地震数据X分量记录图；

图10为图8的数据模型中采用本发明的叠前深度偏移方法得到的纵波偏移深度剖面结果；

图11为图8的数据模型中采用本发明的叠前深度偏移方法得到的转换横波偏移深度剖面结果；

图12为图8的数据模型中采用声波方程偏移并随倾角开孔径得到的纵波偏移深度剖面结果；

图13为图8的数据模型中采用声波方程偏移并随倾角开孔径得到的转换横波偏移深度剖面结果；

图14为图8的数据模型中采用声波方程偏移并固定孔径得到的纵波偏移深度剖面结果；

图15为图8的数据模型中采用声波方程偏移并固定孔径得到的转换横波偏移深度剖面结果。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

方法原理

如图1所示，均匀各向同性介质的三维矢量波场延拓方程为：

$u_n(\stackrel{\→}{x},t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\text{\∫}}}\underset{-\∞}{\overset{\∞}{\∫}}\left\{\begin{array}{c}[(1-{2v}^2){\mathrm{\δ}}_{m3}{\hat{r}}_n+{2v}^2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3]\×{\stackrel{\·}{u}}_m({\stackrel{\→}{x}}^{\′},t+t_P)/\left({\mathrm{rV}}_P\right)\\ +({\mathrm{\δ}}_{\mathrm{mn}}{\hat{r}}_3+{\mathrm{\δ}}_{n3}{\hat{r}}_m-2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3)\×{\stackrel{\·}{u}}_m({\stackrel{\→}{x}}^{\′},t+t_S)/\left({\mathrm{rV}}_S\right)\\ +\left({\hat{r}}_m{\hat{r}}_n\right)T_{3m}({\stackrel{\→}{x}}^{\′},t+t_P)/\left({\mathrm{\ρrV}}_P^2\right)\\ +({\mathrm{\δ}}_{\mathrm{mn}}-{\hat{r}}_m{\hat{r}}_n)T_{3m}({\stackrel{\→}{x}}^{\′},t+t_S)/\left({\mathrm{\ρrV}}_S^2\right)\end{array}\right\}{\mathrm{dx}}_1^{\′}{\mathrm{dx}}_2^{\′}---\left(1\right)$

其中

表示质点位移，单位为m/s；

表示质点速度，单位m；表示应力，单位为N/m²；m，n＝1，2，3分别表示x，y，z坐标，单位m；

为单位矢径，

分别表示射线出射角的单位向量在m、n坐标上的投影；t＝t(x|x′)表示震源到地下绕射点的走时，单位为s；t_p、t_s分别代表绕射点到接收点的PP波、PS波走时，单位为s；r表示绕射点到接收点的距离，单位为m；V_p、V_s分别表示纵、横波速度，单位为m/s；v＝V_s/V_p表示纵横波速度比。

由于自由表面边界条件，地表应力为零，方程(1)的后两项可以舍去，有

$u_n(\stackrel{\→}{x},t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}\underset{-\∞}{\overset{\∞}{\∫}}\left\{\begin{array}{c}[(1-{2v}^2){\mathrm{\δ}}_{m3}{\hat{r}}_n+{2v}^2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3]\×{\stackrel{\·}{u}}_m({\stackrel{\→}{x}}^{\′},t+t_P)/\left({\mathrm{rV}}_P\right)\\ +({\mathrm{\δ}}_{\mathrm{mn}}{\hat{r}}_3+{\mathrm{\δ}}_{n3}{\hat{r}}_m-2{\hat{r}}_m{\hat{r}}_n{\hat{r}}_3)\×{\stackrel{\·}{u}}_m({\stackrel{\→}{x}}^{\′},t+t_S)/\left({\mathrm{rV}}_S\right)\end{array}\right\}{\mathrm{dx}}_1^{\′}{\mathrm{dx}}_2^{\′}---\left(2\right)$

将此方程展开，将三个分量上的能量分别投影到PP分量(3)式和PS分量(4)式上，如图2所示

u_PP＝u₁cosθsinφ+u₂sinθsinφ+u₃cosφ(3)

u_PS＝u₁cosθcosφ+u₂sinθcosφ+u₃sinφ(4)

将(2)式展开后带入(3)式、(4)式，得

$u_{\mathrm{PP}}(x,y,z,t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}\underset{-\∞}{\overset{\∞}{\∫}}\left[\begin{array}{c}\frac{{2v}^2{\hat{r}}_x{\hat{r}}_z}{{\mathrm{rV}}_P}\·{\stackrel{\·}{u}}_x(x^{\′},y^{\′},0,t+t_P)+\frac{{2v}^2{\hat{r}}_y{\hat{r}}_z}{{\mathrm{rV}}_P}\·{\stackrel{\·}{u}}_y(x^{\′},y^{\′},0,t+t_P)\\ +\frac{1-{2v}^2+{2v}^2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_P}\·{\stackrel{\·}{u}}_z(x^{\′},y^{\′},0,t+t_P)\end{array}\right]{\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}---\left(5\right)$

$u_{\mathrm{PS}}(x,y,z,t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}\underset{-\∞}{\overset{\∞}{\∫}}\left[\begin{array}{c}\frac{1-2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_S}\·{\stackrel{\·}{u}}_x(x^{\′},y^{\′},0,t+t_S)+\frac{1-2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_S}\·{\stackrel{\·}{u}}_y(x^{\′},y^{\′},0,t+t_S)\\ +\frac{2({\hat{r}}_x{\hat{r}}_z+{\hat{r}}_y{\hat{r}}_z)}{{\mathrm{rV}}_S}\·{\stackrel{\·}{u}}_z(x^{\′},y^{\′},0,t+t_S)\end{array}\right]{\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}---\left(6\right)$

在二维多分量勘探中，通常使用涨缩震源，因此理论地震记录上只存在P-P波和P-SV波。若令垂向为z轴方向，测线方向为x轴方向，则y分量上无波场信息。将(5)式、(6)式简化为二维弹性波Kirchhoff积分法矢量偏移公式。

$u_{\mathrm{PP}}(x,z,t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}[\frac{{2v}^2{\hat{r}}_x{\hat{r}}_z}{{\mathrm{rV}}_P}\·{\stackrel{\·}{u}}_x(x^{\′},t+t_P)+\frac{1-{2v}^2{\hat{r}}_x{\hat{r}}_x}{{\mathrm{rV}}_P}\·{\stackrel{\·}{u}}_z(x^{\′},t+t_P)]{\mathrm{dx}}^{\′}---\left(7\right)$

$u_{\mathrm{PS}}(x,z,t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}[\frac{1-2{\hat{r}}_z{\hat{r}}_z}{{\mathrm{rV}}_S}\·{\stackrel{\·}{u}}_x(x^{\′},t+t_S)+\frac{2{\hat{r}}_x{\hat{r}}_z}{{\mathrm{rV}}_S}\·{\stackrel{\·}{u}}_z(x^{\′},t+t_S)]{\mathrm{dx}}^{\′}---\left(8\right)$

然而，实际介质是非均匀和各向异性的，地下绕射点到地表接收点的PP波与PS波的射线路径不再是直线，两者的射线路径、走时和出射角均不相同。如图3所示，其中虚线表示S波射线路径，实线表示P波射线路径。因此需要对均匀各向同性介质偏移公式进行改进。走时t、t_p、t_s、路径l_p、l_s、出射角θ_p、θ_s和等效纵横波速度V_p\V_s通过非均匀各向异性介质的射线追踪和速度分析分别求得。

因此可以将(7)式、(8)式进行改进得

$u_{\mathrm{PP}}(x,z,t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}[\frac{{1-2v}^2{\sin }^2{\mathrm{\θ}}_P}{{l_{P}V}_P}\·{\stackrel{\·}{u}}_z(x^{\′},t+t_P)+\frac{{2v}^2\sin {\mathrm{\θ}}_P{\cos \mathrm{\θ}}_P}{l_P{V}_P}\·{\stackrel{\·}{u}}_x(x^{\′},t+t_P)]{\mathrm{dx}}^{\′}---\left(9\right)$

$u_{\mathrm{PS}}(x,z,t)=\frac{1}{4\mathrm{\π}}\underset{-\∞}{\overset{\∞}{\∫}}[\frac{1-2\mathrm{co}{s}^2{\mathrm{\θ}}_S}{{l_{S}V}_S}\·{\stackrel{\·}{u}}_x(x^{\′},t+t_S)+\frac{2\cos {\mathrm{\θ}}_S\sin {\mathrm{\θ}}_S}{{l_{S}V}_S}\·{\stackrel{\·}{u}}_z(x^{\′},t+t_S)]{\mathrm{dx}}^{\′}---\left(10\right)$

以上(9)式、(10)式两个方程就是适用于非均匀各向异性介质的弹性波Kirchhoff积分法矢量偏移公式。

对比目前基于声波方程的均匀各向同性Kirchhoff积分偏移算法(11)式：

$u(x,y,t)=-\frac{1}{2\mathrm{\π}}\underset{A}{\text{\∫\∫}}\frac{\cos \mathrm{\θ}}{{\mathrm{lv}}_{\mathrm{mig}}}\stackrel{\·}{u}(x^{\′},y^{\′},0,t+\frac{1}{v_{\min }}){\mathrm{dx}}^{\′}{\mathrm{dy}}^{\′}---\left(11\right)$

(9)式、(10)式的非均匀各项异性介质中的矢量弹性波场Kirchhoff积分偏移算法具有以下几个优势：

1)具有波场分离的特性，可以做到偏移的同时分离PP、PS波；

2)偏移结果分别输出PP波剖面和PS波剖面，利于分别分析PP波和PS波的AVO特征；

3)使用了弹性波动方程的格林函数，提高算法的保幅性；

4)适用于非均匀各向异性介质。

根据上述方法原理，本发明叠前深度偏移方法实现步骤如下：

步骤1、对地震数据进行偏移速度分析；

步骤2、通过射线追踪计算走时、弧长、出射角、入射角；

步骤3、进行偏移孔径计算；

步骤4、利用Kirchhoff积分法矢量偏移公式进行叠前深度偏移。

其中，步骤1中偏移速度分析采用基于深度聚焦的共成像点速度分析方法；步骤2中，采用直射线追踪法或者采用基于费马原理的最小走时射线追踪法；步骤4中采用适用于非均匀各向异性介质的弹性波Kirchhoff积分法矢量偏移公式进行叠前深度偏移。

本发明叠前深度偏移方法实现步骤的流程参照图4所示。

其中矢量叠前深度偏移算法的流程如下：

A.对SEG-Y数据(地震数据一般以地震道为单位进行组织，采用SEG-Y文件格式存储。SEG-Y格式是由SEG(Society of Exploration Geophysicists，地球物理学勘探协会)提出的标准磁带数据格式之一，它是石油勘探行业地震数据的最为普遍的格式)的所有道循环；

a)输入一道SEG-Y数据，并通过道头读出炮点检波点坐标(输入时判断是X分量还是Z分量，不同的分量将采用不同的走时和格林函数)；

b)将粗糙的走时表插值为与偏移输出剖面相同的网格；

c)对于地下所有成像点循环；

●判断是否在PP波偏移孔径内；

●如果在孔径范围内，则计算该点的格林函数，读出该点PP波走时所对应的能量叠加到成像点的PP分量；

●判断是否在PS波偏移孔径内；

●如果在孔径范围内，则计算该点的格林函数，读出该点PS波走时所对应的能量叠加到成像点的PS分量；

d)对成像点循环结束；

B.SEG-Y数据体循环结束，将偏移结果写入文件，偏移结束。

其总的算法实现流程如表1所示。

表1矢量叠前深度偏移算法的流程表

在本发明叠前深度偏移方法的实现过程中，偏移孔径的选择和偏移速度的获得是关系到偏移效果和计算效率的两个主要控制因素，本发明对其进行了如下改进。

偏移孔径

选择合适的偏移孔径可以提高Kirchhoff积分法矢量偏移公式进行叠前深度偏移的成像质量。偏移孔径过大会降低信噪比；反之则会遗漏有效信息。PP波、PS波在介质中的射线路径不同，因此需要采用不同的孔径。

孔径的选择是对于成像结果(保幅性、精确性)影响最大的因素。现有的文献中所提出的Fresnel带偏移、高斯束偏移、有限孔径偏移事实上都是通过射线理论寻找界面反射能量的真实位置，在这个点的周围开辟孔径。在这个过程中需要重复进行三角函数和反三角函数的运算，需要很大的计算量。而且偏移过程中，需要一次性读入单炮数据。因此在实际生产中，都是基于水平界面假设，在CMP(Common Middle Point，共中心点)附近开孔径，遇到倾角时，只是简单的将孔径扩大。

本发明采用的方法原理与上述方法相同，但是在算法上进行了改进，用中心点坐标X_m和偏移距h分别代替炮点和接收点的坐标，避免了三角函数和反三角函数的计算，提高了运算效率，并且可以适用于基于输入道的偏移方法。

偏移速度分析

本发明中的偏移速度分析采用基于深度聚焦的共成像点速度分析方法，其基础是RCA方法(Liu Z.and Bleistein N.Migration velocity analysisi：Theory and an iterative algorithm，Geophysics，1995，60：142-153)，该方法适用于任何炮检距、倾角和速度变化。在共成像点道集中，成像深度z是炮检距h的函数。当偏移速度等于地下真实速度时，成像深度与炮检距无关，以此为准则进行速度更新，更新公式如下：

$\mathrm{\δ\λ}=-\frac{{\mathrm{\Σ}}_{k=1}^K{\mathrm{\Σ}}_{j=1}^m(g_j^{\left(k\right)}-\stackrel{\‾}{{\hat{g}}^{\left(k\right)}})(z_j^{\left(k\right)}-\stackrel{\‾}{{\hat{z}}^{\left(k\right)}})}{{\mathrm{\Σ}}_{k=1}^K{\mathrm{\Σ}}_{j=1}^m{(g_j^{\left(k\right)}-\stackrel{\‾}{{\hat{g}}^{\left(k\right)}})}^2}$

其中，δλ为速度更新步长，新的迭代速度为V+δλ；K代表共成像点，M表示炮检距；其中，

$g(x,h)=\frac{t_s+t_r}{\cos {\mathrm{\θ}}_s+{\cos \mathrm{\θ}}_r}$

通过射线追踪求得，t为层间旅行时；

表示不同炮检距的成像深度；z通过在共成像点道集上拾取获得。

如图5所示是偏移速度分析的实现流程，其步骤包括：

步骤a1、初始地质模型，并进入步骤a2；

步骤a2、建立初始速度模型，并进入步骤a3；

步骤a3、进行基于目标层位底界面的PSDM，并进入步骤a4；

步骤a4、进行偏移速度分析，并进入步骤a5；

步骤a5、判断是否建立了完整的速度模型，如果为“是”则进入步骤2，如果为“否”则进入步骤a6；

步骤a6、进行基于目标层位底界面的PSDM，并进入步骤a7；

步骤a7、拾取目标层底界面，并进入步骤a8；

步骤a8、更新速度模型，将目标层移至下一层位，并进入步骤a3。

其中，步骤a5中所述完整的速度模型为：从浅到深各层位齐全的以速度、密度、层位埋深作为参数的地震地质模型。

这种方法收敛速度快，适用于各种速度模型和炮检距。

图6A所示为采用基于共CCP道集的叠后时间偏移方法对大庆一条2D3C地震数据进行偏移之后的PP波效果图；图6B所示为采用基于共CCP道集的叠后时间偏移方法对大庆一条2D3C地震数据进行偏移之后的PS波效果图；图7A所示为采用本发明方法对大庆的一条2D3C地震数据进行的偏移处理之后的PP波效果图；图7B所示为采用本发明方法对大庆的一条2D3C地震数据进行的偏移处理之后的PS波效果图。由图6A、图6B、图7A、图7B对比可以看出，本发明的叠前深度偏移方法与以前方法相比，具有更好的效果：

第一、近地表(0-500m)成像效果好。

第二、成像分辨率提高。图6A和图6B所示，采用基于共CCP道集的叠后时间偏移方法所得到的便宜之后的PP波和PS波的效果中，时间偏移剖面上只能看到是多条连续的同相轴。而采用本发明的叠前深度偏移方法后，同相轴的分辨率提高：如图7A和图7B所示的浅层(1000-2000m)同相轴更加清晰。从图7A和图7B中还可以看出同相轴并不是连续的，而是存在很多小的断裂，但在图6A和图6B中并没有该效果；图7A和图7B的深部3000米处的同相轴主频明显提高，并且可看到两条时间偏移剖面上看不到的断层。

第三、转换波成像剖面上3000米处的同相轴的连续性提高，构造更加清晰。

如图8所示为本发明叠前深度偏移方法的模拟数据模型。定义该模型参数为：1000m×1000m，四层；速度(m/s)Vp＝2500，3000，3500，4000，4500；Vs＝1500，1800，2100，2400，2700；密度(g/cm³)＝2.250，2.250，2.250，2.250，2.250。采用弹性波有限差分正演多炮记录；观测***：道间距5m，5炮，每炮200道，初始炮点坐标x＝100m，炮间距dx＝200m；采样间隔0.5ms，采样点数：2401。图9A为图8的数据模型的地震数据Z分量记录图，图9B为图8的数据模型的地震数据X分量记录图。采用本发明的叠前深度偏移方法的偏移结果如图10(纵波偏移深度剖面)、图11(转换横波偏移深度剖面)所示；图12(纵波偏移深度剖面)和图13(转换横波偏移深度剖面)所示是采用随倾角开孔径的声波方程偏移结果；图14(纵波偏移深度剖面)和图15(转换横波偏移深度剖面)所示是采用固定孔径声波方程偏移的结果。成像点网格100*100，dx＝10m，dz＝10m；偏移孔径aper＝0.4，即孔径半径为成像点深度z乘0.4。

从图10-图15对比可以看出，本发明方法的偏移效果最好，构造归位准确、清晰；其次为按倾角开孔径的声波方程偏移方法；最差的是固定孔径声波方程偏移方法。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种叠前深度偏移方法，包括以下步骤：

步骤1、对地震数据进行偏移速度分析；

步骤2、通过射线追踪计算走时、弧长、出射角、入射角；

步骤3、进行偏移孔径计算；

步骤4、利用Kirchhoff积分法矢量偏移公式进行叠前深度偏移。

2.根据权利要求1所述的叠前深度偏移方法，其特征在于：步骤1中所述偏移速度分析采用基于深度聚焦的共成像点速度分析方法。

3.根据权利要求2所述的叠前深度偏移方法，其特征在于，所述步骤1包括如下步骤：

步骤a1、初始地质模型，并进入步骤a2；

步骤a2、建立初始速度模型，并进入步骤a3；

步骤a3、进行基于目标层位底界面的叠前深度偏移，并进入步骤a4；

步骤a4、进行偏移速度分析，并进入步骤a5；

步骤a5、判断是否建立了完整的速度模型，如果为“是”则进入步骤2，如果为“否”则进入步骤a6；

步骤a6、进行基于目标层位底界面的叠前深度偏移，并进入步骤a7；

步骤a7、拾取目标层底界面，并进入步骤a8；

步骤a8、更新速度模型，将目标层移至下一层位，并进入步骤a3。

4.根据权利要求1所述的叠前深度偏移方法，其特征在于：步骤2中，采用直射线追踪法或者采用基于费马原理的最小走时射线追踪法。

5.根据权利要求1所述的叠前深度偏移方法，其特征在于：步骤4中采用适用于非均匀各向异性介质的弹性波Kirchhoff积分法矢量偏移公式进行叠前深度偏移。

Images