一種基于增強型最大二階循環(huán)平穩(wěn)盲解卷積的齒輪箱復(fù)合故障診斷

齊詠生 單成成 賈舜宇 劉利強 董朝軼

1.內(nèi)蒙古工業(yè)大學(xué)電力學(xué)院，呼和浩特，0100802.內(nèi)蒙古自治區(qū)機電控制重點實驗室，呼和浩特，010051

0 引言

齒輪箱工作環(huán)境惡劣，內(nèi)部結(jié)構(gòu)復(fù)雜，當(dāng)其中某部位發(fā)生局部損傷時，產(chǎn)生非平穩(wěn)的周期性沖擊信號會激起齒輪及其相鄰部件的共振，常常會出現(xiàn)多種故障并存的復(fù)合故障狀態(tài)，而復(fù)合故障特征之間存在相互干擾、交叉混疊，且常被干擾噪聲所淹沒，難以提取準(zhǔn)確的特征信息[1]，因此，如何從齒輪箱復(fù)合故障信號中分離出各故障特征仍然是一個亟待解決的難題[2-3]。

為有效分離提取齒輪箱復(fù)合故障特征，研究人員提出了一些診斷的方法，如經(jīng)驗?zāi)B(tài)分解[4](EMD)、集合經(jīng)驗?zāi)B(tài)分解[5](EEMD)、局部均值分解[6](LMD)等，盡管上述方法在齒輪箱復(fù)合故障診斷中取得了一定的成效，但在干擾噪聲的影響下，仍面臨較大挑戰(zhàn)?？紤]到齒輪箱復(fù)合故障振動信號傳遞路徑復(fù)雜多變，可以看作是源信號與信道線性卷積混合的結(jié)果，因此提取故障源信號的沖擊特征便可理解為一個解卷積過程。為此，WIGGINS[7]提出最小熵解卷積(minimum entropy deconvolution, MED)算法，旨在消除傳輸路徑的影響，以峭度值最大化為目標(biāo)，實現(xiàn)對原始沖擊信號的恢復(fù)。冷軍發(fā)等[8]利用MED對齒輪箱振動信號進行解卷積濾波降噪，然后對濾波后的信號作包絡(luò)解調(diào)分析提取故障信號的沖擊特征，盡管該方法實現(xiàn)了齒輪箱早期故障診斷，但MED算法僅能提取少數(shù)突出的尖脈沖，而且存在一些偽脈沖成分[9]。為了更好地從故障振動信號中提取出周期性脈沖信號，MCDONALD等[10]提出最大相關(guān)峭度解卷積(maximum correlate kurtosis deconvolution, MCKD)算法，該方法能夠從振動故障信號中分離出周期性的沖擊特征，并已在旋轉(zhuǎn)機械復(fù)合故障診斷中得到有效應(yīng)用[11]。但MCKD的降噪能力受濾波器長度、解卷積周期等參數(shù)的共同影響，濾波器長度的增大會嚴(yán)重影響其運行速度，且周期為非整數(shù)時需要重采樣或?qū)ζ鋱A整?？傮w來看，只有當(dāng)所有參數(shù)均設(shè)置合適時，MCKD才能發(fā)揮其在故障診斷中的優(yōu)越性[12]。此外，MCKD相比于MED能夠提取更多沖擊脈沖個數(shù)，但也只限于提取局部有限個沖擊。

考慮到齒輪振動信號具有循環(huán)平穩(wěn)特性，近年來，BUZZONI等[13]提出了最大二階循環(huán)平穩(wěn)盲解卷積(maximum second-order cyclostationarity blind deconvolution, CYCBD)算法，該方法能夠提取出包含旋轉(zhuǎn)機械故障沖擊特征的源信號，適用于齒輪局部損傷的故障診斷。CYCBD以最大二階循環(huán)平穩(wěn)指標(biāo)(second-order indicators of cyclostationary,ICS2)為依據(jù)，能夠有效地從振動信號中提取出故障沖擊成分，但CYCBD算法中需預(yù)先根據(jù)先驗知識設(shè)定循環(huán)頻率集，且有效濾波器的長度也必須根據(jù)經(jīng)驗選擇合適，才有可能取得較好的診斷效果，條件比較苛刻，這使得CYCBD算法缺乏自適應(yīng)性而難以應(yīng)用推廣。

鑒于上述情況，本文針對齒輪箱復(fù)合故障難以有效提取的問題，提出一種自適應(yīng)CYCBD(adaptive CYCBD,ACYCBD)與1.5維導(dǎo)數(shù)增強譜相結(jié)合的復(fù)合故障特征提取方法。

1 CYCBD方法

一般來說，盲解卷積旨在從含有噪聲的振動信號X中分離提取出故障沖擊信號y，寫成矩陣形式為

(1)

式中，g為逆濾波器系數(shù)矩陣；y[]、x[]、g[]的中括號表示其為離散信號；N為濾波器長度值；L為振動離散信號的長度。

二階循環(huán)平穩(wěn)性的一般表達(dá)式為

(2)

(3)

(4)

式中，k取正整數(shù)，表示循環(huán)頻率的整數(shù)倍；Ts為故障沖擊的周期，本文將離散信號的循環(huán)頻率定義為1/Ts，此時，α稱為循環(huán)頻率集。

將IICS2寫成矩陣表達(dá)式為

(5)

式中，上標(biāo)H表示矩陣的共軛轉(zhuǎn)置;W為加權(quán)矩陣。

由式(1)和式(5)可以得出：

(6)

式中，RXWX、RXX分別為加權(quán)相關(guān)矩陣和相關(guān)矩陣。

式(6)被稱為廣義Rayleigh熵，由其性質(zhì)可知，要求IICS2最大值，即轉(zhuǎn)換成求解廣義Rayleigh熵的最大特征值λ，即

RXWXg=RXXgλ

(7)

2 自適應(yīng)CYCBD方法(ACYCBD)

2.1 循環(huán)頻率集的選取方法

當(dāng)齒輪箱發(fā)生局部故障時，其信號具有調(diào)幅特性，調(diào)幅信號的二階統(tǒng)計特征具有周期性，因此齒輪箱的故障信號為二階循環(huán)平穩(wěn)信號，其二階循環(huán)統(tǒng)計量Rx(t,τ)表示為

Rx(t,τ)=Rx(t+T0,τ)

(8)

式中，T0、τ分別為周期和延遲量。

式(8)可用傅里葉級數(shù)進行展開：

(9)

其中，Rx(α,τ)為傅里葉系數(shù)，也稱其為循環(huán)自相關(guān)函數(shù)。對于循環(huán)平穩(wěn)信號，其循環(huán)譜密度函數(shù)與循環(huán)自相關(guān)函數(shù)是一對傅里葉變換[14]，因此，循環(huán)譜密度函數(shù)可以表示為

(10)

這里將循環(huán)譜密度函數(shù)S(α,f)簡稱為循環(huán)譜。通過循環(huán)譜分析后，可以將一維振動信號x(t)轉(zhuǎn)換為橫坐標(biāo)為循環(huán)頻率α、縱坐標(biāo)為載波頻率f的二維譜圖。

由于CYCBD方法在旋轉(zhuǎn)機械復(fù)合故障診斷中循環(huán)頻率大都預(yù)先采用公式進行計算，然后在已知條件下進行故障特征的提取,而實際工程應(yīng)用中，傳感器采集的振動信號事先并不知道其中是否存在故障，利用傳統(tǒng)方法進行故障診斷時必須將可能存在的故障情況逐一排查，較為繁瑣,而且，由于理論計算和實際情況常常存在較大誤差，很容易造成誤診和漏診,因此，本文引入循環(huán)譜分析來解決上述問題，通過循環(huán)譜分析可自動識別振動信號中與故障特征相關(guān)的循環(huán)頻率成分，無需預(yù)先采用公式進行計算，避免循環(huán)頻率集選擇上的先驗性和盲目性，使CYCBD 在接下來故障提取過程中具有更好的自適應(yīng)性和針對性。

2.2 濾波器長度的自適應(yīng)選取方法

濾波器長度過大會使計算量增大，導(dǎo)致故障診斷時間過長，濾波器長度過小，CYCBD算法不能起到合理的降噪效果，而較優(yōu)濾波器長度的選取往往依靠經(jīng)驗或者根據(jù)預(yù)先確定的指標(biāo)在一定濾波器長度區(qū)間內(nèi)自適應(yīng)選擇，為此，本文需要確定恰當(dāng)?shù)闹笜?biāo)對濾波器長度進行自適應(yīng)選擇。

當(dāng)振動信號的幅值分布比較均勻、相差不大時，信號的稀疏性會表現(xiàn)出較為微弱的狀態(tài)，這時稀疏度值較小。一旦信號中出現(xiàn)比較明顯的周期性連續(xù)沖擊，信號呈現(xiàn)出較強的稀疏特性，此刻稀疏度指標(biāo)也隨之增大，但齒輪箱振動信號往往受隨機沖擊、干擾噪聲等因素的干擾，會對稀疏度結(jié)果產(chǎn)生影響。前期的研究發(fā)現(xiàn)，高階累積量可以用來測量時間序列的高階互相關(guān)，同時又可體現(xiàn)出隨機過程相對高斯分布的偏移程度，且高斯噪聲的高階累積量為零，因而高階累積量具有抑制噪聲的作用。為此，本文提出采用三階累積量稀疏度(third-order cumulative sparsity，TCS)指標(biāo)對濾波器長度進行自適應(yīng)選取，該指標(biāo)定義如下：

(11)

式中，cn為xn的三階累積量。

sTCS值可以綜合考量故障沖擊信息與環(huán)境噪聲的影響，sTCS越大，說明濾波信號中包含的瞬態(tài)周期性沖擊成分越明顯?？紤]到濾波器的長度值若太小，則ACYCBD不能達(dá)到有效的濾波效果，若濾波器長度值太大，則會增加計算成本且可能導(dǎo)致降噪后的信號失去原始信號的特性等因素，為了兼顧ACYCBD算法的有效性和計算效率，濾波器的長度應(yīng)該能夠基本覆蓋齒輪故障所產(chǎn)生沖擊信號的衰減周期，且長度不宜過大。衰減周期記為D，則D=fs/α*, 其中，fs、α*分別為采樣頻率與循環(huán)頻率值。鑒于此，本文將濾波器長度的搜索范圍設(shè)定在0.5D～1.5D之間，以TCS指標(biāo)最大化為原則，選取該搜索區(qū)間內(nèi)TCS曲線極大值點處對應(yīng)的濾波器長度值作為最優(yōu)值。

2.3 一種新的復(fù)合式1.5維導(dǎo)數(shù)增強譜的提出

由于背景噪聲的存在，利用ACYCBD分離出的最優(yōu)濾波信號仍會受到噪聲的一定干擾，導(dǎo)致故障沖擊特征不夠明晰，診斷效果不能達(dá)到最佳，為此，本文提出了一種新的1.5維導(dǎo)數(shù)增強譜。

包絡(luò)導(dǎo)數(shù)能量算子是一種非線性差分算子[15]，它通過求導(dǎo)引入瞬時頻率的權(quán)重，對噪聲有著較好的抗干擾性，并且計算信號導(dǎo)數(shù)的包絡(luò)來追蹤信號的瞬態(tài)能量特征，能夠?qū)X輪箱的故障沖擊成分進行增強。

同時考慮到1.5維譜能夠有效地抑制高斯白噪聲，強化基頻成分，提高信噪比，是分析非線性、非平穩(wěn)信號的有力工具[16]，為了最大限度地降低干擾噪聲對ACYCBD濾波信號的影響，本文提出將包絡(luò)導(dǎo)數(shù)能量算子與1.5維譜進行融合，形成一種新的1.5維導(dǎo)數(shù)能量算子，將其稱為1.5維導(dǎo)數(shù)(增強)譜。

(1)對于一般的解調(diào)信號，通常使用解析信號的形式來對其進行定義，即

X(t)=x(t)+jH[x(t)]=Aejφ(t)

(12)

式中，H[x(t)]為信號x(t)的希爾伯特變換；A(t)為信號的瞬時幅值；φ(t)為信號的瞬時頻率。

(13)

令x(t)=Acos(ωt+φ)，則

(14)

(15)

將式(14)和式(15)代入式(13)得出包絡(luò)導(dǎo)數(shù)能量算子的時域、頻域表達(dá)式如下：

(16)

式中，*表示卷積運算。

(2)設(shè)解調(diào)后的信號為

(17)

則1.5維導(dǎo)數(shù)譜定義如下：

(18)

其中，Y(ω)為y(t)的傅里葉變換，即

(19)

式中，δ函數(shù)為沖擊函數(shù)。

將式(19)代入式(18)中化簡可得

(20)

由沖擊函數(shù)的性質(zhì)可知，當(dāng)ω=ωk時，式(20)不為零。且對于沖擊函數(shù)，有如下等式成立：

δ(ω-ωu)*δ(ω-ωv)=δ(ω-ωu-ωv)

(21)

因此，可進一步得出1.5維導(dǎo)數(shù)增強譜的表達(dá)式為

(22)

式(22)表明，當(dāng)ωk-ωu-ωv=0時，即ωk=ωu+ωv時,δ(ωk-ωu-ωv)=1，此時譜圖上有幅值，顯示頻率耦合之和的分量。當(dāng)ωk=ωu-ωv或ωk=ωv-ωu時，所對應(yīng)的沖擊函數(shù)幅值為1，譜圖上顯示頻率耦合之差的分量。即融合后的1.5維導(dǎo)數(shù)增強譜在計算信號的瞬態(tài)能量時，通過導(dǎo)數(shù)引入信號的瞬時頻率，克服了包絡(luò)譜的能量僅由幅值決定而忽略了瞬時頻率信息的問題。同時，1.5維導(dǎo)數(shù)增強譜能夠有針對性地保留或增強二次頻率耦合成分，剔除或抑制非耦合的頻率分量。進一步研究發(fā)現(xiàn)，當(dāng)齒輪受到局部損傷時，它產(chǎn)生的故障振動信號為多分量調(diào)幅調(diào)頻信號，通過解調(diào)得到包含故障特征頻率及其若干倍頻的調(diào)制成分，而這些特征頻率及其倍頻之間存在二次頻率耦合關(guān)系，因此利用1.5維導(dǎo)數(shù)譜可有效增強故障特征頻率的譜峰。對于隨機脈沖，背景白噪聲等干擾信號多數(shù)屬于非二次耦合的頻率分量，在導(dǎo)數(shù)增強譜圖中會很容易被剔除或抑制。

綜上，當(dāng)齒輪箱發(fā)生局部故障損傷時，各部件之間摩擦碰撞可使振動過程的能量瞬間發(fā)生改變，利用1.5維導(dǎo)數(shù)增強譜可以剔除與故障特征頻率之間無耦合關(guān)系的噪聲分量，計算得到信號的瞬態(tài)能量能夠更加全面地表征振動信號中周期性故障沖擊成分，因此融合后的1.5維導(dǎo)數(shù)增強譜能夠有效抑制故障非沖擊成分，增強故障沖擊成分。

3 基于ACYCBD的復(fù)合故障診斷

當(dāng)齒輪發(fā)生局部損傷如出現(xiàn)裂紋、點蝕、斷齒等故障時，產(chǎn)生的故障特征往往蘊含在瞬態(tài)沖擊信號中，因此，根據(jù)本文所提方法，首先利用循環(huán)譜分析自適應(yīng)檢測出齒輪箱復(fù)合故障振動信號中與故障特征相關(guān)的循環(huán)頻率值，構(gòu)建不同目標(biāo)類型的循環(huán)頻率集；然后，針對各個故障目標(biāo)類型的循環(huán)頻率集，以TCS值(sTCS)最大化為原則，自適應(yīng)選取CYCBD算法的最優(yōu)濾波器長度值;接下來利用該最優(yōu)參數(shù)，得到包含不同故障沖擊成分的最優(yōu)CYCBD濾波信號。為了進一步消除噪聲因素的干擾，增強故障沖擊特征，最后采用提出的1.5維導(dǎo)數(shù)增強譜分析完成故障診斷。該方法總體流程如圖1所示。具體步驟如下：

(1)使用振動傳感器采集齒輪箱復(fù)合故障信號，并通過循環(huán)譜檢測與故障相關(guān)的循環(huán)頻率，構(gòu)建不同目標(biāo)類型的循環(huán)頻率集。

(2)根據(jù)步驟(1)中不同類型的循環(huán)頻率集，分別以TCS值為指標(biāo)，自適應(yīng)選取最優(yōu)濾波器長度。

(3)基于最優(yōu)濾波器長度，得到不同目標(biāo)類型下的CYCBD最優(yōu)濾波信號。

(4) 對分離出的各個信號分別采用1.5維導(dǎo)數(shù)譜進行特征增強，提高頻譜的可讀性。

(5) 通過分析譜圖中的主導(dǎo)故障特征頻率及其倍頻完成復(fù)合故障診斷。

圖1 本文算法的流程圖

4 仿真信號分析

為了驗證本文所提算法在齒輪箱復(fù)合故障診斷方面的有效性，本節(jié)構(gòu)造了包含齒輪故障產(chǎn)生的周期性沖擊成分，以及諧波信號和隨機噪聲等干擾成分的齒輪箱復(fù)合故障仿真信號來進行算法驗證，仿真信號模型如下：

y(t)=x1(t)+x2(t)+B(t)+n(t)

(23)

式中，y(t)為復(fù)合故障信號；x1(t)和x2(t)為齒輪故障的瞬態(tài)沖擊成分，其故障頻率分別設(shè)定為f1=55 Hz ,f2=90 Hz；B(t)為用于模擬信號中存在的諧波成分；n(t)為隨機噪聲。

仿真過程中，設(shè)置采樣頻率為2.56 kHz，采樣點數(shù)為6000。圖2為仿真信號的時域波形圖，受到諧波信號和隨機噪聲因素的影響，從圖2d復(fù)合故障仿真信號中很難發(fā)現(xiàn)明顯的周期性故障沖擊成分，無法給出有效的故障成因。

利用本文所提方法對仿真信號進行故障診斷，首先對復(fù)合故障信號進行循環(huán)譜分析，結(jié)果如圖3所示。分析可知,在循環(huán)頻率α為55 Hz、90 Hz、110 Hz、165 Hz以及180 Hz處循環(huán)譜密度相對突顯，這表明該仿真信號中存在兩類循環(huán)頻率成分，即分別為以55 Hz及其倍頻構(gòu)成的循環(huán)頻率集(記為A1)，以及90 Hz及其倍頻組成的循環(huán)頻率集(記為A2)。盡管利用循環(huán)譜分析可以初步得出各類故障的局部有限個循環(huán)頻率值，但對復(fù)合故障來說，從圖3中可以觀察到兩類故障循環(huán)頻率之間仍有相互交叉，僅利用循環(huán)譜分析無法將不同故障進行有效分離，因此需要進一步使用ACYCBD對不同故障特征進行分離提取。

(a)沖擊信號1(f1=55 Hz)

(b)沖擊信號2(f2=90 Hz)

(c)諧波信號

(d)復(fù)合故障仿真信號圖2 仿真信號的時域波形圖

圖3 仿真信號的循環(huán)譜分析結(jié)果

由上述分析可知，不同故障類型的循環(huán)頻率值有所不同，而ACYCBD方法正是根據(jù)不同目標(biāo)類型的循環(huán)頻率集，通過選取合適的濾波器長度來提取分離出不同的故障沖擊信號的。首先設(shè)定好循環(huán)頻率集A1，此時可計算得出待提取沖擊信號的衰減周期D為465，濾波器長度搜索區(qū)間為0.5D～1.5D，即[230, 700]。再利用TCS指標(biāo)最大化為原則去篩選CYCBD的有效濾波器長度值，TCS曲線如圖4a所示。不難看出在該搜索區(qū)間內(nèi)TCS峰值處濾波器取值為398，因此將398作為最優(yōu)濾波器長度值對復(fù)合故障仿真信號進行CYCBD運算來提取分離相關(guān)的故障沖擊成分。同理，當(dāng)目標(biāo)類型的循環(huán)頻率集為A2時，濾波器長度的搜索范圍設(shè)定為[140, 430]，TCS曲線如圖4b所示，由圖4b可知，在運用CYCBD提取另一種沖擊成分時，濾波器長度值設(shè)置為286。

(a)循環(huán)頻率集為A1時TCS值曲線

(b)循環(huán)頻率集為A2時TCS值曲線圖4 仿真信號的TCS值曲線

圖5所示為循環(huán)頻率集為A1，濾波器長度取398情況下的ACYCBD盲解卷積后的分析結(jié)果，通過對比可以發(fā)現(xiàn)，圖5a包絡(luò)譜僅提取出若干f1的倍頻，且周圍干擾譜線較多，診斷結(jié)果不夠明晰。同樣在圖5b所示的1.5維譜中f1的部分倍頻被噪聲所淹沒，頻譜可讀性欠佳。在圖5c所示的1.5維導(dǎo)數(shù)譜中f1及其全部倍頻清晰明了，幅值突顯，能夠更充分地說明發(fā)生了與f1故障頻率相關(guān)的故障類型。

通過調(diào)整目標(biāo)故障的循環(huán)頻率集為A2，將濾波器長度最優(yōu)值設(shè)定為286，得到ACYCBD盲解卷積結(jié)果如圖6所示。不難看出盡管圖6a所示包絡(luò)譜中f2的倍頻均有體現(xiàn)，但受到周圍干擾譜線的影響較為嚴(yán)重。而圖6b所示為1.5維譜，圖中f2倍頻不夠明顯，無法給出確切的故障成因。由圖6c所示 1.5維導(dǎo)數(shù)譜可知，f2故障特征頻率及其倍頻處峰值明顯增強，干擾譜線較少，可有效診斷出與f2故障頻率相關(guān)的故障類型。

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖5 循環(huán)頻率集為A1時ACYCBD盲解卷積結(jié)果

為了進一步表明本文所提方法在齒輪箱復(fù)合故障診斷方面的優(yōu)越性，與最小熵解卷積(MED)結(jié)合1.5維導(dǎo)數(shù)譜方法、最大相關(guān)峭度解卷積(MCKD)結(jié)合1.5維導(dǎo)數(shù)譜方法，以及局部均值分解(LMD)結(jié)合1.5維導(dǎo)數(shù)譜方法進行對比分析。用MED對仿真信號作解卷積處理，降噪后信號的1.5維導(dǎo)數(shù)譜如圖7a所示。分析發(fā)現(xiàn)f1與f2故障特征頻率互相干擾，并未實現(xiàn)故障分離。利用MCKD方法對仿真信號進行處理，利用公式T=fs/fi(fs為采樣頻率，fi為故障特征頻率)計算得出不同故障類型的解卷積周期值，進而針對性地進行故障特征提取，提取出的f1故障特征結(jié)果如圖7b所示，可以看到僅提取出少量與f1相關(guān)的特征頻率值，診斷效果不佳。圖7c所示為f2故障特征頻率及其若干倍頻的分離結(jié)果，但圖7c中也發(fā)現(xiàn)f1的相關(guān)倍頻，對分析具體的故障成因造成干擾。

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖6 循環(huán)頻率集為A2時ACYCBD盲解卷積結(jié)果

相比EMD、EEMD方法，LMD具有迭代次數(shù)少，計算速度快，可以更好地抑制端點效應(yīng)等優(yōu)勢[17]。利用LMD對仿真信號進行分解，本文選擇與原始信號相關(guān)系數(shù)較大的前3個PF(productfunction)分量進行分析，其余分量多為噪聲等干擾成分，分析結(jié)果如圖8所示。圖中由上至下依次為PF1～PF3分量的時域波形及其1.5維導(dǎo)數(shù)譜。由圖8b可以看出，在PF1分量中主要將f2特征頻率及其倍頻提取出來，但也存在幅值較弱的f1頻率成分。而PF2、PF3中僅分離出少量與f2相關(guān)的頻率成分，且干擾譜線較多，難以準(zhǔn)確診斷出故障成因。

5 試驗平臺驗證

5.1 QPZZ-Ⅱ試驗?zāi)M平臺

本文首先采用來自QPZZ-Ⅱ旋轉(zhuǎn)機械振動分析及故障診斷試驗平臺數(shù)據(jù)進行驗證，該平臺通過更換有缺陷的齒輪進行齒輪故障模擬，如圖9所示，主要由試驗臺基座、變速驅(qū)動電機、齒輪箱、磁粉制動器、聯(lián)軸器構(gòu)成。其中齒輪箱為單級傳動，主動軸上為小齒輪，從動軸上為大齒輪。該系統(tǒng)中大小齒輪均為圓柱型齒輪，大齒輪齒數(shù)為75，小齒輪齒數(shù)為55，模數(shù)均為2 mm。該試驗平臺通過加速度傳感器獲取齒輪的復(fù)合故障振動信號，采集齒輪數(shù)據(jù)時，設(shè)定采樣頻率fs為5120 Hz。

(a)MED處理結(jié)果

(b)MCKD分離結(jié)果(沖擊信號1)

(c)MCKD分離結(jié)果(沖擊信號2)圖7 MED、MCKD結(jié)合1.5維導(dǎo)數(shù)譜對仿真信號的分析結(jié)果

(a)LMD分解后的信號時域圖(PF1～PF3)

(b)分解信號的1.5維導(dǎo)數(shù)譜(PF1～PF3)圖8 LMD結(jié)合1.5維導(dǎo)數(shù)譜對仿真信號的分析結(jié)果

齒輪故障類型的故障特征頻率根據(jù)齒輪參數(shù)可通過下式計算得到：

式中，f1、f2、fm分別為小大齒輪故障特征頻率和嚙合頻率；z1、z2分別為小齒輪的齒數(shù)和大齒輪的齒數(shù)；n為旋轉(zhuǎn)速度。

圖9 QPZZ-Ⅱ機械故障模擬試驗臺

5.1.1大齒輪點蝕和小齒輪磨損復(fù)合故障診斷

通過在大齒輪的某個齒上制造出若干小圓槽以及將小齒輪的某個齒磨平來模擬齒輪箱的點蝕-磨損復(fù)合故障。實測驅(qū)動電機轉(zhuǎn)速為825 r/min，根據(jù)式(24)可得小齒輪與大齒輪的故障特征頻率理論值分別為13.75 Hz和10.08 Hz，齒輪的嚙合頻率為756.25 Hz。

首先，通過傳感器采集復(fù)合故障振動信號，對該信號進行包絡(luò)解調(diào)，得到包絡(luò)譜如圖10所示，可以發(fā)現(xiàn)大齒輪和小齒輪故障特征頻率交叉混疊，頻譜不夠清晰，無法提取確切的故障特征，給診斷帶來一定難度。

圖10 復(fù)合故障振動信號頻域圖

接下來使用本文方法進行分析。根據(jù)齒輪發(fā)生局部損傷時產(chǎn)生的故障沖擊信號具有循環(huán)平穩(wěn)性的特點，針對齒輪箱點蝕-磨損復(fù)合故障振動信號進行循環(huán)譜分析，結(jié)果如圖11所示?？梢园l(fā)現(xiàn)在循環(huán)頻率α為10.1 Hz、13.75 Hz、20.19 Hz以及27.5 Hz處循環(huán)譜密度值較高，這表明振動信號中存在循環(huán)頻率為10.1 Hz及其倍頻的循環(huán)平穩(wěn)信號，其循環(huán)頻率集記作B1，以及循環(huán)頻率為13.75 Hz及其倍頻的循環(huán)平穩(wěn)信號，其循環(huán)頻率集記作B2。

圖11 大齒輪點蝕和小齒輪磨損復(fù)合故障的循環(huán)譜分析結(jié)果

設(shè)定好CYCBD所需要的循環(huán)頻率集B1，利用TCS指標(biāo)來篩選出最優(yōu)濾波器長度值，TCS值曲線如圖12a所示。濾波器長度值為435時，此刻TCS值達(dá)到區(qū)間極大值，即當(dāng)提取相關(guān)故障特征時，濾波器長度設(shè)定為435。同理可獲得圖12b循環(huán)頻率集為B2時TCS值曲線，該情況下濾波器長度值應(yīng)選擇726進行故障特征提取。

(a)循環(huán)頻率集為B1時TCS值曲線

(b)循環(huán)頻率集為B2時TCS值曲線圖12 大齒輪點蝕和小齒輪磨損復(fù)合信號的TCS值曲線

圖13所示為循環(huán)頻率集為B1，濾波器長度取435參數(shù)下的ACYCBD盲解卷積結(jié)果。通過分析發(fā)現(xiàn)，圖13a包絡(luò)譜基本提取了f1故障特征頻率及其若干倍頻，但周圍干擾譜線相對較多，而13b所示1.5維譜則存在f1部分倍頻被噪聲淹沒的情況，在圖13c所示1.5維導(dǎo)數(shù)譜中，可以看出干擾噪聲得到進一步降低，同時f1故障特征頻率及其倍頻處峰值個數(shù)增多，幅值更加突顯，解調(diào)出的特征頻率值與小齒輪磨損故障理論值基本吻合，且周圍譜線影響很小。

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖13 循環(huán)頻率集為B1時ACYCBD盲解卷積結(jié)果

改變目標(biāo)故障的循環(huán)頻率集為B2，濾波器長度取726進行ACYCBD運算來提取相應(yīng)的故障成分，結(jié)果如圖14所示。圖14a包絡(luò)譜存在較多干擾譜線，部分倍頻峰值不夠突出，診斷效果不佳。圖14b所示為1.5維譜分析結(jié)果，由圖14b可知，僅提取出少量f2的倍頻成分，幅值不夠明晰。圖14c為利用1.5維導(dǎo)數(shù)譜特征增強后的譜圖，其倍頻清晰，幅值突顯，解調(diào)出的特征頻率值與大齒輪點蝕故障理論值相一致，由此可說明大齒輪點蝕故障的發(fā)生。

同理，與MED、MCKD方法進行比較分析，結(jié)果如圖15所示。從圖15a所示的MED處理結(jié)果中發(fā)現(xiàn)，f1、f2特征頻率互相干擾，無法有效分離。此外，通過試驗發(fā)現(xiàn)僅憑MED解卷積很難提取出齒輪箱復(fù)合故障特征。圖15b利用MCKD僅分離出局部有限個小齒輪磨損故障的倍頻特征，其余故障信息幾乎被噪聲所掩蓋，由此診斷容易造成誤診或漏診。圖15c分離出的大齒輪點蝕故障干擾譜線相對較多，少量倍頻模糊不清，影響識別效果。盡管MCKD比MED識別效果有所改善，但在存在噪聲等干擾因素的情況下，仍不具備免疫性。

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖14 循環(huán)頻率集為B2時ACYCBD盲解卷積結(jié)果

(a)MED處理結(jié)果

(b)MCKD分離結(jié)果(小齒輪磨損故障)

(c)MCKD分離結(jié)果(大齒輪點蝕故障)圖15 MED、MCKD結(jié)合1.5維導(dǎo)數(shù)譜對齒輪點蝕與磨損信號的分析結(jié)果

圖16所示為利用LMD結(jié)合1.5維導(dǎo)數(shù)譜對齒輪箱復(fù)合故障信號進行分解結(jié)果，同樣選擇前3個PF分量進行分析，由圖16b可以看出，大小齒輪故障特征頻率及其倍頻在PF1和PF2分量中同時被提取出來，并未將其有效分離，且噪聲干擾較為明顯，而PF3分量僅僅提取出小齒輪故障及其二倍特征頻率，同時受到周圍譜線的嚴(yán)重干擾。因此利用LMD方法對齒輪箱復(fù)合故障進行診斷，仍難以克服噪聲等干擾因素的影響，綜上，可進一步表明本文所提方法的有效性和優(yōu)越性。

(a)LMD分解后的信號時域圖(PF1～PF3)

(b)分解信號的1.5維導(dǎo)數(shù)譜(PF1～PF3)圖16 LMD結(jié)合1.5維導(dǎo)數(shù)譜對齒輪點蝕與磨損信號分析結(jié)果

5.1.2大齒輪斷齒和小齒輪磨損復(fù)合故障診斷

采用本文方法對大齒輪斷齒和小齒輪磨損復(fù)合故障進行診斷，實測驅(qū)動電機轉(zhuǎn)速為812 r/min，計算小齒輪與大齒輪故障特征頻率理論值分別為13.53 Hz和9.92 Hz。首先，采用循環(huán)譜分析檢測復(fù)合故障振動信號中可能存在的具有循環(huán)平穩(wěn)性的故障信號，結(jié)果如圖17所示。從圖17中可獲得兩類不同目標(biāo)故障類型的循環(huán)頻率集，分別是以13.56 Hz及其倍頻組成的循環(huán)頻率集，記作C1，以及9.9 Hz及其倍頻組成的循環(huán)頻率集，記作C2。接下來設(shè)定好ACYCBD算法所需要的故障特征循環(huán)頻率集C1、C2，并利用圖18所示TCS曲線選擇出的最佳濾波器長度值分別為542、570進行解卷積運算。由此可得到解卷積結(jié)果分別如圖19、圖20所示。從圖19c所示的1.5維導(dǎo)數(shù)譜中可以清晰發(fā)現(xiàn)小齒輪磨損主導(dǎo)的故障特征頻率及其倍頻處出現(xiàn)明顯的突出峰值成分，以及由圖20c所示的1.5維導(dǎo)數(shù)譜可知，大齒輪故障特征頻率及其倍頻處峰值突出，且解調(diào)效果均優(yōu)于傳統(tǒng)的包絡(luò)解調(diào)和1.5維譜方法。因此，利用本文方法同樣實現(xiàn)了大小齒輪斷齒與磨損復(fù)合故障診斷。

圖17 大齒輪斷齒和小齒輪磨損復(fù)合故障的循環(huán)譜分析結(jié)果

(a)循環(huán)頻率集為C1時TCS值曲線

(b)循環(huán)頻率集為C2時TCS值曲線圖18 大齒輪斷齒和小齒輪磨損復(fù)合信號的TCS值曲線

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖19 循環(huán)頻率集為C1時ACYCBD盲解卷積結(jié)果

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖20 循環(huán)頻率集為C2時ACYCBD盲解卷積結(jié)果

(a)MED處理結(jié)果

(b)MCKD分離結(jié)果(小齒輪磨損故障)

(c)MCKD分離結(jié)果(大齒輪斷齒故障)圖21 MED、MCKD結(jié)合1.5維導(dǎo)數(shù)譜對齒輪斷齒與磨損信號的分析結(jié)果

與MED、MCKD解卷積方法比較結(jié)果如圖21所示。分析圖21a可知，僅能在f1和f2及其少量倍頻處發(fā)現(xiàn)頻率峰值，且幅值不夠突顯，噪聲干擾比較嚴(yán)重，顯然這樣的解卷積結(jié)果無法給出確切的故障成因。圖21b所示為MCKD解卷積后故障信號的1.5維導(dǎo)數(shù)譜，可以看到譜圖中僅提取出小齒輪磨損故障的局部有限個沖擊特征，且譜峰周圍受到非故障特征頻率的干擾。同樣，在圖21c中雖然可以找出f2～4f2特征頻率值，但之后的若干倍頻被噪聲所覆蓋，不能夠充分說明大齒輪斷齒故障的發(fā)生。

圖22所示為LMD分析結(jié)果，從圖22b中可以發(fā)現(xiàn)PF1、PF3分量中僅提取出大小齒輪故障特征頻率的1倍頻，而PF2分量中兩種故障特征頻率的倍頻值均有體現(xiàn)，說明LMD在分解過程中發(fā)生了模態(tài)混疊，且各分量受干擾噪聲的影響嚴(yán)重，頻譜可讀性較差，無法有效實現(xiàn)齒輪箱復(fù)合故障診斷。

(a)LMD分解后的信號時域圖(PF1～PF3)

5.2 DDS行星齒輪箱復(fù)雜故障平臺試驗驗證

為了更好地表明本文方法的通用性，采用更為復(fù)雜的SpectraQuest公司設(shè)計的傳動故障診斷綜合試驗臺(DDS)進行試驗驗證，試驗平臺如圖23所示。該平臺為一個三級傳動齒輪箱，其中包括一級行星齒輪傳動結(jié)構(gòu)和兩級平行軸齒輪傳動結(jié)構(gòu)，行星齒輪箱中的各級齒輪參數(shù)如表1、表2所示。

圖23 DDS行星齒輪箱試驗臺

表1 行星齒輪的齒輪參數(shù)

表2 平行軸齒輪的齒輪參數(shù)

試驗過程中將第一級行星齒輪中的太陽輪設(shè)置為斷齒故障(局部故障)，第三級平行軸齒輪中主動輪設(shè)置為磨損故障(均勻磨損)，其他齒輪均為正常齒輪。實測驅(qū)動電機的轉(zhuǎn)速為3062 r/min，振動傳感器的采樣頻率為16 kHz，負(fù)載電流設(shè)定為1.5 A。根據(jù)行星齒輪和平行軸齒輪的結(jié)構(gòu)和轉(zhuǎn)動速度，利用公式可計算出故障齒輪的故障特征頻率值。行星齒輪中太陽輪局部故障特征頻率計算公式如下：

(24)

式中，f1為太陽輪的局部故障特征頻率；nd為驅(qū)動電機的轉(zhuǎn)速；Zs為太陽輪齒數(shù)；Zr為內(nèi)齒圈齒數(shù)；Ng為行星輪的個數(shù)。

對于平行軸齒輪來說，當(dāng)有均布故障發(fā)生時，其故障特征頻率為故障齒輪的嚙合頻率，即

(25)

式中，f2為平行軸齒輪的均布磨損故障特征頻率；f0為行星架的轉(zhuǎn)動頻率；Z1、Z2分別為二級平行軸齒輪的主動輪、從動輪齒數(shù)；Z3為三級平行軸主動輪齒數(shù)。

圖24 行星齒輪箱復(fù)合故障信號的循環(huán)譜分析結(jié)果

(a)循環(huán)頻率集為D1時TCS值曲線

(b)循環(huán)頻率集為D2時TCS值曲線圖25 太陽輪斷齒和平行軸齒輪磨損復(fù)合信號的TCS值曲線

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖26 循環(huán)頻率集為D1時ACYCBD盲解卷積結(jié)果

因此，將齒輪相關(guān)參數(shù)代入式(24)、式(25)可計算出太陽輪斷齒故障的特征頻率為159.48 Hz，三級平行軸齒輪中主動輪的磨損故障特征頻率為116.55 Hz。利用本文方法進行故障診斷，得到復(fù)合故障振動信號的循環(huán)譜分析結(jié)果如圖24所示。從圖24中可以發(fā)現(xiàn)存在兩類不同目標(biāo)故障類型的循環(huán)頻率成分，即以159 Hz及其倍頻組成的循環(huán)頻率集，記作D1；116.5 Hz及其倍頻組成的循環(huán)頻率集，記作D2。進一步確定ACYCBD算法所需要的目標(biāo)特征循環(huán)頻率集D1、D2，并根據(jù)TCS曲線選擇出濾波器長度值分別為106、155，如圖25所示。圖26和圖27所示為ACYCBD盲解卷積結(jié)果，可以看出圖26c所示的1.5維導(dǎo)數(shù)譜中將太陽輪斷齒故障及其若干倍頻提取出來，幅值較為突出完整。同理，由圖27c可知，平行軸齒輪磨損故障特征頻率及其倍頻處峰值突出，且解調(diào)效果同樣優(yōu)于傳統(tǒng)的包絡(luò)解調(diào)和1.5維譜方法。接下來與MED、MCKD解卷積方法進行對比，結(jié)果如圖28所示，得到了與前述案例相似的診斷結(jié)果。

(a)包絡(luò)譜

(b)1.5維譜

(c)1.5維導(dǎo)數(shù)譜圖27 循環(huán)頻率集為D2時ACYCBD盲解卷積結(jié)果

(a)MED處理結(jié)果

(b)MCKD分離結(jié)果(太陽輪斷齒故障)

(c)MCKD分離結(jié)果(平行軸齒輪磨損)圖28 MED、MCKD結(jié)合1.5維導(dǎo)數(shù)譜對行星齒輪箱故障信號的分析結(jié)果

最后與LMD進行比較，結(jié)果如圖29所示。分析可知，在PF1分量提取出兩種故障特征頻率的少量倍頻值，表明LMD在分解過程中產(chǎn)生了模態(tài)混疊，且各分量倍頻處的幅值被噪聲干擾。在PF2、PF3分量提取出與平行軸齒輪磨損故障相關(guān)的頻率成分，但有效的特征信息很少，均無法有效實現(xiàn)齒輪箱復(fù)合故障診斷。

為進一步定量論證本文方法的優(yōu)勢，引入裕度因子(margin factor)定量指標(biāo)對不同方法進行性能評價。周期性沖擊信號是齒輪局部故障的主要特征，裕度因子能夠反映振動信號的沖擊程度與離散程度，且穩(wěn)定性較好，其表達(dá)式如下：

(26)

式(26)中的x(n)為不同方法中對原始振動信號進行處理后的濾波信號，裕度因子值越高說明濾波效果越好。不同方法中的裕度因子值如表3所示。分析可知，ACYCBD方法在處理仿真信號時，得到信號1(55 Hz)與信號2(90 Hz)濾波后信號的裕度因子值分別為24.94和15.18，均高于由MCKD方法、MED方法和LMD方法中的PF1～PF3前3個分量計算出的裕度因子值。同樣，在試驗平臺上也定量驗證了本文方法的優(yōu)越性。

(a)LMD分解后的信號時域圖(PF1～PF3)

(b)分解信號的1.5維導(dǎo)數(shù)譜(PF1～PF3)圖29 LMD結(jié)合1.5維導(dǎo)數(shù)譜對行星齒輪箱故障信號分析結(jié)果

表3 不同方法下的裕度因子值

綜上可以發(fā)現(xiàn)，MED方法在進行齒輪箱復(fù)合故障診斷時并不完全適用，這主要是因為MED將故障沖擊信號首個采樣點之前的幅值定義為零，導(dǎo)致其解卷積信號中產(chǎn)生了偽脈沖，且MED方法是基于迭代算法得出濾波器的局部最優(yōu)解，并不能保證是全局最優(yōu)，使得MED方法失效。而MCKD方法主要是受到較多參數(shù)量的影響，同時它也是一種以迭代過程求解最優(yōu)濾波器的方法,盡管能夠?qū)X輪箱復(fù)合故障進行解耦分離，但在噪聲影響下，也只是提取出有限個局部故障沖擊特征。LMD方法在處理齒輪箱多分量耦合調(diào)制復(fù)合故障信號時容易產(chǎn)生模態(tài)混疊，即某個PF分量中同時出現(xiàn)原始信號中不同的尺度特征，導(dǎo)致故障特征提取困難,同時，LMD的局部均值函數(shù)和包絡(luò)估計函數(shù)經(jīng)過多次滑動處理后會產(chǎn)生包絡(luò)誤差，造成部分分解信息失真。而ACYCBD方法約束參數(shù)少，與峭度、相關(guān)峭度等時域指標(biāo)相比，頻域指標(biāo)二階循環(huán)平穩(wěn)性(ICS2)更能夠刻畫故障信號的周期性沖擊特征。綜合以上分析可知，本文方法在降噪、分離和提取齒輪箱復(fù)合故障診斷中具有優(yōu)越性，不過本文方法也存在算法的整體框架相對復(fù)雜，整個算法的運算時間較長，對運算平臺的要求更高等不足。

6 結(jié)論

針對齒輪箱復(fù)合故障特征難以有效提取的問題，提出了一種ACYCBD結(jié)合1.5維導(dǎo)數(shù)譜的復(fù)合故障特征提取方法。該方法能夠有效提取齒輪箱因故障產(chǎn)生的周期性脈沖特征，通過利用循環(huán)譜分析解決了循環(huán)頻率設(shè)定需要先驗知識的問題。并引入三階累積量稀疏度(TCS)指標(biāo)對CYCBD中的濾波器長度進行自適應(yīng)選取，可避免參數(shù)選擇的盲目性，保證了算法處理結(jié)果的有效性。為進一步增強故障沖擊特征，對分離的各個信號分別采用一種1.5維導(dǎo)數(shù)譜來解調(diào)信號中的故障特征頻率，其效果優(yōu)于傳統(tǒng)包絡(luò)解調(diào)與1.5維譜方法。與經(jīng)典的MED、MCKD、LMD方法相比，本文算法有更強的噪聲抑制能力，能夠有效分離和提取齒輪箱復(fù)合故障特征，具有較高的實用價值。