孔轶艳

摘 要:盲源分离是在传输信道模型和信源均未知情况下,仅利用观测信号恢复出源信号的技术。超完备盲分离是指源信号数目大于观测信号数目情况下的盲源分离,由于混合过程中信息的丢失,使得超完备盲分离问题更加具有挑战性。详尽地阐述了超完备盲分离在瞬时混合情况下的数学模型、解决方案及研究现状,并将各算法归为两类,即两步估计法和交替估计法,同时总结了各算法的优缺点及其适用条件。

关键词:盲源分离;独立分量分析;超完备;两步估计法;交替估计法

中图分类号: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.