基于张量正则分解的欠定盲源分离中的混合矩阵识别方法

摘要

本发明公开了一种基于张量正则分解的欠定盲源分离中的混合矩阵识别方法，主要解决现有技术在估计混合矩阵时受特定条件限制的问题。其实现步骤是：(1)对源信号进行采样得到观测数据；(2)利用观测数据的四阶累积量计算在不同时延下的四阶协方差矩阵；(3)将不同时延下的四阶协方差矩阵扩展成三阶张量的形式；(4)对三阶张量进行张量正则分解得到待识别混合矩阵的Khatri-Rao乘积矩阵；(5)利用特征值分解的方法对该乘积矩阵进行处理，得到混合矩阵的估计值。本发明具有识别精度高的优点，可用于语音、通信、雷达及生物医学领域源信号在时频混叠条件下的欠定盲源分离。

说明书

技术领域

本发明属于通信技术领域，特别涉及一种混合矩阵的识别方法，可用于语音、通信、雷 达及生物医学领域源信号在时频混叠条件下的欠定盲源分离中。

背景技术

盲源分离BSS是指在未知的传输通道和源信号的条件下，仅通过传感器收到的观测信号 来达到分离源信号的目的，该方法已经广泛地应用于语音信号处理、图像处理、雷达、通信 及生物医学等各个领域。作为盲源分离的经典算法，独立分量分析ICA及其扩展算法大多用 于解决观测信号数量等于或大于源信号数量条件下的问题，这种盲源分离称为正定或超定的 盲源分离，但在实际过程中，往往需要解决源信号数量小于观测信号数量的问题，即欠定盲 源分离UBSS。欠定盲源分离系统的线性瞬时模型为X(t)＝AS(t)+W(t)，其中X(t)∈C^M表 示观测信号，M为观测信号个数；S(t)∈C^P为未知的源信号，P为源信号个数；W(t)∈C^M表 示加性噪声；未知混合矩阵A＝[a₁,a₂,…,a_P]∈C^M×P。在欠定盲分离系统中，观测信号个数 小于源信号个数，即M<P。欠定条件下的混合矩阵盲识别，就是在未知的混合矩阵A和源 信号S(t)的条件下从观测信号X(t)中识别出混合矩阵，是盲源分离问题的一个难点。

目前，稀疏成分分析SCA是解决欠定盲源分离的问题的主要方法，大多数算法都是通 过“两步法”来完成盲源分离，第一步则是估计出未知的传输通道，即混合矩阵模型，然后 利用识别出的混合矩阵和稀疏分解的方法完成源信号的恢复，因此混合矩阵的识别在盲源分 离问题中非常关键，其精度影响到后续源信号的恢复。部分学者利用信号的稀疏性，采用聚 类的方法进行混合矩阵的识别，当源信号在时域上不满足稀疏性时，则利用傅立叶变换或小 波变换等工具将信号变换到稀疏的频域上，然后利用聚类或势函数的方法识别混合矩阵，例 如，NgutyenLin-Trung,ABelouchrani,KarimA-M.Separatingmoresourcesthansensorsusing  time-frequencydistributions.EURASIPJournalonAppliedSignalProcessing,2005,17,pp.2828- 2847，当源信号在时频域均混叠的情况下，该方法的性能则不理想。部分学者利用时频的方 法，例如，陆凤波,黄知涛,彭耿,等,“基于时频分布的欠定混叠盲分离”，电子学报,2011, 39(9),pp.2067-2072，该方法对观测信号进行时频处理，然后提取信号的自源时频点，利用自 源时频点构造张量模型并对该模型进行张量正则分解，从而完成混合矩阵的识别，但是在频 域交叠比较严重的情况下自源时频点的提取并不理想，因此会影响混合矩阵的识别性能；还 有一部分学者利用信号的统计特性，例如DeLathauwerL,CastaingJ,CardosoJ,“Fourth-order  cumulant-basedblindidentificationofunderdeterminedmixtures”,IEEETransactionsonSignal  Processing,2007,55(6),pp.2965-2973，该方法无需源信号满足稀疏的特性，只需要源信号是统 计独立的非高斯信号，在实际过程中，这个条件往往是容易满足的，但该算法在求解的过程 中，需要假设源信号具有相同符号的峭度，即反映振动信号分布特性的数值统计量，是归一 化的四阶中心矩，而在源信号先验知识不足的情况下该条件往往很难满足，从而影响混合矩 阵的识别精度。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种基于张量正则分解的欠定盲源 分离中的混合矩阵识别方法，以在无需特定条件的前提下，提高识别精度。

本发明的技术方案是：在接收传感器处对经过未知通道的非高斯独立统计源信号进行采 样，得到观测信号；利用观测信号的四阶累积量构建四阶协方差矩阵，并将四阶协方差矩阵 表示成三阶张量模型；采用张量正则分解对三阶张量模型进行求解，对求解后得到的矩阵进 行特征值分解，完成混合矩阵的识别。其实现步骤如下:

(1)在接收端对源信号进行采样，得到观测信号；

(2)计算观测信号的四阶协方差矩阵其中，τ₁＝0，τ₂＝0，τ₃为 整数且τ₃∈[0,R-1]，M为观测信号个数，R是大于P的正整数，取值为2*P，P为源信号的 个数；

(3)将四阶协方差矩阵Q(0,0,0)，Q(0,0,1)，…,Q(0,0,R-1)扩展成三阶张量 T_i,j,k＝[Q(0,0,k)]_i,j，1≤i,j≤M²，0≤k<R；

(4)对三阶张量T进行张量正则分解，得到待识别混合矩阵的Khatri-Rao乘积 源信号的四阶统计特性矩阵D∈C^R×P及A_Q的共轭矩阵

(5)将待识别混合矩阵的Khatri-Rao乘积A_Q的第e列元素b_e表示为矩阵的形式B_e， 其中，B_e的每个元素为：B_e[i,j]＝b_e((i-1)M+j),1≤i,j≤M，1≤e≤P，然后对B_e进行特征值 分解，其中最大的特征值对应的特征向量即为识别出的混合矩阵的第e列。

与现有技术相比，本发明具有以下优点：

第一，本发明利用观测信号的四阶统计特性解决欠定盲源分离中的混合矩阵识别的问题， 克服现有技术在混合矩阵识别时要求源信号满足稀疏性的缺点，使得本发明可以解决源信号 在时频均混叠的条件下的欠定盲源分离中混合矩阵的识别问题。

第二，本发明引入张量正则分解的方法对观测信号的四阶协方差矩阵扩展的三阶张量模 型进行求解，克服现有技术中对观测信号的自源时频点提取困难的缺点，使得本发明提高了 欠定盲源分离中混合矩阵的识别精度。

第三，本发明将不同时延下的四阶协方差矩阵扩展成三阶张量模型，然后再进行求解， 克服现有技术中需要源信号具有相同符号的峭度的缺点，从而可以解决源信号先验知识不足 的条件下欠定盲源分离中混合矩阵的识别。

附图说明

图1为本发明的实现流程图；

图2为本发明仿真实验中设置的4路源信号时域波形图；

图3为本发明仿真试验中设置的4路源信号时频图；

图4为用3个接收传感器对4路源信号的线性混合信号采样得到的3路观测信号；

图5为本发明和现有方法对图4仿真得到的混合矩阵估计性能随信噪比变化的曲线图。

具体实施方式

以下参照附图对本发明作进一步详细的描述。

参照图1，本发明实现步骤如下：

步骤1：在接收端对源信号进行采样得到观测信号。

M个传感器在t时刻对源信号进行等间隔采样，得到观测信号x_i(t)，其中，1≤i≤M， t∈[1,2,…,N]，N为采样数据长度。

步骤2：计算观测信号的四阶协方差矩阵。

(2.1)计算观测信号的四阶矩：

${\hat{m}}_{i,j,k,l}({\tau }_1,{\tau }_2,{\tau }_3)=\frac{1}{T}\underset{t=1}{\overset{N}{\Sigma }}{x}_i\left(t\right)x_j^{*}(t+{\tau }_1)x_k^{*}(t+{\tau }_2)x_l(t+{\tau }_3),$

其中，1≤i,j,k,l≤M，τ₁,τ₂,τ₃分别为第j路，第k路以及第l路观测信号的时延；

(2.2)计算观测信号的互相关：

计算第i路观测信号x_i(t)与第j路观测信号x_j(t)在时延τ₁下的互相关为：

$B_x\left({\tau }_1\right)=\frac{1}{T}\underset{t}{\overset{N}{\Sigma }}{x}_i\left(t\right)x_j^{*}(t+{\tau }_1),$

计算第j路观测信号x_j(t)与第l路观测信号x_l(t)在时延τ₃-τ₂下的互相关为：

$B_x({\tau }_3-{\tau }_2)=\frac{1}{T}\underset{t}{\overset{N}{\Sigma }}{x}_j^{*}(t+{\tau }_2)x_l(t+{\tau }_3),$

计算第i路观测信号x_i(t)与第k路观测信号x_k(t)在时延τ₂下的互相关为：

$R_x\left({\tau }_2\right)=\frac{1}{T}\underset{t}{\overset{N}{\Sigma }}{x}_i\left(t\right)x_k^{*}(t+{\tau }_2),$

第j路观测信号x_j(t)与第l路观测信号x_l(t)在时延τ₃-τ₁下的互相关为：

$R_x({\tau }_3-{\tau }_2)=\frac{1}{T}\underset{t}{\overset{N}{\Sigma }}{x}_j^{*}(t+{\tau }_1)x_l(t+{\tau }_3),$

第i路观测信号x_i(t)与第l路观测信号x_l(t)在时延τ₃下的互相关为：

$B_x\left({\tau }_3\right)=\frac{1}{T}\underset{t}{\overset{N}{\Sigma }}{x}_i\left(t\right)x_l(t+{\tau }_3),$

第j路观测信号x_j(t)与第k路观测信号x_k(t)在时延τ₂-τ₁下的互相关为：

$R_x({\tau }_2-{\tau }_1)=\frac{1}{T}\underset{t}{\overset{N}{\Sigma }}{x}_j^{*}(t+{\tau }_1)x_k^{*}(t+{\tau }_2),$

其中，*表示复数共轭；

(2.3)计算观测信号的四阶累积量：

$\left(\begin{array}{c}{C}_{i,j,k,l}({\tau }_1,{\tau }_2,{\tau }_3)={\hat{m}}_{i,j,k,l}({\tau }_1,{\tau }_2,{\tau }_3)-R_x\left({\tau }_1\right)R_x({\tau }_3-{\tau }_2)-\\ R_x\left({\tau }_2\right)R_x({\tau }_3-{\tau }_1)-R_x\left({\tau }_3\right)R_x({\tau }_2-{\tau }_1)\end{array}\right);$

(2.4)通过四阶累积量C_i,j,k,l(τ₁,τ₂,τ₃)计算观测信号的四阶协方差矩阵Q(τ₁,τ₂,τ₃)：

Q(τ₁,τ₂,τ₃)[M(i-1)+j，M(k-1)+l]＝C_i,j,k,l(τ₁,τ₂,τ₃)，

(2.5)四阶协方差矩阵Q(τ₁,τ₂,τ₃)中取τ₁＝0，τ₂＝0，τ₃∈[0,R-1]，获得不同时延 下的四阶协方差矩阵：Q(0,0,0)，Q(0,0,1)，…,Q(0,0,R-1)，其中R是大于P的正整数，取值为2*P， P为源信号的个数。

步骤3：将观测信号的四阶协方差矩阵Q(0,0,0)，Q(0,0,1)，…,Q(0,0,R-1)扩展成三阶张量 T_i,j,k＝[Q(0,0,k)]_i,j(1≤i,j≤M²,0≤k<R)；

(3.1)根据源信号是统计独立的特性，将四阶协方差矩阵Q(τ₁,τ₂,τ₃)表示为：

$Q({\tau }_1,{\tau }_2,{\tau }_3)=A_Q{C}_S({\tau }_1,{\tau }_2,{\tau }_3)A_Q^H,$

其中，a_p表示混合矩阵的第p列，为Kronecker乘积， C_S(τ₁,τ₂,τ₃)是P×P的对角矩阵；

(3.2)将四阶协方差矩阵Q(0,0,0)，Q(0,0,1)，…,Q(0,0,R-1)表示为：

$\left(\begin{array}{c}Q\left(\mathrm{0,0,0}\right)=A_Q{C}_S\left(\mathrm{0,0,0}\right)A_Q^H,\\ Q\left(\mathrm{0,0,1}\right)=A_Q{C}_S\left(\mathrm{0,0,1}\right)Q_Q^H,\\ .\\ .\\ .\\ Q(\mathrm{0,0},R-1)=A_Q{C}_S(\mathrm{0,0},R-1)A_Q^H,\end{array}\right)$

其中，表示A_Q的酉矩阵；

(3.3)将四阶协方差矩阵Q(0,0,0)，Q(0,0,1)，Q(0,0,2)，…,Q(0,0,R-1)扩展成三阶张量T，其 中，T的第(i,j,k)个元素为T_i,j,k＝[Q(0,0,k)]_i,j，1≤i,j≤M²，0≤k<R-1，Q(0,0,0)为T的第一 维切片，Q(0,0,1)为T的第二维切片，Q(0,0,2)为T的第三维切片，以此类推，Q(0,0,R-1)为T的 第R维切片。

步骤4：对三阶张量T进行张量正则分解：

(4.1)根据三阶张量的定义，求解目标矩阵1≤r≤3，使得代价函数 f(U⁽¹⁾,U⁽²⁾,U⁽³⁾)最小，其中代价函数f(U⁽¹⁾,U⁽²⁾,U⁽³⁾)的定义如下式：

其中，1≤i,j≤M²，0≤k<R-1，和分别为U⁽¹⁾、U⁽²⁾和U⁽³⁾的第r个列向量，和分别为和的第i、j、k个元素；

4.2)采用线性搜索迭代最小二乘LS_ALS算法，对上述代价函数f(U⁽¹⁾,U⁽²⁾,U⁽³⁾)进行优 化，得到目标矩阵U^(r)，1≤r≤3，其中U⁽¹⁾为待识别混合矩阵的Khatri-Rao乘积U⁽²⁾为A_Q的共轭矩阵U⁽³⁾为源信号的四阶统计特性矩阵D∈C^R×P。

步骤5：对待识别混合矩阵的Khatri-Rao乘积A_Q进行特征值分解，得到混合矩阵估计值：

(5.1)将待识别混合矩阵的Khatri-Rao乘积A_Q的每一列b_e表示为矩阵的形式B_e， 其中，B_e的每个元素表示为：B_e[i,j]＝b_e((i-1)M+j),1≤i,j≤M，1≤e≤P；

(5.2)对B_e进行特征值分解，其中最大的特征值对应的特征向量即为识别的混合矩阵的 第e列，1≤e≤P。

本发明的效果可通过以下仿真进一步说明：

为验证本发明的有效性和正确性，采用三种已有的混合矩阵识别方法与本发明方法同时 对某混合矩阵进行识别。所有仿真实验均在Windows8.1操作系统下采用Matlab2012b软件 实现。

1)仿真参数

采样率200MHz，采样点数为1024。4个源信号是时域上和频域上混叠的LFM信号，各 个源信号归一化频率范围为[0.5,0]、[0,0.4]、[0.5,0.24]和[0.35,0.15]，其时域波形图 和时频分布如图2和图3所示。

根据空间自由传输模型，混合矩阵A中每一个元素的定义为：

a_mp＝exp(2πj(α_mcos(θ_p)cos(φ_p)+β_mcos(θ_p)sin(φ_p)))，

其中，α_m＝(R_a/λ)cos(2π(m-1)/M)，β_m＝(R_a/λ)sin(2π(m-1)/M)，R_a/λ＝0.55。信号 入射方位角为θ₁＝3π/10、θ₂＝3π/10、θ₃＝2π/5和θ₄＝0；俯仰角φ₁＝7π/10、φ₂＝9π/10、φ₃＝3π/5 和φ₄＝4π/5。

为了评价算法对混合矩阵的识别性能，采用平均相对误差E_A作为评价因子：

$E_A=E\left\{\frac{{\left|\right|A-\hat{A}\left|\right|}_F}{{\left|\right|A\left|\right|}_F}\right\},$

式中，表示识别的混合矩阵，||·||_F这里表示的是F范数。

2)仿真内容

在信噪比变化为0-35dB的范围，采用基于聚类的混合矩阵识别算法Cluster_based、基于 时频分布的混合矩阵识别算法TFDs_based、基于联合对角化的四阶累积量混合矩阵估计算法 FOOBI以及本发明方法对图4中的3路观测信号做100次MonteCarlo仿真实验，得到混合矩 阵的性能随信噪比的变化曲线，如图5所示。

从图5可见，本发明方法在混合矩阵识别精度上要优于其他三种算法，从而验证了本发 明方法对欠定盲源分离中的混合矩阵识别的有效性和正确性。