CN112541415B - 基于符号传递熵和图论的脑肌功能网络运动疲劳检测方法 - Google Patents

Abstract

本发明公开了一种基于符号传递熵和图论的脑肌功能网络运动疲劳检测方法；首先，通过64导联脑电帽和肌电采集仪采集了10通道脑电信号和尺侧腕屈肌、屈指浅肌、桡侧腕屈肌的表面肌电信号。结合符号传递熵和图论知识建立了脑肌功能网络，并提取了网络特征向量用于K最近邻算法训练分类器，最终用于检测运动疲劳。该方法克服了传统的基于生物电信号的运动疲劳检测方法中未全面考虑脑肌协作进行运动控制的缺点，将脑电信号和肌电信号结合进行运动疲劳检测，结果显示该方法在运动疲劳检测中具有较高的准确率。

Description

基于符号传递熵和图论的脑肌功能网络运动疲劳检测方法

技术领域

本发明属于模式识别领域，涉及一种人体运动疲劳的检测，特别地，涉及一种基于脑电和肌电信号对运动疲劳的检测方法。

背景技术

肌肉疲劳的定义是在连续收缩过程中保持一定强度的能力下降，或在间歇性收缩过程中未能达到初始强度水平，这反映在肌肉的电活动中。该术语表示感觉比预期的难度更大或需要更多的努力。对于健康人来说，反复或持续的肌肉激活会导致疲劳，这在运动或日常生活中很常见。最直观的效果是降低了肌肉产生力量的能力，这增加了肌肉受伤的可能性。这种表现也与运动功能障碍的患者如中风和癌症密切相关。因此，准确有效地检测肌肉疲劳至关重要。

自从首次确定运动过程中的脑电图(electroencephalogram，EEG)和肌电图(electromyogram，EMG)之间的关系以来，大脑和肌肉的交流已成为重要的研究课题，皮层肌肉连贯性(Cortical-muscle coherence，CMC)与身体运动密切相关。同样的研究兴趣也已经应用于肌肉疲劳，即在最大程度的自愿性收缩引起的疲劳后，CMC下降。随后，他们基于疲劳是一个渐进过程这一事实优化了实验范式，并得出了相同的结论。相反，最近的一项研究为志愿者设定了个体疲劳阈值，发现CMC可能增加也可能减少，这可能反映出受试者在疲劳发作后的维持能力。同时，EEG和EMG的功率降低，EEG相干性降低以及EMG循环平稳性提高也相继被报道。结论的不一致可能是因为脑肌网络的信息流方向有关，以及大脑和肌肉都参与了这一过程的事实。

为了解决上述方法的局限性，我们提出了一种结合符号传递熵(symbolictransfer entropy，STE)和图论把EEG和EMG当作一个整体建立脑肌网络的方法。由于不对称(定向)和过程(动态)概率计算，STE方法在检测非线性相互作用方面特别有效，已被应用于EEG和EMG信号。图论为脑网络组织的定量研究提供了一种有力的方法。平均度数、全局效率和全局聚类被计算作为特征，用于之后的K最近邻算法训练肌肉疲劳分类器。

本发明考虑了运动控制时皮层和肌肉双向传递以及皮层和肌肉同时参与的事实，构建了运动控制中的脑肌网络，研究中枢神经系统的调节规律，从神经系统的内源角度上为检测人体运动疲劳提供了一种新的方法。

发明内容

本发明针对为了填补运动疲劳检测方式的不足，提出了一种基于符号传递熵和图论的脑肌功能网络运动疲劳检测方法。从大脑皮层和肌肉共从参与运动控制的事实出发，建立了脑肌功能网络，计算了平均度数、全局效率和全局聚类作为特征进行检测。本方法相比传统方法有着不受电极数量影响、准确率高的特点。由于运动控制过程中，大脑皮层和肌肉共同参与，且互相有信息流入以及流出，单纯的只考虑脑电信号或者只考肌电信号都不能准确的表征神经系统的特征，同时皮层和肌肉间的信息交流具有方向性，无向的耦合计算方法遗失了信息交流中的方向信息。本发明提出了一种结合符号传递熵和图论计算网络特征的方法，根据采集的不同状态下脑电信号和肌电信号判别是否进入肌肉疲劳状态。首先，采集了被试在正常状态下和肌肉疲劳状态下的EEG和EMG信号，其次对获得的原始信号进行了预处理以提高信号质量。根据脑-肌耦合相干规律，选择γ信号频段用于建立脑肌功能网络，并提取了脑肌功能网络相关特征平均度数、全局效率和全局聚类系数作为运动疲劳检测的特征，最后根据K最近邻算法，训练了一个肌肉疲劳检测分类器。

为了实现以上内容，本发明方法主要包括以下步骤：

步骤(1)，通过肌电传感器和脑电帽以1000Hz采样频率同步实时采集正常状态下以及疲劳状态下的表面肌电信号和64通道脑电信号，脑电电极按照国际10-20系统放置。

步骤(2)，对原始信号进行去噪滤波处理：

1)对肌电信号进行小波去噪：首先对原始信号进行小波分解。

给定信号f(k)，k为采样点。其小波变换表达式为：

其中j为尺度参数，J为最佳尺度，N为时间序列长度，ψ代表小波函数。

本发明中使用Mallet算法来实现小波变换，表达式为:

其中h和g分别是尺度函数和小波函数对应的低通和高通滤波器；Sf(0,k)指原始信号；Sf(j,k)指尺度系数；Wf(j,k)指小波系数，简记为w_j,k。

对小波系数进行阈值量化，选择阈值为硬阈值，并对这一层的小波系数进行软阈值量化处理，软阈值操作表达式为：

其中T为阈值，T选取为0.2。

最后对软阈值处理后的小波分量进行重构得到去噪后的肌电信号。重构公式表示为：

其中和/>分别为h和g的共轭。

2)对脑电信号进行独立成分分析去除脑电信号中的眼动、肌电伪迹。独立成分分析的

数学表达式为：

x＝As (5)

其中x为原始数据，A称为混淆矩阵，s则为独立成分分析后的独立成分分量。去除眼动、肌电分量后对独立成分分量进行重构，最后得到预处理后的脑电信号。

3)最后肌电信号和脑电信号分别通过12-30Hz和30-45Hz的带通滤波器被分为beta频段和gamma频段。

步骤(3),应用符号传递熵计算每对信号之间的耦合强度值，并建立符号传递熵邻接矩阵。符号传递熵计算过程如下：

1)给定时间序列信号X＝{x₁,x₂,...x_i}和Y＝{y₁,y₂,...y_i}，i代表采样点。首先进行符号化运算，符号化运算被定义为：

其中S(i)代表符号化后的信号，min和max分别代表时间序列信号的最小值和最大值，delta定义为p被设置为50。

符号化运算后的信号分别为：

2)对符号化后的信号进行传递熵运算，表达式如下：

其中P()代表概率分布，Y指向X的符号传递熵可由相同方式得到。在这种情况下，STE的值表示X和Y之间的耦合关系的强度。请注意，当X和Y相同时，STE＝0。

3)为每对脑电信号、每对肌电信号以及每对脑电-肌电信号计算符号传递熵值，构建加权脑肌功能网络。

步骤(4),根据成本效率Ce确定网络连接传递熵阈值，将步骤(3)中得到的加权脑肌功能网络二值化，成本效率定义为：

Ce＝E_g-D (8)

其中D代表网络密度，是网络实际边数与最大边数的闭值。Eg代表全局效率，计算公式如下：

其中L_i,j代表两个节点直接的最短路径，N是网络中的节点数。

在保证成本效率大的条件下选择网络连接传递熵阈值Th＝0.24，大于阈值置为1，小于阈值置为0，构建二值脑肌功能网络。

步骤(5),基于图论知识对二值脑肌功能网络进行特征提取，分别提取了平均度数、全局效率、全局聚类系数。

度数分为出度和入度，入度代表由其它节点指向此节点的边的总数，出度代表由此节点指向其它节点的边的总数。平均度数De的计算公式如下：

其中De_i表示i节点的度。

全局效率是信息流的标量度量，被定义为给定网络中所有最短路径长度的逆，计算公式如公式(9)。

聚类系数表示节点的邻居成为彼此的邻居的可能性，若聚类系数低，则表明网络连接不紧密，计算公式如下：

其中B_i表示节点i的邻居节点之间的连接数，K_i表示节点i的邻居节点的数目。

步骤(6),针对运动疲劳检测问题，我们基于KNN算法提出了扩展的K最近邻算法来设计分类器，分类器的设计如下：

1)首先从每个类别中找出m个代表向量，计算平均向量构成类中心向量组；类中心向量组的四元表示形式为：{C(d_i),R(d_i),N(d_i),X(d_i)}，其中d_i为样本类别的中心向量，C(d_i)表示类别；R(d_i)为类别覆盖区域半径；N(d_i)为类别覆盖样本点的数目；X(d_i)则为该类别所有样本数据串联后的向量；

2)为了使分类器具有逐步学习的能力，在分类过程中，对于类别得到验证的样本点，距离最近的样本加到已有的模型类中，距离度量为上确界距离，表示为：

d(x,y)＝max(|x₁-x₂|,|y₁-y₂|,|z₁-z₂|) (13)

不能识别的样本点，则在模型类中重新计算四元组。

3)计算待分类样本与各类中心点的距离；在确定待分类样本的类别时，让权重系数乘每一个中心点；由加权后的数据确定待识别样本的类别；

最后将步骤(5)中计算得到的3维特征向量输入扩展的K最近邻算法进行训练分类器对运动疲劳分类。

本发明与已有的运动疲劳检测相比，具有如下特点：

将EEG和EMG信号用于运动疲劳检测，这种方法直接提取人体的表面生物电信号来传递人体的机能信息，提高了疲劳检测的实时性和主动性，从而提高了操作的准确性和便捷性，降低了自我体能评估失误造成的风险，使人体可以更准确低指定运动计划并保护肌肉不受危害。使用K最近邻算法设计分类器可以降低时间复杂度，可满足实时性要求较高的运动疲劳检测要求。

附图说明

图1为实验进程示意图；

图2为发明工作流程示意图；

图3为建立的加权邻接矩阵示意图；

图4为建立的二值脑肌功能网络示意图；

具体实施方式

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

实施范式和实验环境说明：

E1非疲劳状态范式：被试者自然地坐在椅子上，小臂自然水平放在受试者正前方桌子上，同时手中握住握力器。根据电脑音频信号提示开始，在1s内完成30％最大握力，随后保持5s，根据电脑音频信号提示结束。每次实验之间间隔2分钟休息时间以避免发生肌肉疲劳。

E2疲劳状态范式：被试者自然地坐在椅子上，小臂自然水平放在受试者正前方桌子上，同时手中握住握力器。根据电脑音频信号提示开始，在1s内完成30％最大握力，随后保持5s，根据电脑音频信号提示结束。每次实验之间间隔2s，同时根据Borg疲劳自量表记录疲劳指数，疲劳指数大于14的实验数据被记为疲劳数据。

被试组成：青年组10人，年龄22±3，男性受试者8人，女性受试者2人。

为降低环境因素对实验设备精度的影响，本实验采集数据的时段均选择在下午，实验采集环境温度为25±4℃。附图1为实验进程示意图。

附图2为发明实施工作流程示意图，实施过程包括如下步骤：

步骤一：分别采集在2种实验范式下被试的EEG和EMG信号。每个受试者按照E1、E2范式进行实验，每个范式下每个受试者分别采集十组时长为6s的数据。E1范式中，每次实验结束后休息2分钟再进行下一次实验避免发生疲劳。E2范式中，每次实验结束后间隔2s记录疲劳指数后继续实验，目的维持疲劳状态。测试过程中被试被要求避免除握力需要的其它动作，如：咬牙、面部表情、头动等，以防止伪迹的引入，若发生此类动作或者实验人员交谈等干扰事件，则终止该次实验，在休息后重测。

步骤二：对步骤一中获得的原始EEG和EMG信号进行预处理以得到更高质量的生物电信号，并分别使用12-30Hz和30-45Hz带通滤波器提取Btea和Gamma频段数据；

步骤三：对预处理后的EEG和EMG信号进行符号传递熵运算得到13×13的加权邻接矩阵构建加权脑肌功能网络，附图3为加权邻接矩阵示意图，其中10个脑电通道信号(F3、F4、FC3、FC4、C3、C4、CP3、CP4、P3、P4)和3个肌电通道信号(尺侧腕屈肌、屈指浅肌、桡侧腕屈肌)。选取阈值水平为0.24，将低于阈值的边界设为0，高于阈值的边界设为1，最终得到二值脑肌功能网络，附图4为二值脑肌功能网络示意图。

步骤四：按照国际10-20系统64导联图的划分方式，本文分别提取了二值脑肌功能网络特征向量。网络特征提取是所有步骤中最为关键的一步，首先由步骤三中得到的二值脑肌功能网络拓扑图计算计算平均度数、全局效率、全局聚类系数。由该参数来对被试是否进入肌肉疲劳状态进行检测，检测结果表明该指标确实具有良好的区分作用。

表1检测准确率

由表1中结果可知，本文提出的脑肌功能网络特征在运动疲劳检测中具有优秀的检测准确率，尤其在Gamma频段准确率达到了93.3％，这也符合现存的研究发现的运动时EEG和EMG信号的Beta和Gamma频段作为特征频段，而疲劳后Gamma频段更为活跃的事实。总体看来，脑肌功能网络特征能较好地刻画不同肌肉状态的人群神经系统活动规律，对人体运动疲劳的检测精度不亚于传统方法，因此本方法可适用于不同人群和个体的运动疲劳检测。

Claims (1)

1.基于符号传递熵和图论的脑肌功能网络运动疲劳检测方法，其特征在于：把大脑皮质和肌肉群当作一个整体系统，结合脑电信号和肌电信号建立了脑肌功能网络，并提取网络特征对肌肉是否疲劳进行识别；该方法包括如下步骤：

步骤(1)，通过肌电传感器和脑电帽以1000Hz采样频率同步实时采集正常状态下以及疲劳状态下的表面肌电信号和64通道脑电信号，脑电电极按照国际10-20系统放置；

步骤(2)，对原始信号进行去噪滤波处理：

1)对肌电信号进行小波去噪：首先对原始信号进行小波分解；

给定信号f(k)，k为采样点；其小波变换表达式为：

其中j为尺度参数，J为最佳尺度，N为时间序列长度，ψ代表小波函数；

使用Mallet算法来实现小波变换，表达式为:

其中h和g分别是尺度函数和小波函数对应的低通和高通滤波器；Sf(0,k)指原始信号；Sf(j,k)指尺度系数；Wf(j,k)指小波系数，简记为w_j,k；

对小波系数进行阈值量化，选择阈值为硬阈值，并对这一层的小波系数进行软阈值量化处理，软阈值操作表达式为：

其中T为阈值，T选取为0.2；

最后对软阈值处理后的小波分量进行重构得到去噪后的肌电信号；重构公式表示为：

其中和/>分别为h和g的共轭；

2)对脑电信号进行独立成分分析去除脑电信号中的眼动、肌电伪迹；独立成分分析的数学表达式为：

x＝As (5)

其中x为原始数据，A称为混淆矩阵，s则为独立成分分析后的独立成分分量；去除眼动、肌电分量后对独立成分分量进行重构，最后得到预处理后的脑电信号；

3)最后肌电信号和脑电信号分别通过12-30Hz和30-45Hz的带通滤波器被分为beta频段和gamma频段；

步骤(3),应用符号传递熵计算每对信号之间的耦合强度值，并建立符号传递熵邻接矩阵；符号传递熵计算过程如下：

1)给定时间序列信号X＝{x₁,x₂,...x_i}和Y＝{y₁,y₂,...y_i}，i代表采样点；首先进行符号化运算，符号化运算被定义为：

其中S(i)代表符号化后的信号，min和max分别代表时间序列信号的最小值和最大值，delta定义为p被设置为50；

符号化运算后的信号分别为：

2)对符号化后的信号进行传递熵运算，表达式如下：

其中p()代表概率分布，Y指向X的符号传递熵可由相同方式得到；STE的值表示X和Y之间的耦合关系的强度；当X和Y相同时，STE＝0；

3)为每对脑电信号、每对肌电信号以及每对脑电-肌电信号计算符号传递熵值，构建加权脑肌功能网络；

步骤(4),根据成本效率Ce确定网络连接传递熵阈值，将步骤(3)中得到的加权脑肌功能网络二值化，成本效率定义为：

Ce＝E_g-D (8)

其中D代表网络密度，是网络实际边数与最大边数的闭值；E_g代表全局效率，计算公式如下：

其中L_i,j代表两个节点直接的最短路径，N是网络中的节点数；

在保证成本效率大的条件下选择网络连接传递熵阈值Th＝0.24，大于阈值置为1，小于阈值置为0，构建二值脑肌功能网络；

步骤(5),基于图论知识对二值脑肌功能网络进行特征提取，分别提取了平均度数、全局效率、全局聚类系数；

度数分为出度和入度，入度代表由其它节点指向此节点的边的总数，出度代表由此节点指向其它节点的边的总数；平均度数De的计算公式如下：

其中De_i表示i节点的度；

全局效率是信息流的标量度量，被定义为给定网络中所有最短路径长度的逆，计算公式如公式(9)；

聚类系数表示节点的邻居成为彼此的邻居的可能性，若聚类系数低，则表明网络连接不紧密，计算公式如下：

其中B_i表示节点i的邻居节点之间的连接数，K_i表示节点i的邻居节点的数目；

步骤(6),针对运动疲劳检测问题，基于KNN算法提出了扩展的K最近邻算法来设计分类器，分类器的设计如下：

1)首先从每个类别中找出m个代表向量，计算平均向量构成类中心向量组；类中心向量组的四元表示形式为：{C(d_i),R(d_i),N(d_i),X(d_i)}，其中d_i为样本类别的中心向量，C(d_i)表示类别；R(d_i)为类别覆盖区域半径；N(d_i)为类别覆盖样本点的数目；X(d_i)则为该类别所有样本数据串联后的向量；

2)为了使分类器具有逐步学习的能力，在分类过程中，对于类别得到验证的样本点，距离最近的样本加到已有的模型类中，距离度量为上确界距离，表示为：

d(x,y)＝max(|x₁-x₂|,|y₁-y₂|,|z₁-z₂|) (13)

不能识别的样本点，则在模型类中重新计算四元组；

3)计算待分类样本与各类中心点的距离；在确定待分类样本的类别时，让权重系数乘每一个中心点；由加权后的数据确定待识别样本的类别；

最后将步骤(5)中计算得到的3维特征向量输入扩展的K最近邻算法进行训练分类器对运动疲劳分类。