Kriging方法公式推导

1、地统计（Geostatistics）又称地质统计，它是以区域化变量为基础，借助变异函数，研究既具有随机性又具有结构性，或空间相关性和依赖性的自然现象的一门科学。凡是与空间数据的结构性和随机性，或空间相关性和依赖性，或空间格局与变异有关的研究，并对这些数据进行最优无偏内插估计，或模拟这些数据的离散性、波动性时，皆可应用地统计学的理论与方法。地统计方法的条件：(1)随机过程 (2)正态分布 （3)平稳性： 均值平稳 协方差函数有关的二阶平稳和与半变异函数有关的内蕴平稳 Kriging 插值的条件：变异函数和相关分析的结果表明区域化变量存在空间相关性。 首先假设区域化变量 满足二阶平稳假设和本征假设

2、，其数学期望为m，协方差函数 及变异函数 存在。 即 假设在待估计点（x）的临域内共有n个实测点，即x1，x2，xn，其样本值为 。那么，普通克里格法的插值公式为 )(xZc(h)(hmxZE)(2)()()(mhxZxZEhc2)()(21)(hxZxZEhniiixZxZ1*)()()(ixZ(4.2.22) 其中 为权重系数，表示各空间样本点 处的观测值 对估计值 的贡献程度。 可见，克立格插值的关键就是计算权重系数 。显然，权重系数的求取必须满足两个条件： 一是使 的估计是无偏的，即偏差的数学期望为零； 二是最优的，即使估计值 和实际值 之差的平方和最小。 为此，需要满足以下两个条件：

3、iix)(ixZ)(*xZi)(*xZ)(*xZ)(ixZ (1)无偏性。要使 成为 的无偏估计量，即 当 时，也就是当 时，则有 这时， 为 的无偏估计量。 （2）最优性。在满足无偏性条件下，估计方差为mxZExZEniiiniii11)()()()(*xZExZE)(*xZ)(ixZmxZE)(nii11)(*xZ)(ixZ212*2 )()()()(niiiExZxZExZxZE 使用协方差函数表达，它可以进一步写为 (4.2.24) 为使估计方差最小，根据拉格朗日乘数原理，令 (4.2.25) 求F对 和 的偏导数，并令其为0，得克立格方程组 (4.2.26) ninjniiijiji

4、Exxcxxcxxc1112),(2),(),() 1(212niiEFiniiijinjjiFxxcxxcF110) 1(202),(2),(2 niiiExxcxxc12),(),(niiijinjjxxcxxc111),(),(4.2.27)(4.2.28)整理后得 解线性方程组(4.2.27)式,求出权重系数i和拉格朗日系数,代入公式(4.2.24),可得克立格估计方差上述过程也可用矩阵形式表示，令 则普通克立格方程组为 其估计方差为 1),(),(),(,01111112121212222111211xxcxxcxxcDcccccccccKnnnnnnnnDKDK1DxxcTK),(

5、2 在变异函数存在的条件下，根据协方差与变异函数的关系： ，也可以用变异函数表示普通克立格方程组和克立格估计方差，即 (4.2.29) 解线性方程组(4.2.27)式，求出权重系数 和拉格朗日乘数，代入公式(4.2.24)，可得克立格估计方差 ，即 )()0()(hcchniiijinjjxxxx111),(),(),(),(12xxxxniiiK(4.2.30) 也可以将克立格方程组和估计方差用变异函数写成上述矩阵形式。令 在以上的介绍中，区域化变量 的数学期望 可以是已知或未知的。如果m是已知常数，称为简单克立格法；如果m是未知常数，称为普通克立格法。不管是哪一种方法，均可根据方法计算权重

6、系数和克立格估计量。 1),(),(),(,01111112121212222111211xxxxxxDKnnnnnnnnDKDK1),(2xxDTK)(xZmxZE)(4.2.34)(4.2.35) (4.2.36)克里格插值基础克里格方法概述 克里格方法（Kriging）又称空间局部插值法，是以变异函数理论和结构分析为基础，在有限区域内对区域化变量进行无偏最优估计的一种方法，是地统计学的主要内容之一。其实质是利用区域化变量的原始数据和变异函数的结构特点，对未知样点进行线性无偏、最优估计。无偏是指偏差的数学期望为0，最优是指估计值与实际值之差的平方和最小。 克里格插值基础2.克里格方法的具体

7、步骤导入数据数据分析是否服从正态分布是否是否存在趋势否是数据变换泛克里格方法根据数据选择合适的方法计算样点间的距离矩阵计算样点间的属性方差按距离分组按组统计平均距离及对应的平均方差绘制方差变异云图绘制经验半变异函数图拟合理论半变异函数图计算克里格系数进行预测图10.26 克里格方法流程图例如:某地区降水量是一个区域化变量，其变异函数 的实测值及距离h的关系见下表，下面我们试用回归分析方法建立其球状变异函数模型。)(h实测值(h)距离h实测值(h)距离h2.10.69.24.94.31.110.35.15.72.210.56.26.52.510.97.57.83.111.29.58.83.812

8、.49.8 从上面的介绍和讨论，我们知道，球状变异函数的一般形式为 当 时，有ahccahahahcchh03300)223(00)(ah 0330)2()23()(hachacch 如果记 ，则可以得到线性模型 根据表中的数据，对上式进行最小二乘拟合它是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。，得到 (4.2.20) 32132100,21,23,),(hxhxacbacbcbhy22110 xbxbby2192007. 0731. 1048. 2xxy(4.2.19) 比较(4.2.20)式与(4.2.19)式，并做简单计

9、算可知：c0=2.048，c=1.154，a=8.353，所以，球状变异函数模型为535. 8202. 3535. 80)535. 821535. 823(154. 1048. 200)(33*hhhhhh(4.2.21)4个观测点x1，x2，x3，x4的观测值分别为Z(x1)=37、Z(x2)=42、Z(x3)=36、Z(x4)=35，如果假设降水量的变异函数是向同性（即变异函数在各个方向的变化都相同）的二维球状模型，其具体形式为（4.2.21）式。现在，我们用普通克立格法估计观测点x0的降水量值Z(x0)。 当 时， 根据克立格矩阵的对称性，当 时， ，由此计算可得 ji 202. 315

10、4. 1048. 2)0(044332211cccccccji jijiijxxxxcc202. 3)(870. 0)535. 8)2(21535. 8223(154. 1048. 2202. 3)2(202. 3)11(202. 33322042112ccc542. 0)13(202. 3223113cc711. 0)12(202. 322024114ccc601. 0)22(202. 3223223cc383. 0)14(202. 3224334cc466. 0)23(202. 3224224cc 将以上计算结果代入克立格方程组（4.2.31 ），得 952. 0)1(202. 3201c571. 0)3(202. 3203c473. 0301. 0202. 0210. 0287. 01870. 0571. 0711. 0952. 0011111202. 3383. 0466. 0711. 01383. 0202. 3601. 0542. 01466. 0601. 0202. 3870. 01711. 0542. 0870. 0202. 314321 即克立格权重系数分别为：1=0.287，2=0.210,3=0.202