非拓扑一致性三维地质块体建模算法研究

罗开云， 隆振海， 何永清， 蒋先艺， 吴 恒， 宋卫锋

(1.中国石油集团 东方地球物理公司采集技术中心，保定 072751； 2.成都理工大学 信息科学与技术学院，成都 610059)

0 引言

三维块状地质建模[1]以具有近似地质属性的层块作为建模的基本单元，首先利用曲面构建方法建立地质层位和断层曲面，再利用块体追踪技术提取由层位和断层曲面片所形成的封闭地质块体，建模所得到的块状地质模型通常可以用于射线、射线束类正演模拟[2]，也可以通过离散化算法转换为规则或者非规则网格模型，支持有限差分或有限元方法[3]。三维块状建模方法符合人们探索和认知地下结构的规律，利用块体单元可以描述断裂、尖灭、透镜体、侵入体等多种地质结构，因而是目前主流的建模方法，被广泛应用于GoCAD、Petrel等著名地质建模软件中。

三维块状建模算法的历史可以追溯到20世纪80年代，Warbuton[4]提出一种基于无限大非连续面来识别块体的算法；Lin等[5]提出将拓扑概念引入到块体建模中；Lu[6]将方向边和方向环的概念引入了块体建模过程;蒋先艺等[7]提出将具有相同属性的地质层面联系在一起，并且给出了一种基于三角形的块体追踪算法，从而得出地质块体；侯卫生[8]在Lu的基础上进一步改进，使得算法可以在曲面上进行搜索从而得出块体；杨洋等[9]简化并改进了传统的线框架模型，提出了线框架拓扑模型来加速与简化拓扑搜索的过程。

经过数十年的发展三维块状建模算法日渐成熟，在块体追踪方面从简单的逐三角形追踪发展为基于线框拓扑的追踪，不仅在效率上有了很大改进，而且也具有一定的容错性。但本质上，现有的方法都是基于模型交线拓扑结构进行块体追踪，仍需要利用交线的拓扑信息。然而前期的曲面建模，特别是复杂模型的曲面建模，由于存在结构复杂、三角网尺寸差异过大、密度不均匀、浮点误差等原因，很容易造成曲面交线拓扑不一致。而当地质模型存在拓扑不一致的现象时，基于交线拓扑结构的块体追踪方法就会出现不可预期的错误，造成块体建模失败。为了提高三维块状地质建模的成功率和效率，本文提出了一种不依赖交线拓扑信息的块体追踪算法——非拓扑一致的块体追踪算法，该方法模拟人眼视觉观察，完全不依赖交线拓扑信息，仅依靠曲面自身形态进行块体三角网提取，因而具有很强的鲁棒性，对存在拓扑错误，甚至存在漏缝和交越现象的模型都能正常完成块体建模。

1 非拓扑一致性块体建模算法原理

传统的块体追踪算法以交线拓扑结构为基础，利用交线处三角形或者三角网的邻接关系，建立模型的拓扑关系图，并在此基础上通过图的广度或者深度遍历算法进行块体外边界三角网的追踪。对于拓扑关系正确的模型，利用该方法进行块体追踪效率很高，但是当模型的交线拓扑结构存在一定问题，地层或者断层三角网在交线处的拓扑划分不完全一致，或者出现局部漏缝或者交越现象时，追踪得到的块体模型就可能会出现块体缺失、相互粘连等错误。

然而面对有少量拓扑问题的模型，人们会从整体上理解模型结构，忽略局部交线处的细微拓扑问题，即使在模型存在一定问题时也能正确分辨出模型中的各个地质块体。笔者提出的非拓扑一致性块体追踪算法，模拟了人类观察和处理模型的方式，首先从整体上确定块体的大致范围，再利用局部观察的方式确定块体的边界三角网。

图1 非拓扑一致块体追踪算法流程图Fig.1 Flow chart of non-topological consistency blocking algorithm

非拓扑一致性块体追踪算法的基本流程如图1所示，算法首先读入经曲面建模后的三维曲面模型(允许模型存在一定的拓扑错误，并且无需读取任何交线拓扑信息)；接下来对地质模型空间进行均匀网格划分，将曲面三角网中的三角形信息记录在其对应位置的网格单元中；然后使用漫水网格单元追踪方法得出地质块体在空间中所占据的网格单元集合，得到地质块体的基本轮廓；最后再使用类似于人眼视觉观察的方式对块体边界区域的网格单元进行处理，从中提取出地质块体所包含的外表面三角形，从而构造地质块体外表面三角网。

为了直观地说明非拓扑一致性块体追踪算法的关键步骤和执行过程，图2给出了算法执行的具体过程，如图2(a)所示算法首先读取模型的曲面三角网，确定模型的空间范围，将模型空间划分成均匀网格单元，然后用漫水法追踪块体的网格单元集合得到如图2(b)所示的由网格单元所表示的基本块体模型，最后再通过模拟视觉观察方式提取块体的外边界三角网，得到最终的块体模型如图2(c)所示。

2 漫水法追踪块体网格单元

非拓扑一致性块体追踪算法首先要确定块体的基本形态，并大致确定块体可能包含的三角形，为之后精确提取地质块体外表面三角形做好准备。为此可以先将模型空间划分成等距均匀的结构化网格，并将曲面模型中的三角形根据坐标位置放置到与之有交的网格单元中(网格单元是一个长方体)，然后就可以从一个不含有任何三角形的网格单元出发，利用漫水法追踪块体所包含的网格单元集合。

漫水法是图形学中常用的一种区域填充法[10]，其实质是一种广度优先遍历的种子填充法，适用于对内定义区域的填充。而图形学中的内定义区域，是指区域内部的像素具有类似或相同的属性，而区域外的所有像素具有不相同或差别较大的属性，该算法从内区域中的某个像素出发(通常称起始像素为种子)，利用广度优先搜索完成内区域填充。可以对该方法进行改造，用于三维空间中地质块体网格单元集合的追踪，三维空间中地质块体的漫水网格单元追踪算法的具体步骤如下：

图2 非拓扑一致块体追踪算法效果示意图Fig.2 Result of non-topological consistency blocking algorithm(a)A处细节;(b)B处细节;(c)C处细节

1)尝试从三维地质空间中取得一个未被访问过且不含有三角形的网格体作为起始网格体，若不能找到一个所述不含有三角形且未被访问过的网格体则块体追踪过程结束，否则将所述起始网格体加入到遍历队列中并进行步骤2)。

2)若遍历队列为空则重新进行步骤1)，否则从遍历队列中取出一个网格体进行步骤3)。

3)对所述取出的网格体，若其已经被访问过，则重新进行步骤2)，否则进行步骤4)。

4)对所述网格体，找到其周围的六个邻接网格体，对每一个邻接网格体，更新所述邻接网格体的投影方向。

5)对所述邻接网格体，若其中不含有三角形，则将其加入到遍历列表中，并将其标记为已访问，之后返回进行步骤2)，否则将其加入到当前块体的边界区域网格体队列中，再进行步骤2)。

利用上述漫水追踪算法，可以提取模型中各个地质块体所包含的网格单元集合，这些网格单元确定了地质块体的基本形态，而且网格单元集合中包含位于块体边界上的网格单元，这些边界网格单元包含曲面三角形，而其中部分三角形是属于某个特定地质块体的，这为最终提取块体外边界三角形奠定了基础。

3 基于视觉观察的块体三角网建模

非拓扑一致的块体追踪算法的最终目标是提取地质块体的外边界三角形，构造地质块体表面三角网，通过漫水追踪算法获取了地质块体的网格单元集合，为了进一步提取块体的外边界三角形，需要遍历块体的网格单元，从含有三角形的边界网格单元中找出属于与地质块体的三角形，并且将它们组合起来构造块体外表面三角网。

图3 基于视觉观察的块体三角网建模流程图Fig.3 Flow chart of blocking algorithm based on visual observation

从边界网格单元中提取三角形这一过程需要基于模仿人眼视觉观察的方式来进行，方能达到类似于人眼“所见即所得”的观察效果。对于含有三角形的边界区域网格单元，通常它所含的三角形一部分属于当前地质块体，另一部分其他块体，考虑到一个边界网格单元的6个面中至少有一个是与内部网格单元相邻的(内部网格单元指不含有任何三角形的网格单元)，可以将与内部网格单元相邻的面作为观察面，通过该面观察边界网格单元中所包含的三角形，属于当前块体的三角形必然能够被观察到，而不属于当前块体的三角形则会被遮挡而无法观察到，通过“视觉观察”就可以准确提取边界网格单元中属于当前地质块体的三角形了。

当然计算机并没有眼睛，为了实现上述“视觉观察”算法，可以利用OpenGL进行三维图形渲染[11]，并对渲染得到的图像进行判定，从而实现所谓的视觉观察效果。首先确定观察方向，并建立一个与其同位置，同大小的三维平行视景体[11-13]；然后将该网格单元所包含的三角形赋予不同颜色，再利用OpenGL三维图形库将场景渲染到帧缓存对象中[14]；最后从帧缓存中将渲染结果图像取出，遍历图像中每个像素点的颜色，选择在渲染结果中占据了一定像素比例的三角形作为当前地质块体的外边界三角形之一。

对上述算法思路进行总结，可以得到如图3所示的基于视觉观察的块体三角网建模流程图。在地质块体外边界三角形的提取过程中，每一个边界区域中的网格体均会被施以如下处理：

图4 三维地质模型示意图Fig.4 3D geologic model

1)遍历一次所述网格体中已保存的所有三角形信息，遍历过程中检查所述网格体中的三角形是否来自同一个地层面三角网。

2)若所述网格体中的三角形均来自于同一个地层面三角网，则直接将这些三角形判定为当前块体外边界，否则进行步骤3)。

3)创建一个颜色表，所述颜色表中不存在两个相同的颜色，并使颜色表中的颜色与三角形一一对应。

4)使用一个与所述网格体同大小且同位置的视景体进行平行投影绘制得到绘制结果，其中每个三角形均使用其对应所述颜色表中的颜色。

5)利用所述视景体对所述网格体中含有的三角形进行平行投影绘制得出绘制结果，其中每个三角形均使用颜色表中已记录的一个颜色，三角形与颜色表中的颜色一一对应。

6)统计出绘制结果中每种颜色绘制结果中的面积占比，所述颜色对应面积占比大于2%的颜色认定为有效颜色，有效颜色对应的三角形认定为有效三角形。

7)将有效三角形判定为当前块体的表面三角形。

4 算法应用实例

为了验证非拓扑一致性算法的有效性，我们选择了一个13.5 km×13.5 km×8 km较为复杂的地质模型进行算法测试，如图4所示，该模型含有44个层位/断层曲面，存在断裂、尖灭、透镜体等多种地质现象。为了验证算法处理非拓扑一致性模型的能力，在图4中(a)、(b)、(c)三处人为地制造了一些曲面漏洞和交越问题，图5给出了上述三处交线错误的放大图，从图5中可以看到模型存在比较明显的拓扑问题。

图5 三维地质模型细节示意图Fig.5 Details of 3D geologic model(a)A处细节;(b)B处细节;(c)C处细节

图6 传统追踪方法结果图Fig.6 Result of traditional blocking method(a)非拓扑一致性模型结果;(b)拓扑一致性模型结果

图7 本文追踪方法结果图Fig.7 Result of our blocking algorithm(a)整体结果; (b)追踪结果C处细节图

图6(a)是利用传统的拓扑一致性追踪方法建立的块状模型，图6(b)则是用传统方法对原始没有拓扑问题的模型进行追踪建立的块状模型。从图6中可以看出，当模型存在拓扑问题时，传统方法得到的块体模型结构存在问题，局部存在拓扑错误的块体粘连在一起，得到的块体与原始模型追踪得到块体结构和总数存在差异。

图7(a)是利用本文提出的非拓扑一致性算法建立的块状模型，在Intel i7 4790K主频4.0 GHz的测试电脑上，建模时间为23 s。对比图7(a)和图6(b)可以看出，当模型存在一定拓扑问题时非拓扑一致性算法得到的块状模型整体结构正确，块体结构和数目与图6(b)一致。图7(b)给出了非拓扑一致性算法在图7(c)处得到的块体模型，因为原始模型在交线处存在漏缝，所以得到的块体模型也存在一定的缺漏，但是这并没有影响模型构建的整体结果。说明非拓扑一致性算法具有很好的容错能力，可以处理存在一定拓扑问题的模型。

5 结论

三维块状地质建模目前主流的地质建模方法，块体追踪技术是该方法的核心技术，现有的块体追踪技术均依赖于模型的交线拓扑结构，对模型拓扑的一致性和正确性要求较高。笔者提出的非拓扑一致性块体建模算法，首先利用漫水法确定地质块体的基本形态，然后利用基于OpenGL的视觉观察技术实现块体外边界三角网的提取与构建。该算法摆脱了对交线拓扑信息的依赖，仅根据曲面模型自身形态即可完成块体构建，因而对曲面模型的正确性要求大幅降低，当曲面模型中存在交越、漏缝等现象时依然可以正确完成块体模型的构建，与传统块体追踪技术相比，该方法降低了输入模型的要求，具有很强的鲁棒性。虽然该算法运行的效率低于传统的拓扑一致性追踪方法，但是降低了输入曲面模型的质量要求，因而降低了建模的整体难度，能够更好地辅助用户完成复杂三维地质模型的构建，提高建模效率，具有比传统方法更好的实际推广使用价值。