CN111616701A - 基于多元变分模态分解的脑电多域特征提取方法 - Google Patents

Abstract

本发明公开了一种基于多元变分模态分解的脑电多域特征提取方法，本发明首先利用MVMD对原始脑电多通道数据进行自适应分解，然后从分解得到的固有模态函数(IMF)分量提取信号的时域特征以及非线性动力学特征，同时将IMF分量合并构造新的信号矩阵，并采用共空间模式(CSP)法对重构信号提取空间特征，进行时域、非线性动力学以及空域特征的结合，最后通过支持向量机(SVM)对此特征集分类。本发明能有效判别与特定任务相关的重要信息成分，为脑电信号的特征提取提供了新的思路。

Description

基于多元变分模态分解的脑电多域特征提取方法

技术领域

本发明属于模式识别领域，是针对运动想象脑电信号(EEG)，提出利用多元变分模态分解(MVMD)对原始脑电多通道数据进行自适应分解，然后从分解得到的固有模态函数(IMF)分量提取信号的时域特征以及非线性动力学特征，同时将IMF分量合并构造新的信号矩阵，并采用共空间模式(CSP)法对重构信号提取空间特征，来进行时域、非线性动力学以及空域特征的结合的特征及分类的方法。

背景技术

脑机接口技术(BCI)是一种不依赖于人体神经和肌肉组织的正常传输通路，而直接进行人脑与外界之间信息交流的人机交互系统，对于活动能力缺失患者的能力恢复和功能训练具有重要意义，可以为意识障碍和中风患者的康复提供很大的帮助，患者可以通过该技术实现机械设备的控制，完成所需的动作。基于运动想象的BCI是研究最为广泛的一类，运动想象型BCI通过捕捉和识别不同运动想象任务下的EEG信号来实现大脑与外界的信息交换和控制。然而，EEG信号具有非平稳性、非线性等特点，单一角度的特征不能有效体现出信号特点的全面性，因此如何从EEG信号中提取出能全面且有效的识别运动任务的特征，对BCI系统的性能至关重要。

判别运动想象EEG信号通常需要三个步骤：预处理、特征提取、分类器分类。而特征提取往往是最重要的环节，提取的特征好坏直接决定了分类的效果，因此如何提取有利分类的特征是研究的重要课题。

近年来，针对传统特征提取算法在运动想象脑电信号时频特性分析的不足，1998年Huang等人提出了经典的由数据驱动进行分解的时频分析法，称为经验模态分解(Empirical Mode Decomposition,EMD)算法，并且定义了固有模态函数(Intrinsic ModeFunctions,IMF)的概念，这种方法具有很强的自适应性，在信号分析领域受到广泛关注。但EMD存在模态混叠现象，Wu和Huang针对此问题，提出利用噪声辅助EMD的算法EEMD(Ensemble EMD)，因为白噪声在频域内其功率谱密度分布均匀，被分析信号与白噪声混合可以改变极值点特性，补充一些缺失的尺度，在信号分解中有较好的表现。但是，EEMD和EMD不适用于同时处理多通道脑电数据，造成分解得到的IMFs个数以及频率的不匹配现象，2010年，Rehman等人提出多元经验模态分解MEMD(Multivariate EMD)，改进了EMD算法，它能够同时分解多通道的脑电信号，但是仍然没有解决模态混叠的问题。2014年，Dragomiretskiy等提出变分模态分解(variational mode decomposition,VMD)，将信号的分解引入到变分模型中，利用寻找约束变分模型最优解的过程实现信号的分解，同样定义分解得到的分量为固有模态函数(IMF)分量，其中每个分量都具有平稳性，是有限带宽的调频调幅信号。这种方法避免了因局部突变而引起的模态混叠现象，并且在计算速度上具有一定优势，但与EMD算法一样，存在不能同时处理多通道脑电信号的限制，2019年Naveed urRehman等提出的多元变分模态分解MVMD(multivariate variational modedecomposition,MVMD)算法，将VMD算法从单通道拓展到多通道，可以同时处理多通道数据，避免了频率不匹配问题，又解决了MEMD的模式混叠问题，通过MVMD算法对多通道信号自适应地分解成多个具有物理意义的IMF，接着对各阶IMF作Hilbert变换，可以获得高分辨率的时频特征。

此外，通过非线性动力学特征分析脑电信号也是一种有效的处理方法，其中近似熵(Approximate Entropy,AE)、样本熵(Sample Entropy,SE)是近年来脑电信号在非线性动力学研究中的特征代表，Richman和Moorman等提出样本熵，解决了近似熵的计算过程中“自我比较”引起偏差的问题，但没有考虑到时间序列中可能存在的不同时间尺度特征，为了计算不同时间尺度下信号的复杂性，Costa等提出了多尺度熵(MultiScale Entropy)的算法，将样本熵扩展到多个时间尺度，有效的避免原数据进行线性平滑的过程中丢失有效的特征信息。尽管时频特征以及多尺度熵特征已表现出良好的效果，但随着大脑状态的变化特征会有波动，如果可以同步利用EEG信号的空间相关性，便能在一定程度上消除信号的噪声，并实现局部皮层神经活动的定位，共空间模式(Common spatial pattern,CSP)便是一种应用于分析多维脑电信号的空间滤波算法，可以补充EEG信号的空域特征。

近年来相关研究者综合考虑多种特征方法进行特征提取。Chen S等提出了融合香农熵、小波熵和样本熵进行特征提取；杨默涵等基于总体经验模态分解和近似熵提出一种多特征提取方法。以上方法都表现出很好的自适应性和较高的识别准确率，但都忽略了每一类运动想象脑电信号中重要的时域特征和空间分布成分。为了得到更全面有效的特征，本发明提出基于MVMD的时域特征，非线性动力学特征以及空域特征结合的多域特征提取方法。

发明内容

针对运动想象脑电信号具有非平稳性、非线性等特点，单一角度的特征不能有效体现出信号特点的全面性，本发明提出了基于MVMD的时域特征，非线性动力学特征以及空域特征结合的多域特征提取方法，并基于此完成运动想象脑电信号分类的方法。

本发明包括以下步骤：

步骤(1)：多通道脑电信号采集并进行预处理。

采集受试者的多通道运动想象脑电信号；接着利用带通滤波器对脑电信号进行8～30HZ的低通滤波，用于脑电信号特征提取。

步骤(2)：对多通道脑电信号进行多元变分模态分解(MVMD)，得到每个通道的K个本征模态函数(IMF)。

步骤(3)：对各IMF分量采用在时域具有很高分辨率的希尔伯特-黄变换(HHT)来提取瞬时幅值，并求瞬时能量，以此作为区分两类运动想象脑电信号的时域特征。

步骤(4)：对各IMF分量提取多尺度样本熵作为区分两类运动想象脑电信号的非线性动力学特征。

步骤(5)：将每个IMF分量看作为单个通道的脑电信号，并以原始信号的采样频率对其采样，然后合并分量构造两类运动想象任务下新的信号矩阵，对其进行CSP分解得到区分两类的方差向量。

步骤(6)：将步骤(3)～(5)中求得的特征值结合，构造输入特征向量，并对输入向量进行归一化。将训练集特征输入SVM分类器训练分类模型后，对测试集特征进行分类得到测试结果。

本发明的有益效果：运动想象脑电信号通常具有非平稳性、非线性等特点，而传统的方法采用的都是单个角度特征进行特征提取，不能较好地体现有效特征的全面性，本发明从多个角度获取信号的更完整描述，且信号特征都是基于MVMD算法，保证了信号稳定地分解为多个具有物理意义的IMF分量，为获得有效的多特征创造了条件。

附图说明

图1为本发明的流程图；

图2(a)为C3通道运动想象脑电信号MVMD分解后的分量波形；

图2(b)为C4通道运动想象脑电信号MVMD分解后的分量波形；

图3(a)为C3通道特征向量F1的特征值；

图3(b)为C3通道特征向量F1的特征值统计结果。

具体实施方式

下面结合附图详细描述本发明基于MVMD的脑电多域特征提取，如图1，本发明的实施主要包括6个步骤：(1)采集多通道脑电信号及预处理，(2)对多通道脑电信号进行MVMD分解，得到若干本征模态函数IMF分量，(3)根据各IMF分量的Hilbert谱分析求出脑电信号的瞬时能量均值特征，(4)对各IMF分量提取多尺度样本熵特征，(5)将由分量构成的新的信号矩阵通过CSP算法提取方差向量，(6)结合(3)(4)(5)所求特征输入分类器分类得到结果。

下面逐一对各步骤进行详细说明。

步骤(1)：本实例采用BCI Competition II Dataset III数据集进行验证，数据按照以下方式进行采集。实验数据全部来源于同一受试者(女性，25岁)，由位于C3、C4和Cz三个电极前后各2.5cm位置的双级导联记录EEG信号。实验内容是根据提示来想象左右手运动，共包括7组，每组40次实验，每次实验持续9s，前两秒钟被测试者处于安静状态，不做任何思维想象，在第二秒是计算机发出由低到高的提示音，计算机屏幕中央出现“+”图像。第三秒时屏幕出现向左或向右的箭头，提示被测试者想象用左(右)手按箭头指向去移动屏幕上出现的条状棒。采样率为128HZ，并进行了0.5～30HZ带通滤波。在进行的280次实验中，其中140次做训练数据，另外140次为测试数据。

数据预处理：对所有的数据集进行相同的预处理步骤。每位受试者分别选择在视觉提示后的3s之后到8s作为实验数据，然后利用6阶巴特沃兹滤波器对数据进行8～30HZ的带通滤波处理以滤除不相关的低频噪声。

步骤(2)：对多通道数据进行MVMD分解。具体为：

MVMD方法获取IMF分量的方法不同于EMD的循环筛选分解的方法，而是基于输入数据的所有通道之间存在的共同频率分量，在信号的分解过程引入了变分模型，构造一个求变分模型最优解的问题，在这个过程中，各通道的IMF分量同时迭代更新其中心频率和带宽，从而自适应地得到K个窄带IMF分量。对于包含C个数据通道的输入数据X(t)，记为[x₁(t),x₂(t),…x_C(t)]，MVMD算法的具体步骤如下：

①、预先定义k个窄带IMF分量u_k(t)，使

其中，u_k(t)＝[u₁(t),u₂(t),…u_C(t)]。

②、利用希尔伯特变换的算法分别得到向量u_k(t)中每一个元素的解析表示，记为

以计算单边频谱，然后与指数项

相乘来调整其对应的中心频率ω_k(t)，使每个模态的频谱被调制到与它相对应的基频带上，通过谐波转换后

的梯度函数的L₂范数来估计出各模态u_k(t)的带宽。

由于单个的频率分量ω_k被用作为整个向量

的谐波混频，因此需要在多元振荡u_k(t)中找到在多个通道的共同频率分量ω_k，考虑

的所有通道，使每个通道所分解的IMF分量之和能够再现该输入信号，且最小化模态函数带宽之和，受约束的相关优化问题变为：

其中，u_k,c(t)表示相应的通道数c和向量u_k(t)中的元素，

表示对应的u_k,c(t)中每一个元素的解析表示，

表示与时间相关的偏导数。

③、对变分问题求解，构造增广的拉格朗日表示为：

其中，α为二次惩罚因子，其保障存在高斯噪声时信号的重构精度，λ_c(t)为拉格朗日乘子，使得约束条件保持严格性。

④、为了解决这个转变后的非约束性的变分问题，应用乘子交替方向法来更新

和

以寻找增广拉格朗日表达式的鞍点。具体做法为：

1)模态更新

模态更新的最小化问题相当于：

利用Parseval傅里叶变换，将其变换到频域，与VMD算法的模式更新方法类似，得到模式更新为：

2)中心频率更新

由于增广拉格朗日函数的后两项不依赖于ω_k，则中心频率更新所对应的优化问题可以简化为：

同样利用Parseval傅里叶变换，将该优化转换到频域，得到在傅里叶域中的等效问题为：

得到中心频率的更新方法：

通过更新关系自适应地分解信号的频带，得到k个窄带IMF分量。此外，由于MVMD方法能同时计算多个通道数据，保证了通道间频率的一致性，使信号的分析更趋于稳定；

步骤(3)：将式(1)的各分量u_k(t)进行希尔伯特谱分析，定义为u_k(t)与1/πt的卷积，即：

其中，μ为积分变量。

由u_k(t)作为实部，它的希尔伯特变换作为虚部，得到其解析信号和瞬时幅值为：

U_k(t)＝u_k(t)+j₀H[u_k(t)] (10)

其中，j₀为虚数单位。

然后计算每个通道的平均瞬时能量值为：

其中n为采样点数目，

为第i个采样点对应的瞬时能量值。取瞬时能量均值作为时域特征，记为F₁；

步骤(4)：

①根据采样点数对IMF分量作一维离散时间序列{u₁,u₂,…,u_i,…,u_n}，然后进行多尺度化处理来构造粗粒化形式的时间序列{y^(τ)}：

其中，

τ为尺度因子，序列长度为M＝int(n/τ)。

②构造一组m维的向量：

X_m(i)＝{y_i+b:0≤b≤m-1} (14)

③计算X(i)和X(j)间距离:

d[X(i),X(j)]＝max|y_(i+b)-y_(i+b)| (15)

其中，0≤b≤m-1；i,j＝1,2...M-m+1；i≠j，M为序列长度。

④计算d[X(i),X(j)]中小于阈值r的个数，也是模板匹配的个数，记为B^m(i)，在计算B^m(i)与距离总数的比值，记作：

⑤计算

的平均值C^m(r)：

⑥将维数增加至m+1，重复②～⑤步骤，计算

和C^m+1(r)。

得到时间序列长度为M时的样本熵值：

SampEn(m,r,M)＝-ln[C^m+1(r)/C^m(r)] (18)

重复以上步骤得到多尺度下的样本熵特征，记为F₂；

步骤(5)：将每个分量看作为单个通道的脑电信号，并以原始信号的采样频率对其采样，然后合并分量构造新的信号矩阵，矩阵维数为k×n，k表示IMF分量个数，n表示采样点数目。以C3,C4通道为例，矩阵表示为：

记X₁和X₂分别为按式(16)构造的两类运动想象任务下的信号矩阵，构造空间滤波器，首先求信号矩阵协方差，并归一化，得到：

其中，trace(·)为矩阵的迹，然后求出混合空间的协方差矩阵R,并进行特征值分解：

其中，

为两类任务实验的平均协方差矩阵，U为矩阵R的特征向量，Λ为对应的特征值矩阵。进一步可求出白化矩阵P：

接着对R₁与R₂做白化处理以及主成分分解：

则矩阵S₁与矩阵S₂的特征向量相等，即B₁＝B₂，且Λ₁+Λ₂＝I，其中，I为单位矩阵。两类矩阵的特征值之和为1，因此当一类的特征值最大时另一类最小，从而可以实现分类的作用。则构造空间滤波器：

W＝B^TP (24)

将上述由IMF分量构造的运动想象信号矩阵经W滤波，得到特征矩阵，接着取其前q个和后q个特征值构成信号Z_P(P＝1,…,2q)，求取特征向量：

得到表示空域的特征，记为F₃；

步骤(6)：分别求出时域瞬时能量均值特征F₁，多尺度熵特征F₂，以及重构信号的CSP特征F₃，最后得到多域的结合的特征向量F＝{F₁,F₂,F₃}，并进行归一化。将训练集特征输入SVM分类器训练分类模型后，对测试集特征进行分类得到测试结果。

为验证本发明在脑电信号特征提取及分类研究中的真实分类性能，首先对单次实验中C3和C4通道运动想象脑电信号进行多元变分模态分解，分解结果如图2所示，其中包括原始波形和分解得到的4阶IMF分量。在图2中可以看出C3,C4通道的原始信号以及MVMD分解得到的4阶IMF分量信号，每一个IMF分量中都包含了与脑电信号相关的信息，通过对这些分量的处理得到需要的特征，具体处理过程在上述步骤中已有详细的介绍。为了验证所提方法的有效性，就时域特征向量F₁，分别做出所选数据集的C3通道测试数据中想象执行的两类各70次的特征向量，特征对比如图3(a)所示，对特征值统计如图3(b)所示。在图3(a)中，红色星号和蓝色方块分别表示测试数据中想象左手和右手任务时的特征值。根据图中结果，第一类的特征值较大于第二类；图3(b)的统计结果显示了测试数据集的140次实验中，特征值较小时，第一类运动想象任务的个数少于第二类，相反，特征值较大时，第一类运动想象任务的个数多于第二类。因此可以起到分类的作用，但是由于缺少多域特征的结果，在图3中可以看出，个别实验中没有明显的区分效果，因此多角度特征提取是有必要的。

就所选数据集的单一特征和多域的结合特征分别进行分类实验，得到的分类准确率如表1所示。可以看出，利用不同的单一特征向量区分左右手想象运动任务时，时频特征F1表现最好，识别准确率达到83.93％。但是利用本发明所提的多域特征F，可以获得89.64％的准确率，与单一的时频特征，非线性动力学特征和空域特征相比，分别提高了5.71、13.57和10.35个百分点。这是因为与单一特征相比较，多域的结合特征实现了信号不同特征之间的互补，可以获得更全面的信号表达，进而可以获得更好的识别效果。

表1不同特征集分类结果

最后，在表2中给出了其它的多特征结合的运动想象脑电信号分类方法与本发明所提方法的对比，所有方法均是在同一数据集，采用相同的通道数，通过测试分类精度来衡量本发明在脑电信号分类领域的可行性。

表2不同方法分类准确率的比较

可以看出本发明所提方法与HCHT方法中基于希尔伯特-黄变换和CSP的多特征融合提取方法相比，准确率高6.78个百分点，这是因为本发明方法在分解脑电信号时采用MVMD算法，解决了EMD中模态混叠问题，同时避免了通道间固有模态函数个数以及频率不匹配的问题。Linyan Wu的方法结合了MEMD和CSP两种提取特征的方法，就计算过程而言，该方法在MEMD分解之后选择的有效IMF分量为IMF分量间的组合形式，虽然在一定程度上提高了准确率，但是增加了计算的复杂性，且本发明算法识别率与之相比，提高了5.47个百分点。杨默涵和Chen S采用多种非线性特征结合的方法，并取得不错的效果，但都忽略了每一类运动想象脑电信号中重要的时域特征和空间分布成分，而本发明在得到特征向量F时，特征提取方法结合了时域特征，非线性动力学特征以及空域特征，对脑电信号的分析更为全面，因此识别率与之相比，分别提高5.41和3.93个百分点。从整体来看，应用本发明所提方法进行两类运动想象任务的分类，得到的准确率与以上所提方法相比均有所提高，对脑电信号的特征提取提供了新的思路。

Claims (2)

1.基于多元变分模态分解的脑电多域特征提取方法，其特征在于：包括以下步骤：

步骤(1)：多通道脑电信号采集并进行预处理；

步骤(2)：对多通道数据进行MVMD分解；具体为：

对于包含C个数据通道的输入数据X(t)，记为[x₁(t),x₂(t),…x_C(t)]，MVMD算法的具体步骤如下：

①、预先定义k个窄带IMF分量u_k(t)，使

其中，u_k(t)＝[u₁(t),u₂(t),…u_C(t)]；

利用希尔伯特变换的算法分别得到向量u_k(t)中每一个元素的解析表示，记为

以计算单边频谱，然后与指数项

相乘来调整其对应的中心频率ω_k(t)，使每个模态的频谱被调制到与它相对应的基频带上，通过谐波转换后

的梯度函数的L₂范数来估计出各模态u_k(t)的带宽；

由于单个的频率分量ω_k被用作为整个向量

的谐波混频，因此需要在多元振荡u_k(t)中找到在多个通道的共同频率分量ω_k，考虑

的所有通道，使每个通道所分解的IMF分量之和能够再现该输入信号，且最小化模态函数带宽之和，受约束的相关优化问题变为：

其中，u_k,c(t)表示相应的通道数c和向量u_k(t)中的元素，

表示对应的u_k,c(t)中每一个元素的解析表示，

表示与时间相关的偏导数；

③、对变分问题求解，构造增广的拉格朗日表示为：

其中，α为二次惩罚因子，其保障存在高斯噪声时信号的重构精度，λ_c(t)为拉格朗日乘子，使得约束条件保持严格性；

④、为了解决这个转变后的非约束性的变分问题，应用乘子交替方向法来更新

和

以寻找增广拉格朗日表达式的鞍点；具体做法为：

1)模态更新

模态更新的最小化问题相当于：

利用Parseval傅里叶变换，将其变换到频域，与VMD算法的模式更新方法类似，得到模式更新为：

2)中心频率更新

由于增广拉格朗日函数的后两项不依赖于ω_k，则中心频率更新所对应的优化问题可以简化为：

同样利用Parseval傅里叶变换，将该优化转换到频域，得到在傅里叶域中的等效问题为：

得到中心频率的更新方法：

通过更新关系自适应地分解信号的频带，得到k个窄带IMF分量；此外，由于MVMD方法能同时计算多个通道数据，保证了通道间频率的一致性，使信号的分析更趋于稳定；

步骤(3)：将式(1)的各分量u_k(t)进行希尔伯特谱分析，定义为u_k(t)与1/πt的卷积，即：

其中，μ为积分变量；

由u_k(t)作为实部，它的希尔伯特变换作为虚部，得到其解析信号和瞬时幅值为：

U_k(t)＝u_k(t)+j₀H[u_k(t)] (10)

其中，j₀为虚数单位；

然后计算每个通道的平均瞬时能量值为：

其中n为采样点数目，

为第i个采样点对应的瞬时能量值；取瞬时能量均值作为时域特征，记为F₁；

步骤(4)：

①根据采样点数对IMF分量作一维离散时间序列{u₁,u₂,…,u_i,…,u_n}，然后进行多尺度化处理来构造粗粒化形式的时间序列{y^(τ)}：

其中，

τ为尺度因子，序列长度为M＝int(n/τ)；

②构造一组m维的向量：

X_m(i)＝{y_i+b:0≤b≤m-1} (14)

③计算X(i)和X(j)间距离:

d[X(i),X(j)]＝max|y_(i+b)-y_(i+b)| (15)

其中，0≤b≤m-1；i,j＝1,2…M-m+1；i≠j，M为序列长度；

④计算d[X(i),X(j)]中小于阈值r的个数，也是模板匹配的个数，记为B^m(i)，在计算B^m(i)与距离总数的比值，记作：

⑤计算

的平均值C^m(r)：

⑥将维数增加至m+1，重复②～⑤步骤，计算

和C^m+1(r)；

得到时间序列长度为M时的样本熵值：

SampEn(m,r,M)＝-ln[C^m+1(r)/C^m(r)] (18)

重复以上步骤得到多尺度下的样本熵特征，记为F₂；

步骤(5)：将每个分量看作为单个通道的脑电信号，并以原始信号的采样频率对其采样，然后合并分量构造新的信号矩阵，矩阵维数为k×n，k表示IMF分量个数，n表示采样点数目；以C3,C4通道为例，矩阵表示为：

记X₁和X₂分别为按式(16)构造的两类运动想象任务下的信号矩阵，构造空间滤波器，首先求信号矩阵协方差，并归一化，得到：

其中，trace(·)为矩阵的迹，然后求出混合空间的协方差矩阵R,并进行特征值分解：

其中，

为两类任务实验的平均协方差矩阵，U为矩阵R的特征向量，Λ为对应的特征值矩阵；进一步可求出白化矩阵P：

接着对R₁与R₂做白化处理以及主成分分解：

则矩阵S₁与矩阵S₂的特征向量相等，即B₁＝B₂，且Λ₁+Λ₂＝I，其中，I为单位矩阵；两类矩阵的特征值之和为1，因此当一类的特征值最大时另一类最小，从而可以实现分类的作用；则构造空间滤波器：

W＝B^TP (24)

将上述由IMF分量构造的运动想象信号矩阵经W滤波，得到特征矩阵，接着取其前q个和后q个特征值构成信号Z_P(P＝1,…,2q)，求取特征向量：

得到表示空域的特征，记为F₃；

步骤(6)：分别求出时域瞬时能量均值特征F₁，多尺度熵特征F₂，以及重构信号的CSP特征F₃，最后得到多域的结合的特征向量F＝{F₁,F₂,F₃}，并进行归一化；将训练集特征输入SVM分类器训练分类模型后，对测试集特征进行分类得到测试结果。

2.根据权利要求1所述的基于多元变分模态分解的脑电多域特征提取方法，其特征在于：步骤一中多通道脑电信号采集并进行预处理，具体为：

采集受试者的多通道运动想象脑点信号，接着利用带通滤波器对脑电信号进行8-30Hz的低通滤波，用于脑电信号特征提取。

Images