基于VC++的SEGY数据格式地震剖面图绘制

基于VC++的SEGY数据格式地震剖面图绘制

关键词：地震数据，地震剖面，显示波形图，变密度图

1 引言

地震勘探采集的地震数据一般以地震道为单位进行组织,采用SEG-Y文件格式存储。SEG-Y格式是由SEG (Society of Exploration Geophysicists)提出的标准磁带数据格式之一,它是石油勘探行业地震数据的最为普遍的格式。读取数据后经过处理可以得到有关地下构造和岩性等大量有用信息信息。但是原始数据数据量庞大而且杂乱无章,为了使大量的数据变得直观、易于理解,揭示数据间的内在联系,必须对原始数据经过一系列变换转换成图形图像,将信息以更容易理解的形式展现出来。VC++功能强大、界面友好且操作方便,因此目前PC机上通常在VC++环境下对这种格式存放地震勘探数据进行处理并实现可视化显示。笔者在做地震数据可视化时利用VC++6.0编写了地震数据的读取处理程序及绘图程序。 2 读取SEG-Y数据

2.1 SEG-Y格式地震数据的文件结构[1,2] 标准的SEG-Y文件结构如表1所示:前3600字

节为文件头通常被看做一个整体,紧跟其后的是数据体。非标准的SEG-Y文件可能没有前3600字节的文件头。

文件头有3600字节分为两部分,第一部分长度为3200字节,以EBCDIC编码记录地震数据采集系统的有关信息。共有40行,每行长度为80个字节,每个字节表示一个字符,EBCDIC码需先转化成ASCII码,才能在微机中正确显示。第二部分为二进制格式,长度400字节,数据类型为32位或16位的整型,记录数据体信息,如每道的采样点数、采样率、数据格式等[3]。

数据体由多个数据道组成。每道数据分道头、道数据两部分。道头是二进制数据,长度240字节,数据类型:32位或16位的整型,记录采样点数、采样间隔、共深度采集点CDP(Common Depth Point)号、XLine号、YLine号以及坐标信息等。地震道的数据长度为采样点数×4个字节。工作站SEG-Y数据的存储格式一般为IBM格式,微机处理的SEG-Y数据存储格式一般为IEEE格式。通过读取二进制文件头中第3225～3226字节存储的信息可以判定数据的存储类型。 2.2 读取地震数据

地震数据的正确读取是数据处理的前提,任何一位参数的偏差都可能导致整SEG-Y数据文件出现内部参数

矛盾而无法进行数据解读。了解SEG-Y格式地震数据的文件结构后对读取位置定位,读取完成存储为后面的处理作准备。微机版和工作站版SEG-Y数据结构一样,但是数据的存储格式不同,微机中的数据是低字节在前,高字节在后。在微机的操作系统环境下,对SEG-Y格式的数据进行处理时,首先必须对数据进行高位、低位互换处理。

读取SEG-Y数据用到的类如下:在程序中定义了2个结构体:ST_real_info和ST_trace_info,前者用来定义SEGY卷头文件中400字节数据信息代表的变量,后者用来定义道头文件中头240个字节数据信息代表的变量。定义了读取SEG-Y数据的类CReadSegyData,主要功能有:Ebcdic2Ascll ()完成EBCDIC码向ASCLL码的转换;ibm _to_ float ()完成地震数据IBM格式向IEEE数据格式转换;GetRealVar()完成卷头信息的读出;GetOneTraceInfo()完成道头数据的读出;GetOneTraceData()完成道数据的读出等。

读道头和道数据时先要判断文件按有无文件头,如果有则第N道道头和道数据读取位置分别为: info_pos=3600+(N-1)(240 +SampleNum*DataSize) data_pos=3600+(N-1)( 240+SampleNum*DataSize)+240

其中SampleNum是采样点数,DataSize是一个数据点所占的字节数。如果没有卷头,则上面的式子中就没有3600,读文件时用fseek()定位,用fread()读取。读完后判断SEG-Y格式,如果是IBM格式则必须要转化。 3 地震的剖面显示

常见的地震剖面显示方式有波形显示、变面积显示和变密度显示,其中变面积显示又有三种方式:波形变面积显示、单极性显示和双极性显示[3]。图1是绘图参数设置对话框,可以根据显示的需要对各种参数进行设置,比如道与道之间的间隔nTraceInterval和采样点之间的间隔nTimeInterval、起始道m_StartTrace和终止道m_EndTrace、起始时间m_StartTime和终止时间m_EndTime等。为类型设定组合框中的四个单选按钮添加变量m_DrawOption,项目创建后系统自动为我们创建了主框架CMainFrame类、应用程序App类,文档Doc类和视图View类,在视图View类的实现文件中添加单击消息处理函数,单击波形图单选框则将

m_DrawOption赋值1,单击波形加面积则将其赋值2,绘制波形图的函数DrawWaveArea()、绘制波形加面积的函数DrawWaveArea()、绘制双极性的函数DrawChromatoGraph()和绘制变密度的函数DrawVarDensity()也在View类中,在绘图时对m_Draw Option值做判断,如果为1,则调用DrawGraphWave()函数绘制波形图,如果为2则调用DrawWaveArea()函数绘制波形加面积图,依次类推。

在Windows操作系统下,绝大多数具备图形界面的应用程序都离不开图形设备接口GDI(Graphics Device Interface),我们利用GDI所提供的众多函数就可以方便的在屏幕、打印机及其它输出设备上输出图形。在GDI+模式下编程要gdiplus.dll动态链接库的支持,将其复制到程序当前目录,在应用程序App类中添加头文件gdiplus.h,库文件#pragram comment(lib,“Gdiplus.lib”),增加成员变量保存GDI+初始化后在应用程序中的标识,重载初始化函数和修改退出函数。

3.1 绘制波形图

地震数据值记录了反射波在检波点离平衡位置的振动情况,它的大小代表了振幅的强弱,正负表示了振动的方向。波形图是将多道地震数据在一个平面内以曲线的形式表示出来。波形图能够清晰的反应反射波的动力学特征(如振幅、频率和波形等)。

图形区域与离客户区域左右边界的距离sx,sy。道数据储存在数组one_trace_data中,则第i道第j个采样点的横坐标位置X=sx+i*nTraceInterval+one_trace_data[j],纵坐标位置Y=sy+t*nTimeInterval,下一个采样点的横坐标位置X= sx+i*nTraceInte rval+one_trace_data[j+1],纵坐标位置Y= sy+(t+1)*nTime Interval,先根据起始时间和终止时间循环读取一道数据,确定坐标(X,Y)后依次连接起来,得到一道数据的波形图,再根据起始道和终止道循环读取并绘制每一道。

如果时间间隔变大,直接将点与点相连曲线就变得不光滑看起来像折线,可以绘制贝塞尔曲线来改善。贝赛尔曲线是一种以逼近为基础的参数多项式曲线,曲线逼近实际上是一种插值行为。

GDI+的Graphics类提供了绘制贝塞尔曲线的函数:Graphics::DrawBeziers(Pen* pen,PointF*points, INT count),该函数有三个参数:指向画笔的指针、存储控制点坐标的数组和坐标的个数。我们只需将地震数据分组,将该组数据点对应的坐标作为控制点存放在数组中,实例化一个Graphics对象和Pen对象后,调用该函数并传递参数,则由内部封闭的算法来自己完成贝赛尔(Bezier)曲线的绘制。

图2是绘制出来的地震剖面波形图,从图中可以看出地震数据分布情况以及偏离平衡位置的程度,从而了解不同深度地震波的振动情况。

3.2 绘制变面积图

在一系列地震道中,同一反射波的相同相位在相邻地震道上的到达时间是相近的,每道记录的振幅也是相似的,因此相邻地震道之间波形相似,波峰比较靠近。波峰在地震剖面上相互叠套成串,一连串的波峰组成一条线,形成同相轴。同相轴具有一些特点,如振幅显著增强、相邻道波形特征相似等。可以根据同相轴的这些特点来识别地震层位[4]。变面积图就是为了更好地表现同相轴,在绘制时通常将地震数据正值区域或负值区域用特定的颜色来填充。

变面积显示有波形加面积、正极性、双极性显示这几种方式。假设平衡位置的右侧为正值区域,左侧为负值区域。波形加面积是对正值区域进行填充,负值区域以曲线表示;正极性是将正值区域填充,负值和零值区域不做处理;双极性显示是将波峰和波谷分别填充不同的颜色。