，，，，，，，，，，

1. 錢學(xué)森空間技術(shù)實(shí)驗(yàn)室，北京 100094 2. 中國(guó)空間技術(shù)研究院，北京 100094

脈沖星是存在于宇宙中的一種高速運(yùn)轉(zhuǎn)的中子星，因自身極其穩(wěn)定的周期性而被譽(yù)為宇宙中最穩(wěn)定的天文時(shí)鐘，尤其是毫秒脈沖星的周期穩(wěn)定度可達(dá)到10-19～10-21s/s[1，2]。X射線脈沖星導(dǎo)航(X-ray Pulsar-based Navigation，XPNAV)是利用脈沖星輻射的X射線進(jìn)行導(dǎo)航[3-5]的一種新概念導(dǎo)航技術(shù)，它是利用X射線頻段的計(jì)時(shí)觀測(cè)數(shù)據(jù)，通過(guò)差分方程自主確定航天器的導(dǎo)航參數(shù)。X射線脈沖星導(dǎo)航系統(tǒng)對(duì)于未來(lái)導(dǎo)航系統(tǒng)的發(fā)展具有重要的研究意義，近年來(lái)備受國(guó)內(nèi)外學(xué)者的廣泛關(guān)注[6-8]。

2016年11月10日，中國(guó)首顆X射線脈沖星導(dǎo)航試驗(yàn)衛(wèi)星-1(X-ray Pulsar-based Navigation-1,XPNAV-1)成功發(fā)射。XPNAV-1衛(wèi)星工作在500 km的LEO軌道，是一顆質(zhì)量為270 kg的小衛(wèi)星。該衛(wèi)星搭載了兩種不同體制的X射線探測(cè)器——掠入射Wolter-I聚焦型X射線探測(cè)器(簡(jiǎn)稱Wolter-I探測(cè)器)和準(zhǔn)直型微通道板(Microchannel Plates，MCP)探測(cè)器(簡(jiǎn)稱MCP探測(cè)器)。XPNAV-1衛(wèi)星主要任務(wù)是在軌開(kāi)展X射線脈沖星的探測(cè)并進(jìn)行脈沖星導(dǎo)航體制的驗(yàn)證[9，10]。在前期的觀測(cè)中，選擇了蟹狀星云(Crab)脈沖星(PSR B0531+21)[11]進(jìn)行觀測(cè)，目前已經(jīng)完成了Wolter-I探測(cè)器對(duì)Crab脈沖星的多次觀測(cè)，得到了Crab脈沖星輻射的一系列的X射線光子。

在X射線脈沖星導(dǎo)航系統(tǒng)中，脈沖星的脈沖到達(dá)時(shí)間(TOA)作為導(dǎo)航系統(tǒng)的基本觀測(cè)量，其測(cè)量精度是后續(xù)的導(dǎo)航定位精度的決定性因素，因此對(duì)脈沖TOA進(jìn)行高精度的求解是XPNAV系統(tǒng)中一個(gè)非常關(guān)鍵的問(wèn)題[12]。多重信號(hào)子空間分類(Multiple Signal Classification，MUSIC)算法是陣列信號(hào)處理領(lǐng)域的一種高精度算法，主要是用于信號(hào)波達(dá)方向的估計(jì)，為了進(jìn)行高精度的到達(dá)時(shí)間(Time Of Arrival，TOA)估計(jì)，本文提出了利用MUSIC方法，通過(guò)對(duì)接收數(shù)據(jù)的預(yù)處理，完成脈沖星導(dǎo)航系統(tǒng)中的TOA估計(jì)。為了驗(yàn)證方法的性能，對(duì)提出的方法進(jìn)行試驗(yàn)仿真，即對(duì)Wolter-I探測(cè)器觀測(cè)得到的X射線光子進(jìn)行了光子折疊輪廓處理得到折疊輪廓，利用本文的方法，成功得到了折疊輪廓的TOA估計(jì)值，并且估計(jì)精度與傳統(tǒng)互相關(guān)法相當(dāng)。

1 脈沖輪廓

標(biāo)準(zhǔn)輪廓是脈沖星的一個(gè)重要特征，它是通過(guò)對(duì)脈沖星進(jìn)行長(zhǎng)時(shí)間觀測(cè)，將得到的X射線光子TOA通過(guò)光子歷元折疊方法得到的。不同脈沖星的標(biāo)準(zhǔn)脈沖輪廓不同，因此可以將脈沖星的標(biāo)準(zhǔn)脈沖輪廓作為XPNAV系統(tǒng)的一個(gè)必備數(shù)據(jù)。

圖1是Crab脈沖星的標(biāo)準(zhǔn)輪廓[13]，表明了Crab脈沖星的一個(gè)重要特征，它是利用長(zhǎng)時(shí)間多次觀測(cè)的數(shù)據(jù)建立的。在XPNAV-1的觀測(cè)中，截取了一段時(shí)間內(nèi)的觀測(cè)數(shù)據(jù)進(jìn)行光子歷元折疊。觀測(cè)時(shí)間選為UTC 57 727.0約簡(jiǎn)儒略日(Modified Julian Day，MJD)到UTC 57 812.0 MJD，在0.5～10 keV頻段內(nèi)Wolter-I探測(cè)器共接收到4 853 983個(gè)光子，得到了119段PSR B0531+21的光子數(shù)據(jù)，將每段光子數(shù)據(jù)進(jìn)行輪廓折疊[14，15]，可以得到一系列的脈沖折疊輪廓。隨機(jī)選取8段時(shí)間內(nèi)光子數(shù)據(jù)進(jìn)行脈沖折疊輪廓，得到圖2。

圖1 Crab脈沖星標(biāo)準(zhǔn)輪廓Fig.1 Standard profile of Crab pulsar

圖2 Crab脈沖星折疊輪廓Fig.2 Epoch folding profiles of Crab pulsar

對(duì)比圖1和圖2，圖2八段時(shí)間的觀測(cè)輪廓與標(biāo)準(zhǔn)脈沖輪廓形狀相同，但是在橫坐標(biāo)(即脈沖相位)存在著時(shí)間的左右偏離，時(shí)間偏離值稱為脈沖TOA。脈沖TOA與相位之間的關(guān)系為TOA=φT0，T0為脈沖周期。

因此對(duì)于Crab脈沖星，可以定義其標(biāo)準(zhǔn)輪廓為s，觀測(cè)得到的折疊輪廓為r，兩者之間的關(guān)系可以建立為：

r(φi)=βs(φi-φ0)+ω(φi)

(1)

式中：β為幅度因子；φ0為r相對(duì)于s的相位偏離，即為待估計(jì)的TOA；φi(i=1,2,…,Nb)為第i個(gè)格對(duì)應(yīng)的相位；Nb為在輪廓折疊過(guò)程中，將一個(gè)相位周期均勻劃分的格數(shù)；ω(φi)為等效的折疊噪聲，在折疊格中的光子數(shù)多于20個(gè)時(shí)，該噪聲可以認(rèn)為是均值為零的高斯噪聲[16]。

2 脈沖TOA估計(jì)

在X射線脈沖星導(dǎo)航中，利用航天器上的X射線探測(cè)器探測(cè)到脈沖星X射線光子數(shù)據(jù)，通過(guò)對(duì)比脈沖星的同一個(gè)脈沖到達(dá)航天器與基準(zhǔn)點(diǎn)之間的相位差，得到后續(xù)的導(dǎo)航參數(shù)。基準(zhǔn)點(diǎn)可以設(shè)置在太陽(yáng)系質(zhì)心，通過(guò)脈沖星的計(jì)時(shí)觀測(cè)模型可以得到基準(zhǔn)點(diǎn)的TOA，而航天器處脈沖TOA需要利用互相關(guān)法、本文方法等獲取。

2.1 互相關(guān)法

互相關(guān)法是TOA估計(jì)的常用方法，其基本原理在于對(duì)互相關(guān)函數(shù)最大值所對(duì)應(yīng)的自變量的求解。

在用于XPNAV系統(tǒng)中進(jìn)行脈沖TOA估計(jì)時(shí)，首先需要計(jì)算標(biāo)準(zhǔn)輪廓與折疊輪廓的互相關(guān)函數(shù)，即

(2)

由自相關(guān)函數(shù)的性質(zhì)：

|Rss(φ-φ0)|≤Rss(0)

(3)

因此，在Rsr(φ)在φ=φ0處有一個(gè)峰值，峰值對(duì)應(yīng)處即為脈沖TOA的估計(jì)：

(4)

式中：arg{·}表示取函數(shù)的自變量；max[·]表示取函數(shù)的最大值?；ハ嚓P(guān)函數(shù)的峰值點(diǎn)對(duì)應(yīng)的延遲量即為脈沖TOA的估計(jì)。

2.2 本文方法

互相關(guān)法的計(jì)算簡(jiǎn)單，但XPNAV系統(tǒng)中，估計(jì)精度受噪聲及信號(hào)形式的影響較大，因此針對(duì)XPNAV系統(tǒng)，本文提出了利用MUSIC算法對(duì)TOA進(jìn)行估計(jì)。

2.2.1 MUSIC算法原理

MUSIC算法[17，18]具有“超分辨、高精度”的性能，基本原理是將接收的數(shù)據(jù)進(jìn)行特征值分解，分解后得到兩個(gè)相互正交的空間，利用兩個(gè)空間的正交性計(jì)算出尖銳的譜峰，譜峰所對(duì)應(yīng)的參數(shù)即為需要估計(jì)的參數(shù)。

首先可以將MUSIC算法的接收信號(hào)模型建立為：

X(t)=A(θ)S(t)+N(t)

(5)

計(jì)算協(xié)方差矩陣:

R=E[XXH]=AE[SSH]AH+σ2I=

ARSAH+σ2I

(6)

對(duì)協(xié)方差矩陣進(jìn)行特征值分解，可以得到特征值和對(duì)應(yīng)的特征向量。將特征值由大到小排序λ1>λ2>…>λM，將特征值分為兩部分，D個(gè)較大的特征值λ1,…,λD對(duì)應(yīng)著信號(hào)，相應(yīng)的特征向量u1,…,uD構(gòu)成了信號(hào)子空間US=[u1,…,uD]；(M-D)個(gè)較小的特征值對(duì)應(yīng)著噪聲，相應(yīng)的特征向量uD+1,…,uM構(gòu)成了噪聲子空間UN=[uD+1,…,uM]，假設(shè)噪聲為高斯白噪聲，那么有λD+1≈…≈λM≈0。

由MUSIC算法的性質(zhì)[17]可知，US與UN正交，并且導(dǎo)向矢量陣A(θ)與US張成了同一空間，因此有：

(7)

式(7)等號(hào)右側(cè)為一個(gè)零向量，表明陣列的導(dǎo)向矢量陣與噪聲子空間是相互正交的。但是在實(shí)際環(huán)境中，式(7)的右側(cè)并不是一個(gè)嚴(yán)格的零向量，而是近似為零向量。因此，MUSIC算法可以利用此性質(zhì)來(lái)進(jìn)行參數(shù)估計(jì)，即譜估計(jì)公式為：

(8)

最后，搜索PMUSIC(θk)的極大值點(diǎn)即為對(duì)應(yīng)的參數(shù)值。

2.2.2 TOA估計(jì)

在波達(dá)方向(Direction-Of-Arrival，DOA)估計(jì)領(lǐng)域，MUSIC算法因其“高精度”性能備受學(xué)者的青睞。為了利用這個(gè)性質(zhì)，本文將該算法引入到脈沖星導(dǎo)航系統(tǒng)中的TOA估計(jì)。

觀察MUSIC算法的步驟可知，MUSIC算法進(jìn)行DOA估計(jì)時(shí)，利用的是同一個(gè)接收時(shí)間下的不同陣元之間的接收數(shù)據(jù)進(jìn)行信號(hào)參數(shù)的估計(jì)。而利用X射線光子探測(cè)器(Wolter-I探測(cè)器)接收到的光子數(shù)據(jù)并不是同一個(gè)接收時(shí)間下的不同陣元之間的接收數(shù)據(jù)，而是在同一個(gè)陣元的不同接收時(shí)間下的光子數(shù)據(jù)，并不符合MUSIC算法的信號(hào)模型，為此需要將光子數(shù)據(jù)進(jìn)行整合處理。

利用一個(gè)X射線光子探測(cè)器(Wolter-I探測(cè)器)，在不同時(shí)間下得到光子到達(dá)時(shí)間的標(biāo)記，然后進(jìn)行脈沖輪廓折疊，可以得到式(1)的形式，將其寫成矩陣形式為：

r= [r(φ1)r(φ2) …r(φNb)]T=

A(φ0)·β+ω

(9)

式中：r為觀測(cè)的折疊輪廓矢量；A(φ0)=[s(φ1-φ0)s(φ2-φ0)…s(φNb-φ0)]T為導(dǎo)向矢量，是一個(gè)與φ0(即脈沖TOA)有關(guān)的矢量；ω=[ω(φ1)ω(φ2)…ω(φNb)]T為噪聲矢量。

式(9)與式(5)具有相同的組成形式，因此可以考慮通過(guò)MUSIC算法來(lái)求解TOA。具體步驟為：

1)計(jì)算r的協(xié)方差矩陣Rr=E[rrH]；

2)對(duì)Rr進(jìn)行特征值分解，UN=[u2,…,uNb]；

3)計(jì)算譜函數(shù)：

(10)

其中，搜索向量構(gòu)造為：

A(φi)=[s(φ1-φi)s(φ2-φi)s(φNb-φi)]T

(11)

4)對(duì)f=[f(φ1),f(φ2),…,f(φNb)]的最大值進(jìn)行搜索，最大值所對(duì)應(yīng)的φk即為脈沖相位的估計(jì)，即：

(12)

3 TOA估計(jì)結(jié)果及分析

3.1 仿真及結(jié)果分析

首先對(duì)算法的有效性進(jìn)行仿真，選取PSR B0531+21(即Crab脈沖星)為觀測(cè)脈沖星，Crab脈沖星的周期為0.033 1 s，背景流量為5×10-3ph/s/cm2，來(lái)自脈沖星的流量為1.54 ph/s/cm2(流量參數(shù)統(tǒng)計(jì)頻段為2 keV至10 keV內(nèi)單位時(shí)間、單位面積探測(cè)器接收的光子數(shù))。

(1)仿真1 TOA估計(jì)譜

選取Crab脈沖星為觀測(cè)脈沖星，探測(cè)器的面積為0.5 m2，時(shí)間分辨率為100 μs，觀測(cè)周期起點(diǎn)為57 727.0 MJD，觀測(cè)周期長(zhǎng)為0.5 s，起點(diǎn)脈沖相位設(shè)置為φ0=0.7，分別利用互相關(guān)和本文方法進(jìn)行一次脈沖相位估計(jì)，并進(jìn)行歸一化處理，得到圖3。

圖3 TOA估計(jì)譜Fig.3 TOA estimation spectrum

由圖3可知，本文方法的估計(jì)譜要比互相關(guān)方法的估計(jì)譜尖銳，對(duì)譜峰進(jìn)行搜索，得到了互相關(guān)法與本文方法的TOA估計(jì)值分別為0.683 1和0.691 2，即本文方法的估計(jì)值較為準(zhǔn)確。

(2)仿真2 不同觀測(cè)時(shí)間下TOA估計(jì)誤差

為了進(jìn)一步驗(yàn)證算法的性能，觀測(cè)時(shí)間選取0.1 s，0.2 s，0.5 s，1 s，5 s，10 s，20 s，50 s，在不同的觀測(cè)時(shí)間下分別進(jìn)行100次Monte-Carlo仿真，計(jì)算TOA估計(jì)的均方根誤差，其他參數(shù)與仿真1相同。

圖4為不同觀測(cè)時(shí)間下的TOA估計(jì)誤差，可以看到，觀測(cè)時(shí)間越長(zhǎng)，估計(jì)誤差越小，即估計(jì)精度變高，而在相同的觀測(cè)時(shí)間下，本文方法的估計(jì)精度要高于互相關(guān)法。

圖4 不同觀測(cè)時(shí)間下的TOA估計(jì)誤差Fig.4 RMSE in different observation time

3.2 實(shí)測(cè)數(shù)據(jù)分析

圖2是隨機(jī)選取的8段光子數(shù)據(jù)進(jìn)行脈沖折疊輪廓圖，下面對(duì)這8段數(shù)據(jù)分別利用互相關(guān)法和本文方法進(jìn)行脈沖相位的估計(jì)，將兩種方法得到的估計(jì)譜進(jìn)行歸一化處理并進(jìn)行比較，得到圖5。

由圖5可知，互相關(guān)法和本文方法均能得到TOA估計(jì)的譜峰，即兩種方法都成功進(jìn)行了TOA估計(jì)。對(duì)比兩者的估計(jì)譜，可以明顯看出前者的譜峰要比后者的尖銳，因此本文算法的估計(jì)結(jié)果較為精確。

最后統(tǒng)計(jì)8段數(shù)據(jù)的光子數(shù)和利用兩種方法估計(jì)得到的脈沖相位，如表1所示。由表1的結(jié)果可以看出，本文方法與互相關(guān)法的估計(jì)結(jié)果相當(dāng)。

序號(hào)時(shí)間起點(diǎn)/MJD光子數(shù)/個(gè)TOA估計(jì)結(jié)果/周互相關(guān)法MUSIC方法157727.0416180.2197260.217773257732.6597500.4990230.500976357738.3392550.1865230.184570457749.4403250.9365230.940429557758.1417720.4111320.413085657780.0458320.6962890.696289757793.5421220.0869140.086913857810.0444280.3818360.381835

4 結(jié)束語(yǔ)

為了進(jìn)行TOA估計(jì)，本文提出了利用MSUIC算法的高精度性質(zhì)進(jìn)行脈沖星導(dǎo)航系統(tǒng)中的TOA估計(jì)，并且對(duì)互相關(guān)法和本文方法進(jìn)行了對(duì)比。仿真結(jié)果表明，本文方法的TOA估計(jì)的譜峰要比互相關(guān)法尖銳，在相同的情況下，本文方法的估計(jì)精度要高于互相關(guān)法。在XPNAV-1衛(wèi)星前期的觀測(cè)中，選取了光子流量較強(qiáng)的Crab脈沖星，衛(wèi)星所搭載的探測(cè)器成功接收到了Crab脈沖星的光子數(shù)據(jù)，本文中選取了57 727.0 MJD到57812.0 MJD內(nèi)0.5～10 keV頻段內(nèi)的光子觀測(cè)數(shù)據(jù)進(jìn)行處理，得到了脈沖輪廓折疊。XPNAV-1衛(wèi)星的成功發(fā)射，對(duì)建立X射線脈沖星導(dǎo)航系統(tǒng)具有突破性的意義。在XPNAV-1衛(wèi)星后續(xù)的任務(wù)中，我們也規(guī)劃了對(duì)其他光子強(qiáng)度較弱的脈沖星的觀測(cè)。

References)

[1] 帥平.脈沖星:宇宙航行的燈塔[D]. 北京：國(guó)防工業(yè)出版社，2016.

SHUAI P. Pulsar are lighthouses for space flight[M]. Beijing： National Defense Industry Press，2016(in Chinese).

[2] HEWISH A，BELL S J，PILKINGTON D H，et al. Obser-vation of a rapidly pulsating radio source[J]. Nature，1968，217：709-713.

[3] HANSON J E. Principles of X-ray navigation[D]. CA：Stanford University，1996：1-32.

[4] SHEIKH S I. The use of variable celestial X-ray sources for spacecraft navigation[D]. MD：Maryland University，2005：375-403.

[5] 帥平，陳紹龍，吳一帆，等. X射線脈沖星導(dǎo)航原理[J].宇航學(xué)報(bào)，2007，28(6): 1538-1543.

SHUAI P，CHEN S L，WU Y F，et al. Navigation principles using X-ray pulsars[J]. Journal of Astronautics，2007，28(6)：1538-1543(in Chinese).

[6] The ATNF Pulsar Catalogue[EB/OL].[2017-03-17].http:∥www.atnf.csiro.au/research/pulsar/psrcat.

[7] XIONG Z，QIAO L，KIU J Y，et al. GEO satellite autonomous navigation using X-ray pulsar navigation and GNSS measurements[J]. International Journal of Innovative Computing，Information and Control，2012，8(5): 2965-2977.

[8] ANDESON K D，PINES D J. Method of pulse phase tracking for X-ray pulsar based spacecraft navigation using low flux pulsars[C]. Space Ops. Conferences，Pasadena，CA，5-9 May，2014.

[9] 黃良偉，帥平，林晴晴，等. X 射線脈沖星導(dǎo)航標(biāo)稱數(shù)據(jù)庫(kù)構(gòu)建[J].中國(guó)空間科學(xué)技術(shù)，2015，32(3)：66-74.

HUANG L W，SHUAI P，LIN Q Q，et al. Research of the nominal database for X-ray pulsar navigation[J]. Chinese Space Science andTechnology，2015，35(3)：66-74(in Chinese).

[10] 黃良偉，帥平，張新源，等. 脈沖星導(dǎo)航試驗(yàn)衛(wèi)星時(shí)間數(shù)據(jù)分析與脈沖輪廓恢復(fù)[J]. 中國(guó)空間科學(xué)技術(shù)，2017，37(3)：1-10.

HUANG L W，SHUAI P，ZHANG X Y，et al. XPNAV-1 Satellite timing data analysis and pulse profile recovery[J]. Chinese Space Science and Technology，2017，37(3)：1-10(in Chinese).

[11] WINTERNITZ L M，HASSOUNEH M，MITCELL J W，et al. X-ray navigation algorithms and testbed for SEXTANT[C]. IEEE Aerospace conference，2015:1-14.

[12] 林晴晴，帥平，黃良偉，等. 一種高精度的X射線脈沖星導(dǎo)航TOA估計(jì)方法[J].宇航學(xué)報(bào)，2016，37(6):604-701.

LIN Q Q，SHUAI P，HUANG L W，et al. A high precision TOA estimation method for X-ray pulsar-based navigation system[J]. Journal of Astronautic，2016，37(6):604-701(in Chinese).

[13] MPIfR EPN Pulsar Profiles Database，2017，available: http:∥rian.kharkov.ua/decameter/EPN/browser.html.

[14] EMADZADEH A A，SPEYER J L. X-ray pulsar-based relative navigation using epoch folding[J]. IEEE Transactions on Aerospace and Electronic Systems，2011，47(4): 2317-2328.

[15] EMADZADEH A A，SPEYER J L. On modeling and pulse phase estimation of X-ray pulsars[J]. IEEE Transactions on Signal Processing，2010，58(9)：4484-4495.

[16] HANSON J E，SHEIKH S I，GRAVEN P，et al. Noise analysis for X-ray navigation systems[C]. 2008 IEEE/ION Position Location and Navigation Symposium，Monterey，CA，May，2008.

[17] SCHMIDT R O. Multiple emitter location and spectrum estimation[J]. IEEE Transactions on Signal Processing，1986，34(3): 276-280.

[18] JIN H，SWAMY M N S，AHMAD M O. Efficient appli-cation of MUSIC algorithm under the coexistence of far-field and near-field sources[J]. IEEE Transactions on Signal Processing，2012，58(1): 2066-2070.