SEGY格式地震数据的三维可视化

周文辉 石 敏 朱登明 王兆其

1(中国科学院计算技术研究所前瞻研究实验室 北京 100190)2(华北电力大学控制与计算机工程学院 北京 102206)

0 引 言

地震数据三维可视化技术是近年来地质解释的一个重要发展方向[1]，其利用计算机图形学等领域的技术，通过对原始地震数据进行预处理，将其转化为能够反映地质信息的三维地震体数据，并通过一些绘制算法，将地震数据绘制成三维地震体模型[2]，显示在计算机屏幕上。利用交互技术查看各个维度上地质切片信息，挖掘数据间的关系，为断层识别等一些地质勘探活动提供线索[3]。地震数据三维可视化技术的广泛应用，很大程度地降低了地震勘探的成本以及风险，可视化技术高质量地还原了地质的构造信息，提高地震勘探效率的同时提高了勘探的准确性。

本文以SEGY格式的地震数据为例，针对三维数据场的地质解释，研究地震数据的可视化技术。并基于OSG体绘制引擎，实现了SEGY文件的读取、解析以及三维地震数据的可视化。

1 SEGY地震数据文件处理

在地震勘探中，采集到的野外地震勘探数据在大型的工作站中得到进一步的净化，经过一系列处理得到地震数据[4]，将其以道序方式排列，存储为SEGY数据格式文件，SEGY格式存储的地震数据文件目前在地震勘探领域使用得最为广泛[5]。SEGY格式文件按照二进制方式组织数据，其中包含着多种格式的道头，且工作站采集到的数据按照大端格式存储，而地震勘探研究人员通常使用的个人计算机是按照小端方式来存储数据。因此，在PC端正确的加载并解析SEGY格式文件，是地震数据预处理的工作重点，将会影响到人们对地震数据的含义理解及地震数据的可视化效果。

1.1 SEGY地震数据获取

地震勘探通常指反射波地震勘探，目前使用最多的3D地震勘探技术通过人工激发地震波，地震波传播到地层分层处将产生不同强度的反射波。产生的反射波信号将由地面测线上的检波器进行接收并记录在磁带上，得到原始的地震数据，如图1所示。在图1中，横轴表示测线，纵轴表示深度，图中曲线表示每个检波器所采集到的反射波强度数据，一条测线上所有检波器采集到的数据集合又称为一道地震数据。

图1 一道地震数据

由于原始数据中存在着大量的干扰波，因此还需要对地震数据进行预处理(空、废道删除替换等)、常规处理(矫正、振幅补偿等)、目的层特殊处理(属性提取、时频分析等)三个阶段的数据净化，最终得到高精度、高分辨率、高保真的SEGY格式地震数据。

1.2 SEGY数据格式

标准的SEGY文件由两部分数据组成：第一部分数据为记录数据，即两类文件头数据：3 200字节的文本文件头；400字节的二进制文件头。不同文件头的数据均采用不同编码，为了获取地震数据体信息，需要对它们进行格式转换[6]。第二部分数据为地震道数据，即成对出现的240字节地震道头数据与地震数据，其存储着能够反映地质信息的反射波数据，其结构如图2所示。

图2 SEGY格式文件结构图

1.3 SEGY数据类型转换

1.3.1 3 200字节文本文件头格式转换

SEGY文件中，3 200字节文本文件头共包含40个参数，记录了地震数据采集时的外部环境、设备信息等(如设备参数，采集时间)，其数据格式为EBCDIC编码，将其转化为常用的ASCII编码有助于地震解析人员理解文件信息。EBCDIC编码与ASCII编码间可通过一个长度为256字节的映射表进行相互转化，映射表如图3所示。

图3 EBCDIC编码与ASCII编码对应表

1.3.2 400字节二进制文件头格式转换

SEGY文件中400字节二进制文件头共包含32个项参数，记录了地震体的一些测量信息，包括地震数据的数据存储类型、采样数、卷数等。为了正确存储400字节二进制文件头中的数据，需要根据文件头格式定义一个结构体，随后将二进制流的内容拷贝读入结构体中。由于地震采集到的数据是按照工作站大端格式存储，因此对直接读出的数据进行大小端互换得到小端格式，方法如下：在一个循环中，将每个大端格式的字符数据按照其中每个字节的位置，进行第一个字节和最后一个字节互换、第二个字节和倒数第二个字节互换，以此类推，直到循环次数等于字符长度的二分之一时，完成大端格式向小端格式的转换。

1.3.3 地震道数据格式转换

地震道数据分为前240字节的道头数据，与其后的地震数据。位于每个地震道数据前的道头数据，记录了地震数据的采样信息，如采样间隔、震源点号、道集号等。对于地震数据来说，比较常用的数据编码有IBM编码以及IEEE编码。IEEE float类型数据是PC常用的存储格式，具有较好的可读性，在许多领域有着很好的应用。由于两种格式的浮点数位数存在较大的差异，其计算公式如下：

V=(-1)s×M×A(E-B)
M=C+F

式中：V是最终数值，S表示符号，M为尾数，A为基数，E是指数，B是指数偏差，C是尾数常数，F是分数。A、B和C是常量，数值由浮点体系结构决定。地震道数据中，32位IEEE浮点数，A=2，B=127，C=1；32位IBM浮点数，A=16，B=64，C=0。将IBM格式数据转化IEEE格式可根据上述公式进行大量的移位处理，使得IBM格式数据恢复成IEEE格式，存入地震数据内存。

2 三维地震数据体绘制

体绘制指根据三维体数据产生二维图片的一种技术，其在地震数据可视化中应用广泛[7]。在进行三维地震数据体绘制时，由于地震数据规模较大，且数据结构复杂具有空间性，绘制过程涉及大量数据的读取解析、实时绘制等技术[8]。地震数据的存储结构、绘制算法的选择都将会影响到最后的绘制效果及性能。

2.1 三维地震数据结构及特征

来自野外勘探采集得到的三维地震数据，通常是在空间中均匀取样而形成的一个离散的数据体。在绘制过程中需要通过一些数值求解的方式(如插值)进行离散点间的数据恢复，选择数据插值算法时需要注意离散数据间的关联性。

地震体数据采集方式保留了数据的空间性,每一位置点数据都有一个唯一的行号、层号以及列号。SEGY格式文件中数据的空间位置由对应的X_Line纵轴方向、Y_Line横轴方向、Time时间轴方向构成。可以通过三维网格的空间位置，快速访问三维地震数据体中某一位置上的地震强度数据。

2.2 三维地震数据存储结构

地震数据的大小随着采样的精度增加而增加，一般来说，SEGY格式文件存储着大范围区域的地质信息，其大小可达到GB以上。因此在对地震数据进行体绘制、切片处理等操作时，地震数据的载入、提取将会耗费大量的时间，设计一种合适的存储结构将大大改善绘制效率。

由于SEGY格式地震数据按照地震道顺序依道进行存储，在读取地震数据文件时，将数据体进行纵向(X_Line方向)切割，即按照地震道对地震体数据划分。在数据读取过程中，每次读取一个地震道数据，且仅读取需要可视化的地震道数据，将可以大大降低一次性读入内存中的数据量，有利于文件的并行读取。将一道地震数据存储在一个连续的数组中，通过定义缓存可以获得单片地震数据存储空间，随后使用一个指针数组存储所有单片地震数据的起始地址，实现对地震数据的读取。

2.3 应 用

2.3.1 体绘制算法比较

近年来，地震数据体绘制的研究主要集中在体绘制算法及算法的GPU硬件加速上。体绘制算法方面，经典的地震可视化算法有光线投射算法，光线追踪算法以及图元绘制法[9]。对于光线投射算法与光线追踪算法而言，要求体素存在透明度属性，根据光线穿过的体素的透明度进行颜色叠加，最终显示在屏幕上。对于大规模的可视化数据而言，光线投射算法、光线追踪算法均可通过提前终止、跳过空白空间、分布式[10]等方式进行绘制的加速。但是在地震解释上，仅需要查看数据体表面结构，具有透明属性的地震数据将会叠加不同深度的颜色信息，导致无法正确描述当前位置的地质信息，同时因为三维地震数据无需包含透明度信息，不适合上述两种算法。GPU硬件加速方面，GPU强大的并行运算能力能够很大程度减少绘制的时间，提高绘制的效率。GPU体绘制加速算法相比于传统算法，优势在于算法同时有着许多计算单元在进行并行计算。在这个基础上，许多学者将传统算法迁移到GPU上，同时产生了如空白空间跳过(empty space skipping)算法、光线提前终止(early ray termination)算法等进一步提高了体绘制的效率。

在进行地震数据绘制时，本文使用图元绘制算法，图元绘制算法利用几何体图元，将体素相连成片，体素间像素按照插值的方式进行填充。地震切片含有大量的体素，使用图元法进行绘制可选取连续填充四边形串图元(GL_QUAD_ STRIP)，连续填充四边形图元由无限多个点连接成一个密闭的几何体，几何体内部根据顶点颜色进行插值形成。在绘制地震几何体时，将几何体体素组织在一个连续填充四边形串图元中，既可方便找到需要绘制的体素，又可避免对体素的重复绘制,可大大减小体素的重复计算，提升绘制效率。选取OpenGL中的连续填充四边形图元对象，据图元绘制次序将体素面片填充入图元对象中，最后将图元对象添加到绘制场景中，在内存中进行绘制。

2.3.2 地震强度与颜色空间映射

常用的地震解释技术通常使用颜色空间作为地震强度的表达媒介：颜色传递函数将地震数据中的反射波强度数据映射到颜色空间中，为地震数据中每一个体素上色，然后通过体绘制算法构建出几何体模型。具体的映射方式如下：地震数据颜色传递函数[11]可将SEGY格式文件中地震强度数据映射到地震解释常用的蓝、白、红颜色空间中，实现体素的颜色赋值。地震解释使用的蓝、白、红颜色空间三个颜色通道分别代表了地震波正向振幅、零振幅、负向振幅。且颜色数值的大小范围在[-1,1]之间，因而需要对地震数据进行归一化处理，使得地震数据数值与颜色数值一一对应。

SEGY文件中保存的地震数据代表反射波强度，正负值表示地震波振幅方向。地震波的振动特性与数值大小无关，将地震数据进行归一化处理，映射到[-1,1]之间，保留了数据值的相对大小与正负号将不会改变地震波所包含的地质信息。归一化的实现方法为找出地震数据中绝对值最大的值，将其映射为数值1，其他数据按照相同的映射方式将地震数据映射到[-1,1]范围间，实现归一化处理。归一化处理在保证数据不失真的前提下，能够缩小数值范围，减少数据量，提高绘制效率。将单片地震数据归一化后，得到的数据直接作为颜色值通过绘制函数直观地显示到屏幕上，得到的结果如图4所示(颜色不能显示)。

图4 归一化数据绘制

颜色传递函数将归一化后的地震数据通过与颜色空间进行映射达到给体素上色的目的。归一化处理及RGB颜色赋值可通过以下操作一次性完成：

R=intense/maxValue×step(0,intense)
G=intense/maxValue×(2×step(0,intense)-1)
B=intense/maxValue×(step(0,intense)-1)

式中：intense表示地震强度数据，step(0,intense)函数指当intense>= 0(正向振幅)时step() 函数返回1，否则(反向振幅)返回零。

对得到的RGB值进行一次取反操作，并存入一个三维向量中完成对体素的颜色赋值：

pixColor=vec3(1.0-R,1.0-G,1.0-B)

对于单片地震数据，经过颜色传递函数的颜色赋值，将噪音置成白色，得到含有颜色信息的数据体，将其通过绘制函数绘制得到的结果如图5所示。

图5 含颜色信息的地震数据绘制图

3 基于GPU的绘制算法优化

传统的绘制算法在CPU上进行地震数据的可视化绘制，由于地震数据规模较大，且每个像素点在CPU上串行绘制，硬件资源未得到充分利用，绘制效率低下。GPU在设计上采用了数量众多的计算单元及流水线，采用这种结构具有极大的吞吐量，适用于类型统一，无耦合的大规模并行计算。对于地震数据可视化这种计算密集且易于并行计算的程序，迁移到GPU上运行能够发挥出GPU巨大的算力。

GPU图形绘制管线专门描述了GPU在数据可视化上的整个渲染过程，如图6所示，主要包含以下三个阶段：应用程序阶段，主要在CPU内进行图元设置等操作；几何阶段，主要在GPU上执行顶点、坐标计算；光栅阶段，主要为体素进行上色等操作。经过整个图形渲染过程，将三维数据体绘制到二维屏幕上。

图6 GPU图形绘制管线

GPU在进行图形绘制时通过并行计算对绘制过程进行加速[12]，当程序具有以下性质时，将其迁移到GPU上将能够最大发挥出GPU计算上的优势：(1) 算法易于并行执行，算法间耦合度低，数据耦合度低，可以同时进行计算。(2) 计算密集型程序，程序中涉及大量同类型数据的重复计算。

在地震数据体绘制算法中，将每个体素按照图元顶点顺序组织成一个地震切面。图元内部体素与体素间不存在相互依赖关系，顶点之间像素可通过插值得到。图元与图元之间同样不存在依赖关系，因此该算法具有很好的并行条件，易于并行计算。在进行图元绘制时，需要进行体素顶点坐标计算、体素地震强度的颜色转换、体素光栅化以及体素间的颜色插值等大量的重复计算，参与计算的数值类型统一，均为IEEE float类型，属于密集型程序，具有很好的数据并行性。基于此，将图元绘制算法移植到GPU上能够大大提升图元绘制算法的性能。

基于GPU的地震数据图元绘制算法的执行流程如图7所示。

图7 地震数据体绘制流程

首先将三维地震数据体(SEGY文件解析后得到的数据体)作为输入，通过颜色传递函数将体数据中地震强度信息映射到颜色空间，得到体素的颜色信息;同时，通过三维地震体素坐标信息，得到体素在空间中的相对位置坐标。选择连续填充四边形串图元(GL_QUAD_ STRIP)，将体素按照图元前后顺序依次载入内存中。定义着色器，将地震颜色信息以及坐标信息传入GPU中，GPU对所有体素并行执行一遍着色器程序。由于体素按照连续填充四边形图元(GL_QUAD_ STRIP)的定点顺序的方式传入GPU，在绘制地震几何体时，将几何体体素组织在一个连续填充四边形串图元中，四边形串内部根据顶点颜色进行插值生成，这种方式既可方便找到需要绘制的体素，又可避免对体素的重复绘制,可大大减小体素的重复计算；此外不同面片间可进行并行加速绘制，提升了绘制效率。

GPU在进行并行绘制时，对每个体素均需要执行一段shader程序，因此优化shader程序也有利于提升GPU绘制的效率。shader(着色器)程序是能够在GPU上执行的一段绘制程序，使用shader进行GPU编程时，需要编写顶点着色器程序，以及片段着色器程序。其中顶点着色器用于处理顶点数据，将顶点数据载入顶点缓存中，顶点着色器对缓存中的每一个顶点进行渲染；片段着色器则是对片段操作的一种通用功能的可编程方法，其输出是片段的颜色和深度值。shader程序进行的优化如下：在顶点着色器中直接传入地震强度数据、体素坐标数据以及场景变换计算矩阵(用于场景变换时实时更新场景)，减少shader进行数据处理而花费的时间，然后通过当前体素坐标位置与场景变换计算矩阵的乘积，得到更新后的体素位置；简化传递函数将颜色空间归一化在[-1,1]之间，通过传递函数将地震强度数据映射到颜色空间的方式获得体素的颜色，减少计算颜色值所需的时间，随后等待在GPU图形管线光栅阶段进行渲染。

4 基于OSG的地震数据三维可视化

4.1 基于OSG的三维体绘制

OSG是目前地震数据可视化领域中应用较为广泛的一种3D图形开发引擎，它有着丰富的第三方插件，在体绘制效率、三维交互上均有相应的优化。

地震体绘制过程涉及到大量数据的读取解析、实时绘制等技术，OSG三维可视化引擎对体绘制算法进行了大量的优化操作，简化绘制过程，基于OSG的三维体绘制流程图如图8所示。

图8 地震数据体绘制流程

首先将SEGY格式文件载入计算机内存，提取出图元绘制所需的三维地震体表面地震数据。根据控制场景状态控制图形的绘制，因此需要将体素颜色信息以及坐标信息通过传递给场景中的叶子节点。同样，需要将shader属性以及设置绘制图元均加入到同一个叶子节点。OSG中shader对象均通过一个专门的类管理，通过设置子节点的shader状态属性来定义子节点，最后将子节点添加到OSG场景中进行绘制。

此外，当用户对三维可视化得到的地震体模型进行三维交互(旋转、切片等)时，OSG需要完成绘制场景的更新、体素重新绘制等操作，该内容在4.2节中进行介绍。

4.2 三维交互及场景控制

OSG中定义了视景器，实现对三维交互的控制，视景器中包含漫游器类、事件处理器类、场景类以及摄像头类。其中漫游器，用于实现交互的场景漫游(包括坐标变换、体素更新等)；事件处理器，处理鼠标，键盘响应事件；场景，视景器对应的场景根节点；摄像机：实时展示绘制场景。

对于地震体数据可视化而言，交互过程中，场景中三维物体的自由旋转、放大缩小等操作，均会引起相机视口的变换。随着视角变换，场景更新需进行重采样，以及采样点间的线性插值操作。重采样需根据相机的视口变换得到变换矩阵，利用变换矩阵确定采样点在三维体数据中的位置，随后完成从世界坐标到物体坐标的投影，实现原理如下：通过camera类获得相机的视点矩阵；通过将本地坐标转换为世界坐标的方式，获得当前模型的位置矩阵；最后将相机视点矩阵、模型位置矩阵以及相机类中的映射矩阵相乘，其结果即为交互前后顶点位置变换的转换矩阵。

4.3 三维地震体交互

通过体绘制将SEGY格式文件转化为三维地震体，对三维地震体模型进行观察时，通常使用的三维交互方式有局部地震体数据放大，观察局部数据间的关联关系。对地震体的放大缩小操作通过鼠标滚轮控制，其放大的效果图如图9所示。

图9 三维地震体局部放大图

对地震体模型进行切片展示同样是三维交互中的一个重要交互方式。对地震体在X_LINE、Y_LINE、TIME三个方向上进行切片[13]，以便地震解释人员观察地震体内部数据间关系。最后根据地震解释人员的专业背景知识对地质构造信息、石油信息进行还原。对地震体进行三个方向的切片的效果图分别如图10(从左到右分别为X_LINE，Y_LINE，TIME方向)所示。

图10 三维地震体切片

5 实验结果与分析

本文使用一组地震体数据对加速前的体绘制算法与改进的体绘制算法进行对比分析。实验中采用的计算机配置如下：CPU为Intel i5-4210M、显卡为Intel HD Graphics 4600、实现模式为1 920×1 080、屏幕刷新率为60 Hz、内存容量8 GB、操作系统为Window 7 64位旗舰版。针对这组参数，本文设置了三组数据量大小分别为150 MB、350 MB、550 MB的SEGY数据文件，针对体绘制过程中三组重要的性能参数，即数据文件解析时间、模型绘制时间以及三维交互时的每秒所绘制的帧数，分别对加速前的体绘制算法与改进后的体绘制算法进行比较分析，并计算得出性能提升比。实验结果如表1-表3所示。

表1 150 MB文件体绘制算法性能指标

表2 350 MB文件体绘制算法性能指标

表3 550 MB文件体绘制算法性能指标

从表中可以看出，在处理大小为150 MB的地震数据时，本文算法文件解析时间与模型绘制时间均有较明显提升，每秒所绘制的帧数提升较小，原因是数据量较小，仍在显卡承受范围内。150 MB地震数据文件绘制结果如图11所示。

图11 150 MB三维地震体

对于350 MB的地震数据来说，本文算法相比较于传统算法在文件解析时间、模型绘制时间以及每秒所绘制帧数上均有较大提升，整个绘制过程大约需要2 880 ms，相比传统算法绘制需要的6 222 ms，有了比较大的提升，同时FPS保存60 Hz，保证交互时的流畅性。350 MB地震数据文件绘制结果如图12所示。

图12 350 MB三维地震体

当数据量增加到550 MB时，本文算法性能得到最大的发挥，算法的并行性在处理大数据时具有较大的优势。在每秒所绘制帧数这个指标上，本文算法提升了335.71%，有效保障了在进行三维交互时操作的流畅性。550 MB地震数据文件绘制结果如图13所示。

图13 550 MB三维地震体

6 结 语

随着计算机硬件的提升，地震数据三维可视化技术发展迅速，在很多领域已有较为成熟的应用。本文以SEGY格式地震数据为例，针对三维数据场地震强度数据的性质，研究了地震数据的可视化技术。分别对体数据加载、地震数据格式转换、体绘制算法、GPU并行可视化、OSG场景控制进行了深入研究，给出了基于OSG框架下的SEGY格式地震数据的三维可视化实现，为地震解释人员了解地震数据、分析油藏情况提供了一个直观、便捷、高效的途径。