基于實(shí)測(cè)星載原子鐘數(shù)據(jù)的脈沖星守時(shí)方法

王禹淞,王奕迪,鄭 偉

(國(guó)防科技大學(xué)空天科學(xué)學(xué)院,長(zhǎng)沙 410073)

0 引言

衛(wèi)星導(dǎo)航系統(tǒng)是以時(shí)間測(cè)量為基礎(chǔ)的系統(tǒng),精確的位置測(cè)量本質(zhì)上依賴(lài)于精確的時(shí)間測(cè)量[1]。導(dǎo)航衛(wèi)星的時(shí)間參考由其攜帶的星載原子鐘提供,因此,星載原子鐘的鐘差直接影響了衛(wèi)星導(dǎo)航系統(tǒng)的定位與導(dǎo)航精度[2-4]。當(dāng)前,星載原子鐘的鐘差修正主要依靠地面測(cè)控系統(tǒng)測(cè)量并將修正量上傳至衛(wèi)星[5]。該方法精度較高且十分穩(wěn)定,但也存在地面測(cè)控系統(tǒng)的負(fù)擔(dān)大、衛(wèi)星導(dǎo)航系統(tǒng)對(duì)地面測(cè)控系統(tǒng)依賴(lài)程度高的問(wèn)題[5]。

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

在上述研究中,星載原子鐘的鐘差、頻率等數(shù)據(jù)均為仿真生成,并不能完全反映脈沖星守時(shí)系統(tǒng)的性能。本文利用武漢大學(xué)數(shù)據(jù)中心發(fā)布的星載原子鐘鐘差數(shù)據(jù),研究和分析了利用脈沖星對(duì)星載原子鐘鐘差的估計(jì)方法,為驗(yàn)證脈沖星守時(shí)系統(tǒng)的可行性提供了參考。

1 星載原子鐘數(shù)據(jù)處理與分析

1.1 數(shù)據(jù)預(yù)處理

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

式中,xi表示第i個(gè)鐘差數(shù)據(jù);yi表示i個(gè)頻率數(shù)據(jù);Δt為兩個(gè)鐘差數(shù)據(jù)之間的時(shí)間間隔。

將鐘差數(shù)據(jù)轉(zhuǎn)化為頻率數(shù)據(jù)后,將滿足以下條件的數(shù)據(jù)認(rèn)為是粗差并加以剔除

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

圖1 G07的鐘差和頻差數(shù)據(jù)Fig.1 Clock deviation and relative frequency deviation data of G07

1.2 原子鐘鐘差模型

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

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

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

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

式(4)可簡(jiǎn)化為

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

式中,q1為對(duì)應(yīng)εx的過(guò)程噪聲參數(shù),表現(xiàn)為調(diào)相隨機(jī)游走噪聲;q2為對(duì)應(yīng)于εy的過(guò)程噪聲參數(shù),表現(xiàn)為調(diào)頻隨機(jī)游走噪聲;q3為對(duì)應(yīng)于εz的過(guò)程噪聲參數(shù),表現(xiàn)為調(diào)頻隨機(jī)奔跑噪聲。

1.3 頻率穩(wěn)定度分析

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

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

圖2 基于Allan標(biāo)準(zhǔn)差的頻率穩(wěn)定度評(píng)估結(jié)果Fig.2 Results of frequency stability based on Allan deviation

阿倫方差還可用于原子鐘過(guò)程噪聲參數(shù)的估計(jì),具體方法為[16]

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

1.4 原子鐘誤差特性分析

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

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

為對(duì)比脈沖星與原子鐘的性能差異,利用NICER 探測(cè)器的PSR B1937+21 脈沖星在2017年至2019年的實(shí)測(cè)數(shù)據(jù),計(jì)算脈沖星的功率譜,并將其與原子鐘的功率譜進(jìn)行對(duì)比,結(jié)果如圖4所示,該圖為頻率和功率譜的雙對(duì)數(shù)坐標(biāo)圖。圖中,藍(lán)色線表示脈沖星的功率譜,紅色線為脈沖星功率譜的擬合曲線,黃色線為原子鐘功率譜的擬合曲線,紫色和綠色的虛線為原子鐘功率譜的上、下包絡(luò)。

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

如圖4所示,脈沖星的擬合曲線(紅色)與銣鐘的擬合曲線(黃色)在約5.3天相交,之后脈沖星的功率譜擬合曲線(紅色)低于銣鐘的功率譜擬合曲線(黃色)。在約33.2天后,脈沖星實(shí)測(cè)功率譜(藍(lán)色線)低于銣鐘功率譜的下包絡(luò)。因此,脈沖星的長(zhǎng)期性能優(yōu)于原子鐘,可利用脈沖星的觀測(cè)數(shù)據(jù)校準(zhǔn)原子鐘。

圖1 給出了G07(銣鐘)的鐘差和頻差數(shù)據(jù)。可以看出,該銣鐘在約3年的時(shí)間內(nèi),鐘差變化了接近800μs。因此,在含系統(tǒng)性誤差影響的情況下,星載原子鐘的鐘差變化很大。

圖5給出了PSR B1937+21 的計(jì)時(shí)殘差變化情況,可以看出,在接近3年的時(shí)間內(nèi),脈沖星的計(jì)時(shí)殘差不存在系統(tǒng)性誤差,始終不超過(guò)15μs。因此,雖然PSR B1937+21的隨機(jī)性誤差在33天才優(yōu)于G07號(hào)衛(wèi)星的原子鐘,但是由于星載原子鐘存在系統(tǒng)性誤差,即使在較短的時(shí)間內(nèi),也可以通過(guò)脈沖星來(lái)估計(jì)原子鐘的鐘差。

圖5 PSR B1937+21的計(jì)時(shí)殘差Fig.5 Residual of PSR B1937+21

2 脈沖星守時(shí)方法

2.1 脈沖星守時(shí)方案設(shè)計(jì)

脈沖星守時(shí)方案如圖6所示。首先,通過(guò)星載原子鐘測(cè)量脈沖星的脈沖到達(dá)時(shí)間。而后,利用該脈沖到達(dá)時(shí)間,對(duì)星載原子鐘的鐘差進(jìn)行估計(jì)。從圖中可以看出,脈沖星守時(shí)的關(guān)鍵,在于如何利用脈沖到達(dá)時(shí)間,實(shí)現(xiàn)原子鐘鐘差的準(zhǔn)確估計(jì)。

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

圖7給出了利用脈沖星估計(jì)原子鐘鐘差的方法示意圖。首先,利用星載原子鐘,可以測(cè)量脈沖到達(dá)時(shí)間。此外,對(duì)原子鐘的鐘差模型進(jìn)行建模。結(jié)合脈沖到達(dá)時(shí)間和原子鐘的鐘差模型,利用卡爾曼濾波算法,即可實(shí)現(xiàn)原子鐘的鐘差估計(jì)。

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

2.2 脈沖星測(cè)量模型

脈沖星的時(shí)間相位模型是實(shí)現(xiàn)脈沖星守時(shí)的基礎(chǔ),時(shí)間相位模型通常建立在太陽(yáng)系質(zhì)心(solar system barycenter,SSB)。假設(shè)由星載原子鐘記錄的脈沖到達(dá)衛(wèi)星的時(shí)間為tSC,則該脈沖到達(dá)SSB的時(shí)間為[17]

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

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

式中,φ表示脈沖相位;φ0表示tSSB,0=g(t0)時(shí)刻的脈沖相位;f(n)為tSSB,0時(shí)刻的脈沖星自轉(zhuǎn)頻率的n階導(dǎo)數(shù)。利用式(10)所示的時(shí)間相位模型,即可預(yù)測(cè)脈沖星信號(hào)到達(dá)SSB的時(shí)間SSB。

令

因此,脈沖到達(dá)時(shí)間的預(yù)測(cè)值與測(cè)量值之差即反映了原子鐘的鐘差。因此,脈沖星測(cè)量方程可表示為

2.3 脈沖到達(dá)時(shí)間精度分析

脈沖到達(dá)時(shí)間的精度可通過(guò)以下公式表示[18]

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

式中,Δt是觀測(cè)持續(xù)時(shí)間,λp和λn分別是脈沖信號(hào)和背景噪聲的平均流量密度,A是探測(cè)器有效面積。

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

3 守時(shí)性能分析

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

圖8(a)給出了脈沖星守時(shí)系統(tǒng)的仿真計(jì)算結(jié)果,可以看出,最初用最小二乘法解算的結(jié)果鐘差較大,之后脈沖星守時(shí)系統(tǒng)精度逐漸收斂至1μs左右。圖8(b)給出了利用批量最小二乘的脈沖星守時(shí)系統(tǒng)的鐘差,可以看出,該守時(shí)系統(tǒng)的鐘差保持在6μs以?xún)?nèi),守時(shí)精度低于卡爾曼濾波方法。這是因?yàn)?在批量最小二乘算法中未考慮鐘差模型的過(guò)程噪聲,導(dǎo)致基于批量最小二乘的脈沖星守時(shí)系統(tǒng)鐘差大于基于卡爾曼濾波的脈沖星守時(shí)系統(tǒng)鐘差。

圖8 脈沖星守時(shí)系統(tǒng)仿真結(jié)果Fig.8 Simulation results of pulsar timing system

對(duì)基于卡爾曼濾波和批量最小二乘的脈沖星守時(shí)系統(tǒng)進(jìn)行了500次蒙特卡羅仿真計(jì)算,得到基于卡爾曼濾波的守時(shí)系統(tǒng)精度為0.78μs,基于批量最小二乘的守時(shí)系統(tǒng)精度為1.86μs。可以看出,基于卡爾曼濾波的脈沖星守時(shí)系統(tǒng)精度優(yōu)于基于批量最小二乘的脈沖星守時(shí)系統(tǒng)。因此,分析并在鐘差估計(jì)過(guò)程中考慮星載原子鐘模型的過(guò)程噪聲,對(duì)提高守時(shí)系統(tǒng)的性能是有意義的。

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

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

此外,需要說(shuō)明的是,在本文的仿真分析中,脈沖星觀測(cè)策略及需要的探測(cè)器有效面積均是根據(jù)理論公式計(jì)算的,與實(shí)際空間環(huán)境下真實(shí)的脈沖星觀測(cè)性能可能存在一定差距。在后續(xù)的脈沖星守時(shí)方法研究中,應(yīng)該結(jié)合X 射線探測(cè)器的實(shí)際觀測(cè)能力或國(guó)內(nèi)外脈沖星實(shí)測(cè)數(shù)據(jù)開(kāi)展分析與研究。

4 結(jié)論

本文利用實(shí)測(cè)星載原子鐘鐘差數(shù)據(jù)和脈沖星觀測(cè)數(shù)據(jù),開(kāi)展了脈沖星守時(shí)方法研究,得到如下結(jié)論:

1)通過(guò)對(duì)比脈沖星計(jì)時(shí)殘差和星載原子鐘鐘差的功率譜,可以看出脈沖星的長(zhǎng)期性能優(yōu)于星載原子鐘。

2)由于包含系統(tǒng)誤差的影響,星載原子鐘的鐘差變化劇烈,G07號(hào)衛(wèi)星的原子鐘在約3年的時(shí)間內(nèi),鐘差變化了接近800μs。而脈沖星的計(jì)時(shí)殘差幾乎不存在系統(tǒng)誤差,PSR B1937+21在接近3年的時(shí)間內(nèi)計(jì)時(shí)殘差始終不超過(guò)15μs。因此,雖然PSR B1937+21 的隨機(jī)性誤差在33 天后才優(yōu)于G07號(hào)衛(wèi)星的原子鐘,但是由于星載原子鐘存在系統(tǒng)性誤差,即使在較短的時(shí)間內(nèi),也可以通過(guò)脈沖星來(lái)估計(jì)原子鐘的鐘差。

3)基于實(shí)測(cè)星載原子鐘鐘差數(shù)據(jù),開(kāi)展了脈沖星守時(shí)性能分析。計(jì)算結(jié)果表明,若脈沖星的脈沖到達(dá)時(shí)間解算精度為1μs/30 d,則原子鐘鐘差估計(jì)精度可達(dá)到優(yōu)于1μs的水平,初步驗(yàn)證了脈沖星守時(shí)系統(tǒng)的可行性。