1、參賽（由填寫）第九屆“杯”研究生數學建模競賽學校工程大學參賽隊號900520161.李思佳隊員2.鹿傳國3.劉玉磊參賽（由填寫）第九屆“杯”研究生數學建模競賽基于無源探測的空間飛行器主動段軌道估計與誤差分析題 目摘要：本文主要研究了無源探測空間飛行器主動段的軌道估計及誤差分析問題，根據 4 個問題的不同要求，建立了相應的數學模型并求解，實現了軌道估計、飛行器軌道估計、系統誤差估計以及軌道估計和誤差分析等。主要工作包括：（I）問題 1：對描述觀測軌道的簡化運動方程這個二階微分方程進行了降階，繼而提出了一種自適應定點的四階龍格庫塔算法求解運動微分方程，估計了 09 號軌跡，該星在 50s、100s

2、、150s、200s、250s 時刻的位置分別為（1.77381´106 ， 8.16138´106 ， 4.51669 ´106 ）、（ 1.50163´106 ， 8.12676 ´106 ， 4.68468´106 ）、（ 1.22770 ´106 ， 8.08270 ´106 ， 4.84722´106 ）、（ 9.52349 ´105 ， 8.02925´106 ，5.00413´106 ）、（ 6.75894 ´105 ， 7.96650 ´1

3、06 ， 5.15524 ´106 ）。（II）問題 2：合的數據同步算法。因觀測的數據時間異步問題，提出了一種基于樣條插值擬對飛行器的量測所在的觀測坐標系與數據處理所在的基礎坐標系不一致，推導了觀測坐標系與基礎坐標系的轉換關系，將數據信息統一到基礎坐標系下。根據空間飛行器在觀測坐標系下的約束關系，建立了雙星無源逐點交匯的算法模型，并利用同步后的數據解算出雙星觀測條件下的飛行器空G度。根據vr (t) 和m(t) 間軌跡，繼而基于中心差分理論得到飛行器各時刻的速度的物理意義和變化規律，對二者分別建立了含有未知參量的近似模型，完善了主動段空間飛行器的運動方程，直接將飛行器的軌跡、速度和

4、度數據代入運動方程可實現模型中未知參數的估計，為簡化參數估計，提出了一種基于級數理論的多項式模型來擬合vr (t)m(t) / m(t) 這一整體，給出了對多項式各階次系數的估計方法，并對不同階次多項式模型的性能進行了分析對比。空間飛行器的運動數據同（I）用四G 階龍格庫塔法解算運動方程得出。基于上述模型解算出了 0 號飛行器在各個采樣點的位置和速度，各位置分量的殘差分別為：1.60681´102 、1.19288´102 、1.97124 ´102 ，各速度分量的殘差分別為： 6.78000 、4.03933 、8.33612 。（III）問題 3：根據坐標轉換

5、模型和觀測坐標系下的量測方程，確定系統誤差、飛行器空間位置和量測數據的關系，建立方程組，通過對方程組未知量個數與自由度關系的分析，發現不能利用逐點交匯方法估計系統誤差；進而提出了一種系統誤差估計的滑窗四點交匯方法，該方法通過相鄰四點的滑動，實現了系統誤差動態估計。估計出的系統誤差為 da 06 = 1.35166 ´10， db 06 = 1.30491´10，-3-3dq 06 = 8.96365´10 ，da= 8.38309´10 ，d= 5.15881´10 ，d= 1.17204 ´10 。-4-4-4-309b 09q 0

6、9通過求解的系統誤差對附件的觀測數據進行修正，利用修正后的數據，采用（II）的方法對飛行器進行軌道估計和誤差分析。（IV）問題 4： 運動特征，建立了單星單飛行器無源軌道估計問題，結合短時間內飛行器的序貫多點軌道估計模型、無約束最優化軌道估計模型和基于軌道平面約束的軌道估計模型。利用（III）求得的系統誤差對觀測數據進行修正，繼而基于序貫多點軌道估計模型解算出飛行器的軌道參數，實現了單對 1 號飛行器的軌道估計和誤差分析。根據坐標轉換模型、觀測坐標系下的量測方程，通過討論關于軌道參數和系統誤差值的方程組是否可解來確定多星多飛行器聯合軌道與誤差估計的可行性。通過比較方程組自由度和待估計的軌道參數

7、和系統誤差的個數，得出了可估計系統誤差的條件：在觀測滿足2MP ³ 3P + 2M 時，系統誤差可估計。個數 M 與空間飛行器數目 P：軌道估計，系統誤差，四階龍格庫塔算法，滑窗四點交匯，基于級數理論的多項式模型，序貫多點目錄一、問題重述2二、問題假設3三、符號說明4四、問題分析5五、模型的建立與求解65.1問題一的建模與求解75.1.1軌道模型及算法75.1.2計算結果及分析9問題二的建模與求解105.25.2.15.2.25.2.35.2.4觀測數據同步11雙星逐點交匯的建模與求解16G (t) 和m(t) 的模型建立與參數估計21vr飛行器軌道估計及殘差計算265.3問題三的建

8、模與求解315.3.15.3.25.3.3逐點交匯法系統誤差估計可行性分析31系統誤差求解32飛行器軌道估計與殘差計算335.4問題四的建模與求解385.4.15.4.2單星單飛行器無源觀測38多目標系統誤差估計46多六、模型優缺點分析48一、問題重述對他國發射具有敵意的空間飛行器實施并做出快速反應，對于維護安全具有重要的戰略意義。發現發射和探測飛行器軌道參數是實現和做出反應的第一步， 探測空間飛行器發射沒有觀測，后續的與反應都無從談起。居高臨下，是與軌道參數的重要平臺。裝置于上的探測器包括有源和無源兩類：有源探測器采用主動方式（如，激光）搜尋目標，同時具備定向和測距兩種能力；無源探測器則接收

9、目標輻射。采用無源探測器的觀測常采用紅外光學探測器，只接收目標的紅外輻射信息，可定向但不能測距。對于火箭尾部噴焰的高度敏感性是紅外技術的長處，但易受氣候影響與云層干擾則是其缺點。探測的目的是為了推斷空間飛行器的軌道參數，推斷是基于觀測數據并通過數學模型與計算方法做出的。當觀測飛行一段時間，探測器測得目標相對于運動的觀測數據，以觀測和空間飛行器的運動模型和觀測模型為基礎，對空間飛行器的軌道參數（包括軌道位置、速度初值和其他模型參數）進行數學推斷，為飛行器類別、飛行意圖的提供信息基礎。空間飛行器軌道一般可分為三段，依次為：靠火箭推進的主動段、在地球外層空間的慣性飛行段和再入大氣層后的段。主動段通常

10、由多級火箭相繼推進，前一級火箭完成推進后脫落，由后一級火箭接力。慣性飛行段在空氣阻力極小的大氣層外，靠末級火箭關機前獲得的速度在橢圓軌道上作無動力慣性飛行。段則根據任務需求，受后再入大氣層，飛向目標。對于段。而言，在其結束前一直繞地飛行，故無主動段又可細分為若干子段：垂直上升段，程序拐彎段和重力斜飛段。按最優軌道設計，為節約，箭體應盡快穿過稠密大氣層，故火箭一般先垂直發射。程序拐彎段連接垂直上升段與重力斜飛段，在外力矩下使箭體轉過一定角度，該段完成后外加力矩撤銷，進入斜飛狀態。第一級火箭通常負擔“垂直段+程序拐彎段（加外力矩）+重力斜飛段的前段”的推進（視發的特性），重力斜飛段的后程則靠第二、

11、第三級火箭相繼完成。由于斜飛狀態下地球引力與推力不在同一直線，所以箭體質心的運動軌跡為帶一定弧度的光滑曲線。單個紅外光學探測器不具備測距能力，但借助多顆（含兩顆）觀測的同步觀測能夠進行逐點，再結合空間飛行器的運動模型，可以進行軌道參數估計。在單星觀測條件下，利用空間飛行器軌道的特殊性，結合較強的模型約束也可得到一定精度軌道參數估計。由于受大氣影響，垂直上升段的火箭尾焰不易觀測，程序拐彎段的運動方程又較為復雜，所以本題重點關注重力斜飛段的后程段，本題的研究集中于此段。觀測數據不可避免地帶有各種誤差，觀測誤差包括隨機誤差和系統誤差。假設隨機誤差為直接疊加在觀測數據上的白噪聲，可能產生于背景輻射干擾

12、與信息處理等多個方面。經由適當的簡化模型，各種系統誤差最終可以折合為觀測坐標系的原點位置誤差和三軸指向誤差。根據工程經驗，原點位置誤差影響較小，而三軸指向誤差影響較大，對三軸指向誤差進行估計對于提高估計精度很有幫助，本題只考慮三軸指向誤差。三軸指向誤差在二維觀測數據平面上表現為兩個平移誤差和一個旋轉誤差。現以中低軌近圓軌道為觀測星座對假想的空間飛行器進行觀測，生成觀測數據，要求利用觀測數據，對假想空間飛行器的軌道參數進行估計。在僅考慮隨機誤差的條件下，研究下列問題：問題1：觀測在任意時刻的位置計算是估計的前提，請根據附件satinfo.txt和觀2測的簡化運動方程（2），計算09號觀測五個時刻

13、的三維位置。結果保留6位有效數字。在50.0s、100.0s、150.0s、200.0s、250.0s問題2：在本題給定的數據下，06號和09號觀測對0號空間飛行器形成了立體交疊觀測，請結合立體幾何知識按照逐點交匯的思路，給出0號空間飛行器在公式（1）框架下的軌道估計，注意選取適當的 G (t) 和m(t) 的表示模型。按照從50.0s到170.0svr間隔10.0s進行采樣，計算并列表給出0號空間飛行器在各個采樣點的位置和速度，并給出估計殘差。結果保留6位有效數字。同時繪制0號空間飛行器的三個位置t-x、t-y、t-z和三個速度t-vx、t-vy、t-vz曲線示意圖。問題3：若06和09號兩

14、顆觀測均有可能帶有一定的系統誤差，對系統誤差進行正確的估計能夠有效提高精度。利用上述的逐點交匯方法能否同時對系統誤差進行估計？若不能，是否還有其他的思路能夠同時估計系統誤差與軌道？給出你的解決方案與估計結果。在報告中除給出與第二問要求相同的結果外，還應分別給出兩顆觀測的系統誤差估計結果，共六個數值d a , d b , d q ，分別是兩顆的 。問題4：對只有09號觀測單星觀測的01號空間飛行器進行軌道估計，結果形式要求同第三問，注意參考第三問的系統誤差估計結果。并進一步考慮在同時有多顆觀測觀測多個空間飛行器的情況下能否進行系統誤差估計？二、問題假設根據題意，可以進行如下假設：1、地球公轉周期

15、遠大于空間飛行器的觀測弧段時長，故本題在短時間內認定基礎坐標系為慣性坐標系，該基礎坐標系不隨地球旋轉；2、目標運動具有連續性和漸進性；3、假設隨機誤差為直接疊加在觀測數據上的白噪聲，可能產生于背景輻射干擾與信息處理等多個方面；4、假設利用兩個平移誤差可以代表三軸誤差；5、m(t) 一般而言應為嚴格單調遞減的非負函數， Gvr (t) 的方向一般應與飛行器的速度方向相同，其大小穩定；6、題目僅關注飛行器的重力斜飛段的后程段；7、不考慮空間飛行器各點的速度差異；8、不考慮空間飛行器因火箭消耗而產生的位移；9、作用于飛行器上的所有外力都通過質心，即不計空間飛行器在飛行過程中繞質心的轉動。3三、符號說

16、明中涉及到的符號說明如下：為方便理解，其余符號將在文中陸續引出。4符號含義備注G地球引力常數G = 3.986005 *1014m3 / s2mmG (t)為空間飛行器在基礎坐標系下的位置矢量mrc G (t)空間飛行器在基礎坐標系下的度m s2rcm(t )瞬時質量kgm (t)質量變化率kg sa(t)t 時刻 m 號對某一飛行器的估計比值a = ymcxb(t)t 時刻 m 號對某一飛行器的估計比值b = zmcxrG (t)t 時刻 m 號位置r (t) 在 X O Y 面得投影向量mc c c G (t)空間飛行器在基礎坐標系下位置矢量在O X 方mrxss G (t)空間飛行器在基

17、礎坐標系下位置矢量在O Y 方向m rys s G (t)空間飛行器在基礎坐標系下位置矢量在O Z 方mrzs sGRi (t)i 號在基礎坐標系下的位置矢量mG (t)i 號位置矢量 G (t) 在O X 方向上的分量m xiRissG (t)i 號位置矢量 G (t) 在O Y 方向上的分量m yiRis sG (t)i 號位置矢量 G (t) 在O Z 方向上的分量m ziRis s G (t)i 號觀測的飛行器度m s2risG (t)i 號觀測的飛行器速度m srisGJJJJGeXs觀測坐標系的方向向量Oc XcGJJJGeXy觀測坐標系的方向向量OcYcGJJJJGeXy觀測坐標

18、系的方向向量Oc Zc四、問題分析現階段，太空戰略已成為各個競相發展的一個方面，而無源探測系統是太空探測系統的重要組成部分，它不僅關系到導航、氣象、測繪等民用的各個方面，還關系到了領空安全，更是一個國防能力的體現。與有源探測相比，無源探測具有探測距離遠、頻率分配受限少、抗電子干擾和抗毀能力強等優勢。在太空成為新戰場的今天，無源探測技術已成為取得太空優勢的之一。無源探測器就是不主動輻射電磁波，只是接收目標輻射。這使其可定向但不能測距，從而導致比較。目前無源探測是探測研究的熱點，難點問題之一。因為無源探測常采用紅外光學探測器，只接收目標的紅外輻射信息，對于火箭尾部噴焰的高度敏感性是紅外技術的長處，

19、但易受氣候影響與云層干擾則是其缺點，所以無源探測的主要任務之一就是對彈道等空間飛行器的軌道進行估計，推斷該類空間飛行器的軌道參數。當觀測飛行一段時間，探測器測得目標相對于運動的觀測數據，以觀測和空間飛行器的運動模型和觀測模型為基礎，對空間飛行器的軌道參數（包括軌道位置、速度初值和其他模型參數）進行數學推斷，為飛行器類別、飛行意圖的提供信息基礎。該類空間飛行器軌道一般包含三個階段：火箭推進的主動段、地球外層空間的慣性飛行段和再入大氣層后的段。主動段又可細分為若干子段：垂直上升段，程序拐彎段和重力斜飛段。由于重力斜飛段的特殊性，因此對其研究成為了基于空間飛行器軌道估計的重點。該類空間飛行器在基礎坐

20、標系下的主動段的簡化運動方程如下：無源探測的GG G GmGGm (t) r (t) = Fe + FT = -3 r (t) + vr (t) (4-1)G| r (t) |m(t)觀測的簡化運動方程為：G G GmGr(t) = Fe =-3 r(t) G| r(t) |(4-2)通過兩顆進行無源探測的空間飛行器主動段軌道估計過程如圖 4.1 所示。紅色陰影（模塊一）完成了的軌跡估計，綠色陰影（模塊二）表示了通過逐點交匯方法對空間飛行器，并對所建的Gv (t) 和m(t) 模型進行參數估計，得出飛行器的運動方r程進一步估計飛行器軌道，并估計殘差。陰影（模塊三）代表通過時間段交匯方法求解系統

21、誤差，然后在提高精度的前提下，對空間飛行器的軌道進行估計。5G (t) vrm(t) 圖 4.1基于無源探測的空間飛行器軌道估計問題一需要解決的問題是：根據題目給出的觀測簡化運動方程，結合初始條件的位置矢量和速度矢量，利用常微分方程組的數值方法計算出的軌跡，同時確定幾個整數時刻的位置矢量及速度矢量。需注意的一點是微分方程，還需要進行降階。的運動方程是二階常問題二需要解決的問題可以分為四個方面：一是觀測數據同步問題（對應圖4.1 中的雙星同步模塊），即對于不同表示同一目標位置的時刻集合擴大為一致，使得所有觀測軌跡的采樣數據同步可利用樣條插值的方法，將各觀測數據進行同步；二是通過坐標轉換和雙星無源

22、逐點交匯方法，計算得到雙星觀測條件下的飛行器實際空間軌跡；三是通過對所給Gv (t) 和m(t) 的建模，利星無源逐點交匯方法得到的軌跡r對模型中的參數進行估計；四是結合所得到的估計參數，利用飛行器運動方程，通過微分方程組的解出估計出空間飛行器的軌跡，比較通過雙星觀測的飛行器空間軌跡和通過微分方程組解出空間飛行器的軌跡，求解殘差，形成對比。問題三需要解決的問題是：確定雙星無源逐點交匯方法的系統誤差估計的可行性，使用滑窗四點交匯法估計出系統誤差，通過求出的系統誤差對附件給出的 06 與 09 號衛星觀測數據進行修正，利用修正后的數據，采用第二問的解法對飛行器進行逐點交匯定位和軌道估計，并計算出相

23、應的殘差。問題四需要解決的問題可以分為兩個方面：一是研究單星無源探測的特殊性，尋求建立綜合多時刻觀測集合的單星模型，通過利用第三問求出的系統誤差對 09號的觀測數據進行修正，利用修正后的觀測數據對 01 號飛行器進行和軌道估計。二是對多多飛行器無源跟蹤問題進行研究，將系統誤差估計轉化到方程是否可解的問題上來。通過對方程個數與未知量的個數的分析比較，得出可以估計系統誤差的條件。6五、模型的建立與求解5.1問題一的建模與求解軌道的精確估計是空間飛行器軌道估計的基礎，該問題的目的便是對的軌道參數進行估計。在已知初始時刻的軌道參數及其運動方程的條件下，通過解常微分方程即可得到軌道的參數。在對常微分方程

24、進行降階之后，采用了一種自適應定點的四階龍格庫塔算法，得到軌道的數值解，為空間飛行器的軌道估計奠定基礎。5.1.1軌道模型及算法根據題中所給的觀測程組：的簡化運動方程及其初始條件建立軌跡的常微分方2 G ìGGï d R (t) Gm= F= -09R (t) Ge09dt 23| r (t) |ïïG ï dR (t) = v(t) (5.1-1)í09dtï Gï R(0) = ( x)T, y , zï09000ïîv(0) = (vx (0), vy (0), vz (0)GG

25、()代入 R09 (t) = ( x(t), y(t), z(t)TT= v (t), v (t), v (t)于公式(5.1-1)得到模型：和v (t) xyzì d 2 x(t) Gdx(t) ï=-mx(t) = vx (t) ï dt 2ï3/ 2éë x2 (t) + y2 (t) + z2 (t)ùûGdtï d 2 y(t) dy(t) =- = vy (t) my(t) (5.1-2)íï3/ 2éë x2 (t) + y2 (t) + z2 (t)

26、ùûGdt 2dtïï d 2 z(t) =- dz(t) =mz(t) vz (t) ïïî23/ 2dtdtéë x2 (t) + y2 (t) + z2 (t)ùû式(5.1-2)了基本的運動軌道模型，進一步初始時刻的初值，采用自適應定點的四階龍格庫塔算法1(Adaptive Fixed Point-4-Optical Step Runge-Kutta, AFP-4-OSRK)求解出軌跡的數值解。表 5.109 號的初始時刻位置與速度7位置/m數值速度/m/s數值x(t ) 20

27、43922.166765vx (t) -5379.544693y(t ) 8186504.631471vx (t) -407.095342z(t ) 4343461.714791vx (t) -3516.052656由于軌道的確定是空間飛行器軌道估計的基礎，因此很有必要在軌道估計中保持較高的精度，同時確定采樣點時刻軌道估計值。為此在傳統定步長的基礎上提出了 AFP-4-OSRK 算法，該算法既保持了步長的自適應能力，同時提高了算法的精度。具體的 AFP-4-OSRK 算法如下：(t) = dx(t) v= dy(t) v= dz(t) Step1：根據動力學規律有v, y ( )t, z (

28、) t則xdtdtdtd (vy (t)d 2 x(t)d (v (t) d 2 y(t) d 2 z(t)d (v (t)=x,dt=zdt,dt 2dt 2dt 2dtd 2 x(t)dv (t) dx根據動力學規律有v(t) =，則dt為如下式的常系數微分方程組;=x。進而將模型(公式(5.1-2)降階變換dtdt 2æ dvx (t) öæöGç÷- mx(t) dtç÷÷ç÷é x2 (t) + y2 (t) + z2 (t)ù3/ 2ç dvy

29、 (t) ÷çëûç÷ç÷dtG÷ç -my(t) ÷çç÷÷÷ç dvz (t) ÷3/ 2éë x2 (t) + y2 (t) + z2 (t)ùûGm÷ = çG(t) = çdtç -(5.1-3)çdx(t) ÷z(t) ÷÷÷÷÷çç&#

30、231;ççç÷é x2 (t) + y2 (t) + z2 (t)ù3/ 2çdt÷ëûççççdy(t) ÷v (t) ÷xdt÷v (t) yvz (t) dz(t) ÷ç÷èøç÷èdtøStep2：在步長h 條件下建立如下的模型公式：ìG(t +1) = G(t) + h (K + 2K + 2K + K )ï

31、;ïï12346K = f (G(t)1ïK = f (G(t) + h K )(5.1-4)íïïïïïî212K = f (G(t) + h K )322K4 = f (G(t) + hK3 )Nt +1/ 2Step3：嘗試以h /2 為步長一次，結果記為 Y，再將步長連續計算三步i,h/ 2i=1Nt +1/ 2得到同一點的值記為 Y；i,h/6i=1Y t +1/ 2 - Y t +1/ 2，利用 St.1 和 St.2 計算出 Er ，其中e 為誤差限，Step4：記 Er = ma

32、x i,h/6i,h/ 2 e ×Yal(i) 8+ h × f (G ,Y t )Yal(i) = Y t；iikEr £ 1是否成立，若成立則表示滿足精度要求，此步Step5：；否則縮s × h × Er -0 25小步長h =，重新帶入進行 St.3 計算。2Step6：在截斷中考慮到四階方法的截斷誤差，利用兩次的計算結果做周期外插，üNì1()t +1/ 2tt +1/ 2t +1/ 2t +1/ 2(t +1/ 2)= Y+Y- Yt +1/2 t +1/ 2t +1/2 得到Y和Y 與Y之間得Y、Y值。Y

33、7;ýi,h/3i,h/6i,h/6i,h/3i,h/6î15þi=1115()Step7：計算Y (t+1)= Y (t +1/ 2) + Y t +1/ 2 +Y- Yt +1/ 2t +1/ 2。這一結果能達到 5 階的精度要i,h/6i,h/3i,h/6求。為了選取下一步的嘗試步長，算法中要選取誤差界 Er0 。s × h × Er-0 25Er > Er0 是否成立，成立則取Step8：，否則步長為3h 2 。2Step9：通過表 5.1 的初值；依次通過 St.1- St.7 迭代解得微分方程的解。利用 AFP-4-OSRK

34、算法保持了采樣時刻的取得更精確的值，從而通過更精準的估計為空間飛行器的無源探測做好較高精度的鋪墊。軌跡的5.1.2計算結果及分析根據初始數據表 5.1 和步長h = 0.2 ， Er £ 1´10-6 條件， Er 的取值來源于題目所要00最終結果保留 6 位有效數字。利用三維位置，具體結果如表 5.2 所示。編程求得在 50s、100s、150s、200s、250s 的表 5.250s、100s、150s、200s、250s 的三維位置利用 AFP-4-OSRK 算法求解運動微分方程模型，得到 09 號的軌跡（從 t=0s至 t=250s）如圖 5.1(a)所示，其對應時

35、刻的速度如圖 5.1(b)所示。從圖 5.1(a)可以看出 09號的空間軌跡近似一空間連續曲線。而從圖 5.1(b)和圖 5.1(d)可以看出 09 號的速度隨時間是連續變化的，速度沿 X 軸方向上變化最小，沿 Y 軸變化的幅度最大。圖5.1(c)表明在Y 軸方向上的分量幾乎不變，即的運動軌道平面近乎平行于 XOZ9時刻/sX 方向大小/mY 方向大小/mZ 方向大小/m501.77381´1068.16138´1064.51669 ´1061001.50163´1068.12676 ´1064.68468´1061501.22770

36、´1068.08270 ´1064.84722 ´1062009.52349 ´1058.02925´1065.00413´1062506.75894 ´1057.96650 ´1065.15524 ´106面。x 10636005.2340054.832004.630004.4280004.28.2-5350-50028.1-540061.5-5450x 10-1000816x 10-55000.5-1500-5550vy/m/s7.9vx/m/sY/mX/m(a)(b)空間軌跡速度曲線x 106隨時間

37、的位置變化隨時間的速度變化10500084000x軸方向分量y軸方向分量z軸方向分量3000x軸方向分量 y軸方向分量 z軸方向分量620004100002-10000-2000-3000-2-4000-4050100150t/s200250300-5000050100150t/s200250300(c)(d)速度隨時間的變化位置隨時間的變化圖5.109號從t=0s至t=250s間隔1s的空間軌跡和速度曲線5.2問題二的建模與求解利用 5.1 節中的附件給出的不同軌跡估計方法求出 06 號的軌跡，結果如圖 5.2 所示。由于觀測數據時間上不同步，如 meadata_06_00.txt 中觀測的

38、初始時刻為 50.175444807276s，而 meadata_06_00.txt 中觀測的初始時刻為 50.181218909269s，所以對 09 與 06 號的觀測數據要進行數據處理，使得所有觀測數據實現同步，從而確G定飛行器同一時刻在不同中的觀測數據，完成雙星逐點交匯。其次，根據vr (t) 和m(t) 的物理意義和變化規律，對二者分別建立了含有未知參量的近似模型，完善了主動段空間飛行器的運動方程，直接將飛行器的軌跡、速度和度數據代入運動方程可10Z/m分量/mvz/m /s分量/m/s實現模型中未知參數的估計，為簡化參數估計，提出了一種基于級數理論的多項式模型G來擬合v (t)m(

39、t) / m(t) 這一整體，給出了對多項式各階次系數的估計方法，并對不同階 r次多項式模型的性能進行了分析對比。再次，通過雙星逐點交匯得出的飛行器位置對建立的模型中各個參數進行估計，得到基礎坐標系下的主動段的簡化運動方程的明確形式，然后采用四階龍格庫塔來估計出飛行器主動段的軌道。利用估計出的飛行器軌跡與觀測到的軌跡進行比較，計算出殘差。該問題可以分為以下步驟：觀測數據同步、雙星逐點交匯x 10606號5、參量建模與估計、飛行器軌跡估計建模和殘差計算。軌跡45004.54400430048.242008.156x 104100-15002.52-41008.1-20006-4200x 101.

40、5-4300-2500Y/m-44008.05X/m1-3000Vy(t)-4500Vx(t)(a)(b)空間軌跡速度曲線圖5.206號從t=0s至t=250s間隔1s的空間軌跡和速度曲線5.2.1觀測數據同步5.2.1.1樣條插值對異步數據的同步算法主要有最小二乘擬合、多項式擬合、非線性擬合、拉格朗日插值、樣條插值等，各種算法都有相應的優缺點。采用樣條插2值繪制的曲線不僅有很好的光滑度，而且當節點逐漸加密時，其函數值在整體上能很好地逼近函數，相應的導數值也收斂于測數據進行同步。函數的導數，發生龍格現象。因此采用樣條插值法對觀根據附件 meadata_i_j.txt 給出的 06 與 09 號

41、的不同時刻觀測數據，分別畫出 06與 09 號各觀測時刻的觀測量a 、b 的散點圖，如圖 5.3 所示。11Z/mVz(t)：m/s0.120.160.110.140.10.120.090.10.080.080.070.060.06406080100120140160180406080100120140160180時間/s觀測量a 的散點圖時間/s觀測量 b 的散點圖(a) 09 號(b) 09 號-0.60.60.58-0.6050.56-0.610.54-0.6150.52-0.620.5-0.6250.480.46-0.63406080100120時間/s1401601804060801

42、00120140160180時間/s觀測量a 散點圖圖 5.306 與 09 號觀測量 b 的散點圖(c) 06 號(d) 06 號觀測量a 、 b 的散點圖采用樣條插值法對驟如下：觀測數據進行同步，這里以三次樣條插值的為例，具體的步Step1：給定區間a, b 的一個劃分a = t0 < t1 <" < tn-1 < tn = b ，a (t) 是a 的實際函數fa (t) 在該區間上關于該劃分的一個三次樣條函數。設a (t) 在節點ti 處的二階導數為：a ''(ti ) = Mi(i = 0,1,", n) ， 因為在子區間ti

43、-1ti 上a (t) = ai (t) 是三次多項式,所以a ''(t) 在此小區間上是t 的線性函數,且因為用線性插值,可知其表達式為：t Îti-1, ti t - ti-1a ''(ti-1 ) = Mi-1 ,a ''(ti ) = Miti - ta '' (t) = M+ Mii-1 ti t- t- ti-1iii-1ti - t + Mt - ti-1 連續兩次記h = t - t，則有a '' (t) = M得：iii-1i-1iihhii(t - t)3(t - t)3a (t) =

44、 M i+ Mi-1 + A (t- t) + B (t - t)(5.2-1)i-1ii-1iii i6h6hii12其中, Ai , Bi 為常數,可利用插值條件確定,即要求 Ai , Bi 滿足a (t) = 1 M1a (t ) =M h + B h = f (t )h2 + Ah = f (t),2i-1i-1 ii ia i-1i ii ia ii66并記 fa (ti-1 ) = yi-1, fa (ti ) = yi ，則得1 æ yö- 1 MhA =2ç i-1i-1 i ÷ihi è6ø1 æ y

45、46;1B =- M h2ç ii i ÷ihi è6ø(ti - t)3(t - ti-1 )32 ö (ti - t)æ2 ö (t - ti-1)æMi-1Miai (t) = Mi-1+ Mi+ ç yi-1 -+ ç yi -hihi÷÷6h6hè6øhè6øhiiii(t Îti-1, ti , i = 1, 2", n)Step2：由上討論可知,只要確定這 n+1 個值, 就可定出三樣條插值函數a (t

46、) 。為了求出Mi (i = 0,1,", n) ,利用一階導數在子區間連接點上連續的條件a '(ti - 0) = a '(ti + 0) ，求導一次,得在區間ti-1ti 上的表達式為：(t - t)2(t - t)2y - yha ' (t) = -M i+ Mi-1+ ii-1 - i (M- M)(5.2-2)ii-1iii-12h2hh6iii也就是在右端點ti 上有：yi - yi-1yi - yi-1- 0) = hi M - hi (M - M= hi6+ hi3a ' ) +M +(tMi iiii-1i-1i26hhii在左端點t

47、i+1 上有：yi - yi-1yi - yi-1+ 0) = - hi M- hi (M - M= - hi M3- hi6a ' (t) +M +i i-1i-1ii-1i-1i26hhiiStep3：將上式中的 i-1 改為 i,即得在子區間ti-1ti 上的表達式a 'i+1 (t) ,并由此得yi+1 - yia '(t + 0) = - hi+1 M- hi+1 M6+i+1 iii+13hi+1利用a '(t) 在內接點的連續性,即a 'i (ti - 0) = a 'i+1 (ti + 0) 就可得到關于參數的hi + hi+1y

48、i+1 - yiyi - yi-1hihi+16M, M , Mi+1 一個方程，6 Mi-1 +Mi +Mi+1 =-上式兩邊同乘i-1i3hhi+1i136以,即得方程hi + hi+1æ yi+1 - yiy - yöhihi+16(i = 1, 2,", n -1)M+ 2M +M=-ii-1ç÷i-1i+1h + hih + hh + hhhii+1 èi+1iøii+1ii+1ìmhi=ï ih + hïii+1ïhi+1若記 l = 1- mí則所得方程可簡寫成

49、：ih + hiïii+16ï( fa ti , ti+1 - fa ti-1 , ti )mi Mi-1 + 2Mi + li Mi+1 = giìm1M 0 + 2M1 + l1M 2 = g1ïgi =îhi + hi+1(i = 1, 2,", n -1)ïmM + 2M + l M = gï212232(5.2-3)即：íï"ïîmn-1M n-2 + 2Mn-1 + ln-1Mn = gn-1這是一個含有 n+1 個未知數、 n-1 個方程的線性方程組.

50、 要完全確定Mi (i = 0,1,", n) 的值還需要補充兩個條件,這兩個條件通常根據實際問題的需要，根據插值區間a,b的兩個端點處的邊界條件來補充。即已知插值區間兩端的一階導數值a '(t0 ) = fa¢ (t0 ),a '(tn ) = fa¢(tn ) ，則可得到包含ai 的兩個線性方程, a (t) 在子區間t0 , t1 上(t - t)2(t - ty - yh)2的 導 數 為 a ' (a ) = -M 1+ M0+ 10 - 1 (M- M )由 條 件101102h2hh6111h1 + y1 - y06 ( y1

51、 - y0- h1 (M - M ) 即2M + M =a '(t ) = f ¢(x ) = y¢ 得 y¢ = -M- y¢ )0a0000100102h6hh111yn - yn-1 )6Step4：同理由條件a '(t ) = f ¢(t ) = y¢ 得M+ 2M =( y¢ -nannn-1nnhhnn即得確定M 0 , M1,", Mn 的線性方程組é 212%ù é M 0 ùéùg0êmú ê

52、úêúlMgêêêêêëú êúêú1111ú ê#ú = ê #ú%21(5.2-4)ú êMúê gúmln-1 ú ên-1 ún-1 úên-1úû êë Mnúûêëúû2gn14ìg

53、39;= 6 ( f ht , t - y¢ )a0 1001其中， í6ïgt)=( y¢ - f, tan-1 nïnnhîn利用線性代數知識,可以證明方程組的系數矩陣都是非奇異的，因此有惟一解。求解確定a (t) 。同理可求b (t) 。對 06 與 09 號數據同步。的觀測數據分別進行樣條插值，完成5.2.1.2計算結果及分析對附件中的兩個觀測數據進行同步處理時，分別選取不同次數 n 的樣條插值對其插值，并對不同次數 n 樣條插值得到的數值與真實數值的偏差即進行比較，最后確認較優的 n 值。表 5.309 號觀測數據不同次數插值結果分析根據表 5.3 中的數據可以看出當次數取到 4 時插值偏差的均值的數量級就保持不變，所以綜合考慮精度與復雜度對觀測