基于层面结构的地质块体拓扑关系的自动构建

安聪荣，刘 展，王心众

中国石油大学（华东）地球科学与技术学院，山东青岛266555

1 引 言

以地震解释得到的地质界面（如地层界面、断层面等）数据为基础，采用一定的曲面重建算法，在三维空间中可得到拓扑一致的层面模型［1－2］，这是地质构造建模的基本成果，但此时的模型是没有描述块体结构的地质模型，无法在其基础上进行实体（属性）建模［3－4］，这在一定程度成为导致目前油气藏构造模型和属性模型之间难以融合的重要原因之一。

因此，鉴于地质块体结构模型的提取在地质体建模中的重要地位，许多研究者对此进行了深入地研究。文献［5］于1988年提出一种用块体生成语言（block generation language，BGL）来完成地质块体识别的方法，这种方法对简单体的识别非常适合，而对于地质建模中经常遇到的复杂地质体的识别能力有限。之后，文献［6］通过引入拓扑的概念实现地质块体识别。文献［7］提出基于层面结构的地质块体构建方法，该方法采用一种改进的半边数据结构，然后依据三角形间的拓扑关系构建地质块体，由于需要遍历层面模型中的所有三角形，因此在地层三角形数量极大时，方法的效率不高。文献［4］是基于简化面片的闭合块体构造方法，由于在简化后的层面模型包含全局的拓扑信息，因此该方法无需遍历层面模型中所有的三角形，就能实现地质块体的识别与提取，极大提高算法的运行效率，但该方法并没有直接建立地质块体的拓扑信息。文献［8］利用上下地层中三角网的重叠次数分离不同类型的子面，然后依此建立子面与地质块体之间的对应关系，这种方法同样没有直接建立和保存地质块体的拓扑信息。由此可见，目前地质块体结构的构建方法在算法效率提高方面取得较大进展，然而却对地质块体与面片之间以及地质块体与地质块体之间拓扑关系自动建立的研究存在明显的不足，从而在一定程度上限定三维地质体拓扑查询和分析。

为此，本文提出一种地质块体拓扑关系自动构建技术，其核心是建立地质块体与地质面片之间组合关系，并间接表达地质块体之间的邻接关系。为了提高算法的效率，引入方向边和方向三角形的概念以简化层面模型，使得算法在运算中无需遍历层面模型的所有三角形，就能实现地质块体拓扑关系的自动构建。

2 面片的概念及其构建关键算法

地层面上广泛分布的断层线、地层相交线与曲面的边界（外边界、内边界）共同作用，可将地层面划分成若干个相邻的区域，这些区域即为面片。在当前的许多地质数据模型中都将面片作为模型构建元素之一，如3DFDS（formal data structure）模型［9－10］、矢量与栅格集成的三维空间数据模型［11－13］等。在这些模型中，面片一般都为曲面片，由边界和以不规则三角网（irregular triangle net，TIN）表示的内部区域两部分构成，其内的每个三角形都只能被两个不同的地层闭合体所共享。在实际的地质建模中所涉及的面片主要有断层面片、地层面片和边界面片，其边界也可能由断层线、地层相交线、研究区边界线等闭合而构成。面片可以将复杂体分割为左右两个部分，即正体（位于面片正法向量方向一侧的闭合体）和负体（位于面片负法向量方向一侧的闭合体）。

面片的构建主要分为两个步骤进行：

（1）面片边界的确定。正确建立面片边界不仅是确定地层面边界位置的需要，也是正确构建面片之间拓扑关系的需要。在实际建模时，如何从大量弧段中（主要为断层线、地层相交线、边界线等）中提取面片闭合边界，是构建地质面片面临的首要问题。通常采用的方法有两种：一是人工构建，二是自动构建。前者显然已不能满足复杂地质建模的时代要求，后者的难点在于边界弧段与地质面片拓扑关系的自动建立，可以采用2D GIS中多边形拓扑关系自动构建技术来实现地层面片边界弧段的自动成区［14］。

（2）面片内部区域的构建。不同类型的面片所采用的技术也是不同的，如地层面片、断层面片构建常用的方法是将上述步骤（1）所提取的面片边界线作为约束线，利用限定Delaunay算法实现它的构建［15］，这样可以确保边界线在其中的存在性，而地层边界面片则可以采用文献［16］提出的相邻轮廓线同步前进法来实现。

3 地质块体拓扑关系的自动构建算法

3.1 算法的基本思路

导入地质层面模型，提取所有地质面片，并对其进行预处理，然后依据3DFDS模型、矢量与栅格集成模型等的原理，对面片设置正体和负体字段，并从任意面片开始，按照一定的遍历规则建立所有面片的正体和负体，从而实现地质块体拓扑关系的自动构建。

3.2 基本概念

为了描述算法原理的方便，引入以下概念。悬挂面片：正体和负体是同一个体的面片，类似二维空间中的悬挂弧段。

面片排序表：算法运行中的一个变量，记录边界弧段ID和共享该边界弧段并且按照右手法则排序的地质面片集。

活动面片集：是算法运行中的一个变量，命名为PS，记录激活面片的ID和标记，在算法运行起初是空集。

边界弧段的重度：指共享边界弧段的面片个数，是算法运行中的一个变量，利用边界弧段的重度可以判断一个面片是否为悬挂面片。

方向边：边界弧段上的一条矢量线段，仅有一个起点和一个终点，没有中间点。

方向三角形：在层面数据的三角形集中共享方向边的三角形。方向三角形的顶点中有两个是方向边的端点，而第3个是来自各个不同的面片，因此，方向三角形隐含着层面数据中面片的拓扑信息。

3.3 数据预处理

3.3.1 悬挂面片的消除

悬挂面片往往处于一个最小闭合体内部，但又不是这个闭合体的内边界，它是悬挂于这个闭合体内。它的存在会破坏地质块体拓扑关系的自动建立，对于悬挂面片的处理有两种策略，一是预先除去模型内的所有悬挂面片，另一是在闭合体构造过程中除去悬挂面片。本文选择第1种策略，遍历所有面片，求出每条边界弧段的重度，凡是含有重度为1的边界弧段的面片都认为是悬挂面片，给予删除。

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

面片排序表用于记录边界弧段与共享每一个边界弧段面片，并且按照右手法则对面片进行排序。但是由于边界弧段是复杂的折线，直接利用边界弧段和共享边界弧段的面片来判断面片的顺序关系变得异常复杂。为了简化算法设计的难度，采用方向边和方向三角形的概念。

如图1所示，假设V1V2是边界弧段L中的一线段（起点为V1，终点为V2），三角形a、b、c分别来自各面片且共享线段V1V2，则称线段V1V2为方向边，三角形a、b、c为方向三角形，利用右手法则可以得到方向三角形的排序，从而得到共享边界线L的所有面片排序。

图1 方向边和方向三角形Fig.1 Direction－edge and direction－triangle

按照上述思路，在算法运行前，对每条边界弧段，提取方向边和方向三角形，对共享该边的面片进行排序，排序后的结果保存在面片排序表中，值得一提的是，该排序是循环排序。

3.3.3 面片边界弧段方向一致性处理

由于面片生成时，边界是由若干个边界弧段相连构成的闭合线，各个边界弧段的方向可能出现不一致情况，导致无法定义一个面片的正体和负体属性。因此在算法运行前，需要通过遍历面片边界组成弧段，调整首尾节点以及中间点的ID顺序，并保存处理结果。

3.4 算法详细步骤

算法输入：输入层面模型的边界弧段和经预处理后的面片集。

算法输出：输出每个面片的正体和负体。

具体步骤如下：

（1）得到第1面片A，并设置为当前面片，转到步骤（2）。

（2）判断面片A的正体和负体是否为空，如果都非空，转到步骤（1）。当所有面片处理完毕后，算法结束。若正体为空，转到步骤（3）。若负体为空，转到步骤（4）。

（3）创建一个新的正体B，转到步骤（5）。

（4）创建一个新的负体B，转到步骤（5）。

（5）将面片A加入到PS中，并标记面片A，转到步骤（6）。

（6）获取面片A的一条未标记的边界弧段，命名为L，并作标记，转到步骤（7）。

（7）判断B是面片A的正体还是负体，如果是正体转到步骤（8），如果是负体转到步骤（10）。

（8）检查L的方向是否和方向边的方向一致，如果一致，选择面片A之后的下一个面片A′作为扩展面片；如果不一致，选择面片A之前的面片A′作为扩展面片，判断面片A′在PS中是否存在。如果不存在，则加入到其中，并转到步骤（9）；如果存在，转到步骤（12）。

（9）判断边界弧段L的起点和终点的前后顺序是否和面片A′的一致，如果一致，将B设置为面片A′的负体，如果否，设置为面片A′的正体，转到步骤（12）。

（10）检查L的方向是否和方向边的方向一致，如果一致，选择面片A之前的面片A′作为扩展面片；如果不一致，选择面片A之后下一个面片A′作为扩展面片，然后判断面片A′在PS中是否存在。如果否，则加入到其中，并转到步骤（11）；如果存在，转到步骤（12）。

（11）判断边界弧段L的起点和终点的前后顺序是否和面片A′一致，如果一致，将B设置为面片A′的正体；如果否，设置为面片A′的负体，转到步骤（12）。

（12）判断面片A的边界弧段是否都已标记过，如果是，转到步骤（13）；如果否，转到步骤（6）。

（13）判断PS中所有面片是否标记过，如果是，清空PS中的所有记录并转到步骤（1）；如果否，选择一个没有标记的面片作为当前面片A，转到步骤（5）。

3.5 算法涉及的部分数据结构

struct thirdverd ∥方向三角形

4 试 验

4.1 数据准备

试验采用的数据是基于钻井数据构建的层面数据和利用相邻轮廓线同步前进法得到的地层边界面数据，这些数据以三角网格的形式存放。

4.2 算法的实现与测试

利用C＃语言，笔者实现了本文所提出的方法，并对其进行测试。在试验中，使用某地质区域的层面数据（共8层，包括1个地表面和7个地下层面数据）以及边界面数据，共有36个地质面片（如图2所示）。从表1可以看出，三角形数量达19 554的层面模型数据，在引入方向三角形后，在内存中仅保存142个三角形（即方向三角形，包含有全局拓扑信息），三角形的数量大幅度减少，极大提高了算法的运行效率（本算例的运行时间为256ms）。图2是算法运行前的输入数据，图3是根据算法运行后所得的地质块体与面片拓扑信息而建立的地质块体结构模型。结果表明，该算法运行高效、准确。

表1 试验数据预处理结果Tab.1 Results of pre－processing experimental data

图2 层面模型数据Fig.2 Experimental data of stratified geological structure

图3 算法运行后的三维块体模型Fig.3 3DGeological blocks model after running algorithm

5 结 论

通过以上研究以及试验证明，算法设计合理，结果可靠、准确，并得出以下几点结论：

（1）本文提出的地质块体模型构建方法的核心是构建地质面片的正体和负体属性，该方法不仅建立地质块体与其组成面片的拓扑关系，同时也间接建立地质块体与地质块体之间的拓扑信息，有效弥补当前地质块体模型构建过程中未能就地质体拓扑关系信息进行建立的不足。

（2）该算法同时也兼顾了效率问题，引入方向边和方向三角形的概念，简化地质层面数据模型，使得算法无需遍历所有三角形，就能完成地质面片正体与负体属性的计算。

（3）当前的许多地质数据模型中都将面片作为模型构建元素之一，在建模过程都面临着地质体与面片拓扑关系的自动构建，因此，该算法在这些模型的构建过程具有重要意义。

［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.（熊祖强，贺怀建，夏艳华.基于TIN的三维地层建模及可视化技术研究［J］.岩土力学，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.（王纯祥，白世伟，贺怀建.三维地层可视化中地质建模研究［J］.岩石力学与工程学报，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.（魏嘉，唐杰，岳承祺，等.三维地质构造建模技术研究［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.（朱炼，唐杰，袁春风.3维地质体模型中闭合结构的提取［J］.中国图象图形学报，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.（孟宪海，杨钦，李吉刚.基于层面结构的三维闭合地质区块构造算法［J］.北京航空航天大学学报，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模式下面向拓扑的三维地质块体构建方法［J］.岩土力学，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.（李清泉，李德仁.三维空间数据模型集成的概念框架研究［J］.测绘学报，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维空间数据模型及其数据结构设计［J］.测绘学报，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.（龚健雅，夏宗国.矢量与栅格集成的三维数据模型［J］.武汉测绘科技大学学报，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.（齐华.自动建立多边形拓扑关系算法步骤的优化与改进［J］.测绘学报，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.