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

1、勘探线剖面品位分析勘探线剖面品位分析品位吨位曲线分析品位吨位曲线分析组合样品组合样品分析样品分析样品确定矿床块确定矿床块体模型参数体模型参数选择选择插值类型插值类型设置插值参数等设置插值参数等确定搜索邻域确定搜索邻域精度验证精度验证满意满意估值估值矿床品位模型矿床品位模型否是 地理学第一定律地理学第一定律： 距离越近，两点的地理现象相似性越大距离越近，两点的地理现象相似性越大逐点移面内插逐点移面内插：以待插点为中心，确定一个邻域范围，用该邻域内的采样点计算内插点的高程值。反距离加权平均法反距离加权平均法 地质统计学简介地质统计学简介 区域化变量 变差函数建模 克里格插值算法 矿体储量估算应用

2、为解决矿床从普查勘探、矿山设计到矿山开发整为解决矿床从普查勘探、矿山设计到矿山开发整个过程中各种储量计算和误差估计问题发展起来个过程中各种储量计算和误差估计问题发展起来的。的。 地质统计学是数学地质的重要分支，它首先由地质统计学是数学地质的重要分支，它首先由 D DG G克立格克立格（KrigeKrige）工程师在南非的金属矿产）工程师在南非的金属矿产储量计算中使用，后由法国储量计算中使用，后由法国马特隆马特隆（G GMathreonMathreon）教授领导的小组对此作了深入的研究并系统地总教授领导的小组对此作了深入的研究并系统地总结出地质统计学的理论和方法。结出地质统计学的理论和方法。 地

3、质统计学地质统计学(Geostatistics)是以是以区域化变量区域化变量理论理论作为理论基础，以作为理论基础，以变差函数变差函数作为主要工具，对既作为主要工具，对既具有随机性又具有结构性的变量（如品位值）进具有随机性又具有结构性的变量（如品位值）进行研究的科学。其核心即行研究的科学。其核心即“克里格法克里格法”，它是一，它是一种无偏的最小误差的储量计算方法。种无偏的最小误差的储量计算方法。 区域化变量区域化变量 变差函数变差函数 克里格估值克里格估值 从传统方法把部分钻孔品位当作一个块段的品位，从传统方法把部分钻孔品位当作一个块段的品位，从而使高品位估计偏高，低品位估计偏低，而且从而使高品

4、位估计偏高，低品位估计偏低，而且没有考虑矿石品位的空间变异性没有考虑矿石品位的空间变异性，在计算块段平，在计算块段平均品位时，每个样品的贡献仅仅是若干个几何因均品位时，每个样品的贡献仅仅是若干个几何因素。素。 地质统计学方法避免了传统方法的两个缺陷。其地质统计学方法避免了传统方法的两个缺陷。其加权因子是以加权因子是以矿床的各个方向变差函数的参数矿床的各个方向变差函数的参数为为基础计算出来的基础计算出来的, , 这种加权方法充分考虑了矿体这种加权方法充分考虑了矿体形态的空间变化及其品位空间变化特征形态的空间变化及其品位空间变化特征, , 并且采并且采用用了无偏的了无偏的、误差最小误差最小的数理统

5、计方法计算样品的数理统计方法计算样品的加权因子和块段的品位。的加权因子和块段的品位。 完善的理论基础完善的理论基础 基本概念区域化变量 基本工具变差函数 基本假设本征假设 基本方法克里格法 方法与技巧不断涌出方法与技巧不断涌出 析取克里格、多元高斯克里格和各种条件模拟技术的应用和发展 地质统计学的软件包及应用软件不断推出地质统计学的软件包及应用软件不断推出 美国斯坦福大学的GSLIB软件包 挪威ODEN公司的STORM随机建模软件 加拿大的Geostat地质统计学软件 澳大利亚的Surpac VisionMicromine矿山工程软件 地质统计学简介 区域化变量区域化变量 变差函数建模 克里格

6、品位估值 矿体储量估算应用 G.G.马特隆定义区域化变量是：一种在空间上具有马特隆定义区域化变量是：一种在空间上具有数值的实函数，它在空间的每一个点取一个确定数值的实函数，它在空间的每一个点取一个确定的数值，即当由一个点移到下一个点时，函数值的数值，即当由一个点移到下一个点时，函数值是变化的是变化的. . 特征：特征：随机性随机性和和结构性结构性 随机性随机性 结构性结构性从地质及矿业角度来看，区域化变量具有如下性质：（1 1）空间局限性）空间局限性：即它被限制在一个特定的空间（如一个矿体内）；该空间称为区域化的几何域；区域化变量是按几何支撑定义的。（2 2）连续性）连续性：不同的区域化变量具

7、有不同的连续性，这种连续性是通过相邻样品之间的变差函数来描述的。（3 3）异向性）异向性：当区域化变量在各个方向上具有相同的性质时称各向同性，否则称各向异性。（4 4）相关性）相关性：一定范围内、一定程度上的空间相关性，当超出这一范围后相关性减弱以至消失。（5）对于任一区域化变量而言，特殊的变异性是叠加在一般特殊的变异性是叠加在一般规律之上规律之上。 地质统计学简介 区域化变量 变差函数建模变差函数建模 克里格插值算法 矿体储量估算应用 为表征一个矿床金属品位等特征量的变化，经典统计学通常采用均值、方差均值、方差等一类参数，这些统计量只能概括该矿床中金属品位等特征量的全貌，却无法反映局部范围和

8、特定方向局部范围和特定方向上地质特征的变化。地质统计学引入变差函数这一工具，它能够反映区域化变量的空间变化特征相关性和随机性，特别是透过随机性反映区域化变量的结构性，故变差函数又称结构函数结构函数。我们可以把一个矿床看成是空间中的一个域，如图中 为沿 方向被矢量 分割的两个点，其观测值分别为 及 ，该两者的差值差值 就是一个有明确物理意义的结构信息，因而可以看成是一个变量。区域化变量 在空间相距 的任意两点 和 处的值 与 差的方差之半定义为区域化变量区域化变量 的变差函数的变差函数，记为( )Z x()Z xhx, x xhh ( )()Z xZ xhVxxhhx( )Z xhxxh( )Z

9、 x()Z xh( )Z x( , )x h1( , ) ( )()2x hVar Z xZ xh 定义：在任一方向，相距的两个区域化变量和的增量的方差的一半。 公式： 变差函数值与区域化变量位置无关 二阶平稳假设二阶平稳假设和本征假设本征假设a|h( )Z x()Z xhx21( ) ( )()2hE Z xZ xh当区域化变量满足下列两个条件时，称该区域化变量满足二阶平稳：（）在整个研究区内，区域化变量 的期望存在且等于常数： （常数）（）在整个研究区内，区域化变量的空间协方差函数存在且平稳：当时，上式变成：即它有有限先验方差。( )Z x0h ( )E Z xmx2 ( ),() ( )

10、 ()( )Cov Z x Z xhE Z x Z xhmC h, xh x ( )(0)Var Z xC当区域化变量的增量满足下列两个条件时，称该区域化变量满足本征假设：（）在整个研究区内，区域化变量 的增量的期望为： （）对于所有区域化变量的增量的方差函数存在且平稳：即要求的变差函数存在且平稳( )Z x( )()Z xZ xh( )Z x( )()Z xZ xh( )()Z xZ xh( )Z x( )h ( )()0E Z xZ xh, xh 2 ( )() ( )()Var Z xZ xhE Z xZ xh2 ( , )2 ( )x hh, xh 其中: = 两个样本点间的距离 =

11、样本点属性值(位置 ) = 样本点属性值(位置 ) = 样本点数 变差函数计算公式：h( )N h( )iZ x()iZ xhixixh( )211( ) ( )()2( )N hiiihZ xZ xhN h 某地区规则采样数据，数据为属性值，样本间距为100米。 图中表示的是东西方向，相距为100米的样本点对。 通过变差函数计算公式得到东西方向上，滞后距为100米的变差函数值。 变差函数图：滞后距100米的变差函数点024681012141618200100200300400500滞后距变差函数 相距为200米的样本点对。 滞后距为200米的变差函数值。 变差函数图：滞后距200米的变差函数

12、点024681012141618200100200300400500滞后距变差函数 变差函数图：滞后距300米、400米的变差函数点024681012141618200100200300400500滞后距变差函数 计算南北方向滞后距为100米、200米和300米的变差函数。 南北方向400m点数过少，不参与计算。滞后距东西方向南北方向1001.465.352003.39.873004.3118.884006.7变差函数值 变差函数图：东西方向和南北方向024681012141618200100200300400500滞后距变差函数东西方向南北方向对于不规则采样点不规则采样点： 沿某一特定方向和

13、特定滞后距上并没有足够的样本点 采用距离和角度容差解决该问题步长：步长：4m4m步长容差：步长容差：2m2m方位角：方位角：6060倾角：倾角：0 0方位容差：方位容差：22.522.5倾角容差：倾角容差：22.522.5水平带宽：水平带宽：5m5m垂直带宽：垂直带宽：5m5m 变差函数的计算过程是由系统自行完变差函数的计算过程是由系统自行完成的，而合适的参数大小将直接影响计算成的，而合适的参数大小将直接影响计算结果的结果的好坏。好坏。 关于关于参数参数的选取的选取步长大小的选择步长大小的选择：步长间距太小步长间距太小步长间距较合适步长间距较合适步长个数的选择步长个数的选择：原则: 步长大小*

14、步长个数=研究区域长度的一半步长总间距步长总间距 实验变差函数并不能定量的反映数据空间相关性，需要对实验变差函数进行拟合得到理论变差函数。 理论变差函数三参数：块金值块金值/基台值基台值/变程变程(基台值=先验方差)Samples not spatially correlatedSamplesSpatially Correlated基台值基台值变程变程块金值块金值0.(h)h样本空间相关样本空间不相关Samples not spatially correlatedSamplesSpatially Correlated 球状模型 线性模型 指数模型 高斯模型 球状模型公式： 接近原点处，变差函数

15、呈线性形状线性形状，在变程处达到基台值。原点处变差函数的切线在变程的2/3处与基台值相交。 实验变差函数在大多数情况下可以拟合成球状模型。因此，球状模型是应用最广的一种变差函数模型球状模型是应用最广的一种变差函数模型。30300031( )()022hhhhcchaaaccha 指数模型公式： 变差函数渐近地逼近基台值，在实际变程 处，变差函数为0.95 ，模型在原点处为直线。在原点处连续性最好，是一种较稳定的模型。000( )(1)0hahhccehac 高斯模型公式： 变差函数渐近地逼近基台值，在实际变程 处，变差函数为0.95 ，模型在原点处为抛物线。为一种连续性好但稳定性较差的模型。2

16、2000( )(1)0hahhccehac 用球状模型、指示模型或高斯模型对实验变差函数进行拟合。 得到块金值、基台值和变程三个参数。 球状模型变程为4141m，指数模型变程为5823m，高斯模型变程为2884m 观察图形：高斯模型拟合最好，其次是球状模型。 根据实际情况确定变差函数类型，结果因人而异。 基台值相同 变程不同 在不同的方向具有相同的变异程在不同的方向具有相同的变异程度度（基台值相同）（基台值相同）但具有不同的但具有不同的连续程度连续程度（变程不同）（变程不同）为几何各为几何各向异性。向异性。 基台值不同 变程可同可不同 在一些不同的方向上具有不同的变异在一些不同的方向上具有不同

17、的变异程度程度（基台值不同）（基台值不同）连续程度连续程度（变程）（变程）可以相同也可不同为带状各向异性可以相同也可不同为带状各向异性。不同方向结构套合不同方向结构套合 几何各向异性基本思路为通过线性变换将各向异性的坐标向量 转化为各向同性的新坐标向量 设这个线性变换为 ，其中对于各向同性模型， ，其中 对于几何各向异性变差函数 ，变化为矩阵形式(,)Tuvwhh h h(,)TuvwhhhhhAh111213212223313233aaaAaaaaaa11( )(,)( )uvwhh h hh2221uvwhhhh2221()()()uvwuvwhhhhaaa100010001vuwaAaa

18、不同方向结构套合不同方向结构套合 带状各向异性对于带状各向异性，采用分块处理的方法。具体的变差函数模型公式为 ，其中对于 做和几何各向异性相同的处理，对于 做如下处理 ，对于 做如下处理： 总的来说，对于带状各向异性的处理方法是将其看作是几何各向异性进行坐标变换后，再分别对次轴和垂直轴方向上多出的基台值进行叠加处理。1 112123 13( )( )()()hwhwhwh12()h2vvhha13()h3wwhha11( )h各向异性椭球：各向异性椭球：主轴变程次轴变程垂直轴变程方位角倾角旋转角度All points that fall in the block are paired with

19、 the point at (x,y) to create the variogram maps. The size of the block is the lag size. 变差函数是区域化变量空间变异性的一种度量 反映了空间变异程度随距离而变化的特征 可定量的描述区域化变量的空间相关性 地理学第一定律：地理学第一定律： 距离越近的点相似性越大距离越近的点相似性越大 地质统计学简介 区域化变量 变差函数建模 克里格插值算法克里格插值算法 矿体储量估算应用 如果要估算 的值，一般情况下 的值应该是 的平均值，并且 随着距离 的增大而减小。 duu 克里格插值算法建立在变差函数及结构分析理论克

20、里格插值算法建立在变差函数及结构分析理论之上之上 适用条件是变差函数及相关分析的结果表明样品适用条件是变差函数及相关分析的结果表明样品间存在空间相关性间存在空间相关性 其实质是利用其实质是利用区域化变量的原始数据区域化变量的原始数据和和变差函数变差函数的结构特点的结构特点，对未采样点的区域化变量的取值进，对未采样点的区域化变量的取值进行线性、无偏、最优估计。行线性、无偏、最优估计。组合样品分析样品确定块体模型参数选择克里格类型计算实验变差函数并拟合确定搜索邻域交叉验证满意估值品位模型否是 将整个研究区域划分为多个规则小块，分别对小块属性进行估值。 起始点坐标起始点坐标 块大小块大小 块个数块个

21、数 搜索椭圆 直接定义点数 交叉验证的原理为将原始的样品点去除，然后采用原始样品点周围的点来进行克里格估值得到原始样品的估计值，最后做出原始样品和估值样品的散点图，并对估值误差进行统计。 从矿业上的术语具体来说，它是根据一个从矿业上的术语具体来说，它是根据一个待估块待估块段邻域内的若干信息样品的品位数据段邻域内的若干信息样品的品位数据，在考虑了，在考虑了这些这些样品的形状样品的形状、大小大小及及相互位置关系相互位置关系，它们与，它们与待估块段相互之间的空间位置等几何特征，以及待估块段相互之间的空间位置等几何特征，以及品位的变差函数模型所提供的结构信息之后，为品位的变差函数模型所提供的结构信息之

22、后，为了对该块段品位作出一种了对该块段品位作出一种线性、无偏、最小估计线性、无偏、最小估计方差的估计方差的估计而对每个样品值分别赋予一定的权系而对每个样品值分别赋予一定的权系数，最后进行数，最后进行加权平均加权平均来估计该块段品位的方法。来估计该块段品位的方法。 克里格插值算法克里格插值算法: : B.L.U.EB.L.U.E - - best, linear, unbiased estimator best = 最小估计误差 linear = 线性估值方式（同距离反比估值） unbiased = 无偏估计，估计误差之和为0 estimator = 估值方法其中 = 待估点位置和其中一个邻接点

23、 位置 = 估算未知点 用到的邻接点个数 = 和 对应的预测平均值 = 对应的克里格权重 克里格插值公式, u u( )n u*( )Zu( ),()m u m u( )Z u( )u()Z u()Z u( )*1( )( ) ()()n uZum uZ um u关键在于确定邻接权重 最小方差限制条件 无偏估计限制条件 克里格插值公式( )*1( )( ) ()()n uZum uZ um u2*( )( )( )EuVar ZuZ u*( )( )0E ZuZ u 最常用克里格的三种类型最常用克里格的三种类型 简单克里格简单克里格 普通克里格普通克里格 泛克里格泛克里格 其区别在于 的确定方

24、式不同 非线性克里格非线性克里格 指示克里格（Indicator Kriging） 多元高斯克里格（Multi-Gauss kriging) 协克里格（ Cokriging） 块克里格块克里格(Block (Block krigingkriging) )( )mu 六个样本点数据，给出样本点间距、样本点属性值和变差函数变差函数模型：球状模型，块金值0，基台值0.78，变程4141m 3( )(0)( )0.78 (1 1.5 ( /4141)0.5 ( /4141) )C hChhhPnt1Pnt2Pnt3Pnt4Pnt5Pnt6Pnt1018973130244114001265Pnt2189

25、701281145619702280Pnt3313012810152328003206Pnt4244114561523015231970Pnt514401970280015230447Pnt612652280320619704470 由简单克里格插值公式得 0.38,0.56,0.32,0.49,0.46,0.37Tk 0.780.280.060.170.400.430.280.780.430.390.270.200.060.430.780.370.110.060.170.390.370.780.370.270.400.270.110.370.780.650.430.200.060.270.6

26、50.78K( )SKKuk111( )( )1( )( )()()()()n un un un uC uuC uuKC uuC uu1( )()()n uC uukC uu 由简单克里格插值公式得 12314560.14750.45640.02050.27090.25340.0266K k 已知该采样数据平均值为14.70 六个采样点数据的属性值分别为 13.84, 12.15, 12.87, 12.68, 14.41, 14.59 由 得 0.862.551.83( )0.150.460.020.270.250.031.872.010.280.10R u *( )( )ZuR um*( )

27、1.87 14.7012.83Zu 普通克里格插值： 未知 普通克里格估值公式为 由无偏最优估计限制条件无偏最优估计限制条件，构建拉格朗日函数( )( )*11( )( ) ()1( )( )n un uZ uu Z uu m u( )m u( )*1( )( ) ()n uOKOKZuu Z u其中( )1( )1n uOKu( )L u( )21( ),1, ( );2( )( )2( )( ) 1n uOKOKOKEOKLun uuuuu( )*1( )( ) ()()n uZum uZ um u 用“拉格朗日乘数法”求函数f(x,y,z)在条件(x,y,z)=0下的极值方法（步骤）是：

28、 1.做拉格朗日函数L=f(x,y,z)+(x,y,z),称拉格朗日乘数 2.求L分别对x,y,z,求偏导,得方程组,求出驻点P(x,y,z) 如果这个实际问题的最大或最小值存在,一般说来驻点唯一,于是最值可求. 求偏导分别得到下列公式 得到求得到求 的方程组的方程组 估值方差计算公式估值方差计算公式( )11( )( )()()( )02( )n uOKRROKOKL uu CuuCuuuu( )11( )( ) 102( )n uOKOKL uuu 1, ( )n u( )21( )(0)( ) ()( )n uOKOKOKuCu C uuu( )1( )1( ) ()( )()( )1n

29、 uOKOKn uOKu C uuuC uuu1, ( )n u 在地质、物化探数据处理及矿产储量计算中影响计算精度的因素有很多，但主要有以下几个问题： (1)特异值的出现，所谓特异值是指那些比全部数值的平均值或中位数高得多的数值，它既非分析误差所致，也非采样方法等人为误差引起。而是实际存在于所研究的母体之中。这些特异值只占全部数据的极少部分，但却控制了总金属资源量的很大比例。 (2)在一个研究区域或一个矿床中存在几个不同类型的矿化作用，这也影响了品位和储量的精确估计。 为解决上述问题，指示克里格法应运而生，它是在不必去掉重要而实际存在的高值数据的条件下来处理不同的现象，而且给出在一定风险概率

30、条件下未知量的估计值及空间分布。( )Z x 指示克里格是一种非参数地质统计学方法。它是指示克里格是一种非参数地质统计学方法。它是根据一系列的临界值，例如边界品位，先对原根据一系列的临界值，例如边界品位，先对原始数据如下进行转换始数据如下进行转换 然后对转换后的数值求变差函数、进行克里格估然后对转换后的数值求变差函数、进行克里格估值。值。z( )Z x1,( )( , )0,( )Z xzi x zZ xz 简单克里格（整体平稳）简单克里格（整体平稳） 普通克里格（局部平稳）普通克里格（局部平稳） 泛克里格（整体存在趋势）泛克里格（整体存在趋势） 指示克里格（数据不平稳存在极值）指示克里格（数

31、据不平稳存在极值） 协克里格（存在辅助变量）协克里格（存在辅助变量）( ),m umu 已知( )m uu( )W u常量但未知0( )( )( )Kkkkm ua ufu其中 常量但未知( )kka uau( )W u 变差函数参数变差函数参数 块金值块金值：块金值越小，距离越近的点越重要，这样会导致权值的变化范围变大（从负值到大于1的值变化），使数据出现异常。块金值越大，估值结果越平滑。 变程：变程：变程小于任意两点之间的距离，变差函数为纯块金模型，随着变程增大，已知点的位置以及丛聚性变的重要。如果变程很大存在基台值，则变差函数相当于纯块金模型，如果变程很大但不存在基台值，则变差函数为线性

32、模型。 比例比例：比例大小只与克里格估值方差有关，与估值结果无关。 形状形状：球状模型与指数模型在原点位置接近线性关系，指数模型在原点处更加陡峭一些，与变程较小的球状模型相似。高斯模型在原点处为抛物线形式，这种模型要求原始变量要有很高的连续性，否则会出现大量负的权值情况，使估值结果极不稳定。 各向异性各向异性：以各个方向变程不同反映出来，可以理解为对原始数据坐标位置的变化，通过几何校正转化为各向同性情况。 屏蔽效应屏蔽效应第一个点在任何情况下获得的权值更大一些。第二个点的权值小于第第一个点在任何情况下获得的权值更大一些。第二个点的权值小于第一个点，在克里格估值过程中，第二个点的权值可能变为负值

33、，这样一个点，在克里格估值过程中，第二个点的权值可能变为负值，这样第二个点的信息变得重要。例如：当一个负值权重可能隐含某种趋势第二个点的信息变得重要。例如：当一个负值权重可能隐含某种趋势（如果第一个点比第二个点要小，那么估值点可能小与这两个点。（如果第一个点比第二个点要小，那么估值点可能小与这两个点。随着块金值的增大屏蔽效应减弱。所用样品点对估值结果作用基本相随着块金值的增大屏蔽效应减弱。所用样品点对估值结果作用基本相同。同。0.22 0.2 0.18 搜索策略搜索策略 最小值：最小值：作用是限制克里格估值过程中数据个数，如果在搜索椭球范围内，数据个数过少则不进行估值，否则估值结果可能存在很大误差。 最大值最大值：估值过程中用到的数据越多，估值结果越平滑。由于数据的局部平稳性，数据个数不能太多，否则会使不存在相关性的数据仍然参与到估值过程中。 搜索半径搜索半径：搜索半径越大估值结果越平稳，搜索半径与最大值相互影响。搜索半径设置与变差函数变程相关。 克里格估值类型克里格估值类型 简单克里格简单克里格：简单克里格需要一个固定的平均值。在估值过程中平均值被赋予一个相当大的权值，从而导致估值结果过渡平滑。 普通克里格：普通克里格：与简单克里格类似，但平滑性减弱。普通克里格估值权值之和为1，可以合理的减少丛聚效应。 泛克里格：泛克里格：与普通克里格类似，但估值结果