基于時空融合的缺失值填補算法

蘇思凡，竹 翠，朱文軍，趙楓朝

(北京工業大學 信息學部，北京 100022)

0 引 言

空氣污染物預報的方法可分為機理模型和機器學習模型。前者基于空氣動力學知識對污染物的產生和擴散建模。代表算法有CMAQ和WRF-Chem。其理論性強，能對污染成因進行診斷，但需要使用者有豐富的專業知識，時間與計算資源開銷大，還需詳細的污染源資料。為尋找更加簡便的方式，自回歸移動平均[1]、線性回歸[2]、支持向量機[3]、神經網絡[4,5]等機器學習模型被廣泛嘗試。長短時記憶神經網絡[6]、深度置信網絡[7]等深度學習模型也用于該領域。此類方法不用詳述特征與目標間的關系，預測結果可精確到點位，但無法對污染成因給出原理上的解釋。

從數值角度看，PM2.5濃度是一組時間序列。此外，其它地區的污染物也在不斷擴散，任何地區的PM2.5濃度都不能“獨善其身”，勢必會受到相鄰地區的影響[8]。本文分別從時、空兩個維度進行預測，再將預測結果融合，對不同氣象條件下的PM2.5有針對性地回歸。

由于傳感器故障，數據集往往存在大量缺失數據。通過填補缺失值彌補其對預測帶來損失也是本文研究重點之一。盡管缺失值填補已有一些常用的方法[9,10]，但這些方法很少能充分利用數據的時空相關性。因此，本文結合時空背景提出了一種缺失值填補算法，輔助PM2.5濃度預測模型進一步精提升精度。

1 模型原理簡述

1.1 人工神經網絡(ANN)原理簡述

ANN一般由3類結構組成：輸入層、隱藏層和輸出層。每一層包含一個或多個神經元。神經元將前一層的輸出作為當前層的輸入，將這些輸入乘以不同的權值，可以實現對輸入信號的放大與抑制。然后利用激活函數將計算后的值壓縮到某一區間后傳遞給下一層的神經元。數據經過層與層之間連接最終傳導到輸出層，由輸出層神經元將結果送出。

對ANN的訓練就是根據輸出結果與真實值計算誤差，然后根據誤差從后向前不斷調整每層神經元之間權值的過程。這個過程被稱為“誤差的反向傳播”。該過程利用梯度下降的方法對參數進行調整，使模型的輸出值與真實值之間的誤差盡可能減小，從而不斷優化模型，獲得理想的輸出結果。

1.2 長短時記憶神經網絡(LSTM)原理簡述

LSTM[11]是循環神經網絡的變體。與ANN不同的是，LSTM中的神經元除了可以把輸出傳遞到下層神經元，還引入了時序的概念，可以把當前的輸出作為自身下一時刻的輸入，能夠讓信息在網絡中循環。LSTM與ANN相比，最大的優勢在于它能將網絡之前獲得的信息應用在當前神經元的計算中。LSTM的神經元細胞中包含“輸入門”、“忘記門”和“輸出門”，這種結構能有選擇性地記住和忘記曾經獲得的信息，并解決長期依賴問題。正是由于這些特點，LSTM可以從以前的數據中提取知識，從而為當前的預測提供幫助。

1.3 模型樹(MT)原理簡述

MT[12]是將決策樹和線性回歸算法進行融合的一種模型。它首先利用樹構建算法尋找當前最佳的切分屬性與最佳的切分點，不斷將數據集進行二元切分。與普通的決策樹不同的是，MT中的每個葉子結點并不是所在分支對應的標簽，而是基于劃分后的子樣本集所建立的線性回歸模型。與線性回歸算法不同的是，MT的結果是分段線性的，而線性回歸算法是全局線性的。當數據集較大且問題的影響機制復雜時，使用普通線性回歸算法顯然不夠靈巧，不能很好地對全體數據進行擬合。而MT將全體數據切分成多組易建模的數據后再回歸，能夠更好地處理非線性問題，同時提升了算法的可解釋性和預測準確性。

2 基于時空融合的預測模型

一個地區的氣象數據和空氣污染物數據通常由分布在該地區的多個氣象和空氣質量監測站獲取并發布的。本模型也以站點為粒度對PM2.5濃度進行預測。針對每一個站點，分別建立一套預測模型。每套模型的結構相同，但分別依賴各自站點的數據進行訓練與預測。

該模型由3個子模型構成：基于時間維度的LSTM預測器、基于空間維度的ANN預測器和基于時空融合的MT預測器，如圖1所示。其中LSTM負責從時間序列中獲取PM2.5的變化規律，從時間維度進行預測；ANN負責學習站點間PM2.5的影響機制，從空間維度進行預測；而MT負責將時、空間維度的結果整合起來，獲取最終的預測結果。

圖1 基于時空融合的PM2.5濃度預測總體架構

2.1 基于時間維度的LSTM預測器

LSTM可以從歷史數據中獲取變量在時間維度上呈現出的趨勢性和周期性規律。本文利用LSTM通過輸入某一站點前12小時的逐小時PM2.5濃度以及氣象數據(包含6個變量：天氣、濕度、氣壓、風向、風速、溫度)，對該站點下一小時的PM2.5濃度進行多因子的時間序列預測。

在LSTM模型的結構設置上，本文使用了一層LSTM疊加dropout，再疊加一層全連接層的結構來進行訓練與測試，如圖2所示。模型使用mae作為損失函數，adam為優化算法。在訓練前，首先要將每12個相鄰小時的氣象和PM2.5濃度數據處理為一組時間序列，作為訓練數據輸入模型。輸入數據首先通過LSTM層，提取出氣象和PM2.5隨時間變化的時序信息與規律。由于LSTM的學習能力十分強大，為防止過擬合而出現模型泛化性不好的問題，本文在LSTM層結束后加入了dropout，利用以一定比例隨機斷開神經元間連接的策略對網絡進行正則化。經過多次實驗，本文選擇以0.01作為dropout的參數，即在模型訓練過程中，會有1%的神經元的連接會被隨機斷開。經過dropout處理過的數據將通過全連接層進行特征提取與特征映射，根據輸入的氣象和PM2.5數據計算出下一小時的PM2.5濃度值。

圖2 基于時間維度的LSTM預測器結構

2.2 基于空間維度的ANN預測器

某地區的當前的PM2.5濃度值除了會受到本地氣象擴散條件和前期PM2.5水平的影響，周圍地區的PM2.5也會通過大氣運動傳輸到本地。因此，對某一地區PM2.5濃度的預測還要考慮空間維度上存在的相關性。由于PM2.5所呈現的空間相關性除了和距離有關，還會受到空氣對流情況的影響，度量起來較為復雜。為了降低特征工程的難度并獲得一種可靠的空間預測模型，本文利用ANN來模擬周圍地區對本地PM2.5濃度的影響。

基于空間維度的ANN預測器使用了與“基于時間維度的LSTM預測器”類似的序貫模型。該網絡由3個全連接層和dropout組成。dropout的參數設為0.01，全連接層的激活函數選用relu。該模型的輸入數據由本地與其鄰居站點前一小時的PM2.5濃度、本地站點前一小時的風速與風向數據組成。

對于圖3而言，星型代表當前研究的站點，我們以它作為中心，向四周發散出8條射線，每條線間呈45度角。這8條射線將研究站點的周圍劃分成8個區域。對于每一塊區域，我們取離研究站點最近的站點為鄰居站點。按照上述規則，每一個站點最多有8個鄰居站點，分別來自于不同的方向。圖中圓圈代表鄰居站點，方塊代表非鄰居站點。當研究站點前一小時的風向為東北風時，我們通常認為鄰居站點A和鄰居站點B對當前時刻研究站點的PM2.5濃度水平有著較為直接的影響。我們把當前研究站點及其所有鄰居站點前一小時的PM2.5濃度數據以及當前研究站點前一小時的風速風向數據輸入ANN中，試圖讓網絡學習到類似以上風對站點之間PM2.5濃度值的影響機制，從而從空間角度對PM2.5水平進行回歸。如果直接將所有站點的PM2.5濃度數據全部輸入到神經網絡中，輸入的數據量過大，為了得到相對穩定和準確的效果，勢必要增大神經網絡的規模，這將會提升調整參數的難度，也會導致訓練的時間過長，甚至模型震蕩難以收斂。

圖3 鄰居站點

2.3 基于時空融合的MT預測器

有文獻指出，如果各個神經網絡的誤差互不相關，將這些神經網絡作為子預測器進行集成后，總體的預測效果與泛化能力將獲得明顯提升。為了盡量提升子預測器的獨立性，本文分別利用不同的特征構造不同的樣本，并利用不同類型的神經網絡從不同的角度進行預測。現在我們已分別獲得了從時間和空間角度的PM2.5濃度預測結果。當對這兩個結果進行融合時，需要考慮如何進行權重分配的問題。從問題本身來看，當某一站點所處區域的空氣流動性較差時，鄰居站點的PM2.5很難通過大氣運動傳輸到本地。此時，利用ANN所得的結果可能不如利用LSTM所得的結果影響力大。相反，當該站點所處區域的風速較大時，LSTM所得的局部預測結果就應該弱化，而ANN所得的預測結果應該適當得到加強。因此對于這種較為復雜的問題，時、空預測結果的權值不應一概而論，應由當時的氣象擴散條件決定。時空融合預測器所使用的MT算法可以很好地解決這個問題，其示意圖如圖4所示。MT中每個線性回歸方程所確定的參數都適應于當前的氣象條件，比對整體樣本統一進行一次回歸有更好的針對性與適應性。

圖4 基于時空融合的MT預測器結構

3 基于時空融合的缺失值填補算法

缺失數據會使數據中蘊含的規律更難把握，阻礙模型發現與學習特征之間的影響機制。在實際建模中，缺失值的危害主要體現在數據量不足而導致的模型的泛化性差，輸出結果不可靠等方面。此外，由于傳感器硬件故障或通信失敗等情況，數據集經常會存在缺失數據。且本文所使用的數據集是精確到小時粒度的，如果以上故障排除不能及時排除將導致連續缺失值的出現，嚴重降低了數據質量。為了解決傳感器網絡缺失的問題，并進一步提高模型的預測性能，本文提出了一種基于時空融合的缺失值填補算法。該算法的整體思路是利用站點自身的時間關系和站點間的空間關系找到與缺失值最相近的4個參考值，然后將相關系數作為這4個參考值對應的權重。對于連續性變量，將4個參考值進行加權平均后的結果作為填充值。對于非連續性變量，將4個參考值中權重最大的值作為填充值。本方法是對熱卡填充法進行優化后的一種算法，與普通熱卡填充法相比可以有更多的參考值，并從時間和空間兩個維度尋找相似的樣本，提高了填補算法的魯棒性與準確性。

3.1 數據集結構說明

對于n個數據采集站點在相同時間粒度下獲得的觀測值，可以以圖5的方式進行組織排列。其中，每一列代表一個數據采集站點獲得的觀測值的時間序列，每一行代表某一時刻所有站點的觀測量的取值情況。對于記錄中的缺失值，用符號“na”表示。

圖5 數據集結構

對于圖5，我們以3號站點t時刻的缺失值作為研究對象，對其進行填補。

3.2 時間維度

對于時間序列，變量歷史的情況會對當前和未來時刻的情況產生影響，當前時刻的變量是在前一時刻的基礎上發展而來的。特別是對于本文使用的小時粒度的氣溫、氣壓等氣象數據和PM2.5濃度數據，這些變量在一小時內很少會出現大幅度跳變的情況。所以對于時間間隔較小的時間序列，可以簡單地認為與當前時刻與前后兩個時刻的指標水平比較相似。因此，對于所研究的t時刻站點的取值序列，我們認為與t-1時刻和t+1時刻的取值序列最為相近，其相似程度選用皮爾遜相關系數進行量化，公式如式(1)所示。R>0表示變量X和Y間存在正相關關系。R越大，表明兩組變量的正相關程度越強。當R=1時，X和Y呈完全正相關線性關系。反之亦然。R=0表示兩者間無線性關系。本文利用皮爾遜相關系數分別計算t-1時刻序列與t時刻序列的相似度cor1、t時刻序列與t+1時刻序列的相似度cor2

(1)

val1和val2可以看作是與當前缺失值最相似的兩個參考值。本算法利用序列的相似度cor1和cor2分別代表val1與當前缺失值、val2與當前缺失值的相似度，并以此分別作為val1、val2的權值，參與3.4部分時空融合填補的算法。

3.3 空間維度

在空間維度，由于數據采集站點散布在研究區域的不同位置，站點間的距離與位置關系比較復雜，每個站點所處的環境也各不相同，因此很難通過量化這些外部因素的方式來判斷空間中與哪些站點與當前站點最相似。為了找到衡量空間相似度的標準，與時間維度不同，本文直接從各站點采集到的觀測值序列的角度出發，通過計算序列間的相關系數矩陣，從中找到與當前站點在觀測值表現上最相似的其它站點。

對于圖5，通過計算n列兩兩之間的相關系數，我們可得到一個大小為n*n的相關系數矩陣。易知該矩陣是對角線為1的對稱矩陣。由于皮爾遜相關系數的絕對值越大表明兩組變量相關程度越強。因此，遍歷相關系數矩陣中的每一列，排除對角線元素，取其中絕對值最大的前兩個元素，這兩個元素對應的行標簽就是與當前站點相似度最高的兩個站點。對于圖5中當前研究的缺失值，假設通過上述方法找到與其所在的3號站點在相似度最高的兩個站點分別為1和n-1號站點，那么與當前缺失值同一時刻下的val3和val4被認為是空間維度下與當前缺失值相似度最高的兩個參考值。通過相關系數矩陣找到的當前列絕對值最大的前兩個元素cor3和cor4分別對應為val3和val4的權值，參與3.4部分時空融合填補的算法。

3.4 時空融合

經過3.2和3.3兩個部分，對于當前缺失值，我們已從時間和空間維度找到了與其相似度較高的4個參考值，并確定了它們各自的權重。本算法從時間和空間維度共選取4個參考值的原因在于，本課題所用數據中的缺失值較多，如果在時、空維度分別僅選取一個參考值，經常會出現兩個參考值同時都是缺失值的情況，導致參考值不能發揮作用，缺失值依然不能被填補的情況發生。為了在保證填補精度的前提下盡量提高填補效率，本算法分別從時空維度各選取兩個參考值，以降低參考值同時為缺失值的概率。

在對缺失值進行填補時，需要區分兩種情況：一是待填補值是連續型變量，二是待填補值是非連續性變量。對于前者，可使用加權平均的思想。對于后者，由于待填補值非數值，且只能是該特征屬性集合中的某一個，因此加權平均的方式并不適用，這里選擇權值最大的參考值作為填補值。需要說明的是，當某一個或多個參考值為缺失值時，需要將它和它所對應的權值分別從參考值集合和權重集合中剔除，使其不參與時空融合缺失值填補的計算。此時，參考值的總個數從4個退化為不足4個。當且僅當4個參考值全缺失時，當前缺失值無法填補。

填補值fill_val的計算方式如式(2)所示。并利用同樣的方法對整個數據集缺失值進行填補

(2)

4 實驗及分析

本文共設計了3組實驗：時空融合與單維度預測的效果對比、不同時空融合預測器的效果對比、兩種缺失值填補方法的效果對比。第1組實驗是為了驗證本文提出的基于時空融合的預測框架比從時間或空間某一維度進行獨立預測的優越性。第2組實驗是為了驗證在時空融合預測框架的前提下，MT比其它主流回歸算法更適合作為當前問題下的時空融合預測器。由于本文所使用的原始數據集是精確到小時粒度的，且數據缺失程度較高，我們也未能找到其它可以準確參考的小時級數據源，因此本文利用天級粒度的氣象后報與PM2.5濃度后報對所用數據集進行手工近似填補，試圖盡量排除缺失值帶來的不良影響。前兩組實驗均是在手工近似填補后的數據集的基礎上進行的。第3組實驗利用相同的預測模型，分別基于手工近似填補缺失值的數據集和時空融合填補缺失值的數據集進行預測，旨在驗證本文所提出的基于時空融合填補缺失值的有效性。

4.1 實驗數據說明

本文以北京地區為研究對象對PM2.5濃度進行預測。北京共有16個行政區縣和35個空氣質量監測站。本文所使用的數據源可以分為3類：各區縣氣象數據、各空氣質量監測站PM2.5濃度數據、各空氣質量監測站經緯度數據。根據第3類數據可以得到站點間的鄰居關系。前兩組數據的時間跨度為2014年5月1日至2015年4月30日，共計12個月，均以小時為單位進行記錄。但由于缺失數據的存在，前兩組數據的樣本量無法達到小時級粒度下應有的數量。對于以下實驗，我們把前11個月(2014年5月至2015年3月)的數據作為訓練集，把第12個月(2015年4月)的數據作為測試集，對35個站點未來一小時的PM2.5濃度進行預測。

4.2 評價指標

本課題屬于回歸類問題。因此，本文選用了回歸問題常用的兩種評價指標：均方根誤差(RMSE)和決定系數(R2)，對實驗效果進行衡量。這兩個評價指標的公式如式(3)、式(4)所示

(3)

(4)

4.3 時空融合與單維度預測的效果對比

本文第2部分提出了基于時空融合的PM2.5預測模型。該模型分別使用LSTM和ANN從時空兩個角度獨立預測，再將預測結果融合起來用模型樹回歸。為論證該方法的有效性，本實驗將時空融合預測(LSTM+ANN+MT)的效果與時間維度單模型預測(LSTM)和空間維度單模型預測(ANN)的效果進行對比。

表1對時空融合預測模型的整體表現進行了統計。可以看出經過時空融合預測后35個站點的平均指標均好于時空單獨預測。且利用LSTM+ANN+MT進行預測后，大部分站點比時間或空間單維度預測的最好效果還有提升。可見，本文提出的時空融合預測模型是有效的。

4.4 不同時空融合預測器的效果對比

本文2.3部分介紹了基于時空融合的MT預測器。為證明MT作為時空融合預測器的有效性，對于相同的數據集，本文選擇了回歸算法中較為常用的線性回歸(LR)和支持向量回歸(SVR)分別替代MT作為時空融合預測器。

表1 融合預測與單獨預測效果統計

表2對3個模型的整體表現進行了統計。可以看出MT作為時空融合預測模型時，35個站點的平均RMSE和R2均好于LR和SVR。可見，MT在回歸前先對數據集分組，然后分別在每個類中進行回歸，這種方式可以更好地適應樣本間差異較大的數據集，使每次回歸考慮的樣本量縮小，提升回歸算法的針對性和擬合效果。因此利用MT作為時空融合預測模型是有效的。

表2 不同時空融合預測器效果統計

另外，即使是使用LR或SVR作為時空融合預測器，與表1對比，35個站點的平均指標也都優于時空單獨預測，進一步驗證了本文提出的基于時空融合的預測框架的有效性。

4.5 兩種缺失值填補方法的效果對比

本文利用第3部分提出的算法，對數據集中各屬性分別進行缺失值填補。填補前后的數據集缺失程度對比可見表3。從中可以看出，缺失值填補后，缺失值大量減少，可為模型訓練提供更加充足的樣本。其中，填補后數據集依然不完整是由于3.4中所提到的4個參考值均為缺失值的原因所導致的。

此外，對于相同的原始數據集，本文選用兩種方式對缺失值進行填補：利用天級后報數據進行手工近似填補(下稱“手工填補”)、本文提出的基于時空融合進行缺失值填補(下稱“時空填補”)。對于這兩個填補后的數據集，本文利用相同的時空融合預測模型(LSTM+ANN+MT)進行預測，并從RMSE和R2的角度對預測效果進行對比，如圖6、圖7所示。其中，灰色標記代表手工近似填補，黑色標記代表時空融合填補。并以不同形狀對不同預測方式進行區分。

表3 缺失值填補前后數據集缺失程度對比

圖6 兩種缺失值填補方法在RMSE上的預測效果對比

圖7 兩種缺失值填補方法在R2上的預測效果對比

對于圖7，可以看出絕大部分的黑色標記都處于上方，表明利用時空填補方式可以使R2升高。此外，“×”所代表的基于時空融合的預測框架分別處于圖6和圖7中同類填補算法的最下方與最上方，進一步驗證了本文提出的基于時空融合的預測框架的有效性。

表4給出了35個站點整體預測情況的統計表。從中可以看出絕大部分站點在填補后預測水平都有提升。RMSE在原有基礎上進一步降低，R2進一步提高。可見，本文提出的基于時空融合的缺失值填補算法對時空融合預測模型起到了很好的輔助作用。

表4 兩種缺失值填補方法的預測效果對比

5 結束語

本文提出了一種基于時空融合和缺失值填補的PM2.5濃度預測算法。該算法首先從時間和空間維度入手，分別利用LSTM和ANN進行獨立預測，再利用MT將時空預測結構進行融合。此外，還結合時空相關性對缺失值進行填補，輔助上述模型進一步提高預測精度。經過實驗可以得出以下結論：

(1)利用合適的模型，從時空角度獨立預測再融合的預測效果好于單模型單維度的預測效果；

(2)MT作為一種分段線性的回歸模型，比LR和SVR等全局線性的回歸模型更適合作為該場景下的時空融合預測器；

(3)本文提出的基于時空融合的缺失值填補法可以有效地提升數據集的完整性。并且利用該算法填補過的數據集進行PM2.5預測，準確度有進一步提高。