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

地质统计学原理及其在矿床建模与储量估算中的应用中国地质调查局发展研究中心矿床品位建模及其应用需求矿体表面模型（矿化边界）矿床品位模型勘探线剖面品位分析品位－吨位曲线分析矿床品位建模及储量估算流程组合样品分析样品确定矿床块体模型参数选择插值类型设置插值参数等确定搜索邻域精度验证满意估值矿床品位模型否是回顾：地理学第一定律及应用•地理学第一定律：距离越近，两点的地理现象相似性越大逐点移面内插：以待插点为中心，确定一个邻域范围，用该邻域内的采样点计算内插点的高程值。反距离加权平均法内容介绍•地质统计学简介•区域化变量•变差函数建模•克里格插值算法•矿体储量估算应用历史背景与产生•为解决矿床从普查勘探、矿山设计到矿山开发整个过程中各种储量计算和误差估计问题发展起来的。•地质统计学是数学地质的重要分支，它首先由D·G·克立格（Krige）工程师在南非的金属矿产储量计算中使用，后由法国马特隆（G·Mathreon）教授领导的小组对此作了深入的研究并系统地总结出地质统计学的理论和方法。地质统计学定义•地质统计学(Geostatistics)是以区域化变量理论作为理论基础，以变差函数作为主要工具，对既具有随机性又具有结构性的变量（如品位值）进行研究的科学。其核心即“克里格法”，它是一种无偏的最小误差的储量计算方法。•区域化变量•变差函数•克里格估值与传统储量估算方法相比•从传统方法把部分钻孔品位当作一个块段的品位，从而使高品位估计偏高，低品位估计偏低，而且没有考虑矿石品位的空间变异性，在计算块段平均品位时，每个样品的贡献仅仅是若干个几何因素。•地质统计学方法避免了传统方法的两个缺陷。其加权因子是以矿床的各个方向变差函数的参数为基础计算出来的,这种加权方法充分考虑了矿体形态的空间变化及其品位空间变化特征,并且采用了无偏的、误差最小的数理统计方法计算样品的加权因子和块段的品位。地质统计学的发展•完善的理论基础•基本概念—区域化变量基本工具—变差函数•基本假设—本征假设基本方法—克里格法•方法与技巧不断涌出•析取克里格、多元高斯克里格和各种条件模拟技术的应用和发展•地质统计学的软件包及应用软件不断推出•美国斯坦福大学的GSLIB软件包•挪威ODEN公司的STORM随机建模软件•加拿大的Geostat地质统计学软件•澳大利亚的SurpacVision\Micromine矿山工程软件内容介绍•地质统计学简介•区域化变量•变差函数建模•克里格品位估值•矿体储量估算应用区域化变量•G.马特隆定义区域化变量是：一种在空间上具有数值的实函数，它在空间的每一个点取一个确定的数值，即当由一个点移到下一个点时，函数值是变化的.•特征：随机性和结构性•随机性•结构性区域化变量从地质及矿业角度来看，区域化变量具有如下性质：（1）空间局限性：即它被限制在一个特定的空间（如一个矿体内）；该空间称为区域化的几何域；区域化变量是按几何支撑定义的。（2）连续性：不同的区域化变量具有不同的连续性，这种连续性是通过相邻样品之间的变差函数来描述的。（3）异向性：当区域化变量在各个方向上具有相同的性质时称各向同性，否则称各向异性。（4）相关性：一定范围内、一定程度上的空间相关性，当超出这一范围后相关性减弱以至消失。（5）对于任一区域化变量而言，特殊的变异性是叠加在一般规律之上。内容介绍•地质统计学简介•区域化变量•变差函数建模•克里格插值算法•矿体储量估算应用变差函数建模•为表征一个矿床金属品位等特征量的变化，经典统计学通常采用均值、方差等一类参数，这些统计量只能概括该矿床中金属品位等特征量的全貌，却无法反映局部范围和特定方向上地质特征的变化。地质统计学引入变差函数这一工具，它能够反映区域化变量的空间变化特征——相关性和随机性，特别是透过随机性反映区域化变量的结构性，故变差函数又称结构函数。变差函数定义•我们可以把一个矿床看成是空间中的一个域，如图中为沿方向被矢量分割的两个点，其观测值分别为及，该两者的差值就是一个有明确物理意义的结构信息，因而可以看成是一个变量。•区域化变量在空间相距的任意两点和处的值与差的方差之半定义为区域化变量的变差函数，记为()Zx()Zxhx,xxhh[()()]ZxZxhVxxhhx()Zxhxxh()Zx()Zxh()Zx(,)xh1(,)[()()]2xhVarZxZxh变差函数定义•定义：在任一方向，相距的两个区域化变量和的增量的方差的一半。•公式：•变差函数值与区域化变量位置无关•二阶平稳假设和本征假设a||h()Zx()Zxhx21()[()()]2hEZxZxh二阶平稳假设当区域化变量满足下列两个条件时，称该区域化变量满足二阶平稳：（１）在整个研究区内，区域化变量的期望存在且等于常数：（常数）（２）在整个研究区内，区域化变量的空间协方差函数存在且平稳：当时，上式变成：即它有有限先验方差。()Zx0h[()]EZxmx2[(),()][()()]()CovZxZxhEZxZxhmCh,xhx[()](0)VarZxC本征假设当区域化变量的增量满足下列两个条件时，称该区域化变量满足本征假设：（１）在整个研究区内，区域化变量的增量的期望为０：（２）对于所有区域化变量的增量的方差函数存在且平稳：即要求的变差函数存在且平稳()Zx()()ZxZxh()Zx()()ZxZxh()()ZxZxh()Zx()h[()()]0EZxZxh,xh2[()()][()()]VarZxZxhEZxZxh2(,)2()xhh,xh实验变差函数计算其中:=两个样本点间的距离=样本点属性值(位置)=样本点属性值(位置)=样本点数•变差函数计算公式：h()Nh()iZx()iZxhixixh()211()[()()]2()NhiiihZxZxhNh变差函数计算实例•某地区规则采样数据，数据为属性值，样本间距为100米。实验变差函数计算实例•图中表示的是东西方向，相距为100米的样本点对。实验变差函数计算实例•通过变差函数计算公式得到东西方向上，滞后距为100米的变差函数值。实验变差函数计算实例•变差函数图：滞后距100米的变差函数点024681012141618200100200300400500滞后距变差函数实验变差函数计算实例•相距为200米的样本点对。实验变差函数计算实例•滞后距为200米的变差函数值。变差函数计算实例•变差函数图：滞后距200米的变差函数点024681012141618200100200300400500滞后距变差函数变差函数计算实例•变差函数图：滞后距300米、400米的变差函数点024681012141618200100200300400500滞后距变差函数变差函数计算实例•计算南北方向滞后距为100米、200米和300米的变差函数。实验变差函数计算实例•南北方向400m点数过少，不参与计算。滞后距东西方向南北方向1001.465.352003.39.873004.3118.884006.7变差函数值实验变差函数计算实例•变差函数图：东西方向和南北方向024681012141618200100200300400500滞后距变差函数东西方向南北方向实验变差函数计算--距离和角度容差对于不规则采样点：•沿某一特定方向和特定滞后距上并没有足够的样本点•采用距离和角度容差解决该问题520.351MineralResources3511stsemester2005SavethisslidetoEMF,SavethisslidetoEMF,theninsertandcroptheninsertandcrop实验变差函数计算•步长：4m•步长容差：2m•方位角：60•倾角：0•方位容差：22.5•倾角容差：22.5•水平带宽：5m•垂直带宽：5m实验变差函数计算(3D)变差函数的计算过程是由系统自行完成的，而合适的参数大小将直接影响计算结果的好坏。关于参数的选取实验变差函数参数选择步长大小的选择：步长间距太小步长间距较合适实验变差函数参数选择步长个数的选择：原则:步长大小*步长个数=研究区域长度的一半步长总间距理论变差函数•实验变差函数并不能定量的反映数据空间相关性，需要对实验变差函数进行拟合得到理论变差函数。•理论变差函数三参数：块金值/基台值/变程(基台值=先验方差)SamplesnotspatiallycorrelatedSamplesSpatiallyCorrelated基台值变程块金值0.......(h)h样本空间相关样本空间不相关理论变差函数模型SamplesnotspatiallycorrelatedSamplesSpatiallyCorrelated•球状模型•线性模型•指数模型•高斯模型球状模型•球状模型公式：•接近原点处，变差函数呈线性形状，在变程处达到基台值。原点处变差函数的切线在变程的2/3处与基台值相交。•实验变差函数在大多数情况下可以拟合成球状模型。因此，球状模型是应用最广的一种变差函数模型。30300031()()022hhhhcchaaaccha指数模型•指数模型公式：•变差函数渐近地逼近基台值，在实际变程处，变差函数为0.95，模型在原点处为直线。在原点处连续性最好，是一种较稳定的模型。000()(1)0hahhccehac高斯模型•高斯模型公式：•变差函数渐近地逼近基台值，在实际变程处，变差函数为0.95，模型在原点处为抛物线。为一种连续性好但稳定性较差的模型。22000()(1)0hahhccehac变差函数拟合•用球状模型、指示模型或高斯模型对实验变差函数进行拟合。•得到块金值、基台值和变程三个参数。变差函数拟合•球状模型变程为4141m，指数模型变程为5823m，高斯模型变程为2884m•观察图形：高斯模型拟合最好，其次是球状模型。•根据实际情况确定变差函数类型，结果因人而异。变差函数拟合—过程几何各向异性•基台值相同•变程不同在不同的方向具有相同的变异程度（基台值相同）但具有不同的连续程度（变程不同）为几何各向异性。带状各向异性•基台值不同•变程可同可不同在一些不同的方向上具有不同的变异程度（基台值不同）连续程度（变程）可以相同也可不同为带状各向异性。变差函数结构套合不同方向结构套合•几何各向异性•基本思路为通过线性变换将各向异性的坐标向量转化为各向同性的新坐标向量设这个线性变换为，其中对于各向同性模型，，其中对于几何各向异性变差函数，变化为矩阵形式(,,)Tuvwhhhh''''(,,)Tuvwhhhh'hAh111213212223313233aaaAaaaaaa11()(,,)()uvwhhhhh2221uvwhhhh2221()()()uvwuvwhhhhaaa100010001vuwaAaa变差函数结构套合不同方向结构套合•带状各向异性•对于带状各向异性，采用分块处理的方法。具体的变差函数模型公式为，其中对于做和几何各向异性相同的处理，对于做如下处理，对于做如下处理：总的来说，对于带状各向异性的处理方法是将其看作是几何各向异性进行坐标变换后，再分别对次轴和垂直轴方向上多出的基台值进行叠加处理。111212313()()()()hwhwhwh12()h2vvhha13()h3wwhha11()h各向异性椭球各向异性椭球：•主轴变程•次轴变程•垂直轴变程•方位角•倾角•旋转角度几何各向异性结构套合变差函数表面图Allpointsthatfallintheblockarepairedwiththepointat(x,y)tocreatethevariogrammaps.Thesizeoftheblockisthelagsize.变差函数建模•变差函数是区域化变量空间变异性的一种度量•反映了空间变异程度随