CN103439743B - 一种地震数据剖面图形绘制方法 - Google Patents

Abstract

本发明公开了一种地震数据剖面图形绘制方法，采用纵向插值的方法计算出毎道地震波形曲线图对应的设置为线条色的像素点的个数和位值，绘制出地震数据波形曲线，并在波形曲线的基础上，将地震数据的正值区域或者负值区域填充成填充色。相对于现有的方法，本发明的方法简单，且能保证图形绘制效果，适合GPU处理，提高了图形绘制效率；解决了现有的地震数据图形绘制方法对输入的地震数据需要进行特别的处理，不能直接实现绘制的问题。

Description

一种地震数据剖面图形绘制方法

技术领域

本发明涉及地震数据图形绘制领域,特别是一种地震数据剖面图形绘制方法。

背景技术

数据可视化技术已经成为科学研究中的一种流行的技术,这种技术能够提高工程的效率和准确性。数据可视化利用计算机图形学等技术,对数据进行处理,将其转换为可显示在屏幕上的图形或者图像,通过交互处理来挖掘数据相互间的规律。数据可视化技术已经被广泛用于多种领域,特别是在地质勘探工程上的应用,使得地质勘探减少了风险和成本,同时也为其提高了对地震数据处理的效率和质量。

地震数据可视化在国内外的地质解释系统中都有实现，是地质解释软件系统中的基础而又重要的模块。地震数据可视化分为二维和三维地震数据的可视化。二维地震数据图形绘制是二维地震数据可视化的基础，也是三维地震数据分剖面可视化的基础，因此二维地震数据图形绘制是地质解释软件的基础功能。随着数据可视化的不断完善，人们针对二维地震数据的图形绘制提出比较成熟的绘制方法，在“肖汉.地震数据的可视化研究技术研究.长沙：湖南大学，2007,23-25”中提出了基于OpenGL绘制引擎的二维地震剖面B+W图的绘制方法。在“黄丹.基于QT的地震数据可视化的研究及应用[D].成都：电子科技大学，2011,45-50”中提出了基于Qt绘制引擎的二维地震剖面的绘制方法。

随着测量技术的发展，地震数据规模越来越大，目前的地震数据图形绘制方法中地震数据处理的时间比较长，图形光栅化处理效率低，导致了基于Qt和OpenGL等通用绘制引擎的地震数据图形绘制难以支撑应用的需求。

随着高性能计算、图像处理和计算几何等技术的不断发展，尤其随着GPU硬件的发展和GPU可编程性的不断提高，利用GPU完成通用计算的研究渐渐活跃起来。为了提高地震数据可视化的效率，GPGPU技术也逐渐应用到地质勘探软件中。在“K.Xie,P.Wu,S.Yang.GPUand CPU coorperation parallel Visualisation for large seismicdata.Electronics Letters，2010”中提出了GPU和CPU协同并行实现海量地震数据可视化；在“Daniel Patel,Stefan Bruckner,Ivan Viola,Eduard M.Groller.Seismic VolumeVisualization for Horizon Extraction.IEEE Pacific Visualization Symposium，2010”中提出了基于GPU的层位提取的层位生长算法。但是目前还没有适合GPU的地震数据图形绘制方法，并且现有的地震数据图形绘制方法对输入的地震数据需要进行特别的处理，不能直接实现绘制。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种适合GPU处理的地震数据剖面图形绘制方法，提高图形绘制效率；解决现有的地震数据图形绘制方法对输入的地震数据需要进行特别的处理，不能直接实现绘制的问题。

为解决上述技术问题，本发明所采用的技术方案是：一种地震数据剖面图形绘制方法，具体包括如下步骤：

1).设每一道地震数据X方向的间距为traceDistance，Y方向的间距为sampleHeight，输入的地震数据为Date[sampleCount]，输入的地震数据的绝对值最大值为MaxDate，纵向插值的起始点的索引值为sIndex，纵向插值的终止点的索引值为eIndex，则待插值点插值后的地震数据值，即地震数据波形曲线经过的像素点对应的地震数据值InterpData[Index]为：

InterpData[Index]=Date[sIndex]+α*(Date[eIndex]-Date[sIndex])，

其中，α=(Index-sIndex)/(eIndex-sIndex)；Date[sIndex]和Date[eIndex]分别为线性插值起始地震数据值和终止地震数据值；Index为纵向插值起始点和终止点之间的点的索引值；

2).根据下式计算每个所述像素点在特定区域的位置：

其中，Location[Index]为像素点偏离中心位置的像素宽度，Location[Index]的正负分别表示该像素点位于中心位置的左侧或右侧；PixelDate[Index]是地震数据波形曲线经过的所有像素点对应的地震数据值；其中所述特定区域是指毎道地震数据绘制后所占的像素区域，该像素区域的宽度是traceDistance；

3).根据上述步骤1）和步骤2）的计算公式绘制地震数据波形曲线,即根据步骤1）和步骤2）像素点及各像素点的位置确定地震数据波形曲线。

4).将毎道地震数据纵向插值生成的像素点和地震数据对应的像素点标记为边界点，记为BoundaryPoint；判断在所述特定区域内的其他非边界点的像素点与边界点的相对位置关系，将满足填充条件的像素点设置为填充色，其他像素点设置为背景色，从而完成地震数据剖面图形绘制；其中所述插值生成的像素点即地震数据波形曲线经过的像素点。。

所述步骤4）中，对于地震数据的正值区域，填充条件是一道地震数据波形曲线所在的像素区域内的所有的像素点位于地震数据对应的像素点和插值生成的像素点的左边部分；对于地震数据的负值区域，填充条件是一道地震数据波形曲线所在的像素区域内的所有的像素点位于地震数据对应的像素点和插值生成的像素点的右边部分。

与现有技术相比，本发明所具有的有益效果为：本发明的方法不但能应用到CPU上，还适合GPU处理，从而大大提高了地震数据图形绘制效率；传统的地震数据图形绘制方法是调用绘制引擎的API函数实现绘制的，是串行实现的，而本发明的方法在保证绘制效果的前提下，不依赖具体的绘制引擎，对输入的地震数据不需要进行特别的处理，能直接实现绘制，从而简化了绘制过程，使绘制效率得到进一步提高，具有更强的实用性；本发明大大减少了地震数据图形绘制中地震数据处理的时间，可应用于地质解释软件中，提高地震数据的可视化效率。

附图说明

图1为本发明一道地震数据波形曲线绘制示意图；

图2为本发明地震波形曲线图填充示意图；

图3为cuda openGL互操作实现流程图；

图4为使用传统方法绘制的地震数据剖面图形；

图5为使用本发明方法绘制的地震数据剖面图形。

具体实施方式

下面结合附图对本发明的实施例做进一步的说明。

基于通用绘制引擎的地震剖面绘制算法分析可知，地震波形曲线是根据地震数据值、X方向的显示间距和Y方向的间距分布的，因此毎道的地震波形曲线图对应的像素点集固定在特定的区域内。因为波形曲线是由很多的小线段构成的，根据图形显示光栅化的原理可知任何一条线段最终都对应一系列连续的像素点集。因此，相对传统的绘制方法中首先计算地震数据对应的坐标位置，然后再首尾相连，最后由绘制引擎实现线的光栅化处理的流程，本发明用B+W图并行算法直接实现线段的光栅化处理，采用纵向插值的方法计算出毎道地震波形曲线图对应的设置为线条色的像素点的个数和位值。

以一道地震数据的波形曲线图绘制为例，假设X方向的间距为traceDistance，Y方向的间距为sampleHeight，输入的地震数据为Date[sampleCount]，输入地震数据的绝对值最大值为MaxDate，纵向插值的起始点的索引值为sIndex，终止点的索引值为eIndex，则纵向插值计算像素点对应的地震数据值公式如式（2-1）所示。每个像素点在特定区域的位置计算公式如式（2-2）所示。一道地震数据的波形曲线绘制示意图如图1所示。

α=(Index-sIndex)/(eIndex-sIndex)

InterpData[Index]=Date[sIndex]+α*(Date[eIndex]-Date[sIndex]) （2-1）

其中：InterpData[Index]是待插值点插值后的值，Date[sIndex]和Date[eIndex]分别线性插值起始值和终止值。

其中：PixelDate[Index]是地震数据的波形曲线经过的所有像素点对应的地震数据值，traceDistance是X方向的间距，Location[Index]像素点偏离中心位置的像素宽度，该值的正负表示该像素点位于中心位置的左侧还是右侧。

区域填充是在波形曲线的基础上，将地震数据的正值区域或负值区域填充成填充色，因此区域填充光栅化处理就是对毎道所占的固定的像素区域的满足要求的像素点设置成填充色。正值区域的填充和负值区域的填充实现的原理是完全一样的，唯一的区别是填充判断的条件不同。正值区域填充时的填充条件是判断的像素点位于原始地震数据对应的像素点的左边部分，而负值区域填充时的填充条件是判断的像素点位于原始地震数据对应的像素点的右边部分。

以一道地震波形曲线的正值区域填充为例，具体实现方法为：将毎道地震数据纵向插值生成的像素点和原始地震数据对应的像素点标记为边界点，记为BoundaryPoint；判断在固定区域内的其他非边界点的像素点(JudgePoint)与边界点的相对位置关系，在参考点的右边部分的像素点设置为填充色，其他像素点设置为背景色。填充示意图如图2所示。

以下介绍本发明方法在GPU上的实现过程：

在CUDA框架下，一个支持CUDA的GPU作为CPU的协处理器，适用于可以分解为SIMD并行模式的算法。它有自己的静态存储器、片上共享存储器以及由一系列流处理单元(streaming processor)整合而成的流式多处理器(streamming multiprocessors，MP)。静态存储器可以分为只读纹理存储器(texture memory)、读写全局存储器(global memory)、以及只读常量存储器(Constant Memory)，存储于其中的数据可以被所有的流式多处理器访问处理，但速度较慢。而每一个流式多处理器中的读写共享存储器(shared memory)的访问速度较快，但它只能被绑定其上的固定流式多处理器中的处理单元访问。就编程模式而言，GPU执行的所有任务都由一系列用户为GPU分配的线程模块(Block)中的线程(Threads)完成。各个线程以阵列的形式组成了线程模块，而所有的线程模块又以阵列的形式组成—个最终并行计算的网格(Grid)。一旦通过一系列CUDA API函数指定网格和线程模块的维数大小，CUDA API中的内嵌变量就会据此来为这张网格中的每一个线程指定—个索引号。用户通过索引号来控制线程对存储于GPU上的数据读写或运算。而同时，GPU内部根线程模块的分配情况以及存储器的使用情况，使得尽可能多的流式多处理器同时工作，顺序完成程序上分配好的并行线程。MP资源并行的效率受线程模块大小以及位于MP上的共享存储器的制约。当一个Grid中所需并行处理的总线程数量被给定时，每一个线程模块的大小都必须以最大化可用运算资源(在同一个时钟周期内，让尽可能多的流处理单元参与运算)为目标，合理指定线程模块以及共享存储器的大小，同时尽量避免内存访问冲突(既在并行计算时，不同的线程同时访问某内存中的同—个地址而造成的冲突)。

因为本发明提出的并行实现地震剖面图形绘制方法中计算每一个像素值的颜色值具有很好的独立性，其求解过程以及所利用的像素资源之间均不存在相关性，GPU要求的处理的数据之间没有相干性。因为计算像素点的颜色值的方法一样，又满足GPU要求的数据处理的方法一致的设计准则,所以可以利用GPU的并行计算特性进行加速处理。在英特尔Pentium(奔腾)双核E6600@3.06GHz，4G内存，Nvidia GeForce GTX560显卡，1G显存以及在CUDA Tookit v4.2环境下采用CUDA技术实现该算法。GPU计算的关键流程主要分一下几个步骤：

1）：检测GPU性能参数。对于CUDA开发来说，GPU的硬件指标非常重要。针对硬件特点的优化能大幅改善程序的性能。检测系统硬件对CUDA程序来说是必不可少的，以便在程序选定GPU并对其初始化，只有这样，才能开始以后的计算。

2）：将源图像数据复制到CUDA数组，用于GPU并行快速计算。CUDA平台有别于经典GPGPU的重要之处在于它的全局存储器和共享存储器，用户像使用主机内存一样的方便地使用它们。纹理存储器是图形处理器上最重要的存储器，对符合纹理存储器特性的算法，用纹理存储器实现时能有比用全局存储器跟更高的效率。

纹理缓存特别对邻近的数据表现出较大的数据带宽。如果对二维的邻近数据访问得频繁，使用纹理存储器能获得性能的提升。图像上的某一像素块就是这样的数据。通过绑定纹理，一块纹理可以是线性存储器也可以是CUDA数组。因为CUDA数组时NVIDIA针对纹理拾取而特别优化的，所以从CUDA数组中拾取纹理相比从线性存储器中拾取纹理有更多的好处。诸如CUDA数组可以是一维、二维或者三维的，其中的元素可以是标量，或者二维、四维的向量。但线性存储器只能是一维的等等。这就是我们要选择CUDA数组的原因。

3）：指定CUDA纹理属性，绑定纹理到CUDA数组。纹理在纹理存储器上的位置有纹理引用确定，而纹理元在纹理上的位置由纹理坐标确定。由于在实际应用中纹理一般都与CUDA数组绑定。纹理引用是纹理绑定数组的前提，它必须是作用域为文件的静态全局变量，且不允许作为传递给其他函数的参数。纹理引用按如下方式被声明为texture类型：Texture<Type,Dim,ReadMode>texRef;

本发明中的纹理声明为texture<uchar4,2,cudaReadModeNormalizedFloat>rgbaTex；指定了CUDA纹理的属性后，接下来需要分配CUDA数组，然后将数据从主机传输至CUDA数组，再将纹理存储器与CUDA数组绑定。这时，就可以调用使用纹理存储器的内核了。

4）：创建像素缓冲区对象与显示纹理，CUDA和OpenGL互操作。在创建和绑定CUDA纹理后，因为CUDA数组时只读的，我们还需要分配另一个全局存储器来存放结果，本发明用的是unsigned int*dptr。

为了方便和快速显示，本发明利用OpenGL显示结果，利用了cuda和OpenGL互操作实现显示。CUDA和OpenGL互操作的基本方式是使用CUDA生成数据，然后使用OpenGL在屏幕上绘制出数据所表示的图形。两者的结合可以通过两种方式来实现：

i.使用OpenGL的PBO(像素缓冲区对象)。在该方式下，CUDA直接生成像素数据，OpenGL显示这些像素；

ii.使用OpenGL的VBO(顶点缓冲区对象)。在该方式下，CUDA生成顶点网格数据，OpenGL可以根据需要绘制出平滑的表面图或线框图或一系列顶点。

这两种方式的核心都是利用cudaGLMapBufferObject函数将OpenGL的缓冲区映射到CUDA的内存空间上，这样，程序员就可以充分利用CUDA的优点写出性能高的程序在该内存空间上生成数据，这些数据不需要传送，OpenGL可以直接使用。如果不使用CUDA，这些数据需要由CPU来计算产生。一方面，CPU的计算速度通常比GPU慢；另一方面，这些数据需要传送到GPU上以供OpenGL显示使用。鉴于此，当数据量很大时，CUDA和OpenGL的混合使用效果明显。

本发明的CUDA和OpenGL之间的互操作是通过OpenGL的fBO对象实现的。在OpenGL中首先需要创建fBO对象并把数据传送到对应的fBO中。CUDA使用fBO数据的过程其实就是三部分：获得数据的控制权->修改数据->将控制权交回OpenGL进行绘制。具体实现步骤如图3所示。

图4和图5显示了使用传统方法绘制的地震数据剖面图形和使用本发明方法绘制的地震数据剖面图形，比较图4和图5可以看出，本发明的方法保证了海量地震数据的地震剖面的实时绘制。

本发明提出的方法和基于传统绘制引擎的绘制方法的绘制时间对比如表1所示：

表1

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (2)

1.一种地震数据剖面图形绘制方法，其特征在于，具体包括如下步骤：

1).设每一道地震数据X方向的间距为traceDistance，Y方向的间距为sampleHeight，输入的地震数据为Date[sampleCount]，输入的地震数据的绝对值最大值为MaxDate，纵向插值的起始点的索引值为sIndex，纵向插值的终止点的索引值为eIndex，则待插值点插值后的地震数据值，即地震数据波形曲线经过的像素点对应的地震数据值InterpData[Index]为：

InterpData[Index]＝Date[sIndex]+α*(Date[eIndex]-Date[sIndex])，

其中，α＝(Index-sIndex)/(eIndex-sIndex)；Date[sIndex]和Date[eIndex]分别为线性插值起始地震数据值和终止地震数据值；Index为纵向插值起始点和终止点之间的点的索引值；

2).根据下式计算每个所述像素点在特定区域的位置：

$Location\&lsqb;Index\&rsqb;=\frac{PixelDate\&lsqb;Index\&rsqb;}{MaxDate}*\frac{traceDis\tan ce}{2},$

其中，Location[Index]为像素点偏离中心位置的像素宽度，Location[Index]的正负分别表示该像素点位于中心位置的左侧或右侧；PixelDate[Index]是地震数据波形曲线经过的所有像素点对应的地震数据值；其中,所述特定区域是指毎道地震数据绘制后所占的像素区域，该像素区域的宽度是traceDistance；

3).根据上述步骤1)和步骤2)的计算公式绘制地震数据波形曲线；

4).将毎道地震数据纵向插值生成的像素点和地震数据对应的像素点标记为边界点；判断在所述特定区域内的其他非边界点的像素点与边界点的相对位置关系，将满足填充条件的像素点设置为填充色，其他像素点设置为背景色，从而完成地震数据剖面图形绘制；其中所述插值生成的像素点即地震数据波形曲线经过的像素点；

地震数据剖面图形绘制方法在GPU上的实现过程为：

1)：检测GPU性能参数；

2)：将源图像数据复制到CUDA数组；

3)：指定CUDA纹理属性，绑定纹理到CUDA数组；

4)：创建像素缓冲区对象与显示纹理，CUDA和OpenGL互操作，实现显示。

2.根据权利要求1所述的地震数据剖面图形绘制方法，其特征在于，所述步骤4)中，对于地震数据的正值区域，填充条件是一道地震数据波形曲线所在的像素区域内的所有的像素点位于地震数据对应的像素点和插值生成的像素点的左边部分；对于地震数据的负值区域，填充条件是一道地震数据波形曲线所在的像素区域内的所有的像素点位于地震数据对应的像素点和插值生成的像素点的右边部分。