CN104934041B - 基于多目标优化联合块对角化的卷积盲信号分离方法 - Google Patents

Abstract

本发明公开一种基于多目标优化联合块对角化的卷积盲信号分离方法，主要解决现有技术无法从卷积混叠信号中准确分离出所有源信号的问题。其实现步骤为：(1)获取观测数据；(2)计算观测数据的二阶时延相关矩阵；(3)构造块对角化矩阵，并进行子阵划分；(4)建立关于块对角化矩阵的多目标优化模型；(5)根据多目标优化模型，估计块对角化矩阵；(6)判断前后两次块对角化矩阵估计误差的差值绝对值是否大于迭代终止阈值，若是，则输出块对角化矩阵，否则返回步骤(5)；(7)利用块对角化矩阵从观测信号中分离出源信号。本发明能从卷积混叠信号中准确分离出所有源信号，且复杂度低，分离效率高，可用于对语音信号、通信信号的处理。

Description

基于多目标优化联合块对角化的卷积盲信号分离方法

技术领域

本发明涉及信号处理技术领域，更进一步涉及卷积盲信号分离方法，可用于实现语音信号、通信信号的处理。

背景技术

盲信号处理在无线通信、雷达和声呐、图像恢复、语音增强等各个领域具有广泛应用。盲信号分离是盲信号处理的一个重要组成部分，它的唯一假设条件是源信号相互统计独立，这对于自然界的绝大部分情况都是近似成立的，而且盲信号分离技术不需要任何关于源信号和传输信道的先验信息，使得该项技术成为一种非常实用的信号处理方法。尤其当很难或根本无法建立从多个信源到传感器之间的传输信道模型时，盲信号分离就成为唯一可行的信号处理方法。

Buchner H,Aichner R和Kellermann W在论文“A generalization of blindsource separation algorithms for convolutive mixtures based on second-orderstatistics”(IEEE transaction speech and audio process,2005,13(1):120-134)中提出了一种非正交联合块对角化算法。该方法使用最速下降法优化一个针对块对角化的信息论代价函数，提高了信号分离的性能，但是仍然存在收敛速度慢的缺陷。

西安电子科技大学所提出的专利申请“基于三因子迭代联合块对角化的时域混叠盲信号分离方法”(专利申请号：201410008337.7，公布号CN 103871422A)公开了一种基于三因子迭代联合块对角化的时域混叠盲信号分离方法。该方法首先将从传感器接收的时域卷积混叠信号重新排列成观测数据矢量，然后计算观测数据矢量在指定时延下的二阶时延相关矩阵组，然后利用基于三因子迭代的非正交联合块对角化方法计算出真实块混叠矩阵的块对角矩阵，求其伪逆矩阵，计算出源信号估计信号，从而实现时域卷积混叠信号盲源分离。该方法虽然利用构造代价函数求解三组矩阵因子，降低了计算复杂度，克服了对噪声敏感、易产生奇异解等缺点，但是仍然存在的不足之处是，该方法为二次型优化方法，对块对角化矩阵没有任何约束，很容易得到退化解，导致无法实现对所有源信号的分离。

发明内容

本发明的目的在于克服上述现有方法的不足，提出一种基于多目标优化联合块对角化的卷积盲信号分离方法，以实现从卷积混叠信号中准确的分离出所有源信号。

本发明的基本思想是：通过传感器接收到的卷积混叠信号重新排列得到观测数据，计算出指定时延下观测数据的二阶时延相关矩阵；通过对构造的块对角矩阵进行子阵划分建立多目标优化模型；通过对块对角矩阵迭代得到块对角矩阵的最优解，从而从观测信号中分离出源信号。

本发明的具体实现步骤包括如下：

(1)从传感器接收卷积混叠信号数据，将该卷积混叠信号数据根据传感器编号顺序重新排列成观测数据矢量；

(2)计算观测数据矢量在指定时延下的二阶时延相关矩阵R(τ_k)，τ_k表示指定的第k个时延，k＝1,…,Z，Z表示指定时延的总数目；

(3)构造块对角化矩阵，并进行子阵划分：

(3a)构造M×N维块对角化矩阵B，其中，M表示二阶时延相关矩阵R(τ_k)的行数，N表示待分离源信号向量的行数；

(3b)令构造的块对角化矩阵B＝[I,0]^T，得到初始块对角化矩阵，其中I表示单位矩阵，(·)^T表示转置；

(3c)将初始块对角化矩阵按列分成r个子矩阵，r等于待分离源信号的个数q，每个子矩阵的维数为M×n，其中N表示待分离源信号向量的行数；

(4)按照下式，建立非正交联合块对角化多目标优化模型：

其中，表示块对角化矩阵B的估计误差，B表示块对角化矩阵，Z表示指定时延的总数目，r表示子矩阵的总数目，n表示每个子矩阵的总列数，b_lm表示第l个子矩阵的第m列，b_l′m′表示第l′个子矩阵的第m′列，R(τ_k)表示第k个时延对应的二阶时延相关矩阵，b_j表示块对角化矩阵B的第j个列向量，min表示取最小化操作，max表示取最大化操作，det(·)表示求矩阵行列式操作，|·|²表示取平方操作，(·)^H表示复共轭操作，s.t.表示约束条件；

(5)利用非正交联合块对角化多目标优化模型，估计块对角化矩阵：

(5a)令m＝1，m表示子矩阵的列数；

(5b)令l＝1，l表示子矩阵的个数；

(5c)按照下式计算海森矩阵Q_l：

其中，B_l表示块对角化矩阵B的第l个子矩阵，Z表示指定时延的总数目，R(τ_k)表示第k个时延对应的二阶时延相关矩阵，B_(l)表示块对角化矩阵B删除第l个子矩阵后剩余元素构成的矩阵，(·)^H表示复共轭操作；

(5d)按照下式计算正交投影矩阵P_lm：

其中，b_lm表示第l个子矩阵的第m列，B_(lm)表示块对角化矩阵B删除第l个子矩阵的第m列元素后剩余元素构成的矩阵，I表示单位矩阵，[·]^-1表示矩阵的求逆操作；

(5e)计算海森矩阵Q_l的条件数ξ；

(5f)判断海森矩阵Q_l的条件数ξ是否大于设定条件数阈值θ：若是，则执行步骤(5g)-(5h)，否则，执行步骤(5i)；

(5g)按照下式计算中间矩阵C：

C＝U^HP_lmU

其中，U表示海森矩阵Q_l的M-N+n个最小的特征值对应的特征向量矩阵，M表示块对角化矩阵B的行数，N表示块对角化矩阵B的列数，n表示每个子矩阵的总列数；

(5h)计算第l个子矩阵的第m列，并执行步骤(5j)：

b_lm＝Uw

其中，w表示中间矩阵C的最大特征值对应的特征向量；

(5i)计算矩阵对(P_lm,Q_l)最大的广义特征值对应的特征向量，作为第l个子矩阵的第m列，并执行步骤(5j)；

(5j)将l的值增加1，判断此时l的值是否小于子矩阵的总数目r：若是，则返回步骤(5c)，否则，执行步骤(5k)；

(5k)将m的值增加1，判断此时m的值是否小于每个子矩阵的总列数n：若是，则返回步骤(5b)，否则，执行步骤(6)；

(6)判断当前块对角化矩阵B的估计误差与前一次的块对角化矩阵B的估计误差的差值绝对值是否大于设定迭代终止阈值若是，则返回步骤(5)，否则，将得到的块对角化矩阵B输出，并执行步骤(7)；

(7)使用块对角化矩阵B从观测信号中分离出源信号。

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

第一，由于本发明建立的非正交联合块对角化多目标优化模型，考虑了块对角化矩阵的条件数，克服了现有技术对块对角化矩阵没有任何约束，容易得到退化解，导致无法实现对所有源信号的分离的不足，使得本发明能够从卷积混叠信号中准确的分离出所有源信号。

第二，由于本发明首先对块对角化矩阵进行子阵划分，然后根据建立的非正交联合块对角化多目标优化模型对子矩阵按列更新，克服了现有技术使用最速下降法优化构造的代价函数，收敛速度慢的缺点，使得本发明具有计算时间少，分离效率高的优点。

附图说明

图1为本发明的流程图；

图2为用本发明和现有一种信号分离方法的成功执行概率随信噪比变化对比图；

图3为用本发明和现有两种信号分离方法在两种信噪比下的平均收敛性对比图；

图4为本发明仿真使用的三个语音源信号的波形图；

图5为用本发明和现有一种信号分离方法恢复的三个语音信号的波形图，其中，图5(a)为用现有信号分离方法恢复的三个语音信号的波形图，图5(b)为用本发明方法恢复的三个语音信号的波形图。

具体实施方式

下面结合附图对本发明做进一步的描述：

本实施例是基于一个酒会场景，对麦克风接收到的客人的谈话内容和一大堆背景噪声，利用本发明分离出客人的谈话语音。在本实施例中，传感器为麦克风，接收到的卷积混叠信号为语音信号。

参照图1，本发明的具体实施步骤如下：

步骤1，获取观测数据。

从传感器接收卷积混叠信号数据，为了有序处理数据，对使用的传感器进行顺序编号，将该卷积混叠信号数据根据传感器编号顺序重新排列成观测数据矢量X_i(t)。

步骤2，计算二阶时延相关矩阵。

计算观测数据矢量X_i(t)在指定时延下的二阶时延相关矩阵R(τ_k)：

其中，G表示观测信号样本数据段的总数目，τ_k表示设定的第k个时延，k＝1,…,Z，Z表示指定时延的总数目，考虑到人类语音平稳时间约为几十微秒，设定的时延要保证语音信号的短时相关性，同时，时延间隔应该尽量大，以保证在各个时延下的相关矩阵有所差别。

步骤3，构造块对角化矩阵，并进行子阵划分。

3a)构造M×N维块对角化矩阵B，其中，M表示二阶时延相关矩阵R(τ_k)的行数，N表示待分离源信号向量的行数；

3b)令构造的块对角化矩阵B＝[I,0]^T，得到初始块对角化矩阵，其中I表示单位矩阵，(·)^T表示转置；

3c)将初始块对角化矩阵按列分成r个子矩阵，r的取值为待分离源信号的个数，每个子矩阵的维数为M×n，其中

步骤4，建立非正交联合块对角化多目标优化模型。

4a)计算块对角化矩阵的估计误差：

其中，表示块对角化矩阵B的估计误差，B表示块对角化矩阵，Z表示指定时延的总数目，r表示子矩阵的总数目，n表示每个子矩阵的总列数，b_lm表示第l个子矩阵的第m列，b_l′m′表示第l′个子矩阵的第m′列，R(τ_k)表示第k个时延对应的二阶时延相关矩阵，｜·｜²表示取平方操作，(·)^H表示复共轭操作；

4b)计算块对角化矩阵的条件数η_B：

η_B＝det(B^HB)，

其中，det(·)表示求矩阵行列式操作，B表示块对角化矩阵，(·)^H表示复共轭操作；

4c)根据块对角化矩阵的估计误差和块对角化矩阵的条件数η_B，建立如下非正交联合块对角化多目标优化模型：

其中，b_j表示块对角化矩阵B的第j个列向量，min表示取最小化操作，max表示取最大化操作，s.t.表示约束条件。

步骤5，估计块对角化矩阵。

利用步骤4c)建立的非正交联合块对角化多目标优化模型，对块对角化矩阵B每个子矩阵的列向量依次进行更新，估计块对角化矩阵，具体步骤包括：

5a)令m＝1，m表示子矩阵的列数；

5b)令l＝1，l表示子矩阵的个数；

5c)按照下式计算海森矩阵Q_l：

其中，B_l表示块对角化矩阵B的第l个子矩阵，Z表示指定时延的总数目，R(τ_k)表示第k个时延对应的二阶时延相关矩阵，B_(l)表示块对角化矩阵B删除第l个子矩阵后剩余元素构成的矩阵；

5d)按照下式计算正交投影矩阵P_lm：

其中，b_lm表示第l个子矩阵的第m列，B_(lm)表示块对角化矩阵B删除第l个子矩阵的第m列元素后剩余元素构成的矩阵，I表示单位矩阵，[·]^-1表示矩阵的求逆操作；

5e)按照下式，计算海森矩阵Q_l的条件数ξ：

ξ＝||Q₁||·||Q_l ^-1||

其中，[·]^-1表示矩阵的求逆操作，||·||表示求范数操作；

5f)设定条件数的阈值θ，θ的数量级一般为10³，将海森矩阵Q_l的条件数ξ与设定条件数的阈值θ进行比较：若ξ＞θ，则海森矩阵Q_l是病态的，执行步骤5g)-5h)，否则，执行步骤5i)；

5g)按照下式计算中间矩阵C：

C＝U^HP_lmU

其中，U表示海森矩阵Q_l的M-N+n个最小的特征值对应的特征向量矩阵，M表示块对角化矩阵B的行数，N表示块对角化矩阵B的列数，n表示每个子矩阵的总列数；

5h)计算第l个子矩阵的第m列向量：b_lm＝Uw，并执行5j)，其中，w表示中间矩阵C的最大特征值对应的特征向量；

5i)按如下步骤计算矩阵对(P_lm,Q_l)最大的广义特征值对应的特征向量，作为第l个子矩阵的第m列，并执行5j)：

5i.1)计算矩阵对(P_lm,Q_l)的广义特征值和对应的特征向量：

[VD]＝eig(P_lm,Q_l)

其中，D表示矩阵对(P_lm,Q_l)的广义特征值构成的对角矩阵，V表示矩阵对(P_lm,Q_l)的广义特征值对应的特征向量构成的矩阵，eig表示求广义特征值操作；

5i.2)将所述V的第一列作为矩阵对(P_lm,Q_l)最大的广义特征值对应的特征向量；

5j)将l的值增加1，判断此时l的值是否小于子矩阵的总数目r：若是，则返回5c)，否则，执行5k)；

5k)将m的值增加1，判断此时m的值是否小于每个子矩阵的总列数n：若是，则返回5b)，否则，执行步骤6。

步骤6根据需求设定迭代终止阈值 的数量级一般为10^-2，计算当前块对角化矩阵B的估计误差和前一次的块对角化矩阵B的估计误差，判断前后两次估计误差的差值绝对值是否大于设定迭代终止阈值若是，则返回步骤5，否则，将得到的块对角化矩阵B输出，并执行步骤7。

步骤7，分离源信号。

将块对角化矩阵B的伪逆矩阵与观测数据矢量X_i(t)相乘，得到从传感器接收的卷积混叠信号中分离出的源信号：

其中，y(t)表示t时刻的源信号，表示求矩阵的伪逆操作。

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

1.仿真条件：

本发明的仿真运行系统为Intel(R)Core(TM)i7-2600CPU 650@3.40GHz，32位Windows操作系统，仿真软件采用MATLAB(R2013a)。

2.仿真内容与结果分析：

仿真一：人工合成数据仿真

按照下式产生二阶时延相关矩阵R(τ_k)：

其中，块混叠矩阵A，噪声矩阵E_k，以及块对角矩阵D_k的r个对角块中的元素都是随机产生的复数，且服从均值为0，方差为1的复高斯分布，k＝1,…,Z，Z表示二阶时延相关矩阵的总数目，σ²表示噪声方差，(·)^H表示复共轭操作。

定义干扰信号比ISR为：

其中，G＝B^HA表示全局变换矩阵，具有r×r个子块，B表示块对角化矩阵，A表示块混叠矩阵，G_ij表示全局变换矩阵G的第(i,j)个子块，G_ip表示全局变换矩阵G的第(i,p)个子块，G_pj表示全局变换矩阵G的第(p,j)个子块，表示求取得最大值时的p的取值，表示求F范数操作。

设定M＝15，N＝12，Z＝20，r＝3，信噪比用本发明和现有信号分离方法一：JBD-NU各进行500次独立实验，当得到的干扰信号比ISR小于设定阈值μ＝-20dB，认为该次分离是成功的，得到用本发明和现有一种信号分离方法的成功执行概率随信噪比变化对比图如附图2所示。图2中的纵坐标表示成功执行概率，横坐标表示信噪比。

设定M＝N＝12，Z＝20，r＝3，信噪比分别设为SNR＝8dB和SNR＝16dB，用本发明和现有方法一：JBD-NU和方法二：JBD-S两种信号分离方法，各进行500次独立实验，得到在两种信噪比下的平均收敛性对比图如附图3所示。图3中的纵坐标表示干扰信号比，横坐标表示迭代次数。

由图2可以看出无论在怎样的噪声条件下，本发明方法都具有更高的成功执行概率。

由图3可以看出采用本发明方法在两种信噪比下都具有最快的收敛速度，且具有更好的干扰信号比性能。

仿真二：语音源信号分离仿真

对附图4中所示的三个语音源信号，将它们通过模拟的多径效应传输，利用6个传感器接收混合信号，其中多径效应的模拟采用3阶FIR滤波器，滤波器系数随机产生，且服从标准正态分布。将每个接收信号进行5阶重排列，然后利用重排列后的接收信号计算21个时延相关矩阵。

设置对应的仿真参数为：M＝5×6＝30，N＝3×(5+3)＝24，Z＝21，r＝3，设定信噪比SNR＝20dB。用本发明和现有信号分离方法一：JBD-NU恢复三个语音信号，结果如图5，其中5(a)为用现有信号分离方法一：JBD-NU恢复的三个语音信号的波形图；5(b)为用本发明方法恢复的三个语音信号的波形图。图5(a)和图5(b)中的纵坐标表示信号幅值，横坐标表示样本数。

由图5(a)可以看出，采用JBD-NU方法恢复的3个信号都与图4中第2个源信号波形相似，这说明JBD-NU方法事实上收敛到了退化解，它不能分离出所有源信号。

由图5(b)可以看出采用本发明方法恢复的3个语音信号分别与图4中的3个源信号波形相似，这说明本发明方法成功消除了退化解，它能够收敛到良态解，从而分离出所有源信号。

Claims (5)

1.一种基于多目标优化联合块对角化的卷积盲信号分离方法，包括如下步骤：

(1)从传感器接收卷积混叠信号数据，将该卷积混叠信号数据根据传感器编号顺序重新排列成观测数据矢量；

(2)计算观测数据矢量在指定时延下的二阶时延相关矩阵R(τ_k)，τ_k表示指定的第k个时延，k＝1,...,Z，Z表示指定时延的总数目；

(3)构造块对角化矩阵，并进行子阵划分：

(3a)构造M×N维块对角化矩阵B，其中，M表示二阶时延相关矩阵R(τ_k)的行数，N表示待分离源信号向量的行数；

(3b)令构造的块对角化矩阵B＝[I,0]^T，得到初始块对角化矩阵，其中I表示单位矩阵，(·)^T表示转置；

(3c)将初始块对角化矩阵按列分成r个子矩阵，r等于待分离源信号的个数q，每个子矩阵的维数为M×n，其中N表示待分离源信号向量的行数；

(4)按照下式，建立非正交联合块对角化多目标优化模型：

其中，表示块对角化矩阵B的估计误差，B表示块对角化矩阵，Z表示指定时延的总数目，r表示子矩阵的总数目，n表示每个子矩阵的总列数，b_lm表示第l个子矩阵的第m列，b_l′m′表示第l′个子矩阵的第m′列，R(τ_k)表示第k个时延对应的二阶时延相关矩阵，b_j表示块对角化矩阵B的第j个列向量，min表示取最小化操作，max表示取最大化操作，det(·)表示求矩阵行列式操作，|·|²表示取平方操作，(·)^H表示复共轭操作，s.t.表示约束条件；

(5)利用非正交联合块对角化多目标优化模型，估计块对角化矩阵：

(5a)令m＝1，m表示子矩阵的列数；

(5b)令l＝1，l表示子矩阵的个数；

(5c)按照下式计算海森矩阵Q_l：

其中，B_l表示块对角化矩阵B的第l个子矩阵，Z表示指定时延的总数目，R(τ_k)表示第k个时延对应的二阶时延相关矩阵，B_(l)表示块对角化矩阵B删除第l个子矩阵后剩余元素构成的矩阵，(·)^H表示复共轭操作；

(5d)按照下式计算正交投影矩阵P_lm：

其中，b_lm表示第l个子矩阵的第m列，B_(lm)表示块对角化矩阵B删除第l个子矩阵的第m列元素后剩余元素构成的矩阵，I表示单位矩阵，[·]^-1表示矩阵的求逆操作；

(5e)计算海森矩阵Q_l的条件数ξ；

(5f)判断海森矩阵Q_l的条件数ξ是否大于设定条件数阈值θ：若是，则执行步骤(5g)-(5h)，否则，执行步骤(5i)；

(5g)按照下式计算中间矩阵C：

C＝U^HP_lmU

其中，U表示海森矩阵Q_l的M-N+n个最小的特征值对应的特征向量矩阵，M表示块对角化矩阵B的行数，N表示块对角化矩阵B的列数，n表示每个子矩阵的总列数；

(5h)计算第l个子矩阵的第m列，并执行步骤(5j)：

b_lm＝Uw

其中，w表示中间矩阵C的最大特征值对应的特征向量；

(5i)计算矩阵对(P_lm,Q_l)最大的广义特征值对应的特征向量，作为第l个子矩阵的第m列，并执行步骤(5j)；

(5j)将l的值增加1，判断此时l的值是否小于子矩阵的总数目r：若是，则返回步骤(5c)，否则，执行步骤(5k)；

(5k)将m的值增加1，判断此时m的值是否小于每个子矩阵的总列数n：若是，则返回步骤(5b)，否则，执行步骤(6)；

(6)判断当前块对角化矩阵B的估计误差与前一次的块对角化矩阵B的估计误差的差值绝对值是否大于设定迭代终止阈值若是，则返回步骤(5)，否则，将得到的块对角化矩阵B输出，并执行步骤(7)；

(7)使用块对角化矩阵B从观测信号中分离出源信号。

2.根据权利要求1所述的基于多目标优化联合块对角化的卷积盲信号分离方法，其特征在于，所述步骤(2)中的二阶时延相关矩阵R(τ_k)，按照下式计算：

其中，X_i(t)为观测数据，X_i(t+τ_k)为第k个时延τ_k的观测数据，矢量G表示观测信号样本数据段的总数目，τ_k表示设定的第k个时延。

3.根据权利要求1所述的基于多目标优化联合块对角化的卷积盲信号分离方法，其特征在于，所述步骤(5e)中的海森矩阵Q_l的条件数ξ，按照下式计算：

ξ＝||Q_l||·||Q_l ^-1||

其中，[·]^-1表示矩阵的求逆操作，||·||表示求范数操作。

4.根据权利要求1所述的基于多目标优化联合块对角化的卷积盲信号分离方法，其特征在于，所述步骤(5i)中计算矩阵对(P_lm,Q_l)最大的广义特征值对应的特征向量，按如下步骤计算：

第一步，按照下式，计算矩阵对(P_lm,Q_l)的广义特征值和对应的特征向量：

[V D]＝eig(P_lm,Q_l)

其中，D表示矩阵对(P_lm,Q_l)的广义特征值构成的对角矩阵，V表示矩阵对(P_lm,Q_l)的广义特征值对应的特征向量构成的矩阵，eig表示求广义特征值操作；

第二步，将所述V的第一列作为矩阵对(P_lm,Q_l)最大的广义特征值对应的特征向量。

5.根据权利要求1所述的基于多目标优化联合块对角化的卷积盲信号分离方法，其特征在于，所述步骤(7)中的使用块对角化矩阵B从观测信号中分离出源信号，其公式表示如下：

其中，y(t)表示t时刻的源信号，B表示块对角化矩阵，X_i(t)表示观测数据矢量，表示求矩阵的伪逆操作。