CN103971124A - 一种基于相位同步的多类别运动想象脑电信号分类方法 - Google Patents

Abstract

本发明涉及一种基于相位同步的多类别运动想象脑电信号分类方法。本发明首先利用相位锁定值分别计算训练样本和测试样本的相位同步特征，然后计算两者的相关系数并对其去平均取绝对值后按从大到小的顺序排列，接着根据排序后的相关系数进行脑电信号的粗分类，再将粗分类后的脑电信号进行细分类，此过程涉及共空域模式特征提取方法和线性判别分析分类方法。本发明包括脑电信号采集、数据预处理、滤波、相位同步特征相关系数计算、特征提取和分类以及分类准确率计算。分类结果表明：采用基于相位同步的脑电信号分类方法，得到了较好的分类效果，其中基于相位同步的脑电信号的粗分类可以高效率、高准确率地确定测试样本大致所属的类别，减少计算量。

Description

一种基于相位同步的多类别运动想象脑电信号分类方法

技术领域

本发明属于生物特征识别领域中的脑电信号分类领域，具体涉及一种基于脑电信号相位同步获取相位同步特征相关系数并运用于共空域模式算法进行多类别运动想象脑电信号分类的方法。

背景技术

对于神经肌肉系统功能严重丧失的患者来说,一种新的与外界交流的手段是急需的。脑机接口(BCI)正是这样一种技术,它允许使用者直接通过意念而不是传统的神经肌肉通道来实现与外界环境的自主交流。近年来脑机接口(BCI)研究已经逐渐从两类别的模式识别发展为多类别的模式识别，从而实现多种意念与外界环境的自主交流，提升脑机接口技术的实用性。

如何有效地对多类别脑电信号进行分类，提升分类的正确率和分类的速率是目前运动想象多类别脑电信号分类研究中必须面临并解决的一个关键问题。而传统的多类别脑电信号分类方法如基于“一对一”共空域模式方法、基于Jacobi算法的近似联合对角化方法等有计算量较大等缺点。

发明内容

本发明针对现有技术的不足和改善需要，提供了一种基于相位同步的多类别运动想象脑电信号分类方法，该方法是利用相位锁定值分别计算训练样本和测试样本的相位同步特征并根据两者的相关系数依次将多对相应两个类别的脑电信号运用共空域模式算法通过线性判别分析进行多次分类来实现的。

本发明所采用的技术方案：

1.脑电信号采集：

设计实验，使用多通道脑电采集设备采集实验过程中被试者的脑电信号，完成被试者信息录入和脑电信号采集。

2.数据预处理：

对原始脑电信号进行数据预处理，目的是减少伪迹的干扰，提高信噪比，从而提高特征提取的准确性。预处理过程包括带通滤波。

3.滤波：

创建滤波器，将预处理后的脑电信号滤波到计算利用相位同步特征进行分类所需的频段。

4.相位同步特征相关系数计算：

本发明采用相位锁定值(PLV)来表示各次实验所得脑电信号每两个通道之间的相位关系，具体的PLV计算公式如下：

PLV＝|<exp(j{Φ_x(t)-Φ_y(t)})>|

其中Φ_x(t)和Φ_y(t)分别为脑电信号x(t)和y(t)的瞬时相位。

采用希尔伯特变换来计算信号的相位值，信号x(t)的希尔伯特变换定义如下：

$\tilde{x}\left(t\right)=\frac{1}{\mathrm{\&pi;}}P{\&Integral;}_{-\&infin;}^{+\&infin;}\frac{x\left(\mathrm{\&tau;}\right)}{t-\mathrm{\&tau;}}\mathrm{d\&tau;}$

公式中的P为柯西主值。由此可以定义x(t)的解析信号为：

$Z_x\left(t\right)=x\left(t\right)+j\tilde{x}\left(t\right)=A_x\left(t\right)e^{j{\mathrm{\&Phi;}}_x\left(t\right)}$

其中A_x(t)和Φ_x(t)分别为信号x(t)的瞬时振幅以及瞬时相位。

同样地，可以定义信号y(t)的解析信号，并计算瞬时相位Φ_y(t)。

本发明中，根据单次实验的时间流程和样本采样频率等采用合适的时间窗口来计算单次实验特定频段的PLV值。对于训练数据集而言，分别计算各类运动想象脑电信号PLV均值：假设训练数据集中存在N_i段属于第i类运动想象脑电信号样本的不重叠时间片段，则这N_i段片段的PLV均值计算如下：

${\mathrm{PLVavg}}_i=\frac{1}{N_i}|\underset{n=1}{\overset{N_i}{\mathrm{\&Sigma;}}}<\exp \left(\mathrm{j\&Delta;\&Phi;}\right)>|$

其中ΔΦ为信号x(t)和_y(t)的瞬时相位差。

假设本发明选定的脑电通道数为M，利用两两通道构建不同的通道对，计算某一类运动想象单次实验样本所有通道对的PLV值，此时得到一个M×M的上三角矩阵Κ

$K=\left[\begin{array}{ccccc}0& k_{12}& k_{13}& \&CenterDot;\&CenterDot;\&CenterDot;& k_{1M}\\ 0& 0& k_{23}& \&CenterDot;\&CenterDot;\&CenterDot;& k_{2M}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ 0& 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& k_{(M-1)M}\\ 0& 0& 0& 0& 0\end{array}\right]$

该矩阵不仅包含了不同脑电通道两两之间的相位关系，还包含了脑电通道的空间信息。

接下来，对于训练数据集而言，按上述方法计算属于这一类运动想象的其他单次实验样本所有通道对的PLV值，从而得到这一类别运动想象数据集的PLV均值，用同样的方法计算其他各类别运动想象训练数据集的PLV均值，而测试数据集只需计算单次实验样本所有通道对的PLV值。

最后，分别计算以上测试样本所有通道对PLV值所组成的矩阵Α与各类别运动想象数据集PLV均值所组成的矩阵Β的相位同步特征相关系数，具体的相关系数(r)计算公式如下：

$r=\frac{\underset{m}{\mathrm{\&Sigma;}}\underset{n}{\mathrm{\&Sigma;}}(A_{\mathrm{mn}}-\stackrel{\&OverBar;}{A})(B_{\mathrm{mn}}-\stackrel{\&OverBar;}{B})}{\sqrt{\left(\underset{m}{\mathrm{\&Sigma;}}\underset{n}{\mathrm{\&Sigma;}}{(A_{\mathrm{mn}}-\stackrel{\&OverBar;}{A})}^2\right)\left(\underset{m}{\mathrm{\&Sigma;}}\underset{n}{\mathrm{\&Sigma;}}{(B_{\mathrm{mn}}-\stackrel{\&OverBar;}{B})}^2\right)}}$

其中m,n为各矩阵的行列数，从而形成相关系数列向量R

$R=\left[\begin{array}{c}{r}_1\\ r_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ r_i\end{array}\right]$

其中r_i是相关系数，i是运动想象脑电信号类别数。

5.特征提取和分类：

本发明中，特征提取和分类的步骤大致如下：

第一步，基于相位同步特征相关系数的测试样本粗分类：将上述过程计算所得列向量R中的相关系数去平均后取其绝对值按从大到小排序形成列向量C，即

$C=\left[\begin{array}{c}{r}_1^{\&prime;}\\ r_2^{\&prime;}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ r_i^{\&prime;}\end{array}\right],r_1^{\&prime;}>r_2^{\&prime;}>...>r_i^{\&prime;}$

其中r_i'＝abs(r_i-mean(R)),abs()表示相关系数去平均后的绝对值，mean()表示相关系数的平均值，i是运动想象脑电信号类别数。

若原始脑电信号类别数为偶数，则把C中相关系数按顺序依次进行无重叠的两两组合，即每一个相关系数只能被组合一次，并把各两两组合相关系数所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行粗分类,将其归为两个类别中的某一类；若原始脑电信号类别数为奇数，则把C中相关系数按顺序依次进行无重叠的两两组合，此时C中最小的相关系数不参与组合，即最小的相关系数所对应的脑电信号类别不参与组合，并把各两两组合相关系数所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行粗分类,将其归为两个类别中的某一类。

第二步，利用上述相关系数两两组合方法所得到的测试样本可能会属于的新的脑电信号类别对测试样本进行细分类，按不同的情况具体操作如下：①若原始脑电信号类别数为偶数，且新的脑电信号类别数(即为下一次分类的原始脑电信号类别数)也为偶数，则把新的脑电信号类别按顺序依次进行无重叠的两两组合，即每一个新的脑电信号类别只能被组合一次，并把各两两组合所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行分类,将其归为两个类别中的某一类，形成新的脑电信号类别；若新的脑电信号类别数(即为下一次分类的原始脑电信号类别数)为奇数，则把新的脑电信号类别按顺序依次进行无重叠的两两组合，此时最后的脑电信号类别不参与组合，并把各两两组合所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行分类,将其归为两个类别中的某一类，形成新的脑电信号类别。②若原始脑电信号类别数为奇数，则把新的脑电信号类别和上一步骤中不参与组合的脑电信号类别组合成又一新的脑电信号类别(即为下一次分类的原始脑电信号类别)，然后根据这一新的脑电信号类别数的奇偶性按本步骤①中对新的脑电信号类别的处理方法进行操作。

第三步，根据第二步中形成的原始脑电信号类别数和新的脑电信号类别数的奇偶性，循环第二步的操作，直到测试样本被最终分类。

对上述步骤中涉及的CSP特征提取方法和LDA分类方法做如下介绍：

1)特征提取

本发明中，采用共空域模式(CSP)方法对两两组合的两个类别脑电数据集进行特征提取。假设Χ₁和Χ₂分别为两类想象运动任务下维数为N×T的多导联脑电信号,其中N为EEG信号的通道数,T为每一个通道采集点数。Χ₁和Χ₂归一化后的协方差矩阵R₁和R₂分别为:

$\begin{array}{cc}{R}_1=\frac{X_1{X}_1^T}{\mathrm{trace}\left(X_1{X}_1^T\right)}& R_2=\frac{X_2{X}_2^T}{\mathrm{trace}\left(X_2{X}_2^T\right)}\end{array}$

其中Χ^T表示矩阵Χ的转置,trace(Χ)为矩阵Χ对角线上元素的总和。对多次实验得到的空间协方差R₁和Rx取平均值,得到平均协方差矩阵和构造混合协方差矩阵R,并对R做主成分分析可得:

$R=\stackrel{\&OverBar;}{R_1}+\stackrel{\&OverBar;}{R_2}={\mathrm{U\&Sigma;U}}^T$

其中U是混合协方差矩阵R的特征向量组成的矩阵,Σ是R的对角线元为相应特征值的对角矩阵。定义白化矩阵Ρ为:

Ρ＝Σ^-1/2U^Τ

将矩阵和分别通过白化矩阵Ρ进行白化变换:

$\begin{array}{cc}{S}_1=P\stackrel{\&OverBar;}{R_1}{P}^T& S_2=P\stackrel{\&OverBar;}{R_2}{P}^T\end{array}$

然后再对S₁和S₂做主成分分析:

$\begin{array}{cc}{S}_1=U_1{\mathrm{\&Sigma;}}_1{U}_1^T& S_2=U_2{\mathrm{\&Sigma;}}_2{U}_2^T\end{array}$

可以证明矩阵S₁的特征向量矩阵U₁和矩阵S₂的特征向量矩阵U₂是相等的,即U₁＝U₂＝V；同时两个特征值的对角阵Σ₁和Σ₂之和为单位矩阵,即

Σ₁+Σ₂＝Ι

可知,当矩阵S₁的特征值最大时,对应的S₂的特征值最小；同理当S₂的特征值最大时,对应的S₁的特征值最小。根据这一特点,取V中前a行和后a行组成最佳空间滤波器SF,并得到投影矩阵W

W＝SF^ΤΡ

然后将原始EEG信号通过投影矩阵进行投影分别得到两个类别的新信号：

Ζ₁＝WΧ₁ Ζ₂＝WΧ₂

另外，假设Χ为测试样本，则测试样本的新信号为：

Ζ＝WΧ

最后对各类别新信号各维度上的方差进行取对数处理，作为特征：

$f_l=\log \left(\frac{\mathrm{var}\left(Z_l\right)}{{\mathrm{\&Sigma;}}_{l=1}^{2a}\left(\mathrm{var}\left(Z_l\right)\right)}\right),l=1,\&CenterDot;\&CenterDot;\&CenterDot;,2a.$

其中var()表示方差。

2)分类：

本发明中，采用线性判别分析(LDA)对所提取的特征进行分类。

对已知类别的各两两组合的两个类别脑电数据集按上述方法进行特征提取，形成的特征向量作为训练样本，通过训练样本的训练可以得到分类所需的判别函数：

F(x_J)＝w^Τx_J+ω₀

其中{x_J|J＝1,2,...,D}为第J个测试样本的特征向量,D为测试样本的个数。参数w和ω₀的计算准则为：使不同类别间的距离最大，同一类别内的距离最小，即使不同的类别尽可能的分开。

假设X＝{x₁,x₂,...,x_d}为训练样本的数据集，数据集中的每个特征向量分别对应类别{X₁,...,X_c,...,X_C}中的一个类别，则类间离散度矩阵S_b和类内离散度矩阵S_w分别为：

$S_b=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{M}_c(m_c-m){(m_c-m)}^T$

$S_w=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}\underset{x\&Element;X_c}{\mathrm{\&Sigma;}}(x-m_c){(x-m_c)}^T$

其中M_c为类别X_c的向量数量，m_c为类均值向量，m为所有样本的整体均值。

计算矩阵W使类间离散度矩阵和类内离散度矩阵的比值最大：

$W_{\mathrm{opt}}\text{=arg}\underset{w}{\max }\frac{\left|W^T{S}_{b}W\right|}{\left|W^T{S}_{w}W\right|}=[w_1,w_2,...,w_N]$

即可求得所需的判别函数。

利用这些训练样本所得到的判别函数就可以对未知类别的测试样本进行分类。

本发明的有益效果是：提供了一种基于相位同步的能够有效地对多类别运动想象脑电信号进行分类的方法。该方法中基于相位同步的脑电信号的粗分类可以高效率、高准确率地确定测试样本大致所属的类别，减少计算量，为下一步的基于CSP方法的特征提取与分类提供有效的保障，具有更高的实时性和稳定性。与其他传统的“一对一”等CSP方法相比，基于相位同步的脑电信号的分类方法更加高效和实用。

附图说明

图1为本发明流程图；

图2为图1的具体实施流程图，其中图(2)是图(1)中“特征提取和分类”步骤的具体实施流程图；

图3为脑电通道图。

具体实施方式

以下结合附图对本发明作进一步说明。

如图1所示：首先采集所需的脑电信号，接着对脑电信号进行滤波等预处理，然后将预处理后的数据滤波到特定频段计算相位锁定值。通过利用相位锁定值计算所得到的相位同步特征相关系数和共空域模式(CSP)特征提取方法以及线性判别分析分类方法三者的结合运用对测试样本进行分类。

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

步骤S1：通过多通道脑电采集设备采集所需的脑电信号。本实施例中，采样频率为250Hz，电极帽采用国际10/20系统电极放置法，22个电极分别为Fz、FC3、FC1、FCz、FC2、FC4、C5、C3、C1、Cz、C2、C4、C6、CP3、CP1、CPz、CP2、CP4、P1、Pz、P2、POz

参考电极置于左耳垂，右耳垂作为地。22个脑电通道按以上顺序排列的位置如图3所示。

步骤S2：对采集的脑电信号进行预处理，以减少伪迹的干扰。具体过程如下：

带通滤波：所有采集到的脑电信号经过250Hz的采样率采样后通过0.5-100Hz的带通滤波器进行滤波,同时放大器的灵敏度设置为100μV。

步骤S3：利用3阶巴特沃斯带通滤波器将预处理后的脑电信号滤波到8-30Hz频段。

步骤S4：计算相位同步特征相关系数，具体包括：

1)利用相位锁定值(PLV)计算各次实验所得脑电信号每两个通道之间的相位关系，计算公式如下：

PLV＝|<exp(j{Φ_x(t)-Φ_y(t)})>|

其中Φ_x(t)和Φ_y(t)分别为脑电信号x(t)和y(t)的瞬时相位。本实施例中，采用3秒的时间窗口来计算选定通道的PLV值。

2)通过希尔伯特变化计算脑电信号的瞬时相位值，信号x(t)的希尔伯特变换为：

$\tilde{x}\left(t\right)=\frac{1}{\mathrm{\&pi;}}P{\&Integral;}_{-\&infin;}^{+\&infin;}\frac{x\left(\mathrm{\&tau;}\right)}{t-\mathrm{\&tau;}}\mathrm{d\&tau;}$

其中P为柯西主值。由此可得x(t)的解析信号为：

$Z_x\left(t\right)=x\left(t\right)+j\tilde{x}\left(t\right)=A_x\left(t\right)e^{j{\mathrm{\&Phi;}}_x\left(t\right)}$

其中A_x(t)和Φ_x(t)分别为信号x(t)的瞬时振幅以及瞬时相位。

同样地，可以定义信号y(t)的解析信号，并计算瞬时相位Φ_y(t)。

3)本实施例中，共有4个类别的运动想象脑电信号。对于训练数据集而言，分别计算这4个类别运动想象脑电信号各自的PLV均值：假设训练数据集中存在N_i段属于第i类运动想象脑电信号样本的不重叠时间片段，则这N_i段片段的PLV均值计算如下：

${\mathrm{PLVavg}}_i=\frac{1}{N_i}|\underset{n=1}{\overset{N_i}{\mathrm{\&Sigma;}}}<\exp \left(\mathrm{j\&Delta;\&Phi;}\right)>|$

本实施例中，每个类别不重叠时间片段的段数随着训练样本选取的不同而不同。

利用两两通道构建不同的通道对，本实施例中共有22个不同的脑电通道，可以构建231个不同的通道对。计算各类别运动想象训练数据集的PLV均值和测试数据集单次实验样本所有通道对的PLV值，均可放置为一个22×22的上三角矩阵Ρ，该矩阵不仅包含了不同脑电通道两两之间的相位关系，还包含了脑电通道的空间信息。

4)分别计算以上测试样本所有通道对PLV值所组成的矩阵Α与各类别运动想象数据集PLV均值所组成矩阵Β的相关系数。本实施例中，共有4个类别的运动想象脑电信号，因此对于每个测试样本会形成4个相关系数。

步骤S5：特征提取和分类。

本实施例中，特征提取和分类的步骤大致如下：

第一步，基于相位同步特征相关系数的测试样本粗分类：将上述过程计算所得的4个相关系数去平均后取其绝对值按从大到小排序形成列向量C，把C中前两个相关系数所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行粗分类，将其归为两个类别中的某一类；把C中后两个相关系数所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行粗分类，将其归为两个类别中的某一类。

第二步，利用上述相关系数两两组合方法所得到的测试样本可能会属于的两个新的脑电信号类别对测试样本进行细分类，把这两个类别所对应的脑电数据集作为训练样本进行CSP特征提取，并再次对测试样本用LDA分类方法进行分类，得到最终的一个分类结果。

对上述步骤中涉及的CSP特征提取方法和LDA分类方法做如下介绍：

1)特征提取

采用共空域模式(CSP)方法对两两组合的两个类别脑电数据集进行特征提取。本实施例中，假设Χ₁和Χ₂分别为两类想象运动任务下维数为22×750的多导联脑电信号。Χ₁和Χ₂归一化后的协方差矩阵R₁和R₂分别为:

$\begin{array}{cc}{R}_1=\frac{X_1{X}_1^T}{\mathrm{trace}\left(X_1{X}_1^T\right)}& R_2=\frac{X_2{X}_2^T}{\mathrm{trace}\left(X_2{X}_2^T\right)}\end{array}$

其中Χ^T表示矩阵Χ的转置,trace(Χ)为矩阵Χ对角线上元素的总和。对多次实验得到的空间协方差R₁和R₂取平均值,得到平均协方差矩阵和构造混合协方差矩阵R,并对R做主成分分析可得:

$R=\stackrel{\&OverBar;}{R_1}+\stackrel{\&OverBar;}{R_2}={\mathrm{U\&Sigma;U}}^T$

其中U是混合协方差矩阵R的特征向量组成的矩阵,Σ是R的对角线元为相应特征值的对角矩阵。定义白化矩阵Ρ为:

Ρ＝Σ^-1/2U^Τ

将矩阵和分别通过白化矩阵Ρ进行白化变换:

$\begin{array}{cc}{S}_1=P\stackrel{\&OverBar;}{R_1}{P}^T& S_2=P\stackrel{\&OverBar;}{R_2}{P}^T\end{array}$

然后再对S₁和S₂做主成分分析:

$\begin{array}{cc}{S}_1=U_1{\mathrm{\&Sigma;}}_1{U}_1^T& S_2=U_2{\mathrm{\&Sigma;}}_2{U}_2^T\end{array}$

可以证明矩阵S₁的特征向量矩阵U₁和矩阵S₂的特征向量矩阵U₂是相等的,即U₁＝U₂＝V；同时两个特征值的对角阵Σ₁和Σ₂之和为单位矩阵,即

Σ₁+Σ₂＝Ι

可知,当矩阵S₁的特征值最大时,对应的S₂的特征值最小；同理当S₂的特征值最大时,对应的S₁的特征值最小。根据这一特点,本实施例中，取V中前3行和后3行组成最佳空间滤波器SF,并得到投影矩阵W

W＝SF^ΤΡ

然后将原始EEG信号通过投影矩阵进行投影分别得到两个类别的新信号：

Ζ₁＝WΧ₁ Ζ₂＝WΧ₂

另外，假设Χ为测试样本，则测试样本的新信号为：

Ζ＝WΧ

最后对各类别新信号各维度上的方差进行取对数处理，作为特征：

$f_l=\log \left(\frac{\mathrm{var}\left(Z_l\right)}{{\mathrm{\&Sigma;}}_{l=1}^6\left(\mathrm{var}\left(Z_l\right)\right)}\right),l=1,\&CenterDot;\&CenterDot;\&CenterDot;,6.$

其中var()表示方差。

2)分类：

采用线性判别分析(LDA)对所提取的特征进行分类。

本实施例中，共有288个实验样本，随机选取144个作为训练样本，其余144个作为测试样本。把144个训练样本按类别分成4个不同类别的数据集，按上述方法对各两两组合成的两个类别脑电数据集进行特征提取，形成6×I(I为两类别脑电数据集样本总数)维的特征向量作为训练样本，从而得到所需的判别函数：

F(x_J)＝w^Τx_J+ω₀

其中{x_J|J＝1,2,...,D}为第J个测试样本的特征向量,D为测试样本的个数。参数w和ω₀的计算准则为：使不同类别间的距离最大，同一类别内的距离最小，即使不同的类别尽可能的分开。

假设X＝{x₁,x₂,...,x_d}为训练样本的数据集，数据集中的每个特征向量分别对应类别{X₁,...,X_c,...,X_C}中的一个类别，则类间离散度矩阵S_b和类内离散度矩阵S_w分别为：

$S_b=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{M}_c(m_c-m){(m_c-m)}^T$

$S_w=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}\underset{x\&Element;X_c}{\mathrm{\&Sigma;}}(x-m_c){(x-m_c)}^T$

其中M_c为类别X_c的向量数量，m_c为类均值向量，m为所有样本的整体均值。

计算矩阵W使类间离散度矩阵和类内离散度矩阵的比值最大：

$W_{\mathrm{opt}}\text{=arg}\underset{w}{\max }\frac{\left|W^T{S}_{b}W\right|}{\left|W^T{S}_{w}W\right|}=[w_1,w_2,...,w_N]$

即可求得所需的判别函数。

利用这些训练样本所得到的判别函数就可以对未知类别的测试样本进行分类。最后将分类结果与这些特征向量的实际类别进行比较，得到该脑电信号分类方法的分类准确率。

Claims (1)

1.一种基于相位同步的多类别运动想象脑电信号分类方法，其特征在于该方法包括以下步骤：

步骤1、脑电信号采集：使用多通道脑电采集设备采集实验过程中被试者的脑电信号，完成被试者信息录入和脑电信号采集；

步骤2、数据预处理：对原始脑电信号进行数据预处理，具体为带通滤波；

步骤3、滤波：创建滤波器，将预处理后的脑电信号滤波到计算利用相位同步特征进行分类所需的频段；

步骤4、相位同步特征相关系数计算：采用相位锁定值PLV来计算各次实验所得脑电信号每两个通道之间的相位关系，具体的PLV计算公式如下：

PLV＝|<exp(j{Φ_x(t)-Φ_y(t)})>|

其中Φ_x(t)和Φ_y(t)分别为脑电信号x(t)和y(t)的瞬时相位；

采用希尔伯特变换来计算信号的相位值，信号x(t)的希尔伯特变换x定义如下：

$\tilde{x}\left(t\right)=\frac{1}{\mathrm{\&pi;}}P{\&Integral;}_{-\&infin;}^{+\&infin;}\frac{x\left(\mathrm{\&tau;}\right)}{t-\mathrm{\&tau;}}\mathrm{d\&tau;}$

公式中的P为柯西主值；由此可以定义x(t)的解析信号为：

$Z_x\left(t\right)=x\left(t\right)+j\tilde{x}\left(t\right)=A_x\left(t\right)e^{j{\mathrm{\&Phi;}}_x\left(t\right)}$

其中A_x(t)和Φ_x(t)分别为信号x(t)的瞬时振幅以及瞬时相位；

同样地，可以定义信号y(t)的解析信号，并计算瞬时相位Φ_y(t)；

根据单次实验的时间流程和样本采样频率采用时间窗口来计算单次实验特定频段的PLV值；对于训练数据集而言，分别计算各类运动想象脑电信号PLV均值：假设训练数据集中存在N_i段属于第i类运动想象脑电信号样本的不重叠时间片段，则这N_i段片段的PLV均值计算如下：

${\mathrm{PLVavg}}_i=\frac{1}{N_i}|\underset{n=1}{\overset{N_i}{\mathrm{\&Sigma;}}}<\exp \left(\mathrm{j\&Delta;\&Phi;}\right)>|$

其中ΔΦ为信号x(t)和y(t)的瞬时相位差；

设选定的脑电通道数为M，利用两两通道构建不同的通道对，计算某一类运动想象单次实验样本所有通道对的PLV值，此时得到一个M×M的上三角矩阵Κ

$K=\left[\begin{array}{ccccc}0& k_{12}& k_{13}& \&CenterDot;\&CenterDot;\&CenterDot;& k_{1M}\\ 0& 0& k_{23}& \&CenterDot;\&CenterDot;\&CenterDot;& k_{2M}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ 0& 0& 0& \&CenterDot;\&CenterDot;\&CenterDot;& k_{(M-1)M}\\ 0& 0& 0& 0& 0\end{array}\right]$

该矩阵不仅包含了不同脑电通道两两之间的相位关系，还包含了脑电通道的空间信息；

接下来，对于训练数据集而言，按上述方法计算属于这一类运动想象的其他单次实验样本所有通道对的PLV值，从而得到这一类别运动想象数据集的PLV均值，用同样的方法计算其他各类别运动想象训练数据集的PLV均值，而测试数据集只需计算单次实验样本所有通道对的PLV值；

最后，分别计算以上测试样本所有通道对PLV值所组成的矩阵Α与各类别运动想象数据集PLV均值所组成的矩阵Β的相关系数，具体的相关系数(r)计算公式如下：

$r=\frac{\underset{m}{\mathrm{\&Sigma;}}\underset{n}{\mathrm{\&Sigma;}}(A_{\mathrm{mn}}-\stackrel{\&OverBar;}{A})(B_{\mathrm{mn}}-\stackrel{\&OverBar;}{B})}{\sqrt{\left(\underset{m}{\mathrm{\&Sigma;}}\underset{n}{\mathrm{\&Sigma;}}{(A_{\mathrm{mn}}-\stackrel{\&OverBar;}{A})}^2\right)\left(\underset{m}{\mathrm{\&Sigma;}}\underset{n}{\mathrm{\&Sigma;}}{(B_{\mathrm{mn}}-\stackrel{\&OverBar;}{B})}^2\right)}}$

其中m,n为各矩阵的行列数，从而形成相关系数列向量R

$R=\left[\begin{array}{c}{r}_1\\ r_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ r_i\end{array}\right]$

其中r_i是相关系数，i是运动想象脑电信号类别数；

步骤5、特征提取和分类：特征提取和分类的步骤大致如下：

第一步，基于相位同步特征相关系数的测试样本粗分类：将上述过程计算所得列向量R中的相关系数去平均后取其绝对值按从大到小排序形成列向量C，即

$C=\left[\begin{array}{c}{r}_1^{\&prime;}\\ r_2^{\&prime;}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ r_i^{\&prime;}\end{array}\right],r_1^{\&prime;}>r_2^{\&prime;}>...>r_i^{\&prime;}$

其中r_i'＝abs(r_i-mean(R)),abs()表示相关系数去平均后的绝对值，mean()表示相关系数的平均值,i是运动想象脑电信号类别数；

若原始脑电信号类别数为偶数，则把C中相关系数按顺序依次进行无重叠的两两组合，即每一个相关系数只能被组合一次，并把各两两组合相关系数所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行粗分类,将其归为两个类别中的某一类；若原始脑电信号类别数为奇数，则把C中相关系数按顺序依次进行无重叠的两两组合，此时C中最小的相关系数不参与组合，即最小的相关系数所对应的脑电信号类别不参与组合，并把各两两组合相关系数所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行粗分类,将其归为两个类别中的某一类；

第二步，利用上述相关系数两两组合方法所得到的测试样本可能会属于的新的脑电信号类别对测试样本进行细分类，按不同的情况具体操作如下：①若原始脑电信号类别数为偶数，且新的脑电信号类别数也为偶数，则把新的脑电信号类别按顺序依次进行无重叠的两两组合，即每一个新的脑电信号类别只能被组合一次，并把各两两组合所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行分类,将其归为两个类别中的某一类，形成新的脑电信号类别；若新的脑电信号类别数为奇数，则把新的脑电信号类别按顺序依次进行无重叠的两两组合，此时最后的脑电信号类别不参与组合，并把各两两组合所对应的两个类别脑电数据集作为训练样本进行CSP特征提取并对测试样本用LDA分类方法进行分类,将其归为两个类别中的某一类，形成新的脑电信号类别；②若原始脑电信号类别数为奇数，则把新的脑电信号类别和上一步骤中不参与组合的脑电信号类别组合成又一新的脑电信号类别，然后根据这一新的脑电信号类别数的奇偶性按本步骤①中对新的脑电信号类别的处理方法进行操作；

第三步，根据第二步中形成的原始脑电信号类别数和新的脑电信号类别数的奇偶性，循环第二步的操作，直到测试样本被最终分类；

对上述步骤中涉及的CSP特征提取方法和LDA分类方法具体是：

1)CSP特征提取

采用共空域模式(CSP)方法对两两组合的两个类别脑电数据集进行特征提取；假设Χ₁和Χ₂分别为两类想象运动任务下维数为N×T的多导联脑电信号,其中N为EEG信号的通道数,T为每一个通道采集点数；Χ₁和Χ₂归一化后的协方差矩阵R₁和R₂分别为：

$\begin{array}{cc}{R}_1=\frac{X_1{X}_1^T}{\mathrm{trace}\left(X_1{X}_1^T\right)}& R_2=\frac{X_2{X}_2^T}{\mathrm{trace}\left(X_2{X}_2^T\right)}\end{array}$

其中Χ^T表示矩阵Χ的转置,trace(Χ)为矩阵Χ对角线上元素的总和；对多次实验得到的空间协方差R₁和R₂取平均值,得到平均协方差矩阵R₁和构造混合协方差矩阵R,并对R做主成分分析可得:

$R=\stackrel{\&OverBar;}{R_1}+\stackrel{\&OverBar;}{R_2}={\mathrm{U\&Sigma;U}}^T$

其中U是混合协方差矩阵R的特征向量组成的矩阵,Σ是R的对角线元为相应特征值的对角矩阵；定义白化矩阵Ρ为:

Ρ＝Σ^-1/2U^Τ

将矩阵和分别通过白化矩阵Ρ进行白化变换:

$\begin{array}{cc}{S}_1=P\stackrel{\&OverBar;}{R_1}{P}^T& S_2=P\stackrel{\&OverBar;}{R_2}{P}^T\end{array}$

然后再对S₁和S₂做主成分分析:

$\begin{array}{cc}{S}_1=U_1{\mathrm{\&Sigma;}}_1{U}_1^T& S_2=U_2{\mathrm{\&Sigma;}}_2{U}_2^T\end{array}$

可以证明矩阵S₁的特征向量矩阵U₁和矩阵S₂的特征向量矩阵U₂是相等的,即U₁＝U₂＝V；同时两个特征值的对角阵Σ₁和Σ₂之和为单位矩阵,即

Σ₁+Σ₂＝Ι

可知,当矩阵S₁的特征值最大时,对应的S₂的特征值最小；同理当S₂的特征值最大时,对应的S₁的特征值最小；根据这一特点,取V中前a行和后a行组成最佳空间滤波器SF,并得到投影矩阵W

W＝SF^ΤΡ

然后将原始EEG信号通过投影矩阵进行投影分别得到两个类别的新信号：

Ζ₁＝WΧ₁ Ζ₂＝WΧ₂

另外，假设Χ为测试样本，则测试样本的新信号为：

Ζ＝WΧ

最后对各类别新信号各维度上的方差进行取对数处理，作为特征：

$f_l=\log \left(\frac{\mathrm{var}\left(Z_l\right)}{{\mathrm{\&Sigma;}}_{l=1}^{2a}\left(\mathrm{var}\left(Z_l\right)\right)}\right),l=1,\&CenterDot;\&CenterDot;\&CenterDot;,2a;$

其中var()表示方差；

2)LDA分类

对已知类别的各两两组合的两个类别脑电数据集进行特征提取，形成的特征向量作为训练样本，通过训练样本的训练可以得到分类所需的判别函数：

F(x_J)＝w^Τx_J+ω₀

其中{x_J|J＝1,2,...,D}为第J个测试样本的特征向量,D为测试样本的个数；参数w和ω₀的计算准则为：使不同类别间的距离最大，同一类别内的距离最小，即使不同的类别尽可能的分开；

假设X＝{x₁,x₂,...,x_d}为训练样本的数据集，数据集中的每个特征向量分别对应类别{X₁,...,X_c,...,X_C}中的一个类别，则类间离散度矩阵S_b和类内离散度矩阵S_w分别为：

$S_b=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{M}_c(m_c-m){(m_c-m)}^T$

$S_w=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}\underset{x\&Element;X_c}{\mathrm{\&Sigma;}}(x-m_c){(x-m_c)}^T$

其中M_c为类别X_c的向量数量，m_c为类均值向量，m为所有样本的整体均值；

计算矩阵W使类间离散度矩阵和类内离散度矩阵的比值最大：

$W_{\mathrm{opt}}\text{=arg}\underset{w}{\max }\frac{\left|W^T{S}_{b}W\right|}{\left|W^T{S}_{w}W\right|}=[w_1,w_2,...,w_N]$

即可求得所需的判别函数；

利用这些训练样本所得到的判别函数就可以对未知类别的测试样本进行分类。