2012年全國數(shù)學建模競賽優(yōu)秀選-b葛范湯

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

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

|r(t)

r

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

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

n rT，aa

roan

x

zaP在觀測坐標系中的坐標為SP

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

x

xx

2.52.6所示。x11x

6-X(m)-

6

x2.6上述求出的目標飛行器的運動軌跡是根據(jù)06號和09號兩顆的觀測數(shù)x6xx的運動規(guī)律。根據(jù)題目所給信息，在重力斜飛段，目標飛行器僅受重力G和推力T的作用，并且推力T的方向始終與目標飛行器的速度方向相一致。根據(jù)所上述性質(zhì)，首先在中擬合出一個平x6xx-1 -1.162.7接下來根據(jù)投影到平面上的目標飛行器的軌跡點擬合出三個位置分量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三中的運動軌n1n1n?21根據(jù)如下公式計算估計殘差：Qexn1n?n1n?211n nz-i21

，Qey Qez 。其中，ii 中編寫程序，求出三個方向上的估計殘差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é)所擬合出的位中分別繪制三個方程求出的kvf(t2.10③②③②①2.10三個方程分別求出的kvf(t接下來要求與上述方程①、②、③吻合得最好的一條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 中編寫程序，計算50s到70s每隔0.1s的kvf值，并繪制出kvf(t2.11①②③①②③最終曲2.11kvf(t然后對上一步求出的kvf(tkvf(t)=-4.93504e-18t8+4.92301e-15t7-2.08501e-12t6+4.67116e--1.68335e-08t4-8.15737e-06t3+0.00193023t2-

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

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

xx估計曲

xy估計曲

xz估計曲

Vx估計曲

Vy估計曲

Vz估計曲

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

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

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

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

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

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

-

--------

------

ay擬合曲 54

xxx030 3.56x-1x

6-1--

6x圖3.6三中的運動軌03.33.30-------------n1n1n1nn1n1n1n1n1n?21x,yn1n?21x,y,z 中編寫程序，求出三個方向上的估計殘差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)的解析表達式可以看出，kvf(t)m(t)=Cmexp-2.03895e-06t4+6.4385e-04t3-0.0125865t2其中Cm為一常單觀測問題空間飛行物相對地心坐標系的軌跡rt推導(dǎo)出的三個mt曲線應(yīng)該重合，如圖2.10所示。4.1所示，飛行器的合力