CN105974479B - Gpu空间网格体的层析2d/3d各向异性深度域速度建模方法 - Google Patents

Abstract

本发明提出了一种GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，包括：对GPU上采用四阶精度龙格－库塔算法，在预设孔径的数据体内计算每个深度点到炮点和检波点所需要的走时时间，形成走时表；读取所述走时表和时间域道集，进行叠前深度偏移形成CRP道集，并叠加形成深度域剖面；在叠前时间偏移剖面上，解释层位用作沿层速度分析和建模，采用以下三种方式之一进行建模包括：利用扫描剖面实时解释速度；三维沿层自动速度分析；三维射线追踪网格层析建模。本发明可以解决深度域高精度成像问题。

Description

GPU空间网格体的层析2D/3D各向异性深度域速度建模方法

技术领域

本发明涉及地震勘探数据处理技术领域，特别涉及一种GPU空间网格体的层析2D/3D各向异性深度域速度建模方法。

背景技术

深度域速度建模不同于叠前时间偏移建立速度场。叠前时间偏移建立速度场，只需要考虑调整速度，使CRP道集(共反射点道集)能够拉平，由于时间域的特性，无论用什么速度，零偏移距的时间不变，因此比较容易调整。而深度域速度建模，既要考虑CRP道集是否拉平，又要考虑深度是否正确，而且上层的速度变化后，会影响到下覆CRP道集形态。由于偏移是在三维空间进行，侧向速度又会影响到速度分析值。因此三维速度建模是个非常复杂的过程。而且三维速度建模计算量非常大，是否能快速计算也是在实现过程中需要重点考虑的。

发明内容

本发明的目的旨在至少解决所述技术缺陷之一。

为此，本发明的目的在于提出一种GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，可以解决深度域高精度成像问题。

为了实现上述目的，本发明的实施例提供一种GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，包括如下步骤：

步骤S1，对GPU上采用四阶精度龙格－库塔算法，在预设孔径的数据体内计算每个深度点到炮点和检波点所需要的走时时间，形成走时表，其中，根据叠前时间偏移得到的均方根速度，转化成深度域层速度的初始层速度，进行第一次的走时时间计算；

根据前次速度建模的深度与层速度，进行第二次及其后的走时时间的计算；

步骤S2，读取所述走时表和时间域道集，进行叠前深度偏移形成共反射点道集CRP道集，并叠加形成深度域剖面；

步骤S3，在叠前时间偏移剖面上，解释层位用作沿层速度分析和建模，其中，采用以下三种方式之一进行建模包括：

(1)利用扫描剖面实时解释速度，其中，所述利用扫描剖面实时解释速度，包括如下步骤：利用扫描剖面交互速度分析：将二进制速度场按一定的间隔转换成ASCII码深度域层速度曲线，对每个速度点进行垂直方向从上到下的交互速度解释，所述交互速度解释，包括如下步骤：

形成扫描叠加段剖面，在所述扫描叠加段剖面上对该点的上覆地层的深度与层速度进行调高或调低，并在调整后对整个扫描剖面重新更新；

将CRP道集转到时间域，将层速度转成均方根速度，再用均方根速度对剖面扫描，形成速度由小到大的扫描段，将上述扫描段转到深度域，以保证在速度扫描时深度不变；

速度解释是在GPU实时形成的交互分析剖面上进行，每解释一个速度点，GPU刷新扫描剖面，使上覆地层速度变化对下覆地层造成的影响实时体现出来并形成叠加剖面，根据地层产状判断解释的速度是否合理；

解释完成后，将ASCII码速度曲线转换成二进制速度场；

速度分析的方法分三步：第一步是按一定的网格间距将速度场转换成速度层速度曲线；第二步是用交互的方法改变层速度曲线的值；第三步是利用GPU快速刷新全部扫描剖面，观察层速度的改变是否使下覆地层在叠加段上得到最强的叠加能量，如果叠加能量最强，说明选择的层速度正确，得到对应的新的层速度；

(2)三维沿层自动速度分析；

(3)三维射线追踪网格层析建模。

步骤S4，根据步骤S3进行初步的深度域速度分析和建模结果，对需要加载的井深数据进行井深校正。

进一步，在所述步骤S3中，所述三维沿层自动速度分析，包括如下步骤：

三维沿层自动拾取，从上到下逐层计算速度模型，包括：每层之间按照最大扫描能量自动计算每点的最佳层速度，并且自动拾取。

进一步，在所述步骤S3中，所述三维射线追踪网格层析建模，包括如下步骤：

采用三维射线追踪，对每个需要计算的网格做由网格节点向地表的射线追踪，根据每个网格节点的射线追踪的时差建立矩阵，求解该矩阵得到新的速度场。

进一步，在所述步骤S4中，采用倾斜各向异性TTI对需要加载的井深数据进行井深校正。

根据本发明实施例的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，可以解决深度域高精度成像问题。

随着计算机硬件的快速发展，GPU(图形处理器)并行计算因其高效的计算能力，越来越受到人们的重视。GPU最早用来做图形显示，由于一个处理器包含上千个计算核心，善于做大规模的并行计算，被逐步应用到各个领域中。本发明采用GPU并行计算，CUDA语言编程，大大提高了计算效率。

对陆地石油勘探来说，地震信号采集的地表往往不是一个平面，可能会有很大的起伏，深度建模还需要考虑到起伏地表情况。一般来说，速度场要在基准面上建立，而输入的时间域的地震数据是在浮动面上。浮动面要用比较小的平滑参数。既能解决低降速带静校正问题，又对偏移成像精度没有太大影响。浮动面的值保存在数据道头里，计算走时要以这个面的值为起点。本发明的采用的三种建模方法可以满足上述要求，即使三维速度建模计算量非常大，也能实现快速计算。

本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1为根据本发明实施例的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法的流程图；

图2为根据本发明实施例的扫描段交互解释速度的示意图；

图3为在叠前时间偏移剖面上画出的层位，计算三维沿层扫描的示意图；

图4为三维沿层扫描，自动拾取和异常点编辑后的示意图；

图5为TTI之前的深度剖面的示意图；

图6为应用TTI之后的深度剖面的示意图；

图7为速度体V的示意图；

图8为参数体δ的示意图；

图9为参数体θ的示意图。

具体实施方式

下面详细描述本发明的实施例，实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

下面参考图1至图9对本发明实施例的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法进行说明。

如图1所示，本发明实施例的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，包括如下步骤：

步骤S1，对GPU上采用四阶精度龙格－库塔算法，在预设孔径的数据体内计算每个深度点到炮点和检波点所需要的走时时间，形成走时表。

具体点，在GPU上采用四阶精度龙格－库塔算法，在一定孔径的数据体内计算每个深度点到炮点和检波点所需要的时间，并储存到磁盘上。

第一次计算走时的时候，采用的速度可以是根据叠前时间偏移得到的均方根速度，转化成深度域层速度的初始层速。第二次及以后的速度，可采用前次速度建模的深度与层速度。

步骤S2，读取走时表和时间域道集，进行叠前深度偏移PSDM形成共反射点道集CRP道集，并叠加形成深度域剖面。

步骤S3，在叠前时间偏移剖面上，解释层位用作沿层速度分析和建模，其中，采用以下三种方式之一进行建模。

需要说明的是，下面三种方法都是采用GPU作为基础计算硬件，编程是采用英伟达公司提供的CUDA编程语言，采用MPI、QT、CUDA混合编程，其中GPU并行计算用CUDA，交互解释用QT，多节点并行运算用MPI。

1)利用扫描剖面实时解释速度

具体地，对每一个速度点将深度域层速度和深度域CRP道集都转到时间域，并按一定的速度百分比扫描并叠加，再转回到深度域，形成图2所示的百分比扫描剖面。如图2所示，将层速度转成均方根速度，将深度道集转成时间道集，用不同的速度百分比做速度扫描，扫描完成后转回到深度域。

利用扫描剖面交互速度分析：将二进制速度场，按一定的间隔转换成ASCII码深度域层速度曲线，对每个速度点进行垂直方向从上到下的交互速度解释。

为了在速度扫描时深度不变，需要将CRP道集转到时间域，将层速度转成均方根速度，再用均方根速度对剖面扫描，形成速度由小到大的扫描段，再把这些扫描段转到深度域。图2是速度为百分之百时的剖面。如果能量集中在中部，说明道集表示道集已拉平；如果最大能量在右边，说明道集向上弯曲，速度偏小，需要调大上覆地层的层速度；如果最大能量在左边，说明道集向下弯曲，速度偏大，要调小上覆地层的层速度。每当上覆层速度曲线做了调整，GPU重新做以上的各种转换，刷新整个扫描版面。通过反复调层速度曲线，使最大能量集中在中部，就完成了该点的交互深度域层速度分析。

解释的方法是首先形成扫描叠加段剖面，在扫描叠加段剖面上对该点上覆地层的深度与层速度进行调高或调低。由于改变的是层速度，会对下覆地层的剩余速度产生影响，需要对整个扫描剖面重新更新。解释是在GPU实时形成的交互分析剖面上进行，每解释一个速度点，GPU刷新扫描剖面；这种做法，可以使上覆地层速度变化对下覆地层造成的影响实时体现出来，并且形成叠加剖面，根据地层产状，可用来判断解释的速度是否合理；由于实时交互解释计算层速度是垂直方向，重新偏移后速度还会出现不精确之处，需要多轮速度分析以取得精确的速度场。解释完成后，再将ASCII码速度曲线转换成二进制速度场。

本建模方法是在一个速度点上，利用扫描剖面实时解释速度。由于层速度的变化会引起深度变化和下覆地层的剩余速度变化，不利于靠扫描剖面解释剩余速度。本方法采用把深度域层速度转成时间域均方根速度，深度域道集转到时间域，在时间域扫描，扫描剖面就是初始速度按一定百分比形成不同的速度剖面，再将这些剖面按一定的间隔形成叠加段。

由于对于不同的扫描速度，T0时刻不变，将这些扫描段转到深度域，就可以在深度域观察剩余速度的变换趋势，并且每个反射段的深度不变，方便速度解释。

速度分析的方法分三步：第一步是按一定的网格间距将速度场转换成速度层速度曲线；第二步是用交互的方法改变层速度曲线的值；第三步是利用GPU快速刷新全部扫描剖面，观察层速度的改变是否使下覆地层在叠加段上得到最强的叠加能量，如果叠加能量最强，说明选择的层速度正确，得到对应的新的层速度。

这种层速度解释方法，不是直接拾取剩余速度的能量团，而是改变其上部的层速度，GPU即时刷新扫描版面，观察其剩余速度的变化，是最大能量剖面处在100％(扫描段中间)。以往速度分析是在速度谱上拾取最大能量团的方法。这种方法对于信噪比很低的资料有很大的弊端，速度谱上能量团很散无法正确拾取速度。采用该方法，对于信噪比很低的资料，也可以得到正确速度。

深度域拾取速度和时间域不同，每改变一个点，下覆地层的道集曲率都会发生变化，需要重新扫描。如果用传统的CPU(中央处理器)，每扫描一次需要花费数分钟，耗时太多无法达到实时交互。而使用GPU加速后，可实现实时刷新。该方法由于只考虑了垂向速度变化，如果地层倾斜或者横向速度变化大，用调整后的速度重新偏移后，会出现新的误差，但是比之前的误差小。重复解释速度2－3次可逐步消除误差。

从图2所示的剖面上可以观察到CRP道集的反射层是否校平。如果最大叠加段能量在中间(百分之百速度)，说明CRP道集的反射层是平的，如果最大叠加段能量不在中间，说明反射层不平，改变该速度点上方的层速度，采用GPU快速计算更新全部扫描剖面，并观察最大叠加段能量是否在中间。通过交互反复调整层速度曲线，使最大能量集中在中段，可得到一个新的速度场。由于该速度调整只考虑了速度垂向变化，没有考虑反射层的横向移动，因此还是一个好于前一个速度的近似速度场。要想提高精度需要重新偏移，再次调整，逐步逼近正确的速度，一般需要反复3－5次，可得到精确速度场。该方法的特点是简单、实用，不用画速度模型，采用GPU快速更新扫面剖面，使得实时交互解释得以实现。由于采用了扫描叠加段，可观察到低信噪比同向轴，因此适合低信噪比地区的深度建模。解释时，还可用叠加剖面监控地层深度，通过扫描段和叠加剖面的比较，可确定哪些能量段是可靠的。

2)三维沿层自动速度分析

如图3所示，三维沿层自动拾取，手动编辑从上到下，逐层计算速度模型。由于三维数据量大，每个速度点拾取耗费大量时间。可先在时间剖面上画一些层，每个层之间可按最大扫描能量，自动计算每点的最佳层速度，并且自动拾取，对于异常点可进行手动编辑去掉。这种方法效率高，适合大规模三维解释。每计算一层能量时，要考虑到上覆已修正速度的影响。由于是靠能量团解释剩余速度，遇到很低信噪比时能量团不集中，可以与上述利用扫描剖面实时解释速度结合使用。

三维沿层自动速度分析的建模方法是三维沿层自动速度分析。如图4所示，第一步，先在叠前时间偏移剖面上画一些层，画这些层的原则是尽可能沿同一套地层，并且层和层之间要包含几个反射层，以便计算剩余速度能量变化。第二步是由上向下，一层层解释计算，计算下一层时，会用到上覆地层的计算结果。由于速度不断变化，地层深度也会随着变化。由于解释的层位是在时间域，用新的速度转换到深度域，因此，解释的上下两层的地层始终不变。计算可沿层自动拾取最大能量，手动剔除异常值。对于比较大的三维工区，速度点比较多，逐个解释耗时比较长。采用沿层自动速度分析，既能快速拾取沿层的最大能量，又能保证质量，特别适合大区块三维数据。整个过程为第一，解释层位；第二，沿层扫描计算并自动拾取最大能量；第三，手动剔除一些异常值；第四，利用这一层的解释结果，计算新的层速度，再回到第二步，计算下一层，直到完成全部解释的层。

在叠前时间偏移剖面上画出若干套计算用的时窗，根据时窗进行计算；时窗的画法是大致沿地层的走向，可穿过断层，时窗的垂向高度一般要包含2－3套轴；浅层轴多，层位可画薄一些，深层轴少则画厚一些。在计算每一层时，采用的是按一定的层速度百分比层速度扫描，并统计每个变化的百分比速度的在该层的累加能量。

每计算完一层，会根据这层的最大能量进行自动拾取，将自动拾取的点和统计的能量在屏幕上显示出，用交互的方法对一些异常值进行剔除，再计算该层的层速度；由上到下完成全部层的计算。由于该方法只考虑速度垂向变化，没有考虑偏移后地层的横向移动，因此每调完一次需要重新偏速度线，再调速度，反复几次即完成速度建模；

3)三维射线追踪网格层析建模。

拾取多层剩余速度进行三维层析计算，可同时考虑垂向和横向速度变化。采用前两种方法得到速度场，在构造不太复杂的情况下，可取得比较好的效果。对于复杂构造，采取三维层析算法更加精确。该方法采用三维射线追踪，对于每个需要计算的网格做由网格节点向地表的射线追踪，然后根据每个网格节点的射线追踪的时差，建立一个矩阵，求解这个矩阵，得到一个新的速度场。由于每个节点要做射线追踪，计算量非常大，采用GPU可大大提高速度。求解矩阵方程时，也可采用迭代算法，使求解更加稳定。

三维射线追踪网格层析建模，不同于三维沿层自动速度分析的是，不是解释一层，计算一层，而是解释多层剩余速度，从浅到深全部，或者部分解释出来，在计算时考虑浅层速度对深层和侧面速度的影响，影响范围由射线路径确定。对于每个速度网格节点，由这一点向地表的射线追踪，利用射线走时和剩余速度的差异，计算出新的速度。计算深层速度时要考虑到上覆地层剩余速度计算对其产生的影响。对于复杂地区，横向速度变化大，该方向利用剩余速度计算时，不仅考虑垂向速度变化，还考虑横向速度变化，而前面两种建模方法只考虑垂向变化，横向变化要通过重新偏移产生，因此要反复分析速度来取得正确的速度。这种三维层析方法，可以同时考虑垂向和横向的速度变化，使建模更加合理。但是如果剩余速度变化过大，也会影响这种方法的计算精度。一般在前两种建模方法的基础上，使用该方法完成高精度建模。

具体地，利用剩余速度层析计算：首先形成剩余速度扫描剖面，在该剖面上解释CRP道集剩余速度，进行三维网格层析速度计算。对于每个需要计算的网格做由网格节点向地表的射线追踪，然后给这个网格一个速度的扰动，会引起这个网格节点处的道集曲率变化；如果这点的曲率变化和这点剩余速度造成的曲率变化相同，那么这点的速度变化就可以确定。该方法的网格节点计算由上到下，计算下覆地层的节点速度值时，射线追踪已用到了上覆地表新计算的速度值，可以扣除上覆地层速度变化造成的对下覆地层剩余速度的影响，所以可以精确计算。该方法可一次拾取多层剩余速度，统一计算，不用一层层计算，效率高。

步骤S4，根据步骤S3进行初步的深度域速度分析和建模结果，对需要加载的井深数据进行井深校正。

在初步完成建模并做叠前深度偏移，形成深度域剖面后，要加载井资料核查井深，如果井深与剖面不一致可通过调节各向异性参数来解决。对于构造变化比较平缓，可用VTI解决，如果地层变化大、倾角大，采用倾斜各向异性TTI对需要加载的井深数据进行井深校正。对一套地层来说，各向异性参数相对稳定，在过井点求出的各向异性参数，可沿层外推到其区域。通过与井深的误差，可以直接算出δ、ε参数。使得地层和井深相吻合，同时保持反射同向轴平。θ、φ(地层倾角和方位角)的求取在偏移剖面上直接通过相关法求出。可先求出地层主测线和相交线的倾角DIPX和DIPY，再算出θ、φ。图7至图9分别给出了速度体V、参数体δ和参数体θ的示意图。

图5和图6分别给出了TTI前后的深度剖面的示意图，通过TTI，可以消除剖面与井深误差。

如果CRP道集从浅层到深层都已拉平，但叠加剖面的深度和井深不一致，说明存在各向异性问题。各向异性分为VTI和TTI。VTI只考虑δ、ε两个参数，而TTI需要用到δ、ε、θ、φ四个参数。其中δ主要影响地层的深度，ε参数主要影响CRP道集是否拉平，θ是地层的倾角，φ是地层的方位角。θ可以从深度叠加数据体直接计算出来。通过不断调整δ和ε参数，可以得到正确的井深位置，同时也拉平道集。由于δ、ε对一套地层来说变化不大，因此得到一组δ、ε参数后，可以沿层外推到构造的其他部分。

根据本发明实施例的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，可以解决深度域高精度成像问题。

随着计算机硬件的快速发展，GPU(图形处理器)并行计算因其高效的计算能力，越来越受到人们的重视。GPU最早用来做图形显示，由于一个处理器包含上千个计算核心，善于做大规模的并行计算，被逐步应用到各个领域中。本发明采用GPU并行计算，CUDA语言编程，大大提高了计算效率。

对陆地石油勘探来说，地震信号采集的地表往往不是一个平面，可能会有很大的起伏，深度建模还需要考虑到起伏地表情况。一般来说，速度场要在基准面上建立，而输入的时间域的地震数据是在浮动面上。浮动面要用比较小的平滑参数。既能解决低降速带静校正问题，又对偏移成像精度没有太大影响。浮动面的值保存在数据道头里，计算走时要以这个面的值为起点。本发明的采用的三种建模方法可以满足上述要求，即使三维速度建模计算量非常大，也能实现快速计算。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。本发明的范围由所附权利要求及其等同限定。

Claims (4)

1.一种GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，其特征在于，包括如下步骤：

步骤S1，对GPU上采用四阶精度龙格－库塔算法，在预设孔径的数据体内计算每个深度点到炮点和检波点所需要的走时时间，形成走时表，其中，根据叠前时间偏移得到的均方根速度，转化成深度域层速度的初始层速度，进行第一次的走时时间计算；

根据前次速度建模的深度与层速度，进行第二次及其后的走时时间的计算；

步骤S2，读取所述走时表和时间域道集，进行叠前深度偏移形成共反射点道集CRP道集，并叠加形成深度域剖面；

步骤S3，在叠前时间偏移剖面上，解释层位用作沿层速度分析和建模，其中，采用以下三种方式之一进行建模包括：

(1)利用扫描剖面实时解释速度，其中，所述利用扫描剖面实时解释速度，包括如下步骤：利用扫描剖面交互速度分析：将二进制速度场按一定的间隔转换成ASCII码深度域层速度曲线，对每个速度点进行垂直方向从上到下的交互速度解释，所述交互速度解释，包括如下步骤：

形成扫描叠加段剖面，在所述扫描叠加段剖面上对该点的上覆地层的深度与层速度进行调高或调低，并在调整后对整个扫描剖面重新更新；

将CRP道集转到时间域，将层速度转成均方根速度，再用均方根速度对剖面扫描，形成速度由小到大的扫描段，将上述扫描段转到深度域，以保证在速度扫描时深度不变；

速度解释是在GPU实时形成的交互分析剖面上进行，每解释一个速度点，GPU刷新扫描剖面，使上覆地层速度变化对下覆地层造成的影响实时体现出来并形成叠加剖面，根据地层产状判断解释的速度是否合理；

解释完成后，将ASCII码速度曲线转换成二进制速度场；

速度分析的方法分三步：第一步是按一定的网格间距将速度场转换成速度层速度曲线；第二步是用交互的方法改变层速度曲线的值；第三步是利用GPU快速刷新全部扫描剖面，观察层速度的改变是否使下覆地层在叠加段上得到最强的叠加能量，如果叠加能量最强，说明选择的层速度正确，得到对应的新的层速度；

(2)三维沿层自动速度分析；

(3)三维射线追踪网格层析建模。

步骤S4，根据步骤S3进行初步的深度域速度分析和建模结果，对需要加载的井深数据进行井深校正。

2.如权利要求1所述的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，其特征在于，在所述步骤S3中，所述三维沿层自动速度分析，包括如下步骤：

三维沿层自动拾取，从上到下逐层计算速度模型，包括：每层之间按照最大扫描能量自动计算每点的最佳层速度，并且自动拾取。

3.如权利要求1所述的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，其特征在于，在所述步骤S3中，所述三维射线追踪网格层析建模，包括如下步骤：

采用三维射线追踪，对每个需要计算的网格做由网格节点向地表的射线追踪，根据每个网格节点的射线追踪的时差建立矩阵，求解该矩阵得到新的速度场。

4.如权利要求1所述的GPU空间网格体的层析2D/3D各向异性深度域速度建模方法，其特征在于，在所述步骤S4中，采用倾斜各向异性TTI对需要加载的井深数据进行井深校正。