1、第三章 克里格插值内容内容 1. 克里格法概述 2. 简单克里格 3. 普通克里格 4. 泛克里格 5. 协同克里格 6. 指示克里格一、概述一、概述1. 克里格法的定义克里格法的定义 矿业定义：根据一个块段（盘区）内外的若干信息样品的矿业定义：根据一个块段（盘区）内外的若干信息样品的某特征值（品位）数据，对该块段（盘区）品位（特征值）作某特征值（品位）数据，对该块段（盘区）品位（特征值）作出一种线性、无偏、最小估计方差的估计方法。出一种线性、无偏、最小估计方差的估计方法。 数学定义：一种求最优、线性、无偏内插估计量的方法。数学定义：一种求最优、线性、无偏内插估计量的方法。 具体说是：考虑了信
2、息样品的形状、大小及其待估块段相互具体说是：考虑了信息样品的形状、大小及其待估块段相互之间的空间分布位置等几何特征，以及变量的空间结构信息后，之间的空间分布位置等几何特征，以及变量的空间结构信息后，为了达到线性、无偏、最小估计方差的估计，而对每个样品值为了达到线性、无偏、最小估计方差的估计，而对每个样品值分别赋予一定的权系数，最后用加权平均来对待估块段的未知分别赋予一定的权系数，最后用加权平均来对待估块段的未知量进行估计的方法。量进行估计的方法。“特定的滑动平均特定的滑动平均”2. 克里格法的种类克里格法的种类（1）线性平稳地质统计学）线性平稳地质统计学 （2）线性非平稳地质统计学）线性非平稳
3、地质统计学 （3）平稳非线性地质统计学）平稳非线性地质统计学 （4）非参数地质统计学）非参数地质统计学简单克里格、普通克里格法简单克里格、普通克里格法泛克里格法和泛克里格法和K K阶本征函数法等阶本征函数法等 条件数学期望、析取克里格法等条件数学期望、析取克里格法等 指示克立格法等指示克立格法等 （5）多元地质统计学、分形地质统计学、时空地质统计学等等。）多元地质统计学、分形地质统计学、时空地质统计学等等。 3. 克里格法的使用信息及应用条件克里格法的使用信息及应用条件 信息：信息： 一组数据；一组数据； 空间构形（坐标）；空间构形（坐标）； 结构信息（变差函数模型）。结构信息（变差函数模型）
4、。 条件：条件： 二阶平稳（本征）假设、线性估计量二阶平稳（本征）假设、线性估计量普通克里格普通克里格 平稳条件不满足，仍采用线性估计量平稳条件不满足，仍采用线性估计量泛克里格泛克里格 信息不仅包括二阶矩知识，还包括更多知识（二维分信息不仅包括二阶矩知识，还包括更多知识（二维分 布）布）析取克里格析取克里格线性线性平稳平稳线性非线性非平稳平稳非线性非线性平稳平稳二、简单克里格法二、简单克里格法 1. 问题的提出问题的提出 设设Z(x)为为点承载点承载区域化变量，且是二阶平稳（或本征）的。今要对以区域化变量，且是二阶平稳（或本征）的。今要对以x0为中心的盘区为中心的盘区V(x0)的平均值的平均值
5、 （简记为（简记为ZV）进行估计。）进行估计。 d)(1)(0VVxxZVxZV x0v2 x2v1 x1v3 x3v4 x4 2. 线性估计量的构造线性估计量的构造 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，又各不相同。，又各不相同。 其线性估计量为：其线性估计量
6、为： 它是它是n个数值的线性组合。个数值的线性组合。 克里格估值的原则克里格估值的原则：就是在保证这个估值：就是在保证这个估值ZV*是无偏的，且估计是无偏的，且估计方差最小的前提下，求出方差最小的前提下，求出n个权系数个权系数i 。在这样的条件下求得的。在这样的条件下求得的i 所构所构成的估计量成的估计量ZV*称为称为ZV的克里格估计量，记为的克里格估计量，记为ZK* 。这时的估计方差。这时的估计方差称为克里格方差，记为称为克里格方差，记为K*。 当当Z (x)的期望已知时：为的期望已知时：为简单克里格简单克里格；未知时：为；未知时：为普通克里格普通克里格 niiiVZZ1* 3. 简单克里格
7、方程组和简单克里格方差简单克里格方程组和简单克里格方差（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 * 。mZmxxZVxxYVYVVVVd)(1d)(1), 2 , 1( 1nimZYYYiniiiV
8、其中： 显然：显然： YV*是是YV的无偏估计量，且不需要任何条件。因为：的无偏估计量，且不需要任何条件。因为： 为了求出为了求出i，使得使得 最小，首先需求出最小，首先需求出 的表达式：的表达式： 所以：所以： （1） 其中其中 表示协方差函数在待估域表示协方差函数在待估域V上的平均值。上的平均值。)(E)(E0d)(E1)(E0)(E)(E)(E11VVVVniiiniiiVYYxxYVYmZYY22VVEYYE2E ninjjijiniViiVVninjjijiniViiVVVVVVVVExxCxxxCVyxyxCVxYxYExxYxYEVyxyYxYEVYYYYYY1112111222
9、22),(d),(12dd),(1 )()( d)()(12dd)()(1 EE2EEnininjjijiiiExxCVxCVVC1112),(),(2),(),(VVC 为了使为了使 达到最小，按照求极值原理，将前述的达到最小，按照求极值原理，将前述的 公式（公式（1）对）对诸诸i求偏导数，并令其为求偏导数，并令其为0，则有：，则有： （2） 于是得到简单克里格方程组：于是得到简单克里格方程组： 从这个方程组中解出从这个方程组中解出i （i=1，2， ，n），即为所求的简单克里），即为所求的简单克里格系数，它必定满足最小方差无偏估计的要求。格系数，它必定满足最小方差无偏估计的要求。 将克里格
10、方程组两端均乘以将克里格方程组两端均乘以i ，并对，并对i 从从1到到n求和，则有：求和，则有： （3） 将（将（3）式代入公式（）式代入公式（1），则得到简单克里格方差的计算公式：），则得到简单克里格方差的计算公式： （4）0),(2 ),(212njjijiiExxCVxCninjniiijijiVxCxxC11),(),(2E2EniiiKVxCVVC12),(),(), 2 , 1( ),(),(1niVxCxxCinjjii 公式（公式（1）与公式（）与公式（4）中，所用的估计方差符号不一样，（）中，所用的估计方差符号不一样，（1）式表）式表示无偏估计量的估计方差，不能保证估计方差最
11、小，故用记号示无偏估计量的估计方差，不能保证估计方差最小，故用记号 。公。公式（式（4）是在确保估计方差最小的前提下推导出来的，它是克里格方差，）是在确保估计方差最小的前提下推导出来的，它是克里格方差，故记号为故记号为 。其中关键的区别在于。其中关键的区别在于i （i=1，2， ，n）在两个式中的）在两个式中的意义不一样。意义不一样。 从克里格方程组解出从克里格方程组解出i 后，即得到后，即得到YV的简单克里格估计量：的简单克里格估计量： 所以：所以：2E2KnjjjnjjjKKmZmYmYmZ11)(njjnjjjKmZZ111三、普通克里格法三、普通克里格法1. 普通克里格方程组和普通克里
12、格方差普通克里格方程组和普通克里格方差（E(Z(x)=m 未知）未知） 要使估计量要使估计量 是无偏的，就必须增加限制条件：是无偏的，就必须增加限制条件： （1）无偏性条件）无偏性条件 若要使若要使ZV*为为ZV的无偏估计量，即要求的无偏估计量，即要求 EZV*- ZV =0 因为：因为： 又因为：又因为： 所以得无偏性条件为：所以得无偏性条件为：niiiVZZ1*mxxZEVZEVVd)(1)(niiniiiniiiVmZEZEZE111)()(11nii （2）普通克里格方程组）普通克里格方程组 在区域化变量在区域化变量Z(x) 满足二阶平稳的条件下类似于简单克里格方法满足二阶平稳的条件下
13、类似于简单克里格方法的估计方差的推导，同样可以得到估计方差：的估计方差的推导，同样可以得到估计方差： 在无偏性条件在无偏性条件 下，要使得估计方差最小，从而求得诸权下，要使得估计方差最小，从而求得诸权系数系数 i , (i=1,2,n)，这是一个求条件极值的问题，要用拉格朗日乘这是一个求条件极值的问题，要用拉格朗日乘数法。数法。 令：令： ，为，为n个权系数个权系数 i和和 的的(n+1)元函数。元函数。-2 是拉格朗日乘数。求出是拉格朗日乘数。求出F对对 i , (i=1,2,n)以及以及F对对的偏导数，并的偏导数，并令其为零，得到普通克里格方程组。令其为零，得到普通克里格方程组。ninin
14、jjijiiiExxCVxCVVC1112),(),(2),(11nii1212niiEF普通克里格方程组：普通克里格方程组： 整理得：整理得： 这这n+1个方程的方程组，称为普通克里格方程组。个方程的方程组，称为普通克里格方程组。 012), 2 , 1( 02),(2),(211niinijiiiiFnixxCVxCF1), 2 , 1( ),(),(11niiinijiiniVxCxxC普通克里格方差：普通克里格方差： 将上式克里格方程组中的第一式（前将上式克里格方程组中的第一式（前n个方程）两边乘以个方程）两边乘以 i ，再，再对对i 从从1到到n求和得：求和得： 将此式代入到普通克里
15、格估计方差公式中得：将此式代入到普通克里格估计方差公式中得：1), 2 , 1( ),(),(11niiinijiiniVxCxxC),(),(111VxCxxCiniininjjiji),(),(12VxCVVCiniiE （3）用变差函数表示的普通克里格方程组与普通克里格方差）用变差函数表示的普通克里格方程组与普通克里格方差 若若Z(x)只满足本征假设，而不满足二阶平稳假设时，则利用协方只满足本征假设，而不满足二阶平稳假设时，则利用协方差函数与变差函数的关系差函数与变差函数的关系C(h)=C(0) - (h) 可得用变差函数可得用变差函数(h)表示的表示的普通克里格方程组与普通克里格方差：
16、普通克里格方程组与普通克里格方差：),(),(12VVVxniiiE1), 2 , 1( ),(),(11niinjijijniVxxx 用变差函数表示时，普通克里格方程组的矩阵表示形式为：用变差函数表示时，普通克里格方程组的矩阵表示形式为： 1),(),(),( ,2121VvVvVvMnn 01111),(),(),(1),(),(),(1),(),(),(212221212111nnnnnnvvvvvvvvvvvvvvvvvvK MK ),( 2VVMTK （4）信息样品为）信息样品为非点承载非点承载时的普通克里格方程组与普通克里格方差时的普通克里格方程组与普通克里格方差 若样品的承载不
17、能看作是点承载，而是以若样品的承载不能看作是点承载，而是以x i为中心，其体积为为中心，其体积为v i的承载时，样点之间的协方差的承载时，样点之间的协方差C(xi ,x j )，就变为样品域之间的平均协，就变为样品域之间的平均协方差方差 ，相应的普通克里格方程组与普通克里格方差分别写成：，相应的普通克里格方程组与普通克里格方差分别写成： 用变差函数用变差函数(h)表示的普通克里格方程组与普通克里格方差：表示的普通克里格方程组与普通克里格方差：niiiKVvCVVC12),(),(1), 2 , 1( ),(),(11niinjijijniVvCvvC),(jivvC),(),(12VVVvni
18、iiK1), 2 , 1( ),(),(11niinjijijniVvvv （5）普通克里格方程组及其方差的矩阵的表示法）普通克里格方程组及其方差的矩阵的表示法 为简单起见，我们仅给出样品点为非点承载下的普通克里格方程为简单起见，我们仅给出样品点为非点承载下的普通克里格方程组及其方差的矩阵表示形式：组及其方差的矩阵表示形式： 其中：其中： K称为普通克里格矩阵，它是一个对称矩阵，因为有：称为普通克里格矩阵，它是一个对称矩阵，因为有： 估计方差表示为：估计方差表示为： 1),(),(),( ,2121VvCVvCVvCMnnjivvCvvCijji , ),(),( MK MVVCTK),(2
19、01111),(),(),(1),(),(),(1),(),(),(212221212111nnnnnnvvCvvCvvCvvCvvCvvCvvCvvCvvCK 用变差函数表示时，普通克里格方程组的矩阵表示形式为：用变差函数表示时，普通克里格方程组的矩阵表示形式为： 1),(),(),( ,2121VvVvVvMnn01111),(),(),(1),(),(),(1),(),(),(212221212111nnnnnnvvvvvvvvvvvvvvvvvvK MK ),( 2VVMTK进一步的说明进一步的说明 1. 只有当协方差矩阵只有当协方差矩阵C(vi,vj)nn（即矩阵（即矩阵K的左上角的
20、左上角nn阶方阵）是严格正定的，克里格方程组才有唯一解。阶方阵）是严格正定的，克里格方程组才有唯一解。因为此时其系数矩阵的行列式严格大于零。因此，要求所因为此时其系数矩阵的行列式严格大于零。因此，要求所用的点协方差函数用的点协方差函数C(h)是正定的。（若用变差函数是正定的。（若用变差函数(h)表表示，则要求示，则要求- (h)是条件正定的），且数据承载无一重合。是条件正定的），且数据承载无一重合。因为若有因为若有vk=vj，则，则C(vi,vk)= C(vi,vj) (i=1,2,n)，从而矩阵，从而矩阵C(vi,vj)nn中有两列（行）完全相等，故其行列式的值为中有两列（行）完全相等，故其
21、行列式的值为零。零。进一步的说明进一步的说明 2. 克里格估值是一种无偏的内插估值。即若待估块段克里格估值是一种无偏的内插估值。即若待估块段（承载）（承载）V与有效数据的任意承载与有效数据的任意承载vi重合，则由克里格方程重合，则由克里格方程组给出组给出ZK*=Z(vi)及及K2=0。这在制图学中称为这在制图学中称为“克里格估克里格估值曲面通过实测点值曲面通过实测点”。传统的估计方法并没有这种性质。传统的估计方法并没有这种性质。这也说明了克里格估值精度高于其它估值方法。这也说明了克里格估值精度高于其它估值方法。进一步的说明进一步的说明 3. 对于克里格方程组所用到的协方差函数对于克里格方程组所
22、用到的协方差函数C(h)和变差和变差函数函数(h)的模型，不论它们所表征的基本结构如何均可，的模型，不论它们所表征的基本结构如何均可，它们可以是各向同性的，也可以是各向异性；既可以是单它们可以是各向同性的，也可以是各向异性；既可以是单一结构，也可以是套合结构。一结构，也可以是套合结构。进一步的说明进一步的说明 4. 普通克里格方程组和方差只取决于结构模型普通克里格方程组和方差只取决于结构模型C(h)或或(h) ，以及，以及各承载的相对几何特征（或说相对空间位置），而不依赖于数据各承载的相对几何特征（或说相对空间位置），而不依赖于数据Zi 的的具体数值。具体数值。因此，只要知道结构函数因此，只要
23、知道结构函数C(h)或或(h)以及样品的空间位置以及样品的空间位置（数据构形），在开钻前就可得普通克里格方程组及其方差。这样，（数据构形），在开钻前就可得普通克里格方程组及其方差。这样，就可以根据钻孔的空间位置不同，得出不同的克里格方差，从而选择就可以根据钻孔的空间位置不同，得出不同的克里格方差，从而选择较小的克里格方差所对应的钻孔位置构形，较小的克里格方差所对应的钻孔位置构形，在已知结构函数前提下确在已知结构函数前提下确定最优的布孔方案。定最优的布孔方案。进一步的说明进一步的说明 5. 普通克里格矩阵普通克里格矩阵K ，只取决样品承载，只取决样品承载vi (i=1,2,n)的几何特征的几何特
24、征（空间位置），而完全不依赖于待估块段的承载（空间位置），而完全不依赖于待估块段的承载V。因此，只要所用。因此，只要所用的信息样品相同，即使对不同的待估块进行估值，克里格方程组的系的信息样品相同，即使对不同的待估块进行估值，克里格方程组的系数矩阵数矩阵K 也相同。从而只需求一次逆矩阵也相同。从而只需求一次逆矩阵K -1。若估计构形（待估。若估计构形（待估承载与全体样品承载的构形）也相同，则矩阵承载与全体样品承载的构形）也相同，则矩阵M也不变。即只需解也不变。即只需解一次克里格方程组，就可得到线性估计量中的权系数一次克里格方程组，就可得到线性估计量中的权系数i (i=1,2,n)，大大地节省计算
25、时间。（规则勘探网格就满足这一要求）大大地节省计算时间。（规则勘探网格就满足这一要求）进一步的说明进一步的说明 6. 普通克里格方程组及其方差考虑了以下四个方面的因素：普通克里格方程组及其方差考虑了以下四个方面的因素： （1）待估承载）待估承载V 的几何特征（的几何特征（ (V, V) ）；）； （2）数据构形的几何特征（）数据构形的几何特征（ (vi, vj) ）；）； （3）信息样品承载）信息样品承载vi 与待估承载与待估承载V之间的距离（之间的距离（ (vi, V) ） ； （4）反应区域化变量）反应区域化变量Z(x)空间结构特征的变差函数模型（空间结构特征的变差函数模型（ (h) ）
26、。进一步的说明进一步的说明 7. 纯块金效应对普通克里格方程组及其方差的影响纯块金效应对普通克里格方程组及其方差的影响 如果原来变差函数如果原来变差函数1(h) ，后来增加了一个块金常数，后来增加了一个块金常数C0成为：成为：则当所有信息样品承载则当所有信息样品承载vi (i=1,2,n)大小相等，和待估块段大小相等，和待估块段V 彼此都不相彼此都不相交，且交，且V比比vi 大很多（这些条件在实际中常能被满足）时，块金效应（即大很多（这些条件在实际中常能被满足）时，块金效应（即增加了一个块金常数）对普通克里格方程组的影响只是在原普通克里格增加了一个块金常数）对普通克里格方程组的影响只是在原普通
27、克里格矩阵的主对角线上前矩阵的主对角线上前n个元素中减去块金常数。这时方程组变为：个元素中减去块金常数。这时方程组变为：0 0 0)()(01hChhh1),(),(),(),(),(),(),(),(),(),(),(),(111021211121212210212111111211211101niinnnnnnnnnnVvvvCvvvvVvvvvvCvvVvvvvvvvC普通克里格的计算普通克里格的计算克里格计算通常可分两种：克里格计算通常可分两种：1 1、点克里格法、点克里格法: : 若用它们对某一点若用它们对某一点x x0 0处进行克里格估值，处进行克里格估值，则叫点克里格法则叫点克里
28、格法2 2、块段克里格法：若用它们对以某一点、块段克里格法：若用它们对以某一点x x0 0为中心的某一块为中心的某一块段段V V进行克里格估值，则叫块段克里格法进行克里格估值，则叫块段克里格法例例1-点克里格法点克里格法例例2-块段克里格法块段克里格法关于权系数的一点说明关于权系数的一点说明1、克里格权系数从理论上讲可正可负，会导致估值也出现负值。 对于煤层厚度、品位等一些只能取正值的数据，需要进行一些处理： 方法一：人为改变权系数，将负权系数规定为0，其他权系数按比例调整。 方法二：人为改变其估值，权系数不变，当估值出现负数时，规定为0.克里格权系数的特点克里格权系数的特点2、对称性：当区域
29、化变量具有各向同性结构时，在几何形状、位置上待估计承载V具有对称性的信息样品vi,vj，在没有其他效应影响的条件下，具有相同的权系数。 克里格权系数的特点克里格权系数的特点3、减弱丛聚效应（declustering effect）:用克里格法估计时，不会因为有一些样品丛聚在一起而过分增大它们总的权系数 克里格权系数的特点克里格权系数的特点4、屏蔽效应：当块金常数为0或很小时，在待估计块段内部的样品的克里格权值最大，块段周围第一圈样品的权系数也较大，第二圈系数明显减小。这就是克里格权系数的屏蔽效应。随着块金效应的增大，屏蔽效应变弱 四、泛克里格法（四、泛克里格法（UK)普通克立格法，要求区域化变
30、量普通克立格法，要求区域化变量Z(x)是二阶平稳的，至少是准平稳或准内是二阶平稳的，至少是准平稳或准内蕴的，即至少在估计邻域内有蕴的，即至少在估计邻域内有EZ(x) =m (常数常数)成立。成立。然而，在生产实践中，许多区域化变量然而，在生产实践中，许多区域化变量Z（x）在研究区是非平稳的即）在研究区是非平稳的即EZ(x) =m (x)这时就不能用普通克立格法进行估计了。这时就不能用普通克立格法进行估计了。泛克立格法：在漂移的形式泛克立格法：在漂移的形式 EZ(x)=m (x) 和非平稳随机函数和非平稳随机函数Z(x)的协方差的协方差C(h)或变异函数或变异函数(h) 为已知的条件下，一为已知
31、的条件下，一种考虑到有漂移的无偏、线性估计量的地质统计学方法种考虑到有漂移的无偏、线性估计量的地质统计学方法(也称为也称为K阶无偏克阶无偏克立格法立格法)。地质工作中的两种现象地质工作中的两种现象（1）在一个区域中，某区域化变）在一个区域中，某区域化变量量(如金属含量如金属含量)自西向东逐渐升高，自西向东逐渐升高，即所谓有即所谓有“漂移漂移”(drift)存在存在(可以可以理解为理解为“趋势趋势”） ，但在一个小的，但在一个小的局部范围内，又可以认为是平稳的，局部范围内，又可以认为是平稳的，普通克立格法正是利用了这个性质普通克立格法正是利用了这个性质进行估计的。进行估计的。（2）从整体上看，该
32、区的某区域化）从整体上看，该区的某区域化变量是平稳的，但在一个小的局部范变量是平稳的，但在一个小的局部范围，却呈现出存在有漂移的现象，即围，却呈现出存在有漂移的现象，即具有非平稳的数学期望具有非平稳的数学期望 EZ(x) =m (x) 。这时就不能用普通。这时就不能用普通克立格法进行估计。换句话说，当估克立格法进行估计。换句话说，当估计邻域内的有效数据不足以应用普通计邻域内的有效数据不足以应用普通克立格法进行准平稳估计时，就必须克立格法进行准平稳估计时，就必须考虑到漂移的存在。考虑到漂移的存在。对于非平稳的对于非平稳的(或说有漂移的或说有漂移的)的区域化变最的区域化变最z(x)，假设它们可以分
33、解，假设它们可以分解为漂移为漂移(drift)和涨落（和涨落（fluctation）两部分，即）两部分，即 Z(x)=m(x)+R(x)其中，其中， m (x)=EZ(x) 为在点为在点x处处Z(x)的数学期望，称为漂移，的数学期望，称为漂移，R(x)称称为涨落（也称剩余、波动）为涨落（也称剩余、波动） ，它是一个数学期望为零的区域化变量，它是一个数学期望为零的区域化变量.漂移的形式漂移的形式-通常采用多项式形式通常采用多项式形式漂移还有次数的确定漂移还有次数的确定:(1)由于对某一位置品位的估计往往只是根据其有限邻域内的数据来进行的，由于对某一位置品位的估计往往只是根据其有限邻域内的数据来进
34、行的，为了表示该有限邻域内的漂移一般只需用一次或二次多项式就足够了。若为了表示该有限邻域内的漂移一般只需用一次或二次多项式就足够了。若领域小、数据多、递变速度小，则可用一次多项式；反之，则用二次多项式。领域小、数据多、递变速度小，则可用一次多项式；反之，则用二次多项式。在一般情况下，实际矿业工作中在一般情况下，实际矿业工作中90%的漂移问题都可采用线性漂移来解决。的漂移问题都可采用线性漂移来解决。(2)漂移的多项式也不必采用完全的多项式，有时采用局部漂移漂移的多项式也不必采用完全的多项式，有时采用局部漂移(即不完全的多即不完全的多项式项式)就可以了。如只有水平方向的漂移，而无垂向的漂移。就可以
35、了。如只有水平方向的漂移，而无垂向的漂移。（充分考虑地质特征，这样才能准确判断在哪些方向上存在着漂移充分考虑地质特征，这样才能准确判断在哪些方向上存在着漂移）在给定的以点在给定的以点x0为中心的邻域为中心的邻域V内，有任一点内，有任一点xV，其漂移，其漂移m(x)可用多项可用多项式函数来近似表达。式函数来近似表达。（一）已知协方差时，（一）已知协方差时，Z(x)的泛克里格估计的泛克里格估计假设：假设：设设 是非平稳区域化变量，且有是非平稳区域化变量，且有设设 估计量为估计量为klllxfaxmxmxZE0)()()()(),()()()()(yxCymxmyZxZE)(xZ*UKZ)(xZ已知
36、已知n个样品点个样品点 (a=1,2,3,n)，样品值为，样品值为 (a=1,2,3,n) ， 现现在要用这些样品点估计邻域中任一点在要用这些样品点估计邻域中任一点x的品位值的品位值Z(x)，Z(x)可以表为可以表为这这n个个样品值的线性组合样品值的线性组合:系数系数应该遵循无偏性和最优性（估计方差极小）来确定。应该遵循无偏性和最优性（估计方差极小）来确定。xZxZ)(nUKZZ1* x2 x1 x0 x3 x41. 无偏性条件无偏性条件式中：式中： )()(*xmxZEZEUKklxfxfnll.2 , 1 , 0)()(12.最优性条件：最优性条件： 满足无偏性条件下，用满足无偏性条件下，
37、用 估计估计Z(x)的方差是的方差是它必须极小。它必须极小。 *注：推导过程注：推导过程 *UKZnnnUKExxxxCxxCxZZE1112*2),(2),(),()(3.泛克里格方程组：泛克里格方程组：在无偏性条件下使估计方差达到极小。就要利用求条件极值的拉格朗日乘数法，在无偏性条件下使估计方差达到极小。就要利用求条件极值的拉格朗日乘数法，即要求使即要求使达到极小时的达到极小时的泛克立格方程组：泛克立格方程组：上式有上式有n+k+1个未知数个未知数 可由可由n+k+1个方程个方程解得。解得。相应的泛克立格方差：相应的泛克立格方差：kllalnaalExfxfF012)()(2la及klxf
38、xfnxxCxfxxCnllkllln, 2 , 1 , 0)()(, 2 , 1),()(),(101), 2 , 1 , 0;, 2 , 1,(klnlklllnukxfxxCxxC012)(),(),(已知已知n个样品点个样品点 (a=1,2,3,n)，样品值为，样品值为 (a=1,2,3,n) ， 现现在要用这些样品点估计邻域中任一点在要用这些样品点估计邻域中任一点x的的漂移值漂移值m(x)，设设m(x)的估计量为的估计量为m*(x), 它是这它是这n个样品值的线性组合个样品值的线性组合:式中式中应在无偏性和最小估计方差下求得。应在无偏性和最小估计方差下求得。xZxZ)(（二）已知协方
39、差时，（二）已知协方差时，m(x)的泛克里格估计的泛克里格估计naaaZxm1*)(1. 无偏性条件无偏性条件式中：式中： )()(*xmxmE).3 , 2 , 1 , 0()()(1klxfxfnalala2.最优性条件：最优性条件： 满足无偏性条件下，用满足无偏性条件下，用 m*(x) 估计估计m(x)的方差是的方差是 *注：推导过程注：推导过程 nanaExxaC112),(3.泛克里格方程组：泛克里格方程组：在无偏性条件下，求使在无偏性条件下，求使达到极小时的达到极小时的泛克立格方程组：泛克立格方程组：相应的泛克立格方差：相应的泛克立格方差：kllalnaalExfxfF012)()
40、(2la及klxfxfnxfxxCnllkllln, 2 , 1 , 0)()(, 2 , 10)(),(101klllukVfxm02)()( 若欲根据以若欲根据以xa为中心的承载为中心的承载v（xa）=va （a=1,2,n)上的上的n个信息样品的平均品个信息样品的平均品位位 估计以估计以x为中心的块段为中心的块段V(x) 上的平均品位：上的平均品位：则类似地可得到则类似地可得到Zv(x)的泛克里格方程组：的泛克里格方程组：其中：其中：klbbnVvCbvvCnlVlvakllvlnaa, 2 , 1 , 0, 2 , 1),(),(101)(u)du1)(aaxvavvZvxZZ)(u)
41、du1xVvZVZVlvlVlalvdxxfVbdxxfvbaa)(1,)(1V x0v2 x2v1 x1v3 x3v4 x4泛克里格方程组（块估计）泛克里格方程组（块估计）Z(x)的泛克立格方程组：的泛克立格方程组：m(x)的泛克立格方程组：的泛克立格方程组：Z（x）与）与m（x）泛克里格方程组的比较）泛克里格方程组的比较klxfxfnxxCxfxxCnllkllln, 2 , 1 , 0)()(, 2 , 1),()(),(101klxfxfnxfxxCnllkllln, 2 , 1 , 0)()(, 2 , 10)(),(101(1)趋势面分析由于没考虑空趋势面分析由于没考虑空间间数据的
42、结构函数特点，且是对全区数据整体数据的结构函数特点，且是对全区数据整体进行估计的，对进行估计的，对m(x)的的估计也不是最优的估计也不是最优的；而泛克立格法考虑了空间数据而泛克立格法考虑了空间数据的结构函数特点，又采取在邻域内进行局部估计的方的结构函数特点，又采取在邻域内进行局部估计的方j击。故对击。故对m(x)的估计的估计是最优的。是最优的。(2)趋势趋势面面分析虽也能估计分析虽也能估计“趋势趋势”(理解为理解为“漂移漂移”的的“趋势趋势”)m(x)，却不，却不能能估计估计x点处点处Z(x)的其实值，也不能估计以的其实值，也不能估计以x为中心的盘为中心的盘区区V上的平均值上的平均值Zv(x)
43、。而泛克立格法则对此三种值均能估计。而泛克立格法则对此三种值均能估计。(3)趋势面分析无法给出估计趋势面分析无法给出估计m(x) 的误差。的误差。（三）趋势面分析的缺陷（三）趋势面分析的缺陷若若泛克里格方程组泛克里格方程组普通克里格方程组普通克里格方程组普通克里格是泛克里格的特例普通克里格是泛克里格的特例)()()(0常数klllmxfaxmklxfxfnxxCxfxxCnllkllln, 2 , 1 , 0)()(, 2 , 1),()(),(101nnnxxCxxC1011, 2 , 1),(),(对化探数据的处理常用趋势面分析、滑动平均等方法，但都存在一些缺点对化探数据的处理常用趋势面分
44、析、滑动平均等方法，但都存在一些缺点和不足如趋势面次数的确定带有较大的人为性；滑动平均窗口大小的选和不足如趋势面次数的确定带有较大的人为性；滑动平均窗口大小的选择缺少理论依据；而建立在区域化变量理论基础上的地质统计学方法则能择缺少理论依据；而建立在区域化变量理论基础上的地质统计学方法则能弥补这些不足，具体地说，就是可以用泛克立格法或对数泛克立格法来处弥补这些不足，具体地说，就是可以用泛克立格法或对数泛克立格法来处理化探数据。理化探数据。应用实例应用实例 实例：阳新地区实例：阳新地区Cu元素元素l：5万岩石化学测量数据万岩石化学测量数据的泛克里格处理的泛克里格处理接触交代矿床：以铜为主，铁、钼，
45、金、银、硫等皆有产出。接触交代矿床：以铜为主，铁、钼，金、银、硫等皆有产出。（1）数据的统计特征）数据的统计特征数据统计分析数据统计分析（对数（对数直方图直方图）的结果表明，它具有混合总体分布的特征的结果表明，它具有混合总体分布的特征：低值总体代表边缘地区沉积岩中低值总体代表边缘地区沉积岩中Cu元素的分布特征元素的分布特征，高值总体代表本区高值总体代表本区广泛分布的石英闪长岩体中广泛分布的石英闪长岩体中cu元素的分布特征。元素的分布特征。一般情况下，由于化探数据往往取自较大一般情况下，由于化探数据往往取自较大的区域范围的区域范围，不同的地质背景必然会影响不同的地质背景必然会影响变量的分布特征。
46、当有较明显的多峰分布变量的分布特征。当有较明显的多峰分布时，应先对多峰总体筛分后，再分别处理。时，应先对多峰总体筛分后，再分别处理。经筛分后，石英闪长岩体中经筛分后，石英闪长岩体中cu元素含量数元素含量数据近似对据近似对数数正态分布律可对该批数据用正态分布律可对该批数据用对数泛克立格法处理。对数泛克立格法处理。(2)结构分析结构分析 将原始数据取对数后计算出对数变差图将原始数据取对数后计算出对数变差图，呈现出带状异向性，带状方向呈现出带状异向性，带状方向(南北方向南北方向)上存在着线性漂移。这是因为石英闪长岩体在近东西上存在着线性漂移。这是因为石英闪长岩体在近东西(x)方向上方向上Cu元素含量
47、值变化较为稳定、连续，变化幅度较元素含量值变化较为稳定、连续，变化幅度较小小，而近南北，而近南北(y)方向上变方向上变化幅度较大，并存在线性漂移。化幅度较大，并存在线性漂移。Cu元素空问分布的这种东向带状延伸是受元素空问分布的这种东向带状延伸是受岩体及区域性构造控制的岩体及区域性构造控制的。对数泛克立格估值对数泛克立格估值图图不但能较好地反映不但能较好地反映Cu元素在区域上的变化趋势，同时元素在区域上的变化趋势，同时也能较好地反映局部变化特征也能较好地反映局部变化特征。剩余剩余图图上的异常位置、形态和规模与已有矿床、矿点相吻合，而且能为寻上的异常位置、形态和规模与已有矿床、矿点相吻合，而且能为
48、寻找新的化探异常提供信息。找新的化探异常提供信息。五、协同克里格法（五、协同克里格法（CK)概述概述在多金属矿床中，常遇到多个变量协同区域化的现象，如银、铜，铅，在多金属矿床中，常遇到多个变量协同区域化的现象，如银、铜，铅，锌等几个变量的协同区域化，可用锌等几个变量的协同区域化，可用K个区域化变量的集合来表示，即个区域化变量的集合来表示，即Zk（x）。当某一个变量的取样量不足以获得所需精度韵估计量，而其他）。当某一个变量的取样量不足以获得所需精度韵估计量，而其他变量却有较充足的取样量时，研究这个变量与其他变量间的空间相关关变量却有较充足的取样量时，研究这个变量与其他变量间的空间相关关系，借助其
49、他变量的样品信息用协同克立格法就可以提高对这个变量的系，借助其他变量的样品信息用协同克立格法就可以提高对这个变量的估计精度。估计精度。例如，对一个盘区中银的平均品位的估计，除了利用银的品位数据之外，例如，对一个盘区中银的平均品位的估计，除了利用银的品位数据之外，还可以利用与它有协同区域他关系的铜、铅、锌的品位数据。还可以利用与它有协同区域他关系的铜、铅、锌的品位数据。 从理论上来讲，协同克立格法与普通克立格法很类似，没有本质差别，从理论上来讲，协同克立格法与普通克立格法很类似，没有本质差别，主要就是多用了一个下标主要就是多用了一个下标k来表示各个变量。来表示各个变量。协同克立格法（协同克立格法
50、（Cokriging）是以协同区域化理论）是以协同区域化理论为基础，以变异函数为基本工具，研究定义于同一为基础，以变异函数为基本工具，研究定义于同一空间域中，既有统计相关，又有空间位置相关的多空间域中，既有统计相关，又有空间位置相关的多元信息的空间结构的自然现象的科学。元信息的空间结构的自然现象的科学。 KknkkkkkVZZ11*0在无偏、最优条件下得到协同克里格方程组在无偏、最优条件下得到协同克里格方程组协同克里格方差：协同克里格方差：0111,101, 2 , 1;, 2 , 1),()(000;00kkKknvVCvvCkkkkkkkkkkkknnkkKkkkkknkk协同克里格方程组
