CN112130668B - 一种R藤Copula互信息的肌间耦合分析方法 - Google Patents

一种R藤Copula互信息的肌间耦合分析方法 Download PDF

Info

Publication number
CN112130668B
CN112130668B CN202011031460.2A CN202011031460A CN112130668B CN 112130668 B CN112130668 B CN 112130668B CN 202011031460 A CN202011031460 A CN 202011031460A CN 112130668 B CN112130668 B CN 112130668B
Authority
CN
China
Prior art keywords
copula
rattan
tree
channel
mutual information
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
Application number
CN202011031460.2A
Other languages
English (en)
Other versions
CN112130668A (zh
Inventor
佘青山
王洪安
席旭刚
高云园
孔万增
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Dianzi University
Original Assignee
Hangzhou Dianzi University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hangzhou Dianzi University filed Critical Hangzhou Dianzi University
Priority to CN202011031460.2A priority Critical patent/CN112130668B/zh
Publication of CN112130668A publication Critical patent/CN112130668A/zh
Application granted granted Critical
Publication of CN112130668B publication Critical patent/CN112130668B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F3/00Input arrangements for transferring data to be processed into a form capable of being handled by the computer; Output arrangements for transferring data from processing unit to output unit, e.g. interface arrangements
    • G06F3/01Input arrangements or combined input and output arrangements for interaction between user and computer
    • G06F3/011Arrangements for interaction with the human body, e.g. for user immersion in virtual reality
    • G06F3/015Input arrangements based on nervous system activity detection, e.g. brain waves [EEG] detection, electromyograms [EMG] detection, electrodermal response detection
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/72Signal processing specially adapted for physiological signals or for diagnostic purposes
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2203/00Indexing scheme relating to G06F3/00 - G06F3/048
    • G06F2203/01Indexing scheme relating to G06F3/01
    • G06F2203/011Emotion or mood input determined on the basis of sensed human body parameters such as pulse, heart rate or beat, temperature of skin, facial expressions, iris, voice pitch, brain activity patterns

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Biomedical Technology (AREA)
  • General Health & Medical Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • General Engineering & Computer Science (AREA)
  • Physiology (AREA)
  • Pathology (AREA)
  • General Physics & Mathematics (AREA)
  • Neurosurgery (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Neurology (AREA)
  • Psychiatry (AREA)
  • Signal Processing (AREA)
  • Dermatology (AREA)
  • Biophysics (AREA)
  • Human Computer Interaction (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明公开了一种R藤Copula互信息的肌间耦合分析方法。本发明首先进行多通道表面肌电信号的同步采集与预处理,其次利用非参数核密度估计边际分布函数,再进行R藤Copula的简单矩阵表示及参数估计,同时估计R藤Copula互信息和R藤Copula条件互信息;最后进行肌间耦合分析。本发明提出的RVCMI和RVCCMI为肌间耦合分析提供了一种新的研究方法和科学的理论依据,具有良好的应用前景。

Description

一种R藤Copula互信息的肌间耦合分析方法
技术领域
本发明属于神经系统运动控制机制研究领域,涉及R藤Copula、互信息和条件互信息的计算,从而进行肌间耦合分析。
背景技术
神经肌肉系统是高度复杂的,先前的研究表明,位于额叶中央前回的初级运动皮层接受来自几个大脑区域的输入,帮助规划运动,其主要输出刺激脊髓神经元,刺激骨骼肌收缩。然而,中枢神经系统(Central Nervous System,CNS)如何协调由大量肌肉和关节组成的复杂运动以达到行为目标,仍然是一个基本问题。例如,一个平滑的伸展运动可能是肩膀屈肌、肘部伸展肌和支撑姿势的肌肉协调活动的结果。
在中枢神经系统的功能调节和反馈控制下,肌肉之间的耦合现象反映了人体运动过程中多通道肌肉之间的相互作用和相互联系。表面肌电信号(surfaceElectromyographic,sEMG)是由肌肉收缩产生的生物电活动。由于sEMG信号是记录电极附近多个运动单位动作电位(Motor Unit Action Potentials,MUAP)活动总和的测量值,因此可以实时准确地反映肌间耦合的状态和功能。开展肌间耦合特性的研究,将有助于挖掘潜在的中枢神经系统对肌肉的支配调节方式。
目前,分析两个时间序列和多个时间序列耦合的方法很多,包括相干性、互信息(Mutual Information,MI)、S估计和全局同步指数(Global Synchronization Index,GSI)。相干性方法是在频域上量化两个时间序列之间的线性相关性,但该方法未考虑信号内在的非线性特性。在MI方法中,计算两个时间序列的自概率密度分布和联合概率密度分布,并通过计算各种熵来量化两个时间序列之间的线性和非线性统计独立性。然而,由于非高斯分布数据的联合概率密度函数(Probability Density Function,PDF)的估计是一个难题,使得MI实际应用变得十分困难。对于S估计方法,它是一种基于状态空间的同步方法,即通过分析状态空间重构域中多个信号之间的相互依赖关系来计算同步强度。然而,S估计没有充分考虑随机和伪影分量对分析的影响,计算精度有待提高。对于GSI方法,虽然改进了S估计,以分析多维神经序列,但基于协方差矩阵的GSI方法是一种简单的测量多时间序列线性相关性的方法,不估计多时间序列的非线性相关性,同时也在一定程度上受到噪声的干扰。
最近,Ince等将Copula统计理论与高斯变量熵的封闭解相结合,提出了一种实用的高斯Copula互信息(Gaussian Copula Mutual Information,GCMI)估计方法。GCMI为MI提供了一个计算效率高、统计上稳健的下界估计,无需对每个变量的边际分布进行具体假设。该方法使研究人员能够充分利用每个神经信号的特性及其实验设计,从而更好地了解大脑网络的信息处理功能。然而,单一的高斯Copula函数形式单一,只能刻画对称的相关结构,模型拟合容易失真,不能准确描述所有神经信号之间的功能耦合特性。成对Copula(Pair-Copula)结构是多元分布下的依赖关系的灵活表示,近年来已成为多元分析的一个热门话题。成对Copula的思想是将多元分布分解为条件分布,并通过对两个变量同时建模的双变量Copula来描述这些条件分布。特殊的成对Copula结构称为正则藤(R-Vine)Copula结构,它假定分布中特定元素之间的条件独立性,从而使我们能够避开似然估计和抽样中的维数灾难。与其他传统的Copula结构相比,R藤Copula具有更高的灵活性,能够对更广泛的复杂多元相关性进行建模。
受GCMI方法的启发,本发明提出了一种利用R藤Copula估计信息论量MI和条件MI的新方法。它继承了R藤Copula的优点,直接适用于任意数据大小,并将本发明所提方法应用于执行上肢及物任务的多通道sEMG信号,以揭示肌间耦合的线性和非线性特征。
发明内容
本发明的目的在于提供一种能同时准确分析双通道与多通道、线性与非线性、直接与间接的肌间耦合特性的新方法。
为实现上述目的,本发明方法主要包括以下步骤:
步骤(1),多通道表面肌电信号的同步采集与预处理;
具体为:在表面肌电设备的监控下,同步采集上肢上斜方肌(Upper Trapezius,UT)、前三角肌(Anterior Deltoid,AD)、内侧三角肌(Medial Deltoid,MD)、后三角肌(Posterior Deltoid,PD)、胸大肌(Pectoralis Major,PM)、冈下肌(Infraspinatus,IN)、肱二头肌(Biceps,BB)、肱三头肌(Triceps,TB)上的8通道sEMG信号,采样频率为2000Hz,并对采集到的sEMG信号进行简单的预处理。
步骤(2),非参数核密度估计边际分布函数;
具体为:假设各通道sEMG信号是来自连续分布函数Fi(xi)的同分布样本,T为时间序列的长度,i=1,2,...,8,那么Fi(xi)的非参数核密度估计为
其中,为概率密度函数,
其中,h为窗宽,当T→∞,h→0,Th→∞时,非参数核密度估计是真实概率密度分布的一致估计。K(·)为核函数,本发明采用高斯核函数,根据经验法则确定窗宽: 为xi的标准差。
步骤(3),R藤Copula的简单矩阵表示及参数估计;
具体为:一个典型的R藤由树(tree)、节点(node)和边(edge)组成,每一层树上有若干个节点,每一个节点代表一个变量或条件变量,节点之间的连线称为边,每一条边代表相邻两个节点组成的成双Copula(Pair-Copula)。一个R藤结构图代表一个多维Copula密度函数的分解形式。具体说来,一个N通道sEMG信号的R藤结构由N-1层树T1,T2,...,TN-1组成,第i棵树的节点集记为Ni,边集记为Ei,i=1,2,..,N-1,它们满足以下条件:
(1)树T1的节点集N1={1,2,...,N},边集为E1
(2)第i棵树Ti的节点集Ni=Ei,即第i棵树的节点集是第i-1棵树的边集;
(3)如果树Ti中两条边在树Ti+1中用边连接,那么这两条边在树Ti中必须有一个共同的节点;
下面,建立一个N维R藤统计模型。设N通道sEMG信号x1,x2,...,xN构成的随机向量为X={x1,x2,...,xN},其中第i个变量xi的边缘密度函数为fi,i=1,2,...,N,则X的联合概率密度函数可以表示为
其中,Ei中的边e=j(e),k(e)|D(e),j(e)和k(e)是与边e相连接的两个节点,D(e)是条件集,cj(e),k(e)|D(e)表示边e对应的Pair-Copula密度函数,F(xj(e)|xD(e))和F(xk(e)|xD(e))为由条件集D(e)决定的服从[0,1]均匀分布的转化变量。
本发明在定向图模型的基础上,采用一种基于下三角矩阵的R藤矩阵(R-VineMatrices,RVM)在计算机上较为方便地计算和模拟R藤。由于RVM并非是唯一确定的,对于给定的N维R藤,存在2N-1个不同的R藤,为了确定最恰当的R藤结构模型,本发明采用最大遍历树算法(Maximum Spanning Tree,MST),该算法的关键在于保证R藤结构上每一层树节点的Kendall’sτ绝对值之和最大。在RVM确定之后,本发明采用赤池信息准则(AkaikeInformation Criterion,AIC)从众多的Copula函数集中选择最优的Pair-Copula函数,AIC的计算公式为
AIC=-2ln(L)+2k (4)
其中,k为Pair-Copula函数中参数的个数,L为Pair-Copula函数的极大似然函数值。AIC值越小,R藤模型拟合程度越好。在确定最优的Pair-Copula函数后,进一步利用极大似然估计法估计各Pair-Copula函数中的参数。
步骤(4),估计R藤Copula互信息和R藤Copula条件互信息;
具体为:就连续分布而言,两个随机变量X,Y之间的互信息(Mutual Information,MI)定义为
其中,f(x,y)是X和Y的联合概率密度函数,f(x)和f(y)分别是X和Y的边缘概率密度函数。MI度量了两个随机变量之间线性或非线性的相互依赖程度,反映一个随机变量携带另一个随机变量的信息量大小,是一个非负的关联统计,MI(X,Y)≥0,当且仅当X和Y的相互独立时,MI(X,Y)=0,MI越大,则X和Y共享的信息越多。
条件互信息(Conditional Mutual Information,CMI)定义为
CMI表达的是以第三个随机变量Z为条件的两个随机变量X,Y之间的MI,这里,不妨称Z为条件变量,X和Y为观测变量。CMI能够移除其他变量的影响,度量X和Y的直接联系(包含线性和非线性)。CMI(X,Y|Z)≥0,当且仅当在给定Z的条件下,X和Y相互独立。当条件变量Z个数为一时,称为一阶CMI,当条件变量个数大于1时,称为高阶CMI。当对数运算以2为基底时,MI和CMI的单位是比特(Bit)。
由式(5)~(6)可知,MI与CMI的估计严重依赖变量的边际分布函数和联合分布函数的精确表示,边际分布函数可通过直方图法或核密度估计法比较好解决,然而联合分布函数的估计往往十分困难,现有的主要估计方法,如参数化法需要对数据作高斯假设,非参数核估计法和sieve法收敛速度慢,对维数变化敏感,这使得MI与CMI在分析多变量相依程度时十分棘手。
Ma和Sun等证明了MI本质上是一种Copula熵,将Copula密度函数代入MI的表达式中,可得
其中,Hc(F(x),F(y))称为二维Copula熵。推广至高维情形,N维Copula熵被定义为
其中,ui=Fi(xi)。由上式可知,Copula熵不局限于双变量,也适用于多变量的情况,且具有以下推论:
式(9)表明一个N维Copula熵可以分解为两部分:一个N变量联合熵和N个边界熵。当且仅当所有被选择变量相互独立时,Hc(u1,u2,...,uN)=0。
将上述推论代入CMI的熵版本中,可得
式(10)表明CMI可以改写为3项Copula熵的组合:前两项分别为不同观测变量X、Y分别与条件变量Z的Copula熵,后一项为所有变量的Copula熵,这与CMI的原始形式构成了显著的区别。
以上表明,基于Copula熵的MI和CMI的估计关键在于Copula函数的选择。然而,当变量维数大于2时,传统的高维Copula函数在参数估计上面临着维数灾难、模型精度不高、灵活性差等问题。而R藤Copula通过将高维Copula函数分解成一系列Pair-Copula函数的乘积,能有效地规避上述问题。因此,本发明将R藤Copula密度函数代入MI和CMI的表达式中,即可得到R-Vine Copula MI(RVCMI)和R-Vine Copula CMI(RVCCMI)
其中,x1,x2,...,xm表示m个观测变量,z1,z2,...,zn表示n个条件变量,ui=Fi(xi),vi=Fi(zi)。式(11)~(12)表明,R-Vine Copula MI和R-Vine Copula CMI可在统一的理论框架下进行相依程度度量。
步骤(5),肌间耦合分析;
具体为:为体现R藤Copula描述肌间相依结构的优劣,同时将C藤Copula和D藤Copula纳入比较,C藤和D藤是R藤模型中两类最简单的藤,2种藤具有不同的逻辑结构,C藤为星型结构,D藤为平行结构,适合于不同类型的数据集合。进一步,利用Vuong检验R藤Copula是否在统计上显著优于C藤Copula和D藤Copula,若Vuong统计量为正,则倾向于前1个模型,若p-Value大于0.05,则不能拒绝两个模型没有区别的原假设。
在肌间耦合分析时,可简单地分为双通道和多通道分析:1)利用双变量版本的RVCMI(式(11),m=2)和高阶版本的RVCCMI(式12,m=2,n=6)度量双通道肌间间接和直接的非线性耦合强度关系,即将其他6通道sEMG信号视作条件变量;2)利用多变量版本的RVCMI(式(11),m=8)和一阶版本的RVCCMI(式(12),m=7,n=1)度量多通道肌间间接和直接的非线性耦合强度关系,即将某一通道sEMG信号视作条件变量。
本发明与传统的肌间耦合分析方法相比,具有如下优点:
传统的肌间耦合分析方法侧重两两通道之间线性的功能耦合分析,缺乏对多通道信号之间直接的非线性耦合作用的有效描述,针对该问题,本发明提出一种新R藤Copula互信息和R藤Copula条件互信息估计方法,用于准确定量地描述上肢运动过程中多通道肌肉之间的线性和非线性耦合特性,这为探索运动过程肌肉系统耦合特征及中枢神经系统的运动控制机制提供新方法。
附图说明
图1为本发明方法的流程图。
图2为肌肉记录位置。
图3为受试者H1预处理后8通道sEMG信号。
图4(a)为8通道sEMG信号高斯性检验的Q-Q图,图4(b)为边际累积分布函数估计曲线。
图5(a)、图5(b)分别为基于RVCMI和RVCCMI的肌间耦合强度。
图6为多通道肌间耦合作用成分占比。
具体实施方式
下面结合附图对本发明的实施例作详细说明:本实施例在以本发明技术方案为前提下进行实施,给出了详细的实施方案和具体的操作过程。
肌间耦合起源于皮质脊髓通路,可以反映来自脊髓中间神经元的共同驱动,体现为运动过程中肌肉间的相互作用以及中枢神经对肌肉的支配能力。肌间耦合特性的研究对于人体运动理论分析、脑卒中运动功能康复评估以及中枢神经系统运动控制机制的探索等方面具有十分重要的意义。然而,不同的耦合分析方法因其自身的特点不同,将导致所揭示的肌间耦合关系有所差异。
信息论中的MI和CMI是统计学中相关系数和偏相关系数的非线性扩展,不仅可以区分直接联系和间接联系,而且对变量的分布类型没有特殊要求,在很多方面被认为是度量随机变量之间关联程度的完美统计。然而,关于MI和CMI的估计一直是一个具有数值挑战性的问题。本发明在GCMI和GCCMI的基础上,将R藤Copula理论和熵理论相结合提出了RVCMI和RVCCMI方法,不但继承了R藤Copula和MI、CMI的优点,而且还考虑了变量间相依结构。
如图1所示,本发明的实施主要包括六个步骤:(1)同步采集多通道表面肌电信号;(2)预处理;(3)非参数核密度估计边际分布函数;(4)R藤Copula的简单矩阵表示及参数估计;(5)估计R藤Copula互信息和R藤Copula条件互信息;(6)肌间耦合分析。
下面逐一对各步骤进行详细说明。
步骤一:同步采集多通道表面肌电信号
采用表面肌电设备同步采集及物运动过程中上肢上斜方肌(Upper Trapezius,UT)、前三角肌(Anterior Deltoid,AD)、内侧三角肌(Medial Deltoid,MD)、后三角肌(Posterior Deltoid,PD)、胸大肌(Pectoralis Major,PM)、冈下肌(Infraspinatus,IN)、肱二头肌(Biceps,BB)、肱三头肌(Triceps,TB)上的8通道sEMG信号,采样频率为2000Hz,肌肉记录位置如图2所示。
步骤二:预处理
预处理过程,首先手工提取5次有效激活的活动段数据,并通过上下采样的方式维持信号时长为2.5s,然后进行去均值、去基线漂移处理,同时利用二阶IIR陷波滤波器抑制50Hz工频干扰,最后采用四阶巴特沃斯带通滤波进行带通滤波,截止频率分别为5Hz和200Hz,为方便后续分析,将预处理后得到纯净的sEMG信号按实验次序平均。图3给出了受试者H1的结果,左列为各通道sEMG信号的时域波形,右列为对应的最大值归一化后的功率谱。
步骤三:非参数核密度估计边际分布函数
假设各通道sEMG信号是来自连续分布函数Fi(xi)的同分布样本,T为时间序列的长度,i=1,2,...,8,那么Fi(xi)的非参数核密度估计为
其中,为概率密度函数,
其中,h为窗宽,当T→∞,h→0,Th→∞时,非参数核密度估计是真实概率密度分布的一致估计。K(·)为核函数,本发明采用高斯核函数,根据经验法则确定窗宽: 为xi的标准差。
步骤四:R藤Copula的简单矩阵表示及参数估计
一个典型的R藤由树(tree)、节点(node)和边(edge)组成,每一层树上有若干个节点,每一个节点代表一个变量或条件变量,节点之间的连线称为边,每一条边代表相邻两个节点组成的成双Copula(Pair-Copula)。一个R藤结构图代表一个多维Copula密度函数的分解形式。具体说来,一个N通道sEMG信号的R-Vine结构由N-1层树T1,T2,...,TN-1组成,第i棵树的节点集记为Ni,边集记为Ei,i=1,2,..,N-1,它们满足以下条件:
(1)树T1的节点集N1={1,2,...,N},边集为E1
(2)第i棵树Ti的节点集Ni=Ei,i=1,2,..,N-1,即第i棵树的节点集是第i-1棵树的边集;
(3)如果树Ti中两条边在树Ti+1中用边连接,那么这两条边在树Ti中必须有一个共同的节点。
下面,建立一个N维R藤统计模型。设N通道sEMG信号x1,x2,...,xN构成的随机向量为X={x1,x2,...,xN},其中第i个变量xi的边缘密度函数为fi,i=1,2,...,N,则X的联合概率密度函数可以表示为
其中,Ei中的边e=j(e),k(e)|D(e),j(e)和k(e)是与边e相连接的两个节点,D(e)是条件集,cj(e),k(e)|D(e)表示边e对应的Pair-Copula密度函数,F(xj(e)|xD(e))和F(xk(e)|xD(e))为由条件集D(e)决定的服从[0,1]均匀分布的转化变量。
本发明在定向图模型的基础上,采用一种基于下三角矩阵的R藤矩阵(R-VineMatrices,RVM)在计算机上较为方便地计算和模拟R藤。由于RVM并非是唯一确定的,对于给定的N维R藤,存在2N-1个不同的R藤,为了确定最恰当的R藤结构模型,本发明采用最大遍历树算法(Maximum Spanning Tree,MST),该算法的关键在于保证R藤结构上每一层树节点的Kendall’sτ绝对值之和最大。在RVM确定之后,本发明采用赤池信息准则(AkaikeInformation Criterion,AIC)从众多的Copula函数集中选择最优的Pair-Copula函数,AIC的计算公式为
AIC=-2ln(L)+2k (4)
其中,k为Pair-Copula函数中参数的个数,L为Pair-Copula函数的极大似然函数值。AIC值越小,R藤模型拟合程度越好。在确定最优的Pair-Copula函数后,进一步利用极大似然估计法估计各Pair-Copula函数中的参数。
步骤五:估计R藤Copula互信息和R藤Copula条件互信息
就连续分布而言,两个随机变量X,Y之间的互信息(Mutual Information,MI)定义为
其中,f(x,y)是X和Y的联合概率密度函数,f(x)和f(y)分别是X和Y的边缘概率密度函数。MI度量了两个随机变量之间线性或非线性的相互依赖程度,反映一个随机变量携带另一个随机变量的信息量大小,是一个非负的关联统计,MI(X,Y)≥0,当且仅当X和Y的相互独立时,MI(X,Y)=0,MI越大,则X和Y共享的信息越多。
条件互信息(Conditional Mutual Information,CMI)定义为
CMI表达的是以第三个随机变量Z为条件的两个随机变量X,Y之间的MI,这里,不妨称Z为条件变量,X和Y为观测变量。CMI能够移除其他变量的影响,度量X和Y的直接联系(包含线性和非线性)。CMI(X,Y|Z)≥0,当且仅当在给定Z的条件下,X和Y相互独立。当条件变量Z个数为一时,称为一阶CMI,当条件变量个数大于1时,称为高阶CMI。当对数运算以2为基底时,MI和CMI的单位是比特(Bit)。
由式(5)~(6)可知,MI与CMI的估计严重依赖变量的边际分布函数和联合分布函数的精确表示,边际分布函数可通过直方图法或核密度估计法比较好解决,然而联合分布函数的估计往往十分困难,现有的主要估计方法,如参数化法需要对数据作高斯假设,非参数核估计法和sieve法收敛速度慢,对维数变化敏感,这使得MI与CMI在分析多变量相依程度时十分棘手。
Ma和Sun等证明了MI本质上是一种Copula熵,将Copula密度函数代入MI的表达式中,可得
其中,Hc(F(x),F(y))称为二维Copula熵,推广至高维情形,N维Copula熵被定义为
其中,ui=Fi(xi)。由上式可知,Copula熵不局限于双变量,也适用于多变量的情况,且具有以下推论:
式(9)表明一个N维Copula熵可以分解为两部分:一个N变量联合熵和N个边界熵。当且仅当所有被选择变量相互独立时,Hc(u1,u2,...,uN)=0。
将上述推论代入CMI的熵版本中,可得
式(10)表明CMI可以改写为3项Copula熵的组合:前两项分别为不同观测变量X、Y分别与条件变量Z的Copula熵,后一项为所有变量的Copula熵,这与CMI的原始形式构成了显著的区别。
以上表明,基于Copula熵的MI和CMI的估计关键在于Copula函数的选择。然而,当变量维数大于2时,传统的高维Copula函数在参数估计上面临着维数灾难、模型精度不高、灵活性差等问题。而R藤Copula通过将高维Copula函数分解成一系列Pair-Copula函数的乘积,能有效地规避上述问题。因此,本发明将R藤Copula密度函数代入MI和CMI的表达式中,即可得到R-Vine Copula MI(RVCMI)和R-Vine Copula CMI(RVCCMI)
其中,x1,x2,...,xm表示m个观测变量,z1,z2,...,zn表示n个条件变量,ui=Fi(xi),vi=Fi(zi)。式(11)~(12)表明,R-Vine Copula MI和R-Vine Copula CMI可在统一的理论框架下进行相依程度度量。
步骤六:肌间耦合分析
为体现R藤Copula描述肌间相依结构的优劣,同时将C藤Copula和D藤Copula纳入比较,C藤和D藤是R藤模型中两类最简单的藤,2种藤具有不同的逻辑结构,C藤为星型结构,D藤为平行结构,适合于不同类型的数据集合。进一步,利用Vuong检验R藤Copula是否在统计上显著优于C藤Copula和D藤Copula,若Vuong统计量为正,则倾向于前1个模型,若p-Value大于0.05,则不能拒绝两个模型没有区别的原假设。
在肌间耦合分析时,可简单地分为双通道和多通道分析:1)利用双变量版本的RVCMI(式(11),m=2)和高阶版本的RVCCMI(式12,m=2,n=6)度量双通道肌间间接和直接的非线性耦合强度关系,即将其他6通道sEMG信号视作条件变量;2)利用多变量版本的RVCMI(式(11),m=8)和一阶版本的RVCCMI(式(12),m=7,n=1)度量多通道肌间间接和直接的非线性耦合强度关系,即将某一通道sEMG信号视作条件变量。
为了验证本发明方法的性能,实验部分采集了多个受试者的数据,通过性能指标来衡量本发明在肌间耦合分析领域的可行性。实验招募了5名健康受试者参与(H1~H5),平均年龄小于76岁。实验要求每位受试者坐在桌子前,前臂放在一个舒服的位置上,根据每隔10秒激活的语音提示,向着正前方目标进行5次及物动作,目标距桌子高35cm,每次及物后休息10s,受试者在执行运动任务前已知晓实验内容。
为更好理解sEMG信号的分布特点,在估计边际累计分布函数之前,本发明先利用Q-Q图检验了8通道sEMG信号的高斯性,结果如图4(a)所示。由图4(a)可以看出,8通道sEMG信号的分位数与标准高斯分位数不在一条直线上,相距较远,并且在两个对角线外翘,尾部较厚,表明sEMG信号并不遵循高斯分布。边际累计分布函数估计结果如图4(b)所示。由图4(b)可以看出,8通道sEMG信号估计的累计分布曲线单调递增,十分平滑,大部分幅值(-0.5×10-4~0.5×10-4)集中都在零均值两侧,左右分布对称,进一步计算8通道sEMG信号的峰度,结果发现峰度的平均值为5.6780>3,标准差为0.8502,说明sEMG信号存在尖峰厚尾现象。
在估计R藤Copula中的参数以及利用R藤Copula刻画肌间相依结构之前,要求输入的sEMG信号的边际累计分布函数服从(0,1)均匀分布,所以本发明利用单样本Kolmogorov-Smirnov(KS)检验法来检验拟合好的边际累积分布函数,显著性水平取为0.05,KS统计量越小,均匀分布的可能性越大。以图4(b)为例,KS检验结果如表1所示。
表1累积分布函数的KS检验结果(H1)
由表1可见,KS统计量普遍较小,p-Value均大于0.05,表明在0.05的显著性水平下8通道sEMG信号的MCDF不能拒绝(0,1)均匀分布的原假设。在此基础上,本发明根据步骤3所述,借助R软件中的“vinecopula”包构建基于R藤Copula的肌间相依结构,由最大遍历树生成算法得到的R藤矩阵。为体现R藤Copula描述肌间相依结构的优劣,同时将C藤Copula和D藤Copula纳入比较,3种藤Copula模型的分解出的Pair-Copula结构和参数估计结果如表2所示,表2中编号1,2,3,4,5,6,7,8分别表示UT,AD,MD,PD,PM,IN,BB,TB,“|”之后的变量表示条件变量。
表2 C藤Copula的参数估计以及相依结构(H1)
/>
表3 D藤Copula的参数估计以及相依结构(H1)
/>
表4 R藤Copula的参数估计以及相依结构(H1)
/>
注:表2~表4中树的边上的编号1,2,3,4,5,6,7,8分别表示UT,AD,MD,PD,PM,IN,BB,TB。Copula即为Pair-Copula:RC(90)为Rotated Clayton Copula 90degrees,SC为Survival Clayton Copula,SBB8为Survival BB8 Copula,RG(90)为Rotated GumbelCopula 90degrees,SBB8(90)为Rotated BB8 Copula 90degrees,RBB8(90)为Rotated BB8Copula 90degrees,RBB8(270)为Rotated BB8 Copula 270degrees,RT2(90)为RotatedTawn type 2Copula 90degrees,RT2(180)为Rotated Tawn type 2Copula 180degrees,BB7为BB7 Copula,RT1(90)为Rotated Tawn type 1Copula 90degrees,T1为Tawn type1Copula,RC(270)为Rotated Clayton Copula 270degrees。Par1表示Copula中的第1个参数,Par2表示Copula中的第2个参数,如t Copula的Par1表示线性相关参数,Par2表示自由度参数。τ为Kendall’sτ秩相关参数。详细内容请参见R-packages:“vinecopula”和“CDVine”。
由表2~表4可见,在藤Copula建立的肌间相依结构中,T1考虑的是边与边之间的无条件相关性,T2~T6考虑的是1个到多个条件变量下肌间的条件相关性。3种藤Copula在T1~T6的不同边上,根据AIC准则挑选出的最优Pair-Copula既有相同之处,又有所区别。相同点:表中大部分边之间都选择t Copula作为Pair-Copula,在个别边UT和PD(1,4),MD和BB(3,7)上,3种藤Copula分别选择的是SBB8 Copula和RG(90)Copula。不同点:少部分边如PD和PM(4,5),C藤Copula和D藤Copula选择的是Frank,而R藤Copula选择的是RT2(90),再如边IN和BB(6,7),C藤Copula选择的是SBB8(90),D藤Copula和R藤Copula选择的是RBB8(90)。Kendall’sτ秩相关参数常用于衡量节点之间一致性变化程度,不随单调变换的影响,在Copula理论中要优于Pearson线性相关系数。如果节点之间的和谐概率大于不和谐概率,则Kendall’sτ为正,反之,为负。3种藤Copula给出的Kendall’sτ秩相关参数都明显较小,接近0,表明肌间相关性很弱。但仍在一些节点间发现较强的负相关性,AD和TB(2,8)的Kendall’sτ秩相关参数在C藤Copula中为-0.18,在D藤Copula中的为-0.13,在R藤Copula中为-0.17;MD和TB(3,8)的Kendall’sτ秩相关参数在C藤Copula中为-0.24,在D藤Copula中为-0.21,在R藤Copula中为-0.25;PD和TB(4,8)的Kendall’sτ秩相关参数在C藤Copula中为-0.12,在D藤Copula中为-0.18,在R藤Copula中为-0.19。表明TB与三角肌之间存在一定程度上的联系,尤其是TB与MD,在3种藤Copula中的Kendall’sτ秩相关参数的绝对值都是最大的。此外,R藤Copula拥有的AIC值要比C藤Copula和D藤Copula都要低,R藤Copula拥有的对数似然函数值比C藤Copula和D藤Copula都要高,说明作为更一般的R藤Copula对肌间相依结构的拟合程度更高,模型刻画更加准确。
进一步,利用Vuong检验R藤Copula是否在统计上显著优于C藤Copula和D藤Copula,若Vuong统计量为正,则倾向于前1个模型,若p-Value大于0.05,则不能拒绝两个模型没有区别的原假设,检验结果如表3所示。由表3可见,3次Vuong检验都倾向于前者,即在0.05的显著性水平下,R藤Copula建立的肌间相依结构要显著优于C藤Copula和D藤Copula(p-Value<0.05),虽然C藤Copula较之D藤Copula的Vuong统计量为正,但不具有显著性差异(p-Value>0.05)。
表3 3种藤Copula模型Vuong检验结果
注:“*”:p-Value<0.05。
位于图5(a)图5(b)主对角线上的RVCMI值和RVCCMI值即为自信息(SelfInformation,SI),SI衡量了某一通道sEMG信号自身所包含的信息量的多少。由图5(a)或图5(b)可见,8通道sEMG信号的自信息非常相近,在7.36Bit左右。由于高阶的RVCCMI(图6)考虑的是直接的耦合作用,所以自然会低于无条件项的RVCMI(图6)。从图6可以看出,TB与MD的RVCMI值(0.1774Bit)和RVCCMI值(0.1191Bit)都明显要高于其它肌肉之间RVCMI值和RVCCMI值,这说明MD与TB存在强烈的非线性耦合。另外,TB与AD、PD、IN也存在较高的耦合强度,这说明TB与三角肌之间耦合密切。这些与表2中的Kendall’s tau秩相关参数得出的结论一致。除此之外,我们还可以发现,PM和BB与其它肌肉的RVCMI值和RVCCMI值都特别低(<0.01Bit),这说明PM和BB与其它肌肉基本上不存在耦合作用,可能是相互独立的。进一步,利用Mantel test检验H1~H5之间的双通道肌间耦合强度指标RVCMI、RVCCMI的Spearman相关性,显著性水平取为0.05,检验结果如表4所示。
表4受试者之间关于双通道肌间耦合强度的Mantel test相关性检验
注:表中r1,r2分别表示关于RVCMI和RVCCMI的Mantel统计量,p1,p2为相应的p-Value,“*”:p-Value<0.05,“**”:p-Value<0.01。
由表4可知,H1~H5之间的双通道肌间耦合强度指标RVCMI、RVCCMI显著相关(r1>0.62,p1<0.05;r2>0.51,p2<0.05),这说明5名受试者的双通道肌间耦合强度关系是较为一致的。
为体现被移除的某一通道肌肉在整个肌肉系统耦合作用中的相对重要程度,本文进一步计算了一阶RVCCMI与多变量RVCMI的比值,结果以百分数的形式给出,如图6所示。图6中黑色扇形区域表示剩余了多少耦合作用成分,白色色扇形区域则表示了移除了多少耦合作用成分。
由图6可见,当移除MD或TB时,8通道肌间耦合作用强度分别减少了60.13%和69.44%,剩余的7通道肌间耦合作用还不到原来的一半,而在移除PM或BB时,8通道肌间耦合作用强度仅分别减少了4.41%和5.38%,剩余的7通道肌间耦合作用还比较强,并没有受到很大程度的影响(>90%),这说明MD和TB在多通道肌间耦合作用中起着主导作用,而PM与BB以消极的方式参与多通道肌间耦合作用,贡献甚微,这与双通道肌间耦合分析时得出的结论是一致的。进一步,利用Pearson相关系数检验H1~H5之间的多通道肌间耦合强度指标RVCCMI的相关性,显著性水平取为0.05,检验结果如表5所示。
表5受试者之间关于多通道肌间耦合强度的Pearson相关性检验
注:表中r表示关于RVCCMI的Pearson相关系数,p为相应的p-Value,“*”:p-Value<0.05,“**”:p-Value<0.01。
从表5可以看出,H1~H5之间的多通道肌间耦合强度指标RVCCMI几乎是显著相关的(r>0.70,p<0.05),其中,H5与H2、H4之间的Pearson相关系数高达0.9以上,这说明5名受试者的多通道肌间耦合强度关系是非常接近的。
以上所述的实施例仅仅是对本发明的优选实例方式进行描述,并非对本发明的范围进行限定,在不脱离本发明设计精神的前提下,本领域普通技术人员对本发明的技术方案做出的各种变形和改进,均应落入本发明权力要求书确定的保护范围。

Claims (3)

1.一种R藤Copula互信息的肌间耦合分析方法,其特征在于:该方法包括以下主要步骤:
步骤(1),多通道表面肌电信号的同步采集与预处理;
具体为:在表面肌电设备的监控下,同步采集上肢上斜方肌、前三角肌、内侧三角肌、后三角肌、胸大肌、冈下肌、肱二头肌、肱三头肌上的N通道sEMG信号,采样频率为2000Hz,并对采集到的sEMG信号进行预处理;
步骤(2),非参数核密度估计边际分布函数;
具体为:假设各通道sEMG信号是来自连续分布函数Fi(xi)的同分布样本,T为时间序列的长度,i=1,2,...,N,那么Fi(xi)的非参数核密度估计为
其中,为概率密度函数;
步骤(3),R藤Copula的简单矩阵表示及参数估计;
具体为:一个N通道sEMG信号的R藤结构由N-1层树T1,T2,...,TN-1组成,第i棵树的节点集记为Ni,边集记为Ei,i=1,2,..,N-1,它们满足以下条件:
树T1的节点集N1={1,2,...,N},边集为E1
第i棵树Ti的节点集Ni=Ei,即第i棵树的节点集是第i-1棵树的边集;
如果树Ti中两条边在树Ti+1中用边连接,那么这两条边在树Ti中必须有一个共同的节点;
下面,建立一个N维R藤统计模型:设N通道sEMG信号x1,x2,...,x8构成的随机向量为X={x1,x2,...,x8},其中第i个变量xi的边缘密度函数为fi,则随机向量X的联合概率密度函数表示为
其中,Ei中的边e=j(e),k(e)|D(e),j(e)和k(e)是与边e相连接的两个节点,D(e)是条件集,cj(e),k(e)|D(e)表示边e对应的Pair-Copula密度函数,F(xj(e)|xD(e))和F(xk(e)|xD(e))为由条件集D(e)决定的服从[0,1]均匀分布的转化变量;
在定向图模型的基础上,采用基于下三角矩阵的R藤矩阵RVM计算和模拟R藤;在R藤矩阵RVM确定之后,采用赤池信息准则从众多的Copula函数集中选择最优的Pair-Copula函数,在确定最优的Pair-Copula函数后,利用极大似然估计法估计各Pair-Copula函数中的参数;
步骤(4),估计R藤Copula互信息和R藤Copula条件互信息;
具体为:就连续分布而言,两个随机变量X,Y之间的互信息MI定义为
其中,f(x,y)是X和Y的联合概率密度函数,f(x)和f(y)分别是X和Y的边缘概率密度函数;
条件互信息CMI定义为
CMI表达的是以第三个随机变量Z为条件的两个随机变量X,Y之间的MI;
将R藤Copula密度函数代入MI和CMI的表达式中,即可得到RVCMI和RVCCMI
其中,x1,x2,...,xm表示m个观测变量,z1,z2,...,zn表示n个条件变量,ui=Fi(xi),vi=Fi(zi);
步骤(5),肌间耦合分析;
在肌间耦合分析时,分为双通道和多通道分析:
1)令m=2,计算RVCMI;同时令n=6,计算RVCCMI,从而度量双通道肌间间接和直接的非线性耦合强度关系,即将其它6通道sEMG信号视作条件变量;
2)令m=2,计算RVCMI;再令m=7,n=1,计算RVCCMI,从而度量多通道肌间间接和直接的非线性耦合强度关系,即将某一通道sEMG信号视作条件变量。
2.根据权利要求1所述的一种R藤Copula互信息的肌间耦合分析方法,其特征在于:步骤(1)中的预处理具体是:
首先手工提取5次有效激活的活动段数据,并通过上下采样的方式维持信号时长为2.5s;
然后进行去均值、去基线漂移处理,同时利用二阶IIR陷波滤波器抑制50Hz工频干扰;
最后采用四阶巴特沃斯带通滤波进行带通滤波,截止频率分别为5Hz和200Hz。
3.根据权利要求1所述的一种R藤Copula互信息的肌间耦合分析方法,其特征在于:步骤(2)中由于R藤矩阵RVM并非是唯一确定的,对于给定的N维R藤,存在2N-1个不同的R藤,为了确定最恰当的R藤模型,采用最大遍历树算法,该算法的关键在于保证R藤结构上每一层树节点的Kendall’sτ绝对值之和最大。
CN202011031460.2A 2020-09-27 2020-09-27 一种R藤Copula互信息的肌间耦合分析方法 Active CN112130668B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011031460.2A CN112130668B (zh) 2020-09-27 2020-09-27 一种R藤Copula互信息的肌间耦合分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011031460.2A CN112130668B (zh) 2020-09-27 2020-09-27 一种R藤Copula互信息的肌间耦合分析方法

Publications (2)

Publication Number Publication Date
CN112130668A CN112130668A (zh) 2020-12-25
CN112130668B true CN112130668B (zh) 2024-02-02

Family

ID=73840216

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011031460.2A Active CN112130668B (zh) 2020-09-27 2020-09-27 一种R藤Copula互信息的肌间耦合分析方法

Country Status (1)

Country Link
CN (1) CN112130668B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112509689B (zh) * 2021-02-08 2024-05-17 杭州电子科技大学 基于时变Copula互信息的肌间耦合分析方法
CN113440150A (zh) * 2021-06-25 2021-09-28 睿旭康(苏州)智能技术有限责任公司 基于R-vine Copula的皮层肌肉功能网络构建方法
CN113558638A (zh) * 2021-07-05 2021-10-29 杭州电子科技大学 基于Vine Copula的脑肌耦合模型构建方法
CN114298126B (zh) * 2021-10-26 2024-05-31 北京工业大学 一种基于条件互信息与核密度估计的脑功能网络分类方法
CN114137832B (zh) * 2021-10-26 2024-03-29 杭州电子科技大学 基于R藤Copula传递熵的多变量因果关系方法
CN114159081A (zh) * 2021-12-13 2022-03-11 杭州电子科技大学 基于时变混合Copula互信息的肌电耦合方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109165455A (zh) * 2018-09-04 2019-01-08 南京大学 基于互信息和vine copula的水文相依结构建模方法
CN109497999A (zh) * 2018-12-20 2019-03-22 杭州电子科技大学 基于Copula-GC的脑肌电信号时频耦合分析方法
CN110032811A (zh) * 2019-04-17 2019-07-19 电子科技大学 基于Copula函数的工业机器人电气驱动器的可靠性分析方法
CN111258428A (zh) * 2020-01-20 2020-06-09 西安臻泰智能科技有限公司 脑电控制系统及方法
CN111708978A (zh) * 2020-07-23 2020-09-25 杭州电子科技大学 多尺度时频肌间耦合分析方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8311623B2 (en) * 2006-04-15 2012-11-13 The Board Of Trustees Of The Leland Stanford Junior University Systems and methods for estimating surface electromyography
US11406316B2 (en) * 2018-02-14 2022-08-09 Cerenion Oy Apparatus and method for electroencephalographic measurement

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109165455A (zh) * 2018-09-04 2019-01-08 南京大学 基于互信息和vine copula的水文相依结构建模方法
CN109497999A (zh) * 2018-12-20 2019-03-22 杭州电子科技大学 基于Copula-GC的脑肌电信号时频耦合分析方法
CN110032811A (zh) * 2019-04-17 2019-07-19 电子科技大学 基于Copula函数的工业机器人电气驱动器的可靠性分析方法
CN111258428A (zh) * 2020-01-20 2020-06-09 西安臻泰智能科技有限公司 脑电控制系统及方法
CN111708978A (zh) * 2020-07-23 2020-09-25 杭州电子科技大学 多尺度时频肌间耦合分析方法

Also Published As

Publication number Publication date
CN112130668A (zh) 2020-12-25

Similar Documents

Publication Publication Date Title
CN112130668B (zh) 一种R藤Copula互信息的肌间耦合分析方法
CN109262618B (zh) 基于肌肉协同的上肢多关节同步比例肌电控制方法与系统
US20040015894A1 (en) Method and apparatus for extracting low SNR transient signals from noise
Abbaspour et al. A novel approach for removing ECG interferences from surface EMG signals using a combined ANFIS and wavelet
Wang et al. ECG signal denoising based on deep factor analysis
CN111931606B (zh) 基于混合Copula互信息的肌间耦合分析方法
CN109674445B (zh) 一种结合非负矩阵分解和复杂网络的肌间耦合分析方法
CN112232301A (zh) 基于多尺度Copula互信息的肌间耦合网络分析方法
CN113229831B (zh) 基于肌电和肌氧信号的运动功能监测管理方法
Abbaspour et al. Evaluation of wavelet based methods in removing motion artifact from ECG signal
Li et al. ECG denoising method based on an improved VMD algorithm
CN114469124A (zh) 一种运动过程中异常心电信号的识别方法
Liu et al. Multiscale transfer spectral entropy for quantifying corticomuscular interaction
CN112509689A (zh) 基于时变Copula互信息的肌间耦合分析方法
Narmada et al. A novel adaptive artifacts wavelet Denoising for EEG artifacts removal using deep learning with Meta-heuristic approach
Dey et al. Deep learning algorithms for efficient analysis of ecg signals to detect heart disorders
Priyadharsini et al. An Efficient method for the removal of ECG artifact from measured EEG Signal using PSO algorithm
Veer A flexible approach for segregating physiological signals
Jiang et al. Heartbeat classification system based on modified stacked denoising autoencoders and neural networks
CN113967025B (zh) 一种基于肌电信号的多层次运动功能评估方法
Vena et al. Detection of physiological singularities in respiratory dynamics analyzed by recurrence quantification analysis of tracheal sounds
Murthy et al. Design and implementation of hybrid techniques and DA-based reconfigurable FIR filter design for noise removal in EEG signals on FPGA
CN115270847A (zh) 基于小波包分解和卷积神经网络的设计决策脑电识别方法
CN113440150A (zh) 基于R-vine Copula的皮层肌肉功能网络构建方法
N Elbedwehy et al. ECG Denoising using a Single-Node Dynamic Reservoir Computing Architecture.

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