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

胡長濤,張志華,吳智慧

(1. 蘭州交通大學 測繪與地理信息學院,蘭州 730070;2. 地理國情監測技術應用國家地方聯合工程研究中心,蘭州 730070;3. 甘肅省地理國情監測工程實驗室,蘭州 730070)

隨著三維GIS發展及地下工程的需要,地質模型的構建、可視化分析成為地學領域研究熱點[1-3]。三維地質建模是基于現代空間信息理論利用地質調查數據實現地質信息的管理、解譯、空間分析及可視化,可直觀展現地層厚度、巖性及空間形態等屬性[4-5]。高質量的地質模型可有效進行地學空間分析與數值計算,為資源勘探及城市地下空間規劃提供可靠數據支撐[6]。

空間體元模型作為三維地質模型可視化的基礎,在地質資源表達評估及后續定量分析中起著至關重要作用[7]。近年來三維地質建模的方法主要有基于面模型、基于體模型和基于混合模型3大類建模體系[8]。基于面元模型建模方法便于數據更新及可視化,但缺乏內部屬性的表達,無法進行空間查詢分析。混合模型體拓撲關系復雜,技術實現難度大,當前建模方法研究多以不規則體元模型為主。用于地質建模方面的常見不規則體元有廣義三棱柱、六面體和四面體等[9-10]。相比于其它不規則體元,四面體具有體元結構簡單、快速幾何變換、拓撲關系清晰的優點,可以模擬規則與不規則地質體[11]。地質模型以四面體為基本體元,可實現不規則復雜曲面地質體模型的構建。

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

1 建模方法與流程

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

圖1 地質建模流程

2 三維地下空間約束表達

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

圖2 PLC范例

Poly格式文件是PLC的一種線框模型描述,可理解為對分段線性復合體的數字化描述。它由4個部分組成,分別是點列表、面列表、孔洞列表(無空洞即為0)和可選的區域屬性列表。在地質建模中,點列表即為地層邊界面節點,面列表為三角網的三角面,Poly文件適用于地層約束模型的存儲。

3 約束Delaunay四面體剖分算法

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

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

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

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

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

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

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

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

3.2 約束Delaunay四面體剖分

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

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

圖3 約束Delaunay四面體示意圖

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

3.3 特征約束四面體剖分

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

圖4 四面體形狀參數圖

圖5 薄元四面體示意圖

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

4 實驗分析

研究區位于陜西省某地區,寬度為706 m,長度為1 811 m。其中鉆孔數據數量為39個,地層分界面為5個,地層數量為4層。先對鉆孔數據加密,分地層構建多地層表面模型生成約束PLC。利用TetGen進行特征約束Delaunay四面體剖分,向約束邊界和地層內部插入節點,恢復約束線和約束面,構建該地區四面體格網地質模型,并對不同約束下的地質模型進行四面體質量分析。

4.1 地層表面約束模型構建

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

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

3)地層模型側邊界構建。地質體側邊生成即利用三角網縫合上下底層曲面,其中基于平行輪廓線的TIN構建方法應用最為廣泛。縫合的原理為把兩個相鄰的平行輪廓線上的點連接生成滿足一定要求的不相交三角形,實現單個地層模型的構建。基于平行輪廓線的TIN構建方法有最大體積法、最短對角線法和同步前進法[23]。文中研究的地質體對象為正常產狀的巖層,輪廓線大小形狀相近,應用最短對角線法構建的三維表面效果較好。

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

圖6 最短對角線法示意

4.2 三維地質模型可視化

VTK是一個開源的、支持多平臺的軟件系統,主要應用于三維計算機圖形學和三維可視化。其基于三維函數庫OpenGL和面向對象的原理,通過創建.vtk的數據文件格式,為各種數據集類型提供一致的數據表示方案。VTK支持體元繪制,借助于VTK組件庫可實現三維地質模型的可視化及幾何運算,為實際工程施工及地質災害超前預報提供可靠參考。

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

圖7 三維地質模型渲染

圖8 組合剖面圖

4.3 地質模型四面體質量分析

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

圖9 三維地質模型表面格網

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

表1 地質模型數據統計表

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

圖10 四面體格網長寬比統計

圖11 四面體格網線線角統計

5 結 論

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