基于Laplace先驗和稀疏塊相關性的旋轉機械振動信號貝葉斯壓縮重構

馬云飛， 白華軍， 溫亮， 郭馳名， 賈希勝

(1.武警士官學校 軍械系， 浙江 杭州 310023； 2.陸軍工程大學石家莊校區 裝備指揮與管理系， 河北 石家莊 050003)

0 引言

為了實時監控裝備運行狀態，利用分布式無線傳感技術實現機械信號傳輸、監測，具有一定的應用價值和發展潛力。相對于有線的故障預測與健康管理系統，無線監測系統能夠大大增加靈活性、可維護性和可擴展性[1]。而受無線傳輸帶寬限制，想要實時監控必須對信號進行壓縮傳輸。此前，有學者對結構健康監測領域的振動信號壓縮進行了研究[2-3]，但在該領域中進行振動信號監測的采樣頻率較低，信號壓縮實現相對較為容易。例如在橋梁結構監測中僅需要240 Hz的采樣頻率[4]，而要想實現機械振動信號監測，則至少需要5～20 kHz的采樣頻率[5]。由于機械振動信號具有非平穩、非線性等特征，且設備運轉過程中會產生大量噪音，給信號壓縮帶來極大困難。

壓縮感知(CS)理論[6-8]近年來引起學者廣泛關注，其基本原理是利用與壓縮信號線性無關的觀測矩陣對稀疏信號進行采樣，降低采樣信號維數，并在接收端利用重構算法估計原始信號。目前，基于CS的重構算法可以歸為三類：凸優化算法[9]、貪婪算法[10]、貝葉斯CS(BCS)算法[11]。BCS算法將貝葉斯估計原理與CS結合，利用稀疏貝葉斯模型中的相關向量機(RVM)[12]學習估計原始信號的最大后驗概率。BCS算法的優點是充分考慮了信號傳輸中產生的噪音，且并不要求原始信號符合稀疏的性質。此外，Babacan等[13]在BCS框架下，利用Laplace先驗代替高斯先驗約束未知信號，該先驗能夠在對數凹時最大程度地消除局部極小值，從而優選接近于零的信號系數。Ji等[14]將單任務BCS擴展為多任務BCS，通過挖掘多任務之間稀疏信號的統計相關性來實現多個稀疏信號聯合重構。

對機械振動信號CS，文獻[15]提出對變分模態分解分量進行CS，并通過循環模式矩陣定義新的觀測矩陣；文獻[16]將以往頻域稀疏化方法改為相空間稀疏化，并利用主成分分析提取特征，更有效促進信號稀疏性；文獻[17]利用改進m序列原理設計觀測矩陣電路實現，并通過滾珠絲杠動態測試中的振動信號驗證了系統可行性。盡管上述研究取得了較好的效果，但是已有文獻缺乏機械振動信號的BCS相關研究，且沒有考慮到機械信號的內部相關性。

本文提出一種基于Laplace先驗和稀疏塊相關性的BCS(Lap-CBCS)算法。由于Laplace先驗具有更好的稀疏促進作用，特別適合于機械振動信號這一類復雜信號。建立基于該先驗的多任務BCS框架，并通過對模型中各信號塊共享的超參數進行估計，實現信號重構。改進算法充分利用Laplace先驗的稀疏促進作用和類周期信號的內部相關性，以提高重構的準確性，為機械設備無線狀態監測提供技術支撐。

1 CS與貝葉斯重構

1.1 CS原理

CS[8]利用一個低維信號近似表達原始信號。假設原始信號為x(x∈RN，N為原始信號長度)，稀疏變換為Ψ，觀測矩陣Φ∈RM×N(M?N，M為觀測矩陣行數或稱觀測次數)。則信號x在該矩陣上的線性投影y∈RM可以作為壓縮信號，

y=Φx=ΦΨθ=Θθ，

(1)

式中：θ表示稀疏信號；Θ為傳感矩陣，Θ=ΦΨ.

由此，給出壓縮率定義：

(2)

由(2)式可知，原始信號越長，觀測次數越小，則壓縮率越高，也越難實現精確重構。在y和Θ已知情況下，求解θ，可根據信號的稀疏性質，將信號重構問題轉化為0-范數最優化問題，即

min ‖θ‖0，
s.t.y=Θθ=ΦΨx.

(3)

該問題可利用基追蹤(BP)算法[8]、正交匹配追蹤(OMP)算法[10]求解。

1.2 貝葉斯重構

1.2.1 問題轉化

BCS算法[11]從概率角度出發，求解原始信號的最大后驗估計值，以恢復信號。令θs表示保留θ的M個最大分量、剩余分量置0的向量，θe表示保留剩余N-M個分量、其余分量置0的向量，則y可拆分成

y=Θθ=Θθs+Θθe=Θθs+ne，

(4)

式中：ne表示信號次要信息。

信號y傳輸過程中，可能存在噪音nm，則y表示為

y=Θθs+ne+nm=Θθs+n，

(5)

式中：n表示方差為σ2的零均值高斯噪聲。

1.2.2 構造先驗分布

構造y關于參數θs和σ2的先驗密度函數：

(6)

RVM[12]是一個通用的貝葉斯框架，該框架使用參數中的線性模型來獲得回歸和分類任務的稀疏解。傳統模型采用分層先驗模型，并假設θs的先驗分布為高斯分布。與該分布相比，Laplace先驗通過在軸上更多地使用后驗分布來嚴格約束稀疏性，因此信號系數更接近于0，本文假設其服從Laplace分布。

1.2.3 參數迭代求解

在RVM框架下如果y和超參數的先驗分布已知，可以得到θs的均值和方差表達式。進一步，令超參數的邊緣分布最大，利用極大似然估計得到超參數更新表達式。通過反復迭代可估算出所有參數值，最后用θs的均值估計表示θs.

2 Laplace先驗和稀疏塊相關的BCS

2.1 基于Laplace分布的貝葉斯先驗模型

在BCS領域應用較多的是估計值服從高斯先驗的模型。研究表明，在實數域應用高斯分布和指數分布共軛的Laplace先驗模型比高斯先驗具有更好的稀疏促進作用[13]。由于機械振動信號也是實數信號，可以用Laplace先驗約束。假設振動信號θs服從Laplace先驗，參數集γ=[γ1，γ2，…，γj，…，γN]。如圖1所示，Laplace分布可表示為高斯分布和指數分布的混合形式，結合對超參數γj的指數約束，得到θs關于λ的Laplace分布表達式為

圖1 Laplace貝葉斯分級先驗模型Fig.1 Laplace priors-based architectural model

(7)

式中：λ為Laplace模型中描述γj的超參數；p(θs|γ)表示θs關于γ的先驗分布；p(γ|λ)表示γ關于λ的先驗分布；p(θs|λ)表示分層先驗分布。

有關噪聲n的超參數α0=1/σ2仍服從Gamma分布。若給定γ、λ和σ2，則在已知觀測值y下的后驗概率分布可表示為p(θs，γ，λ，σ2|y)=p(θs|y，γ，λ，σ2)p(γ，λ，σ2|y)。其中p(θs|y，γ，λ，σ2)服從關于均值為μ、協方差為Σ的高斯分布，且有

μ=σ-2ΣΘTy，

(8)

Σ=(σ-2ΘTΘ+Λn-λ)-1，

(9)

式中：μ為θs的均值；Σ為θs的方差；Λn-λ=diag(1/γ1， 1/γ2，…，1/γj，…，1/γN)=diag{1/γj}，該矩陣的估計依賴于超參數λ.由于p(γ，λ，σ2|y)∝p(γ，λ，σ2，y)，求p(γ，λ，σ2，y)的解析式如下：

(10)

式中：p(y|θs,σ2)為y關于參數θs和σ2的先驗密度函數；p(θs|γ)為θs關于參數γ的先驗密度函數；p(γ|λ)為γ關于超參數λ的先驗密度函數；p(λ|v)為λ關于超參數v的先驗密度函數；p(σ-2)為σ2的先驗分布函數。

(10)式取對數并略去常數項后，分別對其求γj、λ、σ2偏導并令偏導等于0，根據超參數為正且v=0，解得γj的更新公式為

(11)

2.2 對多稀疏塊的Laplace BCS

(12)

(13)

λn=2(N-1)/(Σjrj)，

(14)

式中：μi，j表示第i塊數據中第j個參數γj的估計均值；Σi，jj表示第i塊數據的方差Σi的第j個對角線元素，其表達式見(9)式。利用(13)式、(14)式迭代計算復雜度較高，本文采用文獻[14]的快速RVM算法。

3 基于Lap-CBCS算法的機械振動信號壓縮重構

類似于機器學習中每次訓練樣本都會對模型參數進行優化，對于多數據塊的壓縮重構也是如此。Lap-CBCS算法參考多任務貝葉斯估計，將多塊數據看作多任務貝葉斯中的任務，構建稀疏塊Laplace先驗分布模型，使得每塊數據都能夠對先驗分布參數進行改進，就可以促進各塊的稀疏重構，大大提高貝葉斯參數估計的準確性。但想要將Lap-CBCS算法應用到實際機械振動信號的壓縮重構中，還需要解決信號分塊和稀疏域選取兩個問題。

3.1 機械振動信號的周期性分塊

由于齒輪箱振動信號具有非線性、非平穩性，其周期性表現得并不明顯，可稱其為類周期性。機械振動信號的類周期確定主要有如下3種方法：

1)根據轉速和采樣頻率計算類周期。假設轉速n=600 r/min，采樣頻率Fs=5 kHz，即采樣點數為5 000個/s，周期T=Fs/n=5 000/10=500個，由此可得每500個采樣點為一個類周期。

2)通過頻譜分析確定?？梢岳酶道锶~變換得到信號頻譜圖，分析信號頻譜圖，得出信號類周期。

3)通過包絡線確定。首先對信號作自相關分析，并求得自相關信號包絡線。然后取包絡線上相鄰兩峰值之間的距離作為類周期。

對機械振動信號進行傅里葉變換后會出現大量的頻譜峰值，難以識別類周期。而通過求包絡線計算誤差也較大，還需要對多個候選值求平均。因此，本文采用通過轉速和采樣頻率計算類周期的方法，轉速信息可通過轉速傳感器測得。

3.2 稀疏域選取

常見稀疏變換方法有離散余弦變換、小波變換等，稀疏域選取可根據信號稀疏度和相似度確定。

定義1稀疏度。假設原始信號x(x∈RN)在某變換域Ψ上是稀疏的，則稀疏信號θ與原始信號關系為x=Ψθ，變換后信號的0-范式‖Ψ-1x‖0，且Ψ-1x∈RM，則稀疏度可表示為K=‖Ψ-1x‖0/M.

定義2相似度。假設函數R(s，t)={δ|如果|si-ti|>ε則δi=1，否則δi=0，i=0，1，2，…，N}，其中ε為閾值，則兩組信號x1，x2∈RN的相似度可表示為D=1-‖R(Ψ-1x1，Ψ-1x2)‖0/M.

3.3 機械振動信號壓縮重構流程

本文提出了基于Laplace先驗和稀疏塊相關性的機械振動信號壓縮重構方法(見圖2)，具體包括以下步驟：

圖2 基于Laplace-CBCS算法的機械振動信號壓縮重構流程Fig.2 Compression and reconstruction processes of mechanical vibration signal based on Laplace-CBCS algorithm

步驟1周期性分塊。根據類周期對采集到的振動信號進行周期性分塊。利用轉速n和采樣頻率Fs計算類周期，計算公式為T=Fs/n.

步驟2信號壓縮。由于機械振動信號是一類非平穩復雜信號，為了提高其稀疏度，首先利用稀疏基Ψ將該信號投影到稀疏域，令原始信號x=Ψθ，得到稀疏域信號θ.利用傳感矩陣Θ=ΦΨ對稀疏信號進行狀態壓縮，得到壓縮后信號y=Θθ.

步驟3信號分塊傳輸。分別傳輸并接收分塊處理后的信號y1，…，yL.

步驟4x信號重構。利用本文Lap-CBCS重構算法以及從前端采集節點傳回的傳感矩陣Θ1，…，ΘL(如果觀測矩陣是確定性矩陣則無需傳回，如果是隨機性矩陣則需要傳回)對y1，…，yL聯合重構，得到稀疏信號θ1，…，θL.對稀疏信號進行逆稀疏變換，得到原始分塊信號x1，…，xL，將原始信號塊按順序拼接在一起，得到完整重構信號x′.重構信號與原始信號相比相位不發生變化。

4 齒輪箱數據實驗驗證

4.1 實驗準備

如圖3所示，實驗齒輪箱型號為JZQ175，三相電磁調速電動機功率4 kW，風冷磁粉制動器為齒輪箱提供載荷。數據采集由美國DYTRAN公司生產的4個3056B4型壓電傳感器組成。利用上述設備進行齒輪箱預制故障實驗，轉速800 r/min，輸入端負載10 N·m，采樣頻率20 kHz，采樣時間為6 s. 將上述實驗環境下采集到的數據進行信號壓縮、傳輸，重構仿真實驗。仿真環境為：MATLAB 2017b，Win10操作系統，i5-8250 CPU.

圖3 實驗臺示意圖Fig.3 Test-rig of gearbox

預置故障實驗主要包括兩種：齒根裂紋故障實驗和斷齒故障實驗。齒根裂紋故障加工在輸出軸大齒輪齒根上，斷齒故障加工在中間軸大齒輪上。本文進行4種預置故障實驗：正常狀態、1 mm齒根裂紋故障、2 mm齒輪斷齒故障、2 mm齒根裂紋故障(1 mm、2 mm表示裂紋深度)，采樣時間為6 s，采樣點數120 000，采集到的數據如圖4所示。

圖4 齒輪箱實測振動信號Fig.4 Measured vibration signals of gearbox

4.2 1 mm齒根裂紋故障仿真實驗

4.2.1 根據轉速和采樣頻率計算類周期

圖5 轉速信號Fig.5 Tachometer signal

4.2.2 稀疏域選擇

針對實測數據，分析離散余弦變換以及小波變換處理后稀疏信號的稀疏度與相似度。觀察處理數據發現，小波分解層數對于變換后信號的稀疏性、相似性影響都不大。而對于稀疏信號，稀疏度越小，相似度越高，越容易進行信號恢復。仿真實驗結果表明：離散余弦變換作為稀疏變換的處理效果是最好的。

4.2.3 利用Lap-CBCS算法對1 mm齒根裂紋故障數據進行仿真實驗

為對比分析，先利用BCS算法對80塊信號分別壓縮重構，再利用Lap-CBCS算法對80塊信號聯合壓縮重構，對比兩次重構。此外，本文擬采用標準均方根誤差(MSE)作為重構效果評價指標，其計算公式如(15)式所示：

MSE=‖u-u′‖2/‖u‖2，

(15)

式中：u表示原始信號；u′表示重構信號。

如圖6所示，隨機選取第1、15、22、36、61塊數據進行比較分析，發現BCS方法對每塊信號單獨壓縮重構MSE大于0.8，而利用Lap-CBCS算法重構可將MSE降低到0.3左右。

圖6 實測齒輪箱1 mm裂紋數據壓縮重構效果對比(80塊信號聯合重構)Fig.6 Comparison results of reconstructing signal blocks for measured 1 mm crack on root of tooth (simultaneous construction of 80 signal blocks)

將Lap-CBCS算法重構的80塊信號拼接起來，可以得到完整重構信號。由于實驗采樣點數較多，無法看出重構細節，選擇采樣點1×104～2×104的數據觀察細節。如圖7所示，Lap-CBCS重構數據點基本能夠覆蓋原數據點。在該范圍內，進一步選擇出1.1×104～1.2×104的數據觀察，發現Lap-CBCS算法重構只有在局部極值點處與原數據略有偏差，其余數據點基本重構成功。

圖7 實測齒輪箱1 mm裂紋局部數據重構細節Fig.7 Details of local data reconstruction for measured 1 mm crack

4.2.4 與其他算法對比

本文利用BP、OMP算法針對上述齒輪箱1 mm齒根裂紋故障數據進行重構實驗。實驗結果表明BP、OMP、BCS算法的標準MSE都在0.5以上，而Lap-CBCS算法的MSE能保持在0.3左右。此外，分析4種算法運行時間發現OMP算法所需時間最少，其次是BP算法，而Lap-CBCS算法運行時間最長。

4.3 齒輪不同故障狀態實驗

選取另外3種狀態：正常運行、2 mm齒輪斷齒故障、2 mm齒根裂紋故障下的實驗數據，由于轉速和采樣頻率不變，分塊周期仍為1 500數據點。對每種狀態數據分別利用BCS、Lap-CBCS、BP、OMP 4種算法進行壓縮重構分析。

選取觀測次數400、700、1 000(對應壓縮率3.75、2.14、1.50)，對每種不同狀態數據計算出MSE和運行時間平均值，得到表1所示結果。

由表1可見4種算法的標準MSE存在如下關系：MSELap-CBCS

表1 不同齒輪箱狀態下重構算法MSE與運行時間對比Tab.1 Calculated MSEs and running times of different algorithms under different gearbox conditions

1)隨著壓縮率降低，CS算法MSE降低，且算法運行時間增加。其中MSELap-CBCS最低，在壓縮率為2.14和1.5時可基本實現重構(MSE不大于0.3)，而BP、BCS算法只有在壓縮率1.5時才能實現重構(齒根裂紋故障狀態下仍不能實現重構)。另外，如果使用OMP算法則MSEOMP較大，表明該算法對于復雜機械振動信號這類非稀疏信號較為敏感，不適用于復雜齒輪箱數據的重構。

2)Lap-CBCS算法運行時間較長，對計算資源要求高。盡管如此，考慮到振動信號傳輸到上位機后后臺可利用的計算資源較多，且在某些特定場景下需要高精度的信號重構技術，因此Lap-CBCS算法具有較好的應用價值。

5 結論

本文將Laplace先驗模型和振動信號周期性稀疏塊相結合，提出一種改進的BCS方法。針對機械振動信號，提出利用轉速和采樣頻率計算類周期進行周期性分塊，并給出振動信號壓縮重構流程。所得主要結論如下：

1)齒輪箱實測信號驗證結果表明，盡管Lap-CBCS運行時間較長，但在標準MSE方面優于BCS、BP、OMP 3種典型算法。

2)Lap-CBCS方法較好地解決了復雜非平穩信號的壓縮重構問題，對于機械設備的無線監測具有一定的促進意義。

但是，Lap-CBCS算法雖然重構精度較高，卻需要較長的運行時間。如何改進算法、提高計算效率，將是下一步可能的研究方向。