CN110367974B - 基于变分模态分解-传递熵的脑肌电耦合的研究方法 - Google Patents

Abstract

本发明提出了信息技术领域一种基于变分模态分解‑传递熵的脑肌电耦合的研究方法，包括以下步骤：步骤1受试者按照实验指令完成实验操作，在实验过程中采集脑电信号和肌电信号；步骤2对上述脑电信号、肌电信号数据的预处理；步骤3对预处理过后的脑电信号和肌电信号进行变分模态分解，将脑电信号和肌电信号分解为若干个不同的本征模态函数；步骤4对使用变分模态分解得到的脑电信号、肌电信号两两之间进行传递熵计算，观察脑肌电耦合强度；步骤5选取不同的时间刻度再次对不同性质、不同分量、不同方向的本征模态函数进行传递熵计算，可以减少频带混叠给后续不同频率波之间进行传递熵计算产生的干扰，提高实验结果准确性。

Description

基于变分模态分解-传递熵的脑肌电耦合的研究方法

技术领域

本发明涉及一种脑肌电耦合研究方法，具体的说是一种基于变分模态分解-传递熵的脑肌电耦合的研究方法，属于认知神经科学和信息技术领域。

背景技术

人体的自主运动是大脑运动皮层发出运动指令，经过中枢神经系统(脑、脊髓)，沿着运动传导通路支配肢体完成一定的动作；同时，人的肢体基于周边环境产生体感，体感沿着神经传导通路反馈给大脑，在大脑的感觉区和运动区加以综合判断，调节运动指令，从而准确的完成动作，这种信号在大脑和肌肉之间的来回传递，表明来自大脑的EEG信号和来自肌肉的EMG信号之间具有一定的相关性，因此我们可以对这种相干性进行研究，这种研究具有一定的实际意义，可以将这种相干性研究用于皮质肌肉耦合、肌肉疲劳、运动障碍、运动康复等一系列的分析中。

生理系统是跨越多尺度时空进行相互作用的复杂系统。运动神经系统通过神经振荡传递运动控制信息，引起相应肌肉的运动单元的同步性振荡活动反映运动响应信息，这种同步振荡活动可以反应多层次的皮质-肌肉功能耦合的连接信息。由于脑电和肌电信号分别包含躯体运动控制信息和肌肉对大脑控制意图的响应信息，因此对运动模式下脑肌电信号间不同时空层次的同步特征分析更有助于理解皮质-肌肉的功能耦合及信息传递。

研究多时空层次的脑肌电耦合信息，关键因素是提取不同尺度上的脑电和肌电信息；首先是粗粒化方法和近似熵结合应用于EEG和 EMG的尺度特性分析中，提取不同时间尺度上的信息；但随着其尺度的增加，序列长度减少使得熵估计不准确；移动均值化的提出克服了这一弊端，使得各个尺度的时间序列长度保持相同；然而粗粒化与移动均值化方法仅是将信号进行时间尺度化，不能刻画EEG和EMG的时频域特性，因此将经验模态分解的方法应用于多尺度划分中，提取 EEG和EMG信号不同时频域尺度特征；EMD方法自适应的将信号分解成一系列的本征模态函数，得到各个本征模态函数的瞬时幅度和频率；但是经过该方法分解之后得到的本征模态函数存在边界效应和模态混叠的现象；因此集总经验模态分解方法应用到脑电信号分析中，消除了模态混叠的现象，相比EMD具有更好的抗噪声性能。但是采用 EEMD方法会造成计算量增大。鉴于以上递归式分解的不足，将变分模态分解引入到尺度化分析中，自使用的提取EEG和EMG中各个具有中心频率的窄带分量；为了定量的研究不同时频尺度之间的脑肌电非线性同步耦合特征及信息传递，本专利结合变分模态分解和传递熵的方法，构建变分模态分解-传递熵模型用于脑肌电频段间耦合性分析。

发明内容

本发明的目的是提供一种基于变分模态分解-传递熵的脑肌电耦合的研究方法，可以减少频带混叠给后续不同频率波之间进行传递熵计算产生的干扰，提高实验结果准确性。

本发明的目的是这样实现的：一种基于变分模态分解-传递熵的脑肌电耦合的研究方法，包括以下步骤：

步骤1)受试者按照实验指令完成实验操作，在实验过程中采集脑电信号和肌电信号，并以数字信号存储在计算机中；

步骤2)对上述脑电信号、肌电信号数据的预处理；

步骤3)对预处理过后的脑电信号和肌电信号进行变分模态分解，将脑电信号和肌电信号分解为若干个不同的本征模态函数；

步骤4)对使用变分模态分解得到的脑电信号、肌电信号两两之间进行传递熵计算，观察脑肌电耦合强度；

步骤5)选取不同的时间刻度再次对不同性质、不同分量、不同方向的本征模态函数进行传递熵计算，观察脑肌电耦合随负荷大小和时间变化的关系。

作为本发明的进一步限定，步骤3)的具体算法如下：

变分模态分解方法是将信号分解过程转移到变分框架内，通过搜寻约束变分模型的最优解实现信号自适应分解；首先将每个模态重新定义为调频-调幅信号，表达为：

u_k(t)＝A_k(t)cos(Φ_k(t)) (1)

式中，A_k(t)为u_k(t)的瞬时幅值，

ω_k(t)为 u_k(t)的瞬时频率，假设每个模态u_k(t)具有中心频率和有限带宽，约束条件为各模态之和等于输入信号f，且每个模态的估计带宽之和最小，模型如下：

st.∑_ku_k＝f (2)

其中，u_k＝{u₁,...,u_k}，表示分解的K有限带宽IMFf分量，为求取(1)、(2)约束变分问题的最优解，引入增广形式的Largrange函数，即：

式中：α为惩罚因子，λ为拉格朗日算子；

采用乘法算子交替方向法解决以上变分问题，通过交替更新u_k ⁿ⁺¹、ω_k ⁿ⁺¹和λⁿ⁺¹，寻找上式的“鞍点”；

其中，

的表达式为：

式中，ω_k等同于

等同于

基于Parseval/Plancherel傅里叶等距变换，将(4)式转变到频域，得到各模态的频域更新；然后将中心频率的取值问题转换到频域，得到中心频率的更新方法；同时更新λ；

具体表达式如下：

每个本征模态函数分量的中心频率和带宽在迭代求解变分模型的过程中不断更新，直至满足迭代停止条件

对于给定的判定精度e＞0，结束整个循环，最终根据实际信号的频域特性得到K个窄带本征模态函数分量。

作为本发明的进一步限定，步骤(4)将经过变分模态分解得到的若干个脑电、肌电信号的本征模态函数两两之间计算传递熵值，具体算法如下：

将使用变分模态分解得到的脑电信号分量和肌电信号分量分别定义为

K₁和K₂分别代表EEG和EMG窄带分量本征模态函数的个数，基于传递熵的定义，构造

到

的变分模态分解-传递熵TE_EEG→EMG，公式如下：

式中，t离散时间指标；u为预测时间；

分别为k₁，k₂的延迟向量，p(·)为变量间的联合概率；

TE_EEG→EMG则表示EEG的k₁分量到EMG的k₂分量间的传递熵值，同理信号

到

的变分模态分解-传递熵TE_EMG→EEG的表达式为：

TE_EMG→EEG则表示EMG的k₂分量到EEG的k₁分量间的传递熵值。

作为本发明的进一步限定，步骤5)的时间刻度分别以1s、2s、 3s、4s为基准，将整个时间段的本征模态函数分割再加以计算，并绘制传递熵变化曲线，观察随着时间变化和负荷逐渐增加时，皮质肌肉耦合的变化情况。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

第一：通过变分模态分解的方法，将经过预处理后的EEG和EMG 信号分别分解成为三个无明显频带混叠的不同模态信号；可以减少频带混叠给后续不同频率波之间进行传递熵计算产生的干扰，提高实验结果准确性。

第二：对经过变分模态分解得到的不同模态信号，采用计算两两不同类模态之间传递熵数值的方法，评估脑肌电信号的耦合程度，同时进一步采用时间分段的方法，将整个实验过程按照不同时间刻度 (1s、2s、3s、4s)划分，可以观察到：极低负荷时，脑肌电耦合性逐渐增加；在刚开始增加负荷时，身体由于不习惯，传递熵数值下降，脑肌电耦合性出现减弱趋势；在负荷逐渐增加的过程和增加结束后，传递熵数值又出现了增加，并且数值超过刚开始极低负荷时的数值，脑肌电耦合性加强，这种多时间尺度的分析方法定量刻画脑肌电信号间的多时频尺度耦合特征及信息传输方向,为进一步探索运动功能障碍的产生机理及康复过程中的运动功能评价方法提供依据。

附图说明

图1为本发明中变分模态分解-传递熵分析脑肌电耦合框架图。

图2为本发明中脑电信号的通道导联示意图。

图3为本发明中肌电信号采集的位置。

图4为本发明中原始脑电信号。

图5为本发明中预处理后的脑电信号。

图6为本发明中原始肌电信号。

图7为本发明中预处理后的肌电信号。

图8为本发明中经过变分模态分解后的脑电模态。

图9为本发明中经过变分模态分解后的脑电各个模态的频谱图。

图10为本发明中经过变分模态分解后的肌电模态。

图11为本发明中经过变分模态分解后的肌电各个模态的频谱图。

图12为整段时间内EEG→EMG传递熵值地形图。

图13为整段时间内EMG→EEG传递熵值地形图。

图14为时间尺度为1s时EEG→EMG传递熵值变化曲线。

图15为时间尺度为2s时EEG→EMG传递熵值变化曲线。

图16为时间尺度为3s时EEG→EMG传递熵值变化曲线。

图17为时间尺度为4s时EEG→EMG传递熵值变化曲线。

图18为时间尺度为1s时EMG→EEG传递熵值变化曲线。

图19为时间尺度为2s时EMG→EEG传递熵值变化曲线。

图20为时间尺度为3s时EMG→EEG传递熵值变化曲线。

图21为时间尺度为4s时EMG→EEG传递熵值变化曲线。

具体实施方式

为下面结合附图对本发明的技术方案做进一步的详细说明：

如图1所示：该图为本实验的整体方案，实验设计部分为：实验共有十名受试者参与(七名男性，三名女性)，受试者均为在校大学生，右利，且均无精神疾病史，无身体运动功能障碍，实验前均告知实验知情同意书，且均未进行过类似的实验。实验开始前，受试者被要求端坐在指定的座椅并被告知尽量不要做规定以外的动作，实验共进行 40s。实验开始时前10s，受试者使用右手握住一个未注水的空杯，实验进行到第10s时，开始将500ml纯净水匀速的注入到受试者的空杯中，控制注水速度，确保500ml纯净水在第20s时全部注入到受试者手中的空杯内，注水结束后，再握持20s，并记录整个过程中的EEG 和EMG信号。之后对原始EEG和EMG信号进行预处理，去除伪迹干扰，工频干扰，基线漂移等，在对EEG和EMG信号进行相同的带通滤波。预处理过后的EEG和EMG信号送入变分分解框架内，进行变分模态分解，得到EEG和EMG分解后的不同频率分量，之后计算这些模态，两两之间的传递熵，通过传递熵数值的大小来判断皮质肌肉的耦合程度——传递熵数值越大，则表明皮质肌肉耦合程度越大，反之亦然。最后，按照不同的时间尺度绘制传递熵的变化曲线，观察在整个实验过程中皮质肌肉耦合的变化情况。

如图2所示：该图为大脑头皮电极放置方法。本发明采用国际脑电图学会标定的10-20电极导联定位标准，，采用双耳垂接法，注射导电膏以增强电极导电性，选择不同脑区的F3、Fz、F4、C3、Cz、 C4、P3、Pz、P4、O1、Oz、O2，12个通道放置电极，以采集不同脑区电信号，通道A1当做参考电极。由前期研究可知，人体肢体的运动，可以引起对侧大脑皮层相应运动区域的信号变化。本实验中，受试者均为右利手，因此仅选择C3导联的数据作为EEG信号的原始数据。

如图3所示：该图为肌肉外置电极的放置方法。本发明测量的是屈腕肌和指浅屈肌两块肌肉的电信号。首先使用去角质膏将两块肌肉表面角质去除，降低阻抗，之后将肌电采集设备的两个采样电极分别黏贴在这两块肌肉处，选择手背作为接地位置。

如图4，5所示：图4 为从脑电设备采集的原始脑电数据，肉眼观察存在较多的伪迹成分，使用自适应带通滤波器去除50Hz工频干扰，再使用0.5-70Hz的带通滤波器对原始数据进行预处理.。图5 为预处理后的脑电数据，肉眼可见已经去除了伪迹。

如图6，7所示：图6 为从肌电设备采集的原始肌电数据，肉眼观察存在较多的伪迹成分，使用和脑电原始信号处理相同的0.5-70Hz 的带通滤波器对原始数据进行预处理.。图7 为预处理后的脑电数据，肉眼可见已经去除了伪迹。

如图8所示：该图为将经过预处理的脑电信号送入变分分解框架内，通过搜寻约束变分模型的最优解将预处理过后的脑电信号和肌电信号自适应分解为三个模态。首先将每个模态重新定义为调频-调幅信号，表达为：

u_k(t)＝A_k(t)cos(Φ_k(t))

式中，A_k(t)为u_k(t)的瞬时幅值。

ω_k(t)为u_k(t) 的瞬时频率，假设每个模态u_k(t)具有中心频率和有限带宽，约束条件为各模态之和等于输入信号f，且每个模态的估计带宽之和最小，模型如下：

其中，u_k＝{u₁,...,u_k}，表示分解的K有限带宽IMFf分量。为求取(1)、(2)约束变分问题的最优解，引入增广形式的Largrange函数，即：

式中：α为惩罚因子，λ为拉格朗日算子。

采用乘法算子交替方向法解决以上变分问题，通过交替更新u_k ⁿ⁺¹、ω_k ⁿ⁺¹和λⁿ⁺¹，寻找上式的“鞍点”。

其中，

的表达式为：

式中，ω_k等同于

等同于

基于Parseval/Plancherel傅里叶等距变换，将(4)式转变到频域，得到各模态的频域更新；然后将中心频率的取值问题转换到频域，得到中心频率的更新方法；同时更新λ。

具体表达式如下：

每个IMF分量的中心频率和带宽在迭代求解变分模型的过程中不断更新，直至满足迭代停止条件

对于给定的判定精度e＞0，结束整个循环，最终根据实际信号的频域特性得到K 个窄带IMF分量，完成信号频带的自适应分割，有效避免模态混叠。其中得到的三个脑电模态的中心频率如下表：

模态1	模态2	模态3
			1.2196Hz	21.2488Hz	10.2086Hz

如图9所示：该图是对分解得到的脑电信号的模态，进行FFT得到，分解后模态的频谱图。

如图10，11所示：图10是对EMG信号进行变分模态分解得到的 EMG信号的三种模态的结果。图11是对EMG信号分解出的三种模态进行FFT变换，得到的频谱图。其中得到的三个脑电模态的中心频率如下表：

模态1	模态2	模态3
			2.1875Hz	22.7196Hz	12.7007Hz

如图12，13所示：该图是计算脑电信号和肌电信号模态之间传递熵的大小来判断皮质肌肉耦合的程度。图12为整段时间内EEG→ EMG传递熵值地形图。图13为整段时间内EMG→EEG传递熵值地形图因为皮质与肌肉之间的信号存在反馈调节机制，在这种作用机制下，皮质肌肉的耦合方式存在EEG→EMG，EMG→EEG两种方向，采用这种的方法可以计算这两种方向下的传递熵值，判断不同方向上的耦合程度。采用地形图的方式来显示不同模态间传递熵数值的大小。图中，颜色越浅的传递熵数值越大。

如图14，15，16，17所示：分别按照时间尺度为1s，2s，3s， 4s连续计算实验过程中的传递熵值的变化情况，并绘制相应的传递熵变化曲线。从四个不同的时间尺度可以看出，随着时间尺度从1s 增加到4s时，传递熵变化的趋势越来越明显，与实验流程进行对照分析可以发现：在第10s开始增加负荷时，传递熵数值均开始呈现下降趋势；第20s负荷增加结束时，传递熵数值开始逐渐上升，并且上升到大于最初的传递熵水平。分析认为，负荷刚开始增加时，大脑与肌肉均未适应负荷变化打来的影响，导致皮质肌肉耦合水平下降；在负荷增加结束时，大脑和肌肉均适应了变化，而且由于负荷变大，需要调动更多的神经元和肌肉细胞来参与到运动控制中，导致皮质肌肉耦合更加强烈。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (1)

1.一种基于变分模态分解-传递熵的脑肌电耦合的研究方法，其特征在于，包括以下步骤：

步骤1)受试者按照实验指令完成实验操作，在实验过程中采集脑电信号和肌电信号，并以数字信号存储在计算机中；采集脑电信号和肌电信号具体过程为：受试者被要求端坐在指定的座椅并不做规定以外的动作，实验共进行40s：实验开始时前10s，受试者使用右手握住一个未注水的空杯，实验进行到第10s时，开始将500ml纯净水匀速的注入到受试者的空杯中，控制注水速度，使得500ml纯净水在第20s时全部注入到受试者手中的空杯内，注水结束后，再握持20s，并记录整个过程中的EEG和EMG信号；

步骤2)对上述脑电信号、肌电信号数据的预处理；

步骤3)对预处理过后的脑电信号和肌电信号进行变分模态分解，将脑电信号和肌电信号分解为若干个不同的本征模态函数；

步骤4)对使用变分模态分解得到的脑电信号、肌电信号两两之间进行传递熵计算，观察脑肌电耦合强度；

步骤5)选取不同的时间刻度再次对不同性质、不同分量、不同方向的本征模态函数进行传递熵计算，获取脑肌电耦合随负荷大小和时间变化的关系；

步骤3)的具体算法如下：

变分模态分解方法是将信号分解过程转移到变分框架内，通过搜寻约束变分模型的最优解实现信号自适应分解；首先将每个模态重新定义为调频-调幅信号，表达为：

u_k(t)＝A_k(t)cos(Ф_k(t)) (1)

式中，A_k(t)为u_k(t)的瞬时幅值，

ω_k(t)为u_k(t)的瞬时频率，假设每个模态u_k(t)具有中心频率和有限带宽，约束条件为各模态之和等于输入信号f，且每个模态的估计带宽之和最小，模型如下：

st.∑_ku_k＝f (2)

其中，u_k＝{u₁,…,u_k}，表示分解的K有限带宽IMFf分量，为求取(1)、(2)约束变分问题的最优解，引入增广形式的Largrange函数，即：

式中：α为惩罚因子，λ为拉格朗日算子；

采用乘法算子交替方向法解决以上变分问题，通过交替更新u_k ⁿ⁺¹、ω_k ⁿ⁺¹和λⁿ⁺¹，寻找上式的“鞍点”；

其中，

的表达式为：

式中，ω_k等同于

等同于

基于Parseval/Plancherel傅里叶等距变换，将(4)式转变到频域，得到各模态的频域更新；然后将中心频率的取值问题转换到频域，得到中心频率的更新方法；同时更新λ；

具体表达式如下：

每个本征模态函数分量的中心频率和带宽在迭代求解变分模型的过程中不断更新，直至满足迭代停止条件

对于给定的判定精度e＞0，结束整个循环，最终根据实际信号的频域特性得到K个窄带本征模态函数分量；

步骤(4)将经过变分模态分解得到的若干个脑电、肌电信号的本征模态函数两两之间计算传递熵值，具体算法如下：

将使用变分模态分解得到的脑电信号分量和肌电信号分量分别定义为

K₁和K₂分别代表EEG和EMG窄带分量本征模态函数的个数，基于传递熵的定义，构造

到

的变分模态分解-传递熵TE_EEG→EMG，公式如下：

式中，t离散时间指标；u为预测时间；

分别为k₁，k₂的延迟向量，p(·)为变量间的联合概率；

TE_EEG→EMG则表示EEG的k₁分量到EMG的k₂分量间的传递熵值，同理信号

到

的变分模态分解-传递熵TE_EMG→EEG的表达式为：

TE_EMG→EEG则表示EMG的k₂分量到EEG的k₁分量间的传递熵值；

步骤5)的时间刻度分别以1s、2s、3s、4s为基准，将整个时间段的本征模态函数分割再加以计算，并绘制传递熵变化曲线，观察随着时间变化和负荷逐渐增加时，皮质肌肉耦合的变化情况。

Images