陳小波，梁書榮，柯 佳，陳 玲，胡 煜

（1.江蘇大學(xué)汽車工程研究院，江蘇 鎮(zhèn)江 212013；2.江蘇大學(xué)管理學(xué)院，江蘇 鎮(zhèn)江 212013）

準(zhǔn)確、完整的交通數(shù)據(jù)是實現(xiàn)交通大數(shù)據(jù)分析與挖掘的基礎(chǔ).然而，實際交通環(huán)境中，由于檢測器故障、惡劣天氣、通信中斷等原因，采集的交通數(shù)據(jù)通常存在缺失值.例如，北京交通流數(shù)據(jù)檢測設(shè)備采集的數(shù)據(jù)中每天約有10%的數(shù)據(jù)缺失[1].傳統(tǒng)的機器學(xué)習(xí)算法，如支持向量機、神經(jīng)網(wǎng)絡(luò)等無法處理不完整數(shù)據(jù)，嚴(yán)重地影響交通預(yù)測與分析智能算法的性能.因此，數(shù)據(jù)缺失問題已成為交通大數(shù)據(jù)分析的關(guān)鍵挑戰(zhàn)之一.

近年來，如何對交通數(shù)據(jù)中的缺失值進(jìn)行準(zhǔn)確恢復(fù)受到國內(nèi)外學(xué)者的廣泛關(guān)注.由于路網(wǎng)的拓?fù)浣Y(jié)構(gòu)、人們出行行為的規(guī)律性和周期性等，同一路網(wǎng)中不同位置的檢測器在不同檢測時間所記錄的數(shù)據(jù)具有較強的時空相關(guān)性，因此，充分挖掘交通數(shù)據(jù)內(nèi)在的時空相關(guān)性是實現(xiàn)缺失數(shù)據(jù)準(zhǔn)確恢復(fù)的關(guān)鍵.目前，交通數(shù)據(jù)缺失值恢復(fù)的主要方法可分為三類：基于預(yù)測的方法、基于插值的方法和基于統(tǒng)計學(xué)習(xí)的方法.基于預(yù)測的方法是利用歷史數(shù)據(jù)的特性將缺失值作為預(yù)測的目標(biāo)進(jìn)行恢復(fù)；例如，Henrickson等[2]針對交通數(shù)據(jù)中的連續(xù)缺失現(xiàn)象，提出一種預(yù)測均值匹配多重恢復(fù)方法，即使是全天數(shù)據(jù)缺失，也能取得良好的效果；但這類方法基于歷史數(shù)據(jù)建模，忽略了缺失值之后的數(shù)據(jù)所提供的信息，導(dǎo)致對信息利用不充分，同時，當(dāng)數(shù)據(jù)變化波動較大時，恢復(fù)數(shù)據(jù)的準(zhǔn)確性將顯著降低.基于插值的方法是利用缺失值前后鄰近的數(shù)據(jù)或者缺失值所在的鄰近檢測器數(shù)據(jù)取加權(quán)平均作為缺失值的恢復(fù)值；孫玲等[3]通過研究缺失值與觀測值的相關(guān)性，對不同相關(guān)性的數(shù)據(jù)給予不同的權(quán)重，通過加權(quán)平均作為插補值；基于插值的方法計算簡單，但要求缺失數(shù)據(jù)前后的相似性較高，并且當(dāng)缺失值周圍數(shù)據(jù)也存在缺失時，這類方法將會出現(xiàn)較大誤差.基于統(tǒng)計學(xué)習(xí)的方法通過引入機器學(xué)習(xí)理論與算法[4-5]建模交通數(shù)據(jù)的時空關(guān)系，實現(xiàn)缺失值恢復(fù)；李林超等[6]提出一種基于多源數(shù)據(jù)融合的異質(zhì)交通流數(shù)據(jù)修復(fù)方法，通過自編碼器提取高維多源交通流數(shù)據(jù)的時空特征，并采用隨機森林估計數(shù)據(jù)中的缺失值；這類方法主要是建立符合數(shù)據(jù)分布的模型，并不斷調(diào)整模型參數(shù)，使之實現(xiàn)對數(shù)據(jù)的最優(yōu)擬合，進(jìn)而實現(xiàn)缺失值的恢復(fù).

低秩矩陣補全（low-rank matrix completion，LRMC）、低秩張量補全（low-rank tensor completion，LRTC）算法被成功應(yīng)用于交通數(shù)據(jù)缺失值恢復(fù)[7].由于交通數(shù)據(jù)的時空相關(guān)性，交通數(shù)據(jù)矩陣具有低秩結(jié)構(gòu)，通過矩陣秩最小化實現(xiàn)對缺失數(shù)據(jù)矩陣的恢復(fù).Chen等[8]將交通量數(shù)據(jù)表示為張量形式，并利用貝葉斯張量分解理論交替優(yōu)化缺失值和概率模型參數(shù).Li等[7]對4種LRTC和4種預(yù)測模型進(jìn)行了研究，表明提高缺失值恢復(fù)精度可以改善交通預(yù)測模型的性能.基于LRMC的缺失值恢復(fù)方法取得了較好的效果，然而，最小化秩函數(shù)是一個NP （nondeterministic polynomial）難的非凸優(yōu)化問題，難以在多項式時間內(nèi)求解.一些方法將秩函數(shù)凸松弛為矩陣核范數(shù)進(jìn)行求解，最小化核范數(shù)是一個凸優(yōu)化問題，具有全局最優(yōu)解.然而，凸松弛問題的最優(yōu)解可能會偏離原問題的最優(yōu)解.同時，大多數(shù)基于LRMC的恢復(fù)方法只利用了數(shù)據(jù)矩陣的全局時空相關(guān)性，而沒有充分利用數(shù)據(jù)的局部結(jié)構(gòu)特性.

基于上述分析，本文提出一種基于圖正則化和Schatten-p范數(shù)（Schatten-pnorm and graph regularization，SPGR）的交通數(shù)據(jù)缺失值恢復(fù)方法.首先，應(yīng)用已有的數(shù)據(jù)缺失值恢復(fù)方法對缺失值進(jìn)行初步估計，基于此，建立刻畫交通數(shù)據(jù)局部近鄰結(jié)構(gòu)的圖模型以反映不同樣本之間的相似度；依據(jù)交通數(shù)據(jù)的低秩結(jié)構(gòu)這一全局先驗信息，應(yīng)用Schatten-p范數(shù)逼近矩陣的秩函數(shù)，提出融合圖正則化和Schattenp范數(shù)最小化的缺失值恢復(fù)模型.然后，基于交替方向乘子法（alternating direction method of multipliers，ADMM）框架設(shè)計優(yōu)化算法，實現(xiàn)對模型的有效求解.最后，通過在真實路網(wǎng)交通流量與速度數(shù)據(jù)上的實驗對比分析驗證該方法的有效性.

1 問題描述

將交通數(shù)據(jù)（如交通流量、速度）表示為變量矩陣X=(X1,X2,···,Xi,···,XN)∈RM×N.N為樣本總數(shù)；M為檢測器一天采集的交通數(shù)據(jù)個數(shù)，即一個樣本中的數(shù)據(jù)量；Xi=(xi(1),xi(2),···,xi(q),···,xi(Q))T，xi(q)為第i個樣本中第q個時刻的交通數(shù)據(jù).設(shè) ?表示X中已知元素的下標(biāo)集合，即X?為觀測值集合，表示X中缺失元素的下標(biāo)集合，因此，交通數(shù)據(jù)恢復(fù)問題就是根據(jù)觀測值集合X?準(zhǔn)確估計缺失值集合.

2 基于SPGR的交通數(shù)據(jù)缺失值恢復(fù)算法

2.1 算法框架

本文提出的基于圖正則化和Schatten-p范數(shù)最小化的交通數(shù)據(jù)缺失值恢復(fù)算法，算法流程如圖1所示，具體包括以下步驟：

圖1 交通數(shù)據(jù)恢復(fù)算法框架Fig.1 Diagram of traffic data imputation algorithm

步驟1應(yīng)用已有的數(shù)據(jù)缺失值恢復(fù)方法（如LRMC）得到缺失值的初始估計，其中，當(dāng)xi(q)已知時，(q)=xi(q).

步驟2根據(jù)初始估計的結(jié)果，計算任意兩個樣本之間的距離，表示為矩陣D=(dij)∈RN×N，其中，dij為第i個樣本與第j個樣本間的距離，j=1,2,···,N，通過加權(quán)歐幾里得公式[9]計算，如式（1）所示.

式中：θq為兩樣本在第q個時刻的差值權(quán)重；θ?q為 θq的歸一化值；α為相對小的正數(shù)，本文取0.1.

步驟3對樣本Xi，根據(jù)距離矩陣D選出與之最接近的K個樣本，表示為集合Ne(i).然后構(gòu)造鄰接矩陣（權(quán)矩陣）S=(sij)∈RN×N.

步驟4基于X與S，建立SPGR模型并求解，得到最終恢復(fù)值Y=(Y1,Y2,···,YN).

2.2 SPGR模型

如前所述，交通數(shù)據(jù)具有較強的時空相關(guān)性，因此，交通數(shù)據(jù)矩陣具有低秩結(jié)構(gòu)[1]，可將缺失值恢復(fù)問題轉(zhuǎn)化為矩陣秩最小化問題，如式（5）所示.

式中：A∈RM×N為觀測值矩陣；P?(?)為 ? 的投影運算，如式（6）所示.

然而，由于秩函數(shù)是非凸且不連續(xù)，所以式（5）是一個NP難問題，難以求解.為解決這一問題，用矩陣的Schatten-p（0

式中：∥X∥sp為X的Schatten-p范數(shù)；σt為X的第t個奇異值.

當(dāng)p= 1 時，Schatten-1范數(shù)即為核范數(shù).當(dāng)p越趨近于0時，X的Schatten-p范數(shù)越接近秩函數(shù).當(dāng)p= 0 時，式（7）就是X的秩函數(shù).總的來說，與核范數(shù)相比，Schatten-p范數(shù)對秩函數(shù)的逼近能力更強，可以更精確刻畫數(shù)據(jù)的低秩結(jié)構(gòu).因此，本文構(gòu)建基于Schatten-p范數(shù)最小化的缺失值恢復(fù)模型，如式（8）所示.

式中：∥?∥F為矩陣的F范數(shù).

同時，鄰近的數(shù)據(jù)樣本具有相似的特征，為更好利用這種信息，在進(jìn)行恢復(fù)時，將所有鄰域樣本之間的距離限制在適當(dāng)?shù)姆秶鷥?nèi)，以防止被恢復(fù)樣本與鄰域樣本的差異過大.基于上述分析，提出圖正則化來刻畫這種局部鄰近結(jié)構(gòu)：

式中：L=E?S，為圖的拉普拉斯矩陣；E為對角元素是的對角矩陣.

結(jié)合式（8）、（9），SPGR的目標(biāo)函數(shù)為

式中：λ>0為控制圖正則化項的權(quán)重常數(shù).

2.3 基于ADMM的優(yōu)化算法

本文采用交替方向乘子法（ADMM）框架尋求式（10）的最優(yōu)解[10].首先，引入輔助變量W和Z，將式（10）轉(zhuǎn)化為式（11）所示等價問題.

構(gòu)造式（11）的增廣拉格朗日函數(shù)為

式中：U、V為拉格朗日乘子；μ1>0，μ2>0，為懲罰系數(shù).

增廣拉格朗日函數(shù)融合了罰函數(shù)法與拉格朗日乘子法的優(yōu)點，依據(jù)ADMM框架對優(yōu)化變量進(jìn)行迭代求解.一個變量進(jìn)行優(yōu)化時，固定其他變量，通過變量交替近似求解，實現(xiàn)算法結(jié)果的優(yōu)化.

1）固定變量Z和X，

式中：l為ADMM算法中的迭代次數(shù)；Wl、Zl、Xl、Ul、Vl分別為W、Z、X、U、V迭代l次后對應(yīng)的數(shù)值.

式（13）的最優(yōu)解通過迭代算法[11]求出.

2）固定變量W和X，

式（14）中對Z求導(dǎo)，并將導(dǎo)數(shù)設(shè)為0，得到

3）固定變量W和Z，

式（16）中對X求導(dǎo)，并將導(dǎo)數(shù)設(shè)為0，得到

式（17）為Sylvester方程，對于方程AX+XB=C，其解表示為X=s(A,B,C)

[12]，s(?)為Sylvester方程的求解運算.因此，

4）更新乘子Ul+ 1和Vl+ 1，

變量W、Z和X按照上述規(guī)則迭代更新，直到算法收斂，收斂條件 ε為

式中：r為很小的正數(shù)，本文取 1×10?5.

當(dāng) ε

綜上，求解式（11）的優(yōu)化算法流程如圖2所示.

圖2 ADMM流程Fig.2 Flow chart of ADMM

3 實驗驗證

3.1 交通數(shù)據(jù)與實驗方案

為了評估SPGR方法的數(shù)據(jù)恢復(fù)性能，在真實的交通流量和交通速度數(shù)據(jù)上進(jìn)行實驗分析.交通流量數(shù)據(jù)來自美國俄勒岡州波特蘭市交通信息中心，從I205和I84州際公路構(gòu)成的路網(wǎng)中選擇40個只有極少缺失數(shù)據(jù)的檢測站采集數(shù)據(jù).由于工作日的交通流量數(shù)據(jù)與周末和節(jié)假日的數(shù)據(jù)存在較大差異，因此，選取2015年中連續(xù)30個工作日的交通流數(shù)據(jù)進(jìn)行研究.傳感器采樣間隔為15 min，最終獲取30 × 40 = 1200個樣本，構(gòu)造為96 × 1200的數(shù)據(jù)矩陣.此外，交通速度數(shù)據(jù)集為中國廣州兩個月（2016年8月1日至2016年9月30日共61 d）內(nèi)以10 min為間隔的7個路段（主要包括城市快速路和干線）的速度信息.因此，得到7 × 61 = 421個樣本，構(gòu)造為144 × 421數(shù)據(jù)矩陣.圖3展示了8個傳感器在同一天的交通數(shù)據(jù)變化情況.

圖3 同一天中不同傳感器交通流和交通速度的變化情況Fig.3 Changes in traffic flow and speed from different sensors over the same day

反映交通數(shù)據(jù)缺失值的復(fù)雜分布，模擬3種常見的數(shù)據(jù)缺失模式：1）完全隨機缺失（missing completely at random，MCAR），缺失值獨立于其他缺失數(shù)據(jù)或已知數(shù)據(jù)，表現(xiàn)為一組隨機分布的孤立點；2）隨機缺失（missing at random，MAR），交通數(shù)據(jù)表現(xiàn)為連續(xù)缺失的現(xiàn)象，即缺失值的恢復(fù)依賴于相鄰的缺失值；3）混合缺失（mixture of MCAR and MAR，MIXED），MAR與MCAR的混合比例各為0.5.圖4為不同缺失模式的示例，圖中每行表示一個流量樣本，每列表示一個變量，黑色表示缺失值.

圖4 數(shù)據(jù)缺失模式模擬示例Fig.4 Simulation examples of data missing modes

為綜合比較不同數(shù)據(jù)恢復(fù)方法的有效性，將提出的SPGR模型、去除圖正則化的SP模型（Schattenp范數(shù)最小化，λ=0）與3種缺失值恢復(fù)方法：LRMC、概率主成分分析（probabilistic principal component analysis，PPCA）、局部最小二乘（local least squares，LLS）進(jìn)行比較[1,13-14].這3種對比方法涵蓋了數(shù)據(jù)缺失值恢復(fù)的主流技術(shù)，包括低秩矩陣補全、概率模型和回歸模型.

實驗中，按照缺失模式和缺失比例模擬缺失值.其中，缺失率 δ 定義為缺失數(shù)據(jù)數(shù)量與總數(shù)據(jù)量之比，以0.1為步長將 δ 從0.1增加到0.5，研究不同缺失率對恢復(fù)性能的影響.為衡量不同算法的恢復(fù)性能，采用缺失項的恢復(fù)值與真實值之間的均方根誤差（RMSE，eRMSE）和平均絕對百分比誤差（MAPE，eMAPE）表示，分別為

式中：C為缺失值的總數(shù)目；和分別為真實值和恢復(fù)值.

RMSE和MAPE越小，算法的恢復(fù)性能越好.為準(zhǔn)確評估5種缺失值恢復(fù)方法在兩種交通數(shù)據(jù)上的性能，減少隨機性對實驗結(jié)果的影響，每種實驗均重復(fù)5次，取實驗結(jié)果的誤差平均值作為評價缺失值恢復(fù)方法性能的依據(jù).

3.2 實驗結(jié)果分析

表1～ 3分別列出了不同算法在MCAR、MAR和MIXED模式下的恢復(fù)誤差，根據(jù)實驗結(jié)果，可以得到以下結(jié)論：1）MCAR缺失模式下，各種算法的恢復(fù)誤差最??；而在MAR缺失模式下，因缺失大量相關(guān)信息，導(dǎo)致數(shù)據(jù)恢復(fù)誤差較大；此外，每種方法的恢復(fù)誤差隨著缺失率的增加而增加.2）在恢復(fù)性能方面，LRMC和PPCA方法不能很好地處理內(nèi)部結(jié)構(gòu)復(fù)雜的數(shù)據(jù)集，導(dǎo)致總體性能比其他方法差；當(dāng)缺失率較低時，PPCA的性能較好，而當(dāng)缺失率增加時，LRMC的性能優(yōu)于PPCA；LLS在低缺失率下恢復(fù)性能較好，然而，當(dāng)缺失率增加時，其恢復(fù)性能會迅速退化.3）本文提出的SPGR算法獲得了更好的恢復(fù)性能；與LLS、PPCA、LRMC方法相比，在MCAR缺失模式下，RMSE降低了3.02% ～ 22.31%，在MAR缺失模式下，RMSE降低了3.23% ～ 28.49%，在MIXED缺失模式下，RMSE降低了3.05% ～ 21.56%；特別是當(dāng)缺失率大于30%時，恢復(fù)誤差降低率越大，相對于LLS算法誤差降低率高達(dá)28.49%，表明該方法可以有效挖掘觀測數(shù)據(jù)的內(nèi)在關(guān)聯(lián)，實現(xiàn)準(zhǔn)確的缺失值恢復(fù).

表1 MCAR模式下不同算法的恢復(fù)誤差Tab.1 Imputation error of different algorithms in MCAR mode

表2 MAR模式下不同算法的恢復(fù)誤差Tab.2 Imputation error of different algorithms in MAR mode

3.3 算法參數(shù)的影響

本文提出的SPGR模型涉及3個參數(shù)：Schattenp范數(shù)的p值、K近鄰（K-nearest neighbor, KNN）方法中的K值、圖正則化的權(quán)重常數(shù) λ.p值決定了對矩陣秩函數(shù)的近似程度，K值決定了用于重建每個樣本的近鄰數(shù)量，λ控制基于局部鄰近的圖正則化的影響.為得到模型參數(shù)的最優(yōu)值，調(diào)整其中一個參數(shù)，固定另外兩個參數(shù)，每次參數(shù)改變時記錄實驗結(jié)果.以交通流量數(shù)據(jù)為例，圖5給出了不同參數(shù)下，缺失值恢復(fù)誤差RMSE的變化.由圖可知：不同缺失率下，p具有不同的最優(yōu)值，由于Schatten-p范數(shù)比核范數(shù)（p=1）更能逼近秩函數(shù)，從而獲得更好的恢復(fù)結(jié)果；如果K值過大或過小，都會導(dǎo)致恢復(fù)精度較差.這是因為K值太小，選擇代表目標(biāo)樣本的相鄰樣本過少，導(dǎo)致用于缺失值恢復(fù)的可用信息不充分，K值過大，遠(yuǎn)離目標(biāo)的樣本將參與重建，也將降低恢復(fù)精度.對于 λ的影響，也可以得出與K值類似的結(jié)論.

表3 MIXED模式下不同算法的恢復(fù)誤差Tab.3 Imputation error of different algorithms in MIXED mode

圖5 SPGR模型在交通流量數(shù)據(jù)上的RMSE隨參數(shù) p、K、λ的變化Fig.5 RMSEs of SPGR model on traffic flow data varied with the parameters ofp、K、λ

3.4 初始化影響

根據(jù)SPGR算法的第一步，為建立表征樣本間相鄰關(guān)系的圖矩陣，需要選擇一種已有的恢復(fù)方法對交通數(shù)據(jù)的缺失值進(jìn)行初始估計.本節(jié)將采用3種不同的初始化方法，即KNN、LLS、LRMC，進(jìn)行缺失值的初始估計，進(jìn)而研究其對SPGR算法性能的影響，實現(xiàn)敏感性分析.以交通流數(shù)據(jù)為例，δ=0.3下的實驗結(jié)果如表4所示.可以看出，不同的初始化方法對SPGR的恢復(fù)性能影響較小，在其他缺失率下也可以觀察到類似現(xiàn)象.這驗證了SPGR對初始值具有較好的魯棒性.

表4 在交通流量數(shù)據(jù)上不同初始化方法對SPGR恢復(fù)誤差的影響Tab.4 Effect of different initialization methods on SPGR imputation error on traffic flow data 輛/15 min

4 結(jié) 論

提出一種融合圖正則化與Schatten-p范數(shù)最小化的交通數(shù)據(jù)缺失值恢復(fù)方法.該方法采用Schatten-p范數(shù)逼近矩陣的秩函數(shù)，對數(shù)據(jù)的低秩先驗信息進(jìn)行約束.通過實驗分析，得到以下結(jié)論：

1）將圖正則化融入到數(shù)據(jù)恢復(fù)框架中，有利于更好地利用數(shù)據(jù)的局部鄰近結(jié)構(gòu).

2）基于真實的高速公路交通量和速度數(shù)據(jù)進(jìn)行仿真實驗表明，提出的方法相對于其他多種方法恢復(fù)誤差降低了3.02%以上，特別是缺失率大于0.3時，誤差降低率達(dá)到20%以上.

3）在未來的工作中，將進(jìn)一步研究交通數(shù)據(jù)在時序上的規(guī)律性，以提升缺失值恢復(fù)的精度.