Surfer網(wǎng)格文件到GOCAD面文件轉(zhuǎn)換方法研究

卞友艷

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

1 引 言

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

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

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

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

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

續(xù)表1

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

length,flag

x1,y1

x2,y2

...

xn,yn

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

3 GOCAD面文件格式

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

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

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

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

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

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

圖2 GOCAD中的面對象及面文件的實例Fig.2 Example of surface object and surface file in GOCAD

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

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

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

圖3 GOCAD面對象的頂點連接方式Fig.3 Vertex connection mode of GOCAD surface objects

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

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

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

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

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

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

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

BinaryReader br = new BinaryReader(fs);//實例化二進制流讀取對象

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

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

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

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

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

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;//計算總網(wǎng)格點數(shù)

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

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

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

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

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

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

//nvert表示邊界點x、y坐標對的數(shù)目

//vertx,verty表示邊界點x、y坐標數(shù)組

//testx,texty表示待判斷的x、y坐標

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ù)索引計算

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

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

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)格點索引

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

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

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)格點坐標及數(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);

//寫入三角面頂點索引

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

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

}

5 應(yīng)用實例

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

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

選取測線上部分測點測量其大地坐標，采用本文介紹的文件格式轉(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中的顯示更加立體形象，可以直觀地分辨測線的位置走向、地形的起伏情況等，證明本文所介紹的方法是可靠有效的。

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

6 結(jié) 論

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

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