地质统计学(6)-普通克里格法-cjg2011

普通克里格法要点1.克里格法的定义2.克里格法的种类3.克里格法的使用信息和应用条件4.普通克里格方程组5.普通克里格方差6.算例与应用实例一、概述1.克里格法的定义矿业定义：根据一个块段（盘区）内外的若干信息样品的某特征值（品位）数据，对该块段（盘区）品位（特征值）作出一种线性、无偏、最小估计方差的估计方法。数学定义：一种求最优、线性、无偏内插估计量的方法。具体说是：考虑了信息样品的形状、大小及其待估块段相互之间的空间分布位置等几何特征，以及变量的空间结构信息后，为了达到线性、无偏、最小估计方差的估计，而对每个样品值分别赋予一定的权系数，最后用加权平均来对待估块段的未知量进行估计的方法。“特定的滑动平均”一、概述2.克里格法的种类（1）普通克里格——满足二阶平稳（或本征）假设的区域化变量（2）泛克里格——非平稳或具有漂移存在的区域化变量（3）析取克里格——计算可采储量时，需要采用非线性估计量时（4）对数克里格——区域化变量符合对数正态分布时（5）随机克里格——数据较少、分布不规则、对估计精度要求不高时（6）因子克里格——（7）指示克里格——处理非参数（类型参数）的区域化变量时3.克里格法的使用信息及应用条件信息：①一组数据；②空间构形（坐标）；③结构信息（变差函数模型）。条件：①二阶平稳（本征）假设、线性估计量——普通克里格②平稳条件不满足，仍采用线性估计量——泛克里格③信息不仅包括二阶矩知识，还包括更多知识（二维分布）——析取克里格线性平稳线性非平稳非线性平稳二、克里格方程组及其方差1.问题的提出设Z(x)为点承载的区域化变量，且是二阶平稳（或本征）的。今要对以x0为中心的盘区V(x0)的平均值（简记为ZV）进行估计。d)(1)(0VVxxZVxZV⊙x0v2⊙x2v1⊙x1v3⊙x3v4⊙x42.线性估计量的构造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)的期望已知时：为简单克里格；未知时：为普通克里格niiiVZZ1*3.简单克里格方程组和简单克里格方差（E(Z(x)=m已知）由于Z(x)的期望为已知，令：Y(x)=Z(x)-m则：E[Y(x)]=E[Z(x)-m]=E[Z(x)]-m=0其协方差函数为：E[Y(x)·Y(y)]=C(x,y)对ZV的估计转化为对YV*的估计，且有：所用的估计量为：只要求得YV的估计值YV*，就能得到ZV的估计值ZV*。mZmxxZVxxYVYVVVVd)(1d)(1),,2,1(1nimZYYYiniiiV其中：显然：YV*是YV的无偏估计量，且不需要任何条件。因为：为了求出λi，使得最小，首先需求出的表达式：所以：（1）其中表示协方差函数在待估域V上的平均值。)(E)(E0d)]([E1)(E0)(E)(E)(E11VVVVniiiniiiVYYxxYVYmZYY22VVEYYE2EninjjijiniViiVVninjjijiniViiVVVVVVVVExxCxxxCVyxyxCVxYxYExxYxYEVyxyYxYEVYYYYYY111211122222),(d),(12dd),(1])()([d)]()([12dd)]()([1][E][E2][EEnininjjijiiiExxCVxCVVC1112),(),(2),(),(VVC为了使达到最小，按照求极值原理，将前述的公式（1）对诸λi求偏导数，并令其为0，则有：（2）于是得到简单克里格方程组：从这个方程组中解出λi（i=1，2，…，n），即为所求的简单克里格系数，它必定满足最小方差无偏估计的要求。将克里格方程组两端均乘以λi，并对i从1到n求和，则有：（3）将（3）式代入公式（1），则得到简单克里格方差的计算公式：（4）0),(2),(212njjijiiExxCVxCninjniiijijiVxCxxC11),(),(2E2EniiiKVxCVVC12),(),(),,2,1(),(),(1niVxCxxCinjjii公式（1）与公式（4）中，所用的估计方差符号不一样，（1）式表示无偏估计量的估计方差，不能保证估计方差最小，故用记号。公式（4）是在确保估计方差最小的前提下推导出来的，它是克里格方差，故记号为。其中关键的区别在于λi（i=1，2，…，n）在两个式中的意义不一样。从克里格方程组解出λi后，即得到YV的简单克里格估计量：所以：2E2KnjjjnjjjKKmZmYmYmZ11)(njjnjjjKmZZ1114.普通克里格方程组和普通克里格方差（E(Z(x)=m未知）要使估计量是无偏的，就必须增加限制条件：（1）无偏性条件若要使ZV*为ZV的无偏估计量，即要求E[ZV*-ZV]=0因为：又因为：所以得无偏性条件为：niiiVZZ1*mxxZEVZEVVd)(1)(niiniiiniiiVmZEZEZE111)()(11nii（2）普通克里格方程组在区域化变量Z(x)满足二阶平稳的条件下类似于简单克里格方法的估计方差的推导，同样可以得到估计方差：在无偏性条件下，要使得估计方差最小，从而求得诸权系数λi,(i=1,2,…,n)，这是一个求条件极值的问题，要用拉格朗日乘数法。令：，为n个权系数λi和μ的(n+1)元函数。-2μ是拉格朗日乘数。求出F对λi,(i=1,2,…,n)以及F对μ的偏导数，并令其为零，得到普通克里格方程组。nininjjijiiiExxCVxCVVC1112),(),(2),(11nii1212niiEF普通克里格方程组：整理得：这n+1个方程的方程组，称为普通克里格方程组。012),,2,1(02),(2),(211niinijiiiiFnixxCVxCF1),,2,1(),(),(11niiinijiiniVxCxxC普通克里格方差：将上式克里格方程组中的第一式（前n个方程）两边乘以λi，再对i从1到n求和得：将此式代入到普通克里格估计方差公式中得：1),,2,1(),(),(11niiinijiiniVxCxxC),(),(111VxCxxCiniininjjiji),(),(12VxCVVCiniiE（3）用变差函数表示的普通克里格方程组与普通克里格方差若Z(x)只满足本征假设，而不满足二阶平稳假设时，则利用协方差函数与变差函数的关系C(h)=C(0)-γ(h)可得用变差函数γ(h)表示的普通克里格方程组与普通克里格方差：),(),(12VVVxniiiE1),,2,1(),(),(11niinjijijniVxxx（4）信息样品为非点承载时的普通克里格方程组与普通克里格方差若样品的承载不能看作是点承载，而是以xi为中心，其体积为vi的承载时，样点之间的协方差C(xi,xj)，就变为样品域之间的平均协方差，相应的普通克里格方程组与普通克里格方差分别写成：用变差函数γ(h)表示的普通克里格方程组与普通克里格方差：niiiKVvCVVC12),(),(1),,2,1(),(),(11niinjijijniVvCvvC),(jivvC),(),(12VVVvniiiK1),,2,1(),(),(11niinjijijniVvvv（5）普通克里格方程组及其方差的矩阵的表示法为简单起见，我们仅给出样品点为非点承载下的普通克里格方程组及其方差的矩阵表示形式：其中：[K]称为普通克里格矩阵，它是一个对称矩阵，因为有：估计方差表示为：1),(),(),(,2121VvCVvCVvCMnnjivvCvvCijji,),(),(MKMVVCTK),(201111),(),(),(1),(),(),(1),(),(),(212221212111nnnnnnvvCvvCvvCvvCvvCvvCvvCvvCvvCK用变差函数表示时，普通克里格方程组的矩阵表示形式为：1),(),(),(,21'21'VvVvVvMnn01111),(),(),(1),(),(),(1),(),(),(212221212111'nnnnnnvvvvvvvvvvvvvvvvvvK'''MK),(''2VVMTK进一步的说明1.只有当协方差矩阵[C(vi,vj)]n×n（即矩阵[K]的左上角n×n阶方阵）是严格正定的，克里格方程组才有唯一解。因为此时其系数矩阵的行列式严格大于零。因此，要求所用的点协方差函数C(h)是正定的。（若用变差函数γ(h)表示，则要求-γ(h)是条件正定的），且数据承载无一重合。因为若有vk=vj，则C(vi,vk)=C(vi,vj)(i=1,2,…,n)，从而矩阵[C(vi,vj)]n×n中有两列（行）完全相等，故其行列式的值为零。进一步的说明2.克里格估值是一种无偏的内插估值。即若待估块段（承载）V与有效数据的任意承载vi重合，则由克里格方程组给出ZK*=Z(vi)及σK2=0。这在制图学中称为“克里格估值曲面通过实测点”。传统的估计方法并没有这种性质。这也说明了克里格估值精度高于其它估值方法。进一步的说明3.对于克里格方程组所用到的协方差函数C(h)和变差函数γ(h)的模型，不论它们所表征的基本结构如何均可，它们可以是各向同性的，也可以是各向异性；既可以是单一结构，也可以是套合结构。进一步的说明4.普通克里格方程组和方差只取决于结构模型C(h)或γ(h)，以及各承载的相对几何特征（或说相对空间位置），而不依赖于数据Zi的具体数值。因此，只要知道结构函数C(h)或γ(h)以及样品的空间位置（数据构形），在开钻前就可得普通克里格方程组及其方差。这样，就可以根据钻孔的空间位置不同，得出不同的克里格方差，从而选择较小的克里格方差所对应的钻孔位置构形，在已知结构函数前提下确定最优的布孔方案。进一步的说明5.普通克里格矩阵[K]，只取决样品承载vi(i=1,2,…,n)的几何特征（空间位置），而完全不依赖于待估块段的承载V。因此，只要所用的信息样品相同，即使对不同的待估块进行估值，克里格方程组的系数矩阵[K]也相同。从而只需求一次逆矩阵[K]-1。