X射線(xiàn)脈沖星導(dǎo)航的最優(yōu)觀測(cè)周期確定

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

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

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

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

1 觀測(cè)相位估計(jì)誤差

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

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

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

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

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

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

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

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

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

2.1 θ的CRLB推導(dǎo)

任意時(shí)刻的相位誤差δ?為

將θ的CRLB表示［13］為

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

其中：

由Iij的表達(dá)式可得Iij=Iji，所以I為對(duì)稱(chēng)矩陣，且

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

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

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

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

所以

綜上可得

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

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

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

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

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

3 仿真實(shí)驗(yàn)驗(yàn)證

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

3.1 θ估計(jì)方差的仿真

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

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

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

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

表3 初始誤差Table 3 Initial error

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

附錄A：

定義

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

式中：

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

將觀測(cè)周期分成N+1段，即

圖A1 觀測(cè)周期分段Fig. A1 Observation period segmentation

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

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

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

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

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

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

下面對(duì)I23、I33的表達(dá)式進(jìn)行簡(jiǎn)化處理。I23的求和項(xiàng)的多項(xiàng)式展開(kāi)式中，去除包含P0的項(xiàng)，共有項(xiàng)非零項(xiàng)，定義Pˉ23為

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

式中：q23為形如的非零項(xiàng)。所以

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

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

附錄B：

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

則θ=中的每一個(gè)元素的CRLB可表示為