觀測輪廓顯著性耦合時的軌道根數(shù)反演方法

蘇劍宇，方海燕，高敬敬

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

0 引 言

作為深空網(wǎng)絡(luò)(Deep space network，DSN)的補(bǔ)充和替代技術(shù)，X射線脈沖星導(dǎo)航(X-ray pulsar-based navigation，XPNAV)是一種新型的自主導(dǎo)航技術(shù)。XPNAV可為在太陽系內(nèi)外飛行的航天器提供位置、時間等導(dǎo)航信息[1-4]，實現(xiàn)航天器高精度自主導(dǎo)航。目前，隨著X射線探測器技術(shù)的不斷進(jìn)步，國內(nèi)外相繼開展了X射線脈沖星導(dǎo)航試驗。如中國空間實驗室天宮二號的γ暴偏振探測科學(xué)實驗[5]，中國首顆X射線脈沖星導(dǎo)航試驗衛(wèi)星(XPNAV-1)在軌開展的X射線脈沖星的探測與脈沖星導(dǎo)航體制的驗證[6-8]，以及國內(nèi)首顆空間X射線天文硬X射線調(diào)制望遠(yuǎn)鏡衛(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導(dǎo)航驗證的主要框架可分為兩種類型：實時和非實時[13]。前者利用脈沖星觀測周期內(nèi)記錄到的一串光子到達(dá)時間(Time of arrivals，TOA)，通過一定的算法提取出脈沖相位以及多普勒頻率，由于航天器在任意時刻觀測的脈沖相位和頻率都可以用該時刻的位置、速度以及太陽系質(zhì)心(Solar system barycenter, SSB)處的脈沖星信號模型準(zhǔn)確地表示，因此通過對相位與多普勒頻率的解算即可得到相應(yīng)時刻的位置與速度信息[14-16]。如SEXTANT任務(wù)通過對多個毫秒脈沖星的觀測數(shù)據(jù)進(jìn)行處理，得到10公里以內(nèi)的導(dǎo)航精度[17]。文獻(xiàn)[18]在XPNAV-1的軌道傳播中引入脈沖星視線相對于地球的測距測量作為控制點(diǎn)，通過對脈沖星PSR B0531+21近三個月的觀測，在控制點(diǎn)進(jìn)行自我定位，平均導(dǎo)航誤差為38.4公里。

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

綜合以上分析，上述文獻(xiàn)方法雖然可在單個軌道根數(shù)存在偏差時，根據(jù)觀測輪廓顯著性反演軌道根數(shù)，但無法在軌道六根數(shù)均存在偏差時，即觀測輪廓顯著性耦合的情況下實現(xiàn)軌道根數(shù)反演。針對該問題，本文研究了顯著性數(shù)據(jù)的分布特征，并根據(jù)其分布特點(diǎn)，提出利用該數(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 軌道動力學(xué)模型

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

(1)

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

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

(2)

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

表1 軌道動力學(xué)模型Table 1 Descriptions of orbital dynamics models

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

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

圖2 坐標(biāo)系示意圖Fig.2 Schematic diagram of the coordinate system

(3)

(4)

(5)

式中：坐標(biāo)變換矩陣Q為[21]

(6)

1.2 觀測脈沖輪廓獲取

1.2.1光子TOA校正

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

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

(7)

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

(8)

式中：第一項為Roemer延遲，第二項為Einstein延遲，第三項為太陽Shapiro延遲。n為脈沖星在太陽系質(zhì)心坐標(biāo)系下的單位方向矢量；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進(jìn)行時間校正后，為獲得觀測脈沖輪廓，還需要將光子TOA轉(zhuǎn)換為光子到達(dá)相位，使用時間相位模型可完成轉(zhuǎn)換，時間相位模型如式(10)所示[15]。

(10)

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

假設(shè)觀測周期為Np個脈沖周期，即光子到達(dá)相位序列可劃分為Np個長度為1的相位區(qū)間。將每一個相位區(qū)間劃分為Nb個等長度的片段，且將獲得的光子到達(dá)相位序列折疊到相位區(qū)間(0,1)內(nèi)的Nb個等長度片段，用第i個片段的中心φi表示該片段內(nèi)的光子到達(dá)相位。通過對相位區(qū)間(0,1)的每一個片段內(nèi)的光子數(shù)進(jìn)行統(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)表示的觀測輪廓進(jìn)行歸一化處理，如式(12)所示。

(12)

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

(13)

(14)

n=1,2,…,Na

(15)

(16)

(17)

將上述過程重復(fù)進(jìn)行，依次可獲得Ω,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)星進(jìn)行導(dǎo)航實驗，脈沖星參數(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，計算機(jī)硬件平臺為DELL T96LKJ4工作站，仿真實驗使用16個線程并行運(yùn)算，計算耗時約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 升交點(diǎn)赤經(jīng)的搜索過程Fig.8 The search process of RAAN

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

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

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

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

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

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

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

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

上一節(jié)在近圓軌道上對文中所提方法進(jìn)行了驗證，本節(jié)在橢圓軌道上對所提軌道根數(shù)反演方法進(jìn)行實驗驗證。橢圓軌道根數(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é) 論