許強(qiáng)， 王宏力， 何貽洋， 由四海， 馮磊

(火箭軍工程大學(xué)導(dǎo)彈工程學(xué)院, 西安 710025)

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

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

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

1 典型脈沖星導(dǎo)航空間方程

1.1 狀態(tài)方程

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

(1)

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

(2)

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

(3)

1.2 觀測方程

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

cΔt(i)=nirsat

(4)

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

(5)

2 觀測模型高階項(xiàng)分析

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

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

2.1 周年視差效應(yīng)

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

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

(6)

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

圖2 周年視差效應(yīng)原理Fig.2 Principle of annual parallax effect

(7)

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

2.2 引力延遲效應(yīng)

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

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

(8)

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

(9)

式中：Msun為太陽質(zhì)量。

圖3 引力延遲效應(yīng)原理[13]Fig.3 Principle of Shapiro effect[13]

2.3 數(shù)值分析

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

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

表1 軌道參數(shù)

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

注：1 kpc=3.08×1019m。

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

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

3 線性觀測方程的改進(jìn)

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

(10)

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

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

(11)

令

(12)

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

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

(13)

式中：Δ為泰勒展開中二階以上高階項(xiàng)。

省略二階以上高階項(xiàng)后可得

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

(14)

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

(15)

4 仿真分析

為驗(yàn)證本文所提出的改進(jìn)觀測方程對導(dǎo)航結(jié)果的提高效果，本文利用表1中的數(shù)據(jù)仿真產(chǎn)生衛(wèi)星運(yùn)行軌道。所用的X射線背景流量BX=0.005 ph·cm-2·s-1，探測器面積A=1 m2，脈沖星導(dǎo)航周期為60 s，導(dǎo)航采用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。本文中脈沖星導(dǎo)航的其他初始條件設(shè)置如下：

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é)對高階項(xiàng)的分析可以得知，假如

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

從衛(wèi)星上對脈沖星進(jìn)行觀測，在一年中有4個(gè)點(diǎn)比較特殊，都是在日衛(wèi)連線(太陽與衛(wèi)星連線)與星日連線相共線或垂直時(shí)。由于導(dǎo)航的截?cái)嗾`差同時(shí)受3顆脈沖星的高階項(xiàng)變化影響，所以3顆脈沖星高階項(xiàng)總和的周年變化在脈沖星導(dǎo)航中是有實(shí)際影響意義的。

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

(16)

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

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

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

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

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

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

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

表4 2種觀測方程的估計(jì)誤差

5 結(jié) 論

本文提出一種適用于脈沖星導(dǎo)航的改進(jìn)線性觀測方程。該方程包含了脈沖星觀測模型中會(huì)對觀測量造成較大高階截?cái)嗾`差的周年視差效應(yīng)和引力延遲效應(yīng)。仿真結(jié)果表明：

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

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