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

肌肉骨骼模型

骨骼的刚体动力学模型

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

首先基于刚体动力学对于人体骨骼系统进行建模。首先确定系统的广义坐标 $q$ ，对于图中所示的二维模型为例，广义坐标包含躯干的水平、数值方向位移和姿态角，以及每条腿三个关节角，一共 $3+2*3=9$ 个自由度。之后通过这些广义变量和广义变量的导数，表示出系统的动能和势能，并构造拉格朗日函数：
$L=T-V$
考虑骨骼系统仅收到地反力和肌肉力的作用，通过拉格朗日方程得到系统的动力学模型：
$\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}$

$\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}$

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

肌肉动力学

肌肉**动力学

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

$f=0.5 \tanh (b(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)$

肌肉收缩动力学

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

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

$\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}$

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

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

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

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

地反力模型

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

$\mathbf{f}_{\text {contact}}=\mathbf{f}_{\text {stiffness}}+\mathbf{f}_{\text {dissipation}}+\mathbf{f}_{\text {friction}}$

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

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

$f_{\text {stiffness}}=f_{H z}=\left(\frac{4}{3} \sigma R^{1 / 2} E^{*}\right) x^{3 / 2}$

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

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

$f_{dissipation}=\frac{3}{2} f_{stiffness} c^{*} \dot{x}$

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

$f_{\text {normal}}=f_{\text {stiffness}}+f_{\text {dissipation}}$

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

$f_{\text {friction}}=\mu(v) f_{\text {normal}}$

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

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

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

综合模型

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

$\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})$

定义状态变量：

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

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

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

肌骨模型的仿真应用

前向仿真

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

逆向分析

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

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

$T_{\mathrm{ID} k}=\sum_{m=1}^{M} r_{m k} F_{\mathrm{T} m}+e_{\mathrm{T} k} T_{\mathrm{max}}$

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

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

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