CN103824329B - 一种地质勘探三维可视化储量估算方法 - Google Patents

Abstract

本申请提供了一种地质勘探三维可视化储量估算方法，该方法包括以下步骤：首先计算矿体的产状方向、走向和倾向，将矿体投影到水平或垂直面上;接着根据投影后的钻孔位置、平均品位和见矿真厚度，将投影后的钻孔进行块段的划分，形成四边形块段；再通过每个钻孔的平均品位、见矿真厚度和四边形块段的面积，计算每个块段对应的矿体的储量；最后对矿体储量进行分级，估算矿体的总储量。

Description

一种地质勘探三维可视化储量估算方法

技术领域

本发明涉及地质勘探三维可视化储量估算方法，属于地球探测与信息技术领域。

背景技术

在地质勘探过程中，通过矿产勘探工作形成了地表、地下大量的三维空间数据。利用计算机三维可视化技术，可以将这些三维空间数据形成三维地表及地下模型，实现三维相关的矿山管理、场景漫游、空间分析、储量估算等功能。

目前，矿产储量估算方法主要分为两大类，一类是传统的以几何计算为基础的常用方法；一类是以统计学为基础的数学地质方法。传统方法在不同条件下的具体估算方法有：等值线法、地质剖面法、算术平均法、地质块段法、多角形法等。数学地质方法主要有距离平方反比法、地质统计学方法（克里格法）、SD储量估算法等。我国矿山在勘查和开采阶段都普遍采用传统矿产储量估算方法，西方主要矿业国家则采用传统方法和地质统计学方法，其中地质统计学方法应用较广。

先进的地质统计学方法引入我国三十多年来，虽然取得了一定的效果，但由于传统观念的约束以及它不太适合我国中小矿多、贫矿多的国情，所以，一直难以推广。相比之下，传统的储量估算方法具有应用广泛、方法简便、易于掌握、能充分发挥地质工作者主观能动性的优点，非常适合我国的实际情况。特别是在早期普查、详查阶段，勘探密度不大等情况下，该方法具有很大的优势，可在仅获得有限的勘探资料的基础上，迅速开展储量估算。因此，在我国目前的经济技术水平下，传统的储量估算方法仍然是矿山矿产储量估算的主要方法，在我国有着广泛的应用前景和其它方法不可替代的优点。

地质块段法在目前传统地质储量估算中占据重要地位，该方法在国内矿产储量估算中应用了几十年，已经成为国家标准。它具有算术平均法的优点，不需做复杂的分析，其估算结果能满足地勘行业对矿产资源的报告要求。当前，地质人员往往使用Excel等表格工具，将大量钻孔岩性数据输入，实现手动计算。在钻孔、剖面、见矿岩性段等各种数据量繁复的情况下，容易产生数据上的错误。调研后发现，能够完成从三维矿体到地质块段法储量报告的实用软件尚未出现。

在申请号为200710053058.2发明专利中公开了一种固体矿床三维可视化储量计算系统及计算方法。它包括以下子系统：（1）矿区点源数据库子系统：该子系统负责包括钻孔、平硐和样品在内的原始数据的采集、入库，为储量计算提供原始数据；（2）矿区各类二维地质图件自动编绘子系统：对入库的原始数据进行交互编辑，生成矿体三维建模和储量计算所需要的剖面图数据和品位图数据；（3）矿体、矿床和复杂地质结构的三维空间分析子系统：在前面两个子系统所提供的数据基础上，完成三维矿体建模和三维地层建模，并对矿体进行进一步的四面体剖分；（4）矿床、矿体储量计算和三维可视化表达子系统：用克吕格或距离幂次反比法对样品品味进行插值计算，从而对矿体中的每个四面体进行赋值并计算出储量，以三维可视化手段表达出来。实际发现，使用该申请的方法计算矿体储量的精确度不高。

发明内容

本发明的目的是提供了一种根据三维实体数据和采样岩性数据完成储量报告的方法，该方法在储量估算的过程中，对每一个中间步骤都可以实时输出，检查钻孔采样值、块段划分、块段平均品位、储量级别等信息，在提高储量估算效率的同时，提高了估算的准确度。

本发明提供一种地质勘探三维可视化储量估算的方法，该方法包括以下步骤：计算矿体的产状方向、走向和倾向，并将矿体投影到水平或垂直面上;根据投影后的钻孔位置、平均品位和见矿真厚度，将投影后的钻孔进行块段的划分，形成四边形块段；其中，矿体边界处有可能形成三角形块段或多边形块段；通过每个钻孔的平均品位、见矿真厚度和四边形块段的面积，计算每个块段对应的矿体的储量；对矿体储量进行分级，估算矿体的总储量。

本发明在分析传统地质块段法计算过程的基础上，将三维可视化技术和传统地质储量估算方法相衔接，利用三维数据建模形成的三维钻孔、剖面、矿体等数据，完成地质块段法储量估算。该方法的优点是在储量估算的过程中，对每一个中间步骤都可以实时输出，检查钻孔采样值、块段划分、块段平均品位、储量级别等信息，在提高储量估算效率的同时，提高了估算的准确度。

附图说明

图1表示在三维视图下，某矿区的矿体图的三维效果。

图2表示经过计算得到的的矿体产状的走向、倾向和倾角。

图3表示图1的矿体在水平投影平面上的投影结果；

图4表示图3的放大效果图，每个钻孔的名称，钻孔的平均品位和见矿总真厚度如图所示。

为了能明确实现本发明的实施例的结构，在图中标注了特定的尺寸、结构和器件，但这仅为示意需要，并非意图将本发明限定在该特定尺寸、结构、器件和环境中，根据具体需要，本领域的普通技术人员可以将这些器件和环境进行调整或者修改，所进行的调整或者修改仍然包括在后附的权利要求的范围中。

具体实施方式

下面结合附图和具体实施例对本发明提供的方法进行详细描述。

在以下的描述中，将描述本发明的多个不同的方面，然而，对于本领域内的普通技术人员而言，可以仅仅利用本发明的一些或者全部结构或者流程来实施本发明。为了解释的明确性而言，阐述了特定的数目、配置和顺序，但是很明显，在没有这些特定细节的情况下也可以实施本发明。在其他情况下，为了不混淆本发明，对于一些众所周知的特征将不再进行详细阐述。

1.准备数据

首先获得三维钻孔数据、剖面数据和矿体数据。其中，钻孔数据包含了参与三维建模的所有钻孔，每个钻孔由若干岩性段组成，每个岩性段又包含了样品分析值，即品位；剖面数据给出了钻孔的排布方式，用于对生成块段的约束；矿体数据是通过三维建模的方法，使用钻孔、剖面、断层、地表等数据生成的，用于对投影后的钻孔位置做有效性判断，同时对储量估算的范围进行约束，图1是在三维视图下，打开钻孔数据、剖面数据和矿体数据的三维效果。与水平面垂直的各条线段表示每个岩性段，与水平面相垂直的各自之间相互平行的平面表示剖面。

2.计算矿体的产状方向、走向和倾向

2.1通过三维矿体数据计算矿体的产状方向、走向和倾向

本实施例借鉴了OBB即：有向包围盒方法，将建立的OBB树的核心过程提取出来，用于计算矿体的产状、走向和倾向，然后进一步定位最合适的投影平面。本发明中的三维矿体模型，使用表面模型的三角网组织形式，一个三维矿体数据，是由N个三角形组成的，其中N≥1。若数据中存在多边形构成的三维矿体表面，可以使用Delaunay三角剖分的方法，将多边形表面分割成若干三角形表面。因此，这样的数据输入形式，增强了数据类型的适应性。使用OBB计算矿体的产状、走向和倾向的过程如下：

以第i个三角形为例，三个顶点分别为pⁱ、qⁱ和rⁱ。则三维矿体的点均值和协方差矩阵分别为如下（1）式和（2）式：

$\mathrm{\μ}=\frac{1}{3N}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{p}^i+q^i+r^i---\left(1\right)$

$C_{\mathrm{jk}}=\frac{1}{3N}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}({\stackrel{-i}{p}}_j{\stackrel{-i}{p}}_k+{\stackrel{-i}{q}}_j{\stackrel{-i}{q}}_k+{\stackrel{-i}{r}}_j{\stackrel{-i}{r}}_k)---\left(2\right)$

其中，1≤j,k≤3， $\stackrel{-i}{p}=p^i-\mathrm{\μ},\stackrel{-i}{q}=q^i-\mathrm{\μ},\stackrel{-i}{r}=r^i-\mathrm{\μ}.$

这里，pⁱ，qⁱ，rⁱ，都是3×1的向量，即μ表示矿体重心，分别表示三角形的三个顶点相对于矿体重心的偏移方向。C_jk是3×3的协方差矩阵。因为C_jk是3×3的对称方阵，所以，它的特征向量为3个，并且是正交的。根据OBB模型的特征，C_jk的三个特征向量分别对应了矿体的产状方向、走向和倾向。若产状方向近似为水平方向，则选择水平投影；若产状方向近似为垂直方向，则计算经过走向方向的垂直面，将该垂直面设定为投影平面，做垂直投影。对于水平投影，投影平面的法线方向为（0,0,1）；对于垂直投影，投影平面的法线方向为（a,b,c），其中，a,b,c为矿体的产状方向，即C_jk的最大特征向量方向。将投影平面记为α，其解析式可记为Ax+By+Cz+D=0，投影平面的绝对位置并不影响矿体储量。因此，取D=0。得到：

$\left\{\begin{array}{c}A=a\\ B=b\\ C=c\\ D=0\end{array}\right.$

结合图2，1表示产状方向，α表示产状方向和水平面之间的夹角。

若产状方向近似为水平方向，则选择水平投影，若产状方向近似为垂直方向，则计算经过走向方向的垂直面，将该垂直面设定为投影平面，做垂直投影，图2的α接近水平方向，所以投影平面为水平平面。

2.2选择合适的投影平面，将三维矿体投影到水平或垂直平面上，形成矿体投影边界从而转到二维图，见图3。

将上一步得到的投影平面记为α，其解析式为Ax+By+Cz+D=0，记矿体中任意一点坐标为p₀=（x₀，y₀，z₀），则p₀在投影面α的投影坐标q₀：

$\left\{\begin{array}{c}{q}_0=(x_0-\mathrm{At},y_0-\mathrm{Bt},z_0-\mathrm{Ct})\\ t=\frac{{\mathrm{Ax}}_0+{\mathrm{By}}_0+{\mathrm{Cz}}_0+D}{A^2+B^2+C^2}\end{array}\right.\}---\left(3\right)$

三维矿体投影到平面上，形成二维多边形，该二维多边形称为矿体投影边界。公式（3）是计算矿体投影边界的基础。按照上述公式，将三角形△的三边投影到平面α上，会得到新的三角形△’。特殊情况下，当Δ⊥α时，△’会退化为一条线段。退化三角形并不影响矿体投影边界的计算。将三维矿体包含的三角形总数记作N，则将所有三角形的边都投影后，在平面α上将得到3·N条二维线段。将这些线段通过结点平差、拓扑重建两个子过程，即可计算得到三维矿体在平面α上的矿体投影边界P。

2.2.1结点平差

经过投影后的线段，是相互离散的。结点平差的目的，是将线段相互关联，对具有公共点的线段分组，为拓扑重建子过程做准备。在数据量不多的情况下，该过程可以将线段端点做两两比较，若相等，则对端点编号，并添加到编号对应的集合中。该方法简单易行，但时间复杂度为O（N²），不适合当数据量大的情况。当数据量大时，可采用将线段端点排序后再比较的方式。端点排序的时间复杂度为O（N·logN），端点比较的时间复杂度为O(N)，总的时间复杂度为O（N·logN）。

2.2.2拓扑重建

拓扑重建是为了找出这些线段的最外侧线段，并将它们首尾连接，形成投影多边形。多边形的搜索方法有多种，常用的有左转算法、右转算法等，本发明采用的是左转算法，它具有计算过程清晰、循环次数恒定等优点。设线段的数目为N，左转算法的具体步骤如下：

（1）计算每个端点相邻接的线段形成的方向角；并按照方向角升序排序；计数器c清零；

（2）计数器c加1，若c大于2·N，转步骤（7）；选择第c条线段(c≤N)或者第c-N条线段(N＜c≤2·N)作为起始线段R_start；

（3）将R_start赋值给R；

（4）若c≤N，查找R尾端点的下一条线段；若N＜c≤2·N，查找R首端点的下一条线段；即：R的后继线段R_next。若R_next等于R_start，则转（5），否则，转（6）；（5）将本次遍历得到的线段形成多边形，若该多边形为逆时针的，则将该多边形输出，转（2）；

（6）将R_next赋值给R，转（4）；

（7）算法结束。

需要说明的是，实际中的三维矿体往往是由多个三维实体构成，因此，矿体投影边界P并不一定是一个多边形，也基于这个原因，在拓扑重建算法的第（5）步，不能在找到一个逆时针的多边形的时候就终止算法。

3）计算投影后的钻孔位置、平均品位、见矿真厚度等信息

以勘探剖面和矿体投影边界为约束，对投影后的钻孔按照一定的地质原则进行划分，形成合理的四边形块段，如图4所示。在每个同心圆旁边都标注有ZKi表示钻孔的名称，钻孔名称旁边的数据分别依次表示钻孔的品位和见矿真厚度。

3.1投影钻孔

以钻孔ZKi为例，从井口位置开始查找ZKi的首次见矿坐标k，使用上述公式，将该坐标投影到平面α，记作k’。根据2.2节的说明得到矿体投影边界P，若k’落到矿体投影边界P中，则k’记为有效钻孔，并添加到有效钻孔集合中。以此处理所有的钻孔，得到有效钻孔集合，参与块段划分。

3.2划分块段

地质块段的划分，是根据矿石的不同工业类型、不同品级、不同矿产资源储量级别等地质特征进行的。在本发明中，依据的地质信息除了钻孔投影位置，还引入了剖面约束。即：计算得到的块段必须沿着勘探剖面排布，并且块段不能跨越勘探剖面。通过读取钻孔所在的剖面信息，按照剖面的空间顺序，对每两个剖面上的钻孔分别进行块段划分，从而得到矿体范围的全体块段。令剖面数为S,投影边界记作P，划分块段的具体算法如下：读取剖面数据，将剖面按照从左到右（或从右到左）的顺序排列。对每个剖面的钻孔，按照左到右（或从右到左）的顺序排列。

（1）读取第i个和第i+1个剖面，以剖面为约束，进行进行Delaunay三角剖分，生成带勘探剖面约束的三角网。其中，1≤i≤S-1。遍历完所有剖面，生成初始三角网，完成钻孔间的初始块段划分。

（2）计算三角网的边界，并计算边界点到矿体投影边界P的最短线段，利用这些最短线段和三角网的边界，完成矿体投影边界P的初始块段划分。

（3）对钻孔间和矿体投影边界上的初始块段划分，使用人机交互的手段，删除多余的线段，将三角块段转换成四边形块段，形成最终的块段划分。

算法第（1）（2）两步，是为了将钻孔按照一定次序排列，从而易于划分块段。这样的处理方式，使得块段划分只局限于剖面之间，既能保证块段划分的正确，也提高了算法速度。算法的第（3）（4）两步，分别得到钻孔间的初始块段和矿体投影边界P上的初始块段。算法第（5）步，引入人机交互手段，使得块段划分更加适合具体的实际要求，该过程完成后形成的块段，大多数是四边形的。同时，在边界钻孔或矿体投影的边界处，有可能形成三角形或者多边形的块段。4）估算每个钻孔的平均品位、见矿总真厚度，计算块段在投影面上的面积、平均真厚度、块段体积、品均品位，估算每个块段的储量。

块段内储量估算的基本依据是：M＝ρ·V·C （4）

其中，M是块段的储量估算值，ρ是平均密度，V是块段体积，C是块段的平均品位。块段平均密度ρ可以通过其他地质资料获得。这里的重点是块段体积V、块段平均品位C的估算。

4.1块段体积V的估算

由于块段已经投影到平面α，其在平面α上的面积容易求得：设多边形T的坐标为(x_i,y_i)，1≤i≤n，为简化公式的表述，定义(x_n+1,y_n+1)＝(x₁,y₁)，则多边形T的面积为：

$S_T=\left|\begin{array}{c}\frac{1}{2}\end{array},\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\right|\left|\begin{array}{cc}{x}_i& y_i\\ x_{i+1}& y_{i+1}\end{array}\right|---\left(5\right)$

其中，为2×2的行列式，其值等于x_i·y_i+1-y_i·x_i+1。

对于有效钻孔集合中的钻孔ZK_i，从钻孔的井口位置到井深位置，统计高于工业品位的见矿总真厚度对于每个块段，设块段包含的有效钻孔数目为m，计算有效钻孔的见矿总真厚度的平均将作为块段真厚度，即：

${\stackrel{\‾}{h}}_{{\mathrm{ZK}}_i}=\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{h}_{{\mathrm{ZK}}_j}---\left(6\right)$

利用基本体积公式V＝S·h （7）

将公式（5）（6）带入公式（7）中，估算块段的体积。块段的体积估算值为：

$V=\mathrm{abs}\left(\begin{array}{c}\frac{1}{2}\end{array}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left|\begin{array}{cc}{x}_i& y_i\\ x_{i+1}& y_{i+1}\end{array}\right|\right)\·\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{h}_{Z{K}_j}---\left(8\right)$

4.2块段平均品位C的估算

将有效钻孔上岩性段的品位以见矿真厚度做加权平均，可以求得钻孔的平均品位。对于有效钻孔ZK_i，从钻孔的井口位置到井深位置进行统计，如果岩性段的品位高于工业品位，则记为(c_i,h_i)，其中，c_i为岩性段的品位，h_i为该岩性段的真厚度。设得到了k个这样的“品位-真厚度”二元组，将岩性段的品位以真厚度做加权平均，作为钻孔的平均品位，即：

$c_{{\mathrm{ZK}}_{\mathrm{ii}}}=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{c}_i\·h_i}{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{h}_i}---\left(9\right)$

块段的品位，可以根据块段包含的有效钻孔的平均品位，以钻孔的见矿总真厚度作为加权平均得到，将块段包含的有效钻孔数目记为m，则块段的品位为：

$C=\frac{\underset{{\mathrm{ZK}}_i=1}{\overset{m}{\mathrm{\Σ}}}{h}_{{\mathrm{ZK}}_i}\·c_{{\mathrm{ZK}}_i}}{\underset{{\mathrm{ZK}}_i}{\overset{m}{\mathrm{\Σ}}}{h}_{{\mathrm{ZK}}_i}}---\left(10\right)$

将公式（8）（10）代入公式（4）中，即求得块段内的储量估算值：

$\mathrm{\ρ}\·\mathrm{abs}\left(\begin{array}{c}\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left|\begin{array}{cc}{x}_i& y_i\\ x_{i+1}& y_{i+1}\end{array}\right|\end{array}\right)\·\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{h}_{Z{K}_j}\·\frac{\underset{{\mathrm{ZK}}_i=1}{\overset{m}{\mathrm{\Σ}}}{h}_{Z{K}_i}\·c_{Z{K}_i}}{\underset{Z{K}_i=1}{\overset{m}{\mathrm{\Σ}}}{h}_{Z{K}_i}}---\left(11\right)$

5）根据块段的投影面积、平均真厚度，平均品位等信息，并结合钻孔数据中提供的伴生有益组份、有效钻孔距离等信息，对块段进行储量分级（分级方法），将同样级别的块段储量相加，估算矿体的储量，并将上述计算结果按照统一标准输出为储量报告。

最后应说明的是，以上实施例仅用以描述本发明的技术方案而不是对本技术方法进行限制，本发明在应用上可以延伸为其他的修改、变化、应用。

Claims (8)

1.一种地质勘探三维可视化储量估算方法，包括：

a.计算矿体的产状方向、走向和倾向，将矿体投影到水平或垂直面上；

b.根据投影后的钻孔位置、平均品位和见矿真厚度，将投影后的钻孔进行块段的划分，形成四边形块段；

c.通过每个钻孔的平均品位、见矿真厚度和四边形块段的面积，计算每个块段对应的矿体的储量；

d.对矿体储量进行分级，估算矿体的总储量；

其特征在于，在所述的步骤a和b之间，还包括建立三维矿体在水平或垂直面上的投影边界，投影边界的建立方法如下：

2.1、设三维矿体的投影平面为α，解析式为Ax+By+Cz+D＝0，其中A、B、C表示投影平面的法线方向，D确定投影平面的空间位置；

2.2、矿体中任一点的坐标为p₀＝(x₀,y₀,z₀),在投影面上的坐标为q₀为：

2.3、当三维矿体投影到平面上，形成二维多边形，该二维多边形为投影边界；

2.4、通过结点平差和拓扑重建方法建立三维矿体的投影边界。

2.根据权利要求1所述的方法，所述的步骤a包括以下步骤：

1.1、首先使用三维建模方法，完成三维钻孔、三维剖面以及三维矿体数据的建模，其中，三维矿体数据是由一组相互邻接的三角形拼接形成的；

1.2、根据三维矿体中各个三角形的顶点相对于矿体重心的偏移向量，计算出整个矿体的产状方向、走向和倾向；

1.3、若矿体的产状方向近似为水平方向，则将矿体投影到水平面，如果矿体的产状方向近似为垂直方向，则将矿体投影到垂直面。

3.根据权利要求2所述的方法，在所述的步骤1.2)中，矿体的重心μ通过以下方法计算得出：将三角形的顶点分别记为Pⁱ、qⁱ和rⁱ， 式中，N表示矿体表面包含N个三角形；第i个三角形的顶点相对于矿体重心的偏移向量分别为：

整个矿体的产状方向、走向和倾向的计算方法如下：

都是3×1的向量，即计算得到的C_jk是3×3协方差矩阵，该协方差矩阵的三个特征向量分别表示矿体的产状方向、走向和倾向。

4.根据权利要求1所述的方法，其特征在于，将矿体上的钻孔投影到投影投影面上，如果钻孔的投影在投影边界外，则为无效钻孔；如果钻孔的投影在投影边界内，则为有效钻孔。

5.根据权利要求1所述的方法，其特征在于，在所述的步骤b包括以下步骤：

3.1、首先将矿体划分为S个相互平行的剖面；

3.2、根据每个有效钻孔所在的剖面信息，对每两个相邻剖面上的有效钻孔进行块段的划分，从而得到矿体范围内的全部块段。

6.根据权利要求5所述的方法，其特征在于，在所述的步骤3.2包括以下步骤：

3.2.1、读取剖面数据，将剖面按照从左到右的顺序排列；

3.2.2、对每个剖面的有效钻孔，按照从左到右的顺序排列；

3.2.3、读取第i和第i+1个剖面，对剖面进行三角剖分，1≤i≤S-1，所有剖面都依次处理，完成钻孔间初始块段的划分；

3.2.4、计算三角网的边界，并计算三角网边界点到矿体投影边界P的最短线段，利用这些最短线段和三角网边界完成对矿体投影边界P上的初始块段划分；

3.2.5、对钻孔间和矿体投影边界的初始块段划分，使用人机交互的手段，删除多余对象线，将三角形块段转换成四边形块段，形成最终的块段划分。

7.根据权利要求1所述的方法，其特征在于，在所述的步骤c中所述的每个块段对应的矿体储量通过下式计算：

M＝ρ·V·C

其中，ρ表示块段平均密度，V表示块段对应的矿体体积，C表示块段对应的矿体平均品位。

8.根据权利要求7所述的方法，其特征在于，所述的块段对应的矿体储量通过以下方法计算：

其中，x_i表示投影面上的块段的在x轴上的一个坐标值，X_i+1表示投影面上的块段的在x轴上相邻的一个坐标值；y_i表示投影面上的块段的在y轴上的一个坐标值，y_i+1表示投影面上的块段的在y轴上相邻的一个坐标值，n表示块段的数量；h_{zk j}表示每个有效钻孔见矿的真厚度，h_zki表示第i个钻孔的真厚度，C_zki表示第i个钻孔的岩性段的品位，m表示每个块段的有效钻孔数目，ρ表示块段的平均密度；M表示块段内的储量估算值。