点参考数据及克里金法

【阅读建议】 点参考数据的空间预测和模拟问题，大致有传统克里金法和目前应用比较广泛的基于似然的方法。本文主要介绍源于地统计学的传统克里金方法，一来掌握空间统计中的基础方法，二来便于与后面几篇文章中提到的高斯过程之间建立联系。克里金方法在对空间随机场作出本征平稳假设的情况下，利用参数化的变异函数对不同位置处随机变量的偏差之间存在的空间结构（相关性）进行建模，利用有限样本点的最大似然求解最优参数，并将其用于预测任务。

【引文信息】

• [1] 史舟， 李艳编， 地统计学在土壤学中的应用. Beijing: Zhong guo nong ye chu ban she， 2006.
• [2] 王政权， 地统计学及在生态学中的应用. Bei jing: Ke xue chu ban she， 1999.

1 引言

空间数据的获取通常具有一定的成本，是进行空间分析的基础与起源。为了提高研究结论的精度，我们希望能够获取研究区域内更多、更全面的精确空间属性数据信息。然而，在实际研究工作中，由于人力成本、资源等外部条件限制，我们无法对全部未知区域加以采样与测量，而往往只能得到研究区域内有限数量的采样点及其相关属性数据。因此，往往可以考虑选取合适的空间采样点，利用一定数学模型，依据已知采样点各自对应属性数据对研究区域所有位置的未知属性信息加以预测（Predict）。

在空间信息技术领域，空间插值（Spatial Interpolation）是实现此类空间预测的另一种称呼，空间插值是一种将离散采样点测量数据转换为连续数据表面的常用方法，包括内插（Interpolation）和外推（Extrapolation）两种形式。一般地，对样本点范围以内（即所有采样点最大外接矩形内部）的空间进行插值，才可称作“内插”（部分文献亦直接使用术语 “插值”）；反之则称 “外推” 或 “预测（Forecast）”，往往认为外推的结果误差较大。

空间插值理论及其方法基于著名的 “地理学第一定律（Tobler’s First Law of Geography）”，即一般地，距离越近的地物具有越高的相关性。这一至今已产生深远影响的地学定律最早由美籍瑞士地理学家 Waldo R. Tobler 教授于 1970 年提出。

在各方法所对应的数学计算原理层面，空间插值一般可以分为确定性插值方法（Deterministic Interpolation）与地统计插值方法（Geostatistics，亦称非确定性插值方法）两种。前者基于研究区域内各信息点之间相似程度或整个表面的平滑程度，创建连续的拟合表面；根据插值计算时所使用样本点范围，又可进一步分为全局插值法与局部插值法。

本文的重点是地统计插值方法，此类方法基于研究区域内各信息点的综合统计学规律，其核心要素有两个：一是变异函数，用于建模空间自相关性以及结构，实现属性的空间自相关性量化；二是克里金法，以变异函数理论与结构分析为基础，通过克里金预测方程创建得出空间过程表面。

2 基础理论

和传统统计技术一样，地统计学也是建立在概率论基础上的。在经典统计学中的一些基本概念例如均值、方差和协方差等，在地统计学中也存在。但经典统计学和地统计学的区别也是真真切切存在的。下面介绍地统计学独有的概念。

2.1 空间过程、区域化变量与随机函数

（1）空间过程

地统计学将研究域内的每一个感兴趣的属性视为一个空间随机过程（或空间随机场），域内任意点处的属性取值都是一个随机变量。直观地理解，就是在空间域内的一个具有一定厚度的表面（Surface），其中每个点的均值一起构成了一个确定性表面，而厚度则代表了各点属性取值的不确定性。

如果对任意点处的随机变量作出高斯分布的假设，则该空间过程是二维或三维空间场景中的一个高斯过程。

与传统统计中通过多次独立实验获得随机变量的统计特性不同，在地统计中，$n$ 个不同空间测站获得的测量值，被视为空间过程某一次随机实验的结果，注意是单次实验结果。

看到表面（Surface）这个词，很容易让人联想到趋势面分析。趋势面分析也是拟合表面的一种形式，但其思路与空间过程的思路存在着本质的不同。趋势面分析源自典型的频率派方法，它认为每个样本点都是来自表面模型的一次独立采样 $(\mathbf{s}, z)$，拟合任务是根据多次重复的观测数据来重建表面模型 $z=f(\mathbf{s})$。在趋势面分析中，任意位置 $\mathbf{s}$ 处的属性值被视为随位置而变化的某个确定值，不同位置处的不确定性主要来源于测量误差，且相互独立。

（2）区域化变量

如果令域内某个位置 $\mathbf{s}$ 处的 $z$ 属性为随机变量，则可以发现 $z$ 会随着 $\mathbf{s}$ 而变化，因此 $z$ 被称为区域化变量，即不同位置是不同的变量。区域化变量可以表示为 $z(\mathbf{s}),\mathbf{s} \in D$ 。

如果假设不同位置处的 $z$ 变量完全独立，则理论上所有位置处可以有自己独立的概率分布形式和均值、方差等统计量。但自然世界的事实恰恰与这种独立性假设相反，根据地理学第一定律，空间事物之间往往体现着依赖性，而地统计学正是基于该定律，利用样本点之间的空间结构和相关性，使得我们能够利用单次实验结果和有限样品就能实现对整个空间过程的建模，并将之用于空间预测（插值）。

（3）随机函数

可以用随机函数 $Z(\mathbf{s})$ 的概念来同时表达空间过程的随机性质和区域化变量的区域性质。其中随机的概念表明在任意给定位置 $\mathbf{s}_i$ 处的属性 $Z(\mathbf{s}_i)$ 是一个随机变量，具有随机特性；而函数的概念表明，这种随机特性会随着 $\mathbf{s}$ 而变化，不同的 $\mathbf{s}$ 值对应不同的随机变量，并且不同位置处的随机变量在一定程度上依赖于所处的区域位置和空间结构。

由此可见， 随机函数可以被视为空间过程的一种形式化表达方式，而区域化变量可以被视为随机函数（空间过程）在某个位置处的一个实现。而实验在 $\mathbf{s}_i$ 处采集的样品，则是区域化变量 $Z(\mathbf{s}_i)$ 的一次采样结果。

注解：

有关贝叶斯分层建模方法的文献中，会对数据、过程、参数进行层次化建模，在此类文献中符号 $Z$ 通常指观测数据，而非空间过程，这与此处的符号可能有些不一致。

2.2 平稳性假设

地统计学认为空间过程（随机函数）中不同点处的随机变量具有不同的分布，其解空间非常庞大，而现实中仅能得到有限位置处随机变量的一个或几个稀少的观测，无法根据这些有限的观测数据求解无穷大的解空间。因此，必须要对空间过程作出一些必要的假设。其中一个重要假设就是平稳性假设，即假设不同位置的随机变量，其概率分布之间存在某种不变性。

（1）严格平稳性假设（强平稳性）

严格平稳假设是指过程中所有随机变量具有完全相同的概率分布，随机性质完全相同。严格平稳假设过于严格了，几乎不适用于现实世界中的任何情况，因此，可以通过放宽一些平稳条件，从而形成一些更有意义的平稳性假设。

（2）二阶平稳性假设（弱平稳性假设）

二阶平稳性假设又称弱平稳假设，该假设不强制要求处处具有相同的概率分布，而是放宽为处处具有相同的二阶统计量，且任意两个随机变量之间的协方差仅仅依赖于其二者之间的距离与方向，而与其具体位置无关。
上述两个条件分别表示为：

条件 1: $E[Z(\mathbf{s})] = m, \quad \mathbf{s} \in D$

条件 2: $\operatorname{Cov}[Z(\mathbf{s}),Z(\mathbf{s+h})] = C(h), \quad \mathbf{s},\mathbf{s+h} \in D$

其中，$E[Z(\mathbf{s})]$ 为区域化变量 $Z(\mathbf{s})$ 的数学期望， $\operatorname{Cov}[Z(\mathbf{s}),Z(\mathbf{s+h})]$ 为区域化变量 $Z(\mathbf{s})$ 与 $Z(\mathbf{s+h})$ 的自协方差函数，$C(\mathbf{h})$ 为仅与 $\mathbf{h}$ 有关的协方差取值，$m$ 为一常数，$\mathbf{h}$ 为步长（或距离、滞后）。条件 2 表明任意两个区域化变量之间的自相关性仅与两者之间的步长有关。

协方差函数依赖于测量随机变量时使用的量纲和尺度，将其进行归一化转换为 自相关函数 可能更容易比较和理解：$\rho(\mathbf{h}) = c(\mathbf{h}/C(0))$。此处的 $C(0)$ 表示 $\mathbf{h}=\mathbf{0}$ 时的协方差，就是先验方差 $\sigma^2$。

上述二阶平稳性假设是针对整个研究区域范围而言。若区域化变量仅仅在某个有限子区域中满足上述条件，则称此区域化变量满足准二阶平稳性假设（Quasi Second Stationary Assumption）。准二阶平稳性假设可以视作一种折中方案，既考虑到平稳范围的大小，又顾及到有效数据的数量。

（3） 本征假设性假设

本征假设（Intrinsic Hypothesis）又称内蕴假设。该假设对二阶平稳性进一步放宽，认为只要区域化变量的增量满足二阶平稳性条件即可。即：

条件 1： $E[Z(\mathbf{s})-Z(\mathbf{s+h})]= 0, \quad \mathbf{s,s+h} \in D$

条件 2: $\operatorname{Var}[Z(\mathbf{s}) - Z(\mathbf{s+h})]=E[\{Z(\mathbf{s}) - Z(\mathbf{s+h})\}^2] = 2 \gamma(\mathbf{h})$

条件 1 表明，在整个研究区域内，区域化随机变量的增量均值为 $0$，或等价的，区域化变量的均值处处相等；条件 2 表明该增量存在方差函数，该方差函数只依赖于步长 $\mathbf{h}$，而与所处位置无关。上式中，$\gamma(\mathbf{h})$ 是步长 $\mathbf{h}$ 上的半方差值，因此也被称为半方差函数，人们更习惯称之为变异函数。

由于二阶平稳假设条件在现实世界中也很难出现，因此本征假设将用户从二阶平稳假设的严格限制中解放了出来，使其适用范围得到了大大扩展。

同样的，本征假设亦是针对整个研究区域范围而言。若区域化变量仅仅某个子区域中满足上述条件，则称此区域化变量满足准本征假设（Quasi Intrinsic Hypothesis）。与准二阶平稳性假设类似，准本征假设亦可视作一种折中方案，其同样既考虑到了本征假设对应范围的大小，又顾及到了有效数据的数量。

可以看出，本征假设运用于区域化变量的增量 $Z(\mathbf{s+h}) - Z(\mathbf{s})]$ 上，而不是运用到 $Z(\mathbf{s})$ 本身。如果一个本征平稳的空间过程在 $Z(\mathbf{s+h})$ 和 $Z(\mathbf{s})$ 之间的协方差函数也存在且仅依赖于步长 $\mathbf{h}$ ，则 $Z(\mathbf{s})$ 也是二阶平稳的。也就是说，二阶平稳性过程肯定是本征平稳的，而本征平稳性过程不一定是二阶平稳的。

本征假设是地统计学中对随机函数的基本假设。

（4） 说明

事实上，本征假设的 条件 1 （即 $E[Z(\mathbf{s})-Z(\mathbf{s+h})]= 0$ 或者 $E[Z(\mathbf{s})] = m$ ）也很少能被一些空间过程（随机函数）满足。例如，春天的平均温度通常是南方高于北方，从北到南温度呈不断增加的趋势。在这种情况下，可以考虑将空间过程分解为一个大尺度趋势（或漂移）和一个中小尺度上的空间过程（只对各位置处的偏差建模）。也就是说，$Z(\mathbf{s})= \mu(\mathbf{s}) + \epsilon(\mathbf{s})$，其中仅 $\epsilon(\mathbf{s})$ 被视为随机过程。

3 变异函数及其拟合

根据上一节的介绍，本征平稳性假设下的变异函数（以及二阶平稳性假设下的协方差函数）决定了空间过程中区域化变量之间的空间结构和相关性，是地统计学的基石。对其进行科学合理的建模，并将模型应用于新位置处区域化变量的预测，是地统计最核心的内容。

3.1 与变异函数相关的统计概念

变异函数与协方差函数、自相关函数之间具有紧密的联系，理解他们之间的关系，有助于对理论方法的理解。

在二阶平稳性假设下，存在以下性质和关系：

（1）协方差函数满足： $|C(h)| \leq C(0) = \operatorname{Var}(Z(x))$

（2）协方差函数与自相关函数的关系：$\rho(\mathbf{h}) = C(\mathbf{h})/C(0)$

（3）变异函数与协方差函数的关系： $\gamma(\mathbf{h}) = C(0) - C(\mathbf{h})$

（4）变异函数与自相关函数的关系：$\gamma(\mathbf{h})=\sigma^2 (1-\rho(\mathbf{h}))$

在本征平稳性假设下：

虽然 $\gamma(\mathbf{h})$ 存在，但协方差函数 $C(\mathbf{h})$ 不一定存在，所以上述等式不一定成立。 显然变异函数能够在更多情况下表现出有效性，因此在地统计学中得到了比协方差函数更为广泛的应用。

二阶平稳性和本征平稳性可以通过变异函数图进行判别：二阶平稳随机过程的半方差会随着步长 $\mathbf{h}$ 增加而增大，到一定程度后该值不再有显著变化，基本保持在一个水平上不变（即基台值）；而本征平稳的随机过程可能不存在基台值。

左图：变异函数与协方差函数的关系；右图：协方差函数不存在时的变异函数

3.2 变异函数及其特征

变异函数，又称为半方差函数，用于描述区域化变量的空间变化特征与强度。根据上一节中的式，变异函数被定义为区域化变量增量平方的数学期望的 $\frac{1}{2}$；变异函数的取值与区域化变量样点所处位置 $\mathbf{s}$ 无关，仅与样点之间的步长（或距离、滞后）$\mathbf{h}$ 有关。

在大多数情况下（并非全部），区域化变量的经验变异函数曲线会呈现出 “先快速上升，再增速减缓，后趋于平稳” 的曲线特征。在地统计学中，此变异函数曲线具有几个十分重要的相关概念：

经验变异函数指根据实验数据获得的变异函数。

（1）块金（Nugget）

块金 代表区域化变量的随机性大小。由理论角度，在步长为 $0$（即滞后距为零）时，区域化变量的样点数值应当相等，因此在步长无限趋近于 $0$ 时，对应变异函数值也应当亦向 $0$ 趋近。但在实际研究中却发现，经验变异函数在步长为 $0$ 时，其取值并不为 $0$ ，而是一个大于 $0$ 的值。这一数值被称为块金常数。

地统计学加将块金效应归因于测量误差（独立的随机白噪声，具有固定的方差）和微尺度变化（由小于样本测量的最小距离的空间变化引起的随机性，由于通过样本无法对其建模，因此只能假设其为独立的、具有固定方差的白噪声）。

（2）基台（Sill）

基台 用于衡量区域化变量变化幅度的大小。当步长无限增大并到达某一程度后，经验变异函数若趋于平稳，则这一平稳水平所对应的变异函数值即为基台值。此情形下的空间过程是一个二阶平稳的随机过程。有些空间过程的半方差值会随着步长无限增长，即不存在基台值，则此时的空间过程可能是本征平稳的，但肯定不是二阶平稳的。

（3）变程（Range）

变程 用于衡量区域化变量自相关范围的大小。当步长无限增大并到达某一程度后，经验变异函数若趋于平稳，则此时对应的步长即为变程。由此可见，变程特指一个阈值，步长小于变程的样本存在空间自相关，而步长大于变程样本则不存在空间自相关。在实际工作中，当变异函数是渐进接近基台时，通常采用 $0.95$ 倍基台值对应的步长作为有效变程。

（4）偏基台（Partial Sill）

基台值与块金常数之间的差值，被称为 偏基台值（Partial Sill），用于衡量空间变异性的范围。

（5）块金系数

块金系数： 块金常数与基台值的比值，用于衡量空间变异性程度。

（6）孔穴效应

有些变异函数的半方差值在达到最大值后开始减小，然后又增大，呈现出某种周期波动的变化。此时变异函数的最大半方差值对应于协方差函数的最小值，两者之间看起来像一个孔穴，因此被称为孔穴效应。孔穴效应的出现，意味着随着距离增加，空间变化呈现一种有规律的重复。

（7）漂移

当变异函数以一个凹向上的形式接近原点时，表明其已经受到了漂移或趋势的影响，导致区域化变量的均值不再是常数，而是一个以位置 $\mathbf{s}$ 为自变量的函数，即 $E[Z(\mathbf{s})] = m(\mathbf{s})$，式中的 $m(\mathbf{s})$ 为漂移。此时，二阶平稳性和本征平稳性假设都不成立。解决此问题的方法是将空间过程分解为两部分：非平稳部分（宏观尺度的变化，用漂移 $m(\mathbf{s})$ 建模） 和平稳部分（微观尺度的变化，用$\varepsilon(\mathbf{s})$ 建模 ），通常漂移 $m(\mathbf{s})$ 是确定性的函数（如全局性的趋势面函数、局部的近邻加权函数等），而 $\varepsilon(\mathbf{s})$ 是其对应的残差。

3.3 变异函数的类型

根据区域化变量的变异函数特征，可以将变异函数分为不同类别加以区分和建模。

依据变异函数基台值的有无，可以将模型分为：

• 有基台值模型
  • 纯块金效应模型
  • 球状模型（Spherical Model）
  • 指数模型（Exponential Model）
  • 高斯模型（Cubic Model或Gaussian Model）
  • 线性有基台模型（Linear with Sill Model）
• 无基台值模型
  • 线性无基台值模型
  • 幂指数模型
  • 对数模型
• 孔穴效应模型
  • 基台值模型
  • 无基台值模型

套合模型：针对某种区域化变量而言，其在不同方向、不同滞后距情况下可能受到不同因素影响；套合模型可以很好解决这一问题。套合结构可以表示为多个变异函数之和，每一个变异函数均代表着某种方向或某一尺度中的变异性，从而对区域化变量的特征加以更好概括。

上述所有模型都可以被视为对变异函数的某种理论上的参数化假设，其参数可以根据观测数据推断得出。因此，对于一次地统计研究过程来说，重要的是根据观测数据给出的结果，科学合理地作出对变异函数的假设，并利用样本数据推断得出变异函数的参数，以完成后续预测任务。

3.4 变异函数的选择

通过对样本数据的探索性分析，研究人员可以尝试上述各种可能的变异函数模型。其大致工作流如下：

（1）获得半方差云图

拿到样品数据后，通过计算空间域 $D$ 内任意位置点对 $(\mathbf{s}_i,\mathbf{s}_j)$ 观测值之间的差异，并量化不同尺度上区域化变量的变异情况，这种差异可以用半方差值来量化计算：

$$\gamma(\mathbf{s}_i,\mathbf{s}_j) = \frac{1}{2}[z(\mathbf{s}_i) - z(\mathbf{s}_j)]^2$$

令半方差值以任意两点之间的向量 $\mathbf{h}_{ij}$ 为自变量，则有:

$$\gamma^*(\mathbf{h}_{ij}) = \frac{1}{2}[z(\mathbf{s}_i) - z(\mathbf{s}_j)]^2$$

如果假设这种变异具有各向同性，则有：

$$\gamma^*(\| \mathbf{h}_{ij} \|) = \frac{1}{2}[z(\mathbf{s}_i) - z(\mathbf{s}_j)]^2$$

此时我们可以根据步长 $\| \mathbf{h}_{ij} \|$ 绘制出所有样品间半方差值构成的云图：

所有样品位置对之间的半方差云图，据此可以大致了解变异函数的宏观情况

（2）生成经验半方差函数图

在实际工作中，通常将步长划分为不同的区间，并在区间内计算半方差的统计特性（如均值、方差等），进而得到经验半方差值和经验半方差图。 一种计算经验半方差值的简单方法是求均值：

$$\gamma^*(h) = \frac{1}{2N(h)} \sum \limits_{i=1}^{N(h)}[z(\mathbf{s}_i + h) - z(\mathbf{s}_i)]^2$$

其中 $N(h)$ 表示落入 $h$ 区间的半方差位置对数量。下图是对上面半方差云图求均值后绘制出的经验半方差函数图。

经验变异函数在原点处的性状非常重要，反映了区域化变量不同尺度的连续型。通常包括五种情况：

• 抛物线型：当 $\|\mathbf{h}\| \rightarrow 0$ 时，$\gamma(\mathbf{h}) \rightarrow A\|\mathbf{h}\|^2$（$A$ 为常数），$\gamma(\mathbf{h})$ 在原点处二次可导，并且随机函数本身在均方意义下可导。其反映出区域化变量具有高度空间连续性特征
• 线性型: 当 $\|\mathbf{h}\| \rightarrow 0$ 时，$\gamma(\mathbf{h}) \rightarrow A\|\mathbf{h}\|$。$\gamma(\mathbf{h})$ 在原点处保持连续。其反映出区域化变量具有平均意义上的连续性。
• 间断型: 当 $\|\mathbf{h}\| \rightarrow 0$ 时，$\gamma(\mathbf{h}) \rightarrow C(0)$。$\gamma(\mathbf{h})$ 在原点处间断，反映块金效应的存在。该性状说明变异函数连续性较差，但当步长增大时，又可逐渐变得比较连续。
• 随机型: 当 $\|\mathbf{h}\| \rightarrow 0$ 时，$\gamma(\mathbf{h}) \rightarrow C(0)$, 但当 $h$ 增大时，$\| \mathbf{h} \| \rightarrow \infty$，$\gamma(\mathbf{h})$ 仍然在 $C(0)$ 附近摆动。无论 $h$ 多么小, 区域化变量 $Z(\mathbf{s})$ 与 $Z(\mathbf{s+h}$ 之间总是不相关。这种类型被称纯块金效应型。它反映了区域化变量完全不存在空间相关的情况,或者说反映了变量是普通的随机变量。此时 $C(0)$ 等于先验方差，$C(0) = \operatorname{Var}(Z(\mathbf{s}))$。
• 过渡型: 当 $\|\mathbf{h}\| \rightarrow 0$ 时，$\gamma(\mathbf{h}) \rightarrow C(0)$；当 $\mathbf{h} = a$ 时，$\gamma(a) = C(0) + C$。这种变异函数是处于抛物线型和纯块金效应型之间的过渡型，即变异函数 $\gamma(h)$ 存在块金常数 $C(0)$ 和基台值 $C(0)+C$。有时也将其中的 $C$ 称为 “拱高”。 当存在块金常数 $C(0)>0$ 时, 基台值等于块金常数加上拱 高, 当块金常数 $C(0)=0$ 时, 基台值就等于“拱高”。过渡型是实际研究工作中最常遇到的一种类型。

（3）选择理论半方差函数

实际上理论变异函数模型 $\gamma(h)$ 是未知的，要从有效的空间取样数据中去估计。对各种不同的 $h$ 值可计算出一系列 $\gamma^*(h)$ 值, 因此, 要用一个理论模型去拟合这一系列的 $\gamma^*(h)$ 值。这本质上是根据上述的经验数据分析，选择 第 4.3 节 中合适的变异函数。有关各类变异函数的详细信息，请参考相关文献。

但有两种现实性的情况是必须考虑的，一是变异函数的结构往往不是单纯的一种，而是多层次结构相互叠加在一起，例如，土壤中某一元素的含量 在 $h \approx 0$ 时，其差异主要由测量误差造成；在 $h \leq 1m$ 时，可能还要加上土壤水分造成的变异成分；在 $h \leq 1m$ 时，可能还要加上地形造成的变异成分。对于这种复杂情况该如何建模？二是不同方向上的变异特征是否相似或相同，如果不同该如何建模？ 这种将出现在不同距离 $h$ 和不同方向上同时起作用的变异性组合模型结构，被称为套合结构。相关建模方法，参见地统计学文献，在此不再赘述。

3.5 变异函数的拟合

变异函数的拟合与具体选择的理论半方差函数有关，但基本思想都是最优曲线拟合。其拟合过程包括三个步骤：确定曲线类型、参数最优估计、最优曲线确定。

确定曲线类型：即 第 3.4 节 所述的经验数据方法，主体依靠专家的判断和假设。

参数最优估计：确定目标函数（通常是最小方差），然后采用最小二乘法做最优化估计。某些过于复杂的模型可以借助泰勒展开转换为线性模型后做近似优化。具体优化算法不再赘述。

最优曲线确定：即在多个曲线模型中选择最适合的那个。可以采用机器学习理论中的模型比较与选择，此外还应当包括对单个模型的检验。

4 克里金法简介

克里金法（Kriging Method）是以上述变异函数理论及其结构分析为基础，在有限区域内对区域化变量进行线性无偏最优估计（Best Linear Unbiased Prediction，BLUP）的一种方法，在地统计学中也被称为空间最优无偏估计器（Spatial BLUP）。

（1）克里金法的内涵

克里金法(Kriging）实质上是利用区域化变量的原始观测数据和变异函数的结构特点， 对未采样点的区域化变量取值进行线性无偏最优估计的一种方法。从数学角度讲，就是一种对空间分布的数据求线性最优无偏内插估计量的一种方法。更具体地讲，它是根据待估样点（或待估块段）有限邻域内若干已测定的样点数据，在认真考虑了样点的分布、数值、空间相互位置关系（体现为变异函数），以及样点与待估点之间的相互空间位置关系，对该待估点值进行的一种线性无偏最优估计。

克里金法与普通的估计方法不同，它最大限度地利用了空间取样所提供的各种信息。在估计未知点数值时，它不仅考虑了落在该点的数据，而且还考虑了邻近样点的数据；不仍考虑了待估点与邻近巳知样点的空间位置，而且还考虑了各邻近样点彼此之间的位置关系。除了上述的几何因素外，还利用了已有观测值空间分布的结构特征，使这种估计比其他传统估计方法更精确，更符合实际，并且避免系统误差的出现，同时该方法能够给出估计误差和精度（ Journel 等，1978；Ripley 1981）。

但需要注意的是：如果变异函数和相关分析结果表明区域化变量的空间相关性不存在，则克里金法不适用。

（2）克里金法的种类

地统计学的诞生与克里金的早期研究工作有着重要联系，以致国际上许多地统计学工作者（ Guarasio 等，1976; Verly 等，1983）将地统计学与克里金法同等对待。但其实更主要的原因是：空间预测（插值） 本身是地统计学最主要的研究任务之一。多年来地统计学家们围绕这一方面进行了大量研究工作，自 70 年代以来发展了多种空间局部估计方法。

克里金法顶层可以分为线性和非线性两类。其中：

线性克里金法：可以根据对变量均值的假设条件不同，分为简单克里金、普通克里金和泛克里金。

非线性克里金扩展：在线性克里金方法基础上，克里金方法朝几个方向进一步做了扩展。

• 一是比泛克里金提供更为宽松的非稳定假设前提，使克里金方法最优估值更佳;
• 二是朝多变量方向发展,如发展为各种协同克里金方法;
• 三是朝非线形克里金方法发展,如指示克里金法、析取克里金法等。

下面是对上述克里金方法的简单介绍:

• ➀ 线性克里金

  • 普通克里金法（Ordinary Kriging）： 是单个变量的局部线形最优无偏估计方法,也是最稳健最常用的一种方法。其假设均值为确定未知的常数。
  • 简单克里金法(Simple Kriging)：它假设空间过程的均值依赖于空间位置，并且是已知的，但在实际中均值一般很难得到，因此很少直接用于估计。但它可以用于其他形式的克里金法，如指示和析取克里金法，在这些方法中，数据进行了转换并且均值已知。
  • 泛克里金法(Universal Kriging)：把一个确定性的趋势模型加人到克里金估值中,将空间过程分解为趋势项和残差项两个部分，有其合理的一面。如果能够很容易地预测残差项的变异函数，那么该方法会得到广泛应用。

• ➁ 非线性克里金

  • 指示克里金法（Indictor Kriging）：将连续型随机变量转换为二值形式,是一种非线性、非参数的克里金预测方法。
  • 对数正态克里金法（Logistic Normal Kriging）：如果样本本身不服从正态分布，但其对数服从正态分布(（土壤的属性常常如此，Webster， 1985），则使用对数正态克里金法。
  • 析取克里金法（Disjunctive Kriging）：是一种有严格参数的非线性克里金方法。这种方法不但可以预测，还能够提供超过或不超过某一阈值的概率，因此从决策非常有用。

• ➂ 协同克里金法

  • 普通协同克里金法（Ordinary Cokriging）：将单个变量的普通克里金法扩展到两个或多个变量，这些变量间要存在一定的协同区域化关系。如果那些测试成本低、样品较多的变量与那些测试成本较高的、样品较少的变量在空间上具有一定的相关性，那么该方法就非常有用。可以利用较密采样得到的数据来提高样品较小数据的预测精度。

• ➃ 其他克里金

  • 概率克里金法(ProbabilityKriging)。由于指示克里金法并没有考虑一个值与阅值的 接近程度而只是它的位置,因此提出了概率克里金法。对每个值它利用rankorder作为辅助变量,利用协同克里金法来预测指示值。

（3）克里金法的主要估计量

对于任何一种估计，实际值和估计值之间总是存在偏差，不能要求估计值和实际值相等。问题是要求的估计方法在估计过程中要满足两点要求:

• 变异的平均数为零，即变异的数学期望等于零，此估计是一种无偏估计;
• 估计结果的方差最小，即估计是一种最优估计

克里金早期的工作中就是采用样品的加权平均值来求估计值的。也就是说，对于任意待估点或块段 $V$ 的实际值 $Z_V(\mathbf{s})$，其估计值 $Z_V^*(\mathbf{s})$ 是通过该待估点或待估块段影响范围内的 $n$ 个有效样品值 $\{ Z(\mathbf{s}_i), i=1,2,\ldots,n\}$ 的线性组合得到的，即

$$Z_V^*(\mathbf{s}) = \sum_{i=1}^{n} \lambda_i Z(\mathbf{s}_i)$$

上式中， $\lambda_i$ 为权重系数，是各已知样品 $Z_V(\mathbf{s})$ 在估计 $Z_V^*(\mathbf{s})$ 时影响大小的系数。由此可见，$Z_V^*(\mathbf{s})$ 的好坏主要取决于怎样计算或选择权重系数 $\lambda_i$，而这正是克里金法的主要内容。

根据式也可以看出，对于不同 $\mathbf{s}$ 显然有单独的一套权重系数，也就是说，克里金法需要对不同位置单独求解一套权重系数。这虽然提高了预测的计算复杂度，但其优点是能够反映不同位置的空间差异特性。

估计量 $Z_V^*(\mathbf{s})$ 被称为实际值 $Z_V(\mathbf{s})$ 的克里金估计量（Kriginge Etimator）。显然， $Z_V^*(\mathbf{s})$ 是 $Z_V(\mathbf{s})$ 的线性、无偏、最优估计量。

5 线性预测克里金法

5.1 普通克里金法

普通克里金法假定均值是未知的确定常数。主要有两种类型： 点状普通克里金 和 块段普通克里金。

5.1.1 点状普通克里金法（ Punctual Kriging ）

（1）预测式

假定 $Z(\mathbf{s})$ 是满足本征假设的一个随机过程, 该随机过程有 $n$ 个观测值 $z(\mathbf{s}_i), i=1,2, \cdots n$。要预测未采样点 $\mathbf{s}_0$ 处的值, 则线性预测值 $Z^*(\mathbf{s}_0)$ 可以表示如下：

$$Z^*\left(\mathbf{s}_0\right)=\sum_{i=1}^n \lambda_i Z\left(\mathbf{s}_i\right) \tag{5.1}$$

（2）模型参数求解

克里金法用无偏、最小方差预测作为准则，来优化权重值。它在满足以下两个条件下,可以实现线性无偏最优估计。

条件 1: 无偏性条件:

$$E\left[Z^*\left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right]=0 \tag{5.2}$$

条件 2: 最优条件:

$$\operatorname{Var}[Z^*(\mathbf{s}_0)-Z(\mathbf{s}_0)]=\min \tag{5.3}$$

利用 式 (5.1) 取代 式 (5.2) 的左边部分, 则有:

$$E\left[Z^*\left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right] = E\left[\sum_{i=1}^n \lambda_i Z\left(\mathbf{s}_i\right)-Z\left(\mathbf{s}_0\right)\right]=\sum_{i=1}^n \lambda_i E\left[Z\left(\mathbf{s}_i\right)\right]-E\left[Z\left(\mathbf{s}_0\right)\right] \tag{5.4}$$

结合本征假设 $E\left[Z\left(\mathbf{s}_i\right)\right]=E\left[Z\left(\mathbf{s}_0\right)\right]=m$ ，可以进一步表示为:

$$\sum_{i=1}^n \lambda_i E\left[Z\left(\mathbf{s}_i\right)\right]-E\left[Z\left(\mathbf{s}_0\right)\right]=\sum_{i=1}^n \lambda_i m-m=m\left(\sum_{i=1}^n \lambda_i-1\right)=0 \tag{5.5}$$

很显然要使 式 (5.5) 成立, 必须满足如下条件:

$$\sum_{i=1}^n \lambda_i=1 \tag{5.6}$$

在本征假设条件下, 式 (5.3) 左边的式子可表示为:

$$\begin{aligned} \operatorname{Var}\left[Z^*\left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right] & =E\left[\left\{Z^*\left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right\}^2\right] \\ & =E\left[\left\{\sum_{i=1}^n \lambda_i Z\left(\mathbf{s}_i\right)-Z\left(\mathbf{s}_0\right)\right\}^2\right] \\ & =2 \sum_{i=1}^n \lambda_i \gamma\left(\mathbf{s}_i, \mathbf{s}_0\right)-\sum_{i=1}^n \sum_{j=1}^n \lambda_i \lambda_j \gamma\left(\mathbf{s}_i, \mathbf{s}_j\right) \tag{5.7} \end{aligned}$$

式中 $\gamma\left(\mathbf{s}_i, \mathbf{s}_j\right)$ 是数据点 $\mathbf{s}_i$ 和 $\mathbf{s}_j$ 之间的半方差值; $\gamma\left(\mathbf{s}_i, \mathbf{s}_0\right)$ 是数据点 $\mathbf{s}_i$ 和预测点 $\mathbf{s}_0$ 之间的半方差值。

由本征平稳假设可知，任何克里金法的核心预测问题，就是计算 式（5.7） 的预测方差。为了得到最优估计，需要令其最小化以求得最优参数值。但根据 式（5.6），这是一个有约束的最优化问题。我们可以通过拉格朗日乘子法来实现优化求解，因此引入拉格朗日乘子 $\varphi$。

定义一个辅助函数 $f\left(\lambda_i, \varphi\right)$, 它的数学表达式为:

$$f\left(\mathbf{s}_i, \varphi\right)=\operatorname{Var}\left[Z^*\left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right]-2 \varphi\left(\sum_{i=1}^n \lambda_i-1\right) \tag{5.8}$$

分别对权重 $\lambda_i$ 和拉格朗日乘子 $\varphi$ 求一阶导数并使其等于 $0$，则对任意 $i=1,2,\ldots,n$，有：

$$\begin{aligned} & \frac{\partial f\left(\lambda_i, \varphi\right)}{\partial \lambda_i}=0 \\ & \frac{\partial f\left(\lambda_i, \varphi\right)}{\partial \varphi}=0 \tag{5.9} \end{aligned}$$

进而得到方程组:

$$\begin{aligned} & \sum_{i=1}^n \lambda_i \gamma\left(\mathbf{s}_i, \mathbf{s}_j\right)+\varphi=\gamma\left(\mathbf{s}_j, \mathbf{s}_0\right) \\ & \sum_{i=1}^n \lambda_i=1 \tag{5.10} \end{aligned}$$

上式中，变异函数 $\gamma(\cdot,\cdot)$ 通过经验变异函数拟合得出，可视为已知（这是克里金法与变异函数之间的主要结合点）。位置 $\mathbf{s}_0, \mathbf{s}_i$ 均为已知。因此，这是一个求 $n+1$ 个未知数的 $n+1$ 阶线性方程组，成熟的解方案，就此我们得到了克里金法所需的所有参数 $\lambda_i$ 和 $\varphi$。

（3）预测方差

将上述系数的解、已观测坐标 $\{\mathbf{s}_i: i =1,\ldots,N \}$ 、待估坐标点 $\mathbf{s}_0$ 直接代入 式（5.7）, 可以得到 $\mathbf{s}_0$ 处的预测方差为:

$$\sigma^2\left(\mathbf{s}_0\right)=\sum_{i=1}^n \lambda_i \gamma\left(\mathbf{s}_i, \mathbf{s}_0\right)+\varphi \tag{5.11}$$

如果预测点 $\mathbf{s}_0$ 恰好是其中的一个实测点 $\mathbf{s}_j$, 则有 $\lambda\left(\mathbf{s}_j\right)=1$ 而其他所有测点的权重等于 $0$，此时 $\sigma^2\left(\mathbf{s}_0\right)$ 最小。事实上将权重代人 式 (5.1) 时，有 $\sigma^2\left(\mathbf{s}_0\right)=0$。也就是说，预测值 $z\left(\mathbf{s}_0\right)$ 就是实测值 $z\left(\mathbf{s}_j\right)$ 。因此，点状克里金法被称为一种精确的插值方法。

（4）向量表示

上面的过程改写成向量形式更容易理解, 对点状克里金预测, 有:

$$A \cdot\left[\begin{array}{l} \boldsymbol{\lambda} \\ \varphi \end{array}\right]=B \tag{5.12}$$

这里

$$\begin{gathered} A=\left[\begin{array}{ccccc} \gamma\left(\mathbf{s}_1, \mathbf{s}_1\right) & \gamma\left(\mathbf{s}_1, \mathbf{s}_2\right) & \cdots & \gamma\left(\mathbf{s}_1, \mathbf{s}_n\right) & 1 \\ \gamma\left(\mathbf{s}_2, \mathbf{s}_1\right) & \gamma\left(\mathbf{s}_2, \mathbf{s}_2\right) & \cdots & \gamma\left(\mathbf{s}_2, \mathbf{s}_n\right) & 1 \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \gamma\left(\mathbf{s}_n, \mathbf{s}_1\right) & \gamma\left(\mathbf{s}_n, \mathbf{s}_2\right) & \cdots & \gamma\left(\mathbf{s}_n, \mathbf{s}_n\right) & 1 \\ 1 & 1 & \cdots & 1 & 0 \end{array}\right], \\ B=\left[\begin{array}{c} \gamma\left(\mathbf{s}_1, \mathbf{s}_0\right) \\ \gamma\left(\mathbf{s}_2, \mathbf{s}_0\right) \\ \cdots \\ \gamma\left(\mathbf{s}_n, \mathbf{s}_0\right) \\ 1 \end{array}\right],\left[\begin{array}{c} \boldsymbol{\lambda} \\ \varphi \end{array}\right]=\left[\begin{array}{c} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_n \\ \varphi \end{array}\right] \end{gathered} \tag{5.13}$$

矩阵 $A$ 中的 $\gamma\left(\mathbf{s}_1, \mathbf{s}_1\right), \ldots, \gamma\left(\mathbf{s}_n, \mathbf{s}_n\right)$ 是实测点之间的半方差值, 矩阵 $B$ 中的 $\gamma\left(\mathbf{s}_1\right.$, $\left.\mathbf{s}_0\right) \cdots \cdots \gamma\left(\mathbf{s}_n, \mathbf{s}_0\right)$ 为实测点 $\mathbf{s}_i$ 和待估点（或内插点） $\mathbf{s}_0$ 之间的半方差值, $\varphi$ 是拉格郎日乘子。矩阵 $A$ 是 可逆的, 因此有克里金参数:

$$\left[\begin{array}{l} \boldsymbol{\lambda} \\ \varphi \end{array}\right]=A^{-1} \cdot B \tag{5.14}$$

根据估计的参数，可得克里金预测方差为:

$$\sigma^2\left(\mathbf{s}_0\right) = B^T\left[\begin{array}{l} \boldsymbol{\lambda} \\ \varphi \end{array}\right] \tag{5.15}$$

5.1.2 块段普通克里金法

块段是空间中的一个子区间，根据其是一维、 二维或者三维的, 块段可以是一条线、一个平面或者一个立方体。令 $B$ 表示待估块段，其中心位置为 $\mathbf{s}$, 其均值为 $Z_B$。则应当有：

$$\begin{align*} Z(B) &= \frac{1}{B} \int_B Z(\mathbf{s}) d \mathbf{s}\\ E[Z(B)] &= \frac{1}{B} \int_B E[Z(\mathbf{s})] d \mathbf{s} = m \end{align*}$$

在待估块段 $B$ 的邻域内，有一组 $n$ 个已知样点 $\mathbf{s}_i, (i=1,2,\ldots,n)$，其观测值为 $Z(\mathbf{s}_i)$。 根据本征假设，其数学期望也为 $m$。

（1）预测公式

如果要预测的是一个块段 $B$ 的值，其实就是求解 $Z(B)$ 的线性估计量 $Z^*(B)$。这种情况下，块段的克里金预测值仍是数据的简单线性加权。

$$Z^*(B)=\sum_{i=1}^n \lambda_i Z\left(x_i\right) \tag{5.16}$$

（2）模型参数求解

与上一节相似，根据无偏性条件可得系数之和为 $1$。根据最优性条件，可得方差为:

$$\begin{aligned} \operatorname{var}\left[Z^*(B)\right] & =E\left[\left\{Z^*(B)-Z(B)\right\}^2\right] \\ & =2 \sum_{i=1}^n \lambda_i \bar{\gamma}\left(x_i, B\right)-\sum_{i=1}^n \sum_{j=1}^n \lambda_i \lambda_j \gamma\left(x_i, x_j\right)-\bar{\gamma}(B, B) \tag{5.17} \end{aligned}$$

式中: $\bar{\gamma}(x_i, B)$ 是第 $i$ 个样点与块段 $B$ 之间的平均半方差值, 可以用数学积分的形式来表示：

$$\bar{\gamma}\left(x_i, B\right)=\frac{1}{|B|} \int_B \gamma\left(x_i, x\right) d x \tag{5.18}$$

式中: $\gamma(x_i, x^{\prime})$ 表示样点 $x_i$ 和用来描述块段 $B$ 的点 $x$ 之间的半方差值。

$$\bar{\gamma}(B, B)=\frac{1}{|B|^2} \int_B \int_B \gamma\left(x_i, x^{\prime}\right) d x d x^{\prime} \tag{5.19}$$

式中: $\gamma\left(x_i, x\right)$ 是点 $x_i$ 和掠过块段 $B$ 的 $x^{\prime}$ 之间的半方差值。在点状克里金法中, $\bar{\gamma}(B$, B）变为 $\gamma\left(x_0, x_0\right)=0$ 。

由此，可得块段克里金方程组为:

$$\begin{aligned} & \sum_{i=1}^n \lambda_i \gamma\left(x_i, x_j\right)+\varphi(B)=\gamma\left(x_j, B\right) \\ & \sum_{i=1}^n \lambda_i=1 \tag{5.20} \end{aligned}$$

（3）预测方差

将结果代入 式 （5.17） ，可得预测方差为:

$$\sigma^2(B)=\sum_{i=1}^n \lambda_i \bar{\gamma}\left(x_i, B\right)+\varphi(B)-\bar{\gamma}(B, B) \tag{5.21}$$

（4）向量形式

改写为向量形式，其中 $B$ 为:

$$B=\left[\begin{array}{c} \gamma\left(x_1, B\right) \\ \gamma\left(x_2, B\right) \\ \cdots \\ \gamma\left(x_n, B\right) \\ 1 \end{array}\right] \tag{5.22}$$

预测方差可表示成:

$$\sigma^2(B)=B^T\left[\begin{array}{l} \lambda \\ \varphi \end{array}\right]-\bar{\gamma}(B, B) \tag{5.23}$$

5.2 简单克里金法

普通克里金法适用于均值未知的应用场景，但有时我们可以从先验知识中知道一个随机变量的均值, 此时可以利用这种先验知识来提高预测精度，这就是简单克里金法。这种克里金预测方法仍然是线性加和, 但将随机过程的均值 $\mu$ 包括了进去, 这种随机过程必须是二阶平稳的。简单克里金预测对那些只满足本征假设的随机过程并不是特别适合。简单克里金的预测公式为:

（1）预测公式

$$Z^*\left(\mathbf{s}_0\right)=\sum_{i=1}^n \lambda_i z\left(\mathbf{s}_i\right)+\left(1-\sum_{i=1}^n \lambda_i\right) \mu \tag{5.24}$$

其中 $\lambda_i$ 是权重，$\mu$ 为已知均值。

（2）模型参数求解

简单克里金法中不再有权重总和等于 $1$ 的条件。权重求解可以利用下面公式计算得来（推导过程参考 式 5.7）:

$$\sum_{i=1}^n \lambda_i \gamma\left(\mathbf{s}_i, \mathbf{s}_j\right)=\gamma\left(\mathbf{s}_0, \mathbf{s}_j\right), \quad j=1,2, \cdots, n \tag{5.25}$$

公式中不再有拉格朗日乘子, 成为一个 $n$ 阶线性方程组。

（3）预测方差

简单克里金的预测方差利用下面公式得到:

$$\sigma_{S K}^2\left(\mathbf{s}_0\right)=\sum_{i=1}^n \lambda_i \gamma\left(\mathbf{s}_i, \mathbf{s}_0\right) \tag{5.26}$$

（4）向量形式

预测方程组可以表示为向量形式为：

$$A \cdot \boldsymbol{\lambda}=B \tag{5.27}$$

其中：

$$\begin{gathered} A=\left[\begin{array}{cccc} \gamma\left(\mathbf{s}_1, \mathbf{s}_1\right) & \gamma\left(\mathbf{s}_1, \mathbf{s}_2\right) & \cdots & \gamma\left(\mathbf{s}_1, \mathbf{s}_n\right) \\ \gamma\left(\mathbf{s}_2, \mathbf{s}_1\right) & \gamma\left(\mathbf{s}_2, \mathbf{s}_2\right) & \cdots & \gamma\left(\mathbf{s}_2, \mathbf{s}_n\right) \\ \cdots & \cdots & \cdots & \cdots \\ \gamma\left(\mathbf{s}_n, \mathbf{s}_1\right) & \gamma\left(\mathbf{s}_n, \mathbf{s}_2\right) & \cdots & \gamma\left(\mathbf{s}_n, \mathbf{s}_n\right) \end{array}\right], \\ B=\left[\begin{array}{c} \gamma\left(\mathbf{s}_1, \mathbf{s}_0\right) \\ \gamma\left(\mathbf{s}_2, \mathbf{s}_0\right) \\ \cdots \\ \gamma\left(\mathbf{s}_n, \mathbf{s}_0\right) \end{array}\right], \boldsymbol{\lambda}=\left[\begin{array}{c} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_n \end{array}\right] \tag{5.28} \end{gathered}$$

矩阵 $A$ 中的 $\gamma\left(\mathbf{s}_1, \mathbf{s}_1\right) \cdots \cdots \gamma\left(\mathbf{s}_n, \mathbf{s}_n\right)$ 是实测点之间的半方差值, 矩阵 $B$ 中的 $\gamma\left(\mathbf{s}_1, \mathbf{s}_0\right)$ $\cdots \gamma\left(\mathbf{s}_n, \mathbf{s}_0\right)$ 为实测点 $\mathbf{s}_i$ 和待估点（内插点） $\mathbf{s}_0$ 之间的半方差值。矩阵 $A$ 是可逆的, 因此有:

$$\boldsymbol{\lambda}=A^{-1} \cdot B \tag{5.29}$$

同普通克里金法一样, 该方法可以推广到块段 $B$ 。通过将 式 (5.25) 中右边部分的 $\gamma\left(\mathbf{s}_0, \mathbf{s}_j\right)$ 用 $\bar{\gamma}\left(B, \mathbf{s}_j\right)$ 替换即可。

一般而言, 简单克里金法的预测方差稍小于由普通克里金预测的方差, 通常认为这是由于引进了数据均值 $\hat{\mu}$ 从而提高了预测的精度。

Wackernagel (1995) 表示, 如果要用到克里金预测的均值，即令 $\hat{\mu}=\hat{Z}(R)$, 则有

$$\sigma_{K K}^2\left(\mathbf{s}_0\right)=\sigma_{S K}^2\left(\mathbf{s}_0\right)+\left\{1-\sum_{i=1}^N \lambda_i^\kappa\right\} \psi(R) \tag{5.30}$$

换句话说, 即普通克里金预测的方差是简单克里金预测的方差与均值预测方差的加和。由于没有更多的信息, 采用这种方法很难额外获得什么。如果均值是利用大量的数据预测得到, 那么与 $\sigma_{S K}^2\left(\mathbf{s}_0\right)$ 相比, $\psi(R)$ 将会是很小的, 即假定简单克里金的权重和接近于 $1$ , 那么上式中右边第二部分就有可能是非常小的。

5.3 泛克里金法

（1）预测公式

泛克里金法假定均值未知, 但不是一个常数, 而是一个由多项式构成的漂移（趋势）:

$$m(\mathbf{s})=\sum_{l=1}^L a_l f_l(\mathbf{s}) \tag{5.31}$$

式中: $a_l$ 为未知系数; $l$ 为多项式阶数。此时, 泛克里金模型可以由随机性和确定性两个部分来构成:

$$Z(\mathbf{s})=m(\mathbf{s})+\boldsymbol{\varepsilon}(\mathbf{s}) \tag{5.32}$$

公式中 $\varepsilon(\mathbf{s})$ 为随机部分, 是一个符合 $[0,1]$ 分布的随机函数, 又叫做残差 (residual)。 确定性部分 $m(\mathbf{s})$ 叫做漂移 (drift)。

（2）模型参数求解

作为一种线性预测克里金方法, 泛克里金的估值结果也是寻求无偏条件下的最优解。但是与之前的简单克里金和普通克里金方法不同，在泛克里金法中，由于漂移的存在，在估计 $Z(\mathbf{s})$ 的估计值 $Z_{uk}^*(\mathbf{s})$ 时，需要先要估计漂移 $m(\mathbf{s})$ 的估计值 $m^*(\mathbf{s})$。也就是说，泛克里金法本质上是对漂移模型参数和残差模型参数同时求解的扩展克里金方法。

首先，泛克里金的预测方差可以分为两个部分（参阅相关文档）：

$$E[Z_{uk}^*(\mathbf{s}_0)-Z(\mathbf{s}_0)]^2 = \operatorname{Var}[Z_{uk}^*(\mathbf{s}_0)-Z(\mathbf{s}_0)] + (E[Z_{uk}^*(\mathbf{s}_0)-Z(\mathbf{s}_0)])^2 \tag{5.33}$$

根据无偏性原则，利用 式（5.31） 的均值表达式，应当有如下期望为 $0$：

$$\begin{aligned} E\left[Z_{uk}^*\left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right] & =\sum_{i=1}^n \lambda_i \sum_{l=1}^L a_l f_l\left(\mathbf{s}_i\right)-\sum_{l=1}^L a_l f_l\left(\mathbf{s}_0\right) \\ & =\sum_{i=1}^L a_l\left(\sum_{i=1}^n \lambda_i f_l\left(\mathbf{s}_i\right)-\sum_{i=1}^L f_l\left(\mathbf{s}_0\right)\right) \tag{5.34} \end{aligned}$$

式中 $a_l$ 为未知系数，要满足此条件，需有:

$$\sum_{i=1}^n \lambda_i f_l\left(\mathbf{s}_i\right)-\sum_{i=1}^L f_l\left(\mathbf{s}_0\right)=0 \text{ 或 } \sum_{i=1}^n \lambda_i f_l\left(\mathbf{s}_i\right)=\sum_{i=1}^L f_l\left(\mathbf{s}_0\right) \tag{5.35}$$

式中 $l=1, \cdots, L$。

可以看出， 式（5.35） 构成了 $L$ 个方程。因此, Matheron (1969) 称之为泛条件, 依此推导的克里金方法称为泛克里金法。

根据方差最小原则，式（5.33） 的均方误差 $E[Z_{uk}^*(\mathbf{s}_0)-Z(\mathbf{s}_0)]^2$ 应当最小化。依据无偏性假设，该式右边的第二项 $(E[Z_i^*(\mathbf{s}_0)-Z(\mathbf{s}_0)])^2$ 趋于 $0$。均方误差可以改写为:

$$\begin{align*} E\left[\left(Z_{uk} \left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right)^2\right]= & \operatorname{Var}\left[Z_{uk}^*\left(\mathbf{s}_0\right)-Z\left(\mathbf{s}_0\right)\right] \\ = & \sum_{i=1}^n \sum_{j=1}^n \lambda_i \lambda_j C_p\left(\mathbf{s}_i, \mathbf{s}_j\right)-2 \sum_{i=1}^n \lambda_i C_p\left(\mathbf{s}_i, \mathbf{s}_0\right) \\ & +C_p\left(\mathbf{s}_0, \mathbf{s}_0\right) \tag{5.36} \end{align*}$$

其中，$C_p\left(\mathbf{s}_i, \mathbf{s}_j\right)$ 表示两个区域化变量 $Z\left(\mathbf{s}_i\right)$ 和 $Z\left(\mathbf{s}_j\right)$ 之间的协方差函数, 即:

$$C_p(\mathbf{s}_i, \mathbf{s}_j) = E[Z(\mathbf{s}_i)Z(\mathbf{s}_j)] - E[Z(\mathbf{s}_i)] E[Z(\mathbf{s}_j)] \tag{5.37}$$

因此，问题最终归结为在 式（5.35） 约束下的 式（5.36） 最优化问题。通过拉格朗日乘子最优化方法, 可以得到如下泛克里金权重计算的方程组:

$$\begin{cases}\sum_{j=1}^n \lambda_j C_p\left(\mathbf{s}_i, \mathbf{s}_j\right)+\sum_{l=1}^L \varphi_k f_l\left(\mathbf{s}_i\right)-C_p\left(\mathbf{s}_i, \mathbf{s}_0\right)=0 & i=1,2, \cdots, n \\ \sum_{i=1}^n \lambda_i f_l\left(\mathbf{s}_i\right)-f_l\left(\mathbf{s}_0\right)=0 & l=0,1, \cdots, L\end{cases} \tag{5.38}$$

式（5.38） 中方程的数量受漂移的多项式假设约束。例如：采用一次多项式时（线性漂移） $L=3$，见 式（5.39）。采用二次多项式（二次漂移）时 $L=6$，见 式（5.40）。

该方程组共有 $n+L$ 个方程组成，需要求解 $n$ 个 $\lambda$ 参数，和 $L$ 个拉格朗日乘子参数 $\varphi$。

$$f_0=1, \quad f_1=\mathbf{s}_1 \quad f_2=\mathbf{s}_2 \tag{5.39}$$

$$f_3=\mathbf{s}_1^2, \quad f_4=\mathbf{s}_1 \mathbf{s}_2, \quad f_5=\mathbf{s}_2^2 \tag{5.40}$$

（3）估计方差

将系数解代入预测方程，可得估计方差为：

$$\sigma^2_{uk} = C_p(\mathbf{s},\mathbf{s})a - \sum_{i=1}^n \lambda_n C_p(\mathbf{s}_0,\mathbf{s}_i) + \sum_{l=1}^L \varphi_l f_l(\mathbf{s})$$

（4）向量形式

泛克里金方程式, 如同普通克里金法一样 [见公式 (5.12)], 是一系列的线性等式, 可 以用下列矩阵表示:

$$A \cdot\left[\begin{array}{l} \lambda \\ \varphi \end{array}\right]=B \tag{5.41}$$

矩阵 $A$ 和矢量 $\lambda$ 和 $B$ 是数据点和目标点的空间位置的函数, 可以用下式表示:

$$A=\left[\begin{array}{cccccc} C_p\left(\mathbf{s}_1, \mathbf{s}_1\right) & \cdots & C_p\left(\mathbf{s}_1, \mathbf{s}_n\right) & f_1\left(\mathbf{s}_1\right) & \cdots & f_L\left(\mathbf{s}_1\right) \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ C_p\left(\mathbf{s}_n, \mathbf{s}_1\right) & \cdots & C_p\left(\mathbf{s}_n, \mathbf{s}_n\right) & f_1\left(\mathbf{s}_n\right) & \cdots & f_L\left(\mathbf{s}_n\right) \\ f_1\left(\mathbf{s}_1\right) & \cdots & f_1\left(\mathbf{s}_n\right) & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ f_L\left(\mathbf{s}_1\right) & \cdots & f_L\left(\mathbf{s}_n\right) & 0 & \cdots & 0 \end{array}\right] \tag{5.42}$$

$${\left[\begin{array}{c} \lambda \\ \varphi \end{array}\right]=\left[\begin{array}{c} \lambda_1 \\ \vdots \\ \lambda_n \\ \varphi_0 \\ \vdots \\ \varphi_L \end{array}\right], B=\left[\begin{array}{c} C_p\left(\mathbf{s}_1, \mathbf{s}_0\right) \\ \vdots \\ C_p\left(\mathbf{s}_n, \mathbf{s}_0\right) \\ f_1\left(\mathbf{s}_0\right) \\ \vdots \\ f_L\left(\mathbf{s}_0\right) \end{array}\right]} \tag{5.43}$$

同普通克里金一样, $A$ 是可逆的, 且权重和朗格拉日乘子可以得到:

$$\left[\begin{array}{l} \lambda \\ \varphi \end{array}\right]=A^{-1} \cdot B \tag{5.44}$$

该权重被揷人到公式 (5.36) 中, 得到泛克里金方差为:

$$\sigma_{uk}^2=C_p\left(\mathbf{s}_0, \mathbf{s}_0\right)-\left[\begin{array}{l} \lambda \\ \varphi \end{array}\right]^T B \tag{5.45}$$

6 非线性预测克里金法

此处待遇到实际问题时再处理

7 协同克里金法

7.1 协同区域化特性与相关函数

自然现象中某一性质 (或变量) 的空间分布常常与其他性质（或变量）密切相关， 因为它们受同样的区域化现象或空间过程影响， 各个变量不但存在着自相关性， 而且与其他变量之间存在着交叉相关性， 这种性质被称为 协同区域化特性 (Co-regionalization)。

协同克里金法可以利用 “同一变量在不同时空” 或 “不同变量在同一时空” 上的协同区域化性质，“用易于测得的变量来对那些难以测得的属性或变量进行估计”，或者 “用样品多的变量対样品少的变量进行估计”。

协同克里金法将其中的一个变量作为主变量，其他的作为辅助变量。“主变量的自相关性” 以及 “主变量与其他变量之间的交叉相关性” 可以更好地进行预测。这种利用其他变量来帮助进行预测的想法非常诱人，但也有一定代价。

协同克里金要比普通克里金进行更多的预测，包括预测每个变量的自相关性和变量之间的交叉相关性。因此理论上，协同克里金法的预测精度总是比普通克里金的预测精度要高，因为它利用了变量间的交叉相关性来帮助进行预测。

不过，每次对未知自相关参数进行预测，又会引进更多的变异。

7.1.1 交叉协方差函数

如果两个区域化变量 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s})$ 是二阶平稳的，且均值分别为 $\mu_1$ 和 $\mu_2$，那么这两者之间的协方差函数分别为:

$$\begin{aligned} & C_{11}(h)=E\left[\left\{Z_1(\mathbf{s})-\mu_1\right\}\left\{Z_1(\mathbf{s}+h)-\mu_1\right\}\right] \\ & C_{22}(h)=E\left[\left\{Z_2(\mathbf{s})-\mu_2\right\}\left\{Z_2(\mathbf{s}+h)-\mu_2\right\}\right] \tag{7.1} \end{aligned}$$

两者之间的交叉协方差函数为:

$$C_{12}(h)=E\left[\left\{Z_1(\mathbf{s})-\mu_1\right\}\left\{Z_2(\mathbf{s}+h)-\mu_2\right\}\right] \tag{7.2}$$

交叉协方差函数公式还隐含另一特性，即不对称性。也就是说：

$$E\left[\left\{Z_1(\mathbf{s})-\mu_1\right\}\left\{Z_2(\mathbf{s}+h)-\mu_2\right\}\right] \neq E\left[\left\{Z_2(\mathbf{s})-\mu_2\right\}\left\{Z_1(\mathbf{s}+h)-\mu_1\right\}\right] \tag{7.3}$$

换言之， 在某一方向上 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s})$ 之间的交叉协方差不同于其反方向上的交叉协方差，函数不对称:

$$C_{12}(h) \neq C_{12}(-h)$$

不对称性非常普遍。例如同一剖面上的表土与底土存在着不对称关系。不过，除非不对称性非常明显或者可以用物理公式合理地解释，否则在预测的大多数场景中，可以认为其是由采样所造成的。

如果两个区域化变量 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s})$ 的实测值分别为 $z_1(\mathbf{s}_i)$ 和 $z_2(\mathbf{s}_i)$， 那么它们之间的经验交叉协方差函数计算如下:

$$C_{12}^{*}(h)=\frac{1}{2 N(h)} \sum_{i=1}^{N(h)} \{z_1(\mathbf{s}_i) z_2(\mathbf{s}_i + h) \} - m_{z_{1_{-h}}} \cdot m_{z_{2_{+h}}} \tag{7.4}$$

这里 $C_{12}(h)$ 是两个变量的经验交叉协方差值; $N(h)$ 是具有相同间距 $h$ 的变量对（ $Z_1(\mathbf{s}), Z_2(\mathbf{s})$ ）的数目; $m_{z_{1_{-h}}}$ 和 $m_{z_{2_{+h}}}$ 为相同间距 $h$ 内 $z_1\left(\mathbf{s}_i\right)$ 和 $z_2\left(\mathbf{s}_i\right)$ 的平均值， 分别计算如 下:

$$\begin{align*} m_{z_{1_{-h}}}=\frac{1}{2 N(h)} \sum_{i=1}^{N(h)} z_1(\mathbf{s}_i) \\ m_{z_{2_{+h}}}=\frac{1}{2 N(h)} \sum_{i=1}^{N(h)} z_2\left(\mathbf{s}_i+h\right) \end{align*}$$

7.1.2 交叉变异函数 (cross variogram)

假设有两个区域化变量 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s})$ 都遵从本征假设。对于变量 $Z_1(\mathbf{s})$ 有:

$$E\left[Z_1(\mathbf{s})-Z_1(\mathbf{s}+h)\right]=0 \tag{7.5}$$

该变量有自身的变异函数，即自变异函数:

$$\gamma_{11}(h)=\frac{1}{2} E\left[\left\{Z_1(\mathbf{s})-Z_1(\mathbf{s}+h)\right\}^2\right] \tag{7.6}$$

同样对区域化变量 $Z_2(\mathbf{s})$， 它增量的数学期望对任意 $\mathbf{s}$ 和 $h$ 存在且等于零， 同时其增量的方差函数对任意 $\mathbf{s}$ 和 $h$ 存在且平稳，也有一个自变异函数:

$$\begin{gather*} E\left[Z_2(\mathbf{s})-Z_2(\mathbf{s}+h)\right]=0 \tag{7.7}\\ \gamma_{22}(h)=\frac{1}{2} E\left[\left\{Z_2(\mathbf{s})-Z_2(\mathbf{s}+h)\right\}^2\right] \tag{7.8} \end{gather*}$$

则两个变量的交叉相关性可以用交叉半方差函数 $\gamma_{12}(h)$ 来表示:

$$\gamma_{12}(h)=\frac{1}{2} E\left[\left\{Z_1(\mathbf{s})-Z_1(\mathbf{s}+h)\right\}\left\{Z_2(\mathbf{s})-Z_2(\mathbf{s}+h)\right\}\right] \tag{7.9}$$

两个区域化变量 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s})$ 的实测值分别为 $z_1\left(\mathbf{s}_i\right)$ 和 $z_2\left(\mathbf{s}_i\right)$，它们之间的经验交叉半方差函数可以采用下式计算:

$$\gamma_{12}^*(h)=\frac{1}{2 N(h)} \sum_{i=1}^{N(h)}\left\{z_1\left(\mathbf{s}_i\right)-z_1\left(\mathbf{s}_i+h\right)\right\}\left\{z_2\left(\mathbf{s}_i\right)-z_2\left(\mathbf{s}_i+h\right)\right\} \tag{7.10}$$

其中 $\gamma_{12}^*(h)$ 是两个变量的经验交叉半方差值; $N(h)$ 是具有相同间距 $h$ 的变量对（ $Z_1(\mathbf{s})，Z_2(\mathbf{s})$ ）的数目。如果两个变量正相关，那么变量 $Z_1(\mathbf{s})$ 从 $\mathbf{s}_i$ 到 $\mathbf{s}_i+h$ 的增加 (或减少) 会相应引起 $Z_2(\mathbf{s})$ 的增加（或减少)， 交叉半方差值就是正值。

交叉变异函数和交叉协方差函数（如果存在的话）有如下关系:

$$\gamma_{12}(h)=C_{12}(0)-\frac{1}{2}\left\{C_{12}(h)+C_{12}(-h)\right\} \tag{7.11}$$

然而这种转化并没有保留所有的信息，我们通过将交叉协方差函数分为偶数部分和奇数部分来说明：

$$C_{12}(h)=\frac{1}{2}\left\{C_{12}(+h)+C_{12}(-h)\right\}+\frac{1}{2}\left\{C_{12}(+h)-C_{12}(-h)\right\} \tag{7.12}$$

从上式可以看出，奇数部分也就是公式右边的第二部分在 式 (7.11) 中并没有出现，因此交叉协方差函数是不对称的函数，而交叉变异函数是对称的偶函数，对任意 $h$ 有:

$$\gamma_{12}(h)=\gamma_{12}(-h) \text { 或 } \gamma_{12}(h)=\gamma_{21}(h)$$

由于交叉变异函数无法表达不对称性，因此不能用于那些具有非常明显不对称性（或不对称性非常重要）的现象或领域中。

7.1.3 交叉相关函数 (cross correlogram)

区域化变量 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s})$ 之间的空间交叉相关函数为 $\rho_{12}(h)$， 公式为:

$$\rho_{12}(h)=\frac{C_{12}(h)}{\sqrt{C_{11}(0) C_{22}(0)}}$$

该公式是变量 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s}+h)$ 在空间域上 Pearson 相关系数的扩展形式。当 $h=0$ 时就是 Pearson 相关系数。但需要注意， $\rho_{12}(0)=0$ 在通常意义下指没有线性相关性，但并不表示在大于零的滞后距离上没有相关性。

交叉相关函数与交叉协方差函数一样，也具有不对称性:

$$\rho_{12}(h) \neq \rho_{12}(-h)$$

两变量之间的经验交叉相关函数计算如下:

$$\rho_{12}^*(h)=\frac{C_{12}^*(h)}{\sigma_{z_{1_{-h}}} \cdot \sigma_{z_{2_{+h}}}}$$

式中: $\rho_{12}^*(h)$ 是两个变量的经验交叉协方差值; $C_{12}^*(h)$ 为经验交叉协方差函数， 可以用 式 (7.4) 计算获得; $\sigma_{z_{1_{-h}}}$ 和 $\sigma_{z_{2_{+h}}}$ 为相同间距 $h$ 内 $z_1\left(\mathbf{s}_i\right)$ 和 $z_2\left(\mathbf{s}_i\right)$ 的标准差值， 分别计算如下:

$$\begin{aligned} & s_{\mathbf{s}_1-h}^2=\frac{1}{N(h)} \sum_{i=1}^{N(h)} z_1^2\left(\mathbf{s}_i\right)-m_{z_1-h}^2， m_{\mathbf{s}_1-h}=2 \frac{1}{N(h)} \sum_{i=1}^{N(h)} z_1\left(\mathbf{s}_i\right) \\ & s_{z_{2+h}}^2=\frac{1}{N(h)} \sum_{i=1}^{N(h)} z_2^2\left(\mathbf{s}_i\right)-m_{z_{++h}}^2， m_{z_{z_{+h}}}=\frac{1}{N(h)} \sum_{i=1}^{N(h)} z_2\left(\mathbf{s}_i\right) \end{aligned}$$

其中 $N(h)$ 是具有相同间距 $h$ 的变量对（ $Z_1(\mathbf{s}),Z_2(\mathbf{s})$ ）的数目。

另外一个表达两个变量之间空间关系的手段是 协离散系数 (codispersion coefficient)， 对任意滞后距离 $h$ 有:

$$v_{12}(h)=\frac{\gamma_{12}(h)}{\sqrt{\gamma_{11}(h) \gamma_{22}(h)}}$$

该系数可以表达两个变量 $Z_1(\mathbf{s})$ 和 $Z_2(\mathbf{s})$ 的 空间离差（spatial differences） 之间的相关性，它是对称的且优于交叉相关性函数 $\rho_{12}(h)$ 。对二阶平稳性而言，当滞后距离 $h$ 的绝对值趋于无穷大时，$v_{12}(h)$ 接近于 $\rho_{12}(k)$。

7.2 普通协同克里金法 (Ordinary Cokriging)

普通协同克里金预测是权重 $\lambda_a^i$ 的一种线性结合。假设待预测变量为 $Z_u$，通常该变量比其他变量的样品少；有 $V$ 个协变量 $Z_j， j=1，2， \cdots， V$。

（1）预测公式

利用普通协同克立格法，有下列预测公式:

$$Z_u^*(B)=\sum_{l=1}^V \sum_{i=1}^{n_{l}} \lambda_{il} \left(z_l\left(\mathbf{s}_i\right)\right) \tag{7.16}$$

式中 $Z_u^*(B)$ 表示待估点或块段处的预测值; 下标 $l$ 表示第 $l$ 个辅助变量， 共有 $V$ 个。下标 $i$ 表示变量 $l$ 实测的样点数目，共有 $n_l$ 个是实测数据; $\lambda$ 是分配给变量 $l$ 的权重。

（2）模型求解

根据无偏性原则，需要满足:

$$\sum_{i=1}^{n_l} \lambda_i=\left\{\begin{array}{l} 1， l=u \\ 0， l \neq u \end{array}\right. \tag{7.17}$$

以两个变量的情况为例。

若只有两个变量 $Z_u$ 和 $Z_v$，其中待预测的变量 $Z_u$ 为主变量，则普通协同克里金的预测 式 (7.16) 可以简化为:

$$Z_u^*(B)=\sum_{i=1}^N \lambda_{u i} z_u\left(\mathbf{s}_i\right)+\sum_{j=1}^p \lambda_{i j} z_v\left(\mathbf{s}_j\right) \tag{7.18}$$

式中: $Z_u^*(B)$ 是待估点或块段处的预测值; $Z_u\left(\mathbf{s}_i\right)$ 和 $Z_v\left(\mathbf{s}_j\right)$ 分别是主变量 $Z_u$ 和辅助变量 $Z_v$ 的实测值; $\lambda_i$ 和 $\lambda_j$ 分别是分配给主变量 $Z_u$ 和辅助变量 $Z_v$ 的实测值的权重， 且 $\sum \lambda_{ui}=1$， $\sum \lambda_{vj}=0$ ; $n$ 和 $p$ 是参与估值的主变量 $Z_u$ 和辅助变量 $Z_v$ 的实测值数目。

根据最优无偏的假设条件，建立普通协同克里金的预测方程组：

$$\begin{aligned} & \sum_{l=1}^v \sum_{i=1}^{n_i} \lambda_{il} \gamma_{lv}\left(\mathbf{s}_i， \mathbf{s}_j\right)+\varphi_v=\bar{\gamma}_{uv}\left(\mathbf{s}_j， B\right) \\ & \sum_{i=1}^{n_l} \lambda_{il}=\left\{\begin{array}{l} 1， l=u \\ 0， l \neq u \end{array}\right. \end{aligned}$$

式中 $\gamma_{lv}(\mathbf{s}_i， \mathbf{s}_j)$ 是变量 $l$ 和 $v$ 在点 $i， j$ 处的交叉半方差值; $\bar{\gamma}_w\left(\mathbf{s}_i， B\right)$ 是点 $i$ 和待预测块段处 $B$ 的交叉半方差值平均值; $\varphi_v$ 是第 $v$ 个变量的拉格朗日乘子。

根据 式 (7.19) 可以计算出权重 $\lambda$， 将其代人普通协同克里金预测 式 (7.16) 中，可以得到预测值。

（3）估计方差

普通协同克里金的预测方差为：

$$\sigma_u^2(B)=\sum_{l=1}^v \sum_{j=1}^{n_l} \lambda_{jl} \bar{\gamma}_{vl}\left(\mathbf{s}_j， B\right)+\varphi_u = \bar{\gamma}_{uv}(B， B)$$

（4）向量表达形式

上述计算可以用矩阵的形式来表达，为了简单起见只考虑两个变量 $u$ 和 $v$，用 $\tau_{uv}$ 表示样点之间的半方差值的矩阵 (包括交叉半方差值，当 $u \neq v$ 时)。假设变量 $u$ 有 $n_u 个点进行了实测，变量 $v$ 在 $n_v$ 个点进行了实测，那么这个矩阵的阶数是 $n_u \times n_v$ :

$$\tau_{u v}=\left[\begin{array}{cccc} \gamma_{uv}\left(\mathbf{s}_1， \mathbf{s}_1\right) & \gamma_{uv}\left(\mathbf{s}_1， \mathbf{s}_2\right) & \cdots & \gamma_{uv}\left(\mathbf{s}_1， \mathbf{s}_{n_v}\right) \\ \gamma_{uv}\left(\mathbf{s}_2， \mathbf{s}_1\right) & \gamma_{uv}\left(\mathbf{s}_2， \mathbf{s}_2\right) & \cdots & \gamma_{uv}\left(\mathbf{s}_2， \mathbf{s}_{n_v}\right) \\ \vdots & \vdots & \cdots & \vdots \\ \gamma_{uv}\left(\mathbf{s}_{n_u}， \mathbf{s}_1\right) & \gamma_{uv}\left(\mathbf{s}_{n_u}， \mathbf{s}_2\right) & \cdots & \gamma_{uv}\left(\mathbf{s}_{n_u}， \mathbf{s}_{n_v}\right) \end{array}\right]$$

变量 $u$ 的自相关向量和变量 $u$ 和 $v$ 之间的交叉相关向量分别用 $b_{uv}$ 和 $b_{vu}$ 表示:

$$b_{vu}=\left[\begin{array}{c} \bar{\gamma}_{vu}\left(\mathbf{s}_1， B\right) \\ \bar{\gamma}_{vu}\left(\mathbf{s}_2， B\right) \\ \vdots \\ \bar{\gamma}_{vu}\left(\mathbf{s}_{n_\alpha}， B\right) \end{array}\right]， \quad b_{uv}=\left[\begin{array}{c} \bar{\gamma}_{uv}\left(\mathbf{s}_1， B\right) \\ \bar{\gamma}_{uv}\left(\mathbf{s}_2， B\right) \\ \vdots \\ \bar{\gamma}_{uv}\left(\mathbf{s}_{n_v}， B\right) \end{array}\right]$$

矩阵等式为:

现用 $G$ 表示等式左边的第一部分， 用 $\lambda$ 表示等式左边的第二部分， 用 $b$ 表示右边的部分， 那么等式可以写为:

$$\lambda=G^{-1} b$$

协同克里金的预测方差为:

$$\hat{\sigma}_u^2(B)=b^T \lambda-\bar{\gamma}_u(B， B)$$

这里块段 $B$ 可以是任何形状和大小，也可以是一个点 $\mathbf{s}_0$ 。这时 $\bar{\gamma}_{w v}\left(\mathbf{s}_{n_v}， B\right)$ 就变为 $\bar{\gamma}_w\left(\mathbf{s}_{n_v}， \mathbf{s}_0\right)$， 且 $\bar{\gamma}_{w t}(B， B)=0$ 。上面的 $b_{w t}$ 和 $b_{w \omega}$ 向量可表示成:

$$b_w=\left[\begin{array}{c} \gamma_\omega\left(\mathbf{s}_1， \mathbf{s}_0\right) \\ \gamma_\omega\left(\mathbf{s}_2， \mathbf{s}_0\right) \\ \vdots \\ \gamma_\omega\left(\mathbf{s}_{n_w}， \mathbf{s}_0\right) \end{array}\right]， \quad b_w=\left[\begin{array}{c} \gamma_w\left(\mathbf{s}_1， \mathbf{s}_0\right) \\ \gamma_w\left(\mathbf{s}_2， \mathbf{s}_0\right) \\ \vdots \\ \gamma_w\left(\mathbf{s}_{n_v}， \mathbf{s}_0\right) \end{array}\right]$$

同时 式 (7.22) 的协同克里金的预测方差可表示为:

$$f_u^2\left(\mathbf{s}_n\right)=b^T \lambda$$

7.3 简单协同克里金法 (Simply Cokriging)

在给定的邻域内，当感兴趣的变量没有样品数据可以用时，普通协同克里金没有均值可用。当情况相反时，可以考虑采用简单协同克里金法。简单协同克里金法要用到变量的经验均值信息，它由感兴趣变量的均值与权重 $\lambda_z$ 的线性结合两个部分组成:

$$Z_u^*(B)=m_u+\sum_{l=1}^V \sum_{i=1}^{n_l} \lambda_{i l}\left(z_l\left(\mathbf{s}_i\right)-m_l\right)$$

式中: $Z_u^*(B)$ 表示所预测的块段（或点）处的预测值; $l$ 表示变量，$l=1，2，\ldots,V$；$i$ 表示变量 $l$ 的实测点序号，共 $n_1$ 个; $\lambda$ 是权重。

简单协同克里金预测系统的矩阵形式为:

$$\left[\begin{array}{ccccc} C_{11} & \cdots & C_{1 j} & \cdots & C_{1 V} \\ \vdots & \ddots & & & \vdots \\ C_{l 1} & & C_{l j} & \ddots & C_{l V} \\ \vdots & & \vdots & \ddots & \vdots \\ C_{V 1} & \cdots & C_{V j} & \cdots & C_{V V} \end{array}\right]\left[\begin{array}{c} \lambda_1 \\ \vdots \\ \lambda_l \\ \vdots \\ \lambda_V \\ \end{array}\right]=\left[\begin{array}{c} c_{1 B} \\ \vdots \\ c_{l B} \\ \vdots \\ c_{V B} \end{array}\right]$$

$$C_{l j}=C_{j l}^T$$

$C_{l j}$ 包含了一定对的变量间样品点之间的协方差， $c_{l B}$ 包括了所要预测的点（或块段）和实测变量 $l$ 的样品点之间的协方差，$\lambda_l$ 是分配给变量 $l$ 的权重。