CN105559777B - 基于小波包和lstm型rnn神经网络的脑电识别方法 - Google Patents

Abstract

基于小波包和LSTM型RNN神经网络的脑电识别方法，属于生物信息技术领域。具体包括：首先，利用平均功率谱法分析脑电信号的时域特性，确定有效的时域范围。其次，对脑电信号进行小波包分解，基于改进的距离准则获取小波包最优子空间，并选取小波包最优子空间的小波包系数构成脑电特征。最后，以LSTM型RNN神经网络作为分类器对脑电信号进行分类。本发明能够自适应地选取小波包最优子空间，并自动确定最佳频段的小波包系数构成脑电时频特征，其保留了脑电信号的时序信息，而LSTM型RNN神经网络能够充分利用基于小波包变换提取的脑电特征的时序信息，从而提高脑电信号的分类准确率，为脑电信号的识别提供一种新思路。

Description

基于小波包和LSTM型RNN神经网络的脑电识别方法

技术领域

本发明涉及脑机接口(brain-computer interface,BCI)系统中脑电信号(Electroencephalograph,EEG)的识别方法，具体涉及利用小波包和LSTM(Long-ShortTerm Memory,LSTM)型RNN(Recurrent Neural Network,RNN)神经网络对运动想象脑电信号(Motor Imagery EEG,MI-EEG)进行识别的方法。

背景技术

近年来随着全球老龄化问题日趋严重，老年人中又多数伴有中风及偏瘫症状，同时因交通事故而造成脊髓受损的人数也越来越多。医学理论和临床医学证明，康复训练对这类患者肢体运动功能的恢复和提高起到非常重要的作用。运动功能康复机器人的出现可以减轻康复医师的工作量，但是缺少患者的主动参与和反馈，治疗效果有待于提高。

BCI技术提供了一个良好的解决该问题的方案。BCI系统通常由五个部分组成：即信号采集系统、信号预处理、特征提取、分类识别和控制外部设备，如图1所示，其关键技术是脑电信号的识别，但由于脑电信号微弱、识别率低、自适应能力差等问题，BCI技术在康复机器人领域的应用还没有得到广泛利用，脑电信号的识别是目前的研究热点，同时也是研究难点。

临床医学研究表明，根据频带特点，人类大脑可以划分成四个脑电节律。人类大脑能够控制人体对侧的运动，比如人体右侧肢体运动(如右手运动)大脑的左侧运动区域的Mu节律(8～12Hz)和Beta节律(14～30Hz)能量明显减小，而同侧的运动感觉区Mu节律和Beta节律能量增大，这种电生理现象也被称为事件相关去同步(Event-RelatedDesynchronization，ERD)和事件相关同步(Event-Related Synchronization，ERS)。脑电信号的ERD/ERS现象是实现BCI技术的神经生理学基础。

脑电信号是一种典型的时间序列信号，而提取的时频脑电特征依然保持其时序信息，现有的脑电分类器，如BP神经网络、支持向量机(Support Vector Machine，SVM)、线性判别式分析(Linear DiscriminantAnalysis，LDA)等没有充分利用脑电信号的时序特征，因此分类效果并不是很理想。

发明内容

针对脑电信号识别率低，脑电特征的时序信息没有被充分利用等问题，提出一种基于小波包和LSTM型RNN神经网络相结合的脑电信号识别方法，记为WPT-LSTM。该方法能够充分利用脑电信号和脑电特征的时序信息，从而提高了脑电信号的分类准确率，为脑电信号的识别提供了一种思路。

本发明为解决所述技术问题，采用如下的技术方案：

首先，利用平均功率谱法分析脑电信号的时域特性，确定有效的时域范围。其次，先利用小波包对脑电信号进行分解，然后采用改进的距离准则获取小波包最优子空间，并选取小波包最优子空间的小波包系数构成脑电特征。最终，利用LSTM神经网络作为分类器对脑电信号进行识别分类。

本发明的具体步骤如下：

步骤1，信号采集及预处理。

首先按照国际标准10-20系统来安放电极，通过脑电采集设备采集n导脑电信号，并将采集到的信号通过有限脉冲响应(Finite Impulse Response,FIR)滤波器进行8-30Hz带通滤波，滤波后信号为其中N为样本点总数，n为脑电导联数目，m为采样点个数，为第i导联第j个采样点滤波后的脑电信号，t＝{1,2,…,N}，i＝1,2,…,n，j＝{1,2,…,m}。

步骤2，平均功率谱法分析脑电时域特性。

EEG信号通过电极帽采集并以电压幅值形式保存，其瞬时功率计算如下：

式中，P_i ^j(t)表示第t个样本第i导联第j个采样点脑电信号的瞬时功率。

假设P_i ^j为N次实验第i导联第j个采样点脑电信号的平均功率，则P_i ^j表示为：

根据式(2)分别计算出每个导联脑电信号的平均功率，并选取区别明显时间段的脑电信号进行特征提取。

步骤3，对步骤2中选取的脑电信号进行小波包分解，然后利用改进的距离准则获取小波包最优子空间，并选取最优子空间的小波包系数构成脑电特征。

步骤3.1，对脑电信号进行小波包分解。

信号在第j层第k点处的小波包分解系数递推公式为：

其中，n为滤波器系数序号，i∈{0,1,2…2^j-1}为小波包子空间序号，j∈{0,1,…}为分解尺度，k∈{1,2,…}为位置指标。g(n)和h(n)为一对正交镜像滤波器。经过小波包分解，原信号按频带划分为若干小波包子空间，各子空间对应的频段为：式中，f_s为信号采样率。

步骤3.2，基于改进的距离准则获取小波包子空间评价值。

设c个类别ω₁,ω₂…,ω_c的样本矢量集为{x^(p,b),p＝1,2,...,c,b＝1,2,...,N_p}，其中x^(p,b)＝(x₁ ^(p,b),x₂ ^(p,b),...,x_m ^(p,b))为ω_p类中的第b个m维样本矢量，为x^(p,b)的第l个分量。N_p为ω_p类中的特征矢量的个数。

首先，计算平均类内距离S_w如下：

式中，表示ω_p类中μ^(p)的第l个分量的均值，l＝1,2,...,m。

其次，设ω_p类样本的平均值矢量为μ^(p)，对c个类别的平均类间距离S_b定义为：

式中，μ_l为样本总体均值矢量的第l个分量。

同时引入权重因子α∈[0,1]表征对类内距离与类间距离的关注程度，得到改进的距离准则如下：

J_Al＝α*S_b-(1-α)*S_w (7)

利用该准则求得的评价值J_Al越大，即表示对应的可分性越好。

步骤3.3，选取最优子空间的小波包系数构成脑电特征。

按式(8)计算第j层第k个小波包子空间对所有样本小波包评价的平均值J_A，以表征对应的特征频段对分类的贡献度。

将小波包子空间对应的评价值按从大到小重新排列，即J^* _A1～J^* _As，即为信号可分性从大到小的排列。将评价值最高的前四个小波包子空间系数抽取出来作为脑电特征。

这里需要说明的是，由于EEG信号的个体差异性，不同实验者的最优子空间可能不同，该方法能够根据不同个体的实际信号自动选择合适的最优子空间。对于同一个个体来说，小波包最优子空间的选取只在训练阶段实施，在测试与使用阶段直接利用选取好的最优子空间的小波包系数构成脑电特征，这样大幅提高计算效率以满足脑电信号在线识别的要求。

步骤4，利用LSTM型RNN神经网络进行脑电特征分类。

RNN神经网络是一种递归神经网络，通过内部反馈存储和利用系统过去时刻的输入输出信息，其网络结构如图2所示。而LSTM单元能有效的解决RNN训练时的梯度爆炸和梯度消失问题，因此，选取LSTM型RNN神经网络对脑电时间序列信号进行分类。LSTM结构如图3所示，其计算过程如下：

(1)计算当前时刻的候选记忆单元值其计算公式如下：

式中，x_t为当前输入数据，h_t-1为上一时刻LSTM单元输出值，W_xc、W_hc分别表示对应输入数据和上一时刻LSTM单元输出的权值，b_c为偏置量。

(2)输入门是用于控制当前数据输入对记忆单元状态值的影响。所有门的计算除了受当前输入数据x_t和上一时刻LSTM单元输出值h_t-1影响外，还受上一时刻记忆单元值c_t-1的影响。输入门i_t值计算公式如下：

i_t＝σ(W_xix_t+W_hih_t-1+W_cic_t-1+b_i) (10)

式中，c_t-1为上一时刻候选记忆单元值，W_xi、W_hi和W_ci分别表示对应输入门数据、上一时刻LSTM单元输入门和上一时刻输入门记忆单元值的权值，b_i为偏置量。

(3)遗忘门是用于控制历史信息对当前记忆单元状态值的影响，其计算公式如下：

f_t＝σ(W_xfx_t+W_hfh_t-1+W_cfc_t-1+b_f) (11)

式中，W_xf、W_hf和W_cf分别表示对应遗忘门、上一时刻LSTM单元遗忘门和上一时刻遗忘门记忆单元值的权值，b_i为偏置量。

(4)计算当前时刻记忆单元状态值c_t，其计算公式如下：

式中，e表示逐点乘积，记忆单元状态更新取决于自身状态c_t-1和当前的候选记忆单元值并且通过输入门和遗忘门分别对这两部分因素进行调节。

(5)输出门用于控制记忆单元状态值的输出，其计算公式如下：

o_t＝σ(W_xox_t+W_hoh_t-1+W_coc_t-1+b_o) (13)

式中，W_xo、W_ho和W_co分别表示对应输出门、上一时刻LSTM单元输出门和上一时刻输出门记忆单元值的权值，b_o为偏置量。

(6)依式(14)计算LSTM单元的输出：

h_t＝o_t ⊙tanh(c_t) (14)

式中，σ一般取logistic sigmoid函数。

通过引入门控的设计，LSTM单元有保存，读取、重置和更新长距离历史信息的能力。LSTM网络采用通过时间的反向传播算法(Back Propagation Through Time,BPTT)进行训练。

本发明具有如下优点:

本发明为采用小波包和LSTM型RNN神经网络相结合的脑电信号识别方法。由于小波包具有任意多尺度、能够避免小波变换固定时频分解的缺陷(如高频部分频率分辨率低)等特点，基于小波包分解并使用改进的距离准则能够获得小波包最优子空间，且可同时获得最优子空间的小波包系数，以其构成脑电时频特征，则保留了原始脑电信号的时序特征。而LSTM型RNN神经网络很好地解决了RNN训练时的梯度爆炸和梯度消失问题，以之为分类器能真正有效利用脑电特征的时序信息，从而提高脑电信号的识别率。

附图说明

图1为BCI系统基本结构示意图；

图2为RNN网络结构图；

图3为LSTM结构图；

图4为脑电采集时序图；

图5为想象右手脑电信号C3、C4导联的平均功率随时间变化图；

图6为想象左手脑电信号C3、C4导联的平均功率随时间变化图；

图7为小波包空间结构分解图；

图8表示小波包子空间评价值排序结果图；

图9为本发明的实施流程图。

具体实施方式

本实施例是在Matlab2013a的仿真环境下进行的。

实施例采用“BCI Competition 2003”的Data set III数据集。实验采集了280次运动想象实验，其中包含想象左右手运动任务各90次的训练数据和包含想象左右手运动任务各50次的测试数据。所有的实验都在同一天完成，每两次实验中间均有几分钟的间隔，实验时序图如图3所示，每次实验过程持续9s。在起初的2s受试者保持休息状态，t＝2s显示器上开始出现1个持续1s的十字光标，同时会伴随1个提示实验开始的声音信号(受试者准备开始想象任务)。在t＝3～9s时，有1个指示左右方向的箭头代替十字光标作为指令，同时受试者按照指令所提示的方向通过想象左右手运动控制进度条向箭头指示方向移动。通过位于C3，C4前后(前“+”后“-”)各2.5cm位置的2对电极记录脑电信号。脑电信号的采样频率是128Hz，已经过0.5～30Hz带通滤波。

1.平均功率谱法分析脑电时域特性

根据公式(1)、(2)分别计算想象左手运动脑电信号C3、C4导联在0～9s内每个采样点对应的平均功率，即PLC₃和PLC₄，以及想象右手运动脑电信号C3、C4导联在0～9s内每个采样点对应的平均功率，即PRC₃和PRC₄，其结果如图2和图4所示。从图中可以看出想象左、右运动的平均功率都在3.5s～8s时间段内存在明显区别。故选取3.5s～8s段内的采样点为512的脑电信号进行特征提取。

2.对步骤1得到的脑电信号首先进行小波包分解，然后利用改进的距离准则获取小波包最优子空间，并选取最优子空间的小波包系数构成分类脑电特征。

步骤2.1，对脑电信号进行小波包分解。

db4小波属于Daubechies小波系，是具有紧支撑集的规范正交小波，能够对EEG进行有效分析，故选取db4小波基函数按照公式(3)、(4)对步骤1得到的脑电信号进行5层小波包分解得到小波包树和小波包系数，小波包空间结构如图6所示。

步骤2.2，基于改进的距离准则获取小波包最优子空间，

按照公式(5)、(6)求得平均类内距离S_w和平均类间距离S_b，然后根据准则(7)求得的评价值J_Al。图7为J^* _A1～J^* _Ad为计算小波包子空间对应的评价值排序。

步骤2.3，选取最优子空间的小波包系数构成分类脑电特征。

每个小波包的评价值J_A根据公式(7)由小波包子空间内所有评价值J_Al计算得到，然后按式(8)计算第j层第k个小波包子空间对所有样本小波包评价平均值J_A，以表征对应的特征频段对分类的贡献度。

将小波包最优子空间对应的评价平均值按从大到小重新排列，即J^* _A1～J^* _As，表征信号的可分性从大到小的排列。小波包子空间个数和权重因子α的取值对分类复杂度有一定影响，实际中应综合考虑。本施例实验数据选取Data set III数据集280组中C3、C4导脑电数据进行5折交叉验证实验。不同的小波包子空间个数和不同的权重因子下的实验结果如表1所示，并在此基础上求出平均分类准确度。

表1权重因子α及特征空间维数对分类正确率的影响

表1中分别列举了α取0.5、0.7、0.9时，分别选取较大可分性的前2～5个小波包子空间对应的小波包系数作为脑电特征，均得到了较高的正确率。当α取0.7且特征维数取64时，得到的结果最好，平均正确率达到94.4％，这说明当α取0.7时，脑电数据的类间距离和类内距离的关注程度设置较合适，故选取4个小波包子空间的小波包系数作为脑电特征最合适。

通过改进的距离准则选取4个小波包子空间分别为第5层的第1、7、8、13个，其大致对应[0.5Hz，1.4Hz]、[6.0Hz，7.9Hz]、[11.6Hz，12.5Hz]几个频段。从这几个频段可以看出选取的频段与想象左、右手运动中分类特征较明显的Mu节律(8～12Hz)频段比较接近。故由这四个频段组成脑电的特征集D_M*N，其中,M＝16，表示一个小波包子空间小波包系数的维度，N＝8，表示分别从选取C3、C4导中第1、7、8、13个小波包子空间的小波包系数。

本次实验数据为同一个人的脑电，所以小波包最优子空间的选取只在训练阶段实施，在测试与使用阶段直接利用选取好的最优子空间的小波包系数构成脑电特征，这样可以较大幅地减少程序运算量，有利于提高程序运行速度。

3.利用LSTM型RNN神经网络对脑电特征进行分类

对步骤2得到的脑电信号特征集D_M*N，其维数为16*8，故构建一个三层的LSTM型RNN网络结果，因同一时刻特征向量维数为8，故输入神经元个数为8，输出为2分类，故输出神经元个数为1。本实施例实验数据选取Data set III数据集280组中C3、C4导脑电数据，基于WPT进行特征提取，并LSTM行RNN神经网络为分类器进行5折交叉验证实验研究，验证该网络的分类性能。

关于隐含层神经元的个数确定，采用了经验+实验的方式，即先根据经验确定大致范围，然后进行多次实验进行确定。

表2隐含层神经元个数对分类正确率的影响

从表2中可知，当隐含层神经元个数为9时，取得了94.4％分类准确率。故本发明隐含层神经元个数设为9。

将上述小波包与LSTM神经网络相结合组成WPT-LSTM脑电信号识别方法，该方法能够充分利用基于小波包变换提取的脑电特征的时序信息，从而提高了脑电信号的分类准确率，为脑电信号的识别提供一种新思路。

如图9所示为本发明的实施流程图。

Claims (1)

1.基于小波包和LSTM型RNN神经网络的脑电识别方法，其特征在于：

首先，利用平均功率谱法分析脑电信号的时域特性，确定有效的时域范围；其次，先利用小波包对脑电信号进行分解，然后采用改进的距离准则获取小波包最优子空间，并选取小波包最优子空间的小波包系数构成脑电特征；最终，利用LSTM神经网络作为分类器对脑电信号进行识别分类；

其特征在于：本方法的具体步骤如下，

步骤1，信号采集及预处理；

首先按照国际标准10-20系统来安放电极，通过脑电采集设备采集n导脑电信号，并将采集到的信号通过有限脉冲响应(Finite Impulse Response,FIR)滤波器进行8-30Hz带通滤波，滤波后信号为其中N为样本点总数，n为脑电导联数目，m为采样点个数，为第i导联第j个采样点滤波后的脑电信号，t＝{1,2,…,N}，i＝1,2,…,n，j＝{1,2,…,m}；

步骤2，平均功率谱法分析脑电时域特性；

EEG信号通过电极帽采集并以电压幅值形式保存，其瞬时功率计算如下：

式中，P_i ^j(t)表示第t个样本第i导联第j个采样点脑电信号的瞬时功率；

假设P_i ^j为N次实验第i导联第j个采样点脑电信号的平均功率，则P_i ^j表示为：

根据式(2)分别计算出每个导联脑电信号的平均功率，并选取区别明显时间段的脑电信号进行特征提取；

步骤3，对步骤2中选取的脑电信号进行小波包分解，然后利用改进的距离准则获取小波包最优子空间，并选取最优子空间的小波包系数构成脑电特征；

步骤3.1，对脑电信号进行小波包分解；

信号在第j层第k点处的小波包分解系数递推公式为：

其中，n为滤波器系数序号，i∈{0,1,2…2^j-1}为小波包子空间序号，j∈{0,1,…}为分解尺度，k∈{1,2,…}为位置指标；g(n)和h(n)为一对正交镜像滤波器；经过小波包分解，原信号按频带划分为若干小波包子空间，各子空间对应的频段为：式中，f_s为信号采样率；

步骤3.2，基于改进的距离准则获取小波包子空间评价值；

设c个类别ω₁,ω₂…,ω_c的样本矢量集为{x^(p,b),p＝1,2,...,c,b＝1,2,...,N_p}，其中为ω_p类中的第b个m维样本矢量，为x^(p,b)的第l个分量；N_p为ω_p类中的特征矢量的个数；

首先，计算平均类内距离S_w如下：

式中，表示ω_p类中μ^(p)的第l个分量的均值，l＝1,2,...,m；

其次，设ω_p类样本的平均值矢量为μ^(p)，对c个类别的平均类间距离S_b定义为：

式中，μ_l为样本总体均值矢量的第l个分量；

同时引入权重因子α∈[0,1]表征对类内距离与类间距离的关注程度，得到改进的距离准则如下：

J_Al＝α*S_b-(1-α)*S_w (7)

利用该准则求得的评价值J_Al越大，即表示对应的可分性越好；

步骤3.3，选取最优子空间的小波包系数构成脑电特征；

按式(8)计算第j层第k个小波包子空间对所有样本小波包评价的平均值J_A，以表征对应的特征频段对分类的贡献度；

将小波包子空间对应的评价值按从大到小重新排列，即J^* _A1～J^* _As，即为信号可分性从大到小的排列；将评价值最高的前四个小波包子空间系数抽取出来作为脑电特征；

这里需要说明的是，由于EEG信号的个体差异性，不同实验者的最优子空间可能不同，该方法能够根据不同个体的实际信号自动选择合适的最优子空间；对于同一个个体来说，小波包最优子空间的选取只在训练阶段实施，在测试与使用阶段直接利用选取好的最优子空间的小波包系数构成脑电特征，这样大幅提高计算效率以满足脑电信号在线识别的要求；

步骤4，利用LSTM型RNN神经网络进行脑电特征分类；

RNN神经网络是一种递归神经网络，通过内部反馈存储和利用系统过去时刻的输入输出信息；而LSTM单元能有效的解决RNN训练时的梯度爆炸和梯度消失问题，因此，选取LSTM型RNN神经网络对脑电时间序列信号进行分类，其计算过程如下：

(1)计算当前时刻的候选记忆单元值其计算公式如下：

式中，x_t为当前输入数据，h_t-1为上一时刻LSTM单元输出值，W_xc、W_hc分别表示对应输入数据和上一时刻LSTM单元输出的权值，b_c为偏置量；

(2)输入门是用于控制当前数据输入对记忆单元状态值的影响；所有门的计算除了受当前输入数据x_t和上一时刻LSTM单元输出值h_t-1影响外，还受上一时刻记忆单元值c_t-1的影响；输入门i_t值计算公式如下：

i_t＝σ(W_xix_t+W_hih_t-1+W_cic_t-1+b_i) (10)

式中，c_t-1为上一时刻候选记忆单元值，W_xi、W_hi和W_ci分别表示对应输入门数据、上一时刻LSTM单元输入门和上一时刻输入门记忆单元值的权值，b_i为偏置量；

(3)遗忘门是用于控制历史信息对当前记忆单元状态值的影响，其计算公式如下：

f_t＝σ(W_xfx_t+W_hfh_t-1+W_cfc_t-1+b_f) (11)

式中，W_xf、W_hf和W_cf分别表示对应遗忘门、上一时刻LSTM单元遗忘门和上一时刻遗忘门记忆单元值的权值，b_f为偏置量；

(4)计算当前时刻记忆单元状态值c_t，其计算公式如下：

式中，⊙表示逐点乘积，记忆单元状态更新取决于自身状态c_t-1和当前的候选记忆单元值并且通过输入门和遗忘门分别对这两部分因素进行调节；

(5)输出门用于控制记忆单元状态值的输出，其计算公式如下：

o_t＝σ(W_xox_t+W_hoh_t-1+W_coc_t-1+b_o) (13)

式中，W_xo、W_ho和W_co分别表示对应输出门、上一时刻LSTM单元输出门和上一时刻输出门记忆单元值的权值，b_o为偏置量；

(6)依式(14)计算LSTM单元的输出：

h_t＝o_t⊙tanh(c_t) (14)

式中，σ取logistic sigmoid函数；

通过引入门控的设计，LSTM单元有保存，读取、重置和更新长距离历史信息的能力；LSTM网络采用通过时间的反向传播算法进行训练。