🔥 大规模空间表面时间序列建模

【摘 要】 在许多现象中观测到的数据都具有空间和时间成分。由于复杂高性能技术的快速发展，现在可以大规模收集时空数据。然而，大型时空数据集的统计建模涉及几个具有挑战性的问题。例如，处理大型数据集和时空非平稳性在计算上具有挑战性。因此，有必要开发新的统计模型。在这里，我们提出了一种新方法来模拟复杂的大型时空数据集。在我们的方法中，在每个时间点估计一个连续的表面，用于捕获空间依赖性（可能是非平稳的）。以这种方式，时空数据产生一系列表面。然后，使用函数型时间序列技术对此表面序列进行建模。函数型时间序列方法使我们能够获得计算上可行的方法，并且还在时间预测方面提供了广泛的灵活性。我们通过蒙特卡罗模拟研究来说明这些优势。我们还使用超过 400 万个值的高分辨率风速模拟数据集测试了方法的性能。总的来说，本方法使用了一种新的数据分析范式，其中随机场被视为一个单一的实体，这在大数据的背景下是一种非常有价值的方法。

【原 文】 I. Martínez-Hernández and M. G. Genton, “Surface time series models for large spatio-temporal datasets”, Spatial Statistics, 2022, doi: https://doi.org/10.1016/j.spasta.2022.100718.

1 简介

在我们所处的时代，数据在许多情况下都是重要资产，包括决策制定。由于复杂、高性能技术的快速发展，现在可以大规模收集数据，从而产生高维和高频数据。统计方法有望适用于这些复杂的大型数据集。然而，复杂和大型数据集的统计分析涉及计算挑战，有时需要高性能计算，这可能是从业者的限制因素。这些挑战将我们引向了一种新的数据分析范式。克服这些挑战的一种方法是假设观测结果具有沿连续体变化的特征，例如曲线或曲面。也就是说，曲线或曲面被视为单点观测。在这里，我们使用这种连续的方法提出了一种方法来建模和预测空间密集且时间密集的复杂和大型时空数据。

人们对时空建模的兴趣越来越大。传统上，时空数据的统计方法侧重于对时空协方差函数 $C(\mathbf{s}_1, t_1, \mathbf{s}_2, t_2)$ 建模，其中 $(\mathbf{s}_i, t_i)，i = 1, 2$ 是两个时空位置。因此，研究人员投入了大量工作和精力来研究协方差函数的有效模型（例如，参见 Chen 等，2021 年）^[4]。一个简单且广泛使用的时空协方差函数是纯空间和纯时间协方差函数的乘积，$C(\mathbf{s}_1, t_1, \mathbf{s}_2, t_2) = CS(\mathbf{s}_1, \mathbf{s}_2) \cdot CT(t_1, t_2)$，其中 $CS$ 和 $CT$ 分别是纯空间协方差函数和纯时间协方差函数。这样的时空协方差函数被称为可分离的。尽管可分离性不允许空间和时间之间的相互作用，但它通常很有吸引力，因为它在某些情况下会导致计算上可行的方法（Genton，2007）^[11]。然而，对时空协方差函数建模的方法对于大型数据集可能不可行。例如，分解密集的 $N × N$（协方差）矩阵的计算成本为 $\mathcal{O}(N^3)$，这对于大 $N$ 具有挑战性。

（1）描述性方法

一些克服计算成本的方法是协方差锥化（Furrer 等，2006 年^[8]；Kaufman 等，2008 年^[20]）和固定秩克里金法（Cressie 和 Johannesson，2008 年^[5]；Nguyen 等，2014 年^[31]）；参见 Sun 等的评论 (2012)^[41]，Heaton 等（2019）^[15]，和 Huang 等 (2021)^[17]。协方差锥化对协方差函数使用紧支撑近似，固定秩克里金法旨在降低参数维数。对于非平稳时空数据，请参阅 Nychka 等 (2018)^[32] 以及 Kuusela 和 Stein (2018)^[22]。

（2）动态模型方法

时空建模的另一种方法是将时间序列和空间统计相结合的技术。这种方法提供了空间描述性和时间动态的统计模型（Wikle 和 Cressie，1999 年 ^[43]；Gelfand 等，2005 年^[10]；Wikle 和 Hooten，2010 年^[44]；Sigrist 等，2012 年^[38]）。时间建模主要考虑时空变化，换句话说，时空数据被认为是一个多元时间序列，其组分在空间上是相关的。按照此想法，可以使用许多种类型的时间序列模型，但可能并不适用于空间上比较丰富的数据（大数据或密集数据）。如果观测数据在空间上是密集的，那么，在每个时间点数据都可以被认为是一个连续的表面（随机场的连续实现），我们按照这个思路，使用时间序列方法对表面进行建模。在面临大 $N$ 问题时，这种连续性方法可以通过在有限维子空间（降维）中逼近连续随机场，来获得计算优势。

（3）函数型数据分析方法

在本文中，我们使用动态模型方法，具体来说，结合了 函数型数据分析 (Functional Data Analysis, FDA) 和空间统计的技术。

函数型数据分析假设观测具有沿连续体变化的特征，例如曲线或曲面（Ramsay 和 Silverman，2005 年^[34]；Kokoszka 和 Reimherr，2017 年^[21]）。因此，函数型数据分析适合于处理本质上在无限维空间中定义的数据。大部分情况下，函数型数据分析关注曲线数据，不过鉴于其对大型复杂数据建模的能力，函数型数据分析方法也已经扩展到了时空数据。一种扩展假设在每个空间位置都观测到了一条曲线（Giraldo 等，2011 年^[12]；Menafoglio 等，2016 年^[30]）。另一个扩展假设空间数据是确定性连续随机场的实现（Sangalli 等，2013 年^[37]；Bernardi 等，2017 年^[2]）。在本文中，我们使用与后一种方法类似的想法。

我们结合了函数型时间序列和连续随机场模型的技术。我们假设每一个时间点的观测是随机场的一次实现（称为表面数据），因此时空数据产生一个表面序列，我们称其为表面时间序列。通过此方法， 空间依赖性可以通过表面的连续性来捕获，而时间依赖性则通过时间序列方法来捕获（例如，参见 Martnez-Hernandez 和 Genton ，2020 年^[26]）。在实际工作中，需要估计连续型表面数据。这种连续型表面的估计通常与经典空间数据分析中的克里金法和高斯过程有关。但在本文中，我们使用了一种类似于一维函数型数据的处理方法来估计连续型表面；具体来说，使用有限基函数来表示连续随机场。

通过有限基函数表示法估计连续曲面是一种常用方法。基函数的一个例子可以是单变量 B 样条的张量积 (Wood, 2006 ^[46]; Xiao 等, 2013 ^[47])。基函数模型的估计可以表示为带有惩罚项的平方误差和的最小化问题，也被称为 基惩罚平滑器（Wood，2017 ^[45]）。基函数的另一种选择是有限元法，例如，参见 Sangalli 等 (2013) ^[37]。

另外一种对基函数模型的估计方法是采用 随机偏微分方程 (SPDE)。在 Lindgren 等 (2011) ^[25]中，研究了高斯随机场与特定随机偏微分方程之间的联系。在那篇论文中，高斯马尔可夫随机场 (GMRF) 被用来近似随机偏微分方程的解。我们在 第 2.1 节 中详细介绍了这种近似。它的优点是：高斯马尔可夫随机场的精度矩阵具有稀疏结构，拥有计算优势，尤其是在大型数据集中。 在本文中，我们使用随机偏微分方程方法和有限元法来估计表面的连续表示。请注意，原则上，可以使用上述任何方法来估计连续表面，因为我们的方法不限于随机偏微分方程方法的特定情况。

（4）本文目标

我们在本文中的目标是：提出一种新方法来分析和预测复杂的大型时空数据，这种方法在计算上是可行的，并且可以使用现有的 R 包实现（R Core Team，2022 ^[33]）。

我们假设数据是在连续空间中取值的时间序列的实现。具体来说：

• 对于每个离散时间，估计其连续随机场。因此，在随机场的连续估计中考虑了（非平稳）空间相关性。
• 假设数据是表面时间序列数据，对时间依赖性进行建模。
• 使用时间序列技术来描述时间依赖性，然后预测时空数据。

本文方法的一个 优点是：它可以处理具有不同数量空间位置集的时空数据，以及在每个时间 $t$ 具有不同位置集的时空数据，因为我们会估计每个离散时间的连续随机场。

本文方法得到高分辨率模拟风数据（来自天气研究和预报 (WRF) 模型）的推动。图 1 说明了我们的想法。在左侧子图中，显示了沙特阿拉伯上空 $\mathbf{s}_i$ 中的 $1500$ 个位置，假设在这些位置观测了风速数据。令 图 1（中） 的 $y_{t}(\mathbf{s}_i)$ 是在固定时间 $t$ 和位置 $\mathbf{s}_i$ 处的观测数据，$t = 1, \ldots , T$。对于每个时间 $t$ 的观测数据 $y_{t}(\mathbf{s}_i)$，我们基于有限基函数 $\{\psi_k\}$ 可以估计一个连续随机场 $x_{t}(\mathbf{s}) = \sum_k b_{k,t} \psi_k(\mathbf{s})$，见 图 1（右）。

然后，我们使用多元函数型时间序列方法，对基函数的系数 $\{b_{k,t}\}$ 进行建模以反映时间依赖性。此环节使用了历史信息和每个的近邻信息，这在第 2.2 节中有详细说明。因此，我们通过系数实现在空间和时间之间进行交互。我们的方法在计算效率上比协方差建模方法更可行。此外，我们的方法可以对空间位置不时变化的数据进行建模，因为随机场 $x_{t}$ 对于每个时间点 $t$ 都被认为是连续的。本文为从业者提供了一种无需高性能计算即可对大型时空数据进行建模的替代方法。

图 1： （左）沙特阿拉伯观测风速的测点位置。我们总共有 $1500$ 个地点。（中） $t$ 时间风速测量值平方根的快照。（右）根据基函数估计的连续随机场 $x_t(\mathbf{s})$。

我们论文的其余部分组织如下：

• 在第 2 节中，介绍了我们提出的方法。本节介绍连续表面估计和表面时间序列建模。
• 在第 3 节中，描述了如何估计方法中涉及的参数。此外，还描述了如何获得预测。
• 在第 4 节中，进行了模拟研究，以评估方法在不同模拟设置下的性能。
• 在第 5 节中，分析了来自天气研究和预报 (WRF) 模型的每小时风速数据集。
• 在第 6 节中提出了一些讨论。

2 表面时间序列模型

我们考虑对复杂的大型时空数据集进行分析，其中 假设空间成分是连续的，而时间成分是离散的。具体来说，我们假设观测到的数据为 $\{y_{1}(\mathbf{s}), y_{2}(\mathbf{s}), . . . , y_{T} (\mathbf{s})\}$，其中对于每个 $t \in \{1, \ldots,T \}$ 来说，$y_{t}: D \subset \mathbb{R}^2 \rightarrow \mathbb{R}$ 都是一个连续函数，$\mathbf{s} \in D$ 代表位置，$D$ 代表研究区域。在实际工作中，数据是在 $D$ 中的一组有限点上观测到的，即对于每个 $t = 1，\ldots, T$, 我们观测到 $m_t$ 个点 $y_{t} = \{y_{t}(\mathbf{s}_1), \ldots, y_{t}(\mathbf{s}_{m_t})\}^T$。

此外，观测到的数据可能存在某种测量误差。因此，假设观测到的数据可以表示为：

$$\begin{align*} y_{t}(\mathbf{s}_i) &= x_{t}(\mathbf{s}_i) + \varepsilon_{t}(\mathbf{s}_i) \\ x_{t}(\mathbf{s}_i) &= \mu_{t}(\mathbf{s}_i) + z_{t}(\mathbf{s}_i) \tag{1} \end{align*}$$

其中 $\mu_{t}(\mathbf{s}_i)$ 是反映宏观尺度变化的均值，$\boldsymbol{\varepsilon}_{t} = \{\varepsilon_{t}(\mathbf{s}_1), \ldots , \varepsilon_{t}(\mathbf{s}_{m_t}) \}^T \sim N_{m_t} (\mathbf{0}, τ^{ −1}_{ε,t} \mathbf{I}_{m_t})$ 是反映测量误差的随机噪声（块金）并假设独立同分布。我们假设 $z_{t}(\mathbf{s}_i)$ 是未知连续高斯随机场的实现，这是对空间（和空间-时间）依赖性建模的最常见选择（Cressie 和 Wikle，2011 ^[6]）。

研究目标是：对连续表面的时间序列 $\{x_{t}; t = 1, \ldots, T\}$ 建模并实现预报。为此，我们使用在函数型时间序列中广泛使用的两步法（参见 Hyndman 和 Ullah (2007)^[19] 以及 Aue 等 (2015)^[1] ）：

• 第一步：估计连续函数 $x_{t}$；
• 第二步：对表面时间序列 $\{x_{t}; t = 1, \ldots , T\}$ 建模。

下面我们详细介绍这两个步骤。

2.1 估计连续随机场

与普通函数型时间序列建模连续函数稍有不同，本文第一步是估计连续随机场。为了简化符号，假设 $t$ 时间的均值是常数，即 $\mu_{t}(\mathbf{s}_i) = \mu_{t}$，因此，我们重点关注随机场 $z_{t}$。

对于每个 $t=1,\ldots,T$ ，连续随机场 $z_{t}(\mathbf{s})$ 解释了空间依赖性。因此，我们使用 Matern 协方差函数估计 $z_{t}(\mathbf{s})$ （参见 Matern，1986 年^[29]；Stein，1999 年^[40] ），对于两个位置 $\mathbf{s}_1$ 和 $\mathbf{s}_2$， 该函数定义为：

$$\operatorname{Cov} \{z_{t}(\mathbf{s}_1), z_{t}(\mathbf{s}_2)\} = \frac{\sigma^2_{z,t}}{ \Gamma(\nu_{t})2^{\nu_{t}−1}} (\kappa_{t} \| \mathbf{s}_1 − \mathbf{s}_2 \| )^{\nu_{t}}K_{\nu_{t}}(\kappa_{t} \| \mathbf{s}_1 − \mathbf{s}_2 \| )$$

其中 $\| \cdot \|$ 为欧氏距离，$K_{\nu_{t}}$ 为第二类 $\nu_{t}$ 阶修正贝塞尔函数，$\sigma^2_{z,t}=1/τ_{z,t}$ 为边缘方差。参数 $\nu_{t}$ 控制底层过程的均方可微性。参数 $\kappa_{t} > 0$ 是与变程 $\rho_t$ 相关的缩放参数。根据经验得出的变程定义为 $\rho_t = \sqrt{8\nu_{t}}/\kappa_{t}$，其中 $\rho_t$ 对应于 $\nu_{t}$ 下空间相关性接近 $0.1$ 的距离。在实际工作中，平滑度参数 $\nu_{t}$ 通常基于我们对底层过程平滑度的先验信念，是固定的。 在这里，我们为所有 $t$ 固定 $\nu_{t} = 1$，类似于 Lindgren 等 (2011) ^[25] 和 Cameletti 等 (2013) ^[3]。 $\nu_{t}$ 的选择足够灵活，可以涵盖一大类空间变化（参见 Rue 等，2009 年^[36]）。

从 式（1） 的模型中，我们有 $\mathbf{y}_{t} | \mathbf{x}_{t}, \tau_{\varepsilon,t} \sim \mathcal{N}_{m_t}(\mathbf{x}_{t}, \tau^{-1}_{\varepsilon,t} \mathbf{I}_{m_t})$，其中 $\mathbf{x}_{t} = \{x_{t}(\mathbf{s}_1), \ldots , x_{t}(\mathbf{s}_{m_t})\}^T$。通常，高斯随机场的推断使用 $\mathbf{z}_{t}$ 的协方差矩阵 $\boldsymbol{\Sigma}_t$ 获得。但这对于大型数据集来说，在计算上是不可行的，因为协方差矩阵过于密集、难以处理。为了降低计算成本，我们注意到 具有 Matern 协方差的高斯空间过程是如下随机偏微分方程的解:

$$\tau_{t} (\kappa_{t}^2 − \Delta)^{\alpha_{t}/2} z_{t}(\mathbf{s}) = W_{t}(\mathbf{s}) \tag{2}$$

其中 $\Delta = \sum^2_{j=1} \partial^2/\partial s_j^2$ 为拉普拉斯算子，$\alpha_{t} = \nu_{t} + 1$ 控制空间过程的平滑度，$\tau_{t} > 0$ 控制方差。过程 $W_{t}$ 是空间高斯白噪声。

对于上述随机偏微分方程，Lindgren 等 (2011)^[25] 使用有限元方法构造了最接近 式 (2) 解的高斯马尔可夫随机场。给定域 $D$ 的三角剖分，式 (2) 的解表示为：

$$z_{t}(\mathbf{s}) = \sum^K_{k=1} b_{k,t} \psi_k(\mathbf{s}) \tag{3}$$

其中 $\{\psi_k\}$ 是在三角网格上定义的分段线性基函数，$\{b_{k,t}\}$ 是高斯分布的随机权重。式 (3) 基于分布原则中的等价原理得到，从而产生高斯马尔可夫随机场。$K$ 的选择取决于应用，由用户指定；有关详细信息，请参阅 Lindgren 等 (2011) ^[25]。因此，我们可以用 式 (3) 的有限元表示法来近似 式(1) 中的连续随机场 $z_{t}$。

在下文中，我们假设数据是一系列具有如 式 (3) 中表示的连续表面，暂时不考虑这些表面是如何估计的。

2.2 时间依赖性建模

第二步是在考虑空间依赖性的情况下，对时间依赖性进行建模。

在本节中，在不失一般性的情况下，假设观测值是随时间观测到的具如 式 (3) 形式的连续随机场（即均值分量 $\mu_{t}$ 为零）。因此，$\mathbf{b}_t = (b_{1,t}, \ldots , b_{K,t})^T$ 是一个 $K$ 维系数向量，用作空间相关基函数的时间相关权重（此处意味着 $\mathbf{b}$ 本身就是一个时间序列，可以用时间序列模型建模）。时间序列 $\{\mathbf{b}_t\}_t$ 可以产生高维和非平稳时间序列。$K$ 将取决于区域 $D$ 的网格有多精细。在我们的数据分析（第 5 节）中，使用了 $K = 364$，$D$ 的相应网格如 图 3 所示。在本文中，我们假设时间序列 $\{\mathbf{b}_t\}_t$ 是高维的，并且区域 $D$ 可能很大或具有很复杂的形状。

我们使用多元函数型时间序列模型对 $\mathbf{b}_t$ 进行建模。为此，对于每个 $k = 1,\ldots,K$，假设 $\{b_{k,t}\}_t$ 是函数型时间序列 $\{f^{(k)}_n(u)； u \in \mathcal{T} , n = 1, \ldots, N \}$ 的实现, 其中:

$$f^{(k)}_n (u) = b_{k,u} \mathbf{1} \{u \in [(n − 1) \delta, n \delta], n \in \mathbb{N}\} \tag{4}$$

其中 $\mathbf{1}(\cdot)$ 为指示函数，$\delta$ 代表包含连续函数 $f^{(k)}_n(·)$ 的时间区间。例如，$\delta$ 可以表示日、月或年，并且由用户根据数据和要解决的问题来定义。在第 5 节中，我们定义 $\delta$ 来表示天数，因此，对于每个 $k$，$\{f^{(k)}_n (u)\}_n$ 是一个日函数型时间序列，其中 $n$ 表示天数的索引，$u$ 表示时间在一天之内。然后，利用 式 (4)，我们将 $\{\mathbf{f}_n\}_n$ 作为 $K$ 维函数型时间序列，其中 $\mathbf{f}_n = \{f^{(1)}_n (u), \ldots , f^{(k)}_n (u)\}^T$。这种方法允许我们预测不同的时间范围，例如，小时、天、周、月或年。

由于维数较高，我们建议使用 函数型动态因子模型（ functional dynamic factor model ） 对 $\{\mathbf{f}_n\}_n$ 建模（Hays 等，2012 年^[14]；Gao 等，2019 年^[9]；Martnez-Hernandez 等，2022 年^[28]）。该模型使用 $\{\mathbf{f}_n\}_n$ 的长期协方差的特征函数。具体来说，我们使用 动态函数型主组分 (FDC) 模型 (Hormann等, 2015 ^[16]) 来近似函数型时间序列，然后使用因子模型对分数的多变量时间序列进行建模：

$$\begin{align*} \mathbf{f}_n(u) &= \boldsymbol{\Lambda}(u) \boldsymbol{\beta}_n + \boldsymbol{\eta}_n(u)\\ \boldsymbol{\beta}_{p,n} &= \mathbf{A}_p \boldsymbol{\theta}_{p,n} + \mathbf{e}_{p,n} \\ \boldsymbol{\theta}_{p,n} − \sum^{q1}_{i=1} \mathbf{G}_i \boldsymbol{\theta}_{p,n−i} &= \boldsymbol{\epsilon}_n + \sum^{q2}_{j=1} \mathbf{M}_j \boldsymbol{\epsilon}_{p,n−j} \tag{5} \end{align*}$$

其中 $\boldsymbol{\boldsymbol{\Lambda}(}u)$ 是 $K×p_0K$ 的块矩阵，其对角块为 $1×p_0$ 的 $\boldsymbol{\lambda}^{(k)}(u) = \{\lambda^{(k)}_1 (u),\ldots , \lambda^{(k)}_{p_0} (u)\}^T$，是 $\{f^{(k)}_n \}$ 的长期协方差的特征函数，$k = 1, \ldots , K$, 且 $p_0$ 是主组分的数量，其他地方为零； $\boldsymbol{\beta}_n = (\beta^{(1)}_{1,n}, \ldots , \beta^{(1)}_{p_0,n}, \beta^{(2)}_{1,n},\ldots , \beta^{(K)}_{p_0,n})^T$ 是包含所有 $K$ 个函数型时间序列的函数型主组分分数构成的 $p_0K$ 维向量；$\boldsymbol{\eta}_n(u)$ 表示从 $p_0 + 1$ 到无穷的剩余主组分项。与普通函数型主组分分析不同，动态函数型主组分分析考虑了函数型时间序列 $\mathbf{f}_n$ 的时间依赖性。

对于每个主组分 $p = 1, \ldots , p_0$, $\boldsymbol{\beta}_{p,n} = (\beta^{(1)}_{p,n}, \beta^{(2)}_{p,n}, \ldots , \beta^{(K)}_{p,n})^T$ 是所有 $K$ 个函数型时间序列的第 $p$ 个函数型主组分分数的向量系列; $\mathbf{A}_p$ 是因子载荷的 $K × L$ 常数矩阵； $\boldsymbol{\theta}_{p,n}$ 为因子时间序列的 $L$ 维向量；误差项 $\mathbf{e}_{p,n}$ 表示因子时间序列无法解释的内容；$\mathbf{G}_i$ 和 $\mathbf{M}_j$ 是代表因子时间序列随时间演化的 $L×L$ 系数矩阵；$\boldsymbol{\epsilon}_{p,n}$ 是白噪声向量。动态函数型主组分分析将函数型时间序列简化为由函数型主组分分数组成的向量时间序列。与经典的函数型主组分分析不同，它还考虑了时间依赖性。

式 (3) 和 式 (5) 中的模型在预测中提供了广泛的灵活性，这对多种应用很有吸引力。例如，风速数据预测的时间尺度从几分钟到几天不等（短期、中期和长期预测）。大多数研究都侧重于短期、中期或长期预测，但没有人仅使用一个模型来预测这些不同时间尺度中的每一个。可以预测不同时间尺度的模型在计算上很方便。它提供了解释的灵活性，无需估计多个模型。

使用我们的 式 (5)，短期和中期预测是通过提前一步预测 $\hat{f}_{N+1}$ 获得的。长期预测是通过提前 $h$ 步的预测 $\hat{f}_{N+h}$ 获得的，其中 $h > 1$（详见第 5 节）。

备注 1：如果均值 $\mu_{t}$ 不为零，则 $\{\mu_{t}\}$ 可以与系数 $\mathbf{b}_t$ 一起建模。也就是说，我们可以不使用 式 (4) 中的 $b_{k,u}$，而是使用 $μ_u + b_{k,u}$。或者，可以与 $\mathbf{b}_t$ 分开建模和预测 $\{\mu_{t}\}$，例如，如果 $\mu_{t}(\mathbf{s}) = \mu_{t}$，则使用单变量函数型时间序列，其中单变量函数型时间序列的构造与 式 (4) 中的类似。

3 推断与预测

3.1 推断

为了估计连续随机场 $x_{t}$，我们需要估计 式（3） 中的系数 $b_{k,t}$ 和均值 $\mu_{t}$（可能取决于 $\mathbf{s}$）。对于此估计，我们使用 R-INLA 框架（Rue 等，2009 年）^[36]，它要求我们在超参数 $\tau_{\varepsilon,t}$、$σ_{x,t}$ 和 $\rho_t$ 上设置先验。先验的选择取决于数据的特征。我们按照 Fuglstad 等（2019）^[7]中描述的相同方式在 $σ_{x,t}$ 和 $\rho_t$ 上分配一个联合先验 。也就是说，联合先验是使用惩罚复杂性 (PC) 的概念间接指定的（参见 Simpson 等，2017 年^[39]）。我们只需要指定尾部概率 $P (\rho_t < \rho_0) = p_{\rho_t}$ 和 $P (σ_{x,t} > σ_0) = p_{σ_{x,t}}$。换句话说，我们只需要为极差设置下尾分位数 $\rho_0$ 和概率 $p_{\rho_t}$，为标准差设置上尾分位数 $σ_0$ 和概率 $p_{σ_{x,t}}$。对于 $\tau_{\varepsilon,t}$，我们假设参数为 $1$ 和 $0.00005$ 的模糊 Gamma 先验。

对于 式（5） 的模型，我们需要估计 $\boldsymbol{\Lambda}(u)$、$\mathbf{A}_p$、$\mathbf{G}_i$ 和 $\mathbf{M}_j$。为了估计 $\boldsymbol{\Lambda}(u)$，我们使用动态函数型主组分分析。一旦估计了 $\boldsymbol{\Lambda}(u)$，我们就得到了 $\boldsymbol{\beta}_n$，然后我们估计因子模型的组成部分。因子模型中的一个常见假设是因子相互独立，我们在此也作同样的假设。分开对每个因子时间序列 $\{θ^{(l)}_{p,n}\}_n, l = 1, . . . , L$ 建模。因此， 式（5） 中的矩阵 $\mathbf{G}_i$ 和 $\mathbf{M}_j$ 是对角矩阵。使用 ftsa 包完成此估算（Hyndman 和 Shang，2020 ^[18]）。下面详细介绍 式（5）的估计过程。

令 $c^{(k)}(u, v) = \sum^{\infty}_{h=−\infty} c^{(k)}_h(u, v)$ 为函数型时间序列 $\{f^{(k)}_n \}$ 的长期协方差，其中 $c^{(k)}_h (u, v) = \operatorname{Cov}\{f^{(k)}_n (u), f^{(k)}_{n+h}(v)\}$ 是滞后 $h$ 处的自协方差函数。我们将长期协方差算子定义为 $C^{(k)}(f )(u) = \int_T c^{(k)}(u, v)f(v)dv$。 $\lambda^{(k)}_p$ 的估计量 $\hat{\lambda}^{(k)}_p$ 被定义为经验长期协方差算子的特征函数。也就是说， $\hat{\lambda}^{(k)}_p$ 使得 $\hat{C}^{(k)} ( \hat{\lambda}^{(k)}_p )(u) = \hat{ζ}^{(k)}_p \hat{\lambda}^{(k)}_p (u)$，其中 $\hat{ζ}^{(k)}_p$ 是降序排列的相应特征值。估计算子 $\hat{C}^{(k)}$ 使用 Rice 和 Shang (2017) ^[35] 中描述的带宽选择方法获得。明确地，$\hat{c}^{(k)}(u, v) := \sum_{|h|≤h} W(h/q) \hat{c}^{(k)}_h (u, v)$，其中 $\hat{c}^{(k)}_h (u, v)$ 是滞后 $h$ 处的经验自协方差函数，$W$ 是对称的权重函数。然后， 式（5） 中的 $\boldsymbol{\beta}_n$ 的分量为 $\tilde{β}^{(k)}_{p,n} = \int_{\mathcal{T}} f^{(k)}_n (u) \hat{\lambda}^{(k)}_p (u)du$。

我们使用估计的 $\tilde{\boldsymbol{\beta}}_{p,n}$ 向量来估计 $\mathbf{A}_p$、$\mathbf{G}_i$ 和 $\mathbf{M}_j$。一种常见方法是通过矩阵的特征分解来估计 $\mathbf{A}_p$ 和 $\boldsymbol{\theta}_{p,n}$，该矩阵涉及 $\boldsymbol{\beta}_{p,n}$、$\boldsymbol{\theta}_{p,n}$ 和 $\mathbf{e}_{p,n}$ 在不同滞后 $h$ 时的协方差和交叉协方差（参见，例如，Lam 等， 2011 年^[23]；Gao 等，2019 年 ^[9]）。

令

• $\boldsymbol{\Sigma}^{(p)}_{\boldsymbol{\beta}} (h) = \operatorname{Cov}(\boldsymbol{\beta}_{p,n+h}, \boldsymbol{\beta}_{p,n})$,
• $\boldsymbol{\Sigma}^{(p)}_{\boldsymbol{\theta}}(h) = \operatorname{Cov}(\boldsymbol{\theta}_{p,n+h}, \boldsymbol{\theta}_{p,n})$,
• $\boldsymbol{\Sigma}^{(p)}_{\mathbf{e}}(h) = \operatorname{Cov}(\mathbf{e}_{p,n+h}, \mathbf{e}_{p,n})$

分别代表各自的协方差矩阵。类似地，令 $\boldsymbol{\Sigma}^{(p)}_{\boldsymbol{\theta},\mathbf{e}}(h)$ 为滞后 $h$ 处 $\boldsymbol{\theta}_{p,n}$ 和 $\mathbf{e}_{p,n}$ 之间的互协方差矩阵。

则从 式（5） 中，我们有：

$$\mathbf{L}^{(p)} = \mathbf{L}^{(p)∗} + \mathbf{E}^{(p)}$$

其中 ：

• $\mathbf{L}^{(p)} = \sum^{h_0}_{h=1} \boldsymbol{\Sigma}^{(p)}_{\boldsymbol{\beta}}(h) \boldsymbol{\Sigma}^{(p)}_{\boldsymbol{\beta}}(h)^T$，
• $\mathbf{L}^{(p)∗} = \mathbf{A}_p \{\sum^{h_0}_{h=1} \mathbf{B}(h) \mathbf{B}(h)^T \} \mathbf{A}_p^T$，
• $\mathbf{E}^{(p)} = \mathbf{A}_p \mathbf{B}(h) \boldsymbol{\Sigma}^{(p)}_{\mathbf{e}}(h)^T + \boldsymbol{\Sigma}^{(p)}_{\mathbf{e}}(h) \mathbf{B}(h) \mathbf{A}_p^T + \boldsymbol{\Sigma}^{(p)}_{\mathbf{e}}(h) \boldsymbol{\Sigma}^{(p)}_{\mathbf{e}}(h)^T$
• $\mathbf{B}(h) = \boldsymbol{\Sigma}^{(p)}_{\boldsymbol{\theta}}(h) \mathbf{A}_p^T + \boldsymbol{\Sigma}^{(p)}_{\boldsymbol{\theta},\mathbf{e}}(h) + \boldsymbol{\Sigma}^{(p)}_{\mathbf{e},\boldsymbol{\theta}} (h)$。

然后，$\mathbf{A}_p$ 的自然估计量定义为 $\hat{\mathbf{A}}_p = (\hat{\mathbf{a}}_{p,1},\ldots,\hat{\mathbf{a}}_{p,L})$，其中 $\hat{\mathbf{a}}_{p,l}$ 是 $\hat{\mathbf{L}}^{(p)}$ 的第 $l$ 个特征向量，并且 $\hat{\mathbf{L}}^{(p)} = \sum^{h_0}_{h=1} \hat{\boldsymbol{\Sigma}}^{(p)}_{\boldsymbol{\beta}}(h) \hat{\boldsymbol{\Sigma}}^{(p)}_{\boldsymbol{\beta}}(h)^T$，其中 $\hat{\boldsymbol{\Sigma}}^{(p)}_{\boldsymbol{\beta}}(h)$ 是 $\tilde{\boldsymbol{\beta}}_{p,n}$ 在滞后 $h$ 处的经验协方差矩阵。最后，得到估计的因子时间序列为 $\hat{\boldsymbol{\theta}}_{p,n} = \hat{\mathbf{A}}_p^T \hat{\boldsymbol{\beta}}_{p,n}$。然后，我们使用最大似然法将自回归移动平均 (ARMA) 模型拟合到 $\hat{\boldsymbol{\theta}}_{p,n}$ 的每个分量。我们使用 Akaike 信息准则 (AIC) 来选择 ARMA 模型的阶。

3.2 预报

本节描述如何预测时间 $T + h_1, h_1 > 0$ 时的连续表面时间序列 $x_{t}(\mathbf{s})$。

再次，不失一般性，让我们假设 $\mu_{t}(\mathbf{s}) = 0$。因此，我们需要预测 $\hat{f}^{(k)}_{N+h_2}$，其中 $h_2 = \lceil h_1/\delta \rceil$。从 式（5） 中，我们得到时间依赖性包含在因子过程 $\{\boldsymbol{\theta}_{p,n}\}_n$ 中。由于这些因子互不相关，我们可以分别为 $\{\boldsymbol{\theta}_{p,n}\}_n, p = 1,\ldots,p_0$ 的每个分量建模。这意味着估计 $L_{p_0}$ 阶的 ARMA 模型。利用 ARMA 模型，我们可以获得预测因子 $\hat{\boldsymbol{\theta}}_{p,N+h_2}$；进而能够得到估计 $\hat{\tilde{\boldsymbol{\beta}}}_{p,N+h_2}$。

由此，函数型时间序列在时间 $N$ 的 $h_2$ 步超前预报为：$\hat{f}^{(k)}_{N+h_2}(u) = \hat{\boldsymbol{\lambda}}^{(k)}(u)^T \hat{\tilde{ \boldsymbol{\beta} }}^{(k)}_{N+h_2}$。

预测的 $K$ 个函数型时间序列 $\{\hat{f}^{(k)}_{N+h_2}(u); k = 1,\ldots , K\}$ 包含获得 $\hat{x}_{T +h_1}(\mathbf{s})$ 所需的系数值 $b_{k,T+h_1}$。预报的系数可以定义为：

$$\hat{b}_{k,T + h_1} = \hat{f}^{(k)}_{N+ \lceil h_1/\delta \rceil}(h_1 − \delta \lfloor h_1/\delta \rfloor), k = 1, \ldots , K \tag{6}$$

最后，连续随机场的平均预测为

$$\hat{x}_{T +h_1} (\mathbf{s}) = \sum^K_{k=1} \hat{b}_{k,T+h_1} \psi_k(\mathbf{s}) \tag{7}$$

可以构造一个表面的预测区间 $\hat{x}_{T+h_1}$。这可以从曲线 $\hat{f}^{(k)}_{N+h_2}$ 的预测区间获得。设 $\hat{f}^{(k)}_{N+h_2,α}$ 和 $\hat{f}^{(k)}_{N+h_2,1−α}$ 为 $\hat{f}^{(k)}_{N+h_2}$ 的预测边界，其中 $α$ 是显着性水平。那么：

• $\hat{x}^α_{T+h_1} (\mathbf{s}) = \sum^K_{k=1} \hat{f}^{(k)}_{N+ \lceil h_1/\delta \rceil,α} (h_1 − \delta \lfloor h_1/\delta \rfloor) \psi_k(\mathbf{s})$
• $\hat{x}^{1−α}_{T+h_1} (\mathbf{s}) = \sum^K_{k=1} \hat{f}^{(k)}_{N+ \lceil h_1/\delta \rceil ,1−α} (h_1 − \delta \lfloor h1/\delta \rfloor )\psi_k(\mathbf{s})$

分别对应于预测的下界和上界，其中 $α$ 为显著性水平。

为了获得 $\hat{f}^{(k)}_{N+h_2,α}$ 和 $\hat{f}^{(k)}_{n +h2,1−α}$，我们使用 Gao 等 （2019） ^[9] 中描述的自举方法。

备注 2： 如果 $\mu_{t}(\mathbf{s}) \neq 0$ 且随时间 $t$ 变化，则 $\mu_{t}(\mathbf{s})$ 的预测类似于 $z_{t}(\mathbf{s})$ 的预测。此外，如果 $\mu_{t}(\mathbf{s})$ 在 $D$ 中是常数，则可以使用单变量函数型时间序列预测方法获得 $\hat{\mu}_{T +h_1}$（Martnez-Hernandez 和 Genton，2021 ^[27]）。

上述过程可通过 R-INLA 和 ftsa 包实现。

4 模拟研究

在这里，我们说明了我们方法的优势。我们提出了两项模拟研究，并将我们的方法与实践中经常使用的可分离时空模型进行了比较（例如，Lenzi 和 Genton，2020 ^[24]）。可分离时空模型被定义为空间域的随机偏微分方程模型和时间维度的 AR(1) 模型。该模型可以明确定义如下：

$$\begin{align*} y_{t}(\mathbf{s}) = μ + ξ_t(\mathbf{s}) + \varepsilon_{t}(\mathbf{s}) \\ ξ_t(\mathbf{s}) = aξ_{t−1}(\mathbf{s}) + W_{t}(\mathbf{s}) \end{align*}$$

其中 $\varepsilon_{t}(\mathbf{s})$ 是测量误差，$ξ_t(\mathbf{s})$ 是遵循一阶自回归动态模型的隐过程，$W_{t}(\mathbf{s})$ 是均值为零的时间独立高斯随机场。空间相关的新项 $W_{t}(\mathbf{s})$ 以空间协方差函数为特征，并且该协方差函数由 Matern 函数定义。我们使用 R-INLA 实现 式（8） 的拟合。在实际工作中，式（8） 因其计算优势而被广泛使用；参见 Cameletti 等 (2013)^[3] 了解更多详情。

在每个模拟研究中，我们将本文方法与 式（8） 结果进行预测准确性比较。对于每项研究，我们模拟了 $1000$ 个蒙特卡罗数据集。

4.1 仿真设置

我们用 式（1） 来模拟时空数据 $\{y_{t}(\mathbf{s}); \mathbf{s} \in D, t = 1, . . . , T \}$。为此，我们指定平均值 $\mu_{t}$ 并模拟 式（3） 中的 $z_{t}(\mathbf{s})$，为此需要指定域 $D$ 和 $D$ 的三角剖分。我们将 $D$ 定义为沙特阿拉伯国家地理边界（见图 1）。然后，得到 $\psi_1,\ldots, \psi_{k=78}$，它们是在 $D$ 的三角网格上定义的基函数。我们模拟测量误差，使得 $\epsilon_{t}(\mathbf{s}_i) \stackrel{i.i.d} \sim \mathcal{N} (0, 0.25)$。我们考虑了 $T = 960$ 个时间点。为了拟合模型，我们在 $100$ 个位置 $\mathbf{s}_1, \ldots , \mathbf{s}_100$ 处计算了 $y_{t}(\mathbf{s}), t = 1,\ldots,T$。这 $100$ 个位置均匀分布在 $D$ 上。上述值的选择是为了能够比较本文模型和 式（8）。尽管本文模型可以处理更大的数据集，但我们无法为更大的 $T$ 值估计 式（8） 。

为了评估模型正确预测的能力，将数据分成 $912$ 个时间点的训练集和其余 $48$ 个时间点的测试集。也就是说，通过预测剩下的 $48$ 个时间点来评估预测能力。我们量化了 $100$ 个位置的预测均方误差 (MSE)：

$$\sum^{100}_{i=1} \{x_{t}(\mathbf{s}_i) − \hat{x}_{t}(\mathbf{s}_i)\}^2/100, t = 913, . . . , 960$$

4.2 模拟 1

首先，我们使用 式（5） 模拟 $z_{t}(\mathbf{s})$ 的系数并设置 $\delta = 24$，这产生 $N = T /\delta = 40$ 条曲线。也就是说，我们模拟多元函数型时间序列 $\{\mathbf{f}_n(u); u \in [0, 1], n = 1, \ldots , N = 40\}$，分量 $f^{(k)}_n (u) = \sum^3_{p=1} β^{(k)}_{p,n} \lambda^{(k)}_p(u)$。特征函数定义为 $\lambda^{(k)}_1(u) = \sin(2πu + πk/2)$、$\lambda^{(k)}_2(u) = \cos(2πu + πk/2)$ 和 $\lambda^{(k)}_3 (u) = \sin(4πu + πk/2)$。系数 $\boldsymbol{\beta}_{p,n} = (\beta^{(1)}_{p,n},\ldots, \beta^{(K)}_{p,n})^T$ 被模拟为 $\boldsymbol{\beta}_{p,n} = \mathbf{A}_p \boldsymbol{\theta}_{p,n}$，其中 $\mathbf{A}_p = (a^{(p)}_{ij} )_{1 ≤i,j≤K}$，其中 $a^{(p)}_{ij} = K^{−1/4}b^{(p)}_{ij}$；对于 $p = 1, 2$，取 $b^{(p)}_{ij} \stackrel{i.i.d} \sim \mathcal{N}(2, 4)$； $b^{(3)}_{ij} \stackrel{i.i.d} \sim \mathcal{N}(0, 0.04)$。

因子时间序列 $\boldsymbol{\theta}_{p,n}$ 是从 1 阶和 2 阶的自回归模型模拟的，即 AR(1) 和 AR(2)。因子时间序列 $\boldsymbol{\theta}_{p,n}$ 的每个分量 $θ^{(l)}_{p,n}$ 被模拟为：

• $θ^{(1)}_{p,n} = 0.5 θ^{(1)}_{p,n−1} + 0.2θ^{(1)}_{p,n− 2} + \epsilon_n$。
• 对于 $l = 2，。 . . , K$, $θ^{(l)}_{p,n} = K^{−1}z^{(l)}_n$，其中 $z^{(l)}_n = 0.2 z^{(l)}_{n−1} + \epsilon_n$
• $\epsilon_n$ 在两种情况下均表示具有单位方差的高斯白噪声。

其次，我们将均值分量 $\mu_{t}$ 定义为均值函数 $μ(u) = 3u \cos(uπ) + 3.6u + 12, u \in [0, 1]$ 的估值结果。

最终，获得的随机场为 $x_{t}(\mathbf{s}) = \mu_{t} + \sum^K_{k=1} b_{k,t} \psi_k(\mathbf{s})$，其中 $b_{k,t}$ 是 $f^{(k)}_n (u)$ 的估值。

请注意，我们模拟了与基数量一样多的因子，$L = K$，其中除了 $l = 1$ 外，第 $k$ 个因子的贡献被降低了 $1/K$。这个模拟提供了一个更复杂的结构来模拟真实数据集。不过，在估计 式（5） 时，我们只估计了三个因子，$L = 3$。

4.3 模拟 2

在本场景中，我们从 式（8） 中模拟 $x_{t}(\mathbf{s})$。也就是说，我们定义 $x_{t}(\mathbf{s}) = μ + ξ_t(\mathbf{s})$，其中 $μ = 10$。过程 $ξ_t(\mathbf{s})$ 被模拟为 $ξ_t(\mathbf{s}) = 0.8ξ_{t−1}(\mathbf{s}) + \sqrt{1 − 0.8^2} W_{t}(\mathbf{s})$，其中 $ξ_1(\mathbf{s}) = W_1(\mathbf{s})$。我们使用方差为 $1$ 且缩放参数为 $0.1$ 的 Whittle 协方差函数生成空间模型 $W_{t}(\mathbf{s})$ 的 $T$ 个独立实现。这大约对应于 $28.2$ 的变程。我们添加项 $\sqrt{1 − 0.8^2}$ 以使过程在时间上平稳。虽然我们的方法不需要时间平稳性，但 式（8） 确实需要。

图 2: 模拟 1 和模拟 2 的 MSE 预测值的箱线图。MSE 值对应于模型 1（我们的建议）和模型 2 预测的 48 个时间点。索引 N+1 代表前 24 个时间点和 N+2代表剩下的。在模拟 1 中，模型 1 呈现最低的 MSE 值；在模拟 2 中，两个模型的性能相似。

4.4 仿真结果

在这里，我们展示了上述两种情况下两种模型的模拟结果。图 2 显示了预测的 $48$ 个时间点预测值的 MSE 。在每个子图（场景）中，我们并排展示了两个模型的性能。总体而言，与 式（8）（标记为模型 2）相比，本文方法（标记为模型 1）具有良好的性能。

在模拟 1（图 2 的第一行）中，我们观测到模型 1（我们的建议）在所有时期（$N + 1$ 和 $N + 2$）都呈现出较低的预测 MSE。这些时期对应于预测的 $48$ 个时间点。模型 1 在 $N + 1$ 和 $N + 2$ 期间的均值分别为 $0.915$ 和 $1.012$。对于模型 2，$N+1$ 和 $N+2$ 期间的 MSE 均值分别为 $1.24$ 和 $1.228$。因此，我们得出模型 1 优于模型 2 的结论。这个结论并不奇怪，因为数据生成过程来自模型 1。

在模拟 2 的场景中（图 2 的第二行），两个模型在所有时间段（$N + 1$ 和 $N + 2$）中的表现相似，即穿越所有 $48$ 个时间点。模型 1 在 $N + 1$ 和 $N + 2$ 期间的均值分别为 $0.26$ 和 $0.259$。对于模型 2，我们获得与模型 1 相同的均值。在此情况下，数据生成过程是模型 2。因此，本文方法与模型 2 得到同样执行结果可以被视为模型的一种优势。

总之，基于这两个模拟，我们的方法优于实践中最常用的时空模型。此外，我们的方法为非平稳、复杂和大型数据集建模提供了广泛的灵活性。

5 风数据分析

5.1 数据说明

我们使用 Yip (2018)^[48] 从天气研究和预报 (WRF) 模型中模拟的风数据。每个测量值对应于 2010 年沙特阿拉伯的每小时风速 (m/s)。最初，风速是在空间中规则的网格点上模拟的，即 $5$ 公里分辨率。

为了达到目标，我们随机选择了沙特阿拉伯的 $1500$ 个地点，并分析了四个月的数据：5 月、6 月、7 月和 8 月。因此，数据集包含 $1500$ 个位置和 $2952$ 个时间点，因此总共有 $4,428,000$ 个值。 $1500$ 个位置显示在 图 1 的左侧子图中。令 $ws_t(\mathbf{s}_i)$ 表示位置 $\mathbf{s}_i$ 和时间 $t$ 处的风速。然后，对于每个 $\mathbf{s}_i，i = 1, \ldots , 1500$, $\{ws_t(\mathbf{s}_i)\}_t$ 是长度为 $T = 2952$ 的小时时间序列。为了使用 式（2） 和 式（3） ，我们考虑了风速的平方根变换。转换后的风速类似于高斯分布（Haslett 和 Raftery，1989）^[13]。我们用 $y_{t}(\mathbf{s}_i) = \sqrt{ws_t(\mathbf{s}_i)}$ 表示平方根变换后的风速。

了解风速行为对于可再生能源非常重要。风能是最有效的可再生能源之一，但这种资源的可用性取决于位置。为确保电力系统的可靠性和质量，开发高精度的风速预测方法非常重要。风速预报根据要预报的时间尺度范围进行分类。这些可以分为短期（提前 30 分钟到 6 小时）、中期（提前 6 小时到 1 天）和长期（提前 1 天到 1 周）预测。文献中已经提出了几种风速预测模型的方法，例如具有时间序列技术、空间统计技术或两者的模型都专注于预测特定的时间尺度范围。在文献中可以找到一些提出用于预测短期和中期时间尺度模型的尝试（例如，参见 Lenzi 和 Genton，2020 年 ^[24] 及其中的参考文献）。相比之下，本文模型无需任何额外修改即可预测短期、中期和长期时间尺度。我们的目标是准确预测短期、中期和长期时间尺度范围内的风速数据。

此处考虑的风速数据集在空间和时间上是庞大而复杂的。因此，数据符合在 第 2 节 中描述的建议方法。本文呢方法（或统计模型）相对于 WRF 模型的优势主要在于计算成本。我们定义 式（4） 中的时间段 $\delta$，使 $n$ 代表天数，即 $\delta = 24$（小时）。 $\delta$ 的这种选择适应了昼夜循环，这是一个外部因素，在模型中是可取的（Weisser 和 Foxon，2003）^[42]。

为了评估预测性能，我们使用了四个月的数据来拟合模型，然后预测 9 月的前三天，即提前 72 小时预测。第一步是定义区域的三角剖分，此处为沙特阿拉伯国家行政区域范围。数据分析中使用的三角剖分显示在 图 3 的第一个子图中。一旦我们定义了域的三角剖分，就可以估计连续随机场。我们将在下面详细说明此估计。

图 3. 具有 $364$ 个节点的沙特阿拉伯网格：$\psi_1,\ldots,\psi_{364}$（左上）和八个连续的连续随机场，对应于 2010 年 6 月 22 日，使用 式（3） 的基础表示法进行估计。每个连续随机场的形式为 $x_{t}(\mathbf{s})=\mu_{t}+\sum_{k=1}^{364} b_{k,t} \psi_k(\mathbf{s})$，其中 $t=1251,\ldots,1258$ 代表小时。

5.2 模型规格

为了估计连续风速随机场，我们使用 $364$ 个基函数，$\psi_1, \ldots , \psi_{364}$, 由三角测量确定。我们将范围先验设置为 $P (\rho_t < 500) = 0.5$，将标准差设置为 $P (σ_{x,t} > 10) = 0.01$。先验范围基于沙特阿拉伯的景观特征。图 3 为 8 个相连的连续随机场的结果：$x_{t}(\mathbf{s}) = \mu_{t} + \sum^{364}_{k=1} b_{k,t} \psi_k(\mathbf{s})$，其中 $t$ 代表小时，对应 2010 年 6 月 22 日 的 $03:00, 04:00，\ldots，10:00 AM$。

在估计连续随机场后，我们将系数 $\{b_{k,t}\}$ 转换为函数型时间序列 $f^k_n$ 。如上所述，我们将时间段定义为 $\delta = 24$（小时）以考虑昼夜循环。因此，函数型时间序列的样本大小为 $N = 123$。图 4 显示了对应于基 $k = 20$ 的一个函数型时间序列，即函数型时间序列 $\{f^{(20)}_n(u); u \in T , n = 1,\ldots, 123 \}$。在左上角的子图中，可以观察到这个函数型时间序列在许多天中显示出一种重复的模式。这种每日的时间依赖模式在所有函数型时间序列 $\{f^{(k)}_n \}_n, k = 1,\ldots, K = 364$ 中都能被观察到。然后，我们继续用函数型时间序列模型描述时间依赖性，并预测接下来的三天：2010 年 9 月 1 日、2 日和 3 日。

我们使用函数型时间序列 $f_n = \{f^{(1)}_n (u), \ldots , f^{(364)}_n (u)\}^T$ 来模拟时间依赖性。在式（5）的动态函数型主组分表示中， 我们设置 $p_0 = 3$ 和因子数 $L = 3$。$p_0$ 的选择基于所解释函数型时间序列的可变性比例。通常，在实际数据应用中，$p_0$ 的值通常小于 $5$。如上所述，ARMA 模型使用似然估计，并通过 AIC 选择阶。使用预报值 $\hat{\tilde{\boldsymbol{\beta}}}_{p,N+1}$, $\hat{\tilde{\boldsymbol{\beta}}}_{p,N+2}$ 和 $\hat{\tilde{\boldsymbol{\beta}}}_{p,N+3}$，我们得到对应的 2010 年 9 月 1 日、2 日和 3 日的预报曲线 $\hat{f}_{N+1}$， $\hat{f}_{N+2}$，$\hat{f}_{N+3}$。

图 4. 函数型时间序列及其预测。左上：函数型时间序列 $\{f_n^{(20)}\}_n$ 对应于基函数 $\psi_{20}$，且 $n=1,\ldots,N=123$（天）。右上角和第二行：对应于 2010 年 9 月 1 日、2 日和 3 日的预报曲线，基分量 $k=20$。这些预报曲线是从式（5）的模型中获得的。虚线曲线表示预报下限和上限。

5.3 结果

图 4 显示了函数型时间序列 $f^{(20)}_n$ 的预报曲线。对于每条预报曲线，我们还绘制了 “真实” 的系数曲线。这些真实曲线是根据 2010 年 9 月 1 日、2 日和 3 日的估计曲面 $x_{T+h}(\mathbf{s})$ 的系数获得的。预报曲线有望接近真实曲线。因此，我们可以评估模型的性能。从 图 4 中，我们观察到预报曲线模仿了真实曲线的主要行为。这意味着预测表面在表示空间依赖性和时间依赖性方面应该是准确的。

我们现在可以在域 $T$ 的任何值处评估预报曲线，以获得预报的连续风速场。我们展示了三个不同时间的预测连续风速场：9 月 1 日 05:00、9 月 1 日 16:00 和 9 月 3 日 19:00。

这意味着我们在 $u = 5$ 和 $u = 16$ 时估值 $\hat{\mathbf{f}}_{N+1}$，在 $u = 19$ 时估值 $\hat{\mathbf{f}}_{N+3}$。这三个时间范围分别对应于短期、中期和长期预测。图 5 显示了对在观测数据的 $1500$ 个位置预报的三个表面的估值。此外，图 5 显示了相应的真实风速数据（未用于拟合模型）。我们观测到预报值呈现出与原始数据相同的变化。例如，9 月 1 日 16:00（第二列）的观测风速值在该国西北部（麦地那地区）较高。这些高值沿着红海沿岸延伸；而在该国中部，风速值较低。我们在预报平均值上获得了相同的变化。我们看到预报的最高值位于西北地区和红海沿岸。该结论适用于所有 $72$ 个预报小时。

图 5： 三个不同时间的预报值结果。第一行显示了在 $1500$ 个位置观测到的风速值的平方根。第二行显示了与观测值保持相同空间模式的预报值。第三行和第四行分别显示预报值的 $5\%$ 和 $95\%$ 分位数。

为了更好地了解方法的性能，我们在 2010 年 9 月 1 日、2 日和 3 日的每个小时计算平方预测误差。也就是说，我们为每个 $t$ 计算 $\{y_{t}(\mathbf{s}_i) − \hat{x}_{t}(\mathbf{s}_i)\}_i^2$，其中 $i$ 是观测数据的位置索引，$\hat{x}_{t}$ 是预测表面。图 6 显示了提前一天、两天和三天的每小时预测误差平方值的箱线图，分别对应于 2010 年 9 月 1 日、2 日和 3 日。

我们观测到，对于 9 月 1 日和 3 日，从这些箱线图中可以知道最高误差值是在中午左右，而对于 9 月 2 日，最高误差值是在上午 6:00 到 9:00 之间。这可能是由于太阳的热量产生对流时的高风速导致观测值和预报值差异更大。尽管存在这些错误值，但数据的一般行为已恢复，如图 5 所示。因此，我们得出结论，本文方法具有良好的性能。

图 6： 预报值的误差值。在 2010 年 9 月 1 日、2 日和 3 日的每个小时，我们都会展示值 $\{y_{t}(\mathbf{s}_i)− \hat{x}_t(\mathbf{s}_i)\}_i^2$ 的箱线图，其中 $\hat{x}_t$ 是预报表面。

6 讨论

许多时空模型的一个重要限制是无法处理大型数据集。由于这种能力缺陷，需要一种新的数据分析范式。沿着这条线，我们提出了一种新方法来模拟大型和复杂的时空数据。这种方法假设在每个时间点，观测是一个连续的随机场。连续随机场是使用有限基函数表示来估计的，其中基函数的系数是假设随机场具有 Matern 协方差结构来估计的。然后，使用多元函数型时间序列方法对表面序列进行建模。

我们的方法可以应用于一大类数据，并不局限于风速数据的具体情况。此外，它可以对空间和时间上不稳定的数据进行建模。这是因为连续表面估计不一定需要空间平稳性，并且基函数的系数是用函数型时间序列方法建模的（平稳性条件跨越索引 $n$ 而不是跨越 $t$）。此外，我们的方法可以处理位置和位置数量随时间的变化。

我们方法的一个局限性是不确定性量化。例如，使用我们的方法们无法获得误差预测方差的封闭形式表达式 $x_{T+h_1} (\mathbf{s}) − \hat{x}_{T +h_1}(\mathbf{s})$；这与克里金法形成对比，克里金法具有封闭形式的表达式。一种可能且合乎逻辑的解决方案是使用自举方法来获得不确定性量化。我们方法的另一个局限性是对要使用的因素 $L$ 的数量的估计。

虽然在实际工作中很常见的是设置 $L = 2$ 或 $3$ 通常效果很好，但需要有一种统计方法来选择 $L$。

我们的结论是，我们的方法提供了一种有价值的替代方法来模拟大型时空数据。

参考文献

• [1] Aue, A., D. D. Norinho, and S. Hormann (2015). On the prediction of stationary functional time series. Journal of the American Statistical Association 110 (509), 378–392.
• [2] Bernardi, M. S., L. M. Sangalli, G. Mazza, and J. O. Ramsay (2017). A penalized regression model for spatial functional data with application to the analysis of the production of waste in Venice province. Stochastic Environmental Research and Risk Assessment 31 (1), 23–38.
• [3] Cameletti, M., F. Lindgren, D. Simpson, and H. Rue (2013). Spatio-temporal modeling of particulate matter concentration through the SPDE approach. AStA. Advances in Statistical Analysis. 97 (2), 109–131.
• [4] Chen, W., M. G. Genton, and Y. Sun (2021). Space-time covariance structures and models. Annual Review of Statistics and Its Application 8 (1), 191–215.
• [5] Cressie, N. and G. Johannesson (2008). Fixed rank kriging for very large spatial data sets. Journal of the Royal Statistical Society. Series B 70 (1), 209–226.
• [6] Cressie, N. and C. K. Wikle (2011). Statistics for Spatio-Temporal Data. Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ.
• [7] Fuglstad, G.-A., D. Simpson, F. Lindgren, and H. Rue (2019). Constructing priors that penalize the complexity of Gaussian random fields. Journal of the American Statistical Association 114 (525), 445–452.
• [8] Furrer, R., M. G. Genton, and D. Nychka (2006). Covariance tapering for interpolation of large spatial datasets. Journal of Computational and Graphical Statistics 15 (3), 502–523.
• [9] Gao, Y., H. L. Shang, and Y. Yang (2019). High-dimensional functional time series forecasting: An application to age-specific mortality rates. Journal of Multivariate Analysis 170, 232–243.
• [10] Gelfand, A. E., S. Banerjee, and D. Gamerman (2005). Spatial process modelling for univariate and multivariate dynamic spatial data. Environmetrics 16 (5), 465–479.
• [11] Genton, M. G. (2007). Separable approximations of space-time covariance matrices. Environmetrics 18 (7), 681–695.
• [12] Giraldo, R., P. Delicado, and J. Mateu (2011). Ordinary kriging for function-valued spatial data. Environmental and Ecological Statistics 18 (3), 411–426.
• [13] Haslett, J. and A. E. Raftery (1989). Space-time modelling with long-memory dependence: Assessing ireland’s wind power resource. Journal of the Royal Statistical Society. Series C (Applied Statistics) 38 (1), 1–50.
• [14] Hays, S., H. Shen, and J. Z. Huang (2012). Functional dynamic factor models with application to yield curve forecasting. The Annals of Applied Statistics 6 (3), 870–894.
• [15] Heaton, M. J., A. Datta, A. O. Finley, R. Furrer, J. Guinness, R. Guhaniyogi, F. Gerber, R. B. Gramacy, D. Hammerling, M. Katzfuss, F. Lindgren, D. W. Nychka, F. Sun, and A. ZammitMangion (2019). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological and Environmental Statistics 24 (3), 398–425.
• [16] Hormann, S., L. u. Kidzi ́ nski, and M. Hallin (2015). Dynamic functional principal components. Journal of the Royal Statistical Society. Series B 77 (2), 319–348.
• [17] Huang, H., S. Abdulah, Y. Sun, H. Ltaief, D. E. Keyes, and M. G. Genton (2021). Competition on spatial statistics for large datasets. Journal of Agricultural, Biological and Environmental Statistics 26 (4), 580–595.
• [18] Hyndman, R. J. and H. L. Shang (2020). ftsa: Functional Time Series Analysis. R package version 5.8.
• [19] Hyndman, R. J. and M. S. Ullah (2007). Robust forecasting of mortality and fertility rates: a functional data approach. Computational Statistics & Data Analysis 51 (10), 4942–4956.
• [20] Kaufman, C. G., M. J. Schervish, and D. W. Nychka (2008). Covariance tapering for likelihoodbased estimation in large spatial data sets. Journal of the American Statistical Association 103 (484), 1545–1555.
• [21] Kokoszka, P. and M. Reimherr (2017). Introduction to functional data analysis. Texts in Statistical Science Series. CRC Press, Boca Raton, FL.
• [22] Kuusela, M. and M. L. Stein (2018). Locally stationary spatio-temporal interpolation of argo profiling float data. Proceedings of the Royal Society A 474 (2220), 20180400.
• [23] Lam, C., Q. Yao, and N. Bathia (2011). Estimation of latent factors for high-dimensional time series. Biometrika 98 (4), 901–918.
• [24] Lenzi, A. and M. G. Genton (2020). Spatiotemporal probabilistic wind vector forecasting over Saudi Arabia. The Annals of Applied Statistics 14 (3), 1359–1378.
• [25] Lindgren, F., H. Rue, and J. Lindstr ̈om (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society. Series B 73 (4), 423–498.
• [26] Martnez-Hernandez, I. and M. G. Genton (2020). Recent developments in complex and spatially correlated functional data. Brazilian Journal of Probability and Statistics 34 (2), 204–229.
• [27] Martnez-Hernandez, I. and M. G. Genton (2021). Nonparametric trend estimation in functional time series with application to annual mortality rates. Biometrics 77 (3), 866–878.
• [28] Martnez-Hernandez, I., J. Gonzalo, and G. Gonz ́ alez-Far ́ıas (2022). Nonparametric estimation of functional dynamic factor model. Journal of Nonparametric Statistics. To appear.
• [29] Matern, B. (1986). Spatial variation (Second ed.), Volume 36 of Lecture Notes in Statistics. Springer-Verlag, Berlin. With a Swedish summary.
• [30] Menafoglio, A., O. Grujic, and J. Caers (2016). Universal kriging of functional data: Tracevariography vs cross-variography? Application to gas forecasting in unconventional shales. Spatial Statistics 15, 39–55.
• [31] Nguyen, H., M. Katzfuss, N. Cressie, and A. Braverman (2014). Spatio-temporal data fusion for very large remote sensing datasets. Technometrics 56 (2), 174–185.
• [32] Nychka, D., D. Hammerling, M. Krock, and A. Wiens (2018). Modeling and emulation of nonstationary Gaussian fields. Spatial Statistics 28, 21–38.
• [33] R Core Team (2022). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.
• [34] Ramsay, J. O. and B. W. Silverman (2005). Functional Data Analysis (Second ed.). Springer Series in Statistics. Springer, New York.
• [35] Rice, G. and H. L. Shang (2017). A plug-in bandwidth selection procedure for long-run covariance estimation with stationary functional time series. Journal of Time Series Analysis 38 (4), 591–609.
• [36] Rue, H., S. Martino, and N. Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society. Series B. 71 (2), 319–392.
• [37] Sangalli, L. M., J. O. Ramsay, and T. O. Ramsay (2013). Spatial spline regression models. Journal of the Royal Statistical Society. Series B 75 (4), 681–703.
• [38] Sigrist, F., H. R. K ̈ unsch, and W. A. Stahel (2012). A dynamic nonstationary spatio-temporal model for short term prediction of precipitation. The Annals of Applied Statistics 6 (4), 14521477.
• [39] Simpson, D., H. Rue, A. Riebler, T. G. Martins, and S. H. Sørbye (2017). Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science 32 (1), 1–28.
• [40] Stein, M. L. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer Series in Statistics. Springer-Verlag, New York.
• [41] Sun, Y., B. Li, and M. G. Genton (2012). Geostatistics for large datasets. In E. Porcu, J.M. Montero, and M. Schlather (Eds.), Advances and Challenges in Space-time Modelling of Natural Events, Berlin, Heidelberg, pp. 55–77. Springer Berlin Heidelberg.
• [42] Weisser, D. and T. Foxon (2003). Implications of seasonal and diurnal variations of wind velocity for power output estimation of a turbine: a case study of grenada. International Journal of Energy Research 27 (13), 1165–1179.
• [43] Wikle, C. K. and N. Cressie (1999). A dimension-reduced approach to space-time Kalman filtering. Biometrika 86 (4), 815–829.
• [44] Wikle, C. K. and M. B. Hooten (2010). A general science-based framework for dynamical spatiotemporal models. TEST 19 (3), 417–451.
• [45] Wood, S. (2017). Generalized Additive Models: An Introduction with R (2 ed.). Chapman and Hall/CRC.
• [46] Wood, S. N. (2006). Low-rank scale-invariant tensor product smooths for generalized additive mixed models. Biometrics 62 (4), 1025–1036.
• [47] Xiao, L., Y. Li, and D. Ruppert (2013). Fast bivariate P-splines: the sandwich smoother. Journal of the Royal Statistical Society. Series B 75 (3), 577–599.
• [48] Yip, C. M. A. (2018). Statistical characteristics and mapping of near-surface and elevated wind resources in the Middle East. Ph.D. Thesis, King Abdullah University of Science and Technology.