论文初稿晚.稿.doc

个人收集整理 勿做商业用途
参赛密码
（由组委会填写）

全
第九届“华为杯”
全国研究生数学建模竞赛

学 校
华北电力大学（保定）
参赛队号
10079017
队员姓名
1. 高成彬
2. 李少岩
3. 刘贺晨

参赛密码
（由组委会填写)

第九届“华为杯"全国研究生数学建模竞赛

题 目 基于遗传算法的无源探测
空间飞行器主动段轨道估计与误差分析
摘 要：

关键词：遗传算法 、参数估计、PID、误差分析、轨道预测
一、问题重述
有些国家会发射特殊目的的空间飞行器,如弹道式导弹、侦察卫星等.对他国发射具有敌意的空间飞行器实施监控并作出快速反应，对于维护国家安全具有重要的战略意义.卫星居高临下，是当今探测空间飞行器发射与轨道参数的重要平台。装置于卫星上无源探测器被动地接收目标辐射,其目的是为了推断空间飞行器的轨道参数
在无源探测中,多使用红外光学探测器，单个红外光学探测器不具备测距能力,但借助多颗（含两颗）观测卫星的同步观测能够进行逐点定位，再结合空间飞行器的运动模型，可以进行轨道参数估计。在单星观测条件下,利用空间飞行器轨道的特殊性，结合较强的模型约束也可得到一定精度轨道参数估计。由于受大气影响，垂直上升段的火箭尾焰不易观测,程序拐弯段的运动方程又较为复杂，所以本题重点关注重力斜飞段的后程段，本题所附仿真数据也集中于此段.
本题以中低轨近圆轨道卫星为观测星座对假想的空间飞行器进行仿真观测，生成仿真观测数据，要求利用仿真观测数据,对假想空间飞行器的轨道参数进行估计。
我们将研究下列问题:
1. 观测卫星在任意时刻的位置计算是估计的前提，根据satinfo.txt和观测卫星的简化运动方程(2），计算09号观测卫星在50.0s、100.0s、150.0s、200.0s、250.0s五个时刻的三维位置。结果保留6位有效数字.
2. 给定的仿真数据下，06号和09号观测卫星对0号空间飞行器形成了立体交叠观测，我们将结合立体几何知识按照逐点交汇定位的思路，给出0号空间飞行器在公式（1）框架下的轨道估计，并选取适当的和的表示模型。按照从50。0s到170。0s间隔10.0s进行采样，计算并列表给出0号空间飞行器在各个采样点的位置和速度，并给出估计残差。结果保留6位有效数字.同时绘制0号空间飞行器的三个位置t—x、t-y、t-z和三个速度t-vx、t—vy、t-vz曲线示意图。
在同时考虑系统误差的条件下，进一步研究下列问题：
3. 若06和09号两颗观测卫星均有可能带有一定的系统误差,对系统误差进行正确的估计能够有效提高精度.我们将考虑利用上述的逐点交汇方法能否同时对系统误差进行估计?若不能，我们将寻找其他的思路能够同时估计系统误差与轨道.并给出与第二问要求相同的结果，以及两颗观测卫星的系统误差估计结果，共六个数值,分别是两颗卫星的。
4. 我们对只有09号观测卫星单星观测的01号空间飞行器进行轨道估计并进一步考虑在同时有多颗观测卫星观测多个空间飞行器的情况下能否联合进行系统误差估计？
5.
二、模型假设
1. 本题假设为形式简单的严格单调递减的非负函数.的方向与飞行器的速度相反方向接近或相同，其大小较为稳定。
2. 不考虑飞行器上各点的速度差异；
3. 不考虑飞行器质心因燃料消耗而产生的位移；
4. 作用于飞行器上的所有外力都通过质心，也就是说不计导弹在飞行过程中绕其质心的转动。
5. 本题假设观测坐标系与基础坐标系之间的转化误差很小，可以忽略不计；
6. 我们仅考虑单级火箭模型，对于多级火箭的模型不予考虑，以简化模型.

三、符号说明

：基础坐标系
：观测坐标系
：空间飞行器在基础坐标系下的位置矢量
:Vr(t)是燃料相对于火箭尾部喷口的喷射速度
：飞行器瞬时质量
：地球引力常数(本题中地球引力常数取）
：飞行器所受的外力加速度之和
：火箭产生的推力加速度
：观测坐标系下的两个无量纲比值
由于本文符号较多，其他符号在论文中第一次出现时,再进行说明.

四、问题分析：
1、问题1：
文件中已经给定基础坐标系下的位置和速度初值，可以利用常微分方程数值解方法计算空间飞行器的运动轨迹。09号观测卫星的运动方程满足题中形式较为简单的微分方程。此外，把09号观测卫星在基础坐标系中的位置和速度作为微分方程的初始条件，通过编程即可得到任意时刻的09号卫星的位置的数值.
2、问题2:
已知06号和09号观测卫星对0号空间飞行器的进行了立体交叠观测,形如meadata_i_j.txt的数据文件中存放着编号为i的卫星对编号为j的飞行器的仿真观测数据信息，按照时间顺序分行，每行分三列，分别是观测时刻以及对应观测数据。在某一观测坐标系中，已知某一时刻的观测数据，显然在空间上确定了一条直线。而两个观测点如果同时观测一个空间飞行器，则是形成了两条直线。如果系统不存在任何误差，当将坐标归算同一坐标系下（不妨为基础坐标系），这两条空间直线将会是共面相交于一点的，这点便是被观测的某一空间飞行器。但由于09号及06号飞行器任意时刻位置是通过插值得到的，以及如题干所说,系统测量值是有误差的，其中影响最大的是三轴指向误差，所以两条观测直线为异面直线。论文【】中提出了一种基于加权最小距离平方和的方法来确定被观测点.本文采用这个思路,对于每一时刻的两条“测量直线”，求其公垂线段中点坐标，便是被观测飞行器的空间位置。选取适当的和的表示模型，设置未知参数的取值范围，引入遗传算法来确定未知参数的值，便可以给出预测出的卫星轨迹.
3、问题3：
为估计系统误差,重新审视空间逐点交汇的思路.如果假设通过逐点交汇法处理的06号和09号观测卫星对于00号空间飞行器的观测数据而得到的00号空间飞行器的位置是准确的.则用逐点交汇得到的数据，通过坐标转换去反算06号09号观测卫星的观测值，便可以得到两个系统分别的系统误差估计结果。
4、问题4：
通过第三问得出的09号空间飞行器的系统误差,对09好空间飞行器对01号空间飞行器的原始观测数据进行处理，得到修正值.根据文献【单行预测的新方法】，利用单一观测卫星对飞行器实施观测和预警，是不满足可观测条件的，因此必须运用
五、问题的求解：
5。1问题1的求解
5.1.1空间飞行器主动段的简化运动方程
根据变质量质点的动力学，空间飞行器在基础坐标系下的主动段的简化运动方程如下：
(1)
其中向量表示飞行器所受的外力加速度之和,表示火箭产生的推力加速度,为瞬时质量；是质量变化率；为空间飞行器在基础坐标系下的位置矢量；表示对时间的二阶导数，即加速度；为地球引力常数（本题中地球引力常数取），为了更明确地表示推力加速度的方向，取的是燃料相对于火箭尾部喷口的喷射速度的逆矢量。
当研究对象为观测卫星时,方程(1）中如果只保留右侧第一项，则可以表示观测卫星的简化运动方程为:
（2)
5.1。2 飞行器常微分方程组的建立
在给定基础坐标系下的位置和速度初值情况下，可以利用常微分方程组数值解方法计算空间飞行器的运动轨迹。
因为位置矢量和速度矢量有如下关系：
（3）
而根据题意有：
（4)
则问题转化为微分方程组（5)：
（5）
初始位置和初始速度如下:

（6）
由于原始数据的时间间隔为0。2s，采用间隔太小，为了得到更高精度的数据，这里对原始数据进行插值，所用方法为三次样条插值,采用时间间隔为0.1s。得到了问题中09号观测卫星在50.0s、100.0s、150。0s、200。0s、250。0s五个时刻的三维位置。
列表如下:
表1 整点时刻09号观测卫星的三维位置
时刻/s
50。0
100.0
150。0
200.0
250。0
Xc/m
9.09204E+06
9。08416E+06
9。07616E+06
9.06806E+06
9.05985E+06
Yc/m
1。73211E+06
1。75437E+06
1.77661E+06
1。79883E+06
1。82102E+06
Zc/m
1。73211E+06
1。75437E+06
1。77661E+06
1.79883E+06
1.82102E+06

直观起见，画示意图如下：

图1整点时刻09号观测卫星的三维位置示意图一
上图中右上角的分布较为集中的五个红点即是整点时刻09号观测卫星的三维位置。

5。2 问题二的求解
5.2.1 坐标系转换
已知09号和06号卫星对00号飞行器的观测角,但是其是在各自的测量坐标系下测得的数据.为了便于分析，我们将两个测量坐标系转换到同一个坐标系下，即基础坐标系。我们采用的方法为通过坐标轴旋转，之后再对测量数据进行相关变换。

图2 基础坐标系与观测坐标系示意图

具体方法如下：
首先将观测坐标系的轴反向，然后绕轴顺时针旋转角，最后绕新的轴顺时针旋转角,则将观测坐标系转换成为了与基础坐标系平行的坐标系系统。转换之后两坐标系向量可进行向量之间的加减。
其中:
(6）
X，Y,Z轴的旋转矩阵及其将X轴方向的转换矩阵如下：

(7)
（8)
（9)
(10)
式中α为旋转角度.转动的角度正负规定为：从原点出发向旋转轴正向观看，顺时针旋转为正，逆时针旋转为负。
按照如上表述,将观测坐标系转换为基础坐标系的表达形式为：
（11）
式中的，，分别在观测坐标系中表示UEN坐标系中的向上(UP）、东（EAST）和北（NORTH)三个方向.
将式（）（）()带入式（)中可得如下关系：
（12)
5.2.2 逐点交汇定位确定飞行轨迹
5.2。2.1 基础坐标系下09号06号卫星观测确定直线的求取
在观测坐标系下，通过09号卫星对00号飞行器的观测角度可确定一条直线：
（13）
同理06号卫星对00号飞行器的观测角度也可确定一条直线,两条直线的交点即为00号飞行器的位置。
为了便于计算,我们将式中的，，通过上述的坐标转换用基础坐标系下的表示。不过值的注意的是转换之后的坐标与基础坐标系下的坐标还存在一个向量差，即：
（14)
将式中的设为K,则,分别表示为，,式中为卫星09在基础坐标系中的坐标。
将其带入上式中化简得到09号卫星确定的直线L1：
(15)
其中：
（16)
同理可得06号卫星确定的直线L2。
5。2.2。2 公垂线的求取
由于材料提供的两卫星的时间不同步，因此不能直接采用材料数据进行计算。为了得到相同时间下的卫星观测数据，我们将材料数据进行三次样条插值拟合，然后求取相同时间两卫星的观测数据.
由于卫星的观测数据存在系统误差.因此直接采用上述方法确定的两条直线可能没用交点，即无法确定00号飞行器的位置。为了解决这个问题，我们采用求取两直线公垂线中点代替00号飞行器位置坐标方法。
设直线L1，L2与其公垂线的交点M1，M2坐标分别为：
（17)
(18）
因为M1M2⊥L1，M1M2⊥L2,从而可求得，两直线与公垂线交点坐标可得，从而公垂线的中点由下式也相应确定。
(19）
5.2.2。4 数据处理的相关问题
原理上，按照上述方法便可得到任一对应测量时刻的0号空间飞行器的三维空间位置，但是由于09号和06号卫星的采样时刻不一致,所以在计算0号空间飞行器的三维空间位置对09号和06号观测卫星的观测值进行三次样条插值处理，得到间隔为0.1秒的数据。使得二者的数据在时刻上保持了一直，进而进行上述步骤以确定三维坐标的观测值。
对得到的0号空间飞行器的三维空间位置的观测值按照从50.0s到170。0s间隔10.0s 进行采样，得到以下数据:
表2 00号空间飞行器的三维空间位置及速度的观测值
时刻/s
50。0
60.0
70.0
80.0
90.0
x/m
—6。84379E+05
—6。85340E+05
—6.88875E+05
-6。94954E+05
-7.03625E+05
y/m
6。63547E+06
6。65952E+06
6.68241E+06
6。70417E+06
6.72469E+06
z/m
2.03978E+06
2。10125E+06
2.16198E+06
2.22192E+06
2。28098E+06
Vx（m/s)
6.35898E+02
-2。58869E+02
-5.42849E+02
—7.28937E+02
-8。40930E+02
Vy(m/s）
3。08770E+03
2.32655E+03
2.30462E+03
2。13288E+03
1。99813E+03
Vz(m/s）
6.65454E+03
6.07775E+03
5。97273E+03
5.94707E+03
5.87309E+03
时刻/s
100。0
110。0
120。0
130。0
140。0
x/m
—7.14927E+05
—7。28978E+05
—7.45842E+05
—7.65642E+05
—7。88506E+05
y/m
6。74394E+06
6.76186E+06
6。77841E+06
6。79349E+06
6.80715E+06
z/m
2。33904E+06
2。39603E+06
2。45179E+06
2。50617E+06
2.55902E+06
Vx（m/s)
—1.30419E+03
-1.53816E+03
-1。81152E+03
—2.14072E+03
—2。33987E+03
Vy(m/s)
1.90028E+03
1.65855E+03
1.52400E+03
1。49723E+03
1.27836E+03
Vz（m/s)
5。80251E+03
5。64522E+03
5。46408E+03
5.35154E+03
5。25209E+03
时刻/s
150.0
160。0
170.0
—
—
x/m
—8.14598E+05
—8。44093E+05
—8。77292E+05
-
-
y/m
6。81929E+06
6.82991E+06
6。83897E+06
-
-
z/m
2。61018E+06
2.65945E+06
2.70662E+06
—
—
Vx（m/s)
—2.73418E+03
—3。15503E+03
—3.42563E+03
—
-
Vy（m/s)
1.19272E+03
1。04478E+03
9。73568E+02
—
—
Vz（m/s)
5.03319E+03
4。85744E+03
4.56032E+03
-
-

下图较为清晰地表示出06号和09号卫星的轨迹和00号空间飞行器观测值确定轨迹。与题干中对于飞行器轨道阶段中主动段的重力斜飞段特点是相符的。

图3 观测卫星轨迹和飞行器50～170s轨迹示意图

5.2.3 对飞行器运动方程模型的求取
一般而言应为严格单调递减的非负函数.的方向一般应与飞行器的速度方向接近或相同，其大小一般较为稳定。
所以假定的模型如下：
（20）
其中是一三阶列向量，其任意元素值接近于零的随机数，k为一常系数。
之后可以根据情况来说明模型改进的问题，这是遗传算法的优势之一。
文献【？？】提出两种常用的变质量模型：
a、质量按线性规律变化
（21）
则：
（22)
（23）
则公式(1)可化为
(24)
b、质量按指数规律变化：
（25）
则：
(26)
(27）
则公式（1）可化为
（28）
在此我们首先选取（20）作为模型，选取(25）作为模型，其中是为调节飞行器的速度方向所加的一个随机小量，以模拟与飞行器速度方向的微小差异。将此代入式(28)中得：
（29）
这里和模型中都含有参数，因此本文将该问题转化为优化问题，目的是得到最优的参数.对于这种双参数的优化问题，如果采用传统遍历形式的暴力搜索方法速度将会非常慢。遗传算法具有较好的全局搜索性能,不容易陷入局部最优，能够同时对搜索空间中的多个解进行评估。因此，本文采用遗传算法进行参数寻优。
最优参数的定义为:使得由公式（*）确定的加速度与通过观测数据得到的加速度的平方和最小。基于此,定义遗传算法的目标函数为:
(30）
其中， 为通过已知数据通过对时间进行微分得到的,见公式（＊＊*)。为通过公式（**）得到的函数值，由于、未知，所以这里预先给定、在[-10,10］之间，利用遗传算法寻找最优的参数使目标函数的数值最小.遗传算法的参数中，最大遗传代数为150，变量的二进制位数为25。最终通过迭代，得到目标函数值最小时对应的参数=-0.0575，=—0。3592。
迭代过程如图（＊）所示

从图中可以看出，遗传150代以后,种群均值基本不变，达到最优的参数

表3 00号空间飞行器的三维空间位置及速度的拟合值
时刻/s
50.0
60。0
70.0
80.0
90。0
x/m
-6。84379E+05
-6.96214E+05
-7.05812E+05
—7。13592E+05
—7。19892E+05
y/m
6。63547E+06
6.69631E+06
6。74041E+06
6.77103E+06
6。79079E+06
z/m
2。03978E+06
2。19647E+06
2.32303E+06
2。42507E+06
2。50715E+06
Vx（m/s）
—3。45140E+03
—2.79987E+03
-2。27010E+03
-1.83931E+03
—1。48900E+03
Vy（m/s）
9.27677E+02
6.82436E+02
4。85148E+02
3。26329E+02
1。98387E+02
Vz(m/s)
4.57684E+03
3。69980E+03
2.98584E+03
2。40464E+03
1.93151E+03
时刻/s
100。0
110.0
120.0
130。0
140。0
x/m
—7.24991E+05
—7.29111E+05
-7.32436E+05
-7。35114E+05
-7.37265E+05
y/m
6.80181E+06
6.80577E+06
6。80404E+06
6。79770E+06
6.78763E+06
z/m
2。57298E+06
2.62558E+06
2.66740E+06
2.70046E+06
2。72637E+06
Vx（m/s)
-1。20411E+03
-9。72418E+02
—7.83976E+02
-6。30698E+02
—5。06008E+02
Vy（m/s)
9。52357E+01
1。19926E+01
-5.52599E+01
-1.09668E+02
—1。53761E+02
Vz（m/s）
1.54634E+03
1.23276E+03
9。77419E+02
7.69474E+02
6。00089E+02
时刻/s
150.0
160.0
170.0
—
—
x/m
-7.38989E+05
—7。40364E+05
-7.41455E+05
-
—
y/m
6。77454E+06
6。75898E+06
6.74142E+06
—
-
z/m
2。74646E+06
2。76181E+06
2。77329E+06
—
-
Vx(m/s)
—4.04563E+02
—3.22014E+02
-2。54830E+02
—
—
Vy（m/s)
—1。89569E+02
-2.18726E+02
—2.42544E+02
—
—
Vz(m/s）
4.62073E+02
3.49576E+02
2。57838E+02
—
-

将空间逐点交汇法得到的0号飞行器的轨迹和通过遗传算法确定参数后的拟合函数对应轨迹放于同一坐标系下对比分析如下：

图4 飞行器拟合轨迹和观测卫星轨迹对比图

图5 3D视图下的飞行器拟合轨迹和观测卫星轨迹对比图

图6 局部放大图

由图可知,首先观测轨迹和拟合轨迹的重合度较高，说明遗传算法参数寻优的优越性以及以及模型的正确性。但是轨迹开始重合度较高，后来差距较大。
以下给出预测以及拟合的0号空间飞行器的三个位置t—x、t—y、t—z和三个速度t—vx、t—vy、t—vz曲线对比示意图。

图7 飞行器观测及拟合轨迹t-x对比图

图8 飞行器观测及拟合轨迹t-y对比图

图9 飞行器观测及拟合轨迹t—z对比图

图10 飞行器观测及拟合轨迹t-vx对比图

图11 飞行器观测及拟合轨迹t—vy对比图

图12 飞行器观测及拟合轨迹t-vz对比图

5。3问题三的求解

本问题主要解决06号及其09号卫星观测带来的系统误差问题。在上一个问题中，我们采用逐点交汇定位的方法确定了00号飞行器的飞行轨迹。为了验证其是否能够估计系统误差，我们作如下推倒:
设为飞行器在基础坐标系的观测位置，为飞行器在基础坐标系中的实际位置，分别为观测卫星的系统误差，则有:
， ，其中的R表示上述的转换矩阵.
为飞行器在基础坐标系中的观测位置坐标,为飞行器在基础坐标系中的实际位置,即为上述公垂线交点坐标.由上式可以推算出二者的坐标对应于基础坐标系中的坐标，.
由公式，可以算出新的α，β值，分别为(），（α,β),则
。
而根据为观测量α，β在平面内的旋转量。

图13 平面示意图
则

图14 的分析
以上图为三个量的表示。可以认为是坐标系通过三个变化叠加得到坐标系，三个变化分别是:第一观测量的平移量，第二观测量的平移量以及观测量在平面内的旋转量即坐标轴在二维平面的转动角度。

5。4 问题四的求解：
5.4。1单星观测定位问题：
5。4.2多对多观测定位及系统误差的确定方法建模:
如果在多颗观测卫星观测多个空间飞行器的情况（假设有n个观测卫星同时观测m个待观测空间飞行器)，如下图所示，可以用以下方法处理数据,联合进行系统误差估计.

图15 多颗观测卫星观测多个空间飞行器示意图

图16 三维平面测位定向示意图
假设所有观测量都转化为统一的基础坐标系下，各观测站的观测噪声为不相关的零均值高斯白噪声，设观测站在基础坐标系中的位置为，观测视线在系统下的向量失量为，则其方程式为：
(30）
又设被观测空间飞行器坐标为.
考虑在空间直角坐标系下的点到观测视线的距离为：
(31)
则，被观测空间飞行器真实位置到各个观测站定位视线的距离平方和为：
（32)
其中,为加权系数。实际过程中可能出现型号不同或工况不同的观测卫星,他们的观测精度不尽相同，在这个情况下，可以用表示不同观测卫星观察结果的可信度。
欲得到被观测空间飞行器的位置,只需要找到使得取得最小值的。于是，只需上式分别对求偏导数即可。如下：
（33)
得到的解便是被测空间飞行器的预测位置（此时得到的是基础坐标系中的坐标值）。
类比第三问系统误差的估计，不难得到多对多观测时，系统误差的估计方法。
首先，由式（）得到的对于的位置估计，以此来反算任一观测卫星在观测坐标系中的“真值”，用解决第三个问题的方法来得到各系统的系统误差的估计值。
之后，通过观测不同的待观察空间飞行器，各个观测系统将会得到m组不尽相同的系统误差估计值。如，对于观测卫星，有以下三个向量：
；;
通过对三组列向量的简单处理，便可得到各个系统的系统误差分量并以此来修正观测值。

6 模型的优缺点及其改进
优点：
1. 我们通过建立相对简单的模型对材料所提出的问题给出了很好的解决方案
2. 运用智能遗传算法，并通过PID算法对对模型进行修改,使得模型预测结果与实际结果非常接近

缺点：
1. 没有考虑坐标系转换带来的转换误差
2. 对系统误差的估计不够准确，未能很好解决数据的白噪声误差
3. 未能求出基于单星观测的轨道模型

模型的改进：
对于原始数据的白噪声误差，可采用卡尔曼滤波对原始数据进行处理,以达到滤波效果.

- 25 -