基于非參數(shù)統(tǒng)計(jì)量的多元控制圖研究

郭佳晟，劉以建

(上海海事大學(xué)物流工程學(xué)院，上海 201306)

0 引言

統(tǒng)計(jì)過(guò)程控制(Statistical Process Control, SPC)作為經(jīng)典質(zhì)量控制方法，可以有效地提升產(chǎn)品質(zhì)量，保障生產(chǎn)過(guò)程。但其針對(duì)于單變量的局限性，對(duì)現(xiàn)代復(fù)雜工藝生產(chǎn)過(guò)程的應(yīng)用效果不佳，生產(chǎn)過(guò)程中，常常存在多個(gè)具有相關(guān)關(guān)系的質(zhì)量特性和過(guò)程參數(shù)，例如零件加工的長(zhǎng)度和直徑、化工過(guò)程的溫度、壓力等。只對(duì)單個(gè)變量監(jiān)控而不考慮變量之間的相關(guān)性會(huì)導(dǎo)致誤報(bào)警率顯著增加，因此需要用多元控制圖進(jìn)行過(guò)程監(jiān)控。Hotelling在1947年首先提出了基于T2統(tǒng)計(jì)量的多元控制圖，用于對(duì)包含多個(gè)質(zhì)量特性的生產(chǎn)過(guò)程實(shí)施統(tǒng)計(jì)監(jiān)控，由此有了多變量統(tǒng)計(jì)控制過(guò)程(Multivariate Statistical Process Control, MSPC)的研究。相繼有了多元累積和(Multivariate Cumulative Sum, MCUSUM)控制圖以及多元指數(shù)加權(quán)移動(dòng)平均(Multivariate Exponentially Weighted Moving Average, MEWMA)控制圖等。

傳統(tǒng)的控制圖可以稱為參數(shù)控制圖，即需要對(duì)總體分布有簡(jiǎn)單假定，例如正態(tài)分布。但實(shí)際所采集到的信息，可能無(wú)法對(duì)總體分布作出簡(jiǎn)單的假設(shè)，例如單邊尺寸線跳動(dòng)，呈現(xiàn)為指數(shù)分布或者 Weibull分布[1]。當(dāng)實(shí)際分布與假定分布有較大偏差時(shí)，基于分布的參數(shù)監(jiān)控會(huì)受到很大影響，監(jiān)控效果大大下降。針對(duì)上述缺陷的解決方法有數(shù)據(jù)變換方法和使用非參數(shù)統(tǒng)計(jì)量?jī)煞N思路，即將過(guò)程數(shù)據(jù)通過(guò)映射函數(shù)轉(zhuǎn)換成符合正態(tài)分布的形式，或使用不依賴于分布的統(tǒng)計(jì)量進(jìn)行監(jiān)控。用非參數(shù)統(tǒng)計(jì)量構(gòu)建的控制圖稱為非參數(shù)控制圖。

在過(guò)去幾年里,非參數(shù)控制圖已經(jīng)引起了很多關(guān)注，文獻(xiàn)[2]針對(duì)再制造過(guò)程的復(fù)雜特性，提出基于Wilcoxon統(tǒng)計(jì)量的EWMA控制圖；文獻(xiàn)[3-4]提出次序秩的非參數(shù)EWMA聯(lián)合控制圖和基于馬爾可夫均值估計(jì)量的自適應(yīng)CUSUM控制圖；文獻(xiàn)[5]采用基于在Logistic分布下尺度參數(shù)的漸近局部最優(yōu)勢(shì)檢驗(yàn)作為統(tǒng)計(jì)量，構(gòu)建了LOG控制圖；文獻(xiàn)[6]考慮同時(shí)對(duì)分布均值和標(biāo)準(zhǔn)差的監(jiān)控，提出結(jié)合Wilcoxon秩和檢驗(yàn)和Ansari-Bradley檢驗(yàn)的非參數(shù)控制圖；文獻(xiàn)[7]使用最小二乘支持向量機(jī)(Least Squares Support Vector Machine，LSSVM)所得到的概率值作為統(tǒng)計(jì)量，提出基于LSSVM的多元非參數(shù)控制圖；文獻(xiàn)[8]結(jié)合漢密爾頓路徑和游程檢驗(yàn)，提出了基于游程檢驗(yàn)的多元非參數(shù)控制圖；文獻(xiàn)[9]結(jié)合多元符號(hào)檢驗(yàn)，提出了非參數(shù)EWMA控制圖；文獻(xiàn)[10]將多元擬合優(yōu)度檢驗(yàn)與最小生成樹結(jié)合，設(shè)計(jì)了SMMST控制圖。

本文提出了一種結(jié)合Wilcoxon秩和檢驗(yàn)和Ansari-Bradley檢驗(yàn)的非參數(shù)控制圖，使用協(xié)方差矩陣構(gòu)建統(tǒng)計(jì)量，實(shí)現(xiàn)對(duì)多變量過(guò)程的監(jiān)控。

1 多元非參數(shù)控制圖

常用的參數(shù)檢驗(yàn)需要對(duì)總體分布有一定的估計(jì)，在對(duì)分布有假設(shè)的基礎(chǔ)上進(jìn)行檢驗(yàn)分析。但在實(shí)際數(shù)據(jù)分析過(guò)程中，可能無(wú)法對(duì)總體分布作簡(jiǎn)單假定，非參數(shù)檢驗(yàn)是不涉及總體分布的參數(shù)的檢驗(yàn)方法，常用秩和作為檢驗(yàn)統(tǒng)計(jì)量。

通過(guò)檢驗(yàn)可以判斷兩分布的參數(shù)是否相同，用于描述分布的參數(shù)有位置參數(shù)、尺度參數(shù)、形狀參數(shù)等。

位置參數(shù)是描述分布集中趨勢(shì)的度量，例如均值和中位數(shù)。常用的對(duì)位置參數(shù)的非參數(shù)檢驗(yàn)有Wilcoxon秩和檢驗(yàn)、游程檢驗(yàn)等。

尺度參數(shù)是描述分布分散程度的參數(shù)。常用的對(duì)尺度參數(shù)的非參數(shù)檢驗(yàn)有Ansari-Bradley檢驗(yàn)、Levene檢驗(yàn)等。

1.1 Wilcoxon秩和檢驗(yàn)

假設(shè)有樣本集X：

X=[X1,X2,…,Xi,…,Xn]

將X排序后得到Xr：

Xr=[X(1),X(2),…,X(k),…,X(n)]

其中，下標(biāo)(k)表示Xi在X中的次序秩，即Ri=k，Ri為樣本X的秩統(tǒng)計(jì)量。Wilcoxon秩和檢驗(yàn)是基于秩統(tǒng)計(jì)量的檢驗(yàn)方法，其構(gòu)建過(guò)程如下：

假設(shè)有樣本集X和Y：

(1)

其中，X服從分布F1(μ1,σ1)，Y服從分布F2(μ2,σ2)。

原假設(shè)為H0:μ1=μ2，備擇假設(shè)為H1:μ1≠μ2，定義Wilcoxon秩和檢驗(yàn)統(tǒng)計(jì)量[2]為：

(2)

其中，N=m+n。統(tǒng)計(jì)量Z的均值和方差[2]為：

1.2 Ansari-Bradley檢驗(yàn)

假設(shè)有樣本集X和Y，如式(1)所示。

原假設(shè)為H0:σ1=σ2，備擇假設(shè)為H1:σ1≠σ2，定義Ansari-Bradley檢驗(yàn)統(tǒng)計(jì)量[6]為：

(3)

其中，N=m+n。統(tǒng)計(jì)量T的均值和方差[6]為：

1.3 設(shè)計(jì)WAB控制圖

假設(shè)有受控?cái)?shù)據(jù)樣本集X0和質(zhì)量樣本X

(4)

其中，X0～F1(μ1,σ1)，X～F2(μ2,σ2)，Xi是d維向量。對(duì)質(zhì)量樣本X增加擾動(dòng)δ，觀察統(tǒng)計(jì)量Z和T的變化情況。取d=2，δ～N(0,1)，混合樣本集的統(tǒng)計(jì)量Z和統(tǒng)計(jì)量T的變化情況如圖1、圖2所示。

圖1 擾動(dòng)為δ時(shí)統(tǒng)計(jì)量的折線圖

圖2 擾動(dòng)為1.5δ時(shí)統(tǒng)計(jì)量的折線圖

觀察圖1、圖2可以發(fā)現(xiàn)，當(dāng)X存在擾動(dòng)時(shí)，統(tǒng)計(jì)量Z和統(tǒng)計(jì)量T的值都有或大或小的漂移。使用樣本協(xié)方差矩陣將統(tǒng)計(jì)量Z和統(tǒng)計(jì)量T結(jié)合，使d維向量轉(zhuǎn)換為單個(gè)值[11]，將其作為統(tǒng)計(jì)量構(gòu)建多元非參數(shù)控制圖，即定義統(tǒng)計(jì)量P為：

(5)

其中，Z是Wilcoxon秩和統(tǒng)計(jì)量向量，S是協(xié)方差矩陣，T是Ansari-Bradley統(tǒng)計(jì)量向量。

使用統(tǒng)計(jì)量P構(gòu)建非參數(shù)控制圖，記為WAB控制圖。WAB控制圖同時(shí)使用了標(biāo)準(zhǔn)化后的Z統(tǒng)計(jì)量和T統(tǒng)計(jì)量，放大了漂移，實(shí)現(xiàn)對(duì)過(guò)程的監(jiān)控。設(shè)置控制限h，當(dāng)P>h時(shí)，說(shuō)明質(zhì)量樣本的分布發(fā)生改變，即生產(chǎn)過(guò)程失控，發(fā)出報(bào)警?？刂葡辢的詳情將在下面部分給出。

1.4 控制限判定

評(píng)價(jià)控制圖的指標(biāo)常用平均運(yùn)行長(zhǎng)度(Average Run Length, ARL)，ARL分為受控ARL(ARL0)和失控ARL(ARL1)[12]。ARL0指受控狀態(tài)下控制圖第一個(gè)虛報(bào)樣本之前的平均樣本數(shù)；ARL1指失控狀態(tài)下控制圖發(fā)現(xiàn)第一個(gè)失控樣本之前的平均樣本數(shù)。本文用蒙特卡羅模擬方法[13]獲取控制限h的值，在確定h值時(shí)，一般設(shè)ARL0=200，計(jì)算步驟為：

(1)選取1組受控?cái)?shù)據(jù)樣本集X0。

(2)給定控制限h的初始值。

(3)生成1000組與X0同分布的樣本集X。

(4)計(jì)算運(yùn)行鏈長(zhǎng)RL(從第一個(gè)移動(dòng)窗口到控制圖首次報(bào)警的移動(dòng)窗口個(gè)數(shù))

(5)重復(fù)過(guò)程(3) ～(5)K次。

(6)計(jì)算K個(gè)RL值的均值。

(7)如果ARL0接近設(shè)定值，則h即為選定的控制限。反之則增大或減小h值并重復(fù)上述過(guò)程。

2 性能分析

使用本文提出的WAB控制圖與T2控制圖進(jìn)行過(guò)程失控狀態(tài)下的性能比較，性能指標(biāo)使用ARL1。樣本數(shù)據(jù)為：

(6)

其中，x·1服從正態(tài)分布，x·2服從Weibull分布。模擬真實(shí)生產(chǎn)過(guò)程先受控后失控的情況，在Xi,i>500處開始加入擾動(dòng)δ，δ～N(0,1)。依據(jù)上文的方法確定控制限h=0.824，選取(m,n)組合為(50,10)。對(duì)不同程度的漂移的ARL1如表1所示。控制圖的虛報(bào)情況如表2所示。

表1 不同δ下ARL1的值

從表1可以看出，WAB控制圖有比T2控制圖更高的敏感度，能更快的發(fā)現(xiàn)失控情況。在小漂移的情況有明顯有優(yōu)勢(shì)，隨著漂移量的增大，逐漸與T2控制圖性能持平。

從表2可以看出WAB控制圖基本沒(méi)有虛報(bào)的現(xiàn)象，比T2控制圖有更高的穩(wěn)定性。

表2 不同δ下誤報(bào)樣本的數(shù)量

3 應(yīng)用實(shí)例

以某工廠汽車發(fā)動(dòng)機(jī)缸蓋生產(chǎn)線為例，對(duì)提出的WAB控制圖進(jìn)行應(yīng)用和驗(yàn)證。選取6個(gè)質(zhì)量參數(shù)進(jìn)行監(jiān)控，采取共100組樣本。其樣本數(shù)據(jù)的數(shù)據(jù)特征，如表3所示。部分質(zhì)量參數(shù)的直方圖如圖3所示，可以看到存在明顯不符合正態(tài)分布的參數(shù)。

表3 質(zhì)量參數(shù)取值范圍表

圖3 部分質(zhì)量參數(shù)直方圖

根據(jù)前文論述的方法，采用WAB控制圖對(duì)樣本數(shù)據(jù)進(jìn)行監(jiān)控，結(jié)果如圖4所示，在第60個(gè)樣本處統(tǒng)計(jì)量超出控制限，觸發(fā)報(bào)警，而樣本實(shí)際在57處存在失控。由此可以表明WAB控制圖的有效性。

圖4 WAB控制圖

4 結(jié)論

當(dāng)樣本數(shù)據(jù)不服從多元正態(tài)分布時(shí)，傳統(tǒng)基于多元正態(tài)分布的控制圖實(shí)用性不佳，針對(duì)這一缺陷，本文提出了WAB控制圖。使用非參數(shù)檢驗(yàn)的方法，使統(tǒng)計(jì)量不依賴于過(guò)程分布，同時(shí)監(jiān)控位置參數(shù)和尺度參數(shù)，通過(guò)協(xié)方差矩陣構(gòu)建多元非參數(shù)控制圖的統(tǒng)計(jì)量P。使用蒙特卡羅方法對(duì)控制圖進(jìn)行性能分析，與T2控制圖進(jìn)行比對(duì)。結(jié)果顯示，WAB控制圖能明顯減少由于非正態(tài)分布所導(dǎo)致的誤報(bào)現(xiàn)象，并且可以更快的發(fā)現(xiàn)分布漂移，有著更高的靈敏度和準(zhǔn)確性。最后，通過(guò)汽車發(fā)動(dòng)機(jī)缸蓋生產(chǎn)過(guò)程中6個(gè)質(zhì)量參數(shù)的測(cè)試數(shù)據(jù)對(duì)WAB控制圖的性能進(jìn)行了驗(yàn)證，結(jié)果表明當(dāng)過(guò)程偏移時(shí)，WAB控制圖能夠有效檢測(cè)出過(guò)程偏移。