基于K-奇異值分解字典學(xué)習(xí)的振動信號壓縮感知方法

何天遠(yuǎn)， 王萬仁2， 吳魯明， 邢亞航， 郝如江

(1. 石家莊鐵道大學(xué)機(jī)械工程學(xué)院， 河北石家莊050043； 2. 中鐵二十二局集團(tuán)軌道工程有限公司， 北京100040)

齒輪箱的狀態(tài)監(jiān)測和故障診斷對預(yù)防突發(fā)事故和維持機(jī)械設(shè)備正常運(yùn)行具有重大的意義。在齒輪箱故障診斷研究中，很少有對信號采集方面的研究。目前被廣泛采用的振動信號采樣方法是以采樣頻率必須大于信號中最高頻率2倍的奈奎斯特采樣定理為基礎(chǔ)的。近年來，隨著大型機(jī)械設(shè)備向綜合、高速、連續(xù)和網(wǎng)絡(luò)化趨勢發(fā)展，如果依然用基于奈奎斯特采樣定理的采集系統(tǒng)采集信號，就會產(chǎn)生較高的采樣頻率，同時會產(chǎn)生海量的監(jiān)測數(shù)據(jù)[1]，因此，齒輪箱狀態(tài)監(jiān)測得到的海量數(shù)據(jù)面臨著難以實(shí)時傳輸和儲存過慢的問題。近幾年出現(xiàn)的壓縮感知(compressed sensing, CS)采樣突破了傳統(tǒng)的奈奎斯特采樣定理[2]的限制。CS理論通過以遠(yuǎn)低于奈奎斯特采樣速率的頻率對信號進(jìn)行壓縮采樣并且直接得到信號的壓縮表示，不需要經(jīng)過信號采樣和壓縮2個單獨(dú)的過程，就達(dá)到采樣和壓縮同時進(jìn)行的目的。最后，利用適當(dāng)?shù)腃S重構(gòu)算法就可以高概率地重構(gòu)出原始信號[3]。

針對上述問題，本文中提出基于K-奇異值分解(K-singular value decomposition, K-SVD)字典學(xué)習(xí)的振動信號CS方法，用以提高重構(gòu)信號與原始信號的相似度，并且最大程度地保留原始信號的頻率特征。實(shí)驗(yàn)分別以基于K-SVD訓(xùn)練生成的2種字典和離散余弦變換(discrete cosine transform, DCT)固定字典為稀疏矩陣，用L1范數(shù)重構(gòu)算法對原始信號進(jìn)行重構(gòu)，并對實(shí)驗(yàn)結(jié)果進(jìn)行對比分析。

1 CS理論

CS利用數(shù)據(jù)的冗余特性，只采集少量有用的樣本來還原真實(shí)的采集數(shù)據(jù)[4]。在采集數(shù)據(jù)時直接采集q個有效的測量值，而不是滿足傳統(tǒng)的奈奎斯特采樣定律的要求，即采樣頻率必須是信號最高頻率的2倍或以上。CS以采集的原始信號在變換域上稀疏作為前提條件，用測量矩陣觀測被測信號，再由觀測值結(jié)合重構(gòu)算法，如L1范數(shù)算法、正交匹配追蹤(orthogonal matching pursuit, OMP)算法等來重構(gòu)原始信號[5]。

CS理論的公式為

y=Φf

，

(1)

式中：y為觀測矢量且y∈m×1，m為觀測矢量的維度；Φ為測量矩陣；f為原始信號且f∈n×1，n為原始信號的維度。f可表示為

f=φα

，

(2)

式中：φ為正交基矩陣；α為稀疏矢量。

CS理論公式中的測量矩陣Φ需要滿足一定的不相關(guān)性，即與φ不相關(guān)。由式(1)、(2)聯(lián)立可得表達(dá)式

y=Φφα

。

(3)

稀疏矢量是r稀疏的，且r

(4)

求解式(4)常用的有L1范數(shù)算法、基追蹤算法和OMP算法等。CS有2個重要的條件，即稀疏性和約束等距性(restricted isometry property, RIP)[6]。原始信號f是稀疏的，或者在某個變換域內(nèi)是稀疏的。測量矩陣Φ則需要滿足RIP條件，

(5)

式中A為感知矩陣。當(dāng)Φ為高斯隨機(jī)矩陣時，雖然完備字典為非正交基，但感知矩陣也具有隨機(jī)性，仍可滿足CS的RIP性質(zhì)。

2 K-SVD方法

Aharon等[7]提出了一種K-SVD的字典訓(xùn)練方法。K-SVD類似于K均值算法的一種泛化形式，其核心思想是通過稀疏約束問題的追蹤和奇異值分解算法的交替進(jìn)行，實(shí)現(xiàn)字典與稀疏系數(shù)的更新。K-SVD算法需要構(gòu)造稀疏字典對數(shù)據(jù)進(jìn)行稀疏表示，在語音信號壓縮、圖像壓縮、體征信號處理等方面已被廣泛應(yīng)用。

2.1 稀疏表示

稀疏表示的理論公式為

Y=dX

，

(6)

式中：Y為樣本矩陣且Y∈w×N，w為字典中每個樣本的維度，N為樣本數(shù)；d為過完備字典且d∈w×k(k>w)，k為一個系數(shù)的長度或訓(xùn)練字典中的原子數(shù)；X為系數(shù)矩陣且X∈k×N。

式(6)的目的是找到讓系數(shù)矩陣X盡量稀疏的字典d，尋找最優(yōu)解(X最稀疏)是非確定性多項(xiàng)式難題(non-deterministic polynomial hard, NP-Hard)，只能逐步逼近最優(yōu)解。

創(chuàng)建一個初始字典，可以是DCT固定字典、小波字典，也可以由采集到的振動信號生成。用該字典對采集的振動信號稀疏表示，即讓數(shù)據(jù)表達(dá)盡可能全面，系數(shù)盡可能小，并得到稀疏表示的系數(shù)矩陣X。將dX分解，也就是用d中的行乘與之對應(yīng)X中的列，即

(7)

式中：di為過完備字典d的行；xi為系數(shù)矩陣X的列。最后將dixi逐個優(yōu)化。

2.2 字典更新

直接把初始字典作為過完備字典會使具有稀疏性的系數(shù)矩陣所表示的數(shù)據(jù)與原始數(shù)據(jù)相比誤差較大。減小誤差的方法是在滿足稀疏度的條件下逐行逐列更新優(yōu)化，逼近可用字典。計(jì)算當(dāng)前誤差矩陣的表達(dá)式為

(8)

產(chǎn)生的誤差為

,

(9)

在分離第k個基分量之后會產(chǎn)生一個空位，而SVD的作用就是找到更適合的基分量來填充這個空位。當(dāng)誤差穩(wěn)定時，更新的字典基本收斂。

2.3 求解流程

K-SVD的求解是一個迭代過程。首先，假設(shè)字典d是固定的，用OMP、基追蹤、L1范數(shù)等算法，都可以得到使Y稀疏表示的系數(shù)矩陣X，然后固定X，根據(jù)X來更新字典d，如此循環(huán)直到收斂為止[8]。

在求解過程中過完備字典d需要逐列更新。首先，假設(shè)系數(shù)矩陣X和過完備字典d是固定的，當(dāng)過完備字典中的第k列dk更新時，dk對應(yīng)系數(shù)矩陣X中第k行為xk，則

(10)

3 仿真分析

3.1 仿真信號

為了驗(yàn)證基于K-SVD字典學(xué)習(xí)的振動信號CS方法是可以提高信號重構(gòu)性能的，設(shè)置如下仿真信號進(jìn)行分析。此仿真信號模擬的是軸承內(nèi)圈局部故障時的振動信號。信號的采樣頻率為16 000 Hz，采樣點(diǎn)數(shù)為4 000。

U(t-hT-τh)2[1-0.5cos(2πfrt)] +n(t)，

(11)

式中：t為時間節(jié)點(diǎn)；沖擊次數(shù)h=1, 2, … ；fr為轉(zhuǎn)頻，值為20 Hz；fn為共振頻率，值為5 000 Hz；C為信號衰減系數(shù)， 值為1 000；τh是第h次沖擊相對于特征周期T的微小波動， 隨機(jī)波動服從均值為0的正態(tài)分布； 標(biāo)準(zhǔn)差為轉(zhuǎn)頻的0.5%；n(t)為高斯噪聲， 標(biāo)準(zhǔn)差為2， 加入噪聲后的信號信噪比為-11.27 dB；U(t)為單位階躍函數(shù)。 計(jì)算所得故障特征頻率fu=3.906 Hz[9]。

沖擊信號波形圖和仿真信號波形圖如圖1所示。從圖中可以明顯看出，噪聲完全掩蓋了周期沖擊成分，沒有規(guī)律可循。

(a)沖擊信號

(b)仿真信號圖1 沖擊信號與仿真信號波形圖

實(shí)驗(yàn)分為3組，在采集的振動信號中隨機(jī)選取500個點(diǎn)作為CS測量數(shù)據(jù)x1，以3組不同的特定字典為基字典，最后用L1范數(shù)算法重構(gòu)信號。

3組特定的字典分別如下：

1)以DCT固定字典為基字典d1；

2)以一部分采集的振動信號為訓(xùn)練樣本，DCT固定字典為初始字典，利用K-SVD模型生成特定基字典d2；

3)以一部分采集的振動信號為訓(xùn)練樣本，用訓(xùn)練樣本生成矩陣作為初始字典，利用K-SVD模型生成特定基字典d3。

3.2 仿真結(jié)果對比分析

當(dāng)壓縮率為60%時， 仿真信號頻域圖如圖2(a)所示； 以d1、d2、d3為基字典的重構(gòu)信號頻域圖如圖2(b)、 (c)、 (d)所示。 對比圖可看出， 經(jīng)過K-SVD訓(xùn)練過的字典在振動信號壓縮采樣中的效果更好，而且圖2(b)中突出頻率被淹沒在毛刺中，突出頻率十分不明顯。

(a)仿真信號頻域圖

(b)基字典為d1，重構(gòu)頻域圖

(c)基字典為d2，重構(gòu)頻域圖

(d)基字典為d3，重構(gòu)頻域圖圖2 仿真信號頻域圖與基字典為d1、 d2、 d3時的重構(gòu)頻域圖

本文中采用弗雷歇距離δ來評價(jià)重構(gòu)信號與原信號的相似度，采用壓縮率來衡量振動信號的壓縮程度。

弗雷歇距離是1906年法國數(shù)學(xué)家Maurice René Fréchet提出的一種基于空間路徑相似度的描述方式，其著重將路徑空間距離考慮進(jìn)去，對于有一定空間時序的曲線相似度的評價(jià)效率更高[10]。

弗雷歇距離就是使最大距離在最小化采樣方式下的比值。弗雷歇距離越小，說明重構(gòu)信號與原始信號的差異越小，2條曲線相似度越高。其數(shù)學(xué)表達(dá)式為

(12)

式中：j為原信號軌跡PF的長度；u為重構(gòu)信號軌跡QF的長度；α(b)為原信號運(yùn)動位置描述函數(shù)；β(b)為重構(gòu)信號運(yùn)動位置描述函數(shù)； 若將變量b約束到區(qū)間[0, 1]內(nèi)， 則有α(0)=0，α(1)=j，β(0)=0，β(1)=u。

在壓縮率分別為30%、40%、50%、60%、70%、80%、90%時，計(jì)算原始信號與字典d1產(chǎn)生的重構(gòu)信號的弗雷歇距離δ1、 原始信號與字典d2產(chǎn)生的重構(gòu)信號的弗雷歇距離δ2、 原始信號與字典d3產(chǎn)生的重構(gòu)信號的弗雷歇距離δ3。

仿真實(shí)驗(yàn)結(jié)果如圖3所示。從圖中可以很明顯地觀察出，在壓縮率相同時,經(jīng)過K-SVD訓(xùn)練過的字典能更好地還原信號，而且以訓(xùn)練樣本為初始字典，經(jīng)過K-SVD模型訓(xùn)練生成的特定字典d3的效果最好。

圖3 仿真信號在不同壓縮率時的弗雷歇距離

4 實(shí)驗(yàn)分析

4.1 實(shí)驗(yàn)介紹

采用圖4所示的動力傳動故障診斷綜合實(shí)驗(yàn)臺，采集振動信號并對其進(jìn)行實(shí)驗(yàn)驗(yàn)證。實(shí)驗(yàn)中電機(jī)的轉(zhuǎn)頻為40 Hz，采樣頻率為12.5 kHz。

圖4 動力傳動故障診斷綜合實(shí)驗(yàn)臺

實(shí)驗(yàn)同樣分為3組，在采集的振動信號中隨機(jī)選取500個點(diǎn)作為CS測量數(shù)據(jù)x2，以3組不同的特定字典為基字典，最后用L1范數(shù)算法重構(gòu)信號。

3組特定的字典分別如下：

1)以DCT固定字典為基字典d4；

2)以一部分采集的振動信號為訓(xùn)練樣本，DCT固定字典為初始字典，利用K-SVD模型生成特定字典d5；

3)以一部分采集的振動信號為訓(xùn)練樣本，用訓(xùn)練樣本生成的矩陣作為初始字典，利用K-SVD模型生成特定字典d6。

4.2 實(shí)驗(yàn)結(jié)果對比分析

壓縮率為60%時，實(shí)測信號頻域圖如圖5(a)所示；以d4、d5、d6為基字典的重構(gòu)信號頻域圖如圖5(b)、 (c)、 (d)所示。對比圖5可知特定字典d6在振動信號壓縮采樣中的效果更好，而且圖5(b)、 (c)中的部分突出頻率被淹沒在毛刺中，突出頻率極不明顯。

實(shí)驗(yàn)結(jié)果如圖6所示。 從圖中可以明顯地觀察出，在壓縮率相同時經(jīng)過K-SVD訓(xùn)練過的字典能更好地還原信號，且此方法對于實(shí)測信號適用。

5 結(jié)語

為了對振動信號更好地進(jìn)行壓縮傳輸，本文中提出了一種基于K-SVD字典學(xué)習(xí)的振動信號CS方法。實(shí)驗(yàn)中分別利用3 種不同的字典作為基字典對信號進(jìn)行重構(gòu)，實(shí)驗(yàn)結(jié)果表明，經(jīng)過K-SVD訓(xùn)練過的字典能有效地提高重構(gòu)精度。

(a)實(shí)測信號頻域圖

(b)基字典為d4，重構(gòu)頻域圖

(c)基字典為d5，重構(gòu)頻域圖

(d)基字典為d6，重構(gòu)頻域圖圖5 實(shí)測信號頻域圖與基字典為d4、 d5、 d6時的重構(gòu)頻域圖

圖6 實(shí)測信號在不同壓縮率時的弗雷歇距離