火星巡视器动力学建模研究与仿真平台实现

徐颢，居鹤华，余萌，顾朋; XU Hao; JU Hehua; YU Meng; GU Peng

2025年7月29日 1:04 星期二

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

火星巡视器动力学建模研究与仿真平台实现 PDF

- ORCID：
徐颢 ¹
- ORCID：
居鹤华 ¹
- ORCID：
余萌 ¹
- ORCID：
顾朋 ²

1. 南京航空航天大学航天学院，南京210016； 2. 北京控制工程研究所，北京100190

中图分类号： V476.4

最近更新：2021-03-23

DOI：10.16356/j.1005-2615.2021.01.002

摘要

火星巡视器是执行复杂的火星表面探测任务的重要载体，为保证巡视器的可靠性，本文对火星巡视器动力学原理进行了研究。首先，分析了巡视器的结构组成，介绍了巡视器动力学建模的相关理论基础。然后根据巡视器有向Span树的轴链系统，基于拉格朗日递归动力学分析方法，得到了多轴系统拉格朗日方程，通过遍历巡视器轴链拓扑完成了巡视器动力学模型的建立。接着，将链符号系统引入经典Bekker轮土力学模型，以提高其准确性。最后，在VS平台MFC框架中开发了动力学仿真软件，对本文的动力学系统进行了测试，验证了巡视器动力学建模与控制的准确性与实时性。

关键词

火星巡视器; 动力学建模; 拉格朗日递归法; 轮土力学; 仿真软件

目前，火星已经成为世界上各个国家主要的探索对象^[

1]。火星巡视器上装载了多种科学仪器，是人类执行火星表面探测任务的重要载体。由于火星表面环境具有复杂性与不确定性，巡视器应具备穿越各种崎岖地形的能力，从而保证火表探测工作的可靠性^{[参考文献 2

百度学术}2]。例如：美国的“机遇号”火星巡视器因车轮陷入沙土而无法动弹^{[参考文献 3

百度学术}3]；“好奇号”火星巡视器自2012年8月登陆火星至2013年10月执行科学探测任务期间，车轮本体结构破损严重^{[参考文献 4

百度学术}4]。因此，针对防止巡视器发生结构性破坏以及巡视器的运行状态等问题，国内外学者对行星巡视器的运动学、动力学等问题进行了大量深入的研究。

Yoshida等^[

5]从轮土力学和铰接体动力学两方面出发研究了行星车的动力学问题，建立了悬架机构的动力学仿真模型；以车轮滑移率为状态变量，在行星车试验台上进行了实验研究。为了研究在崎岖地形中的高速移动机器人的行为控制方法，Iagnemma等^{[参考文献 6

百度学术}6]开发了一套动力学软件，该软件采用四阶龙格库塔法求解以“牛顿‑欧拉”分析动力学方法建立的非线性多阶常微分方程，并以四轮移动机器人为例，进行了仿真和实验。Lindemann^{[参考文献 7

百度学术}7]通过动态测试模型（Dynamic test model, DTM）评估了火星巡视器在模拟火星环境中的移动特性，同时在ADAMS软件中建立了巡视器的动力学模型并给出了仿真结果。居鹤华等^{[参考文献 8

百度学术}8]研究了六轮摇臂式月面巡视器在复杂月面环境下的动力学建模问题，采用牛顿‑欧拉法建立了巡视器动力学模型，利用开发的三维仿真平台对建立的动力学模型进行了验证。目前为止，国内针对行星巡视器动力学的相关研究较少，仍依赖于ADAMS、ODE、Vortex等商业或开源动力学软件进行二次开发，缺少自主理论和相关技术研发，难以满足复杂的行星表面环境探测的工程需求。

针对复杂火星表面环境及高自由度的火星巡视器动力学建模问题，本文引入了运动链符号系统，其具有伪代码的功能，易于计算机编程，计算效率高；同时，本文方法可以动态地建立巡视器的动力学方程，其具有完全参数化、分层迭代式的优点。

1 火星巡视器结构

一类典型火星巡视器的三维模型如图1所示。巡视器采用主动式摇臂结构，是在被动式摇臂的基础上，在左、右主臂与差速器输出端之间以及左、右副臂与各自的主臂之间各加装一个控制关节，从而实现主臂与车体、副臂与主臂的角度控制功能。车体上的差速器保证其出轴与本体的角度大小相等，方向相反。巡视器6个驱动轮均具有舵机，则该主动式摇臂结构具有16个独立的自由度。6个驱动轮的舵机、主臂关节、副臂关节及差速器均具有绝对编码器，驱动轮转动速度由相对编码器进行测量，用作里程计。巡视器还具有4个单自由度的太阳翼、1个2自由度的数传天线及1个3自由度的桅杆等。

图1 火星巡视器三维模型

Fig.1 3D model of the Mars rover

在机器人工程中，首先定义坐标系，然后通过工程测量确定坐标系间的关系，最后以该坐标系为参考进行机器人的运动学及动力学分析。巡视器的轴间运动关系同样需要通过与它们固结的坐标系进行度量。

为了便于分析，本文采用自然坐标系描述巡视器的轴间关系。巡视器移动系统自然坐标系如图2所示，当系统处于初始位置时，只需建立惯性系以及各轴上的参考点，其他运动轴坐标系也自然确定，而不必为每一运动轴确定各自的体坐标系。

图2 巡视器移动系统自然坐标系

Fig.2 Natural coordinate system of the rover mobile system

图2巡视器各部件的体坐标系编号与部件英文缩略标志符一一对应，原点均位于部件质量中心，x、y、z表示坐标轴，下标i和c分别表示系统惯性系和车体系；下标lr、rr分别表示左、右主臂，lb、rb分别表示左、右副臂；三位下标编号的首位缩略标志l或r分别表示左或右部件，第二位缩略标志f、m、r分别表示前、中、后部件，第三位缩略标志w表示车轮、d表示方向机，例如lfd表示左前方向机，lfw表示左前轮，以此类推；坐标轴 $x_{c}$ 、 $y_{c}$ 轴及 $z_{c}$ 轴分别指向巡视器正前方、正左侧及正上方。

2 火星巡视器动力学建模研究

2.1　基本符号定义

本文以运动链符号演算系统为基础^[

9]，通过精确的符号定义以解决高自由度系统符号规范的问题。为方便阅读，下面对本文出现的相关符号及公式做简单说明。

基本定义：在运动链符号演算系统中，具有偏序的属性变量或常量，在名称上包含表示偏序的指标，包含左上角及右下角指标或者右上角及右下角指标；而它们的方向总是由左上角指标至右下角指标。

符号定义如下：

（1）由至少两个简单运动副依序连接构成的运动机构为运动链，简单运动副是构成运动链的基本单位，称之为链节。定义 ${}^{\bar{l}} k_{l}$ 为运动副，其具有全序性，即 ${}^{\bar{l}} k_{l} = {}^{l}k k_{\bar{l}}$ 。

借鉴集合论的链理论^[

10]，将运动副

{}^{\bar{l}} k_{l}

对应的简单运动链

{}^{\bar{l}} l_{l}

定义为偏序集合

{}^{\bar{l}} l_{l} = (\bar{l}, l]

，其中

\bar{l}

是

l

的前继。同样定义

l_{l}^{\bar{l}} = |(\bar{l}, l]|

，即开集

(\bar{l}, l]

的基数^{[参考文献 10

百度学术}10]。

如图3所示，将连接运动副 ${}^{\bar{l}} k_{l}$ 的两部分分离，则轴 $\bar{l}$ 的运动链为 ${}^{i} l_{\bar{l}}$ ，即图中蓝色虚线所示；轴 $l$ 的闭子树为 ${}^{l} L$ 或 $l L$ ，即图中红色虚线框所示；绿色虚线代表运动副 ${}^{\bar{l}} k_{l}$ 。

图3 轴 $l$ 的闭子树

Fig.3 Closed subtree of Axis $l$

（2）任何复合的运动副可由两个基本运动副组成，即转动副与平动副，分别用 $R$ 和 $P$ 来表示。

记平动符为 $r$ 及转动符为 $ϕ$ ，分别对 ${}^{\bar{l}} l_{l}$ 及 $l_{l}^{\bar{l}}$ 进行实例化，得度量空间下的符号，例如： ${}^{k} r_{l}$ 是 ${}^{k} l_{l}$ 的实例，表示由 $k$ 至 $l$ 的位置矢量； $r_{l}^{k}$ 是 $l_{l}^{k}$ 的实例，表示由 $k$ 至 $l$ 的线位置。表1给出了其他符号的定义。

表1 符号定义

Table 1 Symbol definition

符号	含义	符号	含义
$A$	运动轴序列	$i L$	多轴刚体系统
${}^{\|}$	投影符，例： ${}^{\bar{l}} r_{l}$ 对 $i$ 系的矢量投影为 ${}^{i \| \bar{l}} r_{l}$ ，表征将 $\bar{l}$ 系下的 ${}^{\bar{l}} r_{l}$ 变换到 $i$ 系下表示	$\tilde{□}$	叉乘矩阵符号，例： ${}^{\bar{l}} {\tilde{n}}_{l}$ 代表 ${}^{\bar{l}} n_{l}$ 的叉乘矩阵
${}^{\bar{l}} n_{l}$	运动副 ${}^{\bar{l}} k_{l}$ 的轴矢量	${}^{\bar{l}} n_{l}$	${}^{\bar{l}} n_{l}$ 在 $\bar{l}$ 系下的坐标矢量
$ϕ_{l}^{\bar{l}}$	绕轴 ${}^{\bar{l}} n_{l}$ 的角位置	$r_{l}^{\bar{l}}$	沿轴 ${}^{\bar{l}} n_{l}$ 的线位置
${}^{\bar{l}} ϕ_{l}$	绕 ${}^{\bar{l}} n_{l}$ 转动 $ϕ_{l}^{\bar{l}}$ 的坐标矢量	${}^{\bar{l}} r_{l}$	沿 ${}^{\bar{l}} n_{l}$ 平动到 $r_{l}^{\bar{l}}$ 的坐标矢量
${}^{i} {\dot{ϕ}}_{l}$	轴 $l$ 在惯性系 $i$ 下的绝对角速度	${}^{i} {\dot{r}}_{l I}$	惯性原点到轴 $l$ 质心的速度矢量
${}^{i} {\ddot{ϕ}}_{l}$	轴 $l$ 在惯性系 $i$ 下的绝对角加速度	${}^{i} {\ddot{r}}_{l I}$	惯性原点到轴 $l$ 质心的加速度矢量

2.2　火星巡视器动力学建模

火星巡视器动力学建模分为两个步骤。首先，采用有向Span树^[

9]表征巡视器的拓扑关系，本质上反映轴与轴的连接关系以及轴间运动量的参考关系与作用关系；随后，遍历火星巡视器的轴链拓扑结构有向Span树，完成巡视器的动力学建模。

(1)建立火星巡视器拓扑关系

火星巡视器的运动学与动力学依赖于其拓扑关系分析，本文采用有向Span树来表征巡视器的拓扑关系，即忽略巡视器各部分结构的尺寸，仅考虑运动副和相互连接关系而组成的结构，用Span树的节点编号表达动力学系统的运动副。火星巡视器的有向Span树如图4所示。

图4 火星巡视器轴链有向Span树

Fig.4 Directed Span-tree of the Mars rover axis-chain

火星巡视器通过运动轴的连接确定了一个运动链，该链是有序的运动轴的集合，记22轴巡视器移动系统的运动链为 $A = (0,1, 2, \dots, 22]$ 。

(2)建立火星巡视器动力学方程

一个多自由度的火星巡视器系统可视为一个多轴刚体系统即 $i L$ ，由转动副与平动副组成，下面将推导其拉格朗日递归方程。一般的拉格朗日方程表示为

\frac{{}^{i} d}{d t} (\frac{\partial e_{i L}^{i}}{\partial {\dot{q}}_{k}^{\bar{k}}}) - \frac{\partial e_{i L}^{i}}{\partial q_{k}^{\bar{k}}} = 0

(1)

式中： $e_{i L}^{i}$ 为系统能量； $q_{k}^{\bar{k}}$ 为系统广义坐标。

惯性空间记为 $i$ ， $\forall l \in A$ ；轴 $l$ 的能量记为 $e_{l}^{i}$ ，其中平动动能为 ${}_{v} e_{l}^{i}$ ，转动动能为 ${}_{ω} e_{l}^{i}$ ，引力势能为 ${}_{g} e_{l}^{i}$ ；轴 $l$ 受除引力外的外部合力及合力矩分别为 ${}^{i L} f_{l}$ 及 ${}^{i L} τ_{l}$ ；轴 $l$ 的惯性中心或质心 $I$ 记为 $l I$ ，表示惯性中心 $I$ 是轴 $l$ 的子，则轴 $l$ 的质量及质心转动惯量分别为 $m_{l}$ 及 ${}^{l I} J_{l I}$ ；轴 $l$ 的单位轴矢量为 ${}^{\bar{l}} n_{l}$ ；环境 $i$ 作用于 $l I$ 的惯性加速度记为 ${}^{i} g_{l I}$ 。则系统 $i L$ 能量表示 $e_{L}^{i}$ 为

e_{L}^{i} = {}_{m}e e_{L}^{i} + {}_{g}e e_{L}^{i}

(2)

式中

\{\begin{array}{l} {}_{m} e_{L}^{i} = \sum_{k}^{i L} {}_{v} e_{k}^{i} + \sum_{k}^{i L} {}_{ω} e_{k}^{i} \\ {}_{g} e_{L}^{i} = \sum_{k}^{i L} {}_{g} e_{k}^{i} \\ {}_{v} e_{l}^{i} = \frac{1}{2} \cdot m_{l} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I} \\ {}_{ω} e_{l}^{i} = \frac{1}{2} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l} \\ {}_{g} e_{l}^{i} = m_{l} \cdot {}^{i}r r_{l I}^{T} \cdot {}^{i}g g_{l I} = - m_{l} \cdot {}^{i}r r_{l I}^{T} \cdot {}^{i | l I}g g_{i} \end{array}

(3)

由式(1)得多轴系统拉格朗日方程

\{\begin{matrix} \frac{{}^{i} d}{d t} (\frac{\partial {}_{m}e e_{L}^{i}}{\partial {\dot{r}}_{l}^{\bar{l}}}) - \frac{\partial e_{L}^{i}}{\partial r_{l}^{\bar{l}}} = {}^{i | \bar{l}}n n_{l}^{T} \cdot {}^{i L}f f_{l} {}^{\bar{l}}k k_{l} \in P \\ \frac{{}^{i} d}{d t} (\frac{\partial {}_{m}e e_{L}^{i}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}}) - \frac{\partial e_{L}^{i}}{\partial ϕ_{l}^{\bar{l}}} = {}^{i | \bar{l}}n n_{l}^{T} \cdot {}^{i L}τ τ_{l} {}^{\bar{l}}k k_{l} \in R \end{matrix}

(4)

式(4)为轴 $l$ 的控制方程，即在轴矢量 ${}^{\bar{l}} n_{l}$ 上的力平衡方程； ${}^{i | \bar{l}} n_{l}^{T} \cdot {}^{i L}f f_{l}$ 是合力 ${}^{i L} f_{l}$ 在 ${}^{\bar{l}} n_{l}$ 上的分量， ${}^{i | \bar{l}} n_{l}^{T} \cdot {}^{i L}τ τ_{l}$ 是合力矩 ${}^{i L} τ_{l}$ 在 ${}^{\bar{l}} n_{l}$ 上的分量。

考虑轴 $l$ 平动动能对关节速度 ${\dot{q}}_{l}^{\bar{l}}$ 的偏速度对时间的导数得

\begin{array}{l} \frac{d}{d t} (\frac{\partial}{\partial {\dot{q}}_{l}^{\bar{l}}} ({}_{v}e e_{l}^{i})) = \frac{d}{d t} (\frac{\partial}{\partial {\dot{q}}_{l}^{\bar{l}}} (\frac{1}{2} \cdot m_{l} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I})) = \\ \frac{d}{d t} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I}) = m_{l} \cdot \frac{d}{d t} (\frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}}) \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I} + \\ m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I} = m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial q_{l}^{\bar{l}}} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I} + \\ m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I} = \frac{\partial}{\partial q_{l}^{\bar{l}}} ({}_{v}e e_{l}^{i}) + m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I} \end{array}

(5)

考虑轴 $l$ 转动动能对关节速度 ${\dot{q}}_{l}^{\bar{l}}$ 的偏速度对时间的导数得

\begin{array}{l} \frac{d}{d t} (\frac{\partial}{\partial {\dot{q}}_{l}^{\bar{l}}} ({}_{ω}e e_{l}^{i})) = \frac{d}{d t} (\frac{\partial}{\partial {\dot{q}}_{l}^{\bar{l}}} (\frac{1}{2} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) = \\ \frac{d}{d t} (\frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}) = \frac{d}{d t} (\frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}}) \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l} + \\ \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot ({}^{i | l I}J J_{l I} \cdot {}^{i}{\ddot{ϕ}} {\ddot{ϕ}}_{l} + {}^{i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{l} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}) = \\ \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial q_{l}^{\bar{l}}} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot ({}^{i | l I}J J_{l I} \cdot {}^{i}{\ddot{ϕ}} {\ddot{ϕ}}_{l} + {}^{i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{l} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}) = \\ \frac{\partial}{\partial q_{l}^{\bar{l}}} ({}_{ω}e e_{l}^{i}) + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot ({}^{i | l I}J J_{l I} \cdot {}^{i}{\ddot{ϕ}} {\ddot{ϕ}}_{l} + {}^{i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{l} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}) \end{array}

(6)

因 ${}_{g} e_{l}^{i}$ 与 ${\dot{q}}_{l}^{\bar{l}}$ 不相关，由式(4~6)得

\begin{array}{l} \frac{d}{d t} (\frac{\partial}{\partial {\dot{q}}_{l}^{\bar{l}}} ({}_{v}e e_{l}^{i} + {}_{ω}e e_{l}^{i} + {}_{g}e e_{l}^{i})) - \frac{\partial}{\partial q_{l}^{\bar{l}}} ({}_{v}e e_{l}^{i} + {}_{ω}e e_{l}^{i} + {}_{g}e e_{l}^{i}) = \\ \frac{d}{d t} (\frac{\partial}{\partial {\dot{q}}_{l}^{\bar{l}}} ({}_{v}e e_{l}^{i} + {}_{ω}e e_{l}^{i})) - \frac{\partial}{\partial q_{l}^{\bar{l}}} ({}_{v}e e_{l}^{i} + {}_{ω}e e_{l}^{i} + {}_{g}e e_{l}^{i}) = \\ m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot \frac{d}{d t} ({}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}) - \\ \frac{\partial}{\partial q_{l}^{\bar{l}}} ({}_{g}e e_{l}^{i}) = m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot \\ ({}^{i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{l} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l} + {}^{i | l I}J J_{l I} \cdot {}^{i}{\ddot{ϕ}} {\ddot{ϕ}}_{l}) - \frac{\partial}{\partial q_{l}^{\bar{l}}} ({}_{g}e e_{l}^{i}) \end{array}

(7)

动力学系统 $i L$ 的平动动能及转动动能分别表示为

\begin{matrix} {}_{v} e_{L}^{i} = \sum_{l}^{i L} {}_{v} e_{l}^{i} & {}_{ω} e_{L}^{i} = \sum_{l}^{i L} {}_{ω} e_{l}^{i} \end{matrix}

(8)

考虑式(2)和式(3)，即有

\begin{array}{l} \frac{d}{d t} (\frac{\partial}{\partial {\dot{q}}_{l}^{\bar{l}}} (e_{L}^{i})) - \frac{\partial}{\partial q_{l}^{\bar{l}}} (e_{L}^{i}) = \sum_{l}^{{}^{i} L} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I}) + \\ \sum_{l}^{i L} (\frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{q}}_{l}^{\bar{l}}} \cdot ({}^{i} {\dot{\tilde{ϕ}}}_{l} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l} + {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) - \sum_{l}^{i L} (\frac{\partial}{\partial q_{l}^{\bar{l}}} ({}_{g}e e_{l}^{i})) \end{array}

(9)

下面求系统能量对关节速度及坐标的偏速度。若 ${}^{\bar{l}} k_{l} \in P$ ，并考虑 $l \in l L \subset i L$ ，由式(2)及式(3)，得

\begin{array}{l} \frac{\partial}{\partial {\dot{r}}_{l}^{\bar{l}}} ({}_{v}e e_{L}^{i}) = \frac{\partial}{\partial {\dot{r}}_{l}^{\bar{l}}} (\sum_{l}^{l L} (\frac{1}{2} \cdot m_{l} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I})) = \\ \sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{r}}_{l}^{\bar{l}}} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I}) \end{array}

(10)

\begin{array}{l} \frac{\partial}{\partial {\dot{r}}_{l}^{\bar{l}}} ({}_{ω}e e_{L}^{i}) = \frac{\partial}{\partial {\dot{r}}_{l}^{\bar{l}}} (\sum_{l}^{l L} (\frac{1}{2} {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) = \\ \sum_{l}^{l L} (\frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{r}}_{l}^{\bar{l}}} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}) \end{array}

(11)

\begin{array}{l} \frac{\partial}{\partial r_{l}^{\bar{l}}} ({}_{g}e e_{L}^{i}) = \frac{\partial}{\partial r_{l}^{\bar{l}}} (\sum_{l}^{l L} (m_{l} \cdot {}^{i}r r_{l I}^{T} \cdot {}^{i}g g_{l I})) = \\ \sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}r r_{l I}^{T}}{\partial r_{l}^{\bar{l}}} \cdot {}^{i}g g_{l I}) \end{array}

(12)

若 ${}^{\bar{l}} k_{l} \in R$ ，并考虑 $l \in l L \subset i L$ ，由式(2)及式(3)，得

\begin{array}{l} \frac{\partial}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} ({}_{v}e e_{L}^{i}) = \frac{\partial}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} (\sum_{l}^{l L} (\frac{1}{2} \cdot m_{l} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I})) = \\ \sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I}) \end{array}

(13)

\begin{array}{l} \frac{\partial}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} ({}_{ω}e e_{L}^{i}) = \frac{\partial}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} (\sum_{k}^{l L} (\frac{1}{2} {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) = \\ \sum_{l}^{l L} (\frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}) \end{array}

(14)

\begin{array}{l} \frac{\partial}{\partial ϕ_{l}^{\bar{l}}} ({}_{g}e e_{L}^{i}) = \frac{\partial}{\partial ϕ_{l}^{\bar{l}}} (\sum_{l}^{l L} (m_{l} \cdot {}^{i}r r_{l I}^{T} \cdot {}^{i}g g_{l I})) = \\ \sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}r r_{l I}^{T}}{\partial ϕ_{l}^{\bar{l}}} \cdot {}^{i}g g_{l I}) \end{array}

(15)

至此，已完成能量对关节速度及坐标的偏速度计算。

接着，对时间求绝对导数。若 ${}^{\bar{l}} k_{l} \in P$ ，由式（9~11），得

\begin{array}{l} \frac{{}^{i} d}{d t} (\frac{\partial}{\partial {\dot{r}}_{l}^{\bar{l}}} ({}_{m}e e_{L}^{i})) = \\ \frac{d}{d t} (\sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{r}}_{l}^{\bar{l}}} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{r}}_{l}^{\bar{l}}} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) = \\ \frac{\partial}{\partial r_{l}^{\bar{l}}} ({}_{m}e e_{L}^{i}) + \\ \sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{r}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{r}}_{l}^{\bar{l}}} \cdot ({}^{i | l I} J_{l I} \cdot {}^{i}{\ddot{ϕ}} {\ddot{ϕ}}_{l} + {}^{i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{l} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) \end{array}

(16)

若 ${}^{\bar{l}} k_{l} \in R$ ，由式(9、13、14)，得

\begin{array}{l} \frac{d}{d t} (\frac{\partial}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} ({}_{m}e e_{L}^{i})) = \\ \frac{d}{d t} (\sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot {}^{i}{\dot{r}} {\dot{r}}_{l I} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) = \\ \frac{\partial}{\partial ϕ_{l}^{\bar{l}}} ({}_{m}e e_{L}^{i}) + \\ \sum_{l}^{l L} (m_{l} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{l I}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{l I} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot ({}^{i | l I} J_{l I} \cdot {}^{i}{\ddot{ϕ}} {\ddot{ϕ}}_{l} + {}^{i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{l} \cdot {}^{i | l I}J J_{l I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{l})) \end{array}

(17)

至此，已完成对时间 $t$ 的求导。

由偏导数链规则可知^[

\begin{array}{l} \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{n S}}{\partial {\dot{ϕ}}_{k}^{\bar{k}}} = {}^{i | \bar{k}}{\tilde{n}} {\tilde{n}}_{k} \cdot {}^{i | k}r r_{n S} = - {}^{i | k}{\tilde{r}} {\tilde{r}}_{n S} \cdot {}^{i | \bar{k}}n n_{k} \\ \begin{matrix} \begin{matrix} \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{n}}{\partial {\dot{ϕ}}_{k}^{\bar{k}}} = {}^{i | \bar{k}}n n_{k}, & \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{n}}{\partial {\dot{r}}_{k}^{\bar{k}}} = 0_{3} \end{matrix} \\ \frac{\partial {}^{i}r r_{n S}}{\partial r_{k}^{\bar{k}}} = \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{n S}}{\partial {\dot{r}}_{k}^{\bar{k}}} = {}^{i | \bar{k}}n n_{k} \end{matrix} \end{array}

(18)

根据式(18)偏导数关系，将式（12，15~17）代入式(9)，可以推导出轴 $l$ 的拉格朗日递归动力学方程如下

\{\begin{array}{l} \sum_{k}^{i L} (m_{k} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{k I}^{T}}{\partial {\dot{r}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{k I} - m_{k} \cdot \frac{\partial {}^{i}r r_{k I}^{T}}{\partial r_{k}^{\bar{k}}} \cdot {}^{i}g g_{k I}) = \\ \begin{matrix} {}^{i | \bar{l}}n n_{l}^{T} \cdot {}^{i L}f f_{l} & \begin{matrix} {}^{\bar{l}} k_{l} \in P \end{matrix} \end{matrix} \\ \sum_{k}^{i L} (m_{k} \cdot \frac{\partial {}^{i}{\dot{r}} {\dot{r}}_{k I}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot {}^{i}{\ddot{r}} {\ddot{r}}_{k I} + \frac{\partial {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{k}^{T}}{\partial {\dot{ϕ}}_{l}^{\bar{l}}} \cdot ({}^{i | k I} J_{k I} \cdot {}^{i}{\ddot{ϕ}} {\ddot{ϕ}}_{k} + {}^{i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{k} \cdot {}^{i | k I}J J_{k I} \cdot {}^{i}{\dot{ϕ}} {\dot{ϕ}}_{k}) - \\ m_{k} \cdot \frac{\partial {}^{i}r r_{k I}^{T}}{\partial ϕ_{l}^{\bar{l}}} \cdot {}^{i}g g_{k I}) = {}^{i | \bar{l}}n n_{l}^{T} \cdot {}^{i L}τ τ_{l} {}^{\bar{l}}k k_{l} \in R \end{array}

(19)

由式(19)可知，对于多轴系统 $i L$ ，遍历火星巡视器的轴链拓扑结构有向Span树即可完成对巡视器的动力学建模。由此可见，针对高自由度的火星巡视器动力学建模问题，本文的动力学方法可以动态地建立基于运动轴的有向Span树，具有完全参数化、分层迭代式的优点，较传统的牛顿‑欧拉法、拉格朗日法、凯恩法^[

11‑13]及目前国际上的6D空间算子代数法^{[参考文献 14‑15}14‑15]具有更低的计算复杂度；同时，该方法引入了运动链符号系统，具备了伪代码的功能，使得计算机编程更加简单，计算效率高于传统的拉格朗日递归法，适用于高自由度、可变拓扑的机器人动力学建模问题。

在建立系统动力学方程后，紧接着就是方程求解的问题，此为本文后续的工作重点，根据文献[

9]，可将式(19)更进一步表达为规范型方程，并进行多轴系统动力学求解，对此本文不再展开讨论。

3 轮土力学与轮式系统移动维度研究

轮土作用力是轮式多轴系统移动的动因，利用轮土力学知识预测火星巡视器的车轮‑土壤相互作用性能，对巡视器的轮土作用力识别、动力学控制具有重要意义。Bekker轮土模型目前应用广泛，但该模型没有定义驱动轮与土壤相互作用过程中相关物理参数的参考系。因此，本文将链符号系统引入经典Bekker轮土模型，以提高其准确性。

3.1　Bekker轮土力学

Bekker^[

16]轮土力学模型是专门为弹塑性土和刚性驱动轮设计的作用力模型，该模型建立基于水平地形及重力方向载荷。

(1)轮土作用几何关系

如图5所示，给定轴链 $(\bar{\bar{l}}, \bar{l}, l]$ ，记舵机轴系为 $\bar{\bar{l}}$ ，轮地接触坐标系记为 $l$ ，它与环境 $i$ 固结， $O_{l}$ 是轮法向力与轮外侧表面的交点。因轮轴无绝对编码器，故以轮轴系为参考。

图5 轮土模型几何关系

Fig.5 Geometrical relationship of terramechanics model

其中： $θ_{m}$ 为最大应力角； $θ_{1}$ 为泥土进入角； $θ_{2}$ 为泥土退出角； $r$ 为轮半径； $κ_{1}$ 为进入剪应力； $κ_{2}$ 为退出剪应力； $z$ 为轮沉陷量； $σ_{1}$ 为进入正应力； $σ_{2}$ 为退出正应力； ${\dot{ϕ}}_{l}^{\bar{l}}$ 为轮转动速度； ${}^{i} {\dot{r}}_{l}^{[1]}$ 为轮前向速度； ${}^{\bar{l}} τ_{l}^{[2]}$ 为轮轴向力矩； ${}^{\bar{l}} f_{l}^{[3]}$ 为轮法向力； ${}^{\bar{l}} f_{l}^{[1]}$ 为轮牵引力。

根据图5所示几何关系，可得

\{\begin{array}{l} {}^{l |} e_{l}^{[x]} = {}^{l | \bar{\bar{l}}}{\tilde{n}} {\tilde{n}}_{\bar{l}} \cdot {}^{l | \bar{l}}n n_{l} \\ {}^{l |} e_{l}^{[z]} = {}^{l | \bar{l}}n n_{l} \end{array}

(20)

{}^{i | \bar{l}} {\dot{r}}_{l} = {}^{i | l}{\dot{r}} {\dot{r}}_{i} = - {}^{i}{\dot{r}} {\dot{r}}_{l}

(21)

则轮速 $v$ 表达为

\begin{array}{l} v = 0 - {({}^{l | \bar{\bar{l}}} {\tilde{n}}_{\bar{l}} \cdot {}^{l | \bar{l}}n n_{l})}^{T} \cdot ({}^{l | i} {\dot{r}}_{l} + {}^{l | i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{\bar{l}} \cdot {}^{l | \bar{l}}n n_{l} \cdot r) = \\ {}^{l | \bar{l}} n_{l}^{T} \cdot {}^{l | \bar{\bar{l}}}{\tilde{n}} {\tilde{n}}_{\bar{l}} \cdot ({}^{l | i} {\dot{r}}_{l} + {}^{l | i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{\bar{l}} \cdot {}^{l | \bar{l}}n n_{l} \cdot r) \end{array}

(22)

轮 $\bar{l}$ 相对接触点（零速点） $O_{l}$ 的滑移率 $s_{\bar{l}}^{l}$ 为

s_{\bar{l}}^{l} = \{\begin{matrix} - \frac{{}^{l | \bar{l}} n_{l}^{T} \cdot {}^{l | \bar{\bar{l}}}{\tilde{n}} {\tilde{n}}_{\bar{l}} \cdot ({}^{l | i} {\dot{r}}_{l} + {}^{l | i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{\bar{l}} \cdot {}^{l | \bar{l}}n n_{l} \cdot r)}{{}^{l | \bar{l}} n_{l}^{T} \cdot {}^{l | \bar{\bar{l}}}{\tilde{n}} {\tilde{n}}_{\bar{l}} \cdot {}^{l | i}{\dot{\tilde{ϕ}}} {\dot{\tilde{ϕ}}}_{\bar{l}} \cdot {}^{l | \bar{l}}n n_{l} \cdot r} \\ \begin{matrix} {}^{i} {\dot{ϕ}}_{\bar{l}} \neq 0_{3} \end{matrix} \\ 0 \\ \begin{matrix} {}^{i} {\dot{ϕ}}_{\bar{l}} = 0_{3} \end{matrix} \end{matrix}

(23)

(2)改进的轮土应力模型

根据文献[

17]，可得正应力

σ (z)

及剪应力

κ (z)

的表达式，Yoshida等^{[参考文献 18

百度学术}18]应用法向最大应力角

θ_{m}

、泥土进入角

θ_{1}

、轮半径

r

、滑移率

s_{\bar{l}}^{l}

将

σ (z)

及

κ (z)

的表达式修改为

\{\begin{array}{l} σ_{1} (θ) = (k_{1} + k_{2} b) {(\frac{r}{b})}^{n} {(C (θ) - C (θ_{1}))}^{n} \\ θ_{m} \leq θ < θ_{1} \\ σ_{2} (θ) = (k_{1} + k_{2} b) {(\frac{r}{b})}^{n} \cdot \\ {[C (θ_{1} - \frac{θ}{θ_{m}} (θ_{1} - θ_{m})) - C (θ_{1})]}^{n} θ_{2} < θ \leq θ_{m} \end{array}

(24)

\begin{array}{l} κ (θ) = (c + σ t a n ϕ) (1 - \begin{matrix}  \end{matrix} \\ e (- \frac{r}{k} (θ_{1} - θ - (1 - s) (S (θ_{1}) - S (θ))))) \end{array}

(25)

\{\begin{array}{l} \begin{matrix} σ_{1} (θ) = \frac{θ_{1} - θ}{θ_{1} - θ_{m}} σ_{m}, & σ_{2} (θ) = \frac{θ}{θ_{m}} σ_{m} \end{matrix} \\ \begin{matrix} κ_{1} (θ) = \frac{θ_{1} - θ}{θ_{1} - θ_{m}} τ_{m}, & κ_{2} (θ) = c + \frac{θ}{θ_{m}} (τ_{m} - c) \end{matrix} \end{array}

(26)

θ_{m} = (c_{1} + c_{2} s_{\bar{l}}^{l}) θ_{1}

(27)

式中： $S (θ)$ 和 $C (θ)$ 分别代表 $θ$ 的正弦与余弦；其他的轮土参数如表2所示

表2 轮土参数

Table 2 Terramechanics parameters

符号	含义（单位）	符号	含义（单位）
$c$	粘滞系数（N/m²）	$n$	深陷指数
$k_{1}$	摩擦参数	$k_{2}$	粘滞参数
$b$	轮宽（m）	$j$	剪位移
$ϕ$	内摩擦角（rad）	$k$	剪变形参数（m）

给定轮土参数及滑移率，由式(27)得到 $θ_{m}$ ，由式(24)及式(25)分别得 $σ_{m}$ 及 $κ_{m}$ ；然后，由式(26)得 $σ_{1} (θ)$ ， $σ_{2} (θ)$ ， $κ_{1} (θ)$ 及 $κ_{2} (θ)$ 。

由正应力分量 $σ (θ) \cdot C (θ)$ 、剪应力分量 $κ (θ) \cdot S (θ)$ 及式(24)至式(26)，得到轮地接触点 $O_{l}$ 对轮 $\bar{l}$ 的正压力 $N_{\bar{l}}^{l}$ 表达式

\begin{array}{l} N_{\bar{l}}^{l} = \frac{r b}{θ_{m} (θ_{1} - θ_{m})} \cdot \\ [σ_{m} (- θ_{m} C (θ_{1}) + θ_{1} C (θ_{m}) - θ_{1} + θ_{m}) - \\ κ_{m} (θ_{m} S (θ_{1}) - θ_{1} S (θ_{m})) - \\ c (θ_{1} S (θ_{m}) - θ_{m} S (θ_{m}) - θ_{m} θ_{1} + θ_{m}^{2})] \end{array}

(28)

可以证明正压力 $N_{\bar{l}}^{l}$ 是关于进入角 $θ_{1}$ 的单调函数，即

\frac{\partial (N_{\bar{l}}^{l})}{\partial θ_{1}} \geq 0

(29)

(3)牵引力、牵引力矩及侧向力

由式(26)得牵引力矩 $τ_{\bar{l}}^{l}$ 及牵引力 $T_{\bar{l}}^{l}$ 表达式

\begin{array}{l} τ_{\bar{l}}^{l} = r^{2} b \int_{θ_{2}}^{θ_{1}} κ (θ) d θ = \\ r^{2} b \cdot [\int_{0}^{θ_{m}} (c + \frac{θ}{θ_{m}} (κ_{m} - c)) d θ + \int_{θ_{m}}^{θ_{1}} \frac{θ_{1} - θ}{θ_{1} - θ_{m}} κ_{m} d θ] = \\ r^{2} b \cdot (c θ_{m} + \frac{κ_{m} - c}{2 θ_{m}} θ_{m}^{2} + \frac{κ_{m}}{2} (θ_{1} - θ_{m})) = \\ \frac{r^{2} b}{2} \cdot (κ_{m} θ_{1} + c θ_{m}) \end{array}

(30)

T_{\bar{l}}^{l} = τ_{\bar{l}}^{l} / r = \frac{r b}{2} \cdot (κ_{m} θ_{1} + c θ_{m})

(31)

如图6所示，滑移角 $α_{\bar{l}}^{l}$ 由轮轴速度 ${}^{l | i} r_{\bar{l}}^{[1]}$ 及轮地接触系坐标轴 $x_{l}$ 确定。

图6 滑移角

Fig.6 Slip angle

轮侧向力 $L_{\bar{l}}^{l}$ 表达为^[

17]

L_{\bar{l}}^{l} = (c b A + N_{\bar{l}}^{l} \cdot t a n (ϕ)) \cdot (1 - e (- B \cdot α_{\bar{l}}^{l}))

(32)

式中： $A$ 为轮地接触面（m²）； $α_{\bar{l}}^{l}$ 为滑移角； $B$ 为常数（1/rad）； $N_{\bar{l}}^{l}$ 为正压力（N）。

3.2　轮地接触约束方程

理想的轮地接触可视为接触副，该接触副仅约束了关节速度而不是位形，故接触副存在单边约束，即轮地正压力总大于等于零。

记 $u^{″}$ 矢量空间为接触副 ${}^{u^{″}} O_{k^{″}}$ 的约束矢量空间，它是 $u$ 的子空间， $u^{″}$ 位形空间除维度不同外，与矢量空间 $u$ 基本一致。因此，可以定义约束副 ${}^{u^{″}} C_{u}$ 。

接触副 ${}^{u^{″}} O_{k^{″}}$ 的约束要求体 $u^{″}$ 不能刺入体 $k^{″}$ 中，故有

{}^{u^{″}} C_{u} \cdot {}^{u}{\ddot{γ}} {\ddot{γ}}_{u S} \geq 0

(33)

式中 ${}^{u} {\ddot{γ}}_{u S}$ 为接触点 $u S$ 的加速度。

同时，接触力 $f_{k^{″} S}^{u^{″} S}$ 非负，即

f_{k^{″} S}^{u^{″} S} \geq 0

(34)

记接触面法向支撑力为 ${}^{u^{″} S} f_{k^{″} S}$ ，接触条件为

{}^{i | u^{″} S} f_{k^{″} S}^{T} \cdot {}^{u^{″}}C C_{u} \cdot {}^{i | u}{\ddot{γ}} {\ddot{γ}}_{u S} = 0

(35)

由式(33)至式(35)表明接触副约束是不等式方程组。

下面建立接触副 ${}^{u^{″} S} O_{k^{″} S}$ 的单边约束方程。考虑接触副 ${}^{u^{″} S} O_{k^{″} S}$ ，接触点记为 $u^{″} S$ ，接触力记为 $f_{u^{″} S}^{u}$ ，接触点位置记为 ${}^{u} r_{u^{″} S}$ 。接触副 ${}^{u^{″} S} O_{k^{″} S}$ 是特殊的运动副 ${}^{u^{″} S} k_{k^{″} S}$ ，其中一个平动轴受单边约束，转动轴受控。记 $u^{″}$ 对 $k^{″}$ 施加的 ${}^{i | u^{″} S} τ_{k^{″} S}$ 或 ${}^{i | u^{″} S} f_{k^{″} S}$ 为广义约束力， ${}^{i | u^{″}} J_{u S}$ 及 ${}^{i | k^{″}} J_{k S}$ 为广义约束雅可比矩阵。考虑单位质量及转动惯量，则有

[\begin{matrix} {\ddot{r}}_{k^{″} S}^{u^{″} S} \\ {\ddot{r}}_{u^{″} S}^{k^{″} S} \end{matrix}] = [\begin{matrix} {}^{i | u^{″}} J_{u S} & {}^{i | k^{″}} J_{k S} \end{matrix}] \cdot [\begin{matrix} {}^{i} {\ddot{γ}}_{u} \\ {}^{i} {\ddot{γ}}_{k} \end{matrix}] = [\begin{matrix} {}_{c} τ_{k^{'} S}^{u^{'} S} \\ {}_{c} f_{k^{'} S}^{u^{'} S} \end{matrix}]

(36)

[\begin{matrix} {}^{i | u^{″} S} W_{k^{″} S} \\ {}^{i | k^{″} S} W_{u^{″} S} \end{matrix}] = {[\begin{matrix} {}^{i | u^{″}} J_{u S} & {}^{i | k^{″}} J_{k S} \end{matrix}]}^{T} \cdot [\begin{matrix} τ_{k^{″} S}^{u^{″} S} \\ f_{k^{″} S}^{u^{″} S} \end{matrix}]

(37)

式中

\{\begin{array}{l} {}^{i | u^{″}} J_{u S} = [\begin{matrix} {}_{R}^{i | u^{″}} J_{u S} & {}_{P}^{i | u^{″}} J_{u S} \end{matrix}] \\ {}^{i | k^{″}} J_{k S} = [\begin{matrix} {}_{R}^{i | k^{″}} J_{k S} & {}_{P}^{i | k^{″}} J_{k S} \end{matrix}] \\ {}_{R}^{i | u^{″}} J_{u S} = - {}^{i | u}n n_{u^{″}}^{T} \cdot {}^{i | u}{\tilde{r}} {\tilde{r}}_{u^{″} S} \\ {}_{R}^{i | k^{″}} J_{k S} = {}^{i | k}n n_{k^{″}}^{T} \cdot {}^{i | k}{\tilde{r}} {\tilde{r}}_{k^{″} S} \\ {}_{P}^{i | u^{″}} J_{u S} = - {}^{i | u}n n_{u^{″}}^{T} \\ {}_{P}^{i | k^{″}} J_{k S} = {}^{i | k}n n_{k^{″}}^{T} \end{array}

(38)

\{\begin{array}{l} {}^{i | u^{″} S} W_{k^{″} S} = [\begin{matrix} {}^{i | u^{″} S} τ_{k^{″} S} \\ {}^{i | u^{″} S} f_{k^{″} S} \end{matrix}] \\ {}^{i | k^{″} S} W_{u^{″} S} = [\begin{matrix} {}^{i | k^{″} S} τ_{u^{″} S} \\ {}^{i | k^{″} S} f_{u^{″} S} \end{matrix}] \\ W_{k^{″} S}^{u^{″} S} = [\begin{matrix} τ_{k^{″} S}^{u^{″} S} \\ f_{k^{″} S}^{u^{″} S} \end{matrix}] \end{array}

(39)

由式(33~39)得

- [\begin{matrix} {}^{i | u^{″}} J_{u S} & {}^{i | k^{″}} J_{k S} \end{matrix}] \cdot [\begin{matrix} {}^{i} {\ddot{γ}}_{u} \\ {}^{i} {\ddot{γ}}_{k} \end{matrix}] - {}_{c}f f_{k^{″} S}^{u^{″} S} \geq 0

(40)

- f_{k^{″} S}^{u^{″} S} \geq 0

(41)

f_{k^{″} S}^{u^{″} S} \cdot {[\begin{matrix} {}^{i | u^{″}} J_{u S} & {}^{i | k^{″}} J_{k S} \end{matrix}]}^{T} \cdot [\begin{matrix} {}^{i} {\ddot{γ}}_{u} \\ {}^{i} {\ddot{γ}}_{k} \end{matrix}] = 0

(42)

至此，式(40~42)建立了接触副 ${}^{u^{″} S} O_{k^{″} S}$ 的单边约束方程。显然，轮地接触的前向牵引力及侧向力是正压力的函数，轮地作用的合力与正压力成正比。因此，该合力也是单边作用力，通过一个单边约束方程表征。

基于上述Bekker轮土力学代数方程及式(19)，可得完整的巡视器动力学模型，其中包含与25个轴对应的25个动力学方程、与一个差速器对应的一个约束方程及与6个轮地接触点对应的6个轮地接触单边约束方程。

4 仿真平台实现

基于VS2013软件中的MFC框架，采用C++语言开发了火星巡视器三维可视化仿真平台，在仿真平台中，建立了模拟地形场景及一类典型火星巡视器三维模型。仿真平台主要动力学参数中，巡视器车体质量127 kg，相对于本体系的三轴主转动惯量为[40,40,40] kg·m²，车轮质量为3.674 kg，主转动惯量为[2,2,2] kg·m²，其余部件的质量和转动惯量较小。以如图7和图8所示的火星表面模拟环境为场景，进行可视化仿真分析，巡视器可根据期望的移动路线、速度及加速度移动，同时在期望路径上完成原地转向、直线移动、圆弧运动等基本行为。

图7 巡视器沿曲线行驶

Fig.7 Rover trajectory: curvature steering

图8 巡视器原地转向

Fig.8 Rover trajectory: circling around

图9~14给出了巡视器在图7模拟场景下的运动状态及对应的方向机力矩、角度，车轮角速度、力矩，巡视器姿态及速度曲线图。其中，每幅图的点划线、虚线及实线分别与火星巡视器部分结构的前、中及后部件对应。

图9 方向机力矩

Fig.9 Steering gear torque

图10 方向机角度

Fig.10 Steering gear angle

图11 车轮角速度

Fig.11 Angular velocity of wheel

图12 驱动轮力矩

Fig.12 Wheel torque

图13 巡视器姿态

Fig.13 Attitude of the Mars rover

图14 巡视器速度

Fig.14 Speed of the Mars rover

火星表面重力加速度设为3.72 $m / s^{2}$ ，图9~14表明巡视器由起始至5.0 s时的动力学过程，包含500仿真步，每步10 ms。图9表明巡视器在行驶过程中方向机力矩变化较大，这是车轮接触崎岖地形导致的结果。图10表明两侧方向机角度变化幅度较小，因为前5 s内火星巡视器刚开始与地面接触。图11表明左侧车轮角速度由-2.0 (º)/s至2.1 (º)/s变化，右侧车轮角速度由-1.2 (º)/s至2.1 (º)/s变化；且右侧车轮角速度变化幅度大于左侧，符合转向特性。图12反映6个车轮所受到轮土力矩随着轮土作用而逐渐变化，且轮土支撑力之和小于巡视器在火表加速度下的总重力。图13和图14反映了巡视器的前5 s姿态坐标是(200,200,0)，且巡视器以约0.1 m/s速度前进时三轴的平动速度及转动速度，其中巡视器的偏航角变化较大，这是由于巡视器在崎岖的火表地形沿曲线行驶，前侧驱动轮较后侧驱动轮的沉陷量大，反映了地形模拟的真实性。

由于火星表面环境具有复杂性与不确定性，传统的被动六轮摇臂-转向架火星巡视器只能被动地适应地形，在不确定强的地表不能进行自我调节，例如：巡视器行驶在松软地形时易发生车轮沉陷等问题。本文巡视器采用主动式摇臂结构，具有主臂与车体、副臂与主臂的角度控制功能，例如：在复杂的火表中，通过巡视器悬架的主动变形，改变了巡视器移动系统的质心相对位置，可以实现车体高度调整、抬轮、“蟹形”移动等，增加了对不确定性地表的适应能力，更加适合巡视器探测崎岖复杂地形。

5 结　　论

针对中国未来行星表面探测任务的复杂性与不确定性，本文对火星巡视器动力学原理进行了研究。首先对火星巡视器的结构进行介绍并采用自然坐标系统描述巡视器移动系统。然后对本文的基本符号进行定义，建立了巡视器的树形拓扑结构，应用链符号系统推导了多轴拉格朗日方程，计算能量对关节速度和坐标的偏速度及对时间的导数，遍历巡视器拓扑结构即可完成巡视器的动力学建模，该动力学方法具有完全参数化、分层迭代式的优点。接着，通过改进的轮土应力模型，得到轮地作用力和力矩，同时建立了轮地作用接触副的约束方程，用于后续巡视器动力学建模。最后，在开发的可视化仿真平台对巡视器模型进行测试。结果表明，巡视器可以根据期望的速度及加速度移动，同时可在平坦或崎岖地形完成原地转向、直线、圆弧移动等基本行为，体现了本文动力学方法的有效性与优越性，对中国未来行星车动力学研究及仿真软件的开发具有借鉴意义。

参考文献

孙泽洲，孟林智. 中国深空探测现状及持续发展趋势［J］. 南京航空航天大学学报， 2015， 47（6）： 785-791. [百度学术]

SUN Zezhou， MENG Linzhi. Current situation and sustainable development trend of deep space exploration in China［J］. Journal of Nanjing University of Aeronautics & Astronautics， 2015， 47（6）： 785-791. [百度学术]

BIESIADECKI J J， MAIMONE M W. The Mars exploration rover surface mobility flight software driving ambition［C］//Proceedings of 2006 IEEE Aerospace Conference. [百度学术]

［S.l.］： IEEE， 2006： 1-15. [百度学术]

CALLAS J L， GOLOMBEK M P， FRAEMAN A A. Mars exploration rover opportunity end of mission report［R］. Pasadena， CA： Jet Propulsion Laboratory， National Aeronautics and Space Administration， 2019. [百度学术]

TOUPET O， BIESIADECKI J， RANKIN A， et al. Terrain-adaptive wheel speed control on the curiosity Mars rover： Algorithm and flight results［J］. Journal of Field Robotics， 2020， 37（5）： 699-728. [百度学术]

YOSHIDA K， HAMANO H. Motion dynamics of a rover with slip-based traction model［C］//Proceedings of 2002 IEEE International Conference on Robotics and Automation. [百度学术]

［S.l.］： IEEE， 2002， 3： 3155-3160. [百度学术]

IAGNEMMA K， GOLDA D， SPENKO M， et al. Experimental study of high-speed rough-terrain mobile robot models for reactive behaviors［M］//Experimental Robotics VIII. Berlin， Heidelberg： Springer， 2003： 654-663. [百度学术]

LINDEMANN R. Dynamic testing and simulation of the Mars exploration rover［C］//Proceedings of ASME 2005 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. California： ASME Digital Collection， 2005， 47438： 99-106. [百度学术]

居鹤华，田小二. 月面巡视器实时动力学建模与牵引控制［J］. 宇航学报， 2014， 35（7）： 743-752. [百度学术]

JU Hehua， TIAN Xiaoer. Real-time dynamics modeling and traction control for lunar rover［J］. Journal of Astronautics， 2014， 35（7）： 743-752. [百度学术]

JU Hehua. Hand book of space robotics： Axis-invariant based MAS modelling， planning and control［M］. New York： Springer Press， 2018： 43-61. [百度学术]

JECH T. Set theory［M］. ［S.l.］： Springer Science & Business Media， 2013： 70-87. [百度学术]

谢传峰，王琪，程耀，等. 动力学［M］. 北京：高等教育出版社， 2006： 137-161. [百度学术]

XIE Chuanfeng， WANG Qi， CHENG Yao， et al. Dynamics［M］. Beijing： Higher Education Press， 2006： 137-161. [百度学术]

ORIN D E， MCGHEE R B， VUKOBRATOVIC M， et al. Kinematic and kinetic analysis of open-chain linkages utilizing Newton-Euler methods［J］. Mathematical Bioences， 1979， 43（1/2）： 107-130. [百度学术]

KANE T R， LEVINSON D A. Dynamics theory and applications［M］. ［S.l.］： The Internet-First University Press， 2005. [百度学术]

RODRIGUEZ G， JAIN A， KREUTZ-DELGADO K. Spatial operator algebra for multibody system dynamics［J］. Journal of the Astronautical Sciences， 1992， 40（1）： 27-50. [百度学术]

JAIN A. Robot and multibody dynamics： Analysis and algorithms［M］. New York： Springer， 2012： 209-227. [百度学术]

BEKKER M G. Introduction to terrain-vehicle systems［M］. Ann Arbor： University of Michigan Press， 1969： 71-92. [百度学术]

CROLLA D A， EL-RAZAZ A S A. A review of the combined lateral and longitudinal force generation of tyres on deformable surfaces［J］. Journal of Terramechanics， 1987， 24（3）： 199-225. [百度学术]

YOSHIDA K， WATANABE T， MIZUNO N， et al. Terramechanics-based analysis and traction control of a lunar/planetary rover［C］//Proceedings of Field and Service Robotics. Berlin， Heidelberg： Springer， 2003： 225-234. [百度学术]

请使用 Firefox、Chrome、IE10、IE11、360极速模式、搜狗极速模式、QQ极速模式等浏览器，其他浏览器不建议使用!

火星巡视器动力学建模研究与仿真平台实现 PDF

摘要

关键词

1 火星巡视器结构

2 火星巡视器动力学建模研究

2.1 基本符号定义

2.2 火星巡视器动力学建模