轮控小行星探测器建模及跳跃移动仿真

轮控小行星探测器建模及跳跃移动仿真

龚宇鹏;谢天;林晓辉

【摘 要】小行星探测器在弱引力环境下无法采用传统的轮式机构进行移动.为解决

该问题,采用反作用飞轮对探测器跳跃移动进行控制,并分析了该方案的可行性.根据

Hertz碰撞定律及简化的Karnopp切向摩擦力模型,建立了探测器与地面的接触力

模型.分析了轮控小行星探测器的起跳过程,给出了探测器静止起跳所需要的最小飞

轮力矩关系.考虑到反作用飞轮存在惯性、粘滞、摩擦等情况,建立了轮控探测器的

姿态动力学模型,并对探测器在均匀重力场下的连续跳跃过程进行了控制策略设计

及仿真.结果 表明:基于飞轮控制的小行星探测器跳跃移动在微重力环境下是可行的,

且可以通过施加合适的控制力矩维持探测器跳跃的方向及跳跃过程的稳定性.

【期刊名称】《上海航天》

【年(卷),期】2019(036)003

【总页数】7页(P61-67)

【关键词】小行星探测;跳跃行走;动力学;微重力;碰撞

【作 者】龚宇鹏;谢天;林晓辉

【作者单位】哈尔滨工业大学卫星技术研究所,黑龙江哈尔滨150080;哈尔滨工业

大学卫星技术研究所,黑龙江哈尔滨150080;哈尔滨工业大学卫星技术研究所,黑龙

江哈尔滨150080

【正文语种】中 文

【中图分类】TP242.3;V448

0 引言

随着航天技术的飞速发展，各航天大国开展了一系列深空探测活动，并将目光投向

了太阳系中的小行星、彗星等小型天体。这些天体数量众多，对人类研究太阳系和

恒星的成因及演化历史，探索生命起源，获取太空资源具有重要作用[1]。在深空

探测中，小型天体本身的特殊形状导致了其引力场的特殊性，不能将它们单纯视作

中心引力场，因此在航天器对其接近和着陆的过程中，可以将小型天体视作多个引

力中心的组合，从而精确描述其引力特性[2]。

目前，只有美国“近地小行星交会”(NEAR)探测器、日本“隼鸟一号”探测器、

欧洲太空局“罗塞塔号”彗星探测器与小行星表面有过直接接触[3-5]。2000年2

月，NEAR探测器在设计任务之外完成了在“爱神”小行星上的软着陆，获得了宝

贵的探测资料和成果。“隼鸟一号”探测器于2005年11月在“丝川”小行星表

面软着陆并采集了土壤样品。2014年11月，由“罗塞塔号”彗星探测器释放的

“菲莱”探测器成功登陆“67P/楚留莫夫-格拉希门克”彗星，但由于探测器着陆

位置不良，导致其在很长时间内处于休眠状态。另外，美国“黎明号”探测器、日

本“隼鸟二号”、中国“嫦娥二号”探测器也已成功开展了对小行星的探测任务。

其中“嫦娥二号”探测器于2012年12月成功近距离飞跃“战神-图塔蒂斯”小行

星，实现了中国首次对小行星的近距离观测[6]，走出了中国深空探测极其重要的

一步。

研究小行星探测器运动及控制方式能为小行星探测方案设计提供依据，对深空探测

具有重要意义。目前国内外成功使用的行星表面探测器均为轮式车型探测器。虽然

轮式车型探测器在月球、火星等大质量天体上的使用获得了成功，但在小行星的微

重力场下，由于探测器车轮与地面之间的接触力非常小，使车轮打滑而移动缓慢，

且在小行星粗糙的表面，一个小扰动就会使探测器脱离地面，因此，轮式车型探测

器并不适用于小行星表面探测，需要寻求其他可行的方案。

美国、日本、欧洲各国和地区提出了一种新的行走方案——跳跃行走，并对其开

展了相关研究[7]。跳跃行走能很好地适应小行星表面的微重力环境。探测器从小

行星表面起跳后，在空中进行弹道运动并落在另一个位置，从而实现了行走的目的。

探测器在跳跃行走过程中能轻易越过几倍甚至几十倍于自身尺寸的障碍物，能适应

小行星复杂和不可预测的地形地貌，且探测器的运动范围大，能耗低。因此，跳跃

行走是小行星表面巡游探测比较可行的方案。

国内外针对探测器在小行星表面跳跃行走的实现开展了大量研究。如

HANAZAWA等[6-7]设计了一种辐条轮机构，可在松软的行星表面进行移动；柏

龙等[9]设计了用于小行星探测的跳跃机构，通过折叠机构将弹性势能转变为探测

器弹跳的动能，具有极强的越障和行走能力；GAJAMOHAN等[10-11]提出了利

用反作用飞轮进行跳跃的立方体模型，可通过飞轮制动，短时间内获得可使探测器

起跳的较大力矩，但其对起跳后探测器的姿态、跳跃过程的控制问题没有深入研究；

PAVONE等[12-13]设计了一款仿刺猬机器人，其通过刚性杆与地面接触，并通过

反作用飞轮实现对起跳和移动的控制，但未对机器人的移动策略进行进一步研究；

SHEN等[14-15]针对探测器在小行星表面的跳跃运动进行了研究，但未考虑探测

器在运动过程中的控制问题。

采用跳跃机构会增加小行星探测器质量及成本，容易降低探测器可靠性，在不增加

额外跳跃执行机构的情况下，基于反作用飞轮的小行星跳跃探测器具有结构简单、

可靠性强的优势。为此，本文针对轮控探测器建模及跳跃过程的控制进行研究。

1 跳跃移动原理

目前，对小行星跳跃探测器跳跃机理的研究主要分为两类：一是通过跳跃机构将挤

压力作用于行星表面，使探测器获得水平速度的弹跳跳跃方式；二是通过反作用飞

轮的加速改变探测器与地面的接触力，使探测器起跳。

1.1 弹跳跳跃原理

弹跳跳跃的灵感主要来源于将跳跃移动作为主要移动方式的跳蚤、青蛙、袋鼠等动

物。它们一般具有强壮的下肢，在跳跃时依靠下肢折叠产生弹力，并利用上肢调整

跳跃姿态，以实现跳跃移动[8-9]。

弹跳探测器的跳跃移动需要通过跳跃机构来执行，其原理为：跳跃机构通过自身折

叠储存弹性势能，积累后迅速释放，使地面支持力增加，同时产生前向的地面摩擦

力，使探测器跃起并向前移动。

图1 弹跳跳跃原理示意Fig.1 Principle of bounce hopping motion

根据不同的工程应用，弹跳探测器所采用的跳跃机构或运动方式各有不同，但总结

其运动特性，可以得知弹跳跳跃方式有以下几个特点：

1)具有极强的越障能力，可以越过相较于自身尺寸数十倍的障碍物，落脚面小，能

够适应复杂地形与未知环境。

2)能更好地适应小行星的微重力环境，弹跳行走过程能耗较小，并能够解决地面摩

擦力不足的问题。

3)弹跳探测器一般具有较轻的质量和较大的底面，可有效减小压强，在松软的行星

表面跳跃效果较好。

尽管弹跳跳跃能很好解决探测器在微重力环境中的移动问题，但由于引入了跳跃机

构，增加了探测器的结构质量，降低了可靠性，且难以在跳跃过程中对探测器姿态、

运动方向进行控制，无法保证其稳定性，容易发生倾翻，造成仪器及设备的损坏。

1.2 基于反作用飞轮的跳跃原理

基于反作用飞轮跳跃的起跳过程如图2所示，主要包括：

1)反作用飞轮从初始状态开始加速，随着飞轮力矩的逐渐增大，为平衡飞轮力矩，

地面等效支持力的作用点将逐渐向探测器起跳支点D处移动。当地面支持力完全

作用于起跳支点且重力力矩无法与飞轮力矩相平衡时，与倒立摆起摆时的原理相似，

探测器将绕起跳支点发生转动。

2)当探测器发生倾斜后，重力相对于支点处的力矩将随倾角的增加而减小，而飞轮

仍处于加速状态，因此探测器将绕其支点加速旋转。此时，探测器质心的水平速度

与竖直速度都将逐渐增加，直到探测器重力不足以作为以支点D为圆心的加速圆

周运动的向心力时，探测器跃起。

3)在探测器跃起后，由于其在旋转过程中已获得了一定的水平及竖直速度，探测器

将脱离小行星表面进行抛体运动。在探测器起跳后，飞轮可以逐渐减速，提供反作

用力矩，使探测器姿态恢复平稳。

图2 探测器起跳过程的二维模型Fig.2 2D model of probe takeoff process

探测器的成功起跳是在小行星表面进行跳跃探测的基础。为了研究起跳过程，需要

探讨能够使探测器产生翻滚或跳跃的飞轮力矩最小值，并且需要知道施加力矩与探

测器跳跃瞬间运动参数的关系。对探测器起跳的二维动力学模型进行分析。根据牛

顿第二定律和动量矩原理，得到飞轮起跳过程中相对于起跳支点D的动力学方程，

即

(1)

式中：ID为探测器相对于旋转支点D的转动惯量，ID=I+mL2；I为探测器绕质心

的转动惯量；α为探测器与地面接触的支点所对应的对角线与地面X轴负方向的

夹角；为坐标x，y及夹角α的二阶导数。相应的 T为飞轮转动时施加在探测器上

的力矩；m为探测器与飞轮的总质量；g为小行星表面的重力加速度；L为探测器

的质心到任意顶点的距离；ω为探测器绕支点D转动时的角速度。

假设探测器两顶点之间的夹角为2β，则由式(1)可得探测器所需最小起跳控制力矩

的公式，即

Tmin>mgLsin β

(2)

式(2)也是跳跃探测器能够起跳的条件。假设探测器模型为质量分布对称的均匀立

方体，则L为立方体半对角线的长度。由式(2)可知：L越大，则飞轮起跳所需的

控制力矩就越大；同时，探测器起跳所需的控制力矩也与探测器2个顶点之间的

夹角有关，当L不变时，该夹角越大，探测器起跳所需要的飞轮输出力矩也就越

大。

2 探测器动力学建模

2.1 探测器姿态动力学模型

将小行星跳跃探测器假设为1个刚性立方体(见图3)，探测器的质量均匀分布，在

质心处有三轴正交安装的反作用飞轮，通过飞轮加速及制动改变探测器的姿态，实

现跳跃移动。

图3 小行星跳跃探测器三维模型Fig.3 3D model of asteroid hopping probe

为明确反作用飞轮对于探测器姿态的控制作用，需要推导飞轮与探测器的姿态动力

学方程。根据欧拉方程可得

(3)

式中：下标“I”代表探测器在惯性坐标系下的时间导数；下标“B”代表探测器在

本体坐标系下的时间导数；H为探测器的角动量。

由欧拉方程可以推导得到在系统质心坐标系Os-xsyszs下带惯性轮的多刚体姿态

动力学方程，即

(4)

式中：I为系统关于其质心的惯性张量在系统质心坐标系下的矩阵；ω为探测器的

角速度矢量在系统质心坐标系下的矩阵；U为各动量轮组成的安装矩阵；Jw为各

动量轮转动惯量组成的对角阵；Ωw为各动量轮的转速组成的列阵；Tb为外力矩

在系统质心坐标系下的分量列阵，Tb=[Tx Ty Tz]T=rOsX×F，rOsX为小行星坐标

系下探测器质心到与地面碰撞点的向量,F为探测器所受的合外力。

探测器内部的3个飞轮是正交安装的，安装矩阵U为单位阵，则

(5)

(6)

式中：Mk为飞轮k的控制输入力矩，k=B，C，D。通过控制Mk可执行不同探

测器控制策略。

2.2 反作用飞轮模型

反作用飞轮的工作原理是动量矩定理，即星体的总动量矩矢量对时间的导数等于作

用于星体上外力矩矢量的总和。反作用飞轮通常作为执行机构对航天器进行姿态控

制，本文中反作用飞轮是探测器姿态控制及跳跃运动的唯一执行机构，即采用飞轮

加速既能使探测器具有运动速度，还可利用其进行姿态调整和控制。

反作用飞轮系统通常可简化为一阶系统进行控制，即

(7)

式中：Tw=RJ/(KeKm)为飞轮电动机的机电时间常数；R为电动机阻抗；J为反作

用飞轮转动惯量；Ke为电势常数；Km为力矩常数；Kw=1/Ke为飞轮电动机的增

益系数。

轮控系统的实际应用一般更关心反作用飞轮产生的实际力矩M对期望控制信号Tc

的响应，因此希望将控制器环境生成的期望控制力矩，作为飞轮的输入(控制)信号

进行输入。在飞轮电机的模型中加入电枢电流反馈，调节电机电枢电流，产生所需

要的电磁驱动力矩，它们之间的关系为

(8)

(9)

式中：Tc(t)为飞轮的期望控制力矩函数；i(t)为电枢电流；T′w=R/KP；KP为比例

系数。式(9)为力矩模式飞轮的基本数学模型，可看成1个惯性环节。对于性能较

好的惯性轮，惯性环节的时间常数可达到毫秒级。本文选取T′w为0.01 s。

本文选用的反作用飞轮性能指标见表1。

表1 反作用飞轮的性能指标Tab.1 Parameter of reaction wheel类别技术指标

转动惯量/(kg·m2)0.002最大转速/(r·min-1)2 000最大输出力矩/(N·m)0.01最大

角动量/(N·m·s)0.4粘性摩擦系数/(N·m·s)2×10-6

2.3 探测器与地面的接触力模型

由于在飞轮加速使探测器起跳的过程中需要计算碰撞时间和碰撞力等信息，而传统

的瞬间碰撞模型无法揭示这些信息，因此采用等效弹簧阻尼模型，将探测器与地面

之间的接触等效为1个弹簧阻尼模型。在等效弹簧阻尼模型中，法向受力由基于

Hertz定律和非线性阻尼函数的非线性弹簧-阻尼模型计算得到[16-17]。切向力由

简化的Karnopp摩擦力模型计算得到[18]。

探测器与地面的法向碰撞力用弹簧接触力来描述，由Hertz定律确定力的大小，

用阻尼系数表征接触过程的能量损失。基于Hertz定律和非线性阻尼函数的非线

性弹簧-阻尼模型表达式为

(10)

(11)

(12)

式中：k为等效接触刚度，单位为N/m2；Ri为两接触面的曲率半径；Ei，νi分别

为两物体的弹性模量、泊松比；c为阻尼系数；e为恢复系数；v0为撞击点法向初

始相对速度；δ为撞击点法向穿透深度；为撞击点法向相对速度。

为描述探测器的切向碰撞力，需要考虑到探测器在跳跃过程中可能发生的粘滞与滑

动两种运动状态，因此引入简化的Karnopp模型对探测器切向摩擦力进行描述。

Karnopp首先定义了一个速度转换容限VD，当切向速度小于VD时，将探测器状

态视为粘滞，否则为滑动状态。因此，针对简化后的切向碰撞力为

(13)

式中：Fs=μ0Fn为静摩擦力。

3 探测器控制过程及控制律设计

跳跃探测器在小行星表面的跳跃运动主要依靠其内部的飞轮来驱动，飞轮的控制力

矩是整个跳跃系统的控制输入。在开环条件下，探测器能在一定程度上适应地形变

化，但探测器的路径会与预定的路径有较大偏差，无法准确到达预定地点进行探测。

通过合理的控制算法控制飞轮的输入力矩，就能够控制飞轮的转速和角加速度，实

现对探测器跳跃的速度、方向等的控制，实现探测器对行星表面多个预定地点的准

确到达。在整个探测过程中，需要控制探测器的速度不能过大，防止由于速度超过

小行星逃逸速度而脱离小行星的引力场，也要注意飞轮的饱和及卸载问题。另外，

跳跃探测器从地面跃起后，由于其具有初始旋转速度，会在空中进行旋转和翻滚。

如果不对其进行控制，则探测器落地时随机的姿态会使下次跳跃变得十分复杂。因

此，在实际应用中期望每次都以固定的姿态及角速度落回地面，这就需要在每次探

测器跳跃的过程中通过飞轮的控制来调整其姿态，以控制探测器落地姿态。

探测器起跳或与地面发生碰撞后，均会产生一定的角速度。若探测器内部角速度具

有耦合时，探测器还可能具有角加速度，这会使探测器的姿态运动变得复杂，不利

于探测器的控制，因此主要在探测器滞空时对探测器姿态进行调整。

为了对探测器进行姿态稳定控制，需要引入探测器本体坐标系与目标姿态坐标系间

的误差四元数Qe=[qe0 qe]T∈R4×1，假设小行星本体系下的探测器当前姿态四

元数为Q，目标姿态四元数Qf=[qf0 qf]T，则误差四元数

Qe=Q-1Qf

(14)

式中：Qe满足的约束条件为

qe02+qe12+qe22+qe32=1

(15)

(16)

目标姿态四元数可根据不同的姿态任务确定，当前姿态四元数可由姿态测量系统测

量得到。假设探测器当前的姿态角速度为ω，目标角速度为ωf，均定义在航天器

本体坐标系中，则定义误差角速度

ωe=ω-ωf

(17)

将式(17)代入式(4)，可得到误差动力学与运动学方程，即

(ωe+ωf)×[I(ωe+ωf)+JwΩw]+

(18)

(19)

式中：

本节中探测器的姿态机动控制问题就是在给定目标姿态四元数下，设计控制律使得

当前姿态四元数在系统存在参数摄动及外界干扰的情况下，能在有限时间内

Qe→[±1 0 0 0]T，ωe→0。

为了达到更好的控制效果，考虑到飞轮惯性环节及克服摩擦力矩的干扰，本文在控

制环节采取的控制率为

(20)

式中：KL为控制器的回路增益矩阵；k1，k2，k3分别为比例、微分、积分控制

项。经调试，k1=0.02，k2=0.09，k3=0.000 1。

4 跳跃移动仿真校验

为验证探测器基于反作用飞轮进行跳跃移动方案的可行性，需要对其跳跃过程进行

仿真。假设探测器飞轮通过施加控制使探测器起跳后保持匀速前进，仿真条件设置

如下。

1)小行星表面环境及物理参数。根据现有的小行星环境，假设小行星表面重力加速

度g=0.001 m/s2，逃逸速度vescape=2.5 m/s，探测器与小行星表面摩擦系数

μ=0.45，静摩擦系数μ0=0.47，恢复系数e=0.95，刚度系数K=0.5 GPa。

2)探测器的基本参数。质量为5 kg，转动惯量I=diag([0.055 0.055

0.055])kg·m2，探测器边长2l=0.4 m。沿探测器3个飞轮转轴施加的最大飞轮力

矩大小Tc=0.01 N·m。

3)仿真初始条件。探测器初始角速度ω0=[0 0 0]T rad/s，反作用飞轮初始转速

Ω=[0 0 0]T rad/s，探测器初始姿态参数Q=[1 0 0 0]T。

图4为通过飞轮加速起跳后的探测器质心运动轨迹。由图可见，对反作用飞轮施

加控制，可实现探测器在小行星表面跳跃行走。行走过程仿真时间为800 s，得到

探测器的起跳时间为3.16 s，跳跃的最大高度为1.33 m。

图4 探测器质心运动轨迹Fig.4 Trajectory of probe centroid

探测器质心高度变化如图5所示。从图中可以清晰看出探测器每一次跳跃的过程。

探测器在起跳后达到了最高跳跃高度，继续跳跃移动后，其滞空时间与跳跃高度逐

渐减小，直至达到稳定的前向跳跃。

图5 探测器质心高度变化Fig.5 Variation of probe’s centroid height

图6 探测器质心速度变化Fig.6 Variation of probe’s centroid velocity

探测器各方向速度分量随时间的变化情况如图6所示。图中：水平方向速度分量

vx与vy的变化曲线形状为阶梯形。由于跳跃探测器在空中飞行时只受到法向重力

作用，切向受力为0，因此速度大小和方向都未发生改变。但当探测器与地面发生

碰撞时，在水平方向上受到了摩擦力的作用，使质心水平速度在短时间内发生了较

大变化；竖直方向速度分量vz的曲线为锯齿形，这是因为在空中飞行时探测器受

到的重力作用、与地面碰撞时受到竖直方向接触力的作用导致的。由于与地面碰撞

的时间很短，因此每次碰撞时质心速度和角速度均会出现突变。

由以上结果可看出，小行星探测器可以仅通过飞轮在小行星表面微重力的环境下实

现跳跃移动。对跳跃时的探测器姿态施加控制，可维持或改变跳跃方向，达到向不

同方向移动的目的。由于每次跳跃的距离及跳跃高度相对稳定，使用反作用飞轮与

使用弹跳机构的跳跃移动相比，可以更好地对探测器的飞行姿态及跳跃过程施加控

制，更适合于微重力环境下的小行星探测任务的执行。

5 结束语

本文针对小行星探测器跳跃行走方案的动力学和控制策略开展研究。分析了国内外

小行星跳跃探测器的发展现状，建立了跳跃探测器的动力学模型，给出了反作用飞

轮在探测器跳跃移动中所需提供的最小力矩，设计了探测器姿态控制律，进行了仿

真检验。结果表明：探测器能进行完整的姿态控制过程，可以执行小行星表面的探

测任务。对跳跃探测器跳跃行走过程进行了仿真，探测器通过反作用飞轮加速提供

的力矩产生跳跃运动，验证了方案的可行性。考虑到小行星表面的地形变化及非均

匀的重力场环境，未来将研究跳跃探测器在不同地形环境及非均匀重力场下的跳跃

运动，针对不同地形下所采取的跳跃移动策略及相应的控制方法进行设计，以满足

不同条件下小行星着陆探测任务的需求。

参考文献

【相关文献】

[1] 杨剑峰, 盛英华, 徐博. 我国小行星探测发展趋势研究[C]//中国宇航学会深空探测技术专业委员

会学术年会. 2012: 1196-1207.

[2] 刘付成, 朱东方, 黄静. 空间飞行器动力学与控制研究综述[J]. 上海航天, 2017，34(2): 1-26.

[3] VEVERKA J, ROBINSON M, THOMAS P, et al. NEAR at Eros: imaging and spectral

results[J]. Science, 2000, 289(5487): 2088.

[4] ULAMEC S, BIELE J, BOUSQUET P W, et al. Landing on small bodies: from the Rotta

lander to MASCOT and beyond[J]. Acta Astronautica, 2014, 93(2013-03): 460-466.

[5] TARDIVEL S, CANALIAS E, DELEUZE M, et al, Landing MASCOT on asteroid 1999 JU3:

solutions for deploying nanosats to small body surfaces[C]// 46th Lunar and Planetary

Science Conference, 2015， 223-226.

[6] 叶培建, 黄江川，张廷新，等. 嫦娥二号卫星技术成就与中国深空探测展望[J]. 中国科学：技术

科学, 2013(5): 467-477.

[7] HANAZAWA Y, ASANO F. High-speed biped walking using swinging-arms bad on

principle of up-and-down wobbling mass[C]// IEEE International Conference on Robotics

and Automation. 2015: 5191-5196.

[8] XIAO X, KIKUCHI Y, ASANO F, et al. Limit cycle walking of underactuated bipedal

humanoid on slippery road surface[C]// IEEE-Ras International Conference on Humanoid

Robots. 2014: 622-627.

[9] 柏龙, 葛文杰, 陈晓红,等. 用于行星探测的跳跃机器人研究[J]. 机器人, 2009, 31(4): 311-319.

[10] GAJAMOHAN M, MERZ M, THOMMEN I, et al. The Cubli: a cube that can jump up and

balance[C]// IEEE/RSJ International Conference on Intelligent Robots and Systems. 2012:

3722-3727.

[11] GAJAMOHAN M, MUEHLEBACH M, WIDMER T, et al. The Cubli: a reaction wheel

bad 3D inverted pendulum[C]// Control Conference. 2013: 268-274.

[12] PAVONE M, CASTILLOROGEZ J C, NESNAS I A D, et al. Spacecraft/rover hybrids for

the exploration of small Solar System bodies[C]// Aerospace Conference. 2013: 1-11.

[13] ALLEN R, PAVONE M, MCQUIN C, et al. Internally-actuated rovers for all-access

surface mobility: theory and experimentation[C]// IEEE International Conference on

Robotics and Automation. 2013: 5481-5488.

[14] SHEN H X, ZHANG T J, LI Z, et al. Multiple-hopping trajectories near a rotating

asteroid[J]. Astrophysics & Space Science, 2017, 362(3): 45.

[15] BELLEROSE J, SCHEERES D. Dynamics and control for surface exploration of small

bodies [C]//Proceedings of AIAA/AAS 2008 Astro dynamics Specialist Conference. [S.l.]：

AIAA， 2008： 18-21.

[16] HOCKMAN B, FRICK A, NESNAS I A D, et al. Design, control, and experimentation of

internally- actuated rovers for the exploration of low-gravity planetary bodies[M]// Field

and Service Robotics. Springer International Publishing, 2016.

[17] YIGIT A S, CHRISTOFOROU A P, MAJEED M A. A nonlinear visco-elastoplastic impact

model and the coefficient of restitution[J]. Nonlinear Dynamics, 2011, 66(4): 509-521.

[18] 贺波勇,顾绍景,黄海兵,等.中秋节载人登月任务窗口与转移轨道设计研究[J].上海航天, 2017,

34(5): 9-15.