脈沖雙星加速度搜索方法及軟件改進(jìn) *

冀昶旭，于徐紅，劉志杰

(貴州師范大學(xué)貴州省信息與計(jì)算科學(xué)重點(diǎn)實(shí)驗(yàn)室，貴州 貴陽(yáng) 550001)

脈沖雙星系統(tǒng)的發(fā)現(xiàn)極具科學(xué)意義。1974年，赫爾斯和泰勒在天鷹座天域探測(cè)到脈沖信號(hào)，經(jīng)計(jì)算發(fā)現(xiàn)是脈沖雙星[1](PSR B1913 + 16)的信號(hào)。兩人通過(guò)深入研究首次發(fā)現(xiàn)引力波存在的間接定量證據(jù)[2]，并因此獲得1993年諾貝爾物理學(xué)獎(jiǎng)。文[3]基于500 m口徑球面射電望遠(yuǎn)鏡數(shù)據(jù)在M13星團(tuán)中發(fā)現(xiàn)一顆處于雙星系統(tǒng)中的毫秒脈沖星M13F，并認(rèn)證星團(tuán)中另一顆脈沖星M13E為掩食雙星；首次測(cè)量了M13星團(tuán)中4個(gè)脈沖雙星系統(tǒng)的軌道橢率，同時(shí)獲得M13星團(tuán)現(xiàn)有脈沖星的最好計(jì)時(shí)結(jié)果。遺憾的是，現(xiàn)有脈沖雙星的搜索較為耗時(shí)，搜索能力亟待提升。一方面是FAST早期科學(xué)數(shù)據(jù)中心對(duì)數(shù)據(jù)存儲(chǔ)和處理能力的要求高。在FAST開(kāi)啟19波束巡天漂移掃描后，數(shù)據(jù)量激增，每天產(chǎn)生壓縮數(shù)據(jù)50 TB，每年(按200天計(jì)算)產(chǎn)生10 PB數(shù)據(jù)，如何利用好這些數(shù)據(jù)，是對(duì)數(shù)據(jù)中心存儲(chǔ)和超算能力的嚴(yán)峻考驗(yàn)。 另一方面，高效的搜索能力對(duì)FAST探測(cè)優(yōu)質(zhì)毫秒脈沖星具有重要意義。在FAST早期科學(xué)數(shù)據(jù)中心脈沖星分布式搜索計(jì)算系統(tǒng)Craber[4]中使用加速度搜索方法[5-7]，當(dāng)參數(shù)zmax設(shè)置較大時(shí)，搜索效率顯著降低。目前，Ransom版本的加速度搜索方法[8]雖然在中央處理器上進(jìn)行了高度優(yōu)化，但速度仍有欠缺。單機(jī)八核處理大小約140 G的巡天漂移掃描數(shù)據(jù)，zmax參數(shù)預(yù)設(shè)200，耗時(shí)約20 h，這僅是10 min巡天且經(jīng)過(guò)處理的數(shù)據(jù)量。顯然，這種數(shù)據(jù)處理能力極大限制了脈沖雙星的有效發(fā)現(xiàn)。

此外，雙星系統(tǒng)中蘊(yùn)含大量毫秒脈沖星。如表1和圖1中數(shù)據(jù)，球狀星團(tuán)中包含大量脈沖雙星，且多為毫秒脈沖星。截至2020年底，ATNF脈沖星數(shù)據(jù)庫(kù)共收錄脈沖雙星149顆[9]，其中，約有130顆在球狀星團(tuán)中發(fā)現(xiàn)，表1給出了部分球狀星團(tuán)中樣本數(shù)據(jù)統(tǒng)計(jì)結(jié)果，圖1統(tǒng)計(jì)了148顆脈沖星的周期分布(有1顆周期為2.7 s，統(tǒng)計(jì)意義不大，未記入)，有113顆屬于毫秒脈沖星，占雙星系統(tǒng)的四分之三，約75.8%，其中，周期為2～4 ms的毫秒脈沖星約占總數(shù)的60%。

根據(jù)圖1的統(tǒng)計(jì)結(jié)果，雙星系統(tǒng)中毫秒脈沖星數(shù)量居多。文[10]首次提出毫秒脈沖星是正常脈沖星通過(guò)在低質(zhì)量X射線雙星階段吸積物質(zhì)而加速到毫秒量級(jí)。文[11]也認(rèn)為當(dāng)?shù)痛艌?chǎng)的中子星吸積物質(zhì)加速到毫秒級(jí)后會(huì)演化為雙星系統(tǒng)。文[12]模擬的P-P˙圖展示了雙星脈沖星群自轉(zhuǎn)周期峰值在5 ms左右，并將雙星脈沖星群稱為毫秒脈沖星群。由此可知，雙星系統(tǒng)與毫秒脈沖星存在密切關(guān)系。目前，針對(duì)雙星系統(tǒng)的搜索較為耗時(shí)，因此，提高脈沖雙星的搜索能力十分必要。

1 脈沖雙星系統(tǒng)與加速度搜索

1.1 雙星系統(tǒng)及搜索技術(shù)

脈沖雙星系統(tǒng)是相互繞行的兩顆中子星，其中一顆是脈沖星，另一顆稱為伴星[13]。通常情況下，當(dāng)星體自轉(zhuǎn)且磁極波束掃過(guò)信號(hào)探測(cè)設(shè)備時(shí)，探測(cè)設(shè)備能接收到一個(gè)脈沖信號(hào)[14]。傳統(tǒng)的周期搜索是對(duì)采樣信號(hào)進(jìn)行快速傅里葉變換，將時(shí)域信號(hào)轉(zhuǎn)換到頻域，再通過(guò)計(jì)算頻譜功率來(lái)尋找周期。在雙星系統(tǒng)中，脈沖星受到伴星的引力作用，產(chǎn)生具有加速度的軌道運(yùn)動(dòng)，受其影響，脈沖星相對(duì)信號(hào)探測(cè)設(shè)備的視向速度發(fā)生變化[8]。由于多普勒效應(yīng)，探測(cè)設(shè)備接收的星體自旋頻率發(fā)生變化[15]，傳統(tǒng)的周期搜索方法不再適用，也難以檢測(cè)到脈沖星。加速度搜索方法[5-7]可以有效消除這種影響，采用恒定加速度近似軌道運(yùn)動(dòng)的加速度，從而消除觀測(cè)數(shù)據(jù)中由于雙星軌道運(yùn)動(dòng)導(dǎo)致的脈沖到達(dá)時(shí)間的變化[16]。在加速度恒定的假設(shè)下，脈沖星的自旋頻率在觀察者參考系中的漂移速率f˙與加速度a的關(guān)系為

(1)

其中，c為光速；f0為脈沖星在自身慣性系中的自旋頻率。根據(jù)信號(hào)處理角度的不同，加速度搜索可以分為時(shí)域加速度搜索(Time Domain Acceleration Search, TDAS)和頻域加速度搜索(Fourier Domain Acceleration Search, FDAS)[17]。

1.2 時(shí)域加速度搜索

時(shí)域加速度搜索通過(guò)線性展寬或壓縮信號(hào)的方式對(duì)數(shù)據(jù)進(jìn)行重新采樣，以補(bǔ)償軌道運(yùn)動(dòng)引起的多普勒效應(yīng)，并在一定范圍內(nèi)不斷嘗試加速度值，找到最接近真實(shí)值的時(shí)間序列，再利用周期搜索技術(shù)尋找周期信號(hào)。重采樣時(shí)信號(hào)的時(shí)間間隔為

(2)

其中，v為視向速度；c為光速；τ0為保證正確采樣的常數(shù)，

(3)

(3)式中，tsamp為采樣時(shí)間間隔；T為積分時(shí)間；a為嘗試的加速度值；c為光速。

脈沖星搜索軟件包SIGPROC[18]和PEASOUP[19]等應(yīng)用了信號(hào)重采樣技術(shù)。利用時(shí)域加速度搜索方法，文[20]在球狀星團(tuán)中找到了4顆毫秒脈沖星，觀測(cè)參數(shù)見(jiàn)表2。

1.3 頻域加速度搜索

頻域加速度搜索基于相關(guān)性技術(shù)[21]，可以提高頻譜功率。受軌道運(yùn)動(dòng)的影響，接收的脈沖星信號(hào)頻率發(fā)生變化，在頻譜中表現(xiàn)為某一頻率的功率擴(kuò)散至周圍相鄰的頻率點(diǎn)，導(dǎo)致普通的周期搜索方法失效。文[21]提到，可以將原始信號(hào)與鄰近信號(hào)進(jìn)行相關(guān)操作，從而確定相似程度以恢復(fù)信號(hào)，這個(gè)過(guò)程可表示為

(4)

通常認(rèn)為，在觀測(cè)時(shí)間遠(yuǎn)小于軌道周期的情況下使用頻域加速度搜索方法效果更好，即滿足

(5)

其中，Tobs為脈沖星觀測(cè)時(shí)間；Porb為脈沖星軌道運(yùn)動(dòng)周期。當(dāng)觀測(cè)時(shí)間和軌道周期滿足(5)式時(shí)，加速度可以視為常數(shù)，結(jié)合多普勒公式有

(6)

其中，c為光速；T為脈沖星觀測(cè)時(shí)間；f為諧波頻率；z為假設(shè)頻率漂移的頻率點(diǎn)數(shù)量。這樣，根據(jù)漂移的頻率點(diǎn)數(shù)z，通過(guò)局部信號(hào)傅里葉響應(yīng)的互相關(guān)恢復(fù)原始信號(hào)。

其他周期搜索方法也可以消除雙星系統(tǒng)中軌道運(yùn)動(dòng)的影響，但應(yīng)用條件不同。相位調(diào)制搜索主要針對(duì)觀測(cè)時(shí)間與雙星軌道周期大致相同的情況，邊帶搜索在觀測(cè)時(shí)間遠(yuǎn)大于軌道周期時(shí)效果更好。而加速度搜索的使用條件決定了它的適用目標(biāo)為周期更小的毫秒脈沖星，在實(shí)際中也搜尋到許多優(yōu)質(zhì)脈沖星。文[8]首次使用頻域加速度搜索方法在球狀星團(tuán)NGC 6544中發(fā)現(xiàn)了脈沖星PSR J1807-2459A。文[22]利用頻域搜索技術(shù)發(fā)現(xiàn)了高度偏心的毫秒脈沖星PSR J1946+3417，偏心率達(dá)到0.13～0.14，這極為罕見(jiàn)，同時(shí)也預(yù)示著脈沖星不同尋常的形成過(guò)程。

2 頻域加速度搜索并行化的必要性及可行性分析

針對(duì)頻域加速度搜索的并行化提速十分必要。文[16]論述了軌道運(yùn)動(dòng)對(duì)脈沖到達(dá)時(shí)間影響很小，在一定條件下不用加速度搜索也可以發(fā)現(xiàn)雙星。但2020年2月FAST在發(fā)現(xiàn)一顆雙星的過(guò)程中，文[3]將zmax參數(shù)設(shè)置為300才搜索到脈沖星PSR J1641+3627F。可見(jiàn)，隨著搜索要求的不斷提高，擴(kuò)大加速度的搜索范圍成為必然，否則數(shù)據(jù)無(wú)法得到有效處理，可能錯(cuò)失優(yōu)質(zhì)的脈沖星。

相比時(shí)域加速度搜索，頻域加速度搜索更快且更適用于并行化。時(shí)域加速度搜索將處理過(guò)的時(shí)間序列全部加載到內(nèi)存，而頻域加速度搜索處理的是局部序列的卷積，更適合小內(nèi)存的圖形處理器進(jìn)行并行化。對(duì)不同加速度值進(jìn)行處理時(shí)，時(shí)域加速度搜索需要對(duì)每個(gè)加速度重采樣后的序列進(jìn)行快速傅里葉變換，而頻域加速度搜索只需做一次快速傅里葉變換，再進(jìn)行相關(guān)操作即可。同時(shí)，頻域加速度搜索可以看作是一個(gè)濾波過(guò)程，本質(zhì)是原始信號(hào)與濾波器的卷積操作。原始信號(hào)長(zhǎng)度遠(yuǎn)大于濾波器長(zhǎng)度，計(jì)算的時(shí)間復(fù)雜度約為O(N2)，N為信號(hào)長(zhǎng)度。為提高計(jì)算效率，本文使用快速傅里葉變換做循環(huán)卷積。首先將原始信號(hào)拆分成與濾波器等長(zhǎng)度的信號(hào)(濾波器中需要補(bǔ)0)，再分別與濾波器做卷積并傅里葉逆變換到時(shí)域，最后將結(jié)果中重疊的部分截去，再合并為輸出信號(hào)。這一過(guò)程也稱重疊保留法，可以將算法的時(shí)間復(fù)雜度降至O(NlgN)，如圖2。濾波器h中間補(bǔ)充一定長(zhǎng)度的0元素，H為其傅里葉變換。將輸入信號(hào)x分段，每段長(zhǎng)為M，每段信號(hào)xi向前多取N/2個(gè)點(diǎn)，則第1段前面需補(bǔ)充N/2個(gè)0，再對(duì)每段信號(hào)做傅里葉變換得到Xi，Xi與H在復(fù)數(shù)域中相乘再進(jìn)行傅里葉逆變換得到A，B和C，將A，B和C中混疊的部分截去得到A′，B′和C′，拼接后得到等價(jià)的線性卷積，流程如圖3(a)。綜上所述，在頻域加速度搜索過(guò)程中使用循環(huán)卷積分段計(jì)算所具備的局部?jī)?nèi)存性是高度可并行的，非常適合并行化計(jì)算。

圖2 頻域加速度搜索中應(yīng)用的循環(huán)卷積Fig.2 Circular convolution in FDAS

3 頻域加速度搜索并行化的實(shí)現(xiàn)

根據(jù)第2節(jié)的介紹，程序耗時(shí)的矛盾主要在于大量可并行的卷積計(jì)算。一方面，我們通過(guò)循環(huán)卷積來(lái)局部縮短時(shí)間序列的長(zhǎng)度，從而降低內(nèi)存壓力。另一方面，循環(huán)卷積中分段后的序列可以利用圖形處理器計(jì)算單元的局部?jī)?nèi)存完成運(yùn)算。本文實(shí)現(xiàn)的運(yùn)算流程見(jiàn)圖3(b)。

圖3 (a)線性卷積運(yùn)算流程；(b)本文中循環(huán)卷積運(yùn)算流程Fig.3 (a) The flow chart of linear convolution; (b) the process of cyclic convolution in this article

在對(duì)循環(huán)分段的時(shí)間序列進(jìn)行傅里葉變換時(shí)，本文使用CUDA庫(kù)封裝好的cuFFT接口，因?yàn)閏uFFT接口已針對(duì)NVIDIA 圖形處理器進(jìn)行了快速傅里葉變換的高度優(yōu)化。利用批處理cuFFT例程，在圖形處理器上并行調(diào)用一次即完成多段序列的傅里葉變換。處理加速度濾波器時(shí)，我們?cè)谝粋€(gè)圖形處理器網(wǎng)格中加載濾波器，其中每個(gè)線程加載全部濾波器，然后分段處理經(jīng)過(guò)傅里葉變換后的信號(hào)，分別與不同的分段頻域信號(hào)進(jìn)行復(fù)數(shù)域的乘法運(yùn)算。這一過(guò)程在局部進(jìn)行，且每個(gè)線程同步處理相同的運(yùn)算，再將運(yùn)算結(jié)果寫(xiě)回全局內(nèi)存進(jìn)行傅里葉逆變換和重疊段的去除，合并后輸出。

4 實(shí)驗(yàn)及結(jié)果

實(shí)驗(yàn)使用AMD Ryzen7 2700X 8核16線程中央處理器，兩張GeForce GTX1080顯卡，對(duì)漂移掃描數(shù)據(jù)的若干文件(PSR20180811_00A01H， PSR20180826_00A01H， PSR20180828_00A01H和PSR20180902_00A01H， 觀測(cè)參數(shù)見(jiàn)表3)消色散后產(chǎn)生的 .dat文件進(jìn)行測(cè)試，使用不同zmax參數(shù)對(duì)上述4個(gè)文件進(jìn)行處理，產(chǎn)生相應(yīng)標(biāo)記周期信號(hào)的ACCEL文本文件，并與CPU版本、GPU版本(https://github.com/jintaoluo/presto_on_gpu)的加速度搜索做了對(duì)比，結(jié)果見(jiàn)圖4。

通過(guò)不同方法的比較，結(jié)合圖4中不同zmax參數(shù)對(duì)4個(gè)超大漂移掃描數(shù)據(jù)文件的實(shí)驗(yàn)處理，結(jié)果表明，采用文中方法對(duì)加速度搜索程序改進(jìn)后，本文并行化處理方法耗時(shí)相比Ransom版本有極大提升，加速比約10～12倍，較GPU版本也有一定進(jìn)步，提高約1.5～1.7倍。綜合來(lái)看，程序的改進(jìn)有一定成效，達(dá)到預(yù)期效果。

表3 本次實(shí)驗(yàn)使用的觀測(cè)文件參數(shù)Table 3 Observed parameters used in this experiment

圖4 不同方法使用不同zmax參數(shù)對(duì)執(zhí)行4個(gè)文件的消耗時(shí)間的對(duì)比結(jié)果。(a) PSR20180811_00A01H;(b) PSR20180826_00A01H; (c) PSR20180828_00A01H; (d) PSR20180902_00A01H