黃文健,黃瑾珉,曹承昊,陳漢新
武漢工程大學(xué)機(jī)電工程學(xué)院,湖北 武漢 430205
齒輪箱是很多機(jī)械設(shè)備中重要的部件,其常被用來(lái)傳遞動(dòng)力、改變轉(zhuǎn)速、分配動(dòng)力等。在實(shí)際應(yīng)用中齒輪箱的工作狀態(tài)與設(shè)備的運(yùn)行狀態(tài)息息相關(guān),而其中的齒輪,軸,軸承等零部件長(zhǎng)期在高負(fù)荷狀態(tài)下運(yùn)行,故障率較高,很容易引發(fā)機(jī)器故障。在齒輪箱部件失效中由齒輪故障引發(fā)的約占60%[1],因此對(duì)齒輪箱中齒輪故障診斷的研究就顯得很有必要了。
齒輪常見的故障有齒輪磨損、齒面膠合和擦傷、齒面接觸疲勞和斷齒等。在本次實(shí)驗(yàn)中,選取齒輪裂紋來(lái)模擬齒輪故障,得到的實(shí)驗(yàn)信號(hào)摻雜著許多噪聲,所以采用小波包變換對(duì)信號(hào)進(jìn)行降噪[2],之后再提取特征。主元分析法(Principal Component Analysis,PCA)[3]是將多個(gè)線性相關(guān)變量壓縮為少數(shù)幾個(gè)不相關(guān)的變量的一種多元統(tǒng)計(jì)方法。葉永偉等[4]用PCA來(lái)對(duì)設(shè)備特征進(jìn)行提??;杜振寧等[5]將小波包分解與PCA相結(jié)合來(lái)對(duì)軸承進(jìn)行故障診斷。由于PCA具有概念簡(jiǎn)單,計(jì)算方便等優(yōu)點(diǎn),因此常被用來(lái)對(duì)數(shù)據(jù)進(jìn)行降維。
Wald[6]在1947年的著作中系統(tǒng)地介紹了序貫概 率 比 檢 驗(yàn)(Sequential Probability Ratio Test,SPRT)。由于其具有簡(jiǎn)單、高效等特性,因此近些年來(lái)在故障診斷方面使用非常普遍。Goodman等[7]將SPRT與自適應(yīng)波形結(jié)合起來(lái)識(shí)別雷達(dá)信號(hào);Ray等[8]將SPRT應(yīng)用在核電站和飛行器的故障檢測(cè)與識(shí)別;CHEN等[9-10]將SPRT用于齒輪裂紋的故障診斷。SPRT與固定樣本檢驗(yàn)不同之處在于這種算法需要更少的平均樣本量,但是卻具有更高的檢驗(yàn)效率。
在齒輪的故障診斷中,為了避免直接從眾多特征參數(shù)中選取SPRT的參數(shù),本文將PCA與SPRT相結(jié)合起來(lái)應(yīng)用于齒輪裂紋的故障診斷。實(shí)驗(yàn)中采集的信號(hào)在降噪之后,利用PCA來(lái)處理信號(hào)的特征數(shù)據(jù)并選擇合適的主元作為檢驗(yàn)參數(shù)以提高檢驗(yàn)的準(zhǔn)確度與可靠性,通過(guò)分析最后得到的結(jié)果圖來(lái)進(jìn)行識(shí)別與分類[11-12],結(jié)果表明該方法是有效的。
特征提取是從信號(hào)里提取有效信息的過(guò)程,適當(dāng)?shù)姆椒ê陀行У奶卣鲗?duì)后續(xù)過(guò)程很有幫助。在很多研究中對(duì)數(shù)據(jù)進(jìn)行降維會(huì)使之后的計(jì)算分析更加的高效與簡(jiǎn)便。在用PCA研究多變量問(wèn)題時(shí),實(shí)際上是通過(guò)線性變化將數(shù)據(jù)從高維空間映射到低維空間,低維空間中較少的新變量能集中地反映原始變量中包含的絕大部分信息。其基本方法是將數(shù)據(jù)矩陣的協(xié)方差矩陣的特征值按大小進(jìn)行排序,以此順序計(jì)算各個(gè)主元及其貢獻(xiàn)率,之后選取滿足要求的主元從而達(dá)到降維的目的。目前PCA已廣泛用于特征提取、數(shù)據(jù)降噪以及故障重構(gòu)等方面。
PCA的基本原理如下:設(shè)原始數(shù)據(jù)樣本包含m個(gè)n維特征向量 Xk=( x1,x2,… xn)T,則PCA具體步驟如下:
計(jì)算樣本的平均值:
求出樣本向量的協(xié)方差矩陣:
計(jì)算S的特征值 λ1,λ2,…,λn及其對(duì)應(yīng)的特征向量 v1,v2,…,vn。
計(jì)算各主元的貢獻(xiàn)率和累計(jì)貢獻(xiàn)率:
選取前p個(gè)累計(jì)貢獻(xiàn)率滿足要求的主元作為映射矩陣,將原始信號(hào)變換到特征空間:
式(5)中,V=(v1,v2,…vp)。
SPRT的基本原理為:設(shè) x1,x2,…,xn是一組獨(dú)立同分布的隨機(jī)變量序列。對(duì)于這樣的總體樣本值,提出兩個(gè)假設(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í),對(duì)兩種假設(shè)的接受與否的判斷顯得不夠合理。而Wald改進(jìn)的地方是當(dāng) λn遠(yuǎn)小于C時(shí)接受 H0,當(dāng)λn遠(yuǎn)大于C時(shí)拒絕 H0,當(dāng) λn與C很接近時(shí)暫且不做出判斷,而是再繼續(xù)檢測(cè)λn+1,直到之后的似然比遠(yuǎn)大于或遠(yuǎn)小于C為止。
在SPRT中,事先給定兩類錯(cuò)誤的概率α與 β的值,則閾值A(chǔ)、B可以按照下列公式來(lái)確定:
檢驗(yàn)判斷準(zhǔn)則為:代入觀測(cè)值x1,計(jì)算λ1(x),若 λ1(x)<B,則停止觀測(cè)并接受零假設(shè) H0;若λ1(x)>A,則停止觀測(cè)并接受備選假設(shè) H1;若A≤λ1(x)≤B,則繼續(xù)進(jìn)行檢驗(yàn),直至似然比值滿足上述兩種情況為止。
選取正常齒輪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之間會(huì)產(chǎn)生振動(dòng),因此可以任選其中之一來(lái)模擬故障。本次實(shí)驗(yàn)選擇用齒輪3來(lái)進(jìn)行模擬實(shí)驗(yàn)。
圖1 齒輪箱結(jié)構(gòu)圖Fig.1 Structure of gearbox
齒輪箱信號(hào)采集系統(tǒng)如圖2所示,將加速度傳感器安裝在齒輪箱的水平垂直兩個(gè)方向上,通過(guò)動(dòng)態(tài)模擬器采集振動(dòng)信號(hào),并通過(guò)信號(hào)分析器將信號(hào)輸入并存儲(chǔ)在計(jì)算機(jī)中,本文僅研究分析水平方向的信號(hào)。分別用S1、S2和S3表示正常情況與兩種故障情況下的信號(hào)。
圖2 齒輪箱信號(hào)采集系統(tǒng)Fig.2 Signal acquisition system of gearbox
小波包變換將同時(shí)分解信號(hào)的低頻和高頻部分,這有利于有效地區(qū)分信號(hào)中的噪聲。對(duì)實(shí)驗(yàn)信號(hào)進(jìn)行3層小波包分解,得到8個(gè)頻帶的分解信號(hào)后再重構(gòu)信號(hào)。這樣重構(gòu)后的信號(hào)會(huì)更準(zhǔn)確。
振動(dòng)信號(hào)的特征信息是通過(guò)一些特征參數(shù)來(lái)表達(dá)的,每個(gè)特征參數(shù)都能反映該信號(hào)的某種特性,平穩(wěn)的振動(dòng)信號(hào)表明齒輪箱運(yùn)行正常,而異常的或是有較大波動(dòng)的振動(dòng)信號(hào)則表明齒輪箱出現(xiàn)了故障。在常用的時(shí)域特征參數(shù)中,均值能反映信號(hào)波動(dòng)能量的大??;有效值是描述動(dòng)態(tài)信號(hào)強(qiáng)度的指標(biāo);標(biāo)準(zhǔn)差是描述信號(hào)的波動(dòng)分量;峭度指標(biāo)能反應(yīng)振動(dòng)信號(hào)中的沖擊特征;波形指標(biāo)是反映信號(hào)偏離高斯分布程度的指標(biāo);脈沖指標(biāo)與峰值指標(biāo)都是用于檢測(cè)信號(hào)中有無(wú)沖擊的指標(biāo);裕度指標(biāo)常被用來(lái)檢測(cè)機(jī)械設(shè)備的磨損狀況[13-15]。本文選取這8個(gè)參數(shù)指標(biāo)進(jìn)行計(jì)算。假設(shè)xi=[x1,x2,…,xN]是一組離散待檢信號(hào),N=8 192。每組取1 024個(gè)檢驗(yàn)點(diǎn),這樣可以得到7196組檢驗(yàn)數(shù)據(jù)。
對(duì)實(shí)驗(yàn)信號(hào)進(jìn)行降噪之后再提取特征值,本文中n=1 024,特征參數(shù)的計(jì)算公式如下:
將經(jīng)過(guò)預(yù)處理的振動(dòng)信號(hào)代入公式中分別計(jì)算,得到齒輪箱在3種狀況下的3組信號(hào)的8種特征參數(shù)。
對(duì)上述3組數(shù)據(jù)進(jìn)行PCA處理,根據(jù)得到的結(jié)果來(lái)選取待檢驗(yàn)的序列。
經(jīng)過(guò)PCA處理之后,根據(jù)貢獻(xiàn)率的高低來(lái)對(duì)處理后的主元進(jìn)行降序排列,排序越靠前的主元保留的原有信息越多。在得到的處理結(jié)果中,通常是選取累計(jì)貢獻(xiàn)率高于85%的前幾個(gè)主元,這樣選取的主元才能包含足夠多的原有信息。
這些特征值經(jīng)過(guò)PCA處理之后,前4個(gè)主元就能得到超過(guò)85%的信息。這里選取貢獻(xiàn)率最大的序列作為檢驗(yàn)序列,記作 yi=[y1…yk]。其均值和標(biāo)準(zhǔn)差為:
經(jīng)過(guò)以上處理之后的待檢序列大致符合正態(tài)分布。正常狀態(tài)下,該序列滿足零假設(shè)H0∶θ=θ0;故障狀態(tài)下,該序列滿足備擇假設(shè)H1∶θ=θ1。標(biāo)準(zhǔn)差保持為σ。當(dāng)兩個(gè)假設(shè)都成立時(shí),該序列的聯(lián)合概密度為:
式(20)和(21)中:P0i為零假設(shè)條件下的概率密度函數(shù);P1i為備擇假設(shè)條件下的概率密度函數(shù)。
SPRT似然比為:
式(22)中:P0為零假設(shè)條件下的先驗(yàn)概率;P1為備擇假設(shè)條件下的先驗(yàn)概率。
對(duì)似然比公式進(jìn)一步推導(dǎo)簡(jiǎn)化可使得計(jì)算更加簡(jiǎn)便:
簡(jiǎn)化之后的閾值a=lnA,b=lnB。假設(shè)檢驗(yàn)范兩類錯(cuò)誤的概率分別取α=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ù)比較下一個(gè)似然比值與a、b的關(guān)系,如此循環(huán),直至得出結(jié)果。運(yùn)用此方法對(duì)齒輪裂紋進(jìn)行診斷的完整流程如圖3所示。
圖3 齒輪裂紋診斷流程圖Fig.3 Flowchart of gear crack fault diagnosis
利用均方根誤差(Root Mean Square Error,RMSE)來(lái)檢驗(yàn)SPRT應(yīng)用于此的可靠性。均方根誤差計(jì)算公式如下:
對(duì)3組不同的齒輪分別進(jìn)行實(shí)驗(yàn),得到圖4(a)中的原始振動(dòng)信號(hào)。其中S1、S2和S3分別為正常齒輪、25%和50%裂紋齒輪的振動(dòng)信號(hào)。3組原始振動(dòng)信號(hào)經(jīng)過(guò)降噪之后,得到圖4(b)中的信號(hào)。
對(duì)比圖 4(a)和圖 4(b)可以看到,與原始振動(dòng)信號(hào)相比,降噪后信號(hào)的幅值明顯有所減小,而且信號(hào)曲線也更為清晰。圖4(b)基本達(dá)到了預(yù)期效果,這也表明了小波包變換能很好地對(duì)原始振動(dòng)信號(hào)進(jìn)行降噪。
對(duì)3組實(shí)驗(yàn)信號(hào)進(jìn)行特征提取之后再用主元分析法處理,發(fā)現(xiàn)前4個(gè)主元包含超過(guò)85%的有效信息,信號(hào)S1、S2和S3的各主元貢獻(xiàn)率排序如表2所示。由PCA的處理結(jié)果可知這3組信號(hào)的第一主元的貢獻(xiàn)率分別達(dá)到了50.9%、50.2%、51.8%。由于排名第一的主元所占的比例最大,因此被選用來(lái)進(jìn)行SPRT。
圖4 實(shí)驗(yàn)信號(hào):(a)原始振動(dòng)信號(hào),(b)去噪后的信號(hào)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 %
由似然比計(jì)算公式可知,平均值的變化對(duì)似然比影響較大。將正常信號(hào)S1的平均值記為μ0,故障信號(hào)S2和S3的平均值記為μ1。對(duì)信號(hào)進(jìn)行檢測(cè),結(jié)果如圖5所示。
當(dāng)輸入的信號(hào)是 S1和S2時(shí),檢驗(yàn)結(jié)果如圖5(a)所示,從圖中可以判定,信號(hào)S1的似然比Δ1,2(S1)<b,即齒輪箱正常工作;信號(hào) S2的似然比Δ1,2(S2)>a,即齒輪箱出現(xiàn)故障。當(dāng)輸入的信號(hào)是S1和S3時(shí),檢驗(yàn)結(jié)果如圖5(a)所示,從圖中可以判定,信號(hào)S1的似然比Δ1,3(S1)<b,即齒輪箱正常工作;信號(hào)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)不同程度裂紋的齒輪也可以通過(guò)這種方法加以區(qū)分,其方法基本相同。就本次實(shí)驗(yàn)而言,將兩種故障信號(hào)S2和S3進(jìn)行檢測(cè),檢驗(yàn)結(jié)果如圖5(b)所示,從圖中可以判定,信號(hào)S2的似然比Δ2,3(S2)<b,即齒輪具有輕度損傷;信號(hào)S3的似然比Δ2,3(S3)>a,即齒輪出現(xiàn)重度故障。所以根據(jù)檢驗(yàn)結(jié)果可以對(duì)25%裂紋和50%裂紋的信號(hào)進(jìn)行清晰地區(qū)分。
在實(shí)驗(yàn)中,正常狀況下獲取兩組數(shù)據(jù),分別進(jìn)行降噪與PCA處理后共得到14 338個(gè)數(shù)據(jù)點(diǎn)。將這些數(shù)據(jù)分為16組,分別記為A1至A16。對(duì)16組數(shù)據(jù)分別進(jìn)行序貫概率比檢驗(yàn),得到16個(gè)似然比值,記為B1至B16。計(jì)算這16個(gè)值之間的均方根誤差,即B1分別其他15個(gè)值比較,B2分別與其他14個(gè)值比較,…,B15與B16比較,得到15組比值,選取每組比值中的最大值形成誤差曲線。將故障狀況下的似然比記為a,計(jì)算故障狀況與正常狀況之間的誤差,即a與B1,B2…B16比較,將得到16個(gè)誤差值形成誤差曲線。通過(guò)比較誤差來(lái)驗(yàn)證序貫概率比檢驗(yàn)診斷齒輪裂紋故障的可靠性。
將正常信號(hào)S1的平均值記為 μ0,故障信號(hào)S2的平均值記為 μ1。將16組正常數(shù)據(jù)信號(hào)的似然比Δ1,2(S1)輸入程序中,得到了 16 條 normal曲線;將2組故障信號(hào)的似然比Δ1,2(S2),Δ1,2(S3)輸入程序中,得到了2條fault曲線,最后的結(jié)果如圖6(a)所示。圖6(a)表明16組正常曲線的趨勢(shì)基本相同,而故障信號(hào)的曲線與之相比則有較大的區(qū)別。計(jì)算16組正常信號(hào)的似然比之間的均方根誤差E1,2(S1),每組誤差中的最大值形成的誤差曲線如圖6(b)所示。分別計(jì)算故障信號(hào)S1和S2與16組正常信號(hào)似然比之間的誤差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)
將正常信號(hào)S1的平均值記為 μ0,故障信號(hào)S3的平均值記為 μ1。將16組正常數(shù)據(jù)信號(hào)的似然比Δ1,3(S1)輸入程序中,得到了 16 條 normal曲線;將兩組故障信號(hào)的似然比Δ1,3(S2),Δ1,3(S3)輸入程序中,得到了2條fault曲線,最后的結(jié)果如圖7(a)所示。圖7(a)說(shuō)明16組正常曲線的趨勢(shì)基本相同,而故障信號(hào)的曲線與之相比則有較大的區(qū)別。計(jì)算16組正常信號(hào)的似然比之間的均方根誤差E1,3(S1),每組誤差中的最大值形成的誤差曲線如圖7(b)所示。分別計(jì)算故障信號(hào)S1和S2與16組正常信號(hào)似然比之間的誤差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)以正常信號(hào)的平均值為參數(shù)μ0時(shí),不論是以哪種故障信號(hào)的平均值為參數(shù) μ1來(lái)建立模型,從最后的結(jié)果中都能對(duì)3種狀況進(jìn)行清晰地區(qū)分。
本文提出了一種將PCA與SPRT相結(jié)合的檢測(cè)方法,并選取不同程度裂紋的故障齒輪與正常齒輪進(jìn)行模擬試驗(yàn),來(lái)檢測(cè)此方法能否有效區(qū)分不同狀態(tài)下的齒輪。在預(yù)處理之后利用PCA對(duì)實(shí)驗(yàn)信號(hào)所選取的特征參數(shù)進(jìn)行降維處理,有效地減少了特征數(shù)量,使得之后的檢驗(yàn)更加高效準(zhǔn)確,最終有效地識(shí)別出了齒輪箱不同的狀態(tài)。此外,經(jīng)過(guò)多次實(shí)驗(yàn)后證明在其他的工作條件下,該方法也能對(duì)齒輪箱不同的狀態(tài)進(jìn)行清晰地區(qū)分。因此將該方法應(yīng)用于齒輪裂紋故障檢測(cè)與分類方面是有效而且可靠的。此方法也可供其他旋轉(zhuǎn)機(jī)械的診斷研究進(jìn)行參考。