前言
对于角色动画，目前的技术主要包括这两种方式：基于运动学模拟的动画，和基于动力学模拟、即基于物理模拟的动画。第一种方式的应用相当广泛，我们接触的游戏大部分都采用了基于运动学的方式去实现动画，这种方式性能优秀，便于实现，然而却缺少了许多真实性。基于物理模拟的方式来实现动画，对性能的要求比较高，且实现起来比较复杂，优点是交互效果更真实。对人体的肌肉骨骼进行模拟，属于物理模拟的动画实现。
肌肉的模拟较为复杂，涉及到生物力学这门交叉学科的内容。目前较为通用的肌肉模拟的模型是希尔模型(Hill model)。接下来我会根据一篇2013年的TOG论文对肌肉模型的模拟进行介绍。
本篇博客参考自《Flexible Muscle-Based Locomotion for Bipedal Creatures》。
肌肉骨骼模型
希尔模型将肌肉表示成三元素结构，它由下列要素组成：
- 一个收缩单元(contractile element, CE)，表示根据肌肉激活状态收缩的肌肉纤维，这个单元可以产生主动力。
- 一个并行弹性单元(parallel elastic element, PEE)，表示肌肉纤维周围的被动弹性组织，这个单元可以产生被动力。
- 一个串行弹性单元(parallel elastic element, SEE)，代表连接肌肉和骨骼的肌腱，这个单元可以产生被动力。
CE产生的力 F C E F_{CE} FCE，取决于肌肉的恒定最大等距力 F m a x F_{max} Fmax、肌肉激活 a a a、肌肉纤维的长度 L C E L_{CE} LCE，以及收缩速度 V C E V_{CE} VCE：
F C E = a F m a x f L ( L C E ) f V ( V C E ) \begin{align} F_{CE} = aF_{max}f_L(L_{CE})f_V(V_{CE}) \end{align} FCE=aFmaxfL(LCE)fV(VCE)
其中， f L f_L fL表示力与肌肉长度之间的关系， f V f_V fV表示力与当前收缩速度之间的关系，下图即为这两个函数：
弹性元件PEE和SEE产生的被动力 F P E E F_{PEE} FPEE和 F S E E F_{SEE} FSEE根据其长度可以建模为非线性弹簧：
F S E E = f S E E ( L M − L C E ) F P E E = f P E E ( L C E ) \begin{align} F_{SEE} &= f_{SEE}(L_M-L_{CE}) \\ F_{PEE} &= f_{PEE}(L_{CE}) \end{align} FSEEFPEE=fSEE(LM−LCE)=fPEE(LCE)
其中 f S E E f_{SEE} fSEE和 f P E E f_{PEE} fPEE是非线性的力-长度关系， L M L_M LM为肌肉的总长度。
由于SEE与CE和PEE串联在一起，因此总肌肉力 F M F_M FM服从力平衡方程：
F M = F C E + F P E E = F S E E \begin{align} F_M = F_{CE} + F_{PEE} = F_{SEE} \end{align} FM=FCE+FPEE=FSEE
SEE的初始化长度为其最优长度， L C E o p t L_{CE}^{opt} LCEopt。结合肌腱的松弛长度 L S E E s l a c k L_{SEE}^{slack} LSEEslack，可以定义肌肉的休息长度 L M r e s t L_M^{rest} LMrest:
L M r e s t = L S E E s l a c k + L C E o p t \begin{align} L_M^{rest} = L_{SEE}^{slack} + L_{CE}^{opt} \end{align} LMrest=LSEEslack+LCEopt
在实际进行模拟的时候，输入参数为激活状态 a a a以及肌肉的总长度 L M L_M LM。
激活状态的动态更新
肌肉的激活状态 a a a是一个相对缓慢的电化学过程的结果，该过程基于控制系统输出的神经兴奋信号 u u u。这个过程称为激活动力学，其模型为：
∂ a ∂ t = c a ( a − u ) \begin{align} \frac{\partial a}{\partial t} = c_a(a-u) \end{align} ∂t∂a=ca(a−u)
其中 c a c_a ca是恒定的激活和失活速率。
肌肉几何以及肌肉与骨骼的相互作用
骨骼和肌肉具有双向的作用：肌肉施加改变骨骼姿势的力量，而骨骼姿势决定肌肉的长度 L M L_M LM，进而影响收缩动力。整个肌肉的路径是由它的肌腱附着在骨头上的位置、肌肉环绕的骨性标志以及随时间变化的肌肉跨越的关节决定的。在原文的模型中，将肌肉路径定义为一组线段。这个模型是一种简化，具有高性能的优势，并且省略了对肌肉的几何形状建模的需要。如下图，长方体代表着肌肉附着的骨骼，红色连线代表着对肌肉的简化：
任何肌肉 M M M的路径都是由 n n n个附着点数组定义的， [ { b 1 , p 1 } . . . , { b n , p n } ] \left[ \left\{ b_1, p_1 \right\}...,\left\{b_n,p_n\right\} \right] [{b1,p1}...,{bn,pn}]，每一个附着点由一个偏移量 p i \mathbf p_i pi以及一个附着在的骨骼 b i b_i bi定义。偏移量 p i \mathbf p_i pi定义为物体 b i b_i bi坐标系中的一个固定偏移量，它和它附着着的骨骼一起移动和旋转。多个附着点可以附着在一个骨骼上。第一个点和最后一个点 p 1 p_1 p1和 p n p_n pn表示肌肉肌腱附着在骨骼上的位置，其他点都是通过点。
肌肉的总长度 L M L_M LM等于 n − 1 n - 1 n−1个肌肉段长度之和， [ s 1 , . . . , s n − 1 ] [s_1,...,s_{n-1}] [s1,...,sn−1]，这些肌肉段的长度由世界坐标系中每个点的位置 p i W p_i^W piW求出：
L M = ∑ i = 1 n − 1 ∣ ∣ s i ∣ ∣ , s i = p i + 1 W − p i W \begin{align} L_M = \sum_{i=1}^{n-1} || s_i ||, s_i = p_{i+1}^W - p_i^W \end{align} LM=i=1∑n−1∣∣si∣∣,si=pi+1W−piW
这些附着点在这篇文章的定义中，可以有四种活动方式：要么是固定的，要么被限制在一条线、一个平面或一个盒子上：
肌肉如何产生作用
由肌肉总长度 L M L_M LM以及激活状态 a a a可以计算肌肉收缩力 F M F_M FM。由于受力平衡，这个力被传递到肌肉在骨骼上的每个附着点上，并在它跨越的每个关节上产生扭矩。对于每个关节 k k k，在力臂 r k r_k rk定义的方向上可以产生一个扭矩 τ k \tau_k τk。这个力臂对应于穿过关节的肌肉段方向 s c s_c sc与关节中心 j k j_k jk到附着点 p c W p_c^W pcW的矢量之间的叉积：
τ k = F M ∣ ∣ r k ∣ ∣ , r k = ( p c W − j k ) × s c ∣ ∣ s c ∣ ∣ \begin{align} \tau_k = F_M||r_k||, r_k = (p_c^W - j_k) \times \frac{s_c}{||s_c||} \end{align} τk=FM∣∣rk∣∣,rk=(pcW−jk)×∣∣sc∣∣sc
在Unity中实现肌肉骨骼模拟
在unity资源商店中有一个资源，实现了上述的肌肉骨骼模型：
https://assetstore.unity.com/packages/tools/physics/kinesis-physical-muscle-model-based-movement-206089
比较贵，$69.99，忍痛割刀买了。接下来讲解一下这个资源里是如何实现肌肉骨骼模型的。
在这个资源里，实现了调节肌肉激活状态，驱动人物的关节运动的物理模拟。
底层模拟基础
为了实现添加扭矩，驱动关节移动和旋转的功能，利用unity提供的Rigidbody以及Character Joint组件实现底层的模拟。
模拟实现
在每一个时步的模拟中，首先，参考公式 ( 6 ) (6) (6)，根据输入的激活值更新激活状态 a a a。然后，更新SEE部分的长度（如果是初始调用，则不更新）。根据SEE的长度计算归一化后CE的长度，以及归一化后CE的收缩速度（如果是初始调用，则为0）。根据激活状态 a a a、归一化后CE的长度、归一化后CE的收缩速度，计算CE部分产生的主动力，参考公式 ( 1 ) (1) (1)。同时，根据CE的长度，计算PEE部分的被动力。CE部分产生的主动力，加上PEE部分产生的被动力，等于SEE部分产生的被动力，也等于肌肉产生的力，参考公式 ( 4 ) (4) (4)。其中，SEE部分产生的被动力，可以用来更新SEE部分的长度。接下来，利用肌肉产生的力，计算肌肉穿过的关节的扭矩，参考公式 ( 8 ) (8) (8)。
补充
之前在实习时，leader让我做过这样一件事情，根据人物的运动序列，反向计算肌肉的激活状态，当时没能实现。现在想想，可以这样实现：
1、从人物的运动学模型转换到动力学模型，通过PD微分，求解人物每个关节所需要的扭矩。
2、根据每块肌肉和关节的关系，计算出每块肌肉应该产生的发力，通过每块肌肉产生的发力，反求出肌肉的激活状态。（根据公式 ( 1 ) (1) (1)，并对激活状态的反推做简化，用 a = F C E / F m a x a = F_{CE} / F_{max} a=FCE/Fmax计算）
参考
Geijtenbeek T, Van De Panne M, Van Der Stappen A F. Flexible muscle-based locomotion for bipedal creatures[J]. ACM Transactions on Graphics (TOG), 2013, 32(6): 1-11.
https://assetstore.unity.com/packages/tools/physics/kinesis-physical-muscle-model-based-movement-206089