卞友艷

(中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司，湖北 武漢 430063)

1 引 言

近年來(lái)，隨著工程建設(shè)的發(fā)展，以城市地下空間開發(fā)為代表的工程三維勘探技術(shù)得到了廣泛的應(yīng)用[1]，基于其成果的三維地質(zhì)建模工作越來(lái)越重要。GOCAD(Geological Object Computer Aided Design)是由法國(guó)Nancy大學(xué)開發(fā)的一款地質(zhì)對(duì)象計(jì)算機(jī)輔助設(shè)計(jì)軟件，具有強(qiáng)大的三維建模、可視化、地質(zhì)解譯和分析功能，是目前國(guó)際上公認(rèn)的主流三維地質(zhì)建模軟件，在地質(zhì)工程、地球物理勘探、礦業(yè)開發(fā)、水利工程中均有廣泛的應(yīng)用[2,3]。Surfer是地質(zhì)勘察領(lǐng)域常用的繪圖軟件，具有強(qiáng)大的數(shù)據(jù)插值與繪圖功能，地球物理勘探數(shù)據(jù)通常會(huì)利用該軟件進(jìn)行網(wǎng)格化，成果數(shù)據(jù)一般以GRD網(wǎng)格文件的形式保存[4,5]。

在三維地質(zhì)建模中，地球物理勘探成果是重要的建模依據(jù)，而Surfer網(wǎng)格文件作為其常用的數(shù)據(jù)保存格式，無(wú)法直接導(dǎo)入到GOCAD中。雖可在提取網(wǎng)格文件數(shù)據(jù)點(diǎn)的基礎(chǔ)上，再以點(diǎn)數(shù)據(jù)的形式進(jìn)行導(dǎo)入，但其過(guò)程比較繁瑣，且還需要通過(guò)GOCAD手動(dòng)生成網(wǎng)格面[6,7]，整個(gè)過(guò)程需要人工不停地切換軟件，效率很低。為解決上述問(wèn)題，本文研究了Surfer網(wǎng)格文件和GOCAD面文件的格式，分析了其數(shù)據(jù)存儲(chǔ)的特點(diǎn)，并通過(guò)C#語(yǔ)言編程實(shí)現(xiàn)了兩者之間格式的快速轉(zhuǎn)換，以便于在GOCAD中快速導(dǎo)入地球物理勘探成果數(shù)據(jù)，提高三維地質(zhì)建模的效率。

2 Surfer網(wǎng)格文件格式

Surfer是一款基于網(wǎng)格化成圖的繪圖軟件，它可以將不規(guī)則間隔的數(shù)據(jù)插入到規(guī)則間隔的網(wǎng)格中，地球物理勘探中測(cè)量或反演的散點(diǎn)數(shù)據(jù)通常會(huì)利用Surfer生成規(guī)則的網(wǎng)格化文件，再進(jìn)行成圖顯示[8-10]。Surfer經(jīng)過(guò)多年的迭代更新，形成了3種網(wǎng)格文件格式，分別為：Surfer 6 Text Grid、Surfer 6 Binary Grid、Surfer 7 Binary Grid。雖然三種數(shù)據(jù)格式分別采用了ASCII文本格式和二進(jìn)制格式，但所包含的數(shù)據(jù)信息大同小異，以最新的Surfer 7 Binary Grid(以下簡(jiǎn)稱“Surfer7網(wǎng)格文件”)為例，其具體文件格式見表1。

表1 Surfer7網(wǎng)格文件數(shù)據(jù)格式(摘自Surfer幫助文檔)

續(xù)表1

此外，對(duì)于無(wú)意義的網(wǎng)格區(qū)域，通常會(huì)使用白化文件(*.bln)“剔除”該區(qū)域，其文件格式為文本文檔，內(nèi)容如下[11]：

length,flag

x1,y1

x2,y2

...

xn,yn

其中，length表示以下幾行的x、x坐標(biāo)對(duì)的數(shù)量，這些坐標(biāo)點(diǎn)依次連接形成一個(gè)閉合區(qū)域，flag指示將該區(qū)域以內(nèi)的部分還是該區(qū)域以外的部分“白化”，前者值為1，后者為0，一般白化內(nèi)部區(qū)域，即flag為1。

3 GOCAD面文件格式

GOCAD中定義了一些基本的幾何對(duì)象，如點(diǎn)(PointsSet)、線(Curve)、面(Surface)、地層網(wǎng)格(SGrid)等，通過(guò)它們來(lái)構(gòu)建復(fù)雜的三維地質(zhì)模型。這些對(duì)象包含了必要的信息，如空間坐標(biāo)、表示物理特性的數(shù)值數(shù)據(jù)、測(cè)量單位的首選項(xiàng)、可視化顏色等，它們可以單獨(dú)導(dǎo)出保存為文件，也可以從文件導(dǎo)入到項(xiàng)目工程中[12-16]。

GOCAD對(duì)象以ASCII文本文件的格式保存，文件可分為三部分：文件頭、數(shù)據(jù)體和結(jié)束標(biāo)記。文件頭定義了對(duì)象的類型、對(duì)象的名稱以及其他可視化選項(xiàng)等；數(shù)據(jù)體包含對(duì)象的幾何和屬性信息，例如節(jié)點(diǎn)坐標(biāo)、節(jié)點(diǎn)的拓?fù)潢P(guān)系、節(jié)點(diǎn)的屬性值、坐標(biāo)系統(tǒng)、坐標(biāo)單位等；結(jié)束標(biāo)記表示當(dāng)前描述對(duì)象的結(jié)束，利用結(jié)束標(biāo)記，多個(gè)對(duì)象也可以存儲(chǔ)在一個(gè)文件中。

基于Surfer中網(wǎng)格文件的呈現(xiàn)特點(diǎn)，類似于GOCAD中帶屬性的面對(duì)象，可將Surfer網(wǎng)格文件轉(zhuǎn)換為GOCAD面對(duì)象的文件，其格式見圖1。

圖1 GOCAD面對(duì)象(帶屬性)基本文件格式Fig.1 Basic file format of GOCAD surface object (with attributes)

其中，TSurf為GOCAD面對(duì)象的文件標(biāo)識(shí)，后面的數(shù)字表示版本號(hào)；name后面表示該對(duì)象的名稱；PROPERTIES后面為屬性值的名稱，PVRTX表示帶屬性的頂點(diǎn)，其后分別為頂點(diǎn)的索引、頂點(diǎn)的空間坐標(biāo)(xyz)以及頂點(diǎn)的屬性值。

在GOCAD中，面對(duì)象可以不帶屬性，也可以帶一個(gè)或多個(gè)屬性，從Surfer網(wǎng)格文件轉(zhuǎn)換的面對(duì)象，只需一個(gè)屬性即可。GOCAD的面對(duì)象由多個(gè)三角面構(gòu)成，每個(gè)三角面包含三個(gè)頂點(diǎn)，每個(gè)頂點(diǎn)可以包含在多個(gè)三角面中，TRGL后面表示的則是構(gòu)成三角面的頂點(diǎn)的索引。以上為面對(duì)象文件最基本的元素，通過(guò)添加其他的關(guān)鍵字以及數(shù)值可以賦予面對(duì)象更多的信息，例如面對(duì)象的顏色、空值、坐標(biāo)系統(tǒng)、坐標(biāo)單位等，圖2為一個(gè)簡(jiǎn)單的GOCAD面文件的實(shí)例。

圖2 GOCAD中的面對(duì)象及面文件的實(shí)例Fig.2 Example of surface object and surface file in GOCAD

4 轉(zhuǎn)換方法的實(shí)現(xiàn)

4.1 文件轉(zhuǎn)換方法

從Surfer網(wǎng)格文件轉(zhuǎn)為GOCAD面文件，首先需要讀取Surfer網(wǎng)格文件里面的網(wǎng)格點(diǎn)數(shù)據(jù)。結(jié)合x、y方向坐標(biāo)最小值以及網(wǎng)格間距，可以計(jì)算出每一個(gè)網(wǎng)格點(diǎn)的坐標(biāo)值，但該坐標(biāo)值是二維坐標(biāo)，無(wú)法表示數(shù)據(jù)點(diǎn)真實(shí)的空間位置，因此需要進(jìn)行坐標(biāo)的轉(zhuǎn)換匹配。地球物理勘探可分為剖面勘探和面積性勘探，生成網(wǎng)格文件后，其數(shù)據(jù)網(wǎng)格的x、y方向代表的含義不同，坐標(biāo)匹配與轉(zhuǎn)換的方法也不相同。對(duì)于剖面勘探數(shù)據(jù)，只需測(cè)量頂部若干測(cè)點(diǎn)的大地坐標(biāo)，進(jìn)行簡(jiǎn)單的線性插值計(jì)算，即可獲取每個(gè)網(wǎng)格點(diǎn)的大地坐標(biāo)，而對(duì)于面積性勘探數(shù)據(jù)，一般會(huì)實(shí)測(cè)每個(gè)網(wǎng)格點(diǎn)的大地坐標(biāo)或進(jìn)行三維空間插值。對(duì)于坐標(biāo)轉(zhuǎn)換匹配方法已經(jīng)比較成熟，本文不做詳細(xì)介紹。

圖3 GOCAD面對(duì)象的頂點(diǎn)連接方式Fig.3 Vertex connection mode of GOCAD surface objects

在Surfer網(wǎng)格文件中，網(wǎng)格點(diǎn)的數(shù)據(jù)值是按照先從x方向坐標(biāo)由小到大，再?gòu)膟方向坐標(biāo)由小到大的順序排列的?！鞍谆眳^(qū)域的設(shè)置對(duì)圖像顯示效果影響很大，Surfer既可以利用網(wǎng)格點(diǎn)數(shù)據(jù)中的空值(Blank Value)加以區(qū)分，也可以利用另外的白化文件來(lái)設(shè)定，一般采用第二種方式進(jìn)行“白化”。在GOCAD中也有類似“白化”的概念，它是通過(guò)直接設(shè)置空值(No Data Value)來(lái)表示空白區(qū)域。因此，在進(jìn)行文件轉(zhuǎn)換時(shí)，如果存在空白區(qū)域，可以同時(shí)讀取白化文件，獲取白化區(qū)域的邊界點(diǎn)，根據(jù)邊界點(diǎn)形成的閉合多邊形來(lái)判斷每個(gè)網(wǎng)格點(diǎn)是否在白化區(qū)域內(nèi)，若在白化區(qū)域內(nèi)，則將數(shù)據(jù)值設(shè)為空值。判斷某點(diǎn)是否落在區(qū)域內(nèi)，可以通過(guò)面積和判別法、夾角和判別法、引射線法等[17,18]，本文采用射線法進(jìn)行判斷[19-21]，以上過(guò)程可在二維坐標(biāo)下完成。

根據(jù)GOCAD面文件的特點(diǎn)，TRGL后面記錄了組成面的若干個(gè)三角面的頂點(diǎn)索引，在讀取的Surfer網(wǎng)格數(shù)據(jù)中，也需要將每個(gè)網(wǎng)格點(diǎn)進(jìn)行類似的組合連接。由于Surfer網(wǎng)格點(diǎn)是規(guī)則分布的，且是順序保存的，因此可以依次連接空間上相鄰的3個(gè)網(wǎng)格點(diǎn)組成三角面，同時(shí)也可以很容易獲取組成三角面的每個(gè)網(wǎng)格點(diǎn)的索引。如圖3所示的連接組合方式，共可組成(nRow-1)*(nCol-1)*2個(gè)三角面。將匹配轉(zhuǎn)換后的網(wǎng)格點(diǎn)坐標(biāo)及對(duì)應(yīng)的網(wǎng)格點(diǎn)數(shù)據(jù)值按讀取的順序依次寫入新的GOCAD面文件中，同時(shí)將上述計(jì)算的所有三角面的頂點(diǎn)的索引也寫入文件中，即可以得到Surfer網(wǎng)格文件轉(zhuǎn)換后的GOCAD面文件。

4.2 轉(zhuǎn)換程序關(guān)鍵代碼

正如前文所述，轉(zhuǎn)換文件時(shí)關(guān)鍵步驟是讀取Surfer網(wǎng)格文件里面的網(wǎng)格點(diǎn)數(shù)據(jù)、設(shè)置“白化”區(qū)域、面文件數(shù)據(jù)索引計(jì)算和生成GOCAD面文件，詳細(xì)代碼如下。

4.2.1 以Surfer7網(wǎng)格文件為例，讀取Surfer網(wǎng)格文件里面的網(wǎng)格點(diǎn)數(shù)據(jù)

Grid grd = new Grid();//實(shí)例化網(wǎng)格類，存儲(chǔ)網(wǎng)格數(shù)據(jù)

FileStream fs = new FileStream(filepath, FileMode.Open);//實(shí)例化文件流對(duì)象

BinaryReader br = new BinaryReader(fs);//實(shí)例化二進(jìn)制流讀取對(duì)象

br.ReadBytes(20);//...

grd.nRow = br.ReadInt32();//讀取網(wǎng)格行數(shù)

grd.nCol = br.ReadInt32();//讀取網(wǎng)格列數(shù)

grd.xLL = br.ReadDouble();//讀取x方向坐標(biāo)最小值

grd.yLL = br.ReadDouble();//讀取y方向坐標(biāo)最小值

grd.xSize = br.ReadDouble();//讀取x方向網(wǎng)格間距

grd.ySize = br.ReadDouble();//讀取y方向網(wǎng)格間距

grd.vMin = br.ReadDouble();//讀取數(shù)據(jù)體的最小值

grd.vMax = br.ReadDouble();//讀取數(shù)據(jù)體的最大值

grd.Rotation = br.ReadDouble();//讀取網(wǎng)格選裝角度

grd.BlankValue = br.ReadDouble();//讀取空值的代替值

br.ReadBytes(8);//...

int N = grd.nRow * grd.nCol;//計(jì)算總網(wǎng)格點(diǎn)數(shù)

grd.Data = new double[N];//實(shí)例化數(shù)據(jù)數(shù)組

for (int i = 0; i < N; i++)

grd.Data[i] = br.ReadDouble();//讀取數(shù)據(jù)體

br.Close();//關(guān)閉二進(jìn)制流

fs.Close();//關(guān)閉文件流

4.2.2 設(shè)置“白化”區(qū)域

//nvert表示邊界點(diǎn)x、y坐標(biāo)對(duì)的數(shù)目

//vertx,verty表示邊界點(diǎn)x、y坐標(biāo)數(shù)組

//testx,texty表示待判斷的x、y坐標(biāo)

int i, j;

int c = 0;//表示是否在多邊形內(nèi)，是為1，否為0

for (i=0, j=nvert-1; i

{

if (((verty[i] > testy) != (verty[j] > testy)) && (testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))

c = !c;

}

4.2.3 面文件數(shù)據(jù)索引計(jì)算

Trgl[] trgl = new Trgl[(row -1) * (col - 1) * 2];

//實(shí)例化Trgl數(shù)組，每個(gè)Trgl包含3個(gè)索引值

int idx = 0;

for (int i = 0; i < row - 1; i++)

{

for (int j = 0; j < col - 1; j++)

{

trgl[idx] = new Trgl(i * Col + j + 1, i * Col + j + 2, (i + 1) * Col + j + 1);

//添加下層三角面的網(wǎng)格點(diǎn)索引

trgl[idx + 1] = new Trgl((i + 1) * Col + j + 2, i * Col + j + 2, (i + 1) * Col + j + 1);

//添加上層三角面的網(wǎng)格點(diǎn)索引

idx += 2;

}

}

4.2.4 生成GOCAD面文件

using (StreamWriter sw = File.CreateText(filepath))//創(chuàng)建*.ts文件

{

sw.Write(String.Format(@"

GOCAD TSurf 1

HEADER {{

name: {0}

}}

PROPERTIES {1}", "surfer", "prop"));//寫入文件頭信息

for (int i = 0; i < grd.nRow; i++)

{

for (int j = 0; j < grd.nCol; j++)

{

sw.WriteLine(string.Format("PVRTX {0} {1} {2} {3} {4} ", i * grd.nCol + j + 1,grd.xLL + j * grd.xSize, 0, grd.yLL + i * grd.ySize, grd.Data[i * grd.nCol + j]));

//寫入網(wǎng)格點(diǎn)坐標(biāo)及數(shù)據(jù)值

}

}

for (int i = 0; i < (grd.nRow - 1) * (grd.nCol - 1) * 2; i++)

sw.WriteLine("TRGL {0} {1} {2} ", trgl[i].p1, trgl[i].p2, trgl[i].p3);

//寫入三角面頂點(diǎn)索引

sw.WriteLine("END");//寫入結(jié)束標(biāo)記

sw.Close();//關(guān)閉文件

}

5 應(yīng)用實(shí)例

以某隧道的微動(dòng)勘探結(jié)果為例，生成Surfer網(wǎng)格文件，網(wǎng)格大小126×281，橫縱網(wǎng)格間距分別為5.0 m、2.0 m，最小坐標(biāo)為(0，137.793)，由于地形起伏較大，使用了白化文件來(lái)進(jìn)行地表以上區(qū)域的去除，如圖4所示。

圖4 某隧道的微動(dòng)勘探成果Fig.4 The microtremor exploration results of a tunnel

選取測(cè)線上部分測(cè)點(diǎn)測(cè)量其大地坐標(biāo)，采用本文介紹的文件格式轉(zhuǎn)換方法將Surfer網(wǎng)格文件轉(zhuǎn)換為GOCAD面文件，生成新的.ts文件，文件內(nèi)容如下：

GOCAD TSurf 1

HEADER {

name: Tunnel

}

PROPERTIES Velocity

NO_DATA_VALUES -99999

PVRTX 1 606889.8540 3436626.8540 137.793-99999

PVRTX 2 606893.7810 3436630.7805 137.793 -99999

PVRTX 3 606897.7080 3436634.7070 137.793 -99999

...

PVRTX 28202 606889.8540 3436626.8540 337.793 210.541

PVRTX 28203 606893.7810 3436630.7805 337.793 210.150

PVRTX 20204 606897.7080 3436634.7070 337.793 220.009

...

TRGL 1 2 282

TRGL 283 2 282

TRGL 2 3 283

...

END

將生成的GOCAD面文件導(dǎo)入到GOCAD軟件中，結(jié)果如圖5所示，可以看到，轉(zhuǎn)換后的結(jié)果與在Surfer中的顯示基本一致，而在GOCAD中的顯示更加立體形象，可以直觀地分辨測(cè)線的位置走向、地形的起伏情況等，證明本文所介紹的方法是可靠有效的。

圖5 轉(zhuǎn)換后的面文件在GOCAD中的顯示Fig.5 Display of the converted surface file in GOCAD

6 結(jié) 論

本文分析研究了Surfer網(wǎng)格文件與GOCAD面文件的格式，采用C#編程語(yǔ)言實(shí)現(xiàn)了Surfer網(wǎng)格文件到GOCAD面文件的轉(zhuǎn)換。在轉(zhuǎn)換的過(guò)程中，首先需要獲取Surfer網(wǎng)格文件的網(wǎng)格坐標(biāo)及數(shù)據(jù)值，再進(jìn)行坐標(biāo)轉(zhuǎn)換，將二維坐標(biāo)轉(zhuǎn)換為三維坐標(biāo)；然后根據(jù)白化文件，將指定區(qū)域的數(shù)據(jù)值設(shè)為空值，再將各網(wǎng)格點(diǎn)按順序依次連接成三角面，獲取每個(gè)三角面的頂點(diǎn)索引；最后將網(wǎng)格點(diǎn)坐標(biāo)及數(shù)值、三角面頂點(diǎn)索引等寫入新的GOCAD面文件中，完成轉(zhuǎn)換。

本文的研究實(shí)現(xiàn)了地球物理勘探成果文件的快速導(dǎo)入，有助于提高利用GOCAD進(jìn)行三維地質(zhì)建模的效率，具有較高的推廣和應(yīng)用價(jià)值。