技术领域
[0001] 本
发明属于
虚拟现实技术领域，尤其涉及一种流体环境下的关节动画建模方法。
背景技术
[0002] 关节动画属于
角色动画，是一直是虚拟现实技术研究的热点之一。它综合利用计算机科学、艺术、数学、物理学和其它相关学科的知识，在计算机上生成绚丽多彩的、连续的、真实的虚拟画面，给人们提供了一个充分展示个人想象
力和艺术才能的新天地。关节动画使用关节骨架来表示人类或者其他骨架动物的身体结构，是动画驱动技术中最主要的思想。虽然计算机动画在许多领域占据着越来越重要的角色，但是许多问题仍未很好解决。
[0003] 基于物理模型的关节动画技术是八十年代后期发展起来的一种新的计算机动画技术。这种建模方法考虑了物体在真实世界中的属性，如它具有
质量、转动惯矩、弹性、
摩擦力等，并采用动力学原理来产
生物体运动。计算机动画设计者不必关心物体运动过程的细节，只需确定物体运动所需的一些物理属性及一些约束关系，如质量、形状、外力等。经过近几年的发展，它已在图形学中成为一种具有潜在优势的三维造型和运动模拟技术。尽管该技术比传统动画技术的计算复杂度要高得多，但它能逼真地模拟各种自然物理现象，能处理诸如重力、
风、
碰撞检测等在内的复杂动力学模型，这是基于几何的传统动画生成技术所无法比拟的。
[0004] 虽然利用
现有技术能够实现一定程度上的交互，但是复杂场景作用下的关节动画却少有人研究，这是由于受到本身动力学复杂的模型与计算机软
硬件条件的限制。流体环境下的关节动画技术在电影动画，生物游泳力学，视频游戏以及潜
水机器人领域等方面都具有重要作用。此外，基于流体环境的关节动画研究也在一定程度上为其他复杂受力场景提供了关节动画的解算依据。
发明内容
[0005] 为了模拟真实的复杂场力作用下的关节动画，本发明提供了一种流体环境下关节动画建模方法。本发明分别对
铰链体的驱动、动力学和受力作用进行建模。在数据驱动上，使用一种自主发明的基于计算力矩
控制器，提高了
跟踪轨迹的准确性和
稳定性。在动力学上，采用拉格朗日动力学进行建模，减少了计算变量的开销，进而提高运算效率。在外力作用上，将流体对铰链体的外力分为法线和切线方向，进行分别讨论求解，获取更加真实的受力效果。
[0006] 本发明一种流体环境下的关节动画建模方法，该方法包括步骤：
[0007] 步骤1，输入上一时刻计算所得的所有广义
位置和广义速度，基于广义
坐标系；
[0008] 步骤2，计算驱动函数，获得期望的广义位置和广义
加速度，基于广义坐标系；
[0009] 步骤3，通过控制器，获取铰链体的广义驱动力，基于广义坐标系；
[0010] 步骤4，根据拉格朗日动力学，计算铰链体的广义加速度，基于广义坐标系；
[0011] 步骤5，使用欧拉法，计算下一时刻广义位置和中间广义速度，基于广义坐标系；
[0012] 步骤6，计算流体对铰链体产生的外力，包括压力和摩擦力，并从笛卡尔坐标系变换至广义坐标系；
[0013] 步骤7，计算外力产生的广义加速度，基于广义坐标系；
[0014] 步骤8，使用欧拉法，更新广义速度；
[0015] 步骤9，变换坐标系，将铰链体从广义坐标系变换至笛卡尔坐标系，并计算铰链体在笛卡尔坐标系下的位置和方位角；
[0016] 步骤10，根据铰链体的位置和方位角进行
渲染，基于笛卡尔坐标系；
[0017] 所述的步骤1中输入上一时刻计算所得的所有广义位置和广义速度，其上一时刻的广义位置来源于步骤5中，使用欧拉法计算得到的广义位置。上一时刻的广义速度来源于步骤8中，使用欧拉法更新广义速度。第一时刻的广义位置和广义速度进行初始化设置。这里的广义位置指铰链体关节的可控角，广义速度指该可控角的
角速度。
[0018] 所述的步骤2中计算驱动函数，获得期望的广义坐标和广义加速度。该驱动函数采用周期函数，如公式(1a)和(1b)，针对每个广义位置和广义速度进行设置。
[0019] (1a)
[0020] (1b)
[0021] 其中， 、 、 和 分别为振幅、周期、
相位和偏移量。
[0022] 所述的步骤3中计算力矩控制器的计算式如公式(1)所示：
[0023] (1)
[0024] 其中， 和 分别为下一时刻期望到达的广义坐标和广义加速度， 、 和分别为当前时刻广义坐标、广义速度和广义加速度；另外， 、 和 为控制系数，并满足且 与 在同一数量级；此外，根据拉格朗日动力学知， 为下一时刻的质量项， 为当前时刻的科氏力项，使用欧拉法化简 项，如公式(2)和公式(3)所示：
[0025] (2)
[0026] (3)
[0027] 其中， 为当前时刻的雅克比矩阵， 为当前时刻雅克比矩阵的转置， 为当前时刻的质量矩阵，包括质量和惯性张量两部分构成， 为角速度 的反对称矩阵；与标准的拉格朗日动力学不同的是，这里使用的是 和 来近似计算下一时刻产生的项。
[0028] 该控制器能在欧拉法的作用下，既能表现出稳定性，又能表现出跟踪性。下面将进行论证。
[0029] 令 和 分别代表部分参数，如公式(a)和公式(b)所示，
[0030] (a)
[0031] (b)
[0033] (c)
[0034] 论证：满足 ，该控制器在欧拉法中能达到稳定性和跟踪性。
[0035] 为了计算方便，我们假设 ，并代入公式(c)，得
[0036] (d)
[0037] 上下同时除以 ，假设 趋于无穷大，化简得，
[0038] (e)
[0039] 使用欧拉方法分别求得下一时刻的速度 ，下一时刻的位置 和下下时刻的位置 ，可得公式(f) (g) (h)：
[0040] (f)
[0041] (g)
[0042] (h)
[0043] 可以看出该控制器在欧拉方法的作用下，若设置 值比 值过小，则在下两个时间步后获得的广义坐标几乎与预先设定的运动轨迹一致，即稳定性，若设置的 值与 值在同一数量级中，则可以在下两个时间步后，获得预先设定的运动轨迹与拉格朗日动力学中速度反馈产生的
叠加效果，即跟踪性。
[0044] 所述的步骤4中根据拉格朗日动力学，计算铰链体的广义加速度，基于广义坐标系,具体为：
[0045] 步骤4.1，计算铰链体的质量矩阵 ，如公式(4)所示。
[0046] (4)
[0047] 步骤4.2，计算铰链体的科氏力矩阵 ，如公式(3)所示。
[0048] 步骤4.3，计算铰链体的广义重力 ，如公式(5)所示。
[0049] (5)
[0050] 步骤4.4，计算铰链体的广义加速度 ，如公式(6)所示。
[0051] (6)
[0052] 其中， 为当前时刻的内力。
[0053] 所述步骤5使用欧拉法，计算下一时刻每个广义位置和中间广义速度，基于广义坐标系，其具体为：
[0054] 步骤5.1，使用欧拉法计算下一时刻铰链体的位置 ，如公式(7)所示：
[0055] (7)
[0056] 步骤5.2，使用欧拉法计算下一时刻铰链体的速度 ，如公式(8)所示：
[0057] (8)。
[0058] 所述步骤6计算流体对铰链体产生的外力，包括压力和摩擦力，并从笛卡尔坐标系变换至广义坐标系，其具体为：
[0059] 步骤6.1，通过内外
体素化的方法，确定流体与铰链体的耦合面；
[0060] 步骤6.2，根据纳维耶斯托克方程，获得流体对铰链体产生的压力；
[0061] 步骤6.3，计算流体对铰链体产生的摩擦力 ，包括流体产生的粘性摩擦力和库伦摩擦力，使用如公式(9)所示的摩擦力模型；
[0062] (9)
[0063] 其中， 为位置信息向量，包括xyz三个坐标， 为该位置上的粘性摩擦力系数， 为 位置上的流体速度， 为 位置上的库伦摩擦力， 为该位置上的切线方向；另外， 函数如公式(10)所示；
[0064] (10)
[0065] 步骤6.4，将流体产生的压力和摩擦力从笛卡尔坐标系变换至广义坐标系，并计算流体产生的广义外力 ，如公式(11)所示；
[0066] (11)
[0067] 其中， 为每个体素上的片面微元， 为笛卡尔坐标系和广义坐标系的坐标系变换矩阵(3*m，m代表
自由度)。
[0068] 所述步骤7计算外力产生的广义加速度，基于广义坐标系，具体为：通过
牛顿第二定律，计算当前时刻流体外力产生的广义加速度 ，如公式(12)所示：
[0069] (12)。
[0070] 所述步骤8使用欧拉法，更新广义速度；具体为：根据当前时刻流体外力产生的广义加速度 ，更新广义速度 ，如公式(13)所示：
[0071] (13)。
[0072] 所述步骤9中将铰链体从广义坐标系变换至笛卡尔坐标系，并计算铰链体在笛卡尔坐标系下的位置和方位角；具体为：位置的变换如公式(14)所示，方位角的变换如公式(15)所示：
[0073] (14)
[0074] 其中， 为当前铰链在笛卡尔坐标系下的位置向量， 为该铰链所连接的父链在笛卡尔坐标系下的位置向量， 为当前铰链基于父链的旋转矩阵， 为当前铰链的关节可控角的集合， 为该铰链相对于父链的位置向量；
[0075] (15)
[0076] 其中， 为当前铰链基于笛卡尔坐标系的旋转矩阵，分别为所有铰链的关节可控角的集合，n代表当前铰链，n-1代表当前铰链相连的父链，依次类推。
[0077] 本发明的有益效果是：在数据驱动方面，使用基于计算力矩控制器，用以获取铰链体的驱动力，这种控制器除了具有调整每个关节沿着预先设定的轨迹运动的作用之外，还加入了上一时刻铰链体的内力对当前时刻惯性的影响及期望的加速度影响，很好地达到了控制器的跟踪作用，并且配合欧拉方法使用时能达到稳定性（见附加论证）。在动力学方面，采用拉格朗日动力学进行建模，这种广义坐标系的建模方式比起牛顿欧拉动力学来说，它通过添加了内部的约束条件，减少了未知变量，提高了计算速度。在外力作用方面，该发明将外力分为法线和切线方向进行讨论，法向为流体对铰链体的压力，切向为流体对铰链体的摩擦力，并根据纳维耶斯托克方程求解外力，这种方法可以获取更加真实的受力效果。
附图说明
[0080] 图3 示出了外部配置文件需要保存的信息；
[0081] 图4 示出了本实施实例所使用的铰链连接系统；
[0082] 图5 示出了本实施实例中铰链的形状和信息；
[0083] 图6 示出了本实施实例使用的体素化。
具体实施方式
[0084] 下面结合附图和
实施例对本发明优先实施方式进一步说明：
[0085] 实施例：
[0086] 图1所示的流程图给出了本发明整个实施的具体过程：
[0087] 根据算法的流程，本实施搭建了如图2所示框架进行解算。其中分别建立了数据
驱动器、铰链体解算器、拉格朗日解算器、计算力矩控制器、欧拉解算器、外力解算器和外力作用解算器，并将这些解算模
块通过一个可以实现上传和下载的数据管理器进行连接，实现数据传输并达到松耦合的效果。在下述步骤1-8中将分别对这些模块进行说明。
[0088] 在算法实施之前，本实施例中需要使用外部配置文件(如.xml)分别保存铰链、驱动函数和模型的信息作为运算的初始信息，如图3所示。这部分文件信息通过xml文件加载器从外部.xml文件中获得，并通过数据管理器的中转，分配到不同的模块中作为初始信息。此外，模型使用.obj文件。
[0089] 本实施需要预先设置时间步长度 、控制器的可控参数 、 和 ，且设置条件需满足 且 与 在同一数量级。
[0090] 本实施实例解算的是一个具有人体形状的铰链体，自由度为28，其中腰部为根
节点（自由度为6），如图4所示，可看作
树形结构。
[0091] 步骤 1，从数据管理器中取出上一
帧计算所得的广义位置和广义速度，作为输入条件。需注意，第一帧的初始条件来自.xml配置文件的信息。
[0092] 步骤 2，数据驱动器从数据管理中获得驱动函数的信息，计算出下一时刻期望的位置和加速度，并上传至数据管理器。本实施实例使用基本的周期函数作为数据驱动函数，如公式(16a)和(16b)所示：
[0093] (16a)
[0094] (16b)
[0095] 其中， 、 、 和 分别为振幅、周期、相位和偏移量。
[0096] 步骤 3，通过铰链体解算器、拉格朗日解算器和计算力矩控制器计算广义驱动力。进一步包括：
[0097] 步骤 3.1，铰链体解算器从数据驱动器中下载当前时刻的位置和速度信息，并计算每根铰链的雅克比矩阵和铰链的信息（如惯性张量和质量等）。具体做法是，从铰链体的根节点开始，按树的宽度遍历依次解算每个铰链体的局部雅克比矩阵和铰链的信息。
[0098] 步骤 3.1.1，计算每根铰链的局部雅克比矩阵,包括基于速度的雅克比矩阵、基于角速度的雅克比矩阵和求基于时间的一阶导的雅克比矩阵。
[0099] 首先，考虑单个铰链在笛卡尔坐标系下产生的速度，将其分解为线速度 和角速度 ，均为3*1的列向量。假设 和 分别为质心的位置（3*1向量）和关节体的旋转矩阵（3*3矩阵）。线速度 可以表示为公式(17)，
[0100] (17)
[0101] 其中， 和 为n*1的列向量，n为整个链上的自由度，其中 ， 为3*n的矩阵。该实施实例中的n即为22。此外，角速度可以表示为公式(18)，
[0102] (18)
[0103] 由于 恒为反对称矩阵，可以被表示为 （ 的反对称形式， 相当于雅克比矩阵 的第j列），所以角速度 可以用雅克比矩阵表示为公式(19)，
[0104] (19)
[0105] 根据公式(17)和公式(19)可知，每根铰链基于笛卡尔坐标系和广义坐标系转换的雅克比矩阵可以表示为公式(20)中的 ， 为6*n的矩阵。该实施实例中 为6*22的矩阵；
[0106] (20)
[0107] 本实施实例除了需要解算出每个铰链当前时间步的 之外，还需解得 和下一时刻的 。由于下一时刻的信息未知，本实施实例采用欧拉方法，如公式(21)所示，近似求得下一时刻的 ；
[0108] (21)
[0109] 因此针对每个关节的每个广义坐标，求解雅克比分量（3*1的列向量），并上载至数据管理器。其中，每个关节的自由度小于等于3。根据每个关节的旋转信息，可以将每个关节分解成基于父链坐标的局部的绕x、绕y和绕z轴的分量，则局部旋转矩阵的基本形式是可知的，如公式(22a) (22b) (22c)所示。
[0110] (22a)
[0111] (22b)
[0112] (22c)
[0113] 雅克比分量的具体求解方法：首先，计算局部雅克比线速度分量 ，如公式(23a) (23b) (23c)所示。其中， 为父链总的旋转矩阵（3*3的矩阵）,为质心基于父链的局部坐标位置（3*1的列向量）， 为基于广义位置的旋转矩阵， 为旋转矩阵基于广义位置 的一阶导数（根据公式(22a) (22b) (22c)易求得）。除此以外，还需计算一组基于整个链长度的线性速度分量，用于填充拉格朗日动力学中雅克比矩阵的线速度分量，方法同上。
[0114] (23a)
[0115] (23b)
[0116] (23c)
[0117] 其次，计算局部雅克比角速度的分量 。首先将分别初始化为： ， 和 ，则某个广义位置 的初始化向量记为 ，
具体向量与广义坐标的信息有关。其次，分三种情况讨论：
[0118] 1）假设每个关节可用自由度为1（即只存在 ），则 的值就等于初始化的值，如公式(24a)所示。
[0119] 2）假设每个关节可用自由度为2（即存在 和 ），则 的值如公式(24a)和公式(24b)所示。
[0120] 3）假设每个关节可用自由度为3（即存在 、 和 ），则的值如公式(24a)、公式(24b)和公式(24c)所示。雅克比角速度分量则为父链的旋转矩阵乘以 ；
[0121] (24a)
[0122] (24b)
[0123] (24c)
[0124] 再次，需要计算雅克比线速度的分量的导数 （基于时间），如公式(25a) (25b) (25c)所示，其中 是对父链总的旋转矩阵进行求导运
算，该项在宽度遍历时，依次累加求解可得；
[0125] (25a)
[0126] (25b)
[0127] (25c)
[0128] 最后，需要计算雅克比角速度的分量的导数 （基于时间），将 依然初始化为： ， 和 ，则某个广义位置
的初始化向量记为 ，具体向量与广义位置的信息有关。依次分三种情况
讨论：
[0129] 1）假设每个关节可用自由度为1（即只存在 ），则 的值，如公式(26a)所示。
[0130] 2）假设每个关节可用自由度为2（即存在 和 ），则 的值如公式(26a)和公式(26b)所示。
[0131] 3）假设每个关节可用自由度为3（即存在 、 和 ），则 的值如公式(26a)、公式(26b)和公式(26c)所示；
[0132] (26a)
[0133] (26b)
[0134] (26c)
[0135] 因为 为标准的旋转矩阵形式，所以可以预先设定其基于广义位置的一阶导数和二阶导数 的计算式。
[0136] 需要注意，这里除了需要求出当前时刻全部的雅克比分量，还应求得下一时刻的雅克比分量用于计算 ，下一时刻的位置和速度使用欧拉法近似求得。
[0137] 步骤 3.1.2，计算铰链的惯性张量和质量。假设本实施实例中使用的铰链为椭球体，如图4所示。根据平行轴定理知，每个铰链体绕其头部旋转时的惯性张量可表示为公式(27)；
[0138] (27)
[0139] 其中，m为椭球体的质量（计算公式为公式(28)所示），a、b和c分别为椭球体的长、宽和高，x、y和z为椭球体质心相对于旋转基轴的笛卡尔坐标；
[0140] (28)
[0141] 步骤 3.2，拉格朗日解算器从数据驱动器中下载当前时刻的位置和速度信息、雅克比矩阵以及铰链的信息，并计算质量项和科氏力项，完成后上传至数据驱动器。
[0142] 步骤 3.2.1，计算每个铰链的质量项 ，其中 为编号1至n的铰链中的某个铰链。将每个铰链信息中的质量和惯性张量排列为一个6*6的方阵，如公式(29)所示。该实施实例中，n为22；
[0143] (29)
[0144] 其中， 为3*3的单位矩阵。将 和 按公式(20)的排列方式进行上下排列，则根据公式(4)可得到每个铰链的质量项。
[0145] (4)
[0146] 步骤 3.2.2，计算每个铰链的科氏力项。将 、 、 和 依然按照公式(20)的排列方式进行排列获得 和 ，则根据公式(3)可得每个铰链的科氏力项。
[0147] (3)
[0148] 步骤 3.3，计算力矩控制器从数据驱动器中获取质量项、科氏力项、当前时刻的位置、当前时刻的速度、下一时刻期望的位置和下一时刻期望的加速度，并根据公式(1)计算每个铰链的驱动力，并上传数据管理器。
[0149] (1)
[0150] 步骤 4，拉格朗日解算器从数据驱动器中下载每个铰链广义驱动力，然后通过拉格朗日动力学公式，解得当前时刻的加速度，并上传至数据管理器。具体实施如下：
[0151] 步骤 4.1，计算铰链体的质量矩阵 。这里的 为整个铰链体的质量矩阵，规格是n*n（n为自由度）的方阵，是步骤 3.2.1的所有铰链的质量项排列的总和。公式(4)中，为6k*6k的方阵，将每个铰链的质量项排列成如公式(30)所示的形式，k是铰链体关节的总数。 为6k*n的矩阵，如公式(31)所示，其中n为自由度。 和 分别为3k*n矩阵， 的排列形式如公式(32)所示。 的矩阵排列与 相同。该实施实例中的n为22，k为10；
[0152] (30)
[0153] (31)
[0154] (32)
[0155] 步骤4.2，计算铰链体的科氏力矩阵 ，如公式(3)所示。这里的 为整个铰链体的科氏力矩阵，规格是n*n的方阵，是步骤 3.2.2的所有铰链的科氏力项排列的总和。其中矩阵为6k*n的矩阵，参考公式(31)的排列方式。 为角速度 的反对称矩阵，规格为6k*6k的方阵，如公式(33)所示。令角速度 的三个分量为 、 和 ，则 写为公式(34)。分别依次排列所有矩阵的分量，根据公式(3)可以解得 。该实施实例中的n为22，k为
10。
[0156] (33)
[0157] (34)
[0158] 步骤 4.3，计算铰链体的广义重力 ，如公式(5)所示。其中 是3*n的雅克比矩阵，如公式(35)所示。 是3*1的列向量，表示重力，该实施实例中的n为22；
[0159] (5)
[0160] (35)
[0161] 步骤4.4，计算铰链体的广义加速度 ，如公式(6)所示。拉格朗日动力学方程，相当于一个求解具有n个未知数的方程组，n为自由度。使用LU分解和回代的方法，可解得线性方程组的值。该实施实例中的n为22。
[0162] (6)
[0163] 步骤5，欧拉解算器从数据管理器中下载当前时刻的位置、速度和加速度，并计算下一时刻的位置和中间速度，并上传至数据管理器，计算如公式(7)和公式(8)所示。
[0164] (7)
[0165] (8)
[0166] 步骤 6，外力解算器从数据管理器上下载当前时刻每根铰链的雅克比矩阵信息，然后计算流体对铰链体产生的外力，包括压力和摩擦力，并从笛卡尔坐标变换至广义坐标系，最后将计算的广义外力上传至计算数据管理器。进一步包括：
[0167] 步骤 6.1，通过“内外体素化”的方法，确定流体与铰链体的耦合面，作为铰链体的受力面。本实施实例在GPU运算下进行，使用矩阵投影把输入的.obj模型的三角形网状结构，在目的3D纹理结构的每一个切片里渲染一次，如图6所示（左为模型，右为渲染切片）。当绘制几何体时，使用模板缓存(与切片的维数相同)，并将模板缓存初始化为0。结果是网状结构内部的所有体素得到了一个非零的模板值。然后，做一次最终遍历，把模板值复制到障碍物的纹理结构中去。这样就能区分三种类型的单元：内部单元（非零模板值）、外部单元（模板值为0）和紧贴边界的内部单元（取离非零模板值最近的0值）。
[0168] 步骤 6.2，根据纳维耶斯托克方程，获得流体对铰链体产生的压力。具体实施步骤如下：
[0169] 步骤 6.2.1，本实施实例的流
体模拟在GPU下进行，采用MAC网格进行模拟。由于压力值和位置信息存储在3D纹理中，根据耦合面的体素信息，可提取该部分压力信息 和对应的位置信息 。
[0170] 步骤6.2.2，铰链体外设
包围盒，用以检测压力作用在铰链的位置。具体方法是检测提取的压力所对应的位置是否在包围盒内，如果在某根铰链所属的包围盒内，则该压力值作用于此铰链，如果不在该包围盒之内，则判定为不会作用于此铰链。
[0171] 步骤6.2.3，因铰链采用标准的椭球体模型，根据对应压力所在铰链体的位置，易求出铰链体当前位置的法向量 ，得到压力的作用分量 。
[0172] 步骤 6.3，计算流体对铰链体产生的摩擦力，其中包括流体产生的粘性摩擦力和库伦摩擦力，我们如公式(9)所示的摩擦力模型。其具体实施步骤如下。
[0173] (9)
[0174] 步骤 6.3.1，在步骤6.2.1的
基础上提取该部分的速度信息 ，并参考6.2.2的步骤检测速度作用在铰链上的位置。
[0175] 步骤6.3.2，根据对应速度所在铰链体的位置，求出铰链体当前位置的切向量 ，根据公式(7)，可得到摩擦力的作用分量 。
[0176] 步骤 6.4，计算广义外力并上传至数据管理器，具体步骤如下：
[0177] 步骤6.4.1，计算雅克比矩阵。铰链体使用的是标准的椭球体模型，因此易求得法向量和切向量。这里的雅克比矩阵 形式为公式(35)所示。
[0178] 步骤6.4.2，进行数值积分。采用龙贝格求积公式，根据公式(11)所示，针对每个广义外力进行数值积分，积分的范围是从0至整个铰链体的面积，逐次分半的细分粒度为体素的面积。若为了提高效率可减少逐次分半的次数，近似求解。
[0179] (11)
[0180] 步骤 7，外力作用解算器从数据管理器上下载每个广义位置的驱动力和质量信息，并计算广义加速度，如公式(12)所示。这里的 即为当前广义位置所在铰链的惯性张量。
[0181] (12)
[0182] 步骤 8，外力作用解算器使用欧拉法，更新广义速度，并上传至数据管理器，作为下一时刻输入的条件，如公式(13)所示。
[0183] (13)
[0184] 步骤 9，变换坐标系，将铰链体从根节点开始，按广度进行遍历，分别求出基于父链的旋转矩阵和位置，并代入公式(14)和公式(15)求解。其中，位置信息为3*1的向量，方位角信息为3*3方阵。
[0185] (14)
[0186] (15)
[0187] 步骤10，将铰链体的位置和方位角改写为齐次矩阵，并进行渲染。
