基于Hill肌肉模型的人体关节力矩智能预测

2021-02-01熊保平史武翔林昱黄美兰杜民

北京生物医学工程 2021年1期

熊保平史武翔林昱黄美兰杜民,4

0 引言

人体关节力矩测量是当前生物医学领域最具挑战性的课题之一[1]。关节力矩值可应用于康复训练[2-3]、运动员的训练评估[4-5]、假肢及矫形器设计[6-7]等领域。

但由于人体结构和运动过程各肌肉协调的复杂性，还难以在活体上直接测量人体关节力矩值。目前有测力法[8-9]、反向生物力学分析法[10-11]和正向生物力学分析法[12-14]这三种主要的间接测量方法。但这些测量方法存在模型参数校正的时间复杂度高和个体误差等问题。而以数据为驱动的人工智能技术正好可以在训练过程中校正模型参数，适应个体特征且操作简便，使其越来越多地应用于人体关节力矩智能预测中。例如，Uchiyama等[15]利用肩与肘之间肌肉的综合肌电信号以及肩和肘关节角作为人工神经网络的输入项预测肘关节屈伸力矩大小。Luh等[16]利用肌电信号以及肘关节屈伸自由度的角度和角速度作为输入预测肘关节屈伸力矩值。Song 和Tong等[17]把肌电信号，以及肘关节屈伸自由度的角度和角速度作为人工神经网络的输入项预测肘关节屈伸力矩。在Hahn[18]的论文中甚至把年龄、性别、身高、体重以及肌电信号作为神经网络的输入预测膝关节的屈伸力矩。Ardestani等[19]提出利用力板数据和肌电信号作为输入训练人工神经网络预测的髋关节屈伸、旋转、收展以及膝关节屈伸等关节自由度力矩值。还有文献把杠铃的重量和运动轨迹作为神经网络的输入预测下肢髋、膝和踝等关节各个自由度力矩值[20]。通过上述文献可知由于人体运动过程中各肌肉协调的复杂性导致目前缺乏有效方法确定人工神经网络预测关节力矩的输入变量，从而使各方法预测关节力矩所使用的输入变量几乎不相同，而输入变量的特征不完整以及信息冗余都会影响人工神经网络预测结果的准确性[21]。

针对上述问题，本文根据Hill肌肉模型原理[22]结合人体几何学模型[23]建立人体关节力矩智能预测输入-输出关系的数学模型，确定了关节力矩智能预测的输入变量。在该数学模型中把Hill肌肉模型中不易测量的输入变量，如肌肉纤维的长度、收缩速度，以及以关节自由度为支点的肌肉力臂等在活体中不易测量的输入量，转化为肌肉所驱动关节自由度角度和角速度的函数，从而确定关节力矩智能预测的输入变量为该关节自由度所关联肌肉的肌电信号(electromyography,EMG)信号，以及这些肌肉所驱动关节自由度的角度和角速度。为了验证本文提出的关节力矩预测的输入变量的准确性和有效性，一种新型的快速的前馈神经网络学习算法[24-27]——极限学习机被应用于关节力矩的训练与预测。

1 材料与方法

1.1 实验数据

为了验证本文提出方法的有效性和可靠性，从网站(https://simtk.org/projects/emgdri-venmodel)获取一个右下肢偏瘫男性(年龄79岁,身高1.7 m,体重80.5 kg，LE Fugl-Meyer运动功能评估值32/34 pts)运动生物力学模型标定和随意步态下的下肢关节动态力矩预测的实验数据。实验中三维运动捕捉设备采用的是维康公司VICON光学运动捕捉系统，地面反力检测设备采用的是Bertec公司的跑步机，EMG信号的采集设备采用的是Motion Lab Systems公司的信号采集设备。实验采集了实验对象在跑步机上以0.4、0.5、0.6、0.7和0.8 m/s等5个速度行走的运动学数据、动力学和EMG数据。人体运动学数据的采集频率是100 Hz，动力学和EMG数据的采集频率是1 000 Hz。每个速度至少采集40个步态周期运动数据。运动捕捉的marker的安置位置参考文献[28]。地面反作用力和marker运动数据采用四阶零相位滞后巴特沃斯滤波器[29]。滤波器截止频率是一个以7除以步态周期的可变截止频率。

实验根据人体解剖学知识，确定16块对下肢关节运动有主要贡献的肌肉并获取其肌电信号，其中臀大肌、臀中肌、半膜肌、股二头肌长头、股直肌、股中肌、股外侧肌、内侧腓肠肌、胫前肌、腓骨长肌、比目鱼肌11块肌肉采集的是表面肌电信号。内收长肌、髂肌或腰大肌、胫后肌、趾长屈肌、趾长伸肌5块肌肉采用插入式的方式获取肌电信号。采集到的EMG信号首先进行截止频率为40 Hz的高通滤波去除直流和低频噪声，其次进行全波整流使信号非负，然后对信号进行归一化处理去除信号的个体差异。最后对所有运动学数据、肌电以及反向动力计算的力矩值等数据进行重取样，对每个步态周期采集101个数据点，以消除速度对步态分析的影响。详细的处理过程参见文献[13]。

1.2 变量的选择

由于输入变量的特征不完整以及信息冗余都会影响人工神经网络预测结果的准确性[21]，因此本文提出利用人体几何学和Hill肌肉模型来建立关节力矩预测的输入输出数学模型,从而确定智能预测的输入变量。其处理流程如图1所示。

图1 基于Hill肌肉模型的数据处理流程

由Hill肌肉模型[12]可知单块肌肉关于其驱动关节自由度的力矩可以表示为：

(1)

l(t)=a0+a1θ(t)+a2θ(t)2+a3θ(t)3=l(θ)

(2)

式中：l为肌肉纤维的长度；θ(t)为肌肉所驱动关节自由度的角度；a0～a3为多项式的常系数。

肌肉纤维的长度对时间t的一介导数即为肌肉的收缩速度v：

(3)

文献[32]发现由肌肉纤维的长度l对肌肉驱动关节自由度角度的一阶偏导数可得肌肉相对于该关节自由度的力臂，其表达式为：

(4)

式中的负号用于确定力矩的方向。当肌肉驱动两个关节自由度的时候，式(2)～式(4)可扩展为：

(5)

(6)

(7)

(8)

同样,当肌肉驱动3个或者更多个关节自由度的时候可做如上类似的扩展。

由式(1)～式(8)可知肌肉驱动的第i个关节自由度的力矩可以表示为：

(9)

由于所有驱动该关节自由度的肌肉力矩的矢量和即为该关节自由度的力矩，所以第j个关节自由度的力矩可表示为：

(10)

式中：m是第j个关节自由度所关联的肌肉总数。

由此可知,关节力矩智能预测的输入变量是关节自由度所关联肌肉的EMG信号,以及这些肌肉所驱动关节自由度的角度和角速度。

1.3 极限学习机

本文选择极限学习机(extreme learning machine,ELM)作为预测关节力矩的算法。极限学习机是一种新型的快速前馈神经网络学习算法[25]，它的结构如图2所示，分为输入层、隐含层和输出层。极限学习机有别于传统前馈神经网络的显著特点在于其输入层到隐含层的权值随机选取且在训练过程中无需调节，在此基础上隐含层到输出层的权值是通过求解矩阵的 Moore-Penrose 广义逆来唯一确定。这一特点使得网络权重系数的确定过程无需迭代,降低了网络权重系数的调整时间，从而大大提高了学习速率。目前ELM算法已广泛应用于模式识别、计算机视觉、数据挖掘、控制论等领域[33-34]。

图2 极限学习机结构

ELM根据式(10)确定的输入变量进行训练，并预测人体下肢的髋关节屈伸(flexion-extention,FE)、髋关节收展(adduction-abduction,AA)、膝关节屈伸、踝关节背屈(plantar-dorsiflexion,PDF)和踝关节内外翻(inversion-eversion,IE)等5个关节自由度的力矩。由于肌电信号只采集了16块肌肉，所以其他肌肉的肌电信号参考文献[14]，通过相邻肌肉的肌电信号可近似获得。由于相同的输入对于人工神经网络而言是信息冗余，所以本文ELM只把16个EMG信号作为肌电的输入。通过表1可以获得每个关节自由度力所关联的肌肉，以及这些肌肉驱动的关节自由度。例如髋关节屈伸根据表1可知由内收长肌、臀大肌、臀中肌、髂肌、半膜肌、股二头肌长头、股直肌7块肌肉驱动，而这7块肌肉同时又驱动了髋关节屈伸、髋关节收展、膝关节屈伸这三个关节自由度，所以根据式(10)的结论可知预测髋关节屈伸力矩的输入变量为收长肌、臀大肌、臀中肌、髂肌、半膜肌、股二头肌长头、股直肌7块肌肉的肌电信号和髋关节屈伸、髋关节内外翻展和膝关节屈伸这三个关节自由度的角度和角速度共计13个变量。其他关节自由度力矩预测的输入变量可以参照上例分析。

表1 肌肉信号源及其驱动关节自由度列表

2 结果

实验在两个不同的泛化水平下分析ELM的预测能力，它们分别是只把前3个低速(0.4、0.5和0.6 m/s)的步态周期数据用于ELM的训练(level 1)和把所有速度下的步态周期数据用于ELM的训练(level 2)并预测所有速度下的关节力矩值。实验随机选取每个速度下的40个步态周期数据中的30个用于训练，余下的10个步态周期的数据用于测试。实验结果的准确性是通过计算ELM预测的力矩值与反向生物力学分析计算的关节力矩值的归一化均方误差(normalized root mean squared error,NRMSE)和相关系数(ρ)的统计值来验证的。

在level 1泛化水平下各个步速关节自由度力矩的预测结果如图3、图4所示(以0.5 m/s步速为例)。由图可知，该泛化水平除右踝关节内外翻外，其他关节自由度的力矩都可以较精确地预测，且NRMSE<15.48%、ρ>0.86。

图3 以0.5 m/s速度行走的平均关节力矩预测结果(level 1)

图4 在5个不同速度下的各个关节自由度预测力矩的NRMSE值(level 1)

在level 2泛化水平下各个步速关节自由度力矩的预测结果如图5、图6所示(以0.5 m/s步速为例)。由图可知，该泛化水平同样除右踝关节内外翻外，其他关节自由度的力矩都可以较精确地预测，且NRMSE<12.93%、ρ>0.89。由于通过图2～图6无法精确地区分这两个泛化水平的预测精度，所以通过图7、图8可以发现level 2除在0.4 m/s和0.5 m/s这两个速度区分度不大外，其他各关节自由度力矩预测精度明显高于level 1。与相比level 1、level 2的力矩预测的平均NRMSE降低了9.4%，ρ也略有提高。

图5 以0.5 m/s速度行走的平均关节力矩预测结果(level 2)

图6 在5个不同速度下的各个关节自由度预测力矩的NRMSE(level 2)

图7 两个不同泛化水平的NRMSE

图8 两个不同泛化水平的相关系数

在一些文献[15-17]中除肌电信号外还分别使用肩和肘关节角、肩和肘关节角速度、肘关节角度和角速度作为输入预测肘关节力矩。为了进一步验证本文分析输入变量(EAV)预测结果的精确性，该输入变量预测的结果还与其他输入变量的预测结果在level 2的泛化水平下进行了比较。这些输入变量分别是：(1)相关肌肉的肌电信号，以及这些肌肉所驱动的关节自由度的角度(EA)；(2)相关肌肉的肌电信号，以及这些肌肉所驱动的关节自由度的角速度(EV)；(3)相关肌肉的肌电信号，以及力矩预测关节自由度的角度和角速度(EOAV)；(4)相关肌肉的肌电信号，以及力矩预测关节自由度的角度(EOA)；(5)相关肌肉的肌电信号(E)。其中相关肌肉是指关节自由度所关联的肌肉。通过实验比较发现以EA(NRMSE <14.25%)、EV(NRMSE<16.10%)、EOAV(NRMSE <14.94%)、EOA(NRMSE <15.53%)、E(NRMSE <19.51%)作为输入变量在不同速度下各关节自由度预测力矩的NRMSE与EAV(NRMSE <12.93%)的预测结果相比分别增加了10.2%、24.5%、15.5%、20.5%、50.9%，其比较结果如图9所示。其各关节自由度预测力矩的ρ比较结果如图10所示，由图可知本文分析的变量作为输入各关节自由度预测力矩的ρ明显优于其他变量输入的预测结果。

图9 不同输入变量的NRMSE比较

图10 不同输入变量的相关系数ρ值比较

3 讨论

实验结果表明以本文分析变量作为输入的ELM算法可以较精确地实现人体下肢各关节各自由度力矩实时预测。本方法与以往人体关节力矩人工智能预测研究[16,18-19,21,35-36]的不同主要体现在以下两点：(1)建立了人体关节力矩预测的数学模型，并获得其可在线测量的输入变量；(2)把极限学习机算法应用于关节力矩预测，加快了训练和预测的速度。

在整个实验过程中发现右踝关节内外翻展的力矩预测结果的归一化均方根误差和相关系数与其他关节自由度的力矩预测结果存在明显差异。通过分析右踝关节内外翻的反向生物力学计算的力矩，发现其符合N(0.0581，2.1669)的正态分布，与左踝关节内外翻力矩值相比可以认为其是噪声，所以不把踝关节内外翻力矩预测值作为参考。

由图7可知level 2和level 1相比NRMSE只降低了9.4%，所以该方法具有比较好的泛化能力。在实际应用中可以适当减少训练数据，以减少数据收集的时间。

由图9、图10可知EA与EAV相比其各个关节自由度预测力矩的最大平均归一化均方根误差只增长10.04%，平均相关系数略有降低(4.4%)，由此可知本文分析获得的输入变量中角速度和关节力矩的相关度并不高。在预测结果要求不严格的情况下可以不考虑把角速度作为输入，以减少传感器的数量降低实验的复杂度。另外，E与EAV其各个关节自由度预测力矩的NRMSE提高50.9%，ρ降低了20%，关节自由度力矩预测结果的准确度降低，这足以说明肌电信号的大小不能代表关节力矩值[37]且关节自由度的角度对关节力矩预测值影响非常大。并用EOAV是除EA以外预测精度最高的。由此可知被预测关节自由度的角度以及角速度对预测精度的影响比较大。

与反向生物力学分析相比，本文提出的方法首先只需要采集肌电信号和相关关节自由度的角度和角速度，不需要价格昂贵且场地要求比较高的3D步态捕捉系统[10-11]，使在普通医院和实验室进行人体关节自由度力矩的预测成为可能；其次本文提出的方法可以在训练的过程中不断适应个体差异，不需针对不同人进行肌肉骨骼模型的匹配，从而降低了个体差异带来的误差；最后本文提出的方法时间复杂度低，可以实现在线实时预测，使关节预测在实时性比较高的康复训练和机器人控制等方面的使用成为可能。

本研究还存在如下缺陷：(1)人工神经网络的输入和输出之间是黑盒子，无法表达输入和输出之间的关系。(2)考虑到多通道获取下肢肌电信号的复杂性比较大，本文只把一个下肢偏瘫患者的不同速度多步态周期数据作为实验数据。如果有足够多的数据其说服力会更强。(3)本文预测的准确性是和反向生物力学计算的力矩进行比较，虽然反向生物力学计算的关节力矩比较精准，但仍然存在误差，所以还无法完全表达本文方法预测关节力矩的精确度。