基于改進S變換和ICA的相關源分離方法*

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

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

引 言

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

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

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

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

1 分離模型

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

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)的線性加權，估計源信號y(t)就轉化為尋找分離矩陣W的過程，即

y(t)=Wx(t)

(2)

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

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

(3)

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

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

(4)

其中：xrec(t)為獨立成分重構的觀測信號；yD(t)為剔除相關項后的估計源。

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

2 算法原理

2.1 改進S變換

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

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

高斯窗函數中標準差被改造為關于頻率f的倒數，然而標準S變換的窗寬變化趨勢固定，時頻譜唯一。為了提高其適應性，修正窗函數模式固定的缺陷，引入尺度因子λ,α,β，得到改進S變換的表達式為

(7)

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

(8)

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

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

2.2 相關成分檢測

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

(9)

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

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

(10)

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

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

(11)

其實部和虛部變換分別為

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

(12)

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

(13)

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

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)成立，相關成分的相位相差0°或180°的情況下，x(tp,fq)實部和虛部的列向量方向一致。

(16)

實際激勵源中，相關成分的相位往往是隨機的，同相位的概率相當低，所對應的觀測信號在復時頻矩陣中的實部和虛部向量會存在一個角度差Δθ。將其擴展為m×n維混合信號，時頻域中實部和虛部向量之間的夾角可表示為

(17)

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

圖2 相關源分離算法流程圖Fig.2 Correlated source separation flow chart

2.3 相關源估計

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

(18)

其中：y為標準化輸入變量；υ為具有零均值和單位方差的高斯隨機變量；G為非二次函數。

采用漸進方差的形式，選取最優非二次函數G(y)

(19)

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

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

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

(20)

(21)

3 仿真分析

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

(22)

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

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

(23)

ρ的絕對值越接近于1，說明兩變量之間波形的相似度越高，相關性也越強。

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

表1 源信號相關系數Tab.1 Source signal correlation coefficient

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

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

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

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

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

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

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

表2 分離信號相關系數Tab.2 Separated signal correlation coefficient

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

(24)

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

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

圖9 Amari性能指數對比Fig.9 Comparison of Amari performance index

4 實例分析

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

圖10 挖掘機測試圖Fig.10 Experiment of excavator

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

圖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頻率成分為主要峰值，對應發動機發火基頻和10階轉頻，y1中600 Hz附近周期性沖擊成分也較為明顯，應為缸體內氣體燃燒產生的振蕩沖擊引起。

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

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

表3 分離信號相關系數Tab.3 Separated signal correlation coefficient

表4 振源貢獻百分比Tab.4 Percentage of vibration source contribution

5 結束語

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