多通道多分量分解方法在變轉速工況齒輪故障特征提取中的應用

張 亢 田澤宇 陳向民 廖力達 吳家騰

1.長沙理工大學能源與動力工程學院，長沙，4101142.湖南大學機械與運載工程學院，長沙，410082

0 引言

旋轉機械的振動信號包含了機械設備中零部件的健康狀態(tài)信息。對振動信號進行分析和處理是從振動信號中提取故障特征的重要手段，而故障特征是對旋轉機械進行故障狀態(tài)判斷與預測的依據(jù)[1]。齒輪是應用廣泛的旋轉機械部件，當其發(fā)生故障時，故障激勵會導致振動信號產(chǎn)生調(diào)制現(xiàn)象，使得振動信號呈現(xiàn)非平穩(wěn)的多分量調(diào)制特征[2-3]。要提取故障特征就需要減小噪聲和其他分量的干擾，從中分離出調(diào)制成分。然而實際中，齒輪所在齒輪箱結構的復雜性以及運行工況的波動性會使傳感器測得的振動信號包含復雜的耦合成分[4]，且振動信號的分量可能跨多個頻帶或同一頻帶包含多個分量，以及分量在支撐域重疊，從而造成調(diào)制成分分離困難[5]。

目前主要通過對振動信號進行分解來分離調(diào)制成分。根據(jù)信號分解實現(xiàn)方式的不同，現(xiàn)有分解方法主要可分為時域分解方法、頻域濾波分解和時頻域重構方法。經(jīng)驗模態(tài)分解[6]、局域均值分解[7]、本征時間尺度分解[8]等遞歸式時域分解方法得到了大量應用，這一類方法以信號的時間尺度特征為基礎，而時頻域重疊的多分量信號在同一時間尺度上分布著多個分量，此類方法會將其識別為單個分量而無法分解；頻域濾波方法以頻帶劃分為基礎，如奇異譜分解[9]、小波包變換[10]、變分模態(tài)分解[11]等，該類方法對頻帶的劃分是機械的格型劃分，缺乏自適應性，另外對于跨多個頻帶的分量信號，會被分解到不同的頻帶之中，從而失去物理意義；時頻域重構方法通過將目標分量所在的時頻支撐域之外的時頻分布系數(shù)置零，達到分離各分量的目的，如同步壓縮變換[12-13]等。這類方法具有銳化的時頻脊線，提高了時頻的可讀性，然而當信號分量在時頻域發(fā)生交叉時，由于交叉部位擁有相同的時頻分布系數(shù)，故無法對其進行分解或錯誤分解。

根據(jù)上述分析，對于分解具有復雜時頻特征的非平穩(wěn)信號，不能僅僅以信號的時頻域局部特征為基礎，而是要找到既能獨立代表分量又不受時頻特征影響的物理量作為分解的依據(jù)。最近，STANKOVIC等[14]提出了一種新的多通道多分量非平穩(wěn)信號分解方法(MMD)，該方法將信號中的分量用一組權重系數(shù)加權的向量線性組合表示，并以時頻聚集性度量值為優(yōu)化目標，通過迭代優(yōu)化得到各個分量的權重系數(shù)組，從而獲得相應的分量信號。MMD方法與前述以時頻局部特征為分解依據(jù)的分解方法相比，其本質不同在于將信號中的分量表示成一個由權重系數(shù)加權的向量，而不同的分量具有不同的系數(shù)組合，因此，理論上，對跨多個頻帶或單個頻帶中有多個分量以及分量信號發(fā)生了交叉的復雜信號也能進行分解，特別適合具有多分量調(diào)制特點且由于多激勵源耦合造成時頻分布復雜的齒輪箱振動信號的分析。目前，MMD方法僅應用于仿真信號的分析[15-16]，而對于高采樣率的實際振動信號，進行時頻聚集性度量的計算量很大，從而導致MMD算法的分解時間長，實用性差。因此，要利用MMD方法分析實際振動信號，必須解決原算法中時頻聚集性度量計算量大的問題。對數(shù)窗能量準則[17]通過折中選取最佳窗長，降低時頻矩陣的規(guī)模，減少算法迭代次數(shù)，從而可間接降低時頻聚集性度量的計算量與時長。

綜上，本文在將對數(shù)窗能量準則融入MMD方法、提高MMD方法計算效率的基礎上，通過仿真信號將MMD方法與傳統(tǒng)的基于時間尺度的分解方法進行了對比，證明了其優(yōu)勢。進一步，利用這種全新的方法，從復雜的變轉速工況齒輪箱故障振動信號中分離含有故障調(diào)制成分的分量信號，并從中提取故障特征，從而為變轉速工況齒輪故障特征提取提供了一條新的途徑。

1 MMD方法

對于任意一組多通道多分量信號x(n)，MMD方法將x(n)中的每一個分量信號看作若干個特征向量的線性組合，通過優(yōu)化迭代獲得該線性組合的權重系數(shù)，便可得到x(n)相應的分量信號，其具體步驟如下：

(1)將x(n)寫成

(1)

其中，x(m)(n)(m=1,2,…,S)代表第m個通道信號(即傳感器測得的數(shù)據(jù))，每一個通道信號實際上都是一組多分量信號；xp(n)=Ap(n)ejφp(n)(p=1,2,…,P)表示第p個分量成分(即系統(tǒng)產(chǎn)生的源信號)，Ap(n)和φp(n)分別為瞬時幅值和瞬時相位；αmp和φmp分別為源信號經(jīng)傳感器改變的幅值和相位。通過式(1)可以看出任意一個傳感器信號可以由源信號分量的線性組合來表示。

將式(1)進一步寫成矩陣形式：

Xsen=AXcom

(2)

其中，xp(p=1,2,…,P)表示分量信號xp(n)的向量形式。傳感器信號矩陣為

Xsen=[x(1)(n)x(2)(n) …x(S)(n)]T

系統(tǒng)產(chǎn)生的分量信號矩陣為

Xcom=[x1(n)x2(n) …xP(n)]T=

[x1x2…xP]T

常系數(shù)轉換矩陣為

將系統(tǒng)產(chǎn)生的信號分量轉換為傳感器測得的信號。

(2)計算Xsen的自相關矩陣R：

(3)

其中，(·)H表示矩陣的Hermitian轉置，自相關矩陣R為一個方陣。

(3)對自相關矩陣R進行特征值分解：

(4)

其中，矩陣U為特征向量矩陣；Λ為特征值對角陣；λp為特征值；up為對應的特征向量。由式(3)和式(4)可以看出，特征向量up可以用傳感器信號的線性組合表示，由上文已知傳感器信號是源信號分量的線性組合，從而可以將特征向量up用源信號分量的線性組合表示，即

(5)

當S≥P時，由式(5)可以知道，每個分量信號xp同樣可以由特征向量的線性組合表示，即

xp=η1pu1+η2pu2+…+ηPpuP

(6)

其中，ηip(i=1,2,…,P)為未知的權重系數(shù)。

(4)因為不同信號成分在時頻面上的時頻聚集性不同，所以可以用時頻聚集性來判斷信號成分。信號在時頻域的支撐區(qū)間的面積可以定量描述信號的時頻聚集性。設每個分量信號xp在時頻域的支撐區(qū)間為DTp，其面積為Dp，如果分量有重疊，那么其支撐區(qū)間必然也會重疊，但實際中支撐區(qū)間一般不會完全重疊，因此不失一般性，假設D1≤D2≤…≤DP。信號時頻分布的Lq范數(shù)可用來表示信號在時頻面上支撐區(qū)間的面積，即Lq范數(shù)可以作為聚集性測度值M{·}來度量信號的時頻聚集性，即

(7)

其中TFR(n,k)表示信號的時頻分布。因為短時傅里葉變換(short time Fourier transform，STFT)具有計算簡單、不會產(chǎn)生交叉項干擾等優(yōu)點，且得到了成熟應用，因此，本文將STFT作為信號的時頻表示。

(8)

i=1,2,…,P

那么將此時的權重系數(shù)ηi1代入式(6)，便可得到第一個分量信號x1。

式(8)實際上是一個優(yōu)化問題，即以y的時頻聚集性測度最小為優(yōu)化目標，在保證‖y‖2恒定的情況下，利用具體的優(yōu)化算法計算相應的未知權重系數(shù)ηi1。本文采用遺傳算法來求解該優(yōu)化問題。

(6)由式(6)可知，x1=η11u1+η21u2+…+ηP1uP，將第一個分量信號x1替換為第一個特征向量u1，然后將u1(即x1)從剩余的特征向量uk(k=2,3,…,P)中去除，即

(9)

此時得到了一組新的特征向量us(s=1,2,…,P)，形成新的線性組合y=β1u1+β2u2+…+βPuP；對y重復步驟(5)，可得到第二小的支撐區(qū)間面積D2對應的權重系數(shù)ηi2(i=1,2,…,P)，將ηi2代入式(6)即可得到第二個分量信號x2；將x2替換特征向量u2，并從剩余的特征向量uk(k=3,4,…,P)中去除，依此類推，整個過程是一個迭代優(yōu)化的過程，直到?jīng)]有新的特征向量us產(chǎn)生，代表信號中所有的分量xp(p=1,2,…,P)被分解出來，此時迭代終止。

從上述分解過程可以看出，每一個分量信號都是不同權重系數(shù)加權的特征向量的線性組合，其中特征向量來自于信號自身構成的矩陣，即整個分解過程是不需要先驗知識的，而表示每個分量的特征向量的線性組合是唯一的，這確保了得到的分量信號是不重復的；此外，以時頻聚集性測度最小為優(yōu)化目標，通過迭代優(yōu)化獲得權重系數(shù)，這保證了每一個特征向量的線性組合只對應了一個單分量信號，而不是交叉分量或組合分量等虛假分量。

2 時頻聚集性度量與仿真信號分析

2.1 基于對數(shù)窗能量準則的時頻聚集性度量

在MMD方法中，一方面，時頻聚集性度量的準確性決定了分解的效果；另一方面，以信號時頻矩陣的范數(shù)值作為時頻聚集性度量值，時頻矩陣的規(guī)模會影響算法的計算量。而對于實際振動信號，因為其較高的采樣頻率會導致大維數(shù)的時頻矩陣，所以會有很大的計算量。因此，如何保證時頻聚集性度量準確的前提下又具有小的計算量，是MMD方法應用于實際振動信號分析的關鍵。

良好的時頻分辨率是時頻聚集性度量準確的前提[18]，STFT的時頻分辨率與窗長N和采樣頻率fs密切相關。對于實測振動信號，采樣頻率fs是預先設定的，窗長N決定了時頻分辨率。本文通過對數(shù)窗能量準則Q(N)來自適應確定最優(yōu)窗長：

(10)

其中，G(n,k)為信號的STFT時頻分布;lgE(N)為窗能量E(N)的對數(shù)值;Mw為窗函數(shù)個數(shù)。因此，Q(N)實際上為信號STFT時頻能量與窗函數(shù)能量的對數(shù)的比值。隨著窗長N的增大，式(10)中分子分母增長速度的快慢不同，那么理論上Q(N)會在某些窗長發(fā)生折中，而當Q(N)在某個折中位置為最小值時，意味著STFT的時頻聚集性在窗長范圍內(nèi)最優(yōu)，那么此時對應的窗長N為最佳窗長[17-18]。

2.2 仿真信號分析

仿真信號s(t)如下：

(11)

它由3個分量組成，其中調(diào)幅分量s1(t)與線調(diào)頻分量s2(t)在時頻域發(fā)生了交叉，s3(t)為跨時頻尺度的調(diào)制分量。

首先通過對數(shù)窗能量準則確定最優(yōu)窗長，為考察準則的魯棒性，設置采樣頻率分別為800 Hz、1000 Hz和1200 Hz，窗函數(shù)選用Hanning窗，窗長范圍設定為8～256。根據(jù)式(10)計算不同采樣頻率下的對數(shù)窗能量準則曲線，如圖1所示，可以看出，800 Hz、1000 Hz、1200 Hz采樣頻率下的最優(yōu)窗長分別為95、119和137。

根據(jù)最優(yōu)窗長對s(t)進行STFT，同時在相同采樣頻率下，選擇一個較小和一個較大的窗長進行對比，結果如圖2所示?？梢悦黠@看出,在不同采樣頻率下，都是根據(jù)最優(yōu)窗長進行STFT時獲得的時頻效果最佳；另外，在采樣頻率為1000 Hz的情況下，分別選擇不同窗長，在同一臺計算機上針對仿真信號將MMD算法各運行了30次。統(tǒng)計循環(huán)迭代次數(shù)和平均計算時間兩個指標，如表1所示。可以看出，相比于其他窗長，當采用對數(shù)窗能量準則確定的窗長119時，MMD算法的平均運行時間最短，迭代次數(shù)最少，說明了引入對數(shù)窗能量準則的必要性。

圖1 s(t)在不同采樣頻率下的對數(shù)窗能量準則曲線Fig.1 Logarithm window energy rule curve of s(t) atdifferent sampling frequencies

(a)fs=800 Hz，N=50 (b)fs=800 Hz,N=95 (c)fs=800 Hz,N=180

(d)fs=1000 Hz,N=95 (e)fs=1000 Hz,N=119 (f)fs=1000 Hz,N=209

(g)fs=1200 Hz,N=119 (h)fs=1200 Hz,N=137 (i)fs=1200 Hz,N=210

表1 不同窗長下MMD算法的迭代次數(shù)與運算時間

基于上述分析，將仿真信號s(t)的采樣頻率設置為1000 Hz，對應最優(yōu)窗長為119，進行MMD分解。得到的自相關矩陣的特征值如圖3所示，可以看出,非零特征值的個數(shù)為3，即P=3，正好等于s(t)的實際分量數(shù)。根據(jù)S≥P，取S=3，經(jīng)過迭代分解，結果如圖4所示?？梢钥闯觯徊娣至康玫搅撕芎玫姆蛛x，得到的第三代特征向量準確地對應了仿真信號的3個分量，說明了MMD方法的有效性。

圖3 自相關矩陣的特征值Fig.3 Eigenvalues of the autocorrelation matrix

(a)初始特征向量

(b)第一代特征向量

(c)第二代特征向量

(d)第三代特征向量圖4 仿真信號s(t)的MMD分解結果Fig.4 MMD decomposition results of simulationsignal s(t)

采用目前應用廣泛的經(jīng)驗模態(tài)分解(empirical mode decomposition，EMD)和變分模態(tài)分解(variational mode decomposition，VMD)方法分解仿真信號s(t)，為直觀對比，分解結果也用STFT時頻圖表示。EMD與VMD分解得到的前3個分量如圖5和圖6所示。由圖5可以看出，EMD將原信號中的交叉分量當成了一個分量，而非交叉分量也由于時間尺度跨度大，被分解為兩個分量；由圖6可知，VMD算法的頻帶劃分清晰，原信號按頻率被機械地分割到3個頻帶之中，得到的3個分量均失去了物理意義。對比結果進一步說明了MMD方法的優(yōu)越性。

(a)IMF1 (b)IMF2 (c)IMF3圖5 仿真信號s(t)的EMD分解結果Fig.5 EMD decomposition results of simulation signal s(t)

(a)BIMF1 (b)BIMF2 (c)BIMF3圖6 仿真信號s(t)的VMD分解結果Fig.6 VMD decomposition results of simulation signal s(t)

3 在變轉速工況齒輪故障特征提取中的應用

齒輪箱包含多種運動零部件，運行時存在多種振動形式，尤其是當工作在非穩(wěn)態(tài)工況時，各個部件的固有頻率振動成分與齒輪、軸承等旋轉部件的嚙合振動、故障沖擊振動等成分在時頻域很可能發(fā)生重疊。傳統(tǒng)的以時間尺度特征為分解依據(jù)的方法無法分離這些成分，而根據(jù)仿真信號的分析MMD方法特別適合處理該類信號。本節(jié)首先利用MMD方法對變轉速工況下的齒輪箱振動信號進行分解，將包含故障成分的分量從其他干擾分量和噪聲中分離出來，然后通過譜分析提取故障特征。

在SpectraQuest公司的旋轉機械故障模擬綜合實驗臺上進行變轉速工況下齒輪故障實驗，振動信號由5個PCB三向加速度傳感器測取，共有S=15個通道，實驗臺以及傳感器安裝位置如圖7a所示，數(shù)據(jù)采集設備為SIEMENS公司的LMS SCADAS Mobile多通道數(shù)據(jù)采集系統(tǒng)，如圖7b所示。故障齒輪所在齒輪箱的結構如圖8所示，轉軸Ⅰ和轉軸Ⅲ分別為輸入軸和輸出軸，齒輪1、2、3、4均為標準直齒輪，齒數(shù)分別為29、100、36和90。實驗前，分別在齒輪1和3上設置了齒面磨損和斷齒故障(實驗1)。為驗證變轉速工況下分析效果，實驗時輸入軸轉速從440 r/min增至3200 r/min，采樣頻率設置為4096 Hz，某通道的時域信號和轉軸I的轉速曲線如圖9所示，其STFT時頻圖見圖10a?？梢钥闯?，分布有包括隨轉速變化的齒輪嚙合振動、故障沖擊振動以及不隨轉速變化的部件固有頻率振動等多種振動成分，且這些成分跨時間尺度特征大，成分之間存在多處交叉，分布復雜。

(a)旋轉機械故障模擬綜合實驗臺

(b)LMS SCADAS Mobile多通道數(shù)據(jù)采集系統(tǒng)圖7 旋轉機械故障實驗臺與數(shù)據(jù)采集系統(tǒng)Fig.7 Rotating machinery fault test bench and dataacquisition system

圖8 齒輪箱結構簡圖Fig.8 Structural diagram of gearbox

圖9 某通道時域信號與轉軸Ⅰ轉速曲線Fig.9 Time domain signal of a channel and rotatingspeed curve of shaft Ⅰ

采用MMD方法進行處理，根據(jù)時頻圖的分布情況，預設分量數(shù)P=5

(a)某通道信號STFT時頻圖 (b)分量1 STFT時頻圖

(c)分量2 STFT時頻圖圖10 變轉速工況齒輪箱振動信號(實驗1)的MMD分解結果Fig.10 MMD decomposition results of variable rotatingspeed condition gearbox vibration signal(experiment 1)

(a)分量1階次譜

(b)分量2階次譜圖11 分量1和分量2的階次譜Fig.11 Order spectrum of the first component andthe second component

(a)分量1包絡階次譜

(b)分量2包絡階次譜圖12 分量1和分量2的Hilbert包絡階次譜Fig.12 The Hilbert envelope order spectrum of the firstcomponent and the second component

對分量1和2進行階次譜分析，結果如圖11所示，兩對齒輪副的理論嚙合階次分別為O1=29和O2=10.44。由圖11a可以看出，在階次29.02、58.04處有清晰的譜線，對應齒輪1和2的一倍和二倍嚙合階次，并且在其兩邊有不對稱的調(diào)制邊頻帶；由圖11b可以看出，在階次10.45、20.89處有清晰的譜線，對應齒輪3和4的一倍和二倍嚙合階次，同樣在其兩旁有不對稱的調(diào)制邊頻帶。為判斷邊頻特征，對分量1和2的Hilbert包絡信號進行階次分析，結果如圖12所示?？梢钥闯?，圖12a中在階次0.9988處有明顯峰值；圖12b中在階次0.2912處有明顯峰值，分別對應了轉軸Ⅰ和轉軸Ⅱ的轉頻階次，這說明齒輪嚙合振動成分被轉頻調(diào)制，符合齒輪出現(xiàn)局部故障時的故障機理，驗證了MMD方法結合階次分析可有效提取變轉速工況下的齒輪故障特征。

(a)IMF1 STFT時頻圖 (b)IMF2 STFT時頻圖

(c)IMF3 STFT時頻圖圖13 變轉速工況齒輪箱振動信號的EMD分解結果Fig.13 EMD decomposition results of variable rotatingspeed condition gearbox vibration signal

為了對比，采用EMD方法對圖9所示信號進行分解，得到的前3個IMF分量如圖13所示，為了直觀比較，同樣采用STFT時頻表示?？梢钥闯觯捎谵D速變化，時間尺度跨度大，導致兩個嚙合振動分量的不同部分被分解到了3個分量之中，失去了原有的物理意義，并且嚙合振動分量與重疊的固有振動成分也未完全分離，同時還存在較為嚴重的背景噪聲。對3個IMF分量及其包絡信號進行階次譜分析，結果如圖14和圖15所示。在圖14中雖然存在嚙合階次，但由于受到其他成分和噪聲的干擾，表現(xiàn)不明顯；而在圖15中與轉頻調(diào)制相關的特征階次也不明顯，尤其是轉軸Ⅱ的轉頻階次。通過對比說明了MMD方法的優(yōu)勢。

(a)IMF1階次譜

(b)IMF2階次譜

(c)IMF3階次譜圖14 前3個IMF分量的階次譜Fig.14 Order spectrum of the first threeIMF components

(a)IMF1包絡次譜

(b)IMF2包絡次譜

(c)IMF3包絡次譜圖15 前3個IMF分量的Hilbert包絡階次譜Fig.15 The Hilbert envelope order spectrum of thefirst three IMF components

將斷齒齒輪3換成同型號正常齒輪進行實驗(實驗2)，實驗工況與分析方式同前。原始信號及分解得到的分量1和2如圖16所示，可以看出,MMD清晰地分解出了兩個嚙合振動分量，分解效果良好。分量1和2及其包絡信號的階次譜如圖17和圖18所示,可以看出，在階次譜中齒輪1和2的嚙合階次兩旁有明顯的邊頻帶，而齒輪3和4的嚙合階次旁邊頻帶不明顯；在包絡階次譜中主要存在轉軸Ⅰ的轉頻及其倍頻階次，而轉軸Ⅱ的轉頻階次不明顯，這符合齒輪箱實際狀況，進一步說明了MMD方法的有效性。

(a)IMF1 STFT時頻圖 (b)IMF2 STFT時頻圖

(c)IMF3 STFT時頻圖圖16 變轉速工況齒輪箱振動信號(實驗2)的MMD分解結果Fig.16 MMD decomposition results of variable rotatingspeed condition gearbox vibration signal(experiment 2)

(a)分量1階次譜

(b)分量2階次譜圖17 分量1和分量2的階次譜Fig.17 Order spectrum of the first component and thesecond component

4 結論

(a)分量1包絡階次譜

(b)分量2包絡階次譜圖18 分量1和分量2的Hilbert包絡階次譜Fig.18 The Hilbert envelope order spectrum of thefirstcomponentand the second component

本文將一種全新的多通道多分量分解(MMD)方法引入變轉速工況齒輪故障特征提取，在采用對數(shù)窗能量準則提高時頻聚集性度量準確性與效率的基礎上，將MMD方法應用于具有復雜時頻特征的齒輪箱振動信號分析，得到了以下結論：

(1)對數(shù)窗能量準則能自適應地確定不同采樣率下的最佳窗長，融入了對數(shù)窗能量準則的MMD方法可減少算法迭代次數(shù)和計算時間，為實際振動信號分析提供參考。

(2)通過仿真信號的分析表明，相比于經(jīng)典的EMD和VMD方法，MMD方法對時間尺度跨度大且時間尺度有重疊的復雜信號具有明顯分解優(yōu)勢。

(3)通過對變轉速工況齒輪箱振動信號的分析表明，MMD方法可從時頻分布復雜的振動信號中有效分解出含有故障特征的分量成分，為齒輪故障特征提取提供了新的有效途徑。

另外值得說明的是，MMD方法是基于多通道數(shù)據(jù)的，但在某些實際工程設備的齒輪箱上安裝多個振動傳感器可能會受到限制。對于此種情況：一是因為信號之間只要存在相位差，MMD方法就會認為其是不同通道的信號，而實際情況下不同傳感器或同一傳感器不同方向測得的信號都存在相位差，因此在安裝受限的條件下，可以考慮將不同傳感器盡可能地緊鄰布置在可安裝的部位；二是雖然實際振動信號可能包含的分量較多，但包含故障特征的分量一般為前幾階分量，即只要分解出前幾個主要分量即可，因此大部分情況下所需的傳感器數(shù)量不多；三是對于極端情況，如果條件允許，可以嘗試在齒輪箱上增加一個放置傳感器的裝置；反之，可以嘗試對采集的單通道振動數(shù)據(jù)人為設置相位差，從而變?yōu)槎嗤ǖ赖臄?shù)據(jù)，這是基于雖然相位不同，但其中的分量成分的頻率特征不會發(fā)生改變，即不會對故障特征提取造成影響，這也是后續(xù)的研究方向之一。