基于肌电信号的肘关节肌肉疲劳算法性能研究

孟庆丰，陶庆，2，来全宝，胡玉舸

（1.新疆大学机械工程学院，新疆 乌鲁木齐 830047；2.新疆大学机械工程博士后研究中心，新疆 乌鲁木齐 830047）

1 引言

肌肉疲劳是人们日常生活中一种常见的生理现象，一般是指由长时间运动或发力引起的肌无力或酸痛的感觉[1]。日常生活中我们或多或少都会遇到肌肉疲劳的现象，但是如果不顾肌肉的疲劳而过度运动或长时间活动的话，很可能会造成肌体的永久性损伤，更甚者会带来生命危险[2]，因此对于肌肉疲劳评价的研究足以引起我们的高度重视。对于肌肉疲劳评价方法的研究与分析，无论是在体育训练还是在康复医学、运动医学等领域都具有重要意义。

如今，已知通过记录表面肌电信号计算疲劳指数的方法可以分为两个主要类别[3]。第一类由经典的基于信号可量化频谱变化的算法组成。例如经典的平均频率以及频谱矩比[4]和基于小波的方法[5]都属于此类。第二类考虑信号的非线性成分。例如考虑信号重复模式或分形测量的递归量化分析和基于熵的方法的模糊近似熵都属于此类[6]。

由于上肢的运动功能远比下肢复杂，上肢在人体日常的工作生活扮演着重要的角色，而肘关节的屈伸动作甚至出现在人体运动的每个瞬间，并且上肢中肘关节的运动减少50%将导致上肢运动功能降低近80%[7]，肘关节的正常运动对人体的工作生活是非常重要的，因此探究更有效的肘关节肌肉疲劳评价方法迫在眉睫。

为此，国内外学者已经针对肘关节肌肉疲劳评价方法展开了大量的研究，不同方法得出的评价指标获得的肌肉疲劳评价效果不同。文献[8]利用AEMG与MPF等线性和Lempel−Ziv 复杂度与DET等非线性指标分析不同负载下肱二头肌疲劳过程中表面肌电信号的变化，结果显示疲劳指标的变化具有非疲劳特异性。

文献[9]采集不同状态的上肢两侧肌电信号，将时域的积分肌电值、均方根振幅和频域的平均功率频率、中位频率作为疲劳评价指标，研究得出上肢一侧近端肌肉疲劳会使另一侧的运动信号产生伴随性变化，即产生中枢共驱动现象。

文献[10]分析时域的积分肌电值、均方根振幅和频域的平均功率频率、中位频率等指标随着疲劳程度加深的变化趋势，提出一种肌肉疲劳出现时刻的估计方法。文献[11]将时域的iEMG、RMS与频域的MPF、Fc 分析评价指标，利用SPSS 统计分析时域与频域之间的参量关系以及上肢肌肉疲劳间的关系，结果与实际情况更贴合。

以上研究主要是对上肢疲劳时域或频域指标变化规律的探索，缺乏对各疲劳指标的评价效果进行分析。为此，研究了肘关节肌肉疲劳评价方法的评价效果，记录测试对象在不同负载下肘关节做屈曲运动的sEMG 信号，计算MNF、SMR、WIRMI1551、fApEn与RQA%DET的指标，并进一步比较5种方法对肘关节肌肉疲劳的评价效果。

从而获得更有效的肘关节肌肉疲劳评价方法。将从肌肉疲劳评价方法的抗干扰性与区分疲劳程度的能力两个方面对5种疲劳评价方法的效果进行比较，分析所得结果，为肘关节肌肉疲劳提供一个更加准确的人因工程定量评价方法。

2 研究方法

2.1 肌肉疲劳评价方法

2.1.1 平均频率

研究表明，随着肌肉疲劳的累积，sEMG 信号功率谱的主频率由高转向低，平均功率值逐渐变小，并且疲劳程度越大，下降趋势越强。基于对功率谱密度（PSD）的估计，计算了谱疲劳评价方法指标。PSD 是利用Welch 方法估算的（其中k=15），并使用Hamming窗函数进行50%的重叠[12]。基于PSD估计的MNF：

式中：f—频率；fs—采样率。

2.1.2 频谱矩

作为平均频率的替代，频谱矩（SMR）方法被Dimitrov等人描述为新的光谱指数，它是对H/L比率方法的改进。文献[13]用平滑函数加权PSD估计，他们建议使用频率的倒数作为低频带的加权函数，高频带的加权函数是频率平方，立方，四次或五次方。基于阶数k的谱矩Mk：

频谱矩为：

式中：f1=5Hz和f2=500Hz—积分的上下限，阶数k=5。

2.1.3 基于离散小波的方法

量化频率压缩的其他增强疲劳评价方法是Gonzalez−Izal等人报道的基于小波的方法[14]。它们是基于固定小波变换（SWT）的频谱估计而不是PSD 估计。WIRM1551 比值是从Gonzalez−Izal等人建议的小波方法的变体中选择的。它被定义为：

式中：D1—使用sym5小波的SWT的第一尺度的功率谱；D5—SWT的第五尺度的功率谱。

正如建议的那样，f1=8Hz和f2=500Hz为积分的上下限。

2.1.4 模糊近似熵

熵估计是用于疲劳评价的另一种非线性方法。样本熵（Sam⁃pEn）是一种计算时间序列熵的方法，为条件概率的负自然对数，该条件概率是从一个纪元开始的两个序列在m维向量相似，而在m+1 维向量仍然相似，当向量的长度m=2，阈值r=0.3 为最佳搭配。为了疲劳检测的目的，文献[15]报告了一种改进的算法，称为模糊近似熵（fApEn）。与样本熵相反，如果两个序列相差任意恒定偏移量d，则模糊近似熵会将它们视为相似。以m=2，r=0.6为最佳组合计算模糊函数来获得模糊近似熵的近似度。模糊函数为：

2.1.5 递归量化分析

递归量化分析（RQA）是一种试图量化分析信号中重复模式的非线性疲劳评价方法[6]。首先构造距离矩阵，从时间信号的角度导出长度为D（嵌入维）的向量，其中时间延迟间隔λ=5，嵌入尺寸D=15。然后矩阵通过其平均值进行归一化处理，研究表明，阈值r=0.75时矩阵会产生仅包含零和一的二进制递归映射。%DET是用于评估复发图的，其值为形成一定链长的对角线与矩阵中所有对角线的数量之比。研究发现，在疲劳任务下，RQA%DET值的上升与运动单元的更长时标同步非常吻合。

2.2 疲劳评价方法的定量化分析

由于疲劳评价方法难以基于测试对象的主观感受来校准，因此本研究使用定量化分析方法比较疲劳评价方法的差异性。定量化分析方法主要包括同一疲劳条件下评价方法的抗干扰性及不同疲劳条件下区分疲劳程度的能力。相关研究表明，肌肉疲劳评估指数呈现具有恒定斜率的线[16]，本研究利用最小二乘法对五种疲劳评估方法进行了线性回归分析，并求解出每个方法的确定系数R2与回归方程的斜率k。

2.2.1 抗干扰性分析

疲劳评价方法预期在疲劳情况下产生平滑变化的疲劳指数。偏离这个平滑过程被认为是干扰。疲劳指数的干扰可以通过确定系数R2来量化[17]。R2表示具有恒定斜率的线反映信号数据的相关性，R2的数值范围从0到1。R2的值为0意味着信号中没有线性趋势。R2的值为1则意味着信号恰好是具有非零斜率且没有额外干扰的直线。因此R2的值越大表明疲劳评估算法的抗干扰性越强。

2.2.2 区分疲劳程度能力分析

人体肘关节肌肉在不同负载下会产生不同的肌肉疲劳程度。柯尔莫哥洛夫−斯米洛夫检验（Kolmogorov−Smirnov test，即K-S检验）是一种基于概率统计理论的方法，它能够计算两个数据经验分布函数间的最大垂直距离Lmax，Lmax越大说明两个数据相似性越小。

本研究利用K-S检验分析五种评价方法区分疲劳程度的能力[18]。即对五种评价方法在不同负载下的线性回归方程斜率k进行K-S检验得到的Lmax值，Lmax越大表明疲劳评估方法对不同负载下区分疲劳程度的能力越强。

3 肌电信号采集与处理

3.1 测试对象与试验准备

3.1.1 测试对象

为了采集数据分析比较不同疲劳评价方法，本研究采集了12名测试对象（6名女性和6名男性）肱二头肌的表面肌电信号。测试对象的年龄在（19～26）岁之间，无皮肤过敏史且肱二头肌皮肤表面无创伤和疤痕，上肢运动能力处于正常状态。在试验开始之前，所有测试对象都书面同意志愿参加本次试验。

3.1.2 试验准备

为了施加恒定的力，要求每个测试对象以同一固定姿势托举哑铃。实验过程中要求测试对象保持直立坐姿，肘部处于放松状态后抓举哑铃进行实验。有关动作设置，如图1所示。

图1 实验设定Fig.1 Setup of the Experiment

数据采集过程中，首先引导测试对象以最大可实现的用力在短时间内抓举哑铃以单独测取每个测试对象的最大自主收缩力（MVC）。然后在随后的三次试运行期间，根据每个测试对象的MVC调整他们要托举的重量。在三次试验中，要求测试对象分别托举他们MVC 的20%，40%和60%的重量。在20%和40%MVC 负荷的情况下，所有测试对象持续3min 的指定时间，在60%MVC负荷水平下持续135s达到疲劳。每次试验后，测试对象休息30min。

3.2 试验设备与数据采集

3.2.1 试验设备

葡萄牙Plux公司生产ErgoPlux的8通道无线表面肌电测量系统、上海励图医疗器材有限公司生产的LT−7型理疗电极片。

3.2.2 数据采集

测试肌肉位置，如图2所示。先用75%的医用酒精擦拭测试对象右臂肱二头肌的皮肤表面，清理皮肤表面污渍与角质，再使用干燥的棉球擦净皮肤表面残留的酒精后将八组电极片以高度密集的形式放置在肱二头肌上方。即电极以彼此相邻的两列连接，且彼此之间没有任何间隔。根据其直径，相邻放置的电极的中心间隔为26mm。测试电极沿着所测试肌肉的肌纤维方向，贴在肱二头肌最隆处，接地电极贴在肘关节处。试验过程中，每名测试对象依次托举他们最大负载（MVC）的20%、40%与60%重量，利用8 通道的表面肌电测量系统同步记录测试对象的sEMG信号。

图2 高密度电极片的贴点位置图Fig.2 Location of the High−Density Electrode Sheet in the Biceps

3.3 数据处理

经多次试验，测试对象在测试结束后均在不同负载下感受到不同程度的肌肉酸痛僵硬，即肌肉处于疲劳状态。为保证所分析的sEMG信号数据均在疲劳状态下，研究采用测试过程中最后的1min进行数据分析。首先同步采集8个通道的sEMG 会产生基线偏移现象，利用sgolay滤波器消除趋势化，然后sEMG信号采集过程中会受到50Hz工频及其他信号干扰，利用50Hz陷波与巴特沃斯带通滤波器对sEMG信号去噪[18]，最后将处理后的8通道信号取平均值。处理后的示例信号，如图3所示。

图3 处理后的示例信号Fig.3 Exemplary Signals After Processing

根据五种疲劳评价方法的定义，使用Matlab R2014b分别编写五种疲劳评价方法的算法，计算出每个测试对象在不同负载下各疲劳评价方法的指标值，然后对各疲劳评价方法进行定量化分析。定量化分析分为两个部分，一是利用最小二乘法对各方法的指标进行线性回归分析（图4 展示了1 名测试对象在不同负载下最后一分钟的MNF最小二乘法线性回归分析结果，其他4 种方法与MNF 做相同的处理），比较线性回归方程的确定系数R2的均值，利用SPSS25对R2进行单因素方差分析，进一步评价疲劳评价方法的抗干扰性；二是对回归方程的斜率k进行K-S检验得到不同负载下的Lmax值，并对Lmax进行单因素方差分析，评价5种疲劳评价方法区分疲劳程度的能力。图4中=0.792，正如预期的那样，随着负载水平的增加，R2值会增加，因为负载水平越高，MNF的线性趋势越强。

图4 MNF的线性回归分析结果Fig.4 Linear Regression Analysis Results of MNF

4 结果

4.1 抗干扰评价结果

为量化不同疲劳算法的抗干扰性，本研究计算了12个测试对象在不同负载下各疲劳评价方法线性回归方程的确定系数R2平均值与标准差，并对R2进行单因素方差分析结果，如图5所示。

图5（a）显示在20%MVC情况下，SMR的R2平均值为5种疲劳评价方法的最大值，最大平均值=0.237，且与MNF（p=8.28e−04）和RQA%DET（p=5.45e−05）方法存在显著性差异。表明SMR在20%MVC下具有最好的抗干扰性。WIRM1551具有较大的R2平均值=0.223，说明WIRM1551在20%MVC下具有较好的抗干扰能力。RQA%DET的R2平均值最小，表明RQA%DET的抗干扰能力最差。

图5（b）显示在40% MVC 情况下，SMR的R2平均值仍然最大=0.328，且与MNF（p=1.38e−05）和RQA%DET（p=3.69e−06）在评价肘关节肌肉疲劳的抗干扰性的差异具有统计学意义。表明SMR 在40%MVC 下具有最强的抗干扰能力。MNF、WIRM1551与fApEn的R2平均值相差很小，分别为=0.289，=0.288，说明这3种方法在40%MVC下具有相同的抗干扰能力。而RQA%DET的抗干扰性最差。

图5（c）显示MNF、WIRM1551与fApEn在60%MVC具有较大的R2平均值，但SMR仍然为最大值=0.528，且与RQA%DET（p=3.08e−07）的差异具有统计学意义。说明MNF、WIRM1551与fApEn在60%MVC下具有较好的抗干扰性，但SMR的抗干扰能力是最强的。RQA%DET依旧是5种方法中抗干扰性最差的。

图5 五种肌肉疲劳评价方法在不同负载下的R2Fig.5 Five Muscle Fatigue Detection Methods for R2 Under Different Loads

综上所述，SMR在不同负载下的R2平均值均为最大值，拥有最好的抗干扰性。WIRM1551 在不同负载下的R2平均值仅比SMR 小一些，同样具有较好的抗干扰性。MNF 与fApEn 仅在大负载下拥有较大的R2平均值，但是在小负载，低疲劳情况下的R2平均值相比SMR还是有较大差距的，因此，MNF与fApEn仅在大负载下拥有较强的抗干扰能力。而RQA%DET在不同负载下的抗干扰能力是最弱的。

4.2 疲劳程度区分能力评价结果

为量化分析5种疲劳评价方法区分疲劳程度的能力，对线性回归方程的斜率k进行K-S检验，得出不同负载下肘关节肌肉疲劳指标经验概率分布函数的最大垂直距离Lmax，并对Lmax进行单因素方差分析。全局结果，如图6所示。局部Lmax量化分析结果，如图7所示。

图6 全局Lmax量化分析Fig.6 Gl obalLmaxQuantitativeAnalysis

不同负载下全局Lmax量化分析结果，如图6所示。结果显示在疲劳程度较大（40%vs60%）的情况下，MNF、WIRM1551与fAp⁃En的Lmax平均值较大，SMR的Lmax平均值为最大值；在疲劳差异大（20%vs60%）的情况下，SMR与WIRM1551的Lmax平均值为最大值；在疲劳程度小（20%vs40%）的情况下，SMR的Lmax平均值呈现明显的尖点。而RQA%DET的Lmax平均值在三种情况均为最小值。全局Lmax量化分析结果表明在5种方法中SMR拥有最强的区分肘关节肌肉疲劳程度的能力，而RQA%DET是最弱的。

20%与60%MVC的Lmax量化分析结果，如图7（a）所示。结果显示SMR 与WIRM1551的Lmax为最大值，均为0.883，表明SMR与WIRM1551 在区分20%与60%MVC 的疲劳程度的能力最出色。而RQA%DET 是5 种方法中评价疲劳程度效果最差的。SMR与RQA%DET的差异具有统计学意义（p=1.09e−02）。表明SMR与WIRM1551在疲劳程度差异较大时有最强的区分疲劳程度能力。40%与60%MVC的Lmax量化分析结果，如图7（b）所示。结果表明MNF、SMR、WIRM1551 与fApEn 均呈现较大的Lmax平均值，说明4 种方法在评价大负载引起的疲劳程度均有很出色的能力。

其中SMR 是最出色的，Lmax=0.766，且与RQA%DET（p=8.93e−03）存在显著性差异。而RQA%DET区分疲劳程度的能力最差。说明在疲劳程度较大时，MNF、SMR、WIRM1551与fApEn均具有较好的区分能力，且SMR 是最好的。20%vs40%MVC的Lmax量化分析结果，如图7（c）所示。

结果表明SMR的Lmax平均值最大，最大值Lmax=0.622，且与RQA%DET（p=1.6e−03）存在显著性差异。RQA%DET的Lmax平均值仍是最小的。说明在疲劳程度较小的情况下，SMR仍然拥有最强的区分肌肉疲劳程度的能力。

综上所述，SMR相较于另外4种方法在不同负载下都具有最大的Lmax平均值，说明SMR在区分不同负载引起的肘关节肌肉疲劳程度的能力最出色。WIRM1551 虽没有SMR 那么出色，但在不同负载下均具有较好的区分疲劳程度的能力。MNF与fApEn仅在疲劳程度大或疲劳差异较大时具有较好的区分能力。而RQA%DET 在不同负载下具有最小的Lmax平均值，说明RQA%DET在不同负载下引起的肌肉疲劳程度区分性最差。

5 结论

通过采集肱二头肌的表面肌电信号计算了5种常见肌肉疲劳评价算法，从评价算法的抗干扰性与区分疲劳程度的能力两个方面量化其疲劳评价性能，以R2大小体现抗干扰能力的强弱，以Lmax展现区分疲劳程度的能力。实验结果显示，频谱距SMR同时拥有最大的R2和Lmax值，即频谱距SMR在评价肘关节肌肉疲劳时具有最好抗干扰和区分疲劳程度的能力。故频谱距SMR能够更加及时有效的检测出肘关节肌肉疲劳状态。因此在今后研究肘关节屈曲运动的肌肉疲劳时，建议使用频谱矩（SMR）进行实时监测。