重構(gòu)X射線脈沖星信號的純數(shù)值模擬新算法*

桂先洲，黃森林，孫 晨

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

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

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

X射線脈沖星信號的模擬有兩種方法，實物模擬和純數(shù)值模擬。文獻(xiàn)［12］采用實物模擬方法，該方法成本高、精度低。文獻(xiàn)［13］提出使用可見光源和單光子探測器搭建實物模擬系統(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ù)如何求解，后者則是對模型進(jìn)行了近似化簡，避免了反函數(shù)的求解，這導(dǎo)致了信號模擬的精度不高。

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

1 信號模擬算法

1．1 信號模型

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

其期望和方差為

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

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

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

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

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

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

由式(2)可得

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

1．3 信號模擬算法

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

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

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

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

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

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

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

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

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

其反函數(shù)為

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

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

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

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

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

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

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

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

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

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

2 實驗與結(jié)果分析

2.1 泊松過程檢驗

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

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

2.2 整段統(tǒng)計

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

可得

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

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

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

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

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

2.3 分段統(tǒng)計

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

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

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

i

＜10

-3

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

1

］，

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

2

均小于

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

2.4 信號模擬算法比較

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

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

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

2．6 觀測脈沖輪廓比較

由于當(dāng)光子流量較小時，方法二不能有效模擬脈沖信號。故取 λs=1500ph/s，λb=500ph/s，觀測時間取100s，歷元折疊時間倉取128，構(gòu)建觀測脈沖輪廓。進(jìn)行1000次實驗，統(tǒng)計觀測脈沖輪廓每個離散點(diǎn)取值的均方差，并與理論值比較，結(jié)果如圖3所示。方法二模擬的脈沖信號其輪廓誤差最大，其次為方法一，而本文提出的方法三誤差趨近于理論值，說明該方法可以重構(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á)時刻通過歷元折疊后得到的累積脈沖輪廓，可發(fā)現(xiàn)當(dāng)時間越長時，脈沖輪廓越接近標(biāo)準(zhǔn)脈沖輪廓函數(shù)。進(jìn)一步證明了方法三可以重構(gòu)的脈沖星信號。

3 結(jié)論

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

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

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

通過歷元折疊檢驗新算法表明其有效可行。將三種時間長度的模擬信號進(jìn)行歷元折疊，獲得其累積脈沖輪廓，顯示當(dāng)累積時間足夠長時，脈沖輪廓趨近于光子流量函數(shù)(即標(biāo)準(zhǔn)脈沖輪廓)。表明該算法可有效模擬X射線脈沖星信號，為X射線脈沖星導(dǎ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é)報，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］ 楊廷高，仲崇霞．脈沖星時穩(wěn)定度及可能應(yīng)用［J］．時間頻率學(xué)報，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射線脈沖星絕對定位中的整周模糊度改進(jìn)算法研究［J］．國防科技大學(xué)學(xué)報，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é)報，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射線脈沖星信號模擬［J］．物理學(xué)報，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)航半物理仿真實驗系統(tǒng)研究［J］．物理學(xué)報，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］．國防科技大學(xué)學(xué)報，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］．中國科學(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．