三維非拓?fù)湟恢滦詨K體構(gòu)建算法優(yōu)化與應(yīng)用

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

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

0 引言

地震數(shù)值正演模擬廣泛應(yīng)用于基于模型的射線追蹤、覆蓋次數(shù)統(tǒng)計和觀測系統(tǒng)照明等方面[1]。然而三維復(fù)雜地質(zhì)模型的建立是正演的關(guān)鍵步驟，復(fù)雜模型的塊體構(gòu)建難度大、耗時長，特別是在采集設(shè)計領(lǐng)域，由于工期短、時間緊，設(shè)計人員迫切希望在短時間內(nèi)構(gòu)建符合設(shè)計區(qū)域?qū)嶋H地質(zhì)情況的三維模型。

目前采集設(shè)計領(lǐng)域中正演模擬經(jīng)常使用的射線和射線束方法[2]，需要設(shè)計三維塊狀模型，三維塊狀模型的構(gòu)建是以具有近似地質(zhì)屬性的層塊作為塊體構(gòu)建的基本單元。首先利用曲面構(gòu)建方法建立地質(zhì)層位和斷層曲面，再利用塊體追蹤技術(shù)提取由層位和斷層曲面片所形成的封閉地質(zhì)塊體[3]。塊體構(gòu)建所得到的塊狀地質(zhì)模型通?？梢杂糜谏渚€、射線束類正演模擬[4]，也可以通過離散化算法轉(zhuǎn)換為規(guī)則或者非規(guī)則網(wǎng)格模型，支持有限差分或有限元類方法[5]。

三維塊體構(gòu)建算法的歷史可以追溯到20世紀(jì)80年代，Warburton[6]提出一種基于無限大非連續(xù)面識別塊體的算法； Lin等[7]提出將拓?fù)涓拍钜雺K體構(gòu)建中； Lu[8]將方向邊和方向環(huán)的概念引入塊體構(gòu)建過程； 蔣先藝等[9]提出一種封閉模型描述方法，將具有相同屬性的地質(zhì)層面聯(lián)系在一起，并且給出了一種基于三角形的塊體追蹤算法，從而得出以地質(zhì)屬性為基礎(chǔ)的塊體； 侯衛(wèi)生等[10]通過對Lu等[7-8]提出的方法進行改進，可以在曲面上進行搜索得出塊體模型； 楊洋等[11]簡化并改進了傳統(tǒng)的線框架模型，提出了線框架拓?fù)淠Ｐ?，加速并簡化了拓?fù)渌阉鬟^程。

目前三維地質(zhì)塊體構(gòu)建算法已經(jīng)相當(dāng)成熟，在塊體追蹤方面從簡單的逐三角形追蹤發(fā)展為基于線框拓?fù)涞淖粉?，?zhí)行效率和容錯性都有很大提高[12]。但本質(zhì)上，這些方法都是基于模型交線拓?fù)浣Y(jié)構(gòu)進行塊體追蹤，對準(zhǔn)確的地質(zhì)數(shù)據(jù)拓?fù)浣Y(jié)構(gòu)依賴性較強[13]。因而，在進行塊體構(gòu)建之前需進行多次曲面求交或裁剪等操作，以保證地質(zhì)模型拓?fù)潢P(guān)系的正確性。較復(fù)雜地質(zhì)模型存在三角網(wǎng)的尺寸差異過大、密度不均勻、浮點誤差等現(xiàn)象，會導(dǎo)致曲面間的拓?fù)潢P(guān)系存在少量錯誤，使傳統(tǒng)的拓?fù)湟恢滦詨K體構(gòu)建方法失效。然而少量的拓?fù)溴e誤并不會影響模型整體的正確性，如果能構(gòu)建出結(jié)構(gòu)正確的塊體模型，則少量的拓?fù)溴e誤和塊體局部漏縫并不會對射線或模型離散化造成明顯的影響，可以繼續(xù)用于正演數(shù)值模擬[14]。

本文提出一種非拓?fù)湟恢碌娜S地質(zhì)塊體構(gòu)建技術(shù)。該方法模擬人眼視覺觀察，完全不依賴交線拓?fù)湫畔?，僅依靠曲面自身形態(tài)提取塊體三角網(wǎng)，在不對原曲面進行求交或裁剪等操作的情況下構(gòu)建三維地質(zhì)塊體。

1 非拓?fù)湟恢滦詨K體構(gòu)建算法原理

以交線拓?fù)浣Y(jié)構(gòu)為基礎(chǔ)的傳統(tǒng)塊體構(gòu)建算法根據(jù)三角網(wǎng)鄰接關(guān)系建立模型的拓?fù)潢P(guān)系圖，并在此基礎(chǔ)上進行廣度或深度遍歷算法進行塊體外邊界三角網(wǎng)的追蹤[15]。因此拓?fù)潢P(guān)系的正確性決定著此方法的成功率。然而在面對有少量拓?fù)鋯栴}的模型時，通常會從整體把握理解模型結(jié)構(gòu)，忽略局部拓?fù)溴e誤?；诖怂悸?，本文提出非拓?fù)湟恢滦詨K體構(gòu)建算法，模擬人類觀察和處理模型的方式，即首先從整體上確定塊體的大致范圍，再利用局部觀察的方式確定塊體的邊界三角網(wǎng)(圖1)。

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

圖1 非拓?fù)湟恢聣K體構(gòu)建算法流程

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

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

非拓?fù)湟恢滦詨K體構(gòu)建算法首先要確定塊體的基本形態(tài)。為此本文利用漫水填充法追蹤塊體所包含的網(wǎng)格單元集合。漫水法是圖形學(xué)中常用的一種區(qū)域填充法，實質(zhì)上是一種廣度優(yōu)先遍歷的種子填充法，適用于對內(nèi)定義區(qū)域的填充。圖形學(xué)中的內(nèi)定義區(qū)域，是指區(qū)域內(nèi)部的像素具有類似或相同的屬性，而區(qū)域外的所有像素具有不同或差別較大的屬性。該算法從內(nèi)區(qū)域中的某個像素出發(fā)(通常稱起始像素為種子)，利用廣度優(yōu)先搜索完成內(nèi)區(qū)域填充。本文對該方法進行改造，用于三維空間中地質(zhì)塊體網(wǎng)格單元集合的追蹤，具體步驟如下：

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

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

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

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

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

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

圖2 非拓?fù)湟恢聣K體構(gòu)建示意圖(a)輸入曲面模型； (b)對模型空間進行網(wǎng)格劃分； (c)漫水法追蹤塊體外邊界網(wǎng)格單元集合； (d)視覺觀察得到的三角網(wǎng)塊體模型

2 非拓?fù)湟恢聣K體構(gòu)建算法優(yōu)化

在上述非拓?fù)湟恢聣K體構(gòu)建算法中，網(wǎng)格劃分是重要步驟。基本的網(wǎng)格劃分方式是均勻型網(wǎng)格劃分，其劃分方式為無區(qū)別地將整個地質(zhì)模型均勻劃分為規(guī)模統(tǒng)一、精度一致的網(wǎng)格單元。這種劃分方式簡單，便于實現(xiàn)，且索引效率較高。其缺陷是既耗費資源又耗費時間。為此，本文提出一種高效方法，即非規(guī)則八叉樹網(wǎng)格劃分方式。該方法將地質(zhì)模型劃分為大小不一的網(wǎng)格單元，即對含有曲面三角網(wǎng)的區(qū)域劃分為規(guī)模更小、精度更高的網(wǎng)格單元，對不含有三角網(wǎng)的區(qū)域不再繼續(xù)劃分。由于使用了精度可變的網(wǎng)格體，所需求的內(nèi)存空間較靈活，可根據(jù)場景需求調(diào)整[19]。圖3對比了均勻網(wǎng)格體劃分與不規(guī)則八叉樹網(wǎng)格體劃分。

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

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

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

圖4 八叉樹網(wǎng)格劃分算法流程

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

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

表1 不規(guī)則八叉樹劃分與規(guī)則網(wǎng)格劃分的網(wǎng)格單元數(shù)量對比

3 效果展示與對比

為了驗證非拓?fù)湟恢滦运惴ǖ挠行?，選擇了一個大小為32km(x)×18km(y)×15km(z)的較復(fù)雜地質(zhì)模型(圖5)進行算法測試。模型含有44個層位/斷層曲面，存在斷裂、尖滅、透鏡體等多種地質(zhì)現(xiàn)象。為了驗證算法處理非拓?fù)湟恢滦阅Ｐ偷哪芰?，在圖中A、B兩處人為地制造了一些曲面漏縫(圖6左)和交越(圖6右)，可以看到模型存在比較明顯的拓?fù)鋯栴}。

圖7是利用傳統(tǒng)曲面縫合性方法建立的塊狀模型?？梢钥闯?，當(dāng)模型存在拓?fù)鋯栴}時，曲面縫合性方法得到的塊體模型結(jié)構(gòu)存在明顯問題，即局部存在拓?fù)溴e誤的塊體粘連在一起，得到的塊體與原始模型追蹤得到塊體的結(jié)構(gòu)和數(shù)目存在差異。

圖8是用網(wǎng)格塊體構(gòu)建方法建立的塊狀模型，圖9是網(wǎng)格塊體構(gòu)建算法細(xì)節(jié)?？梢钥闯?，網(wǎng)格塊體構(gòu)建算法的整體構(gòu)建效果較好，但對存在拓?fù)鋯栴}的模型，塊體模型沒有保留這些細(xì)節(jié)。

圖5 非拓?fù)湟恢履Ｐ?/p>

圖6 圖5模型中細(xì)節(jié)A(左)和B(右)

圖10是利用本文非拓?fù)湟恢滦运惴ń⒌膲K狀模型。塊體構(gòu)建采用Intel i7 4790K、主頻為4.0GHz的電腦，耗時23s。

圖7 傳統(tǒng)曲面縫合性算法對圖5模型的追蹤結(jié)果

圖8 網(wǎng)格塊體構(gòu)建算法對圖5模型的追蹤結(jié)果

圖9 網(wǎng)格塊體構(gòu)建算法對圖5模型的追蹤結(jié)果細(xì)節(jié)

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

對比圖10與圖8可以看出，當(dāng)模型存在一定拓?fù)鋯栴}時，本文非拓?fù)湟恢聣K體構(gòu)建算法得到的塊狀模型整體結(jié)構(gòu)正確，塊體結(jié)構(gòu)和數(shù)目與圖8一致。圖10右給出了非拓?fù)湟恢滦运惴ㄔ趫D中A、B兩處得到的塊體模型細(xì)節(jié)?？梢娫谠寄Ｐ驮诮痪€處的漏縫，在本文非拓?fù)湟恢滦运惴ń⒌膲K狀模型中也呈現(xiàn)出一定的缺漏現(xiàn)象。

實驗結(jié)果說明非拓?fù)湟恢滦运惴ň哂泻芎玫娜蒎e能力，可以較好地處理存在一定拓?fù)鋯栴}的模型。

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

表2 塊體構(gòu)建測試模型

若某次塊體構(gòu)建的結(jié)果中出現(xiàn)本應(yīng)被分離的兩個塊體黏合到一起，則本次塊體構(gòu)建結(jié)果是失敗的，此時不記錄時間和空間資源的消耗情況。

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

表3 五個模型塊體構(gòu)建算法效果對比

注：算法1指曲面縫合型塊體構(gòu)建算法，算法2指線框塊體構(gòu)建算法，算法3指基于視覺觀察效果的非拓?fù)湟恢聣K體構(gòu)建算法，表4同。

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

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

4 結(jié)論

三維塊體構(gòu)建技術(shù)是構(gòu)造塊狀地質(zhì)模型的核心技術(shù)，在此過程中，要保證復(fù)雜模型的拓?fù)湟恢滦苑浅＠щy。傳統(tǒng)的逐三角形追蹤和線框追蹤方法在模型存在極少量局部拓?fù)溴e誤時能夠建立比較正確的塊狀模型，但對交線拓?fù)浣Y(jié)構(gòu)依賴性強，容錯性和塊體構(gòu)建成功率較低。本文提出的非拓?fù)湟恢滦詨K體構(gòu)建算法，完全不依賴模型交線拓?fù)浣Y(jié)構(gòu)，根據(jù)曲面模型自身形態(tài)即可構(gòu)建，極大地提升了塊體構(gòu)建成功率。然而相比于傳統(tǒng)追蹤算法，該算法效率相對較低，為此應(yīng)用基于八叉樹結(jié)構(gòu)提高了算法效率。通過數(shù)個典型三維模型的塊體構(gòu)建實驗，可以得出以下結(jié)論：

(1)本文基于視覺觀察方式的三角形選取算法提出的非拓?fù)湟恢滦詨K體構(gòu)建算法，不依賴于交線拓?fù)浣Y(jié)構(gòu)，算法容錯性較高，能夠在很大程度上解決復(fù)雜模型拓?fù)洳灰恢碌膯栴}。

(2)非拓?fù)湟恢聣K體構(gòu)建算法以網(wǎng)格塊體追蹤為基礎(chǔ)，追蹤效率較低，且消耗資源較大，使用八叉樹結(jié)構(gòu)對算法進行優(yōu)化，可以極大地降低內(nèi)存消耗，一定程度上提高了追蹤效率。

(3)本文提出基于八叉樹優(yōu)化的非拓?fù)湟恢滦詨K體構(gòu)建算法，能夠有效解決拓?fù)洳灰恢履Ｐ偷膲K體構(gòu)建問題，具有較高的容錯性和塊體構(gòu)建成功率，并且算法資源消耗在目前主流硬件平臺上即可實現(xiàn)。由于該算法降低了對曲面模型的要求，因而也降低了三維復(fù)雜地質(zhì)模型構(gòu)建的整體難度，提高了塊體構(gòu)建的效率，比傳統(tǒng)方法更具實用性。