CN114767132B - 一种基于表面肌电分解的运动单位识别方法 - Google Patents

一种基于表面肌电分解的运动单位识别方法 Download PDF

Info

Publication number
CN114767132B
CN114767132B CN202210417737.8A CN202210417737A CN114767132B CN 114767132 B CN114767132 B CN 114767132B CN 202210417737 A CN202210417737 A CN 202210417737A CN 114767132 B CN114767132 B CN 114767132B
Authority
CN
China
Prior art keywords
motion unit
unmixed
vector
source signal
iteration
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
CN202210417737.8A
Other languages
English (en)
Other versions
CN114767132A (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.)
University of Science and Technology of China USTC
Original Assignee
University of Science and Technology of China USTC
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 University of Science and Technology of China USTC filed Critical University of Science and Technology of China USTC
Priority to CN202210417737.8A priority Critical patent/CN114767132B/zh
Publication of CN114767132A publication Critical patent/CN114767132A/zh
Application granted granted Critical
Publication of CN114767132B publication Critical patent/CN114767132B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/24Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316Modalities, i.e. specific diagnostic methods
    • A61B5/389Electromyography [EMG]
    • 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
    • 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
    • A61B5/7203Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Surgery (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Physics & Mathematics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Pathology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Signal Processing (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Psychiatry (AREA)
  • Artificial Intelligence (AREA)
  • Image Analysis (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明公开了一种基于表面肌电分解的运动单位识别方法,包括:1、预采集指定肌肉运动下的表面肌电信号,构建用以初始化的离线数据集;2、采用离线状态的逐步变量剥离算法分解离线数据集,解算每个运动单位对应的解混向量并保存;3、在线处理阶段,实时获取表面肌电信号并对其分窗,利用已获取的解混向量直接得到不同的运动单位源信号;4、利用自适应修正阈值的方法提取分析窗内源信号对应的运动单位发放序列;5、随着实时数据的输入,执行步骤3和4,待没有新的数据输入时,依次连接不同分析窗的分解结果并输出。本发明能以在线的方式识别表面肌电信号中的运动单位活动,并实时获取到运动单位的发放序列。

Description

一种基于表面肌电分解的运动单位识别方法
技术领域
本发明属于肌电信号处理技术领域,具体涉及一种基于肌电分解技术的运动单位识别方法,主要应用于高效智能的神经指令解码和运动意图捕获。
背景技术
表面肌电信号(Surface Electromyography,SEMG)作为一种伴随运动产生的电生理信号,反映驱动肌肉收缩运动的神经指令。其形成过程由多个运动单位(Motor Unit,MU)受中枢调控激活发放产生的运动单位动作电位(Motor UnitAction Potential,MUAP)沿肌纤维传播在检测电极处的时间、空间叠加而成的,直接反应了用户的运动意图和神经驱动信息。特别的,随着肌电采集技术的进步,包含有丰富的时空信息的高密度表面肌电为肌电控制应用注入了新的活力。更细致准确的解析肌电运动信息,特别是从刻画肌肉整体活动的宏观层面深入至探究支配肌肉活动的神经指令的微观层面,是肌电控制技术革新所面临的根本性问题。因此,现有的肌电控制技术不再满足于各项宏观时频域特征的简单应用,转而聚焦于蕴含在神经驱动下由大量的信源产生的丰富活动信息,寻找更符合运动底层生理规律的控制策略。在解析微观神经驱动信息的过程中,需要获取神经肌肉系统最基本的组成单位——运动单位的活动信息,因此肌电分解是不可或缺的一步。肌电分解可以看作肌电形成的逆过程,将混叠的表面肌电还原为其最基本的组成部分——MUAP序列的形式,帮助我们从运动单位的微观层面获取神经驱动信息。MUAP序列中包含了运动单位的发放和波形信息,其中发放信息是最为重要的,直接展示了神经指令下发和传达的时刻,表征为运动单位的0-1发放序列。大量的学者利用肌电分解算法识别出肌肉活动中运动单位的信息,做出了很多深入的研究,也驱动着肌电控制技术朝着解码微观神经指令的方向逐步发展。
在过去的十年内,最具代表性的肌电分解算法包括卷积核抵消算法(ConvolutionKernel Compensation,CKC)和逐步变量剥离算法(Progressive FastICAPeel-off,PFP)。CKC算法粗略地估计运动单位发放序列和表面肌电信号之间的互相关矢量,然后通过估计的互相关矢量和EMG信号的相关矩阵来计算该运动单位的发放模式。PFP算法通过盲源分离经典算法快速独立成分分析(Fast Independent Component Analysis,FastICA)算法求解出一系列运动单位源信号,提取出其中包括的发放序列,通过带时域约束的FastICA对结果进行检测和修正,并通过剥离策略将已得到的运动单位从原始信号中剥离出来,以寻找到更多的运动单位。两种肌电分解策略在离线的过程中取得了较好的成果,其分解得到的运动单位活动信息是解码表面肌电的神经指令的重要依据,直接反映了用户的运动意图和运动状态,提供了重要的人机接口,应用潜力巨大。
然而,目前的以肌电分解技术为基础的运动单位识别方法大都是在离线的过程中处理,鲜有在在线情况的研究,但是在科技飞速发展的现代社会,对信号的实时处理分析有着巨大的需求,离线状态的运动单位识别显然存在很严重的应用局限性。表面肌电信号呈现为动作电位波形的高度混叠,从中分解出代表微观神经指令的单个运动单位信息已然是领域内公认的难题。对于离线的表面肌电分解技术,需要对较长时间的信号段进行复杂的迭代计算求出解混向量,依赖专家检测的过程来提升分解结果的准确性,计算复杂性较高,确保性能的前提下难以满足肌电控制任务所需的实时获取分解结果的要求,这对于在线运动单位识别的是一个巨大的困难与挑战。
发明内容
本发明为了克服现有运动单位识别方法的不足之处,提出一种基于表面肌电分解的运动单位在线识别方法,以期能够在实时输入肌电信号的情形下,以在线的方式识别表面肌电信号中的运动单位活动,并实时获取到运动单位的发放序列,同时保持较高的精度,从而达到与离线分解相当的效果。
本发明为解决技术问题,采用如下技术方案:
本发明一种基于表面肌电分解的运动单位识别方法的特点是按如下步骤进行:
步骤1、利用肌电测量设备和M通道高密度电极在待测肌肉处预采集t时刻肌肉收缩任务的表面肌电信号数据,记为x(t)=[x1(t),x2(t),…,xI(t),…,xM(t)]T,其中,xI(t)表示t时刻第I个通道的肌电数据;从而得到时长为T1的表面肌电信号数据并作为离线数据集;
步骤2、采用逐步变量剥离算法分解所述离线数据集的肌电信号:
步骤2.1、利用式(1)对所述离线数据集进行扩展,得到扩展信号
式(1)中,K是延迟因子;
定义并初始化运动单位发放集合γ及其对应的解混向量集合ε为空集;
定义并初始化残差信号
步骤2.2、执行FastICA算法,设置算法执行的总次数为n:
步骤2.2.1、初始化i=1;
步骤2.2.2、初始化j=0;
步骤2.2.3、定义并初始化残差信号中第i个解混向量在第j次的计算结果wi,j为一个均值为0,方差为1的高斯信号;
步骤2.2.4、利用式(2)计算第i个解混向量在第j+1次的计算结果wi,j+1
式(2)中,G'表示非多项式函数G的一阶导数;G”表示G的二阶导数;E是数学期望;代表第i个解混向量在第j次的计算结果wi,j的转置;
步骤2.2.5、如果表示预设的迭代计算次数上限,则令j+1赋值给j,返回至步骤2.2.3;若i<n且/>则得到第i个解混向量wi=wi,j,则令i+1赋值给i,返回至步骤2.2.2;若i≥n且/>则得到n个解混向量及其组成的解混矩阵W=[w1,w2,…,wi,…,wn]T,并利用式(3)得到正交化后的解混矩阵W':
W'=(WWT)-1/2W (3)
步骤2.2.6、计算运动单位源信号且Y=[y1,y2,…,yi,…,yn]T,其中,yi代表第i个解混向量wi所对应的运动单位源信号;
步骤2.3:采用图像处理的大津法获取第i个运动单位源信号yi的阈值thi
步骤2.3.1、将区间(0,max(yi))划分为M个等分,得到[thi,1,thi,2,thi,3,…,thi,m,…,thi,M];其中,max(yi)表示运动单位源信号yi的最大值,thi,m表示第m个小区间的最大值;
步骤2.3.2、初始化h=1;
步骤2.3.3、令thi=thi,h;求出yi中大于thi的采样点个数N0,h和小于thi的采样点个数N1,h;并求解个数N0,h对应的灰度μ0,h作为yi中所有大于thi的采样点的平均值;求解个数N1,h对应的灰度μ1,h作为yi中所有小于thi的采样点的平均值;
步骤2.3.4、根据式(4)计算得到类间方差Fi,h
Fi,h=N0,h×N1,h×(μ0,h1,h)2 (4)
步骤2.3.5、若h<M,则令h+1赋值给h,返回至步骤2.3.3;否则进入步骤2.3.6;
步骤2.3.6、找到最大的Fi,h对应的下标hmax,并将下标hmax所对应的thi,hmax作为第i个运动单位源信号yi的阈值thi
步骤2.4、将yi中大于阈值thi的波峰峰值位置置为1,其余置为0,从而得到对应的运动单位发放序列ri
步骤2.5、通过带时域约束的FastICA算法检测分解结果的可靠性:
步骤2.5.1、定义迭代计算的次数为d,并初始化d=1;定义并初始化第d次迭代的拉格朗日乘子μd=0、第d次迭代的惩罚因子δd为随机小数;初始化第d次迭代的第i个运动单位源信号yi,d=yi、第d次迭代的第i个解混向量wi,d=wi;定义第i个运动单位源信号yi和待测发放序列ri之间的相关性函数为其中,ξ是一个预设的相关性下界,初始化第d次迭代的相关性下界ξd=1、第d次迭代的相关性函数/>
步骤2.5.2、通过式(5)得到第d+1次迭代的第i个解混向量wi,d+1
式(5)中,g'd为第d次迭代的相关性函数gd的一阶导数;
步骤2.5.3、对第d+1次迭代的第i个解混向量wi,d+1进行标准化:将wi,d+1/||wi,d+1||2赋予wi,d+1,其中||wi,d+1||2是wi,d+1的2-范数;
步骤2.5.4、根据式(6)、式(7)和式(8)分别更新其他参数:
μd+1=max{0,μddgd} (6)
δd+1=α×δd (7)
ξd+1=β×ξd (8)
式(6)、式(7)和式(8)中,μd+1表示第d+1次迭代的拉格朗日乘子,δd+1表示第d+1次迭代的惩罚因子,ξd+1表示第d+1次迭代的相关性下界;
步骤2.5.5、更新第d+1次迭代的第i个运动单位源信号更新第d+1次迭代的相关性函数/>
步骤2.5.6、若表示预设的迭代计算次数上限,则令d+1赋值给d,返回步骤2.5.2;否则,进行步骤2.5.7;
步骤2.5.7、令第i个运动单位源信号yi=yi,d+1,采用步骤2.3到步骤2.4的过程获取新的发放序列r'i;计算ri和r'i的的互相关系数ci
若ci>C,则表示更新后的运动单位源信号yi可靠,并将新的序列r'i和其对应的解混向量wi,d+1分别放入运动单位发放集合γ和解混向量集合ε中,其中,C表示预设的用来衡量相关性的相关度界限;否则,将新的序列r'i和其对应的解混向量wi,d+1直接删除;返回步骤2.5.1对源信号第i+1个运动单位源信号yi+1继续检测更新,直至最后一个运动单位源信号完成更新;
步骤2.6、遍历运动单位发放集合γ并分别计算任意两个不同的发放序列之间的相关系数,若大于C,则表示两个发放序列重复,并随机选取一个从运动单位发放集合γ和解混向量集合ε中删除;否则,两个发放序列均保留;经此操作后,可得到共有q个发放序列的运动单位发放集合γ和共有q个解混向量的解混向量集合ε;
步骤2.7、将运动单位发放集合γ中q个发放序列表示为R=[r1,r2,…,rq]T并计算其在最小二乘意义下的最优波形从而更新残差信号/>
步骤2.8、将更新后的残差信号代入步骤2.2-步骤2.7的过程进行处理,并判断前后两次处理得到的集合γ和ε中是否有新增的结果,若有,则继续返回至步骤2.2顺序执行,否则停止算法,并最终得到由离线数据集经过分解得到后的运动单位发放集合γ和包含Q个解混向量的解混向量集合ε;由Q个解混向量组成解混矩阵其中,/>表示解混矩阵/>中的第i个解混向量;
步骤3、实时采集肌电信号,并将获取的肌电数据流按照时间轴分为一系列窗长为L,步长为D的分析窗;在每个分析窗内,根据式(1)得到扩展后的肌电数据X,并利用卷积模型 求解运动单位源信号所组成的矩阵S=[s1,s2,…,si,…,sQ]T,si代表其中第i个运动单位源信号;
步骤4、采用自适应的阈值选取策略提取分析窗内第i个运动单位源信号si对应的发放序列并保存;
步骤5、每当新的肌电数据输入时,按照步骤3和步骤4的过程进行处理和保存;
当没有新的数据输入时,利用式(9)将所保存的P个分析窗中的发放序列按照时间顺序进行拼接,得到总的分解结果Spike并作为最终的结果:
Spike=[Win1,Win'2,Win'3,…,Win'i,…,Win'P] (9)
式(9)中,Wini为第i个窗的分解的所有发放序列;Win'i表示Wini的结尾长度为D的数据。
本发明所述的基于表面肌电分解的运动单位识别方法的特点也在于,所述步骤4中的自适应阈值选择策略是按如下过程进行:
步骤4.1、按照步骤2.3中的大津法计算第i个运动单位源信号si的初始阈值th0
步骤4.2、将si中所有大于th0的波峰的幅值置于集合β中,从而得到共有l个波峰幅值的集合β,记作β=[peak1,peak2,…,peakk,…,peakl];其中,peakk代表第k个波峰的幅值;
步骤4.3、初始化k=1;
步骤4.4、令中间变量th=peakk,将si中大于阈值th的波峰峰值位置置为1,其余置为0,从而得到对应的运动单位发放序列bk
步骤4.5、设置向量a,用来存储运动单位发放序列bk中每相邻两数值为1的位置之间的位置差;利用式(10)计算第k个发放变异系数cosk
cosk=std(a)/mean(a) (10)
式(10)中,std(a)表示向量a的标准差,mean(a)表示向量a的平均值;
步骤4.6、若k<l,则令k+1赋值给k,返回步骤4.4;否则找到最大的cosk对应的下标kmin,并将下标kmin所对应的发放序列作为第i个运动单位源信号si对应的发放序列。
本发明能够实时识别表面肌电包含的运动单位,且能够以自适应的方式提取运动单位的发放序列,达到离线肌电分解相当的准确率,具体的有益效果体现在:
1、本发明步骤二对离线数据集进行了分解,获取到了不同运动单位对应的解混向量。基于运动单位结构与位置不变性,步骤三直接将解混向量集合作用于卷积模型求解新获取的肌电信号中的运动单位源信号,从而避免了庞大的迭代计算的步骤,将在线处理的部分计算工作转移到离线的过程中,为在线处理加速。
2、本发明步骤四中通过优化运动单位辨别与提取方法,采用自适应的阈值选取策略,克服了离线肌电分解需要专家检测的弊端,保证了所提的方法能够兼顾高性能和实时性。
3、本发明步骤三中通过设置随时间滑动的分析窗,保证了肌电数据实时输入和分解结果的实时输出,满足了在线处理数据的需求;步骤六通过相关系数的大小链接不同分析窗内的分解结果,保证了分解结果在时间上的连续性。
附图说明
图1为本发明实施例提供的方法流程示意图;
图2为本发明实施例提供的实验示意图;
图3为本发明实施例提供的在线处理和离线分解结果对比示意图;
图4为本发明实施例提供的在线处理结果示意图。
具体实施方式
本实施例中,一种基于表面肌电分解的运动单位识别方法,如图1所示,包括如下步骤:
步骤1、利用肌电测量设备和M通道高密度电极在待测肌肉处预采集t时刻肌肉收缩任务的表面肌电信号数据,记为x(t)=[x1(t),x2(t),…,xI(t),…,xM(t)]T,其中,xI(t)表示t时刻第I个通道的肌电数据;从而得到时长为T1的表面肌电信号数据并作为离线数据集;
具体实施中包括:
(1)募集d位受试者,引导每位受试者将待测手放于3D打印模块上固定,拇指套入连接着力度传感器的戒指环,设备采集手掌大鱼际肌的M通道高密度表面肌电信号,如图2中的(a)部分所示。阵列电极阵列中单个电极触点的直径为p,电极中心间距为q。示例性的,可以设置:d=8,M=64,p=2mm,q=4mm。受试者记为S1-S8;
(2)逐一采集受试者执行拇指收缩任务的连续肌电信号。肌肉收缩任务设置为拇指肌力在2s的时间内从0逐渐增加至30%最大收缩力(在每次实验前测量不同受试者的拇指最大收缩力),然后维持这个力度约3s,最后直接下降至0,如图2中的(b)部分所示。每个手势任务采集5-10次,保证有足够的数据用来测试;
(3)将所采集的数据作为离线数据保存,通常选取5次收缩任务的数据,以保证能够获取到足够的解混向量。故离线数据集时长T1为25s;
步骤2、采用逐步变量剥离算法分解所述离线数据集的肌电信号:
步骤2.1、利用式(1)对所述离线数据集进行扩展,得到扩展信号
式(1)中,K是延迟因子;
定义并初始化运动单位发放集合γ及其对应的解混向量集合ε为空集;
定义并初始化残差信号
步骤2.2、执行FastICA算法,设置算法执行的总次数为n:
步骤2.2.1、初始化i=1;
步骤2.2.2、初始化j=0;
步骤2.2.3、定义并初始化残差信号中第i个解混向量在第j次的计算结果wi,j为一个均值为0,方差为1的高斯信号;
步骤2.2.4、利用式(2)计算第i个解混向量在第j+1次的计算结果wi,j+1
式(2)中,G'表示非多项式函数G的一阶导数;G”表示G的二阶导数;E是数学期望;代表第i个解混向量在第j次的计算结果wi,j的转置;
步骤2.2.5、如果表示预设的迭代计算次数上限,则令j+1赋值给j,返回至步骤2.2.3;若i<n且/>则得到第i个解混向量wi=wi,j,则令i+1赋值给i,返回至步骤2.2.2;若i≥n且/>则得到n个解混向量及其组成的解混矩阵W=[w1,w2,…,wi,…,wn]T,并利用式(3)得到正交化后的解混矩阵W':
W'=(WWT)-1/2W (3)
步骤2.2.6、计算运动单位源信号且Y=[y1,y2,…,yi,…,yn]T,其中,yi代表第i个解混向量wi所对应的运动单位源信号;
步骤2.3:采用图像处理的大津法获取第i个运动单位源信号yi的阈值thi
步骤2.3.1、将区间(0,max(yi))划分为M个等分,得到[thi,1,thi,2,thi,3,…,thi,m,…,thi,M];其中,max(yi)表示运动单位源信号yi的最大值,thi,m表示第m个小区间的最大值;
步骤2.3.2、初始化h=1;
步骤2.3.3、令thi=thi,h;求出yi中大于thi的采样点个数N0,h和小于thi的采样点个数N1,h;并求解个数N0,h对应的灰度μ0,h作为yi中所有大于thi的采样点的平均值;求解个数N1,h对应的灰度μ1,h作为yi中所有小于thi的采样点的平均值;
步骤2.3.4、根据式(4)计算得到类间方差Fi,h
Fi,h=N0,h×N1,h×(μ0,h1,h)2 (4)
步骤2.3.5、若h<M,则令h+1赋值给h,返回至步骤2.3.3;否则进入步骤2.3.6;
步骤2.3.6、找到最大的Fi,h对应的下标hmax,并将下标hmax所对应的作为第i个运动单位源信号yi的阈值thi
步骤2.4、将yi中大于阈值thi的波峰峰值位置置为1,其余置为0,从而得到对应的运动单位发放序列ri
步骤2.5、通过带时域约束的FastICA算法检测分解结果的可靠性:
步骤2.5.1、定义迭代计算的次数为d,并初始化d=1;定义并初始化第d次迭代的拉格朗日乘子μd=0、第d次迭代的惩罚因子δd为随机小数;初始化第d次迭代的第i个运动单位源信号yi,d=yi、第d次迭代的第i个解混向量wi,d=wi;定义第i个运动单位源信号yi和待测发放序列ri之间的相关性函数为其中,ξ是一个预设的相关性下界,初始化第d次迭代的相关性下界ξd=1、第d次迭代的相关性函数/>
步骤2.5.2、通过式(5)得到第d+1次迭代的第i个解混向量wi,d+1
式(5)中,g'd为第d次迭代的相关性函数gd的一阶导数;
步骤2.5.3、对第d+1次迭代的第i个解混向量wi,d+1进行标准化:将wi,d+1/||wi,d+1||2赋予wi,d+1,其中||wi,d+1||2是wi,d+1的2-范数;
步骤2.5.4、根据式(6)、式(7)和式(8)分别更新其他参数:
μd+1=max{0,μddgd} (6)
δd+1=α×δd (7)
ξd+1=β×ξd (8)
式(6)、式(7)和式(8)中,μd+1表示第d+1次迭代的拉格朗日乘子,δd+1表示第d+1次迭代的惩罚因子,ξd+1表示第d+1次迭代的相关性下界;
步骤2.5.5、更新第d+1次迭代的第i个运动单位源信号更新第d+1次迭代的相关性函数/>
步骤2.5.6、若表示预设的迭代计算次数上限,则令d+1赋值给d,返回步骤2.5.2;否则,进行步骤2.5.7;
步骤2.5.7、令第i个运动单位源信号yi=yi,d+1,采用步骤2.3到步骤2.4的过程获取新的发放序列r'i;计算ri和r'i的的互相关系数ci
若ci>C,则表示更新后的运动单位源信号yi可靠,并将新的序列r'i和其对应的解混向量wi,d+1分别放入运动单位发放集合γ和解混向量集合ε中,其中C表示预设的用来衡量相关性的相关度界限;否则将新的序列r'i和其对应的解混向量wi,d+1直接删除;返回步骤2.5.1对源信号第i+1个运动单位源信号yi+1继续检测更新,直至最后一个运动单位源信号完成更新;
步骤2.6、遍历运动单位发放集合γ并分别计算任意两个不同的发放序列之间的相关系数,若大于C,则表示两个发放序列重复,并随机选取一个从运动单位发放集合γ和解混向量集合ε中删除;否则,两个发放序列均保留;经此操作后,可得到共有q个发放序列的运动单位发放集合γ和共有q个解混向量的解混向量集合ε;
步骤2.7、将运动单位发放集合γ中q个发放序列表示为R=[r1,r2,…,rq]T并计算其在最小二乘意义下的最优波形从而更新残差信号/>
步骤2.8、将更新后的残差信号代入步骤2.2-步骤2.7的过程进行处理,并判断前后两次处理得到的集合γ和ε中是否有新增的结果,若有,则继续返回至步骤2.2顺序执行,否则停止算法,并最终得到由离线数据集经过分解得到后的运动单位发放集合γ和包含Q个解混向量的解混向量集合ε;由Q个解混向量组成解混矩阵其中,/>表示解混矩阵/>中的第i个解混向量;
本实施例中,利用了较为代表性的传统的离线分解算法——PFP算法,来识别离线数据集中包含的运动单位。其中,步骤2.1的矩阵扩展部分,延迟因子K选择为5,这与实验的采样率有关;步骤2.2为了尽可能迭代得到较多的解混向量,一般n会取得较大,示例性的,可以设置为n=200;步骤2.2.3中的非多项式函数通常设置为G(x)=log(cosh(x));步骤2.2.5中的为求解每个解混向量所需要的迭代次数上限,通常设置为40;步骤2.3的大津法通过遍历区间内所有点寻求使得类间方差最大的阈值,为了使得遍历的结果尽可能可靠,区间需要划分的尽可能多,因此M需要尽可能大一点,示例性的,可以选择M=200;步骤2.5中带时域约束的FastICA算法是为了检测分解的发放序列的可靠性,通常认为,一个可靠的发放序列收敛至尽可能相似于待测序列的输出,这样判定了待测发放序列的可靠性,同时通过更新解混向量更正了可能的错误或遗漏的发放时刻;步骤2.5.1中的惩罚因子通常初始设置为一个接近于0的小数,示例性的,可以设置为10-5;步骤2.5.4中,为保证收敛,α通常设置为靠近1但比1大的数值,β通常设置为靠近1但比1大的数值,示例性的,可以设置为α=1.02,β=0.97;步骤2.5.6中的迭代次数上限/>通常设置为40,同步骤2.2.5;步骤2.5.7和步骤2.6中的C是两个序列之间的相关度界限,通常设置为0.4,若两个发放序列的互相关系数超过了0.4,则认为这两个发放序列来自于同一运动单位,基于此原理,步骤2.5.7用来判断更新前后发放序列是否发生了质的变化,步骤2.6用来去除集合中重复的发放序列;步骤2.7的剥离策略是为了把已得到的运动单位从原肌电信号中去除,避免梯度下降过程带来的局部收敛,提高分解的产量;
步骤3、实时采集肌电信号,并将获取的肌电数据流按照时间轴分为一系列窗长为L,步长为D的分析窗;在每个分析窗内,根据式(1)得到扩展后的肌电数据X,并利用卷积模型 求解运动单位源信号所组成的矩阵S=[s1,s2,…,si,…,sQ]T,si代表其中第i个运动单位源信号;
本实施例中,实时采集肌电信号的实验方案和设备与离线数据集的采集过程保持一致,这里就不再赘述;为了保证窗长设置的足够合理,既能满足在线数据处理的实时性需求,又可以保证分解结果的可靠,示例性的,可以设置为L=5s,D=1s;
步骤4、采用自适应的阈值选取策略提取分析窗内第i个运动单位源信号si对应的发放序列并保存;肌电信号的采集过程难免会受到环境噪声的影响,由步骤4中卷积模型直接求解的结果难免会有部分噪声序列作为干扰,因此需要通过自适应调节阈值的过程,将信号幅度普遍偏小的噪声序列尽可能去除掉,保证提取的发放序列的准确性。本发明通过执行步骤4.1至步骤4.6来对窗内第i个运动单位源信号si提取发放序列,保证尽可能不受噪声干扰:
步骤4.1、按照步骤2.3中的大津法计算第i个运动单位源信号si的初始阈值th0
步骤4.2、将si中所有大于th0的波峰的幅值置于集合β中,从而得到共有l个波峰幅值的集合β,记作β=[peak1,peak2,…,peakk,…,peakl];其中,peakk代表第k个波峰的幅值;
步骤4.3、初始化k=1;
步骤4.4、令中间变量th=peakk,将si中大于阈值th的波峰峰值位置置为1,其余置为0,从而得到对应的运动单位发放序列bk
步骤4.5、设置向量a,用来存储运动单位发放序列bk中每相邻两数值为1的位置之间的位置差;利用式(9)计算第k个发放变异系数cosk
cosk=std(a)/mean(a) (9)
式(9)中,std(a)表示向量a的标准差,mean(a)表示向量a的平均值;
步骤4.6、若k<l,则令k+1赋值给k,返回步骤4.4;否则找到最大的cosk对应的下标kmin,并将下标kmin所对应的发放序列作为第i个运动单位源信号si对应的发放序列;
本实验中,采用了间隔变异系数作为运动单位序列质量的衡量指标,通常认为正常的运动单位是以较为稳定的时间间隔募集和发放,故其对应的发放变异系数应当尽可能的小,因而这里寻找到发放变异系数最小的结果作为最终的结果;
步骤5、每当新的肌电数据输入时,按照步骤3和步骤4的过程进行处理和保存;
当没有新的数据输入时,利用式(10)将所保存的P个分析窗中的发放序列按照时间顺序进行拼接,得到总的分解结果Spike并作为最终的结果:
Spike=[Win1,Win'2,Win'3,…,Win'i,…,Win'P] (10)
式(10)中,Wini为第i个窗的分解的所有发放序列;Win'i表示Wini的结尾长度为D的数据;
如前所述,实时采集的肌电数据以窗长L、步长D划分分析窗,故两个相邻的分析窗中有L-D的重叠段,因此在连接结果的时候需要去除这部分重复段,仅取每个分析窗结尾为D长度的部分进行连接;此外,在肌电分解领域,通常认为一个解混向量对应一个运动单位,故在本发明直接采用矩阵拼接的形式连接最终的结果。
为了量化评估本发明的效果,通过8位不同受试者的实验数据将本发明方法与传统的离线分解方法(即步骤2内容)相比较,所有的数据处理过程在软件MATLAB 2020a上运行。
对比实验中,在线处理了20s的表面肌电信号,采集过程中共执行了4次的拇指收缩任务。在线过程的20s的信号除了完成了在线识别运动单位发放序列,还应用步骤二所详述的离线分解算法进行处理(将数据分为4个5s的数据段分别分解,每个数据段对应于1次收缩任务)。图3是其中一位受试者一个5s数据段的分解结果示意图。从图上可以看出,在线处理的结果和离线的结果很好地匹配上了,存在较少的错误或者遗漏的发放。图4展示了8位受试者的在线处理结果。可以发现,离线分解算法共获取到了10.31±1.46个运动单位,在线处理过程共获取到了11.38±2.50个运动单位,共有6.56±1.83个运动单位被成功匹配上了,这些运动单位称之为离线和在线过程共有的运动单位。在这些共有的运动单位中,匹配度(Matching Rate,MR)平均为(92.65±3.77)%,错误发现率(False DiscoveryRate,FDR)为0.066±0.030,假阴性率(False Negative Rate,FNR)为0.074±0.044。可以发现,本发明在不同的受试者上均能与离线的分解算法达到90%以上的匹配度,错误率维持在一个较小的状态,且运行速度相较于离线分解得到大幅的提升,平均每个分析窗的处理时间从62.44秒缩减到0.53秒,证明了所提出的在线识别运动单位方法的可行性和可靠性。
综上所述,本发明能够解决肌电分解难以在线运行的问题,将离线分解的结果作为先验知识指导在线的处理过程,并且通过自适应的阈值选取策略,对结果进行修正和完善,以提升整体的精度。与传统的离线分解算法相比,本发明能够更加快速地获取到运动单位的活动信息,在肌电控制方面有着更广泛的应用面,对运动控制和医疗康复等领域具有重要意义。
通过以上的实施方式的描述,本领域的技术人员可以清楚地了解到上述实施例并可以通过相应的软件和代码实现。上述实施例的技术方案可以以软件产品的形式体现出来,该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM,U盘,移动硬盘等)中,包括若干指令用以使得一台计算机设备(可以是个人计算机,服务器,或者网络设备等)执行本发明各个实施例的方法。

Claims (2)

1.一种基于表面肌电分解的运动单位识别方法,其特征是按如下步骤进行:
步骤1、利用肌电测量设备和M通道高密度电极在待测肌肉处预采集t时刻肌肉收缩任务的表面肌电信号数据,记为x(t)=[x1(t),x2(t),…,xI(t),…,xM(t)]T,其中,xI(t)表示t时刻第I个通道的肌电数据;从而得到时长为T1的表面肌电信号数据并作为离线数据集;
步骤2、采用逐步变量剥离算法分解所述离线数据集的肌电信号:
步骤2.1、利用式(1)对所述离线数据集进行扩展,得到扩展信号
式(1)中,K是延迟因子;
定义并初始化运动单位发放集合γ及其对应的解混向量集合ε为空集;
定义并初始化残差信号
步骤2.2、执行FastICA算法,设置算法执行的总次数为n:
步骤2.2.1、初始化i=1;
步骤2.2.2、初始化j=0;
步骤2.2.3、定义并初始化残差信号中第i个解混向量在第j次的计算结果wi,j为一个均值为0,方差为1的高斯信号;
步骤2.2.4、利用式(2)计算第i个解混向量在第j+1次的计算结果wi,j+1
式(2)中,G'表示非多项式函数G的一阶导数;G”表示G的二阶导数;E是数学期望;代表第i个解混向量在第j次的计算结果wi,j的转置;
步骤2.2.5、如果 表示预设的迭代计算次数上限,则令j+1赋值给j,返回至步骤2.2.3;若i<n且/>则得到第i个解混向量wi=wi,j,则令i+1赋值给i,返回至步骤2.2.2;若i≥n且/>则得到n个解混向量及其组成的解混矩阵W=[w1,w2,…,wi,…,wn]T,并利用式(3)得到正交化后的解混矩阵W':
W'=(WWT)-1/2W (3)
步骤2.2.6、计算运动单位源信号且Y=[y1,y2,…,yi,…,yn]T,其中,yi代表第i个解混向量wi所对应的运动单位源信号;
步骤2.3:采用图像处理的大津法获取第i个运动单位源信号yi的阈值thi
步骤2.3.1、将区间(0,max(yi))划分为M个等分,得到[thi,1,thi,2,thi,3,…,thi,m,…,thi,M];其中,max(yi)表示运动单位源信号yi的最大值,thi,m表示第m个小区间的最大值;
步骤2.3.2、初始化h=1;
步骤2.3.3、令thi=thi,h;求出yi中大于thi的采样点个数N0,h和小于thi的采样点个数N1,h;并求解个数N0,h对应的灰度μ0,h作为yi中所有大于thi的采样点的平均值;求解个数N1,h对应的灰度μ1,h作为yi中所有小于thi的采样点的平均值;
步骤2.3.4、根据式(4)计算得到类间方差Fi,h
Fi,h=N0,h×N1,h×(μ0,h1,h)2 (4)
步骤2.3.5、若h<M,则令h+1赋值给h,返回至步骤2.3.3;否则进入步骤2.3.6;
步骤2.3.6、找到最大的Fi,h对应的下标hmax,并将下标hmax所对应的thi,hmax作为第i个运动单位源信号yi的阈值thi
步骤2.4、将yi中大于阈值thi的波峰峰值位置置为1,其余置为0,从而得到对应的运动单位发放序列ri
步骤2.5、通过带时域约束的FastICA算法检测分解结果的可靠性:
步骤2.5.1、定义迭代计算的次数为d,并初始化d=1;定义并初始化第d次迭代的拉格朗日乘子μd=0、第d次迭代的惩罚因子δd为随机小数;初始化第d次迭代的第i个运动单位源信号yi,d=yi、第d次迭代的第i个解混向量wi,d=wi;定义第i个运动单位源信号yi和待测发放序列ri之间的相关性函数为其中,ξ是一个预设的相关性下界,初始化第d次迭代的相关性下界ξd=1、第d次迭代的相关性函数/>
步骤2.5.2、通过式(5)得到第d+1次迭代的第i个解混向量wi,d+1
式(5)中,g'd为第d次迭代的相关性函数gd的一阶导数;
步骤2.5.3、对第d+1次迭代的第i个解混向量wi,d+1进行标准化:将wi,d+1/||wi,d+1||2赋予wi,d+1,其中||wi,d+1||2是wi,d+1的2-范数;
步骤2.5.4、根据式(6)、式(7)和式(8)分别更新其他参数:
μd+1=max{0,μddgd} (6)
δd+1=α×δd (7)
ξd+1=β×ξd (8)
式(6)、式(7)和式(8)中,μd+1表示第d+1次迭代的拉格朗日乘子,δd+1表示第d+1次迭代的惩罚因子,ξd+1表示第d+1次迭代的相关性下界;
步骤2.5.5、更新第d+1次迭代的第i个运动单位源信号更新第d+1次迭代的相关性函数/>
步骤2.5.6、若 表示预设的迭代计算次数上限,则令d+1赋值给d,返回步骤2.5.2;否则,进行步骤2.5.7;
步骤2.5.7、令第i个运动单位源信号yi=yi,d+1,采用步骤2.3到步骤2.4的过程获取新的发放序列r'i;计算ri和r'i的的互相关系数ci
若ci>C,则表示更新后的运动单位源信号yi可靠,并将新的序列r'i和其对应的解混向量wi,d+1分别放入运动单位发放集合γ和解混向量集合ε中,其中,C表示预设的用来衡量相关性的相关度界限;否则,将新的序列r'i和其对应的解混向量wi,d+1直接删除;返回步骤2.5.1对源信号第i+1个运动单位源信号yi+1继续检测更新,直至最后一个运动单位源信号完成更新;
步骤2.6、遍历运动单位发放集合γ并分别计算任意两个不同的发放序列之间的相关系数,若大于C,则表示两个发放序列重复,并随机选取一个从运动单位发放集合γ和解混向量集合ε中删除;否则,两个发放序列均保留;经此操作后,可得到共有q个发放序列的运动单位发放集合γ和共有q个解混向量的解混向量集合ε;
步骤2.7、将运动单位发放集合γ中q个发放序列表示为R=[r1,r2,…,rq]T并计算其在最小二乘意义下的最优波形从而更新残差信号/>
步骤2.8、将更新后的残差信号代入步骤2.2-步骤2.7的过程进行处理,并判断前后两次处理得到的集合γ和ε中是否有新增的结果,若有,则继续返回至步骤2.2顺序执行,否则停止算法,并最终得到由离线数据集经过分解得到后的运动单位发放集合γ和包含Q个解混向量的解混向量集合ε;由Q个解混向量组成解混矩阵/>其中,/>表示解混矩阵/>中的第i个解混向量;
步骤3、实时采集肌电信号,并将获取的肌电数据流按照时间轴分为一系列窗长为L,步长为D的分析窗;在每个分析窗内,根据式(1)得到扩展后的肌电数据X,并利用卷积模型 求解运动单位源信号所组成的矩阵S=[s1,s2,…,si,…,sQ]T,si代表其中第i个运动单位源信号;
步骤4、采用自适应的阈值选取策略提取分析窗内第i个运动单位源信号si对应的发放序列并保存;
步骤5、每当新的肌电数据输入时,按照步骤3和步骤4的过程进行处理和保存;
当没有新的数据输入时,利用式(9)将所保存的P个分析窗中的发放序列按照时间顺序进行拼接,得到总的分解结果Spike并作为最终的结果:
Spike=[Win1,Win′2,Win′3,…,Win′i,…,Win′P] (9)
式(9)中,Wini为第i个窗的分解的所有发放序列;Win'i表示Wini的结尾长度为D的数据。
2.根据权利要求1所述的基于表面肌电分解的运动单位识别方法,其特征在于,所述步骤4中的自适应阈值选择策略是按如下过程进行:
步骤4.1、按照步骤2.3中的大津法计算第i个运动单位源信号si的初始阈值th0
步骤4.2、将si中所有大于th0的波峰的幅值置于集合β中,从而得到共有l个波峰幅值的集合β,记作β=[peak1,peak2,…,peakk,…,peakl];其中,peakk代表第k个波峰的幅值;
步骤4.3、初始化k=1;
步骤4.4、令中间变量th=peakk,将si中大于阈值th的波峰峰值位置置为1,其余置为0,从而得到对应的运动单位发放序列bk
步骤4.5、设置向量a,用来存储运动单位发放序列bk中每相邻两数值为1的位置之间的位置差;利用式(10)计算第k个发放变异系数cosk
cosk=std(a)/mean(a) (10)
式(10)中,std(a)表示向量a的标准差,mean(a)表示向量a的平均值;
步骤4.6、若k<l,则令k+1赋值给k,返回步骤4.4;否则找到最大的cosk对应的下标kmin,并将下标kmin所对应的发放序列作为第i个运动单位源信号si对应的发放序列。
CN202210417737.8A 2022-04-20 2022-04-20 一种基于表面肌电分解的运动单位识别方法 Active CN114767132B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210417737.8A CN114767132B (zh) 2022-04-20 2022-04-20 一种基于表面肌电分解的运动单位识别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210417737.8A CN114767132B (zh) 2022-04-20 2022-04-20 一种基于表面肌电分解的运动单位识别方法

Publications (2)

Publication Number Publication Date
CN114767132A CN114767132A (zh) 2022-07-22
CN114767132B true CN114767132B (zh) 2024-05-07

Family

ID=82430460

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210417737.8A Active CN114767132B (zh) 2022-04-20 2022-04-20 一种基于表面肌电分解的运动单位识别方法

Country Status (1)

Country Link
CN (1) CN114767132B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117017323B (zh) * 2023-09-14 2024-03-29 中国科学技术大学 基于盲源分离的高密度表面膈肌肌电采集与预处理方法
CN117562560B (zh) * 2023-11-20 2024-07-02 北京宜善医学科技有限公司 一种康复训练中的运动效果评估方法、装置及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109359619A (zh) * 2018-10-31 2019-02-19 浙江工业大学之江学院 一种基于卷积盲源分离的高密度表面肌电信号分解方法
CN110464348A (zh) * 2019-07-10 2019-11-19 深圳市智能机器人研究院 基于肌电信号的下肢关节连续运动量识别方法及系统
CN113536911A (zh) * 2021-06-08 2021-10-22 西安交通大学 一种基于双线程的肌电在线实时分解方法
CN113827257A (zh) * 2021-07-06 2021-12-24 杭州电子科技大学 基于FastICA和收缩力的表面肌电分解方法
CN116089851A (zh) * 2022-11-07 2023-05-09 中国科学技术大学 一种自适应更新的在线运动单位识别方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111207875B (zh) * 2020-02-25 2021-06-25 青岛理工大学 基于多粒度并联cnn模型的肌电信号-扭矩匹配方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109359619A (zh) * 2018-10-31 2019-02-19 浙江工业大学之江学院 一种基于卷积盲源分离的高密度表面肌电信号分解方法
CN110464348A (zh) * 2019-07-10 2019-11-19 深圳市智能机器人研究院 基于肌电信号的下肢关节连续运动量识别方法及系统
CN113536911A (zh) * 2021-06-08 2021-10-22 西安交通大学 一种基于双线程的肌电在线实时分解方法
CN113827257A (zh) * 2021-07-06 2021-12-24 杭州电子科技大学 基于FastICA和收缩力的表面肌电分解方法
CN116089851A (zh) * 2022-11-07 2023-05-09 中国科学技术大学 一种自适应更新的在线运动单位识别方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于SEONS算法的表面肌电信号分解方法研究;李强;杨基海;陈香;张旭;;航天医学与医学工程;20070430(第02期);全文 *
基于小波分析的表面肌电信号运动单位动作电位模式研究;李强;杨基海;;生物医学工程学杂志;20100825(第04期);全文 *

Also Published As

Publication number Publication date
CN114767132A (zh) 2022-07-22

Similar Documents

Publication Publication Date Title
CN114767132B (zh) 一种基于表面肌电分解的运动单位识别方法
CN109948647B (zh) 一种基于深度残差网络的心电图分类方法及系统
Clarke et al. Deep learning for robust decomposition of high-density surface EMG signals
CN110786850B (zh) 基于多特征稀疏表示的心电信号身份识别方法及系统
CN112861604B (zh) 一种与用户无关的肌电动作识别与控制方法
CN111657941B (zh) 基于肌肉核心激活区域的电极校正及肌电模式识别方法
CN110598676B (zh) 基于置信度得分模型的深度学习手势肌电信号识别方法
Dai et al. Independent component analysis based algorithms for high-density electromyogram decomposition: Systematic evaluation through simulation
Scano et al. Mixed matrix factorization: A novel algorithm for the extraction of kinematic-muscular synergies
CN116089851A (zh) 一种自适应更新的在线运动单位识别方法
CN111329476B (zh) 一种基于微观神经驱动信息进行肌力估计的方法及装置
CN109598219B (zh) 一种用于鲁棒肌电控制的自适应电极配准方法
CN113397571A (zh) 一种基于先验模板的肌电运动单元分解方法
CN115607164B (zh) 心电特征波分割方法、系统、装置及可读存储介质
CN108509823A (zh) Qrs波群的检测方法及装置
Márquez-Figueroa et al. Optimal extraction of EMG signal envelope and artifacts removal assuming colored measurement noise
CN116671919B (zh) 一种基于可穿戴设备的情绪检测提醒方法
CN110974223A (zh) 基于改进ksvd算法的表面肌电信号压缩重构方法
Zhao et al. Online decomposition of surface electromyogram into individual motor unit activities using progressive FastICA peel-off
Wu et al. A new EMG decomposition framework for upper limb prosthetic systems
CN111110268B (zh) 基于随机向量功能连接网络技术的人体肌音信号预测方法
CN116712099A (zh) 基于多模态数据的胎心状态检测方法、电子设备、存储介质
Nawab et al. Improved resolution of pulse superpositions in a knowledge-based system EMG decomposition
CN112932508B (zh) 一种基于手臂肌电网络的手指活动识别系统
Xue et al. SEMG based intention identification of complex hand motion using nonlinear time series analysis

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