基于MODWPT 平方包絡(luò)峭度譜的軸承聲信號故障診斷方法

李方烜

（中國鐵道科學(xué)研究院集團有限公司 機車車輛研究所，北京 100081）

高速動車組軸箱軸承作為走行部關(guān)鍵旋轉(zhuǎn)部件，時常處于高轉(zhuǎn)速、重載荷和多沖擊的工作環(huán)境中，對其進行服役狀態(tài)監(jiān)測和故障診斷具有重要意義。當(dāng)軸承局部發(fā)生缺陷時，無論缺陷位于軸承的外圈、內(nèi)圈還是滾動體，每當(dāng)軸承元件表面接觸到缺陷部分時，均會引起接觸面的彈性沖擊，同時向外發(fā)射聲信號［1-3］，滾動軸承故障軌旁聲學(xué)診斷系統(tǒng)（簡稱TADS）［4］便是利用軌旁聲學(xué)傳感器拾取軸承聲信號并隨后采用一定的信號分析技術(shù)實現(xiàn)對滾動軸承的故障診斷。

Huang 等人最早提出了經(jīng)驗?zāi)B(tài)分解方法（EMD）［5］，EMD 是將原始非平穩(wěn)信號分解為若干個經(jīng)驗?zāi)B(tài)函數(shù)（IMF）分量信號之和，已有諸多領(lǐng)域利用該方法進行濾波、降噪、時頻分析、特征提取等，但是模態(tài)混疊、端點效應(yīng)一直無法避免。隨后Dragomiretskiy 等人提出了變分模態(tài)分解方法（VMD）［6-7］，該方法將非平穩(wěn)信號分解為非遞歸、變分形式的分解模態(tài)，解決了EMD 問題，但是又引入了懲罰因子α且分解的模態(tài)個數(shù)k需要事先給定，k和α的確定取決于其他算法的支持或研究人員的經(jīng)驗，不同參數(shù)選擇組合對分解結(jié)果影響很大。

最大重疊離散小波包變換（MODPWT）［8-9］可以有效解決上述方法在信號分解過程中出現(xiàn)的問題，非常適合對非平穩(wěn)信號的研究。此外由ANTONI 等［10-12］提出的峭度指標(biāo)對信號中的沖擊成分比較靈敏，可以有效檢測故障信息。因此文中將介紹利用基于MODWPT 平方包絡(luò)峭度譜的軸承故障診斷方法，并通過實測軸承聲信號驗證該方法的可行性和有效性。

1 基本理論

1.1 MODWPT 原理

最大重疊小波包變換（MODWPT）作為一種二叉樹結(jié)構(gòu)濾波器組，可將原始非平穩(wěn)信號在不同尺度下劃分為由不同中心頻率（frequency central）和帶寬（band width）組成的子頻帶，具有很好的局部時頻分辨特性，其原理如下：

對于離散時間列向量序列X｛Xt：t=0，…，N-1｝進行離散小波變換，小波濾波器hl（l=0，…，L-1）和尺度濾波器gl（l=0，…，L-1）存在的關(guān)系為式（1）、式（2）：

L為濾波器長度，假設(shè)L≤N，小波濾波器與其偶數(shù)平移長度濾波器正交，為式（3）：

定義MODWPT 的小波濾波器和尺度濾波器為式（4）、式（5）：

給出傳遞函數(shù)式為（6）、式（7）：

其小波濾波器滿足為式（8）：

將序列X用MODWPT 法分解為分解系數(shù)Wj，n（n=0，…，2j-1，j=1，…，J）；j為分解層數(shù)；n為node，隨j變化，如圖1 所示。我們令W0，0≡X，根據(jù)上文給出的尺度濾波器及其傳遞函數(shù)(f)得到式（9）：

圖1 MODWPT 分解方法

同樣用小波濾波器及其傳遞函數(shù)(f)得到式（10）：

注意到分解系數(shù)W1，0和W1，1的長度均為N，為避免抽采樣導(dǎo)致下一層數(shù)據(jù)點減少，對于j≥1，我們向尺度濾波器和小波濾波器分別插入2j-1-1 個0，為式（11）、式（12）：

采樣該循環(huán)，進行塔式計算得到MODWPT 分解系數(shù)，為式（13）：

注：mod 表示兩數(shù)取余。

1.2 平方包絡(luò)峭度譜

軸承的故障信息通常由原信號中的沖擊信號構(gòu)成，而這些沖擊信號并非總位于某一固定頻帶，因此需要在原信號中找到這些沖擊信號才能隨后對軸承進行故障診斷。平方包絡(luò)峭度譜可以有效檢測由MODWPT 分解得到的各個子頻帶中的沖擊成分，從而確定故障頻帶，其原理如下：

對每一個MODWPT 分解系數(shù)做平方包絡(luò)SE(Wj，n，t)可表示為式（14）：

峭度指標(biāo)對瞬態(tài)沖擊比較敏感，其定義為給定信號的四階中心矩與二階中心矩平方之比，對于長度為N的X（t）信號，其峭度可表示為式（15）：

式中：E·[ ]為期望計算符號，因此對于MODWPT分解系數(shù)Wj，n，t的平方包絡(luò)峭度譜可表示為式（19）：

2 算法流程

文中介紹的方法首先將利用MODWPT 對非平穩(wěn)信號進行分解，得到不同尺度下一系列瞬時頻率和瞬時幅值分量之和，再對信號分解的各個頻帶做平方包絡(luò)峭度計算，得到平方包絡(luò)峭度譜，隨后依據(jù)平方包絡(luò)峭度譜中最大值確定原信號中的目標(biāo)頻帶，接下來對目標(biāo)頻帶進行帶通濾波，最后通過包絡(luò)解調(diào)［13］對軸承進行故障診斷，方法如圖2 所示。接下來將采用Python 語言進行算法開發(fā)，并利用Matplotlib 進行繪圖。

圖2 軸承聲學(xué)故障診斷方法

3 實測數(shù)據(jù)驗證

為驗證文中提出的軸承聲學(xué)故障診斷方法，下面將通過2 個在動車組運用過程中產(chǎn)生的故障軸箱軸承樣本進行故障診斷應(yīng)用，由于未找到實際運用產(chǎn)生的軸承內(nèi)圈、滾動體故障樣本，因此2個樣本均為外圈故障軸箱軸承。

3.1 樣本一實測數(shù)據(jù)分析

3.1.1 樣本一概述

樣本一選取A 型動車組外圈故障軸箱軸承，分別對樣本一3 次測量的軸承聲學(xué)數(shù)據(jù)進行分析，相關(guān)信息見表1。

表1 樣本一3 次測量信息

依據(jù)外圈故障特征頻率公式可以對樣本一3 次測量時的外圈故障理論特征頻率進行計算，公式為式（20）：

式中：Z為滾動體數(shù)目；d為滾動體直徑；D為軸承節(jié)徑；fr為轉(zhuǎn)頻；α為滾動體接觸角。

計算結(jié)果見表2。

表2 樣本一理論故障特征頻率 單位：Hz

3.1.2 樣本一第1 次測量數(shù)據(jù)分析

藍色信號為樣本一第1 次測得的原始聲信號如圖3（a）所示，約0.5 s。利用文中介紹的方法計算得到原信號各子頻帶的平方包絡(luò)峭度譜如圖3（b）所示，可以觀察到峭度最大值位于中心頻率fc=9 544 Hz、帶寬bw=707 Hz 處。接著對該目標(biāo)頻帶進行帶通濾波，得到濾波信號，h（l=9 544，h=707）如圖3（a）中的紅色信號所示，可以觀察到濾波信號中存在一系列沖擊成分，說明該頻帶包含軸承故障信息，此外還可以觀察到目標(biāo)頻帶中的沖擊信號幾乎淹沒在原信號的強噪聲中，這表明文中提供的方法在強噪聲干擾環(huán)境下依舊可以有效定位原信號中的故障頻帶。最后對濾波信號進行包絡(luò)解調(diào)分析，在濾波信號包絡(luò)譜中可以找到67.65 Hz 的基頻及其2、3 倍頻，且頻譜中的基頻與計算得到的外圈故障理論特征頻率接近如圖3（c）所示，由此可驗證該軸承發(fā)生了外圈故障。

圖3 樣本一第1 次測量數(shù)據(jù)

3.1.3 樣本一第2 次測量數(shù)據(jù)分析

藍色信號為樣本一第2 次測得的原始聲信號如圖4（a）所示，約0.63 s。同樣對原信號進行計算得到各子頻帶的平方包絡(luò)峭度譜如圖4（b）所示，可以觀察到峭度最大值位于中心頻率fc=10 251 Hz、帶寬bw=707 Hz 處。接著對該目標(biāo)頻帶進行帶通濾波，得到濾波信號，h（l=10 251，h=707）如圖4（a）中的紅色信號所示，可以觀察到目標(biāo)頻帶中的沖擊信號最大幅值相比第1 次已有較大增加，這是由于故障軸箱軸承此時又運行了約19 002 km，隨著軸箱軸承損傷加重，包含軸箱軸承故障信息的沖擊信號也愈發(fā)明顯。進一步對濾波信號進行包絡(luò)解調(diào)分析，在濾波信號包絡(luò)譜中可以找到54.37 Hz 基頻及其2、3 倍頻，且基頻與計算得到的外圈故障理論特征頻率接近，再次驗證該軸承發(fā)生了外圈故障，如圖4（c）所示。此外，在包絡(luò)譜中也可以發(fā)現(xiàn)故障特征頻率的倍頻成分幅值有明顯增加，說明此時軸承外圈故障相較前一次更加嚴重。

圖4 樣本一第2 次測量數(shù)據(jù)

3.1.4 樣本一第3 次測量數(shù)據(jù)分析

藍色信號為樣本一第3 次測得的原始聲信號如圖5（a）所示，約0.57 s。可以觀察到原信號中已出現(xiàn)一系列可見的沖擊成分，對原信號進行平方包絡(luò)峭度譜分析如圖5（b）所示，可以觀察到峭度最大值位于中心頻率fc=10 251 Hz、帶寬bw=707 Hz處。接著對該目標(biāo)頻帶進行帶通濾波，得到濾波信 號，h（l=10 251，h=707）如 圖5（a）中的紅色信號所示，可以觀察到目標(biāo)頻帶中已存在十分清晰的周期性沖擊信號，此時相較該軸箱軸承第1 次測量時又運行了約28 515 km，軸箱軸承損傷已嚴重劣化。繼續(xù)對濾波信號進行包絡(luò)解調(diào)分析，在濾波信號包絡(luò)譜中可以很明顯找到60.16 Hz 基頻及其2、3、4 倍頻，且基頻與計算得到的外圈故障理論特征頻率接近如圖5（c）所示，驗證該軸承發(fā)生了較嚴重的外圈故障。

圖5 樣本一第3 次測量數(shù)據(jù)

3.2 樣本二實測數(shù)據(jù)分析

3.2.1 樣本二概述

樣本二選取B 型動車組外圈故障軸箱軸承，樣本二與樣本一非同型號軸承，由此盡可能保證兩樣本的差異性，從而驗證本方法的普適性。

用相同的方法對樣本二的3 段軸承聲學(xué)數(shù)據(jù)進行分析，相關(guān)信息見表3。

表3 樣本二3 次測量信息

同樣依據(jù)公式（20）可得到樣本二的故障特征頻率見表4。

表4 樣本二理論故障特征頻率 單位：Hz

3.2.2 樣本二第1 次測量數(shù)據(jù)分析

藍色信號為樣本二第1 次測得的原始聲信號如圖6（a）所示，約0.5 s。利用文中介紹的方法計算得到原信號各子頻帶的平方包絡(luò)峭度譜如圖6（b）所示，可以觀察到峭度最大值位于中心頻率fc=10 605 Hz、帶寬bw=1 414 Hz 處。接著對該目標(biāo)頻帶進行帶通濾波，得到濾波信號，h（l=10 605，h=1 414）如圖6（a）中的紅色信號所示，可以觀察到濾波信號中存在沖擊成分，說明該頻帶包含軸承故障信息。接下來對濾波信號進行包絡(luò)解調(diào)分析，在濾波信號包絡(luò)譜中可以找到67.62 Hz 的基頻及其2 倍頻，頻譜中的基頻與計算得到的外圈故障理論特征頻率接近，如圖6（c）所示，由此可驗證該軸承發(fā)生了外圈故障。

圖6 樣本二第1 次測量數(shù)據(jù)

3.2.3 樣本二第2 次測量數(shù)據(jù)分析

藍色信號為樣本二第2 次測得的原始聲信號如圖7（a）所示，約0.82 s。同樣對原信號進行計算得到各子頻帶的平方包絡(luò)峭度譜如圖7（b）所示，可以觀察到峭度最大值位于中心頻率fc=9 544 Hz、帶寬bw=707 Hz 處。接著對該目標(biāo)頻帶進行帶通濾波，得到濾波信號，h（l=9 544，h=707）如 圖7（a）中的紅色信號所示，可以觀察到目標(biāo)頻帶中的存在較明顯的一系列沖擊成分。進一步對濾波信號進行包絡(luò)解調(diào)分析，在濾波信號包絡(luò)譜中可以找到40.86 Hz 基頻及其2、3 倍頻，且基頻與計算得到的外圈故障理論特征頻率接近，如圖7（c）所示，驗證該軸承發(fā)生了外圈故障。相比第1 次得到的包絡(luò)譜，此時包絡(luò)譜中不僅可以找到3 倍頻，也可以發(fā)現(xiàn)倍頻成分更加明顯，說明故障信號特征較前一次愈發(fā)明顯，這與樣本二較前一次繼續(xù)運行30 936 km 后損傷加重的情況符合。

圖7 樣本二第2 次測量數(shù)據(jù)

3.2.4 樣本二第3 次測量數(shù)據(jù)分析

藍色信號為樣本二第3 次測得的原始聲信號如圖8（a）所示，約0.84 s。對原信號進行平方包絡(luò)峭度譜分析如8（b）所示，可以觀察到峭度最大值位于中心頻率fc=10 251 Hz、帶寬bw=707 Hz 處。接著對該目標(biāo)頻帶進行帶通濾波，得到濾波信號，h（l=10 251，h=707）如圖8（a）中的紅色信號所示，可以觀察到目標(biāo)頻帶中已存在十分清晰的周期性沖擊信號，此時相較樣本二第1 次測量時又運行了約49 870 km，軸箱軸承損傷已很嚴重。對濾波信號進行包絡(luò)解調(diào)分析，在濾波信號包絡(luò)譜中可以很明顯找到39.31 Hz 基頻及其2、3、4 倍頻，且基頻與計算得到的外圈故障理論特征頻率接近，如圖8（c）所示，驗證該軸承發(fā)生了較嚴重的外圈故障。

圖8 樣本二第3 次測量數(shù)據(jù)

4 結(jié)論

通過實測數(shù)據(jù)驗證，文中提供的方法可對故障軸承進行有效診斷，研究結(jié)果表明：

（1）基于MODWPT 的平方包絡(luò)峭度法可以有效提取軸承聲信號中的故障頻帶，且能應(yīng)對一定程度的強噪聲干擾。

（2）利用不同故障軸承樣本進行多次算法驗證，發(fā)現(xiàn)該方法具有不錯的普適性，能夠?qū)收陷S承進行穩(wěn)定識別。

（3）對比同一故障軸承樣本在不同服役里程下的故障診斷結(jié)果，可以發(fā)現(xiàn)其能較準確地反映軸承損傷程度變化，說明該方法應(yīng)用于軸承故障診斷時的有效性和準確性。