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

2022-02-01 13:42:02蘇劍宇方海燕高敬敬

宇航學報 2022年12期

蘇劍宇，方海燕，高敬敬

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

0 引言

作為深空網絡(Deep space network，DSN)的補充和替代技術，X射線脈沖星導航(X-ray pulsar-based navigation，XPNAV)是一種新型的自主導航技術。XPNAV可為在太陽系內外飛行的航天器提供位置、時間等導航信息[1-4]，實現航天器高精度自主導航。目前，隨著X射線探測器技術的不斷進步，國內外相繼開展了X射線脈沖星導航試驗。如中國空間實驗室天宮二號的γ暴偏振探測科學實驗[5]，中國首顆X射線脈沖星導航試驗衛星(XPNAV-1)在軌開展的X射線脈沖星的探測與脈沖星導航體制的驗證[6-8]，以及國內首顆空間X射線天文硬X射線調制望遠鏡衛星(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]。

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

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

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

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

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

1.1 軌道動力學模型

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

(1)

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

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

(2)

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

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

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

描述航天器運動最方便的坐標系之一是近焦點坐標系，如圖2中坐標系Oxyz，該坐標系的基準面是航天器的軌道平面，原點為地球質心，x軸指向近地點，z軸沿航天器角動量方向，y軸與x軸、z軸構成右手直角坐標系[21]。在近焦點坐標系中，航天器的位置和速度矢量可通過軌道六根數來計算，表示為式(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進行校正。對脈沖單星來說，只需轉換到太陽系質心(SSB)處即可。文中使用時間轉換公式如下[22]

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

(7)

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

(8)

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

P≈0.0016568sin(35999.37T+357.5)+

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轉換為光子到達相位，使用時間相位模型可完成轉換，時間相位模型如式(10)所示[15]。

(10)

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

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

(11)

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

(12)

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

(13)

(14)

n=1,2,…,Na

(15)

(16)

(17)

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

(18)

(19)

(20)

(21)

2 仿真校驗

2.1 仿真條件

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

表2 脈沖星參數Table 2 Pulsar parameters

表3 初始軌道參數Table 3 Initial orbit parameters

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

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

2.2 仿真結果與分析

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

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

2.2.2耦合情況下的軌道根數反演

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

圖5 分組卡方的方差(根據軌道半長軸分組)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 升交點赤經的搜索過程Fig.8 The search process of RAAN

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

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

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

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

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

圖12 軌道根數搜索誤差的統計結果Fig.12 Statistical results of searching error of orbit elements

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

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