CN113096240B - 一种局部异向性椭球体动态生成方法 - Google Patents

Abstract

本发明公开了一种局部异向性椭球体动态生成方法，包括以下步骤：S1，对待测区矿体进行网格划分，获取网格与待测区矿体的交点坐标，对三维模型的坐标节点进行遍历搜索，储存所有节点的三维空间坐标；S2，从三维空间坐标集中搜索得到与网格与待测区矿体的交点坐标集中的交点坐标最近的数据点生成当前点的搜索椭球体模型，进而生成所有搜索椭球体模型；S3，生成矿体三维模型的中心线并获取坐标，计算网格与中心线的交点，获取交点坐标；S4，获取中心线上所有搜索椭球体模型；S5，合并边界和中心线的搜索椭球体模型，利用插值法即可得到整个矿体内任意位置的搜索椭球体模型。通过动态生成椭球体，保证了储量估算结果的精度和可信度。

Description

一种局部异向性椭球体动态生成方法

技术领域

本发明涉及地质统计学矿产资源储量估算方法的研究与应用领域，特别是涉及一种局部异向性椭球体动态生成方法。

背景技术

矿产资源储量是经过矿产资源勘查和可行性评价工作所获得的矿产资源蕴藏量的总称，是指矿物含量在指定品位边界以上，并集中埋藏的矿产数量。矿产资源储量估算是地质勘探工作中的一项重要环节，它是根据各勘探阶段中所获得的相关的资料和数据，根据一定的计算方法，来确定矿产资源的数量、质量、空间分布、开采技术条件及其研究精度的过程。储量估算贯穿着矿山的整个生命周期，矿产资源储量估算是对矿产资源“质”与“量”的全面评估与总结，是生产建设和企业投资的主要依据。

目前国际上主流的储量估算方法有距离幂次反比法和地质统计学的克里格法等。其中，克里格法是一种以区域化变量为核心、以变异函数为基本工具的无偏的最优估计方法。距离幂次反比法假设矿石品位之间的相关性与其空间距离呈反比的关系。这两种方法的基本实施流程如下：

1.数据分析及预处理

数据分析和预处理过程主要包括：

(1)基本统计分析：对样品数据的品位、样长等进行基本的统计分析；

(2)特高品位处理：一般情况下，特高品位会对最终估算结果产生较大影响。因此通常要其进行识别和替换处理；

(3)样品组合：其目的是使数据依附固定在相同大小的承载上。通常将组合后的样品长度定为大部分样品的长度或开采台阶高度；

(4)空间变异性分析：利用变异函数计算以获取矿石品位在估算域的连续性方向及变化程度。空间变异性计算结果可以为距离幂次反比法的参数设置提供一定的参考依据，它也是克里格法插值计算的必须参数。

2.建立块模型

根据矿化域或矿体空间范围、形态和产状，以一定长、宽、高的块为单位进行填充，形成整个估算区域的块模型，块体模型中每个单元块体将在步骤4中被赋予矿石品位估计值。

3.建立变异函数模型

通过计算各方向的变异函数，确定估算域的最大连续性方向(主轴)、次连续性方向(次轴)和最小连续性方向(短轴)，三个方向在空间上应相互垂直。

4.搜索椭球体的设置

搜索椭球体定义了对待估块进行品位估算时，样品数据的空间搜索范围和相关参数。一般情况下，要求其空间方位与矿体的形态、状态相匹配。例如：图1所示的矿体与对应的搜索椭球体模型就存在较好的一致性。

但矿体形态如果存在较大的局部变化，这种一致性关系在传统搜索椭球体的框架下便很难实现。例如：假设图1所示的矿体实际上是由一个褶皱控制(图2)，这种情况下之前设置的搜索椭球体便不再适合。

复杂矿体形态情况下搜索椭球体的设置问题，目前的解决方式主要有：

(1)将矿体切分成多个相对规则的区域，在每个区域单独搜索椭球体。

(2)手工或者自动设置多个搜索椭球体模型，然后在矿体模型的不同部位动态调用。这些方法在效率和效果方面，均难以令人满意。

5.交叉验证

为了使估算的结果更加可靠，在估算之前需对拟合得到的变异函数模型的最优性进行校验，即交叉验证。交叉验证的目的是检验变异函数模型参数和搜索椭球体参数的合理性。

6.品位估值及储量汇总

基于前面步骤的模型参数，利用距离幂次反比法或格里格法对块体模型中每个块单元进行品位估值，统计汇总即得到矿区储量结果。

可见，目前的储量流程中搜索椭球体的设置有待进一步改进：在遇到的矿体形状不规则时，传统流程中参数固定的搜索椭球体一般难以满足要求。

发明内容

本发明的目的是提供一种局部异向性椭球体动态生成方法，以解决上述现有技术存在的问题，使地质统计学矿产资源储量估算过程中，在矿体的形态、产状存在明显局部变化的情况下，能生成有效搜索椭球体。

为实现上述目的，本发明提供了如下方案：

本发明提供一种局部异向性椭球体动态生成方法,具体包括以下步骤：

S1，对待测区矿体进行网格划分，计算网格与待测区矿体的边界交点，获取边界交点坐标，基于边界交点坐标，构建第一坐标数据集，基于第一坐标数据集构建矿体三维模型，对三维模型的坐标节点进行遍历搜索，获得坐标节点的三维空间坐标，通过三维空间坐标，构建第二坐标数据集；

S2，基于第一坐标数据集与第二坐标数据集，获得边界交点坐标相邻的数据点，通过数据点构建数据点搜索椭球体模型，基于数据点搜索椭球体模型，构建基于交点的第一搜索椭球体模型；

S3，通过计算获得矿体三维模型的中心线并获取中心线的坐标，基于中心线的坐标构建第三坐标数据集，计算网格与中心线的交点，获取中心线交点坐标，基于中心线交点坐标构建第四坐标数据集；

S4，基于第四坐标数据集与第三坐标数据集，获得中心线交点坐标相邻的数据点，通过数据点构建数据点搜索椭球体模型，基于数据点搜索椭球体模型，构建基于中心线的交点的第二搜索椭球体模型；

S5，基于第一搜索椭球体模型、第二搜索椭球体模型，利用插值法获得第三搜索椭球体模型，用于估算待测区矿体的矿产资源储量，保证储量估算结果的精度和可信度。

进一步地，S1中网格划分的方向为：获得待测区矿体的走向，垂直于走向对待测区矿体进行网格划分。

进一步地，S2中获得数据点搜索椭球体模型的方法为：

S2.1，构建若干个坐标点的空间分布关系，即构建协方差矩阵，计算交点坐标集中的X,Y,Z分量之间形成的协方差矩阵；

S2.2，计算协方差矩阵的特征值和特征向量；

S2.3，将特征值按照从小到大进行排序，取最大的三个特征值对应的特征向量作为当前点搜索椭球体模型的方位；

S2.4，基于方位和指定的轴长，获得数据点搜索椭球体模型。

进一步地，三个特征值对应的特征向量经过归一化处理。

进一步地，S3中获得中心线的方法为：

S3.1，构建代价函数，选取若干个点作为代价函数初始顶点；

S3.2，用动态规则求解使代价函数达到最大值的折线段；

S3.3，在求解出的折线段上利用线性插值插入新的顶点；

S3.4，重复S3.2和S3.3，直到折线段不再发生改变，停止重复，获得目标折线段，目标折线段为中心线。

进一步地，代价函数为：

其中，p_ip_i+1表示线段，a_i表示线段p_ip_i+1的方向，|Δs_i|为线段p_ip_i+1的长度，S代表线段p_ip_i+1的像素集合，

代表概率分布图上线段p_ip_i+1的灰度函数。

进一步地，S5中的所用到的插值法为空间插值算法。

本发明公开了以下技术效果：

本发明在矿体形态复杂的情况下，通过生成实际矿体模型边界和中心线上的椭球体，并以它们作为样品数据对矿体其它位置的椭球体模型进行估值，可以保证生成的搜索椭球体模型与复杂矿体模型的形态、结构具有较好的一致性；

通过动态生成椭球体，使每个局部块体在品位估值过程中都可以用到与矿体局部特征相匹配的搜索椭球体模型，从而保证了储量估算结果的精度和可信度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为简单形态的矿体与静态搜索椭球体模型示意图；

图2为复杂形态的矿体与静态搜索椭球体模型示意图；

图3为矿体模型格网化示意图；

图4为生成矿体模型边界与网格交点处的搜索椭球体示意图；

图5为生成矿体中心线示意图；

图6为生成矿体中心线与网格交点处的搜索椭球体示意图；

图7为插值生成任意位置的搜索椭球体示意图；

图8为全矿区搜索椭球体动态生成效果示意图。

具体实施方式

现详细说明本发明的多种示例性实施方式，该详细说明不应认为是对本发明的限制，而应理解为是对本发明的某些方面、特性和实施方案的更详细的描述。

应理解本发明中所述的术语仅仅是为描述特别的实施方式，并非用于限制本发明。另外，对于本发明中的数值范围，应理解为还具体公开了该范围的上限和下限之间的每个中间值。在任何陈述值或陈述范围内的中间值以及任何其他陈述值或在所述范围内的中间值之间的每个较小的范围也包括在本发明内。这些较小范围的上限和下限可独立地包括或排除在范围内。

除非另有说明，否则本文使用的所有技术和科学术语具有本发明所属领域的常规技术人员通常理解的相同含义。虽然本发明仅描述了优选的方法和材料，但是在本发明的实施或测试中也可以使用与本文所述相似或等同的任何方法和材料。本说明书中提到的所有文献通过引用并入，用以公开和描述与所述文献相关的方法和/或材料。在与任何并入的文献冲突时，以本说明书的内容为准。

在不背离本发明的范围或精神的情况下，可对本发明说明书的具体实施方式做多种改进和变化，这对本领域技术人员而言是显而易见的。由本发明的说明书得到的其他实施方式对技术人员而言是显而易见得的。本申请说明书和实施例仅是示例性的。

关于本文中所使用的“包含”、“包括”、“具有”、“含有”等等，均为开放性的用语，即意指包含但不限于。

本发明中所述的“份”如无特别说明，均按质量份计。

实施例1

本发明提供一种局部异向性椭球体动态生成方法,具体包括以下步骤：

S1，对待测区矿体进行网格划分，计算网格与待测区矿体的边界交点，获取边界交点坐标，基于边界交点坐标，构建第一坐标数据集，基于第一坐标数据集构建矿体三维模型，对三维模型的坐标节点进行遍历搜索，获得坐标节点的三维空间坐标，通过三维空间坐标，构建第二坐标数据集；

S2，基于第一坐标数据集与第二坐标数据集，获得边界交点坐标相邻的数据点，通过数据点构建数据点搜索椭球体模型，基于数据点搜索椭球体模型，构建基于交点的第一搜索椭球体模型；

S3，通过计算获得矿体三维模型的中心线并获取中心线的坐标，基于中心线的坐标构建第三坐标数据集，计算网格与中心线的交点，获取中心线交点坐标，基于中心线交点坐标构建第四坐标数据集；

S4，基于第四坐标数据集与第三坐标数据集，获得中心线交点坐标相邻的数据点，通过数据点构建数据点搜索椭球体模型，基于数据点搜索椭球体模型，构建基于中心线的交点的第二搜索椭球体模型；

S5，基于第一搜索椭球体模型、第二搜索椭球体模型，利用插值法获得第三搜索椭球体模型，用于估算待测区矿体的矿产资源储量，保证储量估算结果的精度和可信度。

S1中网格划分的方向为：获得待测区矿体的走向，垂直于走向对待测区矿体进行网格划分。

S2中获得数据点搜索椭球体模型的方法为：

S2.1，构建若干个坐标点的空间分布关系，即构建协方差矩阵，计算交点坐标集中的X,Y,Z分量之间形成的协方差矩阵；

S2.2，计算协方差矩阵的特征值和特征向量；

S2.3，将特征值按照从小到大进行排序，取最大的三个特征值对应的特征向量作为当前点搜索椭球体模型的方位；

S2.4，基于方位和指定的轴长，获得数据点搜索椭球体模型。

三个特征值对应的特征向量经过归一化处理。

S3中获得中心线的方法为：

S3.1，构建代价函数，选取若干个点作为代价函数初始顶点；

S3.2，用动态规则求解使代价函数达到最大值的折线段；

S3.3，在求解出的折线段上利用线性插值插入新的顶点；

S3.4，重复S3.2和S3.3，直到折线段不再发生改变，停止重复，获得目标折线段，目标折线段为中心线。

进一步地，代价函数为：

其中，p_ip_i+1表示线段，a_i表示线段p_ip_i+1的方向，|Δs_i|为线段p_ip_i+1的长度，S代表线段p_ip_i+1的像素集合，

代表概率分布图上线段p_ip_i+1的灰度函数。

S5中的所用到的插值法为空间插值算法。

实施例2

本发明公开了一种局部异向性椭球体动态生成方法，具体包括以下步骤：

(1)对研究区矿体进行网格划分，下述步骤(2)-(12)均是在相应的每个格网中进行的；

(2)计算格网与矿体区域的交点，获取交点坐标，记为点集P₁；

(3)利用下面的步骤(4)-(8)，在集合P₁中所有数据点位置处生成相应的搜索椭球体模型，记录为集合E₁；

(4)获取矿体三维模型的结点坐标：选择需要进行储量估算的矿体三维模型，对其坐标结点进行遍历搜索，存储所有结点的三维空间坐标，记为点集S₁；

(5)对于P₁中的每一个点，从S₁中搜索得到与其最近的M个数据点。利用步骤(6)-(8)生成当前点的椭球体方位；

(6)构建M个坐标点的空间分布关系(协方差矩阵)：计算结点坐标集中的X,Y,Z分量之间形成的协方差矩阵；

(7)计算协方差矩阵的特征值和特征向量；

(8)有效特征值和特征向量的提取：将特征值按照从小到大进行排序；取最大的三个特征值对应的特征向量，作为当前点椭球体的方位；按照指定轴长生成椭球体模型；

(9)生成矿体模型的中心线，获取矿体三维模型的中心线坐标，将其记为点集S₂；

(10)计算区域划分格网与矿体中心线的交点，获取交点坐标，记为点集P₂；

(11)对于P₂中的每一个点，从S₂中搜索得到与其最近的M个坐标点。重复执行步骤(6)-(8)为P₂中所有点生成一个搜索椭球体模型，记录为集合E₂；

(12)合并上述椭球体模型集合E₁和E₂，记录合并得到椭球体集合为E；

(13)以椭球体集合E为已知样品数据，利用插值算法便可得到整个矿体区域内的任意位置的椭球体模型。

所谓步骤(1)，指的是：根据进行估算的研究区域，估算出整个研究区域范围内矿体的走向，然后垂直于该走向方位对研究区域进行格网化，进而将进行估算的整个矿体进行划分(如附图3)。

步骤(3)，指的是：对区域进行格网化以后，会得到格网与矿体边界相交的结点，在相应的结点处得到相应的结点的三维空间坐标，并生成相应的搜索椭球体模型(见附图4)。

步骤(8)，指的是：根据上面步骤中得到的协方差矩阵来计算搜索椭球体的空间方位，协方差矩阵的特征向量也就是所求搜索椭球体三个轴的方向向量，计算出它的特征向量并将其归一化，这些向量就是搜索椭球体的方向向量：v⁰，v¹，v²。(如附图4)

步骤(9)，指的是：模型中心线提取的思路为：模型中心线用一条含n个顶点的折线段描述；折线段上的顶点可以围绕其初始位置(x_i,y_i)移动，构成一个由折线段组成的集合；模型提取过程为运用动态规划求解使代价函数达到最大值的折线段，将这条折线段视为被提取的模型中心线。设折线段的顶点为p＝{p₁,…,p_n},p_i＝(x_i,y_i)，

式(1)中，S为线段p_ip_i+1的像素集合；

为概率分布图上线段p_ip_i+1的灰度函数。

C_g＝|a_i-a_i+1|＜T₁ (3)

式(2)中，a_i为线段p_ip_i+1的方向；|Δs_i|为线段p_ip_i+1的长度。它们的计算式分别为：

综上，模型中心线可用形式如式(6)的代价函数和一个限制条件(3)表示：

分析式(6)可知，代价函数一系列E_i的和，每一个E_i依赖于折线段的3个相连顶点{p_i-1,p_i,p_i+1}，即每一个顶点p_i只和前一个p_i-1及后一个p_i+1相关。因此，代价函数(6)的最大值可以利用动态规划来求解。

模型中心线的具体提取过程为：

首先，人工选取一系列种子点，作为代价函数(6)中的初始顶点；

然后，运用动态规划求解使代价函数达到最大值的折线段；

其次，利用线性插值在求解出的折线段上插入新的顶点，并再次运用动态规划求解；

最后，重复前两步的过程，直到折线段不再发生变化，将该折线段作为待提取的模型中心线(如附图5)。

所谓步骤(13)，指的是：结合上述步骤得到的边界椭球体以及中心椭球体，利用空间插值算法，例如克里格法，可以插值生成研究区域内任意位置搜索椭球体模型参数(见附图7)。

以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (7)

1.一种局部异向性椭球体动态生成方法，其特征在于：包括以下步骤：

S1，对待测区矿体进行网格划分，计算所述网格与所述待测区矿体的边界交点，获取边界交点坐标，基于所述边界交点坐标，构建第一坐标数据集，基于所述第一坐标数据集构建矿体三维模型，对所述三维模型的坐标节点进行遍历搜索，获得所述坐标节点的三维空间坐标，通过所述三维空间坐标，构建第二坐标数据集；

S2，基于所述第一坐标数据集与所述第二坐标数据集，获得所述边界交点坐标相邻的数据点，通过所述数据点构建数据点搜索椭球体模型，基于所述数据点搜索椭球体模型，构建基于所述交点的第一搜索椭球体模型；

S3，通过计算获得所述矿体三维模型的中心线并获取所述中心线的坐标，基于所述中心线的坐标构建第三坐标数据集，计算所述网格与所述中心线的交点，获取中心线交点坐标，基于所述中心线交点坐标构建第四坐标数据集；

S4，基于所述第四坐标数据集与所述第三坐标数据集，获得所述中心线交点坐标相邻的数据点，通过所述数据点构建数据点搜索椭球体模型，基于所述数据点搜索椭球体模型，构建基于所述中心线的交点的第二搜索椭球体模型；

S5，基于所述第一搜索椭球体模型、第二搜索椭球体模型，利用插值法获得第三搜索椭球体模型，用于估算所述待测区矿体的矿产资源储量，保证储量估算结果的精度和可信度。

2.根据权利要求1所述的局部异向性椭球体动态生成方法，其特征在于：所述S1中网格划分的方向为：获得待测区矿体的走向，垂直于所述走向对待测区矿体进行网格划分。

3.根据权利要求1所述的局部异向性椭球体动态生成方法，其特征在于：所述S2中获得数据点搜索椭球体模型的方法为：

S2.1，构建若干个坐标点的空间分布关系，即构建协方差矩阵，计算交点坐标集中的X,Y,Z分量之间形成的协方差矩阵；

S2.2，计算所述协方差矩阵的特征值和特征向量；

S2.3，将所述特征值按照从小到大进行排序，取最大的三个特征值对应的特征向量作为当前点搜索椭球体模型的方位；

S2.4，基于所述方位和指定的轴长，获得所述数据点搜索椭球体模型。

4.根据权利要求3所述的局部异向性椭球体动态生成方法，其特征在于：所述三个特征值对应的特征向量经过归一化处理。

5.根据权利要求1所述的局部异向性椭球体动态生成方法，其特征在于：所述S3中获得所述中心线的方法为：

S3.1，构建代价函数，选取若干个点作为所述代价函数初始顶点；

S3.2，用动态规则求解使所述代价函数达到最大值的折线段；

S3.3，在求解出的所述折线段上利用线性插值插入新的顶点；

S3.4，重复S3.2和S3.3，直到所述折线段不再发生改变，停止重复，获得目标折线段，所述目标折线段为所述中心线。

6.根据权利要求5所述的局部异向性椭球体动态生成方法，其特征在于：所述代价函数为：

其中所述，p_ip_i+1表示线段，a_i表示线段p_ip_i+1的方向，|Δs_i|为线段p_ip_i+1的长度，S代表线段p_ip_i+1的像素集合，

代表概率分布图上线段p_ip_i+1的灰度函数。

7.根据权利要求1所述的局部异向性椭球体动态生成方法，其特征在于：所述S5中的所用到的插值法为空间插值算法。

Images