楊望燦， 張培林， 陳彥龍， 吳定海， 李海平

(1.91404部隊,秦皇島 066004;2.陸軍工程大學(xué)石家莊校區(qū) 七系，石家莊 050003;3.陸軍特種作戰(zhàn)學(xué)院,桂林 541000；4.陸軍工程大學(xué)石家莊校區(qū) 六系，石家莊 050003)

在機械設(shè)備運行過程中，由于其工作條件和工(作環(huán)境較為惡劣，因此機械設(shè)備部件容易出現(xiàn)磨損和發(fā)生故障[1]。為了及時了解機械設(shè)備的運行狀態(tài)和及早發(fā)現(xiàn)機械設(shè)備的故障情況，振動信號監(jiān)測和分析方法是一種行之有效的技術(shù)手段[2]。但是，由于機械設(shè)備的運行工況較為復(fù)雜，振源較多，采集到的振動信號經(jīng)常被強背景噪聲所淹沒，導(dǎo)致機械設(shè)備故障狀態(tài)特征不明顯，降低了對機械設(shè)備運行狀態(tài)判斷的準確程度。因此，如何對采集的機械振動信號進行降噪處理得到了廣大學(xué)者的熱切關(guān)注[3-5]。

由于小波變換具有多分辨率的時頻分析特性，因此被用于多種類型信號的降噪分析處理中[6-7]。其中，部分學(xué)者通過假設(shè)小波分解系數(shù)符合某一概率統(tǒng)計模型，結(jié)合貝葉斯先驗知識，估計噪聲小波系數(shù)方差，得到小波系數(shù)收縮函數(shù)，從而實現(xiàn)信號的降噪。文獻[8]將雙樹復(fù)小波分解信號得到的高頻系數(shù)和低頻系數(shù)分別采用最大后驗估計算法和廣義形態(tài)濾波進行降噪處理，消除可見光近紅外光譜噪聲。文獻[9]應(yīng)用高斯混合模型(Gaussian Mixture Model，GMM)對小波分解的各類系數(shù)建模，通過噪聲方差估計對信號進行降噪。其中，GMM能夠較好地描述小波系數(shù)的分布情況，降噪效果比較理想。但是傳統(tǒng)的GMM使用全局統(tǒng)一的閾值函數(shù)，算法的局部自適應(yīng)較差，使算法的應(yīng)用受到一定的限制[10]。

自從Eldar等[11]提出量子信號處理(Quantum Signal Processing，QSP)理論后，基于量子理論的信號處理技術(shù)逐漸得到了廣泛地應(yīng)用[12-14]。文獻[15]通過建立量子模型，對聲音信號進行量子傅里葉編碼，相比傳統(tǒng)的聲音信號處理方法，該方法在魯棒性和計算復(fù)雜度方面均有改善。文獻[16]將圖像像素的位置信息和灰度值信息分別用比特量子系統(tǒng)的基態(tài)和量子比特描述，實現(xiàn)了量子圖像的水印加注。受量子衍生信號處理技術(shù)的啟發(fā)，本文提出了一種基于量子高斯混合模型(Quantum Gaussian Mixture Model，QGMM)的振動信號降噪方法，將量子疊加態(tài)理論應(yīng)用于雙樹復(fù)小波包系數(shù)建立的高斯混合模型中，改善高斯混合模型的局部自適應(yīng)性，使小波包系數(shù)自適應(yīng)非線性收縮，提升機械振動信號的降噪處理效果。

1 高斯混合模型的建立

1.1 雙樹復(fù)小波包變換

雙樹復(fù)小波變換是一種有限冗余的小波變換方法，具有近似平移不變性，改善了傳統(tǒng)離散小波變換的Gibbs效應(yīng)和平移敏感性。由于雙樹復(fù)小波變換沒有對高頻部分進一步分解，因此為了提高振動信號的頻率分辨率，本文采用雙樹復(fù)小波包變換(Dual-Tree Complex Wavelet Packet Transform，DTCWPT)對信號的高頻部分和低頻部分同時進行雙樹復(fù)小波變換，用兩棵并行的小波樹實現(xiàn)振動信號的分解。雙樹復(fù)小波包變換如下式所示

(1)

式中：l=(1，2，…，J)為變換尺度因子;J為最大變換尺度。

雙樹復(fù)小波包變換第一層采用傳統(tǒng)的非下抽樣小波變換分解，然后對小波系數(shù)進行奇偶分離，分別采用Mallat算法進行分解和重構(gòu)，雙樹復(fù)小波包分解的小波系數(shù)如下式所示

(2)

(3)

式中:h(n)和g(n)為希爾伯特變換對的濾波器組，滿足半采樣延遲，即

g(n)≈h(n-0.5)

(4)

信號重構(gòu)時，根據(jù)下式對樹a和樹b進行聯(lián)合重構(gòu)

(5)

1.2 高斯混合模型

高斯混合模型是一種應(yīng)用較為廣泛的數(shù)理統(tǒng)計模型，其理論基礎(chǔ)為概率統(tǒng)計，具有較強的靈活性.機械振動信號經(jīng)過雙樹復(fù)小波包變換后，其小波分解系數(shù)呈現(xiàn)稀疏和聚集的分布特性，即在零值附近小波系數(shù)較為稀疏，但幅值較大；在分布兩端小波系數(shù)較為聚集，但幅值較小[17-18]。由于高斯混合模型的分布曲線具有“高峰值”和“長拖尾”的特性，所以可以用高斯混合模型來描述機械振動信號小波系數(shù)的分布規(guī)律。

對含有加性高斯白噪聲的信號y(t)(y(t)=x(t)+n(t))進行雙樹復(fù)小波包變換，根據(jù)小波變換的線性性質(zhì)，信號分解得到的小波系數(shù)滿足以下關(guān)系

Y=X+N

(6)

式中:Y=Yr+iYi，X=Xr+iXi，N=Nr+iNi分別為含噪信號，有用信號和噪聲信號的雙樹復(fù)小波包系數(shù)。

對含噪信號進行雙樹復(fù)小波包分解后，將同一層次、同一子帶內(nèi)復(fù)小波系數(shù)的實部系數(shù)和虛部系數(shù)按照一一間隔的方式排列，即Yd(s)=[Yr(1,s)，Yi(1,s)，Yr(2,s)，Yi(2,s)，…，Yr(nk,s,s)]，其中nk,s為第s層第k個子帶中小波系數(shù)的個數(shù)。對小波系數(shù)排列Yd(s)采用兩狀態(tài)的高斯混合模型建立模型，其概率密度函數(shù)如下式所示

F(Yd)=a1Gaussian(Yd,σ1)+a2Gaussian(Yd,σ2)=

(7)

式中：a1，a2分別為各個單高斯分布的權(quán)重;σ1，σ2分別為各個單高斯分布的標(biāo)準差。

對于高斯混合模型中參數(shù)a1，a2，σ1和σ2，采用最大期望(Expectation Maximization，EM)算法計算確定。EM算法首先計算對數(shù)似然函數(shù)的期望，然后搜索使期望值最大的參數(shù)，得到混合高斯模型參數(shù)的最優(yōu)值，具體步驟參見文獻[19-20]。

對于高斯白噪聲信號，噪聲信號的小波系數(shù)近似服從均值為0，標(biāo)準差為σN的高斯分布，其概率密度函數(shù)如下所示

(8)

1.3 貝葉斯最大后驗估計

對于振動信號的降噪，可以通過從觀測信號y(t)中盡可能精確地估計出真實信號x(t)來完成信號的降噪。本文采用貝葉斯最大后驗估計(Maximum A Posterior，MAP)算法對建立的高斯混合模型進行求解，實現(xiàn)信號的降噪處理。

(9)

根據(jù)貝葉斯理論

(10)

則式(9)變換為

(11)

此式等價于

(12)

令f(X)=lnpX(X)，對式(12)求一階導(dǎo)數(shù)，并使其等于0，得到下式

(13)

將式(7)、式(8)代入式(13)，經(jīng)過推導(dǎo)可得X的估計式，即小波系數(shù)收縮函數(shù)，如下所示

(14)

式中：σN為噪聲信號的標(biāo)準差，其魯棒估計為雙樹復(fù)小波包系數(shù)的絕對值的中值[21]，即

σN=median(Yd)/0.674 5

(15)

根據(jù)建立的混合高斯模型和MAP估計，可以從含噪信號的小波系數(shù)中估計出有用信號的小波系數(shù)，然后對估計的小波系數(shù)進行重構(gòu)，實現(xiàn)信號的降噪。但是從式(14)可知，傳統(tǒng)的高斯混合模型為全局概率模型，其參數(shù)均為全局統(tǒng)一的參數(shù)，不具有局部自適應(yīng)性。當(dāng)信號中信噪比較低時，傳統(tǒng)的高斯混合模型對小波系數(shù)的分布不能達到理想的擬合精度，影響信號的降噪效果。

離合詞“A了個B”與網(wǎng)絡(luò)詞“A了(嘞)個B”同屬概念重組的結(jié)果，形式上相同，兩者都不能接賓語，但如果我們深究入語法功能、語體運用等方面，可發(fā)現(xiàn)它們之間的巨大差別。

2 量子高斯混合模型

傳統(tǒng)的高斯混合模型本質(zhì)上是一種層內(nèi)模型，其僅考慮了同一層次內(nèi)小波系數(shù)的分布情況，但實際上信號的小波系數(shù)還存在著高階相關(guān)性，如尺度間的傳遞性。對于含噪信號的小波系數(shù)，其實質(zhì)上是有用信號小波系數(shù)和噪聲信號小波系數(shù)的疊加，與量子信號處理理論中量子疊加態(tài)原理相似，因此，根據(jù)量子疊加態(tài)理論和小波系數(shù)尺度間的相關(guān)性，提出了一種量子高斯混合模型，使高斯混合模型具有局部自適應(yīng)性，改善信號的降噪效果。

2.1 量子疊加態(tài)原理

在量子理論中，量子比特是描述量子世界的基本單位，量子比特有|0〉和|1〉兩個基本狀態(tài)。根據(jù)量子疊加態(tài)原理，量子比特所表達的狀態(tài)如下所示

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

(16)

式中，系數(shù)a和b稱為|0〉和|1〉兩個基態(tài)的量子概率幅，其模的平方稱為量子概率，表示對應(yīng)基態(tài)出現(xiàn)的概率值。量子概率幅需滿足歸一化條件，即

(17)

含噪信號的小波系數(shù)中包含有用信號和噪聲信號兩種小波系數(shù)，如果用基態(tài)|0〉表示噪聲信號的小波系數(shù);基態(tài)|1〉表示有用信號的小波系數(shù)。這樣量子概率幅a和b就可以表示噪聲信號和有用信號小波系數(shù)的概率，從而可以分析處理二者之間的關(guān)系。

2.2 基于量子疊加態(tài)的模型參數(shù)估計

由于小波系數(shù)具有尺度間的傳遞性，當(dāng)父代系數(shù)的模較大時，那么其子代系數(shù)的模也大，反之，當(dāng)父代系數(shù)的模較小時，其子代系數(shù)的模也較小。所以將父代和子代的小波系數(shù)取模相乘得

(18)

(19)

(20)

(21)

(22)

將小波系數(shù)的量子比特疊加態(tài)應(yīng)用于傳統(tǒng)高斯混合模型得到的小波系數(shù)收縮函數(shù)式(14)可得

(23)

從式(23)可知，量子高斯混合模型通過小波系數(shù)的量子疊加態(tài)特性，根據(jù)信號自身的特點，對傳統(tǒng)高斯混合模型中參數(shù)進行了局部自適應(yīng)調(diào)整。當(dāng)尺度s+1下第j個小波系數(shù)中出現(xiàn)有用信號的概率大時，則自適應(yīng)地增大該位置有用信號的標(biāo)準差，反之，當(dāng)尺度s+1下第j個小波系數(shù)中出現(xiàn)噪聲信號的概率大時，則自適應(yīng)地增大該位置噪聲信號的標(biāo)準差，提高高斯混合模型對小波系數(shù)的分布的擬合精度，改善信號的降噪效果。

2.3 基于量子高斯混合模型的降噪流程

基于量子高斯混合模型的信號降噪流程步驟如表1所示。

表1 量子高斯混合模型降噪算法

3 仿真信號分析

為了模擬行星齒輪箱齒輪局部故障，構(gòu)造了如下的仿真信號

(24)

仿真信號的采樣頻率為4 096 Hz，采樣點數(shù)為2 048個。由式(24)中可知，仿真信號由一個正弦諧波分量，一個調(diào)幅分量和一個調(diào)頻分量線性疊加而成。在仿真信號中加入高斯白噪聲，雙樹復(fù)小波包采用4層分解，采用信噪比σSNR作為評價指標(biāo)，信噪比計算公式如下

(25)

在實驗過程中，引入傳統(tǒng)的高斯混合模型、硬閾值和軟閾值的降噪方法作為對比實驗，傳統(tǒng)的硬閾值、軟閾值的閾值計算按下式計算

(26)

式中:nk,s為第s層第k個子帶中小波系數(shù)的總個數(shù)。硬閾值降噪方法將絕對值低于閾值的小波系數(shù)置0，保留其他小波系數(shù)進行重構(gòu)信號；軟閾值降噪方法將絕對值低于閾值的小波系數(shù)置0，將絕對值高于閾值的小波系數(shù)減去閾值后進行重構(gòu)信號。

不同方法對仿真信號的實驗結(jié)果如圖1所示。從圖1可知，對于仿真信號，量子高斯混合模型去噪效果優(yōu)于其他幾種方法，降噪后的信噪比最高。硬閾值和軟閾值降噪方法，在去除噪聲的同時，也消除了部分原始的真實信號，導(dǎo)致信號部分失真，信噪比提升有限。傳統(tǒng)的高斯混合模型采用全局統(tǒng)一的模型參數(shù)，導(dǎo)致噪聲殘留較多，降噪后的信噪比也不高。而量子高斯模型將小波系數(shù)尺度間的傳遞性用量子比特的疊加態(tài)特性描述，自適應(yīng)調(diào)整模型參數(shù)，在保留原始波形特征的基礎(chǔ)上，最大限度的消除噪聲，取得更優(yōu)的降噪效果。

4 實測信號分析

實測信號來自如圖2所示的行星齒輪箱傳動系統(tǒng)實驗臺，加速度傳感器安裝在行星齒輪箱的箱體上，采集實驗過程中的振動信號，采樣頻率為20 kHz。該實驗通過機械加工的手段模擬了太陽輪齒輪的單個齒齒面的輕微磨損故障。實驗過程中，電機轉(zhuǎn)速為1 200 r/min，磁粉制動器提供的負載為1.2 N·m，采樣時長為2 s。行星齒輪箱的結(jié)構(gòu)參數(shù)如表2所示。根據(jù)理論計算，太陽輪齒輪故障的特征頻率為55.5 Hz。

圖1 不同方法對仿真信號的降噪處理

Fig.1 Denoising process for simulated signal by different methods

圖2 行星齒輪箱實驗臺

齒輪太陽輪行星輪(個數(shù))齒圈齒數(shù)1364(3)146

圖3為行星齒輪箱實驗臺采集的太陽輪齒輪故障的原始振動信號的時域波形和包絡(luò)頻譜圖，為了便于圖形展示，時域波形的時間取0～0.5 s，頻譜圖中頻率取0～500 Hz。從圖3可知，由于故障較為輕微，故障齒輪所產(chǎn)生的沖擊信號不明顯，基本淹沒在背景噪聲中，在包絡(luò)譜中也看不到太陽輪故障的特征頻率，無法判斷行星齒輪箱是否發(fā)生故障。

圖4～圖7為分別采用不同方法對實測太陽輪故障振動信號的降噪結(jié)果。實驗過程中，雙樹復(fù)小波包分解都為4層分解。從圖4可知,采用量子高斯混合模型降噪后，振動信號的時域波形出現(xiàn)了沖擊特征，對其進行包絡(luò)解調(diào)，從包絡(luò)譜圖中可以明顯地看到太陽輪故障的特征頻率55.5 Hz及其2倍頻、3倍頻和4倍頻，符合太陽輪故障的頻譜特點，所以可以判斷行星齒輪箱中太陽輪輪齒發(fā)生了故障。圖5為傳統(tǒng)高斯混合模型的降噪結(jié)果，雖然從包絡(luò)譜圖中能看到太陽輪故障的特征頻率55.5 Hz，但是仍有大量的干擾噪聲，特征頻率淹沒在其他干擾頻率中，不能準確判斷行星齒輪箱中太陽輪是否存在故障。圖6和圖7為采用小波硬閾值和小波軟閾值的降噪結(jié)果，噪聲被大量去除，但是由于太陽輪齒輪磨損輕微，太陽輪故障特征信號也被當(dāng)做噪聲信號被去除或是削弱，所以包絡(luò)譜圖中也無法得到太陽輪故障的特征頻率。

(a) 時域波形

(b) 包絡(luò)譜

(a) 時域波形

(b) 包絡(luò)譜

(a) 時域波形

(b) 包絡(luò)譜

(a) 時域波形

(b) 包絡(luò)譜

5 結(jié) 論

本文提出了一種基于量子高斯混合模型的機械振動信號降噪方法。該方法在傳統(tǒng)高斯混合模型的基礎(chǔ)上，根據(jù)量子信號處理理論，將信號小波系數(shù)相鄰尺度間的傳遞特性轉(zhuǎn)化為量子比特疊加態(tài)，使高斯混合模型中的參數(shù)根據(jù)信號特征自適應(yīng)地調(diào)整大小，增強了高斯混合模型的局部自適應(yīng)性。仿真信號分析表明，相比于其他方法，量子高斯混合模型降噪方法能夠有效提高信號的信噪比。實測行星齒輪箱振動信號結(jié)果表明，本文所提方法能夠抑制信號中的噪聲，保留太陽輪齒輪輕微故障的特征信號，準確提取其故障特征頻率，有效地判斷行星齒輪箱的故障狀態(tài)。

(a) 時域波形

(b) 包絡(luò)譜