師耀龍，滕 曼，李 成，吳曉鳳，柴文軒，楊 婧，楚寶臨,付 強(qiáng)

1.中國環(huán)境監(jiān)測總站，國家環(huán)境保護(hù)環(huán)境監(jiān)測質(zhì)量控制重點(diǎn)實(shí)驗(yàn)室，北京 1000122.河北農(nóng)業(yè)大學(xué)科學(xué)技術(shù)研究院,河北 保定 071001

應(yīng)用穩(wěn)健馬氏距離評價(jià)實(shí)驗(yàn)室間比對能力考核

師耀龍1，滕 曼1，李 成2，吳曉鳳1，柴文軒1，楊 婧1，楚寶臨1,付 強(qiáng)1

1.中國環(huán)境監(jiān)測總站，國家環(huán)境保護(hù)環(huán)境監(jiān)測質(zhì)量控制重點(diǎn)實(shí)驗(yàn)室，北京 1000122.河北農(nóng)業(yè)大學(xué)科學(xué)技術(shù)研究院,河北 保定 071001

基于穩(wěn)健馬氏距離的多元統(tǒng)計(jì)算法，綜合評價(jià)了全國環(huán)境監(jiān)測系統(tǒng)96個(gè)實(shí)驗(yàn)室4種有機(jī)氯農(nóng)藥監(jiān)測能力考核的結(jié)果，從中篩選出25個(gè)考核結(jié)果存在異常的實(shí)驗(yàn)室。將多元統(tǒng)計(jì)結(jié)果與穩(wěn)健Z比分?jǐn)?shù)結(jié)果進(jìn)行比較，發(fā)現(xiàn)基于穩(wěn)健馬氏距離篩選出的異常實(shí)驗(yàn)室其穩(wěn)健Z比分?jǐn)?shù)結(jié)果也存在一定的問題。同時(shí)簡單介紹了穩(wěn)健馬氏距離及其相關(guān)的異常值篩選的算法原理，并詳細(xì)介紹了其基于R語言的實(shí)現(xiàn)過程。

多元穩(wěn)健統(tǒng)計(jì)；穩(wěn)健馬氏距離；實(shí)驗(yàn)室間比對；能力考核

環(huán)境監(jiān)測實(shí)驗(yàn)室能力考核是中國環(huán)境監(jiān)測總站(以下簡稱總站)組織的利用實(shí)驗(yàn)室間比對的方式考核環(huán)境監(jiān)測實(shí)驗(yàn)室技術(shù)能力和質(zhì)量管理水平的活動(dòng)。通過實(shí)驗(yàn)室間比對，有助于總站了解領(lǐng)域內(nèi)實(shí)驗(yàn)室相關(guān)監(jiān)測項(xiàng)目的能力狀況，及時(shí)發(fā)現(xiàn)和糾正影響檢測水平的影響因素，促進(jìn)領(lǐng)域內(nèi)實(shí)驗(yàn)室檢測能力的提高[1-2]。

伴隨著能力考核工作的開展，其數(shù)據(jù)評價(jià)工作中出現(xiàn)了2個(gè)主要問題：①一些實(shí)驗(yàn)室在周期內(nèi)參加了多輪能力考核，或在一輪能力考核中進(jìn)行了多個(gè)樣品的測試，傳統(tǒng)的評價(jià)方法主要是對各實(shí)驗(yàn)室逐個(gè)項(xiàng)目或逐個(gè)輪次的評價(jià)[1-2]，如何用一項(xiàng)綜合指標(biāo)整體評價(jià)實(shí)驗(yàn)室能力尚待研究；②由于能力考核結(jié)果不服從或近似服從正態(tài)分布，需要對傳統(tǒng)統(tǒng)計(jì)量進(jìn)行穩(wěn)健化，如果希望整體評價(jià)，則需要對其參與的各個(gè)項(xiàng)目構(gòu)成的多元數(shù)據(jù)進(jìn)行穩(wěn)健統(tǒng)計(jì)，而穩(wěn)健Z比分?jǐn)?shù)和迭代法都是對能力考核的一元數(shù)據(jù)進(jìn)行穩(wěn)健統(tǒng)計(jì)的方法[3]，對能力考核的多元數(shù)據(jù)進(jìn)行穩(wěn)健統(tǒng)計(jì)的方法尚待研究。

鑒于此，研究詳細(xì)討論了穩(wěn)健馬氏距離在有機(jī)氯監(jiān)測能力考核結(jié)果評價(jià)中的應(yīng)用，并與穩(wěn)健Z比分?jǐn)?shù)評價(jià)結(jié)果進(jìn)行了比較[1-4]。此外，由于馬氏距離與多元數(shù)據(jù)的穩(wěn)健統(tǒng)計(jì)方法在中國環(huán)境監(jiān)測領(lǐng)域應(yīng)用較少，研究對其一般原理與基于開源軟件(R)的實(shí)現(xiàn)方法進(jìn)行了詳細(xì)介紹。

1 一般馬氏距離

在能力考核多元數(shù)據(jù)處理中，馬氏距離具有：①不受不同項(xiàng)目之間量綱不同的影響，所得距離為標(biāo)準(zhǔn)化距離；②計(jì)算過程考慮到了各項(xiàng)目之間相關(guān)性的存在，是建立在總體協(xié)方差矩陣上的標(biāo)準(zhǔn)化距離。馬氏距離已應(yīng)用于包括環(huán)境監(jiān)測領(lǐng)域在內(nèi)的各個(gè)領(lǐng)域內(nèi)多元數(shù)據(jù)異常值篩選工作中[6]。但是，傳統(tǒng)的統(tǒng)計(jì)方法中，異常值的存在會顯著影響中心值和協(xié)方差矩陣的估計(jì)，使一般馬氏距離不能正確反映各個(gè)觀測的偏離程度。對于這類數(shù)據(jù)，需要通過穩(wěn)健統(tǒng)計(jì)的方法，構(gòu)建穩(wěn)定的均值和協(xié)方差矩陣統(tǒng)計(jì)量。

2 穩(wěn)健馬氏距離及其相關(guān)R程序包介紹

開源R軟件中的“robustbase”程序包[7-8]中的covMcd程序是基于Fast-MCD算法和PISON等人在2002年的改進(jìn)完成的[7, 9-11]，廣泛用于估計(jì)多元樣本中的穩(wěn)健統(tǒng)計(jì)量。其算法原理較為復(fù)雜，可簡單的概括為從樣本中選擇h個(gè)觀測，通過不斷的迭代計(jì)算h個(gè)觀測最小的協(xié)方差行列式，該協(xié)方差矩陣通過加權(quán)即可估計(jì)出穩(wěn)健的統(tǒng)計(jì)量，詳細(xì)的計(jì)算方法參見文獻(xiàn)[7,10-12]。穩(wěn)健估計(jì)結(jié)束后，以穩(wěn)健的中心值和協(xié)方差，通過馬氏距離計(jì)算各觀測向量偏離中心值的穩(wěn)健馬氏距離，并可根據(jù)DM(x)2符合卡方分布的特點(diǎn)，篩選出數(shù)據(jù)集中的異常值，這一過程可由R軟件中的“mvoutlier”程序包實(shí)現(xiàn)[13]。

pn(δ)=sup(G(u)-Gn(u))+u≥δ

同時(shí)，計(jì)算pcrit(δ,n,p)

計(jì)算完成后，比較pn(δ)和pcrit(δ,n,p)，若pn(δ)pcrit(δ,n,p)，則在其尾部有pn(δ)×100%個(gè)異常值(從穩(wěn)健馬氏距離最大的觀測算起)。

“mvoutlier”程序包同時(shí)支持多種計(jì)算和畫圖功能，并能從單個(gè)變量觀測整個(gè)數(shù)據(jù)集，分析哪些變量更容易導(dǎo)致異常值的出現(xiàn)，可用于能力考核中不合格站點(diǎn)的篩選與考核項(xiàng)目難易的判斷。

研究將通過穩(wěn)健馬氏距離在有機(jī)氯監(jiān)測能力考核中的應(yīng)用實(shí)例來介紹穩(wěn)健馬氏距離在實(shí)驗(yàn)室能力考核中的應(yīng)用。

3 穩(wěn)健馬氏距離在環(huán)境監(jiān)測實(shí)驗(yàn)室能力考核中的應(yīng)用

3.1 數(shù)據(jù)來源與R軟件計(jì)算過程

總站于2013年開展了針對全國各地市環(huán)境監(jiān)測站的有機(jī)氯監(jiān)測能力考核，數(shù)據(jù)選取了96個(gè)實(shí)驗(yàn)室測定的α-六氯環(huán)己烷(簡寫為α-666)樣品A、α-666樣品B、γ-六氯環(huán)己烷(簡寫為γ-666)樣品A、γ-666樣品B、p, p′-雙對氯苯基三氯乙烷(簡寫為p, p′-DDT)樣品A、p, p′-DDT樣品B、o, p′-雙對氯苯基三氯乙烷(簡寫為o, p′-DDT)樣品A和o, p′-DDT樣品B的含量，構(gòu)成一個(gè)96×8的矩陣，R計(jì)算過程如下：

>x<-read.table(file.choose(), header=TRUE, row.names=1)#導(dǎo)入txt格式矩陣

>library(robustbase)#調(diào)用robustbase

>library(mvoutlier)#調(diào)用mvoutlier>covMcd(x, alpha=0.75) #估計(jì)穩(wěn)健統(tǒng)計(jì)量(h=0.75)

>cov(x) #估計(jì)一般協(xié)方差

>apply(x, 2, mean) #估計(jì)一般中心值

>res1<-dd.plot(x, quan=0.75) #計(jì)算一般馬氏距離和穩(wěn)健馬氏距離

>mdc<-res1$md.cla

>write.table(mdc, file="classical_md.txt")#保存一般馬氏距離

>mdr<-res1$md.rob

>write.table(mdr, file="robust_md.txt")#保存穩(wěn)健馬氏距離

>res2<-aq.plot(x, quan=0.75, alpha=0.05)#篩選不合格機(jī)構(gòu)

>outliers<-which(res2$outliers==T)

>write.table(outliers,file="outliers.txt")#保存不合格機(jī)構(gòu)

>uni.plot(x, quan=0.75, alpha=0.05)#從各個(gè)項(xiàng)目觀測整體數(shù)據(jù)偏離情況

>res1<-covMcd(x, alpha=0.75)

>res2<-arw(x,res1$center,res1$cov,0.05)

>sqrt(res2$cn) #計(jì)算判斷點(diǎn)位數(shù)據(jù)異常的臨界值，若某機(jī)構(gòu)馬氏距離≥該值，則表明該機(jī)構(gòu)能力考核結(jié)果偏離中心值較遠(yuǎn)，結(jié)果存在異常。

>q()

計(jì)算完成后，各監(jiān)測機(jī)構(gòu)測定結(jié)果偏離中心值的穩(wěn)健馬氏距離(robust_md.txt)、傳統(tǒng)馬氏距離(classical_md.txt)和基于穩(wěn)健馬氏距離檢測得到不合格站點(diǎn)名稱(outliers.txt)的分別儲存在R軟件默認(rèn)的工作目錄下，將25個(gè)不合格機(jī)構(gòu)信息匯總進(jìn)入原始數(shù)據(jù)集中進(jìn)一步匯總分析。

3.2 結(jié)果與討論

96家機(jī)構(gòu)8種樣品均值、協(xié)方差矩陣、穩(wěn)健中心值和穩(wěn)健協(xié)方差矩陣見表1～表3。比較穩(wěn)健與非穩(wěn)健均值和協(xié)方差矩陣后發(fā)現(xiàn)，由于94、95、96 3個(gè)極值的存在(表4)，嚴(yán)重影響了該數(shù)據(jù)集的均值和協(xié)方差。說明通過穩(wěn)健統(tǒng)計(jì)估計(jì)出的穩(wěn)健中心值和穩(wěn)健協(xié)方差矩陣可以排除極值對正確估計(jì)數(shù)據(jù)集統(tǒng)計(jì)量的影響，更能體現(xiàn)有機(jī)氯能力考核數(shù)據(jù)的正常分布情況。

表2 96家監(jiān)測機(jī)構(gòu)穩(wěn)健協(xié)方差矩陣Table 2 The covariance of 96 labs evaluated by robust multi-statistics

表3 96家監(jiān)測機(jī)構(gòu)非穩(wěn)健協(xié)方差矩陣Table 3 The covariance of 96 labs evaluated by non-robust multi-statistics

表4 數(shù)據(jù)異常的25家監(jiān)測機(jī)構(gòu)穩(wěn)健馬氏距離、傳統(tǒng)馬氏距離、原始濃度和穩(wěn)健Z比分?jǐn)?shù)判定結(jié)果匯總Table 4 The robust mahalanobis distance, classical mehalanobis distance, concentrations and Z-scores of 25 outlier labs

由表4可見，當(dāng)置信水平a=0.05時(shí)，通過“mvoutlier”程序包計(jì)算出該次能力考核判定合格與否的臨界值為4.3，共篩選出25個(gè)穩(wěn)健馬氏距離≥4.3的能力考核結(jié)果存在異常的機(jī)構(gòu)(機(jī)構(gòu)名稱用序號表示)。94、95、96站點(diǎn)其濃度值顯著高于其他站點(diǎn)，為明顯的異常值，由于這些極值對中心值和協(xié)方差矩陣存在影響，其他機(jī)構(gòu)非穩(wěn)健馬氏距離明顯小于穩(wěn)健馬氏距離，表明穩(wěn)健估計(jì)能夠排除極端值對于馬氏距離計(jì)算的影響。

穩(wěn)健馬氏距離越高表明該機(jī)構(gòu)所得結(jié)果偏離此次能力考核的中心值越遠(yuǎn)，可認(rèn)為其在此次能力考核中的表現(xiàn)越差(表4)。當(dāng)根據(jù)上文介紹的基于穩(wěn)健馬氏距離的算法判斷某監(jiān)測機(jī)構(gòu)能力考核結(jié)果為異常值時(shí)，說明由于某些分析環(huán)節(jié)中出現(xiàn)錯(cuò)誤，或質(zhì)量體系中存在的一些問題，導(dǎo)致其與其他合格站點(diǎn)數(shù)據(jù)分布不一致。通過與Z比分?jǐn)?shù)計(jì)算出的α-666、γ-666、p, p′-DDT和o, p′-DDT能力考核結(jié)果(滿意、有問題或不滿意)進(jìn)行比較(表4)，發(fā)現(xiàn)基于穩(wěn)健馬氏距離判定出來能力考核結(jié)果異常的監(jiān)測機(jī)構(gòu)其Z比分?jǐn)?shù)結(jié)果都存在著不滿意或有問題的項(xiàng)目，說明基于穩(wěn)健馬氏距離的異常值篩選方法能夠較好地從多元統(tǒng)計(jì)的角度發(fā)現(xiàn)能力考核中的不合格單位，其結(jié)果與Z比分?jǐn)?shù)方法結(jié)果較為一致。同時(shí)發(fā)現(xiàn)，穩(wěn)健馬氏距離較高的站點(diǎn)，其Z比分?jǐn)?shù)結(jié)果存在較多的不滿意或有問題，可以考慮在今后的能力考核工作中根據(jù)馬氏距離這一綜合指標(biāo)來對各個(gè)監(jiān)測機(jī)構(gòu)的整體分析能力做出評估。

4 結(jié)論

綜上所述，穩(wěn)健馬氏距離方法在兼具了馬氏距離的優(yōu)點(diǎn)的同時(shí)，又能較好地排除異常值對其的影響。與傳統(tǒng)的針對單輪次、單項(xiàng)目的穩(wěn)健統(tǒng)計(jì)方法(如Z比分?jǐn)?shù)方法、迭代法)相比，該方法可通過馬氏距離這個(gè)單一指標(biāo)根據(jù)多輪次、多項(xiàng)目的能力考核結(jié)果從穩(wěn)健多元統(tǒng)計(jì)的角度對實(shí)驗(yàn)室能力或數(shù)據(jù)質(zhì)量進(jìn)行綜合定量評價(jià)，并能有效篩選出結(jié)果存在問題的監(jiān)測機(jī)構(gòu)，可以在針對多項(xiàng)目、多輪次能力考核結(jié)果的綜合評價(jià)工作中加以試用。該方法在實(shí)際工作中也存在著一定的局限性，如某些機(jī)構(gòu)在年內(nèi)只參加了個(gè)別輪次的能力考核，或在某輪次的能力考核中只進(jìn)行部分樣品的分析測試，穩(wěn)健馬氏距離無法對這些機(jī)構(gòu)該年或該輪次的能力考核結(jié)果進(jìn)行綜合評估。

[1] 滕曼, 付強(qiáng), 楊婧, 等. 2011年全國環(huán)境監(jiān)測實(shí)驗(yàn)室地表水揮發(fā)性有機(jī)物檢測能力分析[J]. 環(huán)境與健康雜志,2013,30(12):1 108-1 109.

TENG M, FU Q, YANG J, et al. Results analysis of proficiency assessment of VOCs monitoring in water [J]. Journal of Environment and Health,2013,30(12):1 108-1 109.

[2] 滕曼, 付強(qiáng), 吳曉鳳, 等. 環(huán)境監(jiān)測實(shí)驗(yàn)室水中砷、汞監(jiān)測能力考核結(jié)果評價(jià)[J]. 中國環(huán)境監(jiān)測,2014,30(4):183-187.

TENG M, FU Q, WU X F, et al. Results analysis of proficiency assessment of As and Hg monitoring in ground water [J]. Environmental Monitoring of China,2014,30(4):183-187.

[3] 刑小茹, 馬小爽, 田文,等. 實(shí)驗(yàn)室間比對能力驗(yàn)證中的兩種穩(wěn)健統(tǒng)計(jì)技術(shù)探討[J]. 中國環(huán)境監(jiān)測,2011,27(4):4-8.

XING X R, MA X S, TIAN W, et al. Two robust statistic techniques in proficiency testing by interlaboratory comparisons [J]. Environmental Monitoring of China,2011,27(4):4-8.

[4] 吳忠祥. 實(shí)驗(yàn)室能力驗(yàn)證中的分割水平檢測樣品與穩(wěn)健統(tǒng)計(jì)技術(shù)[J]. 中國環(huán)境監(jiān)測,2003,19(4):8-10.

WU Z X. Split-level test sample and robust statistics techniques in laboratory proficiency testing [J]. Environmental Monitoring of China, 2003, 19(4): 8-10.

[5] MAHALANOBIS P C. On the generalized distance in statistics [J]. Proceedings of the National Institute of Sciences (Calcutta),1936(2):49-55.

[6] NORSHAHIDA S, ABDUL A J, MOHDT L, et al.

Anomaly detection and sssessment of PM10functional data at several locations in the Klang Valley, Malaysia [J]. Atmospheric Pollution Research,2015(6):365-375.

[7] ROUSSEEUW P J, VAN D K. A fast algorithm for the minimum covariance determinant estimator [J]. Technometrics,1999(41):212-223.

[8] ROUSSEEUW P J, CROUX C, TODOROV V, et al. Robustbase: basic robust statistics [M].2015.

[9] TODOROV V, FILZMOSER P. An object-oriented framework for robust multivariate analysis[J]. Journal of Statistical Software,2009,32(3):1-47.

[10] PISON G, VAN A S, WILLEMS G. Small sample corrections for LTS and MCD [J]. Metrika, 2002(55):111-123.

[11] HUBERT M, ROUSSEEUW P J, VERDONCK T. A deterministic algorithm for robust location and scatter [J].Journal of Computational and Graphical Statistics,2012(21):618-637.

[12] 王斌會, 陳一非. 基于穩(wěn)健馬氏距離的多元異常值檢測[J].統(tǒng)計(jì)與決策,2005(3):4-6.

WANG B H, CHEN Y F. Multivariate outlier detection based on the robust Mahalanobis distance [J]. Statistics & Decision,2005(3):4-6.

[13] FILZMOSER P, GARRETT R G, REIMANN C. Multivariate outlier detection in exploration geochemistry [J]. Computers & Geosciences,2005(31):579-587.

[14] GERVINI D. A robust and efficient adaptive reweighted estimator of multivariate location and scatter [J]. Journal of Multivariate Analysis,2003(84):116-144.

The Application of Robust Mahalanobis Distance in Proficiency Testing by Interlaboratory Comparisions

SHI Yaolong1,TENG Man1,LI Cheng2,WU Xiaofeng1,CHAI Wenxuan1,YANG Jing1,CHU Baolin1,FU Qiang1

1.State Environmental Protection Key Laboratory of Quality Control in Environmental Monitoring,China National Environmental Monitoring Centre,Bejing 100012,China2.Institute of Science and Technology,Agricultural University of Hebei,Baoding 071001,China

The result of proficiency testing of 4 organo-chlorine pesticide species among 96 enviromental monitoring labs was evaluated by the robust mahalanobis, and 25 labs were identified as outliers by this multi-statistics method. By the comparision between robust mahalanobis distance and Z-score, the outliers identified by robust mahalanobis distance were also indentified as outliers by Z-score. In addition, the fundamental and R implementation of outliers detection by robust mahalanobis distance were described in this article.

robust multi-statistics;robust Mahalanobis distance;interlaboratory comparisions;proficiency testing

2016-03-22;

2016-05-05

國家環(huán)保公益性行業(yè)科研專項(xiàng)“國家環(huán)境監(jiān)測網(wǎng)環(huán)境空氣自動(dòng)監(jiān)測(PM2.5、O3) 質(zhì)量保證與質(zhì)量控制技術(shù)體系研究與示范” (201409011)

師耀龍(1988-)，男，河北保定人，博士，工程師。

楚寶臨

X830.3

A

1002-6002(2017)02- 0127- 05

10.19316/j.issn.1002-6002.2017.02.20