基于Matlab和Madagascar的地震三維圖形顯示技術(shù)研究

楊志鵬，陳秀清，陳光容，王登偉，王 松

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

1 引 言

地震數(shù)據(jù)三維圖形顯示技術(shù)是利用三維地震數(shù)據(jù)體顯示、描述和解釋地下地質(zhì)現(xiàn)象和特征的一種圖像顯示手段和工具[1]?，F(xiàn)階段國(guó)內(nèi)外商業(yè)地震數(shù)據(jù)三維圖形軟件主要以工作站Unix操作系統(tǒng)下為主，如LandMark公司的OpenVision和EarthCube、GeoQuest公司的CPS-3和Geoview、Paradigm的VoxelGeo、BGP公司的GeoEast-EasyTrack等三維可視化軟件[2]。雖然這些商業(yè)可視化軟件具有較高的集成度和更強(qiáng)大的顯示功能，但對(duì)于部分地球物理科研領(lǐng)域工作者來(lái)言，在進(jìn)行算法開發(fā)和驗(yàn)證時(shí)，當(dāng)僅需顯示小塊三維地震數(shù)據(jù)模型時(shí)，這些商業(yè)軟件就顯得過(guò)于龐大和復(fù)雜，且如果在當(dāng)前主流的科研計(jì)算平臺(tái)Matlab或Madagascar下進(jìn)行算法開發(fā)，利用商業(yè)軟件進(jìn)行三維顯示還涉及到頻繁的數(shù)據(jù)格式轉(zhuǎn)換問(wèn)題。近年來(lái)，基于Matlab和Madagascar的三維圖形顯示技術(shù)在勘探地震科研領(lǐng)域得到了泛的應(yīng)用和發(fā)展，尤其是在地震數(shù)據(jù)去噪與重建方面[3]，相比于二維圖形顯示[4]，三維顯示結(jié)果能更直觀全面地體現(xiàn)不同方向上數(shù)據(jù)品質(zhì)[5]。本文在前人工作的基礎(chǔ)上[6，7]，詳細(xì)對(duì)比研究了在Matlab和Madagascar軟件下實(shí)現(xiàn)地震三維圖形顯示的技術(shù)方法，并給出了立體圖、切片圖、層位圖、連井剖面圖、展開圖等常見三維顯示方式的核心程序代碼，通過(guò)示例分析證明了本文開發(fā)的地震三維圖形顯示技術(shù)能夠多角度、多方面展示數(shù)據(jù)內(nèi)部信息，具有良好的應(yīng)用前景。

2 Matlab三維圖形顯示技術(shù)

Matlab是MathWorks公司提供的科學(xué)計(jì)算軟件平臺(tái)和編程語(yǔ)言，其在數(shù)據(jù)可視化方面提供了全面和強(qiáng)大的功能，基于Matlab的三維可視化技術(shù)已在煤礦TEM數(shù)據(jù)體顯示、砂巖型鈾礦礦體顯示[8]、三維地質(zhì)建模[9]等領(lǐng)域得到了廣泛應(yīng)用。在勘探地震領(lǐng)域，開源Matlab軟件Seislab3.02地震工具箱[10]雖然提供了繪制三維地震數(shù)據(jù)體的函數(shù)s_volume_browser(seismic)，但其依賴的輸入?yún)?shù)seismic為較復(fù)雜的結(jié)構(gòu)體，在輸入端不易構(gòu)造正確，且在功能上只能顯示立體圖和切片圖，實(shí)際使用效果不佳。因此，本文基于Matlab中多維作圖函數(shù)slice，給出了其在繪制三維地震數(shù)據(jù)立體圖、切片圖、層位圖、連井剖面圖的具體方法和核心程序。

2.1 Matlab三維地震數(shù)據(jù)網(wǎng)格化

Matlab在繪制三維數(shù)據(jù)體時(shí)，首先需得到三維網(wǎng)格坐標(biāo)?？紤]已導(dǎo)入mat格式的三維地震數(shù)據(jù)體Dtime(x,y,z)，其維度大小為x= 1，…，Nsamp，y= 1，…，NInline，z= 1，…，NCrossline。假設(shè)tmin為起始時(shí)間，dt為采樣時(shí)間間隔，T為采樣時(shí)間序列，Inline_first和Inline_end分別為主測(cè)線的起始和結(jié)束序號(hào)，Crossline_first和Crossline_end分別為聯(lián)絡(luò)線的起始和結(jié)束序號(hào)，則可以利用網(wǎng)格化函數(shù)meshgrid求得由三個(gè)向量定義的三維網(wǎng)格坐標(biāo)X，Y，Z，以便用于計(jì)算三變量函數(shù)或繪制三維立體圖，程序如下：

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);

需注意此時(shí)求得的網(wǎng)格坐標(biāo)矩陣X，Y，Z的網(wǎng)格大小為NCrossline×NInline×Nsamp，因此為了將其與原始數(shù)據(jù)體Dtime(x,y,z)的維度對(duì)應(yīng)一致，還應(yīng)使用數(shù)組維度置換函數(shù)permute對(duì)Dtime(x,y,z)的維度進(jìn)行重排，得到與X，Y，Z維度一致的重排數(shù)據(jù)體V。程序如下：

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

2.2 Matlab三維繪圖函數(shù)及程序

slice函數(shù)是Matlab中常用的多維作圖函數(shù)，通常與meshgrid函數(shù)配合使用完成三維圖件的繪制，也是本文繪圖的核心工具之一，其一般用法為[11]：

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

s為返回的圖像句柄；X，Y，Z為三維網(wǎng)格坐標(biāo)，V為三維數(shù)據(jù)體，method為插值方法；xslice、yslice、zslice分別為3個(gè)坐標(biāo)軸的切片值，且可以指定具體形式：

1)標(biāo)量——在指定位置繪制一個(gè)與對(duì)應(yīng)坐標(biāo)軸正交的切片平面;

2)向量——在指定位置繪制多個(gè)與對(duì)應(yīng)坐標(biāo)軸正交的切片平面;

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

4)矩陣——沿曲面而不是平面繪制切片，且此時(shí)xslice、yslice、zslice均必須是相同大小的矩陣。

2.2.1 立體圖和切片圖繪制

本文以某工區(qū)三維疊后實(shí)際地震資料為例，該數(shù)據(jù)體D大小為251×72×45(時(shí)間軸采樣點(diǎn)數(shù)×主測(cè)線采樣點(diǎn)數(shù)×聯(lián)絡(luò)線采樣點(diǎn)數(shù))，采樣間隔為2 ms，通過(guò)繪制不同類型圖件進(jìn)行說(shuō)明。

立體圖和切片圖是當(dāng)前最常用的顯示三維地震數(shù)據(jù)的方式[12]，能夠?qū)θ齻€(gè)平面方向的地震剖面進(jìn)行顯示，直觀且全面地展示構(gòu)造的空間展布形態(tài)。圖1(a)立體圖顯示了t=0 ms，Inline=1，Crossline=1時(shí)的三個(gè)地震剖面，展示了數(shù)據(jù)體各個(gè)維度上最外層構(gòu)造信息。圖1(b)切片圖顯示了t= 250 ms，Inline =[25 65]，Crossline =[20 40]時(shí)的五個(gè)地震剖面，展示了數(shù)據(jù)體各維度上指定位置處的內(nèi)部構(gòu)造信息。程序如下：

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 層位圖和連井剖面圖繪制

層位圖提取在地震解釋工作中有更重要的作用，其反映了特定地震地質(zhì)層時(shí)間域的構(gòu)造面[13]，在Matlab中繪制特定層位面需額外有該層位的走時(shí)數(shù)據(jù)horizonTime，其維度大小為NInline×NCrossline，圖2(a)層位圖顯示了N1層面，Inline = 65，Crossline =40時(shí)的構(gòu)造信息。zh.xd、zh.yd、zh.zd均為繪制曲面切片時(shí)必要的參數(shù)矩陣，一般用slice分兩步進(jìn)行曲面層位圖繪制，程序如下：

圖2 三維地震數(shù)據(jù)層位圖和連井剖面

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

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’);

連井剖面圖是將井的相關(guān)信息在三維空間中顯示，層速度在地震資料地質(zhì)解釋過(guò)程中非常有用，聲波測(cè)井曲線可以得到精細(xì)層速度資料，以約束速度反演結(jié)果[14]。圖2(b)連井剖面圖顯示了在Inline = 15和Crossline = 40這兩個(gè)剖面上的測(cè)井信息，其繪制方法與切片圖類似，設(shè)測(cè)井?dāng)?shù)據(jù)wlogs為結(jié)構(gòu)體數(shù)據(jù)，其包含了鉆井位置坐標(biāo)信息wlogs.inline和wlogs.crossline，以及井曲線數(shù)據(jù)wlogs.data。將實(shí)際測(cè)井曲線值在相應(yīng)位置插入替換掉反演理論值，程序如下：

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

end

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

3 Madagascar三維圖形顯示技術(shù)

Madagascar是由Texas大學(xué)Austin分校的Sergey Fomel等人開發(fā)的可重復(fù)計(jì)算勘探地球物理軟件包，其提供了豐富的三維地震數(shù)據(jù)處理與成像功能模塊，并提供了C、C++、Fortran、Matlab等多個(gè)常用編程語(yǔ)言的接口，應(yīng)用程序涵蓋了數(shù)值分析、通用數(shù)據(jù)分析、地震資料處理、圖形圖像及可視化等方面[15]。本文基于Madagascar中的多維作圖函數(shù)sfgrey3，給出了其在繪制立體圖和展開圖的方法和程序。

3.1 Madagascar數(shù)據(jù)格式及轉(zhuǎn)換方法

Madagascar讀取數(shù)據(jù)為規(guī)則采樣格式(Regularly Sampled Format，RSF)，其由記錄數(shù)據(jù)屬性的文件頭和以二進(jìn)制形式存儲(chǔ)數(shù)值的數(shù)據(jù)體兩部分構(gòu)成，在概念上RSF格式是一個(gè)高維數(shù)據(jù)體[16]。Madagascar可以利用Python腳本文件Sconstruct批量執(zhí)行程序命令，其主要使用4種命令[17]：Fetch命令是從服務(wù)器端下載獲取數(shù)據(jù)；Flow命令是對(duì)輸入數(shù)據(jù)利用各種函數(shù)模塊功能進(jìn)行處理以產(chǎn)生輸出數(shù)據(jù)文件；Plot命令和Result命令均與Flow命令類似，但輸出結(jié)果增加了相應(yīng)的圖件形式?？紤]將數(shù)據(jù)從Matlab的mat格式轉(zhuǎn)換到Madagascar的RSF格式，其做法可分為兩步，首先利用Matlab中fopen函數(shù)和fwrite函數(shù)將數(shù)據(jù)從mat格式轉(zhuǎn)換成二進(jìn)制文件bin格式，然后利用Madagascar中sfbin2rsf函數(shù)將其轉(zhuǎn)換成RSF格式。以上一節(jié)演示數(shù)據(jù)D為例，將D.mat數(shù)據(jù)文件轉(zhuǎn)換成Data.bin文件的程序?yàn)椋?/p>

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

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

并進(jìn)一步在Sconstruct腳本中將其轉(zhuǎn)換成RSF格式，其程序?yàn)椋?/p>

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軟件的格式轉(zhuǎn)換函數(shù)，n1,n2,n3表示數(shù)據(jù)的維度，o1,o2,o3分別表示軸的起始點(diǎn)，d1,d2,d3分別表示軸的采樣間隔。

3.2 Madagascar三維繪圖函數(shù)及程序

圖3(a)立體圖顯示了t=300 ms，Inline=31，Crossline=26時(shí)的三個(gè)地震剖面，圖中的藍(lán)色線條指示了各個(gè)方向上所顯示的剖面序號(hào)。圖3(b)展開圖是將立體圖進(jìn)行展開得到，程序如下：

圖3 三維地震數(shù)據(jù)立體圖和展開圖

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為控制立體圖是否展開的關(guān)鍵參數(shù);frame1，frame2，frame3這三個(gè)參數(shù)控制需顯示的切片剖面;point1，point2這兩個(gè)參數(shù)控制立方體的縱橫比。另外，Madagascar產(chǎn)生的圖件為特殊的vpl圖片格式，可以用sfpen函數(shù)臨時(shí)打開查看，也可以通過(guò)sfvpconvert函數(shù)轉(zhuǎn)化為eps及其他圖片格式，程序如下：

vpconvert file.vpl format = eps;

4 Matlab和Madagascar的應(yīng)用比較

綜合上述分析，Matlab和Madagascar均能有效地對(duì)三維地震數(shù)據(jù)進(jìn)行圖形顯示，但二者又各具特點(diǎn)與差異，總結(jié)如下：

1)在繪制立體圖和展開圖方面，Madagascar可以將三個(gè)維度上任意層面信息投影到最外層進(jìn)行顯示，且有指示線條注明層位序號(hào)，且可以手動(dòng)調(diào)整各個(gè)顯示剖面的縱橫比例，因此相比Matlab顯示功能更豐富。

2)在繪制層位圖和連井剖面圖方面，由于需更多地調(diào)用結(jié)構(gòu)體數(shù)據(jù)進(jìn)行交互顯示，Matlab因其高度集成化的計(jì)算和可視化環(huán)境，以及程序斷點(diǎn)調(diào)試的功能，因此更易于交互實(shí)現(xiàn)，而Madagascar采用類似Python的語(yǔ)法風(fēng)格，程序?qū)懭隨construct腳本編譯執(zhí)行，不利于人機(jī)交互。

3)在圖件格式轉(zhuǎn)換和圖上標(biāo)注等方面，Madagascar圖件有時(shí)需要利用其他圖形處理軟件進(jìn)行后期拼接和裁剪，而Matlab則可以一體化集成實(shí)現(xiàn)。

5 結(jié) 論

本文通過(guò)示例及程序詳細(xì)介紹了基于Matlab和Madagascar兩種科研軟件下實(shí)現(xiàn)地震數(shù)據(jù)體三維圖形顯示的技術(shù)方法和核心實(shí)現(xiàn)步驟，對(duì)比分析討論了兩種軟件繪制三維圖形的特點(diǎn)與差異，具有較大科研實(shí)用性。

致 謝

特別感謝電子科技大學(xué)厙斌博士對(duì)本文的詳細(xì)指導(dǎo)與討論支持。