一种基于先验模板的肌电运动单元分解方法

1.本发明属于信号处理及模式识别领域，涉及一种肌电信号的分解方法，特别涉及一种基于先验模板将肌电信号分解为运动单元成分信息的方法。


背景技术：

2.肌肉收缩时，运动单元中的运动神经元开始兴奋，导致所有由运动神经元支配的肌纤维收缩，从而产生运动单元动作电位，多个运动单元发放的动作电位经过复杂的综合叠加过程最后在检测电极处形成肌电信号。肌电信号的分解实际上就是信号生成的逆过程，其根本目的就是分离出各个运动单元发放的动作电位，获取运动单元相关信息。而这些运动单元相关信息，如发放波形、间隔发放模式和募集规律等，可以提供有关对神经系统健康至关重要的肌肉疾病的临床诊断的决定性细节，有助于从根本上探索神经肌肉系统的生理机制，在肌肉疾病的临床诊断、肌肉组织之间运动协调关系的研究及中枢神经元控制机理的探索等方面具有重大价值。因此，研究肌电信号分解对反映运动神经肌肉系统控制机理和临床应用诊断有着重要的意义。
3.近二十年来，许多研究者致力于肌电信号的分解工作，提出了多种分解技术，取得了一定的突破性成果。根据这些成果分解基本思路的不同，肌电信号分解方法可以被总结分为系统辨识法和形态学分解法两类。系统辨识法主要将人体的神经肌肉系统视为多输入多输出系统，并采用盲源分离方法将肌电信号分解为符合运动单元动作电位分布特征的各个分量。不少学者基于此作出了重大贡献，zazula和holobar等人将盲反卷积和卷积核补偿相结合应用于肌电信号信号分解，并获得了可靠的结果。陈茂启等人构建了一个用于高密度肌电信号分解的渐进式fast ica剥离框架，并通过多次试验和仿真验证了其可靠性。然而，目前的系统辨识法都是基于高密度肌电信号，对单通道肌电信号存在一定局限性。
4.而形态学分解法是根据不同运动单元动作电位波形的固有形态特征，通过优化算法直接从肌电信号中提取运动单元动作电位序列。deluca团队自1982年以来一直从事肌电信号形态学分解工作，最近建立了基于形态分解的precision decomposition iii系统，并将分解结果与肌内肌电信号进行比较以评估分解可靠性。随后，利用改系统分解了75％
‑
100％的最大自主收缩下的肌电信号信号，准确率达到了92.5％。这些研究都为肌电信号分解奠定了坚实的基础，为进一步探索运动单元对运动的影响提供了便利。


技术实现要素：

5.针对现有的肌电信号分解方法很少有通用数学表达式，且分解过程繁冗复杂，在实际应用中缺乏灵活性的问题，本发明提出一种基于先验模板的肌电信号分解方法。
6.本发明首先，利用几个带有可变表达式的先验模板来拟合运动单元动作电位的波形；接着提出了一种迭代匹配策略，实现了运动单元的逐层剥离；随后引入尖峰触发平均算法来重新确定运动单元动作电位波形，使其更符合真实的肌电运动单元。最后，建立基于真实信号分解结果的评估指标，用于评价分解算法的性能。
7.本发明包括以下步骤：
8.步骤1：采集志愿者15％最大自主收缩力持续作用下与静息状态下的肌电信号。
9.步骤2：采用小波包去噪处理步骤1中得到的肌电信号，以获得较为纯净的肌电信号。
10.步骤3：对肌电信号进行运动单元分解，包括以下具体步骤：
11.步骤3.1：根据运动单元发放序列波形的形态学特征来确定了四个先验模板。选择n阶hermite
‑
rodriguez函数作为先验模板对运动单元动作电位波形进行模拟，运动单元动作电位波形大多呈双相或三相的结构，n阶hermite
‑
rodriguez函数与正态函数的n阶导数成正比，能够有效地模拟运动单元动作电位波形。
12.步骤3.2：利用静息状态下的肌电信号作为参考信号对15％最大自主收缩力持续作用下的肌电信号进行有效峰值检测：选择参考信号的均值作为检测阈值，然后基于该检测阈值检测出15％最大自主收缩持续作用力下的肌电信号的有效峰值(波峰和波谷)，并按其幅值的绝对值大小进行降序排列。
13.步骤3.3：在有效峰值附近截取信号段与先验模板进行匹配处理，其过程如下：
14.1.在确定有效峰值后，采用带有可变窗长的滑动窗口在各个有效尖峰附近进行截取，获得部分待分析的信号段作为备选运动单元动作电位。
15.2.采用动态时间规整算法计算每个备选运动单元动作电位与四个先验模板的距离d，分别计算了四个先验模板的平均振幅值mean
t
。
16.3.结合动态时间规整算法和先验模板的平均幅值来确定备选运动单元动作电位是否与先验模板匹配。具体的匹配决策规则为：如果距离d小于二分之一的mean
t
，则认为该备选运动单元动作电位与先验模板匹配；否则不匹配。
17.4.决策完成后，根据与先验模板匹配的备选运动单元动作电位，选取备选运动单元动作电位的开始时间作为发放时刻。
18.步骤3.4：采用脉冲触发平均技术进行运动单元动作电位波形的重构，过程如下：
19.将每个先验模板识别出的运动单元动作电位的发放时刻视为原始semg信号的触发时间点，依次触发得到多个触发波形，然后对触发波形进行平均即可得到对应运动单元动作电位的重构波形。这些重构波形便可视为由各个运动单元发放的真实运动单元动作电位波形。
20.步骤3.5：通过比较幅度因子和参考信号尖峰点平均值的大小判断循环迭代过程是否继续。在每次迭代剥离运动单元动作电位的过程中，若幅值高于参考信号尖峰点平均值，则视残差信号为更新后的肌电信号，并返回步骤3.2继续迭代；否则，循环迭代结束。
21.本发明提供的技术方案的有益效果是：利用该分解方法分解肌电信号的优点包括模板的简单的数学表达式，降低了复杂性，提高了算法的鲁棒性。同时，依照系统辨识模型构建了仿真肌电信号，评价提出分解方法的有效性，结果表明该分解方法的精度为90.9％，对低信噪比很敏感。最后根据真实肌电的分解结果，提出了两个评价指标，逆向对分解算法评估，结果显示评价指标的变化是完全符合预期的，再一次证明了分解的可行性。此外，最后的分析表明，在肌肉收缩期间，大型肌肉运动单元激活相比小型肌肉运动单元激活提供更多的信息，显著揭示了肌肉信息的内部生理机制和主要特征。
附图说明
22.图1为本发明的实施流程图；
23.图2为一阶和二阶hermite
‑
rodriguez函数模拟的先验模板图；
24.图3为15％最大自主收缩力作用下采集的肌电信号分解过程及结果图；
25.图4为不同运动单元数目组合的运动单元动作电位序列与原始肌电信号的相关系数变化图。
具体实施方式
26.本发明在形态学分解思想的基础上提出的一种单通道肌电信号分解方法。其主要技术构思是，根据过去的研究结果确定运动单元的几种固定先验模板，并提取实际肌电信号信号的尖峰信息来设置模板参数，提取有效尖峰。接着，截取有效尖峰附近的信号段进行模板匹配，确定运动单元发放时刻。然后，根据脉冲触发平均技术重新确定运动单元发放波形。随后，从肌电信号中删除已识别的运动单元动作电位获得残差信号，再将残差信号视为新的肌电信号重复上述步骤，不断迭代直至满足终止条件，得到最终的分解结果。最后，基于分解结果建立一种新的互相关系数指标来验证分解方法的有效性和可行性。
27.下面将结合附图详细介绍实施流程，如图1所示，本实施例包括以下几个步骤：
28.步骤1，采集志愿者15％最大自主收缩力持续作用力下与静息状态下的肌电信号，具体方式为：
29.肌电信号是通过delsys trigno
tm
无线肌电系统进行采集，采集采样频率为2000hz，受试者为4名健康的右利手男性，肌电采集通道均为屈指浅肌。实验具体流程如下，前2秒屏幕显示“准备”，当播放提示音“ding”时，受试者逐渐握紧握力传感器，动作时间为2秒，当握力强度达到指定的自主收缩力时，会播放提示音“di”，提醒受试者维持稳定力量5秒，当受试者听到“dong”的提示音后，2秒内松手释放，随后休息2分钟，即完成一次实验流程。实验全程都需要进行握力跟踪，以确保5秒持续力作用过程中肌肉能保持稳定收缩状态。记录5秒的肌电信号，并采用小波包去噪处理，以获得较为纯净的信号。
30.步骤2，高质量的肌电信号将会取得到更好的分解效果，其关键在于有效的消除信号中混叠的噪声信号，降低噪声影响，提高信噪比。因此，本发明采用小波包去噪算法对肌电信号进行预处理，主要思想为：通过小波包变换将含噪肌电信号分解成不同频段的分量，而肌电信号成分一般集中于较大的小波系数中，噪声成分则遍布于小波域内，再利用不同尺度的小波包系数重构即可。参考信号即指在静息状态下采集的肌电信号。
31.步骤3：对肌电信号进行运动单元分解，包括以下具体步骤：
32.步骤3.1，根据运动单元发放序列波形的形态学特征来确定四个先验模板：运动单元动作电位波形大多呈双相或三相的结构，而n阶hermite
‑
rodriguez函数与正态函数的n阶导数成正比，能够有效地模拟运动单元动作电位波形。n阶hermite
‑
rodriguez函数的基本表达式为：
[0033][0034]
其中，n＝1,2,...,n为hr函数的阶数，α表示与muap波形幅值有关的时间尺度因
子，幅值越大，时间尺度因子越大。为n阶厄米特多项式。其中，一阶与二阶hermite
‑
rodriguez函数分别对应双相和三相波形，其表达式如下所示：
[0035][0036][0037]
其中，a1,a2为归一化的运动单元动作电位的幅度因子。通过调节先验模板中幅度因子和时间尺度因子等参数得到的多种形态波形来近似模拟由不同大小运动单元所发放的运动单元动作电位序列，其中，设置模板中因变量t的取值范围为
‑
20～20ms，时间尺度因子α的变化范围为5～20，以更有效地逼近运动单元动作电位波形，图2则为一阶和二阶hr函数的模拟波形示意图，模板的简单的数学表达式，降低了复杂性，提高了算法的鲁棒性。
[0038]
步骤3.2，利用静息状态下的肌电信号作为参考信号对15％最大自主收缩力持续作用下的肌电信号进行有效峰值检测，并初始化幅度因子a，具体如下：首先选择参考信号的均值作为检测阈值，然后基于该阈值检测出15％最大自主收缩持续作用力下的肌电信号的有效峰值(波峰和波谷)，并按其幅值的绝对值大小进行降序排列。接着，按如下约束规则对模板的幅度因子a初始化处理：
[0039][0040]
arc max m
[0041]
s.t.|peak1|
‑
|peak
m
|≤2sd(rest_signal)
[0042]
其中，|peak
k
|表示在某个mvc力量水平下采集的肌电信号中第k个尖峰点的幅值绝对值，且满足|peak1|≥|peak2|≥
…
≥|peak
k
|≥
…
，a0表示幅度因子的初始值，sd(rest_signal)为参考信号的标准差。
[0043]
幅度因子a决定着muap波形幅值的大小，本发明中幅度因子以公比为1/2的等比数列形式进行变化，旨在描述不同muap的幅度水平，其约束条件如下所示：
[0044][0045]
arc max i
[0046]
s.t.a
i
≤mean(rest_peaks)
[0047]
式中，a
i
表示第i次迭代的幅度因子，mean(rest_peaks)表示参考信号中所有尖峰幅值的平均值。
[0048]
基于约束条件检测出15％最大自主收缩持续作用力下的肌电信号的有效峰值(波峰和波谷)，并按其幅值的绝对值大小进行降序排列。
[0049]
步骤3.3，在有效峰值附近截取信号段与先验模板进行匹配处理，模板匹配包括以下步骤：
[0050]
1.在确定幅值因子a
i
后，选取肌电信号中所有满足|peak
i
|
‑
a
i
≥0的尖峰脉冲点作为有效尖峰。传统做法都是以有效尖峰位置为对称点，通过一个固定长度的窗口截取信号段，再将该信号段与先验模板进行匹配分析。但是，不同幅度的运动单元发放运动单元动作
电位的持续时长与发放间隔不同，且同一幅度水平下不同形态的运动单元动作电位也并不一定关于尖峰点对称。因此，本实施例采用了带有可变窗长的滑动窗口在各个有效尖峰附近进行截取，获得部分待分析的信号段作为备选运动单元动作电位。
[0051]
2.采用动态时间规整算法计算每个备选运动单元动作电位与四个模板的相似性，并使相似性规范化。基于动态编程的思想，动态时间规整算法可以测量不同长度的两个时间序列的相似度。
[0052]
给定两个时间序列：运动单元动作电位候选x
l
(t)＝x1,x2,...x
i
,...,x
n
和经验模板y
k
(t)＝y1,y2,...y
i
,...,y
m
，它们的长度分别为n和m。动态时间规整算法通常使用动态编程算法，为了对齐这两个序列，构造一个n*m矩阵d，矩阵元素d
ij
表示xi和yj两点之间的振幅距离。该算法可以理解为寻找该矩阵中从左下角到右上角的路径中的最小路径。该路径是所通过的网格点的累计添加，从连续性和单调性来看，每个矩阵元素d
ij
只有三个前进方向：(i+1,j),(i,j+1)or(i+1,j+1)。动态时间规整算法的目标是使以下式子最小：
[0053][0054]
其中d表示备选运动单元动作电位与先验模板的距离,d
ij
是路径上的矩阵元素表示的值，k是路径经过的矩阵元素个数，此外，还分别计算了四个先验模板的平均振幅值mean
t
。
[0055]
3.结合动态时间规整算法和先验模板的平均幅值来确定备选运动单元动作电位是否与先验模板匹配。具体的匹配决策规则如下：
[0056]
如果满足d(x
l
(t),y
k
(t))<mean
t
/2，直接认为备选运动单元动作电位与先验模板匹配；否则直接认为备选运动单元动作电位与先验模板不匹配；
[0057]
4.决策完成后，根据与先验模板匹配的备选运动单元动作电位，选取备选运动单元动作电位的开始时间作为本次迭代获得的各个运动单元的发放时刻信息。
[0058]
步骤3.4，采用脉冲触发平均技术进行运动单元动作电位波形的重构，过程如下：
[0059]
将每个先验模板识别出的运动单元动作电位的发放时刻视为原始肌电信号的触发时间点，依次触发得到多个触发波形，然后对触发波形进行平均即可得到对应运动单元动作电位的重构波形。这些重构波形便可视为由各个运动单元发放的真实运动单元动作电位波形。
[0060]
至此，就得到了一组运动单元的发放信息和波形信息，15％最大自主收缩力作用下屈指浅肌采集的肌电信号的分解过程及结果如图3所示。由于此次迭代是基于相同幅度因子的，故这组运动单元的发放序列波形可被认为处于同一幅度水平。
[0061]
步骤3.5，通过比较幅度因子a
i
和参考信号尖峰点平均值mean(rest_peaks)的大小判断循环迭代过程是否继续。在每次迭代剥离运动单元动作电位的过程中，若a
i
高于mean(rest_peaks)，则视残差信号为更新后的肌电信号，并返回步骤3.2继续迭代；否则，循环迭代结束，进入之后评估步骤。
[0062]
本实施还包括步骤4，基于分解结果建立互相关系数评估指标验证所述肌电分解方法的有效性和可行性：
[0063]
将分解得到的多个运动单元动作电位序列分别按幅度水平从大到小和从小到大两种策略逐步进行线性组合，形成多个动作组合序列，依次计算运动单元动作电位组合序列与原始肌电信号间的互相关系数，通过观察互相关系数变化趋势来评价分解算法的优劣。若运动单元发放时刻与脉冲触发平均技术重构的动作波形符合真实运动单元发放规律，则当更多的动作序列组合在一起时，原始肌电信号与运动单元动作电位组合序列之间的互相关性也会随之增强。反之，在不断添加新的动作序列更新运动单元动作电位组合序列的过程中，一旦发现更新后的动作组合序列与原始肌电信号间的互相关系数降低，则表示新增的动作序列不能与真实运动单元发放信息相匹配。第i次组合的互相关系数r
i
的计算公式如下所示：
[0064][0065]
式中，x(t)表示原始肌电信号，x
mu(i)
(t)表示第i次组合的运动单元动作电位序列，与则分别表示x(t)和x
mu(i)
(t)的平均值，t为信号采样时长，k则表示分解迭代次数。
[0066]
按照分解得到的各个幅度水平的运动单元动作电位序列，分别以从大到小和从小到大两种策略构建运动单元动作电位组合序列，计算组合序列与原始肌电信号信号的互相关系数，根据系数变化趋势评估分解的有效性和可靠性，互相关系数结果图如图4所示。