基于时间序列变分贝叶斯理论的信号盲源分离.doc

基于时间序列变分贝叶斯理论的信号盲源分离
摘　要：研究信号盲源分离中源信号和混合矩阵估计问题。独立分量分析盲源分离的不足之处在于不能估计混合矩阵和源信号的能量及顺序；变分独立因子分析盲源分离的不足之处在于依赖参数初值。将一般变分贝叶斯理论用于时间序列，推导出时间序列的变分贝叶斯期望极大算法。将此算法用于信号盲源分离，同时将传感器噪声逆方差的分布取为Wishart分布，得到了理论上更合理的后验分布参数更新规则。仿真数据和实际语音信号盲源分离结果表明这种方法可以比较准确地估计混合矩阵和源信号，在一定程度上弥补了独立分量分析和变分独立因子分析盲源分离的不足。
关键词：变分贝叶斯理论；盲源分离；时间序列；期望极大算法
Signal blind source separation based on time series variational Bayesians theory
Abstract：We study the issue of estimating the source signal and the mixing matrix in blind source separation of the signals. The deficiency of the blind source separation through independent component analysis is that it can not estimate the energies and orders of the mixing matrix and the source signals, whereas the deficiency of the blind source separation based on variational independent factor analysis is its strong dependence on initial parameters. Here, the general variational Bayesians theory is applied to the time series, and a novel variational Bayesian expectation-maximization algorithm for the time series is deduced accordingly. Moreover, the proposed method is applied to the blind source separation of the signals, and then more reasonable rules to update posterior parameters are deduced, where the distribution of the inverse variance of the sensor noises is assumed as Wishart distribution. The proposed algorithm was applied to the sources of simulation signals and real speech signals, the results show that the adopted method can accurately estimate the mixing matrix and the source signals, and it can reduce the deficiency of the independent component analysis and the independent variational factor analysis to a certain degree.
Key words：variational Bayesian theory; blind source separation; time series; expectation-maximization algorithm

基于时间序列变分贝叶斯理论的信号盲源分离
1　引　　言
由多个信号源和多个传感器构成的系统，传感器上观测到的信号是各个信号源经过线性或者非线性混合后的结果，混合可能是瞬时的，也可能是延迟的，而且传感器本身也有相关噪声。盲源分离就是利用源信号、传感器噪声以及混合过程的某些性质将各个源信号分离出来，同时计算导致混合的因素。对于瞬时线性混合的情况，目前广泛使用的方法是独立分量分析(ICA) [1-2]。当传感器上无噪声时，ICA效果非常好，计算速度极快。然而，一方面，随着传感器噪声方差变大，ICA效果变差；另一方面，ICA无法确定独立分量的能量、符号、顺序和混合矩阵[3]。而在实际问题中，往往需要从含有传感器噪声的观测数据中估计源信号和混合矩阵。为了解决传感器噪声问题，Attias提出独立因子分析(IFA)[4]。但是IFA对多信号源计算速度极慢。为了克服这一困难，Attias同时提出变分独立因子分析(VIFA) [4]。VIFA使用期望极大(EM)算法，估计出来的参数和源信号是局部最优，依赖于参数初值的选取。1999年，Attias将VIFA中“后验概率密度因子化”的思想发展成变分贝叶斯(VB)理论[5-6]。VB理论是以一定的先验分布概率在解空间搜索最优解，对参数的估计依赖于先验分布参数，先验分布参数越准确，估计出来的参数越准确。但在VB理论的实际应用中，文献[7]中的式(19)、文献[8]中的式(A.5,7-11)和文献[9]中的式(33~37)，都将后验参数更新规则表示为对时间序列求和。可以想象，如果采样的时间点数非常多，求和结果会非常大。此外，文献[10]中的式(4)、文献[11]中的式(1)和文献[8]中的式(25)，都假定传感器噪声不相关(方差为对角矩阵)，而传感器噪声一般是相关的。这些都使得VB理论的工程实现和应用不尽如人意。由此，本文考虑观测数据为时间序列时的VB理论，并将传感器噪声逆方差取为Wishart分布，估计源信号和混合矩阵。
2　时间序列变分贝叶斯期望极大算法
根据Attias的VB理论[6]，一个模型由3个部分构成，即观测数据Y、隐层数据H和模型参数Θ，计算H和Θ的后验分布参数是关键环节。变分贝叶斯期望极大(VBEM)算法分为两步：
1)在变分贝叶斯期望步(VBE-Step)计算隐层H的后验分布q(H)：
(1)
式中：< · >Θ表示在后验分布q(Θ)下计算平均。
2)在变分贝叶斯极大步(VBM-Step)计算模型参数Θ的后验分布q(Θ)：
(2)
式中：p(Θ)为模型参数的先验分布，平均< · >H在后验分布q(H)下进行。
当观测数据、隐层数据和模型参数为时间序列时，即Y={y1,y2,…,yT}，H={h1,h2,…,hT}，Θ={θ1,θ2, …,θT}，从式(1)可以推导出t时刻隐层的后验分布，
(3)
式中：平均是在后验分布q(θ t)下进行。若模型参数随时间变化(动态模型参数)，而只有一个或少数几个时间序列时，无法求出时刻t模型参数的后验分布。若模型参数不随时间变化(静态模型参数)，即θ=θ1=θ2=…=θT，通过一个时间序列就可以估计模型参数，此时，可以推导出时间序列VBEM算法的VBE-Step和VBM-Step,
(4)
(5)
值得注意的是，在式(5)中指数部分为对时间序列计算平均，而不是求和。
3　时间序列变分贝叶斯盲源分离算法
对于瞬时线性混合的情况，信号混合过程可以表示为：
, (t=1,2,…,T) (6)
式中：xt表示t时刻源信号向量(L维)，yt表示t时刻传感器上的观测信号向量(M维)，ut表示t时刻传感器上的0均值高斯噪声向量(M维)。A表示混合矩阵(M×L)，并设为不随时间变化(静态的)。在VB理论的实际应用中，必须确定观测数据和隐层的条件概率密度以及选取模型参数的先验和后验分布。在时间序列VBEM算法用于盲源分离时，根据式(6)，观测数据yt的条件概率密度函数为：
(7)
式中：Λ为传感器噪声逆方差。用高斯混合(MOG)模型表示源信号的概率密度函数[4]，隐层取为MOG模型的状态和源信号，即ht={st,xt}，则隐层的条件概率密度函数为：
(8)
式中：MOG模型参数为{ω, μ, ν}(ω为状态概率，μ为期望，ν为逆方差)。相应的隐层后验分布取为：

(9)
式中：{}为隐层{st,xt}的后验分布参数。模型参数的先验分布取为：
(10)
(11)
(12)

(13)
模型参数后验分布取为:
(14)
(15)
(16)
(17)
整个模型的参数为θ={A, Λ, ω, μ, ν}。将式(7)~(17)代入式(4)，进行时间序列VBE-Step计算，可推导出关于隐层后验分布参数的平均场方程：
(18)
(19)
(20)
(21)
式中：(在后验分布q(A)q(Λ)下平均)，其余几个平均在后验分布q(ω)q(μ,ν)下计算。平均场方程式(18)~（21）是非线性方程，可用迭代法求解[4]。注意式(18)，实际上不依赖于时间，这与Attias的VIFA[4]中的式(6.7)一致。将式(7)~（17）代入式(5)，进行时间序列VBM-Step计算，可推导出模型参数的后验分布参数更新规则：
(22)
(23)
(24)
 (25)
(26)
式中：<xt>和<xt xtT >是在隐层后验分布q(st,xt)下平均，而和表示在隐层条件后验分布q(xt|st)下计算平均。在推导式(22)时，使用了平均场近似。整个算法流程见图1。与文献[7-9]不同，式(22)~（26）是对时间序列平均，不会因采样时间点数过多而导致求和结果过大。
图1　时间序列变分贝叶斯盲源分离算法流程图
Fig.1 The flow chart of the time series VB blind source
separation algorithm
4　结果和讨论
为了验证前述算法，在Visual C++ 6环境下编写了实现算法的程序，分别对仿真数据和实际数据进行计算，并利用MATLAB将计算结果绘图。仿真实验真实源信号见图2(a)。根据式(1)，真实源信号混合再加传感器噪声，得到仿真的观测信号(见图2(b))。真实源信号高斯混合模型参数用文献[12]的方法计算。程序经过7次迭代，收验精度为6.694 73×10-5。估计出来的源信号经sym8小波软阀值3级小波包降噪[13](见图2(c))后，与真实源信号对应分量的误差分别为0.039 8、0.065 9、0.044 0、0.239 4。模型参数A、Λ、ω、μ、ν估计值与真实值的误差分别是 0.192 3、0.036 3、0.011 2、0.032 1、2.88×10-5。
(a)真实源信号
(a) True source signals
(b)混合加噪信号
(b) Signals mixed with sensor noises

(c)估计出来的源信号(经小波包去噪)
(c) The estimated source signals denoised with wavelet package
图2　仿真信号盲源分离实验
Fig.2 The blind source separation experiment
to simulation signals
语音数据来自4个真人说话录音，2个男声，2个女声，时间为2 s，单声道，16位，采样率22 050 Hz，如图3(a)所示。混合后加传感器噪声信号见图3(b)。估计出来的源信号见图3(c)。因为信号数据太长(44100个)，图中只绘制了前面1 000个。估计出来的源信号与真实源信号对应分量的误差分别为0.021 6、0.040 2、0.029 2、 0.030 4。模型参数A、Λ、ω、μ、ν估计值与真实值的误差分别是1.04×10-2、1.12×10-3、1.34×10-3、2.72×10-3、1.47×10-5。
(a)真实语音信号
(a) True speech signals

(b)混合加噪信号
(b)Signals mixed with sensor noises

(c)估计出来的源信号
(c)The estimated source signals
图3　语音信号盲源分离实验
Fig.3 The blind source separation experiment
to speech signals
从以上两实验可以看出，算法是有效的。本文进一步对算法的性能进行了全面的测试，对结果作如下讨论：
1)单次迭代时间复杂度=T×，与Attias的VIFA[4]相同，但是循环次数比VIFA少得多。
2)对高斯混合模型参数的先验分布参数和传感器噪声逆方差的先验分布参数的取值不敏感，但是对混合矩阵的先验分布参数的取值很敏感。
3)对混合矩阵估计有较大误差，出现这种偏差不是因为计算隐层后验分布参数时使用了平均场近似(见式(19))，而是估计混合矩阵后验参数时使用了平均场近似(见式(22))；另外，当传感器噪声相关性极强时，算法效果不太理想，也与这有关。
4)在传感器噪声较强时，比独立分量分析有明显优势。在计算速度方面，比不上独立分量分析。在弱噪声或噪声相关性不强时，独立分量分析仍是极其重要的盲源分离方法。
5)算法相当于以先验分布确定的概率在解空间搜索最优解，收敛速度也快，因此最优解也依赖于模型参数的先验分布。
6)算法对于观测信号个数少于源信号个数的欠定盲源分离问题，从理论上可以计算，但实际效果不理想，还须加上一些限制条件。文献[14]在此领域进行了研究。
5　结　　论
本文将Attias的变分贝叶斯理论应用于时间序列，得到了时间序列的变分贝叶斯期望极大算法，并实现信号盲源分离。与独立分量分析相比，可以用于传感器噪声较强的情况，并且能估计出源信号和混合矩阵，不足之处在于计算速度比不上独立分量分析。与Attias的变分独立因子分析相比，利用了先验知识，以先验分布确定的概率在解空间搜索最优解，最优解依赖于模型参数的先验分布，单次迭代计算速度与变分独立因子相同，但是收敛快，所以整体计算速度更快。与一般的贝叶斯理论一样，模型参数的先验分布参数的取值是计算性能最关键的环节，这取决于人们对研究问题先验知识的积累。
参考文献
[1] BELL A J, SEJNOWSKI T J. An information-maximization approach to blind separation and blind deconvolution[J]. Neural Computation, 1995,7:1129-1159.
[2] HYVARINEN A. Fast and robust fixed-point algorithms for independent component analysis[J]. IEEE Trans. on Neural Networks, 1999,10(3):626-634.
[3] HYVARINEN A, KARHUNEN J, OJA E. Independent component analysis[M]. John Wiley & Sons, Inc, 2001: 154.
[4] ATTIAS H. Independent factor analysis[J]. Neural Computation, 1999,11:803-851.
[5] ATTIAS H. Inferring parameters and structure of latent variable models by variational Bayes[C]. In Proc. 15th Conf. on Uncertainty in Artificial Intelligence, 1999.
[6] ATTIAS H. A variational Bayesian framework for graphical models[G]. NIPS, Cambridge, MA: MIT Press, 2000, 12:209-215.
[7] W.D. PENNY, J. KILNER, BLANKENBURG F. Robust Bayesian general linear models[J]. NeuroImage, 2007, 36:661-671.
[8] ROBERTS S., CHOUDREY R. Data decomposition using independent component analysis with prior constraints[J]. Pattern Recognition, 2003,36:1813-1825.
[9] NASIOS N., BORS G.. Variational learning for gaussian mixture models[J]. IEEE Trans. Syst. Man Cybern-B, 2006,36:849-862.
[10] KIEBEL J, DAUNIZEAU J, PHILLIPS C, et al. Variational Bayesian inversion of the equivalent current dipole model in EEG/MEG[J]. NeuroImage. 2008,39:728-741.
[11] ZUMER M, ATTIAS T, SEKIHARA K, et al. Probabilistic algorithms for MEG/EEG source reconstruction using temporal basis functions learned from data[J]. NeuroImage, 2008,41:924-940.
[12] DEMPSTER A, LAIRD N, RUBIN D. Maximum likelihood from incomplete data via the EM algorithm[J]. Royal Statistical Society Series B, 1977(39):1-38.
[13] DONOHO L. De-noising by soft-thresholding[J]. IEEE, Trans. on Inf. Theory, 1995, 41(3):613-627.
[14] 付宁, 彭喜元. K-Hough欠定盲信道估计算法[J]. 电子测量与仪器学报, 2008, 22(5):63-67.
FU N, PENG X Y. K-Hough underdetermined blind mixing model recovery algorithm[J]. Journal of Electronic Measurement and Instrument, 2008,22(5): 63-67.