CN108338787A - 一种多时段多元多尺度锁相值的相位特征提取方法 - Google Patents

Abstract

本发明涉及一种多时段多元多尺度锁相值(MMMPLV)的相位特征提取方法，首先对运动想象脑电信号进行预处理，并采用噪声辅助多元经验模态分解(NAMEMD)将多导脑电信号同时进行分解，得到多元多尺度固有模态函数(IMFs)；然后根据各尺度IMF的频率范围选取有效的IMF分量；接着分时段计算任意两导联同一尺度IMFs的锁相值(PLV)，将任意两导联各尺度分时段的PLV串联，构成多时段多元多尺度锁相值特征，并将其作为脑电信号的相位特征；将相位特征输入到分类器中进行分类，根据分类正确率选取NAMEMD和MMMPLV中参数的最优值，最终获得最优参数下的脑电相位特征。本发明提高了运动想象脑电信号的分类正确率。

Description

一种多时段多元多尺度锁相值的相位特征提取方法

技术领域

本发明属于脑电信号处理技术领域，尤其涉及一种多时段多元多尺度锁相值的相位特征 提取方法，具体涉及脑-机接口(Brain-Computer Interface,BCI)系统中对运动想象脑电信号 特征的提取方法，采用基于噪声辅助多元经验模态分解(Noise-AssistedMultivariate Empirical Mode Decomposition,NAMEMD)的多时段多元多尺度锁相值(Multi-period Multivariate Multi-scale Phase Locking Value,MMMPLV)方法对运动想象脑电信号进行相位特征提取。

背景技术

人们想象肢体运动时会激活与实际运动相似的脑区并产生运动想象皮层脑电信号(Motor Imagery Electrocorticography，MI-ECoG)。根据皮层可塑性理论，运动想象可以促进受损伤的 运动传导通路修复或重建，所以在脑机接口领域，研究人员试图通过识别不同的运动想象任 务实现大脑与外界设备的通讯，来帮助运动功能障碍患者与外界进行沟通和交流。因此，如 何准确提取MI-ECoG的特征在医疗、康复以及BCI领域尤为重要。

人脑是一个具有复杂结构和功能的生物系统，不同的高级脑功能需要多个不同区域神经 系统之间进行不同层次的整合和协调来完成，而各神经网络的同步振荡是大脑信息整合和信 息处理的主要潜在机制。所以，脑电信号的同步化是研究不同脑区间信息交流的关键部分。 针对MI-ECoG信号的非线性以及非平稳性，由于相位同步(PhaseSynchronization，PS)法不受 信号幅值的影响，多被用于研究不同脑区和不同导联间的同步关系。锁相值(Phase Locking Value,PLV)方法是一种相位同步方法，它又称为平均相位相干(Mean Phase Coherence,MPC) 方法。PLV方法是通过对整个脑电信号时间序列求平均来获得两导信号之间的相位同步关系。 目前，PLV方法和基于经验模态分解(EmpiricalMode Decomposition,EMD)的PLV方法在 MI-ECoG信号的同步性研究中得到了一定的应用。

但是PLV方法存在三方面不足：第一，传统的PLV方法通过希尔伯特变换(HilbertTransform，HT)或小波变换(Wavelet Transform,WT)提取每导信号的相位信息，但是基于WT计算PLV的方法需要事先确定小波的中心频率和带宽，所以仅对某个频段的相位同步 敏感。基于HT计算PLV的方法在处理宽带信号时，要先对信号进行滤波，这会造成部分 信息的遗漏并且缺乏自适应性。第二，传统的PLV方法和基于EMD的PLV方法都只在某 一个频段或某一个尺度上通过平均整个脑电信号时间序列或本征模态函数(Intrinsic ModeFunctions,IMFs)序列来计算两导信号之间的相位同步关系，这忽略了MI-ECoG信号的多尺度特性以及时变特性。第三，EMD方法在分解信号时只能对一导信号进行分解，忽略了多 通道脑电信号之间的信息，并且EMD还会产生模态混叠和模态校准问题，这不利于多通道 MI-ECoG信号的研究，并且会导致模式分类正确率不高。

发明内容

针对现有PLV方法的不足，本发明提出了一种多时段多元多尺度锁相值的相位特征提取 方法，基于NAMEMD的MMMPLV的运动想象脑电信号特征提取。此方法不仅保留了多通道信号之间的有用信息，还从多个尺度上提取MI-ECoG信号的细节相位特征，精细地刻画了MI-ECoG的非线性、非平稳性、多尺度以及时变等特性，并且提高了脑电信号的分类正确率。

为实现上述目的，本发明采用如下的技术方案：

一种多时段多元多尺度锁相值的相位特征提取方法，包括以下步骤：

步骤1、对MI-ECoG脑电信号进行预处理；

步骤2、对步骤1得到的多导脑电信号同时进行NAMEMD分解，得到多元多尺度IMFs分量；

步骤3、接着选取符合MI-ECoG频率范围的IMFs分量，采用MMMPLV方法提取脑 电信号的相位特征，使用支持向量机(Support Vector Machine,SVM)对相位特征进行分类， 根据分类正确率对NAMEMD和MMMPLV中的参数进行优化，最终提取出最优参数下的 脑电相位特征。

MI-ECoG信号特征提取的生理学依据是：研究表明，当受试者在进行单侧肢体运动想象 时，对侧大脑皮层运动感觉区的μ节律(8-13Hz)和β节律(13-30Hz)的振幅明显降低，这种现 象称为事件相关去同步(Event-Related Desynchronizations,ERD)；而同侧大脑皮层运动感觉区 的μ节律和β节律的振幅明显增大，称为事件相关同步(Event-RelatedSynchronizations,ERS)。 ERD/ERS现象为运动想象脑电信号特征提取的研究提供了理论基础。

作为优选，步骤1具体为：

对信号进行8-30Hz带通滤波，得到信号X(t)＝[X₁(t),X₂(t),Λ,X_N(t)]^T∈R^N×K，其中 N和K分别表示导联数和采样点数。

作为优选，步骤2具体为：

步骤2.1，将N导脑电信号和M导不相关的并且和原信号等长的高斯白噪声信号组合， 构成复合信号X¹(t)＝[X₁(t),X₂(t),...,X_N+M(t)]^T∈R^(N+M)×K，其中，噪声幅值为A_m；

步骤2.2，对(N+M)元信号X¹(t)进行多元经验模态分解，参数方向向量的个数为d，信 号X¹(t)被分解成一系列分量和余量r(t)的加和形式，即：

其中，cj(t)为(N+M)元信号的(N+M)元IMFs分量r(t)表示 (N+M)个余量q表示(N+M)元IMFs分量的层数，各层IMF包含 了原信号不同频率尺度特征的信息；

步骤2.3，从(N+M)元IMFs中删除M导噪声通道对应的IMFs分量，保留步骤1中脑 电信号X(t)的N元IMFs分量j＝1,2,Λq。

作为优选，步骤3具体包括：

步骤3.1，根据各尺度IMF的频率范围，选取N导信号的前L层IMFs分量 (j＝1,2,ΛL，L≤q)；

步骤3.2，对N导信号前L层IMFs分量分别构造解析信号并求得其各尺度IMF的瞬时相位公式如下：

其中，i为虚数单位，l表示第l层IMF分量，l＝1,2,Λ,L，X_s表示第s导脑电信号， s＝1,2,Λ,N，为对应的希尔伯特变换为解析信号的幅值；

步骤3.3，计算a,b两导信号同一尺度的IMF的瞬时相位差，公式如下：

其中，a,b＝1,2,Λ,N，a≠b；

步骤3.4，计算给定时间窗下a,b两导信号间第l层IMF的第n段的锁相值；

其中，t_L表示选取的时间窗口长度，<·>表示对时间段t_L取平均，在给定时间窗下，各尺度IMF无重叠地被分为w个时段，w＝T_L/t_L，T_L为总的采样时间，n＝1,2,Λ,w，PLV取值范围为[0,1]；

步骤3.5，构建特征向量F；

将a,b两导信号间的同一尺度IMF的w个时段的锁相值串行融合，构成多时段锁相值， 即MPLV：

再将a,b两导信号间的各尺度IMF的MPLV串行融合形成多时段多尺度锁相值，即MMPLV：

最后，将任意两导联的MMPLV串行融合，构成多时段多元多尺度锁相值，即MMMPLV：

其中，表示特征维数，公式如下：

其中，C表示组合运算；

步骤3.6，根据分类结果优化参数；

对两类运动想象任务的脑电数据分别提取相位特征并使用分类器进行分类，选取分类 正确率最高时的参数N，t_L，A_m和d，最终获得最优参数下的相位特征向量。

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

本发明针对传统PLV方法和基于EMD的PLV方法在提取脑电特征时忽略了MI-ECoG信号的多通道、多尺度以及时变特性的不足，采用NAMEMD将多通道MI-ECoG信号同时进 行分解，保留了MI-ECoG信号多个通道之间的有用信息，并且从多个尺度上分时段提取 MI-ECoG信号的细节相位特征，更加适合MI-ECoG信号的本质分析。相比于传统PLV方法 和基于EMD的PLV方法，本发明提取了MI-ECoG信号更加广泛的相位特征信息，有助于分 析MI-ECoG信号多个通道之间的相互关系。本发明使用支持向量机对提取的脑电特征进行分 类，采用9折交叉验证法，最高平均分类正确率达到了83.07％，高于传统的PLV方法和其他 多种基于经验模态分解的PLV方法的分类结果。本发明不仅可以获取脑电信号多个导联之间 的信息，还从多个尺度上提取细节相位特征，更加适合多通道、多尺度、非线性、非平稳性 的脑电信号的分析，相比于传统的PLV方法，本发明进一步提高了脑电信号的分类正确率。

附图说明

图1为本发明所涉及方法的流程框图；

图2为电极位置分布图；

图3为实验时序图；

具体实施方式

如图1所示，本发明提供一种多时段多元多尺度锁相值的相位特征提取方法，包括以下 步骤：

步骤1，信号预处理；

首先，对信号进行8-30Hz带通滤波，得到信号X(t)＝[X₁(t),X₂(t),Λ,X_N(t)]^T∈R^{N ×K}， 其中N和K分别表示导联数和采样点数；

步骤2，将对步骤1得到的多导脑电信号同时进行噪声辅助多元经验模态分解；

步骤2.1，将N导脑电信号和M导不相关的并且和原信号等长的高斯白噪声信号组合， 构成复合信号X¹(t)＝[X₁(t),X₂(t),...,X_N+M(t)]^T∈R^(N+M)×K，其中，噪声幅值为A_m；

步骤2.2，对(N+M)元信号X¹(t)进行多元经验模态分解，参数方向向量的个数为d，信 号X¹(t)被分解成一系列分量和余量r(t)的加和形式，即：

其中，c^j(t)为(N+M)元信号的(N+M)元IMFs分量r(t)表示 (N+M)个余量q表示(N+M)元IMFs分量的层数，各层IMF包含 了原信号不同频率尺度特征的信息；

步骤2.3，从(N+M)元IMFs中删除M导噪声通道对应的IMFs分量，保留步骤1中脑 电信号X(t)的N元IMFs分量j＝1,2,Λq；

步骤3，采用多时段多元多尺度锁相值方法提取脑电信号的相位特征；

步骤3.1，根据各尺度IMF的频率范围，选取N导信号的前L层IMFs分量 (j＝1,2,ΛL，L≤q)；

步骤3.2，对N导信号前L层IMFs分量分别构造解析信号并求得其各尺度IMF的瞬时相位公式如下：

其中，i为虚数单位，l表示第l层IMF分量，l＝1,2,Λ,L，X_s表示第s导脑电信号， s＝1,2,Λ,N，为对应的希尔伯特变换为解析信号的幅值；

步骤3.3，计算a,b两导信号同一尺度的IMF的瞬时相位差，公式如下：

其中，a,b＝1,2,Λ,N，a≠b；

步骤3.4，计算给定时间窗下a,b两导信号间第l层IMF的第n段的锁相值；

其中，t_L表示选取的时间窗口长度，<·>表示对时间段t_L取平均，在给定时间窗下，各尺度 IMF无重叠地被分为w个时段，w＝T_L/t_L，TL为总的采样时间，n＝1,2,Λ,w，PLV取值范 围为[0,1]；

步骤3.5，构建特征向量F；

将a,b两导信号间的同一尺度IMF的w个时段的锁相值串行融合，构成多时段锁相值， 即MPLV：

再将a,b两导信号间的各尺度IMF的MPLV串行融合形成多时段多尺度锁相值，即MMPLV：

最后，将任意两导联的MMPLV串行融合，构成多时段多元多尺度锁相值，即MMMPLV：

其中，表示特征维数，公式如下：

其中，C表示组合运算；

步骤3.6，根据分类结果优化参数；

对两类运动想象任务的脑电数据分别提取相位特征并使用分类器进行分类，选取分类 正确率最高时的参数N，t_L，A_m和d，获得最优参数下的相位特征向量。

实施例1：

本发明中的实验是在Matlab2011a的仿真环境下进行的。

本发明采用的MI-ECoG数据来源于BCI Competition III的Data Set I数据集。该数据集 记录想象左手小指与想象舌头两种任务，数据通过植入受试者大脑右侧运动皮层的一个8×8 (64导)的铂电极阵列获得，如图2所示，信号采样频率为1000Hz，并经过0.016-300Hz带 通滤波。数据集包含378组实验，其中训练集278组，想象左手小指和舌头各139组，测试 集100组，两种想象任务各50组。实验时序图如图3所示，每次实验持续7秒，每两次实验 均间隔几分钟。t＝0s时，显示器上出现“十”字光标，t＝1～5s时，显示器上的“十”字光标 由舌头或左手小指图片代替，此时，受试者根据提示的图片想象舌头或左手小指运动。t＝5～7s 时，受试者休息。然后进行下一次实验。

(1)数据预处理及导联选取

首先，为避免视觉诱发电位的影响，实验截取1.5～4.5s的3s数据段做离线分析。然后， 本文对Data set I中的训练集和测试集进行下采样至100Hz，即将3000采样点数降至300，提 高数据的处理速度。接着，对信号进行8-30Hz频段带通滤波，以获取ERD生理现象较为明 显的数据段。

当受试者进行想象左手小指或舌头的任务时，不同电极处ERD现象的强度不同。所以， 为了提高识别的精度和速度，获得ERD现象较为明显的导联，本发明方法采用相对距离准则 优选导联，对64导脑电信号优选出的前10导联为38、30、39、29、31、46、40、37、12、 21导联。优选导联后的训练集数据为X_train∈R^10×300×278，测试集数据为X_test∈R^10×300×100。

(2)将优选导联的脑电信号进行噪声辅助多元经验模态分解

首先，将N(3≤N≤10)导脑电信号和2导不相关的噪声幅值为A_m的并且和原信号等长 的高斯白噪声信号组合，构成复合信号X¹(t)∈R^(N+2)×300；然后将复合信号同时经过多元经验模态分解，并删除2导噪声通道对应的IMFs分量，保留脑电信号的N元q层IMFs分 量j＝1,2,Λq；此时，脑电信号的训练集各个导联的数据 为X_train∈R^q×300×278，测试集各个导联的数据为X_test∈R^q×300×100。

(3)采用多时段多元多尺度锁相值方法提取脑电信号的相位特征

根据各尺度IMF分量的频率范围可知，前五层IMF分量在8-30Hz内，所以选取前5层IMFs分量进行实验；由公式(2)、(3)分别对N导信号前5层IMFs分量进行Hilbert 变换得到N导信号前5层IMFs分量的瞬时相位；然后根据公式(4)计算任意两导信号同 一尺度IMF分量的瞬时相位差；选取t_L为0.5s，1s，1.5s，和3s的时间窗长，根据公式(5) 计算给定时间窗下两导信号间同层IMF分量的锁相值，其中，尺寸为3秒的时间窗是指整 个采样时间段；根据公式(6)-(8)得到多时段多元多尺度锁相值MMMPLV特征；将想 象左手小指和舌头运动的两类脑电数据的相位特征输入到SVM(采用LibSVM软件包)中 进行分类，其中，SVM的核函数选取径向基函数，对SVM参数c和g进行遍历，根据9 折交叉验证的最高平均分类正确率选取N，t_L，A_m和d的参数值，最终获得最优参数下的 相位特征向量。由于数据集一共包含378次实验，实验次数可以被9整除，所以选择9折 交叉验证法得到最高平均分类正确率。实验中，当N＝5，t_L＝1s，A_m＝22dBW以及d＝512 时，平均分类正确率达到最高83.07％。此时，脑电信号的训练集相位特征为 F_train＝MMMPLV∈R^150×278，测试集相位特征为F_test＝MMMPLV∈R^{150 ×100}。

表1是当N＝5以及d＝512时，采用MMMPLV方法在不同时间窗和不同噪声大情况下得 到的平均分类正确率。

表1不同时间窗和不同噪声大小对平均分类正确率的影响

从表1中可以看出，t_L＝1s、A_m＝22dBW情况下，平均识别率达到最高83.07％，得 到了最好的分类结果。并且，分时段计算的MMMPLV平均识别率几乎均高于以整个采样 时间作为时间窗的MMMPLV的平均识别率,这证明了MMMPLV方法的有效性。

表2是在相同数据集Data Set I下，将本发明方法与传统的PLV方法和基于多种经验模 态分解的PLV方法进行对比的实验结果，实验将分时段计算的单尺度与多尺度相位特征进行 对比。其中，N＝5，t_L＝1s，采用SVM对特征进行分类，9折交叉验证平均识别率如表2所 示。

表2同一数据库多种PLV方法平均分类正确率结果对比

从表2中可知，几乎所有方法的多尺度PLV平均分类正确率都高于单尺度PLV的平均 分类正确率。本发明方法采用NAMEMD的MMMPLV方法，相比单通道信号分解方法EMD、OEMD和EEMD取得了更好的分类效果，并且NAMEMD进一步去除了MEMD中的模态 混叠现象，提升了分类正确率，证明了MMMPLV方法的有效性。

本发明方法最高平均分类正确率为83.07％，高于传统的PLV方法和其他基于多种经验 模态分解的PLV方法，取得了较好的分类性能。本发明方法不仅获取了多元信号各个通道 之间的信息，从多尺度上衡量了MI-ECoG的复杂性，并且考虑到了MI-ECoG的时变特性， 更精确地提取了相位特征。本发明从MI-ECoG信号的本质出发，针对传统的PLV和基于EMD方法的PLV的不足进行了改进，更加适合具有非线性、非平稳、多通道、多尺度等特 点的MI-ECoG信号的分析，并且进一步提升了分类效果，较好地区分了两类运动想象脑电 信号，证明了本发明方法的有效性。

Claims (4)

1.一种多时段多元多尺度锁相值的相位特征提取方法，其特征在于，包括以下步骤：

步骤1、对MI-ECoG脑电信号进行预处理；

步骤2、对步骤1得到的多导脑电信号同时进行NAMEMD分解，得到多元多尺度IMFs分量；

步骤3、接着选取符合MI-ECoG频率范围的IMFs分量，采用MMMPLV方法提取脑电信号的相位特征，使用支持向量机(Support Vector Machine,SVM)对相位特征进行分类，根据分类正确率对NAMEMD和MMMPLV中的参数进行优化，最终提取出最优参数下的脑电相位特征。

2.如权利要求1所述的多时段多元多尺度锁相值的相位特征提取方法，其特征在于，步骤1具体为：

对信号进行8-30Hz带通滤波，得到信号X(t)＝[X₁(t),X₂(t),Λ,X_N(t)]^T∈R^N×K，其中N和K分别表示导联数和采样点数。

3.如权利要求2所述的多时段多元多尺度锁相值的相位特征提取方法，其特征在于，步骤2具体为：

步骤2.1，将N导脑电信号和M导不相关的并且和原信号等长的高斯白噪声信号组合，构成复合信号X¹(t)＝[X₁(t),X₂(t),...,X_N+M(t)]^T∈R^(N+M)×K，其中，噪声幅值为A_m；

步骤2.2，对(N+M)元信号X¹(t)进行多元经验模态分解，参数方向向量的个数为d，信号X¹(t)被分解成一系列分量和余量r(t)的加和形式，即：

其中，c^j(t)为(N+M)元信号的(N+M)元IMFs分量r(t)表示(N+M)个余量q表示(N+M)元IMFs分量的层数，各层IMF包含了原信号不同频率尺度特征的信息；

步骤2.3，从(N+M)元IMFs中删除M导噪声通道对应的IMFs分量，保留步骤1中脑电信号X(t)的N元IMFs分量

4.如权利要求3所述的多时段多元多尺度锁相值的相位特征提取方法，其特征在于，步骤3具体包括：

步骤3.1，根据各尺度IMF的频率范围，选取N导信号的前L层IMFs分量 (j＝1,2,ΛL，L≤q)；

步骤3.2，对N导信号前L层IMFs分量分别构造解析信号并求得其各尺度IMF的瞬时相位公式如下：

其中，i为虚数单位，l表示第l层IMF分量，l＝1,2,Λ,L，X_s表示第s导脑电信号，s＝1,2,Λ,N，为对应的希尔伯特变换为解析信号的幅值；

步骤3.3，计算a,b两导信号同一尺度的IMF的瞬时相位差，公式如下：

其中，a,b＝1,2,Λ,N，a≠b；

步骤3.4，计算给定时间窗下a,b两导信号间第l层IMF的第n段的锁相值；

其中，t_L表示选取的时间窗口长度，<·>表示对时间段t_L取平均，在给定时间窗下，各尺度IMF无重叠地被分为w个时段，w＝T_L/t_L，T_L为总的采样时间，n＝1,2,Λ,w，PLV取值范围为[0,1]；

步骤3.5，构建特征向量F；

将a,b两导信号间的同一尺度IMF的w个时段的锁相值串行融合，构成多时段锁相值，即MPLV：

再将a,b两导信号间的各尺度IMF的MPLV串行融合形成多时段多尺度锁相值，即MMPLV：

最后，将任意两导联的MMPLV串行融合，构成多时段多元多尺度锁相值，即MMMPLV：

其中，表示特征维数，公式如下：

其中，C表示组合运算；

步骤3.6，根据分类结果优化参数；

对两类运动想象任务的脑电数据分别提取相位特征并使用分类器进行分类，选取分类正确率最高时的参数N，t_L，A_m和d，最终获得最优参数下的相位特征向量。