三維復雜斷層建模關鍵算法研究

侯曉琳 強偉帆

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

0 引 言

斷層是一種重要的地質構造,是控制油田剩余油分布的主要因素,也有可能是油氣滲流的通道,斷層模型的準確性直接影響油藏數值模擬的結果[1-2]。斷層建模是地質結構建模和屬性建模的基礎,由斷層、地層等地質界面封閉包圍成地質實體,三維斷層模型的準確性影響地質結構建模和屬性建模質量[3-5]。斷層模型可以描述斷層構造的幾何形態,斷層之間的接觸關系,可以描述斷層之間的關聯關系,構建斷層模型的幾何拓撲關系實現幾何元素之間的快速索引能大大提高斷層建模處理速度[6-7]。

斷層建模主要包括生成斷層網格模型和描述斷層接觸關系?；跀鄬狱c數據能生成并優化三維斷層網格化模型,真實反映斷層展布情況。在復雜斷層網中,當斷層傾向不一致時可能出現斷層交叉、重疊、截斷等情況,因此會形成Y型[8]、X型[9]、半Y型、半λ型斷層,構建這些斷層模型過程復雜[10]。目前復雜斷層建模方法主要是基于初始化斷層面網格模型,采用人工交互的方式描述斷層接觸關系,生成復雜斷層網格化模型,建模效率較低,并且在應用復雜斷層模型建立三維地質體模型時受到一定條件的限制[11-12]。

為改進現有復雜斷層建模存在的問題,提出一種自動判斷斷層接觸關系的算法,基于先驗地質斷層類型,自動生成復雜斷層模型,減少了人工交互工作量,保留了幾何拓撲關系一致性特征,提高了斷層建模和三維地質建模效率。

1 三維復雜斷層建模

國內外學者對三維地質建模研究較多[13],三維地質結構建模方法主要包括鉆孔建模、基于三棱柱建模、結構化網格建模、隱函數建模等[14-16]。根據不同的三維地質建模思路,采用不同的斷層建模方法,但是描述復雜斷層形態存在一定難度。比如,基于Pillar網格的斷層建模無法準確生成復雜斷層模型[17]。因此,三維復雜斷層建模研究具有意義和應用價值。

三維復雜斷層建模流程主要包括斷層網格模型生成和復雜斷層接觸關系描述兩部分(圖1)?；跀鄬狱c集數據構建三維斷層三角化網格模型,保證幾何拓撲關系一致性,作為初始化斷層面模型。進而控制生成的斷層面網格模型質量,優化并加密網格模型。在構建簡單斷層網格模型基礎上,判斷斷層接觸關系,生成斷層接觸關系描述,對具有接觸關系的復雜斷層進行優化和網格加密,改進網格質量?；谙闰灥牡刭|知識斷層類型,自動生成復雜斷層網格模型,提高復雜斷層建模效率。

圖1 三維復雜斷層建模流程

2 三維斷層三角化網格模型構建方法

2.1 三維斷層三角化網格模型

構建斷層三角化網格模型準確描述斷層形態,結合映射法和二維三角化網格剖分算法[18-19],提出一種基于斷層點集數據的三維斷層三角化網格模型算法(圖2)。

圖2 基于斷層點集數據構建三維斷層面網格模型流程

首先基于三維斷層點集的質心和擬合平面法向量計算斷層點集的擬合平面。假設斷層點集S包含n個三維點S={p1,p2,…,pn},點pi坐標為(xi,yi,zi)。計算點集S的質心點pc=(xc,yc,zc)。

(1)

計算點集S擬合三維平面P：w1x+w2y+w3z+w4=0的法向量,其中w表示待求參數。點集S中的所有點構成n×3矩陣An×3。

(2)

采用奇異值分解方法[20],將矩陣An×3分解為An×3=Un×nΣn×3V3×3,Un×n=[u1,u2,…,un],ui表示Un×n的列向量,對角矩陣Σn×3表示為：

(3)

式中：σ1、σ2、σ3為奇異值。V3×3=[v1,v2,v3],vi為V3×3列向量。最小奇異值σ3對應的列向量v3就是平面P的法向量,即v3=[w1,w2,w3]T。根據質心坐標,計算w4,w4=-w1xc-w2yc-w3zc。根據平面P方程和空間坐標轉換將三維平面映射為二維平面。平面P上三維點pi=(xi,yi,zi)映射為二維點的轉換矩陣M。

(4)

將三維斷層面上的點一一映射到二維平面后,采用二維三角化剖分算法,建立二維Delaunay三角化網格模型T′,算法輸入為三維點集P在二維平面上的映射點集P′,算法輸出為三角化網格單元模型T。二維Delaunay三角化網格模型構建算法流程：

1) 計算二維平面上的點集P′的外包圍框abcd,選取P′中的一點p0,p0與外包圍框abcd構成三角形abp0、bcp0、cdp0、adp0,作為初始化三角化網格模型T′。

(a) (b)圖3 加入點生成新三角形更新網格模型

4) 根據三維點與二維點的一一映射關系,構建并復原三維斷層三角化網格模型T(圖4)。

圖4 基于映射關系構建三維斷層三角化網格模型

基于斷層點集數據,結合映射法和二維三角化網格剖分算法建立三維斷層網格模型,準確描繪斷層在三維空間中的展布情況,圖5展示了基于斷層點集生成的三維斷層網格化模型。

2.2 三維斷層面網格模型優化

基于三維斷層網格化模型,維護幾何拓撲關系,優化網格形態,進一步改進三維斷層模型。三維斷層網格模型優化主要包括畸形網格單元調整和網格加密。

幾何拓撲關系模型是三維斷層建模的基礎,能夠快速索引斷層模型中的幾何元素。圖6是三維斷層網格模型基本數據結構,維護幾何拓撲關系,能夠更加快速索引點、線、面數據,加快網格優化等計算速度。

圖6 三角化網格拓撲關系數據結構設計

基于斷層點集數據生成的三維斷層面網格模型可能存在畸形網格單元等問題,采用網格優化方法改進三維斷層面網格模型形態。三角形單元網格質量評價指標有很多種,本文主要采用兩種評價尺度：相鄰三角形最小內角、三角形外接圓半徑和最短邊的比。通過增加數據點控制網格單元質量,優化網格單元形態。

1) 相鄰三角形最小內角：基于網格幾何拓撲關系采用局部調整方法優化兩個相鄰三角面結構形態。圖7中,相鄰三角形t(abc)和t(acd),計算兩個三角形最小角,變換對角線bd,計算兩個三角形t(abd)和t(bcd)最小角,基于最小角最大化原則選擇最小角最大兩個三角形加入網格模型。

圖7 相鄰三角形最小角最大化局部調整

最小角最大化局部調整算法輸入包括斷層三角形網格模型T,算法輸出包括局部優化后的斷層三角形網格模型T,算法流程：

(1) 依次遍歷網格模型T中的每個三角形t,直至相鄰三角形無須進行最小角最大化局部調整為止算法結束。否則,根據幾何拓撲關系,快速查找三角形t共邊鄰接的三角形ts1、ts2、ts3。

(2) 依次判斷三角形ts1、ts2、ts3與t的空間位置關系。假設tsi(bcd)和t(abc)共邊bc,計算兩個三角形的最小內角α,對換對角線ad,生成兩個新三角形t(abd)和t(acd)。

(3) 計算三角形t(abd)所在三維平面S和點c在平面S上的投影點c′。若點c′在三角形t(abd)內(圖8(a)),無法對換對角線ad,繼續執行(1)。若點c′在三角形t(abd)外(圖8(b)),繼續執行(4)。

(a) (b)圖8 相鄰三角形對換對角線示意圖

(4) 計算三角形t(abd)和t(acd)的最小內角α′,若α′>α,將三角形t(abd)和t(acd)加入到網格模型T,刪除三角形tsi(bcd)和t(abc),更新幾何拓撲關系,繼續執行(1)。

2) 三角形t外接圓半徑和最短邊的比ρ：

(5)

(1) 依次遍歷三角化網格模型T中的每個三角形t,計算三角形t的外接圓與最短邊比ρ,最長邊長l和最短邊長l′。

(3) 若ρ>ρmax,t為畸形網格單元,添加三角形t的最長邊中點,構建新的三角形,更新幾何拓撲關系。若三角形t最長邊只屬于三角形t,則生成兩個新三角形;若三角形t最長邊分屬于兩個不同的三角形,則生成四個新三角形,加入到網格模型T,刪除畸形網格單元。繼續執行(1)。

(a) (b) (c)圖9 不同閾值斷層局部網格優化加密效果

優化并加密后的斷層網格模型作為復雜斷層建模的基礎模型,根據斷層接觸關系,進一步進行復雜斷層建模。

3 復雜斷層自動化建模方法

3.1 斷層接觸關系處理

斷層接觸關系判斷和處理是復雜斷層建模的基礎,根據接觸關系能夠判斷斷層類型,結合先驗地質知識,實現自動化復雜斷層建模。

基于三維斷層網格模型,建立三維空間柵格模型,判斷三維斷層網格模型T1和T2接觸關系。假設。斷層面T1和T2的外包圍框為B1(a1b1c1d1-e1f1g1h1)和B2(a2b2c2d2-e2f2g2h2),表示為：

式中：xmin、xmax分別表示斷層面T1中所有斷層點最小x坐標和最大x坐標(圖12)。柵格模型的大小是兩個斷層模型的并外包圍框B：

圖12 兩個斷層模型T1和T2外包圍框示意圖

B=B1∪B2={xmin,xmax,ymin,ymax,zmin,zmax}

圖13 三角形t所在柵格單元示意圖

(6)

(7)

(8)

記錄三角形t(abc)與所在柵格單元的對應關系,作為判斷斷層接觸關系算法輸入。斷層接觸關系判斷方法簡化為判斷兩個三角化網格模型相交關系。

在某一個柵格單元cell(n,m,k)中的任意兩個三角形t1∈T1和t2∈T2可能相交。遍歷N×M×K個柵格模型中的所有柵格單元cell(n,m,k),若在同一個cell(n,m,k)中包含任意兩個三角形t1(a1b1c1)∈T1和t2(a2b2c2)∈T2,進一步判斷兩個三角形是否相交并計算交點。判斷任意兩個三角形t1(a1b1c1)∈T1和t2(a2b2c2)∈T2相交關系和計算交點的算法輸入包括三角形t1(a1b1c1)和t2(a2b2c2),算法輸出包括是否相交判斷結果,若相交,輸出三角形相交交點。算法流程：

1) 計算三角形t1(a1b1c1)所在平面方程S1((x,y,z,1)T)：(w1,w2,w3,w4)·(x,y,z,1)T=0,將t2(a2b2c2)三點代入平面方程S1得到S1(a2)、S1(b2)、S1(c2)。若S1(a2)、S1(b2)、S1(c2)同為正值或者負值表示兩個三角形不相交,算法結束;反之,執行2)。

3) 三角形t1(a1b1c1)和t2(a2b2c2)相交,計算交點。如圖14所示,依次計算t3-i(a3-ib3-ic3-i),i=1,2三邊與三角形ti(aibici)交點。以邊a3-ib3-i為例,若Si(a3-i)×Si(b3-i)<0,點a3-i、b3-i在平面Si兩側,計算線段a3-ib3-i與平面Si的交點I。若交點I在三角形ti(aibici)內,則交點I就是兩個三角面t1和t2的交點,記錄兩個相交三角形以及t3-i(a3-ib3-ic3-i)的一邊a3-ib3-i。保存三角形與另一三角形邊的相交關系。

圖14 兩個相交三角形位置關系及交點示意圖

在計算出兩個斷層網格模型所有交點的基礎上,采用ear clipping算法基于交點交線分割相交三角形,保留交點、交線的同時保證幾何拓撲關系一致性。ear clipping算法流程主要包括構建三角形閉合多邊形和三角形分割算法兩個部分[21]。斷層網格模型中一個三角形ti(aibici)與另一個斷層模型中的三角形相交,假設交點集合I={I1,I2,…},首先需要基于交點構建三角形ti(aibici)的閉合多邊形。三角形ti(aibici)閉合多邊形Po構建算法輸入包括三角形ti(aibici)和ti(aibici)上交點集合I,算法輸出為三角形ti(aibici)閉合多邊形Po。算法流程：

1) 基于坐標軸對交點集合I進行排序,假設排序后I={I1,I2,…,In},那么交線線段集合為L={I1I2,I2I3,…,In-1In}。

2) 計算三角形ti(aibici)三邊與交線線段L構造閉合多邊形Po?；诮痪€線段構建閉合內環為I內={I1,I2,…,In,In-1,…,I3,I2}?；谌切稳厴嫿ㄩ]合外環。若任意交點Ii位于三角形邊上,將交點Ii加入到閉合外環。根據三角形ti(aibici)與交線線段的位置關系不同,三角形ti(aibici)閉合多邊形Po構建方法不同,分成以下三種情況：

(1) 交線兩端與均在三角形ti邊上(圖15(a))：閉合外環I外={ci,I1,ai,bi,In},基于連接點I1連接閉合外環和閉合內環,構建閉合多邊形Po={ci,I1,I2,…,In,In-1,…,I2,I1,ai,bi,In}。

(2) 交線一端在三角形ti一邊上,其余交點在三角形內(圖15(b))：閉合外環I外={ci,I1,ai,bi},基于連接點I1連接閉合外環和閉合內環,構建閉合多邊形Po={ci,I1,I2,…,In,In-1,…,I2,I1,ai,bi}。

(3) 交線點均在三角形ti內(圖15(c))：閉合外環I外={ci,ai,bi},建立內環與外環連接線段aiI1,基于連接線段,閉合外環和閉合內環構建閉合多邊形Po={ci,ai,I1,I2,…,In,In-1,…,I2,I1,ai,bi,In}。

基于三角形ti頂點與交點構成閉合多邊形Po,依次遍歷并刪除閉合多邊形上的點生成新三角形,直至分割生成所有的新三角形?；陂]合多邊形Po分割三角形算法輸入閉合多邊形Po,算法輸出分割后的新三角形。算法流程：

1) 若ti(aibici)的閉合多邊形Po中有且僅有三個點pi-1、pi、pi+1,則pi-1、pi、pi+1構成新三角形;否則,執行2)。

2) 遍歷ti(aibici)的閉合多邊形Po中一點pi∈P其前繼點表示為pi-1∈P,后繼點為pi+1∈P。若pi-1、pi、pi+1共線,無法構成三角形,繼續遍歷Po中下一點,執行2);否則,執行3)。

3) 若pi-1、pi、pi+1不共線,則pi-1、pi、pi+1可以構成三角形。判斷Po中其他點是否在三角形t(pi-1pipi+1)內,若其他點在三角形內,pi-1、pi、pi+1無法構成新三角形,繼續執行2);反之,pi-1、pi、pi+1構成新三角形t(pi-1pipi+1),在閉合多邊形Po中刪除當前點pi。

(a) (b)圖16 ear clipping算法分割三角形

圖17展示了分屬于不同斷層模型的兩個相交三角形,采用ear clipping算法分割新三角形結果。圖18展示了兩個三角化網格模型相交關系和相交處理結果。圖19展示了基于兩個斷層模型接觸關系處理后的網格模型。

圖17 采用ear clipping算法分割兩個相交三角形結果示意圖

圖18 兩個三角化網格模型相交處理示意圖

圖19 基于兩個斷層模型接觸關系處理后的網格模型

3.2 復雜斷層網格模型優化

采用相交線段描述斷層網格模型的接觸關系,構建復雜斷層網格模型,在交線約束基礎上,優化三角化網格模型形態,加密網格模型。

優先優化以交線為邊或交點為頂點的三角形網格單元形態,準確描述斷層模型接觸關系。在任意一個相交三角形t(abc)中,遍歷所有分割生成的新三角形ti(aibici)的每一條邊e,若邊e不屬于交線線段L且共屬于兩個不同三角形,依據最小角最大原則判定是否對換對角線生成并更新新三角形。圖20中,保留交線線段I1I2、I2I3、I3I4,對換對角線I1b、I1I3,生成新三角形t(aI1I2)、t(abI2)、t(I1I2c)、t(I2I3c),刪除原有三角形t(I1ab)、t(bI1I2)、t(I1I2I3)、t(I1I3c)。

(a)

圖21展示了基于相交線段約束,相交網格模型優化前、后的接觸關系的不同描述結果。

(a)

2) 依次遍歷三角化網格模型Ti中的每個三角形t,計算三角形t的外接圓與最短邊比ρ,最長邊長l和最短邊長l′。

4) 若ρ>ρmax,添加三角形t的最長邊中點m,若三角形t的最長邊不屬于相交線段,構建新三角形,更新幾何拓撲關系。若三角形t的最長邊屬于相交線段e(IiIi+1),執行5)。

圖22展示了兩個相交的三角化網格模型經過優化、加密后更新的網格模型。圖23展示了具有接觸關系的復雜斷層網格模型優化加密后的可視化效果。

圖22 兩個相交三角化網格模型優化加密后的網格模型

圖23 具有接觸關系斷層面網格模型優化加密結果示意圖

3.3 基于地質先驗知識生成復雜斷層模型

在現有地質斷層先驗知識的基礎上,根據斷層和分支斷層關系,基于三維斷層網格模型接觸關系,自動生成三維復雜斷層網格模型。假設兩個相交三維斷層網格模型T1、T2,基于先驗知識已知T1為主斷層,T2為分支斷層,主斷層模型完全截斷分支斷層模型,基于兩個斷層網格T1、T2構建復雜斷層網格模型,描繪復雜斷層展布情況。首先基于主斷層T1、分支斷層T2,計算并生成分支斷層T2的分割線段,分割線段將分支斷層T2截斷成兩個子網格模型,其中之一面積較大網格模型表示分支斷層T2實際展布情況。在分支斷層沒有被主斷層完全截斷的情況下,需基于T1、T2的相交線段L={I1I2,I2I3,…,In-1In}計算相交線段的延長線段形成分割線段。根據不同的接觸關系,分割線段主要包括以下三種情況：

1)T1、T2的相交線段L就是截斷T2的分割線段(圖24(a))。

(a) (b) (c)圖24 不同接觸關系的主斷層與分支斷層

2)T1、T2的相交線段L及一端延長線組成分割線段截斷T2(圖24(b))。

3)T1、T2的相交線段L及兩端延長線組成分割線段截斷T2(圖24(c))。

基于圖24所示的不同接觸關系的主斷層截斷分支斷層,構成復雜斷層模型(圖25)。

(a) (b) (c)圖25 不同接觸關系的主斷層截斷分支斷層示意圖

(a)

根據地質斷層先驗知識,判斷是否需要生成斷層T2截斷線段,截斷T2斷層,生成復雜斷層模型。若需要截斷T2斷層,采用先進先出隊列,遍歷三角形的鄰接關系,當三角形的一邊屬于截斷線段時,與其共邊三角形屬于截斷后的兩個不同子網格模型,反之,則與其共邊三角形屬于同一個網格模型,最后斷層T2被截斷成兩個子網格模型。針對于Y型斷層來說,刪除截斷后網格面積較小的子網格模型,構成Y型斷層網格模型。

圖27 分支斷層模型的截斷線段計算過程

(a)

4 結 語

斷層建模是三維地質建模的基礎,基于斷層點集數據構建三維斷層網格化模型,能夠真實反映斷層展布情況。斷層之間的接觸關系復雜,采用人工交互方法構建復雜斷層模型,工作量較大,建模時間較長,本文提出一種基于斷層點集數據,根據已有的地質斷層類型先驗知識,自動構建復雜斷層網格模型的方法,準確還原復雜斷層構造形態,減少人工交互工作量,提高三維斷層建模效率。