蘇劍宇，方海燕，，高敬敬，趙良

1．西安電子科技大學(xué) 空間科學(xué)與技術(shù)學(xué)院，西安 710126 2．北京臨近空間飛行器系統(tǒng)工程研究所 空間物理重點實驗室，北京 100101

X射線脈沖星導(dǎo)航（X-ray Pulsar-based Navigation，XPNAV）是一種新型的自主導(dǎo)航技術(shù)，可為航天器提供位置、時間等導(dǎo)航信息［1-4］，實現(xiàn)航天器高精度自主導(dǎo)航。目前，國內(nèi)外相繼開展了X射線脈沖星導(dǎo)航試驗，如中國空間實驗室天宮二號（Tiangong-2， TG-2）的γ暴偏振探測科學(xué)實驗［5］，中國首顆X射線脈沖星導(dǎo)航試驗衛(wèi)星（XPNAV-1）在軌開展的X射線脈沖星的探測與脈沖星導(dǎo)航體制的驗證［6-7］，以及國內(nèi)首顆空間X射線天文硬X射線調(diào)制望遠鏡衛(wèi)星（Hard X-ray Modulation Telescope，HXMT） Insight-HXMT的脈沖星定軌精度驗證實驗［8］。如美國NICER（Neutron star Interior Composition Ex?plorer）項 目 的SEXTANT （Station Explorer X-ray Timing And Navigation Technology）搭載國際空間站（International Space Station， ISS）開展的定軌精度驗證工作［9］?？梢?，X射線脈沖星觀測與導(dǎo)航驗證工作正進入蓬勃發(fā)展的時期。XPNAV的基本原理可描述為［10］：在脈沖星觀測周期內(nèi)，航天器上安裝的X射線探測器會記錄到一串光子到達時間（Time of Arrivals，TOA），利用觀測的光子TOA，通過一定的算法提取出在這段觀測時間內(nèi)的某一時刻處航天器所接收的脈沖星信號的相位以及多普勒頻率，由于航天器在任意時刻觀測的脈沖信號相位和頻率都可以用該時刻的位置、速度以及太陽系質(zhì)心（Solar System Barycenter， SSB）處的脈沖星信號模型準(zhǔn)確表示，因此以估計得到的航天器在某一時刻的觀測信號相位和多普勒頻率作為導(dǎo)航觀測量，并將其表示為有關(guān)航天器在該時刻位置和速度的關(guān)系式作為導(dǎo)航觀測方程，獲得任意時刻的觀測相位即可得到相應(yīng)時刻的位置信息。高精度的觀測脈沖相位估計是整個導(dǎo)航的難點，為此，國內(nèi)外學(xué)者也相繼提出許多脈沖相位估計算法［11-18］。

在上述導(dǎo)航算法中，為獲得高精度的位置與速度估計，需要實現(xiàn)觀測相位和頻率的高精度估計。文獻［19］提出在星際轉(zhuǎn)移軌道段可將航天器的運動軌跡簡化為勻加速模型，基于此建模了航天器軌道狀態(tài)誤差與其所引起的脈沖TOA偏差之間的關(guān)系，并對導(dǎo)航觀測方程進行了修正以消除TOA偏差的影響。但線性化過程或產(chǎn)生高階截斷誤差，且該截斷誤差隨脈沖星觀測周期的增大而增大。所以通常為了將軌道位置誤差進行線性化，選擇壓縮脈沖星觀測周期或?qū)⒂^測周期分段。如文獻［20］提出將觀測時段劃分為短時段，在每個較短時段內(nèi)利用針對恒定多普勒頻移模型的相位估計算法估計相位，再借助相位跟蹤的方法提高單段相位估計精度并獲得動態(tài)的多普勒頻移估計。文獻［21］提出了由輪廓形變來求解速度誤差引起的多普勒頻移的新思路，定義了輪廓熵以衡量累積輪廓的不同形變程度，推導(dǎo)了存在速度誤差時，頻率偏差與輪廓熵之間的關(guān)系，利用最優(yōu)化方法通過使熵最小（輪廓最為清晰）從而求解出多普勒頻移。該方法本質(zhì)上仍需要假定在觀測時段內(nèi)的觀測頻率恒定，即航天器的速度恒定，所以觀測周期不能過大。Anderson和Pines［22］研究了基于鎖相環(huán)的脈沖相位跟蹤方法，并給出了針對小流量脈沖星的相位跟蹤方法。文獻［23］將較短觀測時段內(nèi)的預(yù)報位置誤差建模為隨時間近似線性變化，且將這種線性變化的位置誤差對光子到達時間修正的影響建模為恒定的相位偏差加上恒定的多普勒頻移，在此基礎(chǔ)上，提出了一種基于相位和多普勒頻移聯(lián)合觀測的X射線脈沖星導(dǎo)航方法，并進行了實驗驗證與分析。以上算法在進行觀測相位與多普勒頻率估計時，未對脈沖星觀測周期的選取進行深入研究。而最優(yōu)觀測周期的確定對脈沖星導(dǎo)航計算及觀測計劃制定等具有重要的意義，本文對該問題進行研究與討論，以期為觀測方案的制定提供參考。

1 觀測相位估計誤差

最優(yōu)脈沖星觀測周期是為得到高精度的觀測相位估計，所以將觀測相位估計誤差最小時對應(yīng)的脈沖星觀測周期稱為最優(yōu)脈沖星觀測周期。航天器在(t0，tf)的脈沖星觀測周期內(nèi)，任意時刻的觀測相位［10］為

根據(jù)式（1），對?sc(t)的估計誤差取決于對 δr的估計誤差。任意時刻的位置誤差δr為

式中：T為脈沖星觀測周期長度，T=tf?t0；δr0為航天器初始時刻t0的位置誤差。則δr引起的相位誤差δ?為

式中：δ?0為對應(yīng)t0時刻的初始預(yù)測相位。將式（3）展開到2階，表示為

根據(jù)上述分析，對δ?的估計將同時存在方差與偏差，采用均方誤差準(zhǔn)則，并用mse(δ?)表示δ?的估計均方誤差，則均方誤差mse(δ?)為

式中：var(δ?)為δ?的估計方差；b2為δ?的估計偏差，。方差是由θ的估計方差引起的，所以求出估計方差，即可得到使mse(δ?)最小的T的表達式。

根據(jù)文獻［12］，利用探測器探測到的光子到達時間序列，使用最大似然估計可實現(xiàn)對脈沖相位延遲的估計，利用最大似然估計法對θ的估計可表示為

式中：λs、λb分別為脈沖星源的流量與背景流量；h(?)為相位?∈[0，1)的脈沖星的標(biāo)準(zhǔn)輪廓，滿足=1，h(n+?)=h(?)（n為整數(shù)）［12］。Emadzadeh證明了利用式（6）的最大似然估計對脈沖相位延遲的估計方差即為克拉美羅界（Cramer-Rao Lower Bound， CRLB）。所以下面推導(dǎo)對θ估計的CRLB。

2 最優(yōu)脈沖星觀測周期確定

2.1 θ的CRLB推導(dǎo)

任意時刻的相位誤差δ?為

將θ的CRLB表示［13］為

式中：I(θ)為Fisher信息矩陣，根據(jù)Emadzadeh和Speyer［13］的結(jié)論，I(θ)的第i行第j列元素Iij的表達式為

其中：

由Iij的表達式可得Iij=Iji，所以I為對稱矩陣，且

為確定I(θ)矩陣，需要求I11、I12、I13、I23、I335個元素。Emadzadeh給出了I11、I12、I13的表達式，附錄A給出了I23、I33的推導(dǎo)過程，結(jié)果如式（A19）所示。

得到I11、I12、I13、I23、I33的表達式，即可求得的CRLB(θ)。附錄B給出了CRLB(θ)的推導(dǎo)過程。直接求解化解CRLB(θ)難以化解，這里首先化解

式中：a為I(θ)的行列式，將I11、I12、I13、I23、I33代入式（12），化解得到

所以CRLB(δ?d0)=。同理得到

所以

綜上可得

式（17）即為對參數(shù)θ估計的CRLB。

2.2 最優(yōu)觀測周期確定

得 到θ=的CRLB后，可 確 定δ?的估計誤差。根據(jù)式（17）、式（4）中δ?的估計方差var(δ?)可表示為

最優(yōu)觀測周期是指對δ?的估計均方誤差最小時的觀測周期，即所選觀測周期應(yīng)使式（19）中mse(δ?)最 小。對 式（19）求 極 值，得 到 使mse(δ?)最小的觀測周期的表達式為

至此得到了最優(yōu)脈沖星觀測周期。需指出，式（20）所定義的最優(yōu)觀測周期是在式（1）所示的航天器處任意時刻觀測相位模型下推導(dǎo)，以該觀測相位估計誤差最小為準(zhǔn)則得到。

3 仿真實驗驗證

為了驗證式（20）的正確性，對同一顆脈沖星在不同觀測周期下進行仿真實驗，并給出仿真的mse(δ?)隨觀測周期變化的曲線，通過對曲線進行擬合，得到最小mse(δ?)對應(yīng)的觀測周期，并與利用式（20）所預(yù)測的最優(yōu)觀測周期進行對比。

3.1 θ估計方差的仿真

脈沖星觀測信號的仿真是后續(xù)實驗的基礎(chǔ)，所以首先需要仿真航天器在脈沖星觀測周期內(nèi)的脈沖星信號。實驗中采用的軌道根數(shù)如表1所示，所使用的脈沖星參數(shù)如表2所示。光子TOA仿真方法參考文獻［24］。

表1 軌道參數(shù)Table 1 Orbit parameters

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

利用上述脈沖星各項參數(shù)與軌道進行光子序列仿真，航天器初始的偏差與噪聲方差設(shè)置如表3所示。

表3 初始誤差Table 3 Initial error

利用預(yù)測位置偏差對光子TOA進行時間校正，得到修正的光子TOA，然后代入式（6）的算法求?。采用利用格點搜索法［25］求式（6），可知求解過程為三維的搜索過程，的3σ搜索 范 圍 分 別 為與。

重復(fù)仿真光子TOA，利用每一次仿真的光子TOA進行一次搜索，得到θ的第i次估計值，進 行800次 重 復(fù) 試 驗，即i=1，2，…，800。然 后 利 用800次 結(jié) 果，可 求 得δ?0、δf0、δf?0的 均 方 誤 差 mse(δ?0)、mse(δf0)、mse(δf?0)。根據(jù)式（18）、式（19），可得mse(δ?)的仿真結(jié)果為

3.2 不同觀測周期下的mse(δ?)

改變觀測周期，重復(fù)3.1節(jié)中計算mse(δ?)的步驟，可得到mse(δ?)隨觀測周期變化的曲線。脈沖星PSR B1821-24仿真結(jié)果如圖1、圖2所示。

圖1給出了對θ估計的均方根誤差（Root Mean Square Error， RMSE）結(jié)果，根據(jù)2.1節(jié)的推導(dǎo)，CRLB即為對θ估計的RMSE的理論值，根據(jù)圖1的仿真結(jié)果，的估計均方根誤差隨觀測周期的增加首先逐漸靠近其CRLB，證明推導(dǎo)的θ的CRLB的正確性。但當(dāng)觀測周期繼續(xù)增大，的估計均方根誤差逐漸偏離了CRLB，這是因為對式（3）進行了截斷，且截斷誤差近似以代替，該誤差會隨觀測周期的增加而增大。

圖1 對θ估計的均方根誤差（PSR B1821-24）Fig. 1 Root mean square error of estimation of θ（PSR B1821-24）

圖2（a）為mse(δ?)隨觀測周期變化的曲線，其橫縱坐標(biāo)都進行了對數(shù)運算。根據(jù)仿真結(jié)果，mse(δ?)隨觀測周期先減小后增大，且存在使mse(δ?)最小的點，在該點附近擬合，得到圖2（b），根據(jù)擬合曲線可明顯看出極小值點的存在，并可得到mse(δ?)最小時對應(yīng)的觀測周期為620 s，利用式（20）得到仿真條件下的最優(yōu)觀測周期Topt=673 s，與實驗結(jié)果相差53 s，可見預(yù)測最優(yōu)觀測周期在實際最優(yōu)觀測周期附近，證明了該公式的正確性。

圖2 mse(δ?)隨觀測周期變化曲線（PSR B1821-24）Fig. 2 Curve of mse(δ?) with observation period （PSR B1821-24）

同理，圖3、圖4給出了脈沖星PSR B0531+21的仿真結(jié)果。圖3為θ估計的均方根誤差結(jié)果，可看出，δ?0、δf0、的估計均方根誤差隨觀測周期的增加首先逐漸靠近其CRLB，證明推導(dǎo)的θ的CRLB的正確性。當(dāng)觀測周期繼續(xù)增大，由于截斷誤差的影響，δ?0、δf0、的估計均方根誤差逐漸偏離了CRLB。

圖3 對θ估計的均方根誤差（PSR B0531+21）Fig. 3 Root mean square error of estimation of θ（PSR B1821-24）

圖4（a）為仿真得到的mse(δ?)隨觀測周期變化的曲線，其橫縱坐標(biāo)都進行了對數(shù)運算。根據(jù)仿真結(jié)果，mse(δ?)隨觀測周期先減小后增大，且存在使mse(δ?)最小的點；圖4（b）為在最小值附近擬合的結(jié)果，曲線的最小值即為仿真得到的最優(yōu)觀測周期。

圖4 mse(δ?)隨觀測周期變化曲線（PSR B0531+21）Fig. 4 Curve of mse(δ?) with observation period （PSR B0531+21）

根據(jù)擬合曲線可明顯看出極小值點的存在，并可得到mse(δ?)最小時對應(yīng)的觀測周期為435 s，利用式（20）得到仿真條件下的最優(yōu)觀測周期391 s，與實驗結(jié)果相差44 s，可見預(yù)測最優(yōu)觀測周期在實際最優(yōu)觀測周期附近，證明了該公式的正確性。

為進一步驗證公式的正確性，使用NICER的實測數(shù)據(jù)進行了驗證。數(shù)據(jù)包為“ni1013010142_0mpu7_cl.enents.gz”［26］。該 數(shù)據(jù)包為PSR B0531+21的實測數(shù)據(jù)，包含了40279473個光子TOA。表2給出了脈沖星參數(shù)信息，軌道根數(shù)與表1數(shù)據(jù)相同。利用該實測數(shù)據(jù)對θ=進行搜索，并改變觀測周期，重復(fù)上述步驟，可得到mse(δ?)隨觀測周期變化的曲線，圖5（a）為mse(δ?)隨觀測周期變化的曲線，其橫縱坐標(biāo)都進行了對數(shù)運算。通過在mse(δ?)最小的點附近擬合，得到擬合結(jié)果如圖5（b）所示。

圖5 mse(δ?)隨觀測周期變化曲線（PSR B0531+21實測數(shù)據(jù)）Fig. 5 Curve of mse(δ?) with observation period （Ob?servational data of PSR B0531+21）

根據(jù)擬合曲線可明顯看出極小值點的存在，并可得到mse(δ?)最小時對應(yīng)的觀測周期為460 s，利用式（20）得到仿最優(yōu)觀測周期Topt=410 s，與實驗結(jié)果相差50 s，可見預(yù)測最優(yōu)觀測周期在實際最優(yōu)觀測周期附近，且將實測數(shù)據(jù)的驗證結(jié)果與圖4的仿真結(jié)果對比，可看出在相同條件下實測數(shù)據(jù)與仿真數(shù)據(jù)的結(jié)果是接近的，證明了式（20）的可靠性。

使用NICER對PSR B1821-24的觀測數(shù)據(jù)，進一步進行實驗驗證。數(shù)據(jù)包名稱為“ni0070010102_0mpu7_cl_t1”“ni0070010103_0mpu7_cl_t1”“ni0070010104_0mpu7_cl_t1”［26］，共19971個光子TOA，結(jié)果如圖6所示。

從圖6（a）可看出，mse(δ?)隨觀測周期的變化規(guī)律與圖5相同，根據(jù)擬合曲線可明顯看出極小值點的存在，并可得到mse(δ?)最小時對應(yīng)的觀測周期為650 s，利用式（20）得到仿真條件下的最優(yōu)觀測周期714 s，與實驗結(jié)果相差64 s，即預(yù)測最優(yōu)觀測周期在實際最優(yōu)觀測周期附近，證明了式（20）的正確性。且將實測數(shù)據(jù)的驗證結(jié)果與圖2的仿真結(jié)果對比，可看出在相同條件下實測數(shù)據(jù)與仿真數(shù)據(jù)的結(jié)果接近，證明了所提公式的可靠性。

圖6 mse(δ?)隨觀測周期變化曲線（PSR B1821-24實測數(shù)據(jù)）Fig. 6 Curve of mse(δ?) with observation period （Ob?servational data of PSR B1821-24）

4 結(jié)論

在X射線脈沖星導(dǎo)航中，觀測周期的選取是獲得高精度觀測相位估計的關(guān)鍵，本文對該問題進行了研究，得到以下主要結(jié)論：

1）以航天器處觀測相位估計的均方誤差最小為準(zhǔn)則，推導(dǎo)了脈沖星觀測周期與觀測相位估計的均方誤差的關(guān)系，給出了最優(yōu)脈沖星觀測周期的近似計算公式。

2）設(shè)計了仿真驗證實驗，仿真結(jié)果表明給出的最優(yōu)觀測周期計算公式對脈沖星PSR B1821-24與PSR B0531+21的預(yù)測誤差約為64、44 s，而由脈沖星PSR B0531+21的實測數(shù)據(jù)得到的預(yù)測誤差為50 s，即由該公式給出的最優(yōu)觀測周期與實驗中對觀測相位估計的均方誤差最小時對應(yīng)的觀測周期接近，證明所給公式的正確性。

附錄A：

定義

根據(jù)Emadzadeh和Speyer［12］的結(jié)論，I11、I12、I13可表示為

式中：

以下將推導(dǎo)I23、I33的表達式。

將觀測周期分成N+1段，即

圖A1 觀測周期分段Fig. A1 Observation period segmentation

當(dāng)n=0，Pn=0，有

將觀測時間分成N+1段，將式（A4）的積分表示為圖A1中每個小區(qū)間的積分求和，I23可表示為

根據(jù)圖A1，在第n+1個區(qū)間內(nèi)，脈沖星周期為Pn+1，則積分變量t可表示時刻對應(yīng)的相位?(t)為

將式（A6）代入式（A5），以?為新的積分變量，并利用h的周期性，可化解得到

由于N?1，所以式（A7）中占據(jù)主導(dǎo)地位，即

同理，I33的結(jié)果為

下面對I23、I33的表達式進行簡化處理。I23的求和項的多項式展開式中，去除包含P0的項，共有項非零項，定義Pˉ23為

相應(yīng)的，I23、I33可表示為

式中：q23為形如的非零項。所以

綜上，由于N?1，所以可近似認(rèn)為=0，即。同理可證。

將I11、I12、I13表 示 為I23、I33的 形 式，則I11、I12、I13、I23、I33表達式為

附錄B：

I?1(θ)的求解：根據(jù)逆矩陣運算法則，I?1(θ)可表示為

則θ=中的每一個元素的CRLB可表示為