克里金法

）。方块表示数据位置。红色所示的克里金插值曲线沿灰色所示正态分布可信区间的均值延伸。虚线显示了一个平滑的样条，但它显著偏离了由这些均值给出的期望值。]]
克里金或克里金法（kriging，Kriging，发音为）也称为高斯过程回归，在统计学中（最初在地理统计中），是一种基于受先验协方差支配的高斯过程的插值方法。在适当的先验假设下，克里金法在未采样位置给出最佳线性无偏预测（best linear unbiased prediction，BLUP）。基于其他准则（例如）的插值方法可能无法产生最佳线性无偏预测。该方法广泛应用于空间分析和领域。由于諾伯特·維納和安德烈·柯爾莫哥洛夫的贡献，该技术也被称为維納-柯爾莫哥洛夫预测（Wiener–Kolmogorov prediction）。

该方法的理论基础由法国数学家（Georges Matheron）于1960年建立，其依据是（Danie G. Krige）的硕士论文。克里金是南非威特沃特斯兰德礁复合体距离加权平均金品位的开拓性绘图者。克里金试图根据少数钻孔的样本来估算金最可能的分布。英文动词为to krige，最常用的名词是kriging。在文献中，该词有时首字母大写为Kriging。

虽然其基本公式的计算量很大，但克里金法可以使用各种扩展到更大的问题。

主要原理
相关术语与技术
克里金通过计算给定点邻域内函数已知值的加权平均来预测该点的函数值。该方法与迴歸分析密切相关。两种理论都基于对协方差的假设推导出最佳线性无偏估计量，利用高斯-马尔可夫定理证明估计值与误差的独立性，并使用非常相似的公式。即便如此，它们分别应用于不同的框架中：克里金法用于估计随机场的一次实现，而回归模型则基于多变量数据集的多次观测。

克里金估计也可以被视为中的样条函数，其再生核由协方差函数给出。与经典克里金方法的区别在于解释方式：样条的动机是基于希尔伯特空间结构的最小范数插值，而克里金法的动机是基于随机模型的期望平方预测误差。

具有“多项式趋势面”的克里金法在数学上等同于广义最小二乘法多项式曲線擬合。

克里金法也可以被理解为贝叶斯优化的一种形式。克里金法从函数上的先验分布开始。该先验采用高斯过程的形式：来自函数的N个样本将服从正态分布，其中任意两个样本之间的协方差是在两点空间位置上评估的高斯过程协方差函数（或核）。随后观测到一组值，每个值都与一个空间位置相关联。现在，通过将高斯先验与每个观测值的高斯似然函数相结合，可以在任何新的空间位置预测新值。产生的后验概率分布也是高斯的，其均值和协方差可以简单地从观测值、它们的方差以及从先验导出的核矩阵中计算出来。

地质统计估计量
（DEM）区块（使用matplotlib可视化）]]
在地质统计模型中，采样数据被解释为随机过程的结果。这些模型在概念化中包含不确定性，这一事实并不意味着现象（如森林、含水层、矿床）是由随机过程产生的，而是允许人们为未观测位置的数量空间推论建立方法论基础，并量化与估计量相关的不确定性。

在该模型背景下，随机过程仅仅是处理从样本中收集的数据集的一种方式。地质统计调制的第一步是创建一个最能描述观测数据集的随机过程。

来自位置x_1（地理坐标集的通用名称）的值被解释为随机变量Z(x_1)的一次实现z(x_1)。在样本分散的空间A中，存在相互关联的随机变量Z(x_1), Z(x_2), \ldots, Z(x_N)的N次实现。

这组随机变量构成了随机函数，其中只有一次实现是已知的，即观测数据集z(x_i)。由于每个随机变量只有一次实现，理论上不可能确定单个变量或函数的任何统计参数。地质统计形式主义中提出的解决方案包括“假设”随机函数具有不同程度的“平稳性”，以便使某些统计值的推论成为可能。

例如，如果基于变量分布区域A中样本的同质性，假设一阶矩是平稳的（即所有随机变量具有相同的均值），那么就假设均值可以通过采样值的算术平均值来估计。

与二阶矩相关的平稳性假设定义如下：两个随机变量之间的相关性仅取决于它们之间的空间距离，而与它们的位置无关。因此，如果\mathbf{h} = x_2 - x_1且h = |\mathbf{h}|，则：C\big(Z(x_1), Z(x_2)\big) = C\big(Z(x_i), Z(x_i + \mathbf{h})\big) = C(h)\gamma\big(Z(x_1), Z(x_2)\big) = \gamma\big(Z(x_i), Z(x_i + \mathbf{h})\big) = \gamma(h)为简单起见，我们定义C(x_i, x_j) = C\big(Z(x_i), Z(x_j)\big)和\gamma(x_i, x_j) = \gamma\big(Z(x_i), Z(x_j)\big)。

该假设允许人们推论出这两种度量——和协方差函数：\gamma(h) = \frac{1}{2|N(h)|} \sum_{(i,j)\in N(h)} \big(Z(x_i) - Z(x_j)\big)^2C(h) = \frac{1} \sum_{(i,j)\in N(h)} \big(Z(x_i) - m(h)\big)\big(Z(x_j) - m(h)\big)其中：m(h)=\frac{1}{2|N(h)|} \sum_{(i,j)\in N(h)} \left( Z(x_i) + Z(x_j) \right)N(h)表示观测对i，j的集合，使得|x_i - x_j| = h，而|N(h)|是集合中对的数量。

在该集合中，(i,;j)和(j,;i)表示相同的元素。通常使用“近似距离”h，并结合一定的容差来实现。

线性估计
数量Z \colon \mathbb{R}^n \to \mathbb{R}在未观测位置x_0的空间推论或估计，是根据观测值z_i = Z(x_i)和权重w_i(x_0)；i = 1, \ldots, N的线性组合计算得出的： \hat{Z}(x_0) =
\begin{bmatrix}
w_1 & w_2 & \cdots & w_N
\end{bmatrix}
\begin{bmatrix}
z_1 \\
z_2 \\
\vdots \\
z_N
\end{bmatrix} =
\sum_{i=1}^N w_i(x_0) Z(x_i).权重w_i旨在总结空间推论过程中的两个极其重要的步骤：

反映样本与估计位置x_0的结构“接近度”；
同时，它们应该具有去聚集效应，以避免由最终的样本“聚类”引起的偏差。

在计算权重w_i时，地质统计形式主义有两个目标：“无偏性”和“最小估计方差”。

如果将实际值Z(x_0)的点云与估计值\hat{Z}(x_0)进行对比，全局无偏性、该场的“本征平稳性”或广义平稳性的标准意味着估计值的均值必须等于实际值的均值。

第二个标准指出平方偏差\big(\hat{Z}(x) - Z(x)\big)的均值必须最小，这意味着当估计值与实际值的点云更加分散时，估计量就更加不精确。

方法
根据随机场的随机性质和假设的不同程度的平稳性，可以推导出计算权重的不同方法，即应用不同类型的克里金法。经典方法包括：

普通克里金法（Ordinary kriging）假设仅在x_0的搜索邻域内存在恒定的未知均值。
简单克里金法（Simple kriging）假设一阶矩在整个域内平稳且均值已知：E{Z(x)} = E{Z(x_0)} = m，其中m是已知均值。
（Universal kriging）假设一个通用的多项式趋势模型，例如线性趋势模型\textstyle E{Z(x)} = \sum_{k=0}^p \beta_k f_k(x)。
IRFk-克里金法假设E{Z(x)}为x的一个未知多項式。
指标克里金法（Indicator kriging）使用指示函数而不是过程本身，以便估计转移概率。

** 多指标克里金法（Multiple-indicator kriging）是指标克里金法的一个版本，处理一系列指标。最初，MIK作为一种能更准确估算整体全球矿床浓度或品位的新方法，展现了相当大的前景。然而，由于所使用的区块尺寸本质上很大，且缺乏采矿规模的分辨率，这些优势已被建模中其他固有的实用性问题所抵消。在这种情况下，条件模拟正迅速成为被接受的替代技术。

析取克里金法（Disjunctive kriging）是克里金法的非线性推广。
对数正态克里金法通过对数对正数据进行插值。
潜在克里金法假设在的潜在水平（第二阶段）进行各种克里金处理，以产生空间功能预测。当分析空间功能数据{ (y_i, x_i, s_i) }{i=1}^n时，该技术非常有用，其中y_i = (y{i1} ,y_{i2}, \cdots, y_{iT_i})^\top是跨越T_i周期的時間序列数据，x_i = (x_{i1}, x_{i2}, \cdots, x_{ip})^\top是包含p个协变量的向量，而s_i = (s_{i1}, s_{i2})^\top是第i个受试者的空间位置（经度、纬度）。
协同克里金法（Co-kriging）表示对来自多个来源且不同数据源之间存在关系的观测数据进行联合克里金处理。协同克里金法也可以在贝叶斯方法中实现。
贝叶斯克里金法（Bayesian kriging）偏离了对未知系数和超参数的优化，从贝叶斯角度来看，优化被理解为最大似然估计。相反，系数和超参数是从它们的期望值中估计出来的。贝叶斯克里金法的一个优点是，它允许量化克里金仿真器的证据和不确定性。如果使用仿真器来传播不确定性，可以通过比较仿真器不确定性与总不确定性来评估克里金仿真器的质量（另见贝叶斯）。贝叶斯克里金法也可以与协同克里金法混合使用。

普通克里金法
未知值Z (x_0)被解释为位于x_0的随机变量，邻近样本的值Z(x_i),\ i = 1, \ldots, N也是如此。估计量\hat{Z}(x_0)也被解释为位于x_0的随机变量，是变量线性组合的结果。

克里金法旨在最小化在无偏条件下估计Z(x_0)的以下误差的均方值：\epsilon(x_0) = \hat{Z}(x_0) - Z(x_0) =
\begin{bmatrix}
W^T & -1
\end{bmatrix} \cdot
\begin{bmatrix}
Z(x_1) & \cdots & Z(x_N) & Z(x_0)
\end{bmatrix}^T =
\sum^N_{i=1} w_i(x_0) \times Z(x_i) - Z(x_0)前面提到的两个质量标准现在可以用新随机变量\epsilon(x_0)的均值和方差来表示：

无偏性
由于随机函数是平稳的，E[Z(x_i)] = E[Z(x_0)] = m，权重之和必须为1，以确保模型是无偏的。可以看作如下：E[\epsilon(x_0)] = 0 \Leftrightarrow \sum^N_{i=1} w_i(x_0) \times E[Z(x_i)] - E[Z(x_0)] = 0\Leftrightarrow m \sum^N_{i=1} w_i(x_0) - m = 0 \Leftrightarrow \sum^N_{i=1} w_i(x_0) = 1 \Leftrightarrow \mathbf{1}^T \cdot W = 1最小方差

两个估计量可以都有E[\epsilon(x_0)] = 0，但围绕其均值的分散程度决定了估计量质量的差异。为了找到具有最小方差的估计量，我们需要最小化E[\epsilon(x_0)^2]。 \begin{align}
\operatorname{Var}(\epsilon(x_0)) &= \operatorname{Var}\left(\begin{bmatrix} W^T & -1 \end{bmatrix} \cdot
\begin{bmatrix} Z(x_1) & \cdots & Z(x_N) & Z(x_0) \end{bmatrix}^T\right) \\
&= \begin{bmatrix} W^T & -1 \end{bmatrix} \cdot
\operatorname{Var}\left(\begin{bmatrix} Z(x_1) & \cdots & Z(x_N) & Z(x_0) \end{bmatrix}^T\right) \cdot
\begin{bmatrix} W \\ -1 \end{bmatrix}
\end{align}有关详细解释，请参阅协方差矩阵。 \operatorname{Var}(\epsilon(x_0)) = \begin{bmatrix} W^T & -1 \end{bmatrix} \cdot
\begin{bmatrix}
\operatorname{Var}_{x_i} & \operatorname{Cov}_{x_ix_0} \\
\operatorname{Cov}_{x_i x_0}^T & \operatorname{Var}_{x_0}
\end{bmatrix} \cdot
\begin{bmatrix} W \\ -1 \end{bmatrix}其中字母\left\{\operatorname{Var}_{x_i}, \operatorname{Var}_{x_0}, \operatorname{Cov}_{x_ix_0}\right\}代表 \left\{\operatorname{Var}\left(\begin{bmatrix} Z(x_1) & \cdots & Z(x_N) \end{bmatrix}^T\right),
\operatorname{Var}\big(Z(x_0)\big),
\operatorname{Cov}\left(\begin{bmatrix} Z(x_1) & \cdots & Z(x_N) \end{bmatrix}^T,
Z(x_0)\right)\right\}

一旦定义了在Z(x)的所有分析领域中有效的协方差模型或变异函数（C(\mathbf{h})或\gamma(\mathbf{h})），我们就可以写出任何估计量的估计方差表达式，它是样本之间协方差以及样本与待估计点之间协方差的函数： \begin{cases}
\operatorname{Var}\big(\epsilon(x_0)\big) =
W^T \cdot \operatorname{Var}_{x_i} \cdot W -
\operatorname{Cov}_{x_ix_0}^T \cdot W -
W^T \cdot \operatorname{Cov}_{x_ix_0} +
\operatorname{Var}_{x_0}, \\
\operatorname{Var}\big(\epsilon(x_0)\big) =
\operatorname{Cov}(0) +
\sum_i \sum_j w_i w_j \operatorname{Cov}(x_i,x_j) -
2 \sum_iw_i C(x_i,x_0).
\end{cases}从该表达式可以得出一些结论。估计方差：

一旦假设了均值和空间协方差（或变异函数）的平稳性，就无法量化任何线性估计量；
当样本与待估计点之间的协方差减小时，估计方差增大。这意味着，当样本距离x_0越远，估计效果越差；
随变量Z(x)的先验方差C(0)增大而增大；当变量的分散程度较小时，区域A内任何点的方差都较低；
不取决于样本的值，这意味着相同的空间配置（样本与待估计点之间具有相同的几何关系）在区域A的任何部分总是产生相同的估计方差；通过这种方式，方差无法衡量由局部变量产生的估计不确定性。

方程组
W = \underset{\mathbf{1}^T \cdot W = 1}{\operatorname{arg,min}}\left( W^T \cdot \operatorname{Var}{x_i} \cdot W - \operatorname{Cov}{x_ix_0}^T \cdot W - W^T \cdot \operatorname{Cov}{x_ix_0} + \operatorname{Var}{x_0} \right)解决此优化问题（参见拉格朗日乘数）会产生“克里金系统”：\begin{bmatrix}\hat{W}\\\mu\end{bmatrix} = \begin{bmatrix}
\operatorname{Var}_{x_i}& \mathbf{1}\\
\mathbf{1}^T& 0
\end{bmatrix}^{-1}\cdot \begin{bmatrix} \operatorname{Cov}_{x_ix_0}\\ 1\end{bmatrix} = \begin{bmatrix}
\gamma(x_1,x_1) & \cdots & \gamma(x_1,x_n) &1 \\
\vdots & \ddots & \vdots & \vdots \\
\gamma(x_n,x_1) & \cdots & \gamma(x_n,x_n) & 1 \\
1 &\cdots& 1 & 0
\end{bmatrix}^{-1}
\begin{bmatrix}\gamma(x_1,x^) \\ \vdots \\ \gamma(x_n,x^) \\ 1\end{bmatrix}附加参数\mu是一个拉格朗日乘子，用于在最小化克里金误差\sigma_k^2(x)时满足无偏条件。 ====
简单克里金法
的均值和包络线。]]
简单克里金法在数学上是最简单的，但通用性最差。它假设随机场期望值已知，并依赖于协方差函数。然而，在大多数应用中，期望值和协方差都无法预先得知。

应用“简单克里金法”的实际假设是：

场的广义平稳性（方差平稳）。
各处期望值为零：\mu(x) = 0。
已知的协方差函数c(x, y) = \operatorname{Cov}\big(Z(x), Z(y)\big)。

协方差函数是一个关键的设计选择，因为它规定了高斯过程的属性，从而决定了模型的行为。协方差函数编码了诸如平滑度和周期性之类的信息，这些信息反映在生成的估计值中。一种非常常见的协方差函数是平方指数函数，它非常利于平滑函数估计。因此，在许多现实世界的应用中，它可能会产生较差的估计，特别是当真实的底层函数包含不连续性和快速变化时。

方程组
“简单克里金法”的克里金权重没有无偏条件，由“简单克里金方程组”给出： \begin{pmatrix} w_1 \\ \vdots \\ w_n \end{pmatrix} =
\begin{pmatrix}
c(x_1, x_1) & \cdots & c(x_1, x_n) \\
\vdots & \ddots & \vdots \\
c(x_n, x_1) & \cdots & c(x_n, x_n)
\end{pmatrix}^{-1}
\begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}这类似于Z(x_0)对其他z_1, \ldots, z_n的线性回归。

估计
通过简单克里金法进行的插值由下式给出： \operatorname{Var}\big(\hat{Z}(x_0) - Z(x_0)\big) =
\underbrace{c(x_0,x_0)}_{\operatorname{Var}\big(Z(x_0)\big)} -
\underbrace{\begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}'
\begin{pmatrix}
c(x_1,x_1) & \cdots & c(x_1,x_n) \\
\vdots & \ddots & \vdots \\
c(x_n,x_1) & \cdots & c(x_n,x_n)
\end{pmatrix}^{-1}
\begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}}_{\operatorname{Var}\big(\hat{Z}(x_0)\big)}克里金误差由下式给出： \operatorname{Var}\big(\hat{Z}(x_0) - Z(x_0)\big) =
\underbrace{c(x_0,x_0)}_{\operatorname{Var}\big(Z(x_0)\big)} -
\underbrace{\begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}'
\begin{pmatrix}
c(x_1,x_1) & \cdots & c(x_1,x_n) \\
\vdots & \ddots & \vdots \\
c(x_n,x_1) & \cdots & c(x_n,x_n)
\end{pmatrix}^{-1}
\begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}}_{\operatorname{Var}\big(\hat{Z}(x_0)\big)},由此得到高斯-马尔可夫定理的广义最小二乘版本（Chiles & Delfiner 1999, p. 159）： \operatorname{Var}\big(Z(x_0)\big) = \operatorname{Var}\big(\hat{Z}(x_0)\big) + \operatorname{Var}\big(\hat{Z}(x_0) - Z(x_0)\big)贝叶斯克里金法

另见贝叶斯

属性

克里金估计是无偏的：E[\hat{Z}(x_i)] = E[Z(x_i)]。
克里金估计尊重实际观测值：\hat{Z}(x_i) = Z(x_i)（假设没有发生测量误差）。
如果假设成立，克里金估计\hat{Z}(x)是Z(x)的最佳线性无偏估计量。然而（例如Cressie 1993）：
与任何方法一样，如果假设不成立，克里金法可能效果不佳。

** 可能存在更好的非线性及/或有偏方法。
** 当使用错误的变异函数时，无法保证任何属性。然而，通常仍能获得“良好”的插值。
** “最佳”并不一定意味着“好”：例如，在没有空间依赖性的情况下，克里金插值仅与算术平均值一样好。

克里金法提供\sigma_k^2作为精确度的衡量标准。然而，该衡量标准依赖于变异函数的正确性。

应用
尽管克里金法最初是为地质统计学应用而开发的，但它是一种通用的统计插值方法，可以应用于满足适当数学假设的随机场采样数据的任何学科。它可以用于已收集空间相关数据（二维或三维）的情况，并且需要估算实际测量值之间位置（空间间隙）的“填充”数据。

迄今为止，克里金法已用于多种学科，包括：

环境科学
水文地质学
采矿业
自然资源
遥感

集成电路分析与优化
微波器件建模
天文學
页岩油井产油曲线预测

计算机实验的设计与分析
在工程学中，另一个非常重要且快速增长的应用领域是对作为确定性计算机模拟响应变量的数据进行插值，例如有限元法（FEM）模拟。在这种情况下，克里金法被用作元建模工具，即在设计的计算机实验集上建立的黑箱模型。在许多实际工程问题中，例如金属材料成形过程的设计，单次FEM模拟可能长达数小时甚至数天。因此，设计并运行有限次数的计算机模拟，然后使用克里金插值器快速预测任何其他设计点的响应，效率会更高。因此，克里金法经常作为所谓的代理模型，实现在优化程序中。在混合整数输入的情况下，也可以使用基于克里金法的代理模型。

评论 (0)