克里金插值法

克里金插值法及其适用范围

29 巴任若测绘学院

克里金插值法又称空间局部插值法，是以变异函数理论和结构分析为基础，在有限区域内对区域化变量进行无偏最优估计的一种方法，是地统计学的主要内容之一，由南非矿产工程师D. Matheron于1951年在寻找金矿时首次提出，法国著名统计学家G. Matheron随后将该方法理论化、系统化，并命名为Kriging，即克里金插值法。

1 克里金插值法原理

克里金插值法的适用范围为区域化变量存在空间相关性，即如果变异函数和结构分析的结果表明区域化变量存在空间相关性，则可以利用克里金插值法进行内插或外推。其实质是利用区域化变量的原始数据和变异函数的结构特点，对未知样点进行线性无偏、最优估计，无偏是指偏差的数学期望为0，最优是指估计值与实际值之差的平方和最小[1]。因此，克里金插值法是根据未知样点有限领域内的若干已知样本点数据，在考虑了样本点的形状、大小和空间方位，与未知样点的相互空间关系，以及变异函数提供的结构信息之后，对未知样点进行的一种线性无偏最优估计。

假设研究区域a上研究变量Z（x），在点xi∈A（i=1，2，……，n）处属性值为Z（xi），则待插点x0∈A处的属性值Z（x0）的克里金插值结果Z*（x0）是已知采样点属性值Z（xi）（i=1，2，……，n）的加权和，即：

)()(10*

i n i i x Z x Z ∑==λ

（1）

式中i λ是待定权重系数。

其中Z(xi)之间存在一定的相关关系，这种相关性除与距离有关外，还与其相对方向变化有关，克里金插值方法将研究的对象称“区域化变量”

针对克里金方法无偏、最小方差条件可得到无偏条件可得待定权系数i λ (i=1，2，……，n)满足关系式：

11=∑=n i i λ

（2）

以无偏为前提，kriging 方差为最小可得到求解待定权系数i λ的方程组：

⎪⎪⎩⎪⎪⎨⎧=⋯⋯==+∑∑==1)n ,2,1)(,(),(101n i i j j i n i i j x x C x x C λμλ， （3）

式中，C （xi ，xj ）是Z(xi)和Z(xj)的协方差函数。

2 国内外研究进展

从克里金方法被提出到现在已有完善的理论，并在很多领域得到了实际的应用，在某些领域的应用又推动了克里金理论的发展[3]。它的发展可归纳为四个时期，每个时期都是以每一届地质统计学大会的召开为标志。第一时期，初次提出了地质统计学理论，将地质统计学与传统的统计学分开，且提出了区域化变量、简单克里金、普通克

里金、泛克里金的概念。第二时期，地质统计学的理论逐步的幵始改进和完善。第三时期，地质统计学克里金在实践应用的发展相对理论发展更快，形成了两种类型的理论体系：一类是有参数的克里金方法，另一类是没有参数的克里金方法，有参数的克里金方法是指所研究的数据必须符合正态分布，如析取克里金；而没有参数的克里金方法对所研究的变量的分布没有特殊要求，如指示克里金和概率克里金。第四时期，克里金方法的应用领域不断扩展壮大，在研究中有很多新的课题产生，克里金所研究对象已经不再局限于空间领域的变量，随着某些领域的需求，正在向时间-空间领域扩展[4]。

从目前来看，克里金技术的发展可以概括如下：

(1)形成了一套完整的理论体系。线性平稳地质统计学是地质统计学的基础部分，包含基本概念：区域化变量理论；基本工具：变差函数；基本假设：二阶平稳假设和本征假设；基本公式：估计反差和普通克里金法；线性非平稳地质统计学包括了泛克里金和K阶本征函数法等。平稳非线性地质统计学包含析取克里金等。

(2)编制了一些实际有效的程序以及软件。例如斯坦福大学的Geostatistical Earth Modeling Software。

(3)地质统计学的提出原本是为了解决矿产储量的估计，但是随着地质统计学的发展，人们发现其研究对象存在于很多种自然现象中。于是，地质统计学不再是研究地质领域的特有方法，而成为研究某类自然现象通用的方法，例如降水量的分布、水文层的渗透率和孔隙度等属性值、在医学上对骨豁的三维重建[5]等等。

目前国内外学者利用克里金插值法做了大量研究。翟进乾应用克里金插值方法对煤层分布监测进行了系统分析研究[6]；张蕾、陈晓宏将克里金插值方法用于珠江三角洲网河区水位空间插值[7]；尚庆生、郭建文等将克里金插值方法用于计算青藏铁路钻孔地温数据，实现了数据的体视化[8]；颜辉武，祝国瑞等采用克里金插值方法建立水文地质层三维模型[9]，并利用体绘制技术进行可视化表达，取得了良好的效果；刘承香、阮双深、伍小芹提出基于克里金插值方法进行水深数据插值形成规则网格数字高程模型的算法，对海底数字地图的模拟具有重要参考价值，数字仿真结果证明该算法可行[10]。

3 方法步骤

克里金插值法的应用步骤如下：

1、输入原始数据，即采样点，下面以输入三个采样点求待估插值为例来进行说明。如图1所示：

图1 采样点图示

2、网格化，选择区域的范围和网格的大小，对区域进行网格化处理。

3、数据检验与分析，根据采样值是否合乎实际情况，剔除明显差异点。

4、直方图的计算，直方图有助于掌握区域变化的分布规律，以

便决定是否对原始数据进行转换。

5、利用变异函数进行变异函数计算，了解变量的空间结构。

6、克里金插值估计

（1）待估点权重系数估计

利用多边形估计的方法，首先确定离待估点最近的采样点的权重，根据公式（4）进行采样点权重估计： ∑=

++

=n i w w i i d c d c 111

λ （4）

（2）根据搜索策略选择合适的参估点，如图2：

图2 参估点图示

（3）根据已经求出的变异函数以及采样点数量，三个采样点列出三个等式，求出方程组的系数，公式为：

⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡)3,0()2,0()1,0()3,3()2,3()1,3()3,2()2,2()1,2()3,1()2,1()1,1(321C C C C C C C C C C C C λλλ

（5）