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

桂先洲，黃森林，孫 晨

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

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

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

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

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

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

1 信號模擬算法

1．1 信號模型

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

其期望和方差為

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

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

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

脈沖輪廓可由歷元折疊方法構建。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ù)

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

由式(2)可得

其中ti為第i個時間段的中點，故

1．3 信號模擬算法

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

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

在(t，s)內(nèi)到達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

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

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

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

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

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

其反函數(shù)為

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

又G為光子到達時間差，t為當前時刻光子的到達時間，故下一個光子的到達時間為

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

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

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

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

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

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

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

圖1 crab脈沖星標準脈沖輪廓Fig．1 Standard template pulse profile of crab pulsar

2 實驗與結果分析

2.1 泊松過程檢驗

本節(jié)分析1．3．2節(jié)中提出的基于分段線性函數(shù)擬合的標準脈沖輪廓的信號模擬算法，采用χ2擬合優(yōu)度檢驗法驗證生成的信號是否服從泊松分布。以某一時間段內(nèi)的X射線光子數(shù)a={0，1，2…ai}，i=1，2，3…作為樣本，記樣本中等于 ai的個數(shù)為ni并設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ù)集中在期望值附近，而大于或小于期望達到某一值后，樣本數(shù)量變?yōu)?0。將 ni代入式(17)求得 χ2=24．263 8，又．338 6 故)，說明在一整段時間內(nèi)產(chǎn)生的光子到達時刻服從泊松分布。

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

2.3 分段統(tǒng)計

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

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

可能會導致極大的計算誤差，因此將p

i

＜10

-3

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

1

］，

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

2

均小于

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

2.4 信號模擬算法比較

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

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ù)均值，結果如表2所示，由表可見方法二的光子數(shù)目與理論值相差較大，特別是當光子流量較小時尤為明顯，而方法一和方法三則趨近于理論值。

2．6 觀測脈沖輪廓比較

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

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

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

3 結論

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

圖4 累積脈沖輪廓與標準脈沖輪廓的比較Fig．4 Comparison of accumulated pulse profile and standard pulse profile

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

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

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］．宇航學報，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)定度及可能應用［J］．時間頻率學報，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］．國防科技大學學報，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射線脈沖星導航的地面模擬系統(tǒng)研究［J］．物理學報，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］．物理學報，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射線脈沖星導航半物理仿真實驗系統(tǒng)研究［J］．物理學報，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］．國防科技大學學報，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］．中國科學 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．