大規(guī)?？諝赓|(zhì)量監(jiān)測數(shù)據(jù)缺失處理方法實(shí)證研究

張 波,宋國君

大規(guī)?？諝赓|(zhì)量監(jiān)測數(shù)據(jù)缺失處理方法實(shí)證研究

張 波,宋國君*

(中國人民大學(xué)環(huán)境學(xué)院,北京 100872)

基于2016年1月至2021年7月的全國1654個(gè)國控監(jiān)測點(diǎn)小時(shí)級的6種污染物空氣質(zhì)量監(jiān)測數(shù)據(jù),研究缺失值處理方法、效果及其影響.模擬實(shí)驗(yàn)表明交替最小二乘下的低秩矩陣插補(bǔ)算法相比于其他缺失值處理方法擁有更小的均方根誤差、平均百分比誤差,更高的相關(guān)系數(shù)和更快的運(yùn)算速度,在大規(guī)模數(shù)據(jù)集上性能更優(yōu).實(shí)證分析表明應(yīng)用文本方法得到的插補(bǔ)值是有效且合理的,缺失值插補(bǔ)前后污染物濃度評估值會有±10%以內(nèi)的變化,插補(bǔ)后的數(shù)據(jù)集更加準(zhǔn)確和完備.本文建議在基于空氣質(zhì)量監(jiān)測數(shù)據(jù)研究時(shí)應(yīng)先采用本文中的缺失數(shù)據(jù)處理方法,對監(jiān)測數(shù)據(jù)中存在的缺失數(shù)據(jù)進(jìn)行插補(bǔ),提高研究所使用監(jiān)測數(shù)據(jù)的完整性,保證相關(guān)計(jì)算結(jié)果的準(zhǔn)確性和有效性.

監(jiān)測數(shù)據(jù)；大數(shù)據(jù)處理；缺失值處理；矩陣填充；實(shí)證研究

我國2016年實(shí)施的《環(huán)境空氣質(zhì)量標(biāo)準(zhǔn) (GB3095)》(簡稱“標(biāo)準(zhǔn)”)中對污染物監(jiān)測數(shù)據(jù)的有效性做出明確規(guī)定,一天24h中有至少20h以上數(shù)據(jù)時(shí),計(jì)算SO2、NO2、CO、PM2.5和PM10的日均濃度值才是有效日均數(shù)據(jù),O3監(jiān)測數(shù)據(jù)則要求每連續(xù)8h中至少有6h數(shù)據(jù)時(shí)才為有效數(shù)據(jù).但在監(jiān)測設(shè)備的日常運(yùn)行中,由于某些不可預(yù)期因素(如服務(wù)器宕機(jī)、網(wǎng)絡(luò)中斷、設(shè)備故障等)導(dǎo)致監(jiān)測數(shù)據(jù)缺失.目前國家公開發(fā)布的空氣質(zhì)量監(jiān)測數(shù)據(jù)以小時(shí)數(shù)據(jù)為最小顆粒度,小時(shí)數(shù)據(jù)缺失會對后續(xù)日均值、月均值、年度均值、超標(biāo)率、二級標(biāo)準(zhǔn)天數(shù)等指標(biāo)的計(jì)算造成直接影響.2018年生態(tài)環(huán)境部印發(fā)的《城市環(huán)境空氣質(zhì)量排名技術(shù)規(guī)定》[1](簡稱“規(guī)定”)中指出SO2、NO2、PM2.5、PM10的評價(jià)濃度為日均濃度,O3的評價(jià)濃度為日最大8h平均值的90%分位數(shù),CO的評價(jià)濃度為日均濃度的95%分位數(shù),不符合有效性規(guī)定的數(shù)據(jù)不能參與計(jì)算,為空氣質(zhì)量監(jiān)測數(shù)據(jù)的計(jì)算提供了標(biāo)準(zhǔn).大規(guī)?？諝赓|(zhì)量監(jiān)測數(shù)據(jù)是進(jìn)行污染物時(shí)空分布規(guī)律研究的基礎(chǔ)數(shù)據(jù),目前大部分研究均基于中國環(huán)境監(jiān)測總站發(fā)布的連續(xù)在線監(jiān)測數(shù)據(jù)[2-8],數(shù)據(jù)的缺失和不同的缺失處理方法會對研究結(jié)果產(chǎn)生一定影響[9],但空氣質(zhì)量監(jiān)測數(shù)據(jù)中數(shù)據(jù)缺失的影響究竟有多大,使用何種缺失處理方法更加可靠,這方面的研究還非常有限,因而研究如何更加客觀、科學(xué)的處理大規(guī)模環(huán)境監(jiān)測數(shù)據(jù)中缺失問題具有重要的現(xiàn)實(shí)應(yīng)用價(jià)值.

目前國內(nèi)有關(guān)環(huán)境監(jiān)測數(shù)據(jù)缺失問題的研究主要集中在制度建設(shè)和管理規(guī)范方面,如建立外場監(jiān)測和實(shí)驗(yàn)室分析的兩級質(zhì)量保證體系[10],通過開展專項(xiàng)質(zhì)控工作提升監(jiān)測數(shù)據(jù)質(zhì)量[11],加強(qiáng)對生態(tài)環(huán)境監(jiān)測機(jī)構(gòu)的監(jiān)管提升數(shù)據(jù)質(zhì)量[12],通過建立從業(yè)人員監(jiān)管制度確保數(shù)據(jù)質(zhì)量[13],以及對國外經(jīng)驗(yàn)的總結(jié)和借鑒[14],這些研究從體制機(jī)制的角度研究如何提升監(jiān)測數(shù)據(jù)質(zhì)量,但缺少對具體現(xiàn)實(shí)具體問題的探討,尤其是存在大規(guī)模數(shù)據(jù)缺失情況下該如何處理,用什么方法處理,處理后的影響如何等問題還鮮有研究.國際上環(huán)境監(jiān)測數(shù)據(jù)的缺失問題一直以來都是研究者非常關(guān)注的問題之一,很多缺失值處理方法也被應(yīng)用于空氣質(zhì)量監(jiān)測數(shù)據(jù)處理和研究中,如最近鄰算法[15]、EM算法[16]、線性回歸[17]、簡單插補(bǔ)[18-20]、多重插補(bǔ)[18,21]、最小核范數(shù)[22]和低秩矩陣插補(bǔ)[23-24]等,多方法對比研究[21,24]發(fā)現(xiàn)多重插補(bǔ)方法和低秩矩陣方法在數(shù)據(jù)缺失較大的情況下更有優(yōu)勢,能夠產(chǎn)生偏差更小的插補(bǔ)值來替代缺失值.以往研究對象面向的是少數(shù)監(jiān)測站點(diǎn)和較短數(shù)據(jù)序列,現(xiàn)實(shí)中由于監(jiān)測手段快速改進(jìn),空氣質(zhì)量監(jiān)測數(shù)據(jù)已經(jīng)出現(xiàn)爆發(fā)式增加,監(jiān)測點(diǎn)數(shù)量和數(shù)量量都達(dá)到很大規(guī)模,如何處理海量環(huán)境數(shù)據(jù)中的缺失問題,已經(jīng)成為越來越需要密切關(guān)注的問題.本文基于全國空氣質(zhì)量國控點(diǎn)的小時(shí)級監(jiān)測數(shù)據(jù),利用低秩矩陣插補(bǔ)方法處理海量監(jiān)測數(shù)據(jù)中大量缺失的問題,該方法在保證插補(bǔ)精度的同時(shí)能夠大幅度降低大規(guī)模數(shù)據(jù)處理所需時(shí)間,得到的插補(bǔ)值服具有較高有效性.

1 模型方法

在缺失值處理方法中,簡單插補(bǔ)和最近鄰算法方法最簡單,但處理效果較差,EM算法只適用于分布為多元正態(tài)分布的數(shù)據(jù),線性回歸、多重插補(bǔ)雖然效果較好但計(jì)算量巨大,無法適用于大規(guī)模缺失處理.本文使用交替最小二乘下的低秩矩陣插補(bǔ)(Low Rank SVD via Alternating Least Square, softImpute- ALS)[25]算法,其是在低秩矩陣插補(bǔ)方法[26-27]基礎(chǔ)上為解決大規(guī)模數(shù)據(jù)集處理而提出的.低秩矩陣插補(bǔ)方法最初在Netflix競賽中被提出,用于插補(bǔ)電影評分矩陣中的缺失值,實(shí)現(xiàn)電影推薦的目的,因?yàn)槠湓诖笠?guī)模數(shù)據(jù)集上的優(yōu)異性能而備受關(guān)注.

基本思想是在原始存在缺失值的矩陣基礎(chǔ)上進(jìn)行重建,即對矩陣中的缺失值進(jìn)行填充,并保證填充后的完整矩陣與原始矩陣的秩相同.給定數(shù)據(jù)集,且中存在缺失值,由已觀測值的下標(biāo)構(gòu)成集合W,那么softImpute-ALS算法的目標(biāo)是求解以下優(yōu)化問題:

步驟1,初始化矩陣=,其中′r為隨機(jī)生成且列正交,′r為單位矩陣.初始化矩陣=,其中′r=0.

步驟2,固定矩陣,通過優(yōu)化式(2)來更新矩陣,

步驟3,固定矩陣,在第二步中交換矩陣和,使用同樣的方法更新矩陣和矩陣.

步驟4,重復(fù)步驟2和步驟3,直到達(dá)到收斂條件.

步驟5,計(jì)算=′,并對矩陣進(jìn)行SVD分解,有=R.最終輸出結(jié)果、?和,=diag[(1-)+,…,(-)+],以及填充后的矩陣*.

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

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

本文研究數(shù)據(jù)來自中國環(huán)境監(jiān)測總站實(shí)時(shí)發(fā)布的全國367個(gè)城市,1654個(gè)國控監(jiān)測點(diǎn)自2016年1月1日到2021年7月31日的6種空氣污染物小時(shí)濃度監(jiān)測數(shù)據(jù),共計(jì)48912h,理論上如果每個(gè)監(jiān)測點(diǎn)的每種污染物在每個(gè)小時(shí)都有監(jiān)測數(shù)據(jù)的話,那么監(jiān)測數(shù)據(jù)總量為4.2億條,但由于各種原因?qū)е聰?shù)據(jù)出現(xiàn)缺失(即空值),缺失數(shù)據(jù)總計(jì)814萬條.本文將缺失數(shù)據(jù)總量除以理論數(shù)據(jù)總量定義為數(shù)據(jù)缺失率,那么國控點(diǎn)監(jiān)測數(shù)據(jù)的總體缺失率為19.3‰.數(shù)據(jù)中不同污染物缺失率并不相同,SO2缺失率為9.6‰,NO2為10.9‰,CO為13.1‰,O3為14.4‰,PM2.5為15.3‰,PM10為52.4‰,PM10缺失率最高,而SO2缺失率最低.

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

不同監(jiān)測點(diǎn)的缺失情況存在很大差異(表1),以PM2.5為例,缺失最嚴(yán)重的監(jiān)測點(diǎn)的缺失率達(dá)到331.7‰,即三分之一的數(shù)據(jù)因各種原因缺失,致使監(jiān)測點(diǎn)數(shù)據(jù)有效性大打折扣,在使用該數(shù)據(jù)進(jìn)行分析時(shí),可能會導(dǎo)致較大誤差.因而有關(guān)部門在使用空氣質(zhì)量監(jiān)測數(shù)據(jù)時(shí),會對監(jiān)測數(shù)據(jù)進(jìn)行校驗(yàn),但當(dāng)前大部分研究中所使用的公開數(shù)據(jù)集往往是未經(jīng)過校驗(yàn)處理的,其中就存在大量缺失值,缺失值對相關(guān)計(jì)算結(jié)果的影信息響究竟有多大,還沒有被深入討論.

表1 不同污染物的監(jiān)測點(diǎn)數(shù)據(jù)缺失統(tǒng)計(jì)(‰)

2.3 缺失機(jī)制分析

數(shù)據(jù)缺失的內(nèi)在機(jī)理可以劃分3種類型[28],分別是完全隨機(jī)缺失(MCAR)、隨機(jī)缺失(MAR)和非隨機(jī)缺失(MNAR).其中MAR相比于MCAR更加常見和符合現(xiàn)實(shí),也是缺失值處理方法最主要的研究對象.MAR指僅在某個(gè)特定的組內(nèi)缺失值是隨機(jī)產(chǎn)生的,而不同組之間不一定是隨機(jī)的.對應(yīng)到空氣質(zhì)量數(shù)據(jù)的缺失機(jī)制,由于不同污染物缺失率不同,不同監(jiān)測點(diǎn)的缺失率也有很大差異,顯然污染物缺失并不是完全隨機(jī)的,而僅在特定監(jiān)測點(diǎn)和污染物條件下是可以視作是隨機(jī)缺失的.此外,空氣質(zhì)量是連續(xù)在線監(jiān)測模式,如果因設(shè)備故障而導(dǎo)致數(shù)據(jù)缺失的話,數(shù)據(jù)往往會連續(xù)缺失一定時(shí)段直到設(shè)備故障排除.表2顯示,三分之二左右的缺失間隔僅為1,說明大部分情況下缺失值是偶然出現(xiàn).但也有三分之一左右的缺失間隔大于1,缺失在10h以上的情況占比達(dá)到3%以上,這很可能是出現(xiàn)明顯的設(shè)備故障,需要一定的維修時(shí)間.其中PM2.5最大連續(xù)缺失200h,PM10最大連續(xù)缺失196h,SO2最大連續(xù)缺失192h,NO2最大連續(xù)缺失197h,O3最大連續(xù)缺失189h,CO最大連續(xù)缺失191h.研究表明隨著缺失間隔的增大,使用簡單缺失值處理方法(如均值填充等)的有效性不斷下降[29-30],更加需要有效方法來處理大規(guī)模缺失的情況,而softImpute-ALS方法在嚴(yán)重缺失的情況也能有很好的插補(bǔ)性能[25],因而適用于處理空氣質(zhì)量監(jiān)測中的缺失值.

表2 不同污染物缺失間隔占比(%)

3 模擬實(shí)驗(yàn)

3.1 實(shí)驗(yàn)設(shè)計(jì)

為考察softImpute-ALS方法在大規(guī)模缺失數(shù)據(jù)集上的優(yōu)劣,需要構(gòu)造相應(yīng)數(shù)據(jù)集,并計(jì)算各種缺失值處理方法在不同評估指標(biāo)上的表現(xiàn).因此,本文通過以下模擬方案生成算法評估數(shù)據(jù)集.方案包括四個(gè)步驟:(1)抽取可能產(chǎn)生缺失的污染物,抽取概率為不同污染物的缺失值數(shù)量占全部缺失值數(shù)量的比例;(2)抽取可能產(chǎn)生缺失的監(jiān)測點(diǎn),抽取概率為污染物下,監(jiān)測點(diǎn)缺失值數(shù)量占全部該污染物全部缺失值數(shù)量的比例;(3)隨機(jī)抽取可能產(chǎn)生缺失的時(shí)間點(diǎn),即該時(shí)刻開始出現(xiàn)缺失值;(4)抽取產(chǎn)生缺失的間隔長度,抽取概率為表2中缺失間隔對應(yīng)的占比;(5)記錄缺失值的下標(biāo)(,,,),并將原始數(shù)據(jù)中相同下標(biāo)的數(shù)據(jù)標(biāo)記為缺失值.通過以上四個(gè)步驟,在模擬真實(shí)缺失機(jī)制前提下,本文生成得到兩份數(shù)據(jù)集,一份為原始數(shù)據(jù)集,一份是標(biāo)記了缺失值的數(shù)據(jù)集.接下來,本文對數(shù)據(jù)集運(yùn)用不同缺失值方法,通過考察其還原數(shù)據(jù)集的程度來評估方法的優(yōu)劣.

3.2 評估指標(biāo)

使用的評估指標(biāo)有均方根誤差(RMSE)、平均絕對百分比誤差(MAPE)和皮爾遜相關(guān)系數(shù)().RMSE指標(biāo)用于測度缺失值和插補(bǔ)值之間絕對偏差的大小,RMSE越小說明插補(bǔ)值約接近原始值.MAPE指標(biāo)用于測度缺失值和插補(bǔ)值之間相對誤差的大小,MAPE可以比較消除量綱后偏差的大小.指標(biāo)用于測度變量間線性相關(guān)關(guān)系,指標(biāo)值介于-1到1之間,越大說明缺失值與插補(bǔ)值之間相關(guān)性越高,二者擬合效果越好.不妨令模擬數(shù)據(jù)中觀測值下標(biāo)集合為W,則缺失值下標(biāo)構(gòu)成的集合是W的補(bǔ)集記為Wc,為原始的模擬數(shù)據(jù),為缺失值填充后的數(shù)據(jù),為缺失值數(shù)量,則評估指標(biāo)的具體計(jì)算公式為:

3.3 實(shí)驗(yàn)結(jié)果

選擇簡單填充、線性插值、向前填充、最近鄰插補(bǔ)、EM插補(bǔ)、多重插補(bǔ)方法作為對比,比較這些方法和softImpute-ALS在相同數(shù)據(jù)集上,模擬真實(shí)缺失機(jī)制下在RMSE、MAPE、和運(yùn)算時(shí)間(time)上的表現(xiàn).重復(fù)實(shí)驗(yàn)10次,最終結(jié)果為10次實(shí)驗(yàn)的平均值,實(shí)驗(yàn)使用平臺配置為2.2GHz的CPU和16G DDR3內(nèi)存,實(shí)驗(yàn)結(jié)果如表3所示(1)softImpute-ALS插補(bǔ)后的值在RMSE、MAPE兩個(gè)指標(biāo)上均表現(xiàn)最小,說明softImpute-ALS能夠更加準(zhǔn)確的插補(bǔ)缺失值,同時(shí)softImpute-ALS的插補(bǔ)值與真實(shí)值有更高的相關(guān)性,其他方法中只有最近鄰填充能夠與softImpute-ALS方法在準(zhǔn)確度方面更為接近.(2)時(shí)間消耗方面,softImpute-ALS對不同污染染處理時(shí)間介于160～180s之間,雖然相比于簡單填充、線性插值和向前填充耗時(shí)顯著增加,但相比最近鄰方法無疑有巨大的優(yōu)勢,最近鄰方法處理每種污染的用時(shí)在3h左右,是softImpute-ALS的70倍以上.實(shí)驗(yàn)表明softImpute-ALS方法相比于其他方法有更好的準(zhǔn)確度,并且處理時(shí)間可控,因而softImpute-ALS方法在處理大規(guī)模缺失數(shù)據(jù)方面具有顯著的優(yōu)勢.

表3 不同缺失值插補(bǔ)方法比較

4 實(shí)證分析

4.1 數(shù)據(jù)有效性分析

“規(guī)定”中指出“SO2、NO2、PM10、PM2.5的評價(jià)濃度為評價(jià)時(shí)段內(nèi)日均濃度的平均值,O3的評價(jià)濃度為評價(jià)時(shí)段內(nèi)日最大8h平均值的第90百分?jǐn)?shù),CO的評價(jià)濃度為評價(jià)時(shí)段內(nèi)日均濃度的第95百分位數(shù)”的計(jì)算方法,同時(shí)指出城市的日均濃度首先要“計(jì)算各監(jiān)測點(diǎn)位日均濃度,然后計(jì)算算數(shù)平均值得到城市日均濃度,再由此計(jì)算你評價(jià)時(shí)段內(nèi)城市均值或特定百分位數(shù)”.各監(jiān)測點(diǎn)位的日均濃度由該日24h監(jiān)測值求平均而得,在“標(biāo)準(zhǔn)”中規(guī)定每日至少有20h的平均濃度值”才能計(jì)算24h均值,否則該日數(shù)據(jù)無效.

依據(jù)以上標(biāo)準(zhǔn),統(tǒng)計(jì)每個(gè)監(jiān)測點(diǎn)每日濃度并判斷其有效性,然后計(jì)算每個(gè)城市平均濃度及其有效性.研究發(fā)現(xiàn),從2016年1月1日至2021年7月31日共計(jì)2038d,監(jiān)測點(diǎn)有效天數(shù)最大是O-3的1892d,意味著該監(jiān)測點(diǎn)有多達(dá)146d未滿足有效性要求.同時(shí),整體上有超過50%左右的監(jiān)測點(diǎn)數(shù)據(jù)有效時(shí)間介于1655～1840d之間,無效時(shí)間在198～383d之間(表4),可見由于存在大量缺失值,導(dǎo)致數(shù)據(jù)有效時(shí)間顯著減少,每日污染濃度的分析結(jié)果可能存在較大偏差.“規(guī)定”中并未指明在城市某個(gè)或某幾個(gè)監(jiān)測點(diǎn)當(dāng)日數(shù)據(jù)無效情況下,如何計(jì)算該城市當(dāng)日污染物平均濃度,而只是籠統(tǒng)的給出“當(dāng)任何一項(xiàng)污染物不滿足上述有效性規(guī)定且任何一項(xiàng)污染物濃度超過二級標(biāo)準(zhǔn)限值時(shí),以城市當(dāng)日污染物濃度最高點(diǎn)位的數(shù)據(jù)”的規(guī)定,這樣以“最大值”替代的處理方式對于研究污染物時(shí)空分布顯然并不適用,更加合理的方式是利用模型算法對缺失值進(jìn)行準(zhǔn)確的估計(jì).

表4 全部國控監(jiān)測點(diǎn)不同污染物監(jiān)測數(shù)據(jù)有效時(shí)間統(tǒng)計(jì)

4.2 缺失值影響分析

由于“規(guī)定”中沒有明確給出監(jiān)測點(diǎn)當(dāng)日數(shù)據(jù)無效的處理方法,當(dāng)監(jiān)測點(diǎn)當(dāng)日數(shù)據(jù)無效時(shí),分2種情況分別討論.第1種情況,如果一個(gè)城市當(dāng)日有一個(gè)或多個(gè)監(jiān)測點(diǎn)數(shù)據(jù)無效,則當(dāng)日該城市數(shù)據(jù)也無效;第2種情況,如果監(jiān)測點(diǎn)當(dāng)日數(shù)據(jù)無效,則取該監(jiān)測點(diǎn)當(dāng)日濃度最高值作為當(dāng)日有效數(shù)據(jù).在這2種情況下,缺失值均對后續(xù)當(dāng)日城市污染物濃度的計(jì)算均會產(chǎn)生重大影響.

在第1種情況下,以生態(tài)環(huán)境部發(fā)布的《2021年8月全國城市空氣質(zhì)量報(bào)告》[31]中空氣質(zhì)量相對較差的20個(gè)城市為例,分析這些城市由于數(shù)據(jù)缺失導(dǎo)致的無效數(shù)據(jù)問題.表5顯示,在過去5年多中,這20個(gè)城市各種污染物平均每年的無效時(shí)間約為60d,即有將近兩個(gè)月的時(shí)間按照“規(guī)定”是無效數(shù)據(jù),以濟(jì)南市PM10為例,無效時(shí)間更是達(dá)到881d,年均160d,無效時(shí)間的比例遠(yuǎn)超按小時(shí)統(tǒng)計(jì)的1%～5%的數(shù)據(jù)缺失率,其主要原因是每個(gè)監(jiān)測點(diǎn)每日須滿足大于等于20個(gè)有效監(jiān)測數(shù)據(jù).即使按照規(guī)定中日最大濃度替代日均值濃度的方法,也意味著這些城市每年有將近2個(gè)月的時(shí)間使用日最高濃度進(jìn)行排名計(jì)算.

表5 空氣質(zhì)量排名后20城市的數(shù)據(jù)無效時(shí)間統(tǒng)計(jì)

在第2種情況下,污染物濃度在一天內(nèi)存在很大的波動(dòng),造成最大濃度與平均濃度之間有很大差異,最大濃度無法反映真實(shí)情況.以濟(jì)南市機(jī)床二廠監(jiān)測點(diǎn)2021年2月3日為例,當(dāng)日有效監(jiān)測數(shù)據(jù)恰為19個(gè),小于20個(gè)有效監(jiān)測數(shù)據(jù)的要求,缺失14:00、16:00、17:00、18:00和19:00這5個(gè)小時(shí)監(jiān)測數(shù)據(jù).圖1顯示一天內(nèi)不同時(shí)點(diǎn)污染物濃度存在很大差異,該監(jiān)測點(diǎn)PM2.5最大濃度123mg/m3,19個(gè)小時(shí)平均值為81mg/m3,相差52%,PM10最大濃度214mg/m3,19小時(shí)平均值153mg/m3,相差40%,SO2最大濃度55mg/m3,19小時(shí)平均值20mg/m3,相差172%,NO2最大濃度69mg/m3,19小時(shí)平均值40mg/m3,相差75%, O3最大濃度87mg/m3,19小時(shí)平均值59mg/m3,相差47%,CO最大濃度1.3mg/m3,19小時(shí)平均濃度0.9mg/m3,相差41%.“規(guī)定”中的“最大濃度替代”與實(shí)際情況存在很大差異,不能反應(yīng)污染物的真實(shí)濃度,僅僅是在計(jì)算層面對缺失數(shù)據(jù)進(jìn)行懲罰性處理.當(dāng)存在某些不可抗拒因素導(dǎo)致數(shù)據(jù)缺失時(shí),這樣的處理方法并不合理,而應(yīng)當(dāng)采取更科學(xué)和可靠的插補(bǔ)方法,最大程度近似或還原真實(shí)值.

4.3 插補(bǔ)值有效性分析

應(yīng)用softImpute-ALS方法對全國國控點(diǎn)空氣質(zhì)量數(shù)據(jù)進(jìn)行缺失值的插補(bǔ)處理,未驗(yàn)證插補(bǔ)值的有效性,為與上文保持一致,選取濟(jì)南市機(jī)床二廠和科干所兩個(gè)監(jiān)測點(diǎn)來驗(yàn)證插補(bǔ)值的有效性.具體方法是分別計(jì)算兩個(gè)監(jiān)測點(diǎn)每種污染物在剔除缺失值后的數(shù)據(jù)、插補(bǔ)值數(shù)據(jù)、插補(bǔ)缺失值后的數(shù)據(jù)這3個(gè)數(shù)據(jù)集上的相關(guān)系數(shù),并且繪制散點(diǎn)圖(圖2)來分析插補(bǔ)值的有效性與合理性,即如果2個(gè)監(jiān)測點(diǎn)在3個(gè)數(shù)據(jù)集的相關(guān)系數(shù)比較接近,則說明插補(bǔ)值也符合原有數(shù)據(jù)規(guī)律,插補(bǔ)方法是有效的.

由圖2可見,插補(bǔ)值與剔除缺失的數(shù)據(jù)在污染物相關(guān)關(guān)系上是一致的,插補(bǔ)前兩個(gè)監(jiān)測點(diǎn)PM2.5的相關(guān)系數(shù)是0.93,插補(bǔ)后的相關(guān)系數(shù)是0.94,缺失值處理方法得到的插補(bǔ)值也服從監(jiān)測點(diǎn)間相關(guān)的規(guī)律.此外計(jì)算濟(jì)南市其他任意2個(gè)監(jiān)測點(diǎn)之間在3個(gè)數(shù)據(jù)集上的相關(guān)系數(shù),均得到大于0.9的相關(guān)系數(shù).

4.4 實(shí)證結(jié)果分析

應(yīng)用softImpute-ALS方法對全國國控點(diǎn)空氣質(zhì)量數(shù)據(jù)進(jìn)行缺失值的插補(bǔ)處理,并依據(jù)“規(guī)定”中城市濃度計(jì)算標(biāo)準(zhǔn),分析缺失值插補(bǔ)后的城市濃度值和排名的變化情況.以2021年7月為例,分別對168個(gè)城市分別計(jì)算PM2.5,PM10,SO2,NO2,CO-95%分位數(shù),O38h-90%分位數(shù)的濃度值及其排名.生態(tài)環(huán)境部在發(fā)布城市空氣質(zhì)量排名時(shí)會對監(jiān)測數(shù)據(jù)進(jìn)行校驗(yàn),校驗(yàn)后數(shù)據(jù)質(zhì)量要優(yōu)于實(shí)時(shí)監(jiān)測數(shù)據(jù),但校驗(yàn)后數(shù)據(jù)并未公開發(fā)布,本文實(shí)證分析部分基于中國環(huán)境監(jiān)測總站實(shí)時(shí)在線監(jiān)測數(shù)據(jù),僅用于分析缺失值處理方法對計(jì)算結(jié)果的影響,說明在使用公開空氣質(zhì)量監(jiān)測數(shù)據(jù)時(shí)對缺失值處理的必要性.

4.4.1 全國城市PM2.5濃度缺失值插補(bǔ)前后計(jì)算結(jié)果 20個(gè)城市中有18個(gè)相同的城市,說明缺失值插補(bǔ)前后排名總體上一致,但局部排名有明顯的變化.主要原因是2021年7月所有城市PM2.5月度均值總體較低,并且不同城市的監(jiān)測值相差非常小,即計(jì)算值的微小變化會導(dǎo)致排名的較大變化.缺失高濃度監(jiān)測值會顯著提升城市排名,而缺失低濃度監(jiān)測值會顯著降低城市排名,此時(shí)缺失值的處理方法會對城市排名產(chǎn)生關(guān)鍵影響.

表6顯示,不同城市缺失率顯著不同,最高缺失率5.2%,最低缺失率0.8%,而“規(guī)定”中并沒有明確監(jiān)測點(diǎn)小時(shí)監(jiān)測數(shù)據(jù)缺失處理方法,只是給出“以城市污染物濃度最高點(diǎn)位的數(shù)據(jù),統(tǒng)計(jì)當(dāng)日污染物濃度排名”,這樣的處理方式是不合理的,也沒有考慮污染物在一天的周期性變化特征.使用合理插補(bǔ)方法后,PM2.5月度濃度均值變化在-7.8%～ 10.1%之間,基于公開在線監(jiān)測數(shù)據(jù)進(jìn)行研究時(shí),應(yīng)當(dāng)采取合理、可靠的缺失值處理方法對缺失值進(jìn)行插補(bǔ),減少缺失值對最終計(jì)算結(jié)果的影響.

表6 缺失值插補(bǔ)前后城市PM2.5排名比較

注:(1)原排序中相同序號按照城市出現(xiàn)前后的倒序排列;(2)插補(bǔ)前后變動(dòng)計(jì)算公式為:(插補(bǔ)后-插補(bǔ)前)/插補(bǔ)前.

4.4.2 全國城市PM10濃度缺失值插補(bǔ)前后計(jì)算結(jié)果 排名后20城市中有16個(gè)相同城市,PM10缺失率從0.9%～6.0%之間,排名變化較大的是石家莊市,PM10的月度均值從47mg/m3下降至44mg/m3,從原排名163名上升至152名,上升11名,而運(yùn)城市的PM10的月度均值從42mg/m3上升至45mg/m3,排名也從150名下降至157名.

表7 缺失值插補(bǔ)前后城市PM10排名比較

注: (1)原排序中相同序號按照城市出現(xiàn)前后的倒序排名;(2)插補(bǔ)前后變動(dòng)計(jì)算公式為:(插補(bǔ)后-插補(bǔ)前)/插補(bǔ)前

4.4.3 全國城市SO2濃度缺失值插補(bǔ)前后計(jì)算結(jié)果 排名的后20城市全部相同,SO2缺失率從1.4%～ 5.5%之間.由于SO2月度均值差異較大,插補(bǔ)后排名與原排名的最后7名次序均相同,但在之后排名次序顯著不同,即當(dāng)污染物月度均值有較大差異時(shí),缺失值插補(bǔ)對排名影響較小,但如果污染物月度均值比較接近時(shí),缺失值插補(bǔ)會顯著影響排名情況.

4.4.4 全國城市NO2濃度缺失值插補(bǔ)前后計(jì)算結(jié)果 排名后20城市中有17個(gè)相同城市,NO2缺失率為1.2%～5.4%之間.與SO2排名類似,NO2月度均值差異較大,插補(bǔ)后排名與原排名的最后5名次序均相同,其他城市的排名情況均有所變化,但排名上升和下降幅度較小.

表8 缺失值插補(bǔ)前后城市SO2排名比較

注:(1)原排序中相同序號按照城市出現(xiàn)前后的倒序排名;(2)插補(bǔ)前后變動(dòng)計(jì)算公式為:(插補(bǔ)后-插補(bǔ)前)/插補(bǔ)前

表9 缺失值插補(bǔ)前后城市NO2排名比較

注:(1)原排序中相同序號按照城市出現(xiàn)前后的倒序排名;(2)插補(bǔ)前后變動(dòng)計(jì)算公式為:(插補(bǔ)后-插補(bǔ)前)/插補(bǔ)前.

表10 缺失值插補(bǔ)前后城市CO日均值95%分位數(shù)排名比較

注:(1)原排序中相同序號按照城市出現(xiàn)前后的倒序排名;(2)插補(bǔ)前后變動(dòng)計(jì)算公式為:(插補(bǔ)后-插補(bǔ)前)/插補(bǔ)前.

表11 缺失值插補(bǔ)前后城市O3日最大8h90%分位數(shù)排名比較

注:(1)原排序中相同序號按照城市出現(xiàn)前后的倒序排名;(2)插補(bǔ)前后變動(dòng)計(jì)算公式為:(插補(bǔ)后-插補(bǔ)前)/插補(bǔ)前.

4.4.5 全國城市CO日均值95%分位數(shù)和O3日最大8h 90%分位數(shù) 通過CO日均值95%分位數(shù)(表10)與O3日最大8h 90%分位數(shù)(表11)的排名比較,在數(shù)值較為接近時(shí),城市排名會受到缺失值的顯著影響,缺失值如何處理將直接影響城市排名前后次序.新余、承德、唐山、運(yùn)城和忻州這5個(gè)城市CO日均值在插補(bǔ)前后變動(dòng)均超過10%,其中運(yùn)城和忻州的分別從插值前的157名和151名,在插值后變?yōu)?38名和121名,提升20和30名,可見缺失值對城市評價(jià)濃度具有顯著影響.

5 結(jié)論

研究表明數(shù)據(jù)缺失會對計(jì)算結(jié)果造成顯著影響程度,尤其是進(jìn)行日尺度和城市尺度的統(tǒng)計(jì)分析,本研究對于提升監(jiān)測數(shù)據(jù)完整率和分析結(jié)果準(zhǔn)確率方面具有一定的現(xiàn)實(shí)意義.基于softImpute-ALS方法,對2016年以來的全國國控點(diǎn)小時(shí)級的大規(guī)模監(jiān)測數(shù)據(jù)進(jìn)行缺失值插補(bǔ)處理,模擬實(shí)驗(yàn)表明本文所使用方法能夠得到較好的估計(jì)值,RMSE、MAPE和相關(guān)系數(shù)均優(yōu)于其他缺失值常規(guī)處理方法,同時(shí)缺失值估計(jì)過程具有較快的處理速度.實(shí)證分析結(jié)果表明,缺失值插補(bǔ)前后6種污染物評價(jià)濃度會有±10%左右的變化,同時(shí)監(jiān)測點(diǎn)的缺失值數(shù)量越多,缺失值插補(bǔ)前后濃度的計(jì)算值變化也越大,越有必要進(jìn)行缺失值處理.建議在基于公開的大規(guī)?？諝赓|(zhì)量監(jiān)測數(shù)據(jù)進(jìn)行空氣質(zhì)量研究時(shí),首先應(yīng)采用softImpute-ALS方法對監(jiān)測數(shù)據(jù)進(jìn)行缺失值處理,得到更加準(zhǔn)確的估計(jì)值,提升監(jiān)測數(shù)據(jù)完整率,最大程度減少缺失值可能導(dǎo)致的結(jié)果有偏問題,提升相關(guān)研究結(jié)果的準(zhǔn)確性和可靠性.

[1] 中華人民共和國生態(tài)環(huán)境部.城市環(huán)境空氣質(zhì)量排名技術(shù)規(guī)定[R]. 2018.

Ministry of Ecology and Environmental of People's Republic of China. Technical regulations for air quality ranking of cities[R]. 2018.

[2] Deng Q, Yang K, Luo Y. Spatiotemporal patterns of PM2.5in the Beijing–Tianjin–Hebei region during 2013～2016 [J]. Geology, Ecology, and Landscapes, 2017,1(2):95-103.

[3] Li L, Wu A H, Cheng I, et al. Spatiotemporal estimation of historical PM2.5concentrations using PM10, meteorological variables, and spatial effect [J]. Atmospheric Environment, 2017,166:182-191.

[4] Hu M, Wang Y, Wang S, et al. Spatial-temporal heterogeneity of air pollution and its relationship with meteorological factors in the Pearl River Delta, China [J]. Atmospheric Environment, 2021,254:118415.

[5] Li L, Zhang J, Meng X, et al. Estimation of PM2.5concentrations at a high spatiotemporal resolution using constrained mixed-effect bagging models with MAIAC aerosol optical depth [J]. Remote Sensing of Environment, 2018,217:573-586.

[6] Shen Y, Zhang l, Fang X, et al. Spatiotemporal patterns of recent PM2.5concentrations over typical urban agglomerations in China [J]. Science of the Total Environment, 2019,655:13-26.

[7] Zhao S, Yin D, Yu Y, et al. PM2.5and O3pollution during 2015～2019 over 367 Chinese cities: Spatiotemporal variations, meteorological and topographical impacts [J]. Environmental Pollution, 2020,264:114694.

[8] Li K, Jacob D J, Liao H, et al. A two-pollutant strategy for improving ozone and particulate air quality in China [J]. Nature Geoscience, 2019, 12(11):906-910.

[9] Liu J, Li W, Wu J. A framework for delineating the regional boundaries of PM2.5pollution: A case study of China [J]. Environmental Pollution, 2018,235:642-651.

[10] 張 烴,董樹屏,滕 曼,等.區(qū)域大型環(huán)境空氣綜合觀測中外場觀測與實(shí)驗(yàn)室分析數(shù)據(jù)質(zhì)量控制研究[J]. 環(huán)境科學(xué)研究, 2019,32(10): 1664-1671.

Zhang T, Dong S P, Teng M, et al. Quality assurance of field observation and laboratory analysis in regional large scale ambient air joint monitoring campaigns [J]. Research of Environmental Sciences, 2019,32(10):1664-1671.

[11] 師耀龍,呂怡兵,肖建軍.夏季重大活動(dòng)期間O3監(jiān)測數(shù)據(jù)質(zhì)量提升方法研究[J]. 中國環(huán)境監(jiān)測, 2020,36(2):10-14.

Shi Y L, Lyu Y B, Xiao J J. Data quality control method of ozone monitoring during the guarantee for major events in summer [J]. Environmental Monitoring in China, 2020,36(2):10-14.

[12] 師耀龍,陳傳忠,魏俊山,等.加強(qiáng)生態(tài)環(huán)境監(jiān)測機(jī)構(gòu)監(jiān)督管理的思考與分析[J]. 環(huán)境保護(hù), 2018,46(23):56-60.

Shi Yao-long, Chen Chuan-zhong, Wei Jun-shan, et al. The current situation and problem analysis of environmental monitoring organizations' supervision and administration [J]. Environmental Protection, 2018,46(23):56-60.

[13] 劉 媛,彭 溶,張 馳,等.環(huán)境監(jiān)測從業(yè)人員監(jiān)管制度研究[J]. 環(huán)境保護(hù), 2018,46(18):33-35.

Liu Y, Peng R, Zhang C, et al. Research on supervision system of environmental monitoring practitioners [J]. Environmental Protection, 2018,46(18):33-35.

[14] 師耀龍,楊 婧,柴文軒,等.美國環(huán)境空氣監(jiān)測數(shù)據(jù)質(zhì)量核查工作的經(jīng)驗(yàn)與啟示[J]. 中國環(huán)境監(jiān)測, 2017,33(3):8-14.

Shi Y L, Yang J, Chai W X, et al. Experience and illumination of data quality assessment system for ambient air monitoring in the United States [J]. Environmental Monitoring in China, 2017,33(3):8-14.

[15] Rumaling M I, Chee F pien, Dayou J, et al. Missing value imputation for PM10concentration in Sabah using nearest neighbour method (NNM) and expectation-maximization (EM) algorithm [J]. Asian Journal of Atmospheric Environment, 2020,14:62-72.

[16] Junger W L, Ponce D E Leon A. Imputation of missing data in time series for air pollutants [J]. Atmospheric Environment, 2015,102:96- 104.

[17] Larsen L C, Shah M. A context-intensive approach to imputation of missing values in data sets from networks of environmental monitors [J]. Journal of the Air & Waste Management Association (1995), 2016,66(1):38-52.

[18] Junninen H, Niska H, Tuppurainen K, et al. Methods for imputation of missing values in air quality data sets [J]. Atmospheric Environment, 2004,38(18):2895-2907.

[19] Hadeed S J, O’rourke M K, Burgess J L, et al. Imputation methods for addressing missing data in short-term monitoring of air pollutants [J]. Science of The Total Environment, 2020,730:139140.

[20] Real C, ángel Fernández J, Aboal J R, et al. Substituting missing data in compositional analysis [J]. Environmental Pollution, 2011,159(10): 2797-2800.

[21] Gómez-Carracedo M P, Andrade J M, López-Mahía P, et al. A practical comparison of single and multiple imputation methods to handle complex missing data in air quality datasets [J]. Chemometrics and Intelligent Laboratory Systems, 2014,134:23-33.

[22] Chen X, Xiao Y. A novel method for air quality data imputation by nuclear norm minimization [J]. Journal of Sensors, 2018,2018: e7465026.

[23] Moshenberg S, Lerner U, Fishbain B. Spectral methods for imputation of missing air quality data [J]. Environmental Systems Research, 2015,4(1):26.

[24] Liu X, Wang X, Zou L, et al. Spatial imputation for air pollutants data sets via low rank matrix completion algorithm [J]. Environment International, 2020,139:105713.

[25] Hastie T, Mazumder R, Lee J D, et al. Matrix completion and low-rank SVD via fast alternating least squares [J]. 36.

[26] Candès E J, Recht B. Exact matrix completion via convex optimization [J]. Foundations of Computational Mathematics, 2009, 9(6):717.

[27] Candès E J, Tao T. The power of convex relaxation: near-optimal matrix completion [J]. IEEE Transactions on Information Theory, 2010,56(5):2053-2080.

[28] Buuren S van. Flexible imputation of missing data [M]. 2nd edition. Boca Raton: Chapman and Hall/CRC, 2018.

[29] Liu Y, Dillon T, Yu W, et al. Missing value imputation for industrial IoT sensor data with large gaps [J]. IEEE Internet of Things Journal, 2020,7(8):6855-6867.

[30] Velasco-Gallego C, Lazakis I. A novel framework for imputing large gaps of missing values from time series sensor data of marine machinery systems [J]. Ships and Offshore Structures, 2021,10.1080/ 17445302.2021.1943850.

[31] 中華人民共和國生態(tài)環(huán)境部.2021年8月全國城市空氣質(zhì)量報(bào)告[R]. 2021.

Ministry of Ecology and Environmental of People's Republic of China. Air quality reports for cities in China[R]. 2021,8

致謝：本研究受中國人民大學(xué)“雙一流”跨學(xué)科重大創(chuàng)新規(guī)劃平臺—生態(tài)文明跨學(xué)科交叉平臺支持.

Research on the missing value methods for large-scale online air quality monitoring data.

ZHANG Bo, SONG Guo-jun*

(School of Environment and Natural Resources, Renmin University of China, Beijing 100872, China)., 2022,42(5)：2078～2087

Large scale online air quality monitoring data is the basis for air quality research, but there were lots of missing data in large scale online data. In this study, we compared several methods that dealing with the missing values and its impact on the city’s ranking of air quality base on the hourly monitoring data of 1654monitoring sites in China from 1Jan, 2016 to 31July, 2021 of 6types of air pollutants. The simulation results showed that Low Rank SVD via Alternating Least Square had smaller mean squared error, mean absolute percentage error and higher correlation coefficient compared with other traditional methods. The empirical results showed there would be 10% difference before imputation and after imputation for the missing value. The ranking would not change due to the imputation when the air quality assessed value vary greatly, and would change a lot when the assessed value was very close. The study suggested to impute missing value by using the method in this study when analysis the large-scale online air quality monitoring data.

monitoring data；big data；missing value；low rank matrix；empirical research

X323

A

1000-6923(2022)05-2078-10

張 波(1986-),男,內(nèi)蒙古包頭人,講師,博士,研究方向?yàn)榄h(huán)境統(tǒng)計(jì)與建模.發(fā)表論文10余篇.

2021-10-12

中國人民大學(xué)科學(xué)研究基金(中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助)項(xiàng)目成果(22XNF016)

* 責(zé)任作者, 教授, songguojun@ruc.edu.cn