3.4 补充缺失值

3.4.1 问题与挑战

许多传感器被部署在物理世界中，产生了大量具有地理标记的时间序列数据。一般来说，如图3.16a所示，每个传感器都与一个部署位置相关联，在每个时间间隔（如每小时）产生一个读数。实际上，由于传感器或通信错误，读数经常在各个意想不到的时刻丢失。这些缺失的读数不仅影响实时监控，还影响了进一步数据分析的性能。

图3.16b展示了 m 个传感器在 n 个连续时间戳的读数，这些读数以矩阵的形式存储，其中一行代表一个传感器，一列代表一个时间戳。元素 v _ij 指的是第 i 个传感器在第 j 个时间戳的读数。如果一个传感器在某个时间间隔内的读数没有成功接收，就会产生一个空白元素。现在，填补传感器缺失值的问题可以转换为在矩阵中填补空白元素的问题。

在一系列地理感知时间序列数据中填补缺失的读数，存在两个挑战。

1.可能缺失任意传感器和时间戳的读数。在一些极端情况下，我们可能会连续丢失来自一个传感器（例如，图3.16a中的 s ₂ ）的读数，或者同时丢失所有传感器在某个（或多个）时间戳（如图3.16a中的 t ₂ 所示）的读数。我们称这些极端情况为块缺失。对于现有模型来说，处理块缺失问题非常困难，因为我们可能无法为模型找到稳定的输入。例如，非负矩阵分解（NMF）无法处理矩阵中某列或某行数据完全缺失的情况。

2.受到多个复杂因素的影响，传感器读数随位置和时间显著且非线性地变化。首先，距离较近的传感器的读数并不总是比距离较远的传感器更相似。如图3.17a所示，对于地理欧几里得距离， s ₁ 比 s ₃ 更接近 s ₂ 。然而，如图3.17b所示， s ₂ 的空气质量读数与 s ₃ 比与 s ₁ 更相似。原因是 s ₂ 和 s ₃ 位于具有相似地理背景的两个区域，例如POI和交通模式相似，而 s ₁ 部署在一个森林中，并且它与 s ₂ 之间有一个湖泊。这些情况违反了地理学的第一定律，降低了某些基于插值的模型的准确性。其次，传感器读数随时间波动极大，有时会出现突然变化。如图3.16b所示， s ₂ 在第31个时间戳的读数在两小时内下降了两百多。这种突然变化实际上对实时监控和进一步数据分析非常重要，但现有的平滑或插值方法无法很好地处理这种情况。

为了解决这个问题，人们从不同角度提出了许多解决方案。例如，在图3.16b所示的矩阵中， v _{2， j} 和 v _{1， j +1} 是缺失的。可以根据其空间邻域（如 s ₁ 和 s ₃ ）估计 v _{2， j} 的读数，我们称之为 空间视角 。 v _{2， j} 也可以根据相邻时间戳的读数（如 t _{j -1} 和 t _{j +1} ）等进行估计，我们称之为 时间视角 。我们还可以使用不同时间长度的数据进行估计，实现局部和全局视角。例如，我们考虑一个局部数据矩阵中从 t _{j -2} 至 t _{j +2} 的 v _{2， j} 的相邻读数，这被视为 局部视角 ，一个非常长的时间周期内（如从 t ₁ 到 t _n ）的读数，这被视为 全局视角 。局部视角捕捉瞬时变化，而全局视角代表长期模式。

3.4.2 空间模型

这类模型通过考虑传感器的邻域来估计其在时间间隔内的缺失值。三种广泛使用的方法是反距离加权（IDW） ^[39] 、线性估计和克里金法。

3.4.2.1 反距离加权

这个模型根据与目标传感器的距离为每个可用的在地理空间相邻的传感器的读数分配权重，并使用公式（3.14）聚合这些权重以生成预测值 。

其中 d _i 是候选传感器 s _i 与目标传感器之间的空间距离， α 是一个正的幂参数，它通过 控制传感器权重的衰减速率。 为距离较近的传感器的读数分配较大权重，较大的 α 表示权重随距离衰减得更快。

图3.18使用两个数据集（某城市2014年5月至2015年5月的空气质量数据和气象数据）的统计信息来阐述IDW。在这里，我们计算了同一时间戳两个任意传感器读数之间的比。在两个数据集中，随着两个传感器之间的距离增加，比下降。这遵循地理学的第一定律 ^[46] ：一切事物都与其他事物相关，但近的事物比远的事物更相关。这可以被视为地理感知数据的经验空间相关性。

在这里，我们使用图3.16所示的运行示例来展示IDW。假设两个传感器 s ₁ 和 s ₃ 分别距离 s ₂ 约6.5km和8.5km。我们的目标是使用 s ₁ 和 s ₃ 在时间间隔 t _j 的读数（即185和59）来填充 s ₂ 在 t _j 的缺失值。如果设置 α =1，那么两个传感器的权重分别是1/6.5和1/8.5。根据公式（3.14），我们计算预测值 v ^ _gs =（185/6.5+59/8.5）/（1/6.5+1/8.5）=130.4。

3.4.2.2 线性估计

在未观测到的位置 x ₀ 处对量 Z : R ⁿ → R 的空间推断或估计，是通过观测值 z _i = Z （ x _i ）和权重 ω _i （ x ₀ ）的线性组合来计算的，其中 i =0，1，2，…， N ：

权重 ω _i （ x ₀ ）可以通过最小化估计值的均方误差及其真实值 来学习。这是机器学习中的一种线性回归模型。我们将在第7章中详细介绍线性回归。

Pan和Li ^[42] 提出了一个基于 k 最近邻的算法，名为AKE，用于估计传感器的缺失值。AKE首先采用线性回归模型根据最近接收的数据计算每对传感器之间的时间依赖空间相关性。在学习过程结束后，生成一个判定系数，描述学习到的线性函数与接收的数据之间的拟合度。判定系数随后被用作对应传感器 Z （ x _i ）的权重 ω _i （ x ₀ ），在类似于公式（3.11）的加权平均函数中。尽管权重在不同时间间隔动态更新，但不同时间间隔之间的依赖性并未被考虑。

3.4.2.3 克里金法

在统计学中， 克里金法 是一种插值方法，插值的值由具有先验协方差的高斯过程建模。基于对先验的合理假设，克里金法为缺失值提供了最佳的线性无偏预测。

假设有 n +1个点，其索引从0到 n 。第 i 个点表示为 p _i ，其值 z _i 由变量 Z _i 生成。然后我们可以通过以下方式用观察到的值 z ₁ ， z ₂ ，…， z _n 的线性组合来估计未观测点 p ₀ 的值 ：

其中 λ ₁ ， λ ₂ ，…， λ _n 是对应于不同点的权重。在确定这些权重时，克里金法不仅考虑了两点之间的成对距离，还考虑了两点随机变量之间的相关性。一对变量（ Z _i ， Z _j ）之间的相关性由协方差矩阵定义，记作Cov（ Z _i ， Z _j ）。基于对这些变量的假设，已经提出了不同类型的克里金法，例如普通克里金法 ^[27] 、通用克里金法 ^[2] 、IRFk克里金法 ^[40] 和指标克里金法 ^[45] 等。

在本节中，我们重点介绍普通克里金法，它假设所有点的变量的期望值共享相同的常数 c ，尽管 c 的值可能是未知的，如公式（3.17）所示：

对 的估计误差定义为公式（3.18）：

然后，算法会找到参数 λ ₁ 、 λ ₂ 、…、 λ _n ，在 E （ ε ₀ ）=0的条件下使Var（ ε ₀ ）最小，其定义如公式（3.19）和公式（3.20）：

那么优化问题就变成了

通过添加一个拉格朗日乘子 ϕ ，可以求解如下：

实际上，计算Cov（ Z _i ， Z _j ）并不容易。另一种方法是计算变异函数 γ （ Z _i ， Z _j ），定义为公式（3.23）：

假设所有变量具有相同的方差 σ ² ，我们推导出以下关系：

γ （ Z _i ， Z _j ）= σ ² -Cov（ Z _i ， Z _j ）

因此，

通过在公式（3.21）中将Cov（ Z _i ， Z _j ）替换为 σ ² -γ （ Z _i ， Z _j ），普通克里金法方法进一步定义为如下：

最后，通过添加一个拉格朗日乘子 ϕ ，可以求解如下：

一旦确定了参数 λ ₁ ，…， λ _n ，我们就可以根据公式（3.16）计算 。

3.4.3 时间模型

这类方法根据传感器在其他时间间隔的自身读数来推断传感器在某个时间间隔的缺失值。这类方法可以分为两组：基于非特征的方法和基于特征的方法。第一组方法，如简单指数平滑（SES） ^[19] 、自回归移动平均（ARMA） ^[4] 和SARIMA ^[26] ，在推断缺失值时只考虑传感器的读数。Ceylan等人 ^[10] 对这些方法在时间序列中补充缺失值进行了比较。第二组方法，如图模型和回归模型 ^[19，34] ，利用特征函数，考虑了随时间变化的读数的相关性。由于第二组方法将在第7章中详细阐述，因此我们在此之后将重点介绍第一组方法。

3.4.3.1 简单指数平滑

SES通常用于时间序列领域，作为一种指数移动平均模型，正式定义为：

其中 t 是候选读数 v _j 和目标读数之间的时间间隔， β 是一个位于（0，1）之间的平滑参数。一般来说 β （1 -β ） ^{t -1} 会给最近的读数比远处的读数更大的权重，而较小的 β 表示权重随时间间隔更大衰减得更慢。

传统的SES只使用目标时间戳的前驱作为输入。Yi等人 ^[51] 通过使用目标时间戳的前驱和后继来扩展SES。给定一个目标时间戳，改进的SES会给同一传感器的每次读数赋予一个权重 β （1 -β ） ^{t -1} ，公式（3.28）通过根据归一化权重来计算 ：

在实际应用中，选择位于时间阈值（如12h）内的读数，因为远期的读数并不是非常有用。SES模型的灵感来自从时间序列数据观察的结果。图3.19展示了使用与图3.18相同的空气质量数据和气象数据，同一个传感器在不同时间戳的两个任意读数之间的比。图3.19中的两条曲线随着时间间隔的增加而下降，显示出时间序列中的经验时间相关性，也就是说近期时间戳的读数比远期时间戳的读数更为相关。

我们通过图3.16所示的例子来不断展示SES的使用。假设我们旨在使用四个相邻时间戳（ t _{j -2} ， t _{j -1} ， t _{j +1} ， t _{j +2} ）的读数来补充时间间隔 t _j 处的缺失值 s ₂ ，这四个时间戳的值分别是173，136，146，185。如果设定 β =0.5，那么这四个时间戳的权重分别是0.25，0.5，0.5，0.25。因此，最终结果是：

3.4.3.2 ARMA

在时间序列的统计分析中，ARMA模型描述了（弱）平稳随机过程，该过程用两个多项式来表示：一个用于自回归AR（ p ），另一个用于移动平均MA（ q ） ^[4] 。

AR（ p ）指的是阶数为 p 的自回归模型，写作：

其中 X _t 表示时间序列在第 t 个时间间隔的读数， φ ₁ ， φ ₂ ，…， φ _i 是参数， c 是一个常数，随机变量 ε _t 是白噪声。

MA（ q ）指的是阶数为 q 的移动平均模型：

其中 θ ₁ ， θ ₂ ，…， θ _i 是模型的参数， μ 是 X _t 的期望值（通常假设等于0）， ε _t ， ε _{t -1} ，…是白噪声误差项。

ARMA（ p ， q ）指的是包含 p 个自回归项和 q 个移动平均项的模型。这个模型包括了AR（ p ）和MA（ q ）模型，

滞后算子多项式记法 L ⁱ X _t = X _t-i 定义了 p 阶AR滞后算子多项式 φ （ L ）=（1 -φ ₁ L-φ ₂ L ² -… -φ _p L ^p ）和 q 阶MA滞后算子多项式 θ （ L ）=（1+ θ ₁ L + θ ₂ L ² +…+ θ _q L ^q ），我们可以将ARMA（ p ， q ）模型写成：

在实现ARMA模型时，我们要么使用默认参数 φ _i 和 θ _i ，要么预先定义这些参数。通常假设 ε _t 是服从独立同分布的随机变量，它们是从均值为0的正态分布中抽取的样本。

3.4.3.3 ARIMA和SARIMA

为了包含更现实的动态，特别是均值非平稳性和季节性行为，已经提出了许多ARMA模型的变体，包括自回归差分移动平均（ARIMA）模型和季节性自回归差分移动平均（SARIMA）模型。ARIMA中的AR部分表示感兴趣的演化变量回归于其先前的值。MA部分表示回归误差实际上是一系列在同一时期和过去不同时间发生的误差项的线性组合。I（代表“差分”）表示数据值已经被它们的值与先前值的差值替代（并且这个差分过程可能已经执行了多次）。这些特征的目的都是使模型尽可能好地拟合数据。

非季节性ARIMA模型通常表示为ARIMA（ p ， d ， q ），其中参数 p 、 d 和 q 是非负整数， p 是自回归模型的阶数（时间滞后数量）， d 是差分的程度（数据过去值被减去的次数）， q 是移动平均模型的阶数。SARIMA模型通常表示为ARIMA（ p ， d ， q ）（ P ， D ， Q ） m ，其中 m 指的是每个季节中的周期数，大写的 P 、 D 、 Q 分别指的是ARIMA模型季节部分的自回归、差分和移动平均项 ^[26] 。

3.4.4 时空模型

这类模型在估算缺失值时同时考虑了不同地点之间的空间相关性以及不同时间间隔之间的时间依赖性，共有三种不同的组合类型。

3.4.4.1 基于协同过滤的方法

这类方法通过矩阵来适应一段时间内的传感器数据，如图3.16b所示，其中一行代表一个传感器，一列代表一个时间间隔，一个元素存储了特定传感器在特定时间间隔的值。现在，传感器数据补充问题可以转换为推断缺失元素值的问题。两行之间的相似性表示两个传感器之间的空间相关性，两列之间的相似性表示两个时间间隔之间的时间相关性。解决推断问题的经典模型称为 协同过滤 ，这将在第8章中详细介绍。

Li等人 ^[35] 提出了两种基于矩阵分解的模型，包括STR-MF和MTR-MF，以补充传感器的缺失值，同时考虑了不同传感器之间的相关性以及不同时间间隔之间的相关性。具体来说，STR-MF通过在矩阵分解模型的损失函数中添加一个时间正则化项和一个空间接近项，来将时空相似性纳入矩阵分解模型中。对于具有多个传感器（例如温度和湿度）的位置，MTR-MF分别用矩阵容纳每种类型传感器的数据。然后，这些矩阵被集体分解，在损失函数中共享相同的潜在时间空间。此外，Li等人将MRT-MF扩展到三维张量，其中三个维度分别是位置、时间和传感器。随后，提出了一种带有时间正则化项的张量补全方法来估计张量中缺失的项。

3.4.4.2 基于多视图的方法

Yi等人 ^[51] 提出了一种时空多视图学习（ST-MVL）方法，用于集体补充地理感知时间序列数据集中的缺失读数，考虑了同一系列中不同时间戳之间的时间相关性和不同时间序列之间的空间相关性。如图3.20所示，ST-MVL包含四个视图：IDW、SES、基于用户的协同过滤（UCF）和基于物品的协同过滤（ICF）。然后，将这四个视图聚合以生成对缺失读数的最终估计。

IDW根据传感器空间邻域的读数估计传感器的缺失读数。SES基于同一传感器在其他时间戳的读数来估计传感器的缺失读数。由于IDW和SES实际上是基于长时间数据得出的经验模型，它们分别表示对缺失读数的全局空间视图和全局时间视图。

UCF仅根据传感器最近读数与其空间邻域读数之间的局部相似性来估计缺失读数，其中传感器被视为一个用户。同样，ICF基于不同时间戳最近读数之间的局部相似性来估计缺失读数，其中时间戳表示一个物品。由于UCF和ICF只考虑空间和时间角度的局部相似性，它们分别代表局部空间视图和局部时间视图。

为了利用不同视图的优势，我们提出了一种多视图学习算法，该算法根据公式（3.33）找到不同视图预测的线性组合，使得平方误差最小：

其中 b 是一个残差，而 w _i （ i =1，2，3，4）是分配给每个视图的权重。算法1展示了ST-MVL的过程。当一个数据集遇到块缺失问题时，ICF和UCF无法很好地工作，ST-MVL利用IDW和SES为那些缺失的项生成一个初始值（见图3.21中的第3行）。然后，ST-MVL分别使用ICF、UCF、IDW和SES预测每个缺失的项（第4～9行），根据线性核多视图学习框架[见第10行和方程（3.29）]组合这四个预测。分别针对每个传感器训练模型，最小化预测值和基准真实值之间的线性最小二乘误差 ^[33] 。