基于变尺度符号传递熵的多通道脑肌电耦合分析方法与流程

本发明属于生物信号处理领域，涉及一种大脑皮层肌肉和运动肌肉功能耦合的分析方法。

背景技术：

脑电信号(Electroencephalogram,EEG)是通过电极记录到的大脑神经群自发或者节律性活动所产生的信号，它反映了大脑功能区神经细胞群活动时的电位变化。肌电信号(Electromyography,EMG)是通过电极记录肌群活动时所产生的信号。人的自主动作过程中，除了不同脑区之间的协同作用，相关的脑区和肌肉组织相关的周围神经之间也自动地发生了同步化，这种大脑皮层与运动神经肌肉组织两大体系内部及相互间的作用构成了皮层肌肉功能耦合关系(Cortical-Muscular Functional Coupling，CMC)。脑肌电耦合分析在近年来成为许多学者研究的热点，为脑卒中和运动障碍患者提供康复方案和康复状态分析的理论依据。

大脑是一个相当复杂的非线性动力系统，脑电信号的变化和神经中枢兴奋的传导都有一定的非线性特征。20世纪90年代起，非线性动力学的有关方法如分维数、Lorenz散点图、Lyapunov指数、复杂度等被应用到脑电信号的研究中，并取得了一定的进展。近年来相干性分析、格兰杰因果分析、传递熵以及符号化传递熵等方法都被用来分析脑肌电耦合关系。其中符号传递熵研究脑肌电耦合性的关键是符号化后的序列对信息动态特性的保留程度。但是，传统的时间序列符号化方法用同一个符号表示绝对值相当的正负两个值，容易造成对原始信号的混淆现象，丢失信号的动态特性；另一方面，符号化参数固定也会损失信号的部分特征。针对这些问题，提出一种可变尺度参数的符号化方法，并结合传递熵对上肢自主动作中不同握力输出时的脑肌电信号进行了多通道分析，进而提出一种耦合强度的表示方法，对不同握力下脑肌电耦合强度进行定量分析，为运动功能康复状态分析提供了客观依据。

技术实现要素：

本发明为了客观、定量、有效地对大脑皮层与运动神经肌肉组织的耦合性进行分析，本发明利用上肢主动动作时的头皮脑电信号和表面肌电信号，提出了一种基于变尺度符号传递熵的多通道脑肌电耦合分析方法。首先同步采集不同握力动作时多通道脑电信号和相关肌肉组上的表面肌电信号。然后选择可变尺度参数对左、右手相同握力下的脑肌电信号进行符号化，并对符号化后的序列进行传递熵的计算。通过综合分析传递熵的平均值和标准差，以及计算所用时间，选择一个合适有效的符号化尺度参数进行后续分析。进一步对左手/右手动作、不同握力下、多个通道的脑电信号和对应的肌电信号进行符号传递熵分析和比较。最后根据脑电到肌电(EEG—>EMG)和肌电到脑电(EMG—>EEG)的传递熵变化情况，提出了脑肌电信号耦合强度的表示方法，客观定量的反映皮层肌肉和运动肌肉之间的耦合强度。

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

步骤一，获取人体上肢动作时的头皮脑电信号和表面肌电信号样本数据。

步骤二，将获取的EEG和EMG信号设置多个尺度参数进行变尺度符号化处理。

可变尺度符号化方法具体如下：

可变尺度参数的符号化方法，它的符号个数与尺度有关，可动态进行调整，具体思路如下：

1.首先求出序列的最大值和最小值；

2.设定一个可变符号化尺度的参数pieces，它表示将序列最小值到最大值划分为pieces+1份，最终的符号个数也就为pieces+1，pieces小于时间序列的长度，它的值大则表示划分精细，取值小则表示划分粗糙；

3.然后对时间序列符号化，落入最小的区间设置符号为-pieces/2，依次为-pieces/2+0.5，依次类推，最大的符号为pieces/2。

具体函数形式如(1)式所示：

其中i表示时间序列的长度，S(i)代表符号化后的序列，min(x)和max(x)分别代表时间序列的最小和最大值，delta代表每次区间的增加量，其取值为

步骤三，将符号化的脑肌电信号符号序列进行传递熵的计算，取多次的平均值、标准差和计算时间。选取EEG—>EMG传递熵的平均值介于0.40～0.60之间，EMG—>EEG传递熵的平均值介于0.15～0.30之间，标准差小于0.008，一次计算时间少于3.0s时的尺度参数进行后续的分析。

传递熵的计算方法如下：

假如给定两个时间序列X＝{x₁,x₂,···,x_T}和Y＝{y₁,y₂,···,y_T}，其中T是时间序列的长度，x₁、y₁分别是第一个观测值，x₂、y₂分别是第二个观测值，依次类推。我们可以得到Y到X的传递熵(TE_Y→X)和X到Y的传递熵(TE_X→Y)如(2)、(3)式：

其中n为离散时间指标，τ为预测时间，p(·)代表概率分布。

步骤四，利用选定的尺度参数对多个通道脑肌电信号进行符号化，并计算传递熵。

步骤五，提出脑肌电信号耦合强度的表示方法，对皮层肌肉耦合性进行定量分析。

耦合强度的表示方法：

在分析脑肌电传递熵时发现大脑左右半球分界线上的FZ和CZ通道脑肌电传递熵在左、右手握拳时均没有明显变化。因此，为了定量表示明显变化通道符号化传递熵的变化特征以及分析脑肌电的耦合强度，本发明提出以FZ和CZ通道脑肌电传递熵的平均值作为基线，以明显变化通道的传递熵值与基线的差值进行表示，其具体定义如式(4)和(5)：

S_{TE(EEG-＞EMG)}＝TE_{(EEG-＞EMG)i}-Mean_(EEG-＞EMG) (4)

S_{TE(EMG-＞EEG)}＝TE_{(EMG-＞EEG)i}-Mean_(EMG-＞EEG) (5)

其中S_{TE(EEG-＞EMG)}表示，脑电到肌电的耦合强度，S_{TE(EMG-＞EEG)}表示，肌电到脑电的耦合强度；TE_{(EEG-＞EMG)i}表示，脑电到肌电变化明显通道i的传递熵，TE_{(EMG-＞EEG)i}表示，肌电到脑电变化明显通道i的传递熵；Mean_(EEG-＞EMG)表示FZ和CZ通道脑电到肌电的传递熵，Mean_(EMG-＞EEG)表示FZ和CZ通道肌电到脑电的传递熵。

本发明与已有的诸多脑肌电耦合分析方法相比，具有如下特点：

在信号源的选择上，脑肌电的耦合强度可以反映皮层和运动肌肉之间的联系，本发明选择了32通道脑电信号和2通道的肌电信号：脑电信号是大脑功能区神经细胞群活动时的电位变化，体现了大脑功能区的活动状态；肌电信号是肌群活动时的电位变化，反映了肌肉的运动情况。并且上肢运动在大脑皮层占了最大的面积，因此，选择脑电信号和上肢的肌电为脑肌电耦合研究提供了综合、客观、可靠的信号源。

在分析方法上，变尺度符号化方法更能准确的捕捉信息的动态特性，根据信号的特点，可以连续的改变符号化尺度参数，从而选择合适的尺度参数进行时间序列的符号化，因此本发明能适应复杂的信号，抗干扰能力强，算法简单。传递熵对信息的传递方向和能量可以定量的表示。耦合强度的表示客观地体现了大脑运动皮层对运动肌肉的控制能力以及肌肉状态对大脑的反馈。这些为运动障碍患者以及脑卒中病人的康复状态分析提供了理论依据。

附图说明

图1(a)所示为本发明的肌电传感器放置位置，1(b)所示为脑电传感器放置位置；

图2(a)所示为本发明对左手10次动作脑肌电信号不同尺度符号化传递熵平均值的计算结果

图2(b)所示为本发明对左手10次动作脑肌电信号不同尺度符号化传递熵标准差的计算结果

图2(c)所示为本发明对右手10次动作脑肌电信号不同尺度符号化传递熵平均值的计算结果

图2(d)所示为本发明对右手10次动作脑肌电信号不同尺度符号化传递熵标准差的计算结果

图2(e)所示为本发明对10次动作脑肌电信号不同尺度符号化传递熵计算所用时间；

图3-1-1所示为本发明对受试者S2左手5kg握力动作脑肌电信号求符号传递熵的分析结果

图3-1-2所示为本发明对受试者S2右手5kg握力动作脑肌电信号求符号传递熵的分析结果

图3-2-1所示为本发明对受试者S2左手10kg握力动作脑肌电信号求符号传递熵的分析结果

图3-2-2所示为本发明对受试者S2右手10kg握力动作脑肌电信号求符号传递熵的分析结果

图3-3-1所示为本发明对受试者S2左手20kg握力动作脑肌电信号求符号传递熵的分析结果

图3-3-2所示为本发明对受试者S2右手20kg握力动作脑肌电信号求符号传递熵的分析结果

图4-1-1所示为本发明对受试者S1左手5kg、10kg、20kg握力下脑肌电耦合强度分析结果

图4-1-2所示为本发明对受试者S1右手5kg、10kg、20kg握力下脑肌电耦合强度分析结果

图4-2-1所示为本发明对受试者S2左手5kg、10kg、20kg握力下脑肌电耦合强度分析结果

图4-2-2所示为本发明对受试者S2右手5kg、10kg、20kg握力下脑肌电耦合强度分析结果

图4-3-1所示为本发明对受试者S3左手5kg、10kg、20kg握力下脑肌电耦合强度的分析结果

图4-3-2所示为本发明对受试者S3右手5kg、10kg、20kg握力下脑肌电耦合强度分析结果

具体实施方式

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

本实施例包括如下步骤：

1.基于变尺度符号传递熵的多通道脑肌电耦合分析，其具体实施方式包括如下步骤：

步骤一，获取人体上肢动作时的头皮脑电信号和表面肌电信号样本数据，如图1(a)和1(b)所示。

脑电信号反映大脑功能区神经细胞群活动时的电位变化，肌电信号反映肌群活动时的电位变化。人的自主动作过程中，除了不同脑区之间的协同作用，相关的脑区和肌肉组织相关的周围神经之间也自动地发生了同步化。这种大脑皮层与运动神经肌肉组织两大体系内部及相互间的作用构成了皮层肌肉功能耦合关系。因此，本发明采集上肢动作时脑肌电信号对这种耦合关系进行定量分析。

同步采集人体上肢动作时的脑肌电信号。通过128导BrainAmp DC脑电采集系统(Brain Products GmbH,Germany)同步拾取了32通道脑电和2通道肌电信号。选用32通道脑电信号(FP1，FP2，F7，F8，F4，F3，FZ，FC5，FC1，FC2，FC6，T7，C3，CZ，C4，T8，CP5，CP1，CP2，CP6，TP9，P7，P3，PZ，P4，P8，TP10，PO9，O1，OZ，O2，PO10)和左右手臂指浅屈肌(flexor digitorum superficialis,FDS)和掌长肌(palmaris longus,PL)的肌电信号表面肌电为信号来源。受试者(S1、S2、S3)左、右手分别使用弹簧握力计进行5kg、10kg和20kg力输出的握拳实验各10次，共180次。实验前头皮需清洗干净，先用酒精分别在受试者的手臂指浅屈肌和掌长肌上擦涂去污，以增强信号拾取能力。

步骤二，将获取的EEG和EMG信号进行变尺度符号化，其尺度参数设置从20到50以5为步长进行变化，进行符号化处理。

符号化是对生理信号粗粒化的过程，这一过程能够捕捉大尺度的动态特征，从而降低动力学噪声和测量噪声的影响。而时间序列的符号化精确程度影响着传递熵计算的精确性和系统的动态耦合指标。因此，引进变尺度符号函数对时间序列进行符号化。根据信号的复杂程度设置尺度参数从20到50以5为步长进行变化，依据符号函数进行符号化。

步骤三，将符号化的脑肌电信号符号序列进行传递熵的计算，取每10次的平均值、标准差和计算时间，综合考虑平均值、标准差和计算时间指标，选取合适的尺度参数进行后续的分析。

结合不同尺度符号化传递熵的平均值、标准差和计算所用时间分析，每10次的平均值如图2(a)和2(c)所示、标准差如图2(b)和2(d)所示、计算时间如图2(e)所示，发现随着尺度增大标准差、计算时间也快速增大，这说明对于同一个动作计算的传递熵波动变大、耗时更多。对其它几组握力的脑肌电信号进行多尺度符号传递熵的分析，也得到类似的结果。因此，综合考虑传递熵的平均值、标准差和计算时间，选取尺度参数为30，对时间序列符号化，并进行脑肌电信号耦合性分析。

步骤四，利用选定的尺度参数对脑肌电信号进行符号化，并计算传递熵。

为了确定不同握力输出时，大脑运动区各通道脑电与上肢肌电信号的耦合情况，利用步骤三中选定的符号化尺度参数，对三位受试者5kg、10kg和20kg等不同握力下、不同通道的脑电信号、左手/右手动作多方位地对头皮脑电信号和表面肌电信号进行符号化传递熵计算和分析，其中S2的计算结果如图3-1-1～3-3-2所示。

步骤五，提出脑肌电信号耦合强度的表示方法，对皮层肌肉耦合性进行定量分析。

受试者S1的脑肌电耦合强度的分析结果如图4-1-1和4-1-2所示，受试者S2的脑肌电耦合强度的分析结果如图4-2-1和4-2-2所示，受试者S3的脑肌电耦合强度的分析结果如图4-3-1和4-3-2所示，通过分析三位受试者脑肌电耦合强度发现，左手握力状态下耦合性最强的是C4通道脑电，右手握力状态下耦合性最强的是C3通道脑电，不同通道间有显著的差异。因此以大脑左右半球分界线上的FZ和CZ通道脑肌电变尺度符号传递熵的平均值作为基线，以明显变化通道的传递熵值与基线的差值，提出左/右手在不同握力下皮层肌肉耦合强度的表示方法，从而客观定量的表示皮层肌肉和运动肌肉之间的耦合强度，反映大脑皮层对运动肌肉的控制能力以及肌肉状态对大脑的反馈。