基于肌肉骨骼模型的预测仿真

肌肉骨骼模型

骨骼的刚体动力学模型

首先基于刚体动力学对于人体骨骼系统进行建模。首先确定系统的广义坐标qqq，对于图中所示的二维模型为例，广义坐标包含躯干的水平、数值方向位移和姿态角，以及每条腿三个关节角，一共3+2∗3=93+2*3=93+2∗3=9个自由度。之后通过这些广义变量和广义变量的导数，表示出系统的动能和势能，并构造拉格朗日函数：
L=T−VL=T-VL=T−V
考虑骨骼系统仅收到地反力和肌肉力的作用，通过拉格朗日方程得到系统的动力学模型：
ddt(∂L∂q˙j)−∂L∂qj=τcontactj+τmusclej\frac{\mathrm{d}}{\mathrm{d} t}\left(\frac{\partial L}{\partial \dot{q}_{j}}\right)-\frac{\partial L}{\partial q_{j}}=\tau^{j}_{contact}+\tau^{j}_{muscle}dtd(∂q˙j∂L)−∂qj∂L=τcontactj+τmusclej

M(q)⋅q¨+B(q,q˙)=J(q)TFcontact+R(q)Fmuscle\mathbf{M}(\mathbf{q}) \cdot \ddot{\mathbf{q}}+\mathbf{B}(\mathbf{q}, \dot{\mathbf{q}})=\mathbf{J(q)}^T\mathbf{F}_{contact}+\mathbf{R(q)}\mathbf{F}_{muscle}M(q)⋅q¨+B(q,q˙)=J(q)TFcontact+R(q)Fmuscle

其中M(q)\mathbf{M(q)}M(q)为质量矩阵，B(q,q˙)\mathbf{B}(\mathbf{q}, \dot{\mathbf{q}})B(q,q˙)包含重力、离心力和科氏力，J(q)\mathbf{J(q)}J(q)表示雅克比矩阵，将地反力转化为对每个广义坐标的等效力矩，R(q)\mathbf{R(q)}R(q)表示肌肉力对广义坐标的等效力矩。下面我们主要讨论等式右侧的肌肉力与地反力。

肌肉动力学

肌肉激活动力学

尽管各个文献里的描述不尽一致，但一般来说excitation用来描述神经信号，是肌肉的控制输入（整流滤波归一化的EMG信号就是这个），而activation用来表示肌肉的活跃度，是肌肉的内部状态。在静态状态下，这两者可能相等，但在动态变化中，两者一般不同，存在由excitation到activation动态过程，即肌肉激活动力学。各方学者建立了不同的激活动力学模型，一般包含一个微分方程和一些非线性变换，其中Winters开发的一种激活动力学如下：

f=0.5tanh⁡(b(e−a))f=0.5 \tanh (b(e-a)) f=0.5tanh(b(e−a))
dadt=[1τa(0.5+1.5a)(f+0.5)+0.5+1.5aτd(−f+0.5)](e−a)\frac{d a}{d t}=\left[\frac{1}{\tau_{a}(0.5+1.5 a)}(f+0.5)+\frac{0.5+1.5 a}{\tau_{d}}(-f+0.5)\right](e-a) dtda=[τa(0.5+1.5a)1(f+0.5)+τd0.5+1.5a(−f+0.5)](e−a)

肌肉收缩动力学

肌肉收缩的特性可以用如上的Hill肌肉模型来描述，它由肌肉收缩单元CE与被动弹性单元PE并联，并与一个肌腱弹性单元T串联，肌肉与肌腱之间还存在一个角度成为羽状角。由于根据上面的模型，肌肉实际产生的力不仅受到收缩单元CE的激活度的影响，还受到肌肉的长度、收缩速度的影响。由肌肉纤维产生的力具有如下关系：

Fmt(t)=Ft=Fmax⁡[a(t)f(l)f(v)+fp(l)]cos⁡(ϕ(t))\begin{aligned} F^{\mathrm{mt}}(t) &=F^{\mathrm{t}} \\ &=F^{\max }\left[a(t)f(l) f(v) +f_{\mathrm{p}}(l)\right] \cos (\phi(t)) \end{aligned}Fmt(t)=Ft=Fmax[a(t)f(l)f(v)+fp(l)]cos(ϕ(t))

Fmax⁡F^{\max }Fmax表示最大等长收缩力，f(l),f(v),f(p)f(l),f(v),f(p)f(l),f(v),f(p)分别表示力-长度、力-速度、被动力-长度的关系。这些函数由下图所示的实验曲线所描述，一般可以通过多项式函数拟合来近似表示。

注意到上式公式中，实际输入为激活度aaa、归一化肌肉长度lll、和归一化收缩速度vvv，而肌肉长度与收缩速度之间存在关系，收缩动力学的状态量仅有aaa和lll，且上述方程描述的是一个微分方程。我们将整个肌肉动力学表示为如下形式，但这种表示并不严谨，实际上每块肌肉包含着两个微分方程：

Fmuscle=Fmuscle(q,a,a˙,LCE,L˙CE,e)\mathbf{F}_{muscle} = \mathbf{F}_{muscle}(\mathbf{q},\mathbf{a} , \mathbf{\dot a},\mathbf{L}_{CE},\mathbf{\dot L}_{CE},\mathbf{e})Fmuscle=Fmuscle(q,a,a˙,LCE,L˙CE,e)

地反力模型

OpenSim中采用Hunt-Crossley模型来表示接触力，接触力由刚度项、阻尼项和摩擦项组合而成，其表达式为：

fcontact=fstiffness+fdissipation+ffriction\mathbf{f}_{\text {contact}}=\mathbf{f}_{\text {stiffness}}+\mathbf{f}_{\text {dissipation}}+\mathbf{f}_{\text {friction}}fcontact=fstiffness+fdissipation+ffriction

接触力模型中的第一项表示弹性力。如上图所示，我们可以用两个弹性小球来描述接触，当小球接触时会发生弹性变形，变形的程度与各自的弹性模量EEE有关，总的变形为xxx。当变形发生时，每个小球都会受到一个弹性力：

fstiffness=fHz=(43σR1/2E∗)x3/2f_{\text {stiffness}}=f_{H z}=\left(\frac{4}{3} \sigma R^{1 / 2} E^{*}\right) x^{3 / 2}fstiffness=fHz=(34σR1/2E∗)x3/2

RRR表示合成相对曲率，E∗E^*E∗表示合成弹性模量，它们都可以通过两个小球的参数计算出来；σ\sigmaσ表示偏心系数，当σ=1\sigma=1σ=1的时候为圆接触，xxx为总的变形量。

接触力模型中的第二项表示损耗力，或是阻尼力。碰撞的发生往往伴随着能量的损耗，单纯的弹性力无法表示这一点，因此引入Hunt-Crossley耗散力：

fdissipation=32fstiffnessc∗x˙f_{dissipation}=\frac{3}{2} f_{stiffness} c^{*} \dot{x}fdissipation=23fstiffnessc∗x˙

c∗c^{*}c∗表示等效耗散系数，而x˙\dot{x}x˙表现出耗散力的阻尼特性。刚度力和耗散力都垂直接触面。

fnormal=fstiffness+fdissipationf_{\text {normal}}=f_{\text {stiffness}}+f_{\text {dissipation}}fnormal=fstiffness+fdissipation

接触力模型中的最后一项表示摩擦力，方向平行接触面：

ffriction=μ(v)fnormalf_{\text {friction}}=\mu(v) f_{\text {normal}}ffriction=μ(v)fnormal

μ(v)\mu(v)μ(v)表示等效摩擦因数，是关于相对运动速度的函数。

在使用上述的Hunt-Crossley模型来计算地反力时，地面的刚度认为无穷大，则相对弹性变形xxx可以通过接触小球圆心相对地面的距离计算得到，即：

Fcontact=Fcontact(q,q˙)\mathbf{F}_{contact} = \mathbf{F}_{contact}(\mathbf{q}, \dot{\mathbf{q}})Fcontact=Fcontact(q,q˙)

综合模型

综上所述，肌骨模型的动力学可以统一表示为：

M(q)⋅q¨+B(q,q˙)=J(q)TFcontact(q,q˙)+R(q)Fmuscle(q,a,a˙,LCE,L˙CE,e)\mathbf{M}(\mathbf{q}) \cdot \ddot{\mathbf{q}}+\mathbf{B}(\mathbf{q}, \dot{\mathbf{q}})=\mathbf{J(q)}^T \mathbf{F}_{contact}(\mathbf{q}, \dot{\mathbf{q}})+\mathbf{R(q)}\mathbf{F}_{muscle}(\mathbf{q},\mathbf{a} , \mathbf{\dot a},\mathbf{L}_{CE},\mathbf{\dot L}_{CE},\mathbf{e})M(q)⋅q¨+B(q,q˙)=J(q)TFcontact(q,q˙)+R(q)Fmuscle(q,a,a˙,LCE,L˙CE,e)

定义状态变量：

x=(q,q˙,LCE,a)T\mathbf{x}=\left(\mathbf{q}, \dot{\mathbf{q}}, \mathbf{L}_{C E}, \mathbf{a}\right)^{T}x=(q,q˙,LCE,a)T

在定义uuu为16块肌肉的控制量（excitations），上面的动力学方程可表示为：

f(x,x˙,u)=0\mathbf{f}(\mathbf{x}, \dot{\mathbf{x}}, \mathbf{u})=0f(x,x˙,u)=0

肌骨模型的仿真应用

前向仿真

逆向分析

逆向分析是指，在已有运动学、地反力数据后，根据实验数据推断出产生相应运动的肌肉活动，也就是OpenSim中的CMC。由于肌肉之间存在和比较严重的耦合问题，多块肌肉连接同一个关节，同一块肌肉连接着不同关节，因此逆向分析最重要的是解决肌肉冗余 - Solving the Muscle Redundancy Problem。

在这个问题中，已有运动学数据qk\mathbf{q}_kqk和地反力FCk\mathbf{F}_{\mathbf{C}k}FCk序列，k=1,…,Kk=1, \ldots, Kk=1,…,K表示离散时间。我们先计算逆动力学得到关节力矩TIDkT_{\mathbf{ID}k}TIDk，关节力矩应当等于有肌肉收缩的等效力矩：

TIDk=∑m=1MrmkFTm+eTkTmaxT_{\mathrm{ID} k}=\sum_{m=1}^{M} r_{m k} F_{\mathrm{T} m}+e_{\mathrm{T} k} T_{\mathrm{max}}TIDk=m=1∑MrmkFTm+eTkTmax

mmm表示肌肉序号，rrr表示肌肉对关节等效力矩。有时候关节力矩和肌肉里橘无法严格相等，因此我们添加了公式中的最后一项来表示力矩的残差，而eTke_{\mathrm{T} k}eTk表示残差系数，在后面我们会设计优化目标是残差力尽可能的小。

由此建立逆向分析的优化框架：