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]。变面积图就是为了更好地表现同相轴,在绘制时通常将地震数据正值区域或负值区域用特定的颜色来填充。
变面积显示有波形加面积、正极性、双极性显示这几种方式。假设平衡位置的右侧为正值区域,左侧为负值区域。波形加面积是对正值区域进行填充,负值区域以曲线表示;正极性是将正值区域填充,负值和零值区域不做处理;双极性显示是将波峰和波谷分别填充不同的颜色。
填充的关键是要判断数据点与平衡位置的关系以及相邻数据点之间的关系。图3波形中10个特殊点详细的表现了它们之间的关系:
第一种情况是相邻两点一个在平衡位置一个不在平衡位置,这样只需从不为0的点向平衡位置作垂线,垂足与这两个点组成一个三角形。第二种情况是一个点为正值一个点为负值,两点连线与平衡位置有个交点P,同样分别作垂线,垂足、数据点和交点组成2个三角形。第三种情况是相邻两点在平衡位置同一侧,分别作垂线后两点与垂足之间围成一个梯形。第四种情况两点都在平衡位置,这样不需要作处理。
在绘制波形变面积图或正极性图时,如果三角形或梯形在平衡位置的右侧则将其填充,在左侧时不处理。在绘制双极性图的时候才将处于平衡位置左侧的三角形或梯形填上与右侧不同的颜色。
GDI+的Graphics类为我们提供了绘制提供了区域填充的方法Graphics::FillPolygon(Brush*brush, PointF*points,INT count,FillMode fillMode),创建一个Graphics和Brush对象,顶点一般只有三个或4个而且摆放是规则的,故ALTERNATE和WINDING两个填充模式效果一样。将上面求得的多边形的顶点存在一个数组中,向函数传递参数,则可完成面积图的绘制。
图4是填充后绘制出来的三种面积图,从图中可以清晰的表现出自同一界面反射波的同相轴,大致反映层位的基本走向。
3.3 绘制变密度图
变密度图是利用不同的颜色来表示地震数据的大小,用一个连续变化的彩色谱表示地震数据的不同振幅。方法是依据显示比例以采样点为中心画小矩形进行区域填充,根据该样点值的大小进行着色填充。
变密度图的关键是对颜色进行分级。可以自己定义最大正振幅、最大负振幅和零值是对应的颜色,数值从最大到零时从一种颜色逐渐过渡到另一种颜色。以正值为例:设地震数据最大正值为MaxRight,如果最大正振幅对应的颜色为(r1,g1,b1),零值对应的颜色为(r2,g2,b2),则当前数据Data[i]所对应颜色的红色分量r=r2-Data[i]*(r2-r1)/MaxRight,绿色分量g=r2-Data[i]*(g2-g1)/MaxRight,蓝色分量b=b2-Data[i]*(b2-b1)/MaxRight,负值时对应的颜色同理。
确定当前数据点所在的道的基线坐标X1,下一道的基线坐标X2,获得当前采样点的屏幕坐标值yl和下一采样点的坐标值y2,以(x1,y1)为左上点和(x2,y2)为右下点绘制矩形,用采样点振幅值换算出的像素值填充矩形区域,这就是绘制变面积图的基本算法原理。
以上面的算法进行填充效果不会很理想,当道间距和时间间隔放大一定的比例时,就会看到一个个小矩形,而且相邻矩形之间的颜色过渡也不好,可以利用双线性插值方法来改善绘图质量,实现的过程如下:
利用上面算法得到的矩形区域以及填充该区域的颜色,求取填充矩形四个顶点的颜色值(R,G,B)。顶点颜色的求取如图5所示相邻四个矩形有一个共同的顶点,该顶点的颜色值为相邻四个区域颜色值的平均值。之后对矩形左右从上到下进行纵向颜色插值,再根据纵向插值得到的点,对矩形上下从左到右进行横向颜色插值。这些插值点将原来的矩形分割成一个个更小的矩形,每个点对应一个小矩形,用这个点代表的颜色值将其对应的矩形填充即完成整个绘制过程。插值点的数目并不是越多越好,数目大则运算量加大,会影响绘图速度,数目少图像的改善效果不明显。可以根据显示需要合理设置。
图6是绘制出来的变密度图。从绘制结果可以看出数据值与颜色之间有着很好的对应关系,相邻矩形之间的颜色值过渡也很平滑自然。从意义上来看变密度图跟面积图相比更能在细节上精确的反映层位的基本走向并展示整个地层的结构。
4 结束语
地震的剖面显示是地震数据可视化的基本形式,它为用户提供精确的图形显示、交互以及多种信息的综合表达。本文解析了地震勘探数据SEG-Y文件格式,利用VC++在正确读取地震数据的基础上,实现了地震剖面几种常见表现形式的绘制,详细的介绍了绘制的原理和绘制过程。从绘制结果来看,比较真实的反映了数据原貌。可以从整体上看出地震数据的分布及走势,而且有助于对地层结构有个整体的了解,进而为后期的解释及油藏的可能存在位置提供依据。
参考文献:
[1] SEG web site:http://www.seg.org/下载SEG Y标准rev 0、rev 1.
[2] 李捷,王辉等.基于微机的三维地震资料可视化系统的数据格式及数据存取[J].大庆石油学院学报,2005,29(6):25-28.
[3] 肖汉.地震数据的可视化研究[D][湖南大学硕士学位论文].湖南:湖南大学,2007:22-25.
[4] 陆基孟.地震勘探原理[M].山东:石油大学出版社,2004:125-131.
作者简介:刘继承(1970-),男,教授,主要从事油气信息处理技术方面的研究。