CN104422953B - 一种提高地震叠前时间偏移计算效率的方法 - Google Patents

Abstract

本发明提供了一种提高地震叠前时间偏移计算效率的方法，属于石油勘探中的地震资料处理领域。所述方法首先计算在时间轴上相隔步长为n的成像点的旅行时间，然后采用白适应线性插值计算方法得到其余成像点的旅行时间。本发明在未降低成像效果的前提下，提高了计算效率，节约了计算成本。

Description

一种提高地震叠前时间偏移计算效率的方法

技术领域

本发明属于石油勘探中的地震资料处理领域，具体涉及一种提高地震叠前时间偏移计算效率的方法。

背景技术

在石油地震勘探等领域，叠前时间偏移(Pre-stack time migration，PSTM)是构造成像的有效方法之一，能有效适应地层速度变化，适用于大倾角和低信噪比的偏移成像方法，被工业界广泛的应用。

叠前时间偏移是常规地震资料处理的一个非常重要处理步骤，根据数据量的不同和计算集群的规模，叠前时间偏移计算时间从几天到几周，数据量大的甚至达到数月。叠前时间偏移的计算量和需要处理的数据量是巨大的，PSTM每输出一个地震道，就是一次海量计算。以现有常规三维采集的情况下，偏移一个地震道的输出至少需要数百万道甚至更多(偏移孔径决定)的输入道数据，每个点需要进行旅行时的运算和振幅补偿，同时数据道要进行频率滤波处理。

目前，工业界使用大规模的服务器或微机集群来进行叠前偏移处理，其基本原理是将数据分配到各个CPU核上，然后由各个CPU核单独进行计算，最后将结果汇总输出；除了利用CPU进行并行计算外，GPU也是一个很好的选择，这主要是GPU拥有非常可观的浮点运算能力。随着GPU硬件和技术的不断发展，利用GPU集群进行叠前偏移计算，使其计算效率得到了大幅度的提高。但是随着处理数据的不断增加，计算成本仍然非常昂贵，因此，从算法方面提高偏移的计算效率是非常重要。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种提高地震叠前时间偏移计算效率的方法，通过分析叠前时间偏移中旅行时计算的特点，采用自适应的旅行时插值计算方法，减少计算量，提高叠前时间偏移计算效率，为地震资料处理节约计算成本。

本发明是通过以下技术方案实现的：

一种提高地震叠前时间偏移计算效率的方法，首先计算在时间轴上相隔步长为n的成像点的旅行时间，然后采用自适应线性插值计算方法得到其余成像点的旅行时间。

所述方法包括以下步骤：

步骤1：设xt1成像点为需要计算的起始时间点，计算步长为n的三个成像点的旅行时间，三个成像点分别为xt1、xt1+n和xt1+2n，计算得到的旅行时间分别为T1、T2和T3；

步骤2：通过插值计算得到xt1+n成像点的旅行时间T2c；

步骤3：通过插值计算得到的旅行时间T2c与实际计算得到的旅行时间的差值为T2c-T2，如果T2c-T2小于第一阀值eposl，则步长n变大为2倍，即n＝n*2了；如果T2c-T2大于第二阀值epos2，则步长n变小为一半，即n＝n/2；如果上述两个条件均不满足，则步长n保持不变；

步骤4：通过插值计算得到成像点xt1至xt1+2n之间的其他所有成像点的旅行时间；

步骤5：判断是否所有需要计算的点都计算完成，如果否，则设xt1＝xt1+2n，然后返回步骤1，如果是，则结束。

所述步骤1中n的初始值为4、8、16或32。

所述步骤2是这样实现的：采用线性插值方法，T3和T1之间的时间差分为2n等分，每一等分为(T3-T1)/2n，距离xt1为n等分的xt1+n点的旅行时间为：T2c＝T1+n*(T3-T1)/2n)＝(T3+T1)/2。

所述步骤3中的第一阀值eposl为采样时间dt的百分之十五，第二阀值 epos2为采样时间的百分之四十。

与现有技术相比，本发明的有益效果是：本发明在未降低成像效果的前提下，提高了计算效率，节约了计算成本。

附图说明

图1是叠前时间偏移计算的步骤框图。

图2是叠前时间偏移旅行时计算与成像代码。

图3是速度曲线。

图4是旅行时间曲线(成像点和炮检中点距离XX为1500米)。

图5是自适应插值计算出的旅行时与实际旅行时的差。

图6是自适应插值计算旅行时。

图7-1是常规计算的结果。

图7-2是采用本发明方法计算的结果。

图8是本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

克希霍夫叠前时间偏移公式为：

其中：V(M)为M点的像，为振幅加权因子，为地震数据，A为成像点的积分范围，即偏移孔径，为绕射双曲线。克希霍夫叠前时间偏移过程实际上是根据观测系统，计算成像点至每组观测点对(即炮点和检波点)的时间，在对应在数据中把该时间点的振幅值按一定加权原则叠加至成像点，即获得该成像点的像，所有像点都计算完成，则偏移完成，实现步骤如下 (见图1，以2D为例)：

步骤1：加载观测系统，并输入成像范围(NX，NT)(成像范围是由用户自己规定的)。

步骤2：加载成像范围内的速度模型：

地震数据的处理中，专门有一项是速度建模，作为整个处理流程的一个非常重要的环节，为包括偏移等其他处理方法提供速度模型使用，地震处理中很多步骤都需要用到速度模型；速度模型加载过程就是根据提供的速度模型的格式，读取速度值，放入内存中备用即可。

步骤3：观测数据循环(观测系统是由很多接收点(即检波点)和炮点组成，这里的观测数据循环是指每个炮检对(炮点和检波点)所对应的数据循环)，加载一道数据，并获得该道数据的炮检坐标。

步骤4：成像点循环开始，计算每一个成像点至该道观测数据炮检坐标的旅行时间，在数据中取相应时间的振幅，乘以振幅加权因子赋给像点。

步骤5：所有成像点计算完成。

步骤6：所有数据计算完成。

步骤7：偏移结束。

从叠前时间偏移计算流程分析，可以看出每一次成像计算需要三重循环，若是3D偏移则是四重循环，而旅行时间的计算是位于循环的最里层，是最主要的计算量，因此，优化旅行时间的计算，能有效的降低计算成本。图2为叠前时间偏移中核心计算部分的C语言代码。

本发明就是对上述步骤4的计算进行改进。本发明的内容和原理：采用自适应线性插值计算方法对上述计算进行优化，即通过在时间轴上采取一定的步长计算旅行时间(即只计算小部分)，其余点的旅行时通过插值方式计算出来。

图3为常规叠前速度分析得到的速度模型，可以看出来速度模型通常呈线性增加的，因为地层越深，地震波传播速度越大，与实际情况相符合。图4为某一成像点所计算的旅行时间(炮检中点距离成像点1500米)曲线，旅行时间为光 滑的曲线，而且是线性增加的，特别是在成像时间大的地方，曲线越平缓，因此，可以采用一段段的直线来拟合该时间曲线，旅行时只需计算直线的两个端点(指步长为n的两个端点)，中间部分通过插值得到，达到节省计算量的目的。图5为自适应插值计算的旅行时间与实际计算结果的误差，其中dt为2ms，从图中可看出，最大误差仅为0.35ms。

如图8所示，本发明的具体实现步骤如下：

步骤1：设xt1点为需要计算起始的时间点，计算步长为n(即n个成像点)的三个点的旅行时间(见图6)，三个点分布为(xt1，xt1+n，xt1+2n)，计算的旅行时间为(T1，T2，T3)(采用通用的旅行时计算公式计算得到的，例如可参考“起伏地表叠前时间偏移技术研究”，朱海波，张兵，2007，Vol30，No5，勘探地球物理进展，P368，见公式(4))，通常n的初始值给为16(这个值可以是4，8，16，32等，值越大，计算效率越高，但误差可能会增大，当速度垂直变化不大的情况下，值可以放大，反之，应该给小；16是一个实际测试中较为通用的值)。

步骤2：通过插值计算xt1+n点的旅行时间设为T2c：这里采用的是线性插值方法，T3和T1之间的时间差分为2n等分，每一等分为(T3-T1)/2n，那么距离Xt1为n等分的xt1+n点的时间：T2c＝T1+n*(T3-T1)/2n)＝(T3+T1)/2。

步骤3：如果插值计算的旅行时与实际计算的旅行时差(T2c-T2)小于某一阀值epos1，步长n变大为2倍，即n＝n*2。如果插值计算的旅行时与实际计算的旅行时差(T2c-T2)大于某一阀值epos2，步长n变小为一半，即n＝n/2。上述两个条件均不满足，则步长n不变。通过大量的测试和试验，本发明采用的阀值选取标准：eposl为采样时间dt的百分之十五，epos2为采样时间的百分之四十(这个是发明人通过大量实际数据总结的一个经验值，选择标准：不降低效果的情况下，最大化提高计算效率。Eposl过小，会降低效率，epos2过大会降低效果)。

步骤4：插值出点xt1至xt1+2n之间的其他所有点的旅行时间。

步骤5：按新的步长n计算下一个点(此时起点是时间为T3的点xt1+2n，步长为新的n，下一个点是xt1+2n+n)，重复步骤1～4，直至所有需要计算的点计算完成。

为了测试改进的效果，采用国际上通用marmousi模型数据进行测试，测试环境为Linux Version2.6.9-55.0.12.ELsmp操作系统，处理器为8核Intel(R)Xeon(R)CPU X5355@2.66GHz，内存为16G，计算一道数据进行偏移，比较改进前后计算所需的时间，由于单道偏移计算时间很小，因此对其进行2000的重复计算后，统计总的时间，单道计算效率比较的计算结果对比如表1所示，从表1可以看出，对一道数据进行处理时，利用本发明方法对核心算法部分进行改进后，效率提高了51.9％：

表1

对整个模型数据(上面表1中只是对数据中的一道进行测试，marmousi模型数据有很多道，整个模型数据是指所有道的数据)进行计算，采用相同的计算环境，模型数据计算效率结果如表2所示，从表2可以看出，对所有道的数据进行处理时，利用本发明方法对核心算法部分进行改进后，效率提高了约30％。表1和表2中提高效率有差别，这主要是因为本次采用测试数据本省较小(72M)，同时时间采样点数较小为750个点，偏移过程中的其他计算量占据了一部分时间(如数据和速度模型的加载、数据的滤波和结果的输出等)，从而降低了加速时间在整个计算时间的比重。图7-1和图7-2为本发明计算结果与未做效率优化后的计算结果对比，在成像质量上一致，因此，本发明在未降低成像效果的前提下，提高了计算效率，较大的节约了计算成本。

	改进后耗时	改进前耗时	提高效率
				模型数据	366.490m	476.381m	30.0％

表2

分析克希霍夫叠前时间偏移实现过程，找出其计算热点，在2D情况下，旅行时间计算是需要进行三重循环，因此，基于旅行时间的光滑且单调递增的规律，开发出一种在时间轴上采用自适应步长的旅行时计算方法，只计算少量点的旅行时间，其他点旅行时通过插值获得。步长的确定是通过成像点的插值计算和直接计算结果之间的误差来确定的，当误差大于某阀值时，步长自动减半，当误差小于某一较小阀值时，步长自动倍增，这样，在保证了旅行时计算精度(成像品质不降低)的情况下，减少了计算量，节约了计算成本。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (4)

1.一种提高地震叠前时间偏移计算效率的方法，其特征在于：所述方法首先计算在时间轴上相隔步长为n的成像点的旅行时间，然后采用自适应线性插值计算方法得到其余成像点的旅行时间,所述方法包括以下步骤：

步骤1：设xt1成像点为需要计算的起始时间点，计算步长为n的三个成像点的旅行时间，三个成像点分别为xt1、xt1+n和xt1+2n,计算得到的旅行时间分别为T1、T2和T3；

步骤2：通过插值计算得到xt1+n成像点的旅行时间T2c；

步骤3：通过插值计算得到的旅行时间T2c与实际计算得到的旅行时间的差值为T2c-T2，如果T2c-T2小于第一阀值epos1，则步长n变大为2倍，即n＝n*2了；如果T2c-T2大于第二阀值epos2，则步长n变小为一半，即n＝n/2；如果上述两个条件均不满足，则步长n保持不变；

步骤4：通过插值计算得到成像点xt1至xt1+2n之间的其他所有成像点的旅行时间；

步骤5：判断是否所有需要计算的点都计算完成，如果否，则设xt1＝xt1+2n，然后返回步骤1，如果是，则结束。

2.根据权利要求1所述的提高地震叠前时间偏移计算效率的方法，其特征在于：所述步骤1中n的初始值为4、8、16或32。

3.根据权利要求2所述的提高地震叠前时间偏移计算效率的方法，其特征在于：所述步骤2是这样实现的：采用线性插值方法，T3和T1之间的时间差分为2n等分，每一等分为(T3-T1)/2n，距离xt1为n等分的xt1+n点的旅行时间为：T2c＝T1+n*(T3-T1)/2n)＝(T3+T1)/2。

4.根据权利要求3所述的提高地震叠前时间偏移计算效率的方法，其特征在于：所述步骤3中的第一阀值epos1为采样时间dt的百分之十五，第二阀值epos2为采样时间的百分之四十。