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

承诺书

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

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

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

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

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

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A

我们的报名参赛队号为（8位数字组成的编号）：27028001

所属学校（请填写完整的全名）：咸阳师范学院

参赛队员(打印并签名) ：1. 栗楠

2. 尚旭

3. 张瑾

指导教师或指导教师组负责人(打印并签名)：张辉

（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。）

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

编号专用页

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

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

评

阅

人

评

分

备

注

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

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

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

摘要

对于嫦娥三号探测器的软着陆问题，我们需要考虑两个关键：着陆轨道和控制策略的设计。着陆轨道为从着陆准备椭圆轨道近月点至预定着陆点的6个阶段；控制策略要求满足各个阶段在关键时刻所处的状态并尽量减少着陆过程中的燃料消耗。

针对问题一，我们需要确定嫦娥三号探测器在其着陆准备轨道上近月点和远月点的位置以及对应的速度大小及方向，首先我们要建立动力学模型来刻画探测器在着陆准备轨道上的运动。然后运用模型适量计算和已知的预定着陆点来确定近月点、远月点在月球表面投影处对应的经纬度，根据已知条件和天体运动学对称规律，能量守恒定律，万有引力定律进一步求出近日点和远日点的具体位置、速度大小及方向；

针对问题二，首先，我们将嫦娥三号软着陆的6个阶段分为平面斜抛和垂直下降2个阶段分别进行着陆轨道的确定；其次，建立轨道动力学模型，并进行轨道优化，以燃料最优为目标函数对动力学模型进行处理，以获取在嫦娥三号6个阶段的最优控制策略。

1.问题重述

纵观人类探月历史，对于月球的探测大致可分为以下几个阶段：一是环月阶段，发射环月轨道器对月面进行较近距离的成像等科学探测；二是着陆阶段，发射无人月球着陆器和月球车，着陆月面并对月面进行巡视勘察；三是返回阶段，着陆月面后采集月球样本并返回地球；四是载人登月，由宇航员登陆月球采样并返回地球。进入21世纪，包括美国、俄罗斯、欧空局、日本、中国以及印度等又重新掀起了月球探测的高潮, 制定和正在实施的计划从环月、无人软着陆直至载人登月，乃至建立永久性月球基地。

我国相继成功发射“嫦娥一号”和“嫦娥二号”已经实现了绕月探测的目标, 嫦娥三号于2013年12月2日1时30分成功发射，12月6日抵达月球轨道。嫦娥三号在着陆准备轨道上的运行质量为2.4t，其安装在下部的主减速发动机能够产生1500N到7500N的可调节推力，其比冲（即单位质量的推进剂产生的推力）为2940m/s，可以满足调整速度的控制要求。在四周安装有姿态调整发动机，在给定主减速发动机的推力方向后，能够自动通过多个发动机的脉冲组合实现各种姿态的调整控制。嫦娥三号的预定着陆点为19.51W，44.12N，海拔为-2641m （见附件1）。

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

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

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

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

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

2.问题分析

对于嫦娥三号探测器的软着陆问题，我们需要考虑两个关键：着陆轨道和控制策略的设计。着陆轨道为从着陆准备椭圆轨道近月点至预定着陆点的6个阶段；控制策略要求满足各个阶段在关键时刻所处的状态并尽量减少着陆过程中的燃料消耗。

针对问题一，我们需要确定嫦娥三号探测器在其着陆准备轨道上近月点和远月点的位置以及对应的速度大小及方向，首先我们要建立动力学模型来刻画探测器在着陆准备轨道上的运动。然后运用模型适量计算和已知的预定着陆点来确定近月点、远月点在月球表面投影处对应的经纬度，根据已知条件和天体运动学对称规律，能量守恒定律，万有引力定律进一步求出近日点和远日点的具体位置、速度大小及方向；

针对问题二，首先，我们将嫦娥三号软着陆的6个阶段分为平面斜抛和垂直下降2个阶段分别进行着陆轨道的确定；其次，建立轨道动力学模型，并进行轨道优化，以燃料最优为目标函数对动力学模型进行处理，以获取在嫦娥三号6个阶段的最优控制策略。

3.符号说明

h

探测器在主减速阶段竖直方向的位移 s 探测器在主减速阶段水平方向的位移 m g

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

x v 探测器在进入主减速阶段的水平瞬时速度 0x v 探测器在进入主减速阶段的水平初速度 x a

探测器在进入主减速阶段的水平加速度

θ

着陆准备轨道近地点，预定着陆点分别与月心连线之间的夹角 G

万有引力常量 M 地球质量

m

嫦娥三号探测器的质量 k

月球的平均半径 H

近月点距月球表面的距离

λ

嫦娥三号进入主减速段时速度与水平方向的夹角

thrust F 发动机的推力，单位是牛顿； e v

以米/秒为单位的比冲；

m

单位时间燃料消耗的公斤数

4.模型假设

1.假设嫦娥三号月球探测器在变轨和软着陆的过程中只在一个固定的铅垂面内运动，不考虑侧向运动；

2.为了简化问题，假设月球无自传；

3.假设月球引力场为平行定常引力场；

4.假设作用于探测器上的所有外力经过其质心。

5.模型的建立与求解

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

5.1.1第I 阶段（主减速和快速调整阶段）平抛模型

嫦娥三号探测器绕月飞行后，经霍曼变轨后到达近月点，准备进行软着陆，根据天体力学知识，在探测器由着陆准备阶段进入主减速阶段的过程中在水平方向有很大的动能，假设在该阶段探测器发动机推力主要作用于水平方向，在水平速度被基本抵消之后进入最终着陆阶段，最后探测器以垂直姿态软着陆到月球表面。

图1：嫦娥三号软着陆过程示意图

所以在第I 阶段的竖直方向

???

???

?==g g t g m m 612

1h 2 水平方向有?????=-+=s a v v t

a v v x x x

x x x 22

020 在第I 阶段已知：???

?

?

??

??====g

g km h v s

km v m x x 616.120/7.10月球重力加速度竖直下降高度末速度初速度

假定在第I 阶段初始时刻00=t ，将已知数据代入（5.1.1）和（5.1.2）联立求解可得S=105.618km 。

在平抛模型的基础上，模型二加入了对嫦娥三号探测器在进入第I 阶段的速度方向的考虑对平抛模型进行改进。

5.1.2 第I 阶段斜抛模型

当飞行器由近月点进入主减速阶段时，由万有引力提供向心力和黄金代换公式

h

R mv h R GMm x +=+2

2

)(,2

h)g(R GM += 其中G 为月球万有引力常量：2211/1067.6kg m N G ??=-,M 为月球质量，m 为嫦娥三号质量，R 为月球平均半径为1737.013km ，H 为近月点距离月球表面的高度为15km ，x v 为嫦娥三号在近月点水平方向的速度分量。

图2：考虑速度方向时的软着陆过程示意图

则：s km x /691.1v =

km

km

7.1691.1cos =

λ ，则近月点嫦娥三号的速度方向与水平方向的夹角为 1.8。

假设嫦娥三号飞行器在第I 阶段（主减速段和快速调整段），所有推力都作用于B 点。即由A 到B 过程中不受推力，则有：

???

???

?==g 61g 2

1h 2m m t g 代入数值可得t=124.2118s 。

则嫦娥三号从A 到B 点（主减速段和快速调整段）过程中所经历的时间为124.2118s.水平方向的位移为：m t v x k 684.209S ==

由于在月球上空绕月球运行的探测器所受的力是月球对它的万有引力。因此探测器绕月球运行的圆心必与月球质心重合。而这样的轨道有多种，其中比较特殊的有与赤道共面的赤道轨道和通过两极点上空的极地轨道，当然也应存在着与赤道平面成某一角度的绕月轨道，只要圆心在地心，就可能是探测器绕月球运行的轨道。

假设嫦娥三号探测器绕月球南北两极进行极地绕月运动，因而探测器的绕月以及软着陆过程皆在xoz 平面内。则AOB ∠的余弦值可在直角三角形AOB 中算出：

OA

S

=

θarctan 由月球参数可得OA=（1737.013+2.4）km 又S=209.684km 可得o 84.6=θ

图3：嫦娥三号软着陆过程立体示意图

已知嫦娥三号的预定着陆点经纬度为19.51W ，44.12N 。近月点在月球表面的投影的经纬度为19.5W ，50.96N 。则近月点在月球19.5W ，50.96N 上方15km 处的位置，又因为绕月轨道在同一平面内则远月点在月球160.5E ，50.96S 上空100km 处。

5.1.3嫦娥三号探测器速度大小和方向的确定

在准着陆轨道上，嫦娥三号探测器从远月点向近月点运动过程中，遵循能量

守恒定律,则：)(mg 2

1)(m 212m 2

2121n l mv w n l mg v m ++=-++；

其中近月点速度1v =1.7km/s,月球重力加速度m g =2/m 68

.9s ，1l 为近月点到月球表面的距离，月球平均半径n=1737.013km,着陆点的海拔w =-2640km,2v 为远月点的速度，远月点到月球表面的距离2l =100km,将上述已知量代入方程中，得：2v =1.613km/s,

图4：着陆准备轨道切面示意图

如上图所示，近月点与远月点的速度方向，即为其相应速度在x 轴与y 轴方向上的投影。

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

5.2.1着陆轨道的确定

着陆轨道的确定可将6个阶段简化为2个阶段，其中第I 阶段为包括主减速阶段和快速调整阶段的斜抛运动。第II 阶段为包括粗避障阶段、精避障阶段，缓速下降阶段的

1.平面斜抛动力学模型

着陆轨道的主减速段的区间是距离月面15km 到3km 。该阶段的主要是减速，实现到距离月面3公里处嫦娥三号的速度降到57m/s 。

快速调整段的主要是调整探测器姿态，需要从距离月面3km 到 2.4km 处将水平速度减为0m/s ，即使主减速发动机的推力竖直向下。

此段可将月球视为平面来建立月球平面直角坐标系如图 4 所示

图 4 所示的月球平面直角坐标系, 原点o 为下降轨道上制动发动机点火点在月球表面的投影,00y x ，为下降轨道参考系纵向平面,着陆器的下降轨迹位于此平面内。该图表示的是符合重力转弯软着陆的情况，即反推力F 的方向与下降速度方向相反。对于这样的情况，沿两坐标轴方向有如下的动力学方程

???

???

?--=--='=''-=-='

=''m y m y x x g

m v Fv g m F v y m v

Fv m F v x θθsin cos （5.2.1）

上式中, m 为飞行器质量，在短时间内可视为常值;m g 为月球表面的重力加速度，始终垂直于月球表面且为常值；α为飞行路径角， 即为下降速度矢量v 与0x 轴的夹角，从0x 轴开始逆时针量起为正,v 为下降速度的模,在下降速度 v 和垂直于下降速度v 两个方向还可建立如下的动力学方程。

????

????

?=-=-=，e m v F

dt dm v g dt

d v dt dh

,)

sin (,cos βαβ （5.2.2） 根据第I 阶段平抛模型计算和已知有：

1.8==βα，

9.81-90==αβ，

2/m 68

.9g s m =

，2940m/s =v e

代入（5.2.2）式中，可得：

斜抛阶段动力学方程：

????

??

??

?=-=-=2940)9.81sin *63.1(,9.81cos F dt dm v dt

d v dt dh

α （5.2.3） 2.垂直下降动力学模型

在竖直下降阶段，着陆器距离月面很近，且着陆器几乎沿竖直方向下降。因此，该段仍可采用平面月球动力学模型，如图4所示，理想情况下，着陆器在着陆段沿竖直方向下降，则可在平面月球二维模型基础上简化为一维垂直动力学模型，即要求其中的飞行路径角α= 90°。因此,（5.2.1）式可简化为

m g m

u

F dt dW dt y d -?==202 （5.2.4）

式中, u 为制动推力 F 的开关控制量。着陆段一维垂直下降过程如图 5 所示.

图 5 着陆段下降过程示意

5.2.2月球软着陆阶段最优控制策略研究 1.系统模型的建立

取月心o 为坐标原点，oy 指向着陆转移轨道的近月点；r ∈ R+ 为探测器到月心的距 离；θ是oy 和or 的夹角；）（t ?为推力方向与or 垂线的夹角；F 为制动火箭的常值推力大 小，F 取max F 或0。嫦娥三号探测器的质心运动学方程为

???

????????????-=+-==+-==，，，，，e thrust v F dt

dm v m F r dt d dt

d r r u m F dt dv v dt

dr

)2cos (1sin 22ω?ωωθ

ω? 假定在初始时刻00=t ，终端时刻 f t 为任意值，软着陆的初始条件由探测器在椭圆轨道近月点处的状态为

0000)0()0(00)0()0(m m v r r =====，，，，ωωθ

其中月心到近月点的距离15)km 1737.013(0+=r ，探测器初始质量t m 4.20=。满足探测器在月面实现软着陆的终端条件为，，，0)()()(===f f f f f t v t v r t r ω对于

嫦娥三号探测器，最优控制策略是),(*)(t u t u =其中T F t u ],[)(?=，使得着陆过程燃耗最小，即使性能指标函数

?=f

t dt m

J 0

取最小值．

将软着陆过程中消耗燃料分为四大阶段，即为（1）：主减速和快速调整阶段；（2）：粗避障阶段；（3）：精避障阶段；（4）：缓速下降阶段。 其中0m 为初始质量，)(i t F 为各阶段所受推力，

与水平方向的夹角，为主减速阶段初始速度θ为各阶段的加速度。i a （1）在主减速段和快速调整阶段，对着陆器进行受力分析，得到：

000sin )0(a m F g m m =-θ

1

21a h t =

则在这阶段消耗燃料质量为e

v t m 1

187500-=

?（1） （2）粗避障阶段，推力方向向下，则 1187500)(t t F -=

210110)()()(a m m t F g m m m ?-=-?-，

1

222a h t =

则在这阶段消耗燃料质量为e

v t t m )

(87500212+-=

?（2）

（3）精避障阶段，推力方向向下，则)(87500)(212t t t F +-= 2

323a h t =

221021210)()()(a m m m t t F g m m m m ?-?-=+-?-?- 则在这阶段消耗燃料质量为e

v t t t m )

(875003213++-=

?（3）

（4）缓速下降阶段，推力方向向下，则 )(87500)(3213t t t t F ++-= 3

424a h t =

332103213210)()()(a m m m m t t t F g m m m m m ?-?-?-=++-?-?-?- 则在这阶段消耗燃料质量为e

v t t t t m )

(8750043214+++-=

?（4）

消耗的燃料总质量为4321m m m m m ?+?+?+?=总（5） 将式（1）（2）（3）（4）代入（5）中，可得：

2940

8162432300004

321t t t t m ----=总 （6）

由上式可知，消耗燃料的总质量总m 与各阶段所经历的时间存在线性关系，要进行优化控制，则改变各阶段的时间即可达到目的。

2.粗避障阶段

粗避障段嫦娥三号对正下方月面2300×2300m 的范围进行拍照，获得数字高程图，并分析数字高程图，避开大的陨石坑，实现在设计着陆点上方100m 处悬停，并初步确定落月地点，并使嫦娥三号在月面的垂直投影位于预定着陆区域的中心位置。

根据附件提供的*.tif 文件，

1）首先，我们使用Matlab 中的“imread ”命令对具有灰度数的图片文件进行读入

2）同时建立扫描区域，排除过高或过低的位置；对扫描区域内的灰度数进行比较得出方差，确定一片相对较为平坦的区域并得出区域中心位置。

3）确定探测器位置，并与选定区域中心位置进行比较得出需要位移的大小。 其Matlab 程序如附录—程序1所示。（注：该程序运行需较长一段时间） 在确定需要位移的大小和方向后，让嫦娥三号向确定方向进行位移避障。 3.精避障阶段

精避障阶段要求嫦娥三号悬停在距离月面100m 处，对着陆点附近区域100m 范围内拍摄图像，并获得三维数字高程图。分析三维数字高程图，避开较大的陨石坑，确定最佳着陆地点，并实现在着陆点上方30m 处水平方向速度为0m/s 。 同粗避障阶段相似，其Matlab 程序如附录—程序2所示，在确定需要位移的大小和方向后，令嫦娥三号向确定方向进行位移精避障。

问题三的求解 5.3误差分析和敏感性分析 误差分析

由于各种假设条件以及忽略多种外界环境，所以模型的建立过程中存在一定的误差，

敏感性分析

对问题二模型的建立式（6）进行化简，得：

43210027.0005.0008.0t 01.0-2.10m t t t ---=总

（1）总m 对1t 的敏感度为：

4

3214

3214321110027.0005.0008.001.02.100027.0005.0008.001.02.10-0027.0005.0008.0%)101.0(2.10m )%,1t )t (s t t t t t t t t t t t t m n m m ----++++---+-=

-+=总

总

总总（，

=0.1%

根据结果可得，当1t 增加1%时，总m 增加0.1%，s<5%,1t 和总m 的微小变化对模

型结果影响不大，即模型可用。

（n 为当未知量432t ,t ,t 取定值1时，未知量所对应的系数和为定值） （2）总m 对2t 的敏感度为： 同理：总

总

总总，m n m m m -)%,1t ()t (s 22+=

=0.1%

根据结果可得，当2t 增加1%时，总m 增加0.1%，s<5%,2t 和总m 的微小变化对模型结果影响不大，即模型可用。

（n 为当未知量431t ,t ,t 取定值1时，未知量所对应的系数和为定值） （3）总m 对3t 的敏感度为： 同理：%1.0m )n %1t )t (s 33=-+=

总

总

总总，（，m m m

根据结果可得，当3t 增加1%时，总m 增加0.1%，s<5%,3t 和总m 的微小变化对模型结果影响不大，即模型可用。

（n 为当未知量421t ,t ,t 取定值1时，未知量所对应的系数和为定值） （4）总m 对4t 的敏感度为： 同理：%1.0m )%,1t )t (s 44=-+=

总

总

总总（，m n m m

根据结果可得，当4t 增加1%时，总m 增加0.1%，s<5%,4t 和总m 的微小变化对模型结果影响不大，即模型可用。

（n 为当未知量321t ,t ,t 取定值1时，未知量所对应的系数和为定值） （5）综合分析：

根据上面四种结果可得，（总与m )t ,t ,t ,t 4321的微小变化对模型结果影响不大，则模型可用。

6.模型分析与评价

本论文是关于嫦娥三号软着陆轨道设计与控制策略的模型建立与全面综合分析，本题要求具有扎实的物理、天文、数学、天体动力学以及航空学基础，问题难度系数高，实践性较强，深度较广，所需的知识背景较为复杂，在全面建立模型的基础上，我们对模型进行以下评价：

模型优点：

（1）本模型是在查阅大量背景知识的基础上，并对各项数据进行全面分析与深刻理解，对各项物理和数学问题进行了模拟与刻画，可靠性高，实践性强。（2）模型进行了必要的、合理的假设，使问题的主要特征更加突出与易懂，并进行了适当的简化，使问题更容易着手。

（3）模型可靠性较强，有重大的实践意义，科学价值较高，比较符合实际的要求。

模型缺点：

（1）本模型对问题的解决还需待优化，方案对一些实际状况进行了忽略，简单处理了外界环境，对参数的刻画不是很明确。

（2）本模型在数值处理方面存在一定的误差，对系数精确度的描述不突出。（3）在处理问题三过程中，进行了文字方面的客观描述，缺少对数值的处理与分析，缺乏解决问题的应用性。

参考文献

[1] 姜启源谢金星叶俊，数学模型（第四版），北京，出版社，2012年12月。

[2]王鹏基张熇曲广吉，月球软着陆飞行动力学和制导控制建模与仿真，中国科学杂志社，第39卷第三期：521~527，2009年。

[3]张建辉张峰，月球软着陆轨道优化方法比较研究，工程数学学报，第29卷第三期:356~363,2012年6月。

[4]单永正段广仁吕世良月球软着陆的最优控制，光学精密工程，第17卷第9期:2154~2158,2009年9月。

附录

程序一：

clc

clrae

close all

s=imread('C:\Users\dell\Desktop\A\u3.tif');%录入附件3的图像的灰度数

x=ones(60,60); %建立单位扫描区域

[m,n]=size(s);

F=+inf;

z=[0 0]; %记录最小方差矩阵的首坐标

for i=1:m-59

for j=1:n-59

x=s(i:i+59,j:j+59);

if max(max(x))>150|min(min(x))<70 %排除过高或过低的位置break;

else

if max(max(x))-min(min(x))>=22 %最高点与最低点不得超过2.2米break;

else

e=var(x(:)); %通过命令函数计算方差

if e

F=e;

z(1)=i;

z(2)=j;

end

end

end

end

d=z+30 %得出所需矩阵的的中心位置

o=[m/2,n/2]; %确定飞行器的位置

if d(1)-o(1)>0&d(2)-o(2)>0

q=180/pi*atan((d(2)-o(2))/(d(1)-o(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2)));

fprintf('x y %f %f.\n',q,l);

end

if d(1)-o(1)>0&d(2)-o(2)<0

q=180/pi*atan((o(2)-d(2))/(d(1)-o(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2)));

fprintf('x -y %f %f.\n',q,l);

end

if d(1)-o(1)<0&d(2)-o(2)>0

q=180/pi*atan((d(2)-o(2))/(o(1)-d(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2)));

fprintf('-x y %f %f.\n',q,l);

end

if d(1)-o(1)<0&d(2)-o(2)<0

q=180/pi*atan((d(2)-o(2))/(d(1)-o(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2)));

fprintf('-x -y %f %f.\n',q,l);

end

%判断位移的大小和方向

程序2：

clear;

clc;

s=imread('C:\Users\dell\Desktop\A\u4.tif'); %录入附件4图像的灰度数x=ones(60,60); %建立单位扫描区域

[m,n]=size(s);

F=+inf; %方差记录%

z=[0 0]; %记录最小方差矩阵的首坐标

for i=1:m-59

for j=1:n-59

x=s(i:i+59,j:j+59);

if max(max(x))>100|min(min(x))<40 %排除过高或过低的位置break;

else

e=var([x]) %计算方差

if e

F=e;

z(1)=i;

z(2)=j;

end

end

end

end

d=z+30 %得出所需矩阵的的中心位置

o=[m/2,n/2]; %飞行器的位置

if d(1)-o(1)>0&d(2)-o(2)>0

q=180/pi*atan((d(2)-o(2))/(d(1)-o(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2))); fprintf('x y %f %f.\n',q,l);

end

if d(1)-o(1)>0&d(2)-o(2)<0

q=180/pi*atan((o(2)-d(2))/(d(1)-o(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2))); fprintf('x -y %f %f.\n',q,l);

end

if d(1)-o(1)<0&d(2)-o(2)>0

q=180/pi*atan((d(2)-o(2))/(o(1)-d(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2))); fprintf('-x y %f %f.\n',q,l);

end

if d(1)-o(1)<0&d(2)-o(2)<0

q=180/pi*atan((d(2)-o(2))/(d(1)-o(1)));

l=sqrt((d(1)-o(1))*(d(1)-o(1))+(d(2)-o(2))*(d(2)-o(2))); fprintf('-x -y %f %f.\n',q,l);

end

%判断位移的大小和方向