克里金(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^{n}_{i=1}{\lambda_iz_i}- z_o)=0\]

又因为对任意的z都有\(E[z] = c\)，则

\[c \sum^{n}_{i=1}{\lambda_i}- c=0\]

即

\[\sum^{n}_{i=1}{\lambda_i} = 1\]

这是\(\lambda_i\)的约束条件之一。

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

再分析估计误差\(Var(\hat{z_o}-z_o)\)。为方便公式推理，用符号\(J\)表示，即

\[J = Var(\hat{z_o}-z_o)\]

则有

\[\begin{array}{r@{\;=\;}l} J &= Var(\sum^{n}_{i=1}{\lambda_iz_i} – z_o) \\&= Var(\sum^{n}_{i=1}{\lambda_iz_i}) – 2 Cov(\sum^{n}_{i=1}{\lambda_iz_i}, z_o) + Cov(z_o, z_o) \\&= \sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_jCov( z_i, z_j)} – 2 \sum^{n}_{i=1}{\lambda_iCov(z_i, z_o)} + Cov(z_o, z_o) \end{array} \]

为简化描述，定义符号 \(C_{ij} = Cov(z_i,z_j) = Cov(R_i,R_j)\)，这里\(R_i = z_i – c\)，即点\((x_i,y_i)\)处的属性值相对于区域平均属性值的偏差。

则有

\[J = \sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_jC_{ij}} – 2 \sum^{n}_{i=1}{\lambda_iC_{io}} + C_{oo} \]

5.代价函数的最优解

再定义半方差函数 \(r_{ij} = \sigma^2 -C_{ij}\)，带入J中，有

\[\begin{array}{r@{\;=\;}l}J & = \sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j(\sigma^2 – r_{ij})} – 2 \sum^{n}_{i=1}{\lambda_i(\sigma^2 – r_{io})} + \sigma^2 – r_{oo} \\ &=\sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j(\sigma^2)} -\sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j( r_{ij})}-2\sum^{n}_{i=1}{\lambda_i(\sigma^2)}+2 \sum^{n}_{i=1}{\lambda_i(r_{io})}+\sigma^2 – r_{oo} \end{array} \]

考虑到\(\sum^{n}_{i=1}{\lambda_i} = 1\)

\[\begin{array}{r@{\;=\;}l}J &= \sigma^2-\sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})}-2 \sigma^2 +2 \sum^{n}_{i=1}{\lambda_i(r_{io})}+ \sigma^2 – r_{oo}\\&=2 \sum^{n}_{i=1}{\lambda_i(r_{io})} -\sum^{n}_{i=1}\sum^{n}_{j=0}{\lambda_i\lambda_j(r_{ij})} – r_{oo} \end{array} \]

我们的目标是寻找使J最小的一组 \(\lambda_i\)，且J是\(\lambda_i\)的函数，因此直接将J对\(\lambda_i\)求偏导数令其为0即可。即

\[\frac{\partial J}{\partial \lambda_i}= 0;i=1,2,\cdots,n\]

但是要注意的是，我们要保证求解出来的最优 \(\lambda_i\) 满足公式\(\sum^{n}_{i=1}{\lambda_i} = 1\)，这是一个带约束条件的最优化问题。使用拉格朗日乘数法求解，求解方法为构造一个新的目标函数

\[J + 2\phi(\sum^{n}_{i=1}{\lambda_i}-1)\]

其中\(\phi\)是拉格朗日乘数。求解使这个代价函数最小的参数集\({\phi,\lambda_1,\lambda_2,\cdots,\lambda_n}\)，则能满足其在\(\sum^{n}_{i=1}{\lambda_i} = 1\)约束下最小化\(J\)。即

\[\left\{\begin{array}{r@{\;=\;}l}\frac{\partial(J + 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \lambda_k} &= 0;k=1,2,\cdots,n\\ \frac{\partial(J + 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \phi} &= 0 \end{array} \right.\]

\[\left\{\begin{array}{r@{\;=\;}l} \frac{\partial (2 \sum^{n}_{i=1}{\lambda_i(r_{io})} – \sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})} – r_{oo}+ 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \lambda_k} & = 0;k=1,2,\cdots,n\\ \frac{\partial ( 2 \sum^{n}_{i=1}{\lambda_i(r_{io})} – \sum^{n}_{i=1}\sum^{n}_{j}{\lambda_i\lambda_j(r_{ij})} – r_{oo}+ 2\phi(\sum^{n}_{i=1}{\lambda_i}-1))}{\partial \phi} &= 0 \end{array} \right.\]

\[\left\{\begin{array}{r@{\;=\;}l} 2r_{ko} – \sum^{n}_{j=1}{(r_{kj}+r_{jk})\lambda_j}+2\phi&=0;k=1,2,\cdots,n\\ \sum^{n}_{i=1}{\lambda_i} &= 1 \end{array} \right.\]

由于\(C_{ij}=Cov(z_i,z_j)=C_{ji}\)，因此同样地\(r_{ij}=r_{ji}\)，那么有

\[\left\{\begin{array}{r@{\;=\;}l} r_{ko} – \sum^{n}_{j=1}{r_{kj}\lambda_j}+\phi&= 0;k=1,2,\cdots,n\\ \sum^{n}_{i=1}{\lambda_i} &= 1 \end{array} \right.\]

式子中半方差函数\(r_{ij}\)十分重要，最后会详细解释其计算与定义

在以上计算中我们得到了对于求解权重系数\(\lambda_j\)的方程组。写成线性方程组的形式就是：

\begin{equation}\left\{\begin{array}{r@{\;=\;}l} r_{11}\lambda_1 + r_{12}\lambda_2 + \cdots + r_{1n}\lambda_n – \phi&= r_{1o}\\r_{21}\lambda_1 + r_{22}\lambda_2 + \cdots + r_{2n}\lambda_n – \phi&= r_{2o}\\&\cdots\\ r_{n1}\lambda_1 + r_{n2}\lambda_2 + \cdots + r_{nn}\lambda_n – \phi&= r_{no}\\ \lambda_1 + \lambda_2 + \cdots + \lambda_n &= 1\\ \end{array} \right.\end{equation}

写成矩阵形式即为

\[\begin{bmatrix}r_{11}&r_{12}&\cdots&r_{1n}&1\\ r_{21}&r_{22}&\cdots&r_{2n}&1\\\cdots&\cdots&\cdots&\cdots&\cdots\\r_{n1}&r_{n2}&\cdots&r_{nn}&1\\1&1&\cdots&1&0\end{bmatrix}\begin{bmatrix} \lambda_1\\ \lambda_2\\\cdots\\\lambda_n\\-\phi\end{bmatrix}=\begin{bmatrix} r_{1o}\\ r_{2o}\\\cdots\\r_{no}\\1\end{bmatrix}\]

对矩阵求逆即可求解。

唯一未知的就是上文中定义的半方差函数\(r_{ij}\)，接下来将详细讨论

6.半方差函数

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

\[r_{ij} = \sigma^2 -C_{ij}\]

其等价形式为

\[r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]\]

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

根据上文定义 \(R_i = z_i – c\)，有\(z_i-z_j = R_i – R_j\)，则

\[\begin{array}{r@{\;=\;}l} r_{ij} &= \frac{1}{2}E[(R_i-R_j)^2]\\&= \frac{1}{2}E[R_i^2-2R_iR_j+R_j^2]\\&= \frac{1}{2}E[R_i^2]+\frac{1}{2}E[R_j^2]-E[R_iR_j] \end{array} \]

又因为：

\[E[R_i^2] =E[R_j^2] = E[(z_i – c)^2] = Var(z_i) = \sigma^2 \]

\[E[R_iR_j] = E[(z_i – c)(z_j-c)] = Cov(z_i,z_j) = C_{ij}\]

于是有

\[\begin{array}{r@{\;=\;}l} r_{ij} &= \frac{1}{2}E[(z_i-z_j)^2]\\&= \frac{1}{2}E[R_i^2]+\frac{1}{2}E[R_j^2]-E[R_iR_j]\\&= \frac{1}{2}\sigma^2+\frac{1}{2}\sigma^2- C_{ij}\\&=\sigma^2 -C_{ij}\end{array}\]

\( \sigma^2 -C_{ij} = \frac{1}{2}E[(z_i-z_j)^2]\)得证，现在的问题就是如何计算

\[r_{ij} = \frac{1}{2}E[(z_i-z_j)^2]\]

这时需要用到地理学第一定律，空间上相近的属性相近。\(r_{ij} = \frac{1}{2}(z_i-z_j)^2\)表达了属性的相似度；空间的相似度就用距离来表达，定义i与j之间的几何距离

\[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}\]

克里金插值假设\(r_{ij}\)与\(d_{ij}\)存在着函数关系，这种函数关系可以是线性、二次函数、指数、对数关系。为了确认这种关系，我们需要首先对观测数据集

\[\{z(x_1,y_1),z(x_2,y_2),z(x_3,y_3),\cdots,z(x_{n-1},y_{n-1}),z(x_n,y_n)\}\]

计算任意两个点的 距离\(d_{ij}= \sqrt{(x_i-x_j)^2 + (y_i – y_j)^2}\)和 半方差 \(\sigma^2 -C_{ij} =\frac{1}{2}E[(z_i-z_j)^2]\)，这时会得到\(n^2\)个\((d_{ij}, r_{ij})\)的数据对。

将所有的\(d\)和\(r\)绘制成散点图，寻找一个最优的拟合曲线拟合\(d\)与\(r\)的关系，得到函数关系式

\[r = r(d)\]

那么对于任意两点\((x_i,y_i), (x_j,y_j)\)，先计算其距离\(d_{ij}\)，然后根据得到的函数关系就可以得到这两点的半方差\(r_{ij}\)

7. 简单克里金（simple kriging）与普通克里金（ordinary kriging）的区别

以上介绍的均为普通克里金（ordinary kriging）的公式与推理。

事实上普通克里金插值还有简化版，即简单克里金（simple kriging）插值。二者的差异就在于如何定义插值形式：

上文讲到，普通克里金插值形式为

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

而简单克里金的形式则为

\[\hat{z_o} – c= \sum^{n}_{i=1}{\lambda_i(z_i-c)}\]

这里的符号\(c\)在上文介绍过了，是属性值的数学期望，即\(E[z] = c\)。也就是说，在普通克里金插值中，认为未知点的属性值是已知点的属性值的加权求和；而在简单克里金插值中，假设未知点的属性值相对于平均值的偏差是已知点的属性值相对于平均值的偏差的加权求和，用公式表达即为：

\[\hat{R_o} = \sum^{n}_{i=1}{\lambda_iR_i}\]

这里的\(R_i\)在上文定义过了：\(R_i = z_i – c\)。

但是为什么这样的克里金插值称为“简单克里金”呢？由于有假设\(E[z] = c\)，也就是说\(E(R_i + c) = c\)，即\(E(R_i) = 0\)。那么上面的公式\(\hat{R_o} = \sum^{n}_{i=1}{\lambda_iR_i}\)两边的期望一定相同，那么在求解未知参数\(\lambda_i\)就不需要有无偏约束条件\(\sum^{n}_{i=1}{\lambda_i} = 1\)。换句话说，这样的估计公式天生就能满足无偏条件。因此它被称为简单克里金。

从在上文（第4节优化目标/代价函数J）中可以知道，优化目标的推理和求解过程是通过对属性值相对于期望的偏差量\(R_i\)进行数学计算而进行的。也就是说这两种克里金插值方法虽然插值形式不一样，求解方法是一样的，重要的区别是简单克里金插值不需要约束条件\(\sum^{n}_{i=1}{\lambda_i} = 1\)，求解方程组为：

\begin{equation}\left\{\begin{array}{r@{\;=\;}l} r_{11}\lambda_1 + r_{12}\lambda_2 + \cdots + r_{1n}\lambda_n + \phi&= r_{1o}\\r_{21}\lambda_1 + r_{22}\lambda_2 + \cdots + r_{2n}\lambda_n + \phi&= r_{2o}\\&\cdots\\ r_{n1}\lambda_1 + r_{n2}\lambda_2 + \cdots + r_{nn}\lambda_n + \phi&= r_{no}\\ \end{array} \right.\end{equation}

还有更重要的一点，简单克里金的插值公式为：

\[\hat{z_o} = \sum^{n}_{i=1}{\lambda_i(z_i-c)}+c\]

换句话说，在计算未知点属性值\(\hat{z_o}\)前，需要知道该地区的属性值期望\(c\)。事实上我们在进行插值前很难知道这个地区的真实属性值期望。有些研究者可能会采用对观测数据简单求平均的方法计算期望值\(c\)，而考虑到空间采样点位置代表性可能有偏差（比如采样点聚集在某一小片地区，没有代表性），简单平均估计的期望也可能是有偏差的。这是简单克里金方法的局限性。

8.小结

总的来说，进行克里金插值分为这几个步骤：

  1. 对于观测数据，两两计算距离与半方差
  2. 寻找一个拟合曲线（或者其他模型）模拟距离与半方差的关系，从而能根据任意距离计算出相应的半方差
  3. 计算出所有已知点之间的半方差\(r_{ij}\)，直接使用公式\(r_{ij} = \frac{1}{2}(z_i-z_j)^2\)
  4. 对于未知点\(z_o\)，计算它到所有已知点\(z_i\)的距离\(d_{io}\)，并通过第2步的拟合曲线，估计半方差\(r_{io}\)
  5. 求解第四节中的方程组，得到最优系数\(\lambda_i\)
  6. 使用最优系数对已知点的属性值进行加权求和，得到未知点\(z_o\)的估计值，即为\(\hat{z_o} = \sum^{n}_{i=1}{\lambda_iz_i}\)

相关文章:

  1. 加权线性回归(Weighted Linear Regression)的公式及其推理
  2. 用WRF模型进行气象模拟入门(1)——简介以及代码编译
  3. 用WRF模型进行气象模拟入门(2)——WPS的配置与使用
Tags:

96 thoughts on “克里金(Kriging)插值的原理与公式推导”

    1. xg1990

      球面模型是拟合 d 与 r 的关系 时所用到的拟合函数（见6.半方差函数）
      它的参数应该是根据样本点估计出来的： 两两计算样本点的距离 d 跟半方差 r，N 个点得到 N(N-1)/2 组 (r,d)对，然后用各种模型模拟 (r,d)关系

      Reply
    3. Teng Li

      GIS用的Cross-Validation就是DL中经典的留一检验（Leave-one-out），比如你用100个数据，每次把n号数据剔除掉用剩下的99个预测第n个然后计算残差；问题是其实我预测一次就做完了，但验证的计算量确实预测的100倍！数据量越大，越不划算，我猜想可能内部有什么优化算法吧。

      Reply
      1. xg1990

        cross validation方法很多的，leave-one-out只是效率比较低的一种。并不一定需要什么优化算法，换个cross validation方法就行了

        Reply
      2. 王振宇

        您好，在您的评论中学到很多东西，想请教一下您了解泊松克里金方法吗？能否留个联系方剂沟通交流一下？我的邮箱是w97rian@163.com

        Reply
    1. xg1990

      负值说明了空间属性存在负相关。举个具体例子：水面波动，在相距半个波长上的两个点相位永远相反。

      Reply
  5. TOM

    楼主，这篇博客的word版有吗？公式输起来太麻烦，有的话发一份呗！email:yadongzhangcau@sina.com

    Reply
  8. jimbo

    计算观测数据集任意两点的距离和半方差来进行r与d的函数关系拟合，想求教楼主两个点之间的半方差怎么计算

    Reply
      2. erica00

        大神，你在第六节写的用0.5 * (zi−zj)^2来计算半方差，这个和用rij=0.5*E[(zi−zj)^2]计算半方差效果一样吗？

        Reply
        1. xg1990

          我们假设半方差\(r_{ij}\)是距离\(d_{ij}\)的函数，函数值的估计表示为数学期望的形式。
          但是具体的函数形状需要一个一个的观测点来计算，因此需要计算0.5*[(zi−zj)^2]作为观测点来拟合半方差函数

        3. 刘涛

          如果有多组同一空间的数据，那么就可以计算rij = 0.5 * E [（zi – zj）^ 2] 了吧？

  9. tyyyy

    请问每种克里金方法的适用条件是
    什么，如果我的数据转换后仍不符合
    正态分布该怎么办？谢谢

    Reply
    1. xg1990

      不好意思，克里金方法变种太多，我还没能一一学习。普通克里金的假设是空间的期望与方差均一分布，似乎不需要正态分布的假设。

      Reply
      2. Teng Li

        上面的推导过程的确看不出需要硬塞进这个假设，ArcGIS文档也没说清楚这个假设从哪里来的……http://desktop.arcgis.com/en/arcmap/latest/extensions/geostatistical-analyst/examine-the-distribution-of-your-data.htm
        我个人感觉有点像中心极限定理（？），虽然每个采样点自己的分布可以是任意的。但如果采样点足够多，总体的分布就呈现正态（只是这里没有求和或求平均）

        Reply
        1. Teng Li

          对不起，我理解错了，和中心极限没有关系。http://desktop.arcgis.com/zh-cn/arcmap/latest/extensions/geostatistical-analyst/understanding-transformations-and-trends.htm
          “作为预测方法的克里金法并不要求数据具有正态分布。但是，正如了解不同的克里金模型中所述，要获取普通克里金法、简单克里金法和泛克里金法的分位数和概率图则要求数据必须处于正态分布。如果只考虑加权平均值预测方法，则无论数据是否正态分布，克里金法都是最好的无偏差预测方法。但是，如果数据处于正态分布，克里金法会是所有无偏差预测方法（不仅是加权平均值预测）中最好的。”

    2. Teng Li

      普通克里金（ordinary）对数据须满足三条最基本的假设：（1）独立变量正态分布；（2）空间格局全局一致；（3）数据无趋势；此外在解方程上还有无偏假设。
      当这些假设被打破的时候就衍生出一系列变种（除非你手动对数据做变换直接插值后再逆变换回去）：
      （1）如果不满足正态分布，考虑log，arcsin变换或ArcGIS自带的box-cox等，一般一次变换QQ图已经比较好了，如果仍不符合可考虑先看histogram后，剔除异常值与连续多次变换迭代进行。
      （2）空间规律局部变化（即半方差图各地相差很大），首先考虑经验贝叶斯（Empirical Bayesian）+Local Radius参数调节；或者加入协变量（CoKriging），或者手动设置断线（一般是河流、悬崖、道路等）。
      （3）首先观察Explore Data中的Trend Analysis侧面的拟合曲线。数据趋势的剔除，自己做的话可能需要首先拟合一阶或二阶的趋势面，ArcGIS提供的是泛克里金（Universal）工具自动去二次趋势。
      （4）如果除了空间（自）相关外，还利用多个变量相关性提高预测精度，此时变量独立假设打破，进入协同克里金（Ordinary/Universal CoKriging），代价就是额外的交叉相关计算和更多的参数估计。
      （5）特殊数据：前面都在说连续变量（continuous）,如果需要估计的是二值binary data，也无所谓正态分布了；此时，如果用指示克里金预测值（indicator）仍是0/1，如果用概率克里金（Probability kriging），则预测出的是连续概率。
      （6）博主所述，Ordinary和Simple的区别在于无偏性对方程构造的影响，但仍是在线性方程组的领域。一旦进入非线性估计，就是析取（disjunctive）克里金，还要满足个二元正态分布（这个我同样不懂）……
      对于初学者，应尽可能使自己的模型简单化，除非有充分的理由，一般还是推荐使用普通克里金的默认配置。

      Reply
      1. 魏强

        请问，第一个假设：独立变量正态分布，指的是测量点间是独立地吗？如果是的话，测量点相互独立，协方差不就都为0了吗？

        Reply
      2. KyleS

        普通克里金、泛克里金和协同克里金给出最优无偏估计（BLUP）的要求均是随机场满足固有平稳假设。如果使用变异函数建模也要求各向同性假设，但对样本的概率分布没有要求。

        如果样本服从正态分布则克里金法的效果通常更好。如果随机场本身服从正态分布，则BLUP可以基于分位函数给出置信区间。

        一些基于贝叶斯方法的克里金模型要求正态分布的先验假设。

        Reply
    1. xg1990

      从解方程的角度讲，3个点就够了。
      但实际上需要多少点没有定论，取决于具体问题。总的原则是越多越好，分布约随机越好

      Reply
  14. Frank Han

    你好！看了你的文章之后受益匪浅。但我有几个问题：
    请问第5节“代价函数的最优解”中，引入拉格朗日乘数之后的第二个大括号中的第一个式子，分子是不是少了一个“+\phi”？同样的问题也出现在接下来的几个式子中。
    第5节的最后一个矩阵方程，中间那个向量的最后一项是不是应该为“\phi”？

    Reply
  15. 林晨曦

    根据楼主的详细推导，克里金插值利用半方差函数实际只是为了满足误差最小，而并没有体现插值目标点与已知样本点之间的空间自相关关系。比如IDW，在空间上相距较远的点则权值较小，但是kriging似乎并不会因为空间自相关较弱而权值较小，因此插值结果只具有数学意义而不具有地理学的意义。不知道楼主对这个问题的看法是什么？最后十分感谢楼主的文档，受益匪浅！

    Reply
    1. xg1990

      半方差函数泛化了IDW所利用的空间自相关。IDW的「相距较远的点则权值较小」只是一种特殊情况，事实上还可能距离越远权值越大（相关性越高）。更一般地讲，距离与权值大小之间可能是任何函数关系，而半方差函数就是在描述这种关系

      Reply
  16. whu_fish

    错误一：
    第5节，代价函数的最优解中，引入拉格朗日乘子后分别求偏导（即该节中第一个花括号），计算后得到第二个花括号结果错误，第二个花括号中第一个式子左边少了一项“减两倍乘子”，自此下面的公式推导皆少了一项，最终的矩阵形式中，系数矩阵（即第二个矩阵）中最后一行不是0，应是乘子。
    错误二：
    通篇的下标数字0与字母O，再加上i的范围标注应从1开始而不是0，所以会有一定的误导。

    Reply
    1. whu_fish

      关于错误一，“减两倍乘子”应改为“加两倍乘子”，最后的矩阵中最后一行应是“减一倍乘子”

      Reply
  17. 烧饼

    数学渣，跪求详细解释：第4标题代价函数J下面公式的第二步到第三步怎么推出来:
    具体的
    1. 第一块是怎么从Var变成了Cov
    2. 第二块Sigma(lambda * Zi)不是一起的吗，lambda是怎么就拿到外面去了QAQ
    求解释..

    Reply
    1. LEOUN

      (1)直接按照方差公式推导就出来了，很简单。
      (2) lamdbda is a constant, so it can be moved out of the sigma.

      Reply
  19. suedez

    醍醐灌顶！谢谢大神，救我狗命。。。
    小白还想弱弱地请教一下，strict stationary， instinct stationary 和 second-order stationary 到底是什么区别和联系啊 TAT

    Reply
  20. madao

    楼主你好，kriging插值和kriging模型有什么异同？半方差函数和变异函数以及相关函数是不是同一个东西？

    Reply
    1. xg1990

      模型是一个更宽泛的名词了，具体指代什么看上下文吧。半方差函数和变异函数以及相关函数不是同一个东西……

      Reply
  21. 张易

    楼主，已知点的r就直接用rij=1/2（zi-zj）^2计算？在确定rij和dij的函数关系的时候，就是把对应点（（d11,r11），（d12,r12）…）直接把所有点画出来进行观察？不需要处理?（例如删点啥的？）谢谢楼主大神分享，还请赐教。

    Reply
  23. da jax

    作者您好，我有一个问题没想明白，rij=1/2*(zi-zj)^2
    但是zi是一个点，有两个坐标（xi,yi）
    如何计算rij？这个问题还是没有想明白

    Reply
  28. wang

    感谢楼主分享，一般讲解中公式推导时，-2μ 是拉格朗日乘数 ，而不是μ ，能请教下为什么是-2μ ,而不是一般的+μ?

    Reply
  29. Lin

    专业的博主！近期又涉及到GIS方面的东西，蒙圈，请问泛克里金跟普通克里金有啥区别，泛克里金是普通克里金和其他回归方法的结合？本质啥区别？

    Reply
  31. 李酱

    请问一下第3节中：
    zi不是某个确定的值吗？那E[zi]怎么是c而不是zi呢？
    同理预测值也是一个确定值，它的期望不是它本身吗？
    真实值z0是指一个确定的值还是一个在（x0,y0）处的随机变量？

    Reply
  32. zoe zhu

    最后算法步骤中，1计算了两两半方差，怎么3里面还要计算半方差
    没看明白半方差怎么计算的？（dij,rij)中距离d已知公式，但r没有公式，依旧是一个期望，那岂不是应该用样本去预测这个期望？E((zi-zj)的平方)应该怎么求呢？其含义是什么？

    Reply
    1. xg1990

      理论上可以，但是需要好好考虑下怎么计算距离，以及半方差函数的形状。这些都是结合具体数据来分析的

      Reply
      1. li

        可以通过多维变量中的每一维的变化分别计算出距离，然后做出变异函数图，这样做出多个变异函数曲线后，再根据变异函数的嵌套拟合来得出最后的变异函数。不知到这么做可以吗？我现在在处理20维的数据，发现用Matlab的DACE工具包不太好用，不知道您有手写的源码吗，感谢，

        Reply
      2. 7777

        为什么三维或者更高维不能用这个半方差函数，后面拟合的时候只是距离与特征值的拟合，与几维坐标没有关系啊

        Reply
  38. Kevin

    第四、第五部分中的二重求和内层的变量j是否也应从1开始？毕竟Var(z^) = Cov(z^,z^)，展开把求和项提出来的话也应该都从1开始吧？

    Reply
  40. 王武森

    楼主您好！您关于克里金插值的介绍很清晰。在阅读的过程中，我感觉有个地方有点小问题，交流一下：4.优化目标/代价函数J部分J那个式子里j应该是从1开始取值，而不是从0开始取值吧，谢谢！

    Reply
  41. 777

    博主你好，在6.半方差函数这一节里，计算任意已知数据对的半方差时，是直接带入观测值得到rij=((zi-zj)^2)/2就行吗？以及为什么n个已知点可以产生n^2个数据对，需要考虑自身和重复吗，若否是不是应该为n*(n-1)/2个？
    以及第5节求解中 最后那个Bx=l中B是对称的不？

    Reply
  42. 白鹿倚青崖

    楼主，你好，我想问下克里金插值的误差怎么解算啊，
    是解算出权重之后再回代求J吗？
    如果是的话这个J需要除样本数吗？
    还有就是这个权重会求出负值是正常的吗？

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.