非拓?fù)湟恢滦匀S地質(zhì)塊體建模算法研究

羅開云， 隆振海， 何永清， 蔣先藝， 吳 恒， 宋衛(wèi)鋒

(1.中國石油集團 東方地球物理公司采集技術(shù)中心，保定 072751； 2.成都理工大學(xué) 信息科學(xué)與技術(shù)學(xué)院，成都 610059)

0 引言

三維塊狀地質(zhì)建模[1]以具有近似地質(zhì)屬性的層塊作為建模的基本單元，首先利用曲面構(gòu)建方法建立地質(zhì)層位和斷層曲面，再利用塊體追蹤技術(shù)提取由層位和斷層曲面片所形成的封閉地質(zhì)塊體，建模所得到的塊狀地質(zhì)模型通常可以用于射線、射線束類正演模擬[2]，也可以通過離散化算法轉(zhuǎn)換為規(guī)則或者非規(guī)則網(wǎng)格模型，支持有限差分或有限元方法[3]。三維塊狀建模方法符合人們探索和認(rèn)知地下結(jié)構(gòu)的規(guī)律，利用塊體單元可以描述斷裂、尖滅、透鏡體、侵入體等多種地質(zhì)結(jié)構(gòu)，因而是目前主流的建模方法，被廣泛應(yīng)用于GoCAD、Petrel等著名地質(zhì)建模軟件中。

三維塊狀建模算法的歷史可以追溯到20世紀(jì)80年代，Warbuton[4]提出一種基于無限大非連續(xù)面來識別塊體的算法；Lin等[5]提出將拓?fù)涓拍钜氲綁K體建模中；Lu[6]將方向邊和方向環(huán)的概念引入了塊體建模過程;蔣先藝等[7]提出將具有相同屬性的地質(zhì)層面聯(lián)系在一起，并且給出了一種基于三角形的塊體追蹤算法，從而得出地質(zhì)塊體；侯衛(wèi)生[8]在Lu的基礎(chǔ)上進一步改進，使得算法可以在曲面上進行搜索從而得出塊體；楊洋等[9]簡化并改進了傳統(tǒng)的線框架模型，提出了線框架拓?fù)淠Ｐ蛠砑铀倥c簡化拓?fù)渌阉鞯倪^程。

經(jīng)過數(shù)十年的發(fā)展三維塊狀建模算法日漸成熟，在塊體追蹤方面從簡單的逐三角形追蹤發(fā)展為基于線框拓?fù)涞淖粉?，不僅在效率上有了很大改進，而且也具有一定的容錯性。但本質(zhì)上，現(xiàn)有的方法都是基于模型交線拓?fù)浣Y(jié)構(gòu)進行塊體追蹤，仍需要利用交線的拓?fù)湫畔?。然而前期的曲面建模，特別是復(fù)雜模型的曲面建模，由于存在結(jié)構(gòu)復(fù)雜、三角網(wǎng)尺寸差異過大、密度不均勻、浮點誤差等原因，很容易造成曲面交線拓?fù)洳灰恢?。而?dāng)?shù)刭|(zhì)模型存在拓?fù)洳灰恢碌默F(xiàn)象時，基于交線拓?fù)浣Y(jié)構(gòu)的塊體追蹤方法就會出現(xiàn)不可預(yù)期的錯誤，造成塊體建模失敗。為了提高三維塊狀地質(zhì)建模的成功率和效率，本文提出了一種不依賴交線拓?fù)湫畔⒌膲K體追蹤算法——非拓?fù)湟恢碌膲K體追蹤算法，該方法模擬人眼視覺觀察，完全不依賴交線拓?fù)湫畔?，僅依靠曲面自身形態(tài)進行塊體三角網(wǎng)提取，因而具有很強的魯棒性，對存在拓?fù)溴e誤，甚至存在漏縫和交越現(xiàn)象的模型都能正常完成塊體建模。

1 非拓?fù)湟恢滦詨K體建模算法原理

傳統(tǒng)的塊體追蹤算法以交線拓?fù)浣Y(jié)構(gòu)為基礎(chǔ)，利用交線處三角形或者三角網(wǎng)的鄰接關(guān)系，建立模型的拓?fù)潢P(guān)系圖，并在此基礎(chǔ)上通過圖的廣度或者深度遍歷算法進行塊體外邊界三角網(wǎng)的追蹤。對于拓?fù)潢P(guān)系正確的模型，利用該方法進行塊體追蹤效率很高，但是當(dāng)模型的交線拓?fù)浣Y(jié)構(gòu)存在一定問題，地層或者斷層三角網(wǎng)在交線處的拓?fù)鋭澐植煌耆恢?，或者出現(xiàn)局部漏縫或者交越現(xiàn)象時，追蹤得到的塊體模型就可能會出現(xiàn)塊體缺失、相互粘連等錯誤。

然而面對有少量拓?fù)鋯栴}的模型，人們會從整體上理解模型結(jié)構(gòu)，忽略局部交線處的細(xì)微拓?fù)鋯栴}，即使在模型存在一定問題時也能正確分辨出模型中的各個地質(zhì)塊體。筆者提出的非拓?fù)湟恢滦詨K體追蹤算法，模擬了人類觀察和處理模型的方式，首先從整體上確定塊體的大致范圍，再利用局部觀察的方式確定塊體的邊界三角網(wǎng)。

圖1 非拓?fù)湟恢聣K體追蹤算法流程圖Fig.1 Flow chart of non-topological consistency blocking algorithm

非拓?fù)湟恢滦詨K體追蹤算法的基本流程如圖1所示，算法首先讀入經(jīng)曲面建模后的三維曲面模型(允許模型存在一定的拓?fù)溴e誤，并且無需讀取任何交線拓?fù)湫畔?；接下來對地質(zhì)模型空間進行均勻網(wǎng)格劃分，將曲面三角網(wǎng)中的三角形信息記錄在其對應(yīng)位置的網(wǎng)格單元中；然后使用漫水網(wǎng)格單元追蹤方法得出地質(zhì)塊體在空間中所占據(jù)的網(wǎng)格單元集合，得到地質(zhì)塊體的基本輪廓；最后再使用類似于人眼視覺觀察的方式對塊體邊界區(qū)域的網(wǎng)格單元進行處理，從中提取出地質(zhì)塊體所包含的外表面三角形，從而構(gòu)造地質(zhì)塊體外表面三角網(wǎng)。

為了直觀地說明非拓?fù)湟恢滦詨K體追蹤算法的關(guān)鍵步驟和執(zhí)行過程，圖2給出了算法執(zhí)行的具體過程，如圖2(a)所示算法首先讀取模型的曲面三角網(wǎng)，確定模型的空間范圍，將模型空間劃分成均勻網(wǎng)格單元，然后用漫水法追蹤塊體的網(wǎng)格單元集合得到如圖2(b)所示的由網(wǎng)格單元所表示的基本塊體模型，最后再通過模擬視覺觀察方式提取塊體的外邊界三角網(wǎng)，得到最終的塊體模型如圖2(c)所示。

2 漫水法追蹤塊體網(wǎng)格單元

非拓?fù)湟恢滦詨K體追蹤算法首先要確定塊體的基本形態(tài)，并大致確定塊體可能包含的三角形，為之后精確提取地質(zhì)塊體外表面三角形做好準(zhǔn)備。為此可以先將模型空間劃分成等距均勻的結(jié)構(gòu)化網(wǎng)格，并將曲面模型中的三角形根據(jù)坐標(biāo)位置放置到與之有交的網(wǎng)格單元中(網(wǎng)格單元是一個長方體)，然后就可以從一個不含有任何三角形的網(wǎng)格單元出發(fā)，利用漫水法追蹤塊體所包含的網(wǎng)格單元集合。

漫水法是圖形學(xué)中常用的一種區(qū)域填充法[10]，其實質(zhì)是一種廣度優(yōu)先遍歷的種子填充法，適用于對內(nèi)定義區(qū)域的填充。而圖形學(xué)中的內(nèi)定義區(qū)域，是指區(qū)域內(nèi)部的像素具有類似或相同的屬性，而區(qū)域外的所有像素具有不相同或差別較大的屬性，該算法從內(nèi)區(qū)域中的某個像素出發(fā)(通常稱起始像素為種子)，利用廣度優(yōu)先搜索完成內(nèi)區(qū)域填充?？梢詫υ摲椒ㄟM行改造，用于三維空間中地質(zhì)塊體網(wǎng)格單元集合的追蹤，三維空間中地質(zhì)塊體的漫水網(wǎng)格單元追蹤算法的具體步驟如下：

圖2 非拓?fù)湟恢聣K體追蹤算法效果示意圖Fig.2 Result of non-topological consistency blocking algorithm(a)A處細(xì)節(jié);(b)B處細(xì)節(jié);(c)C處細(xì)節(jié)

1)嘗試從三維地質(zhì)空間中取得一個未被訪問過且不含有三角形的網(wǎng)格體作為起始網(wǎng)格體，若不能找到一個所述不含有三角形且未被訪問過的網(wǎng)格體則塊體追蹤過程結(jié)束，否則將所述起始網(wǎng)格體加入到遍歷隊列中并進行步驟2)。

2)若遍歷隊列為空則重新進行步驟1)，否則從遍歷隊列中取出一個網(wǎng)格體進行步驟3)。

3)對所述取出的網(wǎng)格體，若其已經(jīng)被訪問過，則重新進行步驟2)，否則進行步驟4)。

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

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

利用上述漫水追蹤算法，可以提取模型中各個地質(zhì)塊體所包含的網(wǎng)格單元集合，這些網(wǎng)格單元確定了地質(zhì)塊體的基本形態(tài)，而且網(wǎng)格單元集合中包含位于塊體邊界上的網(wǎng)格單元，這些邊界網(wǎng)格單元包含曲面三角形，而其中部分三角形是屬于某個特定地質(zhì)塊體的，這為最終提取塊體外邊界三角形奠定了基礎(chǔ)。

3 基于視覺觀察的塊體三角網(wǎng)建模

非拓?fù)湟恢碌膲K體追蹤算法的最終目標(biāo)是提取地質(zhì)塊體的外邊界三角形，構(gòu)造地質(zhì)塊體表面三角網(wǎng)，通過漫水追蹤算法獲取了地質(zhì)塊體的網(wǎng)格單元集合，為了進一步提取塊體的外邊界三角形，需要遍歷塊體的網(wǎng)格單元，從含有三角形的邊界網(wǎng)格單元中找出屬于與地質(zhì)塊體的三角形，并且將它們組合起來構(gòu)造塊體外表面三角網(wǎng)。

圖3 基于視覺觀察的塊體三角網(wǎng)建模流程圖Fig.3 Flow chart of blocking algorithm based on visual observation

從邊界網(wǎng)格單元中提取三角形這一過程需要基于模仿人眼視覺觀察的方式來進行，方能達到類似于人眼“所見即所得”的觀察效果。對于含有三角形的邊界區(qū)域網(wǎng)格單元，通常它所含的三角形一部分屬于當(dāng)前地質(zhì)塊體，另一部分其他塊體，考慮到一個邊界網(wǎng)格單元的6個面中至少有一個是與內(nèi)部網(wǎng)格單元相鄰的(內(nèi)部網(wǎng)格單元指不含有任何三角形的網(wǎng)格單元)，可以將與內(nèi)部網(wǎng)格單元相鄰的面作為觀察面，通過該面觀察邊界網(wǎng)格單元中所包含的三角形，屬于當(dāng)前塊體的三角形必然能夠被觀察到，而不屬于當(dāng)前塊體的三角形則會被遮擋而無法觀察到，通過“視覺觀察”就可以準(zhǔn)確提取邊界網(wǎng)格單元中屬于當(dāng)前地質(zhì)塊體的三角形了。

當(dāng)然計算機并沒有眼睛，為了實現(xiàn)上述“視覺觀察”算法，可以利用OpenGL進行三維圖形渲染[11]，并對渲染得到的圖像進行判定，從而實現(xiàn)所謂的視覺觀察效果。首先確定觀察方向，并建立一個與其同位置，同大小的三維平行視景體[11-13]；然后將該網(wǎng)格單元所包含的三角形賦予不同顏色，再利用OpenGL三維圖形庫將場景渲染到幀緩存對象中[14]；最后從幀緩存中將渲染結(jié)果圖像取出，遍歷圖像中每個像素點的顏色，選擇在渲染結(jié)果中占據(jù)了一定像素比例的三角形作為當(dāng)前地質(zhì)塊體的外邊界三角形之一。

對上述算法思路進行總結(jié)，可以得到如圖3所示的基于視覺觀察的塊體三角網(wǎng)建模流程圖。在地質(zhì)塊體外邊界三角形的提取過程中，每一個邊界區(qū)域中的網(wǎng)格體均會被施以如下處理：

圖4 三維地質(zhì)模型示意圖Fig.4 3D geologic model

1)遍歷一次所述網(wǎng)格體中已保存的所有三角形信息，遍歷過程中檢查所述網(wǎng)格體中的三角形是否來自同一個地層面三角網(wǎng)。

2)若所述網(wǎng)格體中的三角形均來自于同一個地層面三角網(wǎng)，則直接將這些三角形判定為當(dāng)前塊體外邊界，否則進行步驟3)。

3)創(chuàng)建一個顏色表，所述顏色表中不存在兩個相同的顏色，并使顏色表中的顏色與三角形一一對應(yīng)。

4)使用一個與所述網(wǎng)格體同大小且同位置的視景體進行平行投影繪制得到繪制結(jié)果，其中每個三角形均使用其對應(yīng)所述顏色表中的顏色。

5)利用所述視景體對所述網(wǎng)格體中含有的三角形進行平行投影繪制得出繪制結(jié)果，其中每個三角形均使用顏色表中已記錄的一個顏色，三角形與顏色表中的顏色一一對應(yīng)。

6)統(tǒng)計出繪制結(jié)果中每種顏色繪制結(jié)果中的面積占比，所述顏色對應(yīng)面積占比大于2%的顏色認(rèn)定為有效顏色，有效顏色對應(yīng)的三角形認(rèn)定為有效三角形。

7)將有效三角形判定為當(dāng)前塊體的表面三角形。

4 算法應(yīng)用實例

為了驗證非拓?fù)湟恢滦运惴ǖ挠行?，我們選擇了一個13.5 km×13.5 km×8 km較為復(fù)雜的地質(zhì)模型進行算法測試，如圖4所示，該模型含有44個層位/斷層曲面，存在斷裂、尖滅、透鏡體等多種地質(zhì)現(xiàn)象。為了驗證算法處理非拓?fù)湟恢滦阅Ｐ偷哪芰?，在圖4中(a)、(b)、(c)三處人為地制造了一些曲面漏洞和交越問題，圖5給出了上述三處交線錯誤的放大圖，從圖5中可以看到模型存在比較明顯的拓?fù)鋯栴}。

圖5 三維地質(zhì)模型細(xì)節(jié)示意圖Fig.5 Details of 3D geologic model(a)A處細(xì)節(jié);(b)B處細(xì)節(jié);(c)C處細(xì)節(jié)

圖6 傳統(tǒng)追蹤方法結(jié)果圖Fig.6 Result of traditional blocking method(a)非拓?fù)湟恢滦阅Ｐ徒Y(jié)果;(b)拓?fù)湟恢滦阅Ｐ徒Y(jié)果

圖7 本文追蹤方法結(jié)果圖Fig.7 Result of our blocking algorithm(a)整體結(jié)果; (b)追蹤結(jié)果C處細(xì)節(jié)圖

圖6(a)是利用傳統(tǒng)的拓?fù)湟恢滦宰粉櫡椒ń⒌膲K狀模型，圖6(b)則是用傳統(tǒng)方法對原始沒有拓?fù)鋯栴}的模型進行追蹤建立的塊狀模型。從圖6中可以看出，當(dāng)模型存在拓?fù)鋯栴}時，傳統(tǒng)方法得到的塊體模型結(jié)構(gòu)存在問題，局部存在拓?fù)溴e誤的塊體粘連在一起，得到的塊體與原始模型追蹤得到塊體結(jié)構(gòu)和總數(shù)存在差異。

圖7(a)是利用本文提出的非拓?fù)湟恢滦运惴ń⒌膲K狀模型，在Intel i7 4790K主頻4.0 GHz的測試電腦上，建模時間為23 s。對比圖7(a)和圖6(b)可以看出，當(dāng)模型存在一定拓?fù)鋯栴}時非拓?fù)湟恢滦运惴ǖ玫降膲K狀模型整體結(jié)構(gòu)正確，塊體結(jié)構(gòu)和數(shù)目與圖6(b)一致。圖7(b)給出了非拓?fù)湟恢滦运惴ㄔ趫D7(c)處得到的塊體模型，因為原始模型在交線處存在漏縫，所以得到的塊體模型也存在一定的缺漏，但是這并沒有影響模型構(gòu)建的整體結(jié)果。說明非拓?fù)湟恢滦运惴ň哂泻芎玫娜蒎e能力，可以處理存在一定拓?fù)鋯栴}的模型。

5 結(jié)論

三維塊狀地質(zhì)建模目前主流的地質(zhì)建模方法，塊體追蹤技術(shù)是該方法的核心技術(shù)，現(xiàn)有的塊體追蹤技術(shù)均依賴于模型的交線拓?fù)浣Y(jié)構(gòu)，對模型拓?fù)涞囊恢滦院驼_性要求較高。筆者提出的非拓?fù)湟恢滦詨K體建模算法，首先利用漫水法確定地質(zhì)塊體的基本形態(tài)，然后利用基于OpenGL的視覺觀察技術(shù)實現(xiàn)塊體外邊界三角網(wǎng)的提取與構(gòu)建。該算法擺脫了對交線拓?fù)湫畔⒌囊蕾嚕瑑H根據(jù)曲面模型自身形態(tài)即可完成塊體構(gòu)建，因而對曲面模型的正確性要求大幅降低，當(dāng)曲面模型中存在交越、漏縫等現(xiàn)象時依然可以正確完成塊體模型的構(gòu)建，與傳統(tǒng)塊體追蹤技術(shù)相比，該方法降低了輸入模型的要求，具有很強的魯棒性。雖然該算法運行的效率低于傳統(tǒng)的拓?fù)湟恢滦宰粉櫡椒?，但是降低了輸入曲面模型的質(zhì)量要求，因而降低了建模的整體難度，能夠更好地輔助用戶完成復(fù)雜三維地質(zhì)模型的構(gòu)建，提高建模效率，具有比傳統(tǒng)方法更好的實際推廣使用價值。