考慮地球扁率誤差的兩步Kalman濾波星光導航方法

林 夏，林寶軍，劉迎春，白 濤，武國強

（1. 上海微小衛星工程中心，上海 201210；2. 中國科學院空天信息創新研究院，北京 100094；3. 中國科學院大學，北京 100094； 4. 中國科學院微小衛星創新研究院，上海 201210；5. 上海科技大學，上海 201210；）

天文導航是一類重要的衛星自主導航方法。通過對天文信息的連續觀測，衛星可實現自身位置、速度等導航信息的自主確定[1,2]。利用星敏感器與地球敏感器觀測信息的星光導航是其中一種適用于地球衛星的天文導航方法。星光導航法實現方法簡單，穩定可靠，在低、中、高軌衛星中都有著較多的應用[3]。

但由于地球敏感器敏感地平時存在較大誤差，因而相比于其他的自主導航方法，星光導航法精度較低，在使用上也受到了許多限制[4]。地球敏感器敏感地平誤差由季節變化、地球扁率、環境溫度、安裝偏差、隨機噪聲等多種因素引起[5]。而在這些誤差源中，地球扁率造成的影響最為顯著。學者Tekawy[6]、Li Jie[7]、楊博[8]等人均針對地球扁率誤差提出了相應的修正方法，但這些算法均需要已知衛星軌道信息，并不適合在自主導航算法中直接使用。文獻[9]提出了一種不基于軌道先驗信息的地球扁率誤差修正方法，但算法基于雙圓錐掃描地球敏感器模型設計，并不適合在其他類型的地球敏感器中應用。

為實現地球扁率誤差的星上自主修正，從而提高基于星敏感器與地球敏感器觀測信息的星光導航算法精度，本文提出了一種基于兩步Kalman濾波的星光導航實現方法。首先，基于粗略星光角距觀測信息，利用Kalman濾波算法得到含有一定誤差的衛星粗軌道信息；而后，利用粗軌道信息計算地球敏感器的地球扁率誤差以實現誤差自主補償，得到精度較高的星光角距觀測信息；最后，基于修正后的星光觀測信息，再次利用Kalman濾波算法對衛星軌道信息進行估計，獲得精度較高的衛星自主導航信息。

本文首先介紹了星光導航算法的原理，并依據地球敏感器敏感地平原理研究了地球扁率誤差的自主修正方法，而后詳細給出了兩步Kalman濾波星光導航算法的設計方法，最后通過仿真驗證了算法的有效性。

1 星光導航算法原理

星光導航算法利用星敏感器與地球敏感器對衛星星光角距信息進行測量，并應用Kalman濾波算法對衛星動力學方程獲得的一步軌道信息進行修正，最終得到衛星導航信息的最優輸出。本節將從衛星動力學模型，算法觀測模型與算法濾波模型三方面對星光導航算法進行詳述。

1.1 衛星動力學模型

依據牛頓定律，衛星動力學方程可表示為：

式中，r、v分別為衛星的位置、速度向量。a為衛星加速度向量，可由式(2)表示。wr、wv為位置速度向量的過程噪聲信息，可視為白噪聲向量。

式中，μ=Gm為地球引力常數，G為萬有引力常數，m為地球質量。aNS為地球非球形引力，aNB為N體攝動力加速度，主要為日月引力加速度，aSRP為衛星光壓攝動力加速度，aDG為大氣阻力加速度。

1.2 觀測模型

星光導航算法一般通過觀測星光角距信息進行自主導航，觀測模型如圖1所示，具體模型構造方法描述如下。

圖1 星光角距觀測模型示意圖Fig.1 Observation model of starlight angle

首先，衛星可利用星敏感器觀測得到在星敏感器坐標系下的單位恒星矢量結合星敏感器安裝矩陣Rbs及式(3)，可得到衛星本體坐標系下單位恒星矢量表示。

并且，衛星可利用地球敏感器觀測得到地球敏感器坐標系下的單位地心矢量表示，結合地球敏感器安裝矩陣Rbe及式(4)，可得到地心指向在衛星本體坐標系下單位地心矢量表示。

同時，為對衛星導航信息進行修正，建立衛星位置信息r與星光角距信息θ的關系，形成觀測方程如下：

1.3 濾波模型

由于無跡卡爾曼濾波（UKF）[10]算法不需對{ 非線性系統做簡化近似，相比于擴展卡爾曼濾波算法（EKF），無跡卡爾曼濾波（UKF）算法對于非線性系統有更好的表現。因此對于濾波模型的處理，本文采用UKF算法。

2 基于星光導航的地球扁率誤差自主修正

2.1 地球敏感器敏感地平原理

以圓錐掃描地球敏感器為例，地球敏感器通過掃描地球采集掃入地球點Pin與掃出地球點Pout，確定敏感器坐標系下地心矢量指向。地球敏感器敏感地平模型如圖2所示，具體求解方法如下所述。

圖2 地球敏感器敏感地平模型圖Fig.2 Scanning model of the horizon sensor

地球敏感器繞掃描軸勻速轉動，當捕獲地球時，可獲得掃入地球時間tin，掃出地球時間tout。通過與掃過基準點時間tr相減，可獲得地球敏感器敏感地平的掃描相角Uin、Uout，具體計算方法如式(8)所示[11]。

式中，Ts為地球敏感器的掃描周期。

假設地球敏感器掃描圓錐軸心沿衛星本體系+Y方向，且地球敏感器敏感姿態角—俯仰角與滾動角較小，那么衛星的俯仰角?與滾動角φ可通過式(9)計算得到。相應地，地心矢量在敏感器坐標系下的表示可由式(10)計算得到[12]。

式中，ηin和ηout分別為地球敏感器掃入點天底角與掃出點天底角，T為地球敏感器的半錐角。

2.2 地球扁率誤差自主修正方法設計

由式(9)(10)可知，地球敏感器敏感地平精度由掃描相角Uin、Uout，敏感器半錐角T及天底角ηin、ηout精度確定。由于掃描相角及半錐角均可精確測量，因而地球敏感器敏感地平精度取決于天底角精度。地球扁率誤差為影響天底角精度的主要因素。由于地球扁率誤差修正需已知衛星位置信息，本節設計了基于星光導航信息的地球扁率誤差自主修正方法。

地球為一個南北方向較扁的橢球，在零經度輔助慣性坐標系下的地球球面方程如式(11)所示。

式中，e為地球扁率系數，Re為赤道半徑。

顯然，由于地球的扁球特性，地球敏感器掃入掃出點距地球中心的地心距并非為標準的地球赤道半徑Re，因而當利用式(12)計算天底角時，計算結果將不可避免地存在一定偏差。

基于此，本節設計的地球扁率自主修正方法如下：

首先，采用天文導航算法得到連續穩定的衛星位置信息，該信息引入式(13)中，可計算出地球敏感器穿越點坐標。式(13)可利用地球敏感器敏感地平模型得出，一般采用牛頓迭代法對式(13)求解[6]。

式中，P為地球敏感器掃入掃出地平的位置信息，S為衛星位置，O為地心原點。因而，OP為穿越點慣性坐標系位置向量，OS為衛星地心矢量，0SP為衛星到穿越點的慣性坐標系下的單位向量表示。同時，D=diag(1,1,1/(1 -e)2)，Ys為地球敏感器掃描軸在敏感器坐標系的單位向量表示，Rei為慣性坐標系到地球敏感器坐標系轉換矩陣，可通過式(14)計算得到。式(13)中所有符號也均在圖2中加以標注。

式中，Rsi為慣性坐標系到星敏坐標系的轉換矩陣，其可通過星敏感器實時輸出慣性系四元數信息轉換得到。

3 兩步Kalman濾波星光導航算法設計

基于上節所述地球扁率誤差自主修正方法，本節設計了兩步Kalman濾波算法以將地球扁率誤差修正算法嵌入星光導航算法中，從而實現地球扁率誤差的自主修正及星光導航算法精度的提升。兩步Kalman濾波星光導航算法設計如圖3所示。首先，算法利用地球敏感器原始輸出的地心矢量信息與星敏感器輸出的星光指向信息形成粗星光角距信息。將粗星光角距信息與所需軌道位置速度、協方差預報信息引入星光導航濾波器中得到存在一定誤差的衛星粗位置信息。算法將衛星粗位置信息引入2.2節描述的地球扁率誤差修正方法中，修正地球敏感器輸出，從而獲得精確的地心矢量信息。而后，將得到的精地心矢量信息結合星敏感器輸出獲得精確的星光角距觀測信息，并將精星光角距信息再次代入星光導航算法濾波器中對衛星軌道信息進行濾波求解，最終可得到精確的衛星軌道信息及相應的協方差信息的更新輸出。

圖3 兩步Kalman濾波星光導航算法設計圖Fig.3 Design of the starlight navigation algorithm based on the two-step Kalman filter

4 仿真分析

本節對設計的兩步Kalman濾波星光導航算法進行仿真分析，通過與傳統的星光導航算法比較，評估設計算法的自主導航精度以及地球扁率誤差修正效果。

4.1 仿真場景建立

仿真中選用的衛星軌道為一LEO衛星軌道，軌道參數及敏感器參數設置如表1所示。

表1 算法仿真參數Tab.1 Simulation parameters

同時，衛星采用偏航機動模式調整姿態，即其+Z軸指向地心，并且太陽矢量一直保持在衛星本體系+XOZ面內。仿真采樣周期設定為4 s，仿真時長為43200 s。

4.2 仿真結果分析

首先對地球扁率引起的地球敏感器敏感地平誤差進行評估。如圖4所示，仿真得到的地球敏感器測量的俯仰角與滾動角最大誤差分別為0.327 °和0.116 °。根據地球敏感器俯仰角、滾動角與單位地心矢量的關系，可相應計算得到由地球扁率引起的單位地心矢量誤差。

圖4 地球敏感器滾動俯仰誤差圖Fig.4 Roll and pitch errors in the horizon sensor

而后利用仿真得到的單位地心矢量信息，可應用星光導航算法，對傳統星光導航算法的性能進行評估。以精確的軌道外推模型為基準，可精確確定傳統的星光導航算法誤差。如圖5所示，三軸位置誤差分別為26045 m，32553 m，22852 m；三軸速度誤差分別為21.930 m/s，20.475 m/s，24.811 m/s。可以看到，地球扁率誤差嚴重影響了星光導航算法的自主導航精度。

圖5 傳統星光導航算法位置誤差圖Fig.5 Position errors in the traditional celestial navigation method

隨后，應用本文設計的兩步Kalman濾波星光導航算法，對敏感器地球扁率誤差自主修正精度及星光導航算法精度進行評估。圖6、圖7分別為經兩步Kalman濾波算法計算后的修正后地心矢量誤差及星光導航算法誤差。如圖6所示，經兩步Kalman濾波算法，修正后的地球敏感器滾動角及俯仰角最大誤差分別為0.0011o、0.0006o，地心矢量指向被精確確定。同時，利用修正后的單位地心矢量，可得到精確的星光角距觀測信息，因而星光導航算法精度也得到極大提升。如圖7所示，待算法收斂后，三軸最大位置誤差分別為1508.2 m，947.0 m，2836.8 m；三軸最大速度誤差分別為3.721 m/s，3.540 m/s，4.146 m/s。

圖6 傳統星光導航算法速度誤差圖Fig.6 Velocity errors in the traditional celestial navigation method

圖7 修正后地球敏感器滾動俯仰誤差圖Fig.7 Roll and pitch error in the horizon sensor after revised

最后，將本文設計的兩步Kalman濾波星光導航算法與傳統的星光導航算法得到的地心矢量精度及衛星軌道信息精度分別進行比較。表2為兩算法比較結果。相比于傳統的星光導航算法，利用本文設計的兩步Kalman濾波星光導航算法，地心矢量精度及衛星軌道信息精度均提高了80%以上。因而，兩步Kalman濾波星光導航算法的有效性得到了有效地驗證。

圖8 兩步天文導航算法位置誤差圖Fig.8 Position errors in the two-step celestial navigation method

圖9 兩步天文導航算法速度誤差圖Fig.9 Velocity errors in the two-step celestial navigation method

表2 傳統星光導航算法與兩步星光導航算法精度對比Tab.2 Comparison between traditional starlight navigation algorithm and the two-step starlight navigation algorithm

5 結 論

本文設計了一種基于兩步Kalman濾波的星光導航方法，得到了以下結論。

1）依據地球扁率誤差模型，設計了一種基于星光導航的地球敏感器地球扁率誤差修正方法。利用星光導航算法得到的衛星位置信息，實現了地球敏感器地球扁率誤差的星上自主修正。

2）針對星光導航算法精度較低問題，設計了兩步Kalman濾波星光導航算法。通過第一步星光導航算法可得到精確的星光角距觀測信息，而后再次應用星光導航濾波器可得到精度較高的衛星軌道信息。

3）仿真驗證了算法的有效性。利用本文算法，地球敏感器敏感地心矢量精度由0.3 °提升至0.001 °，衛星三軸位置精度由30000 m提升至3000 m，衛星三軸速度精度由20 m/s提升至4 m/s。因而，相比于傳統的無修正星光導航算法，提升了80%以上。