嫦娥三号软着陆轨道设计与控制策略建模论文

2014年高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 2009 年 9 月 12 日赛区评阅编号（由赛区组委会评阅前进行编号）：2009高教社杯全国大学生数学建模培训竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）： 嫦娥三号软着陆轨道设计与控制策略 摘 要本文根据题目的要求建立了合理的嫦娥三号软着陆轨道设计与控制策略模型模型。，我们借助多种数学软件的优势挖掘出大量数据潜在的信息，并将其合理运用，在此基础上，以最优控制策略为最大目标，长远发展为原则，制定出信息不足条件下的量化综合评价体系。 在本文所建立的模型中，我们采取了层次分析法（AHP）、数据统计拟合以及整数线性规划相结合的手段，这样既借鉴了层次分析法综合评价的优势，又克服了该法中主观因素的不确定性，使模型更具有科学性，要确定着陆准备轨道近月点和远日点的位置，以及嫦娥三号相应速度的大小和方向。考虑了月球自转，针对三维空间内精确定点软着陆问题利用参数 化控制解决了变推力软着陆最优控制问题，此外还针对仅知制动初始点到月心距离而具体位置未知的 情况，对初始点（近月点）的选取进行了研究。 关键字： 嫦娥三号 着陆点 最优控制策略 科氏定律 矩阵 1、问题重述1.1 背景资料根据计划，嫦娥三号将在北京时间12月14号在月球表面实施软着陆。嫦娥三号如何实现软着陆以及能否成功成为外界关注焦点。目前，全球仅有美国、前苏联成功实施了13次无人月球表面软着陆。 目前，全球仅有美国、前苏联成功实施了13次无人月球表面软着陆。北京时间12月10日晚，嫦娥三号已经成功降轨进入预定的月面着陆准备轨道，这是嫦娥三号“落月”前最后一次轨道调整。在实施软着陆之前，嫦娥三号还将在这条近月点高度约15公里、远月点高度约100公里的椭圆轨道上继续飞行。期间，将稳定飞行姿态，对着陆敏感器、着陆数据等再次确认，并对软着陆的起始高度、速度、时间点做最后准备。由于月球上没有大气，嫦娥三号无法依靠降落伞着陆，只能靠变推力发动机，才能完成中途修正、近月制动、动力下降、悬停段等软着陆任务。据了解，嫦娥三号主发动机是目前中国航天器上最大推力的发动机，能够产生从1500牛到7500牛的可调节推力，进而对嫦娥三号实现精准控制。在整个“落月”过程中，“动力下降”被业内形容为最惊心动魄的环节。问题一主要是确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。本文同样考虑了月球自 转，针对三维空间内精确定点软着陆问题利用参数 化控制解决了变推力软着陆最优控制问题，此外还 针对仅知制动初始点到月心距离而具体位置未知的情况，对初始点的选取进行了研究。问题二：确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。问题三：对于你们设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。 3、模型假设 1、假设卫星或飞船相对于地球极小可以看做质点2、假设地球是个规则球体，质量集中于地心3、假设外界引力对该系统可忽略不计4、忽略影响测控站布置的地理因素5、不考虑测控站周围地理环境和天气环境对嫦娥三号测控的影响 4、符号说明 坐标原点，代表月心 指向月球赤道相 对于白道的升交点 指向月球自转角速度方向 按右手坐标系确定 月固坐标系 以月球赤道面为参考平面 指向赤道面与起始子午面的交线方向 指向月球自转角速度方向 按右手坐标系确定 指原点在嫦娥三号质心的轨道坐标系 指向从月心到 着陆器的延伸线方向 垂直指向运动方向 按右手坐标系确定 制动发动机推力 与轴的夹角 在平面上的投影与轴负向所成夹角 为与所成夹角 为在平面上的投影与轴正向所成夹角 为月球自转而产生的月固坐标系相对惯性坐标系的转角 嫦娥三号速度矢量在月固坐标系轴上的投影 嫦娥三号速度矢量在月固坐标系轴上的投影 嫦娥三号速度矢量在月固坐标系轴上的投影 嫦娥三号发动机的推力 嫦娥三号的质量 为某一高度月球重力加速度在月固定系轴上的投影 为某一高度月球重力加速度在月固定系轴上的投影 为某一高度月球重力加速度在月固定系轴上的投影 月球自转角速度 月球引力常量 嫦娥三号制动时的比冲，是一个常值 预定着陆点在月固坐标系中轴的坐标 预定着陆点在月固坐标系中轴的坐标 预定着陆点在月固坐标系中轴的坐标 着陆点到月心距离，即月球半径。 调节参数 控制变量 5、模型建立一、软着陆轨道模型建立 探月飞行器嫦娥三号首先进行霍曼变轨，从圆形环月轨 道进入一条近月点高度为15 km的椭圆轨道；当到 达近月点时，制动发动机点火，探测器进人动力下降 段，最终以很小的相对速度(小于6 ms)降落到月 面指定位置。 如图1所示，定义惯性坐标系，原点在月心，参考平面是月球赤道面，轴指向月球赤道相 对于白道的升交点，轴指向月球自转角速度方 向，轴按右手坐标系确定。再定义月固坐标系，以月球赤道面为参考平面，指向赤道面与起始子午面的交线方向，指向月球自转角速度方向，轴按右手坐标系确定。为原点在嫦娥三号质心的轨道坐标系，指向从月心到 着陆器的延伸线方向，垂直指向运动方向，按右手坐标系确定。制动发动机推力的方向与嫦娥三号纵轴重合，为与轴的夹角，为在平面上的投影与轴负向所成夹角。为与所成夹角，为在平面上的投影与轴正向所成夹角。为月球自转而产生的月固坐标系相对惯性坐标系的转角，不妨假设初始时刻月固坐标系与惯性坐标系重合。 图一：坐标示意图显然有轨道坐标系到惯性坐标系转换矩阵 惯性坐标系到月固坐标系的转换矩阵为根据牛顿第二定律，结合科氏定律整理可以得到嫦娥三号在月固坐标系中的运动方程为其中，,为嫦娥三号速度矢量与月固坐标系各轴上的投影，为发动机的推力，为嫦娥三号的质量，,和为该高度月球重力加速度在月固定系各轴上的投影，为月球自转角速度。因此，在月固坐标系中嫦娥三号的运动方程可表示如下： 其中为月球引力常量，为嫦娥三号制动器的比冲，是一个常值。 取为系统状态变量，为控制变量，则式可以简记为 二、燃料节省最优模型建立 1、 按照耗燃最优的要求，取性能指标为 在实际情况下，通常没必要令嫦娥三号着陆速度严格等于零，只要能保证嫦娥三号以很小的相对速度 降落到月面就足可以接受的。因此，考虑到这一点， 本文将软着陆的末速度要求以惩罚冈子的形式加入 到指标中如下式所示，主要目的是降低最优控制问 题求解的复杂度，该惩罚因子可以通过反复的数值 仿真运算，按经验设定。此外，显然有约束条件 其中，为预定着陆点在月固坐标系中的坐标；为着陆点到月心距离，即月球半径。对于含有形如善。这类关于状态变量在连续时 间上都要满足的不等式约束最优化问题，至今还是 最优化领域的一个难点。文献10中给出一种约束变换技术，使得该类问题得到解决。 显然等价于 但上式显然在时不可微，因此用如下不等式去近似上式 其中 是调节参数。文献证明了当足够小的时候，存在使得对任何满足的能够令对达到满足要求的近似。不妨记为用式替换得到的新的约束函数。因此本文所讨论的软着陆耗燃最优问题转化为：在系统(1)满足约束函数G的情况下，求取适当的控制变量“使指标函数(2)达到最小。2 参数化控制求解耗燃最优问题 假定初始时刻为0，终端时刻为待定参数。选取满足的序列和三组参数，构造形如 的参数化分段常数控制器。其中 用控制器替换系统中的，则问题一转变为：问题2 寻找三组参数，来最小化指标函数，并且满足约束函数。 显然，对于每个给定的P，这都是一个有限维的参数优化问题。文献11中第六章已经证明了当时，问题2的最优解收敛于问题1的最优解。 不过文献12已经证明了在数值计算中，求解问题 2的参数梯度时难度很大甚至求不出真实解，因而 本文引入强化技术来解决这一问题。 从到构造如下变换 上式中,，序列为区间上预先给定的分段点，并且满足。将上式两边对求导可得 其中 不妨令 ，则得到如下增广系统 即 其中O,P,Q与分别为O,P,Q与经过变形后的形式 指标函数变为约束条件变为 其中，由于仅仅已知探测器在软着陆起始点到月心的距离和探测器的起始速度，原来质量，而软着陆起始点与另两个空间位置信息角的初始值未知，因而令为系统待定 参数。则系统可以表示为： 那么问题2转化为如下问题：在系统(6)满足约束并且初始条件如式(9)的情况下，求取适当的控制变量使指标函数(7)达到最小。再由 可知，问题3将最初的探月飞行器软着陆最优控制问题转化成了优化静态控制参数，和以及系统参数，的问题，利用经典的参数优化算法即可求出登月飞行器的软着陆最优控制的一组逼近解和软着陆最优初值点位置以及终端时刻。利用此算法，增加时间的分段点个数叮以重新优化，经过多次优化后即可得到满意精度的参数化解。 此外，假如令系统(1)中的推力F为已知的恒 定推力，令控制变量，则本文问题变为恒定推力下软着陆最优控制问题，依然可以利用本文方法解决，而依据极大值原理结合传统的打靶 法则只能解决恒定推力的情况，因I而相比之下本文方法适用性更广。3 数值仿真已知探测器初始质量；制动发动机最大推力为，最大推力,比冲;初始速度，；月球自转角速度；月球引力常数；近月点距月心距离月球半径。登月点选择月面上的雨海，位置为北纬，西经利用最优控制软件，通过计算机仿真运算，令，即可得到符合精度的最优解，最终利用本文的参数化控制得到软着陆末时刻，末时刻探测器质量，燃料消耗为，最后探测器以的对月速度精确降落到指定登月点。此为可得，从而有最优初始点坐标，。若不考虑对初始点位置的优化，文献8利用打靶 法最终得到着陆时探测器质量为，着陆位置 距预定着陆点，相比之下本文方法在燃料消 耗上节省了，同时落点精确，没有偏差。 图2利用本文方法得出的最优控制率，由于最优控制率是分段长值函数因而为梯形图。图3为三个方向上的速度曲线，因而可以看出探测器软着陆时相对月面速度足够小，软着陆成功实现。图4为软着陆最优轨线，显示了角以及探测器距离月心的距离随时间变化的曲线，图5是探测器质量变化曲线。 图三：软着陆速度曲线 图四：软着陆最优曲线 图五：质量变化曲线 3、问题假设 5、模型的建立及求解5.1 问题一模型的建立及求解5.1.1 模型一：假设卫星或飞船运动轨道为圆在不考虑地球自转的条件下，地球自转时该卫星或飞船在运行过程中相继两圈的经度的差异可不予考虑。卫星或飞船从起飞时加速升空后经一系列加速变轨，最终的运行轨迹是一圆周。即最终卫星或飞船绕地球做匀速圆周运动。用卫星或飞船的运动轨迹所在的平面去切地球会得到一圆面。如图l所示：图3 观测站对圆形卫星轨道覆盖范围示意图地球卫星轨道我们只需在如图C点建立一测控站即可测控A至B之间的劣弧区域，最小测控站数目即为需要覆盖卫星轨道的这样的C点的个数，利用正弦定理解三角形按照此模型以神州七号飞船为例：地球半径为6400公里，飞船进入预定轨道运行稳定后距地球表面高度为343公里，相关数据代入，运用MATLAB计算得出，n=12，即此时需要最少测控站的数目为12个。5.1.2 模型二：考虑到实际，按卫星或飞船运动轨道为椭圆由于在实际情况中飞船的运行轨道为椭圆形，如图2或下图，取椭圆近地点旁边的焦点为地球的圆心，椭圆轨道定位很麻烦，因此先估算，然后再精算1、以地心为圆心，地球半径与近地点之和为半径作圆，如图4、由于圆包含在椭圆区域之内，若能监控到圆周及以外空域，则定能监控到椭圆及以外空域，因此，在地球上均匀建站监控整个圆周。图4 观测站对椭圆形卫星轨道覆盖范围示意图1地心具体算法为：其中是如图所示的圆心角，角如图所示，以神舟七号为例，近地点高度=200公里，所以，用MATLAB软件解得 n=16 。2、以地心为圆心，地球半径与远地点之和为半径作圆，如图5所示。图5 观测站对椭圆形卫星轨道覆盖范围示意图2焦点地心由于大圆包含了椭圆区域，因此只要监控到大圆周及以外空域，则未必能监控整个椭圆周。在地球上均匀建站监控整个圆周，其算法和1）中相同：，以神舟七号为例，其远地点公里，令，解得 n=12。综上，椭圆轨道上的监控站应该在12至16个之间。下面我们进行椭圆轨道的监控站数精确计算。1、椭圆轨道的测控算法思想用逼近方法和迭代算法来实现。由于对椭圆监控不能像对待圆一样均匀分布站点，而，因此分布测控站是一件很棘手的事情。可用前面模型中的算法，可借用的算法来算出的近似值。根据物理学、力学知识及开普勒三大定律可得到右侧的向径，a为椭圆的长半轴，e为离心率，f为向径与椭圆水平长轴的夹角，r为球心到椭圆上的向径。所以椭圆轨道的监控算法为：，用逼近的方法可近似计算出。5.2 对于问题二模型的建立与求解5.2.1、模型一：由于一个卫星或飞船的轨道与地球赤道平面有固定夹角，而求地球自转时该卫星或飞船在运行过程中相继两圈的经度有一些差异，因此，对卫星在地球表面的星下点轨迹进行分析。1、当卫星运行角度与地球自转角速度相同时，卫星沿运行轨道运行一周后星下点轨迹又回到起点，星下点轨迹见图5，但其相继两圈的经度无变化，不合题意。2、当卫星运行角速度是地球自转 速度的两倍是卫星沿运行古道运行两圈后星下点轨迹回到起点。3、当卫星运行角速度是地球自转速度的三倍时，卫星沿运行轨道运行三圈星下点轨迹回到起点。图6 卫星或飞船地面轨迹 卫星运行角速度越大，卫星在地球表面上留下的星下点轨迹越密，从而形成了卫星或飞船地面轨迹，如图6所示。由以上图形观察分析可以得出：星下点轨迹均匀地分布在赤道的两边，即北纬a与南纬a之间，因此，在卫星或飞船的星下点轨迹较为简单的情况下可沿着星下点轨迹设立测控站，对于一般情况，特别是卫星或飞船的星下点较密时，只需测控北纬a与南纬a之间的区域，就可以实现对该卫星的全程监控。在实际测控范围与卫星轨道面的相交圆D中取圆内接正方形并用内接正方形覆盖所要测控的区域。首先根据轨道面与赤道面的夹角a以及每个内接正方形边长在地心所对的圆心角，可求出测控区域中正方形的行数j，如图所示；其次根据每一行正方形覆盖的轨道面的圆周长即可确定每一行所需的正方形的个数i，从而计算出全程监控所需要的监控站的最少个数，如图7所示。图7 圆内接正方形覆盖法模型求解过程：1、计算正方形的边长在地心所对的圆心角为，有由于卫星围绕地球运动时，轨道与赤道平面有一个夹角a，测控范围由正方形覆盖，所以正方形的行数：圆内接正方形每条边在卫星轨道面上所对应的圆弧长为：2、当j为偶数时，各行正方形在南北纬之间对称分布，每一行正方形的总长度是正方形所在轨道面的圆的周长，北纬第一行正方形覆盖的轨道面圆周长，北纬第一行排列的正方形个数，第二行正方形覆盖的轨道面圆周长：，第二排排列的正方形个数，依次类推：北纬第k行排列的正方形个数，为实现全程测控需要的测控站总数当j为奇数时，先用一行正方形覆盖赤道，使正方形的中心在赤道上，然后将各行正方形在南北半球对称分布，所以赤道上排列的正方形的个数为，北纬第一行正方形覆盖的轨道面圆周长为，北纬第一行排列的正方形个数，同理可以求出第k行上分布的测控站的个数，实现全程测控的测控站的总数根据神舟七号的运行数据，将轨道与赤道面夹角42.2度，公里代入得j=4，n=60，即如果要实现对神舟七号的全程监控至少需要60个测控站对其进行测控。5.2.2、问题二中，为使计算方便，采用内接正方形覆盖所要测控的区域.，而圆内接正六边形的面积占圆面积的82.74%，因此可以考虑采用圆内接正六边形覆盖，以提高有效覆盖率，我们因此建立了模型二：设卫星或飞船绕地球飞行的倾角为，离地面的高度为,地球半径为R，卫星或飞船飞行的轨迹为一环绕地球半径为的球面去掉两端球帽的曲面，通过公式计算出球帽的面积得到卫星或者飞船飞行的曲面面积，每个测控站测控的范围与曲面的交线为一个半径为r的圆，只有圆与圆之间有重叠才能够保证测控站的测控范围覆盖所有曲面，这样我们就取每个圆的内接正六边形作为每个测控站对卫星或飞船运行曲面的平均有效测控范围，通过圆内接六边形与圆的面积比率计算出圆内接正六边形的面积，用曲面的面积去除以正六边形的面积，就可以得到需要的最少测控站数目，具体计算方法如下： 球帽面积公式： 球面面积计算公式： 卫星运行曲面的面积： 圆内接正六边形与圆的面积之比:测控站测控的范围与曲面交线圆的面积：对应的圆内接正六边形的面积：其中测控站测控的范围与曲面的交线圆半径： 需要的测控站数目至少为：将神舟七号运行数据：离地面高度H=343千米，地球半径R=6371千米以及倾角代入上式，编辑MATLAB程序（见附件4）运算得到n=67，即最少需要建67个测控站才能全程测控神七飞行。图8 圆内接六边形覆盖法5.3 对问题三资料的查取与分析通过上网查阅资料，搜集了我国神州7号飞船运行资料和相关测控站点的分部信息，飞船运行在轨道倾角42.4度、近地点高度200公里、远地点高度350公里的椭圆轨道上，实施变轨后，进入343公里的圆轨道。各个测控站的地理位置位置如表1。表1 神舟七号测控站地理位置卫星测控站经度纬度1东风站98.50E39.70N2酒泉卫星发射中心100.50E41.83N3喀什卫星测控站75.99E39.45N4和田站79.92E37.10N5青岛站120.30E36.20N6陕西渭南卫星测控站109.50E34.50N7厦门卫星测控站118.01E24.07N8纳米比亚站14.52E22.67N9卡拉奇卫星测控站66.99E24.82N10马林迪站40.10E3.22S11圣地亚哥站70.10E33.43S12远望一号测控船77.0W20.0S13远望二号测控船150.0W31S14远望三号测控船0.0E1.0N15远望四号测控船135.0E30.0N16远望五号测控船180.0E32.0N根据世界地图及神舟七号测控通信系统分布图9所示。图9 神舟七号测控站分布图我们对其中的测控站进行拟合，发现所成图形近似正弦函数的曲线。对此。我们试作以下分析：卫星的运行轨道在地球表面上的投影应大致服从正弦曲线分布。测控站的作用是为了能够更清晰更准确的了解卫星运行的具体情况，故测控站分布在卫星的运行轨道在地球表面上的投影曲线上是十分必要的。查找了神舟六号的观测站数据，各观测站具体位置见附件表1，对模型进行了验证，通过对数据的分析显示，监控占的位置同样接近于正弦曲线，增大了对卫星检测的覆盖面。6、模型的评价及其改进6.1、模型的评价优点：模型有特殊到一般，由简单到复杂，充分考虑了卫星或飞船运行轨道为椭圆时的测算方式，灵活巧妙地计算出了需要最少测控站的区间，减小了数据误差。问题二中我们考虑正四边形和正六边形两种情况，使模型具有更高的可靠性对于问题一，我们建立了圆形轨道模型。由于在模型建立之初所作的假设较少，故与真实情况比较接近，因而实用性强，同时把圆形轨道模型推广到椭圆模型，对问题2的解决起着重要作用。在解决问题2的过程中，我们建立了圆锥测控模型。并且应用此模型对神舟六号与神舟七号飞船的测控站的数目进行了验证，所得结果与实际相近。这充分说明了所建模型比较理想。对于问题3，我们给出了测控站分布的拟合曲线，发现卫星的运行曲线大致服从正弦曲线，如图所示。另外，在某些具体问题上，采用内接正方形或内接正六边形、正八边形覆盖所要测控的区域，会出现较多的浪费，这时不能一味追求圆内接图形面积占圆面积的比率，而应该从整体的角度来考虑测控站的最优分布。6.2、模型的改进关于问题二：在地球自转的影响下卫星运行过程中星下线的轨迹是地球表面的一些曲线，计算测控站的数量比较困难。一种估算方法是设置许多测控站，使得其能覆盖卫星飞过的所有空域。计算这个涵盖赤道的球面的立体角，再用一个观测站所能覆盖的立体角去除这个涵盖赤道的立体角，就可以得到得到要覆盖这个区域至少需要的观测站的数目7、参考文献1 姜启源，数学模型M，北京高等教育出版社，19932 叶其孝，数学建模，教育与国际数学建模竞赛中国工业与应用数学学会工科数学杂志社，19943 萧树铁主编，大学数学试验（第二版）M，高等教育出版社，20064 苏金明 阮沈勇，Matlab 实用教程，北京，电子工业出版社，20055 张海基等，卫星和飞船跟踪测控的数学模型，苏州市职业大学学报，2010年第1期，72-756 戴红兵等，卫星和飞船的跟踪测控，思茅师范高等专科学校学报，2009年第6期，28-327 杨徐昕等，卫星或飞船的跟踪测控模型设计，重庆科技学院学报（自然科学版），2010年第1期，183-1878 柳仲贵，卫星跟踪系统的动态范围，飞行器测控学报，2003附件：附件1：卫星或飞船与地球、太阳关系Matlab程序figure(name,卫星或飞船与地球、太阳关系); s1=0:.01:2*pi;hold on;axis equal;%建立坐标系axis off % 除掉Axesr1=10;%地球到太阳的平均距离r2=3;%卫星或飞船到地球的平均距离w1=1;%设置地球公转角速度w2=12%设置卫星或飞船绕地球公转角速度t=0;%初始时刻为pausetime=.002;%设置暂停时间sita1=0;sita2=0;%设置开始它们都在水平线上set(gcf,doublebuffer,on) %消除抖动plot(-20,18,color,r,marker,.,markersize,40);text(-17,18,太阳);%对太阳进行标识p1=plot(-20,16,color,b,marker,.,markersize,20);text(-17,16,地球);%对地球进行标识p1=plot(-20,14,color,k,marker,.,markersize,13);text(-17,14,卫星或飞船);%对卫星或飞船进行标识plot(0,0,color,r,marker,.,markersize,60);%画太阳plot(r1*cos(s1),r1*sin(s1);%画地球公转轨道set(gca,xlim,-20 20,ylim,-20 20);p1=plot(r1*cos(sita1),r1*sin(sita1),color,b,marker,.,markersize,30);%画地球初始位置l1=plot(r1*cos(sita1)+r2*cos(s1),r1*sin(sita1)+r2*sin(s1);%画卫星或飞船绕地球公转轨道p2x=r1*cos(sita1)+r2*cos(sita2);p2y=r1*sin(sita1)+r2*sin(sita2);p2=plot(p2x,p2y,k,marker,.,markersize,20);%画卫星或飞船的初始位置orbit=line(xdata,p2x,ydata,p2y,color,r);%画卫星或飞船的运动轨迹while 1set(p1,xdata,r1*cos(sita1),ydata,r1*sin(sita1);%设置地球的运动过程set(l1,xdata,r1*cos(sita1)+r2*cos(s1),ydata,r1*sin(sita1)+r2*sin(s1);%设置卫星或飞船绕地球的公转轨道的运动过程ptempx=r1*cos(sita1)+r2*cos(sita2);ptempy=r1*sin(sita1)+r2*sin(sita2);set(p2,xdata,ptempx,ydata,ptempy);%设置卫星或飞船的运动过程p2x=p2x ptempx;p2y=p2y ptempy;set(orbit,xdata,p2x,ydata,p2y);%设置卫星或飞船运动轨迹的显示过程sita1=sita1+w1*pausetime;%地球相对太阳转过的角度sita2=sita2+w2*pausetime;%卫星或飞船相对地球转过的角度pause(pausetime); drawnowend附件2：卫星或飞船地面轨迹图Matlab程序load(topo.mat,topo,topomap1);whos topo topomap1contour(0:359,-89:90,topo,0 0,b)axis equalbox onset(gca,XLim,0 360,YLim,-90 90, . XTick,0 60 120 180 240 300 360, .Ytick,-90 -60 -30 0 30 60 90);hold onx=linspace(0,500*pi,1000);y1=58*sin(0.05*x+0.25*pi);y2=58*sin(0.05*x+0.5*pi);y3=58*sin(0.05*x+0.75*pi);y4=58*sin(0.05*x+pi);y5=58*sin(0.05*x+1.25*pi);y6=58*sin(0.05*x+1.5*pi);y7=58*sin(0.05*x+1.75*pi);y8=58*sin(0.05*x+2*pi);plot(x,y1,k,x,y2,g,x,y3,r,x,y4,y,x,y5,b,x,y6,g,x,y7,m,x,y8,r)hold onz1=0*x+58;z2=-z1;plot(x,z1,k,x,z2,k)附件3：圆内接六边形全覆盖图Matlab程序load iris_datasetnet = newsom(irisInputs,16 3);plotsomtop(net);hold onx=-1:0.05:16.5;y1=0.85+1.15*sin(x);y2=0.85+1.15*sin(x+0.25*pi);y3=0.85+1.15*sin(x+0.5*pi);y4=0.85+1.15*sin(x+0.75*pi);y5=0.85+1.15*sin(x+pi);y6=0.85+1.15*sin(x+1.25*pi);y7=0.85+1.15*sin(x+1.5*pi);y8=0.85+1.15*sin(x+1.75*pi);z1=0*x+2;z2=0*x-0.3;plot(x,y1,k,x,y2,g,x,y3,r,x,y4,y,x,y5,b,x,y6,g,x,y7,m,x,y8,r,x,z1,k,x,z2,k)附件4：function y=cekongzhan(R,H,b) b=b./180.*pi;s1=2.*pi.*(R+H).*(1-sin(b);s=4.*pi.*(R+H).2;s2=s-s1;e=0.827;r=(R+H).*sin(87./180.*pi-asin(R.*(sin(93./180.*pi)./(R+H);s3=pi.*r.2;s4=s3.*e; y=s2./s4;y=ceil(y);end 附件表1 神舟七号测控站地理位置卫星观测站经度纬度1北京航天指挥控制中心120.25E48.07N2酒泉卫星发射中心100.50E41.83N3西安卫星控制中心108.91E34.23N4山西兴县卫星测控站111.22E38.47N5陕西渭南卫星测控站109.50E34.50N6厦门卫星测控站118.01E24.07N7喀什卫星测控站75.99E39.45N8远望一号测控船131.83E37.23N9远望三号测控船13.72W09.50N10卡拉奇卫星测控站66.99E24.82N11和田卫星测控站79.94E37.12N12南非卫星测控站18.42E33.92N13远望二号测控船75.84W19.84N14愿望四号测控船151.00E33.00S28