基于平行因子分析的欠定混合矩阵估计算法

2020-01-08 08:22王少波李红光
探测与控制学报 2019年6期
关键词:张量数目信噪比

王少波,郭 英,眭 萍,李红光,杨 鑫

(空军工程大学信息与导航学院, 陕西 西安 710077)

0 引言

盲源分离是指从观测到的多源混合信号中分离并恢复出相对独立的源信号的过程[1]。盲源分离是现代信号处理的重要前沿研究领域之一,已经在通信、语音处理、机械故障分析、图像处理、生物医学、雷达及经济数据分析等领域得到广泛应用[2-4]。

盲源信号分离技术根据源信号数目和混合信号数目之间的关系可以分为:超定、正定和欠定。超定即观测信号数目大于源信号数目,正定即观测信号数目等于源信号数目,欠定即观测信号数目小于源信号数目,其分别对应着不同的解决方法,目前解决超定、正定盲源分离的方法[5]已经相当成熟,而欠定盲分离问题的解决方法还有待进一步的研究。“两步法”是近年来解决欠定盲源分离问题的一个有效方法,即首先估计混合矩阵,然后在混合矩阵已知的基础上恢复源信号。估计混合矩阵的精度将直接影响分离信号的质量,因此混合矩阵估计算法非常关键。

欠定条件下混合矩阵的估计,主要采用稀疏聚类方法[6-8]和张量分解法[9-13]。文献[8]选择时频域做为稀疏变换域,通过计算接收混合信号时频比的方差来检测单源区域,再利用k-均值聚类来完成混合矩阵的估计。基于稀疏聚类的混合矩阵估计方法,计算相对简单,估计精度较高,但是需要满足一定的稀疏条件,这限制了该方法的应用范围。文献[9]将压缩感知理论与平行因子分析相结合,实现了欠定盲源分离及DOA估计,但其只适用于具有特殊表示形式的源信号。文献[10]和文献[12]将平行因子分析应用到盲源分离问题中,将混合矩阵的估计问题转化为张量分解问题,该方法提高了估计的准确性和鲁棒性,但是进行CP分解的ALS算法对初值敏感且容易陷入局部极值,运算时间较长。文献[13]引入塔克分解先把张量压缩为较低维的核张量,然后运用交替最小二乘对核张量进行标准分解,得到混合矩阵的估计,减少了运算时间,但是初值敏感问题依然存在。基于张量分解的估计算法,克服了稀疏聚类算法稀疏性要求,通常需要利用信号的高阶统计量构造张量,将混合矩阵的估计问题转化为张量分解问题,适用于信源相互独立并且非高斯分布的情况,精度较高,但是算法的计算量较大,运算时间长,且对初值敏感,容易陷入局部收敛。目前研究一种有效的欠定条件下混合矩阵估计方法仍是一个技术难点。

针对现有算法在解决非稀疏信号的欠定混合矩阵估计中,存在的计算时间长、初值敏感且容易陷入局部收敛的问题,提出了基于平行因子分析的欠定混合矩阵估计算法。

1 欠定盲源分离模型及张量理论基础

1.1 欠定盲源分离模型

欠定盲源分离是指在观测数目小于源信号数目的条件下,从观测到的多源混合信号中分离并恢复出相对独立的源信号的过程。图1为欠定盲源分离系统原理图。

设M个相互独立且非高斯分布观测信号是来自于N个源信号的线性混合,则欠定盲源分离信号模型描述为:

X(t)=AS(t)+N(t)

(1)

式(1)中,X(t)=[x1(t),x2(t),…,xM(t)]T表示接收的M个观测信号,S(t)=[s1(t),s2(t),…,sN(t)]T表示输入的N个未知源信号,M

图1 欠定盲源分离系统模型Fig.1 Blind source separation system model

1.2 张量理论基础

一维数组,我们称之为向量;二维数组,我们称之为矩阵;三维数组以及多维数组,我们称之为张量[12]。张量由于其强大的数据分析能力,近年来在数据挖掘、机器学习等领域得到广泛应用。

下面介绍与张量有关的定义与定理。

定义1(向量的外积)[14]:向量a∈RI1×1与b∈RI2×1的外积是一个I1×I2维的矩阵C,定义为:

C=a·b=abT

(2)

定义2(秩1张量)[14]:如果三阶张量χ等于三个向量的外积,则它的秩为1。

定义3(CP分解)[14]:任何张量χ都允许分解为秩-1张量的总和。在三阶张量的情况下,CP分解定义为:

(3)

式(3)中,a∈RI1×1,b∈RI2×1,c∈RI3×1,F为张量χ的秩。

定义4(Kronecker积)[14]:I1×I2维矩阵A=[a1,…,aI2]与I3×I4维矩阵B的Kronecker积是一个维的矩阵,并有如下表达式:

(4)

定义5(Khatri-Ra积)[12]:I1×I0维矩阵A=[a1,a2,…,aI0] 和I2×I0维矩阵B=[b1,b2,…,bI0]的Khatri-Rao积A⊙B是一个I1I2×I0维的矩阵,并有如下表达式:

A⊙B=[a1⊗b1,…,aI0⊗bI0]

(5)

平行因子分析(PARAFAC)最重要的一个特征就是CP分解是具有唯一性的。它是采用平行因子模型进行数据分析的首要条件。下面给出PARAFAC模型分解唯一性的一个重要定理。

定理1(PARAFAC模型的分解唯一性)[13]:若χ∈RI1×I2×I3是一个具有如式的CP分解张量,其中A∈RI1×F、B∈RI2×F和C∈RI3×F,如果满足:

k(A)+k(B)+k(C)≥2F+2

(6)

则矩阵A,B和C在存在尺度模糊和列模糊的条件下是唯一的(也称为本质唯一)。

2 基于平行因子分析的欠定混合矩阵估计

本文利用通信信号的协方差矩阵构造三阶张量,将混合矩阵估计问题转化为了张量分解问题,重点针对现有张量分解算法存在的局部收敛和初始矩阵选择的问题提出改进,采用非迭代的方法确定交替最小二乘算法的初始迭代矩阵,为了减少运算时间,在迭代过程中采用标准线搜索加速收敛。

2.1 构造张量

本文采用文献[10]中的方法来构造三阶张量,具体方法如下:

对于零均值且互不相关的非平稳实信号,源信号在t时刻的协方差矩阵可表示为:

(7)

式(7)中,Dk=E{stst+tk}是对角矩阵,k=1,…,K。为了简化,先忽略噪声的影响。解决的问题是在M

可以将矩阵C1,…,Ck按如下方式压缩为张量C∈M×M×K:(C)ijk=(Ck)ij,i=1,…,M,j=1,…,M,k=1,…,K。通过(D)kf=(Dk)ff定义一个矩阵Dkf,得到:

(8)

将其改写为:

(9)

式(9)中,af和df分别代表矩阵A和D的列向量。

张量沿三个方向的切片展开为矩阵,C(1)∈M2×K,C(2)∈KM×M,C(3)∈MK×M:

(C(1))(i-1)M+j,k=(C(2))(k-1)K+i,j=(C(3))(j-1)M+k,i=Cijk

(10)

根据式(10),张量的三个切片展开矩阵可以写为

C(1)=(A⊙A)·DT

(11)

C(2)=(D⊙A)·AT

(12)

C(3)=(A⊙D)·AT

(13)

由定理1知,满足一定条件,张量的CP分解唯一存在。文献[10]给出了传感器数目和源信号数目使CP分解唯一的关系,如表 1所示。

表1 观测数M与所允许的源信号最大数目Nmax的关系Tab.1 The relationship between the number of observations M and the maximum number Nmaxof source signals

至此,将欠定混合矩阵的估计问题转化为三阶张量的CP分解问题。

2.2 ALS算法

作为解决PARAFAC模型的经典算法,ALS算法算法步骤如图 2所示。

图2 ALS算法流程Fig.2 The ALS algorithm flowchart

算法通过交替最小二乘误差γ,来获得混合矩阵。

(14)

式(14)中,‖·‖F代表F-范数,当矩阵A固定,D的估计值可以由最小二乘结果得到:

(15)

式(15)中,(·)+表示广义逆矩阵。

同理,分别利用式(12)与式(13),通过同样的方式可得

(16)

(17)

由于初始矩阵的选择对ALS算法的准确性及收敛路径有比较大的影响,本文采用非迭代的直接三线性分解法确定ALS算法初始矩阵。

2.3 直接三线性分解

直接三线性分解[16]算法是一种非迭代的求解平行因子模型的三维分解方法,但是其精度和可靠性比ALS算法要差,可以作为ALS算法的初始矩阵。

第一步 根据式(11)、式(12)和式(13),分别对C(1)、C(2)和C(3)进行奇异值分解,其中U、V取前F个奇异值矢量,W取前两个左奇异值矢量。

C(1)=UISIVIU=UI(I×R)

(18)

C(2)=UJSJVJV=UJ(J×R)

(19)

C(3)=UKSKVKW=UK(K×2)

(20)

第二步 利用式(21)、式(22)构造为样本矩阵G1、G2。

(21)

(22)

式中,wki为矩阵W中的元素,Ck∈RM×M为张量C的第k个切片矩阵。

第三步 使用QZ分解阵G1、G2组成的特征方程,L、M分别为方程的特征向量,可证明矩阵A和D分别可由式(23)、式(24)得到[14]。

G1Lλ2r=G2Lλ1rA0=UL-1

(23)

G1Mλ2r=G2Mλ1rA0=VM-1

(24)

至此,通过直接三线性分解得到了粗估的加载矩阵,将其作为ALS算法的初始矩阵,开始迭代。

ALS算法运算时间较长,本文使用标准线搜索加速收敛。

2.4 算法流程

本文采用直接三线性分解粗估混合矩阵,将其作为ALS算法初始矩阵,迭代过程中采用标准的线搜索加速收敛。

线搜索方法首先找到使目标函数减少的反向,然后计算沿着下降方向移动的步长。下降方向可以由各种不同的方法计算,如梯度下降法。

例如在t次迭代,根据算法预测一定的迭代步长ρ,对矩阵A进行更新,可以表示为:A=A+ρΔA。

本文采用的线加速方案为只对矩阵D进行线搜索:

D(new)=D(it-2)+RLS(D(it-1)-D(it-2))

(25)

针对矩阵A,利用式(17)采用最小二乘进一步估计。

算法步长的计算应该十分简单,如果它比相应的迭代需要更多的时间则没有意义。最简单的情况是,RLS被赋予固定值(在1.2~1.3之间)[17],或者被设置为(it)1/3[18]。本文的RLS设置为(it)1/3。

具体算法流程:

步骤1 根据直接三线性分解算法粗估加载矩阵A0,D0;

步骤2 利用式(15)和A0估计出D1,利用式(16)和A0及D1估计出A1,令it=2,则A(it-2)=A0,D(it-2)=D0,A(it-1)=A1,D(it-1)=D1;

步骤3 根据线加速公式(25)计算D(new),利用式(17)及A(it-1)D(it-1)计算A(new);

步骤4 分别利用A(new)、D(new)和A(it-1),D(it-1)及式(14)计算误差函数γ(new)和γ(it-1);

步骤5 若γ(new)<γ(it-1),A(it-1)=A(new),D(it-1)=D(new),γ(it-1)=γ(new);否则,直接执行步骤6;

步骤6 利用式(15)和式(16)及A(it-1),D(it-1);计算得A(it)和D(it):

(26)

(27)

3 仿真实验与分析

为说明算法的性能,算法的效果评估采用混合矩阵估计误差和迭代次数作为指标。

3.1 混合矩阵估计误差

本节通过采用不同调制方式信号来完成欠定条件下混合矩阵估计实验。

实验采用4路不同调制方式的辐射源信号:s1为跳频信号;s2为LFM信号;s3为QPSK信号;s4为FM信号。观测信号数目设为3,信号采样点数为4 096,信噪比变化范围为-5~30 dB,步进为5 dB,每个信噪比处蒙特卡洛仿真实验100次。

评价混合矩阵估计性能采用归一化均方误差:

(28)

实验仿真结果如图3所示,横坐标为信噪比,纵坐标为归一化均方误差。参考算法1采用K均值聚类算法[19],直接进行迭代估计出混合矩阵。参考算法2采用ALS算法[10],受初始值的选择及局部收敛的影响,在低信噪比条件下混合矩阵估计性能一般。参考文献[3]采用的AP聚类方法[20],将每个样本都视为潜在的聚类中心,通过迭代确定聚类中心。本文算法先用非迭代的方法确定初始矩阵,再进行迭代确定混合矩阵。从图 3可以看出,直接采用K均值聚类算法,接收阵元信噪比达到5 dB时, 就已达到极限,约为-15.5 dB,信噪比大于5 dB时,信噪比的增加无法提升混合矩阵的估计准确性;采用本文算法、参考算法2以及参考算法3,混合矩阵估计的准确性随着信噪比的增加而增加,在低信噪比条件下,本文算法估计准确性较参考算法2提升约3 dB,与参考算法3性能接近,在信噪比较高的条件下,本文算法的性能要明显优于参考算法3,略优于参考算法2。对100次蒙特卡洛实验得到的归一化均方误差求方差,得到仿真结果,如图4所示,横坐标为信噪比,纵坐标为100次实验得到的均方误差的方差。从图中可以看出,本文算法均方误差的方差位于ALS算法减小了约0.5,说明本文算法在估计的稳定性上有了较大的提升。

图3 均方误差随信噪比变化曲线Fig.3 Mean square error with SNR curve

3.2 迭代次数

本节通过设定不同的观测数目与信源数来完成算法迭代次数估计实验。实验首先设置观测数目及信源数,然后随机生成A和D,进行蒙特卡洛仿真实验100次,记录归一化均方误差ENMSE达到10 dB时,所需的平均迭代次数。

图4 均方误差的方差随信噪比变化曲线Fig.4 The variance of the mean square error varies with SNR curve

表2列出了两种算法运行一次的平均迭代次数,从表中可以看出在(M,N)取不同值条件下,本文算法迭代次数上较ALS算法减少约41.4%~84.3%,且混合矩形规模越小,减少幅度越大。

表2 (M,N)取不同值时算法的迭代次数Tab.2 Number of iterations of the algorithm when (M,N) take different values

4 结论

本文提出了基于平行因子分析的欠定混合矩阵估计算法。该算法利用信号的协方差矩阵构造三阶张量,采用直接三线性分解确定交替最小二乘算法的初始迭代矩阵,然后在迭代过程中采用标准线搜索加速收敛,最终实现张量分解得到混合矩阵。仿真实验表明,该方法不要求信源的稀疏性,较ALS算法估计精度可以提高约3 dB,迭代次数可以减少约41.4%~84.3%,是一种有效的欠定混合矩阵估计算法。

猜你喜欢
张量数目信噪比
两种64排GE CT冠脉成像信噪比与剂量对比分析研究
浅谈张量的通俗解释
基于经验分布函数快速收敛的信噪比估计器
移火柴
非负张量谱半径上下界的估计不等式
支持张量机算法优化研究综述
自跟踪接收机互相关法性能分析
基于原子分解的非局部结构张量
基于深度学习的无人机数据链信噪比估计算法
牧场里的马