王銀玲 尹顯明 李華聰

(1. 西南科技大學(xué)工程技術(shù)中心 四川綿陽(yáng) 621010； 2. 西北工業(yè)大學(xué)動(dòng)力與能源學(xué)院 西安 710072)

滾動(dòng)軸承廣泛應(yīng)用在旋轉(zhuǎn)機(jī)械中，受負(fù)荷、摩擦、沖擊等非線性因素影響會(huì)出現(xiàn)損傷。通常由于滾動(dòng)軸承工作環(huán)境復(fù)雜，傳感器采集到的信號(hào)含有大量的背景噪聲，造成損傷信號(hào)的非線性、非平穩(wěn)、沖擊性特征被淹沒(méi)。因此，如何通過(guò)對(duì)損傷過(guò)程中產(chǎn)生的非線性、非平穩(wěn)隨機(jī)信號(hào)進(jìn)行分析，準(zhǔn)確檢測(cè)出設(shè)備的故障信號(hào)是滾動(dòng)軸承故障檢測(cè)研究的熱點(diǎn)及難點(diǎn)[1-3]。

在軸承故障診斷中，損傷信號(hào)通常被看作是周期性沖擊信號(hào)，其與系統(tǒng)響應(yīng)的卷積即為傳感器測(cè)量的輸出信號(hào)。若不考慮噪聲的影響，可以尋找一個(gè)逆濾波器，通過(guò)傳感器輸出的觀測(cè)信號(hào)，利用解卷積恢復(fù)出原始的故障信號(hào)。Wiggin[4]結(jié)合熵值最小卷積濾波思想，提出最小熵解卷積方法(Minimum Entropy Deconvolution，MED)，在解卷積過(guò)程中使得峭度最大化，該方法又稱為最大峭度解卷積。但MED只能處理單一沖擊信號(hào)，并不適合旋轉(zhuǎn)機(jī)械的周期性沖擊信號(hào)的處理。

故障信號(hào)通過(guò)EMD分解[5]過(guò)程中存在端點(diǎn)效應(yīng)以及模態(tài)混疊現(xiàn)象。為了解決這些問(wèn)題，學(xué)者們提出了不同的解決方案。陳哲等[6]提出一種集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)[7]與多尺度排列熵(Multiscale Permutation Entropy，MPE)的艦船輻射噪聲復(fù)雜度特征提取方法，解決復(fù)雜海洋環(huán)境中艦船輻射噪聲的特征提取問(wèn)題。耿讀艷等[8]為了有效識(shí)別心沖擊(Ballistocardiogram)信號(hào)，利用噪聲自適應(yīng)完備總體平均經(jīng)驗(yàn)?zāi)B(tài)分解(Complete EEMD with Adaptive Noise，CEEMDAN)的自適應(yīng)降噪優(yōu)勢(shì)，提出一種基于自適應(yīng)噪聲的CEEMDAN聯(lián)合排列熵(PE)的BCG降噪方法，降噪后信號(hào)的幅頻特性得到明顯改善且信噪比較傳統(tǒng)方法有明顯提高。趙榮珍[9]提出一種基于經(jīng)驗(yàn)小波變換、多尺度排列熵、GG聚類算法相結(jié)合的故障診斷方法。劉興教[10]通過(guò)EEMD和最大相關(guān)峭度解卷積(Maximum Correlated Kurtosis Deconvolution，MCKD)提取柔性薄壁軸承故障特征。

由于以上方法不能產(chǎn)生適當(dāng)?shù)男D(zhuǎn)分量的余量，且在信號(hào)分解過(guò)程中為了降低噪聲及模態(tài)混疊，在模態(tài)分解過(guò)程引入了隨機(jī)噪聲，但是引入隨機(jī)噪聲后不能從根本上消除模態(tài)混疊，并可能會(huì)產(chǎn)生虛假分量。為了解決噪聲及模態(tài)混疊對(duì)故障信號(hào)的不利影響，本文提出一種改進(jìn)的固有時(shí)間尺度分解(Intrinsic Time-scale Decomposition, ITD)與最大相關(guān)峭度解卷積(MCKD)[11]相結(jié)合的故障信號(hào)特征提取方法。首先通過(guò)ITD將原始信號(hào)分解為一系列固有旋轉(zhuǎn)分量和一個(gè)剩余的單調(diào)趨勢(shì)分量，然后選擇與原信號(hào)相似度高的旋轉(zhuǎn)分量作為濾波信號(hào)的主要成分，利用粒子群算法優(yōu)化MCKD的參數(shù)，使得峭度值最大，最后對(duì)信號(hào)進(jìn)行包絡(luò)解調(diào)，提取滾動(dòng)軸承信號(hào)的故障特征。

1 改進(jìn)的固有時(shí)間尺度分解

固有時(shí)間尺度分解(ITD)是由Frei等[12]提出的一種非平穩(wěn)信號(hào)的時(shí)頻分析方法。ITD可以實(shí)時(shí)提取信號(hào)的瞬時(shí)特征，將信號(hào)分解為一系列適當(dāng)?shù)男D(zhuǎn)分量和一個(gè)趨勢(shì)分量。對(duì)于t>0的待分解信號(hào)Xt，ITD的分解過(guò)程是定義一個(gè)運(yùn)算算子L，從Xt提取出基線信號(hào)Lt，再通過(guò)Xt-Lt得到旋轉(zhuǎn)分量Ht，用公式表示為：

Xt=LXt+(1-L)Xt=Lt+Ht

(1)

ITD分解[13]首先找出Xt上所有的局部極值點(diǎn)所對(duì)應(yīng)的時(shí)間集合{τk,k=1,2,3…}，時(shí)間集合包含信號(hào)起始時(shí)間點(diǎn)。為了表示方便，用Xk和Lk代替X(τk)和L(τk)。假設(shè)Ht和Lt在時(shí)間區(qū)間[0,τk]上有定義，Xt中時(shí)間t∈[0,τk+2]，這樣就可以在時(shí)間(τk,τk+1]上定義一個(gè)分段線性的基提取算子L：

t∈(τk,τk+1]

(2)

其中：

(1-α)Xk+1

(3)

式中，α∈(0,1)，通常取值為0.5。以這種方式構(gòu)造的基線信號(hào)Lt，可以使其在極值之間保持單調(diào)性，并在極值點(diǎn)的包絡(luò)內(nèi)。由式(2)和式(3)得到基線點(diǎn)信號(hào)后，通過(guò)對(duì)其進(jìn)行線性插值，得到基信號(hào)Lt。用待分解信號(hào)減去基信號(hào)Lt，就得到一個(gè)適當(dāng)?shù)男D(zhuǎn)分量Ht。將提取出的基信號(hào)Lt作為原始信號(hào)，繼續(xù)進(jìn)行ITD分解，直到所分解的基信號(hào)小于某設(shè)定值時(shí)，ITD分解結(jié)束。最終ITD分解可以表示為：

H1+H2+…HP+LP

(4)

式中，Hi表示第i(i=1,2,…p)個(gè)適當(dāng)?shù)男D(zhuǎn)分量，而Lp表示P次后的殘余基線分量。

由于ITD分解過(guò)程中采用線性插值，線性插值只顧及其附近兩點(diǎn)的影響，函數(shù)上兩點(diǎn)之間的近似隨著所近似函數(shù)的二階導(dǎo)數(shù)增大而逐漸變差，為了解決這個(gè)問(wèn)題，本文采用Akima插值方法。Akima插值方法是一種5點(diǎn)求導(dǎo)分段三次多項(xiàng)式插值算法，其5點(diǎn)求導(dǎo)法利用所求點(diǎn)與前后連續(xù)兩點(diǎn)組成5個(gè)點(diǎn)，然后利用點(diǎn)之間延長(zhǎng)線得到交點(diǎn)，最后得到所求點(diǎn)的導(dǎo)數(shù)[14]。

為了驗(yàn)證基于Akima插值ITD分解的有效性，給出一個(gè)合成的數(shù)字信號(hào)，其數(shù)學(xué)表達(dá)式為：

y(t)=y1(t)+y2(t)+0.07·random

(5)

其中random為隨機(jī)信號(hào)：

y1(t)=[1+0.5cos(8πt)]·cos[200πt+

2cos(10πt)]

y2(t)=0.8sin(πt)sin(30πt)

圖1為原始合成信號(hào)及其分量波形，圖2為通過(guò)EMD分解得到前3個(gè)主要IMF分量的時(shí)域波形，圖3為通過(guò)ITD分解后得到的特征波形?？梢钥吹絀TD分解方法要優(yōu)于EMD分解方法，驗(yàn)證了該方法的有效性。

圖1 原始合成信號(hào)及其分量Fig.1 Original composite signal and its components

圖2 合成信號(hào)EMD分解后前3個(gè)IMF分量Fig.2 The first three IMF components ofcomposite signal after EMD decomposition

圖3 合成信號(hào)ITD分解后前3個(gè)分量Fig.3 The first three components of the composite signal after ITD decomposition

2 最大相關(guān)峭度解卷積

傳感器采集到的滾動(dòng)軸承故障信號(hào)x(n)可以表示為：

x(n)=h(n)·y(n)+e(n)

(6)

式中，h(n)為系統(tǒng)的沖擊響應(yīng)，y(n)為系統(tǒng)故障信號(hào)，e(n)為噪聲。在不考慮噪聲的情況下，尋找最優(yōu)濾波器f=[f1,f2,…fL]T，對(duì)系統(tǒng)采集的信號(hào)x(n)進(jìn)行解卷積處理，可還原出故障信號(hào)y(n)，即

y(n)=f(n)·x(n)

(7)

熵表示系統(tǒng)內(nèi)信息的混亂程度，熵值越大則系統(tǒng)越難預(yù)測(cè)。MED方法將熵值最小與卷積濾波思想進(jìn)行了結(jié)合，對(duì)傳感器采集到的原始信號(hào)進(jìn)行解卷積，以反卷積信號(hào)x(n)峭度作為目標(biāo)函數(shù)求解最優(yōu)濾波器，以突出信號(hào)的沖擊成分。但是濾波器選擇不當(dāng)將會(huì)提取出虛假成分，因此出現(xiàn)了以相關(guān)峭度(Correlation Kurtosis，CK)為最優(yōu)指標(biāo)的MCKD，相關(guān)峭度的計(jì)算式為：

(8)

(9)

式中L為濾波器長(zhǎng)度，而f=[f1,f2,…fL]T。對(duì)目標(biāo)函數(shù)求導(dǎo)，使得：

(10)

最終得到濾波器為：

(11)

其中：

(12)

r=[0,T,2T,…,mT]

(13)

(14)

從以上分析可知，通過(guò)MCKD將傳感器采集信號(hào)x(n)進(jìn)行解卷積濾波時(shí)，需要設(shè)定周期T、偏移量M、濾波器長(zhǎng)度L以及迭代次數(shù)N。

3 試驗(yàn)驗(yàn)證

實(shí)驗(yàn)數(shù)據(jù)來(lái)源于某變速器軸承的外圈故障，變速器轉(zhuǎn)速為n=923 r/min，軸承節(jié)徑D=57.50 mm，滾動(dòng)體直徑d=14.22 mm，壓力角α=0°，滾動(dòng)體個(gè)數(shù)Z=7。外圈故障的特征頻率滿足以下關(guān)系：

(15)

通過(guò)上式可以計(jì)算出軸承外圈故障頻率fr=17.43 Hz。

采集變速器軸承外圈故障信號(hào)，采樣頻率為40 kHz，采集數(shù)據(jù)長(zhǎng)度為84×103。首先對(duì)信號(hào)做去均值處理，處理后信號(hào)的時(shí)域與頻域特征如圖4所示。

圖4 原始故障信號(hào)時(shí)域及頻域特征Fig.4 Time domain and frequency domain characteristics of the original fault signal

圖5為原始信號(hào)EMD分解后一階IMF分量的包絡(luò)譜，由圖5可以看出，故障信號(hào)的特征頻率在35.4 Hz，105 Hz處都有較大的峰值，但是由于噪聲的影響，不能判斷是否含有沖擊成分，也就是故障信號(hào)的特征包絡(luò)譜。所以通過(guò)對(duì)原始信號(hào)的包絡(luò)譜分析，不能有效提取故障特征頻率。

圖5 信號(hào)EMD分解后一階IMF信號(hào)包絡(luò)譜Fig.5 First-order IMF signal envelope spectrum after signal EMD decomposition

接下來(lái)需要去除系統(tǒng)的噪聲，突出損傷沖擊信號(hào)特征。對(duì)信號(hào)進(jìn)行基于Akima插值ITD分解，分解后前4個(gè)分量的時(shí)域波形和其頻譜如圖6所示。

圖6 ITD分解后各分量波形及頻域特征Fig.6 Waveform and frequency domain characteristics of each component after ITD decomposition

ITD分解后各分量的能量比、峭度、峭度因子和相似性的統(tǒng)計(jì)特征如圖7所示。由圖7可知，第一階次的能量比最大，是信號(hào)的主要成分；第一、二階次的峭度最大，證明其變化最為劇烈；通過(guò)相似性分析可知第一階次與原始信號(hào)最接近，故將第一階分量作為提取信號(hào)。

圖7 ITD分解后各分量的統(tǒng)計(jì)特征Fig.7 Statistical characteristics of each component after ITD decomposition

通過(guò)ITD分解得到信號(hào)的主要成分后，再對(duì)信號(hào)進(jìn)行MCKD處理以得到信號(hào)的沖擊成分。通過(guò)MCKD對(duì)信號(hào)進(jìn)行解卷積時(shí)，需要設(shè)定周期T、偏移量M、濾波器長(zhǎng)度L及迭代次數(shù)N。其中：周期T計(jì)算值并不一定最優(yōu)，與實(shí)際還有偏差；而濾波器長(zhǎng)度L直接影響著整個(gè)計(jì)算存儲(chǔ)數(shù)據(jù)大小；偏移量M一般取1～7之間；迭代次數(shù)太小有可能不是最優(yōu)解，太大影響計(jì)算效率，一般選取30。以上參數(shù)的選取直接影響著峭度，為了獲得最優(yōu)參數(shù)，本文以最大峭度為目標(biāo)函數(shù)，采用粒子群算法(PSO)對(duì)上述參數(shù)進(jìn)行優(yōu)化，選取最優(yōu)組合。最終選擇濾波器尺寸L=100、最終迭代次數(shù)N=30、解卷積周期為T(mén)=30、偏移量M=5時(shí)取得最優(yōu)解。圖8比較了MCKD處理前后的信號(hào)譜峭度，圖9為MCKD處理后信號(hào)的包絡(luò)譜。圖9清晰顯示出故障頻率為17.09，外圈故障頻率及其倍頻上出現(xiàn)峰值，能夠清晰地診斷出設(shè)備的故障特征。

圖8 MCKD處理前后信號(hào)譜峭度Fig.8 Signal spectral kurtosis before and after MCKD processing

圖9 MCKD處理后信號(hào)的包絡(luò)譜Fig.9 Envelope spectrum of signal processed by MCKD

4 結(jié)束語(yǔ)

對(duì)于非平穩(wěn)性、非線性的軸承故障信號(hào)，本文提出了基于固有時(shí)間尺度分解和相關(guān)峭度解卷積相結(jié)合的特征提取方法。借助固有尺度分解，提取滾動(dòng)軸承故障信號(hào)的主要特征，為了使得到的波形更加光滑平穩(wěn)，在尺度分解中采用Akima插值方法代替線性插值。借助MCKD提取信號(hào)的沖擊成分過(guò)程中，為了選擇合適的參數(shù)，通過(guò)粒子群算法對(duì)參數(shù)進(jìn)行優(yōu)化，以最大峭度為目標(biāo)函數(shù)，得到最優(yōu)參數(shù)，最后對(duì)MCKD后信號(hào)進(jìn)行包絡(luò)譜分解，可以清晰地提取出故障特征頻率。實(shí)驗(yàn)結(jié)果與理論一致，有利于對(duì)噪聲環(huán)境下軸承的狀態(tài)進(jìn)行特征提取。