基于Matlab和Madagascar的地震三维图形显示技术研究

2020-09-04杨志鹏陈秀清陈光容王登伟

工程地球物理学报 2020年4期

杨志鹏，陈秀清，陈光容，王登伟，王松

(四川省地震局西昌地震中心站，四川西昌 615022)

1 引言

地震数据三维图形显示技术是利用三维地震数据体显示、描述和解释地下地质现象和特征的一种图像显示手段和工具[1]。现阶段国内外商业地震数据三维图形软件主要以工作站Unix操作系统下为主，如LandMark公司的OpenVision和EarthCube、GeoQuest公司的CPS-3和Geoview、Paradigm的VoxelGeo、BGP公司的GeoEast-EasyTrack等三维可视化软件[2]。虽然这些商业可视化软件具有较高的集成度和更强大的显示功能，但对于部分地球物理科研领域工作者来言，在进行算法开发和验证时，当仅需显示小块三维地震数据模型时，这些商业软件就显得过于庞大和复杂，且如果在当前主流的科研计算平台Matlab或Madagascar下进行算法开发，利用商业软件进行三维显示还涉及到频繁的数据格式转换问题。近年来，基于Matlab和Madagascar的三维图形显示技术在勘探地震科研领域得到了泛的应用和发展，尤其是在地震数据去噪与重建方面[3]，相比于二维图形显示[4]，三维显示结果能更直观全面地体现不同方向上数据品质[5]。本文在前人工作的基础上[6，7]，详细对比研究了在Matlab和Madagascar软件下实现地震三维图形显示的技术方法，并给出了立体图、切片图、层位图、连井剖面图、展开图等常见三维显示方式的核心程序代码，通过示例分析证明了本文开发的地震三维图形显示技术能够多角度、多方面展示数据内部信息，具有良好的应用前景。

2 Matlab三维图形显示技术

Matlab是MathWorks公司提供的科学计算软件平台和编程语言，其在数据可视化方面提供了全面和强大的功能，基于Matlab的三维可视化技术已在煤矿TEM数据体显示、砂岩型铀矿矿体显示[8]、三维地质建模[9]等领域得到了广泛应用。在勘探地震领域，开源Matlab软件Seislab3.02地震工具箱[10]虽然提供了绘制三维地震数据体的函数s_volume_browser(seismic)，但其依赖的输入参数seismic为较复杂的结构体，在输入端不易构造正确，且在功能上只能显示立体图和切片图，实际使用效果不佳。因此，本文基于Matlab中多维作图函数slice，给出了其在绘制三维地震数据立体图、切片图、层位图、连井剖面图的具体方法和核心程序。

2.1 Matlab三维地震数据网格化

Matlab在绘制三维数据体时，首先需得到三维网格坐标。考虑已导入mat格式的三维地震数据体Dtime(x,y,z)，其维度大小为x= 1，…，Nsamp，y= 1，…，NInline，z= 1，…，NCrossline。假设tmin为起始时间，dt为采样时间间隔，T为采样时间序列，Inline_first和Inline_end分别为主测线的起始和结束序号，Crossline_first和Crossline_end分别为联络线的起始和结束序号，则可以利用网格化函数meshgrid求得由三个向量定义的三维网格坐标X，Y，Z，以便用于计算三变量函数或绘制三维立体图，程序如下：

load D;

[NsampNInlineNCrossline]= size(D);

tmin= 0;

dt= 2;

T=tmin: dt:tmin+(Nsamp-1)*dt;

Inline_first = 1;

Crossline_first = 1;

Inline_end = Inline_first +NInline-1;

Crossline_end = Crossline_first +NCrossline-1;

[X,Y,Z]= meshgrid(Inline_first : Inline_end, Crossline_first : Crossline_end,T);

需注意此时求得的网格坐标矩阵X，Y，Z的网格大小为NCrossline×NInline×Nsamp，因此为了将其与原始数据体Dtime(x,y,z)的维度对应一致，还应使用数组维度置换函数permute对Dtime(x,y,z)的维度进行重排，得到与X，Y，Z维度一致的重排数据体V。程序如下：

V= permute(D,[3 2 1]);

2.2 Matlab三维绘图函数及程序

slice函数是Matlab中常用的多维作图函数，通常与meshgrid函数配合使用完成三维图件的绘制，也是本文绘图的核心工具之一，其一般用法为[11]：

s= slice(X,Y,Z,V, xslice, yslice, zslice, method);

s为返回的图像句柄；X，Y，Z为三维网格坐标，V为三维数据体，method为插值方法；xslice、yslice、zslice分别为3个坐标轴的切片值，且可以指定具体形式：

1)标量——在指定位置绘制一个与对应坐标轴正交的切片平面;

2)向量——在指定位置绘制多个与对应坐标轴正交的切片平面;

3)[]——不绘制任何切片平面;

4)矩阵——沿曲面而不是平面绘制切片，且此时xslice、yslice、zslice均必须是相同大小的矩阵。

2.2.1 立体图和切片图绘制

本文以某工区三维叠后实际地震资料为例，该数据体D大小为251×72×45(时间轴采样点数×主测线采样点数×联络线采样点数)，采样间隔为2 ms，通过绘制不同类型图件进行说明。

立体图和切片图是当前最常用的显示三维地震数据的方式[12]，能够对三个平面方向的地震剖面进行显示，直观且全面地展示构造的空间展布形态。图1(a)立体图显示了t=0 ms，Inline=1，Crossline=1时的三个地震剖面，展示了数据体各个维度上最外层构造信息。图1(b)切片图显示了t= 250 ms，Inline =[25 65]，Crossline =[20 40]时的五个地震剖面，展示了数据体各维度上指定位置处的内部构造信息。程序如下：

Fig1a = slice(X,Y,Z,V, 1, 1, 0, ‘cubic’);

Fig1b= slice(X,Y,Z,V,[25 65],[20 40], 250, ‘cubic’);

shading interp;

set(gca, ‘zdir’, ‘reverse’);

2.2.2 层位图和连井剖面图绘制

层位图提取在地震解释工作中有更重要的作用，其反映了特定地震地质层时间域的构造面[13]，在Matlab中绘制特定层位面需额外有该层位的走时数据horizonTime，其维度大小为NInline×NCrossline，图2(a)层位图显示了N1层面，Inline = 65，Crossline =40时的构造信息。zh.xd、zh.yd、zh.zd均为绘制曲面切片时必要的参数矩阵，一般用slice分两步进行曲面层位图绘制，程序如下：

图2 三维地震数据层位图和连井剖面

load horizonTime;

zh.xd = zeros(NInline,NCrossline);

zh.yd = zeros(NInline,NCrossline);

zh.zd = zeros(NInline,NCrossline);

fori= Inline_first : Inline_end

forj=Crossline_first : Crossline_end

zh.xd(i,j)=i;

zh.yd(i,j)=j;

zh.zd(i,j)= horizonTime(i,j)

end

Fig2a = slice(X,Y,Z,V, 65, 40,[], ‘cubic’);

hold on;

Fig2a = slice(X,Y,Z,V, zh.xd, zh.yd, zh.zd, ‘cubic’);

连井剖面图是将井的相关信息在三维空间中显示，层速度在地震资料地质解释过程中非常有用，声波测井曲线可以得到精细层速度资料，以约束速度反演结果[14]。图2(b)连井剖面图显示了在Inline = 15和Crossline = 40这两个剖面上的测井信息，其绘制方法与切片图类似，设测井数据wlogs为结构体数据，其包含了钻井位置坐标信息wlogs.inline和wlogs.crossline，以及井曲线数据wlogs.data。将实际测井曲线值在相应位置插入替换掉反演理论值，程序如下：

load wlogs;

ix =wlogs.inline-Inline_first + 1;

iy = wlogs.crossline-Crossline_first + 1;

form= ix-1 : ix

forn= iy-1 : iy

V(n,m, :)= wlogs.data(:);

end

Fig2b= slice(X,Y,Z,V, 15, 40, 400, ‘cubic’);

3 Madagascar三维图形显示技术

Madagascar是由Texas大学Austin分校的Sergey Fomel等人开发的可重复计算勘探地球物理软件包，其提供了丰富的三维地震数据处理与成像功能模块，并提供了C、C++、Fortran、Matlab等多个常用编程语言的接口，应用程序涵盖了数值分析、通用数据分析、地震资料处理、图形图像及可视化等方面[15]。本文基于Madagascar中的多维作图函数sfgrey3，给出了其在绘制立体图和展开图的方法和程序。

3.1 Madagascar数据格式及转换方法

Madagascar读取数据为规则采样格式(Regularly Sampled Format，RSF)，其由记录数据属性的文件头和以二进制形式存储数值的数据体两部分构成，在概念上RSF格式是一个高维数据体[16]。Madagascar可以利用Python脚本文件Sconstruct批量执行程序命令，其主要使用4种命令[17]：Fetch命令是从服务器端下载获取数据；Flow命令是对输入数据利用各种函数模块功能进行处理以产生输出数据文件；Plot命令和Result命令均与Flow命令类似，但输出结果增加了相应的图件形式。考虑将数据从Matlab的mat格式转换到Madagascar的RSF格式，其做法可分为两步，首先利用Matlab中fopen函数和fwrite函数将数据从mat格式转换成二进制文件bin格式，然后利用Madagascar中sfbin2rsf函数将其转换成RSF格式。以上一节演示数据D为例，将D.mat数据文件转换成Data.bin文件的程序为：

fid = fopen(‘Data.bin’ , ‘w’);

fwrite(fid, reshape(D, 251, 72*45), ’float’);

并进一步在Sconstruct脚本中将其转换成RSF格式，其程序为：

Flow('Data_RSF' , '../Data.bin' , 'bin2rsf bfile=${SOURCES[0]} n1=251 n2=3240 | put n2=72 n3=45 d1=2 d2=1 d3=1 o1=0 o2=1 o3=1');

其中Data_RSF为输出的RSF格式文件，Data.bin为输入的bin格式文件，bin2rsf为Madagascar软件的格式转换函数，n1,n2,n3表示数据的维度，o1,o2,o3分别表示轴的起始点，d1,d2,d3分别表示轴的采样间隔。

3.2 Madagascar三维绘图函数及程序

图3(a)立体图显示了t=300 ms，Inline=31，Crossline=26时的三个地震剖面，图中的蓝色线条指示了各个方向上所显示的剖面序号。图3(b)展开图是将立体图进行展开得到，程序如下：

图3 三维地震数据立体图和展开图

def GreyFig3a(data, other):

Result(data, ''' byte gainpanel=all maxval=1000 minval=-1000 bar=bar.rsf | grey3 wanttitle=yflat=yframe1=150 frame2=30 frame3=25 point1=0.5 point2=0.5 label1="Time" label2=Inline unit2= label3=Crossline unit3= unit1=ms title= screenratio=1 bar=bar.rsf scalebar=n barlabel=Amplitude color=i %s ''' %other)

def GreyFig3b(data, other):

Result(data, ''' byte gainpanel=all maxval=1000 minval=-1000 bar=bar.rsf | grey3 wanttitle=yflat=nframe1=150 frame2=30 frame3=25 point1=0.5 point2=0.5 label1="Time" label2=Inline unit2= label3=Crossline unit3= unit1=ms title= screenratio=1 bar=bar.rsf scalebar=y barlabel=Amplitude color=i %s ''' %other)

GreyFig3a(‘Data_RSF’, ‘title = ’)

GreyFig3b(‘Data_RSF’, ‘title = ’)

其中flat=y/n为控制立体图是否展开的关键参数;frame1，frame2，frame3这三个参数控制需显示的切片剖面;point1，point2这两个参数控制立方体的纵横比。另外，Madagascar产生的图件为特殊的vpl图片格式，可以用sfpen函数临时打开查看，也可以通过sfvpconvert函数转化为eps及其他图片格式，程序如下：

vpconvert file.vpl format = eps;