羅 蒙，樂 群，張 新

(華東師范大學(xué) 地理科學(xué)學(xué)院/地理信息系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，上海 200241)

?

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

羅 蒙，樂 群，張 新

(華東師范大學(xué) 地理科學(xué)學(xué)院/地理信息系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，上海 200241)

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

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

1引言

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

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

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

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

2.1資料

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

2.2方法

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

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

3結(jié)果與分析

3.1檢驗(yàn)效果評(píng)估

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

表2　個(gè)例分析站點(diǎn)及其異常值情況

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

表3　各檢驗(yàn)法對(duì)10個(gè)個(gè)例站點(diǎn)異常值檢出情況

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

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

表4　各檢驗(yàn)法對(duì)全部站點(diǎn)的異常大值檢驗(yàn)結(jié)果

3.2最佳檢驗(yàn)方法

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

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

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

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

4討論與結(jié)論

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

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

參考文獻(xiàn)：

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

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

[3]潘本鋒, 鄭皓皓, 李莉娜,等. 空氣自動(dòng)監(jiān)測(cè)中PM2.5與PM10“倒掛”現(xiàn)象特征及原因[J]. 中國(guó)環(huán)境監(jiān)測(cè), 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

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

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

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

中圖分類號(hào)：X831

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：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