孫黎明 魏迎奇 蔡紅 嚴俊 宋建正 喬蕓蕓

摘要：為處理地質(zhì)界面之間的空間相交關(guān)系，提出一種新的針對三角地質(zhì)曲面的快速求交方法。該方法融合優(yōu)化八叉樹法和OBB搜索樹方法，可以更快速準(zhǔn)確地剔除遠離交線的其他三角形。求交剩余的三角形得到交線，應(yīng)用三角網(wǎng)局部重構(gòu)和網(wǎng)格優(yōu)化算法修正交線附近的三角網(wǎng)，最終分割交線兩側(cè)的地質(zhì)曲面，完成2個地質(zhì)曲面的離散化求交過程。與AABB、OBB和空間分解法相比，該方法在大數(shù)據(jù)量三角曲面求交中效率優(yōu)勢明顯，可以快速準(zhǔn)確處理地質(zhì)模型構(gòu)建和分析中的曲面求交問題，為三維地質(zhì)模型自動化構(gòu)建的實現(xiàn)提供有效支撐。

關(guān)鍵詞：地質(zhì)建模；地質(zhì)曲面；曲面求交；八叉樹；OBB搜索樹；三角網(wǎng)重構(gòu)；模型切割；區(qū)域離散化

中圖分類號：P221.1

文獻標(biāo)志碼：B

0 引 言

三角形的地質(zhì)曲面是構(gòu)成三維地質(zhì)模型的基本幾何元素，在大區(qū)域復(fù)雜地質(zhì)模型的自動化構(gòu)建中，通常需要處理地質(zhì)界面之間的空間相交關(guān)系，這是三維地質(zhì)建模應(yīng)用研究的熱點問題之一。[1-6]在建模過程中，不整合地層邊界計算、斷層面與地層的交線和任意切割面與地層的交線等[7-8]應(yīng)用廣泛。三維地質(zhì)模型的分析應(yīng)用，例如剖面制作、巷道模擬和地下工程等[9]，需要準(zhǔn)確快速得到相交區(qū)域的空間位置、拓撲信息和內(nèi)部構(gòu)造，以對模型區(qū)域有更準(zhǔn)確的理解，同時也需要進行多次曲面求交和計算分析。[10]隨著地質(zhì)建模數(shù)據(jù)量越來越龐大，尤其是地震數(shù)據(jù)重構(gòu)斷層和地層時，三維地質(zhì)建模需要一種快速穩(wěn)定的三角曲面求交算法。地質(zhì)曲面的求交過程根據(jù)2個已有的三角曲面集合，得到相交曲面之間的交線，并重構(gòu)2個三角曲面和交線兩側(cè)的三角形，完成曲面的離散化。如果循環(huán)判斷求交，那么所有的三角形都要參與求交計算，算法的效率為O（n2）。

實際進行求交的三角形對直接影響三角曲面求交速度，因而算法的核心是快速準(zhǔn)確地找到與交線相關(guān)的三角形。目前，曲面求交和物體碰撞檢測的方法主要有空間分解法（八叉樹、BSP樹和KD樹）和層次包圍盒法（AABB、OBB、球體包圍盒和k-DOP）等。地質(zhì)曲面求交通常直接采用物體碰撞檢測的方法，LINDENBECK等[11]開發(fā)TRICUT軟件，用Rapid庫解決三角網(wǎng)求交問題，應(yīng)用廣泛。CGAL計算幾何算法庫改進了AABB樹求交方法 [12]，基于OBB的地質(zhì)曲面切割方法也廣泛應(yīng)用并取得良好效果[13]，但在數(shù)據(jù)量較大時，構(gòu)建層次搜索樹十分耗時。YANG等[8]提出一種基于混合包圍盒的針對真三維模型的切割方法，可解決多種幾何體一體切割的問題。ELSHEIKH等[14]提出一種基于追蹤法的可靠地質(zhì)曲面求交方法，需要利用三角網(wǎng)拓撲追蹤曲面交線，但尋找起始追蹤點較困難。

本文提出一種改進的地質(zhì)曲面快速求交方法，用結(jié)構(gòu)化的八叉樹算法與改進的OBB混合搜索樹相結(jié)合的方法快速求取三角曲面交線，并重構(gòu)交線所影響的三角形，然后將算法應(yīng)用于地質(zhì)建模分析的典型問題中，包括曲面切割地質(zhì)模型，斷層面求交和地質(zhì)區(qū)域離散化等。為將復(fù)雜地質(zhì)曲面求交問題簡化，將地層求交分解為任意兩兩三角曲面求交。

1 理論方法和求交步驟

檢測2個三角曲面的相交部分，核心是快速找到相交的三角形。無論曲面三角形數(shù)量多大，實際相交的三角形都很少，直接構(gòu)建層次搜索樹效率很低，因此快速剔除大部分不相交的三角形，只留下交線附近區(qū)域的三角形，然后再構(gòu)建層次搜索樹，可大大提升求交速度。采用結(jié)構(gòu)化的八叉樹數(shù)據(jù)結(jié)構(gòu)高效剔除求交曲面中的大多數(shù)不相交三角形，構(gòu)建基于改進OBB的搜索樹，快速準(zhǔn)確地求取交線，改進曲面求交方法流程見圖1。采用靜態(tài)結(jié)構(gòu)化數(shù)據(jù)結(jié)構(gòu)[15-16]表達1個三角曲面，即用1個大于三角網(wǎng)的六面體等分后表達地質(zhì)界面和2個曲面相交的三角形，建立特殊的層次二叉搜索樹，根節(jié)點根據(jù)三角形的AABB表示，葉節(jié)點則用OBB代替，最后對三角形求交得到交線。

1.1 優(yōu)化八叉樹的粗略求交

空間八叉樹是二維平面中四叉樹在三維空間的擴展，是用均勻細分空間描述三維場景中幾何對象的一種樹狀數(shù)據(jù)結(jié)構(gòu)，其每個父節(jié)點用一個正方體表示，每個父節(jié)點都有8個子節(jié)點，是將父節(jié)點均勻細分為8個小正方體得到的。八叉樹通常用于三維空間中海量物體的快速查詢和空間索引，目前廣泛應(yīng)用于大數(shù)據(jù)量的三維GIS，包括三維智慧城市、三維數(shù)字管線、三維BIM可視化和3D大場景游戲等。

八叉樹的特點是可以用多層樹快速索引每個節(jié)點，但對于海量數(shù)據(jù)離散點構(gòu)建的地質(zhì)曲面和細分到單一節(jié)點的八叉樹，每次計算的迭代數(shù)量龐大，且八叉樹的層級非常多，計算效率更低，每次循環(huán)計算耗時很長。本文提出一種非對稱的空間八叉樹模型，用于快速查找2個三角曲面交線所在的三角形。與傳統(tǒng)空間八叉樹模型不同，該模型不是每層都按照8n細分，而是達到一定的邊界條件就采用結(jié)合OBB的方法進行求交。設(shè)計八叉樹的邊界條件，計算八叉樹停止的條件，以省略八叉樹細分到葉節(jié)點的過程，提高求交的計算效率。優(yōu)化八叉樹求交過程見圖2。

地質(zhì)曲面所在的三維笛卡爾坐標(biāo)系定義為R，定義三角網(wǎng)中的點VR、三角形TR、地層三角網(wǎng)SR。任意點P（x， y， z）對應(yīng)立方體ER，立方體的邊長為三角網(wǎng)包圍盒在x、y和z這3個方向上長度的最大值l。依次將每個父節(jié)點分為8個子節(jié)點，葉節(jié)點滿足與同行列子節(jié)點內(nèi)的曲面點個數(shù)接近條件，并且當(dāng)個數(shù)最小時八叉樹就不再細分。與交線不相關(guān)的不相交三角形見圖3，剔除遠離交線后的三角網(wǎng)見圖4。

利用八叉樹結(jié)構(gòu)計算1個子節(jié)點內(nèi)的三角形，可以很快找到距離另一個曲面最近的三角形集合（見圖3）。顯然，該方法得到的結(jié)果仍然包括許多不相交的三角形，如果單用此算法求交，效率不高，因此本文結(jié)合改進OBB層次搜索樹的方法進一步剔除不相交的三角形。

1.2 OBB混合搜索樹的精確求交

OBB是1個緊密的包圍盒，是包含該物體且相對于坐標(biāo)軸任意方向的最小六面體，可提高檢測效率。2個相交地質(zhì)曲面的OBB包圍盒見圖5。

構(gòu)建OBB搜索樹和計算OBB比計算AABB更繁瑣耗時。根據(jù)幾何坐標(biāo)計算物體的OBB，因為三角形OBB需要多2個向量計算，所以大數(shù)據(jù)量時直接構(gòu)建OBB搜索樹復(fù)雜低效。GOTTSCHALK等[17]提出一種有效的實現(xiàn)方法，采用AABB構(gòu)建中間所有層次的根節(jié)點，將第1.1節(jié)中計算的六面體包圍盒直接作為AABB，而葉節(jié)點采用OBB方法，OBB混合搜索樹見圖6。

搜索樹采用二叉樹結(jié)構(gòu)存

儲，最終葉節(jié)點采用OBB求交，能最大限度減少精確判斷相交的三角形對。綜合利用2種包圍盒的優(yōu)點，可減少計算量，加快搜索樹的構(gòu)建速度。[18]

根據(jù)可能相交的三角形構(gòu)建層次搜索樹的簡要步驟如下：

（1）選擇根節(jié)點。將所有的三角形最大包圍盒投影到3個坐標(biāo)軸，找到3個方向中的最長軸，選擇其中間值為最大值與最小值之和的1/2。

（2）構(gòu)建搜索樹。根據(jù)平均三角形AABB的大小，選擇樹的深度層次。樹的深度選擇是影響搜索求交速度的重要因素。經(jīng)過第1.1節(jié)剔除工作后，求交的三角形只占曲面三角形總數(shù)的小部分，層次不需太深。

（3）遞歸分割。同時將每個三角形和點的包圍盒自頂向下遍歷并存放指針到每一層合適的位置節(jié)點，直到單一三角形為止，向下一直到葉節(jié)點都是三角形，葉節(jié)點用OBB代替AABB。

（4）求交。自頂而下判斷包圍盒是否相交，如果遞歸到葉節(jié)點，根據(jù)三角形的3個點構(gòu)造OBB，采用SAT方法求交。[19-20]

1.3 交線區(qū)域網(wǎng)格重構(gòu)

對得到的所有三角形對逐一求交，即可得到最終的交線。一對空間三角形相交的結(jié)果有3種情況：點、線段和相交三角形。地質(zhì)曲面求交中最常見的是線段，線段首尾相連形成交線。交線可以有多條，也可能是封閉環(huán)。目前，空間三角形求交的算法主要分為標(biāo)量判別法和矢量判別法[21]，本文選擇適于地質(zhì)曲面求交的MLLER[22]三角形求交法。

加入新的交線后，將交線點分別加入到交線兩側(cè)的三角網(wǎng)中，旋轉(zhuǎn)變換（不是投影）到合適的二維平面。[23]采用局部三角形修復(fù)方法，若交點在三角形的邊上，則在此處添加新點，將這條邊分割；若交點在三角形內(nèi)部，則在內(nèi)部添加新點，連接頂點構(gòu)成新三角形，盡量滿足Delaunay法則，加入新點后的三角網(wǎng)重構(gòu)見圖7。

新的三角網(wǎng)交線區(qū)域必須進行網(wǎng)格質(zhì)量優(yōu)化。若新加入的點形成狹長三角形和拓撲錯誤三角形，在地質(zhì)模型空間分析、體元劃分和數(shù)值評估中都無法滿足要求，則需要進行三角網(wǎng)質(zhì)量優(yōu)化，局部修改

重構(gòu)。根據(jù)交線點所在的三角形的網(wǎng)格拓撲找到周圍三角形[24-25]，優(yōu)化網(wǎng)格質(zhì)量，包括邊界交換、合并三角形、點重構(gòu)和插入新點。[26]

將局部修復(fù)好的三角形重新變換回三維空間，將原來的曲面以交線為邊界分為2個曲面，判斷其他三角形頂點相對交線的位置，若分別屬于不同的曲面，則求取2個曲面的交線。多個地質(zhì)曲面求交見圖8。

2 實驗與應(yīng)用

2.1 算法性能測試對比

在一臺型號為Intel i3-2330M、CPU為2.1 GHz、內(nèi)存為6 GB的計算機上測試算法效率，并與其他2種方法進行對比，見圖9。靜態(tài)結(jié)構(gòu)化數(shù)據(jù)搜索的速度與分割的包圍盒的個數(shù)有關(guān)，與三角形個數(shù)關(guān)系不大，時間復(fù)雜度為O（2n），保證用于構(gòu)建OBB混合搜索樹的三角形個數(shù)隨求交三角形個數(shù)的增加而平穩(wěn)增長。該改進算法在數(shù)據(jù)量大的曲面求交中效率優(yōu)勢明顯。

2.2 改進算法在三維地質(zhì)建模中的應(yīng)用

由于地質(zhì)構(gòu)造的復(fù)雜性，地質(zhì)曲面之間的空間關(guān)系非常復(fù)雜，因此在建模過程中，地層與地層或斷層相交時，通常需要求取交線并離散化，最常用的方法是任意剖面切割法。三維曲面求交方法為基礎(chǔ)算法，可以很好地應(yīng)用到地質(zhì)自動化建模和分析中。

在三維地質(zhì)建模中，構(gòu)建含斷層的地質(zhì)模型是研究熱點。根據(jù)地質(zhì)剖面圖確定斷點、斷距和傾角等基本信息，構(gòu)建初始斷層面，得到現(xiàn)有斷層面與地層面的求交線，切割上、下2個地層面，重構(gòu)地層和斷層面邊界。保證模型和數(shù)據(jù)的一致性，檢查拓撲關(guān)系的正確性是核心。含簡單斷層的地層模型見圖10。

復(fù)雜的地質(zhì)構(gòu)造有許多界面相交，分割線通過一系列線段表示，不同區(qū)域之間用多個線段隔開。為自動化模型重構(gòu)和模型數(shù)值分析[27]，需要分離交線分割的不同區(qū)域，手動方法繁瑣且不可靠。為實現(xiàn)自動分離不同區(qū)域，嘗試改進前文提出的求交算法以處理簡單區(qū)域劃分的問題。先確定交線類型，再根據(jù)交線追蹤區(qū)域邊界，形成內(nèi)部環(huán)。遇到邊界就停止追蹤，形成多個相鄰多邊形封閉區(qū)域。根據(jù)三角形和交線的相對位置，劃分所有三角形到其所屬的區(qū)域內(nèi)，完成地質(zhì)區(qū)域離散，見圖11。

3 結(jié)束語

為解決三維地質(zhì)建模過程和模型分析中的三角曲面求交問題，提出基于優(yōu)化八叉樹的搜索算法和OBB混合搜索樹相結(jié)合的曲面求交方法。通過結(jié)構(gòu)化搜索快速剔除遠離交線的不相交三角形，結(jié)合OBB混合層次搜索樹，快速求取三角曲面的交線，并將該方法應(yīng)用到地質(zhì)建模和分析中，驗證算法的可靠性和實用性。采用優(yōu)化的八叉樹快速剔除大量不相交三角形，可提升求交算法的效率，除三角曲面外還可以擴展到其他類型曲面。將該算法應(yīng)用到地質(zhì)分析中，解決參數(shù)曲面切割地質(zhì)模型、斷層構(gòu)建和區(qū)域離散化問題。將該算法作為基礎(chǔ)算法工具應(yīng)用到三維地質(zhì)模型的構(gòu)建和分析中，為快速自動化建模分析提供一種新的方法。

改進的曲面求交方法可以較好地處理地質(zhì)曲面求交的相關(guān)問題，未來可以嘗試采用GPU并行計算的方法提升運算效率。在解決地質(zhì)曲面求交相關(guān)應(yīng)用問題上，工作流程自動化和減少人工干預(yù)是未來努力的主要方向。

參考文獻：

[1] 鄧飛， 王瑞， 王美平， 等. 復(fù)雜三維地層建模及快速射線追蹤的研究與實現(xiàn)[J]. 大慶石油地質(zhì)與開發(fā)， 2007， 26（1）： 113-118. DOI： 10.3969/j.issn.1000-3754.2007.01.028.

[2] LI N， BAGAS L， LINDSAY M， et al. An irregular triangle mesh buffer analysis method for boundary representation geological object in three-dimension[J]. Earth Science Informatics， 2017， 10（2）： 149-167.

[3] FU Q， WU Z， WANG X， et al. An algorithm for finding intersection between ball B-spline curves[J]. Journal of Computational and Applied Mathematics， 2018， 327： 260-273.

[4] WELLMANN J F， de la VARGA M， MURDIE R E， et al. Uncertainty estimation for a geological model of Sandstone greenstone belt， Western Australia-insights from integrated geological and geophysical inversion in a Bayesian inference framework[J]. Geological Society， London， Special Publications， 2017， 453（1）： 41. DOI： 10.1144/SP453.12.

[5] CALCAGNO P， COURRIOUX G， LOPEZ S， et al. How geological architecture helps 3D modelling[C]// Proceedings of 4th Meeting of European 3D GeoModelling Community. Orleans， 2018.

[6] JIANG Q P， JIE T， YUAN C F. Fast triangle mesh surface intersection algorithm based on uniform grid[J]. Computer Engineering， 2008， 34（21）： 172-174.

[7] BISTACCHI A， MASSIRONI M， dal PIAZ G V， et al. 3D fold and fault reconstruction with an uncertainty model： An example from an Alpine tunnel case study[J]. Computers & Geosciences， 2008， 34（4）： 351-372.

[8] YANG Y， LIU X J， ZHOU B， et al. Research of 3D GIS section for true three-dimensional geo-spatial model[C]// Proceedings of 2010 International Conference on Multimedia Technology. Ningbo： IEEE， 2010： 1-4. DOI： 10.1109/ICMULT.2010.5631185.

[9] 孫黎明， 李青元， 譚海， 等. 基于Netgen的層狀地質(zhì)體四面體網(wǎng)格劃分方法[J]. 計算機輔助工程， 2013， 22（3）： 47-52. DOI： 10.3969/j.issn.1006-0871.2013.03.010.

[10] FRANK T， TERTOIS A L， MALLET J L. 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[J]. Computers & Geosciences， 2007， 33（7）： 932-943. DOI： 10.1016/j.cageo.2006.11.014.

[11] LINDENBECK C H， EBERT H D， ULMER H， et al. TRICUT： A program to clip triangle meshes using rapid and triangle libraries and visualization toolkit[J]. Computers & Geosciences， 2002， 28（7）： 841-850. DOI： 10.1016/S0098-3004（01）00110-8.

[12] FABRI A， PION S. CGAL-computational geometry algorithms library[C]// Proceedings of 17th ACM SIGSPATIAL International Symposium on Advances in Geographic Information Systems. Washington， 2009. DOI： 10.1145/1653771.1653865.

[13] 黃松柏， 徐華. 基于動態(tài)OBB層次結(jié)構(gòu)的曲面相交算法[J]. 計算機應(yīng)用研究， 2011， 28（8）： 3181-3184. DOI： 10.3969/j.issn.1001-3695.2011.08.106.

[14] ELSHEIKH A H， ELSHEIKH M. A reliable triangular mesh intersection algorithm and its application in geological modelling[J]. Engineering with Computers， 2014， 30（1）： 143-157. DOI： 10.1007/s00366-012-0297-3.

[15] SAMET H. Foundations of multidimensional and metric data structures[M]. San Francisco： Morgan Kaufmann， 2006.

[16] DASSI F， PEROTTO S， FORMAGGIA L， et al. Efficient geometric reconstruction of complex geological structures[J]. Mathematics and Computers in Simulation， 2014， 106： 163-184. DOI： 10.1016/j.matcom.2014.01.005.

[17] GOTTSCHALK S， LIN M C， MANOCHA D. OBB tree： A hierarchical structure for rapid interference detection[C]// Proceedings of 23rd Annual Conference on Computer Graphics and Interactive Techniques. New York， 1996： 171-180. DOI： 10.1145/237170.237244.

[18] FANG Z G， XU J， JIANG J X， et al. Efficient collision detection using bounding volume hierarchies of OBB-AABBs and its application[C]//Proceedings of 2010 International Conference on Computer Design and Applications. Qinhuangdao： IEEE， 2010. DOI： 10.1109/ICCDA.2010.5541315.

[19] CHANG J W， KIM M S. Efficient triangle-triangle intersection test for OBB-based collision detection[J]. Computers & Graphics， 2009， 33（3）： 235-240. DOI： 10.1016/j.cag.2009.03.009.

[20] VAN DEN BERGEN G. Efficient collision detection of complex deformable models using AABB trees[J]. Journal of Graphics Tools， 2012，2（4）： 37-41. DOI： 10.1080/10867651.1997.10487480.

[21] 鄒益勝， 丁國富， 何邕， 等. 空間三角形快速相交檢測算法[J]. 計算機應(yīng)用研究， 2008， 25（10）： 2906-2910. DOI： 10.3969/j.issn.1001-3695.2008.10.007.

[22] MLLER T. A fast triangle-triangle intersection test[J]. Journal of Graphics Tools， 1997， 2（2）： 25-30. DOI： 10.1080/10867651.1997.10487472.

[23] MCLAURIN D， MARCUM D， REMOTIGUE M， et al. Repairing unstructured triangular mesh intersections[J]. International Journal for Numerical Methods in Engineering， 2013， 93（3）： 266-275. DOI： 10.1002/nme.4385.

[24] HOPPE H， DEROSE T， DUCHAMP T， et al. Mesh optimization[C]// Proceedings of 20th Annual Conference on Computer Graphics and Interactive Techniques. New York， 1993. DOI： 10.1145/166117.166119.

[25] ATTENE M， FALCIDIENO B. Remesh： An interactive environment to edit and repair triangle meshes[C]// Proceedings of IEEE International Conference on Shape Modeling and Applications. Matsushima： IEEE， 2006. DOI： 10.1109/SMI.2006.29.

[26] WANG D， HASSAN O， MORGAN K， et al. Enhanced remeshing from STL files with applications to surface grid generation[J]. Communications in Numerical Methods in Engineering， 2007， 23（2）： 227-239. DOI： 10.1002/cnm.894.

[27] 季順迎， 趙金鳳， 狄少丞， 等. 面向環(huán)境力學(xué)的離散元分析軟件研發(fā)和工程應(yīng)用[J]. 計算機輔助工程， 2014， 23（1）： 69-75. DOI： 10.13340/j.cae.2014.01.014.

（編輯 付宇靚）