克里格方法(Kriging)完整版

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

2、差之半定义为Z Z（p p）在）在p p方向上的变差函数，记为方向上的变差函数，记为 )()( 2 1 )(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

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

4、最小化非测点 处的估值方差处的估值方差 , ,可推导出克里可推导出克里 格方程组格方程组 )( ii pzz i p 0 p n i ii zz 1 0 n j ijij hchc 1 0) ()( n j j 1 01 02)( 2 00 2 0 zzE 0 p 方程求解后，方程求解后，可得可得 的估值方差的估值方差为为 0 z n i ii hcc 1 0 2 0 )()0( 03 由此可知，估值由此可知，估值 及估值方差及估值方差 完全取决于完全取决于C C（h h） 0 z 2 0 克里格法步骤克里格法步骤 结构分析与变差函数的拟合、运算。结构分析与变差函数的拟合、运算。 利用利用 (

5、h)(h)= C(0) C(h)= C(0) C(h)公式得到公式得到C(h)C(h) 利用克里格方程求出估计量利用克里格方程求出估计量Z(p)Z(p) 01 变差函数：变差函数：几乎所有的几乎所有的变差函数变差函数理论模型都可归纳为以下形式理论模型都可归纳为以下形式 (h)(h)仅取决于测点的样本值仅取决于测点的样本值， (h)(h)则仅取决于测点的空间分布则仅取决于测点的空间分布 02 A(h)A(h)由下式确定：由下式确定：A(h)=A(h)=C(0)C(0) 03 至于至于B(h) B(h) 的参数的参数 利用最大似然法求解，得到利用最大似然法求解，得到 0 )( )( 1 11 k

6、ij n i n j ij hB hB (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)(1-B(h)C(h)=C(0)(1-B(h) 设设 ，则上式可表示为，则上式可表示为 )(1)(hBhce 令令 将上述式子代入克里格方程组可得与将上述式子代入克里格方程组可得与C(0)C(0)无关的克里无关的克里 格方程组和克里格方差，如下格方程组和克里格方差，如下 )()0()(hcchc e e c)0( n j ieejije hchc 1 0 )()( i1,n 1 1 n j j 22

7、0 )0( e c 和和 )(1)0( 1 0 2 0 n i ieie hcc n i ieiee hc 1 0 2 )(1令令 则则 其中，其中， 取决于区域取决于区域上的样本值，上的样本值， 取取决于区域决于区域上测点的空上测点的空 间分布。间分布。上式上式在优化区域在优化区域上测点的空间分布时，只需任意赋予上测点的空间分布时，只需任意赋予C(0) C(0) 一个正数，而无需实际采集的样本值。一个正数，而无需实际采集的样本值。 ) 0 ( c 2 e 上式说明，随机场上估值方差的分布相对大小仅取决于测点的空间分布上式说明，随机场上估值方差的分布相对大小仅取决于测点的空间分布。 )( )(

8、 )( )2( )( )1( k mv k v k v p 将区域网格化，网格单元为边长等于d的正方形；将落在区 域中的m个网格节点依次编号1、2、m,相应的空间坐标 为q 1、q2、qm p 设置区域上n个测点的初始的空间坐标值值 ，取 一变异函数理论模型为B(h),并给c(0)赋一正值 p 假设测点的空间分布调整了k次后，区域中m个网格节点q1、 q2、qm上的估值方差依次为 、 、 ，将这m个估值 方差按由大到小的次序排列，得到 这里，i和 ，且对于任一 ，当 时， 。 )0()0( 2 )0( 1 , n ppp )( 1 k )( 2 k )(k m , 1 )(miv, 1 mjji )()(jviv p 当n个测点的空间分布由 调整为 时，同理可得m个网格节点上的估值方差序列 p 令i=1，判断 是否成立，若成立，则让i=i+1,继续 判断是否成立， p 当 不成立时，分两种情况 情况一： ，表明网格节点上的较大估值方差变小了， 则接受第（k+1）次测点的移动。 情况二： ，表明网格节点上的较大估值方差变大了， 则取消第（k+1）次测点的移动。 )()( 2 )( 1 , k n kk ppp )1()1( 2 )1( 1 , k n kk ppp )1( )( )1( )2( )1( )1 ( k m