基于量子高斯混合模型的振動信號降噪方法

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

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

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

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

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

1 高斯混合模型的建立

1.1 雙樹復小波包變換

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

(1)

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

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

(2)

(3)

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

g(n)≈h(n-0.5)

(4)

信號重構時，根據下式對樹a和樹b進行聯合重構

(5)

1.2 高斯混合模型

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

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

Y=X+N

(6)

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

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

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

(7)

式中：a1，a2分別為各個單高斯分布的權重;σ1，σ2分別為各個單高斯分布的標準差。

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

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

(8)

1.3 貝葉斯最大后驗估計

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

(9)

根據貝葉斯理論

(10)

則式(9)變換為

(11)

此式等價于

(12)

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

(13)

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

(14)

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

σN=median(Yd)/0.674 5

(15)

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

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

2 量子高斯混合模型

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

2.1 量子疊加態原理

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

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

(16)

式中，系數a和b稱為|0〉和|1〉兩個基態的量子概率幅，其模的平方稱為量子概率，表示對應基態出現的概率值。量子概率幅需滿足歸一化條件，即

(17)

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

2.2 基于量子疊加態的模型參數估計

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

(18)

(19)

(20)

(21)

(22)

將小波系數的量子比特疊加態應用于傳統高斯混合模型得到的小波系數收縮函數式(14)可得

(23)

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

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

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

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

3 仿真信號分析

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

(24)

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

(25)

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

(26)

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

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

4 實測信號分析

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

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

Fig.1 Denoising process for simulated signal by different methods

圖2 行星齒輪箱實驗臺

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

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

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

(a) 時域波形

(b) 包絡譜

(a) 時域波形

(b) 包絡譜

(a) 時域波形

(b) 包絡譜

(a) 時域波形

(b) 包絡譜

5 結 論

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

(a) 時域波形

(b) 包絡譜