趙藝珂 王家序,2 張 新,2 吳 磊 劉治汶
1.西南交通大學(xué)機(jī)械工程學(xué)院,成都,610031 2.重慶大學(xué)機(jī)械傳動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,重慶,400044 3.電子科技大學(xué)自動(dòng)化工程學(xué)院,成都,611731
航空發(fā)動(dòng)機(jī)作為飛機(jī)的關(guān)鍵部件,結(jié)構(gòu)復(fù)雜且工作環(huán)境惡劣(高溫、高壓、高速、高強(qiáng)度、變負(fù)荷等),其核心傳動(dòng)件(如軸承)易出現(xiàn)故障失效[1-3],因此,研究航空發(fā)動(dòng)機(jī)核心傳動(dòng)件的狀態(tài)監(jiān)測(cè)與故障診斷,對(duì)及時(shí)發(fā)現(xiàn)故障隱患、保證飛行安全、降低運(yùn)維成本具有重要意義。
振動(dòng)信號(hào)分析是軸承故障檢測(cè)的一種有效方法。振動(dòng)信號(hào)包含豐富的軸承故障信息,但受復(fù)雜裝備中多個(gè)振動(dòng)源、復(fù)雜傳遞路徑、強(qiáng)噪聲干擾的影響,故障特征信號(hào)十分微弱,因此如何從振動(dòng)信號(hào)中獲取故障信息是故障診斷的關(guān)鍵[4]。為有效提取故障特征,出現(xiàn)了諸多方法,如小波變換方法[5]、經(jīng)驗(yàn)?zāi)B(tài)分解[6]、局部均值分解[7]、變分模態(tài)分解[8]等。小波變換適用于分析非平穩(wěn)信號(hào)或奇異性突變信號(hào),但小波基的選擇依賴(lài)于經(jīng)驗(yàn),自適應(yīng)性差[9]。經(jīng)驗(yàn)?zāi)B(tài)分解、變分模態(tài)分解等解調(diào)方法能自適應(yīng)分離出周期性故障沖擊,但易出現(xiàn)模態(tài)混疊、端點(diǎn)效應(yīng)、過(guò)(欠)包絡(luò)或分解個(gè)數(shù)不易確定等問(wèn)題[10-11]。
盲解卷積方法給故障診斷提供了新的思路。該方法可在未知振動(dòng)源、傳輸通道特性、噪聲強(qiáng)度的情況下,自適應(yīng)地設(shè)計(jì)出脈沖響應(yīng)濾波器來(lái)恢復(fù)故障沖擊。最小熵解卷積(minimum entropy deconvolution,MED)[12]通過(guò)不斷迭代使濾波信號(hào)的峭度達(dá)到最大,得到盲解卷積濾波器,具有算法簡(jiǎn)單、調(diào)節(jié)參數(shù)少等優(yōu)點(diǎn)。旋轉(zhuǎn)機(jī)械開(kāi)始信號(hào)與零輸入信號(hào)間可能存在不連續(xù)的現(xiàn)象,導(dǎo)致在這個(gè)位置的偽沖擊由于延遲而被解卷積,針對(duì)這一問(wèn)題,ENDO等[13]對(duì)MED進(jìn)行調(diào)整,提出了最小熵解卷積調(diào)整(minimum entropy deconvolution adjustment,MEDA)方法。在此基礎(chǔ)上,CABRELLI[14]提出了通過(guò)求導(dǎo)、矩陣運(yùn)算直接求解濾波器系數(shù)的最優(yōu)最小熵解卷積調(diào)整(optimal minimum entropy deconvolution adjustment,OMEDA)方法,獲得了比MED更加精確的解。然而,MED、MEDA、OMEDA在用于故障診斷時(shí)均傾向于恢復(fù)單個(gè)或少量主導(dǎo)沖擊而非軸承局部缺陷引起的周期性的故障沖擊。為此,MCDONALD等[15]提出了將相關(guān)峭度作為最優(yōu)逆濾波器求解指標(biāo)的最大相關(guān)峭度解卷積方法;隨后,MCDONALD等[16]提出了多點(diǎn)DMEDA方法,通過(guò)設(shè)定一個(gè)定義沖擊位置及權(quán)重的目標(biāo)向量來(lái)恢復(fù)振動(dòng)信號(hào)中的故障沖擊;朱丹宸等[17]提出通過(guò)最大濾波信號(hào)的二階循環(huán)平穩(wěn)量來(lái)求解濾波器系數(shù)的最大二階循環(huán)平穩(wěn)盲解卷積方法;LIANG等[18]提出了最大平均峭度解卷積方法,將信號(hào)每個(gè)故障周期內(nèi)的局部峭度代替整段信號(hào)的峭度,實(shí)現(xiàn)了故障沖擊序列的提取。然而,上述幾種方法需要預(yù)知故障特征頻率,甚至需要嚴(yán)格預(yù)知每個(gè)故障沖擊的位置,而航空發(fā)動(dòng)機(jī)等結(jié)構(gòu)復(fù)雜機(jī)械裝備的故障特征頻率信息往往無(wú)法預(yù)先準(zhǔn)確獲知,因此上述方法在實(shí)際工程中的適用性較差。
基于上述分析,本文提出了一種基于無(wú)偏自相關(guān)分析的增強(qiáng)最小熵解卷積方法。該方法將無(wú)偏自相關(guān)變換融入MED的濾波器系數(shù)的迭代求解,利用無(wú)偏自相關(guān)變換能有效抑制濾波信號(hào)中與軸承故障特征無(wú)關(guān)的非周期性干擾成分的特性,解決了傳統(tǒng)MED傾向于恢復(fù)少量主導(dǎo)沖擊的問(wèn)題,從而實(shí)現(xiàn)軸承微弱故障沖擊序列的增強(qiáng)檢測(cè)。
航空發(fā)動(dòng)機(jī)工作環(huán)境惡劣、結(jié)構(gòu)復(fù)雜,軸承測(cè)量信號(hào)通常含有高斯背景噪聲、諧波分量、沖擊噪聲等與軸承故障特征無(wú)關(guān)的干擾成分,周期性故障沖擊常被干擾所淹沒(méi)(見(jiàn)圖1)。測(cè)量信號(hào)模型可表示為[19]
圖1 周期性故障沖擊信號(hào)的盲解卷積恢復(fù)過(guò)程Fig.1 Blind deconvolution recovery process of the periodic fault impact signal
x=e*ge+n*gn
(1)
式中,*代表卷積運(yùn)算;x為測(cè)量信號(hào);e為周期性故障沖擊;n為高斯噪聲、諧波分量、強(qiáng)沖擊噪聲等干擾成分;ge、gn分別為e、n對(duì)應(yīng)的振動(dòng)傳遞函數(shù)。
測(cè)量信號(hào)的各成分特性復(fù)雜且未知,導(dǎo)致周期性故障沖擊難以被準(zhǔn)確提取。盲解卷積方法為之提供了一種可能的解決途徑,該類(lèi)方法通過(guò)設(shè)計(jì)濾波器h來(lái)對(duì)測(cè)量信號(hào)進(jìn)行解卷積,從而提取周期性故障沖擊e(圖1),即
y=x*h≈e
(2)
式中,y為濾波信號(hào)。
盲解卷積方法利用周期性故障沖擊與干擾成分在沖擊性、循環(huán)平穩(wěn)性或稀疏性等特性方面可能存在不同的表現(xiàn)來(lái)設(shè)計(jì)濾波器,達(dá)到恢復(fù)故障沖擊特征的目的,因此,方法的難點(diǎn)在于如何合理找尋和利用上述特性。
MED作為經(jīng)典的盲解卷積方法,在諸多領(lǐng)域有著廣泛應(yīng)用。該方法利用目標(biāo)信號(hào)與干擾成分在峭度上的差異,使濾波信號(hào)的峭度達(dá)到最大來(lái)實(shí)現(xiàn)濾波器系數(shù)的求解,即
(3)
(4)
峭度作為一種非高斯度量指標(biāo),能在一定程度上刻畫(huà)軸承局部缺陷引起的故障沖擊,因此基于峭度的MED方法可用于旋轉(zhuǎn)機(jī)械的故障診斷。MED傾向于得到單個(gè)或少量主導(dǎo)沖擊而非由軸承故障產(chǎn)生的周期性故障沖擊序列。這是因?yàn)閱蝹€(gè)沖擊信號(hào)的峭度通常遠(yuǎn)大于同信號(hào)長(zhǎng)度的周期性沖擊序列的峭度。如圖2所示,周期性故障沖擊序列
(a)周期性故障沖擊的峭度K=10.58
(b)單個(gè)強(qiáng)沖擊的峭度K=149.7圖2 周期性故障沖擊與單個(gè)強(qiáng)沖擊信號(hào)的峭度Fig.2 Kurtosis of periodic fault impacts and single strong impact
t≥0.01k
(5)
(信號(hào)采樣頻率為10 kHz)的峭度為10.58,而相同信號(hào)長(zhǎng)度的單個(gè)沖擊信號(hào)
(6)
(信號(hào)采樣頻率為10 kHz)的峭度高達(dá)149.7。實(shí)際測(cè)量信號(hào)可能還含有少量隨機(jī)強(qiáng)沖擊干擾(如突發(fā)性外部沖擊干擾),導(dǎo)致MED傾向于得到單個(gè)或少量主導(dǎo)沖擊的問(wèn)題變得更為明顯。
本文增強(qiáng)最小熵解卷積方法將無(wú)偏自相關(guān)變換融入MED的濾波器系數(shù)的迭代求解中,在計(jì)算濾波信號(hào)峭度前,先對(duì)濾波信號(hào)進(jìn)行無(wú)偏自相關(guān)變換,以抑制信號(hào)中噪聲、強(qiáng)沖擊等非周期成分的干擾,便于準(zhǔn)確提取周期性的故障沖擊。信號(hào)y的無(wú)偏自相關(guān)變換定義如下:
(7)
式中,τ=q/fs為延遲系數(shù),q=0,1,…,N-1;N為信號(hào)長(zhǎng)度;fs為信號(hào)采樣頻率。
由上述定義可見(jiàn),無(wú)偏自相關(guān)變換可用于檢測(cè)軸承故障信號(hào)中的周期性特征成分,與軸承故障不相關(guān)的非周期成分信號(hào)(噪聲、隨機(jī)強(qiáng)沖擊干擾等)幅值在無(wú)偏自相關(guān)變換后會(huì)驟減接近于零。如圖3~圖5所示,無(wú)偏自相關(guān)變換后,高斯噪聲與強(qiáng)沖擊信號(hào)的幅值接近于零,遠(yuǎn)小于軸承周期性故障沖擊信號(hào)的幅值;周期性故障信號(hào)變換前后具有相同的周期性。由此可見(jiàn),將無(wú)偏自相關(guān)變換融入MED的濾波器系數(shù)求解中,能有效解決MED易于得到單個(gè)或少量強(qiáng)沖擊而非周期性的軸承故障沖擊的問(wèn)題。
(a)高斯噪聲
(b)高斯噪聲的無(wú)偏自相關(guān)變換圖3 高斯噪聲及其無(wú)偏自相關(guān)變換Fig.3 Gussian noise and its unbiased autocorrelation transformation
(a)強(qiáng)沖擊
(b)強(qiáng)沖擊的無(wú)偏自相關(guān)變換圖4 強(qiáng)沖擊及其無(wú)偏自相關(guān)變換Fig.4 Strong impact and its unbiased autocorrelation transformation
(a)周期性的故障沖擊
(b)故障沖擊的無(wú)偏自相關(guān)變換圖5 周期性的故障沖擊信號(hào)及其無(wú)偏自相關(guān)變換Fig.5 Periodic fault impact and its unbiased autocorrelation transformation
無(wú)偏自相關(guān)變換的輸出信號(hào)點(diǎn)數(shù)為輸入信號(hào)點(diǎn)數(shù)的兩倍,在時(shí)域上是前后對(duì)稱(chēng)的。由于信號(hào)中的非周期成分(高斯噪聲和強(qiáng)沖擊干擾)在自相關(guān)變換后的輸出信號(hào)始端出現(xiàn)峰值,而信號(hào)末端由于樣本減少會(huì)出現(xiàn)幅值異常,因此有必要去掉自相關(guān)變換輸出信號(hào)兩端少量的數(shù)據(jù)點(diǎn)[20]?;诖?,在所提方法中,筆者選擇自相關(guān)變換輸出信號(hào)的后半部分,并去掉所選擇信號(hào)的前1%個(gè)點(diǎn)和后10%個(gè)點(diǎn)。
基于上述無(wú)偏自相關(guān)變換,增強(qiáng)最小熵解卷積方法在求解濾波器系數(shù)時(shí)可等效于求解下述優(yōu)化問(wèn)題:
(8)
M1=1.01N+1M2=1.9N-1
類(lèi)似于MED,式(8)所示的濾波器系數(shù)的求解可采用直接求解法或迭代求解法[21]。直接求解法通過(guò)令目標(biāo)函數(shù)對(duì)濾波器系數(shù)的偏導(dǎo)數(shù)為零來(lái)求解濾波器系數(shù),由于需要進(jìn)行矩陣運(yùn)算,因此該方法的計(jì)算效率較低。迭代求解法通過(guò)不斷令目標(biāo)函數(shù)最大來(lái)迭代求解濾波器系數(shù),具有算法簡(jiǎn)單、運(yùn)算高效等優(yōu)勢(shì),得到了更廣泛的應(yīng)用,MED、最大二階循環(huán)平穩(wěn)盲解卷積和最大相關(guān)峭度解卷積均采用此法。鑒于此,本文也采用迭代法求解濾波器系數(shù),具體求解步驟如算法1所示。算法1中,tZ為給定的迭代次數(shù),L為濾波器系數(shù)的長(zhǎng)度,第7步的濾波器系數(shù)的更新方式如下:
(9)
(10)
X0=
(11)
算法1:基于無(wú)偏自相關(guān)分析的增強(qiáng)最小熵解卷積方法輸入:測(cè)量信號(hào)x1.初始化:L,初始濾波器h(1)=[01…0]T,tZ=502.t=13.迭代:4. 按式(2)求解濾波信號(hào)y(t)=x*h(t)5. 按式(7)對(duì)濾波信號(hào)進(jìn)行無(wú)偏自相關(guān)變換R^yy6. 按式(4)計(jì)算R^yy(M1,M2)的峭度7. 按式(9)~式(11)更新濾波器系數(shù)h(t)8. t←t+19. 終止迭代:當(dāng)t>tZ10.h^←hopt(hopt為R^yy(M1,M2)峭度最大時(shí)的濾波器)輸出:故障沖擊估計(jì)值:e^=R^yy(M1,M2)(y=x*h^)
為驗(yàn)證所提方法對(duì)提取周期性故障沖擊的有效性,進(jìn)行了仿真分析。仿真信號(hào)x(t)由故障沖擊e(t)、強(qiáng)沖擊干擾分量b(t)、諧波分量c(t)和高斯噪聲n(t)組成,即x(t)=e(t)+b(t)+c(t)+n(t),其中,c(t)=0.2sin(πt/4000),n(t)為噪聲(加入n(t)后的信號(hào)信噪比為-9 dB)。仿真信號(hào)的采樣頻率設(shè)定為10 kHz。由圖6可見(jiàn),故障沖擊序列為有規(guī)律的周期性信號(hào),但混合信號(hào)中的高斯噪聲、諧波分量和強(qiáng)沖擊干擾成分完全破壞了故障沖擊特征。
(a)周期性的故障沖擊e(t)
(b)單個(gè)強(qiáng)沖擊b(t)
(c)諧波分量c(t)
(d)高斯噪聲n(t)
(e)混合信號(hào)x(t)圖6 仿真信號(hào)Fig.6 The simulated signals
由圖7可見(jiàn),所提方法的濾波信號(hào)在時(shí)域上呈現(xiàn)出明顯的周期性沖擊特征(圖7a),從包絡(luò)譜(圖7b)中能準(zhǔn)確識(shí)別出故障特征頻率f=100 Hz及其倍頻。MED(圖8a)、MEDA(圖9a)與OMEDA(圖10a)的濾波信號(hào)在時(shí)域上均表現(xiàn)出單個(gè)主導(dǎo)沖擊,同時(shí)3種方法所得到的包絡(luò)譜(圖8b、圖9b、圖10b)中故障特征頻率及其倍頻成分不凸顯,干擾頻率較多,故障特征辨識(shí)度較差。
(a)時(shí)域波形
(b)包絡(luò)譜圖7 所提方法對(duì)仿真信號(hào)的分析結(jié)果Fig.7 Analysis result of the proposed method for simulated signal
(a)時(shí)域波形
(b)包絡(luò)譜圖8 MED方法對(duì)仿真信號(hào)的分析結(jié)果Fig.8 Analysis result of the MED for simulated signal
(a)時(shí)域波形
(b)包絡(luò)譜圖9 MEDA方法對(duì)仿真信號(hào)的分析結(jié)果Fig.9 Analysis result of the MEDA for simulated signal
(a)時(shí)域波形
(b)包絡(luò)譜圖10 OMEDA方法對(duì)仿真信號(hào)的分析結(jié)果Fig.10 Analysis result of the OMEDA for simulated signal
為進(jìn)一步凸顯所提方法的優(yōu)勢(shì),本文選用基于峭度最大原則的FK(fast kurtogram)方法進(jìn)行對(duì)比分析。FK方法作為周期性沖擊檢測(cè)方法,在故障診斷領(lǐng)域有著非常廣泛的應(yīng)用,是一種衡量其他方法有效性的基準(zhǔn)方法。如圖11所示,F(xiàn)K方法的濾波信號(hào)只提取了仿真信號(hào)中的單個(gè)強(qiáng)沖擊干擾成分,濾波信號(hào)的平方包絡(luò)譜中也無(wú)法觀測(cè)到與故障有關(guān)的特征信息。
(a)譜峭度圖
(b)濾波信號(hào)
(c)平方包絡(luò)譜圖11 FK方法對(duì)仿真信號(hào)的分析結(jié)果Fig.11 Analysis result of the FK for simulated signal
從上述分析可見(jiàn),MED、MEDA、OMEDA、FK在分析成分復(fù)雜的振動(dòng)信號(hào)時(shí),易得到單個(gè)主導(dǎo)沖擊特征而不是軸承故障產(chǎn)生的周期性故障沖擊。所提方法將無(wú)偏自相關(guān)變換融入到MED的濾波器系數(shù)的迭代求解中,能有效抑制濾波信號(hào)中的高斯噪聲、強(qiáng)沖擊干擾等非周期性成分,準(zhǔn)確提取軸承故障周期性沖擊特征,優(yōu)勢(shì)較明顯。
為驗(yàn)證所提方法對(duì)工程實(shí)際中軸承故障診斷的有效性,進(jìn)行了航空發(fā)動(dòng)機(jī)故障診斷的應(yīng)用。工程數(shù)據(jù)為某航空發(fā)動(dòng)機(jī)附件齒輪箱的軸承故障振動(dòng)數(shù)據(jù)[22]。發(fā)動(dòng)機(jī)結(jié)構(gòu)及傳感器安裝位置如圖12所示。軸承故障情況如下:齒輪箱L1軸軸承外圈有嚴(yán)重劃痕損傷,L5軸軸承外圈有局部剝落。
圖12 發(fā)動(dòng)機(jī)示意圖及傳感器位置Fig.12 Schematic diagram of the aero-engine and the location of sensors
L4軸安裝了脈沖轉(zhuǎn)速計(jì)(每轉(zhuǎn)44脈沖),L1軸、L5軸附近均安裝了加速度計(jì)。信號(hào)采樣頻率為44 100 Hz,L4軸上轉(zhuǎn)速計(jì)的實(shí)測(cè)轉(zhuǎn)頻為180 Hz。數(shù)據(jù)采集過(guò)程中,L1軸附近的加速度計(jì)出現(xiàn)故障,僅L5軸附近加速度計(jì)的測(cè)量信號(hào)可用。根據(jù)L4軸、L5軸的齒輪齒數(shù)62和61,可得L5軸的轉(zhuǎn)頻fr=180×62÷61=182.95 Hz,然后根據(jù)L5軸軸承的故障特征頻率相對(duì)于L5軸轉(zhuǎn)頻的比例系數(shù)(表1)[22],可得L5軸承外圈故障特征頻率fo=7.759fr=1419.51 Hz。
表1 軸承故障特征頻率相對(duì)于L5軸轉(zhuǎn)頻的比例系數(shù)Tab.1 Fault frequencies of bearings of shaft L5 referenced to shaft rotations of L5
圖13給出了L5軸附近加速度計(jì)16 384個(gè)測(cè)量信號(hào)的時(shí)域波形及其包絡(luò)譜,可見(jiàn),時(shí)域波形中的周期性故障沖擊信號(hào)被完全淹沒(méi),包絡(luò)譜中只能觀察到轉(zhuǎn)頻fr,無(wú)法觀察到軸承外圈故障特征頻率。
(a)時(shí)域波形
(b)包絡(luò)譜圖13 航空發(fā)動(dòng)機(jī)故障軸承測(cè)量信號(hào)Fig.13 The measured signal of the aero-engine bearing
圖14~圖17所示為4種盲解卷積方法的濾波信號(hào)及其包絡(luò)譜。對(duì)比濾波信號(hào)可見(jiàn),所提方法外的其余方法均只得到單個(gè)主導(dǎo)沖擊。從所提方法濾波信號(hào)的包絡(luò)譜中能準(zhǔn)確發(fā)現(xiàn)L5軸軸承外圈的故障特征頻率fo。MED濾波信號(hào)包絡(luò)譜明顯不含軸承故障的相關(guān)頻率特征信息,從MEDA與OMEDA的濾波信號(hào)包絡(luò)譜中雖能發(fā)現(xiàn)微弱的故障特征頻率fo,但復(fù)雜的干擾頻率成分影響了故障診斷的準(zhǔn)確性。
(a)時(shí)域波形
(b)包絡(luò)譜圖14 所提方法對(duì)航空發(fā)動(dòng)機(jī)軸承故障信號(hào)的分析結(jié)果Fig.14 Analysis result of the proposed method for measured signal of the aero-engine bearing
(a)時(shí)域波形
(b)包絡(luò)譜圖15 MED方法對(duì)航空發(fā)動(dòng)機(jī)軸承故障信號(hào)的分析結(jié)果Fig.15 Analysis result of the MED for measured signal of the aero-engine bearing
(a)時(shí)域波形
(b)包絡(luò)譜圖16 MEDA方法對(duì)航空發(fā)動(dòng)機(jī)軸承故障信號(hào)的分析結(jié)果Fig.16 Analysis result of the MEDA for measured signal of the aero-engine bearing
(a)時(shí)域波形
(b)包絡(luò)譜圖17 OMEDA方法對(duì)航空發(fā)動(dòng)機(jī)軸承故障信號(hào)的分析結(jié)果Fig.17 Analysis result of the OMEDA for measured signal of the aero-engine bearing
如圖18所示,F(xiàn)K的分析結(jié)果中,濾波信號(hào)具有沖擊序列的特征,從包絡(luò)譜能觀察到故障特征頻率fo,但fo附近存在較強(qiáng)的干擾頻率成分,特征識(shí)別效果較差。
(a)譜峭度圖
(b)濾波信號(hào)
(c)平方包絡(luò)譜圖18 FK方法對(duì)航空發(fā)動(dòng)機(jī)軸承故障信號(hào)的分析結(jié)果Fig.18 Analysis result of the FK for measured signal of the aero-engine bearing
由上述分析結(jié)果可見(jiàn),對(duì)于航空發(fā)動(dòng)機(jī)而言,其復(fù)雜結(jié)構(gòu)、惡劣工況等因素導(dǎo)致軸承故障特征極其微弱,傳統(tǒng)盲解卷積方法及經(jīng)典的FK方法難以滿足信號(hào)分析的需求,而所提方法表現(xiàn)出較強(qiáng)的適用性,在軸承故障診斷中優(yōu)勢(shì)明顯。
(1)提出了一種基于無(wú)偏自相關(guān)分析的增強(qiáng)最小熵解卷積方法,有效解決了傳統(tǒng)盲解卷積方法傾向于恢復(fù)少量主導(dǎo)沖擊而非軸承周期性故障沖擊的問(wèn)題。
(2)仿真信號(hào)分析結(jié)果驗(yàn)證了所提方法在高斯噪聲、強(qiáng)沖擊等復(fù)雜干擾下提取軸承周期性故障沖擊的有效性。
(3)航空發(fā)動(dòng)機(jī)軸承故障診斷案例分析進(jìn)一步證實(shí)了所提方法較傳統(tǒng)盲解卷積方法以及經(jīng)典FK方法具有更強(qiáng)的適用性。