航天器處X 射線脈沖星觀測信號模擬方法*

蘇劍宇 方海燕? 包為民 孫海峰2)3) 趙良

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

2)(陜西省空間超限探測重點(diǎn)實(shí)驗(yàn)室,西安 710126)

3)(鵬城實(shí)驗(yàn)室,深圳 518000)

4)(北京臨近空間飛行器系統(tǒng)工程研究所空間物理重點(diǎn)實(shí)驗(yàn)室,北京 100101)

航天器處X 射線脈沖星觀測信號模擬對脈沖星信號處理方法及導(dǎo)航方案的驗(yàn)證具有重要意義.通過建立航天器處光子到達(dá)時(shí)間與相位之間的關(guān)系,可顯著提高脈沖星信號模擬算法的效率.目前建立的航天器處光子到達(dá)時(shí)間與相位之間關(guān)系模型并未考慮接收信號頻率變化,模擬精度較低.針對這一問題,本文建立了考慮頻率一、二階導(dǎo)數(shù)時(shí)的光子到達(dá)時(shí)間與相位關(guān)系模型,并基于該模型給出航天器處觀測信號模擬方法,提高了脈沖星觀測信號模擬算法精度.仿真實(shí)驗(yàn)證明,相比于迭代法或?qū)崟r(shí)計(jì)算航天器處光子到達(dá)速率函數(shù)的方法,利用推導(dǎo)的光子到達(dá)時(shí)間與相位關(guān)系模型,仿真速度最高可提高3 個(gè)數(shù)量級,保證了信號模擬算法的高效率;且相比于未考慮頻率變化的情況,由包含頻率一、二階導(dǎo)數(shù)的模型得到的Pearson 相關(guān)系數(shù)最高可提高350.0%,顯著提高了計(jì)算精度.

1 引言

X 射線脈沖星導(dǎo)航(X-ray pulsar-based navigation,XPNAV)是一種新型的自主導(dǎo)航技術(shù),可為在太陽系內(nèi)外飛行的航天器提供位置、時(shí)間等導(dǎo)航信息[1?2],實(shí)現(xiàn)航天器高精度自主導(dǎo)航.目前,國內(nèi)外相繼開展了X 射線脈沖星導(dǎo)航試驗(yàn).如中國空間實(shí)驗(yàn)室天宮二號(Tiangong-2,TG-2)的γ暴偏振探測科學(xué)實(shí)驗(yàn)[3],中國首顆X射線脈沖星導(dǎo)航試驗(yàn)衛(wèi)星(X-ray pulsar navigation-I,XPNAV-1)在軌開展的X射線脈沖星的探測與脈沖星導(dǎo)航體制的驗(yàn)證[4?6],以及國內(nèi)首顆空間X 射線天文硬X 射線調(diào)制望遠(yuǎn)鏡衛(wèi)星(hard X-ray modulation telescope,HXMT)Insight-HXMT的脈沖星定軌精度驗(yàn)證實(shí)驗(yàn)[7?8].國外的如美國NICER(neutron star interior composition explorer,NICER)項(xiàng)目的SEXTANT(station explorer Xray timing and navigation technology,SEXTANT)搭載國際空間站(international space station,ISS)開展的定軌精度驗(yàn)證工作[9?10].

由于X 射線無法穿透大氣層,在地面無法觀測到脈沖星輻射的X 射線信號,而巡天觀測和空間飛行試驗(yàn)耗時(shí)長、成本高,且算法的評估和驗(yàn)證難以完全通過飛行試驗(yàn)完成,現(xiàn)階段對X 射線脈沖星導(dǎo)航的相關(guān)研究需要借助地面模擬產(chǎn)生的X 射線脈沖星信號.地面模擬方法可分為半物理方法和數(shù)值方法,其中半物理模擬方法用脈沖星標(biāo)準(zhǔn)輪廓對光子信號的強(qiáng)度進(jìn)行調(diào)制來模擬真實(shí)的脈沖星輻射過程,主要包括機(jī)械轉(zhuǎn)盤[11]和電控光源[12]兩種形式.機(jī)械轉(zhuǎn)盤的模擬方式結(jié)構(gòu)較復(fù)雜,脈沖星輪廓特性和頻率特性的模擬精度相對較低,且模擬不同的脈沖星信號需要更換形狀不同的轉(zhuǎn)盤,靈活性較差.電控光源的方式可以解決上述問題,根據(jù)采用的光源不同,可分為X射線光源[13,14]和可見光源模擬[15].

相比于半物理模擬方法,數(shù)值模擬方法具有簡單、高效、靈活的優(yōu)點(diǎn),且是半物理模擬實(shí)現(xiàn)的算法基礎(chǔ).Emadzadeh 提出了基于反函數(shù)法的脈沖星信號模擬方法[16].Zhang等[17]提出了一種利用高斯分布快速生成X 射線脈沖星信號的方法,計(jì)算效率上有顯著提高.Jin等[18]提出了一種基于蒙特卡羅的X 射線脈沖星信號模擬方法,改進(jìn)了概率統(tǒng)計(jì)方法,克服了非齊次泊松分布的約束,具有仿真時(shí)間短、效率高、精度高的優(yōu)點(diǎn).然而,這些方法只能模擬到達(dá)太陽系質(zhì)心(solar system barycenter,SSB)處或靜止航天器處的X 射線脈沖星信號,而無法模擬信號在動態(tài)環(huán)境中,如在運(yùn)動的航天器中探測到的X 射線脈沖星信號,而實(shí)際中探測器隨航天器運(yùn)動,所以需要模擬運(yùn)動的航天器處的X射線脈沖星觀測信號.

目前,航天器處的X 射線脈沖星信號的數(shù)值模擬主要有兩種方式.一種是先模擬SSB 處的光子到達(dá)時(shí)間序列,然后通過迭代的方式得到航天器處的光子到達(dá)時(shí)間序列[19],稱之為方法1.該方法借助航天器的軌道參數(shù)信息,通過對時(shí)間轉(zhuǎn)換公式的迭代求解,計(jì)算航天器處的光子到達(dá)時(shí)間(time of arrivals,TOAs).由于求解過程主要是利用反復(fù)的迭代,所以對于大流量脈沖星(如脈沖星PSR B0531+21),光子數(shù)越多,針對每個(gè)光子進(jìn)行上述的迭代運(yùn)算將非常耗時(shí),效率非常低.另一種方式基于同一個(gè)光子到達(dá)航天器與SSB 處的相位相等的原理,直接在航天器處模擬光子序列,稱之為方法2.該方法首先根據(jù)SSB 處的脈沖星信號到達(dá)模型和航天器的實(shí)時(shí)位置及速度信息,建立航天器處的光子到達(dá)速率函數(shù),再利用尺度變換法[20]直接產(chǎn)生航天器處的光子TOAs[21].該方法準(zhǔn)確建模了脈沖星信號的頻率緩變,避免了方法1的復(fù)雜迭代過程.但由于航天器處的光子到達(dá)速率函數(shù)隨航天器的狀態(tài)變化而變化,需要實(shí)時(shí)計(jì)算更新,這對于頻率高的脈沖星而言,計(jì)算量龐大,仿真過程同樣非常耗時(shí).針對上述兩種航天器處的X 射線脈沖星信號的模擬方法存在的問題,Su等[22]提出了一種新的X 射線脈沖星信號的模擬方法,文中稱之為方法3.該方法通過公式推導(dǎo),給出了光子到達(dá)航天器的時(shí)間與相位之間的關(guān)系式,再結(jié)合光子到達(dá)航天器與SSB 處的相位相等的原理,只需模擬SSB 處的光子到達(dá)相位序列,即可得到航天器處的光子到達(dá)時(shí)間序列,提高了航天器處脈沖星觀測信號的仿真效率.方法3 雖然提高了航天器處脈沖星觀測信號的仿真效率,但在推導(dǎo)光子到達(dá)航天器處的時(shí)間與相位之間的關(guān)系時(shí),并未考慮接收信號頻率的變化,即在給出的時(shí)間與相位關(guān)系中,并未將頻率的變化建模,導(dǎo)致觀測時(shí)間延長時(shí),時(shí)間與相位關(guān)系的誤差增大,仿真的光子TOAs 誤差增大,且對于頻率越大的脈沖星,仿真的光子TOAs誤差越大.所以該方法只適用較短觀測時(shí)間,觀測時(shí)間較長時(shí),精度無法滿足要求.

綜合上述方法,為在保證航天器處脈沖星觀測信號的仿真效率的同時(shí),提高光子TOAs的計(jì)算精度,本文延續(xù)方法3的思路,推導(dǎo)了考慮接收脈沖星信號頻率一、二階導(dǎo)數(shù)時(shí)的光子到達(dá)時(shí)間與相位之間的關(guān)系模型,并基于該模型給出航天器處觀測信號模擬方法,提高了光子TOAs的計(jì)算精度,并通過與上述3 種方法的對比實(shí)驗(yàn),證明所提方法的正確性與有效性.

2 航天器處X 射線脈沖星信號仿真算法

本文所提航天器處X 射線脈沖星信號仿真方法如圖1 所示.

圖1 航天器處X 射線脈沖星信號模擬方法流程圖Fig.1.Flow chart of X-ray pulsar signal simulation method at spacecraft.

首先利用脈沖星星歷,太陽系星歷與航天器狀態(tài)信息,將航天器處模擬時(shí)段轉(zhuǎn)換到SSB 處的模擬時(shí)段;在SSB處,借助脈沖星標(biāo)準(zhǔn)輪廓仿真SSB 處的光子到達(dá)時(shí)間序列,并利用相位模型將計(jì)算對應(yīng)的到達(dá)相位序列,根據(jù)光子到達(dá)航天器處時(shí)的相位與到達(dá)SSB 處時(shí)的相位關(guān)系,此相位序列即為航天器處的光子到達(dá)相位序列;推導(dǎo)得到航天器處光子到達(dá)時(shí)間與相位關(guān)系,即可根據(jù)航天器處觀測到的光子相位序列得到航天器處的光子到達(dá)時(shí)間序列.

根據(jù)上述過程,推導(dǎo)航天器處光子到達(dá)時(shí)間與相位關(guān)系是仿真算法的關(guān)鍵,本文在推導(dǎo)該關(guān)系式時(shí)考慮了接收脈沖星信號頻率一、二階導(dǎo)數(shù),提高了光子到達(dá)時(shí)間的計(jì)算精度,下面給出所提仿真算法的詳細(xì)推導(dǎo)過程.

2.1 航天器處光子到達(dá)相位序列的生成

2.1.1 SSB 處光子到達(dá)相位序列的生成

為仿真航天器處的光子TOAs,首先需仿真SSB 處的光子TOAs.脈沖星輻射的X 射線光子信號在漫長的傳播過程中不斷衰減,到達(dá)SSB 時(shí)已成為微弱的單光子流,可將其建模為強(qiáng)度為λ(t)的非齊次泊松過程[23],即在(t0,t)時(shí)間段內(nèi),接收到的光子數(shù)Nt0,t等于K的概率P的表達(dá)式為

式中,λ(t)為SSB 處的光子到達(dá)速率函數(shù),且:

利用脈沖星標(biāo)準(zhǔn)輪廓與相位模型,可獲得非齊次泊松過程的強(qiáng)度函數(shù),本文采用尺度變換法[25]產(chǎn)生非齊次泊松過程到達(dá)時(shí)間序列.非齊次泊松事件發(fā)生時(shí)刻s1,s2,···,sn與參數(shù)為1的齊次泊松過程事件發(fā)生時(shí)刻z1,z2,···,zn的關(guān)系為

其中,m(t)為累積強(qiáng)度函數(shù),可表示為

所以首先產(chǎn)生參數(shù)為1的齊次泊松過程事件發(fā)生時(shí)刻,然后利用(4)式的變換關(guān)系即可獲得非齊次泊松事件發(fā)生時(shí)刻.通過仿真非齊次泊松過程的到達(dá)時(shí)間序列就獲得了SSB 處的光子TOAs,利用(5)式,將每個(gè)光子的到達(dá)時(shí)間轉(zhuǎn)化為相位,可得到SSB 處的光子到達(dá)相位序列.

2.1.2 航天器處光子到達(dá)相位序列的生成

任意t時(shí)刻,航天器處觀測到的光子相位?SC與該光子在SSB 處的相位?之間關(guān)系可表示為

其中,τ(t)為時(shí)間轉(zhuǎn)換模型,時(shí)間轉(zhuǎn)換模型描述了光子在傳播過程中經(jīng)過航天器和SSB的時(shí)間差,這一過程包含以下時(shí)間校正項(xiàng)[26,27]

其中,?C為時(shí)鐘校正,?P為視差延遲,?R為Roemer延遲,?E為Einstein延遲,?S為Shapiro延遲,?D為太陽系色散延遲,?A為大氣延遲.由于時(shí)間轉(zhuǎn)換項(xiàng)完整的表達(dá)式非常復(fù)雜,本文采用簡化公式:

式中各變量均在太陽系質(zhì)心慣性坐標(biāo)系下表示,該模型忽略了視差、脈沖星自行以及除太陽外的其他天體的Shapiro 延遲等項(xiàng),轉(zhuǎn)換精度為5—8 μs[21].第一項(xiàng)為幾何延遲,第二項(xiàng)為愛因斯坦延遲,第三項(xiàng)為太陽Shapiro 延遲.其中τ(t)為時(shí)間轉(zhuǎn)換各項(xiàng)之和,n=(nx,ny,nz)為脈沖星在太陽系質(zhì)心坐標(biāo)系下的單位方向矢量,rSC/SSB為航天器相對與太陽系質(zhì)心的位置矢量,rSC/E為航天器相對于地球的位置矢量,vE為地球的速度矢量,rSC/Sun為航天器相對于太陽的位置矢量,c為光速,μS為太陽引力常數(shù).

2.2 航天器處光子到達(dá)時(shí)間與相位的關(guān)系

利用方法3的思路,如果推導(dǎo)出航天器處光子到達(dá)時(shí)間與相位的關(guān)系,即可直接由SSB 處的光子到達(dá)相位序列得到航天器光子到達(dá)時(shí)間序列.假設(shè)航天器處任意時(shí)刻t探測到的觀測相位為?SC(t),相對于觀測起始時(shí)刻t0(對應(yīng)相位?SC(t0)),有如下關(guān)系[28]

其中,fSC(τ)為τ時(shí)刻航天器觀測到的脈沖星頻率,將(9)式進(jìn)行二階泰勒展開,可得到

在(10)式中,記 ??SC=?SC(t)??SC(t0),?t=(t ?t0),則可根據(jù)附錄的推導(dǎo)結(jié)果給出光子到達(dá)時(shí)間與相位的關(guān)系.總結(jié)附錄A的推導(dǎo)結(jié)果,?t=存在多種可能,可將 ?t的解概括為

2.3 fSC,SC,SC的計(jì)算

在得到航天器處光子相位序列后,利用(11)式給出的 ?t與??SC的關(guān)系式,可直接得到航天器處光子TOAs.(11)式中公式選擇及計(jì)算需要使用等參數(shù),下面推導(dǎo)的計(jì)算公式.

對(6)式求導(dǎo),可得到航天器處的頻率表達(dá)式:

通過計(jì)算可得到

其中,vSC/SSB為航天器相對于SSB的速度矢量,vSC/E為航天器相對于地球的速度矢量,aE為地球的加速度矢量.將(14)式各項(xiàng)代入(13)式可得到.(13)式中對t求導(dǎo),可得到為

3 仿真驗(yàn)證與分析

本節(jié)從4 個(gè)方面對模擬方法進(jìn)行驗(yàn)證:1)通過模擬光子TOAs 得到的累積脈沖輪廓與標(biāo)準(zhǔn)脈沖輪廓的對比,證明模擬光子TOAs的正確性;2)通過與實(shí)測數(shù)據(jù)的對比實(shí)驗(yàn),進(jìn)一步驗(yàn)證模擬光子TOAs的正確性;3)通過頻率與頻率導(dǎo)數(shù)的檢驗(yàn),驗(yàn)證所提方法考慮了頻率變化,進(jìn)一步驗(yàn)證模擬光子TOAs的正確性;4)相同條件下,比較所提方法與方法1、方法2 與方法3的仿真效率與精度,證明所提方法相對于現(xiàn)有方法的優(yōu)越性.

3.1 累積輪廓證明

為證明所提方法的正確性,利用仿真光子TOAs恢復(fù)觀測脈沖輪廓,比較觀測脈沖輪廓與仿真中使用的標(biāo)準(zhǔn)脈沖輪廓的相關(guān)系數(shù).仿真中使用的航天器軌道參數(shù)如表1 所示.選取3 顆脈沖星,其參數(shù)如表2 所示,觀測時(shí)長分別為104s,105s 與106s,仿真航天器處接收到的光子到達(dá)時(shí)間序列.

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

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

利用(8)式將航天器處的光子到達(dá)時(shí)間序列進(jìn)行時(shí)間校正,將校正后的光子序列進(jìn)行歷元折疊得到觀測輪廓.觀測輪廓與標(biāo)準(zhǔn)輪廓的對比如圖2所示.根據(jù)圖2的輪廓對比,隨著觀測時(shí)長的增加,利用本文方法仿真的航天器處的光子序列可恢復(fù)標(biāo)準(zhǔn)脈沖輪廓,證明所提航天器處光子序列仿真方法是正確的.

采用Pearson 相關(guān)系數(shù)評價(jià)觀測輪廓與標(biāo)準(zhǔn)輪廓的重合度,表達(dá)式為

其中sk和wk分別為觀測脈沖輪廓與標(biāo)準(zhǔn)脈沖輪廓.表3 給出了觀測時(shí)長為106s 時(shí)的相關(guān)系數(shù).

根據(jù)圖2的輪廓對比及表3 給出的相關(guān)系數(shù),證明由本文方法仿真的航天器處的光子序列可恢復(fù)得到標(biāo)準(zhǔn)脈沖輪廓,證明所提航天器處光子序列仿真方法是正確的.

圖2 觀測輪廓與標(biāo)準(zhǔn)輪廓的對比(a)PSR B0531+21(104s);(b)PSR B1821–24(104s);(c)PSR B1937+21(104s);(d)PSR B0531+21(105s);(e)PSR B1821–24(105s);(f)PSR B1937+21(105s);(g)PSR B0531+21(106s);(h)PSR B1821–24(106s);(i)PSR B1937+21(106s)Fig.2.Comparison between observed profile and standard profile:(a)PSR B0531+21(104s);(b)PSR B1821–24(104s);(c)PSR B1937+21(104s);(d)PSR B0531+21(105s);(e)PSR B1821–24(105s);(f)PSR B1937+21(105s);(g)PSR B0531+21(106s);(h)PSR B1821–24(106s);(i)PSR B1937+21(106s).

表3 觀測輪廓的Pearson 相關(guān)系數(shù)Table 3.Pearson coefficients for observation profiles.

3.2 與實(shí)測數(shù)據(jù)的對比

為驗(yàn)證模擬的航天器處的X 射線脈沖星信號具有同實(shí)測數(shù)據(jù)一樣的物理特性,本節(jié)將仿真結(jié)果與RXTE(Rossi X-ray timing Explorer)衛(wèi)星觀測數(shù)據(jù)進(jìn)行對比,采用RXTE 衛(wèi)星對脈沖星PSR B0531+21的觀測數(shù)據(jù),數(shù)據(jù)包觀測號為95802-01-16-02,其相關(guān)參數(shù)見表4,起始觀測時(shí)間使用簡約儒略日(modified Julian Day,MJD)表示.對未進(jìn)行時(shí)間校正的實(shí)測數(shù)據(jù)進(jìn)行周期折疊,得到累積脈沖輪廓,結(jié)果如圖3 所示.

圖3 不同觀測時(shí)間下未進(jìn)行時(shí)間校正的累積脈沖輪廓(a)10 s;(b)40 s;(c)90 sFig.3.Cumulative pulse profiles without time correction for different observation times:(a)10 s;(b)40 s;(c)90 s.

圖3(a)—(c)分別是利用觀測時(shí)長為10 s,40 s和90 s的實(shí)測數(shù)據(jù)與模擬數(shù)據(jù)直接進(jìn)行周期折疊得到的累積脈沖輪廓,根據(jù)圖中的實(shí)驗(yàn)結(jié)果,如果不進(jìn)行時(shí)間校正,由于航天器運(yùn)動產(chǎn)生的多普勒效應(yīng),接收脈沖星信號頻率變化快,采用固定頻率進(jìn)行周期折疊得到的累積脈沖輪廓出現(xiàn)峰值偏移現(xiàn)象,且隨觀測時(shí)長的增大,這種偏移現(xiàn)象越嚴(yán)重.通過對比,由本文模擬的航天器處X 射線脈沖星信號得到的累積脈沖輪廓隨觀測時(shí)長的變化規(guī)律與實(shí)測數(shù)據(jù)的結(jié)果相同,說明模擬信號具有與實(shí)測數(shù)據(jù)相同的物理特性.

同理,將觀測數(shù)據(jù)進(jìn)行時(shí)間校正后,得到累積脈沖輪廓如圖4(a)—(c)所示,從圖4 可看出,經(jīng)過時(shí)間校正后獲得的累積脈沖輪廓隨著觀測時(shí)間的延長,輪廓越來越清晰,且仿真數(shù)據(jù)得到的結(jié)果與實(shí)測數(shù)據(jù)相同.經(jīng)過由仿真數(shù)據(jù)與實(shí)測數(shù)據(jù)在校正前后得到的累積脈沖輪廓對比,證明空間多普勒效應(yīng)、脈沖星自轉(zhuǎn)特性等都正確地反映在光子序列上,即證明本文模擬方法的正確.

圖4 不同觀測時(shí)間下校正后的累積脈沖輪廓(a)10 s;(b)40 s;(c)90 sFig.4.Cumulative pulse profiles with time correction for different observation times:(a)10 s;(b)40 s;(c)90 s.

從上述結(jié)果可看出由實(shí)測數(shù)據(jù)與模擬數(shù)據(jù)得到的累積輪廓的變化規(guī)律是相同的,為定量的說明這一點(diǎn),可計(jì)算由模擬數(shù)據(jù)得到的累積輪廓與實(shí)測數(shù)據(jù)得到的累積輪廓之間的相位偏差.如表5 所示,給出了觀測時(shí)間從10 s 到90 s的由模擬數(shù)據(jù)得到的累積輪廓與實(shí)測數(shù)據(jù)得到的累積輪廓之間的相位偏差.可計(jì)算得到,未校正前,相位偏差最大值為0.0043;校正后的相位偏差最大值為0.00035,充分說明本文方法模擬的數(shù)據(jù)與實(shí)測數(shù)據(jù)接近,進(jìn)一步證明了方法的正確.

表5 相位差計(jì)算結(jié)果Table 5.Phase difference calculation results.

圖5 給出了時(shí)間校正前后Pearson 相關(guān)系數(shù)隨觀測時(shí)長的變化,圖中采用雙對數(shù)坐標(biāo).圖5(a)為實(shí)測數(shù)據(jù)的仿真結(jié)果,根據(jù)實(shí)驗(yàn)結(jié)果,隨著觀測時(shí)間的延長,由未進(jìn)行時(shí)間校正的光子數(shù)據(jù)得到的脈沖輪廓與標(biāo)準(zhǔn)脈沖輪廓的相關(guān)系數(shù)先增大后減小,校正后的光子數(shù)據(jù)得到的脈沖輪廓與標(biāo)準(zhǔn)脈沖輪廓的相關(guān)系數(shù)逐漸增大且趨向于1.圖5(b)為模擬數(shù)據(jù)的實(shí)驗(yàn)結(jié)果,可看出仿真數(shù)據(jù)和實(shí)測數(shù)據(jù)具有相同的變化規(guī)律,說明模擬數(shù)據(jù)具有和實(shí)測數(shù)據(jù)相同的特性,驗(yàn)證了模擬方法的正確性與合理性.

圖5 Pearson 相關(guān)系數(shù)隨觀測時(shí)間的變化曲線(a)實(shí)測數(shù)據(jù)實(shí)驗(yàn)結(jié)果;(b)模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果Fig.5.Variation curves of Pearson correlation coefficient with observation time:(a)Experiment results for observational data;(b)experiment results for simulation data.

3.3 頻率與頻率導(dǎo)數(shù)的檢驗(yàn)

3.3.1 頻率的檢驗(yàn)

在表2的脈沖星參數(shù)條件下,多次模擬觀測數(shù)據(jù),利用模擬數(shù)據(jù)搜索自轉(zhuǎn)頻率,頻率搜索方法采用法[29],并計(jì)算自轉(zhuǎn)頻率誤差.單次模擬數(shù)據(jù)的觀測時(shí)長為3000 s,重復(fù)200 次實(shí)驗(yàn),得到結(jié)果如圖6 所示.

圖6 給出了3 顆脈沖星的實(shí)驗(yàn)結(jié)果,對200 次實(shí)驗(yàn)得到的自轉(zhuǎn)頻率的估計(jì)值求均方根誤差,計(jì)算公式為

圖6 自轉(zhuǎn)頻率搜索實(shí)驗(yàn)結(jié)果(a)PSR B0531+21;(b)PSR B1821–24;(c)PSR B1937+21Fig.6.Rotation frequency search experimental results:(a)PSR B0531+21;(b)PSR B1821–24;(c)PSR B1937+21.

其中fmid為數(shù) 據(jù)中點(diǎn) 時(shí)刻的 自轉(zhuǎn)頻率,為對fmid估計(jì)的均方根誤差,N=200.經(jīng)計(jì)算,得到脈沖星PSR B0531+21,PSR B1821–24與PSR B1937+21的自轉(zhuǎn)頻率估計(jì)值的均方根誤差分別為1.08×10–5,4.67× 10–5與1.29× 10–4,即模擬的自轉(zhuǎn)頻率誤差在10–5與10–4量級,證明所提模擬方法的正確.

3.3.2 頻率一階導(dǎo)數(shù)的檢驗(yàn)

對脈沖星PSR B0531+21的頻率一階導(dǎo)數(shù)進(jìn)行檢驗(yàn),模擬時(shí)間跨度為5 d,每天模擬6 個(gè)數(shù)據(jù)包,數(shù)據(jù)包之間時(shí)間間隔為14400 s,單個(gè)數(shù)據(jù)包觀測時(shí)長為3000 s,共30 個(gè)數(shù)據(jù)包.對每個(gè)數(shù)據(jù)包進(jìn)行頻率搜索,可得到該數(shù)據(jù)包中點(diǎn)時(shí)刻的頻率,即可表示為

同理給出了脈沖星PSR B1821–24的頻率一階導(dǎo)數(shù)檢驗(yàn)結(jié)果,模擬時(shí)間跨度為360 d,每3 天模擬1 個(gè)數(shù)據(jù)包,單個(gè)數(shù)據(jù)包觀測時(shí)長為5000 s,共120 個(gè)數(shù)據(jù)包.檢驗(yàn)結(jié)果如圖8 所示.利用200 次實(shí)驗(yàn)結(jié)果得到的的估計(jì)值求均方根誤差,得到脈沖星PSR B1821–24的頻率一階導(dǎo)數(shù)的估計(jì)值的均方根誤差為3.12× 10–15,而模擬數(shù)據(jù)中使用的的真實(shí)值為–1.73× 10–13,相對誤差約為1.8%,進(jìn)一步證明了所提模擬方法的正確.

圖7 頻率一階導(dǎo)數(shù)檢驗(yàn)結(jié)果(a)模擬數(shù)據(jù)包得到的待擬合數(shù)據(jù);(b)多次重復(fù)實(shí)驗(yàn)結(jié)果Fig.7.First derivative of rotation frequency test results:(a)Data to be fitted from simulated data package;(b)repeated experimental results.

圖8 頻率一階導(dǎo)數(shù)檢驗(yàn)結(jié)果(a)模擬數(shù)據(jù)包得到的待擬合數(shù)據(jù);(b)多次重復(fù)實(shí)驗(yàn)結(jié)果Fig.8.First derivative of rotation frequency test results:(a)Data to be fitted from simulated data package;(b)repeated experimental results.

3.3.3 頻率二階導(dǎo)數(shù)的檢驗(yàn)

類似于頻率一階導(dǎo)數(shù)的檢驗(yàn),對脈沖星PSR B0531+21的頻率二階導(dǎo)數(shù)進(jìn)行檢驗(yàn).模擬時(shí)間跨度為240 d,單個(gè)數(shù)據(jù)包觀測時(shí)長為3000 s,首先利用上述一階導(dǎo)數(shù)估計(jì)方法每3 天獲得1 個(gè)一階導(dǎo)數(shù)估計(jì)值,并將該值作為這3 天的中點(diǎn)時(shí)刻的一階導(dǎo)數(shù)的值,記為則240 d共可獲得80 個(gè)一階導(dǎo)數(shù)估計(jì)值,即可表示為

根據(jù)圖9(b)中的100 次重復(fù)實(shí)驗(yàn)結(jié)果得到的二階導(dǎo)數(shù)估計(jì)值,計(jì)算得到脈沖星PSR B0531+21的頻率二階導(dǎo)數(shù)的估計(jì)值的均方根誤差為2.86×10–22,而模擬數(shù)據(jù)中使用的真實(shí)值為8.18× 10–21,相對誤差約為3.4%,進(jìn)一步證明了所提模擬方法的正確.

圖9 頻率二階導(dǎo)數(shù)檢驗(yàn)結(jié)果(a)模擬數(shù)據(jù)包得到的待擬合數(shù)據(jù);(b)多次重復(fù)實(shí)驗(yàn)結(jié)果Fig.9.Second derivative of rotation frequency test results:(a)Data to be fitted from simulated data package;(b)repeated experimental results.

3.4 不同方法的比較

3.4.1 方法效率對比

本文所提方法與方法3 思路相近,所以所提方法針對高流量、高頻率脈沖星在仿真速度上也具有優(yōu)勢.為證明這點(diǎn),選取PSR B0531+21 與PSR B1937+21 進(jìn)行仿真效率的實(shí)驗(yàn)驗(yàn)證,脈沖星參數(shù)如表2 所示,探測器面積為1 m2,觀測時(shí)間為480 s.將PSR B0531+21的初始流量設(shè)置為0.5 photons/s,且依次增大流量,分別用方法1、方法3 及本文方法仿真航天器處的光子信號,并統(tǒng)計(jì)仿真所用時(shí)間.

在每個(gè)流量參數(shù)下,進(jìn)行600 次重復(fù)實(shí)驗(yàn),得到3 種方法在該流量條件下的箱線圖;改變流量參數(shù),重復(fù)上述過程,得到3 種方法的仿真時(shí)長隨流量變化的箱線圖,如圖10(a)所示,圖10(b)為仿真所用時(shí)間的統(tǒng)計(jì)平均值變化曲線.根據(jù)圖7,方法1的仿真時(shí)長隨脈沖星流量的增大明顯增大,而所提方法與方法3的仿真速度變化非常小,即脈沖星流量變化對所提方法的效率影響很小.

圖10 方法1 與本文方法在不同流量下的仿真速度對比(a)仿真時(shí)間統(tǒng)計(jì)結(jié)果;(b)時(shí)間平均值Fig.10.Comparison of simulation speed between method 1 and the proposed method under different flux conditions:(a)Simulation time statistics results;(b)time average.

同理,為驗(yàn)證所提方法的效率不受脈沖星頻率變化影響,將PSR B1937+21的初始頻率設(shè)置為5 Hz,依次增大其頻率,分別用方法2、方法3 與所提方法仿真航天器處的光子到達(dá)時(shí)間序列,統(tǒng)計(jì)仿真所用時(shí)間.在每個(gè)頻率參數(shù)下,進(jìn)行600 次重復(fù)實(shí)驗(yàn),得到3 種方法在該頻率條件下的箱線圖;依次增大頻率,重復(fù)該過程,得到3 種方法的仿真時(shí)長隨頻率變化的箱線圖,結(jié)果如圖11 所示.從圖11可看出,方法2的仿真時(shí)長隨脈沖星頻率的增大依次增大,而所提方法的仿真速度與方法3 接近,即脈沖星頻率變化對所提方法的效率影響很小.

圖11 方法2 與本文方法在不同頻率下的仿真速度對比(a)仿真時(shí)間統(tǒng)計(jì)結(jié)果;(b)時(shí)間平均值Fig.11.Comparison of simulation speed between method 2 and the proposed method under different frequencies:(a)Simulation time statistics results;(b)time average.

在表2 中給出的3 顆脈沖星參數(shù)下,用4 種方法仿真600 s 時(shí)長的航天器處的光子到達(dá)時(shí)間序列,進(jìn)行600 次蒙特卡洛仿真,計(jì)算所用時(shí)間的平均值,結(jié)果如表6 所示.

表6 3 種模擬方法效率測試Table 6.Efficiency test of three simulation methods.

根據(jù)表6 中的結(jié)果,所提方法仿真速度與方法3 接近,與方法1 相比,本文方法的仿真速度提高了近1 個(gè)數(shù)量級,尤其是對于大流量的脈沖星PSR B0531+21,仿真速度速度提高了2 個(gè)數(shù)量級;與方法2 相比,本文方法的仿真速度提高了近2 個(gè)數(shù)量級,尤其是對于頻率較大的PSR B1821–24與PSR B1937+21,速度提高了3 個(gè)數(shù)量級.

綜上,所提方法的效率與方法3 接近,較其他兩種方法效率更高,且脈沖星流量或頻率的變化對方法效率影響很小,可高效仿真航天器處的任意脈沖星的信號.

3.4.2 方法精度比較

相較于方法1 與方法2,方法3 在效率上具有明顯優(yōu)勢,但方法3的精度會隨著觀測時(shí)間的延長而降低.本文所提方法保留了方法3 相對于其他兩種方法的優(yōu)勢,且所提方法考慮了頻率一階導(dǎo)數(shù)與二階導(dǎo)數(shù),所以相比于方法3 精度更高.

為證明這一點(diǎn),在每個(gè)觀測時(shí)長條件下,仿真光子數(shù)據(jù),并獲得脈沖輪廓,計(jì)算Pearson 相關(guān)系數(shù).在每個(gè)觀測時(shí)長下,進(jìn)行600 次重復(fù)實(shí)驗(yàn),并求平均值;改變觀測時(shí)長,重復(fù)上述過程,得到兩種方法的Pearson 相關(guān)系數(shù)隨觀測時(shí)長變化的曲線,如圖12 所示.

圖12 兩種方法在不同觀測時(shí)間下的Pearson 相關(guān)系數(shù)對比(a)PSR B0531+21;(b)PSR B1821–24;(c)PSR B1937+21Fig.12.Comparison of Pearson correlation coefficients between two methods for different observation times:(a)PSR B0531+21;(b)PSR B1821–24;(c)PSR B1937+21.

圖12 分別給出了脈沖星PSR B0531+21,PSR B1821–24與PSR B1937+21的仿真結(jié)果.根據(jù)圖12,觀測時(shí)間較短時(shí),兩種方法的Pearson 相關(guān)系數(shù)非常接近,隨著觀測的時(shí)間延長,方法3的Pearson 相關(guān)系數(shù)達(dá)到最大值后逐漸減小,而由本文所提方法得到的Pearson 相關(guān)系數(shù)隨觀測時(shí)間的延長而增大,說明觀測時(shí)間增大時(shí),本文所提方法的精度明顯高于方法3的精度.

經(jīng)計(jì)算,觀測時(shí)間為600 s時(shí),相較于由方法3 得到的Pearson 相關(guān)系數(shù),由所提方法得到的上述3 顆脈沖星的Pearson 相關(guān)系數(shù)可分別提高43.7%,59.5%,350.0%.證明所提方法可明顯提高航天器處光子數(shù)據(jù)的仿真精度,且對于頻率越高的脈沖星(如PSR B1937+21),提高效果更顯著.

4 結(jié)論

本文推導(dǎo)了考慮航天器處接收脈沖星信號頻率一、二階導(dǎo)數(shù)時(shí)的光子到達(dá)時(shí)間與相位之間的關(guān)系模型,并基于該模型,提出了一種航天器處X 射線脈沖星觀測信號的模擬方法.該方法避免了對時(shí)間轉(zhuǎn)換公式的迭代計(jì)算及航天器處實(shí)時(shí)光子到達(dá)速率函數(shù)計(jì)算等復(fù)雜過程,保證了算法的高效率,同時(shí)由于給出的光子到達(dá)時(shí)間與相位關(guān)系考慮了頻率一、二階導(dǎo)數(shù),提高了光子TOAs的計(jì)算精度.仿真實(shí)驗(yàn)證明:1)由所提方法得到的模擬數(shù)據(jù)可恢復(fù)標(biāo)準(zhǔn)脈沖輪廓,且頻率與頻率導(dǎo)數(shù)的檢驗(yàn)證明所提方法考慮了頻率的變化,計(jì)算精度更高,充分驗(yàn)證了該方法的正確性;2)所提方法的仿真速度不受脈沖星流量或頻率變化的影響,且相比于迭代法或?qū)崟r(shí)計(jì)算航天器處光子到達(dá)速率函數(shù)的方法,仿真速度最高可提高3 個(gè)數(shù)量級,保證了信號模擬算法的高效率;3)通過脈沖輪廓Pearson 相關(guān)系數(shù)對比實(shí)驗(yàn),所提方法相比于未考慮接收信號頻率變化的情況,Pearson 相關(guān)系數(shù)最高可提高350.0%,顯著提高了計(jì)算精度.該模擬方法可應(yīng)用于X 射線脈沖星信號處理與導(dǎo)航等算法的驗(yàn)證研究.

附錄A 航天器處的時(shí)間相位關(guān)系推導(dǎo)

一元三次方程的一般形式為

(A4)式所示的導(dǎo)出型一元三次方程的根的判別式為?=(q/2)2+(p/3)3[29],代入p,q的表達(dá)式得到判別式:

1)當(dāng)p≥0時(shí),由(A5)式可知,此時(shí) ?≥0,(A3)式的根可分為兩種情況:

①?>0,(A3)式導(dǎo)出型一元三次方程存在1 個(gè)實(shí)根與兩個(gè)共軛復(fù)根[30],由于 ?t與??SC均為實(shí)數(shù),y也為實(shí)數(shù),所以只考慮實(shí)根,唯一實(shí)根為

②?=0,此時(shí)(A3)式所示的導(dǎo)出型一元三次方程為3 個(gè)實(shí)根,其中兩個(gè)為重根.記3 個(gè)根為y1,y2,y3,其表達(dá)式為[30]

由于p≥0,且?=0,可知p=q=0,所以此時(shí)y1=y2=y3=0.

2)當(dāng)p <0時(shí),此時(shí)Δ的符號不確定,可分為3 種情況:

①?>0,此時(shí)存在唯一實(shí)根,如(A6)式所示.

②?=0,方程的根如(A7)式所示,即存在2 個(gè)實(shí)根,但實(shí)際中,?t>0,??SC>0且?t應(yīng)隨 ??SC的增大 而增大,由y=可知,y也隨??SC的增大而增大,所以需要進(jìn)一步判斷(A7)式中的根隨??SC的單調(diào)性.根據(jù)(A7)式,y1,y2,y3的單調(diào)性由p決定.

根據(jù)(A4)式,當(dāng)f¨SC>0時(shí),q隨??SC增大而增大,而由?=(q/2)2+(p/3)3=0得到可知y1是p的單調(diào)遞增函數(shù),所以y1隨著 ??的增大而增大,所以y1是方程的解;同理,當(dāng)時(shí),y2,y3隨著??的增大而增大,即y2,y3是方程的解.

③?<0,此時(shí)(A7)式所示的導(dǎo)出型一元三次方程存在3 個(gè)互異的實(shí)根為[30]

即函數(shù)f(y)過點(diǎn)(b/3a,d/a),因?yàn)閒(y)的零點(diǎn)為(A3)式的根,所以可通過考察f(y)的零點(diǎn)判斷(A3)式的根.根據(jù)b=的符號,f(y)的函數(shù)圖像為如圖A1 所示的兩種情況.

圖A1 f(y)的函數(shù)圖像Fig.A1.Function images of

當(dāng)b=<0時(shí),b/(3a)<0,如 圖A1(b)所示,此時(shí)y1>y2>b/3a,可知無法根據(jù)y >b/s3a判斷方程的根,需要利用其他條件進(jìn)行判斷.根據(jù)(10)式可知,當(dāng)??SC=0,有 ?t=0,即在 ??SC=0時(shí),(A3)式的根應(yīng)當(dāng)滿足y=?t+(b/3a)=b/3a,下面根據(jù)這一條件對y1,y2進(jìn)行判斷.

注意到(A2)式中,??SC=0時(shí),即 ?t1=0為方程(A6)式的一個(gè)根,假設(shè)其他兩個(gè)根為?t2,?t3,根據(jù)根與系數(shù)之間關(guān)系,?t1,?t2,?t3滿足:

當(dāng)d=?6??SC=0時(shí),且 ?t1=0,?t2,?t3滿足: