一種小天體繞飛軌道及目標天體參數確定方法

王衛東，張澤旭，朱圣英，崔平遠

(1．哈爾濱工業大學 深空探測基礎研究中心，150001哈爾濱，zexuzhang@hit．edu．cn;2．北京理工大學宇航學院，100081北京)

一種小天體繞飛軌道及目標天體參數確定方法

王衛東1，張澤旭1，朱圣英2，崔平遠2

(1．哈爾濱工業大學 深空探測基礎研究中心，150001哈爾濱，zexuzhang@hit．edu．cn;2．北京理工大學宇航學院，100081北京)

提出一種基于SRIF濾波器的小天體探測器繞飛軌道及目標天體參數確定方法．考慮繞飛小天體初期待估參數多、動力學環境復雜等特點，利用SRIF濾波器對組合導航數據進行處理，對探測器繞飛軌道以及目標天體引力場模型、自旋狀態、星歷信息等動力學參數進行估計．利用Eros433的觀測數據進行數學仿真表明，本文給出的方法能夠有效地計算繞飛軌道和確定目標天體的物理參數．

小天體探測;繞飛軌道;光學導航;均方根信息濾波

當探測器最終捕獲目標小天體并形成穩定的繞飛軌道以后，需要確定探測器相對小天體的位置、速度以及姿態等參數，同時在此階段，對目標天體參數的評估也是對探測器的1個重要需求，這些參數包括影響探測器軌道的小天體的引力場系數、自旋角速度、自旋軸指向以及小天體形狀和星歷信息等．NEAR、Rosetta 任務［1－2］證實將小天體表面存在的彈坑地表特征作為導航路標，比利用小天體形狀特征具有更好的軌道確定性能．

在小天體探測器繞飛軌道確定方面，日本ISAS的Kawaguchi等研究了繞飛小行星等小天體的自主光學制導與導航技術［3－4］．美國 JPL 的Bhaskaran等研究了低成本的繞飛小天體軌道自主確定技術，提出了一種利用寬視場相機和預處理的目標小天體模型自主確定繞飛小天體軌道的算法［5－6］．J．J．Bordi利用 NEAR 任務中測量數據研究了激光測距儀的測量對Eros形狀評估和探測器軌道確定的影響，并對只借助于激光測據儀導航方式的導航精度進行了分析，結果顯示此時的導航精度在100 m左右［7］．

利用探測器繞飛小天體的觀測數據來確定小天體模型參數的研究，主要是JPL實驗室和ISAS在進行這方面的研究．JPL實驗室一般利用下面兩種技術，第一種技術是使用PCODP(PC Orbit Determination Program)軟件處理大約30 d的探測器斜距和多普勒數據及小天體表面的光學觀測數據［8－9］，第二種技術是利用另 1 個獨立的軟件ODP(Orbit Determination Program)處理整個軌道的斜距和多普勒數據［10］．Kominato 等［11］研究了利用在軌光學導航參數確定小天體相關物理參數的問題;Catherine等［12］利用飛越小天體的觀測數據研究了小天體的質量確定問題．Johnson等［13］人研究了基于視覺分析的火星軟著陸導航方法．

均方根信息濾波SRIF算法在處理包含了光學數據的觀測數據時能有效地克服濾波器的發散，具有較高的數值穩健型和計算高效性［14］．本文將利用SRIF對繞飛小天體的軌道確定濾波器進行設計，考慮未建模干擾加速度的影響，通過測距數據、測速數據、VLBI數據及導航路標光學數據的處理，對探測器繞飛軌道、小天體物理參數進行確定．

1 變分方程與觀測方程

定義小天體固聯坐標系Σa:坐標原點oa位于小天體的質量中心，za軸沿小天體自旋軸方向，xa軸沿小天體最小慣量軸方向，ya軸定義滿足右手系法則.在小天體固聯坐標系下，繞飛探測器的軌道動力學方程可以表示為

式中:x、y、z為探測器的三軸位置;ω為小天體的自旋角速度;Tx、Ty、Tz為建模加速度包括控制加速度、太陽光壓加速度等;nx、ny、nz為未建模的干擾加速度;V為小天體的引力位函數，利用球體調和函數表示如下:

其中:a為小天體名義半徑;θ、φ為探測器所處的經緯度;Pnm為締結勒讓德多項式.

定義小天體慣性坐標系ΣA:坐標原點oA位于小天體的質量中心，zA軸沿小天體自旋軸方向，xA軸位于黃道平面內垂直zA軸，yA軸定義滿足右手系法則.令小天體自旋軸在J2000坐標系下的赤經、赤緯分別為θa、φa，則小天體慣性坐標系相對J2000坐標系的坐標轉換矩陣為小天體固聯坐標系相對于小天體慣性系坐標轉換矩陣為

其中λ為某歷元時刻xa與xA之間夾角，t為當前時刻距歷元時刻的時間.

探測器繞飛軌道定軌數據有地面站測軌數據和星載相機拍攝的路標圖像數據，其中地面站測軌數據提供探測器絕對位置、速度信息，其本質包含了探測器運動距離信息．路標圖像數據提供了探測器相對目標天體的位置信息，其本質是角度測量信息，結合地面測軌數據包含的運動距離信息和星敏感器提供的姿態信息，能夠確定探測器在小天體固聯坐標系下的三維位置、運動狀態，同時建立小天體固聯坐標系與J2000坐標系之間的轉換關系．路標圖像數據為小天體表面彈坑相對應的像素p、像線l坐標值，其表達式為

其中xci、yci、zci是第i個導航路標在相機坐標系下相對探測器位置矢量rci的三軸分量，f為導航相機焦距.其中

式中r，ρi分別為小天體固聯坐標系下探測器和第i個導航路標的位置矢量;Cca是相機坐標系相對小天體固聯坐標系的坐標轉換矩陣，有Cca=CcbCbICIa成立，Ccb為導航相機坐標系在探測器本體坐標系下的安裝矩陣，CbI為探測器本體坐標系相對慣性空間的姿態轉換矩陣，CIa為J2000坐標系相對小天體固聯坐標系姿態轉換矩陣．

基于以上給出的軌道動力學方程與坐標轉換關系，建立變分方程、量測方程，可以求取狀態轉移偏導數與各觀測數據對應的觀測偏導數，進而確定量測偏導數，以對待估參數進行確定．

2 待估參數與SRIF濾波器設計

2.1 待估參數選取

參考小天體繞飛段軌道變分方程和量測方程，可見直接影響觀測量的因素有探測器相對目標天體的位置、速度和姿態、目標天體的運行軌道、測控站的三維位置、導航路標在目標天體固聯坐標系下的位置，間接影響觀測量的因素有太陽光壓系數、目標天體引力場系數、目標天體自旋狀態、隨機加速度等動力學參數．

在探測器繞飛目標天體階段，星敏感器與陀螺聯合定姿系統能夠提供探測器相對慣性空間的姿態信息，而導航路標所對應像點坐標與探測器相對目標天體的姿態相關，因此光學觀測數據對目標天體相對慣性空間姿態可觀．這里采用小天體固聯坐標系相對J2000姿態轉換矩陣表示目標天體的姿態，該姿態可以利用小天體自旋軸指向、子午線歷元時刻指向、小天體旋轉角速度表示．光學觀測數據同時與探測器相對目標天體的位置相關，地面測控數據則與探測器相對日心慣性空間軌道狀態有關，因此，結合光學觀測數據與地面測控數據能夠對目標天體的星歷信息進行估計，這里采用歷元時刻目標天體相對日心的位置、速度表示其星歷狀態．太陽光壓力、小天體引力、向心加速度、苛氏加速度等動力學參數影響探測器的運行軌道，間接影響各觀測量數據，其隱含在變分方程中，可以通過狀態轉移矩陣對這些動力學參數進行估計．綜合如上因素，考慮科學任務和工程任務的需求，本文選取SRIF濾波器待估參數如表1所示．

表1 軌道確定濾波器估計參數

在小天體表面存在大量彈坑地形，后續自主光學導航需要選取彈坑為導航路標，同時對其三維位置信息進行確定．待估參數中的隨機加速度項為探測器所受到的未建模加速度部分，主要是由于太陽光壓力變化、燃料泄漏引起的，參考國外測控經驗，用高斯－馬爾科夫過程對其建模，該加速度包括時間相關的部分和純粹隨機的部分，其在一定時間內變化不大．SRIF濾波器采用批處理與貫序處理相結合的濾波技術，在濾波過程中，每個批處理過程中均假設未建模加速度為常值，在批與批之間的貫序處理過程中時間相關變化．

基于SRIF濾波的基本原理，考慮到實現繞飛小天體軌道確定的復雜性，將待估參數分為三類進行處理:1)P為相關過程噪聲參數，即為上述中的未建模加速度;2)X為隨時間變化的狀態參數，但并不完全依賴于白噪聲的量，如探測器的相對位置和速度、小天體星歷狀態等;3)Y為不隨時間變化的狀態參數，如探測器太陽光壓系數、小天體引力場系數、小天體姿態狀態、彈坑三維位置等常值參數.此時，濾波器的狀態方程可表示為

式中:下標j表示第j批參數;Wj+1為濾波器的過程噪聲;VP(j)為未建模加速度對狀態參數的作用矩陣，可以表示為

其中ΦX和ΦXP代表了相應的狀態轉移矩陣.由于假設探測器所受未建模加速度滿足高斯馬爾科夫過程，因此式(1)中轉移矩陣Mj可以表示為

其中Δtj表示第j批與第j+1批間隔時間，τi表示噪聲的相關時間.

選取每步待估參數的更新值作為遞推狀態，令初始待估參數服從的誤差協方差陣為P0，對P0陣進行Cholesky變換獲取初始均方差陣S0，進而求得初始信息矩陣R0，滿足

式中R0為上三角陣，對應各待估參數可以表示為

式中 zX、zP、zY為初始虛擬觀測值，vX、vP、vY為初始觀測噪聲，由于初始改變值為零值，所以zX、zP、zY和 vX、vP、vY也均為零值.

2.2 時間更新方程

方程(1)可以等價表示為

將上式代入式(2)中，有

由于Wj+1可以認為是一獨立的隨機過程噪聲，并服從N( 0，σ2)分布，可以用下式描述:式中rw為相應過程噪聲的標準差倒數，考慮到Wj+1為零均值，將zW取為零值.綜合式(4)和式(5)并代入式(3)，有

利用Householder變換［15］構造正交陣T左乘上式，使左端動態矩陣轉化為上三角矩陣，上式變為

將式(6)第2行到第4行的子矩陣保留，即為時間更新到第j+1步的信息矩陣等式，上述過程為SRIF濾波器的時間更新過程，主要考慮了未建模加速度對系統估計參數的影響．

2.3 測量更新方程

利用地面測軌數據與光學數據的線性化量測偏導數，可以將測量殘差表示成相應待估參數的形式

考慮測量數據的噪聲特性，將上式與信息矩陣等式合并，有

其中N為觀測噪聲服從的誤差標準差，即R=NNT，R為觀測噪聲服從的誤差協方差陣.

利用Householder變換構造正交陣T左乘上式，使左端動態矩陣轉化為上三角矩陣，則上式變為

求解式(8)，可得到第j+1步的狀態參數以及相應的協方差陣，利用式(7)中量測更新的信息矩陣等式，繼續進行時間更新，以便數據的貫序批處理，完成探測器軌道狀態和目標天體物理參數的確定.

3 數學仿真與分析

為了驗證軌道確定方案的可行性以及SRIF濾波器的性能，以Eros433為目標星，利用從2012年12月28日起的14天觀測數據對其繞飛軌道與目標天體參數進行估計．小行星Eros433自旋角速度1639.4(°/d)，名義半徑16 km，引力常數4.462 1×105m3/s2．設探測器繞飛軌道高度為80 km，在該軌道高度探測器受到小天體引力場高階項的攝動已較為明顯，如圖1所示．在兩周的引力攝動作用下，探測器軌道漂移將會達到3 km以上，因此在該軌道高度上能夠對目標天體的引力場進行較高精度的評估．這里選取待估目標天體引力場6階系數，繞飛軌道偏心率 0，軌道傾角176°，探測器面質比0. 01，光壓系數為0.6．

參考國外測控經驗，SRIF濾波器采用一階高斯馬爾科夫過程模擬未建模加速度，其相關時間選為2 d，量級為10－12km/s．探測器在80 km 軌道上繞飛階段，地面站每天跟蹤測軌1次，每次2 h，導航相機每天拍照60幅．斜距測量系統差、隨機差均為10 m，斜距變化率測量系統差、隨機差為1 mm/s，VLBI測量測角精度為5 nrad，導航相機視場角為3°×3°，分辨率為1 024 ×1 024，圖像處理精度0.1 piexl，導航路標選取100個，星敏感器確定精度為50 nrad．待估參數初始不確定度及經SRIF濾波器迭代30次得到的最終確定結果如表2所示．

圖1 探測器繞飛小天體Eros433軌道

表2 待估參數初始條件及結果

初始探測器相對目標天體的不確定位置精度為三軸各300 m，速度0.01 m/s，經觀測處理后，位置、速度狀態確定誤差均得到減小，由結果可見，探測器在小天體固聯坐標系下z軸的導航精度明顯地要高于其它兩軸，這與探測器繞飛軌道傾角有關．在本仿真中，繞飛軌道近似為赤道面軌道，導航路標觀測幾何關系使自旋軸方向的可觀度大于其它兩軸．由目標天體自旋軸、子午線與旋轉角速度的估計結果可以看出，采用光學路標導航輔助的方法，可以對目標天體在絕對空間下的姿態指向進行較高精度的估計，其估計精度與星敏感器同等量級．小天體引力常數由初始2 000 m3/s2的不確定度，下降到100 m3/s2以內，估計精度達到了0.02%以內，對小天體各階引力系數估計精度也提高了2個數量級，說明濾波器對小天體質量精確確定的同時，也能夠對小天體的質量分布進行很好地估計，結合小天體形狀信息，通過對小天體密度的分析，可以對其內部組成成分及物質構成提供參考．另一項科學考察數據，小天體的星歷狀態可以達到10 km以內，速度在0.1 m/s以內，達到與基于地面站對深空探測器進行測軌的同等精度量級．SRIF濾波器對100個導航路標的三維位置進行估計，表2中誤差表示形式為位置誤差平均值，其表達式為

同時導航相機可以拍攝到其它路標的大量圖像數據，這些路標雖然沒有在SRIF濾波器中對其位置進行估計，但是結合拍照時刻探測器的位置、姿態等信息，利用多幀圖像數據能夠對這些路標的位置進行估計，這些路標與SRIF濾波器中采用的導航路標聯合組成了導航路標庫，以支持后續探測任務的自主光學導航的實現．

4 結論

本文對探測器繞飛軌道及目標天體參數確定方法進行了研究，針對繞飛小天體初期待估參數多，動力學不確定環境復雜等特點，同時考慮未建模加速度的影響，提出了一種基于SRIF濾波器的軌道確定和目標天體參數確定方法．該方法利用地面測控數據、導航路標的光學數據對探測器繞飛軌道確定，同時對目標天體引力場模型、自旋狀態、星歷信息等動力學參數進行估計．仿真實驗表明該方法能夠有效地對繞飛軌道和目標天體參數進行計算．

［1］KONOPLIV A S，MILLER J K，OWEN W M，et al.A global solution for the gravity field，rotation，landmarks，and ephemeris of eros［J］．Icarus， 2001，160(2):289－299．

［2］ BIELE J，ULAMEC S．Capabilities of philae，the rosetta lander［J］．Space Science Reviews， 2008，138(1/2/3/4):275－289．

［3］KAWAGUCHI J，HASHIMOTO T，MISU T，et al.An autonomous optical guidance and navigation around the asteroid［C］//Proceedings of the 47th International Astronautical Congress．Beijing:［s．n．］，1996．

［4］KAWAGUCHI J，HASHIMOTO T，KUBOTA T，et al.Autonomous optical guidance and navigation strategy around a small body［J］．Journal of Guidance，Control and Dynamics， 1997，20(5):1210 －1220．

［5］BHASKARAN S，RIEDEL J E，SYNNOTT S P．Demonstration of autonomous orbit determination around small bodies［J］．Advances in the Astronautical Sciences Series， 1996，90(2):1297 －1308．

［6］TIMOTHY C D．NEAR Laser rangefinder:a tool for the mapping and topologic study of asteroid 433 Eros［J］．Johns Hopkins APL Technical Digest， 1998，19(2):142－157．

［7］BORDI J J，MILLER J K，WILLIAMS B G．Near Earth Asteroid Rendezvous(NEAR)navigation using altimeter range observations［R］．The Interplanetary Network Progress Report， 2001，IPN PR －42－146．

［8］MILLER J K，WEEKS C J，WOOD L J．Orbit determination strategy and accuracy for a comet rendezvous mission［J］．Journal of Guidance，Control and Dynamics， 1990，13:775 －784．

［9］CHESLEY S R，YEOMANS D K．Comet 9P/Tempel 1 ephemeris development for the deep impact mission［J］．Advances in the Astronautical Sciences， 2006， 2:1271－1282．

［10］MILLER J K，KONOPLIV A S，ANTREASIAN P G，et al.Determination of shape，gravity，and rotational state of asteroid 433 Eros［J］．Icarus， 2002，155(1):3 －17．

［11］KOMINATO T，MATSUOKA M，UO M．Optical hybrid navigation in Hayabusa-Approach，station keeping＆hovering［J］．Advances in the Astronautical Sciences， 2006，124(2):1753－1772．

［12］CATHERINE B S，JEREMY B J，JOHN D A．Small body mass determination from a flyby［J］．Advances in the Astronautical Sciences， 1991，75(2):885 －896．

［13］JOHNSON A E，WILLSON R，CHENG Y．Design through operation of an image-based velocity estimation system for mars landing［J］．International Journal of Computer Vision， 2007，74(3):319－341．

［14］BIERMAN G J．Factorization Methods for Discrete Sequential Estimation．New York:Academic Press，1976．

［15］ALSTON S．Householder，Unitary Triangularization of a Nonsymmetric Matrix［J］．Journal ACM， 1958，5(4)，339－342．

An algorithm of around orbit and parameters determination for small body

WANG Wei-dong1，ZHANG Ze-xu1，ZHU Sheng-ying2，CUI Ping-yuan2

(1．Deep Space Exploration Research Center，Harbin Institute of Technology，150001 Harbin，China，zexuzhang@hit．edu．cn;2．Aerospace College，Beijing Institute of Technology，100081 Beijing，China)

The around orbit parameters and small body physical parameters determination methods are studied based on Square Root Information Filter(SRIF)．Due to a great deal of parameters to be estimated and the complexity of the dynamics，the SRIF is designed to process the integrated navigation data and estimate the gravity field，rotation，and ephemeris of small body and orbit parameters of probe．The mathematic simulation experiments using the measuring data of Eros433 small body have validated the algorithm．

small body exploration;around orbit;optical navigation;Square Root Information Filter(SRIF)

文獻標志碼:A 文章編號:0367－6234(2011)09－0019－06

2010－04－09．

航天創新基金資助項目(CASC200902－4)，深空探測著陸與返回控制技術國防重點學科實驗室，開放基金資助項目(HIT．KLOF．2009070)．

王衛東(1969—)，男，博士研究生;

崔平遠(1961—)，男，教授，博士生導師．

(編輯 張 宏)