CN111814544B - 基于局部平衡极限学习机的动态运动想象脑电识别方法 - Google Patents

Abstract

本发明公开了一种基于局部平衡极限学习机的动态运动想象脑电识别方法，本发明为了提高手部动作分类的性能，采用了一种局部平衡ELM和基于SL量化FBN的新方法。与传统思想不同，本发明从上述加权的FBN提取了两种二值子网络的拓扑特征。此外，本发明指出了ELM融合中的两个不足，并提出了相应的局部平衡ELM。为了达到两个互补特征在ELM特征空间的最佳融合，以及二值化的最佳阈值和正则化ELM的最佳正则化参数，本专利改进了LOO方法的计算复杂度，使用LOO来收敛上述参数的经验范围。

Description

基于局部平衡极限学习机的动态运动想象脑电识别方法

技术领域

本发明属于脑电信号特征分析领域，公开了一种基于局部平衡极限学习机的动态运动想象脑电识别方法。

背景技术

脑机接口(Brain-Computer Interface，BCI)技术在通讯和控制、神经物理调节、康复、游戏和娱乐等领域中得到了十分广泛的应用。

在通信控制，康复训练领域中，运动想象(Motor Imagery，MI)的应用最为广泛。研究指出，当前基于MI的BCI中使用的MI信号几乎都是静态MI(Static MI，sMI)。然而，在医学神经康复领域，MI可分为两个子领域：sMI和动态运动想象(Dynamic Motor Imagery，dMI)。sMI是传统的MI，MI主体仅简单地对动作进行心理表征，而无需实际执行，而dMI则在心理表征的同时增加了实际的动作，在康复训练等领域更具应用价值，在此方面的研究将具有极大的意义。

BCI系统应用的效果取决于系统的解码能力。BCI系统通常分为四个部分：数据采集模块，数据预处理和增强模块，特征提取模块和分类模块。

特征提取是最重要的过程，其效果直接影响分类的准确性。相对于EEG信号分析的传统方法，时域分析，例如能量熵或功率、频域分析、Hilbert变换，以及时频分析，例如小波变换(wavelet transform，WT)。时间序列之间的相互依赖性分析成为更前沿的研究领域，用信息论、相位同步、广义同步和格兰杰因果关系等方法量化不同大脑位置之间的功能连通性(functional connectivity，FC)构建脑功能网络(functional brain network，FBN)，然后辅以图论以提取脑功能网络的拓扑特征。似然同步性(synchronization likelihood，SL)是目前最优异的量化方法之一，它可以量化动态非线性系统中相似模式出现的可能性。

二进制FBN通常是通过消除低于某个阈值的弱FC来构造的。近来，用最小生成树(minimum spanning tree，MST)构造二进制FBN的方法能有效地显示驾驶员的警觉状态和困倦状态之间的差异。MST可以连接原始加权网络中的所有节点，同时最小化路径之和且不形成回路。通过这种方法构建的二值FBN具有相同的节点数和FC数，可以直接比较组之间的网络属性，从而避免了上述方法上的偏差。

尽管通过任意阈值构建二值FBN可能会导致偏差，但是阈值也可以理解为调整网络复杂度的参数。两种方法都按照各自的定义消除了弱FC，这意味着它们在消除弱连接上持有不同的“方向”，从而使得两者生成的二进制FBN之间的相关性很低。故而本发明将这两种互补网络的拓扑特征进行融合。

对于分类模块，近期研究通常致力于深度学习网络上，例如卷积神经网络，标准递归神经网络和长短期记忆。dMI脑电的数据样本较少，且深度学习的训练时间太长。因此，本发明选择了适用于少量样本并且训练快速的改进型极限学习机(extreme learningmachine，ELM)——局部平衡极限学习机(locals_balanced ELM，LBELM)。

发明内容

本发明针对现有技术的不足，提出了一种基于局部平衡极限学习机的动态运动想象脑电识别方法。

与MI相比，dMI在训练和神经康复领域优势明显。本发明提出了一种基于FBN、SL和ELM的新方法。与传统方法不同，本发明从加权FBN构建两种二值FBN的拓扑特征(通过阈值法消除弱连接建立二值FBN，通过构造MST构建二值FBN)。此外，本发明针对互补特征ELM的两个不足，提出了相应的解决方案。为了收敛两个互补特征的最佳融合比例、阈值法的最佳阈值和正则化ELM的最佳正则化参数，本发明采用了LOO方法并且降低了其计算复杂度。

为实现以上目的，本发明方法主要包括以下步骤：

步骤(1)采集脑电数据。实验样本频率设置为5000Hz。根据实验要求，受试者按照随机动作顺序进行曲腕、曲肘和握拳动作。

步骤(2)信号预处理。脑电信号质量的保证是准确有效研究的基础。在脑电信号采集过程中，不可避免地会受到各种信号和噪声的干扰。由于EEG信号非常微弱，因此很容易受到许多干扰源的干扰，包括眼球运动和眨眼伪影，心电图和工频干扰。因此，为了保证基于脑电信号的研究结果的可靠性，有必要对记录的信号进行相应的有效降噪步骤，以尽可能减少噪声对研究结果的干扰。

步骤(3)构建加权FBN。FBN分析是以信息论、相位同步、广义同步和Granger因果关系为基础，辅以图论等方法以表征脑区之间功能运作模式的一种重要方法。SL是最有前途的方法之一，它可以量化动态和非线性系统(如神经网络)中相似模式出现的可能性。

步骤(4)构建二值FBN。二值FBN通常是通过在一定程度上偏置的阈值以下消除弱FC来构造的。近年来，人们提出了一种利用MST构造二值FBN的新方法，以探索驾驶员的警觉状态和困倦状态之间的差异。MST是原网络稳定核的简化表示，它以最小的连接代价连接原加权网络中的所有节点，且不形成圆或环。由于MST的节点数和链路数相同，可以直接比较组之间的网络特性，从而避免了阈值方法上的偏差。尽管任意阈值可能导致偏差，但这也可以作为调整网络复杂性的参数。这两种方法都在各自的定义中去除了弱FC，这意味着它们在消除虚假连接时具有不同的“方向”，并且生成的两个二值FBN之间的相关性很低。因此本发明从加权FBN构建两种二值FBN(通过设置适当的阈值消除弱连接建立二值FBN和通过构造MST构建二值FBN)。

步骤(5)提取拓扑特征。在建立二值FBN之后，我们需要对网络特性进行深入研究。为了更深入地描述网络属性，我们需要定义一些适当的网络拓扑指标。比较和分析每个复杂网络拓扑特征，反映出FBN的功能模式。

步骤(6)使用局部平衡极限学习机分类。ELM由单隐层神经网络(single hidden-layer feedforward neural networks，SLFNs)发展而来。SLFNs的所有参数需要使用梯度下降方法进行收敛，所以SLFNs的学习速度缓慢且易于收敛到局部极小值。为克服上述缺点，黄广斌在SLFNs的进一步研究中提出了两种理论^[34]。在他们的理论中，发现输出层的偏差是不必要的，而输入权重和隐藏层偏差可以随机生成，只需要单独确定SLFNs的输出权重。

对于一个n维的，拥有N个样本，m个类别的训练样本L个隐藏单元，激活函数为g_i(s_j)的SLFNs输出为

其中s_j＝[s_j1,s_j2,…,s_jn]^T为第j个输入，t_j＝[t_j1,t_j2,…,t_jm]^T为对应于第j个输入的期望输出；b_i为第i个隐藏单元的激活阈值；c_i＝[c_i1,c_i2,…,c_in]连接了输入层和第i个隐藏单元；β_i＝[β_i1,β_i2,…,β_im]^T连接了第i个隐藏单元和输出层；o_j为SLFNs对应于第j个输入的实际输出。

SLFNs可以零误差的拟合训练样本等价于由于ELM可以随机生成隐藏层参数{c_i,b_i}，所以隐藏层输出矩阵β

其中，

和

H为隐藏层输出矩阵，为其伪逆。

由于伪逆的不稳定性，本发明应用了正则化以防止过拟合，

这里的λ是正则化系数，正则化之后ELM的隐藏层输出矩阵为

其中I为单位矩阵。

特征在输入ELM之前需要归一化。由于篇幅的限制，本文只讨论了方差归一化、两个互补特征融合和高斯分布生成ELM随机参数的情况。

假设样本s由两个互补特征x₁＝[x_1,1,x_1,2,…,x_1,a]^T和x₂＝[x_2,1,x_2,2,…,x_2,b]^T组成，即s＝[x₁；x₂]，他们的模为‖x₁‖₂和‖x₂‖₂。同时假设某个隐藏层单元的输入权重矩阵为c＝[w₁|w₂]，其中w₁＝[w_1,1,w_1,2,…,w_1,a]^T、w₂＝[w_2,1,w_2,2,…,w_2,b]^T。隐藏层参数由高斯分布N(0,1)生成，所以同理(事实上，模长会随机分布在上述值的附近)。

ELM将上述的两个互补特征整合为一个归一化的(a+b)维特征。然后计算矩阵H，观察公式(3)，发现其中只有阴影项是与输入特征有关

是隐藏层对输入特征的“理解”，其中‖x₁‖₂是特征的模长、是输入特征方向与参考方向夹角的cosine值。ELM的两个缺陷可以被发现：(1)两个“理解”之间的比值是一个定值(事实上，该定值会随机分布在上述值得附近)，导致了两个“理解”的融合权重不能调整；(2)‖x₁‖₂和‖x₂‖₂这两个“理解”是非有限区间的，除非‖x₁‖₂/‖x₂‖₂是一个定值，否则两个“理解”的融合权重将难以调整。

为了克服上述的两个缺陷，本发明在归一化时增加了一个变量p，样本归一化后变为归一化后s'仍然满足方差为一，与此同时我们可以发现矩阵H中的输入特征相关项

中两个“理解”之间的比值变为这就意味着通过改变变量p就可以实现融合比例的调节。

步骤(7)使用留一法(leave-one-out cross validation，LOO)进行参数寻优。LOO方法是机器学习中最有效的模型选择和参数优化方法之一。对于一个N个样本的数据集，LOO的基本思想是将数据集划分为N个不同的子集，轮流将其中一个子集作为测试集，其余的作为训练集N次。为了降低计算复杂度，本发明采用了预测残差平方和(predictedresidual sum of squares，PRESS)统计量，该统计量可有效计算线性模型(而不是多次再训练)的LOO方法的均方误差(mean square error，MSE)。MSE^PRESS的定义是

其中t_j和o_j为第j个样本的期望输出和实际输出。HAT_jj是HAT矩阵第j个对角线元素的值。

HAT矩阵定义为

这里的H矩阵为上文中的公式(3)。

阈值σ、比值p和正则化系数λ都是不能通过经验预先得到，为了收敛参数的经验区间，本发明选择了LOO方法。为了再进一步降低正则化ELM的HAT矩阵的计算复杂度(记为HAT_r)，奇异值分解被用于分离与带收敛参数无关的计算过程。假设H＝UDV^T是矩阵H的奇异值分解，其中D为对角矩阵，U和V为酉矩阵。

If L＞N，

HAT_r＝HH^T(HH^T+λI)^-1

＝UDV^TVDU^T(UDV^TVDU^T+λI)^-1

＝UD²(D²+λI)^-1U^T, (10)

O＝UD²(D²+λI)^-1U^TT, (11)

IfL≤N，

HAT_r＝H(H^TH+λI)^-1H^T

＝UDV^T(VDU^TUDV^T+λI)^-1VDU^T

＝UD(D²+λI)^-1DU^T

＝UD²(D²+λI)^-1U^T, (13)

O＝UD²(D²+λI)^-1U^TT, (14)

对经验集中的参数进行遍历，选择均方误差最小的参数对作为最优参数组合。值得注意的是：(1)计算公式(8)只与实际输出和HAT_r有关；(2)在L≤N和L＞N的情况下，实际输出和HAT_r是一致的(对角矩阵的乘法满足交换律，矩阵的乘法满足结合律)，因此将计算分为两种情况是不必要的；(3)实际输出和HAT_r只出现了U和D²，因此奇异值分解时只需要分解HH^T；(4)矩阵HAT_r的对角线可以直接由矩阵UD²(D²+λI)^-1⊙U的行和来计算。式中，“⊙”表示两个矩阵点乘；(5)因为D²和(D²+λI)^-1是对角矩阵，UD²和UD²(D²+λI)^-1的两次矩阵相乘可以转换为U⊙repmat(diag(D²),size(U,1),1)和UD²⊙repmat(diag((D²+λI)^-1),size(UD²,1),1)的两次矩阵点乘；我们将预计算与λ无关的项以防止重复计算：diag(D²),UD²＝U⊙repmat(diag(D²),size(U,1),1)and U^TT。

本发明与已有的诸多脑电信号的方法相比，具有如下特点：

对二值FBN的构造采用两种算法(通过消除阈值错误连接的二进制FBN和通过构造最小生成树的二进制FBN)；证明了ELM在多特征融合中的缺陷，并提出了相应的解决方案；进一步推导了基于MSE^PRESS的LOO优化算法，以求获得更高效的计算。

附图说明

图1为本发明的实施流程图；

图2为本发明所进行实验中使用的电极分布图；

图3为志愿者1做左手屈腕动作的EEG信号；

图4为志愿者1做左手屈腕动作的去噪脑电信号；

图5为志愿者1动作分类的MSE^PRESS曲面；

图6为8位志愿者的动作识别率；

图7为LOO算法的伪代码；

图8为改进版LOO算法的伪代码。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本发明的实施例采集了八位志愿者分别用左右手执行腕曲，肘曲和握拳共计6种动作的EEG信号。在数据预处理阶段使用了WT，采用基于Manhattan距离的SL构造加权FBN。二进制FBN由阈值和MST两种方法构建。使用局部平衡极限学习机融合两种二值FBN的拓扑特征。为了获得最佳参数，本发明还采用留一法(leave-one-out cross validation，LOO)方法。

本文将介绍参与者、实验任务、数据集的信息与这项工作所涉及的理论，介绍方法的配置，并说明所提出方案的性能。最后与已有的相关研究进行比较，提出预期的工作。

如图1所示，本实施例包括如下步骤：

步骤(1)采集脑电数据。在进行实验之前，应告知受试者完整的实验过程和详细信息，以便他们熟悉整个过程。告诉受试者避免实验计划中未包括的所有肢体运动，尤其是身体和头部的晃动，并尽可能减少对数据的干扰，以使收集的数据更加准确。在数据采集之前，将受试者放在已安装的电极盖上，并将导电胶注入电极中以降低阻抗。调整电极的接触条件，直到阻抗下降到阈值以下，以确保收集信号的质量。实验样本频率设置为5000Hz。实验中使用的电极分布如图2所示。根据实验要求，受试者按照步骤进行曲腕、曲肘和握拳动作。图3为志愿者1做左手WF动作的EEG信号。

步骤(2)信号预处理。脑电信号质量的保证是进行准确有效研究的基础。在脑电信号采集过程中，不可避免地会受到各种信号和噪声的干扰。由于EEG信号非常微弱，因此很容易受到许多干扰源的干扰，包括眼电，心电和工频干扰。因此，为了保证基于脑电信号的研究结果的可靠性，有必要对记录的信号进行WT去噪，以尽可能减少噪声对研究结果的干扰。如图4为志愿者1左手WF的去噪脑电信号。

步骤(3)通过SL构建加权FBN。在建立大脑功能网络模型的过程中，存在两个最重要的问题：如何确定功能网络模型中的节点以及如何量化两个节点之间的连接强度。本文研究的网络模型的节点位置是EEG电极的位置。FBN分析是以信息论、相位同步、广义同步和Granger因果关系为基础，辅以图论等方法以表征脑区之间功能运作模式的一种重要方法。SL是最有前途的方法之一，它可以量化动态和非线性系统(如神经网络)中相似模式出现的可能性。本发明选取了采集信号中的40个通道，使用SL分析方法量化两两通道间的广义同步关系。

步骤(4)构建二值FBN。为了消除节点之间的弱链接，本发明从加权FBN构建两种二值FBN(通过设置适当的阈值消除弱连接建立二值FBN，通过构造MST构建二值FBN)。通过前期研究最佳阈值在0～0.2之间，从中选取阈值构建二值FBN；同时我们将1/SL作为节点之间的路径长度构建基于Prim’s MST算法的二值FBN。

步骤(5)提取拓扑特征。

在建立二值FBN之后，我们需要对网络特性进行深入研究。为了更深入地描述网络属性，我们需要定义一些适当的网络拓扑指标。比较和分析每个复杂网络拓扑特征，反映出FBN的功能模式。

本发明选取了度、聚类系数、介数中心度、核数作为候选拓扑特征。经过分析比较后采用了最富有信息的拓扑特征—度。它在通过阈值法二值FBN和通过构造MST法二值FBN中都具有最优的分类性能。表1分别展示了阈值法二值FBN和通过构造MST法二值FBN四种候选拓扑特征的分类性能(统计量采用了ELM分类的MSE^PRESS，数值越小表明该特征在分类中越优异)；

表1

粗体表示最小的MSE^PRESS

步骤(6)使用局部平衡ELM分类。为了克服ELM融合中的两个缺陷，本发明在归一化时增加了一个变量p，样本归一化后变为归一化后s'仍然满足方差为一，与此同时我们可以发现矩阵H中的输入特征相关项——公式(7)中两个“理解”之间的比值变为这就意味着通过改变变量p就可以实现融合比例的调节。我们牺牲了输入特征的一个维度(模是一维“理解”，方向是(a-1)维/(b-1)维“理解”)来获得比例调节的能力。随着特征维数的增加，“理解”损失相对减少。该方法更适合于高维特征的融合。

步骤(7)使用LOO进行参数寻优。基于MSE^PRESS的LOO方法需要遍历所有候选参数对，最后选择MSE^PRESS值最小的参数对作为最优参数对。阈值候选区间为[0：0.01：0.2]；融合比例候选区间为10^{[-1：0.1：1]}；正则化系数候选区间为e^[-4：0.1:6]。

图7展示了原LOO优化算法的伪代码，图8展示了应用在LBELM中的改进版LOO算法的伪代码。改进版LOO算法用步骤8-13实现了原LOO算法的功能。改进主要体现在了，1)其中改进后的留一法为：令留一法中L≤N情况下的实际输出和HAT_r，与L＞N的情况下的实际输出和HAT_r相同；其中L表示局部平衡极限学习机隐藏层单元的个数，N表示输入特征数量；将原留一法中的预计算项中的矩阵相乘项E＝PU，L＞N和E＝HV，L≤N改进为矩阵点乘项E＝U⊙repmat(diag（D²),size(U,1),1)；其中P＝HH^T,H为局部平衡极限学习机的隐藏层输出矩阵，“⊙”表示两个矩阵点乘，假设H＝UDV^T是矩阵H的奇异值分解，其中D为对角矩阵，U和V为酉矩阵。

图5展示了在最优正则化系数的情况下所有参数对[σ,p]的MSE^PRESS。图5中X轴为融合比例；Y轴为阈值法采用的阈值；Z轴为在该参数情况下ELM分类产生的MSE^PRESS。由图5可知，志愿者一号的最优阈值为0.18，两个特征向量的最优融合比例为3.162。在最优参数的情况下，ELM分类的MSE^PRESS为0.1941，显著减少了分类的误差(前文的单一特征的分类MSE^PRESS分别为1.3633和1.9063)

最后8位志愿者的分类准确性如图6所示。由图可知，简单组合两个特征向量为一个向量的识别率一定在较低单一特征的识别率之上，但是不一定会高于较优的单一特征的识别率，并且简单组合的提升非常有限，验证了ELM在融合中存在的缺陷。而通过增加最优的融合变量p，使得ELM可以自由融合两个特征向量，合理的平衡两个特征向量在ELM特征空间的权重。

Claims (4)

1.基于局部平衡极限学习机的动态运动想象脑电识别方法，其特征在于，包括如下步骤：

步骤(1)采集脑电数据；

步骤(2)信号预处理；

步骤(3)通过似然同步性构建加权脑功能网络；

步骤(4)构建二值脑功能网络；

通过步骤(3)的加权脑功能网络构建两种二值脑功能网络；通过设置的最优阈值消除弱连接建立二值脑功能网络，通过构造最小生成树算法构建二值脑功能网络；

步骤(5)提取拓扑特征；

提取步骤(4)中两种二值脑功能网络各自的拓扑特征；

步骤(6)使用局部平衡极限学习机分类；

通过局部平衡极限学习机实现两种拓扑特征在极限学习机特征空间的最优比例融合，以实现分类效果的提升，并通过最优的正则化系数获得最佳的分类泛化性能；

其中局部平衡极限学习机为在极限学习机的基础上，在归一化时增加了一个变量p，样本归一化后变为归一化后s’仍然满足方差为1，得到矩阵H中的输入特征相关项c·s’；

则极限学习机中的两个“理解”之间的比值变为即通过改变变量p就实现融合比例的调节；其中两个“理解”分别为样本s由两个互补特征x₁＝[x_1，1，x_1，2，…，x_1，a]^T和x₂＝[x_2，1，x_2，2，…，x_2,b]^T组成，即s＝[x₁；x₂]，x₁、x₂的模为||x₁||₂和||x₂||₂，其中x_1，a表示第一个局部特征中第a个元素，x_2，b第二个局部特征中第b个元素；隐藏层单元的输入权重矩阵为c＝[w₁|w₂]，其中w₁＝[w_1,1，w_1,2，…，w_1,a]^T、w₂＝[w_2,1，w_2,2，…，w_2,b]^T；隐藏层参数由高斯分布N(0，1)生成，所以 同理

步骤(7)使用改进后的留一法进行参数寻优；

设置二值化阈值σ候选区间为[0：0.01：0.2]；融合比值p候选区间为10^{[-1：0.1：1]}；正则化系数λ候选区间为e^{[-4：0.1：6]}；通过改进后的留一法得到分类效果最佳的二值化阈值σ、融合比值p和正则化系数λ；

其中改进后的留一法为：令留一法中L≤N情况下的实际输出和HAT_r，与L>N的情况下的实际输出和HAT_r相同；其中L表示局部平衡极限学习机隐藏层单元的个数，N表示输入特征数量；将原留一法中的预计算项中的矩阵相乘项E＝PU，L>N和E＝HV，L≤N改进为矩阵点乘项E＝U⊙repmat(diag(D²)，size(U，1)，1)；其中P＝HH^T，H为局部平衡极限学习机的隐藏层输出矩阵，“⊙”表示两个矩阵点乘，假设H＝UDV^T是矩阵H的奇异值分解，其中D为对角矩阵，U和V为酉矩阵。

2.根据权利要求1所述的基于局部平衡极限学习机的动态运动想象脑电识别方法，其特征在于：所述的采集脑电数据；具体为：

在数据采集之前，将电极帽安装在受试者头上，并将导电胶注入电极中以降低阻抗；调整电极的接触条件，直到阻抗下降到阈值以下，以确保收集信号的质量；实验样本频率设置为5000Hz。

3.根据权利要求1所述的基于局部平衡极限学习机的动态运动想象脑电识别方法，其特征在于：信号预处理；具体为：

对记录的信号进行小波变换去噪，以尽可能减少噪声对研究结果的干扰。

4.根据权利要求1所述的基于局部平衡极限学习机的动态运动想象脑电识别方法，其特征在于：通过似然同步性构建加权脑功能网络；具体为：

选取了采集信号中的40个通道，使用似然同步性分析方法量化两两通道间的广义同步关系，构建了基于似然同步性的加权脑功能网络。