Cressie 的最新空间统计论述

【阅读建议】 本文是 Cressie 在 2021 年新撰写的一篇综述类文章，其主要看点包括：（1）用统一的形式化框架实现了点参考数据、面元数据、点模式数据的建模；（2）对多变量空间统计建模的统一形式化；（3）大数据的空间离散化处理方法（此处尚未理解其优势所在，需要进一步阅读引用的论文）；

【摘 要】 空间统计是一个致力于与空间标签相关数据统计分析的研究领域。地理学家通常将 “位置信息” 与 “属性信息” 联系起来，并且定义了一个被称为 “空间分析” 的研究领域。许多操作空间数据的方法都是由算法驱动的，缺少与之相关的不确定性量化。如果空间分析是统计的（即结合了不确定性量化），则它属于空间统计的研究范畴。空间统计模型的主要特征是邻近的属性值比远处的属性值在统计上更相关，这也被称为地理学第一定律。

【原 文】 N. Cressie and M. T. Moores, “Spatial Statistics,” 2021, doi: 10.48550/ARXIV.2105.07216.

【参 考】

1 导言

空间统计提供了一个概率框架，用于回答数据中包含空间位置信息、且所提问题与空间位置信息相关的科学问题。概率论在空间统计中的主要作用是对不确定性建模，这既包含科学理论中的不确定性，也包括空间数据中的不确定性。

在空间统计中，科学理论中的不确定性可以通过 空间随机过程 以概率形式表示，一般性地可以写成：

$$\{Y(\mathbf{s}) : \mathbf{s} \in \mathscr{D} \} \tag{1}$$

其中 $Y(\mathbf{s})$ 是位置 $\mathbf{s}$ 处的随机属性值，$\mathscr{D}$ 是 $d$ 维空间的子集，$\mathbb{R}^d$ 是欧几里德空间，它索引了感兴趣的所有可能空间位置。 $\mathscr{D}$ 中包含一个集合 $D$ （可能是随机的），它索引了 $\mathscr{D}$ 中与科学研究相关的那些部分空间域。$D$ 可以具有不同的集合性质，主要取决于空间过程是地统计过程、格元过程还是点过程（即空间统计中最基础的三类空间数据类型）。

我们可以将由具有随机性的 $\{Y(\mathbf{s}) : \mathbf{s} \in D\}$ 和具有随机性的 $D$ 定义的联合概率模型，方便地表示为简写 $[Y, D]$，亦即 空间过程模型（Spatial Process Model）。根据贝叶斯规则，有：

$$[Y, D] = [Y \mid D][D] \tag{2}$$

对于一般性的随机变量 $A$ 和 $B$，两者的联合概率在本文中用 $[A, B]$ 表示；给定 $B$ 时 $A$ 的条件概率度量表示为 $[A \mid B]$；$B$ 的边缘概率测度则记为 $[B]$。在本文中，式 (2) 给出了 Cressie（1993 年）作出的空间统计模型的一般性形式化定义。

根据对 $D$ 的三种不同假设，式 (2) 涵盖了三个主要空间统计领域，导致了三种不同类型的空间随机过程 $[Y \mid D]$，具体细节见 第 2 节 的 “空间过程模型”。过去，人们习惯于根据空间数据类型 $\mathbf{Z}$ 对空间统计方法进行分类（例如 Ripley，1981；Upton 和 Fingleton，1985；Cressie，1993），而不是根据生成数据的底层空间过程 $Y$ 的类型来划分。

在本文中，我们将根据过程模型 式 (2) 对空间统计建模选择进行了分类；而数据 $\mathbf{Z}$ 的分布则被定义称为数据模型，定义为同时给定 式 (2) 中的 $Y$ 和 $D$ 时，数据 $\mathbf{z}$ 的条件概率：

$$[\mathbf{Z} \mid Y，D] \tag{3}$$

例如，空间数据 $\mathbf{Z}$ 是在给定空间位置 $\{\mathbf{s}_1, . . . , \mathbf{s}_n \} \subset D$ 处 $\mathbf{Y}$ 的不完美测量值，并且假设数据是条件独立的。也就是说，数据模型是：

$$[\mathbf{Z} \mid Y, D] = \prod^n_{ i=1} [Z(\mathbf{s}_i) \mid Y, D] \tag{4}$$

这里需要注意的是：尽管 式 (4) 来自于条件独立性假设，但边缘分布 $[\mathbf{Z} \mid D]$ 并不表现出独立性，$\mathbf{Z}$ 中的空间统计依赖性继承自 $[Y \mid D]$ 和 式 (4)，并且有：

$$[\mathbf{Z} \mid D] = \int [\mathbf{Z} \mid Y, D][Y \mid D]dY$$

另一个例子是 $D$ 存在随机性，即一个点过程（随机集的特例），此时数据 $\mathbf{Z} = \{N, \mathbf{s}_1, \ldots , \mathbf{s}_n\}$，其中 $N$ 是在有界区域 $\mathscr{D}$ 中点的随机计数，并且 D = \{\mathbf{s}_1,\lots , \mathbf{s}_n\} 是点的随机位置。如果存在与 $D$ 中的随机点相关联的测量（ 有时称为“标记” ）$\{Z(\mathbf{s}_1), \ldots, Z(\mathbf{s}_n)\}$，则这些测量应包含在 $\mathbf{Z}$ 中。也就是说，

$$\mathbf{Z} = \{N, (\mathbf{s}_1, Z(\mathbf{s}_1)) , . . . , (\mathbf{s}_n, Z(\mathbf{s}_n))\} \tag{5}$$

式 (2) 和 式 (3) 给出的空间统计描述，捕获了待解决科学问题中的（已知）不确定性，即空间过程模型中（ 式 (2) ）的科学不确定性和数据模型中（ 式 (3)） 的测量不确定性。 式 (2) 和 式 (3) 一起定义了一个适用于空间数据的分层统计模型。不过，这种条件概率分布 $[\mathbf{Z} \mid Y, D]$, $[Y \mid D]$ 和 $[D]$ 的分层形式，其实适用于所有应用统计。

式 (2) 和 式 (3) 隐含了一个条件，即与 过程模型 和 数据模型 相关的所有参数 $\theta$ 都是已知的。但其实参数恰恰是我们要估计的对象，也存在着不确定性。贝叶斯方法会在 $\theta$ 上放置一个概率分布，即令 $[ \theta ]$ 表示捕获参数不确定性的 参数模型（或先验模型）。通过使用显式符号，所研究科学问题中的所有不确定性，都可以通过联合概率形式表示为：

$$\begin{align*} [\mathbf{Z},Y,D, \theta ] &= [\mathbf{Z},Y,D \mid \theta ][ \theta ] \tag{6}\\ &= [\mathbf{Z} \mid Y, D, \theta ][Y \mid D, \theta ][D \mid \theta ][ \theta ] \tag{7} \end{align*}$$

非常出名的 贝叶斯分层模型 就是使用了式 (7) 给出的链式分解形式。不过，也存在一种 经验分层模型，它将 $\theta$ 的点估计 $\hat{\theta}$ 代入 式 (6) 左侧的第一个因子中，形成如下模型：

$$[\mathbf{Z},Y,D \mid \hat{ \theta} ] = [\mathbf{Z} \mid Y, D, \theta ][Y \mid D, \theta ][D \mid \hat{\theta} ] \tag{8}$$

从空间数据 $\mathbf{Z}$ 中找到 $\theta$ 的有效估计是空间统计中的一个重要问题。但在本文中，我们的重点是 $Y$ 的空间预测问题，因此在下文中将假设参数已知或已被估计得出。为方便起见，我们可以删除 式(8) 中的 $\theta$ ，则问题中的不确定性可以简化为以下联合概率表示：

$$[\mathbf{Z},Y,D] = [\mathbf{Z} \mid Y, D][Y \mid D][D] \tag{9}$$

通过贝叶斯法则，可以利用预测分布来推断未知的 $Y$ 和 $D$：

$$[Y, D \mid \mathbf{Z}] = \frac{[\mathbf{Z} \mid Y, D][Y \mid D][D]}{[\mathbf{Z}]} \tag{10}$$

其中 $[\mathbf{Z}]$ 是确保 式 (10) 右侧积分（或总和）为 $1$ 的归一化常数。如果空间索引集 $D$ 固定且已知，则可以从 式 (10) 和贝叶斯公式中删除 $D$ ，进而简化为：

$$[Y \mid \mathbf{Z}] = \frac{[\mathbf{Z} \mid Y][Y]}{[\mathbf{Z}]} \tag{11}$$

上式是 $Y$ 的预测分布（当 $D$ 固定且已知时）。空间统计中经常使用这种表达方式进行预测。例如，在 式 (2) 和 式 (3) 的高斯分布假设下，可以很容易发现简单克里格预测是 式 (11) 的预测均值（Cressie 和 Wikle，2011）。

本文结构如下：

• 第 2 节 介绍 “空间过程模型” ，包括常用空间过程模型和多变量模型。
• 第 3 节 介绍 “空间离散化” ，重点考虑 $D \subset \mathbb{R}^d$ 的离散化问题，因为在实际计算 式 (10) 或 式 (11) 的预测分布时，这是极其重要的因素之一。
• 第 4 节 介绍 “时空过程模型” ，讨论空间过程模型到时空过程模型的扩展。
• 第 5 节 为 “结论”，简要讨论了空间统计领域近期的重要研究课题，但由于篇幅有限，我们无法完整呈现。

2 空间过程模型

在本节中，我们讨论式 (10) 贝叶斯法则中的概率分布 $[Y \mid D]$、$[Y]$ 和 $[D]$ 及其在空间上下文中的表示。注意不要将过程模型与 $[\mathbf{Z} \mid D]$ 和 $[\mathbf{Z}]$ 这些空间数据的概率分布混淆。在许多空间统计文献中，当研究人员直接为 $[\mathbf{Z} \mid D]$ 建立模型时，这种混淆很常见。采用分层方法有助于我们避免这种混淆，我们可以从统计模型开始获取科学过程知识 $[Y \mid D]$ 和 $[D]$，然后通过数据模型 $[\mathbf{Z} \mid Y，D]$ 对测量误差和缺失数据进行建模。最后，式(10) 的贝叶斯法则允许我们通过预测分布 $[Y, D \mid \mathbf{Z}]$ 来推断 $Y$ 和 $D$。

我们提出了三种类型的空间过程模型，它们的区别是根据定义过程 $Y$ 的空间位置索引集 $D$ 进行的。

• 对于地统计过程，$D = D^G$，这是一个已知集，其位置连续变化且其面积（或体积）$> 0$。
• 对于格元过程，$D = D^L$，这是一个已知集，其位置离散变化并且位置的数量是可数的；请注意，$D^L$ 的面积等于零。
• 对于点过程，$D = D^P$，它是由 $\mathbb{R}^d$ 中的随机点组成的随机集合。

2.1 地统计过程

在本节中，我们假设空间位置 $D$ 由 $D^G$ 给出，其中 $D^G$ 已知。因此，可以从 式（10） 中的任何概率分布中删除 $D$，从而得到 式（11）。这使我们能够专注于 $Y$，并且为了体现空间索引，我们将 $Y$ 等价地写为 $\{Y(\mathbf{s}) : \mathbf{s} \in D^G\}$。地统计过程的一个性质是其空间域 $D^G$ 具有正面积，因此 $\mathbf{s}$ 是连续不可数的。

传统上，地统计过程被指定为二阶矩。从最一般性的定义开始，我们有

$$\begin{align*} \mu_Y(\mathbf{s}) &\equiv E (Y(\mathbf{s})) ; &\mathbf{s} \in D^G \tag{12}\\ C_Y (\mathbf{s,u}) &\equiv \operatorname{cov}(Y(\mathbf{s}),Y(\mathbf{u})); \quad &\mathbf{s,u} \in D^G \tag{13} \end{align*}$$

从 式（12） 和 式（13） 中，我们的目标是获得域内任意点 $\mathbf{s}_0$ 处 $Y(\mathbf{s}_0)$ 的最佳空间线性预测 $\hat{Y}(\mathbf{s}_0)$，该预测依赖于空间数据 $\mathbf{Z} \equiv (Z(\mathbf{s}_1), ..., Z(\mathbf{s}_n))^\prime$，这是一个由 $n$ 个已知空间位置 $D^{G∗} \equiv \{\mathbf{s}_1, \ldots , \mathbf{s}_n \} \subset D^G$ 索引的 $n$ 维数据向量。

在实际工作中，由于缺乏重复数据（来自于随机过程的单样本假设），完全估计出 式（12） 和 式（13） 中的参数 $\theta$ 可能会存在困难，因此 Matheron (1963) 做出了 本征平稳 的平稳性假设。也就是说，对于所有 $\mathbf{s，u} \in D^G$，假设

$$\begin{align*} E(Y(\mathbf{s})) &= \mu ^o_Y \tag{14}\\ \operatorname{var} (Y(\mathbf{s}) −Y(\mathbf{u})) &= 2 \gamma^o_Y (\mathbf{s} − \mathbf{u}) \tag{15} \end{align*}$$

其中 式（15） 等于 $C_Y (\mathbf{s,s}) + C_Y (\mathbf{u,u}) − 2C_Y (\mathbf{s,u})$。量 $2 \gamma^o_Y (·)$ 被称为变异函数， $\gamma^o_Y(·)$ 称为半变异函数（有时也称为半方差）。

如果 式（15） 中的假设被替换为

$$\operatorname{cov}(Y(\mathbf{s}),Y (\mathbf{u}) = C^o_Y (\mathbf{s} − \mathbf{u}), \text{ for all } \mathbf{s,u} \in D^G \tag{16}$$

那么 式（16） 和 式（14） 一起称为二阶平稳。 Matheron 选择 式（15） 是因为本征平稳可以在无需知道或估计 $\mu^o_Y$ 的情况下，就推导出 $Y(\mathbf{s}_0)$ 的空间最优线性预测（即克里格法）。这里，“最优” 指空间线性预测 $\hat{Y}(\mathbf{s}_0)$ 能够最小化均方预测误差 (MSPE) ，MSPE 的定义见下式：

$$E \left[ (\hat{Y}(\mathbf{s}_0) −Y(\mathbf{s}_0))^2 \right] , \text{ for any } \mathbf{s}_0 \in D^G \tag{17}$$

其中 $\hat{Y}(\mathbf{s}_0) \equiv \sum^n_{i=1} \lambda_i Z(\mathbf{s}_i)$。

式（17） 的最小化受无偏估计 $E(\hat{Y}(\mathbf{s}_0)) = E (Y (\mathbf{s}_0))$ 无偏约束，或等效地在 $\{\lambda_i\}$ 上受 $\sum^n_{i=1} \lambda_i = 1$ 约束，并且优化是关于系数 $\{\lambda _i : i = 1, ..., n\}$ 的 。通过选择最优的 $\{\lambda_i\}$，可以得到最优线性预测 $Y(\mathbf{s}_0)$ ，也被称为克里金预测。 Matheron 将这种空间预测方法称为普通克里金法，尽管它在其他领域被称为 BLUP（最优线性无偏预测）； Cressie (1990) 介绍了克里金法的历史，并表明它可以被表述为空间最优线性无偏预测。

图 1：2009 年 1 月澳大利亚温度的克里金预测图，叠加在数据的空间位置上。

图 2：图 1 中克里金预测对应的克里金标准差，式（18） 。

式（14） 中的常量均值假设可以扩展到线性回归 $E(Y(\mathbf{s})) \equiv \mathbf{X}(\mathbf{s})^\prime \boldsymbol{\beta}$，其中 $\mathbf{s} \in D^G$，回归系数 $\boldsymbol{\beta}$ 未知，协变量向量 $\mathbf{X}(\mathbf{s})$ 中包括截距项对应的元素 $1$。 在这种 $E(Y(\mathbf{s}))$ 假设下，普通克里金法 被推广为 通用克里金法，同样记为 $\hat{Y}(\mathbf{s}_0)$。图 1 显示了 2009 年 1 月澳大利亚温度的 通用克里金 预测，覆盖了整个澳洲大陆 $D^G$ ，提供数据 $\mathbf{Z}$ 的气象站位置 $D^{G*} = \{\mathbf{s}_1, \ldots, \mathbf{s}_n\}$ 被叠加在图中。读者可以在 Chiles 和 Delfiner（2012，第 3.4 节）中找到 $\hat{Y}(\mathbf{s}_0)$ 的公式。

式（17） 的 MSPE 被称为 克里金方差，其平方根称为 克里金标准差：

$$\boldsymbol{\sigma}_k(\mathbf{s}_0) \equiv \left( E (\hat{Y} (\mathbf{s}_0) −Y (\mathbf{s}_0))^2\right) ^{1/2} ，\text{ for any } \mathbf{s}_0 \in D^G \tag{18}$$

图 2 显示了与 图 1 中克里金预测变量相关的克里金标准差 $D^G$ 地图。可以看出，较小的 $\boldsymbol{\sigma}_k(\mathbf{s}_0)$ 对应于 $\mathbf{s}_0$ 附近较高密度的气象站。虽然 普通克里金法 和 通用克里金法 能够生成最佳线性预测，但还有一个更好的预测，被称为 最佳最优预测 (best optimal predictor, BOP)，它是在附加约束（如线性）下获得的所有最优预测中最好的那个。根据 式（10） 的贝叶斯规则，在没有任何约束的情况下最小化 MSPE 式（16） 得到的最优预测 $Y^∗(\mathbf{s}_0) \equiv E (Y (\mathbf{s}_0) \mid \mathbf{Z})$，应当是预测分布的均值。请注意，BOP $Y^∗(\mathbf{s}_0)$ 自身是无偏的，即 $E(Y^∗(\mathbf{s}_0)) = E(Y(\mathbf{s}_0))$，无需将其约束为无偏。

2.2 格元过程

在本节中，我们假设空间位置 $D$ 由 $D^L$ 给出，$D^L$ 是 $\mathbb{R}^d$ 的已知可数子集。这通常表示一组网格节点、像素或小区域以及与它们相关联的空间位置；我们将所有此类位置的可数集写为 $D^L \equiv \{\mathbf{s}_1, \mathbf{s}_2,\ldots \}$。其中每个 $\mathbf{s}_i$ 都有一组邻居与其相关联，表示为 $\mathscr{N}(\mathbf{s}_i) \subset D^L \ \mathbf{s}_i$，并且邻居的位置在空间上接近该点（注意，不能视自身为邻居）。格元过程中位置之间的空间统计依赖性根据这些邻居关系来定义。

通常，邻居由 空间依赖矩阵 $W$ 表示，如果 $\mathbf{s} j \in \mathscr{N}(\mathbf{s}_i)$，则元素 $w_{i, j}$ 非零，由于不能视自身为邻居，所以 $W$ 的对角线元素为零。$W$ 的非对角线元素可能与两者之间的距离 $\| \mathbf{s}_i − \mathbf{s}_j \|$ 成反比，或者采用一些其他基于空间接近度的渐进依赖性。如果存在邻域关系时为其赋值 $1$，否则赋值 $0$，则此时的 $W$ 被称为 邻接矩阵，而且当 $\mathbf{s}_j \in \mathscr{N}(\mathbf{s}_i)$ 时也有 $\mathbf{s}_i \in \mathscr{N}(s_j)$，则该矩阵是对称的。

考虑 $\mathbb{R}^2$ 中一个定义在 $D^L = \{(x, y) : x, y = 1, \ldots , 5\}$ 上的格元过程。格网区域内部的某个网格节点 $(x, y)$ ，其一阶邻居是相邻的四个节点，即 $\mathscr{N}(x, y) = \{(x − 1, y), (x, y − 1), (x + 1, y), (x, y + 1)\}$，如图：

其中网格节点 $\mathbf{s}_i$ 用 $×$ 表示，它的一阶邻居用 $\bullet$ 表示。位于格网边界上的节点，其邻居数将少于四个。

最常见的格元过程类型是 马尔可夫随机场 (MRF)，它需要在空间域 $\mathbb{R}^d$ 中满足马尔可夫条件概率假设。即如果对于所有 $\mathbf{s}_i \in D^L$，其条件概率满足下式的马尔可夫性，则格元过程 $\{Y(\mathbf{s}) : \mathbf{s} \in D^L\}$ 是一个马尔可夫随机场：

$$[Y(\mathbf{s}_i) \mid \mathbf{Y}(D^L \neg{\mathbf{s}_i})] = [Y(\mathbf{s}_i) \mid \mathbf{Y}(\mathscr{N}(\mathbf{s}_i))] \tag{19}$$

上式中 $\mathbf{Y}(A) \equiv \{Y (\mathbf{s}_j) : \mathbf{s}_j \in A\}$。

马尔可夫随机场是根据 式（19） 的条件概率定义的，这些条件概率表示只在相邻节点之间存在统计依赖性，这与地统计过程中的协方差函数（或变异函数）完全不同。具体来说，格元过程中的空间依赖性可以表示为：

$$[Y(\mathbf{s}_i) \mid \mathbf{Y}(D^L \neg{\mathbf{s}_i})] =\frac{\exp \{− f (Y (\mathbf{s}_i), \mathbf{Y}(\mathscr{N}(\mathbf{s}_i)))\}}{C} \tag{20}$$

其中 $C$ 是确保 式（20） 的右侧积分（或求和）为 $1$ 的归一化常数。式（20） 在统计力学中也被称为吉布斯随机场，因为在正则条件下，Hammersley-Clifford 定理将联合概率分布与吉布斯测度联系在了一起（Besag，1974）。函数 $f (Y (\mathbf{s}_i), \mathbf{Y}(\mathscr{N}(\mathbf{s}_i)))$ 被称为势能，它量化了邻居之间相互作用的强度。我们可以通过选择不同形式的势能函数，来定义多种马尔可夫随机场模型（Winkler，2003，第 3.2 节）。请注意，需要确保通过所有条件概率分布 $\{[Y (\mathbf{s}_i) \mid \mathbf{Y}(\mathscr{N}(\mathbf{s}_i))] : \mathbf{s}_i \in D^L \}$ 来定义模型整体的联合概率分布 $[\{Y (\mathbf{s}_i) : \mathbf{s}_i \in D^L\}]$（Kaiser 和 Cressie，2000）。

让我们重新审视之前在 $\mathbb{R}^2$ 中定义的规则网格，其一阶邻域结构见下图，请注意位于对角方向上的节点之间是条件独立的。因此，$D^L$ 可以分为两个子网格 $D^L_1$ 和 $D^L_2$，这样在给定 $D^L_1$ 中节点值的条件下，$D^L_2$ 中的节点值之间条件独立，反之亦然（Besag，1974 年；Winkler，2003 年，第 8.1 节）：

这形成了一个棋盘图案，其中给定由 $\circ$ 表示的 $D^L_2$ 节点处的值 $\{Y (\mathbf{u}) : \mathbf{u} \in D^L_2 \}$时，由 $\bullet$ 表示的 $D^L_1$ 节点处的 $\{Y(\mathbf{s}) : \mathbf{s} \in D^L_1 \}$ 之间是相互独立的，。

Besag (1974) 引入了 条件自回归 (CAR) 模型，它是一个根据条件均值和方差定义的高斯马尔可夫随机场。我们也建议读者参考 LeSage 和 Pace (2009) 中讨论的另外一种格元过程模型，被称为 同步自回归 (SAR) 模型，并将其与条件自回归模型进行比较。对于 $\mathbf{s}_i \in D^L$，$Y(\mathbf{s}_i)$，条件自回归模型被定义为由一阶矩和二阶矩定义的条件高斯分布：

$$\begin{align*} E(Y (\mathbf{s}_i) \mid \mathbf{Y}(\mathscr{N}(\mathbf{s}_i))) &= \sum_{\mathbf{s}_j \in \mathscr{N}(\mathbf{s}_i)} c_{i,j} Y (\mathbf{s}_) \tag{21} \\ \operatorname{var}(Y (\mathbf{s}_i) \mid \mathbf{Y}(\mathscr{N}(\mathbf{s}_i))) &= {\tau}^2_i \tag{22} \end{align*}$$

其中 $c_{i, j}$ 是空间自回归系数，其对角线元素 $c_{1,1} = \ldots = c_{n,n} = 0$， $\{ \tau^2_i \}$ 是位置 $\{ \mathbf{s}_i \}$ 处的尺度参数。在下面将提到的正则条件下，该定义使得联合概率分布最终是一个多元高斯分布，即：

$$\mathbf{Y} \sim \operatorname{Gau}(\mathbf{0}, (\mathbf{I-C})^{−1}\mathbf{M}) \tag{23}$$

其中 $\operatorname{Gau}( \boldsymbol{\mu }, \boldsymbol{\Sigma})$ 表示具有均值向量 $\boldsymbol{\mu}$ 和协方差矩阵 $\boldsymbol{\Sigma}$ 的高斯分布；矩阵 $\mathbf{M} \equiv \operatorname{diag}(\tau^2_1,\ldots,\tau^2_n)$ 是对角矩阵；前述所谓的正则条件是： 式（21） 中的系数 $\mathbf{C} \equiv \{c_{i, j}\}$ 必须使 $\mathbf{M}^{−1}(\mathbf{I-C})$ 是一个对称正定矩阵。对于一阶邻域结构（如上面 $\mathbb{R}^2$ 中的简单示例所示），精度矩阵是块对角矩阵，这使利用稀疏矩阵方法从高斯马尔可夫随机场中进行有效采样成为可能（Rue 和 Held，2005 年，第 2.4 节）。

格元过程对应的数据向量为 $\mathbf{Z} \equiv (Z(\mathbf{s}_1), \ldots, Z(\mathbf{s}_n))^\prime$，其中 $D^{L∗} \equiv \{\mathbf{s}_1,\ldots , \mathbf{s}_n \} \subset D^L$。根据 第 1 节 定义，数据模型为 $[\mathbf{Z} \mid \{Y(\mathbf{s}) : \mathbf{s} \in D^L \}]$。为了强调对参数 $\theta$ 的依赖性，我们重新启用其符号并将其写为 $[\mathbf{Z} \mid Y, \theta ]$。 现在，如果我们将格元过程模型记为 $[\{Y(\mathbf{s}) : \mathbf{s} \in D^L \} \mid \boldsymbol{\theta}] \equiv [Y \mid \boldsymbol{\theta} ]$，则通过最大化似然 $\mathscr{L} ( \boldsymbol{\theta} ) \equiv \int [\mathbf{Z} \mid Y, \theta ] [Y \mid \theta ] dY$，我们可以估计出 $\theta$。

对于空间预测任务，根据 $\mathbf{s}_0 \in D^L$ 和通过推断得出的 $\theta$，可以得到 $Y(\mathbf{s}_0)$ 的最优预测 $Y^∗(\mathbf{s}_0) \equiv E(Y (\mathbf{s}_0) \mid \mathbf{Z}, \theta)$（例如，Besag 等，1991 )。请注意，$\mathbf{s}_0$ 可能不属于 $D^{L*}$，因此即使在节点 $\mathbf{s}_0$ 处没有观测数据，$Y^*(\mathbf{s}_0)$ 也可以是 $Y(\mathbf{s}_0)$ 的最优预测。对空间随机过程 $Y$ 的未观测部分进行推断，对于格元过程和地统计过程同等重要。

2.3 空间点过程和随机集

空间点过程是随机位置 $D \equiv D^P \subset \mathscr{D}$ 的可数集合，与该随机点集密切相关的是被表示为 $\{ N(A) : A \subset \mathscr{D}\}$ 的计数过程，其中 $\mathscr{D}$ 索引所有可能感兴趣的位置，我们假设它是有界的。例如，如果 $A$ 是 $\mathscr{D}$ 的给定子集，并且 $A$ 中包含两个随机点 $\{\mathbf{s}_i\}$，则 $N(A) = 2$。由于 $D^P = \{\mathbf{s}_i\}$ 是随机的，而 $A$ 是固定的，$N(A)$ 是定义在非负整数上的随机变量。

显然，对于包含在 $\mathscr{D}$ 中任何子集（可能重叠） $\{A_j : j = 1, \ldots, m\}$ ，和任何 $m = 0, 1, 2,\ldots$，联合分布 $[N(A_1), \ldots, N(A_m)]$ 可以明确定义。空间依赖性可以通过 $\{A_j \}$ 之间的空间接近度表示。现在仅考虑两个固定子集 $A_1$ 和 $A_2$（即 $m = 2$），为了避免因为共享点引起的歧义，令 $A_1 \cap A_2$ 为空集，对于任何不相交的 $A_1$ 和 $A_2$，如果存在统计独立性，则两者之间不存在空间依赖性。即：

$$[N(A_1), N(A_2)] = [N(A_1)] [N(A_2)] \tag{24}$$

最基础的点过程被称为 泊松点过程 ，其具有 式（24） 的独立性，其相关的计数过程满足：

$$[N(A)] = \exp\{−\lambda (A)\} \frac{\lambda (A)^{N(A)}}{ N(A)!} ; A \subset \mathscr{D} \tag{25}$$

其中 $\lambda (A) \equiv \int_A \lambda(\mathbf{s})d \mathbf{s}$。在 式（25） 中，$\lambda (·)$ 是一个强度函数，被定义为：

$$\lambda (\mathbf{s}) \equiv \lim_{| \boldsymbol{\delta}_{\mathbf{s}} | \rightarrow 0} \frac{ E (Y (\boldsymbol{\delta}_{\mathbf{s}})) }{| \boldsymbol{\delta}_{\mathbf{s}} | } \tag{26}$$

其中 $\boldsymbol{\delta}_\mathbf{s}$ 是以 $\mathbf{s} \in \mathscr{D}$ 为中心的小集合，$| \boldsymbol{\delta}_{\mathbf{s}} |$ 为其体积。

在 式（25） 中，对于所有 $\mathbf{s} \in \mathscr{D}$，如果存在 $\lambda(\mathbf{s}) \equiv \lambda$ 的特殊情况，则导致齐次泊松点过程，其模拟如 图 3 所示。该模拟是使用与 式（25） 等效的概率表示获得的计数随机变量 $N(\mathscr{D})$（ $\mathscr{D} = [0, 1] × [0, 1]$）。然后，基于 $N(\mathscr{D})$ ，根据如下均匀分布独立、可识别地模拟了 $\{\mathbf{s}_1, \ldots , \mathbf{s}_N(\mathscr{D})\}$。

$$[\mathbf{u}]=\left\{\begin{array}{l} \frac{1}{\lambda(\mathscr{D})} ; \quad \mathbf{u} \in \mathscr{D} \\ 0 ; \text { elsewhere } \end{array}\right . \tag{27}$$

上述表示也解释了为何齐次泊松点过程通常被称为 完全空间随机 (CSR) 过程，以及为什么它被用做点过程中不存在空间依赖性的测试基线。也就是说， 在将空间模型拟合到点模式之前，通常会对点模式源自完全空间随机过程的原假设进行检验。拒绝完全空间随机过程然后证明空间相关点过程与数据的拟合才是合理的（例如，Ripley，1981；Diggle，2013）。

图 3: 在单位正方形 $\mathscr{D}=[0,1]\times[0,1]$ 上的一个同质点过程（式 25），其中参数 $\lambda= 50$，对于此实现 $N(\mathscr{D})=46$。

很多关于点过程的早期研究致力于建立对各种类型 CSR 过程偏离敏感的测试统计数据（例如，Cressie，1993，第 8.2 节）。随后是研究人员定义并估计空间相关性度量，例如二阶强度函数和 $K$ -函数（例如，Ripley，1981 年，第 8 章），其中推断通常是根据这些函数的矩估计方法。后来出现了更有效的基于似然的推断；Baddeley 等 (2015) 对这些点过程方法进行了全面回顾。从建模角度来看，很多研究特别关注对数高斯 Cox 点过程；在该模型中， 式（25） 中的 $\lambda (·)$ 是随机的，因此 $\{\log(\lambda (\mathbf{s})) : \mathbf{s} \in \mathscr{D}\}$ 是一个高斯过程（例如，Møller 和 Waagepetersen，2003）。该模型自然导致 $\lambda (·)$ 及其参数的分层贝叶斯推断（例如，Gelfand 和 Schliep，2018）。

如果属性过程 $\{Y (\mathbf{s}_i) : \mathbf{s}_i \in D^P\}$ 包含在空间点过程 $D^P$ 中，则可以获得 标记点过程（例如，Cressie，1993，第 8.7 节）。例如，对天然森林的研究，其中位置 $\{\mathbf{s}_i\}$ 和树木的大小 $\{Y (\mathbf{s}_i)\}$ 被概率建模在一起，产生了标记点过程，其中 “标记” 过程是一个关于树大小的空间过程 $\{Y (\mathbf{s}_i) : \mathbf{s}_i \in D^P\}$ 。现在，由 式（10） 的贝叶斯规则（其中 $Y$ 和 $D (= D^P)$ 都是随机的），应该被用于对 $Y$ 和 $D^P$ 作出推断（通过预测分布 $[Y, D^P \mid \mathbf{Z}]$）。这里，$\mathbf{Z}$ 由树木数量、树木位置、数目大小等测量值组成，如 式（5） 中所示。通过边缘化，我们可以得到空间点过程 $D^P$ 的预测分布 $[D^P \mid \mathbf{Z}]$。

空间点过程是 随机集（Random Set） 的特例。随机集是欧几里德空间中的随机量，由 Matheron (1975) 严格定义。一些地质过程更自然地建模为 以集合做为值 的现象（例如，矿化相），但是随机集过程的推断落后于空间点过程的推断。很难根据 集合值 数据来定义似然，这阻碍了统计上有效的推断；尽管如此，基本的矩估计方法通常还是可用的。允许从 集合值 数据中进行推断的最著名随机集是 布尔模型（Boolean Model）（例如，Cressie 和 Wikle，2011，第 4.4 节）。

2.4 多元空间过程

前面的小节介绍了单个空间统计过程，但是随着模型成为复杂世界的更真实的表示，需要表达多个过程之间的相互作用。这通过对矢量值 “地统计过程” $\{ \mathbf{Y}(\mathbf{s}) : \mathbf{s} \in D^G \}$ 和矢量值 “格元过程” $\{ \mathbf{Y}(\mathbf{s}_i) : \mathbf{s}_i \in D^L \}$ 建模最直接地看出，其中 $k$ -维向量 $\mathbf{Y}(\mathbf{s}) \equiv (Y_1(\mathbf{s}),\ldots ,Y_k(\mathbf{s}))^\prime$ 表示通用位置 $\mathbf{s} \in D$ 处的多个过程。可以表示第 2.3 节中讨论的向量值空间点过程作为一组 $k$ 个点过程，$\{ { \mathbf{s}_{1,i} },\ldots, { \mathbf{s}_{k,i} } \}$，这些在 Baddeley 等 （2015 年，第 14 章）中有介绍。如果我们采用分层统计建模方法，则可以构建多元空间过程，其组成单变量过程可以来自前三个小节中介绍的三种空间过程中的任何一种。这是因为，在层次结构的更深层次上，核心多元地统计过程可以控制任何类型过程的空间依赖性，这使得混合多元空间统计过程成为可能。

在下文中，我们简要描述了构建多元地统计过程的两种方法，一种基于联合方法，另一种基于条件方法。我们考虑 $k = 2$ 的情况，即双变量空间过程 $\{(Y_1(\mathbf{s}),Y_2(\mathbf{s}))^\prime : \mathbf{s} \in D^G\}$，以供说明。

（1）联合方法

对于 $\mathbf{s} \in D^G$，联合方法涉及从 $\boldsymbol{\mu} (\mathbf{s}) \equiv ( \mu_1(\mathbf{s}), \mu_2(\mathbf{s}))^\prime \equiv (E(Y_1(\mathbf{s})), E(Y_2(\mathbf{s})))^\prime$ 直接构建有效的空间统计模型，并且

$$\operatorname{cov}(Y_l(\mathbf{s}),Y_m(\mathbf{u})) \equiv C_{lm}(\mathbf{s,u}); l, m = 1, 2 \tag{28}$$

对于 $\mathbf{s,u} \in D^G$。双变量过程均值 $\boldsymbol{\mu} (·)$ 通常建模为向量线性回归；因此，一旦选择了适当的回归变量，就可以直接对双变量均值进行建模

类似于单变量情况，协方差和互协方差函数集 $\{C_{11}(·,·), C_{22}(·,·), C_{12}(·,·), C_{21}(·,·)\}$ 必须满足正-双变量地统计模型有效的确定性条件，需要注意的是，一般来说，$C_{12}(\mathbf{s,u}) \neq C_{21}(\mathbf{s,u})$。有几类有效模型表现出对称的交叉依赖性，即 $C_{12}(\mathbf{s,u}) = C_{21}(\mathbf{s,u})$，例如共同区域化的线性模型 (Gelfand et al., 2004)。当矿体中存在优先成矿时，这些模型不是合理的矿石储量估算模型。

（2）条件方法

条件方法 (Cressie and ZammitMangion, 2016) 中 $k$ 个过程中的每一个都是有向无环图的一个节点，该图在给定剩余过程的情况下指导任何过程的条件依赖性。再次考虑双变量情况（即 $k = 2$），其中只有两个节点，使得 $Y_1(·)$ 在节点 1，$Y_2(·)$ 在节点 2，并且从节点 1 到节点 2 声明一条有向边。那么对联合分布建模的合适方法是通过

$$[Y_1(·),Y_2(·)] = [Y_2(·) \mid Y_1(·)][Y_1(·)] \tag{29}$$

其中 $[Y_2(·) \mid Y_1(·)]$ 是 $[Y_2(·) \mid \{Y_1(s) : \mathbf{s} \in D^G \}]$。

$[Y_1(·)]$ 的地统计模型只是一个基于均值函数 $\mu_1(·)$ 和有效协方差函数 $C_{11}(·,·)$ 的单变量模型，这在 “地统计过程” 中进行了讨论。现在假设 $Y_2(·)$ 取决于 $Y_1(·)$ 如下： 对于 $\mathbf{s,u} \in D^G$，

$$\begin{align*} E[Y_2(s) \mid Y_1(·)] &\equiv \mu_2(s) + ∫ D^G b(s, v) (Y_1(v) − \mu_1(v)) dv \tag{30}\\ \operatorname{cov}(Y_2(s),Y_2(u) \mid Y_1(· )) &\equiv C_{2 \mid 1}(s, u) \tag{31} \end{align*}$$

其中 $C_{2 \mid 1}(·,·)$ 是有效的单变量协方差函数，$b(·,·)$ 是可积交互函数。如果假设 $(Y_1(·),Y_2(·))^\prime$ 是双变量高斯过程，则由 式（30） 和 式（31） 给出的条件矩假设会随之而来。

Cressie 和 Zammit-Mangion (2016) 表明，根据 式（30） 和 式（31），

$$\begin{align*} C_{12}(\mathbf{s,u}) &= ∫ D^G C_{11}(\mathbf{s,v})b(\mathbf{u,v}) d \mathbf{v} \tag{32}\\ C_{21}(\mathbf{s,u}) &= \int D^G C_{11}(\mathbf{v,u})b(\mathbf{v,s}) d \mathbf{v} \tag{33} \\ C_{22}(\mathbf{s,u}) &= C_{2 \mid 1}(\mathbf{s,u}) + \int_{D^G} \int_{D^G} b(\mathbf{s,v}, \mathbf{v})C_{11}(\mathbf{v,w})b(\mathbf{u,w}) d \mathbf{v} d \mathbf{w} \tag{34} \end{align*}$$

对于 $\mathbf{s,u} \in D^G$。连同 $\mu_1(·)$、 $\mu_2(·)$ 和 $C_{11}(·,·)$，这些函数 式（32）– 式（34） 定义了一个有效的双变量地统计过程 $[Y_1(·),Y_2(·)]$。

条件方法的一个显著特性是，如果 $b(\mathbf{s,u}) \neq b(\mathbf{u,s})$，则会发生非对称交叉依赖（即 $C_{12}(\mathbf{s,u}) \neq C_{21}(\mathbf{s,u})$ ）。

总之，条件方法允许通过简单地指定 $\mu (·) = ( \mu_1(·), \mu_2(·))^\prime$ 和两个有效的单变量协方差函数 $C_1(·,·)$ 和 $C_{2 \mid 1}(·,·)$ 来有效地执行多变量建模。条件方法的优点是只需要指定单变量协方差函数（对此有大量研究；例如，Cressie 和 Wikle，2011 年，第 4 章），并且只有交互函数 $b(·,·)$ 需要假设可积性（Cressie 和 Zammit-Mangion，2016）。

3 空间离散化

尽管地统计过程是在连续空间域 $D^G$ 上定义的，但由于计算和数学方面的考虑，这可能会限制统计推断的实际可行性。例如，基于 $n$ 维数据向量的克里金法涉及 $n × n$ 协方差矩阵的求逆，这需要 $n^3$ 次浮点运算和可用内存中的 $n^2$ 次存储。对于大型空间数据集来说，这种计算成本可能让人望而却步；因此，利用空间离散化以实现空间模型的可扩展计算已经成为一个活跃的研究领域。

在实际应用中，空间统计推断需要只需要达到有限空间分辨率即可，因此，许多方法希望通过将空间域 $\mathscr{D}$ 划分离散网格的来利用这一特点，如 图 4 所示。作为这种离散化的结果，地统计过程可以转而通过格元过程来近似，例如高斯马尔可夫随机场（例如，Rue 和 Held，2005 年，第 5.1 节），但有时这会导致不希望的离散化误差和伪影。目前，已经开发出了一些更复杂的方法，来获得在不规则网格上评估地统计（即连续索引）空间过程的高精度近似。

令原始域 $\mathscr{D}$ 有界，并假设它可以被细分为 $\{A_j \subset \mathscr{D} : j = 1, . . . , m\}$ 等小的、无缝无叠的基本区域单元 (Nguyen et al., 2012)，即有 $\mathscr{D} = \cup^m_{ j=1} A_j$ ，且对于任何 $j \neq k \in \{1, \ldots , m\}$, $A_j \hat A_k$ 是空集。图 4 给出了三角形基本区域单元的示例。可以在基本区域单元上定义空间基函数 $\{\phi(·) : l = 1, \ldots , r\}$。例如，Lindgren 等(2011) 使用了三角基函数，其中 $r > m$；而固定秩克里金法（FRK；Cressie 和 Johannesson，2008 年；Zammit-Mangion 和 Cressie，2021）可以为 $r < m$ 使用各种不同的基函数，包括多分辨率小波和双平方。

Vecchia 近似（例如，Datta 等人，2016 年；Katzfuss 等人，2020 年）也是使用离散点网格 $D^L \subset D^G \subset \mathscr{D}$ 定义的，其中包括观测数据的位置 $D^{G∗} = \{\mathbf{s}_1, \ldots , \mathbf{s}_n\}$ 和预测位置 $\{\mathbf{s}_{n+1}, \ldots , \mathbf{s}_{n+p}\}$。令 $[\mathbf{X}] \equiv [\mathbf{Z},Y ]$，其中数据 $\mathbf{Z}$ 和空间过程 $Y$ 与网格 $D^L \equiv \{ \mathbf{s}_1, \ldots , \mathbf{s}_n, \mathbf{s}_{n+1}, \ldots , \mathbf{s}_{n+p} \}$ 相关联。其联合分布和条件分布可以被分解为乘积形式：

$$[\mathbf{X}] =\prod^{n+p}_{i= 1} [X(\mathbf{s}_i) \mid X(\mathbf{s}_1),\ldots , X(\mathbf{s}_{i−1})] \tag{35}$$

在上一节中，空间坐标集 $D^L$ 没有固定顺序。但 Vecchia 近似要求对 $\{\mathbf{s}_1,\ldots, \mathbf{s}_{n+p}\}$ 进行人工排序。如果令排序结果为 $\{\mathbf{s}_{(1)}, \dots, \mathbf{s}_{(n+p)}\}$, 其中元素的邻居集合为 $N (\mathbf{s}_{(i)}) \subset \{\mathbf{s}_{(1)}, \dots, \mathbf{s}_{(i−1)}\}$，这类似于 “格元过程”，只是这些邻域关系不是互惠的：如果对于 $j < i$，$\mathbf{s}_{(j)}$ 属于 $N(\mathbf{s}_{(i)})$，则 $\mathbf{s}_{(i)}$ 不能属于 $N(\mathbf{s}_{(j)})$ 。作为 Vecchia 近似的一部分，选择了邻居数量的固定上限 $q \ll n$ 。即 $| N(\mathbf{s}_{(i)}) | \leq q$，因此由 $\{N (\mathbf{s}_i) : i = 1, \ldots, n + p \}$ 形成的格元是有向无环图，这导致 $\mathscr{D}$ 中的偏序（Cressie 和 Davidson，1998）。

然后由 式(35) 给出的联合分布 $[\mathbf{X}]$ 近似为：

$$\prod^{n+ p}_{i= 1} [X (\mathbf{s}_{(i)}) \mid \mathbf{X}(N (\mathbf{s}_{(i)}))] \equiv [\tilde{\mathbf{X}}] \tag{36}$$

这是一个偏序的马尔可夫模型 (POMM; Cressie and Davidson, 1998)。Vecchia 近似 $[\tilde{\mathbf{X}}]$ 是来自原始、不可数、无序索引集 $D^G$ 上的有效空间过程的分布（Datta 等人，2016），这意味着它可以作为地统计过程模型，并且具有相当大的计算优势。例如，它可以用于分层点过程模型中的随机对数强度函数 $\log(\lambda (\mathbf{s}))$，或者，它也可以与其他过程组合以定义 “多元空间过程” 中描述的模型。但在所有这些情况下，由此产生的预测过程 $[\tilde{\mathbf{X}} \mid \mathbf{Z}]$，都是真实预测过程 $[ \mathbf{x} \mid \mathbf{Z}]$ 的近似值.

4 时空过程

标题为“多变量空间过程”的部分介绍了以矢量形式编写的过程，

$$\mathbf{Y}(\mathbf{s}) \equiv (Y_1(\mathbf{s}), \ldots,Y_k(\mathbf{s}))^\prime; \mathbf{s} \in D^G \tag{37}$$

在那一节中，我们区分了多元空间统计建模的联合方法和条件方法，并且在条件方法下，我们使用有向无环图来给出多元空间相关性的蓝图。

现在，考虑一个时空过程，

$$\{ Y (\mathbf{s};t) : \mathbf{s} \in D^G; t \in \mathcal{T} \}$$

其中 $\mathcal{T}$ 是时间索引集。显然，如果 $\mathcal{T} = \{1, 2, \ldots\}$，那么 式（38） 就变成了时间序列的空间过程，$\{Y (\mathbf{s}; 1),Y (\mathbf{s}; 2),\ldots : \mathbf{s} \in D^G \}$。如果 $\mathcal{T} = \{1, 2, \ldots, k \}$，并且我们定义 $Y_j(\mathbf{s}) \equiv Y(\mathbf{s}; j)$，对于 $j = 1,\ldots, k$，则由此产生的时空过程可以表示为 式（37） 给出的多元空间过程.毫不奇怪，与多变量空间过程一样，时空过程也会出现相同的统计依赖性建模方法二分法（即联合与条件）。

描述 $Y$ 在任何时空“位置” $(\mathbf{s};t)$ 和任何其他位置 $(\mathbf{u};v)$ 之间的所有可能协方差，相当于将 “时间” 视为要添加到 $d$ 维欧几里得空间 $\mathbb{R}^d$ 的另一个维度。采用这种方法，时空统计相关性可以通过协方差函数在 $(d + 1)$ 维空间中表示

$$C(\mathbf{s};t, \mathbf{u}; v) \equiv \operatorname{cov}(Y (\mathbf{s};t),Y (\mathbf{u}; v)); \mathbf{s,u} \in D^G, t, v \in \mathcal{T} \tag{39}$$

当然，时间维度与空间维度的单位不同，未来是不可观察的，其解释也不同。因此，基于 式（39） 的空间和时间的联合建模必须小心进行，以在这种时空建模的描述性方法中考虑时间维度的特殊性质

从当前和过去的时空数据 $\mathbf{Z}$，预测 $Y$ 的过去值称为平滑，预测当前 $Y$ 未观察到的值称为过滤，预测未来 $Y$ 的值称为预测。卡尔曼滤波器 (Kalman, 1960) 的开发是为了使用一种识别时间维度顺序的方法提供对当前状态的快速预测。当收到一组新的当前数据时，今天的过滤值在第二天就变“旧”了。使用我们现在将要描述的动态方法，卡尔曼滤波器非常快速地用今天的数据更新昨天的最优过滤值，以获得当前的最优过滤值。

描述动态方法的最佳方式是将空间域离散化。上一节“空间离散化”描述了实现此目的的多种方法；在这里，我们将考虑在计算机内存中存储属性和位置信息最自然的离散化，即像素或体素（“体积元素”的缩写）的精细分辨率格元 $D^L$。替换 $\{Y (\mathbf{s};t) : \mathbf{s} \in D^G, t = 1, 2, \ldots \}$ 与 $\{Y (\mathbf{s};t) : \mathbf{s} \in D^L, t = 1, 2, \}$，其中 D^L \equiv {\mathbf{s}_1, \ldots,\mathbf{s}_m\} 是构成 $D^G$ 的小面积（或小体积）元素的质心。通常这些元素的面积被指定为相等的，由规则的网格定义。正如我们在下面解释的那样，这允许一种动态方法来为离散时空立方体 $\{\mathbf{s}_1, \ldots,\mathbf{s}_m \} × \{1, 2, \ldots\}$ 上的时空过程 $Y$ 构建统计模型。

定义 $Y_t \equiv (Y (\mathbf{s};t) : \mathbf{s} \in D^L)^\prime$，它是一个 $m$ 维向量。由于时间顺序，我们可以将 $\{Y (\mathbf{s};t) : \mathbf{s} \in D^L, t = 1, ..., k\}$ 从 $t = 1$ 到当前时间 $t = k$ 的联合分布写为

$$[\mathbf{Y}_1, \mathbf{Y}_2, \ldots, \mathbf{Y}_k] = [\mathbf{Y}_1][\mathbf{Y}_2 \mid \mathbf{Y}_1] \ldots [\mathbf{Y}_k \mid \mathbf{Y}_{k-1}, \dots, \mathbf{Y}_2, \mathbf{Y}_1]$$

与 式（35） 具有相同的形式。请注意，这种空间和时间的条件建模是一种自然的方法，因为时间是完全有序的。下一步是做出马尔可夫假设，因此 式（40） 可以写成

$$[\mathbf{Y}_1, \mathbf{Y}_2, \ldots, \mathbf{Y}_k] = [\mathbf{Y}_1] \prod^{k}_{j= 2} [\mathbf{Y}_j \mid \mathbf{Y}_{j-1}]$$

这与我们之前在“格元过程”中讨论的马尔可夫属性相同，只是它现在应用于完全有序的一维域，$\\mathcal{T} = \{1, 2, \ldots \}$, $N(j) = j −1$。马尔可夫假设使我们的方法动态化：它说现在以过去为条件，实际上只取决于“最近的过去”。即，由于 $N(j) = j − 1$，因数 $[\mathbf{Y}_j \mid \mathbf{Y}_{j−1}, \ldots, \mathbf{Y}_2, \mathbf{Y}_1] = [\mathbf{Y}_j \mid \mathbf{Y}_{ j−1}]$，这导致模型 式（41）。

有关 式（39） 给出的描述性方法中使用的模型类型和 式（41） 给出的动态方法中使用的模型类型的更多信息，请参见 Cressie 和 Wikle（2011 年，第 6-8 章）和 Wikle 等人. （2019 年，第 4 章和第 5 章）。对来自这些过程的观察结果的统计分析称为时空统计。使用 R 软件从时空数据推断（估计和预测）可以在 Wikle 等人中找到。 （2019）。

5 结论

空间统计方法对估计或预测中的不确定性提供了经过良好校准的量化，从而使其与地理和环境科学中的其他空间分析方法有了明显区别。

（1）总结

在本文中，空间科学现象中的不确定性由空间过程模型 $\{Y(\mathbf{s}) : \mathbf{s} \in D\}$ 表示，该空间过程定义在 $\mathbb{R}^d$ 中可能随机的空间域 $D$ 上；而观测 $\mathbf{Z}$ 中存在的测量不确定性在数据模型中表示。

在 第 1 节 中，我们了解了如何使用贝叶斯规则来组合上述两个模型来计算统计推断所需的总体不确定性。

除了一些例外情况（例如，Cressie 和 Kornak，2003 年），空间统计模型很少考虑 $\mathscr{D}$ 中的位置测量误差情况。这里我们关注位置误差的空间统计模型：将观察到的位置写为 $D^∗ \equiv \{\mathbf{u}_i : i = 1, \ldots, n\}$；在这种情况下，数据模型的一部分是 $[D^* \mid D]$，过程模型的一部分是 $[D]$。最后，数据由位置和属性组成，并且是 $\mathbf{Z} \equiv \{(\mathbf{u}_i, Z(\mathbf{u}_i)) : i = 1, \ldots , n\}$，空间过程模型为 $[Y, D]$，数据模型为 $[\mathbf{Z} \mid Y，D]$。然后使用 式（10） 给出的贝叶斯规则从预测分布 $[Y, D \mid \mathbf{Z}]$。

空间过程模型主要分为三种类型：

• 地统计过程：过程 $Y$ 中存在不确定性，且在 $D = D^G$ 中被连续索引；
• 格元过程：过程 $Y$ 中同样存在不确定性，但 $Y$ 在可数空间域 $D = D^L$ 上被索引；
• 点过程：在空间域 $D = D^P$ 中存在不确定性。

多个空间过程可以相互影响，形成多元空间过程。重要的是，过程可以随时间和空间变化，形成时空过程。

随着空间数据集的规模急剧增加，越来越多的注意力集中在空间统计模型的可扩展计算上。特别受关注的是使用 “空间离散化” 来近似连续空间域 $D^G$ 的方法。

（2）最新进展

我们认为值得一提的是空间统计方面的其他最新进展，这里进行简短的介绍：

• 关于非平稳性：物理障碍有时会中断空间上非常接近的位置之间的统计关联。已经开发了障碍模型（Bakka 等人，2019 年）来解释空间相关函数中的这些类型的不连续性。对空间过程模型中的非平稳性、各向异性和异方差性进行建模的其他方法是一个活跃的研究领域。

• 关于先验的选择：通常很难为平稳空间过程的参数选择合适的先验分布，例如其相关长度尺度。惩罚复杂性先验（Simpson 等人，2017 年）是一种通过支持参数值来鼓励简约的方法，这些参数值会产生与数据一致的最简单模型。点过程或非高斯格元模型的似然函数在分析和计算上都是难以处理的。已经开发了替代模型、仿真器和准似然来近似这些棘手的似然（Moores 等人，2020 年）。

• 关于多变量建模：Copula （Krupskii 和 Genton，2019 年）是对多变量数据中的空间依赖性进行建模的替代方法，尤其是当数据是非高斯数据时。可能出现非高斯性的一个领域是对极端事件之间的空间关联进行建模，例如温度或降水（Tawn 等，2018 年；Bacro 等，2020 年）。

• 关于小样本建模：来自有限数量观测的可数空间随机变量，可以通过离散化过程来精细化。在现代计算环境中，这是进行空间统计推断（包括克里金法）的关键。