牛露佳,王雙威,曾義文,朱晨媛,王占剛
中國礦業(yè)大學(北京)地球科學與測繪工程學院,北京,100083
內(nèi)容提要:斷層建模是三維地質結構建模中的主要過程之一。斷層面建模過程中需要根據(jù)斷層間的空間接觸或者切割關系進行幾何曲面的裁剪,目前方法利用三角網(wǎng)求交算法進行曲面裁剪,但是該算法處理復雜斷層面切割關系時往往不穩(wěn)定。筆者等提出了基于規(guī)則網(wǎng)格的復雜斷層網(wǎng)絡處理與自動化建模的方法和流程,詳細討論了基于網(wǎng)格化的斷層網(wǎng)絡模型形式化理論表達、建模流程中的斷層網(wǎng)絡空間關系構建以及相交裁剪處理算法等核心步驟。利用測試數(shù)據(jù)和煤礦三維地震構造解釋數(shù)據(jù)進行了驗證,表明該方法可以有效處理多條互相切割、主輔關系復雜的斷層網(wǎng)絡,具有較好的算法穩(wěn)定性;與SKUA—GOCAD斷層建模方法對比,能夠減少交互過程,提高斷層建模的自動化程度。
三維地質結構建模(Structural Geological Modeling)用于構建結構模型來表達地下地質對象的幾何形態(tài)及對象間的空間關系(Myers, 2003;Caumon et al., 2009;Wang Zhangang et al., 2016;Bi Zhengfa et al., 2022)。結構建模需要處理地質構造的接觸關系以及切割關系,其中,斷層建模主要解決斷層面的構建以及切割關系處理問題。復雜斷層網(wǎng)絡處理切割關系主要涉及X型,Y型、λ型、T型、交叉型、半Y型和半λ型等斷層間的組合方式(朱良峰等,2008;李兆亮等,2015)。
斷層建模的一般流程是利用斷層控制點或者斷層棱邊數(shù)據(jù)進行空間插值生成斷層面,然后根據(jù)不同斷層間的空間接觸或者切割關系進行幾何曲面的裁剪。斷層面采用不規(guī)則三角網(wǎng)(TIN,Triangulated Irregular Network)表達,利用三角形求交進行曲面裁剪。這些算法需要顯式處理斷層網(wǎng)絡中主輔斷層的削截關系來提高斷層建模的自動化程度,如基于測地線的路徑切割算法(李兆亮,2015),斷層的錯動和錯層處理(劉光偉等,2018)。在隱式地質建模(Renaudeau et al., 2019;Irakarama et al., 2021;de Kemp, 2021;Grose et al., 2021)中,斷層的處理不依賴于三角網(wǎng)的求交,比如GemPy(de la Varga et al., 2019)以及Grose等(2021)根據(jù)斷層兩盤的運動學特征處理斷層,這類方法需要斷層的幾何形態(tài)和運動學的數(shù)學描述,相關文獻中呈現(xiàn)的斷層以及切割關系相對簡單。從以上研究可以看出,基于不規(guī)則三角網(wǎng)的斷層面求交裁剪是復雜斷層網(wǎng)絡間切割處理與建模的主要手段。不規(guī)則三角網(wǎng)求交一般采用計算幾何中的三角形求交算法以及改進算法(Moller, 1997;趙景昌等,2014),但是在實際應用中三角網(wǎng)求交算法穩(wěn)定性差,直接影響著斷層建模的自動化處理能力,尤其針對多期構造或者斷層網(wǎng)絡的處理(蔣錢平等,2008;Pellerin et al., 2014)。
不采用三角網(wǎng)求交算法進行復雜斷層面切割裁剪處理是斷層建模仍需研究和解決的問題。為此,表達三維地質結構的網(wǎng)格化模型與建模方法被提出以解決地質對象切割、剖切分析等方面的技術難點,比如GTP(吳立新等,2003),規(guī)則網(wǎng)格(Graciano et al., 2018),Cut-Cell網(wǎng)格(Mallison et al., 2014 ),角點網(wǎng)格(于瑞濤,2018)和廣義棱柱網(wǎng)格(GPG)(Li Xin et al., 2018 )?;诰W(wǎng)格化模型的切割算法可以有效實現(xiàn)地質體的剖切處理并取得了較好的效果,如李長春等(2008)基于GTP實現(xiàn)三維地質體的剖切,Zhou Cuiying等(2020)提出垂直投影三角形網(wǎng)格(VPTN)實現(xiàn)地質體的剖切。這類方法的特點是根據(jù)地質分層特性,地質界面的分片構成單元(四邊形或者三角形)在空間XY方向上具有相同的離散結構,曲面分片單元頂點都在三棱柱或者四棱柱邊上,分片單元的切割裁剪只需要在Z方向上進行求交處理,消除了幾何圖形求交算法的誤差影響,保證了算法的穩(wěn)定性(Zhou Cuiying et al., 2020)。研究人員也在探索利用網(wǎng)格化方法進行斷層建模,比如基于GTP的斷層交互建模(車德福等,2008),基于角點網(wǎng)格的局部多級加密(于瑞濤,2018)和六面體元角點移位(劉少華等,2022)。然而,這些方法以交互為主,能自動處理的斷層比較簡單。因此,自動化處理復雜斷層網(wǎng)絡的切割裁剪以及斷層建模仍是一個問題。
由于網(wǎng)格化模型在切割算法穩(wěn)定性方面的優(yōu)勢,筆者等提出基于網(wǎng)格化的斷層網(wǎng)絡模型形式化理論表達,基于此提出基于規(guī)則網(wǎng)格的斷層網(wǎng)絡建模流程,詳細討論了建模流程中斷層空間關系構建以及斷層相交裁剪處理方法。
基于網(wǎng)格化的斷層網(wǎng)絡模型形式化表達為FaultModel={Ω,F,R},其中,定義域Ω={Ωj}表示模型整體空間XY方向的網(wǎng)格化離散域,Ωj是離散單元,類型可以是三角形、四邊形等,j=1...m,m是離散單元數(shù)量;斷層集合F={Fi},i=1...l,l表示斷層的數(shù)量;斷層空間關系集合R={rij},rij=
(1)斷層面Fi是單值曲面,Fi的分段函數(shù)表達為:
u=(x,y),Fi的定義域Ω1∪Ω2∪…∪Ωn?Ω。顯函數(shù)zi(x,y)為離散單元上的插值函數(shù),描述斷層Fi在(x,y)位置處的z值。
(2)斷層空間關系rij定義兩個斷層面是否相交的拓撲關系以及兩個斷層面Z方向上的方向關系。
可見,由于關系rij定義了空間上下關系,不能處理X型斷層組合。根據(jù)此關系可以得到如下推論:
推論1:Fi位于Fj的上方且Fj位于Fk的上方,則Fi位于Fk的上方。即:
zi(x,y)≥zj(x,y)且zj(x,y)≥zk(x,y),則zi(x,y)≥zk(x,y)
推論2:Fi位于Fj的下方且Fj位于Fk的下方,則Fi位于Fk的下方。即:
zi(x,y)≤zj(x,y)且zj(x,y)≤zk(x,y),則zi(x,y)≤zk(x,y)
Fi與Fj的關系rij根據(jù)實際地質情況,利用Fi與Fj的削截主輔關系以及形成先后順序進行設置。
若Fj是主斷層,Fi是輔斷層,則Fj切割Fi,Fi和Fj相交且Fi位于Fj上方或Fi位于Fj下方。
如果離散單元采用三角形單元對定義域Ω進行網(wǎng)格化離散,地質模型表達和建模類似主TIN法(李元亨等,2010)或者GTP建模(吳立新等,2003;車德福等,2008)。
將定義域Ω在XY方向上按一定分辨率均勻離散形成的規(guī)則網(wǎng)格稱為基網(wǎng)格,基網(wǎng)格單元稱為單元柱。如圖1,斷層面在基網(wǎng)格控制下進行網(wǎng)格化,在單元柱內(nèi)中形成四邊形對象稱為Interface。斷層面F由Interface拼合而成,即F={Interfaceij}={
圖1 斷層F網(wǎng)格化 (a)和單元柱(b)Fig.1 Fault F in regular grid(a) and Pillar(b)
通過基于網(wǎng)格化的斷層網(wǎng)絡模型形式化表達和Interface定義,斷層建模實質上將斷層面按網(wǎng)格化生成后進行切割處理,核心過程是在基網(wǎng)格的單元柱內(nèi)根據(jù)斷層空間關系進行Interface的求交和裁剪處理。由于單元柱內(nèi)不同Interface的差異是z值,求交和裁剪主要根據(jù)z值大小關系進行算法設計,有利于算法穩(wěn)定性。斷層空間上下關系決定了Interface間的上下關系,通過空間關系二叉樹(見3.1)管理所有斷層關系,達到自動化處理的目的。
復雜斷層網(wǎng)絡建模處理流程如圖2所示,主要步驟包括:
(1)基網(wǎng)格設置,根據(jù)建模源數(shù)據(jù)區(qū)域范圍和建模網(wǎng)格分辨率將定義域在XY方向上均勻離散形成的基網(wǎng)格。
(2)斷層網(wǎng)絡空間關系構建,根據(jù)實際地質情況,利用斷層的削截主輔關系構建斷層空間關系,包括空間上下關系和相交關系。通過空間關系二叉樹序列(見3.1節(jié))記錄所有斷層空間上下關系,相交關系表記錄斷層相交關系。
(3)網(wǎng)格化斷層面生成,逐個斷層處理,根據(jù)當前斷層F的區(qū)域范圍,利用薄板樣條插值算法或者有限差分插值算法(Irakarama et al., 2021)生成Interface拼合的網(wǎng)格化斷層面,即F={Interfaceij}。同時,使用主成分分析法計算斷層輸入控制點的最佳投影面,在最佳投影面上提取斷層邊界的真實最小凸包并進行調整作為斷層真實邊界條件。
(4)斷層邊界延伸預處理,對空間關系二叉樹中存在相交關系的斷層將切割主斷層面的區(qū)域范圍外延至所有被切割輔斷層區(qū)域范圍的最小公共外包,并對延伸部分進行標記(見3.2節(jié))。
(5)在單元柱內(nèi)實現(xiàn)Interface相交裁剪,在單元柱內(nèi)實現(xiàn)兩個或者多個Interface四邊形求交裁剪算法(見3.3節(jié))。將隸屬同一斷層的Interface處理后的有效曲面拼接起來,進行Delaunay三角剖分生成相交裁剪后的三角網(wǎng)化斷層面。
(6)后處理,利用斷層真實邊界條件對相交裁剪后的三角網(wǎng)化斷層面進行邊界裁剪,生成最終的斷層面。后處理會消除斷層網(wǎng)格化造成的 “假相交”切割關系(見3.1節(jié))和斷層延伸標記邊界(見3.2節(jié)),保證斷層網(wǎng)絡模型空間拓撲的合理性。
下文將重點對斷層間空間關系構建、斷層邊界延伸預處理和Interface相交裁剪3個核心步驟進行詳細介紹。
根據(jù)斷層間空間關系rij的定義,使用空間關系二叉樹和相交關系表記錄斷層空間位置上下和是否相交的信息。空間關系二叉樹的數(shù)據(jù)結構定義如下:
FaultRelationBinaryTree {
樹結點集合D:每個結點一一關聯(lián)斷層集合F={Fi}中對應斷層;
結點關系集合H:結點關系為二元關系:
(1)D中存在唯一的根結點root,在關系H中無前驅;除根結點,其他結點必有唯一前驅結點,即父結點。
(2)除葉結點,每個結點至多有兩個孩子結點,即左子結點和右子結點。
根據(jù)斷層空間關系集合R={rij}構建結點關系,規(guī)定:父結點斷層與左子結點斷層和右子結點斷層都相交;左子結點斷層位于父結點斷層上方;右子結點斷層位于父結點斷層下方。從斷層主輔關系,相對于子結點,父結點為主斷層結點。
空間關系二叉樹描述了有切割聯(lián)系的斷層關系,當斷層數(shù)量過多需要降低樹深度,空間關系二叉樹應滿足平衡二叉樹的結構要求(圖3)。在整個建模工區(qū),相互有切割聯(lián)系的斷層關系會形成多個空間關系二叉樹,整體構成空間關系森林??臻g關系二叉樹按中序序列進行存儲,每個斷層記錄按中序遍歷的序列位置編號,形成斷層關系存儲表(圖3c)。根據(jù)存儲序列編號,無需遍歷二叉樹,通過編號順序可判斷出兩個斷層的空間上下關系,如判斷斷層F2與F4的空間上下關系,F2的中序遍歷位置索引號(5)大于F4的索引號(1),故F2位于F4的下方。
圖3 斷層模型剖面圖(a)、空間關系二叉樹(b)和存儲表(c)Fig.3 Fault model section(a)、 Space relationship binary tree(b) and Storage table(c)
相交關系表記錄任意兩個斷層間是否存在切割關系。相交關系表的意義在于處理斷層網(wǎng)格化后的“假相交”關系,保證斷層建模的正確性。斷層“假相交”切割由規(guī)則網(wǎng)格分辨率引起。如圖4(a)中斷層F3切割F2,F2切割F1,而F1和F3實際中并沒有相交。但是按斷層網(wǎng)格化處理,由于網(wǎng)格分辨率,F3網(wǎng)格化后與F1相交(圖4b)?!凹傧嘟弧鼻懈铌P系需要在相交關系表中設置F1和F3不相交,Interface相交裁剪步驟不處理F1和F3的切割關系,后處理步驟根據(jù)斷層真實邊界進行裁剪。
圖4 斷層切割示意圖:(a) 斷層實際分布、(b) 圖框網(wǎng)格內(nèi)斷層網(wǎng)格化后結果和相交關系表(c)Fig.4 Faults cutting: (a) actual faults distribution, (b) faults in regular grid and intersecting relationship table (c)
由于斷層面裁剪在單個單元柱內(nèi)進行處理,而裁剪前的斷層所覆蓋的單元柱會超出斷層實際區(qū)域。如圖5a所示,斷層F1與F2相交,斷層所覆蓋單元柱有重疊域,也有非重疊域。若F1切割F2且F2在F1下方,在非重疊域上F2在裁剪處理后會出現(xiàn)錯誤(圖5b),其原因在于F1為主斷層切割輔斷層F2,但是F2的單元柱超出了F1的覆蓋范圍。若將斷層F1延伸到F1和F2邊界的最小公共外包(圖5c),可以保證斷層裁剪的正確處理,得到如圖5d所示結果。因此,斷層邊界延伸預處理就是對存在相交關系的主輔斷層,將切割主斷層的區(qū)域范圍外延至所有被切割輔斷層區(qū)域范圍的最小公共外包,對主斷層延伸部分進行標記。斷層面后處理階段應舍棄標記部分。
圖5 斷層延伸處理:(a)待處理斷層;(b)不合理情況;(c)延伸處理;(d)最終結果Fig.5 Fault extension preprocess:(a) original faults; (b) unreasonable case; (c) fault extension; (d) final result
斷層相交裁剪只需在單元柱內(nèi)實現(xiàn)Interface相交裁剪。在單元柱內(nèi),斷層面Interface的4個頂點都在單元柱的棱上(圖1),相鄰頂點構成的邊位于側柱面上,稱為面邊。位于Interface內(nèi)且與單元柱內(nèi)部有交集的邊或線稱為內(nèi)邊,內(nèi)邊的頂點可能在面邊上,也可能位于Interface內(nèi)。face是Interface經(jīng)過相交裁剪后保留的有效部分。最終Interface的有效部分由若干face構成。
根據(jù)交點個數(shù)和位置差異,單元柱內(nèi)兩個不全等Interface的相交可以枚舉出10種情況(圖6)。
圖6 單元柱內(nèi)Interface間相交示意圖Fig.6 Intersection between two Interfaces in the cell column
從圖6可以得到兩個Interface相交的特點:①交點一定在單元柱的側面或者棱上;②交點是兩個Interface在單元柱同一側面上的面邊交點;③交線是交點順次連接而成;④兩個順次連接的交點形成的邊可以是內(nèi)邊或者面邊。根據(jù)以上特點,單元柱內(nèi)Interface求交裁剪步驟如下:
(1)根據(jù)單元柱內(nèi)斷層的空間關系二叉樹中序遍歷,新遍歷的Interface(記為newInterface)和所有已經(jīng)完成遍歷處理的Interface(記為oldInterface)進行面邊求交(圖7a,b),newInterface每個面邊被oldInterface遞歸打斷成多段;
圖7 單元柱內(nèi)Interface相交裁剪(F1切割F2和F3,F2在F1下方,F3在F1上方): (a) F2為newInterface, F1為oldInterface; (b) F3為newInterface,L1為F1和F2的有效內(nèi)邊interiorEdge,L2為F2的有效面邊段;(c) L3為F3和F1相交的內(nèi)邊,虛線為無效的內(nèi)邊和面邊段,P為L1和F3交點Fig.7 Interface cutting in a cell column (F1 cut F2 and F3, F2 below F1, F3 above F1): (a) F2 is the newinterface, F1 is oldinterface; (b) F3 is newinterface, L1 is the valid interior edge of F1 and F2, L2 is the valid face edge segment of F2; (c) L3 is the interior edge of the intersection of F3 and F1, the dashed line is the invalid interior and face edge segments, and P is the intersection of L1 and F3
(2)判斷newInterface面邊的每個分段與oldInterface的空間上下關系是否滿足newInterface和oldInterface的空間上下關系,不滿足的面邊段被剔除,滿足的段為有效面邊段,如圖7bF2的面邊處理;(3)將newInterface與oldInterface的面邊交點依次連接形成內(nèi)邊,并剔除與面邊重合的內(nèi)邊,如圖7a內(nèi)邊L1。內(nèi)邊實際為newInterface和oldInterface的交線;
(4)將oldInterface上有效內(nèi)邊與新的內(nèi)邊在XY投影面上計算交點,記為P(圖7c),內(nèi)邊相應被打斷為多段;
(5)判斷內(nèi)邊的每個分段與oldInterface的空間上下關系是否滿足newInterface和oldInterface的空間上下關系,不滿足的內(nèi)邊段被剔除,滿足的段為有效內(nèi)邊;
(6)使用深度優(yōu)先遍歷將所有有效的面邊段和內(nèi)邊段構成封閉的最小環(huán),即為newInterface相交處理后的有效曲面(圖7d)。
單元柱內(nèi)所有Interface重復以上步驟進行。處理完成后,將隸屬同一斷層的Interface處理后的有效曲面拼接起來,進行Delaunay三角剖分生成相交裁剪后的三角網(wǎng)化斷層面,進入后處理流程。
4 實驗分析
利用Visual C++2019實現(xiàn)了本文方法,測試環(huán)境為CPU R7 5800H,16G內(nèi)存和64位微軟Windows 10操作系統(tǒng)。首先利用花狀構造斷層數(shù)據(jù)進行測試,并與SKUA—GOCAD 19進行了對比;然后利用煤礦三維地震斷層解釋數(shù)據(jù)進行驗證。
花狀構造的斷層網(wǎng)絡測試
數(shù)據(jù)中有8個斷層(圖8a),其中,F1為主斷層,4個分支斷層F2、F3、F4、F5與F1相交,F2切割斷層F3、F8,F5切割其他與之相交斷層?;W(wǎng)格分辨率為1002和20002,利用二維薄板樣條插值算法逐一構建網(wǎng)格化斷層曲面,并根據(jù)斷層建??傮w流程進行自動化處理。得到建模結果如圖8所示,斷層網(wǎng)絡切割處理結果符合斷層組合關系,基網(wǎng)格分辨率越高,曲面光滑性越好。
圖8 花狀構造斷層建模:(a) 源數(shù)據(jù):斷層棱邊;(b) 單個斷層面生成結果;(c) 基網(wǎng)格分辨率1002的結果;(d) 基網(wǎng)格分辨率20002的結果Fig.8 Flower structure fault modeling: (a) source data: fault edges; (b) fault surface; (c) results of base grid resolution 1002 and (d) results of base grid resolution 20002
以此花狀構造數(shù)據(jù)為測試數(shù)據(jù),對比本文方法與SKUA—GOCAD 2019在建模流程和處理結果上的差異。在SKUA—GOCAD 19中利用Structural Modeling功能,輸入圖8a中8個斷層的源數(shù)據(jù),采用DSI空間插值算法生成斷層面,設置斷層間的主輔關系,得到建模結果如圖9a所示,結果中斷層F1和F2以及F2和F3的斷層關系未能正確處理。利用GOCAD的Edit Contacts交互式對建模結果進行修正,最終得到斷層切割關系正確的模型效果如圖9b所示。該測試表明,相比SKUA—GOCAD斷層建模方法,本文方法具有較好的算法穩(wěn)定性,能夠減少交互過程,提高斷層建模的自動化程度。
圖9 SKUA—GOCAD2019花狀構造斷層建模結果: (a)直接建模結果; (b)交互修正斷層邊界Fig.9 Flower structure fault modeling result in SKUA—GOCAD2019: (a) direct modeling results; (b) interactive corrected fault boundaries
仍以花狀構造數(shù)據(jù)為測試數(shù)據(jù)分析基網(wǎng)格分辨率對建模性能的影響。在相同測試環(huán)境下記錄不同分辨率下的斷層建模所需時間,結果如表1所示。從測試結果可以看出,對于地質模型基網(wǎng)格分辨率5002,本文方法在180 s內(nèi)運行,可以滿足一般建模性能需求,但對于高分辨基網(wǎng)格,建模耗時在分鐘和小時數(shù)量級,還需進一步提升建模效率。
表1 不同基網(wǎng)格分辨率下的斷層建模運行時間Table 1 Running time of fault cross-cropping at different base grid resolution
根據(jù)某煤礦三維地震斷層解釋數(shù)據(jù),選取了具有復雜切割關系的9個斷層,斷層輸入數(shù)據(jù)如圖10a所示,其中,斷層DF1切割斷層DF0、DF2,斷層DF5切割斷層DF4、DF6、DF7、DF8?;W(wǎng)格分辨率為5002,利用二維薄板樣條插值算法逐一構建網(wǎng)格化斷層曲面,并根據(jù)斷層建模總體流程進行自動化處理。建模結果如圖10c所示,斷層網(wǎng)絡切割處理結果符合斷層關系,如圖10d中斷層DF5切割斷層DF4、DF6、DF7、DF8的處理結果正確。
圖10 某煤礦斷層網(wǎng)絡建模:(a) 源數(shù)據(jù):斷層棱邊;(b) 單個斷層面生成結果;(c) 斷層建模結果;(d) 局部放大圖Fig.10 Fault network modeling in a coalmine : (a) source data: fault edges; (b) fault surface; (c) fault modeling results; (d) local magnification plot
針對復雜斷層面切割裁剪處理算法穩(wěn)定性以及斷層建模流程自動化等問題,筆者等提出了基于規(guī)則網(wǎng)格的復雜斷層網(wǎng)絡處理與自動化建模的方法和流程,詳細討論了基于網(wǎng)格化的斷層網(wǎng)絡模型形式化理論表達、建模流程中的斷層網(wǎng)絡空間關系構建以及相交裁剪處理算法等核心步驟。通過基于網(wǎng)格化的斷層網(wǎng)絡模型形式化表達和Interface定義,將斷層面建模轉化為在基網(wǎng)格單元柱內(nèi)根據(jù)斷層空間關系進行Interface的求交和裁剪處理,該過程具有較好的算法穩(wěn)定性。利用斷層的削截主輔關系構建斷層間空間關系,通過空間關系二叉樹和相交關系表管理所有斷層的空間關系,按照統(tǒng)一的流程進行斷層切割裁剪和后處理,最終得到符合斷層關系的建模結果。
測試驗證表明,本文方法可以有效處理多條互相切割、主輔關系復雜的斷層網(wǎng)絡,具有較好的算法穩(wěn)定性;通過與SKUA—GOCAD斷層建模方法對比,能夠減少交互過程,提高斷層建模的自動化程度。
本文方法缺點和不足為:① 不能處理多值斷層曲面,如S型斷層曲面和和直立或近直立斷層;② 斷層空間關系需要嚴格的空間上下關系, 故不能處理X型斷層組合;③ 基于基網(wǎng)格的斷層表達是一種擬合方法,斷層面的精度依賴基網(wǎng)格分辨率;④ 由于采用規(guī)則網(wǎng)格,在高分辨率基網(wǎng)格上,建模性能低;⑤ 最終構建的斷層不規(guī)則三角形網(wǎng)格質量還比較差,沒有進行網(wǎng)格優(yōu)化處理。
本文方法的后續(xù)研究可以采用自適應基網(wǎng)格和并行處理提升建模性能,也可以將基網(wǎng)格由規(guī)則網(wǎng)格推廣到三角形網(wǎng)格,提高基于GTP和主TIN等地質建模方法的自動化水平。