基于擴展卡爾曼濾波的XPNAV-1衛(wèi)星自主定軌算法研究

丁陶偉，帥平，黃良偉，張新源

中國空間技術(shù)研究院 錢學(xué)森空間技術(shù)實驗室，北京 100094

X射線脈沖星屬于高速旋轉(zhuǎn)的中子星，具有極穩(wěn)定的周期性。利用X射線脈沖星導(dǎo)航能夠為航天器提供包括時間、位置、速度和姿態(tài)在內(nèi)的全面的導(dǎo)航參考信息，為解決近地軌道和深空探測航天器的持續(xù)高精度自主導(dǎo)航難題提供了一種新思路。近十余年來，利用X射線脈沖星的航天器自主導(dǎo)航一直是國內(nèi)外航天前沿技術(shù)研究的熱點領(lǐng)域[1-2]。2004年，美國國防部國防預(yù)先研究計劃局(DARPA)提出“X射線導(dǎo)航與自主定位驗證(XNAV)”計劃[3]，其目標(biāo)是建立一個脈沖星網(wǎng)絡(luò)，利用脈沖星輻射的X射線信號實現(xiàn)航天器長時間自主導(dǎo)航。自2005年以來，中國持續(xù)開展X射線脈沖星導(dǎo)航系統(tǒng)理論與方法研究，研制了多種類型的X射線探測器產(chǎn)品、脈沖星導(dǎo)航地面試驗系統(tǒng)和脈沖星導(dǎo)航專用試驗衛(wèi)星[4]。2016年11月10日，中國成功發(fā)射全球首顆脈沖星導(dǎo)航專用試驗衛(wèi)星——X射線脈沖星導(dǎo)航1號(XPNAV-1)，率先開展脈沖星導(dǎo)航空間飛行試驗[5-6]。2017年6月3日，美國宇航局(NASA)將“空間站X射線計時與導(dǎo)航技術(shù)試驗(SEXTANT)”設(shè)備發(fā)射到國際空間站(ISS)上，開展脈沖星導(dǎo)航在軌演示驗證試驗[7]。該試驗項目是“中子星內(nèi)部結(jié)構(gòu)探測器(NICER)”任務(wù)之一，其X射線計時儀的有效探測面積大于1 800 cm2，探測能譜范圍為0.2～1.2 keV。

在利用X射線脈沖星的航天器自主軌道確定算法研究方面，前人已經(jīng)進行了大量研究工作，主要包括針對模型不確定性的魯棒濾波算法[8]、為降低非線性影響的非線性預(yù)測濾波(NPF)算法[9]、針對脈沖星方向誤差的擴展?fàn)顟B(tài)無跡卡爾曼濾波(ASUKF)算法[10]、基于多模自適應(yīng)估計的無跡卡爾曼濾波(UKF)算法[11]等，但這些算法都是利用仿真數(shù)據(jù)來實現(xiàn)的。目前，美國戈達德太空飛行中心(GSFC)利用SEXTANT設(shè)備的觀測數(shù)據(jù)開展軌道確定試驗，但尚未檢索到其具體定軌算法。XPNAV-1衛(wèi)星在軌運行4年來，獲取了大量觀測數(shù)據(jù)，按計劃完成了X射線探測器性能測試、典型目標(biāo)脈沖星觀測及脈沖星導(dǎo)航系統(tǒng)體制驗證等試驗任務(wù)[12-14]。利用XPNAV-1衛(wèi)星對蟹狀星云(Crab)脈沖星的完整觀測數(shù)據(jù)，采用基于軌道法平面幾何約束方法，驗證了該衛(wèi)星軌道改進的有效性和脈沖星導(dǎo)航系統(tǒng)體制[15]，但利用單顆脈沖星的觀測數(shù)據(jù)長時間定軌過程存在發(fā)散問題。本文針對XPNAV-1衛(wèi)星拓展任務(wù)及后續(xù)發(fā)展需求，建立其軌道力學(xué)模型和觀測方程，并論述擴展卡爾曼濾波(EKF)算法和分段式定常系統(tǒng)(PWCS)的可觀測性分析方法。

1 XPNAV-1衛(wèi)星任務(wù)

XPNAV-1衛(wèi)星是一顆整星質(zhì)量為243 kg的小衛(wèi)星。其科學(xué)試驗?zāi)繕?biāo)是：1)在空間環(huán)境下實測驗證X射線探測器性能，研究宇宙背景噪聲對探測器作用機理；2)探測Crab脈沖星輻射的X射線光子，提取脈沖輪廓曲線；3)長時間累積探測脈沖星，積累觀測數(shù)據(jù)，探索驗證脈沖星導(dǎo)航體制。

XPNAV-1衛(wèi)星采用降交點地方時為6:00時的太陽同步晨昏軌道，其軌道半長軸為6 878.137 km，軌道傾角為97.4°，偏心率為0.010 3。該衛(wèi)星的有效載荷為兩種不同類型的X射線探測器。其中一種是掠入射式聚焦型探測器，其有效面積為30 cm2，探測能譜范圍為0.5～10 keV，時間分辨率為1.5 μs，視場為15′；另一種是微通道板準(zhǔn)直型探測器，其有效探測面積為1 200 cm2，探測能譜范圍為1～10 keV，時間分辨率為100 ns，視場為2°。針對核心試驗任務(wù)需求，優(yōu)選了8顆目標(biāo)脈沖星進行重點觀測研究，其脈沖星編號、J2000.0國際天球參考系下的赤經(jīng)和赤緯、脈沖周期及光子流量如表1所示。表1中，序號1～4為獨立旋轉(zhuǎn)供能脈沖星(IRP)，其中第1顆脈沖星B0531+21為Crab脈沖星，序號5～8為X射線雙星(XB)。

XPNAV-1衛(wèi)星在軌運行近4年來，對8顆目標(biāo)脈沖星、超新星遺跡以及銀河系X射線等天體進行了觀測研究。由表1可見，XB的光子流量比IRP高2～3個數(shù)量級，觀測XB的主要目的是進行探測器性能測試?？紤]到XB系統(tǒng)本身軌道動力學(xué)的復(fù)雜性，本文主要利用上述4顆IRP開展自主定軌算法研究。

表1 8個候選觀測對象參數(shù)列表

2 XPNAV-1衛(wèi)星自主定軌數(shù)學(xué)模型

XPNAV-1衛(wèi)星軌道為太陽同步軌道，在軌道可見性分析過程中要求脈沖星方向與太陽夾角大于45°，與月球夾角大于5°，并且要考慮地球遮擋、規(guī)避南大西洋異常區(qū)和高緯度極區(qū)輻射等約束條件。針對該衛(wèi)星軌道特點，建立軌道力學(xué)模型和自主定軌觀測方程。

2.1 軌道力學(xué)模型

XPNAV-1衛(wèi)星軌道力學(xué)方程可表示為:

(1)

式中:X(t)為狀態(tài)變量；f[X(t),t]為軌道動力學(xué)方程；w(t)為系統(tǒng)噪聲項。

在該衛(wèi)星自主定軌研究中，采用只考慮地球中心引力加速度和J2項非球形攝動加速度的衛(wèi)星軌道力學(xué)模型，代入式(1)得到衛(wèi)星軌道力學(xué)方程:

(2)

式中:r=[rx,ry,rz]T為衛(wèi)星位置矢量；v=[vx,vy,vz]T為衛(wèi)星速度矢量；X=[rx,ry,rz,vx,vy,vz]T為衛(wèi)星狀態(tài)矢量；w=[wx,wy,wz,wvx,wvy,wvz]T為系統(tǒng)噪聲；μ為地球引力常數(shù)；J2為帶諧項系數(shù)；Re為地球半徑;r為衛(wèi)星到地心的距離。

2.2 自主定軌觀測方程

在太陽系質(zhì)心(solar system barycenter,SSB)坐標(biāo)系下，根據(jù)脈沖星導(dǎo)航的基本原理，可以得到脈沖星導(dǎo)航的基本測量方程[16]：

(3)

式中：tSSB和tsat為脈沖星同一脈沖到達SSB和衛(wèi)星的時間；n為觀測脈沖星方向向量；RSC為衛(wèi)星相對于SSB的位置矢量；c為光速；δt為星載原子時鐘偏差。

相較于脈沖星測距精度和觀測數(shù)目等因素，衛(wèi)星鐘差對定軌精度的影響較小。不失一般性，針對本文開展的自主定軌研究，衛(wèi)星鐘差忽略不計，并考慮地球、衛(wèi)星和SSB三者位置關(guān)系，將式(3)改寫如下：

c·(tSSB-tsat)-n·RE=n·R

式中：RE為地球相對于SSB的位置矢量；R為衛(wèi)星相對于地球質(zhì)心的位置矢量。

在tk時刻衛(wèi)星觀測1顆脈沖星時，測量方程可以表示為：

Zk=HkXk+Vk

(4)

在X射線脈沖星定軌中決定定軌精度的主要因素是脈沖到達時間(TOA)的精度，TOA觀測量的測量噪聲方程可以估算為[17]：

(5)

式中：FB為X射線背景輻射流強；Fx為X射線輻射流強；tobs為信號觀測時間；P為脈沖周期；W為脈沖寬度；pf為脈沖調(diào)制度；Ae為探測器有效面積。

3 EKF算法及PWCS分析方法

標(biāo)準(zhǔn)卡爾曼濾波主要針對隨機離散線性系統(tǒng)模型，針對軌道動力學(xué)非線性的特性，采用EKF來解決隨機連續(xù)非線性系統(tǒng)問題。并在對系統(tǒng)進行濾波處理之前，用PWCS方法定性地分析系統(tǒng)狀態(tài)的可觀測性，這是反映濾波收斂精度和速度的重要指標(biāo)。

3.1 EKF導(dǎo)航濾波算法

EKF的基本思想是將導(dǎo)航系統(tǒng)的狀態(tài)方程進行泰勒展開，忽視二階及二階以上的高階項，將展開式的一階線性項作為原狀態(tài)方程的近似表達式，再采用線性問題的標(biāo)準(zhǔn)卡爾曼濾波方法對近似的狀態(tài)量進行估計，從而得到最終狀態(tài)的次優(yōu)近似值。

對式(2)和式(4)相應(yīng)簡化后可以得到如下系統(tǒng)模型：

(6)

式中：t為時間參數(shù)；w(t)和v(t)分別為系統(tǒng)噪聲和觀測噪聲，且彼此不相關(guān)。

式中：δX(t)為兩者的狀態(tài)估計誤差。

進而得到誤差狀態(tài)的動力學(xué)方程為：

式中：F(t)為雅可比矩陣。進一步將狀態(tài)方程進行離散化處理，得到：

δXk=Φk,k-1δXk-1+Wk-1

(7)

式中：Φk,k-1=I+F(tk-1)·T，為k-1時刻到k時刻的狀態(tài)轉(zhuǎn)移陣；T為濾波周期。

由式(4)和式(7)可以建立系統(tǒng)的離散型線性干擾方程：

擴展卡爾曼濾波流程如下：

(1)濾波初始化

根據(jù)估計狀態(tài)信息，給出狀態(tài)變量及其相應(yīng)誤差方差陣初始值：

(2)狀態(tài)更新

(3)量測更新

Pk=(I-KkHk)Pk/k-1(I-KkHk)T+

3.2 PWCS可觀測性分析方法

不論采用卡爾曼濾波或其他濾波方式，系統(tǒng)狀態(tài)變量的可觀測性是濾波收斂的前提條件。針對X射線脈沖星自主定軌系統(tǒng)是典型的非線性時變系統(tǒng)，提出基于分段式線性定常系統(tǒng)(PWCS)的可觀測度分析方法。在足夠小的時間區(qū)間內(nèi)，可以忽略線性時變系統(tǒng)的系數(shù)矩陣變化，在該時間區(qū)間內(nèi)就可以把時變系統(tǒng)近似為一個分段線性定常系統(tǒng)進行處理，這往往可以大大簡化分析過程。

離散型PWCS系統(tǒng)的一般形式為：

式中：j=1,2,…,q，為時變系統(tǒng)分段間隔序號。于是，系統(tǒng)總的可觀測性矩陣(TOM)可以表示為：

式中：

根據(jù)PWCS分析定理，如果系統(tǒng)滿足

null(Qj)?null(Φj), 1≤j≤q

則有

因此，可以用QS來代替QT分析系統(tǒng)的可觀測性，使問題得到簡化。又因系統(tǒng)狀態(tài)呈遞推關(guān)系，每個時間段的系統(tǒng)狀態(tài)初值就是前一時間段的系統(tǒng)狀態(tài)終值，則可觀測性矩陣仍具有累積繼承的性質(zhì)。這樣可以用某時段j的可觀測性矩陣Qj進一步代替QS，使系統(tǒng)的可觀測性分析進一步簡化。

系統(tǒng)的可觀測矩陣Qj可以表示為：

(8)

式中：

(9)

S中每個元素的量級為μ/r3，對于軌道高度為500 km的衛(wèi)星，其數(shù)量級為10-6，近似為零，將式(9)其代入式(8)，并進行初等行變換，得到：

通過提取可觀測矩陣，進一步進行可觀測度的計算，用可觀測矩陣的條件數(shù)來對可觀測性進行量化，對Qj進行奇異值分解：

Qj=UΛVT

定義可觀測矩陣的條件數(shù)如下：

式中：σmax和σmin分別為M的最大奇異值和最小奇異值。矩陣的條件數(shù)η越大，說明矩陣越接近奇異。

系統(tǒng)可觀測性矩陣的條件數(shù)能夠反映系統(tǒng)狀態(tài)的可觀測程度，即系統(tǒng)可觀測性矩陣的條件數(shù)越大，系統(tǒng)的可觀測度越差。如果系統(tǒng)可觀測性矩陣的條件數(shù)是無窮大，則系統(tǒng)不可觀測；如果系統(tǒng)可觀測性矩陣的條件數(shù)等于1，系統(tǒng)的可觀測性最好。

4 數(shù)值分析試驗

4.1 XPNAV-1衛(wèi)星的觀測數(shù)據(jù)分析

目前，XPNAV-1衛(wèi)星仍正常在軌運行，對Crab脈沖星和3顆低流量IRP進行了詳細觀測，積累了大量觀測數(shù)據(jù)。從實測數(shù)據(jù)分析來看，掠入射式聚焦型探測器觀測的數(shù)據(jù)有較高的信噪比，因此在本文定軌算法研究中，主要利用該探測器采集的數(shù)據(jù)進行數(shù)值分析試驗。

自2016年11月18日以來，XPNAV-1衛(wèi)星一直對Crab脈沖星進行重點觀測，除每年4月30日至8月1日期間Crab方向和太陽夾角小于45°，以及考慮與月球夾角小于5°不可觀測外，其余每天進行觀測，每個觀測弧段持續(xù)觀測時間為30～50 min。目前，通過累積觀測獲取了比較完整的Crab脈沖星觀測數(shù)據(jù)，并提取得到精化的脈沖輪廓，通過計時擬合建立了精度為55.14 μs的計時模型。

此外，該衛(wèi)星在2017年2月21日至2017年8月31日，分別對3顆低流量IRP進行了長期觀測，每個軌道周期持續(xù)開展5～15 min低流量脈沖星探測，累積探測數(shù)據(jù)。經(jīng)統(tǒng)計分析得到了3顆低流量IRP的流量統(tǒng)計以及能譜分析結(jié)果。但由于探測器有效面積較小，脈沖星光子流量較低，目前未能提取到低流量脈沖星完整的脈沖輪廓。

XPNAV-1衛(wèi)星對Crab脈沖星的進行了完整的觀測。其中，時間跨度從2016年11月18日至2017年3月26日的觀測效果較好，共179段觀測數(shù)據(jù)，每個觀測弧段平均觀測時間約45 min。對該段實測數(shù)據(jù)做如下處理：1)在0.5～9.0 keV的能帶內(nèi)選擇光子；2)篩除有效觀測時間短于40 min的觀測數(shù)據(jù)；3)篩除通量小于14個光子/s的觀測數(shù)據(jù)。最后從中選取時間跨度為30 d的數(shù)據(jù)用作數(shù)值試驗分析，共10組觀測值，其TOA誤差及相應(yīng)的測距誤差如表2所示，實測數(shù)據(jù)的起止時間為2017年1月5日至2017年2月5日。

但是由于實測數(shù)據(jù)的信號積分時間較長，觀測量不充足以及測量序列不規(guī)則,所以僅僅只利用以上Crab的實測數(shù)據(jù)會導(dǎo)致濾波結(jié)果的發(fā)散。因此,基于Crab脈沖星實測數(shù)據(jù)的測距精度來模擬生成額外的觀測值，對30 d內(nèi)的實測數(shù)據(jù)進行補充，對觀測周期進行壓縮，每500 s生成一個觀測值。

考慮到單顆脈沖星定軌系統(tǒng)的可觀測性差，再結(jié)合B0540-69，B1509-58，J1846-0258這3顆低流量IRP，利用1顆、2顆、3顆、4顆脈沖星的觀測數(shù)據(jù)組成數(shù)據(jù)組，分別進行軌道確定，分析觀測脈沖星數(shù)目對定軌精度的影響。利用Crab脈沖星實測數(shù)據(jù)，由Crab脈沖星觀測TOA與Crab脈沖星模型預(yù)報TOA的差，可以得到每個觀測時間段的TOA誤差(見表2)，計算得到TOA均方根誤差為99.05 μs，對應(yīng)29.70 km的測距均方根誤差。由于3顆低流脈沖星光子流量較低，探測器有效面積較小，未能得到其完整的脈沖輪廓。由Crab脈沖星的完整觀測數(shù)據(jù)反算出背景噪聲，再根據(jù)3顆低流量脈沖星的實測光子流量統(tǒng)計以及探測器實際有效面積，利用式(5)的估計方程推算得到其測距精度，并由測距精度模擬生成觀測數(shù)據(jù)。表3給出了4顆脈沖星的赤經(jīng)、赤緯和測距精度。

表2 基于Crab脈沖星觀測結(jié)果的TOA誤差及測距誤差

表3 4顆IRP的測距精度估計

4.2 顆脈沖星的可見性分析

在XPNAV-1衛(wèi)星軌道覆蓋性分析過程中，影響其對目標(biāo)脈沖星觀測的約束條件是：1)地球處于觀測目標(biāo)與衛(wèi)星之間，遮擋觀測信號；2)觀測目標(biāo)與太陽之間的夾角小于45°，或者觀測目標(biāo)與月球之間的夾角小于5°，高能粒子對探測器造成影響。XPNAV-1衛(wèi)星的軌道周期約為94 min，圖1和圖2給出了4顆脈沖星在一個軌道周期內(nèi)的可見性狀況。

圖1 4顆脈沖星的可見性分析Fig.1 Visibility analysis of four pulsars

圖2 一個軌道周期內(nèi)可見星個數(shù)Fig.2 Number of visible pulsars in an orbit period

由圖2可以看出，一個軌道周期內(nèi)的任意時刻均最少有2顆脈沖星可被觀測到，脈沖星可見的時間百分比為100%，依靠4顆目標(biāo)脈沖星的觀測信息，可以實現(xiàn)全時段定軌，從每個觀測弧段的可見星中選取單星進行導(dǎo)航，可以保證觀測信息的完備性。

4.3 濾波系統(tǒng)狀態(tài)可觀測性分析

系統(tǒng)可觀測度的數(shù)值試驗結(jié)果如圖3所示，縱坐標(biāo)表示可觀測矩陣條件數(shù)，且取底數(shù)為10的對數(shù)。

表4給出了觀測不同脈沖星數(shù)目時，系統(tǒng)可觀測矩陣的條件數(shù)。當(dāng)全程只觀測1顆或2顆脈沖星時，可觀測矩陣Qj條件數(shù)的數(shù)量級較大，矩陣接近奇異，定軌系統(tǒng)可觀測性差，無法實現(xiàn)濾波過程收斂。當(dāng)觀測3顆脈沖星時，可觀測矩陣條件數(shù)明顯變小，系統(tǒng)可觀測度提高。當(dāng)觀測4顆脈沖星時，可觀測矩陣Qj的條件數(shù)均值約為57，系統(tǒng)可觀測度較好，明顯優(yōu)于前3組結(jié)果?？梢钥偨Y(jié)出當(dāng)觀測脈沖星數(shù)目增加時，可觀測矩陣條件數(shù)變小，系統(tǒng)可觀測度提高，數(shù)值試驗結(jié)果與上節(jié)理論分析結(jié)果一致。

圖3 基于4顆脈沖星的可觀測矩陣條件數(shù)Fig.3 Conditional number of observable matrix based on observing four pulsars

表4 脈沖星數(shù)量對系統(tǒng)可觀測度影響

4.4 自主定軌結(jié)果分析

利用第3組脈沖星數(shù)據(jù)組進行定軌計算，得到軌道誤差和速度誤差序列分別如圖4和圖5所示，進一步統(tǒng)計分析得到位置誤差均方根為56.65 km，速度誤差均方根為0.055 1 km/s。利用第4組脈沖星數(shù)據(jù)進行定軌計算的軌道誤差和速度誤差序列分別如圖6和圖7所示，進一步統(tǒng)計分析得到位置誤差均方根為40.24 km，速度誤差均方根為0.043 3 km/s。

圖4 基于觀測3顆脈沖星的軌道誤差Fig.4 Position error based on observing three pulsars

圖5 基于觀測3顆脈沖星的速度誤差Fig.5 Speed error based on observing three pulsars

圖6 基于觀測4顆脈沖星的軌道誤差Fig.6 Position error based on observing four pulsars

圖7 基于觀測4顆脈沖星的速度誤差Fig.7 Speed error based on observing four pulsars

由圖4～圖7可見，當(dāng)觀測3顆或4顆脈沖星時，濾波過程收斂，并且隨著觀測脈沖星數(shù)目的增加，軌道誤差減小，定軌精度提高。

5 結(jié)束語

針對XPNAV-1衛(wèi)星拓展試驗任務(wù)及脈沖星導(dǎo)航后續(xù)發(fā)展需求，利用多顆脈沖星實測數(shù)據(jù)，對EKF自主定軌算法進行數(shù)值試驗驗證。研究結(jié)果表明：當(dāng)僅觀測1顆脈沖星時，系統(tǒng)的可觀測矩陣條件數(shù)極大，濾波結(jié)果很快發(fā)散；當(dāng)觀測2顆脈沖星時，系統(tǒng)可觀測度提高，但條件數(shù)依然較大，濾波結(jié)果經(jīng)過長時間依舊發(fā)散；當(dāng)觀測3顆脈沖星時，系統(tǒng)可觀測度明顯提高，脈沖星對軌道完全覆蓋，濾波過程收斂，軌道確定精度為56.65 km；當(dāng)觀測4顆脈沖星時，系統(tǒng)可觀測度進一步提高，脈沖星對軌道覆蓋性能好，濾波過程收斂，軌道確定精度為40.24 km；在所選脈沖星滿足該衛(wèi)星整個軌道周期覆蓋的條件下，每個弧段僅需要觀測1顆脈沖星就可以進行定軌計算，其定軌精度取決于所選脈沖星的空間分布和數(shù)量。本文提出的EKF自主定軌算法過程收斂，試驗結(jié)果驗證了算法的合理性和有效性。EKF算法適用于線性高斯白噪聲的系統(tǒng)狀態(tài)估計，而實際的航天器軌道確定數(shù)學(xué)模型往往是非線性或非高斯的，存在線性化誤差、建模誤差和有色噪聲等問題。下一步將利用脈沖星實測數(shù)據(jù)，開展基于UKF、H∞濾波和粒子濾波等算法研究，對比分析各種算法性能，提出最優(yōu)的定軌算法解決方案，為脈沖星導(dǎo)航算法研究和后續(xù)空間飛行試驗提供技術(shù)儲備。