基于實測星載原子鐘數據的脈沖星守時方法

王禹淞,王奕迪,鄭 偉

(國防科技大學空天科學學院,長沙 410073)

0 引言

衛星導航系統是以時間測量為基礎的系統,精確的位置測量本質上依賴于精確的時間測量[1]。導航衛星的時間參考由其攜帶的星載原子鐘提供,因此,星載原子鐘的鐘差直接影響了衛星導航系統的定位與導航精度[2-4]。當前,星載原子鐘的鐘差修正主要依靠地面測控系統測量并將修正量上傳至衛星[5]。該方法精度較高且十分穩定,但也存在地面測控系統的負擔大、衛星導航系統對地面測控系統依賴程度高的問題[5]。

脈沖星是一類高速自轉的中子星,其自轉的長期穩定性極佳,被譽為自然界最穩定的天然時鐘[6]。由于脈沖星具有良好的長期穩定性,眾多學者在脈沖星長期穩定度分析以及脈沖星計時等方面開展了許多研究工作[7-8]。Hobbs等[9]利用國際脈沖星計時陣(international pulsar timing array,IPTA)第一次釋放的數據,綜合利用48顆毫秒脈沖星長期計時資料,構建了綜合脈沖星時。Li等[10]使用昆明40 m 射電望遠鏡的脈沖星計時觀測數據,實現了80 ns量級的本地時鐘鐘跳檢測。Deneva等[11]利用NICER 探測器1 年左右的PSR B1937+21 和9 個月左右的PSR B1821-24觀測數據,開展了X 射線脈沖星計時分析。Sun等[12]利用NICER 探測器5年的PSR B1937+21脈沖星的X 射線觀測數據進行計時分析,并與利用射電波段觀測數據的計時結果進行了對比。除了脈沖星計時外,還可利用脈沖星觀測數據進行星載原子鐘鐘差的估計,從而實現導航衛星自主守時[13]。Sun等[13]分析了脈沖星守時可行性,提出了基于卡爾曼濾波的原子鐘鐘差估計方法。王奕迪等[14]在考慮脈沖星到達時間測量存在系統誤差的情況下,提出了二步卡爾曼濾波算法,有效降低了系統誤差對守時性能的影響。童明雷等[15]利用脈沖星實測數據,提出了利用脈沖星校準星載原子鐘頻率的方法,為脈沖星守時提供了新的思路。

在上述研究中,星載原子鐘的鐘差、頻率等數據均為仿真生成,并不能完全反映脈沖星守時系統的性能。本文利用武漢大學數據中心發布的星載原子鐘鐘差數據,研究和分析了利用脈沖星對星載原子鐘鐘差的估計方法,為驗證脈沖星守時系統的可行性提供了參考。

1 星載原子鐘數據處理與分析

1.1 數據預處理

由于星載原子鐘在運行過程中受到環境變化、設備老化等因素的影響,其鐘差和頻率數據中可能存在粗差或中斷[1-2]。因此,在進行星載原子鐘性能評估前,需要對原始鐘差數據進行預處理。對于數據中斷的情況,可采用線性插值進行擬合[2]。對于數據粗差,需要對其進行檢測和剔除[1]。通常,粗差的檢測一般采用中位數法(median absolute deviation,MAD),檢測對象通常為由原始鐘差數據轉化得到的頻率數據。鐘差數據轉化為頻率數據的方法為

式中,xi表示第i個鐘差數據;yi表示i個頻率數據;Δt為兩個鐘差數據之間的時間間隔。

將鐘差數據轉化為頻率數據后,將滿足以下條件的數據認為是粗差并加以剔除

式中,m為y序列的中位數;n為經驗取值,一般為3～5;M=Median)/0.674 5。 此外,對于存在較大頻率漂移的數據,需要先去除頻率漂移的趨勢項,再進行粗差檢測。圖1為G07衛星(銣鐘)剔除粗差后的頻率數據和鐘差數據。

圖1 G07的鐘差和頻差數據Fig.1 Clock deviation and relative frequency deviation data of G07

1.2 原子鐘鐘差模型

多項式模型是目前最常用的鐘差模型之一,對于一組鐘差序列xi(i=1,2,…,N),xi可用多項式模型表示為[16]

式中,a0,a1,a2,…,an為n個模型參數;n為多項式階數;ei為模型誤差。通常,銣鐘選用2階模型,銫鐘和氫鐘選用1階模型,本文的研究主要考慮2階模型的情況。

根據式(3),銣鐘的狀態方程可表示為[16]

式中,τ為采樣時間;x(t),y(t),z(t)分別為原子鐘的鐘差、鐘速和鐘飄;εx,εy,εz為均值為0的模型隨機誤差。

式(4)可簡化為

式中,Xk=[x(t+τ)y(t+τ)z(t+τ)]為tk時刻的系統狀態;Wk為系統模型誤差,其協方差陣Qk可表示為[16]

式中,q1為對應εx的過程噪聲參數,表現為調相隨機游走噪聲;q2為對應于εy的過程噪聲參數,表現為調頻隨機游走噪聲;q3為對應于εz的過程噪聲參數,表現為調頻隨機奔跑噪聲。

1.3 頻率穩定度分析

頻率穩定度是衡量原子鐘性能的重要指標之一,反映了原子鐘的輸出頻率受噪聲的影響而產生的隨機起伏的情況[1]。阿倫方差(Allan variance)是一種常用的在時域評估穩定度的方法。對于采樣間隔為τ0的頻率數據序列{yn,n=1,2,…,L},該序列的阿倫方差表達式為[16]

式中,τ=lτ0表示平滑時間;ˉyi(l)=表示第i個平滑時間內的l個頻率數據的均值;L'為的個數。本文利用武漢大學分析中心發布的鐘差數據(采樣間隔為900 s)根據G07衛星自2016年7月1日至2020年6月26日的數據,計算了其阿倫標準差,結果如圖2所示,該圖為平滑時間和阿倫標準差的雙對數坐標圖。從圖中可以看出,隨著平滑時間的增長,阿倫標準差先逐漸減小,后逐漸增大。

圖2 基于Allan標準差的頻率穩定度評估結果Fig.2 Results of frequency stability based on Allan deviation

阿倫方差還可用于原子鐘過程噪聲參數的估計,具體方法為[16]

根據式(8),采用最小二乘法,即可計算出過程噪聲參數q1,q2,q3。

1.4 原子鐘誤差特性分析

若真實的原子鐘鐘差系數已知,則原子鐘的鐘差僅受隨機性誤差的影響。因此,首先分析鐘差的隨機性誤差。本文以鐘差的功率譜表征鐘差的隨機性誤差水平。根據原子鐘的鐘差數據,計算了其功率譜,如圖3所示,該圖為頻率和功率譜的雙對數坐標圖。圖中,藍色線為根據鐘差數據計算的功率譜,紅色直線為對功率譜進行擬合的曲線,黃色和紫色曲線分別表示其功率譜的上、下包絡。

圖3 G07的功率譜Fig.3 Power spectrum of G07

為對比脈沖星與原子鐘的性能差異,利用NICER 探測器的PSR B1937+21 脈沖星在2017年至2019年的實測數據,計算脈沖星的功率譜,并將其與原子鐘的功率譜進行對比,結果如圖4所示,該圖為頻率和功率譜的雙對數坐標圖。圖中,藍色線表示脈沖星的功率譜,紅色線為脈沖星功率譜的擬合曲線,黃色線為原子鐘功率譜的擬合曲線,紫色和綠色的虛線為原子鐘功率譜的上、下包絡。

圖4 G07與脈沖星的功率譜對比Fig.4 Compare of the power spectrums of G07 and PSR B1937+21

如圖4所示,脈沖星的擬合曲線(紅色)與銣鐘的擬合曲線(黃色)在約5.3天相交,之后脈沖星的功率譜擬合曲線(紅色)低于銣鐘的功率譜擬合曲線(黃色)。在約33.2天后,脈沖星實測功率譜(藍色線)低于銣鐘功率譜的下包絡。因此,脈沖星的長期性能優于原子鐘,可利用脈沖星的觀測數據校準原子鐘。

圖1 給出了G07(銣鐘)的鐘差和頻差數據?？梢钥闯?該銣鐘在約3年的時間內,鐘差變化了接近800μs。因此,在含系統性誤差影響的情況下,星載原子鐘的鐘差變化很大。

圖5給出了PSR B1937+21 的計時殘差變化情況,可以看出,在接近3年的時間內,脈沖星的計時殘差不存在系統性誤差,始終不超過15μs。因此,雖然PSR B1937+21的隨機性誤差在33天才優于G07號衛星的原子鐘,但是由于星載原子鐘存在系統性誤差,即使在較短的時間內,也可以通過脈沖星來估計原子鐘的鐘差。

圖5 PSR B1937+21的計時殘差Fig.5 Residual of PSR B1937+21

2 脈沖星守時方法

2.1 脈沖星守時方案設計

脈沖星守時方案如圖6所示。首先,通過星載原子鐘測量脈沖星的脈沖到達時間。而后,利用該脈沖到達時間,對星載原子鐘的鐘差進行估計。從圖中可以看出,脈沖星守時的關鍵,在于如何利用脈沖到達時間,實現原子鐘鐘差的準確估計。

圖6 脈沖星守時總體方案Fig.6 The overall plan for pulsar timing

圖7給出了利用脈沖星估計原子鐘鐘差的方法示意圖。首先,利用星載原子鐘,可以測量脈沖到達時間。此外,對原子鐘的鐘差模型進行建模。結合脈沖到達時間和原子鐘的鐘差模型,利用卡爾曼濾波算法,即可實現原子鐘的鐘差估計。

圖7 利用脈沖星估計原子鐘鐘差方法示意圖Fig.7 Schematic diagram of using pulsar to estimate atomic clock deviation method

2.2 脈沖星測量模型

脈沖星的時間相位模型是實現脈沖星守時的基礎,時間相位模型通常建立在太陽系質心(solar system barycenter,SSB)。假設由星載原子鐘記錄的脈沖到達衛星的時間為tSC,則該脈沖到達SSB的時間為[17]

式中,n表示脈沖星方向矢量;r(tSC)表示衛星相對于地球的位置矢量;rE(tSC)表示地球相對于SSB的位置矢量;c表示光速;Mk表示第k個天體的質量;rk(tSC)表示第k個天體相對于衛星的位置矢量;H.O.T 表示可忽略的高階項。

在SSB處的脈沖星時間相位模型為[17]

式中,φ表示脈沖相位;φ0表示tSSB,0=g(t0)時刻的脈沖相位;f(n)為tSSB,0時刻的脈沖星自轉頻率的n階導數。利用式(10)所示的時間相位模型,即可預測脈沖星信號到達SSB的時間SSB。

令

因此,脈沖到達時間的預測值與測量值之差即反映了原子鐘的鐘差。因此,脈沖星測量方程可表示為

2.3 脈沖到達時間精度分析

脈沖到達時間的精度可通過以下公式表示[18]

式中,T50是半流量密度持續時間,Tb是敏感器的時間分辨率,B是背景噪聲的光子數,S是來自脈沖星的光子數。

式中,Δt是觀測持續時間,λp和λn分別是脈沖信號和背景噪聲的平均流量密度,A是探測器有效面積。

將式(15)(16)代入式(14),最后可得脈沖TOA測量精度的表達式為

3 守時性能分析

本章以武漢大學發布的G07衛星自2016年7月1日至2020年6月26日的鐘差數據為基礎,分析脈沖星守時系統性能。假設脈沖到達時間的測量精度為1μs,觀測周期為30天。原子鐘鐘差的過程噪聲參數通過式(8)進行擬合。在進行脈沖星守時計算時,首先使用前3個測量值,利用最小二乘法計算原子鐘的狀態及其狀態協方差矩陣。而后,將最小二乘計算的原子鐘狀態和狀態協方差矩陣作為卡爾曼濾波的初值進行解算。

圖8(a)給出了脈沖星守時系統的仿真計算結果,可以看出,最初用最小二乘法解算的結果鐘差較大,之后脈沖星守時系統精度逐漸收斂至1μs左右。圖8(b)給出了利用批量最小二乘的脈沖星守時系統的鐘差,可以看出,該守時系統的鐘差保持在6μs以內,守時精度低于卡爾曼濾波方法。這是因為,在批量最小二乘算法中未考慮鐘差模型的過程噪聲,導致基于批量最小二乘的脈沖星守時系統鐘差大于基于卡爾曼濾波的脈沖星守時系統鐘差。

圖8 脈沖星守時系統仿真結果Fig.8 Simulation results of pulsar timing system

對基于卡爾曼濾波和批量最小二乘的脈沖星守時系統進行了500次蒙特卡羅仿真計算,得到基于卡爾曼濾波的守時系統精度為0.78μs,基于批量最小二乘的守時系統精度為1.86μs?？梢钥闯?基于卡爾曼濾波的脈沖星守時系統精度優于基于批量最小二乘的脈沖星守時系統。因此,分析并在鐘差估計過程中考慮星載原子鐘模型的過程噪聲,對提高守時系統的性能是有意義的。

分析脈沖星觀測策略和探測器有效面積對脈沖星守時系統性能的影響。表1給出了不同觀測策略下的脈沖星守時系統精度和對應需要的探測器有效面積。根據表中前3列可以看出,相同的觀測周期下,脈沖到達時間解算精度越高,則脈沖星守時系統精度越高。脈沖到達時間解算精度相同的情況下,觀測周期越短則脈沖星守時系統精度越高。在計算探測器有效面積時,假設觀測的脈沖星為PSR B1937+21,其半流量密度持續時間為21μs,λp=0.145 photon/(cm2·s),λn=12 photon/(cm2·s),敏感器的時間分辨率為0.1μs。根據表中第3列和第4列可以看出,對于觀測同一顆脈沖星的情況,探測器有效面積越大,則脈沖星守時系統的性能越好。增大探測器的有效面積對提高守時系統的性能具有積極作用。然而,探測器的有效面積增大意味著航天器的搭載成本提高。因此,需要根據實際的工程應用場景和任務需求,對探測器有效面積進行設計。

表1 不同觀測策略下的脈沖星守時系統性能及探測器有效面積Tab.1 Performance of pulsar timing systems and effective detector areas under different observation strategies

此外,需要說明的是,在本文的仿真分析中,脈沖星觀測策略及需要的探測器有效面積均是根據理論公式計算的,與實際空間環境下真實的脈沖星觀測性能可能存在一定差距。在后續的脈沖星守時方法研究中,應該結合X 射線探測器的實際觀測能力或國內外脈沖星實測數據開展分析與研究。

4 結論

本文利用實測星載原子鐘鐘差數據和脈沖星觀測數據,開展了脈沖星守時方法研究,得到如下結論:

1)通過對比脈沖星計時殘差和星載原子鐘鐘差的功率譜,可以看出脈沖星的長期性能優于星載原子鐘。

2)由于包含系統誤差的影響,星載原子鐘的鐘差變化劇烈,G07號衛星的原子鐘在約3年的時間內,鐘差變化了接近800μs。而脈沖星的計時殘差幾乎不存在系統誤差,PSR B1937+21在接近3年的時間內計時殘差始終不超過15μs。因此,雖然PSR B1937+21 的隨機性誤差在33 天后才優于G07號衛星的原子鐘,但是由于星載原子鐘存在系統性誤差,即使在較短的時間內,也可以通過脈沖星來估計原子鐘的鐘差。

3)基于實測星載原子鐘鐘差數據,開展了脈沖星守時性能分析。計算結果表明,若脈沖星的脈沖到達時間解算精度為1μs/30 d,則原子鐘鐘差估計精度可達到優于1μs的水平,初步驗證了脈沖星守時系統的可行性。