基于零空間表示和最大似然的欠定盲源分離方法

王榮杰，詹宜巨，周海峰

(1.中山大學(xué) 信息科學(xué)與技術(shù)學(xué)院，廣東 廣州 510006；2.中山大學(xué) 工學(xué)院，廣東 廣州 510006；3.集美大學(xué) 輪機(jī)工程學(xué)院，福建 廈門 361021)

1 引言

近年來，鑒于盲源分離優(yōu)越的假設(shè)前提，其模型被廣泛用于數(shù)字通信、語音、圖像處理、船舶工程、機(jī)器人導(dǎo)航和生物醫(yī)學(xué)信號處理等領(lǐng)域[1～6]。所謂的盲源分離，就是在源信號和混合系統(tǒng)(或傳輸通道)等未知的情況下，僅根據(jù)源信號的統(tǒng)計特性，由觀測到的混疊信號恢復(fù)出源信號。根據(jù)觀測信號和源信號的個數(shù)，盲源分離可分為超(正)定盲源分離和欠定盲源分離。超正定盲源分離有很多成熟的方法，如 ICA獨(dú)立分量分析[7,8]、修正的Stone法[9]、階統(tǒng)計代數(shù)法[10]等。但這些方法不適用于觀測信號的個數(shù)少于源信號的個數(shù)的情況，欠定盲源分離(UBSS, underdetermined blind source separation)是本文的研究重點(diǎn)。傳統(tǒng)的欠定盲源分離算法主要可分為2類：一類是利用源信號的稀疏性來處理，文獻(xiàn)[11]提出利用聚類技術(shù)和最小化 l1范數(shù)相結(jié)合的方法來恢復(fù)時域稀疏的源信號，文獻(xiàn)[12,13]提出利用時頻變換方法將非稀疏信號變換到時頻域，以達(dá)稀疏的目的，它們認(rèn)為在每個時頻區(qū)域內(nèi)只有一個源信號不為零；這些方法都選擇服從超高斯分布的語音為源信號。另一類是利用廣義分布模型作為源信號的概率密度的貝葉斯估計法[14,15]，這類方法的主要缺點(diǎn)是計算復(fù)雜，在一定程度上降低了源信號的恢復(fù)質(zhì)量。這些傳統(tǒng)的算法都假設(shè)源數(shù)為已知。本文提出一種源數(shù)未知的欠定情況下任意分布源信號的盲分離方法，該方法首先利用S變換和聚類技術(shù)相結(jié)合來估算源數(shù)和混疊矩陣，然后將源信號以零空間形式表示，最后通過最大似然 (ML,maxima likelihood) 法估計關(guān)于它的后驗(yàn)概率以達(dá)到恢復(fù)源信號的目的。

2 問題的描述

假設(shè)n個彼此相互獨(dú)立的未知源信號s(t)=[s1(t),s2(t), …, sn(t)]T，通過一未知瞬時線性混合系統(tǒng)后，得到m個觀測信號x(t)=[x1(t), x2(t), …, xm(t)]T。觀測信號x(t)與源信號s(t)的關(guān)系可由式(1)表示。

其中，A∈Rm×n為混合矩陣，它反映了混合系統(tǒng)或信道的傳輸特性，要求它為行滿秩，t=0, … ,N-1為時域采樣點(diǎn)。為了便于分析，在推導(dǎo)過程中將s(t)和x(t)分別記為s和x。

在m ≥ n的情況下，給定A，源信號s可由式(2)估計得到。

其中，A*為A的廣義逆矩陣，A*= AT(AAT)。

當(dāng)m＜ n時，即為欠定情況，既便A已知，對于源信號s的恢復(fù)也不是唯一，只能通過估計方法估計出s的最優(yōu)估計值。

3 源數(shù)和混疊矩陣的估計

在第4節(jié)中分析的前提條件是假設(shè)混疊矩陣A和n為已知，特別在傳統(tǒng)算法中的n都設(shè)為已知，但在實(shí)際中A和n都是未知的。本文采用了一種基于S變換的單源測試(SSD, single source detention)方法，首先以此為準(zhǔn)則選擇出時頻域的單源觀測信號，然后采用聚類技術(shù)分別估計源信號個數(shù)n和混疊矩陣A。

3.1 基于S變換的單源測試方法

S變換是Stockwell于1996年提出的一種可逆的局部時頻分析方法，其思想是對連續(xù)小波變換和短時傅里葉變換的發(fā)展。它不同于短時傅里葉變換之處在于采用的高斯窗口的高度和寬度隨頻率而變化，這樣就克服了短時傅里葉變換窗口高度和寬度固定的缺陷[16,17]。觀測信號的 S變換可由式(3)和式(4)計算得到。

式(3)和式(4)中，i=1,…,m，代表時間的 t =0, …,N-1，代表頻率的k=0, …, N-1；式(4)中的(·)為式(3)中信號的傅里葉變換。

表1 不同單源信號測試方法計算復(fù)雜度的比較

式(5)中，║·║，Re[·]和 Im[·]分別為 Frobenius 范數(shù)，取實(shí)部和取虛部運(yùn)算。如果時頻點(diǎn)(t, k)為單源點(diǎn)時，且存在(t, k) ≠ 0(i=1, …, m)，那么 XS(t, k)中所有的元素均由某一個相同的源信號與混合矩陣 A中的元素相乘獲得，因此為純實(shí)數(shù)向量，即；如果(t, k)不為單源點(diǎn)時，那么為復(fù)數(shù)向量，而不為零且為一正數(shù)。式(5)中的ε為一個正數(shù)。由于具有單源特性的觀測信號XS(t,k)的實(shí)部或虛部與對應(yīng)的A中的列向量之間具有相同的絕對方向[19]，因此本文利用與傳統(tǒng)的基于稀疏特性的欠定混合矩陣估計算法中的聚類技術(shù)對具有單源特性的觀測信號進(jìn)行聚類估計混合矩陣，同時源信號的個數(shù)也可利用單源信號的特性和聚類技術(shù)來估計；式(5)中的Ψs為具有單源特性觀測信號XS(t,k)實(shí)部的集合。基于S變換的單源測試方法不僅克服了其他傳統(tǒng)單源測試方法中的時頻變換窗函數(shù)和長度重疊選擇盲目性的缺陷；同時，與傳統(tǒng)單源測試方法相比較，本文方法的計算簡單，通過表1比較了3種不同單源測試方法的計算復(fù)雜度可以證明。

表1中時頻變換的點(diǎn)數(shù)均設(shè)為N；利用表中的3種方法將每一個單源性觀測信號變?yōu)橛糜诠浪慊旌暇仃嚨木垲愒仨?xiàng)還需如下計算量：本文的方法只需 m次取實(shí)部運(yùn)算；而文獻(xiàn)[13]的方法需3m+2次乘法運(yùn)算、3(m-1)次加法運(yùn)算和 2m 次取實(shí)(虛)部運(yùn)算；文獻(xiàn)[19]的方法需m+1次乘法運(yùn)算、(m-1)次加法運(yùn)算和m次取實(shí)(虛)部運(yùn)算。

3.2 源數(shù)和混疊矩陣的估計

由于源數(shù)是未知的，并且欠定情況下的源數(shù)也不能用PCA或SVD進(jìn)行估算。為了估算源信號的個數(shù)，本節(jié)采用一種基于聚類驗(yàn)證技術(shù)確定源數(shù)，該方法根據(jù)不同的聚類數(shù)(聚類數(shù)c=2,…, cmax)對式(5)中的單源性元素進(jìn)行模糊 c均值聚類(fuzzy c-means clustering)，并以式(6)為準(zhǔn)則對不同c的聚類結(jié)果進(jìn)行評估得到最優(yōu)聚類數(shù)，即為源數(shù)。

式(6)中的 Scat(c)和 Sep(c)不僅與聚類數(shù) c有關(guān)，還與單源觀測信號集合 Ψs中的元素和聚類中心相關(guān)，它的具體計算可參考文獻(xiàn)[20]。Scat(c)表示聚類數(shù)為c的聚類緊湊性，當(dāng)聚類緊湊性越好時Scat(c)的值將越小，它的值在0～1之間；Sep(c)用于描述聚類數(shù)為c時聚類中心分布情況，稱為聚類分離度；當(dāng)聚類中心分布越合理時，它的值也將越小。因此，最小化式(6)，便可獲得 Ψs的最佳聚類數(shù)，即確定源信號個數(shù)。當(dāng)源數(shù)確定后，混疊矩陣將根據(jù)最優(yōu)聚類下的元素來估算。設(shè)ai為A的第i列向量，具體由式(7)估算得到。

4 源信號的恢復(fù)

在A和n都已知的欠定情況下，盲源分離算法就是要從x得到s的最優(yōu)估計。本節(jié)在假設(shè)A和n已知情況下，首先介紹了源信號的零空間表示[21]，在此基礎(chǔ)上采用了一種基于ML估計的方法來恢復(fù)源信號。

4.1 源信號的零空間表示

為了便于分析，先設(shè)式(1)中的A由一個m階的對角陣Λ和m×(n-m)維的0矩陣組成，記Λ=diag[λ1, …, λm]且 λi＜ λi-1, λi( i=1, … ,m)＞0，那么式(1)可改寫成式(8)。

由式(8)可知，s的前m個值可由式(9)得到。

其余的(n-m)個s用一個(n-m)維列向量的r表示，r =[r1,…, rn-m]T。由此式(8)中s的解可用式(10)來描述。

當(dāng)混合矩陣A不滿足式(8)中的特殊形式時，它的s解也就不能直接用式(10)來描述。但由于A為行滿秩，它的SVD奇異值分解中含m個不為零的特征值，可寫成A=U(Λ, 0)V，將其代入式(1)，并根據(jù)式(10)，可得到s的解如式(11)所示。

由式(12)可知，當(dāng) A已知，那么估計了(n-m)個的r就相當(dāng)于估計n個的s。

4.2 源信號的恢復(fù)

為了能同時分離服從任意分布的源信號，本文選取具有對稱單峰分布特性的廣義高斯模型 (GGM,generalized Gaussian model)[18,22]作為源信號的概率密度，式(13)為它的通用數(shù)學(xué)表達(dá)式。

其中，α為信號z的方差，()Γ·為Gamma函數(shù)；調(diào)節(jié)β的值可得出不同的分布函數(shù)模型，β=2表示標(biāo)準(zhǔn)高斯分布，β＜2時為超高斯分布，β＞2時為亞高斯分布。

設(shè) X=[ x(0), …, x(N-1)]，R=[r(0), … ,r(N-1)]，S=[ s(0), …, s(N-1)]，q=[α1, …, αn, β1, …,βn]T，并記S的概率密度為f(S|q)。那么可得X和R的聯(lián)合分布函數(shù)如式(14)所示。

為了估計每個源信號分布函數(shù)模型中的α和β，構(gòu)造關(guān)于參數(shù)向量q的先驗(yàn)似然函數(shù)如式(15)所示。

由最大似然估計法可知

即

式(14)～式(17)中，由于A獨(dú)于立于q，且分布p(X|A,q)和 p(A)是非負(fù)的，因此 q的最大似然估計可簡化為式(18)。

式(18)中估計出q后，根據(jù)關(guān)于R的后驗(yàn)概率P(R | X q; A)利用Metropolis-Hastings算法[23]將產(chǎn)生R的抽樣序列(k=K0, …, K; K0＞0為burn in周期)。基于這抽樣序列，可定義關(guān)于后驗(yàn)概率均值的二次型損失期望函數(shù)如式(19)所示。

最小化式(19)中的 J(R)可得到后驗(yàn)概率均值的估計值，它可由式(20)抽樣序列的經(jīng)驗(yàn)均值來逼近。

從上述的分析過程可知，由式(18)和式(20)結(jié)合可估計出r(t)，然后將其代入式(12)中就可以得到s(t)的估計值。

5 仿真實(shí)驗(yàn)分析

本文提出的源數(shù)未知的欠定盲源分離算法實(shí)現(xiàn)步驟如下。

step 1 源數(shù)n和混疊矩陣A的估計

1) 由式(3)和式(4)對 X進(jìn)行 S變換，得到XS(t,k)，t=0,…,N-1, k=0,…,N-1;

2) 由式(5)選擇XS(t,k)中具有單源性的信號，得到Ψs；

3) 利用文獻(xiàn)[20]中的FCM對Ψs進(jìn)行c(c=2,…,cmax)聚類，并由式(6)確定n；

4) 由式(7)估算混疊矩陣A；

step 2 源信號s的恢復(fù)

1) 利用最大似然法估計式(18)中的參數(shù)向量q;

2) 利用文獻(xiàn)[23]中的 Metropolis-Hastings算法和式(20)相結(jié)合估算r；

3) 將step 1中的A和r代入式(12)，即可得到s的估計值。

在仿真驗(yàn)證實(shí)驗(yàn)中選取的源信號時域波形如圖1所示，本文的仿真平臺為MATLAB R2010b。

為了更好地檢驗(yàn)基于S變換和聚類驗(yàn)證技術(shù)相結(jié)合的源數(shù)估計法，利用例1、例2和例3這3個例子加以驗(yàn)證。例1的源信號選{s1, s3}，例2的源信號選{s1, s2, s3}，例3的源信號選{s1,s2, s3, s4}，在這3個例子的仿真實(shí)驗(yàn)中A都是在[-1 1]之間隨機(jī)產(chǎn)生的，且 m=3。它們的仿真結(jié)果如圖2所示。由圖2可知，聚類數(shù)越接近于最優(yōu)值時它的聚類驗(yàn)證指數(shù)越小，由式(6)可以準(zhǔn)確地確定例 1、例 2和例 3的源數(shù)；由此可表明基于S變換和聚類驗(yàn)證技術(shù)相結(jié)合的估計法不僅能準(zhǔn)確估計欠定情況下的源數(shù)，也適用于超定情況。

圖1 源信號

圖2 聚類驗(yàn)證指數(shù)曲線

為了全面驗(yàn)證本文提出的UBSS方法的有效性，還將該方法與其他文獻(xiàn)的方法進(jìn)行比較加以驗(yàn)證。在這個仿真實(shí)驗(yàn)中，式(1)中的 A同樣也是在[-1 1]之間隨機(jī)產(chǎn)生的，它的值見式(21)，本文欠定混合矩陣估計方法與文獻(xiàn)[13]、文獻(xiàn)[16]、文獻(xiàn)[19]等不同方法估計的結(jié)果如表 2所示。利用本文、文獻(xiàn)[13]和文獻(xiàn)[24]等 3種不同UBSS方法實(shí)現(xiàn)的分離信號如圖3所示；源信號估計方法的性能由式(22)評價，利用它對源信號估計性能的比較如表3所示。

表2 不同混合矩陣估計方法的估計結(jié)果的比較

表3 不同UBSS方法分離結(jié)果的比較

由表2和表3的比較結(jié)果，以及圖3的波形可知，本文的欠定盲源分離方法不僅能較好地分離出超高斯和亞高斯2種不同分布的信號，同時它在A和源信號恢復(fù)的性能上也體現(xiàn)了它比其他方法更具有優(yōu)越性。另外，在仿真實(shí)驗(yàn)中，由式(18)得到源信號(s1,s2,s3,s4)分布函數(shù)中的 α的估計值分別為1.059 3，1.010 7，0.091 0，0.083 5；β的估計值分別 70.010 9，172.850 5，0.795 7，1.888 7，由 β 的值可知s1, s2為亞高斯分布信號，s3, s4為超高斯分布信號；為驗(yàn)證這一結(jié)論，直接計算源信號的峭度，(s1,s2,s3,s4)的峭度分別為1.500 1，1.500 0，6.417 0，5.830 0，當(dāng)峭度大于3的信號為超高斯分布，小于3的信號為亞高斯分布，由此說明上述結(jié)論是正確的。

圖4 不同源信號估計算法的恢復(fù)信號

為了進(jìn)一步評價第3節(jié)提出的源信號估計算法的性能，分別利用該算法與文獻(xiàn)[14,15]的算法分離出欠定情況下的EEG(腦電圖)信號進(jìn)行比較分析，3種源信號恢復(fù)算法的復(fù)雜度運(yùn)算量為：本文的算法需估計2n個參數(shù)，2n次最大化運(yùn)算，1次Metropolis-Hastings抽樣運(yùn)算；文獻(xiàn)[14]的算法需要估計2n個參數(shù)，6n次最大化運(yùn)算，6n次Metropolis-Hastings抽樣運(yùn)算；文獻(xiàn)[15]的算法需要估計2nL個參數(shù)，2nL次最大化運(yùn)算，L次Metropolis-Hastings抽樣運(yùn)算，L為該算法的收斂迭代步長。圖 4(a)中的源信號s1和s2為服從超高斯分布的EEG信號，而s3為服從亞高斯分布的噪聲信號，它們都取自文獻(xiàn)[25]；而混疊矩陣A的值見式(25)，它也是在[-1 1]之間隨機(jī)產(chǎn)生，m=2。在這個仿真實(shí)驗(yàn)中，假設(shè)混合矩陣A和源數(shù)n均為已知，圖4(b)為觀測信號，圖4(c)～圖 4(e)分別為 3種不同源信號估計算法的分離信號，它們對源信號估計的性能比較如表4所示。

表4 不同源信號估計算法實(shí)驗(yàn)結(jié)果的比較

由圖4和表4的比較結(jié)果可知，雖然文獻(xiàn)[15]中概率分布函數(shù)對信號建模具有較高的自由度，但大量參數(shù)的估計不僅使得該算法的計算復(fù)雜，且源信號恢復(fù)性能差；利用本文的算法恢復(fù)的源信號s1和s2的β估計值分別為1.116 8和1.559 1；與文獻(xiàn)[14,15]的算法比較，該算法具有計算簡單和估計性能優(yōu)越等特點(diǎn)。

6 結(jié)束語

針對欠定情況下的源數(shù)估計、混疊矩陣和源信號恢復(fù)等關(guān)鍵技術(shù)，本文提出了一種源數(shù)未知的欠定盲源分離算法。該方法首先利用S變換和聚類技術(shù)相結(jié)合的方法來估算源數(shù)和混疊矩陣；然后將源信號以零空間形式表示，這種表示形式將求解n個未知數(shù)的問題變換成求解(n-m)個的未知，再通過最大似然法估計關(guān)于它的后驗(yàn)概率以達(dá)到恢復(fù)源信號的目的。仿真結(jié)果表明了該方法不僅能較好地分離出服從不同分布的源信號，同時它比其他方法具有更好的估計性能。本文提出的源信號個數(shù)和混合矩陣的估計是利用聚類技術(shù)對 SSD的元素進(jìn)行聚類獲得的，而SSD元素的選取依據(jù)是實(shí)數(shù)的觀測信號在S變換時頻域上滿足提出的SSD條件；在源信號恢復(fù)方面，通過最大似然法估計關(guān)于它的后驗(yàn)概率只適用于恢復(fù)混合矩陣為實(shí)數(shù)矩陣情況下的源信號，但源信號的零空間表示法也適用于混合矩陣為復(fù)數(shù)矩陣情況，這為復(fù)數(shù)混合矩陣的欠定盲源分離方法的研究提供了重要的思路，同時復(fù)數(shù)混合矩陣的欠定盲源分離方法的研究將成為我們下一個研究目標(biāo)。

[1] LIU Y D, ZHOU Z T, HU D W. A novel method for spatio-temporal pattern analysis of brain fMRI data[J]. Science in China Series F, Information Sciences, 2005, 48(2): 151-160.

[2] ARAKI S, MAKINO S, BLIN A. Underdetermined blind separation for speech in real environment with sparseness and ICA[A]. Processings of ICASSP ’04[C]. Montreal, Canada, 2004. 881-884.

[3] 王榮杰, 周海峰, 詹宜巨. 船舶噪聲的自適應(yīng)分離技術(shù)[J]. 中國航海, 2011, 34(3): 10-15.WANG R J, ZHOU H F, ZHAN Y J. Adaptive separation technology of ship noise[J]. Navigation of China, 2011, 34(3): 10-15.

[4] OHNISHI N Y, IMIYA A. Independent component analysis of optical flow for robot navigation[J]. Neurocomputing, 2008, 171(10-12):2140-2163.

[5] ANNA T, LUIGI B, EMANUELE S. A markov model for blind image separation by a mean-field EM algorithm[J]. IEEE Transactions on Image Processing, 2006, 15(2): 473-482.

[6] BAI E W, LI Q Y, ZHANG Z Y. Blind source separation/channel equalization of nonlinear channels with binary inputs[J]. IEEE Transactions on Signal Processing, 2005, 53(7): 2315-2323.

[7] CICHOCKI A, AMARI S. Adaptive Blind Signal and Image Processing[M]. New York: Wiley, 2002.

[8] DOUGLAS S C, GUPTA M, SAWADA H. Spatio-temporal fastICA algorithms for the blind separation of convolutive mixtures[J]. IEEE Trans Audio, Speech, Lang Process, 2007, 15(5): 1540-1550.

[9] XIE S L, HE Z S, FU Y L. A note on stone’s conjecture of blind separation[J]. Neural Computation, 2005, 117(2): 321-330.

[10] EVEN J, MOISAN E. Blind source separation using order statistics[J].Signal Processing, 2005, 85(9): 1744-1758.

[11] FANG Y, ZHANG Y. A robust clustering algorithm for underdetermined blind separation of sparse sources[J]. Journal of Shanghai University(English Edition), 2008, 12(3): 228-234.

[12] BOFILL P, ZIBULEVSKY M. Underdetermined source separation using sparse representation[J]. Signal Process, 2001, 81(11):2353-2362.

[13] ABDELDJALIL A, NGUYEN L, KARIM A. Underdetermined blind separation of nondisjoint sources in the time-frequency do-main[J]. IEEE Transactions on Signal Processing, 2007, 55(3):897-907.

[14] CEMGIL A T, FEVOTTE C, GODSIIL S J. Variational and stochastic inference for bayesian source separation[J]. Digital Signal Processing,2007, 17(5): 891-913.

[15] SNOUSSI H C, IDIER J. Bayesian blind separation of generalized hyperbolic processes in noisy and underdeterminate mixtures[J]. IEEE Transactions on Signal Processing, 2006, 54(9): 3257-3269.

[16] WANG R J, ZHAN Y J. Application of similarity in fault diagnosis of power electronics circuits[J]. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, 2010, E93-A(6):1190-1195.

[17] STOCKWELL R G, MANSINA L, LOWER R P. Localization of the complex spectrum: the s transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001.

[18] KIM S Y, CHANG D Y. Underdetermined blind source separation based on subspace representation[J]. IEEE Transactions on Signal Processing, 2009, 57(7): 2604-2614.

[19] REJU V G, KOH S N, SOON I Y. An algorithm for mixing matrix estimation in instantaneous blind source separation[J]. Signal Processing, 2009, 89(9): 1762-1773.

[20] SUN H J, WNNG S R, JIANG Q S. FCM-based model selection algorithms for determining the number of cluster[J]. Pattern Recognition, 2004, 37(10): 2027-2037.

[21] CHEN R B, WU Y N. A null sapce method for over-complete blind source separateion[J]. Computational Statistics & Data Analysis, 2007,51(12): 5519-5536.

[22] 史習(xí)智. 盲信號處理[M]. 上海：上海交通大學(xué)出版社,2006.SHI X Z. Blind Signal Processing[M]. Shanghai: Shanghai Jiaotong University Press, 2006.

[23] ROBERT C, CASELLA G. Monte Carlo Statistical Methods[M].New York: Springer-verlag, 1999.

[24] KHOR L C. Robust adaptive blind signal estimation algorithm for underdetermined mixture[J]. IEE Proceedings-Circuits, Devices and Systems, 2006, 153(4): 320-331.

[25] CICHOCKI A, AMARI S, SIWEK K. ICALAB toolboxes [EB/OL].http://www.bsp.brain.riken.jp/ICALAB/ICALABSignalProc/benchmar ks, 2007.