基于量子Hadamard變換的滾動軸承振動信號分析方法

王懷光，陳彥龍，楊望燦，王 強

(1. 軍械工程學(xué)院 七系，石家莊 050003; 2. 陸軍特種作戰(zhàn)學(xué)院，廣西 桂林 541000)

滾動軸承是機械設(shè)備的重要組成部件，是機械設(shè)備正常運行的保證[1]。滾動軸承振動信號中包含了豐富的運行狀態(tài)信息，因此在滾動軸承故障診斷與狀態(tài)監(jiān)測過程中，發(fā)揮著重要作用[2]。但現(xiàn)場采集的振動信號通常受到較強的噪聲干擾，故障信息往往淹沒在強烈的背景噪聲當(dāng)中，難以提取反映滾動軸承運行的關(guān)鍵特征，使得振動信號的利用率極大降低，在很大程度上制約了滾動軸承效能的有效發(fā)揮。對于這一問題，振動信號降噪[3]是一種有效的解決方式。

小波分析[4]與形態(tài)學(xué)分析[5]是振動信號降噪中兩種常用的方法。在不同應(yīng)用條件下都取得了良好的降噪效果。

小波變換是目前運用最廣泛的降噪方法。以小波分析理論為基礎(chǔ)，文獻(xiàn)[6]采用了雙樹復(fù)小波對信號進(jìn)行分析，通過塊閾值的方法，結(jié)合領(lǐng)域與空域，有效降低背景噪聲。文獻(xiàn)[7]對不同分解層次的系數(shù)求模，然后對系數(shù)采用非線性時間序列分析，消除隨機噪聲的干擾，最后再對系數(shù)采取軟閾值去除直流信息，順利提取軸承中的微弱故障成分，等等。分析小波降噪的原理，其充分利用了有用信號與噪聲信號在小波域變換系數(shù)的特性，然后利用軟閾值或者硬閾值的方式將部分小波系數(shù)濾除，從而保留下信號的有用信息。但是這種方式建立在信號整體特征的基礎(chǔ)上提取有用信息，信號的降噪處理更多的是對信號進(jìn)行一體化處理，因此在降噪過程中容易忽略每個采樣點中所包含的細(xì)節(jié)信息，容易造成噪聲信號與非噪聲信號混雜，有時會抑制非噪聲信號，嚴(yán)重影響了信號的降噪效果。

形態(tài)學(xué)分析是一種非線性分析工具，近年來在信號降噪領(lǐng)域得到廣泛關(guān)注。以形態(tài)學(xué)分析為基礎(chǔ)，文獻(xiàn)[8]提出多尺度平均組合形態(tài)濾波方法，對原始振動信號進(jìn)行降噪處理，實現(xiàn)對模態(tài)混疊的抑制。文獻(xiàn)[9]研究了基于廣義形態(tài)分量分析的降噪方法，首先把一維觀測信號轉(zhuǎn)換為多維虛擬觀測信號，再借助廣義形態(tài)分量分析方法，通過對觀測信號的有效分離實現(xiàn)降噪。基于形態(tài)學(xué)分析的振動信號降噪方法主要利用腐蝕、膨脹、形態(tài)開運算、形態(tài)閉運算四種基本算子，并利用信號與結(jié)構(gòu)元素的匹配相似度，實現(xiàn)有用信號與噪聲信號的分離，達(dá)到降噪的目的。但是這種分析方式在信號處理過程中并沒有充分考慮到信號的隨機性對降噪效果的影響，并且對于某一特定振動信號而言，結(jié)構(gòu)元素的選取對其信號處理過程的影響較大，導(dǎo)致結(jié)構(gòu)元素的適應(yīng)性弱，在一定程度上影響了形態(tài)學(xué)分析方法的降噪效果。

為進(jìn)一步改進(jìn)信號的降噪效果，本文改變傳統(tǒng)的信號降噪理念，結(jié)合量子理論[10]，提出基于量子Hadamard變換[11]的降噪方法。該方法充分發(fā)掘量子理論在振動信號處理中的潛能，通過量子化的方法將振動信號單個采樣點進(jìn)行表達(dá)，賦予了單個采樣點信號量子特性。通過量子化，單個采樣點的幅值特性被弱化，而是轉(zhuǎn)化成具有特定分布特性的概率分布模型。通過分析處理單個采樣點的概率幅，判斷信號是否為噪聲信號，并對故障信號進(jìn)行有效凸顯。與基于小波變換的信號降噪方式相比，本文算法通過對單個采樣點信號進(jìn)行降噪處理，能夠克服信號一體化處理帶來的不足，避免了閾值處理過程中有用信號與噪聲信號的混雜。與基于形態(tài)學(xué)分析的信號降噪方法相比，本文算法對單個采樣點的降噪過程合理利用了采樣點的概率分布特性，避免了降噪過程中信號隨機性、結(jié)構(gòu)元素等因素對信號降噪過程的影響，對分布參數(shù)的處理方法具有較強的適應(yīng)性。實驗結(jié)果表明，該方法能夠有效改善信號的降噪效果，同時機械振動信號的故障特征得到凸顯。

1 量子Hadamard變換

Hadamard變換(Hadamard Transform，HT)是量子理論框架下一種普遍使用的變換方式。HT與量子比特(Quantum Bit，QB)息息相關(guān)，在一個量子系統(tǒng)中，假設(shè)其QB長度為n，則該系統(tǒng)的態(tài)矢總數(shù)為N=2n。構(gòu)造N×N的酉矩陣H(H=HT，HH=I)，將系統(tǒng)在二維空間下展開，那么HT通?？杀硎緸閇12-13]

(1)

量子比特表達(dá)式為

|Ψ>=a|0>+b|1>

(2)

對量子比特|Ψ>=a|0>+b|1>應(yīng)用HT可得

(3)

2 基于HT的振動信號量子化

2.1 基于HT的量子概率幅的參數(shù)范圍

根據(jù)量子理論的相關(guān)約束條件，設(shè)定a，b∈[0, 1]，振動信號的狀態(tài)包括噪聲信號以及非噪聲信號(故障信號)，為便于量子理論的應(yīng)用，將故障信號設(shè)為基態(tài)|0>，將噪聲信號設(shè)為基態(tài)|1>。對于實際采集的振動信號，其信號狀態(tài)通常在故障信號以及噪聲信號之間變化，故a，b可設(shè)定為

a,b∈(0, 1)

(4)

根據(jù)量子理論的基本知識，結(jié)合式(3)中可知，狀態(tài)|0>出現(xiàn)的概率為

(5)

同理可得，狀態(tài)|1>出現(xiàn)的概率為

(6)

式(5)和式(6)稱為Hadamard量子概率(Hadamard Quantum Probability，HQP)。

從式(5)和式(6)中可以發(fā)現(xiàn)，對于不同的基本狀態(tài)而言，其HQP取值取決于a×b，由于a，b為量子比特|Ψ>=a|0>+b|1>的概率幅，根據(jù)文獻(xiàn)[14]，振動信號QB概率幅可以表示為

(7)

(8)

結(jié)合前文可知，p滿足p∈(0,1)，故

(9)

對其進(jìn)行一階求導(dǎo)可得

(10)

圖1 HQP變化Fig.1 Change of HQP

為進(jìn)一步簡化過程，設(shè)定式(5)和式(6)在其可行域內(nèi)單調(diào)增或者單調(diào)減，當(dāng)p∈(0,0.5]時，a，b可以表示為

(11)

(12)

2.2 基于HT的振動信號量子化

結(jié)合前文可知，對于采集的振動信號，需要進(jìn)行歸一化處理

(13)

式中：s(k)為信號采樣點，歸一化處理后信號幅值z(k)∈[0,1]，量子化后的概率幅由z(k)折半生成，定義為

(14)

由于每一個信號點狀態(tài)介于故障信號與噪聲信號之間，對于QB，其概率幅可表示為

(15)

(16)

ε>0，為取值趨近于0的正值。進(jìn)而進(jìn)行量子化處理，可以表示為

(17)

其中，

(18)

(19)

則根據(jù)陳彥龍等的研究，經(jīng)過HT后，QB可以表示為

(20)

分析得知，經(jīng)過HT后，式(5)、式(6)可表示為

(21)

(22)

從式中可以看出，基本狀態(tài)|0>出現(xiàn)的概率與abs(s(k))呈正比，基本狀態(tài)|1>出現(xiàn)的概率與abs(s(k)呈反比，因此采樣點為故障信號的概率與abs(s(k))呈正比，為噪聲信號的概率與abs(s(k)呈反比。

3 基于HT的降噪方法

3.1 可行性分析

3.1.1 衡量算子

對于滾動軸承的振動信號而言，若其不含有噪聲，則信號間的關(guān)聯(lián)性較強，而噪聲信號會破壞這種關(guān)聯(lián)性，因此本文借助量子理論提出衡量算子(Measurement Operator，MO)的概念。圖2為振動信號1×3的鄰近信號窗口。

sn(k-1)sn(k)sn(k+1)圖2 1×3鄰域位置關(guān)系Fig.2 1×3 neighborhood

衡量算子可以表達(dá)為

mo(k)=a(k-1)×a(k)×a(k+1)

(23)

式中：mo(k)為衡量算子；a(k)為振動信號量子化后的概率幅。

mo(k)可具體表達(dá)為

(24)

經(jīng)過上述分析得知，mo(k)的大小不僅取決于k點sn(k)的大小，還取決于左右相鄰點sn(k-1)，sn(k+1)的大小，因此在sn(k-1)，sn(k)，sn(k+1)均為較大值時，mo(k)才能取得較大值，結(jié)合振動信號的特點，通常故障信號在峰值附近sn能夠連續(xù)取得較大值，因此故障信號在峰值附近取得較大值的概率較大。這種變化使得故障信號能夠得到有效突出，有利于信號的降噪處理。

3.1.2 閾值確定

滾動軸承振動信號本身具有隨機性，使得正向脈沖信號P和負(fù)向脈沖信號N的產(chǎn)生具有隨機性，針對這種情況，可以從每一個采樣點的自身特性出發(fā)，設(shè)定一個不同的閾值T(n)，并以閾值作為確定不同點處理方式的依據(jù)。

閾值的選取，對降噪效果有著至關(guān)重要的影響，當(dāng)閾值選取過大時，可能會抑制信號中的非噪聲信號，而當(dāng)閾值較小時，噪聲信號可能無法得到有效抑制，無法達(dá)到理想的降噪效果?？紤]到脈沖信號P和負(fù)向脈沖信號N的隨機性，本文提出與單個采樣點信號特征相適應(yīng)的閾值確定方式。

單個采樣點閾值的確定通過中值濾波器(Median Filter，MF)實現(xiàn)，中值濾波器是一種平滑濾波器，在振動信號去噪中應(yīng)用較多，通常情況下，無論是噪聲還是故障信號都會導(dǎo)致振動信號發(fā)生突變，針對這一特點，本文設(shè)計合適的MF來選取合適的閾值，達(dá)到故障信號突出，噪聲信號消除的目的，其過程可以表示為

T(k)=med(mo(sn(k-3)),,mo(sn(k+3)))

(25)

從式(25)中可以看出，采用MF濾波器，有效保證了閾值T(k)低于峰值點mo值。

3.2 降噪方法

3.2.1 正向故障脈沖信號處理

針對正向信號P，假設(shè)P在k點附近表現(xiàn)為故障，則故障特征越明顯，a×b計算結(jié)果值越大，兩種基本的狀態(tài)在經(jīng)過HT后概率變化情況為：式(5)取值變大，式(6)取值變小。根據(jù)分析得知，要達(dá)到故障突出表達(dá)的目的，對于mo(k)>T(k)的情況，該點信號應(yīng)加強，對于mo(k)≤T(k)的情況，該點信號應(yīng)減弱。因此可以利用式(5)、式(6)中的HQP實現(xiàn)信號降噪。在mo(k)>T(k)的條件下，該采樣點很大概率上為故障信號，進(jìn)而利用|0>的HQP增強信號

(26)

在mo(k)的HQP降噪

(27)

由式(26)、式(27)的計算結(jié)果可知，當(dāng)mo(k)>T(k)時，信號加強的程度與信號幅值s(k)成正比；當(dāng)mo(k)≤T(k)時，減弱程度與s(k)成反比。

由式(26)、式(27)分析得到

sn(k)+0.5≤0.5+sn(k)+a(k)×b(k)≤sn(k)+1

(28)

sn(k)-0.5≤-0.5+sn(k)+a(k)×b(k)≤sn(k)

(29)

從式(28)、式(29)中可以看出，對于不同采樣點信號，若mo(k)>T(k)，sn(k)加強，加強的范圍為0.5～1.0，且當(dāng)sn(k)值較大時，sn(k)增加較多，最多增加1；若mo(k)≤T(k)，sn(k)數(shù)值減弱，減弱的范圍為0～0.5，且當(dāng)sn(k)值越小時，sn(k)減少的越多，最多減少0.5。

3.2.2 負(fù)向故障脈沖信號處理

針對負(fù)向信號N，假設(shè)N在k點附近表現(xiàn)為故障，在mo(k)>T(k)的條件下，該采樣點很大概率上為故障信號，進(jìn)而利用|0>的HQP進(jìn)行負(fù)向增強

(30)

在mo(k)的HQP降噪

(31)

由式(30)、式(31)的計算結(jié)果可知，當(dāng)mo(k)>T(k)時，信號減少的程度(故障信號加強)與信號絕對值成正比；當(dāng)mo(k)≤T(k)時，信號幅值增加程度(故障信號減弱)與信號絕對值成反比。

由式(30)、式(31)可得

-sn(k)-1≤-0.5-sn(k)-a(k)×b(k)≤ -sn(k)-0.5

(32)

-sn(k)-0.5≤-0.5-sn(k)+a(k)×b(k)≤ -sn(k)

(33)

從式(32)、式(33)中可以看出，對于不同采樣點信號，若mo(k)>T(k)，-sn(k)數(shù)值得到減少，減少的范圍為0.5～1.0，且當(dāng)-sn(k)值較小時，-sn(k)減少較多，最多減少1；若mo(k)≤T(k)，-sn(k)數(shù)值增加，增加的范圍為0～0.5，且當(dāng)sn(k)值越大時，-sn(k)增加的越多，最多增加0.5。

4 降噪算法步驟

根據(jù)上述分析，通過HT得到了滾動軸承振動信號時域處理算法，該算法對振動信號進(jìn)行量子化處理，進(jìn)而實現(xiàn)噪聲的消除與故障信號的增強，并且在整個過程中參數(shù)變化具有適應(yīng)性?；诹孔親T的分析方法(Analysis Method Based on Quantum Hadamard Transform，AMQHT)步驟可以描述為：

步驟1基于HT，利用式(13)、式(14)實現(xiàn)振動信號的量子化，處理方法如式(20)所示；

步驟2利用HT，對每一個信號點進(jìn)行處理，根據(jù)式(21)和式(22)計算出不同點對應(yīng)的HQP；

步驟3利用式(23)，計算不同信號點的MO值；

步驟4利用式(25)確定不同信號點的T(k)大?。?/p>

步驟5若采樣點滿足s(k)≥0，設(shè)定為P(正向)，利用式(26)和式(27)進(jìn)行降噪處理。

步驟6若采樣點滿足s(k)<0，設(shè)定為N(負(fù)向)利用式(30)和式(31)進(jìn)行降噪處理。

另一個方面，上述步驟5、步驟6，根據(jù)分析可以看出，出現(xiàn)在s(k)≥0區(qū)域的波谷，處理完之后同樣為波谷；出現(xiàn)在s(k)≤0區(qū)域的波峰，處理完之后同樣為波峰。因此，經(jīng)過本文算法進(jìn)行處理后，波峰波谷區(qū)域在實質(zhì)上并不會改變。

5 實驗結(jié)果與分析

5.1 仿真信號實驗

為驗證本文算法的有效性，設(shè)計仿真信號，進(jìn)行信號的降噪實驗。

仿真信號采樣頻率為960 Hz，采樣時間為1 s。仿真信號中成分主要包括兩部分

y=y1+y2

(34)

式中：y為仿真信號；y1為本底信號，其表達(dá)式為

y1(t)=sin(2π×100×t)+sin(2π×200×t)

(35)

y2仿真噪聲信號，實驗中，為分析不同類噪聲對本文算法降噪效果的影響，分別設(shè)置了不同種噪聲信號，包括白噪聲、有色噪聲、脈沖干擾。在滾動軸承振動信號采集過程中，真實信號通常會淹沒在背景噪聲之中，因此設(shè)置白噪聲、有色噪聲標(biāo)準(zhǔn)差為4，有色噪聲表達(dá)式可以表示為

e(k)=x(k)+0.5×x(k-1)

(36)

式中：x(k)為白噪聲信號的第k個點。脈沖干擾為周期性的沖擊信號，單個脈沖信號可以表示為

p(t)=A0e-40×t·cos(2π×150×t)

(37)

式中：A0為幅值，仿真實驗中設(shè)置A0=4，載波信號的頻率為150 Hz，沖擊頻率為10 Hz。仿真信號時域如圖3所示(圖中時間為0.5 s)，頻譜如圖4所示。從時域圖中可以看出，本底信號完全淹沒在噪聲信號中，而在頻域內(nèi)，本底信號特征頻率也受到嚴(yán)重污染。

圖3 仿真信號時域圖Fig.3 The simulated signals in time domain

圖4 仿真信號頻域圖Fig.4 The simulated signals in frequency domain

針對上述仿真信號，利用本文所提的AMQHT算法進(jìn)行降噪，降噪信號的頻譜如圖5所示，從圖中可以看出，對于白噪聲和有色噪聲，經(jīng)過AMQHT降噪后，本底信號的特征頻率得到了充分保留，同時噪聲信號在頻域內(nèi)受到了有效抑制，這就使得頻域內(nèi)特征頻率得到了有效凸顯。實驗結(jié)果表明，當(dāng)原始信號淹沒在白噪聲、有色噪聲信號中時，本文算法能夠降低噪聲的干擾，實現(xiàn)有用信號的凸顯。而對于脈沖干擾而言，本文算法很難將其與本底信號相區(qū)分，甚至將脈沖干擾當(dāng)做有用信號進(jìn)行凸顯，因此在圖5(d)中，本底信號、脈沖載波信號、脈沖頻率均得到保留。實驗結(jié)果表明，本文算法對脈沖干擾的抑制效果不足。

圖5 降噪效果Fig.5 The de-noising effect

5.2 實測信號實驗

為進(jìn)一步驗證本文算法的有效性，利用實測信號進(jìn)行降噪實驗。故障信號來自于某型滾動軸承變速箱。在變速箱的軸承上設(shè)置內(nèi)圈故障。在軸承內(nèi)圈上加工出1 mm×0.2 mm(長×深)的劃痕。采集加速度振動信號，實測轉(zhuǎn)動速度1 748 r/min(29.1 Hz)，傳感器安裝在對應(yīng)軸承位置的箱蓋上方。

滾動軸承內(nèi)圈故障的特征頻率計算公式為

(38)

式中：n為軸承轉(zhuǎn)速，r/min； 轉(zhuǎn)頻fr與轉(zhuǎn)速n的關(guān)系為fr=n/60；D為軸承中滾動體的直徑；D為軸承的節(jié)徑；β為軸承的接觸角；z為滾動體的個數(shù)，經(jīng)過計算，內(nèi)圈故障理論頻率應(yīng)為f=157.7 Hz。

振動信號的采樣頻率fs為12 kHz，采樣時間0.6 s，圖6顯示了變速箱軸承內(nèi)圈發(fā)生故障時的振動信號時域波形及其頻譜。時域波形中雖然出現(xiàn)了沖擊的信號，但是振動信號的頻譜卻難以觀測到故障特征頻率f。

圖6 滾動軸承內(nèi)圈故障的振動信號Fig.6 Vibration signals from rolling bearing with fault in inner ring

為驗證本文所提算法的效果，利用軸承內(nèi)圈故障信號進(jìn)行降噪實驗，降噪效果如圖7所示。

對降噪后的軸承內(nèi)圈故障信號進(jìn)行傅里葉變換，頻域特點如圖7(b)所示。從圖中可以看出，利用AMQHT降噪后，故障信號頻率f和及其二倍頻2f得到有效表示。

圖7 AMQHT處理結(jié)果Fig.7 Result of AMQHT

為進(jìn)一步驗證本文算法的降噪效果，利用不同算法對上述內(nèi)圈故障信號進(jìn)行降噪處理，對比不同算法的降噪效果，對比的算法為數(shù)學(xué)形態(tài)學(xué)(Analysis Method Based on Mathematical Morphology，AMMM)。降噪后時域波形與頻譜，如圖8所示。

圖8 AMMM處理結(jié)果Fig.8 Result of AMMM

為量化降噪算法的降噪效果，引入三個振動信號降噪的量化指標(biāo)

(1) 降噪指標(biāo)

(39)

(2) 增強指標(biāo)

(40)

(3) 頻率指標(biāo)

頻率指標(biāo)表示故障特征頻率f的幅度。

不同算法的降噪效果如表1所示，可以看出，在頻率指標(biāo)上，AMQHT稍差于AMMM，但從降噪指標(biāo)及增強指標(biāo)來看，AMQHT要好于AMMM。分析其原因，AMMM對故障信號的增強效果要稍好于AMQHT，但由于AMMM對信號進(jìn)行處理過程中，很難實現(xiàn)故障信號與非故障信號的有效分離，因此容易將非故障信號進(jìn)行增強，表現(xiàn)在降噪指標(biāo)與增強指標(biāo)上，AMMM要差于AMQHT。綜合分析，AMQHT降噪效果要好于AMMM，實驗結(jié)果充分驗證本文算法的有效性。

表1 降噪效果對比Tab.1 Comparison of de-noising effect

6 結(jié) 論

針對滾動軸承信號的降噪問題，引入量子Hadamard變換，建立起一種用于滾動軸承故障狀態(tài)下的振動信號分析方法。該方法以量子理論為基礎(chǔ)，深入考慮每一個采樣點中噪聲和故障信息的變化。將每一個采樣點進(jìn)行量子化后，獨立分析每一個采樣點的信息，有效克服了現(xiàn)有降噪算法對細(xì)節(jié)信息考慮不充分的問題。最后，將本算法應(yīng)用于實測軸承內(nèi)圈故障振動信號的降噪，驗證了本文算法的有效性。