残障人士运动状态下脑电信号的高鲁棒解析方法与流程

本发明涉及人机交互技术领域，尤其涉及一种残障人士运动状态下脑电信号的高鲁棒解析方法。

背景技术：

脑机接口(brain-computer interface，BCI)是大脑与外部设备间建立的直接连接通路。利用大脑中的各种信号实现人脑与计算机或其他电子设备的通讯和控制，是一种全新的人机接口方式。当前获取脑信号的方式主要有脑电图EEG、功能磁共振成像fMRI和脑磁图MEG等多种形式。EEG信号具有时间分辨率高和可以无创获取等特点，使得利用大脑无时无刻不在产生的脑电波信号来解读大脑工作状态和思维活动成为当前BCI研究的热门和主流方向。

BCI设备还可以帮助残障人士。例如，美国研制了由肌肉电信号和脑电信号联合控制的义肢。在2012年11月，美国男子扎克·沃特利用脑控义肢，成功攀登了103层的芝加哥威利斯大厦，在这个过程中，沃特只要想着触动电机和链条，使假肢的脚踝和膝盖的动作保持一致，这个仿生假肢就会对腿筋部位的肌肉产生电脉冲，实现行走的功能。

目前BCI技术发展总体上还是处于初级的研究阶段，BCI系统的实用场景仍然处于实验室的理想条件下，距离真正实际应用还有相当长的路要走。首先，当前的脑电波解析技术还不够完善，把情绪与认知状态作为固定模式看待，缺少动态、系统的研究视点，还有待进一步研究发展。其次，国内外脑认知研究中都要求使用者静止不动并尽可能地屏蔽掉各种干扰信号，而在实际应用情况下这些要求都是不能满足的。第三，探索在工作人员运动的状态下检测和解读脑电信号，目前的国内外研究都没有涉及这一点。第四，脑认知从来就是一个复杂的动态过程，需要时序建模技术来分析处理，这也是目前的国内外研究没有深入的领域。此外，目前研究中大多数都只是分析单一动作所引起的脑电波变化特征，通常还会要求实验人员要处于心情平和的状态，从而排除情绪带来影响，这样通常不适用于实际应用的复杂环境和状态。所以要想将BCI研究进一步深入和实用化，就必须要攻克运动状态下情绪和脑认知状态的高鲁棒解析的问题。

技术实现要素：

本发明的目的在于克服突破在安静条件下研究脑认知活动的传统做法，针对运动协变脑电波，提出了一种残障人士运动状态下脑电信号的高鲁棒解析方法。

为达上述目的，本发明通过以下技术方案实现：

一种残障人士运动状态下脑电信号的高鲁棒解析方法，包括：

获取运动状态下的脑电样本信号，建立具有变化运动特征的信号集；获取残障人士运动状态下的待认知脑电信号；以及

对EEG信号进行基于自适应CQT的脑电信号timbre分析，包括：

首先，对脑电图EEG信号进行预处理，去除噪声与伪迹干扰；其次，根据脑电信号timbre的谐波规律，进行自适应差异的提取，自适应地找到各频带中心位置；最后，根据各频带中心位置，计算CQT频谱，提取各频带频谱特征参数构造特异性timbre；其中，将传统CQT公式修改为：

其中相邻谱线的间隔B_k由B_k＝f_k+1-f_k计算得到，各频带采样带宽N_k由N_k＝f_s/B_k计算得到；

使用基于希尔伯特-黄变换HHT的多尺度时频分析算法对脑电信号的谐波成分和精细结构进行分析，提取脑电信号时-频-空域特征：所述多尺度时频分析算法首先对脑电信号进行经验模态分解EMD分解，得到若干本征模态函数IMF，对这些IMF进行希尔伯特谱分析HSA，计算各IMF对应的瞬时频率，然后进行空间电极的选择，最后提取0-30Hz频段的IMF分量，基于这些IMF分量重建原信号；

基于分层递阶的可视化特征融合与降维算法对上述步骤获得的高维的特征进行降维处理：首先，对高维数据进行特征排序，并进行分层处理，之后对各层数据进行可视化特征融合，得到能够反映此层数据特性的特征，并将各层提取的特征组合在一起作为下一级处理的输入，重复上述过程，得到能够表示高维数据的低维特征。

进一步地，自适应差异的提取具体为：使用各频带能量最大值的频率位置来表示个体间的差异，计算方法如下：1)样本选取：从信号中随机取出一维，并取出长为L的一段来计算功率谱和频率序列；2)插值计算：使用样条插值法在频率序列基础上提高频率分辨率；3)频带能量提取：将原始EEG信号划分为Delta(0.5-4Hz)，Theta(4-8Hz)，Alpha(8-13Hz)，Beta(13-20Hz)，以及Gamma(30-50Hz)五个频带，分别计算这五个频带能量值；4)差异计算：求出各段频率能量最大值时频率的位置。

进一步地，空间电极选择为：针对特定被试者，首先选定一个最优准则J，计算所有电极的最优准则J，并且选择J最大的那个电极作为初始电极；任取一个其它的电极与初始电极组合在一起并且计算准则J，如果结果小于初始电极的J，则丢弃这个电极数据；如果大于初始电极的J，则取这个新的电极组合为初始电极；继续上一步，直到取完所有电极，此时得到的初始电极即为针对该被试者的最优的采集电极组合。

进一步地，设原始特征Feature为(x₁,x₂,…,x_i,…,x_n)，x_i表示第i维特征，n为特征维数，基于分层递阶的可视化特征融合与降维算法流程如下：

输入：原始特征Feature(x₁,x₂,…,x_i,…,x_n)

输出：特征融合降维后特征Re_fea

(1)F-score值计算：对每维特征x_i计算其F-Score值F_i；

(2)特征排序：对所有特征(x₁,x₂,…,x_i,…,x_n)按F-Score值的大小降序排列，得到排列后的特征(x₁',x'₂,…,x_i',…,x'_n)及F-Score值排列(F₁',F₂',…,F_i',…,F_n')，其中F′₁≥F′₂≥…≥F′_i≥…≥F′_n；

(3)分层：对所有特征的F-Score值求平均，得到F_Mean，求得第一个小于F_Mean的点F_first，即F₁'≥F₂'≥…≥F'_{F_first-1}≥F_Mean，F_Mean＜F'_{F_first}，则将(x₁',x'₂,…,x'_{F_first-1})作为第一层，计算F_first至n的F-Score的平均值F_Mean2，同理求得第一个小于F_mean2的点F_second，将(x'_{F_first},x'₂,…,x'_{F_second-1})作为第二层，(x'_{F_second},x'_{F_second+1},…,x'_n)作为第三层；

(4)可视化特征融合：对第一层数据构成雷达图，m＝2求取雷达图的三角形的重心特征，以完成可视化特征融合，对第二层数据构成雷达图，m＝3求取雷达图的四边形的重心特征，以完成可视化特征融合，对第二层数据构成雷达图，m＝4求取雷达图的四边形的重心特征，以完成可视化特征融合，将三层特征融合后的特征组合构造出初始特征融合集S’；

(5)对初始特征融合集S’，重复(1)，(2)，(3)步，得到降维后特征Re_fea。

本发明的有益效果是：本发明提出了一种残障人士运动状态下脑电信号的高鲁棒解析方法，包括获取运动状态下的脑电样本信号，建立具有变化运动特征的信号集；获取残障人士运动状态下的待认知脑电信号；以及对EEG信号进行基于自适应CQT的脑电信号timbre分析，使用基于希尔伯特-黄变换HHT的多尺度时频分析算法对脑电信号的谐波成分和精细结构进行分析，提取脑电信号时-频-空域特征：基于分层递阶的可视化特征融合与降维算法对上述高维的特征进行降维处理。运动协变脑电波样本数据和残障人士的脑电波数据通过以上步骤可构造出时域-频域-空域相配合的认知信号图谱，进而分类识别出残障人士的脑电波数据所对应的含义。本发明突破在安静条件下研究脑认知活动的传统做法，提出运动协变脑波的新概念，研究人体运动对脑电信号的影响规律，满足实用化产业需求；摒弃必须分离干扰信号与目标信号的传统思想，研究协变脑电的根本规律和直接解析方法。

附图说明

图1是基于CQT变换的脑电信号频带分析过程示意图；

图2是基于希尔伯特黄变换的脑电信号时频分析过程示意图；

图3是分层递阶特征融合与降维算法示意图。

具体实施方案

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

首先，本发明以从静止、步行到奔跑的一系列运动状况为实验条件，用平均法、PCA、ICA等脑电(ERP)提取算法分别实现运动状态下脑静息态时的脑电波信息的提取，根据提取结果，建立具有变化运动特征的信号集，如简单运动信号集，复杂运动信号集等。用信号集内的所有信号作为刺激源，通过实验得到对应的被试的ERP信号集合。计算ERP信号的潜伏期、波幅、波形变化等时域特征和频谱特征等变换域特征，使用统计学的方法和数据挖掘的方法，寻找ERP信号特征变化与运动信号特征变化的规律，以及引起ERP信号改变的最低刺激信号改变的阈值。从而得出运动状态下脑静息时脑波的基本规律，建立其电生理模型。

在音频信号中，不同的发声体由于材料、结构不同，发出声音的音质(timbre)也就不同，因此，可以通过音质的不同规律去分辨不同的声源。根据timbre的特点，即使在同一音高和同一声音强度的情况下，也能区分出是不同乐器或人发出的。同样，以脑电信号为代表的各种认知信号，也具有规律的谐波成分和精细结构，本发明借用音乐领域的音质概念(timbre)来对此进行表达。不同类型的认知信号其timbre各不相同但具有一定的规律性，不同人同一类信号的timbre也具有其个人特点。本发明根据认知信号timbre的个人特点和不同采集位置的特异性(如：频带分布特点，成分能量分布特点，谐波结构模式等)，提出基于自适应CQT的信号timbre分析方法，能够快速找到脑电信号timbre的个人特点，实现过程如附图1所示。

常数Q变换(CQT：Constant Q Transform)由Brown和Pluckette在1991年提出并成功应用于音乐信号处理领域。鉴于音阶频率不是等间隔分布，而是按指数规律分布的。因此用传统的等间隔频率分析方法(如离散傅里叶变换DFT)来估计音阶频率时，会产生非常大的估计误差。从滤波器组的观点来解释，DFT就是一个中心频率等间隔分布、带宽相同的滤波器组。而CQT却具有多分辨率与中心频率自适应的特点，通过计算音乐信号的CQT谱，可直接得到该信号在各音阶频率的频谱值，因此CQT在信号处理中正逐渐获得广泛应用。

有限长序列x(n)的CQT表达式如下：

式中N_k是计算第k条频率f_k的CQT变换时所对应的窗长度，是长度为N_k的窗函数，Q是CQT变换中的常数因子，k是序列CQT谱的频率下标。

本发明的基于自适应CQT的脑电信号timbre分析过程主要包含三个部分。首先，对各被试EEG信号进行预处理，去除噪声与伪迹干扰；其次，根据脑电信号timbre的谐波规律，自适应地找到各频带中心位置；最后，根据各频带中心位置，计算CQT频谱，提取各频带频谱特征参数构造特异性timbre。

然后开始自适应差异的提取：使用各频带能量最大值的频率位置来表示个体间的差异。计算方法如下：1)样本选取：从信号中随机取出一维，并取出长为L的一段来计算功率谱和频率序列；2)插值计算：使用样条插值法在频率序列基础上提高频率分辨率；3)频带能量提取：将原始EEG信号划分为Delta(0.5-4Hz)，Theta(4-8Hz)，Alpha(8-13Hz)，Beta(13-20Hz)，以及Gamma(30-50Hz)五个频带，分别计算这五个频带能量值；4)差异计算：求出各段频率能量最大值时频率的位置。

接下来进行中心频率的估计：不同于音阶频率指数分布的规律，脑电信号频率分布规律在自适应差异提取步骤中自适应地找到。为此，将传统CQT公式修改为：

其中相邻谱线的间隔(带宽)B_k可由B_k＝f_k+1-f_k计算得到，因此各频带采样带宽N_k就可以由计算得到。

这样针对认知信号之间的个体差异问题，提出的能够自适应找到信号差异的方法，可以比传统的等间隔频带划分方法更能揭示脑电信号频谱特点。通过CQT可以快速找到脑电信号timbre的个人特点后，就可以使用基于HHT的多尺度时频分析技术对脑电信号的谐波成分和精细结构进行分析，提取脑电信号时-频-空域特征。

希尔伯特-黄变换(Hilbert Huang Transform，HHT)由美国航空航天局的黄锷教授于1998年提出。HHT是一种基于经验的信号分析方法，可以自适应地对信号进行分解，得到非线性非平稳信号的具有物理意义的表示。它由两部分构成，经验模态分解(Empirical Mode Decomposition，EMD)和希尔伯特谱分析(Hilbert Spectrum Analysis，HSA)。对于一个非线性、非平稳的信号，如脑电信号，直接进行希尔伯特变换得到的结果在很大程度上失去了原有的物理意义。EMD可以将信号分解成为有限并少量的具有不同特征尺度的本征模态函数(Intrinsic Mode Function，IMF)，它是从信号本身的特征时间尺度出发对信号进行分解，没有固定的先验基，是一种自适应的信号分解方法，因此得到的IMF分量一般具有明显的物理意义，这些IMF反映了信号内部固有的波动性，在每个周期上仅包含一个波动模态。对IMF进行HSA，可以得到各IMF的瞬时频率和瞬时振幅具有清晰的物理意义，能够表征信号的局部特征，继而得到的希尔伯特谱能够准确地反映出该物理过程中能量在各种尺度的时频空间上的分布规律。正好适应了我们对信号分析时在不同频率范围所需要不同的分辨率这一基本要求。

目前主流的EEG信号分析方法有两种：一种是传统时频特征组合法，将时域均值、频域功率谱组合作为特征矢量，主要是利用多种类别信息提供更多的特征，但较多的特征使得建立分类器的模型比较复杂，不利于实际系统中的应用；另一种是小波变化系数法，依据先验知识、抽取感兴趣频段的小波系数作为特征，但选择不同的小波对分类结果有一定的影响。这些脑电信号预处理方法存在一定的缺陷，即忽略了空间电极的选择而导致的大数据计算量，同时还有依据经验选择的电极无法根据不同被试者自身脑电信号的特点动态变化。针对这些不足，本发明提出了一种基于多域的脑电信号处理算法，不仅在时频域对脑电信号进行适当的预处理操作，而且进一步把预处理应用到脑电信号的空间分布中，有效克服传统预处理方法中系统数据计算量大的缺点。算法主要分为三个步骤，如附图2所示。

在对脑电信号的时频分析阶段，使用的是前面提到的HHT。有研究表明，脑电信号一般存在于0-30Hz的低频段，所以首先对脑电信号进行EMD分解，得到若干IMF，对这些IMF进行HSA，计算各IMF对应的瞬时频率，然后提取0-30Hz频段的IMF分量，最后基于这些IMF分量重建原信号。

A.经验模态分解

首先，对脑电信号进行EMD分解，其核心过程是通过“筛选”对信号进行分解的过程，对于给定的信号x(t)，具体的处理过程如下所示：

1)分别用光滑的三次样条函数连接信号x(t)的所有极大值点组成上包络x_up(t)，以及所有极小值点组成下包络x_low(t)，所有的数据都应该被上下包络包含在内；

2)对数据的上下包络求取均值，并从原始数据中减去进而得到差值函数d₁(t)，公式描述如下：

d₁(t)＝x(t)-(x_up(t)+x_low(t))/2 (3)

3)将d₁(t)作为新的x(t)，重复执行步骤1)和2)，直到d(t)满足IMF的两个条件，即：1)对于整个时间序列，数据中的极值点数目与过零率数相等或者最多相差一个；2)在任一个数据点位置，由局部极大值点确定的上包络以及由局部极小值点确定的下包络，它们的均值为0。此时记d₁(t)为c₁(t)，c₁(t)便是从原始信号中提取的第一个IMF，它包含信号的最高频成分；

4)计算差值信号r₁(t)，计算公式如下所示：

r₁(t)＝x(t)-c₁(t) (4)

将r₁(t)视为新的初始信号x(t)，重复执行1)-3)的筛选过程，直到第N阶的差值信号成为单调函数，不能再筛选出IMF分量。

5)最后，获得了一系列固有模态分量，这些简单而有着良好特性的分量就组成了原始信号x(t)的基函数，x(t)的扩展表达式如下：

其中，N表示IMF的数量，r_N(t)是最后的差值信号，c_j(t)即是求得的IMF。

B.希尔伯特谱分析

通过HSA计算每个IMF的瞬时频率，来判断每个IMF所处的频带。每个IMF被表示为x_imf(t)，它的希尔伯特变换被表示为y_imf(t)，与希尔伯特变换相关的分析函数z_imf(t)如下所示。

其中，

其中，a_imf(t)表示瞬时幅值函数，θ_imf(t)表示瞬时相位函数，瞬时频率可以按下式给出

C.空间电极选择

针对特定被试者，首先选定一个最优准则J。计算所有电极的最优准则J，并且选择J最大的那个电极作为初始电极；任取一个其它的电极与初始电极组合在一起并且计算准则J，如果结果小于初始电极的J，则丢弃这个电极数据；如果大于初始电极的J，则取这个新的电极组合为初始电极。继续上一步，直到取完所有电极。此时得到的初始电极即为针对该被试者的最优的采集电极组合。对于不同的被试者，可以得到不同的最优采集电极组合。

D.信号重建

如式(10)所示，通过将处于0-30Hz频带内的IMF进行叠加，便可得到重构的脑电信号。

其中，x^*(t)表示重构后的脑电信号，M表示处于0-30Hz频带内的IMF个数。

为了分析以上提取的特征，本发明采用基于分层递阶的可视化特征融合与降维算法对特征进行分析，并最终建立模型。在完成了脑电波信号的特征提取之后，得到一组高维的数据。本发明将可视化技术作为将高维数据降维的处理方法，并采用多层信息融合的思想，对特征空间按照特征重要性进行分层，分段，并通过递阶的方式运用雷达图挖掘特征间的特性，对数据化繁为简，达到信息融合与降维的目的，最终完成建模。

雷达图又可以称为蜘蛛图或者星点图，雷达图的图形信息称为图形特征。雷达图的图形特征包括面积特征、重心特征、相邻幅值比、位置特征、分区面积比等。其中重心特征能够较好的反应雷达图间的各个维之间的内部关系，是重要的可视化特征。

一个M维的样本数据构成的雷达图r₁,r₂,…,r_i,…,r_M，任意连续相邻的m维变量构成的m多边形(r_i,r_i+1,…,r_i+m-1)，由几何代数能够推出，此m多边形的重心特征为

其中，w为雷达图中相邻两维特征所构成的夹角，w＝2π/M，abs_im表示雷达图中(r_i,r_i+1,…,r_i+m-1)构成的m多边形的重心的幅值大小，angle_i表示雷达图中(r_i,r_i+1,…,r_i+m-1)构成m多边形的重心的角度。从上述重心特征的计算公式能够看出，雷达图重心特征的提取方法其实是一种无监督的非线性映射。

F-Score原本是一种衡量单个特征对两类样本分辨能力的评价指标，由于实际中多分类问题更加普遍，因而将F-Score进行扩展至多分类问题上。F-Score的定义为：给定训练样本集x_k∈R^m，k＝1,2,…,n，l(l≥2)为样本类别数，n_j为第j类的样本个数，其中j＝1,2,…,l。则训练样本第i个特征的F-Score定义为：

其中表示的为第i个特征在整个训练集上的平均值，表示的为第i个特征在第j类数据集上的平均值；表示的是第j类的第k个样本点的第i维特征的特征值。

基于分层递阶的可视化特征融合与降维算法的主要思想是首先对高维数据进行特征排序，并进行分层处理，之后对各层数据进行特征融合，得到能够反映此层数据特性的特征，并将各层提取的特征组合在一起作为下一级处理的输入，重复上述过程，得到能够表示高维数据的低维特征。算法流程如附图3所示。

设原始特征Feature为(x₁,x₂,…,x_i,…,x_n)，x_i表示第i维特征，n为特征维数，则基于分层递阶模型的可视化特征融合与降维方法的算法如下：

输入：原始特征Feature(x₁,x₂,…,x_i,…,x_n)

输出：特征融合降维后特征Re_fea

(1)F-score值计算：对每维特征x_i按式(12)计算其F-Score值F_i；

(2)特征排序：对所有特征(x₁,x₂,…,x_i,…,x_n)按F-Score值的大小降序排列，得到排列后的特征(x₁',x'₂,…,x_i',…,x'_n)及F-Score值排列(F₁',F₂',…,F_i',…,F_n')，其中F′₁≥F′₂≥…≥F′_i≥…≥F′_n；

(3)分层：对所有特征的F-Score值求平均，得到F_Mean，求得第一个小于F_Mean的点F_first，即F₁'≥F₂'≥…≥F'_{F_first-1}≥F_Mean，F_Mean＜F'_{F_first}，则将(x₁',x'₂,…,x'_{F_first-1})作为第一层，计算F_first至n的F-Score的平均值F_Mean2，同理求得第一个小于F_mean2的点F_second，将(x'_{F_first},x'₂,…,x'_{F_second-1})作为第二层，(x'_{F_sec}ond,x'_{F_second+1},…,x'_n)作为第三层；

(4)可视化特征融合：对第一层数据构成雷达图，并按式(11)，m＝2求取雷达图的三角形的重心特征，以完成可视化特征融合，对第二层数据构成雷达图，并按式(11)，m＝3求取雷达图的四边形的重心特征，以完成可视化特征融合，对第二层数据构成雷达图，并按式(11)，m＝4求取雷达图的四边形的重心特征，以完成可视化特征融合，将三层特征融合后的特征组合构造出初始特征融合集S’；

(5)对初始特征融合集S’，重复(1)，(2)，(3)步，得到降维后特征Re_fea。

运动协变脑电波样本数据和残障人士的脑电波数据通过以上步骤可构造出时域-频域-空域相配合的认知信号图谱，进而分类识别出残障人士的脑电波数据所对应的含义。本发明突破在安静条件下研究脑认知活动的传统做法，提出运动协变脑波的新概念，研究人体运动对脑电信号的影响规律，满足实用化产业需求；摒弃必须分离干扰信号与目标信号的传统思想，研究协变脑电的根本规律和直接解析方法。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。