張 陽，秦幸幸

(天津商業(yè)大學 1.管理學院；2.管理創(chuàng)新與評價研究中心，天津 300134)

在統(tǒng)計過程控制中，產(chǎn)品或過程的質(zhì)量特性可以通過某種函數(shù)關系來描述，這種函數(shù)關系被稱為輪廓(profile)[1]。輪廓控制(profile monitoring)即基于統(tǒng)計過程控制理論對質(zhì)量特性表現(xiàn)為輪廓的產(chǎn)品或過程質(zhì)量進行監(jiān)控，通常包括兩個階段。第1階段基于歷史輪廓數(shù)據(jù)分析來判斷過程的穩(wěn)定狀況，并識別受控輪廓集和確定輪廓控制的參數(shù)，可稱為輪廓控制參數(shù)識別階段；第2階段則基于所確定的控制參數(shù)，建立控制圖，實施過程監(jiān)控，稱為輪廓控制監(jiān)控階段。輪廓控制參數(shù)識別的關鍵在于受控輪廓集的識別，它是確定輪廓控制參數(shù)的基礎。因此，本文將研究重點聚焦于受控輪廓集的識別。傳統(tǒng)的輪廓控制參數(shù)識別方法大多假設輪廓內(nèi)測量數(shù)據(jù)之間相互獨立[2-3]。然而，在實際情況中，因連續(xù)測量時間短、數(shù)據(jù)量大等原因，會使輪廓內(nèi)測量點間存在相關關系。此外，實際生產(chǎn)加工過程中存在的系統(tǒng)性誤差，也會導致輪廓內(nèi)測量數(shù)據(jù)存在相關性[4-7]。因此，忽略輪廓內(nèi)部相關性的輪廓控制方法已不再適用。

線性混合模型通常用來對內(nèi)部存在相關性的輪廓進行建模[8-9]，并在此基礎上構建輪廓控制參數(shù)識別方法[10-13]。該類方法首先采用線性混合模型進行輪廓建模，并估計樣本輪廓參數(shù)。然后，基于輪廓參數(shù)建立T2控制圖，用于識別受控輪廓集，進而確定輪廓控制參數(shù)。構建T2控制圖時，可采用連續(xù)差分估計量(SDE)[8]、最佳線性無偏估計量(eblups)[10]和基于最小體積橢球(MVE)的估計[11-12]等方法估計T2統(tǒng)計量的方差?協(xié)方差矩陣。然而，基于上述方法的T2控制方法是基于全部樣本輪廓數(shù)據(jù)來估計方差?協(xié)方差矩陣，進而確定控制限。這將會使得T2控制圖的識別性能因異常輪廓的存在而降低，尤其是基于最佳線性無偏估計量的方法受到的影響極大。

為降低因樣本輪廓中存在異常輪廓所帶來的影響，Chen等[13]提出基于層次聚類(cluster based，CB)算法的受控輪廓識別方法，并在此基礎上建立T2控制圖，一定程度上降低異常輪廓的影響，且與其他方法相比基于層次聚類的T2控制圖的性能更好。然而，此方法需要首先選定一半以上的樣本輪廓作為初始受控輪廓，但其中往往會存在較多異常輪廓。由于聚類結果依賴于距離度量的選取，異常輪廓的存在，尤其是較少離群值的存在，會影響最終聚類結果，進而影響最終確定的受控輪廓集。這會使得基于層次聚類算法的T2控制圖識別性能會降低。綜上所述，降低離群值和異常輪廓對輪廓控制參數(shù)識別階段的影響，是構建穩(wěn)健輪廓控制參數(shù)識別方法的重點和難點。

為了更準確地從樣本輪廓數(shù)據(jù)中識別出受控輪廓集，本文提出基于密度(density based，DB)的輪廓控制參數(shù)識別方法?；诿芏鹊乃惴軌蛴行У乜朔哟尉垲愃惴ǖ娜毕?，為構建穩(wěn)健的T2控制圖提供保證。所提方法使用樣本輪廓參數(shù)向量計算樣本輪廓的密度；再基于樣本輪廓的密度確定一部分初始受控輪廓，排除異常輪廓，尤其是離群值的干擾；隨后經(jīng)過多次迭代獲得最終受控輪廓集；在此基礎上采用參數(shù)均值方法確定輪廓控制參數(shù)?；诿芏鹊妮喞刂茀?shù)識別方法得到的輪廓控制參數(shù)，可用于構建輪廓控制圖。

1 模型假設

在生產(chǎn)加工過程中，假定已獲取m條樣本輪廓數(shù)據(jù)，且第i條樣本輪廓的觀測數(shù)據(jù)為(xit,yit) ，i=1,2, ···,m；t=1, 2, ···,n。本文假設在不同輪廓中測量點位置固定，即當t不變時，xit=xt?；贘ensen等[8]所使用的線性混合模型，建立輪廓模型為

其中，Yi=(yi1,yi2, ···,yin)T為 輪廓測量值；Xi是n×o階向量，表示與固定效應相關的回歸變量，o為回歸變量數(shù)目； β為o×1階向量，表示輪廓固定效應參數(shù)，且在不同的輪廓中 β不 變；Zi為n×s階矩陣；bi為s×1階向量，表示輪廓隨機效應參數(shù)，且bi服從均值為0、協(xié)方差矩陣為G的多元正態(tài)分布 (MN)，記為bi～ MN(0,G)，G為s×s階 矩陣； εi為隨機誤差，且εi～MN(0,Ri)，Ri為n×n階矩陣[14]。

Ruppert等[15]給出模型(1)的矩陣表達式為

在實際情況中，V在一般情況下是未知的，常用最大似然估計或最小二乘法來得到B和R[15]，從而得出V。在3.1節(jié)仿真設置中，假設V已知。

基于式(3)和(4)可估計輪廓參數(shù)向量為

其中，b?i為 矩陣b?中的第i行元素所組成的向量。

2 基于密度的參數(shù)識別方法

2.1 密度概念

本文將數(shù)據(jù)向量密度定義為以空間內(nèi)某個數(shù)據(jù)向量為中心的最小半徑的倒數(shù)，且此最小半徑需包含周圍P個數(shù)據(jù)向量[16]。其中，P表示密度參數(shù)，可根據(jù)具體情況來確定。不難看出，根據(jù)密度大小可以確定該數(shù)據(jù)向量是離群點的可能性。如果某個數(shù)據(jù)向量的密度越小，其離周圍P個數(shù)據(jù)向量的最大距離就越大，則該數(shù)據(jù)向量為離群點的可能性就越大。具體來講，假設空間內(nèi)有一組數(shù)據(jù)向量集Uω={γ1, γ2, ···, γω}， 其中，γl=(γ1l, γ2l, ···, γφl)，l=1,2, ···, ω ， 為 φ維向量。對于任意數(shù)據(jù)向量，分別計算γl與其最近的P個數(shù)據(jù)向量的距離Dp,l={d1,l,d2,l, ···,dp,l}。 則數(shù)據(jù)向量γl的 密度ρl=1/Dmax,l，其中，Dmax,l=max{d1,l,d2,l, ···,dp,l}。數(shù)據(jù)向量密度的矩陣式算法如下。

2.2 基于密度的輪廓控制參數(shù)識別方法

基于密度的輪廓控制參數(shù)識別方法包括建立輪廓模型、確定初始受控輪廓集、識別受控輪廓集和確定輪廓控制參數(shù)等4個階段。其中，初始受控輪廓確定階段是整個輪廓控制參數(shù)識別方法的重點，是根據(jù)輪廓模型參數(shù)向量的密度來選定初始受控輪廓。與基于層次聚類的方法相比，基于密度的方法在確定初始受控輪廓時可以事先排除離群輪廓。因此，基于密度的方法可以有效降低離群輪廓對控制參數(shù)識別的影響?；诿芏鹊妮喞刂茀?shù)識別流程如圖1所示，具體步驟如下。

圖1 基于密度的輪廓參數(shù)識別流程圖Figure 1 The flow chart of profile parameters identification method based on data-density

階段1 建立輪廓模型。

基于式(2)對m條樣本輪廓建模，然后可得到輪廓參數(shù)向量的估計 β?i。記C={β?1, β?2, ···, β?m}，表示所有樣本輪廓集合。

階段2 確定初始受控輪廓集。

本階段將采用輪廓密度來選定初始受控輪廓。首先，基于2.1節(jié)中密度的矩陣式算法，計算并獲得參數(shù)向量集C的密度矩陣DC。 之后，選取DC中第P列的m個元素，并將其按照從大到小的順序排列。最后，選取元素排列中前Q個元素所對應的輪廓作為初始受控輪廓集C0。此處，參數(shù)Q和P的取值影響所提基于密度的輪廓控制參數(shù)識別的性能。參數(shù)Q和P設置對所提方法的影響將在仿真分析部分討論。

階段3 受控輪廓集識別。

1) 設定受控輪廓集為Cq，并先設定其初始值為初始受控輪廓集，即Cq=C0。

2) 基于Cq中Q條輪廓數(shù)據(jù)，利用式(3) ～ (5)更新Cq中每個輪廓的參數(shù)，并計算均值 βˉ。 記Cm?q=C?Cq為不包含在Cq中 的輪廓集，并計算其中第 τ個輪廓的T2統(tǒng)計量為

4)如果r= 0，那么Cq則為基于密度方法得到的最終受控輪廓集 CI， 且 CI=Cq。此時可以終止程序。如果r≠ 0，則令q=q+r，Cq=Cq+r，并返至2)，繼續(xù)受控輪廓集的識別。

階段4 確定輪廓控制參數(shù)。

利用所確定的受控輪廓集 CI，根據(jù)線性混合模型使用極大似然估計或最小二乘法，重新估計每條輪廓的參數(shù)，并求輪廓參數(shù)的均值，即可得輪廓控制參數(shù)。

3 仿真分析

輪廓控制參數(shù)識別方法的關鍵是正確地識別出受控輪廓集。仿真分析主要通過蒙特卡洛模擬來分析正確識別受控輪廓集的能力，并與已有方法進行比較分析。

3.1 仿真設置

3.2 識別性能指標確定

輪廓控制參數(shù)識別方法正確識別受控輪廓集的能力，可以通過分類精度、敏感性和特異性等性能指標來分析。基于表1的混淆矩陣可根據(jù)表2的公式來計算分類精度、敏感性和特異性等性能指標[18]。其中，a為表示正確識別異常輪廓的數(shù)目；e為把異常輪廓識別為受控輪廓的數(shù)目；f為把受控輪廓識別為異常輪廓的數(shù)目；h為正確識別受控輪廓的數(shù)目。容易得知，分類精度、敏感性和特異性的取值范圍均為 [ 0, 1]。而且，當這些指標的值越大時，則表明所提方法正確識別受控輪廓集的能力越好。

表1 混淆矩陣Table 1 Confusion matrix

表2 性能指標公式Table 2 Performance index formula

然而，特異性和敏感性呈現(xiàn)的負相關關系，即敏感性取值變大(小)時特異性取值則變小(大)，會為同時采用敏感性和特異性比較不同方法性能時帶來一定困難。因此，為方便比較，可綜合考慮敏感性和特異性指標，采用約登指數(shù)(Youden index)[19-22]來分析受控輪廓集正確識別的能力，其定義為

同樣，約登指數(shù)取值范圍也為 [ 0, 1]，且取值越大時識別性能越好。本文將基于分類精度和約登指數(shù)2個指標來分析參數(shù)取值和置信水平的影響。

3.3 參數(shù)P和Q的影響分析

在分析參數(shù)P和Q對識別性能的影響時，基于3.1節(jié)中輪廓模型，首先通過多次蒙特卡洛模擬仿真初步分析參數(shù)P和Q的取值及其對應性能指標取值。根據(jù)初步分析結果得知，當P和Q的值介于區(qū)間[[m/3]+1, [m/2]+1]中時，性能指標相對較好。因此，為進一步分析所提方法的識別性能，將針對不同偏移量進行以下4種組合試驗，即

在蒙特卡洛仿真分析中，不同偏移量下各組合模擬仿真次數(shù)均為5 000次，其中，α設定為0.025。最終仿真結果如表3所示，其中加粗結果為特定偏移量下不同參數(shù)組合中分類精度和約登指數(shù)的最優(yōu)結果。

表3 不同參數(shù)組合下識別性能比較Table 3 Performance comparisons on different combinations of parameters P and Q

由表3可知，當P一定時，對于任意偏移量，采用較小的Q值可獲得分類精度和約登指數(shù)的較優(yōu)結果。另外，當Q一定時，對于小于等于0.25的偏移量，P較小時分類精度和約登指數(shù)結果較優(yōu)；而對于偏移量大于0.25的情況，P較大時可獲得較好的性能。此外，對于特定參數(shù)組合，隨著偏移量的不斷增加，分類精度和約登指數(shù)也在逐步變大；而且當偏移量大于等于0.25時，分類精度和約登指數(shù)的取值大都超過0.95，并接近于1。

通過比較參數(shù)P和Q的不同組合下分類精度和約登指數(shù)的值，可以發(fā)現(xiàn)，無論P取何值，當Q=[m/3]+1時總能獲得最優(yōu)識別性能。因此，在實際應用中，可以設定參數(shù)Q為 [m/3]+1，然后再根據(jù)偏移量大小設定參數(shù)P的取值。具體而言，當偏移量為中小偏移時，設定P=[m/3]+1且Q=[m/3]+1；當偏移量為大偏移時，設定P=[m/3]+1且Q=[m/3]+1。

然而，上述參數(shù)P和Q的設定需要假設偏移量根據(jù)以往經(jīng)驗預先已知。當偏移量未知時，建議選取參數(shù)組合C4，即設定參數(shù)P和Q分別為 [m/3]+1和[m/3]+1。因為，雖然對于大偏移的情況，組合C4下分類精度和約登指數(shù)的取值比組合C3下較低，但是也超過0.95，并逐漸接近于1。

3.4 置信水平α的影響分析

基于前述分析結果，本節(jié)將在參數(shù)P和Q分別取[m/3]+1和 [m/3]+1的情況下研究置信水平α對所提方法性能的影響。在不同α取值下，分類精度和約登指數(shù)隨偏移變化情況如圖2(a)和(b)所示。

由圖2可知，對于偏移量為中小偏移的情況，α= 0.05時，分類精度和約登指數(shù)均能最優(yōu)，除了分類精度在偏移量為0.05時取值最小之外。而當偏移量較大時，α= 0.01時性能較優(yōu)。因此，若偏移量已知時，可以為中小偏移給定相對較大的α值，此時識別準確性不高，較大的α值能夠有效控制第2類錯誤的發(fā)生；而對于較大偏移，可以設定相對較小的α值，來降低識別準確性很高時容易發(fā)生的第1類錯誤概率。

圖2 偏移量對分類精度和約登指數(shù)的影Figure 2 The effect of offset on classification accuracy and Yoden index

然而，當實際應用中偏移量未知時，建議α設定為0.05。盡管較大的α值會削弱偏移量較大時所提方法的識別性能，但是偏移量較大時分類精度和約登指數(shù)都非常接近于1，也能得到良好的分類效果。

3.5 性能比較研究

本節(jié)將比較分析所提基于密度的分析方法(DB)與已有方法的識別性能。此處選取P=[m/3]+1，Q=[m/3]+1， α=0.05。已有方法將選取基于層次聚類的分析方法(CB)[13]、基于最小體積橢球的分析方法(MVE)[11-12]和基于連續(xù)差分的分析方法(SDE)[8]。性能比較分析采用用蒙特卡洛模擬仿真，各種方法的性能指標結果均基于5 000次重復計算，并列于表4中，其中加粗結果為特定偏移量下不同方法中各性能指標中的最優(yōu)結果。

表4 不同分析方法的性能比較Table 4 Performance of different methods

由表4易知，根據(jù)約登指數(shù)可以得出所提基于密度的方法要優(yōu)于所有已有方法，而且根據(jù)分類精度也可以看出所提方法在大多數(shù)情況下也表現(xiàn)出最優(yōu)的識別效果。當偏移量為0.05或0.3時，根據(jù)分類精度的結果，雖然基于密度的方法表現(xiàn)較弱，但是，與最優(yōu)方法相比，在分類精度的結果上相差不大。綜合在分類精度與約登指數(shù)上的整體表現(xiàn)，所提基于密度的分析方法的綜合性能最好，表現(xiàn)出較優(yōu)的識別性能。

4 結論

目前已有的輪廓控制參數(shù)識別方法受異常輪廓的干擾較大。本文所提出的基于密度的受控輪廓識別方法，改進了已有方法在初始受控輪廓確定時因異常輪廓存在造成敏感性降低的不足之處，進而提高了估計輪廓參數(shù) βˉ的準確性，獲得更為穩(wěn)健的T2統(tǒng)計量，從而可以提高識別性能。另外，本文建議使用約登指數(shù)作為性能指標，可以綜合分析識別性能，從而在比較各種方法的性能時更加方便。此外，基于蒙特卡洛模擬，分析所提方法中初始受控輪廓數(shù)目Q、密度參數(shù)P和置信水平α對識別性能的影響，并給出當偏移量未知時Q、P和α的設定建議。最后，比較分析所提方法與已有方法的識別性能。結果顯示，綜合分類精度和約登指數(shù)的結果，所提基于密度的輪廓參數(shù)控制方法的性能要優(yōu)于其他方法。因此，基于密度的輪廓控制參數(shù)識別方法在很大程度上削弱了異常輪廓對T2控制圖的干擾，對于識別產(chǎn)品或過程中的受控輪廓集非常有效。

由于所提基于密度的方法依賴于初始受控輪廓數(shù)目Q和密度參數(shù)P，從而影響初始受控輪廓集的確定。因此，未來研究中，可以綜合運用其他方法，如K均值聚類等，與所提方法來改進初始受控輪廓集的確定。如首先基于不同方法各自確定初始受控輪廓，然后通過比較分析，取不同方法中都包含的初始受控輪廓作為初始受控輪廓集，可以進一步降低異常輪廓的干擾，從而提升參數(shù)估計精度，進而提高識別性能和輪廓控制參數(shù)的合理性。