基于PCA與SPRT的機(jī)械故障診斷方法研究

黃文健，黃瑾珉，曹承昊，陳漢新

武漢工程大學(xué)機(jī)電工程學(xué)院，湖北 武漢 430205

齒輪箱是很多機(jī)械設(shè)備中重要的部件，其常被用來傳遞動力、改變轉(zhuǎn)速、分配動力等。在實(shí)際應(yīng)用中齒輪箱的工作狀態(tài)與設(shè)備的運(yùn)行狀態(tài)息息相關(guān)，而其中的齒輪，軸，軸承等零部件長期在高負(fù)荷狀態(tài)下運(yùn)行，故障率較高，很容易引發(fā)機(jī)器故障。在齒輪箱部件失效中由齒輪故障引發(fā)的約占60%［1］，因此對齒輪箱中齒輪故障診斷的研究就顯得很有必要了。

齒輪常見的故障有齒輪磨損、齒面膠合和擦傷、齒面接觸疲勞和斷齒等。在本次實(shí)驗(yàn)中，選取齒輪裂紋來模擬齒輪故障，得到的實(shí)驗(yàn)信號摻雜著許多噪聲，所以采用小波包變換對信號進(jìn)行降噪［2］，之后再提取特征。主元分析法（Principal Component Analysis，PCA）［3］是將多個線性相關(guān)變量壓縮為少數(shù)幾個不相關(guān)的變量的一種多元統(tǒng)計(jì)方法。葉永偉等［4］用PCA來對設(shè)備特征進(jìn)行提?。欢耪駥幍龋?］將小波包分解與PCA相結(jié)合來對軸承進(jìn)行故障診斷。由于PCA具有概念簡單，計(jì)算方便等優(yōu)點(diǎn)，因此常被用來對數(shù)據(jù)進(jìn)行降維。

Wald［6］在1947年的著作中系統(tǒng)地介紹了序貫概 率 比 檢 驗(yàn)（Sequential Probability Ratio Test，SPRT）。由于其具有簡單、高效等特性，因此近些年來在故障診斷方面使用非常普遍。Goodman等［7］將SPRT與自適應(yīng)波形結(jié)合起來識別雷達(dá)信號；Ray等［8］將SPRT應(yīng)用在核電站和飛行器的故障檢測與識別；CHEN等［9-10］將SPRT用于齒輪裂紋的故障診斷。SPRT與固定樣本檢驗(yàn)不同之處在于這種算法需要更少的平均樣本量，但是卻具有更高的檢驗(yàn)效率。

在齒輪的故障診斷中，為了避免直接從眾多特征參數(shù)中選取SPRT的參數(shù)，本文將PCA與SPRT相結(jié)合起來應(yīng)用于齒輪裂紋的故障診斷。實(shí)驗(yàn)中采集的信號在降噪之后，利用PCA來處理信號的特征數(shù)據(jù)并選擇合適的主元作為檢驗(yàn)參數(shù)以提高檢驗(yàn)的準(zhǔn)確度與可靠性，通過分析最后得到的結(jié)果圖來進(jìn)行識別與分類［11-12］，結(jié)果表明該方法是有效的。

1 PCA

特征提取是從信號里提取有效信息的過程，適當(dāng)?shù)姆椒ê陀行У奶卣鲗罄m(xù)過程很有幫助。在很多研究中對數(shù)據(jù)進(jìn)行降維會使之后的計(jì)算分析更加的高效與簡便。在用PCA研究多變量問題時，實(shí)際上是通過線性變化將數(shù)據(jù)從高維空間映射到低維空間，低維空間中較少的新變量能集中地反映原始變量中包含的絕大部分信息。其基本方法是將數(shù)據(jù)矩陣的協(xié)方差矩陣的特征值按大小進(jìn)行排序，以此順序計(jì)算各個主元及其貢獻(xiàn)率，之后選取滿足要求的主元從而達(dá)到降維的目的。目前PCA已廣泛用于特征提取、數(shù)據(jù)降噪以及故障重構(gòu)等方面。

PCA的基本原理如下：設(shè)原始數(shù)據(jù)樣本包含m個n維特征向量 Xk=（ x1，x2，… xn）T，則PCA具體步驟如下：

計(jì)算樣本的平均值：

求出樣本向量的協(xié)方差矩陣：

計(jì)算S的特征值 λ1，λ2，…，λn及其對應(yīng)的特征向量 v1，v2，…，vn。

計(jì)算各主元的貢獻(xiàn)率和累計(jì)貢獻(xiàn)率：

選取前p個累計(jì)貢獻(xiàn)率滿足要求的主元作為映射矩陣，將原始信號變換到特征空間：

式（5）中，V=（v1，v2，…vp）。

2 SPRT理論

SPRT的基本原理為：設(shè) x1，x2，…，xn是一組獨(dú)立同分布的隨機(jī)變量序列。對于這樣的總體樣本值，提出兩個假設(shè)，零假設(shè)（H0∶θ=θ0）和備選假設(shè)（H1∶θ=θ1）。它們的聯(lián)合概率密度可以定義為：

其中j=0，1，檢驗(yàn)的似然比計(jì)算如下：

式（1）中，f(x/θ)為條件概率分布；θ為分布參數(shù)。

在以前的檢驗(yàn)中，當(dāng)λn的值在閾值C附近時，對兩種假設(shè)的接受與否的判斷顯得不夠合理。而Wald改進(jìn)的地方是當(dāng) λn遠(yuǎn)小于C時接受 H0，當(dāng)λn遠(yuǎn)大于C時拒絕 H0，當(dāng) λn與C很接近時暫且不做出判斷，而是再繼續(xù)檢測λn+1，直到之后的似然比遠(yuǎn)大于或遠(yuǎn)小于C為止。

在SPRT中，事先給定兩類錯誤的概率α與 β的值，則閾值A(chǔ)、B可以按照下列公式來確定：

檢驗(yàn)判斷準(zhǔn)則為：代入觀測值x1，計(jì)算λ1(x)，若 λ1(x)＜B，則停止觀測并接受零假設(shè) H0；若λ1(x)＞A，則停止觀測并接受備選假設(shè) H1；若A≤λ1(x)≤B，則繼續(xù)進(jìn)行檢驗(yàn)，直至似然比值滿足上述兩種情況為止。

3 實(shí)驗(yàn)系統(tǒng)

選取正常齒輪F1和有裂紋的故障齒輪F2、F3進(jìn)行實(shí)驗(yàn)。齒輪裂紋全深度a=2.4 mm，裂紋全寬度b=25 mm，厚度為0.4 mm，齒輪裂紋的角度為45°。實(shí)驗(yàn)中齒輪箱為空載狀態(tài)，轉(zhuǎn)速為800 r/min。3組齒輪的參數(shù)如表1所示。

表1 三種故障模式Tab.1 Three fault modes

圖1是齒輪箱的結(jié)構(gòu)圖。在運(yùn)行中齒輪3和4之間會產(chǎn)生振動，因此可以任選其中之一來模擬故障。本次實(shí)驗(yàn)選擇用齒輪3來進(jìn)行模擬實(shí)驗(yàn)。

圖1 齒輪箱結(jié)構(gòu)圖Fig.1 Structure of gearbox

齒輪箱信號采集系統(tǒng)如圖2所示，將加速度傳感器安裝在齒輪箱的水平垂直兩個方向上，通過動態(tài)模擬器采集振動信號，并通過信號分析器將信號輸入并存儲在計(jì)算機(jī)中，本文僅研究分析水平方向的信號。分別用S1、S2和S3表示正常情況與兩種故障情況下的信號。

圖2 齒輪箱信號采集系統(tǒng)Fig.2 Signal acquisition system of gearbox

4 齒輪裂紋的故障診斷

4.1 小波包降噪

小波包變換將同時分解信號的低頻和高頻部分，這有利于有效地區(qū)分信號中的噪聲。對實(shí)驗(yàn)信號進(jìn)行3層小波包分解，得到8個頻帶的分解信號后再重構(gòu)信號。這樣重構(gòu)后的信號會更準(zhǔn)確。

4.2 選取與計(jì)算特征參數(shù)

振動信號的特征信息是通過一些特征參數(shù)來表達(dá)的，每個特征參數(shù)都能反映該信號的某種特性，平穩(wěn)的振動信號表明齒輪箱運(yùn)行正常，而異常的或是有較大波動的振動信號則表明齒輪箱出現(xiàn)了故障。在常用的時域特征參數(shù)中，均值能反映信號波動能量的大??；有效值是描述動態(tài)信號強(qiáng)度的指標(biāo)；標(biāo)準(zhǔn)差是描述信號的波動分量；峭度指標(biāo)能反應(yīng)振動信號中的沖擊特征；波形指標(biāo)是反映信號偏離高斯分布程度的指標(biāo)；脈沖指標(biāo)與峰值指標(biāo)都是用于檢測信號中有無沖擊的指標(biāo)；裕度指標(biāo)常被用來檢測機(jī)械設(shè)備的磨損狀況［13-15］。本文選取這8個參數(shù)指標(biāo)進(jìn)行計(jì)算。假設(shè)xi=[x1，x2，…，xN]是一組離散待檢信號，N=8 192。每組取1 024個檢驗(yàn)點(diǎn)，這樣可以得到7196組檢驗(yàn)數(shù)據(jù)。

對實(shí)驗(yàn)信號進(jìn)行降噪之后再提取特征值，本文中n=1 024，特征參數(shù)的計(jì)算公式如下：

將經(jīng)過預(yù)處理的振動信號代入公式中分別計(jì)算，得到齒輪箱在3種狀況下的3組信號的8種特征參數(shù)。

4.3 特征參數(shù)的降維

對上述3組數(shù)據(jù)進(jìn)行PCA處理，根據(jù)得到的結(jié)果來選取待檢驗(yàn)的序列。

經(jīng)過PCA處理之后，根據(jù)貢獻(xiàn)率的高低來對處理后的主元進(jìn)行降序排列，排序越靠前的主元保留的原有信息越多。在得到的處理結(jié)果中，通常是選取累計(jì)貢獻(xiàn)率高于85%的前幾個主元，這樣選取的主元才能包含足夠多的原有信息。

這些特征值經(jīng)過PCA處理之后，前4個主元就能得到超過85%的信息。這里選取貢獻(xiàn)率最大的序列作為檢驗(yàn)序列，記作 yi=[y1…yk]。其均值和標(biāo)準(zhǔn)差為：

4.4 SPRT

經(jīng)過以上處理之后的待檢序列大致符合正態(tài)分布。正常狀態(tài)下，該序列滿足零假設(shè)H0∶θ=θ0；故障狀態(tài)下，該序列滿足備擇假設(shè)H1∶θ=θ1。標(biāo)準(zhǔn)差保持為σ。當(dāng)兩個假設(shè)都成立時，該序列的聯(lián)合概密度為：

式（20）和（21）中：P0i為零假設(shè)條件下的概率密度函數(shù)；P1i為備擇假設(shè)條件下的概率密度函數(shù)。

SPRT似然比為：

式（22）中：P0為零假設(shè)條件下的先驗(yàn)概率；P1為備擇假設(shè)條件下的先驗(yàn)概率。

對似然比公式進(jìn)一步推導(dǎo)簡化可使得計(jì)算更加簡便：

簡化之后的閾值a=lnA，b=lnB。假設(shè)檢驗(yàn)范兩類錯誤的概率分別取α=0.001，β=0.001，則計(jì)算可得a=6.91，b=-6.91。計(jì)算出似然比值后根據(jù)檢驗(yàn)準(zhǔn)則進(jìn)行判斷。如果Δ＜b，則停止檢驗(yàn)并接受H0，即為正常狀態(tài)；如果Δ＞a，停止檢驗(yàn)并接受H1，即為故障狀態(tài)；如果b＜Δ＜a，則繼續(xù)比較下一個似然比值與a、b的關(guān)系，如此循環(huán)，直至得出結(jié)果。運(yùn)用此方法對齒輪裂紋進(jìn)行診斷的完整流程如圖3所示。

圖3 齒輪裂紋診斷流程圖Fig.3 Flowchart of gear crack fault diagnosis

4.5 誤差分析

利用均方根誤差（Root Mean Square Error，RMSE）來檢驗(yàn)SPRT應(yīng)用于此的可靠性。均方根誤差計(jì)算公式如下：

5 結(jié)果與分析

5.1 小波包降噪結(jié)果

對3組不同的齒輪分別進(jìn)行實(shí)驗(yàn)，得到圖4（a）中的原始振動信號。其中S1、S2和S3分別為正常齒輪、25%和50%裂紋齒輪的振動信號。3組原始振動信號經(jīng)過降噪之后，得到圖4（b）中的信號。

對比圖 4（a）和圖 4（b）可以看到，與原始振動信號相比，降噪后信號的幅值明顯有所減小，而且信號曲線也更為清晰。圖4（b）基本達(dá)到了預(yù)期效果，這也表明了小波包變換能很好地對原始振動信號進(jìn)行降噪。

5.2 PCA處理結(jié)果

對3組實(shí)驗(yàn)信號進(jìn)行特征提取之后再用主元分析法處理，發(fā)現(xiàn)前4個主元包含超過85%的有效信息，信號S1、S2和S3的各主元貢獻(xiàn)率排序如表2所示。由PCA的處理結(jié)果可知這3組信號的第一主元的貢獻(xiàn)率分別達(dá)到了50.9%、50.2%、51.8%。由于排名第一的主元所占的比例最大，因此被選用來進(jìn)行SPRT。

圖4 實(shí)驗(yàn)信號：（a）原始振動信號，（b）去噪后的信號Fig.4 Experimental signals：（a）original vibration signals，（b）de-noised signals

表2 主元貢獻(xiàn)率及累計(jì)貢獻(xiàn)率Tab.2 Principal component contribution rate and accumulative contribution rate %

5.3 SPRT結(jié)果

由似然比計(jì)算公式可知，平均值的變化對似然比影響較大。將正常信號S1的平均值記為μ0，故障信號S2和S3的平均值記為μ1。對信號進(jìn)行檢測，結(jié)果如圖5所示。

當(dāng)輸入的信號是 S1和S2時，檢驗(yàn)結(jié)果如圖5（a）所示，從圖中可以判定，信號S1的似然比Δ1，2（S1）＜b，即齒輪箱正常工作；信號 S2的似然比Δ1，2（S2）＞a，即齒輪箱出現(xiàn)故障。當(dāng)輸入的信號是S1和S3時，檢驗(yàn)結(jié)果如圖5（a）所示，從圖中可以判定，信號S1的似然比Δ1，3（S1）＜b，即齒輪箱正常工作；信號S3的似然比Δ1，3（S3）＞a，即齒輪箱出現(xiàn)故障。

圖5 SPRT結(jié)果：（a）似然比（Δ1，2（S1），Δ1，2（S2），Δ1，3（S1），Δ1，3（S3））與序貫檢驗(yàn)迭代次數(shù)的關(guān)系，（b）似然比（Δ2，3（S2），Δ2，3（S3））與序貫檢驗(yàn)迭代次數(shù)的關(guān)系Fig.5 Results of SPRT：（a）relationship between SPRT likelihood ratios（Δ1，2（S1），Δ1，2（S2），Δ1，3（S1），Δ1，3（S3））and the number of sequential testing iteration，（b）relationship between SPRT likelihood ratios（Δ2，3（S2），Δ2，3（S3））and the number of sequential testing iteration

在實(shí)驗(yàn)中發(fā)現(xiàn)不同程度裂紋的齒輪也可以通過這種方法加以區(qū)分，其方法基本相同。就本次實(shí)驗(yàn)而言，將兩種故障信號S2和S3進(jìn)行檢測，檢驗(yàn)結(jié)果如圖5（b）所示，從圖中可以判定，信號S2的似然比Δ2，3（S2）＜b，即齒輪具有輕度損傷；信號S3的似然比Δ2，3（S3）＞a，即齒輪出現(xiàn)重度故障。所以根據(jù)檢驗(yàn)結(jié)果可以對25%裂紋和50%裂紋的信號進(jìn)行清晰地區(qū)分。

5.4 誤差結(jié)果分析

在實(shí)驗(yàn)中，正常狀況下獲取兩組數(shù)據(jù)，分別進(jìn)行降噪與PCA處理后共得到14 338個數(shù)據(jù)點(diǎn)。將這些數(shù)據(jù)分為16組，分別記為A1至A16。對16組數(shù)據(jù)分別進(jìn)行序貫概率比檢驗(yàn)，得到16個似然比值，記為B1至B16。計(jì)算這16個值之間的均方根誤差，即B1分別其他15個值比較，B2分別與其他14個值比較，…，B15與B16比較，得到15組比值，選取每組比值中的最大值形成誤差曲線。將故障狀況下的似然比記為a，計(jì)算故障狀況與正常狀況之間的誤差，即a與B1，B2…B16比較，將得到16個誤差值形成誤差曲線。通過比較誤差來驗(yàn)證序貫概率比檢驗(yàn)診斷齒輪裂紋故障的可靠性。

將正常信號S1的平均值記為 μ0，故障信號S2的平均值記為 μ1。將16組正常數(shù)據(jù)信號的似然比Δ1，2（S1）輸入程序中，得到了 16 條 normal曲線；將2組故障信號的似然比Δ1，2（S2），Δ1，2（S3）輸入程序中，得到了2條fault曲線，最后的結(jié)果如圖6（a）所示。圖6（a）表明16組正常曲線的趨勢基本相同，而故障信號的曲線與之相比則有較大的區(qū)別。計(jì)算16組正常信號的似然比之間的均方根誤差E1，2（S1），每組誤差中的最大值形成的誤差曲線如圖6（b）所示。分別計(jì)算故障信號S1和S2與16組正常信號似然比之間的誤差E1，2（S2）和E1，2（S3），得到圖6（b）中的另外兩條誤差曲線。在誤差方面，從圖6（b）中可以看出不同故障之間的誤差遠(yuǎn)大于同種故障之間的，這表明序貫概率比檢驗(yàn)可以很好地區(qū)分不同故障狀況。

圖6 SPRT及誤差結(jié)果：（a）似然比（Δ1，（2S1），Δ1，（2S2），Δ1，（2S3））與序貫檢驗(yàn)迭代次數(shù)的關(guān)系，（b）均方根誤差值E1，（2S1），E1，（2S2）和E1，（2S3）Fig.6 Results of SPRT and RMSE：（a）relationship between SPRT likelihood ratios（Δ1，（2S1），Δ1，（2S2），Δ1，（2S3））and thenumber of sequential testing iteration，（b）RMSEs E1，（2S1），E1，（2S2）and E1，（2S3）

將正常信號S1的平均值記為 μ0，故障信號S3的平均值記為 μ1。將16組正常數(shù)據(jù)信號的似然比Δ1，3（S1）輸入程序中，得到了 16 條 normal曲線；將兩組故障信號的似然比Δ1，3（S2），Δ1，3（S3）輸入程序中，得到了2條fault曲線，最后的結(jié)果如圖7（a）所示。圖7（a）說明16組正常曲線的趨勢基本相同，而故障信號的曲線與之相比則有較大的區(qū)別。計(jì)算16組正常信號的似然比之間的均方根誤差E1，3（S1），每組誤差中的最大值形成的誤差曲線如圖7（b）所示。分別計(jì)算故障信號S1和S2與16組正常信號似然比之間的誤差E1，2（S2）和E1，3（S3），得到圖7（b）中的另外2條誤差曲線。在誤差方面，從圖7（b）中可以看出不同故障之間的誤差遠(yuǎn)大于同種故障之間的，這表明SPRT可以很好地區(qū)分不同故障狀況。

圖7 SPRT及均方根誤差結(jié)果：（a）似然比（Δ1，（3S1），Δ1，（3S2），Δ1，（3S3））與序貫檢驗(yàn)迭代次數(shù)的關(guān)系，（b）均方根誤差值E1，（3S1），E1，（3S2）和E1，（3S3）Fig.7 Results of SPRT and RMSE：（a）relationship between SPRT likelihood ratios（Δ1，（3S1），Δ1，（3S2），Δ1，（3S3））and thenumber of sequential testing iteration；（b）RMSEs E1，（3S1），E1，（3S2）and E1，（3S3）

由圖6和圖7可知，當(dāng)以正常信號的平均值為參數(shù)μ0時，不論是以哪種故障信號的平均值為參數(shù) μ1來建立模型，從最后的結(jié)果中都能對3種狀況進(jìn)行清晰地區(qū)分。

6 結(jié) 語

本文提出了一種將PCA與SPRT相結(jié)合的檢測方法，并選取不同程度裂紋的故障齒輪與正常齒輪進(jìn)行模擬試驗(yàn)，來檢測此方法能否有效區(qū)分不同狀態(tài)下的齒輪。在預(yù)處理之后利用PCA對實(shí)驗(yàn)信號所選取的特征參數(shù)進(jìn)行降維處理，有效地減少了特征數(shù)量，使得之后的檢驗(yàn)更加高效準(zhǔn)確，最終有效地識別出了齒輪箱不同的狀態(tài)。此外，經(jīng)過多次實(shí)驗(yàn)后證明在其他的工作條件下，該方法也能對齒輪箱不同的狀態(tài)進(jìn)行清晰地區(qū)分。因此將該方法應(yīng)用于齒輪裂紋故障檢測與分類方面是有效而且可靠的。此方法也可供其他旋轉(zhuǎn)機(jī)械的診斷研究進(jìn)行參考。