脈沖星信號(hào)相干消色散與非相干消色散的比較研究*

黃玉祥，汪 敏，郝龍飛，李志玄，徐永華

(1.中國科學(xué)院云南天文臺(tái)，云南昆明 650011；2.中國科學(xué)院大學(xué)，北京 100049)

脈沖星是一種具有超高溫、超高壓、超高密度、超強(qiáng)磁場(chǎng)、超強(qiáng)電場(chǎng)和超強(qiáng)引力場(chǎng)等極端物理?xiàng)l件的天體，一般認(rèn)為是超新星爆發(fā)的產(chǎn)物，典型半徑約為10 km，而質(zhì)量卻與太陽相當(dāng)，核心密度高達(dá)1014g/cm3[1]。目前，大多數(shù)人認(rèn)為，脈沖星是高速旋轉(zhuǎn)的磁中子星，其自轉(zhuǎn)軸與磁軸之間有一個(gè)夾角，兩個(gè)磁極各有一個(gè)輻射束。脈沖星的輻射幾乎占據(jù)了整個(gè)電磁波譜，包括射電、紅外、可見光、紫外、X射線和γ射線等電磁波頻段。其輻射具有較高的周期性，輻射周期和自轉(zhuǎn)周期一致，且周期變化率較小。

脈沖星極端的物理特性和輻射的高周期性具有較大的研究意義和應(yīng)用前景[2]。從科學(xué)研究來講，脈沖星可以用于驗(yàn)證和發(fā)展極端條件下的物理學(xué)和天文學(xué)理論，毫秒脈沖雙星的研究已經(jīng)檢驗(yàn)了廣義相對(duì)論和引力波輻射[3-4]；從工程應(yīng)用來講，毫秒脈沖星可以用于精確計(jì)時(shí)和星際導(dǎo)航[1，5]。由于毫秒脈沖星的周期變化率小于10-19，周期噪聲極小，周期穩(wěn)定度比較高，且隨時(shí)間的增加會(huì)持續(xù)改善，甚至能超過氫原子鐘，多個(gè)不同天區(qū)的毫秒脈沖星計(jì)時(shí)觀測(cè)可以組合成一個(gè)平均的脈沖星鐘，這可以用于現(xiàn)有原子鐘長期穩(wěn)定性的修正。

為了進(jìn)行上述研究，就要獲得精確的脈沖星脈沖到達(dá)時(shí)間[1](Time Of Arrival，TOA)，由于星際介質(zhì)的存在，觀測(cè)設(shè)備接收到的脈沖星信號(hào)發(fā)生色散。為了簡化理論，星際介質(zhì)被看成是低溫、低密度等離子體，它對(duì)在其中傳播的頻率為f的電磁波的折射率為

其中，fp為等離子體頻率ne取 0.03 cm-3[1]。

對(duì)于射電波段觀測(cè)，接收到的脈沖信號(hào)是一維時(shí)域上的脈沖輪廓，不同頻率的脈沖之間發(fā)生色散延遲，頻譜起始頻率信號(hào)相對(duì)于其截止頻率信號(hào)延遲時(shí)間[1]為

其中，DM(Dispersion Measure)為傳播路徑上的星際介質(zhì)的色散量，色散量單位為為傳播路徑上的平均電子數(shù)密度；D為星際介質(zhì)的色散常量，MHz2pc-1cm3s；flow和fhigh分別為接收頻帶的起始頻率和截止頻率，單位為MHz。

由(2)式可知，一定帶寬脈沖星信號(hào)的色散延遲會(huì)造成時(shí)域脈沖的相位移動(dòng)、幅度降低和脈沖輪廓展寬等色散現(xiàn)象，這限制了到達(dá)時(shí)間分辨率。因此需要對(duì)接收到的脈沖星信號(hào)進(jìn)行消色散并獲得精確的到達(dá)時(shí)間，通過對(duì)脈沖到達(dá)時(shí)間的擬合可以得到脈沖星的相關(guān)物理參數(shù)，進(jìn)而研究其自身的輻射機(jī)制和其他相關(guān)特性。

消色散的宗旨是消除信號(hào)頻帶內(nèi)的色散效應(yīng)，目前主要有非相干消色散和相干消色散[1]。非相干消色散是基于時(shí)間序列，用濾波器組分通道，按相對(duì)于參考通道的延遲時(shí)間單個(gè)通道進(jìn)行平移從而消除色散；相干消色散是在頻域上利用星際介質(zhì)的傳輸函數(shù)(chirp函數(shù))實(shí)現(xiàn)通帶內(nèi)的消色散，如圖1。

圖1 消色散算法流程圖；上半部分是相干消色散的流程圖，下半部分是非相干消色散的流程圖Fig.1 Flow chart of coherent de-dispersion algorithm；the above belongs to coherent de-dispersion，the below belongs to incoherent de-dispersion

目前已經(jīng)有學(xué)者在比較兩種消色散方法方面進(jìn)行了一些研究。文[6]在430 MHz處對(duì)PSRJ2322＋2057進(jìn)行了觀測(cè)，獲得了相干和非相干消色散后的平均脈沖輪廓，并計(jì)算出脈沖到達(dá)時(shí)間的精度，后者與前者到達(dá)時(shí)間的均方根誤差比值為0.3，如圖2。文[7]用數(shù)字信號(hào)處理(Digital Signal Processing，DSP)算法在327 MHz處比較了這兩種方法得到的一些脈沖星信號(hào)的積分輪廓，相干消色散后數(shù)據(jù)的時(shí)域序列脈沖幅度和時(shí)間分辨率是非相干消色散后數(shù)據(jù)的10～15倍[7]。國內(nèi)，中國科學(xué)院新疆天文臺(tái)的脈沖星觀測(cè)團(tuán)隊(duì)利用基于MK5A的VLBI記錄系統(tǒng)建立一套帶寬為64 MHz的相干消色散系統(tǒng)，得到了一些脈沖星的平均脈沖輪廓，與以前使用的2×128×2.5 MHz多通道非相干消色散系統(tǒng)的結(jié)果進(jìn)行了比較[8]，結(jié)論是相干消色散優(yōu)于非相干消色散，但沒有進(jìn)行系統(tǒng)和定量的分析研究。本文旨在用相關(guān)系數(shù)的方法定量地分析兩種方法消色散效果的差別。

1 消色散原理

兩種消色散方法基于不同的消色散原理。非相干消色散是基于時(shí)域上各個(gè)頻段的信號(hào)相位移動(dòng)進(jìn)行的不連續(xù)消色散，相干消色散是基于頻域乘積的連續(xù)消色散，具體描述如下。

1.1 非相干消色散

非相干消色散的過程是首先把觀測(cè)帶寬為BW的時(shí)域信號(hào)由濾波器組分成若干狹窄通道，然后在時(shí)域上把每個(gè)通道信號(hào)按適當(dāng)?shù)臅r(shí)延量進(jìn)行平移，最后將脈沖對(duì)齊后的所有通道信號(hào)幅度累加得到消色散后的時(shí)域信號(hào)序列，如圖3。非相干消色散的優(yōu)點(diǎn)是數(shù)據(jù)量小，計(jì)算量少。這種方法的缺點(diǎn)：(1)單個(gè)通道內(nèi)的色散效應(yīng)沒有消除，增大通道的數(shù)量可以減少剩余的色散效應(yīng)，但是由不確定性原理tres∝1/(2×bw)可知，增大通道數(shù)量，單個(gè)通道信號(hào)的時(shí)間分辨率會(huì)增大，所以要根據(jù)需求，適當(dāng)調(diào)整通道數(shù)的數(shù)值；(2)由于非相干消色散過程是單通道時(shí)域序列的平移，該方法會(huì)損失信號(hào)的相位信息。

圖2 脈沖星PSR J2322＋2057的脈沖輪廓。實(shí)線是相干消色散后的輪廓，虛線是非相干消色散后的輪廓[6]Fig.2 Pulse profile of PSR J2322＋2057.Solid line： pulse profile observed by coherent de-dispersion.Dashed line： The same profile， filtered by incoherent dedispersion

圖3 消色散原理圖[9]。上半部分是非相干消色散的原理，下半部分是相干消色散的原理Fig.3 Schematic diagram of incoherent de-dispersion.The upper part shows the principle of incoherent dedispersion，the lower part shows the Principle of coherent de-dispersion

1.2 相干消色散

星際介質(zhì)的色散效應(yīng)相當(dāng)于移相器的作用，因此接收到的脈沖星信號(hào)相當(dāng)于原始信號(hào)加上移相器得到的結(jié)果，移相器可以用傳遞函數(shù)H(f)表征，具體表達(dá)[1]為

其中，fo為本振頻率；f1為中頻頻率。

由傳遞函數(shù)H(f)可知，原始信號(hào)可以通過乘以其復(fù)共軛H(-f)進(jìn)行完全消色散，消除整個(gè)觀測(cè)帶寬內(nèi)的色散效應(yīng)，如圖3，最大的時(shí)間分辨率為1/2BW( )。這種消色散方法可用于高精度毫秒脈沖星計(jì)時(shí)和脈沖星射電輻射過程的研究。

理論上，兩種消色散結(jié)果的差異主要由各個(gè)通道的剩余色散量造成的累加效應(yīng)，因此分析兩種消色散效果可以通過分析各個(gè)通道的剩余色散量實(shí)現(xiàn)。

對(duì)于上邊帶信號(hào)的非相干消色散，由(2)式可知，整個(gè)帶寬BW內(nèi)的最大時(shí)延：

其中，BW為觀測(cè)帶寬。

單個(gè)通道的時(shí)間延遲：

2 脈沖星信號(hào)模擬、消色散算法及比較方法實(shí)現(xiàn)

本文用MATLAB程序模擬了射電望遠(yuǎn)鏡接收到的脈沖星信號(hào)，實(shí)現(xiàn)了兩種消色散過程和兩種消色散方法效果的定量比較過程，具體過程如下。

2.1 原始信號(hào)模擬

脈沖星輻射的原始信號(hào)很弱，為了提高信噪比，需要長時(shí)間的積分。為了簡化計(jì)算，此處模擬一個(gè)周期的基帶信號(hào)。為了便于比較兩種方法的消色散效果，采用短周期0.003 906 25 s，同時(shí)為了提高程序的運(yùn)算效率，采樣率(2 Gbps)、采樣時(shí)長和信號(hào)帶寬等都取2的冪指數(shù)。加色散前的基帶信號(hào)和加色散后基帶信號(hào)如圖4。色散的基帶信號(hào)模擬過程如下：

圖4 (a)加色散前的基帶模擬信號(hào)，(b)加色散后的基帶模擬信號(hào)；信號(hào)頻譜帶寬BW為1 GHz，起始頻率fo為1 GHz，色散量DM為40 cm-3·pc，圖5也取同樣的值Fig.4 The left image is the baseband analog signal before adding dispersion，the right image is the baseband signal after adding dispersion； Signal spectrum bandwidth BW： 1GHz， the starting frequency fo： 1GHz， the Dispersion Measure DM：40cm-3·pc.In the figure.5 those variables also take same value

(1)模擬帶寬為BW，起始頻率為fo，譜指數(shù)為-2.0的基帶信號(hào)[1]。

(2)把頻域的模擬基帶信號(hào)通過逆傅里葉變換到時(shí)域，進(jìn)行高斯脈沖調(diào)制。

(3)把調(diào)制后信號(hào)再通過傅里葉變換到頻域，乘以傳輸函數(shù)H(f)后進(jìn)行加色散。

(4)把加色散的頻域信號(hào)再變換到時(shí)域，加高斯白噪聲成模擬的脈沖星時(shí)域信號(hào)。

2.2 消色散算法實(shí)現(xiàn)

根據(jù)模擬長度為N的脈沖星信號(hào)數(shù)據(jù)，利用MATLAB程序?qū)崿F(xiàn)信號(hào)的消色散過程，并實(shí)現(xiàn)兩種消色散后的序列與加色散前序列的相關(guān)過程。

2.2.1 非相干消色散步驟

(1)數(shù)據(jù)分箱：把模擬的基帶數(shù)據(jù)按順序平均分箱。

(2)多相位濾波①https：//casper.berkeley.edu/wiki/The＿Polyphase＿Filter＿Bank＿Technique：對(duì)每個(gè)箱中的數(shù)據(jù)進(jìn)行多相位濾波，數(shù)據(jù)序列長度減小原始的這可以降低傅里葉變換的頻譜泄露[10]。多相位濾波過程為：利用相同長度的sinc函數(shù)序列依次乘以每個(gè)箱中的數(shù)據(jù)，把得到的序列按順序平均分成4段，并逐次對(duì)應(yīng)相加得到原始箱內(nèi)序列長度數(shù)據(jù)。

(3)傅里葉變換：對(duì)每個(gè)箱中的數(shù)據(jù)進(jìn)行傅里葉變換。

(4)分通道并平移數(shù)據(jù)：取參考頻率為截止頻率fhigh，任意中心頻率f的通道按延遲時(shí)間平移。由(2)式可知，相對(duì)延遲時(shí)間為

(5)逆傅里葉變換并求和：把每個(gè)箱中數(shù)據(jù)再進(jìn)行逆傅里葉變換，然后對(duì)箱內(nèi)的數(shù)據(jù)求和得到分箱后的數(shù)據(jù)序列，非相干消色散后信號(hào)的二維頻譜見圖5(a)。

圖5 (a)非相干消色散后信號(hào)的二維頻譜；(b)相干消色散后信號(hào)的二維頻譜Fig.5 The two-dimensional spectrum of the signal after incoherent de-dispersion，the two-dimensional spectrum of the signal after coherent de-dispersion

2.2.2 相干消色散步驟

基帶信號(hào)只有一個(gè)周期，為了簡化運(yùn)算，將每次處理的數(shù)據(jù)長度也定為一個(gè)周期。具體步驟如下：

(1)直接進(jìn)行傅里葉變換：把模擬的基帶數(shù)據(jù)進(jìn)行傅里葉變換，得到頻域數(shù)據(jù)。

(2)頻域消色散：把頻域數(shù)據(jù)乘以消色散函數(shù)H(-f)得到消色散后的頻域數(shù)據(jù)。

(3)逆傅里葉變換：把消色散后的頻域數(shù)據(jù)進(jìn)行逆傅里葉變換得到原始序列長度的消色散后的時(shí)域數(shù)據(jù)。

(4)折疊成分箱序列：把消色散后的數(shù)據(jù)依次折疊成與非相干消色散分箱后序列長度相同的箱，然后把每個(gè)箱中的數(shù)據(jù)累加得到序列，相干消色散后信號(hào)的二維頻譜見圖5(b)。

2.3 消色散效果的定量比較

為了定量分析兩種消色散，讓兩種消色散后的結(jié)果與加色散前的數(shù)據(jù)進(jìn)行互相關(guān)，求它們的互相關(guān)系數(shù)：

其中，xi屬于長度為N的序列X；yi屬于長度為N的序列Y。目前普遍認(rèn)為，相關(guān)系數(shù)的絕對(duì)值在0.3以下代表無直線相關(guān)；0.3以上代表直線相關(guān)；0.3～0.5代表低度相關(guān)；0.5～0.8代表顯著相關(guān)(中等程度相關(guān))；0.8以上代表高度相關(guān)；當(dāng)且僅當(dāng)其絕對(duì)值為1時(shí)，代表線性相關(guān)。對(duì)于研究的信號(hào)，表明兩個(gè)信號(hào)完全相同，如相同的脈寬和信噪比。具體步驟為：

(1)把加色散前的模擬數(shù)據(jù)按上面的方法折疊累加成長度與非相干消色散分箱后序列長度相同的數(shù)據(jù)序列。

(2)把兩種消色散后的數(shù)據(jù)與上一步的分箱數(shù)據(jù)進(jìn)行互相關(guān)得到互相關(guān)系數(shù)。

3 兩種消色散的定量比較

首先，研究用上述兩組相關(guān)系數(shù)定量比較兩組消色散的效果。然后，估計(jì)這兩種相關(guān)系數(shù)小于等于0.01時(shí)信號(hào)頻譜的起始頻率。

在信號(hào)模擬和消色散過程中，色散量的取值為1～40 cm-3·pc，步長為1 cm-3·pc。信號(hào)頻譜帶寬BW取值為32 MHz、64 MHz、128 MHz、256 MHz、512 MHz和1 GHz，單通道帶寬為1 MHz。信號(hào)頻譜起始頻率的取值范圍為256 MHz～8 GHz，比較運(yùn)算耗時(shí)和消色散效果，其步長為64 MHz，對(duì)于確定兩種消色散效果相同時(shí)的頻譜起始頻率，其步長為16 MHz。

3.1 兩種消色散結(jié)果得到的互相關(guān)系數(shù)的比較

通過計(jì)算得到兩種消色散后的數(shù)據(jù)與加色散前數(shù)據(jù)的相關(guān)系數(shù)如圖6、圖7。首先分析兩種消色散結(jié)果得到互相關(guān)系數(shù)隨各變量的變化特點(diǎn)。

圖6 相干消色散后的數(shù)據(jù)與加色散前的數(shù)據(jù)的相關(guān)系數(shù)；為了便于顯示，此處相關(guān)系數(shù)乘了1 000倍Fig.6 The coherence coefficient of the data after the coherent de-dispersion with the data before the adding dispersion.To facilitate the display，the coherence coefficient is multiplied by 1000

由圖6可知，相干消色散得到的相關(guān)系數(shù)接近于1，這說明，在不考慮離散化運(yùn)算過程中引入的誤差時(shí)，相干消色散可以完全恢復(fù)原始加色散之前的信號(hào)。結(jié)合兩圖可以發(fā)現(xiàn)，當(dāng)BW減小時(shí)，相關(guān)系數(shù)的波動(dòng)變大，這與(6)式有關(guān)；其次，數(shù)據(jù)處理過程存在頻譜泄露，這是由于數(shù)據(jù)在周期性分箱過程中產(chǎn)生一些截?cái)囝l點(diǎn)，有限長離散數(shù)據(jù)的傅里葉變換有頻譜泄露，這些效應(yīng)的影響隨著BW的減小而增大。

兩種消色散得到的相關(guān)系數(shù)之差見圖7。對(duì)于相同的色散量DM，頻譜起始頻率fo的取值在一定范圍內(nèi)，相關(guān)系數(shù)之差隨著變量fo指數(shù)減小，但到某值f1之上，該相關(guān)系數(shù)之差接近于1且波動(dòng)逐漸變小，這說明在頻譜起始頻率fo低于f1時(shí)，非相干消色散的效果比相干消色散的差，且兩種方法的消色散效果的差距隨著fo的增加而縮小。在fo高于f1時(shí)，兩種方法的消色散效果相同。結(jié)合兩圖可以看出，相同fo隨DM變化對(duì)應(yīng)的相關(guān)系數(shù)的變化曲線隨著BW的減小起伏變大，這是數(shù)據(jù)處理時(shí)頻譜泄露的誤差造成的。相同fo、相同DM對(duì)應(yīng)的相關(guān)系數(shù)之差隨BW的增大而減小，且f1隨BW的減小而增大。結(jié)合(6)式和圖7可知，在不考慮誤差范圍時(shí)，在一定觀測(cè)頻率內(nèi)，相干消色散效果優(yōu)于非相干消色散效果，兩種方法的消色散效果的差距隨觀測(cè)帶寬、觀測(cè)頻率的增加而近似按指數(shù)關(guān)系減小，隨著色散量的增大而近似按冪指數(shù)關(guān)系增大，同時(shí)能夠確定，可以找到兩種方法消色散效果相同時(shí)的觀測(cè)起始頻率。

圖7 兩種消色散方法得到的相關(guān)系數(shù)之差Fig.7 The coherence coefficient of the data after the incoherent de-dispersion with the data before the adding dispersion

3.2 兩種消色散效果相同對(duì)應(yīng)的通帶起始頻率

兩種消色散效果相同時(shí)，這對(duì)應(yīng)于非相干消色散單通道內(nèi)的時(shí)延tdely取某一接近零的常數(shù)。非相干消色散過程，若單個(gè)通道帶寬bw為常數(shù)，由(6)式可得

由上式，tdely為常數(shù)，有

經(jīng)過計(jì)算，在兩種消色散后的數(shù)據(jù)與加色散前的數(shù)據(jù)的相關(guān)系數(shù)小于某一定值條件下，得到頻譜起始頻率fo如表1，考慮到如上所述數(shù)據(jù)處理中引入的誤差，這個(gè)值取0.01。

表1 不同BW、DM對(duì)應(yīng)消色散效果相同時(shí)的f o，其單位是MHz，這只是部分?jǐn)?shù)據(jù)Table 1 The fo to different BW and DM， which corresponds to the de-dispersion effect，the unit is MHz，This is just part of the data

4 結(jié) 論

運(yùn)用MATLAB程序?qū)崿F(xiàn)了兩種消色散過程的定量比較，并確定了兩種消色散方法相同時(shí)的頻譜起始頻率。在實(shí)現(xiàn)比較過程中，色散量的取值為1～40 cm-3·pc，步長為1 cm-3·pc；信號(hào)頻譜帶寬BW取值為32 MHz、64 MHz、128 MHz、256 MHz、512 MHz和1 GHz，單通道帶寬為1 MHz；信號(hào)頻譜起始頻率的取值范圍為256 MHz～8 GHz。對(duì)于比較運(yùn)算耗時(shí)和消色散效果，其步長為64 MHz，對(duì)于確定兩種消色散效果相同時(shí)的頻譜起始頻率，其步長為16 MHz。定量的比較結(jié)果說明，在一定的觀測(cè)頻率之內(nèi)，相干消色散的效果優(yōu)于非相干消色散，實(shí)際觀測(cè)中可以根據(jù)觀測(cè)要求，參考這些結(jié)果，選擇消色散方法。最后研究了在不同色散量和不同頻譜帶寬下，兩種消色散技術(shù)消色散效果相同時(shí)的起始頻率，非相干消色散脈沖星觀測(cè)系統(tǒng)在實(shí)際的觀測(cè)中可以參考這些頻率制定觀測(cè)頻段。非相干消色散過程引用了多相位濾波器，一定程度上緩解了傅里葉變換的頻譜泄露問題。但是頻譜泄露問題依然存在，尤其在低頻端和高頻端。下一步將擴(kuò)展色散量的變化范圍和變化的單個(gè)通道帶寬，同時(shí)繼續(xù)改善頻譜泄漏問題。

致謝：由衷感謝北京大學(xué)李柯伽老師和國家天文臺(tái)盧吉光博士對(duì)脈沖星信號(hào)模擬的指導(dǎo)。