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

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

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

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

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

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

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

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

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

消色散的宗旨是消除信號頻帶內的色散效應，目前主要有非相干消色散和相干消色散[1]。非相干消色散是基于時間序列，用濾波器組分通道，按相對于參考通道的延遲時間單個通道進行平移從而消除色散；相干消色散是在頻域上利用星際介質的傳輸函數(shù)(chirp函數(shù))實現(xiàn)通帶內的消色散，如圖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

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

1 消色散原理

兩種消色散方法基于不同的消色散原理。非相干消色散是基于時域上各個頻段的信號相位移動進行的不連續(xù)消色散，相干消色散是基于頻域乘積的連續(xù)消色散，具體描述如下。

1.1 非相干消色散

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

圖2 脈沖星PSR J2322＋2057的脈沖輪廓。實線是相干消色散后的輪廓，虛線是非相干消色散后的輪廓[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 相干消色散

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

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

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

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

對于上邊帶信號的非相干消色散，由(2)式可知，整個帶寬BW內的最大時延：

其中，BW為觀測帶寬。

單個通道的時間延遲：

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

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

2.1 原始信號模擬

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

圖4 (a)加色散前的基帶模擬信號，(b)加色散后的基帶模擬信號；信號頻譜帶寬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的基帶信號[1]。

(2)把頻域的模擬基帶信號通過逆傅里葉變換到時域，進行高斯脈沖調制。

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

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

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

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

2.2.1 非相干消色散步驟

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

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

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

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

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

圖5 (a)非相干消色散后信號的二維頻譜；(b)相干消色散后信號的二維頻譜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 相干消色散步驟

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

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

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

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

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

2.3 消色散效果的定量比較

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

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

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

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

3 兩種消色散的定量比較

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

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

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

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

圖6 相干消色散后的數(shù)據(jù)與加色散前的數(shù)據(jù)的相關系數(shù)；為了便于顯示，此處相關系數(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可知，相干消色散得到的相關系數(shù)接近于1，這說明，在不考慮離散化運算過程中引入的誤差時，相干消色散可以完全恢復原始加色散之前的信號。結合兩圖可以發(fā)現(xiàn)，當BW減小時，相關系數(shù)的波動變大，這與(6)式有關；其次，數(shù)據(jù)處理過程存在頻譜泄露，這是由于數(shù)據(jù)在周期性分箱過程中產生一些截斷頻點，有限長離散數(shù)據(jù)的傅里葉變換有頻譜泄露，這些效應的影響隨著BW的減小而增大。

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

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

3.2 兩種消色散效果相同對應的通帶起始頻率

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

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

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

表1 不同BW、DM對應消色散效果相同時的f o，其單位是MHz，這只是部分數(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 結 論

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

致謝：由衷感謝北京大學李柯伽老師和國家天文臺盧吉光博士對脈沖星信號模擬的指導。