CN111931606B - 基于混合Copula互信息的肌间耦合分析方法 - Google Patents
基于混合Copula互信息的肌间耦合分析方法 Download PDFInfo
- Publication number
- CN111931606B CN111931606B CN202010715642.5A CN202010715642A CN111931606B CN 111931606 B CN111931606 B CN 111931606B CN 202010715642 A CN202010715642 A CN 202010715642A CN 111931606 B CN111931606 B CN 111931606B
- Authority
- CN
- China
- Prior art keywords
- copula
- inter
- mutual information
- function
- mixed
- 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
Links
- 241000039077 Copula Species 0.000 title claims abstract description 125
- 230000008878 coupling Effects 0.000 title claims abstract description 48
- 238000010168 coupling process Methods 0.000 title claims abstract description 48
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 48
- 238000004458 analytical method Methods 0.000 title claims abstract description 18
- 230000006870 function Effects 0.000 claims abstract description 83
- 210000001364 upper extremity Anatomy 0.000 claims abstract description 10
- 238000000034 method Methods 0.000 claims description 44
- 238000005315 distribution function Methods 0.000 claims description 11
- 238000012360 testing method Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 7
- 238000005070 sampling Methods 0.000 claims description 6
- 230000001419 dependent effect Effects 0.000 claims description 5
- 238000001914 filtration Methods 0.000 claims description 4
- 230000001360 synchronised effect Effects 0.000 claims description 4
- 241000932075 Priacanthus hamrur Species 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 238000012544 monitoring process Methods 0.000 claims description 2
- 238000001228 spectrum Methods 0.000 claims description 2
- 230000033001 locomotion Effects 0.000 abstract description 12
- 208000006011 Stroke Diseases 0.000 description 16
- 210000003205 muscle Anatomy 0.000 description 10
- 230000003183 myoelectrical effect Effects 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 206010008190 Cerebrovascular accident Diseases 0.000 description 3
- 208000012661 Dyskinesia Diseases 0.000 description 3
- 238000000540 analysis of variance Methods 0.000 description 3
- 210000003169 central nervous system Anatomy 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 238000011156 evaluation Methods 0.000 description 3
- 230000003993 interaction Effects 0.000 description 3
- 230000007246 mechanism Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000004913 activation Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000002490 cerebral effect Effects 0.000 description 2
- 208000026106 cerebrovascular disease Diseases 0.000 description 2
- 230000001186 cumulative effect Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 208000035475 disorder Diseases 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 238000007918 intramuscular administration Methods 0.000 description 2
- 210000000653 nervous system Anatomy 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 206010019468 Hemiplegia Diseases 0.000 description 1
- 206010061296 Motor dysfunction Diseases 0.000 description 1
- 208000007101 Muscle Cramp Diseases 0.000 description 1
- 208000005392 Spasm Diseases 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000017531 blood circulation Effects 0.000 description 1
- 210000004556 brain Anatomy 0.000 description 1
- 230000005978 brain dysfunction Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 210000003710 cerebral cortex Anatomy 0.000 description 1
- 230000003920 cognitive function Effects 0.000 description 1
- 230000001447 compensatory effect Effects 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 230000001054 cortical effect Effects 0.000 description 1
- 230000002594 corticospinal effect Effects 0.000 description 1
- 238000002567 electromyography Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 210000000245 forearm Anatomy 0.000 description 1
- 230000001771 impaired effect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 210000000629 knee joint Anatomy 0.000 description 1
- 210000003141 lower extremity Anatomy 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 210000002161 motor neuron Anatomy 0.000 description 1
- 230000003387 muscular Effects 0.000 description 1
- 230000007433 nerve pathway Effects 0.000 description 1
- 230000001537 neural effect Effects 0.000 description 1
- 230000007230 neural mechanism Effects 0.000 description 1
- 210000002569 neuron Anatomy 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 210000000976 primary motor cortex Anatomy 0.000 description 1
- 210000003314 quadriceps muscle Anatomy 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 210000000278 spinal cord Anatomy 0.000 description 1
- 238000013179 statistical model Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H20/00—ICT specially adapted for therapies or health-improving plans, e.g. for handling prescriptions, for steering therapy or for monitoring patient compliance
- G16H20/30—ICT specially adapted for therapies or health-improving plans, e.g. for handling prescriptions, for steering therapy or for monitoring patient compliance relating to physical therapies or activities, e.g. physiotherapy, acupressure or exercising
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Physical Education & Sports Medicine (AREA)
- General Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Biophysics (AREA)
- General Engineering & Computer Science (AREA)
- Epidemiology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Primary Health Care (AREA)
- Public Health (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开一种基于混合Copula互信息的肌间耦合分析方法。本发明首先采用4种典型的单参数Copula函数建立混合Copula模型，然后利用惩罚似然函数的期望最大化(EM)和拟牛顿(BFGS)算法估计模型参数，再根据信息熵理论构造混合Copula互信息，最后将混合Copula互信息用于研究卒中后上肢及物运动时肱二头肌(BB)和肱三头肌(TB)在α、β和γ频段上的肌间耦合强度。本发明提出的混合Copula互信息不仅能度量肌间非线性耦合强度大小，而且还能全面刻画肌间依赖结构关系，对上肢运动功能评价具有良好的应用价值。
Description
技术领域
本发明属于神经系统运动控制机制研究领域，涉及混合Copula和互信息的计算，从而进行肌间功能耦合分析。
背景技术
脑卒中又称为中风，脑血管意外，是在脑血管病变或血流障碍基础上发生的局限性或弥漫性脑功能障碍，业已成为成人获得性运动障碍的主要原因。卒中后最常见的缺陷是对侧上/下肢偏瘫，由于上肢功能精细且灵巧，在康复治疗过程中，恢复较为困难，目前尚未找到有效便利的康复治疗技术。由于肌肉活动代表了神经系统的输出，对肌肉活动状态的检查能够了解运动障碍患者神经机制的灵活性和适应性的差异。现阶段在脑卒中患者康复过程中，运动功能评价的主要依据是康复医师的主观经验或依靠监测患者的肌力状态，难以客观、准确、定量评定康复效果。
肌间耦合的概念来源于皮层肌肉功能耦合研究，指的是运动过程中肌肉间的相互作用。肌间耦合在人体运动中起着重要的作用，决定了人体运动的整体肌肉模式，能够反映肌肉系统在中枢神经系统控制下的运动功能状态和信息交互方式。由于表面肌电信号(surface Electromyography,sEMG)非平稳、频域特征突出，在不同频段上耦合特征存在明显差异。肌间耦合被发现主要集中在α(8～15Hz)、β(15～30Hz)和γ(30～60Hz)频段，α频段的肌间耦合与运动过程中姿势的维持和非自主收缩有关，β频段的肌间耦合代表了从初级运动皮层到运动神经元的传递过程，γ频段的肌间耦合体现与认知功能的脑皮层信息整合过程。
近年来，时域的互信息(Mutual Information,MI)方法和频域的相干性(Coherence)方法，被广泛应用于脑肌电同步耦合分析之中。Chang等计算了行走过程中股四头肌间在α、β、高/低γ频段的交叉相关系数和相干性区域，实验结果表明，由于皮质脊髓传导的疲劳相关障碍，肌间交叉相关系数和β频带的一致性增强可能是在疲劳期间维持膝关节协同肌肉协调的一种代偿策略。Ouyang等提出了一种MI测度来描述手张开和手闭合状态下sEMG信号之间的相关性，实验结果表明，MI测量法能够有效地提取sEMG信号之间的相互关系，成功地区分不同的人手打开和关闭状态，提示MI能够从肌电信号中揭示肌间相互作用的本质特征。然而，相关系数(Correlation)和相干性只是表征两序列间是否线性的一个度量，并不能定量地描述两序列间相互依赖程度有多大，是否独立。基于概率分布的MI，对二阶和所有高阶相关系数都很敏感，能定量地表示出两序列间的相互依赖关系。然而，由于非高斯分布数据的联合概率密度函数(Joint Probability Density Function,JPDF)的估计是一个困难的问题，因此MI的估计往往十分困难。
Sklar定理指出，一个N维分量的联合分布函数可以由N个变量的边缘分布和1个Copula函数来描述。Copula函数称为“连接函数”或“相依函数”，它是把多个随机变量的联合分布与各自的边缘分布相连接的函数。Copula函数描述了变量间的相关性，具有诸多优点：1)构造灵活的多元分布，对边际分布没有特定要求；2)适用于严格单调增变换；3)刻画变量间的非线性、非对称相关性；4)如果边际分布函数连续，则Copula函数唯一确定；5)理论上不受维度限制。近年来，Copula理论作为一种灵活稳健的相关性分析工具，被广泛应用在金融风险度量、保险精算学、生物统计等领域。
Copula函数主要分为椭圆函数族(Ellipse Copula)和阿基米德函数(Archimedean Copula)两种类型。其中，椭圆Copula包括高斯Copula和t Copula，阿基米德Copula主要有Gumbel Copula，Clayton Copula和Frank Copula。不同类型的Copula函数具有不同的函数结构，不同形式的Copula模型导致不同的分析结果。现有研究在应用Copula函数时，大多采用单个Copula函数进行数据拟合，考虑到数据蕴含的复杂相关结构，仅用一种Copula函数可能会出现失真的情况。Hu等提出混合Copula函数(Mixed-Copula，M-Copula)的概念，建立了一个混合的Copula模型，使其能够成功捕获几个主要金融市场依赖结构的模式。M-Copula函数就是利用已知的一些Copula函数进行线性组合，形成新的Copula函数，能够克服单一Copula函数的缺点，具有多样性。Ma等证明了负Copula熵等价为MI，提出利用Copula函数理解和估计MI的新方法，通过Copula熵估计互信息方法简单，计算量小，能有效避免对联合密度函数地估计。
为了准确、定量描述不同特征功能频段上肌间线性和非线性耦合强度大小，本发明在Copula熵的基础上，将混合Copula与互信息有机结合，提出一种新的混合Copula互信息方法，并应用肌间耦合分析，旨在为运动功能康复评估提供新的研究方法和科学依据。
发明内容
本发明的目的在于提供一种可准确得到上肢肌间功能耦合特性的分析方法。
为实现上述目的，本发明方法主要包括以下步骤：
步骤(1)，表面肌电信号的同步采集以及预处理；
具体为：在表面肌电设备的监控下，同步采集上肢肱二头肌(Biceps，BB)和肱三头肌(Triceps，TB)上的2通道sEMG信号，采样频率为2000Hz。鉴于sEMG信号在低频段耦合特征显著，所以本发明重点关注α(8～15Hz)、β(15～30Hz)和γ(30～60Hz)3个特征功能频段。
步骤(2)，非参数核密度估计边际分布；
具体为：假设各频段上的sEMG信号是来自连续分布函数Fi(xi)的同分布样本，T为时间序列的长度，i＝1,2，那么Fi(xi)的非参数核密度估计为
其中，为概率密度函数，
其中，h为窗宽，当T→∞，h→0，Th→∞时，非参数核密度估计是真实概率密度分布的一致估计。K(·)为核函数，本发明采用高斯核函数，根据经验法则确定窗宽： 为xi的标准差。
步骤(3)，构建混合Copula模型；
具体为：混合Copula密度函数定义为：
其中，s为选择的Copula的个数，ωk为权重参数或称为形状参数，θk为相依参数，ck为第k个Copula密度函数，u,v～U[0,1]分别表示变量的边际分布函数。本发明采用4种常见的单参数的二元Copula函数用于构建混合Copula模型，分别为
(1)二元高斯Copula密度函数
其中，参数θGausssian∈(-1,1)为相关系数，二元高斯Copula函数适合对称的尾部特征，尾部渐进独立。
(2)二元Clayton Copula密度函数
其中，参数θClayton∈[-1,∞)\{0}，二元Clayton Copula函数适合不对称的尾部特征，对下尾部变化敏感。
(3)二元Frank Copula密度函数
其中，参数θFrank∈(-∞,∞)\{0}，二元Frank Copula函数适合对称的尾部特征，尾部渐进独立，对上下尾部变化都不敏感。
(4)二元Gumbel Copula密度函数
其中，参数θGumbel∈[0,∞)，二元Gumbel Copula函数适合不对称的尾部特征，对上尾部变化敏感。
步骤(4)，EM-BFGS算法估计混合Copula模型参数；
具体为：在对边际分布进行估计后，可得到估计值序列将其代入关于参数φ＝{ωT,θT}的对数似然函数L(φ)中
关于上式的求解，可在其后添加一个惩罚项，即
其中，γT为光滑参数，控制模型复杂度。为SCAD惩罚函数，满足无偏性、稀疏性和连续性。再利用期望最大化(EM)算法和拟牛顿(BFGS)算法估计ω*和θ*。
步骤(5)，模型拟合优度检验；
具体为：构建的混合Copula模型对肌间耦合关系的拟合优度的衡量标准可采用广泛使用的Akaike信息准则(AIC)和贝叶斯信息准则(BIC)：
AIC＝-2ln(L)+2k (10)
BIC＝-2ln(L)+k·ln(N) (11)
其中，k为模型中参数的个数，L为似然函数值。AIC、BIC值越小，模型拟合程度越好。
步骤(6)，计算混合Copula互信息；
具体为：根据互信息与Copula熵的等价关系，混合Copula互信息定义为
(12)
其中，称为混合Copula熵，可通过下式计算得到
步骤(7)，肌间耦合分析。
具体为：采用步骤(6)计算的混合Copula互信息度量特征功能频段上的肌间线性和非线性耦合强度。
本发明与传统的肌间耦合分析方法相比，具有如下优点：
传统的肌间耦合分析方法无法准确描述肌间线性耦合强度大小，针对该问题，本发明利提出一种新的混合Copula互信息方法，采用4种典型的单参数Copula函数建立混合Copula模型，利用惩罚似然函数的期望最大化(EM)和拟牛顿(BFGS)算法估计模型参数，再根据信息熵理论构造混合Copula互信息，这为探索运动过程肌间耦合特征及中枢神经系统的运动控制机制提供新方法。
附图说明
图1为本发明方法的流程图。
图2(a)、(b)、(c)为受试者H1的BB与TB分离出的α、β和γ频段波形。图2(d)、(e)、(f)为受试者S1的BB与TB分离出的α、β和γ频段波形。
图3(a)、(b)、(c)为受试者H1的3个频段波形的边际分布函数曲线。图3(d)、(e)、(f)为受试者S1的3个频段波形的边际分布函数曲线。
图4(a)、(b)、(c)为健康组和卒中组在α、β和γ频段的肌间耦合强度对比。
具体实施方式
下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方案和具体的操作过程。
卒中后，由于运动神经通路受损阻碍了神经振荡的传导，造成了运动功能障碍。卒中患者常常由于肌力弱、肌肉痉挛等原因无法完成及物动作，严重影响了人在日常生活中大部分时间的功能需求。肌间耦合可以反映来自脊髓中间神经元的共同驱动，研究运动过程中肌间耦合特性，挖掘潜在的中枢神经系统运动控制机制，对于上肢康复运动的功能评价具有十分重要的意义。
针对现有肌间耦合分析方法存在的问题，本发明将统计学中的Copula理论引入到肌间耦合分析之中。如图1所示，本发明的实施主要包括七个步骤：(1)表面肌电信号的同步采集以及预处理；(2)非参数核密度估计边际分布；(3)构建混合Copula模型；(4)EM-BFGS算法估计混合Copula模型参数；(5)模型拟合优度检验；(6)计算混合Copula互信息；(7)肌间耦合分析。最后，采用单因素一元方差分析方法进行组间显著性检验。
下面逐一对各步骤进行详细说明。
步骤一：表面肌电信号的同步采集以及预处理
采用表面肌电设备同步采集上肢2块肌肉上的sEMG信号，采样频率为2000Hz。预处理过程，首先手工提取5次有效活动段数据，然后通过上下采样的方式维持信号时长2.5s，接着进行去均值、去基线漂移，利用谱插值算法抑制50Hz工频干扰，最后采用4阶巴特沃斯带通滤波进行5-200Hz带通滤波，得到纯净的3通道sEMG信号。由于sEMG信号非平稳、频域特征突出，在低频段耦合特征存在明显差异，所以本发明重点关注α、β和γ这3个特征功能频段，通过快速傅里叶变换提取的频段波形如图2所示。
步骤二：非参数核密度估计边际分布
假设各频段上的sEMG信号是来自连续分布函数Fi(xi)的同分布样本，T为时间序列的长度，i＝1,2，那么Fi(xi)的非参数核密度估计为
其中，为概率密度函数，
其中，h为窗宽，当T→∞，h→0，Th→∞时，非参数核密度估计是真实概率密度分布的一致估计。K(·)为核函数，本发明采用高斯核函数，根据经验法则确定窗宽： 为xi的标准差。
步骤三：构建混合Copula模型
不同Copula密度函数的凸线性组合仍然是Copula密度函数，通过形状参数的变化可以对变量间复杂的相关结构进行更加细致的刻画，具有更大的灵活性和实用性。
混合Copula密度函数定义为：
其中，s为选择的Copula的个数，ωk为权重参数或称为形状参数，θk为相依参数，u,v～U[0,1]分别表示变量的边际分布函数。本发明采用4种常见的单参数的二元Copula函数用于构建混合Copula模型，分别为
(1)二元高斯Copula密度函数
其中，参数θGausssian∈(-1,1)为相关系数，二元高斯Copula函数适合对称的尾部特征，尾部渐进独立。
(2)二元Clayton Copula密度函数
其中，参数θClayton∈[-1,∞)\{0}，二元Clayton Copula函数适合不对称的尾部特征，对下尾部变化敏感。
(3)二元Frank Copula密度函数
其中，参数θFrank∈(-∞,∞)\{0}，二元Frank Copula函数适合对称的尾部特征，尾部渐进独立，对上下尾部变化都不敏感。
(4)二元Gumbel Copula密度函数
其中，参数θGumbel∈[0,∞)，二元Gumbel Copula函数适合不对称的尾部特征，对上尾部变化敏感。
步骤四：EM-BFGS算法估计混合Copula模型参数
在对边际分布进行估计后，可得到估计值序列将其代入关于参数φ＝{ωT,θT}的对数似然函数L(φ)中
关于上式的求解，可在其后添加一个惩罚项，即
其中，γT为光滑参数，控制模型复杂度。为SCAD惩罚函数，满足无偏性、稀疏性和连续性。利用期望最大化(EM)算法和拟牛顿(BFGS)算法估计ω*和θ*。
步骤五：模型拟合优度检验
在Copula的应用中，一个重要的问题是如何从一组给定的候选Copula中选择最合适的Copula来描述依赖结构。一种常用的方法是Akaike信息准则(AIC)。Akaike信息准则是衡量统计模型拟合优度的标准，其定义为：
AIC＝-2ln(L)+2k (23)
另一个常用的准则是贝叶斯信息准则(BIC)，BIC定义为：
BIC＝-2ln(L)+k·ln(N) (24)
其中，k为模型中参数的个数，L为似然函数值。AIC、BIC值越小，模型拟合程度越好。
步骤六：计算混合Copula互信息
二元Copula熵定义为
Hc(u,v)＝-∫∫c(u,v)×log(c(u,v))dudv (25)
具有以下性质
Hc(u,v)＝H(x,y)-H(x)-H(y)＝-I(x,y) (26)
其中H(*)表示Shannon熵，H(*,*)表示联合熵，I(*,*)表示互信息。式(26)表明负的Copula熵等价于互信息。显然，不同的Copula函数估计出的互信息是不同的。将混合Copula密度函数代入式(25)或式(26)中，可得到混合Copula熵，取负后即为混合Copula互信息
关于混合Copula熵的估计可按式(28)进行计算[21]
其中，E[*]表示求期望操作。当对数运算以2为基底，单位为比特(Bit)。
步骤七：肌间耦合分析
采用步骤六计算的混合Copula互信息度量特征功能频段上的肌间线性和非线性耦合强度，采用单因素一元方差分析方法进行组间显著性检验。
为了验证本发明方法的性能，实验部分采集了多个受试者的数据，通过性能指标来衡量本发明在肌间耦合分析领域的可行性。实验招募了共10名受试者参与实验，对照组：5名健康被试(H1～H5)，平均年龄小于76岁；实验组：5名卒中患者(S1～S5)，中风后少于33.7天，显示轻度运动障碍，Fugl-Meyer平均评分高于43.732分(总分为66分)。要求每位受试者坐在桌子前，前臂放在一个舒服的位置上，根据每隔10秒激活的语音提示，要求向着正前方目标进行5次及物动作，目标距桌子高35cm，每次及物后休息10s，受试者在执行运动任务前已了解实验内容。采用表面肌电设备同步采集上肢2通道肱二头肌(Biceps，BB)、肱三头肌(Triceps，TB)上的sEMG信号，采样频率为2000Hz。
利用快速傅里叶变换，提取BB和TB的α、β和γ三个频段上的子波。其中，健康组被试H1和卒中组被试S1结果如图2所示。由图2可见，H1和S1的TB与BB呈现出错峰活动，BB的激活程度要低于TB。H1的TB在0.5～1s附近显著激活，S1的TB在α、β和γ频段上激活时间要长于H1，S1的BB在α频段的激活程度要高于H1。
图3给出了非参数核密度估计边际分布函数的估计结果，从图3可以看出，H1和S1的BB与TB在不同频段子波的累计概率密度函数(Cumulative Density Function,CDF)都具有一定对称性，与高斯分布的CDF比较接近，其中BB较TB的尾部占比更小，方差更小。
为凸显混合Copula模型的灵活性，在通过IFM算法估计单个Copula函数的参数后，利用EM算法和BFGS算法估计混合Copula函数的相关参数和权重参数，按式(10)和(11)计算AIC、BIC拟合优度检验指标，健康组被试H1～H5和卒中组被试S1～S5的5次平均结果如表1所示。
表1 Copula、混合Copula参数估计以及AIC、BIC指标(平均值)
注：Sub.表示Subjects，Ban.表示Bands，Ga.表示Gaussian，Cl.表示Clayton，Fr.表示Frank，Gu.表示Gumbel。
由表1可知，在健康组中，Gaussian Copula和Frank Copula在单独运用时表现效果很好{α：H1,H4,H5；β：H1～H5；γ：H2～H5}，具有较低的AIC值和BIC值。在α频段混合Copula拟合效果是最好的(H1～H5)，Gaussian Copula的权重最大(H1,H3,H4)，影响最大，但在β频段和γ频段上混合Copula的性能则不如Gaussian Copula和Frank Copula(H1,H3,H4,H5)，其中，β频段上Gumbel Copula的权重最大(H1,H2,H4,H5)，γ频段Gaussian Copula的权重最大(H1,H2,H4,H5)。在卒中组中，Gaussian Copula和Frank Copula在单独运用时表现效果较好{α：S2,S3；β：S1,S3,S5；γ：S2,S3,S5}，但都不如混合Copula{α：S1～S5；β：S1～S4；γ：S1～S5}，AIC值和BIC值更低，在混合Copula中，Gaussian Copula的权重最大{α：S3,S5；β：S1,S4,S5；γ：S1～S5}。
为表明混合Copula互信息的优良估计性能，本发明与以上4种单一的Copula模型获得的Copula互信息相比较，均用于表征BB和TB分离出的3个频段信号间的耦合强度，结果以误差柱状图的形式给出，如图4所示。由图4可见，健康组和卒中组在α频段的平均肌间耦合强度是较为接近的。混合Copula互信息值最大(>0.04Bit)，其次是Gaussian Copula。β、γ频段的平均肌间耦合强度都较低，卒中组比健康组的混合Copula互信息值更高(>0.02Bit)，肌间耦合更强。混合Copula互信息较单个Copula互信息更容易区分组间肌间耦合强度的差异。
利用单因素一元方差分析方法检验组间肌间耦合强度差异是否具有统计学意义，结果如表2所示。由表2可知，卒中组在β和γ频段平均肌间耦合强度高于健康组(β:0.0242Bit>0.0104Bit,γ:0.0268Bit>0.0116Bit)，但具有较大的方差，3个频段的肌间耦合强度不具有显著性差异(p>0.05)。
表2各频段组间单因素一元方差分析(平均值±标准差，Bit)
以上所述的实施例仅仅是对本发明的优选实例方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。
Claims (2)
1.基于混合Copula互信息的肌间耦合分析方法，其特征在于：该方法包括以下步骤：
步骤(1)，表面肌电信号的同步采集以及预处理；
具体为：在表面肌电设备的监控下，同步采集上肢肱二头肌和肱三头肌上的2通道表面肌电信号，采样频率为2000Hz；
步骤(2)，非参数核密度估计边际分布；
具体为：假设各频段上的表面肌电信号是来自连续分布函数Fi(xi)的同分布样本，T为时间序列的长度，i＝1,2，那么Fi(xi)的非参数核密度估计为
其中，为概率密度函数；
其中，h为窗宽，当T→∞，h→0，Th→∞时，非参数核密度估计是真实概率密度分布的一致估计，K(·)为核函数；
步骤(3)，构建混合Copula模型；
具体为：定义混合Copula密度函数：
s.t.∑ωk＝1,ωk≥0
其中，s为选择的Copula的个数，ωk为权重参数，θk为相依参数，ck为第k个Copula密度函数，u,v～U[0,1]分别表示变量的边际分布函数，采用单参数的二元Copula函数用于构建混合Copula模型；
步骤(4)，EM-BFGS算法估计混合Copula模型参数；
具体为：在对边际分布进行估计后，可得到估计值序列将其代入关于参数φ＝{ωT,θT}的对数似然函数L(φ)中
关于上式的求解，其后添加一个惩罚项，即
其中，γT为光滑参数，控制模型复杂度；为SCAD惩罚函数，满足无偏性、稀疏性和连续性；再利用期望最大化算法和拟牛顿算法估计ω*和θ*；
步骤(5)，模型拟合优度检验；
具体为：构建的混合Copula模型对肌间耦合关系的拟合优度的衡量标准采用广泛使用的Akaike信息准则和贝叶斯信息准则：
步骤(6)，计算混合Copula互信息；
具体为：根据互信息与Copula熵的等价关系，混合Copula互信息定义为：
其中，称为混合Copula熵，可通过下式计算得到
步骤(7)，肌间耦合分析；
具体为：采用步骤(6)计算的混合Copula互信息度量特征功能频段上的肌间线性和非线性耦合强度。
2.根据权利要求1所述的基于混合Copula互信息的肌间耦合分析方法，其特征在于：
步骤(1)所述的预处理过程具体为：首先手工提取5次有效活动段数据，然后通过上下采样的方式维持信号时长2.5s，接着进行去均值、去基线漂移，利用谱插值算法抑制50Hz工频干扰，最后采用4阶巴特沃斯带通滤波进行5-200Hz带通滤波，得到纯净的3通道表面肌电信号。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010715642.5A CN111931606B (zh) | 2020-07-23 | 2020-07-23 | 基于混合Copula互信息的肌间耦合分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010715642.5A CN111931606B (zh) | 2020-07-23 | 2020-07-23 | 基于混合Copula互信息的肌间耦合分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111931606A CN111931606A (zh) | 2020-11-13 |
| CN111931606B true CN111931606B (zh) | 2023-10-20 |
Family
ID=73315284
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010715642.5A Active CN111931606B (zh) | 2020-07-23 | 2020-07-23 | 基于混合Copula互信息的肌间耦合分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111931606B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112232301A (zh) * | 2020-11-16 | 2021-01-15 | 杭州电子科技大学 | 基于多尺度Copula互信息的肌间耦合网络分析方法 |
| CN112509689B (zh) * | 2021-02-08 | 2024-05-17 | 杭州电子科技大学 | 基于时变Copula互信息的肌间耦合分析方法 |
| CN113440150A (zh) * | 2021-06-25 | 2021-09-28 | 睿旭康(苏州)智能技术有限责任公司 | 基于R-vine Copula的皮层肌肉功能网络构建方法 |
| CN114137832B (zh) * | 2021-10-26 | 2024-03-29 | 杭州电子科技大学 | 基于R藤Copula传递熵的多变量因果关系方法 |
| CN114159081A (zh) * | 2021-12-13 | 2022-03-11 | 杭州电子科技大学 | 基于时变混合Copula互信息的肌电耦合方法 |
| CN115474945B (zh) * | 2022-09-15 | 2024-04-12 | 燕山大学 | 一种面向多通道脑肌电耦合分析的多元全局同步指数方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108345887A (zh) * | 2018-01-29 | 2018-07-31 | 清华大学深圳研究生院 | 图像语义分割模型的训练方法及图像语义分割方法 |
| CN108742613A (zh) * | 2018-05-30 | 2018-11-06 | 杭州电子科技大学 | 基于传递熵与广义偏定向相干性的肌间耦合分析方法 |
| CN109497999A (zh) * | 2018-12-20 | 2019-03-22 | 杭州电子科技大学 | 基于Copula-GC的脑肌电信号时频耦合分析方法 |
| US10453197B1 (en) * | 2019-02-18 | 2019-10-22 | Inception Institute of Artificial Intelligence, Ltd. | Object counting and instance segmentation using neural network architectures with image-level supervision |
| CN110638444A (zh) * | 2019-09-05 | 2020-01-03 | 杭州电子科技大学 | 基于MEMD-rTVgPDC的皮层肌肉耦合分析方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9672193B2 (en) * | 2013-03-15 | 2017-06-06 | Sas Institute Inc. | Compact representation of multivariate posterior probability distribution from simulated samples |
-
2020
- 2020-07-23 CN CN202010715642.5A patent/CN111931606B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108345887A (zh) * | 2018-01-29 | 2018-07-31 | 清华大学深圳研究生院 | 图像语义分割模型的训练方法及图像语义分割方法 |
| CN108742613A (zh) * | 2018-05-30 | 2018-11-06 | 杭州电子科技大学 | 基于传递熵与广义偏定向相干性的肌间耦合分析方法 |
| CN109497999A (zh) * | 2018-12-20 | 2019-03-22 | 杭州电子科技大学 | 基于Copula-GC的脑肌电信号时频耦合分析方法 |
| US10453197B1 (en) * | 2019-02-18 | 2019-10-22 | Inception Institute of Artificial Intelligence, Ltd. | Object counting and instance segmentation using neural network architectures with image-level supervision |
| CN110638444A (zh) * | 2019-09-05 | 2020-01-03 | 杭州电子科技大学 | 基于MEMD-rTVgPDC的皮层肌肉耦合分析方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111931606A (zh) | 2020-11-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111931606B (zh) | 基于混合Copula互信息的肌间耦合分析方法 | |
| Wang et al. | Recognition of muscle fatigue status based on improved wavelet threshold and CNN-SVM | |
| Chen et al. | Adaptive real-time identification of motor unit discharges from non-stationary high-density surface electromyographic signals | |
| CN204931634U (zh) | 基于生理信息的抑郁症评估系统 | |
| CN111708978B (zh) | 多尺度时频肌间耦合分析方法 | |
| CN112130668B (zh) | 一种R藤Copula互信息的肌间耦合分析方法 | |
| CN112509689B (zh) | 基于时变Copula互信息的肌间耦合分析方法 | |
| CN109674445B (zh) | 一种结合非负矩阵分解和复杂网络的肌间耦合分析方法 | |
| CN110638444A (zh) | 基于MEMD-rTVgPDC的皮层肌肉耦合分析方法 | |
| CN111931129A (zh) | 基于高斯Copula传递熵的肌间耦合网络分析方法 | |
| Majid et al. | EMG feature extractions for upper-limb functional movement during rehabilitation | |
| Liu et al. | The progression of muscle fatigue during exercise estimation with the aid of high-frequency component parameters derived from ensemble empirical mode decomposition | |
| CN114648040A (zh) | 生命体征多生理信号提取、融合分析方法 | |
| Al-Taee et al. | Optimal feature set for finger movement classification based on sEMG | |
| Arjunan et al. | Fractal properties of surface electromyogram for classification of low-level hand movements from single-channel forearm muscle activity | |
| CN113558638A (zh) | 基于Vine Copula的脑肌耦合模型构建方法 | |
| CN117281479A (zh) | 一种基于表面肌电信号多维度特征融合的人体下肢慢性疼痛判别方法、存储介质及装置 | |
| CN114052750B (zh) | 基于标准模板肌电分解的脑肌信息传递规律提取方法 | |
| Zhao et al. | Online Decomposition of Surface Electromyogram into Individual Motor Unit Activities Using Progressive FastICA Peel-off | |
| CN111543990B (zh) | 肌肉协同的提取方法及基于肌肉协同的运动功能分析评估方法 | |
| Cahyadi et al. | Analysis of EMG based arm movement sequence using mean and median frequency | |
| Diez et al. | A comparative study of the performance of different spectral estimation methods for classification of mental tasks | |
| Sburlea et al. | Predicting EMG envelopes of grasping movements from EEG recordings using unscented kalman filtering | |
| Long et al. | Estimation of fatigue status by sEMG signal using SVM algorithm in massage assessment | |
| Zhao et al. | Evaluation of lower-limb rehabilitation based on muscle fatigue |
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 |