基于特征約束的四面體剖分方法研究

胡長濤,張志華,吳智慧

(1. 蘭州交通大學(xué) 測繪與地理信息學(xué)院,蘭州 730070;2. 地理國情監(jiān)測技術(shù)應(yīng)用國家地方聯(lián)合工程研究中心,蘭州 730070;3. 甘肅省地理國情監(jiān)測工程實驗室,蘭州 730070)

隨著三維GIS發(fā)展及地下工程的需要,地質(zhì)模型的構(gòu)建、可視化分析成為地學(xué)領(lǐng)域研究熱點[1-3]。三維地質(zhì)建模是基于現(xiàn)代空間信息理論利用地質(zhì)調(diào)查數(shù)據(jù)實現(xiàn)地質(zhì)信息的管理、解譯、空間分析及可視化,可直觀展現(xiàn)地層厚度、巖性及空間形態(tài)等屬性[4-5]。高質(zhì)量的地質(zhì)模型可有效進行地學(xué)空間分析與數(shù)值計算,為資源勘探及城市地下空間規(guī)劃提供可靠數(shù)據(jù)支撐[6]。

空間體元模型作為三維地質(zhì)模型可視化的基礎(chǔ),在地質(zhì)資源表達評估及后續(xù)定量分析中起著至關(guān)重要作用[7]。近年來三維地質(zhì)建模的方法主要有基于面模型、基于體模型和基于混合模型3大類建模體系[8]?；诿嬖Ｐ徒７椒ū阌跀?shù)據(jù)更新及可視化,但缺乏內(nèi)部屬性的表達,無法進行空間查詢分析?；旌夏Ｐ腕w拓?fù)潢P(guān)系復(fù)雜,技術(shù)實現(xiàn)難度大,當(dāng)前建模方法研究多以不規(guī)則體元模型為主。用于地質(zhì)建模方面的常見不規(guī)則體元有廣義三棱柱、六面體和四面體等[9-10]。相比于其它不規(guī)則體元,四面體具有體元結(jié)構(gòu)簡單、快速幾何變換、拓?fù)潢P(guān)系清晰的優(yōu)點,可以模擬規(guī)則與不規(guī)則地質(zhì)體[11]。地質(zhì)模型以四面體為基本體元,可實現(xiàn)不規(guī)則復(fù)雜曲面地質(zhì)體模型的構(gòu)建。

Delaunay三角剖分具有嚴(yán)密的理論基礎(chǔ)和數(shù)學(xué)特性[12],限定條件下的三維剖分研究比較成熟,但針對地學(xué)領(lǐng)域復(fù)雜的邊界面及內(nèi)部特征約束較多的地理空間對象方面的研究較少。因此,文中提出一種基于地質(zhì)鉆孔數(shù)據(jù),以地層邊界面不規(guī)則三角網(wǎng)(triangulate irregular network, TIN)為約束條件,對其進行四面體剖分的模型構(gòu)建方法。該方法得到可視化效果良好的地質(zhì)模型,可快速實現(xiàn)三維剖面模型及礦體產(chǎn)量估算等空間分析操作,為地下巖礦體開采等施工提供決策數(shù)據(jù)支撐。

1 建模方法與流程

考慮Delaunay三角剖分的優(yōu)良邊界約束特性,文中利用克里金插值法生成虛擬鉆孔數(shù)據(jù)并建立地層表面模型。地層側(cè)面利用最短對角線法構(gòu)建側(cè)邊三角網(wǎng),地層面三角網(wǎng)與側(cè)面三角網(wǎng)構(gòu)成多地層表面模型。以地層表面模型轉(zhuǎn)換為分段線性復(fù)合體作為輸入,進行約束四面體剖分得到以四面體為體元的三維地質(zhì)模型,通過設(shè)置不同半徑邊長比及二面角等特征約束控制格網(wǎng)質(zhì)量,實現(xiàn)地質(zhì)模型的更新與重構(gòu)。模型的基本幾何元素有節(jié)點、邊、三角形面片和簡單塊體4種。建模過程包括地層表面模型構(gòu)建、約束四面體剖分兩部分,地質(zhì)模型的可視化借助于VTK(visualization toolkit)組件庫實現(xiàn),技術(shù)路線如圖1所示。

圖1 地質(zhì)建模流程

2 三維地下空間約束表達

地學(xué)領(lǐng)域地層的邊界面及內(nèi)部約束特征通常較多,選擇適宜的三維空間表達方式有利于構(gòu)建精確完整的地質(zhì)模型。文中地質(zhì)體三維空間的約束表達為分段線性復(fù)合體(piecewise linear complex, PLC)。PLC由Miller等[13]提出,三維PLC中X是一個集頂點、邊、面和多面體的集合,如圖2所示。PLC定義中不允許它的單體非法相交,兩個邊的相交只能在一個公共的頂點處,并且該點位于X內(nèi);四面體T的兩個面相交只能是頂點和線段的結(jié)合,并且位于X內(nèi)。另外,PLC允許點、線、面浮動于多面體區(qū)域內(nèi),在三維表達上靈活,適用于地學(xué)建模領(lǐng)域。

圖2 PLC范例

Poly格式文件是PLC的一種線框模型描述,可理解為對分段線性復(fù)合體的數(shù)字化描述。它由4個部分組成,分別是點列表、面列表、孔洞列表(無空洞即為0)和可選的區(qū)域?qū)傩粤斜怼Ｔ诘刭|(zhì)建模中,點列表即為地層邊界面節(jié)點,面列表為三角網(wǎng)的三角面,Poly文件適用于地層約束模型的存儲。

3 約束Delaunay四面體剖分算法

TetGen是德國柏林魏爾斯特拉斯應(yīng)用分析和隨機研究所開發(fā)的開源高質(zhì)量的四面體格網(wǎng)生成程序[14]。它可生成精確的約束Delaunay四面體化、邊界一致的Delaunay網(wǎng)格和Voronoi分區(qū),可以穩(wěn)健地處理任意復(fù)雜的三維幾何形狀。

3.1 基于頂點插入的增量Delaunay算法

三維Delaunay剖分算法中使用最廣泛的是Bowyer-Watson算法[15],其原理如下:

1)建立一個包含所有點的初始格網(wǎng),格網(wǎng)通常為點集的最小長方體包圍盒;

2)向初始格網(wǎng)中插入點,遍歷尋找包含該點的四面體單元;

3)通過鄰接關(guān)系找出外接球包含的所有四面體單元,刪除這些單元構(gòu)成Delaunay空腔,連接插入點與空腔中的每個頂點生成新的網(wǎng)格單元;

4)重復(fù)2)和3)操作直至所有頂點插入完畢,刪除格網(wǎng)頂點,得到點集Delaunay四面體格網(wǎng)。

另外,TetGen使用空間排序體制法[16]對插入點進行預(yù)排序處理,通過簡單的隨機行走算法[17]進行點定位,提高了算法的運行速度。通過使用orient 3D(用于點定位)和in sphere(用于更新Delaunay四面體)約束方法保證了算法的健壯性。

3.2 約束Delaunay四面體剖分

約束Delaunay四面體剖分(constrained delaunay tetrahedralization, CDT)由Shewchuk[18]提出,是Delaunay四面體剖分的變體。約束Delaunay四面體化的定義如下:

設(shè)X是一個三維PLC,平面f∈X使得P和Q兩點位于這個平面的兩側(cè),線段PQ與這個切面相交,PQ兩個點也位于X內(nèi),但相互不可見。線段不影響可見性,AB是線段,但C和D之間可見。點及平面的空間位置見圖3。對于X中頂點形成的四面體PABC,由于它的外接球不包括X中對其可見的點,則四面體PABC為約束Delaunay四面體。其中f相當(dāng)于約束面。X的四面體剖分是一個三維單純復(fù)形T,T和X的頂點一致且X中的每個體元都是T中單純形的并集X。對于X的約束四面體剖分T,如果滿足每個四面體都符合約束Delaunay,則T就為X的約束Delaunay四面體剖分。

圖3 約束Delaunay四面體示意圖

Delaunay四面體剖分和約束Delaunay四面體剖分的不同之處在于,對邊界多邊形內(nèi)的三角形去掉了局部Delaunay的要求。因此,約束Delaunay四面體剖分保留了很多Delaunay四面體剖分的有利屬性[19]。與Delaunay四面體剖分相比,CDT生成所需的Steiner點通常會更少。尤其是當(dāng)PLC包含一些銳化特征,幾個單體會在一個很小的角度(或是二面角)相交時,Delaunay四面體剖分是難以生成的,而且可能需要大量的Steiner點,而CDT可以很容易地處理這種銳化特征。

3.3 特征約束四面體剖分

特征約束是對四面體體元的形狀參數(shù)化約束控制,實現(xiàn)四面體格網(wǎng)的質(zhì)量優(yōu)化。四面體形狀參數(shù)沒有唯一的定義,一般的含義為避免內(nèi)部存在非常小或非常大的二面角。對于一個單純形四面體質(zhì)量判定的最普遍方法是長寬比。圖4為單個四面體形狀參數(shù)圖,其中r為四面體外接球半徑,h為點到對應(yīng)三角面的距離。四面體PABC的長寬比是最長邊的長度|AC|與最小高度h之間的比值,其中最小高度h是四面體4個頂點到對應(yīng)平面距離的最小值。TetGen并不是用長寬比,而是選擇半徑邊長比和二面角(兩個面的夾角)作為四面體的形狀控制,二者相互結(jié)合可實現(xiàn)對四面體的特征約束。四面體T的半徑-邊長比(radius-edge ratio)是外接球半徑r與四面體最短邊長度d的比值,即為r與|PC|的比值。大多數(shù)形狀不好的四面體會有一個比較大的半徑-邊長比值,只有沒有短邊且體積接近0的扁平四面體薄片[20]會存在較小的半徑邊長比。薄元四面體形狀如圖5所示。TetGen采用四面體的最小二面角作為Delaunay細(xì)化算法的第二個約束方法,通過Delaunay細(xì)化迭代次數(shù)來移除薄片。由于具有銳化特征的四面體是永遠(yuǎn)不會被移除的,因此只有少數(shù)畸形四面體會保留下來。

圖4 四面體形狀參數(shù)圖

圖5 薄元四面體示意圖

TetGen采用了一個簡單的“爬山”方案來控制四面體形狀特征,即只有當(dāng)生成的新四面體參數(shù)值優(yōu)于當(dāng)前最差的四面體時才進行局部優(yōu)化運算。初始化形狀較差的四面體列表,若其參數(shù)值小于給定的目標(biāo)值,對其使用局部操作:使用邊/面翻轉(zhuǎn)和頂點平滑來刪除它們。結(jié)合這兩種操作,不斷迭代實現(xiàn)優(yōu)化的四面體替換劣質(zhì)四面體。

4 實驗分析

研究區(qū)位于陜西省某地區(qū),寬度為706 m,長度為1 811 m。其中鉆孔數(shù)據(jù)數(shù)量為39個,地層分界面為5個,地層數(shù)量為4層。先對鉆孔數(shù)據(jù)加密,分地層構(gòu)建多地層表面模型生成約束PLC。利用TetGen進行特征約束Delaunay四面體剖分,向約束邊界和地層內(nèi)部插入節(jié)點,恢復(fù)約束線和約束面,構(gòu)建該地區(qū)四面體格網(wǎng)地質(zhì)模型,并對不同約束下的地質(zhì)模型進行四面體質(zhì)量分析。

4.1 地層表面約束模型構(gòu)建

1)鉆孔數(shù)據(jù)預(yù)處理。原始鉆孔數(shù)據(jù)一般包含坐標(biāo)、巖性、深度及方位角等信息。建立地層面Delaunay三角網(wǎng)需要獲得鉆孔軌跡線上巖性分界點的三維點坐標(biāo),故需要對原始數(shù)據(jù)進行預(yù)處理。由于研究區(qū)常存在鉆孔資料分布不均,采用克里金插值算法對原始鉆孔數(shù)據(jù)進行加密處理。通過插值處理計算新增巖性分界點坐標(biāo),構(gòu)建相對精確完整的地層分界面。

2)地層模型上下邊界構(gòu)建。地質(zhì)界面的建模本質(zhì)上是利用面模型實現(xiàn)對地層起伏形態(tài)的模擬。目前地層邊界面的建模主要采用不規(guī)則三角網(wǎng)法[21]。TIN模型結(jié)構(gòu)簡單,建模理論成熟,精度高,可視化效果好,適用于地質(zhì)體表面的幾何建模[22]。文中采用Delaunay三角剖分算法構(gòu)建TIN,TIN由節(jié)點、三角形邊和三角形面構(gòu)成。在地質(zhì)建模中,地層分界點即為節(jié)點,對節(jié)點進行Delaunay三角剖分,依次建立各地層上下邊界面的三維模型。

3)地層模型側(cè)邊界構(gòu)建。地質(zhì)體側(cè)邊生成即利用三角網(wǎng)縫合上下底層曲面,其中基于平行輪廓線的TIN構(gòu)建方法應(yīng)用最為廣泛?？p合的原理為把兩個相鄰的平行輪廓線上的點連接生成滿足一定要求的不相交三角形,實現(xiàn)單個地層模型的構(gòu)建?；谄叫休喞€的TIN構(gòu)建方法有最大體積法、最短對角線法和同步前進法[23]。文中研究的地質(zhì)體對象為正常產(chǎn)狀的巖層,輪廓線大小形狀相近,應(yīng)用最短對角線法構(gòu)建的三維表面效果較好。

最短對角線法具體步驟為:獲取上下輪廓線,判斷上下輪廓線頂點數(shù)量是否相等,若不相等,在數(shù)量少的輪廓線中最長距離的兩個頂點之間外推一個頂點,重復(fù)以上操作直至數(shù)量相等。最短對角線法的空間示意圖見圖6,在上輪廓線上找到與Pi點距離最近的下輪廓點Qi,判斷Pi,Pi+1,Qi,Qi+1組成的四邊形對角線PiQi+1和Pi+1Qi長度大小;若|PiQi+1|>|Pi+1Qi|,則構(gòu)建△PiPi+1Qi和△Pi+1QiQi+1,否則構(gòu)建△PiPi+1Qi+1和△PiQiQi+1;直到完成所有頂點判斷,即實現(xiàn)地層側(cè)邊的生成。

圖6 最短對角線法示意

4.2 三維地質(zhì)模型可視化

VTK是一個開源的、支持多平臺的軟件系統(tǒng),主要應(yīng)用于三維計算機圖形學(xué)和三維可視化。其基于三維函數(shù)庫OpenGL和面向?qū)ο蟮脑?通過創(chuàng)建.vtk的數(shù)據(jù)文件格式,為各種數(shù)據(jù)集類型提供一致的數(shù)據(jù)表示方案。VTK支持體元繪制,借助于VTK組件庫可實現(xiàn)三維地質(zhì)模型的可視化及幾何運算,為實際工程施工及地質(zhì)災(zāi)害超前預(yù)報提供可靠參考。

對多地層表面模型進行約束四面體剖分輸出.vtk文件,使用VTK組件庫實現(xiàn)以四面體為基本體元的地質(zhì)模型的渲染,地質(zhì)體模型可視化如圖7所示。另外在地學(xué)分析中,三維地質(zhì)切割是三維建模的重要信息交互手段,能夠分析礦體的內(nèi)部結(jié)構(gòu),揭示地質(zhì)體內(nèi)部隱藏的地質(zhì)信息[24]。三維地質(zhì)模型的切割主要分為任意平面切割和組合切割。任意平面切割是最廣泛的切割方式,組合切割是利用多個平面組合構(gòu)建出三維地質(zhì)剖面模型。地質(zhì)模型組合剖面圖如圖8所示。借助VTK組件庫可快速實現(xiàn)模型切割及剖面可視化。

圖7 三維地質(zhì)模型渲染

圖8 組合剖面圖

4.3 地質(zhì)模型四面體質(zhì)量分析

三維地質(zhì)模型質(zhì)量主要與四面體質(zhì)量相關(guān)。在相同地層表面模型條件下,不同特征參數(shù)約束所生成的四面體格網(wǎng)質(zhì)量不同。通過設(shè)置3種不同特征參數(shù)約束條件,分析四面體格網(wǎng)的質(zhì)量情況。3種具體參數(shù)約束為半徑邊長比為4.0,同時最小二面角為8°;半徑邊長比為2.0,同時最小二面角為8°;以及半徑邊長比為2.0,同時最小二面角為20°的約束,在此分別簡稱為參約1、參約2和參約3。3種參數(shù)約束下的地質(zhì)模型表面格網(wǎng)如圖9所示。

圖9 三維地質(zhì)模型表面格網(wǎng)

地質(zhì)模型以四面體為基本體元,可根據(jù)四面體頂點及地層屬性等信息快速實現(xiàn)地層體積計算。通過對實例模型統(tǒng)計計算,可得底層地層的體積為2.01×108m3,第二層地層的體積較小為1.18×107m3,第三層地層的體積為1.99×108m3,最頂層地層的體積為2.61×108m3,總體積為6.73×108m3。對3種參約剖分后的四面體格網(wǎng)的節(jié)點數(shù)量、四面體數(shù)量和模型內(nèi)存大小統(tǒng)計得到表1。由表1可知,在參約1條件下模型的四面體數(shù)量最少,模型節(jié)點數(shù)量最少,占用儲存空間最小。另外通過圖9(a)可知,參約1模型表面三角網(wǎng)稀疏,間接反映四面體格網(wǎng)質(zhì)量較差;參約3條件下,模型四面體數(shù)量最多,模型節(jié)點數(shù)量多,表面格網(wǎng)三角形分布均勻,模型質(zhì)量較高,但占用儲存空間最大。

表1 地質(zhì)模型數(shù)據(jù)統(tǒng)計表

為分析不同剖分參數(shù)下四面體格網(wǎng)的整體質(zhì)量,文中以四面體長寬比和線線角兩個評價標(biāo)準(zhǔn)定量分析四面體格網(wǎng)質(zhì)量。四面體長寬比越小代表四面體質(zhì)量越好。線線角為四面體三角面中兩條線的夾角,范圍為0°～180°,可間接反映四面體質(zhì)量情況。由于不同參約下模型四面體總數(shù)不同,故將長寬比和線線角劃范圍按百分比進行統(tǒng)計,結(jié)果如圖10和圖11所示。通過圖10可知,參約3條件下在小于4的長寬比范圍內(nèi),經(jīng)統(tǒng)計四面體數(shù)占量比達85%,明顯大于另外兩種參數(shù),四面體格網(wǎng)質(zhì)量更高。在參約3線線角集中在40°～80°,占比達79%,遠(yuǎn)高于其它兩參數(shù)的42%和57%。線線角集中在0°或180°附近意味著位于薄元四面體上。由圖11可知,參約3下在0°或180°附近的線線角數(shù)量很少,薄元四面體數(shù)量最少。通過對比3種參約剖分結(jié)果,參約3下模型四面體整體質(zhì)量最高,參約2次之,參約1下模型四面體整體質(zhì)量最差。本方法利用特征約束可生成不同質(zhì)量要求的四面體模型,可滿足不同地質(zhì)分析條件的需要。

圖10 四面體格網(wǎng)長寬比統(tǒng)計

圖11 四面體格網(wǎng)線線角統(tǒng)計

5 結(jié) 論

基于鉆孔數(shù)據(jù)建立Delaunay三角網(wǎng)生成多地層地質(zhì)體邊界模型,并對其進行約束Delaunay四面體剖分,通過半徑邊長比和最小二面角大小等特征約束生成網(wǎng)格質(zhì)量較高的三維地質(zhì)模型。該模型以四面體為基本體元,可進行數(shù)據(jù)存儲與快速可視化,并能夠進行拓?fù)潢P(guān)系查詢、體積計算、剖切面生成、模型切割及有限元分析等空間分析操作。通過實驗分析證明了建模方法的有效性和可行性,可為地下礦體開采等工程規(guī)劃提供數(shù)據(jù)支持?？紤]到地下空間點線面及空洞的復(fù)雜性,還需對斷層、尖滅等復(fù)雜地質(zhì)體建模進行研究,并將模型應(yīng)用于地學(xué)數(shù)值計算分析,為地質(zhì)災(zāi)害模擬及預(yù)測提供借鑒。