经验模态分解方法在脉冲星计时残差分析及预报中的应用∗

高 峰高玉平童明雷杨廷高赵成仕

(1中国科学院国家授时中心西安710600)

(2中国科学院时间频率基准重点实验室西安710600)

(3中国科学院大学北京100049)

(4西安科技大学应用物理系西安710054)

1 引言

脉冲星计时分析是整个脉冲星领域研究的基石,在诸如强引力场中的广义相对论效应、检验引力波、探测太阳系外行星以及提供计时标准等研究中[1−4],都需要高精度的脉冲星计时分析的结果.当前,由于计时模型不完善及一些未能消除的误差源的存在,限制了脉冲星计时分析精度进一步提高,其影响主要体现在脉冲星计时残差的波动性上.脉冲星计时残差是脉冲星时与原子时比对的产物,其波动性主要来源于脉冲星计时观测误差、脉冲到达时间(Time of Arrival,TOA)转换误差、脉冲星本身内部噪声等不确定因素的影响.不同的误差源对计时残差的影响各不相同[5−7],分析计时残差诸多波动特征,将有利于解读上述各误差的来源.参考文献[8−10]已报道了利用指数模型、频域方法、功率谱方法等实现计时残差中低频噪声的度量及分类.频域分析的最大优点是通过傅立叶或小波变换能将波动信号的某些特征清楚地显现出来.然而,上述变换所用的基底函数系和“窗口”分析都具有明显的主观选择性.Huang等人于1998年提出了经验模态分解(Empirical Mode Decomposition,EMD)方法用于分析非线性非平稳时间序列[11].EMD方法具有良好的功能:对数据具有自适应性,能用几个内在的本征模态函数(Intrinsic Mode Function,IMF)和一个剩余趋势揭示时间序列的振荡结构和非平稳性.EMD方法一经提出,便在振动故障诊断、光学信号、滤波及噪声消除等多个领域得到了广泛应用.

本文基于EMD方法分析时间序列的良好功能,将脉冲星计时残差看作时间序列并采用EMD方法进行分析,旨在揭示计时残差内在的振荡结构,这有助于理解不同误差源对脉冲星计时的影响;其次,脉冲星计时残差的预报能反推各类误差源影响的趋势,本文进一步利用EMD方法分解后的IMF分量或其组合,探讨EMD方法在脉冲星计时残差预报方面所具有的优势.

2 EMD方法的数学原理

EMD方法是分解信号序列的一种有效工具,它把时间序列分解为有限个IMF分量和一个剩余趋势项[12].各分量应满足如下约束条件:(1)在整个资料集里,极值点的数目与穿零点的数目必须相等或者最多相差1个;(2)由局部极大值所构成的包络线以及由局部极小值所构成的包络线的平均值为零.

EMD方法的关键技术问题是如何把一个非线性非平稳时间序列(假定初始时间序列为X(t))分解为有限个IMF分量和一个趋势项.解决问题的方法是通过筛选,即:首先,通过3次样条曲线联结所在局部极大值得到上包络线,用同样的方法连接局部极小值形成下包络线,上下两条包络线应当分别涵盖全部的极大值和极小值;其次,取上下包络线的均值线为m1(t).原序列X(t)与均值m1(t)之差为h1(t),即

这样得到的h1(t)未必满足上述IMF分量的约束条件.把h1(t)当作初始时间序列,重复筛选过程

直到第k次筛选的结果符合要求,这样就得到了第1个IMF分量C1.C1应是信号的特征时间尺度最小的高频IMF分量.接下来,把C1从X(t)中分离出去,得到剩余时间序列.剩余时间序列重复(1)–(2)式的筛选过程,直到剩余项成为一个单调函数.于是得到一系列的IMF分量,即

(3)式中,C1,···,Cn是IMF的n个分量,r是剩余趋势项.

3 数据来源及分析

本文选取4颗毫秒脉冲星B1937+21、J0613–0200、J1012+5307、J1643–1224进行分析,其相关TOA数据和脉冲星参数来源于Arzoumanian等[13]2015年公布的北美计时阵(North American Nanohertz Observatory for Gravitational Waves,NANOGrav)观测的37颗毫秒脉冲星数据.上述数据经标准脉冲星数据处理流程拟合脉冲星天文参数,随后获得的计时残差数据作为本文分析的原始数据.4颗毫秒脉冲星观测跨度最长为9.2 yr,释放的TOA数据计数多达11995点(部分基本参数见表1,P为脉冲星自转周期).

表1 4颗脉冲星基本参数Table 1 The basic parameters of 4 pulsars

本文选取脉冲星数据的标准是基于以下考虑:首先,脉冲星计时模型参数精度的提高虽然有赖于多种因素,但观测TOA计数越多,越有利于脉冲星计时模型参数的拟合,从而降低拟合误差的引入;其次,采样时间跨度较长,将有利于脉冲星计时残差中红噪声成分的测量.所选的4颗脉冲星在文献[13]中皆给出了残差的红噪声模型,在本文分析中,其约9 yr观测跨度将非常有利于解析低频噪声对计时残差的影响.目前,脉冲星计时观测并不是规则采样,为了避免或减小因数据不规则采样引入新的统计误差,本文所选脉冲星计时残差数据均做了规则化预处理.其处理过程以脉冲星J0613–0200为例进行说明:根据表1第3列和第5列可知,脉冲星J0613–0200在约8.6 yr的观测跨度内有7651个计数点,对应的约化儒略日(Modified Julian Day,MJD)有较高的分辨率,即同一天可能对应数个观测点.令同一天有观测点时取样1个,其值为同一天内数个点求平均,对应MJD数值向下取整.这样处理过程等同于1 d之内脉冲星时与原子时比对1次取计时残差,经过处理后脉冲星J0613–0200计时残差为171个点.随后进行线性插值取样,获得均匀分布的同样数目计时残差点数,即171个点.对其他3颗脉冲星残差做相同规则化预处理,预处理后的4颗脉冲星计时残差作为本文框架内分析的初始数据,其残差分布如图1所示.

将脉冲星计时残差数据作为信号序列应用EMD方法进行分解,将讨论以下问题:

(1)各脉冲星计时残差的波动由哪些时间尺度的振荡模态变化构成;

(2)哪些振荡模态对整个计时残差的变化起主要作用,不同脉冲星计时残差波动有何共性及个性;

(3)不同振荡模态时间尺度的变化与脉冲星计时测量误差来源有何联系.

图2(a)–(h)分别表示脉冲星计时残差数据X(t)、应用EMD方法将X(t)进行分解后得到的6个IMF分量C1–C6和剩余趋势项r的分布.实线、虚线、点线、点虚线分别表示数据来源于脉冲星B1937+21、J0613–0200、J1012+5307、J1643–1224的残差.整体来看,对比X(t)可以发现:图2(a)中4颗脉冲星残差数据X(t)幅值大小涨落很不规则且有明显的低频特征.图2(b)–(g)中分解之后C1–C6分量呈现出围绕零均值线振荡且局部极大值和极小值基本对称的振荡形式.它们的均值都为零,不随时间变化,频率变化较小,波形具有比X(t)规整、简单、非平稳性减弱等特点.图2(h)中剩余趋势项r是近似单调变化的曲线.除此之外,分解之后通过不同IMF分量可以清晰地表明X(t)的多尺度变化,C1–C6分量时间尺度逐渐增大,亦即从高频到低频.上述特点都非常符合EMD方法分解时序信号的特征.

图1 4颗脉冲星计时残差分布.纵坐标是计时残差,横坐标是MJD的跨度.实线、虚线、点线、点虚线分别代表脉冲星B1937+21、J0613–0200、J1012+5307、J1643–1224的残差数据.Fig.1 The timing residuals of 4 pulsars.The vertical axis stands for the timing residual,and the horizontal axis presents span of MJD.The solid,dashed,dotted,and dot-dashed lines represent B1937+21,J0613–0200,J1012+5307,and J1643–1224,respectively.

考察IMF分量对于残差数据X(t)的相对重要性时,可根据各分量振幅变化的量级(如图2的纵坐标数值)比较;也可以参照谐波分析,用方差贡献表现各分量对原序列的重要性.考虑到标准差与数据的量纲一致,在描述数据波动性时更明显,且方差与标准差很容易转换,本文采用标准差贡献来表征各分量对原序列的重要性.表2第2列是各脉冲星初始残差序列的标准差值,其中脉冲星J1643–1224标准差最大,为2.5658µs. 脉冲星J0613–0200标准差值最小,为0.7575 µs. 第3–9列是IMF各分量C1–C6及剩余趋势项r对各脉冲星残差序列波动性的标准差贡献.相对于后3颗脉冲星各分量贡献而言,脉冲星B1937+21各分量中标准差最大的3个依次为C6,r和C4,它们的贡献分别为0.7482µs,0.5321µs和0.4748µs. 尤其是C6,比C1–C3大了数倍,说明C6的贡献对X(t)振荡起决定性作用,即X(t)中最显著的低频特征由C6反映出来.脉冲星J0613–0200除C3–C4分量偏小外,其他分量及趋势项标准差值相当.脉冲星J1012+5307中C2最大,为1.0623µs,除C6外,其他各项随时间尺度的增大,标准差贡献都是下降的.而脉冲星J1643–1224各IMF分量和趋势项标准差值相差不大.

图2 图(a)–(h)分别是4颗脉冲星计时残差X(t)、IMF分量C1–C6和剩余趋势项r的分布.实线、虚线、点线、点虚线分别表示数据来源于脉冲星B1937+21、J0613–0200、J1012+5307、J1643–1224的残差.Fig.2 The panels(a)–(h)indicate pulsar timing residuals X(t),C1–C6of the intrinsic mode function,and trend term r,respectively.The solid,dashed,dotted,and dot-dashed lines represent B1937+21,J0613–0200,J1012+5307,and J1643–1224,respectively.

表2 各IMF分量对4颗脉冲星计时残差波动性的标准差贡献Table 2 The contributions of IMF components to the standard deviations of the volatility for 4 pulsars timing residuals

综上所述,我们可以看到:采用EMD方法分析能较好地用几个内在IMF分量和一个剩余趋势揭示各脉冲星残差序列的振荡结构和非平稳趋势.尤其是结合各分量的标准差贡献可以清楚地看出某个或几个IMF分量在原残差序列中所起的作用.本文所选的4颗脉冲星其共性都具有明显的低频波动(红噪声),这与文献[13]能用红噪声模型描述残差中的红噪声成分是一致的.但通过EMD方法分解后分析可以看出,4颗脉冲星残差的波动特征仍有较大的差异.以脉冲星B1937+21为例,其低频成份贡献非常明显,这表明其残差序列的频谱分析在低频端必定有明显的能量包存在.而其他3颗脉冲星残差低频成分则没有这么显著.另外,4颗脉冲星的差异也表现在各自的剩余趋势项r上,结合图2(h)可以看出:脉冲星B1937+21和J1643–1224近10 yr来残差序列总体上是上升趋势,而脉冲星J0613–0200和J1012+5307是下降趋势.最后,如果简单地将脉冲星计时残差区分为白噪声和红噪声两种成分,白噪声主要来源于计时观测的随机误差,具有高频、随机性强等特征.而红噪声来源广泛,脉冲星本身内部噪声、传播误差、脉冲星到达时间转换误差等都可累积形成低频噪声.单个或多个因素组合形成周期性特征并且对残差波动贡献较强时,理论上EMD方法都能在时域分析中做到较好的辨识.所以,EMD方法是研究不同误差来源对脉冲星计时影响的一种重要的工具,可与计时残差的频谱分析形成交叉证认.

4 ARMA模型及其实现

自回归滑动平均模型(Auto-Regressive Moving Average Model,ARMA)是时间序列分析与前推预测领域中一类常用的重要方法[14].模型形式ARMA(p,q)中包含了p阶自回归项和q阶移动平均项.ARMA(p,q)模型可以表示为

其中,Xt是平稳或准平稳的时间序列,φi和θj分别是自回归参数和滑动平均参数,εt为均值为零、方差为σ2的正态白噪声,可认为是随机干扰项.

(4)式中,若Xt仅与t时刻的白噪声εt和以前的变化Xt−i有关,即

此时,(5)式称为p阶自回归(Auto-Regressive,AR)模型,可记作AR(p)或者ARMA(p,0),其中p为模型阶数.

若Xt仅与t时刻的白噪声εt及以前的白噪声εt−j存在关联,即

则(6)式称为q阶滑动平均(Moving Average,MA)模型,记作MA(q)或是ARMA(0,q).所以,ARMA(p,q)模型可看作是AR模型与MA模型的有机综合.

利用ARMA模型进行时间序列的预报通常由以下4步确认或计算:

(1)模型识别 对于一个平稳时间序列,可以通过该序列的自相关函数(Auto-Correlation Function,ACF)和偏相关函数(Partial Correlation Function,PCF)的“拖尾”(Trailing Lag,相关函数随延迟量增加逐渐归于零)或“截尾”(Truncation Lag,延迟量大于某值后相关函数恒等于零)特性来合理判断序列更符合哪个模型,其判断准则可见表3.

(2)模型定阶模型类别确立后,模型阶次可由信息论准则、传递函数准则或是最终预测误差准则等多种方式确定模型的最优阶次[15].本文采用Akaike信息论法则(Akaike Information Criterion,AIC)来选择模型的最佳阶次.

(3)模型参数估计模型参数估计方法有多种,本文采用Yule-Walker方程估计法.

(4)预报模式选择预报模式同样多样且各具优点,本文采用递推模式进行预报,其特点是只需一次建模计算、计算量小、可减少建模的复杂性可能引入的多种误差.

表3 ARMA模型识别准则Table 3 ARMA identification criterion

5 脉冲星计时残差的预报

在第4节中,我们讨论了ARMA模型相关原理及其对平稳时间序列进行预报的步骤.本节中将进一步讨论应用EMD方法对脉冲星计时残差预报的优势.需要说明的是,时序数据的回归预报对数据的平稳性有较高的要求,脉冲星J0613–0200残差数据相较于其他3颗脉冲星有更好的准平稳性.所以,本节中仅对脉冲星J0613–0200应用ARMA模型进行计时残差序列的预报,各参与预报的数据均通过了用于单位根检验的扩展迪基-福勒检验(Augmented Dickey-Fuller test)[16−17].

开展计时残差序列预报的思路是:首先截取初始数据前部采用直接建模预报和对截取数据应用EMD方法分解后建模间接预报两种模式;随后两种预报结果与初始数据后部(真实观测数据)进行比对,以检验两种预报方式的优劣及其精度.其具体步骤如下:

(1)直接预报 取脉冲星J0613-0200残差数据(共171个点)的前161个点建模并预报后10个点;

(2)间接预报 取残差数据的前161个点并应用EMD方法分解,利用分解后的IMF分量或其组合建模并预报后10个点.

采用直接预报模式时,根据表3中ARMA模型识别准则,发现建模数据适合AR(p)模型.随后逐次进行模型定阶、参数估计,并采用递推模式进行预报.采用EMD方法分解进行间接预报有两点考虑:一是在比对预报结果时为体现两种预报模式的优劣,本文中将控制间接预报模型、参数估计和预报模式与直接预报流程中选择尽量相同或类似;二是EMD方法分解后的IMF高频分量会出现随机性太强不能提取有益信息进行预报,或是低频分量趋势性太强不符合平稳时间序列的要求.本文采用各分量组合重构的方式解决上述问题.如图3所示,上图为C1与C3相加重构,中图为C2加C5重构,下图为C4与趋势项r重构.这类组合方式并没有检验是否最优,但各重构时序均通过了单位根检验,并根据表3中ARMA模型识别准则,均适合AR(p)模型.类似于直接预报流程,3个重构时序分别给出各自的预报结果,最后叠加给出间接预报的结果.

图4中给出了直接预报数据、间接预报数据与初始数据(后10个点,下同)比对的结果,以此检验两种预报模式的优劣及其精度.4个标注参量X(t)、AR(X)、AR(C)和X∗(t)分别表示初始数据、直接预报数据、间接预报数据和初始数据平滑后的值.

图3 不同IMF分量的组合重构.上图是C1和C3累加,中图是C2和C5累加,下图是C4和r累加.Fig.3 The reconstruction of different IMF components.The top panel stands for the sum of C1and C3,the middle panel is the sum of C2and C5,and the bottom panel presents the sum of C4and r.

图4 脉冲星计时残差的预报数据和初始数据的比对.虚线和点线分别表示直接预报数据和间接预报数据.上下图实线分别表示初始数据和平滑后的初始数据.Fig.4 The comparison between the forecast data and the initial data of pulsar timing residuals.The dashed and dotted lines indicate the direct and indirect forecast data,respectively.The solid lines in the upper and lower panels represent the initial data and the smoothed initial data,respectively.

分析图4可发现:上图显示初始数据X(t)随机起伏较大,直接预报数据AR(X)和间接预报数据AR(C)较为平滑.检验预报结果好坏最直接的方法是比对预报数据与初始数据接近的程度.为消除初始数据随机性因素的干扰,下图将初始数据进行平滑后与预报数据进行比对.图5中显示的是两种预报数据AR(X)、AR(C)分别与平滑数据X∗(t)作差后的结果.从接近程度上看,AR(C)在短期(约4个月)范围内要明显优于AR(X).图5中更长期(大于4个月)预报时,直接预报数据AR(X)更接近平滑数据X∗(t),但考虑初始数据末尾存在端部效应等因素,上述结果不能形成明显结论.除此之外,初始数据的1倍标准差约为0.37µs,短期预报时,间接预报数据明显落在初始数据1倍标准差之内,具有较高的精度.所以,综上所述,采用EMD方法分解后对脉冲星J0613–0200计时残差进行短期预报,其预报模式合理有效,且预报结果明显优于直接预报.

图5 脉冲星计时残差的预报数据和平滑后初始数据的比对.实线和短划线分别来自直接预报数据和间接预报数据与平滑后初始数据作差的绝对值.Fig.5 The comparison between the forecast data and smoothed initial data of pulsar timing residuals.The solid line indicates the absolute value of the difference between direct forecast data and the smoothed initial data,and the dashed line represents the absolute value of the difference between the indirect forecast data and the smoothed initial data.

6 结论

分析脉冲星计时残差的波动特征有利于研究各类误差源对脉冲星计时精度的影响.EMD方法分析具有对时序数据较强的自适应性,能用几个内在的本征模态函数和一个剩余趋势揭示时间序列的振荡结构和非平稳性等特点,非常适合脉冲星计时残差数据的分析.本文采用EMD方法对脉冲星B1937+21、J0613–0200、J1012+5307、J1643–1224的计时残差振荡模式进行了分析,并进一步探究了应用EMD方法在计时残差预报方面所具有的优势.通过以上研究,我们发现:

(1)脉冲星的计时残差的波动可以分解为多个振荡模态和一个趋势项,不同的振荡模态清晰地表明了计时残差波动的多时间尺度的变化.通过比较4颗脉冲星各振荡模态的方差贡献,发现脉冲星B1937+21低频成分对残差波动贡献最大,其余3颗脉冲星则不是很显著.这种脉冲星残差序列波动的差异性也表现在分解后的剩余趋势项上.

(2)EMD方法具有在时域上辨识残差序列高频成分与低频成分的能力,在研究不同误差来源对脉冲星计时残差影响领域,可与残差的频谱分析形成交叉证认.

(3)通过分析单颗脉冲星J0613–0200计时残差序列回归预报,揭示在短期预报范围内采用EMD方法将计时残差序列分解后再预报要优于直接预报结果.究其原因,可认为EMD方法是将时间序列平稳化的一种过程,这将有利于提高模型预报精度.

致谢此研究利用了北美计时阵的观测数据,并特别感谢中国科学院国家授时中心雷雨博士对自回归模型的有益讨论.