0 引言
       假设盲源分离的源信号个数为J，接收传感器的个数为R，盲源分离可以分为超定盲源分离（J≥R）和欠定盲源分离（J＜R）两种情况。大多数盲分离算法都假设接收传感器个数不少于源信号个数，然而在实际应用中，接收传感器个数往往有限，有时会出现接收传感器小于接收源信号个数的欠定混合情况（Underdetermined Blind Source Separation，UBSS）。欠定盲源分离一般分为两步：(1)分离混合矩阵；(2)恢复源信号。本文只考虑对混合矩阵的估计。
       针对欠定盲源分离问题，大多数文献提出的算法是将观测信号在时域或频域稀疏化，这势必会产生庞大的计算量，并且应用范围局限于观测信号和源信号数量较少的情况。考虑到源信号一般均满足相互独立和具有时间结构等特性，L.De Lathauwer提出了二阶欠定盲辨识算法(Second-order Blind Identification of Underdetermined Mixtures，SOBIUM)[1]，该方法不要求源信号在时域或变换域是稀疏的，通过对观测信号的时延协方差矩阵组成三阶张量直接进行平行因子分解实现对混合矩阵的估计。TICHAVSKY P在SOBIUM的基础上提出了加权欠定混合矩阵盲分离算法[2]，该算法通过加权张量分解来完成混合矩阵的估计，提高了分离信号的信干比，但是SOBIUM的迭代收敛时间较长，而且可能产生局部收敛。
       针对以上问题，本文在SOBIUM方法的基础上，加入了增强线搜索算法（Enhanced Line Search，ELS) 。ELS可以显著改善最小二乘法的性能，降低局部收敛的风险，更重要的是减少了迭代次数,并且复杂度不高。
1 欠定盲分离与PARAFAC分解
1.1 瞬时欠定盲源分离模型
       瞬时混合模型下的欠定盲源分离，其含噪混合模型为：
       X(t)=AS(t)+N(t)    t=1，…，T                                                  (1)
       其中，X(t)=[x1(t)，x2(t)，…，xR(t)]T为R维接收信号矢量，A表示一个未知的J×R的混合矩阵，S(t)=[s1(t)，s2(t)，…，sR(t)]T为R维源信号矢量，N(t)=[n1(t)，…，nR(t)]T为R维噪声矢量，(*)T代表转置。在噪声不存在或者可以忽略不计的情况下，式(1)可以化简为：
       X(t)=AS(t)(2)
1.2 PARAFAC分解
       平行因子(Parallel Factor，PARAFAC)分析又叫标准分解，是三面或更高面阵低秩分解的总称。平行因子分解在多个应用领域发挥着有广泛的作用，远远超出了化学计量学。平行因子分析在信号处理和通信领域中的数据域和子空间域[1-2]表现出良好的实用性,观测数据被转换为张量形式进行运算。下面给出关于平行因子的定义：
       定义1：若矩阵A的任意kN个列线性独立,则最大kN的值称之为矩阵A的Kruskal秩,简称k秩。
       定义2：如果一个张量X∈RI×J×K等于三个向量a，b，c的外积，则这个张量的秩为1。
       定义3：一个三阶张量X∈RI×J×K可以写成秩为1的张量的最小数量的线性组合，叫作平行因子分解。这一最小数量(源信号数N)等于张量X的秩(可用于对源信号数的估计)即：

       式中ar、br、cr分别代表矩阵A∈CI×R、B∈CJ×R和C∈CK×R的第r列。其中xijk=ai bj ck，i=1，…，I，j=1，…，J，k=1，…，K。
平行因子分解的矩阵模式可以写为：
       X(1)=(B⊙C)AT，X(1)∈CJK×I
       X(2)=(C⊙A)BT，X(2)∈CKI×J
       X(3)=(A⊙B)CT，X(3)∈CIJ×K(4)
       平行因子的唯一性在文献[3-5]中进行了研究，可以总结为下面的定理：
       定理1：如果满足
       kA+kB+kC≥2R+2(5)
       则平行因子分解唯一。kA、kB、kC分别代表矩阵A、B、C的秩。R为源信号的个数。
       用平行因子分解解决欠定盲分离混合矩阵问题时，源信号个数J与接收传感器最大个数Rmax的关系如表1所示。

1.3 PARAFAC分解估计欠定混合矩阵
       若源信号为零均值且互不相关的非平稳信号，那么源信号在t时刻的二阶自相关矩阵可表示为：

       式中DN=Est s■■是块对角阵，n=1，…，N，N是分块的个数，矩阵A称为分块成型矩阵。式(3)中时间延时?子n可以为零，上标T表示转置。
    将矩阵组{RN}转为(d，d，M)维的三阶张量形式：

       其中θ代表一个影响矩阵A和D的所有元素的参数向量。

       式中R代表张量，R为源信号的数目，⊙表示张量的外积，{an}和{dn}分别为A和D的列向量，上标*代表共轭转置。
2 加权增强最小二乘法
2.1 交替最小二乘法算法
    张量的标准分解通常使用三线性交替最小二乘（Alternating Least Squares，ALS）算法实现。迭代过程中的代价函数为：

       ||·||_F表示Frobenius矩阵范数。ALS的目标是在每一步迭代中，使张量R与它的当前估计值的差的范数最小。用于平行因子分析模型拟合的 ALS 过程即在固定上次迭代获取的部分矩阵估计值基础上, 估计其他矩阵, 该交错映射形式的最小二乘回归过程循环下去, 直至收敛。矩阵A、B和C的估计可以表示为：

其中上标“+”代表Moore-Penrose伪逆。
2.2 增强的线搜索（ELS）
       通常在数据量非常大，或当两个因素几乎共线时，ALS的收敛性是非常缓慢的[6]。压缩和线搜索是应对收敛慢问题的两种解决方案。本文采用增强线搜索来加快ALS：

上式中上标(k)、(k-1)、(k-2)分别代表第k次，第k-1次，第k-2次迭代。令：

其中代表迭代的方向。松弛因子?籽表示迭代的步长。?籽的选择十分重要。同一个算法只改变?籽的值，迭代收敛的速度变化如图1所示。

图1  松弛因子?籽对迭代收缩的影响

2.3 平行因子的增强加权目标函数
       平行因子分解模型同独立分量分析模型一样，具有置换不确定性和排列不确定性。为了有效地解决这个问题而不牺牲算法的收敛性，本文采用如下收敛函数：

       其中随着迭代次数的增加?着趋于0。I表示与XJK×I相同维数的单位矩阵。式(15)可以写为：

       其中JK×I维矩阵T3、T2、T1和T0分别表示为：

       上标k和k-2为了简便已省略。定义Vec为矩阵矢量化符号，例如有矩阵A∈CI×J，则：

       则等式（16）等效于

       其中4×1维矢量代表共轭转置，式(18)对于复数和实数都适用[7]。IJK×4维矩阵T分别由T3、T2、T1和T0的列矢量组合而成，可以表示为：

3 数据实验：分离混合语音的混合矩阵
       本节通过数值仿真实验对本文算法(SO-WALS-ELS)与SOBIUM算法的性能做比较。混合矩阵估计的相对误差公式为：

其中，。假设它们的列向量均已单位化且消除了排列顺序的不确定性。
       实验中采用4路独立源信号混合成3路观测信号为例，4路语音从语音库中随机选取，采样率为16 kHz，取160 000点，混合方式为瞬时混合，H为混合矩阵。图2为源信号，图3为混合信号，图4为分离信号，图5为本文算法与SOBIUM算法的对比图。

图2  源信号

图3  混合信号

图4  分离信号

    从图2～图5可以看出，改进的算法分离出了大概原始信号，但分离信号的顺序和极性都发生了变化，这也是目前平行因子分解尚无法解决的问题。
   

图5  本文算法与SOBIUM算法的对比图

       根据公式，改进前的算法相对误差为0.055 9，改进后的相对误差为0.029 8。经过多次实验，改进后的方法比原方法具有更快的收敛速度，并且更精确。
4 结论
       本文提出了一种基于增强加权最小二乘法的欠定混合矩阵分离的新算法，适用于非平稳信号。首先，该算法将接收信号的空间协方差矩阵叠加成三阶张量，然后再对此三阶张量进行平行因子分解，最后利用增强加权最小二乘法完成混合矩阵估计。仿真实验结果表明：本文提出的算法具有比SOBIUM算法更好的分离效果和更好的鲁棒性，而且实现简单，可满足实际应用的要求。
参考文献
[1] LIEVEN De LATHAUWER.Blind identification of underde-termined mixtures by simultaneous matrix diagonalization[J].IEEE Transactions on Signal Processing，2008，56(3)：1096-1105.
[2] PETR Tichavsky,ZBYNEK Koldovsky.Weight adjusted tensor method for blind separation of underdetermined mixtures of nonstationary sources[J].IEEE Transactions on Signal Processing，2011，59：1037-1047.
[3] KRUSKAL J B.Three-way arrays：Rank and uniqueness of trilinear decompositions，with application to arithmetic complexity and statistics，Linear Algebra Appl.，1977(18)：95-138.
[4] STEGEMAN  A，SIDIROPOULOS N D.On kruskal’s uniqueness condition for the Candecomp/ PARAFAC decomposition[C].Linear Algebra and its Applications，2007
(420)：540-552.
[5] NICHOLAS D.SIDIROPOULOS，RASMUSBRO.On the Uni-queness of Multilinear Decomposition of N-way Arrays[J].J.Chemometrics，2000，14(3)：229-239.
[6] RAJIH M，COMON P，HARSHMAN R A.Enhanced line search：A novel method to accelerate PARAFAC[J].SIAM Journal on Matrix Analysis，2008，30(3)：1148-1171.
[7] NION D，DE LATHAUWER L.Line Search computation of the Block Factor Model for blind multi-user access in wireless communications[C].Advances in Wireless Communi-cations(SPAWC) France，Cannes，.IEEE Workshop on Signal Processing，2006.
[8] NION D，DE LATHAUWER I.An enhanced line Search scheme for complex-valued tensor decompositions.Application in DS-CDMA[J].Signal Processing，2008，3(88)：749-755.