CN108282424B - 用于四个数据集联合盲源分离的四阶张量联合对角化算法 - Google Patents

Abstract

本发明公开了一种用于四个数据集联合盲源分离的四阶张量联合对角化算法，包括以下步骤：S1、观测信号：对四个数据集观测信号分别进行预白化；S2、目标张量：预白化之后，构造一组互四阶累积量张量；S3、初始化因子矩阵；S4、代价函数收敛计算：若一次扫描过后，算法收敛，则计算结束，若算法仍未收敛，则以此次更新所得的因子矩阵作为初始值，进行下一次扫描，遍历更新雅克比旋转矩阵，更新因子矩阵，直至收敛为止。本发明所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，该算法基于正交旋转变换，因而能够获得最小二乘意义上的最优解，是一种针对不多于四个数据集信号的J‑BSS方法。

Description

用于四个数据集联合盲源分离的四阶张量联合对角化算法

技术领域

本发明涉及生物医学、通信、语音、阵列等信号处理领域，一种可用于四个数据集联合盲源分离的四阶张量联合对角化算法。

背景技术

联合盲源分离(Joint Blind Source Separation，J-BSS)作为一种新兴的、基于数据驱动的多数据集融合处理技术，在生物医学、通信、语音、阵列等信号处理领域获得了大量关注[1]–[14]。

同传统的盲源分离(BSS)不同[14]–[24]，J-BSS通常假设如下多数据集信号模型：

x^(m)(t)＝A^(m)s^(m)(t),m＝1,2,...,M, (1A)

其中，x(t)^(m)∈C^N和s(t)^(m)∈C^R分别表示第m个数据集的观测信号和源信号， A^(m)∈C^N×R表示第m个数据集的混合矩阵。

参数M、N、R分别表示数据集个数、观测信号通道数及源信号个数。

J-BSS旨在仅已知多数据集观测信号x^(m)(t)的前提下，对多数据集混合矩阵 A^(m)及源信号s^(m)(t)进行联合辨识。

在上述模型中，多数据集信号往往是同一目标或物理现象，经由不同的观测方式所获得的多组数据。

在这里，“不同的观测方式”是不同的设备或实验方案，如脑电(EEG)与功能磁共振(fMRI)[2]；或是同一种设备和方案下，不同的被试个体，如多被试fMRI[3],[7]；或者是相同的设备、方案及被试个体，在不同变换域的观测 (如频域、时域、空域、统计域等[3],[10])。

在处理上述多数据集信号时，由于J-BSS利用了不同集合间的相关性及相异性，经融合处理后的信息能够更为全面地描述观测对象在多层次、多剖面、多模态意义上的特征，因此在辨识能力、分离精度等方面之BSS具有性能优势。

关于J-BSS的研究最早可追溯至上世纪30年代，文献[1]所提出的典范相关分析方法(Canonical Correlation Analysis:CCA)。在后来的较长一段时间内， CCA被广泛用于两个数据集信号的融合预处理之中。

近些年来，CCA被拓展至多于两个数据集的情形，即多集合CCA(Multiset CCA:M-CCA)[6]，并在多被试fMRI数据融合中得以应用。随后，M-CCA 的耦合矩阵分解模型进一步演变为矩阵广义联合对角化(Generalized joint diagonalization:GJD)模型[8]。

此外，由经典的独立成分分析(Independent Component Analysis:ICA)发展而来的独立矢量分析(Independent Vector Analysis:IVA)也是J-BSS的主要方法之一[26]–[28]。

最近一两年，面向多集合、多维源信号分离问题的联合独立子空间分析 (jointindependent subspace analysis:J-ISA)逐步被人们所关注[29]。

在上述研究中，GJD是J-BSS的一类重要方法。这类方法通过计算互累积量，将多数据集信号转化为一组相互耦合，且具有可联合对角化结构的矩阵切片，从而通过对上述结构进行代数拟合实现J-BSS。

目前，有关GJD的研究工作大多基于信号二阶互协方差。而仅有的基于四阶互累积量的方法[8]，也只是考虑了四阶互累积量在矩阵层面的低阶可联合对角化结构，未能充分挖掘其高阶可联合对角化结构。此外，该方法并未对信号的时域特性(如非平稳性)加以利用，进而导致算法的精准度不足。

发明内容

根据上述提出的技术问题，而提供一种用于四个数据集联合盲源分离的四阶张量联合对角化算法，用以弥补现有技术未利用信号联合高阶统计特性的缺陷，进而导致算法的精准度不足的缺点。本发明采用的技术手段如下：

一种用于四个数据集联合盲源分离的四阶张量联合对角化算法，包括以下步骤：

S1、观测信号：

对四个数据集观测信号分别进行预白化：

记第m个数据集的白化矩阵为W^(m)∈C^R×N,m＝1,…,4；

则预白化处理之后，第m个数据集观测信号为：

y^(m)(t)＝W^(m)x^(m)(t)＝U^(m)s^(m)(t)∈C^R, (1)

其中，酉矩阵U^(m)＠W^(m)A^(m)∈C^R×R是经预白化之后，第m个数据集信号的混合矩阵；

S2、目标张量：

预白化之后，构造一组互四阶累积量张量如下：

T_k＝＝cum(y⁽¹⁾(t_k),y⁽²⁾(t_k),y⁽³⁾(t_k),y⁽⁴⁾(t_k)),k＝1,...,K, (2)

其中，“cum(·)”用于计算互四阶累积量，其定义为：

其中，“E(·)”用于计算数学期望，上标“*”表示共轭，根据定义可知，每一个张量T_k∈C^R×R×R×R表示四组观测信号在时刻t_k的互四阶累积量张量，k＝ 1,…,K；

其公式表示如下：

T_k＝D_k×₁U⁽¹⁾×₂U^(2)*×₃U^(3)*×₄U⁽⁴⁾, (4)

其中×_n表示张量和矩阵的n模乘积，n＝1,…,4；一个四阶张量

与矩阵

的n模乘积定义为：

公式(4)中D_k＝cum(s⁽¹⁾(t_k),s⁽²⁾(t_k),s⁽³⁾(t_k),s⁽⁴⁾(t_k))∈C^R×R×R×R为四组数据集的源信号在时刻t_k的互四阶累积量张量；

S3、初始化因子矩阵；

S4、代价函数收敛计算：

首先，优化准则：

定义目标函数λ为数据张量T_k与拟合张量D_k×₁U⁽¹⁾×₂U^(2)*×₃U^(3)*×₄U⁽⁴⁾间的均方误差：

其中||·||_F表示弗罗贝尼乌斯范数，通过最小化λ，计算混合矩阵在最小二乘意义下的最优估计值；

注意到U^(m)为酉矩阵，m＝1,2,3,4，根据酉变换的保范性，公式(7)改写为：

其中，off(·)用于将其输入张量的超对角元置0，diag(·)用于对其输入张量的非对角元置0，

为常数。因此，对λ进行最小化等价于对

进行最大化，因此本发明通过如下准则估计混合矩阵U^(m),m＝1,2,3,4：

其中

分别表示U⁽¹⁾,U⁽²⁾,U⁽³⁾,U⁽⁴⁾的估计值；

公式(7)–(9)表明，通过最大化T_k×₁U^(1)H×₂U^(2)T×₃U^(3)T×₄U^(4)H的超对角元范数平方和，即对T₁,...,T_K进行联合对角化，获得混合矩阵U⁽¹⁾,U⁽²⁾,U⁽³⁾,U⁽⁴⁾在最小二乘意义下的最优估计值；

其次，雅克比迭代：

雅克比迭代将待求解的酉矩阵写成一系列雅克比旋转矩阵的乘积，进而通过分别优化每一个雅克比旋转矩阵，最终实现目标函数的最大化；

具体而言，记数据张量T_k及混合矩阵U^(m)在前一次迭代的更新值分别为

和

其在当前迭代的更新值分别为

和

每一次迭代通过雅克比旋转矩阵

对T_k和U^(m)进行更新，即：

其中雅克比旋转矩阵

定义如下：

其中

‘i’表示虚部单位；

根据定义，

除了(i,i),(i,j),(j,i),(j,j)四个位置的元素不为0，对角线上元素为1之外，其余位置皆为0；

令坐标i的值从1取至R，j的值由i取至R，对于某一对固定的坐标值(i,j)，由(9)和(10)可知，

的最优解通过求解下述优化问题获得，m＝1,2,3,4：

获得

之后，按照式(10)对T_k和U^(m)进行更新；坐标(i,j)遍历所有取值时所包含的全部迭代称为一次扫描；

若一次扫描过后，算法收敛，则计算结束，若算法仍未收敛，则以此次更新所得的因子矩阵作为初始值，进行下一次扫描，遍历更新雅克比旋转矩阵，更新因子矩阵，直至收敛为止。

作为优选在步骤S1之前，做如下假设：

(A1)组内独立性：当1≤r≠u≤R时，

和

统计独立；

(A2)组间相关性：当1≤r＝u≤R时，

和

统计相关；

(A3)源信号为非高斯、非平稳信号；

(A4)所考虑的多数据集模型为正定的，即N≥R；

(A5)数据集的数目M＝4。

作为优选在实际中，在假设条件(A1)和(A4)下，白化矩阵W^(m)∈C^R×N可通过观测信号x^(m)(t)的二阶协方差矩阵的奇异值分解而获得。

作为优选由假设(A1)，(A2)及(A3)可知，每一个张量D_k均为超对角张量，即：

定义U′^(m)@U^(m)D^(m)P,，m＝1,2,3,4，其中D⁽¹⁾,D⁽²⁾,D⁽³⁾,D⁽⁴⁾∈C^R×R为对角矩阵，且满足D⁽¹⁾D^(2)*D^(3)*D⁽⁴⁾＝I，I为单位矩阵，P∈C^R×R为排序矩阵；

不难得知，若将公式(4)中的U^(m)替换为U′^(m)，等式依然成立；

因此，U′^(m)在不考虑幅度/相位模糊以及顺序模糊时，和混合矩阵U^(m)等价；

通过对方程(4)进行求解，可对混合矩阵U^(m)进行估计。

作为优选遍历更新雅克比旋转矩阵中，每一步迭代中所涉及

的闭式最优解；

雅克比旋转矩阵的闭式最优解

根据雅克比旋转矩阵的性质，酉变换

仅改变张量

在如下八部分的取值：

这里我们采用MATLAB符号来表示

的子张量，例如

表示固定张量

的第一个索引值为i，其余索引值不固定时所获得的子张量；

能够看出，最大化张量

超对角元素的范数平方和，等价于最大化上述八个子张量相交所构成的大小为2×2×2×2的子张量的超对角元素范数平方和；

因此，优化准则(12)可改写为：

为了进一步简化优化过程，我们采用交替更新的方式，即用如下四步来代替(13)，其中每一步用于交替更新

其中

且

接下来，我们以

为例，解释(14)的求解步骤；

首先根据定义可知：

其中，

根据公式(15)，有：

其中M_i,j＝[m_1,i,j,...,m_K,i,j,q_1,i,j,...,q_K,i,j]∈C^3×2K；

上述推导表明，对公式(16)进行最大化时，

对应于矩阵

的主特征向量，记

的估计值为

则有

进而可由(11) 获得

接下来，进行更新

并按照同求解

类似的步骤获得

同理，计算获得

和

作为优选迭代停止条件的判定为：

计算相邻两次扫描所获得张量的对角元范数和的相对变化值：

当其小于阈值ε时，则认为算法收敛，迭代终止；

或，判断当前扫描所获得的雅克比旋转矩阵是否接近于单位矩阵，具体而言，定义

当ξ⁽¹⁾、ξ⁽²⁾、ξ⁽³⁾和ξ⁽⁴⁾同时小于阈值ε时，则说明

和

已经近似于单位阵，不能更新目标张量，可认为算法收敛，迭代终止。

与现有技术相比较，本发明所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，是一种针对不多于四个数据集信号的J-BSS方法。该方法对多数据集信号计算其在不同时间段的互四阶累积量张量，并通过对其进行张量联合对角化，实现多数据集信号的J-BSS。

此外，我们提出了一种基于Givens连续旋转的张量联合对角化算法。同现有技术[30]中提出的仅适用于单数据集BSS的张量对角化算法相比，该算法适用于张量的加载因子矩阵不相同的情况，从而适用于J-BSS。同现有技术[8]中提出的非正交张量对角化算法相比，该算法基于正交旋转变换，因而能够获得最小二乘意义上的最优解。

本发明所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，应用了信号的时域特性，本发明对多数据集信号计算其在不同时间段的互四阶累积量张量，提高了算法的精准度。

本发明所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，利用四个因子矩阵构建了张量，在算法的过程中，同时交替更新四个因子矩阵，算法精准，同现有技术中的算法只是交替更新前三个因子矩阵相比，第四个因子矩阵在得出前三个因子矩阵后，自动求出。

本发明的优化准则直接对四个因子矩阵进行优化，现有技术的算法只优化前三个因子矩阵。由于优化准则不同，且被优化的张量也不同，所以后续的优化过程都是不同的，现有技术中由于前三个因子矩阵已被求出，故而第四个因子矩阵通过还原目标张量而获得，但是这样获得的因子矩阵精度低。本发明直接优化四个因子矩阵，提高了第四个因子矩阵的精度。

本发明在计算张量时，将信号分为k段，共计算出k个张量，同现有技术中算法信号不分段，计算出n个张量，也就是同因子矩阵的维度个张量相比，利用了更多的信号信息，从而提高算法精度。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明算法流程图。

图2-1是本发明实施例1中无噪声条件下JTD和GOJD算法的A-PI值随扫描次数的变化曲线(K＝3，N＝3)。

图2-2是本发明实施例1中无噪声条件下JTD和GOJD算法的A-PI值随扫描次数的变化曲线(K＝3，N＝50)。

图2-3是本发明实施例1中无噪声条件下JTD和GOJD算法的A-PI值随扫描次数的变化曲线(K＝50，N＝3)。

图2-4是本发明实施例1中无噪声条件下JTD和GOJD算法的A-PI值随扫描次数的变化曲线(K＝50，N＝50)。

图3-1是本发明实施例2中A-PI随SNR变化曲线示意图(快拍数＝100000)。

图3-2是本发明实施例2中A-PI随快拍数变化曲线示意图(SNR＝5dB)。

图4本发明实施例3中八路ECG数据图。

图5是本发明实施例3中实际FECG信号分离结果图。

具体实施方式

如图1所示，一种用于四个数据集联合盲源分离的四阶张量联合对角化算法，首先吗，做如下假设：

(A1)组内独立性：当1≤r≠u≤R时，

和

统计独立；

(A2)组间相关性：当1≤r＝u≤R时，

和

统计相关；

(A3)源信号为非高斯、非平稳信号；

(A4)所考虑的多数据集模型为正定的，即N≥R；

(A5)数据集的数目M＝4。

在上述假设中，(A1)和(A2)是现有的联合盲源分离方法中普遍用到的假设条件；条件(A3)和(A4)则表明本发明所述方法适用于正定情况下，非高斯及非平稳信号的联合盲源分离；条件(A5)则限定本发明方法所能处理的数据集数目不超过四个，本发明以四个数据集为例阐述所提出的联合盲源分离算法，低于四个数据集的情况通过类似的推导获得。

具体包括以下步骤：

S1、观测信号：

对四个数据集观测信号分别进行预白化：

记第m个数据集的白化矩阵为W^(m)∈C^R×N,m＝1,…,4；

则预白化处理之后，第m个数据集观测信号为：

y^(m)(t)＝W^(m)x^(m)(t)＝U^(m)s^(m)(t)∈C^R, (1)

其中，酉矩阵U^(m)＠W^(m)A^(m)∈C^R×R是经预白化之后，第m个数据集信号的混合矩阵；在实际中，在假设条件(A1)和(A4)下，白化矩阵W^(m)∈C^R×N可通过观测信号x^(m)(t)的二阶协方差矩阵的奇异值分解而获得[8]。

S2、目标张量：

预白化之后，构造一组互四阶累积量张量如下：

T_k＝＝cum(y⁽¹⁾(t_k),y⁽²⁾(t_k),y⁽³⁾(t_k),y⁽⁴⁾(t_k)),k＝1,...,K, (2)

其中，“cum(·)”用于计算互四阶累积量，其定义为：

其中，“E(·)”用于计算数学期望，上标“*”表示共轭，根据定义可知，每一个张量T_k∈C^R×R×R×R表示四组观测信号在时刻t_k的互四阶累积量张量，k＝ 1,…,K；

其公式表示如下：

T_k＝D_k×₁U⁽¹⁾×₂U^(2)*×₃U^(3)*×₄U⁽⁴⁾, (4)

其中×_n表示张量和矩阵的n模乘积，n＝1,…,4；一个四阶张量

与矩阵

的n模乘积定义为：

公式(4)中D_k＝cum(s⁽¹⁾(t_k),s⁽²⁾(t_k),s⁽³⁾(t_k),s⁽⁴⁾(t_k))∈C^R×R×R×R为四组数据集的源信号在时刻t_k的互四阶累积量张量；

由假设(A1)，(A2)及(A3)可知，每一个张量D_k均为超对角张量，即：

定义U′^(m)@U^(m)D^(m)P,，m＝1,2,3,4，其中D⁽¹⁾,D⁽²⁾,D⁽³⁾,D⁽⁴⁾∈C^R×R为对角矩阵，且满足D⁽¹⁾D^(2)*D^(3)*D⁽⁴⁾＝I，I为单位矩阵，P∈C^R×R为排序矩阵；

不难得知，若将公式(4)中的U^(m)替换为U′^(m)，等式依然成立；

因此，U′^(m)在不考虑幅度/相位模糊(由对角矩阵D^(m)表征)以及顺序模糊 (由排序矩阵P表征)时，和混合矩阵U^(m)等价。通过对方程(4)进行求解，可对混合矩阵U^(m)进行估计。

联合盲源分离的目的，正是在不考虑上述两类平凡模糊的前提条件下，对混合矩阵进行辨识，进而对源信号进行估计。

同传统的BSS有所不同的是，当m取不同索引值时，不同的U′^(m)相对于真实混合矩阵U^(m)的排序矩阵P是完全相同的，这意味着经由J-BSS所估计的混合矩阵及源信号成分，在顺序上是自然对齐的，这也是J-BSS相较于BSS的优势之一。

本发明中利用四个因子矩阵构建了张量，在算法的过程中，同时交替更新四个因子矩阵，同现有技术中的算法只是交替更新前三个因子矩阵相比，第四个因子矩阵在得出前三个因子矩阵后，自动求出。

本发明在计算张量时，将信号分为k段，共计算出k个张量，同现有技术中算法信号不分段，计算出n个张量，也就是因子矩阵的维度个张量相比，利用了更多的信号信息，从而提高算法精度。

S3、初始化因子矩阵；

S4、代价函数收敛计算：

首先，优化准则：

定义目标函数λ为数据张量T_k与拟合张量D_k×₁U⁽¹⁾×₂U^(2)*×₃U^(3)*×₄U⁽⁴⁾间的均方误差：

其中||·||_F表示弗罗贝尼乌斯范数(Frobenius norm)，通过最小化λ，计算混合矩阵在最小二乘意义下的最优估计值；

注意到U^(m)为酉矩阵，m＝1,2,3,4，根据酉变换的保范性，公式(7)改写为：

其中，off(·)用于将其输入张量的超对角元置0，diag(·)用于对其输入张量的非对角元置0，

为常数。因此，对λ进行最小化等价于对

进行最大化，因此本发明通过如下准则估计混合矩阵U^(m),m＝1,2,3,4：

其中

分别表示U⁽¹⁾,U⁽²⁾,U⁽³⁾,U⁽⁴⁾的估计值；

公式(7)–(9)表明，通过最大化T_k×₁U^(1)H×₂U^(2)T×₃U^(3)T×₄U^(4)H的超对角元范数平方和，即对T₁,...,T_K进行联合对角化，获得混合矩阵U⁽¹⁾,U⁽²⁾,U⁽³⁾,U⁽⁴⁾在最小二乘意义下的最优估计值；接下来，我们提出一种基于雅克比迭代的张量联合对角化算法。

本发明的优化准则直接对四个因子矩阵进行优化，现有技术中算法只优化前三个因子矩阵。由于优化准则不同，且被优化的张量也不同，所以后续的优化过程都是不同的，现有技术中由于前三个因子矩阵已被求出，故而第四个因子矩阵通过还原目标张量而获得，但是这样获得的因子矩阵精度低。本发明直接优化四个因子矩阵，提高了第四个因子矩阵的精度。

其次，雅克比迭代：

雅克比迭代将待求解的酉矩阵写成一系列雅克比旋转矩阵的乘积，进而通过分别优化每一个雅克比旋转矩阵，最终实现目标函数的最大化；

具体而言，记数据张量T_k及混合矩阵U^(m)在前一次迭代的更新值分别为

和

其在当前迭代的更新值分别为

和

每一次迭代通过雅克比旋转矩阵

(又称为Givens旋转矩阵)对T_k和U^(m)进行更新，即：

其中雅克比旋转矩阵

定义如下：

其中

‘i’表示虚部单位。根据定义，

除了(i,i),(i,j),(j,i),(j,j)四个位置的元素不为0，对角线上元素为1之外，其余位置皆为0；

令坐标i的值从1取至R，j的值由i取至R，对于某一对固定的坐标值(i,j)，由(9)和(10)可知，

的最优解通过求解下述优化问题获得，m＝1,2,3,4：

获得

之后，按照式(10)对T_k和U^(m)进行更新；坐标(i,j)遍历所有取值时所包含的全部迭代称为一次扫描(Sweep)；

若一次扫描过后，算法收敛，则计算结束，若算法仍未收敛，则以此次更新所得的因子矩阵作为初始值，进行下一次扫描，遍历更新雅克比旋转矩阵，更新因子矩阵，直至收敛为止。

迭代停止条件的判定为：

计算相邻两次扫描所获得张量的对角元范数和的相对变化值：

当其小于阈值ε时，则认为算法收敛，迭代终止；

或，判断当前扫描所获得的雅克比旋转矩阵是否接近于单位矩阵，具体而言，定义

当ξ⁽¹⁾、ξ⁽²⁾、ξ⁽³⁾和ξ⁽⁴⁾同时小于阈值ε时，则说明

和

已经近似于单位阵，不能更新目标张量，可认为算法收敛，迭代终止。

从上述描述可知，雅克比迭代将优化问题(9)，化解为一系列子优化问题 (12)。由于每一个雅克比旋转矩阵仅有两个未知参数，其最优解有望具有简单的闭式解形式。接下来，我们将推导每一步迭代中所涉及

的闭式最优解。

遍历更新雅克比旋转矩阵中，每一步迭代中所涉及

的闭式最优解；

雅克比旋转矩阵的闭式最优解

根据雅克比旋转矩阵的性质，酉变换

仅改变张量

在如下八部分的取值:

这里我们采用MATLAB符号来表示

的子张量，例如

表示固定张量

的第一个索引值为i，其余索引值不固定时所获得的子张量；

不难看出，最大化张量

超对角元素的范数平方和，等价于最大化上述八个子张量相交所构成的大小为2×2×2×2的子张量的超对角元素范数平方和；

因此，优化准则(12)可改写为：

为了进一步简化优化过程，我们采用交替更新的方式，即用如下四步来代替(13)，其中每一步用于交替更新

其中

且

接下来，我们以

为例，解释(14)的求解步骤；

首先根据定义可知：

其中

根据公式(15)，有：

其中M_i,j＝[m_1,i,j,...,m_K,i,j,q_1,i,j,...,q_K,i,j]∈C^3×2K；

上述推导表明，对公式(16)进行最大化时，

对应于矩阵

的主特征向量，记

的估计值为

则有

进而可由(11) 获得

接下来，进行更新

并按照同求解

类似的步骤获得

同理，计算获得

和

上述推导假设J-BSS问题为复值，当J-BSS问题为实值时，亦通过类似的计算将其转化为实值四阶张量的联合对角化问题，并通过实值雅克比旋转对之进行求解，计算步骤与复值情况类似。

需要注意的是，此时实值雅克比旋转矩阵定义为：

且公式(15)中

m_k,i,j,q_k,i,j的表达式由下式给出：

实施例1，随机生成超对角张量D_k，和酉矩阵U⁽¹⁾,U⁽²⁾,U⁽³⁾,U⁽⁴⁾,并按照公式 (4)直接构造张量

算法性能通过如下式定义的性能指标(Performance Index:PI)进行评价：

其中，

上标

表示Moore-Penrose伪逆。对分离结果分别计算其PI之后，再进行平均。以此平均PI(Average PI:A-PI)作为算法性能的评价指标。

通过计算每一次扫描结束时算法所对应的PI值，可反映算法的收敛性。图 2-1至图2-4绘制了在K和N取不同值时，JTD算法在10次独立实验中，其 A-PI值随扫描次数的变化曲线。

作为参考，我们绘制在相同条件下，广义正交矩阵对角化算法(GOJD，[8]) 的A-PI曲线(其数据矩阵由T_k＝固定后三维索引获得)。看到，所提出的JTD 算法呈现出单调下降、线性收敛的收敛模式。同GOJD算法相比，JTD算法平均5次扫描即能达到收敛，具有更优的收敛性能。

图2-1至图2-4，为无噪声条件下JTD和GOJD算法的A-PI值随扫描次数的变化曲线。这里展示四种条件下的实验结果：(a)如图2-1所示，张量的尺寸N和数量K均取较小值；(b)如图2-2所示，张量的尺寸N取较大值，数量 K取较小值；(c)如图2-3所示，张量的尺寸N取较小值，数量K取较大值； (d)如图2-4所示，张量的尺寸N和数量K均取较大值。

实施例2，本实验中，根据公式(1A)产生四组观测信号x^(m)(t),m＝1,2,3,4：

其中混合矩阵A^(m)∈C^N×R由程序随机产生，噪声项n^(m)(t)∈C^N×T为零均值单位方差的高斯白噪声。

源信号由下式构造：

s_r(t)＝Q_rs′_r(t), (22)

其中

s′_r(t)∈C⁴为分段调幅的随机相位信号，其相位及调幅系数由程序随机产生。矩阵Q_r∈R^4×4由程序随机产生，用于引入集合间的相关性。不难得知，所构造的源信号为非平稳非高斯信号(假设 A3)，且不同数据集的源信号满足组间相关和组内独立的假设条件(假设A1， A2)。

在公式(21)中，σ_s和σ_n分别代表信号强度和噪声强度。进一步可定义信噪比(SNR)如下：SNR@20lgσ_s/σ_n。

我们用所提出的基于JTD的J-BSS算法对上述多数据集信号进行分离。并将其同下述同类型的BSS及J-BSS算法进行对比：

(a)特征矩阵近似联合对角化算法(Joint Approximate Diagonalization ofEigenmatrices:JADE)[16]；

(b)三阶张量联合对角化算法(Simultaneous Third-Order TensorDiaognalization:STOTD)[30]；

(c)多集合典范相关分析算法(Multiset Canonical Correlation Analysis:MCCA)[6]；

(d)广义正交矩阵联合对角化算法(Generalized Orthogonal JointDiagonalization:GOJD)[8]。

其中JADE和STOTD为基于四阶累积量的BSS算法，MCCA和GOJD为 J-BSS算法。在GOJD算法的实验中，我们分别采用了基于四阶累积量和二阶协方差的算法，分别记为GOJD4和GOJD2。

作我们用采样互四阶累积量为(2)所定义的互四阶累积量的估计值(2)，其计算公式如下：

其中

表示预白化以后观测信号在时间维度上的第k个分段，m＝ 1,…,4。这里分段长度为L，相邻段之间的重叠率α∈[0,1)，定义为相邻段的重合点数与L的比值。

在本实验中，我们固定各数据集的观测通道数N＝10，源信号数为R＝5， A-PI曲线中的每个点由200次蒙特卡洛独立实验的结果平均所获得。首先，固定时间采样点数T＝100000，L＝T/10，α＝0.3。

令SNR从0dB变化至10dB。接下来，固定SNR＝5dB，令时间采样点数由10000变化至100000。在上述两种情况下，所比较算法的A-PI曲线如图3-1 和图3-2所示。看到，所提出的JTD算法性能明显优于除GOJD2之外的其它算法。

同GOJD2相比，由图3-1看出，在SNR为2–8dB时，JTD算法的分离精度更高。图3-2则表明，当快拍数较高时(80000–100000)所提出的JTD算法精度更高，反之则是GOJD2更具优势。

这主要是因为GOJD2所采用的二阶协方差较之JTD所采用的的四阶累积量，在短快拍的情况下，具有更小的有限采样误差。

如图3-1和图3-2所示，JTD，GJD，MCCA，JADE，STOTD分离计算机合成信号时的性能对比。

实施例3，实验旨在从实际采集的，混合了孕妇心电信号(Mother ECG, MECG)，胎儿心电信号(Fetal ECG，FECG)，及其它干扰信号的观测数据中，提取FECG信号。

采用的八路ECG数据取自鲁汶大学的DAISY数据库0，信号波形如图4 所示。采样频率为250Hz，采样点数为2500，采样时间是10s，如图4所示，是实际采集的八路ECG数据，数据取自鲁汶大学的DAISY数据库。

按照下式，将八路观测信号划分为四个数据集，每个数据集中含有五路观测信号：

将四个数据集的观测信号预白化后分段，每段长度为L＝500，重叠率为α＝0.5，不难得知此时K＝10。对每一段信号计算互四阶累积量以构建K个目标张量。用JTD算法对其进行联合对角化，获得解混矩阵后，由公式(1)可获得各源信号成分，结果如图4所示，

表示第m个数据集的源信号的估计值，如图5所示，能够看到MECG信号和FECG信号被成功分离。

从图中看出，每个数据集第1路信号和第2路信号为MECG，前三个数据集中第3路信号为FECG。结果表明，基于JTD的J-BSS算法能够成功分离实际采集的FECG信号。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

参考文献

[1]Hotelling H.Relations between two sets of variates[J].Biometrika,1936,28(3/4):321-377.

[2]Sui J,Adali T,Yu Q,et al.A review of multivariate methods formultimodal fusion of brain imaging data[J].Journal of Neuroscience Methods,2012,204(1):68-81.

[3]Kim T,Attias H T,Lee S Y,et al.Blind source separation exploitinghigher-order frequency dependencies[J].IEEE Trans on Audio,Speech andLanguage Processing,2007,15(1): 70-79.

[4]Lee J H,Lee T W,Jolesz F A,et al.Independent vector analysis(IVA):multivariate approach for fMRI group study[J].Neuroimage,2008,40(1):86-109.

[5]Anderson M,Adali T,Li X L.Joint blind source separation withmultivariate Gaussian model:Algorithms and performance analysis[J].IEEE Transon Signal Processing,2012, 60(4):1672-1683.

[6]Li Y O,Adali T,Wang W,et al.Joint blind source separation bymultiset canonical correlation analysis[J].IEEE Trans on Signal Processing,2009,57(10):3918-3929.

[7]Correa N M,Adali T,Li Y O,et al.Canonical correlation analysis fordata fusion and group inferences[J].IEEE Signal Processing Magazine,2010,27(4):39-50.

[8]Li X L,Adal₁ T,Anderson M.Joint blind source separation bygeneralized joint diagonalization of cumulant matrices[J].Signal Processing,2011,91(10):2314-2322.

[9]Congedo M,Phlypo R,Chatel-Goldman J.Orthogonal and non-orthogonaljoint blind source separation in the least-squares sense[A].Signal ProcessingConference(EUSIPCO)[C]. Bucharest,Romania,2012:1885-1889.

[10]Gong X F,Wang X L,Lin Q H.Generalized Non-Orthogonal JointDiagonalization With LU Decomposition and Successive Rotations[J].IEEE Transon Signal Processing,2015, 63(5):1322-1334.

[11]Lahat D,Jutten C.Joint independent subspace analysis:a quasi-Newton algorithm[A]. International Conference on Latent Variable Analysis andSignal Separation[C].Liberec, Czech Republic,2015:111-118.

[12]Lahat D,Adali T,Jutten C.Multimodal data fusion:an overview ofmethods,challenges,and prospects[J].Proceedings of the IEEE,2015,103(9):1449-1477.

[13]Adali T,Levin-Schwartz Y,Calhoun V D.Multimodal data fusion usingsource separation: Two effective models based on ICA and IVA and theirproperties[J].Proceedings of the IEEE, 2015,103(9):1478-1493.

[14]Hérault J,Jutten C,Ans B.Détection de grandeurs primitives dansun message composite par une architecture de calcul neuromimétique enapprentissage non supervisé[A]. Proceedings of the 10th Workshop Traitementdu signal et ses applications(GRETSI)[C]. Nice,France,1985:1017–1022.

[15]Comon P.Independent component analysis,a new concept？[J].SignalProcessing,1994, 36(3):287-314.

[16]Cardoso J F,Souloumiac A.Blind beamforming for non-Gaussiansignals[A].IEE Proceedings F(Radar and Signal Processing)[C].1993,140(6):362-370.

[17]Oja E,Hyvarinen A.A fast fixed-point algorithm for independentcomponent analysis[J]. Neural Computation,1997,9(7):1483-1492.

[18]Belouchrani A,Abed-Meraim K,Cardoso J F,et al.A blind sourceseparation technique using second-order statistics[J].IEEE Trans on SignalProcessing,1997,45(2):434-444.

[19]Novey M,Adali T.Complex ICA by negentropy maximization[J].IEEETrans on Neural Networks,2008,19(4):596-609.

[20]Sawada H,Araki S,Makino S.Underdetermined convolutive blindsource separation via frequency bin-wise clustering and permutation alignment[J].IEEE Trans on Audio,Speech and Language Processing,2011,19(3):516-527.

[21]Yu Y,Petropulu A P.PARAFAC-based blind estimation of possiblyunderdetermined convolutive MIMO systems[J].IEEE Trans on Signal Processing,2008,56(1):111-124.

[22]Lin Q H,Liu J,Zheng Y R,et al.Semiblind spatial ICA of fMRI usingspatial constraints[J].Human Brain Mapping,2010,31(7):1076-1088.

[23]Cichocki A,Zdunek R,Phan A H,et al.Nonnegative matrix and tensorfactorizations: applications to exploratory multi-way data analysis and blindsource separation[M].John Wiley&Sons,2009.

[24]Handbook of Blind Source Separation:Independent componentanalysis and applications[M].Academic Press,2010.

[25]Li Y O,Wang W,Adali T,et al.CCA for joint blind source separationof multiple datasets with application to group fMRI analysis[A].IEEEInternational Conference on Acoustics, Speech and Signal Processing[C].2008:1837-1840.

[26]Kim T,Eltoft T,Lee T W.Independent vector analysis:An extensionof ICA to multivariate components[A].International Conference on IndependentComponent Analysis and Signal Separation[C].Springer Berlin Heidelberg,2006:165-172.

[27]Anderson M,Li X L,Adal₁ T.Nonorthogonal independent vectoranalysis using multivariate Gaussian model[A].International Conference onLatent Variable Analysis and Signal Separation[C].Springer Berlin Heidelberg,2010:354-361.

[28]Anderson M,Adali T,Li X L.Joint blind source separation withmultivariate Gaussian model:Algorithms and performance analysis[J].IEEE Transon Signal Processing,2012, 60(4):1672-1683.

[29]Lahat D,Jutten C.Joint independent subspace analysis:a quasi-Newton algorithm[C]. International Conference on Latent Variable Analysis andSignal Separation[C].Liberec, Czech Republic,2015:111-118.

[30]De Lathauwer L,De Moor B,Vandewalle J.Independent componentanalysis and (simultaneous)third-order tensor diagonalization[J].IEEE Transon Signal Processing,2001, 49(10):2262-2271.

Moor D.Database for the identification of systems(DaISy)[Z].http://www.esat.kuleuven.ac. be/sista/daisy,2010。

Claims (6)

1.一种用于四个数据集联合盲源分离的四阶张量联合对角化算法，其特征在于包括以下步骤：

S1、观测信号：

对四个数据集观测信号分别进行预白化：

记第m个数据集的白化矩阵为

则预白化处理之后，第m个数据集观测信号为：

其中，酉矩阵U^(m)@W^(m)A^(m)∈C^R×R是经预白化之后，第m个数据集信号的混合矩阵，x^(m)(t)为第m个数据集的观测信号，s^(m)(t)为第m个数据集的源信号，A^(m)为第m个数据集的混合矩阵；

S2、目标张量：

预白化之后，构造一组互四阶累积量张量如下：

其中，“cum(·)”用于计算互四阶累积量，其定义为：

其中，“E(·)”用于计算数学期望，上标“*”表示共轭，根据定义可知，每一个张量

表示四组观测信号在时刻t_k的互四阶累积量张量，k＝1,…,K；

其公式表示如下：

其中×_n表示张量和矩阵的n模乘积，n＝1,…,4；一个四阶张量

与矩阵

的n模乘积定义为：

公式(4)中

为四组数据集的源信号在时刻t_k的互四阶累积量张量；

S3、初始化因子矩阵；

S4、代价函数收敛计算：

首先，优化准则：

定义目标函数λ为数据张量

与拟合张量

间的均方误差：

其中||·||_F表示弗罗贝尼乌斯范数，通过最小化λ，计算混合矩阵在最小二乘意义下的最优估计值；

注意到U^(m)为酉矩阵，m＝1,2,3,4，根据酉变换的保范性，公式(7)改写为：

其中，off(·)用于将其输入张量的超对角元置0，diag(·)用于对其输入张量的非对角元置0，

为常数；

因此，对λ进行最小化等价于对

进行最大化，因此，本发明通过如下准则估计混合矩阵U^(m),m＝1,2,3,4：

其中

分别表示U⁽¹⁾,U⁽²⁾,U⁽³⁾,U⁽⁴⁾的估计值；

公式(7)–(9)表明，通过最大化

的超对角元范数平方和，即对

进行联合对角化，获得混合矩阵U⁽¹⁾,U⁽²⁾,U⁽³⁾,U⁽⁴⁾在最小二乘意义下的最优估计值；

其次，雅克比迭代：

雅克比迭代将待求解的酉矩阵写成一系列雅克比旋转矩阵的乘积，进而通过分别优化每一个雅克比旋转矩阵，最终实现目标函数的最大化；

具体而言，记数据张量

及混合矩阵U^(m)在前一次迭代的更新值分别为

和

其在当前迭代的更新值分别为

和

每一次迭代通过雅克比旋转矩阵

对

和U^(m)进行更新，即：

其中雅克比旋转矩阵

定义如下：

其中

‘i’表示虚部单位；

根据定义，

除了(i,i),(i,j),(j,i),(j,j)四个位置的元素不为0，对角线上元素为1之外，其余位置皆为0；

令坐标i的值从1取至R，j的值由i取至R，对于某一对固定的坐标值(i,j)，由(9)和(10)可知，

的最优解通过求解下述优化问题获得，m＝1,2,3,4：

获得

之后，按照式(10)对

和U^(m)进行更新；坐标(i,j)遍历所有取值时所包含的全部迭代称为一次扫描；

若一次扫描过后，算法收敛，则计算结束，若算法仍未收敛，则以此次更新所得的因子矩阵作为初始值，进行下一次扫描，遍历更新雅克比旋转矩阵，更新因子矩阵，直至收敛为止。

2.根据权利要求1所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，其特征在于：

在步骤S1之前，做如下假设：

(A1)组内独立性：当1≤r≠u≤R时，

和

统计独立；

(A2)组间相关性：当1≤r＝u≤R时，

和

统计相关；

(A3)源信号为非高斯、非平稳信号；

(A4)所考虑的多数据集模型为正定的，即N≥R；

(A5)数据集的数目M＝4。

3.根据权利要求2所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，其特征在于：

在假设条件(A1)和(A4)下，白化矩阵

可通过观测信号x^(m)(t)的二阶协方差矩阵的奇异值分解而获得。

4.根据权利要求3所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，其特征在于：

由假设(A1)，(A2)及(A3)可知，每一个张量

均为超对角张量，即：

定义

其中D⁽¹⁾,D⁽²⁾,D⁽³⁾,

为对角矩阵，且满足D⁽¹⁾D^(2)*D^(3)*D⁽⁴⁾＝I，I为单位矩阵，

为排序矩阵；

若将公式(4)中的U^(m)替换为U′^(m)，等式依然成立；

因此，U′^(m)在不考虑幅度/相位模糊以及顺序模糊时，和混合矩阵U^(m)等价；

通过对方程(4)进行求解，可对混合矩阵U^(m)进行估计。

5.根据权利要求1所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，其特征在于：

遍历更新雅克比旋转矩阵中，每一步迭代中所涉及

的闭式最优解；

雅克比旋转矩阵的闭式最优解求解如下：

根据雅克比旋转矩阵的性质，酉变换

仅改变张量

在如下八部分的取值：

采用Matlab符号来表示

的子张量，

表示固定张量

的第一个索引值为i，其余索引值不固定时所获得的子张量；

最大化张量

超对角元素的范数平方和，等价于最大化上述八个子张量相交所构成的大小为2×2×2×2的子张量的超对角元素范数平方和；

因此，优化准则(12)可改写为：

为了进一步简化优化过程，我们采用交替更新的方式，即用如下四步来代替(13)，其中每一步用于交替更新

其中

且

接下来，我们以

为例，解释(14)的求解步骤；

首先根据定义可知：

其中，

根据公式(15)，有：

其中

上述推导表明，对公式(16)进行最大化时，

对应于矩阵

的主特征向量，记

的估计值为

则有

进而可由(11)获得

接下来，进行更新

并按照同求解

类似的步骤获得

同理，计算获得

和

6.根据权利要求1或5所述的用于四个数据集联合盲源分离的四阶张量联合对角化算法，其特征在于：

迭代停止条件的判定为：

计算相邻两次扫描所获得张量的对角元范数和的相对变化值：

当其小于阈值ε时，则认为算法收敛，迭代终止；

或，判断当前扫描所获得的雅克比旋转矩阵是否接近于单位矩阵，具体而言，定义

当ξ⁽¹⁾、ξ⁽²⁾、ξ⁽³⁾和ξ⁽⁴⁾同时小于阈值ε时，则说明

和

已经近似于单位阵，不能更新目标张量，算法收敛，迭代终止。

Images