約束Delaunay四面體剖分在三維地質(zhì)建模中的應(yīng)用

余淑娟，郭飛，李想，徐峰

（南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210046）

0 引言

分析與解決城市規(guī)劃建設(shè)中的環(huán)境地質(zhì)問(wèn)題，已成為城市發(fā)展的當(dāng)務(wù)之急。近年來(lái)，隨著計(jì)算機(jī)技術(shù)的發(fā)展和城市規(guī)劃建設(shè)的需要，利用計(jì)算機(jī)建立科學(xué)、精細(xì)、可靠的三維淺層空間地質(zhì)體模型，指導(dǎo)城市規(guī)劃建設(shè)，逐漸成為計(jì)算機(jī)和地學(xué)領(lǐng)域內(nèi)研究的重要內(nèi)容［1］。

三角剖分，尤其是限定的三角剖分一直是近年來(lái)計(jì)算幾何及其相關(guān)領(lǐng)域的主要研究課題之一。Delaunay三角網(wǎng)具有優(yōu)良的幾何特性，如空外接圓性質(zhì)、最小角最大的性質(zhì)、局部修改性以及良好的邊界約束適應(yīng)能力等，是公認(rèn)的最優(yōu)三角網(wǎng)。Delaunay三角化具有很好的理論基礎(chǔ)和數(shù)學(xué)特性，一直在網(wǎng)格剖分中占有重要地位。對(duì)于沒(méi)有限定條件的離散點(diǎn)而言，Delaunay三角剖分的算法已十分成熟，但在引入限定條件后，其剖分算法還是一個(gè)研究熱點(diǎn)，并不是很成熟。對(duì)約束三角網(wǎng)在三維地質(zhì)建模應(yīng)用中而言，郝海森等利用含強(qiáng)約束條件的Delaunay不規(guī)則三角網(wǎng)構(gòu)建地層底板的三維空間實(shí)體模型［2］；蔡強(qiáng)等采用區(qū)域子分和聯(lián)動(dòng)剖分算法實(shí)現(xiàn)了地質(zhì)結(jié)構(gòu)重疊域的限定Delaunay三角剖分［3］；陳永鋒等對(duì)傳統(tǒng)凸包生成、三角剖分和空外接圓檢測(cè)等算法進(jìn)行了改進(jìn)，對(duì)地質(zhì)面進(jìn)行Delaunay三角剖分［4］；孟永東等依據(jù)帶約束Delaunay三角剖分算法分別生成地層、斷裂、界線類(lèi)地質(zhì)結(jié)構(gòu)的TIN模型，拼合形成工程整體三維地質(zhì)模型［5］。綜合考察現(xiàn)有的研究成果［2－10］，大多數(shù)集中于通過(guò)Delaunay三角剖分構(gòu)建地質(zhì)結(jié)構(gòu)TIN模型，對(duì)于利用地質(zhì)層面Delaunay三角網(wǎng)為約束條件直接構(gòu)建四面體地質(zhì)體模型方面，目前還需要進(jìn)一步研究。

本文考慮到Delaunay三角剖分良好的邊界約束適應(yīng)能力，基于鉆孔數(shù)據(jù)插值處理，提出了一種以地質(zhì)層面Delaunay三角網(wǎng)為約束條件的四面體剖分地質(zhì)建模方法。首先，將研究區(qū)域邊界離散，根據(jù)研究區(qū)域內(nèi)鉆孔勘察資料的空間分布，自適應(yīng)的在數(shù)據(jù)相對(duì)稀疏的地區(qū)添加用于插值處理的虛擬鉆孔，采用克里金插值算法計(jì)算出虛擬鉆孔中和邊界點(diǎn)的地層底板埋深數(shù)據(jù)；其次，對(duì)原始鉆孔和虛擬鉆孔數(shù)據(jù)以及限定的研究區(qū)域邊界構(gòu)建每個(gè)地層面的Delaunay三角網(wǎng)；再次，以地層面三角網(wǎng)和鉆孔數(shù)據(jù)及其空間分布信息為約束，構(gòu)造地質(zhì)幾何模型；然后，以該幾何模型作為輸入的分段連續(xù)線性組合物（Piecewise Linear Complexes，PLC）［11］限定條件，設(shè)置四面體的體積或者半徑邊長(zhǎng)比進(jìn)行網(wǎng)格質(zhì)量控制，構(gòu)建三維地質(zhì)四面體實(shí)體模型，并借助VTK（Visualization Toolkit）組件庫(kù)實(shí)現(xiàn)了該地質(zhì)模型的可視化。

1 三維地質(zhì)模型

三維地質(zhì)建模（3DGeological Modeling）也可稱(chēng)為三維地學(xué)建模、三維地質(zhì)數(shù)字化建模等，指在獲取原始地質(zhì)勘測(cè)數(shù)據(jù)（如地質(zhì)點(diǎn)、鉆孔、平硐、航拍片等）的基礎(chǔ)上，通過(guò)對(duì)地質(zhì)實(shí)體對(duì)象拓?fù)湫畔⒌慕庾g和幾何形態(tài)的確定，將各種信息綜合形成一個(gè)復(fù)雜整體三維模型的過(guò)程。三維地質(zhì)建模技術(shù)就是為了解決地學(xué)領(lǐng)域中遇到的三維問(wèn)題而提出的，如三維地層、斷裂、礦體和巷道的顯示以及三維巷道的空間拓?fù)浞治觥⑷S礦體的體積與儲(chǔ)量計(jì)算等問(wèn)題。

三維地質(zhì)建模常用方法包括塊段構(gòu)模法、線框構(gòu)模法、表面構(gòu)模法及實(shí)體構(gòu)模法［12］。有很多方法可以用來(lái)表達(dá)表面，如等高線模型、網(wǎng)格模型等，最常用的是TIN模型。實(shí)體構(gòu)模法包括四面體法、六面體法及三棱柱法等，目前利用三棱柱體元構(gòu)建地質(zhì)體模型的研究比較成熟［13－15］。本文針對(duì)城市地下淺層空間，以地質(zhì)層面Delaunay三角網(wǎng)及鉆孔點(diǎn)為約束條件進(jìn)行四面體剖分實(shí)現(xiàn)層狀地質(zhì)實(shí)體的三維地質(zhì)體模型方法屬于實(shí)體構(gòu)模法，該方法能夠得到可視化效果較好的地質(zhì)模型，并可高效地支撐地學(xué)空間分析與計(jì)算，為城市規(guī)劃提供更好的決策支持。

2 算法描述

2.1 鉆孔數(shù)據(jù)預(yù)處理與插值處理

本文以地質(zhì)層面Delaunay三角網(wǎng)為約束條件進(jìn)行四面體地質(zhì)建模的方法，對(duì)鉆孔數(shù)據(jù)的格式有一定的要求，必須對(duì)原始鉆孔數(shù)據(jù)進(jìn)行處理。原始的鉆孔數(shù)據(jù)點(diǎn)在縱向是深度遞增值，根據(jù)孔口坐標(biāo)與地層分界點(diǎn)的水平偏移量（或是方位角），計(jì)算鉆孔數(shù)據(jù)點(diǎn)（地層分界點(diǎn)或切穿點(diǎn)）的水平坐標(biāo)，同時(shí)給這些點(diǎn)賦予地層屬性，為后續(xù)三維構(gòu)模提供方便。

由于建模研究區(qū)域內(nèi)鉆孔資料分布極不均衡，直接利用原始的鉆孔資料難以構(gòu)建相對(duì)精確、完整且具有良好可視化效果的三維地質(zhì)體模型。因此，需要對(duì)原始分布不均勻的鉆孔數(shù)據(jù)進(jìn)行適當(dāng)?shù)目臻g插值處理以獲取地層層面信息。本文根據(jù)用戶(hù)指定的添加插值點(diǎn)的間距，離散研究區(qū)域邊界，添加邊界插值點(diǎn)；規(guī)則離散研究區(qū)域內(nèi)部，添加內(nèi)部插值點(diǎn)；根據(jù)原始鉆孔勘察資料分布情況，將靠近（距離小于用戶(hù)指定的間距值）原始鉆孔點(diǎn)的插值點(diǎn)移除；最后在邊界插值點(diǎn)和內(nèi)部插值點(diǎn)處進(jìn)行插值處理，采用克里金插值算法，利用原始鉆孔勘察資料計(jì)算新增插值點(diǎn)的地層參數(shù)。

2.2 約束Delaunay三角／四面體剖分算法

Delaunay三角網(wǎng)具有優(yōu)良的幾何特性，是公認(rèn)的最優(yōu)三角網(wǎng)。Delaunay三角剖分有兩條重要的性質(zhì)：圓形準(zhǔn)則和最大－最小角準(zhǔn)則，后者保證了Delaunay三角剖分是當(dāng)前點(diǎn)集的最優(yōu)剖分。

在實(shí)際應(yīng)用中，往往要求生成的Delaunay三角網(wǎng)滿足一些約束（限定）條件：在二維情況下有約束點(diǎn)、約束線，在三維情況下還包括約束面。在上述約束條件下生成Delaunay三角網(wǎng)格的過(guò)程稱(chēng)作約束Delaunay三角剖分，也稱(chēng)為限定Delaunay三角剖分。在本文地質(zhì)建模中，地質(zhì)結(jié)構(gòu)中的二維約束條件的輸入格式為平面直線圖（Plane Straight Line Graph，PSLG），此格式為點(diǎn)與線段的集合［11］，如圖1的PSLG范例所示。

圖1 PSLG范例Fig.1 PSLG example

三維約束條件的輸入格式為分段連續(xù)線性組合物（PLC），此格式是節(jié)點(diǎn)、線段與刻面的集合［11］。節(jié)點(diǎn)和線段的定義與PSLG相同，刻面可以是復(fù)雜或不規(guī)則的多邊形，也可能有數(shù)條邊、可能不是凸面狀，刻面上可能有洞、裂縫或節(jié)點(diǎn)，但刻面永遠(yuǎn)是平面。三維空間中任何由平面片、直線段和點(diǎn)構(gòu)成的限定條件經(jīng)過(guò)有限次的運(yùn)算后都可以轉(zhuǎn)換成規(guī)范的PLC限定條件，如圖2PLC范例所示。

圖2 PLC范例Fig.2 PLC example

約束Delaunay三角剖分（CDT），允許在域中加入額外的點(diǎn)從而使得域的所有邊界存在于剖分結(jié)果中，剖分后的所有單元均不包含其他可視節(jié)點(diǎn)，但并不是所有單元均符合Delaunay準(zhǔn)則。由于對(duì)邊界的細(xì)化處理加入了新的節(jié)點(diǎn)，導(dǎo)致初始剖分點(diǎn)集是最后的剖分點(diǎn)集的子集。

CDT的可視性：如圖3a、圖3b所示，粗線為輸入的線段，三角形t與邊e外接圓內(nèi)都有兩個(gè)侵入的節(jié)點(diǎn)，這兩節(jié)點(diǎn)對(duì)t或e而言都是“不可視的”，因?yàn)槎急惠斎刖€段擋住。三維CDT則增加了刻面對(duì)輸入節(jié)點(diǎn)的可視限制。如圖3c所示，粗線為三維的CDT四面體，而此四面體的外接球內(nèi)有一個(gè)節(jié)點(diǎn)存在，但被其陰影的刻面遮蔽了，則是不可視的。

圖3 節(jié)點(diǎn)可視性Fig.3 The visibility of vertex

約束Delaunay三角／四面體剖分新增點(diǎn)原則：1）侵入線段：當(dāng)輸入線段的徑圓／徑球包含其他節(jié)點(diǎn)即被其他節(jié)點(diǎn)侵入時(shí)，則新增此線段的中點(diǎn)。2）侵入刻面：當(dāng)輸入刻面的赤道球被其他節(jié)點(diǎn)侵入時(shí)，則新增此赤道球的球心即刻面的外接圓圓心。但若新增的圓心侵入了某條線段的徑球時(shí)，則刪除此圓心，改為新添此線段的中點(diǎn)。

約束Delaunay三角／四面體剖分方法如下：1）邊界離散。將輸入的PSLG或PLC中的線段和刻面離散成小的線段或者小三角形。2）點(diǎn)集Delaunay三角／四面體剖分。將離散后的所有頂點(diǎn)集合及約束點(diǎn)集合采用Bowyer－Watson算法［16］進(jìn)行點(diǎn)集的Delaunay三角／四面體剖分。3）邊界恢復(fù)。對(duì)每個(gè)小線段和三角形檢查其在三角／四面體網(wǎng)格中的存在性，如果一個(gè)小線段和小三角形在網(wǎng)格中不存在，則利用CDT的可視性對(duì)其進(jìn)行恢復(fù)。4）新增節(jié)點(diǎn)。若邊界線段和三角形被其他節(jié)點(diǎn)侵入則按照對(duì)應(yīng)的新增點(diǎn)原則對(duì)其進(jìn)行細(xì)分以減小線段和三角形的尺度，直到所有小三角形和小線段都不被其他節(jié)點(diǎn)侵入。5）網(wǎng)格質(zhì)量控制。若網(wǎng)格單元超出給定的網(wǎng)格質(zhì)量限制參數(shù)（三角形的最小角或最大面積／四面體的半徑邊長(zhǎng)比或體積），則插入其外接圓／外接球的圓／球心，并遵守Delaunay外接圓／球定理［17］。但是若新添的圓／球心侵入了某小線段或者小三角形，則刪除此圓／球心，按照步驟4添加點(diǎn)。6）刪除邊界以外的網(wǎng)格單元，算法結(jié)束。

2.3 構(gòu)建地質(zhì)體模型算法

三維地質(zhì)建模方法主要有模塊法、線框法、實(shí)體法、表面法4種，本文針對(duì)地下淺層空間提出的以地質(zhì)層面Delaunay三角網(wǎng)為約束條件進(jìn)行約束四面體剖分構(gòu)建地質(zhì)體模型的算法屬于實(shí)體法。該算法在對(duì)原始鉆孔數(shù)據(jù)進(jìn)行預(yù)處理和插值處理后，提取區(qū)域邊界線，刪除邊界外的鉆孔點(diǎn)（包括虛擬鉆孔點(diǎn)），利用剩余的鉆孔點(diǎn)以及邊界構(gòu)建各層面約束條件PSLG，依照前述三角剖分算法進(jìn)行層面Delaunay三角剖分。在層面三角網(wǎng)的基礎(chǔ)上，結(jié)合鉆孔數(shù)據(jù)及其空間分布信息構(gòu)建地質(zhì)體約束條件PLC，根據(jù)設(shè)定的體積或者半徑邊長(zhǎng)比控制網(wǎng)格質(zhì)量，對(duì)PLC進(jìn)行約束Delaunay四面體剖分，最后利用VTK組件庫(kù)對(duì)建模結(jié)果進(jìn)行可視化。

3 實(shí)例分析

采用本文的地質(zhì)體建模方法，以南京市河西地區(qū)為研究區(qū)域，構(gòu)建了該區(qū)域的四面體三維地質(zhì)模型。圖4a為對(duì)鉆孔數(shù)據(jù)進(jìn)行插值處理且刪除區(qū)域外鉆孔點(diǎn)后的鉆孔點(diǎn)及區(qū)域邊界線即為各地質(zhì)層面的約束Delaunay三角剖分的限定條件PSLG。圖4b為地質(zhì)體四面體剖分約束條件PLC，共有5層約束三角網(wǎng)。根據(jù)上述PLC及用戶(hù)指定的網(wǎng)格質(zhì)量控制參數(shù)，運(yùn)用前述算法進(jìn)行地質(zhì)體約束Delaunay四面體剖分，構(gòu)建地質(zhì)體模型，其可視化結(jié)果如圖4c和圖4d所示?？梢钥闯?，本文的算法能夠成功地對(duì)層狀地質(zhì)體進(jìn)行約束Delaunay四面體剖分，并且可以設(shè)置四面體的體積或者半徑變長(zhǎng)比進(jìn)行網(wǎng)格質(zhì)量控制，改善網(wǎng)格單元的質(zhì)量，構(gòu)建良好的三維地質(zhì)體模型。

圖4 構(gòu)建地質(zhì)體模型實(shí)例Fig.4 The example of constructing geological model

算法實(shí)現(xiàn)環(huán)境為微機(jī)環(huán)境，CPU為Pentium IV 2.0GHz，內(nèi)存為2GB，軟件環(huán)境為 Windows XP Professional版，算法使用C＃語(yǔ)言實(shí)現(xiàn)，結(jié)合VTK組件庫(kù)進(jìn)行可視化，在Visual Studio 2008編譯環(huán)境下編譯通過(guò)。

4 結(jié)論

本文考慮到Delaunay三角剖分良好的邊界約束適應(yīng)能力，提出了以地質(zhì)層面Delaunay三角網(wǎng)為約束條件進(jìn)行四面體剖分的地質(zhì)建模算法。該算法針對(duì)地質(zhì)鉆孔數(shù)據(jù)稀疏且分布不均勻的問(wèn)題，采用空間插值算法進(jìn)行了適當(dāng)?shù)膬?nèi)插，利用邊界約束和鉆孔點(diǎn)約束構(gòu)建了各地質(zhì)層面的約束Delaunay三角網(wǎng)，在此基礎(chǔ)上，結(jié)合鉆孔數(shù)據(jù)及其空間分布信息構(gòu)建地質(zhì)體約束，根據(jù)設(shè)定的最大體積或者半徑邊長(zhǎng)比控制網(wǎng)格質(zhì)量，進(jìn)行四面體剖分構(gòu)建三維地質(zhì)實(shí)體模型，借助VTK組件庫(kù)進(jìn)行模型可視化。實(shí)例證明本文提出的算法是可行的，可得到可視化效果良好的四面體地質(zhì)模型，且可高效地支撐地學(xué)空間分析與計(jì)算，為城市規(guī)劃提供更好的決策支持。

［1］ 金江軍，潘懋，曲紅剛，等.三維地質(zhì)建模及其在地下空間開(kāi)發(fā)中的應(yīng)用［J］.信息技術(shù)，2007（3）：26－29.

［2］ 郝海森，吳立新.基于強(qiáng)約束Delaunay－TIN的三維地學(xué)模擬與可視化［J］.地理與地理信息科學(xué)，2003，19（2）：15－18.

［3］ 蔡強(qiáng)，楊欽，陳其明.地質(zhì)結(jié)構(gòu)重疊域的限定Delaunay三角剖分研究［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2004，16（6）：766－771.

［4］ 陳永鋒，陳曉鵬，魏麗英.基于鉆孔數(shù)據(jù)的地質(zhì)面Delaunay三角剖分［J］.金屬礦山，2009（9）：120－122.

［5］ 孟永東，徐衛(wèi)亞，田斌，等.基于帶約束三角剖分的三維地質(zhì)建模方法及應(yīng)用［J］.系統(tǒng)仿真學(xué)報(bào)，2009，21（19）：5985－5989.

［6］ 李培軍.層狀地質(zhì)體的三維模擬與可視化［J］.地學(xué)前緣，2000（7）：271－277.

［7］ 孟憲海，楊欽，李吉?jiǎng)?基于層面結(jié)構(gòu)的三維閉合地質(zhì)區(qū)塊構(gòu)造算法［J］.北京航空航天大學(xué)學(xué)報(bào)，2005，31（2）：182－186.

［8］ 紅雄，張建勛.基于TIN和STP的3D地學(xué)混合構(gòu)模［J］.工程圖學(xué)學(xué)報(bào)，2006（4）：47－71.

［9］ 劉少華，吳東勝，羅小龍，等.復(fù)雜地質(zhì)體的三維建模與可視化方法的研究［J］.礦業(yè)研究與開(kāi)發(fā)，2007，27（2）：56－58.

［10］ 蔡強(qiáng)，李海生，左敏，等.基于Delaunay三角剖分的復(fù)雜地質(zhì)結(jié)構(gòu)建模［J］.金屬礦山，2010（4）：126－130.

［11］ PREPARATA F P，SHAMOS M I.Computational Geometry［M］.New York／Berlin：Springer－Verlag，1985.

［12］ HOULDING S W.3DGeoscience Modelling：Computer Techniques for Geological Characterization［M］.London：Spring－Vetlag，1994.

［13］ 張煜，白世偉.一種基于三棱柱體體元的三維地層建模方法與應(yīng)用［J］.中國(guó)圖象圖形學(xué)報(bào)，2001，6（3）：285－290.

［14］ 齊安文，吳立新.基于類(lèi)三棱柱的三維地質(zhì)模擬與拓?fù)溲芯浚跩］.礦山測(cè)量，2003（3）：65－67.

［15］ 車(chē)德福，陳學(xué)習(xí)，吳立新，等.基于廣義三棱柱體元的三維地層建模方法［J］.遼寧工程技術(shù)大學(xué)學(xué)報(bào)，2006，25（1）：36－38.

［16］ EDELSBRUNNER H.Geometry and Topology of Mesh Generation［M］.England：Cambridge Univ.Press，2001.

［17］ SI H，G?RTNER K.Meshing Piecewise Linear Complexes by Constrained Delaunay Tetrahedralizations［A］.Proceedings of the 14th International Meshing Roundtable［C］.California／San Diego，2005.147－163.