固定秩克里金法

【摘 要】 非常大的空间数据集的空间统计具有挑战性。数据集的大小 $n$ 会导致计算最优空间预测变量（例如克里金法）出现问题，因为其计算成本为 $n$ 的三次方。此外，大型数据集通常是在大型空间域上定义，因此感兴趣的空间过程通常在该域上表现出非平稳行为。 通过使用一组固定数量的基函数，可以定义一个灵活的非平稳协方差函数族，这产生了我们称为 “固定秩克里金法” 的空间预测方法。具体来说，固定秩克里金法就是此类非平稳协方差函数支撑下的克里金法。当 $n$ 非常大时，它依赖于计算简化，以获得隐空间过程的空间最佳线性无偏预测器及其均方预测误差。基于最小化加权 Frobenius 范数的方法产生协方差函数参数的最佳估计量，然后将其代入固定秩克里金方程。新方法适用于在整个地球上观测到的非常大的臭氧数据集，其中 $n$ 约为数十万。

【原 文】 Cressie, N. and Johannesson, G. (2008) ‘Fixed rank kriging for very large spatial data sets: Fixed Rank Kriging’, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(1), pp. 209–226. Available at: https://doi.org/10.1111/j.1467-9868.2007.00633.x.

1 引言

克里金法或空间最佳线性无偏预测 (BLUP) 在地球和环境科学中非常流行，有时被称为最佳插值。 Matheron (1962)^[26] 创造了术语 “克里金法” 以纪念南非采矿工程师 D. G. Krige (Cressie, 1990) ^[6]。通过协方差函数（或变异函数）对空间可变性进行内部量化，克里金法可以根据不完整和含噪声的空间数据生成最佳预测图和相关预测标准差（例如 Cressie (1993) ^[7]，第 3 章）。有时空间数据的获取成本很高（例如，为估算石油储量而钻井），在这种情况下，样本量 $n$ 通常较小，克里金法可以直接执行。但随着卫星遥感平台的普及，数据库范式从小到大，通常是每天千兆字节的量级。求解克里金方程直接涉及 $n \times n$ 协方差矩阵 $\boldsymbol{\Sigma}$ 的求逆，其中 $n$ 个数据可能需要 $\mathcal{O}(n^3)$ 计算才能获得 $\boldsymbol{\Sigma}^{-1}$。在这些情况下，无法直接对海量数据集进行克里金法处理。

本文的目标是开发将克里金法的计算成本降低到 $\mathcal{O}(n)$ 的方法。

（1）移动窗口法和稀疏化法方法

即使是几千数量级的空间数据集也会导致计算速度减慢，因此：

• 移动窗口方法：数据子集方法由 Haas (1995) ^[13] 的移动窗口方法形式化，不过在窗口内拟合的局部协方差函数，在较大空间滞后时会产生不相容的协方差。
• 协方差矩阵稀疏化：当协方差函数的变程有限时，协方差矩阵 $\boldsymbol{\Sigma}$ 通常是稀疏的，因此 $\boldsymbol{\Sigma}^{-1}$ 可以通过使用稀疏矩阵技术获得。
• 精度矩阵稀疏化：Rue 和 Tjelmeland (2002) ^[33] 直接将空间过程近似为包裹在环面上的高斯马尔可夫随机场，并将精度矩阵 $\boldsymbol{\Sigma}^{-1}$ 近似为稀疏矩阵。

（2）局部克里金方法

当数据集很大（数万量级）到非常大（数十万量级）时，直接克里金法可能会失效，通常会使用局部克里金邻域（例如 Cressie (1993) ^[7]，第 131–134页)。已有成果包括 根据正交基给出等效表示然后采取截断基的形式、 协方差锥化、 近似迭代方法（如共轭梯度）、 使用归纳变量实现稀疏近似 或 用更小的空间内填充位置集替换数据位置。

• 最近研究的一个途径是近似克里金方程（Nychka 等，1996 年 ^[30]，2002 年^[31]；Nychka，2000 年 ^[29]；Billings 等，2002 年^[2][3]；Furrer 等，2006 年^[12]；Quinonero-Candela 和 Rasmussen，2005 年 ^[32]）。 Kammann 和 Wand (2003)^[23] 在拟合一类他们称为地理加性模型的空间模型时采用了最后一个想法。

（3）多分辨率克里金法

另一种方法是选择协方差函数的类别，即使空间数据集很大，也可以对其进行精确的克里金法（例如 Huang 等 (2002) ^[18]、Johannesson 和 Cressie (2004a) ^[19] 以及 Johannesson 等 (2007) ^[21]）)。在这些论文中，构建了一个多分辨率空间（和时空）过程，以便克里金法可以迭代和极快地计算，计算复杂度与数据大小成线性关系。在空间案例中，Johannesson 和 Cressie (2004a) ^[19] 实现了 $10^8$ 数量级的加速比直接求解克里金方程。他们可以在大约 $3$ 分钟内为 $n=160 000$ 计算出整个地球的最佳空间预测变量及其相关的均方预测误差。拥有允许精确计算的空间模型的一个优点是，无需担心近似克里金预测变量有多接近近似均方预测误差为相应的理论值。此类精确方法有两个重要问题需要确定：一是用于克里金法的空间协方差函数有多灵活？二是如何拟合该协方差函数？

（4）本文途径

上面提到的多分辨率模型隐含的空间协方差是 “非平稳的” 和 “块状的” 的。而在本文中，我们将使用一种不同的方法来实现最佳空间预测的量级加速，该方法既可以使用非常灵活的协方差函数，也可以选择是否平滑，这取决于通过空间数据所展示的空间依赖类型（与 Tzeng 等（2005）^[38] 的方法相反）。我们将证明存在一类非常丰富的空间协方差函数，从中可以精确地执行大型空间数据集的克里金法，其计算成本仅为 $\mathcal{O}(n)$。

我们将考虑一类 $n \times n$ 的协方差矩阵 $\boldsymbol{\Sigma}$，通过我们的方法使得 $\boldsymbol{\Sigma}^{-1}$ 可以通过对 $r \times r$ 矩阵求逆得到。由于其中 $r$ 是固定的，所以我们将该方法称为固定秩方法。例如：在第 4 节中给出的臭氧 (TCO) 数据应用中，$n$ 为 $173 405$，而 $r$ 被选择为 $396$。在 第 2.3 节 中，我们给出了该方法的计算复杂度为 $\mathcal{O}(nr^2)$，它仅随样本量线性增加。

此外， 假设数据集是覆盖全球的卫星遥感数据，那么数据中的任何空间依赖性几乎肯定会在全球范围内的异质性。

本文方法的创新之处在于：同时解决了 “数据集大小” 和 “空间异质性” 这两个问题。研究成果是我们称之为固定秩克里金法 (FRK) 的空间 BLUP 过程，它依赖于对 $r$ 固定且独立于 $n$ 的 $r \times r$ 矩阵求逆。

为了完整起见，我们提到了另一种基于平滑样条的空间预测方法。与克里金法相比，平滑样条不依赖于空间随机过程（其协方差函数必须建模、拟合并用于计算最佳预测变量），但有结点（knots）和平滑（smoothing）等参数需要确定，而且空间数据集大小也会导致计算困难。Hastie (1996) ^[14] 和 Johannesson 和 Cressie (2004b) ^[20] 为海量数据集开发了低阶样条平滑器。

（5）本文贡献

（1） 在协方差结构建模方面：要进行 FRK，必须指定（非平稳）协方差函数的形式；我们提出了一类足够灵活、可以对多个尺度的空间变化进行建模的协方差函数，它能够产生 $n \times n$ 协方差矩阵 $\boldsymbol{\Sigma}$，且其对应的精度矩阵可直接计算。

（2） 在空间预测方面：最小化均方预测误差的空间 BLUP 在各种矩阵计算中涉及 $\boldsymbol{\Sigma}^{-1}$ 的计算问题；我们表明，FRK 的计算成本与 $n$ 呈线性关系，远远优于 传统的 $\mathcal{O}(n^3)$。空间协方差函数通过最小化加权 Frobenius 范数来拟合经验协方差。

本文结构如下：

• 第 2 节 介绍了克里金法并给出了定义 FRK 的方程。
• 在第 3 节中，研究了 FRK 中使用的一类非平稳协方差函数，包括如何找到最适合数据的协方差函数。
• 第 4 节将该方法应用于 TCO 数据，其中 $n = 173 405$；通过直接求逆数据的 $n \times n$ 理论协方差矩阵来进行克里金法是不可能的。
• 第 5 节包含讨论和结论，后面是技术附录 A。

2 Kriging：最优线性空间预测

在本节中，我们将介绍克里金法的表示法，并将其等同于空间设置中的 BLUP。当空间数据集很大时，通常不可能精确计算克里金法。在本节的后半部分，我们展示了选择特定类别的非平稳空间协方差如何允许快速计算克里金预测变量（即空间 BLUP）和克里金标准差（即均方根预测误差）。

2.1 克里金方程

设 $Y(\mathbf{s}) : \mathbf{s} \in \mathcal{D} \subset \mathbb{R}^d\}$ 为实值空间过程。我们希望在包含测量误差的数据基础上，对真实的 $Y$ 过程进行推断；考虑如下实际观测的过程 $Z(\cdot)$：

$$Z(\mathbf{s}) \equiv Y(\mathbf{s}) + \varepsilon(\mathbf{s}), \mathbf{s} \in \mathcal{D} \tag{2.1}$$

其中 $\{\varepsilon(\mathbf{s}) : \mathbf{s} \in \mathcal{D} \}$ 是零均值的空间白噪声过程，过程方差为 $\operatorname{var}{\varepsilon(\mathbf{s})} = \sigma^2 v(\mathbf{s}) \in (0, \infty)$，其中 $\mathbf{s} \in \mathcal{D}$，$\sigma^2 > 0$ 且 $v(\cdot)$ 已知。事实上，过程 $Z(\cdot)$ 仅在有限数量的空间位置 $\{\mathbf{s}_1,\ldots ,\mathbf{s}_n\}$ 已知；将可用数据的向量定义为

$$\mathbf{Z} \equiv (Z(\mathbf{s}_1),\dots,Z(\mathbf{s}_n))^{\prime} \tag{2.2}$$

假设隐过程 $Y(\cdot)$ 具有线性均值结构，

$$Y(\mathbf{s}) = \mathbf{t}(\mathbf{s})^{\prime} \boldsymbol{\alpha} + ν(\mathbf{s}), \qquad \mathbf{s} \in \mathcal{D} \tag{2.3}$$

其中 $\mathbf{t}(\cdot) \equiv (\mathbf{t}_1(\cdot),\ldots ,\mathbf{t}_p(\cdot))^{\prime}$ 表示已知协变量的向量；$\boldsymbol{\alpha} \equiv (\alpha_1,\ldots ,\alpha_p)^{\prime}$ 表示未知的系数 ，$ν(\cdot)$ 为零均值的过程，其方差为 $0 < \operatorname{var}\{ν(\mathbf{s})\} < \infty$，所有的 $\mathbf{s} \in \mathcal{D}$，并且具有不一定平稳的空间协方差函数：

$$\operatorname{cov}\{ν(\mathbf{u}), ν(\mathbf{v})\} \equiv C(\mathbf{u,v}),\qquad \mathbf{u, v} \in \mathcal{D} \tag{2.4}$$

暂时未指定其矩。

如果我们用类似于 $Z$ 的方式定义 $\boldsymbol{\varepsilon}$、$\mathbf{Y}$ 和 $\boldsymbol{ν}$，则 式（2.1） – 式（2.4） 表示了一个通用的线性混合模型，

$$\mathbf{Z} = \mathbf{T} \boldsymbol{\alpha} + \boldsymbol{\delta}, \qquad \boldsymbol{δ = ν + \varepsilon} \tag{2.5}$$

其中 $\mathbf{T}$ 是协变量 $(\mathbf{t}(\mathbf{s}_1),\ldots ,\mathbf{t}(\mathbf{s}_n))^{\prime}$ 的 $n \times p$ 矩阵。观察 式（2.5），误差项 $\boldsymbol{\delta}$ 由两个独立的零均值分量组成，导致 $\mathbb{E}(\boldsymbol{\delta}) = 0$ 和 $\operatorname{var}(\boldsymbol{\delta}) = \boldsymbol{\Sigma} \equiv (σ_{ij})$，其中

$$\sigma_{ij} = \begin{cases} C(\mathbf{s}_j,\mathbf{s}_j) + \sigma^2 v(\mathbf{s}_j), &i = j,\\ C(\mathbf{s}_i,\mathbf{s}_j), &i \neq j \end{cases}$$

如果记 $\mathbf{C} \equiv (C(\mathbf{s}_i,\mathbf{s}_j))$ 和 $\mathbf{V} \equiv \operatorname{diag}\{v(\mathbf{s}_1),\ldots,v(\mathbf{s}_n)\}$，则显而易见有：

$$\boldsymbol{\Sigma} = \mathbf{C} + \sigma^2 \mathbf{V} \tag{2.6}$$

注意：这里没有假设协方差函数的平稳性或各向同性。

我们感兴趣对 $Y$ 过程的推断，而不是含噪声的 $Z$ 过程。对于点预测，我们希望预测位置 $\mathbf{s}_0$ 处的 $Y$ 过程值（$\mathbf{s}_0 \in \mathcal{D}$），而不关心 $\mathbf{s}_0$ 是否是已观测位置。Cressie (1993 ^[7]) 的 第 3.4.5 节 给出了新位置点处过程值 $Y(\mathbf{s}_0)$ 的 克里金预测公式：

$$\hat{Y}(\mathbf{s}_0) = \mathbf{t}(\mathbf{s}_0)^{\prime} \hat{\boldsymbol{\alpha}} + \mathbf{k}(\mathbf{s}_0)^{\prime}(\mathbf{Z} − \mathbf{T} \hat{\boldsymbol{\alpha}}) \tag{2.7}$$

其中，回归系数的估计为：

$$\hat{\boldsymbol{\alpha}} = \mathbf{(T^{\prime} \boldsymbol{\Sigma}^{-1} T)^{-1} T^{\prime} \boldsymbol{\Sigma}^{-1}Z} \tag{2.8}$$

$$\mathbf{k(\mathbf{s}_0)^{\prime} = c(\mathbf{s}_0)^{\prime} \boldsymbol{\Sigma}^{-1}} \tag{2.9}$$

其中，$\mathbf{c}(\mathbf{s}_0) \equiv (C(\mathbf{s}_0, \mathbf{s}_1),\ldots,C(\mathbf{s}_0, \mathbf{s}_n))^{\prime}$ 为新位置 $\mathbf{s_0}$ 处的过程值变量与观测点处过程值变量之间的协方差向量。

式（2.7） 与克里金法的等价性可能不会立即显现出来，因为克里金法的传统推导是根据变异函数得出的，并且没有测量误差（即 式（2.1） 中的 $\varepsilon$ -过程等于 $0$）；参见 Journel 和 Huijbregts (1978 ^[22]) 的 第五章。 克里金预测标准差 为 $\hat{Y}(\mathbf{s}_0)$ 的均方根预测误差，即 $[\mathbb{E}\{Y(\mathbf{s}_0) − \hat{Y}(\mathbf{s}_0)\}^2]^{1/2}$，可以由下式解析给出：

$$\sigma_k(\mathbf{s}_0) = \{C(\mathbf{s}_0, \mathbf{s}_0) − \mathbf{k}(\mathbf{s}_0)^{\prime} \boldsymbol{\Sigma} \mathbf{k}(\mathbf{s}_0) + (\mathbf{t}(\mathbf{s}_0) − \mathbf{T}^{\prime} \mathbf{k}(\mathbf{s}_0))^{\prime} (\mathbf{T}^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{T})^{-1}(\mathbf{t}(\mathbf{s}_0) − \mathbf{T}^{\prime} \mathbf{k}(\mathbf{s}_0)) \}^{1/2} \tag{2.10}$$

当 式（2.7） 和 式（2.10） 中的预测位置 $\mathbf{s}_0$ 在域 $\mathcal{D}$ 中变化时，就会分别生成 kriging 预测图和 kriging 标准差图。 不过，在实际工作中，预测位置的数量通常有限，并且经常被取为叠加在 $\mathcal{D}$ 上的某个精细分辨率网格上的节点。

仔细观察 式（2.7） 和 式（2.10），可以发现协方差矩阵的求逆运算 $\boldsymbol{\Sigma}^{-1}$ 是一个重要组成部分，也是出现计算瓶颈最显著的地方。传统 $n \times n$ 对称正定矩阵的求逆计算成本为 $\mathcal{O}(n^3)$。当 $n$ 为数万及以上时， 式（2.7） 和 式（2.10） 通常无法在任何可接受的时间内完成计算。

在下一小节中，我们将展示如何通过选择协方差函数类，来加速最佳空间预测（即克里金法）。

2.2 空间协方差函数

一般而言， 式（2.4） 定义的协方差函数 $C(\mathbf{u,v})$ 在 $\mathbb{R}^d \times \mathbb{R}^d$ 上必须是正定的。 $C(\mathbf{u,v})$ 通常被建模为平稳的，此时 $C$ 是 $\mathbf{u-v}$ 的非负定函数。

在本文中，我们将尝试通过 $r$ 个（不一定正交的）基函数来捕获不同尺度的空间依赖性（见 式 2.11）：

$$\mathbf{S(\mathbf{u})} \equiv (S_1(\mathbf{u}),\ldots,S_r(\mathbf{u}))^{\prime}, \qquad \mathbf{u} \in \mathbb{R}^d \tag{2.11}$$

式中的 $r$ 值是固定的，基函数的具体形式我们将 在 3.1 节 中给出。如果令 $\mathbf{K}$ 是任意 $r \times r$ 大小的正定矩阵，则协方差 $\operatorname{cov} \{ Y(\mathbf{u}), Y(\mathbf{v}) \}$ 可以建模为：

$$C(\mathbf{u,v}) = \mathbf{S}(\mathbf{u})^{\prime} \mathbf{KS}(\mathbf{v}), \qquad \mathbf{u, v} \in \mathbb{R}^d \tag{2.12}$$

上式已被证明是一个非负定函数（第 3.1 节），因此是有效的协方差函数（Cressie 和 Johannesson，2006 年 ^[8] ）。当然，也可以进一步将噪声方差 $\tau^2 \mathbb{I}(\mathbf{u=v})$ 添加到 式（2.12） 中，不过我们在本文中没有这样做；见第 5 节。

很容易看出， 式 (2.12) 是将空间过程 $ν(\mathbf{s})$ 改写为 $ν(\mathbf{s})= \mathbf{S}(\mathbf{s})^{\prime} \boldsymbol{\eta}$ 的结果（ $\mathbf{s} \in \mathcal{D}$ ），其中 $\boldsymbol{\eta}$ 是一个的 $r$ 维的随机向量，其协方差矩阵为 $\operatorname{cov}(\boldsymbol{\eta}) = \mathbf{K}$。我们常常将 $ν(\cdot)$ 的模型称为 空间随机效应。 因此，式（2.3） 中的模型 $Y(\mathbf{s}) = \mathbf{t}(\mathbf{s})^{\prime} \boldsymbol{\beta} + \mathbf{S}(\mathbf{s})^{\prime} \boldsymbol{\eta}$，是一个混合效应的线性模型，通常被称为 空间混合效应模型。

2.3 固定秩克里金法

（1）训练点的协方差矩阵

由 式（2.12） 可得 $\mathbf{Y}$ (或 $\boldsymbol{ν}$) 的 $n \times n$ 的理论协方差矩阵为 $\mathbf{C=SKS^{\prime}}$，因此

$$\boldsymbol{\Sigma} = \mathbf{SKS^{\prime}} + \sigma^2 \mathbf{V} \tag{2.13}$$

其中未知参数是正定的 $r \times r$ 矩阵 $\mathbf{K}$，并且 $\sigma^2 > 0$。矩阵 $\mathbf{S}$ 的大小为 $n \times r$，其第 $(i，l)$ 个元素表示第 $i$ 个观测的第 $l$ 个基函数值 $S_l(\mathbf{s}_i)$，矩阵 $\mathbf{V}$ 是一个对角矩阵，其元素值由测量误差的方差给出。通常 $\mathbf{S}$ 和 $\mathbf{V}$ 都假设已知（注： $\mathbf{S}$ 已知表明基函数的形式已知）。

（2）测试点的协方差

进一步地，与新位置 $\mathbf{s_0}$ 相关的协方差矩阵为：

$$\operatorname{cov} \{Y(\mathbf{s}_0), \mathbf{Z}\} = \mathbf{c(s_0)^{\prime}} = \mathbf{S(\mathbf{s}_0)^{\prime}KS^{\prime}} \tag{2.14}$$

即在由 式（2.1） 、 式（2.3） 和 式（2.12） 指定的模型中，我们可以找到 式（2.7） 和 式（2.10） 中需要的所有分量的表达式。

（3）精度矩阵

当 $n$ 非常大时，仍然存在计算效率问题，但正如我们将要展示的那样， 式（2.12） 的协方差函数形式使得克里金预测方程中的计算仅涉及 $r \times r$ 矩阵的求逆。回顾 式（2.13），则协方差矩阵 $\boldsymbol{\Sigma}$ 的逆矩阵（即精度矩阵）为：

$$\boldsymbol{\Sigma}^{-1} = \sigma^{-1} \mathbf{V}^{-1/2} \{\mathbf{I} + (\sigma^{−1} \mathbf{V}^{-1/2} \mathbf{S}) \mathbf{K} (\sigma^{-1} \mathbf{V}^{-1/2}\mathbf{S})^{\prime}\}^{-1} \sigma^{-1} \mathbf{V}^{-1/2} \tag{2.15}$$

注意其中大括号内的项，具有 $\mathbf{I + PKP^{\prime}}$ 的形式，其中 $\mathbf{P}$ 是任意 $n \times r$ 矩阵。

让我们首先分析 $\mathbf{I + PKP^{\prime}}$ 的求逆问题，然后将分析的结果应用于 式 2.15。

根据 Sherman–Morrison–Woodbury 公式， $\mathbf{I + PKP^{\prime}}$ 可以重写为：

$$\mathbf{I + PKP^{\prime} = I + (I + PKP^{\prime})PK(I + P^{\prime}PK)^{−1}P^{\prime}}$$

在公式两边乘以 $\mathbf{(I + PKP^{\prime})^{-1}}$，可以得到该形式的通用逆矩阵解：

$$\mathbf{(I + PKP^{\prime})^{-1} = I − P(\mathbf{K}^{-1} + P^{\prime}P)^{-1}P^{\prime}}$$

（注：此处为 Sherman–Morrison–Woodbury 公式所涵盖的结果，参见 Henderson 和 Searle (1981 ^[16])）。

将上式用于 式（2.15），可以得到如下精度矩阵的简化计算公式：

$$\boldsymbol{\Sigma}^{-1} = (\sigma^2 \mathbf{V})^{-1} − (\sigma^2 \mathbf{V})^{-1} \mathbf{S} \{\mathbf{K}^{-1} + \mathbf{S}^{\prime}(\sigma^2 \mathbf{V})^{-1} \mathbf{S}\}^{-1} \mathbf{S}^{\prime}(\sigma^2 \mathbf{V})^{-1} \tag{2.16}$$

请注意， 式（2.16） 中原始协方差矩阵的求逆运算仅涉及 $r \times r$ 正定矩阵 $\mathbf{K}$ 和 $n \times n$ 对角矩阵 $\mathbf{V}$ 的求逆，如果 $r \ll n$ 的话，则会大大提升计算效率。

（4）空间预测

最后， 式（2.7） 的克里金法预测器（空间 BLUP）可以改写为：

$$\hat{Y}(\mathbf{s}_0) = \mathbf{t}(\mathbf{s}_0)^{\prime} \hat{\boldsymbol{α}} + \mathbf{S}(\mathbf{s}_0)^{\prime} \mathbf{KS}^{\prime} \boldsymbol{\Sigma}^{-1} (\mathbf{Z} − \mathbf{T} \hat{\boldsymbol{\alpha}}) \tag{2.17}$$

其中 $\hat{\boldsymbol{\alpha}} = (\mathbf{T}^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{T})^{-1} \mathbf{T}^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{Z}$ 和 $\boldsymbol{\Sigma}^{-1}$ 均可由 式（2.16） 计算给出。

（5）预测方差

式（2.10） 的克里金标准差可以改写为：

$$\begin{align*} \sigma_k(\mathbf{s}_0) = \{ &\mathbf{S}(\mathbf{s}_0)^{\prime} \mathbf{KS(s_0)} − \mathbf{S(s_0)}^{\prime} \mathbf{KS}^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{SK S(s_0)} \\ &+ ( \mathbf{\mathbf{t}(s_0) − T^{\prime} \boldsymbol{\Sigma}^{-1} SK S(\mathbf{s}_0)})^{\prime} ( \mathbf{T^{\prime} \boldsymbol{\Sigma}^{-1} T})^{-1} (\mathbf{\mathbf{t}(s_0) − T^{\prime} \boldsymbol{\Sigma}^{-1} SK S(s_0)}) \}^{1/2} \end{align*} \tag{2.18}$$

其中 $\boldsymbol{\Sigma}^{-1}$ 仍然由 式（2.16） 给出。

固定秩克里金法（Fixed Rank Kriging, FRK) 是 式（2.16） – 式（2.18） 方法的总称（Cressie 和 Johannesson，2006 ^[8]）。当 式（2.17） 和 式（2.18） 中的预测位置 $\mathbf{s}_0$ 随 $\mathcal{D}$ 变化时，可以分别生成克里金预测图和克里金标准差图。

（4）计算复杂度分析

仔细分析 式（2.16） – 式（2.18），对于固定的回归变量数量 $p$ 和 $\mathbf{K}$ 的固定秩 $r$，在由 式（2.12） 定义的协方差模型中，FRK 的计算负担仅关于 $n$ 呈线性，而不是三次方。为了看到这一点（为便于理解，假设 $\sigma^2 \mathbf{V} = \mathbf{I}$）：

• 与 式（2.17） 和 式（2.18） 相关的计算涉及 $\mathbf{S}^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{S}$、$\mathbf{S}^{\prime} \boldsymbol{\Sigma}^{-1} \mathbf{a}$ 和 $\boldsymbol{\Sigma}^{-1} \mathbf{a}$ ，其中向量 $\mathbf{a}$ 的长度为 $n$。
• 为了执行上述计算，可以预计算 $\mathbf{A} \equiv \mathbf{S}^{\prime} \mathbf{S}$ 和 $\mathbf{B} \equiv (\mathbf{K}^{-1} + \mathbf{S^{\prime}S})^{1/2}$，其最大计算成本为 $\mathcal{O}(nr^2)$（我们假设 $n > r$）。
• 由 式（2.16） 可知，$\mathbf{S^{\prime} \boldsymbol{\Sigma}^{-1} S = A − ABA}$，$\mathbf{ABA}$ 需要 $\mathcal{O}(r^3)$ 的计算成本。
• $\mathbf{S^{\prime} \boldsymbol{\Sigma}^{-1}a}$ 和 $\boldsymbol{\Sigma}^{-1} \mathbf{a}$ 的计算成本不会超过 $\mathcal{O}(nr^2)$。
• 式（2.17） 的预测方程计算成本为 $\mathcal{O}(r)$。
• 对于固定的 $\mathbf{s}_0$，式（2.18） 的预测方差计算成本为 $\mathcal{O}(r^2)$ （假设 $p \ll r$）。

因此，预测的总计算规模为 $\mathcal{O}(nr^2)$。

正如 第 4 节 所证实的那样，FRK 方法使得构建基于非常大的空间数据集的克里金预测变量图和克里金标准差图成为可能。

克里金方法和平滑方法之间的关系已经被很好地建立起来了（参见如 Cressie (1993 ^[7]) 的 第 5.9 节 和 Nychka (2000 ^[29])）。事实上，依赖于协方差模型 式（2.12） 的 FRK 是由固定秩平滑技术所驱动的，该技术从正则化和岭回归演化而来（Johannesson 和 Cressie，2004b ^[20]）。FRK 的新颖之处在于结合了 固定秩正定矩阵 $K$（待估计的参数）和 基函数 $\{S_l(\cdot)\}$（待指定），能够产生 式 2.12 形式的空间协方差函数，进而产生高计算效率的克里金预测和克里金标准差。

在下一节中，我们将详细考虑 $\mathbf{K}$ 和测量误差方差 $\sigma^2$ 的估计。

3 协方差函数的选择

在本文中考虑的协方差函数族，需要满足如下基函数表达形式：

$$C(\mathbf{u,v}) = \mathbf{S(\mathbf{u})^{\prime} K S(\mathbf{v})}, \qquad \mathbf{u, v} \in \mathbb{R}^d$$

其中 $\mathbf{K}$ 是一个 $r \times r$ 正定矩阵，$\mathbf{S}(\cdot)$ 是由基函数值 $S_1(\cdot),\ldots ,S_r(\cdot)$ 组成的 $r \times 1$ 向量，其中 $r$ 是固定值。最终结果是一个标量。

上述表示形式类似于 Stroud 等 (2001 ^[37] )，不过他们仅使用方法进行时空建模，但没有考虑为克里金预测求 $\boldsymbol{\Sigma}$ 的逆。

下面我们将给出此类协方差函数的一些性质。在经典地统计学意义上，我们还将展示数据是如何被二次使用（例如 Cressie (1989)，注：此处指参数估计和空间预测两个过程都需要观测数据的参与）：它们不仅（线性地）存在于克里金预测器（ 式（2.7））中；还可（非线性地）被用于估计空间协方差参数 $\mathbf{K}$ 和 $\sigma^2$。

3.1 基本性质

（1）协方差函数的非负定性

最重要的性质是：函数 $C(\mathbf{u,v})$ 是非负定的，其证明很简单：对于 $\mathbb{R}^d$ 中的任何位置 $\{\mathbf{s}_i : i = 1,\ldots ,m\}$，任何实数 $\{b_i : i = 1,\ldots ,m \}$ 和任何整数 $m$。因为 $\mathbf{K}$ 是正定的，所以有：

$$\sum^{m}_{i=1} \sum^{m}_{j=1} b_i b_j C(\mathbf{s}_i,\mathbf{s}_j) = \mathbf{b^{\prime}_m (S_m K S^{\prime}_m)b_m} = \mathbf{(S^{\prime}_m b_m)^{\prime} K (S^{\prime}_mb_m)} \geq 0$$

（2）与 KL 展开的关系

一个与 式（2.12） 相关但不同的模型是 Karhunen–Loéve 展开（例如 Adler (1981 ^[1] 的 第 3.3 节）。

注：

KL 展开： 将一个随机变量展开为多个随机系数与正交基函数的线性组合，即： $y(\mathbf{s}) = \sum^m_i \lambda_i \phi_i(\mathbf{s})$。其中，$\phi_i(\mathbf{s})$ 为相互正交的特征函数（或称基函数），$\lambda_i$ 为随机的非负特征值（或称基系数）。

根据 KL 展开，可以将随机过程的协方差函数改写为：

$$C_1(\mathbf{u,v}) \equiv \sum^{\infty}_{l=1} λ_l \phi_l(\mathbf{u}) \phi_l(\mathbf{v}) \tag{3.1}$$

其中 $\{λ_l\}$ 为非负特征值，$\{\phi_l(\cdot)\}$ 为正交特征函数。该式来源于如下积分方程：

$$\int C_1(\mathbf{u,v}) \phi(\mathbf{v}) d \mathbf{v} = λ \phi(\mathbf{u})$$

如果将 式（3.1） 在第 $k$ 项处截断，则可以得到随机过程的近似协方差函数：

$$C_2(\mathbf{u，v}) = \sum^{k}_{l=1} λ_l \phi_l(\mathbf{u}) \phi_l(\mathbf{v}) \equiv \boldsymbol{\phi(\mathbf{u})^{\prime} \Lambda \phi(\mathbf{v})}$$

上式最右侧为矩阵-向量表示，其中 $\boldsymbol{\Lambda}$ 是一个元素值非负的 $k \times k$ 对角矩阵。

不失一般性：

• 假设截断只保留具有正特征值的项；那么截断后的 Karhunen–Loéve 展开显然是 式（2.12） 的一个特例。

• 反之，如果将 式（2.12） 中的 $\mathbf{K}$ 改写成谱形式（或二次型形式），即 $\mathbf{K = P \boldsymbol{\Lambda} P^{\prime}}$，则我们可以看到 $C(\mathbf{u,v}) = \mathbf{(P^{\prime} S(u))^{\prime} \boldsymbol{\Lambda} (P^{\prime} S(v))}$，它看起来像一个截断的 Karhunen–Loéve 展开，只是采用了的基函数 $\{\phi_l(\cdot)\}$ 并非正交函数。

综上所述，实现 式（2.12） 的模型需要做如下工作：

• （1）选择一组有限的（可能非正交的）基函数 $\boldsymbol{\phi}(\mathbf{s})$，见 第 3.2 节。
• （2）估计一个固定秩的协方差矩阵 $\mathbf{K}$ （ 通常非单位矩阵 $\mathbf{I}$ ），见 第 3.3 节；

3.2 基函数的选择

（1）基函数的可选范围

固定秩的协方差矩阵 $\mathbf{K}$ 可以通过数据进行估计，但基函数 $\mathbf{S}(\cdot) \equiv (S_1(\cdot),\ldots ,S_n(\cdot))^{\prime}$ 通常必须提前指定。但由于不要求满足正交性约束，所以 $S_1(\cdot),\ldots ,S_r(\cdot)$ 的选择不受限制，可以包括平滑样条基函数（例如 Wahba (1990 ^[40])）、小波基函数（例如 Vidakovic (1999 ^[39])）和径向基函数（例如 Hastie 等 (2001 ^[15])，第 186-187 页）。

Nychka (2000 ^[29]) 和 Nychka 等 (2002 ^[31]) 汇集了各种类型的基函数，其中假设 $r = n$（即秩 $r$ 等于样本大小，因此可能不固定；或者如果 $r$ 固定，则 $\mathbf{K}$ 是对角的）。 Nychka 等 (2002 ^[31]) 令人信服地证明了 式（2.12） 协方差模型的固定秩基函数分解形式，能够逼近地统计学中使用的其他协方差函数，例如各向同性的指数模型。事实上，在第 4 节中，我们将在地球上进行克里金法，而选择的基函数为多分辨率局部双平方函数。

（2）建议考虑多分辨率

在基函数选择方面，我们的主要建议是要考虑多分辨率，以使 式（2.12） 的协方差模型能够捕获多个尺度的可变性。

在 第 2.2 节 中，我们曾提到 式（2.12） 可以被等效地认为是空间随机效应模型 $\nu(\cdot)=\mathbf{S}(\cdot)^{\prime} \boldsymbol{\eta}$，其中随机效应 $\boldsymbol{\eta}$ 具有由 $\mathbf{K}$ 指定的依赖结构（注：效应是在某些领域中对回归系数的特定称呼）。因此，如果 $\mathbf{S}(\cdot)$ 中包含多分辨率的分量，会捕获更多空间尺度的变化。

实际上， 式（2.3） 中的均值函数 $\mathbf{t}(\cdot)^{\prime} \boldsymbol{\alpha}$ 遗漏的一些较大尺度变化，有可能会被 $\mathbf{S}(\cdot)^{\prime} \boldsymbol{\eta}$ 的一些大尺度空间随机效应分量恢复。

最显著的多分辨率函数类别是各种类型的小波；同样，第 4 节 中使用的局部双平方函数类也是多分辨率的（但不正交）。事实上，在对气溶胶卫星数据的分析中，Shi 和 Cressie (2007 ^[35]) 选择了（非正交的）$W$ 小波作为均值函数向量 $\mathbf{t}(\cdot)$ 和基函数向量 $\mathbf{S}(\cdot)$ 的分量。此时的问题可以简化为在 $\mathbf{t}(\cdot)$ 和 $\mathbf{S}(\cdot)$ 选择使用哪些小波的问题，Shi 和 Cressie (2007 ^[35]) 给出了解决方案。

（3）小结

从多种基函数中选择哪一类基函数的问题尚在研究之中。但如果希望比较两个 FRK 图以检测变化异常，我们建议两者采用相同的 $\mathbf{t}(\cdot)$ 组分和 $\mathbf{S}(\cdot)$ 组分。

从计算角度来看，使用一类可以快速计算 $\mathbf{S}^{\prime} \mathbf{V}^{-1} \mathbf{S}$ 和 $\mathbf{S^{\prime}a}$ 的基函数是有益的（ $\mathbf{a}$ 为任意向量）。尽管我们在第 2.3 节中看到此类计算通常为 $\mathcal{O}(nr^2)$，但通过使用第 4 节中的双平方基或小波基以及稀疏矩阵库，计算成本实际上可以降低到 $\mathcal{O}(kr^2)$，其中 $k < n$。

3.3 协方差函数的拟合

本文中拟合空间协方差函数的策略与经典地统计学方法一致，如 Matheron (1963 ^[27]) 和 Journel 和 Huijbregts (1978 ^[22]) 的方法。其基本过程为：

• 首先，基于矩量法获得协方差矩阵的经验估计量 $\hat{\boldsymbol{\Sigma}}$，该估计量包含有噪声，而且可能不是正定的。
• 然后，选择一个能够保证正定的参数化矩阵类 $\{\boldsymbol{\Sigma(\theta) : \theta} \in \Theta\}$ ，通过优化方法估计能够使 $\hat{\boldsymbol{\Sigma}}$ 和 $\boldsymbol{\Sigma(\hat{\theta})}$ 最接近的 $\hat{\boldsymbol{\theta}} \in \Theta$，并将其作为最优参数值。
• 最后，将 $\boldsymbol{\Sigma(\hat{\theta})}$ 作为固定秩协方差矩阵的估计，代入克里金方程 式（2.17） 和 式（2.18） 。

注： 在 Zammit-Mangion 等（2021）中，提出了对 FRK 的改进和基于 EM 算法的估计。见 “Zammit-Mangion, Andrew, and Noel Cressie. “FRK : An R Package for Spatial and Spatio-Temporal Prediction with Large Datasets.” Journal of Statistical Software 98, no. 4 (2021). https://doi.org/10.18637/jss.v098.i04” 。

在 FRK 中，参数 $\boldsymbol{\theta}$ 包含两个部分：一是 $r \times r$ 的固定秩正定矩阵 $\mathbf{K}$ ，二是方差分量 $\sigma^2 \in (0, \infty )$。$\hat{\mathbf{K}}$ 和 $\hat{\sigma}^2$ 的估计可以将表示矩阵之间相似性的 Frobenius 范数作为目标函数，也就是说，通过最小化 “经验协方差矩阵” 和 “理论协方差矩阵” 之间 Frobenius 范数得到两部分参数的估计。

如果矩阵 $\mathbf{B}$ 是矩阵 $\mathbf{A}$ 的一个近似，那么误差矩阵就是 $\mathbf{B−A}$。最相似的矩阵应当使得误差矩阵有最小范数，而 Frobenius 范数就是一种最常用范数。该范数被定义为矩阵中各元素的平方和的二次方根，表达式为

$$\|A\|_{F}=\sqrt{\sum_{i=2}^{m} \sum_{j=1}^{n}\left|a_{i j}\right|^{2}}=\sqrt{\operatorname{trace}\left(A^{*} A\right)}=\sqrt{\sum_{i=1}^{\min \{m, n\}} \sigma_{i}^{2}}$$

可以将矩阵看做一个由各列向量拼接起来的长向量，而 Frobenius 范数正是该向量的欧式模。

（1） 去除趋势分量

如果要定义方差和协方差的经验估计，首先需要去除观测数据中的趋势分量（即均值函数部分）。由于缺乏空间依赖性的初始知识，同时考虑到计算效率，我们使用 $\boldsymbol{\alpha}$ 的普通最小二乘估计：

$$\boldsymbol{\bar{\alpha}} \equiv \mathbf{(T^{\prime} T)^{-1} T^{\prime} Z} \tag{3.2}$$

然后，计算细节残差：

$$D(\mathbf{s}_i) \equiv Z(\mathbf{s}_i) − \mathbf{t}(\mathbf{s}_i)^{\prime} \bar{\boldsymbol{\alpha}},\qquad i = 1,\ldots ,n \tag{3.3}$$

（2）协方差函数拟合

与经典地统计学一样，我们将数据 “装箱” 以计算空间依赖性的矩量法估计量。箱（bin）的数量 $M$ 是固定的且大于基函数的数量 $r$。因此，对于协方差函数的估计和拟合，一旦数据被装箱，计算复杂度就和 $n$ 无关了。

装箱后拟合的矩量法，常见于传统的克里金方法。其本质上就是拟合一条与箱均值最接近的协方差函数曲线。

假设 $\{\mathbf{u}_j : j = 1,\ldots ,M\}$ （ $r \leq M < n$ ）是一组位置（或箱中心位置），能够提供对 $\mathcal{D}$ 的良好覆盖。我们在第 $j$ 个箱中心 $\mathbf{u}_j$ 周围定义一个邻域（表示为 $N(\mathbf{u}_j)$ ），并定义如下 0-1 权重：

$$w_{ji} \equiv \begin{cases} 1, & \text{ if } \mathbf{s}_i \in N(\mathbf{u}_j), \\ 0, & \text{otherwise.} \end{cases} \tag{3.4}$$

其中 $i = 1,\ldots ,n$, $j = 1,\ldots ,M$。

令 $\mathbf{w}_j \equiv (w_{j1},\ldots ,w_{jn})^{\prime}$ ，并定义

$$\bar{Z_j} \equiv \mathbf{w}^{\prime}_j (\mathbf{Z} − \mathbf{T} \bar{\boldsymbol{\alpha}}) / \mathbf{w}^{\prime}_j \mathbf{1}_n \tag{3.5}$$

其中 $\mathbf{1}_n$ 是由 $1$ 构成的 $n \times 1$ 向量，是与箱中心 $\mathbf{u}_j, j = 1,\ldots ,M$ 相关联的平均去趋势数据。

在 附录 A 的近似 (A.4) 中，我们展示了 $\boldsymbol{\Sigma}_M \equiv \operatorname{var}（\bar{Z_1},\ldots ,\bar{Z_M})$ 可以近似为 $\bar{\boldsymbol{\Sigma}}_M(\mathbf{K}, \sigma^2) \equiv \bar{\mathbf{S}} \mathbf{K} \bar{\mathbf{S}}^{\prime} + \sigma^2 \bar{\mathbf{V}}$，其中 $\bar{\mathbf{S}}$ 和 $\bar{\mathbf{V}}$ 是 $\mathbf{S}$ 和 $\mathbf{V}$ 的易于计算的合并版本；与经典地统计学一样，近似是由分箱引起的小偏差引起的。此外，在附录 A 的表达式 (A.2) 中，我们给出了基于细节残差 式（3.3） 的经验正定估计 $\hat{\boldsymbol{\Sigma}}_M$。然后我们选择 $\mathbf{K}$ 正定和 $\sigma^2 \in (0, \infty )$ 使得 $\bar{\boldsymbol{\Sigma}}_M(\mathbf{K}, \sigma^2)$ 尽可能 “接近” $\hat{\boldsymbol{\Sigma}}_M$。

同阶的两个矩阵 $\mathbf{A}$ 和 $\mathbf{B}$ 之间存在各种矩阵范数。我们将使用的是 Frobenius 范数，

$$\|A − B\|^2 \equiv \operatorname{tr}\{(\mathbf{A − B})^{\prime}( \mathbf{A − B})\} = \sum_{j,k}(A_{jk} − B_{jk})^2 \tag{3.6}$$

这也被 Hastie (1996 ^[14]) 用于推导伪样条曲线，以及 Donoho 等 (1998 ^[10]) 估计协方差。本着与变异函数的加权最小二乘估计相同的精神 (Cressie, 1985 ^[4])，我们最终将使用加权 Frobenius 范数，但目前仅考虑 式（3.6） 的未加权版本。

当 $\sigma^2 = 0$ 时，$\hat{\boldsymbol{\Sigma}}_M(\mathbf{K},0) = \mathbf{ \bar{S} K \bar{S}^{\prime}}$。使用 Frobenius 范数，最小化 $\|\hat{\boldsymbol{\Sigma}}_M − \bar{\boldsymbol{\Sigma}}_M(\mathbf{K},0)\|$ 的 $\mathbf{K}$ 由下式给出（参见附录 A）

$$\mathbf{ \hat{K} = R^{−1} Q^{\prime} \hat{\boldsymbol{\Sigma}}_M Q (R^{−1})^{\prime}} \tag{3.7}$$

对应的拟合协方差矩阵为 $\bar{\boldsymbol{\Sigma}}_M ( \hat{\mathbf{K}} ,0) = \mathbf{QQ^{\prime} \hat{\boldsymbol{\Sigma}}_M QQ^{\prime}}$，其中 $\mathbf{\bar{S} = QR}$ 是 $\bar{\mathbf{S}}$ 的 Q–R 分解（即 $\mathbf{Q}$ 是 $M \times r$ 正交矩阵，$\mathbf{R}$ 是非奇异 $r \times r$ 上三角矩阵）。 Q-R 分解的计算成本为 $\mathcal{O}(r^3)$。由于 $\hat{\boldsymbol{\Sigma}}_M$ 是正定的，因此 $\hat{\mathbf{K}}$ 也是正定的。

当 $\sigma^2 \in (0, \infty )$ 时：

$$\|\hat{\boldsymbol{\Sigma}}_M − \bar{\boldsymbol{\Sigma}}_M(\mathbf{K}, \sigma^2)\| = \| \hat{\boldsymbol{\Sigma}}_M − \sigma^2 \bar{\mathbf{V}} − \bar{\mathbf{S}} \mathbf{K} \bar{\mathbf{S}}^{\prime}\|$$

产生的最佳参数估计（假设 $\sigma^2$ 为已知） 为：

$$\mathbf{ \hat{K} = R^{−1} Q^{\prime} (\hat{\boldsymbol{\Sigma}}_M − \sigma^2 \bar{\mathbf{V}}) Q(R^{−1})^{\prime}} \tag{3.8}$$

相应的拟合后协方差矩阵是 $\bar{\boldsymbol{\Sigma}}_M (\hat{\mathbf{K}}, \sigma^2 )$ 为：

$$\hat{\boldsymbol{\Sigma}}_M(\hat{\mathbf{K}}, \sigma^2 ) = \mathbf{QQ^{\prime}(\hat{\boldsymbol{\Sigma}}_M − \sigma^2 \bar{\mathbf{V}}) QQ^{\prime} + \sigma^2 \bar{\mathbf{V}}}\\ = \mathbf{QQ^{\prime} \hat{\boldsymbol{\Sigma}}_M QQ^{\prime} + \sigma^2 (\bar{\mathbf{V}} − QQ^{\prime} \bar{\mathbf{V}} QQ^{\prime})} \tag{3.9}$$

此时，可以通过关于 $\sigma^2 \in (0, \infty )$ 最小化下式，来获得 $\sigma^2$ 的估计：

$$\|\hat{\boldsymbol{\Sigma}}_M − \bar{\boldsymbol{\Sigma}}_M (\hat{\mathbf{K}}, \sigma^2)\|^2 = \sum_{j,k} \{( \hat{\boldsymbol{\Sigma}}_M − \mathbf{P}(\hat{\boldsymbol{\Sigma}}_M ))_{jk} − \sigma^2 (\bar{\mathbf{V}} − \mathbf{P} (\bar{\mathbf{V}}))_{jk} \}^2$$

其中对于任何 $M \times M$ 矩阵 $\mathbf{A}$，$\mathbf{P(A) \equiv QQ^{\prime} A QQ^{\prime}}$。计算成本为 $\mathcal{O}(M^3)$。请注意，这只是一个具有斜率 $\sigma^2$ 和零截距的简单线性回归。因此，最小化过程很容易实现。固定秩协方差矩阵的最终估计为：

$$\mathbf{K = R^{−1} Q^{\prime} (\hat{\boldsymbol{\Sigma}}_M − \hat{\sigma}^2 \bar{\mathbf{V}}) Q (R^{−1})^{\prime}} \tag{3.10}$$

对于 $r<M$，上述参数估计的计算成本为 $\mathcal{O}(M^3)$，这意味着 $M$ 的一个好的选择是独立于 $n$，因此不支配克里金法的计算成本 $\mathcal{O}(nr^2)$；参见第 4 节的时序数据比较。

最后，利用 $\hat{\mathbf{K}}$ 和 $\hat{\sigma}^2$ 的估计值，我们可以实施克里金预测；只需要将估计值代入 式（2.16） – 式（2.18） 即可。由此产生的 FRK 仅涉及 $r \times r$ 的固定秩矩阵和 $n \times n$ 的对角矩阵求逆，也就是说，我们实现了与数据量成线性关系的计算复杂度。

（3）使用加权范数作为目标函数

应当为变化较小或数据较多的箱赋予更多权重。因此考虑如下加权 Frobenius 范数：

$$\|\hat{\boldsymbol{\Sigma}}_M −\bar{\boldsymbol{\Sigma}}_M (\mathbf{K}, \sigma^2 )\|^2_a \equiv \sum_{j,k} a_ja_k \{(\hat{\boldsymbol{\Sigma}}_M)_{jk} − (\bar{\boldsymbol{\Sigma}}_M(K, \sigma^2 ))_{jk}\}^2 \tag{3.11}$$

其中 $a_1,\ldots ,a_M$ 是已知的正权重。等价地有，

$$\|\hat{\boldsymbol{\Sigma}}_M −\bar{\boldsymbol{\Sigma}}_M (\mathbf{K}, \sigma^2)\|^2_a = \|\bar{\mathbf{A}}^{1/2} \hat{\boldsymbol{\Sigma}}_M \bar{\mathbf{A}}^{1/2} − \bar{\mathbf{A}}^{1/2} \bar{\boldsymbol{\Sigma}}_M (\mathbf{K}, \sigma^2) \bar{\mathbf{A}}^{1/2}\|^2$$

其中 $\bar{\mathbf{A}} \equiv \operatorname{diag}(a_1,\ldots ,a_M)$，即 Frobenius 范数的加权版本仅涉及按 $\{a_j : j = 1,\ldots ,M\}$ 的比例缩放矩阵 $\hat{\boldsymbol{\Sigma}}_M$ 和 $\bar{\boldsymbol{\Sigma}}_M (\mathbf{K}, \sigma^2)$ 的行和列，因此在计算上并不比未加权版本更繁重。从附录 A 中的表达式 (A.5)，我们从统计上的选择是

$$a_j \propto (\mathbf{w}^{\prime}_j \mathbf{1}_n)^{1/2} / \mathbf{V}_D(\mathbf{u}_j), j = 1,\ldots ,M$$

这是一个基于数据的权重，其中 $\mathbf{V}_D(\mathbf{u}_j)$ 是第 $j$ 个箱中的经验方差，由 附录 A 中的表达式 (A.1) 给出。

总之，$\mathbf{K}$ 和 $\sigma^2$ 的估计问题基于最小化加权 Frobenius 范数。这是一个加权最小二乘准则，类似于 Cressie (1985 ^[4]) 给出的变异函数估计方法。它是基于矩的方法，而不是基于似然的方法。

在高斯假设下，$\mathbf{K}$ 和 $\sigma^2$ 的似然取决于 $\boldsymbol{\Sigma}^{-1}$ 和 $|\boldsymbol{\Sigma}|$。从 式（2.16） 的 Sherman–Morrison–Woodbury 公式 ，我们可以得到协方差矩阵的逆 $\boldsymbol{\Sigma}^{-1}$。

类似的公式可以得到原始的协方差矩阵行列式：

$$|\boldsymbol{\Sigma}|=|\sigma^2 \mathbf{V}| |K| |\mathbf{K}^{-1} + \mathbf{S}^{\prime}(\sigma^2 \mathbf{V})^{−1} \mathbf{S}|$$

可以看出，该计算仅涉及 $r \times r$ 矩阵的行列式。也就是说，计算 $\mathbf{K}$ 和 $\sigma^2$ 的似然是可行的。不过，最大化似然依然存在问题，需要进一步对 $\mathbf{K}$ 进行参数化 (Stein, 2008 ^[36])。Fuentes (2007 ^[11]) 在高斯性和协方差平稳性假设下，给出了大型空间数据集的近似似然。

4 臭氧卫星数据固定秩克里金法

5 讨论

本文介绍了大规模空间数据集的精确克里金法（空间 BLUP）方法。从计算成本来看，FRK 是线性可扩展的，可以处理海量数据集（千兆字节量级）。我们的结果依赖于 使用一类源自空间随机效应模型的非平稳协方差函数。回想一下 式（2.12）：

$$C(\mathbf{u,v}) \equiv \mathbf{S(\mathbf{u})^{\prime} K S(\mathbf{v})}, \qquad \mathbf{u, v} \in \mathbb{R}^d,$$

其中 $\mathbf{S}(\cdot) \equiv (S_1(\cdot),\ldots ,S_r(\cdot))^{\prime}$ 是基函数向量。 $r \times r$ 的正定矩阵 $\mathbf{K}$ 是空间相关参数，我们使用加权最小二乘法以经典地统计学方式对其进行估计； $\mathbf{K}$ 的最大似然估计是未来研究的主题。正在开发中的贝叶斯方法在 $\mathbf{K}$ 上放置先验（例如 Wishart 分布），而后贝叶斯模型定义可以基于该假设完成，例如，$Y(\cdot)$ 是独立于高斯分布的高斯过程白噪声过程 $\varepsilon(\cdot)$，$\sigma^2$ 的先验是伽玛分布。

贝叶斯分析还允许在非线性地统计模型（Diggle 等 (1998 ^[9])）中进行最佳空间预测。 尽管此类模型对于大规模空间数据集的计算量可能很大，但由于在分析中使用了马尔可夫链蒙特卡罗计算，仍然有机会通过使用 式（2.12） 来实现计算加速。作为这方面的证据，Hrafnkelsson 和 Cressie (2003 ^[17]) 将标准地统计协方差模型与直接对逆协方差矩阵建模的模型进行了比较，他们报告说后者导致计算效率提高了 $5$ 倍以上。

（隐）$Y$ 过程中的微尺度变化可以通过在 式（2.6） 中包含另一个对角矩阵来建模。当两个对角矩阵彼此成比例时，测量误差参数 $\sigma^2$ 和微尺度参数 $\tau^2$ 无法单独识别，尽管两者之和 $\tau^2 + \sigma^2$ 可以识别。该和在地统计学文献中通常被称为 “块金效应”。在本文中，我们假设 $Y(\cdot)$ 是平滑的（即 $\tau^2 = 0$），其余的可变性则是由测量误差引起的。

由 式（2.12） 给出的非平稳协方差函数具有显著的支持度变化特性。令 $B \subset \mathbb{R}^d$ 并定义 $Y(B) \equiv \int_B Y(\mathbf{s}) d \mathbf{s}=|B|$，其中 $|B|$ 是 $B$ 的 $d$ 维体积。那么

$$\operatorname{cov}\{Y(B_1), Y(B_2)\} = \mathbf{S}(B_1)^{\prime} \mathbf{KS}(B_2), B_1, B_2 \subset \mathbb{R}^d,$$

其中 $\mathbf{S}(B) \equiv (S_1(B),\ldots,S_r(B))^{\prime}$ 和 $S_l(B) \equiv \int_B S_l(\mathbf{s}) d \mathbf{s}= |B|$，其中 $B \subset \mathbb{R}^d$。

因此，无论数据和预测变量的支撑如何，克里金方程都采用与 式（2.16） – （2.18） 相同的形式。实际上，基本功能将离线集成。

最后，从第 3.2 节中给出的空间随机效应模型 $ν(\mathbf{s}) = \mathbf{S}(\mathbf{s})^{\prime} \boldsymbol{\eta}$ 可以自然推广到时空随机效应模型 $ν(\mathbf{s}, t) = \mathbf{S}(\mathbf{s})^{\prime} \boldsymbol{\eta}(t)$，其中 $\{\boldsymbol{\eta}(t) : t = 0, 1, 2,\ldots \}$ 是均值为 $0$ 且 $\operatorname{cov} \{\boldsymbol{\eta}(t_1), \boldsymbol{\eta}(t_2)\} \equiv \mathbf{K} (t_1, t_2)$, $t_1, t_2 = 0, 1, 2,\ldots$。 基于此模型的时空克里金法和时空卡尔曼滤波，以及涉及趋势的混合模型版本，目前正在研究中。

参考文献

• [1] Adler, R. J. (1981) The Geometry of Random Fields. Chichester: Wiley.
• [2] Billings, S. D., Beatson, R. K. and Newsam, G. N. (2002a) Interpolation of geophysical data using continuous global surfaces. Geophysics, 67, 1810–1822.
• [3] Billings, S. D., Newsam, G. N. and Beatson, R. K. (2002b) Smooth fitting of geophysical data using continuous global surfaces. Geophysics, 67, 1823–1834.
• [4] Cressie, N. (1985) Fitting variogram models by weighted least squares. J. Int. Ass. Math. Geol., 17, 563–586.
• [5] Cressie, N. (1989) Geostatistics. Am. Statistn, 43, 197–202.
• [6] Cressie, N. (1990) The origins of kriging. Math. Geol., 22, 239–252.
• [7] Cressie, N. (1993) Statistics for Spatial Data, revised edn. New York: Wiley.
• [8] Cressie, N. and Johannesson, G. (2006) Spatial prediction of massive datasets. In Proc. Australian Academy of Science Elizabeth and Frederick White Conf ., pp. 1–11. Canberra: Australian Academy of Science.
• [9] Diggle, P. J., Tawn, J. A. and Moyeed, R. A. (1998) Model-based geostatistics. Appl. Statist., 47, 299–326.
• [10] Donoho, D. L., Mallet, S. and von Sachs, R. (1998) Estimating covariances of locally stationary processes: rates of convergence of best basis methods. Technical Report 517. Stanford University, Stanford.
• [11] Fuentes, M. (2007) Approximate likelihoods for large irregularly spaced spatial data. J. Am. Statist. Ass., 102, 321–331.
• [12] Furrer, R., Genton, M. G. and Nychka, D. (2006) Covariance tapering for interpolation of large spatial datasets. J. Computnl Graph. Statist., 15, 502–523.
• [13] Haas, T. C. (1995) Local prediction of a spatio-temporal process with an application to wet sulfate deposition. J. Am. Statist. Ass., 90, 1189–1199.
• [14] Hastie, T. (1996) Pseudosplines. J. R. Statist. Soc. B, 58, 379–396.
• [15] Hastie, T., Tibshirani, R. and Friedman, J. (2001) Elements of Statistical Learning: Data Mining, Inference, and Prediction. New York: Springer.
• [16] Henderson, H. V. and Searle, S. R. (1981) On deriving the inverse of a sum of matrices. SIAM Rev., 23, 53–60.
• [17] Hrafnkelsson, B. and Cressie, N. (2003) Hierarchical modeling of count data with application to nuclear fall-out. Environ. Ecol. Statist., 10, 179–200.
• [18] Huang, H.-C., Cressie, N. and Gabrosek, J. (2002) Fast, resolution-consistent spatial prediction of global processes from satellite data. J. Computnl Graph. Statist., 11, 63–88.
• [19] Johannesson, G. and Cressie, N. (2004a) Variance-covariance modeling and estimation for multi-resolution spatial models. In geoENV IV—Geostatistics for Environmental Applications (eds X. Sanchez-Vila, J. Carrera and J. Gomez-Hernandez), pp. 319–330. Dordrecht: Kluwer.
• [20] Johannesson, G. and Cressie, N. (2004b) Finding large-scale spatial trends in massive, global, environmental datasets. Environmetrics, 15, 1–44.
• [21] Johannesson, G., Cressie, N. and Huang, H.-C. (2007) Dynamic multi-resolution spatial models. Environ. Ecol. Statist., 14, 5–25.
• [22] Journel, A. G. and Huijbregts, C. (1978) Mining Geostatistics. London: Academic Press.
• [23] Kammann, E. E. and Wand, M. P. (2003) Geoadditive models. Appl. Statist., 52, 1–18.
• [24] London, J. (1985) The observed distribution of atmospheric ozone and its variations. In Ozone in the Free Atmosphere (eds R. C. Whitten and S. S. Prasad), pp. 11–80. New York: Van Nostrand Reinhold.
• [25] Madrid, C. R. (1978) The Nimbus-7 User’s Guide. Greenbelt: NASA.
• [26] Matheron, G. (1962) Traite de Geostatistique Appliqueé, vol. I. Paris: Technip.
• [27] Matheron, G. (1963) Principles of geostatistics. Econ. Geol., 58, 1246–1266.
• [28] McPeters, R. D., Bhartia, P. K., Krueger, A. J., Herman, J. R., Schlesinger, B. M., Wellemeyer, C. G., Seftor, C. J., Jaross, G., Taylor, S. L., Swissler, T., Torres, O., Labow, G., Byerly, W. and Cebula, R. P. (1996) The Nimbus-7 Total Ozone Mapping Spectrometer (TOMS) Data Products User’s Guide. Greenbelt: NASA.
• [29] Nychka, D. (2000) Spatial-process estimates as smoothers. In Smoothing and Regression: Approaches, Computation, and Application (ed. M. G. Schimek), pp. 393–424. New York: Wiley.
• [30] Nychka, D., Bailey, B., Ellner, S., Haaland, P. and O’Connell, M. (1996) FUNFITS: Data Analysis and Statistical Tools for Estimating Functions. Raleigh: North Carolina State University.
• [31] Nychka, D., Wikle, C. and Royle, J. A. (2002) Multiresolution models for nonstationary spatial covariance functions. Statist. Modllng, 2, 315–331.
• [32] Quiñonero-Candela, J. and Rasmussen, C. E. (2005) A unifying view of sparse approximate Gaussian process regression. J. Mach. Learn. Res., 6, 1939–1959.
• [33] Rue, H. and Tjelmeland, H. (2002) Fitting Gaussian Markov random fields to Gaussian fields. Scand. J. Statist., 29, 31–49.
• [34] Sahr, K. (2001) DGGRID Version 3.1b: User Documentation for Discrete Global Grid Generation Software. Ashland: Southern Oregon University. (Available from http://www.sou.edu/cs/sahr/dgg/.)
• [35] Shi, T. and Cressie, N. (2007) Global statistical analysis of MISR aerosol data: a massive data product from NASA’s Terra satellite. Environmetrics, 18, 665–680.
• [36] Stein, M. (2008) A modeling approach for large spatial datasets. J. Kor. Statist. Soc., 37, in the press.
• [37] Stroud, J. R., Müller, P. and Sansó, B. (2001) Dynamic models for spatiotemporal data. J. R. Statist. Soc. B, 63, 673–689.
• [38] Tzeng, S., Huang, H.-C. and Cressie, N. (2005) A fast, optimal spatial-prediction method for massive datasets. J. Am. Statist. Ass., 100, 1343–1357.
• [39] Vidakovic, B. (1999) Statistical Modeling by Wavelets. New York: Wiley.
• [40] Wahba, G. (1990) Spline Models for Observational Data. Philadelphia: Society for Industrial and Applied Mathematics.