六面體體元網格三維地質模型剖切算法①

張文東, 明志強, 劉培剛

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

六面體體元網格三維地質模型剖切算法①

張文東, 明志強, 劉培剛

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

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

三維地質模型; 六面體體元; 剖面切割; 投影法; 不規則三角網

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

圖1 三維坐標系下的地質模型

1 算法說明

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

將地質體模型置于三維笛卡爾坐標系X-Y-Z中, 如圖1所示, 用X方向表示地質體長, Y方向表示地質體寬,Z方向表示地質體深. 因為生成地質體模型的數據個數是已知的, 所以每個具體的地質體模型X, Y, Z三個方向上數據點的個數是確定的, 根據數據點個數可算出每個方向上六面體體元數量. 假定在Z軸方向有k個有效數據點, 說明在深度上有(k-1)層六面體體元, X軸方向若有i個數據點, 地質體長度上就有(i-1)層六面體體元, Y軸方向若有j個數據點, 寬度方向就有(j-1)層六面體體元.

圖2 三維地質模型剖切流程圖

2 算法步驟

實際的三維地質體模型多由非規則六面體體元組成, 圖3為P1、P2、P3、P4、P5、P6、P7、P8八個數據點組成的單個體元剖切情況圖, 六個面的大小不同.當剖切面為四邊形A’B’C’D’時, 切面只與體元的頂面或底面相切, 使用投影求交點. 當剖切面為四邊形ABCD, 出現剖切面與體元側面相切的情況, 圖中點E為側面交點無法通過投影獲取, 需要單獨計算.

2.1 分層投影求交

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

圖3 單個體元剖切圖

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

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

投影求交具體步驟為:

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

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

因為交點坐標過空間直線方程P1P2, 將二維空間求得的交點坐標(X, Y)帶入公式(1)可求出三維交點坐標(X, Y, Z). P2P3、P3P4、P4P1依次與線段L判交, 并記錄每次的交點坐標.

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

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

2.2 投影求交優化

如果線段L依次與網格上每個邊進行判交, 判交次數過多. 在實際的剖切過程中, 切割面不會與地質模型的所有邊都有交點, 找出可能存在交點的區域以縮小判交范圍. 因此在傳統判交基礎上進行了改進, 減少不必要的求交判斷以提高效率.

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

圖5 投影后地質模型與切割面3種關系圖

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

除去地質模型最后一層, 其他各層六面體體元底層與相鄰下一層的頂層由相同的頂點組成, 因此每次投影后選擇只對每個體元結構的頂面判交. 再加上最后一層體元結構的底面判交結果.

2.3 計算側面交點

在每層投影求交時記錄各交點出自該層哪一個四邊形, 四邊形用索引值區分. 圖3中頂面M與底面M’有相同索引值, 頂面M中有交點A、B、A’、B’, 底面M’有交點C’、D’, 交點C、D落在M’之外, 此時體元出現與側面相切.

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

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

2.4 生成剖切面

判交結束得到投影層及側面交點坐標, 將各交點坐標按一定規則相連得到要求的剖面. 如圖6所示, 第t(1<=t

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

圖7 四邊形三角化轉換圖

3 算法實驗

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

圖8 實驗效果圖

由表1的實驗結果可知, 六面體體元網格地質模型在剖切時間上明顯比四面體格網模型地質體剖切時短,且數據點個數越多, 效果越明顯.

表1 兩種地質模型組織方式剖切時間比較

4 結語

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

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

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

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

4劉光偉, 白潤才, 呂進國, 等. 基于三維地質實體模型生成地質剖面圖的應用. 遼寧工程技術大學學報(自然科學版),2010, 29(4): 557–559.

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

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

7趙龍, 閔世平, 代強玲. 基于八叉樹的三維地質剖面生成算法. 計算機工程, 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任輝龍, 鞏生龍, 蔡永昌. 基于投影法三維模型全六面體網格的劃分方法. 計算機輔助工程, 2013, 22(5): 122–128.

10鄒新龍, 石丹, 劉茂, 等. 簡單多邊形的三角化算法. 環境技術, 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

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

2016-10-27; 收到修改稿時間: 2016-12-12