超完備瞬時盲分離算法研究

孔軼艷

摘 要:盲源分離是在傳輸信道模型和信源均未知情況下,僅利用觀測信號恢復出源信號的技術。超完備盲分離是指源信號數目大于觀測信號數目情況下的盲源分離,由于混合過程中信息的丟失,使得超完備盲分離問題更加具有挑戰性。詳盡地闡述了超完備盲分離在瞬時混合情況下的數學模型、解決方案及研究現狀,并將各算法歸為兩類,即兩步估計法和交替估計法,同時總結了各算法的優缺點及其適用條件。

關鍵詞:盲源分離;獨立分量分析;超完備;兩步估計法;交替估計法

中圖分類號:TP911文獻標識碼:A

文章編號:1004-373X(2009)19-030-05

Research about Algorithms of Instantaneous Over-complete BSS

KONG Yiyan

(Liuzhou Vocational & Technical College,Liuzhou,545006,China)

Abstract:The Blind Source Separation (BSS) is a technique to separate the source only by the observed sinals,without any information of the channel and source.The over-complete blind source separation is the BSS when the number of the original signals is larger than that of the mixtures.The over-complete BSS is more challenging because some information are lost when mixing.In this paper,mathematics model,solution and current research of the instantaneous over-complete BSS are introduced in detail.Moreover,the algorithms are classified into two kinds,the two-step estimation approach and the alternately estimate approach.Finally,the advantage and disadvantage of these algorithms,also the application condition are concluded.

Keywords:blind source separation;independent component analysis;over-complete;two-step estination;alternate estimation

0 引 言

獨立分量分析[1](Independent Component Analysis,ICA)是一種通過最大化多維觀測向量的統計獨立性來尋找合適的線性變換的統計方法。其目標是在源信號向量和傳輸信道均未知的情況下,從觀測數據中重建和恢復各獨立源。獨立分量分析在盲源分離(Blind Source Separation,BSS)方面的應用,使得BSS的研究取得了突破性的進展,發展了很多優秀的算法。然而傳統的ICA是要求觀測信號數目大于或等于源信號數目的,現實中存在很多源信號數目大于觀測信號數目的情況,這種情況下的盲源分離被稱為超完備盲分離[1-13](over-complete BSS)。

由于觀測信號數目少于源信號數目,相當于信源在經過混合信道后,發生了有損壓縮,因此采用傳統的ICA通過對混合系統求偽逆的過程已無法恢復出源信號,這些丟失的信息無法僅通過數學手段來恢復,只能通過一些先驗、假設或限制條件(如:獨立性、稀疏性等)進行彌補。

本文就超完備瞬時盲分離的問題描述和目前的主要解決方案做了一個詳細的介紹,并依據不同算法的特點,將之歸納為兩類,總結了各自的優缺點及其適用條件。

1 問題描述

盲源分離的數據模型如下:

x=As+ε(1)

式中:s=[s1,s2,…,sN]T是N維未知獨立源信號矢量,經過線性系統A混合,再疊加噪聲信號ε=[ε1,ε2,…,εM]后,得到M維的觀測信號矢量x=[x1,x2,…,xM]T,線性系統A為M×N未知混合矩陣,在超完備情況下,假設M

完備情況下的盲分離一般可以同時估計出混合矩陣和源信號,但在超完備情況下,同時估計源和混合矩陣十分困難,所以分離過程一般分為兩步,先估計混合矩陣,再由估計出的混合矩陣和觀測信號去估計源信號。也可以在同一迭代過程中交替地估計源信號和混合矩陣。因此根據各算法的估計特點, 可以將各算法分為兩步估計法和交替估計法兩類,下面將分別介紹這兩類算法。

2 兩步估計法[7-10]

2002年,Theis F J在文獻[10]及其他學者研究成果的基礎上,提出了兩步實現超完備盲分離的思想[7],即BMMR和BSR,先由觀測信號估計混合矩陣,再在觀測信號和估計出的混合矩陣基礎上估計源信號。這種方法多基于無噪數據模型,即x=As,各變量定義同式(1)。

2.1 BMMR

混合矩陣的盲恢復可以采用聚類[10,12],幾何ICA[8,9]等方法來實現。聚類算法一般要求源信號比較稀疏,可以近似地認為在每一時刻僅有一個信號非零,而采用幾何ICA算法對稀疏性并不做特別要求。

2.1.1 聚類算法

采用無噪數據模型x=As。以M=2為例介紹該算法,系統模型可以寫為:

x1x2=a11a12…a1N

a21a21…a2Ns1s2髎N(2)

一般假設s1,s2,…,sN具有稀疏分布,那么在x1,x2中它們常常只是單獨出現。如果某一時刻t只有s1有值,則有xt1=a11st1,xt2=a21st1,因而此時xt1/xt2=a11/a21。可見,屬于這種情況的觀測數據x在x1~x2散點圖上將聚集在斜率為a11/a21的斜線上(見圖1)。同理,當只有si出現的時刻,x1~x2散點圖上將聚集在斜率為a1i/a2i的斜線上。那么在其散點圖上將會聚集N條線,根據這些線的斜率,再結合混合矩陣A各列的歸一化條件|ai|2=1(i=1,2,…,N),就可以把矩陣A的各列向量估計出來。

圖1 x1~x2散點圖(N=3)

推而廣之,當觀測信號數目M>2時,觀測數據的散點圖將聚集在各ai方向上,因此求各ai相當于找數據在M維空間的聚集位置(這可以通過神經網絡的C均值聚類/K均值聚類來實現)。這些ai即為混合矩陣A各列的矢量。

2.1.2 幾何ICA算法

該算法是采用了“贏家通吃”(winner-takes-all)規則,可以描述如下:在單位球SM-1糝M上選擇2N個單位根ω1,ω′1,…,ωN,ω′N,使得ωi=-ω′i(i=1,2,…,N),且每對根彼此線性獨立,這些根稱為神經元。取固定的學習率η:N→R,使得滿足η(n)>0,∑n∈Nη(n)=∞和∑n∈Nη(n)2<∞。然后按照下面的步驟進行迭代,直到遇到收斂條件。

選擇隨機變量x中的一個樣本x(t)∈R,若x(t)=0,則重新選擇一個非零的樣本(當x的聯合概率密度函數ρx為連續時將不會出現這種情況),將x(t)投影到單位球,得到:y(t)=x(t)/|x(t)|。假設ωi或ω′i是所有根中根據歐氏距離算出的離y(t)最近的點,則進行如下迭代:

ωi(t+1)=π{ωi(t)+η(t)sgn[y(t)-ωi(t)]}(3)

式中:π:RM\{0}→SM-1表示在RM中的M-1單位球SM-1上的投影。另外有:

ω′i(t+1)=-ωi(t+1)(4)

其他的神經元在本次迭代中則不進行更新。

所有的神經元均按照這種迭代方式進行更新。迭代的最終結果將獲得一組ω1,ω′1,…,ωN,ω′N,將ω1,ω2,…,ωN組成向量矩陣A,即為待恢復的混合矩陣。

2.1.3 算法比較

聚類算法要求源信號滿足一定的稀疏性,但隨著信號維數的增加,計算代價不會迅速增長。

幾何ICA算法從信號的幾何特性出發,便于理解,所以得到了更廣闊的應用。它對信號的稀疏特性沒有特別的要求。且Theis F J在文獻[8]中給出了矩陣等價的概念和算法收斂的條件,并嚴格證明了分離的混合矩陣的惟一性。幾何ICA的主要缺點就是隨著觀測信號維數的增長,樣本和收斂次數呈指數增長,因此幾何ICA一般都限制在低維情況。

2.2 BSR

在超完備基條件下,就算己知混合矩陣A,源信號的分離也是不惟一的,因此常常假設源信號滿足一定的稀疏性,然后利用稀疏分量分析理論,根據最大稀疏性準則來分離源信號,最短路徑算法和廣義FOCUSS算法都是以尋求信號的最大稀疏性表示為目標的源恢復方法,最短路徑算法是基于單次觀測樣本進行處理的,需要的樣本數據較多,但單次計算的代價較低,而FOCUSS則是對信號的所有樣本進行批處理的,效率較高,計算代價相對較高。

2.2.1 最短路徑算法[7,10,12]

最短路徑算法是以“L1范數最小”為目標的一種分解方法,此時要求:

xt=∑Nj=1ajstj, 且∑Nj=1|stj|=min(5)

上式前一部分說明,每個可行解都是各矢量ajstj的矢量和,而且其總和等于xt,上式的后一部分說明,在全部可行解中尋找路徑(矢量長度之和),最短的一條就是所求答案。

以M=2為例,說明最短路徑尋優方法,如圖2所示。

(1) 計算在某時刻t時xt的矢量方向θt=arctan(xt2/xt1)。

(2) 選擇各aj(j=1,2,…,N)中方向居于θt兩側,且與θt最接近的兩個矢量作為基矢量(圖中的aA與aB)。

(3) 以選定的兩個基矢量為方向構造平行四邊形,使其合成矢量為xt,從而確定矢量的長度stA和stB。也就是說此時的分解結果為xt=aAstA+aBstB,從而得到t時刻的源信號表示:

(st)i=[(aA|aB)-1xt]A,i=A

[(aA|aB)-1xt]B,i=B

0,otherwise(6)

從而把t時刻的觀測信號xt分解成只含有兩個信源的稀疏組合。

(4) 依次對t=1~T逐點進行,便將x分解成若干稀疏信源的組合。

圖2 最短路徑法的說明

aA和aB是A中最靠近xt的單位矢量,xt=aAstA+aBstB,而且|stA|+|stB|為最小。

2.2.2 似p范數代價函數和廣義FOCUSS算法[11,12]

分散性可以看作稀疏性的反義詞,很多文獻認為“似p范數”可以作為分散性的度量:

Jp(s)=∑Ni=1sip,0≤p≤1(7)

當p=0時,J0(s)將等于s中的非零元素的個數;p=1時,J1(s)將等于s中所有元素的絕對值之和。Jp(s)愈大說明能量越分散,所以最小化Jp(s)將會使s的能量越集中,稀疏性越大。因此為了在等式約束x=As下最小化式(7),定義Lagrange函數L(s,λ)。

L(s,λ)=Jp(s)+λ(x-As)(8)

式中:λ∈RN,為Lagrange乘子矢量,那么式(8)的平衡點可以估計如下:

齭L(s*,λ*)=齭Jp(s)-ATλ*=0

λL(s*,λ*)=x-As*=0(9)

式中:齭Jp(s)=pD-1(s)s,其中D(s)=diag[sip-2]是對角矩陣,解方程(9)可得:

λ*=p[AD(s*)AT]-1x(10)

s*=p-1D(s*)ATλ*=D(s*)AT[AD(s*)AT]-1x(11)

從而可以得到估計最優矢量s*的廣義FOCUSS迭代算法為:

s(k+1)=D[s(k)]AT{AD[s(k)]AT}-1x(12)

3 交替估計法

這類算法的特點是源信號和混合矩陣的估計過程交替進行,直到收斂。它通常以最大化似然函數為目標,通過梯度尋優的方式實現對混合矩陣的估計,再在某一假設源信號概率分布的情況下采用最大化后驗概率的方法估計源信號,常采用的源信號模型有拉普拉斯分布[3-5]、高斯混合分布[6]、廣義高斯分布等。

3.1 最大化似然函數法估計混合矩陣

對于式(1)定義的模型,可以得到觀測數據x的似然函數為:

L=p(x|A)=∫p(x|A,s)p(s)ds(13)

采用最大化對數似然函數的方法來估計混合矩陣,那么目標函數就變為:

=argmaxA[ln p(x|A)]

=argmaxA[ln∫p(x|A,s)p(s)ds](14)

式中:p(x|A,s)根據系統假設的噪聲模型的不同而不同,當噪聲滿足高斯分布時,有:

p(x|A,s)=12πM1σexp-‖x-As‖22σ2(15)

當噪聲假設為具有更寬泛分布的廣義高斯分布時,有:

p(x|A,s)=ρ2λΓ(1ρ)exp-|x-As|λρ(16)

式中:λ=Γ(1/ρ)Γ(3/ρ);Γ(?)表示伽瑪分布。

由于通信系統中很多情況下的噪聲都可以近似看作高斯噪聲,所以在后面的討論中,p(x|A,s)的分布一般采取式(15)。然而不管噪聲是何種分布,在過完備的情況下,由于A不可逆,使得式(13)中的積分計算相當復雜,只能對上述積分進行數值近似,常采用的近似方法有最大值近似 [2]和高斯積分法則近似[3-5]。最大值近似法一般要求被積函數p(x|A,s)p(s)具有相當尖銳的分布,近似的過程忽略了大量的概率信息,且為了保證算法的收斂,要對A的列矢量的長度加以限制,不然將會造成較大的近似誤差。然而高斯積分法則相對準確,對被積函數的分布沒有特別要求。下面將分別介紹這兩種近似方法。

3.1.1 最大值近似

1997年,Olshausen和Field提出了最大值近似方法[1],它假設式(13)的被積函數p(x|A,s)p(s)具有相當尖銳的分布,即:有一個最大值,其他的都很小,于是可以僅取此最大值作為式(13)的積分近似值,那么混合矩陣的估計式(14)就變為:

=argmaxA(17)

式中:arg(?)表示求值。記E(x,s|A)=-ln[p(x|A,s)p(s)],則由式(15)可知:

E(x,s|A)=12σ2‖x-As‖2-∑Ni=1ln pi(si)+Const.(18)

那么式(17)變為:

=argminA(19)

最小化E(x,s|A)的過程相當于最小化式(18)的前兩項,最小化第一項相當于減小重建向量的誤差,從而使張成信號子空間,最小化第二項使s趨于最大稀疏分布,這兩項在最小化過程中相互影響,從而使的估計達到最優值。

式(19)的優化過程一般分兩步:首先,固定A,最小化E(x,s|A),求得,然后再固定,用梯度下降法最小化式(18)可得估計矩陣。

3.1.2 高斯積分近似

由于高斯積分近似中有:

∫f(x)dx靋onst?f()?-dd2xlog f(x)|x=-12(20)

所以式(13)可以近似為:

∫p(x|A,s)p(s)ds=const?p(x|A,s)?

p()|H|-(1/2)(21)

式中:H是p(x|A,s)p(s)在s=的Hessian矩陣:

H=-ddssTln p(x|A,s)p(s)|s=(22)

將式(15)代入式(21)并對之兩邊取對數得:

ln∫p(x|A,s)p(s)ds

∝ln p()-12σ2|x-A|2-12ln|H|(23)

H=1σ2ATA-ddssTln p(s)|s=(24)

將式(23),式(24)代入式(14),然后通過梯度上升法尋找其最大值點,從而得到矩陣A自然梯度迭代估計式如下:

ΔA∝AAT氮礎ln P(x|A)-A[φ()T+I](25)

式中:φ()=[φ(1),φ(2),…,φ(N)]T,φ(i)=祃og P(i)礽稱為激活函數。

3.2 最大化后驗概率法估計源信號[3-6]

在估計出混合矩陣后,可以在假設源信號滿足一定的概率分布下,通過最大化信源的后驗概率來實現對源信號的估計:

=maxs P(s|x,A)=maxsP(x|A,s)P(s)(26)

采用對數表示后可得:

=argmaxsln p(x|A,s)+ln p(s)(27)

對式(27)的計算可以采用梯度上升法或擬牛頓法求極值,兩者相比,擬牛頓法計算復雜度大幅度增加,但卻加快了收斂速度。

定義目標函數J=ln p(x|A,s)+ln p(s),則梯度法的源信號迭代過程為:

s(k+1)=s(k)+ηs礘祍

=s(k)+ηs[-ATΨ()+φ()](28)

式中:=x-As,Ψ()=祃n P(1)1,祃n P(2)2,…,祃n P(M)礛T。

擬牛頓法的迭代過程[6]為:

s(k+1)=s(k)+ηs2J祍sT-1礘祍(29)

3.3 源信號模型

前面提到,超完備條件下的盲分離,由于混合過程中丟失的信息無法用傳統的數學手段來恢復,而只能通過一些先驗或假設來彌補這些丟失的信息,而在交替估計算法中,無論是對混合矩陣的估計,還是對源信號的估計,都需要知道源信號的概率密度函數,因此需要建立源信號模型,一般可以采用拉普拉斯分布或高斯混合分布。

3.3.1 拉普拉斯分布模型[3-5]

對于語音信號,大多具有超高斯分布,所以可以假設源信號服從拉普拉斯分布:P(si)∝exp(-α|si|),這時比較適合做語音信號的盲分離。通常還需假設源信號向量相互獨立,則:

P(s)=∏Ni=1P(si)∝∏Ni=1exp(-α|si|)(30)

此即為N維源信號的聯合概率分布。

3.3.2 高斯混合分布模型[6]

對于很多通信信號或圖像信號,其概率密度并不滿足拉普拉斯分布。于是提出了采用更寬泛的信號概率模型:高斯混合模型,它既可以描述單峰分布,又可以描述多峰分布;既可以描述亞高斯信號,又可以描述超高斯信號。即每一個信號都可以用一個Q維的高斯信號按照一定的比例混合得到:

p()=∑Qq=1κqpq(|λq)=

∑Qq=1κq(2π)N2|Σq|12exp-12(-μq)TΣ-1q(-μq)(31)

式中:N是s的維數;κq,μq和Σq分別是混合權重,均值和協方差矩陣,它們均可以通過最大期望算法(Expectation Maximization,EM)[6]或交替條件期望最大(ACEM)算法[13]進行估計,在此不做具體介紹,可參看相關文獻。

3.3.3 信號模型比較

拉普拉斯分布模型,參數較少,因此計算代價相對較低,但是一般只適用于超高斯稀疏分布的語音信號,一旦源信號與假設不同,將會造成較大的誤差。

高斯混合模型能更好地模擬不同的源信號分布,適用的范圍比較寬,但是參數很多,且會隨著源信號數目的增加呈指數增長,加上EM算法本身運算復雜度高,所以這種模型的運算量大,計算代價很高。

4 結 語

本文總結了大部分的超完備盲分離算法,將各算法歸納為兩步估計法和交替估計法兩類,這兩類算法各有許多不同的處理方法,文中較詳細地介紹了各算法的思想,且對各算法的適用條件和優缺點進行了對比和總結。

從文中介紹的各算法來看,目前超完備盲分離的研究基本上還都限于低維情況的研究,各算法的復雜度都較高,且都對系統模型加了許多假設和限制,所以今后超完備盲分離的研究重點是在非稀疏情況下的分離,在非對稱信號的情況下的盲分離,在高維觀測數據情況下的盲分離,以及計算復雜度的降低等方面的研究。

參考文獻

[1]Hyvarinen A,Karhunen J,Oja E.Independent Component Analysis[M].New York:John Wieley &Sons;,Inc.,2001.

[2]Olshausen B A,Field D J.Sparse Coding with an Overcomplete Basis Set:A Strategy Employed by VI[J].Vision Research,1997,37(23):3 311-3 325.

[3]Te-Won Lee,Lewicki M S.Mark Girolami,et al.Blind Source Separation of more Source than Mixtures using Overcomplete Representations[J].IEEE Signal Processing Letters,1999,6(4):87-90.

[4]Lewicki M S,Olshausen B A.Probabilistic Framework for the Adaptation and Comparison of Image Codes[J].Journal of the Optical Society of America: Optics,Image Science and Vision,1999,16(7):1 587-1 601.

[5]Lewicki M S,Sejnowski T J.Learning Overcomplete Representations[J].Neural Computation,2000,12(2):337-365.

[6]Khor L C,Woo W L,Dlay S S.Non-sparse Approach to Underdetermined Blind Signal Estimation[A].IEEE International Conference on Acoustic,System and Signal Processing[C].2005(5):309-312.

[7]Theis F J,Lang E W.Formalization of the Two-step App-roach to Overcomplete BSS[A].Proc.of SIP[C].Hawaii,2002.

[8]Theis F J,Jung A,Lang E W,et al.A Theoretic Model for Geometric Linear ICA[A].Proc.of ICA[C].2001:349-354.

[9]Theis F J,Lang E W.Geometric Overcomplete ICA[A].Proc.of ESANN[C].2002:217-223.

[10]Bofill P,Zibulevsky M.Underdetermined Blind Source Separation Using Sparse Representations[J].Signal Process,2001,81(11):2 353-2 362.

[11]Amari S,Cichocki A.自適應盲信號與圖像處理[M].北京:電子工業出版社,2005.

[12]楊福生,洪波.獨立分量分析的原理與應用[M].北京:清華大學出版社,2006.

[13]檀志斌.線性盲信號分離中多信號源少觀測源問題的研究[D].合肥:中國科技大學,2004.