趙藝珂 王家序，2 張 新，2 吳 磊 劉治汶

1.西南交通大學機械工程學院，成都，610031 2.重慶大學機械傳動國家重點實驗室，重慶，400044 3.電子科技大學自動化工程學院，成都，611731

0 引言

航空發(fā)動機作為飛機的關鍵部件，結構復雜且工作環(huán)境惡劣(高溫、高壓、高速、高強度、變負荷等)，其核心傳動件(如軸承)易出現(xiàn)故障失效[1-3]，因此，研究航空發(fā)動機核心傳動件的狀態(tài)監(jiān)測與故障診斷，對及時發(fā)現(xiàn)故障隱患、保證飛行安全、降低運維成本具有重要意義。

振動信號分析是軸承故障檢測的一種有效方法。振動信號包含豐富的軸承故障信息，但受復雜裝備中多個振動源、復雜傳遞路徑、強噪聲干擾的影響，故障特征信號十分微弱，因此如何從振動信號中獲取故障信息是故障診斷的關鍵[4]。為有效提取故障特征，出現(xiàn)了諸多方法，如小波變換方法[5]、經(jīng)驗模態(tài)分解[6]、局部均值分解[7]、變分模態(tài)分解[8]等。小波變換適用于分析非平穩(wěn)信號或奇異性突變信號，但小波基的選擇依賴于經(jīng)驗，自適應性差[9]。經(jīng)驗模態(tài)分解、變分模態(tài)分解等解調(diào)方法能自適應分離出周期性故障沖擊，但易出現(xiàn)模態(tài)混疊、端點效應、過(欠)包絡或分解個數(shù)不易確定等問題[10-11]。

盲解卷積方法給故障診斷提供了新的思路。該方法可在未知振動源、傳輸通道特性、噪聲強度的情況下，自適應地設計出脈沖響應濾波器來恢復故障沖擊。最小熵解卷積(minimum entropy deconvolution，MED)[12]通過不斷迭代使濾波信號的峭度達到最大，得到盲解卷積濾波器，具有算法簡單、調(diào)節(jié)參數(shù)少等優(yōu)點。旋轉機械開始信號與零輸入信號間可能存在不連續(xù)的現(xiàn)象，導致在這個位置的偽沖擊由于延遲而被解卷積，針對這一問題，ENDO等[13]對MED進行調(diào)整，提出了最小熵解卷積調(diào)整(minimum entropy deconvolution adjustment，MEDA)方法。在此基礎上，CABRELLI[14]提出了通過求導、矩陣運算直接求解濾波器系數(shù)的最優(yōu)最小熵解卷積調(diào)整(optimal minimum entropy deconvolution adjustment，OMEDA)方法，獲得了比MED更加精確的解。然而，MED、MEDA、OMEDA在用于故障診斷時均傾向于恢復單個或少量主導沖擊而非軸承局部缺陷引起的周期性的故障沖擊。為此，MCDONALD等[15]提出了將相關峭度作為最優(yōu)逆濾波器求解指標的最大相關峭度解卷積方法；隨后，MCDONALD等[16]提出了多點DMEDA方法，通過設定一個定義沖擊位置及權重的目標向量來恢復振動信號中的故障沖擊；朱丹宸等[17]提出通過最大濾波信號的二階循環(huán)平穩(wěn)量來求解濾波器系數(shù)的最大二階循環(huán)平穩(wěn)盲解卷積方法；LIANG等[18]提出了最大平均峭度解卷積方法，將信號每個故障周期內(nèi)的局部峭度代替整段信號的峭度，實現(xiàn)了故障沖擊序列的提取。然而，上述幾種方法需要預知故障特征頻率，甚至需要嚴格預知每個故障沖擊的位置，而航空發(fā)動機等結構復雜機械裝備的故障特征頻率信息往往無法預先準確獲知，因此上述方法在實際工程中的適用性較差。

基于上述分析，本文提出了一種基于無偏自相關分析的增強最小熵解卷積方法。該方法將無偏自相關變換融入MED的濾波器系數(shù)的迭代求解，利用無偏自相關變換能有效抑制濾波信號中與軸承故障特征無關的非周期性干擾成分的特性，解決了傳統(tǒng)MED傾向于恢復少量主導沖擊的問題，從而實現(xiàn)軸承微弱故障沖擊序列的增強檢測。

1 問題描述

1.1 故障信號盲解卷積模型

航空發(fā)動機工作環(huán)境惡劣、結構復雜，軸承測量信號通常含有高斯背景噪聲、諧波分量、沖擊噪聲等與軸承故障特征無關的干擾成分，周期性故障沖擊常被干擾所淹沒(見圖1)。測量信號模型可表示為[19]

圖1 周期性故障沖擊信號的盲解卷積恢復過程Fig.1 Blind deconvolution recovery process of the periodic fault impact signal

x=e*ge+n*gn

(1)

式中，*代表卷積運算；x為測量信號；e為周期性故障沖擊；n為高斯噪聲、諧波分量、強沖擊噪聲等干擾成分；ge、gn分別為e、n對應的振動傳遞函數(shù)。

測量信號的各成分特性復雜且未知，導致周期性故障沖擊難以被準確提取。盲解卷積方法為之提供了一種可能的解決途徑，該類方法通過設計濾波器h來對測量信號進行解卷積，從而提取周期性故障沖擊e(圖1)，即

y=x*h≈e

(2)

式中，y為濾波信號。

盲解卷積方法利用周期性故障沖擊與干擾成分在沖擊性、循環(huán)平穩(wěn)性或稀疏性等特性方面可能存在不同的表現(xiàn)來設計濾波器，達到恢復故障沖擊特征的目的，因此，方法的難點在于如何合理找尋和利用上述特性。

1.2 MED方法及其主要問題

MED作為經(jīng)典的盲解卷積方法，在諸多領域有著廣泛應用。該方法利用目標信號與干擾成分在峭度上的差異，使濾波信號的峭度達到最大來實現(xiàn)濾波器系數(shù)的求解，即

(3)

(4)

峭度作為一種非高斯度量指標，能在一定程度上刻畫軸承局部缺陷引起的故障沖擊，因此基于峭度的MED方法可用于旋轉機械的故障診斷。MED傾向于得到單個或少量主導沖擊而非由軸承故障產(chǎn)生的周期性故障沖擊序列。這是因為單個沖擊信號的峭度通常遠大于同信號長度的周期性沖擊序列的峭度。如圖2所示，周期性故障沖擊序列

(a)周期性故障沖擊的峭度K=10.58

(b)單個強沖擊的峭度K=149.7圖2 周期性故障沖擊與單個強沖擊信號的峭度Fig.2 Kurtosis of periodic fault impacts and single strong impact

t≥0.01k

(5)

(信號采樣頻率為10 kHz)的峭度為10.58，而相同信號長度的單個沖擊信號

(6)

(信號采樣頻率為10 kHz)的峭度高達149.7。實際測量信號可能還含有少量隨機強沖擊干擾(如突發(fā)性外部沖擊干擾)，導致MED傾向于得到單個或少量主導沖擊的問題變得更為明顯。

2 增強最小熵解卷積方法

本文增強最小熵解卷積方法將無偏自相關變換融入MED的濾波器系數(shù)的迭代求解中，在計算濾波信號峭度前，先對濾波信號進行無偏自相關變換，以抑制信號中噪聲、強沖擊等非周期成分的干擾，便于準確提取周期性的故障沖擊。信號y的無偏自相關變換定義如下：

(7)

式中，τ=q/fs為延遲系數(shù)，q=0，1，…，N-1；N為信號長度；fs為信號采樣頻率。

由上述定義可見，無偏自相關變換可用于檢測軸承故障信號中的周期性特征成分，與軸承故障不相關的非周期成分信號(噪聲、隨機強沖擊干擾等)幅值在無偏自相關變換后會驟減接近于零。如圖3～圖5所示，無偏自相關變換后，高斯噪聲與強沖擊信號的幅值接近于零，遠小于軸承周期性故障沖擊信號的幅值；周期性故障信號變換前后具有相同的周期性。由此可見，將無偏自相關變換融入MED的濾波器系數(shù)求解中，能有效解決MED易于得到單個或少量強沖擊而非周期性的軸承故障沖擊的問題。

(a)高斯噪聲

(b)高斯噪聲的無偏自相關變換圖3 高斯噪聲及其無偏自相關變換Fig.3 Gussian noise and its unbiased autocorrelation transformation

(a)強沖擊

(b)強沖擊的無偏自相關變換圖4 強沖擊及其無偏自相關變換Fig.4 Strong impact and its unbiased autocorrelation transformation

(a)周期性的故障沖擊

(b)故障沖擊的無偏自相關變換圖5 周期性的故障沖擊信號及其無偏自相關變換Fig.5 Periodic fault impact and its unbiased autocorrelation transformation

無偏自相關變換的輸出信號點數(shù)為輸入信號點數(shù)的兩倍，在時域上是前后對稱的。由于信號中的非周期成分(高斯噪聲和強沖擊干擾)在自相關變換后的輸出信號始端出現(xiàn)峰值，而信號末端由于樣本減少會出現(xiàn)幅值異常，因此有必要去掉自相關變換輸出信號兩端少量的數(shù)據(jù)點[20]?；诖耍谒岱椒ㄖ?，筆者選擇自相關變換輸出信號的后半部分，并去掉所選擇信號的前1%個點和后10%個點。

基于上述無偏自相關變換，增強最小熵解卷積方法在求解濾波器系數(shù)時可等效于求解下述優(yōu)化問題：

(8)

M1=1.01N+1M2=1.9N-1

類似于MED，式(8)所示的濾波器系數(shù)的求解可采用直接求解法或迭代求解法[21]。直接求解法通過令目標函數(shù)對濾波器系數(shù)的偏導數(shù)為零來求解濾波器系數(shù)，由于需要進行矩陣運算，因此該方法的計算效率較低。迭代求解法通過不斷令目標函數(shù)最大來迭代求解濾波器系數(shù)，具有算法簡單、運算高效等優(yōu)勢，得到了更廣泛的應用，MED、最大二階循環(huán)平穩(wěn)盲解卷積和最大相關峭度解卷積均采用此法。鑒于此，本文也采用迭代法求解濾波器系數(shù)，具體求解步驟如算法1所示。算法1中，tZ為給定的迭代次數(shù)，L為濾波器系數(shù)的長度，第7步的濾波器系數(shù)的更新方式如下：

(9)

(10)

X0=

(11)

算法1:基于無偏自相關分析的增強最小熵解卷積方法輸入:測量信號x1.初始化:L,初始濾波器h(1)=[01…0]T,tZ=502.t=13.迭代:4. 按式(2)求解濾波信號y(t)=x*h(t)5. 按式(7)對濾波信號進行無偏自相關變換R^yy6. 按式(4)計算R^yy(M1,M2)的峭度7. 按式(9)～式(11)更新濾波器系數(shù)h(t)8. t←t+19. 終止迭代:當t>tZ10.h^←hopt(hopt為R^yy(M1,M2)峭度最大時的濾波器)輸出:故障沖擊估計值:e^=R^yy(M1,M2)(y=x*h^)

3 仿真信號分析

為驗證所提方法對提取周期性故障沖擊的有效性，進行了仿真分析。仿真信號x(t)由故障沖擊e(t)、強沖擊干擾分量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)后的信號信噪比為-9 dB)。仿真信號的采樣頻率設定為10 kHz。由圖6可見，故障沖擊序列為有規(guī)律的周期性信號，但混合信號中的高斯噪聲、諧波分量和強沖擊干擾成分完全破壞了故障沖擊特征。

(a)周期性的故障沖擊e(t)

(b)單個強沖擊b(t)

(c)諧波分量c(t)

(d)高斯噪聲n(t)

(e)混合信號x(t)圖6 仿真信號Fig.6 The simulated signals

由圖7可見，所提方法的濾波信號在時域上呈現(xiàn)出明顯的周期性沖擊特征(圖7a)，從包絡譜(圖7b)中能準確識別出故障特征頻率f=100 Hz及其倍頻。MED(圖8a)、MEDA(圖9a)與OMEDA(圖10a)的濾波信號在時域上均表現(xiàn)出單個主導沖擊，同時3種方法所得到的包絡譜(圖8b、圖9b、圖10b)中故障特征頻率及其倍頻成分不凸顯，干擾頻率較多，故障特征辨識度較差。

(a)時域波形

(b)包絡譜圖7 所提方法對仿真信號的分析結果Fig.7 Analysis result of the proposed method for simulated signal

(a)時域波形

(b)包絡譜圖8 MED方法對仿真信號的分析結果Fig.8 Analysis result of the MED for simulated signal

(a)時域波形

(b)包絡譜圖9 MEDA方法對仿真信號的分析結果Fig.9 Analysis result of the MEDA for simulated signal

(a)時域波形

(b)包絡譜圖10 OMEDA方法對仿真信號的分析結果Fig.10 Analysis result of the OMEDA for simulated signal

為進一步凸顯所提方法的優(yōu)勢，本文選用基于峭度最大原則的FK(fast kurtogram)方法進行對比分析。FK方法作為周期性沖擊檢測方法，在故障診斷領域有著非常廣泛的應用，是一種衡量其他方法有效性的基準方法。如圖11所示，F(xiàn)K方法的濾波信號只提取了仿真信號中的單個強沖擊干擾成分，濾波信號的平方包絡譜中也無法觀測到與故障有關的特征信息。

(a)譜峭度圖

(b)濾波信號

(c)平方包絡譜圖11 FK方法對仿真信號的分析結果Fig.11 Analysis result of the FK for simulated signal

從上述分析可見，MED、MEDA、OMEDA、FK在分析成分復雜的振動信號時，易得到單個主導沖擊特征而不是軸承故障產(chǎn)生的周期性故障沖擊。所提方法將無偏自相關變換融入到MED的濾波器系數(shù)的迭代求解中，能有效抑制濾波信號中的高斯噪聲、強沖擊干擾等非周期性成分，準確提取軸承故障周期性沖擊特征，優(yōu)勢較明顯。

4 航空發(fā)動機軸承故障診斷

為驗證所提方法對工程實際中軸承故障診斷的有效性，進行了航空發(fā)動機故障診斷的應用。工程數(shù)據(jù)為某航空發(fā)動機附件齒輪箱的軸承故障振動數(shù)據(jù)[22]。發(fā)動機結構及傳感器安裝位置如圖12所示。軸承故障情況如下：齒輪箱L1軸軸承外圈有嚴重劃痕損傷，L5軸軸承外圈有局部剝落。

圖12 發(fā)動機示意圖及傳感器位置Fig.12 Schematic diagram of the aero-engine and the location of sensors

L4軸安裝了脈沖轉速計(每轉44脈沖)，L1軸、L5軸附近均安裝了加速度計。信號采樣頻率為44 100 Hz，L4軸上轉速計的實測轉頻為180 Hz。數(shù)據(jù)采集過程中，L1軸附近的加速度計出現(xiàn)故障，僅L5軸附近加速度計的測量信號可用。根據(jù)L4軸、L5軸的齒輪齒數(shù)62和61，可得L5軸的轉頻fr=180×62÷61=182.95 Hz，然后根據(jù)L5軸軸承的故障特征頻率相對于L5軸轉頻的比例系數(shù)(表1)[22]，可得L5軸承外圈故障特征頻率fo=7.759fr=1419.51 Hz。

表1 軸承故障特征頻率相對于L5軸轉頻的比例系數(shù)Tab.1 Fault frequencies of bearings of shaft L5 referenced to shaft rotations of L5

圖13給出了L5軸附近加速度計16 384個測量信號的時域波形及其包絡譜，可見，時域波形中的周期性故障沖擊信號被完全淹沒，包絡譜中只能觀察到轉頻fr，無法觀察到軸承外圈故障特征頻率。

(a)時域波形

(b)包絡譜圖13 航空發(fā)動機故障軸承測量信號Fig.13 The measured signal of the aero-engine bearing

圖14～圖17所示為4種盲解卷積方法的濾波信號及其包絡譜。對比濾波信號可見，所提方法外的其余方法均只得到單個主導沖擊。從所提方法濾波信號的包絡譜中能準確發(fā)現(xiàn)L5軸軸承外圈的故障特征頻率fo。MED濾波信號包絡譜明顯不含軸承故障的相關頻率特征信息，從MEDA與OMEDA的濾波信號包絡譜中雖能發(fā)現(xiàn)微弱的故障特征頻率fo，但復雜的干擾頻率成分影響了故障診斷的準確性。

(a)時域波形

(b)包絡譜圖14 所提方法對航空發(fā)動機軸承故障信號的分析結果Fig.14 Analysis result of the proposed method for measured signal of the aero-engine bearing

(a)時域波形

(b)包絡譜圖15 MED方法對航空發(fā)動機軸承故障信號的分析結果Fig.15 Analysis result of the MED for measured signal of the aero-engine bearing

(a)時域波形

(b)包絡譜圖16 MEDA方法對航空發(fā)動機軸承故障信號的分析結果Fig.16 Analysis result of the MEDA for measured signal of the aero-engine bearing

(a)時域波形

(b)包絡譜圖17 OMEDA方法對航空發(fā)動機軸承故障信號的分析結果Fig.17 Analysis result of the OMEDA for measured signal of the aero-engine bearing

如圖18所示，F(xiàn)K的分析結果中，濾波信號具有沖擊序列的特征，從包絡譜能觀察到故障特征頻率fo，但fo附近存在較強的干擾頻率成分，特征識別效果較差。

(a)譜峭度圖

(b)濾波信號

(c)平方包絡譜圖18 FK方法對航空發(fā)動機軸承故障信號的分析結果Fig.18 Analysis result of the FK for measured signal of the aero-engine bearing

由上述分析結果可見，對于航空發(fā)動機而言，其復雜結構、惡劣工況等因素導致軸承故障特征極其微弱，傳統(tǒng)盲解卷積方法及經(jīng)典的FK方法難以滿足信號分析的需求，而所提方法表現(xiàn)出較強的適用性，在軸承故障診斷中優(yōu)勢明顯。

5 結論

(1)提出了一種基于無偏自相關分析的增強最小熵解卷積方法，有效解決了傳統(tǒng)盲解卷積方法傾向于恢復少量主導沖擊而非軸承周期性故障沖擊的問題。

(2)仿真信號分析結果驗證了所提方法在高斯噪聲、強沖擊等復雜干擾下提取軸承周期性故障沖擊的有效性。

(3)航空發(fā)動機軸承故障診斷案例分析進一步證實了所提方法較傳統(tǒng)盲解卷積方法以及經(jīng)典FK方法具有更強的適用性。