三維非拓撲一致性塊體構建算法優化與應用

羅開云 房 鵬 何永清 蔣先藝 吳 恒

(①東方地球物理公司采集技術中心,河北涿州 072751； ②成都理工大學信息科學與技術學院，四川成都 610059)

0 引言

地震數值正演模擬廣泛應用于基于模型的射線追蹤、覆蓋次數統計和觀測系統照明等方面[1]。然而三維復雜地質模型的建立是正演的關鍵步驟，復雜模型的塊體構建難度大、耗時長，特別是在采集設計領域，由于工期短、時間緊，設計人員迫切希望在短時間內構建符合設計區域實際地質情況的三維模型。

目前采集設計領域中正演模擬經常使用的射線和射線束方法[2]，需要設計三維塊狀模型，三維塊狀模型的構建是以具有近似地質屬性的層塊作為塊體構建的基本單元。首先利用曲面構建方法建立地質層位和斷層曲面，再利用塊體追蹤技術提取由層位和斷層曲面片所形成的封閉地質塊體[3]。塊體構建所得到的塊狀地質模型通常可以用于射線、射線束類正演模擬[4]，也可以通過離散化算法轉換為規則或者非規則網格模型，支持有限差分或有限元類方法[5]。

三維塊體構建算法的歷史可以追溯到20世紀80年代，Warburton[6]提出一種基于無限大非連續面識別塊體的算法； Lin等[7]提出將拓撲概念引入塊體構建中； Lu[8]將方向邊和方向環的概念引入塊體構建過程； 蔣先藝等[9]提出一種封閉模型描述方法，將具有相同屬性的地質層面聯系在一起，并且給出了一種基于三角形的塊體追蹤算法，從而得出以地質屬性為基礎的塊體； 侯衛生等[10]通過對Lu等[7-8]提出的方法進行改進，可以在曲面上進行搜索得出塊體模型； 楊洋等[11]簡化并改進了傳統的線框架模型，提出了線框架拓撲模型，加速并簡化了拓撲搜索過程。

目前三維地質塊體構建算法已經相當成熟，在塊體追蹤方面從簡單的逐三角形追蹤發展為基于線框拓撲的追蹤，執行效率和容錯性都有很大提高[12]。但本質上，這些方法都是基于模型交線拓撲結構進行塊體追蹤，對準確的地質數據拓撲結構依賴性較強[13]。因而，在進行塊體構建之前需進行多次曲面求交或裁剪等操作，以保證地質模型拓撲關系的正確性。較復雜地質模型存在三角網的尺寸差異過大、密度不均勻、浮點誤差等現象，會導致曲面間的拓撲關系存在少量錯誤，使傳統的拓撲一致性塊體構建方法失效。然而少量的拓撲錯誤并不會影響模型整體的正確性，如果能構建出結構正確的塊體模型，則少量的拓撲錯誤和塊體局部漏縫并不會對射線或模型離散化造成明顯的影響，可以繼續用于正演數值模擬[14]。

本文提出一種非拓撲一致的三維地質塊體構建技術。該方法模擬人眼視覺觀察，完全不依賴交線拓撲信息，僅依靠曲面自身形態提取塊體三角網，在不對原曲面進行求交或裁剪等操作的情況下構建三維地質塊體。

1 非拓撲一致性塊體構建算法原理

以交線拓撲結構為基礎的傳統塊體構建算法根據三角網鄰接關系建立模型的拓撲關系圖，并在此基礎上進行廣度或深度遍歷算法進行塊體外邊界三角網的追蹤[15]。因此拓撲關系的正確性決定著此方法的成功率。然而在面對有少量拓撲問題的模型時，通常會從整體把握理解模型結構，忽略局部拓撲錯誤。基于此思路，本文提出非拓撲一致性塊體構建算法，模擬人類觀察和處理模型的方式，即首先從整體上確定塊體的大致范圍，再利用局部觀察的方式確定塊體的邊界三角網(圖1)。

算法首先讀取模型的曲面三角網，確定模型空間范圍。如圖2a所示，假如模型尺寸為7500m(x)×6400m(y)×7400m(z)， 包含3條斷層、 15個層位曲面。

圖1 非拓撲一致塊體構建算法流程

然后對當前的地質模型進行網格劃分。如圖2b所示，以任一不包含三角形的網格單元作為起始種子(圖2b中紅點所標明的網格單元)，使用漫水填充法[16]對網格單元模型進行邊界網格單元追蹤，得到近似表示當前地質塊體的網格塊體——包含塊體界面三角形的網格單元集合(圖2c)。

最后以基于視覺觀察方式的三角形選取算法逐個對求得的網格塊體選取對應外表面的三角形，如圖2d所示。

非拓撲一致性塊體構建算法首先要確定塊體的基本形態。為此本文利用漫水填充法追蹤塊體所包含的網格單元集合。漫水法是圖形學中常用的一種區域填充法，實質上是一種廣度優先遍歷的種子填充法，適用于對內定義區域的填充。圖形學中的內定義區域，是指區域內部的像素具有類似或相同的屬性，而區域外的所有像素具有不同或差別較大的屬性。該算法從內區域中的某個像素出發(通常稱起始像素為種子)，利用廣度優先搜索完成內區域填充。本文對該方法進行改造，用于三維空間中地質塊體網格單元集合的追蹤，具體步驟如下：

(1)從內存中取得起始種子網格單元，加入到遍歷隊列；

(2)若遍歷隊列為空，則返回步驟(1)，否則從遍歷隊列中取出一個網格單元進行步驟(3)；

(3)若步驟(2)得到的網格單元已被訪問過，則重新進行步驟(2)，否則進行步驟(4)；

(4)找到所述網格單元周圍的六個鄰接網格單元，對每一個鄰接網格單元，更新所述鄰接網格體的投影方向；

(5)對所述鄰接網格單元，若其中不含三角形，則將其加入到遍歷列表中，并將其標記為“已訪問”，并返回步驟(2)； 否則，將其加入到當前塊體的邊界區域網格塊體隊列中，并進行步驟(2)，直至遍歷所有網格單元。

通過漫水追蹤算法獲取了地質塊體的網格單元集合，為精確提取地質塊體外表面三角形做好了準備。非拓撲一致的塊體構建算法的最終目標是提取地質塊體的外邊界三角形，構造地質塊體表面三角網。為此需要遍歷塊體的網格單元，從含有三角形的邊界網格單元中找出屬于地質塊體的三角形，并將它們組合起來，構造塊體外表面三角網。在具體實現上，為達到“篩選出被觀察到的三角形”這一目的，算法使用了圖形庫OpenGL對各個網格體進行渲染繪制，選取三角形單元，并采用平行投影的方式繪制渲染網格體相關三角形[17]。將渲染繪制的結果輸出到OpenGL的幀緩存對象(FBO)中。FBO是OpenGL內置的對象之一，利用該機制，編碼人員可以直接將繪制結果輸出到內存中而不是屏幕上[18]。在本文算法中，這一機制恰好符合從繪制結果直接篩選三角形的需求。

圖2 非拓撲一致塊體構建示意圖(a)輸入曲面模型； (b)對模型空間進行網格劃分； (c)漫水法追蹤塊體外邊界網格單元集合； (d)視覺觀察得到的三角網塊體模型

2 非拓撲一致塊體構建算法優化

在上述非拓撲一致塊體構建算法中，網格劃分是重要步驟。基本的網格劃分方式是均勻型網格劃分，其劃分方式為無區別地將整個地質模型均勻劃分為規模統一、精度一致的網格單元。這種劃分方式簡單，便于實現，且索引效率較高。其缺陷是既耗費資源又耗費時間。為此，本文提出一種高效方法，即非規則八叉樹網格劃分方式。該方法將地質模型劃分為大小不一的網格單元，即對含有曲面三角網的區域劃分為規模更小、精度更高的網格單元，對不含有三角網的區域不再繼續劃分。由于使用了精度可變的網格體，所需求的內存空間較靈活，可根據場景需求調整[19]。圖3對比了均勻網格體劃分與不規則八叉樹網格體劃分。

圖3 均勻網格體劃分(左)與不規則八叉樹網格體劃分(右)對比

八叉樹網格剖分將立體空間從三個維度上同時進行二分，即將原空間劃分為8個幾何相似的子空間，之后再根據實際情況決定是否需要對某個或某幾個子空間進一步劃分。同時，任何一個被劃分出的子空間都將作為一個節點被保存到八叉樹數據結構中，表示整個原始地質空間的數據結構對象將作為八叉樹的根節點，而沒有被繼續劃分子空間的數據結構對象將作為八叉樹的葉子節點。為了實現上下級節點的索引，表示節點的數據結構應該保留指向上下級數據結構的指針以及相對的位置關系。圖4為八叉樹網格剖分算法流程圖。

通過遞歸方式建立不規則八叉樹結構，首先從輸入的根節點開始，檢查當節點的劃分層數是否小于最大劃分層數，并判斷節點包含的三角形數目是否大于1。如果滿足條件則從x、y、z三個維度上同時進行二分，將節點劃分成8個子節點，為每個子節點設定空間范圍，并根據空間范圍置入所包含的三角形。然后依次對8個子節點進行前述的檢查和遞歸劃分。最終遞歸調用結束后，不規則八叉樹即建立完成。

圖4 八叉樹網格劃分算法流程

由于使用八叉樹網格剖分將地質空間劃分為大小不同、精度不一的網格體，相比均勻型網格塊體，其搜索遍歷方法更復雜。因此本文使用一種基于局部坐標系機制的溢出方向搜索算法[20]解決尋找指定方向上的鄰接網格體問題，完善了對鄰接網格體搜索算法。改進后的算法分為三個階段。第一階段：搜索目標網格體可能存在的上級網格體，記錄由上級網格體到目標網格體的搜索路徑；第二階段：通過第一階段所給出的搜索路徑進行目標網格體的搜索；第三階段：對目標網格體進行嘗試性的細分，以保證作為算法結果的網格體一定是地質空間網格體八叉樹的葉子節點。

為了對比八叉樹與非八叉樹方式的存儲效率，對圖2a所示的曲面模型進行實驗，結果如表1所示。若模型包含薄層或透鏡體等地質構造，為了精細體現這些構造的細節特征，不規則八叉樹應具有較深的劃分層數。當劃分層數不大于5時，不規則八叉樹與規則網格劃分規模一致。由表1中網格規模和劃分時間可以看出，隨著劃分層數的增加，規則網格劃分方式的網格單元數目以8倍的速率增長，而八叉樹劃分方式的網格單元數目和劃分耗時增長幅度遠小于8，并且呈遞減趨勢趨于2。因此，使用不規則八叉樹方式對地質空間進行網格劃分極大地節省了存儲空間，提高了塊體追蹤效率。

表1 不規則八叉樹劃分與規則網格劃分的網格單元數量對比

3 效果展示與對比

為了驗證非拓撲一致性算法的有效性，選擇了一個大小為32km(x)×18km(y)×15km(z)的較復雜地質模型(圖5)進行算法測試。模型含有44個層位/斷層曲面，存在斷裂、尖滅、透鏡體等多種地質現象。為了驗證算法處理非拓撲一致性模型的能力，在圖中A、B兩處人為地制造了一些曲面漏縫(圖6左)和交越(圖6右)，可以看到模型存在比較明顯的拓撲問題。

圖7是利用傳統曲面縫合性方法建立的塊狀模型。可以看出，當模型存在拓撲問題時，曲面縫合性方法得到的塊體模型結構存在明顯問題，即局部存在拓撲錯誤的塊體粘連在一起，得到的塊體與原始模型追蹤得到塊體的結構和數目存在差異。

圖8是用網格塊體構建方法建立的塊狀模型，圖9是網格塊體構建算法細節。可以看出，網格塊體構建算法的整體構建效果較好，但對存在拓撲問題的模型，塊體模型沒有保留這些細節。

圖5 非拓撲一致模型

圖6 圖5模型中細節A(左)和B(右)

圖10是利用本文非拓撲一致性算法建立的塊狀模型。塊體構建采用Intel i7 4790K、主頻為4.0GHz的電腦，耗時23s。

圖7 傳統曲面縫合性算法對圖5模型的追蹤結果

圖8 網格塊體構建算法對圖5模型的追蹤結果

圖9 網格塊體構建算法對圖5模型的追蹤結果細節

圖10 本文算法對圖5模型的追蹤結果(左)及細節(右)

對比圖10與圖8可以看出，當模型存在一定拓撲問題時，本文非拓撲一致塊體構建算法得到的塊狀模型整體結構正確，塊體結構和數目與圖8一致。圖10右給出了非拓撲一致性算法在圖中A、B兩處得到的塊體模型細節。可見在原始模型在交線處的漏縫，在本文非拓撲一致性算法建立的塊狀模型中也呈現出一定的缺漏現象。

實驗結果說明非拓撲一致性算法具有很好的容錯能力，可以較好地處理存在一定拓撲問題的模型。

為更好地測試基于視覺觀察的非拓撲一致塊體構建算法的適用性，利用不同形狀和復雜程度的曲面模型模擬真實的地層，如表2所示。與基于三角形追蹤的封閉模型描述方法，即曲面縫合型塊體構建算法[6]與線框塊體構建法[8]作對比。塊體構建數據和效果對比情況見表3和表4。實驗使用的計算機操作系統為Windows10專業版x64、處理器為Intel(R) Core(TM) i7-6800K@3.4GHz(C6T12)、內存為DDR416GB3000MHz。

表2 塊體構建測試模型

若某次塊體構建的結果中出現本應被分離的兩個塊體黏合到一起，則本次塊體構建結果是失敗的，此時不記錄時間和空間資源的消耗情況。

由表4可見，本文基于視覺觀察效果的塊體構建算法雖然在時間資源和空間資源消耗上都數倍于曲面縫合型塊體構建算法和線框塊體構建算法，但就算法本身而言，對于表2中最復雜的地質模型b62，塊體構建的空間占用不超過2GB，算法執行時間不超過25s，因此一般配置的計算機都可滿足其時間和空間需求。該實驗證明，即使在地質模型存在交線等拓撲結構的情況下，亦可得到較客觀的地質情況。換言之，文中算法具有較強的魯棒性，且資源消耗在可接受范圍內，是一種實用且可靠的地質塊體構建算法。

表3 五個模型塊體構建算法效果對比

注：算法1指曲面縫合型塊體構建算法，算法2指線框塊體構建算法，算法3指基于視覺觀察效果的非拓撲一致塊體構建算法，表4同。

表4 五個模型不同算法時間空間資源消耗對比

注：算法名后括號中的數字為指定的網格體劃分層數，如“算法3(5)”即為使用算法3，并指定網格劃分層數為5的條件下進行地質塊體構建。

4 結論

三維塊體構建技術是構造塊狀地質模型的核心技術，在此過程中，要保證復雜模型的拓撲一致性非常困難。傳統的逐三角形追蹤和線框追蹤方法在模型存在極少量局部拓撲錯誤時能夠建立比較正確的塊狀模型，但對交線拓撲結構依賴性強，容錯性和塊體構建成功率較低。本文提出的非拓撲一致性塊體構建算法，完全不依賴模型交線拓撲結構，根據曲面模型自身形態即可構建，極大地提升了塊體構建成功率。然而相比于傳統追蹤算法，該算法效率相對較低，為此應用基于八叉樹結構提高了算法效率。通過數個典型三維模型的塊體構建實驗，可以得出以下結論：

(1)本文基于視覺觀察方式的三角形選取算法提出的非拓撲一致性塊體構建算法，不依賴于交線拓撲結構，算法容錯性較高，能夠在很大程度上解決復雜模型拓撲不一致的問題。

(2)非拓撲一致塊體構建算法以網格塊體追蹤為基礎，追蹤效率較低，且消耗資源較大，使用八叉樹結構對算法進行優化，可以極大地降低內存消耗，一定程度上提高了追蹤效率。

(3)本文提出基于八叉樹優化的非拓撲一致性塊體構建算法，能夠有效解決拓撲不一致模型的塊體構建問題，具有較高的容錯性和塊體構建成功率，并且算法資源消耗在目前主流硬件平臺上即可實現。由于該算法降低了對曲面模型的要求，因而也降低了三維復雜地質模型構建的整體難度，提高了塊體構建的效率，比傳統方法更具實用性。