CN109674445B - 一种结合非负矩阵分解和复杂网络的肌间耦合分析方法 - Google Patents

Abstract

本发明提出了一种结合非负矩阵分解和复杂网络的肌间耦合分析方法；本发明采集了不同握力输出时，7通道表面肌电信号，首先采用广义偏定向相干计算多通道之间的相干性；然后用非负矩阵分解算法对相干性值进行分解；最后用复杂网络建立不同条件下的肌肉功能网络。利用图论的特征指标，评估肌肉功能网络信息流的传递效率。结果显示：不同握力下肌肉的激活程度存在显著性差异；肌肉间的耦合在10～20Hz波段上较为显著；并且在10～20Hz波段上其耦合程度随着握力水平不同呈现显著性变化。本发明通过对多通道肌间耦合分析，体现了中枢神经系统对不同输出握力的控制模式，为运动功能障碍的诊断和康复效果的评价提供依据。

Description

一种结合非负矩阵分解和复杂网络的肌间耦合分析方法

技术领域

本发明属于信号处理领域，涉及一种脑电、肌电信号的耦合性方法，特别涉及一种应用于肌间耦合状态分析方法。

背景技术

表面肌电信号(surface electromyography，sEMG)是肌肉运动过程中，肌肉收缩从而所产生的动作电位。其在皮肤表面叠加而成，与肢体的运动具有直接关系。肌间耦合(Intermuscular Coupling，IMC)是指肌肉之间的相互作用或相干性、时间同步性。肌间耦合可以间接揭示皮质下共同输入到肌肉的信息，有利于我们探究中枢神经系统对肌肉的协调控制。人类可以根据不同的需求完成不同的动作，其中抓握是人体上肢运动最为基础的动作之一，在完成抓握过程中需要多块肌肉的密切配合，尤其是在调节握力的大小中，肌肉之间的协调关系更为复杂。本发明对不同握力状态下的肌间耦合进行分析研究。近年来，利用肌间耦合研究运动过程中的肌肉协调状况特征，挖掘潜在的神经中枢系统运动控制机制的研究相继展开。有研究表明，肌间耦合存在频段显著特征，运动过程中肌间耦合主要表现在beta(15-30Hz)频段，beta频段的肌间耦合代表了从初级运动皮层到运动神经元的信息传递过程。Charissou Camille等通过基于小波的时频分析，研究在重复最大疲劳收缩前后膝关节最大等长伸展时beta频段肌间相互作用(IM)的调节，发现疲劳状态前后beta频段耦合特征存在差异。Dannados Santos A等采用肌间耦合对不同握力下三位抓握过程中手部肌肉之间协调影响进行分析，不同握力下的肌间一致性没有显著性差异，但是在不同的肌肉对之间的存在的显著性差异。杜文浩等构建变分模态分解-相干分析模型，分析20％的最大自主收缩力(maximum voluntary contraction,MVC)下特征频段的肌间耦合特性，并发现beta与gamma频段耦合强度随运动时间呈增强趋势。上述研究都只分析两块肌肉间的相干性，没有对多块肌肉之间的耦合关系进行深入探究。

针对这一问题，本发明提出了一种基于非负矩阵分解和复杂网络的肌间耦合分析方法，可以利用图论的知识，从整体来探索肌肉之间耦合关系及信息流动性变化。

发明内容

为了能够快速准确的分析上肢握力输出时，不同肌肉之间的耦合关系，本发明提出了一种基于非负矩阵分解和复杂功能网络的肌间耦合分析方法。采集了不同握力输出时，右手臂7通道的表面肌电信号(surface electromyography，sEMG)，首先利用采用广义偏定向相干(Generalized Partial Directed Coherence，GPDC)对预处理后的信号计算多通道之间的相干性；然后运用非负矩阵分解算法(Nonnegative Matrix Factorization，NMF)对计算得到的相干性值进行分解；最后运用复杂网络(Complex network，CN)建立不同条件下的肌肉功能网络；以确定其耦合关系。

本发明方法主要包括以下步骤：

一种结合非负矩阵分解和复杂网络的肌间耦合分析方法，该方法包括如下步骤：

步骤1、使用生物电采集设备采集5％、20％、40％、60％、80％最大自主收缩力抓握实验下的肌电信号；具体为：实验开始前先测量每位受试者的MVC，我们使用弹簧握力计分别进行5％、20％、40％、60％、80％MVC力输出的抓握实验，每种不同握力大小状态下恒力输出各3次，每次抓握5s，然后再放松20s；每种握力为一组，休息5分钟，然后做下一组握力实验。

步骤2、对采集到的肌电信号进行预处理，包括降频、降噪；

步骤3、采用广义偏定向相干方法对预处理后的信号计算多通道之间的相干性；

步骤4、运用非负矩阵分解算法对计算得到的相干性值进行分解，得到四个不同频段的相干性的值；

步骤5、根据步骤4的结果运用复杂网络建立不同条件下的肌肉功能网络。

作为优选，步骤3中的具体算法如下为：

广义偏定向相干方法是一种基于多变量自回归处理时间序列模型的格兰杰因果性频域测量方法；该函数分析多个通道序列在频域上的相关程度，用数值大小在0到1之间的归一化反映其相关程度大小；归一化的广义偏定向相干值

在[0，1]之间，表示x_g流向x_h的信号占所有从x_g流出信号的比例，其中x_g表示第g个通道，x_h表示第h个通道；具体研究算法，对于多通道的AR模型，表示为：

式中，x_m代表表面肌电信号sEMG的第m个通道，m表示通道数；P表示多变量自回归模型阶数，在处理肌电信号时选择P＝5，C_i为第i阶两两通道之间的相关系数：

然后对式(2)进行傅里叶变换，得式

式中I表示维数为m＝7的单位矩阵；这样

定义为式

式中，a_k,h(f)表示

中第h列的第k个元素，

表示

中第h列的第g个元素。

作为优选，步骤4中的具体算法如下为：

非负矩阵分解算法属于盲源分离方法的一种，由于NMF分解的结果具有非负的约束，其分解结果更具有实际意义，采用非负矩阵对相干性频率分解中，分解出相干性的四个主要作用频段；非负矩阵分解的基础形式为：

式中，W_ηω为基矩阵，C_ημ为系数矩阵，n为基矩阵分解列数；NMF的求解写成标准的优化问题：

求解优化问题(15)的迭代法则如下：

式(7)和(8)中的W和C交替优化，最终得到全局最优解；原矩阵M_ημ的列向量解释为C_ωμ中所有列向量的加权和，而权重系数为W_ij中对应列向量中的元素；将所有数据串联在一起进行非负矩阵分解，得到权重矩阵和系数矩阵，其中权重矩阵为所有受试者的相干性大小；系数矩阵为分解出来的主要频率段。

作为优选，步骤5中的具体算法如下为：

利用复杂网络，计算节点度、平均度、聚类系数、平均路径长度特征参数，进一步对肌肉功能网络的连接规律进行评价；通过对比相同握力不同频段间、特定频段不同握力下的肌肉功能网络拓扑结构的不同，研究探索大脑对肌肉的协调控制方式和调节机理；一个复杂网络图G(V,E)包括节点集合V和边的集合E；网络中各个节点之间的关系用邻接矩阵描述，矩阵中的元素A_gh为两个节点g,h之间的关系；在二值化网络中A_gh＝0时表示两个节点之间无连接边；A_gh＝1，表示有连接边，其特征参数如下：

节点度k_g表示与节点g相关联的边的条数，反映这个节点在网络中的活跃度，和节点对其相邻节点的影响力；节点度越大，节点对其相邻节点的重要度贡献越大，在网络中起着更重要的作用：

式中,m为网络中节点数目，A_gh为邻接矩阵的元素，平均度<k>表示整个网络中所有节点度的平均值：

聚类系数CC_g表示节点g聚集程度的系数,是描述网络集团化的重要指标；肌肉功能网络的聚类系数，反映了肌肉之间的功能性分化机理，即肌肉功能区域内的神经元之间连接紧密程度：

CC_g＝2e_g/k_g(k_g-1) (12)

式中,e_g为节点g与相连的k_g个节点间实际存在的边数；对网络中所有节点的聚类系数取平均值得到整个网络的平均聚类系数即：

最短路径长度表示两个节点i,j之间边数最少的一条通路的长度，平均路径长度L表示任意两个节点之间最短路径长度的平均值，不包括自身之间的连接，是研究复杂特性的关键参数：

式中,d_gh为从g到h之间的最短路径长度。

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

探究肌肉与肌肉之间的耦合关系时，传统的方法主要分析两块肌肉间的相干性，没有对多块肌肉之间的耦合关系进行深入探究。对于多通道的时间序列，不同的通道可能会受到其他通道相互作用，难以寻找较为准确耦合关系。这对这一问题，本发明提出了基于非负矩阵分解和复杂网络的肌间耦合分析，在探究耦合关系之前，先计算多通道肌电之间的相干性，然后运用非负矩阵分解相干性值，最后建立肌肉功能网络，分析不同握力下的网络特征。这将会从图论的角度分析，多块肌肉之间的耦合关系，反映中枢神经系统对不同输出握力的控制模式。

附图说明

图1为本发明的实施流程图。

图2为非负矩阵分解出的相干性及频率段示意图

图3为相同握力、不同频段下肌肉功能网络特征参数对比分析

图4为特定频段、不同握力下肌肉功能网络特征参数对比分析

具体实施方式

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

如图1所示，本实施专利包括如下步骤：

步骤1，采集了不同握力输出时，7通道表面肌电信号(EMG)。具体的采集操作如下：

选取6名无任何病史的右利手健康受试者，在安静的实验室内，受试者端坐于木制座椅上，将握力计放置于受试者手心。实验开始前我们先测量每位受试者的MVC，我们使用弹簧握力计分别进行5％、20％、40％、60％、80％MVC力输出的抓握实验，每种不同握力大小状态下恒力输出各3次，每次抓握5s，然后再放松20s；每种握力为一组，休息5分钟，然后做下一组握力实验。数据采集所用的是生物电采集系统(Brain Products GmbH,Germany)。主要采集右手臂右手肱二头肌(biceps brachi i，BB)、肱桡肌(brachioradialis，B)、桡侧腕长伸肌(Extensor carpi radialis，ECR)、指伸肌(extensor digitorum，ED)、尺侧腕伸肌(extensor carpi ulnaris，ECU)、指浅屈肌(flexor digitorum superficialis,FDS)、尺侧腕屈肌(Flexor carpi ulnaris，FCU)7通道的肌电信号(EMG)。在电极安放前，皮肤表面也经过酒精清理，采样频率为1000Hz。

步骤2，对步骤一中获取7通道表面肌电信号，进行消噪及滤波处理，得到分析所需的肌电信号。

步骤3，使用GPDC算法，计算多通道肌肉之间的相干性，得到相干性随频率变化的值。

步骤4，使用NMF分解前面得到的相干性的值，分解后得到系数矩阵和基矩阵，如图2所示，左侧为非负矩阵分解出来的主要频段范围(即系数矩阵)，右侧为三种不同握力下的相干性矩阵(即基矩阵)，每个方格代表两块肌肉之间的耦合性大小，白色代表耦合性较差，黑色代表耦合性较好。

步骤5，首先建立相同握力、不同频段下肌肉功能网络，通过分析比较复杂功能网络的3个重要特征参数(平均度、平均聚类系数、平均路径长度)，结果如图3所示，在特征频段(10～20Hz)内，肌间的耦合性更加密切，肌肉功能网络的节点间的连线数增加，网络的平均度也随之增加；相邻节点之间的联系变得更加紧密,平均聚类系数增加；节点之间的距离随之缩短，网络的平均路径长度减小从而使得全局效率有所提高。

然后以特征频段(10～20Hz)为主要分析对象，建立特征频段内、不同握力下的肌肉功能网络，利用图论分析，结果如图4所示，肌间耦合随着握力的增大而增大，肌肉功能网络的3个特征参数具有显著性变化，即肌肉之间的耦合程度在不同的握力状态下存在显著性不同。握力越大肌肉之间的协调作用越强烈，肌肉功能网络间的内部连通性越好,有利于维持较高的握力状态。

这些肌肉功能网络特征参数的不同，表明肌肉的耦合程度在相同握力不同频段间、特征频段不同握力间有显著性不同。在10～20Hz频段内以及维持较高握力时肌肉功能网络间的内部连通性变强,信息交流更为密切。

Claims (5)

1.一种结合非负矩阵分解和复杂网络的肌间耦合分析方法，其特征在于，该方法包括如下步骤：

步骤1、使用生物电采集设备采集5％、20％、40％、60％、80％最大自主收缩力抓握实验下的肌电信号；

步骤2、对采集到的肌电信号进行预处理，包括降频、降噪；

步骤3、采用广义偏定向相干方法对预处理后的信号计算多通道之间的相干性；

步骤4、运用非负矩阵分解算法对计算得到的相干性值进行分解，得到四个不同频段的相干性的值；

步骤5、根据步骤4的结果运用复杂网络建立相同握力不同频段以及特定频段不同握力下的肌肉功能网络。

2.根据权利要求1所述的结合非负矩阵分解和复杂网络的肌间耦合分析方法，步骤1具体为：实验开始前先测量每位受试者的MVC，我们使用弹簧握力计分别进行5％、20％、40％、60％、80％MVC力输出的抓握实验，每种不同握力大小状态下恒力输出各3次，每次抓握5s，然后再放松20s；每种握力为一组，休息5分钟，然后做下一组握力实验。

3.根据权利要求1所述的一种结合非负矩阵分解和复杂网络的肌间耦合分析方法，步骤3中的具体算法如下为：

广义偏定向相干方法是一种基于多变量自回归处理时间序列模型的格兰杰因果性频域测量方法；利用广义偏向相干方法分析多个通道序列在频域上的相关程度，用数值大小在0到1之间的归一化反映其相关程度大小；归一化的广义偏定向相干值

在[0，1]之间，表示x_g流向x_h的信号占所有从x_g流出信号的比例，其中x_g表示第g个通道，x_h表示第h个通道；具体研究算法，对于多通道的AR模型，表示为：

式中，x_m代表表面肌电信号sEMG的第m个通道，m表示通道数；P表示多变量自回归模型阶数，在处理肌电信号时选择P＝5，C_i为第i阶两两通道之间的相关系数：

然后对式(2)进行傅里叶变换，得式

式中I表示维数为m＝7的单位矩阵；这样

定义为式

式中，a_k,h(f)表示

中第h列的第k个元素，

表示

中第h列的第g个元素。

4.根据权利要求1所述的一种结合非负矩阵分解和复杂网络的肌间耦合分析方法，步骤4中的具体算法如下为：

非负矩阵分解算法属于盲源分离方法的一种，由于NMF分解的结果具有非负的约束，其分解结果更具有实际意义，采用非负矩阵对相干性频率分解中，分解出相干性的四个主要作用频段；非负矩阵分解的基础形式为：

式中，W_ηω为基矩阵，C_ημ为系数矩阵，n为基矩阵分解列数；NMF的求解写成标准的优化问题：

求解优化问题(15)的迭代法则如下：

式(7)和(8)中的W和C交替优化，最终得到全局最优解；原矩阵M_ημ的列向量解释为C_ωμ中所有列向量的加权和，而权重系数为W_ηω中对应列向量中的元素；将所有数据串联在一起进行非负矩阵分解，得到权重矩阵和系数矩阵，其中权重矩阵为所有受试者的相干性大小；系数矩阵为分解出来的主要频率段。

5.根据权利要求1所述的一种结合非负矩阵分解和复杂网络的肌间耦合分析方法，步骤5中的具体算法如下为：

利用复杂网络，计算节点度、平均度、聚类系数、平均路径长度特征参数，进一步对肌肉功能网络的连接规律进行评价；通过对比相同握力不同频段间、特定频段不同握力下的肌肉功能网络拓扑结构的不同，研究探索大脑对肌肉的协调控制方式和调节机理；一个复杂网络图G(V,E)包括节点集合V和边的集合E；网络中各个节点之间的关系用邻接矩阵描述，矩阵中的元素A_gh为两个节点g,h之间的关系；在二值化网络中A_gh＝0时表示两个节点之间无连接边；A_gh＝1，表示有连接边，其特征参数如下：

节点度k_g表示与节点g相关联的边的条数，反映这个节点在网络中的活跃度，和节点对其相邻节点的影响力；节点度越大，节点对其相邻节点的重要度贡献越大，在网络中起着更重要的作用：

式中,m为网络中节点数目，A_gh为邻接矩阵的元素，平均度<k>表示整个网络中所有节点度的平均值：

聚类系数CC_g表示节点g聚集程度的系数,是描述网络集团化的重要指标；肌肉功能网络的聚类系数，反映了肌肉之间的功能性分化机理，即肌肉功能区域内的神经元之间连接紧密程度：

CC_g＝2e_g/k_g(k_g-1) (12)

式中,e_g为节点g与相连的k_g个节点间实际存在的边数；对网络中所有节点的聚类系数取平均值得到整个网络的平均聚类系数即：

最短路径长度表示两个节点i,j之间边数最少的一条通路的长度，平均路径长度L表示任意两个节点之间最短路径长度的平均值，不包括自身之间的连接，是研究复杂特性的关键参数：

式中,d_gh为从g到h之间的最短路径长度。

Images