CN106503304A - 一种在有限元人体模型中实现肌肉主动力的模拟方法及控制模型 - Google Patents

Abstract

本发明公开了一种在有限元人体模型中实现肌肉主动力的模拟方法及控制模型，将人体中枢神经系统进行运动控制过程中感觉系统的动态特性用传递函数模拟，建立基于PID闭环控制的肌肉主动控制模型；在有限元模拟平台采用关键字进行控制程序的编写，实现控制模型与力学模型的有机结合，不仅可实现肌肉主动控制功能，还使控制模型具有随研究对象的变化而调整的适应性。在有限元模型中采用了将三维实体单元与一维梁单元以共节点的方式并联，实现了被动力模拟和主动力模拟；采用控制理论计算出与碰撞情况相对应的激活等级，同时，控制模型具有随研究对象的变化而调整的适应性，而不是现有技术中采用预先设定激活等级，更加真实地还原了实际碰撞情况。

Description

一种在有限元人体模型中实现肌肉主动力的模拟方法及控制 模型

技术领域

本发明涉及人-机-环境的仿生技术领域，特别是涉及一种在有限元人体模型中实现肌肉主动力的模拟方法及控制模型。

背景技术

在汽车碰撞安全研究领域，人体有限元模型为损伤机理和损伤防护研究提供了有效的手段，许多研究工作者建立了基于人体解剖学结构的人体三维有限元模型，研究人体颈部在汽车碰撞事故中的生物力学特性、损伤机制及损伤准则，但这些研究通常只考虑了肌肉的被动效应。肌肉主动力的作用的缺失，将制约着损伤机理被完全重现和认知；此外，随着主动安全技术越来越得到重视和发展，主、被动安全技术的结合发展趋势，碰撞发生前的载荷施加过程及载荷等级所引起的肌肉主动反应对驾乘人员的运动学响应有着巨大的影响，这就要求有一种工具即可评估汽车安全系统在碰撞阶段的性能，又可评估其在碰撞前阶段的性能，即在人体有限元模型中实现肌肉的主动控制模拟。

目前，肌肉主动力建模方法主要是通过定义激活方程式等方法设定肌肉的激活等级来模拟神经系统对肌肉的刺激，并在仿真过程中定义某一个特定时刻作为激活初始时刻。这种方法定义的肌肉模型，其反射响应取决于所定义的时间常数，或预先定义的激活等级曲线。通过这种方法建立的肌肉主动力模型存在以下缺点：一是模型无法重现人体姿态的维持特性，如果预碰撞时间较长，有限元模型将可能倾倒、塌陷，而无法维持现实中乘员身体的正常姿态；二是无法模拟肌肉的条件反射响应，即不能产生与实际碰撞事故发生时的汽车加速度、安全带预紧等因素相匹配的肌肉响应。

发明内容

为了克服上述现有技术的不足，本发明提供了一种在有限元人体模型中实现肌肉主动力的模拟方法及控制模型，该方法将人体中枢神经系统进行运动控制过程中感觉系统的动态特性用传递函数模拟，在有限元软件中采用关键字编写控制程序，建立基于PID闭环控制的肌肉主动控制模型，将控制模型与肌肉力学模型有机结合，在有限元模拟平台中实现人体肌肉主动控制功能。

一种在有限元人体模型中实现肌肉主动力的模拟方法，通过有限元实体肌肉单元的材料特性定义肌肉被动力特性，采用控制模型与力学模型相耦合的方式，根据Hill肌肉本构模型，对肌肉采用三维实体单元与一维梁单元以共节点的方式并联，定义肌肉的主动力响应特性；根据肌肉的收缩速度和收缩量采用PID反馈控制，生成模拟肌肉刺激信号，利用肌肉刺激信号进行肌肉激活等级的调整，利用得到的肌肉激活等级产生对应的肌肉主动力，将肌肉主动力作用于人体有限元模型，实现肌肉主动力模拟。

所述肌肉的主动力F_CE根据Hill肌肉材料模型，按照公式获得：F_CE＝σ_ce*A_pcsa，且采用*MAT_MUSCLE定义肌肉的主动力属性；

其中，σ_ce＝A(t)f_l(l)f_v(v)σ_max，A(t)为肌肉激活等级，依据肌肉受到碰撞产生的刺激信号取值；f_l(l)和f_v(v)分别为为肌肉力－肌肉长度型函数以及肌肉力－收缩速度型函数，σ_max为最大收缩应力，A_pcsa为肌肉的生理横截面积；σ_max、A_pcsa、f_l(l)以及f_v(v)均为肌肉的生物力学参数，根据人体生物力学计算获得；

肌肉激活等级依据肌肉受到碰撞产生的刺激信号过程采用两个串联的低通滤波器模拟，第一个低通滤波器模拟由肌肉刺激信号u(t)获得神经刺激信号Ne(t)的过程，由神经刺激时间常数T_ne的大小调节神经刺激信号的输出速率；第二个滤波器模拟肌肉激活/钝化的过程，通过调整激活时间常数T_naa和钝化时间常数T_nad的大小获得肌肉激活等级A(t)：

其中，T_na表示激活时间常数T_naa或钝化时间常数T_nad，当模拟过程为激活时，T_na表示激活时间常数T_naa；当模拟过程为钝化时，T_na表示钝化时间常数T_nad。

所述模拟肌肉刺激信号按照公式计算获得：

且通过物理建模的方法，利用LS-DYNA求解器实现积分和微分运算，获取公式中误差信号的积分和微分曲线；

其中，e(t)＝r(t)-y_d(t)，e(t)表示肌肉内骨骼关节弯曲偏差信号，r(t)表示肌肉内骨骼关节碰撞前弯曲程度参考信号，y_d(t)表示肌肉内骨骼关节碰撞后的弯曲程度经过神经系统延迟后的反馈信号，k_p、k_i以及k_d分别表示PID控制算法中的比例系数、积分系数以及微分系数；表示对肌肉内骨骼关节弯曲偏差信号e(t)进行积分运算；de(t)/dt表示对肌肉内骨骼关节弯曲偏差信号e(t)进行微分运算。

PID控制是工程中非常成熟的控制方法，三个参数的整定是根据稳定性、快速性、准确性的原则，以消除偏差为目的，按照先比例后积分再微分的步骤，逐步凑试获得；

所述误差信号的积分和微分曲线通过*BOUNDARY_PRESCRIBED_MOTION_NODE定义某骨骼关节节点的速度曲线为偏差信号e(t)，并通过*DEFINE_CURVE_FUNCTION的指令dx与accx获取该骨骼关节节点的位移和加速度曲线，实现积分ie(t)和微分运算de(t)。

采用两个串联的滤波器对PID控制模块生成模拟肌肉刺激信号进行整定，使得利用刺激信号u(t)调整后的肌肉激活等级在最小极限0.5％和最大极限70％之间。

根据Janet Brelin-Fornari的研究可知，肌肉激活等级的最大极限为70％，最小极限为0.5％^[17]，激活等级越大，肌肉的主动收缩力越大。激活等级整定模块对刺激信号u(t)进行规范并将其限制在最大最小极限之间，得到可用的激活等级A(t)。

所述神经系统延迟利用一阶惯性环节的相位滞后特性模拟，具体过程如下：

在有限元模型中定义两个参考点ni和nj，用纯粘性属性的离散单元连接这两个参考点，C为离散单元的粘性参数值，设置参考点nj质量为m，并使用关键字*BOUNDARY_PRESCRIBED_MOTION_NODE设置参考点ni的输入利用关键字*HISTORY_NODE获取参考点nj的输出则神经系统延迟的传递函数为：ω表示神经系统延迟的角频率；

采用Vlugt实验获得的神经系统延迟时间T_de，令一阶惯性系统的群时延等于T_de，按以下公式确认m和C的取值：

神经延迟取决于刺激信号往返中枢神经所需要的时间，根据Vlugt实验^[15]获得的神经延迟时间Tde，令一阶惯性系统的时延等于Tde，根据一阶惯性系统的群时延表达式确定合适的m，C值，从而获得需要的延迟时间。

所述通过有限元实体肌肉单元的材料特性定义肌肉被动力特性是指三维实体单元采用超弹性Ogden橡胶材料模型，且肌肉设置为具有高度非线性粘弹性的被动力属性。

有益效果

本发明提供了一种在有限元人体模型中实现肌肉主动力的模拟方法及控制模型，通过对人体骨骼肌系统复杂的动力学特性研究，以及对肌肉主动力产生的影响因素挖掘的基础上，将人体中枢神经系统进行运动控制过程中感觉系统的动态特性用传递函数模拟，建立基于PID闭环控制的肌肉主动控制模型；在有限元模拟平台采用关键字进行控制程序的编写，实现控制模型与力学模型的有机结合，不仅可实现肌肉主动控制功能，还使控制模型具有随研究对象的变化而调整的适应性。在本方法中尤其在有限元模型中采用了将三维实体单元与一维梁单元以共节点的方式并联，实现了被动力模拟和主动力模拟；而且采用控制理论计算出与碰撞情况相对应的激活等级，同时，控制模型具有随研究对象的变化而调整的适应性，而不是现有技术中采用预先设定激活等级，更加真实地模拟了实际碰撞情况，为安全产品的研发提供了理论基础。

附图说明

图1为基于PID的颈部肌肉反馈控制模型；

图2为颈部关节角计算参考点设置；

图3一阶惯性环节的物理建模示意图。

具体实施方式

下面将结合附图和实施例对本发明做进一步的说明。

一种在有限元人体模型中实现肌肉主动力的模拟方法，包括以下几个步骤：

步骤(1)对人体有限元模型的肌肉模型进行重建，定义肌肉的被动属性和主动属性。

步骤(2)、通过MRI影像资料获得人体骨骼肌三维模型的几何，根据人体解剖学结构中肌肉的起止端，将肌肉与骨骼模型进行匹配。

步骤(3)、根据Hill肌肉本构模型，采用三维实体单元与一维梁单元以共节点的方式并联，完成肌肉模型重建，模拟骨骼肌的主被动力响应特性。

其中三维实体单元采用超弹性Ogden橡胶材料模型(*MAT_OGDEN_RUBBER)，定义肌肉具有高度非线性粘弹性的被动力属性；一维梁单元采用Hill肌肉材料模型(*MAT_MUSCLE)，定义肌肉的主动力属性。

步骤(4)、根据Hill肌肉材料模型，每条肌肉的主动力F_CE均可按以下公式计算得到：

σ_ce＝A(t)f_l(l)f_v(v)σ_max (1)

F_CE＝σ_ce*A_pcsa (2)

式中A(t)为肌肉激活等级，f_l(l)和f_v(v)分别为为肌肉力－肌肉长度型函数以及肌肉力－收缩速度型函数，σ_max为最大收缩应力，A_pcsa为肌肉的生理横截面积。

其中σ_max A_pcsa f_l(l)f_v(v)为肌肉的生物力学参数，已有大量研究结论，在材料模型*MAT_MUSCLE中定义即可，肌肉激活等级A(t)与碰撞事故发生时的状况有关，即与碰撞事故发生时的肌肉收缩量和收缩速度有关；

肌肉激活等级依据肌肉受到碰撞产生的刺激信号过程采用两个串联的低通滤波器模拟，第一个低通滤波器模拟由肌肉刺激信号u(t)获得神经刺激信号Ne(t)的过程，由神经刺激时间常数T_ne的大小调节神经刺激信号的输出速率；第二个滤波器模拟肌肉激活/钝化的过程，通过调整激活时间常数T_naa和钝化时间常数T_nad的大小获得肌肉激活等级A(t)：

其中，T_na表示激活时间常数T_naa或钝化时间常数T_nad，当模拟过程为激活时，T_na表示激活时间常数T_naa；当模拟过程为钝化时，T_na表示钝化时间常数T_nad。

以颈部为例，颈部肌肉的生物力学参数参考已发表研究结果，如表1所示，

表1.颈部肌肉参数

以颈部为例，颈部肌肉的控制参数如表2所示；

表2肌肉主动力控制系统参数

步骤(5)、采用PID反馈闭环控制技术根据碰撞时刻的实际情况调整肌肉激活等级A(t)。肌肉的主动作用反映在预碰撞中的姿势维持以及紧急情况下的条件反射；

采用keyword关键字编程建立控制模型，根据碰撞情况调整肌肉主动属性的激活等级Na(t)，产生与碰撞强度相适应的肌肉主动力F_ce，再作用于人体有限元模型，实现模型在预碰撞中姿势维持以及乘员在紧急情况下的条件反射模拟。

本发明正是依据这一思路来进行肌肉激活等级A(t)的调整。人体姿势的保持及改变主要通过关节角的变化来实现，以头颈部模型为例，如图1所示，颈部实时弯曲角度为输出y(t)，并作为反馈信号y_d(t)参与颈部运动的调节；g(t)为头颈部有限元模型碰撞仿真的运动输入，通常为碰撞加速度信号；参考信号r(t)为参考关节角，即碰撞前颈部的弯曲程度；反馈信号y_d(t)经神经延迟作用后与参考信号r(t)进行比较，得到偏差信号e(t)；PID控制模块根据偏差信号的大小作出相应调控，得到肌肉刺激信号u(t),模拟中枢神经对肌肉收缩进行调控的电位信号；激活等级生成模块根据肌肉激活理论对刺激信号u(t)进行调整和规范，从而获得肌肉激活等级A(t)，控制肌肉的激活程度；A(t)作用于头颈部有限元模型的Hill肌肉单元，控制头颈部关节弯曲角度的变化，实时角度的输出反馈又将产生新的误差，这样反复调整直到误差为零，从而实现了肌肉主动力对头颈部姿势的调节。

步骤(6)、在有限元平台上实现步骤(5)所述的肌肉主动力模拟。

在实现上述方法中，所述的肌肉内骨骼关节角度信号通过矢量计算完成测量和计算。

矢量计算是成熟的方法，以头颈部为例来阐述通过矢量计算获得颈部关节角。

以头颈部模型为例，如图2所示分别取人体模型中的第一胸椎T1中心点和第一颈椎C1中心点为基准点n1、n2，根据式3～5可计算两点所成的矢量与竖直方向的夹角，从而获得颈部关节角。

U_x＝n_1x-n_2x (3)

Uz＝n_1z-n_2z (4)

n_1x和n_1z分别为某时刻T1基准点n1的x、z坐标值，n_2x和n_2z为C1基准点n2的x、z坐标值，U_x和U_z分别为x向和z向矢量，y(t)则为该时刻的颈部关节角。

通过LS-DYNA函数库中的指令CX(节点ID)、CY(节点ID)、CZ(节点ID)获取基准点n1、n2的实时坐标，并通过*DEFINE_CURVE_FUNCTION完成公式计算，实时颈部弯曲角的计算。

本方法中神经延迟的模拟利用简单一阶惯性环节的相位滞后特性来实现神经延迟的模拟，在有限元环境中构建一阶惯性环节的思路如图3所示，在有限元模型中定义两个参考点(ni,nj)，用纯粘性属性的离散单元连接这两个参考点，C为离散单元的粘性参数值，定义nj质量为m，并使用关键字*BOUNDARY_PRESCRIBED_MOTION_NODE定义点ni的输入利用关键字*HISTORY_NODE获取nj的输出则其传递函数为：

神经延迟取决于刺激信号往返中枢神经所需要的时间，本文采用Vlugt实验获得的神经延迟时间T_de，令一阶惯性系统的群时延：

选定合适的m，C值，即可获得需要的延迟时间T_de。

源代码如下：

基于PID控制生成肌肉刺激信号的具体过程如下：

PID控制器对颈部实时弯曲角度与碰撞前的弯曲角度之间的误差信号进行比例、积分、微分运算，得到姿势保持所需的相应肌肉刺激信号u(t)，计算过程如式7～8所示。LS-DYNA函数库中没有积分和微分函数，本发明通过物理建模的方法，利用LS-DYNA求解器实现积分和微分运算，获取误差信号的积分和微分曲线。具体过程如下：通过*BOUNDARY_PRESCRIBED_MOTION_NODE定义某节点的速度曲线为偏差信号e(t)，并通过*DEFINE_CURVE_FUNCTION的指令dx(节点id)与accx(节点id)获取该节点的位移和加速度曲线，即分别完成了积分和微分运算ie(t)、de(t)，然后利用LS-DYNA关键字编写PID计算程序，完成PID控制模块的建立。

e(t)＝r(t)-y_d(t) (7)

源代码如下：

对激活等级整定模块的具体实现过程如下：

采用两个串联的滤波器实现刺激信号的整定。PID控制模块生成模拟肌肉用以维持姿势所受的电位信号u(t)，但该信号不能直接作用于力学模型，肌肉激活等级的最大极限为70％，最小极限为0.5％，需采用激活等级整定模块对刺激信号u(t)进行规范将其限制在最大最小极限之间才能得到可用的激活等级A(t)。源代码如下：

已有的肌肉主动力建模通常通过预先设定固定的激活等级曲线来实现。通过这种方法建立的肌肉主动力模型无法模拟肌肉的条件发射响应，即不能产生与实际碰撞事故发生时的汽车加速度、安全带预紧等因素相匹配的肌肉响应。

本发明在Hill肌肉模型的基础上加入PID反馈闭环控制技术，将控制技术与肌肉力学模型有机结合，实现神经系统对肌肉主动作用的模拟，并可通过PID参数的变化使肌肉模型具有随研究对象的变化而调整的特性。

以上应用了具体个例对本发明进行阐述，只是为了帮助本领域中的普通技术人员很好的理解。在不偏离本发明的精神和范围的情况下，还可以对本发明的具体实施方式作各种推演、变形和替换。这些变更和替换都将落在本发明权利要求书所限定的范围内。

Claims (8)

1.一种在有限元人体模型中实现肌肉主动力的模拟方法，其特征在于，通过有限元实体肌肉单元的材料特性定义肌肉被动力特性，采用控制模型与力学模型相耦合的方式，根据Hill肌肉本构模型，对肌肉采用三维实体单元与一维梁单元以共节点的方式并联，定义肌肉的主动力响应特性，根据肌肉的收缩速度和收缩量采用PID反馈控制，生成模拟肌肉刺激信号，利用肌肉刺激信号进行肌肉激活等级的调整，利用得到的肌肉激活等级产生对应的肌肉主动力，将肌肉主动力作用于人体有限元模型，实现肌肉主动力模拟。

2.根据权利要求1所述的方法，其特征在于，所述肌肉的主动力F_CE根据Hill肌肉材料模型，按照公式获得：F_CE＝σ_ce*A_pcsa，且采用*MAT_MUSCLE定义肌肉的主动力属性；

其中，σ_ce＝A(t)f_l(l)f_v(v)σ_max，A(t)为肌肉激活等级，依据肌肉受到碰撞产生的刺激信号取值；f_l(l)和f_v(v)分别为为肌肉力－肌肉长度型函数以及肌肉力－收缩速度型函数，σ_max为最大收缩应力，A_pcsa为肌肉的生理横截面积；σ_max、A_pcsa、f_l(l)以及f_v(v)均为肌肉的生物力学参数，根据人体生物力学计算获得；

肌肉激活等级依据肌肉受到碰撞产生的刺激信号过程采用两个串联的低通滤波器模拟，第一个低通滤波器模拟由肌肉刺激信号u(t)获得神经刺激信号Ne(t)的过程，由神经刺激时间常数T_ne的大小调节神经刺激信号的输出速率；第二个滤波器模拟肌肉激活/钝化的过程，通过调整激活时间常数T_naa和钝化时间常数T_nad的大小获得肌肉激活等级A(t)：

$T_{ne}\frac{dNe\left(t\right)}{dt}=u\left(t\right)-Ne\left(t\right)$

$T_{na}\frac{dA\left(t\right)}{dt}=Ne\left(t\right)-A\left(t\right)$

其中，T_na表示激活时间常数T_naa或钝化时间常数T_nad，当模拟过程为激活时，T_na表示激活时间常数T_naa；当模拟过程为钝化时，T_na表示钝化时间常数T_nad。

3.根据权利要求2所述的方法，其特征在于，所述模拟肌肉刺激信号按照公式计算获得：且通过物理建模的方法，利用LS-DYNA求解器实现积分和微分运算，获取公式中误差信号的积分和微分曲线；

其中，e(t)＝r(t)-y_d(t)，e(t)表示肌肉内骨骼关节弯曲偏差信号，r(t)表示肌肉内骨骼关节碰撞前弯曲程度参考信号，y_d(t)表示肌肉内骨骼关节碰撞后的弯曲程度经过神经系统延迟后的反馈信号，k_p、k_i以及k_d分别表示PID控制算法中的比例系数、积分系数以及微分系数；表示对肌肉内骨骼关节弯曲偏差信号e(t)进行积分运算；de(t)/dt表示对肌肉内骨骼关节弯曲偏差信号e(t)进行微分运算。

4.根据权利要求3所述的方法，其特征在于，所述误差信号的积分和微分曲线通过*BOUNDARY_PRESCRIBED_MOTION_NODE定义某骨骼关节节点的速度曲线为偏差信号e(t)，并通过*DEFINE_CURVE_FUNCTION的指令dx与accx获取该骨骼关节节点的位移和加速度曲线，实现积分ie(t)和微分运算de(t)。

5.根据权利要求4所述的方法，其特征在于，采用两个串联的滤波器对PID控制模块生成模拟肌肉刺激信号进行整定，使得利用刺激信号u(t)调整后的肌肉激活等级在最小极限0.5％和最大极限70％之间。

6.根据权利要求3-5任一项所述的方法，其特征在于，所述神经系统延迟利用一阶惯性环节的相位滞后特性模拟，具体过程如下：

在有限元模型中定义两个参考点ni和nj，用纯粘性属性的离散单元连接这两个参考点，C为离散单元的粘性参数值，设置参考点nj质量为m，并使用关键字*BOUNDARY_PRESCRIBED_MOTION_NODE设置参考点ni的输入利用关键字*HISTORY_NODE获取参考点nj的输出则神经系统延迟的传递函数为：ω表示神经系统延迟的角频率；

采用Vlugt实验获得的神经系统延迟时间T_de，令一阶惯性系统的群时延等于T_de，按以下公式确认m和C的取值：

7.根据权利要求6所述的方法，其特征在于，所述通过有限元实体肌肉单元的材料特性定义肌肉被动力特性是指三维实体单元采用超弹性Ogden橡胶材料模型，且肌肉设置为具有高度非线性粘弹性的被动力属性。

8.一种在有限元人体模型中实现肌肉主动力的模拟控制模型，其特征在于，基于权利要求1-7任一项所述的方法，该控制模型包括PID控制模块、激活等级生成模块、Hill肌肉单元、肌肉骨骼模型以及神经系统延迟模块；

所述PID控制模块、激活等级生成模块、Hill肌肉单元以及肌肉骨骼模型依次相连，肌肉骨骼模型的输出端与神经系统延迟模块的输入端相连，将实时输出的肌肉内骨骼关节碰撞后的弯曲程度信号输入至神经系统延迟模块；神经系统延迟模块的输出端与肌肉内骨骼关节碰撞前参考信号端相连，两者相减后得到的差值信号端与PID控制模块的输入端相连。