蘇 哲,王 勇,許錄平,羅 楠

(西安電子科技大學電子工程學院,西安710071)

0 引 言

X射線脈沖星導航(XPNAV:X-ray Pulsar-based Navigation)是一種新的天文自主導航技術,可為近地軌道、深空和星際空間飛行的航天器提供精確的位置、速度、姿態(tài)和時間等導航信息。利用觀測信號辨識脈沖星是XPNAV的重要環(huán)節(jié)[1]。根據辨識結果從星載導航數據庫中提取該脈沖星的輻射方向,可推出航天器的姿態(tài);已知航天器姿態(tài)定向觀測某顆脈沖星時,也需要根據辨識結果來驗證姿態(tài)的正確性和準確度[2];根據脈沖星類型從導航數據庫中提取其發(fā)射頻率,與實測脈沖頻率比較,測定多普勒頻移量,可計算航天器在脈沖星視線上的運動速率;根據脈沖星類型提取導航數據庫中的標準輪廓和星歷等參數,測量脈沖到達時間,可解算出航天器在三維空間中的位置[3]。因此,在有限的觀測時間內,快速、準確的辨識脈沖星對于航天器定姿、測速和定位都具有重要意義。迄今為止發(fā)現的脈沖星中,沒有哪兩顆具有完全相同的形成過程和輻射機理,所以可通過周期和累積脈沖輪廓唯一確定一顆脈沖星。但是,在軌航天器觀測的脈沖星信號十分微弱,且受到多普勒效應和噪聲的影響,測量所得周期值存在誤差[4],同時,許多適于XPNAV的毫秒脈沖星的周期非常接近,故采用周期辨識難以取得良好的辨識效果。累積脈沖輪廓被認為是脈沖星極冠區(qū)的輻射窗口,這種由物理性質所決定的信號特征極其穩(wěn)定。因此,XPNAV系統(tǒng)中可通過比較累積脈沖輪廓和導航數據庫中的標準輪廓辨識脈沖星編號。

通過累積脈沖輪廓辨識脈沖星的基礎是特征向量的提取。由于疊加起始時刻不同,累積脈沖輪廓和標準輪廓間存在相位差,這要求特征向量具有平移不變性。同時,在XPNAV系統(tǒng)中,由于多普勒效應、相位間隔的不同和脈沖周期測量的不準確性,累積脈沖輪廓相對于標準輪廓還存在尺度伸縮。綜上所述,在XPNAV系統(tǒng)中通過累積脈沖輪廓辨識脈沖星,提取累積脈沖輪廓中具有平移和尺度伸縮不變性的特征信息,并抑制噪聲,是解決問題的關鍵。

現有的射電脈沖星累積脈沖輪廓辨識算法主要有一維選擇線譜算法[5]。該算法在選擇雙譜辨識算法[6]的基礎上,改進了特征向量的降維方法。這兩種辨識算法均利用雙譜的平移不變性提取特征信息,并抑制高斯噪聲,辨識效果相差不大。由于雙譜不具有尺度伸縮不變性,這兩種算法均無法有效辨識具有尺度伸縮效應的累積脈沖輪廓,不適于XPNAV系統(tǒng)。

針對XPNAV系統(tǒng)中累積脈沖輪廓的辨識問題,本文在雙譜的基礎上,引入Mellin變換,提出了一種基于選擇Bispectra-Mellin(BM)譜的累積脈沖輪廓辨識算法。該算法通過BM變換提取累積脈沖輪廓中的平移和尺度伸縮不變特征信息,并濾除噪聲,然后利用簡化的Fisher可分離度對BM幅度譜降維,提取特征向量,用于相關辨識。該算法能有效抑制累積脈沖輪廓的平移、尺度伸縮和噪聲對辨識的影響,運算量小,適用于XPNAV系統(tǒng)。

1 XPNAV系統(tǒng)中累積脈沖輪廓的特點分析

RXTE(Rossi X-ray Timing Explorer)是美國NASA探索者系列衛(wèi)星中具有高時間分辨率(1μs),大有效面積(6250 cm2)和寬能譜范圍(2-200keV)的X射線探測器,具有對突發(fā)現象迅速反應并實施探測的能力,適于對X射線脈沖星進行觀測研究。本文從美國高能天文數據中心(HEASARC:High Energy Astrophysics Science Archive Research Center)獲取RXTE的脈沖星觀測數據進行研究。數據的提取和初步處理采用HEASARC提供的專用工具包FTOOLS。

以脈沖星B0531+21為例,濾除RXTE的無效觀測數據后,進行時間轉換,周期疊加,得到該脈沖星的標準輪廓(累積時間略大于7小時,相位間隔設為2.8125度),如圖1所示。其脈沖寬度比(半高寬度與周期的比值)為43.12%。

圖1 脈沖星B0531+21的標準輪廓,相位間隔為2.8125度Fig.1 Thestandard profile of pulsar B0531+21,the phasebin is 2.8125 degree

模擬XPNAV系統(tǒng)中的情況,假設未知觀測信號的脈沖星類型,利用FTOOLS工具包中的Powspec工具估計其自轉頻率,并設自轉頻率的1、2階導數均為零。截取一段10秒種的觀測數據,進行光子到達時間的轉換后,將相位間隔設為2.8125度,進行周期疊加。由于RXTE的軌道數據更新時間為60秒,所以時間轉換無法消除10秒內衛(wèi)星運動引起的多普勒效應。所得累積脈沖輪廓如圖2(a)所示。計算其脈沖寬度比為43.58%。利用標準輪廓計算其信噪比為9.43dB。對比圖1和2(a)可看出,航天器高速運動引起的多普勒效應和對信號疊加時采用的周期不準確,導致累積脈沖輪廓略微展寬,相對于標準輪廓的尺度伸縮幅度為0.989。

圖2 脈沖星B0531+21累積脈沖輪廓,累積時間10秒,相位間隔分別為2.8125度(a)和7.2度(b)Fig.2 Theintegrated pulse profile of pulsar B0531+21,theintegrated time is 10 seconds,and thephase bin is 2.8125 degree(a)and 7.2 degree(b)respectively

將相位間隔設為7.2度,其他條件不變,進行周期疊加,所得累積脈沖輪廓如圖2(b)所示。其信噪比為11.21dB,相對于標準輪廓的尺度伸縮幅度為2.53。對比圖2(a)和2(b),在相同的累積時間下,2(b)具有較長的相位間隔,單個相位間隔內的累積時間較長,其信噪比高于2(a)。因此,適度延長相位間隔雖然會導致累積脈沖輪廓產生尺度伸縮,但卻能大幅度提高其信噪比。換句話說,如果能找到一種可屏蔽尺度伸縮的方法,就能夠極大地縮短辨識所需累積脈沖輪廓的疊加時間,從而縮短整個XPNAV系統(tǒng)導航參數的輸出時間間隔。

Taylor所建立的射電脈沖星累積脈沖輪廓模型[7]中,將累積脈沖輪廓看作是標準輪廓時延和加噪后的結果。但是,圖2所示累積脈沖輪廓除受到上述因素影響外,也產生了尺度伸縮,伸縮幅度分別為0.98和2.53。為準確描述累積脈沖輪廓相對于標準輪廓的尺度伸縮,引入尺度因子,將其表達式改為:

上式中,a為直流偏差,b為幅度因子,f(t)為標準輪廓,g(t)為噪聲,服從Poisson分布,d為時間延遲,s為尺度因子。

由(1)式可知,通過比較累積脈沖輪廓p(t)和導航數據庫中的標準輪廓f(t)來辨識脈沖星,辨識算法必須能夠屏蔽時間延遲(d)和尺度(s)上的差異,提取具有平移和伸縮不變性的特征信息,并抑制噪聲g(t)。

2 基于選擇BM譜的脈沖星累積脈沖輪廓辨識算法

2.1 Bispectra-Mellin變換及其平移和尺度伸縮不變性

實信號x(t)同其伸縮量x(st)的Mellin變換為:

可見,信號x(t)尺度的變化僅僅導致其Mellin變換的相位發(fā)生變化,而對幅值沒有影響。本文對兩維信號x(t1,t2)的Mellin變換做出如下定義:

尺度伸縮后的信號x(s1 t1,s2 t2)的Mellin變換為:

可見,(4)式所定義的兩維Mellin變換的幅值同樣具有伸縮不變性,可用于提取信號中的尺度伸縮不變特征量。

借鑒一維Mellin變換的實現方法——DMT算法[8],本文給出一種兩維Mellin變換的實現方法。假設x(t1,t2)在很小的積分區(qū)域{(t1,t1+Δt),(t2,t2+Δt)}內為常量,則可將(4)式展開為:

為提取平移和尺度伸縮不變特征量,并抑制噪聲的影響,本文將兩維Mellin變換和雙譜相結合,給出實信號x(t)的Bispectra-Mellin變換的定義:

其中,max(·)為取極大值運算符,用于對雙譜幅度進行歸一化,消除幅度因子 b對辨識的影響。B(ω1,ω2)是信號 x(t)的雙譜。

文獻[9]指出,當每個相位間隔中的光子數均大于20時,Poisson噪聲g(t)近似服從零均值Gaussian分布。理論上,Gaussian過程的雙譜為零。因此,假設p(t)中直流分量 a已經濾除,且累積時間足夠長使得每個相位間隔中光子數均大于20,易知信號f(t)和p(t)的雙譜間有如下關系:

結合(7)式,可得f(t)和p(t)的BM譜間關系為:

可知,BM幅度譜具有平移和尺度伸縮不變性,且繼承了雙譜抑制噪聲的能力,適用于在XPNAV系統(tǒng)中提取脈沖星累積脈沖輪廓的特征信息。

2.2 BM幅度譜的降維

由于BM幅度譜維數較高,不利于辨識,需對其進行將維處理。降維的過程本質上是從高維的BM幅度譜中選擇最能代表不同類型脈沖星間差異的點,組成特征向量?，F有辨識算法中,降維的方法主要有Fisher可分離度方法[6]和一維選擇線譜方法[5]。在實驗中發(fā)現,一維選擇線譜方法的運算量極大,星載計算機運算能力有限,難以利用該算法實時地辨識脈沖星信號。因此,本文仍采用Fisher可分離度方法對BM幅度譜降維,并針對脈沖星的特點對Fisher可分離度進行了適當簡化,使之更適于XPNAV系統(tǒng)。

脈沖星標準輪廓是脈沖星極冠區(qū)的輻射窗口,在軟X射線頻段標準輪廓具有唯一性。因此,在提取特征向量時,每類脈沖星信號僅需一個訓練樣本——由長時間觀測數據疊加所得標準輪廓。此時,僅需尋找能使不同類別間差異最大的若干點組成特征向量,即可取得良好的辨識效果。

假定 {M(i)(ρ1,ρ2)}i=1,2,…,I是訓練樣本集合的BM變換幅度譜,其中上標 i表示不同編號的脈沖星,I表示訓練樣本集合中脈沖星的種類數目。該訓練樣本集合在(ρ1,ρ2)點的可分離度定義為:

根據上式,計算BM譜域內所有點的 m(ρ1,ρ2),選擇具有最大 m(ρ1,ρ2)值的 N個點組成特征向量,即可有效辨識該訓練樣本內的脈沖星。

實驗中發(fā)現,N的取值非常重要,N太小,特征向量不足以涵蓋該訓練樣本集合中所有脈沖星標準輪廓的特征信息,導致辨識效果不佳;N太大,特征向量包含了不同類標準輪廓間的共有信息,不僅不能提高辨識效果,還會增加運算負擔。

為確定N的取值,定義一種衡量辨識效果的性能指標——可辨識度。累積脈沖輪廓的可辨識度r定義為:

其中,corc為脈沖星累積脈沖輪廓和其標準輪廓特征向量的相關系數,corm為累積脈沖輪廓和訓練樣本中其它類型脈沖星標準輪廓的特征向量相關系數的最大值。易知,r＞1時,根據該累積脈沖輪廓可得出正確的辨識結果,且 r越大,辨識效果越好,抗干擾能力越強;r≤1時,無法得到正確的辨識結果。

實驗中,對于選定的訓練樣本,首先令 N=2,計算脈沖星的可辨識度;增加N的取值,可辨識度也不斷提高;當N大于N0時,繼續(xù)增加N的值時,可辨識度不再有明顯提高,有時還會略微減小。此時,對于該訓練樣本,可取N=N 0。

3 實驗與結果分析

實驗數據由RXTE觀測得到,取自HEASARC數據庫,數據格式為FITS(Flexible Image Transport System)。實驗工作在Matlab7.04環(huán)境下完成。

3.1 累積脈沖輪廓的BM變換

分別計算圖1所示標準輪廓和圖2(b)所示累積脈沖輪廓的雙譜和BM幅度譜,如圖3、4所示。在計算雙譜變換時,為抑制噪聲,在雙譜域加入Rao-Gabr平滑窗[10],窗口長度設為5倍于采樣周期。在計算兩維Mellin變換時,根據雙譜的對稱性,僅取雙譜域右上角1/4區(qū)域的點參與計算,可在降低運算量的同時不丟失特征信息。

對比圖3(a)和4(a)可以看出,標準輪廓和累積脈沖輪廓(尺度因子2.53,信噪比11.21)的雙譜具有相似性,但尺度不同,說明雙譜變換不具有伸縮不變性。對比圖3(b)和4(b)可以看出,BM幅度譜圖形相似,且矯正了雙譜尺度上的差異,說明BM變換具有平移和尺度伸縮不變性。

圖3 脈沖星B0531+21標準輪廓的雙譜幅度譜(a)和BM幅度譜(b)的等高線圖Fig.3 The contour line of bispectra(a)and BM spectra(b)amplitude of the standard profile of pulsar B0531+21

圖4 脈沖星B0531+21累積脈沖輪廓的雙譜幅度譜(a)和BM幅度譜(b)的等高線圖Fig.4 The contour line of bispectra(a)and BM spectra(b)amplitude of the integrated pulse profile of pulsar B0531+21

3.2 辨識實驗

選擇雙譜算法和一維選擇線譜算法均是利用雙譜的平移不變性提取特征信息。通過圖3(a)和4(a)可知,無論采用什么樣的降維方法,均很難通過雙譜域上的特征向量辨識不同尺度的累積脈沖輪廓。下面將進一步通過比較選擇雙譜算法和本文算法的辨識效果說明該問題。

從文獻[11]列出的適于導航的脈沖星(具有較高品質因子)中選擇12顆已被RXTE觀測過的脈沖星 :B1929+10、B1055-52、B1046-58、B1323-62、B0531+21 、B1713-40 、B1719-37 、B1742-30 、B1749-28、B1821-24、J1730-2304 和 J1617-5055,按照先后順序,其編號分別為1～12。將相位間隔均設為2.8125度,長時間累積得到它們的標準輪廓,作為訓練樣本集合。

分別計算訓練樣本集合中各標準輪廓的BM幅度譜。然后,利用(10)式計算BM譜域內每個點的可分離度。對于該訓練樣本,取N=10。提取具有最大可分離度的10個點作為標準輪廓的特征向量。待辨識信號仍選用B0531+21的累積脈沖輪廓,累積時間10秒,相位間隔7.2度,相對于標準輪廓的尺度因子2.53,每個相位間隔中的光子數均大于300,噪聲服從高斯分布。計算其BM幅度譜,并從中提取特征向量。分別計算該特征向量與訓練樣本集合中所有標準輪廓特征向量的相關系數,結果如圖5中實線所示。其中,橫坐標表示脈沖星編號。圖5中虛線為通過選擇雙譜算法提取脈沖星標準輪廓和累積脈沖輪廓的特征向量,進行相關運算的結果。

圖5 脈沖星B0531+21累積脈沖輪廓的特征向量與訓練樣本集合中標準輪廓特征向量的相關系數Fig.5 The feature vectors correlation coefficients of theintegrated pulseprofileof B0531+21 and standard profiles in training set

可以看出,當尺度因子 s≠1時,BM幅度譜所提取的累積脈沖輪廓的特征向量和其標準輪廓的特征向量具有較高的相關性,可用于辨識脈沖星類型;反觀選擇雙譜算法,由于雙譜不具有尺度伸縮不變性,相同類型脈沖星的累積脈沖輪廓和標準輪廓的特征向量的相關性較低,難以用于辨識。

3.3 相位間隔長度對辨識效果的影響

為討論相位間隔長度與辨識效果的關系,分別計算具有不同相位間隔的累積脈沖輪廓的可辨識度。訓練樣本集合仍采用前一節(jié)中列舉的12顆脈沖星的標準輪廓。對B0531+21的同一段觀測數據(10秒鐘)進行周期疊加,設置不同的相位間隔,可得不同尺度因子和信噪比的累積脈沖輪廓,分別計算它們的可辨識度,如表1所示。

可以看出,選擇雙譜算法僅當尺度因子近似為1時,可辨識度略大于1;而本文算法對于不同的尺度因子,均可有效辨識。而且隨尺度的增加,可辨識度總體上有提高的趨勢。這與增加相位間隔長度使累積脈沖輪廓的信噪比增加相吻合。但是,當相位間隔長度過大時(大于等于9),由于累積脈沖輪廓無法提供足夠多的可供辨識的信息,可辨識度小于1。當相位間隔為4.8度時,可辨識度最高,此時累積脈沖輪廓不僅能提供足夠多的辨識信息而且信噪比較高。

表1 不同相位間隔的累積脈沖輪廓的可辨識度Table 1 The recognition rateof the accumulated pulse profileswith different phase bin sizes

3.4 運算量分析

下面從運算量的角度考察本文算法。實驗用計算機基本配置為Intel雙核2.2GHz CPU,2GDDR2內存。表2對比了選擇雙譜算法和本文算法所耗時間。

表2 選擇雙譜和本文算法耗時對比Table2 The comparison of consuming time between selected bispectra algorithm and the proposed algorithm

在標準輪廓特征向量的提取階段,需分別計算訓練樣本在BM譜域上每個點的可分離度,然后從中選擇具有最大可分離度的若干點作為特征向量。該階段運算量較大,但為前期準備階段,并不影響算法的實時性,因此可選擇經長時間累積、具有高信噪比、一個周期內有較多相位間隔的標準輪廓作為訓練樣本集合。對累積脈沖輪廓的辨識階段,僅需計算累積脈沖輪廓BM變換的若干點,組成特征向量,并依次與已經存儲于計算機內存中的標準輪廓特征向量進行相關匹配,運算量較小。本文算法與選擇雙譜算法相比,需另外計算若干點的兩維Mellin變換,故其辨識階段的運算量略大于選擇雙譜算法,但毫秒級的辨識時間同脈沖星信號的累積時間(本實驗中為10 s)相比,幾乎可忽略不計。

4 結 論

在有限的觀測時間內,快速、準確的辨識脈沖星類型對于利用X射線脈沖星定姿、測速和定位具有重要意義。本文研究了XPNAV系統(tǒng)中累積脈沖輪廓的特點,針對XPNAV系統(tǒng)提出了一種基于選擇BM譜的脈沖星累積脈沖輪廓辨識算法。理論分析和實驗結果表明:(1)BM幅度譜具有平移和尺度伸縮不變性,且能抑制噪聲;(2)本文算法可有效辨識具有平移和尺度伸縮效應的累積脈沖輪廓,辨識效果優(yōu)于現有辨識算法;(3)利用相同的觀測數據,通過合理選擇相位間隔的長度(4.8度),本文算法能夠達到最佳可辨識度;(4)本文算法運算量小,適用于XPNAV系統(tǒng)。

[1] SHEIKH S I.The use of variable celestial X-ray sources for spacecraft navigation[D].Maryland:Maryland University,2005.

[2] SHEIKE SI,PINES D J.Spacecraft navigation using X-ray pulsars[J].Journal of Guidancecontrol and dynamics,2006,29(1):49-63.

[3] 帥平,陳紹龍,吳一帆,等.X射線脈沖星導航原理[J].宇航學報,2007,28(6):1538-1543.[SHUAI Ping,CHENG Shaolong,WU Yi-fan,et al.Navigation principles using X-ray pulsars[J].Journal of Astronautics,2007,28(6):1538-1543.]

[4] 蘇哲,許錄平,王光耀,等.基于離散方波變換的脈沖星微弱信號周期性檢測[J].宇航學報,2009,30(6):2243-2248.[SUZhe,XULu-ping,WANGGuang-yao,et al.Pulsar weak signal periodicity detection based on discrete square wave transform[J].Journal of Astronautics,2009,30(6):2243-2248.]

[5] XIEZheng-hua,XU Lu-ping,NIGuang-ren,et al.A new feature vector using selected line spectra for pulsar signal bispectrum characteristic analysis and recognition[J].Chinese Journal of Astronomy and Astrophysics,2007,7(4):565-571.

[6] ZHANG Xian-da,YUShi,BAOZheng.A new feature vector using selected bispectra for signal classification with application in radar target recognition[J].IEEE Trans.on Signal Processing,2001,49(9):1875-1885.

[7] Taylor JH.Pulsar timing and relativistic gravity[J].Class.Quantum Grav,1993,10:S167-S174.

[8] Zwicke P E,Kiss I.A new implementation of the mellin transform and its application to radar classification of ships[J].IEEE Trans.on Pattern Analysis and Machine Intelligence,1983,5(2):191-198.

[9] Hanson J,Sheikh S,Graven P,et al.Noise analysisfor X-ray navigation systems[C].2008 IEEE/ION Position,Location and Navigation Symposium,Monterey,USA:IEEE,2008:704-713.

[10] Nikias C L,Petropulu A P.Higher-order spectra analysis[M].New Jersey:PTRPrentice-Hall,1993.

[11] Josep Sala,Andreu Urruela,Xavier Villares.Feasibility study for a spacecraft navigation system relying on pulsar timing information final report[R].Technical University of Catalonia Final Report for ARIADNA project 03/4202,2004.