姜 宇，陳 曉，王龍奇，金 晶，馬家辰

（1.哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱 150001；2.上海衛(wèi)星工程研究所，上海 201109）

0 引言

脈沖星是一種高速自轉(zhuǎn)的強(qiáng)磁化中子星，具有自轉(zhuǎn)周期穩(wěn)定、輻射光子能量集中的特點(diǎn)。它能夠?yàn)榻剀壍?、深空和星際空間飛行的航天器提供位置、速度、時(shí)間和姿態(tài)等導(dǎo)航信息，實(shí)現(xiàn)對(duì)長(zhǎng)時(shí)間在軌航天器的高精度自主導(dǎo)航與控制，具有廣闊的工程應(yīng)用前景［1-2］。

脈沖星的自轉(zhuǎn)軸與磁軸存在固定夾角，又因?yàn)槊}沖星的自轉(zhuǎn)周期非常穩(wěn)定，航天器接收到周期性的光子信號(hào)。

記錄航天器單位時(shí)間內(nèi)接收的光子數(shù)可以獲得脈沖星觀測(cè)輪廓，對(duì)比脈沖星的觀測(cè)輪廓和太陽系質(zhì)心（Solar System Barycenter，SSB）的標(biāo)準(zhǔn)輪廓可以得到光子到達(dá)航天器與SSB 的時(shí)間差。把時(shí)間差作為觀測(cè)量建立三維觀測(cè)模型，結(jié)合軌道動(dòng)力學(xué)模型，利用相關(guān)濾波算法可以得到航天器的位置和速度信息［3-7］。

隨著空間探測(cè)向大尺度時(shí)空發(fā)展，計(jì)時(shí)精度向納秒量級(jí)提高，三維歐式空間下的脈沖星光子傳播模型已經(jīng)不能滿足精度要求。光子在傳播過程中，不僅需要考慮其周年視差效應(yīng)、多普勒頻移以及等離子體的色散延緩，還要考慮太陽系引力的彎曲效應(yīng)［8-10］。

文獻(xiàn)［11］基于DSX 體系下的后牛頓近似方法，將SSB 坐標(biāo)系作為航天器導(dǎo)航的參考系，推導(dǎo)了1PN、2PN 度規(guī)下的高精度觀測(cè)模型，可以滿足納秒量級(jí)的時(shí)間測(cè)量精度。

文獻(xiàn)［12-13］提出在四維光坐標(biāo)框架下討論脈沖星光子傳播問題，將4 個(gè)脈沖星的固有時(shí)作為四維平直時(shí)空的光坐標(biāo)，可以為太陽系內(nèi)的任意觀測(cè)者導(dǎo)航，由此建立相對(duì)論定位系統(tǒng)理論框架。

文獻(xiàn)［14］在四維黎曼時(shí)空中引入相對(duì)論定位系統(tǒng)，建立四維觀測(cè)模型，與傳統(tǒng)的三維觀測(cè)模型相比引入時(shí)間狀態(tài)量，提高了時(shí)間估計(jì)的準(zhǔn)確性，且不必推算標(biāo)準(zhǔn)輪廓的相位，減小了標(biāo)準(zhǔn)輪廓的擬合誤差。

X 射線脈沖星自主導(dǎo)航依賴于精確的時(shí)間測(cè)量，航天器時(shí)鐘記錄的固有時(shí)與SSB 坐標(biāo)系下的坐標(biāo)時(shí)不同步，在500 s 的觀測(cè)窗口下，固有時(shí)和坐標(biāo)時(shí)會(huì)產(chǎn)生微妙量級(jí)的誤差，代入迭代過程會(huì)導(dǎo)致濾波器發(fā)散，所以需要引入時(shí)間轉(zhuǎn)換模型。利用光子運(yùn)動(dòng)的四維時(shí)空間隔為零，可以得到固有時(shí)與坐標(biāo)時(shí)的基本關(guān)系。

文獻(xiàn)［15-16］研究環(huán)火星探測(cè)器星載時(shí)鐘在1PN 度規(guī)下的時(shí)間轉(zhuǎn)換關(guān)系，通過數(shù)值分析得出時(shí)間不同步主要受太陽引力和航天器速度影響，可達(dá)亞秒/年量級(jí)。上述關(guān)系式以坐標(biāo)時(shí)作為自變量進(jìn)行積分運(yùn)算，并且積分式中涉及航天器的位置、速度信息，很難實(shí)現(xiàn)從固有時(shí)到坐標(biāo)時(shí)的轉(zhuǎn)換。

文獻(xiàn)［17］根據(jù)廣義相對(duì)論的后牛頓近似理論，在忽略太陽系天體自轉(zhuǎn)和扁率的情況下，利用能量守恒定律推導(dǎo)出航天器繞飛的時(shí)間轉(zhuǎn)換模型，模型精度在納秒量級(jí)。文獻(xiàn)［18］針對(duì)航天器繞地變軌時(shí)可能做拋物線或雙曲線運(yùn)動(dòng)推導(dǎo)了更加全面的時(shí)間轉(zhuǎn)換模型，對(duì)文獻(xiàn)［19］進(jìn)行了補(bǔ)充分析。

本文基于繞飛軌道研究高精度脈沖星導(dǎo)航模型，為脈沖星導(dǎo)航的深空探測(cè)提供理論基礎(chǔ)，所做工作集中在以下兩方面：在四維黎曼時(shí)空下建立脈沖星導(dǎo)航系統(tǒng)狀態(tài)模型和四維觀測(cè)模型，并基于無跡卡爾曼濾波（UKF）方法設(shè)計(jì)脈沖星導(dǎo)航算法；對(duì)引入的四維觀測(cè)模型和時(shí)間轉(zhuǎn)換模型進(jìn)行誤差分析，研究影響導(dǎo)航精度的主要因素。

1 四維黎曼時(shí)空下的X 射線脈沖星導(dǎo)航模型

1.1 航天器軌道力學(xué)模型

選取太陽系質(zhì)心天球參考系（Barycentric Celestial Reference System，BCRS）的時(shí)空度規(guī)，建立帶有時(shí)間轉(zhuǎn)換項(xiàng)的航天器軌道運(yùn)動(dòng)學(xué)模型。X 射線脈沖星導(dǎo)航原理圖如圖1 所示。

圖1 X 射線脈沖星導(dǎo)航原理圖Fig.1 Schematic diagram of X ray pulsar-based navigation

四維黎曼時(shí)空下航天器的狀態(tài)x定義為

式中：t為航天器的坐標(biāo)時(shí)；rsa=rse+rea、vsa=vse+vea分別為航天器相對(duì)于SSB 的位置和速度矢量；rse、vse分別為繞飛行星相對(duì)于SSB 的位置和速度矢量 ；rea=[rea，x rea，y rea，z]T、vea=[vea，x vea，y vea，z]T分別為航天器相對(duì)于行星的位置和速度矢量。

相對(duì)繞飛行星建立航天器軌道運(yùn)動(dòng)學(xué)方程為［5］

式中：Re為繞飛行星的半徑；μe為繞飛行星的引力常數(shù)；J2為攝動(dòng)項(xiàng)系數(shù)；ΔF=[ΔFx，ΔFy，ΔFz]T為攝動(dòng)項(xiàng)矩陣，指的是包括月球引力攝動(dòng)、潮汐攝動(dòng)在內(nèi)的其他攝動(dòng)力影響。將式（2）簡(jiǎn)寫為

在四維黎曼時(shí)空下，航天器的坐標(biāo)時(shí)可由時(shí)鐘記錄的固有時(shí)轉(zhuǎn)換得到。在1PN 度規(guī)下，忽略星體的自轉(zhuǎn)、扁率以及太陽與SSB 間距，可得到從航天器固有時(shí)至SSB 滿足納秒量級(jí)精度的時(shí)間轉(zhuǎn)換模型［16］：

式中：Δτ為星載時(shí)鐘走過的固有時(shí)；Δt為Δτ時(shí)間內(nèi)時(shí)鐘經(jīng)歷的坐標(biāo)時(shí)；c為光速；等號(hào)后的第1 項(xiàng)為Δt時(shí)間內(nèi)運(yùn)動(dòng)產(chǎn)生的綜合時(shí)延；為引力修正因子；μs為太陽引力常數(shù)；μe為行星的引力常數(shù)；ase為行星繞太陽作橢圓運(yùn)動(dòng)的長(zhǎng)半軸；aea為航天器繞行星運(yùn)動(dòng)的長(zhǎng)半軸。

將式（4）簡(jiǎn)寫為

式中：Δt0=LΔτ包含了行星引力產(chǎn)生的頻移，稱為引力時(shí)延；包含了航天器運(yùn)動(dòng)以及行星繞太陽運(yùn)動(dòng)偏離圓軌道所產(chǎn)生的頻移，稱為偏軌時(shí)延；包含了航天器軌道曲率變化所產(chǎn)生的頻移，稱為航天器曲率時(shí)延；包含了行星軌道曲率變化所產(chǎn)生的頻移，稱為行星曲率時(shí)延；Δt0為長(zhǎng)期影響項(xiàng)；Δt1、Δt2、Δt3為周期影響項(xiàng)。

根據(jù)式（1）～式（4），航天器的狀態(tài)模型最終描述為

式中：f(x，τ)=；τ為星載時(shí)鐘的固有時(shí)刻，這是實(shí)際飛行任務(wù)中可獲取的時(shí)間狀態(tài)量；w(τ)為系統(tǒng)的傳遞噪聲矩陣。

1.2 四維觀測(cè)模型

通常意義下的時(shí)間概念是太陽系質(zhì)心坐標(biāo)系下的航天器坐標(biāo)時(shí)，而實(shí)際中可獲取的時(shí)間變量是星載時(shí)鐘原子時(shí)，它們之間存在一種轉(zhuǎn)換關(guān)系。在高精度要求的情況下，坐標(biāo)時(shí)成為了一個(gè)關(guān)于航天器狀態(tài)和原子時(shí)變化的新觀測(cè)量，同時(shí)單顆脈沖星僅可以在組合導(dǎo)航中應(yīng)用，或者作投影方向上的導(dǎo)航算法精度驗(yàn)證，不能實(shí)現(xiàn)自主導(dǎo)航。因此，需要引入四維觀測(cè)模型。

根據(jù)相對(duì)論定位系統(tǒng)的基本思想和后牛頓引力理論，將脈沖星P（P=1，2，3，4）走過的固有時(shí)TP作為直接觀測(cè)量［13］，建立四維觀測(cè)模型如下：

根據(jù)式（7）和式（8），脈沖星導(dǎo)航的四維觀測(cè)模型簡(jiǎn)寫為

四維觀測(cè)模型相比傳統(tǒng)的三維觀測(cè)模型，引入了時(shí)間狀態(tài)量，且不必推算標(biāo)準(zhǔn)輪廓的相位，具有更精確的時(shí)間估計(jì)，更小的相位擬合誤差。

2 基于UKF 的脈沖星導(dǎo)航算法

根據(jù)式（6）和式（9）得到脈沖星導(dǎo)航系統(tǒng)的離散模型為

式中：wk、vk的協(xié)方差矩陣分別為Qk和Rk，且滿足

利用UKF 算法進(jìn)行脈沖星導(dǎo)航濾波的步驟如下：

步驟1時(shí)間轉(zhuǎn)換。上一時(shí)刻到當(dāng)前時(shí)刻走過的固有時(shí)為Δτk?1，根據(jù)式（4）可將航天器記錄的固有時(shí)轉(zhuǎn)換為坐標(biāo)時(shí)，即

步驟2構(gòu)造采樣點(diǎn)?？柭鼮V波傳遞狀態(tài)量的一階矩和二階矩信息，對(duì)于狀態(tài)均值及其協(xié)方差矩陣Pk?1，選取比例對(duì)稱采樣策略對(duì)其進(jìn)行采樣，但在采樣前需要保證算法的穩(wěn)定性。由于協(xié)方差矩陣Pk?1存在時(shí)間項(xiàng)的極小對(duì)角值，在迭代過程中不能保證其對(duì)稱半正定，需對(duì)Pk?1作對(duì)稱正定化處理：

式中：α為比例縮放因子，可通過調(diào)整α取值來調(diào)節(jié)采樣點(diǎn)與中心點(diǎn)的距離；β為函數(shù)f(?)的高階項(xiàng)參數(shù)；κ涉及高 階矩的選取，一般 取為0 或3?n，但應(yīng)確保協(xié)方差矩陣傳遞的半正定性。

式中：(·)i算子為取矩陣的第i列信息。

將采樣點(diǎn)展開為

步驟3預(yù)測(cè)。利用狀態(tài)模型對(duì)每個(gè)采樣點(diǎn)進(jìn)行預(yù)測(cè)，得到

計(jì)算狀態(tài)均值及其協(xié)方差矩陣為

計(jì)算觀測(cè)量及其均值為

步驟4更新。觀測(cè)量與狀態(tài)量之間的預(yù)測(cè)協(xié)方差陣分別為

濾波器增益為

更新時(shí)刻的狀態(tài)均值及其協(xié)方差矩陣為

3 仿真分析

3.1 仿真條件

仿真采用DE405 星歷，初始時(shí)間為1990 年7 月7 日12 時(shí)0 分0 秒，時(shí)長(zhǎng)為35 d，觀測(cè)時(shí)間間隔為500 s。

本文的導(dǎo)航模型適用于太陽系任一行星的繞飛軌道，這里以繞飛地球高軌為例，其軌道參數(shù)見表1。

表1 軌道參數(shù)Tab.1 Orbit parameters

選用Arecibo 天文臺(tái)長(zhǎng)期觀測(cè)的4 顆毫秒級(jí)脈沖星，其參數(shù)見表2。

表2 脈沖星參數(shù)Tab.2 Pulsar parameters

計(jì)算機(jī)仿真環(huán)境為CPU 4 核2.5 GHz，內(nèi)存8 GB，Matlab 2012b。初始狀態(tài)誤差δx0、初始協(xié)方差陣P0、系統(tǒng)噪聲協(xié)方差陣Qk、觀測(cè)噪聲協(xié)方差陣Rk設(shè)置如下：

3.2 仿真結(jié)果

UKF 濾波后的導(dǎo)航結(jié)果如圖2所示，導(dǎo)航的位置誤差穩(wěn)定在10 m 量級(jí)，時(shí)間誤差穩(wěn)定在nm 量級(jí)，速度誤差穩(wěn)定在0.01 m/s 量級(jí)。下面分別討論時(shí)間轉(zhuǎn)換模型和四維觀測(cè)模型中的時(shí)延項(xiàng)對(duì)導(dǎo)航精度的影響。

圖2 四維黎曼時(shí)空下脈沖星導(dǎo)航模型的估計(jì)誤差Fig.2 Estimation errors of the pulsar-based navigation model in 4D Riemann space-time

3.2.1 時(shí)間轉(zhuǎn)換模型精度分析

時(shí)間轉(zhuǎn)換模型的各個(gè)時(shí)延項(xiàng)在迭代中所占的比重如圖3 所示。其中，引力時(shí)延Δt0數(shù)值最大，達(dá)10?5s 量級(jí)；偏軌時(shí)延Δt1、行星曲率時(shí)延Δt3比重相對(duì)較小，為10?7s 量級(jí)；航天器曲率時(shí)延Δt2數(shù)值最小，為10?10s 量級(jí)。忽略時(shí)間轉(zhuǎn)換模型中的任一時(shí)延項(xiàng)進(jìn)行導(dǎo)航，結(jié)果如圖4 所示?？梢钥闯?，忽略Δt0對(duì)導(dǎo)航精度的影響最大，達(dá)105m 量級(jí)；忽略Δt1、Δt3對(duì)導(dǎo)航精度的影響相對(duì)較小，為103m 量級(jí)；忽略Δt2對(duì)導(dǎo)航精度的影響最小，為10 m 量級(jí)。這表明在時(shí)間轉(zhuǎn)換過程中，引力時(shí)延、航天器偏軌時(shí)延以及繞飛行星曲率時(shí)延等因素對(duì)導(dǎo)航精度的影響較大，航天器繞行星飛行的曲率變化對(duì)導(dǎo)航精度的影響較小。上述分析驗(yàn)證了引入時(shí)間轉(zhuǎn)換模型的必要性。

圖3 時(shí)間轉(zhuǎn)換模型中的時(shí)延項(xiàng)比重Fig.3 Proportion of delay items in the time transformation model

圖4 忽略時(shí)間轉(zhuǎn)換模型時(shí)延項(xiàng)的導(dǎo)航定位誤差Fig.4 Positioning errors without a time transformation term

3.2.2 四維觀測(cè)模型精度分析

四維觀測(cè)模型中的各個(gè)時(shí)延項(xiàng)在迭代中所占比重如圖5 所示。圖5 中，引力時(shí)延ρ3(r)數(shù)值較大，達(dá)10?4s 量級(jí)，周年視差時(shí)延ρ1(r)數(shù)值較小，達(dá)10?7s量級(jí)，脈沖星自行時(shí)延ρ2(r)總體上呈增長(zhǎng)趨勢(shì)，在仿真時(shí)間內(nèi)與ρ1(r)達(dá)到同一量級(jí)。若忽略其中一項(xiàng)進(jìn)行導(dǎo)航，結(jié)果如圖6 所示。

圖6 忽略四維觀測(cè)模型時(shí)延項(xiàng)的導(dǎo)航定位誤差Fig.6 Positioning errors without a time delay term of the 4D observation model

可以看出，忽略ρ3(r)對(duì)導(dǎo)航精度的影響很大，達(dá)104m 量級(jí)；忽略ρ1(r)、ρ2(r)對(duì)導(dǎo)航精度的影響較小，為103m級(jí)，并且由于脈沖星自行時(shí)延有累積增長(zhǎng)趨勢(shì)；忽略ρ2(r)引起的誤差按其趨勢(shì)增長(zhǎng)，與圖5的結(jié)果吻合。上述分析表明，在量測(cè)過程中太陽系的引力遠(yuǎn)大于周年視差效應(yīng)、脈沖星自行對(duì)導(dǎo)航精度的影響。

圖5 四維觀測(cè)模型中時(shí)延項(xiàng)比重Fig.5 Proportions of delay items in the 4D observation model

4 結(jié)束語

本文建立的導(dǎo)航模型考慮了廣義相對(duì)論的影響，通過引入時(shí)間轉(zhuǎn)換模型可適應(yīng)ns 量級(jí)的計(jì)時(shí)精度，在此基礎(chǔ)上分析了導(dǎo)航模型的精度影響因素，為深空探測(cè)的工程應(yīng)用提供了理論參考。同時(shí)使用本文建立的導(dǎo)航模型可以簡(jiǎn)化脈沖星信號(hào)處理過程，為脈沖星導(dǎo)航算法提供更精準(zhǔn)的模型基礎(chǔ)。