CN111708978B - 多尺度时频肌间耦合分析方法 - Google Patents

Abstract

本发明公开了一种多尺度时频肌间耦合分析方法。本发明首先进行多通道表面肌电信号同步采集并进行预处理；并对预处理后的数据进行噪声辅助的多元经验模态分解，得到有用的IMF尺度分量。其次对IMF尺度分量进行同步提取变换；具体为：对每个IMF尺度分量进行短时傅里叶变换，再乘以一个相位因子后进行同步压缩变换。然后计算时频互信息、时频归一化互信息和时频条件互信息；最后将上述计算结果进行多尺度时频肌间耦合统计分析。本发明为定量研究脑卒中患者上肢康复运动过程中不同时频尺度下的肌间非线性耦合强度特性提供了一种新方法。

Description

多尺度时频肌间耦合分析方法

技术领域

本发明属于神经系统运动控制机制研究领域，涉及噪声辅助的多元经验模态分解、时频互信息、时频归一化互信息和时频条件互信息的计算，从而进行多尺度时频肌间耦合特性分析。

背景技术

如何促进脑卒中患者上肢运动功能康复是目前来康复医学界关注的热点和难点。康复治疗过程中运动功能评估具有重要作用，一方面评估结果能够实时地反映卒中患者运动功能障碍程度，另一方面还可以辅助临床诊断并指导个性化康复训练方案的制定。

从人体骨骼肌表面通过电极引导、记录下来的表面肌电信号(surfaceElectromyography,sEMG)与肌肉的活动状态和功能状态存在着不同程度的关联，能够在一定程度上反映神经、肌肉的活动，经常被用于定性评价与定量分析脑卒中患者肢体肌肉运动功能的恢复情况。

由于运动神经系统是通过神经振荡的方式来传递运动控制信息，与肌间各区域的相互作用来实现功能调节和整合，体现为不同层次的同步耦合，而脑卒中患者正是由于运动中枢神经受损阻碍了神经振荡的传导才导致了运动功能障碍，所以通过肌间耦合分析能够反映运动控制过程中神经肌肉的异常功能联系，揭示神经网络的协同工作方式，理解运动功能障碍的病理机制。

近年来，肌间耦合分析算法层出不穷，相干性(Coherence)因其算法原理简单、易于实现被广泛应用到肌间耦合分析中。Xie等利用Coherence对比分析了中风患者运动过程中健侧、患侧上肢拮抗肌间一致性特征，发现患侧运动过程中beta频段肌间一致性相对于健侧存在明显缺少。然而，Coherence仅能描述线性耦合关系，Faes等指出Coherence所测得的耦合关系包含了直接和间接的影响，会过度估计肌间耦合强度。

Summers等研究表明，肌间耦合存在线性与非线性的关系，且非线性耦合在神经运动控制中起着更加重要的作用。为此，Du等构建了交叉频率相干-非负矩阵分解(CFC-NMF)模型，研究上肢肘关节屈伸运动中肌间线性和非线性耦合特性，发现不同频率比下中风患者肌间耦合强度明显低于健康人。然而，生理系统需要跨越多个时空尺度进行控制，导致sEMG信号具有多尺度特性，传统的单尺度分析只能从单一尺度上反映信号的特性，不能反映尺度变化导致的信号特性的变化。为此，Ma等将多元经验模态分解(MultivariateEmpirical Mode Decomposition,MEMD)与传递熵(Transfer Entropy,TE)结合，应用于健康人和病人静态握力输出条件下脑肌电同步耦合分析，研究发现皮层肌肉功能耦合是双向的，在高频段(40～75Hz)脑电对肌电的耦合强度大于肌电对脑电的耦合强度。但MEMD易受噪声影响，且存在模态混叠等问题。并且，sEMG信号是一种非线性、非平稳的混沌信号，单从时域或频域角度剖析肌间耦合关系是比较片面的，为此，有学者根据Coherence原理结合时频分析方法，提出小波相干性(Wavelet Coherence,WC)应用于电生理信号的耦合分析，但其仍遗留了线性时频变换的缺点，且无法给出定量分析指标。

最近，Rehman等人结合EEMD和MEMD的优点，提出了噪声辅助的多元经验模式分解(Noise Assisted Multivariate Empirical Mode Decomposition,NA-MEMD)方法，不但充分利用了MEMD处理白噪声时具有的二元滤波器组特性，而且加入额外的独立白噪声保证了MEMD分解后信号的本征模态函数(Intrinsic Mode Function，IMF)分量与噪声的IMF分量之间的可分离性。对分解得到的多个平稳的IMF尺度分量做时频分析，既考虑了信号的多尺度特性，又能兼顾时频局部特性。

受限于Heisenberg测不准原理的传统时频分析方法，时频分辨率较低，Daubechies等于2011年提出了一种新的时频变换方法，即同步压缩变换(Synchrosqueezing Transform,SST)，该方法在小波尺度方向上利用一种重分配方法将时间-尺度平面转化为时间-频率平面，在提高时频聚集性的同时支持信号重构，但SST无法处理多分量、强调频信号，具有较差的噪声鲁棒性。最近，Yu等提出一种新的时频分析方法—同步提取变换(Synchroextracting Transform,SET)。SET也是一种时频分析后处理技术，与SST不同的是，SET通过仅保留在时频脊线位置处的时频系数，剔除其余发散能量，以此达到时频聚焦、提高时频分辨率的目的。考虑到现有的同步耦合分析方法无法同时从时频联合域的角度出发，定量地反映时频耦合特性。基于信息论的互信息(Mutual Information,MI)、归一化互信息(Normalized Mutual Information,NMI)和条件互信息(ConditionalMutual Information,CMI)方法不依赖既定模型，能够度量两个或多个随机变量之间的非线性依赖程度，被广泛应用在电生理信号的同步耦合分析之中。直观地，MI衡量了变量间的共享信息，即当知道其中一个变量时减少了多少另一个变量的不确定性。NMI是MI的归一化版本，消除了量纲影响。当已知第三方变量的情况下，变量间直接的MI即为CMI。

为定量研究脑卒中患者上肢康复运动过程中不同时频尺度下的肌间非线性耦合强度特性，本发明提出一种新的时频耦合(Time Frequency Coupling,TFC)分析方法，并将TFC与NA-MEMD相结合，构建多尺度时频耦合(Multi-scale TFC,Ms-TFC)分析模型，来进行肌间耦合分析。

发明内容

本发明的目的在于提供一种可得到上肢肌间不同时频尺度下的功能耦合特性的分析方法。

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

步骤(1)，多通道表面肌电信号同步采集并进行预处理；

具体为：同时采集5名健康被试和5名卒中患者执行及物运动任务时上肢3通道表面肌电信号；对获取的数据进行去均值，去基线漂移，同时利用谱插值算法抑制50Hz工频干扰。

步骤(2)，对预处理后的数据进行噪声辅助的多元经验模态分解；

具体为：向预处理后的m通道表面肌电信号X＝{x₁,x₂,…,x_i,…x_m}中添加两路互不相关的高斯白噪声，组成多元复合信号X_c＝{x₁,x₂,…,x_m,x_m+1,x_m+2}，长度为T，经过多次多元经验模态分解(MEMD)，并将与噪声相关的两路残余分量剔除，得到m×q个有用的IMF尺度分量。

步骤(3)，对IMF尺度分量进行同步提取变换。

具体为：首先对每个IMF尺度分量imf_i ^k(t)进行短时傅里叶变换(STFT)，然后再乘以一个相位因子e^jωt，可得到

其中，k∈{1,2,...,m}，i∈{1,2,...,q}。

根据同步压缩变换(SST)的思想，同步提取变换(SET)可表示为

SET(t,ω)＝STFT_e(t,ω)·δ(ω-ω₀(t,ω))   (2)

其中，δ(ω-ω₀(t,ω))称为同步提取算子(SEO)，δ(*)为delta函数，ω₀(t,ω)为STFT谱中时频系数的瞬时频率。

步骤(4)，计算时频互信息、时频归一化互信息和时频条件互信息。

具体为：在同一IMF尺度分量上，计算第k个通道和第l个通道的sEMG信号的时频互信息(TFMI)、时频归一化互信息(TFNMI)和时频条件互信息(TFCMI)，即

TFMI(U,V)＝H(U)+H(V)-H(U,V)

TFCMI(U,V|W)＝H(U,W)+H(U,W)-H(W)-H(U,V,W)   (3)

其中，U、V、W分别为

其中，k≠l≠s∈{1,2,...,m},i∈{1,2,...,q}，H(*)表示时频熵，H(*,*)和H(*,*,*)表示时频联合熵，对数运算时基底取为2，单位为比特(Bit)。

步骤(5)，采用步骤(4)的计算结果进行多尺度时频肌间耦合统计分析。

本发明与传统的肌间耦合分析方法相比，具有如下有益效果：

传统的肌间耦合分析方法并不能有效地描述sEMG信号在不同时频尺度上线性与非线性、直接与间接的功能耦合特性，针对该问题，本发明首先提出一种新的时频耦合(Time Frequency Coupling,TFC)分析方法。该方法计算SET变换后时频能量谱间的MI、NMI和CMI，记为时频互信息(TFMI)、时频归一化互信息(TFNMI)和时频条件互信息(TFCMI)，利用TFMI、TFNMI和TFCMI指标度量时频能量谱间的耦合强度。然后将TFC与NA-MEMD相结合，构建多尺度时频耦合(Multi-scale TFC,Ms-TFC)分析模型，为定量研究脑卒中患者上肢康复运动过程中不同时频尺度下的肌间非线性耦合强度特性提供新方法。

附图说明

图1为本发明方法的流程图。

图2(a)和图2(b)为分别为受试者H1和S1预处理后的3通道sEMG信号。

图3(a)和图3(b)分别为受试者H1和S1的3通道sEMG信号经NA-MEMD分解后的IMF分量和残余分量。

图4(a)和图4(b)分别为受试者H1和S1前4个IMF分量的SET平均时频表示。

图5(a)和图5(b)分别为受试者H1和S1在4个IMF尺度上肌间的时频耦合强度值。

具体实施方式

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

有研究表明，脑卒中患者正是由于上运动神经元受损，运动系统失去高位中枢的控制，导致运动控制和感觉神经通路部分或完全阻断，从而引起异常运动模式。肌间耦合，即在运动中枢神经的控制下不同肌肉间的相互关联和相互作用，能够反映大脑皮层与肌肉间不同层次的功能联系。因而，通过对肌间耦合特性的研究可以探索脑卒中患者运动功能障碍的产生机制，为康复训练进展评估提供理论依据，具有十分重要的意义。

本发明提出一种新的肌间时频耦合(TFC)分析方法，包括TFMI、TFNMI和TFCMI，其本质是计算时频能量谱间的MI、NMI和CMI。如图1所示，本发明的实施主要包括六个步骤：(1)多通道表面肌电信号同步采集并进行预处理；(2)对预处理后的数据进行噪声辅助的多元经验模态分解；(3)对IMF尺度分量进行同步提取变换；(4)计算时频互信息、时频归一化互信息和时频条件互信息；(5)采用步骤(4)的计算结果进行多尺度时频肌间耦合统计分析。

下面逐一对各步骤进行详细说明。

步骤一：多通道表面肌电信号同步采集并进行预处理

采用表面肌电设备同步采集肩膀3块肌肉上的sEMG信号，采样频率为2000Hz。预处理过程，首先手工提取5次有效活动段数据，然后通过上下采样的方式维持信号时长2.5s，接着进行去均值、去基线漂移，利用谱插值算法抑制50Hz工频干扰，最后采用4阶巴特沃斯带通滤波进行5-200Hz带通滤波，得到纯净的3通道sEMG信号，预处理后时域和频域效果如图2(a)和图2(b)所示。

步骤二：对预处理后的数据进行噪声辅助的多元经验模态分解

根据NA-MEMD算法，对m通道的sEMG信号X(t)添加2通道互不相关的高斯白噪声，组成多元复合信号

长度为T，经过多次MEMD分解，获得一系列IMF分量和余量的加和：

其中，q表示IMF尺度个数，

与

分别对应于m组信号和2组噪声的IMF矢量和残余矢量，将与噪声相关的2通道剔除，保留有用的m组IMF。

步骤三：对IMF尺度分量进行同步提取变换

对每个IMF分量

进行短时傅里叶变换(Short Time Fourier Transform,STFT)后再乘以一个相位因子e^jωt，可得

其中，k∈{1,2,...,m}，i∈{1,2,...,q}。利用式(3)，估计出STFT谱中时频系数的瞬时频率，

其中，

是STFT_e(t,ω)对时间的一阶偏导数。根据SST的思想，SET可以写为，

SET(t,ω)＝STFT_e(t,ω)·δ(ω-ω₀(t,ω))   (4)

其中，δ(ω-ω₀(t,ω))称为同步提取算子(SEO)，δ(*)为delta函数。

步骤四：计算时频互信息、时频归一化互信息和时频条件互信息

在同一IMF尺度分量上，计算第k个通道和第l个通道的sEMG信号的时频互信息(TFMI)、时频归一化互信息(TFNMI)和时频条件互信息(TFCMI)，即

TFMI(U,V)＝H(U)+H(V)-H(U,V)

TFCMI(U,V|W)＝H(U,W)+H(U,W)-H(W)-H(U,V,W)   (3)

其中，U、V、W分别为

其中，k≠l≠s∈{1,2,...,m},i∈{1,2,...,q}，H(*)表示时频熵，H(*,*)和H(*,*,*)表示时频联合熵，对数运算时基底取为2，单位为比特(Bit)。当U与V完全相关时，TFMI(U,V)最大，TFNMI(U,V)等于1，若W与U，W与V统计独立，则TFCMI(U,V|W)等于TFMI(U,V)；当U与V不相关时，TFMI(U,V)＝TFCMI(U,V|W)＝0Bit，TFNMI(U,V)＝0。在进行熵估计时，采用简单的直方图做法，即将时频能量值映射到不同数值区间(能量等级)上，再统计边际概率和联合概率。

为了验证本发明方法的性能，实验部分采集了多个受试者的数据，通过性能指标来衡量本发明在肌间耦合分析领域的可行性。实验招募了共10名受试者参与实验，对照组：5名健康被试(H1～H5)，平均年龄小于76岁；实验组：5名卒中患者(S1～S5)，中风后少于33.7天，显示轻度运动障碍，Fugl-Meyer平均评分高于43.732分(总分为66分)。要求每位受试者坐在桌子前，前臂放在一个舒服的位置上，根据每隔10秒激活的语音提示，要求向着正前方目标进行5次及物动作，目标距桌子高35cm，每次及物后休息10s，受试者在执行运动任务前已了解实验内容。采用表面肌电设备同步采集上肢前三角肌(Anterior Deltoid，AD)、内侧三角肌(Medial Deltoid，MD)、后三角肌(Posterior Deltoid，PD)上的sEMG信号，采样频率为2000Hz。

图3(a)、(b)为NA-MEMD算法分解获得的多个IMF尺度分量和残余分量。由图3(a)、(b)可知H1的3通道sEMG信号分解出6组IMF分量和1组残余分量，S1分解出5组IMF分量和1组残余分量，在实验时需要将残余分量进行丢弃。各IMF分量的频率分布从高到低依次排列，主频错开，可见NA-MEMD算法有效地避免了模态混叠和模式分离等问题。

由于个体差异性和sEMG信号的时变性，不同的被试或实验(trials)在NA-MEMD分解时可能存在IMF分量个数不一致的情况，造成分析困难，因此利用方差贡献率(VarianceContribution Rate,VCR)判定IMF分量的重要程度，挑选出累计方差贡献率大于90％的前几个IMF分量用于时频耦合分析。H1和S1上各个IMF分量的平均VCR如表1所示。

表1 IMF分量VCR(平均值±标准差)

由表1可见，AD、MD和PD的IM1～IMF4分量的平均VCR出现相同的分布规律，平均累积VCR已超过90％，其中，IMF2分量平均VCR最大，认为前4组IMF分量足够代表(解释)源信号绝大部分信息。

前4组IMF分量经SET时频变换后得到的平均时频能量谱如图4(a)和图4(b)所示，图中颜色深浅代表能量谱密度大小(dB)。由图可见，H1和S1的AD、MD和PD分解出的各IMF分量时频能量集中的位置接近，前0.5s能量值几乎为0。IMF1分量带宽范围大致在100～250Hz，IMF2分量带宽范围大致在50～150Hz，IMF3分量带宽范围大致在40～80Hz，IMF4分量带宽范围大致在20～60Hz。注意，这里为方便显示，对时频能量谱做对数变换后，幅值较低的频率成分有相对拉高。

取能量等级L为8，即将时频能量值映射到0～2⁷，AD、MD和PD在4个IMF分量间的平均时频耦合强度TFMI、TFNMI和TFCMI结果如图5(a)和图5(b)所示。由图可见，H1的肌间耦合强度TFMI、TFNMI和TFCMI值在IMF1～IMF4分量上都比较低且比较接近，TFMI和TFCMI在0.005～0.008Bit之间，TFNMI平均值从IMF1到IMF4有逐渐上升的趋势，在0.04～0.15之间；S1的肌间耦合强度和H1有所不同，TFMI和TFCMI平均值在0.003～0.017Bit之间，TFNMI平均值从IMF1到IMF4也有逐渐上升的趋势，在0.05～0.27之间，直观上要高于H1，尤其是MD与PD。并且，S1中MD与PD的TFMI、TFNMI和TFCMI平均值均要明显高于AD与MD、AD与PD。

更加具体地，H1～H5、S1～S5在不同IMF尺度分量上的肌间耦合强度TFMI、TFNMI和TFCMI结果如表2所示，单因素多元方差分析用于组间显著性检验。

表2不同IMF分量上TFMI、TFNMI和TFCMI结果(平均值±标准差，×10^-2)

由表2可见，在对照组中，除H1外，H2～H5的肌间耦合强度情况是相似的，从IMF1到IMF4，TFMI和TFCMI平均值有逐渐的下降趋势。TFMI范围在0.0022～0.0147Bit之间，TFCMI本质上是要低于或等于TFMI，范围在0.0017～0.0119Bit之间。然而，H1～H5的TFNMI从IMF1到IMF4都是呈现出逐渐上升的特点，范围在0.03～0.1845之间，此外，AD、MD和PD的肌间耦合强度TFMI、TFNMI和TFCMI在同一IMF分量上都十分接近。实验组的肌间耦合强度情况与对照组有所区别，TFMI范围在0.0030～0.0214Bit之间，TFCMI范围在0.0030～0.0182Bit之间。除S3外，TFNMI也是随着IMF1到IMF逐渐上升，范围在0.0367～0.3025之间，除S1外，MD与PD在IMF1～IM4分量上的肌间耦合强度TFNMI高于AD与MD、AD与PD的情况并不是很明显。单因素多元方差分析结果显示，组间AD、MD和PD的肌间耦合强度TFMI、TFNMI和TFCMI不具有显著性差异(p>0.05)。

以上所述的实施例仅仅是对本发明的优选实例方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (4)

1.多尺度时频肌间耦合分析方法，其特征在于：该方法包括以下步骤：

步骤(1)，多通道表面肌电信号同步采集并进行预处理；

步骤(2)，对预处理后的数据进行噪声辅助的多元经验模态分解，得到有用的IMF尺度分量；

步骤(3)，对IMF尺度分量进行同步提取变换；具体为：首先对每个IMF尺度分量进行短时傅里叶变换，再乘以一个相位因子，最后进行同步压缩变换；

步骤(4)，计算时频互信息、时频归一化互信息和时频条件互信息；

具体为：在同一IMF尺度分量上，计算第k个通道和第l个通道的表面肌电信号的时频互信息TFMI(U,V)、时频归一化互信息TFNMI(U,V)和时频条件互信息TFCMI(U,V|W)，即

TFMI(U,V)＝H(U)+H(V)-H(U,V)

TFCMI(U,V|W)＝H(U,W)+H(V,W)-H(W)

-H(U,V,W)

其中，U、V、W分别为：

其中，k≠l，k≠s，l≠s，k∈{1,2,...,m}，l∈{1,2,...,m}，s∈{1,2,...,m}，i∈{1,2,...,q}，SET(t,ω)表示同步压缩变换后的IMF尺度分量，H(*)表示时频熵，H(*,*)和H(*,*,*)表示时频联合熵，对数运算时基底取为2，单位为比特，q为IMF尺度分量数，m为表面肌电信号数；

步骤(5)，采用步骤(4)的计算结果进行多尺度时频肌间耦合统计分析。

2.根据权利要求1所述的多尺度时频肌间耦合分析方法，其特征在于：步骤(1)具体是：

采用表面肌电设备同步采集肩膀三块肌肉上的表面肌电信号，采样频率为2000Hz；

预处理过程：首先手工提取5次有效活动段数据，然后通过上下采样的方式维持信号时长2.5s，接着进行去均值、去基线漂移，利用谱插值算法抑制50Hz工频干扰，最后采用4阶巴特沃斯带通滤波进行5-200Hz带通滤波，得到纯净的3通道表面肌电信号。

3.根据权利要求1所述的多尺度时频肌间耦合分析方法，其特征在于：步骤(2)具体为：向预处理后的m通道表面肌电信号X＝{x₁,x₂,…,x_i,…x_m}中添加两路互不相关的高斯白噪声，组成多元复合信号X_c＝{x₁,x₂,…,x_m,x_m+1,x_m+2}，长度为T，经过多次多元经验模态分解，将与噪声相关的两路残余分量剔除，得到m×q个有用的IMF尺度分量。

4.根据权利要求1所述的多尺度时频肌间耦合分析方法，其特征在于：在进行熵估计时，采用直方图，即将时频能量值映射到不同数值区间上，再统计边际概率和联合概率。

Images