三维地质界面网格模型构建算法研究

侯 晓 琳

(北京大学地球与空间科学学院 北京 100871)

0 引 言

在地质学相关工程中通过多种技术手段获取的多类型地质数据不能够直观展示真实地质情况。三维地质建模技术利用地质数据建立三维地质模型，采用计算机可视化方法更加直观地模拟出地下构造情况，成为地质学研究重点[1-2]。目前针对于三维地质建模技术的研究成果已经逐步应用于生产实践项目中，多种三维地质建模软件尝试为石油勘探开发、矿产资源预测等项目提供三维定量化模型，为实际开发过程提供理论依据，降低开发成本[3]。

根据不同的地质数据类型[4]，三维地质模型可以分成两类[5]：地质构造模型，通过解释原始地质构造数据，描述地层、断层等地质界面空间形态展布和拓扑关系[6]；地质属性模型，采用分析预测方法，定量化描述地质变量空间变化规律[7]。地质构造模型是三维地质模型基础[8]，其研究内容包括地质体模型建模方法、地质层面建模方法以及地质要素拓扑关系描述方法。从数据结构角度，将地质构造模型分成面模型、体模型以及混合模型。通过建立三维地质界面模型描述复杂地质构造，以地质界面模型为约束建立三维地质体模型。三维地质层界面模型是描述地质要素之间拓扑关系的基础[9]。

本文在现有研究的基础上，提出一种改进的三维地质界面三角化网格模型构建算法，采用实际的断层面数据建立并可视化展示断层面三角网格模型。与现有的基于映射关系建立三维三角化曲面模型算法相比，改进后的算法减少了映射关系计算过程，避免了可能出现的三维空间中多点映射到二维平面上一点而导致无法生成三维曲面模型的问题，改进的局部映射方法可在三维空间中直接生成曲面网格模型，该算法同样也适用于解决其他领域三维三角化曲面构建问题。

1 相关研究

三维地质构造模型的几何模型包括点模型、线模型、面模型和体模型。面模型用于表示三维地质体界面，反映界面间构造关系，例如地层面、断层面等。面模型主要类型包括不规则三角面模型、规则网格模型、边界表示模型和线框模型等[10]。

在三维地质体构造建模中，通常采用规则网格和非规则网格描述地质界面和地质体。规则网格优点在于网格剖分算法较为简单，便于进行数值模拟计算。非规则网格则能够更好地模拟地质界面状态和地质体构造形态，反映出地质构造的非均质性[11]。与规则的四边形网格相比，非规则三角化网格在描述曲面和不规则几何体时具有一定优势。在三维地质体建模中，表示地质界面的原始数据为点集数据，例如，断点数据描述地质断层界面，分层点数据表示地质层面。不失一般性，三维地质界面三角化网格模型构造方法同样可以应用于基于点数据集构建三维曲面模型的相关问题。

现有研究中主要采用映射法建立三维曲面网格模型，将三维数据点投影到二维平面上，在平面上进行二维三角剖分，根据点映射关系，将二维平面剖分结果映射到三维空间，生成三维曲面三角化网格模型[12]。在映射法中，映射函数和映射平面的选取方法具有一定不确定性，映射函数的影响生成三角化网格模型的形态[13]。很多非规则网格生成算法均采用映射法实现思路。经典的三角化剖分算法有Lawson算法和Watson算法等[14-15]，这些算法实现思路清晰，并且能够保证三角形网格质量。但是，映射法在某些特定情况下也无法处理。三维点集中任意点p(x,y,z)，必须保证其在物理域和映射域中是一一映射关系。假设三维空间上的两点p1(x,y,z1)和p2(x,y,z2)，映射到某一二维平面上的点为同一点p′(x,y)，在这种情况下映射法就很难处理，确定映射函数和选择映射平面是基于映射法建立三维曲面模型的关键问题。当三维空间中的不同点映射到二维平面的同一点时，主要解决方法是将曲面点集划分多个子集，基于不同子集建立三角化网格后再进行拼接，可能导致出现各部分曲面边界的拓扑关系不一致问题。确定映射平面、获取点集包围盒等子问题增加映射法实现难度。

本文提出一种改进的基于边界控制-局部映射的曲面三角化网格生成算法，在三维空间中基于点集构造三角化网格模型，直接剖分算法减少映射算法中数据预处理等计算过程，通过分析三维空间点位置关系生成三角化剖分模型。在存在边界条件约束和无边界条件约束下，均可以建立三维地质界面三角化网格模型，客观反映地质界面展布形态，避免了三维空间中多点映射二维平面同一点的多对一映射问题。

2 三维空间曲面三角化网格模型

2.1 Delaunay三角化

Delaunay三角化特征能够保证三角化网格模型具有良好的构网形态[16]。Delaunay三角化的三个重要性质描述如下[17]：

(1) 空外接圆(球)性质。在二维空间平面Delaunay三角网中每个三角形的外接圆均不包含其他节点。三维空间中，每个四面体的外接球不包含其他节点。

(2) 最小角最大化性质。满足所有三角形中最小角的最大化性质。

(3) 局部优化性质。三角网局部优化过程能够保证全局优化。

在二维平面上可以证明能够构建Delaunay三角化网格模型，满足Delaunay三角化特征。基于映射法，在二维平面构建的Delaunay三角化网格模型映射到三维空间时，Delaunay三角化特征无法描述。但其几何定义是保证三角网格质量理论依据。因此可以将Delaunay三角化特征应用于建立三维曲面三角化网格模型，控制三角网格形态。

2.2 曲面三角化模型构建算法

基于Delaunay三角化特征提出一种边界控制-局部映射的曲面三角化网格模型构建算法，算法流程见图1，主要包括三个步骤：(1) 三角化网格模型初始化；(2) 三角化网格模型更新；(3) 网格模型边界约束更新。

图1 曲面三角化网格模型构建算法流程

断层的断点点集S={p1,p2,…,p25}作为实验数据，详细说明边界控制-局部映射曲面三角化网格模型构建算法步骤。假设算法输入为点集S，算法输出为边集合E以及三角形网格模型T。

2.2.1三角化网格模型初始化

算法首先从三维点集合中挑选三个点构成网格模型的初始化三角形。以初始三角形为中心，不断加入新的三维点更新网格模型。曲面空间中三角形由三边构成，将边定义为两种类型：内部边和外界边。在曲面三角网格模型T中：(1) 内部边：一条边属于两个不同的三角形；(2) 边界边：一条边属于且只属于一个三角形。算法需要区分并记录边类型，边界边集合表示为borderset，边界边集合内所有边构成闭环。基于边界边拓展网格模型，基于内部边更新网格模型，同时保证三角化网格质量。

采用随机方式在点集中确定第一个三角形，该三角形的三条边均属于边界边，记录边界边信息，图2显示了断层断点数据集中随机初始三角形结果。

图2 基于点集S初始化种子三角形示意图

三角化网格模型随机初始化三角形思路：(1) 在点集S中，获取随机点p9，根据最近邻原则选取最近点p4，构成第一条三角形边e1=(p4，p9)，将p4、p9标记成已访问状态，以标记点集S′=(p4，p9)，点集S更新为S-S′。(2) 根据点到边e1中点的最小距离确定三角形第三点。点p到边e中点距离表示为dist(p,e)=dist(p,pm)，pm是e的中点。因此，p20作为距离边e1中点最近点构成第一个三角形，并将其标记为已访问，更新边界边集合borderset和三角化网格T：

borderset={e1=(p9,p4)，e2=(p4,p20)，e3=(p20,p9)}

T={t1=(p9,p4,p20)}

2.2.2三角化网格模型更新

三角网格模型更新过程为从点集S中选取未标记的三维点加入到三角化网格模型中，更新边界边集合borderset和三角网模型T，直至点集S中所有点都加入到三角化网格模型，如图3所示。

图3 三角化网格模型更新流程

计算点集S中的点到borderset集合中边界边中点的最小距离。以图2中断点三角化网格为例，计算可得dist(p14,e1)≤dist(p′,e1)，dist(p15,e2)≤dist(p′,e2)，dist(p14,e3)≤dist(p′,e3)，∀p′∈S，且dist(p14,e1)

1) 局部映射。假设曲面连续光滑，基于Delaunay三角化特征采用局部映射方法判断点与三角形的位置关系，更新三角化网格模型，确保三角网格良好形态。

采用局部映射方法判断空间一点与三角形t位置关系。假设任意三角形三点坐标为p1={x1,y1,z1}，p2={x2,y2,z2}，p3={x3,y3,z3}，以其中一边(p2，p3)为标准，过边(p2，p3)垂直于三角形的平面方程A(x,y,z)，即：

式中：

通过垂直平面方程，以边(p2，p3)为标准，将点p0={x0,y0,z0}和三角形第三点p1坐标代入A计算可得A1和A2。任意一点p0与三角形t的空间位置关系定义为：

① 若A1×A2<0，p0位于以边(p2，p3)为标准三角形的负向空间；

② 若A1×A2>0，p0位于以边(p2，p3)为标准正向空间；

③ 若A1×A2>0，p0在三角形所在平面的垂直平面上，边(p2，p3)也在该垂直平面上。

判断一点p0在三角形p1p2p3平面的局部投影是否在该三角形范围内，以三边(p1，p2)、(p2，p3)、(p1，p3)为标准，若p0均在各边的正向空间中，则表示该点局部映射在三角形内，反之，则表示局部映射不在三角形范围内。

图4 三角形与点空间位置关系示意图

根据三种空间位置关系，采用三种不同方法更新三角化网格模型：

(1) 点pt位于边界边e负向空间时，点pt与边界边e构成新的三角形加入网格，边界边e变为内部边，并增加边界边e1=(p2,pt1)和e2=(pt1,p3)。

在该种情况下，点pt在平面A上的投影与边界边关系可为三种，见图5，局部投影点可能位于边界边e上、边界边e沿p2p3方向的延长线上、边沿p3p2方向的延长线上。

图5 边界边e与在平面A上点的局部映射点位置关系

在图5中根据不同的位置关系生成不同三角形，更新边界边集合，如图6所示。

图6 根据边界边上映射点位置关系更新网格模型

图7 边界边与其正向空间映射点关系

根据正向空间内不同的位置关系采用不同方法生成新三角形：

图8 边界边与前序边界边最近点相同情况下生成新三角形

图9 边界边与前序边界边最近点不相同情况下生成新三角形

2) 最小角最大化原则。基于Delaunay三角化特征保证三角化剖分质量，引入最小角最大化对三角网格模型进行适当变化调整。在图10中，点p4与边界边e组成新三角形t2=(p4,p3,p2)，共边e的另一三角形为t1=(p1,p2,p3)，计算两个三角形的最小角。对换对角线(p1,p4)，计算两个三角形最小角值，与变换前最小角值比较，若变换后最小角值较大，则进行对角线变换，变换后三角形为t2=(p4,p3,p1)、t1=(p1,p2,p4)。当三角形对角线产生对换后，同样检验其相邻三角形是否需要进行对角线兑换，保证最小角最大原则。

图10 三维空间三角形对角线变换

保证最小角最大化原则，更新三角形网格模型主要流程为：

(1) 在网格模型中每加入一个新三角形，采用队列q处理三角形对角线变换。对于每个新三角形，将其内部边放入队列构成队列初始状态。

(2) 依次从队列取出一条内部边e，直至队列q为空结束，否则，继续执行(3)。

(3) 计算内部边e所属两个三角形t1、t2，根据最小角最大化原则进行判断。若需要对角线变换，则交换两个三角形对角线，并将t1、t2中其他内部边加入到队列q中，继续执行(2)。

以图2中点集S为例，不断加入未标记点更新网格模型，当网格模型中三角形个数为12时，三角网格模型如图11所示。

图11 点集S三角网格模型更新过程中示意图

此时边界边集合borderset和三角网格模型T更新为：

borderset={(p9,p4)，(p4,p5)，(p5,p10)，(p10,p15)，(p15,p20)，(p20,p25)，(p25,p24)，(p24,p18)，(p18,p13)，(p13,p8)，(p8,p14)，(p14,p9)}

T={(p9,p4,p20),(p4,p15,p20),(p4,p10,p15),(p4,p5,p10),(p9,p20,p14)，(p14,p20,p25),(p14,p19,p25)，(p19,p24,p25),(p13,p24,p18),(p13,p19,p24)，(p13,p19,p8),(p3,p19,p14)}

当点集S为空时结束三角化网格更新，记录下边界边集合，下一步进行网格边界约束更新。图12展示图2中所有点加入到三角化网格模型，显然三角化网格模型构成凹多边形，进一步更新网格模型边界，采用凸多边形模拟三维曲面边界。

图12 所有点加入三角化网格模型示意

2.2.3网格模型边界约束更新

根据边界约束更新三角化网格模型，主要有两种方法：(1) 无边界数据约束情况下，自动填充边界，原则上自动填充为凸多边形边界。其基本思路是按照边界边集合中边界连接顺序，沿着边界边顺序逐一更新模型，直至边界边无需更新操作为止。(2) 根据已有的点集S曲面边界数据约束更新三角化网格模型。

在边界边集合中任意边界边e都有与其相连的前边界边ebefore与后边界边enext。在有曲面边界数据约束的情况下，若边界边e被确定为曲面边界，则无需处理；反之，则对当前边界边进行更新，判断边界边e能否与ebefore或是enext构成新三角形，更新模型边界。根据三边关系更新网格模型边界基本思路如下：

① 若ebefore和e同属于一个三角形：ebefore和e无法构成新三角形，同样，若e和enext同属于同一个三角形，则e和enext无法构成新三角形。

② 若ebefore和e不同属于一个三角形：如果没有点的局部投影在边ebefore和e所在平面的三角形内，则ebefore和e可以构成三角形，反之，则无法构成三角形。采用同样方法判断enext和e是否能够构成三角形。

③ 若只有边ebefore和e能构成三角形，将新三角形加入到网格模型，删除边界边ebefore和e,将新三角形第三边加入边界边集合继续边界更新操作；若只有enext和e能够构成新三角形，加入网格模型，删除边界边enext和e，第三边加入到边界边集合；若ebefore和e，enext和e都能构成新三角形，根据两个三角形的最小角，选择最小角最大三角形加入到网格模型中。每增加一个新三角形到网格模型时，基于最小角最大化原则，局部更新网格模型，保证网格形态。

在无曲面边界约束数据的情况下，图2中点集S进行模型边界更新后，三角化模型展示如图13所示。

图13 点集合S曲面三角网格模型

3 地质界面三角化网格模型实例

地质界面数据常采用点集合数据表示，例如断层采用断点数据描述，地层采用地层点数据描述。本文提出一种基于点集数据的局部映射-边界控制三角化网格模型构建算法。

以断点数据集为例说明算法实现流程与效果。在三角化网格模型构建算法流程中，随机初始化三角形网格模型，依次加入断点集合中其他点数据，更新三角化网格模型，直至断点数据均加入到网格模型结束更新，其构成的三角化网格模型如图14所示。

图14 加入所有点数据后生成的三角化网格模型

基于图14中的三角化网格模型和边界边集合，更新曲面三角化网格模型边界。在无曲面边界数据约束的情况下，填充边界三角形，根据边界边集合自动生成凸多边形曲面，如图15(a)所示。当加入边界约束数据时，恢复模型边界，如图15(b)所示。

(a) 无边界约束边界自动更新 (b) 加入边界约束数据模型图15 某断层曲面三角化网格模型

图15所显示的断层三角化曲面网格模型中，存在两点距离较近和两点映射到某平面共线等情况。若采用映射法，选择不合适的映射平面会导致三角化结果错误，映射平面直接影响网格模型生成质量。直接在三维空间中通过空间位置判断生成三维曲面三角化网格模型，能够较好地还原曲面的真实形态。虽然空间位置判断较为复杂，实际上，在算法实现过程中减少数据预处理和其他的计算过程。

图16展示多个断层面三角化网格模型形态，在生成多断层面三角化网格模型时，算法优点明显，与映射法不同，不需要计算每个断层面的二维映射平面，而是直接在三维空间中构建多个断层面的网格模型，减少算法其他辅助计算步骤。

图16 多个断层面三角化网格模型

与经典映射法相比，本文提出的曲面三角化网格模型构建算法能够直接应用于构造较复杂的三维曲面模型。图17、图18展示了基于点集数据采用局部映射-边界控制算法生成的较为复杂三维曲面模型，经典映射法较难实现，原因在于很难确定一个映射平面，将三维数据点映射到二维平面上，进而较难生成准确的三维曲面三角化模型。

(a) 斜视图

(b) 俯视图

(c) 侧视图图17 复杂三维曲面三角化网格模型示意图

(a) 斜视图

(b) 侧视图

(c) 俯视图图18 复杂三维曲面三角化网格模型示意图

4 结 语

在三维地质建模中，基于地质构造模型研究地质界面和地质体构造形态具有重要的意义，构造三维地质界面三角化网格模型作为地质体构造约束，是研究地质构造的基础。本文基于现有曲面三角化网格构建算法研究分析，结合地质数据特点，提出一种局部映射-边界控制的三维三角化网格模型构建算法。基于点集数据构建地质界面三角网格化曲面模型，模拟地质界面特征。采用断层点集数据验证算法有效性，可视化展示断层三角化网格模型。与经典映射法相比，局部映射-边界控制三角化网格模型构建算法能够在三维空间中直接基于点集数据建立曲面三角化网格模型，无需对点集数据进行预处理，避免映射法三维空间中多点映射到二维平面一点所引起问题，并且能够保证较好的网格模型形态。该算法可以拓展到其他领域用于建立三维曲面网格模型。