这篇文章是转载的一个大神的，因为那个大神的知乎回答的公式坏了，因此整理了一下公式，分享一下，讲的真的挺好的，大神的博客链接：克里金(Kriging)插值的原理与公式推导 - xg1990
0. 引言——从反距离插值(IDW)说起
空间插值问题，就是在已知空间上若干离散点 \(\left(x_i, y_i\right)\) 的某一属性(如气温，海拔)的观测值 \(z_i=z\left(x_i, y_i\right)\) 的条 件下，估计空间上任意一点 \((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_{i=1}^n \frac{1}{d^\alpha} z_i \]其中 \(\alpha\) 通常取1或者2。
即，用空间上所有已知点的数据加权求和来估计末知点的值，权重取决于距离的倒数（或者倒数的平方）。 那么，距离近的点，权重就大；距离远的点，权重就小。
反距离插值可以有效的基于地理学第一定律估计属性值空间分布，但仍然存在很多问题：
-
\(\alpha\) 的值不确定
-
用倒数函数来描述空间关联程度不够准确
因此更加准确的克里金插值方法被提出来了
1. 克里金插值的定义
相比反距离插值，克里金插值公式更加抽象
\[\hat{z_o}=\sum_{i=1}^n \lambda_i z_i \]这里的 \(\lambda_i\) 是权重系数。它同样是用空间上所有已知点的数据加权求和来估计末知点的值。但权重系数并非距离的倒数，而是能够满足点 \(\left(x_o, y_o\right)\) 处的估计值 \(\hat{z}_o\) 与真实值 \(z_o\) 的差最小的一套最优系数，即
\[\min _{\lambda_i} \operatorname{Var}\left(\hat{z}_o-z_o\right) \]同时满足无偏估计的条件
\[E\left(\hat{z}_o-z_o\right)=0 \]2. 假设条件
不同的克里金插值方法的主要差异就是假设条件不同。本文仅介绍普通克里金插值的假设条件与应用。
普通克里金插值的假设条件为，空间属性 \(z\) 是均一的。对于空间任意一点 \((x, y)\) ，都有同样的期望 \(c\) 与方差 \(\sigma^2\) 。即对任意点 \((x, y)\) 都有
\[\begin{gathered} E[z(x, y)]=E[z]=c \\ \operatorname{Var}[z(x, y)]=\sigma^2 \end{gathered} \]换一种说法：任意一点处的值 \(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)\) 处的偏差，其方差均为常数
\[\operatorname{Var}[R(x, y)]=\sigma^2 \]3. 无偏约束条件
先分析无偏估计条件 \(E\left(\hat{z_o}-z_o\right)=0\) ，将 \(\hat{z_o}=\sum_{i=1}^n \lambda_i z_i\) 带入则有
\[E\left(\sum_{i=1}^n \lambda_i z_i-z_o\right)=0 \]又因为对任意的都有 \(E[z]=c\) ，则
\[c \sum_{i=1}^n \lambda_i-c=0 \]即
\[\sum_{i=1}^n \lambda_i=1 \]这是 \(\lambda_i\) 的约束条件之一。
4. 优化目标/代价函数J
再分析估计误差 \(\operatorname{Var}\left(\hat{z_o}-z_o\right)\) 。为方便公式推理，用符号 \(J\) 表示，即
\[J=\operatorname{Var}\left(\hat{z_o}-z_o\right) \]则有
\[\begin{aligned} J & =\operatorname{Var}\left(\sum_{i=1}^n \lambda_i z_i-z_o\right) \\ & =\operatorname{Var}\left(\sum_{i=1}^n \lambda_i z_i\right)-2 \operatorname{Cov}\left(\sum_{i=1}^n \lambda_i z_i, z_o\right)+\operatorname{Cov}\left(z_o, z_o\right) \\ & =\sum_{i=1}^n \sum_{j=0}^n \lambda_i \lambda_j \operatorname{Cov}\left(z_i, z_j\right)-2 \sum_{i=1}^n \lambda_i \operatorname{Cov}\left(z_i, z_o\right)+\operatorname{Cov}\left(z_o, z_o\right) \end{aligned} \]为简化描述，定义符号 \(C_{i j}=\operatorname{Cov}\left(z_i, z_j\right)=\operatorname{Cov}\left(R_i, R_j\right)\) ，这里 \(R_i=z_i-c\) ， 即点 \(\left(x_i, y_i\right)\) 处的属性值相对于区域平均属性值的偏差。则有
\[J=\sum_{i=1}^n \sum_j^n \lambda_i \lambda_j C_{i j}-2 \sum_{i=1}^n \lambda_i C_{i o}+C_{o o} \]5. 代价函数的最优解
再定义半方差函数 \(r_{i j}=\sigma^2-C_{i j}\) ，带入J中，有
\[\begin{aligned} J & =\sum_{i=1}^n \sum_{j=0}^n \lambda_i \lambda_j\left(\sigma^2-r_{i j}\right)-2 \sum_{i=1}^n \lambda_i\left(\sigma^2-r_{i o}\right)+\sigma^2-r_{o o} \\ & =\sum_{i=1}^n \sum_{j=0}^n \lambda_i \lambda_j\left(\sigma^2\right)-\sum_{i=1}^n \sum_{j=0}^n \lambda_i \lambda_j\left(r_{i j}\right)-2 \sum_{i=1}^n \lambda_i\left(\sigma^2\right)+2 \sum_{i=1}^n \lambda_i\left(r_{i o}\right)+\sigma^2-r_{o o} \end{aligned} \]考虑到 \(\sum_{i=1}^n \lambda_i=1\)
\[\begin{aligned} J & =\sigma^2-\sum_{i=1}^n \sum_j^n \lambda_i \lambda_j\left(r_{i j}\right)-2 \sigma^2+2 \sum_{i=1}^n \lambda_i\left(r_{i o}\right)+\sigma^2-r_{o o} \\ & =2 \sum_{i=1}^n \lambda_i\left(r_{i o}\right)-\sum_{i=1}^n \sum_{j=0}^n \lambda_i \lambda_j\left(r_{i j}\right)-r_{o o} \end{aligned} \]我们的目标是寻找使 最小的一组 \(\lambda_i\) ，且 \(J\)是 \(\lambda_i\) 的函数，因此直接将对 \(\lambda_i\) 求偏导数令其为 0 即可。即
\[\frac{\partial J}{\partial \lambda_i}=0 ; i=1,2, \cdots, n \]但是要注意的是，我们要保证求解出来的最优 \(\lambda_i\) 满足公式 \(\sum_{i=1}^n \lambda_i=1\) ，这是一个带约束条件的最优化问 题。使用拉格朗日乘数法求解，求解方法为构造一个新的目标函数
\[J+2 \phi\left(\sum_{i=1}^n \lambda_i-1\right) \]其中 \(\phi\) 是拉格朗日乘数。求解使这个代价函数最小的参数集 \(\phi, \lambda_1, \lambda_2, \cdots, \lambda_n\) ，则能满足其在 \(\sum_{i=1}^n \lambda_i=1\) 约 束下最小化 \(J\) 。即
\[\begin{gathered} \left\{\begin{array}{l} \frac{\partial\left(J+2 \phi\left(\sum_{i=1}^n \lambda_{i-1}-1\right)\right)}{\partial \lambda_k}=0 ; k=1,2, \cdots, n \\ \frac{\partial\left(J+2 \phi\left(\sum_{i=1}^n \lambda_i-1\right)\right)}{\partial \phi}=0 \end{array}\right. \\ \left\{\begin{array}{l} \frac{\partial\left(2 \sum_{i=1}^n \lambda_i\left(r_{i o}\right)-\sum_{i=1}^n \sum_j^n \lambda_{i j} \lambda_j\left(r_{i j}\right)-r_{\infty}+2 \phi\left(\sum_{i=1}^n \lambda_i-1\right)\right)}{\partial \lambda_k}=0 ; k=1,2, \cdots, n \\ \frac{\partial\left(2 \sum_{i=1}^n \lambda_i\left(r_{i o}\right)-\sum_{i=1}^n \sum_j^n \lambda_i \lambda_j\left(r_{i j}\right)-r_{\infty}+2 \phi\left(\sum_{i=1}^n \lambda_i-1\right)\right)}{\partial \phi}=0 \end{array}\right. \\ \left\{\begin{array}{c} 2 r_{k o}-\sum_{j=1}^n\left(r_{k j}+r_{j k}\right) \lambda_j+2 \phi=0 ; k=1,2, \cdots, n \\ \sum_{i=1}^n \lambda_i=1 \end{array}\right. \end{gathered} \]由于 \(C_{i j}=\operatorname{Cov}\left(z_i, z_j\right)=C_{j i}\) ，因此同样地 \(r_{i j}=r_{j i}\) ，那么有
\[\left\{\begin{aligned} r_{k o}-\sum_{j=1}^n r_{k j} \lambda_j+\phi & =0 ; k=1,2, \cdots, n \\ \sum_{i=1}^n \lambda_i & =1 \end{aligned}\right. \]式子中半方差函数 \(r_{i j}\) 十分重要，最后会详细解释其计算与定义
在以上计算中我们得到了对于求解权重系数 \(\lambda_j\) 的方程组。写成线性方程组的形式就是：
写成矩阵形式即为
\[\left[\begin{array}{ccccc} r_{11} & r_{12} & \cdots & r_{1 n} & 1 \\ r_{21} & r_{22} & \cdots & r_{2 n} & 1 \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ r_{n 1} & r_{n 2} & \cdots & r_{n n} & 1 \\ 1 & 1 & \cdots & 1 & 0 \end{array}\right]\left[\begin{array}{c} \lambda_1 \\ \lambda_2 \\ \cdots \\ \lambda_n \\ -\phi \end{array}\right]=\left[\begin{array}{c} r_{1 o} \\ r_{2 o} \\ \cdots \\ r_{n o} \\ 1 \end{array}\right] \]对矩阵求逆即可求解。
唯一末知的就是上文中定义的半方差函数 \(r_{i j}\) ，接下来将详细讨论上
6. 半方差函数
文中对半方差函数的定义为
\[r_{i j}=\sigma^2-C_{i j} \]其等价形式为
\[r_{i j}=\frac{1}{2} E\left[\left(z_i-z_j\right)^2\right] \]这也是半方差函数名称的来由，接下来证明这二者是等价的：
根据上文定义 \(R_i=z_i-c\) ，有 \(z_i-z_j=R_i-R_j\) ，则
又因为：
\[\begin{aligned} E\left[R_i^2\right] & =E\left[R_j^2\right]=E\left[\left(z_i-c\right)^2\right]=\operatorname{Var}\left(z_i\right)=\sigma^2 \\ E\left[R_i R_j\right] & =E\left[\left(z_i-c\right)\left(z_j-c\right)\right]=\operatorname{Cov}\left(z_i, z_j\right)=C_{i j} \end{aligned} \]于是有
\[\begin{aligned} r_{i j} & =\frac{1}{2} E\left[\left(z_i-z_j\right)^2\right] \\ & =\frac{1}{2} E\left[R_i^2\right]+\frac{1}{2} E\left[R_j^2\right]-E\left[R_i R_j\right] \\ & =\frac{1}{2} \sigma^2+\frac{1}{2} \sigma^2-C_{i j} \\ & =\sigma^2-C_{i j} \end{aligned} \]\(\sigma^2-C_{i j}=\frac{1}{2} E\left[\left(z_i-z_j\right)^2\right]\) 得证，现在的问题就是如何计算
\[r_{i j}=\frac{1}{2} E\left[\left(z_i-z_j\right)^2\right] \]这时需要用到地理学第一定律，空间上相近的属性相近。 \(r_{i j}=\frac{1}{2}\left(z_i-z_j\right)^2\) 表达了属性的相似度；空间的相似度就用距离来表达，定义 \(\mathrm{i} \mathrm{j}\) 之间的几何距离
\[d_{i j}=d\left(z_i, z_j\right)=d\left(\left(x_i, y_i\right),\left(x_j, y_j\right)\right)=\sqrt{\left(x_i-x_j\right)^2+\left(y_i-y_j\right)^2} \]克里金插值假设 \(r_{i j}\) 与 \(d_{i j}\) 存在着函数关系，这种函数关系可以是线性、二次函数、指数、对数关系。为了确认 这种关系，我们需要首先对观测数据集
\[\left\{z\left(x_1, y_1\right), z\left(x_2, y_2\right), z\left(x_3, y_3\right), \cdots, z\left(x_{n-1}, y_{n-1}\right), z\left(x_n, y_n\right)\right\} \]计算任意两个点的距离 \(d_{i j}=\sqrt{\left(x_i-x_j\right)^2+\left(y_i-y_j\right)^2}\) 和半方差 \(\sigma^2-C_{i j}=\frac{1}{2} E\left[\left(z_i-z_j\right)^2\right]\) ，这时会得到 \(n^2\) 个 \(\left(d_{i j}, r_{i j}\right)\) 的数据对。
将所有的 \(d\) 和 \(r\) 绘制成散点图，寻找一个最优的拟合曲线拟合 \(d\) 与 \(r\) 的关系，得到函数关系式
\[r=r(d) \]那么对于任意两点 \(i\) 和 \(j\) ，先计算其距离 \(d_{i j}\) ，然后根据得到的函数关系就可以得到这两点的半方差 \(r_{i j}\)
7. 简单克里金 (simple kriging) 与普通克里金 (ordinary kriging) 的区别
以上介绍的均为普通克里金 (ordinary kriging) 的公式与推理。
事实上普通克里金插值还有简化版，即简单克里金 (simple kriging) 插值。二者的差异就在于如何定义插值 形式：
上文讲到，普通克里金插值形式为
而简单克里金的形式则为
\[\hat{z_o}-c=\sum_{i=1}^n \lambda_i\left(z_i-c\right) \]这里的符号c在上文介绍过了，是属性值的数学期望，即 \(E[z]=c_{\mathrm{o}}\) 也就是说，在普通克里金插值中，认为末知点的属性值是已知点的属性值的加权求和；而在简单克里金插值中，假设末知点的属性值相对于平均值的偏差是已知点的属性值相对于平均值的偏差的加权求和，用公式表达即为：
\[\hat{R}_o=\sum_{i=1}^n \lambda_i R_i \]这里的 \(R_i\) 在上文定义过了: \(R_i=z_i-c_{\text {。 }}\)
但是为什么这样的克里金揷值称为"简单克里金“呢? 由于有假设 \(E[z]=c\) ，也就是说 \(E\left(R_i+c\right)=c\) ，即 \(E\left(R_i\right)=0\) 。那么上面的公式 \(\hat{R}_o=\sum_{i=1}^n \lambda_i R_i\) 两边的期望一定相同，那么在求解末知参数 \(\lambda_i\) 就不需要有无 偏约束条件 \(\sum_{i=1}^n \lambda_i=1\) 。换句话说，这样的估计公式天生就能满足无偏条件。因此它被称为简单克里金。
从在上文 (第4节优化目标/代价函数) 中可以知道，优化目标的推理和求解过程是通过对属性值相对于期望 的偏差量 \(R_i\) 进行数学计算而进行的。也就是说这两种克里金揷值方法虽然揷值形式不一样，求解方法是一样 的，重要的区别是简单克里金揷值不需要约束条件 \(\sum_{i=1}^n \lambda_i=1\) ，求解方程组为：
\[\left\{\begin{aligned} r_{11} \lambda_1+r_{12} \lambda_2+\cdots+r_{1 n} \lambda_n+\phi & =r_{1 o} \\ r_{21} \lambda_1+r_{22} \lambda_2+\cdots+r_{2 n} \lambda_n+\phi & =r_{2 o} \\ & \cdots \\ r_{n 1} \lambda_1+r_{n 2} \lambda_2+\cdots+r_{n n} \lambda_n+\phi & =r_{n o} \end{aligned}\right. \]还有更重要的一点，简单克里金的揷值公式为：
\[\hat{z_o}=\sum_{i=1}^n \lambda_i\left(z_i-c\right)+c \]换句话说，在计算末知点属性值 \(\hat{z}_o\) 前，需要知道该地区的属性值期望 \(c_0\) 事实上我们在进行揷值前很难知道这个地区的真实属性值期望。有些研究者可能会采用对观测数据简单求平均的方法计算期望值 \(c\) ，而考虑到空间采样点位置代表性可能有偏差 (比如采样点聚集在某一小片地区，没有代表性)，简单平均估计的期望也可 能是有偏差的。这是简单克里金方法的局限性。
8. 小结
总的来说，进行克里金揷值分为这几个步聚：
- 对于观测数据，两两计算距离与半方差
- 寻找一个拟合曲线 (或者其他模型) 模拟距离与半方差的关系，从而能根据任意距离计算出相应的半方差
- 计算出所有已知点之间的半方差 \(r_{i j}\) ，直接使用公式 \(r_{i j}=\frac{1}{2}\left(z_i-z_j\right)^2\)
- 对于末知点 \(z_o\) ，计算它到所有已知点 \(z_i\) 的距离 \(d_{i o}\) ，并通过第2步的拟合曲线，估计半方差 \(r_{i o}\)
- 求解第四节中的方程组，得到最优系数 \(\lambda_i\)
- 使用最优系数对已知点的属性值进行加权求和，得到末知点 \(z_o\) 的估计值，即为 \(\hat{z_o}=\sum_{i=1}^n \lambda_i z_i\)