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

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

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

2)(陜西省空間超限探測重點實驗室,西安 710126)

3)(鵬城實驗室,深圳 518000)

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

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

1 引言

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

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

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

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

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

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

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

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

首先利用脈沖星星歷,太陽系星歷與航天器狀態(tài)信息,將航天器處模擬時段轉換到SSB 處的模擬時段;在SSB處,借助脈沖星標準輪廓仿真SSB 處的光子到達時間序列,并利用相位模型將計算對應的到達相位序列,根據光子到達航天器處時的相位與到達SSB 處時的相位關系,此相位序列即為航天器處的光子到達相位序列;推導得到航天器處光子到達時間與相位關系,即可根據航天器處觀測到的光子相位序列得到航天器處的光子到達時間序列.

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

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

2.1.1 SSB 處光子到達相位序列的生成

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

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

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

其中,m(t)為累積強度函數,可表示為

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

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

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

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

其中,?C為時鐘校正,?P為視差延遲,?R為Roemer延遲,?E為Einstein延遲,?S為Shapiro延遲,?D為太陽系色散延遲,?A為大氣延遲.由于時間轉換項完整的表達式非常復雜,本文采用簡化公式:

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

2.2 航天器處光子到達時間與相位的關系

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

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

在(10)式中,記 ??SC=?SC(t)??SC(t0),?t=(t ?t0),則可根據附錄的推導結果給出光子到達時間與相位的關系.總結附錄A的推導結果,?t=存在多種可能,可將 ?t的解概括為

2.3 fSC,SC,SC的計算

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

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

通過計算可得到

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

3 仿真驗證與分析

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

3.1 累積輪廓證明

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

表1 初始軌道參數Table 1.Initial orbit parameters.

表2 脈沖星參數Table 2.Pulsar parameters.

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

采用Pearson 相關系數評價觀測輪廓與標準輪廓的重合度,表達式為

其中sk和wk分別為觀測脈沖輪廓與標準脈沖輪廓.表3 給出了觀測時長為106s 時的相關系數.

根據圖2的輪廓對比及表3 給出的相關系數,證明由本文方法仿真的航天器處的光子序列可恢復得到標準脈沖輪廓,證明所提航天器處光子序列仿真方法是正確的.

圖2 觀測輪廓與標準輪廓的對比(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 相關系數Table 3.Pearson coefficients for observation profiles.

3.2 與實測數據的對比

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

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

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

圖4 不同觀測時間下校正后的累積脈沖輪廓(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.

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

表5 相位差計算結果Table 5.Phase difference calculation results.

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

圖5 Pearson 相關系數隨觀測時間的變化曲線(a)實測數據實驗結果;(b)模擬數據實驗結果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 頻率與頻率導數的檢驗

3.3.1 頻率的檢驗

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

圖6 給出了3 顆脈沖星的實驗結果,對200 次實驗得到的自轉頻率的估計值求均方根誤差,計算公式為

圖6 自轉頻率搜索實驗結果(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為數 據中點 時刻的 自轉頻率,為對fmid估計的均方根誤差,N=200.經計算,得到脈沖星PSR B0531+21,PSR B1821–24與PSR B1937+21的自轉頻率估計值的均方根誤差分別為1.08×10–5,4.67× 10–5與1.29× 10–4,即模擬的自轉頻率誤差在10–5與10–4量級,證明所提模擬方法的正確.

3.3.2 頻率一階導數的檢驗

對脈沖星PSR B0531+21的頻率一階導數進行檢驗,模擬時間跨度為5 d,每天模擬6 個數據包,數據包之間時間間隔為14400 s,單個數據包觀測時長為3000 s,共30 個數據包.對每個數據包進行頻率搜索,可得到該數據包中點時刻的頻率,即可表示為

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

圖7 頻率一階導數檢驗結果(a)模擬數據包得到的待擬合數據;(b)多次重復實驗結果Fig.7.First derivative of rotation frequency test results:(a)Data to be fitted from simulated data package;(b)repeated experimental results.

圖8 頻率一階導數檢驗結果(a)模擬數據包得到的待擬合數據;(b)多次重復實驗結果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 頻率二階導數的檢驗

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

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

圖9 頻率二階導數檢驗結果(a)模擬數據包得到的待擬合數據;(b)多次重復實驗結果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)勢.為證明這點,選取PSR B0531+21 與PSR B1937+21 進行仿真效率的實驗驗證,脈沖星參數如表2 所示,探測器面積為1 m2,觀測時間為480 s.將PSR B0531+21的初始流量設置為0.5 photons/s,且依次增大流量,分別用方法1、方法3 及本文方法仿真航天器處的光子信號,并統(tǒng)計仿真所用時間.

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

圖10 方法1 與本文方法在不同流量下的仿真速度對比(a)仿真時間統(tǒng)計結果;(b)時間平均值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.

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

圖11 方法2 與本文方法在不同頻率下的仿真速度對比(a)仿真時間統(tǒng)計結果;(b)時間平均值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 顆脈沖星參數下,用4 種方法仿真600 s 時長的航天器處的光子到達時間序列,進行600 次蒙特卡洛仿真,計算所用時間的平均值,結果如表6 所示.

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

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

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

3.4.2 方法精度比較

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

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

圖12 兩種方法在不同觀測時間下的Pearson 相關系數對比(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的仿真結果.根據圖12,觀測時間較短時,兩種方法的Pearson 相關系數非常接近,隨著觀測的時間延長,方法3的Pearson 相關系數達到最大值后逐漸減小,而由本文所提方法得到的Pearson 相關系數隨觀測時間的延長而增大,說明觀測時間增大時,本文所提方法的精度明顯高于方法3的精度.

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

4 結論

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

附錄A 航天器處的時間相位關系推導

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

當d=?6??SC=0時,且 ?t1=0,?t2,?t3滿足: