非正態(tài)變異下的非線性輪廓異常點(diǎn)識(shí)別方法研究

聶 斌， 王 曦， 胡 雪

(天津大學(xué) 管理與經(jīng)濟(jì)學(xué)部，天津 300072)

0 引言

統(tǒng)計(jì)過程控制(Statistical Process Control, SPC)是質(zhì)量管理領(lǐng)域的常用方法之一。此方法將統(tǒng)計(jì)技術(shù)應(yīng)用于過程監(jiān)控，通過分析工具觀察過程的運(yùn)行狀態(tài)，對(duì)識(shí)別出的異常進(jìn)行處理，以確保過程處于穩(wěn)定、可接受的狀態(tài)從而保證產(chǎn)品和服務(wù)的質(zhì)量。

產(chǎn)品的質(zhì)量特征由一元逐漸發(fā)展為多元，而且近年來輪廓型數(shù)據(jù)的監(jiān)控問題已成為質(zhì)量領(lǐng)域關(guān)注的熱點(diǎn)問題。產(chǎn)品或過程的質(zhì)量被描述為一個(gè)函數(shù)關(guān)系，其響應(yīng)變量對(duì)應(yīng)一個(gè)或多個(gè)解釋變量。這種函數(shù)關(guān)系稱為輪廓(profile)，或函數(shù)數(shù)據(jù)(Functional Data)。Woodall等[1]將輪廓監(jiān)控分為兩個(gè)階段：第一階段是以采集到的歷史樣本數(shù)據(jù)為基礎(chǔ)，判斷失控樣本點(diǎn)并移除，建立穩(wěn)定的狀態(tài)模型；第二階段是在第一階段建立受控模型的基礎(chǔ)上對(duì)實(shí)時(shí)數(shù)據(jù)進(jìn)行監(jiān)控，以及時(shí)發(fā)現(xiàn)過程中出現(xiàn)的變異。

本文研究的異常點(diǎn)識(shí)別問題屬于輪廓監(jiān)控第一階段。由于第二階段是根據(jù)第一階段建立的模型進(jìn)行實(shí)時(shí)判斷，因此在第一階段的模型中有效識(shí)別隱藏在所有觀測(cè)點(diǎn)中的小部分異常數(shù)據(jù)點(diǎn)是非常重要的。異常點(diǎn)(outliers)通常預(yù)示著系統(tǒng)中的非正常狀態(tài)，這些異常點(diǎn)可以提示一些有用的信息，甚至是能夠有一些重大發(fā)現(xiàn)。關(guān)于異常點(diǎn)識(shí)別問題至今在各領(lǐng)域已有較多應(yīng)用，如網(wǎng)絡(luò)安全[2～5]，視覺監(jiān)控[6,7]，遙感技術(shù)[8,9]，內(nèi)科診斷學(xué)[10]等。

很多學(xué)者致力于線性輪廓異常點(diǎn)監(jiān)控問題的研究[11,12]，而目前在醫(yī)療、工程、信息技術(shù)等領(lǐng)域的理論研究與工程應(yīng)用方面，輪廓數(shù)據(jù)更多表現(xiàn)為非線性，即輪廓線的自變量與因變量之間的函數(shù)關(guān)系是非線性的。因此需要建立更加復(fù)雜的模型加以描述。不同學(xué)者在解決非線性輪廓異常點(diǎn)識(shí)別問題時(shí)提出了不同的方法[13,14]，很多數(shù)據(jù)處理技術(shù)被應(yīng)用其中。當(dāng)前研究多以隨機(jī)變異服從正態(tài)分布為基本假設(shè)，然而在實(shí)際應(yīng)用中大量存在非正態(tài)分布的情況，如經(jīng)濟(jì)指標(biāo)多為偏態(tài)分布、產(chǎn)品壽命通常服從威布爾分布等。當(dāng)正態(tài)分布假設(shè)不滿足時(shí)，異常點(diǎn)識(shí)別性能會(huì)受到影響。

本文提出將小波分析、數(shù)據(jù)深度、聚類分析等技術(shù)相結(jié)合的非線性輪廓異常點(diǎn)識(shí)別方法，簡稱為WDC法。小波分析(Wavelet Analysis)是數(shù)據(jù)降噪的處理方法。數(shù)據(jù)深度(Data Depth)是一種通過對(duì)高維數(shù)據(jù)進(jìn)行排序來描述數(shù)據(jù)的離心程度的方法。深度值越大，表示數(shù)據(jù)點(diǎn)越接近數(shù)據(jù)中心。聚類分析是依據(jù)數(shù)據(jù)自身的特征值將集合內(nèi)的所有數(shù)據(jù)聚集為若干具有相似特征的數(shù)據(jù)類的方法。在以數(shù)據(jù)深度識(shí)別異常點(diǎn)的方法中均采用設(shè)定閾值的方法進(jìn)行判定[15,16]，然而本文運(yùn)用K-means聚類區(qū)分正常輪廓和異常點(diǎn)，經(jīng)驗(yàn)證表明準(zhǔn)確性更高。選擇以上技術(shù)的另一個(gè)原因是其對(duì)正態(tài)分布假設(shè)的依賴性不強(qiáng)，應(yīng)用范圍廣，經(jīng)仿真驗(yàn)證在非正態(tài)分布變異下具有良好的性能。本文使用Matlab仿真的方法將新方法與Zhang等[20]提出的χ2控制圖方法進(jìn)行對(duì)比，結(jié)果顯示新方法的性能更優(yōu)。最后，采用木板垂直密度輪廓經(jīng)典案例證明新方法在實(shí)際應(yīng)用中是有效的。

1 數(shù)據(jù)模型及流程研究

假設(shè)研究對(duì)象是N條非線性輪廓，每一條輪廓線由M個(gè)點(diǎn)組成，每個(gè)點(diǎn)(xij,yij)是一個(gè)觀測(cè)值，因此可以把這N條非線性輪廓看成N×M的矩陣。N條輪廓線所構(gòu)成的樣本均為隨機(jī)抽樣得到。輪廓線數(shù)據(jù)中包含正常輪廓線和異常輪廓線。每條輪廓線的一般函數(shù)表示為：

yij=f(xij)+εij

(1)

其中i=1,2,…,N；j=1,…,M；εij為噪聲因子，εij∈N(0,1)。

當(dāng)M的取值較大時(shí)，輪廓表現(xiàn)為高維特性。同時(shí)，非線性函數(shù)f具有相對(duì)復(fù)雜的特點(diǎn)。基于以上數(shù)據(jù)模型及非線性輪廓特點(diǎn)，本文提出的WDC法基本流程步驟如下所示：

第一步: 隨機(jī)生成500條包含異常點(diǎn)的輪廓線，通過小波降噪的方法去除輪廓線中的噪聲因子。

第二步：將降噪后的輪廓線通過馬氏深度的方法轉(zhuǎn)化為數(shù)據(jù)深度值。

第三步：對(duì)輪廓線對(duì)應(yīng)的數(shù)據(jù)深度值進(jìn)行K-Means聚類分析，取聚類數(shù)量為2，將數(shù)量較少的類做為異常點(diǎn)的集合。

當(dāng)每條輪廓所包含的點(diǎn)的數(shù)量較大時(shí)，噪聲因子將在很大程度上影響輪廓形狀的變異性。因此，本文采用小波分析對(duì)原始輪廓數(shù)據(jù)進(jìn)行降噪處理，盡可能降低噪聲因子對(duì)異常輪廓線識(shí)別效果的影響。由于輪廓非線性函數(shù)的復(fù)雜性，無法以精確的函數(shù)模型描述其形狀特征。因此，本文將輪廓視為高維空間中的點(diǎn)，并以其離心程度作為判斷異常點(diǎn)的依據(jù)，從而實(shí)現(xiàn)在不需要對(duì)輪廓進(jìn)行函數(shù)建模的前提下識(shí)別異常點(diǎn)。經(jīng)典的多元統(tǒng)計(jì)方法在高維度的情況下效果受到一定的限制。當(dāng)前普遍采用降維的方法進(jìn)行預(yù)理，如變量選擇法和投影法。然而，降維方法會(huì)造成一定程度上的信息損失，不能全面反映每個(gè)數(shù)據(jù)的變異。對(duì)于輪廓型數(shù)據(jù)，某些部分?jǐn)?shù)據(jù)的變化表現(xiàn)為局部形狀的變化。如果在降維過程中忽視了這些部分的變化，則不能準(zhǔn)確的識(shí)別出表現(xiàn)為局部變化的異常輪廓。然而對(duì)于高維的非線性輪廓，局部變化是不能被忽視的，它是異常輪廓類型的重要組成部分。因此，本文利用數(shù)據(jù)深度中的馬氏深度對(duì)輪廓的離心程度進(jìn)行度量，能夠有效避免信息損失。為了識(shí)別異常點(diǎn)，對(duì)轉(zhuǎn)化后的數(shù)據(jù)深度值進(jìn)行K-Means聚類分析。從高維空間上將異常點(diǎn)與占有絕對(duì)多數(shù)的正常點(diǎn)區(qū)分開來。最終實(shí)現(xiàn)對(duì)非線輪廓異常點(diǎn)的準(zhǔn)確識(shí)別。

2 方法技術(shù)

2.1 小波降噪

小波降噪是小波變換在信號(hào)處理領(lǐng)域的重要應(yīng)用之一，本質(zhì)是抑制信號(hào)中的無用部分而增強(qiáng)信號(hào)中有用部分的過程。通常情況下實(shí)際測(cè)試中含有大量的噪聲，這些噪聲會(huì)影響到分析結(jié)果和分析精度，進(jìn)而對(duì)得出結(jié)論的準(zhǔn)確性產(chǎn)生影響，因此小波降噪處理是非常有必要的。

含有噪聲的一維信號(hào)模型：

y(t)=f(t)+σe(t)

(2)

其中t=0,…,n-1，f(t)為真實(shí)信號(hào)，σ為噪聲強(qiáng)度，e(t)為噪聲，y(t)為含噪聲的信號(hào)。信號(hào)的噪聲一般情況下可以分為有色噪聲和白色噪聲兩類。對(duì)于有色噪聲的消除可以先采用對(duì)有色噪聲白色化再按照白色噪聲的方式進(jìn)行處理。因此，小波降噪的重點(diǎn)在于消除信號(hào)中的白噪聲[17]。

假設(shè)e(t)為高斯白噪聲，則噪聲強(qiáng)度為1。通常情況下有用信號(hào)表現(xiàn)為一些較平穩(wěn)的信號(hào)或低頻信號(hào)，而噪聲信號(hào)表現(xiàn)為高頻信號(hào)。

對(duì)混合信號(hào)y(t)消噪的目的在于消除噪聲部分σe(t)，從y(t)中恢復(fù)出真實(shí)信號(hào)f(t)。小波降噪過程可分為三個(gè)步驟進(jìn)行處理[18]：

(1)對(duì)信號(hào)進(jìn)行小波分解，通過分解計(jì)算，分解為低頻段和高頻段；

(2)噪聲因子通常在高頻段信號(hào)中，因此對(duì)高頻段信號(hào)以閾值量化等形式對(duì)小波系數(shù)進(jìn)行處理；

(3)對(duì)信號(hào)進(jìn)行小波重構(gòu)以實(shí)現(xiàn)降噪。根據(jù)小波分解的低頻段系數(shù)和高頻段系數(shù)進(jìn)行小波重構(gòu)。

小波降噪步驟圖如圖1所示：

圖1 小波降噪步驟圖

2.2 馬氏深度

數(shù)據(jù)深度通過對(duì)高維數(shù)據(jù)排序來描述數(shù)據(jù)的離心程度，以數(shù)據(jù)yi的深度D(yi)來描述yi在數(shù)據(jù)群中位置。數(shù)據(jù)深度越大，表示數(shù)據(jù)點(diǎn)越靠近數(shù)據(jù)群中心，數(shù)據(jù)點(diǎn)的信息量越大、重要性越高。目前，數(shù)據(jù)深度包含馬氏深度(Mahalanobis depth)、半空間深度(Half-space depth)、空間秩深度(Spatial rank depth)、單純性深度(The simplicial depth)。本文采用馬氏深度法經(jīng)測(cè)試有較好的性能。

馬氏深度的提出源于馬氏距離。設(shè)x,y∈Rd，則x和y關(guān)于正定矩陣Md×d的馬氏距離為[19]：

(3)

定義x關(guān)于分布F的馬氏深度(Mahalanobis)為：

(4)

亦可寫作：

(5)

其中μF和∑F分別是對(duì)應(yīng)于分布F的均值向量和協(xié)方差矩陣。

2.3 K-Means聚類

K-means算法作為經(jīng)典的聚類分析算法中的一種，最早由MacQueen在1967年提出，目前是在眾多領(lǐng)域應(yīng)用較多的有效算法之一。K-means算法是根據(jù)樣本的相似度將樣本聚集到K個(gè)簇當(dāng)中，使簇內(nèi)樣本相似度高而簇間相似度低。

給定一個(gè)包含n個(gè)d維數(shù)據(jù)點(diǎn)的數(shù)據(jù)集X={x1,x2,…,xi,…,xn}，其中xi∈Rd。K-means算法將該數(shù)據(jù)集劃分為K類，C={ck,i=1,…,K}，每個(gè)類ck的質(zhì)心為mk，選取歐式距離作為相似性和距離判斷準(zhǔn)則，計(jì)算各類總的距離平方和E(C)，聚類目標(biāo)是使E(C最小。

(6)

圖2 K-means算法步驟流程圖

3 性能分析

本文使用MATLAB對(duì)非線性輪廓異常點(diǎn)識(shí)別問題進(jìn)行仿真并檢測(cè)本文所提方法識(shí)別異常點(diǎn)的效果，并與Zhang等[20]所提出的用χ2方法進(jìn)行對(duì)比。在本節(jié)依次介紹仿真過程、效果評(píng)價(jià)指標(biāo)、仿真性能結(jié)果等。

3.1 仿真過程

在仿真過程中，首先生成隨機(jī)的N條輪廓線，即N×M矩陣，然后進(jìn)行小波降噪，減少噪聲因子，進(jìn)而運(yùn)用馬氏深度的方法以及K-means聚類將異常點(diǎn)識(shí)別出來。為保證方法的有效性，仿真生成500條輪廓線(包含正常輪廓和異常輪廓)用所提方法進(jìn)行識(shí)別異常點(diǎn)。本仿真驗(yàn)證中以對(duì)數(shù)正態(tài)分布為輪廓的隨機(jī)變異分布特征，目的是驗(yàn)證新方法對(duì)非正態(tài)變異的識(shí)別效果。以式(7)表示生成的正常輪廓曲線：

Y=3cos(x)+5sin(x)+εN

(7)

其中x∈(0, 2π)，εN為噪聲因子，服從對(duì)數(shù)正態(tài)分布。εN的數(shù)學(xué)期望和方差分別為μN(yùn)=0,σN=0.7，即ε服從期望為0，標(biāo)準(zhǔn)差為0.7的對(duì)數(shù)正態(tài)分布。

以式(8)表示生成的異常輪廓曲線：

Y=3cos(x)+5sin(x)+εA

(8)

其中x∈(0,2π)，εA為噪聲因子，服從對(duì)數(shù)正態(tài)分布，σA=0.7。與正常輪廓線不同之處在于εA的期望μA是一個(gè)變化值，μA=0.1∶0.1∶2，即μA從0.1開始以0.1為步長，均勻取20個(gè)點(diǎn)至2止。圖3展示出5條正常輪廓曲線和5條μA為0.5的異常輪廓曲線。

圖3 輪廓線示意圖

3.2 效果評(píng)價(jià)指標(biāo)

在本研究中采用第一錯(cuò)誤率(Type Error I)和第二錯(cuò)誤率(Type Error II)作為異常點(diǎn)識(shí)別的效果評(píng)價(jià)指標(biāo)。

第一錯(cuò)誤率為生產(chǎn)過程中穩(wěn)定過程被判定為不穩(wěn)定的概率，即正常輪廓判定為異常點(diǎn)的數(shù)量與所有正常輪廓線數(shù)量的比值。

(9)

其中，M1為正常輪廓錯(cuò)判為異常點(diǎn)的數(shù)量；PN為正常輪廓數(shù)量。

相反，第二錯(cuò)誤率為不穩(wěn)定誤判為穩(wěn)定的概率，即異常點(diǎn)誤判為正常輪廓的數(shù)量與所有異常點(diǎn)數(shù)量的比值。具體計(jì)算如下：

(10)

其中，M2為異常點(diǎn)誤判為正常輪廓的數(shù)量；PA為異常點(diǎn)數(shù)量。

前文提到異常輪廓線的誤差項(xiàng)期望μA是一個(gè)變化值，下節(jié)仿真結(jié)果將討論隨μA值逐漸增大，即異常輪廓線與正常輪廓線差異逐漸增大的時(shí)候，TypeErrorI、TypeErrorII如何變化。

3.3 仿真結(jié)果對(duì)比

為了探究在異常點(diǎn)數(shù)量不同的情況下，本文所提方法的性能效果的差異，本文引入異常點(diǎn)占比(記作Perc)指標(biāo)，研究在異常點(diǎn)占總輪廓線比例(Perc)不同的情況下，異常輪廓線的誤差項(xiàng)期望μA與TypeErrorI、TypeErrorII的函數(shù)關(guān)系。

仿真共隨機(jī)生成500條輪廓線，考慮異常點(diǎn)占總輪廓線比例(Perc)分別為0.10、0.15、0.20這三種情況下，即異常輪廓線數(shù)量分別為50、75、100。按照本文所提方法在仿真中重復(fù)循環(huán)1000次，將本文所提方法性能與Zhang等[20]方法性能進(jìn)行對(duì)比。圖4展示的異常點(diǎn)占比(Perc)分別為0.10、0.15、0.20時(shí)μA與TypeErrorI對(duì)比圖，實(shí)線表示本文所提方法，虛線表示對(duì)比方法。

圖4 第一錯(cuò)誤率性能效果對(duì)比圖

表1為第一類錯(cuò)誤的對(duì)比數(shù)據(jù)，由數(shù)據(jù)和圖像可以得到以下幾點(diǎn)結(jié)論：

(1)本文所提方法的第一類錯(cuò)誤均明顯小于對(duì)比方法；

(2)在μA相同情況下，對(duì)比方法的第一類錯(cuò)誤和異常點(diǎn)占比(Perc)成明顯反比，即TypeErrorI隨Perc升高而減小。而本文所提方法在μA相同情況下，TypeErrorI隨Perc升高而略有減小甚至變化不大；

(3)隨μA增長，本文所提方法性能的差異性小于對(duì)比方法，說明本方法的性能在不同情況下具有一定的穩(wěn)定性。

表1 第一錯(cuò)誤率(Type Error I)性能效果數(shù)據(jù)表

圖5 第二錯(cuò)誤率性能效果對(duì)比圖

圖5展示的異常點(diǎn)占比(Perc)分別為0.10、0.15、0.20時(shí)μA與TypeErrorII對(duì)比圖，實(shí)線表示本文所提方法，虛線表示對(duì)比方法。

表2為第二類錯(cuò)誤的對(duì)比數(shù)據(jù)，由數(shù)據(jù)和圖像可以得到以下幾點(diǎn)結(jié)論：

(1)本文所提方法的第二類錯(cuò)誤均小于對(duì)比方法；(2)在μA相同情況下，本文所提方法和對(duì)比方法的第二類錯(cuò)誤和異常點(diǎn)占比(Perc)均成正比，即TypeErrorII隨Perc升高而增大。

(3)異常點(diǎn)占比(Perc)不同時(shí)，本文所提方法的μA-TypeErrorII函數(shù)曲線變化幅度小于對(duì)比方法。即Perc不同時(shí)，本文所提方法性能穩(wěn)定性優(yōu)于對(duì)比方法。

表2 第二錯(cuò)誤率(Type Error II)性能效果數(shù)據(jù)表

綜上所述，經(jīng)過仿真分析得到的性能對(duì)比結(jié)果可知，對(duì)于非正態(tài)變異的非線性輪廓線本文所提方法的第一錯(cuò)誤率明顯小于對(duì)比方法，第二錯(cuò)誤率略小于對(duì)比方法。由此可知，本文所提方法在解決異常點(diǎn)識(shí)別問題時(shí)效果優(yōu)于Zhang等[20]所提出的用χ2方法。

4 實(shí)例分析

在本節(jié)中，引入垂直密度輪廓(Vertical Density Profile, VDP)案例對(duì)本文所提方法的實(shí)際應(yīng)用效果進(jìn)行驗(yàn)證。本案例最先由Walker和Wright[21]提出，進(jìn)而很多質(zhì)量監(jiān)控領(lǐng)域?qū)W者在研究過程中均對(duì)其進(jìn)行分析，主要研究木板厚度方向上的密度變化。通常情況下，纖維木板為中間芯層斷面密度高，兩側(cè)斷面密度低。通常使用激光準(zhǔn)確掃描模板各個(gè)位置的垂直密度來得到不同輪廓曲線。每張木板的垂直密度輪廓線體現(xiàn)了其產(chǎn)品質(zhì)量水平。VDP數(shù)據(jù)共包含24條輪廓，每條輪廓線包含314個(gè)觀測(cè)值。該24條輪廓線呈現(xiàn)出浴盆形態(tài)，如圖6所示。

圖6 VDP輪廓數(shù)據(jù)圖

采用本文所提出的WDC方法對(duì)VDP數(shù)據(jù)進(jìn)行異常點(diǎn)識(shí)別，利用Matlab軟件進(jìn)行分析計(jì)算，得到的結(jié)果為：輪廓線3、6、16、11、21為異常點(diǎn)。圖7展示了VDP數(shù)據(jù)的異常點(diǎn)識(shí)別結(jié)果。由圖像可看出，第3條、第6條輪廓線明顯偏移整體輪廓線，第16條輪廓線整體位于偏下方的位置，偏離于主體輪廓線，故這三條輪廓線符合異常輪廓的特征。對(duì)于輪廓線11、21，可以看出這兩條輪廓線的兩端偏離于主體輪廓線的兩端位置，輪廓線21兩端高于主體輪廓線的兩端而輪廓線11兩端低于主體輪廓線兩端，但是由于和主體輪廓形態(tài)沒有明顯差異，因此判定輪廓線11、21為疑似異常點(diǎn)。

圖7 WDC方法識(shí)別出的異常點(diǎn)

5 結(jié)論及展望

針對(duì)非線性輪廓異常點(diǎn)識(shí)別問題，本文提出了基于小波降噪、馬氏深度和聚類分析的WDC方法，此方法首次使用聚類分析識(shí)別轉(zhuǎn)化為數(shù)據(jù)深度值的輪廓線中存在的異常點(diǎn)。通過MATLAB仿真測(cè)試新方法的性能并與χ2控制圖方法對(duì)比。結(jié)果顯示對(duì)于非正態(tài)變異的非線性輪廓新方法的第一錯(cuò)誤率明顯優(yōu)于對(duì)比方法并且第一錯(cuò)誤率隨異常點(diǎn)占比改變變化不明顯，同時(shí)第二錯(cuò)誤率也優(yōu)于對(duì)比方法。WDC方法在性能與穩(wěn)定性方面均優(yōu)于對(duì)比方法。

由于小波降噪、數(shù)據(jù)深度等技術(shù)仍有很多不同的方法，因此如何合理選擇相應(yīng)的技術(shù)，針對(duì)不同對(duì)象獲得最佳的分析效果仍有待深入研究。另外，非正態(tài)分布的分布類型也有很多種，在對(duì)數(shù)正態(tài)分布以外的其他類型過程中本方法的性能表現(xiàn)仍有待研究。非線性輪廓的形態(tài)也會(huì)影響識(shí)別效果，因此仍有必要對(duì)新方法的適用性進(jìn)行全面探討。