羅 蒙，樂 群，張 新

(華東師范大學 地理科學學院/地理信息系統(tǒng)教育部重點實驗室，上海 200241)

?

不一致性檢驗法對細顆粒物監(jiān)測數(shù)據(jù)中異常大值的檢驗

羅 蒙，樂 群，張 新

(華東師范大學 地理科學學院/地理信息系統(tǒng)教育部重點實驗室，上海 200241)

摘要：根據(jù)目前國內(nèi)關(guān)于細顆粒物監(jiān)測數(shù)據(jù)質(zhì)量控制方面研究很少的現(xiàn)狀，使用Barnett總結(jié)出的正態(tài)樣本不一致性檢驗法以及Verma模擬出的超大樣本臨界值，對可吸入顆粒物小時監(jiān)測數(shù)據(jù)進行了異常值檢驗。所用數(shù)據(jù)為環(huán)保部發(fā)布的可吸入顆粒物小時監(jiān)測數(shù)據(jù)，數(shù)據(jù)時長為一年半左右，結(jié)果表明：各種檢驗法中以N1和N4為代表的偏差/尺度型統(tǒng)計量和平方和型統(tǒng)計量較適合PM2.5監(jiān)測數(shù)據(jù)的異常值檢驗；N1檢驗法比N4檢驗法更不容易受臨界值插值誤差的影響。

關(guān)鍵詞：PM2.5；小時監(jiān)測數(shù)據(jù)；異常值檢驗；不一致性檢驗

1引言

在諸多空氣污染物中，細顆粒物(PM2.5)因其直徑甚微，吸入后可直接進入肺部，對人類身體健康造成較大威脅，近年來受到社會各界的廣泛關(guān)注。目前我國已在多數(shù)城市設(shè)立了PM10和PM2.5的監(jiān)測站點，并計劃將監(jiān)測網(wǎng)絡(luò)覆蓋至所有地級以上城市，監(jiān)測方式由過去的人工稱重測量(振蕩天平法)逐漸演變?yōu)楝F(xiàn)在的自動監(jiān)測(β射線法)[1]。然而，我國目前的監(jiān)測網(wǎng)絡(luò)還不夠完善，數(shù)據(jù)積累時長較短，站點少且分布不均勻，數(shù)據(jù)質(zhì)量較差。在監(jiān)測過程中，由于受儀器故障、人為測量失誤、惡劣天氣、數(shù)據(jù)傳輸故障等因素影響，監(jiān)測數(shù)據(jù)會出現(xiàn)異常值，往往表現(xiàn)為異常大值、負值，或數(shù)值異常起伏、長時間平緩監(jiān)測結(jié)果[2]或PM2.5濃度大于PM10[3]，若在分析過程中不加以仔細甄別，勢必會對研究結(jié)果造成較大影響，甚至得出錯誤的結(jié)論。

在統(tǒng)計學領(lǐng)域，異常值的檢驗始終是一個重要而復雜的問題。Bendre給出了指數(shù)分布樣本中屏蔽效應(yīng)(masking effect)的范例，其中屏蔽效應(yīng)是指異常值檢驗中常常發(fā)生的一種因為有其他異常值存在而導致某些異常值無法被識別的現(xiàn)象[4]。關(guān)于正態(tài)樣本異常點的研究成果最為詳盡，如Barnett在處理符合正態(tài)分布的樣本時，認為明顯偏離樣本分布特征的值是異常值，檢驗這些異常值的方法被稱為不一致性檢驗(discordancy test)[5]。為了拓展檢驗法的應(yīng)用范圍，Verma利用蒙特卡洛方法對15種不一致性檢驗法(以N1-N15指代)在7個顯著水平下(α=0.3，0.2，0.1,0.05,0.02，0.01，0.005)的臨界值表，以及臨界值隨樣本量n變化的插值公式，將檢驗法的應(yīng)用范圍擴大到容量最多為30000的超大樣本，極大地拓展了檢驗法的應(yīng)用范圍[6～10]。

由于PM2.5小時監(jiān)測數(shù)據(jù)樣本容量很大，許多只適用于小樣本的異常值檢驗法難以應(yīng)用于長時段的監(jiān)測數(shù)據(jù)，而Verma的工作剛好拓展了不一致性檢驗的應(yīng)用范圍，使其可以應(yīng)用于最多長達3年左右的小時監(jiān)測數(shù)據(jù)時間序列。針對目前國內(nèi)細顆粒物監(jiān)測數(shù)據(jù)質(zhì)量較差、關(guān)于異常值檢驗方面研究較少的現(xiàn)狀，以Verma改進的正態(tài)樣本不一致性檢驗方法為基礎(chǔ)，通過對比選出最適合的PM2.5小時監(jiān)測數(shù)據(jù)特征的檢驗方法，主要檢驗數(shù)據(jù)中存在的異常大值錯誤數(shù)據(jù)。

2數(shù)據(jù)來源和方法

2.1資料

選用數(shù)據(jù)為中國環(huán)境保護部公布的PM2.5和PM10小時監(jiān)測數(shù)據(jù)，共有363個城市的1575個國家空氣質(zhì)量自動監(jiān)測站點，數(shù)據(jù)收集起止時間為2013年11月至2015年5月，其中有兩段較長時段的數(shù)據(jù)因網(wǎng)絡(luò)故障沒有收集(2013年11月26日12時至2013年12月5日15時，2015年2月7日11時至2015年3月2日6時)，作為缺測時段處理。因PM10缺測率較高，僅對PM2.5監(jiān)測數(shù)據(jù)進行異常值檢驗，而將PM10數(shù)據(jù)作為異常值人工檢視的參考數(shù)據(jù)。

2.2方法

重點探討監(jiān)測序列中的出現(xiàn)的異常大值，故略去了各方法檢測最小值一側(cè)的統(tǒng)計量。在檢驗之前首先剔除所有超過測量范圍的錯誤值?？紤]到PM2.5監(jiān)測數(shù)據(jù)往往不符合正態(tài)分布[11，12]，先對原始序列進行自然對數(shù)轉(zhuǎn)換，再進行異常值檢驗。

表1　 6種Barnett正態(tài)檢驗法的概況

3結(jié)果與分析

3.1檢驗效果評估

為更好地評估檢驗法的檢驗功效，我們從1575個時間序列中選出有代表性的8個含異常值站點進行個例分析(表2)。這些異常值通常數(shù)值很大，或前后一段時間內(nèi)序列比較平穩(wěn)，只在異常點處有尖峰，或異常大值持續(xù)數(shù)個至數(shù)十個小時，與異常值相鄰的時次往往存在缺測，這些特征都與正常序列有較大差別。為比較異常值與正確大值記錄的區(qū)別，同時也選取了兩個由污染過程導致的PM2.5高濃度序列作為對比，其一是喀什地區(qū)市環(huán)境監(jiān)測站在2015年5月的一段監(jiān)測序列，期間發(fā)生一次強沙塵暴；其二是株洲天臺山莊在2014年1月31日凌晨發(fā)生的一次強污染過程，1月30日為大年除夕夜，在全國其他很多站點都發(fā)現(xiàn)該時段出現(xiàn)大監(jiān)測值，因此認定這是大規(guī)模燃放煙花爆竹導致的污染。

表2　個例分析站點及其異常值情況

注：*該時段內(nèi)除異常值外，其余時次為缺測

表3　各檢驗法對10個個例站點異常值檢出情況

注：括號中為檢出右側(cè)異常大值的個數(shù)

將6種方法應(yīng)用于所有站點，檢驗結(jié)果見表4，可見，仍然是N1與N4的檢驗結(jié)果最接近人工檢視發(fā)現(xiàn)的異常值數(shù)；N9對多異常值站點漏檢嚴重；N6存在較多誤檢，但對于一些存在多個異常值的站點又存在漏檢；N14檢出了大量的異常小值，雖然也檢出了許多異常大值，然而大部分是誤檢，許多真正的異常值并沒有被檢出；N15則存在較多的誤檢，檢驗效果介于幾種方法之間。因此認為N1和N4是最適合PM2.5小時監(jiān)測數(shù)據(jù)的異常值檢驗方法。N1和N4的異常值檢出數(shù)比人工檢視少，這是因為異常值中有一部分處于次序樣本的中部，即這些異常值不在樣本最大值一側(cè)，無法被不一致性檢驗法識別。這部分異常值多屬于樣本局部的跳變，可用濾波等其他方法予以排除，不予詳細討論。

表4　各檢驗法對全部站點的異常大值檢驗結(jié)果

3.2最佳檢驗方法

排除4個不適合的檢驗法后，再來詳細對比N1與N4的檢驗效果。圖1為1%顯著水平下對各站應(yīng)用N1方法檢驗最大值項是否為異常值的情況。圖2與圖1類似，為相同條件下N4方法的情況?？梢?，雖然N1和N4的檢驗結(jié)果較為接近，但TN1與臨界值的距離要遠遠大于TN4與臨界值的距離，而當統(tǒng)計量與臨界值十分接近時，容易受檢驗方法誤差和臨界值插值誤差等因素影響，被檢出的異常值可能難與正常的大值記錄區(qū)分，從而發(fā)生誤檢或漏檢情況。TN4與臨界值十分接近是由其算法決定的，當樣本容量較大時，去掉一項對整體樣本的離差平方和影響較小，因此TN4一般是小于1且非常接近于1的分數(shù)，且與臨界值的偏差通常小于0.01。

圖1　各站點應(yīng)用N1檢驗法時最大值x(n)是否為異常值的檢驗結(jié)果與相應(yīng)的統(tǒng)計量TN1和臨界值

通過人工檢視N1與N4的檢驗結(jié)果發(fā)現(xiàn)，N1檢出值比N4稍多，誤檢正確值的情況多于N4，但N1漏檢明顯異常值的情況遠遠少于N4，由于誤檢正確情況可以經(jīng)由人工檢視進行排除，而漏檢則無法進行補救，故N1比N4更好。綜合來看，N1是最佳檢驗方法。

圖2　各站點應(yīng)用N4檢驗法時最大值x(n)是否為異常值的檢驗結(jié)果與相應(yīng)的統(tǒng)計量TN4和臨界值

4討論與結(jié)論

N1方法雖然能檢驗出大部分異常值，但會受到樣本容量、統(tǒng)計量和臨界值的制約，對異常大值的檢驗能力存在最低閾值，低于閾值的異常值無法被識別。例如，若樣本均值和標準差過大，則由表1可知，其對應(yīng)的能識別的最小x(n)也偏大，因此在實際應(yīng)用中，需要對原始數(shù)據(jù)的基本統(tǒng)計特征有一定了解后才能應(yīng)用該方法，對于樣本均值和標準差過大的站點可能會因檢驗閾值過高而漏掉部分異常值，此時不宜使用N1檢驗法。

6種Barnett總結(jié)的正態(tài)樣本不一致性檢驗方法中，N1和N4檢驗效果與實際情況最為接近；N6方法本身無法判斷異常值出現(xiàn)在哪一側(cè)，且易受屏蔽效應(yīng)影響；N9能檢驗出只存在一個異常值的站點，但易受屏蔽效應(yīng)影響無法檢出存在多個異常值的站點；N14的統(tǒng)計量形式?jīng)Q定其會將左側(cè)異常小值識別為異常值，與研究目的不符，不適合PM2.5的異常值檢驗；N15檢驗效果介于N1、N4和N6之間。N1與N4相比，N1誤檢情況多于N4，漏檢情況少于N4，且統(tǒng)計量TN1與臨界值的距離大于TN4與臨界值的距離，檢驗結(jié)果不容易受到臨界值誤差的影響。綜上所述，N1為最適合PM2.5小時監(jiān)測數(shù)據(jù)的異常大值檢驗方法。

參考文獻：

[1]潘本鋒,汪巍,王瑞斌,等.我國PM2.5監(jiān)測網(wǎng)絡(luò)布局與監(jiān)測方法體系構(gòu)建策略分析[J].環(huán)境與可持續(xù)發(fā)展,2013,38(3):9～13.

[2]師建中,陳丹青.PM2.5監(jiān)測數(shù)據(jù)質(zhì)量主要影響因素和控制方法探討[J].綠色科技,2012(5):179～180.

[3]潘本鋒, 鄭皓皓, 李莉娜,等. 空氣自動監(jiān)測中PM2.5與PM10“倒掛”現(xiàn)象特征及原因[J]. 中國環(huán)境監(jiān)測, 2014,30(5):90～95.

[4]Bendre S M, Kale B K. Masking effect on tests for outliers in exponential models[J]. Journal of the American Statistical Association, 1985, 80(392): 1020～1025.

[5]Barnett V, Lewis T. Outliers in Statistical Data[M]. Chichester: John Wiley, 1978.

[6]Verma S P. Sixteen statistical tests for outlier detection and rejection in evaluation of International Geochemical Reference Materials: Example of microgabbro PM‐S[J]. Geostandards Newsletter, 1997, 21(1): 59～75.

[7]Verma S P, Quiroz Ruiz A. Critical values for six Dixon tests for outliers in normal samples up to sizes 100, and applications in science and engineering[J]. Revista Mexicana de CienciasGeológicas, 2006, 23(2): 133～161.

[8]Verma S P, Quiroz Ruiz A. Critical values for 22 discordancy test variants for outliers in normal samples up to sizes 100, and applications in science and engineering[J]. Revistamexicana de CienciasGeológicas, 2006, 23(3): 302～319.

[9]Verma S P, Quiroz Ruiz A, Díaz-González L. Critical values for 33 discordancy test variants for outliers in normal samples up to sizes 1000, and applications in quality control in Earth Sciences[J]. Revista Mexicana de CienciasGeológicas, 2008, 25(1): 82～96.

[10]Verma S P, Quiroz-Ruiz A. Critical values for 33 discordancy test variants for outliers in normal samples of very large sizes from 1000 to 30000 and evaluation of different regression models for the interpolation and extrapolation of critical values[J]. Revista Mexicana de CienciasGeológicas, 2008, 25(3): 369～381.

[11]Karaca F, Alagha O, Ertürk F. Statistical characterization of atmospheric PM 10 and PM 2.5 concentrations at a non-impacted suburban site of Istanbul, Turkey[J]. Chemosphere, 2005, 59(8): 1183～1190.

[12]Lu H C, Fang G C. Estimating the frequency distributions of PM 10 and PM 2.5 by the statistics of wind speed at Sha-Lu, Taiwan[J]. Science of the total environment, 2002, 298(1): 119～130.

收稿日期：2016-05-09

基金項目：國家重點基礎(chǔ)研究發(fā)展規(guī)劃973項目(編號：2012CB955803)；華東師范大學大型儀器設(shè)備開放基金；華東師范大學研究生科研創(chuàng)新實踐資助項目(編號：YJSKC2015-15)

作者簡介：羅蒙(1989—)，男，華東師范大學地理科學學院碩士研究生。

通訊作者：樂群(1967—)，男，副教授，博士，主要從事氣候數(shù)值模擬及城市大氣環(huán)境方面的教學與研究工作。

中圖分類號：X831

文獻標識碼：A

文章編號：1674-9944(2016)12-0129-04

Upper Outlier Detection of Fine Particulate Matter Monitoring DataUsing the Discordancy Tests

Luo Meng, Yue Qun, Zhang Xin

(SchoolofGeographicSciences/KeyLaboratoryofGeographicInformationScience,MinistryofEducation,EastChinaNormalUniversity,Shanghai200241,China)

Abstract:The quality control studies of fine particulate matter monitoring data in China is very limited at present.Based on the Normal Sample discordance tests summarized by Barnett and the critical values for super-large sized samples simulated by Verma, the outliers in fine particulate matter monitoring data had been detected.The hourly monitoring data of inhalable particles used here was released by Ministry of Environmental Protection of the People’s Republic of China.The data length is about one and a half year. Among the various detection method,the N1 and N4 test had the best results,which representeddeviation/spread statistics and sums of squares statistics,respectively.N1 was less vulnerable to the interpolation error of critical values than N4.

Key words:PM2.5; hourly monitoring data;outlier detection;discordance test