1. 引言
表面肌电(sEMG)信号是人体的生物电信号之一,在肌肉接收到神经系统发出的运动指令后会执行相应的肢体运动 [1]。肌肉收缩时,细胞内外的带电离子发生相应的传输,运动神经元细胞产生一个动作电位由轴突传至肌纤维,sEMG信号就是由多个单纤维动作电位传到皮肤表面形成的 [2]。分析sEMG信号的特征值是临床医学、运动医学和人机工程学领域评价肌肉功能的重要的手段,也是研究肌肉疲劳状态的重要依据。康复病人在恢复期间容易运动过度导致肌肉损伤,运动员训练过程中训练过量也会导致肌肉拉伤。据统计,在大学生运动损伤案例中,肌肉损伤的病例就高达40% [3],因此,设计一款可穿戴的运动疲劳监测装置预防肌肉损伤十分有必要。近年来,sEMG信号较多地应用于肌肉疲劳评估中 [4] [5] [6]。李柽安等设计了一款可穿戴腰部运动损伤防护智能服装 [7];叶科学等研制了一款基于sEMG信号的椎旁肌监测系统 [8]。这些研究中,sEMG信号的采集和特征值分析是其主要基础。本文设计了一款基于肌电特征K值的穿戴式膝关节运动疲劳的监测系统。该系统通过实时采集人体大腿股四头肌部位的sEMG信号,并对其进行时域和非线性特征参数分析,提出了肌电特征 值来表征肌肉疲劳程度,以实现实时监测人体膝关节部位肌肉疲劳状态。该装置的研发对运动员和康复病人预防肌肉损伤具有重要意义。
2. 系统硬件设计
sEMG信号的幅值范围是0~5 mV,频谱主要分布在10~500 Hz,是一种十分微弱的生理电信号,易受皮肤表面静电、外界无线电磁的干扰,因此,要设计合理的sEMG信号采集电路来采集该信号,可穿戴式膝关节运动疲劳的监测系统的总体框图如图1所示,主要包括:无线穿戴式信号采集装置、手持式数据处理与预警装置。sEMG信号随电极进入无线穿戴式信号采集装置,主要包括:sEMG信号采集与预处理模块、微处理器、无线数据传输模块、智能电源管理模块。手持式数据处理与预警装置包括:无线数据传输模块、基于特征参数K值的信号处理模块、声光预警电路、电源模块。

Figure 1. Overall block diagram of wearable knee joint exercise fatigue monitoring system
图1. 穿戴式膝关节运动疲劳的监测系统的总体框图
本设计中,采集电极选择氯化银贴片电极,相比针电极,贴片电极适于测量运动时的肌电变化,具有无创性、简便性等特点。电极布局采用三点式差分方式。
sEMG信号采集与预处理模块完成对sEMG信号的拾取、放大、滤波和模数转换。主要功能模块包括:前置差分放大电路、带通滤波电路、后级放大电路、模数转换电路。其中,前置差分放大电路,去除sEMG信号中的共模噪声干扰,放大倍数为20倍;带通滤波电路,主要滤除10~500 Hz之外的噪声;后级可编程放大电路,可根据信号强度自适应调节放大倍数;模数转换电路采用双极性的16位采集系统,与微处理器之间通过SPI串行通信 [9],采样频率为2000 Hz。
智能电源管理模块,由锂电池供电,为整个系统提供±3.3V、±5V电压。为提高无线穿戴式信号采集装置的续航能力,该模块采用模拟开关切换控制低功耗待机与工作模式,手持式数据处理与预警装置端按下开始按钮时,无线传输模块将该信息发送到无线穿戴式信号采集装置的微处理器中,控制连接MOS管的IO口推挽输出高电平,无线穿戴式信号采集装置由低功耗模式进入工作模式;手持式数据处理与预警装置端按下结束按钮时,通过同样方法控制连接MOS管的IO口推免输出低电平,无线穿戴式信号采集装置恢复低功耗模式。
无线数据传输 [10] 模块采用CC101芯片,将肌电数据高效稳定地无线传输到基于特征参数K值的信号处理模块,进行滤波降噪处理和特征参数计算,并智能识别肌肉疲劳状态,在疲劳发生时控制声光预警电路进行提示。CC1101芯片是一款为低功耗无线应用而设计的低成本UHF收发器,射频收发器集成了一个高度可配置的基带调制解调器。该调制解调器支持各种调制格式,数据传输率可达500 kbps。微处理器的SPI接口连接CC1101无线模块来实现sEMG数据的无线传输,电路设计部分如图2所示。发送端接收无线穿戴式信号采集装置采集到的肌电数据,并将数据无线传输至手持式数据处理模块,供微处理器进一步处理。经过测试,采用CC1101无线模块进行数据传输与普通串口数据线相比,不存在多收、误收和数据丢失的情况。

Figure 2. Circuit design of CC1101 wireless transmission module
图2. CC1101无线传输模块电路设计
3. 基于样本熵与经验模态分解的去噪算法设计
3.1. 算法原理分析
相比于静止状态下的肌电信号采集而言,在运动过程中所采集的肌电信号更容易受到电极移位、动作摆动的干扰,因此,需要对采集的信号进行去噪预处理。
经验模态分解(EMD)算法是一种能够将原始信号分解为多个本征模态函数(IMF),适用于处理非线性、非平稳信号的新方法。EMD算法的基本思想是把复合信号分解为一系列IMF信号的组合,然后对各分量IMF进行希尔伯特变换进而得到希尔伯特谱,这种方法只与信号的本质特征相关。EMD方法具有自适应时频分析特征。
样本熵是近似熵的改进算法 [11],其大小反映了信号的复杂程度,其值越大表明信号越复杂,所含噪声越多,样本熵越小则表明信号的自相关程度高,所含噪声少 [12]。样本熵非常适用于生物信号的分析处理中。
小波变换去噪方法是分析非线性信号的常用方法之一,具有时频局部化与多分辨率分析的特点。其去噪基本原理是将含有噪声的原始信号进行多尺度小波变换得到小波系数和尺度系数,对不同尺度上的信号和噪声,根据小波系数性质进行去噪。
在肌电信号预处理去噪中,采用了基于样本熵的EMD算法 [13],其基本原理是将原始信号首先进行经验模态分解得到多个IMF分量,再对IMF分量计算样本熵,对样本熵较大的IMF分量进行小波去噪后重构得到去噪信号。
3.2. 算法实现流程
本文中应用基于样本熵与EMD算法对sEMG信号去噪流程如图3所示。

Figure 3. Denoising algorithm flow based on sample entropy and EMD
图3. 基于样本熵与经验模态分解的去噪算法流程
第一步,输入原始sEMG信号序列
,确定肌电信号序列的所有局部极值点,用三次样条线分别将所有的局部极大值点、极小值点连接起来形成上、下包络线,上下包络线平均值组成的序列记为
,得到:
(1)
式中,
为原始肌电信号序列,y为不断更新的原始肌电信号序列的分量序列,
为y序列与其上下包络线平均值组成序列
之间差值。
第二步,如果
满足如下条件:整个y序列的极点个数和过零点数相同或相差不超过一个,且整个序列任意一点处,由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值必须为零,则
为
的本征模态函数(IMF)的第一个分量。
如果
不满足上述条件,则将
的值赋给y,重复第一步,得到新的
;若
仍不满足条件,重复循环i次,得到
使得
满足IMF条件,记
,则
为原始信号
的第一个满足IMF条件的分量。
将
从y中分离出来,得到残差函数r:
(2)
第三步,对r进行上述一、二步计算,重复循环j次,直到r满足终止条件:即不断更新的r成为一个单调函数或常函数,或者只包含一个极值点,无法进一步分解时结束循环。即:
(3)
此时得到j个IMF序列
和一个残余分量
;其中,j个IMF序列代表原始肌电信号在不同频段的成分。然后计算每一个IMF序列
对应的样本熵
,对样本熵较大的IMF序列进行去噪,即满足:
去噪条件,则对相应的
进行小波去噪得到新的序列
,
为满足小波去噪条件的IMF序列个数;将不满足条件的IMF序列保留组成新的序列
。
最后,将残余分量、
序列和
序列重构得到去噪后的信号
:
(4)
其中,
为原始肌电序列去噪后的信号,
为经过小波去噪的IMF分量之和,
为保留未处理的IMF分量之和,
为经过EMD算法的残余分量。
4. 肌电特征K值的特征提取算法
sEMG信号的时域特征参数主要包括肌电积分值(IEMG)、均方根(RMS)。IEMG能够在时间维度上反映sEMG信号幅值的变化,当肌肉收缩时,IEMG的变化与肌力、肌张力呈正相关。其计算方法如式(5)所示。
(5)
其中,
代表sEMG信号序列的每一个数据,N为采样点数。
RMS主要反映了肌肉活动时运动单位激活的数量以及同步程度,它的变化主要表征sEMG信号的有效值。其计算方法如式(6)。
(6)
其中,
为sEMG信号序列的平均值。
样本熵S是sEMG信号的非线性特征参数,它除了反映人体肌肉活动的内在非线性特性外,还可以表征特定的肌肉活动,对于研究肌疲劳有重要意义 [14]。具体计算流程如下:
对于一个长度为N的肌电信号时间序列
,取一个正整数m满足
,则在U中存在
个长度为m的连续片段,其中第i个片段为:
(7)
对于给定的m,计算所有的m长连续片段之间的距离为:
(8)
给定阈值r,对于每个i计算其余
个m长连续片段与
的距离小于r的数目以及在所有片段中所占的比重,分别记为
,
,然后对所有的i,计算平均值:
(9)
所以,此肌电信号的样本熵为:
(10)
在本设计中,结合非线性特征参数样本熵S与时域特征参数均方根RMS来表征肌肉疲劳状态,引入了肌电特征参数K值。对于一个长度为N的肌电信号序列,考虑对所有数据处理,设定滑动窗口的宽度为与步长为
,其中N为
的倍数,得到
,其K值计算公式如式:
, (
) (11)
肌肉疲劳程度的评估与个体间的差异,肌肉的质量、锻炼程度相关,为了具体定义疲劳等级,我们引入疲劳因子
,计算公式如下:
(12)
其中,
为受试者正常初始状态下的特征参数,
为受试者随疲劳程度加深的实时特征参数。通过受试者与自身初始状态相比,随肌肉疲劳程度不断加深,
不同取值对应了不同疲劳等级,当疲劳因子等于0时,对应受试者非疲劳状态;当疲劳因子等于1时,对应受试者疲劳状态;当疲劳因子等于2时,对应受试者达到了极度疲劳状态,这时应当提示受试者停止运动、调整休息。
5. 实验验证与数据分析
5.1. 膝关节运动疲劳实验的方案设计
sEMG信号采集设备实物图如图4,图4(a)为sEMG信号采集硬件实物图;图4(b)为sEMG信号采集整体装配图,图4(b)中“测量电极”分别是正极、负极和参考接地电极。“电源开关”为整个设备的供电部位,内置2颗可充电纽扣电池,电池型号为常用的LIR2032,每颗电池容量约为40 mAh,额定电压3.6V。“采集按钮”的k0控制装置开始采集数据,k1控制装置复位操作。
选择(23 ± 3)岁、身高(170 ± 8) cm、体重在(60 ± 10) kg、24 h内无剧烈运动的男性大学生进行肌疲劳实验。在膝关节肌肉疲劳实验中,设计并完成了不同运动方式下的膝关节肌疲劳实验,包括:深蹲运动疲劳实验、坐姿下反复抬腿运动疲劳实验、模拟骑行阻尼自行车运动疲劳实验。具体实验动作如图5所示。贴电极片之前用酒精棉等清洁和湿润皮肤表面,预防皮肤表面静电干扰。考虑到股四头肌位于大腿前侧,是维持人体直立姿势和行走的重要肌肉之一,并且股直肌面积较大,方便采集数据,因此实验选取人体右腿的股直肌为信号采集部位。目前临床上尚未有定量指标对肌疲劳进行评估,主观疲劳量表(RPE)是一种在国内外被广泛采用的肌疲劳程度评估手段 [15] [16] [17] [18] [19]。本文参考RPE量表将疲劳程度等级分为三级,当肌肉处于放松状态时定义为正常状态;当肌肉开始疲劳时定义为疲劳状态;当肌肉非常疲劳时定义为极度疲劳状态。实验开始前设置采样频率2000 Hz,依次采集受试者在正常、疲劳和极度疲劳三种状态下的连续sEMG数据,每个状态下采集7000个采样点,采集时长为3.5 s。深蹲肌肉疲劳实验具体方案如下:首先让受试者保持站立姿势,此时定义为初始未疲劳状态,采集并记录正常状态下的sEMG数据;之后让受试者进行深蹲训练,规定3 s完成一组深蹲动作,当受试者刚开始感受到疲劳时定义为疲劳状态,采集并记录这一时刻的sEMG数据;让受试者继续深蹲训练,直至受试者感到极度疲劳并无法坚持时定义为极度疲劳状态,采集并记录这一时刻的sEMG数据。让受试者休息2小时,确保受试者恢复正常状态,重复以上步骤。每位受试者重复进行三次肌疲劳实验,共采集了10位受试者的实验数据。
(a)
(b)
Figure 4. Physical map of sEMG signal acquisition equipment: (a) Hardware physical map; (b) Whole assembly drawing
图4. sEMG信号采集设备实物图。(a) 硬件实物图;(b) 整体装配图
(a)
(b)
(c)
Figure 5. Different exercise muscle fatigue experiment. (a) Full squat repeatedly raises the leg; (b) Sitting posture; (c) Simulated riding damped bikes
图5. 不同运动的肌肉疲劳实验过程:(a) 深蹲;(b) 坐姿下反复抬腿;(c) 模拟骑行阻尼自行车
5.2. 肌电去噪算法的实验验证
通过sEMG硬件采集电路收集的肌电数据仍然存在工频噪声、基线漂移和高斯白噪声,即
,其中
为包含噪声的肌电信号,
为去噪后信号,
为工频干扰,
为基线漂移噪声,
为高斯白噪声。采用了24阶的巴特沃斯型数字滤波器进一步去除50 Hz工频干扰,其去噪效果如图6所示。与原始信号相比,去除了50 Hz左右的工频干扰。

Figure 6. The original waveform amd the power frequency denoising waveform of sEMG: (a) Time domain: Raw EMG signal; (b) Frequency domain: Original EMG signal; (c) Time domain: Filtered EMG signal; (d) Frequency domain: Filtered EMG signal
图6. 肌电信号的原始波形及工频去噪后波形:(a) 时域:原始肌电信号;(b) 频域:原始肌电信号;(c) 时域:工频、带通滤波后的肌电信号;(d) 频域:工频、带通滤波后的肌电信号
经过工频去噪后,使用了两种方法对sEMG信号进行降噪,去除基线漂移和高斯白噪声。第一种方法是小波阈值去噪方法,基本原理是将信号通过Matlab算法处理后,选择产生的小波系数 [20]。信号经过小波分解的小波系数大于噪声经过小波分解的小波系数,选择一个小波阈值,保留大于小波阈值的小波系数信号,将小于小波阈值的小波系数信号置零达到去噪的效果 [21]。第二种方法是基于样本熵的经验模态分解去噪方法。由于实际采集的sEMG数据本身包含各种噪声,无法作为“干净”信号计算信噪比,因此选择NinaPro肌电数据库2中的5位健康志愿者手部运动时采集到的sEMG信号进行去噪效果的分析。为了定量分析去噪效果,添加信噪比为10 db模拟高斯白噪声和基线漂移噪声,并分别计算两种去噪算法后获得的信噪比(SNR)和均方误差(RMSE),如式(13)、(14)所示:
(13)
(14)
其中,
是NinaPro肌电数据库2中的“干净”sEMG信号,
为去噪后的sEMG信号,N为数据长度;SNR表示去噪信号与噪声的能量之比,反映了信号噪声的去除程度,RMSE反映了信号降噪前后的偏离程度。当SNR越大,RMSE越小时,则表明一个算法的去噪效果越好。
去噪效果如表1,结果表明:经过样本熵的经验模态分解去噪方法,去噪后的sEMG信号SNR更大,RMSE更小,去噪效果更好。最终选择基于样本熵的经验模态分解去噪方法对实验采集的sEMG数据去噪。

Table 1. Ninapro EMG database two data used two methods to compare the denoising effects
表1. NinaPro肌电数据库2数据使用两种方法去噪效果对比
5.3. 肌电特征提取算法的实验验证
本系统按照sEMG信号的基本分析方法,首先对深蹲实验采集的sEMG信号进行时域指标与非线性指标分析。本文以500个采样点作为步长,500个点为窗口函数,进行50%重叠率的滑动计算,计算得到RMS、样本熵S和K值。受试者a肌疲劳前后的特征值变化结果如图7所示,横坐标代表窗口数,三条分割线划分了正常、疲劳和极度疲劳三种状态下的特征值变化结果,时域特征值RMS随着疲劳程度的加深而递增,非线性特征值样本熵S随着疲劳程度的加深而递减,特征参数K值随着疲劳程度的加深而递增,且递增趋势相较于RMS更加明显。

Figure 7. Results before and after fatigue characteristic value for subject
图7. 受试者a疲劳前后特征值变化结果
对受试者a的三种特征参数进行配对t检验 [22],结果如表2,对比正常状态–疲劳状态、疲劳状态–极度疲劳状态,可知:K值肌疲劳前后的显著性P值最小,说明K值肌疲劳前后的显著差异性最大,更容易区分肌肉不同疲劳状态,适合作为评估肌疲劳程度的特征参数。

Table 2. Paired t-test results of the characteristic parameters for subject
表2. 受试者a的特征参数的配对t检验
将10位受试者深蹲实验的特征指标参数取平均值,整理肌疲劳前后所有特征参数变化柱状图如图8,从正常–疲劳–极度疲劳这一过程中,时域特征指标RMS呈递增变化,在疲劳–极度疲劳段变化趋势较小;从正常–疲劳–极度疲劳这一过程中,非线性特征指标样本熵S呈递减变化,但整体变化趋势较小;从正常–疲劳–极度疲劳这一过程中,K值呈递增变化,变化趋势最为明显,衡量肌肉疲劳程度的效果优于样本熵S和RMS,更适合作为评估膝关节运动疲劳的特征指标参数。
(a)
(b)
(c)
Figure 8. The variation of time, nonlinear and K-value before and after squat muscle fatigue of 10 subjects: (a) Time domain characteristic index RMS changes before and after muscle fatigue; (b) The nonlinear characteristic index S changes before and after fatigue; (c) Changes of K value before and after muscle fatigue
图8. 10位受试者时、非线性和K值特征指标深蹲肌疲劳前后变化情况:(a) 时域特征指标RMS肌疲劳前后变化情况;(b) 非线性特征指标样本熵S肌疲劳前后变化情况;(c) K值特征指标肌疲劳前后变化情况
为了验证算法的有效性和泛化性,我们针对受试者不同动作,如深蹲、坐姿下反复抬腿和模拟骑行阻尼自行车。进一步分析肌电特征K值在不同疲劳状态下的变化情况如图9所示,三种实验动作下的K
(a)
(b)
(c)
Figure 9. K-value of 10 subjects before and after muscle fatigue caused by different movements: (a) Change of K-value before and after muscle fatigue caused by squat; (b) Changes in K-value before and after muscle fatigue caused by repeated leg lifting under sitting posture; (c) Changes of K-value before and after muscle fatigue caused by simulated cycling damping bicycle
图9. 10位受试者不同动作方式导致肌肉疲劳前后的K值变化:(a) 深蹲导致肌肉疲劳前后的K值变化;(b)坐姿下反复抬腿导致肌肉疲劳前后的K值变化;(c) 模拟骑行阻尼自行车导致肌肉疲劳前后的K值变化
值变化表现出一致规律性,即随着疲劳程度的加深,K值呈递增趋势。其中,坐姿下反复抬腿实验组的K值整体幅值变化相较于其他两组数据较小,对其分析发现:坐姿下反复抬腿过程中,运动强度低于其他两种运动方式,使得运动肌群的负载相对较低,导致受试者在整个实验过程中时域特征值RMS变化趋势较小,根据公式
可知肌电特征K值的变化幅值也相对较小。综上所述,不同运动的膝关节肌肉疲劳对应的肌电特征K值变化呈现出相似的变化规律,我们提出的算法可以适用不同运动方式下的肌肉疲劳识别。
5.4. 系统实时性分析
将Matlab软件实现的滤波和特征值处理算法移植到微处理器上,编程控制系统实时处理sEMG信号,并在原设备基础上增加HMI串口显示屏,系统的功能样机如图10所示。整个设备工作流程:首先进行
500个采样点的数据采集,采样频率为2000 Hz;在每次数据采集完成之后,微处理器进行预处理去噪、特征值计算和HMI串口屏显示任务,在此期间不进行数据采集,相邻数据采集的时间间隔设置为1.20 s。通过实验测试了微处理器每次处理500个肌电数据所需时间,包括数据采集、数据预处理、特征值计算和HMI串口屏显示任务,共需0.98 s时长,如表3所示。当前的数据采集间隔能满足运动过程中膝关节肌疲劳监测的实时性需求。

Table 3. The time consumed by the system for one round of data sampling and processing
表3. 系统每轮采样与数据处理所需时间
6. 结论
本文提出了一种基于肌电特征K值的肌疲劳评估方法,设计了基于sEMG信号的穿戴式膝关节运动疲劳监测系统的样机。通过采集人体膝关节股四头肌在疲劳前后的sEMG信号,分析sEMG信号的时域、非线性特征值,设计肌电特征K值作为评估肌疲劳程度的指标参数。通过对疲劳前后的三个特征参数:样本熵、均方根和肌电特征K值分别进行配对t检验,结果表明:肌电特征K值在肌疲劳前后的显著差异性P值最小,区分肌肉疲劳与非疲劳状态效果更加明显,适合作为评估膝关节运动疲劳状态的特征参数。通过实验验证,我们提出的肌电特征K值算法能够适用于三种不同运动方式下的膝关节肌肉疲劳识别,K值在肌疲劳前后表现出相似的变化趋势,该算法的泛化能力良好。所研发的穿戴式膝关节运动疲劳的监测系统可以实现运动过程中的肌电采集与肌疲劳评估,进而为预防肌肉运动损伤提供参考。
后期,拟采用近红外肌氧监测仪监测待测肌群的肌氧饱和度(SaO2) [23],并结合主观疲劳量表RPE评分对肌肉状态进行评估,增加实验数据提高肌疲劳评估的准确度。并且,进一步使系统样机小型化、集成化,将其应用于临床,积累更多的临床数据,为实现人工智能评估肌疲劳提供依据。
基金项目
上海市自然科学基金面上项目(20ZR1437700);
上海市2022年度生物医药科技支撑专项(22S31902200);
上海市市级科技重大专项Shanghai Municipal Science and Technology Major Project(2021SHZDZX);
上海市“科技创新行动计划”(22S31902200);
上海市普陀区特色专科项目(2019tszk01)。
NOTES
*通讯作者。