六面體體元網(wǎng)格三維地質(zhì)模型剖切算法①

張文東, 明志強(qiáng), 劉培剛

(中國石油大學(xué) 計(jì)算機(jī)與通信工程學(xué)院, 青島 266580)

六面體體元網(wǎng)格三維地質(zhì)模型剖切算法①

張文東, 明志強(qiáng), 劉培剛

(中國石油大學(xué) 計(jì)算機(jī)與通信工程學(xué)院, 青島 266580)

針對常用的六面體體元網(wǎng)格三維地質(zhì)模型, 提出了一種求剖切面的算法. 首先, 采用分層投影求交點(diǎn)的方式, 將地質(zhì)體模型與切割面投影到同一平面, 三維空間下的地質(zhì)體模型與切割面的剖切轉(zhuǎn)化為二維平面上的四邊形網(wǎng)格與切線段求交點(diǎn)的運(yùn)算. 為減少判交次數(shù), 先根據(jù)切線走勢判斷可能存在交點(diǎn)的區(qū)域, 再對可能區(qū)域進(jìn)行精確判交. 其次, 找到并求出不能通過投影方式得到的交點(diǎn). 然后, 將得到的所有交點(diǎn)按規(guī)則組成四邊形網(wǎng)格, 對每個(gè)四邊形三角化處理得到TIN形式的剖切面. 最后, 對該TIN面進(jìn)行顯示. 實(shí)驗(yàn)證明了對六面體體元網(wǎng)格三維地質(zhì)模型剖切的可行性.

三維地質(zhì)模型; 六面體體元; 剖面切割; 投影法; 不規(guī)則三角網(wǎng)

在三維地質(zhì)體的可視化研究領(lǐng)域, 對地質(zhì)模型進(jìn)行剖面切割是一個(gè)重要的組成部分. 剖切后生成的剖切面方便分析不透明地質(zhì)體基本內(nèi)部構(gòu)造及屬性. 常用的三維地質(zhì)建模方式有表面模型(主要是TIN模型)、體元模型和混合模型[1]. 不同的建模形式也決定了不同的剖切方法. 文獻(xiàn)[2]提出基于TIN數(shù)據(jù)三維地質(zhì)體的折剖面切割算法, 提高了算法效率, 但地質(zhì)體模型基于邊界表達(dá), 剖切后形成的剖面不能直觀有效的地表現(xiàn)內(nèi)部構(gòu)造. 文獻(xiàn)[3]提出了基于四面體格網(wǎng)模型的地質(zhì)體剖面生成算法, 該方法在求邊與切面的交點(diǎn)時(shí), 四面體的六個(gè)邊都要與切面進(jìn)行判斷, 并有重復(fù)計(jì)算的情況, 降低了執(zhí)行效率. 文獻(xiàn)[4]提出一種基于AutoCAD平臺的三維地質(zhì)實(shí)體模型自動生成地質(zhì)剖面的方法,該方法借助于現(xiàn)有的三維實(shí)體模型, 省去大量計(jì)算, 但剖切的三維實(shí)體模型是已存在的, 缺乏靈活性. 文獻(xiàn)[5,6]引入虛擬現(xiàn)實(shí)中碰撞檢測包圍盒的方法, 可以有效的減少求交次數(shù), 但包圍盒樹的構(gòu)建與存儲相對復(fù)雜. 文獻(xiàn)[7,8]提出基于八叉樹的快速構(gòu)建三維地質(zhì)剖面的算法, 該算法在一定程度上提高了切割效率, 降低時(shí)間復(fù)雜度. 本文提出六面體體元網(wǎng)格三維地質(zhì)模型剖切方法, 被剖切的地質(zhì)體模型是由多個(gè)六面體體元組成, 每個(gè)體元又由8個(gè)相鄰原始地質(zhì)數(shù)據(jù)點(diǎn)按一定規(guī)則組成,如圖1所示, 剖切時(shí)切割面會與各六面體體元的邊相交,采用分層投影求交點(diǎn)的方式求出各交點(diǎn), 所有交點(diǎn)所構(gòu)成的平面就是所求剖切面. 通過此方法得到的剖面可更準(zhǔn)確客觀地表達(dá)模型內(nèi)部結(jié)構(gòu).

圖1 三維坐標(biāo)系下的地質(zhì)模型

1 算法說明

實(shí)現(xiàn)三維地質(zhì)模型的剖切, 首先要確定三維地質(zhì)模型, 選擇合適的切割面并進(jìn)行剖切處理, 最后將剖切面三維顯示. 本文使用六面體體元網(wǎng)格的地質(zhì)體模型,該模型是一種結(jié)構(gòu)化的體元模型, 基本結(jié)構(gòu)是一系列六面體, 每個(gè)體元由8個(gè)數(shù)據(jù)點(diǎn)組成且有獨(dú)立的描述和存儲. 與四面體模型相比, 六面體在達(dá)到模型要求精度所用的節(jié)點(diǎn)數(shù)及網(wǎng)格數(shù)明顯少于四面體模型[9]. 圖2為三維地質(zhì)剖切整體流程圖. 本文將重點(diǎn)放在剖切面的生成上.

將地質(zhì)體模型置于三維笛卡爾坐標(biāo)系X-Y-Z中, 如圖1所示, 用X方向表示地質(zhì)體長, Y方向表示地質(zhì)體寬,Z方向表示地質(zhì)體深. 因?yàn)樯傻刭|(zhì)體模型的數(shù)據(jù)個(gè)數(shù)是已知的, 所以每個(gè)具體的地質(zhì)體模型X, Y, Z三個(gè)方向上數(shù)據(jù)點(diǎn)的個(gè)數(shù)是確定的, 根據(jù)數(shù)據(jù)點(diǎn)個(gè)數(shù)可算出每個(gè)方向上六面體體元數(shù)量. 假定在Z軸方向有k個(gè)有效數(shù)據(jù)點(diǎn), 說明在深度上有(k-1)層六面體體元, X軸方向若有i個(gè)數(shù)據(jù)點(diǎn), 地質(zhì)體長度上就有(i-1)層六面體體元, Y軸方向若有j個(gè)數(shù)據(jù)點(diǎn), 寬度方向就有(j-1)層六面體體元.

圖2 三維地質(zhì)模型剖切流程圖

2 算法步驟

實(shí)際的三維地質(zhì)體模型多由非規(guī)則六面體體元組成, 圖3為P1、P2、P3、P4、P5、P6、P7、P8八個(gè)數(shù)據(jù)點(diǎn)組成的單個(gè)體元剖切情況圖, 六個(gè)面的大小不同.當(dāng)剖切面為四邊形A’B’C’D’時(shí), 切面只與體元的頂面或底面相切, 使用投影求交點(diǎn). 當(dāng)剖切面為四邊形ABCD, 出現(xiàn)剖切面與體元側(cè)面相切的情況, 圖中點(diǎn)E為側(cè)面交點(diǎn)無法通過投影獲取, 需要單獨(dú)計(jì)算.

2.1 分層投影求交

進(jìn)行分層投影求交, 分別將切割面與被切割地質(zhì)體投影到二維X-0-Y面上(算法中用到的切割面要求不與平面X-0-Y平行), 分層投影次數(shù)由地質(zhì)模型在Z軸方向的六面體體元層數(shù)決定, 每次投影后三維地質(zhì)體變成只保留長和寬的四邊形網(wǎng)格結(jié)構(gòu), 切割面變成一條線段, 把該線段記為L. 求三維剖切面變成二維上求線段與四邊形各邊交點(diǎn)的運(yùn)算.

圖3 單個(gè)體元剖切圖

圖4為一次投影結(jié)束后四邊形網(wǎng)格與線段L相交情況. 假設(shè)地質(zhì)體在Z軸方向有(k-1)層六面體體元, 每次投影后的平面是(i*j)個(gè)頂點(diǎn)組成的四邊形網(wǎng)格, 就會有(i-1)*(j-1)個(gè)四邊形, 每個(gè)四邊形四條邊分別記為:P1P2、P2P3、P3P4、P4P1.

圖4 地質(zhì)體模型與切割面投影后的二維形式

投影求交具體步驟為:

(1) 地質(zhì)體模型與切割面完成一次投影, 分別得到四邊形網(wǎng)格與線段L.

(2) 從網(wǎng)格坐標(biāo)原點(diǎn)起始位置開始每次選擇一個(gè)四邊形, 讓其四條邊分別與線段L判交. 如果邊P1P2與線段L有交點(diǎn), 首先在二維空間下求出平面坐標(biāo)(X, Y),然后映射回三維空間下求出三維交點(diǎn)坐標(biāo)(X, Y, Z). 已知三維空間兩點(diǎn)坐標(biāo)P1、P2, 設(shè)點(diǎn)P1為(Xp1, Yp1, Zp1),點(diǎn)P2(Xp2, Yp2, Zp2), 由空間直線兩點(diǎn)式方程可求出過P1、P2的直線方程.

因?yàn)榻稽c(diǎn)坐標(biāo)過空間直線方程P1P2, 將二維空間求得的交點(diǎn)坐標(biāo)(X, Y)帶入公式(1)可求出三維交點(diǎn)坐標(biāo)(X, Y, Z). P2P3、P3P4、P4P1依次與線段L判交, 并記錄每次的交點(diǎn)坐標(biāo).

(3) 選擇下一個(gè)四邊形并判斷其合法性, 若沒有超出網(wǎng)格四邊形總個(gè)數(shù)(i-1)*(j-1), 執(zhí)行該四邊形與線段L的判交操作, 超出就退出本次投影求交.

(4) 對地質(zhì)模型體元結(jié)構(gòu)進(jìn)行下一層的投影, 判斷是否超出Z軸方向總層數(shù)(k-1), 若沒有超出繼續(xù)執(zhí)行步驟(2)(3), 否則就退出全部投影求交.

2.2 投影求交優(yōu)化

如果線段L依次與網(wǎng)格上每個(gè)邊進(jìn)行判交, 判交次數(shù)過多. 在實(shí)際的剖切過程中, 切割面不會與地質(zhì)模型的所有邊都有交點(diǎn), 找出可能存在交點(diǎn)的區(qū)域以縮小判交范圍. 因此在傳統(tǒng)判交基礎(chǔ)上進(jìn)行了改進(jìn), 減少不必要的求交判斷以提高效率.

在切割面與平面X-0-Y不平行的前提下, 投影后地質(zhì)體模型與切割面的關(guān)系會有圖5(a)(b)(c)三種情況,切割面投影后形成的線段L的斜率依次為負(fù)值、正值和零. 可在每次判交前確定線段L的斜率, 根據(jù)斜率確定判交的方向走勢.

圖5 投影后地質(zhì)模型與切割面3種關(guān)系圖

對圖5(a)的情況, 先確定第一個(gè)交點(diǎn)所在四邊形位置, 標(biāo)記該交點(diǎn)出現(xiàn)在X軸方向本行第m1個(gè)四邊形上;計(jì)算線段L斜率, 這里為負(fù)值, 可初步判斷下次出現(xiàn)交點(diǎn)的四邊形只能是下一行第m1個(gè)或m1之后的位置;Y軸方向行數(shù)加1, 新行的判交從該行第m1個(gè)四邊形開始, 忽略m1之前的判交操作. 標(biāo)記新行出現(xiàn)交點(diǎn)的四邊形為第m2個(gè); 依次類推執(zhí)行至結(jié)束. 圖5(b)的情況, 首先找到第一個(gè)交點(diǎn)出現(xiàn)的四邊形, 標(biāo)記該交點(diǎn)出現(xiàn)在X軸方向本行第n1個(gè)四邊形上; 計(jì)算線段L斜率, 這里為正值, 可初步判斷下次出現(xiàn)交點(diǎn)的四邊形只能是第n1個(gè)或n1之前的位置; 行數(shù)加1, 新行的判交操作從第n1個(gè)開始依次往前遞減執(zhí)行, 忽略第n1個(gè)之后的判交.標(biāo)記新行出現(xiàn)交點(diǎn)的四邊形為第n2個(gè); 依次類推執(zhí)行至結(jié)束. 圖5(c)的情況, 首先找到第一個(gè)交點(diǎn)同樣標(biāo)記位置; 計(jì)算線段L斜率, 這里為1; 此種情況所有交點(diǎn)坐標(biāo)的X值相同, 按順序求出對應(yīng)不同Y值與Z值即可.

除去地質(zhì)模型最后一層, 其他各層六面體體元底層與相鄰下一層的頂層由相同的頂點(diǎn)組成, 因此每次投影后選擇只對每個(gè)體元結(jié)構(gòu)的頂面判交. 再加上最后一層體元結(jié)構(gòu)的底面判交結(jié)果.

2.3 計(jì)算側(cè)面交點(diǎn)

在每層投影求交時(shí)記錄各交點(diǎn)出自該層哪一個(gè)四邊形, 四邊形用索引值區(qū)分. 圖3中頂面M與底面M’有相同索引值, 頂面M中有交點(diǎn)A、B、A’、B’, 底面M’有交點(diǎn)C’、D’, 交點(diǎn)C、D落在M’之外, 此時(shí)體元出現(xiàn)與側(cè)面相切.

P2、P6為已知數(shù)據(jù)點(diǎn), 設(shè)點(diǎn)P2為(Xp2, Yp2, Zp2),P6(Xp6, Yp6, Zp6), 由空間直線兩點(diǎn)式方程(2)確定過P2、P3的直線方程.

投影求交時(shí)得到A、B、C、D坐標(biāo)值, 取其任意三點(diǎn)唯一確定平面ABCD. 設(shè)點(diǎn)A(Xa, Ya, Za), B(Xb, Yb, Zb),C(Xc, Yc, Zc). 已知三點(diǎn)不在同一直線, 由空間平面三點(diǎn)式方程(3)求出過A、B、C三點(diǎn)的平面方程. 聯(lián)立方程(2)與方程(3)求出交點(diǎn)E.

2.4 生成剖切面

判交結(jié)束得到投影層及側(cè)面交點(diǎn)坐標(biāo), 將各交點(diǎn)坐標(biāo)按一定規(guī)則相連得到要求的剖面. 如圖6所示, 第t(1<=t

對四邊形網(wǎng)格中每個(gè)四邊形進(jìn)行三角化處理[10],每個(gè)四邊形分割為兩個(gè)三角形, 如圖7, 四邊形P1’P2’P3’P4’可轉(zhuǎn)化為三角形P1’P3’P2’和三角形P1’P4’P3’. 最終將四邊形網(wǎng)格形式的剖切面轉(zhuǎn)化成三角面片網(wǎng)格形式的剖切面, 方便剖切面三維顯示.

圖7 四邊形三角化轉(zhuǎn)換圖

3 算法實(shí)驗(yàn)

為驗(yàn)證算法的可行性及剖切效率, 使用本算法中六面體體元模型和文獻(xiàn)[3]中四面體格網(wǎng)模型分別實(shí)現(xiàn)某地區(qū)地質(zhì)體的三維顯示并剖切, 切面投影后切線斜率為負(fù)值. 實(shí)驗(yàn)涉及到的硬件配置為: 操作系統(tǒng)windows7旗艦版, CPU主頻2.2 GHz, 安裝內(nèi)存4 G(3.5 G可用). 兩種方式都能獲得近似一致的地質(zhì)模型及剖切面. 圖8(a)為六面體體元網(wǎng)格生成的三維地質(zhì)模型,8(b)為生成的地質(zhì)剖切面. 由實(shí)驗(yàn)效果圖可知利用該方法可以準(zhǔn)確的得到剖切面.

圖8 實(shí)驗(yàn)效果圖

由表1的實(shí)驗(yàn)結(jié)果可知, 六面體體元網(wǎng)格地質(zhì)模型在剖切時(shí)間上明顯比四面體格網(wǎng)模型地質(zhì)體剖切時(shí)短,且數(shù)據(jù)點(diǎn)個(gè)數(shù)越多, 效果越明顯.

表1 兩種地質(zhì)模型組織方式剖切時(shí)間比較

4 結(jié)語

本文針對六面體體元網(wǎng)格的實(shí)體地質(zhì)模型求剖切面. 結(jié)合該模型的特點(diǎn)在求交時(shí)主要使用了分層投影求交點(diǎn)的方法, 并對特殊切點(diǎn)進(jìn)行了專門處理, 通過實(shí)驗(yàn)證明了此方法的可行性. 另外本算法在剖切時(shí)通過對切線走勢進(jìn)行預(yù)判斷以及投影求交的方式與文獻(xiàn)[3]中四面體格網(wǎng)模型的地質(zhì)體剖切相比, 減少了判交次數(shù), 提高了整體剖切效率. 本次求出的剖切面只以同一顏色進(jìn)行了三維顯示. 由于每個(gè)體元的屬性可以被獨(dú)立描述和存儲, 每個(gè)原始數(shù)據(jù)都有自己的屬性值, 而判交得到的新數(shù)據(jù)點(diǎn)屬性值可根據(jù)鄰近原始數(shù)據(jù)點(diǎn)獲得, 不同屬性值又可用不同顏色區(qū)分. 后期為更有效地展現(xiàn)體元網(wǎng)格剖切的優(yōu)勢, 會將每個(gè)數(shù)據(jù)點(diǎn)屬性值應(yīng)用到三維顯示上, 得到一個(gè)用不同顏色值區(qū)分不同屬性區(qū)的三維地質(zhì)剖切面.

1劉蘇. 三維地質(zhì)模型可視化方法及應(yīng)用. 中國高新技術(shù)企業(yè), 2015, (30): 133–134.

2明鏡, 潘懋, 屈紅剛, 等. 基于TIN數(shù)據(jù)三維地質(zhì)體的折剖面切割算法. 地理與地理信息科學(xué), 2008, 24(3): 37–40.

3沈敬偉, 周廷剛, 劉德兒, 等. 一種基于四面體格網(wǎng)模型的地質(zhì)體剖面生成算法. 西南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2014,36(8): 123–129.

4劉光偉, 白潤才, 呂進(jìn)國, 等. 基于三維地質(zhì)實(shí)體模型生成地質(zhì)剖面圖的應(yīng)用. 遼寧工程技術(shù)大學(xué)學(xué)報(bào)(自然科學(xué)版),2010, 29(4): 557–559.

5李運(yùn)鋒, 劉修國. 基于方向包圍盒投影轉(zhuǎn)換的輪廓線拼接算法. 計(jì)算機(jī)應(yīng)用, 2011, 31(12): 3353–3356.

6李偉波, 劉嘉, 陳耀華. 三維地層Tin模型剖切的改進(jìn)算法.計(jì)算機(jī)應(yīng)用與軟件, 2013, 30(8): 158–161,169.

7趙龍, 閔世平, 代強(qiáng)玲. 基于八叉樹的三維地質(zhì)剖面生成算法. 計(jì)算機(jī)工程, 2014, 40(2): 250–255.

8Delamé T, Roudet C, Faudot D. From a medial surface to a mesh. Computer Graphics Forum, 2012, 31(5): 1637–1646.[doi: 10.1111/cgf.2012.31.issue-5]

9任輝龍, 鞏生龍, 蔡永昌. 基于投影法三維模型全六面體網(wǎng)格的劃分方法. 計(jì)算機(jī)輔助工程, 2013, 22(5): 122–128.

10鄒新龍, 石丹, 劉茂, 等. 簡單多邊形的三角化算法. 環(huán)境技術(shù), 2014, (增刊): 137–140.

Hexahedral Voxel Grid 3D Geological Model Partitioning Algorithm

ZHANG Wen-Dong, MING Zhi-Qiang, LIU Pei-Gang
(College of Computer and Communication Engineering, China University of Petroleum, Qingdao 266580, China)

In view of the common hexahedral voxel grid 3D geological model, we introduce a 3D geological model partitioning algorithm. In the process of cutting, the cut surface will intersect with the edge of the hexahedral element.Firstly, the geological model and the cut plane are projected into the same 2D plane at the same time, the process of portioning between geological model and cut plane in the 3D is converted to an operation that looks for the intersection point between the line and the quadrilateral grid. In order to reduce the number of judgement intersections, we would find out the possible intersection area by the Line slope direction of cutting plane, then judge whether there is an intersection carefully. Next, we find other points of intersection that could not be got by projection, and connect the node coordinates according to certain rules form the quadrilateral grid and triangulation. Finally, we display the TIN. Experimental results prove that this method is feasible.

three-dimensional geological model; hexahedral voxel; profile cutting; projection method; TIN

張文東,明志強(qiáng),劉培剛.六面體體元網(wǎng)格三維地質(zhì)模型剖切算法.計(jì)算機(jī)系統(tǒng)應(yīng)用,2017,26(7):195–199. http://www.c-s-a.org.cn/1003-3254/5858.html

2016-10-27; 收到修改稿時(shí)間: 2016-12-12