非正態樣本下的多元控制圖的性能優化

郭佳晟，劉以建

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

0 引言

統計過程控制(Statistical Process Control, SPC)是經典質量控制方法，其對于單一變量的監控可以有效的提高生產質量。但在現代復雜工藝的生產過程的生產過程中，常常存在多個具有相關關系的質量特性和過程參數，例如零件加工的長度和直徑、化工過程的溫度、壓力等。Hotelling H在1947年首先提出了基于T2統計量的多元控制圖，用于對包含多個質量特性的生產過程實施統計監控，由此有了多變量統計控制過程(Multivariate Statistical Process Control, MSPC)的研究[1]。相繼有了多元累積和(Multivariate Cumulative Sum, MCUSUM)控制圖[2]以及多元指數加權移動平均(Multivariate Exponentially Weighted Moving Average, MEWMA)控制圖[3]等。

普通多元控制圖都是基于多元正態的假設上構建的。然而實際生產過程中，樣本數據存在非正態現象[4]。當過程變量較多時，過程數據也很難滿足多元正態假設，多元控制圖的效果會大大下降，導致誤報和漏報。因此針對非參數控制圖的研究是必要的，即不依賴總體分布，依舊能對生產過程達到監控效果的控制圖。

在過去幾年里,非參數控制圖已經引起了很多關注，有了許多關于非正態數據的質量控制和MSPC的研究，如喬新穎[5]提出基于偏度校正的非正態過程質量控制圖和基于比例加權方差方法的非正態過程質量控制圖；魏明月[6]提出對Box-Cox數據轉換函數的改進方法；姜興宇等[7]應用秩統計量分割數據與其分布的關系，結合滑動窗口實現了動態再制造過程質量的自適應監控；石文華等[8]考慮到軸承徑向間隙參數可能存在的非正態情況，用Johnson轉換的方法進行非正態過程能力分析；劉瀏等[9]使用基于次序秩的非參數控制圖，引入馬爾可夫均值估計量預測未知的漂移大小；裴德昭等[10]提出了基于游程檢驗的多元非參數控制圖。

本文提出了一種基于改進高斯混合模型的T2控制圖，使用高斯混合模型對樣本數據擬合，用多個服從正態的數據簇進行T2值計算，實現統計量與數據分布之間關系的解耦。針對高斯混合模型的初始化方法，提出數據密度的初始化方法。為了方便，將本文構建的控制圖記為G-T2控制圖。

1 HotellingT2控制圖

HotellingT2控制圖是基于T2統計量的多元控制圖。T2控制圖的構建過程如下：假設有n組樣本的樣本數據集X:

(1)

X服從p元正態分布Np(μ,Σ)，其中,μ為均值向量，Σ為協方差矩陣。T2統計量的定義為：

(2)

(3)

(4)

T2控制圖的上、下控制限(UCL、LCL)的計算公式為：

(5)

其中，Fα(p,n-p)是自由度為p和n-p的F分布在α分位點的值。當T2>UCL時，表明生產過程失控，控制圖發出警報。

2 高斯混合模型

聚類算法是數據挖掘時常用的處理手段，是非監督學習的一種經典算法[12]。常見的聚類算法有K-Means，DBSCAN，GMM等。

高斯混合模型(Gaussian Mixture Model, GMM)是一種常見的聚類算法，通過最大期望算法(Expectation-Maximization Algorithm, EM算法)將已知數據分解為多個高斯分布的疊加，對原始數據完成分類的操作。理論上，如果某個高斯混合模型包含的高斯模型個數足夠多，它們之間的權重設定得足夠合理，這個混合模型可以擬合任意分布的樣本。

2.1 高斯模型

高斯分布(單高斯模型)是最常見的分布，其概率密度函數為:

(6)

其中，X=(x1,x2,…,xn)是n維向量，Σ是X的協方差矩陣。

高斯混合模型是單高斯模型的拓展，假設有k個高斯模型，則混合高斯模型的概率密度函數為：

(7)

2.2 EM算法

2.2.1 似然函數

設有樣本集Y=y1,y2,…,yn，樣本y服從高斯分布。p(yn|μ,Σ)是表示變量Y=yn的概率。假設樣本之間符合獨立同分布(IID)，則樣本集Y的似然函數：

(8)

抽取任意樣本y的概率為：

(9)

在等式(8)中將權重ξi作為條件，得到更新后的似然函數[13-14]：

(10)

2.2.2 E-Step

對樣本集Y=y1,y2,…,yn引入隱變量γ，將非完全數據集Y轉變為完全數據集：

Z=(Y,γ)=(yi|γ1,γ2,…,γk)

(11)

當yi由第k個高斯模型采樣得來時，γ1=0,γ2=0,…,γk-1=0,γk=1。此時完全數據集的似然函數為：

(12)

其中,θ表示高斯分布的參數。

定義對數似然函數的期望為：

Q(θ,θm)=E[lnL(θ)|Y,θm]

(13)

其中,θ表示參數，θm表示已知參數估計值。

2.2.3 M-Step

最大化期望值Q(θ,θm)，即找到θmax，滿足：

θmax=argmax[Q(θ,θm)]

(14)

2.3 K-Means初始化

GMM的初始化值會影響算法的效果。不良的初始化值可能使算法陷入局部最優，且增加計算時間。若初始化值在真實聚類中心附近，則能大大地降低計算速度，且避免算法得出局部最優解。除了隨即初始化，最常見的初始化方法為K-Means初始化。

K-Means方法是一種帶參數的算法，需要預先確定K的值，即確定分成幾類。算法的流程圖如圖1所示。

圖1 K-Means算法流程

3 基于改進GMM的非參數T2控制圖

G-T2控制圖的監控流程圖如圖2所示。監控過程分為離線學習模型和在線監控兩部分。

(15)

計算Ni內數據點的T2值，根據每個樣本數據計算得出的T2值繪制T2控制圖，根據受控數據得出的T2值設定控制限的值。

圖2 非參數控制圖監控流程圖

在線監控過程中，對于新的樣本數據Xn+1，使用GMM模型預測其對應的子類Nk，使用μk和σk計算Xn+1的T2值，并判斷是否受控。

采用K-Means初始化的GMM算法需要根據指定類別個數，沒有良好的先驗知識可能導致效果很差。有學者提出了基于密度的聚類算法[15-16]，根據概率和人們的直觀感受，聚類中心應該有更多的數據點。在此基礎上，本文提出一種基于密度的非參數GMM初始化方法。

計算初始聚類中心采用廣度優先的搜索，即選取某一個數據點Xi=[x1x2…xp]為圓心，ε為半徑得到一個超球體Ci，Ci內包裹有其他數據點Xj為Ci所捕獲的點。每次對Ci內所有點判斷完畢后，再依次向外擴張獲取Cj，直至超球體不再包裹其他數據點。此時所選擇的所有點若大于設定的最小樣本數m，則判斷為一個類，該類中最密集的區域，即被包裹最多次的點為估計聚類中心。計算過程中判斷不成類的點則為離散點。

通常ε和m是根據經驗設定，本文使用統計特征中的十分位數，確定ε和m的值。避免先驗知識確定參數，實現非參數的目的。

當p=2時，基于密度的初始化方法的聚類中心選擇結果示意圖如圖3所示。初始化方法流程圖如圖4所示。

圖3 聚類中心選擇結果示意圖

圖4 基于密度的初始化方法流程圖

4 G-T2控制圖的數據仿真驗證

本文選用蒙特卡羅方法(Monte Carlo method)驗證G-T2控制圖的可行性。采用Weibull分布產生非正態數據。Weibull的概率密度函數為：

其中，λ>0是比例參數，k>0是形狀參數。

生成200組樣本測試數據，每組樣本包含兩個特征參數，即樣本集為：

其中，x·1服從正態分布，x·2服從Weibull分布，數據整體不服從二元正態分布。

4.1 無擾動測試

先使用不加入擾動的數據，即受控數據進行測試。由于隨機生成的測試數據服從一定的分布，所以只存在隨機因素對數據產生影響，即由這組數據所模擬的生產過程是處于統計受控狀態[3]。

樣本數據的直方圖如圖5所示，從圖中可以看到，x·1服從正態分布，x·2不服從正態分布。

圖5 數據直方圖

通過多次測試，發現混合權重系數α=0.8時為最佳。使用G-T2值繪制的控制圖與T2值繪制的控制圖進行比對，如圖6所示，可以發現改進后的T2值曲線更平滑，離散程度也更低。T2值與G-T2值的分布情況如表1所示。

圖6 T2控制圖與G-T2控制圖的比較

表1 統計量值分布情況

4.2 帶擾動測試

對樣本數據加入擾動，測試控制圖的靈敏性。在Xi,i>100處開始加入擾動δ，δ～N(0,1)。數據直方圖如圖7所示，控制圖如圖8所示。

圖7 數據直方圖

圖8 T2控制圖與G-T2控制圖的比較

在圖8中可以看到，相比于T2控制圖，G-T2控制圖異常值更加明顯，大約在第100個樣本處出現了異常值。T2控制圖整體受控值和異常值區別不大，容易造成對生產過程的誤判。

4.3 重復隨機試驗測試

評價控制圖的性能常用的指標是平均運行鏈長(ARL)，其分為受控平均鏈長ARL(0)和失控平均鏈長ARL(1)。此處使用ARL(1)對控制圖的進行測試。

同4.2節，對樣本數據加入擾動，對擾動發現的越早即表示控制圖的性能越好。控制限的值由受控狀態下的統計量最大值決定。

本節對文中所提方法進行了1000次重復隨機試驗，計算平均值以獲取穩定的結果。結果如表2所示。

表2 重復試驗結果

可以從表中看到，G-T2控制圖的UCL比T2控制圖低，這表明了在受控狀態下G-T2控制圖中統計量值的離散程度更低，曲線更加平滑；G-T2控制圖的ARL(1)比T2控制圖低，這表明了當出現失控情況時，G-T2控制圖能更快地在發現過程失控。

5 結論

當樣本數據不服從多元正態分布時，傳統基于多元正態分布的控制圖實用性不佳，針對這一缺陷，本文提出了G-T2控制圖。使用GMM方法對樣本數據進行擬合，采取數據變換的思路，將非正態轉化為多元正態分布，并對GMM算法進行了優化。通過蒙特卡羅方法進行數據仿真實驗，分別對T2控制圖和G-T2控制圖進行性能評估。結果表明：當過程受控時，G-T2控制圖能明顯減少由于非正態特性所導致的誤報現象，較T2控制圖有更高的監控準確率；當過程出現失控時，G-T2控制圖能及時發現失控現象，較T2控制圖有更小的失控鏈長。因此，本文提出的G-T2控制圖可以在過程分布未知時，對生產過程依然保持良好的監控效果。