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

时间 2019-12-10

标签克里 kriging 插值原理公式推导繁體版

原文原文链接

转自：http://xg1990.com/blog/archives/222ide

学过空间插值的人都知道克里金插值，可是它的变种繁多、公式复杂，还有个半方差函数让人不知所云
本文讲简单介绍基本克里金插值的原理，及其推理过程。函数

0.引言——从反距离插值(IDW)提及

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

直观来说，根据地理学第必定律，atom

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

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

z^= \sum i = 0 n 1 d α z i

其中 α 一般取1或者2。orm

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

反距离插值能够有效的基于地理学第必定律估计属性值空间分布，但仍然存在不少问题：ci

α 的值不肯定
用倒数函数来描述空间关联程度不够准确

所以更加准确的克里金插值方法被提出来了数学

1.克里金插值的定义

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

z o^= \sum i = 0 n λ i z i

其中 zo^ 是点 (xo,yo) 处的估计值，即 zo=z(xo,yo) 。

这里的 λi 是权重系数。它一样是用空间上全部已知点的数据加权求和来估计未知点的值。但权重系数并不是距离的倒数，而是可以知足点 (xo,yo) 处的估计值 zo^ 与真实值 zo 的差最小的一套最优系数，即

min λ i V a r (z o^- z o)

同时知足无偏估计的条件

E (z o^- z o) = 0

2.假设条件

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

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

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

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

V a r [z (x, y)] = σ 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) 处的误差，其方差均为常数

V a r [R (x, y)] = σ 2

3.无偏约束条件

先分析无偏估计条件 E(zo^−zo)=0 ，将 zo^=∑ni=0λizi 带入则有

E (\sum i = 0 n λ i z i - z o) = 0

又由于对任意的z都有 E[z]=c ，则

c \sum i = 0 n λ i - c = 0

即

\sum i = 0 n λ i = 1

这是 λi 的约束条件之一。

4.优化目标/代价函数J

再分析估计偏差 Var(zo^−zo) 。为方便公式推理，用符号 J 表示，即

J = V a r (z o^- z o)

则有

J = V a r (\sum n i = 0 λ i z i - z o) = V a r (\sum n i = 0 λ i z i) - 2 C o v (\sum n i = 0 λ i z i, z o) + C o v (z o, z o) = \sum n i = 0 \sum n j = 0 λ i λ j C o v (z i, z j) - 2 \sum n i = 0 λ i C o v (z i, z o) + C o v (z o, z o)

为简化描述，定义符号 Cij=Cov(zi,zj)=Cov(Ri,Rj) ，这里 Ri=zi−c ，即点 (xi,yi) 处的属性值相对于区域平均属性值的误差。

则有

J = \sum i = 0 n \sum j n λ i λ j C i j - 2 \sum i = 0 n λ i C i o + C o o

5.代价函数的最优解

再定义半方差函数 rij=σ2−Cij ，带入J中，有

J = \sum n i = 0 \sum n j = 0 λ i λ j (σ 2 - r i j) - 2 \sum n i = 0 λ i (σ 2 - r i o) + σ 2 - r o o = \sum n i = 0 \sum n j = 0 λ i λ j (σ 2) - \sum n i = 0 \sum n j = 0 λ i λ j (r i j) - 2 \sum n i = 0 λ i (σ 2) + 2 \sum n i = 0 λ i (r i o) + σ 2 - r o o

考虑到 ∑ni=0λi=1

J = σ 2 - \sum n i = 0 \sum n j λ i λ j (r i j) - 2 σ 2 + 2 \sum n i = 0 λ i (r i o) + σ 2 - r o o = 2 \sum n i = 0 λ i (r i o) - \sum n i = 0 \sum n j = 0 λ i λ j (r i j) - r o o

咱们的目标是寻找使J最小的一组 λi ，且J是 λi 的函数，所以直接将J对 λi 求偏导数令其为0便可。即

\partial J \partial λ i = 0; i = 1, 2, \dots, n

可是要注意的是，咱们要保证求解出来的最优 λi 知足公式 ∑ni=0λi=1 ，这是一个带约束条件的最优化问题。使用拉格朗日乘数法求解，求解方法为构造一个新的目标函数

J + ϕ (\sum i = 0 n λ i - 1)

其中 ϕ 是拉格朗日乘数。求解使这个代价函数最小的参数集 ϕ,λ1,λ2,⋯,λn ，则能知足其在 ∑ni=0λi=1 约束下最小化 J 。即

⎧ ⎩ ⎨ ⎪ ⎪ \partial ( J + ϕ ( \sum n i = 0 λ i - 1 ) ) \partial λ i \partial ( J + ϕ ( \sum n i = 0 λ i - 1 ) ) \partial ϕ = 0; i = 1, 2, \dots, n = 0

⎧ ⎩ ⎨ ⎪ ⎪ \partial ( 2 \sum n i = 0 λ i ( r i o ) - \sum n i = 0 \sum n j λ i λ j ( r i j ) - r o o ) \partial λ i \partial ( \partial 2 \sum n i = 0 λ i ( r i o ) - \sum n i = 0 \sum n j λ i λ j ( r i j ) - r o o ) \partial ϕ = 0; i = 1, 2, \dots, n = 0

{2 r i o - \sum n j = 1 (r i j + r j i) λ j \sum n i = 0 λ i = 0; i = 1, 2, \dots, n = 1

因为 Cij=Cov(zi,zj)=Cji ，所以一样地 rij=rji ，那么有

{r i o - \sum n j = 1 r i j λ j \sum n i = 0 λ i = 0; i = 1, 2, \dots, n = 1

式子中半方差函数 rij 十分重要，最后会详细解释其计算与定义

在以上计算中咱们获得了对于求解权重系数 λj 的方程组。写成线性方程组的形式就是：

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ r 11 λ 1 + r 12 λ 2 + \dots + r 1 n λ n + ϕ r 21 λ 1 + r 22 λ 2 + \dots + r 2 n λ n + ϕ r n 1 λ 1 + r n 2 λ 2 + \dots + r n n λ n + ϕ λ 1 + λ 2 + \dots + λ n = r 1 o = r 2 o \dots = r n o = 1 (1)

写成矩阵形式即为

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ r 11 r 21 \dots r n 1 1 r 12 r 22 \dots r n 2 1 \dots \dots \dots \dots \dots r 1 n r 2 n \dots r n n 1 11 \dots 10 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ λ 1 λ 2 \dots λ n 0 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ r 1 o r 2 o \dots r n o 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

对矩阵求逆便可求解。

惟一未知的就是上文中定义的半方差函数 rij ，接下来将详细讨论

6.半方差函数

上文中对半方差函数的定义为

r i j = σ 2 - C i j

其等价形式为

r i j = 1 2 E [(z i - z j) 2]

这也是半方差函数名称的来由，接下来证实这两者是等价的：

根据上文定义 Ri=zi−c ，有 zi−zj=Ri−Rj ，则

r i j = 1 2 E [(R i - R j) 2] = 1 2 E [R 2 i - 2 R i R j + R 2 j] = 1 2 E [R 2 i] + 1 2 E [R 2 j] - E [R i R j]

又由于：

E [R 2 i] = E [R 2 j] = E [(z i - c) 2] = V a r (z i) = σ 2

E [R i R j] = E [(z i - c) (z j - c)] = C o v (z i, z j) = C i j

因而有

r i j = 1 2 E [(z i - z j) 2] = 1 2 E [R 2 i] + 1 2 E [R 2 j] - E [R i R j] = 1 2 σ 2 + 1 2 σ 2 - C i j = σ 2 - C i j

σ2−Cij=12E[(zi−zj)2] 得证，如今的问题就是如何计算

r i j = 1 2 E [(z i - z j) 2]

这时须要用到地理学第必定律，空间上相近的属性相近。 rij=12(zi−zj)2 表达了属性的类似度；空间的类似度就用距离来表达，定义i与j之间的几何距离

d i j = d (z i, z j) = d ((x i, y i), (x j, y j)) = (x i - x j) 2 + (y i - y j) 2 - - - - - - - - - - - - - - - - - - \sqrt

克里金插值假设 rij 与 dij 存在着函数关系，这种函数关系能够是线性、二次函数、指数、对数关系。为了确认这种关系，咱们须要首先对观测数据集

{z (x 1, y 1), z (x 2, y 2), z (x 3, y 3), \dots, z (x n - 1, y n - 1), z (x n, y n)}

计算任意两个点的距离 dij=(xi−xj)2+(yi−yj)2−−−−−−−−−−−−−−−−−−√ 和半方差 σ2−Cij=12E[(zi−zj)2] ，这时会获得 n2 个