譚正華 ，王李管，熊書敏，劉任任

(1.湘潭大學(xué) 湖南省高等學(xué)校智能制造重點實驗室，湖南 湘潭，411105；2.中南大學(xué) 數(shù)字礦山研究中心，湖南 長沙，410083)

地質(zhì)體工程剖面圖通常是沿勘探線方向切繪，反映出該剖面上的礦(巖)體界線、地質(zhì)體巖性，工程地質(zhì)特征、水文特征、斷層位置和構(gòu)造形態(tài)等，它是分析地質(zhì)構(gòu)造、儲量估算和指導(dǎo)采礦設(shè)計的重要圖件之一[1]，是工程師、設(shè)計管理人員等交流的共同語言。地質(zhì)剖面圖的生成要依賴于特定的地質(zhì)體模型，由于地質(zhì)體數(shù)據(jù)量大、結(jié)構(gòu)復(fù)雜(多含夾石、斷層)，且當(dāng)?shù)V體內(nèi)含有孔隙或巖漿侵入另一種礦體，切割后形成的地質(zhì)區(qū)域通常含有孔洞或孤島，不同區(qū)域的部分邊界相互粘連。同時，為了表示不同的地質(zhì)構(gòu)造，需要對不同地質(zhì)體內(nèi)部區(qū)域進(jìn)行識別，并進(jìn)行巖性圖案的填充。國外學(xué)者較早地提出了一些三維實體切割、開挖算法，主要算法有Shostko等[2]提出并實現(xiàn)的算法、gOcad項目組實現(xiàn)的切割算法[3]和三角網(wǎng)切割算法TRICUT[4]，主要思想是將復(fù)雜地質(zhì)模型的剖面生成問題歸結(jié)為空間三角網(wǎng)的相互切割問題，基本操作是三角形與三角形求交和平面與三角形求交，存在大量的幾何運算，時間復(fù)雜度較高。另外，國內(nèi)學(xué)者在實體切割和地質(zhì)體剖面圖生成問題上提出了一些改進(jìn)算法：陳國良等[5]采用一種基于OBB樹的三角網(wǎng)切割算法來實現(xiàn)空間兩個三角網(wǎng)格的切割；陳建宏等[6]提出了基于鉆孔數(shù)據(jù)的剖面圖生成算法；朱瑩等[7?8]提出采用GIS組件技術(shù)實現(xiàn)地質(zhì)體剖面圖的生成等。本文作者提出一種新的基于平面隱函數(shù)的實體切割思想和區(qū)域自動識別技術(shù)的地質(zhì)體剖面圖自動生成方法，將三角形和切割面之間的幾何運算轉(zhuǎn)換為數(shù)學(xué)運算，同時通過樹結(jié)構(gòu)形象表達(dá)地質(zhì)體輪廓線形成的復(fù)雜區(qū)域(內(nèi)含孔、島)，最后通過對復(fù)雜區(qū)域的巖性圖案填充，從而實現(xiàn)地質(zhì)體剖面圖自動生成。實例驗證了該算法有效性和可行性。

1 基本原理及流程

由于地質(zhì)體的形態(tài)復(fù)雜，直接基于原始地質(zhì)模型的剖切較為復(fù)雜。因此，本文作者提出基于隱函數(shù)平面切割算法，加快幾何運算，生成離散的交線；采用基于KD樹的空間索引方法，確定交線的鄰接關(guān)系，生成地質(zhì)體輪廓線；采用樹結(jié)構(gòu)形象表達(dá)地質(zhì)體輪廓線形成的復(fù)雜區(qū)域（簡稱“區(qū)域樹”）；為了表達(dá)地質(zhì)體巖性信息，對剖面圖地質(zhì)體區(qū)域按照標(biāo)準(zhǔn)巖性圖例進(jìn)行圖案填充，完成剖面圖的自動繪制。算法流程見圖1。

圖1 地質(zhì)體剖切算法流程Fig.1 Process of 3D entities cutting

2 三維地質(zhì)結(jié)構(gòu)模型的組織

地質(zhì)體剖切算法依賴特定的數(shù)據(jù)模型。表達(dá)地質(zhì)體基本數(shù)據(jù)模型有表征實體外表輪廓的表面模型、基于體元的塊段模型以及混合模型[9]。其中，表面模型側(cè)重于三維空間實體的表面表示，能精確表達(dá)地質(zhì)體的形態(tài)。這些被模擬的實體可以是封閉的，如礦體、硐室和采場等；也可以是非封閉的，如地層界面、斷層和地表等。

數(shù)字采礦系統(tǒng)多采用不規(guī)則三角格網(wǎng)(TIN)來模擬各種復(fù)雜地質(zhì)體邊界及人工工程的表面，本文研究的切剖對象為TIN表達(dá)的地質(zhì)體模型，如圖2所示。

圖2 基于TIN表達(dá)的地質(zhì)體模型Fig.2 Geological model expressed by TIN

3 算法的關(guān)鍵過程

3.1 隱函數(shù)的平面切割算法

3.1.1 隱函數(shù)的性質(zhì)

用平面切割多面體形成一系列的輪廓線，是將三維問題轉(zhuǎn)化為二維問題關(guān)鍵的一步，速度和正確性是首先要關(guān)注的問題。地質(zhì)剖面圖一般是沿勘探線上的平面切割，切割面上的所有三角面片，它們的平面方程是相同的，這樣的平面可以采用一個方程來表示，如：

根據(jù)這一重要特點，采用隱函數(shù)實現(xiàn)平面切割算法，省去了大量幾何運算，提高了運算速度。

隱函數(shù)具有以下3個重要性質(zhì)。

性質(zhì)1：可以描述一些基本幾何形態(tài)，包括平面、球、圓柱、圓錐、橢球體和二次曲面；

性質(zhì)2： 將三維歐式空間分為3個不同區(qū)域，即在隱函數(shù)里、上和外3個區(qū)域，這些區(qū)域分別被定義為F(x,y,z)＜0,F(x,y,z)=0,F(x,y,z)＞0。

性質(zhì)3： 可以將空間位置化為一個數(shù)值，如將空間的坐標(biāo)(xi,yi,zi)轉(zhuǎn)變?yōu)橐粋€數(shù)值ti。

利用隱函數(shù)的性質(zhì) 3，并結(jié)合等值面(線)算法(如Marching cubes等算法)，很容易將切割線與被切實體之間的交線求出。

3.1.2 基于平面隱函數(shù)的實體切割

基本思想是：通過隱函數(shù)為被切實體所有的頂點坐標(biāo)產(chǎn)生一個數(shù)值，然后求F(x,y,z)=0的等值線(根據(jù)隱函數(shù)性質(zhì)2)，即為所求的交線?；静襟E如下：

步驟1: 通過隱函數(shù)(1)為被切實體的所有頂點產(chǎn)生數(shù)值；

步驟2: 求F(x,y,z)=0的等值線，即為交線，形成一個或多個多邊形區(qū)域。

對于每一個三角形面片，如果它的每一個三角形都為正或為負(fù)，則該三角形不會被切割面切割，也不會形成交線，然而，如果這些數(shù)值既有正也有負(fù)，則該三角形被切割，必然有交線通過該三角形，交點位于三角形的邊上，通過沿邊線線性插值計算交點坐標(biāo)。三角形被切割可以分為如圖3所示8種情況。對于每個三角形，根據(jù)其3個頂點的數(shù)值找到對應(yīng)8種情況的一種情況，找到交點所在的邊。如果其中一個交點所在邊的2個頂點對應(yīng)的F(x,y,z)數(shù)值為v1，v2。則交點的數(shù)值為：

則交點坐標(biāo)為：

圖3 切剖三角形的8種情形Fig.3 Eight cases of triangular surface cut

兩交點連線構(gòu)成該三角形和平面的切割線。

3.2 基于KD樹的輪廓線生成算法

如果一個平面切割一個封閉的實體，如存在交線，則所有交線合在一起將形成1個或多個封閉的交線圈(多邊形)。三角形和平面的交線是一系列空間離散的直線段，不是閉合的輪廓線，必須將前后相鄰的線段串聯(lián)為1個或多個多邊形區(qū)域，以便后續(xù)處理。

由n條孤立的直線段形成閉合的多段線，最簡單的方式是：取出1條直線段的1個端點與其他剩下的直線段端點比較，如果相等，則將被找到的直線和前面的直線段連接，然后取出被找到直線的另一個端點，繼續(xù)查找，直到形成一閉合多邊形或搜索完畢。該方法速度比較慢，最壞情況下時間復(fù)雜度為O(n2)。

本文基于KD樹查找直線段的鄰接關(guān)系，快速實現(xiàn)閉合多段線的建立。

3.2.1 KD樹簡介

KD樹[10?11]是二叉搜索樹的擴(kuò)展，K表示空間維數(shù)，故又稱K維搜索樹，特別適合空間點狀目標(biāo)的索引，其示意圖見圖4。在每個內(nèi)部節(jié)點中，用一個K?1維與坐標(biāo)軸平行的超平面(如二維空間中的線)將節(jié)點所表示的K維空間分成2個部分，存儲在子樹中的點大約一半落入一側(cè)，而另一半落入另一側(cè)。當(dāng)一個節(jié)點中的點數(shù)少于給定的最大點數(shù)時，劃分結(jié)束。這些超平面在K個可能的方向上交替出現(xiàn)，每個超平面中至少包括一個數(shù)據(jù)點，這就保證了KD樹不存在“無點空間”的浪費。

圖4 KD樹示意圖(K=2)Fig.4 KD tree diagram

3.2.2 KD樹構(gòu)建及索引

為了構(gòu)建平衡樹，本文給出交點數(shù)為n時KD樹構(gòu)建算法。

步驟 1：初始化一個新結(jié)點 TopNode，將所有點作為其點數(shù)據(jù)；

步驟 2：判斷點數(shù)據(jù)個數(shù)是否小于最小個數(shù)，如果小于最小個數(shù)，則結(jié)束；

步驟 3：沿數(shù)據(jù)的最大延伸方向?qū)?shù)據(jù)排序，采用過中間的位置的點且垂直于最大延伸方向的超平面分割數(shù)據(jù)，形成2個新左子結(jié)點LeftNode和右子結(jié)點RightNode；

步驟 4：對LeftNode執(zhí)行步驟2；

步驟 5：對RightNode執(zhí)行步驟2。

為了便于查詢，KD樹記錄其覆蓋范圍最大、最小值，其查詢算法和二叉樹查找相同，在此不再贅述，其時間復(fù)雜度為O(nlog2n)。

3.3 基于樹結(jié)構(gòu)的封閉區(qū)域組織和索引

基本原理是采用樹結(jié)構(gòu)來表達(dá)區(qū)域的存儲和組織方式，以下稱為區(qū)域樹[11]。區(qū)域樹有如下特點：

(1)一個封閉區(qū)域映射一個區(qū)域樹。

(2)區(qū)域樹的結(jié)點表示區(qū)域的內(nèi)、外邊界，其根結(jié)點表示該區(qū)域的最外圍邊界，子結(jié)點表示區(qū)域內(nèi)部的孔、島邊界。

(3)結(jié)點之間的上下層次關(guān)系表示邊界的包含關(guān)系，同層結(jié)點表示邊界的并列關(guān)系。

(4)邊界是有方向性的，并對邊界方向作如下規(guī)定：區(qū)域最外圍邊界的走向為逆時針方向，內(nèi)部孔洞邊界沿順時針方向，島邊界沿逆時針方向。

某地質(zhì)體切割輪廓線示意圖如圖5所示。該地質(zhì)體輪廓線將平面劃分為 3個區(qū)域(包含最外的無窮大區(qū)域∞)。目標(biāo)區(qū)域及其對應(yīng)的樹如圖6所示。

圖5 某地質(zhì)體輪廓線示意圖Fig.5 Schematic outline of a geological

圖6 區(qū)域及其對應(yīng)的樹結(jié)構(gòu)Fig.6 Regions expressed by tree

區(qū)域的提取方法有多種[12?14]，其中文獻(xiàn)[13]給出了定向極大、極小閉環(huán)的概念，通過建立結(jié)點?路徑網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)圖，采用內(nèi)層路徑優(yōu)先查找算法，提取所有定向閉環(huán)，最后建立定向極小閉環(huán)和定向極大閉環(huán)的隸屬關(guān)系，建立區(qū)域樹，有效地實現(xiàn)了復(fù)雜區(qū)域的識別和提取，具有普遍適用性。本文采用文獻(xiàn)[13]的算法實現(xiàn)了復(fù)雜輪廓線的封閉區(qū)域提取和識別，當(dāng)存在l個極大閉環(huán)，k個極小閉環(huán)時，該算法時間復(fù)雜度為O(l×k)。

區(qū)域樹的構(gòu)造過程見參考文獻(xiàn)[13]。區(qū)域樹的構(gòu)造為地質(zhì)體剖面圖矢量圖案填充提供了數(shù)據(jù)準(zhǔn)備。

4 算法分析和應(yīng)用實例

4.1 算法復(fù)雜度分析

該算法時間復(fù)雜度主要由3個部分決定：基于平面隱函數(shù)的實體切割，基于KD樹的閉合輪廓線構(gòu)造和區(qū)域樹的構(gòu)造。假如構(gòu)成復(fù)雜地質(zhì)體三角面片為m，平面切割后生成n條孤立直線段，所有直線段構(gòu)成l個極大閉環(huán)，k個極小閉環(huán)時，則總的時間復(fù)雜度為：。

4.2 應(yīng)用實例

通過 VC+開發(fā)工具和可視化工具包實現(xiàn)了該算法。運行環(huán)境為Windows 2000操作系統(tǒng)，CPU為Inter Pentium4 3.00 GHz，內(nèi)存512 MB。當(dāng)復(fù)雜地質(zhì)體三角面片數(shù)為1 589，3 264，5 470，7 412，10 345時，根據(jù)勘探線剖切形成的離散直線段、定向閉環(huán)和切割三角面片耗時見表1，并做出三角形個數(shù)?切割時間曲線，如圖7所示。由圖7可以看出，該三角形切割算法時間復(fù)雜度近似線性。

表1 不同三角面片個數(shù)算法耗時Table 1 Consuming time for different number of triangles

圖7 三角面片數(shù)?時間曲線Fig.7 Time changes with different number of triangles

“DIMINE”數(shù)字礦山軟件是中南大學(xué)數(shù)字礦山研究中心開發(fā)的一套礦床地質(zhì)建模、礦山工程測量、數(shù)據(jù)處理和地下、露天礦床開采設(shè)計等功能于一體的集成化、智能化、可視化軟件系統(tǒng)。該系統(tǒng)中地質(zhì)平剖面投影圖的生成、縱剖面圖的繪制等模塊中采用了該算法。該算法可有效地用于復(fù)雜地質(zhì)體剖面圖自動生成操作，如圖8所示。

圖8 地質(zhì)體剖面圖自動生成示例Fig.8 An example of automatic generation of Geological Profile

5 結(jié)論

(1)提出平面隱函數(shù)的實體切割算法，該算法將三角形和平面間的幾何運算轉(zhuǎn)換為數(shù)學(xué)運算，加快了運算速度，實現(xiàn)了地質(zhì)實體的快速剖切。

(2)提出基于 KD樹封閉輪廓線生成方法，同過對鄰近點對的快速索引，實現(xiàn)了封閉輪廓線的快速生成。該方法同樣適用于開口實體(如地層、地表)的非閉合輪廓線生成。

(3)采用圖論中的樹結(jié)構(gòu)形式化表達(dá)區(qū)域的空間組織結(jié)構(gòu)，為巖性圖案填充提供了數(shù)據(jù)來源。

(4)為進(jìn)一步提高運算速度，可建立三角形之間的鄰接關(guān)系，從而在生成離散線段的同時，快速建立線段間的相鄰關(guān)系，是下一步研究的難點和方向。

[1]曹燮明.采礦手冊[M].北京: 冶金工業(yè)出版社,1999:172?175.CAO Xie-ming.Mining guide[M].Beijing: Metallurgical Industry Press,1999: 172?175.

[2]Shostko A A,Lohner R,Sandberg W C.Surface triangulation over intersecting geometries[J].International Journal for Numerical Methods in Engineering,1999,44(9): 1359?1376.

[3]Coelho L C G,Gattass M,Fiqueiredo L H.Intersecting and trimming parametric meshes on finite-elements shells[J].International Journal for Numerical Methods in Engineering,2000,47(4): 777?800.

[4]Lindenbecka C H,Ebertb H D,Ulmera H.A program to clip triangle meshes using the rapid and triangle libraries and the visualization toolkit[J].Computers & Geosciences,2002,28(6):841?850.

[5]陳國量,劉修國,商建嘎,等.三維地質(zhì)結(jié)構(gòu)模型的切割分析技術(shù)及方法[J].計算機(jī)工程,2007,33(20): 184?186.CHEN Guo-liang,LIU Xiu-guo,SHANG Jian-ga,et al.Incision analysis technology and method for 3D geological structure model[J].Computer Engineering,2007,33(20): 184?186.

[6]陳建宏,周智勇,陳綱,等.基于鉆孔數(shù)據(jù)的勘探線剖面圖自動生成方法[J].中南大學(xué)學(xué)報: 自然科學(xué)版,2005,36(3):487?490.CHEN Jian-hong,ZHOU Zhi-yong,CHEN Gang,et al.Automatic formation method of prospecting line profile map based on drill hole database[J].Journal of Central South University: Science and Technology,2005,36(3): 487?490.

[7]朱瑩,劉學(xué)軍,陳鎖忠,等.基于GIS的地質(zhì)剖面圖自動繪制軟件的研究[J].南京師大學(xué)報: 自然科學(xué)版,2007,30(4):104?108.ZHU Ying,LIU Xue-jun,CHEN Suo-zhong,et al.A research on automatic generating software of geologic section based on GIS[J].Journal of Nanjing Normal University: Natural Science Edition,2007,30(4): 104?108.

[8]王建芳,包世泰,余應(yīng)剛,等.基于GIS模板的地質(zhì)剖面圖模型及其實現(xiàn)[J].測繪科學(xué),2008,33(5): 184?186.WANG Jian-fang,BAO Shi-tai,YU Ying-gang,et al.Realization of geological section map model based GIS template[J].Science of Surveying and Mapping,2008,33(5): 184?186.

[9]王家耀.空間信息系統(tǒng)原理[M].北京: 科學(xué)出版社,2003:100?105.WANG Jia-yao.Principles of spatial information systems[M].Beijing: Science Press,2003: 100?105.

[10]Bentley J L.Multidimensional binary search trees used for associative searching[J].Commun ACM,1975,18(9): 509?517.

[11]Bentley J L.K-d trees for semi-dynamic point sets[C]//Proc 6th Annual Symposium on Computational Geometry.New York:Association for Computing Machinery,1990: 187?197.

[12]譚正華,王李管,畢林,等.參數(shù)曲線集復(fù)雜區(qū)域的全自動識別算法[J].計算機(jī)工程,2010,36(8): 23?26.TAN Zheng-hua,WANG Li-guan,BI Lin,et al.Automatic identification algorithm for complicated regions of parameterized curve set[J].Computer Engineering,2010,36(8):23?26.

[13]宋懷波,路長厚,王富春.一種新的復(fù)雜區(qū)域孔洞填充算法[J].桂林電子科技大學(xué)學(xué)報,2006,26(6): 451?454.SONG Huai-bo,LU Chang-hou,WANG Fu-chun.A new hole-filling algorithm for complicated connecting regions[J].Journal of Guilin University of Electronic Technology,2006,26(6): 451?454.