朱亞靜，李志農(nóng)，，楊曉飛，李澤東，陶俊勇

(1.南昌航空大學(xué) 無損檢測教育部重點實驗室，南昌 330063;2.國防科學(xué)技術(shù)大學(xué) 裝備綜合保障技術(shù)重點實驗室，長沙 410073)

1 引言

盲分離(blind source separation，BSS)可以在源信號和信號傳輸?shù)幕旌贤ǖ赖认闰炛R未知的情況下分離觀測信號獲取有用源信號。傳統(tǒng)機械故障盲分離方法中常采用二階統(tǒng)計量的方法進行混合矩陣和源信號的估計，例如，劉歡[1]在風(fēng)力機主軸承的故障診斷中引入二階統(tǒng)計量BSS算法該應(yīng)用取得了良好的效果。文獻[2]中將EMD和ICA相結(jié)合用于去除風(fēng)電機機組的強干擾信息，分離混合信號得到源信號估計。但是這些傳統(tǒng)的源數(shù)估計方法和故障源信號識別方法均屬于二維數(shù)據(jù)處理方法，采用矩陣分解，需要添加一些約束或者限制條件(比如恒模約束以及正交約束)才能使其分解具有唯一性。實際應(yīng)用中并不能滿足這些苛刻的約束條件，因此需要研究新的方法來解決這些不足。

PARAFA因分解具有唯一性而備受關(guān)注。PARAFAC概念開始于心理學(xué)領(lǐng)域，后來快速的在許多領(lǐng)域得到發(fā)展和應(yīng)用[3]。尤其，PARAFAC能夠充分利用信號的代數(shù)性質(zhì)使其在通訊和機械故障信號處理等領(lǐng)域發(fā)展迅速，其主要通過多維數(shù)據(jù)的擬合迭代求得各種參數(shù)估計[4-5]。楊若男等[6]為了解決SIMO-OFDM系統(tǒng)下的信道估計和符號檢測問題，建立了PARAFAC模型，利用此模型中的行滿秩特性，并且結(jié)合奇異值分解，提出了一種聯(lián)合盲估計的求解方法，成功減少了計算的復(fù)雜度，同時也避免了性能下降的問題。但是實際的應(yīng)用中，使用現(xiàn)有的三線性平行因子方法實現(xiàn)機械故障的診斷時，只能對采集信號的傳感器數(shù)目，振動信號的分段數(shù)以及每個數(shù)據(jù)段包含的數(shù)據(jù)點數(shù)信息進行建模，振動信號中的時間序列信息被完全忽略，振動信號包含的信息不完整，所以有必要將三維平行因子方法擴展至四維。目前，三線性平行因子向四線性平行因子的擴展有許多的進展。Shang F[7]在三線性平行因子模型的基礎(chǔ)上，加入一個維度使其成為四線性平行因子模型，進行四維數(shù)據(jù)的獲取與分析，結(jié)果表明四線性平行因子比三線性平行因子效果更好。文獻[8]中通過使用模擬和實驗數(shù)據(jù)集比較PARAFAC和四線性PARAFAC，結(jié)果表明四線性PARAFAC可以從數(shù)據(jù)中提取更多的信息，而且證實了與三線性平行因子分解算法相比四線性平行因子具有更高的參數(shù)估計精度。隨著研究的深入諸多文獻均表明四線性PARAFAC比三線性PARAFAC有更好的應(yīng)用價值。

基于以上理念，本文中提出了一種基于四線性平行因子的機械故障源盲分離方法，四線性盲分離模型相比于傳統(tǒng)的盲分離優(yōu)勢明顯，主要表現(xiàn)在：分解在寬松條件下具有唯一性，可擴展范圍廣，利用數(shù)學(xué)分離的思想代替?zhèn)鹘y(tǒng)的信號分離的思想。該方法相比于在傳統(tǒng)平行因子分解算法主要優(yōu)勢主要表現(xiàn)在：① 在三維的基礎(chǔ)引入第四維解決振動信號時序信息被忽略的問題，數(shù)據(jù)包容量大，建模時包含的信息更加完整；② 在寬松條件下分解具有唯一性，可擴展范圍廣。利用四線性交替最小二乘進行迭代擬合，收斂更加平穩(wěn)，預(yù)測能力更強。

2 四線性平行因子理論基礎(chǔ)

四維矩陣X∈CI×J×K×L的四線性成分模型分解示意圖如圖1。

圖1 四線性成分模型分解直觀分解示意圖Fig.1 Intuitive decomposition diagram of quadrilinear component model decomposition

其中，X是一個四維矩陣，將X分解為大小分別為I×R，J×R，K×R，N×R的相關(guān)矩陣A，B，C，D。G是一個四維核對角矩陣，維數(shù)為I×J×K×N，并且超對角元素為1，其他元素為0。E為一個I×J×K×N四維殘差矩陣。

矩陣X∈CI×J×K×N的四線性成分模型標(biāo)量形式通常可以表示為：

(1)

其中，A∈CI×R，B∈CJ×R，C∈CK×R，D∈CN×R分別為四線性成分模型的4個承載矩陣。k=1，2,…,K，n=1，2,…,N，i=1，2,…,I，j=1，2,…,J，E∈CI×J×K×N為四維殘差矩陣。

四線性平行因子模型是將四線性成分模型中的四維矩陣X鋪展成二維矩陣形式，沿4個方向進行交替分解得到，有4種表述形式：

(2)

其中，Di(A)，Dj(B)，Dk(C)，Dn(D)為對角化算子，Di(A)表示矩陣其余的位置均為零，主對角線元素以矩陣A的第i行作為元素的新生成矩陣。其他對角化算子表示的含義與Di(A)表示的含義相似。

3 四線性平行因子盲分離模型

線性瞬時混疊模型可以表示為：

X(t)=Α(t)·S(t)+E(t)

(3)

式中：X(t)=[x1(t)，x2(t)，…,xN(t)]T為N維觀測信號；S(t)=[s1(t)，s2(t)，…，sM(t)]T為M個源信號；Α(t)是維度為N×M的混合矩陣；E(t)是噪聲信號。假定一個時間段內(nèi)每個觀測信號的數(shù)據(jù)點為L個(共N個傳感器)，將觀測信號中心化處理，接著將處理后的觀測信號平均截成不會重復(fù)的數(shù)據(jù)段J個，每一個數(shù)據(jù)段里面含有的數(shù)據(jù)點數(shù)為L/J個，令L/J等于K。采樣時間為t，每一段數(shù)據(jù)對應(yīng)的采集時間為h，則共分為t/h個時間段記其為I。則式(3)可以表示為：

X(i，h)=A·S(i，h)+E(i，h)

(4)

計算各段時滯協(xié)方差矩陣Rk：

Rk(tk，h)=ARx(tk，h)AT

(5)

其中，k=1，2,…,K，序號tk表示t內(nèi)第k個數(shù)據(jù)點，式中：

(6)

假設(shè)每個源信號均不相關(guān)，則塊時滯協(xié)方差矩陣化可以化為對角矩陣，令這個對角矩陣為B，則式(5)化為

Rk(tk，h)=Adiag(B)AT

(7)

接著計算一個時間段內(nèi)N個觀測信號的源時滯協(xié)方差Rn：

Rn(tn，h)=ARS(tn，h)AT

(8)

式中tn表示t內(nèi)第n個觀測信號，式(8)中RS(tn，h)可以表示為：

(9)

將源時滯協(xié)方差矩陣化為對角矩陣，并令這個對角矩陣為C，則式(7)化為：

Rn(tn，h)=Adiag(C)AT

(10)

所有的源時滯協(xié)方差矩陣和塊時滯協(xié)方差矩陣疊加為三階張量，并且將這個三階張量記為Rk，n：

(11)

其中，diag(A)，diag(B)，diag(C)，diag(D)是載荷矩陣A，B，C，D的對角陣。k=1，2，…，K，n=1，2，…，N。K×N個切片累積成N×J×K×I的四維形式，R為最佳組分?jǐn)?shù)，則其四線性平行因子標(biāo)量形式可表示為：

(12)

式中，r=1，2，…，R，并且n=1，2，…，N，j=1，2，…，J，i=1，2，…，I，k=1，2，…，K，式(12)即為四線性平行因子盲分離模型。

4 四線性平行因子的分解

kA+kB+kC+kD≥2R+3

(13)

則在寬松條件下，A，B，C，D可唯一確定，也就是X的K秩分解具有唯一可以確定，即四線性分解具有唯一性。其中，kA，kB，kC，kD分別是載荷矩陣A，B，C，D的秩，R表示四線性平行因子模型中的成分?jǐn)?shù)。

四線性平行因子用四線性交替最小二乘方法(quadrilinear alternating least square，QALS)完成四線性平行因子模型的擬合迭代，其基本思路是每一步更新一個矩陣，更新辦法是：對余下的矩陣，依據(jù)前一次估計的結(jié)果，利用最小二乘法來更新，該形式的更新迭代過程一致重復(fù)，直至算法收斂。相比于三線性迭代方法，該方法擬合精度更高，表現(xiàn)更平穩(wěn)，收斂速度更快。一般來講，迭代過程是通過初始化承載矩陣B，C，D，然后根據(jù)式(14)求出A，根據(jù)求得的A，由式(15)、式(16)、式(17)依次分別求出矩陣B，C，D。循環(huán)該過程直至算法收斂。

A=XI×JKN[(B⊙D⊙C)T]+

(14)

B=XJ×KLI[(A⊙D⊙C)T]+

(15)

C=XK×NIJ[(B⊙D⊙A)T]+

(16)

D=XN×IJK[(B⊙A⊙C)T]+

(17)

四線平行因子分析的步驟為：

1) 估計最佳組分?jǐn)?shù)

2) 初始化3個載荷矩陣B，C，D

3) 由式(14)迭代計算A

4) 依據(jù)迭代計算的A由式(15)迭代計算B，依次計算C，D。

5) 重復(fù)步驟3)和4)，到收斂至滿足收斂準(zhǔn)則。

該算法可通過優(yōu)化式(18)所示的迭代函數(shù)實現(xiàn)X的分解:

(18)

四線性平行因子迭代的停止準(zhǔn)則是：

(19)

其中：上標(biāo)符號“+”含義為Moore-Penrose偽逆；⊙代表矩陣的Khatri-Rao積；ε殘余誤差的平方和，m為迭代次數(shù)。一般認(rèn)為ε為1×10-6，當(dāng)滿足式(19)就可以認(rèn)為算法收斂。

5 仿真研究

模擬兩組調(diào)制信號s1(t)，s2(t)來產(chǎn)生數(shù)據(jù)仿真說明算法有效性。信號s1(t)由基頻為70 Hz，調(diào)頻為30 Hz的信號和一頻率為150 Hz的正弦信號疊加而成。信號s2(t)由基頻為120 Hz，調(diào)頻為10 Hz的信號和一頻率為145 Hz的正弦信號疊加而成。

(20)

2組源信號的采樣頻率為Fs=3 000 Hz，采樣點數(shù)為N=10 240，源信號s1(t)，s2(t)的時域波形圖如圖2，幅值譜圖如圖3。

圖2 源信號時域波形圖Fig.2 Time domain waveform of Source signal

圖3 源信號幅值譜圖Fig.3 Amplitude spectrum of Source signal

源信號幅值譜圖如圖3(a)以及圖3 (b)所示，圖3(a)中源信號s1(t)和圖3(b)源信號s2(t)的所顯示的頻率均為基頻頻率和調(diào)頻頻率的疊加，s1(t)其特征頻率分別為40 Hz、70 Hz、100 Hz、150 Hz。s2(t)的所顯示的頻率也為基頻頻率和調(diào)頻頻率的疊加，其特征頻率分別為120 Hz、110 Hz、130 Hz、145 Hz。

為了得到虛擬的觀測信號將兩路觀測信號混疊為一路，這里選擇一個隨機混合矩陣A并且添加噪聲干擾得到混疊信號時域波形圖如圖4，得到混疊信號幅值譜圖如圖5。由圖4以及圖5可知中兩路仿真信號的源信號時域波形圖完全混疊在了一起，并且兩源信幅值均在混合幅值譜上顯示，其特征頻率分別為40 Hz、70 Hz、100 Hz、110 Hz、120 Hz、145 Hz、150 Hz，兩源信號幅值譜圖完全混疊，并且無法辨識兩源信號的特征頻率，表明兩路觀測信號已經(jīng)完全混疊。

圖4 混疊信號時域波形圖Fig.4 Time domain waveform of aliased signal

圖5 混疊信號幅值譜圖Fig.5 Amplitude spectrum of aliased signal

為了比較三線性和四線性的迭代平穩(wěn)性和以及迭代收斂性，比較仿真信號滿足三線性和四線性平行因子迭代的停止準(zhǔn)則時三線性和四線性平行因子迭代收斂過程中迭代函數(shù)值的變化以及迭代結(jié)束時的迭代次數(shù)，比較結(jié)果如圖6所示。由圖6可知，迭代結(jié)束時四線性平行因子的迭代次數(shù)明顯低于三線性平行因子，并且在迭代過程中，四線性迭代函數(shù)的值明顯比三線性低，表明四線性相比于三線性平行因子收斂性更好，意味著四線性在的估計精度更高估計誤差更小。

圖6 迭代次數(shù)曲線Fig.6 Comparison of iteration times of different methods

用四線性平行因子盲分離方法分離混疊信號波形圖得到的時域波形圖，如圖7(a)和圖7(b)。用四線性平行因子盲分離方法分離混疊信號波形圖恢復(fù)得到的幅值譜圖，如圖8(a)與圖8(b)。用三線性平行因子分離得到的時域波形圖，如圖9(a)和圖9(b);得到幅值譜圖，如圖10(a)與圖10(b)。

圖7 恢復(fù)信號時域波形圖(四線性平行因子)Fig.7 Time-domain waveform of recovered signal (quadrilinear parallel factors)

圖8 恢復(fù)信號幅值譜圖(四線性平行因子)Fig.8 Altitude spectrum of recovered signal (quadrilinear parallel factors)

由圖8可知，用四線性平行因子盲分離方法能夠分離分離混疊信號，由圖7(a)可知，其主要特征頻率主要為40 Hz、70 Hz、100 Hz、150 Hz，其與源信號s1(t)的特征頻率基本一致，所以用四線性平行因子盲分離算法分離混疊信號其可以識別源信號s1(t)。由圖8(b)可知，特征頻率主要為110 Hz、120 Hz、130 Hz、145 Hz，其與源信號s2(t)的特征頻率基本一致，所以用四線性平行因子盲分離算法分離混疊信號其可以識別源信號s2(t)。由圖7與圖2以及圖4的對比(即四線性平行因子盲分離方法分離得到的時域波形圖與源信號時域波形圖的對比以及混疊信號時域波形圖對比)可以發(fā)現(xiàn)四線性平行因子盲分離方法很好地分離了混疊信號，分離得到的時域波形圖與源信號時域波形圖相似度很高。

圖9 恢復(fù)信號時域波形圖(三線性平行因子)Fig.9 Time-domain waveform of recovered signal (trilinear parallel factors)

圖10 恢復(fù)信號幅值譜圖(三線性平行因子)Fig.10 Altitude spectrum of recovered signal (trilinear parallel factors)

由圖8與圖10的對比(即四線性平行因子盲分離方法分離得到的幅值譜圖與三線性平行因子盲分離方法分離得到的幅值譜圖對比)可知三線在分離時效果沒有四線性分離效果好。由圖8與圖3以及圖5的對比(即四線性平行因子盲分離方法分離得到的幅值譜圖與源信號幅值譜圖以及混疊信號幅值譜圖對比)可知，四線性平行因子盲分離方法成功分離了混疊信號，并且分離得到的信號特征頻率與兩源信號特征頻率一致，表明四線性盲分離算法可以有效地分離混合信號得到源信號。

盲源分離中通常用性能指數(shù)以及相似系數(shù)來評判分離效果的好壞，這里選擇相似系數(shù)來分析評定。將混疊信號分別用三線性平行因子盲分離算法以及4-PARAFAC-BASS分離算法的分離，然后計算相似系數(shù)。式(21)是相似系數(shù)值的計算方法，式中cov(·)為方差。一般認(rèn)為恢復(fù)信號與源信號越相似，則相似系數(shù)|ρij|的值越接近為1。若|ρij|值越大則表示分離效果越好，|ρij|值越小分離效果越不好。

(21)

用ρij在信噪比為20 dB時，PARAFAC-BSS分離算法與4-PARAFAC-BASS分離算法的分離性能，如表1所示。

表1 PARAFAC-BSS與4-PARAFAC-BSS的分離性能Table 1 Comparison of separation performance between PARAFAC-BSS and 4-PARAFAC-BSS

由表1可知，三線性平行因子盲源分離(PARAFAC-BSS)得到的相似系數(shù)分別是0.667 7與0.304 1與四線性平行因子盲源分離(4-PARAFAC-BSS)得到的相似系數(shù)0.932 0與0.746 4相比，四線性得到的相似系數(shù)的值更接近于1，由此可知四線性的分離效果更好。并且三線性分離中s2(t)的相似系數(shù)的值較低，分離效果很不好。綜上可知，表明四線性盲分離算法可以有效地分離混合信號得到源信號。并且四線性平行因子盲分離算法得到的分離效果比三線性平行盲分離算法得到的分離效果更好，證實了所提算法的有效性。

6 實驗研究

用如圖11所示的實驗裝置，驗證所提四線性平行因子盲分離算法的有效性。用電火花加工技術(shù)分別在驅(qū)動端軸承內(nèi)圈中央位置以及軸承外圈中央位置加工點蝕故障，故障直徑為0.007英寸，深度為0.011英寸。將加工過的故障軸承重新裝入測試電機中，在電機負(fù)載為0馬力的工況條件下工作，電機轉(zhuǎn)速1 796 r/min，相應(yīng)的轉(zhuǎn)頻為fr=29.17 Hz，采樣頻率為12 000 Hz，采樣點數(shù)為24 000，軸承內(nèi)圈的故障頻率為104.31 Hz，軸承外圈的故障頻率為143.95 Hz。當(dāng)電機1單獨運行時用加速度傳感器采集到的振動信號數(shù)據(jù)，截取0.4 s的振動信號數(shù)據(jù)得到時域波形圖如圖11中(a)所示，當(dāng)電機2單獨運行時用加速度傳感器2采集振動信號數(shù)據(jù)，截取0.4 s的振動信號數(shù)據(jù)得到時域波形圖如圖11中(b)所示。

圖11 實驗裝置示意圖Fig.11 Test bed

當(dāng)兩電機分別單獨轉(zhuǎn)動時，采集數(shù)據(jù)，得到的兩源信號時域波形圖分別如圖12(a)與圖12(b)所示，幅值譜圖13(a)對應(yīng)特征頻率為電機轉(zhuǎn)頻29.17 Hz，軸承故障特征頻104.31 Hz，以及轉(zhuǎn)頻的二倍頻等以及故障特征頻率和轉(zhuǎn)頻的疊加影響頻率。幅值譜圖13(a)與圖13(b)所示分別對應(yīng)特征頻率為電機轉(zhuǎn)頻29.17 Hz，軸承故障特征頻率143.95 Hz，以及轉(zhuǎn)頻的二倍頻等以及故障特征頻率和轉(zhuǎn)頻的疊加影響頻率。

圖12 源信號時域波形圖Fig.12 Time-domain waveform of source signal

當(dāng)電機1和電機2同時運行時得到雙通道混合信號時域波形圖如圖14(a)與圖14(b)所示，得到的幅值譜圖如圖15(a)與圖15 (b)所示。

圖13 源信號幅值譜圖Fig.13 Altitude spectrum of source signal

圖14 混疊信號時域波形圖Fig.14 Time domain waveform of aliased signal

圖15 混疊信號幅值譜圖Fig.15 A ltitude spectrum of aliased signal

由圖14的混疊信號時域波形圖與圖12源信號時域波形圖對比可知，每一個通道采集得到的振動信號都是兩個電機同時轉(zhuǎn)動時信號的疊加。由圖13的混疊信號的幅值譜圖與圖13的源信號幅值譜圖對比可知，圖15中兩源信號的特征頻率均同時在混合信號的幅值譜上顯示，無法辨識各源信號的特征頻率，所以兩源信號完全混合。

用四線性平行因子盲分離算法分解得到的時域波形圖如圖16(a)與圖16 (b)示，以及得到的幅值譜圖如圖17(a)與圖17(b)所示。

圖16 恢復(fù)信號的時域波形圖(四線性平行因子)Fig.16 Time domain waveform of recovered signal separated by quadrilinear parallel factor

圖17 恢復(fù)信號幅值譜圖(四線性平行因子)Fig.17 Altitude spectrum of recovered signal separated by quadrilinear parallel factor

由圖16可知，由圖17所示的四線性平行因子盲分離方法能夠分離分離混疊信號。由圖17中(a)可知，特征頻率主要為29.17 Hz、58.34 Hz、85.61 Hz、143.95 Hz，其分別對應(yīng)著電機轉(zhuǎn)頻、轉(zhuǎn)頻的二倍頻、二倍頻與故障頻率的疊加頻率、以及故障特征頻率。由其故障特征頻率為143.95 Hz可知，該頻譜為軸承外圈的故障。由圖17中(b)可知，特征頻率主要為29.17 Hz、104.31 Hz、208.61 Hz，其分別對應(yīng)著電機轉(zhuǎn)頻、故障特征頻率以及兩倍故障特征頻率。由其故障特征頻率為104.31 Hz可知，該頻譜為軸承內(nèi)圈的故障。由圖17與圖12以及圖14的對比(即四線性平行因子盲分離方法分離得到的時域波形圖與源信號時域波形圖的對比以及混疊信號時域波形圖對比)可以發(fā)現(xiàn)，分離得到的時域波形圖與源信號時域波形圖有很高的相似度，分離效果良好。由圖17與圖13以及圖15的對比(即四線性平行因子盲分離方法分離得到的幅值譜圖與源信號幅值譜圖以及混疊信號幅值譜圖對比)可知，四線性平行因子盲分離方法成功分離了混疊信號，并且分離得到的信號特征頻率與兩源信號特征頻率一致，但是由于盲分離本身模糊性的特性，會導(dǎo)致了盲分離得到的幅值譜圖排序次序出現(xiàn)顛倒或者互換，這種排序及波形幅值的變換并不影響分離結(jié)果的準(zhǔn)確性。從圖17可以看出分離得到的幅值譜圖除了排序及波形幅值有些變換外，恢復(fù)以及分離效果良好。四線性平行因子盲分離方法有效并且分離效果良好。

7 結(jié)論

本文提出基于四線性平行因子的機械故障盲分離方法，將三線性平行因子擴展到四線性，彌補了三線性建模時只對振動信號部分信息建模的缺陷，使數(shù)據(jù)包含的信息更完整。該方法在盲分離無需通道信息和源信號等先驗知識的條件下能對分離出未知源信號，并且將盲分離與四線性平行因子在寬松條件下分解唯一性結(jié)合，在振動信號特征提取和信號分離時能充分發(fā)揮2種算法的優(yōu)勢。并且該方法使用四線性最小二乘法進行迭代擬合，相比于三線性平行因子迭代速度更快，迭代誤差更小，迭代精確度更高。最后利用仿真和實驗證實了所提的四線性平行因子的機械故障盲分離方法可以很好地提取故障特征并分離混合源信號。