基于層面結(jié)構(gòu)的地質(zhì)塊體拓?fù)潢P(guān)系的自動(dòng)構(gòu)建

安聰榮，劉 展，王心眾

中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，山東青島266555

1 引 言

以地震解釋得到的地質(zhì)界面（如地層界面、斷層面等）數(shù)據(jù)為基礎(chǔ)，采用一定的曲面重建算法，在三維空間中可得到拓?fù)湟恢碌膶用婺Ｐ停?－2］，這是地質(zhì)構(gòu)造建模的基本成果，但此時(shí)的模型是沒有描述塊體結(jié)構(gòu)的地質(zhì)模型，無法在其基礎(chǔ)上進(jìn)行實(shí)體（屬性）建模［3－4］，這在一定程度成為導(dǎo)致目前油氣藏構(gòu)造模型和屬性模型之間難以融合的重要原因之一。

因此，鑒于地質(zhì)塊體結(jié)構(gòu)模型的提取在地質(zhì)體建模中的重要地位，許多研究者對此進(jìn)行了深入地研究。文獻(xiàn)［5］于1988年提出一種用塊體生成語言（block generation language，BGL）來完成地質(zhì)塊體識別的方法，這種方法對簡單體的識別非常適合，而對于地質(zhì)建模中經(jīng)常遇到的復(fù)雜地質(zhì)體的識別能力有限。之后，文獻(xiàn)［6］通過引入拓?fù)涞母拍顚?shí)現(xiàn)地質(zhì)塊體識別。文獻(xiàn)［7］提出基于層面結(jié)構(gòu)的地質(zhì)塊體構(gòu)建方法，該方法采用一種改進(jìn)的半邊數(shù)據(jù)結(jié)構(gòu)，然后依據(jù)三角形間的拓?fù)潢P(guān)系構(gòu)建地質(zhì)塊體，由于需要遍歷層面模型中的所有三角形，因此在地層三角形數(shù)量極大時(shí)，方法的效率不高。文獻(xiàn)［4］是基于簡化面片的閉合塊體構(gòu)造方法，由于在簡化后的層面模型包含全局的拓?fù)湫畔ⅲ虼嗽摲椒o需遍歷層面模型中所有的三角形，就能實(shí)現(xiàn)地質(zhì)塊體的識別與提取，極大提高算法的運(yùn)行效率，但該方法并沒有直接建立地質(zhì)塊體的拓?fù)湫畔?。文獻(xiàn)［8］利用上下地層中三角網(wǎng)的重疊次數(shù)分離不同類型的子面，然后依此建立子面與地質(zhì)塊體之間的對應(yīng)關(guān)系，這種方法同樣沒有直接建立和保存地質(zhì)塊體的拓?fù)湫畔ⅰＳ纱丝梢?，目前地質(zhì)塊體結(jié)構(gòu)的構(gòu)建方法在算法效率提高方面取得較大進(jìn)展，然而卻對地質(zhì)塊體與面片之間以及地質(zhì)塊體與地質(zhì)塊體之間拓?fù)潢P(guān)系自動(dòng)建立的研究存在明顯的不足，從而在一定程度上限定三維地質(zhì)體拓?fù)洳樵兒头治觥?/p>

為此，本文提出一種地質(zhì)塊體拓?fù)潢P(guān)系自動(dòng)構(gòu)建技術(shù)，其核心是建立地質(zhì)塊體與地質(zhì)面片之間組合關(guān)系，并間接表達(dá)地質(zhì)塊體之間的鄰接關(guān)系。為了提高算法的效率，引入方向邊和方向三角形的概念以簡化層面模型，使得算法在運(yùn)算中無需遍歷層面模型的所有三角形，就能實(shí)現(xiàn)地質(zhì)塊體拓?fù)潢P(guān)系的自動(dòng)構(gòu)建。

2 面片的概念及其構(gòu)建關(guān)鍵算法

地層面上廣泛分布的斷層線、地層相交線與曲面的邊界（外邊界、內(nèi)邊界）共同作用，可將地層面劃分成若干個(gè)相鄰的區(qū)域，這些區(qū)域即為面片。在當(dāng)前的許多地質(zhì)數(shù)據(jù)模型中都將面片作為模型構(gòu)建元素之一，如3DFDS（formal data structure）模型［9－10］、矢量與柵格集成的三維空間數(shù)據(jù)模型［11－13］等。在這些模型中，面片一般都為曲面片，由邊界和以不規(guī)則三角網(wǎng)（irregular triangle net，TIN）表示的內(nèi)部區(qū)域兩部分構(gòu)成，其內(nèi)的每個(gè)三角形都只能被兩個(gè)不同的地層閉合體所共享。在實(shí)際的地質(zhì)建模中所涉及的面片主要有斷層面片、地層面片和邊界面片，其邊界也可能由斷層線、地層相交線、研究區(qū)邊界線等閉合而構(gòu)成。面片可以將復(fù)雜體分割為左右兩個(gè)部分，即正體（位于面片正法向量方向一側(cè)的閉合體）和負(fù)體（位于面片負(fù)法向量方向一側(cè)的閉合體）。

面片的構(gòu)建主要分為兩個(gè)步驟進(jìn)行：

（1）面片邊界的確定。正確建立面片邊界不僅是確定地層面邊界位置的需要，也是正確構(gòu)建面片之間拓?fù)潢P(guān)系的需要。在實(shí)際建模時(shí)，如何從大量弧段中（主要為斷層線、地層相交線、邊界線等）中提取面片閉合邊界，是構(gòu)建地質(zhì)面片面臨的首要問題。通常采用的方法有兩種：一是人工構(gòu)建，二是自動(dòng)構(gòu)建。前者顯然已不能滿足復(fù)雜地質(zhì)建模的時(shí)代要求，后者的難點(diǎn)在于邊界弧段與地質(zhì)面片拓?fù)潢P(guān)系的自動(dòng)建立，可以采用2D GIS中多邊形拓?fù)潢P(guān)系自動(dòng)構(gòu)建技術(shù)來實(shí)現(xiàn)地層面片邊界弧段的自動(dòng)成區(qū)［14］。

（2）面片內(nèi)部區(qū)域的構(gòu)建。不同類型的面片所采用的技術(shù)也是不同的，如地層面片、斷層面片構(gòu)建常用的方法是將上述步驟（1）所提取的面片邊界線作為約束線，利用限定Delaunay算法實(shí)現(xiàn)它的構(gòu)建［15］，這樣可以確保邊界線在其中的存在性，而地層邊界面片則可以采用文獻(xiàn)［16］提出的相鄰輪廓線同步前進(jìn)法來實(shí)現(xiàn)。

3 地質(zhì)塊體拓?fù)潢P(guān)系的自動(dòng)構(gòu)建算法

3.1 算法的基本思路

導(dǎo)入地質(zhì)層面模型，提取所有地質(zhì)面片，并對其進(jìn)行預(yù)處理，然后依據(jù)3DFDS模型、矢量與柵格集成模型等的原理，對面片設(shè)置正體和負(fù)體字段，并從任意面片開始，按照一定的遍歷規(guī)則建立所有面片的正體和負(fù)體，從而實(shí)現(xiàn)地質(zhì)塊體拓?fù)潢P(guān)系的自動(dòng)構(gòu)建。

3.2 基本概念

為了描述算法原理的方便，引入以下概念。懸掛面片：正體和負(fù)體是同一個(gè)體的面片，類似二維空間中的懸掛弧段。

面片排序表：算法運(yùn)行中的一個(gè)變量，記錄邊界弧段ID和共享該邊界弧段并且按照右手法則排序的地質(zhì)面片集。

活動(dòng)面片集：是算法運(yùn)行中的一個(gè)變量，命名為PS，記錄激活面片的ID和標(biāo)記，在算法運(yùn)行起初是空集。

邊界弧段的重度：指共享邊界弧段的面片個(gè)數(shù)，是算法運(yùn)行中的一個(gè)變量，利用邊界弧段的重度可以判斷一個(gè)面片是否為懸掛面片。

方向邊：邊界弧段上的一條矢量線段，僅有一個(gè)起點(diǎn)和一個(gè)終點(diǎn)，沒有中間點(diǎn)。

方向三角形：在層面數(shù)據(jù)的三角形集中共享方向邊的三角形。方向三角形的頂點(diǎn)中有兩個(gè)是方向邊的端點(diǎn)，而第3個(gè)是來自各個(gè)不同的面片，因此，方向三角形隱含著層面數(shù)據(jù)中面片的拓?fù)湫畔ⅰ?/p>

3.3 數(shù)據(jù)預(yù)處理

3.3.1 懸掛面片的消除

懸掛面片往往處于一個(gè)最小閉合體內(nèi)部，但又不是這個(gè)閉合體的內(nèi)邊界，它是懸掛于這個(gè)閉合體內(nèi)。它的存在會(huì)破壞地質(zhì)塊體拓?fù)潢P(guān)系的自動(dòng)建立，對于懸掛面片的處理有兩種策略，一是預(yù)先除去模型內(nèi)的所有懸掛面片，另一是在閉合體構(gòu)造過程中除去懸掛面片。本文選擇第1種策略，遍歷所有面片，求出每條邊界弧段的重度，凡是含有重度為1的邊界弧段的面片都認(rèn)為是懸掛面片，給予刪除。

3.3.2 共享邊界弧段面片的排序

面片排序表用于記錄邊界弧段與共享每一個(gè)邊界弧段面片，并且按照右手法則對面片進(jìn)行排序。但是由于邊界弧段是復(fù)雜的折線，直接利用邊界弧段和共享邊界弧段的面片來判斷面片的順序關(guān)系變得異常復(fù)雜。為了簡化算法設(shè)計(jì)的難度，采用方向邊和方向三角形的概念。

如圖1所示，假設(shè)V1V2是邊界弧段L中的一線段（起點(diǎn)為V1，終點(diǎn)為V2），三角形a、b、c分別來自各面片且共享線段V1V2，則稱線段V1V2為方向邊，三角形a、b、c為方向三角形，利用右手法則可以得到方向三角形的排序，從而得到共享邊界線L的所有面片排序。

圖1 方向邊和方向三角形Fig.1 Direction－edge and direction－triangle

按照上述思路，在算法運(yùn)行前，對每條邊界弧段，提取方向邊和方向三角形，對共享該邊的面片進(jìn)行排序，排序后的結(jié)果保存在面片排序表中，值得一提的是，該排序是循環(huán)排序。

3.3.3 面片邊界弧段方向一致性處理

由于面片生成時(shí)，邊界是由若干個(gè)邊界弧段相連構(gòu)成的閉合線，各個(gè)邊界弧段的方向可能出現(xiàn)不一致情況，導(dǎo)致無法定義一個(gè)面片的正體和負(fù)體屬性。因此在算法運(yùn)行前，需要通過遍歷面片邊界組成弧段，調(diào)整首尾節(jié)點(diǎn)以及中間點(diǎn)的ID順序，并保存處理結(jié)果。

3.4 算法詳細(xì)步驟

算法輸入：輸入層面模型的邊界弧段和經(jīng)預(yù)處理后的面片集。

算法輸出：輸出每個(gè)面片的正體和負(fù)體。

具體步驟如下：

（1）得到第1面片A，并設(shè)置為當(dāng)前面片，轉(zhuǎn)到步驟（2）。

（2）判斷面片A的正體和負(fù)體是否為空，如果都非空，轉(zhuǎn)到步驟（1）。當(dāng)所有面片處理完畢后，算法結(jié)束。若正體為空，轉(zhuǎn)到步驟（3）。若負(fù)體為空，轉(zhuǎn)到步驟（4）。

（3）創(chuàng)建一個(gè)新的正體B，轉(zhuǎn)到步驟（5）。

（4）創(chuàng)建一個(gè)新的負(fù)體B，轉(zhuǎn)到步驟（5）。

（5）將面片A加入到PS中，并標(biāo)記面片A，轉(zhuǎn)到步驟（6）。

（6）獲取面片A的一條未標(biāo)記的邊界弧段，命名為L，并作標(biāo)記，轉(zhuǎn)到步驟（7）。

（7）判斷B是面片A的正體還是負(fù)體，如果是正體轉(zhuǎn)到步驟（8），如果是負(fù)體轉(zhuǎn)到步驟（10）。

（8）檢查L的方向是否和方向邊的方向一致，如果一致，選擇面片A之后的下一個(gè)面片A′作為擴(kuò)展面片；如果不一致，選擇面片A之前的面片A′作為擴(kuò)展面片，判斷面片A′在PS中是否存在。如果不存在，則加入到其中，并轉(zhuǎn)到步驟（9）；如果存在，轉(zhuǎn)到步驟（12）。

（9）判斷邊界弧段L的起點(diǎn)和終點(diǎn)的前后順序是否和面片A′的一致，如果一致，將B設(shè)置為面片A′的負(fù)體，如果否，設(shè)置為面片A′的正體，轉(zhuǎn)到步驟（12）。

（10）檢查L的方向是否和方向邊的方向一致，如果一致，選擇面片A之前的面片A′作為擴(kuò)展面片；如果不一致，選擇面片A之后下一個(gè)面片A′作為擴(kuò)展面片，然后判斷面片A′在PS中是否存在。如果否，則加入到其中，并轉(zhuǎn)到步驟（11）；如果存在，轉(zhuǎn)到步驟（12）。

（11）判斷邊界弧段L的起點(diǎn)和終點(diǎn)的前后順序是否和面片A′一致，如果一致，將B設(shè)置為面片A′的正體；如果否，設(shè)置為面片A′的負(fù)體，轉(zhuǎn)到步驟（12）。

（12）判斷面片A的邊界弧段是否都已標(biāo)記過，如果是，轉(zhuǎn)到步驟（13）；如果否，轉(zhuǎn)到步驟（6）。

（13）判斷PS中所有面片是否標(biāo)記過，如果是，清空PS中的所有記錄并轉(zhuǎn)到步驟（1）；如果否，選擇一個(gè)沒有標(biāo)記的面片作為當(dāng)前面片A，轉(zhuǎn)到步驟（5）。

3.5 算法涉及的部分?jǐn)?shù)據(jù)結(jié)構(gòu)

struct thirdverd ∥方向三角形

4 試 驗(yàn)

4.1 數(shù)據(jù)準(zhǔn)備

試驗(yàn)采用的數(shù)據(jù)是基于鉆井?dāng)?shù)據(jù)構(gòu)建的層面數(shù)據(jù)和利用相鄰輪廓線同步前進(jìn)法得到的地層邊界面數(shù)據(jù)，這些數(shù)據(jù)以三角網(wǎng)格的形式存放。

4.2 算法的實(shí)現(xiàn)與測試

利用C＃語言，筆者實(shí)現(xiàn)了本文所提出的方法，并對其進(jìn)行測試。在試驗(yàn)中，使用某地質(zhì)區(qū)域的層面數(shù)據(jù)（共8層，包括1個(gè)地表面和7個(gè)地下層面數(shù)據(jù)）以及邊界面數(shù)據(jù)，共有36個(gè)地質(zhì)面片（如圖2所示）。從表1可以看出，三角形數(shù)量達(dá)19 554的層面模型數(shù)據(jù)，在引入方向三角形后，在內(nèi)存中僅保存142個(gè)三角形（即方向三角形，包含有全局拓?fù)湫畔ⅲ?，三角形的?shù)量大幅度減少，極大提高了算法的運(yùn)行效率（本算例的運(yùn)行時(shí)間為256ms）。圖2是算法運(yùn)行前的輸入數(shù)據(jù)，圖3是根據(jù)算法運(yùn)行后所得的地質(zhì)塊體與面片拓?fù)湫畔⒍⒌牡刭|(zhì)塊體結(jié)構(gòu)模型。結(jié)果表明，該算法運(yùn)行高效、準(zhǔn)確。

表1 試驗(yàn)數(shù)據(jù)預(yù)處理結(jié)果Tab.1 Results of pre－processing experimental data

圖2 層面模型數(shù)據(jù)Fig.2 Experimental data of stratified geological structure

圖3 算法運(yùn)行后的三維塊體模型Fig.3 3DGeological blocks model after running algorithm

5 結(jié) 論

通過以上研究以及試驗(yàn)證明，算法設(shè)計(jì)合理，結(jié)果可靠、準(zhǔn)確，并得出以下幾點(diǎn)結(jié)論：

（1）本文提出的地質(zhì)塊體模型構(gòu)建方法的核心是構(gòu)建地質(zhì)面片的正體和負(fù)體屬性，該方法不僅建立地質(zhì)塊體與其組成面片的拓?fù)潢P(guān)系，同時(shí)也間接建立地質(zhì)塊體與地質(zhì)塊體之間的拓?fù)湫畔?，有效彌補(bǔ)當(dāng)前地質(zhì)塊體模型構(gòu)建過程中未能就地質(zhì)體拓?fù)潢P(guān)系信息進(jìn)行建立的不足。

（2）該算法同時(shí)也兼顧了效率問題，引入方向邊和方向三角形的概念，簡化地質(zhì)層面數(shù)據(jù)模型，使得算法無需遍歷所有三角形，就能完成地質(zhì)面片正體與負(fù)體屬性的計(jì)算。

（3）當(dāng)前的許多地質(zhì)數(shù)據(jù)模型中都將面片作為模型構(gòu)建元素之一，在建模過程都面臨著地質(zhì)體與面片拓?fù)潢P(guān)系的自動(dòng)構(gòu)建，因此，該算法在這些模型的構(gòu)建過程具有重要意義。

［1］ XIONG Zuqiang，HE Huaijian，XIA Yanhua.Study on Technology of 3DStratum Modeling and Visualization Based on TIN［J］.Rock and Soil Mechanics，2007，28（9）：1954－1958.（熊祖強(qiáng)，賀懷建，夏艷華.基于TIN的三維地層建模及可視化技術(shù)研究［J］.巖土力學(xué)，2007，28（9）：1954－1958.）

［2］ WANG Chunxiang，BAI Shiwei，HE Huaijian.Study on Geological Modeling in 3DStratum Visualization［J］.Chinese Journal of Rock Mechanics and Engineering，2003，22（10）：1722－1726.（王純祥，白世偉，賀懷建.三維地層可視化中地質(zhì)建模研究［J］.巖石力學(xué)與工程學(xué)報(bào)，2003，22（10）：1722－1726.）

［3］ WEI Jia，TANG Jie，YUE Chengqi，et al.Study of 3D Geological Structure Model Building［J］.Geophysical Prospecting for Petroleum，2008，47（4）：319－327.（魏嘉，唐杰，岳承祺，等.三維地質(zhì)構(gòu)造建模技術(shù)研究［J］.石油物探，2008，47（4）：319－327.）

［4］ ZHU Lian，TANG Jie，YUAN Chunfeng.Construction of Coherent Structure from 3DGeological Model［J］.Journal of Image and Graphics，2009，14（12）：2582－2587.（朱煉，唐杰，袁春風(fēng).3維地質(zhì)體模型中閉合結(jié)構(gòu)的提?。跩］.中國圖象圖形學(xué)報(bào)，2009，14（12）：2582－2587.）

［5］ HELIOT D.Generating a Blocky Rock Mass［J］.International Journal of Rock Mechanics and Mining Science and Geomech Abstracts，1988，25（3）：127－138.

［6］ JING L，STEPHANSSON O.Topological Identification of Blocky Assemblages for Jointed Rock Masses［J］.International Journal of Rock Mechanics and Mining Science and Geomechanics Abstracts，1994，31（2）：163－172.

［7］ MENG Xianhai，YANG Qin，LI Jigang.Construction of Coherent 3DGeological Blocks from Stratified Geological Structure［J］.Journal of Beijing University of Aeronautics and Astronautic，2005，31（2）：182－186.（孟憲海，楊欽，李吉?jiǎng)?基于層面結(jié)構(gòu)的三維閉合地質(zhì)區(qū)塊構(gòu)造算法［J］.北京航空航天大學(xué)學(xué)報(bào)，2005，31（2）：182－186.）

［8］ LI Yuanheng，CHEN Guoliang，LIU Xiuguo，et al.The Topology－oriented Method of Building 3DGeological Block Model Based on Primary TIN［J］.Rock and Soil Mechanics，2010，31（6）：1902－1906.（李元亨，陳國良，劉修國，等.主TIN模式下面向拓?fù)涞娜S地質(zhì)塊體構(gòu)建方法［J］.巖土力學(xué)，2010，31（6）：1902－1906.）

［9］ MOLENAAR M.A Formal Data Structure for Three Dimensional Vector Maps［C］∥Proceedings of 4th International Symposium on Spatial Data Handling.Zurich：［s.n.］，1990：830－843.

［10］ MOLENAAR M.A Topology for 3DVector Maps［J］.ITC International Technology Communication Journal，1992，1：25－33.

［11］ LI Qingquan，LI Deren.Research on the Conceptual Frame of the Integration of 3DSpatial Datamodel［J］.Acta Geodaetica et Cartographica Sinica，1998，27（4）：325－330.（李清泉，李德仁.三維空間數(shù)據(jù)模型集成的概念框架研究［J］.測繪學(xué)報(bào)，1998，27（4）：325－330.）

［12］ CHENG Penggen，GONG Jianya.Design of Threedimensional Spatial Data Model and Its Data Structure in Geological Exploration Engineering［J］.Acta Geodaetica et Cartographica Sinica，2001，30（1）：74－81.（程朋根，龔健雅.地勘工程3維空間數(shù)據(jù)模型及其數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)［J］.測繪學(xué)報(bào)，2001，30（1）：74－81.）

［13］ GONG Jianya，XIA Zongguo.An Integrated Data Model in Three Dimensional GIS［J］.Journal of Wuhan Technical University of Surveying and Mapping，1997，22（1）：7－15.（龔健雅，夏宗國.矢量與柵格集成的三維數(shù)據(jù)模型［J］.武漢測繪科技大學(xué)學(xué)報(bào)，1997，22（1）：7－15.）

［14］ QI Hua.The Optimization and Improvement for the Algorithm Steps on the Automatic Creation of Topological Relation of Polygons［J］.Acta Geodaetica et Cartographica Sinica，1997，26（3）：254－260.（齊華.自動(dòng)建立多邊形拓?fù)潢P(guān)系算法步驟的優(yōu)化與改進(jìn)［J］.測繪學(xué)報(bào)，1997，26（3）：254－260.）

［15］ SLOAN S W.A Fast Algorithm for Constructing Delaunay Triangulations in the Plane［J］.Advances in Engineering Soft Ware，1987，9（1）：34－55.

［16］ GANAPATHY S，DENEHY T G.A New General Triangulation Method for Planar Contours［J］.Computer Graphics，1982，16（3）：69－75.