克里金法

统计学专业术语

克里金法（Kriging）是依据协方差函数对随机过程/随机场进行空间建模和预测（插值）的回归算法。在特定的随机过程，例如固有平稳过程中，克里金法能够给出最优线性无偏估计（Best Linear Unbiased Prediction, BLUP），因此在地统计学中也被称为空间最优无偏估计器（spatial BLUP）。

历史

克里金法的命名来自于南非金矿工程师丹尼·克里格（Danie G. Krige），以纪念其使用回归方法对空间场进行预测的开创性研究。克里金法（普通克里金）的提出者为法国统计学家乔治斯·马瑟伦（Georges Matheron），在其1963年发表的著作Principles of geostatistics中，马瑟伦将克里金法定义为“对已知样本加权平均以估计平面上的未知点，并使得估计值与真实值的数学期望相同且方差最小的地统计学过程”

（原文）“ … and must have the same average value within the whole large field … the (weights) have such values that estimation variance of by ,… should take the smallest possible value.”

马瑟伦在提出克里金法的同时使用了BLUP理论，为克里金法的发展奠定了基础。但在这一时期，不同领域的诸多研究都得到了等价或类似于克里金的估计方法。

在大气科学领域，前苏联气象学家列夫·舍米诺维奇·加丁（Лев Семено́вич Гандин）在其1963年发表的著作“气象场的客观分析（Objective analysis of meteorological field）”中开展了与马瑟伦相同的工作。在该著作1965年的英语译本中，简单克里金被称为“最优插值（optimum interpolation）”、普通克里金被称为“归一化权重最优插值（optimum interpolation with normalization of weighting fuctors ）”、协同克里金被称为“最优空间场匹配（optimum matching of flelds）”。加丁在其著作中同样发展了克里金法的BLUP理论并讨论了克里金法在气象领域的应用。

在高斯过程领域，安德雷·柯尔莫哥洛夫（Андрей Николаевич Колмогоров）和罗伯特·维纳（Robert Wiener）在二十世纪40年代开展了包括线性预测和BLUP理论在内的工作，并将其应用于时间序列数据中，被后世称为维纳-柯尔莫哥洛夫滤波器（Wiener-Kolmogorov filter）的信号处理技术是与克里金法十分接近的求解系统。随后在70年代，高斯过程理论和贝叶斯推断被逐步应用于空间场的研究，促进了克里金法的发展。

理论

固有平稳过程

克里金法在对空间场进行预测时，将空间场视为随机过程（stochastic process）的推广，即随机场（random field）。具体而言，若对指数集有随机过程且为一拓扑空间（topological space），则被称为随机场。克里金法中随机场所对应的指数集通常为地理坐标，而随机场内每一个点的测度都是一个随机变量，服从特定的概率分布。

根据地理学第一定律（the first law of geography），地理空间上的所有值都是互相联系的，且距离近的值具有更强的联系。随机场使用协方差函数对上述结论进行刻画，对应高斯过程回归（GPR）理论中的“核函数（kernel function）”。使用克里金法需要随机场满足两个假设：

满足上述假设的随机过程被称为固有平稳过程（intrinsically stationary process），二阶平稳过程（second-order stationary process）是其特例。平稳高斯过程不仅是严格的平稳过程，而且在许多问题中是一个理由充分的假设，因此克里金法有部分研究完全基于高斯随机场展开。

克里金法通常假设固有平稳过程是各向同性（isotropy）的，即其协方差函数仅是点间欧氏距离（Euclidean distance）的函数。各向同性随机场可以直接使用变异函数（variogram）进行建模，简化了克里金法的求解步骤。此外，一些特定类型的各项异性，例如几何各向异性（geometric anisotropy）可以通过坐标变换转化为各向同性。

类似于高斯过程回归，克里金法不要求样本服从特定的概率分布，但实践中当样本没有偏斜数据时，克里金法往往有好的效果。

BLUP

克里金法是对已知函数进行加权平均估计未知函数的方法，其预测理论接近于线性回归，因此也有类似高斯-马尔可夫定理（Gauss-Markov theorem）的BLUP理论。

将随机场中变量的估计表示为包含随机误差的线性系统，则BLUP可表示为选择线性系统参数使估计值和真实值的方差最小：

式中为未知点，为随机场的样本，为权重系数，通常被称为克里金权重（Kriging weights）。由方差定义可知，当估计值和真实值的数学期望相同时，两者的方差最小：

使用上述BLUP条件求解的权重系数包含样本点与未知点间的协方差函数（或变异函数），因此只有随机场的协方差函数已知，或随机场为固有平稳过程时，克里金法才能给出BLUP。

算法

普通克里金

普通克里金（Ordinary Kriging, OK）是最早被提出和系统研究的克里金法，并随着地统计学的发展衍生出一系列变体和改进算法。普通克里金是一个线性估计系统，适用于任何满足各向同性假设的固有平稳随机场。

各向同性假设下的固有平稳随机场中，数学期望与其位置无关，且协方差仅是点间距离的函数。通常随机场的协方差函数是未知的，需要使用变异函数作为近似，此时变异函数也仅与点间距离有关：

定义为n个样本的对应值，则普通克里金问题和克里金方差有如下表示：

式中是在未知点处的估计，是点和间的协方差函数。由BLUP理论易证明普通克里金的无偏估计条件是所有权重系数之和为1：。由此可使用拉格朗日乘数法构造如下求解函数并得到普通克里金问题的方程组：

该方程组有时也被称为普通克里金系统（ordinary Kriging system），包含n+1个方程以求解n个权重系数。常见的求解方式是将克里金系统写为矩阵形式并对矩阵求逆。求解后以矩阵形式表示的克里金权重为：

式中为协方差矩阵，为未知点和样本间协方差组成的列向量，为n个1组成的列向量。将所得权重带入先前公式中可以得到普通克里金的无偏估计：

由上式易知，随机场的数学期望为。若随机场的数学期望已知（通常假设为0），则普通克里金退化为简单克里金（simple Kriging）。由于固有平稳随机场的数学期望处处相等，因此简单克里金自身满足BLUP条件，容易解得其克里金权重为。

在高斯随机场中，普通克里金和简单克里金的估计值有如下分布的置信区间：

式中为标准正态分布的分位函数，即累积分布函数的反函数。

改进算法

泛克里金（Universal Kriging, UK）

在已知随机场是各项同性的固有平稳过程和漂移量（drift）的叠加时，可以使用泛克里金对随机场进行建模。泛克里金假设漂移量是一系列已知解析函数的线性组合：

式中是随机场的漂移，是满足各向同性假设的固有平稳随机场。最常见的漂移是线性函数，对应具有线性趋势的随机场。泛克里金的问题表述和克里金方差与普通克里金相同，其无偏估计条件有如下表述：

当且时，泛克里金与普通克里金的无偏估计条件等价，因此普通克里金可以视为泛克里金的一个特例。将上式作为拉格朗日乘子可得泛克里金的求解系统：

将求解所的权重带入先前公式中可以得到泛克里金在的估计值。在高斯随机场中，泛克里金和普通克里金以相同的方法估计置信区间。

协同克里金（Co-Kriging, CK）

协同克里金是克里金法在处理多变量问题时的改进，其中需要建模的随机场称为主变量，参与建模的其它随机场称为协变量。协同克里金可以有任意个数的协变量，但主变量和协变量必须具有相关性（correlation）且均是满足各向同性假设的固有平稳随机场。协同克里金可以使用互变异函数（cross-covariogram）估计随机场间的互协方差函数（cross-covariance），但要求后者是对称的，即。

若主变量拥有M个协变量且每个变量拥有个样本，则协同克里金的优化问题和克里金方差有如下表示：

式中为主变量的协方差函数，为主变量与协变量间的互协方差函数。协同克里金的无偏估计条件为主变量权重系数之和为1,协变量权重系数之和为0：

类似普通克里金法，使用拉格朗日乘数法可以构造协同克里金系统求解系数 :

协同克里金的估计结果有如下分布的置信区间：

上述算法为基于普通克里金发展的协同克里金，泛克里金也有对应的协同算法，被称为“泛协同克里金（universal co-Kriging）”。协同克里金和泛协同克里金通常要求协变量拥有充足的样本。若协变量的样本量不足以计算互变异函数。可计算伪互变异函数（pseudo cross-variogram）作为近似。

析取克里金（Disjunctive Kriging, DK）

普通克里金是随机场的BLUP，但在随机场和指数集之间存在非线性关系时，线性估计结果往往不是最优的。析取克里金将普通克里金中的权重系数推广为函数，从而实现了对随机场的非线性估计。析取克里金可以定义为如下的优化问题：

式中函数为预先给定的非线性函数。最常见地，给定指示函数时，析取克里金又被称为指示克里金（Indicator Kriging, IK）。析取克里金求解给定函数使得成为真实值在由构成的向量空间中的正交投影，由此可得其问题表述为：

析取克里金要求样本集是同因子的（isofactorial），但应用中通常直接假设样本集服从联合正态分布，此时使用阶埃尔米特多项式（Hermite polynomial）对函数进行展开可将上式转化为：

式中为k阶埃尔米特多项式，为满足如下关系的待定参数：

为和的相关系数，为k阶埃尔米特多展开系数：

由埃尔米特多展开后的析取克里金问题，可得其克里金方差为：

析取克里金在每个埃尔米特展开上使用普通克里金求解和具体的求解过程依赖于数值计算，考虑计算的复杂程度，通常将埃尔米特展开级数限制在100以内。

混合算法

回归克里金（regression-Kriging）

回归克里金是广义线性模型（Generalized Linear Model, GLM）和克里金法相结合的算法，也是最常见的混合算法。回归克里金首先使用GLM估计空间场中的系统性效应（determinstic effect），随后使用克里金法估计由回归残差构成的随机场：

回归克里金考虑了空间场的趋势，因此和泛克里金较为相似，但后者是基于随机场假设的空间估计，而回归克里金则将线性趋势和随机过程完全分开。

神经网络克里金（neural Kriging）

神经网络克里金可见于大气科学研究中，一些文献也将其称为“直接神经网络残差克里金（Direct Neural Network Residual Kriging, DNNRK）”。神经网络克里金与回归克里金在逻辑上类似，即首先使用各类人工神经网络（Artificial Neural Network, ANN）算法对空间场进行建模，随后使用克里金法估计残差。

贝叶斯克里金（Bayesian Kriging）

贝叶斯克里金是使用贝叶斯推断（Bayesian inference）的克里金法的统称。贝叶斯克里金使用由超参数（hyper-parameter）的先验（prior）定义克里金系统的参数（权重系数）并估计其后验（posterior）。贝叶斯克里金的先验通常为正态分布和Gamma分布。

常见的贝叶斯克里金框架包括Kitanidis算法和Handcock算法。其中Kitanidis算法要求随机场的协方差除尺度外已知，在实际应用中难以满足；Handcock算法不要求协方差函数已知，而是将其表示为马顿核（Matérn kernel）形式的高斯过程先验并使用极大似然估计（Maximum Likelihood Estimation, MLE）求解超参数。两种算法都以高斯随机场为前提，且后者是与高斯过程回归等价的方法。对非高斯随机场，Handcock算法有贝叶斯高斯变换（Bayesian Transformed Gaussian, BTG）算法。而对更一般的各向异性随机场，以及有同时包含空间和时间观测样本的情形，可以使用等级贝叶斯克里金（hierarchical bayesian Kriging）。该方法可以对非平稳随机场使用，其中协方差函数完全由超参数进行模拟。

性质

克里金法是一个确切估计器（exact estimator），由克里金法的求解系统可知，其估计的随机场在样本点的取值与对应观测值一致。这一性质的优势在于克里金法的估计总是与观测值接近，不会与实际情况相距太远；但估计值总是介于观测值之间，不能模拟剧烈的的变化。

由于克里金法在应用中使用拟合经验变异函数的方式估计随机场的协方差，而变异函数模型除块金（原点）外都是连续函数，因此克里金法对随机场的估计是平滑的。这意味着在偏斜样本中克里金法对异常值敏感，因此对偏斜样本进行正态变换，例如Box-Cox变换是常见的做法。

克里金法仅在随机场为固有平稳过程时提供BLUP，因此在使用克里金法前有必要通过样本考察随机场的稳定性，一个常见的方法是绘制沃罗诺伊图（Voronoi diagram），计算标准差并观察是否有局地异常值。

在克里金法中，所有样本都会参与对未知点的估计，且在变异函数随距离递增的情形下，与未知点距离近的样本影响更大。实际应用中为突出随机场的局地特征，可以人为设定点的影响范围，不考虑范围外样本的影响。这一性质也意味着若随机场某处没有临近的样本，则其克里金方差会增大，因此克里金法仅适用于格点数据的内插值。普通克里金的外插结果随点距离的增大趋近于样本平均值；泛克里金会在外插时将漂移量无限外推。

与高斯过程回归（GPR）的关系：GPR是与克里金法相近的非参数模型，常见于机器学习领域。若协方差函数的形式等价，简单克里金和普通克里金的输出与GPR在正态似然下输出的数学期望相同。GPR与克里金法的不同点在于，前者假设随机场为高斯过程并给出其对测试样本后验的完整分布；后者假设随机场为固有平稳过程并给出其对测试样本的最优无偏估计。此外，一些贝叶斯克里金方法，例如Handcock算法与GPR是等价的，可参见Handock and Stein (1993)。

应用

克里金法被广泛用于各类观测的空间插值，例如地质学中的地下水位和土壤湿度的采样；环境科学研究中的大气污染（例如臭氧）和土壤污染物的研究；以及大气科学中的近地面风场、气温、降水等的单点观测。

克里金法在工程问题的数值试验中可作为代理模型（surrogate model）对有限的模拟结果进行插值。具体而言，若对问题全局使用确定性模拟方法（deterministic computer simulations），例如有限元方法会占用大量计算资源而无法（快速）实现时，可以仅模拟局部个别点的结果并使用克里金法插值到全局。

包含克里金法的编程模块

作为一种典型的地统计学算法，克里金法被广泛收录于各类代码库中，以下对其进行部分归纳：

全国各地天气预报查询

上海市

云南省

云南省

云南省

云南省