矿山三维模型无缝集成方法与研究

李　江, 刘修国

(1.中国地质大学(武汉)信息工程学院,湖北 武汉　430074; 2.湖北省国土资源厅 信息中心,湖北 武汉　430071)

0　引言

近年来,随着数字矿山技术(Digital Mine,DM)的蓬勃发展,三维地质建模理论在采矿领域得到了广泛应用,尤其在地下空间三维建模方面,三维模型的无缝集成更是三维GIS、计算机可视化与城市、矿山及岩土工程等领域研究的交叉热点。数字化三维建模技术及采矿软件以精细、可靠的空间三维集成模型为采矿设计、智能化矿山爆破过程分析、生产调度、虚拟仿真、决策支持等方面矿山生产提供了强有力的支撑平台。目前,针对矿山地下空间地质构造、矿体、地层、矿山巷道、地下工程开挖体等对象进行的三维建模研究基本都是独立建模,各实体的几何模型完全相互独立,实体模型的边界交界处势必形成边界缝隙,而多层界面结构的缝隙也必将影响应用分析的精确度,难以满足矿山应用的要求。因此,开展面向多层次应用需求的多分辨率无缝集成建模,为地下空间的认知活动提供合理、科学、高效的多分辨率无缝集成三维模型,对矿山的可持续发展具有重大意义。

矿山三维模型的建立集成了勘探地质学、数学地质、地球物理、矿山测量、地理信息系统、图形图像学等相关知识,其起源于因矿山管理和采矿的需要,由计算机辅助设计技术发展而来。国外在该领域的研究起步较早,Carlson Eric于20世纪80年代末从地质学角度提出了地下空间三维结构模型,并由此提出了单纯复形模型(Simplicial Complex),构建了反映地质矿体立体概念的地下空间结构理论[1]。三维地学模型由Raper J.F于1989年提出[2]。Fritsch于1990年提出三维GIS数据结构[3]。90年代初,Joe创立了基于三维点集的局部变换构建三维Delaunay三角形算法[4],即基于空球准则的伪局部优化法(Pseudo loeally optima)。Molenar以点、弧、边、面四元素为基础构建了形式化三维数据模型——三维矢量数据模型(3D Formal Data Structure)[5]。Houlding S .W[6-7]提出了规则三维格网(Regular 3D grid)、非规则块(Irregular block)、断面(Sectional)和体(Volume)的数据结构。Vietor和Pilout M在基于Delaunay四面体的三维矢量数据模型(TEN)方面进行了相应研究[8-9]。Lattuada and Raper在1998年对三维GIS环境下的Delaunay三角形在地质领域的应用做了相关研究工作,证明了四面体网格具有较高精度、方便可视化并易于实现关系查询的优点[10]。中国在该领域的研究近年来同样得到了发展和突破,获得了较多实用成果。李清泉、孙敏等提出了以TIN与CSG集成建模作为城市3D GIS和3D CM(City Model)建模的主要方式[11-12]。王彦兵[13]针对城市空间(地上空间、地表空间和地下空间)实体数据来源及分布特征,提出了TIN与广义三棱柱(GTP)集成建模方法。在矿山建模方面,龚健雅等[14]提出了一种矢栅集成的建模方式,将空间对象定义为0-3D及复杂对象,将数字表面模型、数字立体模型、像元和体元、柱状实体等引入到模型中。程朋根[15]等以地勘工程为研究背景,发展了以地勘工程为应用环境的矢栅集成数据模型。侯恩科[16]在其博士论文中也提出了涉及多种空间对象类的三维地质模拟的集成数据模型。

综上所述,目前针对矿山地下空间实体的独立建模研究较多,并形成了较为成熟的理论和方法,但针对地下空间实体间相关关联和影响的集成化建模研究较少,尤其是对复杂地质形体三维模型构建技术的相关理论和技术很不完善,不能满足实际问题的需要。本文以四面体网格的数据结构和生成算法为基础,针对地下空间不同实体的不同结构特征和多种建模数据源,研究实现了不同数据源之间的融合机制,建立了集成式、统一结构的空间数据模型,研究了三维模型无缝集成中的关键算法,提出了矿山地下空间实体无缝集成建模的体系结构和建模方法,实现了矿山空间实体完整、连续、无缝的模型化表达。

1　无缝集成建模关键算法及研究

1.1　Delaunay四面体网格构造算法

地质三维空间数据模型按单元维数可分为面元模型(Facial Model)和体元模型(Volumetric Model )两大类,按建模方式可分为单一建模模型(single modeling)、混合建模(compound modeling)和集成建模(integrated modeling)三大类[17]。

四面体网格(TEtrahedral Network,TEN)由三维Delaunay三角化基础上提出,是不规则三角网模型(TIN)在三维空间上的扩展,其实质是利用紧密排列但不重叠的不规则四面体作为基本体元来完成对空间实体的表达[4,8-9,18]。四面体模型以空间离散点为顶点,离散点集利用互不相交的直线连接形成三角面片,所有互不穿越的三角面片就构成了四面体网格,每个四面体内不包含点集中的任意一点。四面体模型是面最少的体元结构,表达精度高,可通过插值运算实现实体的内部表达,在布尔计算和体积面积等属性计算方面效率较高,便于对空间实体实现三维分析和显示,能较好地应用于地质矿山领域,实现对复杂地质体的描述。本文采用四面体网格为数据模型,实现三维模型的无缝集成。

四面体网格的生成方法中以Delaunay法的使用最为广泛,其中,局部变换法和Bowyer-Waster算法是Delaunay三角剖分的两个经典算法,随着三角剖分研究的深入,越来越多的计算方法开始出现。具有代表性的如:多种Incremental[19-20]算法,Gift-Wrappingl算法,Quick Hull[21]算法和Divide-and-Conquer算法等。现将Divide-and-Conquer(分治-合并)法介绍如下。

首先将数据点集(Point Set,PS)按升序排序,使Pi(Xi,Yi,Zi)

① 数据点集P以X坐标为主,Y、Z坐标为辅排序;② 将P分为近似相等的Pi和Pr;③ 判断Pi/Pr中是否包含4-7个点。如果是,合并Pi和Pr的两个四面体网格,执行LOP算法优化四面体网格,否则,调用Lee(Pi)或 Lee(Pr);④ 递归上述过程,直至四面体网格形成。

上述LOP即局部优化算法(Local Optimization Procedure)采用是Lawson提出的换边、换面算法。该算法在加入点后,对于共边或共面的三角形四面体进行空外接球或空外接圆测试,如外接球或外接圆内包含其他点,则进行换边或换面的优化操作,LOP是运用四面体网格性质对有两个公共面的四面体组成的六面体进行判断,如果四面体的外接球面中包含有第五个顶点,则将这个六面体的公共面相互交换(如图1所示)。

图1　四面体优化示意图

1.2　约束四面体网格(CD-TEN)生成算法

矿山空间对象普遍具有边界复杂、空间特征约束较多等特点,除外部边界外,可能存在的特征约束以两种方式存在:线约束和面约束。线约束是空间一系列特征线形成的约束,如山脊与谷底线、道路与河流等;面约束是空间分界面形成的约束,具有不同的属性,如地质分层、褶皱、断层,采空区、地下工程体外围轮廓等。这些约束的存在增加了用四面体网格模型构建建立矿山空间实体的复杂度,因此,针对约束四面体网格(Constrained Delaunay TEtrahedral Network,CD-TEN)生成算法研究也越来越得以关注。

实现四面体网格的约束化是在普通四面体网格的基础上分三部分进行:恢复约束边,恢复约束面,清除多余的四面体网格。现以四面体网格中的Flip局部变换为例介绍约束边的恢复过程(图2)。

设ViVJ是四面体中待恢复的边,管道P(ViVJ)穿过所有与ViVJ将相交的三角面(阴影部分为相交三角面)。

管道P(ViVJ)形成的关键是快速搜索到所有被约束边穿越三角面。其过程为:以Vi为起点向VJ延伸,穿越过的第一个四面体面为F(该四面体句柄),将F指向面赋予P(ViVJ)后,对于F的邻接四面体F=F.sym,重复以下步骤,最终形成P(ViVJ):

① 如VJ在F内,则P(ViVJ)完成;② 在F的其他三个面中寻找与VJ相交的面,将其赋给F并添加P(ViVJ),取得F的邻接四面体F=F.sym,返回步骤①。

在以上步骤中,如果在查找F的过程中,没有三角面与VJ相交,则该约束边与四面体相交,此时应采用边交换方法替换该约束边并继续以上步骤,如果交换方法不成功,则采用插入点方式完成[22]。

图2　四面体约束边恢复过程

1.3　空间插值算法

建立矿山实体三维模型时,数据的来源一般分为地震资料数据、直接测量数据和钻孔数据三个方面,其中钻探过程是地质勘探过程中获取地下矿体、岩层、地质构造等信息的主要手段,过程中获取的钻孔数据是建立三维矿山模型的主要数据来源,钻孔数据具有不规则性、分散性、空间性等特点,其质量的好坏将影响到建模质量的好坏。

钻孔原始资料信息只局限于描述一定范围内的地质属性值,如果不对数据进行预处理和补充,在其基础上建立的空间实体模型将难有较好的精度保证。因此,在矿山实体建模过程中,利用有限的钻孔采样点数据,采用一定的计算方法,对未知区域的特征数据进行推理和估计,提高整体数据密度,就是空间插值的概念,即通过已知的空间数据按照一定规则估算或预测未知空间数据值[23]。空间插值过程也是一个将非规则分布的地下矿山空间数据转变为规则分布,从而构建出一个连续的地下空间数据曲面的过程。在矿山应用领域,成熟的插值方法包括最近邻点插值法、克里金法(Kriging)和距离幂次反比法,下面以克里金法为例介绍空间插值的过程。

克里金插值方法是在地质统计学方法(Geostatistics)基础上发展的空间插值计算方法。该算法的特点是在插值过程中能够根据特定的优化准则动态决定变量的数值,因此在插值计算中能够有效解决其它插值方法中误差难以分析的问题,并能逐点对误差进行理论估计,避免了“边界效应”的产生[24]。该算法与传统的加权插值算法类似,不同之处在于权重的赋予方面,距离加权反比插法仅考虑了插值点与采样点的距离,而克里金法除距离因素外,通过结构分析和变差函数,同时将采样点的空间分布信息、采样点与插值点的空间方位关系、采样点的总体趋势变化等相关因素统一纳入考虑范围。克里金插值方法范例如下:

设对于插值点0,其邻域内共有n个有效的离散采样点对0产生影响,那么插值点的估算值Z*(X0)如式(1)所示:

(1)

式中:Z(Xi)为已知采样点Xi(i=1,2，…,n)的属性值;λi为加权因子,表示各产生影响的已知点的属性值Z(Xi)对插入点的估计值Z*(Xi)的影响大小。该方法需满足以下两个前提,一是插值点的实际值与估算值的误差期望为0,如式(2)所示:

E[Z(X0)-Z*(Xi)]=0

(2)

二是估计方差最小,如式(3)所示:

(3)

对于无偏条件,应该满足各权系数的和等于1,如式(4):

(4)

(5)

因此,由式(5)可得式(6):

(6)

克里金插值法有很多分支,如普通克里金法、范克里金法、协同克里金法和指示克里金法等。不同的克里金法适用于不同条件,可以满足不同的研究需要。

1.4　曲面求交运算

矿山空间领域存在大量层面、层状现象,当出现断层错断岩层、地层不整合、地层尖灭等情形时普遍会遇到曲面间的求交问题,同时,在进行复杂地质构造的三维无缝建模时,针对生成一致且封闭块体模型的要求,褶皱、断层以及地层界面相互之间的切割和调整处理同样涉及曲面求交,因此,曲面求交(Surface-Surface Intersection,SSI)运算是实体建模的关键,也是三维模型无缝集成过程中的关键算法。曲面求交运算类型共分四类:隐式曲面—隐式曲面求交,隐式曲面—参数曲面求交,参数曲面—参数曲面求交,多边形曲面—多边形曲面求交。下面以曲面求交四个步骤(曲面求交判断、交线跟踪、局部区域重构、优化处理)中曲面求交判断为例对该算法作简要介绍。

曲面求交判断又称作碰撞检测,包括空间分解法(Space Decomposition)和层次包围盒法(Hierarchical Bounding Volumes)。空间分解法是将虚拟空间划分为相等体积的独立单元格,对占据了同一单元格或相邻单元格的几何对象进行相交测试。层次包围盒法是采用体积略大且几何特性简单的包围盒近似描述复杂的几何对象,从而仅对包围盒重叠的对象进行进一步的相交测试,同时通过构造树状层次结构逼近对象几何模型,直至基本完全获得对象的几何特性。比较典型的曲面求交算法包括轴向方向包围盒法(Axis-Aligend Bounding Boxes,AABB)和方向包围盒法(Oriented Bounding Boxes,OBB)。AABB是最简单的包围盒,但其紧密性较差,OBB紧密性最好,但其计算和相交测试的算法较为复杂。

OBB树是利用OBB 建立的树结构,它的基本描述参数为:原点ori (orix,oriy,oriz)、方向dir (dirx,diry,dirz)和大小size ( sizex,sizey,sizez)。其建立过程首先建立整个曲面的OBB,然后将其分解为2 个节点OBB,并将这2 个节点的孩子节点重复以上步骤,直至其最小单元无法继续分割的OBB,最终形成一个仅包含一个三角形的OBB,如图3所示。

图3　OBB包围盒示意图

2　矿山三维模型无缝集成方法与过程

对矿山地下空间实体进行有效、完整的统一描述和无缝集成建模是目前三维建模领域的难点所在,单一的建模方式无法完成对地下复杂实体的空间查询和可视化分析需求,体元模型能够描述地质体的内部品位变化,却无法精确描述地质体边界、矿体与断层褶皱、矿体与矿井巷道等空间关系,面元构模方式能够精确描述矿体或地层的层间边界轮廓,却无法针对地质体内部的属性变化进行精确表达,混合模型存在的数据存储量大、计算复杂等缺点也同样不适用于地质体三维模型的集成。因此,在保证空间地下实体三维模型无缝和拓扑一致的前提下,利用某种合适的关联集成不同模型之间的公共界面,构建一种建立在多模型基础上的新型空间数据模型是建立复杂地质体模型的唯一途径。

2.1　三维模型无缝集成原理

矿山地下三维空间对象涵盖地质体和井巷开挖工程等实体,针对地层、矿体、断层、褶皱、巷道、竖井等实体进行三维无缝建模不仅能实现地下空间的数字化、可视化表现,同时也能对矿山资源评价、储量计算、矿山生产的规划和过程模拟等提供有力支撑。

三维模型无缝集成的原理是基于模型关联和数据关联集成机制,采用一体化的三维空间数据模型和空间数据结构,将地下三维空间对象相互集成。针对地质体模型间的无缝集成原理为:采用CD-TIN模型对地质分界面、断层、褶皱等进行构模,采用TEN模型对地层、矿体等地质体进行构建,利用不同实体间的相交检测、实体切割拆分、单元重构、约束模型重建等方式,完成地质体模型之间的无缝集成。针对地质体模型与地下开挖工程的无缝集成原理是:首先利用GTP模型完成地下开挖工程建模,然后采用模型转换算法,实现GTP模型向TEN模型的转换,将一个GTP单元剖分为三个四面体模型,最后通用体-体模型间相交检测、实体切割拆分、单元重构、约束模型重建方式,完成地质体模型与地下开挖工程的无缝集成。集成原理如图4。

图4　矿山三维模型无缝集成原理图

2.2　断层与地层的无缝构模过程

断层与地层无缝构模的过程是:利用地层数据建立初始四面体网格模型,然后将断层与地层拟合后的面模型进行曲面相交计算,经过曲面重构和局部优化,将产生的集合三角面作为约束条件,最后生成断层和地层的一体化无缝三维模型。断层与地层统一无缝构模流程图如图5所示。建模过程可分为七个步骤:

(1)数据处理。针对不同类型和不同精度的原始数据按照模型预先设计的数据库格式进行整理,将整理后的各类数据如钻孔断点数据、地质剖面数据、断层表格数据、地层分层数据等进行分类保存。

(2)数据加工。对原始采样点数据进行加工处理。针对地层数据,利用加密、插值计算对分布离散的钻孔数据进行补充,使之满足地层层面拟合需要。对断层数据,需对稀疏的断点数据加密、插值,同时根据地质工作经验对数据补充修正,对地质剖面图反映的各地层与断层的交切关系以及各个断点的准确坐标进行完善补充。

(3)生成初始三维地层模型。利用完善后的地层层面钻孔数据,采用四面体生成算法,完成三维地层实体模型构建。

(4)拟合断层面。利用断层参数数据、断点数据、剖面数据等断层数据,采用约束Delaunay三角剖分算法拟合合理的断层面网格。如断层之间存在切错关系,需首先按多条断层的切错关系进行优先排列,然后再进行切割处理。

(5)拟合地层层面。该步骤以加密插值后的钻孔数据为基础对每个地层层面分别进行拟合。生成地层层面三角网模型时需考虑断层线的约束作用。

(6)断层面与地层面曲面相交计算。利用曲面求交运算,对断层面与地层面进行切割、调整、约束处理。当地层层面与断层面完成切割调整后,还要根据断层之间的断距对交线附近的地层层位高程进行调整,并将调整后的交线约束到断层面上。

(7)地层与断层无缝建模。利用上一步骤生成的三角网结合作为约束条件,针对初始的地层三维模型,采用约束四面体网格的生成算法,而生成无缝且一致的断层与地层一体化模型,最终实现地层与断层的无缝建模。

图5　断层与地层统一无缝构模基本流程图

2.3　矿体与地下开挖工程间无缝构模方法

矿山地下开挖工程是矿山空间领域内根据一定几何、属性、拓扑结构进行一致表达的特定空间实体,包括斜井、竖井、巷道等。矿山地下开挖工程在地下空间占有一定的空间和范围,具有一定的形态和属性,具备与地质体一样的三维特性并与自然地质对象紧密联系。

矿体与矿山开挖工程无缝构模的基本思想是:分别利用TEN和GTP模型建立矿体与工程体的三维模型,然后,利用转换算法将工程体GTP模型转换为TEN模型,将转换后的工程体模型与矿体三维模型(TEN)进行体—体模型布尔运算,最终实现矿体与矿山地下开挖工程的三维模型无缝集成,图6以矿山巷道和矿体的无缝建模为例,展示其基本建模流程。

(1)数据处理。即钻孔原始数据和矿山巷道勘测数据的预处理过程。包括针对单工程矿段截取过程所需的数据补充和修正,针对矿体初始三维模型建立的钻孔数据加密插值以及针对矿山巷道横断面的辅助数据增加等;

(2)矿体初始四面体网格建立。将加密插值后的钻孔数据根据矿体品位划分为多个离散数据段,建立矿体的初始三维模型;

(3)勘察剖面图的生成。包括单工程矿段截取和矿体圈定两个过程,最终形成包含多个矿体边界的多边形轮廓线;

(4)矿体三维模型的建立。以勘察剖面图为约束条件,在矿体的初始三维模型的基础上,采用约束四面体网格的生成算法,最终形成矿体三维模型;

(5)矿山巷道GTP模型建立。利用矿山巷道勘测数据,针对横断面形状固定和变化的两种情况,生成巷道GTP三维模型;

(6)巷道GTP模型向四面体模型转换。利用前述的转化算法,实现GTP模型转换为UC-TEN模型的转换,此时,矿体和巷道的三维模型都以四面体网格形式构建;

(7)矿体和巷道三维模型集成。利用四面体之间的布尔运算,按相交检测、单元拆分、单元重构的体相交原则,最终实现两者之间的无缝构模。

图6　矿体与矿山巷道无缝构模基本流程图

3　结语

矿山地下空间的无缝集成三维建模是数字矿山和三维GIS领域的难点。本文针对复杂地质体和矿山地下开挖工程的建模过程,提出一种基于四面体网格模型的无缝集成建模方法,其总体思路为:利用基础数据经插值加密后构建基础三维模型,然后利用多重约束和曲面求交运算等方法重构各种界面,最后通过块体相交算法重构三维实体模型。该方法已经在部分矿山储量计算过程中实现了复杂地质体的无缝建模与可视化,并提供了相应的计算模型,下一步将完善无缝集成模型的空间分析、数据挖掘和地学统计等相应功能,使之真正成为数字化矿山的核心技术和平台。

参考文献：

[1]Erie C.Three-dimensional conceptual modeling of subsurface struetures[C].Baltimore ed.ASPRS,ACSM,Annual Convention,1987,1987(4):188-200.

[2]Raper J.F. The-dimensional geo-seientific mapping and modeling systems : aconceptualdesign[M].London: Taylor & franeis,1989:11-19.

[3]Fritseh D. Towards three-dimensional data struetures in Geographic Information Systems[A].EGIS Foundation,Utrecht,1990:335-345.

[4]JoeB.Constrution of three-dimensional Delaunay triangulatios using local transformations[J].Computer Aided Geometric Design,1991,8(2):123-142.

[5]Molenaar M. A topology for 3D vector maps[J].ITC Journal，1992(1):25-33.

[6]Houlding.S.W.3D Geoscience Modeling: Computer Techniquesfor Geological Characterization [M].New York and Heidelburg:Springon-verlag,1994:121-125.

[7]Houlding.S.W. Practical geostatistics,modeling and spatial analysis [M]. NewYork and Heidelburg: Springer-Verlag，2000.

[8]Vietor. J.D.,Alan P. Delaunay tetrahedral Data Modeling for 3D GIS Application[A].Proc .GIS/LIS.1993,93:671-678.

[9]Pilout M.,Tempfli K.,Molenaar M. A Tetrahedron-based 3D Vector Data Model for Geoinformation[C]//Netherlands Geodetic Commision,Publications on Geodesy.1994(40):129-140.

[10]Roberte Lattuada,Jonathan Raper. Application of 3D Delaunay triangulation algorithms in geoscientific modelling[C]//The Third International Conference/Workshop on Integrating GIS and Environmental Modeling CD-ROM,1998：1-9.

[11]李清泉,杨必胜,史文中,等.三维空间数据的获取、建模与可视化[M].武汉:武汉大学出版社,2003.

[12]孙敏.三维城市模型的数据模型研究[D].长沙:中南工业大学,2000.

[13]王彦兵.基于TIN耦合的城市地上下三维空间无缝集成建模研究[D].北京:中国矿业大学,2005.

[14]龚健雅,夏宗国.矢量与栅格集成的三维数据模型[J].武汉测绘科技大学学报,1997,22(l):7-15.

[15]程朋根,龚健雅.地勘工程三维空间数据模型及其数据结构设计[J].测绘学报,2001,30(1):74-81.

[16]侯恩科.三维地质模拟的若干关键问题研究[D].北京：中国矿业大学,2002.

[17]吴立新,车德福,郭甲腾.面向地上下无缝集成建模的新一代三维地理信息系统[J].测绘工程,2006,15(5):1-6.

[18]Pilouk M.Integrated modelling for 3D GIS [D].The Netherlands: ITC，1996.

[19]Barry Joe.Construction of Three-Dimensional Triangulations using Local Transformations[J].Computer Aided Geometric Design,1991,8:123-142.

[20]Herbert Edelsbrunner and Nimish Shah. Incremental topological flipping works for regular triangulations[J].Algorithmica,1996,15:223-241.

[21]Qhull,The Geometry Center,University of Minneapolis,Brad Barber,bradb@geom.umn.edu,http://www.geom.umn.edu/locate/qhull.

[22]Steven J.Owen,Constrained Triangulation:Application to Hex-Domaint Mesh Generation[C]//Proceedings,8th International Meshing Roundtable,South Lake Tahoe:CA,1999:31-41.

[23]李翠平,李仲学,胡乃联.面向地矿工程体系化3种空间插值方法的应用 [J]. 金属矿山,2004(4):56-58,77.

[24]Jorunel A. G.,Huijbregts CH.矿业地质统计学[M].侯景儒,黄競先,译. 北京:冶金工业出版社,1982.