高海燕,馬文娟,李唯欣,張 悅

(蘭州財經(jīng)大學(xué) 統(tǒng)計學(xué)院,甘肅 蘭州 730020)

隨著物聯(lián)網(wǎng)傳感、智能終端監(jiān)測技術(shù)的日益成熟,環(huán)境空氣污染物監(jiān)測站點采集的數(shù)據(jù)呈現(xiàn)“樣本多、采樣頻率高”等特點。然而,相關(guān)研究表明,與環(huán)境空氣質(zhì)量監(jiān)測數(shù)據(jù)的數(shù)量相比,數(shù)據(jù)的有效性和準確性更為重要[1]。2016 年1 月1日,在全國范圍內(nèi)實施的《環(huán)境空氣質(zhì)量標準》(GB 3095—2012)(簡稱“標準”)中對環(huán)境空氣質(zhì)量數(shù)據(jù)的有效性做出了明確的規(guī)定,如SO2、NO2、CO 等計算24 h 平均濃度時,要求每日至少有20 h 平均濃度值或采樣時間,O3計算8 h 平均濃度時,要求每8 h 至少有6 h 平均濃度值。但在監(jiān)測設(shè)備正常運行時,通信故障、傳輸期間數(shù)據(jù)包丟失等使得空氣質(zhì)量監(jiān)測數(shù)據(jù)存在大量缺失。一方面,數(shù)據(jù)的大規(guī)模缺失導(dǎo)致污染物數(shù)據(jù)不符合有效性規(guī)定,直接影響后續(xù)日均值、月均值等指標的計算[2];另一方面,完整的環(huán)境空氣質(zhì)量數(shù)據(jù)是分析污染物特征的基礎(chǔ),目前大多研究均基于中國環(huán)境監(jiān)測總站或當?shù)厣鷳B(tài)環(huán)境局發(fā)布的在線實時監(jiān)測數(shù)據(jù)[3-6],其中缺失值的存在會不同程度地增大統(tǒng)計分析的復(fù)雜性和難度,降低統(tǒng)計推斷的精度,最終導(dǎo)致統(tǒng)計分析結(jié)果存在偏誤。因此,如何科學(xué)、有效地處理環(huán)境空氣質(zhì)量監(jiān)測數(shù)據(jù)缺失問題具有重要的應(yīng)用價值。

目前,對于環(huán)境空氣質(zhì)量缺失數(shù)據(jù)的部分研究主要從體制機制的角度出發(fā),通過完善監(jiān)測制度規(guī)范及監(jiān)測環(huán)境,從根本上減少缺失值的出現(xiàn),從而提高監(jiān)測數(shù)據(jù)的質(zhì)量[7-8]。同時,另一部分研究者利用缺失插補方法填充環(huán)境空氣質(zhì)量缺失數(shù)據(jù),但當數(shù)據(jù)缺失類型有別于傳統(tǒng)類型且缺失比例較高,甚至出現(xiàn)大量“條狀、塊狀”缺失時,傳統(tǒng)插補方法失效。常規(guī)缺失數(shù)據(jù)處理策略可歸納為兩步法和一步法。當缺失比例較高時,兩步法的插補方法失效,直接刪除法[9]會造成數(shù)據(jù)信息的二次損失;“均值類”插補方法[10-11](如均值插補、中位數(shù)插補等)不能反映缺失數(shù)據(jù)的變異性;當數(shù)據(jù)高度集中時,熱卡插補法填充效果差,并在模擬數(shù)據(jù)的分布特征時缺乏準確性[12];涉及“距離度量”的插補方法也會導(dǎo)致處理效果不理想。如當缺失比例較大或缺失數(shù)據(jù)點大量連續(xù)時K-近鄰算法[13]插補效果較差;EM 插補要求數(shù)據(jù)分布服從正態(tài)分布,且當數(shù)據(jù)量大時計算過程較為繁瑣、算法收斂速度過于依賴初值的選擇[1];多重插補考慮了缺失數(shù)據(jù)的不確定性,插補效果雖好但計算量巨大,不適用于處理大規(guī)模缺失數(shù)據(jù)[15]。支持向量機的參數(shù)選擇較為困難,當訓(xùn)練樣本數(shù)目大時計算速度慢[16]。插補與分析并行的一步法中,貝葉斯模型[17]適用數(shù)據(jù)量較小的情況,當缺失率較高時,估計參數(shù)增加,使得估計方差增大、填充效果降低;神經(jīng)網(wǎng)絡(luò)[18]對于缺失數(shù)據(jù)敏感性較差,同時要求大樣本。此外,考慮到環(huán)境空氣質(zhì)量數(shù)據(jù)的時序特征及空間特征,大量時空數(shù)據(jù)插補方法被提出并應(yīng)用于處理缺失值[19-20]。上述傳統(tǒng)的插補方法大都沒有考慮空氣質(zhì)量數(shù)據(jù)相鄰時間、相鄰站點之間的相關(guān)特征信息以及變化軌跡的潛在模式,而函數(shù)型數(shù)據(jù)分析[21](Functional Data Analysis,FDA)是一種旨在克服這些挑戰(zhàn)的方法。

大規(guī)模空氣質(zhì)量數(shù)據(jù)的外在表現(xiàn)形式雖是離散、稀疏的片段點集,但內(nèi)在結(jié)構(gòu)卻呈現(xiàn)出連續(xù)、動態(tài)的函數(shù)曲線(或曲面)特征。因此,FDA 方法被廣泛應(yīng)用于環(huán)境空氣質(zhì)量數(shù)據(jù)的研究中[23-26],這些研究表明FDA 方法用于空氣質(zhì)量數(shù)據(jù)分析是有效可行的,且從函數(shù)的角度分析可以得到更有效的結(jié)果。同時,利用數(shù)據(jù)的函數(shù)性特征,可以從函數(shù)的視角對缺失數(shù)據(jù)進行處理[26]。故本研究從空氣質(zhì)量數(shù)據(jù)的“函數(shù)”特性出發(fā)插補其缺失值,利用一種函數(shù)型數(shù)據(jù)插補方法——概率函數(shù)聚類插補(Probabilistic Functional Clustering Imputation,PFCI)方法[27]處理北京市35 個監(jiān)測站點的空氣質(zhì)量數(shù)據(jù)。相較于常規(guī)插補方法,本研究所用插補方法效果更好,這為大規(guī)模稀疏空氣質(zhì)量數(shù)據(jù)補全提供了一種有效的思路與方法,從而保障后續(xù)統(tǒng)計分析的可靠性。

1 模型方法

盡管FDA 在空氣質(zhì)量研究中應(yīng)用廣泛,但學(xué)者們在處理缺失數(shù)據(jù)時仍采用傳統(tǒng)插補方法,沒有充分利用空氣質(zhì)量數(shù)據(jù)函數(shù)曲線所含的信息。同時,考慮到結(jié)合聚類過程可以進一步提高缺失數(shù)據(jù)的插補精度,采用PFCI 方法處理空氣質(zhì)量數(shù)據(jù)的缺失問題。該方法對已知數(shù)據(jù)進行聚類分析,在分析過程中預(yù)測缺失數(shù)據(jù),是一種聚類與插補并行完成的方法。設(shè)每個站點監(jiān)測到的空氣質(zhì)量數(shù)據(jù)變動軌跡為受測量誤差污染的隨機函數(shù)的實現(xiàn),并利用子空間投影技術(shù)[28]識別不完整函數(shù)數(shù)據(jù)的潛在模式。假設(shè)每個函數(shù)實現(xiàn)來自隨機過程的混合,其中每個子過程代表一個聚類子空間,該子空間由均值函數(shù)和基于函數(shù)型主成分分析(Functional Principal Component Analysis,FPCA)的協(xié)方差算子的特征函數(shù)所張成,PFCI 方法是非參數(shù)的、不需要任何分布假設(shè),能夠自適應(yīng)地擬合每個類別中的均值和協(xié)方差函數(shù)。因此,PFCI 方法能有效識別空氣質(zhì)量數(shù)據(jù)隨機變化模式的類別,并借助每類子空間中的信息插補缺失值,能提高插補精度。

其中,Γ(s,t)可以通過特征分解表示為Γ(s,t)=∑rλrφr(s)φr(t),λ1≥λ2≥…≥0 是滿足∑rλr＜∞的特征值,{φr}是對應(yīng)于特征值{λr}的特征函數(shù)。在不考慮潛在分組結(jié)構(gòu)的情況下,隨機軌跡X用FPCA 模型[29]表示為:

其中,隨機效應(yīng)ξr=(X-μ,φr)是零均值,方差為λr的不相關(guān)的函數(shù)型主成分得分。此外,根據(jù)式(3):

觀察在時間點Tij收集的被測量誤差污染的第i條隨機軌跡Xi(t)。其中,測量誤差εij是具有零均值、有限方差Var{εij}=σ2的不相關(guān)隨機誤差且獨立于ξir=(Xi-μ,φr)。

基于FPCA 和子空間投影技術(shù)的PFCI 方法,能并行實現(xiàn)空氣質(zhì)量數(shù)據(jù)的自適應(yīng)聚類及缺失值插補,具體執(zhí)行步驟如表1 所示:

表1 PFCI 方法過程

2 數(shù)據(jù)介紹與分析

2.1 數(shù)據(jù)介紹

本研究所用污染物濃度數(shù)據(jù)來自中國環(huán)境監(jiān)測總站(http://www.cnemc.cn/),選取北京市35個環(huán)境空氣污染物監(jiān)測站點2021 年1 月1 日—12 月31 日的6 種污染物(PM2.5、PM10、SO2、NO2、CO、O3)日均數(shù)據(jù)作為研究對象,共計365 d,缺失數(shù)據(jù)共計18 531 條,總體缺失率為24.17%。而不同空氣污染物的缺失率各不相同,PM10為61%,CO 為20%,O3為19%,NO2為17%,SO2為17%,PM2.5為11%,其中PM10缺失率最高,PM2.5缺失率最低。

2.2 監(jiān)測點數(shù)據(jù)缺失情況

如表2 所示,不同污染物35 個監(jiān)測站點的缺失情況存在較大差異,以PM10為例,缺失最嚴重的監(jiān)測站點(京東南區(qū)域點)缺失率高達78%,大規(guī)模缺失導(dǎo)致京東南區(qū)域點監(jiān)測到的數(shù)據(jù)有效性不足,若直接使用該數(shù)據(jù),可能會降低統(tǒng)計推斷的精度,最終導(dǎo)致統(tǒng)計分析結(jié)果偏誤。

表2 各監(jiān)測站點6 種污染物數(shù)據(jù)缺失統(tǒng)計 %

2.3 缺失機制分析

Rubin[30]研究了數(shù)據(jù)的不同缺失機制,并劃分為完全隨機缺失(MCAR)、隨機缺失(MAR)以及非隨機缺失(NMAR)三種類型。在實踐中,空氣質(zhì)量數(shù)據(jù)的缺失模式可能是MCAR 和MAR 模式的組合,由于無法根據(jù)數(shù)據(jù)對MAR、MCAR 與NMAR 進行區(qū)分,本研究將其缺失模式分為點缺失(Point missing,PM)、區(qū)間缺失(Interval missing,IM)和混合缺失(PM/IM)[27]。PM 集合中缺失點完全獨立、孤立且隨機分散,而IM 集合中包含隨機分布的未觀測區(qū)間,PM/IM 集合是PM 和IM 的混合。以缺失率最高的PM10為例,北京市各監(jiān)測站點PM10濃度數(shù)據(jù)的缺失特征如圖1 所示。圖1 中利用灰度的深淺程度近似表示PM10濃度觀測數(shù)據(jù)取值的大小,顏色越淺表示值越小,黑色方塊表示數(shù)據(jù)缺失點,每一行為各監(jiān)測站點365 d 的缺失情況,每一列為一天中35 個站點的缺失情況。從圖1 可直觀看出,黑色方塊集中表現(xiàn)為“條狀、塊狀”分布,故PM10濃度數(shù)據(jù)存在大規(guī)模“條狀、塊狀”缺失、缺失比例較高的特點,為混合PM/IM 缺失模式。PM 和IM 頻率分布情況如表3 所示。每個數(shù)據(jù)集中大部分缺失模式為PM,其中缺失比率最高達該數(shù)據(jù)集總?cè)笔У?8.58%。缺失間隔長度從2～17 d 不等。因此,傳統(tǒng)插補方法處理上述缺失數(shù)據(jù)集會存在局限性。

圖1 北京市35 個監(jiān)測站點PM10 濃度數(shù)據(jù)缺失情況

表3 點缺失(零長度)和不同長度缺失區(qū)間的頻數(shù)分布

3 模擬實驗

3.1 實驗設(shè)計

通過模擬缺失值插補實驗,驗證PFCI 方法插補大規(guī)模稀疏空氣質(zhì)量函數(shù)型數(shù)據(jù)的有效性,并在不同缺失比例的設(shè)定下,與三種典型插補方法進行對比,評價其插補性能。因此,本研究通過以下三個步驟完成方案設(shè)計:

同時，材料的力學(xué)性能也是重中之中，比如屈服極限、彈性模量、彎曲強度、表面硬度等都應(yīng)該在合適的條件內(nèi)。如下表所示：

3.1.1 數(shù)據(jù)準備

選取空氣質(zhì)量在線監(jiān)測平臺(http://www.aqistudy.cn/historydata)發(fā)布的2014—2021 年京津冀地區(qū)14 個城市污染物的完整月度數(shù)據(jù)?？紤]到6 種污染物之間數(shù)據(jù)特征以及單位的不同,選取PM2.5、CO 濃度數(shù)據(jù)作為實驗數(shù)據(jù)集。

3.1.2 人工生成帶有PM 和IM 缺失的條目

具體工作如下:

(1) 點缺失(PM):以逐點方式對缺失的條目進行統(tǒng)一隨機采樣。

(2) 區(qū)間缺失(IM):生成不同大小的缺失區(qū)間,這些區(qū)間在每條曲線軌跡中的起始位置以及長度遵循均勻分布。

為了模擬第2.3 節(jié)中介紹的缺失模式,本研究采用PM 以及混合PM/IM 模式,同時為每種模式生成缺失點,其中是曲線軌跡數(shù),是時間點數(shù),是缺失率并分別設(shè)定為5%,10%,20%,30%,40%,50%和60%。

3.1.3 確定PFCI 方法的對比方法以及評估指標

(1) 對比方法:選取常用的均值填充方法,以及兩種以矩陣填充技術(shù)為核心思想的典型函數(shù)型數(shù)據(jù)插補方法[31]:SFI(Soft Functional Impute,SFI)、HFI(Hard Functional Impute,HFI)。

(2) 評估指標:利用均方根誤差(RMSE)、平均絕對百分比誤差(MAPE)評估每種方法的插補性能。RMSE、MAPE 的值越小,表明插補值與真實值之間的誤差越小。兩個評估指標的計算公式分別為:

式中:n0——包含人工缺失條目的曲線軌跡數(shù);mi——第i條軌跡的缺失點數(shù);——第i條曲線軌跡第j個生成缺失的時間點。

3.2 實驗結(jié)果

在上述實驗設(shè)計的基礎(chǔ)上,針對不同的缺失率,分別采用均值填充、SFI、HFI 和PFCI 方法修復(fù)兩個數(shù)據(jù)集并進行插補性能的比較。實驗結(jié)果如圖2、圖3 所示:

圖3 不同缺失率下CO 濃度數(shù)據(jù)的插補誤差RMSE 和MAPE 對比(10 次重復(fù)模擬結(jié)果均值)

由圖2 可以發(fā)現(xiàn),在插補PM2.5濃度缺失數(shù)據(jù)時,常規(guī)插補方法中的均值填充在兩種缺失模式以及不同缺失率的情況下,插補誤差RMSE 均大于函數(shù)型數(shù)據(jù)缺失值插補方法,并且當缺失率增大到40%,50%,60%時,RMSE 和MAPE 顯著增大。此外,相較于SFI 和HFI,不同缺失率下PFCI 方法插補PM2.5濃度數(shù)據(jù)的RMSE、MAPE更小,插補性能更優(yōu)。當缺失率達到50%,60%時,SFI 和HFI 的插補誤差顯著增大,而PFCI 方法的插補誤差變化較小、趨于穩(wěn)定。

觀察圖3 可以發(fā)現(xiàn),與插補PM2.5濃度缺失數(shù)據(jù)類似,四種方法在插補CO 濃度缺失數(shù)據(jù)時,均值填充的插補誤差均大于SFI、HFI 和PFCI,其中,當缺失率大于等于10%時,均值填充的MAPE達到100%,說明插補值與真實值之間的誤差較大。此外,當缺失率為10%,20%,30%時,SFI、HFI 和PFCI 三種方法的插補誤差無明顯差異,但當缺失率為40%,50%以及60%時,SFI 和HFI 方法的插補誤差顯著增大,相反,PFCI 方法的插補誤差卻無明顯變化,趨于平穩(wěn)。綜上所述,相較于其他三種插補方法,PFCI 方法當缺失率較高時具有明顯的插補優(yōu)勢。

圖4 為在不同缺失率下均值填充、SFI、HFI和PFCI 方法插補PM2.5、CO 濃度缺失數(shù)據(jù)10 次模擬實驗運行時間對比圖。由圖4 可知,PFCI 方法在插補兩個數(shù)據(jù)集時運行時間都高于均值填充、SFI 以及HFI,這是PFCI 方法在執(zhí)行的過程中加入Logit 模型迭代分類所導(dǎo)致的,而迭代分類步驟是PFCI 方法的重要子過程,收斂的聚類結(jié)果提高了PFCI 方法插補缺失值的有效性。

圖4 4 種方法插補PM2.5 濃度缺失數(shù)據(jù)和CO 濃度缺失數(shù)據(jù)的運行時間對比

4 PFCI 方法在稀疏空氣質(zhì)量數(shù)據(jù)中的實證應(yīng)用

在上述模擬實驗驗證PFCI 方法插補環(huán)境空氣質(zhì)量數(shù)據(jù)有效性的基礎(chǔ)上,針對大規(guī)模稀疏空氣質(zhì)量數(shù)據(jù)開展插補應(yīng)用。利用PFCI 方法對2021 年北京市35 個監(jiān)測站點的6 種空氣污染物濃度數(shù)據(jù)進行缺失插補。北京市空氣質(zhì)量35 個監(jiān)測站點的站點名稱、站點類型等基本信息如表4 所示,依據(jù)位置信息可將站點分為5 類(為描述方便起見,分別定義為第1 類、第2 類、第3 類、第4 類、第5 類)。

表4 北京市空氣質(zhì)量監(jiān)測站點基本信息

為檢驗PFCI 方法聚類的準確性,以缺失率最高的PM10為例,將PFCI 方法中Logit 模型的收斂結(jié)果作為35 個監(jiān)測站點的最終聚類結(jié)果,聚類結(jié)果如表5 所示。觀察表5 可知,PM10缺失率高達61%,導(dǎo)致PFCI 方法的聚類結(jié)果與實際站點的布設(shè)情況存在一些差異,但整體上聚類結(jié)果較好,其中第3 類的聚類準確率達到75%。借助最終聚類結(jié)果,插補PM10濃度數(shù)據(jù)缺失值,以京東南區(qū)域點(PM10缺失率最高的監(jiān)測站點)為例,6 種污染物插補結(jié)果如圖5 所示,其中觀測值、預(yù)測軌跡和插補的缺失值分別用黑色的點、灰色曲線以及黑色小圓圈標記。

圖5 京東南區(qū)域點6 種污染物濃度數(shù)據(jù)插補效果對比

續(xù)圖5

表5 北京市空氣質(zhì)量監(jiān)測站點聚類結(jié)果

圖5 中京東南區(qū)域點各項污染物的預(yù)測軌跡由預(yù)測函數(shù)式(6)生成,

圖6 京東南區(qū)域點6 種污染物真實觀測值關(guān)于預(yù)測值的散點圖

上述分析結(jié)果表明,PFCI 方法作為“聚類+插補”并行的一步法,不僅能夠較好地識別監(jiān)測站點的空間布局,在一定程度上降低缺失值對聚類分析的影響。同時,PFCI 方法在處理缺失值時,大規(guī)模點狀、區(qū)間缺失或者“條狀、塊狀”缺失以及高缺失率等特征對該方法的影響較小,且相較于其他插補方法,PFCI 方法對高缺失率的敏感度較低。

5 結(jié)論

智能終端監(jiān)測技術(shù)使得環(huán)境空氣質(zhì)量在線監(jiān)測系統(tǒng)覆蓋全國,系統(tǒng)在運行期間受眾多因素的影響,導(dǎo)致出現(xiàn)數(shù)據(jù)缺失問題,而數(shù)據(jù)的缺失通常會對后續(xù)統(tǒng)計分析產(chǎn)生顯著影響。本研究首先選取2014—2021 年京津冀地區(qū)14 個城市PM2.5、CO的完整數(shù)據(jù)進行插補方法性能評估實驗,實驗發(fā)現(xiàn),當出現(xiàn)大規(guī)?！皸l狀、塊狀”缺失、缺失率較高時,PFCI 方法的RMSE、MAPE 均小于其他方法,并且插補誤差變化較小。基于此,利用PFCI 方法對2021 年北京市35 個監(jiān)測站點的6 種空氣污染物缺失數(shù)據(jù)進行插補應(yīng)用,結(jié)果表明,PFCI 方法不僅能準確識別空氣污染物變化模式的類別,還利用每一類中的信息處理缺失值,插補性能較好。