基于EMD-AR谱的3C 273光变周期分析*

蔡剑华 刘备 杨江河， 庹满先 聂建军 汪胜辉 樊军辉

(1 湖南文理学院数理学院 常德 415000)

(2 广州大学天体物理中心 广州 510006)

1 引言

耀变体是活动星系核中观测性质极其特殊的子类， 而光变是它们最突出的观测特征， 大多数表现为长期的无规则快速光变， 有的光变表现出周期性[1-7]. Blazar的光变表现出光谱随亮度的变化， 光变时标不但带来了辐射区大小信息， 甚至带来了Blazar喷流强弱的信息[8-16]. 周期性光变可能是Blazar内部结构的反映[8]， 但是由于观测数据分布的不均匀性， 使得从光变曲线中分析周期性会有一定的不确定性. 在长期的研究中， 人们寻找天体光变周期的常用方法有:周期图谱法[17-18]、相位分析法[19]、相关函数方法[20-21]、Jurkevich方法[22]、模型谱估计法[23-25]、小波分析法[26-28]和Hilbert-Huang变换[29]等. 后来， 为了更好地处理天体观测中的非等间隔数据， 人们又提出了CLEAN方法[30]、CLEANEST方法[31-32]和基于小波的向量投影方法[32]. 这些方法在天文观测和理论研究中得到了广泛应用， 并有作者不断将最新发展起来的现代信号处理方法应用到天体物理领域. 3C 273是最早发现的类星体， 是研究最多的天体之一. 对类星体3C 273光变周期的分析， 已有学者开展了相关的工作. Smith等[33]研究表明3C 273在1887-1963年的光变曲线中有5-6个不同周期的波动， 并称其有12.7-15.2 yr的周期(P). Babadzhanyants等[34]分析了3C 273在1887-1991年的B波段光变曲线， 发现P= 13.4 yr; 2001年樊军辉等[35-36]研究了3C 273的110 yr光变数据，发现了其P=2.0、(13.65±0.20)和(22.5±0.2)yr的周期.2007年Fan等[37]还研究了它的射电光变周期性， 发现其14.5 GHz光变曲线的周期为(8.2±0.2)yr，8.0 GHz光变曲线的周期为(8.3±0.2) yr和(19.3±1.7) yr， 4.8 GHz光变曲线的周期为(8.8±0.3) yr. 2013年Vol’vach等[38]根据3C 273在1963-2011年的光变曲线， 得出了P=(11.2±2.3)、(7.2±0.8)、(4.9±0.3)和(2.8±0.3) yr的周期. 2014年Fan等[39]综合自己以前的工作和对1998-2008年观测数据的分析， 进一步得出3C 273的光变周期为P=(21.10±0.14)、(10.90±0.14)、(13.20±0.09)、(7.30±0.10)、(2.10±0.06)和(0.68±0.05)yr.2014年唐洁[40]基于密歇根大学射电天文台数据库(Michigan Radio Astronomy Observatory， UMRAO)从1965年到2012年所收集的3C 273天体47 yr的8.0 GHz射电波段观测数据， 得到3C 273的平均周期是16.9、7.5和2.5 yr.

本文试图在前人工作的基础上， 结合近期发展起来的非平稳信号处理方法-经验模态分解和传统自回归(Auto-regressive，AR)模型谱估计的优点，提出一种组合的功率谱估计方法， 并将他们应用于类星体3C 273的光变周期估算， 旨在为天体的周期光变研究探索一条新途径.

2 方法原理

2.1 经验模态分解

经验模态分解(Empirical Mode Decomposition， EMD)， 又叫Huang变换， 它是由美籍华人科学家N. E. Huang于1998年提出的一种新型信号分解方法. 与被广泛使用的小波分解不同， 经验模态分解无须选择基函数， 特别适合于对非平稳、非线性信号的处理. 经验模态分解采用“筛”选的方法， 将信号x(t)分解为一组固有模态函数(Intrinsic Mode Function， IMF)的和[41-42]，t表示时间.

其中n表示分解得到的最大阶数，i表示取其中的各阶，Ci(t)即为各阶模态函数(为作图和描述方便常简为Ci)，R(t)为剩余分量， 具体分解过程见文献[41]. 固有模态函数满足信号关于时间轴局部对称和零点数与极点数相等或至多相差1的两个条件， 包含了原信号不同时间尺度的局部特征信息.EMD分解实质上是对非平稳数据进行平稳化处理的一个过程， 在空间域表现为从小尺度到大尺度的层层滤波， 在频域则表现为从高频到低频的层层滤波. 同时， Huang等[41-42]也证明了经验模态分解的完备性， 即分解得到的各固有模态以及剩余分量相加能够恢复原始信号. 这些特征为信号的相关性分析和构造滤波器抑制噪声提供了条件[43-44].

若要消除高频部分的噪声， 一个包含多个IMF的低通时空滤波器可构造为

k为要保留的IMF的最小阶数. 若噪声仅存在低频部分， 高通时空滤波器为

h为要保留的IMF的最大阶数. 若噪声在高低频都存在， 则带通时空滤波器为

本文就是在EMD分解后得到各阶分量， 并计算其与原始光变数据的相关度， 以此构造滤波器，取相关系数大的系列IMF分量来进行进一步谱估计的.

2.2 AR谱估计方法

随机过程s(m)的AR模型可以表示为[45-46]

m为离散数据点，al(l=1，2，···，p)为AR模型的系数，l为各系数序列标号，p为系数序列标号的最大值.u(m)是一个方差为δ2的白噪声序列. 求解AR模型系数常用的方法有Yule-Walker法、Burg法和最小二乘法， 本文采用的是Yule-Walker法. AR模型常用的阶数判定准则有3种: FPE (Final Predietion Error)准则、AIC (Akaike Information Criterion)准则和BIC (Bayesian Information Criterion)准则， 本文采用的是AIC准则. 随机过程的功率谱估计定义为[47-48]

j为虚数单位，ω为角变量. 以上(5)式和(6)式给出的模型被称为自回归模型， 简称AR模型， 获得的功率谱估计称为AR谱. 求得s(m)的前p+1个自相关函数， 代入公式即可求得s(m)的功率谱.

2.3 方法步骤

基于EMD-AR谱的天体光变周期分析流程如图1所示， 描述如下:

图1 基于EMD-AR谱的天体光变周期分析流程Fig.1 Analysis flow of variability periodicity based on EMD-AR spectrum

(1)首先对观测数据进行EMD分解， 得到各阶模态函数“C1”-“Cn”;

(2)计算各阶IMF分量与原始观测数据的相关度， 去除相关系数小的分量， 用余下的IMF分量“Ch”-“Ck”重构数据，k和h的取值由(2)-(4)式确定;

(3)最后估算其AR谱， 进一步分析天体的光变周期.

3 数据样本

本文分析的数据来自广州大学樊军辉教授课题组在其系列论文[16， 18， 35-37， 39， 49]中收集和观测的类星体3C 273光变数据， 时间跨度为1887-2016年， 长达130 yr， 共有4342个观测数据. 由于历史原因和天体观测的特点， 有些月份甚至年份的观测数据缺失. 但在后期， 尤其是在2000年后， 数据观测记录的频次很高， 出现了多段连续几个小时观测， 且以1 min为时间间隔记录的观测数据， 这些数据为3C 273短时标周期的分析提供了条件.

用于长时标周期分析的数据如下. 在已有的1887-2016年的共4342个观测数据中， 根据时标每年随机提取出一个数据， 共112个观测数据， 为便于分析我们将此设为L1光变曲线. 在4342个观测数据中， 按时标每月随机提取出一个数据， 共469个观测数据， 设为L2光变曲线. 同上， 在观测数据中，按时标随机每天提取出一个数据， 共1254个观测数据， 设为L3光变曲线. L1-L3光变曲线的时域波形如图2 (a)-(c)所示.

图2 类星体3C 273的观测数据: (a)-(c)用于长时标周期计算; (d)-(f)用于短时标周期计算.Fig.2 Observational data of quasar 3C 273: (a)-(c) for long time scale period calculation， (d)-(f) for short time scale period calculation.

用于短时标周期分析的数据如下. 分析4342个观测数据， 发现有很多段较连续的以“min”为时间间隔的观测数据， 为我们进行天内光变的研究提供了条件. 选取较长的3段: Julian day (JD)2454459.4701-JD2454459.6516数据段中， 在连续的4 h内有192个观测数据， 我们将此设为G1光变曲线; JD2454460.4704-JD2454460.6459数据段中， 在连续的4 h内有198个观测数据， 设为G2光变曲线; JD2454461.4713-JD2454461.6165数据段中，在连续的3 h内有165个观测数据， 设为G3光变曲线. G1-G3光变曲线的时域波形如图2 (d)-(f)所示. 这3段观测数据连续性很好， 几乎是等间隔的以“min”为单位的采样观测， 分别对它们进行周期分析， 尝试得到以“min”为单位的3C 273短时标周期.

4 仿真分析

为验证采用本文方法进行功率谱估计的有效性， 作者首先对一组仿真信号进行了处理. 设计两个频率为30 Hz和70 Hz正弦分量的叠加信号并对其采样， 采样频率为1000 Hz、采样点数为1459以构成仿真信号. 采用本文提出的EMD-AR谱估计方法， 对其进行功率谱的估计.

图3给出了仿真信号EMD分解后的各IMF分量及其对应的功率谱. 从图3可以清楚地看到EMD分解在空间域表现为从小尺度到大尺度， 在频域则表现为从高频到低频的分解过程. 分析图3 (a)“C1”分量和“C2”分量的时域波形图及其图3 (b)对应的频谱， 30 Hz和70 Hz的2个频率成分已被精确地分解出来. 计算“C1”-“R”分量与原始信号的相关度分别为: 0.7548、0.7075、0.0006、0.0025、0.0023、0.0016、0.0011. 显然， “C1”-“C2”分量与原始数据的相关系数比“C3”-“R”分量的相关系数大得多， 他们不在同一数量级. 选择相关系数大的“C1”-“C2”分量重构数据， 并进行AR谱估计. 图4所示为采用本文提出的EMD-AR谱方法得到的功率谱估计. 可以看出， 功率谱的精度很高，30 Hz和70 Hz的频率成分已被精确地分离和表征出来.

图3 仿真信号EMD分解后的各IMF分量及其对应的功率谱Fig.3 IMF components and its power spectrum of simulation signal after EMD decomposition

图4 EMD-AR谱方法估算的仿真信号的功率谱Fig.4 Power spectrum of simulation signal estimated by EMD-AR method

为进一步验证本文方法可有效过滤非周期性信号且不会得到虚假的周期信号， 作者在仿真上加入了一些强度随机、宽度随机的“Bumps”脉冲噪声信号， 再进一步作谱估计. 如图5 (a)所示， “星划”线为原始的两个正弦函数叠加的仿真信号， “点划”线为加入的“Bumps”脉冲噪声. 加噪后的仿真信号如图5 (b)所示. 依照本文提出的方法， 先对信号进行EMD分解，得到的各级IMF分量如图6(a)所示. 分别计算他们与原始仿真信号的相关系数为:0.7255、0.6946、0.0006、0.0006、0.0001、0.0119、0.0039. 取相关系数大的“C1”-“C2”分量重构数据，并进一步作AR谱估计， 得到的功率谱如图6 (b)所示. 显然， 在图6 (b)中30 Hz和70 Hz的频率成分也非常清晰， 本文方法能正确地将周期信息从含噪信号中辨识、提取出来， 可有效过滤这些非周期性信号且不会得到虚假的周期信号， 为正确估算光变周期提供了条件.

图5 加噪前后的仿真信号: (a)原始仿真信号和“Bumps”脉冲噪声; (b)加入“Bumps”脉冲噪声后的仿真信号.Fig.5 Simulation signals before and after adding noise: (a) original simulation signal and “Bumps” impulse noise; (b)simulation signal after adding “Bumps” impulse noise.

图6 加噪仿真信号的EMD分解(a)和本文方法估算的加噪仿真信号的功率谱(b)Fig.6 EMD decomposition of noisy simulation signal (a) and power spectrum of noisy simulation signal estimated by this method (b)

5 光变周期分析

5.1 长时标周期计算

L1光变曲线是以年为单位的数据样本， 在1887-2016年的观测数据中， 每年随机保留1个数据. 对提取出的数据进行三次样条插值预处理(本文三次样条插值采用的边界条件是非结点边界条件， 由于本数据曲线两端数据差的绝对值也较小， 非结点边界条件有利于对周期的估算， 下同.)， 然后对其进行EMD分解， 自适应地得到6阶模态函数“C1”-“R”， 如图7 (a)所示. 计算他们与原始信号的相关系数分别为: 0.5455、0.5611、0.3500、0.2371、0.0232、0.2338. 去除相关度很小的“C5”分量，用其余分量重构信号， 进一步估算AR谱如图7 (b)所示.从图7的功率谱中， 可以计算出L1光变曲线的光变周期为11.02和4.69 yr.

图7 L1光变曲线的EMD分解(a)及其光变周期估计(b)Fig.7 EMD decomposition of L1 light curve (a) and estimation of variability periodicity (b)

类比L1的处理， 对L2光变曲线进行EMD分解后得到8阶模态函数“C1”-“R”， 如图8 (a)所示. 分别计算他们与原始信号的相关系数为: 0.4979、0.3705、0.4546、0.4411、0.2935、0.0509、0.1208、0.1593. 去除相关度很小的“C6”分量， 用其余分量重构信号， 再进一步估算AR谱如图8 (b)所示. 从图8的功率谱中， 计算得到L2光变曲线的光变周期为13.51、2.76和1.46 yr.

图8 L2光变曲线的EMD分解(a)及其光变周期估计(b)Fig.8 EMD decomposition of L2 light curve (a) and estimation of variability periodicity (b)

同上处理方法， 对L3光变曲线进行EMD分解后得到11阶模态函数“C1”-“R”， 如图9 (a)所示. 计算他们与原始信号的相关系数为: 0.0404、0.0493、0.0822、0.1269、0.1341、0.3268、0.4112、0.4100、0.2474、0.0939、0.4152. 去除相关度很小的“C1”、“C2”、“C3”和“C10”分量， 用其余分量重构信号，再进一步估算AR谱如图9 (b)所示. 从图9的功率谱中， 可以计算出L3光变曲线的光变周期为: 21.23、11.02、5.51、4.69、3.79、2.04 yr.

图9 L3光变曲线的EMD分解(a)及其光变周期估计(b)Fig.9 EMD decomposition of L3 light curve (a) and estimation of variability periodicity (b)

5.2 短时标周期计算

G1段数据为4 h内连续记录的192个观测数据，采用三次样条插值的方法补充完整数据为240采样点， 然后对其进行EMD分解. 如图10 (a)所示， 光变曲线被自适应地分解， 得到7阶模态函数“C1”-“R”，计算他们与原始信号的相关系数分别为: 0.5864、0.5097、0.4326、0.2836、0.3491、0.0515、0.3796.去除相关度很小的“C6”分量， 用其余分量重构信号， 再进一步估算AR谱如图10 (b)所示. 从图10的功率谱中， 可以计算出G1光变曲线的周期为: 240、30、15.3、10、7.3、6.4、4.2 min.

图10 G1光变曲线的EMD分解(a)及其光变周期估计(b)Fig.10 EMD decomposition of G1 light curve (a) and estimation of variability periodicity (b)

类比G1段数据的处理， 对G2光变曲线进行EMD分解后得到8阶模态函数“C1”-“R”，如图11(a)所示. 分别计算他们与原始信号的相关系数为:0.7559、0.4842、0.2828， 0.2951、0.1245、0.0329、0.1922、0.1800. 去除相关度很小的“C6”分量，用其余分量重构信号， 再进一步估算AR谱如图11 (b)所示. 从图11的功率谱中， 可以计算出G2光变曲线的周期为: 240、30、14.9、7.5、6.2、4.4、3.5、3.02 min.

图11 G2光变曲线的EMD分解(a)及其光变周期估计(b)Fig.11 EMD decomposition of G2 light curve (a) and estimation of variability periodicity (b)

同上处理方法， 将G3段数据进行EMD分解后得到7阶模态函数“C1”-“R”， 如图12 (a)所示. 计算他们与原始信号的相关系数分别为: 0.6786、0.3648、0.3090、0.3922、0.1689、0.2518、0.1924.7个相关度系数都在同一个数量级， 所以保留所有分量重构信号， 再进一步估算AR谱如图12 (b)所示. 从图12的功率谱中， 可计算出G3光变曲线的周期为: 180、29.1、15.1、10、7.2、5、4.4、3.5、3.03 min.

图12 G3光变曲线的EMD分解(a)及其光变周期估计(b)Fig.12 EMD decomposition of G3 light curve (a) and estimation of variability periodicity (b)

5.3 结果与讨论

为便于分析与讨论， 将3条长时标周期光变曲线估算的光变周期列表， 并与一些文献中得到的类星体3C 273的光变周期作对比， 如表1所示(P1-P6表示6个估算的周期). 本文得到的21.23 yr的周期与Fan等[39]中给出(21.10±0.14) yr的结论一致. 本文得到的13.51 yr的周期与Smith等[33]给出的12.7-15.2 yr、Babadzhanyants等[34]给出的13.4 yr和Fan等[39]给出的(13.20±0.09) yr吻合得很好. 本文得到的11.02 yr的周期与Vol’vach等[38]给出的(11.2±2.3）yr和Fan等[39]给出的(10.90±0.14) yr一致性很好， 且本文得到了5.51 yr的周期， 与11.02 yr和21.23 yr的周期呈较好的倍周期关系. 本文得到的4.69 yr、2.76 yr的周期与Vol’vach等[38]给出的(4.9±0.3) yr和(2.8±0.3)yr的周期一致. 综合上述分析结果， 本文得到的类星体3C 273的长时标周期为: 21.23、11.02、5.51、13.51、4.69、3.79、2.76 yr.值得一提的是，因为原始数据跨距了130 yr， 由于历史原因， 早期数据的可靠性无法保证， 我们只用近半个世纪的数据重复做了分析， 得到的结果是一致的.

表1 不同方法估算的类星体3C 273的长时标周期Table 1 Long time scale period of quasar 3C 273 estimated by different methods

同样的， 为便于分析， 将3段短时标周期数据估算的光变周期列表， 如表2所示(P1-P9表示9个估算的周期). G1和G2光变曲线计算所得的240 min周期及G3光变曲线所得的180 min周期， 分别与他们的数据等长， 不能作为周期判定. 3条光变曲线分别在30、15和7.5 min， 10和5 min以及6和3 min附近都出现了峰值， 且表现为倍周期关系. 综合3段光变曲线的分析， 可预测类星体3C 273的短时标周期为: (30±1)、(15±0.3)、(7.5±0.2)、(10±0.1)、(5±0.6)和(6±0.4)、(3±0.5) min.

表2 不同光变曲线估算的类星体3C 273的短时标周期Table 2 Short time period of quasar 3C 273 estimated by different light curves

6 结论

本文利用类星体3C 273长达130 yr的观测数据， 提取得到年、月和日的光变数据， 作为长时标周期估算的光变曲线; 同时选取3段连续性很好， 几乎是等间隔的以“min”为单位的观测数据， 作为对该星体进行短时标周期估算的数据样本. 采用最新发展起来的经验模态分解与传统的AR模型谱估计的组合方法， 将具有非线性、非平稳特性的光变数据分解为具有单分量信号特征的模态函数， 再选取相关系数高的分量重构信号后估算AR谱， 以此计算光变曲线的周期. 通过对类星体3C 273的光变资料分析， 我们得到类星体3C 273的长时标周期为:21.23、13.51、11.02、5.51、4.69、3.79、2.76 yr，与此前文献[33-34， 38-40]计算的结果一致; 短时标周期为: (30±1)、(15±0.3)、(7.5±0.2)、(10±0.1)、(5±0.6)和(6±0.4)、(3±0.5) min， 这些可能的周期还存在3组大致倍周期的关系. 分析结果也表明， 基于EMD-AR谱的方法对计算Blazar天体的光变周期是可靠的， 为天体的周期光变研究提供了一条新思路.