2014年全国数学建模论文

承 诺 书

我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的报名参赛队号为（8位数字组成的编号）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)：

（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以

上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。）

日期： 年 月 日 赛区评阅编号（由赛区组委会评阅前进行编号）：

2014高教社杯全国大学生数学建模竞赛

编 号 专 用 页

评 阅 人 评 分 备 注 赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

A题 嫦娥三号软着陆轨道设计与控制策略

摘 要

月球软着陆轨道研究是我国开展登月活动所必须先期解决的关键问题之一。针对这个问题，本文提出了基于景象匹配技术并采用分段控制来实现对月面某已知点进行精确软着陆的方法[1]。

就问题一而言，我们采用二体模型对着陆准备轨道近月点和远月点的位置进行了确定，利用直角坐标之间的转换关系来求解出嫦娥三号在该位置的相应速度的大小和方向，不足之处是：忽略的影响因素比较多，使模型中的误差较大。

根据问题二的要求，以消耗燃料最低为目的建立了一个模型，计算出着陆轨道，并针对着陆轨道的不同阶段分别采用分段控制的方法，在着陆准备轨道、主减速、快速调整、粗避障段飞行、精避障段飞行、缓速下降月面等过程进行了研究，得到了着陆轨道的最优控制策略，此模型的缺点是缺少仿真分析进行对比检验，验证模型的正确性和结果的最优性，需要改进的是：对自己分段控制最优策略加以验证。

在问题三的解决中，我们利用matlab7.0软件绘出嫦娥三号的轨迹曲线以及相对应的高度、速度等的变化曲线，然后找出误差源，利用误差源中主要因素的相关数据进行分析，判断影响误差的因素及其严重性。缺憾的是：对于敏感性分析，考虑到嫦娥三号的初始质量、推力、比冲等多个不确定因素的影响，但是没有运用Excel的模拟运算表，对嫦娥三号安全着陆的概率进行敏感性分析。

关键字：景象匹配；分段控制；燃料最优；精确软着陆

1

一、 问题重述

嫦娥三号准备在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆轨道，着陆轨道为近月点至着陆点，其软着陆过程为6个阶段，要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。需解决如下问题：

问题一：确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小和方向。

问题二：确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。

问题三：对于我们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。

二、问题分析

嫦娥三号发射成功开启了我国航天事业的新篇章，它的软着陆首次实现我国航天器在地外天体巡视探测，首次实现对月面探测器的遥操作[2]。在嫦娥三号高速飞行的情况下，确保它可以在月球预定区域内实现软着陆，我们先使用简单的二体模型确定着陆轨道近月点和远月点的位置并求解出对应的速度大小和方向。接着采用最优控制燃料的模型确定着陆轨迹，运用不同的控制法和动力方程建立模型，求解出6个阶段的最优控制策略。在第三个问题的解答中，运用matlab软件做出不同的误差影响因素下嫦娥三号着陆的高度和速度变化曲线，分析出误差源中的重要影响因素。

三、模型假设

1、假设嫦娥三号可以实现月球表面的精确软着陆，探测器准确地降落在月面预定的区域；

2、假设嫦娥三号软着陆只考虑发动机的推力和发射角度，其他因素都可以忽略； 3、假设嫦娥三号着陆准备轨道不考虑质量的损失； 4、月球引力均匀，忽略月球自传；

5、嫦娥三号的所有电子设备均处于正常工作状态。

四、符号说明

2

符号 意义 初始轨道半径 月球半径 轨道高度 初始横向速度 单位 r0 R h vx0 vy 0 初始径向速度 嫦娥三号在远月点的月心距 推力方向与当地水平线的夹角 月球引力常数 行器发动机的推力 ra   F

/ 五、模型的建立与求解

5.1 问题1的模型建立与求解

考虑从15公里左右的轨道高度下降到接近月面这一阶段的飞行器轨道控制方法。由于月球表面附近没有大气，所以在飞行器的动力学模型中没有大气阻力项。而且从15公里左右的轨道高度软着陆到非常接近月球表面的时间比较短，一般在几百秒的范围内，所以诸如月球引力非球项、日地引力摄动等影响均可忽略不计。因此，使用较为简单的二体模型就可以很好的描述这一问题。

以月球飞行旗某时刻的星下点为原点的直角坐标系形式受控飞行器动力学方程为

vx vyFcos mr3Fsin(yR) mr3 xvx yvy

3

mQ

式中

—月球引力常数，=0.4902802627/;

F飞行器发动机的推力，其赋值一般认为是常值T，和分别是可提供选择的发动机

推力赋值的上下限。

Q－飞行器发动机的燃料秒消耗率；

－推力方向角（操纵角），即推力方向与当地水平线的夹角；

R－月球半径；

r－嫦娥三号月心距，r= 飞行器的初始条件

vx0vpr0Rh

2ra()r0r0ra 其中初始轨道半径r0为开普勒轨道近月点月心距，R为月球半径，h为轨道高度；初始横向速度vx并非当地环绕速度，而是在开普勒轨道运动的飞行器从位于较高停泊轨道的

0远月点运动到近月点的速度，这一速度大于当地环绕速度，是飞行器在远月点的月心距；由于初始时刻飞行器在近月点，所以初始径向速度vy=0.

0终端约束条件为

vx(tf)0 ；vy(tf)0r(tf)R

；

式中tf软着陆完毕时刻其物理意义是飞行器降落到月球表面时，速度为零。 经过求解得到，vpvx01.7km/s

根据已知材料得到嫦娥三号着陆过程中经度基本不变，纬度改变，月球赤纬和地球纬度一样，也分为南北各90个分度，又因为月球极区半径R=1735.843km,所以每个纬度的竖直高度差为19.28714km。

即近月点位置为（19.51W，67.55N，15km），远月点位置为（160.49E，22.45S，100km）

5.2 问题2的模型建立与求解

问题2中，确定嫦娥三号软着陆轨迹，最佳轨道即是达到燃料最省的性能指标，此时发

4

动机的比冲较大，推力会取为最大值，使之在尽量短的时间内完成软着陆轨道任务[3]。需要求解最优控制，对于最优控制问题，需要引入如下变量：

共轭变量：

=[r,v,,m,]T

Tf状态变量： x=[r m]T 构造哈密顿函数为 H(x,,u)= H=T(x,u)

f=rf1vf2f3mf4f5

L(yRL)LFF) =r(cos3x)V(sin3mrmr(rcosvsin)m]+H2 +Tmax=F[mIspg =H1(,F)+H2

在上式中，H2是与控制量无关的部分。根据Pontryagin极大值原理，最优控制使H最大，也就是H1最大，最优控制规律为：

F

(rcosvsin)m式中，函数s(t)=

mIspg可得到控制方程：

H0 最优推力控制角有 * =-

关于解决软着陆六个阶段的最优控制策略，采用了不同的控制方法和动力学模型，进行分段研究，控制方法主要包括控制策略（指对推力的大小，方向和发动机开机时间进行控制的方法。由于发动机的推力由发动机自身性能决定，于是本文主要针对给定推力大小的发动机的推力方向和开关机点进行研究）和对控制过程计算方法的研究。以下为各个过程的详细介绍。如图所示

5

图5.2-1

5.2.1 着陆准备轨道(H0-H1)

着陆准备段的主要任务是在动力下降点火前，修正地面注入的点火时刻机对应轨道，计算点火目标姿态并调整到位[3]。此时嫦娥三号应该在近月点做好软着陆的准备。忽略月球形状近视为球体，以月球中心为坐标建立平面直角坐标系，由近月点和远月点可得椭圆方程

x2y21 22(ar)(br)约束条件： a=15km b=110km r=1737.64km

嫦娥三号在次轨道进行软着陆准备。

5.2.2 主减速过程（H1-H2）

1、 主减速段控制策略 推力方向控制：

6

如图5.2-1所示，在主减速阶段嫦娥三号的速度方向可分解为横向速度v和径向速度

vr，在减速阶段主要的任务是减小嫦娥三号的横向速度为其后的粗避障飞行做准备。根据

此任务要求设计如下的控制过程：

嫦娥三号在近月点H1的初始轨道高度进行降轨,此时发动机的推力方向与嫦娥三号的速度相反,近月点轨道速度为1.7km/s，速度方向与当地水平面平行。此时嫦娥三号纵轴应与水平面平行，即安装在嫦娥三号纵轴方向的发动机喷口与当地水平面平行，推力方向与嫦娥三号横向速度方向相反（图5.2-1状态）所示。

主减速阶段的任务是减小嫦娥三号横向速度即嫦娥三号相对地平面的速度，为尽量减小水平方向速度并尽量降低控制的复杂程度，提高可靠性。决定采用定推力方向的控制策略，即要求推力方向始终保持在嫦娥三号纵轴的方向，并与嫦娥三号的飞行方向相反。在嫦娥三号减速初期推力方向与当地水平面的方向平行并指向速度的相反方向（推力方向角）。但由于月球是近似的球体，所以月球表面存在弧度。嫦娥三号在横向飞行阶段纵轴会与月面产生夹角（与当地水平面夹角图5.2-1状态所示），此时嫦娥三号的推力可分解为一个与当地水平面平行的分力用于减小水平速度，一个垂直方向的分力将增加垂直方向速度（由于水平飞行距离相对于月球半径是小量所以发动机推力方向相对于当地水平面的偏角极小即垂直分力极小）。

约束条件（H2点的约束条件）：

根据任务需要，在横向减速完成时（到达H2点），嫦娥三号的横向距离将减小到零即

v=0，所以终端约束条件为v=0，其飞行过程如图5.2-1.

2、 控制过程计算方法：首先将H1点初始状态参量带入到横向减速动力学方程中

v2vr2asin

rrvrvvacos

rrvr

7

v

r由于求解微分方程，所以计算机将对参数进行迭代计算得到下一个迭代时间点的各个状态参量。并以此量为初始参量返回带入动力学模型中进行下一次迭代计算，直到满足嫦娥三号终端约束条件，计算停止并输出各个状态值作为下一阶段控制过程的初始值。

计算初值:

（1）初始条件：v=1.7km/s

vt=0

h=15km

m0=2400kg

=0

（2）约束条件： v=0

将以上条件带入方程，可得到嫦娥三号主发动机的工作时间（即减速发动机关机点时刻）并输出H2点的各个状态值作为下一阶段计算的初始值。控制过程流程如图5.2-1所示。下图5.2-1中所示的X轴位嫦娥三号纵轴，Y轴位嫦娥三号横轴，为嫦娥三号纵轴与当地水平面的夹角。

5.2.3 快速调整阶段 (H2-H3)

在H2-H3段，以一个阶段的终点高度H2为起点，以满足相机开机要求的高度H3为终点，选择一个适当的发动机开启时刻以满足H3点降落相机时机速度要求。

1、快速调整阶段控制策略 推力方向控制：

由上一个飞行阶段终点H2作为起点，根据上一个控制过程产生的H1点状态值，即有嫦娥三号纵横相对当地水平面的偏角得到姿态调整信息，对嫦娥三号的进行姿态调整（如图

5.2-1 状态）；是嫦娥三号纵轴沿垂直平面旋转，最终将主发动机的喷口方向

28

（嫦娥三号纵轴）垂直当地月表面（如图5.2-1状态），选择主发动机的开机点（如图5.2-1 状态）利用主发动机推力抵消径向速度。 约束条件（H3点的约束条件）：

H3点的约束条件：由嫦娥三号携带的降落相机的开机高度和速度要求为基准，得到嫦娥三号进行径向减速后所需的径向速度和距离月面的高度。

vy=20m/s

vx=0

h=2.4km

2、控制过程的计算方法：

首先将下落时间t计入径向减速动力学方程，设自由落体的时间为（主发动机开机点其中包括姿态调整时间），主发动机开机时间为，可得公式

t=+

(5-1)

设初始速度为，自由下落末速度（开机点速度）为，终点速度为。自由下落距离为，开机减速距离为，下落总距离为S。可将式5-1代入动力学方程

Tvyg (5-2)

mxvx yvy

mQ

可得到重力场影响下径向速度公式

vygvy0gt （5-3）

vygvy1t20Fdtgt2 （5-4） mmt方程5-2可得到重力场影响下径向位移公式

12 （5-5） S1vy0gtdtvy0t1gt202t29

S2t20F(vyggtdt)dtmmt(5-6a)

mmmm

ln(1mt)dtln(1mt)ttm(1mt) （5-6b）

代入式5-4可得

12Fmmm S2vygt2gt2+[ln(1t2)t2t2ln(1t2)]dt （5 -6c）

2mmmm

112FmmmSvy0t1gt12+gt2[ln(1t2)t2t2ln(1t2)]22mmmm

vygvy0gt1 （5-7）

vygvy1t20FFmdtgt2=vy1[ln(1t2)]gt

mmtmmg 由式中5-7可知，和vy为未知项，可由已知条件vy及约束条件S和vy代入到方程中解

01得，其中时间需大于等于姿态调整所需要的时间。具体计算过程为：

将H2点初始状态参量和嫦娥三号姿态调整所需要的时间代入到径向减速动力学方程5-7中，对方程加入H3点的终端约束，方程所得未知数为主发动机开机时间和开机高度H，由于方程有解析解，未知数可求。从而得到了径向减速段的控制过程。即嫦娥三号由H2高度开始进行自由下落，下落期间完成嫦娥三号的姿态调整过程，是主发动机喷口垂直朝向月面。其后经过时间t，嫦娥三号到达高度H，主发动机开机进行径向减速，最终在H3点满足降落相机的开机条件，径向减速控制过程结束。输出H3点各个状态值，作为其后的粗避障嫦娥三号的初值。

5.2.4 粗避障飞行过程（H3-H4）

以上一个阶段降落相机开启点H3为起点，再导致信息形成并完成对准目标点的调姿工作后，横向辅助发动机起动飞向目标点，经向减速发动机同时开启以减小垂直方向上的速度，到达目标点上空后，主发动机起动对垂直方向进行主减速最终保证嫦娥三号到达距100m的H4点时满足对月速度要求此时粗避障飞行结束。 1、粗避障飞行控制策略 推力方向控制：

10

由于上一个控制段结束时嫦娥三号各状态值已经达到降落相机开机条件，此时嫦娥三号所拍摄的图像和机器内存储的环月阶段所拍摄的基准图像进行匹配计算，形成目标点位置信息。由此信息可得到辅助发动机喷口方向与嫦娥三号和目标的角度偏差和距离r。

以此信息作为粗避障飞行条件，首先进行嫦娥三号的水平姿态调整，使嫦娥三号沿本体纵轴旋转角，以使辅助发动机的喷口线方向朝向嫦娥三号和目标连线方向。其后嫦娥三号在水平面将沿直线飞向目标点，在水平飞行过程中，背向目标点的喷口先起动，加速向目标点飞行。接着在点进行喷口切换，背向目标点喷口关闭，面向目标点的发动机喷口开启，减速飞翔目标点并最终在目标点上空横向减速为0.此时横向制导飞行段结束。

在横向粗避障飞行段辅助发动机工作的同时，主发动机也开启工作，由于采用推力阶梯可调式发动机，比冲一定时，推力越小的发动机燃料消耗越少，可选择较小的推力进行径向减速以保证最终满足约束条件。

横向粗避障飞行段结束后，由于嫦娥三号的横向的速度减为零，为了在H4点满足约束条件，嫦娥三号将再进入一个径向减速阶段。采用先自由下落后主发动机开启减速的过程，知道在点满足约束条件，粗避障飞行段结束。 约束条件：

H4点的约束条件由嫦娥三号精度要求和飞行仍无要求为准，得到嫦娥三号进行垂直减速后所要达到的径向速度以及距离月表面的高度。

vy=20m/s

vx=0

r=100m

2、控制过程的计算方法

首先将粗避障飞行段分为横向粗避障段和径向减速段两部分进行计算，与径向减速过程方程建立方法相似。设横向粗避障飞行总时间为t，自由下降时间为，从主发动机和辅助发动机开季到辅助发动机喷口切换点的时间为，从喷口切换点到横向粗避障段结束的时间为，可得方程

tt1tHxt2 （5--8）

（1）粗避障段径向运动方程建立：由于初始横向速度为0，设初始径向速度为vy，自

0由下落末速度为vyg，横向粗避障终点径向速度为vy.自由下落距离为Sy，开机减速距离

1111

为Sy，下落总距离为Sy。由于主发动机开机减速时间为tH2xt2，由5-3，5-4可得垂直

方向速度公式

vygvy0gt1 （5—9）

vygvy1tHxt20Fdtg(tHxt2)mmtFmln(1(tHxt2))g(tHxt2)mm (5—10)

由径向减速方程式可得：

Syvy0t1(tHx121Fmgt1+vyg(tHxt2)g(tHxt2)2[ln(1(tHxt2))22mmmmt2)(tHxt2)ln(1(tHxt2))]mm(5-11)

vygvy0gt1

vygvy1式中

tHxt20FFmdtg(tHxt2)vy1ln(1(tHxt2))g(tHxt2)mmtmmSy--------垂直下落总位移；

vy0-------点初始速度

vyg自由落体末速度；

vy1点的终端速度；

自由下落时间（发动机开机点）；

tHxt2----主发动机开机时间；

由式5-11可知，tHxt2，和S为未知项，可由已知条件和横向飞行时间代入方程中

解，其中时间需大于等于姿态调整状态所需的时间，最后解得位移Sy项作为其后的径向减速段初始条件。

（2）横向粗避障段横向运动方程的建立：横向运动与径向减速段类似可得到横向速度和位移公式，设辅助发动机开机点到喷口切换点的横向位移为SH,喷口切换点到横向粗避障

x12

段结束（横向速度减为0）的横向位置Sx,喷口切换点横向速度为vH，横向总位移为

2x，喷口切换点嫦娥三号质量为mHmmtH代SxSx2SHx（由粗避障信息给出）

xx入动力学模型

TyTx，则横向速vx vyg xvx yvy mvQmmFFmdtln(1mmtmm度公式可表示为

vHxt20Hxt ) （5-12）

vHx横向位移公式为

t20FFmdtvHxln(1t2)0 （5-13）

mHxmtmmHxFmmmSHx[ln(1tHx)tHxtHxln(1tHx)]mmmmSx2vHxt2 联立方程可得

(5-14)

Fmmm[ln(1t2)t2t2ln(1t2)] (5-15) mmHxmmHxSxvHxt2

Fmm[ln(1t2t)t22mmHxmmln(1tmHx)]2)]



Fmm[ln(1tHxtH)tHxxmmmFmln(1mm)HtxFmmln(1tHxmmln(1mHx2vHx t ) （5---16）

mHxmmtHx

式中 --------反向喷口工作时间

--------正向喷口工作时间

由式5---16可知，，vH为未知项，可由已知条件及约束条件带入到方程中解得，tH为

xx喷口切换时间，tH，t2可作为径向减速初始条件。

x（3）粗避障飞行结束后的径向减速飞行方程建立：由于径向减速过程与H2-H3的径向减速过程基本相同，由5—7可得到主发动机的开机点和开机时间.

13

将H3点初始状态参量和嫦娥三号姿态调整及景象匹配计算时间约束代入到横向粗避障飞行段方程5—16中，并加入粗避障信息产生的横向惟一约束条件1，解得的方程未知数为辅助发动机喷口切换时间点和辅助发动机开机总时间tHxt2。

由径向运动方程可得到开机高度，方程解得的未知数又可作为横向粗避障飞行结束后的径向减速段初始值代入径向减速动力学方程5—7.代入H4点的约束条件2解得径向减速的开机点和开机时间。从而得到了径向段的减速段的控制过程。即嫦娥三号由横向减速终点高度开始进行自由下落，经过时间，主发动机开机进行径向减速，最终在H4满足任务要求。输出H4点各状态值，作为其后垂直降落月面的初值。

5.2.5 精避障飞行过程（H4-H5）

以上一个阶段精避障结束点H4为起点，在精避障信息形成并完成对准目标的调资工作后，减速发动机同事开启以减小垂直方向上的速度，到达目标点上空后，主发动机启动对垂直方向进行主减速最终保证嫦娥三号到达距30m的H5点时满足对月速度要求此时精避障飞行结束。 约束条件：

H5点的约束条件由嫦娥三号精度要求和飞行任务要求为准，得到嫦娥三号进行垂直减速后所要达到的径向速度以及距离月球表面的高度。

vy20m/svx0

r=30

控制过程的计算方法:

精避障的控制分过程与粗避障径向减速过程相似，所以方程建立方法相似。 径向减速方程为：

11Syvy0t1gt12vyg(tHxt2)g(tHxt2)222

Fmmm[ln(1(tHxt2))(tHxt2)(tHxt2)ln(1(tHxt2))]mmmmvygvy0gt1

vygvy1

tHxt20FFmdtg(tHxt2)vy1ln(1(tHxt2))g(tHxt2)mmtmm14

式中

Sy-------垂直下落总位移；

vy0------点初始速度

vyg自由落体末速度；

vy1 点的终端速度；

； t1自由下落时间（发动机开机点）

tHxt2主发动机开机时间；

代入H5点的约束条件解得径向减速的开机点和开机时间。从而得到了径向段的减速段的控制过程。即嫦娥三号由横向减速终点高度开始进行自由下落，经过时间，主发动机开机进行径向减速，最终在H5满足任务要求。

精避障过程结束的H4点开始，首先是嫦娥三号的自由下落阶段，主发动机选择合适的开机点以保证到达月球表面H5时相对月球的速度为0。

5.2.6 缓速下降阶段

1、垂直降落过程的控制策略

垂直降落过程与镜像减速过程基本一致，但是少了嫦娥三号进行姿态调整的过程，控制策略更加简化。

约束条件（H5点的约束条件）：

H5点的约束条件即软着陆的基本要求为已接近0的速度降落在月球表面，即

vy=0

vx=0

r=

2、控制过程的计算方法

由于径向减速过程与H2-H3的径向减速过程基本相同（只缺省了姿态调整的约束，不影响方程的计算），可将方程5-11中由（抵达月面速度）代替，可得方程

112Svy0t1gt12vygt21gt222 (5-17)

Fmmm[ln(1t2)t2t2ln(1t2)]mmmm15

vygvy0gt1

vygv月+

t20FFmdtgt2v月ln(1t2)gt2

mmtmmS垂直下落总位移；

vy0H3点的初始速度；

vyg自由落体末速度；

v月月球软着陆速度；

自由下落时间； 主发动机开机时间；

与径向减速段控制方法相同，将H4点初始状态参量代入到径向减速动力学方程5-17中，并代入H5点的终端约束，方程所得未知数为主发动机开机时间开机点，由于方程有解析解，从而得到了径向减速段的控制方法。

5.3 问题3 着陆轨道和控制策略的误差分析和敏感性分析

1、误差分析

理想情况下不存在各种系统误差，采用的各种制导过程，各个着陆指标均能达到设计要求。我们利用matlab7.0软件分析做出近月点到月面的高度曲线图、各个阶段速度变化的曲线图、质量变化曲线图（具体编程过程见附件），图形如下所示5-3-1。

14121086420h/m050100150t/s200250300 16

高度随时间变化曲线

1800160014001200V/(m/s)10008006004002000050100150t/s200250300 速度随时间变化曲线

图5-3-1

软着陆飞行轨迹的影响因素有很多，误差源主要包括导航设备测量误差和执行机构误差2部分，此外还有，月球非球形引力摄动和日、地引力摄动等环境干扰引起的的误差。其中测量测量误差包括地面测轨误差、惯性装置（IMU）测量误差、雷达高度计以及多普勒测速雷达误差；执行机构误差主要包括初始质量的估计误差和推力误差[4]。 （1）初始质量估计误差

初始估计的质量增加3%的误差，高度、速度变化图如下5-3-2所示。由于无法确定质量偏差带来的航程偏差，可以移动初始质量估计正偏差上界的带来的航程偏差，再引入推力开关函数来调整航程，进而能够达到招录区域[5]。

17

高度随时间变化曲线1

速度随时间变化曲线1

图5-3-2

（2）推力大小误差

推力增加3%的误差，高度和速度的变化形式就有所变动，这样看来，推力偏差影响比较大，这样就需要对制动发动机的推力提出更高的相对精度要求。高度变化和速度变化如下图5-3-3。

18

高度随时间变化曲线2

速度随时间变化曲线2

2、敏感性分析

敏感性分析是系统分析中分析系统稳定性的一种方法，对变动参数分析的一种不可或缺的工具，可以用Excel中的模拟运算表来实现，可分为单变量模拟运算和双变量模拟运算。我们分析得到不确定因素有：推力大小、比冲、初始质量等，经过与理想状态下的曲线相比，可以看出，嫦娥三号的轨迹会受以上因素的影响而发生明显改变，也进一步说明误差分析的重要性。

19

六、结果分析

通过简单的二体模型确定了着陆轨道上的近月点和远月点，从而我们可以解决嫦娥三号在这两点的运动情况，接着为了达到消耗燃料最低的目的建立轨道模型，并采用分段控制的方法和动力学方程得到嫦娥三号着陆轨道六个阶段的最优控制策略，最后我们用matlab7.0软件做出嫦娥三号的高度变化和速度变化的曲线，进行了简单的数据的误差分析，是我们的模型更加完整。

七、模型评价与改进

（1）本文所建模型的优点在于以燃料消耗最优性为出发点，设计分段控制和动力方程建立模型并求解，无需迭代计算。

（2）本文对软着陆过程中的六个阶段分别进行了分析和建模求解，体现出了阶段最优控制策略。

（3）如果可以进行仿真计算，模型的建立和求解将会更加完善，充分体现最优控制的软着陆导律的时效性和优越性。

参考文献：

[1]刘畅.基于景象匹配的月球探测器精确软着陆轨道的研究.哈尔滨工业大学硕士论文.2006:26-43

[2]新华网.嫦娥三号成功发射的意义:中国科技发展又一成就.http://chuzhong.eol.cn/focus_9144/20131202/t20131202_1047065.shtml .2014.9.13. [3] 孙泽洲,张廷新,张熇,贾阳,张洪华,陈建新,吴学英,申振荣等. 嫦娥三号探测器的技术设计与成就[J]. 中国科学:技术科学,2014,04:331-343

[4]单永正，段广仁，吕世良等.月球探测器软着陆的最优控制[J],2009,09:2153-2158 [5]关轶峰,王大轶,黄翔宇. 月球软着陆主制动过程的数学仿真分析[A]. 中国宇航学会深空探测技术专业委员会、国防科学技术工业委员会科技与质量司.中国宇航学会深空探测技术专业委员会第三届学术会议论文集[C].中国宇航学会深空探测技术专业委员会、国防科学技术工业委员会科技与质量司:,2006:7.

[6] 单永正.月球探测器软着陆的制导问题研究[D].哈尔滨工业大学 2009

20

附件

MATLAB程序： 子程序

a=imread(‘附件4 距月面100m处的数字高程图.tif’); A=im2double(a); xa=1:1000; xb=1:1000;

[xx,yy]=meshgrid(xa,xb); mesh(xx,yy,A); figure(gcf);

21

子程序：

a=imread('附件3 距2400m处的数字高程图.tif'); b=imread('附件4 距月面100m处的数字高程图.tif'); [x,y]=size(a); [X,Y]=meshgrid(1:x,1:y); Z=double(a); mesh(X,Y,Z); colormap

MATLAB程序2：

高度变化曲线的MATLAB程序： a=xlsread('C:\\MATLAB7\\work\\Book1.xls'); n=size(a,1); x=a(:,1); y=a(:,2); plot(x,y,'rs-')

set(gca,'xtick',[0:50:300],'ytick',[0:2:15]) grid on

速度变化曲线的MATLAB程序：

a=xlsread('C:\\MATLAB7\\work\\Book1.xls'); n=size(a,1); x=a(:,1); y=a(:,4); plot(x,y,'rs-')

22

set(gca,'xtick',[0:50:300],'ytick',[0:200:1800]) grid on

时间/s

0 13 26 39 52 高度/km 极角

15 14 13 12.1 11.3 0 1.3 2.1 2.7 3.5 速度 推力攻m/s 角

1700 182 1650 178.1 1600 176.9 1520 177.2 1480 176.8 65 10.5 78 9.8 91 8.8 104 8.1 117 7.4 130 6.4 143 5.8 156 5 169 4 182 3.3 195 2.5 208 1.9 221 1 234 0.5 247 0.2 270

0

4 1400 4.6 1300 5.1 1220 5.6 1180 6 1100 6.5 1000 6.9 900 7.3 800 7.5 710 7.8 610 8 540 8.1 420 8.15 350 8.2 220 8.23 130 8.3

10 23

176.2 176 175.5 175 174.5 174.1 173.9 173.6 173 172.5 172.3 172.1 172.2 172.8 174.8 167