三维矿产资源预测评价中信息综合处理装置及其方法
【专利摘要】本发明公开了一种三维矿产资源预测评价中信息综合处理装置及其方法，其中该方法包括：步骤一，构建三维成矿预测的初始预测体；步骤二，对初始预测体快速提取空间离散点的预测要素；步骤三，对由表面模型表达的预测要素进行栅格化，并将栅格化结果赋值到初始预测体的属性表中；步骤四，对初始预测体的属性表进行筛选与优选，获取预测目标体。该方法用于辅助地质专家解决已知矿山深部及外围找矿工作，满足实现基于地质、地球物理、地球化学等数据的多元成矿信息的信息自动模拟、综合分析与可视化等方面的需求。
【专利说明】三维矿产资源预测评价中信息综合处理装置及其方法
【技术领域】
[0001]本发明涉及数学地质与地理信息系统地图制图领域，特别是涉及一种大比例尺三维成矿预测多元信息综合模型快速构建、信息自动综合与可视化分析的技术及装置。
【背景技术】
[0002]随着地表矿、浅部矿及易识别矿的日益减少，找矿难度日益增大，找矿效果日益降低，进入21世纪80年代以来，重点找寻深部隐伏矿、新类型矿及新领域矿已经成为世界各国所关注的找矿方向，其中使用三维立体模型开展深部隐伏复杂地质体(即，第二找矿空间)的寻找，已经成为许多国家和地区找矿的主要对象，因此大比例尺(通常大于1:50000)的成矿预测的作用更加突出，它已经成为矿床勘查工作的一个重要组成部分(赵鹏大等，1992)。
[0003]目前，国内外还没有开发出一套成熟的专门用于深部矿调、大比例尺三维成矿预测工作的应用软件。大量的三维成矿预测研究首先通过矿山勘探软件完成成矿信息提取与三维可视化分析等工作，其次，在其他软件中完成信息综合分析。国内外矿区勘探三维建模与可视化，地球化学、地球物理数据三维可视化分析等应用领域开发了不少成熟实用的三维地学软件。从软件应用角度分析，国内的三维软件市场主要被国外软件厂商垄断。因此更加需要大比例尺成矿预测研究与应用等实际要求的软件工具。
[0004]综上所述，形成一套集快速建模，高精度三维可视化分析，定位、定量以及定概率预测方法为一体的具有独立自主知识产权的专业三维成矿预测方法，已经成为找矿工作中急需的技术支撑手段。为达成这一目的，需要应用到如下技术:
[0005]I)矿产资源评价中独立目标图层信息综合模型
[0006]在矿产资源评价领域，相互独立的图层是指控矿证据层之间都是独立的。如果一个控矿存在取决于另一个控矿因素，且把这两个都作为证据层，将会产生过高或者过低的成矿概率，将导致最终的预测结果受到影响。因此，在实际应用中，必须使用数学模型对进入成矿预测的图层进行条件独立性检验，剔除非独立图层；目标图层是指从原始变量中挑选与预测矿种有较密切关系的控矿证据层，这一挑选的过程是基于数学模型实现的；信息综合是指在控矿证据图层具有独立且与成矿有利的基础上，将在同一比例尺、同一数学基础的不同空间相同位置的信息叠加在一起，以每一个图层作为一个输入特征的属性，通过一定的数学模型将上述特征属性组合起来用于描述新的特征，即矿产预测靶区。
[0007]2) 二维扫面线算法
[0008]二维空间建立索引结构的方法有很多，例如:AVL树、红黑树、BSP树、均匀条带分害I]、等格网分割、四叉树分割和自适应分块等方法。当数据量很大时，为了能够快速完成查询检索，同时考虑到时间与空间复杂度问题，一般使用网格索引技术。二维扫描线最初是应用于二维图像处理领域。其主要思想是通过分别向待判定像素相邻的4个方向进行某种条件搜索来确定该像素的状态。为了实现快速、准确地高精度三维多元信息综合模型的构建技术，将上述两种算法思想应用到三维空间，解决复杂地质体的高精度立方格模型(即块段模型)快速建模。
[0009]3)三维空间K阶近邻搜索技术
[0010]K-近邻查询是一类空间检索算法，主要解决针对空间离散点的搜索问题，是地理信息系统的基本算法之一。随着三维离散点的获取技术不断完善，对于大量三维离散点的分类、插值和曲面重建等需求不断增加，快速、高效地实现给定条件下的K-近邻查询已经成为当前三维地理信息系统研究的热点问题。K-近邻查询是在二维或三维离散数据集S中查找K个与待查点欧氏距离最近的点。
[0011]在大比例尺成矿预测研究中，存在大量地质、地球化学、地球物理原始点采样数据。针对这些原始数据提取的方法是基于预测单元格中心点的空间插值算法，例如:以单元格的中心为待插值点，使用克立格、反距离加权法等方法获取单元格的属性信息。因此，K阶近邻的快速搜索技术是实现预测要素快速提取的重要步骤之一。
[0012]4)基于证据权重算法与三维可视化技术的隐伏靶区优选技术
[0013]近年来，矿产资源勘查评价中的GIS技术应用的通常做法是采用所谓证据权重法。证据权法是定量分析和综合多源专题地学数据层的方法，生成量化的含矿潜力图。在采用传统方法的矿产勘查工作中，常常需要对比分析地质、地球物理和地球化学等专题图件，圈出值得进一步工作的靶区，证据权法实际上就是传统勘查方法的数字化“版本”。
[0014]证据权重法中的“证据”由一系列勘查数据集(即地学专题图件)组成，“权重值”是根据已知矿床或特定勘查区勘查模型估算出来的(Basab Mukhopadhyay等,2000)。其基本思路是将分别赋以不同权值的专题图数据(证据)在GIS环境下进行空间叠加和综合分析，生成具有不同含矿概率值的网格专题图，这种定量分析结果反映了不同的含矿潜力分区在空间上的分布，从而为进一步的详细勘查工作提供区域上的超前评价结果。
[0015]5 )基于OpenGL的三维可视化技术
[0016]本发明中提及的三维可视化技术实现主要指基于OpenGL技术的三维可视化程序的开发。OpenGL是OpenGraphicsLib的缩写,是一套三维图形处理库,也是该领域的工业标准。计算机三维图形是指将用数据描述的三维空间通过计算转换成二维图像并显示或打印出来的技术。OpenGL被设计成独立于硬件，独立于窗口系统的，在运行各种操作系统的各种计算机上都可用，并能在网络环境下以客户/服务器模式工作，是专业图形处理、科学计算等高端应用领域的标准图形库。
【发明内容】
[0017]本发明的目的在于提供一种三维成矿预测多元信息综合处理装置及其方法，用于辅助地质专家解决已知矿山深部及外围找矿工作，满足实现基于地质、地球物理、地球化学等数据的多元成矿信息的信息自动模拟、综合分析与可视化等方面的需求。
[0018]为了实现上述目的，本发明提供一种三维矿产资源预测评价中信息综合处理装置，其特征在于，包括:
[0019]预测体构建模块，用于构建三维成矿预测的初始预测体；
[0020]空间搜索模块，连接所述预测体构建模块，用于对初始预测体快速提取空间离散点的预测要素；
[0021]变量提取模块，连接所述空间搜索模块，用于对由表面模型表达的预测要素进行栅格化，并将栅格化结果赋值到初始预测体的属性表中；
[0022]靶区优选模块，连接所述预测体构建模块、所述变量提取模块，用于对初始预测体的属性表进行筛选与优选，获取预测目标体。
[0023]所述的装置，其中，所述预测体构建模块进一步包括:
[0024]外包围盒提取模块，用于提取体表面或曲面数据的外包围盒；
[0025]外包围盒网格化模块，用于按照设定的网格粒度将外包围盒分解为六面体；
[0026]面体求交模块，用于根据外包围盒进行三角面片与六面体求交计算；
[0027]六面体分类模块，基于种子填充法和求交计算结果,对构成三维成矿预测体的六面体进行快速筛选，生成初始预测体以及六面体模型的属性表。
[0028]所述的装置，其中，所述空间搜索模块进一步包括:
[0029]第一索引建立模块，用于建立三维离散点的八叉树空间索引；
[0030]第二索引建立模块，用于建立八叉树空间索引的三维R-Tree树索引；
[0031]搜索范围获取模块，用于获取待查点的搜索范围；
[0032]空间插值实现模块，用于根据搜索范围和第一索引建立模块、第二索引建立模块的结果，并应用克立格、反距离加权公式进行空间插值；
[0033]搜索赋值模块，用于将空间插值结果赋值给初始预测体，并保留在初始预测体的属性表中。
[0034]所述的装置，其中，所述变量快速提取模块进一步包括:
[0035]栅格化处理模块，用于栅格化表面模型体或曲面数据；
[0036]结果赋值模块，用于将栅格化结果按照空间位置赋值到初始预测体的属性表中。
[0037]所述的装置，其中，所述靶区优选模块进一步包括:
[0038]成矿概率计算模块，用于以初始预测体的属性表为初始数据，基于证据权重法计算三维预测体中每个六面体的成矿概率；
[0039]优选体模块，用于以初始预测体模型为基础，并结合成矿概率计算结果，提取成矿预测目标体。
[0040]为了实现上述目的，本发明提供一种三维矿产资源预测评价中信息综合处理方法，其特征在于，包括:
[0041]步骤一，构建三维成矿预测的初始预测体；
[0042]步骤二，对初始预测体快速提取空间离散点的预测要素；
[0043]步骤三，对由表面模型表达的预测要素进行栅格化，并将栅格化结果赋值到初始预测体的属性表中；
[0044]步骤四，对初始预测体的属性表进行筛选与优选，获取预测目标体。
[0045]所述的方法，其中，所述步骤一，进一步包括:
[0046]提取体表面或曲面数据的外包围盒；
[0047]按照设定的网格粒度将外包围盒分解为六面体；
[0048]根据外包围盒进行三角面片与六面体求交计算；
[0049]基于种子填充法和求交计算结果，对构成三维成矿预测体的六面体进行快速筛选，生成初始预测体以及六面体模型的属性表。
[0050]所述的方法，其中，所述步骤二，进一步包括:[0051]建立三维离散点的八叉树空间索引；
[0052]建立八叉树空间索引的三维R-Tree树索引；
[0053]获取待查点的搜索范围；
[0054]根据搜索范围和第一索引建立模块、第二索引建立模块的结果，并应用克立格、反距离加权公式进行空间插值；
[0055]将空间插值结果赋值给初始预测体，并保留在初始预测体的属性表中。
[0056]所述的方法，其中，所述步骤三，进一步包括:
[0057]栅格化表面模型体或曲面数据；
[0058]将栅格化结果按照空间位置赋值到初始预测体的属性表中。
[0059]所述的方法，其中，所述步骤四，进一步包括:
[0060]以初始预测体的属性表为初始数据，基于证据权重法计算三维预测体中每个六面体的成矿概率；
[0061]以初始预测体模型为基础，并结合成矿概率计算结果，提取成矿预测目标体。
[0062]本发明提出的方法，适用于大比例尺三维成矿预测中预测要素模拟与多元找矿信息综合分析模型快速构建及其可视化分析等领域，主要用于针对大比例尺矿产勘查过程中获得的各种数据，例如采样点数据、地质地形数据等，利用快速空间插值、曲面自动拟合、栅格模型快速构建等三维计算机技术，实现地质、地球物理、地球化学等多元数据的三维信息自动模拟、综合分析与可视化显示。
【专利附图】
【附图说明】
[0063]图1是三维矿产资源预测多元信息综合处理装置结构图；
[0064]图2是三维矿产资源预测多元信息综合处理方法流程图；
[0065]图3A、3B是R-Tree空间索引模型示意图；
[0066]图4是本发明大比例尺矿产资源评价多元采样点的快速空间插值方法流程图；
[0067]图5A、5B、5C是基于Flood-Fill思想的栅格快速筛选算法的基础——体素之间的相邻关系；
[0068]图6是本发明基于复杂地质体模型与地表曲面模型约束的找矿信息综合分析模型快速构建方法流程图；
[0069]图7是本发明Flood-Fill算法二维空间实现示意图；
[0070]图8是本发明Flood-Fill算法二维空间实现不意图；
[0071]图9是本发明中变量提取模块的示意图；
[0072]图10是本发明中变量提取模块的示意图。
【具体实施方式】
[0073]以下结合附图和具体实施例对本发明进行详细描述，但不作为对本发明的限定。
[0074]如图1所示，是三维矿产资源预测多元信息综合处理装置结构图。该装置100包括:预测体构建模块10、空间搜索模块20、变量提取模块30、靶区优选模块40。
[0075]预测体构建模块10、空间搜索模块20、变量提取模块30、靶区优选模块40是相互关联的。[0076]开展大比例尺三维成矿预测，首先需要构建高精度三维多元信息综合预测模型体(栅格体)。该模型体由大量的六面体构成，用于表示初始预测体，之后可以使用基于六面体的定量与定位方法，计算出预测目标体的位置。预测目标体是以初始预测体为基础，通过单元格搜索与变量提取功能得到的。在进行三维成矿预测时，首先需要将不同图层、不同类型的数据提取到三维预测体的属性表中，因此，需要空间搜索模块20和变量提取模块30两个模块实现对初始预测体的筛选。
[0077]预测体构建模块10，用于构建三维成矿预测的初始预测体。该初始预测体可以由表面模型体通过算法转换过来，也可以直接计算构建。该表面模型体由大量的六面体构成，用于表示初始预测体。
[0078]空间搜索模块20，用于以三维空间离散点K阶近邻快速搜索算法为基础，实现快速提取空间离散点等的预测要素。
[0079]变量提取模块30，主要是对基于表面模型体表达的预测要素进行改进的快速栅格化，实现快速提取曲面和预测体等的预测要素变量。
[0080]靶区优选模块40，用于根据证据权重法对三维成矿预测体的初始变量属性表进行筛选与优选，获取预测目标体。
[0081]三维可视化贯穿于整个大比例尺成矿预测过程之中，为预测体构建模块10、空间搜索模块20、快速提取模块30、靶区优选模块40提供逼真的计算机三维可视化表达。
[0082]进一步地，预测体构建模块10包括:
[0083]外包围盒提取模块，用于提取表面模型体或曲面数据的外包围盒；
[0084]外包围盒网格化模块，用于按照一定的网格粒度，将外包围盒分解为一组六面体；
[0085]面体求交模块，用于根据外包围盒进行三角面片与六面体求交；此处的三角面片指构造成体表面的基本图元，即三角形。
[0086]六面体分类模块，基于Flood-Fill算法和求交计算结果，对构成三维成矿预测体的六面体进行快速筛选，生成初始预测体，同时生成六面体模型的属性表，即初始预测体的
属性表。
[0087]进一步地，空间搜索模块20包括:
[0088]第一索引建立模块，用于建立三维离散点的八叉树空间索引，此索引为第一索弓I ;
[0089]第二索引建立模块，用于建立八叉树空间索引的三维R-Tree树索引，此索引为第二索引；
[0090]搜索范围获取模块，用于通过输入与自动计算得到待查点的搜索范围；
[0091]空间插值实现模块，用于根据搜索范围和第一索引建立模块、第二索引建立模块的结果，并应用普通克立格、反距离加权公式进行快速空间插值；
[0092]第二索引是建立在第一索引上的索引。即，基于原始数据建立第一索引。之后，在第一索引的基础上建立第二索引。在使用克立格或者反距离加权计算时，首先使用第二索引；之后，在第二索引计算的基础上，使用第一索引，得到最终的计算结果。
[0093]筛选赋值模块，用于将空间插值结果赋值给初始预测体，并保留在初始预测体的属性表中。[0094]进一步地，变量快速提取模块30包括:
[0095]栅格化处理模块，用于根据预测体构建模块10的算法转换功能，栅格化表面模型体或者曲面数据，且栅格的长、宽、高等参数与预测体构建模块10中的一致；
[0096]结果赋值模块，用于将栅格化结果按照空间位置直接赋值到初始预测体的属性表中。
[0097]筛选赋值模块所实现的赋值是指赋值过程中需要先进行搜索计算，再进行赋值。所谓搜索计算是指插值结果与初始预测体中的六面体并不是一一对应的，因此，需要根据筛选算法确定哪个或者哪些插值结果被如何保留到初始预测体中。而在结果赋值模块中，这一赋值过程是——对应的。
[0098]初始预测体的属性表指预测模型的属性库，它的存在贯穿于整个装置。它的主要目的是保存筛选赋值模块和结果赋值模块计算的结果。
[0099]进一步地，靶区优选模块40包括:
[0100]成矿概率计算模块，用于以初始预测体的属性表为初始数据，基于证据权重法计算三维成矿预测体中每个六面体的成矿概率；
[0101]证据权重法(Bonham-Carter, 1989)是现有技术,是经常用于矿产资源评价工作的方法。
[0102]优选体模块，用于以初始预测体模型为基础，基于OpenGL实现三维交互式编辑功能设定空间筛选范围，并结合成矿概率计算结果，提取成矿预测目标体。
[0103]如图2所示，是三维矿产资源预测多元信息综合处理方法流程图。
[0104]该方法适用于矿山深部和外围找矿的多元信息综合分析模型快速构建与可视化分析等领域，属于数学地质与地理信息系统地图制图范畴，主要用途针对大比例尺矿产勘查过程中获得的地质、地球化学、地球物理等信息，利用计算机三维技术，实现基于多元成矿信息的自动模拟、综合分析与可视化，辅助地质专家解决已知矿山深部及外围找矿的模型及其装置。
[0105]此方法的主要步骤如下:
[0106]步骤201，构建三维成矿预测的初始预测体。
[0107]进一步地，步骤201是下述步骤202至步骤204的基础。步骤203和步骤204都是为由步骤201生成的初始预测体提取初始变量的步骤。步骤201的核心内容是基于Flood-Fill思想，实现矢量体到栅格体的快速转换；步骤204中的基于证据权重算法与三维可视化技术的隐伏靶区优选是在步骤201的基础上进行初始预测体的分割的。
[0108]体的栅格模型又称为块段模型，该模型比较适合表示非均匀体数据的空间分布，在大比例尺深边部成矿预测、三维反演、矿山储量估算等领域得到广泛应用。在大比例尺成矿预测研究中，首先需要建立初始预测栅格体。初始预测栅格体是指在成矿规律与研究区地质背景研究的基础上，初步圈定的一个用于成矿预测的范围。对于这个范围的计算机表达，通常是使用一个体来表示。体数据的来源，是利用体建模功能新建一个，更加普遍的方法是由地表和一个向下延伸距离形成一个体。建立高精度三维成矿预测体需要解决由体的表面模型到栅格模型的快速转换问题。
[0109]步骤202，地质采样点数据快速空间搜索。
[0110]进一步地，步骤202以步骤201的结果为基础，提取地质采样点数据。步骤202的核心内容是基于八叉树与R-Tree树的K-近邻查询算法，可以作为任意空间插值算法的基础；步骤204中将使用由步骤202提取属性的步骤201的结果模型，划分深边部找矿靶区。
[0111]使用地质、地球物理、地球化学等多元离散采样点信息时，一般首先需要进行空间插值。该步骤主要解决在三维空间中待插值点较多，同时存在任意方向搜索需求的前提下，快速实现空间插值的解决方案。本发明主要使用普通克立格与反距离加权法实现快速空间插值；
[0112]步骤203，离散点图层到独立目标图层的成矿信息的提取。
[0113]在步骤201和步骤202的基础上，步骤203将以六面体单元格的中心点为空间待插值点并根据搜索到的成矿要素点，计算每个六面体单元格所对应的成矿要素值。成矿要素值的计算主要包括6种方法。即:最大值、最小值、算数平均值、算法求和、普通克立格插值以及反距离加权法插值等方法获得。离散点成矿信息的提取的核心是离散点的K-近邻查询。
[0114]步骤204，复杂体或曲面图层到独立目标图层的成矿信息提取。
[0115]在本发明中，体和曲面都是由表面模型表示的。复杂体是指包含多个部分但没有三角面片相交的三角网。其中，多个部分的位置关系既可以是相离的，也可以是相互包含的。
[0116]本步骤的核心是步骤201中的算法，将复杂体或者曲面转换为栅格体。唯一的区别在于:栅格体中所包含的六面体单元格不是由用户设定的，而一定是自动确定下来的，参数与步骤201中用户设置的参数相同，即其标准为预测栅格体的单元格大小，同时新生成单元格的外包围盒与初始预测体的外包围盒进行“求并”计算，以确保新生成单元格在空间位置上一定与某一个预测栅格体的单元格重叠，如图9、10所示。同时，将与新生成单元格重叠的预测栅格体单元格对应的属性值设置为1，否则设置为O。
[0117]在图9中，描述了本发明中变量提取模块的示意图，其主要目的为了说明提取体或者曲面属性的过程。首先将体或者曲面栅格化，其次将其与初始预测体的栅格模型进行求交计算；最后保存体或者曲面的存在属性。
[0118]在图10中，描述了本发明中变量提取模块的示意图。主要是针对图9中的方法进行改进。采用了“对齐”策略。主要解决图9中所述的第二步，即“与初始预测体的栅格模型进行求交计算”，这一步将会导致大量计算。当采用“对齐”策略后，不需要进行求交计算，仅需要比较相同位置上是否同时存在两个栅格化的立方格。
[0119]步骤205，基于证据权重算法筛选成矿靶区。
[0120]首先，将每一种地质标志图层都用二态变量来表不，用I表不地质标志存在，O表示地质标志不存在；其次，每一个地质标志都计算一对权系数，一个表示该标志存在时的权，另一个表示该标志不存在时的权。当无法确定该标志存在与否时，令权系数为O ;第三，预测矿种矿床产出的后验概率比(odds)的对数值等于先验概率比的对数值与各种地质标志的权系数之和。该方法源于概率论中的贝叶斯关系式，在计算成矿后验概率时，需要用到有限个独立随机事件的概率乘法公式，因此，每一种控矿地质因素相对于矿床产出这一概率事件来说，都必须是条件独立的。
[0121]步骤206，在地质背景与成矿规律研究的基础上，应用三维可视化，进行交互式编辑，提取成矿靶区。[0122]根据步骤205，地质学家可以根据后验概率对初始预测栅格体进行分类，例如，后验概率大于等于0.95的是A类靶区，在区间[0.8，0.95)上的是B类靶区，在区间[0.6，0.8)上的是C类靶区。但是，在一般情况下，每一类靶区的分布规律可能与地质学家对研究区的地质背景与成矿规律研究不符，此时需要进行一定的交互式编辑工作。即，将一些地质专家认为重要的区域单独圈定出来。在此，需要应用基于OpenGL的三维编辑功能，绘制一个规则的筛选六面体，将此筛选六面体与希望进行筛选的A或B或C类靶区求交，提取最终的矿区深边部矿体。
[0123]采用以上6个步骤，可在一定程度上满足大比例尺三维成矿找矿信息提取与综合分析模型快速构建与三维可视化等方面的需求。下面结合附图对如上方法流程进行具体描述。
[0124]如图3A-3B所示，是R-Tree空间索引模型示意图；图4是本发明大比例尺矿产资源评价多元采样点的快速空间插值方法流程图。
[0125]三维空间离散点的K-近邻查询是三维地理信息系统研究的重要问题。与二维空间离散点的K-近邻查询相比，三维空间离散点对K-近邻查询的速度要求更高，其实现也更加复杂。目前，关于三维空间中K-近邻查询的研究大多建立在某一种空间索引结构的基础上。但是，由于每一种空间索引都存在着不足，因此，使用单一的空间索引结构进行三维空间离散点的K-近邻查询的效果并不十分理想。针对上述问题，本发明提出一种三维空间离散点的K-近邻查询的改进方法。该方法将八叉树和R-Tree树两种空间索引结构复合使用，完成三维空间离散点的K-近邻查询，提高了 K-近邻查询的效率。实验证明，该方法可以快速、准确地进行三维空间离散点的K-近邻查询。在三维空间中，随着点的数量大量增加，K-近 邻查询是非常耗时的。例如，在三维空间中对12000个点进行反距离加权空间插值，设待插值点达到1000万，如果使用线性方法进行插值，其时间复杂度将达到0(12000X10000000)。可以看出上述方法的时间复杂度是用户无法接受的。
[0126]目前，针对此问题开展了一系列快速算法研究。这些算法主要可以分为两大类。第一类，利用点集VOTonoi图来进行K个最近点的搜索，但点集的VOTonoi图的计算量非常大。第二类，利用辅助空间索引结构进行K个最近点的搜索。Gumhold等人提出空间分块策略，但是该方法既不能保证它的空间分块具有最佳或接近于最佳的搜索速度，也不能保证每个数据点都能找到K个最近领域。熊邦书等对该类方法进行了改进，但仍然需要估算分块的大小。刘宇等应用构建三维R-Tree空间索引，提高查找速度，但是当初始采样点很多(^ IO4)时，R-Tree树所花费的时间与空间开销不能忽略。同时，在进行K-近邻查询时，其搜索半径影响域不一定是水平摆放的矩形或椭圆形。特别是在矿产资源评价领域，空间待插值点受到已知三维模型空间形态的约束，例如:矿体产状信息等。因此，在进行空间插值时，还需要考虑搜索半径影响域在空间的形态。
[0127]首先介绍一些相关的概念及术语以便后文论述使用。
[0128]定义1.R-Tree空间索引结构
[0129]R树是一种高度平衡树，是B树在K维(K>=2)空间上的自然扩展。R树由中间节点和叶节点组成，实际数据对象的最小外接矩形存储在叶节点中，中间节点通过聚集其低层节点的外接矩形形成，包含所有这些外接矩形。同时，R树是一种动态索引结构，即:它的查询可与插入或删除同时进行，而且不需要定期地对树结构进行重新组织。如图3Α、3Β所/Jn ο
[0130]图3A表示R-Tree是一棵索引树，可以使用该结果进行空间查询；图3B表示图3A所表述的R-Tree的原始数据的空间位置分布。具体而言，矩形块R8~R19未包含任何其他矩形块，因此，这些矩形块在图3A中为叶子节点。以此类推，由图3B可以得到图3A。
[0131]定义2.八叉树
[0132]将整个空间区域按照2nX2nX2n的方式在空间八个象限内递归分割区域，逐步分解成为被单一类型区域内含的立方体区域。
[0133]变量η是一个确定值，且在一个具体计算中是唯一的值。具体而言，η由初始预测体在三维空间中的长、宽、高中计算2的对数最大的那个值。
[0134]综上所述，总结本部分算法的总体流程如下，如图4所示，描述了大比例尺矿产资源评价多元采样点的快速空间插值方法。
[0135]步骤401，建立三维点云的八叉树空间索引；
[0136]该步骤中，首先，计算点云的外包围盒；其次，将外包围盒按照2η幂方式分割，直到达到默认或者用户指定标准为止。
[0137]变量η是一个确定值，且在一个具体计算中是唯一的值。具体而言，η由初始预测体在三维空间中的长、宽、高中计算2的对数最大的那个值。
[0138]步骤402，建立八叉树空间索引的三维R-Tree树索引；
[0139]该步骤中，是将由八叉树分解空间获得的子区域按照R-Tree树建模规则，构建空间子区域的R-Tree查找树。
[0140]空间八叉树将整个空间区域按照2nX2nX2l^方式在空间八个象限内递归分割区域，获得一组子区域。首先，这些子区域是规则的，但是包含数据的子区域并不是规则的。因此，使用R-Tree管理这些包含数据的子区域。
[0141]步骤403，通过矩阵计算获得待插值点的搜索范围；
[0142]该步骤中，在空间插值计算过程中，首先需要确定待插值点的搜索范围。这一范围，一般使用六面体、球体或椭球体表示。在矿产资源预测评价领域，一般使用六面体或椭球体。在三维空间中，确定六面体和椭球体的三个轴不一定与笛卡尔坐标系的X，Y, Z轴平行；
[0143]此处的矩阵计算指三维空间坐标变换。由上面一段可知，点的搜索使用六面体和椭球体。其搜索有方向，即六面体或者椭球体的空间形态，由倾向、倾角与侧伏角决定。在进行空间数据搜索时，首先根据搜索六面体或者椭球体的空间形态进行空间坐标变换。其次，判断待插值点是否在六面体或者是椭球体中。具体的变换包括平移和旋转。
[0144]步骤404，应用普通克立格或者反距离加权法进行空间插值；
[0145]该步骤中，执行该步骤的主要目的是依据点云数据得到待插值点的属性值。
[0146]根据步骤403中的搜索结果，使用普通克立格或反距离加权法计算待插值点处的值。具体而言，当存在一 个空间待插值点，首先使用步骤403中的搜索算法在点云数据中进行搜索，其次，使用普通克立格或反距离加权法估计待插点的值。这里对普通克立格进行简单介绍。普通克里格(Ordinary Kriging)是区域化变量的线性估计,它假设数据变化成正态分布，认为区域化变量Z的期望值是未知的。插值过程类似于加权滑动平均，权重值的确定来自于空间数据分析。[0147]如图5A、5B、5C所示，是本发明体素之间的相邻关系；图6是本发明基于复杂地质体模型与地表曲面模型约束的找矿信息综合分析模型快速构建方法流程图；图7是本发明Flood-Fill算法二维空间实现示意图；图8是本发明Flood-Fill算法三维空间实现示意图。
[0148]目前，对于三维体形数据模型及其建模方法，国内外已有大量研究。其建模方法主要包括:块段构模法、线框构模法、表面构模法、实体构模法和断面构模法。块段模型的主要特点是对复杂地质体内的不均匀性具有较强的表达能力，并且易于实现复杂地质体的定量计算和空间分析，因此广泛应用于三维反演、矿山储量估算和大比例尺三维成矿预测等矿山实际生产中。
[0149]建立复杂地质体的块段模型，首先，需要完成复杂地质体的线框或表面构模。其次，在此基础上通过空间查询运算完成复杂地质体的块段构模。块段模型在矿山应用中比较适合于预先给定体素(本发明主要指长方体)尺寸的应用，如储量估算和三维隐伏复杂地质体预测工作。即便如此，由于需要保证定量计算的精度，上述应用还是会产生大量的体素单元(>106)。为了解决上述转换问题，国内外学者进行了深入的研究。其核心思想是建立辅助空间索引，加速长方体与复杂地质体表面多边形空间查询计算的速度。Tomaminen提出使用射线法实现表面模型到块段模型的转换，由于需要对每个体素引出的射线进行处理，因此很难满足大量立方格快速计算的需求，同时该方法存在精度问题。在Tomaminen工作的基础上，Amanatides和Jiang提出基于BSP树空间索引结构进一步改进转换算法。对于达到IO6数量级的体素，使用BSP树结构，会造成树的深度过深导致查找效率下降等问题。毕林等和荆永滨等提出使用线性八叉树结构快速构建复杂地质体的块段模型，该方法能够满足快速生成复杂地质体表面体素。但是，在判断复杂地质体内部体素时，仍然延续使用了算法Tomaminen的思想，在数据量很大时，整个算法效率仍然会明显降低。此外，研究人员在实际开发中发现八叉树模型更加适合对已经建立的块段模型进行管理，而在建立块段模型过程中，使用八叉树结构无疑将增加系统开销。
[0150]综上所述，实现复杂地质体块段构模的关键是快速、精确完成由线框或表面模型到块段模型转换。目前研究中存在的主要问题是，当需要处理的长方体Cuboid O IO6)数量较大时，上述算法都没有给出令人满意的建模方法，因此需要提出解决方案。
[0151 ] 首先介绍一些相关的概念及术语以便后文使用。
[0152]定义1.复杂地质体
[0153]体是三维几何元素，使用封闭表面围成的维数一致的有效空间。一般满足上述描述的体称为正则体。即，任何面必须是表面的一部分，不能悬面；每条边有且仅有两个邻面，不能悬空；每个顶点至少与三条边相邻，不能孤立。在此基础上，复杂地质体是由一组体构成的，这些体互不相交，但可以相互包含或者分离。
[0154]定理I复杂地质体模型三角网格边的条数=1.5X三角形的个数
[0155]证明:已知三角形网格所包含的边与三角形存在如下关系:
[0156]3T=2Ei+Eb
[0157]其中，T是三角形的个数，Ei为所有公共边的条数，Eb为非公共边的条数。此公式中的1、b可任意给定。
[0158]又已知在空间中具有拓扑结构的复杂地质体模型为封闭三角网格，则Eb=O,则存在如下关系:3T=2Ei；即Ei=L 5XT
[0159]证明完毕。
[0160]定义2体素之间的相邻关系
[0161]体素可以理解为二维像素在三维空间中的扩展，是一组分布在正交网格中心的立方体单元。两个体素之间分别按照共点、共边和共面等限定条件分为26-邻接、18-邻接和6-邻接，如图5A、5B、5C所示。
[0162]图5A-5C表不空间某一立方体与其它立方体的位置关系。图5A中表不26-邻接；图5B表示18-邻接；图5C表示6-邻接。
[0163]定义3.复杂地质体边界长方体
[0164]设存在长方体Cell，Cell的6_邻接长方体中至少有I个长方体不存在，则称Cell为复杂地质体的边界长方体。
[0165]本发明提出应用均匀网格分割与Flood-fill思想实现复杂地质体表面模型到块段模型的转换。在本发明提出算法中等网格分割算法主要用于剖分复杂地质体的外包围盒。这样做的好处是剖分过程的时间复杂度为0(1)，且剖分后的长方体是有序的，其定位查询时间复杂度也为0(1)。由等网格分割得到的长方体根据其与复杂地质体的关系可以分为3类。第一类，在复杂地质体的外部；第二类，在复杂地质体的内部；第三类，在复杂地质体的表面上。根据体形块段模型的定义可知，块段模型由第二类和第三类长方体组成。第三类长方体中的绝大部分可以在等网格分割的基础上由求交计算快速得到。而第三类的一部分和第二类的长方体则需要应用Flood-fill算法计算获得。
[0166]综上所述，总结本发明提出算法的总体流程如下，如图6所示:
[0167]步骤601，建立复杂地质体的外包围盒；
[0168]该步骤中，取原有地质体在X, Y, Z三个坐标轴的范围[(xmin, ymin, zmin),(xmax, ymax, zmax)]，构建一个六面体。
[0169]步骤602，按照用户输入的长方体尺寸均匀分割复杂地质体的外包围盒；
[0170]该步骤中，采用人机交互操作模式，由用户指定分割外包围盒的精度，即在X，Y，Z轴三个方向上的步长。
[0171]步骤603，按照均匀分割后长方体之间的6-邻接关系判断面元(三角形)与长方体是否相交；
[0172]面元是图形学中描述空间几何对象的基本单位。说明“面元即指三角形”，避免歧义。在其他不同的算法中，有可能面元是矩形。
[0173]该步骤中，6-邻接关系如图5C所示。关于空间线段的求交计算，可以在计算几何书中找到。
[0174]步骤604，在已知复杂地质体表面长方体的基础上，应用扫描线思想检测剩余长方体是否在复杂地质体的内部；
[0175]该步骤中，如下述算法3所示。
[0176]步骤605，保留复杂地质体表面和内部长方体并且加以区分，得到最终的块段模型。
[0177]该步骤中，如下述算法3所示。
[0178]综上所述，根据算法主要思想和总体流程，总结算法的流程如图6所示。[0179]根据如上描述，本发明算法主要包括两大部分。第一部分是查找复杂地质体表面上的长方体，第二部分是查找复杂地质体内部的长方体。
[0180]第一部分主要包括两步:第一步是查找复杂地质体表面边界边上的长方体；第二步是查找复杂地质体表面三角面片上的长方体。
[0181]本发明将使用下述算法I和算法2完成第一部分长方体的查找。
[0182]算法I基于外包围盒等网格分割思想查找边经过的长方体
[0183]设BM是一组由等网格分割得到长方体；
[0184]及是复杂地质体上的边。
[0185]输入:I和 BM
[0186]输出:AB经过的一组长方体CubeArray
[0187]步骤(I)计算点A和点B所在BM中的长方体CellA和CellB ；
[0188]步骤(2)判断长方体CellA和CellB是否相等。若不相等，则计算与长方体CellA的相交矩形P并且转到步骤(3)，否则，记录长方体CellA并转到步骤(4)；
[0189]步骤(3)根据6-邻接定义和相交矩形P，获得与长方体CellA相邻的长方体Cell’，标记长方体Cell’的属性m_bCellStatus=l。若Cell’ ！=CellB则返回步骤(2),否则转到步骤(4)。
[0190]步骤⑷结束。
[0191]长方体Cell’是一个新得到的长方体；其依据与CellA的6_邻接和矩形P获得。
[0192]算法2判断三角面片与长方体是否相交
[0193]输入:CubeArray和三角面片t
[0194]输出:三角面片上的一组长方体cTArray
[0195]步骤I,在CubeArray中查找一个长方体c,c满足至少存在一条边e与三角面片t相交。如果存在这样的长方体C，转到步骤2，否则转到步骤3。
[0196]步骤2，按照6-邻接查找规则，查找一个长方体c’，c’满足至少存在一条边e’(e’古e)与三角面片t相交。如果存在这样的长方体c’，继续查找计算。否则转到步骤
O O
[0197]步骤3，结束。
[0198]在获得与复杂地质体表面上的长方体体素之后，还需要获得其他在块段模型上的长方体。在对第二部分长方体的查找过程中，本发明提出将Flood-fill算法应用到三维空间，即通过查找待判断长方体的6-邻接方向上长方体的内外状态，确定当前长方体是否在复杂地质体内部。由于这个判断过程是基于等网格分割的结果，因此算法所需开销相对较低。算法3的思想如图7和8所示。在二维空间中，如图7，应用Flood-fill算法确定当前网格单元在多边形内外的状态时，需要分别向网格单元格的4个方向进行搜索，依据搜索结果确定待判断网格单元的内外状态。同理，在三维空间中，如图8所示，待判定长方体可以通过其6-邻接方向上长方体的内外状态，确定当前其是否在复杂地质体内部。
[0199]最后，基于Flood-fill算法获得块段模型的剩余长方体。使用扫描线方法标记部分外部长方体。沿射线方向跟踪其经过的长方体Cell，若Cell为复杂地质体表面的长方体，则停止跟踪并且标记Cell为definedCube ;否则，标记Cell为外部的长方体并且继续跟踪直到遇到另一个边界长方体；遍历所有未标记的长方体。设其中任意一个未标记长方体undefinedCube,依据其6_邻接规则沿6个方向做射线,取最近遇到外部长方体的射线r的方向作为undefinedCube判断方向；计算射线r所经过的长方体的数量η判断立方格在地质体的内外。如下面算法3所述。
[0200]算法3基于扫描线算法获得块段模型的剩余长方体
[0201]输入:边界长方体集合S
[0202]输出:块段模型BlockModel
[0203]步骤I由S向外包围盒内部做射线，这些射线平行于X轴或Y轴或Z轴；
[0204]步骤2使用扫描线方法标记部分外部长方体。沿射线方向跟踪其经过的长方体Cell，若Cell为矿体表面的长方体，则停止跟踪，并且标记Cell为definedCube ;否则，标记Cell为外部的长方体并且继续跟踪直到遇到另一个边界长方体；
[0205]步骤3遍历所有未标记的长方体。设其中任意一个未标记长方体undefinedCube,依据其6-邻接规则沿6个方向做射线，取最近遇到外部长方体的射线r的方向作为undefinedCube 判断方向；
[0206]步骤4计算射线r所经过的长方体的数量η,如果η为奇数,则undefineCube为内部长方体；否则，UndefinedCell为外部长方体；
[0207]步骤5结束。
[0208]本发明提出的方法，适用于大比例尺三维成矿预测中预测要素模拟与多元找矿信息综合分析模型快速构建及其可视化分析等领域，主要用于针对大比例尺矿产勘查过程中获得的各种数据，例如采样点数据、地质地形数据等，利用快速空间插值、曲面自动拟合、栅格模型快速构建等三维计算机技术，实现地质、地球物理、地球化学等多元数据的三维信息自动模拟、综合分析与可视化显示。
【权利要求】
1.一种三维矿产资源预测评价中信息综合处理装置，其特征在于，包括: 预测体构建模块，用于构建三维成矿预测的初始预测体； 空间搜索模块，连接所述预测体构建模块，用于对初始预测体快速提取空间离散点的预测要素，将结果赋值给初始预测体并保留在初始预测体的属性表中； 变量提取模块，连接所述预测体构建模块，用于对由表面模型表达的预测要素进行栅格化，并将栅格化结果赋值到初始预测体的属性表中； 靶区优选模块，连接所述空间搜索模块、所述变量提取模块，用于对初始预测体的属性表进行筛选与优选，获取预测目标体。
2.根据权利要求1所述的装置，其特征在于，所述预测体构建模块进一步包括: 外包围盒提取模块，用于提取体表面或曲面数据的外包围盒； 外包围盒网格化模块，用于按照设定的网格粒度将外包围盒分解为六面体； 面体求交模块，用于根据外包围盒进行三角面片与六面体求交计算； 六面体分类模块，用于基于种子填充法和求交计算结果，对构成三维成矿预测体的六面体进行快速筛选，生成初始预测体以及初始预测体的属性表。
3.根据权利要求1所述的装置，其特征在于，所述空间搜索模块进一步包括: 第一索引建立模块，用于建立三维离散点的八叉树空间索引； 第二索引建立模块，用于建立八叉树空间索引的三维R-Tree树索引； 搜索范围获取模块，用于获取待查点的搜索范围； 空间插值实现模块，用于根据搜索范围和第一索引建立模块、第二索引建立模块的结果，并应用克立格、反距离加权公式进行空间插值； 搜索赋值模块，用于将空间插值结果赋值给初始预测体，并保留在初始预测体的属性表中。
4.根据权利要求1、2或3所述的装置，其特征在于，所述变量快速提取模块进一步包括: 栅格化处理模块，用于栅格化表面模型体或曲面数据； 结果赋值模块，用于将栅格化结果按照空间位置赋值到初始预测体的属性表中。
5.根据权利要求1、2或3所述的装置，其特征在于，所述靶区优选模块进一步包括: 成矿概率计算模块，用于以初始预测体的属性表为初始数据，基于证据权重法计算三维预测体中每个六面体的成矿概率； 优选体模块，用于以初始预测体模型为基础，并结合成矿概率计算结果，提取成矿预测目标体。
6.一种三维矿产资源预测评价中信息综合处理方法，其特征在于，包括: 步骤一，构建三维成矿预测的初始预测体； 步骤二，对初始预测体快速提取空间离散点的预测要素，将结果赋值给初始预测体并保留在初始预测体的属性表中； 步骤三，对由表面模型表达的预测要素进行栅格化，并将栅格化结果赋值到初始预测体的属性表中； 步骤四，对初始预测体的属性表进行筛选与优选，获取预测目标体。
7.根据权利要求6所述 的方法，其特征在于，所述步骤一，进一步包括:提取体表面或曲面数据的外包围盒； 按照设定的网格粒度将外包围盒分解为六面体； 根据外包围盒进行三角面片与六面体求交计算； 根据种子填充法和求交计算结果，对构成三维成矿预测体的六面体进行快速筛选，生成初始预测体以及初始预测体的属性表。
8.根据权利要求6所述的方法，其特征在于，所述步骤二，进一步包括: 建立三维离散点的八叉树空间索引； 建立八叉树空间索引的三维R-Tree树索引； 获取待查点的搜索范围； 根据搜索范围和第一索引建立模块、第二索引建立模块的结果，并应用克立格、反距离加权公式进行空间插值； 将空间插值结果赋值给初始预测体，并保留在初始预测体的属性表中。
9.根据权利要求6、7或8所述的方法，其特征在于，所述步骤三，进一步包括: 栅格化表面模型体或曲面数据； 将栅格化结果按照空间位置赋值到初始预测体的属性表中。
10.根据 权利要求6、7或8所述的方法，其特征在于，所述步骤四，进一步包括: 以初始预测体的属性表为初始数据，基于证据权重法计算三维预测体中每个六面体的成矿概率； 以初始预测体模型为基础，并结合成矿概率计算结果，提取成矿预测目标体。
【文档编号】G06Q50/02GK103903061SQ201410077136
【公开日】2014年7月2日 申请日期:2014年3月4日 优先权日:2014年3月4日
【发明者】李楠, 肖克炎, 邹伟, 孙莉, 丁建华, 娄德波, 阴江宁, 范建福 申请人:中国地质科学院矿产资源研究所, 李楠