Traffic Data Imputation Based on Graph Regularization and Schatten-p Norm Minimization
-
摘要:
为充分利用交通数据低秩特性与局部近邻关系，准确恢复交通数据采集系统中的缺失数据，首先，应用基于核范数的低秩矩阵补全模型对交通数据矩阵进行预插补，以获得缺失值的初始估计，基于此，构建表征数据局部近邻结构的图模型；然后，提出融合图正则化和Schatten-
p 范数最小化的交通数据缺失值恢复模型；进一步，提出基于交替方向乘子框架的优化算法，求解缺失值恢复的最优化问题，得到最终的数据恢复结果；最后，用实际的高速公路交通流量和速度数据比较多种方法的恢复误差，同时给出所提方法的参数敏感性分析. 实验结果表明：在完全随机缺失、随机缺失和混合缺失模式下，缺失率为10% ~ 50%时，相比于局部最小二乘、概率主成分分析和低秩矩阵补全等方法，基于图正则化和Schatten-p 范数最小化的算法恢复误差降低了3.02% ~ 28.49%.-
关键词:
- 智能交通 /
- 数据恢复 /
- Schatten-p范数 /
- 交通数据 /
- 图正则化
Abstract:To make full use of the low-rank characteristics and local neighbor relationship of the traffic data, and accurately recover the missing data in traffic data acquisition system, firstly, the traffic data matrix is pre-interpolated by the low-rank matrix completion model based on kernel norm to obtain the initial estimate of the missing data. Based on this, a graph model that characterizes the local neighbor structure of the data is constructed. Then, a missing value imputation model combining graph regularization and Schatten-
p norm minimization is proposed. Furthermore, an optimization algorithm based on alternating direction multiplier framework is proposed to solve the optimization of missing value imputation, so as to obtain the final imputation result. Finally, the real expressway traffic volume and speed data are used to compare the imputation errors of several methods, and the parameter sensitivity of the proposed method is analyzed. The experimental results show that compared with local least squares, probabilistic principal component analysis and low-rank matrix completion, the proposed method reduces the error of traffic data imputation by 3.02%−28.49% when the missing rate is 10%−50% in missing completely at random mode, missing at random mode and mixed missing mode.-
Key words:
- intelligent transportation /
- data imputation /
- Schatten-p norm /
- traffic data /
- graph regularization
-
准确、完整的交通数据是实现交通大数据分析与挖掘的基础. 然而，实际交通环境中，由于检测器故障、恶劣天气、通信中断等原因，采集的交通数据通常存在缺失值. 例如，北京交通流数据检测设备采集的数据中每天约有10%的数据缺失[1]. 传统的机器学习算法，如支持向量机、神经网络等无法处理不完整数据，严重地影响交通预测与分析智能算法的性能. 因此，数据缺失问题已成为交通大数据分析的关键挑战之一.
近年来，如何对交通数据中的缺失值进行准确恢复受到国内外学者的广泛关注. 由于路网的拓扑结构、人们出行行为的规律性和周期性等，同一路网中不同位置的检测器在不同检测时间所记录的数据具有较强的时空相关性，因此，充分挖掘交通数据内在的时空相关性是实现缺失数据准确恢复的关键. 目前，交通数据缺失值恢复的主要方法可分为三类：基于预测的方法、基于插值的方法和基于统计学习的方法. 基于预测的方法是利用历史数据的特性将缺失值作为预测的目标进行恢复；例如，Henrickson等[2]针对交通数据中的连续缺失现象，提出一种预测均值匹配多重恢复方法，即使是全天数据缺失，也能取得良好的效果；但这类方法基于历史数据建模，忽略了缺失值之后的数据所提供的信息，导致对信息利用不充分，同时，当数据变化波动较大时，恢复数据的准确性将显著降低. 基于插值的方法是利用缺失值前后邻近的数据或者缺失值所在的邻近检测器数据取加权平均作为缺失值的恢复值；孙玲等[3]通过研究缺失值与观测值的相关性，对不同相关性的数据给予不同的权重，通过加权平均作为插补值；基于插值的方法计算简单，但要求缺失数据前后的相似性较高，并且当缺失值周围数据也存在缺失时，这类方法将会出现较大误差. 基于统计学习的方法通过引入机器学习理论与算法[4-5]建模交通数据的时空关系，实现缺失值恢复；李林超等[6]提出一种基于多源数据融合的异质交通流数据修复方法，通过自编码器提取高维多源交通流数据的时空特征，并采用随机森林估计数据中的缺失值；这类方法主要是建立符合数据分布的模型，并不断调整模型参数，使之实现对数据的最优拟合，进而实现缺失值的恢复.
低秩矩阵补全（low-rank matrix completion，LRMC）、低秩张量补全（low-rank tensor completion，LRTC）算法被成功应用于交通数据缺失值恢复[7]. 由于交通数据的时空相关性，交通数据矩阵具有低秩结构，通过矩阵秩最小化实现对缺失数据矩阵的恢复. Chen等[8]将交通量数据表示为张量形式，并利用贝叶斯张量分解理论交替优化缺失值和概率模型参数. Li等[7]对4种LRTC和4种预测模型进行了研究，表明提高缺失值恢复精度可以改善交通预测模型的性能. 基于LRMC的缺失值恢复方法取得了较好的效果，然而，最小化秩函数是一个NP （non-deterministic polynomial） 难的非凸优化问题，难以在多项式时间内求解. 一些方法将秩函数凸松弛为矩阵核范数进行求解，最小化核范数是一个凸优化问题，具有全局最优解. 然而，凸松弛问题的最优解可能会偏离原问题的最优解. 同时，大多数基于LRMC的恢复方法只利用了数据矩阵的全局时空相关性，而没有充分利用数据的局部结构特性.
基于上述分析，本文提出一种基于图正则化和Schatten-p范数（Schatten-p norm and graph regularization，SPGR）的交通数据缺失值恢复方法. 首先，应用已有的数据缺失值恢复方法对缺失值进行初步估计，基于此，建立刻画交通数据局部近邻结构的图模型以反映不同样本之间的相似度；依据交通数据的低秩结构这一全局先验信息，应用Schatten-p范数逼近矩阵的秩函数，提出融合图正则化和Schatten-p范数最小化的缺失值恢复模型. 然后，基于交替方向乘子法（alternating direction method of multipliers，ADMM）框架设计优化算法，实现对模型的有效求解. 最后，通过在真实路网交通流量与速度数据上的实验对比分析验证该方法的有效性.
1. 问题描述
将交通数据（如交通流量、速度）表示为变量矩阵
${\boldsymbol{X}} = \left( {{{\boldsymbol{X}}_1},{{\boldsymbol{X}}_2}, \cdots ,{{\boldsymbol{X}}_i}, \cdots ,{{\boldsymbol{X}}_N}} \right) \in {{\bf{R}}^{M \times N}}$ . N为样本总数；M为检测器一天采集的交通数据个数，即一个样本中的数据量；${{\boldsymbol{X}}_i} = {\left( {\mathop x\nolimits_i \left( 1 \right),\mathop x\nolimits_i \left( 2 \right), \cdots,{x_i(q)},\cdots,\mathop x\nolimits_i \left( Q \right)} \right)^{\rm{T}}}$ ，$\mathop x\nolimits_i \left( q \right)$ 为第i个样本中第q个时刻的交通数据. 设$ \varOmega $ 表示X中已知元素的下标集合，即${{\boldsymbol{X}}_\varOmega }$ 为观测值集合，$ \overline \varOmega $ 表示X中缺失元素的下标集合，因此，交通数据恢复问题就是根据观测值集合${{\boldsymbol{X}}_\varOmega }$ 准确估计缺失值集合${{\boldsymbol{X}}_{\overline \varOmega }}$ .2. 基于SPGR的交通数据缺失值恢复算法
2.1 算法框架
本文提出的基于图正则化和Schatten-p范数最小化的交通数据缺失值恢复算法，算法流程如图1所示，具体包括以下步骤：
步骤1 应用已有的数据缺失值恢复方法（如LRMC）得到缺失值的初始估计
$\widetilde {\boldsymbol{X}} = ( {{\widetilde {\boldsymbol{X}}}_1},{{\widetilde {\boldsymbol{X}}}_2}, \cdots ,{{\widetilde {\boldsymbol{X}}}_i}, \cdots , {{\widetilde {\boldsymbol{X}}}_N} )$ ，其中，$\widetilde {\boldsymbol{X}}_i = \left( {{{\tilde x}_i}(q)} \right)$ ，当$\mathop x\nolimits_i \left( q \right)$ 已知时，$\mathop {\tilde x}\nolimits_i \left( q \right){\text{ = }}\mathop x\nolimits_i \left( q \right)$ .步骤2 根据初始估计的结果，计算任意两个样本之间的距离，表示为矩阵
${\boldsymbol{D}}{\text{ = }}\left( {{d_{ij}}} \right) \in \mathop {\bf{R}}\nolimits^{N \times N}$ ，其中，$ \mathop d\nolimits_{ij} $ 为第$ i $ 个样本与第$ j $ 个样本间的距离，$j = 1, 2, \cdots ,N$ ，通过加权欧几里得公式[9]计算，如式（1）所示.$$ \mathop d\nolimits_{ij} = \sqrt {\sum\nolimits_{q = 1}^Q {\mathop {\tilde \theta }\nolimits_q } \left( {\mathop {\tilde x}\nolimits_i \left( q \right) - \mathop {\tilde x}\nolimits_j \left( q \right)} \right)} , $$ (1) $$ \mathop {\tilde \theta }\nolimits_q {\text{ = }}\frac{{\mathop \theta \nolimits_q }}{{\displaystyle\sum\nolimits_{q = 1}^Q {\mathop \theta \nolimits_q } }} , $$ (2) $$ {{\displaystyle \theta }}_{q}\text=\left\{\begin{array}{l} 1\text{，}\quad{{\displaystyle x}}_{i}\left(q\right)\text{、}{{\displaystyle x}}_{j}\left(q\right)已知,\\ \alpha \text{，}\quad 其他,\end{array}\right. $$ (3) 式中：
$ {\theta _q} $ 为两样本在第q个时刻的差值权重；$ \mathop {\tilde \theta }\nolimits_q $ 为$ {\theta _q} $ 的归一化值；$ \alpha $ 为相对小的正数，本文取0.1.步骤3 对样本
${{\boldsymbol{X}}_i}$ ，根据距离矩阵${\boldsymbol{D}}$ 选出与之最接近的K个样本，表示为集合${N_{\rm{e}}}( i )$ . 然后构造邻接矩阵（权矩阵）${\boldsymbol{S}} = ( {{s_{ij}}} ) \in \mathop {\bf{R}}\nolimits^{N \times N}$ .$$ {s}_{ij}=\left\{\begin{array}{l}1,\quad j\in {{\displaystyle N}}_{{\rm{e}}}\left(i\right)\;或\;i\in {{\displaystyle N}}_{{\rm{e}}}\left(j\right),\\ 0,\quad 其他.\end{array}\right. $$ (4) 步骤4 基于X与
$ {\boldsymbol{S}} $ ，建立SPGR模型并求解，得到最终恢复值${\boldsymbol{Y}} = ( {{{\boldsymbol{Y}}_1},{{\boldsymbol{Y}}_2}, \cdots ,{{\boldsymbol{Y}}_N}} )$ .2.2 SPGR模型
如前所述，交通数据具有较强的时空相关性，因此，交通数据矩阵具有低秩结构[1]，可将缺失值恢复问题转化为矩阵秩最小化问题，如式（5）所示.
$$ \left\{ {\begin{array}{*{20}{l}} {\mathop {\min }\limits_{\boldsymbol{X}} } & {\rm{rank}}\;\boldsymbol{X}, \\ {\rm { s.t. }} & {\boldsymbol{P}}_{\varOmega}(\boldsymbol{X})={\boldsymbol{P}}_{\varOmega}(\boldsymbol{A}), \end{array}} \right. $$ (5) 式中：
${\boldsymbol{A}} \in {{\bf{R}}^{M \times N}}$ 为观测值矩阵；${{\boldsymbol{P}}_\varOmega }( {\text{•}} )$ 为$\varOmega$ 的投影运算，如式（6）所示.$$ {\left( {{P_\varOmega }\left( {\boldsymbol{X}} \right)} \right)_{ij}} = \left\{ {\begin{array}{*{20}{l}} {{X_{ij}},\quad i,j \in \varOmega ,} \\ {0,\quad\;\;\; i,j \notin \varOmega .} \end{array}} \right. $$ (6) 然而，由于秩函数是非凸且不连续，所以式（5）是一个NP难问题，难以求解. 为解决这一问题，用矩阵的Schatten-p （
$ {0 \lt p \lt \infty } $ ）范数逼近式（5）中的秩函数，具体表示为$$ {\left\| {\boldsymbol{X}} \right\|_{{\rm{sp}}}} = {\left( {\sum\nolimits_{t = 1}^{\min \left\{ {N,M} \right\}} {\sigma _t^p} } \right)^{\frac{1}{p}}} , $$ (7) 式中：
${\left\| {\boldsymbol{X}} \right\|_{{\rm{sp}}}}$ 为${\boldsymbol{X}}$ 的Schatten-p范数；$ {\sigma _t} $ 为${\boldsymbol{X}}$ 的第$ t $ 个奇异值.当 p = 1 时，Schatten-1范数即为核范数. 当p越趋近于0时，
${\boldsymbol{X}}$ 的Schatten-p范数越接近秩函数. 当p = 0 时，式（7）就是${\boldsymbol{X}}$ 的秩函数. 总的来说，与核范数相比，Schatten-p范数对秩函数的逼近能力更强，可以更精确刻画数据的低秩结构. 因此，本文构建基于Schatten-p范数最小化的缺失值恢复模型，如式（8）所示.$$ \mathop {\min }\limits_{\boldsymbol{X}} \left\| {\boldsymbol{X}} \right\|_{_{{\rm{sp}}}}^p + \left\| {{P_\varOmega }\left( {\boldsymbol{X}} \right) - {P_\varOmega }\left( {\boldsymbol{A}} \right)} \right\|_{\rm{F}}^2 , $$ (8) 式中：
$ \left\| {\text{•}} \right\|_{\rm{F}}^{} $ 为矩阵的F范数.同时，邻近的数据样本具有相似的特征，为更好利用这种信息，在进行恢复时，将所有邻域样本之间的距离限制在适当的范围内，以防止被恢复样本与邻域样本的差异过大. 基于上述分析，提出图正则化来刻画这种局部邻近结构：
$$ \sum\limits_{{i,j}} {{s_{ij}}} {\left\| {{x_i} - {x_j}} \right\|^2} = {\rm{tr}}\left( {{{\boldsymbol{X}}^{\rm{T}}}{\boldsymbol{LX}}} \right) , $$ (9) 式中：
${\boldsymbol{L}} = {\boldsymbol{E}} - {\boldsymbol{S}}$ ，为图的拉普拉斯矩阵；E 为对角元素是${E_{ii}} = \displaystyle\sum\nolimits_{j = 1}^N {{s_{ij}}}$ 的对角矩阵.结合式（8）、（9），SPGR的目标函数为
$$ \mathop {\min }\limits_{\boldsymbol{X}} \left\| {\boldsymbol{X}} \right\|_{{\rm{sp}}}^p + \left\| {{{\boldsymbol{P}}_\varOmega }\left( {\boldsymbol{X}} \right) - {{\boldsymbol{P}}_\varOmega }\left( {\boldsymbol{A}} \right)} \right\|_{\rm{F}}^2 + \lambda\; {\text{tr}}\left( {{{\boldsymbol{X}}^{\rm{T}}}{\boldsymbol{LX}}} \right) , $$ (10) 式中：
$ \lambda \gt 0 $ 为控制图正则化项的权重常数.2.3 基于ADMM的优化算法
本文采用交替方向乘子法（ADMM）框架寻求式（10）的最优解[10]. 首先，引入辅助变量
${\boldsymbol{W}}$ 和${\boldsymbol{Z}}$ ，将式（10）转化为式（11）所示等价问题.$$ \left\{ {\begin{array}{*{20}{l}} {\mathop {\min }\limits_{\boldsymbol{X}} \left\| {\boldsymbol{W}} \right\|_{{\text{sp}}}^p + \left\| {{P_\varOmega }\left( {\boldsymbol{Z}} \right) - {P_\varOmega }\left( A \right)} \right\|_{\text{F}}^2 + \lambda {\text{tr}}\left( {{{\boldsymbol{X}}^{\text{T}}}{\boldsymbol{LX}}} \right),} \\ {{\rm{s.t}}.\quad {\boldsymbol{X}} = {\boldsymbol{W}}{\text{ = }}{\boldsymbol{Z}}. } \end{array}} \right. $$ (11) 构造式（11）的增广拉格朗日函数为
$$ \begin{split} &\mathcal{L}(\boldsymbol{W}, \boldsymbol{Z}, \boldsymbol{X}, \boldsymbol{U}, \boldsymbol{V})=\|\boldsymbol{W}\|^{p}_{{\rm{sp}}}+\left\|P_{\varOmega}(\boldsymbol{Z})-P_{\varOmega}(\boldsymbol{A})\right\|_{\mathrm{F}}^{2} \;+\\ &\quad \lambda \operatorname{tr}\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{L} \boldsymbol{X}\right)+\operatorname{tr}\left(\boldsymbol{U}^{\mathrm{T}}(\boldsymbol{X}-\boldsymbol{W})\right) +\dfrac{\mu_{1}}{2}\|\boldsymbol{X}-\boldsymbol{W}\|_{\mathrm{F}}^{2}\;+\\ &\quad \operatorname{tr}\left(\boldsymbol{V}^{\mathrm{T}}(\boldsymbol{X}-\boldsymbol{Z})\right) +\dfrac{\mu_{2}}{2}\|\boldsymbol{X}-\boldsymbol{Z}\|_{\mathrm{F}}^{2},\\[-15pt] \end{split} $$ (12) 式中：
${\boldsymbol{U}}$ 、${\boldsymbol{V}}$ 为拉格朗日乘子；$ {\mu _1} \gt 0 $ ，$ {\mu _2} \gt 0 $ ，为惩罚系数.增广拉格朗日函数融合了罚函数法与拉格朗日乘子法的优点，依据ADMM框架对优化变量进行迭代求解. 一个变量进行优化时，固定其他变量，通过变量交替近似求解，实现算法结果的优化.
1） 固定变量Z 和X，
$$ \begin{split} &{{\boldsymbol{W}}_{l + 1}} = \arg \mathop {\min }\limits_{\boldsymbol{W}} \mathcal{L}\left( {{\boldsymbol{W}},{{\boldsymbol{Z}}_l},{{\boldsymbol{X}}_l},{{\boldsymbol{U}}_l},{{\boldsymbol{V}}_l}} \right)= \\ &\quad \arg \mathop {\min }\limits_{\boldsymbol{W}} \left\| W \right\|_{{\text{sp}}}^p +{\text{tr}}\left( {{\boldsymbol{U}}_l^{\rm{T}}\left( {{{\boldsymbol{X}}_l} - {\boldsymbol{W}}} \right)} \right)+ \dfrac{{{\mu _1}}}{2}\left\| {{{\boldsymbol{X}}_l} - {\boldsymbol{W}}} \right\|_{\rm{F}}^2 = \\ &\quad \arg \mathop {\min }\limits_{\boldsymbol{W}} \left\| {\boldsymbol{W}} \right\|_{{\rm{sp}}}^p+ \dfrac{1}{{{\mu _1}}}\left\| {{\boldsymbol{W}} - \left( {{{\boldsymbol{X}}_l} + \dfrac{1}{{{\mu _1}}}{{\boldsymbol{U}}_l}} \right)} \right\|_{\rm{F}}^2, \\[-12pt] \end{split} $$ (13) 式中：l为ADMM算法中的迭代次数；
${{\boldsymbol{W}}_{l}}$ 、$ {{\boldsymbol{Z}}_l} $ 、$ {{\boldsymbol{X}}_l} $ 、$ {{\boldsymbol{U}}_l} $ 、$ {{\boldsymbol{V}}_l} $ 分别为W、$ {\boldsymbol{Z}} $ 、$ {\boldsymbol{X}} $ 、$ {\boldsymbol{U}} $ 、$ {\boldsymbol{V}} $ 迭代l次后对应的数值.式（13）的最优解通过迭代算法[11]求出.
2） 固定变量W 和X，
$$ \begin{split} & {{\boldsymbol{Z}}_{l + 1}} = \arg \mathop {\min }\limits_{\boldsymbol{Z}} \mathcal{L}\left( {{{\boldsymbol{W}}_{l + 1}},{\boldsymbol{Z}},{{\boldsymbol{X}}_l},{{\boldsymbol{U}}_l},{{\boldsymbol{V}}_l}} \right)= \\ &\quad \arg \mathop {\min }\limits_{\boldsymbol{Z}}\left\| {{P_\varOmega }\left( {\boldsymbol{Z}} \right) - {P_\varOmega }\left( A{\boldsymbol{}} \right)} \right\|_{\rm{F}}^2\;+\\ &\quad {\text{tr}}\left( {{\boldsymbol{V}}_l^{\rm{T}}\left( {{{\boldsymbol{X}}_l} - {\boldsymbol{Z}}} \right)} \right) + \dfrac{{{\mu _2}}}{2}\left\| {{{\boldsymbol{X}}_l} - {\boldsymbol{Z}}} \right\|_{\rm{F}}^2. \\[-6pt] \end{split} $$ (14) 式（14）中对Z求导，并将导数设为0，得到
$$ \begin{split} &{{\boldsymbol{Z}}_{l + 1}}{\text{ = }}\frac{1}{{{\mu _2}}}{{\boldsymbol{V}}_l}{\text{ + }}\frac{2}{{{\mu _2}}}{P_\varOmega }\left( {\boldsymbol{A}} \right){\text{ + }}{{\boldsymbol{X}}_l}\;- \\ &\quad \frac{2}{{2 + {\mu _2}}}{P_\varOmega }\left( {\frac{1}{{{\mu _2}}}{{\boldsymbol{V}}_l}{\text{ + }}\frac{2}{{{\mu _2}}}{P_\varOmega }\left( {\boldsymbol{A}} \right){\text{ + }}{{\boldsymbol{X}}_l}} \right). \\[-10pt] \end{split} $$ (15) 3） 固定变量W 和Z，
$$ \begin{split} & {{\boldsymbol{X}}_{l + 1}} = \arg \mathop {\min }\limits_{\boldsymbol{X}} \mathcal{L}\left( {{{\boldsymbol{W}}_{l + 1}},{{\boldsymbol{Z}}_{l + 1}},{\boldsymbol{X}},{{\boldsymbol{U}}_l},{{\boldsymbol{V}}_l}} \right)= \\ &\quad \arg \mathop {\min }\limits_{\boldsymbol{X}} \lambda {\text{tr}}\left( {{{\boldsymbol{X}}^{\rm{T}}}{\boldsymbol{LX}}} \right)+ {\text{tr}}\left( {{\boldsymbol{U}}_l^{\rm{T}}\left( {{\boldsymbol{X}} - {{\boldsymbol{W}}_{l + 1}}} \right)} \right) + \\ &\quad \frac{{{\mu _1}}}{2}\left\| {{\boldsymbol{X}} - {{\boldsymbol{W}}_{l + 1}}} \right\|_{\rm{F}}^2+ {\text{tr}}\left( {{\boldsymbol{V}}_l^{\rm{T}}\left( {{\boldsymbol{X}} - {{\boldsymbol{Z}}_{l + 1}}} \right)} \right) + \\ &\quad \frac{{{\mu _2}}}{2}\left\| {{\boldsymbol{X}} - {{\boldsymbol{Z}}_{l + 1}}} \right\|_{\rm{F}}^2. \\[-10pt] \end{split} $$ (16) 式（16）中对X求导，并将导数设为0，得到
$$ \left( {2\lambda {\boldsymbol{L}} + {\mu _1}{\boldsymbol{I}}} \right){\boldsymbol{X}} + {\mu _2}{\boldsymbol{X}} = {\mu _1}{{\boldsymbol{W}}_{l + 1}} + {\mu _2}{{\boldsymbol{Z}}_{l + 1}} - {{\boldsymbol{U}}_l} - {{\boldsymbol{V}}_l}. $$ (17) 式（17）为Sylvester方程，对于方程
${\boldsymbol{AX }} + {\boldsymbol{XB}} = {\boldsymbol{C}}$ ，其解表示为${\boldsymbol{X}} = s\left( {{\boldsymbol{A}},{\boldsymbol{B}},{\boldsymbol{C}}} \right)$ [12]，$s\left( {\text{•}} \right)$ 为Sylvester方程的求解运算. 因此，$$ {{\boldsymbol{X}}_{l + 1}} = s\left( \begin{gathered} 2\lambda {\boldsymbol{L}} + {\mu _1}{\boldsymbol{I}},\;{\mu _2}{\boldsymbol{I}},\;{\mu _1}{{\boldsymbol{W}}_{l + 1}}{\text{ + }} {\mu _2}{{\boldsymbol{Z}}_{l + 1}} - {{\boldsymbol{U}}_l} - {{\boldsymbol{V}}_l} \\ \end{gathered} \right) . $$ (18) 4） 更新乘子Ul + 1 和Vl + 1，
$$ {{\boldsymbol{U}}_{l + 1}} = {{\boldsymbol{U}}_l} + {\mu _1}\left( {{{\boldsymbol{X}}_{l + 1}} - {{\boldsymbol{W}}_{l + 1}}} \right) , $$ (19) $$ {{\boldsymbol{V}}_{l + 1}} = {{\boldsymbol{V}}_l} + {\mu _2}\left( {{{\boldsymbol{X}}_{l + 1}} - {{\boldsymbol{Z}}_{l + 1}}} \right) . $$ (20) 变量W、Z和X按照上述规则迭代更新，直到算法收敛，收敛条件
$ \varepsilon $ 为$$ \varepsilon = \frac{{{{\left\| {{{\boldsymbol{X}}_{l + 1}} - {{\boldsymbol{X}}_l}} \right\|}_{\rm{F}}}}}{{\max \{ {{{\left\| {{X_{l + 1}}} \right\|}_{\rm{F}}},1\}} }} < r , $$ (21) 式中：r为很小的正数，本文取
$ 1 \times {10^{ - 5}} $ .当
$ \varepsilon < r $ 时，该算法收敛，停止迭代，输出结果.综上，求解式（11）的优化算法流程如图2所示.
3. 实验验证
3.1 交通数据与实验方案
为了评估SPGR方法的数据恢复性能，在真实的交通流量和交通速度数据上进行实验分析. 交通流量数据来自美国俄勒冈州波特兰市交通信息中心，从I205和I84州际公路构成的路网中选择40个只有极少缺失数据的检测站采集数据. 由于工作日的交通流量数据与周末和节假日的数据存在较大差异，因此，选取2015年中连续30个工作日的交通流数据进行研究. 传感器采样间隔为15 min，最终获取30 × 40 = 1200个样本，构造为96 × 1200的数据矩阵. 此外，交通速度数据集为中国广州两个月（2016年8月1日至2016年9月30日共61 d）内以10 min为间隔的7个路段（主要包括城市快速路和干线）的速度信息. 因此，得到7 × 61 = 421个样本，构造为144 × 421数据矩阵. 图3展示了8个传感器在同一天的交通数据变化情况.
图 3 同一天中不同传感器交通流和交通速度的变化情况Figure 3. Changes in traffic flow and speed from different sensors over the same day为反映交通数据缺失值的复杂分布，模拟3种常见的数据缺失模式：1） 完全随机缺失（missing completely at random，MCAR），缺失值独立于其他缺失数据或已知数据，表现为一组随机分布的孤立点；2） 随机缺失（missing at random，MAR），交通数据表现为连续缺失的现象，即缺失值的恢复依赖于相邻的缺失值；3） 混合缺失（mixture of MCAR and MAR，MIXED），MAR与MCAR的混合比例各为0.5. 图4为不同缺失模式的示例，图中每行表示一个流量样本，每列表示一个变量，黑色表示缺失值.
为综合比较不同数据恢复方法的有效性，将提出的SPGR模型、去除图正则化的SP模型（Schatten-p范数最小化，
$ \lambda = 0 $ ）与3种缺失值恢复方法：LRMC、概率主成分分析（probabilistic principal component analysis，PPCA）、局部最小二乘（local least squares，LLS）进行比较[1,13-14]. 这3种对比方法涵盖了数据缺失值恢复的主流技术，包括低秩矩阵补全、概率模型和回归模型.实验中，按照缺失模式和缺失比例模拟缺失值. 其中，缺失率
$ \delta $ 定义为缺失数据数量与总数据量之比，以0.1为步长将$ \delta $ 从0.1增加到0.5，研究不同缺失率对恢复性能的影响. 为衡量不同算法的恢复性能，采用缺失项的恢复值与真实值之间的均方根误差（RMSE，eRMSE）和平均绝对百分比误差（MAPE，eMAPE）表示，分别为$$ \mathop {{e}}\nolimits_{{\rm{RMSE}}} = \sqrt {\frac{1}{C}\mathop {\left| {{\boldsymbol{X}}_{\overline \varOmega }^{} - \widetilde {\boldsymbol{X}}_{\overline \varOmega }^{}} \right|}\nolimits^2 } , $$ (22) $$ \mathop {{e}}\nolimits_{{\rm{MAPE}}} = \frac{1}{C}\left| {\frac{{{\boldsymbol{X}}_{\overline \varOmega }^{} - \widetilde {\boldsymbol{X}}_{\overline \Omega }^{}}}{{{\boldsymbol{X}}_{\overline \varOmega }^{}}}} \right| \times 100{\text{%}} , $$ (23) 式中：
$ C $ 为缺失值的总数目；${\boldsymbol{X}}_{\overline \varOmega }^{}$ 和$\widetilde {\boldsymbol{X}}_{\overline \varOmega }^{}$ 分别为真实值和恢复值.RMSE和MAPE越小，算法的恢复性能越好. 为准确评估5种缺失值恢复方法在两种交通数据上的性能，减少随机性对实验结果的影响，每种实验均重复5次，取实验结果的误差平均值作为评价缺失值恢复方法性能的依据.
3.2 实验结果分析
表1 ~ 3分别列出了不同算法在MCAR、MAR和MIXED模式下的恢复误差，根据实验结果，可以得到以下结论：1） MCAR缺失模式下，各种算法的恢复误差最小；而在MAR缺失模式下，因缺失大量相关信息，导致数据恢复误差较大；此外，每种方法的恢复误差随着缺失率的增加而增加. 2） 在恢复性能方面，LRMC和PPCA方法不能很好地处理内部结构复杂的数据集，导致总体性能比其他方法差；当缺失率较低时，PPCA的性能较好，而当缺失率增加时，LRMC的性能优于PPCA；LLS在低缺失率下恢复性能较好，然而，当缺失率增加时，其恢复性能会迅速退化. 3） 本文提出的SPGR算法获得了更好的恢复性能；与LLS、PPCA、LRMC方法相比，在MCAR缺失模式下，RMSE降低了3.02% ~ 22.31%，在MAR缺失模式下，RMSE降低了3.23% ~ 28.49%，在MIXED缺失模式下，RMSE降低了3.05% ~ 21.56%；特别是当缺失率大于30%时，恢复误差降低率越大，相对于LLS算法误差降低率高达28.49%，表明该方法可以有效挖掘观测数据的内在关联，实现准确的缺失值恢复.
表 1 MCAR模式下不同算法的恢复误差Table 1. Imputation error of different algorithms in MCAR mode算法 δ 交通流量 交通速度 RMSE/（辆·
（15 min）−1）MAPE/% RMSE/
（km·h−1）MAPE/% LRMC 0.1 80.63 13.29 2.53 5.76 0.2 83.80 13.75 2.64 6.07 0.3 86.58 14.29 2.65 6.40 0.4 90.94 14.86 2.79 6.83 0.5 95.91 15.74 2.99 7.37 PPCA 0.1 79.06 13.00 2.41 5.82 0.2 82.98 13.80 2.53 6.08 0.3 85.47 14.67 2.72 6.27 0.4 89.69 15.55 2.81 6.85 0.5 95.83 17.11 2.98 6.95 LLS 0.1 76.62 13.26 3.15 7.40 0.2 82.26 14.21 3.28 7.74 0.3 89.94 15.81 3.43 8.13 0.4 99.01 17.61 3.67 8.82 0.5 113.01 19.92 4.06 9.76 SP 0.1 75.77 12.73 2.39 5.65 0.2 78.36 13.18 2.50 5.95 0.3 82.43 13.73 2.63 6.31 0.4 88.06 14.51 2.78 6.74 0.5 93.66 15.56 2.98 7.29 SPGR 0.1 74.56 12.67 2.38 5.64 0.2 76.96 13.17 2.50 5.95 0.3 80.93 13.59 2.63 6.30 0.4 86.24 14.38 2.76 6.70 0.5 92.40 15.44 2.95 7.19 表 2 MAR模式下不同算法的恢复误差Table 2. Imputation error of different algorithms in MAR mode算法 δ 交通流量 交通速度 RMSE/（辆·
（15 min）−1）MAPE/% RMSE/
（km·h−1）MAPE/% LRMC 0.1 98.31 14.62 3.38 8.42 0.2 100.92 15.14 3.45 8.68 0.3 103.11 15.72 3.48 8.77 0.4 106.22 16.30 3.55 8.98 0.5 109.81 17.03 3.63 9.15 PPCA 0.1 90.14 13.71 3.60 8.48 0.2 95.13 14.70 3.75 8.83 0.3 97.45 15.49 3.77 8.87 0.4 101.55 18.49 3.80 8.93 0.5 105.51 17.72 3.85 9.02 LLS 0.1 84.31 13.40 3.75 8.89 0.2 91.89 14.70 3.86 9.30 0.3 100.10 16.38 3.91 9.40 0.4 112.55 16.44 4.11 10.02 0.5 130.83 21.11 4.51 10.69 SP 0.1 88.08 13.64 3.26 8.17 0.2 90.99 14.09 3.37 8.42 0.3 94.50 14.96 3.43 8.60 0.4 99.86 15.83 3.53 8.86 0.5 105.22 16.77 3.60 9.00 SPGR 0.1 86.31 13.51 3.22 7.97 0.2 89.01 13.93 3.31 8.23 0.3 92.76 14.79 3.36 8.35 0.4 97.51 15.63 3.48 8.64 0.5 101.82 16.62 3.56 8.81 3.3 算法参数的影响
本文提出的SPGR模型涉及3个参数：Schatten-p范数的
$ p $ 值、K近邻（K-nearest neighbor, KNN）方法中的K值、图正则化的权重常数$ \lambda $ . p值决定了对矩阵秩函数的近似程度，$ K $ 值决定了用于重建每个样本的近邻数量，$ \lambda $ 控制基于局部邻近的图正则化的影响. 为得到模型参数的最优值，调整其中一个参数，固定另外两个参数，每次参数改变时记录实验结果. 以交通流量数据为例，图5给出了不同参数下，缺失值恢复误差RMSE的变化. 由图可知：不同缺失率下，p具有不同的最优值，由于Schatten-p范数比核范数（$ p = 1 $ ）更能逼近秩函数，从而获得更好的恢复结果；如果K值过大或过小，都会导致恢复精度较差. 这是因为K值太小，选择代表目标样本的相邻样本过少，导致用于缺失值恢复的可用信息不充分，K值过大，远离目标的样本将参与重建，也将降低恢复精度. 对于$ \lambda $ 的影响，也可以得出与K值类似的结论.表 3 MIXED模式下不同算法的恢复误差Table 3. Imputation error of different algorithms in MIXED mode算法 δ 交通流量 交通速度 RMSE/（辆·
（15 min）−1）MAPE/% RMSE/
（km·h−1）MAPE/% LRMC 0.1 89.57 13.56 2.97 7.10 0.2 92.73 14.31 2.99 7.28 0.3 95.18 14.70 3.10 7.59 0.4 99.00 15.61 3.19 7.86 0.5 103.12 16.24 2.99 8.24 PPCA 0.1 84.40 12.95 3.13 7.24 0.2 89.37 14.08 3.10 7.33 0.3 91.37 14.76 3.18 7.57 0.4 95.26 16.10 3.24 7.73 0.5 100.80 17.33 2.98 7.91 LLS 0.1 79.48 13.11 3.65 8.05 0.2 86.60 14.34 3.54 8.39 0.3 93.75 15.77 3.66 8.74 0.4 104.31 17.83 3.82 9.20 0.5 118.90 20.15 4.06 10.09 SP 0.1 82.10 12.85 2.92 6.95 0.2 85.54 13.62 3.02 7.10 0.3 88.72 14.34 3.10 7.49 0.4 94.10 15.35 3.18 7.79 0.5 99.60 16.18 3.31 8.15 SPGR 0.1 80.51 12.74 2.86 6.81 0.2 83.82 13.49 2.89 7.03 0.3 86.62 14.18 3.02 7.36 0.4 92.01 15.21 3.10 7.62 0.5 97.81 16.04 3.24 7.98 
图 5 SPGR模型在交通流量数据上的RMSE随参数$ p $ 、$ K $ 、$ \lambda $ 的变化Figure 5. RMSEs of SPGR model on traffic flow data varied with the parameters of$ p、K、\lambda $ 3.4 初始化影响
根据SPGR算法的第一步，为建立表征样本间相邻关系的图矩阵，需要选择一种已有的恢复方法对交通数据的缺失值进行初始估计. 本节将采用3种不同的初始化方法，即KNN、LLS、LRMC，进行缺失值的初始估计，进而研究其对SPGR算法性能的影响，实现敏感性分析. 以交通流数据为例，
$\delta \;{\text{ = }}\; 0.3$ 下的实验结果如表4所示. 可以看出，不同的初始化方法对SPGR的恢复性能影响较小，在其他缺失率下也可以观察到类似现象. 这验证了SPGR对初始值具有较好的鲁棒性.表 4 在交通流量数据上不同初始化方法对SPGR恢复误差的影响Table 4. Effect of different initialization methods on SPGR imputation error on traffic flow data辆/15 min 数据模式 KNN LLS LRMC MCAR 81.21 80.99 80.93 MAR 93.68 93.14 92.76 MIXED 87.74 87.48 86.62 4. 结 论
提出一种融合图正则化与Schatten-p范数最小化的交通数据缺失值恢复方法. 该方法采用Schatten-p范数逼近矩阵的秩函数，对数据的低秩先验信息进行约束. 通过实验分析，得到以下结论：
1） 将图正则化融入到数据恢复框架中，有利于更好地利用数据的局部邻近结构.
2） 基于真实的高速公路交通量和速度数据进行仿真实验表明，提出的方法相对于其他多种方法恢复误差降低了3.02%以上，特别是缺失率大于0.3时，误差降低率达到20%以上.
3） 在未来的工作中，将进一步研究交通数据在时序上的规律性，以提升缺失值恢复的精度.
-
图 3 同一天中不同传感器交通流和交通速度的变化情况
Figure 3. Changes in traffic flow and speed from different sensors over the same day
图 5 SPGR模型在交通流量数据上的RMSE随参数
$ p $ 、$ K $ 、$ \lambda $ 的变化Figure 5. RMSEs of SPGR model on traffic flow data varied with the parameters of
$ p、K、\lambda $ 表 1 MCAR模式下不同算法的恢复误差
Table 1. Imputation error of different algorithms in MCAR mode
算法 δ 交通流量 交通速度 RMSE/（辆·
（15 min）−1）MAPE/% RMSE/
（km·h−1）MAPE/% LRMC 0.1 80.63 13.29 2.53 5.76 0.2 83.80 13.75 2.64 6.07 0.3 86.58 14.29 2.65 6.40 0.4 90.94 14.86 2.79 6.83 0.5 95.91 15.74 2.99 7.37 PPCA 0.1 79.06 13.00 2.41 5.82 0.2 82.98 13.80 2.53 6.08 0.3 85.47 14.67 2.72 6.27 0.4 89.69 15.55 2.81 6.85 0.5 95.83 17.11 2.98 6.95 LLS 0.1 76.62 13.26 3.15 7.40 0.2 82.26 14.21 3.28 7.74 0.3 89.94 15.81 3.43 8.13 0.4 99.01 17.61 3.67 8.82 0.5 113.01 19.92 4.06 9.76 SP 0.1 75.77 12.73 2.39 5.65 0.2 78.36 13.18 2.50 5.95 0.3 82.43 13.73 2.63 6.31 0.4 88.06 14.51 2.78 6.74 0.5 93.66 15.56 2.98 7.29 SPGR 0.1 74.56 12.67 2.38 5.64 0.2 76.96 13.17 2.50 5.95 0.3 80.93 13.59 2.63 6.30 0.4 86.24 14.38 2.76 6.70 0.5 92.40 15.44 2.95 7.19 
下载: 导出CSV
表 2 MAR模式下不同算法的恢复误差
Table 2. Imputation error of different algorithms in MAR mode
算法 δ 交通流量 交通速度 RMSE/（辆·
（15 min）−1）MAPE/% RMSE/
（km·h−1）MAPE/% LRMC 0.1 98.31 14.62 3.38 8.42 0.2 100.92 15.14 3.45 8.68 0.3 103.11 15.72 3.48 8.77 0.4 106.22 16.30 3.55 8.98 0.5 109.81 17.03 3.63 9.15 PPCA 0.1 90.14 13.71 3.60 8.48 0.2 95.13 14.70 3.75 8.83 0.3 97.45 15.49 3.77 8.87 0.4 101.55 18.49 3.80 8.93 0.5 105.51 17.72 3.85 9.02 LLS 0.1 84.31 13.40 3.75 8.89 0.2 91.89 14.70 3.86 9.30 0.3 100.10 16.38 3.91 9.40 0.4 112.55 16.44 4.11 10.02 0.5 130.83 21.11 4.51 10.69 SP 0.1 88.08 13.64 3.26 8.17 0.2 90.99 14.09 3.37 8.42 0.3 94.50 14.96 3.43 8.60 0.4 99.86 15.83 3.53 8.86 0.5 105.22 16.77 3.60 9.00 SPGR 0.1 86.31 13.51 3.22 7.97 0.2 89.01 13.93 3.31 8.23 0.3 92.76 14.79 3.36 8.35 0.4 97.51 15.63 3.48 8.64 0.5 101.82 16.62 3.56 8.81
下载: 导出CSV
表 3 MIXED模式下不同算法的恢复误差
Table 3. Imputation error of different algorithms in MIXED mode
算法 δ 交通流量 交通速度 RMSE/（辆·
（15 min）−1）MAPE/% RMSE/
（km·h−1）MAPE/% LRMC 0.1 89.57 13.56 2.97 7.10 0.2 92.73 14.31 2.99 7.28 0.3 95.18 14.70 3.10 7.59 0.4 99.00 15.61 3.19 7.86 0.5 103.12 16.24 2.99 8.24 PPCA 0.1 84.40 12.95 3.13 7.24 0.2 89.37 14.08 3.10 7.33 0.3 91.37 14.76 3.18 7.57 0.4 95.26 16.10 3.24 7.73 0.5 100.80 17.33 2.98 7.91 LLS 0.1 79.48 13.11 3.65 8.05 0.2 86.60 14.34 3.54 8.39 0.3 93.75 15.77 3.66 8.74 0.4 104.31 17.83 3.82 9.20 0.5 118.90 20.15 4.06 10.09 SP 0.1 82.10 12.85 2.92 6.95 0.2 85.54 13.62 3.02 7.10 0.3 88.72 14.34 3.10 7.49 0.4 94.10 15.35 3.18 7.79 0.5 99.60 16.18 3.31 8.15 SPGR 0.1 80.51 12.74 2.86 6.81 0.2 83.82 13.49 2.89 7.03 0.3 86.62 14.18 3.02 7.36 0.4 92.01 15.21 3.10 7.62 0.5 97.81 16.04 3.24 7.98
下载: 导出CSV
表 4 在交通流量数据上不同初始化方法对SPGR恢复误差的影响
Table 4. Effect of different initialization methods on SPGR imputation error on traffic flow data
辆/15 min 数据模式 KNN LLS LRMC MCAR 81.21 80.99 80.93 MAR 93.68 93.14 92.76 MIXED 87.74 87.48 86.62
下载: 导出CSV
-
[1] 陈程. 基于稀疏表示的交通数据缺失值恢复方法研究[D]. 镇江: 江苏大学, 2020. [2] HENRICKSON K, ZOU Y J, WANG Y H. Flexible and robust method for missing loop detector data imputation[J]. Transportation Research Record: Journal of the Transportation Research Board, 2015, 2527(1): 29-36. doi: 10.3141/2527-04 [3] 孙玲,刘浩,牛树云. 考虑时空相关性的固定检测缺失数据重构算法[J]. 交通运输工程学报,2010,10(5): 121-126. doi: 10.3969/j.issn.1671-1637.2010.05.021SUN Ling, LIU Hao, NIU Shuyun. Reconstructive method of missing data for location-specific detector considering spatio-temporal relationship[J]. Journal of Traffic and Transportation Engineering, 2010, 10(5): 121-126. doi: 10.3969/j.issn.1671-1637.2010.05.021 [4] CHEN Y Y, LV Y S, WANG F Y. Traffic flow imputation using parallel data and generative adversarial networks[J]. IEEE Transactions on Intelligent Transportation Systems, 2020, 21(4): 1624-1630. doi: 10.1109/TITS.2019.2910295 [5] RODRIGUES F, HENRICKSON K, PEREIRA F C. Multi-output Gaussian processes for crowdsourced traffic data imputation[J]. IEEE Transactions on Intelligent Transportation Systems, 2019, 20(2): 594-603. doi: 10.1109/TITS.2018.2817879 [6] 李林超,曲栩,张健,等. 基于特征级融合的高速公路异质交通流数据修复方法[J]. 东南大学学报(自然科学版),2018,48(5): 972-978. doi: 10.3969/j.issn.1001-0505.2018.05.029LI Linchao, QU Xu, ZHANG Jian, et al. Missing value imputation method for heterogeneous traffic flow data based on feature fusion[J]. Journal of Southeast University (Natural Science Edition), 2018, 48(5): 972-978. doi: 10.3969/j.issn.1001-0505.2018.05.029 [7] LI Q, TAN H C, WU Y K, et al. Traffic flow prediction with missing data imputed by tensor completion methods[J]. IEEE Access, 2020, 8: 63188-63201. doi: 10.1109/ACCESS.2020.2984588 [8] CHEN X Y, HE Z C, CHEN Y X, et al. Missing traffic data imputation and pattern discovery with a Bayesian augmented tensor factorization model[J]. Transportation Research Part C: Emerging Technologies, 2019, 104: 66-77. doi: 10.1016/j.trc.2019.03.003 [9] CHEN X B, CAI Y F, YE Q L, et al. Graph regularized local self-representation for missing value imputation with applications to on-road traffic sensor data[J]. Neurocomputing, 2018, 303: 47-59. doi: 10.1016/j.neucom.2018.04.029 [10] BOYD S, PARIKH N, CHU E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends, 2011, 3(1): 1-122. [11] NIE F P, WANG H, HUANG H, et al. Joint Schatten p-norm and ℓ p-norm robust matrix completion for missing value recovery[J]. Knowledge and Information Systems, 2015, 42(3): 525-544. doi: 10.1007/s10115-013-0713-z [12] BARTELS R H, STEWART G W. Solution of the matrix equation AX + XB = C[J]. Communications of the ACM, 1972, 15(9): 820-826. doi: 10.1145/361573.361582 [13] LI L C, ZHANG J, WANG Y G, et al. Missing value imputation for traffic-related time series data based on a multi-view learning method[J]. IEEE Transactions on Intelligent Transportation Systems, 2019, 20(8): 2933-2943. doi: 10.1109/TITS.2018.2869768 [14] ZHANG Y, LIU Y C. Data imputation using least squares support vector machines in urban arterial streets[J]. IEEE Signal Processing Letters, 2009, 16(5): 414-417. doi: 10.1109/LSP.2009.2016451 期刊类型引用(4)
1. 安少聪. 热耦合映射下的光伏电池组升温异常态势识别. 自动化与仪器仪表. 2024(07): 284-287+292 . 
百度学术2. 汤伟，漆苏应，杨晓东，李国强. 基于随机森林和最近邻插值法的交通流量数据修复方法. 科学技术与工程. 2024(32): 14056-14065 . 百度学术3. 武江南，张红梅，赵永梅，曾航，胡钢. 基于张量加权与截断核范数的交通数据修复方法. 计算机科学. 2023(08): 45-51 . 百度学术4. 武江南，张红梅，赵永梅，曾航. 基于张量奇异值理论的交通数据重构方法. 计算机应用研究. 2022(05): 1449-1453+1459 . 百度学术其他类型引用(2)
-
下载:
下载:
百度学术
点击查看大图
计量
- 文章访问数: 502
- HTML全文浏览量: 214
- PDF下载量: 38
- 被引次数: 6