2012年全國(guó)數(shù)學(xué)建模競(jìng)賽優(yōu)秀選-b葛范湯

基于無(wú)源探測(cè)的空間飛行器一、問題的提有些國(guó)家會(huì)發(fā)射特殊目的的空間飛行器，如彈道式、偵察等。對(duì)他國(guó)發(fā)射具有敵意的空間飛行器實(shí)施并作出快速反應(yīng)，對(duì)于國(guó)家安全具有重要的意義。發(fā)現(xiàn)發(fā)射和探測(cè)其軌道參數(shù)是實(shí)現(xiàn)和作出反應(yīng)的第一步，沒有觀測(cè)，后續(xù)的判斷與反應(yīng)都無(wú)從談起。居高臨下，是探測(cè)過(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è)模型為本題以中低軌近圓軌道為觀測(cè)星座對(duì)假想的空間飛行器進(jìn)行仿真觀測(cè)，生成仿真觀測(cè)數(shù)據(jù)，要求利用仿真觀測(cè)數(shù)據(jù)，對(duì)假想空間飛行器的軌道參的281、根據(jù)觀測(cè)信息和觀測(cè)的簡(jiǎn)化運(yùn)動(dòng)方程，計(jì)算09號(hào)觀測(cè)成了交疊觀測(cè)，請(qǐng)結(jié)合幾何知識(shí)按照逐點(diǎn)交匯定位的思路，給出0號(hào)空間飛行器在公式（1）框架下的軌道估計(jì)，注意選取適當(dāng)?shù)膙r(t)和m(t)的表同時(shí)繪0號(hào)空間飛行器的t-x、t-y、t-z和三個(gè)速t-vx、t-vy、t-vz3若06和09測(cè)帶系對(duì)誤差進(jìn)行正確的估計(jì)能夠有效提高精度。利用上述的逐點(diǎn)交匯方法能否同時(shí)對(duì)系統(tǒng)誤差進(jìn)行估計(jì)？若不能，是否還有其他的思路能夠同時(shí)估計(jì)系統(tǒng)誤差與軌道？給出你的解決方案與估計(jì)結(jié)果。在報(bào)告中除給出與第二問要求相同的結(jié)果外，還應(yīng)分別給出兩顆觀測(cè)的系統(tǒng)誤差估計(jì)結(jié)果，共六個(gè)數(shù)值，分別是兩顆的d,d,d。觀測(cè)為中低軌近圓軌道r(t)FF r(t)v(t) |r(t) r(t)Fer(t|r(r(t)Fer(t|r(t)基礎(chǔ)坐標(biāo)系OcXcYcZc基礎(chǔ)坐標(biāo)系OcXcYcZc系和隨運(yùn)動(dòng)的觀測(cè)坐標(biāo)系OsXsYsZs如下圖所Fe：表示飛行器所受的外力加速度之和r(t)：為空間飛行器在基礎(chǔ)坐標(biāo)系下的位置矢量r(t)：表示r(t)對(duì)時(shí)間tmGm：為地球引力常數(shù)（本題中地球引力常數(shù)取G3.9860051014m3/s2mvr(t)：是相對(duì)于火箭尾部噴口的噴射速度ys

,三軸指向誤差用三個(gè)常值小量dddd：第一觀測(cè)量d：觀測(cè)量在觀測(cè)軌道的計(jì)算r(t)Fe

|r(t)

r

值解方法計(jì)算空間飛行器的運(yùn)動(dòng)軌跡。不同空間飛行器的本質(zhì)差異就在于vr(t)和m(t)的模型不同，m(t)一般而言應(yīng)為嚴(yán)格單調(diào)遞減的非負(fù)函數(shù)。vr(t的方向傳統(tǒng)的常微分方程組的數(shù)值方法是Runge-Kutta法。但由于四階Runge-Kutta法不收斂；本文中采用了Runge-Kutta-Fehlberg法，通過(guò)控制本文在中編寫程序，計(jì)算出09號(hào)的運(yùn)行軌道，如圖1.1所250.0s五個(gè)時(shí)刻的三維位置，如表1.1所示，結(jié)果保留6位有效數(shù)字。xx642y0--5-0xx--42z0z---8x圖1.109號(hào)的運(yùn)行軌表1.109號(hào)觀測(cè)在五個(gè)離散時(shí)刻的三維位時(shí)刻利星定位估計(jì)空間飛行器的軌道雙星觀測(cè)數(shù)據(jù)時(shí)間不對(duì)齊問題的解決通過(guò)觀察仿真數(shù)據(jù)文件meadata_06_00.txt和meadata_09_00.txt中的仿真觀測(cè)數(shù)據(jù)信息可以發(fā)現(xiàn)雖然06號(hào)和09號(hào)均對(duì)00號(hào)空間飛行器進(jìn)行觀測(cè)，但是它們的觀測(cè)時(shí)間并不嚴(yán)格對(duì)齊。因此，首先需要對(duì)其中一顆的觀測(cè)數(shù)據(jù)進(jìn)行擬合，求出其在另一顆觀測(cè)時(shí)刻上的觀測(cè)值。二乘法對(duì)觀測(cè)數(shù)據(jù)進(jìn)行多項(xiàng)式擬合。本節(jié)以09號(hào)的觀測(cè)時(shí)刻為準(zhǔn)，對(duì)06時(shí)刻上。在中編寫程序，繪制出擬合曲線以及原始數(shù)據(jù)點(diǎn)如圖2.1所原原始觀測(cè)數(shù)據(jù) （a）（a）觀測(cè)量 擬（b）觀測(cè)量圖2.106號(hào)觀測(cè)數(shù)據(jù)的擬合同時(shí)得出觀測(cè)量、的擬合曲線表達(dá)式為：-4.17448e-16t61.82007e-13t5-3.87875e-，。-8.56294e-09t3-1.15533e-06t2然后將09號(hào)的觀測(cè)時(shí)刻代入上述多項(xiàng)式表達(dá)式求出06號(hào)在09號(hào)觀測(cè)時(shí)刻上的觀測(cè)值，在中編寫程序，繪制出擬合曲線以及對(duì)齊后的始數(shù)據(jù)點(diǎn)如圖2.2所示。對(duì)對(duì) 擬合曲擬合曲對(duì)齊后測(cè)數(shù)據(jù)（a）0000 （b）觀測(cè)量對(duì)齊后的數(shù)據(jù)點(diǎn)圖2.206號(hào)對(duì)齊后的觀測(cè)數(shù)據(jù)計(jì)算雙星在觀測(cè)時(shí)刻上的位置1Runge-Kutta-Fehlberg方法，在06號(hào)和09號(hào)的運(yùn)動(dòng)方程，并繪制出軌道以及兩在09號(hào)的觀測(cè)時(shí)刻上的位置，如圖2.3所示，圖中粗線表示06號(hào)和09號(hào)兩在09號(hào)xx09星軌06號(hào)軌864206240420x8x圖2.306號(hào)和09號(hào)的運(yùn)行軌基礎(chǔ)坐標(biāo)系與觀測(cè)坐標(biāo)系之間的轉(zhuǎn)換基礎(chǔ)坐標(biāo)系OcXcYcZc系和隨運(yùn)動(dòng)的觀測(cè)坐標(biāo)系OsXsYsZs如下圖所2.4設(shè)觀測(cè)坐標(biāo)系原點(diǎn)O在基礎(chǔ)坐標(biāo)系中的坐標(biāo)為 zT。則坐標(biāo) Xs在基礎(chǔ)坐標(biāo)系X,Y,Z三軸上的投影分別為xosrn，yosrn和zosrn，其 prn x2y2z2；坐標(biāo)軸Zs在基礎(chǔ)坐標(biāo)系X,Y,Z p

zr2zoszossgnzosa2a2 ra 。注意，由于觀測(cè)坐標(biāo)系是，UEN坐標(biāo)系，a2a2 設(shè)觀測(cè)坐標(biāo)系相對(duì)于基礎(chǔ)坐標(biāo)系的旋轉(zhuǎn)矩陣為CRn

n rT，aa

roan

x

zaP在觀測(cè)坐標(biāo)系中的坐標(biāo)為SP

S該點(diǎn)在基礎(chǔ)坐標(biāo)系中的坐標(biāo)應(yīng)為CPCTSPS同理，設(shè)第i號(hào) 的觀測(cè)數(shù)據(jù)為 i，即從06號(hào) 坐標(biāo)系中的坐標(biāo)為1 ，那么該方向向量在世界坐 iTiCT1 T 觀測(cè)誤差的存在，在同一觀測(cè)時(shí)刻06號(hào)的觀測(cè)射線與09號(hào)的觀測(cè)射P1P2Pmid2.5所示。因此，Pmid該觀測(cè)時(shí)

x

xx

2.52.6所示。x11x

6-X(m)-

6

x2.6上述求出的目標(biāo)飛行器的運(yùn)動(dòng)軌跡是根據(jù)06號(hào)和09號(hào)兩顆的觀測(cè)數(shù)x6xx的運(yùn)動(dòng)規(guī)律。根據(jù)題目所給信息，在重力斜飛段，目標(biāo)飛行器僅受重力G和推力T的作用，并且推力T的方向始終與目標(biāo)飛行器的速度方向相一致。根據(jù)所上述性質(zhì)，首先在中擬合出一個(gè)平x6xx-1 -1.162.7接下來(lái)根據(jù)投影到平面上的目標(biāo)飛行器的軌跡點(diǎn)擬合出三個(gè)位置分量p-7.18812e-10t62.44420e-07t55.13831e-155.495t155.495t p2.24400e-10t61.54686e-07t57.03955e- -0.0217185t35.42187t2317.048tp8.70490e-10t62.93240e-07t56.03146e- +0.0150115t36.13428t2178.866t并且求出位置量的一階導(dǎo)數(shù)

dp,

d

,vd

和二階導(dǎo)數(shù)

adv,adv,adv。 dt dt dt

x

66666666

x

x

vx擬合

vy擬合曲 vz擬合曲

0

ax擬合曲

ay擬合曲6543

az擬合曲

2.8xx

x圖2.9三中的運(yùn)動(dòng)軌n1n1n?21根據(jù)如下公式計(jì)算估計(jì)殘差：Qexn1n?n1n?211n nz-i21

，Qey Qez 。其中，ii 中編寫程序，求出三個(gè)方向上的估計(jì)殘差Qex=Qey=8.58835(m)Qez5.48184(m)vr(t)FF r(t) mv

|r(t) 因此，可設(shè)v(tkv(t)k<0m(t)f(t v vv(t)r(t)|r(t)r(t)kf vv(t)r(t)rx|r(t)r(t) ①v(t)ry(t)G|r(t)②v(t) G|r(t) kf(t)kvf(t)kf(t)可以看出，方程①、②、③均可以求出kvf(t2.5節(jié)所擬合出的位中分別繪制三個(gè)方程求出的kvf(t2.10③②③②①2.10三個(gè)方程分別求出的kvf(t接下來(lái)要求與上述方程①、②、③吻合得最好的一條kvf(t變化曲線?？蓃(t) r(t) |r(t)|3 rry(t)kvf(t)=v(t)ry(t)|r(t)|3ry(r zr(t) r(t)z |r(t)|3 r(t)r(t)x|r(t)r(t)xrx A=r(t)B=v r

r(t) r(t) |r(t)|3 vv 中編寫程序，計(jì)算50s到70s每隔0.1s的kvf值，并繪制出kvf(t2.11①②③①②③最終曲2.11kvf(t然后對(duì)上一步求出的kvf(tkvf(t)=-4.93504e-18t8+4.92301e-15t7-2.08501e-12t6+4.67116e--1.68335e-08t4-8.15737e-06t3+0.00193023t2-

從kvf(t)的解析表達(dá)式可以看出，kvf(t)m(t)=Cmexp-2.03934e-06t4+6.43412e-04t3-0.0125231t2其中Cm為一常擬合擬合線2.12kvf(t因此，將kvf(t)的表達(dá)式代入方程（2.1可得空間飛行器的運(yùn)動(dòng)學(xué)方rr(t)FF |r(t)r(t)kf(t)rvr中編寫程序，計(jì)算出微分方程（2.7）50s170s區(qū)間內(nèi)的數(shù)值解，并繪制2.1350.0s170.0s10.0s進(jìn)行采樣，計(jì)同時(shí)0號(hào)空間飛行器t-x、t-y、t-z和三個(gè)速t-vx、t-vy、t-vz2.146x6-xX(m)-

6x2.132.10V（m/s）-------------

xx估計(jì)曲

xy估計(jì)曲

xz估計(jì)曲

Vx估計(jì)曲

Vy估計(jì)曲

Vz估計(jì)曲

2.14估計(jì)對(duì)空間飛行器觀測(cè)的系統(tǒng)誤差用三個(gè)常值小量dd,d來(lái)表示，分別表示第一觀測(cè)量的平移量、第二觀測(cè)量的平移量以及觀測(cè)量在平面內(nèi)的旋轉(zhuǎn)量。的x(t),y(t),z(t)滿足目標(biāo)飛行器的簡(jiǎn)化運(yùn)動(dòng)方程。x(ty(tz(t進(jìn)行反推計(jì)算，與觀測(cè)的進(jìn)行比較，06號(hào)和 的觀測(cè)數(shù)據(jù)與估計(jì)數(shù)據(jù)分別如表3.1和表3.2所示，其中=估計(jì)-測(cè)量 ，均為微小量，且隨時(shí)間變化不明顯，因此這一偏差主要不是由于隨機(jī)誤差引表3.106號(hào)觀測(cè)值與估計(jì)值比時(shí)刻測(cè)測(cè)估估表3.209號(hào)觀測(cè)值與估計(jì)值比時(shí)刻測(cè)量測(cè)量計(jì)算計(jì)算------------------------------------------------接由觀測(cè)得出的觀測(cè)量。在-平面中，點(diǎn)，繞坐標(biāo)原點(diǎn)旋轉(zhuǎn)角度d，得到點(diǎn)的坐標(biāo)為cosdsind，cosdsind，再將該點(diǎn)平移向量d，d，得到點(diǎn)的坐標(biāo)為cosdsind+d，cosd+sind+d, 3.1由于系統(tǒng)誤差全部折算到三軸指向誤差ddd?，?=cosd-sindd，cosdsindd由于為一個(gè)微小的旋轉(zhuǎn)量，所以有cosd=1,sindd。將該關(guān)系代入上

由公式(3.1)、(3.2)可知，?-與存性關(guān)系，?-與也 性據(jù)，得到與?-之間的關(guān)系為：?-0.00014050.0002417，從而可以得出d0.0001405，d0.00024173.2所示，采用一次多項(xiàng)式擬合數(shù)據(jù)，得到與?-?-0.00020620.0003062，從而可以得出d0.0002062d0.0003062

圖3.206 與?-之間的關(guān)

圖3.306 與?-之間的關(guān)將（3.2） ，從而可得d000017dd6.542e05綜上，06號(hào)的三軸指向誤差為d0.0001627，d0.0002417d0.0003062然后對(duì)09號(hào)的數(shù)據(jù)進(jìn)行擬合。如圖3.4所示，采用一次多項(xiàng)式擬合數(shù)出d0.0003108，d0.00027133.5所示，采用一次多項(xiàng)式擬合數(shù)據(jù)，得到?-?-0.00036100.0005022，從而可以得出d0.0003610d0.0005022。

圖3.209 與?-之間的關(guān)

圖3.309 與?-之間的關(guān)綜上，09號(hào)的三軸指向誤差為d0.0003359，d0.0002713d0.0005022考慮了系統(tǒng)誤差，對(duì)測(cè)量值進(jìn)行系統(tǒng)誤差修正，利用修正后的結(jié)果重新估計(jì)空間飛行器的運(yùn)動(dòng)學(xué)方程，在中分別繪制三個(gè)方程求出的kvf(t)變化情況，如圖3.4所示。3.4中可以看出：三條曲線幾乎重合在一起，這說(shuō)明本節(jié)求出的擬合3.4系統(tǒng)誤差修正后的kvf(t考慮了系統(tǒng)誤差，對(duì)測(cè)量值進(jìn)行系統(tǒng)誤差修正，利用修正后的結(jié)果重3.5所示。xx---1--1

-

--------

------

ay擬合曲 54

xxx030 3.56x-1x

6-1--

6x圖3.6三中的運(yùn)動(dòng)軌03.33.30-------------n1n1n1nn1n1n1n1n1n?21x,yn1n?21x,y,z 中編寫程序，求出三個(gè)方向上的估計(jì)殘差Qex=Qey=8.58727(m)Qez5.48488(m)3.3節(jié)求出的kvf(tkvf(t)=-4.95254e-18t8+4.9396e-15t7-2.0916e-12t6+4.6852e-從kvf(t)的解析表達(dá)式可以看出，kvf(t)m(t)=Cmexp-2.03895e-06t4+6.4385e-04t3-0.0125865t2其中Cm為一常單觀測(cè)問題空間飛行物相對(duì)地心坐標(biāo)系的軌跡rt推導(dǎo)出的三個(gè)mt曲線應(yīng)該重合，如圖2.10所示。4.1所示，飛行器的合力