林 夏，林寶軍，劉迎春，白 濤，武國(guó)強(qiáng)

（1. 上海微小衛(wèi)星工程中心，上海 201210；2. 中國(guó)科學(xué)院空天信息創(chuàng)新研究院，北京 100094；3. 中國(guó)科學(xué)院大學(xué)，北京 100094； 4. 中國(guó)科學(xué)院微小衛(wèi)星創(chuàng)新研究院，上海 201210；5. 上?？萍即髮W(xué)，上海 201210；）

天文導(dǎo)航是一類重要的衛(wèi)星自主導(dǎo)航方法。通過對(duì)天文信息的連續(xù)觀測(cè)，衛(wèi)星可實(shí)現(xiàn)自身位置、速度等導(dǎo)航信息的自主確定[1,2]。利用星敏感器與地球敏感器觀測(cè)信息的星光導(dǎo)航是其中一種適用于地球衛(wèi)星的天文導(dǎo)航方法。星光導(dǎo)航法實(shí)現(xiàn)方法簡(jiǎn)單，穩(wěn)定可靠，在低、中、高軌衛(wèi)星中都有著較多的應(yīng)用[3]。

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

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

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

1 星光導(dǎo)航算法原理

星光導(dǎo)航算法利用星敏感器與地球敏感器對(duì)衛(wèi)星星光角距信息進(jìn)行測(cè)量，并應(yīng)用Kalman濾波算法對(duì)衛(wèi)星動(dòng)力學(xué)方程獲得的一步軌道信息進(jìn)行修正，最終得到衛(wèi)星導(dǎo)航信息的最優(yōu)輸出。本節(jié)將從衛(wèi)星動(dòng)力學(xué)模型，算法觀測(cè)模型與算法濾波模型三方面對(duì)星光導(dǎo)航算法進(jìn)行詳述。

1.1 衛(wèi)星動(dòng)力學(xué)模型

依據(jù)牛頓定律，衛(wèi)星動(dòng)力學(xué)方程可表示為：

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

式中，μ=Gm為地球引力常數(shù)，G為萬(wàn)有引力常數(shù)，m為地球質(zhì)量。aNS為地球非球形引力，aNB為N體攝動(dòng)力加速度，主要為日月引力加速度，aSRP為衛(wèi)星光壓攝動(dòng)力加速度，aDG為大氣阻力加速度。

1.2 觀測(cè)模型

星光導(dǎo)航算法一般通過觀測(cè)星光角距信息進(jìn)行自主導(dǎo)航，觀測(cè)模型如圖1所示，具體模型構(gòu)造方法描述如下。

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

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

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

同時(shí)，為對(duì)衛(wèi)星導(dǎo)航信息進(jìn)行修正，建立衛(wèi)星位置信息r與星光角距信息θ的關(guān)系，形成觀測(cè)方程如下：

1.3 濾波模型

由于無(wú)跡卡爾曼濾波（UKF）[10]算法不需對(duì){ 非線性系統(tǒng)做簡(jiǎn)化近似，相比于擴(kuò)展卡爾曼濾波算法（EKF），無(wú)跡卡爾曼濾波（UKF）算法對(duì)于非線性系統(tǒng)有更好的表現(xiàn)。因此對(duì)于濾波模型的處理，本文采用UKF算法。

2 基于星光導(dǎo)航的地球扁率誤差自主修正

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

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

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

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

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

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

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

2.2 地球扁率誤差自主修正方法設(shè)計(jì)

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

地球?yàn)橐粋€(gè)南北方向較扁的橢球，在零經(jīng)度輔助慣性坐標(biāo)系下的地球球面方程如式(11)所示。

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

顯然，由于地球的扁球特性，地球敏感器掃入掃出點(diǎn)距地球中心的地心距并非為標(biāo)準(zhǔn)的地球赤道半徑Re，因而當(dāng)利用式(12)計(jì)算天底角時(shí)，計(jì)算結(jié)果將不可避免地存在一定偏差。

基于此，本節(jié)設(shè)計(jì)的地球扁率自主修正方法如下：

首先，采用天文導(dǎo)航算法得到連續(xù)穩(wěn)定的衛(wèi)星位置信息，該信息引入式(13)中，可計(jì)算出地球敏感器穿越點(diǎn)坐標(biāo)。式(13)可利用地球敏感器敏感地平模型得出，一般采用牛頓迭代法對(duì)式(13)求解[6]。

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

式中，Rsi為慣性坐標(biāo)系到星敏坐標(biāo)系的轉(zhuǎn)換矩陣，其可通過星敏感器實(shí)時(shí)輸出慣性系四元數(shù)信息轉(zhuǎn)換得到。

3 兩步Kalman濾波星光導(dǎo)航算法設(shè)計(jì)

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

圖3 兩步Kalman濾波星光導(dǎo)航算法設(shè)計(jì)圖Fig.3 Design of the starlight navigation algorithm based on the two-step Kalman filter

4 仿真分析

本節(jié)對(duì)設(shè)計(jì)的兩步Kalman濾波星光導(dǎo)航算法進(jìn)行仿真分析，通過與傳統(tǒng)的星光導(dǎo)航算法比較，評(píng)估設(shè)計(jì)算法的自主導(dǎo)航精度以及地球扁率誤差修正效果。

4.1 仿真場(chǎng)景建立

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

表1 算法仿真參數(shù)Tab.1 Simulation parameters

同時(shí)，衛(wèi)星采用偏航機(jī)動(dòng)模式調(diào)整姿態(tài)，即其+Z軸指向地心，并且太陽(yáng)矢量一直保持在衛(wèi)星本體系+XOZ面內(nèi)。仿真采樣周期設(shè)定為4 s，仿真時(shí)長(zhǎng)為43200 s。

4.2 仿真結(jié)果分析

首先對(duì)地球扁率引起的地球敏感器敏感地平誤差進(jìn)行評(píng)估。如圖4所示，仿真得到的地球敏感器測(cè)量的俯仰角與滾動(dòng)角最大誤差分別為0.327 °和0.116 °。根據(jù)地球敏感器俯仰角、滾動(dòng)角與單位地心矢量的關(guān)系，可相應(yīng)計(jì)算得到由地球扁率引起的單位地心矢量誤差。

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

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

圖5 傳統(tǒng)星光導(dǎo)航算法位置誤差圖Fig.5 Position errors in the traditional celestial navigation method

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

圖6 傳統(tǒng)星光導(dǎo)航算法速度誤差圖Fig.6 Velocity errors in the traditional celestial navigation method

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

最后，將本文設(shè)計(jì)的兩步Kalman濾波星光導(dǎo)航算法與傳統(tǒng)的星光導(dǎo)航算法得到的地心矢量精度及衛(wèi)星軌道信息精度分別進(jìn)行比較。表2為兩算法比較結(jié)果。相比于傳統(tǒng)的星光導(dǎo)航算法，利用本文設(shè)計(jì)的兩步Kalman濾波星光導(dǎo)航算法，地心矢量精度及衛(wèi)星軌道信息精度均提高了80%以上。因而，兩步Kalman濾波星光導(dǎo)航算法的有效性得到了有效地驗(yàn)證。

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

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

表2 傳統(tǒng)星光導(dǎo)航算法與兩步星光導(dǎo)航算法精度對(duì)比Tab.2 Comparison between traditional starlight navigation algorithm and the two-step starlight navigation algorithm

5 結(jié) 論

本文設(shè)計(jì)了一種基于兩步Kalman濾波的星光導(dǎo)航方法，得到了以下結(jié)論。

1）依據(jù)地球扁率誤差模型，設(shè)計(jì)了一種基于星光導(dǎo)航的地球敏感器地球扁率誤差修正方法。利用星光導(dǎo)航算法得到的衛(wèi)星位置信息，實(shí)現(xiàn)了地球敏感器地球扁率誤差的星上自主修正。

2）針對(duì)星光導(dǎo)航算法精度較低問題，設(shè)計(jì)了兩步Kalman濾波星光導(dǎo)航算法。通過第一步星光導(dǎo)航算法可得到精確的星光角距觀測(cè)信息，而后再次應(yīng)用星光導(dǎo)航濾波器可得到精度較高的衛(wèi)星軌道信息。

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