据统计，现代人平均每天约行走一万步，一生中约行走上亿步^[1]。人体双足行走时包含非常复杂的肌肉-骨骼-神经系统运动^[2-3]，其外在运动形式体现为自然的摆臂动作和类倒立摆^[4]的腿部动作，2种动作间的相互协调使行走成为一种比较省能的运动形式。这个特点使得行走成为人类日常运动的主要形式，研究其助力方式对提高人体活动能力具有重要意义。

髋关节和踝关节是为行走提供动力的主要部位。髋关节周围肌肉厚壮粗大，分布着众多双关节肌群，如股二头肌长头、股直肌和髂腰肌等，可进行跨关节的能量转移，对髋、膝、踝关节均能提供动力。从做功和能量代谢的角度来看，行走过程中髋关节所做的功为42%，大于踝关节的41%；髋关节的代谢能为45%，大于踝关节的32%^[5-6]。由此可见，髋关节肌群的作用要大于踝关节肌群。

鉴于髋关节在行走过程中的重要作用，许多研究团队对髋关节助力外骨骼进行了研究。德兰黑大学的Barazesh等^[7]设计的髋关节被动外骨骼通过利用与腘绳肌平行的弹性单元储存和释放能量的方法来对伸髋肌群助力，仿真结果表明助力后代谢能下降了10.0%，经实验验证穿戴外骨骼比不穿戴时代谢能下降8.9%。苏黎世联邦理工学院的Haufe等^[8]设计了一款髋关节被动助力外骨骼，其弹簧装置在步态周期内储存和释放能量，从而降低代谢能，实验结果表明：髋关节代谢功率没有增加，踝关节的腓肠肌与比目鱼肌的激活程度降低了28%。澳大利亚亚科廷大学的Shen等^[9]通过MATLAB软件和OpenSim软件的联合仿真，在髋、膝、踝关节处分别添加刚度为4.15，7.35，8.20 kN/m的弹簧，使得代谢能分别降低37.63%，4.85%，6.38%。上海理工大学康复研究院的胡冰山等^[10]设计了一款基于凸轮机构的髋关节变刚度储能外骨骼，并利用OpenSim软件分析了外骨骼刚度对代谢能的影响，结果表明，当刚度为171.9 Nm/rad时，代谢能降低效果最佳，降低了9%。

通过分析发现，多数学者在研究下肢外骨骼时仅考虑了对关节力矩的分析，而忽略了对肌肉-肌腱运动代谢能的分析。运动代谢能是评价下肢外骨骼助力效能的重要指标，若在设计阶段就考虑代谢能，则对助力设计的益处极大，但这很难以常用运动学、动力学分析方法来实现。首先，人体运动机械能与代谢能之间的映射关系较为复杂；其次，跨关节肌肉和弹性肌腱的存在，使得能量可在关节之间转移^[11]。因此，通过逆动力学求解得到的关节机械功可能会出现关节处输出的功与关节处肌肉所做的功不相等的情况。

人体代谢能消耗模型须根据肌肉-肌腱-骨骼特性建立，针对该问题，笔者提出基于Anybody软件与Hill肌肉模型建立髋关节肌肉-肌腱-骨骼简化模型，并对人体行走过程中髋关节的能耗进行分析，以此确定较好的助力方法和助力外骨骼机构设计方案。同时，通过建立人体-外骨骼一体化模型，对被动助力外骨骼机构的设计参数进行优化，并通过行走实验来验证所构建简化模型的正确性、外骨骼机构的助力效果以及研究方法的可行性，旨在为运动助力外骨骼的开发提供一种基于人体动力学分析和代谢能优化的研究方法。

根据人体解剖学知识，与矢状面内行走相关的下肢骨骼肌肉可分为伸展肌群和屈曲肌群。髋关节的伸展肌群主要由腘绳肌（包括半腱肌、半膜肌和股二头肌，以股二头肌为主）和臀大肌组成，其主要为膝关节屈曲和髋关节伸展提供动力，并保持膝关节的稳定性，是防止胫骨过度屈曲和膝关节过度伸展的重要动态稳定结构；屈曲肌群主要由髂腰肌和股前肌群（主要包括股直肌和缝匠肌，称为股匠肌群）组成，用于步态周期中支撑相到摆动相的过渡区间的屈髋抬腿。

以Anybody软件对1.5 m/s正常步态下的行走实验数据进行运动学、动力学仿真模拟。人体行走时单个步态周期的划分如图1所示，单个步态周期内髋关节的动力学仿真结果如图2所示。与以往行走步态研究文章中的结果相比，肌肉长度、肌肉激活程度、关节角度和关节力矩等数据的变化基本相符^[12-13]。

髋关节肌肉是支持行走运动的主要肌群。据研究，在行走过程中，踝关节力矩最大，但其部分能量由髋关节的跨关节肌群提供。同时，在以往研究中发现，仅对髋关节助力时，膝关节和踝关节的肌肉激活程度明显下降^[14]；与踝关节肌肉相比，髋关节肌肉-肌腱的刚性明显较大，即顺应度较低^[15]。因此，对踝关节的助力大部分仅作用于腿部而不是被内部肌腱吸收，故对髋关节助力效果更佳。

由图2可知，髋关节的伸展、屈曲肌群在单个完整的步态周期内交替发力。以股匠肌、髂腰肌为主的屈髋肌群的发力主要集中在步态周期的摆动相阶段，屈曲阶段髋关节角度为负。以腘绳肌、臀大肌为主的伸髋肌群的发力阶段分为两部分：在支撑相初期做正功加速伸髋，在摆动相做负功减速屈髋，从而完成支撑相到摆动相的抬腿迈步，伸展阶段髋关节角度为正。

助力决策的关键在于助力肌群的选择。分析图2可知：1）总体来看，伸髋肌群的总肌力与屈髋肌群相近，但伸髋肌群的功率峰值较大；2）从肌力来看，伸髋肌群的肌力曲线与关节角度具有更高的相关度，在整个步态周期中呈接近对称的形式，有利于被动助力外骨骼的布置；3）人体行走时的动力来源于足底蹬地获得的地面反力，而蹬地动作是通过伸展肌群发力使关节伸展来实现的。综上，若对髋关节伸展肌群助力，则可以加强前行动力或减小动力肌肉的负担，从而提高行走效率。因此，确定以伸髋肌群为助力对象，通过在髋关节处对目标肌群施加外辅助力矩来实现助力，根据髋关节代谢能最低原理确定助力强度。

图3所示为经典的Hill肌肉模型^[16]，该模型将肌肉简化为3组弹性元素相结合的收缩单元，包括串联主动收缩单元（contractile element, CE）、非线性串联弹性单元（series elastic element, SEE）和并联弹性单元（parallel elastic element, PEE），分别对应肌肉、肌腱与肌肉周围的结缔组织。虽然Hill肌肉模型没有直接模拟肌肉收缩的微观机制，但可以很好地再现肌肉的许多特性，被广泛应用于人体下肢运动模型的构建。图3中：Lmt为肌肉-肌腱单元总长度，Lt为肌腱单元长度，Lm为肌肉单元长度；α为羽状角，其大小决定了肌肉经肌腱作用于骨骼上力的大小；Ft为肌腱力。

考虑到个体差异，对Hill肌肉模型的参数进行归一化处理：L˜m=Lm/Lm0、v˜m=vm/vmax、L˜t=(Lt-Lt0)/Lt0。其中：L˜m、v˜m和L˜t分别为归一化处理后肌肉单元的长度、收缩速度和肌腱单元的长度，L_m0为肌肉单元在达到最大等距力Fmax时的长度，即最佳肌肉长度，v_m、vmax为肌肉的收缩速度和最大收缩速度，L_t0为肌腱松弛长度。

由图3可知，肌肉-肌腱单元总长度L_mt等于肌肉单元长度L_m和肌腱单元长度L_t之和，即：

三元素之间力的关系为：

式中：FMTU为总肌力；FM为肌肉单元及其周围结缔组织的收缩力。

肌肉单元及其周围结缔组织的力—长度关系分别如下：

式中：f˜mv˜m、f˜mL˜m分别为肌肉收缩速度影响因子和肌肉长度影响因子的归一化值^[5,17]；a为肌肉的激活程度。

结合上述公式，利用式(5)可求得肌肉运动的代谢功率P_met：

式中：fmetv˜m为代谢能影响因子。

利用Anybody软件仿真得到髋关节主要作用肌肉的基本参数，如表1所示。表中：PCSA（physiological cross-section area）表示肌肉生理横截面积。

为简化模型，在髋关节肌肉-肌腱-骨骼建模时作如下约定：

1）仅研究人体下肢行走过程中在矢状面内的运动。根据上文仿真结果，将髋关节运动视作由几部分主要肌肉（伸展肌群：臀大肌、股二头肌长头、股二头肌短头、半腱肌、半膜肌；屈曲肌群：股直肌、缝匠肌、髂腰肌）共同作用的结果^[18]，并将其简化为4个肌群：臀大肌群、腘绳肌群、髂腰肌群和股匠肌群。

2）基于Hill肌肉模型建立肌肉-肌腱-骨骼简化模型。

3）以肌群中最具代表性的肌束等效代替肌群的整体作用。每个肌群由多根肌束共同组成，而各肌束的Hill肌肉模型参数并不统一。为简化计算，选择各肌群中的一根肌束，将其作用力放大后等效代替整个肌群的作用力。

为获得下肢行走过程中髋关节各肌群的肌力、激活程度以及肌肉的基本参数，将利用AMTI生物力学测力板与Vicon运动捕捉系统获得的行走实验数据导入Anybody软件中的下肢模型，作为关节驱动参数，以对行走过程进行仿真分析。建立的髋关节肌肉-肌腱-骨骼简化模型如图4所示。将表1数据代入2.1小节中的数学模型，可分别计算得到各肌群的肌力和代谢能等，但各肌群的力臂不能直接计算得到，需要通过优化计算获得。

在Anybody软件的下肢模型中，以等效肌束替代肌群。例如：在计算腘绳肌群肌力时，选择肌力与肌群总肌力变化趋势和大小相近的股二头肌长头作为等效肌束，放大后的等效肌力按式(6)计算。设臀大肌群、腘绳肌群、髂腰肌群、股匠肌群的等效肌束的放大系数分别为y₁、y₂、y₃、y₄，利用MATLAB软件中的粒子群算法求解放大系数，优化目标函数设为各采样周期内肌群总肌力与等效肌力差值的均方根误差最小，其中各肌群总肌力与肌束力均由Anybody软件仿真获得。肌群肌力与等效肌力差值的均方根误差按式(7)计算。

式中：F_等效为等效肌力；F_肌束为肌束力；y为放大系数；R_EMS为均方根误差；F_肌群为仿真获得的肌力。

通过优化得到，臀大肌群、腘绳肌群、髂腰肌群、股匠肌群的等效肌束放大系数分别为6.56，1.75，5.24，1.40。图5所示为单个步态周期内髋关节的4个肌群的肌力与等效肌力的对比。为使所建立的肌肉-肌腱-骨骼简化模型能够复现人体下肢行走时髋关节的生物力学特性，建立一个约束优化函数来求解各肌群的作用力臂，以使模型输出的关节力矩与实验结果相匹配。

考虑到作用力臂会随运动发生变化，除肌肉起止点明显的股二头肌外，难以准确获取其他肌肉的作用力臂。因此，对于髂腰肌群、股匠肌群和臀大肌群，以定值力臂近似代表变值力臂；对于肌力作用线明显的腘绳肌群，根据几何关系计算不同时刻的力臂。根据实验中的标记点与骨骼之间几何关系以及刚体结构，计算得到股二头肌长头的实时力臂以及髋关节回转中心位置。设髂腰肌群、股匠肌群、臀大肌群的力臂长度分别为d₁、d₂、d₃，以关节力矩的均方根误差最小为优化目标，利用粒子群算法求解未知力臂，在计算髋关节的实时力矩时，规定顺时针为负，逆时针为正。根据文献[3]，设置各力臂的约束区间分别为：0.02<d₁<0.05 m、0.05<d₂<0.065 m、0.065<d₃<0.09 m。经优化得到各力臂的最优值分别为0.024 9，0.062 4，0.083 6 m。为验证优化后关节力矩的准确性，对基于优化力臂的髋关节力矩进行标准化计算（力矩/体重），并与实验结果进行比较，计算得到两者的均方根误差R_MSE=0.103 8 Nm/kg、相关度r²=0.95。根据文献[14,18]中关于人体简化模型准确度的判断条件（R_MSE<0.2 Nm/kg、r²>0.8），认为所构建的简化模型可以反映人体生物力学的真实情况。图6所示为利用优化力臂计算得到的髋关节力矩与仿真髋关节力矩的对比。由图可知，两者符合程度较高^[19]。

人体与外骨骼助力机构的结合并非简单的机构叠加，而是形成了一个相互影响的耦合系统，共同影响行走运动。因此，在对外骨骼助力机构进行效果评估和优化时，应该将两者作为一个整体来分析。为了在肌肉-肌腱-骨骼层面分析助力机构的作用，以便找到最佳的助力区间和助力曲线、优化机构结构参数以及被动助力所选用的弹性单元的刚度等，在肌肉-肌腱-骨骼简化模型中添加虚拟扭簧，以等效髋关节被动助力外骨骼的作用，如图7(a)所示。图7(b)所示为被动助力外骨骼机构原理。

根据图7，求得在该外骨骼助力下的髋关节各肌群的发力情况，再以Hill肌肉模型进行逆向求解，可分析肌肉-肌腱单元的肌力和代谢能，计算框架如图8所示。

在外辅助力矩作用下，求解得到髋关节伸展、屈曲肌群的残余作用力矩，结合上文的优化力臂，得到各肌束作用在髋关节上的力。由Hill肌肉模型可知，肌腱受力等于各肌束在羽状角影响下作用于关节上的力，利用肌腱力-长度计算公式可求得步态周期内不同时刻的肌腱长度。而整个肌束的长度只与运动状态下髋关节的角度有关，在已知肌腱长度与关节角度的情况下，可求得肌束长度，对其微分即可得到肌肉收缩速度。而代谢能大小只与串联主动收缩单元（CE）的参数有关，可先计算并联弹性单元（PEE）的收缩力，再由三元素之间的力学关系求得串联主动收缩单元的实时受力，结合计算得到的肌肉收缩速度等参数，求得代谢能和代谢功率。

助力外骨骼的设计目标为使人体关节总体能量消耗最少，因此以髋关节肌群的整体代谢能最低作为目标函数来优化所设计髋关节助力外骨骼的刚度。规定外骨骼大腿杆与人体躯干共线时的髋关节角度为0º，大腿杆向前为负，向后为正。图7(a)所示的模拟外骨骼的虚拟扭簧的作用如下：当髋关节角度小于0º时，大腿向前运动时扭簧在屈髋肌群牵引作用下伸长，做负功以储存能量，大腿向后运动时扭簧收缩，辅助伸髋肌群做正功以释放能量；当髋关节角度大于0º时，扭簧松弛，不做功。髋关节外骨骼助力过程如图9所示。

设虚拟扭簧的刚度为k，髋关节角度为β，则扭簧所提供的力矩M可表示为：

在理想情况下，应使外骨骼作用于髋关节的辅助力矩与目标助力肌群所提供的关节力矩相等，即肌肉不发力，仅通过外骨骼助力。因选用被动式助力，助力目标的拮抗肌群是助力机构的能量来源，且扭簧储存、释放能量随角度变化呈对称性，故无法实现理想助力矩曲线。因此，以扭簧辅助力矩与无助力时伸髋力矩差值的均方根误差最小为目标函数，利用优化算法求解最佳扭簧刚度k，这样可在满足助力最大的同时使拮抗肌群（即屈髋肌群）为扭簧储存能量而额外做功最小。最终的优化结果为k=39.500 1 Nm/rad，R_EMS=6.001 1 N·m，据此进行仿真，得到原始伸髋力矩与辅助力矩和有、无助力下的髋关节力矩对比曲线，结果如图10所示。从图10(a)中可以看出，外骨骼为伸髋动作提供了同步的辅助力矩，但其幅值基本不超过原始伸髋力矩，在助力阶段不会产生额外的拮抗肌力，其能量来源于步态周期中间阶段屈髋肌肉所做的功。由图10(b)可以看出，在外骨骼助力下，伸髋肌群提供的力矩下降，达到助力效果。

根据分配力矩计算得到肌力大小，由Hill肌肉模型逆向求解髋关节各功能肌群在单个步态周期内的代谢能和代谢功率。以代谢能最低为目标，优化分配各肌群占比系数，伸展肌群中腘绳肌群和臀大肌群的分配系数为0.24和0.76，屈曲肌群按发力区间计算。图11所示为基于Hill肌肉模型逆向计算得到的有、无助力下髋关节所有肌群的总代谢功率曲线。由图11可见，在步态周期的起止阶段，即伸髋阶段，代谢功率明显下降，对其积分可计算得到相较于无助力时髋关节肌群的代谢能下降了12.4%。

基于上文对外骨骼辅助力矩的优化以及代谢能的对比，设计可穿戴式髋关节助力外骨骼（弹性元件选用拉簧）。为验证模型及辅助力矩曲线的准确性，考虑到在行走过程中腰部在矢状面和冠状面内相对平稳这一特性，选择腰部为髋关节助力外骨骼的固定位。

髋关节屈伸角度、大腿和小腿的长度为已知的人体参数，行走过程中的髋关节角度可由动作捕捉系统获得。需要确定髋关节助力外骨骼机构下连接杆BC与大腿杆AB间夹角θ、下绑缚点到膝关节的距离L_BE 、上固定点D到髋关节转动中心的距离L_AD 和拉簧刚度K。设优化参数为 X =[x₁x₂x₃x₄]，其中x_i 分别表示上述4个未知参数，单位分别为rad、mm、mm、N/m。

考虑实际情况，设置上述4个参数的优化范围，使用MATLAB软件中的粒子群算法进行优化计算，目标函数设置为外骨骼机构产生的实际辅助力矩与图10(a)所示的仿真辅助力矩差值的均方根误差最小。经优化计算，得到 X =[1.7, 100, 200, 1 200]。根据该优化参数，计算外骨骼机构对髋关节的辅助力矩，并与图10(a)中辅助力矩仿真曲线进行比较，结果如图12所示。由图可见，两者在整个步态周期中比较接近。

为验证髋关节肌肉-肌腱-骨骼简化模型的正确性及所设计髋关节助力外骨骼的可行性，搭建图13所示的实验平台。选择身体健康的实验人员（男，26岁，身高1.76 m，体重67 kg）以1.5 m/s的速度在跑步机上行走。为验证所设计外骨骼的助力效果，选用肌电信号采集设备（Sichiray，双通道，1 000 Hz）来采集髋关节伸展、屈曲肌群的肌电信号，以评估肌肉的激活程度。因深层肌群的肌电信号无法直接测量，选择测量较为方便、位置明显的股直肌和股二头肌的肌电信号，这2块肌肉互为拮抗肌群，分别为髋关节伸展、屈曲提供动力。为减少穿戴外骨骼的不适应性，在正式实验前进行20 min的适应训练。

肌电信号与人体运动时的代谢能消耗呈正相关，采用肌电信号来间接表示代谢能的影响。实验者穿戴不同拉簧刚度（0.8，1.0，1.3 kN/m）的髋关节助力外骨骼，开展行走实验，并对行走过程中股直肌和股二头肌的肌电信号进行采集。选择在行走中间段的步态稳定区间内采集肌电信号，对采集到的肌电信号进行四阶巴特沃兹带通滤波器（20~400 Hz）滤波处理，随后选取合适的参考值进行归一化处理^[20-21]，结果如图14所示。

由图14可知，随着拉簧刚度的增大，股二头肌的平均肌电信号减小，在拉簧刚度为1.0~1.3 kN/m时，肌电信号减小趋势减缓，这是因为在外骨骼助力下肌肉自身发力减小，肌肉收缩力仅由肌肉组织中的结缔组织代为补偿。由Hill肌肉模型可知，结缔组织的收缩力只与肌肉单元变化长度相关，不会影响肌肉的激活程度，故导致该肌肉总体收缩力减小到一定程度后变小趋势减缓；作为目标助力肌群的拮抗肌群，股直肌在步态周期的中间阶段起屈髋作用并为外骨骼提供能量。因此，若外骨骼拉簧刚度增大，则会导致股直肌额外发力，即致使平均肌电信号增大。

由上文理论推导可知，外骨骼助力机构弹性元件的刚度存在最佳值，其结果表现为该刚度下伸、屈髋肌群的肌电信号数值求和最小。不同拉簧刚度下髋关节伸屈肌群的总肌电信号的统计结果如图15所示。由图可以看出，在无助力（刚度为0 kN/m）至拉簧刚度为1.0 kN/m范围内，髋关节肌群的总肌电信号呈现下降趋势，而在拉簧刚度为1.3 kN/m时表现为上升且高于无助力时。由于拉簧刚度的选择较少，无法明确取何值时整体肌电信号达到最小值，但可以看出最佳刚度存在于1.0 kN/m附近，接近于基于简化模型优化得到的刚度。

本文利用人体行走实验的运动学和力学数据驱动Anybody仿真，对行走过程中髋关节主要肌群的动力学表现进行了研究，并以此为依据提出了髋关节被动助力外骨骼的优化设计方案。同时，利用Hill肌肉模型、髋关节肌肉-肌腱-骨骼简化模型和髋关节-外骨骼耦合模型优化得到了外骨骼助力机构的刚度和结构参数，并分析了外骨骼在行走过程中对髋关节肌肉在力学性能、代谢能等方面的影响。实验结果表明，外骨骼的助力可有效地降低髋关节肌群的代谢能，验证了所提出的基于人体动力学分析的外骨骼助力设计方法的合理性，可为其他助力外骨骼的设计提供参考。

未来需要进一步开展的工作包括：1）以轻量化为目标，对现有外骨骼进行优化；2）精进模型，分别对矢状面、冠状面、额状面内的肌群发力情况进行分析；3）以主、被动方式对屈髋、伸髋肌群助力，以达到屈髋主动助力的同时为伸髋被动助力储存能量的效果；4）以采集耗氧量的方式，计算外骨骼助力下新陈代谢的降低效果。