高階局部投影算法及其在滾動(dòng)軸承故障診斷中的應(yīng)用

呂 勇， 施 威， 易燦燦

(武漢科技大學(xué) 機(jī)械自動(dòng)化學(xué)院 冶金裝備及其控制教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430081)

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

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

1 局部投影算法基本原理

一個(gè)從混沌系統(tǒng)測(cè)得的長(zhǎng)度為N的時(shí)間序列s(1),s(2),….,s(N)∈R,選擇合適的嵌入維數(shù)m和延遲時(shí)間τ就可以將其重構(gòu)到m維的相空間中,則重構(gòu)后相空間中的每個(gè)相點(diǎn)可以表示為

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

(1)

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

步驟1選擇合適的參數(shù)m和τ將含噪時(shí)間序列重構(gòu)到m維相空間。

步驟2對(duì)于每個(gè)相點(diǎn)，確定其鄰域Un。鄰域的確定有固定鄰域數(shù)目和固定鄰域半徑兩種方式，本文采用的是前者。

步驟3計(jì)算每個(gè)鄰域相點(diǎn)的質(zhì)心。

步驟4計(jì)算每個(gè)相點(diǎn)鄰域的協(xié)方差矩陣C并對(duì)其進(jìn)行特征值分解，求出m-m0個(gè)較小特征值對(duì)應(yīng)的特征向量aq，其中q=1,2,…,m-m0。

步驟5減去相點(diǎn)在噪聲子空間上的投影，即

(2)

步驟6返回步驟2，直到所有數(shù)據(jù)處理完畢。

2 高階局部投影算法

針對(duì)局部投影降噪步驟3，最初Cawley等將鄰域相點(diǎn)的平均值

(3)

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

(4)

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

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

圖1 質(zhì)心估計(jì)方法示意圖

Fig.1 Schematic diagram of centroid estimation methods

Moore等使用更高階的多項(xiàng)式對(duì)鄰域質(zhì)心進(jìn)行更準(zhǔn)確地估計(jì)，并將其運(yùn)用于醫(yī)學(xué)信號(hào)處理中，取得了較好的效果。本文將其稱為高階局部投影算法，基本原理如下：

對(duì)于任意實(shí)數(shù)δ>0，定義連續(xù)移動(dòng)平均算子Iδ

(5)

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

(6)

適用于所有的單項(xiàng)式

(7)

這個(gè)公式是用來證明對(duì)于k=1,2,…,20，由k個(gè)線性方程組成的系統(tǒng)

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

(8)

的唯一解是

(9)

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

綜上所述，高階局部投影算法的實(shí)現(xiàn)步驟如下：

步驟1確定合適的嵌入維數(shù)m和延遲時(shí)間τ，將一維含噪信號(hào)重構(gòu)到m維相空間。

步驟2給定鄰域數(shù)目k，確定每個(gè)相點(diǎn)的鄰域Un。

步驟3利用前面提到的高階質(zhì)心估計(jì)方法取代傳統(tǒng)的質(zhì)心估計(jì)方法，計(jì)算出每個(gè)鄰域的質(zhì)心。

步驟4計(jì)算每個(gè)相點(diǎn)的協(xié)方差矩陣C，并求得C的特征值和對(duì)應(yīng)的特征向量。

步驟5確定m0個(gè)較大的特征值和剩余的m-m0個(gè)較小的特征值對(duì)應(yīng)的特征向量。

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

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

3 仿真計(jì)算與工程應(yīng)用

為了驗(yàn)證本文所提方法的有效性和實(shí)用性，本文首先采用調(diào)頻信號(hào)和混沌信號(hào)進(jìn)行數(shù)值仿真，并與奇異譜降噪、平滑局部子空間投影[14]和標(biāo)準(zhǔn)局部投影三種基于相空間重構(gòu)的非線性時(shí)間序列降噪方法做了對(duì)比，最后對(duì)從工業(yè)現(xiàn)場(chǎng)測(cè)得的風(fēng)機(jī)軸承數(shù)據(jù)進(jìn)行了分析，并成功地提取出了風(fēng)機(jī)軸承的故障特征。

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

3.1 對(duì)調(diào)頻信號(hào)的數(shù)值仿真

采用的仿真信號(hào)為

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

(10)

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

(a)含噪信號(hào)的時(shí)域圖(b)降噪后的時(shí)域圖

圖2 正弦信號(hào)降噪前后時(shí)域圖

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

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

(11)

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

圖3 四種降噪方法的對(duì)比結(jié)果

3.2 對(duì)混沌信號(hào)的數(shù)值仿真

本文采用Lorenz混沌時(shí)間序列進(jìn)行仿真分析。其方程為

(12)

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

圖4 不含噪Lorenz信號(hào)的相圖

圖5 含噪Lorenz信號(hào)的相圖

圖6 降噪后Lorenz信號(hào)的相圖

圖4是未添加噪聲前的相圖，圖5是添加一定噪聲后信噪比為3 dB的相圖，圖6是使用高階局部投影降噪后的相圖?？梢钥吹剑词乖谠肼暠容^大的情況下，高階局部投影仍然具有較好的降噪效果，并且能夠很好地保持吸引子的幾何形狀。

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

3.3 工程應(yīng)用

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

圖7 四種降噪方法的對(duì)比結(jié)果

圖8 大脫硫風(fēng)機(jī)的結(jié)構(gòu)簡(jiǎn)圖

轉(zhuǎn)頻f/Hz內(nèi)圈fi/Hz外圈fo/Hz滾動(dòng)體fr/Hz保持架fc/Hz12.3396.2064.1328.864.93

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

圖9 軸承振動(dòng)信號(hào)降噪前后時(shí)域圖

圖10 軸承振動(dòng)信號(hào)降噪前后時(shí)域圖

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

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

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

圖11 軸承振動(dòng)信號(hào)降噪前后頻譜放大圖

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

圖12 標(biāo)準(zhǔn)局部投影降噪后的頻譜放大圖

通過對(duì)以上數(shù)值仿真信號(hào)和實(shí)測(cè)的風(fēng)機(jī)軸承振動(dòng)信號(hào)的分析表明：使用高階多項(xiàng)式對(duì)鄰域質(zhì)心進(jìn)行估計(jì)能夠更好地逼近吸引子的幾何結(jié)構(gòu)，減少了吸引子的扭曲變形，進(jìn)一步抑制了噪聲，提升了局部投影算法的降噪效果。

4 結(jié) 論

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

[1] 陽建宏, 徐金梧, 呂勇, 等. 基于相重構(gòu)和主流形識(shí)別的非線性時(shí)間序列降噪方法[J]. 北京科技大學(xué)學(xué)報(bào), 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] 呂志民, 張武軍, 徐金梧. 基于奇異譜的降噪方法及其在故障診斷技術(shù)中的應(yīng)用[J]. 機(jī)械工程學(xué)報(bào), 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] 徐金梧, 呂勇, 王海峰. 局部投影算法及其在非線性時(shí)間序列分析中的應(yīng)用[J]. 機(jī)械工程學(xué)報(bào), 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] 陽建宏, 徐金梧, 楊德斌, 等. 鄰域自適應(yīng)選取的局部投影非線性降噪方法[J]. 振動(dòng)與沖擊, 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]. 電子與信息學(xué)報(bào), 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] 余成義. 流形學(xué)習(xí)在混沌時(shí)間序列降噪中的應(yīng)用[D]. 武漢: 武漢科技大學(xué), 2012.