桂先洲，黃森林，孫 晨

(國(guó)防科技大學(xué)航天科學(xué)與工程學(xué)院，湖南長(zhǎng)沙 410073)

空間探測(cè)中如何實(shí)現(xiàn)深空探測(cè)器的完全自主導(dǎo)航是一個(gè)頭等難題。目前應(yīng)用的導(dǎo)航方式主要有天文導(dǎo)航、衛(wèi)星導(dǎo)航、慣性導(dǎo)航和地基導(dǎo)航。傳統(tǒng)的天文導(dǎo)航依賴近天體進(jìn)行定位，定位精度不高，且操作復(fù)雜［1-3］。近地空間航天器可以利用衛(wèi)星導(dǎo)航系統(tǒng)實(shí)現(xiàn)較精確導(dǎo)航，但衛(wèi)星導(dǎo)航系統(tǒng)本身需要定位，自主性較差。慣性導(dǎo)航系統(tǒng)的誤差累積效應(yīng)嚴(yán)重，不能單獨(dú)進(jìn)行長(zhǎng)時(shí)間的導(dǎo)航，一般作為組合導(dǎo)航的輔助導(dǎo)航手段［4］。地基導(dǎo)航系統(tǒng)通常使用雷達(dá)測(cè)距和光學(xué)跟蹤進(jìn)行定位，但對(duì)于深空導(dǎo)航，其定位的橫向精度差，且完全依賴地面站，地面站要進(jìn)行大量的運(yùn)算［5］。X射線脈沖星是一類高速旋轉(zhuǎn)，并發(fā)射X射線的中子星，探測(cè)器接收到的X射線信號(hào)體現(xiàn)為周期性的脈沖，其周期穩(wěn)定性可以媲美現(xiàn)在精度最高的原子鐘［6］。X射線脈沖星是天然的導(dǎo)航信標(biāo)，脈沖星導(dǎo)航是目前最具潛力實(shí)現(xiàn)深空自主導(dǎo)航的技術(shù)。

目前對(duì)X射線脈沖星導(dǎo)航的研究主要集中在定位方面［7-8］，關(guān)鍵之一是脈沖相位估計(jì)［9-10］。脈沖相位估計(jì)算法的仿真分析需要以X射線脈沖星光子的到達(dá)時(shí)間作為實(shí)驗(yàn)數(shù)據(jù)［11］。由于觀測(cè)數(shù)據(jù)的精確相位未知，不能進(jìn)行誤差分析;而重構(gòu)的X射線脈沖星信號(hào)其相位已知，因此利用模擬算法重構(gòu)的脈沖星信號(hào)進(jìn)行相位估計(jì)算法仿真十分必要。

X射線脈沖星信號(hào)的模擬有兩種方法，實(shí)物模擬和純數(shù)值模擬。文獻(xiàn)［12］采用實(shí)物模擬方法，該方法成本高、精度低。文獻(xiàn)［13］提出使用可見光源和單光子探測(cè)器搭建實(shí)物模擬系統(tǒng)，提高系統(tǒng)的靈活性和精度。

文獻(xiàn)［14］提出的模擬方法為純數(shù)值模擬，該方法成本低、靈活性強(qiáng)，但模擬過程將非線性泊松過程的強(qiáng)度函數(shù)分段近似為恒定值，極大地降低了模擬的精度。文獻(xiàn)［15-16］都基于反函數(shù)構(gòu)建了數(shù)值模擬算法的數(shù)學(xué)模型，但前者并沒有進(jìn)一步闡述反函數(shù)如何求解，后者則是對(duì)模型進(jìn)行了近似化簡(jiǎn)，避免了反函數(shù)的求解，這導(dǎo)致了信號(hào)模擬的精度不高。

本文提出一種基于精確光子流量函數(shù)的純數(shù)值X射線脈沖星信號(hào)的模擬算法。該算法的數(shù)學(xué)模型基于反函數(shù)建立，脈沖輪廓通過分段線性函數(shù)擬合，從而反函數(shù)易求。然后以χ2擬合優(yōu)度檢驗(yàn)作為驗(yàn)證方法，對(duì)信號(hào)模擬算法的有效性進(jìn)行定量驗(yàn)證。最后比較了三種信號(hào)模擬算法，并將本文提出的算法生成的觀測(cè)脈沖輪廓與標(biāo)準(zhǔn)輪廓比較，對(duì)信號(hào)模擬算法的有效性進(jìn)行驗(yàn)證。

1 信號(hào)模擬算法

1．1 信號(hào)模型

X 射線脈沖星發(fā)射X射線光子是一個(gè)隨機(jī)過程，探測(cè)器記錄的光子到達(dá)時(shí)間則服從非齊次泊松分布。設(shè)光子流量函數(shù)為λ(t)，該函數(shù)表征單位時(shí)間到達(dá)探測(cè)器的光子數(shù)目，標(biāo)準(zhǔn)脈沖輪廓函數(shù)的表達(dá)式與光子流量函數(shù)相同。觀測(cè)時(shí)段為(t0，t)內(nèi)到達(dá)的光子數(shù)目為(t，s)內(nèi)到達(dá)的光子數(shù)目，則Nts服從參數(shù)為的泊松分布，在(t，s)內(nèi)到達(dá) k 個(gè)光子的概率為:

其期望和方差為

其中，φ+φ0為t時(shí)刻對(duì)應(yīng)的相位，φ0為信號(hào)的初始相位和 λs分別是背景輻射密度和脈沖星輻射平均密度，單位與λ相同(ph/s，即每秒接收的光子數(shù)目)，這兩個(gè)量可從天文觀測(cè)數(shù)據(jù)得知。h(φ)以1為周期，且滿足歸一化條件

當(dāng) φ ＞1時(shí)，h(φ)=h(φ -n)，n=floor(φ)，n是不大于φ的最大整數(shù)。

1．2 脈沖輪廓函數(shù)

脈沖輪廓可由歷元折疊方法構(gòu)建。X射線脈沖星的信號(hào)由時(shí)間序列組成，在信號(hào)周期已知的前提下，可以將若干個(gè)周期的信號(hào)折疊到一個(gè)周期內(nèi)。時(shí)間序列可表示為{ti}m1，轉(zhuǎn)換方法如式(5)所示。

其中，T為信號(hào)周期，floor(x)表示不大于x的最大整數(shù)。將一個(gè)周期劃分為Nb個(gè)足夠小的區(qū)間，統(tǒng)計(jì)落在每個(gè)區(qū)間的光子數(shù)目ni，每個(gè)區(qū)間長(zhǎng)度為L(zhǎng)=T/Nb，共有 Np個(gè)周期的數(shù)據(jù)，則第 i個(gè)區(qū)間的平均光子流量函數(shù)

當(dāng)觀測(cè)時(shí)間足夠長(zhǎng)時(shí)Np→∞，Nb→∞每個(gè)區(qū)間的長(zhǎng)度L→0，可以得到光滑、近似連續(xù)的標(biāo)準(zhǔn)脈沖輪廓。

由式(2)可得

其中ti為第i個(gè)時(shí)間段的中點(diǎn)，故

1．3 信號(hào)模擬算法

X射線脈沖星信號(hào)的數(shù)值模擬算法可分為兩種，一種基于光子數(shù)分布率，一種基于分布函數(shù)的反函數(shù)。

1．3．1 基于光子數(shù)分布率的信號(hào)模擬算法

在(t，s)內(nèi)到達(dá)k個(gè)光子的概率P{Nts=k}可由式(1)計(jì)算，將一個(gè)周期劃分為Nb個(gè)足夠小的區(qū)間，每個(gè)區(qū)間稱為一個(gè)時(shí)間倉(cāng)(bin)，第i個(gè)時(shí)間倉(cāng)內(nèi)的光子數(shù)記為Ni，其分布率可記為

k 0 1 … ＞p P{Ni=k} P0 P1 … 1-∑p i=1 Pi

當(dāng)時(shí)間倉(cāng)足夠小時(shí)，每個(gè)時(shí)間倉(cāng)中的光子數(shù)目取值基本為0或1。對(duì)于某一特定長(zhǎng)度的時(shí)間倉(cāng)，大于p個(gè)光子數(shù)的事件為小概率事件，可以將其光子數(shù)近似為p。

根據(jù)分布律可以利用隨機(jī)數(shù)生成每個(gè)時(shí)間倉(cāng)內(nèi)的光子數(shù)目，以時(shí)間倉(cāng)的中點(diǎn)作為該時(shí)間倉(cāng)內(nèi)光子的到達(dá)時(shí)刻。

1．3．2 基于反函數(shù)的信號(hào)模擬算法

X射線脈沖星信號(hào)用光子到達(dá)時(shí)間表示，因此，其模擬算法可按如下方法導(dǎo)出。

設(shè)X射線脈沖星光子的到達(dá)時(shí)間間隔為x，則x的分布函數(shù)與具體的時(shí)刻t有關(guān)，可用條件分布函數(shù)表示為

其反函數(shù)為

由于R是(0，1)區(qū)間上的隨機(jī)變量，故(1-R)也是(0，1)區(qū)間上的隨機(jī)變量，式(12)可簡(jiǎn)化為

又G為光子到達(dá)時(shí)間差，t為當(dāng)前時(shí)刻光子的到達(dá)時(shí)間，故下一個(gè)光子的到達(dá)時(shí)間為

脈沖輪廓λ(φ)的表達(dá)式一般并不規(guī)則，故使用解析方法不易求出Λ(t)和Λ-1。基于反函數(shù)的信號(hào)模擬算法之間的區(qū)別就在于如何處理Λ(t)和 Λ-1。

1)針對(duì)這一問題文獻(xiàn)［16］采用了近似模型，當(dāng)光子到達(dá)時(shí)刻間隔很小，即δt→0時(shí)，

但是當(dāng)脈沖星光子流量較小時(shí)，并不滿足δt→0的條件，不能采用上式進(jìn)行信號(hào)模擬。

2)本文采用分段線性函數(shù)擬合脈沖輪廓函數(shù)λ(φ)。脈沖標(biāo)準(zhǔn)輪廓也是通過歷元折疊構(gòu)建，因此由離散點(diǎn)表示。當(dāng)利用指數(shù)函數(shù)或高斯函數(shù)等擬合時(shí)［17-18］，雖然可以得到平滑的輪廓曲線，但相應(yīng)的求輪廓函數(shù)積分及積分函數(shù)的反函數(shù)時(shí)十分困難。

信號(hào)模擬主要是為了分析相位估計(jì)算法，輪廓形狀會(huì)影響相位估計(jì)的精度，因此必須針對(duì)不同的脈沖星采用其標(biāo)準(zhǔn)脈沖輪廓。但是輪廓是否光滑并不影響結(jié)果。因此采用分段線性函數(shù)擬合脈沖輪廓是可行的，而且λ(φ)分段可積，Λ(t)的反函數(shù)也易求。

設(shè)歷元折疊的時(shí)間倉(cāng)數(shù)目為M，則標(biāo)準(zhǔn)脈沖輪廓的離散點(diǎn)數(shù)目為M+1，為了保證離散點(diǎn)可以表示整個(gè)周期的輪廓，使第一個(gè)點(diǎn)的相位為0，最后一個(gè)點(diǎn)的相位為1。第 i個(gè)離散點(diǎn)表示為，括號(hào)中第一項(xiàng)為相位，第二項(xiàng)為光子流量。線性函數(shù)擬合的脈沖輪廓可表示為

由于式(16)是線性的，所以Λ(t)易求，而Λ(t)是關(guān)于相位φ的二次函數(shù)，故Λ-1可以通過求解二次方程得到，為了簡(jiǎn)化計(jì)算，還可以通過內(nèi)插值得到近似解。值得一提的是，這里的近似解是利用精確模型確定光子到達(dá)時(shí)刻范圍后取近似值，基本不會(huì)影響后續(xù)信號(hào)的模擬，而式(15)則是對(duì)模型取近似，然后用近似的模型生成光子到達(dá)時(shí)刻，該時(shí)刻可能偏離精確模型計(jì)算值較遠(yuǎn)，從而影響后續(xù)信號(hào)的生成。

圖1 crab脈沖星標(biāo)準(zhǔn)脈沖輪廓Fig．1 Standard template pulse profile of crab pulsar

2 實(shí)驗(yàn)與結(jié)果分析

2.1 泊松過程檢驗(yàn)

本節(jié)分析1．3．2節(jié)中提出的基于分段線性函數(shù)擬合的標(biāo)準(zhǔn)脈沖輪廓的信號(hào)模擬算法，采用χ2擬合優(yōu)度檢驗(yàn)法驗(yàn)證生成的信號(hào)是否服從泊松分布。以某一時(shí)間段內(nèi)的X射線光子數(shù)a={0，1，2…ai}，i=1，2，3…作為樣本，記樣本中等于 ai的個(gè)數(shù)為ni并設(shè)a共有k個(gè)不同的取值，由式(18)可求得a=ai=i-1的概率pi，樣本總數(shù)為n=，則皮爾遜χ2統(tǒng)計(jì)量為

對(duì)給定的顯著性水平α，檢驗(yàn)規(guī)則為:

2.2 整段統(tǒng)計(jì)

生成25個(gè)周期，即1s的觀測(cè)數(shù)據(jù)，即t=0，s=1。然后統(tǒng)計(jì)1s內(nèi)觀測(cè)到的光子數(shù)目。皮爾遜χ2統(tǒng)計(jì)量的計(jì)算為:

可得

其中，i=1，2，3，…，取 a={0，1，…，48，［49，∞)}，則 k=50，i的最大值為50，記

n50為大于光子數(shù)，大于等于49的樣本數(shù)。

根據(jù)式(17)即可計(jì)算出χ2。

由信號(hào)模擬算法生成5000個(gè)樣本數(shù)據(jù)，ni的分布如圖2所示，由圖可見，在光子數(shù)集中在期望值附近，而大于或小于期望達(dá)到某一值后，樣本數(shù)量變?yōu)?0。將 ni代入式(17)求得 χ2=24．263 8，又．338 6 故)，說明在一整段時(shí)間內(nèi)產(chǎn)生的光子到達(dá)時(shí)刻服從泊松分布。

圖2 光子數(shù)目分布圖Fig．2 Distribution of the number of photons

2.3 分段統(tǒng)計(jì)

生成250個(gè)周期，即10s的觀測(cè)數(shù)據(jù)，然后進(jìn)行歷元折疊后，將一個(gè)周期均分為8個(gè)部分，對(duì)每部分的光子數(shù)進(jìn)行χ2擬合優(yōu)度檢驗(yàn)，皮爾遜χ2統(tǒng)計(jì)量計(jì)算方法與2．1節(jié)相同。各段的光子個(gè)數(shù)期望E(Nts)值計(jì)算如表1所示?？紤]到pi過小

表1 分段統(tǒng)計(jì)各段的參數(shù)及計(jì)算結(jié)果Tab．1 Parameters and results of each section for the section-divided method

可能會(huì)導(dǎo)致極大的計(jì)算誤差，因此將p

i

＜10

-3

的光子數(shù)目按區(qū)間統(tǒng)計(jì)，a的取值集合為(［0，a

1

］，

，實(shí)驗(yàn)數(shù)據(jù)及結(jié)果如表 1 所示。由表中數(shù)據(jù)可知每一段的統(tǒng)計(jì)量χ

2

均小于

，說明在每一段上產(chǎn)生的光子數(shù)目也服從泊松分布。

2.4 信號(hào)模擬算法比較

前面介紹了兩種信號(hào)模擬算法，并提出了一種基于分段線性函數(shù)擬合脈沖輪廓函數(shù)的模擬算法，將基于光子數(shù)分布率的信號(hào)模擬算法記為方法一，將采用式(15)由反函數(shù)推導(dǎo)的信號(hào)模擬算法記為方法二，基于分段線性函數(shù)擬合脈沖輪廓函數(shù)的模擬算法也利用反函數(shù)導(dǎo)出，與方法二的區(qū)別在于計(jì)算過程不采用近似，將其記為方法三。下面分別比較三種算法生成的光子數(shù)期望和由模擬信號(hào)構(gòu)建的觀測(cè)脈沖輪廓。

2.5 光子數(shù)目比較

根據(jù)X射線脈沖星的光子流量分布規(guī)律，取三組不同的λs和λb，

1)λs=15ph/s，λb=5ph/s

2)λs=150ph/s，λb=50ph/s

3)λs=1500ph/s，λb=500ph/s

表2 光子數(shù)目比較Tab．2 Comparison of photon counts

觀測(cè)時(shí)間取100s，統(tǒng)計(jì)1000次實(shí)驗(yàn)的光子數(shù)均值，結(jié)果如表2所示，由表可見方法二的光子數(shù)目與理論值相差較大，特別是當(dāng)光子流量較小時(shí)尤為明顯，而方法一和方法三則趨近于理論值。

2．6 觀測(cè)脈沖輪廓比較

由于當(dāng)光子流量較小時(shí)，方法二不能有效模擬脈沖信號(hào)。故取 λs=1500ph/s，λb=500ph/s，觀測(cè)時(shí)間取100s，歷元折疊時(shí)間倉(cāng)取128，構(gòu)建觀測(cè)脈沖輪廓。進(jìn)行1000次實(shí)驗(yàn)，統(tǒng)計(jì)觀測(cè)脈沖輪廓每個(gè)離散點(diǎn)取值的均方差，并與理論值比較，結(jié)果如圖3所示。方法二模擬的脈沖信號(hào)其輪廓誤差最大，其次為方法一，而本文提出的方法三誤差趨近于理論值，說明該方法可以重構(gòu)標(biāo)準(zhǔn)脈沖輪廓。

圖3 脈沖輪廓誤差與理論值的比較Fig．3 Comparison between the error of pulse profile and analysis result

圖4 顯示的是將模擬產(chǎn)生的 10s，100s，1000s的光子到達(dá)時(shí)刻通過歷元折疊后得到的累積脈沖輪廓，可發(fā)現(xiàn)當(dāng)時(shí)間越長(zhǎng)時(shí)，脈沖輪廓越接近標(biāo)準(zhǔn)脈沖輪廓函數(shù)。進(jìn)一步證明了方法三可以重構(gòu)的脈沖星信號(hào)。

3 結(jié)論

本文提出一種基于精確光子流量函數(shù)的純數(shù)值重構(gòu)X射線脈沖星信號(hào)的模擬算法，該算法依據(jù)X射線脈沖星信號(hào)的光子到達(dá)數(shù)目服從泊松分布性質(zhì)，用概率論的相關(guān)方法模擬光子到達(dá)的隨機(jī)事件，并采用分段線性函數(shù)擬合脈沖輪廓解決了輪廓函數(shù)積分及其反函數(shù)不易求得的問題。該方法不需要借助其他硬件設(shè)備即可生成模擬的X射線脈沖星信號(hào)，成本低、靈活性高，可以作為X射線脈沖星導(dǎo)航實(shí)驗(yàn)的信號(hào)源。通過擬合優(yōu)度檢驗(yàn)驗(yàn)證，本算法重構(gòu)的脈沖星信號(hào)服從泊松分布，符合脈沖星信號(hào)特征。

圖4 累積脈沖輪廓與標(biāo)準(zhǔn)脈沖輪廓的比較Fig．4 Comparison of accumulated pulse profile and standard pulse profile

將提出的新算法與兩種已有的算法比較，結(jié)果顯示新算法在光子數(shù)和脈沖輪廓誤差上都最接近理論值，說明其最接近實(shí)際信號(hào)。

通過歷元折疊檢驗(yàn)新算法表明其有效可行。將三種時(shí)間長(zhǎng)度的模擬信號(hào)進(jìn)行歷元折疊，獲得其累積脈沖輪廓，顯示當(dāng)累積時(shí)間足夠長(zhǎng)時(shí)，脈沖輪廓趨近于光子流量函數(shù)(即標(biāo)準(zhǔn)脈沖輪廓)。表明該算法可有效模擬X射線脈沖星信號(hào)，為X射線脈沖星導(dǎo)航的研究提供了信號(hào)源基礎(chǔ)。

References)

［1］ Hosken R W，Wertz JR．Microcosm autonomous navigation system on-orbitoperation［C］//Proceedings of the Annual AAS Guidance and Control Conference，Keystone，Colorado，1995:491-506．

［2］ 楊博．航天器星敏感器自主定位方法及精度分析［J］．宇航學(xué)報(bào)，2002，23(3):81-84．YANG Bo．The method of autonomous position determination using star sensor and analysis of precision for spacecraft［J］．Journal of Astronautics，2002，23(3):81-84．

［3］ Counley R，White R，Gai E．Autonomous satellite navigation by stellar refraction［J］．Journal of Guidance，Control，and Dynamics，1984，7(2):129-134．

［4］ Wertz JR．Spacecraft attitude determination and control［M］．Kluwer Academic Publishers，1978．

［5］ Downs G，Reichley P E．Techniques formeasuring arrival times of pulsar signals I:DSN observations from 1968 to 1980［J］．NASA Jet Propulsion Laboratory，California Institute of Technology，Pasadena CA，NASA Technical Reports NASACR-163564，1980．

［6］ 楊廷高，仲崇霞．脈沖星時(shí)穩(wěn)定度及可能應(yīng)用［J］．時(shí)間頻率學(xué)報(bào)，2004，27(2):129-137．YANG Tinggao，ZHONG Chongxia．Frequency stability of pulsar time scale and possible application［J］．Journal of time and frequency，2004，27(2):129-137．

［7］ 桂先洲，黎勝亮，李志豪．基于X射線脈沖星絕對(duì)定位中的整周模糊度改進(jìn)算法研究［J］．國(guó)防科技大學(xué)學(xué)報(bào)，2010，32(3):33-36．Gui Xianzhou，Li Shengliang，Li Zhihao．Study of betterment arithmetic of ambiguity in absolution position based on X-ray pulsar［J］． Journal of National University of Defense Technology，2010，32(3):33-36．

［8］ Sheikh S I，Hanson J E．Spacecraft navigation and timing using X-ray pulsar［J］．Navigation，2011，58(2):165-186．

［9］ 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．

［10］ Rinauro S，Colonnese S，Scarano G．Fast near-maximum likelihood phase estimation of X-ray pulsars［J］．Signal Processing，2013，93(1):326-331．

［11］ Sheikh S I．The use of variable celestial X-ray sources for spacecraft navigation［D］．University of Maryland，2005．

［12］ 胡慧君，趙寶升，盛立志，等．基于X射線脈沖星導(dǎo)航的地面模擬系統(tǒng)研究［J］．物理學(xué)報(bào)，2011，60(2):853-861．Hu Huijun，Zhao Baosheng，Sheng Lizhi，etal．A simulation experiment system for X-ray pulsar based navigation［J］．Acta Physica Sinica，2011，60(2):853-861．

［13］ 孫海峰，謝楷，李小平，等．高穩(wěn)定度X射線脈沖星信號(hào)模擬［J］．物理學(xué)報(bào)，2013，62(10):109701．Sun Haifeng，Xie Kai，Li Xiaoping，et al．A simulation technique of X-ray pulsar signals with high timing stability［J］．Acta Physica Sinica，2013，62(10):109701．

［14］ 蘇哲，許錄平，王婷．X射線脈沖星導(dǎo)航半物理仿真實(shí)驗(yàn)系統(tǒng)研究［J］．物理學(xué)報(bào)，2011，60(11):819-826．Su Zhe， Xu Luping， Wang Ting． X-ray pulsar-based navigation semi-physical simulation experiment system［J］．Acta Physica Sinica，2011，60(11):819-826．

［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］ 黎勝亮，劉昆，吳錦杰．脈沖輪廓建模仿真與去噪研究［J］．國(guó)防科技大學(xué)學(xué)報(bào)，2013，35(4):26-29．Li Shengliang，Liu Kun，Wu Jinjie． On modeling and denoising of X-ray pulsar profile［J］．Journal of National University of Defense Technology，2013，35(4):26-29．

［17］ 徐軒彬，吳鑫基．脈沖星PSR2111+46平均脈沖分析和譜特性［J］．中國(guó)科學(xué) A 輯，2002，32(12):1134-1141．Xu Xuanbing， Wu Xing Ji． Spectrum characteristic and analysis ofmean pulse of pulsar PSR2111+46［J］．Science China(A)，2002，32(12):1134-1141．

［18］ Rankin JM．Toward an empirical theory of pulsar emission．IV-Geometry of the core emission region［J］．The Astrophysical Journal，1990，352(3):247-257．