基于相干性—非负矩阵分解的多通道肌电耦合特性分析方法与流程

本发明涉及神经系统运动控制机制研究领域，尤其是一种基于相干性—非负矩阵分解的多通道肌电耦合特性分析方法。

背景技术：

肌间耦合是肢体在运动过程中不同肌肉间的相互关联与相互协调作用。通过研究多通道表面肌电信号(surfaceelectromyography，semg)间各特征频段的耦合特性，可以获得多通道肌肉间的功能联系及中枢神经系统支配肢体运动的执行与协调方式机理。近年来，基于传统的相干性分析方法研究肢体运动过程中肌肉间的耦合特性相继展开。有学者利用一致性分析方法计算两个肌电信号的互谱密度对信号自谱密度函数的归一化，以反映肌电信号在频域内的耦合关系。但是，传统的肌间一致性分析方法只能反映频域下的相干性，无法提取肌电信号在不同时频尺度下的特征信息，并且无法反映各频段的肌肉间功能连接强度。另外，在肢体运动过程中，多块肌肉同时动作，导致单个通道或者两个通道肌电信号无法全面反映肢体运动过程中肌肉间的功能耦合关系。

技术实现要素：

本发明目的在于提供一种可得到多通道肌电信号间的耦合特性、还能反映各个频段的肌肉间功能连接强度的基于相干性—非负矩阵分解的多通道肌电耦合特性分析方法。

为实现上述目的，采用了以下技术方案：本发明所述方法将相干性分析和非负矩阵分解相结合，具体步骤如下：

步骤1，同步采集多通道肌电信号并对其进行预处理；

步骤2，计算多通道肌电信号相干性；

步骤3，通过非负矩阵分解获得各个频段肌间功能连接强度。

进一步的，步骤1中，采集多通道肌电信号时，利用美国delsys公司trigno^tmwirelessemg采集设备，分辨率设为16bit，采样率为2000hz；采集信号前，被试静坐使上臂自然下垂，肘部用绷带固定在支架上，保证实验过程中姿势不变，调节支架使前臂与地面平行，前臂与上臂的夹角约为90^°，同时采集前臂旋前动作下的多块肌肉的肌电信号。

进一步的，步骤1中，对肌电信号进行预处理时，利用自适应50hz工频陷波滤波器对肌电信号进行处理，去除工频干扰；选用巴特沃斯三阶带通fir滤波器对肌电信号进行处理，使肌电信号主要集中在5-200hz之间。

进一步的，步骤3中，各频段肌间耦合特性分析方法的具体内容如下：

首先将多通道肌电信号进行相干性分析，然后利用非负矩阵分解方法将多通道肌电信号间的相干性结果进行分解，得到各个频段的相干性，进而定量分析多通道肌肉间耦合特性；

相干性能体现两个时间序列在频域上的相关程度，设x和y为两组时间序列，两信号相干性计算公式如下：

式中，sxy(f)为x和y在频率f上的互谱密度，sxx(f)、syy(f)分别为x和y的自谱密度；cxy为x和y的相干性，其取值范围为0-1；如果cxy(f)＝1，说明x和y在频率f上完全线性相关；如果cxy(f)＝0，说明x和y在频率f上完全独立；如果cxy(f)取值在0到1之间，说明x和y在频率f上部分线性相关，可能存在非线性关系；

非负矩阵分解方法的基本思想为：对于任意给定的非负矩阵vi×j，非负矩阵分解方法能够寻找一个非负矩阵wi×p和一个非负矩阵hp×j，使得满足

v≈wh(2)

或者

式中，矩阵vi×j为连接矩阵，矩阵wi×p为基矩阵，矩阵hp×j为系数矩阵；

在非负矩阵分解中，目标函数用来衡量分解结果的逼近程度；

对于欧几里得距离的目标函数的迭代的规则

目标函数||v-wh||²是单调的，但不是增函数，且||v-wh||²保持不变的条件是矩阵w和h固定；

对于k-l散度的目标函数的迭代的规则

目标函数d(v||wh)是单调的，但不是增函数，且d(v||wh)保持不变的条件是矩阵w和h固定；

协同个数p的值和肌电信号通道个数i、信号的时间序列长度j满足

(i+j)×p＜i×j(8)

由于分解前后的矩阵中仅包含非负元素，因此，原矩阵v的列向量可以解释为对基矩阵w中所有列向量的加权和，而权重系数为h中对应的列向量中的元素。

与现有技术相比，本发明具有如下优点：提出了一种用于多通道肌电耦合特性研究的新方法，定量刻画肌电信号在不同时频尺度上的功能耦合特征，进而定量分析多通道肌电信号间的耦合特性以及各个频段肌肉间功能连接强度，是进行特征频段肌间耦合特性分析的有效方法，为深入探索中枢神经系统运动控制机制提供有效的观察手段，同时也提供一种多通道肌电信号的研究方法。

附图说明

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

图2为本发明方法的肌电信号采集位置图。

图3为肌电信号预处理前后的对比图。

图4为被试的多通道肌电信号相干性—非负矩阵分解结果图。

具体实施方式

下面结合附图对本发明做进一步说明：

人体肌电信号随着肢体运动状态变化而改变，且每一个动作同时有多块肌肉参与其中。因此，分析多通道肌电信号间的耦合特性对中枢神经系统的运动控制机制研究尤为重要。本发明将传统相干分析方法和非负矩阵分解相结合，如图1所示，具体包括多通道肌电信号同步采集、肌电信号预处理、多通道肌电信号相干性计算、非负矩阵分解、各频段肌间耦合特性分析、功能耦合特性评价。具体方法包括肌电信号采集部分和信号处理部分：

肌电信号采集，利用美国delsys公司trigno^tmwirelessemg采集设备，分辨率设为16bit，采样率为2000hz。采集信号前，被试静坐使上臂自然下垂，肘部用绷带固定在支架上，以保证实验过程中姿势不变，调节支架使前臂与地面平行，前臂与上臂的夹角约为90°，放置电极前用75％酒精擦拭被测皮肤的表面，去除皮肤表面油脂和皮屑，同步采集前臂旋前动作下的多块肌肉的肌电信号，如图2所示，具体包括被试右臂的指浅屈肌(fds)、尺侧腕伸肌(ecu)、指伸肌(ed)、桡侧腕伸肌(ecr)、桡侧腕屈肌(fcr)、掌长肌(pl)、肱二头肌(bb)和肱桡肌(b)。

信号处理部分包括预处理和各频段肌间耦合特性分析

信号预处理：肌电信号是一种微弱信号，容易受到噪声的干扰，采集到的多通道原始肌电信号需要进行预处理，利用自适应50hz工频陷波滤波器对肌电信号进行处理，去除工频干扰；选用巴特沃斯三阶带通fir滤波器对肌电信号进行处理，使肌电信号主要集中在5-200hz之间。信号预处理前后对比如图3所示。从图3可以看出，预处理有效滤除了原始肌电信号中的50hz工频干扰及其倍频干扰。

各频段肌间耦合特性分析

本发明首先对各个通道肌电信号进行相干性分析，然后利用非负矩阵分解方法将肌间相干性结果进行分解，得到各个频段的相干性，进而定量分析肌肉间耦合特性。

相干性可以体现两个时间序列在频域上的相关程度，设x和y为两组时间序列，两信号相干性计算公式如下：

式中，sxy(f)为x和y在频率f上的互谱密度，sxx(f)、syy(f)分别为x和y的自谱密度。cxy为x和y的相干性，其取值范围为0-1；如果cxy(f)＝1，说明x和y在频率f上完全线性相关；如果cxy(f)＝0，说明x和y在频率f上完全独立；如果cxy(f)取值在0到1之间，说明x和y在频率f上部分线性相关，可能存在非线性关系。

进行相干性分析之后，要对相干性结果进行非负矩阵分解，进而分析多通道各频段的肌间耦合特性。

nmf的基本思想可以简单的描述为：对于任意给定的非负矩阵vi×j，nmf能够寻找一个非负矩阵wi×p和一个非负矩阵hp×j，使得满足

vi×j≈wi×php×j(2)

或者

式中，矩阵vi×j为连接矩阵，矩阵wi×p为基矩阵，矩阵hp×j为系数矩阵。

在非负矩阵分解中，目标函数的选择至关重要，利用它来衡量分解结果的逼近程度。非负矩阵分解算法目标函数的选取有许多，其中较为常用的是基于k-l散度(kullback-leiblerdivergence)和欧几里得距离(eulideandistance)的目标函数。

基于欧几里得距离的目标函数为：

式中，目标函数||v-wh||²取得最小值的条件是v＝wh，且最小值为0。

基于k-l散度的目标函数为：

式中，目标函数d(v||wh)取得最小值的条件是v＝wh，且最小值为0。

在最优化问题中，矩阵w和h均是变量，不管选择哪种目标函数，矩阵w和h均不是凸函数，因而求得其最优解比较困难。为此，采取以下迭代规则，既能保证运算速度，又能方便运算。

对于欧几里得距离的目标函数的迭代规则

目标函数||v-wh||²是单调的，但不是增函数，且||v-wh||²保持不变的条件是矩阵w和h固定。

对于k-l散度的目标函数的迭代规则

目标函数d(v||wh)是单调的，但不是增函数，且d(v||wh)保持不变的条件是矩阵w和h固定。

协同个数p的值一般要用适合的方法进行严格的选择，且和肌电信号通道个数i、信号的时间序列长度j满足

(i+j)×p＜i×j(10)

由于分解前后的矩阵中仅包含非负元素，因此，原矩阵v的列向量可以解释为对基矩阵w中所有列向量的加权和，而权重系数为h中对应的列向量中的元素。这种基于基向量组合的表示形式具有很直观的语义解释，它反映了人类思维中“局部构成整体”的概念。

为验证本发明所述的一种基于相干性—非负矩阵分解的多通道肌电耦合特性分析方法，采集6名健康被试(年龄为(25±3)岁)的上肢semg，被试相关信息如表1所示。

表1被试相关信息

要求被试实验前无肌肉疲劳现象、精神状态良好，且熟悉实验流程。按照本发明所述的多通道肌电信号采集与处理过程，采集健康被试的上肢表面肌电信号，并分析多通道肌电的耦合特性，进而研究中枢神经系统的运动控制机制。

图4为被试多通道肌电信号相干性-非负矩阵分解结果图(不同颜色代表不同的耦合强度)。图4中，左侧图横坐标为肌电信号频率、纵坐标为信号频谱，右侧图横、纵坐标为被试上肢8块肌肉；由图4(a)～(d)图可以直观看出，被试上肢的semg信号被分解到4个频段内，且通过方格颜色的不同代表不同肌肉间的耦合强度大小，体现出肌间耦合强度的差异。因此，多通道肌电信号经过相干性—非负矩阵分解可以得到不同频段肌电信号对应的肌间耦合特性强度，为深入探索中枢神经系统运动控制机制提供有效的观察数据。

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