多成分脈沖輪廓模版及其在脈沖星測(cè)時(shí)中的應(yīng)用

盧吉光，劉雨蘭，姜鵬

多成分脈沖輪廓模版及其在脈沖星測(cè)時(shí)中的應(yīng)用

盧吉光1,2，劉雨蘭1,2，姜鵬1,2

（1. 中國(guó)科學(xué)院 國(guó)家天文臺(tái)，北京 100101；2. 貴州射電天文臺(tái)，貴陽(yáng) 550000）

部分射電脈沖星存在模式變換現(xiàn)象，射電脈沖輪廓的變化給脈沖星測(cè)時(shí)工作帶來(lái)了一定的困難。本文介紹了對(duì)存在多種穩(wěn)定成分的脈沖輪廓進(jìn)行成分分離的方案，并通過(guò)模擬數(shù)據(jù)實(shí)驗(yàn)說(shuō)明了利用多成分脈沖標(biāo)準(zhǔn)輪廓進(jìn)行測(cè)時(shí)優(yōu)于使用單成分標(biāo)準(zhǔn)輪廓進(jìn)行測(cè)時(shí)。

脈沖星；脈沖星測(cè)時(shí)；模式變換；脈沖輪廓模版；主成分分析

0 引言

射電脈沖星的平均脈沖輪廓通常是穩(wěn)定的[1]，因此可以使用多次觀測(cè)獲得的脈沖輪廓進(jìn)行疊加獲得最終收斂的穩(wěn)定結(jié)果。利用獲得的穩(wěn)定結(jié)果或提取其形狀特征信息作為標(biāo)準(zhǔn)脈沖輪廓模版，可與單次觀測(cè)的平均脈沖輪廓進(jìn)行比對(duì)來(lái)獲取脈沖輻射到達(dá)觀測(cè)站的時(shí)間，從而進(jìn)行脈沖星測(cè)時(shí)。

然而某些脈沖星的射電脈沖輪廓是不穩(wěn)定的，會(huì)在多種模式中切換[2-6]。如果不同模式之間的特征差別較大，且可以根據(jù)單脈沖形狀劃分到不同的模式中，則不同模式的脈沖輪廓可以對(duì)不同模式下的單脈沖進(jìn)行分別平均計(jì)算獲得，同時(shí)也可以按照不同模式下的脈沖星輪廓分別進(jìn)行測(cè)時(shí)以獲得好的測(cè)時(shí)結(jié)果。但對(duì)于絕大部分發(fā)生模式變換的脈沖星，要么輻射強(qiáng)度不足以分辨單脈沖形狀，要么單脈沖本身變化較大難以區(qū)分其屬于哪種輻射模式。因此對(duì)于存在模式變換的脈沖星，特定模式下的高精度的標(biāo)準(zhǔn)脈沖輪廓與觀測(cè)脈沖輪廓難以獲得，而對(duì)于單個(gè)模式持續(xù)時(shí)間較短的脈沖星這一問(wèn)題尤其明顯。這為高精度脈沖星測(cè)時(shí)工作帶來(lái)了極大挑戰(zhàn)。

本文將介紹一種通過(guò)數(shù)學(xué)方法獲得含有多種成分的脈沖輪廓模版，并使用多成分輪廓模版對(duì)脈沖星進(jìn)行測(cè)時(shí)的手段，從而有望提高具有多種穩(wěn)定輪廓成分的脈沖星的測(cè)時(shí)結(jié)果的精度。

1 脈沖輪廓成分分離

本節(jié)將介紹對(duì)含有多種輻射模式的脈沖星觀測(cè)輪廓進(jìn)行成分分離的相關(guān)方案。

1.1 脈沖輪廓成分介紹

1.2 經(jīng)典主成分分析方法

經(jīng)典主成分分析方法[7]常被用在天文數(shù)據(jù)分析中，比如斯隆望遠(yuǎn)鏡數(shù)字巡天項(xiàng)目（Sloan Digital Sky Survey，SDSS）中的類星體光譜分類[8]和弱引力透鏡現(xiàn)象中的點(diǎn)擴(kuò)散函數(shù)變化分析[9]等都使用了經(jīng)典主成分分析方法。

如果單次脈沖星觀測(cè)中包含了多種輻射模式，且無(wú)法將多種輻射模式精確分離，則可用自助抽樣法（bootstrap方法）從觀測(cè)數(shù)據(jù)中抽選指定數(shù)目的單脈沖疊加形成多個(gè)脈沖平均輪廓，而每個(gè)平均輪廓都可以視為不同模式的輪廓成分的線性組合。一般情況下，利用長(zhǎng)期測(cè)時(shí)獲得的脈沖星星歷，在用于對(duì)單次脈沖星觀測(cè)數(shù)據(jù)進(jìn)行折疊時(shí)，不會(huì)在不同周期中導(dǎo)致大的相位偏差，因此可以將此時(shí)獲得的不同脈沖輪廓的脈沖相位視為已對(duì)齊的。即公式（1）可簡(jiǎn)化為

。（5）

1.3 一般情況

經(jīng)典主成分分析方法對(duì)應(yīng)的求解方式簡(jiǎn)潔明了，但在很多情況下的數(shù)據(jù)并不滿足其適用條件。比如蓋亞探測(cè)器（Gaia）的光譜數(shù)據(jù)中，不同數(shù)據(jù)的誤差不同，因此在進(jìn)行成分分析時(shí)需要考慮不同數(shù)據(jù)有著不同的權(quán)重，即帶權(quán)重的主成分分析（weighted PCA）[10]；在聲學(xué)領(lǐng)域，很多情況下需要分析的信號(hào)不存在確定的起始時(shí)間，因此需要采用時(shí)移性主成分分析（time-shifted PCA）[11-12]。在實(shí)際的成分分析過(guò)程中，需要根據(jù)實(shí)際情況來(lái)選取成分分析方法。

2 數(shù)據(jù)模擬多成分脈沖輪廓模版在脈沖星測(cè)時(shí)中的應(yīng)用

本節(jié)將使用數(shù)值模擬的方式展示如何在具有模式變換現(xiàn)象的脈沖星數(shù)據(jù)中提取多成分脈沖輪廓模版，并將它們應(yīng)用在脈沖星測(cè)時(shí)中。

本文中所進(jìn)行的各種數(shù)據(jù)模擬以及脈沖星測(cè)時(shí)等操作均為使用python3編程進(jìn)行，其中繪圖部分采用matplotlib模塊中的相關(guān)函數(shù)完成，涉及到函數(shù)擬合的部分采用的是scipy.optimize模塊中的leastsq函數(shù)，以Levenberg-Marquardt算法完成。詳細(xì)的數(shù)據(jù)模擬過(guò)程如下。

首先，需要模擬存在多種輻射模式的脈沖輪廓。在本文的模擬中，假定脈沖星存在3種不同的輻射模式，生成模擬脈沖輪廓數(shù)據(jù)的方式如下：① 假定數(shù)據(jù)中的脈沖星存在3種輻射模式，不同模式下的脈沖輪廓如圖1所示；② 考慮不同觀測(cè)數(shù)據(jù)中3種輻射的強(qiáng)度在0～1中均勻隨機(jī)分布，即每次模擬觀測(cè)數(shù)據(jù)中的脈沖輪廓為3種模式下脈沖輪廓乘以0～1之間的均勻分布隨機(jī)浮點(diǎn)數(shù)再相加；③ 設(shè)定不同位置處的系統(tǒng)噪聲強(qiáng)度為0.015，由脈沖星輻射造成的流量起伏為3種輻射模式脈沖輪廓強(qiáng)度之和的0.02倍，利用高斯分布的隨機(jī)浮點(diǎn)數(shù)生成器生成不同相位處的系統(tǒng)噪聲與流量起伏；④ 設(shè)定不同觀測(cè)中的相位移動(dòng)為0～1之間的均勻分布隨機(jī)浮點(diǎn)數(shù)，將模擬數(shù)據(jù)與系統(tǒng)噪聲和流量起伏相加，然后將相位按設(shè)定值進(jìn)行平移，獲得模擬脈沖輪廓數(shù)據(jù)。圖2中展示的是模擬數(shù)據(jù)中的16個(gè)脈沖輪廓?？梢钥吹?，不同觀測(cè)中的脈沖輪廓具有不同的形狀與脈沖相位。

圖1 模擬數(shù)據(jù)中3種輻射模式下的平均脈沖輪廓

圖2 模擬數(shù)據(jù)中的16個(gè)觀測(cè)脈沖輪廓示例

為了檢驗(yàn)新的多成分脈沖標(biāo)準(zhǔn)輪廓在脈沖星測(cè)時(shí)過(guò)程中的應(yīng)用，我們需要比較這一標(biāo)準(zhǔn)輪廓組在測(cè)時(shí)時(shí)生成的結(jié)果與使用單成分標(biāo)準(zhǔn)輪廓對(duì)脈沖星進(jìn)行測(cè)時(shí)的結(jié)果。首先，將模擬觀測(cè)脈沖輪廓數(shù)據(jù)對(duì)齊后疊加生成單成分的脈沖標(biāo)準(zhǔn)輪廓（對(duì)照組，如圖3中點(diǎn)線所示）；接下來(lái)，我們利用本節(jié)第一段中的方法重新生成了一組模擬脈沖輪廓，然后考察利用單成分脈沖標(biāo)準(zhǔn)輪廓與多成分脈沖輪廓分別對(duì)這一組新生成的模擬脈沖輪廓進(jìn)行比對(duì)，獲得這一組脈沖輪廓與兩組脈沖標(biāo)準(zhǔn)輪廓之間的相位差，并將求得的相位差與生成脈沖輪廓時(shí)加入的相位移動(dòng)進(jìn)行對(duì)比。單成分脈沖標(biāo)準(zhǔn)輪廓測(cè)時(shí)使用的是經(jīng)典的相位梯度法計(jì)算標(biāo)準(zhǔn)輪廓與模擬脈沖輪廓的相位差；多成分脈沖標(biāo)準(zhǔn)輪廓的測(cè)時(shí)稍微復(fù)雜一些，需要將標(biāo)準(zhǔn)輪廓中的多個(gè)成分以不同的線性組合方式進(jìn)行疊加，并以相位梯度法計(jì)算每個(gè)疊加后的輪廓與模擬脈沖輪廓的相位差，最終保留計(jì)算誤差最小的一組結(jié)果。圖4中繪制的即為兩種不同測(cè)時(shí)方式獲得的測(cè)時(shí)結(jié)果。在圖4所示的多次試驗(yàn)中，單輻射成分進(jìn)行測(cè)時(shí)給出的相位誤差的標(biāo)準(zhǔn)差為0.00653，而多輻射成分測(cè)時(shí)給出的相位誤差的標(biāo)準(zhǔn)差為0.00271。顯然，當(dāng)脈沖輪廓中存在多種輻射成分時(shí)，使用多成分脈沖標(biāo)準(zhǔn)輪廓進(jìn)行測(cè)時(shí)明顯優(yōu)于使用傳統(tǒng)的單成分標(biāo)準(zhǔn)輪廓進(jìn)行測(cè)時(shí)。

圖3 成分分離獲得的3種不同脈沖輪廓成分（即多成分脈沖標(biāo)準(zhǔn)輪廓）與單成分脈沖標(biāo)準(zhǔn)輪廓對(duì)比圖

圖4 分別利用單成分與多成分脈沖標(biāo)準(zhǔn)輪廓進(jìn)行測(cè)時(shí)的測(cè)時(shí)殘差

3 關(guān)于多成分脈沖輪廓模版獲取方式的討論

在上一節(jié)中采用了數(shù)值模擬的方式展示了多成分脈沖輪廓模版的應(yīng)用，然而模擬數(shù)據(jù)與真實(shí)數(shù)據(jù)存在一定的差別，因此需要針對(duì)性地進(jìn)行相應(yīng)討論。

在以上多成分脈沖輪廓分離過(guò)程的模擬中，由于事先假定了輪廓中存在3個(gè)成分，在擬合時(shí)同樣采用了3成分?jǐn)M合，這樣才得到圖3中所示的結(jié)果。在使用實(shí)際的脈沖星觀測(cè)輪廓進(jìn)行以上操作時(shí)，由于無(wú)法事先判斷脈沖輪廓中存在幾個(gè)輻射成分，需要依次對(duì)成分?jǐn)?shù)進(jìn)行試驗(yàn)，然后考察分離輪廓成分的線性相關(guān)性，然后采用分離出的所有脈沖輪廓成分的最小線性無(wú)關(guān)組作為最終的成分分離結(jié)果。

在脈沖星模式變換時(shí)，有可能會(huì)存在過(guò)渡態(tài)[14]，這一狀態(tài)的脈沖形狀不屬于任何一種穩(wěn)定模式。如果存在過(guò)渡態(tài)，那么將觀測(cè)數(shù)據(jù)直接折疊積分獲得的脈沖輪廓可能并非是簡(jiǎn)單的多種模式的線性組合，而是有可能存在非齊次或非線性效應(yīng)。此時(shí)再使用以上成分分離方法時(shí)可能導(dǎo)致結(jié)果有較大誤差。因此，在實(shí)際的成分分離過(guò)程中，如果采用不同的脈沖輪廓組進(jìn)行成分分離實(shí)驗(yàn)時(shí)獲得了相差較大的結(jié)果，需要參考核主成分分析（kernel PCA，KPCA）[15]方法考慮成分疊加時(shí)的非線性效應(yīng)后再進(jìn)行成分分離。

4 結(jié)語(yǔ)

通過(guò)對(duì)存在多種輻射模式的脈沖星的模擬觀測(cè)輪廓進(jìn)行成分分離實(shí)驗(yàn)，發(fā)現(xiàn)對(duì)該類脈沖星進(jìn)行測(cè)時(shí)時(shí)，使用多成分標(biāo)準(zhǔn)輪廓給出的測(cè)時(shí)結(jié)果的精度要高于單成分標(biāo)準(zhǔn)輪廓的測(cè)時(shí)結(jié)果精度。這意味著在進(jìn)行脈沖星測(cè)時(shí)時(shí)，對(duì)存在多種輻射成分的脈沖星采用多成分脈沖標(biāo)準(zhǔn)輪廓是有必要的，且對(duì)于無(wú)法判斷是否存在多種輻射成分的脈沖星，也有必要采用成分分離法對(duì)其輪廓是否含有多種線性無(wú)關(guān)輪廓成分進(jìn)行判斷，以便于獲得更高精度的測(cè)時(shí)結(jié)果。

[1] HELFAND D J, MANCHESTER R N, TAYLOR J H. Observations of pulsar radio emission. III-stability of integrated profiles[J]. The Astrophysical Journal, 1975(198): 661-670.

[2] TAYLOR J H, MANCHESTER R N, HUGUENIN G R. Observations of pulsar radio emission. I-total intensity measurements of individual pulses[J]. The Astrophysical Journal, 1975(195): 513-528.

[3] BARTEL N, MORRIS D, SIEBER W, et al. The mode-switching phenomenon in pulsars[J]. The Astrophysical Journal, 1982, 258: 776-789.

[4] RANKIN J M. Toward an empirical theory of pulsar emission. III-mode changing, drifting subpulses, and pulse nulling[J]. The Astrophysical Journal, 1986(301): 901-922.

[5] REDMAN S L, WRIGHT G A E, RANKIN J M. Pulsar PSR B2303+30: a single system of drifting subpulses, moding and nulling[J]. Monthly Notices of the Royal Astronomical Society, 2005, 357(3): 859-872.

[6] WANG N, MANCHESTER R N, JOHNSTON S. Pulsar nulling and mode changing[J]. Monthly Notices of the Royal Astronomical Society, 2007, 377(3): 1383-1392.

[7] TIPPING M E, BISHOP C M. Probabilistic principal component analysis[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 1999, 61(3): 611-622.

[8] YIP C W. Spectral classification of quasars in the sloan digital sky survey: eigen spectra, redshift, and luminosity effects[J]. The Astronomical Journal, 2004, 128(6): 2603-2630.

[9] JARVIS M, JAIN B. Principal component analysis of PSF variation in weak lensing surveys[J/OL]. (2004-12-19)[2023-03-10]. https://arxiv.org/abs/astro-ph/0412234.

[10] DELCHAMBRE L. Weighted principal component analysis: a weighted covariance eigen decomposition approach[J]. Monthly Notices of the Royal Astronomical Society, 2015, 446(4): 3545-3555.

[11] HE J J, TAN E L, GAN W S. Time-shifted principal component analysis based cue extraction for stereo audio signals[C] // IEEE International Conference on Acoustics, Speech and Signal Processing, 2013.

[12] CHEVEIGNéA D, SIMON JONATHAN Z. Denoising based on time-shift PCA[J]. Journal of Neuroscience Methods, 2013, 165(2): 297-305.

[13] HYVARINEN A, OJA E. Fast independent component analysis, algorithms and applications[J]. Neural Networks, 2000, 13(4-5): 411-430.

[14] SULEYMANOVA S, BILOUS A. The discovery of transitive phenomenon in the radio emission of the mode-switcher PSR B0943+10[J/OL]. (2022-10-9)[2023-03-10].https://arxiv.org/abs/2210.04215.

[15] SCH?LKOPF B, SMOLA A, MüLLER K R. Kernel principal component analysis[C]//International conference on artificial neural networks, Berlin, Heidelberg: Springer, 1997.

Multi-component pulse profile template and the application in pulsar timing

LU Ji-guang1,2, LIU Yu-lan1,2, JIANG Peng1,2

(1. National Astronomical Observatory, Chinese Academy of Sciences, Beijing 100101, China;2.Guizhou Radio Astronomical Observatory, Guizhou University, Guiyang 550000, China)

The mode change phenomena exist in the radio observation data of some pulsars, and the change of mean pulse profile introduces difficulties into pulsar timing. This paper introduces a scheme for decomposing the observed pulse profiles with multiple components, and shows that multi-component pulse profile template is better than single-component template in measuring the pulse arrival time.

pulsar; pulsar timing; mode switching; pulse profile template; PCA(principal component analysis)

盧吉光, 劉雨蘭, 姜鵬. 多成分脈沖輪廓模版及其在脈沖星測(cè)時(shí)中的應(yīng)用[J]. 時(shí)間頻率學(xué)報(bào), 2023, 46(3): 171-177.

10.13875/j.issn.1674-0637.2023-03-0171-07

2023-03-10；

2023-03-29

國(guó)家自然科學(xué)基金資助項(xiàng)目（12003047）