数学建模嫦娥三号运行轨迹及着陆点分析

2014高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：

嫦娥三号软着陆轨道设计与控制策略摘要：根据题目附录和文献[4]中提供的嫦娥三号的运行参数，利用角动量守恒及向量几何的方法，分别确定了近日点、远日点的位置向量和速度向量。与文献[4]的真实数据比较发现吻合良好。本文重点关注优化减速控制与着陆点避障两方面：前者燃耗最大，后者决定着陆成败。首先，在多重坐标变换基础上，建立了飞行器制动的动力学方程。并以燃耗为最优化性能指标、近月点状态为初始条件、着陆点状态为终端条件，利用极值原理求解飞行器的着陆轨迹，及其最优控制参数。其次，对避障阶段采集的高程图采取水平剖分、比较高程方差的方法，解出最优降落点。关键词：软着陆；最优轨道；避障问题重述嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道，于北京时间12月14号在月球表面实施软着陆。嫦娥三号在着陆准备轨道上的运行质量为2.4t，安装在其下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。嫦娥三号四周安装了姿态调整的发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m。嫦娥三号在高速飞行的情况下，为了保证嫦娥三号能准确地在月球预定区域内实现软着陆，关键的问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求如下：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段，要求满足每个阶段在关键点所处的状态；尽量减少软着陆过程的燃料消耗。

根据上述的基本要求，建立数学模型解决下面的问题：

（1）计算其着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。

（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。（3）对于设计的着陆轨道和控制策略进行相应的误差分析和敏感性分析。2、问题分析2.1技术背景月球软着陆的方式有两种，分别是：1）从月球转移轨道直至月球轨道，从而实现软着陆；2）从月球的停泊轨道变轨到近月点，然后实现软着陆。[2]嫦娥三号要求准确的在月球预定区域内实现软着陆，它降落的过程为如下几个阶段：嫦娥三号先进行霍曼变轨，由高度约100km的停泊轨道进入到近月点高度约15km的椭圆轨道；到达近月点时，探测器的制动发动机点火，进入动力下降阶段，即上图中的DA段；在距离月面约2km时，水平速度减为0，即上图中的AF段，调整姿态后，探测器以自由落体的方式降落到月面，即上图中的FL段。2.2着陆过程的主要考虑因素着陆期间主要考虑两大因素：一、主减速段是用时最长、推进剂消耗最多，因此最小化推进剂消耗是该段制导律的主要设计目标。二、避障阶段决定飞行器能否顺利着陆，因此统计分析高程数据并采取相应横向推力最关重要。2.3力学分析如上图示，以月心为坐标原点，月球自传轴为OZ轴，OX轴过零经度线，利用右手准则确定OY轴，建立空间直角坐标系。嫦娥三号探测器绕月运动中，受到指向月心的引力F引力，以及相背于速度方向的由发动机产生的推力F引力是在月心坐标下考虑。F其参数为：G为万有引力常量且G=6.67×10M为月球的质量，M=7.349×10m为航天器的质量r=（发动机推力是在飞行器轨道坐标系下考虑的：FVe为发动机的冲力，是由发动机构造、性质、燃料共同决定。飞行过程中是不变的。mdmdtF是推力方向，由飞行器的飞行姿态决定。F因此运动的加速度为：d2.4坐标系2.4.1直角坐标系的旋转我们假设月球是一个不旋球体，不考虑侧向运动，我们建立的月球探测器的数学模型，利用燃耗最优原则，依据庞特里亚金最大值原理设计了最优轨道。首先，我们建立了以OXYZ为原点的月心惯性坐标系，OX轴是月球的0经度线；OZ是月球的自转方向；OY是根据右手坐标来确立的。要在月心惯性坐标下计算推力，需进行坐标旋转变换。a根据欧拉姿态角的变化过程，OXYZ经过两次有顺序的坐标轴旋转，得到两个矩阵：先将轨道坐标OX1Y1ZM再将之沿绕OX方向旋转角度-βM因此得到总的旋转矩阵：TAX1Y1Z1为原点在探测器的轨道坐标系,AZ1指向从月心到着陆器的延伸线方向，AY1垂直指向运动方向,AX1按右手坐标系确定。制动发动机推力F的方向与探测器纵轴重合，θ为在月心惯性坐标中，推力为：

FF是推力方向，如上图所示。AZϑ推力方向与AZφ是推力在AX1YΘ和φ随飞行器飞行状态而变，那么：F同理，在月心惯性坐标下，总的加速度为：a2.4.2直角坐标系与经纬坐标系的变换如上图所示，以月心为坐标原点，月球自转轴为竖轴建立空间直角坐标系OXYZ，假设近月点的空间直角坐标为(x,y,z)，用μ和σ分别表示近月点位置的经度和纬度，可由以下经纬度转换成为直角坐标系坐标的公式：Z=r∙sin3.模型假设不考虑月球的自转。因为飞行器的速度相比月球自转的速度快、着陆时间快得多。不考虑除月球以外其他天体引力对飞行器的扰动。不考虑着陆区域附近因地形不同而造成的质量不同，进而影响引力场。即假设引力场在同样高度是均匀的。假设月球引力场是从月心点发出的，引力指向月心点。假设飞行器燃料足够使用。4、符号及字母描述5.近月点与远月点参数的确定5.1近月点的参数根据文献资料数据[4]，由北京工程研究所和空间智能控制技术国家级重点实验室的研究员们在这篇文献中对嫦娥三号的发射与着陆的设计做出了详细的说明,参与嫦娥三号设计与制造的他们在这篇文献中对嫦娥三号在登月过程中的各个阶段的各项参数都做出了详细而精确的推到与证明，通过这篇文献我们可以得到探测器近月点的位置信息为：经度：19.0464°W纬度：28.9989°N高度：15km俯视姿态角：85°资料中所给的近月点的速度V=1.6957km/s 5.1.1计算近月点直角坐标由2.4.2中构建的方程组，设点N为探测器在近月点处的位置，且该点所处在0度经线的位置上；已知近月点的位置为（19.0464°W，28.9989°N），高度为月球半径+15km，着陆点的位置为（19.51°W，44.12°N）可得到近月点的直角坐标：XYZ5.2远月点的参数5.2.1用近月点计算远月点位置如上图所示，由于近月点、月心、远月点在同一直线上，已知近月点的位置为经度：19.0464°W，纬度：28.9989°N。运用投影的方法可得到远月点的位置，计算方法如下：经度=180°-近月点经度=160.9536°W纬度=28.9989°S5.2.2远月点的直角坐标由2.4.2中的方程组结合远月点的经纬度数据可以得到远月点的直角坐标为X1Y15.3飞行器在近月点和远月点的速度5.3.1推导、计算月球的第一宇宙速率根据万有引力定律：GMmr2=mv(月球质量M=7.349*1022kg，G=6.67*10-11，月球的半径r=1.737013由此可以得到月球的第一宇宙速度为：V1=1.68km/s由于远月点处于嫦娥三号探测卫星的变轨点也就是变轨后的轨道与环月运动轨道的交点处，因此，远月点的速度即为探测器绕月运行的速度，即月球的第一宇宙速度，既得远月点的速度为V1=1.68km/s。5.3.2利用角动量守恒估算近月点速度首先由角动量的定义可知：L其中L为角动量，r为运动轨迹向量，p为动量，v为速度向量，m为探测器的重量。然后，刚根据角动量守恒定律：r×v=常量。已知，远日点（初次变轨时的位置）的路径长度：r1=100km+月球半径rr代入数据可得v2与资料中所给的v2=1.6957km/s相比较，计算所得的数据与资料所给数据产生的误差（3.8%）<5.3.3计算近月点速度方向由5.1.1中图所示，设点N为探测器在近月点处的位置，且该点所处在0度经线的位置上；点L为探测器着陆点的位置，由此构建方向向量OL，ON以及探测器的速度方向向量V，已知近月点的位置为（19.0464°W，28.9989°N），高度为月球半径r+15km，着陆点的位置为（19.51°W，44.12°N）根据向量的原理整理总结可得如下方程组：ON∙将ON，OL数据代入结合5.1.1中所计算出的近月点N的空进直角坐标，运用matlab进行运算后整理可得如下线性方程组：1719.5综上所述及计算我们可以得到近月点的速度方向向量为：V6、主减速阶段的轨道方程6.1动力学方程根据2.3中牛顿力学分析；经过2.4.1坐标方程变换后，飞行器的动力学方程如下：axaya(其中Q=dmd6.2控制方程整理6.1中的方程，将其打开引入简记符号后，系统状态方程可表示为：dX其中，a=b=c=cosαsinθsinφ-sinαsinθcosφggg以上式子中α，由上式运动方程可知，决定飞行器状态的变量为：sr=（X,Y,Z）6.3燃耗最优控制（1）性能指标α，β由以上结论可得出性能指标为：J=既降落过程中消耗的燃料最少的情况。（2）初始状态有上文数据可得，系统的初始状态为：v（3）控制变量Q（θ，φ（4）终止状态v注：动力下降过程中推进剂消耗约为1.4t。7.调整与控制阶段的控制策略7.1选择最优着陆地点在粗避障阶段，将2.4km时的高程图用Matlab进行剖分，将其剖分100m*100m的方格。根据高程图海拔数据，计算以方格为单位的每块区域的方差，方差最小的区域即为最平整的区域，即探测器选择的着陆区域。同理，在精避障阶段将高程图剖分为5m*5m的方格，根据海拔数据资料，计算以方格为单位的每块区域的方差，方差最小的区域即为最平整的区域，即探测器的着陆位置。数据结果如下：（1）2400米避障定位：图示区域：2300m*2300m

剖分块尺寸：100m*100m，共529个块

最优剖分块的高程方差：2.6882

最优剖分块中点的相对位置：(250m,1950m)（2）100米避障定位：图示区域：100m*100m剖分块尺寸：5m*5m，共400个块最优剖分块的高程方差：13.0012最优剖分块中点的相对位置：(52.5m,27.5m)参考文献[1]孙刚.测绘卫星的发展及技术现状.测绘科学与工程；Vo1．27，No．1Mar．2007[2]软着陆./view/63739.htm?fr=aladdin#3；Sep.14[3]杨峰，曹麒麟，段海星．基于DNSBlocklist的反垃圾邮件系统的设计与实现[J]．计算机工程与应用；2003(7)：11·12．45．[4]张洪华,关轶峰,黄翔宇,李骥,赵宇,于萍,张晓文,杨巍,梁俊,王大轶.嫦娥三号着陆器动力下降的制导导航与控制.中国科学：技术科学；Vol.44,377-384,2014[5]叶培建,彭兢.深空探测与我国深空探测展望.中国工程科学,2006,8:13–18[6]KlumppAR.ApolloLunarDescentGuidance,R-695.Cambridge:CharlesStarkDraperLaboratoryofMassachusettsInstituteofTechnology,1971.1–4[7]LeeAY,ElyT,SostaricR,etal.Preliminarydesignoftheguidance,navigation,andcontrolsystemofthealtairlunarlander.AIAA2010-7717[8]McHenryRL,BrandTJ,LongAD,etal.SpaceShuttleascentguidance,navigation,andcontrol.JournaloftheAstronauticalS,1979,27:1–38[9]KlumppAR.Apollolunardescentguidance.Automatica,1974,10:133–146[10]SavagePG.Strapdowninertialnavigationintegrationalgorithmdesignpart1:Attitudealgorithms.JGuidControlDynam,1998,21:19–28[11]王寨,李铁寿,王大轶.探月卫星变轨时的姿态控制研究.航天控制,2005,23:11–14[13]王大轶，李铁寿，马兴瑞．月球最优软着陆两点边值问题的数值解法[J]．航天控制，2000，3：44—49．[14]周净扬，周获.月球探测器软着陆精确建模及最优轨道设计.宇航学报.V01．28No．6.Nov.2007附录部分MatLab代码经纬坐标与直角坐标的转换---------------------------functionp=polar2orth(alpha,beta,radius)%alpha：经度，正负号（E:+,W:-）%beta：纬度，正负号（N:+,S:-）p(1)=radius*cos(beta)*cos(alpha);p(2)=radius*cos(beta)*sin(alpha);p(3)=radius*sin(beta);end（2）2400M处高程图最优落点位置算法------------------------------clear;clcM_img100_src=imread('A2.tif');M_img100=double(M_img100_src);N=50;l=length(M_img100)/N;flat_var=zeros(l);flat_mean=zeros(l);%求方格内的均值fori=1:lforj=1:lflat_mean(i,j)=mean(mean(M_img100((i-1)*N+1:i*N,(j-1)*N+1:j*N)));endendfori=1:lforj=1:lsum2=0;forp=(i-1)*N+1:i*Nforq=(j-1)*N+1:j*Nsum2=sum2+(M_img100(p,q)-flat_mean(i,j))^2;endendflat_var(i,j)=sum2/N^2;endendmin_var=min(min(flat_var));[x_m,y_m]=find(flat_var==min_var);x_m=x_m(1)y_m=y_m(1)min_varpos_w=x_m*N-N/2pos_h=y_m*N-N/2imshow(M_img100_src);holdon;plot(pos_w,pos_h,'^w','MarkerSize',10,'MarkerFaceColor','w')plot(pos_w,1:length(M_img100),'w')plot(1:length(M_img100),pos_h,'w')plot(length(M_img100)/2,length(M_img100)/2,'+b','MarkerSize',10)（3）100M处高程图最优落点算法-------------------------clear;clcM_img2400_src=imread('A1.tif');M_img2