第一作者代士超男，碩士，1989年生

通信作者郭瑜男，博士,教授，博士生導(dǎo)師，1971年12月生

基于子頻帶譜峭度平均的快速譜峭度圖算法改進(jìn)

代士超，郭瑜，伍星，那靖

(昆明理工大學(xué)機(jī)電工程學(xué)院，昆明650500)

摘要：快速譜峭度圖(Fast Kurtogram)算法具有能自適應(yīng)選取共振解調(diào)頻帶并實(shí)現(xiàn)包絡(luò)解調(diào)提取的優(yōu)點(diǎn)，在滾動(dòng)軸承包絡(luò)分析中有廣闊的應(yīng)用前景，但其在實(shí)際應(yīng)用中，若被采集信號(hào)中包含有較高峰值的脈沖干擾時(shí)，將可能導(dǎo)致譜峭度圖的自適應(yīng)共振帶確定失效，最終導(dǎo)致無(wú)法獲得包含有效滾動(dòng)軸承的故障特征信息的包絡(luò)信號(hào)。為解決快速譜峭度圖算法的上述不穩(wěn)定問(wèn)題，本文提出了一種基于子頻帶譜峭度平均的改進(jìn)快速譜峭度圖算法，其可有效消除或削弱脈沖干擾成分對(duì)譜峭度圖結(jié)果的影響，提高了共振解調(diào)頻帶確定的魯棒性，實(shí)現(xiàn)了基于快速譜峭度圖算法的滾動(dòng)軸承故障特征準(zhǔn)確提取。仿真和試驗(yàn)結(jié)果驗(yàn)證了本方法的有效性。

關(guān)鍵詞：譜峭度，快速譜峭度圖算法，包絡(luò)分析，譜峭度平均

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(51365023)

收稿日期：2013-12-19

中圖分類(lèi)號(hào):TH133.33

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2015.07.016

Abstract:The fast kurtogram algorithm has the advantages of choosing resonance demodulation frequency band adaptively and being able to realize the envelope demodulation extraction, so, it has broad application prospect in the rolling bearing envelope analysis. But in the fast kurtogram calculation, the existense of a high peak pulse interference in the vibration signal may result in that, the fast kurtogram calculation fails to determine the adaptive resonance frequency band, eventually the envelope may not be effectively extracted and used as the rolling bearing fault characteristic information. In order to solve the instability problem of fast kurtogram algorithm, an improvement fast kurtogram algorithm based on the average of sub-band spectral kurtosis was presented. The method can effectively eliminate or weaken the influence of pulse interference on the result of fast kurtogram, improve the robustness of the algorithm and make better the stability of the frequency band selection for resonant demodulation. The method proposed can extract the rolling bearing fault feature accurately. The simulation and experimental results veried its effectiveness.

Improvement on fast kurtogram algorithm based on sub-frequency-band spectral kurtosis average

DAIShi-chao,GUOYu,WUXing,NAJing(Faculty of Mechanical and Electrical Engineering, Kunming University of Science and Technology, Kunming 650500, China)

Key words:spectral kurtosis; fast kurtogram; envelope analysis; spectral kurtosis on average

包絡(luò)分析(又稱(chēng)為共振解調(diào)分析)是目前被廣泛應(yīng)用的滾動(dòng)軸承初期故障特征提取技術(shù)之一。包絡(luò)分析的關(guān)鍵在于確定合適的共振解調(diào)頻帶參數(shù)(濾波中心頻率和帶寬)。近來(lái)提出的基于譜峭度的包絡(luò)分析方法[1]，能夠根據(jù)譜峭度值自適應(yīng)地獲取優(yōu)化的濾波中心頻率和帶寬，并實(shí)現(xiàn)對(duì)故障特征的有效提取，其已被證明是一個(gè)非常有效的滾動(dòng)軸承故障包絡(luò)分析方法[2]。但研究[3]表明，當(dāng)拾取信號(hào)中包含有較高峰值的脈沖干擾時(shí)，常規(guī)的快速譜峭度圖算法(Fast Kurtogram)[4]將可能無(wú)法確定合理的共振解調(diào)頻帶[3]，進(jìn)而無(wú)法有效獲得滾動(dòng)軸承故障信息。文獻(xiàn)[3]中提出了一種無(wú)該問(wèn)題的Protrugram算法，但該算法存在需要預(yù)先確定共振頻帶大概寬度的不足。本文在常規(guī)快速譜峭度圖算法的基礎(chǔ)之上，提出了一種基于子頻帶譜峭度平均的改進(jìn)快速譜峭度圖算法，其可有效消除或削弱脈沖干擾成分對(duì)原譜峭度圖結(jié)果的影響，可實(shí)現(xiàn)存在較高峰值脈沖干擾的條件下滾動(dòng)軸承故障振動(dòng)共振解調(diào)優(yōu)化頻帶的準(zhǔn)確確定，進(jìn)而實(shí)現(xiàn)故障特征的有效提取。仿真和試驗(yàn)結(jié)果驗(yàn)證了本方法的有效性。

1基于Fast Kurtogram的復(fù)包絡(luò)提取簡(jiǎn)介

1.1譜峭度介紹

峭度(Kurtosis)對(duì)沖擊信號(hào)敏感，用于描述振動(dòng)信號(hào)中沖擊的強(qiáng)弱，但易受噪聲干擾而難以表征沖擊的特征。近來(lái)提出的譜峭度 (Spectral Kurtosis，SK)，可反映信號(hào)峭度值隨頻率的變化。計(jì)算公式為[1,5]:

(1)

1.2Fast Kurtogram算法簡(jiǎn)介

Fast Kurtogram 算法的主要思想為通過(guò)構(gòu)建一系列具有不同頻帶的1/3-二叉樹(shù)帶通濾波器組實(shí)現(xiàn)各子(頻)帶的譜峭度計(jì)算，并通過(guò)對(duì)子帶譜峭度值的比較確定優(yōu)化解調(diào)頻帶，進(jìn)而實(shí)現(xiàn)解調(diào)頻帶對(duì)應(yīng)的包絡(luò)的提取[4]。其主要實(shí)現(xiàn)步驟包括：

(1)濾波器的構(gòu)建及復(fù)包絡(luò)的提取。與小波包分解算法相似，首先構(gòu)建兩個(gè)準(zhǔn)解析濾波器，分別為低通濾波器h0(j)和高通濾波器h1(j)：

(2)

式中:h(j)為截止頻率fc=1/8+，≥0的低通濾波器(此處頻率歸一化，即fs=1，一般fc取0.3)；h0(j)是由h(j)頻移1/8后得到的準(zhǔn)解析低通濾波器，其帶寬為[0,1/4]；h1(j)是由h(j)頻移3/8后得到的準(zhǔn)解析高通濾波器，其帶寬為[1/4,1/2]。利用以上兩個(gè)濾波器分別對(duì)振動(dòng)信號(hào)進(jìn)行濾波，并采用2倍降采樣方法進(jìn)行迭代，將原信號(hào)進(jìn)行L級(jí)分解，每級(jí)濾波后信號(hào)的個(gè)數(shù)為2l,把原信號(hào)分解成不同子頻帶下的復(fù)包絡(luò)信號(hào)[4]cli(n)，其中分解級(jí)數(shù)l=0,1,…,L-1，子頻帶信號(hào)位置系數(shù)i=0,1,…,2l-1。分解算法可用以下公式實(shí)現(xiàn)：

(3)

式中:*為卷積運(yùn)算。當(dāng)l=0時(shí)，c0(n)≡x(n)。cli(n)可認(rèn)為是原始振動(dòng)信號(hào)x(t)經(jīng)不同濾波器濾波后得到的子頻帶信號(hào)，濾波器的中心頻率和帶寬分別為fi=(i+1)×2-l-1和fl=2-l-1。根據(jù)包絡(luò)檢波原理，(n) 為第l級(jí)第i個(gè)子頻帶信號(hào)，是信號(hào)的包絡(luò)信號(hào)。為了提高分析的精度，對(duì)頻帶進(jìn)一步細(xì)分，在l+1和l+2級(jí)之間插入32l個(gè)濾波器，構(gòu)建三個(gè)頻帶分別為[0,1/6]、[1/6,1/3]和[1/3,1/2]的準(zhǔn)解析帶通濾波器gj(n)(j=0,1,2)，對(duì)每個(gè)子頻帶信號(hào)進(jìn)行再分解[4]。

(4)

(3)優(yōu)化包絡(luò)參數(shù)及復(fù)包絡(luò)的獲取。當(dāng)式(4)中譜峭度K值最大時(shí)，對(duì)應(yīng)于優(yōu)化的中心頻率fo、帶寬Δfo及復(fù)包絡(luò)信號(hào)c0(n)。即：

(fo,Δfo,co(n))=argmax{K(fi,Δfl)}

(5)

式中：argmax表示取最大值所對(duì)應(yīng)的參數(shù)，即獲得譜峭度值為最大時(shí)對(duì)應(yīng)的優(yōu)化中心頻率和帶寬，以及對(duì)應(yīng)的優(yōu)化復(fù)包絡(luò)信號(hào)。

2Fast Kurtogram算法改進(jìn)

2.1傳統(tǒng)快速譜峭度圖算法的不足

譜峭度對(duì)沖擊敏感，其能表明振動(dòng)信號(hào)中一系列瞬態(tài)的存在并指出它們?cè)陬l域中的位置[7]。因此，利用滾動(dòng)軸承局部故障對(duì)應(yīng)振動(dòng)信號(hào)的沖擊性特征，由Fast Kurtogram算法可以自適應(yīng)地確定優(yōu)化的滾動(dòng)軸承故障引起沖擊的共振解調(diào)頻帶。但應(yīng)注意軸承早期故障對(duì)應(yīng)的振動(dòng)信號(hào)通常為弱信號(hào)，易受其他高峰值脈沖成分的干擾。文獻(xiàn)[3]指出，當(dāng)采集的滾動(dòng)軸承振動(dòng)信號(hào)中包含高峰值脈沖干擾時(shí)，其會(huì)導(dǎo)致某些頻帶的峭度值變大，影響滾動(dòng)軸承故障振動(dòng)所對(duì)應(yīng)共振解調(diào)頻帶的優(yōu)化選擇，進(jìn)而導(dǎo)致包絡(luò)解調(diào)失效。因此，直接用原Fast Kurtogram算法確定優(yōu)化解調(diào)頻帶并不穩(wěn)定，易受沖擊噪聲影響而導(dǎo)致后續(xù)的包絡(luò)分析失敗，本文后面的實(shí)驗(yàn)分析中也將驗(yàn)證這一點(diǎn)。

2.2基于子頻帶譜峭度平均的Fast Kurtogram算法

圖1　基于子頻帶譜峭度平均的快速譜峭度圖算法流程 Fig.1 Process of fast kurtogram based on the average of sub-band spectral kurtosis

眾所周知，頻域平均是常用的提高信噪比的有效手段，其可有效地去除頻譜中非周期的隨機(jī)干擾[8]。其在研究中提示我們可以將平均的概念引入到譜峭度圖中作為降噪手段。為了解決上述脈沖干擾問(wèn)題，在原Fast Kurtogram算法的基礎(chǔ)之上，本文提出了基于子頻帶譜峭度平均的Fast Kurtogram算法。其算法流程如圖1，主要步驟見(jiàn)圖1。

(1)原始采集信號(hào)分段

將原始采集振動(dòng)序列x(n)分為等長(zhǎng)度的M段，得到分段信號(hào)xm(n)，(m=1,2,3,4,…M)。

圖2　第m段數(shù)據(jù)經(jīng)樹(shù)狀濾波器組濾波后 各級(jí)子頻帶峭度值 Fig.2 Kurtosis values by the tree filter group for the mth data section

(2)各段信號(hào)譜峭度計(jì)算

(3)子頻帶譜峭度平均

將步驟(2)中得到的m組譜峭度值數(shù)組進(jìn)行對(duì)應(yīng)位置求和取平均?？杀硎緸椋?/p>

(6)

3仿真試驗(yàn)分析

為驗(yàn)證本方法的有效性，進(jìn)行了仿真試驗(yàn)，其過(guò)程如下：仿真一個(gè)系統(tǒng)固有頻率fr=4 000 Hz幅值為0.1的單一微弱沖擊衰減的振動(dòng)作為軸承故障每轉(zhuǎn)內(nèi)的一次沖擊：

(7)

軸承故障頻率為124 Hz。

仿真由三個(gè)強(qiáng)低頻正弦干擾組成的低頻干擾信號(hào)y(t)，頻率分別為445 Hz、35 Hz和11 Hz，對(duì)應(yīng)的幅值分別為0.3、1.0和0.8。

仿照文獻(xiàn)[3]，仿真一固有頻率f=6 000 Hz的振蕩衰減脈沖p(t)作為脈沖干擾成分：

p(t)=0.75e(-Dt)cos(2πfgt)

(8)

式中：D表示震蕩衰減系數(shù)。

最后仿真一個(gè)幅值為0.05的高斯白噪聲干擾n(t)。最終得到仿真的軸承振動(dòng)信號(hào)為：

x(t)=R(t)+y(t)+p(t)+n(t)

(9)

式中：R(t)表示為仿真的軸承衰減沖擊振動(dòng)信號(hào)。仿真過(guò)程中取采樣頻率fs=25 kHz，采樣時(shí)間為10 s，仿真的滾動(dòng)軸承振動(dòng)信號(hào)及其各分量(圖中顯示時(shí)間為6.7 s~6.77 s)如圖3、圖4，從圖中可以發(fā)現(xiàn)，脈沖干擾成分在整個(gè)信號(hào)中比較明顯。對(duì)仿真的10 s軸承故障振動(dòng)信號(hào)進(jìn)行分析，截取1秒時(shí)長(zhǎng)的信號(hào)段進(jìn)行計(jì)算，有脈沖干擾存在的數(shù)據(jù)段經(jīng)Fast Kurtogram計(jì)算得到濾波中心頻率和帶寬會(huì)失效，如圖5所示。

圖3　仿真軸承時(shí)域信號(hào)(虛線圈內(nèi)為干擾) Fig.3 Time domain simulation of bearing signal (interference in the dotted circles)

圖4　仿真軸承信號(hào)組成分量 Fig.4 Simulation of bearing signal components

圖5　脈沖干擾時(shí)段數(shù)據(jù)計(jì)算得到的譜峭度圖 Fig.5 Kurtogram of the signal with an impulsive interference

圖6　仿真信號(hào)有干擾成分平方包絡(luò) Fig.6 Square envelope of simulation signal with an impulsive interference

圖7　仿真信號(hào)包含沖擊干擾時(shí)的平方包絡(luò)譜 Fig.7 Squared envelope spectrum of the simulation signal with an impulsive interference component

從圖5可以看出譜峭度計(jì)算得到的濾波中心頻率和帶寬分別為fc=6 054.687 5 Hz、B=390.625 Hz，顯然中心頻率與脈沖干擾的固有頻率6 000 Hz相吻合，用該參數(shù)提取平方包絡(luò)[9]的時(shí)間波形如圖6所示，脈沖干擾十分明顯，對(duì)該包絡(luò)曲線進(jìn)行FFT運(yùn)算，得到的平方包絡(luò)譜如圖7所示，不難看出無(wú)法獲得軸承的故障特征頻率譜線，包絡(luò)分析失效。按論文中介紹的改進(jìn)方法進(jìn)行譜峭度平均計(jì)算，按照2.2節(jié)中所述，將10 s信號(hào)平均分為10段，分別計(jì)算譜峭度后再進(jìn)行子頻帶譜峭度平均計(jì)算，如圖8所示，得到的包絡(luò)中心頻率fc=3 906.25 Hz、帶寬B=1 562.5 Hz，與軸承仿真數(shù)據(jù)4 000 Hz的固有頻率基本吻合，用該參數(shù)進(jìn)行平方包絡(luò)提取獲得的時(shí)間波形如圖9所示，對(duì)應(yīng)的平方包絡(luò)譜如圖10所示，與原Fast Kurtogram算法相比，采用本文所提方法計(jì)算得到的包絡(luò)參數(shù)更為合理和準(zhǔn)確。仿真結(jié)果表明本方法的可行性和有效性。

圖8　子頻帶平均后得到的譜峭度圖 Fig.8 Kurtogram after sub-frequency-band average

圖9　仿真信號(hào)譜峭度平均后得到的平方包絡(luò) Fig.9 Square envelope of simulation signal after spectral kurtosis average

圖10　仿真信號(hào)譜峭度平均后包絡(luò)譜 Fig.10 FFT of simulation signal after spectral kurtosis average

4測(cè)試試驗(yàn)分析

測(cè)試研究中以QPZZ-Ⅱ旋轉(zhuǎn)機(jī)械振動(dòng)分析及故障模擬實(shí)驗(yàn)系統(tǒng)為測(cè)試對(duì)象，如圖11。

圖11　QPZZ-Ⅱ旋轉(zhuǎn)機(jī)械振動(dòng)分析及故障模擬實(shí)驗(yàn)系統(tǒng) Fig.11 QPZZ-Ⅱ rotating machinery vibration analysis and fault simulation system

在故障軸承座內(nèi)安裝有一故障滾動(dòng)軸承，進(jìn)行滾動(dòng)軸承故障實(shí)驗(yàn)，以外圈故障為例，試驗(yàn)參數(shù)如下：故障軸承的型號(hào)為NU205EM；安裝軸為恒定轉(zhuǎn)速1 200 r/min(對(duì)應(yīng)轉(zhuǎn)頻為20 Hz)；數(shù)據(jù)采集設(shè)備為NI PXI-1042Q動(dòng)態(tài)信號(hào)采集系統(tǒng)；采樣頻率為20 kHz；采樣時(shí)間為10 s；加速度傳感器為PCB加速度傳感器PCB603C01；安裝位置為軸承座垂直位置和兩側(cè)水平位置；轉(zhuǎn)速脈沖測(cè)量傳感器為電渦流位移傳感器，型號(hào)為4842。圖12為采集的原始振動(dòng)信號(hào)時(shí)間波形，在數(shù)據(jù)段中存在某一未知源的脈沖干擾。

圖12　原始測(cè)試振動(dòng)信號(hào) Fig.12 Original testing vibration signa

將數(shù)據(jù)信號(hào)平均分10段(分段可根據(jù)分析需要或數(shù)據(jù)長(zhǎng)度而定)，分別對(duì)每一段數(shù)據(jù)進(jìn)行Fast Kurtogram計(jì)算，獲取自適應(yīng)濾波中心頻率和帶寬。含有未知源脈沖干擾成分段(如圖13)計(jì)算結(jié)果如圖14所示，得到的濾波中心頻率和帶寬分別為fc=3 281.25 Hz、B=312.5 Hz。用該參數(shù)進(jìn)行平方包絡(luò)提取，其包絡(luò)后時(shí)間信號(hào)如圖15所示，脈沖干擾十分明顯，對(duì)該包絡(luò)進(jìn)行FFT運(yùn)算，如圖16所示，譜線很雜亂，很難獲得軸承的故障特征頻率譜線，難以做出準(zhǔn)確的故障診斷。

圖13　有干擾成分的原始信號(hào)段 Fig.13 Original signal with an impulsive interference

圖14　存在干擾成分信號(hào)段的譜峭度圖 Fig.14 Kurtogram of signal with an impulsive interference

圖15　有干擾成分譜峭度計(jì)算后的平方包絡(luò) Fig.15 Squared envelope of signal with an interference

圖16　原方法在有沖擊干擾成分下得到的平方包絡(luò)譜 Fig.16 Squared envelope spectrum by original method for signal with an impulsive interference

采用本文所提出子頻帶譜峭度平均的方法，將10段信號(hào)的譜峭度值進(jìn)行平均計(jì)算，得到的優(yōu)化濾波中心頻率和帶寬分別為fc=4 375 Hz、B=1 250 Hz，如圖17所示。用該參數(shù)對(duì)原始軸承振動(dòng)信號(hào)進(jìn)行平方包絡(luò)提取，其時(shí)間波形如圖18所示。

圖17　子頻帶平均后得到的譜峭度圖 Fig.17 Kurtogram after sub-frequency-band average

圖18　譜峭度平均后得到的平方包絡(luò) Fig.18 Squared envelope after spectral kurtosis average

從圖18和圖15包絡(luò)曲線的對(duì)比可以看出，采用子頻帶譜峭度平均的方法得到的參數(shù)在進(jìn)行平方包絡(luò)時(shí)可有效削弱脈沖干擾造成的包絡(luò)不理想情況。

對(duì)該包絡(luò)進(jìn)行FFT運(yùn)算，結(jié)果如圖19所示。對(duì)比原Fast Kurtogram算法結(jié)果，采用本文所提出的方法可以在存在脈沖干擾情況下準(zhǔn)確獲得滾動(dòng)軸承的故障特征頻率成分及其倍頻。

圖19　譜峭度平均后得到的平方包絡(luò)譜 Fig.19 Squared envelope spectrum after spectral kurtosis average

5結(jié)論

本文提出的基于子頻帶譜峭度平均的改進(jìn)快速譜峭度圖算法，可有效消除或削弱高峰值干擾成分對(duì)譜峭度計(jì)算的影響，避免信號(hào)中脈沖干擾導(dǎo)致滾動(dòng)軸承優(yōu)化解調(diào)頻帶確定的失效，提高了包絡(luò)分析中濾波中心頻率和帶寬選擇的準(zhǔn)確性，可實(shí)現(xiàn)對(duì)滾動(dòng)軸承故障特征的準(zhǔn)確提取。相比文獻(xiàn)[3]中提出的方法而言，本文所提出的方法不需要確定故障頻率帶寬，可直接對(duì)故障頻率未知的含有高峰值干擾成分振動(dòng)信號(hào)進(jìn)行有效的故障特征提取和包絡(luò)分析。仿真和試驗(yàn)結(jié)果驗(yàn)證了本方法的有效性。

參考文獻(xiàn)

[1]Sawalhi N，Randall R B，Endo H. The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis [J]. Mechanical Systems and Signal Processing, 2007,21：2616-2633.

[2]Antoni J，Randall R B. The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines [J]. Mechanical Systems and Signal Processing, 2006,20(2)：308-331.

[4]Antoni J. Fast computation of the kurtogram for the detection of transient faults [J]. Mechanical Systems and Signal Processing, 2007,21(1)：108-124.

[5]郭瑜，鄭華文，高艷，等. 基于譜峭度的滾動(dòng)軸承包絡(luò)分析研究 [J]. 振動(dòng).測(cè)試與診斷, 2011,31(4)：517-521.

GUO Yu，ZHENG Hua-wen，GAO Yan，et al. Envelope analysis of rolling bearings based on spectral kurtosis [J]. Journal of Vibration，Measurement & Diagnosis，2011，31(4)：517-521.

[6]劉亭偉，郭瑜，李斌，等. 基于譜峭度的滾動(dòng)軸承故障包絡(luò)階比跟蹤分析 [J]. 振動(dòng)與沖擊，2012,31(17)：149-153.

LIU Ting-wei，GUO Yu，LI Bin，et al. Envelope order tracking analysis for rolling element bearing faults based on spectral kurtosis [J]. Journal of Vibration and Shock，2012，31(17)：149-153.

[7]Antoni J. The spectral kurtosis: a useful tool for characterising non-stationary signals[J]. Mechanical Systems and Signal Processing, 2006，20(2)：282-307.

[8]鞠萍華，秦樹(shù)人，秦毅，等. 多分辨 EMD 方法與頻域平均在齒輪早期故障診斷中的研究 [J]. 振動(dòng)與沖擊, 2009，28(5)：97-101.

JU Ping-hua，QIN Shu-ren，QIN yi，et al. Research on gear earlier fault diagnosis by method of multi-resolution empirical mode decomposition and frequency domain averaging[J]. Journal of Vibration and Shock，2009，28(5)：97-101.

[9]Ho D，Randall R B. Optimisation of bearing diagnostic techniques using simulated and actual bearing fault signals [J]. Mechanical Systems and Signal Processing，2000，14(5)：763-788.