張 彬，梁 楠,2，于徐紅，陳宗浩，劉志杰，游善平*

(1. 貴州師范大學(xué)貴州省信息與計算科學(xué)重點實驗室，貴州 貴陽 550001；2. 中國科學(xué)院云南天文臺，云南 昆明 650216)

脈沖星是一種快速旋轉(zhuǎn)、高磁化的中子星，是理解超密度物質(zhì)的關(guān)鍵，可以用于各種基礎(chǔ)物理實驗[1]。目前，脈沖星搜索方法主要分為周期性搜索和單脈沖搜索兩大類。其中，周期性搜索通過應(yīng)用快速傅里葉變換(Fast Fourier Transform, FFT)將時間序列轉(zhuǎn)化到頻域，以識別周期性信號，然后在確定的周期內(nèi)對原始時間序列數(shù)據(jù)進行折疊，以提高周期性信號的信噪比[2]。而單脈沖搜索尋找強的、非周期的脈沖，并未使用快速傅里葉變換和折疊。單脈沖搜索非常適合發(fā)現(xiàn)周期性搜索中無法發(fā)現(xiàn)的孤立爆發(fā)。研究人員應(yīng)用單脈沖搜索發(fā)現(xiàn)了旋轉(zhuǎn)射電暫現(xiàn)源(Rotating Radio Transients, RRAT)[3]和快速射電暴的(Fast Radio Bursts, FRB)[4]。

2003年，文[5]首次提出了一個理論框架，運用單脈沖搜索方法探測射電天文觀測數(shù)據(jù)中的脈沖星信號。2009年，文[6]把單脈沖搜索應(yīng)用于阿雷西博L波段饋源陣列脈沖星巡天(Pulsar Arecibo L-band Feed Array, PALFA)七波束觀測數(shù)據(jù)，發(fā)現(xiàn)了7顆脈沖星。2010年，文[7]使用單脈沖搜索方法對帕克斯多波束脈沖星巡天(Parkes Multibeam Pulsar Survey, PMPS)數(shù)據(jù)進行重新分析，發(fā)現(xiàn)10顆旋轉(zhuǎn)射電暫現(xiàn)源。2015年，文[8]在單脈沖理論的基礎(chǔ)上，設(shè)計針對旋轉(zhuǎn)射電暫現(xiàn)源的搜索工具RRATtrap，應(yīng)用在綠岸射電望遠鏡(Green Bank Telescope, GBT)350 MHz觀測數(shù)據(jù)，最終探測到18顆旋轉(zhuǎn)射電暫現(xiàn)源。2016年，文[9]首次將單脈沖搜索與機器學(xué)習(xí)相結(jié)合，提出遞歸峰值識別算法(Recursive Algorithm for Peak IDentification, RAPID)，實現(xiàn)了自動化篩選脈沖星候選體，并探測到6顆脈沖星。

2016年9月,位于貴州喀斯特山區(qū)的500 m口徑球面射電望遠鏡進入科學(xué)運行階段，望遠鏡先后配置了超寬帶接收機和L波段19波束接收機[10]。為充分利用FAST優(yōu)異的探測能力，2018年，文[11]設(shè)計并實施了FAST多科學(xué)目標(biāo)同時掃描巡天，即同時使用多個數(shù)字終端采集脈沖星、中子氫、分子譜線、旋轉(zhuǎn)射電暫現(xiàn)源、快速射電暴等多個科學(xué)目標(biāo)的觀測數(shù)據(jù)。鑒于單脈沖搜索方法在各大射電望遠鏡的成功應(yīng)用，本文對CRAFTS超寬帶數(shù)據(jù)文件的單脈沖搜索結(jié)果進行研究，發(fā)現(xiàn)CRAFTS單脈沖搜索篩選得到脈沖星候選體存在數(shù)以萬計的假陽性樣本。因此，我們根據(jù)脈沖星信號在色散量信噪比曲線大致是高斯曲線以及信噪比數(shù)據(jù)分布是高斯分布的假設(shè)，構(gòu)造3個顯著區(qū)分脈沖星信號和干擾的特征，旨在緩解出現(xiàn)大量假陽性樣本的問題。試驗結(jié)果表明，添加特征值判斷的單脈沖篩選工具與傳統(tǒng)方法相比，假陽性樣本數(shù)量減少了20%。

1 CRAFTS超寬帶數(shù)據(jù)單脈沖搜索實驗

在射電天文觀測數(shù)據(jù)中發(fā)現(xiàn)脈沖星通常分為4個階段：收集、消色散、周期性搜索或單脈沖搜索和人工檢查[1]。第1階段，原始數(shù)據(jù)以電壓時間序列的形式由射電望遠鏡收集；第2階段，脈沖輻射通過星際介質(zhì)(Interstellar Medium, ISM)的色散效應(yīng)，導(dǎo)致較低頻率的脈沖比較高頻率晚到達[12]，消色散可以去除這些與頻率有關(guān)的延遲效應(yīng)的影響；第3階段，使用周期性搜索或單脈沖搜索找出脈沖星候選體；第4階段，對判斷為脈沖星候選體的數(shù)據(jù)進行人工檢查。

我們先使用PRESTO(Pulsar Exploration and Search Toolkit)脈沖星搜索工具處理2017年8月至2018年5月CRAFTS的超寬帶數(shù)據(jù)，處理流程包括去干擾、消色散和單脈沖搜索等步驟；然后，應(yīng)用單脈沖篩選工具RRATtrap從PRESTO的Single_pulse_search.py輸出文件中篩選出脈沖星候選體。

單脈沖篩選工具RRATtrap根據(jù)脈沖星信號的特性區(qū)分脈沖星信號和干擾。(1)脈沖星信號出現(xiàn)在一定色散量(Dispersion Measure, DM)范圍內(nèi)，在最佳色散量時檢測到峰值信噪比，而在該色散量上下信噪比下降，這是由于脈沖以不準(zhǔn)確的色散量消色散時導(dǎo)致脈沖展寬造成的。(2)由于信號在最佳色散量檢測到峰值信噪比，預(yù)計來自地面的信號(即射頻干擾)在色散量為0 pc·cm-3時達到峰值。同時射頻干擾(Radio Frequency Interference, RFI)不受色散效應(yīng)的影響，因此，脈沖星信號出現(xiàn)在一個非常大的色散量范圍內(nèi)[8]。

圖1是單脈沖篩選工具RRATtrap探測到已知脈沖星PSR B0540+23的診斷圖。左側(cè)子圖縱坐標(biāo)為試驗色散量范圍，橫坐標(biāo)為觀測時間，CRAFTS超寬帶數(shù)據(jù)每個觀測文件為52 s，右側(cè)子圖描述試驗色散量對應(yīng)的信噪比大小。從圖1可以看到呈紡錘狀的單脈沖事件組(在相鄰的色散量與時間窗口內(nèi)所有信噪比大于5的事件)在最佳色散量78 pc·cm-3時取得峰值信噪比68，而在該色散量上下，信噪比下降。其中，強的脈沖星信號標(biāo)記為Excellent；弱的脈沖星信號表標(biāo)記為Very good；不規(guī)則的脈沖星信號標(biāo)記為Good或Ok；射頻干擾標(biāo)記為RFI；宇宙噪聲標(biāo)記為Noise。在圖2中，我們看到在較寬的色散量范圍上極強的射頻干擾的實例。此外，在t=20～40 s，極強的射頻干擾錯誤標(biāo)記為脈沖星信號。

圖1 RRATtrap探測到脈沖星PSR B0540+23的診斷圖

圖2 RRATtrap篩選出假陽性實例的診斷圖

本文試驗選取2017年8月至2018年5月共計約32萬個超寬帶CRAFTS數(shù)據(jù)文件[1]進行單脈沖搜索試驗。由于處理的數(shù)據(jù)規(guī)模較大，對所有數(shù)據(jù)進行全面的人工檢查工作量非常大。因此，我們先使用RRATtrap初步篩選脈沖星候選體。然后，我們對數(shù)據(jù)文件中標(biāo)記為Excellent，Very good，Good和Ok共計約5萬顆脈沖星候選體的診斷圖進行手工檢查，發(fā)現(xiàn)僅有772個診斷圖真正具有脈沖星信號，對應(yīng)101顆已知脈沖星。通過

(1)

計算的假陽性率(False Positive Rate, FPR)達到98.5%，其中，F(xiàn)P是沒有脈沖星信號的數(shù)據(jù)文件標(biāo)記為候選脈沖星的數(shù)量；TP是含有脈沖星信號的數(shù)據(jù)文件標(biāo)記為候選脈沖星的文件數(shù)量。每100個脈沖星候選體中，僅有1～2個包含脈沖星信號(已知或新發(fā)現(xiàn)的)。由于CRAFTS超寬帶數(shù)據(jù)文件中包含的具體脈沖星數(shù)量未知，本文未對假陰性情況進行分析。

2 區(qū)分脈沖星信號和干擾的特征

CRAFTS超寬帶數(shù)據(jù)單脈沖搜索結(jié)果表明，脈沖星候選體僅有1.5%真正具有脈沖星信號，這主要是未能成功區(qū)分脈沖星信號和干擾造成的，如圖2。本節(jié)首先探討不準(zhǔn)確的色散量對信噪比的影響。然后，構(gòu)造3個顯著區(qū)分脈沖星信號和干擾的特征，并在探測到的101顆脈沖星樣本中選取同時具有脈沖星信號、射頻干擾和宇宙噪聲的79顆脈沖星作為樣本，計算它們在3個特征的取值情況。最后，根據(jù)脈沖星信號、射頻干擾和宇宙噪聲在3個特征取值分布的差異，提出合理閾值應(yīng)用在單脈沖篩選工具，對脈沖星候選體進行進一步篩選。

2.1 色散量與信噪比曲線

文[5]探討了試驗與真實色散量之間的偏差對信噪比的影響。通常，對一個特定色散量的時間序列進行多次下采樣并重新搜索，當(dāng)有效的采樣時間最接近脈沖寬度時，得到的信噪比最高。在不同的色散量信道中，隨著試驗色散量與真實色散量的偏差越大，信噪比越小，如圖1。測量信噪比S(δDM)與真實信噪比S的比值與色散量的偏差δDM滿足[5]

(2)

這里，

(3)

其中，erf為誤差函數(shù)；δDM為試驗色散量與真實色散量的偏差；Δv為總帶寬；v為中心頻率；W為脈沖寬度，單位ms。

利用(2)式和(3)式計算天體物理脈沖的信噪比和寬度，我們可以計算預(yù)期的色散量偏差δDM。圖3展示了在中心頻率546 MHz，時間分辨率為100 μs和匹配濾波使用30的PSR B2000+40一組單脈沖事件預(yù)期和擬合信噪比下降曲線。在這兩種情況下，信噪比峰值處的色散量是真實的色散量。我們計算預(yù)期的信噪比下降，利用峰值信噪比和匹配濾波得到脈沖寬度，和實際脈沖寬度很接近。首先，我們使用非線性最小二乘法(Non-linear Least Squares, NLS)對原始數(shù)據(jù)進行回歸，得到擬合的峰值信噪比和脈沖寬度。然后，根據(jù)擬合結(jié)果，代入(2)式得到擬合的色散與信噪比曲線。由圖3可以看出，預(yù)期和擬合的色散量與信噪比曲線大致是高斯曲線。使用觀測到的峰值信噪比和匹配濾波得到脈沖寬度，預(yù)期的信噪比下降與擬合值在峰值右側(cè)非常接近，但在左側(cè)出現(xiàn)偏差。這是匹配濾波得到的脈沖寬度與實際脈沖寬度的偏差造成的。

2.2 3個顯著特征

文[13]指出，脈沖星的輪廓可以簡化為高斯型(對大多數(shù)脈沖星來說是一個合理的近似)。按照文[5]提出的單脈沖搜索理論，天體物理脈沖的色散量和信噪比曲線大致是高斯曲線，如圖3。由于高斯曲線通常是對稱的，文[14]提出了兩種對稱特征(SIDM和SIS/N)分別表征單脈沖事件組的色散量和信噪比的對稱性。SIDM的計算公式為

(4)

其中，DMpeak為信噪比峰值對應(yīng)的色散量；DMmax為單脈沖事件組最大的色散量；DMmin為單脈沖事件組最小的色散量。SIS/N的計算公式為

(5)

其中，∑S/Nleft為信噪比峰值左側(cè)所有單脈沖事件的信噪比之和；∑S/Nright為信噪比峰值右側(cè)所有單脈沖事件的信噪比之和。

由(4)式和(5)式定義的SIDM和SIS/N的取值范圍在0～1之間，這兩個特征值越高，色散量與信噪比曲線越對稱，相反不遵循(2)式描述規(guī)律的射頻干擾，通常在色散量與信噪比空間是單調(diào)遞減(或遞增)。所以，它們通常具有接近于0的對稱值，可以與脈沖星信號進行區(qū)分。圖4和圖5分別展示了脈沖星樣本SIDM和SIS/N的對稱值，并對脈沖星信號、射頻干擾和噪聲進行了對比。從圖4和圖5可以看出，脈沖星信號通常具有較大的對稱值，而射頻干擾和噪聲的對稱值相對較小，這表明脈沖星信號的色散量和信噪比曲線比干擾更對稱。此外，圖4和圖5標(biāo)記一些對稱值比較大的噪聲，這些通常是高斯噪聲。

圖4 脈沖星樣本SIDM對稱值

圖5 脈沖星樣本SIS/N對稱值

文[5]預(yù)測在沒有任何宇宙噪聲和射頻干擾的情況下，信噪比的數(shù)據(jù)分布是高斯分布。峰度是描述總體數(shù)據(jù)分布與高斯分布陡緩程度的統(tǒng)計量。為了判斷脈沖星樣本信噪比的數(shù)據(jù)分布是否符合高斯分布，我們使用峰度值表征單脈沖事件組內(nèi)信噪比的數(shù)據(jù)分布與高斯分布的符合程度。峰度值的計算公式為

(6)

其中，S/N為單脈沖事件的信噪比；μ為單脈沖事件組內(nèi)平均信噪比；σ為單脈沖事件組內(nèi)信噪比的方差。

當(dāng)(6)式定義的峰度值是3時，信噪比的數(shù)據(jù)分布服從高斯分布，隨著峰度值與3的差距越大，其分布形態(tài)的陡緩程度與高斯分布的差異程度越大。圖6用箱線圖展示了脈沖星樣本的峰度值分布情況，并對脈沖星信號、射頻干擾和噪聲進行了對比。從圖6可以看出，脈沖星信號峰度值集中在3附近，而射頻干擾和噪聲的峰度值分散，且與3差距較大，表明脈沖星信號信噪比分布近似服從高斯分布，而大部分干擾與噪聲不具有此規(guī)律。

圖4～圖6可以明顯看出脈沖星信號、射頻干擾和宇宙噪聲在3個特征分布的差異。因此，我們根據(jù)脈沖星信號和干擾特征值的分布區(qū)間(表1)，選擇合適的特征閾值，剔除不滿足閾值的脈沖星候">選體，從而達到降低假陽性率的目的(比如設(shè)置脈沖星信號特征SIDM最小值0.1作為閾值，對低于閾值的候選體視為干擾，那么，特征值在0.02～0.1區(qū)間的射頻干擾將不會錯誤標(biāo)記為脈沖星候選體)。同時，我們選取不同閾值進行試驗，并統(tǒng)計其對應(yīng)的假陽性率和遺漏脈沖星信號的結(jié)果。由圖7各個特征在不同閾值假陽性及遺漏脈沖星的情況可以看出，隨著閾值變大或者閾值區(qū)間縮小，假陽性率不斷減小，與此同時，遺漏發(fā)現(xiàn)的脈沖星數(shù)量在不斷增加。脈沖星搜尋的首要前提是保證觀測數(shù)據(jù)中所有脈沖星信號不遺漏。因此，我們選取表1中脈沖星信號特征SIDM最小值0.1和SIS/N最小值0.29作為閾值，特征峰度-1.54～27.11作為閾值區(qū)間，對低于閾值或不在閾值區(qū)間的脈沖星候選體視為射頻干擾或宇宙噪聲。

圖6 脈沖星樣本峰度值

表1 3個特征在3種類別的分布區(qū)間

圖7 各個特征在不同閾值假陽性和遺漏脈沖星的情況，其中紅色豎線代表所選擇的閾值

此外，為了證實所述特征是否在其他觀測數(shù)據(jù)有效，我們對帕克斯多波束數(shù)據(jù)進行了試驗。根據(jù)文[15]公開的帕克斯單脈沖數(shù)據(jù)庫，我們構(gòu)建了一個PMPS(Parkes Multibeam Pulsar Survey)數(shù)據(jù)庫，其中包含帕克斯望遠鏡發(fā)現(xiàn)的部分脈沖星信號。我們計算得到PMPS數(shù)據(jù)集中脈沖星信號在3個特征值的分布情況，并與CRAFTS進行對比，如表1。由表1中CRAFTS及PMPS的脈沖星信號在3個特征值分布情況可以發(fā)現(xiàn)，它們的特征值分布十分接近，盡管它們在消色散網(wǎng)格、周圍射頻干擾的環(huán)境等存在顯著差異。因此，本文所述的特征閾值同樣能夠區(qū)分帕克斯數(shù)據(jù)庫中脈沖星信號和干擾。

我們將上述閾值和閾值區(qū)間應(yīng)用于單脈沖篩選工具，對脈沖星候選體進行重新篩選，試驗結(jié)果如表2。由表2可以看出，添加單個特征閾值判斷的單脈沖篩選工具能在一定程度上降低假陽性率，并且3個特征閾值組合判斷假陽性率從98.5%降低到78.4%，表明添加特征閾值判斷的單脈沖篩選工具提高了脈沖星搜索效率。

表2 添加特征后假陽性率結(jié)果

3 總結(jié)與展望

本文使用PRESTO的單脈沖搜索方法對CRAFTS超寬帶數(shù)據(jù)文件進行的試驗表明，現(xiàn)有的單脈沖篩選工具難以區(qū)分真實信號與噪聲或射頻干擾信號，造成數(shù)以萬計的假陽性樣本出現(xiàn)，顯著增加人工篩選單脈沖候選體的時間開銷以及候選數(shù)據(jù)存儲壓力。因此，我們提出3個顯著區(qū)分脈沖星信號和干擾的特征，并選取同時具有脈沖星信號、射頻干擾和宇宙噪聲的79顆脈沖星樣本，計算它們在3個特征取值情況。然后，根據(jù)脈沖星信號、射頻干擾和宇宙噪聲在3個特征取值分布的差異，提出合理閾值應(yīng)用于單脈沖篩選工具，用于對脈沖星候選體進行嚴(yán)格的判斷。最后，使用添加特征閾值判斷的篩選工具對脈沖星候選體進行重新試驗。結(jié)果表明，在保證所有脈沖星信號不遺漏的情況下，假陽性率從98.5%降低到78.4%。因此，本文所述3個特征具有實用性和有效性，有助于單脈沖搜索在CRAFTS巡天數(shù)據(jù)的應(yīng)用。

致謝：本文在500 m口徑球面射電望遠鏡(FAST)數(shù)據(jù)基礎(chǔ)上完成。FAST是由中國科學(xué)院國家天文臺運行和管理的國家大科學(xué)裝置。感謝中國科學(xué)院天文大科學(xué)研究中心FAST重大成果培育項目對本文工作的資助。