克里金(Kriging)插值

本文为克里金插值的原理及推导过程的学习笔记。

反距离插值(IDW)

空间插值问题，就是在已知空间上若干离散点的某一属性的条件下，估计空间上任意一点的属性值问题。
根据地理学第一定律：
All attribute values on a geographic surface are related to each other, but closer values are more strongly related than are more distant ones.
对于空间上任意一点的属性，定义反距离插值公式:

$$\hat z = \sum _{i=1}^n \frac 1{d ^\alpha} z _i$$

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

1. 的值不确定
2. 倒数函数描述空间关联程度并不准确

为此，有人提出了更为准确的克里金插值方法。

克里金插值定义

克里金插值公式与反距离插值公式相类似：

$$\hat z _0 = \sum _{i=1}^n \lambda _i z _i$$

式中的为权重系数，是能够满足点处的估计值与真实值之差的方差最小的一套最优系数，即

$$\underset {\lambda_i}{min} Var(\hat z _0 - z _0)$$

同时满足无偏估计条件

$$E(\hat z _0 - z _0) = 0$$

假设条件

不同克里金插值方法的主要差异就是假设条件不同。
普通克里金插值的假设条件为：空间属性z是均一的，对于空间任意一点(x,y)，都有相同的期望和方差。
即对任意一点(x,y)都有

$$E[Z(X,Y)] = E[Z] = c$$ $$Var[z(x,y)] = \sigma _2$$

或者说：任意一点(x,y)处的空间属性值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$$

无偏约束条件

要满足无偏估计条件，将带入则有

$$E(\sum _{i=1}^n \lambda _i z _i - z _0) = 0$$

又因为对于任意的z都有，则

$$c \sum _{i=1}^n \lambda _i - c = 0$$

即

$$\sum _{i=1}^n \lambda _i = 1$$

上式为的约束条件之一。

优化目标/代价函数J

令符号J表示估计误差

$$J = Var(\hat z _0 - z _0)$$

则有

$$\begin{aligned} J&= Var(\sum _{i=1}^n \lambda _i z _i - z _0)\\\\ &= Var(\sum _{i=1}^n \lambda _i z _i) + Var(z _0) - 2Cov(\sum _{i=1}^n \lambda _i z _i, z _0)\\\\ &= \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j Cov(z _i,z _j) -2 \sum _{i=1}^n \lambda _i Cov(z _i,z _0) + Var(z _0)\\\\ \end{aligned}$$

为简化符号描述，定义,这里，即点处的属性值相对于区域平均属性值的偏差。

则上式可写为

$$J = \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j C _{ij} - 2 \sum _{i=1}^n \lambda _i C _{i0} + C _{00}$$

代价函数的最优解

定义半方差函数，带入J中，有

$$\begin{aligned} J&= \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j (\sigma^2-r _{ij}) - 2 \sum _{i=1}^n \lambda _i (\sigma^2-r _{i0}) + \sigma^2-r _{00}\\\\ &= \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j \sigma^2 - \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j r _{ij} - 2 \sum _{i=1}^n \lambda _i \sigma^2 + 2 \sum _{i=1}^n \lambda _i r _{i0} + \sigma^2 - r _{00} \end{aligned}$$

考虑到，则

$$\begin{aligned} J&= \sigma^2 - \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j r _{ij} - 2 \sigma^2 + 2 \sum _{i=1}^n \lambda _i r _{i0} + \sigma^2 - r _{00}\\\\ &= 2 \sum _{i=1}^n \lambda _i r _{i0} - \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j r _{ij} - r _{00} \end{aligned}$$

为找到使J最小的一组(J是的函数)，将J对求导并令其为0。即

$$\frac {\partial J}{\partial \lambda _i} = 0; i=1,2,\cdots,n$$

需要注意，最优解还需要满足条件，这是一个带约束条件的最优化问题。
使用拉格朗日乘数法求解，首先构造一个新的目标函数：

$$J + 2\phi(\sum _{i=1}^n \lambda _i -1)$$

其中是拉格朗日乘数。求解使这个代价函数最小的参数集，则能在满足约束下最小化J。即

$$\begin{cases} \frac {\partial(J + 2\phi(\sum _{i=1}^n \lambda _i -1))}{\partial \lambda _k} &=0; k=1,2,\cdots,n\\\\ \frac {\partial(J + 2\phi(\sum _{i=1}^n \lambda _i -1))}{\partial \phi} &=0 \end{cases}$$ $$\begin{cases} \frac {\partial(2 \sum _{i=1}^n \lambda _i r _{i0} - \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j r _{ij} - r _{00} + 2\phi(\sum _{i=1}^n \lambda _i -1))}{\partial \lambda _k} &=0; k=1,2,\cdots,n\\\\ \frac {\partial(2 \sum _{i=1}^n \lambda _i r _{i0} - \sum _{i=1}^n \sum _{j=1}^n \lambda _i \lambda _j r _{ij} - r _{00} + 2\phi(\sum _{i=1}^n \lambda _i -1))}{\partial \phi} &=0 \end{cases}$$

化简得

$$\begin{cases} 2r _{k0} - \sum _{j=1}^n (r _{kj}+r _{jk}) \lambda _j + 2 \phi &= 0; k=1,2,\cdots,n\\\\ \sum _{i=1}^n \lambda _i &= 1 \end{cases}$$

由于，因此,同样地，那么有

$$\begin{cases} r _{k0} - \sum _{j=1}^n r _{kj} \lambda _j + \phi &= 0; k=1,2,\cdots,n\\\\ \sum _{i=1}^n \lambda _i &= 1 \end{cases}$$

将上式写为线性方程组形式：

$$\begin{cases} r _{11} \lambda _1 + r _{12} \lambda _2 + \cdots + r _{1n} \lambda _n - \phi &= r _10\\\\ r _{21} \lambda _1 + r _{22} \lambda _2 + \cdots + r _{2n} \lambda _n - \phi &= r _10\\\\ &\vdots\\\\ r _{n1} \lambda _1 + r _{n2} \lambda _2 + \cdots + r _{nn} \lambda _n - \phi &= r _10\\\\ \lambda _1 + \lambda _2 + cdots + \lambda _n &= 1 \end{cases}$$

写为矩阵形式：

$$\begin{bmatrix} r _{11} & r _{12} & \cdots & r _{1n} & 1\\\\ r _{21} & r _{22} & \cdots & r _{2n} & 1\\\\ \vdots & \vdots & \vdots & \vdots & \vdots\\\\ r _{n1} & r _{n2} & \cdots & r _{nn} & 1\\\\ 1 & 1 & \cdots & 1 & 0 \end{bmatrix} \begin{bmatrix} \lambda _1\\\\ \lambda _2\\\\ \vdots\\\\ \lambda _n\\\\ - \phi \end{bmatrix}= \begin{bmatrix} r _{10}\\\\ r _{20}\\\\ \vdots\\\\ r _{n0}\\\\ 1 \end{bmatrix}$$

求解该线性方程组即可得到所需参数集。

半方差函数

因为：

$$\begin{aligned} \sigma _2 &= Var(z _i) = E(z _i^2) - E^2(z _i) = E(z _i^2) - 2 E^2(z _i) + E^2(z _i)\\\\ &= E(z _i^2) - 2 c E(z _i) + c^2 = E(z _i^2 -2 c z _i + c^2) = E[(z _i - c)^2]\\\\ &= E(R _i^2) \end{aligned}$$ $$\begin{aligned} C _{ij} &= Cov(z _i,z _j) = E(z _i z _j) - E(z _i)E(z _j)\\\\ &= E(z _i z _j) - c E(z _i) - c E(z _j) + c^2 = E(z _i z _j - c z _i - c _j + c^2)\\\\ &= E[(z _i - c)(z _j - c)]\\\\ &= E(R _i R _j) \end{aligned}$$

于是有：

$$\begin{aligned} r _{ij} &= \sigma^2 - C _{ij}\\\\ &= \frac {1}{2} \sigma^2 + \frac {1}{2} \sigma^2 - C _{ij}\\\\ &= \frac {1}{2} E(R _i^2) + \frac {1}{2} E(R _j^2) - E(R _i R _j)\\\\ &= \frac {1}{2} E(R _i^2 - 2 R _i R _j + R _j^2)\\\\ &= \frac {1}{2} E[(R _i^2 - R _j^2)^2]\\\\ &= \frac {1}{2} E[(R _i^2 + c - R _j^2 + c)^2]\\\\ &= \frac {1}{2} E[(z _i^2 - z _j^2)^2] \end{aligned}$$

那么，计算的问题就转换为计算。
根据地理学第一定律，空间上距离相近的点属性相似。表达了两点属性的相似度；空间上的距离定义为：

$$d _{ij} = d(z _i,z _j) = d[(x _i,y _i),(x _j,y _j)] = \sqrt{(x _i - x _j)^2 + (y _i - y _j)^2}$$

克里金插值假设空间距离和属性的相似度之间存在着函数关系。

为了确认这种关系，首先需要对观测数据集计算任意两点之间的距离和半方差/属性相似度。将所有的和绘制成散点图，并寻找一个最优的拟合曲线来拟合d与r的关系。

得到该拟合曲线之后，对于任意两点，只要先计算得两点间的空间距离，便可得到两点的半方差。

一般情况下，如果地理属性之间存在空间自相关，那么 r 会随着 d 的增大而逐渐增大，d 增大到一定程度(range)后 r 会收敛于某个值(sill)，表示距离 d 超过 range 的观测值对预测值的自相关性已经可以忽略不计，只需要使用 range 范围内的观测数据来计算预测值。

补充证明1

UMVUE(一致最小方差无偏估计量、最优无偏估计量)

为找到一个“最好”的无偏估计量。
对无偏估计可使用估计偏差的方差作为判据。
对有偏估计使用均方误差作为判据。

$$\begin{aligned} MSE&= E[(\hat \theta -\theta)^2]\\\\ &= E[\epsilon^2]\\\\ &= Var(\epsilon) + E^2(\epsilon)\\\\ &= Var(\epsilon) + E^2(\hat \theta -\theta) \end{aligned}$$

均方误差MSE由估计偏差的方差与估计偏差的平方两部分组成。
如果满足无偏估计，则，此时可用估计偏差的方差来代替均方误差；若不满足无偏估计，则需要考虑均方误差MSE。

补充证明2

期望、方差基本运算

运算1：方差基本运算

$$\begin{aligned} Var(x)&= E[(x-E(x))^2]\\\\ &= E(x^2) - E^2(x) \end{aligned}$$

运算2：方差线性组合运算1

$$\begin{aligned} Var(x+y)&= E[(x+y-E(x+y))^2]\\\\ &= E[(x-E(x)+y-E(y))^2]\\\\ &= E[(x-E(x))^2] + E[(y-E(y))^2] + 2E[(x-E(x))(y-E(y))]\\\\ &= Var(x) + Var(y) + 2Cov(x,y) \end{aligned}$$

运算3：方差线性组合运算2

$$\begin{aligned} Var(\sum _{i=1}^n x _i)& = \sum _{i,j=1}^n Cov(x _i,x _j)\\\\ &= \sum _{i=1}^n Var(x _i) + \sum _{i \neq j} Cov(x _i,x _j) \end{aligned}$$

运算4：协方差基本运算1

$$\begin{aligned} Cov(x,y)&= E[(x-E(x))(y-E(y))]\\\\ &= E(xy) - E(x)E(y) - E(y)E(x) + E(x)E(y)\\\\ &= E(xy) - E(x)E(y) \end{aligned}$$

运算5：协方差基本运算2

$$\begin{aligned} Cov(x,x)&= E[(x-E(x))(x-E(x))]\\\\ &= E(x^2) - E^2(x)\\\\ &= Var(x) \end{aligned}$$

运算6：协方差基本运算3

$$Cov(\sum _{i=1}^n x _i,\sum _{j=1}^m y _j) = \sum _{i=1}^n \sum _{j=1}^m Cov(x _i,y _j)$$