2012年优秀参赛集题

参 （ 填写题目 基于无源探测的空间飞行器主动段轨道估计 要空间飞行器的快速探测和实时对于国家安全具有重要的意义，探测器发射和轨道参数是进行实时并作出后续反应的基础。在实际侦测中，常采用多颗中低轨近圆实现组网探测，以快速发现空间飞行器。构建的空间飞行器的轨道方程有利于实时估计和其运行状态，本文主要考虑空间飞行器在重力斜飞段的后对于问题一，首先根据题中给定的观测的简化运动方程；再采用常用的常微分方程组的数值解法-龙格·库塔积分方法计算出了09号观测的运动轨迹，并50.0100.0150.0200.0250.0s（1·库塔积分方法的优缺点。（3（9；接着使用0050.0s170.0s27，8，够同时估计系统误差和轨道的数学模型计算出了6号9号观测的系统误（见4（23）59最后了扩展卡尔曼滤波模型的优缺点。当同时有多颗观测观测多个飞行器的情况下，可以利用半参数估计法联合计算出空间飞行器的快速探测和实时对于国家安全具有重要的意义，探测器发射和轨道参数是进行实时并作出后续反应的基础。中低轨近圆是一种探1、2观测在任意时刻的位置计算是估计的前提，请根据satinfo.txt和观（2，6在本题给定的仿真数据下，06号和09号观测对0号空间飞行器形成了立体交叠观测，请结合几何知识按照逐点交汇定位的思路，给出0号空间飞行器在公式（1）框架下的轨道估计，注意选取适当的vr(t和m(t50.0s170.0s10.0s0号空间飞行器在各个采样点的位置位置t-x、t-y、t-z和三个速度t-vx、t-vy、t-vz曲线示意图。若06和09号两颗观测均有可能带有一定的系统误差对系统误差进行确的估计能够有效提高精度。利用上述的逐点交汇方法能否同时对系统误差进行估计？若不能，是否还有其他的思路能够同时估计系统误差与轨道？给出你的解决方案与估计结果。在报告中除给出与第二问要求相同的结果外，还应分别给出两颗观测卫星的系统误差估计结果，共六个数值，分别是两颗的d,d,d。对只有09号观测单星观测的01号空间飞行器进行轨道估计，结果形式要问题分析与坐标系构程组的数值解法计算出09号观测的运动轨迹。050.0s内采样点（10s一次采样）方程即可。当同时有多颗观测观测多个飞行器的情况下，可以通过半参数估计法（1）基础坐标系与观察坐标系的建为描述观测和空间飞行器的运动，需要建立适当的坐标系。本题基础坐标系为随地心平移的坐标系，取地球中心Ocz轴，指向北极为正xOc指向零时刻的0经度线，再按右手系确定y轴，建立直角坐标系OcXcYcZc。地心Oc在绕日椭圆轨道上运动，所以理论上OcXcYcZc系是非惯性系。1第二个坐标系是随运动的观测坐标系OsXsYsZs2，原点取为中心OsXs轴沿OcOsZsXsYs轴按右手系确定。由于一般测量的轨道都不会严格经过南北极上空，所以这种坐标系的定义是明确的。如此定义的观测坐标系也叫做UEN（UP图 观测坐标系示意

飞行器的瞬时质量m(t)3相对于火箭尾部喷口的喷射速度的矢量vr(t)的方向与飞行器的速度方向反 附件中给定的观测初始参数中可以不考虑误差名符单名符单m(tm(tkgOs-XsYsr(tmm3/r(t对时间r(tvrv(t问题一求根据题意，观 的运动方程可简化如下r(t)r(t)Fe|r(t)r其中向量Fe表 在基础坐标系下的位置矢量r(t表示r(t对时间tGm为地球引力常数（常数取Gm3.98600510m/svr(t14 设观测在t时刻的观测位置为(x(t),y(t),z(t))，则有r(t)=(x(t),y(t),z(t))。于是，drr

，初值为

r(t

r

|r(t)步法。针对运动的特点，在大多数轨道计算程序中通常选择龙格-库塔（Runge-Kutta，简称为RK）单步法、阿达姆斯（Adams）和科威尔（Cowell）多步AdamsCowell方法用于解二阶常微分方程[1,2]。本文选择龙格-库塔（RK）求解观测在给定时刻的三维位置。对于一个时

观测的速度向量r(t0+nh)；再将求解结果带入方程组（2）中的第一行方程，求解可得观测的位置向量r(t0+nh)，即观测的三维坐标。 进行计算调用其预定义的ode45命令求解，积分步长取为10s。最终计算出的观测 在50.0s、100.0s、150.0s、200.0s、250.0s五个时刻的三维位置的结果如表1所示（结果保留6位有效数字。表1：观测在5个时刻的三维位置（基础坐标系时刻观测的三维位置（106XYZ模型分对于求解常微分方程组，可以采用解析法或数值法。数值方法通常较解析法具有对于问题一，采用了龙格-库塔（）积分器进行计算，该方法的优点是其可以自起步积分，这是多步积分器所不具有的独有特性。积分器的积分步长与积分阶数的选择依赖于所需要的精度和微分方程右函数的性质，通常将积分步长与积分阶数作为问题二求问题二要求基于附件中给定的06号和09号观测对0号空间飞行器形成了立的vr(t和m(t)050.0s170.0s（10s一次采样0号空间飞行器的三个位置t-xt-yt-z和三个速度t-vx、t-vy、t-vz曲线示意图。vr(tm(t)的表示模m(t)量变化模型[2]为m(t)=m(1m(t)tm0m00m0行器的初始质量。由于题中给出了相对于火箭尾部喷口的喷射速度的矢量vr(t)与飞行器的方向反向线接近或相同，其大小一般较为稳定，因为可以将vr(t的大小视为常数，在任意时刻仅是方向发生变化。于是vr(t)和m(t)的表示模型可设为下列表达0m(t)=m(1-m(t)t0

v(t)=-|v(t)|

其中v(tv(t)|v基于角度信息的定位模为了确定飞行器在基础坐标系和观测坐标系中相对于两颗观测的方位角和仰角的关系，现将基础坐标系的原点分别6号、9号观测的中心得出O1c-X1cY1cZ1cO2c-X2cY2cZ2c两个坐标系。设空间飞行器在基础坐标系中的坐标为(x,y,z)，其对应的在6号、9号观测坐标系中坐标记为(x1,y1,z1)和(x2,y2,z2)，而对应于O1c-X1cY1cZ1c和O2c-X2cY2cZ2c两个坐标系中的坐标记为(x1c,y1c,z1c和(x2c,y2c,z2c。6号、9号观测在基础坐标系中的位置分别为(x'1,y'1,z'1)、(x'2,y'2,z'2)，可根据第一问建立的模型求解得到。飞行器相对于6号观测的方位角和仰角分别为1、1， 对应的观测值为m、m，上述仰角均基于基础坐标系。坐标系之间的变换关系如下 x1 x1c x2 x2cy1A1y1c，y2A2y2c

z z 1 1c 2 2c其中变换矩阵A1A2cosB1Scos

sinB1S sin

cos sinBcos

sinBsin cosB

1ScosB2Scos

cosB2Ssin sinB2S sin

cos sin cos

sin sin cos

2S式中的BS,LS分别为观测相对于基础坐标系的方位角和仰角。变换矩阵A1、A2可由通过先先求出的三维位置，再求出其对应的方位角、仰角即可计算出。6号、9号观测的观测方程[3]如下：m=+n=arctanx-x'i

m=+n=arctan +n(x'-x)2+(y'-(x'-x)2+(y'-ii式中ni1和ni2cos x2x2+y z

sin= x2x2+y2+z 根据题中给定的无量纲比值i=i xi

i=i sinsin cossin

x

，X=y

sin2

z sin2sin cos2sin cos2 x'sinsin+y'cossin-z' Y=

1x'2cos2-y'2sin x'sinsin+y'cossin-z'cos 2在（8）式中，只有X为未知数，其余参数均可利用方。于是，如果矩B和YBm和Ym（即用观测值mm代替B 和Y中的i、i。若仅考虑Y的噪声干扰，并假设噪声的各个分量统计独立，则可BTBT

X=（BT

飞行器的运动轨迹方程（最小二乘法r(t)r(t)FF |r(t)r(t)v(t)r结合选取的vr(t和m(t)的表示模型（3）和飞行器角度定位方程（9）rr(t)|r(t)r(t)v(t)rm(t)=m(1-m(t)t0 0v(t)=-|v(t)|v

x

|vr(t)=X=y=（BT

zz x

rr|r(t)r|vr(t)|m(t)v|v(t)|(mr(t)=y=（BTzz 由于的喷射速度vr(t)大小相对稳定，故|vr(t)|可视为常数，并根据所选取m(t)m(t)、m0模型求察附件中的数据，发现两颗观测的采样时间并不同步，但采样间隔均为 的采样时间进行同步本文采用了曲线拟合方法计算出无量纲比值和随时间的变化方程，再根据拟合方程求出在给定时刻和的观测值，拟合的图3:6 拟合得出的和图4:9 拟合得出的和从拟合结果中可以看出，采用三次曲线拟合已经可以很好地近原始数据曲线，故在本文中采用了三次方程计算给定时刻的6号、9号观测得到值与时间t6 ：1=-8.2263*10-9*t3-7.5502*10-7*t2+9.0069*10-19 ：

=-1.1348*10-8*t3-3.7316*10-9*t2+5.3769*10-4*t-=-3.7190*10-9*t3-4.3284*10-7*t2+8.7754*10-（12，即r(t)=

|r(t)

r

|vr(t)|m(t)v|v(t)|(m0- xr(t)=y=（BT|vr(t)|

zzm(t、m0为常数。从方程组中可知，可将m(t、m0为一个整体来计算，即只需计算|vr(t)|和m0/m(t)。根据两 的观测数据，通过（x,y,z)；再根据公式（12）程，可进一步计算出飞行器的速度。由于|vr(t)|、m0m(t)具体表达式，需要先求出|vr(t)|m(t、m0的三维位置可以拟合出飞行器的在对应时间段内的位置曲线通过拟合发现，方程，需要根据拟合方程求出位置的常数|vr(t)|和m0m(t)，选取t=50s和t=60s两个时|vr(t)|= 0r(t)=

|r(t)

r

上式（13）50s170s（10s）之2所示。飞行器的速度XYZXYZt-x、t-y、t-zt-vx、t-vy、t-vz5-7所示，从图中可以看出，x、y、z的大小均随时间而增大，其对应的速度大小也随时 5：t-x、t-vx6：t-y、t-vy7：t-z、t-vz残差和标准差分由于残差是指观测值与值（拟合值）之间的差[3]，即是实际观察值与回归估而标准差（即均方根差）

eiXi-ˆi

(Xi-ˆ

)2/在本文中Xi对应于根据观测值，采用角度交汇定位和最小二乘法计算得出的飞行器三维位置。由于飞行器的速度不能通过观测直接得出，故在此只50s170s（10s）的残差如833（106（106XYZXYZ000a)x的残差曲 b)y的残差曲c)z的残差曲线 d)位置残差随时间的变化关系图8：残差曲线飞行器轨道估计的标准差（即均方根差）

(Xi-ˆ

)2/在求解飞行器动力学轨迹方程时，需要给定初值。首先根据附件的数据，初始轨道信息，但其不能提供预报轨道，而动力学法定轨虽考虑了飞行器的受力情况，能提供高精度的轨道和预报轨道，但是要求作用于的力学模型十分精确，则任何力学模型误差都将带人历元状态参数估值中，一般来说观测量离解算历元越远，动力模型的误差影响越大。本文在计算飞行器动力学轨迹方程时，通过最小二乘法计算得出其初值，即通过XBT1BTY来计算初值。然而，该方法仅考虑Y的噪声干扰，并假设噪声的各个分量统计独立。若同时考虑Bm和Ym的噪声干扰，并且假设噪声的各个分量统计独立，则可以得到总体最小二乘意义下的解为

v42/v4v(3)/v(1)

v(4)/v(1)

Bm进行奇异值分解后最小奇异值所对应的又奇异向量，v4i向量v4的第i个元素YmBm中的噪声分量并不统计独立，为了得到它们之间的

i1

m m并考虑如下Taylor

i sinm

i sin

cosmo(n

coscosm cosm

sinmo(n

sin i

sinmn i isin i i

cosmo(n

icoscosm icosm

sinmo(n i

iBBmFNFNFN

44F1diag[f11F2diag[

f1Nf2N ，f

sin iiFdiag[ ]，

cosmsin sinmcos

iF4diag[

cos

f4N

f=

，

= sinmsin cosmcos sin i i

xsinmycos =(ysinmxcosm

(xsinmycosm)cosmzsinm iN N从（14）式中可以看出，YmBm所受到的噪声干扰均来源于噪声向量N，于是可 NR2N,XBm

YmX

XF1NF2NF3NF4N1 式中W是个对称矩阵，显然（19）式是一个有约束的优化问题，该问题可转化为TXT

1T)1

X

0 XX

X(1)FX

X(3)FX(4)F，

YmX(iX的第i 0（ 0问题三求问题三要求基于问题二的求解结果，并考虑观 的系统误差，判断能否利果，共六个数值，分别是两颗的d,d,d。 d,d,d来表示，分别表示第一观测量的平移量、第二观测量的平移量以及观测量在平面内的旋转量。同时估计系统误差与轨道的模设飞行器的状态向量（包括位置和速度）XX=(x,y,z,vx,vy,vzX(t)f(X(t))Z(t)h(X Xkk1Xk-1Wk-

ZHX h(X其中{Wk1k=1,2}、{Nkk=1,2}是随机过程噪声序列，其方差分别为Qh(X为测量值向量。

k

Hk

Xk

可比矩阵在k-1k时刻的值。WkNk扩展卡尔曼滤波[6]算法的执行过程分为两个阶段：和更新。在阶段，利用当前状态和协方差计算下一时刻目标状态的先验估计；在更新阶段，利用到的测量值更正先验估计以活动改进的后验估计。扩展卡尔曼滤波算法以“①利用目标的当前状 下一时刻的状ˆ TXkXk1f(tk1)Tk1f(tk1)②下一状态的误差协方PPT kk PHK k HPHk XˆXˆK(ZHXˆ P(IKH)P(IKH)TKRKT k XZH分别为状态方程和QRK尔曼增益，I为单位矩阵，P为估计的误差协方差。

(X)I(X)T

XT222Xkk

0Xk1 I k

k1

XZk Z

kkkkkk式中 Tkk参数对应的矩阵，其余符号的表示意义同前文。状态方程的雅可比XHXEk的具体表达式如下（推导过程由于过状态方程的雅可比X(X)f(X) I 0I

3G 3G 3Gxz A 3G 3A 3G 3G 3Gz2 HX 0 H(X)= 0 0 观测值中的三维坐标值可以通过(xyz)T（BT1BTY由于ddd分别表示第一观测量的平移量以及观测量在平面内的旋转量。在本文中，将旋转量理解为(点平面坐标系原点连线与坐标轴的夹角。于是有dim dim,i1, dim- 式中m、m、m表示观测值，、、表示真实值。同时，给定的无量纲比值

i=i xi=模型求解及残差分

i i通过编程实现扩展卡尔曼滤波算法，计算得出6号和9号的系统4所示。表4:两颗的系统误d（10-4d（10-4d（10-469号飞行器轨道估计的标准差（即均方根差）

(Xi-ˆ

)2/ 50s170s（10s采样间隔）5所示，其运行状态示意9所示。5飞行器的三维位置（106飞行器的速度XYZXYZa)t-x曲线示意 b)t-vx曲线示意c)t-y曲线示意 d)t-vy曲线示意e)t-z曲线示意图 f)t-vz曲线示意图图9：飞行器的运动状态示意图模型分析与态参数的更新。所以扩展卡尔曼滤波算法一般只需前一个历元的状念参数估值，无须所有历史观测信息，具有很高的计算效率，并可以进行实时估计。应用扩展卡尔曼滤波算法可以充分考虑飞行器的运动信息，在确定轨道的同时还可以估计出的运行速度，可以提供短期的预报轨道。当某个历元缺失观测值时，几何法定问题四求问题四基于问题三求出的9号观测的系统误差估计结果，利用09号观测01由于问题三已经求解出了9号观测的系统误差，在求解01号空间飞行器道的过程中，只需考虑随机误差了。附件给出了9号观测所观测01号空间飞行器的和，可以计算出相应飞行器相对9号观测的方位角和仰角；根据和角度变化率建立观测模型。观测量为被观测相对于观测的方位角和仰角。2个角度实际就是在地固坐标系下观测到空间飞行器的位置矢量的方位角和仰角，无须做较为复杂的转换。设(t)和(t)分别为方位角(t)和仰角(t)的角度变化率，因(t)则

，

(t)x2(t)y2(t)(t)

xsinycos

(t)xsincosysinsinz x2y令

h(X(t))

(t)(t)

x2y2h(X(t))

(t) h4(X(t)) (t)43X(t)(x,yzx,yz)43H

H

1

sincos

(cos2sin2)y2cossinx sinsin (cos2sin2)x2cossinyH2

HH 3

/r/