地质统计学原理及其在矿床建模与储量估算中的应用

会计学1地质统计学原理及其在矿床建模与储量估算中的应用矿床品位建模及其应用需求矿体表面模型（矿化边界）矿床品位模型勘探线剖面品位分析品位－吨位曲线分析第1页/共98页矿床品位建模及储量估算流程组合样品分析样品确定矿床块体模型参数选择插值类型设置插值参数等确定搜索邻域精度验证满意估值矿床品位模型否是第2页/共98页回顾：地理学第一定律及应用地理学第一定律：距离越近，两点的地理现象相似性越大逐点移面内插：以待插点为中心，确定一个邻域范围，用该邻域内的采样点计算内插点的高程值。反距离加权平均法第3页/共98页内容介绍地质统计学简介区域化变量变差函数建模克里格插值算法矿体储量估算应用第4页/共98页历史背景与产生为解决矿床从普查勘探、矿山设计到矿山开发整个过程中各种储量计算和误差估计问题发展起来的。地质统计学是数学地质的重要分支，它首先由

D·G·克立格（Krige）工程师在南非的金属矿产储量计算中使用，后由法国马特隆（G·Mathreon）教授领导的小组对此作了深入的研究并系统地总结出地质统计学的理论和方法。第5页/共98页地质统计学定义地质统计学(Geostatistics)是以区域化变量理论作为理论基础，以变差函数作为主要工具，对既具有随机性又具有结构性的变量（如品位值）进行研究的科学。其核心即“克里格法”，它是一种无偏的最小误差的储量计算方法。区域化变量变差函数克里格估值第6页/共98页与传统储量估算方法相比从传统方法把部分钻孔品位当作一个块段的品位，从而使高品位估计偏高，低品位估计偏低，而且没有考虑矿石品位的空间变异性，在计算块段平均品位时，每个样品的贡献仅仅是若干个几何因素。地质统计学方法避免了传统方法的两个缺陷。其加权因子是以矿床的各个方向变差函数的参数为基础计算出来的,这种加权方法充分考虑了矿体形态的空间变化及其品位空间变化特征,并且采用了无偏的、误差最小的数理统计方法计算样品的加权因子和块段的品位。第7页/共98页地质统计学的发展完善的理论基础基本概念—区域化变量基本工具—变差函数基本假设—本征假设基本方法—克里格法方法与技巧不断涌出析取克里格、多元高斯克里格和各种条件模拟技术的应用和发展地质统计学的软件包及应用软件不断推出

美国斯坦福大学的GSLIB软件包

挪威ODEN公司的STORM随机建模软件

加拿大的Geostat地质统计学软件

澳大利亚的SurpacVision\Micromine矿山工程软件第8页/共98页内容介绍地质统计学简介区域化变量变差函数建模克里格品位估值矿体储量估算应用第9页/共98页区域化变量G.马特隆定义区域化变量是：一种在空间上具有数值的实函数，它在空间的每一个点取一个确定的数值，即当由一个点移到下一个点时，函数值是变化的.特征：随机性和结构性

随机性

结构性第10页/共98页区域化变量从地质及矿业角度来看，区域化变量具有如下性质：

（1）空间局限性：即它被限制在一个特定的空间（如一个矿体内）；该空间称为区域化的几何域；区域化变量是按几何支撑定义的。

（2）连续性：不同的区域化变量具有不同的连续性，这种连续性是通过相邻样品之间的变差函数来描述的。

（3）异向性：当区域化变量在各个方向上具有相同的性质时称各向同性，否则称各向异性。

（4）相关性：一定范围内、一定程度上的空间相关性，当超出这一范围后相关性减弱以至消失。 （5）对于任一区域化变量而言，特殊的变异性是叠加在一般规律之上。第11页/共98页内容介绍地质统计学简介区域化变量变差函数建模克里格插值算法矿体储量估算应用第12页/共98页变差函数建模为表征一个矿床金属品位等特征量的变化，经典统计学通常采用均值、方差等一类参数，这些统计量只能概括该矿床中金属品位等特征量的全貌，却无法反映局部范围和特定方向上地质特征的变化。地质统计学引入变差函数这一工具，它能够反映区域化变量的空间变化特征——相关性和随机性，特别是透过随机性反映区域化变量的结构性，故变差函数又称结构函数。第13页/共98页变差函数定义我们可以把一个矿床看成是空间中的一个域，如图中

为沿方向被矢量分割的两个点，其观测值分别为及，该两者的差值就是一个有明确物理意义的结构信息，因而可以看成是一个变量。区域化变量在空间相距的任意两点和处的值与差的方差之半定义为区域化变量的变差函数，记为

第14页/共98页变差函数定义定义：在任一方向，相距的两个区域化变量和的增量的方差的一半。公式：变差函数值与区域化变量位置无关二阶平稳假设和本征假设第15页/共98页二阶平稳假设当区域化变量满足下列两个条件时，称该区域化变量满足二阶平稳：（１）在整个研究区内，区域化变量的期望存在且等于常数：

（常数）（２）在整个研究区内，区域化变量的空间协方差函数存在且平稳：

当时，上式变成：

即它有有限先验方差。第16页/共98页本征假设当区域化变量的增量满足下列两个条件时，称该区域化变量满足本征假设：（１）在整个研究区内，区域化变量的增量的期望为０：

（２）对于所有区域化变量的增量的方差函数存在且平稳：

即要求的变差函数存在且平稳第17页/共98页实验变差函数计算其中:

=两个样本点间的距离

=样本点属性值(位置)=样本点属性值(位置)=样本点数变差函数计算公式：第18页/共98页变差函数计算实例某地区规则采样数据，数据为属性值，样本间距为100米。第19页/共98页实验变差函数计算实例图中表示的是东西方向，相距为100米的样本点对。第20页/共98页实验变差函数计算实例通过变差函数计算公式得到东西方向上，滞后距为100米的变差函数值。第21页/共98页实验变差函数计算实例变差函数图：滞后距100米的变差函数点024681012141618200100200300400500滞后距变差函数第22页/共98页实验变差函数计算实例相距为200米的样本点对。第23页/共98页实验变差函数计算实例滞后距为200米的变差函数值。第24页/共98页变差函数计算实例变差函数图：滞后距200米的变差函数点024681012141618200100200300400500滞后距变差函数第25页/共98页变差函数计算实例变差函数图：滞后距300米、400米的变差函数点024681012141618200100200300400500滞后距变差函数第26页/共98页变差函数计算实例计算南北方向滞后距为100米、200米和300米的变差函数。第27页/共98页实验变差函数计算实例南北方向400m点数过少，不参与计算。滞后距东西方向南北方向1001.465.352003.39.873004.3118.884006.7变差函数值第28页/共98页实验变差函数计算实例变差函数图：东西方向和南北方向024681012141618200100200300400500滞后距变差函数东西方向南北方向第29页/共98页实验变差函数计算--距离和角度容差对于不规则采样点：沿某一特定方向和特定滞后距上并没有足够的样本点采用距离和角度容差解决该问题第30页/共98页实验变差函数计算步长：4m步长容差：2m方位角：60倾角：0方位容差：22.5倾角容差：22.5水平带宽：5m垂直带宽：5m第31页/共98页实验变差函数计算(3D)第32页/共98页变差函数的计算过程是由系统自行完成的，而合适的参数大小将直接影响计算结果的好坏。关于参数的选取第33页/共98页实验变差函数参数选择步长大小的选择：步长间距太小步长间距较合适第34页/共98页实验变差函数参数选择步长个数的选择：原则:步长大小*步长个数=研究区域长度的一半步长总间距第35页/共98页理论变差函数实验变差函数并不能定量的反映数据空间相关性，需要对实验变差函数进行拟合得到理论变差函数。理论变差函数三参数：块金值/基台值/变程(基台值=先验方差)SamplesnotspatiallycorrelatedSamplesSpatiallyCorrelated基台值变程块金值0.......(h)gh样本空间相关样本空间不相关第36页/共98页理论变差函数模型SamplesnotspatiallycorrelatedSamplesSpatiallyCorrelated球状模型线性模型指数模型高斯模型第37页/共98页球状模型球状模型公式：接近原点处，变差函数呈线性形状，在变程处达到基台值。原点处变差函数的切线在变程的2/3处与基台值相交。实验变差函数在大多数情况下可以拟合成球状模型。因此，球状模型是应用最广的一种变差函数模型。第38页/共98页指数模型指数模型公式：变差函数渐近地逼近基台值，在实际变程处，变差函数为0.95，模型在原点处为直线。在原点处连续性最好，是一种较稳定的模型。第39页/共98页高斯模型高斯模型公式：变差函数渐近地逼近基台值，在实际变程处，变差函数为0.95，模型在原点处为抛物线。为一种连续性好但稳定性较差的模型。第40页/共98页变差函数拟合用球状模型、指示模型或高斯模型对实验变差函数进行拟合。得到块金值、基台值和变程三个参数。第41页/共98页变差函数拟合球状模型变程为4141m，指数模型变程为5823m，高斯模型变程为2884m观察图形：高斯模型拟合最好，其次是球状模型。根据实际情况确定变差函数类型，结果因人而异。第42页/共98页变差函数拟合—过程第43页/共98页几何各向异性基台值相同变程不同在不同的方向具有相同的变异程度（基台值相同）但具有不同的连续程度（变程不同）为几何各向异性。第44页/共98页带状各向异性基台值不同变程可同可不同

在一些不同的方向上具有不同的变异程度（基台值不同）连续程度（变程）可以相同也可不同为带状各向异性。第45页/共98页变差函数结构套合不同方向结构套合几何各向异性基本思路为通过线性变换将各向异性的坐标向量

转化为各向同性的新坐标向量设这个线性变换为

，其中对于各向同性模型，，其中对于几何各向异性变差函数，变化为矩阵形式第46页/共98页变差函数结构套合不同方向结构套合带状各向异性对于带状各向异性，采用分块处理的方法。具体的变差函数模型公式为，其中对于

做和几何各向异性相同的处理，对于

做如下处理

，对于做如下处理：

总的来说，对于带状各向异性的处理方法是将其看作是几何各向异性进行坐标变换后，再分别对次轴和垂直轴方向上多出的基台值进行叠加处理。第47页/共98页各向异性椭球各向异性椭球：主轴变程次轴变程垂直轴变程方位角倾角旋转角度第48页/共98页几何各向异性结构套合第49页/共98页变差函数表面图Allpointsthatfallintheblockarepairedwiththepointat(x,y)tocreatethevariogrammaps.Thesizeoftheblockisthelagsize.第50页/共98页变差函数建模变差函数是区域化变量空间变异性的一种度量反映了空间变异程度随距离而变化的特征可定量的描述区域化变量的空间相关性地理学第一定律：距离越近的点相似性越大第51页/共98页内容介绍地质统计学简介区域化变量变差函数建模克里格插值算法矿体储量估算应用第52页/共98页如果要估算的值，一般情况下的值应该是的平均值，并且随着距离的增大而减小。第53页/共98页克里格插值算法克里格插值算法建立在变差函数及结构分析理论之上适用条件是变差函数及相关分析的结果表明样品间存在空间相关性其实质是利用区域化变量的原始数据和变差函数的结构特点，对未采样点的区域化变量的取值进行线性、无偏、最优估计。第54页/共98页克里格插值过程组合样品分析样品确定块体模型参数选择克里格类型计算实验变差函数并拟合确定搜索邻域交叉验证满意估值品位模型否是第55页/共98页块体模型定义将整个研究区域划分为多个规则小块，分别对小块属性进行估值。起始点坐标块大小块个数第56页/共98页搜索邻域确定搜索椭圆直接定义点数第57页/共98页交叉验证交叉验证的原理为将原始的样品点去除，然后采用原始样品点周围的点来进行克里格估值得到原始样品的估计值，最后做出原始样品和估值样品的散点图，并对估值误差进行统计。第58页/共98页克里格插值算法从矿业上的术语具体来说，它是根据一个待估块段邻域内的若干信息样品的品位数据，在考虑了这些样品的形状、大小及相互位置关系，它们与待估块段相互之间的空间位置等几何特征，以及品位的变差函数模型所提供的结构信息之后，为了对该块段品位作出一种线性、无偏、最小估计方差的估计而对每个样品值分别赋予一定的权系数，最后进行加权平均来估计该块段品位的方法。第59页/共98页克里格插值算法克里格插值算法:

B.L.U.E

-best,linear,unbiasedestimator

best=

最小估计误差

linear=

线性估值方式（同距离反比估值）

unbiased=无偏估计，估计误差之和为0

estimator=

估值方法第60页/共98页克里格插值算法定义其中

=待估点位置和其中一个邻接点位置

=估算未知点用到的邻接点个数

=和对应的预测平均值

=对应的克里格权重

克里格插值公式第61页/共98页克里格插值算法基础关键在于确定邻接权重最小方差限制条件无偏估计限制条件克里格插值公式第62页/共98页克里格插值类型最常用克里格的三种类型简单克里格普通克里格泛克里格其区别在于的确定方式不同

非线性克里格

指示克里格（IndicatorKriging）

多元高斯克里格（Multi-Gausskriging)协克里格（Cokriging）块克里格(Blockkriging)第63页/共98页简单克里格插值应用实例六个样本点数据，给出样本点间距、样本点属性值和变差函数变差函数模型：球状模型，块金值0，基台值0.78，变程4141m

Pnt1Pnt2Pnt3Pnt4Pnt5Pnt6Pnt1018973130244114001265Pnt2189701281145619702280Pnt3313012810152328003206Pnt4244114561523015231970Pnt514401970280015230447Pnt612652280320619704470第64页/共98页简单克里格插值应用实例由简单克里格插值公式得

第65页/共98页简单克里格插值应用实例由简单克里格插值公式得

第66页/共98页简单克里格插值应用实例已知该采样数据平均值为14.70六个采样点数据的属性值分别为

13.84,

12.15,

12.87,

12.68,

14.41,

14.59

由得

第67页/共98页简单克里格插值应用实例第68页/共98页普通克里格插值普通克里格插值：

未知普通克里格估值公式为由无偏最优估计限制条件，构建拉格朗日函数其中第69页/共98页拉格朗日乘数法用“拉格朗日乘数法”求函数f(x,y,z)在条件φ(x,y,z)=0下的极值方法（步骤）是：

1.做拉格朗日函数L=f(x,y,z)+λφ(x,y,z),λ称拉格朗日乘数

2.求L分别对x,y,z,λ求偏导,得方程组,求出驻点P(x,y,z)

如果这个实际问题的最大或最小值存在,一般说来驻点唯一,于是最值可求.第70页/共98页普通克里格插值求偏导分别得到下列公式得到求的方程组估值方差计算公式第71页/共98页指示克里格估值在地质、物化探数据处理及矿产储量计算中影响计算精度的因素有很多，但主要有以下几个问题：(1)特异值的出现，所谓特异值是指那些比全部数值的平均值或中位数高得多的数值，它既非分析误差所致，也非采样方法等人为误差引起。而是实际存在于所研究的母体之中。这些特异值只占全部数据的极少部分，但却控制了总金属资源量的很大比例。(2)在一个研究区域或一个矿床中存在几个不同类型的矿化作用，这也影响了品位和储量的精确估计。为解决上述问题，指示克里格法应运而生，它是在不必去掉重要而实际存在的高值数据的条件下来处理不同的现象，而且给出在一定风险概率条件下未知量的估计值及空间分布。第72页/共98页指示克里格估值指示克里格是一种非参数地质统计学方法。它是根据一系列的临界值，例如边界品位，先对原始数据如下进行转换然后对转换后的数值求变差函数、进行克里格估值。第73页/共98页总结简单克里格（整体平稳）普通克里格（局部平稳）泛克里格（整体存在趋势）指示克里格（数据不平稳存在极值）协克里格（存在辅助变量）已知常量但未知其中常量但未知第74页/共98页几点注意内容变差函数参数

块金值：块金值越小，距离越近的点越重要，这样会导致权值的变化范围变大（从负值到大于1的值变化），使数据出现异常。块金值越大，估值结果越平滑。

变程：变程小于任意两点之间的距离，变差函数为纯块金模型，随着变程增大，已知点的位置以及丛聚性变的重要。如果变程很大存在基台值，则变差函数相当于纯块金模型，如果变程很大但不存在基台值，则变差函数为线性模型。

比例：比例大小只与克里格估值方差有关，与估值结果无关。

形状：球状模型与指数模型在原点位置接近线性关系，指数模型在原点处更加陡峭一些，与变程较小的球状模型相似。高斯模型在原点处为抛物线形式，这种模型要求原始变量要有很高的连续性，否则会出现大量负的权值情况，使估值结果极不稳定。

各向异性：以各个方向变程不同反映出来，可以理解为对原始数据坐标位置的变化，通过几何校正转化为各向同性情况。第75页/共98页几点注意内容屏蔽效应第一个点在任何情况下获得的权值更大一些。第二个点的权值小于第一个点，在克里格估值过程中，第二个点的权值可能变为负值，这样第二个点的信息变得重要。例如：当一个负值权重可能隐含某种趋势（如果第一个点比第二个点要小，那么估值点可能小与这两个点。随着块金值的增大屏蔽效应减弱。所用样品点对估值结果作用基本相同。0.220.20.18第76页/共98页几点注意内容搜索策略最小值：作用是限制克里格估值过程中数据个数，如果在搜索椭球范围内，数据个数过少则不进行估值，否则估值结果可能存在很大误差。

最大值：估值过程中用到的数据越多，估值结果越平滑。由于数据的局部平稳性，数据个数不能太多，否则会使不存在相关性的数据仍然参与到估值过程中。

搜索半径：搜索半径越大估值结果越平稳，搜索半径与最大值相互影响。搜索半径设置与变差函数变程相关。第77页/共98页几点注意内容克里格估值类型

简单克里格：简单克里格需要一个固定的平均值。在估值过程中平均值被赋予一个相当大的权值，从而导致估值结果过渡平滑。

普通克里格：与简单克里格类似，但平滑性减弱。普通克里格估值权值之和为1，可以合理的减少丛聚效应。

泛克里格：与普通克里格类似，但估值结果边缘会存在一个奇异值，这是由趋势面构造的不合理造成的。

块克里格：由于块克里格估值的平滑效应会使克里格估值结果方差减小。随着块的增大，估值结果越来越平滑。第78页/共98页内容介绍地质统计学简介区域化变量变差函数建模克里格品位估值矿体储量估算应用第79页/共98页矿床品位建模及储量估算流程组合样品分析样品确定矿床块体模型参数选择克里格类型计算变差函数并拟合确定搜索邻域精度验证满意估值矿床品位模型否是第80页/共98页克里格法资源储量估算流程－关键步骤否否是是数据分析数据变换选择泛克里格估值选择其它克里格估值剔出趋势导入组合样品数据是否服从正态分布是否存在趋势数据分析与变换第81页/共98页克里格法资源储量估算流程－关键步骤计算实验变差函数并绘制变差函数表面图是否服从各向异性单方向变差函数拟合确定最大变程方向多方向变差函数拟合是否球状模型拟合变差函数拟合单方向变差函数拟合多方向变差函数拟合指数模型拟合高斯模型拟合几何各向异性带状

各向异性变差函数计算及拟合数据分析与变换结