0 概述

时空数据是指同时具有时间和空间维度的数据^[1]，传感器、移动电话、射频识别（RFID）和智能电网等智能设备的发展促进了实时时空数据流的采集。考虑一个时空随机过程，时空建模的目标是基于时空数据构建时空模型以对给定时刻所有位置的行为进行预测^[2]。若不考虑时间因素，可以采用单纯的空间模型进行建模，如kriging方法^[3]，但是其准确性较低，添加时间维度可以提高预测的准确性。随着时空数据流采集难度的降低，时空序列建模逐渐成为学者们的研究热点之一。

目前，关于时空序列预测主要分为基于物理模型、基于统计模型和基于机器学习的3种方法。基于物理模型的方法首先对时空序列的机理进行研究，寻找其内在的系统动力学规律并构建系统动力学模型，然后对时空序列进行表达进而预测。其中，JONES等人^[4]提出随机偏微分方程以描述连续型时空随机过程。基于统计模型的方法主要分为描述型时空模型和动态时空模型两类。前者利用统计学中的描述型统计量表达时空模型的性质，并对统计量进行建模以消除随机误差；后者考虑时间和空间的自相关性并利用过去和其他地区的数据对当前数据进行建模，然后实现迭代更新预测。近年来，随着机器学习和深度学习的不断发展，人工智能技术被广泛应用于时空序列建模和预测任务。通过机器学习模型能够提取时空序列中复杂的特征模式，也可以对高维时空序列进行降维和聚类从而使得分析更简便。

本文提出一种分离时空数据中的时间模式和空间模式并分别建模的方法。对原始数据进行平稳性检验并中心化，利用奇异值分解（SVD）分解中心化的数据集，通过时间序列模型中经典的ARIMA模型对时间模式建模并检验其有效性，然后利用ARIMA模型预测时间序列，将预测结果与空间模式相结合并对真实时空序列进行重建，以得到各个地理观测点的预测值。

1 相关工作

描述型时空模型较早以时空协方差为研究对象，通过对样本协方差值进行曲面拟合获得协方差函数，然后利用协方差函数分析时空模式的演变。时空kriging方法^[5-7]基于时空过程的协方差函数给出未知地区给定时刻的最优线性无偏估计。但是，由于时空协方差是一种描述型统计量，很难解释时空模式的内在动态变化。描述型时空模型在数学上更通用，但是动态时空模型在科学上有更强的解释性^[8]。动态时空模型基于条件概率分布进行建模，其中最主要的动态时空模型为层次时空模型。层次时空模型可分为2个主要类别：一类是经验层次模型，其认为观测到的时空过程是真实时空过程的演变以及真实过程通过某种函数作用产生观测过程，其机制类似于隐马尔可夫过程；另一类是贝叶斯层次模型，和经验层次模型的主要区别在于，贝叶斯层次模型认为真实过程中的参数也是动态变化的，其在经验层次模型的基础上增加了底层的参数过程，因此，贝叶斯层次时空模型将时空序列过程分解为参数过程、真实过程和数据过程3个层次并分别建模^[8-10]。

无论经验层次模型还是贝叶斯层次模型，真实过程都是最重要的，其对理解时空动态变化模式具有重要意义。因此，时空模型的一个研究重点在于真实过程的模型构建。统计时空模型的构建主要来源于时间模型和空间模型的结合。CLIFF和ORD较早将时间序列模型应用于空间分析中，提出空间自回归模型（SAR）、空间移动平均模型（SMA）和空间回归模型（SR）等^[11]。MARTIN和OEPPEN将空间信息整合到传统的ARIMA模型^[12]中，提出STARMA模型^[12]。STARMA定义了空间阶次的概念并在真实应用中产生了良好效果^[13-15]。但是，随着时空数据的概念外延，STARMA模型中关于欧式距离越小则空间阶次越低的假设越来越难以满足，使得其在一些未知空间相关性结构的数据集中表现较差。BESSA等人结合其他地区的历史数据和待预测地区的数据，构建向量自回归模型VAR^[16]以对时空序列进行建模描述。但是，VAR模型中的待估计参数空间较大，一方面需要消耗极大的计算资源，另一方面可能由于样本量不足而引起过拟合问题。因此，基于Lasso的VAR（Lasso-VAR）模型被广泛应用^[17]，尽管Lasso-VAR在一定程度上解决了模型过拟合问题，但是其优化模型变得更难求解，计算成本过高。

BAHADORI等人^[18]通过将时空数据作为张量进行处理，提出一种低秩张量学习框架以进行多元时空序列分析。BAROCIO等人^[19]通过动态模式分解的方式对时空数据进行降维并提取时空特征。LI^[20]利用梯度提升回归树（Gradient Boosting Regression Tree，GBRT）算法对城市共享单车的时空数据进行建模并预测数量。在深度学习方法中，递归神经网络（Recurrent Neutral Network，RNN）和深度神经网络（Deep Neutral Network，DNN）被广泛应用于时空序列模型构建任务。SHI^[21]利用RNN模型的一个变体，即长短时记忆（Long and Short Term Memory，LSTM）网络对地区的降雨量进行预测。CHE等人^[22]将传统的RNN拓展到时空领域，提出时空递归神经网络（Spatio-Temporal Recurrent Neural Network，ST-RNN），以对时空序列进行建模预测。类似地，在深度学习方面，ZHANG等人^[23]将深度残差网络拓展到时空领域，提出时空深度残差模型（ST-ResNet）以对人流量进行预测。

2 算法描述 2.1 时空数据的奇异值分解

SVD是一种矩阵分解技术，其在信号处理和统计学中有很多应用^[24]。给定一个秩为$ l $的时空数据矩阵$ {\boldsymbol{Y}}_{D\times T} $，其中，$ D $表示空间中观测点的个数，$ T $表示采样时间点的个数。时空数据矩阵$ {\boldsymbol{Y}}_{D\times T} $的奇异值分解如下：

$ \boldsymbol{Y}=\boldsymbol{U}\boldsymbol{S}{\boldsymbol{V}}^{\text{'}}={s}_{1}{\boldsymbol{u}}_{1}{\boldsymbol{v}}_{1}^{\text{'}}+{s}_{2}{\boldsymbol{u}}_{2}{\boldsymbol{v}}_{2}^{\text{'}}+\cdots +{s}_{l}{\boldsymbol{u}}_{l}{\boldsymbol{v}}_{l}^{\text{'}} $

(1)

其中，$ \boldsymbol{U}=\left({\boldsymbol{u}}_{1}, {\boldsymbol{u}}_{2}, \cdots , {\boldsymbol{u}}_{l}\right) $，$ \boldsymbol{V}=\left({\boldsymbol{v}}_{1}, {\boldsymbol{v}}_{2}, \cdots , {\boldsymbol{v}}_{l}\right) $，$ \boldsymbol{S}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{s}_{1}, \right. $ $ \left.{s}_{2}, \cdots , {s}_{l}\right\} $，且$ {s}_{1}\ge {s}_{2}\ge \cdots \ge {s}_{l}\succ 0 $。向量$ {\boldsymbol{u}}_{i} $（$ i=1, 2, \cdots , l $）是左奇异矩阵的列向量，向量$ {\boldsymbol{v}}_{i}^{\text{'}} $（$ i=1, 2, \cdots , l $）是右奇异矩阵的行向量，标量$ {s}_{i} $称为奇异值。

假设$ \left\{{\boldsymbol{c}}_{m}:m=\mathrm{1, 2}, \cdots , T\right\} $是矩阵$ {\boldsymbol{Y}}_{D\times T} $的列向量，$ {\boldsymbol{c}}_{m} $代表给定的$ m $时刻$ D $中所有空间单元的观测值，$ \boldsymbol{Y}{\boldsymbol{Y}}^{\text{'}} $表达了$ D $空间单元之间的相关性，这里假定$ {\boldsymbol{Y}}_{D\times T} $已经去中心化为零均值矩阵。矩阵$ {\boldsymbol{Y}}_{D\times T} $的$ {\boldsymbol{u}}_{i} $事实上是相关矩阵$ \boldsymbol{Y}{\boldsymbol{Y}}^{\text{'}} $的特征向量，$ {\boldsymbol{u}}_{1} $表示相关矩阵$ \boldsymbol{Y}{\boldsymbol{Y}}^{\text{'}} $对应特征值最大的特征向量，包含了空间相关性最多的信息量，或被称为“空间模态”，表征了空间相关性的模式。$ {\boldsymbol{u}}_{i} $的第$ j $个分量$ {u}_{ij} $表示第$ j $地区对第$ i $空间模态的“贡献”。类似地，假定$ \left\{{\boldsymbol{r}}_{n}^{\mathrm{\text{'}}}:n=\mathrm{1, 2}, \cdots , D\mathrm{ }\right\} $是矩阵$ {\boldsymbol{Y}}_{D\times T} $的行向量，$ {\boldsymbol{r}}_{n}^{\mathrm{\text{'}}} $代表$ n $位置在整个时间段的观测值向量，$ {\boldsymbol{Y}}^{\text{'}}\boldsymbol{Y} $表达了不同时刻之间的相关性，矩阵$ {\boldsymbol{Y}}_{D\times T} $的$ {\boldsymbol{v}}_{i}^{\text{'}} $事实上是相关矩阵$ {\boldsymbol{Y}}^{\text{'}}\boldsymbol{Y} $的特征向量，$ {\boldsymbol{r}}_{1} $表示相关矩阵$ {\boldsymbol{Y}}^{\text{'}}\boldsymbol{Y} $对应特征值最大的特征向量，包含了时间相关性最多的信息量，或被称为“时间模态”，表征了时间相关性的模式。$ {\boldsymbol{r}}_{i}^{\mathrm{\text{'}}} $的第$ j $个分量$ {r}_{ij}^{\mathrm{\text{'}}} $表示第$ j $时刻对第$ i $时间模态的“贡献”。时空矩阵$ {\boldsymbol{Y}}_{D\times T} $分解后的S是奇异值矩阵，$ {s}_{i} $表示模式$ i $的重要程度，例如，若$ {s}_{1} $是最大的奇异值，则$ {s}_{1} $对应的模式1具有表征空间模式的最重要的特征。

2.2 基于SVD的时空序列模型

给定历史时空数据矩阵$ {\boldsymbol{Y}}_{D\times (t-1)} $，对当前时刻t的各个地理观测点进行预测的具体步骤如下：

步骤1 通过SVD对中心化后的时空矩阵进行分解。

假定历史时空数据矩阵$ {\boldsymbol{Y}}_{D\times (t-1)} $的秩为$ l $，可利用式（1）得到如下分解：

$ {\boldsymbol{Y}}_{D\times \left(t-1\right)}={s}_{1}{\boldsymbol{u}}_{1}{\boldsymbol{v}}_{1}^{\mathrm{\text{'}}}+{s}_{2}{\boldsymbol{u}}_{2}{\boldsymbol{v}}_{2}^{\mathrm{\text{'}}}+\cdots +{s}_{l}{\boldsymbol{u}}_{l}{\boldsymbol{v}}_{l}^{\mathrm{\text{'}}} $

(2)

SVD有一个重要的性质，定义奇异值占比$ {E}_{r}=\frac{\sum\limits_{i=1}^{r}{s}_{i}}{\sum\limits_{i=1}^{l}{s}_{i}} $，则当$ r <<l $时，$ {E}_{r} $可达到85%以上的水平，剩余的可认为是噪声。因此，通过选取前几个奇异值与对应的左奇异向量和右奇异向量进行重建，可以对矩阵实现降噪，如下：

$ {\tilde{\boldsymbol{Y}}}_{D\times \left(t-1\right)}={s}_{1}{\boldsymbol{u}}_{1}{\boldsymbol{v}}_{1}^{\mathrm{\text{'}}}+{s}_{2}{\boldsymbol{u}}_{2}{\boldsymbol{v}}_{2}^{\mathrm{\text{'}}}+\cdots +{s}_{r}{\boldsymbol{u}}_{r}{\boldsymbol{v}}_{r}^{\mathrm{\text{'}}} $

(3)

步骤2 通过ARIMA^[25]对时间模式进行建模预测。

由于分解之后得到的右奇异向量$ {\boldsymbol{v}}_{i}^{\text{'}} $可以看作时间序列，因此本文利用时间序列中应用最广泛、效果最好的ARIMA模型进行建模。ARIMA的标准模型如下：

$ \begin{array}{l}{\nabla }^{d}{v}_{i, t}=\mu +{\varphi }_{1}{\nabla }^{d}{v}_{i, t-1}+\cdots +{\varphi }_{p}{\nabla }^{d}{v}_{i, t-p}+{\varepsilon }_{t}-\\ {\theta }_{1}{\varepsilon }_{t-1}-\cdots -{\theta }_{q}{\varepsilon }_{t-q}\end{array} $

(4)

其中，$ {\nabla }^{d}{v}_{i, t} $代表t时刻第i个向量的d阶差分，$ {\varepsilon }_{t} $是t时刻均值为0的随机误差，$ \mu $、$ {\varphi }_{i} $（$ i=1, 2, \cdots , p $）、$ {\theta }_{i} $（$ i=1, 2, \cdots , q $）为待估计参数。利用AIC、BIC信息准则和最大似然法进行模型选择和估计，当得到估计好的模型后，利用该模型进行h步向前预测，如下：

$ \begin{array}{l}{\widehat{v}}_{i, t+h}={\varphi }_{1}{\nabla }^{d}{v}_{i, t+h-1}+\cdots +{\varphi }_{p}{\nabla }^{d}{v}_{i, t+h-p}+{\varepsilon }_{t}-\\ {\theta }_{1}{\varepsilon }_{t+h-1}-\cdots -{\theta }_{q}{\varepsilon }_{t+h-q}\end{array} $

(5)

步骤3 利用SVD进行重建得到h步向前预测结果。

当得到时间模式的估计值后，利用已经存储的奇异值和对应的左奇异向量重建时空矩阵，得到最终预测结果：

$ {\widehat{\boldsymbol{Y}}}_{t+h}={s}_{1}{\boldsymbol{u}}_{1}{\widehat{\boldsymbol{v}}}_{1, t+h}^{\mathrm{\text{'}}}+{s}_{2}{\boldsymbol{u}}_{2}{\widehat{\boldsymbol{v}}}_{2, t+h}^{\mathrm{\text{'}}}+\cdots +{s}_{r}{\boldsymbol{u}}_{r}{\widehat{\boldsymbol{v}}}_{r, t+h}^{\mathrm{\text{'}}} $

(6)

2.3 模型优化

模型优化包括奇异值选择和ARIMA模型参数选择过程。针对奇异值选择，不同的奇异值个数重建的矩阵精度不同，通常情况下，利用前几个较大奇异值即可基本重构原始时空矩阵，剩余奇异值可理解为由数据波动形成的噪音。本文通过遍历的方式验证了不同的奇异值个数对最后效果的影响，最终设定对前2个奇异值对应的时间模式进行建模。针对时间模式ARIMA模型的构建，首先需要对时间模式的平稳性进行检验，若不平稳，需要将其转化为平稳模式并在后续模型中逆推回真实预测结果；当数据平稳性检验通过后，利用ARIMA模型对平稳时间模式进行建模，并利用ACF和PACF图^[25]确定ARIMA模型中的p、d和q参数取值；最后通过交叉验证以及信息准则AIC、BIC^[26]对模型有效性进行检验并选择最优模型，在检验通过后，利用得到的ARIMA模型完成预测。本文所提ST-SVD算法描述如算法1所示，算法流程如图 1所示。

算法1 ST-SVD算法

输入  时空矩阵$ {\boldsymbol{Y}}_{D\times (t-1)} $，预测步长l

输出 D个基站的l步预测值

1.给定奇异值分解的阈值thre（本文设定为85%）

2.对$ {\mathrm{Y}}_{\mathrm{D}\times (\mathrm{t}-1)} $进行奇异值分解，得到U、S和V

3.利用阈值thre（85%）选择U、S和V中的有用成分

4.通过ARIMA模型对选定的V进行建模并根据AIC、BIC选定最优ARIMA模型

5.利用最优ARIMA模型进行l步预测

6.将预测结果与选定的U、S进行乘积，重建D个基站的预测值

3 案例分析 3.1 数据集描述

本文利用中国某大型城市的2 333个基站在216 h（共9天）内的流量数据对所提ST-SVD算法进行验证，数据的采集频率为1次/h。图 2所示为2 333个基站的相对位置布局，经纬度已经过处理，表 1所示为其中5个基站在13 h内的流量数据示例。图 3所示为3个基站在216 h内的流量变化情况，从图 3可以看出，基站3具有较明显的9个峰，表明基站流量的变化基本以一天为周期，虽然另外2个基站中基站1也存在较类似的峰值，但是两者的整体变化有较大差异。本文将216 h内的流量数据拆分成训练集和测试集，训练集包含前160 h的数据，测试集包含剩余56 h的数据。

3.2 ST-SVD模型构建

ST-SVD模型构建步骤如下：

步骤1 通过奇异值分解对中心化后的时空矩阵进行分解。在本案例中，时空矩阵$ \boldsymbol{Y} $的大小是2 333$ \times $160，在进行数据预处理（异常值处理、平稳性处理）之后，通过对处理后的$ \boldsymbol{Y} $进行奇异值分解，得到左奇异矩阵、右奇异矩阵和奇异值。图 4所示为截取的时空矩阵$ \boldsymbol{Y} $的右奇异矩阵，即时空序列的时间模式，此处截取了一天内每个小时之间的相关性情况，黄色区域表明相关性较强（彩色效果见《计算机工程》官网HTML版），从图 4可以看出，时间相关性具有较明显的周期模式，并且可预测性较强。图 5所示为时空矩阵$ \boldsymbol{Y} $中不同地点的皮尔逊相关系数与距离之间的关系，从图 5可以看出，针对该区域基站流量的时空数据，距离越近相关性越大的假设并不成立。

图 6所示为排序后的奇异值，一般而言，前几个奇异值即可涵盖大部分信息。从图 6可以看出，前2个奇异值占据了奇异值之和的89%，因此，本文分别构建一个奇异值的重建算法ST-SVD（1）和两个奇异值的重建算法ST-SVD（2）。

图 7、图 8分别对应前2个奇异值的左奇异矩阵（空间模式）和右奇异矩阵（时间模式）。从中可以看出，空间模式较为复杂，没有明显规律，但是时间模式显示出明显的周期性和可预测性。因此，本文利用ARIMA模型分别对2个时间序列进行建模并预测。

步骤2 通过AIC和BIC信息准则（表 2）选择ARIMA模型中的p、q、d参数并得到下述结果：

$ \begin{array}{l}{v}_{1, t}=-0.072\mathrm{ }6+1.635\mathrm{ }1\mathrm{ }{v}_{1, t-1}-0.531\mathrm{ }9\mathrm{ }{v}_{1, t-2}-\\ 0.196\mathrm{ }4\mathrm{ }{v}_{1, t-3}+0.150\mathrm{ }9\mathrm{ }{\varepsilon }_{t-1}-\\ 0.545\mathrm{ }0\mathrm{ }{\varepsilon }_{t-2}-0.369\mathrm{ }1\mathrm{ }{\varepsilon }_{t-3}\end{array} $

(7)

$ \begin{array}{l}{v}_{2, t}=1.263\mathrm{ }5\mathrm{ }{v}_{2, t-1}+0.263\mathrm{ }5\mathrm{ }{v}_{2, t-2}+\\ 0.721\mathrm{ }1\mathrm{ }{\varepsilon }_{t-1}-0.721\mathrm{ }1\mathrm{ }{\varepsilon }_{t-2}\end{array} $

(8)

RMSE和MAE^[27]的计算公式分别如下：

$ \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{DT}\sum\limits_{i=1}^{D}\sum\limits_{t=1}^{T}{\left({y}_{it}-{\widehat{y}}_{it}\right)}^{2}} $

(9)

$ \mathrm{M}\mathrm{A}\mathrm{E}=\frac{1}{DT}\sum\limits_{i=1}^{D}\sum\limits_{t=1}^{T}\left|{y}_{it}-{\widehat{y}}_{it}\right| $

(10)

步骤3 利用奇异值分解进行重建得到h步向前预测结果。在得到t时刻时间模式的预测值$ {v}_{1, t} $和$ {v}_{2, t} $后，即可利用式（1）结合左奇异矩阵和奇异值重建时空矩阵，得到t时刻空间各个位置的预测值。考虑到流量的周期性一般为一天，因此，本文利用ARIMA模型分别向前1步、向前6步、向前12步和向前24步进行预测，并利用Bootstrap^[28]从2 333个基站中抽取不同的样本量，从100次实验中取均值作为最终结果，以评估算法在整个周期内不同预测长度下的准确度和预测性能。

3.3 模型性能比较

本文将ST-SVD（1）、ST-SVD（2）与现有常用的ARIMA、Lasso-VAR、LSTM和STARMA 4种模型进行对比。其中，ARIMA模型并不是时空序列模型，但是在不考虑空间观测点的相关性时时空序列变成独立的多个时间序列，可以分别利用ARIMA进行建模预测。ARIMA模型时间成本极高，但是可作为一种基线模型进行对比。Lasso-VAR是带有Lasso正则化约束的VAR模型，其认为时空模型是时间序列模型加空间维度，即增加一维，然后通过传统的VAR模型并添加Lasso正则化来降低过拟合风险。LSTM是递归神经网络的变体，适用于时间序列，其与VAR类似，将时空数据集的空间维度叠加到时间序列中进行训练预测。STARMA模型是经典的时空分析模型，本文采用欧氏距离定义模型中的空间权重矩阵。实验过程中使用的软件、软件依赖包信息以及模型关键参数如表 3所示。

利用10个、20个、50个和100个基站160 h内的数据分别对上述6种模型进行训练，并给出向前1步、6步、12步和24步的预测结果，利用常见的预测精度指标——均方根误差RMSE和绝对值误差MAE对预测性能进行评估。由于本文案例中共有2 333个基站，为了提高性能评估的准确性并降低方差，通过Bootstrap在2 333个基站中随机选取上述10个、20个、50个和100个基站100次，并对100次的实验结果取平均值以作为最终的性能评估结果。

表 4所示为上述6种模型向前1步的部分预测结果，加粗数字为最优预测结果，括号中的百分数表示预测百分比误差，计算公式如式（11）所示，$ \widehat{y} $表示预测值，$ y $表示真实值。

$ \mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}=\frac{\left|\widehat{y}-y\right|}{y} $

(11)

从表 4可以看出：LSTM模型的预测精度最差，原因是其数据量过少，模型欠拟合，这表明神经网络模型需要足够多的样本来提高精度；STARMA模型优于不添加空间信息的ARIMA模型；ST-SVD的2种模型相较于其他4种模型预测准确率更优。具体地，利用2个奇异值的ST-SVD（2）模型的预测误差约为0.13，ARIMA、Lasso-VAR、LSTM和STARMA的误差分别约为0.22、0.21、0.92和0.19。ST-SVD（1）和ST-SVD（2）明显优于其他4种对比模型且ST-SVD（2）优于ST-SVD（1）。

表 5所示为上述6种模型在4种不同基站个数以及4种不同预测步长情况下的RMSE，括号中为MAE。从表 5可以看出，ST-SVD模型的性能明显优于其余4种对比模型，而且ST-SVD（2）的重构结果稍优于ST-SVD（1）的重构结果。从图 9、图 10可以直观地看出，2种ST-SVD模型的误差低于其余4种对比模型。

4 结束语

时空序列模型STARMA通过构建空间权重矩阵来表征数据的空间相关模式，但是空间权重的构建大多依赖距离等主观性因素，导致STARMA难以适用于多数数据集。本文建立一种新的时空序列模型ST-SVD，其利用SVD技术对时空数据集的时间模式和空间模式进行自动分解，通过ARIMA模型拟合时间模式并建模预测，最终重建出时空预测结果。ST-SVD模型不需要对数据集的空间结构进行假设，只需对时间序列实现建模，大幅降低了问题复杂度和模型训练成本。实验结果表明，ST-SVD模型的预测效果优于LSTM、STARMA等时空序列模型。但是，本文研究尚存在一定不足，一是ST-SVD认为空间模式是时不变的，即空间作用和时间作用相互独立，二是在奇异值分解后的时间序列建模中利用了较为传统的ARIMA模型，该模型是一种线性模型，无法捕捉到时间序列中的非线性模式。下一步将利用机器学习、深度学习等技术对时间模式进行建模，然后通过奇异值分解重建时空序列，以解决上述问题。