一种新的复杂地质体采矿工程剖面图自动生成方法

谭正华 ，王李管，熊书敏，刘任任

(1.湘潭大学 湖南省高等学校智能制造重点实验室，湖南 湘潭，411105；2.中南大学 数字矿山研究中心，湖南 长沙，410083)

地质体工程剖面图通常是沿勘探线方向切绘，反映出该剖面上的矿(岩)体界线、地质体岩性，工程地质特征、水文特征、断层位置和构造形态等，它是分析地质构造、储量估算和指导采矿设计的重要图件之一[1]，是工程师、设计管理人员等交流的共同语言。地质剖面图的生成要依赖于特定的地质体模型，由于地质体数据量大、结构复杂(多含夹石、断层)，且当矿体内含有孔隙或岩浆侵入另一种矿体，切割后形成的地质区域通常含有孔洞或孤岛，不同区域的部分边界相互粘连。同时，为了表示不同的地质构造，需要对不同地质体内部区域进行识别，并进行岩性图案的填充。国外学者较早地提出了一些三维实体切割、开挖算法，主要算法有Shostko等[2]提出并实现的算法、gOcad项目组实现的切割算法[3]和三角网切割算法TRICUT[4]，主要思想是将复杂地质模型的剖面生成问题归结为空间三角网的相互切割问题，基本操作是三角形与三角形求交和平面与三角形求交，存在大量的几何运算，时间复杂度较高。另外，国内学者在实体切割和地质体剖面图生成问题上提出了一些改进算法：陈国良等[5]采用一种基于OBB树的三角网切割算法来实现空间两个三角网格的切割；陈建宏等[6]提出了基于钻孔数据的剖面图生成算法；朱莹等[7−8]提出采用GIS组件技术实现地质体剖面图的生成等。本文作者提出一种新的基于平面隐函数的实体切割思想和区域自动识别技术的地质体剖面图自动生成方法，将三角形和切割面之间的几何运算转换为数学运算，同时通过树结构形象表达地质体轮廓线形成的复杂区域(内含孔、岛)，最后通过对复杂区域的岩性图案填充，从而实现地质体剖面图自动生成。实例验证了该算法有效性和可行性。

1 基本原理及流程

由于地质体的形态复杂，直接基于原始地质模型的剖切较为复杂。因此，本文作者提出基于隐函数平面切割算法，加快几何运算，生成离散的交线；采用基于KD树的空间索引方法，确定交线的邻接关系，生成地质体轮廓线；采用树结构形象表达地质体轮廓线形成的复杂区域（简称“区域树”）；为了表达地质体岩性信息，对剖面图地质体区域按照标准岩性图例进行图案填充，完成剖面图的自动绘制。算法流程见图1。

图1 地质体剖切算法流程Fig.1 Process of 3D entities cutting

2 三维地质结构模型的组织

地质体剖切算法依赖特定的数据模型。表达地质体基本数据模型有表征实体外表轮廓的表面模型、基于体元的块段模型以及混合模型[9]。其中，表面模型侧重于三维空间实体的表面表示，能精确表达地质体的形态。这些被模拟的实体可以是封闭的，如矿体、硐室和采场等；也可以是非封闭的，如地层界面、断层和地表等。

数字采矿系统多采用不规则三角格网(TIN)来模拟各种复杂地质体边界及人工工程的表面，本文研究的切剖对象为TIN表达的地质体模型，如图2所示。

图2 基于TIN表达的地质体模型Fig.2 Geological model expressed by TIN

3 算法的关键过程

3.1 隐函数的平面切割算法

3.1.1 隐函数的性质

用平面切割多面体形成一系列的轮廓线，是将三维问题转化为二维问题关键的一步，速度和正确性是首先要关注的问题。地质剖面图一般是沿勘探线上的平面切割，切割面上的所有三角面片，它们的平面方程是相同的，这样的平面可以采用一个方程来表示，如：

根据这一重要特点，采用隐函数实现平面切割算法，省去了大量几何运算，提高了运算速度。

隐函数具有以下3个重要性质。

性质1：可以描述一些基本几何形态，包括平面、球、圆柱、圆锥、椭球体和二次曲面；

性质2： 将三维欧式空间分为3个不同区域，即在隐函数里、上和外3个区域，这些区域分别被定义为F(x,y,z)＜0,F(x,y,z)=0,F(x,y,z)＞0。

性质3： 可以将空间位置化为一个数值，如将空间的坐标(xi,yi,zi)转变为一个数值ti。

利用隐函数的性质 3，并结合等值面(线)算法(如Marching cubes等算法)，很容易将切割线与被切实体之间的交线求出。

3.1.2 基于平面隐函数的实体切割

基本思想是：通过隐函数为被切实体所有的顶点坐标产生一个数值，然后求F(x,y,z)=0的等值线(根据隐函数性质2)，即为所求的交线。基本步骤如下：

步骤1: 通过隐函数(1)为被切实体的所有顶点产生数值；

步骤2: 求F(x,y,z)=0的等值线，即为交线，形成一个或多个多边形区域。

对于每一个三角形面片，如果它的每一个三角形都为正或为负，则该三角形不会被切割面切割，也不会形成交线，然而，如果这些数值既有正也有负，则该三角形被切割，必然有交线通过该三角形，交点位于三角形的边上，通过沿边线线性插值计算交点坐标。三角形被切割可以分为如图3所示8种情况。对于每个三角形，根据其3个顶点的数值找到对应8种情况的一种情况，找到交点所在的边。如果其中一个交点所在边的2个顶点对应的F(x,y,z)数值为v1，v2。则交点的数值为：

则交点坐标为：

图3 切剖三角形的8种情形Fig.3 Eight cases of triangular surface cut

两交点连线构成该三角形和平面的切割线。

3.2 基于KD树的轮廓线生成算法

如果一个平面切割一个封闭的实体，如存在交线，则所有交线合在一起将形成1个或多个封闭的交线圈(多边形)。三角形和平面的交线是一系列空间离散的直线段，不是闭合的轮廓线，必须将前后相邻的线段串联为1个或多个多边形区域，以便后续处理。

由n条孤立的直线段形成闭合的多段线，最简单的方式是：取出1条直线段的1个端点与其他剩下的直线段端点比较，如果相等，则将被找到的直线和前面的直线段连接，然后取出被找到直线的另一个端点，继续查找，直到形成一闭合多边形或搜索完毕。该方法速度比较慢，最坏情况下时间复杂度为O(n2)。

本文基于KD树查找直线段的邻接关系，快速实现闭合多段线的建立。

3.2.1 KD树简介

KD树[10−11]是二叉搜索树的扩展，K表示空间维数，故又称K维搜索树，特别适合空间点状目标的索引，其示意图见图4。在每个内部节点中，用一个K−1维与坐标轴平行的超平面(如二维空间中的线)将节点所表示的K维空间分成2个部分，存储在子树中的点大约一半落入一侧，而另一半落入另一侧。当一个节点中的点数少于给定的最大点数时，划分结束。这些超平面在K个可能的方向上交替出现，每个超平面中至少包括一个数据点，这就保证了KD树不存在“无点空间”的浪费。

图4 KD树示意图(K=2)Fig.4 KD tree diagram

3.2.2 KD树构建及索引

为了构建平衡树，本文给出交点数为n时KD树构建算法。

步骤 1：初始化一个新结点 TopNode，将所有点作为其点数据；

步骤 2：判断点数据个数是否小于最小个数，如果小于最小个数，则结束；

步骤 3：沿数据的最大延伸方向对数据排序，采用过中间的位置的点且垂直于最大延伸方向的超平面分割数据，形成2个新左子结点LeftNode和右子结点RightNode；

步骤 4：对LeftNode执行步骤2；

步骤 5：对RightNode执行步骤2。

为了便于查询，KD树记录其覆盖范围最大、最小值，其查询算法和二叉树查找相同，在此不再赘述，其时间复杂度为O(nlog2n)。

3.3 基于树结构的封闭区域组织和索引

基本原理是采用树结构来表达区域的存储和组织方式，以下称为区域树[11]。区域树有如下特点：

(1)一个封闭区域映射一个区域树。

(2)区域树的结点表示区域的内、外边界，其根结点表示该区域的最外围边界，子结点表示区域内部的孔、岛边界。

(3)结点之间的上下层次关系表示边界的包含关系，同层结点表示边界的并列关系。

(4)边界是有方向性的，并对边界方向作如下规定：区域最外围边界的走向为逆时针方向，内部孔洞边界沿顺时针方向，岛边界沿逆时针方向。

某地质体切割轮廓线示意图如图5所示。该地质体轮廓线将平面划分为 3个区域(包含最外的无穷大区域∞)。目标区域及其对应的树如图6所示。

图5 某地质体轮廓线示意图Fig.5 Schematic outline of a geological

图6 区域及其对应的树结构Fig.6 Regions expressed by tree

区域的提取方法有多种[12−14]，其中文献[13]给出了定向极大、极小闭环的概念，通过建立结点−路径网络拓扑结构图，采用内层路径优先查找算法，提取所有定向闭环，最后建立定向极小闭环和定向极大闭环的隶属关系，建立区域树，有效地实现了复杂区域的识别和提取，具有普遍适用性。本文采用文献[13]的算法实现了复杂轮廓线的封闭区域提取和识别，当存在l个极大闭环，k个极小闭环时，该算法时间复杂度为O(l×k)。

区域树的构造过程见参考文献[13]。区域树的构造为地质体剖面图矢量图案填充提供了数据准备。

4 算法分析和应用实例

4.1 算法复杂度分析

该算法时间复杂度主要由3个部分决定：基于平面隐函数的实体切割，基于KD树的闭合轮廓线构造和区域树的构造。假如构成复杂地质体三角面片为m，平面切割后生成n条孤立直线段，所有直线段构成l个极大闭环，k个极小闭环时，则总的时间复杂度为：。

4.2 应用实例

通过 VC+开发工具和可视化工具包实现了该算法。运行环境为Windows 2000操作系统，CPU为Inter Pentium4 3.00 GHz，内存512 MB。当复杂地质体三角面片数为1 589，3 264，5 470，7 412，10 345时，根据勘探线剖切形成的离散直线段、定向闭环和切割三角面片耗时见表1，并做出三角形个数−切割时间曲线，如图7所示。由图7可以看出，该三角形切割算法时间复杂度近似线性。

表1 不同三角面片个数算法耗时Table 1 Consuming time for different number of triangles

图7 三角面片数−时间曲线Fig.7 Time changes with different number of triangles

“DIMINE”数字矿山软件是中南大学数字矿山研究中心开发的一套矿床地质建模、矿山工程测量、数据处理和地下、露天矿床开采设计等功能于一体的集成化、智能化、可视化软件系统。该系统中地质平剖面投影图的生成、纵剖面图的绘制等模块中采用了该算法。该算法可有效地用于复杂地质体剖面图自动生成操作，如图8所示。

图8 地质体剖面图自动生成示例Fig.8 An example of automatic generation of Geological Profile

5 结论

(1)提出平面隐函数的实体切割算法，该算法将三角形和平面间的几何运算转换为数学运算，加快了运算速度，实现了地质实体的快速剖切。

(2)提出基于 KD树封闭轮廓线生成方法，同过对邻近点对的快速索引，实现了封闭轮廓线的快速生成。该方法同样适用于开口实体(如地层、地表)的非闭合轮廓线生成。

(3)采用图论中的树结构形式化表达区域的空间组织结构，为岩性图案填充提供了数据来源。

(4)为进一步提高运算速度，可建立三角形之间的邻接关系，从而在生成离散线段的同时，快速建立线段间的相邻关系，是下一步研究的难点和方向。

[1]曹燮明.采矿手册[M].北京: 冶金工业出版社,1999:172−175.CAO Xie-ming.Mining guide[M].Beijing: Metallurgical Industry Press,1999: 172−175.

[2]Shostko A A,Lohner R,Sandberg W C.Surface triangulation over intersecting geometries[J].International Journal for Numerical Methods in Engineering,1999,44(9): 1359−1376.

[3]Coelho L C G,Gattass M,Fiqueiredo L H.Intersecting and trimming parametric meshes on finite-elements shells[J].International Journal for Numerical Methods in Engineering,2000,47(4): 777−800.

[4]Lindenbecka C H,Ebertb H D,Ulmera H.A program to clip triangle meshes using the rapid and triangle libraries and the visualization toolkit[J].Computers & Geosciences,2002,28(6):841−850.

[5]陈国量,刘修国,商建嘎,等.三维地质结构模型的切割分析技术及方法[J].计算机工程,2007,33(20): 184−186.CHEN Guo-liang,LIU Xiu-guo,SHANG Jian-ga,et al.Incision analysis technology and method for 3D geological structure model[J].Computer Engineering,2007,33(20): 184−186.

[6]陈建宏,周智勇,陈纲,等.基于钻孔数据的勘探线剖面图自动生成方法[J].中南大学学报: 自然科学版,2005,36(3):487−490.CHEN Jian-hong,ZHOU Zhi-yong,CHEN Gang,et al.Automatic formation method of prospecting line profile map based on drill hole database[J].Journal of Central South University: Science and Technology,2005,36(3): 487−490.

[7]朱莹,刘学军,陈锁忠,等.基于GIS的地质剖面图自动绘制软件的研究[J].南京师大学报: 自然科学版,2007,30(4):104−108.ZHU Ying,LIU Xue-jun,CHEN Suo-zhong,et al.A research on automatic generating software of geologic section based on GIS[J].Journal of Nanjing Normal University: Natural Science Edition,2007,30(4): 104−108.

[8]王建芳,包世泰,余应刚,等.基于GIS模板的地质剖面图模型及其实现[J].测绘科学,2008,33(5): 184−186.WANG Jian-fang,BAO Shi-tai,YU Ying-gang,et al.Realization of geological section map model based GIS template[J].Science of Surveying and Mapping,2008,33(5): 184−186.

[9]王家耀.空间信息系统原理[M].北京: 科学出版社,2003:100−105.WANG Jia-yao.Principles of spatial information systems[M].Beijing: Science Press,2003: 100−105.

[10]Bentley J L.Multidimensional binary search trees used for associative searching[J].Commun ACM,1975,18(9): 509−517.

[11]Bentley J L.K-d trees for semi-dynamic point sets[C]//Proc 6th Annual Symposium on Computational Geometry.New York:Association for Computing Machinery,1990: 187−197.

[12]谭正华,王李管,毕林,等.参数曲线集复杂区域的全自动识别算法[J].计算机工程,2010,36(8): 23−26.TAN Zheng-hua,WANG Li-guan,BI Lin,et al.Automatic identification algorithm for complicated regions of parameterized curve set[J].Computer Engineering,2010,36(8):23−26.

[13]宋怀波,路长厚,王富春.一种新的复杂区域孔洞填充算法[J].桂林电子科技大学学报,2006,26(6): 451−454.SONG Huai-bo,LU Chang-hou,WANG Fu-chun.A new hole-filling algorithm for complicated connecting regions[J].Journal of Guilin University of Electronic Technology,2006,26(6): 451−454.