2012年全国数学建模竞赛优秀选

（由填写第九届“杯”研究生数学建模竞题目 基于无源探测的空间飞行器 几何知识以及逐点交汇定位的思路，通过向量运算方法逐点定位目标，仿真06号和09号观测对0号空间飞行器的立体交叠观测轨道；然后选取线性质量变化方程m(t)m0kt，建立飞行器运动方程模型，通过分析和处理观测数据，拟合确定方程参重新建立选点定位模型，并设定指标w评价系统误差的稳定性，建立优化模型并求解。向交叉定位无源探测微分方程数值解轨道估计误差分目 .-3----------一、绪问题背景观测是探测空间飞行器发射与轨道参数的重要平台，按轨道特点，可分为高轨地球同步轨道和中低轨近圆轨道。其中，中低轨近圆轨道距地球表面数百在地球外层空间的惯性飞行段和再入大气层后的段。主动段通常由多级火箭相继推气层外，靠末级火箭关机前获得的速度在椭圆轨道上作无动力惯性飞行。段则根据任务需求，受控制后再入大气层，飞向目标。对于而言，在其结束前一直绕地再次证明，辐射源无源是现代防空系统、机载对地、对海以及对隐颗）观测的同步观测能够进行逐点定位，再结合空间飞行器的运动模型，便可以进发展技术与研究现状国内外现有的无源主要是利用目标自身辐射信号来对目标进行定位，关于表系统是洛克希德.马丁公司的“沉默哨兵”系统、英国防御研究局（DERA）研究备自主无源导航定位能力的国家。、等凭借本国的科技实力及军事战略需求也在极力地发展本国的无源导航，以求在未来的战场上掌握信息主导权。近年来，国内有关单位如国防、电子科技5129所、航天局、北航、西安电子科技大学、电子科技大学等都对单站、双站无源进行了和研究，并随着计算机科学技术、定位和算法、测量技术等的进步，无源定位和取问题提出主要工作建立模型首先，确立基础坐标系，即选取地心坐标系OcXcYcZc为基础坐标系。然后，利用建立改进的模型首先，利用空间立体几何知识进行坐标变换，即将观测坐标系OsXsYsZs下的位置矢量转化为基础坐标系OcXcYcZc下的位置矢量。其次，基于测向线公垂线的中点，利用双站测向定位原理，并根据仿真数据进行逐点定位，生成06号和09号观测对0首先，选择适当的vr(t和m(t来表示模型。其次，并利用最小二乘法对飞行器主动件，进一步导出飞行器的轨道估计方程。最后，再利用数值法按照从50.0s170.0s通过建立模型计算求出偏移量d、d，通过分析d、d，断定选取公垂线的中点对重新选取测向线公垂线上的点F，通过逐点交叉对飞行器的轨迹进行重新定位。考虑到测向线公垂线长度相对于测向线长度极小，公垂线长度视为固定值，在公垂线上选取该F，使线段MF与线段FN成线性比例。选定该点后，进行偏移量以建立误差参数求解模型，确定误差函数中的两个参数μ和，确定误差函数，二 、符号说明及基本假2.1基础坐标系为随地心平移的坐标系，取地球中心Oczx轴由Oc0y轴，建立直角坐标系OcXcYcZc。地心Oc在绕日椭圆轨道上运动，所以理论上OcXcYcZc系是非惯性系。随运动的观测坐标系OsXsYsZs2，原点取为中心OsXs轴沿OcOsZsXs垂直指向正北，Ys轴按右手系确定。由于一般测的观测坐标系也叫做UEN坐标系，因为三个坐标轴分别指向上（UP、东（EAST）和

2观测坐标系示意新坐标轴OX，OY，OZ的方向余弦分OX——l1m1n1OY—l2m2n2OZ——l3m3方向交汇点进行计算无源定位。2.2符号说明FT：m(t)：质量变化r(t：r(t对时间tGm：地球引力常数（本题中地球引力常数取Gm m3/s2vr(t)： Ys(t0)：零时刻探测在基础坐标系下y轴方向的位Ys(t0)：零时刻探测在基础坐标系下y轴方向的速度Zs(t0：零时刻探测在基础坐标系下z轴方向的速度2.3基本假设观测的运行轨道相对于地球是不变由于空间飞行器的主动段主要依靠第一级火箭的“垂直段+程序拐弯段（加外三、问题分问题一分析 ：在仅考虑随机误差的条件下，根据数据文件satinfo.txt中的观测

|r(t)

r 计算09号观测在50.0s、100.0s、150.0s、200.0s、250.0s五个时刻的三维位置。本题需要建立轨道模型，推导出观测和空间飞行器的轨道运动方程，然后进行求体，则它对人造的吸引可等效为一个质点，相当于全部质量集中于地球质心这一点上，于是地球和就构成一个简单的二体系统。因此我们选择地心为坐标系的原点，来讨论人造相对地心的运动，这是很自然的。但严格来说，该地心坐标系并非惯性问题二分析题目：在本题给定的仿真数据下，06号和09号观测对0号空间飞行器形r(t)r(t)FFGmr(t)v(t) |r(t)rys

xsyszs标变换法则，将观测坐标系OsXsYsZs下的位置矢量转化为基础坐标系OcXcYcZc下的中给的仿真数据（编号为i的在观测时刻t对编号为j的飞行器以及对应观测数然后，选择适当的vr(t)和m(t)来建立轨道模型，并利用最小二乘法对飞行器问题三分析统误差影响比较大的三轴指向误差，即三个常值小量d,d,d，分别表示第一观测量的平移量、第二观测量的平移量以及观测量在平面内的旋转量。问题二中在仅考虑随机误差的情况下，我们基于S1、S2的测向线公垂线的中模型计算求出d、d重新在公垂线上选取点，通过逐点交叉对飞行器进行重新定位，保证d、d（这是由系统误差的稳定性决定的又已知是由观测坐标系下的两个无量纲比值确定的值（ys

zsxsys,zs为空间飞行器在观测坐标系中的坐标。也就是说，系统误差d、d的确定问题四分析题目中要求对只有09号观测单星观测的01号空间飞行器进行轨道估计。该题涉及到单星无源。单星对星的无源定位，就是把无源侦察设备安装到平台上，利用单颗平台截获目标发射的通讯、测控和导航信号并进行处理，进四、问题一的建模、求的简化运动方程r(t)

|r(t)

r(t)以及零时刻观测的位置和速度x,y,z,x,y,z，立观测的轨道模型进行求解4.1模r(t)r(r(t)r(

(t)

GmX|r(t) Ys(t)

|r(t)G

Ys

Zs(t) mZs|r(t)r(t)X(t)Y(t)Z 其中（Xs(t)Ys(t)、Zs(t)为观 在地心坐标系Oc-XcYcZc下的坐标分量Xs(t) Ys(t)、Zs(t)分别为在三个坐标轴方向的加速度。Xs(t0、Ys(t0、Zs(t0、Xs(t0、Ys(t0、Zs(t0为微分方程组的初始状态参数，也就是题目中已给出的初始位置和速度，rr(t Xs(t)Ys(t)Zs(t222.3 A(t) mXs |r(t) B(t) mYs |r(t) C(t) |r(t)

Zs Xs(t0、Ys(t0、Zs(t0Xs(t0、Ys(t0、Zs(t0为微分方程组的初始状态参数，这1取步长h=tt0t2t1t3t2...tn1tn（h为常数，根 法的计算，则1

)GmX(t)A(t ||r(t) X2(t)Y2(t)Z2(t

|r(t) s

根据以上数值分析方法，我们选取以下步长h，分别计算得到09号观测在xyz1.77380E8.16139-1.501628.126784.68469-1.227708.082724.84723-9.523538.029295.00415-6.758987.966545.15526- 图4-1-1轨 测的初始位置（.166765，.631471，.714791）运行轨迹误差，见表4-1-2。步运行一圈后和初始位置的坐标差值xyzyz4.34325E2.03606E-4.71839E--4.86518E-1.84732E-4.67887E-7.92706E-5.33559E-4.43129E-4.13038E-

图4-1-2轨道误差显示好于步长为0.050.5的。因此，步长的选取不易过大，亦不易过小。道，如图4-1-8所示。误差产生的原因，如图4-1-9所示。 误差产改进的模使计算产生的误差尽量减小，如图4-2-1所示。1具体分析：将运动轨迹的变化微元AB转化为微元AB'（AB'为运动在微元时间△t内的变化微元轨迹，该弧线的切线方向为r，所在圆心为O，O为地球圆心1这样' 即为减小的矢量差，从而减小误差，从而提高观测轨道的精确度，减小计算结果误。

rA0BrA0BH1在改进轨道模型中，向量OB'为观测 的位置矢量r， 只要推导出向量OB'的坐标，后续求解过程便类似于上述 模型，此处不再详细进行模型的求解。下面，我们主要进行向量OB'的推导。1已知飞行器加速度r0同位置矢量r0rrGm0|r00根据圆周运动规律，得到AB'所在圆的半径R，即B'O'2r02r0/

根据向量叉乘原理，易得AO'的方向向量nnnr(r令BOO'，根据向量三角函数，得的余弦

cos OBO'B'2 ' OB OO O OBOBOHO'H为△OBO'的垂线，求得该垂线的高

O'H OO'2OH在△B'O'H OO'2OHHB' R2HB' R2O'H'R2 22OB'OB'OHOB'7OB'OB根据以上数值分析方法，我们选取以下步长h，分别计算得到09号观测在xyz1.77380E8.16138E4.51669-1.50162E8.12676E4.68468-1.22769E8.08269E4.84721-9.52346E8.02925E5.00412-6.75890E7.96650E5.15524-xyzyz-2.47084E-2.36085E--6.63349E-模型小结通过模型一和模型二分别求出观测的运行轨迹，并给出各模型条件下09号观测在50.0s、100.0s、150.0s、200.0s、250.0s五个时刻的三维位置。通过对模型的五、问题二的建模、求基于测向线公垂线中点的两站测向交叉定位模型标系OsXsYsZs下观测对空间飞行器进行观测的数据，而空间飞行器的主动段的运 r(tFeFT r(tvr 是建立在基础坐标系OcXcYcZc下的方程，因|r(t) 现根据问题二，我们已知飞行器在观测坐标系OsXsYsZs下的位置，下面将对飞行

OcOs，不妨设观测坐标系OsXsYsZs的原点为（X0YoZo，则平移向量OcOs=（X0，YoZoOSYS、OSZS根据三维坐标系下的向量运算以及右手法则，分别求出坐标轴OSXS、OSYS、OSZS的Xs轴沿OcOs连线，离开地球方向为正，可得OSXSOX( zs=(l,m,n, ,

'OSZsOcOs(OcZcOcOs将向量OY'进行单位化，得到OY的余弦向量（不妨设为（lmnS S

OSZs(l3,m3,n3'OSYsOsZsOs进而，求得OSYS的余弦向量OSYs(l1,m1,n1将新坐标轴OSXS、OSYS、OSZS的方向余弦，带入旋转变换，得到坐标旋转变Xcl1Xsm1YsYlXmYn 2 2 2ZlXmYn 3 3 3Xcl1Xsm1Ysn1ZsXYlXmYnZ 3 3 3 无源探测器的观测常采用红外光学探测器，单个红外光学探测器不具备测距能探测和09号探测的两条测向线可能不会交汇于被探测空间飞行器这一点，大逐点交叉定位，模型如图5-2-1所示。其中，S1、S2分别为06号、09号 n其 n0

MN(S1S2n)nn

OS'OS sinS'NS sinNS'S

f1f1f1f112sinNSS'2

1f 1f2S fS

SMSSM

1sinNSS 112 sinS' 2112SNS'SN12 sinS'12

sinNS 347的长度，进而得到向量S1MS2NSMSMSM11f1SNSN22f2又根据向量和，分别得到向 、OMOS1ONOS2中垂线中心向 2飞行器运动方程模型选取火箭喷气速度vr(t)和火箭质量瞬时质量m(t)为一随时间变化的线性方程：m(t)=m0-kt 则火箭质量变化率m(t)为一常系数k，即

m(t)

（h=1s（ABBC）进行研究，则我们可以假设火箭在每个微元ds段内进行匀加速直线运FT为火箭产生的推力加速度大小，Fe为飞行器所受的外力加速度大小，方向是中线指向O，大小为固定值。设A点处速度为v1，B点处速度为v2FGm A r rA

FTvr

2

FvrTm0Tk

ABv1v22

根据力学几何关系得出v2与v1的关系v2(FTFecosABh+v1 cosAB(rArBAB。联立式（5-2-7（5-2-8）消除 ，v(FGm )h+ TT

Fˆ

A rABFˆ(2v22BCGmcosBT r 进行计算仿真，得到火箭的推力加速度拟合曲线，如图5-2-1.5-2-1s：F 308.6

r(t)Fr(t)FFGmr(t)v(t) |r(t)rr(t)FFr(t)FFGmr(t) |r(t)7863308.6t5.2.1..

XX (t)X

f|r(t)f

(t) 308.6.

|r(t)Yf(t)

GmY(t)7863

Yf

|r(t)

308.6

|r(t). fZf(t)

Z(t)7863Zf

|r(t) 308.6t|r(t)(t)Yf(t)Zf Xft0、Yf(t0、Zf(t0、Xf(t0、Yf(t0、Zf(t0为微分方程组的初始状态参数 就是题目中已给出 初始位置和速度，且 Xf(t)Yf(t)Zf(t).AXfBYfCZ

Af Af(t) mXf(t) |r(t)

308.6t|r(t) 7863Bf(t) B(t) mY |r(t)

308.6t|r(t)

7863C

3ffC(t) m3ff |r(t)

(t)

308.6

|r(t)这样将二阶微分方程转换为一阶微分方程。取步长h=t1t0t2t1t3t2tn1（h为常数，根据法的计算，则X(tn1)hA(tn)X(tnAA(t)h(GmX(t)|r(t)n308.6 rA(tn))A(tnn时间xyzvxvyvz------------- 速度矢量修正模型元向量FR，其单位向量为 ，我们的目的就是要将速度微元向量FR修正到 FR0f0FR易得速度微元向量FR、FR0的单位向量 fff0fiff那么，速度方向的修正因子向量dddf0i从而得到速度修正的单位向量vv1iv0i进而得到修正后的速度向量v1vv1vf容易得到速度大小的修正因子大小sfssv1v1v0/测数据计算出的运行轨迹进行比较，如图5-3-2所示。

六、问题三的建模、求统误差影响比较大的三轴指向误差，即三个常值小量d,d,d，分别表示第一观测量的平移量、第二观测量的平移量以及观测量在平面内的旋转量。问题二中在仅考虑随机误差的情况下，我们基于S1、S2的测向线公垂线的中模型计算求出d、d重新在公垂线上选取点，通过逐点交叉对飞行器进行重新定位，保证d、d（这是由系统误差的稳定性决定的又已知是由观测坐标系下的两个无量纲比值确定的值（ys

zsxsys,zs为空间飞行器在观测坐标系中的坐标。也就是说，系统误差d、d的确定坐标反变换模型坐标，我们可以将其简化为矩阵形式的变换（6-1-1） X0

YY c 3s 0 其中，（Xc，YcZc）为地心坐标系下的坐标（Xs，YsZs）为观测坐标系下的坐标（X0X0X0）2 Xs l31 Y s 3 n3ZcZ0设观测S1在观测坐标系下的坐标为S1（Xs1，Ys1，Zs1），S1的测向线与公垂线的交点M，M在地心坐标系下的M（XMYMZM，公垂线中F的坐标为F（XFYFZFM在观测坐标系下的坐标为Ms1（XMS1YMS1ZMS1，又可表示成Ms1（1，0 ）则根据坐标反变换，得出式（6-1-Ms1TQ1(MTST Q

l3m m m 。m m m 。 n3系统误差判断模型三轴指向误差中的两个偏移量d、d分别为式（6-2-1（6-2-d0d0

通过对偏移量d、d及旋转量进d行线性拟合，得到如下仿真图（1）06号 06号d的分布 06号d的分布（2）09号

06号d的分布 09号d的分布 09号d的分布 09号d的分布观察到两个偏移量d、d以及旋转量d重新选取测向线公垂线上的点F，通过逐点交叉对飞行器的轨迹进行重新定位。公垂线选点定位模型时间变化的曲线图，发现公垂线长度的变化范围为2200m-2400m，而飞行器运行轨道400000m2400<<400000，因此可以近似看作为固定值。这就对以后公根据上述判断，我们认定测向线公垂线的长度不变，不妨将三内公垂线与测向线组成的曲面展开成二内的矩形平面，如图6-3-1所示。从50s开始，每隔一每条公垂线上选取该F，使线段MF与线段FN成线性比例，该比例从μ变化到。

时 119+1 （1-n）+n 根据以上推导，容易得到（50+n）s时刻下观测坐标系下的向量OF50nOF50n=（1-'）OM（1-n）+n 联合上述两个模型，得知观测S1、S2在二维观测平面上的偏移量d1、d1、d2、d2以及所在平面的旋转量d是关于自变量t的函数，其中包含参数μ、，形如d1(t;u)。发现d1d1d2d2d1d2近似围绕常数上下波动，符合系统误差的性质（1）06号 06号d的分布 06号d的分布（2）09号

06号d的分布 09号d的分布 09号d的分布 09号d的分布误差参数求解模型接下来，进行参数u和的求解，根据第三个模型得到d1d1d2d2均是关于u和的函数。令k1d1u、k2d1u,)、k3d2u,)、k4d2u,)，为使系统误差尽量减小，即偏