1、2012 年 3 月 Journal on Communications March 2012第 33 卷第 3 期 通 信 学 报 Vol.33 No. 3基于零空间表示和最大似然的欠定盲源分离方法王荣杰 1,3,詹宜巨 2,周海峰 3(1. 中山大学 信息科学与技术学院,广东 广州 510006;2. 中山大学 工学院,广东 广州 510006;2. 集美大学 轮机工程学院,福建 厦门 361021)摘 要:针对欠定情况下源数的估计、解混叠矩阵和源信号恢复关键技术,提出一种源数未知的欠定盲源分离算法,首先利用 S 变换和聚类技术相结合来估算源数和混叠矩阵,然后将源信号以零空间形式表示,再通
2、过最大似然估计关于其后验概率以达到恢复源信号的目的。仿真实验结果表明了该方法不仅能同时分离服从超高斯分布和亚高斯分布的源信号,且它比其他传统的方法具有更优越的估计性能。关键词:欠定盲源分离;S 变换;模糊 c 均值聚类;零空间表示;最大似然中图分类号:TN911.7 文献标识码:A 文章编号:1000-436X(2011)03-Underdetermined blind source separation based onnull-space representation and maximum likelihoodWANG Rong-jie1,3, ZHAN Yi-ju2, ZHOU Hai
3、-feng3( 1. School of Information Science and Technology, Sun Yat-sen University, Guangzhou 510006, China;2. School of Engineering, Sun Yat-sen University, Guangzhou 510006, China;3. Marine Engineering Institute, Jimei University, Xiamen 361021, China. )Abstract: Aiming to the estimation of source nu
4、mbers, mixing matrix and source signals under underdetermined case, and a method of underdetermined blind source separation with an unknown number of sources was proposed. Firstly, an algorithm based on S transforms and clustering technology was introduced to estimate the number of sources and mixin
5、g mixtures. Then sources were represented as null space form, the recovery of source signals using a method based on maximum likelihood. The simulation results show that the proposed method can separate sources of super-Gaussian distrisutron and sub-Gaussian distribution, and compared to other conve
6、ntional algorithms, estimated mixing matrix and separated sources with higher accuracy.Key words: underdetermined blind source separation; S transforms; fuzzy c-mean clustering; null-space representation; maximum likelihood 1 引言近年来,鉴于盲源分离优越的假设前提,其模型被广泛用于数字通信、语音、图像处理、船舶工程、机器人导航和生物医学信号处理等领域16。所谓的盲源分离,
7、就是在源信号和混合系统(或传输通道) 等未知的情况下,仅根据源信号的统计特性,由观测到的混叠信号恢复出源信号。根据观测信号和源信号的个数,盲源分离可分为超(正)定盲源分离和欠定盲源分离。超正定盲源分离有很多成熟的方法,如 ICA 独立分量分析 7,8、修正的 Stone 法 9、阶统计代数法 10等。但这些方收稿日期:2010-12-20;修回日期:2011-12-17基金项目:国家自然科学基金资助项目(51179074)Foundational Items: The National Natural Science Foundation of China (51179074)2 通 信 学
8、报 第 33 卷法不适用于观测信号的个数少于源信号的个数的情况,欠定盲源分离(UBSS, underdetermined blind source separation)是本文的研究重点。传统的欠定盲源分离算法主要可分为 2 类:一类是利用源信号的稀疏性来处理,文献11提出利用聚类技术和最小化 l1 范数相结合的方法来恢复时域稀疏的源信号,文献12,13提出利用时频变换方法将非稀疏信号变换到时频域,以达稀疏的目的,它们认为在每个时频区域内只有一个源信号不为零;这些方法都选择服从超高斯分布的语音为源信号。另一类是利用广义分布模型作为源信号的概率密度的贝叶斯估计法 14,15,这类方法的主要缺点是
9、计算复杂,在一定程度上降低了源信号的 恢 复 质 量 。这 些 传 统 的 算 法 都 假 设 源 数 为 已 知 。 本 文 提 出 一 种源 数 未 知 的 欠 定 情 况 下 任 意 分 布 源 信 号 的 盲 分 离 方法 , 该 方 法 首 先 利 用 S 变 换 和 聚 类 技 术 相 结 合 来 估算 源 数 和 混 叠 矩 阵 , 然 后 将 源 信 号 以 零 空 间 形 式 表示,最后通过最大似然法(ML, maxima likelihood)估计关于它的后验概率以达到恢复源信号的目的。2 问题的描述假设 n 个彼此相互独立的未知源信号 s(t)=s1(t),s2(t),
10、,sn(t)T,通过一未知瞬时线性混合系统后,得到 m 个观测信号 x(t)=x1(t), x2(t), , xm(t)T。观测信号 x(t)与源信号 s(t)的关系可由式 (1)表示。x(t)=As(t) (1)其中,AR mn 为混合矩阵,它反映了混合系统或信道的传输特性,要求它为行满秩,t=0, ,N1 为时域采样点。为了便于分析,在推导过程中将 s(t)和 x(t)分别记为 s 和 x。在 m n 的情况下,给定 A,源信号 s 可由式(2)估计得到。(2)*s其中,A *为 A 的广义逆矩阵,A *= AT(AAT)。当 m0,那么式(1)可改写成式(8) 。x=( , 0)s (8
11、)由式(8)可知,s 的前 m 个值可由式(9)得到。i=1, ,m (9)/ii其余的(nm) 个 s 用一个(nm) 维列向量的 r 表示,r =r1, , rnmT。由此式(8)中 s 的解可用式(10) 来描述。(10)1()()mnnmx00srI式(10)中的 I 为 (nm)阶单位矩阵;为混合矩阵 A 的广义逆矩阵。1()nm0当混合矩阵 A 不满足式(8)中的特殊形式时,4 通 信 学 报 第 33 卷它的 s 解也就不能直接用式(10)来描述。但由于 A为行满秩,它的 SVD 奇异值分解中含 m 个不为零的特征值,可写成 A=U( , 0)V,将其代入式(1),并根据式(10
12、),可得到 s 的解如式 (11)所示。(11)1()T() mnnm0sVxrI其中, 为 A 的广义逆变换矩阵,记1()nA*;V 的后( nm)个列向量为 A 特征值为零对应的特征向量,它们组成的矩阵记为 V2,即为 A 的零空间。那么式(11)可简写成式(12)。s = A*x+V2r (12)由式(12)可知,当 A 已知,那么估计了(nm)个的 r 就相当于估计 n 个的 s。4.2 源信号的恢复 为了能同时分离服从任意分布的源信号,本文选取具有对称单峰分布特性的广义高斯模型 (GGM, generalized gaussian model) 18, 22作为源信号的概率密度,式(
13、13)为它的通用数学表达式。(13)(|,)exp2(1/)zpz其中, 为信号 z 的方差, 为 Gamma 函数;调节 的值可得出不同的分布函数模型,=2 表示标准高斯分布, 2 时为亚高斯分布。设 X= x(0), , x(N1),R= r(0), , r(N1),S= s(0), , s(N1),q= 1, , n, 1, , nT,并记S 的概率密度为 f(S|q)。那么可得 X 和 R 的联合分布函数如式(14) 所示。(14)(,|)(|)det)NpfASq为了估计每个源信号分布函数模型中的 和,构造关于参数向量 q 的先验似然函数如式(15)所示。(15)(|)(|,)(dA
14、ppXA由最大似然估计法可知(16)argmx(|,)(qAq即(17)(|,)(d0pAXqAq式(14)式(17)中,由于 A 独于立于 q,且分布p(X|A,q)和 p(A)是非负的,因此 q 的最大似然估计可简化为式(18) 所示。(18)argmx(|,)pqX式(18)中估计出 q 后,根据关于 R 的后验概率P(R | X q; A)利用 Metropolis-Hastings 算法 23将产生R 的抽样序列 (k=K0, , K; K00 为 burn in 周期) 。基于这抽样序列,可定义关于后验概率均值的二次型损失期望函数如式(19)所示。(19)2()| ; ()dJEp
15、RXqAR最小化式(19)中的 J(R)可得到后验概率均值的估计值 ,它可由式(20)抽样序列 的经验均值 k来逼近。(20)01| ;()KkEXqA从上述的分析过程可知,由式(18)和式(20)结合可估计出 r(t),然后将其代入式 (12)中就可以得到 s(t)的估计值。5 仿真实验分析本文提出的源数未知的欠定盲源分离算法实现步骤为:step 1 源数 n 和混叠矩阵 A 的估计1) 由式(3)和式(4)对 X 进行 S 变换,得到 XS(t,k),t=0,N-1, k=0,N-1;2) 由式 (5)选择 XS(t,k)中具有单源性的信号,得到s ;3) 利用文献20 中的 FCM 对
16、s 进行 c(c=2, ,cmax)聚类,并由式(6)确定 n;4) 由式(7)估算混叠矩阵 A;step 2 源信号 s 的恢复1) 利用最大似然法估计式(18)中的参数向量 q;2) 利用文献 23中的 Metropolis-Hastings 算法和式(20)相结合估算 r;3) 将 step 1 中的 A 和 r 代入式(12) ,即可得到 s 的估计值。在仿真验证实验中选取的源信号时域波形如第 3 期 王荣杰等:基于零空间表示和最大似然的欠定盲源分离方法 5图 1 所示,本文的仿真平台为 MATLAB R2010b。为了更好地检验基于 S 变换和聚类验证技术相结合的源数估计法,利用例
17、1、例 2 和例 3 这 3个例子加以验证。例 1 的源信号选s 1, s3,例 2 的源信号选s 1, s2, s3,例 3 的源信号选 s1, s2, s3, s4,在这 3 个例子的仿真实验中 A 都是在1 1之间随机产生的,且 m=3。它们的仿真结果如图 2 所示。由图 2 可知,聚类数越接近于最优值时它的聚类验证指数越小,由式(6)可以准确地确定例 1、例 2和例 3 的源数;由此可表明基于 S 变换和聚类验证技术相结合的估计法不仅能准确估计欠定情况下的源数,也适用于超定情况。图 1 源信号图 2 聚类验证指数曲线为 了 全 面 验 证 本 文 提 出 的 UBSS 方 法 的 有
18、效性 , 还 将 该 方 法 与 其 它 文 献 的 方 法 进 行 比 较 加以 验 证 。 在 这 个 仿 真 实 验 中 , 式 (1)中 的 A 同 样也 是 在 1 1之 间 随 机 产 生 的 , 它 的 值 见 式 (21),本 文 欠 定 混 合 矩 阵 估 计 方 法 与 文 献 13、 文 献16、 文 献 19等 不 同 方 法 估 计 的 结 果 如 表 2 所示 。 利 用 本 文 、 文 献 13和 文 献 24等 3 种 不 同UBSS 方 法 实 现 的 分 离 信 号 如 图 3 所 示 ; 源 信 号估 计 方 法 的 性 能 由 式 (22)评 价 , 利
19、 用 它 对 源 信 号估 计 性 能 的 比 较 如 表 3 所 示 。(21)0.649 .2- 0.657- .4887-1. . . .A表 2 不同混合矩阵估计方法的估计结果的比较估计值真值 本文提出的方法文献13 的方法文献16的方法文献19的方法A(1,1) 0.363 7 0.262 1 0.342 4 0.289 7A(1,2) 0.193 4 0.072 5 0.107 5 0.149 4A(1,3) 0.330 3 0.273 5 0.400 6 0.410 5A(1,4) 0.951 5 1.035 2 0.819 0 0.961 5A(2,1) 0.480 5 0.5
20、79 9 0.485 2 0.410 0A(2,2) 0.366 7 0.156 7 0.239 9 0.187 0A(2,3) 0.983 0 0.777 0 1.122 0 0.949 0A(2,4) 0.817 7 1.051 1 0.816 3 0.708 1A(3,1) 0.407 4 0.301 3 0.470 6 0.581 0A(3,2) 0.672 9 0.407 2 0.543 7 0.592 1A(3,3) 0.515 4 0.633 0 0.575 8 0.538 2A(3,4) 0.6523 0.6926 0.5080 0.7132,i=1, ,n (22)21log
21、iiiSSIR式中的 为与 对应的估计。isi6 通 信 学 报 第 33 卷图 3 不同 UBSS 方法的分离信号表 3 不同 UBSS 方法分离结果的比较UBSS 方法源信号 本文提出的方法文献24的方法文献13的方法s1 19.626 8 17.098 8 15.860 2s2 26.230 7 19.116 6 12.721 5s3 25.680 4 23.038 8 23.950 0s4 23.901 4 19.294 4 22.380 1由表 2 和表 3 的比较结果,以及图 3 的波形可知,本文的欠定盲源分离方法不仅能较好地分离出超高斯和亚高斯 2 种不同分布的信号,同时它在 A
22、 和源信号恢复的性能上也体现了它比其他方法更具有优越性。另外,在仿真实验中,由式(18)得到源信号 (s1,s2,s3,s4)分布函数中的 的估计值分别为 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 为 超 高 斯 分 布 信 号 ; 为 验 证 这 一 结 论 , 直 接 计算 源 信 号 的 峭 度 , (s1,s2,s3,s4)的 峭 度 分 别 为 1.500 1, 1.500 0, 6.417 0, 5.8
23、30 0, 当 峭 度 大 于 3 的 信 号为 超 高 斯 分 布 , 小 于 3 的 信 号 为 亚 高 斯 分 布 , 由 此说 明 上 述 结 论 是 正 确 的 。第 3 期 王荣杰等:基于零空间表示和最大似然的欠定盲源分离方法 7图 4 不同源信号估计算法的恢复信号为了进一步评价第 3 节提出的源信号估计算法的性能,分别利用该算法与文献14,15的算法分离出欠定情况下的 EEG(脑电图) 信号进行比较分析,3 种源信号恢复算法的复杂度运算量为:本文的算法需估计 2n 个参数,2n 次最大化运算,1 次Metropolis- Hastings 抽样运算;文献14的算法需要估计 2n
24、个参数,6n 次最大化运算,6n 次Metropolis-Hastings 抽样运算;文献15的算法需要估计 2nL 个参数,2nL 次最大化运算,L 次Metropolis-Hastings 抽样运算,L 为该算法的收敛迭代步长。图 4(a)中的源信号 s1 和 s2 为服从超高斯分布的 EEG 信号,而 s3 为服从亚高斯分布的噪声信号,它们都取自文献25;而混叠矩阵 A 的值见式(23),它也是在1 1之间随机产生,m =2。在这个仿真实验中,假设混合矩阵 A 和源数 n 均为已知,图 4(b)为观测信号,图 4(c)图 4(e)分别为 3种不同源信号估计算法的分离信号,它们对源信号估计
25、的性能比较如表 4 所示。(23) 0.357 .86 30.1 12576 A表 4 不同源信号估计算法实验结果的比较源信号恢复算法源信号 本文提出的算法文献 14的算法文献15的算法s1 17.397 7 15.544 2 3.510 1s2 14.363 6 11.555 5 7.504 4s3 26.232 1 25.383 6 17.474 1由图 4 和表 4 的比较结果可知,虽然文献15中概率分布函数对信号建模具有较高的自由度,但大量参数的估计不仅使得该算法的计算复杂,且源信号恢复性能差;利用本文的算法恢复的源信号 s1 和 s2 的 估计值分别为 1.116 8 和 1.559
26、 1;与文献14,15 的算法比较,该算法具有计算简单和估计性能优越等特点。6 结束语针对欠定情况下的源数估计、混叠矩阵和源信号恢复等关键技术,本文提出了一种源数未知的欠定盲源分离算法。该方法首先利用 S 变换和聚类技术相结合的方法来估算源数和混叠矩阵;然后将源信号以零空间形式表示,这种表示形式将求解 n 个未知数的问题变换成求解(nm) 个的未知,再通过最大似然法估计关于它的后验概率以达到恢复源信号的目的。仿真结果表明了该方法不仅能较好地分离出服从不同分布的源信号,同时它比其他方法具有更好的估计性能。本文提出的源信号个数和混合矩阵的估计是利用聚类技术对 SSD 的元素进行聚类获得的,而 SS
27、D 元素的选取依据是实数的观测信号在 S 变换时频域上满足提出的 SSD 条件;在源信号恢复方面,通过最大似然法估计关于它的后验概率只适用于恢复混合矩阵为实数矩阵情况下的源信号,但源信号的零空间表示法也适用于混合矩阵为复数矩阵情况,这为复数混合矩阵的欠定盲源分离方法的研究提供了重要的思路,同时复数混合矩阵的欠定盲源分离方法的研究将成为我们下一个研究目标。参考文献:1 LIU Y D, ZHOU Z T, HU D W. A novel method for spatio-temporal pattern analysis of brain fMRI dataJ. Science in Chin
28、a 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 ICAA. Processing of ICASSP 04C. Montreal, Canada, 2004. 881-884.3 王荣杰, 周海峰, 詹宜巨. 船舶噪声的自适应分离技术J. 中国航海, 2011, 34(3): 10-15.WANG R J, Z
29、HOU H F, ZHAN Y J. Adaptive separation technology of ship noiseJ. Navigation of China, 2011, 34(3): 10-15. 4 OHNISHI N Y, IMIYA A. Independent component analysis of optical flow for robot navigationJ. Neurocomputing, 2008, 171(10-12): 2140-2163.5 ANNA T, LUIGI B, EMANUELE S. A markov model for blind
30、 image separation by a mean-field EM algorithmJ. 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 inputsJ. IEEE Transaction on Signal Processing, 2005, 53(7): 2315-2323.7 CICHOCKI
31、A, AMARI S. Adaptive blind signal and image processingM. New York: Wiley, 2002.8 通 信 学 报 第 33 卷8 DOUGLAS S C, GUPTA M, SAWADA H. Spatio-temporal fastICA algorithms for the blind separation of convolutive mixturesJ. IEEE Trans Audio, Speech, Lang Process, 2007, 15(5): 1540-1550.9 XIE S L, HE Z S, FU
32、Y L. A note on stones conjecture of blind separationJ. Neural Computation, 2005, 117(2): 321-330.10 EVEN J, MOISAN E. Blind source separation using order statisticsJ. Signal Processing, 2005, 85(9): 1744-1758.11 FANG Y, ZHANG Y. A robust clustering algorithm for underdetermined blind separation of s
33、parse sourcesJ. Journal of Shanghai University(English Edition), 2008, 12(3): 228-234.12 BOFILL P, ZIBULEVSKY M. Underdetermined source separation using sparse representationJ. Signal Process, 2001, 81(11): 2353-2362.13 ABDELDJALIL A, NGUYEN L, KARIM A. Underdetermined blind separation of nondisjoin
34、t sources in the time-frequency domainJ. 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 separationJ. Digital Signal Processing, 2007, 17(5): 891-913.15 SNOUSSI H C, IDIER J. Bayesian blind s
35、eparation of generalized hyperbolic processes in noisy and underdeterminate mixturesJ. 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 circuitsJ. IEICE Transactions on Fundamentals of Electronics,
36、 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 transformJ. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001.18 KIM S Y, CHANG D Y. Underdetermined blind source separation based on subsp
37、ace representationJ. 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 separationJ. Signal Processing, 2009, 89(9): 1762-1773.20 SUN H J, WNNG S R, JIANG Q S. FCM-based model selectio
38、n algorithms for determining the number of clusterJ. Pattern Recognition, 2004, 37(10): 2027-2037.21 CHEN R B, WU Y N. A null sapce method for over-complete blind source separateionJ. Computational Statistics & Data Analysis, 2007, 51(12): 5519-5536.22 史习智. 盲信号处理M. 上海:上海交通大学出版社,2006.SHI X Z. Blind S
39、ignal ProcessingM. 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 mixtureJ. 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/benchmarks, 2007.作者简介:王荣杰(1981-),男,福建晋江人,博士研究生,主要研究方向为智能信息处理和电力电子电路故障诊断。詹宜巨(1955-),男,江苏南京人,中山大学教授、博士生导师,主要研究方向为电子信息技术应用。周海峰(1970-),男,福建莆田人,集美大学副教授,主要研究方向为智能信息处理。