余 路， 曲建嶺， 高 峰， 田沿平， 郭超然， 李 儼

(1.海軍航空大學(xué)青島校區(qū) 航空儀電控制工程與指揮系，青島 266041；2.西北工業(yè)大學(xué) 自動(dòng)化學(xué)院，西安 710072)

盲源分離(Blind Source Separation,BSS)是機(jī)械設(shè)備進(jìn)行復(fù)合故障診斷的一種有效方法。目前在機(jī)械設(shè)備的盲分離方法中，通常假設(shè)觀測(cè)信號(hào)的個(gè)數(shù)不小于源信號(hào)的個(gè)數(shù)，然而這在實(shí)際的故障診斷實(shí)踐中并不總是成立。受限于機(jī)械設(shè)備關(guān)鍵部件的安裝位置和傳感器的布置成本，觀測(cè)信號(hào)個(gè)數(shù)小于源信號(hào)個(gè)數(shù)的欠定盲源分離問題更加普遍，而單通道盲源分離(Single Channel Blind Source Separation，SCBSS)問題可以看作是欠定盲源分離問題中的特例。由于觀測(cè)信號(hào)個(gè)數(shù)越少意味著獲得的可用信息越少，對(duì)于源信號(hào)的恢復(fù)越困難，因此，單通道盲源分離問題一直以來是盲源分離問題研究的熱點(diǎn)和難點(diǎn)問題[1-2]。

多年來，已有眾多學(xué)者對(duì)此問題展開了研究。一些學(xué)者采用經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)及其變種理論的方法展開研究。例如，文獻(xiàn)[3]利用EMD算法增加觀測(cè)信號(hào)個(gè)數(shù)，將SCBSS問題轉(zhuǎn)化為正定或超定BSS問題，再使用成熟的FastICA算法求解。為了抑制EMD算法的端點(diǎn)效應(yīng)問題，文獻(xiàn)[4]采用EEMD算法(Ensemble EMD,EEMD)和基于四階累積量的特征矩陣聯(lián)合近似對(duì)角化算法(JADE，Joint Approximative Diagonalization of Eigenmatrix)并成功應(yīng)用于地聲信號(hào)的分離。這些算法雖然在一定條件下都可以較準(zhǔn)確的恢復(fù)源信號(hào)，但EMD算法本質(zhì)的端點(diǎn)效應(yīng)及模態(tài)混疊問題并沒有得到根本的解決，同時(shí)，這些方法改變了信號(hào)的結(jié)構(gòu)特征，受人工經(jīng)驗(yàn)的因素影響較大，因此給這類算法的廣泛應(yīng)用帶來了限制。還有部分學(xué)者考慮觀測(cè)信號(hào)的稀疏性問題，文獻(xiàn)[5]和文獻(xiàn)[6]采用時(shí)頻分析和聚類的方法解決欠定盲源分離問題，但對(duì)于SCBSS問題，線性稀疏性的條件常常不能滿足，而更復(fù)雜的非線性混合問題則需要更大的運(yùn)算量和新的求解算法，目前在這方面的研究還尚在起步階段。

基于以上分析，本文采用了另外一種思路求解SCBSS問題，即將移不變字典學(xué)習(xí)算法加入到SCBSS求解問題的框架中。

具有移不變特性的稀疏編碼算法最早由Smith在2006年的《Nature》中提出[7]，其本質(zhì)是將信號(hào)以多個(gè)基函數(shù)及其相應(yīng)系數(shù)卷積的形式表示，自適應(yīng)的學(xué)習(xí)蘊(yùn)含在信號(hào)中的移不變基函數(shù)，從而實(shí)現(xiàn)對(duì)原信號(hào)的稀疏表達(dá)。對(duì)于旋轉(zhuǎn)機(jī)械來說，反應(yīng)信號(hào)特征的分量往往都是重復(fù)出現(xiàn)的，比如故障軸承的沖擊，齒輪的磨損及旋翼轉(zhuǎn)速異常等[8]。在恒定工況下，不同激振源的信號(hào)特征相似，不同激振源的信號(hào)有所不同。利用移不變字典學(xué)習(xí)算法(Shift Invariant Dictionary Learning，SIDL)求解觀測(cè)信號(hào)中的移不變基函數(shù)，引入改進(jìn)的模糊C均值聚類算法自適應(yīng)地將估計(jì)源信號(hào)個(gè)數(shù)的同時(shí)，實(shí)現(xiàn)對(duì)各個(gè)移不變基函數(shù)時(shí)域重構(gòu)移不變分量進(jìn)行聚類，實(shí)現(xiàn)同類源信號(hào)的聚類和不同類源信號(hào)的區(qū)分，從而實(shí)現(xiàn)單通道振動(dòng)信號(hào)的盲源分離。

1 SIDL學(xué)習(xí)算法

1.1 SIDL學(xué)習(xí)算法簡(jiǎn)介

SIDL與標(biāo)準(zhǔn)稀疏編碼不同，其將輸入信號(hào)y表示成各個(gè)基函數(shù)與其系數(shù)的卷積

*xj+ε

(1)

式中，“*”表示卷積算子，dj∈Rm為基函數(shù)，M為基函數(shù)個(gè)數(shù)，xj∈Rl為稀疏系數(shù)，m為基函數(shù)長(zhǎng)度，l為信號(hào)的分段長(zhǎng)度，n為原始信號(hào)長(zhǎng)度，且三者滿足m

1.2 系數(shù)求解

求解稀疏系數(shù)的過程相當(dāng)于把所有基函數(shù)固定，求解每個(gè)基函數(shù)對(duì)應(yīng)的系數(shù)

(2)

1.3 字典學(xué)習(xí)

在該階段中，固定系數(shù)x求解字典矩陣D，則原優(yōu)化問題轉(zhuǎn)化為l2約束優(yōu)化問題

(3)

式中，N為輸入信號(hào)個(gè)數(shù)。直接求解上式計(jì)算量巨大，本文在字典學(xué)習(xí)階段采用時(shí)頻域轉(zhuǎn)化的思想，將時(shí)域計(jì)算量較大的卷積運(yùn)算轉(zhuǎn)換到頻域，由頻域的乘積運(yùn)算代替時(shí)域的卷積運(yùn)算，使運(yùn)算復(fù)雜度由O(N2)降低為O(NlgN)，提高了運(yùn)算效率。由Parseval定理可知，式(3)可轉(zhuǎn)化為如下優(yōu)化問題[13]

(4)

(5)

其中，λ∈RM為對(duì)偶變量，單位向量I∈RM，并且：

(6)

2 模糊C均值聚類算法

模糊C均值(Fuzzy C-Means，F(xiàn)CM)聚類算法已被廣泛應(yīng)用于數(shù)據(jù)挖掘和信號(hào)分析[14]，與傳統(tǒng)的基于硬閾值聚類的算法不同，F(xiàn)CM聚類算法利用模糊集理論實(shí)現(xiàn)對(duì)數(shù)據(jù)的軟劃分。文獻(xiàn)[15]指出，對(duì)于振動(dòng)信號(hào)這類非平穩(wěn)信號(hào)，采用軟閾值比硬閾值劃分往往能取得更佳的效果，因此本文采用模糊C均值聚類算法實(shí)現(xiàn)移不變分量的聚類。

給定數(shù)據(jù)集X=[x1,x2,…,xn]為s維空間Rs的有限觀測(cè)樣本，采用FCM算法針對(duì)每個(gè)數(shù)據(jù)點(diǎn)產(chǎn)生模糊隸屬關(guān)系uij=ui(xj)∈[0,1]。聚類中心由采樣點(diǎn)的數(shù)值及其隸屬化目標(biāo)函數(shù)的最小值確定

(7)

(8)

(9)

3 基于改進(jìn)字典學(xué)習(xí)的單通道振動(dòng)信號(hào)盲源分離算法

利用SIDL自學(xué)習(xí)觀測(cè)信號(hào)中的移不變基函數(shù)，將各個(gè)基函數(shù)進(jìn)行時(shí)域重構(gòu)得到移不變分量(Shift Invariant Component，SIC)，將各個(gè)SIC進(jìn)行分段并提取關(guān)鍵點(diǎn)，將FCM算法對(duì)各SIC進(jìn)行自適應(yīng)聚類，以聚類后源信號(hào)的相關(guān)度最小作為對(duì)源信號(hào)數(shù)目的最佳估計(jì)。輸出聚類后的SIC作為最終分離的源信號(hào)。

算法流程圖如圖1所示。本文算法的具體實(shí)現(xiàn)步驟如下。

圖1 本文算法簡(jiǎn)要流程圖

步驟1 將觀測(cè)信號(hào)進(jìn)行分段和平滑預(yù)處理。

步驟2 利用SIDL算法自學(xué)習(xí)信號(hào)中的移不變基函數(shù)。

步驟3 將各移不變基函數(shù)進(jìn)行時(shí)域重構(gòu)得到各個(gè)移不變分量(SIC)。

步驟4 求解各SIC的包絡(luò)譜。

步驟5 將各SIC的包絡(luò)譜進(jìn)行分段并執(zhí)行局部最大值檢測(cè)，得到反應(yīng)包絡(luò)譜概略波形特征的峰值點(diǎn)。

步驟6 將反映波形特征的峰值進(jìn)行模糊C均值聚類，聚類個(gè)數(shù)依次從2到SIC的總個(gè)數(shù)，依次計(jì)算不同聚類數(shù)下的恢復(fù)的源信號(hào)相關(guān)系數(shù)平均值。

步驟7 以相關(guān)系數(shù)平均值最低作為對(duì)源信號(hào)個(gè)數(shù)的估計(jì)，并對(duì)此時(shí)的SIC進(jìn)行合并，得到分離的源信號(hào)。

4 實(shí)驗(yàn)驗(yàn)證

4.1 實(shí)驗(yàn)設(shè)置

在大型設(shè)備中，齒輪和軸承是重要的機(jī)械部件。當(dāng)齒輪和軸承發(fā)生故障時(shí)，往往給機(jī)械設(shè)備的健康性能帶來?yè)p失。因此，本文根據(jù)齒輪與軸承的故障模型，分別構(gòu)造了故障齒輪信號(hào)和故障軸承信號(hào)的仿真模型[17-18]，同時(shí)加以白噪聲。

當(dāng)齒輪發(fā)生故障時(shí)，往往表現(xiàn)出幅值調(diào)制與相位調(diào)制的現(xiàn)象，其振動(dòng)信號(hào)模型為

(10)

式中，Am(t)和θm分別為第m次諧波的幅值和相位，f1為齒輪嚙合頻率。在本例中取M=5，Am(t)=2.5[sin(60πt)+2]/m，θm=mπ/4，f1=340 Hz。

軸承發(fā)生故障時(shí)往往伴隨著振動(dòng)信號(hào)的沖擊響應(yīng)，因此故障軸承的模型可以簡(jiǎn)化為

y2=a(t)e-ptsin(2πf2t+φ)

(11)

式中，a(t)為調(diào)幅幅值，p為衰減系數(shù)，f2為軸承固有頻率，φ為相位。在本例中取a(t)=20[sin(60πt)+2]，φ=0，f2=3 400 Hz，T=0.01 s，采樣頻率fs=12 kHz。觀測(cè)信號(hào)可以表示為y=y1+y2+yn，其中yn為高斯白噪聲，使觀測(cè)信號(hào)與原始信號(hào)的信噪比為10 dB，源信號(hào)與觀測(cè)信號(hào)如圖2所示(為了更好的展示細(xì)節(jié)，各個(gè)信號(hào)取1 000個(gè)點(diǎn)，下同)。應(yīng)用本文提出的字典學(xué)習(xí)算法，設(shè)置參數(shù)如表1所示。

圖2 源信號(hào)及加噪混合信號(hào)

參數(shù)數(shù)值基函數(shù)長(zhǎng)度50基函數(shù)個(gè)數(shù)6子信號(hào)長(zhǎng)度256子信號(hào)重疊率50%稀疏度160迭代次數(shù)60

4.2 仿真實(shí)驗(yàn)

經(jīng)過60次迭代后算法收斂，學(xué)習(xí)到的基函數(shù)及其對(duì)應(yīng)的移不變分量(SIC)如圖3所示，從基函數(shù)中可以隱約看出，故障軸承的沖擊響應(yīng)與故障齒輪的調(diào)制特性，但從對(duì)應(yīng)的SIC中可以比較明顯的看出兩種故障已經(jīng)實(shí)現(xiàn)了區(qū)分。為了以數(shù)值計(jì)算的方式而不是主觀臆斷作為兩種故障區(qū)分的依據(jù)，考慮將各個(gè)SIC轉(zhuǎn)化至頻域，采用包絡(luò)譜對(duì)其進(jìn)行進(jìn)一步分析。

由于采樣點(diǎn)數(shù)以及采樣頻率的原因，原始包絡(luò)譜呈現(xiàn)出凌亂復(fù)雜的波形，不利于后期的處理分析。因此本文提出了局部最大值檢測(cè)的處理策略，即將包絡(luò)譜進(jìn)行不重疊分段處理，以每一區(qū)間段內(nèi)的最大值點(diǎn)作為關(guān)鍵點(diǎn)，將所有關(guān)鍵點(diǎn)平滑連接形成了易于理解與處理的包絡(luò)譜波形特征，同時(shí)波形的峰值關(guān)鍵點(diǎn)可以作為接下來模糊C均值聚類的特征向量加以分析處理，如圖4所示?？紤]到信號(hào)源個(gè)數(shù)的未知性，本文采取自適應(yīng)的源信號(hào)數(shù)目估計(jì)，即將聚類個(gè)數(shù)從2開始逐一增加到預(yù)先設(shè)定的基函數(shù)個(gè)數(shù)，在每個(gè)聚類數(shù)目下執(zhí)行模糊C均值聚類得到源信號(hào)的估計(jì)，若達(dá)到理想的分離效果，則分離的源信號(hào)之間的皮爾遜相關(guān)系數(shù)的均值達(dá)到最小。在本例中，聚類數(shù)目與分離的源信號(hào)相關(guān)系數(shù)均值(corrmean)的關(guān)系如圖5所示，可以看出，當(dāng)取值為2時(shí)達(dá)到最小值，這也符合本例的實(shí)際情況。

圖3 學(xué)習(xí)到的基函數(shù)及其對(duì)應(yīng)的時(shí)域移不變分量

圖4 各個(gè)移不變分量的包絡(luò)譜(實(shí)線)及局部最大值檢測(cè)(虛線)的結(jié)果

Fig.4 Envelope spectrums of SICs(line) and results of local maximum detection(dotted line)

圖5 聚類數(shù)與corrmean關(guān)系

聚類數(shù)目為2時(shí)分離的源信號(hào)的時(shí)域波形如圖6所示。從圖中可以看出，分離效果良好，保留了故障軸承和故障齒輪的主要特征信息，只是在波形的幅值上與原始信號(hào)存在一定的出入。

圖6 本文算法分離結(jié)果

作為對(duì)比，分別以文獻(xiàn)[3]和文獻(xiàn)[4]的方法進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果分別如圖7和圖8所示，可以看出并沒有達(dá)到良好的分離效果。為了定量衡量本文算法與其他算法的分離效果，采用峰值信噪比(SNR)和均方差(MSE)作為指標(biāo)參數(shù)，不同算法的實(shí)驗(yàn)結(jié)果如表2所示?？梢钥闯?，本文算法在三者中保持最高信噪比的同時(shí)具有最低的均方差，表明本文算法較好的保持了源信號(hào)的特征。

圖7 EMD+ICA法分離結(jié)果

圖8 EEMD+JADE法分離結(jié)果

4.3 實(shí)測(cè)航空發(fā)動(dòng)機(jī)振動(dòng)信號(hào)實(shí)驗(yàn)

為驗(yàn)證本文算法在實(shí)際裝備中的可行性，將其應(yīng)用于某型直升機(jī)振動(dòng)數(shù)據(jù)采集與處理。眾所周知，直升機(jī)的飛行安全至關(guān)重要，一旦機(jī)體發(fā)生故障將會(huì)給機(jī)組人員的生命和裝備造成損失，因此常利用振動(dòng)信號(hào)傳感器采集機(jī)體各位置處的振動(dòng)信號(hào)作為故障診斷的依據(jù)，但受限于直升機(jī)各部件的緊湊性，單個(gè)傳感器采集到的信號(hào)可能是多個(gè)機(jī)體振源的疊加，因此有必要對(duì)其進(jìn)行分離作進(jìn)一步的判斷。直升機(jī)部分振動(dòng)傳感器的安裝如圖9所示。

表2不同算法分離效果對(duì)比

Tab.2Comparisonofseparatingresultsofdifferentmethods

方法y1y2SNRMSESNRMSEEMD+ICA(文獻(xiàn)[3])2.88.353.55.38EEMD+JADE(文獻(xiàn)[4])5.23.954.34.46本文算法13.51.0514.21.21

將位于右發(fā)動(dòng)機(jī)前部的振動(dòng)傳感器采集到的信號(hào)作為研究對(duì)象，其時(shí)域波形和頻域包絡(luò)譜如圖10所示，從圖中無法看出明顯的特征。設(shè)置基函數(shù)個(gè)數(shù)為4，其余參數(shù)與表1相同，應(yīng)用本文的字典學(xué)習(xí)算法學(xué)習(xí)到的基函數(shù)如圖11所示，對(duì)其時(shí)域重構(gòu)的移不變分量(SIC)進(jìn)行包絡(luò)譜分析及局域最大值檢測(cè)，實(shí)驗(yàn)結(jié)果如圖12所示。從圖中可以看出SIC1與SIC3的包絡(luò)譜略有相似，SIC2與SIC4的包絡(luò)譜也有相似之處，為了驗(yàn)證猜想，應(yīng)用模糊C均值聚類算法，聚類個(gè)數(shù)從2到4，發(fā)現(xiàn)在聚類個(gè)數(shù)為2時(shí)達(dá)到corrmean最小值，因此判斷采集到的振動(dòng)信號(hào)包含2個(gè)信號(hào)源。

(a) 觀測(cè)信號(hào)時(shí)域波形

(b) 觀測(cè)信號(hào)包絡(luò)圖

圖11 學(xué)習(xí)到的基函數(shù)

圖12 各移不變分量的包絡(luò)譜及局部最大值檢測(cè)

Fig.12 Envelope spectrums of SICs and local maximum detection

分離的時(shí)域信號(hào)如圖13所示，對(duì)各個(gè)分離信號(hào)求取包絡(luò)譜如圖14所示，可以看出分離信號(hào)1為慢車時(shí)右發(fā)動(dòng)機(jī)轉(zhuǎn)子的轉(zhuǎn)速fm(293 Hz)，分離信號(hào)2為減速器齒輪振動(dòng)信號(hào)(fc，2fc，3fc，4fc)。因此該振動(dòng)傳感器采集到的信號(hào)已成功得到分離，一旦某部位發(fā)生故障異常振動(dòng)時(shí)可以準(zhǔn)確定位到具體的位置。

(a) 分離信號(hào)1時(shí)域波形

(b) 分離信號(hào)2時(shí)域波形

(a) 分離信號(hào)1包絡(luò)圖

(b) 分離信號(hào)2包絡(luò)圖

作為對(duì)比，分別將文獻(xiàn)[3]和文獻(xiàn)[4]的方法應(yīng)用于上述實(shí)驗(yàn)數(shù)據(jù)，分離后信號(hào)頻譜分別如圖15和16所示。

(a) 分離信號(hào)1包絡(luò)圖

(b) 分離信號(hào)2包絡(luò)圖

Fig.15 Envelope spectrum of separated signals by paper[3]

(a) 分離信號(hào)1包絡(luò)圖

(b) 分離信號(hào)2包絡(luò)圖

Fig.16 Envelope spectrum of separated signals by paper[4]

通過分離信號(hào)的包絡(luò)譜對(duì)比可以看出，文獻(xiàn)[3]和文獻(xiàn)[4]方法均一定程度上存在兩種信號(hào)源混疊的現(xiàn)象。例如，文獻(xiàn)[3]中的分離信號(hào)1中含有減速器齒輪振動(dòng)信號(hào)的各階譜，但也部分含有右發(fā)動(dòng)機(jī)轉(zhuǎn)子轉(zhuǎn)速的頻率，給分離信號(hào)的識(shí)別增加了困難。同理，文獻(xiàn)[4]方法的分離結(jié)果也存在類似的問題，表明分離效果較本文算法存在差距。

取三種算法各運(yùn)行20次的平均運(yùn)算時(shí)間結(jié)果進(jìn)行對(duì)比，實(shí)驗(yàn)軟件環(huán)境配置采用Linux系統(tǒng)下Matlab2014a，硬件環(huán)境采用i5HQ6300處理器，8 GB內(nèi)存。實(shí)驗(yàn)結(jié)果如表3所示。

表3 不同算法運(yùn)算時(shí)間對(duì)比

不可否認(rèn)，文獻(xiàn)[3]和文獻(xiàn)[4]算法相比于本文算法具有更快的運(yùn)算速度，本文算法的時(shí)間花費(fèi)主要在于移不變基函數(shù)的迭代求解，但通過迭代求解的移不變基函數(shù)可以更好地學(xué)習(xí)到觀測(cè)信號(hào)中存在的本質(zhì)特征，因此可以獲得更佳的信號(hào)分離效果。

5 結(jié) 論

本文提出了一種采用改進(jìn)字典學(xué)習(xí)對(duì)單通道振動(dòng)信號(hào)進(jìn)行特征自學(xué)習(xí)和盲源分離的解決方案。移不變字典學(xué)習(xí)算法能夠自學(xué)習(xí)反映信號(hào)本質(zhì)特征的移不變基函數(shù)，解決了特征提取依賴專家知識(shí)與人工提取的困難；以局部最大值檢測(cè)算法提取各個(gè)移不變分量的包絡(luò)譜關(guān)鍵點(diǎn)，降低了后期聚類運(yùn)算的效率和準(zhǔn)確性，而且更有利于反映頻域信號(hào)的波形特征；利用自適應(yīng)模糊C均值聚類算法實(shí)現(xiàn)對(duì)源信號(hào)的恢復(fù)，解決了源數(shù)估計(jì)和信號(hào)恢復(fù)的問題。由仿真和實(shí)測(cè)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果可知，本文提出的算法具有較高的源信號(hào)恢復(fù)精度，有利于對(duì)單通道振動(dòng)信號(hào)的進(jìn)一步分析，具有較強(qiáng)的工程應(yīng)用推廣價(jià)值。