基于八叉树和混合搜索树的地质曲面快速求交方法

孙黎明 魏迎奇 蔡红 严俊 宋建正 乔芸芸

摘要：为处理地质界面之间的空间相交关系，提出一种新的针对三角地质曲面的快速求交方法。该方法融合优化八叉树法和OBB搜索树方法，可以更快速准确地剔除远离交线的其他三角形。求交剩余的三角形得到交线，应用三角网局部重构和网格优化算法修正交线附近的三角网，最终分割交线两侧的地质曲面，完成2个地质曲面的离散化求交过程。与AABB、OBB和空间分解法相比，该方法在大数据量三角曲面求交中效率优势明显，可以快速准确处理地质模型构建和分析中的曲面求交问题，为三维地质模型自动化构建的实现提供有效支撑。

关键词：地质建模；地质曲面；曲面求交；八叉树；OBB搜索树；三角网重构；模型切割；区域离散化

中图分类号：P221.1

文献标志码：B

0 引 言

三角形的地质曲面是构成三维地质模型的基本几何元素，在大区域复杂地质模型的自动化构建中，通常需要处理地质界面之间的空间相交关系，这是三维地质建模应用研究的热点问题之一。[1-6]在建模过程中，不整合地层边界计算、断层面与地层的交线和任意切割面与地层的交线等[7-8]应用广泛。三维地质模型的分析应用，例如剖面制作、巷道模拟和地下工程等[9]，需要准确快速得到相交区域的空间位置、拓扑信息和内部构造，以对模型区域有更准确的理解，同时也需要进行多次曲面求交和计算分析。[10]随着地质建模数据量越来越庞大，尤其是地震数据重构断层和地层时，三维地质建模需要一种快速稳定的三角曲面求交算法。地质曲面的求交过程根据2个已有的三角曲面集合，得到相交曲面之间的交线，并重构2个三角曲面和交线两侧的三角形，完成曲面的离散化。如果循环判断求交，那么所有的三角形都要参与求交计算，算法的效率为O（n2）。

实际进行求交的三角形对直接影响三角曲面求交速度，因而算法的核心是快速准确地找到与交线相关的三角形。目前，曲面求交和物体碰撞检测的方法主要有空间分解法（八叉树、BSP树和KD树）和层次包围盒法（AABB、OBB、球體包围盒和k-DOP）等。地质曲面求交通常直接采用物体碰撞检测的方法，LINDENBECK等[11]开发TRICUT软件，用Rapid库解决三角网求交问题，应用广泛。CGAL计算几何算法库改进了AABB树求交方法 [12]，基于OBB的地质曲面切割方法也广泛应用并取得良好效果[13]，但在数据量较大时，构建层次搜索树十分耗时。YANG等[8]提出一种基于混合包围盒的针对真三维模型的切割方法，可解决多种几何体一体切割的问题。ELSHEIKH等[14]提出一种基于追踪法的可靠地质曲面求交方法，需要利用三角网拓扑追踪曲面交线，但寻找起始追踪点较困难。

本文提出一种改进的地质曲面快速求交方法，用结构化的八叉树算法与改进的OBB混合搜索树相结合的方法快速求取三角曲面交线，并重构交线所影响的三角形，然后将算法应用于地质建模分析的典型问题中，包括曲面切割地质模型，断层面求交和地质区域离散化等。为将复杂地质曲面求交问题简化，将地层求交分解为任意两两三角曲面求交。

1 理论方法和求交步骤

检测2个三角曲面的相交部分，核心是快速找到相交的三角形。无论曲面三角形数量多大，实际相交的三角形都很少，直接构建层次搜索树效率很低，因此快速剔除大部分不相交的三角形，只留下交线附近区域的三角形，然后再构建层次搜索树，可大大提升求交速度。采用结构化的八叉树数据结构高效剔除求交曲面中的大多数不相交三角形，构建基于改进OBB的搜索树，快速准确地求取交线，改进曲面求交方法流程见图1。采用静态结构化数据结构[15-16]表达1个三角曲面，即用1个大于三角网的六面体等分后表达地质界面和2个曲面相交的三角形，建立特殊的层次二叉搜索树，根节点根据三角形的AABB表示，叶节点则用OBB代替，最后对三角形求交得到交线。

1.1 优化八叉树的粗略求交

空间八叉树是二维平面中四叉树在三维空间的扩展，是用均匀细分空间描述三维场景中几何对象的一种树状数据结构，其每个父节点用一个正方体表示，每个父节点都有8个子节点，是将父节点均匀细分为8个小正方体得到的。八叉树通常用于三维空间中海量物体的快速查询和空间索引，目前广泛应用于大数据量的三维GIS，包括三维智慧城市、三维数字管线、三维BIM可视化和3D大场景游戏等。

八叉树的特点是可以用多层树快速索引每个节点，但对于海量数据离散点构建的地质曲面和细分到单一节点的八叉树，每次计算的迭代数量庞大，且八叉树的层级非常多，计算效率更低，每次循环计算耗时很长。本文提出一种非对称的空间八叉树模型，用于快速查找2个三角曲面交线所在的三角形。与传统空间八叉树模型不同，该模型不是每层都按照8n细分，而是达到一定的边界条件就采用结合OBB的方法进行求交。设计八叉树的边界条件，计算八叉树停止的条件，以省略八叉树细分到叶节点的过程，提高求交的计算效率。优化八叉树求交过程见图2。

地质曲面所在的三维笛卡尔坐标系定义为R，定义三角网中的点VR、三角形TR、地层三角网SR。任意点P（x， y， z）对应立方体ER，立方体的边长为三角网包围盒在x、y和z这3个方向上长度的最大值l。依次将每个父节点分为8个子节点，叶节点满足与同行列子节点内的曲面点个数接近条件，并且当个数最小时八叉树就不再细分。与交线不相关的不相交三角形见图3，剔除远离交线后的三角网见图4。

利用八叉树结构计算1个子节点内的三角形，可以很快找到距离另一个曲面最近的三角形集合（见图3）。显然，该方法得到的结果仍然包括许多不相交的三角形，如果单用此算法求交，效率不高，因此本文结合改进OBB层次搜索树的方法进一步剔除不相交的三角形。

1.2 OBB混合搜索树的精确求交

OBB是1个紧密的包围盒，是包含该物体且相对于坐标轴任意方向的最小六面体，可提高检测效率。2个相交地质曲面的OBB包围盒见图5。

构建OBB搜索树和计算OBB比计算AABB更繁琐耗时。根据几何坐标计算物体的OBB，因为三角形OBB需要多2个向量计算，所以大数据量时直接构建OBB搜索树复杂低效。GOTTSCHALK等[17]提出一种有效的实现方法，采用AABB构建中间所有层次的根节点，将第1.1节中计算的六面体包围盒直接作为AABB，而叶节点采用OBB方法，OBB混合搜索树见图6。

搜索树采用二叉树结构存

储，最终叶节点采用OBB求交，能最大限度减少精确判断相交的三角形对。综合利用2种包围盒的优点，可减少计算量，加快搜索树的构建速度。[18]

根据可能相交的三角形构建层次搜索树的简要步骤如下：

（1）选择根节点。将所有的三角形最大包围盒投影到3个坐标轴，找到3个方向中的最长轴，选择其中间值为最大值与最小值之和的1/2。

（2）构建搜索树。根据平均三角形AABB的大小，选择树的深度层次。树的深度选择是影响搜索求交速度的重要因素。经过第1.1节剔除工作后，求交的三角形只占曲面三角形总数的小部分，层次不需太深。

（3）递归分割。同时将每个三角形和点的包围盒自顶向下遍历并存放指针到每一层合适的位置节点，直到单一三角形为止，向下一直到叶节点都是三角形，叶节点用OBB代替AABB。

（4）求交。自顶而下判断包围盒是否相交，如果递归到叶节点，根据三角形的3个点构造OBB，采用SAT方法求交。[19-20]

1.3 交线区域网格重构

对得到的所有三角形对逐一求交，即可得到最终的交线。一对空间三角形相交的结果有3种情况：点、线段和相交三角形。地质曲面求交中最常见的是线段，线段首尾相连形成交线。交线可以有多条，也可能是封闭环。目前，空间三角形求交的算法主要分为标量判别法和矢量判别法[21]，本文选择适于地质曲面求交的MLLER[22]三角形求交法。

加入新的交线后，将交线点分别加入到交线两侧的三角网中，旋转变换（不是投影）到合适的二维平面。[23]采用局部三角形修复方法，若交点在三角形的边上，则在此处添加新点，将这条边分割；若交点在三角形内部，则在内部添加新点，连接顶点构成新三角形，尽量满足Delaunay法则，加入新点后的三角网重构见图7。

新的三角网交线区域必须进行网格质量优化。若新加入的点形成狭长三角形和拓扑错误三角形，在地质模型空间分析、体元划分和数值评估中都无法满足要求，则需要进行三角网质量优化，局部修改

重构。根据交线点所在的三角形的网格拓扑找到周围三角形[24-25]，优化网格质量，包括边界交换、合并三角形、点重构和插入新点。[26]

将局部修复好的三角形重新变换回三维空间，将原来的曲面以交线为边界分为2个曲面，判断其他三角形顶点相对交线的位置，若分别属于不同的曲面，则求取2个曲面的交线。多个地质曲面求交见图8。

2 实验与应用

2.1 算法性能测试对比

在一台型号为Intel i3-2330M、CPU为2.1 GHz、内存为6 GB的计算机上测试算法效率，并与其他2种方法进行对比，见图9。静态结构化数据搜索的速度与分割的包围盒的个数有关，与三角形个数关系不大，时间复杂度为O（2n），保证用于构建OBB混合搜索树的三角形个数随求交三角形个数的增加而平稳增长。该改进算法在数据量大的曲面求交中效率优势明显。

2.2 改进算法在三维地质建模中的应用

由于地质构造的复杂性，地质曲面之间的空间关系非常复杂，因此在建模过程中，地层与地层或断层相交时，通常需要求取交线并离散化，最常用的方法是任意剖面切割法。三维曲面求交方法为基础算法，可以很好地应用到地质自动化建模和分析中。

在三维地质建模中，构建含断层的地质模型是研究热点。根据地质剖面图确定断点、断距和倾角等基本信息，构建初始断层面，得到现有断层面与地层面的求交线，切割上、下2个地层面，重构地层和断层面边界。保证模型和数据的一致性，检查拓扑关系的正确性是核心。含简单断层的地层模型见图10。

复杂的地质構造有许多界面相交，分割线通过一系列线段表示，不同区域之间用多个线段隔开。为自动化模型重构和模型数值分析[27]，需要分离交线分割的不同区域，手动方法繁琐且不可靠。为实现自动分离不同区域，尝试改进前文提出的求交算法以处理简单区域划分的问题。先确定交线类型，再根据交线追踪区域边界，形成内部环。遇到边界就停止追踪，形成多个相邻多边形封闭区域。根据三角形和交线的相对位置，划分所有三角形到其所属的区域内，完成地质区域离散，见图11。

3 结束语

为解决三维地质建模过程和模型分析中的三角曲面求交问题，提出基于优化八叉树的搜索算法和OBB混合搜索树相结合的曲面求交方法。通过结构化搜索快速剔除远离交线的不相交三角形，结合OBB混合层次搜索树，快速求取三角曲面的交线，并将该方法应用到地质建模和分析中，验证算法的可靠性和实用性。采用优化的八叉树快速剔除大量不相交三角形，可提升求交算法的效率，除三角曲面外还可以扩展到其他类型曲面。将该算法应用到地质分析中，解决参数曲面切割地质模型、断层构建和区域离散化问题。将该算法作为基础算法工具应用到三维地质模型的构建和分析中，为快速自动化建模分析提供一种新的方法。

改进的曲面求交方法可以较好地处理地质曲面求交的相关问题，未来可以尝试采用GPU并行计算的方法提升运算效率。在解决地质曲面求交相关应用问题上，工作流程自动化和减少人工干预是未来努力的主要方向。

参考文献：

[1] 邓飞， 王瑞， 王美平， 等. 复杂三维地层建模及快速射线追踪的研究与实现[J]. 大庆石油地质与开发， 2007， 26（1）： 113-118. DOI： 10.3969/j.issn.1000-3754.2007.01.028.

[2] LI N， BAGAS L， LINDSAY M， et al. An irregular triangle mesh buffer analysis method for boundary representation geological object in three-dimension[J]. Earth Science Informatics， 2017， 10（2）： 149-167.

[3] FU Q， WU Z， WANG X， et al. An algorithm for finding intersection between ball B-spline curves[J]. Journal of Computational and Applied Mathematics， 2018， 327： 260-273.

[4] WELLMANN J F， de la VARGA M， MURDIE R E， et al. Uncertainty estimation for a geological model of Sandstone greenstone belt， Western Australia-insights from integrated geological and geophysical inversion in a Bayesian inference framework[J]. Geological Society， London， Special Publications， 2017， 453（1）： 41. DOI： 10.1144/SP453.12.

[5] CALCAGNO P， COURRIOUX G， LOPEZ S， et al. How geological architecture helps 3D modelling[C]// Proceedings of 4th Meeting of European 3D GeoModelling Community. Orleans， 2018.

[6] JIANG Q P， JIE T， YUAN C F. Fast triangle mesh surface intersection algorithm based on uniform grid[J]. Computer Engineering， 2008， 34（21）： 172-174.

[7] BISTACCHI A， MASSIRONI M， dal PIAZ G V， et al. 3D fold and fault reconstruction with an uncertainty model： An example from an Alpine tunnel case study[J]. Computers & Geosciences， 2008， 34（4）： 351-372.

[8] YANG Y， LIU X J， ZHOU B， et al. Research of 3D GIS section for true three-dimensional geo-spatial model[C]// Proceedings of 2010 International Conference on Multimedia Technology. Ningbo： IEEE， 2010： 1-4. DOI： 10.1109/ICMULT.2010.5631185.

[9] 孙黎明， 李青元， 谭海， 等. 基于Netgen的层状地质体四面体网格划分方法[J]. 计算机辅助工程， 2013， 22（3）： 47-52. DOI： 10.3969/j.issn.1006-0871.2013.03.010.

[10] FRANK T， TERTOIS A L， MALLET J L. 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[J]. Computers & Geosciences， 2007， 33（7）： 932-943. DOI： 10.1016/j.cageo.2006.11.014.

[11] LINDENBECK C H， EBERT H D， ULMER H， et al. TRICUT： A program to clip triangle meshes using rapid and triangle libraries and visualization toolkit[J]. Computers & Geosciences， 2002， 28（7）： 841-850. DOI： 10.1016/S0098-3004（01）00110-8.

[12] FABRI A， PION S. CGAL-computational geometry algorithms library[C]// Proceedings of 17th ACM SIGSPATIAL International Symposium on Advances in Geographic Information Systems. Washington， 2009. DOI： 10.1145/1653771.1653865.

[13] 黃松柏， 徐华. 基于动态OBB层次结构的曲面相交算法[J]. 计算机应用研究， 2011， 28（8）： 3181-3184. DOI： 10.3969/j.issn.1001-3695.2011.08.106.

[14] ELSHEIKH A H， ELSHEIKH M. A reliable triangular mesh intersection algorithm and its application in geological modelling[J]. Engineering with Computers， 2014， 30（1）： 143-157. DOI： 10.1007/s00366-012-0297-3.

[15] SAMET H. Foundations of multidimensional and metric data structures[M]. San Francisco： Morgan Kaufmann， 2006.

[16] DASSI F， PEROTTO S， FORMAGGIA L， et al. Efficient geometric reconstruction of complex geological structures[J]. Mathematics and Computers in Simulation， 2014， 106： 163-184. DOI： 10.1016/j.matcom.2014.01.005.

[17] GOTTSCHALK S， LIN M C， MANOCHA D. OBB tree： A hierarchical structure for rapid interference detection[C]// Proceedings of 23rd Annual Conference on Computer Graphics and Interactive Techniques. New York， 1996： 171-180. DOI： 10.1145/237170.237244.

[18] FANG Z G， XU J， JIANG J X， et al. Efficient collision detection using bounding volume hierarchies of OBB-AABBs and its application[C]//Proceedings of 2010 International Conference on Computer Design and Applications. Qinhuangdao： IEEE， 2010. DOI： 10.1109/ICCDA.2010.5541315.

[19] CHANG J W， KIM M S. Efficient triangle-triangle intersection test for OBB-based collision detection[J]. Computers & Graphics， 2009， 33（3）： 235-240. DOI： 10.1016/j.cag.2009.03.009.

[20] VAN DEN BERGEN G. Efficient collision detection of complex deformable models using AABB trees[J]. Journal of Graphics Tools， 2012，2（4）： 37-41. DOI： 10.1080/10867651.1997.10487480.

[21] 鄒益胜， 丁国富， 何邕， 等. 空间三角形快速相交检测算法[J]. 计算机应用研究， 2008， 25（10）： 2906-2910. DOI： 10.3969/j.issn.1001-3695.2008.10.007.

[22] MLLER T. A fast triangle-triangle intersection test[J]. Journal of Graphics Tools， 1997， 2（2）： 25-30. DOI： 10.1080/10867651.1997.10487472.

[23] MCLAURIN D， MARCUM D， REMOTIGUE M， et al. Repairing unstructured triangular mesh intersections[J]. International Journal for Numerical Methods in Engineering， 2013， 93（3）： 266-275. DOI： 10.1002/nme.4385.

[24] HOPPE H， DEROSE T， DUCHAMP T， et al. Mesh optimization[C]// Proceedings of 20th Annual Conference on Computer Graphics and Interactive Techniques. New York， 1993. DOI： 10.1145/166117.166119.

[25] ATTENE M， FALCIDIENO B. Remesh： An interactive environment to edit and repair triangle meshes[C]// Proceedings of IEEE International Conference on Shape Modeling and Applications. Matsushima： IEEE， 2006. DOI： 10.1109/SMI.2006.29.

[26] WANG D， HASSAN O， MORGAN K， et al. Enhanced remeshing from STL files with applications to surface grid generation[J]. Communications in Numerical Methods in Engineering， 2007， 23（2）： 227-239. DOI： 10.1002/cnm.894.

[27] 季顺迎， 赵金凤， 狄少丞， 等. 面向环境力学的离散元分析软件研发和工程应用[J]. 计算机辅助工程， 2014， 23（1）： 69-75. DOI： 10.13340/j.cae.2014.01.014.

（编辑 付宇靓）