基于Bayesian估計的X射線脈沖星微弱信號檢測*

張 華 許錄平 謝 強 羅 楠

(西安電子科技大學電子工程學院，西安 710071)(2010年4月8日收到;2010年7月10日收到修改稿)

基于Bayesian估計的X射線脈沖星微弱信號檢測*

張 華 許錄平 謝 強 羅 楠

(西安電子科技大學電子工程學院，西安 710071)(2010年4月8日收到;2010年7月10日收到修改稿)

累積輪廓、流量和周期是X射線脈沖星輻射信號的三個重要特征，將其應用于X射線脈沖星信號檢測中，提出了一種基于Bayesian估計的X射線脈沖星周期輻射信號時域檢測方法.該方法以非脈沖區噪聲觀測為先驗知識，利用X射線脈沖星輻射信號的泊松分布模型推導了信號概率密度分布函數，以該函數的累積分布函數為判據，對X射線脈沖星微弱信號進行檢測，并提取位相偏移量.利用仿真數據和RXTE衛星的實測數據進行實驗驗證，結果表明:本文方法性能優于同類的基于高斯分布模型的檢測方法，在檢測信號的同時能在一定精度下給出信號位相偏移值.

脈沖星，Bayesian估計，位相測量，時域檢測

PACS:97.60.Gb，98.70.Qy，96.60.tk

1.引 言

脈沖星是一類高速穩定自旋的中子星，其輻射信號覆蓋了射電、紅外、可見光、紫外、X射線及γ射線波段［1］.其中X射線波段集中了輻射信號的大部分能量，易于小型化設備探測和處理，使基于X射線脈沖星的航天器定時、定姿和定位成為可能［2］.脈沖星距太陽系通常都有幾千光年，空間X射線探測器探測到的X射線脈沖星輻射信號非常微弱且湮沒在背景噪聲中，因此X射線脈沖星信號檢測就成為應用中的一大難題.常用的檢測方法以基于FFT的頻域周期檢測法為主［3］，有時還需要將脈沖星信號進行周期累加得到脈沖星累積脈沖輪廓模型，再提取平移不變特征量進行檢測［4，5］.但基于FFT的頻域檢測方法僅在高斯白噪聲前提下適用，而且當脈沖星的輻射脈沖較窄時，信號能量分散于諧波中，若噪聲包含頻率成分與諧波接近，則諧波能量得到增強，這將導致信號被誤檢測，比如，利用該方法對ROSAT望遠鏡觀測的光變曲線數據進行周期檢測就存在失效情況［6］;此外，這種頻域方法計算量大，這就限制了其在計算能力有限的星載計算機中的實現.相對于這種頻域方法而言，時域方法可以避免頻域變換時能量分散于諧波中，因此對信息利用更為充分，可以得到更優的檢測性能［6］.文獻［7］最早提出了利用統計推理方法對感興趣區域脈沖信號進行檢測的時域方法，但其前提是信號和噪聲呈高斯分布.本文對這一時域方法進行拓展，根據X射線脈沖星輻射信號的泊松分布模型，提出一種新的基于 Bayesian估計［8］的 X射線脈沖星輻射信號時域檢測方法，討論了其在位相測量中的作用，并通過實驗說明了該方法的有效性.

2.基于Bayesian估計的X射線脈沖星輻射信號檢測

2.1.脈沖星輻射信號劃分及其建模

脈沖星的脈沖周期和輪廓特征在X射線探測器探測的原始信號上表現為X射線光子流的平均強度變化，為了同時利用觀測數據所包含的脈沖星輻射信號的輪廓信息、強度信息和周期信息，根據累積輪廓的結構，將脈沖星信號分為信號段和噪聲段.以噪聲段的噪聲觀測作為先驗知識，利用Bayesian估計推導出信號段信號強度的后驗概率分布.不失一般性，以脈沖星B1937+21為例，如圖1，該脈沖星一段觀測信號累積出的脈沖輪廓包含一個主脈沖和一個次脈沖，組成信號段，其在一個脈沖周期內持續的時間分別設為 WP1，WP2，令 WP=WP1+WP2;其他部分組成噪聲段，在一個脈沖周期內持續的時間分別設為 WN1，WN2和 WN3，令 WN=WN1+WN2+WN3;WP+WN=T.

圖1 脈沖星信號劃分

空間X射線光子到達探測器具有隨機性，若將單位時間內到達探測器的光子數視為隨機事件，可用泊松分布對該事件進行建模［9，10］.由于 X射線脈沖星輻射信號強度周期性變化，光子到達探測器的密度隨時間分布是不均勻的，定義落在區間(t1，t2)內的光子數等于k的概率為［9］

式中，λ(t)表示 t時刻光子流量密度，λ(t)=λs(t)+λn，其中λs(t)為信號輻射強度，當脈沖周期為T時，有λs(t)=λs(t+T)，λn為各種噪聲強度的總和，包括背景輻射噪聲、系統熱噪聲等，簡單起見視噪聲平均強度為常數;t2－t1=ΔT為觀測時間.

2.2.X射線脈沖星信號后驗概率表示

為了簡化分析，假設:1)背景輻射噪聲平穩，信號段的背景噪聲分布與噪聲段一致，且可用噪聲段觀測數據建模;2)一次分析時間為一個脈沖周期，在整個觀測期間假設脈沖周期不變;3)在背景輻射噪聲強度未知時假設其服從均勻分布.定義:Cg，Cb和Cp分別為信號段光子數、噪聲段光子數和一個脈沖周期時間觀測的光子總數，則Cp=Cb+Cg，Cg=定義:μp，μb和μg分別為信號光子平均計數率、背景光子平均計數率和信號段總的光子計數率，有μg=μp+μb.對X射線脈沖星信號分析的目的即在已知Cb和Cg的前提下，估計后驗概率 P(μp/Cb，Cg).由查普曼-科爾莫格洛夫方程有

在平均背景輻射強度已知時，即μb已知時，對于泊松分布模型而言，可認為 P(μp/μb)=P(μg= μp+μb)，得到

其中 μg僅與 Cg有關，有

由X射線脈沖星輻射的泊松分布模型可得

根 據 Bayesian 原 理 有 P(μg/Cg) =，假設 P(Cg)服從均勻分布，有E［P(μg)］WP=E［P(Cg)］，近似得到

在假設1)前提下，用噪聲段的觀測估計μb，

與(6)式同理得到

將(3)，(4)，(6)，(7)和(8)式代入(2)式得到

其中μg=μp+μb，根據二項式定理有(μgWP)Cg=，代入(9)式得到

再根據 Gamma積分［11］

對(10)式積分運算得到

(12)式中 P(μp/Cb，Cg)表示在已知觀測量Cb，Cg時該觀測中存在信號的可能性，當觀測到的信號段總計數Cg趨近無窮大時，P(μp/Cb，Cg)趨近于1.

2.3.基于后驗概率累積分布函數的X射線脈沖星信號檢測算法

通過對(12)式的后驗概率密度函數積分得到其累積分布函數轉4);

4)若 j≤k，s(t)=s(t+jΔt)轉2)，否則轉5);

5)求 R(n)(0≤n≤k)的方差δ，搜索是否存在R(n)>3δ，若存在說明檢測到信號.

2.4.基于信號檢出率的位相偏移測量算法

在圖1所示的脈沖星信號劃分中，噪聲段 WN

從2.1節和2.2節的分析知道，Cb，Cg分別反映了噪聲和含噪聲的信號的先驗知識，(13)式中 PCb，Cg則反映了從先驗知識Cb，Cg中獲取的信號的信息.先驗知識越準確，Cg中包含的信號信息越多，PCb，Cg越大，當PCb，Cg大于或等于某一置信水平時，就認為探測到了信號.對于X射線探測器接收到的信號而言，認為其是以采樣時間為單位的光子計數率的時間序列，用 s(nΔt)表示該時間序列，s(n1Δt，n2Δt)表示序列中的某一段，T表示脈沖星周期，Δt表示采樣間隔，T=kΔt，k 為每周期采樣數，以 PCb，Cg為判據，檢測序列中是否包含特定脈沖星信號的步驟為

1)對某顆脈沖星，按圖1所示的方法將其輪廓劃分為信號段WP和噪聲段WN，設閾值水平為L，令i=0，j=0，HS=0，HB=0，待檢測序列長度為 mT(m為整周期數);

2)對 s(ikΔt，(i+1)kΔt)計算 PCb，Cg，若 PCb，Cg≥ L，HS=HS+1，否則 HB=HB+1;

3)若 i

式中s∈{0，k－1}表示起始位相，P表示觀測整周期數.具體實現方法為在2.3節的步驟5)中提取r=max

n{R(n)}.由于時域采樣的離散性，時域位相提取的精度受最小采樣間隔的限制.Taylor FFT等頻域方法［5］可以擺脫這種限制，使位相測量精度僅與累積輪廓信噪比有關，但頻域方法中諧波分量會導致多解情況，造成真實解判定的混亂.因此，位相時域解最少可以作為頻域解搜索過程中的初值，為更精確的位相提取提供有力的參考.

3.脈沖星特征參數對檢測性能影響的理論分析

為了分析各參數對新方法檢測性能影響，以(13)式所表示的概率分布函數 PCb，Cg為分析對象，分析結果如圖2所示:Cg，Cb和 β分別表示每脈沖周期信號段光子計數、噪聲段光子計數和累積輪廓占空比;圖 2(a)以 Cg=0，1，…，5 ph(光子)和 Cb=3 ph為例，說明 PCb，Cg隨 β 變化情況;圖 2(b)以Cb=0，1，…，5 ph 和 Cg=3ph 為例，說明 PCb，Cg隨 β變化情況;圖 2(c)以 β =0.5，Cg=0，1，…，5 ph 為例，說明 PCb，Cg隨 Cb變化情況;圖 2(d)以 β =0.5，Cb=0，1，…，5 ph 為例，說明 PCb，Cg隨 Cg變化情況.

圖2 概率分布函數與各參數的關系 (a)Cb=3 ph，PCb，Cg與 β 關系;(b)Cg=3ph，PCb，Cg與 β 關系;(c)PCb，Cg與 Cb關系;(d)PCb，Cg與 Cg關系

利用累積分布 PCb，Cg作為判據，對于圖2(a)而言，β 一定時，PCb，Cg在 Cg=0，1，…，5 ph 時的各曲線之間的差值越大，說明檢測靈敏度越高，檢測性能越好;通過觀察該差值隨β分布可以發現，其僅在β∈［0.2，0.8］時較大，說明本文方法不適用于占空比過大或過小的脈沖星信號檢測，而且該差值在β∈［0.2，0.5］時的值大于 β∈［0.5，0.8］的值，說明本文方法對信號段較窄的脈沖星檢測性能更好些;對于圖 2(b)而言，β 一定時 PCb，Cg在 Cb=0，1，…，5 ph的各曲線之間的差值越大，說明檢測靈敏度受噪聲影響越大，與圖2(a)的分析方法相同，可以發現本文方法檢測信號段較窄的脈沖星信號時對噪聲的抑制能力強.從圖2(a)和圖2(b)中還發現，當 β一定時，Cg越大檢測性能越好;相反，Cb越大檢測性能越差，這與一般物理規律是符合的.圖2(c)和圖2(d)進一步說明了這一結論.如圖2(c)所示，隨著 Cb增加，PCb，Cg減小，Cg=0，1，…，5 ph 的 PCb，Cg各曲線間的差也在減小，說明檢測靈敏度降低;反之，如圖2(d)所示，對Cb一定的信號而言，隨著Cg增強，PCb，Cg增大，則檢測性能得到增強，當 Cb越小時，檢測性能隨Cg增強越快.

4.實驗與分析

4.1.仿真實驗及分析

以輻射強度較弱的 B1937+21和 B1821－24兩顆脈沖星為例，主要參數如表1.設X射線探測器有效探測面積為1 m2，平均背景輻射強度為5×102ph/cm2.s－1，根 據 脈 沖 星 周 期、標 準 累 積 輪廓［12—14］、平均流量和每周期采樣次數計算每兩次采樣之間信號流量δ(t)和背景輻射流量α，單周期仿真信號表示為:nΔt))+ α}，其中 S為采樣點數，poissrnd(·)為泊松隨機數生成函數.設探測置信水平為95%(虛警率水平5%)，做1000次Monte Carlo實驗，分別統計兩顆脈沖星信號檢出率和虛警率與累積周期數關系，結果如圖3(a).可以看到，隨著累積周期數的增加，檢測性能逐漸提高，且對 B1821-24檢測性能提高更快.需要說明的是，由于周期累加過程中，噪聲也增強了，根據圖2(c)的結果，為增加靈敏度，累加后的信號應消除部分累加噪聲，這里減去最小值并乘以系數0.02.為了說明輪廓、流量和周期參數差異對檢出率的影響，將B1821-24的各個參數分別修改為與B1937+21一致，再進行仿真，結果如圖3(b).從圖3(b)分析，當脈沖星形狀、信號流量和周期改變時，本文方法在檢測性能上存在差異，對于脈沖星B1821－24而言，輪廓占空比影響比較大，但這僅是個例，不同的脈沖星受各參數的影響在不同條件下是不同的，圖2的分析可以說明這一點.

為了評價本文方法檢測性能，與同類的基于高斯分布的檢測方法對比.基于高斯分布的檢測方法假設信號和背景噪聲服從高斯分布，其噪聲抽取方法為［7］

圖3 檢出率與累積周期的關系 (a)檢出率比較;(b)參數差異影響

此時判決閾值定義為

其中，δ0為信號方差，K為系數，這里令 K=2，當μgauss≤Lc，表示檢測到信號，否則表示出現虛警.設累積周期為10000，以B1937+21為例，信號檢出率與虛警率的差和背景輻射強度的關系如圖4所示.當背景噪聲較小時，基于高斯分布檢測方法的性能與本文方法的性能相當，但隨著背景計數率的增加，基于高斯分布檢測方法性能迅速下降，本文方法性能優于基于高斯分布的檢測方法.

圖4 與高斯檢測方法性能的比較

對于觀測時間內可能出現多個信號的情況，由于新方法引入了周期和輪廓特征信息，對多個信號間的干擾具有一定的免疫能力.如圖6(a)所示，設信號和背景噪聲輻射強度均為 5×10－2ph/cm2·s－1，在僅包含 B1821-24信號的一段觀測中滑動檢測B1937-21信號，然后疊加 B1937-21信號采用同樣的方法檢測，包含B1937-21信號可以觀察到明顯的尖峰，說明檢測到信號.同理，如圖6(b)所示，在B1937-21信號中疊加 B1821-24信號，得到類似的檢測結果.

圖5 位相測量的均方差

表1 X射線脈沖星參數:周期、X射線流量密度、每周期采樣次數、占空比

本文方法利用了X射線脈沖星的輪廓信息，可以提取出觀測信號的位相，不過所提取位相的精度受到時域采樣間隔的限制，精度上不如頻域方法［5］，但頻域方法受諧波分量的影響，不可避免出現多解情況，因此時域位相解至少可以作為頻域方法初值，提高搜索真實位相解的速度.對于新方法在位相測量的性能，本文利用B1821-24和B1937+21的仿真數據做了仿真分析.設背景噪聲輻射強度均為 5 ×10－2ph/cm2·s－1，仿真中先隨機移動仿真數據的初始位相，然后利用本文方法滑動檢測初始位相，記錄下真實位相與檢測位相的差，多次測量求均方差，該均方差與累積周期數的關系如圖5所示.可見隨著疊加次數的增加，整體上位相測量方差在減小，在信噪比相同時，B1937+21位相測量均方差較小，這與B1937+21采樣率高的原因有直接關系.

圖6 混合多脈沖星信號情況下的檢測 (a)檢測B1937-21;(b)檢測B1821-24

4.2.基于RXTE實測數據的實驗及分析

從HEASARC數據庫的Public Data Access連接進入，以B1855+09為關鍵字搜索出以該脈沖星為觀測對象的觀測號，然后從其FTP［15］上下載與該觀測號對應的數據，利用FTOOLS工具處理得到該脈沖星X射線光變曲線.經過修正，剔除掉不良數據，得到256個周期長的信號，對該信號經累加去均值歸一化后的輪廓如圖7所示.由于該段數據背景計數率較大，為了應用本文方法，先逐周期疊加，再減去最小值，乘以系數 0.02，求 PCb，Cg，然后將信號循環右移1rad位相，用同樣的方法再計算PCb，Cg，結果如圖8所示，可見:位相無誤時，隨著累加周期數的增加，累積分布PCb，Cg逐漸提高，而當位相偏移1rad時，PCb，Cg隨著累加周期數的變化并不明顯，說明了本方法在檢測信號和辨別位相時的有效性.為了進一步說明本文方法在位相測量上作用，將原信號位相右移πrad，然后從零位相開始滑動檢測，PCb，Cg隨滑動位相變化情況如圖9所示，可見位相對 PCb，Cg有直接影響，準確的位相能最大化 PCb，Cg，從而通過PCb，Cg可以提取出位相.

圖7 歸一化累積輪廓

圖8 分布與累積周期數的關系

圖9 PCb，Cg隨位相變化情況

5.結 論

基于Bayesian估計的X射線脈沖星信號檢測方法利用X射線脈沖星信號輪廓特征、周期特征和強度等信息作為信號檢測的基本參數，在時域對X射線脈沖星信號檢測，避免了頻域方法中諧波分量引起的能量分散導致的誤檢測，在檢測性能上，優于同類的基于高斯分布模型的檢測方法.利用RXTE的實測數據也說明了該方法的有效性.在位相測量方面，該方法能在一定精度要求下給出位相偏移值，雖然受時域采樣精度的限制，性能上不如Taylor FFT等頻域方法，但至少能作為頻域方法的初值，減少頻域方法多解情況引起的誤判，提高頻域方法位相搜索的速度.

［1］ Pei Y J，Jun Z 2008 Chin.Phys.B 17 1

［2］ Sheikh S I 2005 Ph.D.Dissertation(Maryland:University of Maryland)

［3］ Seward F D，Harnden F R，Helfand D J 1984 Astrophysical Journal 287 L19

［4］ Xie Z H，Xu L P，Ni G R 2007 J.Infrared Millim.Waves 26 187(in Chinese)［謝振華、許錄平、倪廣仁 2008紅外與毫米波學報 26 187］

［5］ Xie Z H，Xu L P，Ni G R 2007 Acta Phys.Sin.26 187(in Chinese)［謝振華、許錄平、倪廣仁 2007物理學報 26 187］

［6］ Gregory P C 1996 Astrophysical Journal 10 1

［7］ Currie L A 1968 Analytical Chemistry 40 586

［8］ Gregory P 2003 Bayesian Logical Data Analysis for the Physical Sciences(U.K.:Cambridge Univ)p5

［9］ Sala J，Urruela A，Villares X 2004 European Space Agency Advanced Concepts Team ARIADNA Study June 4202 23

［10］ Li J X，Ke X Z 2010 Acta Phys.Sin.59 8304(in Chinese)［李建勛、柯熙政 2010物理學報 59 8304］

［11］ Young J M K A 2009 IEEE Trans.Nucl.Sci.56 1278

［12］ Epn.http://www.jb.man.ac.uk/research/pulsar/resources/epn/browser.html.

［13］ Kuiper L， Hermsen W 2003 proceedings of AGILE Science Workshop Rome，Italy，Dec，2003 p11

［14］ Nicastro L，Cusumano G，Lohmer O，Kramer M，Kuiper L，Hermsen W，Mineo T，Becker W 2004 Astron.Astrophys 413 1065

［15］ Rxte.FTP:legacy.gsfc.nasa.gov.

PACS:97.60.Gb，98.70.Qy，96.60.tk

X-ray pulsar weak signal detection based on Bayesian estimation*

Zhang Hua Xu Lu-Ping Xie Qiang Luo Nan
(School of Electronic Engineering，Xidian University Xi’an 710071，China)(Received 8 April 2010;revised manuscript received 10 July 2010)

Integrated profile，flux and period are three major characteristics of pulsar.Taking advantage of these characteristics，a time domain detecting method is presented based on Bayesian estimation to detect pulsar periodic radiation signals.The signal probability distribution is deduced based on the Poisson distribution model of X-ray pulsar，with noise observation of non-pulse region used as the priori knowledge.The cumulative distribution function of the signal probability is used as a criterion to detect weak signals of X-ray pulsars and extract the phase offset.The RXTE data and the simulation data are used for experimental verification，and the results show that the method outperforms a similar method which is based on the Gauss distribution model，in addition，it can give the phase offset in a certain accuracy.

pulsars，Bayesian estimation，phase measurement，time domain detection

*國家高技術研究發展計劃(批準號:2007AA12Z323)，國家自然科學基金(批準號:60772139)資助的課題.

E-mail:zhanghua@mail.xidian.edu.cn

*Project supported by the National High Technology Research and Development Program of China(Grant No.2007AA12Z323)，the National Natural Science Foundation of China(Grant No.60772139).

E-mail:zhanghua@mail.xidian.edu.cn