1、普通克里格法,要 点 1. 克里格法的定义 2. 克里格法的种类 3. 克里格法的使用信息和应用条件 4. 普通克里格方程组 5. 普通克里格方差 6. 算例与应用实例,一、概述,1. 克里格法的定义 矿业定义:根据一个块段(盘区)内外的若干信息样品的某特征值(品位)数据,对该块段(盘区)品位(特征值)作出一种线性、无偏、最小估计方差的估计方法。 数学定义:一种求最优、线性、无偏内插估计量的方法。 具体说是:考虑了信息样品的形状、大小及其待估块段相互之间的空间分布位置等几何特征,以及变量的空间结构信息后,为了达到线性、无偏、最小估计方差的估计,而对每个样品值分别赋予一定的权系数,最后用加权平均
2、来对待估块段的未知量进行估计的方法。“特定的滑动平均”,一、概述,2. 克里格法的种类 (1)普通克里格满足二阶平稳(或本征)假设的区域化变量 (2)泛克里格非平稳或具有漂移存在的区域化变量 (3)析取克里格计算可采储量时,需要采用非线性估计量时 (4)对数克里格区域化变量符合对数正态分布时 (5)随机克里格数据较少、分布不规则、对估计精度要求不高时 (6)因子克里格 (7)指示克里格处理非参数(类型参数)的区域化变量时,二、克里格方程组及其方差,1. 问题的提出 设Z(x)为点承载的区域化变量,且是二阶平稳(或本征)的。今要对以x0为中心的盘区V(x0)的平均值 (简记为ZV)进行估计。,2
3、. 线性估计量的构造 Zi (i=1,2, ,n)是一组离散的信息样品数据,它们定义在点承载xi (i=1,2, ,n)上的;或是确定在以xi 点为中心的承载vi (i=1,2, ,n)上的平均值Zvi (xi) (简记Zi )。且这n个承载vi (i=1,2, ,n)既不同于V,又各不相同。所采用的线性估计量为: 它是n个数值的线性组合。 克里格估值的原则:就是在保证这个估值ZV*是无偏的,且估计方差最小的前提下,求出n个权系数i 。在这样的条件下求得的i 所构成的估计量ZV*称为ZV的克里格估计量,记为ZK* 。这时的估计方差称为克里格方差,记为K*。 当Z (x)的期望已知时:为简单克里
4、格;未知时:为普通克里格,3. 简单克里格方程组和简单克里格方差(E(Z(x)=m 已知) 由于Z (x)的期望为已知,令:Y(x)=Z(x)-m 则:EY(x)=EZ(x)-m=EZ(x)-m=0 其协方差函数为: EY(x) Y( y) =C(x, y ) 对ZV的估计转化为对YV* 的估计,且有: 所用的估计量为: 只要求得YV的估计值YV* ,就能得到ZV的估计值ZV * 。,显然: YV*是YV的无偏估计量,且不需要任何条件。因为: 为了求出i,使得 最小,首先需求出 的表达式: 所以: (1) 其中 表示协方差函数在待估域V上的平均值。,公式(1)与公式(4)中,所用的估计方差符号
5、不一样,(1)式表示无偏估计量的估计方差,不能保证估计方差最小,故用记号 。公式(4)是在确保估计方差最小的前提下推导出来的,它是克里格方差,故记号为 。其中关键的区别在于i (i=1,2, ,n)在两个式中的意义不一样。 从克里格方程组解出i 后,即得到YV的简单克里格估计量: 所以:,4. 普通克里格方程组和普通克里格方差(E(Z(x)=m 未知) 要使估计量 是无偏的,就必须增加限制条件: (1)无偏性条件 若要使ZV*为ZV的无偏估计量,即要求 EZV*- ZV =0 因为: 又因为: 所以得无偏性条件为:,(2)普通克里格方程组 在区域化变量Z(x) 满足二阶平稳的条件下类似于简单克
6、里格方法的估计方差的推导,同样可以得到估计方差: 在无偏性条件 下,要使得估计方差最小,从而求得诸权系数 i , (i=1,2,n),这是一个求条件极值的问题,要用拉格朗日乘数法。 令: ,为n个权系数 i和 的(n+1)元函数。-2 是拉格朗日乘数。求出F对 i , (i=1,2,n)以及F对的偏导数,并令其为零,得到普通克里格方程组。,普通克里格方程组: 整理得: 这n+1个方程的方程组,称为普通克里格方程组。,普通克里格方差: 将上式克里格方程组中的第一式(前n个方程)两边乘以 i ,再对i 从1到n求和得: 将此式代入到普通克里格估计方差公式中得:,(3)用变差函数表示的普通克里格方程
7、组与普通克里格方差 若Z(x)只满足本征假设,而不满足二阶平稳假设时,则利用协方差函数与变差函数的关系C(h)=C(0) - (h) 可得用变差函数(h)表示的普通克里格方程组与普通克里格方差:,(4)信息样品为非点承载时的普通克里格方程组与普通克里格方差 若样品的承载不能看作是点承载,而是以x i为中心,其体积为v i的承载时,样点之间的协方差C(xi ,x j ),就变为样品域之间的平均协方差 ,相应的普通克里格方程组与普通克里格方差分别写成: 用变差函数(h)表示的普通克里格方程组与普通克里格方差:,(5)普通克里格方程组及其方差的矩阵的表示法 为简单起见,我们仅给出样品点为非点承载下的
8、普通克里格方程组及其方差的矩阵表示形式: 其中: K称为普通克里格矩阵,它是一个对称矩阵,因为有: 估计方差表示为:,用变差函数表示时,普通克里格方程组的矩阵表示形式为:,进一步的说明,进一步的说明,2. 克里格估值是一种无偏的内插估值。即若待估块段(承载)V与有效数据的任意承载vi重合,则由克里格方程组给出ZK*=Z(vi)及K2=0。这在制图学中称为“克里格估值曲面通过实测点”。传统的估计方法并没有这种性质。这也说明了克里格估值精度高于其它估值方法。,进一步的说明,3. 对于克里格方程组所用到的协方差函数C(h)和变差函数(h)的模型,不论它们所表征的基本结构如何均可,它们可以是各向同性的
9、,也可以是各向异性;既可以是单一结构,也可以是套合结构。,进一步的说明,4. 普通克里格方程组和方差只取决于结构模型C(h)或(h) ,以及各承载的相对几何特征(或说相对空间位置),而不依赖于数据Zi 的具体数值。因此,只要知道结构函数C(h)或(h)以及样品的空间位置(数据构形),在开钻前就可得普通克里格方程组及其方差。这样,就可以根据钻孔的空间位置不同,得出不同的克里格方差,从而选择较小的克里格方差所对应的钻孔位置构形,在已知结构函数前提下确定最优的布孔方案。,进一步的说明,5. 普通克里格矩阵K ,只取决样品承载vi (i=1,2,n)的几何特征(空间位置),而完全不依赖于待估块段的承载
10、V。因此,只要所用的信息样品相同,即使对不同的待估块进行估值,克里格方程组的系数矩阵K 也相同。从而只需求一次逆矩阵K -1。若估计构形(待估承载与全体样品承载的构形)也相同,则矩阵M也不变。即只需解一次克里格方程组,就可得到线性估计量中的权系数i (i=1,2,n),大大地节省计算时间。(规则勘探网格就满足这一要求),进一步的说明,6. 普通克里格方程组及其方差考虑了以下四个方面的因素: (1)待估承载V 的几何特征( (V, V) ); (2)数据构形的几何特征( (vi, vj) ); (3)信息样品承载vi 与待估承载V之间的距离( (vi, V) ) ; (4)反应区域化变量Z(x)
11、空间结构特征的变差函数模型( (h) ) 。,进一步的说明,7. 纯块金效应对普通克里格方程组及其方差的影响 如果原来变差函数1(h) ,后来增加了一个块金常数C0成为:则当所有信息样品承载vi (i=1,2,n)大小相等,和待估块段V 彼此都不相交,且V比vi 大很多(这些条件在实际中常能被满足)时,块金效应(即增加了一个块金常数)对普通克里格方程组的影响只是在原普通克里格矩阵的主对角线上前n个元素中减去块金常数。这时方程组变为:,普通克里格的计算,普通克里格的计算分两种:1、点克里格法: 样品承载很小(可看成点承载),若用它们对某一点x0处进行克里格估值,则叫点克里格法2、块段克里格法:样品承载很小(可看成点承载),若用它们对以某一点x0为中心的某一块段V进行克里格估值,则叫块段克里格法,例1-点克里格法,例2-块段克里格法,关于权系数的一点说明,克里格权系数的特点,克里格权系数的特点,克里格权系数的特点,