克里格方法Kriging

1、n克里格法是地质统计学的核心。克里格法是地质统计学的核心。n解决问题：主要对矿产资源储量进行估计，解决问题：主要对矿产资源储量进行估计，现已推广运用到各领域。现已推广运用到各领域。n方法概要：根据已知样品的空间位置和相关方法概要：根据已知样品的空间位置和相关程度，求出未知区域线性无偏、估计误差最小程度，求出未知区域线性无偏、估计误差最小的储量。的储量。n优点：考虑到样品的空间变异性特征。优点：考虑到样品的空间变异性特征。01变差函数：变差函数：z(p)z(p)为一随机过程，为一随机过程，z(p)z(p)在在p p，p+hp+h两点处的值之差两点处的值之差的方差之半定义为的方差之半定义为z z（

2、p p）在）在p p方向上的变差函数，记为方向上的变差函数，记为)()(21)(hpzpzvarh )(h变差函数描述了区域化变量的变差函数描述了区域化变量的空间结构性。性。 只依赖于只依赖于h h。协方差函数：随机过程协方差函数：随机过程z(p) z(p) 在在p p1 1、p p2 2处的两个随机变量处的两个随机变量z(pz(p1 1) )和和z(pz(p2 2) )的二阶混合中心矩，即的二阶混合中心矩，即covz(pcovz(p1 1), z(p), z(p2 2)=ez(p)=ez(p1 1) )* *z(pz(p2 2)-ez(p)-ez(p1 1)* *ez(pez(p2 2)，记

3、，记为为 c(pc(p1 1, p, p2 2) )02整个区域中，整个区域中，z(p)z(p)的协方差函数存在且相同，即只依赖于的协方差函数存在且相同，即只依赖于h hcovz(p),z(p+h) c(h)；当当h=0h=0时，时，c(0)=varz(x)c(0)=varz(x)， x x03 (h)(h)= c(0) c(h)= c(0) c(h)z(p)z(p)为为区域区域上上随机过程，随机过程，p p； 上上有有n n个测点（样本点），个测点（样本点）， 在在 处的测值处的测值，则，则 处的最优线性估计为处的最优线性估计为 最小化非测点最小化非测点 处的估值方差处的估值方差 , ,可推

4、导出克里可推导出克里格方程组格方程组)(iipzz ip0pniiizz10njijijhchc10)()(njj10102)(20020zze0p方程求解后，方程求解后，可得可得 的估值方差的估值方差为为0 zniiihcc1020)()0(03由此可知，估值由此可知，估值 及估值方差及估值方差 完全取决于完全取决于c c（h h）0 z20克里格法步骤克里格法步骤结构分析与变差函数的拟合、运算。结构分析与变差函数的拟合、运算。利用利用 (h)(h)= c(0) c(h)= c(0) c(h)公式得到公式得到c(h)c(h)利用克里格方程求出估计量利用克里格方程求出估计量z(p)z(p)01

5、变差函数：变差函数：几乎所有的几乎所有的变差函数变差函数理论模型都可归纳为以下形式理论模型都可归纳为以下形式 (h)(h)仅取决于测点的样本值仅取决于测点的样本值， (h)(h)则仅取决于测点的空间分布则仅取决于测点的空间分布02a(h)a(h)由下式确定：由下式确定：a(h)=a(h)=c(0)c(0)03至于至于b(h) b(h) 的参数的参数 利用最大似然法求解，得到利用最大似然法求解，得到0)( )( 111 kijninjijhbhb (h)(h)=a(h)=a(h)* *b(h)b(h)s,21由由 (h)(h)=c(0)b(h)=c(0)b(h)，可得，可得 c(h)=c(0)(

6、1-b(h)c(h)=c(0)(1-b(h)设设 ，则上式可表示为，则上式可表示为)(1)(hbhce令令将上述式子代入克里格方程组可得与将上述式子代入克里格方程组可得与c(0)c(0)无关的克里无关的克里格方程组和克里格方差，如下格方程组和克里格方差，如下)()0()(hcchceec)0(njieejijehchc10)()( i1,n 11njj220)0(ec和和 )(1)0(1020niieiehccniieieehc102)(1令令则则其中，其中， 取决于区域取决于区域上的样本值，上的样本值， 取取决于区域决于区域上测点的空上测点的空间分布。间分布。上式上式在优化区域在优化区域上测

7、点的空间分布时，只需任意赋予上测点的空间分布时，只需任意赋予c(0) c(0) 一个正数，而无需实际采集的样本值。一个正数，而无需实际采集的样本值。) 0 ( c2e上式说明，随机场上估值方差的分布相对大小仅取决于测点的空间分布上式说明，随机场上估值方差的分布相对大小仅取决于测点的空间分布。)()()()2()()1(kmvkvkvp 将区域网格化，网格单元为边长等于d的正方形；将落在区域中的m个网格节点依次编号1、2、m,相应的空间坐标为q 1、q2、qmp 设置区域上n个测点的初始的空间坐标值值 ，取一变异函数理论模型为b(h),并给c(0)赋一正值p 假设测点的空间分布调整了k次后，区域

8、中m个网格节点q1、q2、qm上的估值方差依次为 、 、 ，将这m个估值方差按由大到小的次序排列，得到 这里，i和 ，且对于任一 ，当 时， 。)0()0(2)0(1,nppp)(1k)(2k)(km, 1 )(miv, 1 mjji )()(jvivp 当n个测点的空间分布由 调整为 时，同理可得m个网格节点上的估值方差序列p 令i=1，判断 是否成立，若成立，则让i=i+1,继续判断是否成立，p 当 不成立时，分两种情况 情况一： ，表明网格节点上的较大估值方差变小了，则接受第（k+1）次测点的移动。情况二： ，表明网格节点上的较大估值方差变大了，则取消第（k+1）次测点的移动。 )()(2)(1,knkkppp)1()1(2)1(1,k