matlab克里金插值法,克里金(Kriging)插值的原理与公式推导

学过空间插值的人都知道克里金插值，但是它的变种繁多、公式复杂，还有个半方差函数让人不知所云

本文讲简单介绍基本克里金插值的原理，及其推理过程，全文分为九个部分：

0.引言－从反距离插值说起

1.克里金插值的定义

2.假设条件

3.无偏约束条件

4.优化目标／代价函数

5.代价函数的最优解

6.半方差函数

7.普通克里金与简单克里金

8.小结

0.引言——从反距离插值(IDW)说起

空间插值问题，就是在已知空间上若干离散点 \((x_i,y_i)\) 的某一属性(如气温，海拔)的观测值\(z_i=z(x_i,y_i)\)的条件下，估计空间上任意一点\((x,y)\)的属性值的问题。

直观来讲，根据地理学第一定律，

All attribute values on a geographic surface are related to each other, but closer values are more strongly related than are more distant ones.

大意就是，地理属性有空间相关性，相近的事物会更相似。由此人们发明了反距离插值，对于空间上任意一点\((x,y)\)的属性\(z=z(x,y)\)，定义反距离插值公式估计量

\[\hat{z} = \sum^{n}_{i=1}{\frac{1}{d^\alpha}z_i}\]

其中\(\alpha\)通常取1或者2。

即，用空间上所有已知点的数据加权求和来估计未知点的值，权重取决于距离的倒数(或者倒数的平方)。那么，距离近的点，权重就大；距离远的点，权重就小。

反距离插值可以有效的基于地理学第一定律估计属性值空间分布，但仍然存在很多问题：

\(\alpha\)的值不确定

用倒数函数来描述空间关联程度不够准确

因此更加准确的克里金插值方法被提出来了

1.克里金插值的定义

相比反距离插值，克里金插值公式更加抽象

\[\hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i}\]

其中\(\hat{z_o}\)是点\((x_o,y_o)\)处的估计值，即\(z_o=z(x_o,y_o)\) 。

这里的\(\lambda_i\)是权重系数。它同样是用空间上所有已知点的数据加权求和来估计未知点的值。但权重系数并非距离的倒数，而是能够满足点\((x_o,y_o)\)处的估计值\(\hat{z_o}\)与真实值\(z_o\)的差最小的一套最优系数，即

\[\min_{\lambda_i} Var(\hat{z_o}-z_o)\]

同时满足无偏估计的条件

\[E(\hat{z_o}-z_o)=0\]

2.假设条件

不同的克里金插值方法的主要差异就是假设条件不同。本文仅介绍普通克里金插值的假设条件与应用。

普通克里金插值的假设条件为，空间属性\(z\)是均一的。对于空间任意一点\((x,y)\)，都有同样的期望c与方差\(\sigma^2\)。

即对任意点\((x,y)\)都有

\[E[z(x,y)] = E[z] = c\]

\[Var[z(x,y)] = \sigma^2\]

换一种说法：任意一点处的值\(z(x,y)\)，都由区域平均值\(c\)和该点的随机偏差\(R(x,y)\)组成，即

\[z(x,y)=E[z(x,y)] + R(x,y)] = c + R(x,y)\]

其中\(R(x,y)\)表示点\((x,y)\)处的偏差，其方差均为常数

\[Var[R(x,y)] = \sigma^2\]

3.无偏约束条件

先分析无偏估计条件\(E(\hat{z_o}-z_o)=0\)，将\(\hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i}\)带入则有

\[E(\sum^{