CN108595749A - 一种使用变异函数单一方向结构分析的资源储量评估方法 - Google Patents

Abstract

本发明涉及使用变异函数单一方向结构分析的资源储量估算方法，该方法包括如下步骤：对矿段做等长组合采样，进行信息统计，得出样品数据的品位值的分布范围，均值、方差、偏度；针对品位值，构建分形模型，获得各组分的概率分布及组分个数，所述组分是指在分形范围内拟合出一条或多条线段，每一个线段代表了一个组分；确定块金方差及部分基台值；利用组分个数及部分基台值信息，结合实验变异函数曲线特征，建立各组分变程与其部分基台值的联系；建立变异函数结构分析与地质意义的联系，以及对资源储量估算的参数约束，进一步进行资源储量评估。

Description

一种使用变异函数单一方向结构分析的资源储量评估方法

技术领域

本发明涉及地质资源勘探领域，特别涉及一种使用变异函数单一方向结构分析的资源储量估算方法。

背景技术

地质资源储量是矿山设计和开发的基础。采用先进的资源储量估算方法来提高估算的精度，对于降低矿山设计风险，增加矿山开发的盈利能力，促进矿山可持续发展有着重要意义。近半个世纪以来，资源储量估算从传统的块段法、平行剖面法、多边形法、距离幂次反比法为主，逐渐发展成为现代主流的地质统计学方法(克里格法)。地质统计学估值具有最优无偏估计特征，它考虑了地质变量的空间关系，具有比传统方法估算更精确、可以给出估算精度、能够进行整体和局部估计等诸多优点(侯景儒和黄竞先，1982；王仁铎和胡光道，1987)。尤其是随着近年来计算机科学的发展，地质统计学方法也在不断的发展完善，逐渐成为国际上进行矿产资源储量估算的指定方法。我国从上世纪90年代开始，也开始引进地质统计学方法，并将其作为资源储量评审备案认可的估算方法之一。然而，国内应用地质统计学方法提交资源储量报告的案例还很少，尚处于实践摸索阶段。

地质统计学以变异函数为基本工具。变异函数既能描述区域化变量的结构性变化，又能描述其随机性变化，是地质统计学估算的基础(刘爱利等，2015)。以普通克里格法为例，其估计量为：

各样品点的权系数λ_i通过普通克里格方程组来求解：

普通克里格估计方差为：

式中，γ(x_i,x_j)、均为与空间矿化结构有关的变异函数γ(h)。可以看出，变异函数在资源储量估算的克里格权重系数和估计方差计算中，都起着基础性、决定性的作用。而且，变异函数提供了矿化空间变异性的程度(基台值)和范围(变程)信息，为后续资源储量估算参数的确定提供了重要依据。因此，优化变异函数结构分析是提高资源储量估算准确度和精度的有效途径。

变异函数通过研究相隔距离为h(h为向量)的两个区域化变量Z(x)和Z(x+h)数值的离差平方的期望，来研究任意两点Z(x)和Z(x+h)的空间相关性，也就是结构性。其一般表达式为：

区域化变量的随机性用变异函数在原点处的性状——块金方差(γ₀)来表示。根据变异函数的性质有，

γ(h)＝γ₀+γ₁(h)+γ₂(h)+γ₃(h)….+γ_i(h)，i＝1,2,3……

式中，γ_i(h)表征变异函数单一方向不同尺度或者不同方向不同尺度的结构(structure)或者组分(components)。变异函数这种由多个结构或者组分叠加的形式叫做套合结构(nest structure)。套合结构中，每一个变异函数代表一种特定尺度上的变异性。

变异函数结构分析的一般流程是：选取区域化变量、数据审计、基本统计分析、计算实验变异函数、运用理论变异函数进行结构最优拟合、模型检验以及结构的地质解译。其中，结构分析的最核心内容是通过实验变异函数的最优拟合，来确定资源储量估算的块金方差、部分基台值以及变程等重要参数。一般的做法是先进行单一方向不同尺度的结构分析，然后运用线性变换等方法进行不同方向变异函数的套合。

目前，单一方向变异函数结构分析的主要内容是(多)组分识别、理论变异函数拟合和地质解译与成果应用。

(1)单一方向结构分析的基本原理

区域化变量的每一个采样点数据可以当做它的一个实现。利用已知空间采样数据可以计算出实验变异函数。它的数学表达式如下：

其中，h为滞后距(lag)，Z(x_i)和Z(x_i+h)为区域化变量在x_i和x_i+h点处的值，二者相隔h，N(h)为数据对个数，γ*(h)为计算的实验变异函数值(experimental semivariogram)。

在某一方向上，不断的改变h，就可以得出一系列h，N(h)和γ*(h)值(或称Gamma*值)。以滞后距离h为横轴，γ(h)为纵轴，以散点图的形式把变异函数表现出来，如图1所示。

实验变异函数的拟合可以用不同的数学模型，包括球状模型、高斯模型、指数模型等，最常用的是球状模型。利用数学模型对实验变异函数进行拟合，可以得到理论变异函数曲线及结构模型。主要数学模型的公式为：

球状模型

高斯模型

指数模型

例如，可以利用球状模型对实验变异函数曲线进行拟合(图1)。

其中，C₀表示块金方差(nugget)或块金值，C₁为部分基台值(partial sill)，C为基台值(sill)，有C＝C₀+C₁。h为滞后距或步长(lag)，a为变程(range)，超过这一尺度，认为区域化变量不具有空间相关性或结构性，也即h＝a时，γ(h)趋近于基台值。

以上是单一结构组分的情形。单一方向上，如果存在不同尺度的套合结构，则理论变异函数可以写成单个变异函数的加和形式。套合结构的每个组分可以是不同数学模型的变异函数。以2组分的球状模型为例进行说明。

其块金方差模型为：

第一个组分的球状模型为：

第二个组分的球状模型为：

那么，总的套合结构为γ(h)＝γ₀+γ₁(h)+γ₂(h)，当a₁＜a₂时，γ(h)可以写成分段函数形式：

其2组分结构拟合结果示意图见图2。

(2)单一方向套合结构的识别

传统方法主要是通过观察实验变异函数曲线特征来识别结构组分。一般来讲，如果存在多个结构的套合，那么实验变异函数曲线就存在拐点。具体的，若发现靠近原点的前区曲线很陡，后区曲线较缓时，就可以考虑用套合结构模型来拟合(王仁铎和胡光道，1987)。实验变异函数曲线上明显的肩部或者曲线斜率变化处，可能指示了有多少个组分。具体包含几个组分，取决于研究区特征，但是要尽可能的用最少的组分来拟合实验变异函数。

(3)单一方向套合结构的拟合

用理论变异函数拟合实验变异函数时，如果存在多个组分，各个组分之间是存在相互作用的。改变大尺度组分的部分基台值，总基台值一定的情况下，小尺度组分的基台值也相应发生改变。因此，拟合时需要关注组分个数、部分基台值和相应变程等参数变化。常用的拟合方法主要有经验法和最小二乘法。

①经验法。

首先模拟块金方差，根据井向方向实验变异函数头2-3个点进行拟合，拟合曲线与纵轴的交点即为块金方差值。然后根据拐点特征，组分从少向多，变程由小到大进行逐个拟合。由于理论模型的总基台值为块金方差与各组分的部分基台值之和(C＝C₀+C₁+C₂+…)，那么计算某一组分的部分基台值，可以用总基台值减去块金方差以及之前所有组分的部分基台值。通过反复调整参数，直至理论变异函数曲线能够反映出实验变异函数的趋势特征。

②最小二乘法。

块金方差用一元一次线性拟合来实现，这是一种计算机自动拟合的方法。令y＝γ*(h)，x＝h，则有y＝ax+b，b值即为块金方差C₀。

套合结构的拟合。以2组分球状模型的套合结构为例，根据拐点信息，将全部选取的数据点分为前区和后区两部分，前区和后区的分界点选在曲线由陡变缓的转折点处为宜。分界的数据点，既作为前区的最后一个数据点，也作为后区的第一个数据点。前区的数据点主要用来拟合第一段二级套合球状模型理论变异函数图，后区的数据点则用来拟合第二段二级套合球状模型变异函数图。其套合结构数学模型为：

设γ(h)＝y，h＝x₁，h₃＝x₂，C₀＝b₀

则有，y＝b₀+b₁x₁+b₂x₂

因此，实验变异函数的拟合就成了二元线性回归问题了，需要根据数理统计学中二元线性回归的理论和方法进行求解(王仁铎和胡光道，1987)。3组分以上套合结构的拟合方法也类似。

当然，也有大量学者认为采用计算机自动拟合是认识上的错误，是不可取的，应该用人工拟合(尹镇南，2012)。申请者也持这种观点，因为每个矿床的地质规律不一样，每个变异函数点的地质意义也不一样，有些点应该参与变异函数拟合，而有些点则应该果断的舍去。

传统方法进行单一方向的变异函数结构分析主要存在以下不足之处。

(1)结构识别：主要通过变异函数图进行识别。存在2种极端情形，如果曲线没有明显的肩部或者拐点，组分将难以通过曲线特征进行肉眼识别；或者存在多个拐点(参见图6)，识别的组分太多，容易造成对函数的过度拟合。因此，传统方法在变异函数结构识别上具有很强的人为性。

(2)结构拟合：以2组分套合结构为例，经验法由于存在5个未知参数(块金方差、2组分的部分基台值和变程)，在确定了组分个数和块金方差的基础上，反复调整部分基台值和变程，拟合是一个比较耗时的过程。最小二乘法同样也依赖于组分的识别，组分分界点的识别不合理，最小二乘法的拟合结果差别会很大；另外，实验变异函数不同于理论模型，它或多或少存在噪音，这些结构中的扰动可能会影响拟合结果。不论用哪种方法，都要求实验变异函数曲线上存在明显的拐点或肩部，但经验表明这一要求常常难以满足。

(3)地质解释：基于实验变异函数形态的结构分析，可以定量地确定存在于不同尺度的多个组分，而组分的地质意义会具有比较强的多解性。因此，传统变异函数结构分析的结论应用具有一定局限性。

因此亟需一种利用变异函数结构分析的资源储量估算方法，具有更高的精度和速度，以达到优化资源储量估算的目的。

发明内容

本发明的目的是针对现有技术存在的问题，提供一种使用变异函数单一方向结构分析的资源储量估算方法，克服传统的变异函数结构分析方法具有结构识别存在主观性、拟合过程繁琐、地质解译多解性等不足。

上述目的是通过下述方案实现的：

一种使用变异函数单一方向结构分析的资源储量评估方法，其特征在于包括如下步骤：

a.对矿段做等长组合采样，进行信息统计，得出样品数据的品位值r的分布范围，均值m、方差S²、偏度k；

b.针对品位值r，构建分形模型，获得各组分的概率分布P_i及组分个数，所述组分是指在分形范围内拟合出一条或多条线段，每一个线段代表了一个组分；

c.确定块金方差C₀及部分基台值C_i，C_i＝P_i(C-C₀)，其中C₀为块金方差，C为总基台值，i为组分数；

d.利用组分个数及部分基台值C_i信息，结合实验变异函数曲线特征，建立各组分变程a_i与其部分基台值的联系；

e.建立变异函数结构分析与地质意义的联系，以及对资源储量估算的约束进一步进行资源储量评估。

根据上述评估方法，其特征在于，所述步骤b中，针对不同的品位值r，统计大于r的数据个数N(r)，据此计算logr和logN(r)，在双对数坐标轴上，做logr～logN(r)散点图；用最小二乘法进行拟合，拟合优度R²≥0.95，拟合得到一条或多条线段，每一条线段代表一个组分；利用该线段两端点的纵坐标的反对数求差，得到该组分中的样品点个数N_i(r)，据此计算各组分的概率，

P_i＝N_i/N

其中，N为分形范围内总样品数，N_i为某组分的样品个数。

根据上述评估方法，其特征在于，所述步骤c中，改变滞后距h大小，计算实验变异函数γ*(h)值，做h-γ*(h)散点图；利用实验变异函数的头2-3个点拟合块金方差C₀；根据公式C_i＝P_i(C-C₀)，计算出各组分的部分基台值；在满足采样空间范围大于3倍变程的条件下，用样品方差估计基台值，可得C_i＝P_i(S²-C₀)。

根据上述评估方法，其特征在于，所述步骤d中，先初步确定各组分变程，再通过变程微小的调整，使实验变异函数的拟合达到最优。

根据上述评估方法，其特征在于，所述步骤e中，根据所确定的理论变异函数的C₀、C_i和a_i结构参数，确定变异函数套合结构的数学模型；根据该模型对资源储量估算的条件和参数进行约束。

根据上述评估方法，其特征在于，所述约束为确定矿化背景值和矿化结构，确定资源储量估算邻域。

根据上述评估方法，其特征在于，块金效应越小，屏蔽作用越强，估计邻域范围可以缩小；相反若块金效应大，屏蔽作用弱，估计邻域范围可以扩大。

本发明通过区域化变量的分形模型，可以定量给出变异函数结构的组分个数及其部分基台值等拟合关键参数，不仅仅依赖于实验变异函数曲线特征，减少了拟合的人为性；

本发明结合组分部分基台值和实验变异函数曲线拐点和肩部信息，可以大致确定变程，实现对曲线的快速拟合；

本发明有着明确的物理意义，变异函数结构的地质意义与构造分形模型时的区域化变量具体属性有关，同时对资源储量估算的条件和参数有明显的约束。

附图说明

图1单一方向实验变异函数曲线及其拟合示意图；

图2单一方向变异函数套合结构的拟合示意图；

图3使用单一方向变异函数结构分析的资源储量评估流程图；

图4某铜钴矿床等长组合样品数据频率直方图；

图5某铜钴矿床N-S分形模型拟合图；

图6某铜钴矿床实验变异函数曲线特征及块金方差的拟合；

图7某铜钴矿床各组分变程的拟合。

具体实施方式

下面结合本发明说明书实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述。

在本实施方式中，以钻探和槽探工程的采样分析数据为区域化变量，基于Micromine软件平台，探讨具体的操作流程以及实施案例。参阅图3的使用单一方向变异函数结构分析的资源储量评估流程图。

1.基本统计分析。

在本实施方式中，用等长组合样品数据做频率直方图，并进行基本信息统计，主要包括样品数据的分布范围(min和max)、均值(m)、方差(S²)、偏度(k)等特征。

2.运用分形统计方法获得各组分的概率分布P_i。

分形方法是近30多年来被广泛应用于地质地球化学数据处理的一种新方法。大量实践证明，矿化元素在特定的标度区间内具有分形分布特征，指示了特定矿化作用或矿化类型，并且利用分形方法可以有效区分这种多期次或不同类型的矿化叠加作用(Blenkinsop,1991；Cheng,1996；Afzal etal.,2013；Heidari et al.,2013)。

分形模型的构建。根据Number-Size(N-S)模型(Turcotte,1992)，若样品品位数据满足分形分布则有，N(≥r)＝Cr-D，其中r为品位值，N(r)为品位大于r的样品个数，C为常数，D为分维值。对于给定不同的品位r，统计大于或等于r的品位数据个数N(r)，在双对数坐标轴上，做logr～logN(r)散点图。用最小二乘法拟合(拟合优度≥0.950)，拟合直线斜率的绝对值就是分形维数D(Heidari et al.，2013)。并不是所有数据点都符合分形分布特征(Blenkinsop,1991)，把具有分形分布特征的品位区间称作分形标度区间。如果品位分布只在一个特定标度区间内具有分形分布特征，叫做单标度分形；如果在多个标度区间内进行拟合，叫做多标度分形，对应多个分维值。所有分形标度区间加起来的变化范围称为分形范围(Fractal range)。分形范围之外的数据点属于特异值，应该在做变异函数分析前处理掉。在分形范围内，可能拟合出一条或多条线段，每一个线段代表了一个组分。每一组分中样品点个数可以用该线段端点的纵坐标的反对数求差得来，其分布概率(累计频率)P_i为该组分的样品点个数(累计频数)除以分形范围内总样品数。

P_i＝N_i/N

其中，N为分形范围内总样品数，N_i为某组分的样品个数，P_i为某组分的累计频率，当样品数据足够多的情况下，可以作为该组分的累计概率。

在本实施方式中，构建N-S分形模型，对于不同的品位r，计算大于r的数据个数N(r)。在双对数坐标轴上，做logr～logN(r)散点图。用最小二乘法进行拟合，形成连续不同斜率的线段，拟合优度R²＝0.95。拟合时，尽量将拟合的线段条数降到最低，这样得到最少的组分数。通过线段端点信息，计算出每一组分的数据频数(N_i)，当数据量足够大时，则各组分的概率可以用频率来表示，即P_i＝N_i/N。

3.利用组分个数及其部分基台值(C_i)信息，结合实验变异函数曲线特征，建立各组分变程(a_i)与其部分基台值(C_i)的联系。

实验变异函数的计算一般在软件中进行。以不同的h，计算出γ*(h)值，做h-γ*(h)散点图。

由于块金方差反映的是系统误差(如测量误差)或者小于最小取样样长的微观结构，因此，可以利用井向变异函数来拟合块金方差。利用实验变异函数的头2-3个点，很容易通过手动方式，找到理论变异函数曲线与纵轴的交点，即为块金方差值C₀。

再确定各个组分的部分基台值，利用公式C_i＝P_i(C-C₀)，计算出各组分的部分基台值。

各个组分在理论变异函数曲线中的位置，与组分的变程有关。一般地，变程越小的组分，越位于理论变异函数曲线的底部，变程越大的组分，越位于曲线的顶部，但也存在很多特殊的情形。为了便于说明，我们讨论最一般的情况，就是品位越高的组分，其连续性越差，具有较小的变程，位于理论变异函数曲线的底部；品位越低的组分，连续性相对较好，具有较大的变程，位于理论变异函数曲线的顶部。对于特殊的情形，因为组分的个数和部分基台值已经确定，也很容易从曲线上找到各自的位置。

如果实验变异函数曲线存在肩部或者拐点，那么变程的拟合就很简单。从图6上可以看出，曲线中明显的肩部或者拐点的横坐标直接指示了不同组分的变程。但是它们的纵坐标并不能直接反映各组分的部分基台值(图2)，这是因为较小滞后距的情况下，拐点处的实验变异函数值为多个组分变异函数值加和的结果。因此，变程较小组分的部分基台值要比拐点对应实验变异函数值低。也正是利用这一关系，可以定性确定各组分在理论变异函数曲线中的具体位置。通过上述方法，先初步确定各组分变程，再通过变程微小的调整，来使实验变异函数的拟合达到最优。

在本实施方式中，改变滞后距大小，计算实验变异函数值γ*(h)，做出h-γ*(h)散点图，连接各点，形成实验变异函数曲线。先用头2-3个实验变异函数点拟合块金方差C₀，然后根据各组分概率计算各组分的部分基台值C_i＝P_i(C-C₀)。当采样空间范围都大于3倍该方向变程，可以用样品方差估计基台值，则有C_i＝P_i(S²-C₀)，当采样空间范围小于3倍的变程，基台值直接从实验变异函数图中读取。

4.建立结构分析与地质意义的联系，以及对资源储量估算的约束。

确定了理论变异函数的C₀、C_i和a_i等结构参数后，可以确定变异函数套合结构的数学模型。变异函数结构分析的最后一步是地质解释，也就是说明变异函数结构分析的实际物理意义。

N-S分形模型的不同组分代表了不同的品位域。不同组分有着不同的变程，指示了不同品位区间样品空间上分布的自相关性或连续性差异，变程越大，该品位域的空间连续性越好，反之连续性越差。因而，变异函数结构的地质意义与构造分形模型时的区域化变量具体属性有关。推而广之，需要了解什么样的地质属性，就将这种属性作为区域化变量，来构建特定的分形模型，这样通过变异函数的结构分析，可以很好的揭示这种地质属性的趋势特征。

更为重要的是，基于上述变异函数的结构分析，可以为资源储量估算提供更多的参数约束，以提高估算精度。

(1)对γ(x_i,x_j)的约束

γ(x_i,x_j)是克里格方程组中的重要条件，其物理意义为估计邻域内已知样品点对的变异函数。γ(x_i,x_j)取决于哪些已知样品点要参与资源储量估算。资源储量估算的品位域应遵循变异函数的结构特征，不破坏矿化空间结构，换言之要与区域化变量的地质统计学特征一致。采用本申请的变异函数结构分析方法，可以很容易确定哪些是指示了矿化背景值，哪些是资源储量估算真正关心的矿化结构。

(2)对的约束

是克里格方程组和估计方差的计算内容，其物理意义为样品点到待估矿块之间的所有点对的平均变异函数值。显然与已知样品点到待估矿块之间的距离有关。变异函数的结构分析显示不同的结构有着不同的变程，资源储量估算应采用矿化域所对应结构的变程信息。因此，在确定资源储量估算邻域(样品搜索半径)时，应该充分考虑结构分析的成果。此外，还应考虑块金方差对估计邻域的影响。一般而言，块金效应(块金方差/基台值)越小，屏蔽作用越强，估计邻域范围可以适当缩小；相反若块金效应大，屏蔽作用弱，估计邻域范围可以适当扩大。

在本实施方式中，结合实验变异函数曲线的肩部和拐点信息，先确定各个组分大致的变程(a_i)，然后根据实际情况进行变程的微调，达到实验变异函数的最佳拟合。

在本实施方式中，在确定块金方差、各组分部分基台值及变程后，即可将套合结构的理论变异函数写成各组分加和的形式，完成最终的拟合。由于理论变异函数不同品位域结构有着各自空间连续性或者地质统计学规律，结合特定矿床的经济边界品位，可以很容易确定资源储量估算的边界。设置资源储量估算的估计邻域时，应充分考虑矿化结构的变程，同时结合块金效应的屏蔽作用对估计邻域的影响。

在一个场景示例中，研究区位于中非铜矿带刚果(金)矿段某铜钴矿床。矿体产出受岩性控制，经历了热水沉积、热液改造和表生氧化等多期矿化作用，为层控型矿床。本案例选取该矿床的其中一个矿段进行研究。矿体呈EW向展布，向南倾斜，倾角50-70°，矿体长约350m，宽约200m，厚度约100m。所研究的矿段由地表8个浅井、8个探槽和11个钻孔控制，勘查工程间距为50×50m。根据矿化规律圈定了矿化体，矿化体中探槽和浅井采样基本样长为1m，共823个数据；钻孔采样样长整体在1-2m之间共785个数据。为了消除样品权重的影响，用1m等长组合样后(共1886个数据)进行变异函数分析。

(1)基本统计分析

用1m等长组合样品数据做频率直方图(图4)，并进行基本信息统计，主要包括样品数据的分布范围0.01-22.17％，样品均值(m)为1.69％，样品方差(S²)为7.26％²。

(2)分形统计分析

构建N-S分形模型，对于不同的品位r＝0.01％，0.02％，...，19.9％，统计大于r的数据个数N(r)，据此计算logr和logN(r)，见表1。在双对数坐标轴上，做logr～logN(r)散点图。用最小二乘法进行拟合，形成连续不同斜率的线段，拟合优度R²≥0.95，见图5。拟合共得到3个组分AB、BC和CD，分别对应于低品位域(0.01-0.3％)、中品位域(0.3-6％)和高品位域(6-22.17％)，不存在特异值，可以直接计算各组分概率。通过线段端点B和C信息，计算出各组分的数据个数，据此计算各组分的概率(表2)，分别为0.41、0.58和0.8。

表1 N-S分形模型计算

表2各组分数据个数及概率

组分	品位区间	数据个数(Ni)	概率(Pi)
				低品位域	0.01-0.3％	N＝1886-1112＝774	P＝774/1886＝0.41
中品位域	0.3-6％	N＝1112-153＝959	P＝959/1886＝0.51
				高品位域	6-22.17％	N＝153	P＝153/1886＝0.08

(3)块金方差拟合和各组分部分基台值计算

在Micromine中，改变滞后距大小，计算实验变异函数值γ*(h)，见表3。做出h-γ*(h)散点图，连接各点，形成实验变异函数曲线。

表3实验变异函数计算结果

h	N(h)	γ*(h)	h	N(h)	γ*(h)	h	N(h)	γ*(h)	h	N(h)	γ*(h)
												1	1861	1.267	20	1472	6.122	39	1143	8.079	58	820	7.536
2	1837	2.349	21	1455	6.240	40	1126	8.275	59	804	7.564
												3	1813	2.975	22	1438	6.327	41	1108	8.454	60	788	7.545
4	1789	3.408	23	1421	6.576	42	1090	8.455	61	772	7.455
												5	1765	3.692	24	1404	6.758	43	1072	8.423	62	756	7.261
6	1742	3.958	25	1386	6.841	44	1055	8.349	63	740	7.054
												7	1721	4.125	26	1369	7.014	45	1038	8.336	64	725	7.087
8	1700	4.203	27	1352	7.070	46	1021	8.319	65	710	6.998
												9	1680	4.613	28	1335	7.123	47	1004	8.230	66	695	6.966
10	1660	4.985	29	1318	7.243	48	987	8.150	67	681	6.854
												11	1640	5.140	30	1301	7.439	49	970	8.067	68	667	6.940
12	1621	5.235	31	1284	7.656	50	953	8.000	69	654	6.991
												13	1601	5.341	32	1265	7.857	51	936	7.807	70	641	7.024
14	1581	5.461	33	1247	7.902	52	919	7.671	71	628	6.992
												15	1561	5.627	34	1229	7.824	53	902	7.588	72	615	7.010
16	1542	5.675	35	1211	7.842	54	885	7.361	73	602	6.974
												17	1524	5.841	36	1194	7.925	55	868	7.520
18	1506	6.001	37	1177	8.012	56	852	7.573
												19	1489	6.088	38	1160	7.987	57	836	7.580

先用头2-3个实验变异函数点拟合块金方差C₀。Micromine中提供了一般线性拟合工具，利用头2个点进行拟合，拟合结果为C₀＝0(图6)。

然后根据公式C_i＝P_i(C-C₀)，计算各组分的部分基台值。从实验变异函数曲线图上看，γ*(h)有明显高出样品方差7.26％²的部分，且在γ*(h)＝8.0％²左右时，存在一个明显的肩部，部分学者会选择使用8.0％²作为理论变异函数的基台值C。但从实验变异函数曲线可以发现，井口方向上最大变程在30m左右，而钻孔的采样深度在100-150m左右，满足采样空间范围大于3倍变程的条件。因此，仍可以用样品方差估计基台值，也即C_i＝P_i(S²-C₀)＝P_i*(7.26-0)，结果见表4。

表4某铜钴矿床各组分的部分基台值计算

组分	品位区间	概率(Pi)	部分基台值(Ci)
				低品位域	0.01-0.3％	0.41	2.98％2
中品位域	0.3-6％	0.51	3.70％2
				高品位域	6-22.17％	0.08	0.58％2

(4)各组分变程拟合

本案例中实验变异函数曲线中有多个“可疑”的肩部和拐点(图6)。在手工拟合的情况下，需要结合各组分的部分基台值和曲线特点，逐个来确定各组分在拟合的理论变异函数中的位置。曲线上a点是一个明显的肩部，其对应的h值(变程)是7m左右，γ*(h)值在4.2％²左右，该值与中品位域的部分基台值3.7％²最为接近，据此可以大致第一个组分变程是7m，对应的部分基台值3.7％²。

然后再分别确定第二组分和第三组分的位置，只存在两种可能，分别尝试一下哪种拟合结果更优。结合曲线看，a-b段比b-f段曲线更为陡立，应该代表了不同的组分。计算发现a、b两点的γ*(h)值之差在0.8％²左右，略高于高品位域的部分基台值0.58％²，因此第二组分的变程应该是b点指示的10m，相应的部分基台值为0.58％²。剩下的就是第三组分，变程在30m左右，部分基台值为2.98％²。将以上参数输入软件，稍微调整变程，直至观察到的拟合结果最佳为止。实际上，在Micromine平台中，在输入已知的各组分部分基台值后，软件会自动给出了组分在理论变异函数中的位置。各组分在变异函数中的位置与输入数据时组分在软件中的具体位置无关。为了便于理解，我们把最终结果整理成了从第一组分到第三组分，组分变程越来越大的基本形式，而且选取了指数模型和高斯模型等不同的数学模型，以使拟合结果最优。

通过调整各组分变程，最终的拟合结果是，C₀＝0％²；C₁＝3.7％²，a₁＝6m；C₂＝0.58％²，a₂＝14m；C₃＝2.98％²，a₃＝33m。从图7可以看出，变异函数的小滞后距部分得到了很好的拟合，大滞后距部分理论变异函数曲线与实验变异函数曲线整体形态接近，因此拟合的整体效果较好。

(5)结构套合与地质解释

拟合的井向方向理论变异函数的最终套合结构形式为：

表5各组分的部分基台值计算

从变异函数结构分析(表5)可以看出，刚果(金)某铜矿床低品位域组分可能对应了矿化背景值，而中、高品位域则指示了矿化体。经矿山设计论证，该矿床铜的工业指标为边界品位＝1％，因此资源储量估算的品位范围可能是Cu≥1％，或者Cu≥0.3％，这取决于Cu＝1％的矿化边界性质。如果该边界为硬边界，资源储量估算利用的数据为根据Cu＝1％圈定的矿化体内的所有数据；如果Cu＝1％圈定矿化边界为软边界，也即矿化体内的数据点估计还需要利用边界之外的数据，这样资源储量估算的品位范围最佳范围为用Cu＝0.3％圈定的矿化体。因为，0.3％-6％范围内数据具有统一的空间自相关性，这样资源储量估算的范围确定充分考虑了地质统计学规律，也结合了矿床的经济开发价值。当然，我们还可以考察6％-22.17％范围内的数据的空间分布。如果空间上可以与0.3％-6％的中品位域数据分离，则可以单独圈定富矿体(富矿体的存在会大大提升矿山的盈利能力)，在确定边界性质后进行资源储量估算。

此外，在确定资源储量估算邻域时，要根据中、高品位域组分的变程来设置搜索半径。因此，井向方向上搜索半径设置在6-14m为最佳。反之，如果利用矿化背景值结构的变程33m作为搜索半径，可能会夸大矿化体的空间连续性，导致资源储量估算的偏差。同时，由于理论变异函数块金方差C₀＝0，说明数据的屏蔽作用很强，估计邻域的设置尽可能小一些。

以上所述的仅是本发明的较佳实施例，并不局限发明。应当指出对于本领域的普通技术人员来说，在本发明所提供的技术启示下，还可以做出其它等同改进，均可以实现本发明的目的，都应视为本发明的保护范围。

Claims (7)

1.一种使用变异函数单一方向结构分析的资源储量评估方法，其特征在于包括如下步骤：

a.对矿段做等长组合采样，进行信息统计，得出样品数据的品位值r的分布范围，均值m、方差S²、偏度k；

b.针对品位值r，构建分形模型，获得各组分的概率分布P_i及组分个数，所述组分是指在分形范围内拟合出一条或多条线段，每一个线段代表了一个组分；

c.确定块金方差C₀及部分基台值C_i，C_i＝P_i(C-C₀)，其中C₀为块金方差，C为总基台值，i为组分数；

d.利用组分个数及部分基台值C_i信息，结合实验变异函数曲线特征，建立各组分变程a_i与其部分基台值的联系；

e.建立变异函数结构分析与地质意义的联系，以及对资源储量估算的参数约束，进一步进行资源储量评估。

2.根据权利要求1所述的评估方法，其特征在于，所述步骤b中，针对不同的品位值r，统计大于r的数据个数N(r)，据此计算logr和logN(r)，在双对数坐标轴上，做logr～logN(r)散点图；用最小二乘法进行拟合，拟合优度R²≥0.95，拟合得到一条或多条线段，每一条线段代表一个组分；利用该线段两端点的纵坐标的反对数求差，得到该组分中的样品点个数N_i(r)，据此计算各组分的概率，

P_i＝N_i/N

其中，N为分形范围内总样品数，N_i为某组分的样品个数。

3.根据权利要求1所述的评估方法，其特征在于，所述步骤c中，改变滞后距h大小，计算实验变异函数γ*(h)值，做h-γ*(h)散点图；利用实验变异函数的头2-3个点拟合块金方差C₀；根据公式C_i＝P_i(C-C₀)，计算出各组分的部分基台值；在满足采样空间范围大于3倍变程的条件下，用样品方差估计基台值，可得C_i＝P_i(S²-C₀)。

4.根据权利要求1所述的评估方法，其特征在于，所述步骤d中，先初步确定各组分变程，再通过变程微小的调整，使实验变异函数的拟合达到最优。

5.根据权利要求1所述的评估方法，其特征在于，所述步骤e中，根据所确定的理论变异函数的C₀、C_i和a_i结构参数，确定变异函数套合结构的数学模型；根据该模型对资源储量估算的条件和参数进行约束。

6.根据权利要求5所述的评估方法，其特征在于，所述约束为确定矿化背景值和矿化结构，确定资源储量估算邻域。

7.根据权利要求6所述的评估方法，其特征在于，块金效应越小，屏蔽作用越强，估计邻域范围可以缩小；相反若块金效应大，屏蔽作用弱，估计邻域范围可以扩大。