CN103584855B - 脑肌电同步采集及信息传递特性分析方法 - Google Patents

Abstract

本发明涉及康复工程领域、神经工程及运动障碍机制研究领域，具体涉及一种脑肌电同步采集及信息传递特性分析方法，所述方法包括脑肌电信号同步采集部分和信号处理部分，脑肌电信号同步采集部分包括脑电信号采集和肌电信号采集；信号处理部分包括信号预处理和脑肌电信息传递特性分析方法。本发明具有可应用性、可采纳性，在康复医学领域具有重要的应用价值。

Description

脑肌电同步采集及信息传递特性分析方法

技术领域

本发明涉及康复工程领域、神经工程及运动障碍机制研究领域，具体涉及一种在脑肌电同步采集下基于状态转移率的脑肌电信息传递特征分析方法。

背景技术

头皮脑电（EEG)是脑皮层神经电活动的表现，表面肌电（sEMG)是肌肉电活动的表现，二者分别作为人体指令中心和执行末端的代表性信号，在机体运动功能领域有着非常重要的研究意义。自从1995年Conway发现了运动过程中肌肉电信号和脑磁信号之间具有相关性，基于脑皮层和肌肉的相关同步特征探索大脑感觉运功皮层与肌肉功能联系的研究相继开展。2004年日本科学家Tomohiro Shibata等同步测量正常人和镜像运动病人在手臂动作时的sEMG和EEG信号，通过分析脑肌电信号的相干性，验证了镜像运动病人控制神经分布与正常人的差异特征；2009年美国勒纳研究所的Yin Fang等比较健康人和中风病人在手部运动中的脑肌电同步情况，发现中风病人在gamma频段同步缺失。相干分析已成为神经科学领域估计神经功能联系的经典方法，它能够有效的估计信号频域之间的线性相关性，而神经之间有可能存在非线性联系。为此，一些学者将信息论的方法引入脑肌电同步分析，最初是基于互信息表征两个变量之间的关系，可以发现二阶或更高阶的相关性。由于传统的互信息方法无方向性且不能排除同源信号的影响。本发明作者曾提出了基于复杂度信息熵的信息传递指数，表征序列在不同特征空间中的能量分布不确定性。该方法最初提出是用于机械故障诊断，对于运动神经系统控制反馈机理的研究尚存在以下不足：通过计算两个信源之间的共有信息量来描述信息传递关系，不考虑特定信息出现的顺序和时间，本质上没有方向分辨力，且不能排除同源信息的干扰。

发明内容

本发明的目的在于提出基于时序局部互信息的状态转移率，能发现线性和非线性的信息相关性、有助于深入探索大脑皮层与肌肉之间的信息传递特征的脑肌电同步采集及信息传递特性分析方法。

本发明所述的方法包括信号同步采集部分和信号处理部分；

（1）信号同步采集部分包括：脑电信号采集和肌电信号采集；

（1-1）脑电信号采集，采用国际标准10-20电极放置标准，通过电极帽实现电极与头皮接触，采用单级导联法，A1、A2导联分别连接到左右耳后乳突作为参考电极，接地电极布置在头顶正中；

（1-2）肌电信号采集，利用三点式差分输入，两个信号为肌电的差分输入端，另一个为参考地；差分输入电极顺着肌纤维方向，放置在肌腹位置处；

利用32通道脑电采集仪采集脑电信号，修改其中两个差分通道通带宽度用于采集肌电信号，进而实现脑肌电同步采集；

（2）信号处理部分包括信号预处理和脑肌电信息传递特性分析方法；

（2-1）信号预处理：利用自适应高通滤波器和自适应50Hz工频陷波滤波器分别对脑电信号和肌电信号进行滤波处理，去除信号中的基线漂移和工频干扰；然后选用巴特沃斯三阶带通FIR滤波器对脑电信号和肌电信号进行处理，根据信号的有效频段特征，选取脑电信号的截止频率为：2Hz和50Hz；肌电信号的截止频率为：2Hz和200Hz；

（2-2）脑肌电信息传递特性分析方法

利用状态转移率分析脑肌电信号信息传递特性，定量描述脑电信号与肌电信号之间线性和非线性的信息传递；

状态转移率的定义如下：

${\mathrm{STR}}_{x\→y}=\frac{H\left(y_{t+u}\right|y_t^n)-H\left(y_{t+u}\right|x_t^m,y_t^n)}{H\left(y_{t+u}\right|y_t^n)}\∈\left[\mathrm{0,1}\right]---\left(1\right)$

式中，x为源时间序列；y为目标时间序列；STR_x→y为x到y的状态转移率； 分别为序列x和y在t时刻的状态向量；y_t+u为序列y在t+u时刻的状态；

为目标序列y在自身历史状态和源序列x历史状态下的条件熵；根据条件熵的定义可以由（2）式计算

$H\left(y_{t+u}\right|y_t^n)=H(y_{t+u},y_t^n)-H\left(y_t^n\right)---\left(2\right)$

式中，为目标序列y相对于自身历史状态新增加的信息量；为序列y在t+u时刻状态和自身历史状态的总信息量；为序列y历史状态包含的信息量；

同理计算如式（3）

$H\left(y_{t+u}\right|y_t^n{,x}_t^m)=H(y_{t+u},y_t^n,x_t^m)-H(y_t^n,x_t^m)---\left(3\right)$

为目标序列y相对于源序列x与自身历史状态新增加的信息量，两者之差为目标序列新增信息量中来自于源序列历史状态的信息量，即为源序列到目标序列的信息转移量；该信息转移量与目标序列新增信息量的比值即为状态转移率；

所述的状态转移率可由下式得出：

${\mathrm{STR}}_{x\→y}=\frac{H\left(y_{t+u}{,y}_t^n\right)+H(x_t^m,y_t^n)-H(y_{t+u},x_t^m,y_t^n)}{H(y_{t+u},y_t^n)-H\left(y_t^n\right)}---\left(4\right)$

式中，为序列y在t+u时刻状态和自身历史状态的联合熵；为x和y的历史状态的联合熵；为序列y历史状态的熵；为序列y在t+u时刻状态和x、y的历史状态的联合熵；

如果序列y未来某个时刻的状态完全由自身历史状态决定，即则x到y的状态转移率等于0；

如果序列y不能完全由自身的历史状态确定，但可由x和y的历史状态完全确定，则而这种情况下x到y的状态转移率为1(前提是假设y不能由自身的历史状态完全确定)。

为进一步体现不同测点之间的信息传递特征，以不同通道之间的STR作为矩阵元素构建状态转移矩阵，矩阵如下：

与已有技术相比，本发明的有益效果为：

本发明同步采集人体运动中相关脑肌电信号，利用状态转移率分析脑肌电信号信息传递特性，定量描述脑电信号与肌电信号之间线性和非线性的信息传递，有助于探索大脑皮层与肌肉之间的功能联系，研究运动控制反馈机制及运动障碍病理机制，建立基于脑肌电信号的康复状态评价指标，构建康复机器人运动状态及患者生理状态评价机制，可以获得可观的社会效益和经济效益。

附图说明

图1为本发明的脑肌电同步采集系统结构简图。

图2为本发明的工作流程图。

图3为本发明同步采集脑肌电信号的变化结果。

图4为中风患者与健康人膝关节屈伸运动中脑肌电状态转移率的结果对比图。

图5为中风患者与健康人双侧肢体运动时脑肌电状态转移率的差值对比图。

附图标号：1-脑电电极、2-电极帽、3-脑肌电采集仪、4-肌电导联线、5-肌电电极。

具体实施方式

人体生理电信号随着外界刺激和自身状态的变化而改变，导致了生理电信号的非平稳性，使得脑电信号和肌电信号具有非线性的特点。由于信息熵能够刻画信号间线性和非线性的信息相关性，本发明通过计算信号新增状态与源信号历史状态之间的共有信息量，获得不同运动状态下大脑皮层与肌肉之间信息传递关系，进而研究运动功能障碍产生的生理机制。具体方法包括信号同步采集部分和信号处理部分；

（1）信号同步采集部分包括：脑电信号采集和肌电信号采集；

（1-1）脑电信号采集，采用国际标准10-20电极放置标准，通过电极帽2实现脑电电极1与头皮接触。进行膝关节屈伸运动下脑肌电同步采集实验，选择C1、C2以及F3、F4区进行采集，采用单级导联法，A1、A2导联分别连接到左右耳后乳突作为参考电极，接地电极布置在头顶正中；

（1-2）肌电信号采集，利用三点式差分输入，两个信号为肌电的差分输入端，另一个为参考地；差分输入电极顺着股四头肌肌纤维方向，放置在肌腹位置处，使用酒精擦拭被测部位的皮肤表面，去除皮肤表面油脂和皮屑，粘贴肌电电极5，并将肌电导联线4适当固定尽量减少运动过程中导联线晃动的干扰。

利用32通道脑肌电采集仪3采集脑电信号，修改其中两个差分通道通带宽度用于采集肌电信号，进而实现脑肌电同步采集；

信号预处理：脑肌电信号均属于微弱信号，容易受到噪声干扰，采集到的原始数据需要进行预处理，利用自适应高通滤波器和自适应50Hz工频陷波滤波器分别对脑电信号和肌电信号进行滤波处理，去除信号中的基线漂移和工频干扰；然后选用巴特沃斯三阶带通FIR滤波器对脑电信号和肌电信号进行处理，根据信号的有效频段特征，选取脑电信号的截止频率为：2Hz和50Hz；肌电信号的截止频率为：2Hz和200Hz；

脑肌电信号状态转移率计算

由于采集到的信号只是大脑皮层或肌肉这种复杂系统的一个观测变量输出，而有意义的因果性假设是基于系统状态变化，而不是实际采集的信号。因此，对每个通道采集的时间序列进行状态空间重构，实际采集到的长度为N的时间序列x₁,x₂,…,x_n-1,x_n，采用Takens延迟坐标法构造状态向量：

$x_t^d=(x_t,x_{t-\mathrm{\τ}},x_{t-2\mathrm{\τ}},...,x_{t-(d-1)\mathrm{\τ}})$

式中，t=1,2,…,N-(d-1)τ，延迟时间τ由自相关法确定，取自相关函数下降到初始值的1-1/e时的τ为延迟时间，同时嵌入维数d由Cao方法确定。

按下式计算脑肌电通道间的状态转移率

${\mathrm{STR}}_{x\→y}=\frac{H\left(y_{t+u}\right|y_t^n)-H\left(y_{t+u}\right|x_t^m,y_t^n)}{H\left(y_{t+u}\right|y_t^n)}\∈\left[\mathrm{0,1}\right]---\left(1\right)$

式中，x为源时间序列；y为目标时间序列；STR_x→y为x到y的状态转移率； 分别为序列x和y在t时刻的状态向量；y_t+u为序列y在t+u时刻的状态；为目标序列y在自身历史状态下的条件熵；为目标序列y在自身历史状态和源序列x历史状态下的条件熵；根据条件熵的定义，可以由下式计算

$H\left(y_{t+u}\right|y_t^n)=H(y_{t+u},y_t^n)-H\left(y_t^n\right)---\left(2\right)$

同理计算 $H\left(y_{t+u}\right|y_t^n,x_t^m),$

$H\left(y_{t+u}\right|y_t^n{,x}_t^m)=H(y_{t+u},y_t^n,x_t^m)-H(y_t^n,x_t^m)---\left(3\right)$

则状态转移率可由下式计算得到

${\mathrm{STR}}_{x\→y}=\frac{H\left(y_{t+u}{,y}_t^n\right)+H(x_t^m,y_t^n)-H(y_{t+u},x_t^m,y_t^n)}{H(y_{t+u},y_t^n)-H\left(y_t^n\right)}---\left(4\right)$

为进一步体现不同测点之间的信息传递特征，以不同通道之间的STR作为矩阵元素构建状态转移矩阵，矩阵如下：

为验证本发明所述的脑肌电信号信息传递特性分析方法的可行性和有效性，与秦皇岛市第一医院合作，进行6名中风引起的右侧运动功能障碍患者的测试实验，如表1所示。按照本发明所述的脑肌电采集与分析过程，同步采集患者膝关节屈伸运动下脑肌电信号，并分析与健康人信息传递特性的差异，得到患者运动神经环路状态的评价。

表1患者相关信息

本次实验采集股四头肌表面肌电信号与对侧两通道脑电信号，并计算状态转移率。

图4为中风患者与健康人膝关节屈伸运动中脑肌电状态转移率的结果对比图，图中SR代表患者右侧肢体运动时脑肌电间状态转移率，SL代表患者左侧，CR和CL分别代表健康人的右侧和左侧。对比分析多组数据，中风患者右侧（患侧)运动过程中脑肌电状态转移率平均在0.04左右明显低于左侧的0.06左右、低于健康人（平均在0.07到0.08之间)。

图5为中风患者与健康人双侧肢体运动时脑肌电状态转移率的差值对比图。从图中可清晰判断分析出患者的两侧肢体运动情况，为医治提供有效数据。

本发明提供了用于脑肌电信息传递特性分析的新方法，能够定量描述脑电信号与肌电信号之间线性和非线性的信息传递，为深入探索运动功能障碍产生机制提供有效的观察手段，同时为运动功能康复状态和康复手段的有效性提供新的评价指标，可以广泛应用于康复工程、神经工程和运动机制研究领域。

Claims (2)

1.一种脑肌电同步采集及信息传递特性分析方法，其特征在于：包括信号同步采集部分和信号处理部分；

(1)信号同步采集部分包括：脑电信号采集和肌电信号采集；

(1-1)脑电信号采集，采用国际标准10-20电极放置标准，通过电极帽实现电极与头皮接触，采用单级导联法，A1、A2导联分别连接到左右耳后乳突作为参考电极，接地电极布置在头顶正中；

(1-2)肌电信号采集，利用三点式差分输入，两个信号为肌电的差分输入端，另一个为参考地；差分输入电极顺着肌纤维方向，放置在肌腹位置处；

利用32通道脑电采集仪采集脑电信号，修改其中两个差分通道通带宽度用于采集肌电信号，进而实现脑肌电同步采集；

(2)信号处理部分包括信号预处理和脑肌电信息传递特性分析方法；

(2-1)信号预处理：利用自适应高通滤波器和自适应50Hz工频陷波滤波器分别对脑电信号和肌电信号进行滤波处理，去除信号中的基线漂移和工频干扰；然后选用巴特沃斯三阶带通FIR滤波器对脑电信号和肌电信号进行处理，根据信号的有效频段特征，设定高通截止频率为2Hz及低通截止频率为50Hz的FIR带通滤波器对脑电信号进行有效滤波；设定高通截止频率为2Hz及低通截止频率为200Hz的FIR带通滤波器对肌电信号进行滤波；

(2-2)脑肌电信息传递特性分析方法

利用状态转移率分析脑肌电信号信息传递特性，定量描述脑电信号与肌电信号之间线性和非线性的信息传递；

状态转移率的定义如下：

${\mathrm{STR}}_{x\→y}=\frac{H\left(y_{t+u}\right|y_t^n)-H\left(y_{t+u}\right|x_t^m,y_t^n)}{H\left(y_{t+u}\right|y_t^n)}\∈\left[\mathrm{0,1}\right]---\left(1\right)$

式中，x为源时间序列；y为目标时间序列；STR_x→y为x到y的状态转移率；分别为序列x和y在t时刻的状态向量；y_t+u为序列y在t+u时刻的状态；为目标序列y在自身历史状态下的条件熵；为目标序列y在自身历史状态和源序列x历史状态下的条件熵；根据条件熵的定义，可以由(2)式计算

$H\left(y_{t+u}\right|y_t^n)=H(y_{t+u},y_t^n)-H\left(y_t^n\right)---\left(2\right)$

式中，为目标序列y相对于自身历史状态新增加的信息量；为序列y在t+u时刻状态和自身历史状态的总信息量；为序列y历史状态包含的信息量；

同理计算如式(3)

$H\left(y_{t+u}\right|y_t^n,x_t^m)=H(y_{t+u},y_t^n,x_t^m)-H(y_t^n,x_t^m)---\left(3\right)$

为目标序列y相对于源序列x与自身历史状态新增加的信息量，两者之差为目标序列新增信息量中来自于源序列历史状态的信息量，即为源序列到目标序列的信息转移量；该信息转移量与目标序列新增信息量的比值即为状态转移率；

所述的状态转移率可由下式得出：

${\mathrm{STR}}_{x\→y}=\frac{H(y_{t+u},y_t^n)+H(x_t^m,y_t^n)-H\left(y_t^n\right)-H(y_{t+u},x_t^m,y_t^n)}{H(y_{t+u},y_t^n)-H\left(y_t^n\right)}---\left(4\right)$

式中，为序列y在t+u时刻状态和自身历史状态的联合熵；为x和y的历史状态的联合熵；为序列y历史状态的熵；为序列y在t+u时刻状态和x、y的历史状态的联合熵；

如果序列y未来某个时刻的状态完全由自身历史状态决定，即则x到y的状态转移率等于0；

如果序列y不能完全由自身的历史状态确定，但可由x和y的历史状态完全确定，则而这种情况下x到y的状态转移率为1。

2.根据权利要求1所述的脑肌电同步采集及信息传递特性分析方法，其特征在于：为进一步体现不同测点之间的信息传递特征，以不同通道之间的STR作为矩阵元素构建状态转移矩阵，矩阵如下：