2014数学建模A题资料

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

摘 要：

嫦娥三号卫星采用的是软着陆方式登陆月球，在卫星高速飞行的情况下，我们要精确地在月球预定区域内实现软着陆，需要对其运行轨道进行设计并制定相应控制策略。由于天体的运动均满足开普勒三大定律以及总能量守恒定律，我们据此建立一系列的方程，最终求得卫星在近月点处的速度大小s km a GM c a b Va /6922.1/*)/(≈-= 求得卫星在远月点处的速度大小s km a GM c a b Vb /6139.1/*)/(≈+=。其速度方向均为当前运动轨道的切方向。

嫦娥三号的软着陆过程又分为6个阶段，对每个阶段都要进行一定的控制。通过对数值的分析，我们制定如下策略：在主减速阶段，首先卫星的发动机应全部用于水平方向的减速，竖直方向以自由落体状态加速下降，直至竖直方向速度达到56.8505m/s ，水平方向速度为4.9522m/s ，耗费燃料90kg ，卫星下落高度约为14695m 。之后发动机推力仍以最大输出工作，其推力在竖直方向上瞬间产生约3752N 的向上阻力，在水平方向上产生6494N 的阻力，使得卫星在主减速的后半阶段在竖直方向开始匀速运动，水平方向继续做减速运动水平速度最终降为0，直至卫星降落到距月3000m 高度处；在快速调整阶段，姿态发动机调整卫星的运动方向为竖直向下，主发动机继续工作，发动机在水平方向上产生1097N 推力，竖直方向上产生3752N 推力，使得卫星水平方向的速度降为0m/s ，竖直方向继续以57m/s 的速度做匀速运动；在粗避障阶段，卫星对其正下方月面进行拍照获得数字高程图，利用matlab 软件对图像进行灰度色差分析，进而初略确定了一个着陆点的像素点坐标(193，1169），过程中，姿态发动机根据图像进行卫星的位置调整，主发动机用于竖直方向上的减速工作，使得卫星在距月100m 处达到悬停状态；精避障阶段，进行更精细的月面拍摄，采用同粗避障段类似方法得到更精确的像素点坐标为(318，651),同时利用姿态发动机调整运动方向的同时主发动机工作对水平方向进行速度控制，使其在距月30m 处水平速度为0m/s ；缓慢下降阶段，发动对卫星竖直方向进行减速控制，使其在距月4m 处合速度变为0m/s ；在最后阶段，关闭发动机，使其做自由落体运动，最终成功着陆。

关键词：

软着陆 开普勒定理 优化控制 基于C 语言编程

一、问题重述

嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。此次嫦娥三号采用的是软着陆方法。

但嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是对着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道；着陆轨道为从近月点至着陆点，其软着陆过程共分为6个阶段，并要求满足每个阶段在关键点所处的状态，尽量减少软着陆过程的燃料消耗。

根据题目给出的已知量建立合适的模型求解分析下列几个问题：

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

（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。

（3）对设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。

二、问题分析

本问题的一个关键点在于求解出嫦娥三号卫星的绕月运行轨道。由于已知卫星的运行轨迹为近月点15km，远月点100km的椭圆形轨道。

对于第一问的关键就是求解出该椭圆轨道与月球赤平面的夹角。由于已知落点的经纬度坐标和软着陆过程中6个状态的数值，我们就可以使用逆推法求解出着陆准备轨道近月点和远月点的经纬度坐标位置。由于卫星的运动满足开普勒三大定律以及能量守恒定律，我们可以根据这几大定律确定几个等式联立方程组从而求得嫦娥三号在近月点和远月点的相应速度大小与方向。

对于第二问，要设计一种方案使得在软着陆的过程中耗能最少，并达到预设的各项指标。我们需要不断的去设计、计算、调整，在不断的尝试摸索中寻找出一个比较不错的软着陆方案。

三、模型假设及符号约定

3.1模型假设

（1） 假设卫星在主减速阶段的运动轨迹为抛物线

（2） 假设卫星在3000m 处基本就处于目标地点的竖直上方 （3） 假设月球为球体

（4） 假设在近月点即将开始软着陆状态的前后速度基本不变

3.2 符号约定 符号 含义 G 引力常量 M 月球的质量

)(t m t 时刻卫星的质量 r

月球的半径

'g 月球表面的重力加速度

Vi

i 位置处卫星的速度大小 ai i 方向的加速度

Si i 位置处相对应的面积大小 Ei i 位置处卫星的能量 Ti 到i 位置所用时间 Ft t 时刻发动机的推力 Ve

以m/s 为单位的比冲 ?

m

单位时间燃料消耗的质量 a

椭圆的半长轴 b 椭圆的半短轴 c

椭圆的焦点

四、 模型建立与求解

4.1 模型建立

4.1.1 问题1卫星轨道模型的建立

我们将嫦娥三号卫星绕月运行的椭圆轨道抽象出如下图所示的一个简单的几何图形，月球的月心位于椭圆的一个焦点F 上，椭圆的半长轴为a ，半短轴为b ，半焦距为c 。A 点为近月点，速度为Va ，B 点为远月点，速度为Vb 。易

知A 、B 两点距月心的距离c a La -=，c a Lb +=。在一个极小的时间段t ?内，卫星与月心连线扫过的面积分别为La t Va Sa **2/?=，Lb t Vb Sb **2/?=。由开普勒第二定律可知，卫星与中心天体连线在单位时间内扫过的面积相等，所以Sb Sa =，代入化简后可得公式(1)Va c a c a Vb *)/()(+-=；由于卫星运动的总机械能等于其动能和引力势能之和，所以在A 点，卫星的总机械

能)/(*2/12La GMm mVa Ea -+= （公式2），同理B 点的总机械能)/(*2/12Lb GMm mVb Eb -+= （公式3）

。卫星在运行过程中只有动能和引力势能之间的转化，机械能守恒，所以Eb Ea = （公式4）。[1]

由上面(1)(2)(3)(4)四个公式我们可以解得a GM c a b Va /*)/(-=，

a GM c a

b Vb /*)/(+=。

4.1.2 问题2卫星软着陆各阶段控制策略模型的建立

主减速段控制策略：

第一阶段，主减速发动机全力用于水平方向减速，竖直方向呈自由落体加速下落，至下落速度达到y V 1 (m/s)为止，用时记为1T ，质量减少到1M ；

1'*1T g y V = (1) 当10T t <<时，由Ve m F *?

=可得出下式，

t dt

m t m t

5510.224002400)(0

-=-

=??

(2)

三号卫星运行速度V 水平方向分量为：

(3)

在10T t <<内，卫星下落的高度2^'**5.0t g ht =, 水平方向位移

?-+-=

t

dt t t s 0

))2210.22400ln(2940)2400ln(29401700

()(

(4)

第二阶段，竖直方向以y V 1匀速下降，主减速发动机仍以最大推力工作，)(t F 的竖直方向分力为'*)(g t m -，水平方向分力为22)'*)((7500)(g t m t Fx -=

。

记第二阶段结束时刻，即距月面3000m 时为2T ，当21T t T <<时，

t t t m 5510.224002940/*75002400)(-=-= (5)

下落高度

)1(*11'**5.0)(2T t y V T g t h -+= (6)

水平方向加速度

)

5510.22400/()'*)*5510.22400(7500()(/)()(2

22t g t t m t fx t ax ---== (7)

水平方向速度

??

----=-

=t

T t

T ds s g s T Vx ds s ax T Vx t Vx 1

1

2

22)5510.22400/()')5510.22400(7500()1()()1()(

=?-----

)

5510.22400(')

15510.22400('22)/7500(5510

.21

)1(t g T g

du u u T Vx

=)15510.22400(')5510.22400('2

222)75007500ln

75007500(5510

.21)1(T g t g u

u u T Vx ----

+--(8)

水平方向位移 dt t Vx T s t s t

T

?+=1

)()1()(

(9)

快速调整段控制策略：

在快速调整阶段，卫星仍然以s m T g Vy /571'*≈=的速度下落，水平方向推力逐渐减小到0，这样卫星的方向就正对目标下方了。

在2T t >时，发动机给与卫星水平方向一个力Fx ，使得其水平方向的速度最终降为0m/s ；同时主发动机在竖直方向上也对卫星施加一个力，使其在竖直方向上做匀速运动。水平方向位移2^3**5.03)3(0T a T V T Sx -=。

水平方向推力为)(*)()(t ax t m t Fx =

，合力为

?

+=

t

T e

m dr

r m v g 2

)(12

两边对t 求导得：

)

(1)(2

't m v g t m e

m +=

再由2T t =时，)()(2T m t m =得，当96.422+<

)

(1exp(

)()(22

2T t v g T m t m e

m -+=

水平方向速度

2

22)()()()(2

T t T v dr r a T v t v x t

T x x x +-=-=?

水平方向位移

)

(5.0))()(()()()()(2

22222222

T t T t T T v T s dr r v T s t s x t

T x ---++=+=?

()96.42(+T s ，15000)就是近月点坐标。

粗避障段控制策略：

以下步骤的实施均在matlab 上实践，源代码见附录2

将抓拍到的图片利用matlab 读入，求出图像平均灰度，以此来表示理想降落地点的高度。但通过比较发现平均灰度并不能代表理想地点高度，于是求出灰度众数average ，作为理想降落高度。

然后新建图像，使新图像上的每点的灰度等于原图像对应点的灰度与灰度众数average 的差的绝对值。于是新图像上灰度越高的点表示越不理想的点。（matlab 运行效果图见附录3）

选取10 作为分界线，灰度值大于10的改为150，表示不理想的降落地点，灰度值小于10的改为0，表示理想的降落地点。将图片边界四周都像素点全部改为灰度150，表示不理想降落地点，因为越靠近图片边缘则越缺少信息。(matlab 运行效果图见附录4)

接下去通过将不理想的点向四周扩散来逐渐覆盖全图，留下0.001的空隙时停止扩散，则此时剩余的空隙则为相对理想的降落地点，求出这些降落地点中离中心飞船最近的地点，作为最终降落地点。

最后我们分别绘制了覆盖率为0.7时的图像、覆盖率为0.99的图像、覆盖率为0.999的图像(分别件附录5、6、7)。

精避障段控制策略：

此阶段只是在粗壁障阶段的基础上对着落点做了一个更精确的定位，所以选点的原理同上。

4.2 模型求解 4.2.1 问题1

由题目已知，卫星近月点为15km ，远月点为100km ，据此我们可以得出下面一个等式100)(15=-++r c c ，求解得2/9r c += km 。据此又可以导出椭圆的半长轴的大小r r c a *2/32415+=++=，半短轴22c a b -=

。然后将c b a 、、等数值

代入模型中，便可以直接求出卫星在近月点的速度大小s km V /6922.1a ≈，卫星在远月点的速度大小s km V /6139.1b ≈。

4.2.2 问题2 主减速阶段：

由天体运动黄金代换式'2g R GM =，得出22/6243.1/'s m R GM g ≈= 利用C 语言编程使用迭代法搜索T1的值： 具体编程算法见附录1 Step 1 令t=T1,

下落高度 2^'**5.0)(t g t h =, deltt=1 竖直方向速度t g t Vy '*)(=,

水平方向速度 （初始速度有第一问可知为1692m/s ）

)2210.22400ln(2940)2400ln(29401692)(t t Vx -+-= 水平方向位移

t

u ulunu t t S 2210.224002400

)(2210

.22940))2400ln(29401692()(---+-= Step 2 令t=t+ deltt

下落高度 )1(*2^1'**5.0)(T t Vy T g t h -+= 竖直方向速度1'*)(T g t Vy = 水平方向速度由（8）式计算可得 水平方向位移

Step 3 如果264112000)(+, then go to step

2

如果264112000)(+

step 1

如果264112000)(+>t h , then end

最后利用C 语言编程解得

s T 351=,m h 85.14695=,s m Vy s m Vx /8505.56,/9522.4==。

从而我们可以得出，在主减速阶段，首先卫星的发动机的全部动力要用于卫星水平方向上的减速，竖直方向以自由落体状态加速下降，直至竖直方向速度达到56.8505m/s ，水平方向速度为4.9522m/s ，此时卫星的质量约为2310kg ，耗费燃料90kg ，卫星下落高度约为14695m 。之后发动机推力仍以最大输出工作，其推力在竖直方向上瞬间产生约3752N 的向上阻力，在水平方向上产生6494N 的阻力，使得卫星在主减速的后半阶段在

竖直方向开始匀速运动，水平方向继续做减速运动水平速度最终降为0，直至卫星降落到距月3000m 高度处。

快速调整阶段：

根据分运动的等时性，可以利用卫星竖直方向上的运动方程t V h 0=，求得卫星在快速调整阶段的最大用时s V h t 52.10/0≈=，从而根据水平方向的运动公式

t a V V o *+=求得水平方向的最小加速度s m a /475.0=。所以可以得出该过程中发

动机在水平方向上的最小推动力N Fx 1097min =，竖直方向上的推力

N mg Fxy 3752'==，并求得水平位移m S 15234≈。

五、 模型评价与推广

5.1 模型评价

优点：本模型采用了化繁为简的方法，将复杂问题简单化，把卫星绕月运行抽象出一个椭圆的几何问题再结合物理模型利用数学方法求解相应的运动数值

缺点：本模型很多地方理想化，比如忽视了月球的扁率，卫星运行中姿态是变化的某些部分不应该用质点来处理

5.2 模型推广 利用本模型的算法可以求出任何天体在每个时刻的运行速度及其他物理量，方便研究天体运动，方便针对性的对需要探测的星球设计相应的卫星轨道。

六、 参考文献

[1] 王健伟 李兴，近日点和远日点速度的两种典型求法，物理教师，第34卷第6期：58，2013

七、 附录

附录1：

下面采用的是C 语言程序来摸索求解1T 的值（运行环境VC++ 6.0） 文件夹内find_t1.cpp 文件

#include #include

#include

#define g 1.6243

int main()

{

int t,t1,deltt=1;

t=t1=4;

double h,vx,vy,compare,u1,u2;

h=0.5*g*t*t;

vx=1692-2940*log(2400)+2940*log(2400-2.2210*t);

vy=g*t;

while(1)

{

compare=sqrt(57*57-vy*vy);

if(h<12000+2641 && vx>compare)

{

t=t+deltt;

u1=g*(2400-2.5510*t1);

u2=g*(2400-2.5510*t);

vy=g*t1;

h=0.5*g*t1*t1+(t-t1)*vy;

vx=1692-2940*log(2400)+2940*log(2400-2.2210*t)-(1.0/2.5510)*(sqrt(7500*7500-u1*u1)+7 500*log((7500-sqrt(7500*7500-u1*u1))/u1)-(sqrt(7500*7500-u2*u2)+7500*log((7500-sqrt(7 500*7500-u2*u2))/u2)));

}

else if(h<12000+2641 && vx

{

t1++;

t=t1;

h=0.5*g*t*t;

vx=1692-2940*log(2400)+2940*log(2400-2.2210*t);

vy=g*t;

}

else if(h>12000+2641)

{

printf("t1=%d\n",t1);

printf("h=%lf\n",h);

printf("vx=%lf\n",vx);

printf("vy=%lf\n",vy);

break;

}

}

return 0;

}

附录2：

big=imread('附件3 距2400m处的数字高程图');

[m,n]=size(big);

notenum=[1:256];

for i=1:m;

for j=1:n;

notenum(big(i,j)+1)=notenum(big(i,j)+1)+1;

end

end

max=0;

maxnum=0;

for i=1:256;

if max

max=notenum(i);

maxnum=i;

end

end

average=maxnum;

big1=big;

for i=1:m

for j=1:n

big1(i,j)=abs(int8(big1(i,j))-average);

end

end

big2=big1;

for i=1:m

for j=1:n

if big2(i,j)>15

big2(i,j)=150;

else big2(i,j)=0;

end

end

end

for i=1:m;

big2(i,1)=150;

end

for i=1:m;

big2(i,m)=150;

end

for j=2:n-1;

big2(1,j)=150;

big2(m,j)=150;

end

for k=1:200

num=[m:n];

numsum=0;

for i=1:m;

for j=1:n;

if big2(i,j)>20;

num(i,j)=1;

numsum=numsum+1;

else num(i,j)=0;

end

end

end

if numsum/(m*n)>0.999

break

end

for i=2:m-1;

for j=2:n-1;

sum=0;

if num(i,j)==1;

continue

else

for i1=-1:1;

for i2=-1:1;

if num(i+i1,j+i2)==1;

sum=sum+1;

end

end

end

if rand(1)*3

big2(i,j)=150;

end

end

end

end

end

distance=[m,n];

for i=1:m;

for j=1:n;

if num(i,j)==0;

distance(i,j)=(i-m/2)*(i-m/2)+(j-n/2)*(j-n/2);

else

distance(i,j)=m*m+n*n;

end

end

end

min=m*m+n*n;

mini=0;

minj=0;

for i=1:m;

for j=1:n;

if num(i,j) == 0;

if min > double(distance(i,j));

min=double(distance(i,j))

mini=i;

minj=j;

end

end

end

end

附录3

附录4

附录5 覆盖率为0.7的图像

附录6 覆盖率为0.99的图像

附录7 覆盖率为0.999的图像