CN111931129A - 基于高斯Copula传递熵的肌间耦合网络分析方法 - Google Patents

Abstract

本发明公开了一种基于高斯Copula传递熵的肌间耦合网络分析方法。本发明首先采用高斯Copula函数估计传递熵，然后采集了8名健康被试和5名卒中患者执行及物运动任务时5通道的表面肌电信号，利用快速傅里叶变换提取出α、β和γ频段信息，进而采用高斯Copula传递熵作为测度建立肌间耦合网络，最后通过网络拓扑参数分析肌间特征频段的功能耦合和信息流向特征。本发明避免了对联合概率密度的估计，能更好地推断出复杂网络中因果关系，对于挖掘潜在的运动控制机制和运动功能康复评估具有良好的应用价值。

Description

基于高斯Copula传递熵的肌间耦合网络分析方法

技术领域

本发明属于神经系统运动控制机制研究领域，涉及高斯Copula和传递熵的计算，从而进行肌间因果耦合网络分析。

背景技术

脑卒中又称为中风，脑血管意外，是在脑血管病变或血流障碍基础上发生的局限性或弥漫性脑功能障碍，业已成为成人获得性运动障碍的主要原因。目前在脑卒中患者康复过程中，运动功能评价的主要依据是康复医师的主观经验或依靠监测患者的肌力状态，难以客观、准确、定量评定康复效果。肌间耦合的概念来源于皮层肌肉功能耦合研究，指的是运动过程中肌肉间的相互作用。肌间耦合在人体运动中起着重要的作用，决定了人体运动的整体肌肉模式，能够反映肌肉系统在中枢神经系统控制下的运动功能状态和信息交互方式。由于肌肉活动代表了神经系统的输出，对肌间耦合状态的检查能够了解运动障碍患者神经机制的灵活性和适应性的差异。

复杂网络理论作为一种研究复杂系统动力学的新方法最早由Watts和Barabasi等人提出，复杂网络分析是把复杂系统的内部单元抽象为网络的节点，单元之间以边的形式相连接，通过量化分析网络的拓扑结构特性及动力学行为，可以揭示复杂系统的内在属性及运行规律等重要信息。复杂网络的局部和全局特性能够清晰地刻画组成复杂系统的不同元素之间的相互关系和信息流动过程。近年来，复杂网络吸引了相关研究人员的广泛关注，尤其是复杂脑网络的研究，已成为神经科学的热点。大量研究表明，复杂网络既不是规则网络，也不是随机网络，具有小世界、无标度特性。目前，复杂网络的研究视角主要包括结构性网络、功能性网络和因效性网络。不同于无向连接的功能性网络，因效性网络重点强调网络中各种连接的方向性，着重分析各网络节点之间的因果关系以及统计趋势，并根据信息在节点之间的传播方向来分析网络演化的工作过程。因效性网络和功能性网络的差别在于如何量化测度网络节点之间的关系，一般采用因果关系分析来对节点的连接强度进行量化。

传统的因果性分析方法中，大量研究者使用格兰杰因果性(Granger Causality,GC)和传递熵(Transfer Entropy,TE)来分析复杂的数据，特别是具有非线性、多维度特性的数据间的因果关系，往往能得到问题的解释。GC和TE都是基于相同的因果关系定义的，Barnett等证明了在高斯假设下，GC和TE等价，建立在线性回归模型上的GC分析是从预测的角度来检测变量间的因果性，可描述双向因果传递、计算相对简单。Samani等利用GC研究了6块手部肌肉对之间功能连接随时间变化情况，发现只有三角肌前束至肱二头肌的因果关系随时间增加(三角肌前肌束→肱二头肌)，提示这可能与疲劳发展的共收缩效应相关。但GC存在两个问题：一是，模型阶数难以确定，频率分辨率受限；二是，难以定量描述高阶复杂的因果关系。TE是基于转移概率定义的非对称性测量方法，包含着方向性和动态信息，具有不依赖既定模型、非线性定量分析的特点，本质上等价于条件互信息(Conditional MutualInformation,CMI)。Xie等提出了一种基于TE的多尺度传递熵(MSTE)方法，利用MSTE方法计算不同频段、不同尺度上的脑肌电功能耦合特征，发现MSTE可以定量估计感觉运动皮层与肌肉之间的非线性连接。但是，传统的TE需要估计变量的概率密度函数，计算复杂度高且不稳定。

Sklar定理指出，一个N维分量的联合分布函数可以由N个变量的边缘分布和1个Copula函数来描述。这个Copula函数译为“连接函数”或“相依函数”，是把多个随机变量的联合分布与各自的边缘分布相连接的函数，既能够完整地刻画变量间的相关结构，又可以将单个随机变量的边际分布与变量间的相关结构拆开来处理，具有不受变量维度限制、对边际分布没有特定要求等优点。在Copula理论中，Copula熵被定义为包含随机变量之间所有依赖信息的Copula函数的统计依赖度的度量，具有多元、对称、非负、不受单调变换影响、在高斯情形下等价于相关系数等性质，被认为是理想统计独立测度应具有的性质。Ma等证明了负Copula熵等价为互信息(Mutual Information,MI)，提出利用Copula函数理解和估计MI的新方法，发现通过Copula熵估计互信息方法简单，计算量小，能有效避免对联合密度函数地估计，这为估计TE提供了一种新的参考。实际上，Copula也已被用来估计因果关系，Hu和Liang等提出用Copula来计算GC，将GC的定义写成条件Copula密度的形式，然后根据经验估计条件Copula密度，但该方法缺乏收敛性保证。

针对复杂网络中TE估计困难等问题，本发明提出了一种高斯Copula传递熵的新方法，即利用高斯Copula函数去估计TE，并应用多特征功能频段的肌间耦合网络分析，旨在探索中枢神经系统的运动控制机制以及为运动功能康复评估提供新的研究方法和理论依据。

发明内容

本发明的目的在于从复杂网络的角度提供一种可准确得到上肢肌间因果耦合特性的分析方法。

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

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

具体为：在表面肌电设备的监控下，同步采集上肢前三角肌(Anterior Deltoid，AD)、内侧三角肌(Medial Deltoid，MD)、后三角肌(Posterior Deltoid，PD)、肱二头肌(Biceps，BB)、肱三头肌(Triceps，TB)上的5通道sEMG信号，采样频率为2000Hz。鉴于sEMG信号在低频段耦合特征显著，所以本发明重点关注α(8～15Hz)、β(15～30Hz)和γ(30～60Hz)三个特征功能频段。

步骤(2)，高斯Copula传递熵的定义与计算；

具体为：传递熵是基于转移概率定义的非对称性测量方法，包含着方向性和动态信息，具有不依赖既定模型、非线性定量分析的特点，本质上等价于条件互信息(Conditional Mutual Information,CMI)。然而，传递熵的估计一直是一个具有数值挑战性的问题，一般来说，这种估计依赖于相关变量的联合概率分布的精确表示。本发明提出一种新的高斯Copula传递熵方法，即利用高斯Copula函数去估计传递熵。假设两通道sEMG信号X＝{x₁,x₂,...,x_T}，Y＝{y₁,y₂,...,y_T}分别为一阶马尔科夫过程，T为时间序列的长度，X→Y的高斯Copula传递熵定义为

其中，τ和t分别为X和Y的滞后时间，

k和l分别为X和Y的嵌入维度，ρ₁，ρ₂和ρ₃分别为高斯Copula函数参数，H_c(*)表示变量间的高斯Copula熵。由高斯Copula传递熵的表达式可知，其关键在于3项高斯Copula熵的估计，为此，本发明提出一种简单可行的高斯Copula熵估计方法，包含以下3个步骤：

(a)非参数核密度估计边际分布函数

假设各通道sEMG信号

是来自连续分布函数F_i(x_i)的同分布样本，T为时间序列的长度，i＝1,2,...,5，那么F_i(x_i)的非参数核密度估计为

其中，

为概率密度函数，

其中，h为窗宽，当T→∞，h→0，Th→∞时，非参数核密度估计是真实概率密度分布的一致估计。K(·)为核函数，本发明采用高斯核函数，根据经验法则确定窗宽：

为x_i的标准差。

(b)典型极大似然法估计高斯Copula参数

常用的Copula参数估计方法有最大似然估计(ML)、两阶段极大似然估计(IFM)和典型极大似然估计(CML)。由于CML估计在边际分布未知时，相对ML估计和IFM估计稳定性更好，所以本发明选择CML估计高斯Copula参数。根据CML估计方法，在利用非参数核密度估计边际分布函数后，可以无需估计边际分布中的参数，只要估计Copula函数中的参数ρ

(c)蒙特卡洛法估计高斯Copula熵

高斯Copula熵的计算可直接通过多重积分法求得，但当变量较多时，多重积分法较为困难，本发明采用蒙特卡洛法模拟高斯Copula熵，表示为

H_c(x₁,x₂,...,x_d)＝-E[logc(u₁,u₂,...,u_d)] (5)

其中，E[·]表示数学期望，当对数运算以2为基底时，高斯Copula传递熵的单位为比特(Bit)。

步骤(3)，肌间耦合网络的搭建与分析；

具体为：复杂网络是研究自然界各种复杂系统的有力工具，肌间耦合网络的建立有助于研究肌肉系统之间的相互作用、网络拓扑信息以及与肌肉功能、疾病产生的联系。一般而言，肌间耦合网络的搭建及分析主要包含以下3个步骤：

(a)定义网络的节点和边

肌间耦合网络本质上是有向加权复杂网络，可用G表示，G＝(V,E)。V＝{v₁,v₂,...,v_n}为肌肉节点集合，

为边集合。肌间耦合网络的节点数目为n＝|V|，边数为m＝|E|。v_i∈V,(i＝1,2,...,n)表示网络中的一个节点，v_i→v_j∈E表示节点v_i到节点v_j的一条有向边。

表示有向边v_i→v_j的权值(或耦合强度)，一般情况下，

(b)确定网络连接阈值

为尽可能满足肌间耦合网络的连通性，避免孤立节点的出现，本文设置阈值在0.1～0.9*max(GCTE)变化，其中max(GCTE)表示网络中GCTE的最大值，将肌间耦合网络的拉普拉斯矩阵的第二小非零特征值λ_min所对应的阈值确定为最佳阈值。简言之，设第k+1个阈值对应的邻接矩阵为GCTE^k+1，k∈[1,9]，肌间耦合网络的拉普拉斯矩阵L的第二小特征值λ_min等于0

L＝D-GCTE^k+1 (6)

其中，D为肌间耦合网络的度矩阵，对角元素D_ii为GCTE^k+1的第i列元素之和，其余元素为零。则我们将取第k个阈值为最佳阈值，阈值化后的邻接矩阵为

(c)分析静态网络参数

进一步，利用图论中度(degree)和聚类系数(clustering coefficient)等静态网络参数，来分析各频段上肌间耦合网络的拓扑属性。具体分为两个方面：一是面向节点的节点度和节点聚类系数，反映网络的局部特性；二是面向网络的网络平均度和聚类系数，反映网络的全局特征。

1)节点度和聚类系数

有向网络中节点度可分为入度和出度，节点度越大意味着这个节点越“重要”。节点的聚类系数表示某一结点的邻居间互为邻居的可能，衡量的是网络集团程度。依据肌间耦合网络，定义每个肌肉节点v_i的入度和出度分别为D_in(v_i)和D_out(v_i)，聚类系数为CC(v_i)

其中，D_in(v_i)越大，表明该节点在网络中是多个节点共同作用的结果，影响该节点的连接越多。D_out(v_i)越大，表明该节点在网络中所起作用越大，是多个节点的“因”。当节点之间具有相同的入度和出度时，可认为这些节点是同质的。CC(v_i)越大，表明该节点与相邻节点之间存在规则且频繁的转换。

2)网络的平均度和聚类系数

网络的平均度描述了网络整体的连通程度，而网络的聚类系数度量了网络的积聚情况，定义肌间耦合网络的平均度AD为所有节点的入度和出度的平均值，网络的聚类系数NCC为所有节点的聚类系数的平均值

其中，AD越大，网络中连边越多，连通性越好。NCC∈[0,1]，当NCC＝0，表示网络中所有节点都是孤立节点，当NCC＝1，表示网络中任意节点间都有边相连。

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

传统的肌间耦合分析方法侧重两两通道之间功能耦合分析，缺乏对多通道信号之间信息定向传递的有效描述，针对该问题，本发明提出一种新的高斯Copula传递熵估计方法，并基于高斯Copula传递熵构建肌间耦合网络模型，定量描述上肢运动过程中多通道肌肉之间的因果耦合特性，这为探索运动过程肌肉系统耦合特征及中枢神经系统的运动控制机制提供新方法。

附图说明

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

图2(a)、图2(b)分别为受试者H1和S1的5通道sEMG信号分离出的α、β和γ频段波形。

图3为时延τ对GCTE性能的影响曲线。

图4(a)为不同阈值下的肌间耦合网络，图4(b)为肌间耦合网络在不同阈值下的拉普拉斯矩阵的第二小特征值，圆点对应的横坐标即为选取的最佳阈值。

图5(a)、图5(b)分别为健康组和卒中组在α、β和γ频段下肌间耦合网络对比。

图6(a)为健康组和卒中组分别在α、β和γ频段下肌间耦合网络肌肉节点的入度；图6(b)为健康组和卒中组分别在α、β和γ频段下肌间耦合网络肌肉节点的出度；图6(c)为健康组和卒中组分别在α、β和γ频段下肌间耦合网络肌肉节点的聚类系数。

图7(a)为健康组和卒中组分别在α、β和γ频段下肌间耦合网络的平均度；图7(b)为健康组和卒中组分别在α、β和γ频段下肌间耦合网络的聚类系数。

具体实施方式

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

卒中后，中枢神经系统受损，大脑对低级中枢(脊髓、脑干)调节和控制丧失，导致原始反射出现，运动功能模式异常。表现为肌张力增高，肌群间协调紊乱，甚至痉挛，严重影响患者的日常生活活动。肌间耦合可以反映来自脊髓中间神经元的共同驱动，体现运动过程中肌肉间的相互作用以及中枢神经对肌肉的支配能力。探讨肌肉系统的同步耦合特性和功能联系，对于指导卒中患者上肢运动功能康复评定具有十分重要的意义。

针对现有肌间耦合分析方法存在的问题，本发明将统计学中的Copula理论引入到肌间耦合网络分析之中。如图1所示，本发明的实施主要包括三个步骤：(1)多通道表面肌电信号的同步采集与预处理；(2)高斯Copula传递熵的定义与计算；(3)肌间耦合网络的搭建与分析。

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

步骤一：多通道表面肌电信号的同步采集与预处理

采用表面肌电设备同步采集及物运动过程中上肢前三角肌(Anterior Deltoid，AD)、内侧三角肌(Medial Deltoid，MD)、后三角肌(Posterior Deltoid，PD)、肱二头肌(Biceps，BB)、肱三头肌(Triceps，TB)上的5通道sEMG信号，采样频率为2000Hz。预处理过程，首先手工提取5次有效活动段数据，然后通过上下采样的方式维持信号时长2.5s，接着进行去均值、去基线漂移，利用谱插值算法抑制50Hz工频干扰，最后采用4阶巴特沃斯带通滤波进行5-200Hz带通滤波，得到纯净的3通道sEMG信号。由于sEMG信号非平稳、频域特征突出，在低频段耦合特征存在明显差异，所以本发明重点关注α、β和γ这3个特征功能频段，通过快速傅里叶变换提取的频段波形如图2(a)和图2(b)所示。

步骤二：高斯Copula传递熵的定义与计算

传递熵是基于转移概率定义的非对称性测量方法，包含着方向性和动态信息，具有不依赖既定模型、非线性定量分析的特点，本质上等价于条件互信息(ConditionalMutual Information,CMI)。然而，传递熵的估计一直是一个具有数值挑战性的问题，一般来说，这种估计依赖于相关变量的联合概率分布的精确表示。实际上，传递熵可理解为互信息的不对称形式，而Copula熵则被认为是互信息的一个扩展。受到这一思想的启发，本发明提出一种新的高斯Copula传递熵方法，即利用高斯Copula函数去估计传递熵。假设两通道sEMG信号X＝{x₁,x₂,...,x_T}，Y＝{y₁,y₂,...,y_T}分别为一阶马尔科夫过程，T为时间序列的长度，X→Y的高斯Copula传递熵定义为

其中，τ和t分别为X和Y的滞后时间，

k和l分别为X和Y的嵌入维度，ρ₁，ρ₂和ρ₃分别为高斯Copula函数参数，H_c(*)表示变量间的高斯Copula熵。由高斯Copula传递熵的表达式可知，其关键在于3项高斯Copula熵的估计，为此，本发明提出一种简单可行的高斯Copula熵估计方法，包含以下3个步骤：

(a)非参数核密度估计边际分布函数

假设各通道sEMG信号

是来自连续分布函数F_i(x_i)的同分布样本，T为时间序列的长度，i＝1,2,...,5，那么F_i(x_i)的非参数核密度估计为

其中，

为概率密度函数，

其中，h为窗宽，当T→∞，h→0，Th→∞时，非参数核密度估计是真实概率密度分布的一致估计。K(·)为核函数，本发明采用高斯核函数，根据经验法则确定窗宽：

为x_i的标准差。

(b)典型极大似然法估计高斯Copula参数

常用的Copula参数估计方法有最大似然估计(ML)、两阶段极大似然估计(IFM)和典型极大似然估计(CML)。由于CML估计在边际分布未知时，相对ML估计和IFM估计稳定性更好，所以本发明选择CML估计高斯Copula参数。根据CML估计方法，在利用非参数核密度估计边际分布函数后，可以无需估计边际分布中的参数，只要估计Copula函数中的参数ρ

其中c()为Copula密度函数，u₁～u_d表示d个产量(d维)的累积分布函数。

(c)蒙特卡洛法估计高斯Copula熵

高斯Copula熵的计算可直接通过多重积分法求得，但当变量较多时，多重积分法较为困难，本发明采用蒙特卡洛法模拟高斯Copula熵，表示为

H_c(x₁,x₂,...,x_d)＝-E[logc(u₁,u₂,...,u_d)] (17)

其中，E[·]表示数学期望，当对数运算以2为基底时，高斯Copula传递熵的单位为比特(Bit)。

步骤三：肌间耦合网络的搭建与分析

复杂网络是研究自然界各种复杂系统的有力工具，肌间耦合网络的建立有助于研究肌肉系统之间的相互作用、网络拓扑信息以及与肌肉功能、疾病产生的联系。一般而言，肌间耦合网络的搭建及分析主要包含以下3个步骤：

(a)定义网络的节点和边

肌间耦合网络本质上是有向加权复杂网络，可用G表示，G＝(V,E)。V＝{v₁,v₂,...,v_n}为肌肉节点集合，

为边集合。肌间耦合网络的节点数目为n＝|V|，边数为m＝|E|。v_i∈V,(i＝1,2,...,n)表示网络中的一个节点，v_i→v_j∈E表示节点v_i到节点v_j的一条有向边。

表示有向边v_i→v_j的权值(或耦合强度)，一般情况下，

(b)确定网络连接阈值

为尽可能满足肌间耦合网络的连通性，避免孤立节点的出现，本文设置阈值在0.1～0.9*max(GCTE)变化，其中max(GCTE)表示网络中GCTE的最大值，将肌间耦合网络的拉普拉斯矩阵的第二小非零特征值λ_min所对应的阈值确定为最佳阈值。简言之，设第k+1个阈值对应的邻接矩阵为GCTE^k+1，k∈[1,9]，肌间耦合网络的拉普拉斯矩阵L的第二小特征值λ_min等于0

L＝D-GCTE^k+1 (18)

其中，D为肌间耦合网络的度矩阵，对角元素D_ii为GCTE^k+1的第i列元素之和，其余元素为零。则我们将取第k个阈值为最佳阈值，阈值化后的邻接矩阵为

(c)分析静态网络参数

进一步，利用图论中度(degree)和聚类系数(clustering coefficient)等静态网络参数，来分析各频段上肌间耦合网络的拓扑属性。具体分为两个方面：一是面向节点的节点度和节点聚类系数，反映网络的局部特性；二是面向网络的网络平均度和聚类系数，反映网络的全局特征。

1)节点度和聚类系数

有向网络中节点度可分为入度和出度，节点度越大意味着这个节点越“重要”。节点的聚类系数表示某一结点的邻居间互为邻居的可能，衡量的是网络集团程度。依据肌间耦合网络，定义每个肌肉节点v_i的入度和出度分别为D_in(v_i)和D_out(v_i)，聚类系数为CC(v_i)

其中，D_in(v_i)越大，表明该节点在网络中是多个节点共同作用的结果，影响该节点的连接越多。D_out(v_i)越大，表明该节点在网络中所起作用越大，是多个节点的“因”。当节点之间具有相同的入度和出度时，可认为这些节点是同质的。CC(v_i)越大，表明该节点与相邻节点之间存在规则且频繁的转换。

2)网络的平均度和聚类系数

网络的平均度描述了网络整体的连通程度，而网络的聚类系数度量了网络的积聚情况，定义肌间耦合网络的平均度AD为所有节点的入度和出度的平均值，网络的聚类系数NCC为所有节点的聚类系数的平均值

其中，AD越大，网络中连边越多，连通性越好。NCC∈[0,1]，当NCC＝0，表示网络中所有节点都是孤立节点，当NCC＝1，表示网络中任意节点间都有边相连。

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

根据GCTE的表达式可知，GCTE是以时延τ为自变量，时延τ的选取直接影响着GCTE度量肌间因果耦合关系的准确性。有文献指出在TE取最大时，对应于真正的延迟时间，这可以作为确定时延τ的参考标准。本发明在计算两两通道之间的GCTE值时，设置时延τ在0～50ms变化，并取GCTE在此范围内的最大值，作为肌间因果耦合强度指标。关于时延τ的确定，以受试者H1的MD与TB在α频段的GCTE为例，如图3所示。从图3可以看出，MD→TB的GCTE在MD取时延τ为21ms时，达到最大值0.0950Bit，TB→MD的GCTE在TB取时延τ为19.5ms时，达到最大值0.0697Bit。这说明在α频段当MD与TB取相近的时延τ时，互有最大的信息传递，且MD对TB的因果耦合强度要高于TB对MD，当不同因果耦合方向上GCTE达到最大值时，对应的时延τ为最优参数。其他受试者在α、β和γ频段关于GCTE的时延τ的确定规则与之类似。

阈值的选取直接影响着以GCTE所表示肌间耦合网络的拓扑性质，为充分考虑肌间耦合网络的代数连通性，本发明通过判断肌间耦合网络的拉普拉斯矩阵的第二小特征值λ_min是否为零来确定最佳网络阈值，仍以受试者H1在不同阈值下α频段所构建的肌间耦合网络为例进行说明，由图4(a)可以看出，随着阈值的增加，肌间耦合网络的边数、节点的入度与出度均出现了明显的下降，阈值越大，肌间耦合网络中节点的连边越少，连通性越差。当阈值高于0.5*max(GCTE)时，肌间耦合网络开始出现孤立节点BB和PD，虽然较强的肌间因果耦合关系(如AD与MD)得以保留，但较弱的肌间因果耦合关系却难以察觉。在图4(b)中，随着阈值的增加，肌间耦合网络的拉普拉斯矩阵的第二小特征值λ_min从非零值(>0)逐渐下降，当阈值高于0.6*max(GCTE)以后，λ_min始终取零，刚好与图4(a)中的肌间耦合网络的连通情况相对应，这说明拉普拉斯矩阵的第二小特征值λ_min可以作为肌间耦合网络阈值选取的指标，最佳阈值可选择为0.5*max(GCTE)。其他受试者在α、β和γ频段关于肌间耦合网络阈值的选取规则与上述类似。

在确定好GCTE的时延和网络的阈值后，即可得到肌间耦合网络的邻接矩阵。为便于组间比较，将每位受试者的邻接矩阵进行组内平均，构建3个频段下的肌间耦合网络如图5(a)和图5(b)所示。肌间耦合网络中各肌肉节点间的信息流动是双向的，健康组和卒中组在α、β和γ频段下的肌间耦合网络不能发现明显的区别。直观上，MD、PD和TB拥有更高的入度和出度，并且MD与PD，PD与TB在3个频段上均存在着较强的信息交互。AD在α频段上的入度和出度要高于β和γ频段，同时AD与MD，AD与PD，AD与TB在α频段上的信息交互也要强于β和γ频段。BB在3个频段上的入度和出度都明显低于其他肌肉节点，一旦网络的阈值增加，BB势必被认为是网络中的孤立节点，并且BB与其他肌肉节点间几乎不存在任何信息交互，表现为肌间耦合网络中较细的连线。这或许说明MD、PD和TB在上肢及物运动过程中占据着主导作用，彼此间协同工作，信息交流密切，AD与其他肌肉间的耦合主要发生在α频段，在β和γ频段耦合程度较低，而BB可能独立于其他肌肉，其功能比较特殊。

进一步，计算3个频段下肌间耦合网络中各肌肉节点的入度、出度和聚类系数，并利用配对样本t检验进行组间显著性检验，显著性水平取为0.05，结果如图6(a)、图6(b)、图6(c)和表1所示。

表1α、β和γ频段下肌间耦合网络中节点的入度D_in、出度D_out和聚类系数CC对比(平均值±标准差)

3个频段下各肌肉节点的入度、出度和聚类系数的分布十分相似，MD在α频段拥有更高的入度和出度，PD在β和γ频段拥有更高的入度和出度。配对样本t检验结果显示，除BB在γ频段的出度存在显著性差异外(p＝0.026<0.05)，其他情况下组间各肌肉节点的入度、出度和聚类系数都不具有显著性差异(p>0.05)，其中，健康组的BB在γ频段的出度为0.03±0.03Bit，卒中组的BB在γ频段的出度为0.08±0.05Bit。这说明卒中组的AD、MD、PD和TB与健康人的肌间耦合情况相似，而其运动功能缺陷可能是由于BB在γ频段出现了异常耦合，表现为过多的信息释放。

利用单因素一元方差分析检验肌肉节点之间网络参数是否具有显著性差异，结果显示健康组和卒中组BB的入度和出度在3个频段上都要显著低于其他肌肉(p<0.05)，这更加证实了BB不同于其他肌肉，在肌间耦合网络中耦合程度非常低，发挥作用小。同时发现AD的入度在β和γ频段要显著低于PD(p＝0.048<0.05,p＝0.006<0.05)，AD的出度在γ频段要显著低于PD和TB(p＝0.018<0.05,p＝0.044<0.05)，这也表明了不同于α频段，AD在β和γ频段的肌间耦合网络中耦合程度较低，但仍强于BB。此外，其他肌肉节点之间的入度、出度和聚类系数在各对应频段均不具有显著性差异(p>0.05)，这同样说明了MD、PD和TB在肌间耦合网络中具有同质性。

进一步，利用单因素多元方差分析检验各肌肉节点的入度和出度之间是否具有显著性差异，结果显示在各频段上，每组肌肉节点的入度和出度不具有显著性差异(α:p＝0.711>0.05,β:p＝0.881>0.05,γ:p＝0.654>0.05)。这说明各肌肉在运动控制过程中相互制约，处于动态平衡。

最后，利用Kruskal Wallis检验各频段间网络参数是否具有显著性差异，结果显示BB的入度在α频段要显著高于β频段(p＝0.019<0.05)和γ频段(p＝0.005<0.05)，AD的出度在α频段要显著高于β频段(p＝0.022<0.05)和γ频段(p＝0.003<0.05)，除此之外，其他肌肉节点间的入度、出度和聚类系数在频段间不具有显著性差异(p>0.05)。这说明AD在α频段与其他肌肉间的耦合主要体现为信息的输出，BB在α频段接收到了更高的信息输入，结合卒中组在α频段的肌间耦合网络(图5(b))可以看出，AD对BB具有一定的驱动作用，MD、PD和TB在3个频段上耦合情况基本一致。

在对节点的度和聚类系数进行平均后，可以得到网络的平均度和聚类系数，同样利用配对样本t检验进行组间显著性检验，结果如图7(a)和图7(b)所示。由配对样本t检验结果可见，健康组和卒中组的肌间耦合网络的平均度和聚类系数都不具有显著性差异(p>0.05)，健康组和卒中组在α、β和γ频段上肌间耦合网络的平均度分别为0.23±0.15Bit、0.19±0.12Bit和0.20±0.13Bit，聚类系数分别为0.04±0.04、0.04±0.03和0.04±0.03。同时，利用独立样本t检验3个频段间肌间耦合网络的平均度和聚类系数是否具有显著性差异，结果显示3个频段间肌间耦合网络的平均度和聚类系数同样也不具有显著性差异(p>0.05)。这说明健康组和卒中组的肌间耦合网络随机性强、紧凑程度低，存在一定的相似性，在不同频段上、不同频段间没能发现明显区别，这可能与中枢神经系统的运动控制策略有关。

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

Claims (2)

1.基于高斯Copula传递熵的肌间耦合网络分析方法，其特征在于：该方法包括以下主要步骤：

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

具体为：在表面肌电设备的监控下，同步采集上肢前三角肌、内侧三角肌、后三角肌、肱二头肌、肱三头肌上的5通道表面肌电信号，采样频率为2000Hz；

步骤(2)，高斯Copula传递熵的定义与计算；

具体为：假设两通道表面肌电信号X＝{x₁,x₂,...,x_T}，Y＝{y₁,y₂,...,y_T}分别为一阶马尔科夫过程，T为时间序列的长度，X→Y的高斯Copula传递熵定义为

其中，τ和t分别为X和Y的滞后时间，

k和l分别为X和Y的嵌入维度，ρ₁，ρ₂和ρ₃分别为高斯Copula函数参数，H_c(*)表示变量间的高斯Copula熵；

由高斯Copula传递熵的表达式可知，其关键在于3项高斯Copula熵的估计，高斯Copula熵的估计采用以下3个步骤：

(a)非参数核密度估计边际分布函数

假设各通道表面肌电信号

是来自连续分布函数F_i(x_i)的同分布样本，那么F_i(x_i)的非参数核密度估计为

其中，

为概率密度函数，

其中，h为窗宽，当T→∞，h→0，Th→∞时，非参数核密度估计是真实概率密度分布的一致估计，K(·)为核函数；

(b)典型极大似然法估计高斯Copula参数

在利用非参数核密度估计边际分布函数后，可以无需估计边际分布中的参数，只要估计Copula函数中的参数ρ

(c)蒙特卡洛法估计高斯Copula熵

H_c(x₁,x₂,...,x_d)＝-E[logc(u₁,u₂,...,u_d)] (5)

其中，E[·]表示数学期望；

步骤(3)，肌间耦合网络的搭建与分析；

具体为：

(a)定义网络的节点和边

肌间耦合网络本质上是有向加权复杂网络，可用G表示，G＝(V,E)；V＝{v₁,v₂,...,v_n}为肌肉节点集合，

为边集合；肌间耦合网络的节点数目为n＝|V|，边数为m＝|E|；v_i∈V表示网络中的一个节点，v_i→v_j∈E表示节点v_i到节点v_j的一条有向边；

表示有向边v_i→v_j的权值；

(b)确定网络连接阈值

设置阈值在0.1～0.9*max(GCTE)变化，其中max(GCTE)表示网络中GCTE的最大值，将肌间耦合网络的拉普拉斯矩阵的第二小非零特征值λ_min所对应的阈值确定为最佳阈值；

设第k+1个阈值对应的邻接矩阵为GCTE^k+1，k∈[1,9]，肌间耦合网络的拉普拉斯矩阵L的第二小特征值λ_min等于0

L＝D-GCTE^k+1 (6)

其中，D为肌间耦合网络的度矩阵，对角元素D_ii为GCTE^k+1的第i列元素之和，其余元素为零；则将取第k个阈值为最佳阈值，阈值化后的邻接矩阵为

(c)分析静态网络参数

利用图论中度和聚类系数来分析各频段上肌间耦合网络的拓扑属性；具体分为两个方面：

一是面向节点的节点度和节点聚类系数，反映网络的局部特性；

二是面向网络的网络平均度和聚类系数，反映网络的全局特征；

节点度和节点聚类系数

依据肌间耦合网络，定义每个肌肉节点v_i的入度和出度分别为D_in(v_i)和D_out(v_i)，聚类系数为CC(v_i)

其中，D_in(v_i)越大，表明该节点在网络中是多个节点共同作用的结果，影响该节点的连接越多；D_out(v_i)越大，表明该节点在网络中所起作用越大；当节点之间具有相同的入度和出度时，认为这些节点是同质的；CC(v_i)越大，表明该节点与相邻节点之间存在规则且频繁的转换；

网络的平均度和聚类系数

定义肌间耦合网络的平均度AD为所有节点的入度和出度的平均值，网络的聚类系数NCC为所有节点的聚类系数的平均值

其中，AD越大，网络中连边越多，连通性越好；NCC∈[0,1]，当NCC＝0，表示网络中所有节点都是孤立节点，当NCC＝1，表示网络中任意节点间都有边相连。

2.根据权利要求1所述的基于高斯Copula传递熵的肌间耦合网络分析方法，其特征在于：步骤(1)中的预处理过程为：首先手工提取5次有效活动段数据，然后通过上下采样的方式维持信号时长2.5s，接着进行去均值、去基线漂移，利用谱插值算法抑制50Hz工频干扰，最后采用4阶巴特沃斯带通滤波进行5-200Hz带通滤波，得到纯净的3通道表面肌电信号。

Images