基于截斷誤差的改進脈沖星導航觀測方程

許強， 王宏力， 何貽洋， 由四海， 馮磊

(火箭軍工程大學導彈工程學院, 西安 710025)

X射線脈沖星是宇宙中高速旋轉并周期性地輻射電磁脈沖的中子星。由于具有高度穩(wěn)定的旋轉周期，使得脈沖星所輻射脈沖的周期及輪廓具有較強的不變性和較高的可識別性。X射線脈沖星導航系統(tǒng)就是利用其穩(wěn)定的X射線脈沖信號作為觀測量進行航天器導航的系統(tǒng)，它具有隱蔽性高、自主性強、抗干擾性好等優(yōu)點。近年來關于X射線脈沖星導航的相關研究也日益深入。美國國防高等研究計劃局早在2004年就已經(jīng)正式提了X射線脈沖星導航計劃——XNAV計劃，并于2005年開始研制[1]，歐洲航空局也在2005年制定了類似的“深空探測器脈沖星導航研究(deep space vessel pulsar navigation study)”計劃。后續(xù)俄羅斯、日本等也都相繼開展了相關的研究項目。中國在2016年11月10日于酒泉衛(wèi)星發(fā)射中心成功發(fā)射了第一顆脈沖星導航試驗衛(wèi)星XPNAV-1[2]。

X射線脈沖星導航的測量原理是根據(jù)航天器與太陽系質心(SSB)在脈沖星方向的多普勒延遲得到航天器在SSB坐標系中的狀態(tài)[3]。觀測量是脈沖到達航天器與SSB的時間差或其與光速的乘積。但在宇宙中，脈沖信號的傳播會受天體運動、相對論效應等的影響而產(chǎn)生變化，也就在觀測方程之中引入了高階項[4]。在當前研究中，擴展卡爾曼濾波(EKF)算法因具有工程實現(xiàn)簡單、算法成熟、實時性好等特點而得到廣泛應用。但是在X射線脈沖星導航中應用EKF算法時卻常常因為線性化需要而忽略觀測模型高階項的影響[5-8]，這實際會給導航結果帶來較大截斷誤差甚至導致發(fā)散。這些高階項的變化大多是慢時變的，短時間內(nèi)可以認為是常值并進行補償，長時間來看卻有著較大的變化波動。對于X射線脈沖星導航在深空探測器、長期在軌衛(wèi)星上的應用影響較大，有必要通過在線估計的方法予以補償。

針對以上問題，本文以地球長期在軌衛(wèi)星為背景，通過對高階項的產(chǎn)生機理及數(shù)值量級進行仿真分析，提出一種改進的線性觀測方程。該方程既不需要額外增加狀態(tài)量，也滿足以EKF為基礎的各種拓展濾波算法應用于脈沖星導航時對觀測方程線性化的要求。并通過對比不同時段X射線脈沖星導航系統(tǒng)的仿真結果,證明了該方程的有效性。

1 典型脈沖星導航空間方程

1.1 狀態(tài)方程

X射線脈沖星導航的狀態(tài)方程是基于航天器在地球質心慣性系中的動力學方程建立起來的。狀態(tài)量選為航天器在地球質心慣性系中的位置和速度。在考慮二階帶諧項時，其狀態(tài)方程為[9]

(1)

則可將其狀態(tài)方程簡寫為

(2)

式中：w(t)為系統(tǒng)噪聲矩陣，且

(3)

1.2 觀測方程

本文X射線脈沖星導航中，觀測量選用X射線脈沖信號到達航天器和SSB的時間差Δt與光速c的乘積。如果不考慮其他物理效應的影響，其簡化的觀測模型如下[10]：

cΔt(i)=nirsat

(4)

式中：cΔt(i)為第i顆脈沖星對應的觀測量;ni為第i顆脈沖星的方向矢量，下標i為脈沖星的編號；rsat為衛(wèi)星在SSB坐標系中的位置矢量。由于狀態(tài)方程中使用的狀態(tài)量是航天器在地球質心慣性坐標系中的位置和速度，所以可將rsat表示為rsat=rE+r，rE為SSB坐標系中地球的位置，可以通過地球星歷預報得到,r為衛(wèi)星在地球質心坐標系的位置矢量。所以，X射線脈沖星的觀測方程可以簡寫為

(5)

2 觀測模型高階項分析

如圖1所示，簡化的線性觀測模型中僅有幾何延遲項，即nrsat=cΔt，n為脈沖星的方向矢量。脈沖信號在宇宙中傳播時理論上能夠對傳播時間產(chǎn)生影響的還有相對論效應、色散延遲等[4]。但是考慮到當前光子到達時間的測量精度在微秒量級，理論精度應達到10-8s，所以能夠產(chǎn)生截斷誤差的因素主要考慮地球周年運動導致的周年視差效應和恒星引力場作用下光線彎曲產(chǎn)生的引力延遲效應[11]。

OS-XSYSZS—太陽質心坐標系; OSSB-XSSBYSSBZSSB—太陽系質心坐標系; b—太陽質心到SSB的矢量； rE—地球在太陽系質心坐標系中的位置。圖1 脈沖星導航基本原理Fig.1 Basic principle of pulsar navigation

2.1 周年視差效應

假設觀測者在地球上對恒星進行觀測，由于地球的公轉，恒星在天球上的位置也會發(fā)生改變。以太陽上觀測的恒星在天球上位置作為其平均位置，從地球上觀測的位置為其實際位置。由于地球圍繞太陽公轉，2個位置會存在一定的偏差，也就是所謂的恒星周年視差。周年視差在日地連線(太陽與地球連線)同星日連線(脈沖星與太陽連線)垂直時達到極大值，如圖2所示。

理論上而言，被觀測恒星距離太陽的距離越遠，周年視差效應會越不明顯，所以脈沖星導航中此項的時間量也是一個很小的數(shù)值。但是在X射線脈沖星導航的解算過程中，該項會直接與光速相乘，而且由幾何原理可知，該項產(chǎn)生的誤差投影到3個坐標方向上之后，各方向誤差的代數(shù)和會變大，如式(6)所示:

(6)

在短時間內(nèi)周年視差效應引起的偏差變化不明顯，但是一年內(nèi)的變化幅值卻比較大。如果X射線脈沖星導航系統(tǒng)長時間的運行，那該項的影響必須予以考慮。關于周年視差效應在X射線脈沖星觀測中的影響的具體推導屬于天文學范疇，在此直接由式(7)給出[11-12]：

圖2 周年視差效應原理Fig.2 Principle of annual parallax effect

(7)

式中：D0為脈沖星距離SSB的距離。

2.2 引力延遲效應

引力延遲效應也稱Shapiro效應，是在太陽系內(nèi)驗證廣義相對論的4個經(jīng)典試驗之一。它是指當射電或者光信號經(jīng)過大質量的天體時，受天體重力場的影響其速度方向會發(fā)生變化，傳播的行程增加，傳播到目的地或往返的時間也會因此增加，如圖3所示[13]。所以如果觀測者同樣在地球上，當?shù)厍蛭挥谛侨者B線上時，引力延遲效應最明顯；當?shù)厍蛭挥谛侨者B線的延長線上時，引力延遲效應相互抵消，影響效果最小。

引力延遲效應對導航結果的影響同周年視差效應相同，在此不再贅述。根據(jù)天文學知識，可以得到太陽系內(nèi)天體引力延遲效應的總和為[14]

(8)

由于太陽系內(nèi)太陽質量最大，產(chǎn)生的引力延遲也就最大，約為微秒量級[15]。而其他天體的引力延遲均不大于10-8s，所以受光子測量精度影響可以不予考慮。若僅考慮太陽引力場的作用，則式(8)可以簡化為

(9)

式中：Msun為太陽質量。

圖3 引力延遲效應原理[13]Fig.3 Principle of Shapiro effect[13]

2.3 數(shù)值分析

通過以上分析，能夠產(chǎn)生截斷誤差的主要為地球運動的周年視差效應和太陽的引力延遲效應。利用STK9.0版本軟件按照表1中的軌道參數(shù)分別仿真產(chǎn)生一天及一年內(nèi)的高精度HPOP衛(wèi)星運行軌道，并以表2[16]中B1821-24脈沖星的參數(shù)為例計算式(7)和式(9)的結果，太陽質量Msun=1.989 1×1030kg，萬有引力常數(shù)G=6.67×10-11m3·kg-1·s-2。一天和一年內(nèi)高階項變化情況結果分別如圖4和圖5所示。

通過分析可以發(fā)現(xiàn)：周年視差效應和太陽引力延遲效應的變化幅度明顯，最大與最小值能相差10倍之多，一年內(nèi)周期變化明顯。最大值最小點計算的結果與第2.1節(jié)及2.2節(jié)的理論分析相吻合。

表1 軌道參數(shù)

表2 B1821-24脈沖星的參數(shù)[16]

注：1 kpc=3.08×1019m。

圖4 一天內(nèi)高階項變化情況Fig.4 Variation of higher order terms in a day

圖5 一年內(nèi)高階項變化情況Fig.5 Variation of higher order terms in a year

3 線性觀測方程的改進

由于濾波周期內(nèi)衛(wèi)星繞地球的運動相對于衛(wèi)星到SSB的距離|rsat|來說非常小，所以結合式(7)，在k時刻可對脈沖星導航中的周年視差效應做如下化簡：

(10)

式中：rsat(k-1)為k-1時刻衛(wèi)星在SSB中的位置；rsat(k)為k時刻衛(wèi)星在SSB中的位置。

引力延遲效應的影響量級最大，同理也可對k時刻的式(9)進行如下簡化處理：

(11)

令

(12)

則將式(11)在k-1時刻對rsat進行一階泰勒展開：

a1ln(a2+a3rsat(k))=a1ln(a2+a3rsat(k-1))+

(13)

式中：Δ為泰勒展開中二階以上高階項。

省略二階以上高階項后可得

Δt2≈a1ln(a2+a3rsat(k-1))-

(14)

綜上所述，結合式(5)、式(10)、式(12)、式(14)，可將改進的帶有高階項的線性觀測方程寫為

(15)

4 仿真分析

為驗證本文所提出的改進觀測方程對導航結果的提高效果，本文利用表1中的數(shù)據(jù)仿真產(chǎn)生衛(wèi)星運行軌道。所用的X射線背景流量BX=0.005 ph·cm-2·s-1，探測器面積A=1 m2，脈沖星導航周期為60 s，導航采用B0531+21、B1821-24及B1937+21三顆脈沖星，具體參數(shù)見表3[16]。地球引力常數(shù)為3.986 004 418×1014N·m2/kg，光速為3×108m/s，重力二階帶諧系數(shù)為0.001 082 63，地球赤道半徑為6.378 137×106m。本文中脈沖星導航的其他初始條件設置如下：

1) 初始誤差為

δx=(1 000 m,1 000 m,1 000 m,20 m/s,20 m/s,

20 m/s)

2)初始誤差協(xié)方差為

P(0)=diag[1 0002m2,1 0002m2,1 0002m2,

202m2/s2,202m2/s2,202m2/s2]

3) 系統(tǒng)噪聲協(xié)方差為

式中:q1=8 cm;q2=0.05 mm/s。

4) 量測噪聲協(xié)方差為

Rk=diag[1092m2,3252m2,3442m2]

通過第2節(jié)對高階項的分析可以得知，假如

表3 導航用脈沖星參數(shù)[16]

從衛(wèi)星上對脈沖星進行觀測，在一年中有4個點比較特殊，都是在日衛(wèi)連線(太陽與衛(wèi)星連線)與星日連線相共線或垂直時。由于導航的截斷誤差同時受3顆脈沖星的高階項變化影響，所以3顆脈沖星高階項總和的周年變化在脈沖星導航中是有實際影響意義的。

3顆脈沖星高階項總和的周年變化情況如圖6所示。為了驗證改進的線性觀測方程的有效性，本文選取其最大值和最小值點，即2015年10月17日與2015年12月29日進行仿真運算,仿真結果見圖7～圖10。同時為了更好地評價算法的估計性能，本文采用均方根誤差(RMSE)作為導航誤差的計算公式[17]，統(tǒng)計區(qū)間為一天中的(900,1 440) min區(qū)間，RMSE的表達式為

(16)

式中：Δrk為第k時刻真實的軌道位置與濾波器估計位置之間的距離。

圖6 全年3顆脈沖星高階項總和變化Fig.6 Whole year variation of summation of three pulsars’ higher order terms

圖7 2015年10月17日簡化線性觀測方程估計結果Fig.7 Estimation results of simplified linear measurement equation on October 17, 2015

圖8 2015年10月17日改進線性觀測方程估計結果Fig.8 Estimation results of improved linear measurement equation on October 17, 2015

圖9 2015年12月29日簡化線性觀測方程估計結果Fig.9 Estimation results of simplified linear measurement equation on December 29, 2015

圖10 2015年12月29日改進線性觀測方程估計結果Fig.10 Estimation results of improved linear measurement equation on December 29, 2015

為了降低隨機因素對結果的影響，將2個線性觀測方程分別獨立運算50次并對導航結果求平均值，結果信息統(tǒng)計如表4所示。通過運行結果的圖表來看，在考慮高階項的情況下使用簡化線性觀測方程進行EKF解算時，忽略高階項影響而產(chǎn)生的較大截斷誤差已經(jīng)導致濾波發(fā)散，而改進的線性觀測方程卻仍然能夠較好較快地將位置估計誤差收斂到250 m以內(nèi)，速度估計誤差2 m/s以內(nèi)，并且在高階項的最大值點與最小值點的位置估計誤差僅相差2.6 m，速度誤差僅相差0.038 2 m/s。方程對高階項的變化表現(xiàn)出較好的魯棒性。通過運算時間反應運算量，在60 s更新周期的條件下，用配置i3處理器的臺式計算機運行2012A版本MATLAB進行導航解算，2種方程進行一天內(nèi)的導航解算，總運行時間分別為2.054 4 s和1.749 0 s,平均每次計算僅相差0.000 2 s。

表4 2種觀測方程的估計誤差

5 結 論

本文提出一種適用于脈沖星導航的改進線性觀測方程。該方程包含了脈沖星觀測模型中會對觀測量造成較大高階截斷誤差的周年視差效應和引力延遲效應。仿真結果表明：

1) 在考慮周年視差效應和引力延遲效應影響的情況下，簡化的線性觀測方程在EKF算法中會產(chǎn)生發(fā)散，而改進的線性觀測方程卻能較快的得到250 m和2 m/s以內(nèi)的位置和速度估計誤差。

2) 改進的線性觀測方程對高階項的變化具有一定的魯棒性。在高階項變化的最大值和最小值點位置估計誤差僅相差2.6 m，速度估計誤差僅相差0.038 2 m/s，適合長時在軌航天器的導航解算。