2014高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
嫦娥三号软着路轨道设计与控制策略
摘要
本文主要为分阶段研究嫦娥三号的软着陆轨道设计与最优控制策略。 建立模型一确定近月点和远月点的位置，以及嫦娥三号速度大小与方向。首先以月球中心为坐标原点建立空间坐标系，根据计算的作用力可知地球影响较小，故忽略不计。然后将嫦娥三号软着陆看作抛物线的运动过程，计算在最大推力下的减速运动，求得月面偏移距离为，由此计算出偏移角度为15.25°。从而得出近月点和远月点的经纬度分别为（34.76°W，44.12°N）和（34.76°E，44.12°S）。最后在软着陆的椭圆轨道上，由动力势能和重力势能的变化，计算出嫦娥三号在远月点和近月点的速度分别为，沿轨道切线方向。
建立模型二和模型三确定着陆轨道和在6个阶段的最优控制策略。模型二主要对主减速阶段和快速调整阶段进行初步分析。模型三分六个阶段确定轨道和最优控制策略，主减速阶段建立目标函数燃料，假设推力最大,将最优燃耗软着陆问题转化为最短时间控制问题，然后采用拟牛顿法和四阶Admas预测-校正得到；快速调整阶段采用重力转弯制导，在假设条件下对嫦娥三号进行受力分析，得到嫦娥三号的动力学模型，然后通过开关控制得到燃耗最优控制，并画出仿真图；粗避障阶段采用多项式制导，通过初始状态和末端状态反解多项式系数进而求取标称轨迹，然后将避障区域网格化，比较网格内的方差大小确定最优区域范围；精避障阶段需在满足本文提出的避障原则式下搜索全局最优解，以网格区域总体得分作为目标函数，得到最优区域为坐标附近，并以螺旋搜索法搜索安全半径的个数。其余阶段仅对其做简单物理分析后绘制出六个阶段的着陆轨道。
建立模型四做相应的误差分析和敏感性分析。首先以模型二为基础进行误差分析，当主减速阶段的推力、初始质量变化时，计算嫦娥三号质量和燃料消耗速率的变化趋势。再以模型三为基础进行分析，对初始高度变化前后主减速阶段的的偏角和和着陆轨道进行对比分析并计算误差。然后进行敏感性分析，主要利用蒙特卡洛分析着陆轨道的粗避障阶段和精避障阶段月面不同地形高度，对嫦娥三号降落时所需调整概率大小的影响，接着分析嫦娥三号着陆占地面积大小对着陆调整概率的影响。
关键字：抛物线、燃料、拟牛顿法、Admas、网格化、蒙特卡洛模拟
1.问题重述
嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m。
嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段，要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。
根据上述的基本要求，请你们建立数学模型解决下面的问题：
（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。
（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
（3）对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
2.问题的分析
本文所研究的问题一主要为基础计算和物理知识,首先我们需要根据预定的着陆点的经纬度确定轨道，然后通过抛物线的运动计算出在月球着陆时的水平路程，然后计算出偏移角度，据此确定近月点的经纬度，而嫦娥三号的着陆轨道为过月球中心点的椭圆轨道，所以远月点的经纬度和近月点对称，则可以由近月点计算出远月点的经纬度。最后因为在着陆轨道上卫星的能量守恒，则可以通过势能和动能的转换来计算嫦娥三号的速度和方向。
本文所研究的问题二主要为过程的最优控制和建立嫦娥三号软着陆轨道。因为嫦娥三号的软着陆主要分为六个阶段，所以此问应分为六个阶段来求解。主减速阶段采用燃料最优制导律来分析，建立着陆坐标系，将最优燃耗软着陆问题转化为最短时间控制问题，然后得到目标函数；快速调整阶段采用重力转弯制导，对嫦娥三号进行受力分析，得到嫦娥三号的动力学模型，然后计算出燃耗最优控制，并画出仿真图；粗避障阶段采用多项式制导，首先列出加速度、速度、位移的多项式，然后通过初始状态和末端状态反解多项式系数进而求取标称轨迹；精避障阶段首先设定嫦娥三号的体型大小，然后处理数据的数量级不同，最后在整个降落区域的范围内搜索最优着陆点；由于在缓速下降和自由落体阶段中，发动机已经关闭，故仅对其做简单物理分析。最后通过整个分析得出总的着陆轨道。
本文所研究的问题三主要为着陆轨道和控制策略做误差分析和敏感度分析，
需要对问题二所设计的着陆轨道和控制策略中的发动机推力、初始速度、初始高度进行误差分析。然后进行敏感度分析，即对着陆轨道的粗避障阶段和精避障阶段月面不同地形高度对嫦娥三号降落时所需调整概率大小的影响，最后分析嫦娥三号着陆占地面积大小对着陆调整概率的影响。
3.模型的假设
假设一：嫦娥三号与月球均不受其他行星及卫星的影响
假设二：不考虑月球绕地及其他星球的公转和月球的自转
假设三：将月球近似的看做标准球体
假设四：嫦娥卫星的燃料消耗主要是在着陆的主减速阶段
假设五：软着陆的四、五、六阶段着陆轨迹基本在同一平面内
4.符号与公式的约定和说明
： G=为引力常量，m、M分别为两物体质量，R为两物体距离，为两
物体间的作用力
： 为物体质量，为物体在作用下产生的加速度
： 软着陆起始速度
： 加速度
： 平抛产生的距离
： 物体的动能（
： 物体的重力势能（
： 嫦娥三号的推力
： 偏好系数
： 降落地点总体得分
： 第段离散段的平均加速度
由于本文使用参数和公式较多，其他公式和符号在具体模型中再做说明。
5.模型的建立与求解
5.1模型一的建立
5.1.1模型的假设
由万有引力公式计算，再由牛顿第二定律计算地球和月球在近月点和远月点处的重力加速度。
表1 地球及月球在近月点和远月点的重力加速度（单位：）
三号与月球影响很小，故可忽略不计。所以本模型只考虑月球对嫦娥三号的影响。
5.1.2模型的分析
根据附件2给出的软着陆过程示意图，即嫦娥三号将在近月点15公里处以抛物线下降，相对速度从每秒1.7公里逐渐降为零。整个过程大概需要750秒，
我们将其看作匀减速运动过程。利用matlab绘制嫦娥三号绕月飞行的三维动态图，更直观的反应嫦娥三号的环月飞行，如图3（源程序见附录）：
图2 嫦娥三号绕月轨道坐标图 图3 嫦娥三号环月飞行
同时由附件二所给的嫦娥三号着陆区域和着陆点示意图可知，只要保证嫦娥三号的着陆区域在虹湾着陆区，则认为着陆成功。
为保证嫦娥三号以最大概率降落到精准的着陆点和虹湾着陆区，经分析后得出，选择以北纬44.12°作为软着陆的绕月轨道。在这种确定纬度的绕月轨道中，月球对嫦娥三号的万有引力，可以分解为两个方向。一个是绕月的向心力，一个是与绕行面相切的力，则选择最终状态为绕赤道运行更为准确。故根据实际分析，嫦娥三号的绕月平面应与南北极轴重合。
图4 嫦娥三号绕月飞行轨道分析
5.1.3模型的建立与计算
据了解，嫦娥三号主发动机是目前中国航天器上最大推力的发动机，能够产生从1500牛到7500牛的可调节推力，故可根据推力范围求取嫦娥三号的加速度范围。并用最大的加速度计算平抛产生的距离。
主减速段看作平抛运动：
起始速度
加速度的取值范围
平抛产生的距离 （
图5 嫦娥三号抛物示意图
由上图，并结合计算所得的抛物距离，得到准备着陆的点与软着陆点相差15.25°，即可算出近月点的经纬度，同时根据对称性，又可求得远月点的经纬度。
由附件所给条件可知距离月球表面15km时，速度的大小为，则此速度看作近月点速度，在稳定的轨道下，从近月点到远月点可看作重力势能和动能相互转换的过程，而远月点距离地球表面为100km，可以计算重力势能的变化，即可算出远月点的速度：
(1)
根据以上公式可得出近月点与远月点的速度（速度方向沿轨道切线方向），连同经纬度，如下表所示：
表6 近月点、远月点位置与速度
5.2模型二的建立
5.2.1模型的分析
本模型主要对主减速阶段和快速调整阶段进行初步分析
首先分析嫦娥三号在此阶段的的受力情况，假设受力与竖直方向的夹角为：
图7主减速阶段受力分析图 图8 不考虑质量变化时的受力分析
利用动量守恒定律可得：
(2)
(3)
由题目和附件可知，嫦娥三号在运行过程中有燃料的消耗，本模型分为两种情况考虑，一种为考虑质量变化，另一种为不考虑质量变化。由于主减速阶段燃料消耗很大，故作为质量变化考虑；而快速调整阶段速度很小，质量变化很小，故作为质量不变考虑。
考虑质量变化（主减速阶段），推力大小
此阶段的燃料的消耗量为
不考虑质量变化（快速调整阶段）：由于值较小，可以通过姿态调整发动机进行微调,假设此阶段质量的变化较小，则可以假设质量基本保持不变。
通过受力分析，可得到以下分析式：
最后得到燃料消耗为
(4)
5.1.2模型的建立
建立目标规划函数，计算最少的燃料消耗。由分析阶段的计算可以得出总燃料消耗量：
(5)
由表达式可以画出总燃料消耗量与质量和时间的关系
图9 总燃料消耗量与时间的关系
由图可以看出，嫦娥三号的质量随时间递增而减少，而燃料的消耗随着时间递增而增加。
5.3模型三的建立
本模型为分阶段深入分析嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
5.3.1主减速阶段制导控制律（燃料最优率制导[2]）
模型的准备
拟牛顿法是求解非线性优化问题最有效的方法之一。拟牛顿法只要求每一步迭代时知道目标函数的梯度。通过测量梯度的变化，构造一个目标函数的模型使之足产生超线性收敛性。构造目标函数在当前迭代的二次模型和割线公式
预估—校正算法的方法包括三步四阶Adams外插法和三步四阶Adams内插法为了保证计算得精度，本文采用内插法
模型的分析与建立 嫦娥三号主减速阶段从距离月球表面15km开始，由初
速度为 开始主减速。建立二维模型描述嫦娥三号在此阶段的运动。令月心O为坐标原点, y 指向动力下降段的开始制动点, x 向着陆器的开始运动方向，见下图：
图10 着陆坐标系
由坐标系可建立嫦娥三号的质心动力学方程，描述如方程组（6）：
（6）
式中: ,,和分别为嫦娥三号的月心距、极角、角速度和质量;
为嫦娥三号沿方向上的速度;
为制动发动机的推力(固定的常值或0) ;
为其比冲;
为月球引力常数;
为发动机推力与当地水平线的夹角即推力方向角。
动力下降的初始条件由霍曼变轨后的椭圆轨道的近月点确定,终端条件
为嫦娥三号在月面实现软着陆。令初始时刻,终端时刻不定,则此过程的约束条件可以表示为方程组（7）：
（7）
对的求解 月球软着陆的最优轨道设计就是要在满足上述初始条件和终端
约束条件的前提下, 调整推力大小和方向,使得嫦娥三号实现燃料最优软着陆,则设燃料最优目标函数为表达式（8）：
（8）
在无奇异情况下,推力应为开关控制。要么以最大推力工作,要么以最小推力工作。但为了简化问题,采用常值推力假设，即认为制动发动机一直以最大推力工作。这一方法一方面有利于优化，另一方面可降低发动机复杂性。采用常值推力假设后,月球最优燃耗软着陆问题转化为最短时间控制问题,即寻找实现软着陆的最短时间，求解步骤如下:
：确定一终端时间,满足条件
：求解无约束最优控制问题状态方程式,终端时间为，性能指标为：
（9）
其中下标表示在时刻的取值。
： 根据终端能量特性修正，然后返回，直到。
终端时刻的初始值估计，由于软着陆时着陆器能量为零,可知推力作用主要是抵消能量，将该能量等效为动能,则可推出等效速度为
假设采用脉冲推力模式,将该速度抵消需要消耗的燃料量为
而对于实际的有限推力模式，与相对应的时间为
（10）
式中为发动机燃料秒流量
最终得计算结果为：
因脉冲推力比有限推力消耗的燃料量少，所以使得该计算结果偏小。 目标函数的求解 第二阶段垂直方向上的减速最大值为
由文献可知,为使卫星在第六阶段自由落体，则快速调整阶段的速度范围为：
假设主减速阶段卫星以一定角度提供向上的推动力，则等效速度为
由于值较小，故可以忽略不计。
此问题为终端时间固定型无约束最优控制问题,本模型将其转化为非线性规划问题,然后借助于拟牛顿法和四阶Admas 预测-校正积分格式快速求解。为保证优化精度,转化方法采用计算量稍大但精度较高的直接离散化方法。
直接离散化方法将整个最优控制过程分成若干个时间段,时间段之间的端点称为节点;选择节点处的控制变量作为未知参数,通过插值得到整个最优控制过程的控制变量积分状态方程;根据这些控制变量积分状态方程形成目标函数,得到一个无约束数学规划问题。具体如下：
(1) 将整个飞行时间分为N 个时间段,形成N+ 1 个时间节点 ( i = 0 ,1 , ⋯,
N) ,取时刻的控制量为优化变量,共有N + 1 个变量；
(2) 整个飞行过程的控制量可以通过在各时间节点处线性插值得到；
(3) 采用拟牛顿法和四阶Admas预测-校正积分，得到从到积分状态方程(6)
和目标函数(9)。
图11 偏角和垂直速度随时间变化的趋势
5.3.2快速调整段制导律（重力转弯制导[4]）
模型的分析 由于在最终着陆段中，嫦娥三号的距月面距离只有 2 千米左
右，远远小于月球的半径 1738 千米，因此在建模时可以忽略月球的曲率，将月面近似看为水平面；且考虑到在最终着陆段中嫦娥三号的切向速度只有几十米每秒，设切向速度给嫦娥三号所带来的离心加速度为，月球半径为。因为嫦娥三号的切向速度为，则计算切向速度给嫦娥三号所带来的离心加速度公式为：
因此可以忽略嫦娥三号的离心加速度，只考虑重力加速度。
模型的建立 假设嫦娥三号的下降轨迹在一个平面内，设制动发动机的比冲
为，秒耗量为，嫦娥三号的垂直高度为，切向速度为，质量为，制动发动机的推力方向与垂直方向夹角为。在以上假设条件下，我们对嫦娥三号进行受力分析，可以得到嫦娥三号的动力学模型为：
（12）
模型的最优解 为了使嫦娥三号在最终着陆段中的燃料消耗达到最小，则设
嫦娥三号软着陆燃料消耗为：
（13）
对于重力转弯制导法下的软着陆模型，推力的燃耗最优控制是开关控制，而且开关次数最多不会超过 1 次。要实现嫦娥三号的终端状态约束，嫦娥三号只能先进行自由落体，直到开关切换函数为 0 时，制动发动机工作，嫦娥三号进行制动减速，直至在到达月面时减速为 0,仿真图如下所示：
图12 快速调整阶段运动状态
5.3.3粗避障段制导律（多项式制导[5] ）
模型的分析 嫦娥三号软着陆粗避障阶段持续时间较短，所以需要设计有效
的制导律使探测器能在有限的时间内跟踪上标称轨迹，外部环境的干扰是影响着陆精度的主要因素。所以，本模型首先给出了多项式，然后通过初始状态和末端状态反解多项式系数进而求取标称轨迹，然后设计终端滑模制导律跟踪标称轨迹。
模型的建立 多项式形式的标称轨迹规划一般假设系统状态变量为多项式，
基于边界条件和着陆时间解相关系数。对于嫦娥三号粗避障阶段，首先可以将着陆器的加速度表示为二次多项式的形式：
(14)
其中，和分别为待定常数矢量。对式（14）等式两边积分可以得到嫦娥三号的速度矢量和位置矢量的表达式为：
+ (15)
+ (16)
给定着陆时间和初末端状态的情况下，可以解出：
=
模型的计算和分析 生成标称轨道的仿真参数为着陆器在着陆点平移坐标
下的初始位置矢量 ,初始速度矢量，着陆时间为，将参数代入到式（17）可得常矢量为：
基于光学图像的粗障碍检测就是利用月球岩石和坑的图像特征识别大障碍, 确定安全区域。根据岩石和坑的特征，本文选取避障原则如下式：
图13 粗避障阶段的等高线
将此区域图片看做的矩阵，进一步分割为个的矩阵。根据组成地面高度的矩阵，利用var函数求解计算每一个矩阵的方差。方差的大小代表地面的平坦程度。
图14 粗避障阶段最优着陆点
图中白色区域为方差最小点，即为不考虑避障阶段速度增量的值时，需要搜寻的最优着陆区域。
5.3.4精避障阶段
精避障阶段，推力和姿态发动机的比冲较小且时间短，不将比冲燃料消耗计算在内。为了在整个降落区域的范围内搜索最优着陆点，将图片区域网格化处理为的矩阵，选择最优区域的准则为总高度和总平坦度值的大小。
用Min-Max标准化方法消除数量级的不同
设置偏好系数表示区域总高度对降落点得分的影响，表示区域总平坦度对降落点得分的影响。则降落地点总体得分。
图15 精避障阶段的等高线
对着陆所占用的不同区域下的计算，得出结论在占用区域面积时，最优点为的附近区域。
表28 精避障阶段最优降落点
根据需要着陆的大小，对整个各个区域进行搜索满足的点，即为可选择的降落点
5.3.5轨道的确定
上文对着陆轨道的六个阶段进行分析，主减速阶段嫦娥卫星的速度和质量变化最大，对轨道的计算也最为重要。对于缓速下降和自由落体阶段，由于发动机已经关闭，则对于最优控制和轨道设计不必过多分析。通过前面四个阶段的分析和自由落体的规律，得出最终的着陆轨道如下图：
图16 最终着陆轨道的设计
5.4误差分析与敏感度分析
主要对模型三设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
5.4.1误差分析
本模型主要分析发动机推力误差、初始速度误差、初始高度误差等。 发动机推力误差：主要分析为主减速阶段推力变化和嫦娥三号初始质量变化对嫦娥三号质量和燃料消耗的影响。
首先设定嫦娥三号的推力为最大推力7500N,然后将分别乘以1.1、0.9，观察的变化对嫦娥三号质量和燃料消耗的影响，如下图：
图17 推力改变时的误差分析
由图可以看出，嫦娥三号的推力变化会引起嫦娥卫星的质量和燃料消耗的变化，推力越大，质量改变越小，燃料消耗越少。
由题目所给条件可知嫦娥三号的初始质量为 =2400kg，然后将嫦娥三号初始质量乘以1.1和0.9，观察此时嫦娥三号的质量和燃料消耗的变化。
图18嫦娥三号质量改变时的误差分析
由图可知，嫦娥三号的初始质量的变化会引起嫦娥三号的质量和燃料消耗的变化，初始质量越大嫦娥三号的质量变化越大，燃料消耗的越多。
对主减速阶段的初始速度和初始高度进行误差分析，
嫦娥三号的预定着陆点
海拔为-2641m，则将主减速阶段的高度设置为15Km至17.641Km之间。将其与原有状态下的运动状态相互比较。仅考虑切向速度变化，根据燃料最优制导模型的计算方法，利用四阶龙格-库塔公式和拟牛顿法将主减速的30个阶段嫦娥三号偏角的变化与原变化进行比较，如下图：
图19 偏角的变化
上图蓝线表示原的变化，绿线为改变切向速度时的变化，红线为两者的误差，可以看出前期原偏角大于改变后的偏角，后期则相反。误差也随着时间变少。由误差计算公式 ，计算偏角总误差为-9.49% 。
根据已求得的偏角的的值，将主减速段运动路径分割为30个阶段，并将轨道离散化
图20 初始高度变化时轨道的变化
图21 初始高度变化的轨道离散化
图18的红线为原高度时轨道变化，粉红线为改变原高度时的轨道变化。由误差公式可得，在主减速阶段的误差为,误差率为 。
已知嫦娥三号的初始比冲量为2940，将其分别乘以0.9、1.1，即改变比冲量，观察嫦娥三号质量和燃料消耗的变化。
图22 比冲量变化时轨道特性的变化
由图可以看出，比冲量的值越大，嫦娥三号的质量变化越大，燃料消耗越大。
5.4.2 敏感度分析
粗避障段：粗避障段的范围是距离月面2.4km到100m区间，其主要是要求避开大的陨石坑，实现在设计着陆点上方100m处悬停，并初步确定落月地点。 将附件所给图片网格化为2300×2300的矩阵，本文根据处的月球高度，得到避障原则： （19）
使用matlab软件并采用用蒙特卡洛的方法进行1000次仿真（源程序见附录），模拟分析月面不同地形高度对嫦娥三号降落时所需调整概率大小的影响。
图23 粗避障阶段地形 图24 粗避障阶段不同降落高度所需调整 精避障段：精细避障段的区间是距离月面100m到30m。要求嫦娥三号悬停在距离月面100m处，对着陆点附近区域100m范围内拍摄图像，并获得三维数字高程图。分析三维数字高程图，避开较大的陨石坑，确定最佳着陆地点，实现在着陆点上方30m处水平方向速度为0m/s。
与粗避障一样，在满足同样的避障原则下，分析月面不同降落地形高度对嫦娥三号降落时所需调整概率大小的影响。
图25 精避障阶段地形 图26 精避障阶段不同降落高度所需调整
在精确避障阶段的避障原则下，为了研究嫦娥三号在降落时占地面积大小对轨道降落的敏感度的影响，选择和这两个数据作为嫦娥三号降落时的占地面积。并用这两个数据在matlab软件中做50次模拟比较（源程序见附录）。
图27 两种不同着陆占地面积着陆比较
由图可以看出，的占地面积的非调整降落次数高于。，由此推测着陆占地面积越大，可直接顺利着陆的概率越小。
分别选用六个不同的占地面积，对其进行1000次模拟，计算出1000次模拟中无需调整即可顺利着陆的次数，如下表：
表29 不同占地面积无需调整即可顺利着陆的次数
的概率越大。
6.模型的评价与改进方向
6.1模型的评价
6.1.1模型一的评价
模型一分别以着陆点的经度和纬度作为准备着陆轨道，选取经度不变的轨道处于稳定状态，不需要产生推力，此种轨道保证了燃料消耗的最优。选取纬度不变的轨道则保证了当平抛距离较大且难以精确确定时以最大概率降落在着陆区域内。
6.1.2 模型三的评价
在模型三中建立的主减速阶段燃料最优目标规划函数，利用时间逼近法快速求解月球最优软着陆问题。对于终端时间固定型最优控制问题,将其直接离散化为非线性规划问题,采用拟牛顿法和四阶Admas预测-校正积分方法快速求解。此方法优化精度较高,收敛速度快，比近年较为流行的智能算法(如遗传算法等)，减少了计算量且更符合实际需求和精度要求。
对嫦娥三号软着陆的其他阶段也分别建立了动力学模型。并且分段建立了最优目标函数。确定了着陆轨道。在精避障阶段，综合考虑了着陆位置的总高度和总平坦度，对不同数量级的数据标准化，设置偏好系数后对所有点进行全局搜索，得出了最优降落策略和最优降落点。
6.1.3 模型四的评价
模型四对误差进行了多方面的分析，包括对的最大推力、初始速度的变化的轨道特性分析、进一步对主减速阶段的偏角的趋势分析、计算误差和灵敏度。
6.1.4 模型的不足与改进方向
模型的不足：由于轨道的复杂多变性，本文简化了模型的数学推导，将高度变化引起的轨道路径长度变化忽略，只重点考虑和计算了主减速阶段的轨道特性，造成了设计的轨道系统的误差。且未对轨道路径做出明显的全局优化。
改进方向：
1. 将着陆轨道的六个阶段燃料量作为规划函数，将自适应遗传算法与模拟退火
算法相结合,形成一种自适应模拟退火遗传算法，增强轨道路径设计的整体搜索能力。
2. 根据月球岩石和坑的特征, 设计了粗障碍识别和安全着陆区选取算法: 1)
图像直方图分析; 2)K 均值聚类; 3) 过亮障碍识别; 4) 过暗障碍识别; 5)
纹理障碍识别; 6) 采用螺旋搜索算法确定每个单元格的安全半径(图3); 7) 根据安全半径, 选取候选安全着陆点; 8) 评估候选安全着陆点避障所需的速度增量; 9) 根据安全半径和速度增量评价值, 综合确定安全着陆点.
7.参考文献
[1] 张德丰，MATLAB数值分析，北京：机械工业出版社，2012。
[2] 赵吉松，谷良贤，高原，月球软着陆轨道的时间逼近法快速优化设计[J]，宇航学报，第29卷第5期：1-5，2008. 9。
[3] 朱建丰，徐世杰，基于自适应模拟退火遗传算法的月球软着陆轨道优化[J]，航空学报，第28卷第4期：2-3，2007.7。
[4] 于彦波，火星探测器动力下降段制导律研究[D]，哈尔滨，哈尔滨工业大学， 2013。
[5] 张仲满，月球软着陆的制导算法研究[D]，哈尔滨，哈尔滨工业大学， 2009。
[6]田青，常微分方程初值问题数值解的实现与分析，http://www.doc88.com/p-[1**********]3.html,2014.09.13。
[7]张洪华，梁俊等，嫦娥三号自主避障软着陆控制技术[J]，中国科学：科学计算，第44卷第6期：2-4,2014。
附录
蒙特卡洛分析不同降落地形高度选取调整概率
I2=imread ('data.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
c=zeros(200,1);
d=zeros(200,1);
for j=1:1:200
for i=1:1000
aa=floor(2299*rand(1)+1);
bb=floor(2299*rand(1)+1);
if pp2(aa,bb)
c(j,1)=c(j,1)+1;
else
d(j,1)=d(j,1)+1;
end
end
end
plot(1:200,c/1000,1:200,d/1000);
title('不同降落地形高度选取调整概率');
xlabel('地形高度');
ylabel('调整与非调整概率');
legend('不需调整','需要调整');
精避障阶段50次模拟分析嫦娥安好占地面积大小对着陆的影响： d1=zeros(50,1);
for j=1:50
for i=1:1000
aa=floor(979*rand(1)+11);
bb=floor(979*rand(1)+11);
aa1=aa-7;
aa2=aa+7;
bb1=bb-7;
bb2=bb+7;
c=0;
for cc=aa1:aa2
for dd=bb1:bb2
if pp2(cc,dd)>=50&pp2(cc,dd)
c=c+1;
end
if c==15*15;
d(j,1)=d(j,1)+1;
end
end
end
end
end
plot(1:50,d)
快速调整阶段运动状态
v1=-90;
a=(7500-1200*1.6)/1200;
d=0;
for t=0.1:0.1:14
d=d+1;
x(1,d)=3000+v1*t+0.5*a*t^2;
end
plot(0.1:0.1:14,x')
title('快速调整段制导律下运动状态');
xlabel('时间/t');
ylabel('高度/m');
总燃料消耗量与时间的关系
F=7500;
m=2400;
ve=2940;
v=22;
g=1.63;
a=3.125;
m1=zeros(1,487);
for i=1:487
m1=F/ve;
m=m-m1-m*g/ve;
A(1,i)=m;
F=m*a;
end
subplot(1,2,1)
plot(1:487,A)
title('速度');
xlabel('时间/t');
ylabel('嫦娥卫星质量/Kg');
subplot(1,2,2)
plot(1:486,A(2:487)-A(1:486))
title('能料消耗的变化');
xlabel('时间/t');
ylabel('质量/Kg');
快速调整阶段最优降落点
I2=imread ('data1.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
a=input('请输入嫦娥号的大小a=');
for aa=a+1:1:1000-a
for bb=a+1:1:1000-a
AA=pp2(aa-a+1:aa+a,bb-a+1:bb+a);
A(aa-a,bb-a)=sum(sum(AA));
if A(aa-a,bb-a)==2389
aa
bb
end
end
end
初始高度变化时轨道的变化
H1=15000;
H2=15000+2641;
X=0;
X2=0;
i=0:1:510;
I=0:1:600;
X=1700.*i-0.5*3.125*(i.^2);
X2=1700.*I-0.5*2.878*(I.^2);
Y=H1-22.*i;
Y7=min(H2-22*1.1*I);
Y5=H2-22.*I*1.1;
i1=510:559;
I1=600:640;
X1=57*(i1-510)+460590-0.5*2.*((i1-510).^2);
Y1=min(Y)-22.*(i1-510);
Y8=Y7-22*1.1.*(I1-600);
X3=460982:0.1:460990;
Y3=0:2700/80:2700;
Y_1=0:100:3500;
X_1=linspace(501967,500720,36)
plot(X,Y,'r',X2,Y5,'m',X1,Y1,'k',X3,Y3,'g',X_1,Y_1); title('初始高度变化时轨道的变化');
xlabel('和近月点的距离');
ylabel('距地面高度/m');
text(250000,11000,'主减速阶段');
text(420000,3500,'快速调整阶段');
text(380000,2000,'粗避障阶段');
text(380000,1000,'其他阶段');
网格化矩阵方差
I2=imread ('data.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
A=zeros(100,100);
for m=1:22
for n=1:22
for i=(1:100)+m*100
for j=(1:100)+m*100
A((i-100)*m,(j-100)*n)=pp2(i,j);
c=var(A(:))
end
end
end
end
搜索可选择的着陆点
I2=imread ('data1.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
a=zeros(1000,1000);
for i=1:1000
for j=1:1000
if pp2(i,j)>=60&pp2(i,j)
a(i,j)=1;
else
a(i,j)=0;
end
end
end
bb=input('输入bb=');
for m1=bb+1:999-bb
for m2=bb+1:999-bb
if sum(sum(a(m1-bb:m1+bb,m2-bb:m2+bb)))==(2*bb+1)^2 m1
m2
else
continue;
end
end
end
2014高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
嫦娥三号软着路轨道设计与控制策略
摘要
本文主要为分阶段研究嫦娥三号的软着陆轨道设计与最优控制策略。 建立模型一确定近月点和远月点的位置，以及嫦娥三号速度大小与方向。首先以月球中心为坐标原点建立空间坐标系，根据计算的作用力可知地球影响较小，故忽略不计。然后将嫦娥三号软着陆看作抛物线的运动过程，计算在最大推力下的减速运动，求得月面偏移距离为，由此计算出偏移角度为15.25°。从而得出近月点和远月点的经纬度分别为（34.76°W，44.12°N）和（34.76°E，44.12°S）。最后在软着陆的椭圆轨道上，由动力势能和重力势能的变化，计算出嫦娥三号在远月点和近月点的速度分别为，沿轨道切线方向。
建立模型二和模型三确定着陆轨道和在6个阶段的最优控制策略。模型二主要对主减速阶段和快速调整阶段进行初步分析。模型三分六个阶段确定轨道和最优控制策略，主减速阶段建立目标函数燃料，假设推力最大,将最优燃耗软着陆问题转化为最短时间控制问题，然后采用拟牛顿法和四阶Admas预测-校正得到；快速调整阶段采用重力转弯制导，在假设条件下对嫦娥三号进行受力分析，得到嫦娥三号的动力学模型，然后通过开关控制得到燃耗最优控制，并画出仿真图；粗避障阶段采用多项式制导，通过初始状态和末端状态反解多项式系数进而求取标称轨迹，然后将避障区域网格化，比较网格内的方差大小确定最优区域范围；精避障阶段需在满足本文提出的避障原则式下搜索全局最优解，以网格区域总体得分作为目标函数，得到最优区域为坐标附近，并以螺旋搜索法搜索安全半径的个数。其余阶段仅对其做简单物理分析后绘制出六个阶段的着陆轨道。
建立模型四做相应的误差分析和敏感性分析。首先以模型二为基础进行误差分析，当主减速阶段的推力、初始质量变化时，计算嫦娥三号质量和燃料消耗速率的变化趋势。再以模型三为基础进行分析，对初始高度变化前后主减速阶段的的偏角和和着陆轨道进行对比分析并计算误差。然后进行敏感性分析，主要利用蒙特卡洛分析着陆轨道的粗避障阶段和精避障阶段月面不同地形高度，对嫦娥三号降落时所需调整概率大小的影响，接着分析嫦娥三号着陆占地面积大小对着陆调整概率的影响。
关键字：抛物线、燃料、拟牛顿法、Admas、网格化、蒙特卡洛模拟
1.问题重述
嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m。
嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段，要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。
根据上述的基本要求，请你们建立数学模型解决下面的问题：
（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。
（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
（3）对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
2.问题的分析
本文所研究的问题一主要为基础计算和物理知识,首先我们需要根据预定的着陆点的经纬度确定轨道，然后通过抛物线的运动计算出在月球着陆时的水平路程，然后计算出偏移角度，据此确定近月点的经纬度，而嫦娥三号的着陆轨道为过月球中心点的椭圆轨道，所以远月点的经纬度和近月点对称，则可以由近月点计算出远月点的经纬度。最后因为在着陆轨道上卫星的能量守恒，则可以通过势能和动能的转换来计算嫦娥三号的速度和方向。
本文所研究的问题二主要为过程的最优控制和建立嫦娥三号软着陆轨道。因为嫦娥三号的软着陆主要分为六个阶段，所以此问应分为六个阶段来求解。主减速阶段采用燃料最优制导律来分析，建立着陆坐标系，将最优燃耗软着陆问题转化为最短时间控制问题，然后得到目标函数；快速调整阶段采用重力转弯制导，对嫦娥三号进行受力分析，得到嫦娥三号的动力学模型，然后计算出燃耗最优控制，并画出仿真图；粗避障阶段采用多项式制导，首先列出加速度、速度、位移的多项式，然后通过初始状态和末端状态反解多项式系数进而求取标称轨迹；精避障阶段首先设定嫦娥三号的体型大小，然后处理数据的数量级不同，最后在整个降落区域的范围内搜索最优着陆点；由于在缓速下降和自由落体阶段中，发动机已经关闭，故仅对其做简单物理分析。最后通过整个分析得出总的着陆轨道。
本文所研究的问题三主要为着陆轨道和控制策略做误差分析和敏感度分析，
需要对问题二所设计的着陆轨道和控制策略中的发动机推力、初始速度、初始高度进行误差分析。然后进行敏感度分析，即对着陆轨道的粗避障阶段和精避障阶段月面不同地形高度对嫦娥三号降落时所需调整概率大小的影响，最后分析嫦娥三号着陆占地面积大小对着陆调整概率的影响。
3.模型的假设
假设一：嫦娥三号与月球均不受其他行星及卫星的影响
假设二：不考虑月球绕地及其他星球的公转和月球的自转
假设三：将月球近似的看做标准球体
假设四：嫦娥卫星的燃料消耗主要是在着陆的主减速阶段
假设五：软着陆的四、五、六阶段着陆轨迹基本在同一平面内
4.符号与公式的约定和说明
： G=为引力常量，m、M分别为两物体质量，R为两物体距离，为两
物体间的作用力
： 为物体质量，为物体在作用下产生的加速度
： 软着陆起始速度
： 加速度
： 平抛产生的距离
： 物体的动能（
： 物体的重力势能（
： 嫦娥三号的推力
： 偏好系数
： 降落地点总体得分
： 第段离散段的平均加速度
由于本文使用参数和公式较多，其他公式和符号在具体模型中再做说明。
5.模型的建立与求解
5.1模型一的建立
5.1.1模型的假设
由万有引力公式计算，再由牛顿第二定律计算地球和月球在近月点和远月点处的重力加速度。
表1 地球及月球在近月点和远月点的重力加速度（单位：）
三号与月球影响很小，故可忽略不计。所以本模型只考虑月球对嫦娥三号的影响。
5.1.2模型的分析
根据附件2给出的软着陆过程示意图，即嫦娥三号将在近月点15公里处以抛物线下降，相对速度从每秒1.7公里逐渐降为零。整个过程大概需要750秒，
我们将其看作匀减速运动过程。利用matlab绘制嫦娥三号绕月飞行的三维动态图，更直观的反应嫦娥三号的环月飞行，如图3（源程序见附录）：
图2 嫦娥三号绕月轨道坐标图 图3 嫦娥三号环月飞行
同时由附件二所给的嫦娥三号着陆区域和着陆点示意图可知，只要保证嫦娥三号的着陆区域在虹湾着陆区，则认为着陆成功。
为保证嫦娥三号以最大概率降落到精准的着陆点和虹湾着陆区，经分析后得出，选择以北纬44.12°作为软着陆的绕月轨道。在这种确定纬度的绕月轨道中，月球对嫦娥三号的万有引力，可以分解为两个方向。一个是绕月的向心力，一个是与绕行面相切的力，则选择最终状态为绕赤道运行更为准确。故根据实际分析，嫦娥三号的绕月平面应与南北极轴重合。
图4 嫦娥三号绕月飞行轨道分析
5.1.3模型的建立与计算
据了解，嫦娥三号主发动机是目前中国航天器上最大推力的发动机，能够产生从1500牛到7500牛的可调节推力，故可根据推力范围求取嫦娥三号的加速度范围。并用最大的加速度计算平抛产生的距离。
主减速段看作平抛运动：
起始速度
加速度的取值范围
平抛产生的距离 （
图5 嫦娥三号抛物示意图
由上图，并结合计算所得的抛物距离，得到准备着陆的点与软着陆点相差15.25°，即可算出近月点的经纬度，同时根据对称性，又可求得远月点的经纬度。
由附件所给条件可知距离月球表面15km时，速度的大小为，则此速度看作近月点速度，在稳定的轨道下，从近月点到远月点可看作重力势能和动能相互转换的过程，而远月点距离地球表面为100km，可以计算重力势能的变化，即可算出远月点的速度：
(1)
根据以上公式可得出近月点与远月点的速度（速度方向沿轨道切线方向），连同经纬度，如下表所示：
表6 近月点、远月点位置与速度
5.2模型二的建立
5.2.1模型的分析
本模型主要对主减速阶段和快速调整阶段进行初步分析
首先分析嫦娥三号在此阶段的的受力情况，假设受力与竖直方向的夹角为：
图7主减速阶段受力分析图 图8 不考虑质量变化时的受力分析
利用动量守恒定律可得：
(2)
(3)
由题目和附件可知，嫦娥三号在运行过程中有燃料的消耗，本模型分为两种情况考虑，一种为考虑质量变化，另一种为不考虑质量变化。由于主减速阶段燃料消耗很大，故作为质量变化考虑；而快速调整阶段速度很小，质量变化很小，故作为质量不变考虑。
考虑质量变化（主减速阶段），推力大小
此阶段的燃料的消耗量为
不考虑质量变化（快速调整阶段）：由于值较小，可以通过姿态调整发动机进行微调,假设此阶段质量的变化较小，则可以假设质量基本保持不变。
通过受力分析，可得到以下分析式：
最后得到燃料消耗为
(4)
5.1.2模型的建立
建立目标规划函数，计算最少的燃料消耗。由分析阶段的计算可以得出总燃料消耗量：
(5)
由表达式可以画出总燃料消耗量与质量和时间的关系
图9 总燃料消耗量与时间的关系
由图可以看出，嫦娥三号的质量随时间递增而减少，而燃料的消耗随着时间递增而增加。
5.3模型三的建立
本模型为分阶段深入分析嫦娥三号的着陆轨道和在6个阶段的最优控制策略。
5.3.1主减速阶段制导控制律（燃料最优率制导[2]）
模型的准备
拟牛顿法是求解非线性优化问题最有效的方法之一。拟牛顿法只要求每一步迭代时知道目标函数的梯度。通过测量梯度的变化，构造一个目标函数的模型使之足产生超线性收敛性。构造目标函数在当前迭代的二次模型和割线公式
预估—校正算法的方法包括三步四阶Adams外插法和三步四阶Adams内插法为了保证计算得精度，本文采用内插法
模型的分析与建立 嫦娥三号主减速阶段从距离月球表面15km开始，由初
速度为 开始主减速。建立二维模型描述嫦娥三号在此阶段的运动。令月心O为坐标原点, y 指向动力下降段的开始制动点, x 向着陆器的开始运动方向，见下图：
图10 着陆坐标系
由坐标系可建立嫦娥三号的质心动力学方程，描述如方程组（6）：
（6）
式中: ,,和分别为嫦娥三号的月心距、极角、角速度和质量;
为嫦娥三号沿方向上的速度;
为制动发动机的推力(固定的常值或0) ;
为其比冲;
为月球引力常数;
为发动机推力与当地水平线的夹角即推力方向角。
动力下降的初始条件由霍曼变轨后的椭圆轨道的近月点确定,终端条件
为嫦娥三号在月面实现软着陆。令初始时刻,终端时刻不定,则此过程的约束条件可以表示为方程组（7）：
（7）
对的求解 月球软着陆的最优轨道设计就是要在满足上述初始条件和终端
约束条件的前提下, 调整推力大小和方向,使得嫦娥三号实现燃料最优软着陆,则设燃料最优目标函数为表达式（8）：
（8）
在无奇异情况下,推力应为开关控制。要么以最大推力工作,要么以最小推力工作。但为了简化问题,采用常值推力假设，即认为制动发动机一直以最大推力工作。这一方法一方面有利于优化，另一方面可降低发动机复杂性。采用常值推力假设后,月球最优燃耗软着陆问题转化为最短时间控制问题,即寻找实现软着陆的最短时间，求解步骤如下:
：确定一终端时间,满足条件
：求解无约束最优控制问题状态方程式,终端时间为，性能指标为：
（9）
其中下标表示在时刻的取值。
： 根据终端能量特性修正，然后返回，直到。
终端时刻的初始值估计，由于软着陆时着陆器能量为零,可知推力作用主要是抵消能量，将该能量等效为动能,则可推出等效速度为
假设采用脉冲推力模式,将该速度抵消需要消耗的燃料量为
而对于实际的有限推力模式，与相对应的时间为
（10）
式中为发动机燃料秒流量
最终得计算结果为：
因脉冲推力比有限推力消耗的燃料量少，所以使得该计算结果偏小。 目标函数的求解 第二阶段垂直方向上的减速最大值为
由文献可知,为使卫星在第六阶段自由落体，则快速调整阶段的速度范围为：
假设主减速阶段卫星以一定角度提供向上的推动力，则等效速度为
由于值较小，故可以忽略不计。
此问题为终端时间固定型无约束最优控制问题,本模型将其转化为非线性规划问题,然后借助于拟牛顿法和四阶Admas 预测-校正积分格式快速求解。为保证优化精度,转化方法采用计算量稍大但精度较高的直接离散化方法。
直接离散化方法将整个最优控制过程分成若干个时间段,时间段之间的端点称为节点;选择节点处的控制变量作为未知参数,通过插值得到整个最优控制过程的控制变量积分状态方程;根据这些控制变量积分状态方程形成目标函数,得到一个无约束数学规划问题。具体如下：
(1) 将整个飞行时间分为N 个时间段,形成N+ 1 个时间节点 ( i = 0 ,1 , ⋯,
N) ,取时刻的控制量为优化变量,共有N + 1 个变量；
(2) 整个飞行过程的控制量可以通过在各时间节点处线性插值得到；
(3) 采用拟牛顿法和四阶Admas预测-校正积分，得到从到积分状态方程(6)
和目标函数(9)。
图11 偏角和垂直速度随时间变化的趋势
5.3.2快速调整段制导律（重力转弯制导[4]）
模型的分析 由于在最终着陆段中，嫦娥三号的距月面距离只有 2 千米左
右，远远小于月球的半径 1738 千米，因此在建模时可以忽略月球的曲率，将月面近似看为水平面；且考虑到在最终着陆段中嫦娥三号的切向速度只有几十米每秒，设切向速度给嫦娥三号所带来的离心加速度为，月球半径为。因为嫦娥三号的切向速度为，则计算切向速度给嫦娥三号所带来的离心加速度公式为：
因此可以忽略嫦娥三号的离心加速度，只考虑重力加速度。
模型的建立 假设嫦娥三号的下降轨迹在一个平面内，设制动发动机的比冲
为，秒耗量为，嫦娥三号的垂直高度为，切向速度为，质量为，制动发动机的推力方向与垂直方向夹角为。在以上假设条件下，我们对嫦娥三号进行受力分析，可以得到嫦娥三号的动力学模型为：
（12）
模型的最优解 为了使嫦娥三号在最终着陆段中的燃料消耗达到最小，则设
嫦娥三号软着陆燃料消耗为：
（13）
对于重力转弯制导法下的软着陆模型，推力的燃耗最优控制是开关控制，而且开关次数最多不会超过 1 次。要实现嫦娥三号的终端状态约束，嫦娥三号只能先进行自由落体，直到开关切换函数为 0 时，制动发动机工作，嫦娥三号进行制动减速，直至在到达月面时减速为 0,仿真图如下所示：
图12 快速调整阶段运动状态
5.3.3粗避障段制导律（多项式制导[5] ）
模型的分析 嫦娥三号软着陆粗避障阶段持续时间较短，所以需要设计有效
的制导律使探测器能在有限的时间内跟踪上标称轨迹，外部环境的干扰是影响着陆精度的主要因素。所以，本模型首先给出了多项式，然后通过初始状态和末端状态反解多项式系数进而求取标称轨迹，然后设计终端滑模制导律跟踪标称轨迹。
模型的建立 多项式形式的标称轨迹规划一般假设系统状态变量为多项式，
基于边界条件和着陆时间解相关系数。对于嫦娥三号粗避障阶段，首先可以将着陆器的加速度表示为二次多项式的形式：
(14)
其中，和分别为待定常数矢量。对式（14）等式两边积分可以得到嫦娥三号的速度矢量和位置矢量的表达式为：
+ (15)
+ (16)
给定着陆时间和初末端状态的情况下，可以解出：
=
模型的计算和分析 生成标称轨道的仿真参数为着陆器在着陆点平移坐标
下的初始位置矢量 ,初始速度矢量，着陆时间为，将参数代入到式（17）可得常矢量为：
基于光学图像的粗障碍检测就是利用月球岩石和坑的图像特征识别大障碍, 确定安全区域。根据岩石和坑的特征，本文选取避障原则如下式：
图13 粗避障阶段的等高线
将此区域图片看做的矩阵，进一步分割为个的矩阵。根据组成地面高度的矩阵，利用var函数求解计算每一个矩阵的方差。方差的大小代表地面的平坦程度。
图14 粗避障阶段最优着陆点
图中白色区域为方差最小点，即为不考虑避障阶段速度增量的值时，需要搜寻的最优着陆区域。
5.3.4精避障阶段
精避障阶段，推力和姿态发动机的比冲较小且时间短，不将比冲燃料消耗计算在内。为了在整个降落区域的范围内搜索最优着陆点，将图片区域网格化处理为的矩阵，选择最优区域的准则为总高度和总平坦度值的大小。
用Min-Max标准化方法消除数量级的不同
设置偏好系数表示区域总高度对降落点得分的影响，表示区域总平坦度对降落点得分的影响。则降落地点总体得分。
图15 精避障阶段的等高线
对着陆所占用的不同区域下的计算，得出结论在占用区域面积时，最优点为的附近区域。
表28 精避障阶段最优降落点
根据需要着陆的大小，对整个各个区域进行搜索满足的点，即为可选择的降落点
5.3.5轨道的确定
上文对着陆轨道的六个阶段进行分析，主减速阶段嫦娥卫星的速度和质量变化最大，对轨道的计算也最为重要。对于缓速下降和自由落体阶段，由于发动机已经关闭，则对于最优控制和轨道设计不必过多分析。通过前面四个阶段的分析和自由落体的规律，得出最终的着陆轨道如下图：
图16 最终着陆轨道的设计
5.4误差分析与敏感度分析
主要对模型三设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。
5.4.1误差分析
本模型主要分析发动机推力误差、初始速度误差、初始高度误差等。 发动机推力误差：主要分析为主减速阶段推力变化和嫦娥三号初始质量变化对嫦娥三号质量和燃料消耗的影响。
首先设定嫦娥三号的推力为最大推力7500N,然后将分别乘以1.1、0.9，观察的变化对嫦娥三号质量和燃料消耗的影响，如下图：
图17 推力改变时的误差分析
由图可以看出，嫦娥三号的推力变化会引起嫦娥卫星的质量和燃料消耗的变化，推力越大，质量改变越小，燃料消耗越少。
由题目所给条件可知嫦娥三号的初始质量为 =2400kg，然后将嫦娥三号初始质量乘以1.1和0.9，观察此时嫦娥三号的质量和燃料消耗的变化。
图18嫦娥三号质量改变时的误差分析
由图可知，嫦娥三号的初始质量的变化会引起嫦娥三号的质量和燃料消耗的变化，初始质量越大嫦娥三号的质量变化越大，燃料消耗的越多。
对主减速阶段的初始速度和初始高度进行误差分析，
嫦娥三号的预定着陆点
海拔为-2641m，则将主减速阶段的高度设置为15Km至17.641Km之间。将其与原有状态下的运动状态相互比较。仅考虑切向速度变化，根据燃料最优制导模型的计算方法，利用四阶龙格-库塔公式和拟牛顿法将主减速的30个阶段嫦娥三号偏角的变化与原变化进行比较，如下图：
图19 偏角的变化
上图蓝线表示原的变化，绿线为改变切向速度时的变化，红线为两者的误差，可以看出前期原偏角大于改变后的偏角，后期则相反。误差也随着时间变少。由误差计算公式 ，计算偏角总误差为-9.49% 。
根据已求得的偏角的的值，将主减速段运动路径分割为30个阶段，并将轨道离散化
图20 初始高度变化时轨道的变化
图21 初始高度变化的轨道离散化
图18的红线为原高度时轨道变化，粉红线为改变原高度时的轨道变化。由误差公式可得，在主减速阶段的误差为,误差率为 。
已知嫦娥三号的初始比冲量为2940，将其分别乘以0.9、1.1，即改变比冲量，观察嫦娥三号质量和燃料消耗的变化。
图22 比冲量变化时轨道特性的变化
由图可以看出，比冲量的值越大，嫦娥三号的质量变化越大，燃料消耗越大。
5.4.2 敏感度分析
粗避障段：粗避障段的范围是距离月面2.4km到100m区间，其主要是要求避开大的陨石坑，实现在设计着陆点上方100m处悬停，并初步确定落月地点。 将附件所给图片网格化为2300×2300的矩阵，本文根据处的月球高度，得到避障原则： （19）
使用matlab软件并采用用蒙特卡洛的方法进行1000次仿真（源程序见附录），模拟分析月面不同地形高度对嫦娥三号降落时所需调整概率大小的影响。
图23 粗避障阶段地形 图24 粗避障阶段不同降落高度所需调整 精避障段：精细避障段的区间是距离月面100m到30m。要求嫦娥三号悬停在距离月面100m处，对着陆点附近区域100m范围内拍摄图像，并获得三维数字高程图。分析三维数字高程图，避开较大的陨石坑，确定最佳着陆地点，实现在着陆点上方30m处水平方向速度为0m/s。
与粗避障一样，在满足同样的避障原则下，分析月面不同降落地形高度对嫦娥三号降落时所需调整概率大小的影响。
图25 精避障阶段地形 图26 精避障阶段不同降落高度所需调整
在精确避障阶段的避障原则下，为了研究嫦娥三号在降落时占地面积大小对轨道降落的敏感度的影响，选择和这两个数据作为嫦娥三号降落时的占地面积。并用这两个数据在matlab软件中做50次模拟比较（源程序见附录）。
图27 两种不同着陆占地面积着陆比较
由图可以看出，的占地面积的非调整降落次数高于。，由此推测着陆占地面积越大，可直接顺利着陆的概率越小。
分别选用六个不同的占地面积，对其进行1000次模拟，计算出1000次模拟中无需调整即可顺利着陆的次数，如下表：
表29 不同占地面积无需调整即可顺利着陆的次数
的概率越大。
6.模型的评价与改进方向
6.1模型的评价
6.1.1模型一的评价
模型一分别以着陆点的经度和纬度作为准备着陆轨道，选取经度不变的轨道处于稳定状态，不需要产生推力，此种轨道保证了燃料消耗的最优。选取纬度不变的轨道则保证了当平抛距离较大且难以精确确定时以最大概率降落在着陆区域内。
6.1.2 模型三的评价
在模型三中建立的主减速阶段燃料最优目标规划函数，利用时间逼近法快速求解月球最优软着陆问题。对于终端时间固定型最优控制问题,将其直接离散化为非线性规划问题,采用拟牛顿法和四阶Admas预测-校正积分方法快速求解。此方法优化精度较高,收敛速度快，比近年较为流行的智能算法(如遗传算法等)，减少了计算量且更符合实际需求和精度要求。
对嫦娥三号软着陆的其他阶段也分别建立了动力学模型。并且分段建立了最优目标函数。确定了着陆轨道。在精避障阶段，综合考虑了着陆位置的总高度和总平坦度，对不同数量级的数据标准化，设置偏好系数后对所有点进行全局搜索，得出了最优降落策略和最优降落点。
6.1.3 模型四的评价
模型四对误差进行了多方面的分析，包括对的最大推力、初始速度的变化的轨道特性分析、进一步对主减速阶段的偏角的趋势分析、计算误差和灵敏度。
6.1.4 模型的不足与改进方向
模型的不足：由于轨道的复杂多变性，本文简化了模型的数学推导，将高度变化引起的轨道路径长度变化忽略，只重点考虑和计算了主减速阶段的轨道特性，造成了设计的轨道系统的误差。且未对轨道路径做出明显的全局优化。
改进方向：
1. 将着陆轨道的六个阶段燃料量作为规划函数，将自适应遗传算法与模拟退火
算法相结合,形成一种自适应模拟退火遗传算法，增强轨道路径设计的整体搜索能力。
2. 根据月球岩石和坑的特征, 设计了粗障碍识别和安全着陆区选取算法: 1)
图像直方图分析; 2)K 均值聚类; 3) 过亮障碍识别; 4) 过暗障碍识别; 5)
纹理障碍识别; 6) 采用螺旋搜索算法确定每个单元格的安全半径(图3); 7) 根据安全半径, 选取候选安全着陆点; 8) 评估候选安全着陆点避障所需的速度增量; 9) 根据安全半径和速度增量评价值, 综合确定安全着陆点.
7.参考文献
[1] 张德丰，MATLAB数值分析，北京：机械工业出版社，2012。
[2] 赵吉松，谷良贤，高原，月球软着陆轨道的时间逼近法快速优化设计[J]，宇航学报，第29卷第5期：1-5，2008. 9。
[3] 朱建丰，徐世杰，基于自适应模拟退火遗传算法的月球软着陆轨道优化[J]，航空学报，第28卷第4期：2-3，2007.7。
[4] 于彦波，火星探测器动力下降段制导律研究[D]，哈尔滨，哈尔滨工业大学， 2013。
[5] 张仲满，月球软着陆的制导算法研究[D]，哈尔滨，哈尔滨工业大学， 2009。
[6]田青，常微分方程初值问题数值解的实现与分析，http://www.doc88.com/p-[1**********]3.html,2014.09.13。
[7]张洪华，梁俊等，嫦娥三号自主避障软着陆控制技术[J]，中国科学：科学计算，第44卷第6期：2-4,2014。
附录
蒙特卡洛分析不同降落地形高度选取调整概率
I2=imread ('data.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
c=zeros(200,1);
d=zeros(200,1);
for j=1:1:200
for i=1:1000
aa=floor(2299*rand(1)+1);
bb=floor(2299*rand(1)+1);
if pp2(aa,bb)
c(j,1)=c(j,1)+1;
else
d(j,1)=d(j,1)+1;
end
end
end
plot(1:200,c/1000,1:200,d/1000);
title('不同降落地形高度选取调整概率');
xlabel('地形高度');
ylabel('调整与非调整概率');
legend('不需调整','需要调整');
精避障阶段50次模拟分析嫦娥安好占地面积大小对着陆的影响： d1=zeros(50,1);
for j=1:50
for i=1:1000
aa=floor(979*rand(1)+11);
bb=floor(979*rand(1)+11);
aa1=aa-7;
aa2=aa+7;
bb1=bb-7;
bb2=bb+7;
c=0;
for cc=aa1:aa2
for dd=bb1:bb2
if pp2(cc,dd)>=50&pp2(cc,dd)
c=c+1;
end
if c==15*15;
d(j,1)=d(j,1)+1;
end
end
end
end
end
plot(1:50,d)
快速调整阶段运动状态
v1=-90;
a=(7500-1200*1.6)/1200;
d=0;
for t=0.1:0.1:14
d=d+1;
x(1,d)=3000+v1*t+0.5*a*t^2;
end
plot(0.1:0.1:14,x')
title('快速调整段制导律下运动状态');
xlabel('时间/t');
ylabel('高度/m');
总燃料消耗量与时间的关系
F=7500;
m=2400;
ve=2940;
v=22;
g=1.63;
a=3.125;
m1=zeros(1,487);
for i=1:487
m1=F/ve;
m=m-m1-m*g/ve;
A(1,i)=m;
F=m*a;
end
subplot(1,2,1)
plot(1:487,A)
title('速度');
xlabel('时间/t');
ylabel('嫦娥卫星质量/Kg');
subplot(1,2,2)
plot(1:486,A(2:487)-A(1:486))
title('能料消耗的变化');
xlabel('时间/t');
ylabel('质量/Kg');
快速调整阶段最优降落点
I2=imread ('data1.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
a=input('请输入嫦娥号的大小a=');
for aa=a+1:1:1000-a
for bb=a+1:1:1000-a
AA=pp2(aa-a+1:aa+a,bb-a+1:bb+a);
A(aa-a,bb-a)=sum(sum(AA));
if A(aa-a,bb-a)==2389
aa
bb
end
end
end
初始高度变化时轨道的变化
H1=15000;
H2=15000+2641;
X=0;
X2=0;
i=0:1:510;
I=0:1:600;
X=1700.*i-0.5*3.125*(i.^2);
X2=1700.*I-0.5*2.878*(I.^2);
Y=H1-22.*i;
Y7=min(H2-22*1.1*I);
Y5=H2-22.*I*1.1;
i1=510:559;
I1=600:640;
X1=57*(i1-510)+460590-0.5*2.*((i1-510).^2);
Y1=min(Y)-22.*(i1-510);
Y8=Y7-22*1.1.*(I1-600);
X3=460982:0.1:460990;
Y3=0:2700/80:2700;
Y_1=0:100:3500;
X_1=linspace(501967,500720,36)
plot(X,Y,'r',X2,Y5,'m',X1,Y1,'k',X3,Y3,'g',X_1,Y_1); title('初始高度变化时轨道的变化');
xlabel('和近月点的距离');
ylabel('距地面高度/m');
text(250000,11000,'主减速阶段');
text(420000,3500,'快速调整阶段');
text(380000,2000,'粗避障阶段');
text(380000,1000,'其他阶段');
网格化矩阵方差
I2=imread ('data.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
A=zeros(100,100);
for m=1:22
for n=1:22
for i=(1:100)+m*100
for j=(1:100)+m*100
A((i-100)*m,(j-100)*n)=pp2(i,j);
c=var(A(:))
end
end
end
end
搜索可选择的着陆点
I2=imread ('data1.tif');
p2=I2;
[y,x]=size(p2);
[X,Y]=meshgrid(1:x,1:y);
pp2=double(p2);
a=zeros(1000,1000);
for i=1:1000
for j=1:1000
if pp2(i,j)>=60&pp2(i,j)
a(i,j)=1;
else
a(i,j)=0;
end
end
end
bb=input('输入bb=');
for m1=bb+1:999-bb
for m2=bb+1:999-bb
if sum(sum(a(m1-bb:m1+bb,m2-bb:m2+bb)))==(2*bb+1)^2 m1
m2
else
continue;
end
end
end