CN114259242B - 基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法 - Google Patents
基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法 Download PDFInfo
- Publication number
- CN114259242B CN114259242B CN202111611263.2A CN202111611263A CN114259242B CN 114259242 B CN114259242 B CN 114259242B CN 202111611263 A CN202111611263 A CN 202111611263A CN 114259242 B CN114259242 B CN 114259242B
- Authority
- CN
- China
- Prior art keywords
- time
- scale
- frequency
- time scale
- coupling
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000010168 coupling process Methods 0.000 title claims abstract description 67
- 238000012546 transfer Methods 0.000 title claims abstract description 43
- 210000003205 muscle Anatomy 0.000 title claims abstract description 42
- 230000001054 cortical effect Effects 0.000 title claims abstract description 29
- 238000001228 spectrum Methods 0.000 title claims abstract description 27
- 230000008878 coupling Effects 0.000 claims abstract description 58
- 238000005859 coupling reaction Methods 0.000 claims abstract description 58
- 210000004556 brain Anatomy 0.000 claims abstract description 42
- 238000004458 analytical method Methods 0.000 claims abstract description 14
- 238000001914 filtration Methods 0.000 claims abstract description 9
- 238000000034 method Methods 0.000 claims description 44
- 230000009471 action Effects 0.000 claims description 32
- 239000011159 matrix material Substances 0.000 claims description 28
- 230000003595 spectral effect Effects 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 8
- 239000013598 vector Substances 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 5
- 238000012935 Averaging Methods 0.000 claims description 4
- 230000005540 biological transmission Effects 0.000 claims description 4
- 238000000354 decomposition reaction Methods 0.000 claims description 4
- 238000012880 independent component analysis Methods 0.000 claims description 4
- 230000003183 myoelectrical effect Effects 0.000 claims description 4
- 230000001419 dependent effect Effects 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 2
- 238000012360 testing method Methods 0.000 claims 1
- 210000001364 upper extremity Anatomy 0.000 abstract description 5
- 230000002146 bilateral effect Effects 0.000 abstract 1
- 238000002474 experimental method Methods 0.000 abstract 1
- 230000008569 process Effects 0.000 description 9
- 230000007423 decrease Effects 0.000 description 5
- 238000000540 analysis of variance Methods 0.000 description 3
- 210000003710 cerebral cortex Anatomy 0.000 description 3
- 230000001427 coherent effect Effects 0.000 description 3
- 230000005611 electricity Effects 0.000 description 3
- 210000003414 extremity Anatomy 0.000 description 3
- 238000000556 factor analysis Methods 0.000 description 3
- 230000007246 mechanism Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000002490 cerebral effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 239000012634 fragment Substances 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 210000000337 motor cortex Anatomy 0.000 description 2
- 238000005312 nonlinear dynamic Methods 0.000 description 2
- 208000012661 Dyskinesia Diseases 0.000 description 1
- 208000016285 Movement disease Diseases 0.000 description 1
- 206010049565 Muscle fatigue Diseases 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002457 bidirectional effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000001364 causal effect Effects 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 230000004064 dysfunction Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 230000007604 neuronal communication Effects 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 230000007170 pathology Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000001020 rhythmical effect Effects 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 230000002747 voluntary effect Effects 0.000 description 1
Landscapes
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开了一种基于多时间尺度传递谱熵的功能性皮质肌肉耦合分析方法。首先，在实验过程中采集受试者的脑电信号和上肢相关肌肉的肌电信号。然后对脑电信号和肌电信号去除伪影，去噪，滤波处理。计算处理后的脑电信号和肌电信号之间的多时间尺度传递谱熵以及尺度为s时，特定频段内的显著性面积值。然后分析双侧运动区的脑电信号与指浅屈肌的肌电信号之间的功能性皮质肌肉耦合。本发明的多时间尺度传递谱熵不仅抑制了虚假耦合，而且还描述了局部频段内的信息传递。此外，多时间尺度传递谱熵为研究功能性皮质肌肉耦合的多尺度特性提供了新的视角，而且可以更好的分析功能性皮质肌肉耦合的耦合强度、信息流、局部频带等特性。
Description
技术领域
本发明属于生物电信号处理领域，涉及一种基于脑电信号和肌电信号的功能性皮质肌肉耦合方法。
背景技术
大脑运动皮质和肌肉活动之间有节律的功能性皮质肌肉耦合是研究感觉运动系统中神经通讯的有效工具。基于脑电信号和肌电信号可以间接解释感觉运动系统中的多层次和非线性的动态特性。大脑运动皮层的脑电信号和相应的肢体对侧的肌电信号能够分别反映大脑对肢体运动的控制信息及肌肉对大脑的反馈信息，揭示了肢体运动过程中大脑皮层与相应肌肉间的功能性皮质肌肉耦合。功能性皮质肌肉耦合研究在一定程度上能够反映大脑皮层与肌肉间的相互作用关系，有助于理解运动控制与运动障碍病理机制。因此，研究探索关于功能性皮质肌肉耦合的有效方法对感觉运动系统有一定意义。
目前，许多方法被用来分析功能性皮质肌肉耦合的耦合强度、信息流、局部频带等特性。比如Grange因果关系(GC)、相干分析、互信息、传递熵(TE)。虽然相干分析可以反映脑肌电信号在频域上的线性相互作用，但是它可能不能反映皮层肌肉间的因果关系。与相干分析方法相同，Grange因果关系(GC)常用于分析线性模型，所以相干分析和GC不能有效刻画感觉运动系统中的非线性特性。互信息方法是研究脑肌电信号非线性耦合特征的一种有效方法，但是互信息无法体现功能性皮质肌肉耦合的信息流动方向且不能排除共同信号源的影响。传递熵(TE)是基于转移概率定义的非对称性测量方法。TE包含着方向性和动态信息,能够测量有效连通性且不用对复杂系统进行建模。然而，TE方法无法分析功能性皮质肌肉耦合的局部频段特性。最近，陈晓玲等人将TE方法扩展到频域，并且提出传递谱熵(TSE)方法。与TE方法相比，TSE方法不仅抑制了虚假耦合，而且还描述了局部频段内的信息传递。因此，TSE比TE更适合分析功能性皮质肌肉耦合的耦合强度、信息流、局部频带等特性。但是TSE没能进一步挖掘功能性皮质肌肉耦合的多时间尺度特性。
生理系统是跨越多时间尺度时空进行相互作用的复杂结构，多时间尺度分析可能有助于分析脑电信号和肌电信号之间不同时空层次的同步特征。Costa等人通过对时间序列进行粗粒化处理，提出多时间尺度熵(MSE)方法用于分析生理信号中的多时间尺度和动态特性。随后，谢平等人将粗粒化处理引入TE，提出多时间尺度传递熵(MSTE)用于定量研究脑电信号和肌电信号在不同时间尺度上的非线性耦合特征。然而，粗粒化处理存在一些问题，对给定时间序列的粗粒化过程，会随着其尺度的增加，大幅度减少这一时间序列的长度。并且当原始数据长度不是尺度因子的整数倍时，将会造成部分数据的丢失。因此，当粗粒化过程应用到短时间序列时，序列长度的减小会导致TE估计的精确性降低。为了解决这些问题，wu等人提出了复杂多时间尺度化(MA)方法，并且将传统MSE算法中的粗粒化过程替换为复杂多时间尺度化方法(MA)，进一步得到改进的修正多时间尺度熵(MMSE)。与MSE相比，MMSE可以更精确地估计熵，避免引入不确定熵。
因为以上方法的局限性，提出一种合适的方法来分析功能性皮质肌肉耦合的多时间尺度特性以及耦合强度、信息流、局部频带等特性具有一定的意义。
发明内容
考虑到生理信号中固有的多时间尺度特征，为了克服现有方法的局限性，本发明提供了一种基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法。
本发明首先，采集受试者在三级稳态握力任务中的脑电信号和上肢相关肌肉的肌电信号。然后应用多时间尺度传递谱熵计算得到多时间尺度和多频段显著性面积值矩阵，根据目标函数G选择出的26的时间尺度用于接下来的研究。最后，在26的时间尺度下，通过多因素方差方法(ANOVA)得到受试者在不同频段条件下，不同实验动作条件下或者不同方向条件(脑电信号到肌电信号或肌电信号到脑电信号)下的显著性面积值差异。
为了实现以上目的，本发明方法主要包括以下步骤：
步骤1：采集受试者的脑电信号和肌电信号。
(1)以1000Hz的采样频率分别采集受试者在做实验动作时的脑电信号；
(2)以1000Hz的采样频率分别采集受试者在做实验动作时对应肌肉位置的肌电信号；
步骤2：使用EEGLAB工具箱对脑电信号进行预处理，独立分量分析用于去除脑电信号的伪影。然后采用经验模态分解(EMD)和小波阈值相结合的方法对肌电信号进行去噪处理。对脑电信号和肌电信号进行1-60HZ的带通滤波，最后使用固定滞后卡尔曼平滑器滤除脑电信号和肌电信号的50HZ电源频率。
步骤2：处理脑电信号和肌电信号。
2-1.独立分量分析去除脑电信号的伪影；
2-2.通过经验模态分解(EMD)和小波阈值相结合对肌电信号进行去噪处理；
2-3.对脑电信号和肌电信号进行1-60HZ的带通滤波。
2-4.使用固定滞后卡尔曼平滑器滤除脑电信号和肌电信号的50HZ电源频率。
步骤3：通过多时间尺度传递谱熵和显著性面积值计算方法分别构建不同实验动作以及不同方向(脑电到肌电或肌电到脑电)的多时间尺度和多频段显著性面积值矩阵，通过这些矩阵得到受试者功能性皮质肌肉耦合的多时间尺度特性、多频段特性、耦合强度特性和信息流特性。
3-1定义多时间尺度传递谱熵。具体步骤如下：
(1)设时间尺度为s，s取整数值且不超过时间序列长度。
(2)对于具有耦合关系的时间序列，比如脑电和肌电时间序列。X＝{x(t),t＝1,2,…,T}和Y＝{y(t),t＝1,2,…,T}，其中，序列长度均为T(T≥1)。利用复杂多时间尺度化对信号进行多时间尺度化，X和Y经过复杂多时间尺度化后，得到与尺度相关的时间序列分别为XS和YS，公式如下：
式中，s为尺度，j为信号XS和YS的序号，T为序列X和Y的长度。
(3)对于尺度为s的时间序列XS和YS，基于相空间重构方法分别将原始时间序列XS和YS拓展到M维的向量空间中，得到矩阵Ps和Qs：
其中τ为延迟时间,可以由相关积分法确定，m为嵌入维数，可由Cao准则确定，s为当前时间尺度。
(4)然后对Ps和Qs的矩阵进行二维傅里叶变换得到时间尺度为s的矩阵VS和WS：
其中频率分辨率Δf＝fs/N,fs为采样频率，N＝T-s+1是数据的长度。
(5)对于给定的一个频率值f，计算矩阵VS和WS中对应频率值f的列向量vs(f)和ws(f)之间的多时间尺度传递谱熵：
其中f是离散频率值，τ是频率预测值，和分别是向量vs(f)和ws(f)在频域点f处的φ-维和φ-维的延时向量。/>表示频率f处的k+δ时刻序列vs(f)基于k时刻序列vs(f)下的条件熵，/>则表示频率f值处的k+δ时刻序列vs(f)基于k时刻序列vs(f)和ws(f)下的条件熵。同理可以得到从YS到XS的/>
3-2.定义特定频段内的显著性面积值AS。
计算时间尺度为s时，特定频段内的显著性面积值AS：
式中Δf表示频率分辨率，为当前频段的最低频率，/>为当前频段的最高频率。AS值越大，则两个信号之间的耦合强度越强。即AS能够定量描述脑电信号与肌电信号间在多时间尺度上的非线性同步耦合特征，s为当前时间尺度。
3-3.选择需要构建矩阵的实验动作及其信息流传输方向(脑电信号到肌电信号或肌电信号到脑电信号)。
3-4.选择多个时间尺度，时间尺度的可选数值范围为1到脑电及肌电信号时间序列长度。
3-5.分别计算各个时间尺度在多个频段下的多时间尺度传递谱熵。
(1)计算相关脑电信号和肌电信号之间在范围为1～60HZ中的每一个整数频率的当前方向的多时间尺度传递谱熵。
(2)针对得到的每一个频率的多时间尺度传递谱熵，计算不同频段下(α频段：8～14Hz，β1频段：15～25Hz，β2频段：25～35Hz，γ1频段：35～45Hz，γ2频段：45～60Hz)的显著性面积值。对相同频段上得到的显著性面积值进行平均。
(3)对与每个受试者在相同频段上的显著性面积值取平均，得到时间尺度为s时的多个频段下的多时间尺度传递谱熵。
(4)替换时间尺度s为下一个需要计算的时间尺度，重复(1)-(3)步骤。
(5)将计算出来的值按矩阵行设为不同时间尺度，矩阵列设为不同频段，值越高的填充越深的红色，值越低的地方填充越深的蓝色的规则绘图，得到多时间尺度和多频段显著性面积值矩阵。
3-6.选择其他需要构建矩阵的实验动作及其信息流传输方向重复以上步骤。
3-7.通过从不同角度观察以上得到的多个矩阵得到受试者功能性皮质肌肉耦合的多时间尺度特性、多频段特性、耦合强度特性和信息流特性。
(1)沿着每个矩阵横向观察，探寻随着时间尺度增大的耦合强度变化规律。
(2)沿着每个矩阵纵向观察，探寻哪一个频段上的耦合强度相对于其他频段较大。
(3)通过观察矩阵各个区域颜色的深浅可以得知某一个区域耦合强度的大小。
(4)观察同一个实验动作的不同方向(脑电信号到肌电信号或肌电信号到脑电信号)的矩阵，可以得知某一个区域上的不同方向的耦合强度强弱。
步骤4：定义目标函数G，根据目标函数G选择最优时间尺度。从显著性面积值矩阵中选择最优时间尺度以及显著性面积值最大的频段，用多因素方差分析(ANOVA)方法计算这些部分值的显著差异得出受试者左右手之间、不同信息流方向之间、不同实验动作之间的耦合差异。
4-1.定义目标函数G：
其中，M和S分别表示显著性面积值AS的平均值和标准差.最优时间尺度参数确定为G达到峰值时的尺度。
4-2.从不同显著性面积值矩阵中选择最合适的相同实验动作的多个频段在所有尺度下的显著性面积值以及相同频段下的多个实验动作在所有尺度下的显著性面积值。计算他们的目标函数G的值。
4-3.观察目标函数G的值随着时间尺度增长的变化情况，选择最多次使得G达到峰值的尺度值作为最优时间尺度。下面的步骤基于最优时间尺度。
4-4.选取显著性面积值矩阵最优时间尺度上的值，从中选取红色最深的频段，即显著性面积值最强，也即耦合强度最强的频段进行下面的步骤。
4-5.以使用手(左手和右手)、不同实验动作、信息流方向(脑电到肌电、肌电到脑电)显著性面积为自变量，显著性面积值为因变量进行三因素方差分析。
4-6.随后对左手和右手，所有实验动作，方向(脑电到肌电、肌电到脑电)进行多重比较分析。得到相同实验动作及相同信息流方向上的左右手之间、相同手及相同实验动作的不同信息流方向之间、相同手及相同信息流方向上的不同实验动作之间的显著性耦合差异。
本发明具有如下有益效果：
以往的对于时间序列多时间尺度的处理存在一些问题，对给定时间序列的多时间尺度化过程，会随着其尺度的增加，大幅度减少这一时间序列的长度。并且当原始数据长度不是尺度因子的整数倍时，将会造成部分数据的丢失。因此，当多时间尺度化过程应用到短时间序列时，序列长度的减小会导致俩个信号之间耦合强度估计的精确性降低。而本发明提出的多时间尺度传递谱熵方法把传统的多时间尺度化过程替换为复杂多时间尺度化方法。与传统的多时间尺度化方法相比，复杂多时间尺度化方法(MA)保证了各尺度上多时间尺度化后的时间序列与原始时间序列长度一致，并且避免了数据丢失的现象，大大提高了方法的准确性。多时间尺度传递谱熵分析功能性皮质肌肉耦合的多时间尺度特性提供了新的视角，为进一步探索运动功能障碍的产生机理及康复过程中的运动功能评价方法提供了依据。
附图说明
图1为本发明的实施流程图；
图2为右手执行握力任务时，相关的单个通道脑电信号和上肢相关肌肉的肌电信号的多时间尺度和多频段显著性面积值矩阵。
(a)表示15％MVC握力，脑电信号→肌电信号；
(b)表示25％MVC握力，脑电信号→肌电信号；
(c)表示50％MVC握力，脑电信号→肌电信号；
(d)表示15％MVC握力，肌电信号→脑电信号；
(e)表示25％MVC握力，肌电信号→脑电信号；
(f)表示50％MVC握力，肌电信号→脑电信号；
图3为左手执行握力任务时，相关的单个通道脑电信号和上肢相关肌肉的肌电信号的多时间尺度和多频段显著性面积值矩阵。
(a)表示15％MVC握力，脑电信号→肌电信号；
(b)表示25％MVC握力，脑电信号→肌电信号；
(c)表示50％MVC握力，脑电信号→肌电信号；
(d)表示15％MVC握力，肌电信号→脑电信号；
(e)表示25％MVC握力，肌电信号→脑电信号；
(f)表示50％MVC握力，肌电信号→脑电信号；
图4为利用目标函数G来确定最优的时间尺度参数。
图5用不同握力握拳条件下左右手之间显著性面积值的比较以及脑电信号和肌电信号不同信息流方向之间显著性面积值的比较。
图6脑电信号和肌电信号不同信息传递方向条件下不同握力之间显著性面积值的比较。
具体实施方式
本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程。下面结合附图对本发明的实施例作详细说明：
如图1所示，本实施例包括如下步骤：
步骤一，采集受试者的脑电信号和肌电信号。
对于6名健康受试者(四名男性和两名女性)，首先使用弹簧式握力计(中国，LynxMall，EH101)测量每个受试者的最大自愿收缩(MVC)。之后，受试者执行15％MVC、25％MVC、50％MVC握力水平的握力任务。所有受试者分别用左手和右手执行每个握力任务。在同一握力水平的所有任务完成后，受试者休息20min，然后切换到下一个握力水平的握力任务，以防止肌肉疲劳。在连接电极之前，每个受试者都清洁自己的头发和上肢皮肤。
然后以1000Hz的采样频率采集受试者在做三级稳态握力任务时的脑电信号，采集的脑电信号包括C3、C4通道；
以1000Hz的采样频率分别采集受试者在做实验动作时对应肌肉位置的肌电信号，采集的肌电信号包括指浅屈肌(FDS)的信号。
步骤二：处理脑电信号和肌电信号。
首先用独立分量分析去除脑电信号的伪影，如果握力偏离5％MVC的范围，或在1-60Hz范围内出现较大幅度伪影，则丢弃对应脑电信号片段，并排除这些脑电信号片段对应的肌电信号。
采用经验模态分解(EMD)和小波阈值相结合的方法对肌电信号进行去噪处理。
对脑电信号和肌电信号进行1-60Hz范围的带通滤波。
最后使用固定滞后卡尔曼平滑器滤除脑电信号和肌电信号的50Hz电源频率。
步骤三：通过多时间尺度传递谱熵和显著性面积值计算方法分别构建不同实验动作以及不同方向(脑电到肌电或肌电到脑电)的多时间尺度和多频段显著性面积值矩阵，通过这些矩阵得到受试者功能性皮质肌肉耦合的多时间尺度特性、多频段特性、耦合强度特性和信息流特性。
选择每一个受试者的每一不同握力级别的握力任务和每一种信息流传输方向，在尺度s为2,6,10,14,18,22,26,30,34,38,42时，分别进行以下计算：
计算C4通道脑电信号和左手指浅屈肌的肌电信号之间在频率范围为1～60HZ(步长为1HZ，即：1,2,3...60HZ)的双向(脑电信号→肌电信号和肌电信号→脑电信号)多时间尺度传递谱熵。类似地，可以计算出C3通道脑电信号和右手指浅屈肌的肌电信号之间在频率范围为1～60HZ的双向多时间尺度传递谱熵。
然后针对得到的每一组多时间尺度传递谱熵，根据显著性面积公式计算不同频段下(α频段：8～14Hz，β1频段：15～25Hz，β2频段：25～35Hz，γ1频段：35～45Hz，γ2频段：45～60Hz)的显著性面积值。
对5次重复握力任务计算得到的5个显著性面积值进行平均，得到每个受试者在不同握力任务中的不同频段上的显著性面积值。最后对每个组的6名受试者在特定任务的特定频段上的显著性面积值取平均。每个受试者组的显著性面积值数量如下：健康对照组：6个握力任务(左手和右手的15％MVC、25％MVC和50％MVC)×5个频段×2个方向,患者组：4个握力任务(左手和右手的15％MVC、25％MVC)×5个频段×2个方向。
通过以上计算出来的显著性面积值得出多时间尺度和多频段显著性面积值矩阵，如图2和图3所示。
结果表明，大脑运动皮质和指浅屈肌(FDS)之间在不同尺度上存在双向耦合关系，并且在不同方向上(EEG→EMG及EMG→EEG),耦合强度在β1、β2、γ2频段上更高。β1频段的耦合强度在绝大部分尺度上都比其他各个频段的耦合强度要强。随着尺度的增加，各个频段上的耦合强度逐渐递增，在26～30尺度耦合强度达到最大后逐渐降低。
步骤4：定义目标函数G，根据目标函数G选择最优时间尺度。从显著性面积值矩阵中选择最优时间尺度以及显著性面积值最大的频段，用多因素方差分析(ANOVA)方法计算这些部分的值的显著差异得出受试者左右手之间、不同信息流方向之间、不同实验动作之间的耦合差异。
定义目标函数G用于选择合适的尺度参数。
其中，M和S分别表示显著性面积值的平均值和标准差，最优尺度参数确定为G达到峰值时的尺度。
最优时间尺度的具体选择流程如下：在25％MVC握力任务中，步长为4的范围为2～42的每个尺度中：通计算出β1频段上左手和右手在俩个方向上的目标函数G的值；同理通过计算出在右手EEG→EMG上方向上不同频段的目标函数G的值，得到图4。如图4所示，在四种握力任务中目标函数G的最大值所对应的尺度分别为22、26、30、26。在5个频段中目标函数G的最大值所对应的尺度分别为22、26、26、30、34。综上所述，选择26作为最优尺度进行接下来的分析研究。
因为β1频段在所有频段中耦合强度最大，在β1频段上，进一步分析不同握力之间，不同使用手之间和不同方向之间的显著性面积值的差异。在β1频段上，本实施例进一步研究了不同握力之间，不同使用手之间和不同方向之间的差异。以使用手(左手和右手)、握力大小(15％MVC、25％MVC、50％MVC)、方向(脑电信号→肌电信号，肌电信号→脑电信号)为自变量，显著性面积值为因变量进行三因素方差分析。使用手(F(1,60)＝12.256,p<0.001)、握力大小(F(1,60)＝8.4274,p＝0.005)和方向(F(1,60)＝5.124,p＝0.027)均存在主效应。
随后对左手和右手，三种握力大小，脑电信号→肌电信号，肌电信号→脑电信号进行多重比较分析。
从图5中可以看出，在所有的三种握力中，无论是在脑电信号→肌电信号方向还是肌电信号→脑电信号方向上，右手的显著性面积值都要高于左手(p<0.05)。在所有的握力任务中，右手15％MVC、25％MVC和左手50％MVC握力任务脑电信号→肌电信号方向上的显著性面积值高于肌电信号→脑电信号方向(p<0.05)。从图6中可以看出在所有的四个通道中，50％MVC握力的AS值都显著高于25％MVC握力和15％MVC握力的As值(p<0.5)。然而，25％MVC握力下的显著性面积值仅仅在C3->FDS通道高于15％MVC握力的显著性面积值(F(2,60)＝3.568,p＝0.034)。
综上，相对于其他功能性皮质肌肉耦合分析方法，多时间尺度传递谱熵不仅抑制了虚假耦合，而且还描述了局部频段内的信息传递。而且多时间尺度传递谱熵可以更好的得出脑电和肌电之间的多层次和非线性的动态特性，为研究功能性皮质肌肉耦合的多时间尺度特性提供了新的视角。使用此方法可以更好的有助于理解运动控制与运动障碍病理机制。
Claims (6)
1.基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法，其特征在于该方法包括如下步骤：
步骤1：采集受试者的脑电信号和肌电信号；
步骤2：处理脑电信号和肌电信号；
步骤3：通过多时间尺度传递谱熵和显著性面积值计算方法分别构建不同实验动作以及不同方向的多时间尺度和多频段显著性面积值矩阵，通过这些矩阵得到受试者功能性皮质肌肉耦合的多时间尺度特性、多频段特性、耦合强度特性和信息流特性；
步骤4：定义目标函数G，根据目标函数G选择最优时间尺度；从显著性面积值矩阵中选择最优时间尺度以及显著性面积值最大的频段，用多因素方差分析方法计算这些部分的值的显著差异得出受试者左右手之间、不同信息流方向之间、不同实验动作之间的耦合差异；
其中步骤1包括：
以1000Hz的采样频率分别采集受试者在做实验动作时的脑电信号；
以1000Hz的采样频率分别采集受试者在做实验动作时对应肌肉位置的肌电信号。
2.根据权利要求1所述的基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法，其特征在于：所述的步骤2具体包括：
2-1、独立分量分析去除脑电信号的伪影；
2-2、通过经验模态分解和小波阈值相结合对肌电信号进行去噪处理；
2-3、对脑电信号和肌电信号进行带通滤波；
2-4、使用固定滞后卡尔曼平滑器滤除脑电信号和肌电信号的电源频率。
3.根据权利要求1所述的基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法，其特征在于：所述的步骤3中的多时间尺度传递谱熵具体设计如下：
3-1.设时间尺度为s，s取整数值且不超过时间序列长度；
3-2.使用复杂多时间尺度化方法对脑电信号时间序列和肌电信号时间序列进行多时间尺度化；
3-3.基于相空间重构方法分别将步骤3-2中复杂多时间尺度化后的脑电和肌电时间序列拓展到M维的向量空间中，从而使得原来的一维时间序列都变为矩阵；
3-4.对步骤3-3的两个矩阵都进行二维傅里叶变换，得到这两个矩阵的频域表示；
3-5.基于步骤3-4得到的矩阵对于给定的一个频率值f，计算矩阵中对应频率值的两个列向量之间的多时间尺度传递谱熵。
4.根据权利要求1所述的多时间尺度传递谱熵的功能性皮质肌肉耦合方法，其特征在于：步骤3中所述的显著性面积值计算方法设计如下
其中频率分辨率Δf＝fs/N，fs为采样频率，N是数据的长度，为当前频段的最低频率，为当前频段的最高频率，s为当前时间尺度。
5.根据权利要求1所述的多时间尺度传递谱熵的功能性皮质肌肉耦合方法，其特征在于：步骤3中具体步骤如下：
5-1.选择需要构建矩阵的实验动作及其信息流传输方向；
5-2.选择多个时间尺度，时间尺度的可选数值范围为1到脑电及肌电信号时间序列长度；
5-3.分别计算各个时间尺度在多个频段下的多时间尺度传递谱熵；
(1)计算相关脑电信号和肌电信号之间在范围为1～60HZ中的每一个整数频率的当前方向的多时间尺度传递谱熵；
(2)针对得到的每一个频率的多时间尺度传递谱熵，计算不同频段下的显著性面积值；对相同频段上得到的显著性面积值进行平均；
(3)对与每个受试者在相同频段上的显著性面积值取平均，得到时间尺度为s时的多个频段下的多时间尺度传递谱熵；
(4)替换时间尺度s为下一个需要计算的时间尺度，重复(1)-(3)步骤；
(5)将计算出来的值按矩阵行设为不同时间尺度，矩阵列设为不同频段，值越高的填充越深的红色，值越低的地方填充越深的蓝色的规则绘图，得到多时间尺度和多频段显著性面积值矩阵；
5-4.选择其他需要构建矩阵的实验动作及其信息流传输方向重复以上步骤；
5-5.通过从不同角度观察以上得到的多个矩阵得到受试者功能性皮质肌肉耦合的多时间尺度特性、多频段特性、耦合强度特性和信息流特性。
6.根据权利要求1所述的多时间尺度传递谱熵的功能性皮质肌肉耦合方法，其特征在于：步骤4中具体步骤如下：
6-1.定义目标函数G的计算公式；
6-2.从不同显著性面积值矩阵中选择最合适的相同实验动作的多个频段在所有尺度下的显著性面积值以及相同频段下的多个实验动作在所有尺度下的显著性面积值；计算他们的目标函数G的值；
6-3.观察目标函数G的值随着时间尺度增长的变化情况，选择最多次使得G达到峰值的尺度值作为最优时间尺度；下面的步骤基于最优时间尺度；
6-4.观察显著性面积值矩阵最优时间尺度上的值，从中选取红色最深的频段，即耦合强度最强的频段进行下面的步骤；
6-5.以使用手、不同实验动作、信息流方向的显著性面积为自变量，显著性面积值为因变量进行三因素方差分析；
6-6.随后对左手和右手，所有实验动作，信息流方向进行多重比较分析；得到相同实验动作及相同信息流方向上的左右手之间、相同手及相同实验动作的不同信息流方向之间、相同手及相同信息流方向上的不同实验动作之间的耦合差异。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111611263.2A CN114259242B (zh) | 2021-12-27 | 2021-12-27 | 基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111611263.2A CN114259242B (zh) | 2021-12-27 | 2021-12-27 | 基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114259242A CN114259242A (zh) | 2022-04-01 |
| CN114259242B true CN114259242B (zh) | 2024-04-12 |
Family
ID=80830308
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111611263.2A Active CN114259242B (zh) | 2021-12-27 | 2021-12-27 | 基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114259242B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106821318A (zh) * | 2016-08-31 | 2017-06-13 | 天津市人民医院 | 一种脑电信号多尺度定量分析方法 |
| CN109144259A (zh) * | 2018-08-23 | 2019-01-04 | 杭州电子科技大学 | 一种基于多尺度排列传递熵的脑区间同步关系分析方法 |
| CN111227830A (zh) * | 2020-02-14 | 2020-06-05 | 燕山大学 | 一种基于复杂改进多尺度传递熵的脑肌电耦合分析方法 |
| CN113499084A (zh) * | 2021-07-09 | 2021-10-15 | 燕山大学 | 脑卒中运动功能评价模型建立方法、评价方法及系统 |
-
2021
- 2021-12-27 CN CN202111611263.2A patent/CN114259242B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106821318A (zh) * | 2016-08-31 | 2017-06-13 | 天津市人民医院 | 一种脑电信号多尺度定量分析方法 |
| CN109144259A (zh) * | 2018-08-23 | 2019-01-04 | 杭州电子科技大学 | 一种基于多尺度排列传递熵的脑区间同步关系分析方法 |
| CN111227830A (zh) * | 2020-02-14 | 2020-06-05 | 燕山大学 | 一种基于复杂改进多尺度传递熵的脑肌电耦合分析方法 |
| CN113499084A (zh) * | 2021-07-09 | 2021-10-15 | 燕山大学 | 脑卒中运动功能评价模型建立方法、评价方法及系统 |
Non-Patent Citations (2)
| Title |
|---|
| Multiscale Transfer Spectral Entropy for Quantifying Corticomuscular Interaction;Jinbiao Liu 等;《IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS》;20210630;第25卷(第6期);第2281-2292页 * |
| Transfer Spectral Entropy and Application to Functional Corticomuscular Coupling;Xiaoling Chen 等;《IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING》;20190531;第27卷(第5期);第1092-1102页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114259242A (zh) | 2022-04-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Fang et al. | Extracting features from phase space of EEG signals in brain–computer interfaces | |
| Phinyomark et al. | Application of wavelet analysis in EMG feature extraction for pattern classification | |
| Jha et al. | Electrocardiogram data compression using DCT based discrete orthogonal Stockwell transform | |
| CN110367974B (zh) | 基于变分模态分解-传递熵的脑肌电耦合的研究方法 | |
| CN111067514B (zh) | 一种基于多尺度多变量传递熵的多通道脑电耦合分析方法 | |
| Chen et al. | Transfer spectral entropy and application to functional corticomuscular coupling | |
| CN111227830B (zh) | 一种基于复杂改进多尺度传递熵的脑肌电耦合分析方法 | |
| CN113499084B (zh) | 脑卒中运动功能评价模型建立方法、评价方法及系统 | |
| AlOmari et al. | Analysis of extracted forearm sEMG signal using LDA, QDA, K-NN classification algorithms | |
| JP5252432B2 (ja) | 指関節角度推定装置 | |
| Cranstoun et al. | Time-frequency spectral estimation of multichannel EEG using the auto-SLEX method | |
| Azami et al. | Entropy analysis of univariate biomedical signals: Review and comparison of methods | |
| CN109674445B (zh) | 一种结合非负矩阵分解和复杂网络的肌间耦合分析方法 | |
| Mir et al. | ECG denoising and feature extraction techniques–a review | |
| CN111543946A (zh) | 基于改进变分模态分解算法的癫痫脑电信号自动检测方法 | |
| CN105266804A (zh) | 一种基于低秩和稀疏矩阵分解的新型脑电信号处理方法 | |
| CN107320097B (zh) | 利用肌电信号边际谱熵提取肌肉疲劳特征的方法和装置 | |
| Van Putten et al. | Application of a neural complexity measure to multichannel EEG | |
| Rosso et al. | Analysis of wavelet-filtered tonic-clonic electroencephalogram recordings | |
| Xi et al. | Analysis of functional corticomuscular coupling based on multiscale transfer spectral entropy | |
| CN114387668A (zh) | 一种基于多层次神经肌肉耦合特征信息融合的分类方法 | |
| CN114259242B (zh) | 基于多时间尺度传递谱熵的功能性皮质肌肉耦合方法 | |
| Gupta et al. | A novel approach for extracting feature from eeg signal for mental task classification | |
| CN115299960A (zh) | 基于短时变分模态分解的电信号分解方法及脑电信号分解装置 | |
| Shaw et al. | Efficacy of adaptive directed transfer function for neural connectivity estimation of EEG signal during meditation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |