唐振宇，黃凱，楊期江，朱曉彬

(廣州航海學(xué)院船舶與海洋工程學(xué)院，廣東廣州 510725)

0 前言

在旋轉(zhuǎn)機(jī)械轉(zhuǎn)子系統(tǒng)的運(yùn)行過程中，因受到測試環(huán)境、測試儀器和人為因素等影響，傳感器所采集的振動(dòng)信號(hào)易被噪聲污染，且旋轉(zhuǎn)機(jī)械結(jié)構(gòu)復(fù)雜，振動(dòng)信號(hào)往往由多個(gè)部件的振動(dòng)疊加而成。轉(zhuǎn)子振動(dòng)信號(hào)具有非高斯、非平穩(wěn)的特性。傳統(tǒng)的信號(hào)處理方法對(duì)含有噪聲的混疊信號(hào)進(jìn)行分析時(shí)有明顯不足，這也對(duì)機(jī)械的狀態(tài)監(jiān)測和故障診斷帶來了巨大的困難。

傳統(tǒng)的信號(hào)處理方法如小波分解、經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)已被廣泛應(yīng)用于旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)分析，但對(duì)于含噪的混疊信號(hào)，上述方法仍存在明顯不足。小波分解在對(duì)轉(zhuǎn)子振動(dòng)信號(hào)分析時(shí)存在以下問題:(1)小波分解層數(shù)的確定需要信號(hào)的頻率結(jié)構(gòu)和采樣頻率等先驗(yàn)知識(shí);(2)小波基和閾值的選擇困難。經(jīng)驗(yàn)?zāi)B(tài)分解也存在模態(tài)混疊和端點(diǎn)效應(yīng)問題，導(dǎo)致基于EMD的方法效果不佳。改進(jìn)的EMD方法如集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)是將噪聲輔助分析加在經(jīng)驗(yàn)?zāi)B(tài)分解中，有效地抑制了模態(tài)混疊的現(xiàn)象，但無法完全消除混疊現(xiàn)象。因此對(duì)于信號(hào)的分解，應(yīng)該根據(jù)信號(hào)本身的結(jié)構(gòu)特點(diǎn)，自適應(yīng)選擇基函數(shù)完成信號(hào)的分解。1993年MALLAT和ZHANG提出了基于過完備字典對(duì)信號(hào)進(jìn)行稀疏分解的思想，同時(shí)也提出了匹配追蹤算法。2004年有學(xué)者建立了壓縮感知理論(Compressing Sensing，CS)。壓縮感知理論證明了信號(hào)在某個(gè)變換空間內(nèi)具有稀疏性，可僅用少量數(shù)據(jù)保留信號(hào)的大部分特征，這樣可以用遠(yuǎn)低于奈奎斯特采樣定理要求的采樣頻率完成信號(hào)的重建恢復(fù)，這也將稀疏分解帶到了一個(gè)新高度。稀疏分解存在稀疏度過大，重構(gòu)信號(hào)的噪聲分量增加，而稀疏度過小，重構(gòu)信號(hào)又會(huì)丟失重要特征信息的問題。其中，稀疏度為信號(hào)中非零元素的個(gè)數(shù)。針對(duì)稀疏度選取問題的研究，文獻(xiàn)[6]提出了稀疏度自適應(yīng)匹配追蹤算法(Sparsity Adaptive Matching Pursuit,SAMP)，文獻(xiàn)[7]提出變步長自適應(yīng)匹配追蹤算法，文獻(xiàn)[8]提出基于CS的稀疏度變步長自適應(yīng)壓縮采樣匹配追蹤算法。

對(duì)于稀疏度的研究，國內(nèi)學(xué)者大多是對(duì)SAMP算法進(jìn)行改進(jìn)。本文作者提出互相關(guān)稀疏分解的方法，將互相關(guān)函數(shù)引入稀疏分解中以凸顯稀疏分解后各分量與原信號(hào)之間的相關(guān)性，并計(jì)算相關(guān)系數(shù)以解決稀疏度的選取問題。采用所提方法進(jìn)行仿真實(shí)驗(yàn)和轉(zhuǎn)子振動(dòng)信號(hào)實(shí)例分析，成功提取振動(dòng)信號(hào)周期特征。

1 互相關(guān)稀疏分解原理

1.1 信號(hào)的壓縮感知

以轉(zhuǎn)子系統(tǒng)為核心部件的旋轉(zhuǎn)機(jī)械運(yùn)行工況復(fù)雜，傳感器采集到的信號(hào)含有大量噪聲。轉(zhuǎn)子振動(dòng)信號(hào)數(shù)學(xué)模型為

()=()+()

(1)

式中：()為轉(zhuǎn)子振動(dòng)信號(hào)；()為噪聲。對(duì)于實(shí)信號(hào)()，如果它在稀疏字典中是稀疏的，則用觀測矩陣對(duì)信號(hào)進(jìn)行壓縮傳感，得到觀測信號(hào)。壓縮傳感方程表示為

=

(2)

有效的壓縮傳感要求觀測矩陣與稀疏字典互不相關(guān)，否則無法求解。高斯矩陣是壓縮傳感常用的測量矩陣，滿足RIP性質(zhì)。

1.2 稀疏分解理論

(3)

其中:為每個(gè)字典原子相對(duì)于原信號(hào)的權(quán)重，集合為稀疏系數(shù)，且只有個(gè)非零值。

1.3 正交匹配追蹤算法

常見的自然信號(hào)在時(shí)域是不稀疏的，信號(hào)稀疏分解理論指出，信號(hào)可以經(jīng)過變換矩陣進(jìn)行稀疏表示。由式(2)可以看作原信號(hào)在觀測矩陣下的線性投影，顯然的維數(shù)遠(yuǎn)低于的維數(shù)，用進(jìn)行稀疏分解更為容易。結(jié)合式(2)(3)得出下式:

=α

(4)

式中:=為×的矩陣(?)，稱為感知矩陣。此時(shí)滿足約束等距條件，可以通過求解式(5)的最優(yōu)化范數(shù)式重構(gòu)信號(hào)。

s.t.=

(5)

式(5)本質(zhì)上是一個(gè)優(yōu)化問題，即NP-hard問題。本文作者利用正交匹配追蹤算法(Orthogonal Matching Pursuit，OMP)選取多個(gè)字典原子計(jì)算信號(hào)最優(yōu)逼近，使得誤差最小。OMP通過求殘差與感知矩陣中各個(gè)原子之間內(nèi)積值的絕對(duì)值，來計(jì)算稀疏系數(shù)。內(nèi)積的絕對(duì)值越大說明感知矩陣中的這一列原子與信號(hào)越相關(guān)。

={|=|〈,〉|,=1,2,…,}

(6)

并采用最小二乘法對(duì)信號(hào)進(jìn)行逼近和殘差更新：

(7)

(8)

OMP算法具體步驟如下:

(1)初始?xì)埐?span id="caseq0s" class="emphasis_italic">=，迭代次數(shù)=1，索引值集合=、=;

(2)計(jì)算殘差與感知矩陣中各個(gè)原子之間內(nèi)積值的絕對(duì)值，并將中的最大值對(duì)應(yīng)的索引值存入中;

(3)更新感知矩陣，其中=∪;

(4)應(yīng)用公式(7)得到，同時(shí)用式(8)對(duì)殘差進(jìn)行更新;

1.4 過完備字典構(gòu)造

過完備字典是一組具有廣泛時(shí)頻特性的函數(shù)波形,每個(gè)波形作為字典中的一列，這一列稱為原子，原子的數(shù)量是冗余的?；谶^完備字典的稀疏分解能夠根據(jù)信號(hào)的結(jié)構(gòu)特征提取信號(hào)的各種局部特征，實(shí)現(xiàn)信號(hào)的自適應(yīng)分解，獲得信號(hào)的稀疏表示。

針對(duì)轉(zhuǎn)子振動(dòng)信號(hào)的周期特征，構(gòu)造離散傅里葉變換矩陣(Discrete Fourier Transform,DFT)，長度為的轉(zhuǎn)子振動(dòng)信號(hào)()通過DFT矩陣變換后轉(zhuǎn)化為()的形式:

(9)

即=，上述過程表示為矩陣形式為式(10)

(10)

(11)

公式(11)將其歸一化后可得到過完備字典。

1.5 皮爾遜相關(guān)系數(shù)

皮爾遜相關(guān)系數(shù)是用于度量2個(gè)向量之間線性相關(guān)性的變量，計(jì)算公式為

(12)

其中:為信號(hào)的點(diǎn)數(shù);為[-1 1]之間的數(shù)，越接近于1，則2個(gè)變量之間的相關(guān)性越高。

采用皮爾遜系數(shù)作為稀疏分解效果的評(píng)判標(biāo)準(zhǔn)。計(jì)算稀疏分解后的信號(hào)與原信號(hào)的相關(guān)性。一般情況下，當(dāng)>0.7時(shí)，認(rèn)為稀疏分解的信號(hào)與原信號(hào)相關(guān)性好，是真實(shí)的分量，予以保留。

2 仿真實(shí)驗(yàn)及分析

轉(zhuǎn)子振動(dòng)信號(hào)波形復(fù)雜，由其轉(zhuǎn)頻及倍頻和隨機(jī)噪聲等組成。依據(jù)上述特點(diǎn)，構(gòu)造仿真信號(hào)()，設(shè)置采樣頻率=2 000 Hz，采樣點(diǎn)數(shù)=1 200,則有：

()=sin(2π)+sin(2π)+sin(2π)+()+()

(13)

()=e-sin(2π)

(14)

其中：為轉(zhuǎn)子轉(zhuǎn)頻;、分別為2倍頻和3倍頻;()為轉(zhuǎn)子正常運(yùn)行時(shí)產(chǎn)生的沖擊衰減分量;()為均值為0且方差為1的高斯白噪聲。

設(shè)幅值=2 mm、=1.6 mm、=1.2 mm、=0.2 mm、=50 mm;基頻=30 Hz; 倍頻、分別為60、90 Hz;沖擊頻率=50 Hz。仿真信號(hào)()的時(shí)域波形如圖1所示，頻譜如圖2所示。

圖1 仿真信號(hào)X(t)的時(shí)域波形

圖2 仿真信號(hào)X(t)的頻譜

構(gòu)造大小為516×1 200的高斯觀測矩陣，對(duì)仿真信號(hào)進(jìn)行降維處理，經(jīng)式(2)得到大小為1×516的觀測信號(hào)。再構(gòu)造大小為1 200×1 200的離散傅里葉變換字典與觀測矩陣相乘得到感知矩陣?；诟兄仃?span id="so0qgcq" class="emphasis_italic">，采用正交匹配追蹤算法對(duì)觀測信號(hào)()進(jìn)行重構(gòu)，得到稀疏系數(shù)的集合。根據(jù)式(12)分別計(jì)算出不同稀疏度時(shí)，重構(gòu)信號(hào)()與仿真信號(hào)()的皮爾遜相關(guān)系數(shù)如表1所示。仿真信號(hào)()由3個(gè)主要分量組成，稀疏度取值應(yīng)大于3。稀疏度越大，重構(gòu)信號(hào)與仿真信號(hào)的相關(guān)性越強(qiáng)，當(dāng)值大到一定程度時(shí)，相關(guān)性并無明顯提高，甚至?xí)烊朐肼?。根?jù)實(shí)驗(yàn)結(jié)果并考慮時(shí)效性，重構(gòu)信號(hào)的稀疏度選取為6。=6時(shí)，重構(gòu)信號(hào)()如圖3所示，其頻譜如圖4所示。

表1 稀疏度K

圖3 重構(gòu)信號(hào)Sc(t)

圖4 重構(gòu)信號(hào)Sc(t)頻譜

由圖3、圖4可知：原始信號(hào)經(jīng)重構(gòu)后，噪聲和沖擊分量得到了有效抑制且重構(gòu)信號(hào)保留了原始信號(hào)的大部分特征，譜峰與原始信號(hào)的譜峰誤差精確到0.01 Hz?，F(xiàn)將稀疏系數(shù)中不為零的項(xiàng)與感知矩陣對(duì)應(yīng)原子相乘，分離仿真信號(hào)不同的頻率成分。稀疏度=6時(shí)，稀疏分解結(jié)果如圖5—圖7所示。經(jīng)MATLAB仿真計(jì)算，分離出的()～()與仿真信號(hào)()對(duì)應(yīng)的轉(zhuǎn)頻、倍頻成分相關(guān)性均達(dá)到0.99。

圖5 稀疏分解信號(hào)Sc0(t)

圖6 稀疏分解信號(hào)Sc1(t)

圖7 稀疏分解信號(hào)Sc2(t)

3 轉(zhuǎn)子振動(dòng)信號(hào)分析

3.1 軸承轉(zhuǎn)子實(shí)驗(yàn)臺(tái)

實(shí)驗(yàn)對(duì)象為Bently小型軸承轉(zhuǎn)子實(shí)驗(yàn)臺(tái)，該實(shí)驗(yàn)臺(tái)主要由驅(qū)動(dòng)系統(tǒng)、控制系統(tǒng)、采集系統(tǒng)及機(jī)械結(jié)構(gòu)等部分組成。為降低外界振動(dòng)對(duì)軸承轉(zhuǎn)子系統(tǒng)的影響，在實(shí)驗(yàn)平臺(tái)與軸承轉(zhuǎn)子系統(tǒng)之間增加具有良好減振性能的減振支座，最大程度地隔絕外界因素對(duì)實(shí)驗(yàn)臺(tái)的干擾。Bently小型軸承轉(zhuǎn)子實(shí)驗(yàn)臺(tái)如圖8所示，其中振動(dòng)采集系統(tǒng)為LMS多功能數(shù)據(jù)采集器。

圖8 轉(zhuǎn)子實(shí)驗(yàn)臺(tái)實(shí)物圖

3.2 轉(zhuǎn)頻和倍頻信號(hào)提取

使用LMS數(shù)據(jù)采集系統(tǒng)采集振動(dòng)加速度信號(hào)，采樣頻率為16 384 Hz，轉(zhuǎn)速為4 000 r/min。其時(shí)域波形、頻譜及時(shí)頻圖分別如圖9—圖11所示。在4 000 r/min工況下，信號(hào)D1的轉(zhuǎn)頻為65.6 Hz，頻譜內(nèi)存在高頻噪聲和轉(zhuǎn)子系統(tǒng)運(yùn)行時(shí)激發(fā)的諧波分量。

圖9 轉(zhuǎn)子振動(dòng)信號(hào)D1時(shí)域波形

圖10 轉(zhuǎn)子振動(dòng)信號(hào)D1頻譜

圖11 轉(zhuǎn)子振動(dòng)信號(hào)D1小波時(shí)頻圖

由圖10可知：可以提取4個(gè)分量，分別為(65.6 Hz,0.363 9 mm)、(131.2 Hz,0.269 7 mm)、(196.8 Hz,0.135 6 mm)、(326.4 Hz,0.074 51 mm)，濾除其余無用成分。

經(jīng)MATLAB計(jì)算，當(dāng)稀疏度=16時(shí)，重構(gòu)信號(hào)與信號(hào)D1的相關(guān)系數(shù)為0.9，接近于1，相關(guān)性高。文中取=16對(duì)信號(hào)進(jìn)行稀疏分解，結(jié)果如圖12—圖15所示。可以看出：文中方法準(zhǔn)確地提取了轉(zhuǎn)子振動(dòng)信號(hào)的轉(zhuǎn)頻及倍頻成分。

圖12 轉(zhuǎn)子轉(zhuǎn)頻

圖13 轉(zhuǎn)子2倍頻

圖14 轉(zhuǎn)子3倍頻

圖15 轉(zhuǎn)子5倍頻

以上實(shí)際工程應(yīng)用證明文中方法的有效性和可行性，為轉(zhuǎn)子運(yùn)行狀況監(jiān)測和故障診斷提供純凈的信號(hào)，有利于進(jìn)一步信號(hào)分析。

4 結(jié)論

為解決傳統(tǒng)信號(hào)處理方法難以對(duì)含噪混疊信號(hào)進(jìn)行分析的難點(diǎn)，本文作者提出互相關(guān)稀疏分解法。利用該方法對(duì)信號(hào)采用壓縮感知，僅需少量數(shù)據(jù)就可以保留信號(hào)的大部分特征，用遠(yuǎn)低于奈奎斯特采樣定理要求的采樣頻率完成信號(hào)的重建恢復(fù)，進(jìn)而引出稀疏度對(duì)重構(gòu)信號(hào)的影響。稀疏度過大，重構(gòu)信號(hào)的噪聲分量會(huì)增加，過小，重構(gòu)信號(hào)會(huì)丟失特征信息。因此，需選擇最合適的稀疏度完成對(duì)仿真信號(hào)的降噪重構(gòu)。

為了驗(yàn)證該方法的可行性，首先利用皮爾遜系數(shù)計(jì)算稀疏分解后的信號(hào)與原信號(hào)的相關(guān)性，并通過計(jì)算相關(guān)系數(shù)來解決稀疏度的取值問題。當(dāng)稀疏分解的信號(hào)與原信號(hào)相關(guān)性好(>0.7)時(shí)，保留該信號(hào)。將稀疏系數(shù)與字典中對(duì)應(yīng)的原子相乘，分離出仿真信號(hào)中不同的頻率成分，最后對(duì)轉(zhuǎn)子振動(dòng)信號(hào)實(shí)例進(jìn)行分析，成功提取轉(zhuǎn)子的轉(zhuǎn)頻及其倍頻成分。仿真實(shí)驗(yàn)結(jié)果驗(yàn)證了該方法的有效性和可行性，可為轉(zhuǎn)子運(yùn)行狀況監(jiān)測和故障診斷提供純凈的信號(hào)，有利于進(jìn)一步分析信號(hào)。結(jié)果表明：與傳統(tǒng)的信號(hào)處理方法相比，用互相關(guān)稀疏分解法處理的旋轉(zhuǎn)機(jī)械運(yùn)行的含噪混疊信號(hào)更易于分析，有助于旋轉(zhuǎn)機(jī)械的運(yùn)行狀況監(jiān)測和故障診斷。