張貝娜 ，馮震華 ，張豐 ，杜震洪 *，劉仁義 ，周芹

（1.浙江大學(xué)浙江省資源與環(huán)境信息系統(tǒng)重點實驗室，浙江杭州310028；2.浙江大學(xué)地理信息科學(xué)研究所浙江杭州310027；3.北京超圖軟件服務(wù)有限公司，北京100015）

由于城市化的飛速發(fā)展，伴隨而來的城市空氣質(zhì)量問題已成為全球居民普遍關(guān)心的問題。城市物理空間中部署了大量的傳感器，用以持續(xù)和連續(xù)地監(jiān)測環(huán)境信息，生成了海量的帶有地理位置信息的時間序列數(shù)據(jù)，幫助人類更好地理解周邊環(huán)境[1]。在實際監(jiān)測中，由于傳感器故障、數(shù)據(jù)傳輸錯誤等，監(jiān)測數(shù)據(jù)常會出現(xiàn)隨機性甚至連續(xù)性缺失，嚴(yán)重影響實時監(jiān)測數(shù)據(jù)的查看和后續(xù)數(shù)據(jù)挖掘訓(xùn)練集的組織。

經(jīng)典數(shù)據(jù)補全方法包括數(shù)學(xué)統(tǒng)計、時間序列預(yù)測、空間插值等。RUBIN[2]提出用多重補全方法迭代進行數(shù)據(jù)補全；LI等[3]比較了周期性填補法、均值填補法和三次樣條函數(shù)插值法3種方法對時間序列缺失數(shù)據(jù)的補全情況；TURRADO等[4]采用自適應(yīng)的反距離加權(quán)方法對太陽輻射缺失數(shù)據(jù)進行補全。

上述方法僅考慮單視圖，即僅使用一種維度信息來表示一個對象[5]。目前，也有考慮多視圖信息的研究，例如CHEN等[6]采用時空k最近鄰方法對推文標(biāo)記數(shù)據(jù)進行補全，YI等[7]提出結(jié)合經(jīng)驗統(tǒng)計模型和數(shù)據(jù)驅(qū)動算法的線性回歸方法進行時空數(shù)據(jù)補全。因當(dāng)前該類方法普遍立足于線性關(guān)系，不利于多視圖信息的關(guān)聯(lián)挖掘。

隨著機器學(xué)習(xí)的發(fā)展，出現(xiàn)了一些基于非線性關(guān)系的數(shù)據(jù)補全方法，例如XIA等[8]提出采用非負矩陣分解補全社交圖像的缺失標(biāo)記，魏晶茹等[9]采用支持向量機和粒子群優(yōu)化方法補全污水處理廠的污染物測量數(shù)據(jù)。因當(dāng)前該類方法普遍基于單視圖信息，不利于連續(xù)缺失數(shù)據(jù)的補全。

在數(shù)據(jù)預(yù)測問題的解決中，由于未來條件的不可預(yù)知性，其他影響因素（例如氣象條件）往往作為輸入補充條件。而數(shù)據(jù)補全問題是基于歷史和未來數(shù)據(jù)的，本身蘊含較為豐富的信息，因此，本文重點研究監(jiān)測數(shù)據(jù)的缺失數(shù)據(jù)重建。

針對當(dāng)前數(shù)據(jù)補全方法的優(yōu)點和不足，以北京市一年的城市空氣質(zhì)量數(shù)據(jù)為研究對象，基于指數(shù)移動平均、普通克里金、非凸矩陣完備，提出了一種基于時空多視圖BP神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)補全方法，用于補全隨機缺失、時間連續(xù)缺失和空間連續(xù)缺失的數(shù)據(jù)，并用實驗驗證和分析該方法的有效性和優(yōu)越性。

1 時空多視圖

傳統(tǒng)數(shù)據(jù)采用單視圖組織形式。隨著大數(shù)據(jù)時代數(shù)據(jù)獲取和存儲能力的顯著提高，同一對象的數(shù)據(jù)表達形式日益豐富，例如一幅圖像可被形狀、紋理和顏色共同表示[10]。多角度的數(shù)據(jù)組織形式能夠描述同一對象的不同方面，能夠刻畫同一目標(biāo)的不同特性，對信息的描述具有補充價值，被稱為多視圖數(shù)據(jù)[11]。

理論上，多視圖只是一個泛指，它可以是不同特征，也可以是不同的模型架構(gòu)和不同的數(shù)據(jù)集。本文重點關(guān)注多視圖特征，從局部時間信息、局部空間信息和全局時空信息3個角度入手，構(gòu)建3時空多視圖特征，提取不同角度的數(shù)據(jù)內(nèi)容，挖掘數(shù)據(jù)本身的潛在信息。

1.1 基于指數(shù)移動平均的時間視圖

采用指數(shù)移動平均（exponential moving average，EMA）模型表征時間視圖。該模型最早由BROWN[12]提出，對過去的時間序列數(shù)據(jù)進行加權(quán)平均，時間間隔越近的數(shù)據(jù)權(quán)重越大，可用遞推公式表示：

其中，St是t時刻的平滑值，vt是t時刻的實際值，α是平滑常數(shù)，滿足0＜α＜1。

數(shù)據(jù)補全時應(yīng)同時考慮包含過去和未來的雙向數(shù)據(jù)，為充分利用數(shù)據(jù)，采用拓展為雙向的指數(shù)移動平均方法。此外，考慮到過于久遠的數(shù)據(jù)對當(dāng)前數(shù)據(jù)影響較小，故只將最近一定時間范圍內(nèi)的數(shù)據(jù)作為輸入進行計算，公式為

其中，是t時刻的平滑值，vt是t時刻的實際觀察值，threshold是時 間 窗 口，將 [t－threshold,t+threshold]范圍內(nèi)除t外的數(shù)據(jù)作為t時刻數(shù)據(jù)的計算輸入，α是平滑常數(shù)，滿足0＜α＜1。

1.2 基于普通克里金的空間視圖

采用普通克里金（ordinary Kriging，OK）模型表征空間視圖。該模型的雛形來自南非礦業(yè)工程師KRIGE[13]，法國地質(zhì)數(shù)學(xué)家MATHERON[14]對其進行了整理歸納。方法的依據(jù)是地理學(xué)第一定律[15]，即已知空間上的地理屬性具有空間相關(guān)性，且越相近的事物越相似，其基本公式為

克里金插值采用統(tǒng)計規(guī)律信息計^算權(quán)重系數(shù)，目標(biāo)是得到滿足點(xs,ys)處的估計值與真實值vs的差最小的一套最優(yōu)系數(shù)，即滿足最優(yōu)系數(shù)條件：

同時滿足無偏估計條件：

克里金具有不同變種，主要區(qū)別在于假設(shè)條件的不同。本文采用的普通克里金插值方法的假設(shè)條件為：空間屬性v是均一的；對于空間任意一點(x,y)，都有同樣的期望與方差。

基于上述3項條件，經(jīng)過公式推導(dǎo)，可得到權(quán)重系數(shù)λi的求解等價式：

其中，γij是半方差函數(shù)值，φ是拉格朗日乘數(shù)。

普通克里金的求解過程最終轉(zhuǎn)換為半方差函數(shù)的求解問題，其定義為

已知點間的半方差函數(shù)值γij可通過實測值求得；而已知點和未知點間的半方差函數(shù)值γis通過統(tǒng)計信息和空間相似性的關(guān)系求得。

克里金插值假設(shè)觀測點i和j的半方差函數(shù)γij與兩點間的距離dij存在函數(shù)關(guān)系。首先，繪制兩兩已知點間“距離-半方差函數(shù)值”的散點圖，然后尋找一條最優(yōu)的擬合曲線（常用半方差模型，有球面、圓、指數(shù)、對數(shù)等模型）擬合d和γ的關(guān)系，得到半方差函數(shù)。對于任意兩點(xi,yi)和(xj,yj)，先計算其距離dij，然后根據(jù)函數(shù)關(guān)系求得這兩點間的半方差γij，最終求得權(quán)重系數(shù)和待插值點的估計值。

1.3 基于非凸矩陣完備的時空關(guān)聯(lián)視圖

采用非凸矩陣完備（non-convex matrix completion）模型表征時空關(guān)聯(lián)視圖。矩陣完備是在一個矩陣有元素缺失的情況下，通過對缺失值之外的有效元素位置進行采樣，進行一定的計算和分析，最后補全整個矩陣的一種算法[16]。從數(shù)學(xué)角度而言，絕大多數(shù)矩陣補全問題是通過尋找低秩條件下原始矩陣的近似矩陣解決的，矩陣完備問題可形式化表示為

其中，X是目標(biāo)完整矩陣，Ω是對應(yīng)矩陣中所有能觀測到的坐標(biāo)位置集合，MΩ是選取X在Ω位置的掩膜操作，Y是實測樣本。求解式(8)是NP難問題，可將其轉(zhuǎn)換為具有高概率恢復(fù)矩陣能力的凸松弛問題：

其中，‖‖X*是核范數(shù)，被定義為奇異值之和。

其中，‖X‖p*是矩陣奇異值的lp范數(shù)。

稍微放寬式(10)中的限制條件，目標(biāo)問題便轉(zhuǎn)變?yōu)?/p>

其中，‖·‖F(xiàn)是弗羅貝尼烏斯范數(shù)（Frobenius norm），ε為誤差值。

將非凸矩陣完備進一步轉(zhuǎn)換，可數(shù)值計算量化的目標(biāo)問題：

最終可通過轉(zhuǎn)換優(yōu)化和冷卻算法進行迭代，得到補全的目標(biāo)矩陣。

2 時空多視圖BP神經(jīng)網(wǎng)絡(luò)模型

時空多視圖BP神經(jīng)網(wǎng)絡(luò)數(shù)據(jù)補全模型（spatiotemporal multi-view-BP,STMV-BP）框架見圖1，包括時空多視圖和BP神經(jīng)網(wǎng)絡(luò)2個重點內(nèi)容。時空多視圖從時間、空間和時空3個互為補充的視角切入，綜合抽取不同角度的時空數(shù)據(jù)信息；BP神經(jīng)網(wǎng)絡(luò)融合時空多視圖特征，映射時空多視圖和目標(biāo)補全數(shù)據(jù)之間的非線性關(guān)系，構(gòu)成面向城市空氣質(zhì)量數(shù)據(jù)補全的方法架構(gòu)。

圖1 時空多視圖BP神經(jīng)網(wǎng)絡(luò)模型框架Fig.1 Framework of spatio-temporal multi-view-BP model

本文設(shè)計了3個獨立的時間、空間和時空視圖，將其作為特征，沒將所有數(shù)據(jù)作為神經(jīng)網(wǎng)絡(luò)輸入的原因如下：

（1）從特征空間的角度看，時間和空間視圖使用的基礎(chǔ)數(shù)據(jù)無重疊，提供了不同視角的信息；

（2）從模型的角度看，3個視圖分別模擬局部和全局信息，表征局部信息的時間視圖和空間視圖從本質(zhì)上說是回歸問題，表征全局信息的時空關(guān)聯(lián)視圖則是非線性插值問題，需采用不同的模型處理；

（3）從網(wǎng)絡(luò)訓(xùn)練的角度看，將所有數(shù)據(jù)作為特征模型參數(shù)易導(dǎo)致效率降低，基于有限訓(xùn)練數(shù)據(jù)的時空特征提取，有利于提升參數(shù)訓(xùn)練效率，提高模型準(zhǔn)確度。

本文采用BP神經(jīng)網(wǎng)絡(luò)而非線性回歸或其他復(fù)雜神經(jīng)網(wǎng)絡(luò)進行特征融合，主要出于以下兩點考慮：

（1）BP神經(jīng)網(wǎng)絡(luò)能夠映射輸入數(shù)據(jù)和目標(biāo)數(shù)據(jù)間的線性關(guān)系，彌補了線性回歸僅能進行線性關(guān)系映射的不足，能夠更深入地挖掘時空多視圖信息；

（2）BP神經(jīng)網(wǎng)絡(luò)較其他復(fù)雜神經(jīng)網(wǎng)絡(luò)建模簡單，更重要的是，對數(shù)據(jù)補全方法的時空多視圖特征，最終凝練出的3個視圖已能有效映射此類簡單維度的關(guān)系。

3 實驗與分析

硬件環(huán)境：處理器Intel（R）Core（TM）i7-4770 CPU@3.40 Hz，內(nèi)存8.00 GB。

軟件環(huán)境：操作系統(tǒng)Microsoft Windows7（64bit），Matlab版本 2017b。

數(shù)據(jù)來源于“城市計算”（urban computing）項目公開數(shù)據(jù)集[1,19]，采用2014年5月1日至2015年4月30日北京市36個站點的空氣污染物小時濃度數(shù)據(jù)，站點分布見圖2，污染物包括 PM2.5、PM10、NO2、CO、O3和SO26種。

圖2 北京市空氣質(zhì)量監(jiān)測站點分布Fig.2 Distribution of air quality monitoring stations in Beijing

對278 023條原始空氣質(zhì)量記錄數(shù)據(jù)進行預(yù)處理。首先，分站點規(guī)整化數(shù)據(jù)，將各站點數(shù)據(jù)處理為按照時間序列排列的固定1 h間隔的規(guī)整數(shù)據(jù)，并將缺失數(shù)據(jù)用NaN值進行填充。然后，根據(jù)污染物類別劃分?jǐn)?shù)據(jù)，以順序時刻為行，以站點為列，將數(shù)據(jù)劃分為對應(yīng)污染物的6個數(shù)據(jù)矩陣，每個矩陣大小為8 760×36。其中，矩陣行大小代表365 d×24 h·d-1共8 760 h，矩陣列大小代表36個站點。

各空氣污染物數(shù)據(jù)存在隨機缺失、時間連續(xù)缺失和空間連續(xù)缺失3種缺失模式，3種缺失模式可能存在交叉，缺失統(tǒng)計情況如表1所示。

實驗采用均方根誤差（root mean squared error，RMSE）、平均絕對誤差（mean absolute error，MAE）和平均相對誤差（mean relative error，MRE）3種誤差作為評價指標(biāo)。

表1 各空氣污染物缺失統(tǒng)計Table1 Statistics on the missing of air pollutants

3.1 不同缺失率下的數(shù)據(jù)補全

為探究不同缺失率下的表現(xiàn)，采用5%，10%，15%，20%，25%，30%6種缺失率進行3種缺失模式的隨機挖空，訓(xùn)練集和測試集的比例為8:2。其中，時間連續(xù)缺失，定義為某一站點11個時刻連續(xù)缺失，為時間視圖范圍大??；空間連續(xù)缺失，定義為36個站點連續(xù)缺失，為所有空間站點數(shù)量。

考慮數(shù)據(jù)選取的隨機性會導(dǎo)致補全難度不一，本實驗中，每一高缺失率挖空點位集合包含了前一低缺失率挖空點位集合。

由表2可知，本文方法對PM2.5時空數(shù)據(jù)3種缺失模式的數(shù)據(jù)補全是有效且實用的。對應(yīng)PM2.5數(shù)值在 1.000～1 000.000 μg·m-3時，RMSE＜25，MAE＜14，MRE＜0.17。時間或空間連續(xù)缺失的各項誤差指標(biāo)值明顯較大，說明連續(xù)缺失模式的補全難度更大。

表2 不同缺失率下3種缺失模式數(shù)據(jù)補全結(jié)果（PM2.5）Table2 Data completion results for different missing modes under different missing rates（PM2.5）

3.2 不同視圖的數(shù)據(jù)補全

為探究不同視圖的數(shù)據(jù)補全效果，在15%缺失率下，分別進行時間視圖（EMA）、空間視圖（OK）和時空關(guān)聯(lián)視圖（Non-convex）的單視圖數(shù)據(jù)補全效果實驗，并與本文的時空多視圖方法進行對比。

由表3知，本文方法的各項指標(biāo)均優(yōu)于單視圖數(shù)據(jù)補全。此外，單視圖補全方法在極端情況下無法補全數(shù)據(jù)。如：時間視圖基于時間序列數(shù)據(jù)，在時間連續(xù)缺失情況下無法進行補全；空間視圖基于空間數(shù)據(jù)，在空間連續(xù)缺失情況下無法進行補全；時空視圖無法基于全空行或全空列求取相似低秩矩陣，在全空行和全空列下無法進行補全。

表3 不同視圖數(shù)據(jù)補全結(jié)果（PM2.5）Table3 Data completion results for different views（PM2.5）

由圖3可知，時空多視圖方法最優(yōu)；3個單視圖中，空間視圖最差，時間視圖最優(yōu)，時空視圖介于中間，說明各視圖對方法貢獻度由高到低依次為時間視圖、時空視圖和空間視圖。

3.3 與其他方法的比較

圖3 在隨機缺失模式下不同視圖的數(shù)據(jù)補全結(jié)果（PM2.5）Fig.3 Data completion results for different views in random missing mode（PM2.5）

將STMV-BP和以下5種典型數(shù)據(jù)補全方法的表現(xiàn)效果進行了比較，結(jié)果見表4。5種典型數(shù)據(jù)補全方法為：（1）ARIMA，基于過去短期數(shù)據(jù)預(yù)測缺失數(shù)據(jù)，不斷迭代預(yù)測缺失數(shù)據(jù)；（2）STKNN，基于k個最鄰近的時間和空間數(shù)據(jù)，取其平均進行缺失數(shù)據(jù)的填補；（3）NonCVX，基于EMA和OK進行局部初始化，采用非凸矩陣完備算法進行缺失數(shù)據(jù)填補；（4）EMA+OK，基于EMA和OK計算時間視圖和空間視圖，取其平均進行缺失數(shù)據(jù)的填補；（5）STMV-R，采用本文提出的時空多視圖作為特征，采用線性回歸融合時空多視圖特征。

由表4知，在各種缺失場景下，STMV-BP的補全效果明顯優(yōu)于各對比方法；考慮時間和空間的多角度方法，明顯優(yōu)于僅考慮單視圖的方法；僅考慮時間線性關(guān)系的ARIMA補全效果最差，且在時間連續(xù)缺失場景下無法有效進行數(shù)據(jù)補全。

考慮時空多視圖和對多視圖特征進行融合，可明顯提高數(shù)據(jù)的補全效果；BP神經(jīng)網(wǎng)絡(luò)3種缺失模式下的補全效果均明顯優(yōu)于線性回歸。

圖4 不同方法數(shù)據(jù)補全誤差的累積分布函數(shù)（PM2.5）Fig.4 Cumulative distribution function of data completion error for different methods（PM2.5）

各方法數(shù)據(jù)補全誤差的累積分布函數(shù)（cumulative distribution function，CDF）見圖4，從圖4中可發(fā)現(xiàn)，STMV-BP在誤差為0處導(dǎo)數(shù)最大且最先趨近于1，進一步說明本文方法整體表現(xiàn)更優(yōu)。

表4 不同方法數(shù)據(jù)補全結(jié)果（PM2.5）Table4 Data completion results for different methods（PM2.5）

3.4 各空氣污染物的數(shù)據(jù)補全

在15%缺失率下，各污染物數(shù)據(jù)補全效果見表5。各污染物的誤差項均在合理范圍內(nèi)，以隨機缺失為例，MRE普遍在0.16以下，說明本文提出的STMV-BP方法在空氣污染物數(shù)據(jù)補全上具有適用性。

各空氣污染物平均相對誤差和缺失率對照見圖5。通過比較發(fā)現(xiàn)，除SO2外，其他污染物的平均相對誤差與缺失率普遍呈對應(yīng)關(guān)系，即缺失率越高，數(shù)據(jù)補全的效果越差，缺失率越低，數(shù)據(jù)補全的效果相應(yīng)較好。而SO2數(shù)據(jù)缺失率的平均相對誤差明顯偏高，可能原因為SO2數(shù)據(jù)的時空變換規(guī)律較其他污染物更難捕捉。

3.5 模型參數(shù)

圖5 各空氣污染物平均相對誤差和缺失率對照Fig.5 Comparison of mean relative error and missing rate for different air pollutants

本文在BP神經(jīng)網(wǎng)絡(luò)的具體構(gòu)建中，采用4層網(wǎng)絡(luò)結(jié)構(gòu)，其中第1層為包含3個神經(jīng)元的輸入層，第2層為包含6個神經(jīng)元的第1隱含層，第3層為包含2個神經(jīng)元的第2隱含層，第4層為包含1個神經(jīng)元的輸出層；將Sigmoid函數(shù)作為激活函數(shù)，均方誤差（mean square error，MSE）作為損失，Levenberg-Marquardt算法作為優(yōu)化方法，對數(shù)據(jù)集預(yù)先進行歸一化處理。

時空多視圖BP神經(jīng)網(wǎng)絡(luò)模型受不同參數(shù)的影響。在時空多視圖構(gòu)建中，指數(shù)移動平均方法的threshold和α，非凸矩陣完備方法的p取不同值時的補全效果見圖6。

本文最終選擇參數(shù)threshold=5，此時補全效果已較優(yōu)，當(dāng)參數(shù)取更大值時，所需數(shù)據(jù)量更大且會影響數(shù)據(jù)處理效率；α=0.85，此時補全效果最優(yōu)；p=0.5，此時補全效果達到最優(yōu)。

4 結(jié) 論

圖6 時空多視圖參數(shù)影響Fig.6 Impact of parameters for spatio-temporal multi-view

基于指數(shù)移動平均、普通克里金和非凸矩陣完備，提出了面向城市空氣質(zhì)量數(shù)據(jù)補全的時空多視圖BP神經(jīng)網(wǎng)絡(luò)模型。由于集合了從不同視角出發(fā)具有互補性的不同視圖，本方法可在數(shù)據(jù)隨機缺失、時間連續(xù)缺失和空間連續(xù)缺失3種模式下有效補全數(shù)據(jù)，與單視圖補全方法相比，不僅提升了隨機補全的效果，而且能夠在極端缺失情況下補全數(shù)據(jù)。由于采用BP神經(jīng)網(wǎng)絡(luò)映射非線性關(guān)系，本方法能夠更好地擬合時空多視圖關(guān)系，各誤差指標(biāo)明顯較線性方法低。實驗證明，該方法能夠有效補全各空氣污染物數(shù)據(jù)。本文研究成果可為城市空氣質(zhì)量數(shù)據(jù)補全工作提供方法支持，研究思路可供時空數(shù)據(jù)挖掘參考。下一步將重點研究本方法對其他類型時空缺失數(shù)據(jù)的補全效果。

表5 各空氣污染物數(shù)據(jù)補全結(jié)果Table5 Data completion results for different air pollutants