（由填寫(xiě)題目無(wú)源探測(cè)的飛行器主動(dòng)段軌道估：本文首先通過(guò)直接求解的簡(jiǎn)化運(yùn)動(dòng)方程，得到06號(hào)與09號(hào)兩顆標(biāo)軌道的方法。該方法以逐點(diǎn)交匯定位法為基礎(chǔ)，利用測(cè)量值誤差在小角度條件下轉(zhuǎn)化為測(cè)量軸偏轉(zhuǎn)角的方法估計(jì)系統(tǒng)誤差，再利用該系統(tǒng)誤差修：雙星單行彈道估計(jì)最小二乘估一、問(wèn)題重有些國(guó)家會(huì)發(fā)射特殊目的的空間飛行器，如彈道式、偵察等。對(duì)他國(guó)發(fā)射具有敵意的空間飛行器實(shí)施并作出快速反應(yīng)，對(duì)于國(guó)家安全具有重要的戰(zhàn)略意義。發(fā)現(xiàn)發(fā)射和探測(cè)其軌道參數(shù)是實(shí)現(xiàn)和作出反應(yīng)的第一步，沒(méi)有觀測(cè)，后續(xù)的判斷與反應(yīng)都無(wú)從談起。居高臨下，是探測(cè)觀測(cè)按軌道特點(diǎn)，可分為高軌地球同步軌道和中低軌近圓軌道衛(wèi)3.6萬(wàn)千米，軌道平面與地球赤道平面重合，理論上用3顆間隔120度分布的同步軌道可覆蓋地球絕大部分表面。中低道平面交成一定角度，且常由若干顆實(shí)現(xiàn)組網(wǎng)探測(cè)。裝置于上的探測(cè)同時(shí)具備定向和測(cè)距兩種能力；無(wú)源探測(cè)器則接收目標(biāo)輻射。采用無(wú)源探測(cè)器的觀測(cè)常采用紅外光學(xué)探測(cè)器，只接收目標(biāo)的紅外輻射信息，可定向過(guò)數(shù)學(xué)模型與計(jì)算方法作出的。當(dāng)觀測(cè)飛行一段時(shí)間，探測(cè)器測(cè)得目標(biāo)相對(duì)于運(yùn)動(dòng)的觀測(cè)數(shù)據(jù)，以觀測(cè)和空間飛行器的運(yùn)動(dòng)模型和觀測(cè)模型為外層空間的慣性飛行段和再入大氣層后的段。主動(dòng)段通常由多級(jí)火箭相繼飛行。段則根據(jù)任務(wù)需求，受控制后再入大氣層，飛向目標(biāo)。對(duì)于而言，在其結(jié)束前一直繞地飛行，故無(wú)段。料，箭體應(yīng)盡快穿過(guò)稠密大氣層，故火箭一般先垂直發(fā)射。設(shè)A點(diǎn)為地面發(fā)射點(diǎn)，AB為垂直上升段，BC弧段為程序拐彎段，CD弧段為重力斜飛段，DE垂直緯 橢圓軌赤道平1坐標(biāo)系為隨地心平移的坐標(biāo)系，取地球中心Oc為原點(diǎn)，地球自轉(zhuǎn)軸取為z立直角坐標(biāo)系OcXcYcZc。地心Oc在繞日橢圓軌道上運(yùn)動(dòng)，所以理論上2第二個(gè)坐標(biāo)系是隨運(yùn)動(dòng)的觀測(cè)坐標(biāo)系OsXsYsZs，見(jiàn)圖2，原點(diǎn)取為衛(wèi)星中心OsXs軸沿OcOsZsXs垂直指向正北，以這種坐標(biāo)系的定義是明確的。如此定義的觀測(cè)坐標(biāo)系也叫做UEN坐標(biāo)系，因?yàn)槿齻€(gè)坐標(biāo)軸分別指向上（UP、東（EAST）和北（NORTH）三個(gè)方向。

3r(t)vr(t) |r(t) 其中向量Fe表示飛行器所受的外力加速度之和，F(xiàn)T表示火箭產(chǎn)生的推力表示r(t對(duì)時(shí)間tGr常數(shù)（本題中地球引力常數(shù)取

m3/s2

vr

3r |r(t)

vrvrvr和m(t)的模型不同，m(t)一般而言應(yīng)為嚴(yán)格單調(diào)遞減的非負(fù)函數(shù)vrys

,考慮與平臺(tái)相關(guān)的系統(tǒng)誤差，即不同觀測(cè)的系統(tǒng)誤差相互沒(méi)有關(guān)聯(lián)，同一觀測(cè)對(duì)于不同空間飛行器的系統(tǒng)誤差是一樣的。經(jīng)由適當(dāng)?shù)暮?jiǎn)化模型，用三個(gè)常值小量ddd來(lái)表示，分別表示第一觀測(cè)量的平移量、第二觀測(cè)的平移量以及觀測(cè)量在平面內(nèi)的旋轉(zhuǎn)量。單個(gè)紅外光學(xué)探測(cè)器不具備測(cè)距能力，但借助多顆（含兩顆）觀測(cè)的參數(shù)文件satinfo.txt用來(lái)觀測(cè)信息，每行表示一顆，包含六表示編號(hào)為i的對(duì)編號(hào)為j的飛行器的仿真觀測(cè)數(shù)據(jù)信息，按照時(shí)間順序分行，每行分三列，分別是觀測(cè)時(shí)刻t以及對(duì)應(yīng)觀測(cè)數(shù)據(jù)觀測(cè)在任意時(shí)刻的位置計(jì)算是估計(jì)的前提請(qǐng)根據(jù)satinfo.txt和觀測(cè)的簡(jiǎn)化運(yùn)動(dòng)方（2計(jì)算09號(hào)觀測(cè)在50.0s100.0s150.0s200.0s、250.0s五個(gè)時(shí)刻的三維位置。結(jié)果保留6位有效數(shù)字?？臻g飛行器在（1）框架下的軌道估計(jì)，注意選取適當(dāng)?shù)膙r(t)和m(t)的0號(hào)空間飛行器的三個(gè)位置t-x、ty、tz和三個(gè)速度t-vx、t-vy、t-vz若06和09號(hào)兩顆觀測(cè)均有可能帶有一定的系統(tǒng)誤差，對(duì)系統(tǒng)誤差星的ddd如果還有時(shí)間和，還可考慮下列09號(hào)觀測(cè)單星觀測(cè)的01號(hào)空間飛行器進(jìn)行軌道估計(jì)，結(jié)二、問(wèn)題分對(duì)于第一個(gè)問(wèn)題，由于(2)式描述了的運(yùn)動(dòng)，因此以satinfo.txt文件中所提供的位置坐標(biāo)及速度作為初始條件，解算微分方程(2)即可得到軌道，進(jìn)而求得在特定時(shí)刻的位置。求出由測(cè)量數(shù)據(jù)得出的飛行器飛行軌跡。隨后選取參數(shù)vr(t)、m(t)的模型為 r(t)、m(t)mmt，其中v、m為常數(shù)，并利用最小二乘法vr vr 思想將該估值問(wèn)題轉(zhuǎn)化為參數(shù)優(yōu)化問(wèn)題，并利用遺傳算法確定vr(t、m(t模型標(biāo)飛行軌道的方法。本文將三軸指向誤差d、d、d視 坐標(biāo)系繞系繞三個(gè)坐標(biāo)軸的旋轉(zhuǎn)。該偏差由常值的系統(tǒng)誤差d、d、d三、模型假2、測(cè)量的白噪聲與輸入不相關(guān)3、忽略目標(biāo)紅外信號(hào)的時(shí)間及對(duì)紅外信號(hào)的處理時(shí)間4、兩顆測(cè)量的時(shí)鐘始終保持同步四、符號(hào)說(shuō)明及名Rc——目標(biāo)在地心慣性坐標(biāo)系下的位置矢量Rs——目標(biāo)在觀測(cè)坐標(biāo)系下的位置矢量sCc——測(cè)量坐標(biāo)系到地心坐標(biāo)系的轉(zhuǎn)換矩sh0、hf——初始（及終點(diǎn)）——誤差——?dú)埐钗?、模型的建立及?/p>

3r |r(t)

地心慣性坐標(biāo)系中，06號(hào)的初始條件為x0,6 .220573，y0,6 .771852，z0,6 vx0,64453.807606vy0,61566.513180vz0,64453.80760609號(hào)的初始條件為x0,9 .166765，y0,9 .631471，z0,9 vx0,95379.544693vy0,9407.095342vz0,93516.052656該微分方程，得到09號(hào)的軌道如圖3所示：圖 09號(hào)軌09號(hào)在下列時(shí)刻在地心慣性坐標(biāo)系中的位置如下xyzMM圖 逐點(diǎn)交匯法示意 1c R1S1

s111

1R2S2nCcs222

其中R1、R2為目標(biāo)在地心慣性坐標(biāo)系下的位置矢量，S1、S2為在心慣性坐標(biāo)系下的位置矢量，、、、為S1、S2的測(cè)量數(shù)據(jù) 為S1、S2測(cè)量坐標(biāo)系到地心慣性坐標(biāo)系的轉(zhuǎn)換矩陣c s (L)R(B)cosc s (L)R(B)sinLcos cos sinLsin

sis 0 s 其中L、B分別 當(dāng)前的緯度和經(jīng)度。注意到Cc1

Cc1 s11

s22 1 2目標(biāo)與S1、S2的連線在地心慣性坐標(biāo)系下的方向向量，分別記為I1、I2。nI2mI1S1S

設(shè)I 、

，S1S2

1

l2 zl1 2 j2nj1mknkm lnlm (k2j2)n(k1j1)my(lk)n(lk)mz

mx(k2l2)y(l2j2)z(j2k2 (k1l1)(j2k2)(j1 l j k

n (j2k2)(k1l1)(k2l2)(j1k1 xc1 xc2將(9)式代入(5)、(6)兩式，可解出R1y、R2 c1 c2zc1 zc2R1R2xc1xc2 R

y R 2 c2 z c2 將兩 在同一時(shí)刻的觀測(cè)數(shù)據(jù)作為1、1、2、2代入上述方法中決這一問(wèn)題。隨著ti從t0增加到tF，由此形成目標(biāo)飛行器的運(yùn)動(dòng)軌跡，如圖5圖 目標(biāo)飛行器運(yùn)動(dòng)軌 vr v立m(t)的模型為：m(t)m0 其中m為一常數(shù)且m0 r(t) 3r(t)vr m|r(t)

，

m0

r(t)K1

|r(t)

r(t)K2給定任意K、K及初始條件 、r(t)，則微分方程(12)可解，這里 r(t0 r(t0初始條件 r(t0r(r(t0nT0)R(t0nT0N2N其中TTtFt0 z該目標(biāo)量描述了由動(dòng)力學(xué)模型得出的目標(biāo)運(yùn)動(dòng)軌跡r(tym(tz

xc(t) cz

(t)至此，對(duì)K1、K2的估值問(wèn)題轉(zhuǎn)化為一個(gè)優(yōu)化問(wèn)題：即求K1、K2、r(t0的取值，使得M最小。本文采用遺傳算法實(shí)現(xiàn)這一優(yōu)化過(guò)程，遺傳代數(shù)為100K1 K K運(yùn)動(dòng)軌跡r(t。圖 目標(biāo)飛行器t-x曲圖 目標(biāo)飛行器t-y曲圖8 目標(biāo)飛行器t-z曲線圖 目標(biāo)飛行器t-vx曲圖 目標(biāo)飛行器t-vy曲圖 目標(biāo)飛行器t-vz曲txyz--------------------------------圖12為測(cè)量值與彈道估計(jì)值的偏差示意圖，圖中rs表示某一時(shí)刻目標(biāo)相對(duì)于的實(shí)際方位，Rs表示該時(shí)刻所測(cè)得的目標(biāo)的方位。S圖 測(cè)量偏差示意rsRs經(jīng)旋轉(zhuǎn)所得，由于rsRsRs到rs的旋轉(zhuǎn) Rsrs RsrsP arccosRsrs RsrsP在觀測(cè)坐標(biāo)系x、y、z軸上的投影分別為PsxPsyPszPPsxPsyPsx、PsyPszPsxPsy、Psz可表示測(cè)量坐標(biāo)系繞其x、y、z軸的歐拉角轉(zhuǎn)動(dòng)。如圖13x1y1z1繞z軸旋轉(zhuǎn)一個(gè)小角度x2y2z21 1在坐標(biāo)系x2y2z2中表示的M點(diǎn)的值 x2y2y1x x

y1x

d21My2M圖13 d與坐標(biāo)軸轉(zhuǎn)角的示意圖Psx

d

dd W1、W2、W3為零均值白噪聲過(guò)程在該時(shí)刻的樣本，d、d、d為常 按其采樣間隔不斷對(duì)目標(biāo)方位進(jìn)行探測(cè)，因此(17)式中、、dE(d)E(Psx

dd() E(Pszd06-1.69793×10-3.79254×10--3.52652×10-092.31695×10-1.27422×10--3.08118×10-DD

R(d)R(d)R(d

y z

D

D D 圖 目標(biāo)飛行器t-x曲圖 目標(biāo)飛行器t-y曲圖16 目標(biāo)飛行器t-z曲線圖 目標(biāo)飛行器t-vx曲圖 目標(biāo)飛行器t-vy曲圖 目標(biāo)飛行器t-vz曲txyz--------------------------飛行器在地心慣性坐標(biāo)系下的質(zhì)心運(yùn)動(dòng)微分方程如(1)所示，在飛行器的重力轉(zhuǎn)彎段，其攻角數(shù)值很小，為簡(jiǎn)化模型可視為零。即=0。攻角確定之后，飛行器的推力加速度模型FT成為彈道的因素。本文采用如下模型表示 Fv(t) =1 k-2kk

表示與速度r (21)所示；mm

k r(t)FeFT

r(t) rr|r(t) k2-飛行器在地心坐標(biāo)系下的位置速度參數(shù)(xc0yc0zc0vxc0vyc0vzc0，另一部是飛行器推力模型參數(shù)(k1,k2)對(duì)于積分初值(xc0yc0zc0vxc0vyc0vzc0而言，若直接將其作為待估參數(shù)低。本節(jié)通過(guò)對(duì)參數(shù)的等價(jià)變換，將(xc0yc0zc0vxc0vyc0vzc0的估計(jì)轉(zhuǎn)化為對(duì)(0,0,h0v00,hf)的估計(jì)，利用飛行器軌道的特殊性，獲得更加可靠的約(0,0)為初始觀測(cè)時(shí)刻 對(duì)飛行器的真實(shí)觀測(cè)量。 為起始觀測(cè)時(shí)飛行器的真實(shí)飛行高度，這三個(gè)參數(shù)用來(lái)確定位置初值(xc0yc0zc0，其等價(jià)轉(zhuǎn)xc0 1yc0 =xCs ss0 s0c0 x2+y2x2+y2 即(xc0,yc0,zc0)=f(0,0,h0 vx2+vy2 其中：vx2+vy2 度矢量方向落在由目標(biāo)初始位置矢量R0和目標(biāo)終點(diǎn)位置矢量Rf構(gòu)成的平面內(nèi)如圖20所示。即滿足： hMfR M020Rf可根據(jù)式(23)由(f,f,hf三個(gè)速度初值(vxc0,vyc0,vzc0)的等價(jià)轉(zhuǎn)換關(guān)系如下式所示vx2+vy2vx2+vy2 vxc0(RR)V=(RR)*vy c0

顯然，首發(fā)高程 可以由(26)式確定其區(qū)間范圍 x2+y2+z2-Re 其中Re為地球半徑，hminhmax分別為首發(fā)高程的最大值和最(0,0)（記為U0）約束區(qū)間確定。位置積分初值(xc0yc0zc0（記為Rc）與測(cè)量坐標(biāo)系下的坐標(biāo)(xs0,ys0,zs0)（記為Rs）關(guān)系如(27)所示： Rs=Cs(R-S+ 所處的地心坐標(biāo)系位置，表示誤差項(xiàng)，此處只c慮觀測(cè)的系統(tǒng)誤差。Cs表示從地心坐標(biāo)系到觀測(cè)坐標(biāo)系的坐標(biāo)轉(zhuǎn)化矩陣。則有觀測(cè)量的真實(shí)值和測(cè)量值滿足關(guān)系式[3]：c y y xs

y

ysxs

maxxs-

Hmin-其中Emax表示系統(tǒng)誤差上界，參考問(wèn)題三結(jié)果即可得到。Mmax表示目標(biāo)在Mmax=max{,}，根 提供的測(cè)量數(shù)據(jù)即可得到該上界。Hmin表示目標(biāo)同理可和測(cè)量滿足關(guān)系式由此可得(0,0)的約束式

Emax+MmaxEmaxHmin-Emax

-Emax+MmaxEmaxEmax+MmaxEmax Emax 對(duì)于(v0,0)的