融合整体法和Simulink的人体肌肉力矩动态仿真计算

王彦凤，王瑞，蔡玉强，赵欣

(1. 华北理工大学 机械工程学院，河北 唐山 063210；2.唐山市拓又达科技有限公司，河北 唐山 063000)

0 引言

穿戴式下肢外骨骼康复机器人是一种典型的外骨骼助力装置[1]，融合了传感、控制、信息获取、移动计算等机器人技术, 帮助操作者在无意识控制下完成站立、行走等功能[2]。人体肌肉力量作为康复机器人共享控制和自适应控制中常用的一种参量，很多学者对肌肉力量的计算进行了研究。娄佐政[3]根据人体解剖学建立了人体体段肌肉骨骼模型，然后据此模型在 MATLAB 开发平台上设计了一套可视化的人体关节静态力学评估软件 Bio SAFE估计人体的关节力矩和肌肉力。李海龙等[4-6]通过提取患者下肢肌电信号在时域内的特征量来测量下肢肌肉的肌肉力量。张希安等[7-9]使用骨肌仿真技术，选取 Hill 三元素肌肉模型与Lagrange 动力学模型相结合，建立骨肌系统模型计算下肢肌肉力矩。贾惠珍等[9]利用计算机和高速数字摄像技术，通过量化计算产生力的速率(RFD)，测量等长收缩所需的时间来测定人体肌肉力量。

目前，针对肌肉力量测量研究中主要存在三方面的问题:(1)只做静态动力学分析。由于人体步态规律变化复杂，单纯的静力学分析没有考虑肢体惯性力、离心力等动态因素；(2)建立数学模型时将下肢简化为多刚体的系统组合，没有反应出肌肉等柔性组织对下肢关节力矩计算的影响；(3)测量存在的辅助设备多、信号分析复杂、易受环境变化影响和反馈滞后。为了快速、准确地求出下肢肌肉力量的变化规律，该项研究提出了一种融合整体法和 MATLAB /Simulink 的动态仿真计算方法。整体法主要解决2个以上物体组成的系统力学问题[10]，使用整体法对人体-外骨骼系统应用牛顿欧拉法建立动力学方程，并加入了肢体的被动粘弹性扭矩，对模型进行了完整的主动步态训练仿真计算，求解出了髋关节与膝关节的肌肉力矩，验证了所提方法的正确性。

1 建立人体-外骨骼系统动力学方程

1.1 人体肌肉骨骼动态力学分析

人体下肢可看作刚性体和柔性体的组合。其中骨骼为刚性体，肌肉、关节软骨等为柔性体。肌肉骨骼的动作过程等效如图1所示。在使用穿戴式外骨骼机器人进行步态训练时，下肢的角度变化大多集中在矢状面内，并且踝关节多数为被动自由度，这为问题分析提供了便利。在研究中作以下假设：(1)只考虑人体矢状面的下肢情况；(2)股骨和胫骨为刚体，体段间由铰链连接，杆件间无间隙无摩擦；(3)将脚与小腿作为一个整体。人体下肢模型的运动方程表示为：

(1)

图1 肌肉骨骼的动作过程

1.2 外骨骼动态力学分析

外骨骼依据人体仿生学设计，与人体生理结构高度相似度，所以与人体有相同的运动学特性。外骨骼机器人产生的驱动力矩转化为以下两部分：(1)克服外骨骼自身运动中的惯性力、重力、离心力和科里奥利力；(2)补充给人体运动，施加给下肢的τext。由此得到外骨骼机器人的运动方程为

(2)

1.3 人体-外骨骼系统动态力学分析

整体法适用于求解系统内各物体保持相对静止或匀速直线运动的情况，也适用于相对作匀变速直线运动的情况[11]。外骨骼带动人体进行康复训练时二者相对静止，而且为了使患者更舒适、安全，外骨骼作匀变速直线运动，使用整体法分析人体-外骨骼系统更加便捷合理。人体-外骨骼系统模型如图2所示。

图2 人体-外骨骼系统模型

该模型由下肢外骨骼与人体下肢组成，各有5个参数定义:关节长度(L)、体段质量(m)、质心在连杆平行方向上的位置(a)、惯性矩(Iz)、关节角(θ)。系统的牛顿欧拉方程为：

(3)

HR指标代表了人类-外骨骼系统的组合，即MHR=MH+MR,VHR=VH+VR,GHR=GH+GR具体如式(4)-式(21)所示,式中角标1代表大腿，角标2代表小腿。

高技术产业是我国经济的战略性先导产业，对推动产业结构升级和经济增长方式转变具有重要的推动作用。高技术产业的发展水平关系到一个国家或地区的综合竞争力，近年来我国高技术产业得到了较快发展，为我国的经济发展做出了巨大贡献。但高技术产业快速发展的同时也存在许多问题，如高技术产业的区域发展不均衡，在技术水平、技术创新效率等方面与发达国家还存在较大差距。由于高技术产业具备技术和知识集中、资源消耗量少等特点，使得技术创新效率对产业发展起到更为重要的引导作用，也决定了各地区高技术产业的总体发展水平。因此，如何科学有效地对区域高技术产业创新效率进行评价，对于区域高技术产业的健康发展具有重要意义。

(4)

MHR,11=JHR1+2L1XHR2cosθ2

(5)

MHR,11=JHR1+2L1XHR2cosθ2

(6)

MHR,12=MHR,21=JHR2+L1XHR2cosθ2

(7)

MHR,22=JHR2

(8)

(9)

(10)

(11)

(12)

GHR,1=g(XHR1sinθ1+XHR2sinθ12)

(13)

GHR,2=gXHR2V(14)

(15)

(16)

(17)

(18)

(19)

(20)

X2=m2a2(21)

式(16)和式(17)使用Riener的双指数方程[12]对P(θ)进行计算。为了简洁起见，在式(18)～式(21)中省略了角标H、R和HR。人体体段惯性参数为了符合中国人的特点，参照我国的成年人体生物力学参数(GB/T17245-2004)标准[13]，应用多元线性回归方程进行计算。

2 建立Simulink仿真模型

Simulink提供了一个动态系统建模、仿真和综合分析的集成环境，是MATLAB最重要的组件之一[14]。人体-外骨骼的动力学方程结合Simulink，动态仿真计算流程如图3所示。未知量外骨骼机器人关节角度和关节驱动力矩由康复机器人自带的角度编码器和扭矩传感器进行计算，只需将动力学方程模块化嵌入到仿真模型中便可实时动态计算出患者的关节肌肉力矩。建立仿真计算平台如图4所示。

图3 动态仿真计算流程

图4 仿真计算平台

3 下肢主动康复训练实验及结果分析

3.1 肌肉扭矩测量实验过程

以某公司自主研发的自平衡下肢外骨骼康复机器人为实验平台，如图5所示，该款外骨骼康复机器人的髋关节和膝关节内置角度编码器和扭矩传感器。

图5 外骨骼康复机器人

为了验证所提方法的准确性，在每位受试者的右腿臀大肌、股内侧肌、股直肌、半腱肌、腓肠肌上贴有测量电极片，用以获得表面肌电信号，通过等距收缩实验建立表面肌电信号(EMG)数据和关节扭矩之间的关系，以此作为对比实验。选取3名成年男性和2名成年女性作为实验目标，其中1名受试者和与他相符合的外骨骼康复机器人基本参数如表1所示。

表1 人体与外骨骼基本参数

5名受试者因患病程度不同，因此肌肉有力量也不相同。当人产生轻微的迈步动作时就会被外骨骼的传感系统捕捉到，控制系统开始逐步的增加外骨骼的驱动扭矩，直至人体-外骨骼系统整体关节速度达到预先的设定的正常人步态速度；对于下肢完全没有力量的患者由医护人员辅助其迈步。实验过程如图6所示。受试者每次进行10个步幅循环，测量3次，每次间隔30 min,每次使用之前需要重新标定测量位置。角度编码器与外骨骼扭矩传感器的测量数据实时传输至上位机，由MATLAB计算肌肉扭矩。表面肌电信号测量设备将电极采集的原始肌电图信号经Bagnoli 8通道系统放大至增益1 000，然后进行带通滤波，所有信号的采样速率均为1 000 Hz。

图6 实验过程

3.2 实验结果分析

将5名受试者的外骨骼角度编码器和扭矩传感器获得的数据以数值矩阵格式导入MATLAB工作空间，取仿真时间为3 s、固定步长0.000 1。其中1名受试者的simulink仿真结果如图7所示；表面肌电信号测得的模拟量，经过等距肌肉收缩实验换算转换后的肌肉扭矩如图8所示。

图7 MATLAB/simulink计算肌肉扭矩

图8 表面肌电信号测量肌肉扭矩

从图7与图8曲线的变化规律来看，所提方法获得的下肢肌肉扭矩与表面肌电信号测量出的下肢肌肉扭矩十分相似。为了更加直观的比较两种方法,将五位受试者由simulink计算出的肌肉扭矩作为观测值，表面肌电信号测量的肌肉扭矩作为真值，求出两种方法获取的肌肉扭矩的均方根误差，其计算公式如下：

(22)

式中n为受试者人数5；Xs,i为simulink计算的第i位受试者在j时刻的肌肉扭矩；XsEMG,i

使用MATLAB对其进行计算，髋关节肌肉扭矩的RMSE最大值为2.36，膝关节肌肉扭矩的RMSE最大值为2.89。两种方法获取的髋关节肌肉扭矩和膝关节肌肉扭矩偏差很小，证明所提方法能够准确测量穿戴外骨骼康复机器人进行主动康复训练时的人体下肢肌肉扭矩。

4 结论

(1)提出了一种融合整体法和simulink动态仿真计算人体肌肉扭矩的方法。将进行主动康复训练的人体下肢与穿戴式外骨骼康复机器人看做一个整体，并将该整体简化为多刚体系统，用牛顿-欧拉法求解动力学问题，建立动力学方程。

(2)使用simulink搭建仿真计算平台，将角度编码器与扭矩传感器获得的数字量导入MATLAB，进行人体肌肉扭矩计算。通过与表面肌电信号实验结果进行对比，验证了所提方法的可行性。

(3)该方法不需要在受试者身上安装传感器，提高了测量的舒适性；信号的采集、分析和计算也比传统的方法更简便、省时。可作为穿戴式外骨骼康复机器人康复效果测评和主动训练时共享控制及自适应控制的一种方法。