一种基于双线程的肌电在线实时分解方法与流程

1.本发明属于肌电实时分解技术领域，具体涉及一种基于双线程的肌电在线实时分解方法。


背景技术：

2.肌电分解(emg decomposition)是根据肌电信号形成的机理，逆向求解还原出构成肌电信号的各个运动单元动作电位(motor unit action potential，muap)序列或放电事件序列的过程。肌电分解所得各个运动单元放电事件信息可用于估计肌肉收缩强度，进而解码人体复杂运动意图并控制假肢等外部辅助运动设备。与传统基于肌电信号时频域特征的肌肉收缩强度估计方法相比，利用运动单元放电事件信息估计肌肉收缩强度，不仅更具有生理学上的合理性，同时避免了肌电信号形成过程中各种干扰因素对肌肉收缩强度估计的影响。
3.离线状态下肌电分解技术发展较早，早期的肌电分解算法以carlo j.de luca及其团队提出的muap模板匹配法为代表，尽管获得了广泛应用，但是其存在分解过程耗时长和分解所得运动单元个数少等问题。此后，盲源分离技术的出现为肌电分解提供了新的思路，主要包括卷积核补偿法(convolution kernel compensation，ckc)和独立成分分析法(independent component analysis，ica)。ckc放弃了对muap波形(卷积核)的直接估计，转而直接重建各运动单元的放电事件序列，但不同muap重叠增加会导致其分解性能下降，因此ckc法较适用于低肌肉收缩强度情形。ica法假设各个运动单元放电活动的独立性，依据对源成分独立性测量的不同，又可细分为基于fastica(负熵)、informaxica(互信息)和robustica(峰度)的肌电分解方法。此外，一些改进算法也被相继提出，如结合k
‑
means聚类分析和ckc的分解算法，结合fastica和ckc的分解算法，以及增加运动单元分解个数的peel
‑
off
‑
fastica算法等。上述算法仅适用于离线状态下的肌电分解，运动单元放电事件提取过程耗时长，无法满足运动意图实时解码和假肢等外部设备实时控制的需求。
4.为了实现在线实时肌电分解，holobar教授提出一种基于ckc法的肌电在线分解算法，但分解输出运动单元个数较少且无法满足一般外部设备实时控制要求。
5.目前，肌电实时分解过程一般可分为两个步骤，第一个步骤识别运动单元，在初始化离线阶段完成，主要用于提取各个运动单元对应的特征信息。特征信息因分解方法不同而有所差异，如模板匹配法中指的是各个运动单元的muap波形，而盲源分离法中指运动单元对应的分离向量。第二个步骤为运动单元放电事件的实时检测，主要利用初始化阶段提取的各个运动单元特征信息在线实时提取运动单元放电事件。尽管国内外研究人员开展了肌电实时分解技术的初步研究工作，但是现有方法均建立在运动单元募集以及运动单元特征信息稳态的基础上。这一假设虽然简化了肌电实时分解问题，但是其与实际肢体随意运动下运动单元募集以及特征信息具有非稳态特征不符。运动单元募集的非稳态性表现在，正常生理状态下，肌肉持续收缩时运动单元会出现交替放电现象，即开始被募集的运动单元停止放电，然后募集开始未放电的运动单元，这一神经调节机制能够有效缓解肌疲劳。运
动单元特征信息的非稳态性主要由运动单元muap波形变异造成，如运动单元持续放电时通常伴随其muap波形的变异，动态收缩及肢体体位变化等会导致肌纤维与皮肤表面相对位置变化，进而同样引起运动单元的muap波形变异。在现有肌电实时分解方法中，运动单元的识别过程仅在初始化阶段进行，因此运动单元募集的非稳态性会导致后来被募集的运动单元放电事件无法被检测到。此外，现有肌电实时分解方法通常在初始化阶段提取运动单元特征信息后认为其保持不变，而实际情况下特征信息的非稳态性会造成运动单元放电事件检测准确性显著降低。这些共同造成对肌肉收缩强度估计的偏差，进而影响运动意图解码准确性。


技术实现要素：

6.本发明的目的是针对现有肌电实时分解技术的局限性，提出一种基于双线程的肌电在线实时分解方法，能够在线自动识别新募集运动单元，并能够感知运动单元特征信息的变化，自动调整运动单元特征信息，从而使得所提自适应肌电实时分解技术能够更全面、更准确地检测肌肉运动单元集的放电活动，提高对肌肉收缩强度的估计的准确性，最终提高运动意图解码的准确性，满足假肢等外部辅助运动设备精细运动控制的需求。
7.本发明采用的技术方案是：
8.一种基于双线程的肌电在线实时分解方法，肌电实时分解方法包括主线程和副线程两个平行线程，具体包括以下步骤：
9.步骤1：副线程用于运动单元在线识别，即运动单元特征信息在线提取与更新；
10.步骤2：主线程用于运动单元放电事件实时检测，即利用步骤1提取的各个运动单元特征信息在线实时提取运动单元放电事件。
11.优选的，在步骤1中，运动单元特征信息在线提取与更新主要包括以下步骤：
12.步骤101：精炼已识别运动单元的分离向量，从而得到分离矩阵
13.步骤102：识别新运动单元，提取新识别运动单元特征信息，并将新识别运动单元分离矩阵b
new
与精炼后分离矩阵合并，从而更新运动单元分离矩阵为
14.优选的，在步骤101中，给定已识别运动单元的分离向量w0和muap波形，以及25秒肌电数据x，其分离向量w0的精炼过程为：
15.(1)对肌电数据x进行ckc迭代，获得运动单元新的分离向量w和放电事件序列t；
16.(2)若放电事件总次数小于100，则拒绝新的分离向量w，进行步骤(6)，否则，进行步骤(3)；
17.(3)利用放电事件序列t和肌电数据x，通过锁时平均法获得运动单元新的muap波形；
18.(4)计算运动单元新的muap波形和给定的muap波形的muap波形相似指数；
19.(5)若muap波形相似指数大于0.95，则使用新的分离向量w和muap波形更新给定运动单元的特征信息；否则，拒绝新的分离向量w和muap波形，并进行步骤(6)；
20.(6)结束给定运动单元分离向量精炼。
21.优选的，在步骤(1)中，ckc迭代过程如下：
22.a、将肌电数据x进行扩展与白化，得到肌电矩阵
23.b、初始化放电事件序列t0为全零向量，w＝w0；
24.c、计算运动单元对应源信号：
[0025][0026]
d、通过对源信号s进行峰值检测以及峰值二分类聚类分析，获得放电事件序列t；
[0027]
e、重新计算分离向量为其中t(t
j
)＝1，j代表放电总次数；
[0028]
f、若t不等于t0，用t代替t0，并回到步骤c；否则，进行步骤g；
[0029]
g、输出精炼后分离向量w和新的放电事件序列t。
[0030]
优选的，在步骤(2)中，设定最少放电次数100。
[0031]
优选的，在步骤(4)中，若新的muap波形和给定的muap波形分别为x
i
(t)和y
j
(t)，i,j＝1,2,...,m，t＝1,...,t，其中m为高密度肌电信号电极个数，t为muap波形时长，则muap波形相似指数计算如下：
[0032][0033]
其中，其中，式(2)右边第二项可看作各电极处muap波形差异的能量加权平均，s变化范围从0到1，0表明两组muap波形完全相反，1表示两组muap波形完全相同；常数c为尺度调节因子，可调节信号能量高、低电极处的权重分布，本发明中常数c设置为4。
[0034]
优选的，在步骤(5)中，当muap波形相似指数低于设定的阈值时，该阈值设置为0.95，说明ckc迭代过程使得迭代后的分离向量收敛到其他运动单元对应的局部最优解，则舍弃ckc迭代后的分离向量，并放弃对给定运动单元的本次精炼。
[0035]
优选的，在步骤102中，首先通过fastica算法得到一定数量运动单元的特征信息，然后利用muap波形相似指数判断新识别运动单元和已识别运动单元中存在重复的运动单元，并去除重复运动单元，假设还保留n2个运动单元，已识别运动单元集中包含n1个运动单元，则识别去除重复运动单元获得运动单元分离矩阵b
new
的过程如下：
[0036]
a、将n2个保留的运动单元全部初始化为有效，且非已识别运动单元的副本；令i用于标记n1个已识别运动单元，j用于标记n2个保留的运动单元；
[0037]
b、从n1个已识别运动单元中选择运动单元i，从n2个保留的运动单元中选择运动单元j；
[0038]
c、若运动单元j被标记为无效或者某个运动单元的副本，则进行步骤h；否则进行步骤d；
[0039]
d、利用求取运动单元i在副线程滑动窗内的源信号，并获得运动单元i的放电事件序列，计算运动单元i和j同步放电占总放电次数百分比；
[0040]
e、若同步放电占比超过80％，进行步骤f；否则，标记运动单元j为真正的新识别运动单元，并进行步骤h；
[0041]
f、计算运动单元i和j间的muap波形相似指数；
[0042]
g、若muap波形相似指数大于0.95，则标记运动单元j为运动单元i的副本；否则，标记运动单元j为无效运动单元；
[0043]
h、若n1个已识别运动单元中的每个运动单元和n2个保留运动单元中的每个运动单元均进行了步骤c
‑
g，则结束此过程；否则，返回步骤b；
[0044]
经过上述循环后，n2个运动单元被标记为如下三类，即无效运动单元、运动单元集中某一运动单元的副本、真正的新识别运动单元；首先，直接去除无效运动单元，然后，如果多个运动单元是同一运动单元的副本，则选择源信号峰值二分类效果最好的一个替换中对应的运动单元，并去除所有副本，最后，b
new
中仅保留真正的新识别运动单元。
[0045]
优选的，在步骤2中，运动单元放电事件在线实时检测过程如下：
[0046]
步骤201：首先对1秒长的肌电数据进行扩展与白化预处理；
[0047]
步骤202：然后用更新后运动单元分离矩阵b中包含的各运动单元分离向量左乘肌电矩阵即得到各个运动单元对应的源信号s；
[0048]
步骤203：对源信号s进行峰值检测和二分类聚类分析，得到各个运动单元的放电事件；计算各个源信号峰值二分类的加权平均silhouette距离，当其小于一定阈值，则认为所检测放电事件为假阳性事件并舍弃，否则认为真阳性事件并保留，并保留最后0.2秒内的放电事件。
[0049]
本发明的有益效果：本发明提出的一种基于双线程的肌电在线实时分解方法具有以下明显效果：
[0050]
(1)可有效应对肌肉持续收缩时运动单元交替放电带来的新募集运动单元识别问题，从而减少遗漏某些运动单元放电信息的可能性，保证肌肉收缩强度估计的准确性。
[0051]
(2)可有效应对肢体复杂随意运动时各种因素导致的运动单元特征信息非稳态问题，从而提高运动单元放电信息检测的准确性，保证肌肉收缩强度估计的准确性。
附图说明
[0052]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0053]
图1所示为本发明的一种基于双线程的肌电在线实时分解方法的流程图。
[0054]
图2所示为本发明的仿真实验结果图。
具体实施方式
[0055]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
[0056]
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护
的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0057]
本发明具体提供了一种基于双线程的肌电在线实时分解方法，其中，运动单元识别和放电事件检测方法主要基于fastica算法，在fastica算法中，原始高密度肌电数据x首先经过扩展与白化处理，获得新的肌电矩阵然后，对肌电矩阵进行fastica迭代，每次迭代收敛到局部最优解后得到一个运动单元的相关信息，包括分离向量w，源信号s和放电事件序列t，其中源信号s由式(1)得到
[0058][0059]
放电事件序列t由源信号s依次经过峰值检测和k
‑
means聚类(二分类)分析得到。由于肌电预处理中对原始肌电数据进行了扩展处理，因此，不同次迭代可能收敛到同一运动单元对应的最优解。为此，可利用放电事件同步性找出同属一个运动单元的特征信息，并去除多余的特征信息，使得每个已识别运动单元均是唯一的，且仅保留一组特征信息。在本实施例中，当同步放电事件占比总放电事件高于85％时，认为两个已识别运动单元实际同属一个运动单元。此外，有的迭代过程可能收敛至运动伪迹，其突出特征是在较长的一段时间内(本实施例中为25秒)，仅检测到1
‑
2个孤立的放电事件，将此类“运动单元”去除。有些迭代所得源信号峰值二分类效果较差，导致放电事件检测的准确性很低，同样将此类“运动单元”去除，其中二分类效果用源信号峰值二分类的加权平均silhouette距离进行评价。在本实施例中，重复运动单元，运动伪迹对应“运动单元”以及二分类效果差的“运动单元”去除过程统称为无效运动单元去除。最后，利用原始高密度肌电数据x以及放电事件序列t，通过锁时平均法获得各个保留运动单元的muap波形。在本实施例中，运动单元的分离向量以及muap波形称为运动单元的特征信息。将运动单元分离向量用于新获取的肌电数据并进行峰值检测，聚类分析和假阳性放电判别，即可获取新肌电数据对应各个运动单元的放电事件。所有已识别运动单元的分离向量合并构成已知运动单元的分离矩阵，即b＝[w1,w2,
…
w
n
]。
[0060]
通过上述分析，本实施例中，肌电实时分解方法包括主线程和副线程两个平行线程，如图1所示，具体包括以下步骤：
[0061]
步骤1：副线程用于运动单元在线识别，即运动单元特征信息在线提取与更新；
[0062]
步骤2：主线程用于运动单元放电事件实时检测，即利用步骤1提取的各个运动单元特征信息在线实时提取运动单元放电事件。
[0063]
在步骤1中，运动单元特征信息在线提取与更新主要包括以下步骤：
[0064]
步骤101：精炼已识别运动单元的分离向量，从而得到分离矩阵
[0065]
步骤102：识别新运动单元，提取新识别运动单元特征信息，并将新识别运动单元分离矩阵b
new
与精炼后分离矩阵合并，从而更新运动单元分离矩阵为
[0066]
上述步骤101中，给定已识别运动单元的分离向量w0和muap波形，以及25秒肌电数据x，其分离向量w0的精炼过程为：
[0067]
(1)对肌电数据x进行ckc迭代，获得运动单元新的分离向量w和放电事件序列t；
[0068]
(2)若放电事件总次数小于100，则拒绝新的分离向量w，进行步骤(6)，否则，进行步骤(3)；
[0069]
(3)利用放电事件序列t和肌电数据x，通过锁时平均法获得运动单元新的muap波形；
[0070]
(4)计算运动单元新的muap波形和给定的muap波形的muap波形相似指数；
[0071]
(5)若muap波形相似指数大于0.95，则使用新的分离向量w和muap波形更新给定运动单元的特征信息；否则，拒绝新的分离向量w和muap波形，并进行步骤(6)；
[0072]
(6)结束给定运动单元分离向量精炼。
[0073]
步骤(1)中，ckc迭代过程如下：
[0074]
a、将肌电数据x进行扩展与白化，得到肌电矩阵
[0075]
b、初始化放电事件序列t0为全零向量，w＝w0；
[0076]
c、计算运动单元对应源信号：
[0077][0078]
d、通过对源信号s进行峰值检测以及峰值二分类聚类分析，获得放电事件序列t；
[0079]
e、重新计算分离向量为其中t(t
j
)＝1，j代表放电总次数；
[0080]
f、若t不等于t0，用t代替t0，并回到步骤c；否则，进行步骤g；
[0081]
g、输出精炼后分离向量w和新的放电事件序列t。
[0082]
在步骤(2)中，当给定运动单元放电总次数过少时，通过锁时平均法获得的muap波形可能不准确，因此本发明中设定了一个最少放电次数100，用于保证所提取muap波形的准确性。
[0083]
在步骤(3)锁时平均法步骤：
[0084]
①
以放电事件序列t中每次放电时刻为0时刻，前后各20ms截取肌电数据x，得到长度为40ms的数据段，数据段个数等于放电总次数。
[0085]
②
将截取得到的所有40ms数据段进行平均，得到运动单元的muap波形。
[0086]
在步骤(4)中，若新的muap波形和给定的muap波形分别为x
i
(t)和y
j
(t)，i,j＝1,2,...,m，t＝1,...,t，其中m为高密度肌电信号电极个数，t为muap波形时长，则muap波形相似指数计算如下：
[0087][0088]
其中，其中，式(2)右边第二项可看作各电极处muap波形差异的能量加权平均，s变化范围从0到1，0表明两组muap波形完全相反，1表示两组muap波形完全相同；常数c为尺度调节因子，可调节信号能量高、低电极处的权重分布，常数c设置为4。
[0089]
在步骤(5)中，当muap波形相似指数低于某一阈值(本实施例依据经验值将阈值设置为0.95)时，说明ckc迭代过程使得迭代后的分离向量收敛到其他运动单元对应的局部最
优解，因此需要舍弃ckc迭代后的分离向量，并放弃对给定运动单元的本次精炼。这一过程主要用于保证运动单元在每次副线程更新前后的一致性。
[0090]
步骤102中，首先通过fastica算法得到一定数量运动单元的特征信息，考虑本实施例每间隔180秒进行一次新运动单元识别，每次副线程进行80次fastica迭代，从而得到80个运动单元特征信息，将无效运动单元和放电次数少于100次的运动单元去除后，假设还保留n2个运动单元，已识别运动单元集中包含n1个运动单元，则识别去除重复运动单元获得运动单元分离矩阵b
new
的过程如下：
[0091]
a、将n2个保留的运动单元全部初始化为有效，且非已识别运动单元的副本；令i用于标记n1个已识别运动单元，j用于标记n2个保留的运动单元；
[0092]
b、从n1个已识别运动单元中选择运动单元i，从n2个保留的运动单元中选择运动单元j；
[0093]
c、若运动单元j被标记为无效或者某个运动单元的副本，则进行步骤h；否则进行步骤d；
[0094]
d、利用求取运动单元i在副线程滑动窗内的源信号，并获得运动单元i的放电事件序列，计算运动单元i和j同步放电占总放电次数百分比；
[0095]
e、若同步放电占比超过80％，进行步骤f；否则，标记运动单元j为真正的新识别运动单元，并进行步骤h；
[0096]
f、计算运动单元i和j间的muap波形相似指数；
[0097]
g、若muap波形相似指数大于0.95，则标记运动单元j为运动单元i的副本；否则，标记运动单元j为无效运动单元；
[0098]
h、若n1个已识别运动单元中的每个运动单元和n2个保留运动单元中的每个运动单元均进行了步骤c
‑
g，则结束此过程；否则，返回步骤b；
[0099]
经过上述循环后，n2个运动单元被标记为如下三类，即无效运动单元、运动单元集中某一运动单元的副本、真正的新识别运动单元；首先，直接去除无效运动单元，然后，如果多个运动单元是同一运动单元的副本，则选择源信号峰值二分类效果最好的一个替换中对应的运动单元，并去除所有副本，最后，b
new
中仅保留真正的新识别运动单元。
[0100]
由于副线程进行的运动单元特征信息在线提取与更新过程较为耗时，无法做到实时处理，因此副线程所用滑动窗窗长为25秒，而步长为180秒。即，在利用实时获取的肌电数据，每0.2秒提取一次各已知运动单元放电事件的同时，每隔180秒用25秒长的肌电数据进行一次运动单元分离矩阵的更新。
[0101]
在步骤2中，运动单元放电事件在线实时检测过程如下：
[0102]
步骤201：首先对1秒长的肌电数据进行扩展与白化预处理，1秒滑动窗的步长为0.2秒；
[0103]
步骤202：然后用更新后运动单元分离矩阵b中包含的各运动单元分离向量左乘肌电矩阵即得到各个运动单元对应的源信号s；
[0104]
步骤203：对源信号s进行峰值检测和二分类聚类分析，得到各个运动单元的放电事件；计算各个源信号峰值二分类的加权平均silhouette距离，当其小于一定阈值，则认为所检测放电事件为假阳性事件并舍弃，否则认为为真阳性事件并保留，并保留最后0.2秒内
的放电事件。
[0105]
为验证本实施例所提方法的有效性，利用仿真数据验证了所提方法在线识别新募集运动单元，以及更新运动单元特征信息用以提高放电事件实时检测准确率的能力。每段仿真数据时长30分钟，共生成10段数据。为模拟运动单元交替放电现象，利用一随机交替放电仿真模型生成共计80个运动单元的放电事件序列。为模拟运动单元特征信息非稳态性，每次放电事件对应的muap波形幅值随机变动。分别利用传统方法和本发明所提方法进行运动单元放电事件实时提取，其中前25秒数据用于离线初始化，副线程分离矩阵更新每隔3分钟进行一次。10段数据处理后的平均结果如图2所示。
[0106]
其中，图2中上图是每次分离矩阵更新后已识别运动单元个数的变化情况，可以看出识别运动单元个数由初始化时(第一分钟)的10个左右增加到最后(第28分钟)的34个左右，且增加速度随时间逐渐减少，说明本实施例所提方法能够保证所识别运动单元的唯一性，否则识别运动单元个数会持续增加直至超过总运动单元个数。图2中下图比较了本实施例所提方法和传统方法在放电事件检测准确率方面的差异，可以看出，由于本实施例所提方法定期对运动单元特征信息进行更新，改善了其提取放电事件的准确性，仿真数据表明，平均检测准确率从86％提高到91％。
[0107]
本发明能够在线自动识别新募集运动单元，并能够感知运动单元特征信息的变化，自动调整运动单元特征信息，从而使得所提自适应肌电实时分解技术能够更全面、更准确地检测肌肉运动单元集的放电活动，提高对肌肉收缩强度的估计的准确性，最终提高运动意图解码的准确性，满足假肢等外部辅助运动设备精细运动控制的需求。
[0108]
以上所述，仅用以说明本发明的技术方案而非限制，本领域普通技术人员对本发明的技术方案所做的其它修改或者等同替换，只要不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。