基于數學形態學分形理論與PNN的軸承狀態監測與故障診斷

林 懿， 龍 偉， 李炎炎， 張慶華

(四川大學機械工程學院， 成都 610065)

1 引 言

滾動軸承作為旋轉機械中應用最廣泛的關鍵部件，因其惡劣的工作環境，導致故障頻發，給工業安全生產帶來了嚴峻挑戰.《中國制造2025》已明確將提高大型設備的可靠性與安全性列為重要發展方向.因此在工業4.0背景下，對滾動軸承進行基于數據驅動的實時監測與故障診斷具有重要的理論和現實意義.

機電設備運行時產生的振動信號中包含著豐富的軸承狀態信息，且抗干擾能力強，常用于設備故障診斷[1].但由于滾動軸承工況復雜、激勵源多、耦合性強，振動特性表征模糊，診斷效果往往不佳，因此故障特征的精確提取是實現高精度故障診斷的重要前提[2].經典的時頻域處理法存在特征提取混肴、表征性差、自適應性缺乏等缺點，不適用于智能化監測與診斷[3].近年興起的基于深度學習的軸承故障診斷方法，雖然取得了不錯的診斷效果，但因其硬件成本高、運算量大、實時性差等問題，在工業生產方面同樣缺乏使用環境[4].分形幾何學是一門以不規則幾何形態為研究對象的幾何學,通過提取復雜信號的分形維數，定量表征信號的復雜性和非線性.李洋等[5]將分形幾何理論引入軸承故障診斷中，成功提取了不同狀態軸承的特征信息.在多種分形理論中，基于數學形態學的分形維數計算方法以其魯棒性強、計算量少等優點被廣泛應用于信號處理方向[6].李兵等[7]首次使用數學形態學理論計算一維振動信號的分形維數，成功識別軸承四種基本狀態，但由于振動信號中噪聲的影響，加之表征維度低，難以區分不同損傷程度的軸承故障.

在實測振動信號中，大量的背景噪聲因其隨機性，不具備明顯的分形特性，必須對原始信號降噪提純，才能保證分形維數的準確性[8].經驗模態分解(Empirical Mode Decomposition，EMD)通過迭代分解將原始信號重構為一系列本征模態函數(Intrinsic Mode Function， IMF)達到降噪目的，因其自適應性、正交性及非線性被廣泛應用于信號處理領域，但“端點效應”和“模態混疊”的存在影響了其降噪性能.盡管國內外學者對其進行大量的改進與優化，迭代次數多、運算量大及實時性差等問題的存在使其難以適用于軸承狀態的在線監測[9].局部均值分解(Local Mean Decomposition，LMD)于2005年被Jonathan S. Smith提出并首次運用于EEG信號中，它通過自適應的將復雜原始信號分解成一系列純調頻信號與包絡信號的乘積函數(PF分量)降噪提純，不僅延緩了EMD“模態混疊”與“端點效應”的缺陷，并且分解速度更快.張亢等[10]通過LMD分解軸承振動信號降低噪聲影響，并提取相關性系數最大分量的分形維數作為分類依據，基本實現了軸承工作狀態的識別，但軸承工況復雜，使用單一特征表征軸承狀態缺乏空間性與全面性.金成功[11]則通過EMD的改進算法CEEMDAN降低振動信號的噪聲，并提升特征維度，但結果同樣只能識別4種基本狀態，且算法運行時間較長.直到齊詠生等[12]使用LMD算法分解振動信號，通過提取前三個PF分量，將其重構為特征向量進行軸承故障診斷，才對同種故障不同損傷程度狀態的識別取得了較好的效果，不過，因為LMD算法過量的篩選迭代，導致PF虛假分量的存在，使得第三個PF分量與原始信號的相關性較差，難以精確表征，因此，在LMD算法的改進與軸承故障特征的選擇上，仍有進一步提升的可能.

將提取的特征向量與模式識別理論結合，可以實現軸承的在線智能監測與故障診斷[13].其中，聚類分析難以對高維數據分類識別[14]，SVM[15]算法冗余，缺乏效率，且運算時間長.而概率神經網絡(Probabilistic Neural Network，PNN)是一種依據貝葉斯決策理論構建的智能分類器.因其訓練簡潔、收斂快以及容錯高的優點，適用于數據的實時處理[16].

就上述不足，本文通過構建“三次樣條插值”與“迭代終止條件”完善LMD算法的不足，優化算法實時性與PF分量的精確性；同時，結合分形盒維數的基本原理，首次引入物理量“形態學覆蓋面積”作為特征向量，代替相關性系數較差的PF分量，仍舊從多維度表征軸承狀態，提升表征精度.

2 背 景

2.1 數學形態學分形維數

分形維數的計算是用不同尺度表征分形集，以定值其非線性與復雜性，同時量化其結構特征.在一維離散數據的背景下，數學形態學維數的計算是基于盒維數與多尺度數學形態學的覆蓋思想，對一維振動信號在不同尺度下進行腐蝕(Θ)和膨脹(⊕)運算，具體原理如下.

設f(n)與g(m)分別為定義在F={0，1，2，…，N-1}和G={0，1，2，…，M-1}上的離散函數，且N≥M，則數學形態學中腐蝕濾波器(Θ)和膨脹濾波器(⊕)的定義如下.

(fΘg)=min[f(n+m)-g(m)]

(1)

(f⊕g)=max[f(n-m)+g(m)]

(2)

其中，f為輸入信號；g為結構元素.

在尺度ε(整數，且1≤ε≤εmax)下的結構元素定義為ε-1次腐蝕或者膨脹運算，即

gΘε=gΘgΘ…Θg

(3)

g⊕ε=g⊕g⊕…⊕g

(4)

則尺度ε下，對信號的覆蓋面積為

(5)

Maragos于1993年提出，當ε→0時有

(6)

其中，ε=1, 2 …,εmax；DM為Minkowski-Bouligand維數；c為常數；同時，對{ln(Ag(ε)/ε2)，ln(1/ε)}使用最小二乘線性擬合，所得直線斜率則為分形維數DM的極限近似值.

2.2 局部均值分解算法及其改進

2.2.1 基本算法 局部均值分解的本質就是從信號中逐步分離出調頻信號和調幅包絡信號的過程.通過多次篩選，自適應的將PF分量由高頻到低頻逐漸分離，直到剩余信號為常數或不再包含振蕩為止.具體過程的步驟如下.

(1) 計算原始信號x(t)相鄰局部極值點的均值，如式(7).

(7)

將各mi線性連接，再用滑動平滑法處理，得到局部均值函數m11(t).

(2) 局部包絡函數，計算方法如式(8).

(8)

用局部均值函數相同的方法處理，得到包絡估計函數a11(t).

(3) 將m11(t)從x(t)中分離，得到殘差信號h11(t).

h11(t)=x(t)-m11(t)

(9)

(4) 振幅解調.

s11(t)=h11(t)/aii(t)

(10)

我們將s11(t)視為原始信號，重復步驟(1)和步驟(2)得到包絡估計函數a12(t)，若a12=1，則判定s11(t)為純調頻函數.否則，以s11為原始信號，重復上述步驟，直到滿足條件得到純調頻信號s1n(t).不過在實際應用中，通常設定閾值α，使|a1n-1|<α時，終止迭代.

(5) 所有包絡估計函數的乘積為調幅包絡信號a1(t)，即

(11)

(6) 求出乘積函數PF1(t).

PF1=a1t(t)s1n(t)

(12)

(7) 用x(t)減去PF1(t)得到u1(t)函數，重復上述步驟，直到uk(t)變成常數或不再振蕩.

因此，通過局部均值算法的分解，原始信號x(t)如式(13)構成.

(13)

2.2.2 弊端及改進 盡管局部均值在一定程度上改善了經驗模態分解理論支撐缺乏、運算量大以及模式混淆等缺點，但其本身仍存在大量迭代“篩選”，實時性不足.比如，在計算局部均值函數與包絡估計函數的過程中，連續的滑動平移動作會增加運算復雜度，損耗大量的計算內存；其次，過多的篩選次數不僅浪費時間，也容易產生相關性系數極差的虛假PF分量.

針對上述問題，我們提出了以下兩個改善方案.

(1) 采用“三次樣條插值法”代替滑動平移法，只需一次插值，即可獲得平滑的理想結果，減少了算法的復雜度，改善了實時性.

(2) 設定合理的閾值γ，作為PF分量的選擇標準，當其與原始信號的相關性系數符合ρ<γ，則判定為虛假分量，及時終止迭代，不僅降低了算法運行時間，同時保證了分量的合理性，提高了特征提取的準確率.相關性系數ρ計算公式如下.

(14)

其中，x(t)為原始信號；p(t)為PF分量時序信號，并且經過多次調試，選擇0.2作為相關系數閾值，既保證了PF分量的個數，又極大的避免了虛假分量的出現.

2.3 特征提取

對冗余過多的特征向量智能識別，不僅增加硬件負擔，不利于算法實時性，還會影響算法識別的正確率.因此，只提取相關性系數最大的兩個分量的分形維數作為特征向量是合理的.其次，形態學覆蓋面積，作為描述分形集整體空間大小的物理量，同樣可定量且直觀地描述信號的特性.其理論依據基于分形盒維數的覆蓋方法，即假設有一個面積為S的未知平面，用邊長為ε的正方形填充該區域，需要正方形的個數為N(ε)，則有下列推論：

(15)

兩邊取對數：

lnN(ε)=lnS+2 ln (1/ε)

(16)

得到未知平面的維數為：

(17)

因此，形態學覆蓋面積S，相較于分形維數D，同樣可以作為軸承分類依據，此外，在基于數學形態濾波的分形維數提取中，因腐蝕濾波器與膨脹濾波器的特性，形態學覆蓋面積與形態學分形維數同樣具有平移不變性，因此，將式(5)中的形態學覆蓋面積Ag作為特征向量引入軸承故障診斷也具有科學的理論依據.

綜上所述，本文提取相關性系數最大的兩個PF分量的數學形態學分形維數，結合信號的數學形態學覆蓋面積共同構建特征向量，輸入概率神經網絡中以達到狀態識別的目的.

2.4 算法步驟

(1) 信號去噪.使用改進局部均值算法，對采集到的原始信號x(t)去噪，保留相關性系數最大PF1，PF2分量.

(2) 特征提取.提取PF1，PF2的數學形態學分形維數T1，T2并對信號x(t)的形態學覆蓋面積Ag照盒維數覆蓋思想取對數，得到第三個特征向量T3.

(3) 狀態識別.將特征向量T1，T2，T3輸入已經訓練好的概率神經網絡模型中，完成對軸承工作狀態的識別.

3 實 驗

為了驗證該方法在改善實時性與準確率方面的有效性，采用CWRU(美國凱斯西儲大學)軸承數據中心[17]的實測數據.本文選取滾動體、外圈和內圈三種故障進行分析，其中，為了更直觀體現所提方法的實用性，滾動軸承、內、外圈分別采用了損傷直徑為0.007 in和0.021 in(1 in=2.54 cm)的故障信號.因此，包含正常狀態在內，共采用了七種類型共819個樣本的振動信號，每種類型的樣本數量為117，單個樣本中含有1 024個樣本點.數據其他信息為：轉速為1 750 r/min，采樣頻率為12 kHz，故障深度為0.011 in.

3.1 信號分解

滾動軸承工作狀態發生改變時，其振動信號也會產生相應變化.如圖1所示，正常工作狀態的滾動軸承，振動信號隨機性強，呈現出非線性、非平穩的特點；滾動體故障時，有一定頻率的沖擊，但整體隨機性較強；內、外圈故障狀態下，時域信號波形呈現周期性沖擊，特征明顯.因此，軸承四種狀態中，正常狀態與滾動體故障狀態相似，沖擊特性模糊、隨機性強，且振幅較小，特征混雜，兩者不易區分.而內、外圈故障狀態規律強、振幅較大，具有明顯的沖擊特性，但兩者特性相似，仍然不易區分.

圖1 滾動軸承四種狀態時域波形圖Fig.1 Time domain waveform diagram of rolling bearing in four states

振動特征容易淹沒在摩擦、剛度、載荷等復雜工況環境帶來的噪聲中，因此，需要對原始振動信號分解提純.以0.07 in內圈故障為例，分別使用改進局部均值分解(ILMD)與標準局部均值分解(LMD)處理原始信號，分解結果如圖2和圖3，且兩種算法各PF分量的相關系數ρ如表1所示.可以看到，因為閾值γ的緣故，當ILMD分解出的PF分量與原始信號的相關系數ρ小于0.2時，迭代終止，不僅減少了運行時間，也有效的避免了虛假分量的產生；此外，三次樣條插值法的采用，一定程度的改善了“端點效應”的影響，并且顯著提升了前兩階PF分量的相關系數，最大程度保留了原始信號的振動特征，提高了特征提取的精度.同時，為更進一步驗證所提算法的優勢，將其與經驗模態分解(EMD)、自適應噪聲的完全集合經驗模態分解(CEEMDAN)比較.

圖2 ILMD時域圖Fig.2 Temporal waveform of ILMD

圖3 LMD時域圖Fig.3 Temporal waveform of LMD

表1 PF分量相關系數對比

通過對0.007 in軸承內圈故障信號的處理，分別對比前兩階分量相關性系數ρ、程序運行時間t，為避免偶然因素，結果為重復20次實驗求取平均值.具體對比結果如表2所示，顯然，ILMD前兩階分量的相關系數比EMD、CEEMDAN大，對軸承狀態信息表征更好，故障特征提取更精確.此外，CEEMDAN算法復雜，2.94 s的運行時間，遠超ILMD算法與EMD算法.綜合比較，ILMD算法在表征精確性與運行實時性兩方面，都優于EMD及CEEMDAN算法，因此，選擇ILMD對原始振動信號分解是合理且具有優勢的.

表2 各算法比較

3.2 特征提取

在信號多尺度形態學分析的過程中，結構元素G和最大形態學尺度εmax的選取非常關鍵.結構元素G的類型和長度是影響形態學覆蓋面積與分形維數計算的重要參數，將“零高度”的扁平型G作為單位結構元素，一方面保持了信號的原始特性，即其幅值免受結構元素高度的影響；另一方面，削減了部分計算量，增加算法實時性[18].因此，選取的單位結構元素為[0 0 0].而最大形態學尺度εmax的選取也會對結果產生重大影響，一般而言，εmax的標準為[1,N/2]之間的整數，通過文獻參考與多次調試，再權衡算法實時性與準確率之后，選擇εmax=8最為合適.

使用前文提到的ILMD算法分解七種工作狀態下的軸承振動信號，計算第一個PF分量的數學形態學的分形維數，并在一維空間中進行軸承狀態識別.當識別類型為正常狀態及0.007 in損傷直徑的三種故障時，比較結果如圖4所示.可見不同狀態的分形維數呈現出明顯的分形區間，故障診斷效果較好.當識別類型增加了0.021 in損傷直徑的三種故障后，比較結果如圖5所示，可以看到，正常狀態與滾動體故障有一定的分形區間，分形特性明顯，區分容易，但不同損傷程度滾動體故障的分形維數交叉重疊，難以區分；此外，對于0.007 in與0.021 in 的內、外圈故障，沖擊特性模糊，辨識不易，需要進一步的特征提取.綜上可知，利用分形特性對軸承工作狀態進行識別的方法是可行的，并且對不同軸承故障的識別效果良好，但對同種故障但不同損傷程度的狀態分辨效果差.因此，使用單一維度特征表征軸承狀態缺乏空間性與全面性，可以增加特征維數以更好的表征軸承狀態.

圖4 不同故障一維PF分量分形維數Fig.4 Fractal dimension of one-dimensional PF component of different faults

圖5 七種狀態下一維PF分量分形維數Fig.5 Fractal dimension of one-dimensional PF component in seven states

圖6為使用兩個PF分量分形維數的二維散點圖，由于向量維度的增加，盡管“特征混疊”現象有所改善，但0.021 in損傷直徑的內、外圈故障及不同損傷直徑的滾動體故障的特征仍然有部分混疊交聯，不易區分.圖7為三個PF分量的三維散點圖，相較于二維散點圖，“特征重疊”現象幾乎沒有改善.一方面是因為單一尺度特征向量的表征缺少完備性，另一方面則是因為第三階PF分量與原始信號的相關性系數較小，難以精準表征原始信號的分形特性.因此，需要采用其他特征替代第三階PF分量以改善識別效果.

圖6 七種狀態下二維PF分量分形維數Fig.6 Fractal dimension of two-dimensional PF component in seven states

圖7 七種狀態下三維PF分量分形維數Fig.7 Fractal dimension ofthree-dimensional PF component in seven states

作為描述分形集整體空間大小的物理量，“形態學覆蓋面積”也可以定值定量的表征軸承工作狀態，其理論依據已在前文給出.因此，將形態學覆蓋面積的對數作為第三維特征向量引入聚類識別中，其結果如圖8所示.相較于前文三種方法，該方法類聚合度高，“特征混疊”現象改善明顯，不僅在識別軸承不同工作狀態方面表現優異，在區分同種故障的不同損傷程度方面也表現良好.

圖8 引入形態學覆蓋面積的三維特征Fig.8 Three-dimensional features of the morphological coverage area introduced

上述實驗證明了本文提出的特征提取方法具有顯著優勢，并為軸承狀態監測與故障診斷提供了一種有效可行的新指標.此外，對不同損傷程度狀態的精確識別，也有助于軸承壽命的預測.

3.3 智能識別與分類

為減少人工干預，實現在線監測，引入概率神經網絡(PNN)對軸承狀態智能分類.實驗時，將不同狀態的樣本隨機劃分為兩份，第一份約占樣本數3/4，作為訓練樣本，另一份為測試樣本，因此在7種不同狀態下，訓練集共609組，測試集共210組.此外，在PNN網絡訓練中，較小的spread系數會優化分類結果[19]，經過多次調試，該參數設定為0.01.

確定參數后，用本文提出的三維特征向量構建輸入矩陣，并訓練生成神經網絡模型，以實現對軸承狀態的識別分類.為保證分類結果的準確性與科學性，隨機抽取樣本并重復實驗20次，統計其平均準確率.圖9為20次重復實驗中分類效果最優的一次，7種狀態下的識別準確率達到100%，不僅能夠識別不同狀態的軸承，并且成功對同種故障不同損傷程度的軸承狀態進行了分類.為進一步體現本文算法在實時性與準確率方面的提升，將其與不同算法比較，結果如表3所示.

圖9 七種狀態下的識別結果Fig.9 Recognition results in seven states

表3 不同算法各狀態的識別率與識別時間

綜上所述，針對滾動軸承狀態的在線監測與故障的智能診斷問題，本文提出算法在監測實時性與診斷準確度上面都有大幅改善，為滾動軸承的智能管理提供了一種有效的新方法.

4 結 論

本文提出了一種基于ILMD與形態學分形理論的特征提取方法，并結合PNN成功實現對軸承狀態的監測與故障診斷.該方法首先通過ILMD分解原始振動信號，抑制噪聲，提高特征維度；其次，使用數學形態學覆蓋求取前兩階PF分量的分形維數，并結合盒維數分形理論，引入形態學覆蓋面積構成三維特征矩陣；最后，通過CWRU中心軸承真實數據構建PNN神經網絡模型完成識別，并與其余算法比較，驗證了該方法的可行性與優異性.結果表明，本文所提方法不僅能識別不同狀態下的軸承，也可精確表征不同損傷程度的軸承.其平均識別時間只要0.21 s，平均識別準確率高達99.6 %，遠超其余算法，為滾動軸承狀態在線監測與智能故障診斷提供了一種高效、科學的方法.