盲源分离的SVM概率密度函数估计算法

2009，45（17）Computer Engineering and Applications计算机工程与应用
1引言
随着信号处理技术和神经网络技术的日趋成熟，盲源分离技术已经被广泛地应用在通信、语音识别、图像处理等众多领
域。

盲源分离的信号模型为Z=A×S+Noise，其中Z=[z1，z2，…，z N]T
是N维观测矢量，S=[s1，…，s M]T是M维独立分量矢量，A为N×M的混合矩阵，Noise为噪声。

盲源分离就是在混合矩阵A和信号源S未知的情况下，寻找分离矩阵W，使得Y=WZ为信号源S的一个最优估计。

Cardoso对盲分离领域的发展进行了详细地介绍并指出信源自适应是一个比较困难的问题[1]。

Boscolo 提出用投影搜索的核密度方法估计代价函数[2]。

Vlassis用一种加权高斯和的方法估计信源概率密度[3]，该方法能实现对超高斯和亚高斯信源的估计，但误差较大并且具有较高的时间复杂度。

Choi提出用扩展高斯模型估计超高斯与亚高斯分布[4]，实现盲分离，该方法的局限在于只能估计单峰，判别超高斯与亚高斯用经验选择系数不够准确。

盲信号抽取定点（Fixed-Point）算法由Hyvarinen和Oja提出[5]，也称为快速ICA（Fast ICA），它通过多层神经网络顺序从混合信号中抽取源信号。

该算法分别对分离矩阵的每一列进行更新，每一次找到一个独立分量，这种算法可以运用于任何的非平方的代价函数来获得对独立性的估计。

设第i个抽取单元的输出为：
y
i
=w T
i
z（1）式中：w i=[w i1，w i2，…，w in]T，z是去均值和白化处理后的混合信号。

对于任一抽取单元输出y i，取负熵作为其非高斯性度量的代价函数：
Δw
i
∝
坠J（y
i
）
坠w
i
=γE[zf（w T
i
z）]（2）
式中w T
i
是分离矩阵W的第i行，γ=E[F（y
i
）]-E[F（v）]，v为
N（0，1），F（·）是用来逼近概率密度的非线性函数，f（·）是F（·）的导数。

J（·）是负熵，它是任意概率密度函数和具有同样方差的高斯型概率密度函数间的K-L散度。

固定点算法的表达式为[5]：
w
i
（k+1）=E{zf[w T
i
（k）z]}-E[f′（w T
i
（k）z）]w
i
（k）（3a）
w
i
（k+1）←
w
i
（k+1）
‖w
i
（k+1）‖
2
（3b ‖
‖
‖
‖‖
‖
‖
‖
‖‖
‖
）
收敛后的w T
i
是矩阵W中的一行，因此由y
i
（n）=w T
i
z（n）分
离出某一信源s
i
（n）。

盲分离算法中，通常采用非线性函数代替信号的激活函数；而对于统计特性未知的信号源，非线性函数只能凭经验选取。

如
盲源分离的SVM概率密度函数估计算法
胡波平，何选森
HU Bo-ping，HE Xuan-sen
湖南大学计算机与通信学院，长沙410082
School of Computer and Communication，Hunan University，Changsha410082，China
E-mail：hbp1982@
HU Bo-ping，HE Xuan-sen.Novel BSS algorithm for estimating PDF based on puter Engineering and Applications，2009，45（17）：142-144.
Abstract：This paper utilizes radial basis function to construct the kernel function and forms a sparse representation of the probability density of the mixed signal based on SVM regression method of Neural network theory，so the expressions of the PDF of the output signals is gotten and a new method of estimating the active function is put forward.A novel BSS algorithm is pre－sented in this paper by combining the PDF estimation with the fixed-point algorithm.The simulations have verified that this algo－rithm can successfully separate the mixed sub-Gaussian and super-Gaussian source signals.
Key words：support vector machine；probability density function；active function；fixed-point algorithm
摘要：基于神经网络理论中的支持向量机回归方法，利用径向基函数构造核函数，给出类高斯函数的混合信号概率密度稀疏表达，进而得到输出信号的概率密度的显式表达；提出一种估计激活函数的新方法，与盲信号抽取定点算法相结合，形成一种新的盲分离算法。

通过仿真实验，验证了该方法能成功地分离超、亚高斯混合信号。

关键词：支持向量机；概率密度；激活函数；定点抽取算法
DOI：10.3778/j.issn.1002-8331.2009.17.043文章编号：1002－8331（2009）17-0142-03文献标识码：A中图分类号：TN911
基金项目：湖南省湘潭市科技计划项目（No.ZJ20071008）。

作者简介：胡波平（1982－），男，硕士研究生，主要研究领域为盲信号处理、阵列信号处理；何选森（1958-），男，副教授，研究方向为信号与信息处理。

收稿日期：2008-04-09修回日期：2008-06-23
142
2009，45（17）信号服从亚高斯分布，其激活函数为Фi （y i ）=αy i +y i |y 2
i |，而对于超高斯情况，其激活函数为Фi （y i ）=αy i +tan h （βy i ），这里α≥0，β≥2。

但是在实际中并不知道源信号是服从超高斯还是亚高斯分布，如果按照某一种统计特性来选取非线性函数，对源信号的分离效果差，甚至根本无法分离。

提出一种基于支持向量机（Support Vector Machine ）对分离信号进行概率密度函数的回归估计，得到一组稀疏的解系数，配合选定的核函数得出对分离信号概率密度函数的显式表达，并以此构造激活函数应用到FastICA 算法中，有效实现超高斯和亚高斯混合信源的分离。

2基于支持向量机
（SVM ）的概率密度函数回归方法常用的概率密度估计方法分为参数估计和非参数估计两
大类。

参数估计方法具有较大的局限性；而非参数估计具有更
广的应用范围，其中又以Parzen 窗方法最为常用。

但该方法估计概率密度要用到所有的训练样本，当训练样本非常多时，这种方法不实用[6]。

为了解决这个问题，就必须寻找稀疏解的办法。

而支持向量机（SVM ）的解只和训练样本中的支持向量有关，因此它提供了一种良好稀疏解的方法。

若考虑损失函数为ε-的不敏感函数L ε（d ，
y ），其中d 为期望输出，y 为实际输出，
即：L ε（d ，
y ）=0if |d-y|≤ε
|d-y|-εotherwis ≤
e
（4）SVM 的基于结构风险最小化（SRM ）原则为求使经验风险1
n
n
i =1ΣL ε
（d i
，
y i
）与结构风险之和最小的u ，即：min R SVM （C ）=C 1n n
i =1ΣL ε（d i +y i ）+
12
‖u ‖2
（5）
概率密度估计就是求线性算子方程：
+∞
-∞
乙
θ（x-t ）p （t ）d t=F （x ）（6）
的解，其中θ（x ）=1，x >0
0，x ≤乙
0，p （·）是概率密度函数，F
（·）为其未知分布函数。

同时解必需满足如下条件：p （x ）≥0，
+∞
-∞
乙
p （x ）d x =
1。

在上述线性算子方程
（6）中，分布函数F （x ）未知，但是给出了一定的独立同分布样本x 1，…，x n ，利用这些样本，构造经验分布函数F n （x
）=1n n
i =1
Σθ
（x-x i
）。

从概率密度的定义求解线性算子方程Ap （x ）=
+∞
-∞
乙
θ（x-t ）p （t ）d t 寻找该方程的解与在F （x ，u ）=
∞
r =0
Σu r
Φr
（x ）的空间中寻找线性回归函数是等价的。

由于概率密
度函数通常为非线性函数，为使其线性可分，将样本空间X 奂R
n
映射到高维特征空间Φ（X ）奂R N
（Φ（·）为非线性映射，N 可能
为无穷）：X →Φ
（X ）以构造超平面。

因此上述问题就转化为用SVM 方法对（x 1，F n （x 1）），…，（x n ，F n （x n
））在像空间中进行回归估计的问题，其解为
u =n *
i =0
Σα*
i Φ
（x *
i ）（7）
其中n *为支持向量数目，x *
i 为支持向量，α*
i 为其对应的非零拉
格朗日系数，Φ
（x ）是样本空间到高维空间映射的像。

定义像空间中的核函数为[6]：
K （x i ，x j ）=∞
r =0ΣΦr （x i ）Φr （x j
）（8）
交叉核函数为：K （x i ，x ）=∞
r =0ΣΦr （x i ）φr （x
）（9）
其中φ为Φ的微分，利用SVM 方法解线性算子方程可转换为利用核函数和数据对（x 1，F n （x 1）），…，（x n ，F n （x n ））在象空间中进行回归估计，获得u ，并通过交叉核函数得到线性算子方程的解[6]：
p （x ，α）=n *
i =1Σα*
i k （x *
i ，x
）（10）
3基于SVM 概率密度估计的FastICA 算法
在对观测信号x 概率密度的SVM 估计当中可以选择多种
核函数，为了使估计出的概率密度函数曲线更平滑，取径向基函数为核函数，其形式为：K （x ，y ）=exp{-σ‖x-y ‖2
}，σ>0为指
数因子，得到关于观测信号的概率密度函数的形式为[7]：
p （x ，α）=n *
i =1Σα*
i exp{-σ‖x-x *
i ‖2
}
（11）
这样的形式跟高斯函数多项式非常相似，只是差一个常系数，而这一点可以由α*
i 来弥补。

盲信号分离中，源信号是未知的，没有训练样本，所以根据概率密度函数的估计方法并不能直接用来估计源信号的概率密度函数。

要构造激活函数只能用源信号的估计y 来实现。

直接对y 的概率密度进行估计，
看起来似乎可行；但是y 由迭代产生，增大了计算量，而且分离矩阵W 的初值是人为选定的，不能保证算法收敛。

本文采用直接估计白化信号z 的概率密度的方式，Fast-I －CA 中系统输出Y =WZ ，
由于Y 是源的估计S 赞，其概率密度函数为（取σ=1）[8]
：
p （y ）=|det W |-1
p z （W -1
y ）=n *
i =1
Σα*
i N （W i z *
i ，W i W T
i
）（12）
式中N （μ，δ）表示均值为μ，方差为δ的正态过程，将p （y ）=f （y ）代入式（11）参与盲分离计算可得到目标函数的近似最优估计。

若要提取多个信源，只要把上述算法重复运行，每次取不同的w i （0）即可提取出多个信源。

为了保证每次提取出来的都是不同的信号源，必须在重复算法前将分离矩阵进行正交化处理把已经提取的分量去掉。

算法全过程如下：
（1）把原始数据x 去均值并白化，得z ；
（2）用SVM 回归算法估计z 的概率密度函数p （z ）；（3）设M 为待提取独立分量的个数，开始计数c =1；（4）给w c （0）设置初始值，计算向量的2范数，要求‖w c （0）‖2=1；（5）将p （z ）和w c （k ）代入式（3）计算y 的概率密度函数p （y ），并置f （y ）=p
（y ）；（6）迭代：w c （k +1）=E {z f [w T
c （k ）z ]}-E {[f ′（w T
c （k ）z ）]w c （k ）}，其中集合平均可用时间平均代替；
（7）将得到的W 正交归一化处理；（8）若w c 未收敛，回到步骤（5）；
胡波平，何选森：盲源分离的SVM 概率密度函数估计算法
143
（9）令c=c +1，判断c 是否大于M ，是则算法结束，否则转
到步骤（4
）。

4仿真实验
为了验证以上算法性能，在Matlab 平台进行了仿真。

选择6个独立信号作为源信号：符号函数sign （cos （2π155*t ）），高频正弦信号sin （2π800t ），低频正弦信号sin （2π90t ），调制信号sin （2π300t +6×cos （2π60t ）），高斯随机信号randn （1，1500），超高斯信号[（mod （t ，23）-11）/9]5+3；其峭度分别为：0.9998，1.5，
1.504，1.9547，3，5.411，
这里所给的峭度是未经归一化的值，以峭度值3作为超、亚高斯信号的分界点。

混合信道A 服从[-1，1]区间均匀分布的随机数，以速率f =10kHz 对传感器信号采样，
获得观测信号x
（t ），采样长度1500点。

算法的收敛性能用串音误差来评价[8]，即：
CE =m
i=1
Σ[m
j =1Σ
|c ij ，k |ij ，k k
-1]+m j =1Σ[m
i=1
Σ|c ij ，k |ij ，k k
-1]（13）
其中C =[c ij ]=W ΣA ，Σ为白化阵。

仿真中分离矩阵的初值为0.5I ，I 为与W 阶数相同的单位阵。

从分离结果来看，各路信号波形得到了很好的重现，失真很小。

图3反映了分离信号的峭度随迭代次数的变化，在前6次迭代中，算法并未收敛，分离出来的信号没有达到各自的峭度值。

直到第7步各信号峭度基本稳定在各自的数值上，几乎没有波动。

从图4的串音误差图上可以看到在前11步迭代当中算法未收敛，误差起伏很剧烈，而后算法迅速收敛，误差很小。

在仿真实验中由于是超、亚高斯信号共存，该算法与传统FastICA 算法分离精度比较如表1。

从表1可以看出传统的FastICA 算法在只有200个样本点时性能不稳定，甚至可能发散，而且收敛后的精度也不高，在采样800点时算法才能够得到较好的分离精度。

而文中算法的稳定性和精度都有较大提高。

图1源信号
S o u r c e W a v e f o r
m
10-110-110-110-150-51050
050010001500
050010001500
050010001500
050010001500
500
1000
1500
5001000
1
500
Sample Number
图2
分离信号
20-250-520-220-220-250-5
050010001500
050010001500
050010001500
050010001500
050010001500
50010001
500
Sample Number
S e p a r a t e d W a v e f o r m
图3分离信号归一化峭度值曲线
图4
串音误差
表1
算法性能比较
Algorithm Identity
Originnal FastICA Originnal FastICA Algorithm of this paper
Sample number
200points 800points 200points
Item Kurtosis
Kurtosis Error Ratio/
（%）Kurtosis Error Ratio/
（%）Kurtosis Error Ratio/
（%）Source 10.99981.01281.301.0010.121.00160.18
Source 21.51.663410.891.53372.251.54853.23
Source 31.50411.722814.541.54722.871.5714.45
Source 41.95472.495127.652.08536.681.99572.10
Source 5
32.93892.042.92752.422.95241.59
Source 65.4115.78926.995.5011.665.44170.57
（下转166页）
5结束语
通过采用支持向量机对观测信号的概率密度函数进行学
习的方式，间接估计输出信号的概率密度函数，由此构造更为合理的激活函数，对以负熵为代价函数的FastICA 算法进行了改进，获得了较高的分离精度。

仿真结果表明该算法能够对超、亚高斯同时存在的混合独立信号进行分离，其收敛速度快且分离精度高。

该算法可直接对观测信号进行操作，无须判断信号分布的类型，具有较强的自适应能力和稳健性。

参考文献：
[1]Cardoso J F.Blind signal separation[J].Statistical Principles Proc ，
1998，86（10）：2009-2025.
[2]Boscolo R ，Vwani H P.Independent component analysis based on
nonparametric density estimation.[J].IEEE Trans on Neural Net －
works ，2004，15（1）：55-65.
[3]Vlassis N ，Motomura Y.Efficient source adaptivity in independent
component analysis.[J].IEEE Trans on Neural Networks ，2001，12（3）：559-565.
[4]Choi S ，Cichocki A ，Amari S.Flexible independent component anal －
ysis.[C]//IEEE Workshop on Neural Networks for Signal Process －ing ，Cambridge ，UK ，1998：83-92.
[5]Hyvarinen A.Fast and robust fixed-point algorithms for independent
component analysis[J].IEEE Trans on Neural Networks ，1999，10（3）：626-634.
[6]张炤.基于支持向量机的概率密度估计方法[J].系统仿真学报，2005，17（10）：2355-2357.
[7]Vapnik V N.统计学习理论的本质[M].张学工，
译.北京：清华大学出版社，2000.
[8]李小军.基于估计概率密度函数的独立分量分析方法[J].西安电子
科技大学学报，2005，32（4）：574-578.
（上接144页
）
（d ）对图（c ）检测的结果
（c ）增加2001个噪声点
（b ）增加667个噪声点（a ）原图图3
合成图像检测实例
（a ）原图像
（b ）对图（a ）的预处理结果（c ）本文算法对图（b ）的检测结果
图4实际图像检测实例
噪声比/（%
）020406080100120140160
RHT/s 0.01040.06870.30181.27703.43657.409417.594534.881273.4238
文献[5]/s 0.00830.01560.02280.02940.03680.04120.04540.04830.0517
本文算法/s 0.00350.00380.00450.00520.00600.00790.01090.01300.0147
表1
三种算法执行时间比较
圆心横坐标a
圆心纵坐标b 序号12345678
检测值110.8902173.122373.991829.9571111.151588.658950.8182189.0418
精确值11117374301118951189
检测值134.134173.695042.713142.9286158.984843.0166115.8636119.0997
精确值13474434315943116119
检测值10.173264.925416.043520.67459.897620.993843.228350.1089
精确值1065162110214350
圆半径r 表2
本文算法对图3（c ）的检测结果
检测50次、500次的平均时间。

用本文算法对图3（c ）进行50次检测均能准确提取其圆心坐标和半径，结果如表2所示，其中图3的坐标原点为矩形框左下角、坐标轴分别平行于矩形的两条邻边。

实验2对实际图像图4（a ）进行相关预处理后得到图4（b ），
用RHT 与本文算法分别对图4（b ）检测500次，RHT 的平均检测时间为4.4905s ，本文算法只需0.0765s 。

两种算法均能准确提取其圆心坐标和半径，本文算法的检测结果如图4（c ）所示。

5结束语
提出一种基于随机Hough 变换的快速多圆检测算法。

算
法只需随机采样一点，通过该点按照一定方向搜索满足文中性质4的另外两点，从而确定候选圆；采用快速的证据积累方法，
同时确保判断真圆的可靠性。

由性质1~性质3可知，本文算法能快速地检测出图像中的多个圆，也适用于残缺圆的检测，如待检测圆的部分区域出现凸块、凹块、缺块等情形。

参考文献：
[1]Xu L ，Oja E ，Kultanen P.A new curve detection method ：randomized
Hough transform （RHT ）[J].Pattern Recognition Letters ，1990，11
（5）：331-338.
[2]Xu L ，Oja E.Randomized Hough transform （RHT ）：basic mechanisms ，
algorithms ，and computational complexities[J].Computer Vision Graph －ics Image Process ：Image Understanding ，1993，57（2）：131-154.[3]束志林，
戚飞虎．一种新的随机Hough 快速圆检测算法[J]．计算机工程，2003，29（6）：87-88.
[4]王建锋，
吴庆标．一种随机Hough 变换检测圆的改进算法[J]．计算机工程与应用，2005，41（14）.
[5]黎自强，
滕弘飞．广义Hough 变换：多个圆的快速随机检测[J]．计算机辅助设计与图形学学报，2006，18（1）：27-33.。