1 引 言

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

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

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

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

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

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

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

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

为了描述算法原理的方便，引入以下概念。

悬挂面片：正体和负体是同一个体的面片，类似二维空间中的悬挂弧段。

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

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

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

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

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

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

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

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

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

图选项

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

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

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

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

具体步骤如下：

(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)。

struct thirdverd    //方向三角形

{

        public long PointID;//第三点ID，方向三角形的另两个点为方向边上的端点

        public long FaceID; //方向三角形所在面片ID

}

struct diredge    //方向边

{

        public long ID; //方向边所在的边界弧段ID

        public long[] DiredgePointID;//方向边上起点和终点ID

        public thirdverd[] Direction_Tri;//排序后的方向三角形(即面片排序表)

}

struct dBorderLine    //边界弧段

{

        public long ID; //边界弧段ID

        public long startID, endID; //边界起点和终点ID

        public long [] allpointIDs; //边界弧段上所有点ID

        public int N; //边界弧段的重度

}

class Face //面片

{

        public long ID; // 面片ID;

        public dBorderLine [] BorderLines; //边界弧段ID

        public long[] TriangleIDs; //面片包含三角形ID

        public long nomalbody=-1; // 正体ID，初始值为-1

        public long countbody=-1; // 负体ID，初始值为-1

}

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

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

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

图选项

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

图选项

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

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

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

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