參數優化最大二階循環平穩盲解卷積行星輪軸承故障提取

林 云， 郭 瑜， 陳 鑫

(昆明理工大學 機電工程學院，昆明 650500)

行星齒輪箱廣泛用于風力發電、直升機、工程機械等大型復雜機械裝備中，其惡劣的工作環境導致其易產生故障。行星輪軸承是行星齒輪箱中較容易產生故障的零部件，因此研究行星輪軸承故障特征提取方法對機械設備的安全運行具有重要意義。

滾動軸承故障源信號的傳遞過程可以看作是源信號與信道的一個線性卷積混合過程[1]，因此可以通過解卷積運算恢復故障沖擊脈沖。相關解卷積算法中，Mcdonald等[2]最近提出的多點優最小熵反卷積(multipoint optimal minimum entropy deconvolution adjusted, MOMEDA)算法，在滾動軸承故障特征提取上取得較好效果[3]，但Buzzoni等[4]使用的周期準則與軸承故障沖擊的循環平穩性不符，在強噪聲干擾下，算法可能會失效。近期，Buzzoni等提出最大二階循環平穩盲解卷積(maximum second order cyclostationary blind deconvolution，CYCBD)算法。算法利用了二階循環平穩指標(second order indicators of cyclostationary，ICS2)對故障沖擊敏感的特性，具有較好的魯棒性。近年來，針對滾動軸承故障診斷的研究表明，該算法可以有效提取軸承故障沖擊，且相對于Momeda等的算法更具優勢[5]。但在軸承滑移條件下，利用CYCBD技術提取行星輪軸承故障特征，算法中循環頻率和濾波器長度的準確選取仍是一個難點。

近期，羅忠等[6]提出了故障特征比(fault feature ratio，FFR)指標，其旨在反映包絡譜中的故障特征頻率及其諧波的顯著程度，但尚存在如下問題：①未考慮軸承滑移的情況；②指標中諧波均值較大時，不一定反映了各個諧波均顯著；③指標所選的頻率范圍受到干擾較大。為此，本文研究中提出一種改進故障特征比(improved fault feature ratio，IFFR)指標，以該指標作為粒子群算法的適應度函數自動獲取CYCBD算法中循環頻率和濾波器長度。通過行星輪軸承外圈故障特征提取的仿真和試驗對比分析，結果反映了本文方法的有效性和優勢。

1 基于改進FFR指標的適應度函數

在利用智能算法優化參數時，適應度函數決定了整個優化過程的方向和目標，而建立合適的適應度函數是此類問題的重點和難點。近期，羅忠等提出了FFR指標作為粒子群算法的適應度函數。該指標通過計算包絡譜中故障特征頻率分量所占比重來反映故障特征提取效果，表達式如下

(1)

ASi=max[AS(ifm-0.02fm,ifm+0.02fm)]

(2)

式中:ASi為包絡譜中第i階故障特征頻率處的幅值；I為選取的故障諧波數；ASt為包絡譜中1～1 000 Hz頻率幅值之和；E(·)為求均值；fm為故障特征頻率。

盡管該指標能夠一定程度反映故障特征提取效果，但仍存在以下不足：①在式(1)分子中，指標通過故障各階諧波的均值來反映故障特征顯著程度，而當包絡譜中某一階故障特征頻率明顯突出時，即使其他諧波處幅值即使并不顯著，也會使該均值較大，此時FFR指標將主要反映單個諧波的顯著程度；②式(2)中對每一階諧波附近進行峰值搜索的范圍均為理論故障頻率的±2%內，當軸承在運行過程中存在一定滑移時，在較高階的諧波處，實際故障特征頻率將可能超出該峰值搜索范圍，導致指標失效；③平方包絡譜開始的前幾赫茲頻率有可能出現一些大峰值的影響[7]，因此，固定地棄掉前1 Hz頻率成分難以有效消除平方包絡譜具有較大幅值的低頻干擾。限制了該指標在平方包絡譜中的應用。

為解決上述問題，研究中提出IFFR指標，其表達式為

IIFFR(i)=min[Mi/E(MSES)],i=1,…,N

(3)

Mi=max{M[ifω(1-0.02),ifω(1+0.02)]}

(4)

MSES=[M(fω-0.5fω,Nfω+0.5fω)]

(5)

式中:M(·)為取信號幅值；i為包絡譜中第i階故障諧波；N為所取的故障諧波總數；E(·)為求均值；fω為故障特征頻率。

IFFR指標選取了信號前i階故障諧波中故障特征比的最小值，從而保證了當該指標越大，前i階故障諧波越顯著；第i階理論故障諧波處的峰值搜索范圍相應擴大為第一階理論故障諧波的i倍，保證了滑移條件下，實際故障特征頻率及其諧波不會超出搜索范圍；該指標在求包絡譜均值時選取式(5)中較窄的頻率范圍，避免了其他頻率范圍的影響。

行星輪軸承外圈故障仿真和試驗信號處理結果表明，該指標能夠較好地揭示包絡譜中各階諧波的顯著程度，進而保障故障特征提取的準確性和效果。

2 參數優化CYCBD算法

2.1 CYCBD算法原理

CYCBD算法是一種基于廣義瑞利熵和循環平穩指標的新型盲解卷積算法，可以用于提取具有循環平穩特性的故障沖擊信號。在解卷積算法模型中，一般通過測量信號x來估計目標輸入信號s0，此類算法的基本原理可以用式(6)表示

s=x*h≈s0

(6)

式中:s為估計源信號；h為逆濾波器；*為卷積運算。式(6)的矩陣形式為

s=Xh

(7)

(8)

式中:L為離散信號s的點數；N為逆濾波器h的長度。

在廣義瑞利熵的框架下，求解最大ICS2值被轉化為求解廣義瑞利熵的特征值，此時ICS2可定義為

(9)

式中，RXWX和RXX分別為加權相關矩陣和相關矩陣。加權矩陣W可以表示為

(10)

E=(e1…ek…ek)

(12)

(13)

式中:k為循環頻率集的大小；Ts為故障沖擊周期。由于循環頻率被認為是與信號能量波動相關的頻率，可與軸承故障或齒輪故障等相關，因此可將離散信號的循環頻率定義為k/Ts[8]，而離散信號的理論故障頻率為1/Ts，因此可基于理論故障頻率計算循環頻率。

需要說明的是，CYCBD算法目的在于求解最佳逆濾波器，而獲取加權矩陣W是其核心的步驟。在計算W時，需要提供濾波器長度和循環頻率這兩個參數，其設置準確與否對最終求得的逆濾波器會產生較大影響。

2.2 參數優化CYCBD

粒子群算法[9]是一種仿生優化模型，該算法具備較好的全局優化性能。其原理簡述如下：①在D維空間中，N個粒子組成一個種群，則第m個粒子為一個D維向量，粒子的位置和速度分別為Xm=(xm1,xm2,…,xmD)，Vm=(vm1,vm2,…,vmD)，其中，m=1,2,…,N；②第m個粒子搜索到的最優位置和整個粒子群搜索到的最優位置可分別記為Pbest=(pm1,pm2, …,pmD),gbest=(g1,g2, …,gD),其中，m=1,2,…,N;③各個粒子可分別根據式(14)和式(15)來更新自己的速度和位置

vmn(t+1)=λ·vmn(t)+c1r1(t)[pmn(t)-xmn(t)]+

c2r2(t)[pgn(t)-xmn(t)]

(14)

xmn(t+1)=xmn(t)+vmn(t+1)

(15)

式中:c1,c2為學習因子,分別為個體經驗和社會經驗對粒子運動路線的影響；r1,r2為0～1的隨機數；λ為慣性權重, 慣性權重的大小表示粒子下一步速度對當前速度繼承量的百分比。

通過合理設置各項參數，可以使粒子快速而準確到達目標位置，完成參數優化過程。在式(14)中學習因子c1和c2按照文獻[10]設置為c1=c2=1.5；慣性權重設置在0.9～1.2時，算法具有比較好的搜索性能[11]，本文設置為λmax=1.2，λmin=0.9。

3 故障特征提取流程

本文提出一種基于參數優化最大二階循環平穩盲解卷積的行星輪軸承故障特征提取方法，其實現基本步驟為：

步驟1利用鍵相脈沖信號將原始時域信號通過階比分析技術轉換為角域信號，消除轉速波動干擾；

步驟2利用離散隨機分離(discrete random separation，DRS)[12]技術將齒輪相關信號和軸承信號進行分離，抑制齒輪相關信號的影響；

步驟3利用CYCBD技術對軸承信號部分進一步處理，在CYCBD算法中，基于理論故障特征階次確定循環頻率搜索范圍，基于理論故障周期點數確定濾波器長度搜索范圍，粒子的位置坐標由循環頻率和濾波器長度兩個參數構成，根據參數搜索范圍和搜索精度，設定種群大小N和最大迭代次數I；

步驟4設置適應度函數為IFFR，初始化CYCBD參數，其中迭代步數設置為5，迭代精度為10-4；

步驟5計算種群中N個個體的適應度值，記錄最大的適應度值，記錄此時粒子的位置；

步驟6更新種群位置，求種群N個個體的適應度值，記錄此輪迭代最大適應度值，與當前最大適應度值進行對比，取其大者，更新最大適應度值；

步驟7重復步驟5和步驟6，直到達到最大迭代次數I，輸出最大適應度值及其對應的位置坐標；

步驟8根據位置坐標得到的循環頻率和濾波器長度作為CYCBD的優化輸入參數，算法中循環頻率集大小取為k=100。求解卷積信號的平方包絡譜，提取故障特征。

該方法的流程如圖1所示。

圖1 基于參數優化最大二階循環平穩盲解卷積的行星輪軸承故障提取方法Fig.1 Fault extraction method of planet bearing based on parameter optimized CYCBD

4 仿真驗證

4.1 仿真信號構造

滾動軸承的局部故障振動仿真信號，可用式(16)～式(18)表示[13]

(16)

Ai=Aocos(2πQt+φA)+CA

(17)

s(t)=e-βtsin(2πfnt+φw)

(18)

式中:Ai為信號幅值調制部分；M為故障沖擊脈沖個數；T為故障沖擊時間間隔；τi為滾動軸承隨機滑移；fn為系統的固有頻率；Q為軸承轉頻；φw和φA為初始相位；n(t)為高斯白噪聲。

仿真的行星輪軸承外圈局部故障振動信號主要包括：齒輪信號、軸承信號和噪聲信號三部分。

(1)軸承信號。與單一滾動軸承故障不同的是，行星輪軸承發生局部故障時，軸承信號具有更為復雜的調制。由于行星傳動結構，行星輪軸承外圈既有自轉又有公轉。其中，自轉與行星輪繞行星輪軸的旋轉有關，二者轉頻相同，在故障振動信號中表現為行星輪的轉頻調制現象。整個行星輪軸承會隨行星架一起旋轉，產生時變路徑效應，在故障振動信號中表現為與行星架轉頻相關的調幅作用。對于傳感器安裝于齒圈上方的情況，時變路徑調幅效應可以用Hanning窗函數來表示[14]。基于上述分析，行星輪軸承外圈局部故障仿真，軸承信號部分可用式(19)和式(20)表示

[1-cos(2πfct)]×s(t-iTP-τi)(19)

s(t)=e-βtsin(2πfn1t+φw)

(20)

式中:Ai為第i個故障沖擊的幅值；fp為行星輪的旋轉頻率；fc為行星架旋轉頻率；τi為軸承滑移，在0～0.02TP隨機取值；fn1為軸承固有頻率；β為阻尼比。

(2)齒輪信號。行星齒輪箱在運轉過程中，齒輪嚙合包括行星輪分別和太陽輪、齒圈的嚙合。齒輪輪齒嚙合同樣會激起共振，且輪齒嚙合引起的共振大于故障軸承[15]。齒輪信號部分可表示為

(21)

式中:Aj為第j次輪齒嚙合的幅值；N為齒輪嚙合次數；TQ為齒輪信號的嚙合周期；fn2為齒輪固有頻率。

(3)噪聲信號。考慮噪聲影響，在仿真信號中加入信噪比為-10 dB的高斯白噪聲。

行星輪軸承外圈故障仿真信號可用式(22)表示

x(t)=xout,bearing(t)+xgear(t)+xnoise(t)

(22)

式中，xnoise(t)為高斯白噪聲。

仿真信號的具體參數如表1所示。

表1 行星輪軸承外圈故障仿真參數Tab.1 Simulation parameters of planetary gear bearing with fault on outer race

行星輪軸承故障沖擊周期可由式(23)計算

(23)

式中，fbo為行星輪軸承外圈故障特征頻率。

4.2 故障特征提取

仿真信號原始時域信號如圖2所示，其頻譜如圖3所示，圖中輪齒嚙合頻率及其諧波較為顯著，而行星輪軸承外圈故障特征頻率fbo(100 Hz)被噪聲淹沒。

圖2 時域信號Fig.2 Time domain signal

圖3 信號頻譜Fig.3 Spectrum of signal

原始信號中，齒輪信號等強周期信號的干擾較大，難以直接提取行星輪軸承故障沖擊特征。因此，先用DRS技術對原始時域信號進行信號分離，得到周期信號部分和隨機信號部分，其中隨機信號部分包含了軸承故障沖擊信號，用該部分作為CYCBD算法的輸入信號。

由表1可知，行星輪軸承外圈理論故障特征頻率為fbo=100 Hz，考慮軸承存在約2%以內的滑移[16]，粒子群算法中故障循環頻率的搜索區間設置為：[98,102]；而最小濾波器長度要能覆蓋軸承故障所產生的沖擊信號衰減周期且最大長度不宜過大，本文設置濾波器長度的搜索范圍為0.5Tbo～1.5Tbo，即[120,360]。此外，循環頻率搜索區間長度(按循環頻率的搜索精度為0.1 Hz計算)為40，按照黃包裕等設置種群大小為區間長度的1/2，設置迭代次數I=10。

首先利用基于FFR指標的粒子群算法優化CYCBD參數，指標中諧波總數N取為4。在粒子群算法中，循環頻率搜索速度范圍設置為[-0.02,0.02]；濾波器長度搜索速度范圍設置為[-10,10]。算法的迭代圖如圖4(a)所示，圖中可見第6次迭代時，FFR指標值達到最大。對應的優化循環頻率和濾波器長度分別為99.6 Hz和296。將此參數應用于CYCBD算法，解卷積結果如圖4(b)所示。從圖4(b)中可見，相對于圖2所示的原始故障信號，故障沖擊一定程度被增強。結果的平方包絡階比譜如圖4(c)所示，該指標指示的平方包絡譜中，故障特征頻率的二階諧波較為顯著，而故障特征頻率處被噪聲淹沒，說明FFR指標存在主要反映單個諧波顯著程度的情況，不利于故障特征的有效提取。

圖4 基于FFR指標的CYCBD算法處理仿真信號Fig.4 FFR index based CYCBD algorithm to process simulation signal

隨后，采用本文所提出的基于IFFR指標的粒子群算法優化CYCBD參數。圖5(a)所示為利用基于IFFR指標的粒子群算法優化CYCBD參數過程的迭代圖，圖中可見在第7輪迭代時，優化值達到最大，此時對應的優化循環頻率和濾波器長度分別為99.1 Hz和342。將此參數應用于CYCBD算法，解卷積結果如圖5(b)所示。相對于圖4(b)所示的解卷積結果，故障沖擊增強效果更好。結果的平方包絡譜如圖5(c)所示。圖中故障特征階次及其諧波較為明顯的反映出了行星輪軸承外圈故障特征，表明利用IFFR指標能夠較好地反映平方包絡譜中的故障顯著程度。

圖5 基于IFFR指標的CYCBD算法處理仿真信號Fig.5 IFFR index based CYCBD algorithm to process simulation signal

相對于圖4(c)，圖5(c)能夠更好地反映行星輪軸承外圈故障特征，而計算圖5(c)所示結果的FFR值僅為0.31，小于圖4(c)結果的FFR值(2.33)。可知在獲取CYCBD算法優化參數組合時，利用IFFR指標作為適應度函數，其效果要優于相同條件下的FFR指標。

5 試驗驗證

5.1 試驗數據采集

利用行星齒輪箱綜合試驗臺進行試驗，以NGW型單級行星齒輪箱為試驗對象，采集行星輪軸承故障振動數據，輸入轉速設置為1 000 r/min。對本文所提方法進行驗證。加速度傳感器和電渦流傳感器的安裝位置如圖6所示，分別用于拾取振動信號和同步轉速脈沖。其中加速度傳感器為RION PV-86 4527型；電渦流傳感器為DH90型。通過齒圈傳遞到箱體正上方的信號傳遞路徑最短，可以得到信噪比相對較高的信號，試驗中在該處拾取振動信號。此外，信號采集系統中使用NI USB9234型采集卡和DH-5853型電荷放大器，采樣頻率為51.2 kHz，信號放大倍數為3。

圖6 傳感器安裝位置Fig.6 Location of sensor installation

齒輪箱內行星輪軸承、齒輪參數分別如表2、表3 所示。由相關故障理論，利用表2和表3數據計算得到行星輪軸承外圈故障相關頻率如表4所示。

表2 行星輪軸承參數Tab.2 Parameters of planet bearing

表3 齒輪參數Tab.3 Parameters of gears

表4 行星輪軸承外圈故障相關頻率Tab.4 Frequencies related to the planet bearing with faulty outer race 單位：Hz

研究中，選擇輸出軸(對應行星架轉軸)作為參考軸，即其階次為1，其他相關元件可利用對應信號頻率值按式(24)計算階次

(24)

式中:l為階次；f為信號頻率；n為參考軸轉速。

由式(24)可計算得到行星輪軸承外圈故障相關階次如表5所示。

表5 行星輪軸承外圈故障相關階次Tab.5 Orders related to the planet bearing with faulty outer race

圖7所示為軸承故障為電火花加工的人造局部故障，故障寬度約為1 mm，深度約為0.5 mm。

圖7 行星輪軸承外圈人造局部故障Fig.7 Artificial local fault of outer race of planet bearing

5.2 故障特征提取

為消除轉速波動的影響，以行星架轉軸為參考，將原始時域振動信號轉換到角域。設置等角度重采樣率設置為Fs=9 088 r/min，為齒圈齒數的128倍。原始角域信號如圖8所示，其包絡階比譜如圖9所示。從圖9中可見，各部件轉頻成分及其諧波較為顯著，而行星輪軸承外圈故障特征階次lbo(9.57)幾乎被噪聲淹沒。

圖8 角域信號Fig.8 Angular domain signal

圖9 角域信號包絡階比譜Fig.9 Envelope order spectrum of angular domain signal

利用DRS技術對原始角域信號進行信號分離，分離得到的隨機信號部分進行后續的參數自適應CYCBD算法處理。考慮到算法計算效率，對DRS分離后的隨機信號部分進行降采樣，重采樣頻率Fs=4 544點/轉。由表5可知，行星輪軸承外圈理論故障階次為lbo=9.57，根據上文中關于理論故障頻率與循環頻率的關系，對應循環頻率的搜索區間設置為：[9.37,9.76], 區間長度(按循環頻率的搜索精度為0.01階計算)為40；由式(23)計算得，Tbo=474.8點，設置濾波器長度的搜索范圍為0.5Tbo～1.5Tbo，即[237,713]。設置種群大小為區間長度的1/2，設置迭代次數I=10。

首先利用基于FFR指標的粒子群算法優化CYCBD參數，指標中諧波總數N取為4。在粒子群算法中，循環頻率搜索速度范圍設置為[-0.02,0.02]；濾波器長度搜索速度范圍設置為[-10,10]。算法的迭代圖如圖10(a)所示，圖中可見第7次迭代時，FFR指標值達到最大。對應的優化循環頻率(階次)和濾波器長度分別為9.66和395。將此參數應用于CYCBD算法，解卷積結果如圖10(b)所示，對比圖8原始角域信號，圖10(b)中難以觀察到故障特征增強效果。結果的平方包絡階比譜如圖10(c)所示，圖中虛線指示為故障特征階次及其諧波位置，該指標指示的平方包絡譜中，難以觀察到故障特征。

圖10 基于FFR指標的CYCBD算法處理試驗信號Fig.10 FFR index based CYCBD algorithm to process experimental signal

隨后，采用本文所提出的方法。圖11(a)所示為利用基于IFFR指標的粒子群算法優化CYCBD參數過程的迭代圖，圖中可見在第6輪迭代時，優化值達到最大，此時對應的優化循環頻率(階次)和濾波器長度分別為9.42和714。將此參數應用于CYCBD算法，解卷積結果如圖11(b)所示。相對于圖10(b)所示的解卷積結果，故障沖擊增強效果更為明顯。結果的平方包絡譜如圖11(c)所示。圖中行星輪軸承外圈故障特征階次及其諧波較為顯著。表明基于IFFR指標的粒子群算法能夠有效的獲取CYCBD算法中較佳的參數組合，進而利用參數優化的CYCBD算法可以準確提取行星輪軸承外圈故障特征。需要說明的是，試驗信號由于軸承滑移現象，行星輪軸承外圈實際故障階次與理論故障階次(9.57)存在約1.6%的誤差。

圖11 基于IFFR指標的CYCBD算法處理試驗信號Fig.11 IFFR index based CYCBD algorithm to process experimental signal

相對于圖10(c)結果，圖11(c)結果能夠更好地反映行星輪軸承外圈故障特征，而計算圖11(c)所示結果的FFR值僅為3.50，小于圖10(c)結果的FFR值(5.70)。可見在獲取CYCBD算法優化參數組合時，利用IFFR指標作為適應度函數，其效果要優于相同條件下的FFR指標。

6 結 論

本文針對現有CYCBD算法FFR指標存在的問題，提出一種基于參數優化最大二階循環平穩盲解卷積的行星輪軸承故障提取方法。對原FFR指標進行了改進，提出了IFFR指標。將IFFR作為粒子群算法的適應度函數對CYCBD算法的循環頻率和濾波器長度進行優化選取可以取得較好的效果。通過行星輪軸承外圈故障仿真和試驗驗證了方法的有效性和優勢。