一种基于肌动信号的股四头肌收缩力量估计方法研究*

王大庆,郭伟斌,吴海峰,高理富*

(1.中国科学院合肥智能机械研究所,合肥 230031;2.中国科学技术大学,合肥 230026;3.安徽省仿生感知与先进机器人技术重点实验室,合肥 230031)

随着老龄化社会以及机器人与人越来越广泛的共存共融,可穿戴机器人系统日渐成为机器人领域的研究热点。简便实时准确获取人体运动意图及相关参数是实现舒适、自然、灵活的人机协调运动的关键[1]。肌肉收缩力量作为肌肉活动的主要特征,在一定程度上反映了肢体的相关运动信息和疲劳状态,研究其测量方法,可将其用于可穿戴设备的控制或者监测人体肌肉的活动状态,对于可穿戴设备、康复训练以及科学运动等领域具有重要的意义。

由于肌肉被皮肤、脂肪等组织包裹,目前还难以在活体上直接测量肌肉发力大小。传统运动生物力学多借助于三维测力台和等速测试设备(如Cybex测力仪)进行肌肉收缩力量的测量,应用场合有限。而基于运动学和动力学分析方法的肌肉收缩力量估计模型需要采集准确的运动学参数,实验设计复杂,成本较高。为了克服上述缺点,研究人员尝试利用表面肌电信号sEMG(Surface Electromyography)来实现人体肌肉收缩力量的非侵入性测量与评估,取得了很好的效果[2-5]。但是肌电信号较微弱,易受其他电子设备的干扰,对于采集系统和测量环境要求较高;肌电测量电极大多需要放置在皮肤表面,易受皮肤表面的阻抗变化的影响[6],这些限制因素导致sEMG的实际应用不够灵活、方便。

研究表明,肌肉活动时不仅会伴随有电信号sEMG产生,同时还会产生低频的机械振动信号MMG(Mechanomyography),中文译为肌声[7],肌动[8]或者肌音[9]信号。MMG信号记录的是肌肉在主动收缩时横向振动的力学现象,反映了肌肉运动过程中的力学特性。[8-9]。

和sEMG相比,MMG具有采集系统成本低,抗干扰性好的特点。更为重要的,由于其机械振动信号的特性,MMG信号对于传感器的摆放位置不敏感。此外,采集时不需要直接接触皮肤表面,隔着衣物依然可以有效采集,从而不受皮肤表面汗液,温度引起的阻抗变化的影响[10],尤其适合应用于可穿戴设备。

MMG信号与肌肉力之间的关系大致为:在0～80%MVC(MVC是maximum voluntary contraction的缩写,意为最大肌肉随意收缩力。由于每个人的肌肉力量大小有差异,所以在表示肌肉力量时,通常以百分比的形式描述),MMG信号幅度随肌肉力的增大而增大;在80%～100%MVC,根据不同的情况,信号幅度会继续增大或者下降;并且会受到肌肉运动单元的激活策略,肌肉纤维的构成和肌肉强壮程度等因素的影响[11]。因此,MMG信号与肌肉力之间是比较复杂的非线性关系。此外,研究人员针对不同的应用背景,从时域、频域等各个不同的角度提出了不同的特征提取方法[12],但是寻找最能反映原始信号变化趋势的特征参数始终是一个难点。对于MMG-力关系模型的构建,传统的多项式回归算法在面对复杂的非线性关系时就显得不够灵活。因此,研究人员把目光投向了神经网络[13-15],但神经网络也存在着诸如泛化能力差、在训练样本量比较少时容易出现过拟合等问题。

鉴于此,本文针对可穿戴设备中的人机交互过程控制所需的肌肉力/力矩信息测量需求,选取人体下肢最大、膝关节活动最重要的股四头肌作为研究对象,利用加速度传感器采集股四头肌活动时的MMG信号,从时域、频域、熵和信号相关性角度提取若干特征,最后通过相关向量机算法建立了MMG信号-力的关系模型,并通过实验验证了模型对股四头肌收缩力量进行估计的有效性。

1 信号获取与特征提取

1.1 MMG信号的获取

在之前的研究中,许多不同类型的传感器被用来检测MMG信号,其中包括麦克风传感器(MIC)、压电传感器(PIZ)、加速度传感器(ACC)和激光测距传感器等。其中,加速度传感器以其体积小、重量轻,便于安装固定而备受青睐;并且对于人体较大块的肌肉,从一个测量点采集的MMG信号不能代表肌肉整体的活动情况[16],而基于MEMS工艺的ACC传感器成本较低适合组建多通道采集甚至传感器阵列。

图1 有效信号的提取(股内侧肌)

MMG信号本质上是由肌肉纤维收缩所引起的一种机械振动信号,其有用信息主要集中在低频段[17]。通过ACC传感器采集的信号中,除高频噪声外,还会掺杂重力加速度信息和由于肌肉收缩引起的低频运动伪迹如图1所示。

此前,有研究人员仅利用原始的ACC信号进行相关动作的识别[18]。但其信号受传感器的摆放位置影响较大,以致每次测量信号波形差别较大,使得准确度难以保证。因此必须对原始加速度信号进行滤波处理。如图1下图所示,本文采用了通带为5 Hz～40 Hz的4阶巴特沃斯带通滤波器,只保留由肌肉纤维振动所引起的低频MMG信号,提高测量系统的对于传感器摆放位置的鲁棒性。滤波后的MMG信号幅值与力传感器所测力值呈较明显的正相关。

1.2 MMG信号特征提取

MMG信号的特征提取方法因不同的研究背景和目的而异[12]。从图1可知,肌肉力量与MMG信号的幅值密切相关[12]。因此,本文首先计算了描述信号幅值变化的平均绝对值MAV(Mean Absolute Value):

(1)

此外,研究表明肌肉产生的张力与同时参与收缩的肌纤维数目有关。同时兴奋的运动单位越多,参与收缩的纤维数目就越多,产生的力量也越大。那么,随着参与运动单位数量的增多,MMG信号中包含的成分理应增多,换句话说,信号会变得更复杂。本文使用了样本熵SampEn(Sample entropy)来衡量这种复杂度的变化,

SampEn(m,r,N)=-ln[Bm+1(r)/Bm(r)]

(2)

式中:m为嵌入维数,r为阈值,N为序列长度,具体计算过程详见文献[19]。SampEn值越小,代表序列自我相似性越高,值越大,序列越复杂。值得注意的是,阈值r的选择对于SampEn值有着很大的影响。通常,r取序列的0.1至0.2个标准差。但在本文中,r取固定值(例如0.025),即r的取值不随滑动窗口序列改变时,SampEn值与肌肉收缩力量呈现较强的相关性。

另外,对于肌肉收缩力量的估计,频率特征是否有用以及来自不同肌肉的MMG信号之间的串扰是否存在,学术界存在争议。

本文通过计算每个通道的平均功率频率MPF,在特征中加入了频率信息。此外,计算不同2个采集通道信号之间的相关系数CC2Cs(Correlation Coefficient between 2 different Channels)(式(3)),将其作为评价信号之间串扰情况的特征。

(3)

算法表明,在回归建模中,加入频率信息和不同通道信号之间的相关性信息可以提高模型估计的准确性。

2 相关向量机回归算法

tn=y(xn;ω)+εn

(4)

式中:εn服从均值为零,方差为σ2的高斯分布。则

p(tn|x)=N (tn|y(xn),σ2)

(5)

上式均值y(xn)由以下线性模型给出

(6)

式中:ω0是一个偏置参数,其形式与SVM类似,但基函数由核给出,训练集的每个数据点关联着一个核,且没有正定核的限制。

根据tn的独立观测假设,似然函数可表示为

(7)

为了避免出现过拟合,假定参数向量ω服从零均值的高斯先验分布,即

(8)

RVM的关键在于为每个权参数ωi都单独引入一个超参数αi,而不是一个共享的超参数,其中αi即为对应权值ωi的超参数。而α与σ2服从Gamma分布,通过将其分布的形状参数和尺度参数设为零,可以得到一致超先验,可得到p(ωi)。根据贝叶斯公式可得到权值的后验概率分布仍然是高斯分布:

p(ω|t,α,σ2)=N (ω|m,Σ)
m=σ2ΣΦTt
Σ=(A+βΦTΦ)-1

(9)

式中:Φ是N×M的设计矩阵,元素为Φni=φi(xn),i=1,…,N且ΦnM=1,n=1,…,N,A=diag(αi)。

然后通过最大化边缘似然函数来确定α和σ2的值。边缘似然函数通过对权向量积分的方式得到

(10)

由于这表示两个高斯分布的卷积,因此可以计算对数边缘似然函数,形式为

(11)

式中:t=(t1,…,tN)T,并且定义了N×N的矩阵C,形式为

C=(σ2)-1I+ΦA-1ΦT

(12)

接下来令似然函数(9)导数等于零来求解超参数α和σ2,可得

(13)

式中:mi是式(9)定义的后验均值m的第i个分量。γi度量了对应的参数ωi由数据确定的效果,定义为

γi=1-αiΣii

(14)

式中:Σii是式(9)给出的后验协方差Σ的第i个对角元素。学习过程就是利用式(12)不断的调整超参数,然后再利用式(9)重新计算后验概率的均值和协方差,直到满足一个合适的收敛准则。

作为最优化的结果,很大一部分αi都趋于无穷,对应的权参数ωi的后验概率的均值和方差都是零。也就是说,与这些参数关联的基函数φi(x)对于模型的预测不起作用,因此被高效的剔除,从而生成了一个稀疏的模型。对应于剩下的非零权值的输入xn被称为相关向(relevance vector),因为它们是通过自动相关性检测的方法得到的。

找到了最大化边缘似然函数的超参数α*和σ2*之后,对于一个新的输入x*,就可以计算t上的预测分布。由式(5)和式(9)可得:

=N [t|mTφ(x),σ2(x)]

(15)

因此预测均值由式(5)给出,其中ω被设置为后验均值m,预测分布的方差为

σ2(x)=(σ2*)-1+φ(x)TΣφ(x)

(16)

3 基于肌动信号的人体股四头肌收缩力量估计

为了能够准确的得到肌肉收缩过程中的力信息,本文采用了等长收缩的方式进行数据采集。基于肌动信号MMG的人体股四头肌收缩力量估计方法首先记录股四头肌活动时的MMG信号和相应的力信息,然后将滤波、特征提取后MMG信号与以%MVC表示的力信息作为输入,通过相关向量机算法训练得到MMG-力回归模型,最终实现对股四头肌收缩强度的连续估计(如图2)。

图2 肌肉收缩力量估计方法框图

等长收缩实验时,参与者采取坐姿状态,右腿与装置固定在一起(图3)。膝关节弯曲近似90°,小腿前侧与力传感器接触,用于采集实验过程中,膝关节所产生的伸小腿力信息。根据解剖学知识,小腿伸主要依靠的是股四头肌收缩,故此处测得的力可以反映股四头肌收缩力量的相对大小。

图3 实验装置及股四头肌的组成

实验时,记录股直肌RF(rectus femoris)、股内侧肌VM(vastus medialis)和股外侧肌VL(Vastus Lateralis)的MMG信号(股四头肌中的股中间肌位于股直肌下方,无法直接测得)。

MMG传感器采用了由ADXL327芯片(3轴,测量范围±2个重力加速度)制成的加速度传感器(尺寸27 mm×22 mm,重2.89 g)。实验时3个传感器放置于衣物表面,用绑带加以固定(图4)。

图4 加速度传感器及固定示意图

7名肢体健康的测试者(6男1女,身高(171.8±6.23)cm,体重(64.0±13.2)kg)参与实验。实验开始前会测量每名测试者在该姿势下的最大肌肉随意收缩力(MVC),实验中会将测量的力数值转换成其百分比(%MVC)。每名测试者完成4组实验数据的采集,每组实验持续120 s,测试者依据屏幕的实时力反馈,完成近似正弦波形的发力曲线。所有信号的采集使用数据采集卡NI USB-6215,采样频率设为500 Hz。

由于所采用的是一款3轴加速度传感器,从理论上说,垂直于皮肤表面的分量信号会比其他2个分量信号强,但为了消除每次传感器摆放位置带来的差别,以及实验过程中传感器因肌肉收缩引起的角度变化,最终取3个分量信号的合加速度作为MMG信号的原始信号。

(17)

原始MMG信号经4阶巴特沃斯带通(5 Hz～40 Hz)滤波后进行特征提取。设置1 000 ms的时间窗和250 ms的增量窗,用以提取MMG特征。为了增强特征提取方法的鲁棒性,本文尝试从不同的计算角度选取了MAV、MPF、SampEn和CC2Cs 4个特征。选取的特征与MMG原始信号之间的关系如图5所示。

图5 原始波形的特征表达

从图5可以看出MAV和SampEn可以很好地反映MMG信号幅值的变化,而MPF和CC2Cs也似乎蕴含了有用信息。

本文从3个通道的MMG信号中各提取3个特征及他们之间的相关系数,构成1个12维的输入向量,目标值t即为股四头肌收缩力量的相对大小(%MVC),由力传感器测得的数值转化而来。核函数选择高斯径向基(RBF)核函数,最大迭代次数设为1 000次。

图6 回归模型输入构成

每个参与者进行4组测试,所得数据集通过3次10折交叉验证[21]来测试算法的准确性。每次验证计算模型在测试集上的均方根误差RMSE(Root-mean-square error)和决定系数R2(coefficient of determination),用以表征模型估计结果的优劣。

4 实验结果与分析

4.1 不同特征组合所得估计结果对比

图7列出了在将不同特征组合作为输入的情况下,模型对同一个参与者的数据验证估计的结果。均方根误差RMSE反映了估计值同真实值之间的偏差,决定系数R2则反映了回归曲线的拟合程度,值越大,意味着自变量(MMG特征)对因变量(力)的解释程度越高。

图7 不同特征组合输入模型的准确度对比

从图7可以看出,多输入的估计效果一般优于单输入。尽管MPF和CC2Cs与原信号之间的相关性在图5中表现的不高,但如图7所示,通过引入这两个特征,加入频率信息和信号串扰特征,提高了回归估计的准确度。

此外,在信号特征提取时,均方根RMS(Root-Mean-Square)也常用来描述信号的幅值变化,过零点ZC(Zero Crossing)也可以在一定程度上反映信号频率的变化情况。但在本问题中,ZC的引入似乎并不能提高估计的精度,并且RMS和MAV刻画的均为原始信号幅值信息,最终只保留了特征MAV。因此,本文根据图7所得结果,MAV、MPF、SampEn和CC2Cs 4个特征作为估计模型的输入。

4.2 准确度结果统计

7名参与者的数据集单独验证的结果如表1。估计模型的平均RMSE为8.7%MVC,R2为0.817,估计准确度要优于文献[13](RMSE=14.1%MVC),接近文献[15](RMSE=8.6%MVC)。但本文所采用的测量方式比文献[15]更为灵活、简单。

与此同时,表1右侧列出了模型在不同参与者之间验证的准确性,具体的做法就是用6个参与者的数据训练模型,剩下1人的数据作为测试集验证模型的准确性。从表中可以看出,无论是RMSE还是R2都不如图1左侧的准确度高。这可能的原因一方面是样本数据太少,另一方面也是由于不同参与者之间的身体肌肉状况差异比较大。

表1 模型估计的准确度

图8展示了模型估计的结果。黑色曲线代表力传感器测得的真实力值,蓝色曲线是利用不同参与者数据训练构建的模型估计的结果,红色曲线是同一参与者自身数据建立的模型估计的结果。两个模型均能在一定程度上反映肌肉力量的变化趋势,但红色曲线的估计结果明显优于蓝色曲线。

图8 模型估计结果

4.3 肌肉收缩力量估计算法对比

文献[13-15]均采用了神经网络BPNN算法,并且SVM与RVM也有诸多相似之处,因此针对肌肉收缩力量估计问题,有必要将RVM算法与SVM、神经网络(以后向传播BPNN为例)进行对比,以验证RVM算法在此问题的有效性。

BPNN模型采用3层网络结构,其中输入层节点数12,隐含层节点数20,输出层节点数1。目标误差10-5,最大迭代次数1 000。SVM模型选用高斯核函数,利用网格搜索方法确定SVM参数惩罚因子和核函数参数的最优值[22]。所有运算均在MATLAB软件环境下完成。图9展示了SVM、BPNN和RVM针对同一参与者的数据训练建模并估计的结果(3次10折交叉验证)。

图9 BPNN、SVM与RVM估计结果对比

图9中数据显示,SVM与BPNN准确度接近,但BPNN估计结果的标准差稍大。RVM的准确度比BPNN提高4.4%,比SVM提高4.6%。此外,针对1860个样本点的回归问题,RVM得到的相关向量数量为139±2.16个,而SVM的支持向量为467±8.18个,这也体现了RVM在预测函数上的稀疏性。

如表2,RVM训练时间较长,但预测时间优于BPNN,SVM在运算时间上优势明显。由于RVM算法在训练中的每次迭代和预测时还需要矩阵求逆等运算,因此像文献[23]分类问题中描述的RVM训练时间和识别时间均优于SVM的情况并未出现。这可能是由于数据点规模不够大,无法发挥出RVM预测函数稀疏性的优势。

表2 RVM与SVM、BPNN运算时间对比

5 结论

本文通过采集人体股四头肌肌肉活动时的MMG信号,建立了MMG信号与肌肉收缩力量关系模型,实验结果验证了本文方法的有效性,并且相比sEMG,MMG的获取可以不用直接接触皮肤表面,测量方式更为便捷、灵活。同时,基于稀疏贝叶斯思想的相关向量机算法在此回归问题上表现出了比SVM和BP神经网络略优的准确性。下一步的研究,准备把本文的人体肌肉收缩力量的估计方法运用到研制的可穿戴助力机器人的交互控制中,检验本方法的有效性和适用性,实现更为灵活的人机协调运动。