基于先验模板的表面肌电信号渐进分解算法研究

罗万国 侯文生,# 郑小林,# 万小萍 陈海燕 周 平 吴小鹰*

1(重庆大学生物医学工程系， 重庆 400044)2(重庆市医疗电子工程技术研究中心， 重庆 400044)

基于先验模板的表面肌电信号渐进分解算法研究

罗万国1侯文生1,2#郑小林1,2#万小萍1陈海燕2周 平2吴小鹰1*

1(重庆大学生物医学工程系， 重庆 400044)2(重庆市医疗电子工程技术研究中心， 重庆 400044)

本研究通过将表面肌电信号(sEMG)分解为运动单元动作电位序列(MUAPTs)，来研究神经-肌肉系统中运动单元(MU)的募集与发放模式。针对高收缩力情况下MUAP叠加问题，首先采用FastICA算法和小波包去噪算法对信号进行预处理；然后基于先验知识构建了4种形态的可伸缩MUAP模板；最后，采用“先大后小”的渐进识别方式，逐个对MUAP进行自动提取。在此基础上，还将该算法应用于8名受试者(3组/人)不同手指活动模式下的指浅屈肌多通道(12通道)sEMG分解；单通道分解结果显示，高力量水平下sEMG中的主体MUAPTs能够被有效检测和分类；统计结果证实，随着力量水平的增加，MUAP的数目增加；不同大小MUAP比重的变化与活动手指和力量水平具有显著的相关性。本文的实验结果，初步验证了利用先验模板从sEMG中渐进提取MUAP的可行性，为sEMG分解和进一步研究MU发放规律提供于一种新的思路。

表面肌电信号；分解；先验模板；运动单元动作电位

引言

表面肌电信号(surface electromyography, sEMG)，是由肌肉兴奋时募集的多个运动单元(motor unit, MU)产生的动作电位(action potential, AP)，在表面电极处的综合叠加[1]。由于sEMG具有无创的特点，在假肢控制[2]、神经肌肉系统功能异常的辅助诊断[3]、运动康复诊断[4]等方面受到重视。虽然直接提取sEMG的时域、频域和时频域特征，可用于描述肌肉活动及其功能状态[4-5]，但为了进一步认识MU募集的细节信息，通过分解sEMG来还原构成它的主体运动单元动作电位序列(motor unit action potential trains, MUAPTs)，在近年来成为了一个重要的研究方向[1, 3, 6-13]。

基于统计学的盲源分离(BSS)算法，尤其是独立分量分析(ICA)，是比较典型的sEMG分解算法之一。日本坂田大学的Garcia研究小组，采用一系列ICA算法对仿真信号和实测sEMG进行了分解[6-7]，得到了一些与主体MUAPTs相似的分离信号；中国科技大学的神经肌肉控制实验室，将多种BSS算法应用于多通道sEMG分解，得到了4类MUAP信息，并较好地解决了波形叠加问题[11]；重庆大学安媛等利用FastICA算法，从6通道的指浅屈肌sEMG中提取了4种MUAP信息，并初步研究了力量水平对MUAP募集模式的影响[12]。然而，由于ICA算法要求观测信号中有，且仅有一个高斯源信号，观测矢量的个数还不小于源矢量的个数，故基于ICA的多通道sEMG分解，大多在较低力量水平进行。此外，由于需要反复迭代，其分解效率也有较大局限[11, 14]。

另一种典型的sEMG分解思路，是基于波形检测和模式识别的方法对sEMG进行分解，其中最具代表性的是De Luca等在插入式肌电信号(iEMG)分解方法的基础上，提出来新的sEMG分解策略[1]。最新报道显示，该算法可用于不同肌肉、不同力量水平(最高100% MVC)的sEMG分解，所分解的MU类型数高达40，准确度也高达97%[8]。但是，也有学者认为其准确性仍需改进[9]。本文结合一些现有sEMG处理与分解算法的优点，以MUAP波形结构为先验知识，设计了一种渐进式sEMG分解算法：通过模板拟合和逐次迭代提取MUAP信息，以提高MUAP识别效率；并通过不同手指活动模式下采集的指浅屈肌sEMG，对算法的有效性进行了初步验证。

1 渐进分解算法设计

1.1模板建立

MUAP波形被证实大多呈现双相或三相结构[1, 8, 10, 12]，MU遵循“大小原则”进行募集，后募集的MU具有更大的激活阈值且产生尺寸更大的MUAP。为了使所建立的模板波形更具有代表性，提取出更多MU所发放的MUAP信息，依据现有报道的结果[1, 8, 10, 12]，构建了如图1所示的4种基本形态的模板波形(Tp1～Tp4)，且每种形态的模板波形，通过幅度和时间伸缩以代表不同激活阈值MU发放的MUAP。

为了计算方便，将Hermite-Rodriguez函数进行归一化和变形来表示每种形态MUAP的波形，上述4种模板的近似数学表达式为

(1)

式中，Tpi(t)表示第i个MUAP模板为了更好地逼近模板波形，变量t的取值范围取：[-2.5,2.5]，A为MUAP的幅度伸缩因子，采用式(2)所示一个比值为1/2的等比数列来代表不同大小MUAP的幅值A0可通过MVC力量水平下采集到的sEMG由式(3)进行估计，N的大小则取决于式(6)中Thrmin的大小；λ1,λ2,λ3,λ4为与MUAP幅值有关的时间伸缩因子，幅值越大，时长相对越长。参考相关文献[1, 8]，其总变化范围设为：5～20 ms。

(2)

(3)

式中，peakm为MVC力量水平下采集的sEMG峰值，M为满足|peak1|-|peakm|≤2σ的最大m值，σ为基线的方差。

1.2分解步骤

sEMG经常表现为不同幅值MUAP的相互叠加，幅值大的MUAP波形相对容易被检测出来。因此，采用如图2所示的两层循环迭代算法，按照“从大到小”的顺序，自动对sEMG中的MUAP进行逐个剥离。

步骤1：阈值初始化。比较输入信号峰值与MVC力量水平下sEMG估计出的MUAP幅值序列，确定出输入信号中所包含的最大MUAP的幅值，并以此作为备选峰值检测点的初始阈值。

步骤2：备选峰值点检测。所谓的备选峰值点检测，就是在sEMG中检测绝对峰值大于阈值的峰值点的位置，然后按照峰值由大到小的顺序，将其对应的时刻点保存到1个一维数组中。

步骤3：MUAP识别与提取。从峰值索引数组中得到一个峰值索引之后，通常的做法是采用一个固定窗长的窗口，在索引点对称截取出一个待分析信号段与模板波形进行匹配。然而，不同层次的MUAP的持续时间不同，同一层次不同形态MUAP关于峰值点并不都是对称的，而且不同波形其峰值到左右两边比例也不同(如图1所示)。因此，采取移动的可变窗长截取索引点附近的待分析信号段。除此之外，为了方便进行参数计算，将不同窗长截取的MUAP波形都进行了3次样条插值，使其与模板波形具有相同的点数。

待分析信号段与模板波形的匹配程度，主要采用相关系数与残差比进行判定，两个参数分别由式(4)和式(5)进行计算。

(4)

(5)

式中，Tpi(t)表示第i个MUAP模板，xj(t)表示第j个信号片段，cov(·)表示协方差运算，D(·)表示方差运算，C表示相关系数矩阵，R表示残差比矩阵。

得到相关系数矩阵C和残差矩阵R之后，采用如下规则进行MUAP波形判定：

(1)如果最大相关系数大于0.95，则判定具有最大相关系数的模板与信号片段匹配；

(2)如果最大相关系数介于0.7与0.95之间，则判定相关系数大于0.7且残差比最小的模板与信号片段匹配；

(3)如果最大相关系数小于0.7，则判定该峰值索引点附近无MUAP发放。

识别和提取完某一峰值索引点的MUAP之后，则将残余sEMG用于下一峰值索引点的MUAP识别与提取。当遍历所有备选峰值点之后，退出子程序，完成该阈值下的MUAP识别与提取。

从可分离的角度来看，当所发放MUAP的幅值小于基线时，就很难检测和识别到相应的MUAP。因此，利用式(6)计算的最小阈值来判定是否结束分解。

(6)

式中，peakn为基线信号中的峰值点，N为满足|peak1)|-|peakn|≤σ的最大n值。

2 sEMG检测及分解

2.1sEMG采集

为了验证算法的有效性，实验募集了8名(男性4名，女性4名，年龄20～24岁)身体健康的大学生志愿者进行sEMG采集实验，受试者均知情同意。实验中，受试者按如图3(a)所示姿势分别使用食指、中指和环指进行20% MVC，40% MVC和60% MVC指力跟踪实验(详见文献[12])。采用图3(b)所示的粘贴在指浅屈肌(FDS)肌腹25%～76%处的一块6×2(行×列)Ag-AgCl柔性阵列电极进行sEMG检测。检测到的sEMG经过信号放大电路之后，采用RM6280C多道生理记录仪进行记录、显示和存储。同步记录的还有以相同采样率(2 kHz)采集的4个指力传感器的输出电压信号。

2.2预处理

作为一种弱生物电信号，在采集sEMG的过程中不可避免地受环境噪声、设备固有噪声、工频干扰、运动伪迹等噪声的影响。因此，对原始sEMG依次进行了有效段自动提取、带通滤波、FastICA去工频干扰和小波包去噪等预处理来提高信号的质量，具体处理过程与效果参见文献[15]。

2.3MUAP分解

根据去噪后的基线信号和MVC力量水平下采集的sEMG，估计得到可分离MUAP的最大幅值约等于4 mV，最小幅值约等于1 mV。为了计算方便，可分离MUAP的幅值分别设为4、2、1 mV。对于12通道的sEMG，采取的是单通道分别处理的方式。单通道sEMG分解可以概括为4个步骤。

步骤1：输入待分解sEMG，根据信号最大峰值确定其中所包含的最大MUAP幅值，并以该幅值作为初始阈值进行备选峰值检测；

步骤2：检测信号中绝对峰值大于阈值的峰值点的位置，然后按照峰值由大到小的顺序，将所有备选峰值点对应的索引点保存到1个一维数组中；

步骤3：从最大峰值索引点开始，在索引点附近根据4种不同幅值的MUAP模板，选择不同变化范围的窗长(4 mV：15～20 ms；2 mV：10～16 ms；1 mV：5～12 ms)移动截取待分析信号片段，计算特征参数，识别MUAP；若在某个备选峰值点附近识别出有效的MUAP，则从信号中减去相应的模板；若未识别出有效MUAP，则进入下一个索引点进行识别，直到遍历所有备选峰值点；

步骤4：检测残余信号峰值，若残余信号峰值大于最小检测阈值，则以更小MUAP的幅值作为备选峰值点的检测阈值，返回步骤2继续对残余信号进行分解；否则，结束分解。

3 结果

3.1单通道分解结果

图4为一段60%力量水平下，食指恒力收缩时采集到的sEMG，经过预处理和分解的效果图。从图4中可以看到，预处理之后的信号无明显的工频干扰，高频噪声明显降低，曲线更加平滑；另外，渐进逐层分解后得到了幅值约为1、2、4 mV等3层MUAPTs，共12个MUAPTs加上一个幅值和能量较少的残余信号，且残余信号中无明显MUAP波形残存；从不同MUAP的发放时间序列上可以看到，所提取出的MUAP与原始波形在时间轴上具有很好的对应关系；仔细观察右侧的细节放大图，可以清楚地看到，叠加波形中的主体MUAP都被成功地分解出来，并归类到了相应的MUAPTs中，且其发放时刻定位准确。由此可见，所提出的“先大后小”的渐进式分解算法，可以有效地实现叠加波形的分解。

3.2MUAP发放总数随力量水平的变化趋势

分别统计了8名受试者(3组/人)食指、中指和环指在20% MVC、40% MVC和60% MVC等3个力量水平下采集到的时长为2 s的12通道sEMG分解得到的MUAP数目，结果如图5所示。显然，对于所有受试者，3种手指活动模式下，随着力量水平的增加，MUAP发放的数目都呈增加趋势。

3.3不同手指活动模式下不同大小MUAP发放比重变化

从图5可以观察到，MUAP发放总数总是随着力量水平的增加而增加。然而，不同活动模式下，所分解得到的3种不同大小MUAP的变化并不都随着力量水平的增加而增加。由于不同的受试者，不同活动模式下分解得到的MUAP的总数不同，因此以3种幅值不同的MUAP数目占总数目的百分比(即比重)为指标，统计了8名受试者在不同活动模式下(每种活动模式3组重复试验数据，共24组数据)不同MUAP的贡献。结果如图6所示。

从图6所示的统计对比分析结果可以观察到，3种手指活动模式下，大幅值MUAP(A≈4 mV)的发放比重，随着力量水平的增加而增加；对于中等幅值(A≈2 mV)MUAP，随着力量水平的增加，不同手指则表现出不同的变化模式：食指，显著性增加；中指，20% ～ 40% MVC显著性增加，40% ～ 60% MVC则基本保持稳定；环指，20%～40% MVC显著性增加，40% ～ 60% MVC显著性减少。小幅值(A≈1 mV)MUAP的发放比重，随着力量水平的增加而减小。

4 讨论

近年来，随着盲源分离、人工神经网络和人工智能等新兴技术的不断引入，sEMG分解研究取得了很大的进步，但MUAP叠加波形分解仍是sEMG分解中难点问题之一。针对这一问题，本研究设计了一种渐进式sEMG分解算法，基于先验模板匹配和分级量化，实现了对不同力量水平下的MUAP提取。实测sEMG分解结果显示，与其他报道给出的结果[1, 8, 16]一样，随着力量水平的增加，分解得到的MUAP总数呈增加趋势。由此可见，本研究的分解算法提取出的主体MUAPTs，确实能够从整体上反映肌肉活动的强弱变化。

目前，基于MUAP形态的sEMG分解算法，在处理叠加波形时主要有两种思路：一种是将不同的MUAP模板进行不同程度的叠加合成再与sEMG进行比对；另一种则是从sEMG中逐个剥离MUAP[16]。本研究的基本思路与后者相同，但对分解策略进行了优化，以适应高力量水平的sEMG分解。首先，本研究采用了一种简洁的类似二进制的分级量化方法设置不同MUAP的幅值；其次，传统的sEMG分解算法在波形分割时，通常采用对称的固定窗口进行信号分割，忽略了不同MUAP时长不同，对称性不同[17]为了弥补其不足，在同一个峰值点，采用不同窗长，移动截取的方式，截取多个待分析波形，以便更加准确地对主体MUAP进行识别；最后，在波形识别时，相关系数只作为基本判定指标，还考虑了残差的大小，能够很好的减少误判。从图4中60% MVC力量水平下的sEMG分解结果来看，所提出的算法与其他基于波形检测和模式识别的方法[1, 8]一样，所得到的MUAPTs与预处理后的sEMG中的MUAP叠加波形，在时间轴上具有很好的对应关系，同一时刻同时发放的MUAP也得到了有效分解。另外，分解之后得到残余分量较小，无明显MUAP残余。综上，本研究所提算法分解得到的主体MUAPTs中MUAP的发放时刻定位准确，有效地实现了叠加波形的分解，适用于高力量水平下的sEMG进行分解。

由图6可知，同一手指活动模式下，随着力量水平的增加，不同大小MUAP的比重不同，即不同大小MUAP对手指力量的贡献率不同。小幅值(A≈1 mV)、持续时间短的MUAP，在较低力量水平下(20% MVC)的贡献率最大，随着力量水平的增加，贡献率降低；大幅值(A≈4 mV)、持续时间长的MUAP，在较低力量水平下(20% MVC)的贡献率最小，随着力量水平的增加，贡献率增加。产生上述结果的原因可能是，在较低力量水平下，激活阈值较小的MU被激活并发放幅值较小的MUAP，随着力量水平增加，为了产生更大的肌力，具有更高激活阈值的MU被激活，并且主要由后激活的MU来产生肌力。事实上，已有结果证实，MUAP的大小主要是由MU中包含的肌纤维数决定，为了产生更大的肌力，包含有更多肌纤维数的高阈值MU将被激活[16]。

同时，本实验结果还表明，不同手指活动模式下，部分MUAP贡献率随着力量水平的变化具有显著性差异。如图6所示，3种手指活动模式下，中等幅值MUAP(A≈2 mV)贡献率随着力量的变化趋势具有显著性差别。产生这一结果的原因，可能是不同手指活动时，神经-肌肉系统采取了不同的MU募集策略，选择性地激活MU[18]。

5 结论

针对高收缩力情况下sEMG分解中的叠加波形分解问题，本研究基于MU募集的“大小原则”和MUAP形态，提出一种“从大到小”的渐进式MUAP分解算法。为了验证算法的有效性，对不同手指活动模式下的指浅屈肌sEMG进行了分解。结果显示：该算法能有效分解sEMG；随着力量水平的增加，分解得到的MUAP数目呈现增加趋势；不同大小MUAP的比重随着活动手指和力量水平的变化而变化。所提出的分解算法，为深入研究神经-肌肉控制系统中的MU募集与发放信息提供了可能。而所提出算法MUAP基本模板类型较少，只提取了一部分MU所发放的MUAPTs信息，后续研究中需要建立有效的模板更新机制，实现更多形态MUAP的提取。此外，后续研究需增加样本量，并设计实验对自主收缩时MU的募集模式进行研究。

[1] De Luca CJ, Adam A, Wotiz R,etal. Decomposition of surface EMG signals[J]. Journal of Neurophysiology,2006,96(3):1646-1657.

[2] 侯文生，万莎，吴小鹰，等. 手指运动姿态检测及假肢手指动作实时控制[J]. 重庆大学学报,2013(04):103-109.

[3] 姚博. 表面肌电信号分解算法及其在小儿脑瘫评诂中的应用[D]. 合肥：中国科学技术大学,2012.

[4] 王耀兵，季林红，黄靖远. 康复训练过程中上肢运动肌电信号特征的对比分析[J]. 中国临床康复,2004(23):4721-4723.

[5] 王刚，王志中，胡晓，等. 基于最佳小波包的表面肌电信号分类方法[J]. 中国医学物理学杂志,2006(01):45-48.

[6] Garcia GA, Okuno R, Akazawa K. A decomposition algorithm for surface electrode-array electromyogram-A noninvasive, three-step approach to analyze surface EMG signals[J]. IEEE Engineering in Medicine and Biology Magazine,2005,24(4):63-72.

[7] Naik GR, Arjunan S, Kumar D. Applications of ICA and fractal dimension in sEMG signal processing for subtle movement analysis: a review [J]. Australasian Physical & Engineering Sciences in Medicine, 2011,34(2):179-193.

[8] Nawab SH, Chang SS, De Luca CJ. High-yield decomposition of surface EMG signals[J]. Clinical Neurophysiology, 2010,121(10):1602-1615.

[9] Farina D, Enoka RM. Surface EMG decomposition requires an appropriate validation [J]. Journal of Neurophysiology, 2011,105(2):981-982.

[10] 安媛. 手指活动影响前臂多腱肌运动单元募集模式的初步研究[D]. 重庆：重庆大学,2012.

[11] 李强. 表面肌电信号的运动单位动作电位检测[D]. 合肥：中国科学技术大学,2008.

[12] 侯文生，安媛，杨丹丹，等. 食指力量对指浅屈肌运动单元募集模式的影响[J]. 医用生物力学,2012(06):648-654.

[13] Merletti R, Holobar A, Farina D. Analysis of motor units with high-density surface electromyography[J]. Journal of Electro- Myography and Kinesiology,2008,18(6):879-890.

[14] 姚博，杨基海，陈香，等. 基于稀疏分量分析的欠定盲源分离用于表面肌电信号分解[J]. 航天医学与医学工程,2012(02).

[15] 罗万国，吴小鹰，杨丹丹，等. 基于形态学的表面肌电信号分解算法研究[Z]. 成都:电子科技大学生命科学与技术学院,2013.

[16] Stashuk D. EMG signal decomposition: how can it be accomplished and used? [J]. Journal of Electromyography and Kinesiology, 2001,11(3SI):151-173.

[17] Chauvet E, Fokapu O, Hogrel JY,etal. Automatic identification of motor unit action potential trains from electromyographic signals using fuzzy techniques[J]. Medical & Biological Engineering & Computing, 2003,41(6):646-653.

[18] Butler TJ, Kilbreath SL, Gorman RB,etal. Selective recruitment of single motor units in human flexor digitorum superficialis muscle during flexion of individual fingers[J]. Journal of Physiology-London, 2005,567(1):301-309.

StudyonScaledDecompositionofSurfaceEMGBasedonPriorTemplates

LUO Wan-Guo1HOU Wen-Sheng1,2#ZHENG Xiao-Lin1,2#WAN Xiao-Ping1CHEN Hai-Yan2ZHOU Ping2WU Xiao-Ying1*

1(DepartmentofBiomedicalEngineering,ChongqingUniversity,Chongqing400044,China)2(ChongqingEngineering&TechnologyResearchCenterforMedicalElectronics,Chongqing400044,China)

surface electromyography; decomposition; prior templates; motor unit action potentials

10.3969/j.issn.0258-8021. 2014. 03.015

2014-02-05， 录用日期:2014-04-20

国家科技支撑计划(2012BAI16B02)；国家自然科学基金(30970758)；重庆市自然科学基金(cstc2012jjA10103)

R318

D

0258-8021(2014) 03-0366-07

#中国生物医学工程学会会员(Member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)，E-mail: wuxiaoying69@163.com