基于高階圖卷積網(wǎng)絡(luò)的城市空氣質(zhì)量推斷模型*

陳 杰 ，許鎮(zhèn)義 ，2

(1.中國科學(xué)技術(shù)大學(xué) 自動化系，安徽 合肥230026；2.合肥綜合性國家科學(xué)中心人工智能研究院，安徽 合肥230088)

0 引言

近年來，隨著經(jīng)濟(jì)的增長，環(huán)境問題也變得日益突出， 大氣污染問題正受到前所未有的關(guān)注和重視[1]。城市空氣中，如一氧化碳(CO)、碳?xì)浠?HC)、氮氧化物(NOx)、固體顆粒物(PM2.5、PM10)等污染物濃度與人們的身體健康息息相關(guān)[2-3]。 空氣質(zhì)量指數(shù)(Air Quality Index，AQI)是定量描述空氣質(zhì)量狀況的指數(shù)，其數(shù)值越大說明空氣污染狀況越嚴(yán)重，對人體健康的危害也就越大[4]。

為了及時準(zhǔn)確地反映空氣質(zhì)量及發(fā)展趨勢，需要精確的空氣質(zhì)量監(jiān)測設(shè)備。然而這些監(jiān)測設(shè)備的建設(shè)成本以及后期的維護(hù)費用是昂貴的[4]，所以，在一個城市不同區(qū)域安裝大量的空氣監(jiān)測設(shè)備是不現(xiàn)實的。 例如，北京市主城區(qū)面積約為1 381 平方公里，僅僅包含17 個空氣監(jiān)測站點，各個監(jiān)測站位置分布由圖1 所示。

圖1 北京市空氣質(zhì)量監(jiān)測站位置分布

這些年來，關(guān)于空氣質(zhì)量推斷的研究主要分為：基于知識驅(qū)動和基于數(shù)據(jù)驅(qū)動的兩類方法。 基于知識驅(qū)動的方法主要使用數(shù)學(xué)模型[5-6]和物理知識[7]，模擬空氣污染物的實際物理擴(kuò)散過程，通過計算仿真來推斷未觀測區(qū)域的空氣質(zhì)量分布。 這些方法都采用了一系列優(yōu)化技術(shù)，如相關(guān)性分析、聚類分析和多目標(biāo)優(yōu)化[8-10]。 為了達(dá)到穩(wěn)定狀態(tài)，仿真過程不僅需要復(fù)雜的系統(tǒng)編程，而且會消耗大量的計算能力。 模型簡化和穩(wěn)定性假設(shè)，也會使模型效率進(jìn)一步降低。

近幾年來，基于數(shù)據(jù)驅(qū)動的方法也在不斷發(fā)展,并在空氣質(zhì)量推斷方面取得了很好的效果。ZHENG Y 等[4]提出了一種基于Co-training 協(xié)同訓(xùn)練框架[11]的半監(jiān)督學(xué)習(xí)方法來推斷城市空間細(xì)粒度的空氣質(zhì)量。 但是該框架由兩個獨立的分類器組成，很難捕獲不同特征之間的復(fù)雜交互性。CHEN L等[12]采用K-近鄰策略來構(gòu)建外部數(shù)據(jù)源(如路網(wǎng)結(jié)構(gòu)、交通等)與空氣質(zhì)量分布之間的關(guān)系。 然而K的具體數(shù)值沒有確切的理論標(biāo)準(zhǔn)，最近的K個站點的特征不一定是最有效的，并且不同站點的影響程度可能會隨時間而變化，因此有可能導(dǎo)致出現(xiàn)非一致性的問題。 HSIEH H P 等[13]基于圖標(biāo)簽節(jié)點傳播[14]，設(shè)計了一種基于親和關(guān)系的空氣質(zhì)量推斷模型(the Affinity-based AQI Inference Model，AQInf)。 然而由于其自身模型性能的限制，AQInf 很難捕獲圖節(jié)點之間的反映空氣質(zhì)量時空變化趨勢的關(guān)聯(lián)程度。

本文主要想解決的問題是：如何基于城市已建立的有限數(shù)量的空氣質(zhì)量監(jiān)測站，來精確地推斷未監(jiān)測區(qū)域的空氣質(zhì)量分布。 因而將城市空間細(xì)粒度空氣質(zhì)量預(yù)測問題轉(zhuǎn)化為城市時空圖(Spatio-Temporal Graph，STG)預(yù)測問題，其中，時空圖中的每個節(jié)點表示城市中的每個子區(qū)域(1 km×1 km)，并且每個節(jié)點關(guān)聯(lián)隨時間變化的空氣質(zhì)量數(shù)值。 任務(wù)的挑戰(zhàn)性主要基于以下幾個方面：

(1)空氣質(zhì)量分布受到很多復(fù)雜的外部環(huán)境因素的影響(例如天氣、交通、土地用途等)。 對于商業(yè)中心等交通流量大的區(qū)域，空氣質(zhì)量往往比公園湖泊等區(qū)域差。 這樣會導(dǎo)致由于地理因素影響而出現(xiàn)不光滑數(shù)據(jù)值，很難通過基于插值的方法來精確推斷未監(jiān)測區(qū)域的空氣質(zhì)量。 例如，圖1 中的S5、S6監(jiān)測站點和S11、S13 監(jiān)測站點，在地理位置上很接近，但是PM10、AQI 指數(shù)分布卻常年差異很大，具體分布如圖2 所示。 造成這個現(xiàn)象最可能的原因是，空氣質(zhì)量較優(yōu)的站點，靠近公園或者湖泊，而空氣質(zhì)量較差的站點，處于商業(yè)中心或主干道附近，交通量較大。

(2)空氣質(zhì)量分布具有時空交互性。 在同一時刻內(nèi)， 不同空間區(qū)域內(nèi)的空氣質(zhì)量分布相互影響；在同一空間內(nèi)，相同節(jié)點在不同時刻的AQI 數(shù)值相互關(guān)聯(lián)。

1 問題描述

1.1 城市時空圖

基于城市地理信息，將該城市區(qū)域劃分成不相交的網(wǎng)格(1 km×1 km 的子區(qū)域)，并且每一個網(wǎng)格與隨時間變化的AQI(Air Quality Index)相關(guān)聯(lián)，如圖3 所示。

將每個劃分的網(wǎng)格看作一個節(jié)點，構(gòu)建時空圖G=(V；ε；A)，其中，V是節(jié)點的集合，|V|=N；ε 是邊的集合，表示節(jié)點之間的連通性；A∈Rn×n為圖G 的鄰接矩陣。 時空圖中，V可以表示為V=L∪U，其中，L表示為已建立空氣質(zhì)量監(jiān)測站的節(jié)點集合，U表示為未建立空氣質(zhì)量監(jiān)測站的節(jié)點集合。將L中的節(jié)點作為標(biāo)簽節(jié)點，U中的節(jié)點作為未標(biāo)簽節(jié)點。時空圖的空間相關(guān)性通過同一時刻下節(jié)點之間的地理距離特征來反映，時間相關(guān)性通過同一節(jié)點相關(guān)聯(lián)的不同時刻AQI 數(shù)值來反映。

圖2 不同監(jiān)測站點的PM10，AQI 實時分布

圖3 城市區(qū)域網(wǎng)格劃分

1.2 數(shù)據(jù)收集

收集與空氣質(zhì)量分布相關(guān)的外部影響因素，該外部影響因素主要包括氣象數(shù)據(jù)、子區(qū)域POIs(Pointof-Interests)數(shù)目以及區(qū)域路網(wǎng)信息(如各區(qū)域高速路長度，主干道長度等[15]。 這些數(shù)據(jù)特征集合表示為F={f1，f2， …，fm}，fm表示收集到的第m個特征。POIs 類型如表1 所示。

表1 POIs 類型列表

1.3 空氣質(zhì)量分布預(yù)測問題

問題：在{Ti|i=1，2，…，n}時刻，基于時空圖G=(V；ε；A)，利用已觀測區(qū)域節(jié)點L的空氣質(zhì)量分布fl∈R|L|×P來預(yù)測任意未觀測區(qū)域U的空間細(xì)粒度空氣質(zhì)量分布fu∈R|U|×P。 空氣質(zhì)量分布預(yù)測任務(wù)主要是學(xué)習(xí)一個函數(shù)F(·)，將輸入fl∈R|L|×P映射到輸出fu∈R|U|×P：

其中，|L|為標(biāo)簽節(jié)點的數(shù)目，|U|為未標(biāo)簽節(jié)點的數(shù)目，P是每個節(jié)點特征的數(shù)目。

2 研究方法

2.1 圖卷積網(wǎng)絡(luò)

通用化卷積神經(jīng)網(wǎng)絡(luò)以處理任意圖結(jié)構(gòu)數(shù)據(jù)是近些年來研究的熱點之一。BRUNA J 等[16]通過將濾波器和信號轉(zhuǎn)換到傅里葉域，然后將相乘結(jié)果返回離散域從而實現(xiàn)圖卷積。 信號變換是通過與圖拉普拉斯圖的特征向量相乘來實現(xiàn)的，需要對稱拉普拉斯算子的二次特征分解。 然而，特征分解的低秩近似可以使Chebyshev 近似多項式來計算以降低計算復(fù)雜度[17-18]。 KIPF T N 等[19]通過局部頻譜濾波器的一階近似，設(shè)計了一種圖卷積網(wǎng)絡(luò)用于半監(jiān)督學(xué)習(xí)。

作為圖卷積的典型應(yīng)用，SUN J 等[20]設(shè)計一種多視圖圖卷積網(wǎng)絡(luò)用于預(yù)測城市不規(guī)則區(qū)域人群流量。 XU Z 等[21]設(shè)計一種時空圖卷積多融合網(wǎng)絡(luò)用于城市機動車排放預(yù)測。 CUI Z 等[22]設(shè)計了一種交通高階圖卷積循環(huán)神經(jīng)網(wǎng)絡(luò)，用于交通學(xué)習(xí)和預(yù)測。

2.2 高階圖卷積網(wǎng)絡(luò)的設(shè)計

最近的一些研究[23-25]表明，對于圖節(jié)點半監(jiān)督學(xué)習(xí)， 用隨機游走統(tǒng)計來學(xué)習(xí)表示的效果比較好，因為這樣可以保留圖結(jié)構(gòu)的完整性。 在特殊條件下，如果圖卷積網(wǎng)絡(luò)的第一層激勵函數(shù)是一個恒等函數(shù)，那么文獻(xiàn)[19]中的GCN 模型也可以學(xué)習(xí)隨機游走。 假設(shè)第一層的激勵函數(shù)是一個恒等函數(shù)，即：

該GCN 模型可以進(jìn)一步簡化為：

可以發(fā)現(xiàn)，在特定條件下的GCN 模型實際上做了一步隨機游走，可以將節(jié)點的信息傳遞到高階鄰居節(jié)點。 然而在Eq.2 中，信息特征在圖邊之間的傳遞呈縮減狀態(tài)(左乘標(biāo)準(zhǔn)化鄰接矩陣A?，輸出經(jīng)過了非線性激勵函數(shù)σ(·))。 因此，不能在原來的多層GCN 網(wǎng)絡(luò)上直接對A?取高階， 否則會導(dǎo)致網(wǎng)絡(luò)性能下降。 提高圖卷積網(wǎng)絡(luò)的階數(shù)，同時減少網(wǎng)絡(luò)的層數(shù)，并將每階網(wǎng)絡(luò)提取的特征進(jìn)行融合從而得到最終特征表示：

其中，W(k)為的權(quán)重參數(shù)，GCk表示圖卷積網(wǎng)絡(luò)提取到的k階特征。

2.3 空氣質(zhì)量推斷模型

首先定義時空圖的高階鄰居矩陣，時空圖的一階鄰居矩陣即為鄰接矩陣A，K階鄰接矩陣可以通過A 的k次冪得到。 如圖4 右邊所示，一個節(jié)點(五角星)的K階鄰居節(jié)點為圖中的黑色邊可以到達(dá)的節(jié)點。

將高階圖卷積網(wǎng)絡(luò)(High-order Graph Convolution，HGCN)應(yīng)用到環(huán)境監(jiān)測領(lǐng)域，并設(shè)計成空氣質(zhì)量推斷模型(Air quality inference model based on High-order Graph Convolution，HGCNInf)，其中，HGCN 用來捕獲空氣質(zhì)量分布的時空交互性，全連接神經(jīng)網(wǎng)絡(luò)(Fullyconnected Neural Network，F(xiàn)NN)用來提取復(fù)雜的外部影響因素特征，具體框架如圖4 所示。 為了使模型更具有合理性和解釋性，將Eq.7 的K階卷積操作修改為如下形式：

圖4 空氣質(zhì)量推斷模型框架(圖的右邊為高階圖卷積網(wǎng)絡(luò)構(gòu)圖)

將圖卷積提取的K階特征連接在一起，定義如下：

GCk∈RN×(K×P)是圖卷積網(wǎng)絡(luò)提取的K階特征集合，這樣就達(dá)到了提取輸入數(shù)據(jù)高階特征的目的。將GCk∈RN×(K×P)通過全連接網(wǎng)絡(luò)ffc1，通過反向傳播共同訓(xùn)練。 定義如下：

其中，空氣質(zhì)量分布特征GC ∈RN×P，可訓(xùn)練參數(shù)Wfc1∈RP×(K×P)。

然后用全連接網(wǎng)絡(luò)FNN 來捕獲外部影響因素特征，定義如下：

其中，fm∈RN×M為M個外部影響因素，可訓(xùn)練參數(shù)Wfc2∈RP×M，GCf∈RN×P。

將外部影響因素特征GCf∈RN×P與空氣質(zhì)量分布特征GC∈RN×P進(jìn)行融合：

其中，可訓(xùn)練參數(shù)Wfc3∈RP×P，得到總特征GCall∈RN×P。

進(jìn)一步，得到未標(biāo)簽節(jié)點的AQI 分布Fu：

其中，F(xiàn)l∈R|L|×P為標(biāo)簽節(jié)點的空氣質(zhì)量分布輸出，F(xiàn)u∈R|U|×P為未標(biāo)簽節(jié)點的空氣質(zhì)量分布輸出。

2.4 網(wǎng)絡(luò)損失和訓(xùn)練

最小化模型的輸出Fl和已知訓(xùn)練標(biāo)簽Y的交叉熵?fù)p失，交叉熵?fù)p失定義如下：

其中，*表示哈達(dá)瑪積；diag(Vl)表示一個對角陣，該對角陣的元素，當(dāng)(i，i)∈Vl時，設(shè)置為1，其余元素為0。

為了使圖卷積特征更具有穩(wěn)定性，在損失函數(shù)上面增加正則化項。 相鄰節(jié)點對特定節(jié)點的影響必須通過兩節(jié)點之間的所有節(jié)點傳遞，圖卷積操作提取到的不同階節(jié)點特征不應(yīng)該有太大變化。 為了限制相鄰階特征之間的差異，將基于L2 范數(shù)的卷積特征正則化項添加到損失函數(shù)中。 定義如下：

總的損失函數(shù)如下：

其中，參數(shù)α 用于控制正則化項的權(quán)重大小。 最小化模型損失函數(shù)對網(wǎng)絡(luò)進(jìn)行訓(xùn)練：

當(dāng)網(wǎng)絡(luò)訓(xùn)練結(jié)束時，從其未觀測區(qū)域的AQI 分布Fu中找到量化概率最高的值，并定義為未標(biāo)簽節(jié)點的AQI 預(yù)測值：

3 實驗結(jié)果

3.1 數(shù)據(jù)描述

采用北京市空氣質(zhì)量AQI 實時數(shù)據(jù)驗證了提出的模型， 該數(shù)據(jù)收集于北京市生態(tài)環(huán)境監(jiān)測中心，每小時收集一次，時間跨度從01/01/2015 到12/31/2015，包含17 個監(jiān)測站點，總計共17×8 760 條數(shù)據(jù)。 在實驗中，將北京主城區(qū)劃分成30×38 個網(wǎng)格子區(qū)域(1 km×1 km)，其中17 個子區(qū)域建有空氣質(zhì)量監(jiān)測站，剩余的1 123 個子區(qū)域的空氣質(zhì)量分布未知。

3.2 數(shù)據(jù)處理

本文的輸入數(shù)據(jù)主要包含北京市空氣質(zhì)量數(shù)據(jù)和復(fù)雜的外部影響因素數(shù)據(jù)。 采用one-hot 編碼來處理空氣質(zhì)量AQI 數(shù)據(jù)，采用Z-Score Normalization將外部影響因素數(shù)據(jù)進(jìn)行線性變換，使其結(jié)果映射到[0，1]范圍內(nèi)。 北京城市時空圖的鄰接矩陣根據(jù)節(jié)點間的距離，通過閾值化高斯核加權(quán)函數(shù)[26]來計算，計算公式如下：

其中，dist(pi，pj)2表示節(jié)點pi和pj之間的歐式距離，σ2和ε 是閾值用來控制矩陣Wij的分布和稀疏性。本實驗中，σ2和ε分別設(shè)置為10 和0.05。

3.3 實驗設(shè)置

3.3.1 實驗評價指標(biāo)

為了評價和比較不同模型之間的性能，采用兩個常用的指標(biāo)： (1)Mean Absolute Percentage Error(MAPE)；(2)Root Mean Squared Error(RMSE)。

3.3.2 基準(zhǔn)模型

通過以下基準(zhǔn)模型和提出的模型進(jìn)行了比較：

SVR：Support vector regression[27]。 SVR 利 用 監(jiān) 測站點的歷史數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)用來預(yù)測未標(biāo)簽節(jié)點的空氣質(zhì)量數(shù)值。

AQInf：Affinity-based AQI inference model[13]。 對于城市空氣質(zhì)量預(yù)測而言，基于目前的知識，這是最先進(jìn)的方法。 利用標(biāo)簽節(jié)點的歷史數(shù)據(jù)，結(jié)合外部影響因素，用來預(yù)測未標(biāo)簽節(jié)點的空氣質(zhì)量分布。

GCNInf：GCN-Based AQI inference model。 設(shè)計了一個如文獻(xiàn)[19]所描述的兩層圖卷積網(wǎng)絡(luò)。 類似于AQInf，利用標(biāo)簽節(jié)點的歷史數(shù)據(jù)，結(jié)合外部影響因素，用來預(yù)測未標(biāo)簽節(jié)點的空氣質(zhì)量分布。

所有的神經(jīng)網(wǎng)絡(luò)模型實施都是基于Pytorch version：1.1.0， 并且評價和訓(xùn)練都基于CPU Intel?CoreTMi5 -4210U CPU@1.70 GHz，NVIDIA GeForce 840M。

3.4 實驗結(jié)果分析

在實驗中，由于僅僅具有標(biāo)簽節(jié)點(已建立空氣質(zhì)量監(jiān)測站的子區(qū)域)的真實數(shù)據(jù)，對于未標(biāo)簽節(jié)點(沒有建立監(jiān)測站的子區(qū)域)而言，沒有真實數(shù)據(jù)可以驗證空氣質(zhì)量推斷模型的精度。 因此，采用交叉驗證的方法，在17 個標(biāo)簽節(jié)點中隨機選擇12個節(jié)點作為訓(xùn)練節(jié)點，產(chǎn)生12×8 760 個訓(xùn)練數(shù)據(jù)，剩余的5 個標(biāo)簽節(jié)點用來作為驗證節(jié)點，產(chǎn)生5×8 760個驗證數(shù)據(jù)。

3.4.1 訓(xùn)練效率

考慮不同卷積階數(shù)K對所提出來的HGCNInf推斷性能的影響。 圖5(a)、圖5(b)和表2 是HGCNInf模型取不同階數(shù)的性能比較。 對于北京市空氣質(zhì)量分布數(shù)據(jù)而言，在兩個指標(biāo)MAPE 和RMSE 上，當(dāng)高階圖卷積網(wǎng)絡(luò)階數(shù)K=3 時，模型HGCNInf 的推斷性能最好。 圖5(c)和表3 為HGCNInf 網(wǎng)絡(luò)取不同階數(shù)時，在每個時刻訓(xùn)練1 000 次所消耗的平均時間，可以發(fā)現(xiàn)，隨著HGCNInf 卷積網(wǎng)絡(luò)階數(shù)的增加，訓(xùn)練時間也隨之增加；此外，當(dāng)模型HGCNInf 考慮外部影響因素后，訓(xùn)練所消耗的時間將大幅度增加。

3.4.2 空氣質(zhì)量推斷模型結(jié)果對比

綜合模型預(yù)測精度和時間效率，選取K=3 作為HGCNInf 的階數(shù)，其預(yù)測精度如表4 所示。

圖5 不同階數(shù)的HGCNInf 性能和訓(xùn)練時間

表2 不同階數(shù)的HGCNInf 性能比較(%)

表3 不同階數(shù)的HGCNInf 訓(xùn)練時間(s)

表4 空氣質(zhì)量推斷模型的性能比較(%)

圖6 空氣質(zhì)量推斷模型的性能比較

對各個空氣質(zhì)量推斷模型分別訓(xùn)練三次，每次訓(xùn)練10 000 次，然后取平均值作為最終結(jié)果。 表4和圖6 展示了提出的模型HGCNInf 和其他基準(zhǔn)模型在北京市空氣質(zhì)量分布上兩個指標(biāo)的評價對比。 可以發(fā)現(xiàn)， 提出的模型HGCNInf， 無論在指標(biāo)MAPE上,還是指標(biāo)RMSE 上，都遠(yuǎn)遠(yuǎn)優(yōu)于其他基準(zhǔn)模型。另外，對比空氣質(zhì)量推斷模型在未加入外部影響因素和已加入外部影響因素在兩個評價指標(biāo)上的差異，同樣的模型，在考慮了外部影響因素影響后，對于城市空氣質(zhì)量分布預(yù)測精度的提升有著顯著的作用。因此，可以得出結(jié)論，復(fù)雜的外部影響因素與空氣質(zhì)量分布息息相關(guān)，這也很好地解釋了圖1 中地理位置很接近的監(jiān)測站點，在空氣質(zhì)量數(shù)值上卻常年差異較大的現(xiàn)象。

3.4.3 空氣質(zhì)量推斷模型的魯棒性

為了測試各個空氣質(zhì)量推斷模型對標(biāo)簽節(jié)點數(shù)目的敏感性，隨機丟棄相同的標(biāo)簽節(jié)點，觀察各個空氣質(zhì)量推斷模型預(yù)測性能的變化程度，具體結(jié)果如表5 和圖7 所示。

從表5 和圖7 中可以發(fā)現(xiàn)，當(dāng)逐漸減少標(biāo)簽節(jié)點的數(shù)目時，所提出的模型HGCNInf 的預(yù)測精度在評價指標(biāo)MAPE 和RMSE 上總是優(yōu)于各個基準(zhǔn)模型， 并且MAPE 和RMSE 增長速度相對緩慢。 這樣的結(jié)果表明，模型HGCNInf 的魯棒性要遠(yuǎn)遠(yuǎn)優(yōu)于各個基準(zhǔn)模型， 即使當(dāng)只有非常少的標(biāo)簽節(jié)點可用時，模型HGCNInf 也同樣是適用的。

表5 空氣質(zhì)量推斷模型的魯棒性(%)

圖7 空氣質(zhì)量推斷模型的魯棒性

4 結(jié)論

本文將城市空間細(xì)粒度空氣質(zhì)量預(yù)測問題轉(zhuǎn)化為時空圖預(yù)測問題，基于改進(jìn)的高階圖卷積網(wǎng)絡(luò)設(shè)計了一種有效的空氣質(zhì)量推斷模型HGCNInf。HGCNInf 不僅可以捕獲大氣污染物分布的時空相互作用，還可以提取復(fù)雜的外部影響因素特征，從而可以準(zhǔn)確地預(yù)測整個市區(qū)的細(xì)粒度空氣質(zhì)量分布。在接下來的工作中，將專注于通過并行化提高該模型的效率。此外，未來將尋求模型的更多應(yīng)用，特別是在城市地區(qū)的空氣質(zhì)量監(jiān)測站布點智能化選址領(lǐng)域。