高階局部投影算法及其在滾動軸承故障診斷中的應用

呂 勇， 施 威， 易燦燦

(武漢科技大學 機械自動化學院 冶金裝備及其控制教育部重點實驗室，武漢 430081)

滾動軸承已經廣泛地運用于各種旋轉機械設備中，由于軸承等關鍵零部件發生故障導致整個設備停止運行的現象層出不窮，由此造成了巨大的經濟損失。目前，由于一般很難知道機械系統具體的數學模型，因而對滾動軸承運行狀態的監測主要是通過傳感器測得反映其某一運行狀態的振動信號進行的。這些振動信號本質上是由速度、加速度或位移等單一變量組成的一維時間序列。然而，測得的振動信號中含有很多與軸承運動狀態無關的干擾信號，在信號的實際測量過程中，很容易受到環境噪聲的干擾，最終獲得的信號具有非平穩、非線性的特點，即呈現一定的混沌特性。由于在頻域中混沌信號與噪聲具有相似的寬頻特征，因此傳統的線性濾波方法具有很大的局限性。隨著非線性科學的不斷發展，基于相空間重構的降噪方法[1-5]逐漸成為研究非線性系統的重要工具之一。其中局部投影[6-9]是降噪效果比較好的一種方法，而且它不需要對系統模型有任何的先驗知識，因而廣泛地應用于各個領域[10-11]。

鄰域質心的選取是局部投影降噪過程中十分關鍵的一步，它決定著降噪效果的好壞。Cawley等將鄰域相點的平均值作為質心，這種方法由于局部線性化會產生較大的誤差，對降噪效果影響較大。后來Sauer為了進一步地抑制噪聲而利用二階多項式對鄰域質心進行了更精確地估計，使其與投影的超平面近似相切。最近，Moore等[12]提出對鄰域質心進行更高階地估計，提出了一種高階局部投影降噪算法，并將其運用于醫學信號處理中。雖然局部投影算法在機械故障診斷領域已有較為廣泛的應用，但大多是針對齒輪信號而且都是從時域提取特征。本文提出將高階局部投影單獨地應用于滾動軸承的故障診斷中，首先通過仿真信號驗證了該方法相比于傳統的基于相空間重構的非線性時間序列降噪方法的優越性。最后利用高階局部投影算法對從工業現場測得的風機軸承振動信號進行了降噪處理，降噪后從頻域成功地提取出了風機軸承故障特征。

1 局部投影算法基本原理

一個從混沌系統測得的長度為N的時間序列s(1),s(2),….,s(N)∈R,選擇合適的嵌入維數m和延遲時間τ就可以將其重構到m維的相空間中,則重構后相空間中的每個相點可以表示為

Xn=(sn,sn-τ,L,sn-(m-1)τ)

(1)

根據Takens[13]嵌入定理：當系統吸引子的分數維為d時，只需要滿足最小嵌入維數m≥2d+1時，重構后的相空間與原來的動力學系統是微分同胚的，即具有相同的動力學特性。在m維相空間中，反映系統動力學行為的吸引子通常只局限在一個m0的低維子空間(d

步驟1選擇合適的參數m和τ將含噪時間序列重構到m維相空間。

步驟2對于每個相點，確定其鄰域Un。鄰域的確定有固定鄰域數目和固定鄰域半徑兩種方式，本文采用的是前者。

步驟3計算每個鄰域相點的質心。

步驟4計算每個相點鄰域的協方差矩陣C并對其進行特征值分解，求出m-m0個較小特征值對應的特征向量aq，其中q=1,2,…,m-m0。

步驟5減去相點在噪聲子空間上的投影，即

(2)

步驟6返回步驟2，直到所有數據處理完畢。

2 高階局部投影算法

針對局部投影降噪步驟3，最初Cawley等將鄰域相點的平均值

(3)

作為鄰域質心，如圖1(a)所示，這種方法由于局部線性化會產生較大的誤差，各個鄰域的中心與曲線都不是相切的，而向曲線內部發生了一定程度的偏移，對降噪效果影響較大。后來Sauer為了進一步地抑制噪聲而利用二階多項式

(4)

對鄰域質心進行了估計，如圖1(b)所示，本質上是通過將鄰域中心向外移的方法使其與投影的超平面近似相切，本文將其稱為標準局部投影算法。

(a)局部線性化(b)Sauer的改進方法

圖1 質心估計方法示意圖

Fig.1 Schematic diagram of centroid estimation methods

Moore等使用更高階的多項式對鄰域質心進行更準確地估計，并將其運用于醫學信號處理中，取得了較好的效果。本文將其稱為高階局部投影算法，基本原理如下：

對于任意實數δ>0，定義連續移動平均算子Iδ

(5)

Sauer通過考慮將一階和二階的線性組合作為鄰域的質心，提高了局部投影算法的降噪功效，依此類推，如果令

(6)

適用于所有的單項式

(7)

這個公式是用來證明對于k=1,2,…,20，由k個線性方程組成的系統

n=0,2,…,2(k-1)

(8)

的唯一解是

(9)

顯然，由式(9)所確定的系數包括了Cawley等(i=1)和Sauer(i=2)提出的兩種估算鄰域質心的方法，而且還將鄰域質心的估計推廣到更高階(i≥3)。

綜上所述，高階局部投影算法的實現步驟如下：

步驟1確定合適的嵌入維數m和延遲時間τ，將一維含噪信號重構到m維相空間。

步驟2給定鄰域數目k，確定每個相點的鄰域Un。

步驟3利用前面提到的高階質心估計方法取代傳統的質心估計方法，計算出每個鄰域的質心。

步驟4計算每個相點的協方差矩陣C，并求得C的特征值和對應的特征向量。

步驟5確定m0個較大的特征值和剩余的m-m0個較小的特征值對應的特征向量。

步驟6利用式(2)進行修正，從而消除噪聲。

步驟7返回步驟2，直至所有相點處理完后將其還原成一維時間序列，這樣就得到了降噪后的信號。

3 仿真計算與工程應用

為了驗證本文所提方法的有效性和實用性，本文首先采用調頻信號和混沌信號進行數值仿真，并與奇異譜降噪、平滑局部子空間投影[14]和標準局部投影三種基于相空間重構的非線性時間序列降噪方法做了對比，最后對從工業現場測得的風機軸承數據進行了分析，并成功地提取出了風機軸承的故障特征。

由于這四種方法都是基于相空間重構，而嵌入維數m和延遲時間τ是其中兩個非常關鍵的參數。目前針對嵌入維數的選取方法有：偽最小鄰近法、關聯維數法和奇異值分解法等。文獻[15]中提出了一種嵌入維數自適應選取方法，并通過仿真和計算驗證了該方法的有效性，本文將其作為嵌入維數的選擇依據。而針對延遲時間的選取，通常使用平均互信息曲線第一次達到極小值時的τ值作為延遲時間。但是從數學理論的角度上講，τ可以任意取值，本文經過多次仿真實驗證明，τ值的選取對降噪效果影響很微弱。另外，τ取較小值時，鄰域點能夠較好地逼近吸引子的主流形，所以文中延遲時間τ均取值為1。

3.1 對調頻信號的數值仿真

采用的仿真信號為

x=cos(70πt+sin(20πt))

(10)

其中采樣點數N=8 000，采樣頻率fs=8 000 Hz，添加一定水平的高斯白噪聲，使其添加噪聲后(降噪前)的信噪比為3 dB。然后使用高階局部投影算法進行降噪處理。降噪效果如圖2所示，圖2(a)是添加一定噪聲后信噪比為3 dB的時域圖，圖2(b)是使用高階局部投影算法降噪后的時域圖。可以看到，降噪后的時域圖上基本看不出噪聲的影響，說明高階局部投影算法具有較好的降噪效果。

(a)含噪信號的時域圖(b)降噪后的時域圖

圖2 正弦信號降噪前后時域圖

Fig.2 Time waveforms of sinusoidal signal before and after noise reduction

為了進一步說明高階局部投影算法的降噪效果，將高階局部投影算法與奇異譜降噪、平滑局部子空間投影和標準局部投影降噪算法進行對比，使用降噪前后信噪比(Signal-to-Noise Ratio,SNR)的增益作為衡量降噪效果的標準。其中SNR計算式為

(11)

同樣采用之前提到的調頻信號進行仿真，添加不同水平的高斯白噪聲，使其添加噪聲后(降噪前)的信噪比分別為0 dB、3 dB、6 dB、9 dB、12 dB、15 dB和18 dB。然后對這7組數據分別使用這四種方法進行降噪處理。采樣點數N=8 000，采樣頻率fs=8 000 Hz。計算結果均是經過多次計算后取平均值以減小偶然誤差，圖3可以看出高階局部投影算法相比其它幾種非線性時間序列降噪方法具有更好的降噪效果。

圖3 四種降噪方法的對比結果

3.2 對混沌信號的數值仿真

本文采用Lorenz混沌時間序列進行仿真分析。其方程為

(12)

式中：σ=10,r=28,b=8/3，此時Lorenz系統呈現出混沌特性。取Lorenz系統的x分量進行仿真分析，采樣點數N=3 000，采樣頻率fs=300 Hz，設定積分步長為0.05，用四階Runge-Kutta法進行計算。對Lorenz系統的x分量添加一定水平的高斯白噪聲，使其添加噪聲后(降噪前)的信噪比為3 dB。然后使用高階局部投影算法進行降噪處理，降噪效果如圖4～圖6所示。

圖4 不含噪Lorenz信號的相圖

圖5 含噪Lorenz信號的相圖

圖6 降噪后Lorenz信號的相圖

圖4是未添加噪聲前的相圖，圖5是添加一定噪聲后信噪比為3 dB的相圖，圖6是使用高階局部投影降噪后的相圖。可以看到，即使在噪聲比較大的情況下，高階局部投影仍然具有較好的降噪效果，并且能夠很好地保持吸引子的幾何形狀。

為了進一步說明高階局部投影算法的降噪效果，將高階局部投影算法與奇異譜降噪、平滑局部子空間投影和標準局部投影算法對比。同樣采用之前提到的混沌信號進行仿真，添加不同水平的高斯白噪聲，使其添加噪聲后(降噪前)的信噪比分別為0 dB、3 dB、6 dB、9 dB、12 dB、15 dB和18 dB。然后對這7組數據分別使用這四種方法進行降噪處理。采樣點數N=8 000，采樣頻率fs=8 000 Hz，分別計算四種方法降噪后的信噪比。為了減小偶然誤差，計算結果均是經過多次計算后取平均值得到的。圖7可以看到高階局部投影具有最大的信噪比增益，說明了高階局部投影算法在保留吸引子確定性結構的基礎上對非線性信號具有較強的處理能力。

3.3 工程應用

大脫硫風機作為某煉鋼廠的關鍵設備，在環保除塵，治理環境污染方面起著重要的作用。在實際生產過程中一旦發生故障，將導致轉爐內的灰塵無法及時清除，不僅會造成環境的污染，而且有可能造成巨大的經濟損失和安全隱患。設備的結構簡圖如圖8所示，現場檢測人員使用美國CSI2130振動采集儀在軸承3的垂直方向進行信號采集，2006年9月發現該軸承出現異常的振動，將采集的信號分別使用高階局部投影算法和標準局部投影算法進行分析。軸承型號為22344CA，采樣頻率為25 600 Hz，轉速為740 r/min，這里取12 000個點進行分析，計算得到的故障特征頻率如表1所示。

圖7 四種降噪方法的對比結果

圖8 大脫硫風機的結構簡圖

轉頻f/Hz內圈fi/Hz外圈fo/Hz滾動體fr/Hz保持架fc/Hz12.3396.2064.1328.864.93

(a)原始信號時域圖(b)高階局部投影降噪后的時域圖

圖9 軸承振動信號降噪前后時域圖

圖10 軸承振動信號降噪前后時域圖

Fig.10 The spectrum of bearing vibration signal before and after noise reduction

從圖9和圖10，即時域圖和頻域圖中可以看出，經過本文所提出的方法降噪處理后，背景噪聲大大降低了。將頻域圖放大后，在沒有使用高階局部投影降噪處理前，由于受到大量噪聲的干擾，很多故障特征頻率被噪聲所淹沒，無法判斷出軸承的故障類型。而經高階局部投影降噪處理后，由圖11(b)可以清晰地觀察到軸的轉頻12.8 Hz，故障頻率64 Hz及其二倍頻128 Hz、三倍頻192 Hz和四倍頻256 Hz。圖12是經過標準局部投影降噪后的頻域放大圖，可以看到故障頻率64 Hz及其二倍頻128 Hz的幅值明顯低于圖11(b)中相對應頻率的幅值，因此高階局部投影算法能夠更加有效地提取故障特征。雖然與理論計算得到的外圈故障特征頻率值有微小的偏差，但在誤差允許的范圍內，符合要求。因此可以判斷是軸承外圈發生了故障，后來經檢測證實為外圈故障。

(a)原始信號頻譜放大圖(b)高階局部投影降噪后的頻譜放大圖

圖11 軸承振動信號降噪前后頻譜放大圖

Fig.11 The spectrum zoom of bearing vibration signal before and after noise reduction

圖12 標準局部投影降噪后的頻譜放大圖

通過對以上數值仿真信號和實測的風機軸承振動信號的分析表明：使用高階多項式對鄰域質心進行估計能夠更好地逼近吸引子的幾何結構，減少了吸引子的扭曲變形，進一步抑制了噪聲，提升了局部投影算法的降噪效果。

4 結 論

針對局部投影降噪算法過程中鄰域質心的選取對吸引子形狀的影響，使用高階多項式對鄰域質心進行了更精確地估計，減少了吸引子的扭曲變形，進一步抑制了噪聲，提升了局部投影算法的降噪效果。通過使用調頻信號和Lorenz混沌信號進行數值仿真并和其它基于相空間重構的非線性降噪方法對比，凸顯了高階局部投影算法的優越性。最后將其運用于工業現場風機軸承的故障診斷中，有效降低了背景噪聲，并從頻域成功提取出了軸承的故障特征，證明了高階局部投影算法的實用性和可行性。

[1] 陽建宏, 徐金梧, 呂勇, 等. 基于相重構和主流形識別的非線性時間序列降噪方法[J]. 北京科技大學學報, 2005,27(5):631-634.

YANG Jianhong, XU Jinwu, Lü Yong, et al. Nonlinear time series noise reduction method based on phase reconstruction and principal manifold learning[J].Journal of Beijing University of Science and Technology, 2005,27(5):631-634.

[2] CAWLEY R, HSU G H. Local-geometric-projection method for noise reduction in chaotic maps and flows[J]. Physical Review A, 1992,46(6):3057-3082.

[3] SAUER T. A noise reduction method for signals from nonlinear systems[J].Physical D,1992,58(1):193-201.

[4] KANTZ H, SCHREIBER T. Nonlinear time series analysis[M]. Cambridge: Cambridge University Press, 2003.

[5] 呂志民, 張武軍, 徐金梧. 基于奇異譜的降噪方法及其在故障診斷技術中的應用[J]. 機械工程學報, 2000,35(3):95-99.

Lü Zhimin, ZHANG Wujun, XU Jinwu. A noise reduction method based singular spectrum and its application in machine fault diagnosis[J].Journal of Mechanical Engineering, 2000,35(3):95-99.

[6] 徐金梧, 呂勇, 王海峰. 局部投影算法及其在非線性時間序列分析中的應用[J]. 機械工程學報, 2003,39(9):146-150.

XU Jinwu, Lü Yong, WANG Haifeng. Local projective method and its application on nonlinear time series[J].Journal of Mechanical Engineering,2003,39(9): 146-150.

[7] 陽建宏, 徐金梧, 楊德斌, 等. 鄰域自適應選取的局部投影非線性降噪方法[J]. 振動與沖擊, 2006,25(4):64-67.

YANG Jianhong, XU Jinwu, YANG Debin, et al. Nonlinear noise reduction method by local adaptive neighborhood selection[J].Journal of Vibration and Shock, 2006,25(4):64-67.

[8] 韓敏, 劉云俠. 基于非線性約束的局部投影降噪[J]. 電子與信息學報, 2009,31(2):400-404.

HAN Min, LIU Yunxia. Local projection noise reduction based on nonlinear constraints[J].Journal of Electronics and Information Technology,2009,31(2):400-404.

[9] MICHAEL T J, RICHARD J P. Generalized phase space projection for nonlinear noise reduction[J]. Physica D, 2005,201(3/4):306-317.

[11] KOTAS M. Projective filtering of time warped ECG beats[J]. Computers in Biology and Medicine, 2008,38(1):121-137.

[12] MOORE J M, SMALL M, KARRECH A. Improvements to local projective noise reduction through higher order and multiscale refinements[J]. Chaos, 2015,25(6): 653-671.

[13] TAKENS F. Detecting strange attractors in turbulence[J]. Lecture Notes in Mathematics,1981, 898(1):366-381.

[14] CHELIDZE D. Smooth local subspace projection for nonlinear noise reduction[J]. Chaos, 2014,24(1): 274-276.

[15] 余成義. 流形學習在混沌時間序列降噪中的應用[D]. 武漢: 武漢科技大學, 2012.