蘇劍宇，方海燕，高敬敬

(西安電子科技大學空間科學與技術(shù)學院，西安 710126)

0 引 言

作為深空網(wǎng)絡(Deep space network，DSN)的補充和替代技術(shù)，X射線脈沖星導航(X-ray pulsar-based navigation，XPNAV)是一種新型的自主導航技術(shù)。XPNAV可為在太陽系內(nèi)外飛行的航天器提供位置、時間等導航信息[1-4]，實現(xiàn)航天器高精度自主導航。目前，隨著X射線探測器技術(shù)的不斷進步，國內(nèi)外相繼開展了X射線脈沖星導航試驗。如中國空間實驗室天宮二號的γ暴偏振探測科學實驗[5]，中國首顆X射線脈沖星導航試驗衛(wèi)星(XPNAV-1)在軌開展的X射線脈沖星的探測與脈沖星導航體制的驗證[6-8]，以及國內(nèi)首顆空間X射線天文硬X射線調(diào)制望遠鏡衛(wèi)星(Hard X-ray Modulation Telescope，HXMT)Insight-HXMT的脈沖星定軌精度驗證實驗[9-10]。國外的如美國NICER(Neutron star interior composition explorer)項目的SEXTANT(Station explorer X-ray timing and naviga-tion technology)搭載在國際空間站(ISS)上開展的定軌精度驗證工作[11-12]。

現(xiàn)有XPNAV導航驗證的主要框架可分為兩種類型：實時和非實時[13]。前者利用脈沖星觀測周期內(nèi)記錄到的一串光子到達時間(Time of arrivals，TOA)，通過一定的算法提取出脈沖相位以及多普勒頻率，由于航天器在任意時刻觀測的脈沖相位和頻率都可以用該時刻的位置、速度以及太陽系質(zhì)心(Solar system barycenter, SSB)處的脈沖星信號模型準確地表示，因此通過對相位與多普勒頻率的解算即可得到相應時刻的位置與速度信息[14-16]。如SEXTANT任務通過對多個毫秒脈沖星的觀測數(shù)據(jù)進行處理，得到10公里以內(nèi)的導航精度[17]。文獻[18]在XPNAV-1的軌道傳播中引入脈沖星視線相對于地球的測距測量作為控制點，通過對脈沖星PSR B0531+21近三個月的觀測，在控制點進行自我定位，平均導航誤差為38.4公里。

根據(jù)以上分析，SEPO方法是一種離線確定軌道根數(shù)的方法，該方法利用觀測輪廓的形變確定軌道根數(shù)，而非脈沖TOA，因此不依賴脈沖星標準輪廓，有助于應用到不同的任務場景中。但是這兩個導航實驗只考慮了單個軌道根數(shù)存在誤差的情況下軌道根數(shù)的反演，實際中軌道根數(shù)可能都存在誤差，此時得到的觀測脈沖輪廓顯著性會產(chǎn)生耦合，無法直接通過觀測脈沖輪廓顯著性的極大值確定軌道根數(shù)的最優(yōu)估計值，如何在觀測輪廓顯著性耦合情況下進行軌道根數(shù)的反演并未進行討論。文獻[13]在SEPO方法的基礎(chǔ)上提出顯著性卡方分組的思想解耦顯著性卡方，并通過對分組后的顯著性卡方再次計算二次卡方確定軌道根數(shù)的最優(yōu)估計值。該方法的顯著性卡方分組的思想為解耦顯著性卡方提供了思路，但該文獻方法仍然使用二次卡方的最大值確定軌道根數(shù)的最優(yōu)估計是不合理的，因為顯著性卡方反映了輪廓的顯著性，在軌道真值附近取極大值，但二次卡方與觀測輪廓并無直接關(guān)系，所以二次卡方最大值并不一定對應軌道根數(shù)的最優(yōu)估計值，這將導致軌道根數(shù)估計值偏離真值；且在逐一估計得到軌道六根數(shù)的過程中，本次的軌道根數(shù)估計將會影響下一軌道根數(shù)的分組，所以本次軌道根數(shù)的較大估計誤差將導致下一軌道根數(shù)產(chǎn)生更大誤差。

綜合以上分析，上述文獻方法雖然可在單個軌道根數(shù)存在偏差時，根據(jù)觀測輪廓顯著性反演軌道根數(shù)，但無法在軌道六根數(shù)均存在偏差時，即觀測輪廓顯著性耦合的情況下實現(xiàn)軌道根數(shù)反演。針對該問題，本文研究了顯著性數(shù)據(jù)的分布特征，并根據(jù)其分布特點，提出利用該數(shù)據(jù)的方差這一數(shù)字特征作為軌道根數(shù)最優(yōu)估計的依據(jù)，實現(xiàn)觀測輪廓顯著性耦合情況下的軌道六根數(shù)的估計，并利用實測數(shù)據(jù)驗證了方法的有效性。

1 顯著性耦合下的軌道根數(shù)反演方法

圖1 顯著性耦合下的軌道根數(shù)反演方法Fig.1 Flow chart of the orbital elements inversion method for significant coupling

1.1 軌道動力學模型

將航天器視為質(zhì)點，在J2000.0地心慣性坐標系(ECI)下，如圖2所示坐標系OXYZ。航天器的軌道動力學模型可表示為

(1)

式中：r與v分別為航天器的位置和速度；a為航天器的加速度。航天器在軌道上運動要受到各種力的影響，產(chǎn)生的攝動加速度是多方面的。為避免定軌過程中計算量過大，建立如下的近地軌道衛(wèi)星的預測模型，該模型由如下加速度項組成。

a=aTB+aNS+aS+aM+aD+aSP

(2)

式中：aTB為航天器所受到的地球二體引力加速度；aNS為地球非球形攝動加速度；aS與aM分別為太陽引力與月球引力的攝動加速度；aD為大氣阻力攝動加速度；aSP為太陽光壓攝動加速度。表1給出了文中使用的動力學模型中攝動力計算使用的模型名稱，計算公式可參考文獻[20]。

表1 軌道動力學模型Table 1 Descriptions of orbital dynamics models

軌道六根數(shù)分別為半長軸a、偏心率e、升交點赤經(jīng)Ω(RAAN)、軌道傾角i、近地點幅角ω和真近點角θ。

描述航天器運動最方便的坐標系之一是近焦點坐標系，如圖2中坐標系Oxyz，該坐標系的基準面是航天器的軌道平面，原點為地球質(zhì)心，x軸指向近地點，z軸沿航天器角動量方向，y軸與x軸、z軸構(gòu)成右手直角坐標系[21]。在近焦點坐標系中，航天器的位置和速度矢量可通過軌道六根數(shù)來計算，表示為式(3)與式(4)[20]。

圖2 坐標系示意圖Fig.2 Schematic diagram of the coordinate system

(3)

(4)

(5)

式中：坐標變換矩陣Q為[21]

(6)

1.2 觀測脈沖輪廓獲取

1.2.1光子TOA校正

X射線脈沖星信號微弱，只能以光子的形式被探測器收到，由于引力、相對論等效應的影響，在獲得觀測脈沖輪廓之前，需要對接收到的光子TOA進行校正。對脈沖單星來說，只需轉(zhuǎn)換到太陽系質(zhì)心(SSB)處即可。文中使用時間轉(zhuǎn)換公式如下[22]

tSSB-tsc=ΔC+ΔR+ΔE+ΔS

(7)

式中：ΔC為時鐘校正；ΔR為Roemer延遲，表示光子在真空中從航天器到SSB的傳播時間；ΔE為Einstein延遲，是由太陽系內(nèi)天體的運動所引起的引力紅移及時間變慢效應；ΔS為Shapiro延遲，由天體引力引起的時空彎曲矯正。由于時間轉(zhuǎn)換項完整的表達式非常復雜，文中采用如式(8)所示的簡化公式，式中的各變量均在太陽系質(zhì)心慣性坐標系下表示，該模型忽略了視差、脈沖星自行以及除太陽外的其他天體的Shapiro延遲等項，轉(zhuǎn)換精度為5～8 μs[23]。

(8)

式中：第一項為Roemer延遲，第二項為Einstein延遲，第三項為太陽Shapiro延遲。n為脈沖星在太陽系質(zhì)心坐標系下的單位方向矢量；rsc/ssb為航天器相對于太陽系質(zhì)心的位置；rsc/E為航天器相對于地球的位置；vE為地球的速度；rsc/sun為航天器相對于太陽的位置；c為光速；μs為太陽引力常數(shù)；P是周期項,可表示為[24]

P≈0.0016568sin(35999.37T+357.5)+

0.0000224sin(32964.5T+246)+

0.0000224sin(32964.5T+246)+

0.0000138sin(71998.7T+355)+

0.0000048sin(3034.9T+25)+

0.0000048sin(34777.3T+230)

(9)

式中：T=(JD-2452545.0)/36525, JD為儒略日。

1.2.2觀測脈沖輪廓獲取

利用式(8)將探測器探測到的光子TOA進行時間校正后，為獲得觀測脈沖輪廓，還需要將光子TOA轉(zhuǎn)換為光子到達相位，使用時間相位模型可完成轉(zhuǎn)換，時間相位模型如式(10)所示[15]。

(10)

式中：下標0表示在參考歷元t0時刻的自轉(zhuǎn)和相位參數(shù)，φ0即t0時刻的初始相位，自轉(zhuǎn)頻率k-1階偏導數(shù)ν(k-1)是待擬合的參數(shù)，k=1時，ν(0)=ν0，即t0時刻的初始自轉(zhuǎn)頻率。式(10)相當于對脈沖周期的標準化處理，一個脈沖周期對應的相位區(qū)間長度為1。

假設觀測周期為Np個脈沖周期，即光子到達相位序列可劃分為Np個長度為1的相位區(qū)間。將每一個相位區(qū)間劃分為Nb個等長度的片段，且將獲得的光子到達相位序列折疊到相位區(qū)間(0,1)內(nèi)的Nb個等長度片段，用第i個片段的中心φi表示該片段內(nèi)的光子到達相位。通過對相位區(qū)間(0,1)的每一個片段內(nèi)的光子數(shù)進行統(tǒng)計，可得到平均觀測脈沖輪廓，如式(11)所示。

(11)

式中：Ci表示折疊后相位區(qū)間(0,1)內(nèi)的第i個片段內(nèi)的光子總數(shù)；cj(φi)表示第j個相位區(qū)間內(nèi)第i個片段內(nèi)的光子個數(shù)。將式(11)表示的觀測輪廓進行歸一化處理，如式(12)所示。

(12)

1.3 顯著性耦合下的軌道根數(shù)反演方法

(13)

(14)

n=1,2,…,Na

(15)

(16)

(17)

將上述過程重復進行，依次可獲得Ω,i,ω的估計值，如式(18)，(19)與(20)所示。

(18)

(19)

(20)

(21)

2 仿真校驗

2.1 仿真條件

選用RXTE衛(wèi)星的PSR B0531+21脈沖星實測數(shù)據(jù)包95802- 01- 16- 02，95802- 01- 16- 03與95802- 01- 16- 04，對RXTE衛(wèi)星進行導航實驗，脈沖星參數(shù)如表2所示。

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

表3 初始軌道參數(shù)Table 3 Initial orbit parameters

表4 軌道根數(shù)的偏移量Table 4 The bias of the orbital elements

本文仿真實驗使用的計算平臺為MATLAB，計算機硬件平臺為DELL T96LKJ4工作站，仿真實驗使用16個線程并行運算，計算耗時約4 h。

2.2 仿真結(jié)果與分析

2.2.1耦合情況下的顯著性卡方

圖3 軌道六根數(shù)均存在誤差時的顯著性卡方值Fig.3 Significance chi square values for deviations in six orbital elements

2.2.2耦合情況下的軌道根數(shù)反演

圖4 分組卡方值的統(tǒng)計結(jié)果(根據(jù)軌道半長軸分組)Fig.4 Statistical results of grouped chi square values (Grouped according toorbit semi-major axis)

圖5 分組卡方的方差(根據(jù)軌道半長軸分組)Fig.5 Variance of grouped chi square values (Grouped according to orbit semi-major axis)

圖6 軌道偏心率的搜索過程Fig.6 The search process of orbit eccentricity

圖7 確定軌道偏心率后的顯著性卡方Fig.7 Significance chi square after determining orbit eccentricity

圖8 升交點赤經(jīng)的搜索過程Fig.8 The search process of RAAN

圖9 根據(jù)軌道傾角分組的卡方方差Fig.9 Variance of chi square grouped according to orbital inclination

軌道傾角確定后，只有近地點幅角ω與真近點角θ兩個軌道根數(shù)耦合。圖10為根據(jù)近地點幅角分組的卡方方差，通過對方差曲線擬合，并取方差最小值對應的近地點幅角，可得到近地點幅角的估計值。

圖10 根據(jù)近地點幅角分組的卡方方差Fig.10 Variance of chi square grouped according to argument of perigee

圖11 確定近地點幅角后的顯著性卡方Fig.11 Significance chi square after determining argument of perigee

通過上述過程，可逐一搜索得到軌道六根數(shù)，為避免脈沖星觀測數(shù)據(jù)的隨機性對本文方法實驗結(jié)果的影響，在相同條件下，多次仿真光子到達時間序列，并利用每一組光子數(shù)據(jù)反演軌道根數(shù)，計算軌道根數(shù)的搜索誤差。重復進行100次仿真實驗，每一次實驗可獲得一組軌道根數(shù)，并計算得到軌道根數(shù)的估計誤差。圖12給出了每個軌道根數(shù)搜索誤差的統(tǒng)計圖，對每個軌道根數(shù)的所有實驗結(jié)果取平均值，以此作為軌道根數(shù)的估計誤差，如表5所示，給出了軌道根數(shù)的平均估計誤差及對應的位置、速度誤差。

圖12 軌道根數(shù)搜索誤差的統(tǒng)計結(jié)果Fig.12 Statistical results of searching error of orbit elements

表5 軌道根數(shù)估計誤差Table 5 Estimation errors of orbital elements

在相同的實驗條件下，比較本文方法的搜索結(jié)果與利用文獻[13]方法得到的軌道根數(shù)搜索結(jié)果，如表6所示。根據(jù)表中結(jié)果，只有真近點角的搜索誤差與本文的結(jié)果是一樣的，其他軌道根數(shù)的搜索誤差很大，這是因為在搜索真近點角時，其他軌道根數(shù)已經(jīng)搜索出結(jié)果，此時變?yōu)閱蝹€軌道根數(shù)的搜索問題，證明文獻[13]的方法主要適用于單個軌道根數(shù)的搜索問題，不適用于軌道六根數(shù)同時存在誤差的情況。

表6 軌道根數(shù)估計誤差對比Table 6 Comparison of estimation errors of orbital elements

2.2.3橢圓軌道的軌道根數(shù)反演

上一節(jié)在近圓軌道上對文中所提方法進行了驗證，本節(jié)在橢圓軌道上對所提軌道根數(shù)反演方法進行實驗驗證。橢圓軌道根數(shù)如表7所示。

表7 橢圓軌道參數(shù)Table 7 Elliptical orbit parameters

圖13 軌道根數(shù)逐一分組后的顯著性卡方的方差Fig.13 Variance of chi square of significance after grouping the orbital elements one by one

表8 軌道根數(shù)的偏移量Table 8 The bias of the orbital elements

3 結(jié) 論