一种基于灰度图的地震变面积波形绘制方法

摘要： 在地球物理勘探中，地震资料计算机显示方法主要有：变密度彩色显示、波形显示、变面积黑白（填充波形）显示等，其中变面积显示方式在数据采集和资料处理阶段使用最为广泛。现有变面积显示技术在查看小范围数据时波形填充显示效果比较好，当扩大显示范围时由于屏幕分辨率限制，只能对波形进行抽希显示，这种现实方法精度不高，显示效果不美观。本文提出一种基于灰度的变面积绘制方法，将像素位置看成一个单元，首先统计要填充的波形落在单元内部的强弱权值，最后根据每个单元内的波形权值设置像素的颜色深浅，即使用颜色的深浅表示波形的强弱。实验表明，本文提出方法绘制的地震剖面显示精度更高，效果更加美观。

Abstract： Currently in geophysical exploration， there are some kinds of methods for seismic waveform display with computer： variable density color display， waveform display， variable area filled waveform display， etc. The variable area filled waveform display effect is good when displaying small range data. However， when the display range expands， due to the limit of screen resolution， the only way is to resample data. This method's accuracy is not high， showing unsightly. This article describes one method for display waveform which is based on gray variable area. Take every pixel as a unit， add up the power of wave in the unit， then use the power to calculate the blackness， which means using the blackness to show the waveform. Experiment shows that the seismic profile displayed with this method has higher accuracy， and the effect is more beautiful.

关键词： 地震剖面；波形显示；变面积剖面；灰度

Key words： seismic profile；waveform display；variable area profile；gray

中图分类号：TP317.4 文献标识码：A 文章编号：1006-4311（2014）29-0208-03

0 引言

在地球物理勘探[1]中，地震资料在计算机上的显示方法主要有：变密度彩色剖面、波形显示、变面积波形填充三种。因变面积剖面能够反映界面的形态，直观性强，外形与地质剖面接近，而且能够突出反射层次，有利于波形分析和对比，所以在地震资料采集和处理阶段主要使用波形、变面积剖面或波形加变面积显示方式查看资料，而地震资料解释阶段则主要使用变密度彩色显示方式来查看地震资料[2-4]。

对于波形显示和变面积波形填充显示方式[5]，当一个窗口范围内显示的地震资料比较多、显示地震道密度很高时[6]，比如一个像素宽度范围就有一个地震道时，如果采用正常显示方法那么整个视图区域将会是漆黑一片，什么也看不到。目前通常的解决方案是对地震道进行抽希显示，即抽取地震道显示，保证显示的地震道有足够的像素宽度来显示。这种方案将会有可能丢弃掉很多关键数据，比如从十几道或者几十道中抽取一道来显示，大量的数据信息被丢弃，导致无法看到完整的数据信息。图2为抽希方式显示的波形效果图，图（a）和（b）分别为每道地震数据最小显示宽度为5个和3个像素宽度的抽希策略显示效果图。

本文提出使用灰度图方式方式解决这一问题，当视图中显示的数据量很大时，视图中的每个像素使用灰度表示该像素位置的波形的强弱，即用灰度图展示地震波形的强弱变化的地震剖面，使用这种方法实现变面积剖面显示，使绘图精度更高，效果更加美观。

1 灰度图绘图原理

经过分析发现，之所以在地震道足够密时要进行抽希显示的原因是：绘制波形时，如果地震道足够密时，每个像素点要么是被绘制成一个颜色，要么是背景色，那么绘制出的图形视觉效果跟预期的相差太多，所以导致显示精度不高。人们真正需要看到的效果是通过某种形式展现出波形的真实强弱，比如通过不同灰度展现出波形的最接近的形状。如图3，（b）为只有一种填充颜色的显示效果，（c）则为使用灰度方式显示的效果，可以看出，灰度方式的显示效果更能反映真实的数据资料。

这种情况只需要计算出每个像素内的波形面积大小，并根据该面积所占像素整个面积的百分比计算出该波形的灰度值大小，最后根据灰度值大小得出该像素的RGB颜色值，即可以绘制出灰度剖面。

还需要考虑另外一种情况，当地震道显示宽度很小时，相邻的波形可能会重复绘制在一个像素中，即有些像素会被多次绘制，如图4中颜色突出的像素部分将会被多次填充。此时这两个像素的灰度值应该是两次绘制波形灰度值的累加值，即两次绘制过程灰度值的和。

2 灰度图绘制

要实现地震变面积剖面绘制过程，必须首先通过视图坐标变换过程将逻辑坐标（地震道索引和样点索引）变换到窗口视图中，然后根据变换逻辑将数据波形变换到对应的像素位置，最后计算出窗口中每个像素的灰度值，将每个像素颜色修改成对应的RGB灰度颜色。考虑到如果逐点设置窗口像素颜色效率会比较低，实现时在内存中申请一个和窗口大小一致的位图，修改位图的每个颜色值，绘制结束后将位图呈现在窗口中。

2.1 地震道绘制 窗口能显示的数据范围只是整个数据范围的一部分，必须允许用户通过缩放和滚动变换窗口的显示范围。为了实现滚动和缩放控制，这里引入四个变量控制缩放和滚动，水平方向缩放比例（xx），垂直方向缩放比例（yy），水平方向显示位置（dx）和垂直方向显示位置（dy）。这四个变量可以用矩阵T来表示：T=xx yydx dy

那么第x道第y个样点波形的中心点位置在窗口中的位置为：

（x′，y′）=T×（x，y）=xx yydx dy×（x，y）（xx×x+dx，yy×y+dy）

这里矩阵中的xx表示每道地震波形所占的像素宽度，yy表示每个样点间隔所占的像素高度。为了将波形绘制在一道宽度内，还需要将波形归一化到（-0.5，0.5），即一道的宽度内。假设波形最大值为Amax，最小值为Amin，则x道、y样点、波形Ai的像素位置为：

y′=yy×y+dy x′=xx×（■-0.5）+dx

为了填充波形，这里还需要一个参数定义填充的基准线，假设波形填充基线为Abase，则基线的像素位置为：

xbase=xx×（■-0.5）+dx

当x′Xbase时填充正波形。

需要注意的是如果yy参数大于1，即每个样点距离大于一个像素，需要对两个像素之间通过某种插值算法计算出垂直方向每个像素位置对应的振幅值，并通过振幅值计算出对应的像素位置，这里采用线性插值方法计算两个波形控制点之间的振幅。假设上一个振幅值为A0，下一个振幅为A1，yi和yi+1之间养点振幅差值公式为：

Ai=A0+■

其中波形线是否显示、正波形是否填充、负波形是否填充是用户配置的，这些位置是否绘制需要根据配置参数确定。

2.2 显示数据范围 在绘图之前，必须先确定需要绘制的地震道范围和样点范围，然后从文件或内存中读取对应的样点值进行会之后。假设窗口视图的水平显示范围为（x1，x2），垂直显示范围为（y1，y2），总地震道个数为TraceCount，样点个数为SampleCount则通过上面公式我们可以计算出需要绘制的地震道索引范围和样点范围：

Trace1′=xx×x1+dx Trace2′=xx×x2+dx

Sample1′=yy×y1+dy Sample2′=yy×y2+dy

此时计算出的范围是浮点值，而且必须考虑相邻道和相邻样点的绘制，所以最终需要绘制的数据范围为：

Trace1=max（0，integer（Trace1′+1））

Trace2=min（TraceCount-1，integer（Trace2′+1））

Sample1=max（0，integer（Sample1′+1））

Sample2=min（SampleCount-1，integer（Sample2′+1））

这里integer为对浮点数取整函数，max表示两个数中取较大者，min表示两个数中取较小者。

2.3 灰度方式绘制 因为灰度[7]方式计算量比较大，效率相对较低，所以在地震道显示宽度足够宽，比如大于5个像素时，不使用灰度方式绘图，使用常规方式绘制波形即可，而当地震道显示宽度小于5个像素时才使用灰度方式绘图。通过上面的方法可以计算出每个地震道的宽度（即波形的整体宽度）：TraceWidth=abs（xx）

这里abs指取绝对值函数。当TraceWidth大于5个像素时，使用正常波形绘制方式，否则使用灰度方式绘制波形。通过试验发现，当TraceWidth小于3个像素时，波形显示已经没有意义，此时不需要显示波形线，只显示波形填充即可。

2.3.1 正常绘制波形步骤

①在地震道（Trace1，Trace2）范围逐道读取（Sample1，Sample2）范围样点数据绘图，对每道地震数据执行②-④步骤。

②计算Y方向每行像素位置的振幅值，当Y方向每个采样间隔的显示高度大于一个像素时，需要对中间的行位置根据2.1中的公式进行插值。然后对每个像素位置的振幅值执行③-④步骤。

③如果当前视图显示波形线，则将x′位置的像素设置为波形线颜色。

④计算需要填充的像素范围：

x1=round（x′） x2=round（xbase）

x1′=min（x1，x2） x2′=max（x1，x2）

将（x1′，x2′）范围内的像素颜色为波形填充颜色，这里min指两个数中取较小的值，max指两个数中去交大的值，round指对浮点数进行四舍五入取整型值。

2.3.2 灰度方式绘制波形步骤

①在内存中申请一个浮点二维数组，该数组的大小和当前绘图的位图大小一致，该数组中的值代表每个像素的灰度值，0表示白色，1表示黑色。将该数组中每个数据值的初始状态设置为0。

②在地震道（Trace1，Trace2）范围逐道读取（Sample1，Sample2）范围样点数据绘图，对每个样点数据执行③-⑤步骤。

③计算Y方向每行像素位置的振幅值，当Y方向每个采样间隔的显示高度大于一个像素时，需要对中间的行位置根据2.1中的公式进行插值，计算该位置的振幅值。根据每个样点所占像素个数计算其垂直方向灰度权值：

Pxi=abs（yy）？叟1：1abs（yy）

然后对每个像素位置的振幅值执行④-⑥步骤。

④根据2.1中的公式计算每行像素位置的振幅值对应的像素位置x′。

⑤如果绘制波形并且TraceWidth大于3并且当前显示配置参数为显示波形线，则计算的灰度值，因为我们希望波形线的颜色要看起来稍微浅一些，所以将该波形的颜色值乘以0.2，作为波形线的灰度值累加到该像素位置的灰度值中：Grayround（x′）=Grayround（x′）+Pxi×0.2

⑥确定需要绘制的开始和结束像素位置：

x1=min（x′，xbase） x2=max（x′，xbase）

x1′=integer（x1） x2′=integer（x2）

这里min指两个数中取较小的值，max指两个数中去交大的值，integer指对浮点数进行截断取整，变换成整型值。

当x1′=x2即波形宽度小于一个像素时，x1′位置灰度为：Grayx′=Grayx′+[Pxi×abs（x1-x2）]

否则对于（x1′，x2′）范围像素分别计算灰度值：

Grayx′=Grayx′+Px■×xi=x1：1-x1+（integer）x1x1

⑦最后对每个像素根据灰度值计算图片的RGB颜色值：（R，G，B）=（255×（1-gray），255×（1-gray），255×（1-gray））

3 实验结果

为了对比抽希绘图方式和灰度绘图方式的效果，这里使用某地区叠后地震资料做一个显示对比。当前窗口大小为500*400，当水平方向显示道数为250道时使用抽希方式绘图和使用灰度方式绘图效果对比如图5。

而当继续扩大水平方向显示范围（缩小视图），窗口水平方向显示500道数据时，使用抽希方式绘图和使用灰度方式绘图效果对比如图6。

通过试验结果对比可以看出，使用灰度方式绘图精度更高，在窗口显示大范围地震道时能准确呈现出地震资料的细节信息，使人眼睛看到的效果和实际资料更加吻合。

4 结论

本文所介绍的灰度图绘制地震波形和变面积剖面方法相对简单，容易实现，同时不需要特殊显卡支持，主要使用CPU进行计算处理。使用这种方法绘制地震剖面可以显示更大范围地震资料，显示精度高，更加美观，为资料处理和解释人员了解资料的整体详细信息提供了更加方便和精确的方法。

参考文献：

[1]陆基盂.地震勘探原理[M].石油工业出版社，1992.

[2]蔡军，龚幸林.面块切片在三维地震资料解释中的应用[J].中国煤田地质，2004，16（2）：40-41.

[3]杨勤勇.全三维地震解释[J].地球物理学进展，1999，14（3）：128-137.

[4]陈茂山，王云高.油气勘探软件用户界面设计方法探索与实践[J].勘探地球物理进展，2004，27（2）：132-138.

[5]王宏琳.地震软件技术-勘探地球物理计算机软件开发[M].北京：石油工业出版社，2005：4.

[6]刘继承，姚亮.基于VC++的SEGY数据格式地震剖面图绘制[J].自动化技术与应用，2012，31（3）：24-27.

[7]周雪芹，刘晓红.利用VC实现真彩色位图到灰度图的转换[J].自动化技术与应用，2005，24（3）：25-27.