基于ICSA-MCKD方法的滾動軸承聲信號微弱故障診斷

王 樹 杰， 李 宏 坤， 王 朝 閣， 孫 斌， 劉 艾 強

(大連理工大學 機械工程學院，遼寧 大連 116024)

0 引 言

滾動軸承作為旋轉機械的核心部件，通常承載著較大的荷載，當滾動軸承出現故障時，會造成機械傳動故障，進而產生嚴重的經濟損失.因此，對滾動軸承的診斷和監測可以提高機械系統的安全性和可靠性[1].目前，監測和診斷方法主要基于振動信號分析[2-4].然而振動傳感器采用接觸式方式安裝，限制了其使用，而聲音傳感器采用非接觸式，安裝方便.聲音信號包含豐富的信息，當機械設備發生故障時，聲音信號會產生相應的變化.

針對聲音信號低信噪比這一問題，提取出信號中的微弱故障信息是主要的研究方向，一些學者做了相關的研究.文獻[5]以峭度為指標提出了最小熵反卷積(minimum entropy deconvolution，MED)算法，但是該算法對隨機脈沖較為敏感，導致算法失效.為了充分利用軸承故障信號的周期性，文獻[6]以相關峭度為指標提出了最大相關峭度反卷積(maximum correlation kurtosis deconvolution，MCKD)算法，克服了MED算法易受單個脈沖影響的問題.Miao等[7]對MCKD算法進行了改進，通過尋找包絡信號的自相關最大值來估計迭代周期，克服了輸入參數嚴格和復雜重采樣的難題.申博文等[8]提出了基于最大相關峭度反卷積與自適應噪聲完備集合經驗模態分解的聲信號故障特征提取方法，先用MCKD算法增強聲信號中的沖擊，然后再計算每個經驗模態分量的峭度，選出最優分量，提取故障信息.Chen等[9]針對最大相關峭度反卷積、多點最優最小熵反卷積調整、最大二階循環平穩盲反卷積等方法高度依賴測量信號先驗周期的問題，引入了6種周期性檢測技術來自適應地識別重復脈沖的周期.張俊等[10]將變分模態分解(variational mode decomposition，VMD)和MCKD相結合，并采用粒子群算法(particle swarm optimization，PSO)對VMD和MCKD中的參數組合進行尋優.劉尚坤等[11]將Teager能量算子和MCKD用于滾動軸承的故障識別，首先采用MCKD對信號進行降噪，然后使用Teager算子增強信號中的周期性沖擊.Lyu等[12]用量子遺傳算法(quantum genetic algorithm，QGA)自適應地選擇MCKD的濾波器長度和反卷積周期，用于齒輪和軸承復合故障診斷.Wang等[13]提出了一種基于并行雙參數優化共振稀疏信號分解(resonance-based sparse signal decomposition，RSSD)和優化的多點最優最小熵反卷積調整(multipoint optimal minimum entropy deconvolution adjusted，MOMEDA)的復合故障特征提取方法.

MCKD算法可以有效提取出滾動軸承聲信號中的故障周期信息，然而其依賴于濾波器長度和先驗故障周期的選擇.布谷鳥搜索算法(cuckoo search algorithm，CSA)是一種生物啟發式算法，該算法模擬了布谷鳥將自己的卵寄生在其他鳥類巢穴的行為.CSA具有優秀的全局尋優能力和局部尋優能力，但是該算法同其他生物啟發式算法一樣，也存在后期尋優精度低、迭代速度慢、易陷入局部最優的問題.Tsipianitis等[14]在算法中引入靜態和動態懲罰函數，以增強CSA，并且結合鳥群算法(bird swarm algorithm，BSA)的關鍵參數，提出了具有動態懲罰的CS混合方法.Mareli等[15]提出了3種隨迭代次數動態增加的發現概率函數，使全局尋優和局部尋優更加平衡，通過測試函數驗證了該方法優于發現概率固定的布谷鳥搜索算法.李榮雨等[16]通過調整萊維飛行步長和在偏好隨機游走中引入動態慣性權重及記憶策略，提高了算法的穩定性和搜索能力.針對CSA固定步長和發現概率問題，本文提出改進的布谷鳥搜索算法(improved cuckoo search algorithm，ICSA)，在萊維飛行中引入步長尺度因子，在偏好隨機游走中將固定發現概率調整為隨迭代次數增加自適應調整的變發現概率.故本文將ICSA算法與MCKD算法相結合，用ICSA搜索出MCKD算法的濾波器長度和故障周期最佳組合，然后對反卷積后的信號進行包絡譜分析，進行滾動軸承的故障診斷.

1 基本原理

1.1 最大相關峭度反卷積

實測信號xn(n=1，2，…，N)是源振動信號sn(n=1，2，…，N)與傳遞路徑h卷積得到的，MCKD是通過尋找一個FIR濾波器，進行反卷積運算得到源信號sn的近似解yn(n=1，2，…，N)：

(1)

式中：f(f1，f2，…，fL)為濾波器，L為濾波器長度，*為卷積運算.

通過最大化相關峭度求解濾波器f，相關峭度的定義為

(2)

式中：T為故障周期，T=fs/fi，fs為采樣頻率，fi為故障特征頻率；M為位移數，一般取M=7.

最大化相關峭度求解濾波器f(f1，f2，…，fL)：

(3)

相關峭度對濾波器求導：

(4)

由式(1)～(4)可求得濾波器的系數，并表示成矩陣形式：

(5)

式中

綜上所述，MCKD算法的流程如下：

(1)初始化故障周期T、濾波器長度L、位移數M等參數；

(4)根據yn計算αm、β；

(5)更新濾波器f；

(6)判斷ΔKc,M(T)是否小于閾值，若小于閾值，結束迭代，否則重復步驟(3)～(5).

1.2 CSA

CSA由控制全局搜索的萊維飛行和局部搜索的偏好隨機游走組成.在CSA中有3條理想化規則：

(1)每只布谷鳥一次產一個卵，并隨機選擇一個寄生鳥巢.

(2)最高質量的卵可以存活到下一代.

(3)寄生鳥巢數量固定，且寄生卵被發現的概率是pa.當宿主發現寄生卵后，會選擇拋出寄生卵或者放棄鳥巢重新筑巢.

每個寄生鳥巢代表一組解，通過萊維飛行更新鳥巢位置：

(6)

(7)

式中：μ和ν服從標準正態分布，β=1.5.

(8)

萊維飛行的步長

(9)

式中：α0為步長因子，通常α0=0.01；Xb表示當前最優解.

綜合式(6)～(8)通過萊維飛行到達的新位置為

(10)

在萊維飛行后，產生一個隨機數與發現概率pa比較，一般取pa=0.25，當隨機數大于pa，代表寄生卵被發現，進行偏好隨機游走，產生新的鳥巢：

(11)

1.3 ICSA

在標準CSA中，萊維飛行的步長因子α0是固定的，但是步長因子對給定的優化問題會很敏感，若α0一直較大，則算法的全局搜索能力很強，但是會降低搜索精度，若α0一直較小，會導致算法收斂太慢.α0可以是固定的或是變化的[17]，故本文引入步長尺度因子，提出變步長的CSA(ICSA)，將式(9)改為式(12)：

(12)

式中：F為步長尺度因子，服從[0，1]上的均勻分布.

式(10)更新為式(13)：

(13)

在標準CSA中，偏好隨機游走的發現概率pa為固定值0.25，pa是局部最優解搜索的關鍵，pa太大會導致搜索精度不夠，pa太小會導致錯過最優解，使得算法收斂過慢.故本文提出了隨迭代次數更新的發現概率，在算法開始時，pa較小，增加算法的搜索精度，在算法后期，有較大的pa，增加算法的收斂速度.調整固定pa為式(14)：

(14)

1.4 ICSA-MCKD方法

在用ICSA優化參數時，需要確定一個適應度函數.文獻[18]提出了諧波顯著性指標(harmonic significant index，HSI)來識別故障信息，其表達式為

(15)

式中：F(ω)為故障頻率ω處的幅值，K為要計算的諧波數量，N(ω)為ω附近的噪聲量，P(kω)=F(kω)/N(kω)表示故障頻率ω的顯著程度.但是，該指標對于實際信號不能完全適用：

(1)計算ω附近的噪聲量N(ω)時，均采用ω左右各5個頻率點對應的幅值進行平均，沒有考慮內圈故障時轉頻的調制影響，不能反映出故障頻率周圍的整體噪聲.

(2)H(ω)為多階諧波指標相乘，受單階諧波指標影響較大，當只有某一階的諧波指標較大，而其他階諧波指標較小時，也會使H(ω)值較大，導致誤判.

基于以上問題，提出調整的諧波顯著性指標(adjusted harmonic significant index，Ha)，將其作為ICSA-MCKD的適應度函數，其表達式為

(16)

式中：F(kω)為反卷積信號的包絡譜在kω處的幅值，N(kω)為kω左右轉頻范圍內的平均噪聲，a(·)為取均值.調整的諧波顯著性指標計算kω附近以轉頻為間隔的邊頻帶內的噪聲量而不是左右各5個頻點，考慮了轉頻的調制作用；計算多階諧波指標的平均值，克服了只有單階諧波指標較大導致HSI誤判的問題.

綜上所述，ICSA-MCKD故障診斷方法的流程如圖1所示，主要步驟如下：

(2)計算每個鳥巢的目標值，找到當前位置目標值最大的鳥巢，作為最優鳥巢Xb.

(3)進行萊維飛行，更新位置，最優鳥巢保持不動.計算新位置下每個鳥巢的目標值Ha(ni)，并與舊鳥巢的目標值比較，如果新鳥巢目標值大于舊鳥巢目標值，則替換舊鳥巢.

(4)進行偏好隨機游走，計算新鳥巢的目標值，如果大于舊鳥巢的目標值，則替換舊鳥巢.

(5)選擇最優解，判斷是否滿足迭代停止條件，若滿足，則根據最優參數組合進行MCKD濾波，進而進行包絡譜分析，若不滿足則返回步驟(3)繼續循環.

2 仿真研究

為驗證本文提出的方法對滾動軸承故障診斷的有效性，構造了內圈故障軸承的仿真信號：

x(t)=e(t)+r(t)+p(t)+n(t)

其中

δi)+φa]

徐進步一拍腦袋:“啊,想起來了,‘世上本無路’那一句,難怪聽著有印象。可就他,八成沒讀過魯迅的什么書吧?”

式中：e(t)代表故障脈沖信號，W1=128為故障脈沖個數，A(t)=0.5(1-cos(2πfrt))為故障脈沖幅值，fr=10 Hz為轉頻，Ta=1/128 s為脈沖間隔時間，即故障特征頻率為fi=128 Hz，fa=2 600 Hz為固有頻率，ξa、δi、φa分別為衰減系數、滾子滑移時間延遲和初始相位.r(t)為隨機脈沖信號，W2=3為脈沖個數，振幅Bs和出現時間Ts分別服從正態分布和均勻分布，fb=1 700 Hz.p(t)代表諧波信號，W3=2為諧波個數，且諧波頻率f1=10 Hz，f2=20 Hz.n(t)為高斯白噪聲.加入噪聲后信噪比為-11.6 dB，采樣頻率為12 800 Hz，采樣時間為1 s.

圖2(a)為各分量及加噪后的仿真信號時域圖，故障周期信號已經被噪聲淹沒；圖2(b)為仿真信號包絡譜，由于噪聲干擾，圖中難以看出故障特征頻率.

現將本文的方法用于仿真信號處理.利用ICSA優化MCKD參數組合，初始化ICSA的參數，鳥巢個數n=15，最大迭代次數tmax=50，目標函數為Ha(ni).濾波器長度L的搜索區間為[100，1 000]，故障周期T的搜索區間為[98，102].

為驗證ICSA的優越性，用CSA和PSO對同一仿真信號進行參數尋優，圖3(a)為各參數優化算法的迭代圖，為了公平起見，CSA的初始化參數和ICSA的保持一致.PSO的種群個數和最大迭代次數也與ICSA的保持一致，分別為n=15，tmax=50，學習因子c1=c2=1，慣性權重w=1.濾波器長度L和故障周期T的搜索區間都分別設置為[100，1 000]和[98，102].從圖中可以看出，ICSA在迭代30次時達到最優解，Ha為3.434；CSA在迭代37次時達到最優解，Ha為3.205；PSO在迭代33次時達到最優解，Ha為2.832.Ha作為反卷積效果的一個標準，Ha越大，說明反卷積效果越好，故障特征越明顯.由此可見，ICSA在迭代次數最少時有最大的適應度函數值，故ICSA有更快的迭代速度和搜索精度.ICSA的尋優參數組合為L=217，T=100，該參數下MCKD反卷積信號的包絡譜如圖3(b)所示，在圖中可以明顯看到故障特征頻率128 Hz及其倍頻，并且可以看到轉頻調制的邊頻帶，證明了本文方法的有效性.

為了證明ICSA對MCKD參數尋優結果的可靠性，現改變最佳參數組合中的一個，并以更改后的參數組合輸入MCKD中對仿真信號進行處理.圖4(a)為將最佳參數組合[217，100]中的濾波器長度L改為215后的MCKD反卷積信號包絡譜.與圖3(b)相比圖4(a)中內圈故障特征頻率的四倍頻、五倍頻被噪聲干擾嚴重，看不出明顯的邊頻帶，且內圈故障特征頻率及其倍頻的幅值都有所下降.圖4(b)為將最佳參數組合[217，100]中的故障周期T改為101后MCKD反卷積信號的包絡譜.與圖3(b)相比圖4(b)中只能隱約看出內圈故障特征頻率及其二倍頻和三倍頻，且噪聲干擾嚴重，而高階倍頻都被噪聲淹沒.通過圖4(a)與圖4(b)的對比分析可知，MCKD對參數的選擇有很大的依賴性，且對故障周期T的選擇更為敏感.而人為主觀地選擇參數具有偶然性，導致反卷積效果不理想，而通過ICSA自適應地選擇參數，可以達到最佳的分析效果.

為了凸顯MCKD的優越性，采用MED對上文的模擬信號進行處理，MED的濾波器長度和MCKD的一致(L=217)，經MED濾波處理后信號的包絡譜如圖5所示，沒有明顯的故障特征頻率及其倍頻，說明MED對有隨機沖擊和強噪聲信號沒有很好的作用，對比圖3(b)和圖5可知，MCKD有更好的故障信息提取能力.

3 試驗研究

為了進一步驗證該方法的有效性，現對滾動軸承內圈故障的實測信號進行分析.試驗在QPZZ-Ⅱ試驗平臺上進行，如圖6(a)所示.測點1選在遠離軸承座的位置，測點方向為正對軸承座，以模擬聲音信號的微弱故障，另選測點2靠近軸承座的位置做對比，測點方向也為正對軸承座，如圖6(b)所示.軸承內圈故障如圖6(c)所示，加工方式為線切割，所用軸承型號為NU205EM/PS，具體參數如表1所示.使用PCB噪聲麥克風和NI9234采集卡采集聲音信號，設置采樣頻率為12 800 Hz，采樣時間為2 s，轉速為900 r/min，根據理論公式計算fi=116.25 Hz，T=110.

現分別對測點1和測點2處信號進行包絡譜分析，圖7(a)為測點2處信號的包絡譜，在圖中可以明顯看到內圈故障特征頻率及其二倍、三倍故障特征頻率，以及轉頻調制的邊頻帶.雖然，通過測點2處信號的包絡譜可以判斷出故障特征，但是故障特征頻率被轉頻及其倍頻以及噪聲所干擾，對故障診斷造成困難.圖7(b)為測點1處信號包絡譜，在圖中只能看到明顯的轉頻以及二倍轉頻，而故障頻率被噪聲淹沒，無法判斷故障信息.由此可得出結論：聲音信號對距離較為敏感，考慮到實際工程應用中，傳感器距被測對象都有一定的距離，故需要對實測信號進行進一步的處理.

現用本文提出的方法對測點1處的信號進行處理.初始化ICSA的參數，鳥巢個數n=15，最大迭代次數tmax=50，目標函數為Ha(ni).濾波器長度L的搜索區間為[100，1 000]，故障周期的搜索區間為[108，113].

為驗證ICSA對實際信號也有很好的搜索性能，與CSA、PSO對同一試驗信號進行參數尋優.圖8為各參數優化算法的迭代圖，CSA和PSO的初始化參數同仿真信號保持一致，濾波器長度L和故障周期T的搜索區間都分別為[100，1 000]和[108，113].從圖中可以看出，ICSA在迭代33次時達到最優解，Ha為6.596；CSA在迭代39次時達到最優解，Ha為6.482；PSO在迭代47次時達到最優解，Ha為6.41.由此可見，ICSA在迭代次數最少時有最大的適應度函數值，故ICSA 有更高的迭代速度和搜索精度.

ICSA的最佳尋優參數為L=622，T=111.5，該參數下MCKD反卷積信號的包絡譜如圖9(a)所示，在圖中可以明顯看到故障特征頻率115.5 Hz及其倍頻.CSA的最佳尋優參數為L=654，T=111，該參數下MCKD反卷積信號的包絡譜如圖9(b)所示，相比于圖9(a)，整體譜線幅值更小，噪聲干擾更大，證明了ICSA的尋優優越性.

為驗證故障周期的選擇會影響MCKD對實際信號的反卷積效果，現把最佳參數組合[622，111.5]中的故障周期換為理論故障周期110，進行MCKD反卷積，圖10為反卷積后的包絡譜.與圖9(a)相比圖10中故障頻率的高階倍頻被噪聲干擾嚴重，且整體譜線的幅值都有所下降.說明在實際工程問題中，轉速波動以及滾動體滑移會導致實際故障特征頻率與理論計算值有一定的誤差，用理論值作為MCKD的先驗參數會導致故障特征不明顯，對故障診斷造成困難.

為了凸顯MCKD對實際信號處理的優越性，對同一實際信號用MED進行濾波處理，MED的濾波器長度和ICSA-MCKD的一致(L=622)，對MED濾波后的信號進行包絡譜分析，如圖11所示，圖中雖能找到故障特征頻率，但是噪聲干擾嚴重，對比圖9(a)和圖11可知，MCKD有更好的故障特征提取能力.

4 結 語

本文提出了一種基于改進的布谷鳥搜索算法(ICSA)和最大相關峭度反卷積(MCKD)相結合的滾動軸承故障診斷方法，通過仿真信號和實際信號的比較，證明了該方法的有效性和優越性.MCKD可以在強噪聲背景下提取出滾動軸承的故障信息，然而從實測信號分析可知，由于轉速波動和滾動體滑移等因素的影響，滾動軸承的實際故障頻率和理論計算值存在偏差，這就會影響MCKD的效果，故本文提出使用參數優化算法搜索MCKD的最佳參數組合.CSA在算法后期存在尋優精度不足和收斂速度較慢的問題，本文提出了ICSA，在萊維飛行中引入了步長尺度因子和自適應發現概率，提高了算法的搜索精度和迭代速度，并與CSA、PSO算法對比，驗證了ICSA的優越性.