基于改進(jìn)S變換和ICA的相關(guān)源分離方法*

韋成龍 ， 周以齊 ， 李 瑞 ， 于 剛

(1.山東大學(xué)機(jī)械工程學(xué)院 濟(jì)南，250061) (2. 濟(jì)南大學(xué)自動化與電氣工程學(xué)院 濟(jì)南，250061)

引 言

機(jī)械系統(tǒng)運(yùn)行時各個部件同時工作，由于受到不同振動源的激勵，測試采集的振動信號往往是多源混合的結(jié)果。如何從復(fù)雜的觀測信號中快速準(zhǔn)確地分離出源信號，對于振源識別和貢獻(xiàn)量估計具有重要意義。獨(dú)立成分分析作為盲源分離的常用方法，在源未知的情況下，僅利用觀測信號和先驗(yàn)知識分離和提取信號源的過程，廣泛應(yīng)用于噪聲控制、振源識別和故障診斷等領(lǐng)域[1-2]。

ICA以源信號相互統(tǒng)計獨(dú)立作為假設(shè)條件，基于非高斯性最大化原理分離觀測變量，而在實(shí)際工程應(yīng)用中，嚴(yán)格的統(tǒng)計獨(dú)立往往過于理想化。尤其對于動力機(jī)械來說，振源眾多且傳遞路徑復(fù)雜，發(fā)動機(jī)作為主要的動力源，振動信號具有明顯的沖擊和時變特征，且各部件級聯(lián)耦合，觀測信號中獨(dú)立成分和相關(guān)成分混合在一起，難以獲得理想的分離效果[3-4]。當(dāng)源信號存在統(tǒng)計相關(guān)成分時，分離獲得的估計源與真實(shí)源之間將相差一個與源向量協(xié)方差有關(guān)的因子向量[5]。

在解決相關(guān)性源信號分離的問題上，國內(nèi)外學(xué)者開展了相關(guān)研究，主要方法有子帶分解[6-7]、核典型相關(guān)[8]和稀疏成分分析[9-10]等。文獻(xiàn)[7]提出基于小波包分解的相關(guān)源分離方法，利用互信息作為獨(dú)立性測度重構(gòu)觀測信號，獲得分離矩陣，但分解性能受先驗(yàn)知識和小波分解層數(shù)的影響。文獻(xiàn)[8]采用核典型相關(guān)分析的方法將非線性信號映射到核特征空間中，將非線性問題轉(zhuǎn)化為線性，解決非線性相關(guān)信號的分離，而核化過程并不能有效消除相關(guān)項(xiàng)。文獻(xiàn)[9]利用稀疏信號延混合矩陣列方向向量的線性聚類特征，能夠?qū)崿F(xiàn)源數(shù)欠定條件下的估計，當(dāng)信源的稀疏度不夠時，分離誤差較大。

筆者借鑒子帶分解和稀疏分量分析的思想，以相關(guān)源中部分成分統(tǒng)計獨(dú)立為前提假設(shè)，提出了一種基于改進(jìn)S變換和ICA的相關(guān)源分離方法。在S變換[11-12]中引入?yún)?shù)可調(diào)的尺度因子，將時域觀測信號擴(kuò)展到時頻域，利用混合觀測信號在時頻域中實(shí)部和虛部的方向向量，檢測并剔除源信號中的相關(guān)成分，避免了獨(dú)立子帶方法中信號分解層數(shù)的確定和互信息在同頻成分檢測上的失效問題，通過重構(gòu)觀測信號以滿足盲源分離的統(tǒng)計獨(dú)立條件，進(jìn)而采用負(fù)熵最大化的獨(dú)立成分分析獲得分離矩陣，實(shí)現(xiàn)相關(guān)源的有效分離。

1 分離模型

盲源分離是在混合系統(tǒng)未知的情況下僅利用觀測信號恢復(fù)源信號的過程。假設(shè)正定或超定條件下，即n個源信號瞬時混合被m個傳感器接收(n≤m)，將觀測信號中的多變量數(shù)據(jù)簡化為線性混合問題，其數(shù)學(xué)模型可表示為

x(t)=As(t)+n(t)

(1)

其中：A∈Rm×n為列滿秩混合矩陣；x=(x1,x2,…,xm)T為觀測變量；s=(s1,s2,…,sn)T為n維未知的信號源；n(t)為附加噪聲。

觀測變量x(t)作為信號源s(t)的線性加權(quán)，估計源信號y(t)就轉(zhuǎn)化為尋找分離矩陣W的過程，即

y(t)=Wx(t)

(2)

當(dāng)激勵源存在同頻分量時，信號之間統(tǒng)計相關(guān)，無法采用以獨(dú)立源作為假設(shè)前提的標(biāo)準(zhǔn)ICA方法直接分離出來。假設(shè)相關(guān)源由獨(dú)立成分sD(t)和非獨(dú)立成分sID(t)兩部分組成

s(t)=sD(t)+sID(t)

(3)

觀測信號可看作獨(dú)立成分和相關(guān)成分的線性疊加，在瞬時混合模型下，各成分的混合系數(shù)相同，分離矩陣W可通過部分獨(dú)立成分的混合信息獲得。

yD(t)=Wxrec(t)=WAsD(t)

(4)

其中：xrec(t)為獨(dú)立成分重構(gòu)的觀測信號；yD(t)為剔除相關(guān)項(xiàng)后的估計源。

筆者通過時頻預(yù)處理方法，從混合變量中剔除源信號的相關(guān)成分，對剩余信號進(jìn)行重構(gòu)，從而滿足統(tǒng)計獨(dú)立的條件。對重構(gòu)后的信號進(jìn)行盲分離，估計出分離矩陣W，由式(2)恢復(fù)相關(guān)源信號y(t)。

2 算法原理

2.1 改進(jìn)S變換

S變換作為一種短時傅里葉變換的可變窗擴(kuò)展和連續(xù)小波的相位修正，能夠有效克服窗寬固定的缺點(diǎn)，時頻變換后各分量的相位與原始信號保持直接聯(lián)系，對工程測試中非平穩(wěn)信號具有較好的時頻特性[13-15]。標(biāo)準(zhǔn)S變換是在短時傅里葉變換的基礎(chǔ)上引入時寬與頻率倒數(shù)成正比關(guān)系的高斯窗，定義為

其中：ST(τ,f)為信號x(t)的S變換；g(τ-t,f)為窗函數(shù)；τ為平移因子，用于控制高斯窗在時間軸t上的位置；f為信號頻率。

高斯窗函數(shù)中標(biāo)準(zhǔn)差被改造為關(guān)于頻率f的倒數(shù)，然而標(biāo)準(zhǔn)S變換的窗寬變化趨勢固定，時頻譜唯一。為了提高其適應(yīng)性，修正窗函數(shù)模式固定的缺陷，引入尺度因子λ,α,β，得到改進(jìn)S變換的表達(dá)式為

(7)

其逆變換(inverse modified S-transform，簡稱 IMST)表示為

(8)

圖1為引入尺度因子后，窗函數(shù)半高時寬隨頻率的變化曲線。當(dāng)λ=1,α=1,β=0時，可以看作標(biāo)準(zhǔn)S變換。當(dāng)0<α<1時，時窗寬度的減小速率降低，高頻部分頻率分辨率提高，而α>1時，通常用于捕捉高頻瞬態(tài)信號。當(dāng)0<λ<1時，頻率分辨率提高，窗寬呈非線性變化，趨于平緩，λ>1時，則相反。β通常取正值，用于提高低頻區(qū)域的時間分辨率。λ,α,β均具有控制窗口寬度和變化率的作用，可針對實(shí)際需求優(yōu)化三者參數(shù)，以獲得最優(yōu)時頻分辨率。

圖1 窗函數(shù)半高時寬隨頻率變化曲線Fig.1 The curve of temporal full width at half maximum versus frequency

2.2 相關(guān)成分檢測

觀測信號x(t)經(jīng)MST得到一個關(guān)于時間t和頻率f的二維復(fù)時頻矩陣x(t,f),aj=[a1j,a2j,…,amj]T為混合矩陣A的第j列向量，則瞬時混合模型在時頻域中的表達(dá)式為

(9)

僅存在兩個激勵源s1,s2時，式(9)可寫作

x(t,f)=a1s1(t,f)+a2s2(t,f)

(10)

滿足統(tǒng)計獨(dú)立的條件下，混合矩陣中的每個時頻點(diǎn)僅為其中一個獨(dú)立源的響應(yīng)。此時，假設(shè)在時間tp和頻率fq處，s1(tp,fq)≠0，s2(tp,fq)=0，式(10)可表示為

x(tp,fq)=a1s1(tp,fq)

(11)

其實(shí)部和虛部變換分別為

Re{x(tp,fq)}=a1Re{s1(tp,fq)}

(12)

Im{x(tp,fq)}=a1Im{s1(tp,fq)}

(13)

由式(12)，(13)可知，x(tp,fq)實(shí)部和虛部的方向向量與混合矩陣的列向量方向一致。當(dāng)兩激勵源在時間tp和頻率fq處頻率相同，即存在相關(guān)成分時，源信號s1(tp,fq)≠0，s2(tp,fq)≠0，對應(yīng)觀測信號在時頻域中實(shí)部和虛部變換為

Re{x(tp,fq)}=a1Re{s1(tp,fq)}+a2Re{s2(tp,fq)}

(14)

Im{x(tp,fq)}=a1Im{s1(tp,fq)}+a2Im{s2(tp,fq)}

(15)

此時，只有在式(16)成立，相關(guān)成分的相位相差0°或180°的情況下，x(tp,fq)實(shí)部和虛部的列向量方向一致。

(16)

實(shí)際激勵源中，相關(guān)成分的相位往往是隨機(jī)的，同相位的概率相當(dāng)?shù)停鶎?yīng)的觀測信號在復(fù)時頻矩陣中的實(shí)部和虛部向量會存在一個角度差Δθ。將其擴(kuò)展為m×n維混合信號，時頻域中實(shí)部和虛部向量之間的夾角可表示為

(17)

式(17)表明，當(dāng)cosΔθ=1時，實(shí)部和虛部向量方向一致，觀測信號中對應(yīng)成分獨(dú)立，而cosΔθ<1時，對應(yīng)成分相關(guān)。考慮到工程測試中測點(diǎn)位置和噪聲的影響，很難嚴(yán)格滿足條件，通過設(shè)置閾值ε來降低其影響，ε取接近于1的值，當(dāng)cosΔθ<ε時，檢測為相關(guān)成分并將其剔除。閾值的確定需要根據(jù)實(shí)測數(shù)據(jù)來選擇，至少保留10%以上的數(shù)據(jù)用來保證分離矩陣的估計。

圖2 相關(guān)源分離算法流程圖Fig.2 Correlated source separation flow chart

2.3 相關(guān)源估計

基于MST和ICA的相關(guān)源分離算法流程如圖2所示。分離出相關(guān)成分后，利用IMST重構(gòu)剩余信號，從而滿足ICA對信號源統(tǒng)計獨(dú)立的要求，然后以最大化負(fù)熵作為獨(dú)立性判據(jù)，采用快速固定點(diǎn)算法[1](fast fixed-point independent component analysis，簡稱 FastICA)獲得分離矩陣。負(fù)熵作為一種非高斯性的度量，較峭度具有更好的魯棒性和統(tǒng)計特性，以非二次函數(shù)近似的形式估計負(fù)熵[16]

(18)

其中：y為標(biāo)準(zhǔn)化輸入變量；υ為具有零均值和單位方差的高斯隨機(jī)變量；G為非二次函數(shù)。

采用漸進(jìn)方差的形式，選取最優(yōu)非二次函數(shù)G(y)

(19)

其中：1≤α≤2時，取G(y)=log2cosh(y)；α<1，即超高斯程度很高時，取G(y)=-exp(-y2/2)。

通過極大化E{G(WTz)}獲得負(fù)熵的最大近似值，根據(jù)牛頓迭代更新分離矩陣W，如式(20)所示，每輪迭代對W進(jìn)行正交化和標(biāo)準(zhǔn)化，直到收斂為止。

W+=E{zG(WTz)}-E{zG′(WTz)}W

(20)

(21)

3 仿真分析

為了驗(yàn)證提出方法的有效性，構(gòu)造4個仿真信號模擬機(jī)械振源。s1，s2分別為20 Hz和120 Hz諧波信號，模擬中低頻簡諧振動，s3為調(diào)幅信號，模擬幅值調(diào)制特征；s4為指數(shù)衰減信號，模擬周期性沖擊。在s1～s4中分別加入相位隨機(jī)的40 Hz諧波成分sd。設(shè)采樣頻率為1 024 Hz，采樣時間為1 s，信號函數(shù)為

(22)

圖3為加入40 Hz同頻相關(guān)成分后源信號的時域波形。采用皮爾森相關(guān)系數(shù)衡量信號之間相似程度，假設(shè)兩個隨機(jī)變量x和y，相關(guān)系數(shù)ρ定義為

圖3 源信號時域波形Fig.3 Waveform of source signals

(23)

ρ的絕對值越接近于1，說明兩變量之間波形的相似度越高，相關(guān)性也越強(qiáng)。

計算仿真源之間相關(guān)系數(shù)如表1所示。加入同頻成分后，源信號呈現(xiàn)不同程度的相關(guān)性。隨機(jī)生成(0,1)均勻分布的4×4混合矩陣A，以線性疊加的方式混合相關(guān)源得到4個觀測信號x1～x4，其時域波形如圖4所示。

表1 源信號相關(guān)系數(shù)Tab.1 Source signal correlation coefficient

圖4 觀測信號時域波形Fig.4 Waveform of observed signals

為使待分離信號滿足統(tǒng)計獨(dú)立條件，去除相關(guān)成分的干擾，首先利用MST提取時頻系數(shù)，取尺度參數(shù)α=0.3，β=5，λ=1，計算觀測信號在各時頻點(diǎn)上實(shí)部和虛部向量的夾角Δθ，設(shè)置閾值ε=0.9，將cosΔθ<ε成分置零，再經(jīng)IMST重構(gòu)剩余信號。圖5,6分別為觀測信號x1和重構(gòu)信號xrec1的MST時頻譜。對比看出，經(jīng)本方法處理后，觀測信號中的40 Hz相關(guān)成分得到有效去除。然后采用ICA分離重構(gòu)信號，利用分離矩陣得到估計源y1～y4。

圖5 觀測信號x1時頻譜Fig.5 Time-frequency spectrum of observed signal x1

圖6 重構(gòu)信號xrec1時頻譜Fig.6 Time-frequency spectrum of reconstruction signal xxec1

圖7 分離信號獨(dú)立分量時域波形Fig.7 Waveform of separated independent components

圖8 估計源時域波形Fig.8 Waveform of estimated source signals

圖7為分離重構(gòu)信號得到獨(dú)立分量yD1～yD4的時域波形。可以看出，分離源中20 Hz和120 Hz的中低頻諧波成分，調(diào)幅成分以及瞬態(tài)沖擊成分被有效分離出來，可用于相關(guān)源信號中獨(dú)立特征的提取。圖8為本研究方法得到的估計源y1～y4的時域波形。對比圖3，波形上與加入40 Hz相關(guān)成分后的仿真源基本一致，并計算分離信號與仿真源之間的相關(guān)系數(shù)，如表2所示。y1，y2，y3，y4分別與仿真源s4，s2，s3，s1的相關(guān)系數(shù)為1，說明在無噪條件下，幾乎可以實(shí)現(xiàn)相關(guān)源信號的完全分離與估計。

表2 分離信號相關(guān)系數(shù)Tab.2 Separated signal correlation coefficient

相關(guān)系數(shù)僅用來衡量兩兩變量間的相似程度，為了進(jìn)一步評價算法的分離性能，采用Amari性能指數(shù)[18](performance index，簡稱 PI)來度量混合矩陣A和分離矩陣W之間的差異，定義為

(24)

其中：gij作為全局矩陣G=WA的第(i,j)個元素，PI越接近于零，說明算法分離性能越好，當(dāng)超過0.2時，可認(rèn)為分離失效。

在x1～x4中加入信噪比為5～50 dB的高斯噪聲，對比本分離方法與直接采用FastICA分離時的PI值，如圖9所示。可以看出，隨著信噪比的降低，算法的分解性能變差，而本研究方法的PI指數(shù)明顯小于直接使用FastICA的分解結(jié)果，在信噪比20 dB以上的有噪條件下，對相關(guān)源混合信號仍能夠取得較好的分離效果。

圖9 Amari性能指數(shù)對比Fig.9 Comparison of Amari performance index

4 實(shí)例分析

工程機(jī)械中觀測變量通常為多激勵源的綜合響應(yīng)，以某型挖掘機(jī)為測試對象，發(fā)動機(jī)作為主要動力源，振動信號經(jīng)4個減振裝置傳遞到車架上并混合，將每個減振點(diǎn)看作一個激勵，利用車架振動信號估計振源。測試采用8通道NI數(shù)據(jù)采集系統(tǒng)和壓電式加速度振動傳感器，在空擋最大油門工況下，同時測量四點(diǎn)激勵與上車架的振動信號。發(fā)動機(jī)轉(zhuǎn)速為2 280 r/min，采樣頻率為5 120 Hz,實(shí)驗(yàn)樣車和上車架傳感器測點(diǎn)布置如圖10(a)和圖10(b)所示。

圖10 挖掘機(jī)測試圖Fig.10 Experiment of excavator

由于車架振動來源于發(fā)動機(jī)和減振器的共同作用，激勵信號之間存在不同程度的相關(guān)性，四點(diǎn)減振處振動信號的相關(guān)系數(shù)最大為0.36。利用本研究方法分離車架振動信號，取尺度參數(shù)α=0.3，β=56，λ=0.8，設(shè)置閾值ε=0.9，去除相關(guān)項(xiàng)后重構(gòu)獨(dú)立混合信號，最后利用ICA獲得分離矩陣W,估計出振源信號，確定各個振源對車架振動的貢獻(xiàn)量。

圖11 上車架振動信號Fig.11 Upper frame vibration signals

圖12 分離信號Fig.12 Separated signals

圖13 分離信號時頻譜Fig.13 Time-frequency spectrum of separated signals

圖11，圖12分別為車架振動信號x1～x4和本研究方法分離得到的估計源y1～y4的時域波形和頻譜圖。圖13為分離信號的MST時頻譜。可以看出，4個分離信號的能量主要在1 000 Hz以下，以76 Hz和380 Hz頻率成分為主要峰值，對應(yīng)發(fā)動機(jī)發(fā)火基頻和10階轉(zhuǎn)頻，y1中600 Hz附近周期性沖擊成分也較為明顯，應(yīng)為缸體內(nèi)氣體燃燒產(chǎn)生的振蕩沖擊引起。

分離信號與激勵信號之間的相關(guān)系數(shù)如表3所示。估計源y1～y4分別與左后減振、左前減振、右前減振、右后減振處的振動信號相關(guān)系數(shù)達(dá)到0.8以上，為顯著相關(guān)，說明分離出來信號與4個激勵點(diǎn)處振動信號的波形基本一致。本研究方法能夠從車架混合信號中成功分離出激勵源，對于工程測試信號，依然能夠取得較好的分離效果。

利用分離矩陣與估計信號定量計算各個振源對車架振動的貢獻(xiàn)比，如表4所示。可以看出，右側(cè)前后減振處激勵對車架振動貢獻(xiàn)較大，其中，右前減振作為第1主振源，貢獻(xiàn)比為40.4%。結(jié)合車架振動以及分離信號的譜分析結(jié)果推測，右側(cè)車架的380 Hz峰值頻率成分應(yīng)主要由右前減振處振動激勵引起，可通過針對性的調(diào)節(jié)右前減振裝置、優(yōu)化車架結(jié)構(gòu)，來達(dá)到改善車架振動的目的。

表3 分離信號相關(guān)系數(shù)Tab.3 Separated signal correlation coefficient

表4 振源貢獻(xiàn)百分比Tab.4 Percentage of vibration source contribution

5 結(jié)束語

針對傳統(tǒng)盲源分離方法對于統(tǒng)計獨(dú)立條件的限制，提出了一種基于改進(jìn)S變換和獨(dú)立成分分析的相關(guān)源分離方法，用來解決工程應(yīng)用中相關(guān)性機(jī)械振源在混合觀測信號中的分離和估計問題。該方法在線性瞬時混合盲源分離模型的基礎(chǔ)上，針對獨(dú)立成分和相關(guān)成分共存的信號源，以相關(guān)成分相位隨機(jī)為假設(shè)，利用其在時頻域中的方向向量，有效剔除混合信號中的相關(guān)分量，避免了相關(guān)成分對分離結(jié)果的影響。在時頻化處理中，以S變換為基礎(chǔ)引入窗寬可調(diào)的尺度因子，改善時頻分辨率固定的缺陷，提高了復(fù)雜工程測試信號的適應(yīng)性。通過仿真和實(shí)測分析驗(yàn)證了提出方法在相關(guān)源分離中的有效性，分離性能優(yōu)于未對相關(guān)項(xiàng)進(jìn)行預(yù)處理的盲源分離結(jié)果，通過閾值的合理設(shè)置，能夠從混合信號準(zhǔn)確分離出相關(guān)源信號，得到各個振源對混合信號的貢獻(xiàn)比。該方法適用于正定、超定條件下含有交叉頻率成分的相關(guān)源分離，在信源識別和故障診斷等領(lǐng)域具有一定的工程應(yīng)用價值。