基于RPCA低秩稀疏分解的循環頻率檢測方法

王 冉，余龍靖，余 亮，蔣偉康

(1. 上海海事大學 物流工程學院, 上海 201306; 2. 上海交通大學 機械系統與振動國家重點實驗室，上海 200240)

循環平穩是處理隨機信號中隱藏的周期性現象的一個強大工具。循環平穩信號是指一類特殊的非平穩信號，其性質是隨機的，但在統計中表現出周期性[1]。旋轉機械由于其旋轉運動特點，其產生信號具備周期時變特性，是典型的循環平穩信號。目前，循環平穩信號處理方法已廣泛應用于旋轉機械故障診斷、機械系統識別和盲源分離等[2]。明陽等[3]提出了一種基于循環維納濾波器和包絡譜分析的方法對滾動軸承微弱故障特征進行提取。近年來通過改善循環平穩理論進行旋轉機械故障診斷仍是一個熱門話題[4-6]。Chen等提出了一種新的結合循環譜相干和深度學習的故障診斷方法，將滾動軸承故障信號的二維循環譜相干圖與卷積神經網絡(convolutional neural network, CNN)結合完成滾動軸承的故障診斷。

在循環平穩信號處理方法的應用中，循環頻率的檢測至關重要。例如：基于循環平穩分析實現機械源分離的前提是利用循環頻率設計濾波器。在旋轉機械的故障診斷中，循環頻率往往就是旋轉部件的故障特征頻率，對軸承、齒輪等部件的振動信號進行循環平穩分析的主要目的就是提取循環頻率從而確定設備的故障類型。在工程中的許多其他問題也與循環頻率的檢測息息相關。例如：Umut等[7]利用輪船聲信號具備的二階循環平穩特性，通過對水聽器測量的聲信號進行循環頻率檢測來探測輪船；Antoni等[8]利用循環調制相干性從測量信號的循環平穩特征中檢測水面艦船。此外，Li等[9]利用離心泵振動信號的流動誘導效應的循環平穩特性，使用循環平穩方法對其進行檢測。在上述應用中，為了檢測到有效的循環頻率及其高階諧波，常規的基于循環平穩特征的方案傾向于簡單地將信號特征的值與預定閾值進行比較。然而，當信號被噪聲干擾時，這樣的方案將導致檢測性能大幅下降。

主成分分析(principal component analysis, PCA)是數據科學領域中獲取和分析數據的基本工具，其基本思想就是從高維數據中提取出低維的有效數據。但主成分分析對于噪聲干擾非常敏感。魯棒主成分分析(robust principal component analysis, RPCA)是為了提高主成分分析的魯棒性而提出的用于矩陣低秩稀疏分解的方法。魯棒主成分分析的本質是將一個已知的采樣數據矩陣M∈M×N分解成一個低秩矩陣L∈M×N和一個稀疏矩陣S∈M×N，其中，M×N表示大小為M×N的實數矩陣。Candès等[10]通過將一般矩陣分解為低秩和稀疏矩陣解決了魯棒性問題，為信號處理和計算機視覺應用提供了合適的框架，并被廣泛應用于聲學信號處理[11-12]與機器學習[13-14]中。Yu等[15]利用RPCA對軸承振動信號短時傅里葉變換得到的時頻圖進行分解，從嚴重噪聲污染的滾動軸承振動信號中提取瞬態信號以進行故障診斷。

針對強噪聲干擾下難以準確檢測循環頻率的問題，本文提出了一種強干擾下將基于魯棒主成分分析(robust principal component analysis, RPCA)的低秩稀疏分解技術應用于循環譜密度(cyclic spectral density, CSD)矩陣從而進行循環頻率檢測的方法。該方法首先計算強背景噪聲下的循環平穩信號的CSD矩陣，利用RPCA將循環譜密度矩陣分解為代表噪聲干擾的低秩矩陣和包含循環頻率特征的稀疏矩陣，再利用稀疏矩陣構建檢測函數對循環頻率進行自動檢測。仿真結果驗證了該方法在強噪聲干擾下進行循環頻率檢測的優越性，并可根據檢測循環頻率各高階諧波的ROC曲線為不同信噪比條件下選擇檢測階數提供參考。最后將該方法應用在滾動軸承早期故障檢測上，試驗結果與常用的包絡譜和譜峭度(spectral kurtosis, SK)方法[16-18]對比以證明該方法在早期故障診斷應用中的優越性。

1 理論基礎

1.1 循環譜密度

循環平穩信號是一類特殊的非平穩信號，其統計特性隨時間周期性變化。包含額外信息的隱藏周期性表現在循環平穩信號中。

對于循環平穩信號x(t)，其時變自相關函數Rx(t,τ)為

(1)

Rx(t,τ)=Rx(t+T,τ)

(2)

(3)

式中，α為循環頻率，且α=m/T,m=1,2,3,…。對式(3)中的時延τ作Fourie變換即可得到CSD函數，也稱譜相關密度(spectral correlation density, SCD)函數

(4)

式中，f為頻譜頻率。CSD也可看作是時變自相關函數Rx(t,τ)的二階Fourier變換，因此也可寫為信號x(t)左右移頻α/2的兩信號的互譜密度

(5)

循環譜密度函數表示每個離散循環頻率分量α上頻譜頻率f的功率分布。

1.2 魯棒主成分分析

魯棒主成分分析算法是由Candès等提出的實現圖像背景與前景分離的重要技術之一。RPCA采用低秩稀疏模型，將數據矩陣分解為低秩矩陣和稀疏矩陣。設L為低秩矩陣，S為稀疏矩陣，M為原始矩陣，故RPCA 的低秩稀疏分解可以看成優化問題

(6)

(7)

2 基于RPCA的循環頻率檢測方法

對于一個循環頻率為αk(k=1,2,3,…)的循環平穩信號x(t)，其循環譜密度矩陣是由離散分布在對應循環頻率列上的稀疏值組成的，因此CSD中的循環頻率特征具有稀疏特性。且平穩噪聲的CSD具備低秩特性。本文針對噪聲和循環平穩信號的CSD矩陣分別具有的低秩和稀疏特性，提出一種對CSD矩陣低秩稀疏分解的方法。該方法利用RPCA對強干擾下循環平穩信號的CSD矩陣M進行低秩稀疏分解的方法，將噪聲的CSD建模為低秩矩陣L，循環平穩信號的CSD建模為稀疏矩陣S，進而利用稀疏矩陣S構建檢測器進行循環頻率的檢測。方法的具體流程如圖1所示。

圖1 基于RPCA的循環頻率檢測方法流程Fig.1 The flowchart of RPCA-based cyclic frequency detection method

對于輸入的帶有強干擾的循環平穩隨機信號，首先利用式(4)計算其循環譜密度矩陣。其次利用RPCA對CSD矩陣進行低秩稀疏分解。然而值得注意的是，稀疏矩陣S在λ較小時經常被噪聲污染(這表明一定量的噪聲被分解為稀疏分量)。因此，可以進一步使用濾波器來修正最終結果，其形式為

(8)

式中，ε為控制|S(m,n)|和|L(m,n)|幅度增益比的參數，本文中取定值ε=1。通過計算哈達瑪積(Hadamard product)E=F⊙S可以通過濾波減少稀疏矩陣S中的噪聲干擾從而得到稀疏性更強的增強稀疏矩陣E。最后，利用得到的增強稀疏矩陣E構建檢測函數

(9)

式中：TR∈Ck；分子表示增強稀疏矩陣E第k列元素之和；分母表示原始CSD矩陣M第k列的均值。最后通過一個帶有固定閾值的peak-picking算法檢測出檢測函數的峰值點即循環頻率。

3 仿真分析

為驗證方法的有效性，用帶有強噪聲干擾的循環平穩仿真信號輸入對其進行循環頻率檢測。仿真的循環平穩信號以軸承外圈故障模型為參考，其信號模型為

(10)

式中：s(t)為局部點蝕失效引起的靜態振蕩信號；T為振蕩發生周期；n(t)為干擾噪聲；β為共振阻尼衰減系數；fn為共振頻率。仿真時，采樣頻率為10 000 Hz，共有10 000個采樣點。振蕩發生周期T=0.05 s，故循環頻率為20 Hz。共振阻尼衰減系數β=500且共振頻率fn=3 600 Hz。

為驗證在低信噪比條件下的檢測性能，為仿真信號加以信噪比(signal-to-noise ratio，SNR)為-20～0 dB的高斯白噪聲分別進行檢測。圖2為加了RSN=-10 dB 噪聲的仿真信號的波形、循環譜密度圖。

圖2 加噪仿真信號(RSN=-10 dB)Fig.2 The noisy simulated signal (RSN=-10 dB)

由圖2(b)中加噪信號的CSD圖可以看得，在對仿真的循環平穩信號加入了RSN=-10 dB的高斯白噪聲后，循環平穩特征已被干擾成分所掩蓋。對該加噪信號矩陣進行RPCA低秩稀疏分解，其中參數λ=0.06。并進一步利用式(8)所示的濾波器對得到的稀疏分量進行增強，參數ε=1。分解得到的低秩分量和增強稀疏分量如圖3所示。

圖3 加噪仿真信號(RSN=-10 dB)低秩稀疏分解結果Fig.3 The low-rank and sparse decomposition results of noise-added simulation signal (RSN=-10 dB)

從分解結果可以看出在噪聲干擾的信噪比達到-10 dB 時仍可在稀疏分量中將背景噪聲幾乎分離，且在稀疏分量中可以明顯看出循環頻率及其高階諧波的位置。在得到帶有循環平穩特征的增強稀疏矩陣后，利用式(9)構建循環頻率的檢測函數。圖4(a)為直接使用原始CSD構建檢測函數的檢測結果，圖4(b)為檢測函數經過閾值為0.1的peak-picking算法檢測到的峰值點(即循環頻率及其高階諧波點)。由圖4可以看出在RSN=-10 dB時本文提出的檢測方案仍可以準確地檢測到循環頻率(20 Hz)及其前2～5階倍頻且相較于原始CSD直接檢測的效果優越。

圖4 循環頻率檢測結果Fig.4 The result of the cyclic frequency detection

為進一步評估所提出的循環頻率檢測器對于噪聲干擾的魯棒性，本文利用蒙特卡洛模擬，在各SNR噪聲干擾下分別對檢測循環頻率1～5階諧波頻率進行100次檢測模擬。通過Pd=Md/Mtot來計算檢測概率繪制ROC曲線，其中Md和Mtot為檢測到循環頻率及其高階諧波的次數和檢測的總次數。本文提出的檢測器的ROC曲線如圖5所示。結果表明了在強噪聲干擾下檢測概率方面的優越性。該ROC曲線也可為不同信噪比條件下選擇檢測階數提供參考。根據檢測循環頻率各高階諧波的ROC曲線，為保證檢測準確率高于90%的條件下盡可能多地檢測循環頻率階數，可在RSN≥-10 dB 時檢測5階諧波頻率，-11 dB≤RSN<-10 dB時檢測4階諧波頻率，-14 dB≤RSN<-11 dB時檢測3階諧波頻率，-15 dB≤RSN<-14 dB時檢測1階諧波頻率。

圖5 蒙特卡洛模擬ROC曲線Fig.5 The ROC curve of Monte Carlo simulation

4 試驗驗證

滾動軸承故障信號是典型的循環平穩信號，故利用本文提出的方法可以被應用到滾動軸承故障診斷中。為進一步驗證方法的有效性，本文將所提出的方法應用于滾動軸承加速疲勞壽命試驗數據中以實現滾動軸承早期故障的檢測。

4.1 試驗設置

本試驗為在杭州軸承試驗研究中心支持的加速軸承壽命測試儀(ABLT-1A)上進行的加速軸承從運行到故障的全壽命測試。試驗裝置主要由試驗臺、滾動軸承架、傳動系統、加載系統、計算機控制系統和數據采集系統組成。整個試驗設備如圖6所示。試驗臺在由交流電機驅動的單軸上同時承載4個滾動軸承。傳感器安裝簡圖如圖7所示。軸轉速設置為3 000 r/min，軸上的徑向載荷增加為12 744 kN。每分鐘采集一組 0.8 s的數據，采樣率為25.6 kHz。被測軸承型號為6307，根據軸承的結構參數和轉速，可由運動學公式計算出軸承失效的特征頻率，如表1所示。

圖6 ABLT-1A型軸承壽命強化試驗機Fig.6 Accelerated bearing life tester, ABLT-1A

圖7 傳感器安裝簡圖Fig.7 The installation diagram of the sensor

表1 6307滾動軸承特征頻率Tab.1 Characteristic frequencies of rolling bearing 6307

4.2 試驗結果

在加速軸承故障壽命測試中，總共測試了8個滾動軸承。作為第一個故障軸承，軸承B12不易受到其他軸承的干擾。在我們的試驗中，選擇失效滾動軸承B12的數據作為試驗對象。圖8為B12軸承的全壽命有效值(RMS)變化曲線(RMS是工程上使用較為廣泛的設備衰退評估指標)。由圖8可知，軸承RMS值在2 302點發生突變。

圖8 B12軸承全壽命有效值RMSFig.8 The RMS of B12 bearing life cycle

在這種情況下，將考慮整個壽命周期的第2 297 min 的數據。圖9顯示了2 297 min原始信號的時域波形和包絡譜分析。顯然，無法從波形中獲取軸承的故障信息。分析信號的包絡譜如圖9(a)所示，由包絡譜只能看見內圈的故障特征頻率一階基頻且噪聲干擾成分大，影響故障檢測。

圖9 B12軸承全壽命周期的第2 297 min的數據Fig.9 The data at the 2 297th minute of the full life cycle of the B12 bearing

使用所提出的方法對分析信號進行循環頻率檢測，其CSD與分解后的增強稀疏分量如圖10所示。如圖10(a)所示，原始信號的CSD中軸承的的故障特征成分十分微弱且被強干擾成分掩蓋，無法從中檢測故障。經本文所提出方法低秩稀疏分解后并增強后的稀疏分量如圖10(b)所示，其中參數λ=0.025，ε=1。從圖中可以明顯地看出軸承內圈故障特征成分以及轉頻的邊頻帶成分。

圖10 分析信號CSD圖與低秩稀疏分解結果Fig.10 The CSD of analysis signal and low-rank and sparse decomposition results

最后，對于得到的增強稀疏分量，使用式(9)計算檢測函數并采用閾值為0.1的peak-picking算法進行檢測，檢測的結果如圖11所示。

圖11 循環頻率檢測結果Fig.11 The result of the cyclic frequency detection

在尋峰算法找到的峰值中，可以找到內圈故障頻率fbpi及其倍頻成分，且可以明顯地看到轉頻調制的邊頻帶成分，雖有得到一些干擾頻率成分，但已可以充分地檢測到軸承的內圈故障。將本文提出方法的檢測結果對比于圖9(b)的包絡譜結果，其效果有明顯的提升。該方法在2 297 min即可檢測到軸承故障，比利用RMS值檢測故障的時刻早。這充分地說明該方法在滾動軸承早期故障檢測中的優越性。

為進一步證明該方法在滾動軸承早期故障檢測中的優越性，本文使用譜峭度的方法對第2 297 min數據故障特征提取的結果作為對比。

與共振頻帶相關的信息最豐富的頻帶可以從圖12(a)中所示的譜峭度圖中獲得，由虛線圓圈突出顯示。信息最豐富的頻譜頻帶的中心頻率fc=2 333.3 Hz，帶寬為bw=133.3 Hz。隨后利用帶通濾波器將該信息最豐富的頻帶濾出，濾波后信號的時域波形如圖12(b)所示。接著對該濾波信號進行包絡譜分析，如圖12(c)所示。由濾波信號包絡譜可以看到，經譜峭度方法濾波后分析的方法得到的包絡譜中僅可以看到內圈故障特征頻率的一階諧波，在更高階的頻率上并沒有峰值存在。將此結果對比圖11中本文所提出方法的檢測函數，可以進一步反映出本文提出方法的優越性。

圖12 譜峭度分析Fig.12 The results of spectral kurtosis

5 結 論

針對低信噪比條件下循環頻率難以檢測的局限，本文提出了一種低信噪比下將基于魯棒主成分分析的低秩稀疏分解技術應用于循環譜密度矩陣從而進行循環頻率檢測的新方法。利用低秩稀疏分解出的帶有循環平穩特征的稀疏矩陣構建檢測函數對循環頻率進行檢測。仿真結果表明該方法在強噪聲干擾下循環頻率檢測的優越性：在信噪比為-10 dB下仍精準可檢測到前5階循環頻率。根據蒙特卡洛模擬得到ROC曲線可以為不同信噪比下檢測不同階數的循環頻率提供參考。該方法也可進一步應用到滾動軸承的早期故障檢測上，通過軸承加速疲勞壽命試驗數據進行驗證，試驗結果表明：該方法在第2 297 min可以有效檢測到軸承內圈故障，對比傳統包絡譜和譜峭度的方法有明顯的優勢。