曹凯, 潘懋, 孙鹏
北京大学 地球与空间科学学院, 北京 100871
收稿日期：2016-09-18
基金项目：国家重大科技专项 (2016ZX05010-001)。
作者简介：曹凯 (1990-), 男, 江西上伐人, 北京大学博士研究生;
潘懋 (1954-), 男, 内蒙古赤峰人, 北京大学教授，博士生导师。

摘要：在三维结构建模中, 构建地层面时常常因不符合真实地质情况的插值点出现而产生不合理的层面厚度.针对这一问题, 提出一种厚度控制算法.首先对不规则三角网建立空间索引和OBB碰撞模型, 在此基础上进行相交测试找出不合理厚度位置, 然后采用拓扑学中几何图形剖分的思想, 在0维和1维上进行厚度调整.该方法已在石油、矿山领域成功运用, 应用结果表明本文算法求解速度快、运行稳定、结果可靠, 能自动有效地控制层面厚度, 避免了大量人工交互的厚度调整操作.
关键词：地层厚度控制三维地质建模三维结构建模不规则三角网碰撞检测空间索引
A Stratigraphic Thickness Control Algorithm in 3D Structural Modeling
CAO Kai, PAN Mao, SUN Peng
School of Earth and Space Sciences, Peking University, Beijing 100871, China
Corresponding author: CAO Kai, E-mail: caokai0831@163.com
Abstract: In 3D stuctural modeling, during construction of the stratigraphic surface, unreasonable stratigraphic thickness often arises by the interpolation point that is not consistent with the actual geological conditions. To address this problem, a new method for the thickness control is proposed. Firstly, a model of spatial index and OBB collision for triangulated irregular network is established. Secondly, according to the model, the position of the layer with unreasonable thickness is found out by intersection test. Thirdly, with the thought of geometry subdivision in topology, the thickness control is performed on zero and one dimensions. This method has been successfully applied in petroleum and mine field. The results show that the algorithm is fast, stable and reliable. It can be used to control the thickness of the layer automatically, avoiding massive manual interaction operations of thickness adjustment.
Key Words: stratigraphic thickness control3D geological modeling3D structural modelingtriangulated irregular net-work (TIN)collision testspatial index
地质曲面的构建是三维结构建模的关键技术之一.在地学领域中, 地层产状描述其在空间的产出形态, 包含其大小、形态和方位.在建模过程中, 如何约束地层面的产状, 以获取合理的地层面是决定模型准确性的关键问题.三维地质建模数据来源多样, 但是相对于复杂的地质现象其数据量仍然显得稀疏.例如测井数据, 虽然在纵向上分布密集, 但横向上的分布仅限于井点.显然, 仅靠实测的数据不足以表达地质界面形态, 通常需要使用空间插值技术获取更密的数据点^[1].常见的插值方法主要有距离反比、最近邻、克里格等^[2-4]; 然而, 上述这些插值方法仍然不能完全有效地构建符合真实地质规律的地层面产状形态 (尤其是局部产状), 并且它们仅仅关注单个层面的构建, 并不考虑多个层面之间的关系.在实际应用中, 常常导致局部地层厚度不理想, 尤其在薄层和尖灭位置.例如在不该尖灭的地方出现帖层, 产生不合理的地层剥蚀, 甚至会有地层“负”厚度出现.对此, 通常需要地质人员在三维可视化场景中对构建的地层面进行有效性检查, 然后在不合理的区域通过修改或增加约束控制点, 减少不合理的插值点出现, 直到层面符合要求.这种人工迭代建模的方式, 大大降低了效率, 而且存在遗漏的风险.随着三维地质建模技术的发展, 基于隐函数的建模方法被提出.该方法把地质界面作为在三维空间标量场的等值面, 先利用实测数据构造标量场, 然后在其中追踪等值面形成地层面^[5-6].它虽然支持薄层和尖灭, 并且在理论上保证了层面不交叉, 但是它对数据要求严格且在应用上还不成熟, 目前绝大多数的层面构建仍然采用插值的方法.可见, 层面厚度的控制是一个不容易解决的难点.
针对上述建模过程中需要人工迭代的问题, 本文提出一种厚度控制算法, 在保证其时间效率的基础上, 不仅准确检查出不合理的层厚, 而且能自动地调整该处厚度以符合要求.该算法已经成功运用于三维结构建模中.
1 算法概述地质界面用几何曲面来表达, 不规则三角网 (triang ulated irregular net-word, TIN) 能有效表达复杂不规则几何地质体表面形态, 利用不规则三角网建模是当前主流的三维地质结构模型表达方法之一.因此本文提出基于TIN实现地层厚度控制.
1.1 实现思路算法主要思路是:输入待处理的地层面 (以下称为活动面) 和层序上与之相邻接的上一个地层面 (以下称为参考面) 以及最小厚度参数 (以下称为Minthick), 找出厚度不符合此参数的局部三角形面片 (以下称为TinLocal), 保持参考面不变, 修改活动面上TinLocal的点的z值 (以z值代表层面的高程).因此, 算法包含两个关键内容:相交检测和厚度调整.然而, 在TIN中, 仅仅保证点的厚度不能避免地质界面的不合理相交, 还需要控制边的厚度.图 1a为参考面和活动面相交方式的简化模型, 虚线部分代表在参考面下方 (包括C和线段AB).如图 1b所示, 将点C提高到参考面上方C′位置后, 线段AB仍然在参考面下方.
图 1(Fig. 1)

图 1 点的厚度调整Fig.1 Thickness control of point between stratigraphic surfaces (a)—调整前；(b)—调整后.

为解决这个问题, 采用拓扑学中几何图形剖分的思想:几何体可剖分成基本组成部分 (点, 边, 三边形, 四面体, …), 剖分的基本组成成分叫做“单形”, “点”是0维单形, “边”是1维单形, “三边形”(包括内部) 是2维单形等.本文的厚度调整对象实质上是2维单形 (三角形), 要完成对2维对象的操作, 需将问题分解为对0维和1维对象的处理.因此, 厚度调整需包括两个部分:①0维点的厚度调整; ②1维线的厚度调整.当这两者都完成后, 才能实现曲面的厚度控制.
1.2 相交检测在实现厚度调整之前, 先要找出地层面不合理相交的位置.在使用TIN表达曲面的情况下, 该问题实质上是TIN的求交判断.考虑到本文算法应用于复杂大区域构造, 涉及大规模的TIN求交运算, 对算法时间复杂度有一定要求.
对于TIN的求交判断, 蛮力解法是遍历所有三角形分别进行两两求交判断, 不考虑三角形对之间的几何求交, 其时间复杂度为O(N²).减少需要求交的三角形数, 是提高求交效率的主要手段之一.因此需要在进行三角形求交之前, 尽量多地剔除那些肯定不会相交的三角形.本文借鉴空间分解^[7]的思想, 使用空间索引的方法:先建立TIN的空间索引, 然后将三角形映射到空间网格索引单元中, 仅对于同一个索引单元的三角形进行求交判断.在此基础上, 采用碰撞检测, 进一步剔除不相交的三角形对.
空间索引的目的是为了对几何对象快速定位, 以提升空间操作速度以及效率.目前主要的索引方法有R树系列、四叉树、固定格网以及K-D-B树等.本文采用格网索引, 不仅简洁易实现, 而且高效.其基本思想是把空间几何对象的外包矩形划分为若干个立方体空间网格, 然后将待处理的空间几何元素映射到立方体网格单元中, 这样就建立了空间几何元素集的索引, 方便之后操作的快速定位.假设TIN的外包矩形范围是 (X_Min, Y_Min, Z_Min) 到 (X_Max, Y_Max, Z_Max), 其中的网格单元以 (i, j, k) 来标识, 大小为cellsize, 某个三角形的外包矩形范围是 (x_Min, y_Min, z_Min) 到 (x_Max, y_Max, z_Max), 那么其映射到的网格单元范围为^[8]

(1)

(2)

(3)

碰撞检测通过利用体积稍大且特性简单的几何体来近似地代替复杂几何对象的思想来快速求得相交对象的方法.碰撞检测有多种分类, 本文采取其中性能较好的方向包围盒算法OBB (oriented bounding box).OBB最大特点是它方向的任意性, 这使得它可以根据被包围对象的形状特点尽可能紧密地包围对象, 能比较显著地减少包围体的个数, 从而避免了大量包围体之间的相交检测.在实现过程中先建立待求对象的OBB, 以之作为根节点.然后将其分解成2个OBB, 将形成的新的2个节点作为根节点的子节点.以此类推, 直到它的最小单元是不可再分割的OBB.进行相交测试时, 利用轴分离方法, 即把包围盒投影到空间坐标轴上并检查其是否线性相交, 据此判断两包围盒是否相交, 可达到O(NlgN) 的效率.
需要注意的是, 地层面之间本身有地质意义上的相接 (如不整合, 剥蚀造成的尖灭), 在建模的过程中这种地层边界上的厚度为“0”.上述相交检测会把这些合理的接触检查出来作为待处理的TinLocal, 因此需要摘除这种情况, 只留“负”厚度三角形对.Magali在文献[9]中提出了三维结构建模过程中的基本原则, 其中的拓扑一致性原则可以作为区分出“0”厚度的准则.所谓的拓扑一致性原则, 是指地质界面只能在共同边界相交, 并且在交线处必须保证点与线的一致性 (如图 2).
图 2(Fig. 2)

图 2 拓扑一致性原则[9]Fig.2 Topological consistency principle[9]

1.3 厚度调整厚度调整通过调整不合理相交位置处活动面上的点到参考面的竖直距离来实现, 在这过程中需要进行点到面的竖直投影, 即求取点到三角网的竖直方向上的距离以及投影点.目前广泛使用的是基于空间数据结构的方法, 这种方法一般采用空间剖分结构或层次包围盒结构对三角网进行预处理, 以方便快速存取操作并进行计算.对于上述两种结构有很多文献进行研究和改进, 前文所述的基于空间索引的OBB结构就是一种改进的层次包围盒结构^[10-11].
本文提出厚度调整 (即点到面的投影) 仍然采用基于空间索引的OBB碰撞检测方法.一方面, 与相交检测方法统一, 达到简洁高效的目的 (只需一次建立索引结构)；另一方面, 本文算法应用在三维地质建模中, 两个TIN面为邻接地层, 它们在同一个边界限定的区域中, 有XY方向投影重叠的特性.利用这个特性, 使用碰撞检测的方法可以很容易进行距离求取, 其具体实现方式如图 3所示, 待投影点为P, 投影面以A, B两个三角形简化表示.以P为底边中点, P_maxP_min为顶点构造竖直的两个三角形P₁, P₂, 其中P_maxP_min的Z值分别大于和小于投影面的最大和最小Z值, 以确保P₁, P₂能与投影面相交.在图示的例子中, 通过基于空间索引的OBB碰撞检测方法, 可求得P₁与投影面的三角形A, B相交.问题转化为求点P到常数个三角形 (A, B) 的竖直距离, 大大减少了几何计算量, 保证了算法的高效性.
图 3(Fig. 3)

图 3 构造三角形进行碰撞检测求点到三角网的距离Fig.3 Construction of triangle to get the distance between point to TIN by collision detection

此外, 由于地质模型中断层的存在, 投影过程需要考虑断层的影响.在某些情况下 (例如逆断层附近), 如果投影被断层“遮挡”(该方向投影到断层的距离更小), 那么该投影视为无效投影.这是因为在建模之前的步骤中断层已作为分界面约束将地层面错开, 两侧的层面关系无需处理.
2 主要流程依据前述的主要思路及关键技术, 设计算法流程图如图 4所示.
图 4(Fig. 4)

图 4 算法主要流程示意图Fig.4 Explanatory view of main process of the algorithm

步骤1??分别获得活动面和参考面的外包矩形, 按式 (1)~式 (3) 对每个三角形建立空间索引, 再对每个索引单元中的三角形建立OBB碰撞模型，并全局保存, 在后续步骤中多次使用.
步骤2??如前文1.3所述, 求取活动面的每个点到参考面的距离, 即该点厚度.提高厚度小于Minthick的点的Z值, 使其厚度达到Minthick.本步骤调整活动面上所有三角形顶点的厚度, 初步实现厚度控制.
步骤3??如前文1.3所述, 把参考面上的点对活动面进行投影, 同样求取投影距离作为厚度.如果厚度小于Minthick, 把投影点加入活动面并提高其厚度达到Minthick.如何把投影点加入活动面取决于其在三角形的位置:如果投影点在顶点, 直接处理顶点; 如果投影点在边附近, 在该点处打裂边; 如果在三角形内部, 则在该点打裂三角形 (图 5).本步骤仍然是改变活动面上点的厚度, 与步骤2区别在于该点来自于参考面点的投影, 其目的是最大程度地在0维上实现厚度控制.
图 5(Fig. 5)

图 5 投影点打裂三角形Fig.5 Split of triangle by projection point (a)—投影在顶点; (b)—投影在边; (c)—投影在内部.

步骤4??如前文1.2所述对两个面进行相交测试, 选出“负”厚度TinLocal.再对其中的每对三角形, 在XY平面上判断求交.三角形的求交转化为线段的求交, 假设线段AB与CD存在交点, 其坐标分别为A(X₁, Y₁, Z₁), B(X₂, Y₂, Z₂), C(X₃, Y₃, Z₃), D(X₄, Y₄, Z₄), 根据其线性关系可得如下方程组 (此处不考虑坐标相等, 等式无意义的情况):

(4)

求解式 (4) 可得二维交点 (X, Y, Z) 以及 (X, Y, Z′), 从而得出厚度差.若此厚度小于Minthick, 则用活动面三角形上的交点打裂边, 然后提高该点高程值达到Minthick.在某些特殊情况下, 随机取出三角形对进行厚度处理会导致产生新的负厚度, 使算法进入死循环.一个解决方法是求出所有三角形对的厚度差进行降序排序, 按厚度差递减的顺序处理三角形对.实践证明这种方法有效排除了死循环的情况.本步在1维上实现厚度控制, 结合之前的0维厚度控制, 在理论上完全实现了整个TIN面之间的厚度调整, 算法结束.
3 实验分析及应用3.1 算法分析为验证算法的效率, 挑选了3组测试样例对其进行实验, 并与不使用基于空间索引的OBB碰撞方法的普适求解算法进行了比较.实验硬件和软件环境信息如下:CPU为intel (R) Core (TM) i5 3.20 GHz, 内存为12 GB, 操作系统为Windows 8 64位, 实验平台为Visual Studio 2008, 开发语言为C++.实验结果如表 1所示.组别1和2的结果表明, 在同样的输入规模下, 相交三角形对的数量增多对普适算法时间增长幅度较小, 符合普适算法遍历所有三角形求交检测而不考虑无效求交的特点.组别1和3的结果表明, 本文算法对处理大规模TIN面其时间效率显著, 其原因是避免了大量的无效几何求交.OBB碰撞算法把相交检测从O(N²) 降为O(NlgN), 空间索引把其中N的数量级别再次降低 (具体取决于cellsize).算法时间主要消耗在建立空间索引和碰撞模型上, 而无论是对于相交检测还是厚度调整, 它们的建立只需一次.
表 1(Table 1)

表 1 厚度控制算法时间效率对比Table 1 Comparison of time efficiency of thickness control algorithm 组别 三角形数 相交三角形对数 普适算法时间/ms 本文算法时间/ms 1 7 892 8 256 1 203 16 042 5 902 2 7 924 8 301 2 291 17 930 8 723 3 82 804 77 921 1 539 > 300 000 22 943

表 1 厚度控制算法时间效率对比 Table 1 Comparison of time efficiency of thickness control algorithm

3.2 应用分析本文算法已经成功运用在Creatar三维地质建模平台中.图 6为某地区的构造模型建立过程.该地区地形结构复杂, 覆盖面积广, 是大区域复杂地质构造建模.数据来源包括地表DEM, 钻孔以及地质剖面.除地表DEM数据比较充分以外, 其他数据都相对不足, 因此在建模过程中必然产生因插值造成的不合理地层面相交 (图 6a).本文算法自动地检测并调整不合理厚度位置 (图 6b), 在可接受的时间内完成了厚度控制, 并最终形成合理的闭合地质体 (图 6c).
图 6(Fig. 6)

图 6 算法实例——某地区构造模型Fig.6 Algorithm example—a structural model of some area (a)—不合理层厚；(b)—算法调整后的地层；(c)—闭合地质体.

4 结语本文提出一种厚度控制算法, 能自动调整两个层面的局部产状, 达到合理厚度.该算法主要思想是，检查出不合理的厚度位置，然后采用拓扑学中几何图形剖分的思想分别在0维和1维的程度上对该位置的高程值进行调整, 直到厚度符合要求.算法主要包括相交检测和厚度调整, 本文提出了使用基于空间索引的OBB碰撞检测方法的思想解决这两个问题.相比于使用传统的求解算法, 在处理大规模复杂数据时大大提高了效率.该方法已在石油、矿山领域成功运用, 其应用结果表明本文算法求解速度快、运行稳定、结果可靠, 避免了大量人工交互的厚度调整操作.
参考文献

[1]	明镜. 三维地质建模技术研究[J].地理与地理信息科学, 2011, 27(4): 14–18. ( Ming Jing. A study on three-dimensional geological modeling[J].Geography and Geo-Information Science, 2011, 27(4): 14–18.)
[2]	Shepard D.A two-dimensional interpolation function for irregularly-spaced data[C]//Proceedings of the 1968 23rd ACM National Conference.New York:ACM, 1968:517-524.
[3]	Watson D F, Philip G M. Neighborhood-based interpolation[J].Geobyte, 1987, 2(2): 12–16.
[4]	Deutsch C V, Journel A G. GSLIB, geostatistical software library and user's guide[M]. 2nd ed. New York: Oxford University Press, 1998.
[5]	Hoppe H, Derose T, Duchamp T, et al. Surface reconstruction from unorgnized Points[J].Computer Graphics, 1992, 26(2): 71–78.DOI:10.1145/142920
[6]	Frank T, Tertois A L, Mallet J L. 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[J].Computer & Geosciences, 2007, 33(7): 932–943.
[7]	张少丽, 王毅刚, 陈小雕. 基于空间分解的三角网格模型求交方法[J].计算机应用, 2009, 29(10): 2671–2673. ( Zhang Shao-li, Wang Yi-gang, Chen Xiao-diao. Intersection method for triangle mesh model based on space division[J].Journal of Computer Applications, 2009, 29(10): 2671–2673.)
[8]	Lo S H, Wang W X. A fast robust algorithm for the intersection of triangulated surfaces[J].Engineering with Computers, 2004, 20(1): 11–21.DOI:10.1007/s00366-004-0277-3
[9]	Lecour M, Cognot R, Duvinage I, et al. Modelling of stochastic faults and fault networks in a structural uncertainty study[J].Petroleum Geoscience, 2001, 7(sup): 31–42.
[10]	Gottschalk S, Lin M, Manocha D.OBB-Tree:a hierarchical structure for rapid interference detection[C]//Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques 1996.New York:ACM Press, 1996:171-180.
[11]	Johnson D, Cohen E.A framework for efficient minimum distance computations[C]// Proceedings of IEEE International Conference on Robotics and Automation 1998.Leuven:IEEE Computer Society Press, 1998:3678-3683.