随着深空探测技术的不断发展，利用行星车对目标星球进行探测成为获得该星球信息的一个重要途径。如美国已经发射的火星巡视器索杰纳号(1997年)、勇气号(2004年)、机遇号(2004年)及好奇号(2011年)都对火星表面进行了一系列探测及分析，中国成功着陆于月球表面的玉兔号月面巡视器(2013年)也对月球表面进行了一定的科学探测。上述部分行星车如图 1所示。

2020年中国即将把火星车送至火星表面，进行火星探测活动。为了更全面地得到行星表面的各种信息以及更好地适应行星表面复杂的地形地貌，行星车的结构设计变得越来越复杂，其自由度也随之增加，导致行星车的运动控制变得越来越困难，上述问题一直是国内外相关行业专家关注的重点也是难点。为解决行星车运动控制的难题，需要研究行星车的动力学建模与求解。

研究该问题时，通常将行星车看作由多刚体组成的机构。在此基础上，由于行星车具有差速结构，故在拓扑上具有闭链结构(见图 2)，该问题转化为闭链多刚体动力学建模与求解的问题。此外还需研究车轮与地面之间的轮土作用力，该力也是系统的外力。本文针对上述2个方面，介绍国内外研究者在近几十年内取得的主要研究成果，阐述相关算法的研究现状及基本思想。

1 行星车动力学建模

回溯动力学发展的历史，它起源于1686年Newton对于单质点的动力学分析。之后在1776年，Euler提出了刚体的概念，并结合Newton的质点动力学，给出了关于刚体运动的Newton-Euler方程^[1]，表示如下：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{f}} = m\mathit{\boldsymbol{\dot v}},\\ \tau = {\mathit{\boldsymbol{J}}_{\rm{C}}}\mathit{\boldsymbol{\dot \omega }} + \mathit{\boldsymbol{\omega }} \times {\mathit{\boldsymbol{J}}_{\rm{c}}} \cdot \mathit{\boldsymbol{\omega }}. \end{array} \right. $

(1)

其中：f和τ分别代表作用在刚体上的力和力矩，m和J_c分别代表刚体的质量和相对于质心的转动惯量，v和ω分别代表刚体的线速度和角速度，进而$\dot{\mathit{\pmb{ v}}}$和$\dot{\mathit{\pmb{ \omega }}}$分别代表刚体的加速度和角加速度。

而在1788年，Lagrange创立了分析动力学，根据能量守恒原理提出了刚体系统的Lagrange方程^[1]，表示如下：

$ \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial T}}{{\partial {{\dot q}_j}}}} \right) - \frac{{\partial T}}{{\partial {q_j}}} = {Q_j}\;\;\;\;j = 1,2, \cdots ,k. $

(2)

其中：q_j是第j个物体的广义坐标，T是系统动能，Q_j是对应于广义坐标q_j的广义力。在此之后，Hamilton提出了一种改进Lagrange方程的方法^[1]，即通过引入正则动量${p_j} = \partial T/\partial {\dot q_j}$, j=1, 2, …, k，用正则相空间(p₁, p₂, …, p_k；q₁, q₂, …q_k)，代替广义速率和广义坐标变量$\left( {{{\dot q}_1}, {{\dot q}_2}, \ldots {{\dot q}_k};{q_1}, {q_2}, \ldots {q_k}} \right)$。因此，Lagrange方程与式(3)等价。

$ {\dot p_j} = - \frac{{\partial H}}{{\partial {q_j}}},{\dot q_j} = - \frac{{\partial H}}{{\partial {p_j}}},\left( {j = 1,2, \cdots ,k} \right). $

(3)

其中H表示系统总能量。之后Boltzman-Hamel方程也被提出，丰富了分析动力学的内容。很多学者都在Newton-Euler法及Lagrange法的基础上持续研究刚体系统动力学的问题，并取得了一定的成果。

1.1 Kane法

虽然传统的Newton-Euler法及Lagrange法都可以解决刚体动力学建模问题，但在建立多刚体系统动力学模型时，他们都有各自的不足。1965年，Kane等^[2]提出了一种描述广义速度的方法，Singh等^[3]根据D'Alembert原理及文[2]，提出了对多刚体系统进行显式动力学建模的Kane方法，表示如下：

$ \begin{array}{*{20}{c}} {{m_j} \cdot \frac{{\partial \mathit{\boldsymbol{v}}_j^{\rm{T}}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}} \cdot {\mathit{\boldsymbol{a}}_j} + \frac{{\partial \mathit{\boldsymbol{\omega }}_j^{\rm{T}}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}} \cdot \left( {{\mathit{\boldsymbol{J}}_j} \cdot {\mathit{\boldsymbol{\alpha }}_j} + {{\mathit{\boldsymbol{\tilde \omega }}}_j} \cdot {\mathit{\boldsymbol{J}}_j} \cdot {\mathit{\boldsymbol{\omega }}_j}} \right) = }\\ {\frac{{\partial \mathit{\boldsymbol{v}}_j^{\rm{T}}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}} \cdot {\mathit{\boldsymbol{f}}_j} + \frac{{\partial \mathit{\boldsymbol{\omega }}_j^{\rm{T}}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}} \cdot {\mathit{\boldsymbol{\tau }}_j}.} \end{array} $

(4)

该方法结合了Newton-Euler方程与Lagrange方程的优点，并改进了二者的缺点。对二者的改进分别为：1)不考虑Newton-Euler方程中刚体之间的约束力及约束力矩；2)无需像Lagrange方程那样对能量方程求偏导数。

利用Kane方法，多体系统的动力学方程可以显式描述。文[4-5]给出了2种利用Kane方法进行动力学分析的应用。

虽然Kane方程可以给出多体系统动力学的显式方程，但在1965年，偏速度$\partial \mathit{\pmb{v}}_j^{\rm{T}}/\partial {\dot{\mathit{\pmb{ q}}}_j}$及$\partial \mathit{\pmb{\omega }}_j^{\rm{T}}/\partial {\dot{\mathit{\pmb{ q}}}_j}$很难求解，需要通过数值计算的方法对结果进行逼近，因此其计算效率也并未得到显著提升。1985年，Kane对该套方法思想进行了总结，消除了部分缺陷^[6]。

1.2 递归法

20世纪70年代末期至80年代初，随着计算机的发展，如何利用计算机解算多体系统动力学成为研究的焦点。因便于用计算机进行建模及计算，基于Newton-Euler方法及Lagrange法的递归算法受到众多学者关注。

在递归法中，多体系统的连接关系受到了大家的关注，如何用拓扑图反映系统内各零件的连接关系及相互运动关系成为焦点。Spanning树因其具有非环、可分叉、线性遍历等性质而受到学者们的关注。

1) Newton-Euler递归法。

文[7-9]先后提出了基于Newton-Euler的递归动力学方法，其动力学模型表示如下：

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{m_j}{\mathit{\boldsymbol{I}}_3}}&{ - {m_j} \cdot {{\mathit{\boldsymbol{\tilde r}}}_{{j_C}}}}\\ {{m_j} \cdot {{\mathit{\boldsymbol{\tilde r}}}_{{j_C}}}}&{{\mathit{\boldsymbol{J}}_{{j_C}}}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\dot v}}}_j}}\\ {{{\mathit{\boldsymbol{\dot \omega }}}_j}} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\omega }}_j} \times {\mathit{\boldsymbol{\omega }}_j} \times \left( {{m_{\rm{j}}} \cdot {{\mathit{\boldsymbol{\tilde r}}}_{{{\rm{j}}_{\rm{C}}}}}} \right)}\\ {{\mathit{\boldsymbol{\omega }}_j} \times {\mathit{\boldsymbol{J}}_{{{\rm{j}}_{\rm{C}}}}} \cdot {\mathit{\boldsymbol{\omega }}_j}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{}^j{\mathit{\boldsymbol{f}}_j}}\\ {{}^j{\mathit{\boldsymbol{\tau }}_j}} \end{array}} \right].} \end{array} $

(5)

其中：${\tilde {\mathit{\pmb{r}}}_{jc}}$为体j质心位置的叉乘矩阵，I₃为单位矩阵。在此基础上，不同学者针对传统Newton-Euler递归法进行了改进，如表 1所示。

文[10-14]通过一系列文章提出了一种基于Newton-Euler递归法的求解多体动力学的理论-空间算子代数法。该方法的核心为给定空间核心算子ξ_ϕ(k, m)及空间传播算子ϕ(k, m)：

$ {\xi _\phi }\left( {k,m} \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_3}}&{\mathit{\boldsymbol{\tilde r}}\left( {k,m} \right)}\\ {\bf{0}}&{{\mathit{\boldsymbol{I}}_3}} \end{array}} \right]. $

(6)

其中$\tilde {\mathit{\pmb{r}}}$(k, m)表示关节k及关节m之间的距离。该系统给出如下设定：假定系统固接处编号为N，从该固接处依次对关节进行编号直至末端。关节k与杆件(k+1)及杆件k相连，其中杆件k与k相连的部分定义为p^－(k)，而杆件(k+1)与k相连的部分定义为p⁺(k)，p为任意属性符号，包括速度，加速度等。根据空间核心算子，得到式(1)的迭代表示：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{V}}\left( k \right) = \xi _\phi ^{\rm{T}}\left( {k + 1,k} \right)\mathit{\boldsymbol{V}}\left( {k + 1} \right) + \\ {\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right)\dot \theta \left( k \right),\\ \mathit{\boldsymbol{\alpha }}\left( k \right) = \xi _\phi ^{\rm{T}}\left( {k + 1,k} \right)\alpha \left( {k + 1} \right) + \\ {\mathit{\boldsymbol{H}}^{\rm{T}}}\left( k \right)\ddot \theta \left( k \right) + \mathit{\boldsymbol{a}}\left( k \right),\\ \mathit{\boldsymbol{f}}\left( k \right) - {\xi _\phi }\left( {k,k - 1} \right)\mathit{\boldsymbol{f}}\left( {k - 1} \right) = \\ \mathit{\boldsymbol{M}}\left( k \right)\mathit{\boldsymbol{\alpha }}\left( k \right) + \mathit{\boldsymbol{b}}\left( k \right),\\ \mathit{\boldsymbol{\tau }}\left( k \right) = \mathit{\boldsymbol{H}}\left( k \right)\mathit{\boldsymbol{f}}\left( k \right). \end{array} \right. $

(7)

其中：V(k)表示杆件k的空间广义速度，$\mathit{\pmb{\alpha}} \left( k \right) = {\dot{\mathit{\pmb{ V}}}}\left( k \right)$表示杆件k的空间广义加速度，b(k)表示与速度相关的回转力。假设ϕ(k, m)=(I₆－ξ_ϕ(k, m))^－1，则式(7)可改写为

$ \left\{ \begin{array}{l} V = {\phi ^{\rm{T}}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\dot \theta ,\\ \mathit{\boldsymbol{\alpha }} = {\phi ^{\rm{T}}}\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\ddot \theta + a} \right),\\ \mathit{\boldsymbol{f}} = \phi \left( {\mathit{\boldsymbol{M}}\alpha + b} \right),\\ \mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{Hf}}. \end{array} \right. $

(8)

其中ϕ^T(k+1, k)可以看作表示关节两端传递的Jacobi矩阵。该算子为线性算子，利用它可以实现广义力、广义速度及广义加速度的反向或正向传播。利用空间算子代数，系统可从高层框架中描述对机构的控制及轨迹跟踪，文[13]给出了具有多叉树拓扑结构的多体系统的动力学方程。该方程为：

$ \left\{ \begin{gathered} \mathit{\boldsymbol{\tau }} = \mathbb{M}\left( \theta \right)\ddot \theta + \mathit{\boldsymbol{C}}\left( {\theta ,\dot \theta } \right), \hfill \\ \mathbb{M}\left( \theta \right) = \mathit{\boldsymbol{H}}\phi \mathit{\boldsymbol{M}}{\phi ^{\text{T}}}{\mathit{\boldsymbol{H}}^{\text{T}}} \in {{\bf{R}}^{N \times N}}, \hfill \\ \mathit{\boldsymbol{C}}\left( {\theta ,\dot \theta } \right) = \mathit{\boldsymbol{H}}\phi \left( {\mathit{\boldsymbol{M}}{\phi ^{\text{T}}}\mathit{\boldsymbol{a}} + \mathit{\boldsymbol{b}}} \right) \in {{\bf{R}}^N}. \hfill \\ \end{gathered} \right. $

(9)

文[13-14]还给出了闭链系统的解决方案，其方法为选定一个关节，将其从关节处切开，则该关节内部的相互作用力转变为系统两边各增添一个外力，再对系统进行分析。

文[15-20]提出了利用类似于Plücker符号的6-D空间矢量对Euler方程进行了改进，其模型为：

$ {\mathit{\boldsymbol{f}}_i} = {\mathit{\boldsymbol{I}}_i}{\mathit{\boldsymbol{a}}_i} + {\mathit{\boldsymbol{\tilde v}}_i} \cdot {\mathit{\boldsymbol{I}}_i}{\mathit{\boldsymbol{v}}_i}. $

(10)

其中: f_i为作用在体i上的广义空间作用力，v_i和a_i分别为体i的广义速度和广义加速度，I_i为体i的惯性张量。它们都是6-D空间矢量完整地讲述了6-D空间代数的计算规则、运动学的递归公式、动力学的递归公式及闭链系统的解决方法。文[20]给出了便于计算机实现的算法。对于闭链问题，其思路与Jain^[13]提到的思路一致，都是将闭链拆开，将物体之间的相互作用力当作作用于系统的外力，重新对系统动力学进行分析，如图 3所示。

在最近的多体动力学建模理论中，文[21-22]针对式(11)中的闭链系统，提出了一种基于Newton-Euler方程的解决方法。

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{M\ddot q}} + \mathit{\Phi }_q^{\rm{T}}\mathit{\boldsymbol{\lambda }} = \mathit{\boldsymbol{f}},\\ \mathit{\Phi }\left( \mathit{\boldsymbol{q}} \right) = 0. \end{array} \right. $

(11)

其中: λ为Lagrange乘子矢量；Φ为运动学约束方程，表示闭链中相连物体的约束。利用嵌入技术对Lagrange乘子进行消元，得到新的动力学方程如下：

$ \mathit{\boldsymbol{M\ddot q}} + \mathit{\Phi }_q^{\rm{T}}\mathit{\boldsymbol{\rho }}\left( {\mathit{\ddot \Phi } + 2\mathit{\boldsymbol{\omega \zeta }}\mathit{\dot \Phi } + {\mathit{\boldsymbol{\omega }}^2}\mathit{\Phi }} \right) = \mathit{\boldsymbol{f}}. $

(12)

其加速度迭代式表示为

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\ddot q}}}_0} = {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{f}};\\ {{\mathit{\boldsymbol{\ddot q}}}_k} = {\left( {\mathit{\boldsymbol{M}} + \mathit{\Phi }_q^{\rm{T}}\mathit{\boldsymbol{\rho }}{\mathit{\Phi }_q}} \right)^{ - 1}} \cdot \\ \left\{ {\mathit{\boldsymbol{M}}{{\mathit{\boldsymbol{\ddot q}}}_{k - 1}} - \mathit{\Phi }_q^{\rm{T}}\mathit{\boldsymbol{\rho }} \cdot } \right.\\ \left. {\left( {{{\mathit{\dot \Phi }}_q}\mathit{\boldsymbol{\dot q}} + {{\mathit{\dot \Phi }}_t} + 2\mathit{\boldsymbol{\omega }}\zeta \mathit{\dot \Phi } + {\mathit{\boldsymbol{\omega }}^2}\mathit{\Phi }} \right)} \right\},\;\;\;\;\;k \ge 1. \end{array} \right. $

(13)

其中：ρ、ω、ζ分别为惩罚系数矩阵、自然频率矩阵、阻尼率矩阵。因Φ与机构参数相关，故可以利用Gröbner基对Φ进行三角化，避免了对Φ进行迭代求解，保证了系统的实时应用。该方法需要先利用Maple对公式符号进行简化，再利用C++进行编程，因此每次更换系统时，都需要首先使用Maple化简公式，再利用其他软件进行计算。因此，该方法的实用性还需要继续提高。

Newton-Eular法的建模复杂度与关节个数呈线性关系，该特点保证其建模效率很高。但对于每一个刚体，都需要用6个方程描述刚体的力或力矩平衡。对于用计算机实现一个多体系统的动力学解算，使用该方法需要很强的计算机编程功底，也需要对Newton-Euler动力学非常了解。

2) Lagrange递归法。

对于Lagrange递归动力学，Hollerbach^[23]提出了基于Lagrange方程的动力学建模方法，该方法将Newton-Euler法中位于质心的体坐标系变为位于关节上的D-H坐标系，通过沿轴的运动表示刚体运动，得到一个与关节个数呈线性关系的递归动力学方程，该方法在时间复杂度上与Newton-Euler递归动力学方法相同，但在每一步的递归操作中，要执行更多的加法与乘法操作。Hollerbach提出的方法具有运动学的正向递归与力的反向递归，同时他利用3×3维的旋转变换阵代替4×4维的齐次变换阵^[24-25]，降低了每一次递归的计算量。

Silver^[26]选择了用广义速度及广义加速度矢量作为方程的表示，改进了Hollerbach给出的方程，并证明了在新的方程中，Lagrange递归方法与Newton-Euler递归方法具有相同的时间复杂度。方程表示如下：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\tau }}_j} = \frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\frac{{\partial T}}{{\partial {{\dot q}_j}}}} \right) - \frac{{\partial T}}{{\partial {q_j}}} = }\\ {\sum\limits_i {\left( {{m_i}\mathit{\boldsymbol{\dot v}}_i^{\rm{T}} \cdot \frac{{\partial {\mathit{\boldsymbol{v}}_i}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}} + {{\left[ {{{\mathit{\boldsymbol{\dot \omega }}}_i} \cdot {\mathit{\boldsymbol{J}}_i} + {{\mathit{\boldsymbol{\tilde \omega }}}_i} \cdot {\mathit{\boldsymbol{J}}_i} \cdot {\mathit{\boldsymbol{\omega }}_i}} \right]}^{\rm{T}}} \cdot \frac{{\partial {\mathit{\boldsymbol{\omega }}_i}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}}} \right)} ;} \end{array} $

(14)

$ \begin{array}{l} \frac{{\partial {\mathit{\boldsymbol{v}}_i}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}} = \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\tilde z}}}_j} \cdot {\mathit{\boldsymbol{r}}_{j,i}},\;\;\;\;\;j \in {\bf{R}};\\ {{\mathit{\boldsymbol{\tilde z}}}_j},\;\;\;\;\;\;\;\;\;\;j \in {\bf{P}}; \end{array} \right.\\ \frac{{\partial {\mathit{\boldsymbol{\omega }}_i}}}{{\partial {{\mathit{\boldsymbol{\dot q}}}_j}}} = \left\{ \begin{array}{l} {\mathit{\boldsymbol{z}}_j},\;\;\;\;\;j \in {\bf{R}};\\ 0,\;\;\;\;\;\;j \in {\bf{P}}. \end{array} \right. \end{array} $

(15)

其中：z_j表示第j轴的方向矢量，r_{j, i}表示关节j与关节i的距离矢量。

居鹤华等^[27]对关节进行抽象，得到了轴不变量的概念，轴不变量即为关节的运动轴。根据轴不变量，建立以轴系为参考，以Spanning树为拓扑结构的轴链系统。在此基础上根据Lagrange方程推导出Lagrange递归平衡方程，并将其归纳为规范型，表示如下：

$ \left\{ \begin{array}{l} {}^{i\left| {\bar u} \right.}\mathit{\boldsymbol{n}}_u^{\rm{T}} \cdot \mathit{\boldsymbol{M}}_P^{\left[ u \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q + }}{}^{i\left| {\bar u} \right.}\mathit{\boldsymbol{n}}_u^{\rm{T}} \cdot \mathit{\boldsymbol{h}}_P^{\left[ u \right]} = f_u^D,\\ \;\;\;\;\;{}^{\bar u}{\mathit{\boldsymbol{k}}_u} \in \mathit{\boldsymbol{P}};\\ {}^{i\left| {\bar u} \right.}\mathit{\boldsymbol{n}}_u^{\rm{T}} \cdot \mathit{\boldsymbol{M}}_R^{\left[ u \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q + }}{}^{i\left| {\bar u} \right.}\mathit{\boldsymbol{n}}_u^{\rm{T}} \cdot \mathit{\boldsymbol{h}}_R^{\left[ u \right]} = \tau _u^D,\\ \;\;\;\;\;{}^{\bar u}{\mathit{\boldsymbol{k}}_u} \in \mathit{\boldsymbol{R}}. \end{array} \right. $

(16)

其中：M_P^[u][·]和M_R^[u][·]分别表示平动副和转动副的惯性矩阵，$^{i\left| {\bar u} \right.}{\mathit{\pmb{n}}_u}$表示轴不变量，f_u^D和τ_u^D分别表示系统投影到轴u的广义力和广义力矩。式(16)中系统两侧在轴上保持广义力平衡。此外，代表矩阵M_P及M_R的第u行的表达式分别为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_P^{\left[ u \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q}} = \sum\limits_l^{{i_{{l_{\bar u}}}}} {\left( {\begin{array}{*{20}{c}} {\sum\limits_j^{{u_L}} {\left( {{m_j}} \right) \cdot 1} }\\ { - \sum\limits_j^{{u_L}} {\left( {{m_j} \cdot {}^{i\left| l \right.}{\mathit{\boldsymbol{r}}_{jI}}} \right)} } \end{array}\left\| { \cdot {}^{i\left| {\overline l } \right.}{\mathit{\boldsymbol{n}}_l} \cdot \mathit{\boldsymbol{\ddot q}}_l^{\overline l }} \right.} \right)} + }\\ {\sum\limits_k^{{u_L}} {\left( {\begin{array}{*{20}{c}} {\sum\limits_j^{{k_L}} {\left( {{m_j}} \right) \cdot 1} }\\ { - \sum\limits_j^{{k_L}} {\left( {{m_j} \cdot {}^{i\left| k \right.}{\mathit{\boldsymbol{r}}_{jI}}} \right)} } \end{array}\left\| { \cdot {}^{i\left| {\overline k } \right.}{\mathit{\boldsymbol{n}}_k} \cdot \mathit{\boldsymbol{\ddot q}}_k^{\overline k }} \right.} \right)} ,} \end{array} $

(17)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_R^{\left[ u \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q}} = }\\ {\sum\limits_l^{{i_{{l_{\bar u}}}}} {\left( {\begin{array}{*{20}{c}} {\sum\limits_j^{{u_L}} {\left( {{m_j} \cdot {}^{i\left| u \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}}} \right)} }\\ {\sum\limits_j^{{u_L}} {\left( \begin{array}{l} {}^{i\left| {jI} \right.}{\mathit{\boldsymbol{J}}_{jl}} - \\ {m_j} \cdot {}^{i\left| u \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}} \cdot {}^{i\left| l \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}} \end{array} \right)} } \end{array}\left\| { \cdot {}^{i\left| {\overline l } \right.}{\mathit{\boldsymbol{n}}_l} \cdot \mathit{\boldsymbol{\ddot q}}_l^{\overline l }} \right.} \right)} + }\\ {\sum\limits_k^{{u_L}} {\left( {\begin{array}{*{20}{c}} {\sum\limits_j^{{k_L}} {\left( {{m_j} \cdot {}^{i\left| u \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}}} \right)} }\\ {\sum\limits_j^{{k_L}} {\left( \begin{array}{l} {}^{i\left| {jI} \right.}{\mathit{\boldsymbol{J}}_{jl}} - \\ {m_j} \cdot {}^{i\left| u \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}} \cdot {}^{i\left| k \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}} \end{array} \right)} } \end{array}\left\| { \cdot {}^{i\left| {\overline k } \right.}{\mathit{\boldsymbol{n}}_k} \cdot \mathit{\boldsymbol{\ddot q}}_k^{\overline k }} \right.} \right)} ,} \end{array} $

(18)

其中：$^i{l_{\bar u}}$表示运动链i节点到u的父节点，^uL表示u节点的闭子树。作用于节点u的Coriolis力及重力的表示h_P^[u]及h_R^[u]分别为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{h}}_P^{\left[ u \right]} = \sum\limits_k^{{u_L}} {\left( {{m_k} \cdot \sum\limits_l^{{i_{{l_{kI}}}}} {\left( {{}^i\tilde \phi _{\overline l }^{ \vdots 2} \cdot {}^{i\left| {\bar l} \right.}{r_l}} \right)} + } \right.} }\\ {\left. {\left( {2 \cdot {}^i{{\tilde \phi }_{\overline l }} \cdot {}^{i\left| {\bar l} \right.}{{\dot r}_l}} \right)} \right) - \sum\limits_k^{{u_L}} {\left( {{m_k} \cdot {}^i{\mathit{\boldsymbol{g}}_{kl}}} \right)} ,} \end{array} $

(19)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{h}}_R^{\left[ u \right]} = \sum\limits_k^{{u_L}} {\left( {{m_k} \cdot {}^{i\left| u \right.}{{\mathit{\boldsymbol{\tilde r}}}_{kl}} \cdot \sum\limits_l^{{i_{{l_{kI}}}}} {\left( {{}^i\tilde \phi _{\overline l }^{ \vdots 2} \cdot {}^{i\left| {\bar l} \right.}{r_l}} \right)} + } \right.} }\\ {\left. {\left( {2 \cdot {}^i{{\tilde \phi }_{\overline l }} \cdot {}^{i\left| {\bar l} \right.}{{\mathit{\boldsymbol{\dot r}}}_l}} \right)} \right) + \sum\limits_k^{{u_L}} {\left( {{}^i{{\tilde \phi }_k} \cdot {}^{i\left| {kI} \right.}{\mathit{\boldsymbol{J}}_{kI}} \cdot {}^i{\phi _k}} \right)} - }\\ {\sum\limits_k^{{u_L}} {\left( {{m_k} \cdot {}^{i\left| u \right.}{{\mathit{\boldsymbol{\tilde r}}}_{kI}} \cdot {}^i{\mathit{\boldsymbol{g}}_{kI}}} \right)} .} \end{array} $

(20)

作用于节点u的广义合外力f_u^D及τ_u^D表示为

$ \left\{ \begin{array}{l} f_u^\mathit{\boldsymbol{D}} = f_u^c + \dot r_u^{\bar u} \cdot {\rm{G}}\left( {f_u^c} \right) + \\ \;\;\;\;\;\;{}^{i\left| {\bar u} \right.}\mathit{\boldsymbol{n}}_u^{\rm{T}} \cdot \sum\limits_l^{{u_L}} {\left( {{}^{iS}{\mathit{\boldsymbol{f}}_{lS}}} \right)} ,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{}^{u - }{\mathit{\boldsymbol{k}}_u} \in {\bf{P}};\\ \tau _u^\mathit{\boldsymbol{D}} = \tau _u^c + \dot \phi _u^{\bar u} \cdot {\rm{G}}\left( {\tau _u^c} \right) + {}^{i\left| {\bar u} \right.}\mathit{\boldsymbol{n}}_u^{\rm{T}} \cdot \\ \;\;\;\;\;\;\;\;\sum\limits_l^{{u_L}} {\left( {{}^{i\left| {\bar u} \right.}{{\mathit{\boldsymbol{\tilde r}}}_{lS}} \cdot {}^{iS}{\mathit{\boldsymbol{f}}_{lS}} + {}^i{\mathit{\boldsymbol{\tau }}_l}} \right)} ,\;\;\;\;\;\;\;\;\;{}^{u - }{\mathit{\boldsymbol{k}}_u} \in {\bf{R}}. \end{array} \right. $

(21)

其中$f_u^c + \dot r_u^{\bar u}\cdot{\rm{G}}\left( {f_u^c} \right)$和$\tau _u^c + \dot \phi _u^{\bar u}\cdot{\rm{G}}\left( {\tau _u^c} \right)$表示广义共轴驱动力即电机驱动力和力矩。$^{i\left| {\bar u} \right.}\mathit{\pmb{n}}_u^{\rm{T}}\cdot\sum\limits_l^{^uL} {\left( {^{iS}{\mathit{\pmb{f}}_{lS}}} \right)} $及$^{i\left| {\bar u} \right.}\mathit{\pmb{n}}_u^{\rm{T}}\cdot\sum\limits_l^{^uL} {\left( {^{i|u}{{\tilde {\mathit{\pmb{r}}}}_{lS}}{\cdot^{iS}}{\mathit{\pmb{f}}_{lS}}{ + ^i}{\mathit{\pmb{\tau}} _l}} \right)} $表示外力在轴u的投影。

在式(17)—(20)中，

$ \left. {{}_b^a} \right\| \cdot {}^{\bar k}{\mathit{\boldsymbol{n}}_k} \cdot q_k^{\bar k} \buildrel \Delta \over = \left\{ {\begin{array}{*{20}{c}} {a \cdot {}^{\bar k}{\mathit{\boldsymbol{n}}_k} \cdot q_k^{\bar k},}&{{}^{\bar k}{\mathit{\boldsymbol{k}}_k} \in {\bf{P}};}\\ {b \cdot {}^{\bar k}{\mathit{\boldsymbol{n}}_k} \cdot q_k^{\bar k},}&{{}^{\bar k}{\mathit{\boldsymbol{k}}_k} \in {\bf{R}}.} \end{array}} \right. $

(22)

进一步，针对星球车动基座的特点，对车体进行动力学分析，因车体有6个自由度，无须将广义力投影到某个轴上，直接得到基座的动力学模型，表示如下：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{M}}_P^{\left[ c \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q}} + \mathit{\boldsymbol{h}}_P^{\left[ c \right]} = {}^{i\left| \mathit{\boldsymbol{D}} \right.}{\mathit{\boldsymbol{f}}_c},\\ \mathit{\boldsymbol{M}}_R^{\left[ c \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q}} + \mathit{\boldsymbol{h}}_R^{\left[ c \right]} = {}^{i\left| \mathit{\boldsymbol{D}} \right.}{\mathit{\boldsymbol{\tau }}_c}. \end{array} \right. $

(23)

其中

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_P^{\left[ c \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q}} = \sum\limits_k^{{c_L}} {\left( {{m_k}} \right) \cdot {}^i{\mathit{\boldsymbol{Q}}_c} \cdot \mathit{\boldsymbol{\ddot r}}_{\left( {i,c} \right]}^{\rm{T}}} - }\\ {\sum\limits_k^{{c_L}} {\left( {{m_k} \cdot {}^{i\left| c \right.}{{\mathit{\boldsymbol{\tilde r}}}_{kI}}} \right) \cdot {}^i{\mathit{\Theta }_c} \cdot \ddot \phi _{\left( {i,c} \right]}^{\rm{T}}} + }\\ {\sum\limits_k^{{c_{ - L}}} {\left( {\left. {\begin{array}{*{20}{c}} {\sum\limits_j^{{k_L}} {\left( {{m_j}} \right)} }\\ { - \sum\limits_j^{{k_L}} {\left( {{m_j} \cdot {}^{i\left| k \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jl}}} \right)} } \end{array}} \right\| \cdot {}^{i\left| {\bar k} \right.}{\mathit{\boldsymbol{n}}_k} \cdot \ddot q_k^{\bar k}} \right)} ,} \end{array} $

(24)

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}_P^{\left[ c \right]\left[ \cdot \right]} \cdot \mathit{\boldsymbol{\ddot q}} = \sum\limits_k^{{c_L}} {\left( {{}^{i\left| {kI} \right.}{\mathit{\boldsymbol{J}}_{kl}} - {m_k} \cdot {}^{i\left| c \right.}\mathit{\boldsymbol{\tilde r}}_{kl}^{ \vdots 2}} \right)} \cdot {}^i{\mathit{\Theta }_c} \cdot }\\ {\ddot \phi _{\left( {i,c} \right]}^{\rm{T}} + \sum\limits_k^{{c_L}} {\left( {{m_k} \cdot {}^{i\left| c \right.}{{\mathit{\boldsymbol{\tilde r}}}_{kI}}} \right) \cdot {}^i{\mathit{\boldsymbol{Q}}_c} \cdot \mathit{\boldsymbol{\ddot r}}_{\left( {i,c} \right]}^{\rm{T}} + } }\\ {\sum\limits_k^{{c_{ - L}}} {\left( {\left. {\begin{array}{*{20}{c}} {\sum\limits_j^{{k_L}} {\left( {{m_j} \cdot {}^{i\left| c \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}}} \right)} }\\ {\sum\limits_j^{{k_L}} {\left( \begin{array}{l} {}^{i\left| {jI} \right.}{\mathit{\boldsymbol{J}}_{jI}} - \\ {m_j} \cdot {}^{i\left| c \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}} \cdot {}^{i\left| k \right.}{{\mathit{\boldsymbol{\tilde r}}}_{jI}} \end{array} \right)} } \end{array}} \right\| \cdot {}^{i\left| {\bar k} \right.}{\mathit{\boldsymbol{n}}_k} \cdot \ddot q_k^{\bar k}} \right)} } \end{array} $

(25)

而h_P^[u]及h_R^[u]与式(19)和(20)一致。

对于Lagrange迭代法，其正动力学解为

$ \mathit{\boldsymbol{\ddot q}} = {\mathit{\boldsymbol{M}}^{ - 1}} \cdot \left( {\mathit{\boldsymbol{f}} - \mathit{\boldsymbol{h}}} \right). $

(26)

根据上文描述可知，传统Lagrange递归法定义清晰明确，比Newton-Euler法具有更简单的方程描述。对于非动力学的研究人员，利用编程实现Lagrange递归动力学比Newton-Euler递归动力学更加容易。而基于轴不变量的Lagrange递归法将广义力投影到运动轴上，其编程过程更加简单，且计算效率略高于传统Lagrange递归法。

根据文[13, 20, 24, 27]，对各递归方法的建模复杂度及正动力学计算复杂度进行总结，如表 2所示。

由表 2可知，递归法的建模复杂度都为O(n)，因惯性矩阵求逆的计算复杂度为O(n²)，故正动力学计算复杂度不会优于O(n²)。

1.3 性能比较

行星车是一种几乎所有关节都是转动副的多刚体系统，只有车身具有平移能力，文[27]将其称为动基座的多体系统。当前有多种开源软件库或商业软件可以对多体动力学进行建模。ODE、Bullet、SIMBODY、RBDL等都是著名的开源软件。相关的商业软件有ADAMS、RecurDyn、SIMPACK及SAMCEF等。许多实验室也针对各自的课题，开发出相应的动力学仿真软件包。在动力学建模与求解中，每个方法都有自己的特点，各模型之间并没有绝对的好坏之分，很难区分应用不同方法建模后模型精度的高低，一般通过建模速度及求解速度区分算法的优劣。但对于不同的建模方法，在工程运用中受到编程水平及所应用不同的编程语言导致效率不同的约束，其求解时间存在差异，因此极少能在文献中查到对不同建模方法效率的评价。对于上述已知的软件，绝大多数都应用Newton-Euler迭代法对多体系统进行建模，而Lagrange递归法多用于解决机械臂柔性控制的问题^[28-30]。

本实验室利用基于轴不变量的Lagrange递归法通过C++语言对22自由度的月球车进行了动力学建模及解算。根据执行100 000次正动力学循环的求解时间，得到单步的平均求解时间为17.95 ms。仿真软件的示意图如图 4所示。

在文[31]中，RBDL软件的开发者Felis对比了利用C语言开发的RBDL与SIMBODY软件包的仿真结果。RBDL软件中使用的算法主要为文[20]提到的6-D空间矢量法。SIMBODY主要采用文[13]中提到的空间算子操作代数算法对多刚体进行动力学建模。Felis使用2款软件分别求解一个35自由度的多刚体系统的正动力学，平均求解时间表明RBDL软件包具有更快的计算效率。其具体执行结果如表 3所示。

此外，作者还分别对比了针对5~100个自由度的系统，2款软件执行100 000次正动力学函数所消耗的平均时间，结果如图 5所示。

由表 3和图 5可知，RBDL比SIMBODY软件的求解效率更高，是一个优秀的第三方软件库。当今最著名的机器人仿真平台ROS系统就使用了RBDL作为其动力学的支持库。

2 轮土力学

轮土力学是研究车轮与土壤之间作用关系的学科。对于行星车探测，因行星车车轮为钢轮结构，而行星表面土壤多为松软泥土，因此轮土力学中普遍使用钢轮软土模型，如图 6所示，其中p₁~p₆为轮地接触点，f₁~f₆为轮地接触力，f_s为车体牵引力。

Bekker^[33-34]轮土力学是该领域内最初的成果，该模型是针对弹塑性的土壤及刚性驱动轮的作用力模型。

$ \sigma \left( z \right) = \left( {\frac{{{k_1}}}{b} + {k_2}} \right) \cdot {z^n}. $

(27)

Bekker轮土力学模型是基于水平地形及重力方向载荷建立的。但上述方法存在2个缺陷：1)是轮式多轴系统驱动轮接触的地形通常是起伏变化的，相应地载荷随地形起伏而变化；2) Bekker轮土力学模型未给出驱动轮与土壤作用过程中相关物理量的参考系定义。因此，需要以Bekker轮土力学模型为基础，增加相应的指标系统，提高轮土力学模型表达的准确性。

Janosi等^[35]提出了土壤与驱动轮的剪应力模型；Wong等^[36-37]继续研究土壤的剪切变形，最终将剪应力表达式化简为

$ \tau = \left( {c + \sigma \tan \left( \phi \right)} \right)\left( {1 - {\rm{e}}\left( { - j/k} \right)} \right). $

(28)

式(27)和(28)中的轮土基本参数如表 4所示。

文[32]及[38-42]不断对式(27)和(28)进行改进，最后得到了一种轮土几何关系, 如图 7所示。图中补充的轮土参数如表 5所示。

图 7中，有

$ z\left( \theta \right) = r\left( {\cos \theta - \cos {\theta _1}} \right). $

(29)

将式(29)代入式(27)—(28)中，得到：

$ {\sigma _1}\left( \theta \right) = \left( {\frac{{{k_1}}}{b} + {k_2}} \right) \cdot {\left( {r\left( {\cos \theta - \cos {\theta _1}} \right)} \right)^n}, $

(30)

$ \begin{array}{*{20}{c}} {{\sigma _2}\left( \theta \right) = \left( {\frac{{{k_1}}}{b} + {k_2}} \right) \cdot }\\ {{{\left( {r\left( {\cos \left( {{\theta _1} - \theta \cdot \frac{{{\theta _1} - {\theta _{\rm{m}}}}}{{{\theta _{\rm{m}}}}}} \right) - \cos {\theta _1}} \right)} \right)}^n}.} \end{array} $

(31)

此外，法向最大应力角θ_m表示为

$ {\theta _{\rm{m}}} = \left( {{c_1} + {c_2}s} \right) \cdot {\theta _1}. $

(32)

其中s为轮地滑移率。根据θ_m可以得到对应的σ_m及τ_m，进一步得到修正后的正应力及剪应力表达式：

$ \begin{array}{l} {\sigma _1}\left( \theta \right) = \frac{{{\theta _1} - \theta }}{{{\theta _1} - {\theta _{\rm{m}}}}}{\sigma _{\rm{m}}},\;\;\;\;\;\;{\sigma _2}\left( \theta \right) = \frac{\theta }{{{\theta _{\rm{m}}}}}{\sigma _{\rm{m}}},\\ {\tau _1}\left( \theta \right) = \frac{{{\theta _1} - \theta }}{{{\theta _1} - {\theta _{\rm{m}}}}}{\tau _{\rm{m}}},\;\;\;\;\;\;\;\;{\tau _2}\left( \theta \right) = c + \frac{\theta }{{{\theta _{\rm{m}}}}}\left( {{\tau _{\rm{m}}} - c} \right). \end{array} $

(33)

进一步根据图 8所示的车轮坐标系，分析车轮所受的力与力矩，根据滑移角α得到轮地法向力、牵引力、牵引力矩及侧向力。

法向力表示为

$ \begin{array}{*{20}{c}} {N = \frac{{rb}}{{{\theta _{\rm{m}}}\left( {{\theta _{\rm{l}}} - {\theta _{\rm{m}}}} \right)}} \cdot }\\ {\left[ {\begin{array}{*{20}{c}} {{\sigma _{\rm{m}}}\left( { - {\theta _{\rm{m}}}\cos \left( {{\theta _{\rm{1}}}} \right) + {\theta _{\rm{1}}}\cos \left( {{\theta _{\rm{m}}}} \right) - {\theta _{\rm{1}}} + {\theta _{\rm{m}}}} \right)}\\ { - {\kappa _{\rm{m}}}\left( {{\theta _{\rm{m}}}\sin \left( {{\theta _1}} \right) - {\theta _1}\sin \left( {{\theta _{\rm{m}}}} \right)} \right)}\\ { - c\left( { - {\theta _{\rm{1}}}\sin \left( {{\theta _{\rm{m}}}} \right) - {\theta _{\rm{m}}}\sin \left( {{\theta _{\rm{m}}}} \right) - {\theta _{\rm{m}}}{\theta _{\rm{1}}} + \theta _{\rm{m}}^2} \right)} \end{array}} \right].} \end{array} $

(34)

牵引力及牵引力矩表示为

$ \left\{ \begin{array}{l} T = \frac{{rb}}{2}\left( {{\tau _{\rm{m}}}{\theta _1} + c{\theta _{\rm{m}}}} \right),\\ \tau = \frac{{{r^2}b}}{2}\left( {{\tau _{\rm{m}}}{\theta _1} + c{\theta _{\rm{m}}}} \right). \end{array} \right. $

(35)

侧向力表示为

$ L = \left( {cbA + N \cdot \tan \left( \phi \right)} \right) \cdot \left( {1 - {\rm{e}}\left( { - B \cdot \alpha } \right)} \right). $

(36)

其中：A为轮地接触面积，B为常系数。

根据Coulomb摩擦因数μ_x，可得轮土作用的前向Coulomb摩擦力f_x表达式：

$ {f_x} = {\mu _x} \cdot N \cdot {\rm{sign}}\left( {{v_x}} \right). $

(37)

利用上述模型，研究人员分别对比模拟仿真实验及在真实环境的实验，证明了式(34)—(36)轮土力学模型的可行性。

在文[40]中给出了结合轮土力学的行星车动力学建模的方法并进行了相应的仿真实验。仿真实验与实物实验的结果对比表明：无论是平地行驶还是爬坡实验，仿真实验的曲线都能反映实际的车体的运行基本状态。平地实验结果如图 9所示，其中左侧为实物实验数据，右侧为仿真数据。

文[43]中比较了2组利用轮土力学进行的动力学仿真实验与实物实验的位置误差与方向误差的计算结果(见图 10)。实验条件为：左右前轮都分别固定为15°及30°，左右后轮固定为0°。结果表明，在模拟行星车轮与土壤的相互作用时，利用上述提到的一系列方法，可以高精度对轮土力学进行仿真。结果如表 6所示。

文[44-45]介绍了在中国嫦娥三号辅助实验中应用Newton-Euler递归动力学模型及轮土力学模型实现嫦娥三号实时动力学控制的成功经验。其中地面软土模型通过点云数据实现，通过布置在轮轴平面上的虚拟传感器向下发射虚拟射线，可以得到射线与点云数据的交点。从而得到钢轮的泥土进入角和泥土退出角，再进一步计算法向力、牵引力及侧向力，如图 11所示。

虽然该轮土力学模型能够正确地计算钢轮软土的相互作用力，但使用该方法首先需要确定土壤参数。若土壤摩擦参数、黏滞参数未知，则计算结果将存在一定的误差。

3 结论

近些年来，对于多刚体动力学建模的研究逐渐减少，许多的研究者都将目标转向了柔性体的动力学建模。而对于轮土力学，研究人员仍按照文中提到的成熟算法进行仿真或实验。在有关行星车车体的动力学仿真实验中，不同学者利用Newton-Euler递归法或Lagrange递归法进行了各种各样的实验，都取得了不错的结果，二者之间很难明确地划分性能的优劣。由于Kane方法不适于用计算机完成多体系统的动力学建模，故并未见到许多利用该方法完成实验的成果。