基于Python的地震数据可视化

河南省航空物探遥感中心 师天祺

在地震勘探行业中，地震数据的读取与显示是地震处理与解释的关键一环，随着浅层地震勘探、工程地震勘探的兴起，基于Linux的地震处理、解释系统已经无法适应快速发展的局面。本文详细介绍了SEGY数据格式，并重点分析了SEGY数据的难点，借助Python实现了地震数据的可视化，满足了施工现场快速实现地震数据的处理与成像的要求，对现场工作有一定指导意义。

地震数据的读取与显示是地震处理与解释的关键一环，可视化技术利用计算机技术，将数据进行处理，转为显示在屏幕上，通过交互处理挖掘数据之间的规律[1]。可视化技术以被运用到多个领域，特别是在找煤、找油、找气、找断层上取得了重要的意义[2]。但是国内外上通常使用Linux主机作为地震数据的处理、解释的载体，笨重的设备、高昂的生产成本无法适应如今工程物探领域短平快的发展要求。为适应快速发展的工程物探技术，本文借助Python语言内置函数库及Matplotlib绘图库，在详细的分析SEGY格式数据基础上，实现地震数据的读取以及波形的显示，为野外工作起到了指导作用。

1 SEGY数据结构

SEGY是以地震道为单位进行存储，是由SEG协会提出的标准磁带数据格式之一，它是石油勘探行业地震数据的最为普遍的格式之一。标准SEGY数据文件件一般包括三个部分，格式如下：

(1)EBCDIC文件头。由40个卡组成，长度为3200个字节，用来储存对地震数据进行描述的信息，通常需要先转换ASCII码才能进行显示。(2)二进制文件头。长度为400个字节，采用大字节序储存，通常使用32位、16位整型，共包含有32项信息，用来储存描述SEGY数据的关键参数，包括数据格式、采样点数、采样间隔、采样编码等，正确的读取文件头信息将直接影响后续地震道数据的读取。(3)地震道数据。地震道数据由两部分组成，每条地震道数据包含240字节的道头信息和地震数据。道头信息一般保存该道对应的线号、道号、采样点数、坐标等信息。地震数据通常是以IBM浮点形式进行保存的。

2 程序设计及难点分析

Python是一种完全面向对象的语言，其优雅的可读性、可扩展性、快速的开发模式及丰富的第三方库，Python越来越受到开发者的喜爱。本文为快速实现SEGY数据的可视化，缩短开发周期，采用Python作为程序开发语言，仅使用Python内置库及Matplotlib绘图库实现了SEGY数据的读取与显示，同时Python其优秀的第三方库，尤其是地震方向的第三方库，也为未来程序的开发提供了一些优秀的思路和扩展。程序完整流程如图1所示。

图1 程序流程图Fig.1 Program flow chart

想要清晰的显示SEGY数据，除需要理解SEGY数据结构外，还需要解决字节序、IBM浮点读取等一系列问题[3-5]。

(1)字节序。大字节序和小字节序是数据储存规则，小字节序是高位数据储存在内存高位地址，低位数据储存在低位地址，大字节序是高位数据储存在低位地址，低位数据储存在高位地址。字节序示意图如图2所示。

SEGY数据使用大字节序进行编码，这是因为处理器的结构不同造成的，如IBM机、Unix工作站等通常使用大字节序格式，个人电脑通常使用小字节序编码。在实际读取数据过程中如果忽略字节序问题，那么读取的数据可能截然相反。例如采样点数在大字节序为2500，在小字节序则读取为0，所以在读取字节时，应首先将字节序调整至正确编码才按字节依次读取。

(2)IBM浮点。SEGY数据通常使用IBM浮点，偶尔会使用32位IEEE浮点或一些非常用数据格式进行记录。地震数据通过何种格式进行记录，是由卷头中3225-3226字节储存的信息决定的。通过读取这两个字节的内容，可以判断数据的储存格式，从而确定每道地震记录数据所占的字节空间长度，即数据长度=数据类型（字节数）×采样点数[10]Python默认使用IEEE浮点，两种数据结构存在较大的差异，会对地震数据读取与显示造成影响。为实现IBM浮点的读取，需要首先了解IBM与IEEE浮点的区别[6]。

IBM和IEEE浮点数的结构如下：S代表符号、E代表指数部分、F代表位数部分。

浮点数表示如下：

注：32位IEEE浮点数：A=2，B=127，C=1；32位IBM浮点数：A=2，B=64，C=1。

IBM和IEEE浮点的区别在于位数的有效范围，IEEE浮点尾数最小值是1/(2^23)，而IBM浮点尾数最小是为1/16，而不是1/(2^24)。

(3)地震数据绘制。地震原始数据范围较大且杂乱无章，直接使用原始数据成图未必能做出想要的样子。为了使原始数据经过一系列变换转为图形图像，使信息更容易理解，需要对原始数据进行一定程度的处理[7-10]。本文使用生产的单炮数据（500道接收，1ms采样间隔，采样时间为2.5s）作为示例，本文对数据进行归一化处理，使数据范围在（-1，1）之间并设置道间隔为1，保证了道与道之间波形不会发生相交的情况，如图3所示显示了不归一化和归一化数据的成像差异并与使用Seisee软件的成像效果进行对比。示例单炮图像横轴为道序列号，道与道之间间距为设置为1，纵轴为双程时，绘制单炮数据。

图3 成像效果图Fig.3 Imaging rendering

由图3可明显看出，由于未对原始数据进行归一化处理，地震数值随着双程时的增大而迅速减少，由于数据差值过大，波形明显振动主要集中在前0.2s内，后面由于数值相对较小，波形基本属于平直状态。由于横轴设置的道间距过小，振幅数值远超道间距，成像整体挤在了一起，虽能见波形振动，但以无法分辨为哪一道数据，更无法分辨单炮质量的好坏，可见未对原始数据进行归一化处理难以满足成像要求。经过归一化处理的单炮数据整体图像效果较好，炮点位置清晰可见，初至波可明显分辨，初至波之下可见明显层状反射波有效信息，显示效果较好。使用归一化的数据未造成因为距离震源较近，接收能量较大造成的对远处地震波成像影响的效果，与使用Seisee软件的成像效果类似。由于Seisee使用了AGC增益，效果要好于开发的成像程序，但初步做到的Windows平台下的SEGY数据成像，并且为以后工程地震勘探的定制开发打下了基础。

3 结论

随着个人计算机存储技术和计算能力的提升，使用Windows系统进行地震资料处理与解释已经成为行业趋势，尤其是使用个人电脑在施工环境进行野外施工监控、野外数据处理以及工程地震勘探环境需要快速对数据处理及成像。在Windows环境下，本文在充分掌握SEGY数据的基础上，尝试使用Python脚本编写了程序，正确解析了SEGY数据结构，初步完成了成像系统的设计，并应用在实际生产环境中。虽然在成像效果上与主流成像软件还有一定的差距，但是打下了坚实的基础，下一步建议在滤波、AGC上进行研究，增加多种成像模式如灰度、变密度等可视化技术，使之达到主流成像软件的效果。