三維地質(zhì)界面網(wǎng)格模型構(gòu)建算法研究

侯 曉 琳

(北京大學(xué)地球與空間科學(xué)學(xué)院 北京 100871)

0 引 言

在地質(zhì)學(xué)相關(guān)工程中通過多種技術(shù)手段獲取的多類型地質(zhì)數(shù)據(jù)不能夠直觀展示真實地質(zhì)情況。三維地質(zhì)建模技術(shù)利用地質(zhì)數(shù)據(jù)建立三維地質(zhì)模型，采用計算機可視化方法更加直觀地模擬出地下構(gòu)造情況，成為地質(zhì)學(xué)研究重點[1-2]。目前針對于三維地質(zhì)建模技術(shù)的研究成果已經(jīng)逐步應(yīng)用于生產(chǎn)實踐項目中，多種三維地質(zhì)建模軟件嘗試為石油勘探開發(fā)、礦產(chǎn)資源預(yù)測等項目提供三維定量化模型，為實際開發(fā)過程提供理論依據(jù)，降低開發(fā)成本[3]。

根據(jù)不同的地質(zhì)數(shù)據(jù)類型[4]，三維地質(zhì)模型可以分成兩類[5]：地質(zhì)構(gòu)造模型，通過解釋原始地質(zhì)構(gòu)造數(shù)據(jù)，描述地層、斷層等地質(zhì)界面空間形態(tài)展布和拓?fù)潢P(guān)系[6]；地質(zhì)屬性模型，采用分析預(yù)測方法，定量化描述地質(zhì)變量空間變化規(guī)律[7]。地質(zhì)構(gòu)造模型是三維地質(zhì)模型基礎(chǔ)[8]，其研究內(nèi)容包括地質(zhì)體模型建模方法、地質(zhì)層面建模方法以及地質(zhì)要素拓?fù)潢P(guān)系描述方法。從數(shù)據(jù)結(jié)構(gòu)角度，將地質(zhì)構(gòu)造模型分成面模型、體模型以及混合模型。通過建立三維地質(zhì)界面模型描述復(fù)雜地質(zhì)構(gòu)造，以地質(zhì)界面模型為約束建立三維地質(zhì)體模型。三維地質(zhì)層界面模型是描述地質(zhì)要素之間拓?fù)潢P(guān)系的基礎(chǔ)[9]。

本文在現(xiàn)有研究的基礎(chǔ)上，提出一種改進的三維地質(zhì)界面三角化網(wǎng)格模型構(gòu)建算法，采用實際的斷層面數(shù)據(jù)建立并可視化展示斷層面三角網(wǎng)格模型。與現(xiàn)有的基于映射關(guān)系建立三維三角化曲面模型算法相比，改進后的算法減少了映射關(guān)系計算過程，避免了可能出現(xiàn)的三維空間中多點映射到二維平面上一點而導(dǎo)致無法生成三維曲面模型的問題，改進的局部映射方法可在三維空間中直接生成曲面網(wǎng)格模型，該算法同樣也適用于解決其他領(lǐng)域三維三角化曲面構(gòu)建問題。

1 相關(guān)研究

三維地質(zhì)構(gòu)造模型的幾何模型包括點模型、線模型、面模型和體模型。面模型用于表示三維地質(zhì)體界面，反映界面間構(gòu)造關(guān)系，例如地層面、斷層面等。面模型主要類型包括不規(guī)則三角面模型、規(guī)則網(wǎng)格模型、邊界表示模型和線框模型等[10]。

在三維地質(zhì)體構(gòu)造建模中，通常采用規(guī)則網(wǎng)格和非規(guī)則網(wǎng)格描述地質(zhì)界面和地質(zhì)體。規(guī)則網(wǎng)格優(yōu)點在于網(wǎng)格剖分算法較為簡單，便于進行數(shù)值模擬計算。非規(guī)則網(wǎng)格則能夠更好地模擬地質(zhì)界面狀態(tài)和地質(zhì)體構(gòu)造形態(tài)，反映出地質(zhì)構(gòu)造的非均質(zhì)性[11]。與規(guī)則的四邊形網(wǎng)格相比，非規(guī)則三角化網(wǎng)格在描述曲面和不規(guī)則幾何體時具有一定優(yōu)勢。在三維地質(zhì)體建模中，表示地質(zhì)界面的原始數(shù)據(jù)為點集數(shù)據(jù)，例如，斷點數(shù)據(jù)描述地質(zhì)斷層界面，分層點數(shù)據(jù)表示地質(zhì)層面。不失一般性，三維地質(zhì)界面三角化網(wǎng)格模型構(gòu)造方法同樣可以應(yīng)用于基于點數(shù)據(jù)集構(gòu)建三維曲面模型的相關(guān)問題。

現(xiàn)有研究中主要采用映射法建立三維曲面網(wǎng)格模型，將三維數(shù)據(jù)點投影到二維平面上，在平面上進行二維三角剖分，根據(jù)點映射關(guān)系，將二維平面剖分結(jié)果映射到三維空間，生成三維曲面三角化網(wǎng)格模型[12]。在映射法中，映射函數(shù)和映射平面的選取方法具有一定不確定性，映射函數(shù)的影響生成三角化網(wǎng)格模型的形態(tài)[13]。很多非規(guī)則網(wǎng)格生成算法均采用映射法實現(xiàn)思路。經(jīng)典的三角化剖分算法有Lawson算法和Watson算法等[14-15]，這些算法實現(xiàn)思路清晰，并且能夠保證三角形網(wǎng)格質(zhì)量。但是，映射法在某些特定情況下也無法處理。三維點集中任意點p(x,y,z)，必須保證其在物理域和映射域中是一一映射關(guān)系。假設(shè)三維空間上的兩點p1(x,y,z1)和p2(x,y,z2)，映射到某一二維平面上的點為同一點p′(x,y)，在這種情況下映射法就很難處理，確定映射函數(shù)和選擇映射平面是基于映射法建立三維曲面模型的關(guān)鍵問題。當(dāng)三維空間中的不同點映射到二維平面的同一點時，主要解決方法是將曲面點集劃分多個子集，基于不同子集建立三角化網(wǎng)格后再進行拼接，可能導(dǎo)致出現(xiàn)各部分曲面邊界的拓?fù)潢P(guān)系不一致問題。確定映射平面、獲取點集包圍盒等子問題增加映射法實現(xiàn)難度。

本文提出一種改進的基于邊界控制-局部映射的曲面三角化網(wǎng)格生成算法，在三維空間中基于點集構(gòu)造三角化網(wǎng)格模型，直接剖分算法減少映射算法中數(shù)據(jù)預(yù)處理等計算過程，通過分析三維空間點位置關(guān)系生成三角化剖分模型。在存在邊界條件約束和無邊界條件約束下，均可以建立三維地質(zhì)界面三角化網(wǎng)格模型，客觀反映地質(zhì)界面展布形態(tài)，避免了三維空間中多點映射二維平面同一點的多對一映射問題。

2 三維空間曲面三角化網(wǎng)格模型

2.1 Delaunay三角化

Delaunay三角化特征能夠保證三角化網(wǎng)格模型具有良好的構(gòu)網(wǎng)形態(tài)[16]。Delaunay三角化的三個重要性質(zhì)描述如下[17]：

(1) 空外接圓(球)性質(zhì)。在二維空間平面Delaunay三角網(wǎng)中每個三角形的外接圓均不包含其他節(jié)點。三維空間中，每個四面體的外接球不包含其他節(jié)點。

(2) 最小角最大化性質(zhì)。滿足所有三角形中最小角的最大化性質(zhì)。

(3) 局部優(yōu)化性質(zhì)。三角網(wǎng)局部優(yōu)化過程能夠保證全局優(yōu)化。

在二維平面上可以證明能夠構(gòu)建Delaunay三角化網(wǎng)格模型，滿足Delaunay三角化特征。基于映射法，在二維平面構(gòu)建的Delaunay三角化網(wǎng)格模型映射到三維空間時，Delaunay三角化特征無法描述。但其幾何定義是保證三角網(wǎng)格質(zhì)量理論依據(jù)。因此可以將Delaunay三角化特征應(yīng)用于建立三維曲面三角化網(wǎng)格模型，控制三角網(wǎng)格形態(tài)。

2.2 曲面三角化模型構(gòu)建算法

基于Delaunay三角化特征提出一種邊界控制-局部映射的曲面三角化網(wǎng)格模型構(gòu)建算法，算法流程見圖1，主要包括三個步驟：(1) 三角化網(wǎng)格模型初始化；(2) 三角化網(wǎng)格模型更新；(3) 網(wǎng)格模型邊界約束更新。

圖1 曲面三角化網(wǎng)格模型構(gòu)建算法流程

斷層的斷點點集S={p1,p2,…,p25}作為實驗數(shù)據(jù)，詳細說明邊界控制-局部映射曲面三角化網(wǎng)格模型構(gòu)建算法步驟。假設(shè)算法輸入為點集S，算法輸出為邊集合E以及三角形網(wǎng)格模型T。

2.2.1三角化網(wǎng)格模型初始化

算法首先從三維點集合中挑選三個點構(gòu)成網(wǎng)格模型的初始化三角形。以初始三角形為中心，不斷加入新的三維點更新網(wǎng)格模型。曲面空間中三角形由三邊構(gòu)成，將邊定義為兩種類型：內(nèi)部邊和外界邊。在曲面三角網(wǎng)格模型T中：(1) 內(nèi)部邊：一條邊屬于兩個不同的三角形；(2) 邊界邊：一條邊屬于且只屬于一個三角形。算法需要區(qū)分并記錄邊類型，邊界邊集合表示為borderset，邊界邊集合內(nèi)所有邊構(gòu)成閉環(huán)。基于邊界邊拓展網(wǎng)格模型，基于內(nèi)部邊更新網(wǎng)格模型，同時保證三角化網(wǎng)格質(zhì)量。

采用隨機方式在點集中確定第一個三角形，該三角形的三條邊均屬于邊界邊，記錄邊界邊信息，圖2顯示了斷層斷點數(shù)據(jù)集中隨機初始三角形結(jié)果。

圖2 基于點集S初始化種子三角形示意圖

三角化網(wǎng)格模型隨機初始化三角形思路：(1) 在點集S中，獲取隨機點p9，根據(jù)最近鄰原則選取最近點p4，構(gòu)成第一條三角形邊e1=(p4，p9)，將p4、p9標(biāo)記成已訪問狀態(tài)，以標(biāo)記點集S′=(p4，p9)，點集S更新為S-S′。(2) 根據(jù)點到邊e1中點的最小距離確定三角形第三點。點p到邊e中點距離表示為dist(p,e)=dist(p,pm)，pm是e的中點。因此，p20作為距離邊e1中點最近點構(gòu)成第一個三角形，并將其標(biāo)記為已訪問，更新邊界邊集合borderset和三角化網(wǎng)格T：

borderset={e1=(p9,p4)，e2=(p4,p20)，e3=(p20,p9)}

T={t1=(p9,p4,p20)}

2.2.2三角化網(wǎng)格模型更新

三角網(wǎng)格模型更新過程為從點集S中選取未標(biāo)記的三維點加入到三角化網(wǎng)格模型中，更新邊界邊集合borderset和三角網(wǎng)模型T，直至點集S中所有點都加入到三角化網(wǎng)格模型，如圖3所示。

圖3 三角化網(wǎng)格模型更新流程

計算點集S中的點到borderset集合中邊界邊中點的最小距離。以圖2中斷點三角化網(wǎng)格為例，計算可得dist(p14,e1)≤dist(p′,e1)，dist(p15,e2)≤dist(p′,e2)，dist(p14,e3)≤dist(p′,e3)，?p′∈S，且dist(p14,e1)

1) 局部映射。假設(shè)曲面連續(xù)光滑，基于Delaunay三角化特征采用局部映射方法判斷點與三角形的位置關(guān)系，更新三角化網(wǎng)格模型，確保三角網(wǎng)格良好形態(tài)。

采用局部映射方法判斷空間一點與三角形t位置關(guān)系。假設(shè)任意三角形三點坐標(biāo)為p1={x1,y1,z1}，p2={x2,y2,z2}，p3={x3,y3,z3}，以其中一邊(p2，p3)為標(biāo)準(zhǔn)，過邊(p2，p3)垂直于三角形的平面方程A(x,y,z)，即：

式中：

通過垂直平面方程，以邊(p2，p3)為標(biāo)準(zhǔn)，將點p0={x0,y0,z0}和三角形第三點p1坐標(biāo)代入A計算可得A1和A2。任意一點p0與三角形t的空間位置關(guān)系定義為：

① 若A1×A2<0，p0位于以邊(p2，p3)為標(biāo)準(zhǔn)三角形的負(fù)向空間；

② 若A1×A2>0，p0位于以邊(p2，p3)為標(biāo)準(zhǔn)正向空間；

③ 若A1×A2>0，p0在三角形所在平面的垂直平面上，邊(p2，p3)也在該垂直平面上。

判斷一點p0在三角形p1p2p3平面的局部投影是否在該三角形范圍內(nèi)，以三邊(p1，p2)、(p2，p3)、(p1，p3)為標(biāo)準(zhǔn)，若p0均在各邊的正向空間中，則表示該點局部映射在三角形內(nèi)，反之，則表示局部映射不在三角形范圍內(nèi)。

圖4 三角形與點空間位置關(guān)系示意圖

根據(jù)三種空間位置關(guān)系，采用三種不同方法更新三角化網(wǎng)格模型：

(1) 點pt位于邊界邊e負(fù)向空間時，點pt與邊界邊e構(gòu)成新的三角形加入網(wǎng)格，邊界邊e變?yōu)閮?nèi)部邊，并增加邊界邊e1=(p2,pt1)和e2=(pt1,p3)。

在該種情況下，點pt在平面A上的投影與邊界邊關(guān)系可為三種，見圖5，局部投影點可能位于邊界邊e上、邊界邊e沿p2p3方向的延長線上、邊沿p3p2方向的延長線上。

圖5 邊界邊e與在平面A上點的局部映射點位置關(guān)系

在圖5中根據(jù)不同的位置關(guān)系生成不同三角形，更新邊界邊集合，如圖6所示。

圖6 根據(jù)邊界邊上映射點位置關(guān)系更新網(wǎng)格模型

圖7 邊界邊與其正向空間映射點關(guān)系

根據(jù)正向空間內(nèi)不同的位置關(guān)系采用不同方法生成新三角形：

圖8 邊界邊與前序邊界邊最近點相同情況下生成新三角形

圖9 邊界邊與前序邊界邊最近點不相同情況下生成新三角形

2) 最小角最大化原則?；贒elaunay三角化特征保證三角化剖分質(zhì)量，引入最小角最大化對三角網(wǎng)格模型進行適當(dāng)變化調(diào)整。在圖10中，點p4與邊界邊e組成新三角形t2=(p4,p3,p2)，共邊e的另一三角形為t1=(p1,p2,p3)，計算兩個三角形的最小角。對換對角線(p1,p4)，計算兩個三角形最小角值，與變換前最小角值比較，若變換后最小角值較大，則進行對角線變換，變換后三角形為t2=(p4,p3,p1)、t1=(p1,p2,p4)。當(dāng)三角形對角線產(chǎn)生對換后，同樣檢驗其相鄰三角形是否需要進行對角線兌換，保證最小角最大原則。

圖10 三維空間三角形對角線變換

保證最小角最大化原則，更新三角形網(wǎng)格模型主要流程為：

(1) 在網(wǎng)格模型中每加入一個新三角形，采用隊列q處理三角形對角線變換。對于每個新三角形，將其內(nèi)部邊放入隊列構(gòu)成隊列初始狀態(tài)。

(2) 依次從隊列取出一條內(nèi)部邊e，直至隊列q為空結(jié)束，否則，繼續(xù)執(zhí)行(3)。

(3) 計算內(nèi)部邊e所屬兩個三角形t1、t2，根據(jù)最小角最大化原則進行判斷。若需要對角線變換，則交換兩個三角形對角線，并將t1、t2中其他內(nèi)部邊加入到隊列q中，繼續(xù)執(zhí)行(2)。

以圖2中點集S為例，不斷加入未標(biāo)記點更新網(wǎng)格模型，當(dāng)網(wǎng)格模型中三角形個數(shù)為12時，三角網(wǎng)格模型如圖11所示。

圖11 點集S三角網(wǎng)格模型更新過程中示意圖

此時邊界邊集合borderset和三角網(wǎng)格模型T更新為：

borderset={(p9,p4)，(p4,p5)，(p5,p10)，(p10,p15)，(p15,p20)，(p20,p25)，(p25,p24)，(p24,p18)，(p18,p13)，(p13,p8)，(p8,p14)，(p14,p9)}

T={(p9,p4,p20),(p4,p15,p20),(p4,p10,p15),(p4,p5,p10),(p9,p20,p14)，(p14,p20,p25),(p14,p19,p25)，(p19,p24,p25),(p13,p24,p18),(p13,p19,p24)，(p13,p19,p8),(p3,p19,p14)}

當(dāng)點集S為空時結(jié)束三角化網(wǎng)格更新，記錄下邊界邊集合，下一步進行網(wǎng)格邊界約束更新。圖12展示圖2中所有點加入到三角化網(wǎng)格模型，顯然三角化網(wǎng)格模型構(gòu)成凹多邊形，進一步更新網(wǎng)格模型邊界，采用凸多邊形模擬三維曲面邊界。

圖12 所有點加入三角化網(wǎng)格模型示意

2.2.3網(wǎng)格模型邊界約束更新

根據(jù)邊界約束更新三角化網(wǎng)格模型，主要有兩種方法：(1) 無邊界數(shù)據(jù)約束情況下，自動填充邊界，原則上自動填充為凸多邊形邊界。其基本思路是按照邊界邊集合中邊界連接順序，沿著邊界邊順序逐一更新模型，直至邊界邊無需更新操作為止。(2) 根據(jù)已有的點集S曲面邊界數(shù)據(jù)約束更新三角化網(wǎng)格模型。

在邊界邊集合中任意邊界邊e都有與其相連的前邊界邊ebefore與后邊界邊enext。在有曲面邊界數(shù)據(jù)約束的情況下，若邊界邊e被確定為曲面邊界，則無需處理；反之，則對當(dāng)前邊界邊進行更新，判斷邊界邊e能否與ebefore或是enext構(gòu)成新三角形，更新模型邊界。根據(jù)三邊關(guān)系更新網(wǎng)格模型邊界基本思路如下：

① 若ebefore和e同屬于一個三角形：ebefore和e無法構(gòu)成新三角形，同樣，若e和enext同屬于同一個三角形，則e和enext無法構(gòu)成新三角形。

② 若ebefore和e不同屬于一個三角形：如果沒有點的局部投影在邊ebefore和e所在平面的三角形內(nèi)，則ebefore和e可以構(gòu)成三角形，反之，則無法構(gòu)成三角形。采用同樣方法判斷enext和e是否能夠構(gòu)成三角形。

③ 若只有邊ebefore和e能構(gòu)成三角形，將新三角形加入到網(wǎng)格模型，刪除邊界邊ebefore和e,將新三角形第三邊加入邊界邊集合繼續(xù)邊界更新操作；若只有enext和e能夠構(gòu)成新三角形，加入網(wǎng)格模型，刪除邊界邊enext和e，第三邊加入到邊界邊集合；若ebefore和e，enext和e都能構(gòu)成新三角形，根據(jù)兩個三角形的最小角，選擇最小角最大三角形加入到網(wǎng)格模型中。每增加一個新三角形到網(wǎng)格模型時，基于最小角最大化原則，局部更新網(wǎng)格模型，保證網(wǎng)格形態(tài)。

在無曲面邊界約束數(shù)據(jù)的情況下，圖2中點集S進行模型邊界更新后，三角化模型展示如圖13所示。

圖13 點集合S曲面三角網(wǎng)格模型

3 地質(zhì)界面三角化網(wǎng)格模型實例

地質(zhì)界面數(shù)據(jù)常采用點集合數(shù)據(jù)表示，例如斷層采用斷點數(shù)據(jù)描述，地層采用地層點數(shù)據(jù)描述。本文提出一種基于點集數(shù)據(jù)的局部映射-邊界控制三角化網(wǎng)格模型構(gòu)建算法。

以斷點數(shù)據(jù)集為例說明算法實現(xiàn)流程與效果。在三角化網(wǎng)格模型構(gòu)建算法流程中，隨機初始化三角形網(wǎng)格模型，依次加入斷點集合中其他點數(shù)據(jù)，更新三角化網(wǎng)格模型，直至斷點數(shù)據(jù)均加入到網(wǎng)格模型結(jié)束更新，其構(gòu)成的三角化網(wǎng)格模型如圖14所示。

圖14 加入所有點數(shù)據(jù)后生成的三角化網(wǎng)格模型

基于圖14中的三角化網(wǎng)格模型和邊界邊集合，更新曲面三角化網(wǎng)格模型邊界。在無曲面邊界數(shù)據(jù)約束的情況下，填充邊界三角形，根據(jù)邊界邊集合自動生成凸多邊形曲面，如圖15(a)所示。當(dāng)加入邊界約束數(shù)據(jù)時，恢復(fù)模型邊界，如圖15(b)所示。

(a) 無邊界約束邊界自動更新 (b) 加入邊界約束數(shù)據(jù)模型圖15 某斷層曲面三角化網(wǎng)格模型

圖15所顯示的斷層三角化曲面網(wǎng)格模型中，存在兩點距離較近和兩點映射到某平面共線等情況。若采用映射法，選擇不合適的映射平面會導(dǎo)致三角化結(jié)果錯誤，映射平面直接影響網(wǎng)格模型生成質(zhì)量。直接在三維空間中通過空間位置判斷生成三維曲面三角化網(wǎng)格模型，能夠較好地還原曲面的真實形態(tài)。雖然空間位置判斷較為復(fù)雜，實際上，在算法實現(xiàn)過程中減少數(shù)據(jù)預(yù)處理和其他的計算過程。

圖16展示多個斷層面三角化網(wǎng)格模型形態(tài)，在生成多斷層面三角化網(wǎng)格模型時，算法優(yōu)點明顯，與映射法不同，不需要計算每個斷層面的二維映射平面，而是直接在三維空間中構(gòu)建多個斷層面的網(wǎng)格模型，減少算法其他輔助計算步驟。

圖16 多個斷層面三角化網(wǎng)格模型

與經(jīng)典映射法相比，本文提出的曲面三角化網(wǎng)格模型構(gòu)建算法能夠直接應(yīng)用于構(gòu)造較復(fù)雜的三維曲面模型。圖17、圖18展示了基于點集數(shù)據(jù)采用局部映射-邊界控制算法生成的較為復(fù)雜三維曲面模型，經(jīng)典映射法較難實現(xiàn)，原因在于很難確定一個映射平面，將三維數(shù)據(jù)點映射到二維平面上，進而較難生成準(zhǔn)確的三維曲面三角化模型。

(a) 斜視圖

(b) 俯視圖

(c) 側(cè)視圖圖17 復(fù)雜三維曲面三角化網(wǎng)格模型示意圖

(a) 斜視圖

(b) 側(cè)視圖

(c) 俯視圖圖18 復(fù)雜三維曲面三角化網(wǎng)格模型示意圖

4 結(jié) 語

在三維地質(zhì)建模中，基于地質(zhì)構(gòu)造模型研究地質(zhì)界面和地質(zhì)體構(gòu)造形態(tài)具有重要的意義，構(gòu)造三維地質(zhì)界面三角化網(wǎng)格模型作為地質(zhì)體構(gòu)造約束，是研究地質(zhì)構(gòu)造的基礎(chǔ)。本文基于現(xiàn)有曲面三角化網(wǎng)格構(gòu)建算法研究分析，結(jié)合地質(zhì)數(shù)據(jù)特點，提出一種局部映射-邊界控制的三維三角化網(wǎng)格模型構(gòu)建算法。基于點集數(shù)據(jù)構(gòu)建地質(zhì)界面三角網(wǎng)格化曲面模型，模擬地質(zhì)界面特征。采用斷層點集數(shù)據(jù)驗證算法有效性，可視化展示斷層三角化網(wǎng)格模型。與經(jīng)典映射法相比，局部映射-邊界控制三角化網(wǎng)格模型構(gòu)建算法能夠在三維空間中直接基于點集數(shù)據(jù)建立曲面三角化網(wǎng)格模型，無需對點集數(shù)據(jù)進行預(yù)處理，避免映射法三維空間中多點映射到二維平面一點所引起問題，并且能夠保證較好的網(wǎng)格模型形態(tài)。該算法可以拓展到其他領(lǐng)域用于建立三維曲面網(wǎng)格模型。