一种强鲁棒性的盲分离混合矩阵估计方法

陈永强 王宏霞

①(西南交通大学信息科学与技术学院 成都 610031)

②(成都信息工程学院电子实验中心 成都 610225)

1 引言

稀疏成分分析(SCA)[1-7]目前已成为解决欠定盲分离的重要手段，大多数文献通常都是先估计出混合矩阵或假设混合矩阵已知的情况下然后再分离源信号，因此混合矩阵的精确估计对于源信号的分离很重要。文献[1-3]通过K-means和C-means算法来估计混合矩阵，但这两种方法对初始值要求都很高，而如何设置初始值通常是不知道的，不理想的初始值很容易使得算法落入局部极值点。文献[4]通过LOST算法来估计混合矩阵，但也很依赖于初始值，其性能并不比K-means好。文献[5]采用的分层聚类方法虽然不依赖初始值，但是需要设置类与类之间的最小门限值，该门限值非常难以设置，算法容易受到 “野值”的影响，造成错误估计。文献[6]采用系统聚类的方法来估计混合矩阵，但这种方法实质上还是分层聚类。

最近，文献[8,9]提出基于人工蜂群[10](Artificial Bee Colony,ABC)算法的聚类方法，并表明ABC算法相比传统的群智能算法，例如遗传算法[11](GA)、粒子群优化算法[12](PSO)具有更强的全局搜索能力。但是直接将文献[8,9]的聚类方法用于盲混合矩阵的估计，并不能达到很好效果。这是由于盲混合矩阵的估计具有其自身的一些特点，而且文献[8,9]只是简单地通过优化目标函数来达到聚类的目的，并没有充分利用数据聚类的特点，导致算法收敛速度缓慢，搜索精度不够高。本文根据稀疏信号混合矩阵估计的特点，充分利用观测数据聚类的特点，提出一种两阶段采用不同食物源编码方式和搜索策略的蜂群算法，通过新的蜜蜂搜寻行为，明显加快了收敛速度，并通过子蜂群之间的协作搜索方式加强局部搜索能力，进一步提高了混合矩阵的估计精度。本文提出的方法对初始值不敏感，在信噪比很低的情况下依然能够体现出很好的鲁棒性。

2 瞬时混合模型及单源点检测

假设源信号的个数为N，传感器的个数为M。对信号进行短时傅里叶变换(STFT)，在某个时频点(t,f)，瞬时混合模型如下所示。

其中X(t,f)=[x1(t,f),x2(t,f),…,xM(t,f)]T是M维观测数据向量；S(t,f)=[s1(t,f),s2(t,f),…,sN(t,f)]T是N维源信号数据向量；A是M×N维实混合矩阵，且是时不变的；N(t,f)=[n1(t,f),n2(t,f),…,nM(t,f)]T是M维高斯白噪声向量。

对于弱稀疏信号，通过单源点检测后能够呈现出更好的方向性，而且可以大大减小计算量。本文同时利用观测向量的实、虚部[5]以及相邻单源点之间的“相关性”[13]来检测单源点，将相邻频率的两个时频点观测向量的实部Xreal(t,f),Xreal(t,f′)和虚部Ximag(t,f),Ximag(t,f′)构成集合{Yii=1,2,3,4}，将满足

且长度较大的数据归一化并使第1个元素非负，然后构成单源点数据集合这里，代表向量长度，δ是一个正数(为了增加鲁棒性，选的可以大一些，本文仿真中取为0.4)。

3 基于人工蜂群算法的混合矩阵估计方法

文献[10]对ABC算法进行了详尽的描述。该算法模拟蜂群的采蜜行为，将蜂群中蜜蜂的种类分为雇佣蜂(employed bees)、观察蜂(onlooker bees)和侦察蜂(scouts bees)。雇佣蜂和观察蜂在整个蜂群中各占一半，每个食物源处只有一个雇佣蜂，雇佣蜂在特定时候可以转化成侦察蜂。为了提高蜜蜂对食物源的搜索效率，加强局部搜索能力，本文将算法分为两个阶段，第1阶段主要进行全局搜索，第2阶段进行局部搜索，两阶段采用不同的策略。下面分成全局搜索策略和局部搜索策略对两阶段算法进行描述。

3.1全局搜索策略

3.1.1 食物源编码方式及适应度函数本文在第1阶段采用矩阵编码方式，即每个食物源的位置用一个矩阵表示， 即

每个食物源的位置其实就是混合矩阵的一个潜在解。都是归一化的列向量，即

在聚类问题中通常采用欧氏距离作为数据和聚类中心之间“相似度”的度量[8,9]。欧氏度量对点聚类问题更加适合，而混合矩阵的估计是一个线聚类问题。虽然也可以考虑将直线投影到球面上，变成球面上的点，然后用欧氏距离来度量点与点之间的相似性，但是这样做会带来问题。假设某个时频点处只有s1起作用，观测向量X=a1s1(a1是A的第1列)。由于s1可正可负，因此为了衡量距离，需要将X投影到上半部分单位球面上。例如a1=[0.0159- 0 .6828 - 0 .7304]T,X在上半部分单位球面的投影只可能是a1。但实际情况下，源信号并不是完全稀疏的，除了s1的主导作用外，其它源也会起次要作用，导致X投影后的单位向量可能变为X=[0.0170 0.6981 0.7158]T，显然这两个数据的方向很接近，但欧氏距离却很大。为此，本文采用单位向量的内积绝对值来度量相似度，并以此来构建新的目标函数(适应度函数)，该目标函数如式(4)所示。

这里，max{·}表示取集合中的最大值。由于向量都做了归一化，因此两个向量方向越相同，内积绝对值就越大。

3.1.2 通过改变雇佣蜂行为提高收敛速度在基本ABC算法中，每个蜜蜂的搜索行为具有很大的随机性，从而使得算法的收敛速度缓慢。本文在第1阶段提出一种提高蜂群搜索效率的策略，其基本思想是：单源点观测数据会在一些方向上较为集中，如果将单源点数据归属到iU的各个列构成各个集合，那么这些集合的中心附近存在真正聚类中心的可能性一定大于随机产生的候选解。因此，在雇佣蜂寻找食物的阶段，本文将其搜索行为改为如下方式：

首先根据第i个食物源位置对单源点Xssp(l)(l=1,2,…,L)数据进行归类。如果Xssp(l)和的内积绝对值最大，那么候选食物源的位置由下式产生

3.1.3 观察蜂行为保证全局搜索能力观察蜂选择第i个食物源的概率仍然按照文献[10]计算，即

由于基本 ABC算法是在向量编码下产生候选解的，因此这里需要将其拓展到矩阵编码的情况，本文的观察蜂采用式(7)来产生候选食物源。

j,k分别是从集合{1,2,…,N},{1,2,…,I}中按等概率方式随机选取的，且k≠i。这里的,,都是向量，也就是说，只随机改变Ui的某一列，其它列均不变，从而得到候选解。由于只改变矩阵的某一列，仍然是在食物源位置的邻域内搜索。若gf()＜gf(),则取代，否则不变。这里也要进行归一化处理。

3.1.4 算法流程第1阶段算法的完整流程如下：

(1)首先确定参数：食物源个数I和观察蜂循环搜索次数limit；

(2)随机产生I+1个M×N维矩阵(列归一化)，分别作为I个食物源的初始位置和全局最优解G=[g1g2…gN]的初始位置；

(3)雇佣蜂按式(5)产生候选解，根据适应度函数(式(4))决定候选解取舍；

(4)观察蜂先按式(6)对食物源进行选择，然后按式(7)产生候选解，并根据适应度函数决定其取舍，这一过程循环进行limit次；

(5)在所有蜜蜂结束搜索后，将最好食物源Ubest(即gf(Ui)i=1,2,…,I中最大的那个食物源)的适应度函数值gf(Ubest)和gf(G)进行比较，如果gf(Ubest)＞gf(G)，则G=Ubest，否则G不变；

(6)若某个食物源邻域被观察蜂搜索limit次，仍然无法找到更好的食物源，则该食物源处的雇佣蜂转化成侦察蜂，在整个解空间随机产生一个新的食物源；

(7)如果gf(G)已经几乎不变，第1阶段搜索结束，否则回到(3)，再开始下次循环。

3.2局部搜索策略

在第1阶段算法结束后，各个聚类中心的估计值gn(n=1,2,…,N)并不一定非常准确。为此，本文采用子蜂群之间的协作，通过子适应度函数的优化来更新聚类中心的位置，然后通过适应度函数值的好坏判断，对更新后的聚类中心进行取舍，保证估计值向真正的聚类中心一步步靠近。

3.2.1 编码方式和子适应度函数将蜂群分解为N个子蜂群，每个子蜂群只负责对特定的若干食物源进行搜索，这些食物源在当前的聚类中心附近产生。这些由第n个子蜂群负责搜索的食物源位置用Wn,i(i=1,2,…,Qn)表示，Qn是由第n个子蜂群负责搜索的食物源个数，食物源位置采用向量编码方式，即

Wn,i是混合矩阵的部分潜在解，即混合矩阵第n列(各个聚类中心)的潜在解。

每个子蜂群都有各自的适应度函数，本文称为子适应度函数，其定义为

这里，Cgn表示属于当前聚类中心gn的数据集合。显然，如果Wn,i使得式(9)达到最大，Wn,i就是该集合的聚类中心。由于这里采用向量编码，各个子蜂群就可以根据集合Cgn中的数据，直接用基本ABC算法对式(9)进行优化。不过这里删除侦察蜂，因为侦察蜂的作用主要是为了加强ABC算法的全局搜索能力，提高食物源的多样性，而第2阶段的目的是进行局部搜索，因此侦察蜂失去了作用。

3.2.2 加强搜索效率和加快收敛速度的初始化方法第n个子蜂群负责搜索的Qn个食物源的初始化按下式进行

这里，Gj,n表示G的第j行第n列所对应的元素；k是从{1,2,…,M}随机取的数；δ是一个较小的常数；randn是服从正态分布N(0,1)的一个随机数。同样，初始化后要对Wn,i进行归一化。显然，食物源都在聚类中心的附近产生，这样可以加快算法收敛速度。

3.2.3 子蜂群协作算法流程第2阶段其实就是通过各个子蜂群的协同优化作业来达到获得更好聚类中心的目的，下面给出这一算法的完整流程：

(1)首先将蜂群分解为N个子蜂群，对于每个子蜂群、Qn、循环次数MCN都取为一样；

(2)对单源点数据Xssp(l)进行分类，若|Xssp(l)T⋅gj|(j=1,2,…,N)最大，则Xssp(l)∈Cgj；

(3)对食物源按式(10)进行初始化；

(4)每个子蜂群按照基本 ABC算法循环MCN次。每个子蜂群也有各自的最优部分解每次循环根据子适应度函数值(式(9))决定是否要对其进行更新；

(5)当各子蜂群结束搜索后，按式(4)计算P=[p1p2…pN]的适应度函数值gf(P)，如果gf(P)＞gf(G)，则G=P，否则G保持不变；

(6)若gf(G)几乎不变时，则第2阶段搜索结束，最后得到的G就是对混合矩阵A的估计，即=G；否则回到(2)，进行下一次循环。

4 计算机仿真

本文将语音信号作为源信号，语音采样率为16 kHz，持续时间3 s。本文进行混合矩阵估计时，都在时频域里进行，短时傅里叶变换采用窗长为1024点的汉宁窗，重叠率为50%。以下分成几个方面将本文方法与已有方法进行比较。

针对不同的源信号个数，本文的改进ABC算法参数都设置如下：第1阶段的食物源个数I=10,limit=20，当gf(G)的变化量小于0.01就结束；第2阶段每个子蜂群的食物源个数Qn=5,MCN=20,当gf(G)的变化量小于0.001时算法结束。

(1)混合矩阵的估计 这里将本文提出的混合矩阵估计方法和文献[2]，文献[5]，文献[14]采用的混合矩阵估计方法进行比较。文献[2]采用K-means算法、文献[5]采用的是分层聚类(hierarchical)方法、文献[14]则采用PSO算法。在估计混合矩阵之前，各算法都先进行单源点检测。为了公平比较，这里都采用本文的单源点检测方法，并且这里将文献[14]采用的基本PSO算法替换为最近由文献[12]提出的一种改进PSO算法(称为pPSA)，并以本文的式(4)作为目标函数。仿真发现，当pPSA算法的种群规模为50，权因子w=0.2，学习因子c1=c2=2可以取得相对较好的效果，其它参数和文献[12]一样。pPSA一般要到150步以后才能收敛，因此这里迭代次数(循环次数)设置为180次。K-means算法、pPSA算法、本文方法都用随机的方式产生初始值。

混合矩阵均方估计误差MSE(dB)采用下式来计算

图1是不同算法在3个传感器8个源信号时得到的MSE和信噪比(MSE-SNR)关系曲线，这里的估计误差MSE是随机产生50次混合矩阵得到的平均误差(由于pPSA耗时太长，只取20次进行平均)，混合矩阵的各列向量之间的夹角大于15°。可以看到，由于K-means算法对初始值很敏感，经常会收敛到错误位置，造成很大误差，MSE是最高的。pPSA有时同样会落入局部极值点，虽然比K-means稍好，但仍然有较大的MSE。分层聚类算法由于不受初始值影响，因此MSE要小一些。图中，K-means(repeated)表示K-means采用不同初始值重复运行40次，然后取最好的那一次(使式(4)目标函数值最大)作为混合矩阵的估计。可以看到，通过重复运算，K-means的性能得到了很大改善。但是发现，即使在SNR较高时，对50次不同混合矩阵估计中，K-means仍然有2～3次落入了局部极值点，造成较大误差。从图中可以看到，本文提出的方法明显好于其它方法，本文算法虽然也需要设置初始值，但对初始值鲁棒性非常好，因此不需要设置不同的初始值重复运行也可以获得很高的精度。

图1 MSE-SNR关系曲线图

表1是3个传感器8个源信号时，SNR=35 dB时各算法的鲁棒性比较。VAR是随机产生50次混合矩阵得到的估计误差MSE(dB)的方差，它可以衡量算法的鲁棒性。分层聚类和K-means都不是直接优化目标函数，容易造成错误聚类。虽然PSO算法直接优化目标函数，但其全局优化能力有限，容易落入局部极值点。这是因为PSO在进化后期，各个粒子趋于相同，失去多样性，造成“早熟”，这也是目前其它一些群智能算法的缺点。ABC算法在进化后期，虽然各食物源也会收敛到最优解附近，但是食物源位置(i=1,2,…,I)各列的次序是不相同的，因此并不会失去多样性，通过式(7)仍然可以产生新的候选解，保证了其具有很强的全局搜索能力，因此本文方法的鲁棒性比其它方法要好很多。当信噪比降低时，其它算法会产生更多次的较大误差。而本文方法即使在SNR=5 dB，单源点数据呈现出的方向性已经变得很差时，仍然没有产生过一次较大误差。

表1 算法鲁棒性比较

(2)ABC改进算法的优势 图2是3个传感器8个源信号时的性能比较。图中，ABC表示用文献[8]聚类方法估计混合矩阵得到的结果，不过这里增加了归一化处理，并且采用式(4)作为目标函数，食物源个数和limit参数与本文相同；ABC*算法表示采用全局搜索策略，并按本文的方法改变蜜蜂搜索行为得到的结果；ABC**表示采用两阶段不同搜索策略，而且改变蜜蜂搜索行为得到的结果。

图2(a)是在收敛速度方面的比较。由于ABC只有一个循环，而本文有两个循环(内循环和外循环)，为了公平比较，这里的循环次数是指总的循环次数(内循环次数和外循环次数的乘积)。适应度函数值是30次不同混合矩阵下得到的平均值。

可以看到，ABC和ABC*的初始值是完全一样的。ABC*在收敛过程中出现两次跳跃，第1次出现在循环刚开始的时候，这是因为随机产生的iU的某些列和聚类中心比较靠近，通过式(5)产生候选解Vi,Vi的某些列就会迅速接近聚类中心，从而使适应度函数值迅速增加。随着各列在循环过程中都比较接近聚类中心时，Vi的各列都会逼近聚类中心，于是适应度函数出现了第2次跳跃。ABC*算法平均88步就可以收敛，而ABC平均238步才收敛。因此，改变蜜蜂的搜索行为明显加快了收敛速度。

图2(b)是在混合矩阵估计误差方面的比较。可以看到，通过改变蜜蜂搜索行为不但可以加快收敛速度，还可以提高搜索精度。由于第2阶段加强了算法的局部搜索能力，所以采用两阶段搜索策略，可以进一步提高算法的搜索精度，减小估计误差。ABC*比ABC算法的适应度函数值提高了0.12,ABC**则提高了0.23，但估计精度却可以得到明显提高，充分说明了精确搜索的重要性。由于计算子适应度函数(式(9))时，不需要再对Xssp(l)进行归类，因此相比适应度函数(式(4))的计算量要小得多，而且只需要60步(总循环次数)左右就可以收敛，因此第2阶段获得更高估计精度的代价是很小的，这一点可以从表2看出。

表2 各算法运行时间(s)

(3)各算法运行时间比较表2是3个传感器8个源信号时各算法运算时间的比较。可以看到，分层聚类用时最短，pPSA用时最长。改进的ABC算法虽然不是用时最短的，但是运算时间是可以接受的。ABC*的运行时间只是ABC的39%。ABC**只比ABC*多出了4%的计算量，运行时间也只是ABC的41%。

图2 改进ABC算法与基本ABC算法的比较

5 结束语

本文根据稀疏信号的特点提出两阶段不同搜索策略的蜂群算法，并通过改变蜜蜂搜索食物的方式，明显提高了搜索精度，加快了收敛速度，减小了算法运行时间，将其用于盲分离混合矩阵的估计取得了很好的效果。目前已有的混合矩阵估计方法只在源信号很稀疏或者源信号个数较少时才能有相对较好的性能，但对于语音等弱稀疏信号，尤其当源个数比较多时，混合矩阵的估计经常会出现较大误差。本文的混合矩阵估计方法解决了这一问题，体现出了非常好的鲁棒性。而且，本文方法在具有很好性能的同时并不需要太大的计算量，这一点也是非常重要的。

[1]Aissa-El-Bey A,Linh-Trung N,Abed-Meraimet K,et al..Underdetermined blind separation of nondisjoint sources in the time- frequency domain[J].IEEE Transactions on Signal Processing,2007,55(3)：897-907.

[2]Kim S G and Yoo C D.Underdetermined blind source separation based on subspace representation[J].IEEE Transactions on Signal Processing,2009,57(7)：2604-2614.

[3]陈晓军,成昊,唐斌.基于ICA的雷达信号欠定盲分离算法[J].电子与信息学报,2010,32(4)：919-924.Chen X J,Cheng H,and Tang B.Underdetermined blind radar signal separation based on ICA[J].Journal of Electronics＆Information Technology,2010,32(4)：919-924.

[4]O’Grady P D and Pearlmutter B A.The LOST algorithm：finding lines and separating speech mixtures[J].EURASIP Journal on Advances in Signal Processing,2008,784296：1-17.

[5]Reju V G,Koh S N,and Soon I Y.An algorithm for mixing matrix estimation in instantaneous blind source separation[J].Signal Processing,2009,89(9)：1762-1773.

[6]陆凤波,黄知涛,姜文利.基于时频域单源区域的延迟欠定混合非平稳信号盲分离[J].电子学报,2011,39(4)：854-858.Lu F B,Huang Z T,and Jiang W L.Underdetermined blind separation of time-delayed non-stationary signal based on single source region in the time-frequency domain[J].Acta Electronica Sinica,2011,39(4)：854-858.

[7]Lu F B,Huang Z T,and Jiang W L.Underdetermined blind separation of non-disjoint signals in time-frequency domain based on matrix diagonalization[J].Signal Processing,2011,91(7)：1568-1577.

[8]Karaboga D and Ozturk C.A novel clustering approach：Artificial Bee Colony (ABC)algorithm[J].Applied Soft Computing,2011,11(1)：652-657.

[9]Zhang C S,Ouyang D T,and Ning J.An artificial bee colony approach for clustering[J].Expert Systems with Applications,2010,37(7)：4761-4767.

[10]Karaboga D and Basturk B.A powerful and efficient algorithm for numerical function optimization：Artificial Bee Colony (ABC)algorithm[J].Journal of Global Optimization,2007,39(3)：459-471.

[11]李敏强,寇纪凇,等.遗传算法的基本理论与应用[M].北京：科学出版社,2002：17-75.Li M Q,Kou J S,et al..The Basic Theory of GA and the Application[M].Beijing：Science Press,2002：17-75.

[12]Zhao X C.A perturbed particle swarm algorithm for numerical optimization[J].Applied Soft Computing,2010,10(1)：119-124.

[13]Bofill P.Identifying single source data for mixing matrix estimation in instantaneous blind source separation[C].International Conference on Artificial Neural Networks(ICANN),Prague,2008：759-767.

[14]高鹰,谢胜利,许若宁,等.基于粒子群优化算法的稀疏信号盲分离[J].系统仿真学报,2006,18(8)：2264-2266.Gao Y,Xie S L,Xu R N,et al..Blind sparse source separation based on particle swarm optimization[J].Journal of System Simulation,2006,18(8)：2264-2266.