类星体长周期光变分析方法的研究*

余 莲，张 雄，张皓晶，李富婷，徐小林，任国伟，吴月承

(云南师范大学物理与电子信息学院，云南 昆明 650500)

对长周期光变的研究是分析耀变体性质的重要方法之一，它的确定关系到耀变体的结构和辐射等相关的多个物理量的估算，比如辐射区域的半径、喷流的多普勒因子和黑洞质量等，长周期光变为理论模型的建立提供一定的参数。光变可分为长时标光变、中等时标光变、短时标光变[1]。对于长时标光变和中等时标光变来说，由于受到观测仪器、月相和天气等多种因素的影响，很难得到比较完整的光变的观测数据序列[2]。在研究光学剧变类星体光变周期性时，通过研究长周期光变分析方法，可以获取周期光变分析的最佳参数，用有限的实际观测数据序列获得最佳周期估算值。

通常分析周期的方法有：功率谱法、Jurkevich方法、小波、结构函数等，这些方法的特点是可以对高于奈奎斯特采样定律均匀采样的数据进行可靠分析[3-5]。但是实际的数据分析，特别是针对天体观测所获得耀变体的长时标光变周期分析中，这些方法的使用受到许多条件的限制，比如说傅里叶分析[6]要求连续的等间隔采样，缺失数据点的处理引进了一些不真实的信息。所以这些方法应用在天体的光变周期分析中，增大了确定周期的误差。而文[7]提出的Jurkevich方法，是一种建立在期待值均方误差基础上的统计方法，对于处理不等间隔观测数据有很大帮助[7]。本文利用4种方法对这一问题进行讨论，并提出相应的改进，使得在不等间隔时间序列的类星体光变曲线中，寻找周期性的分析计算变得更为简便准确。

1 类星体长周期光变的计算方法

1.1 Jurkevich方法

(7)

显然，V2m

在周期分析中按照不同的试验周期P把数据折叠，然后把折叠好的资料按相位从小到大排列，分成m组。若给定一个观测时间为t，试验周期为P，则其对应的相位定义为[7]

(8)

其中，φ(t)满足0≤φ(t)≤1；t0为时间原点。计算出每组数据的方差V2l和总方差V2m，如果所得的实验周期等于实际周期，则V2m的值将达到最小，由实验周期Pn与V2m的关系曲线中的V2m最小值对应的时间即可得到样本的周期。

1.2 时间补偿离散傅里叶变换分析方法

时间补偿离散傅里叶变换方法是计算光变周期最常用的方法之一，此方法最早由文[9]在1981年提出。在过去的研究中文[10]用该方法分析PKS 1510-089红外光变周期。通过对1，sinωt，cosωt作Gram-Schmidt正交化，得到3个正交向量，将数据投影到3个正交向量上得到了频谱。具体过程如下：

H0=1；H1=cosωt；H2=sinωt.

(9)

正交化后：h0=a0H0，

(10)

h1=a1H1-a1a0(h0,H1)，

(11)

h2=a2H2-a2a0(h0,H2)-a2a1(h1,H2)，

(12)

小括号表示求两个向量的内积。根据以下关系确定h0，h1，h2：

(h0,h0)=(h1,h1)=(h2,h2)=1 .

(13)

在不均匀采样的情况下，加权的时间补偿离散傅里叶变换：

(14)

在许多类星体的观测中，观测数据f(ti)的精度各不相同，考虑到这个问题，引进权重方程ωi=ω(ti)，重新定义内积：

(g1,g2)=∑ωig1(ti)g2(ti).

(15)

(16)

(17)

(18)

在内积中引入权重后得到回归系数：c0=0，

(19)

c1=a1∑ωif(ti)cosxi，

(20)

(21)

(22)

由线性回归理论可知0≤I(ω)≤Q，式中：

Q=(f,f)=∑ωif(ti)2，

(23)

利用这一性质，引进标准化因子：统计量S(ω)=I(ω)/Q，称这个量为谱相关系数，对于所有频率ω，0≤S(ω)≤1。

1.3 离散相关分析方法

离散相关分析方法是Edelson和Krolik用来分析两组离散数据相关性的方法之一，文[11]用此方法分析PKS 2155-304的光变。该方法可以指出两个变量时间序列的相关性和时间延迟，并可应用于周期性分析[11]。具体步骤如下：

首先计算两组数据的离散相关函数值，如数组ai和bi，则离散相关函数值为

(24)

其次，计算DCF(τ)值。通过时间延迟Δtij=ti-tj把两组数联系起来，假如时间延迟为τ，在区间τ± Δτ/2中有M个Δtij，则DCF(τ)值为

(25)

再次离散相关函数的误差为

(26)

对于所得到的离散相关图，如果峰值在0的右边，表明数组ai早于数组bi的变化。反之，数组ai迟于数组bi的变化。

1.4 功率谱密度分析方法

功率谱密度的定义[12]是如果u(t)是一个可以进行傅里叶变换的函数，则

u(ν)=∑u(t)e-2icvt，

(27)

因为u(t)是实函数，u(ν)是一个复函数，它们之间满足Parseval公式

(28)

若u(t)表示光谱，则等式左端表示u(t)在(-,)上的总能量，右端的函数称为能谱密度，它是一个非负实数，表示单位频率所具有的能量。从数学意义上讲，大多数u(t)不能进行傅里叶变换，如果使用截尾函数ur(t)截取u(t)，

(29)

那么对于持续时间有限的截尾函数ur(t)而言可以进行傅里叶变换，变换式：

uT(ν)=∑uT(t)e2icνt=∑u(t)e2icνt，

(30)

(31)

根据截尾函数的定义，将上式两端除以2T并令T→，得到

(32)

与能谱密度的定义相对应，上式右端函数称为功率谱密度，记为

(33)

从表达式中可见功率谱密度是一个非负实数，从整个功率谱密度的推导可以看出，功率谱密度是一个频域中的量，它直接反应了在频域中不同频率对应的值。

2 天文观测数据中周期信号的模拟检验

2.1 天文观测周期信号的模拟检验结果

为检验上述4种研究方法的可靠性，用一个模拟的周期信号作为天文观测数据分析上述4种方法。在这里以2π为周期的正弦函数y=sinθ，检验分析4种研究方法的准确性，使用的单位为弧度。在实验中选择0 rad为起点，步长为0.1 rad，不同的数据点个数共取15组，所有研究方法考虑噪声等因素的影响。

图1为检验Jurkevich方法分析天体周期的可靠性，取80～360个数据点，每组增加20个数据点，通过Jurkevich方法分析后得到15组不同数据点的V2m-P曲线，进而得到的分析结果为(6.30 ± 0.02) rad。

图1 Jurkevich方法对sin函数不同数据点的周期性分析
Fig.1 Periodic analysis of sin functions′ different data points by using Jurkevich method

图2为检验时间补偿离散傅里叶变换分析方法分析天体周期的可靠性，取80～360个数据点，每组增加20个数据点，用时间补偿离散傅里叶变换分析方法分析后得到15组不同数据点的DCDFT-Frequency曲线，进而得到的分析结果为(6.33 ± 0.13) rad。

图2 时间补偿离散傅里叶变换分析方法(DCDFT)对sin函数不同数据点的周期性分析
Fig.2 Periodic analysis of sin functions′ different data points by using DCDFT method

图3为检验离散相关分析方法分析天体周期的可靠性，取80～360个数据点，每组增加20个数据点，用离散相关分析方法分析后得到15组不同数据点的DCF-Delay曲线，进而得到的分析结果为(6.287 ± 0.088) rad。

图4为检验功率谱密度分析方法分析天体周期的可靠性，取80～360个数据点，每组增加20个数据点，用功率谱密度分析方法分析后得到15组不同数据点的Power-Frequency曲线，进而得到的分析结果为(6.287 ± 0.045) rad。

由图5得到，Jurkevich方法、时间补偿离散傅里叶变换分析方法、离散相关分析方法和功率谱密度分析方法周期需要的数据点最低分别为50个、120个、100个和60个。获取最短的连续数据采样后，Jurkevich方法最有效。

2.2 利用天文模拟数据寻找Jurkevich方法的最佳参数

图3 离散相关分析方法(DCF)对sin函数不同数据点的周期性分析
Fig.3 Periodic analysis of sin functions′ different data points by using DCF method

如图6，当取360个数据点，分组数为1～6组时，要求不等间隔的数据样本周期数不少于6个[7]，因此周期性可忽略；当分组数为7～12组时，在第10组、第11组和第12组可能出现了倍周期，即可能刚好是第1个周期的重复出现。即使出现极小微差的伪周期，周期性最好的仍为第9组。

当取720个数据点，如图7，分组数为1～6组时，同样没有周期性；当分组数为7～12组时，在第10组和第12组出现了倍周期，即可能刚好是第1个周期的重复出现；在第11组出现了明显的伪周期，因此第9组周期性的分析结果最佳。

如图8，当取1 500个数据点，分组数为1～6组时，周期性不明显，所以这几组的周期可忽略。当分组数为7～12组时，在第8组、第10组和第12组出现了倍周期，第11组的周期性显然没有第9组明显。总的来说，周期性的最佳分析结果为第9组。下面用具有多个观测数据的两个源验证此结果。

图4 功率谱密度分析方法(PSD)对sin函数不同数据点的周期性分析
Fig.4 Periodic analysis of sin functions′ different data points by using PSD method

图5 sin函数用时间补偿离散傅里叶变换分析方法 (DCDFT)、Jurkevich方法、离散相关分析方法 (DCF) 和功率谱密度分析方法 (PSD) 进行周期分析时要求的最低数据点 (右边的两个图形从上向下分别为Jurkevich方法和功率谱密度分析方法(PSD) 的分析图形)

Fig.5 The minimum data points required for the sin function to perform periodic analysis using the DCDFT method, Jurkevich method, DCF method and PSD method (The two graphs on the right are the Jurkevich method and PSD method analysis graph from top to bottom)

图6 sin函数的V2m-P，左图分组数为1～6组，右图分组数为7～12组
Fig.6 TheV2m-Pcurve of sin, The number of groupings on the left are Group 1-6, and the grouping on the right are Group 7-12

图7 sin函数的V2m-P，左图分组数为1～6组，右图分组数为7～12组
Fig.7 TheV2m-Pcurve of sin, The number of groupings on the left are Group 1-6, and the grouping on the right are Group 7-12

图8 sin函数的V2m-P，左图分组数为1～6组，右图分组数为7～12组

Fig.8 TheV2m-Pcurve of sin, The number of groupings on the left are Group 1-6,and the grouping on the right are Group 7-12

3 计算并分析类星体3C 279和3C 454.3的光变周期

3.1 类星体3C 279的数据点及光变周期

本文研究的类星体3C 279在B波段和3C 454.3在B波段数据点主要是从网站(http://www.astro.yale.edu/smarts/glast/home.php)获取。

从2008年2月5日到2017年7月19日，由以上网站搜寻到类星体3C 279在B波段的数据点有661个，B波段的光变曲线如图9，可以看出有几次大的爆发，最亮时B波段星等为14.9，光变曲线中有6个以上周期存在，这些均满足Jurkevich方法中确认长周期存在的必要条件[13]。

利用Jurkevich方法分组数m=9，图10可得到3C 279在B波段的光变周期(2.81 ± 0.54)年。周期Pn与V2m的关系曲线中有很多的极小值，为了有效地判断周期Pn与V2m关系曲线中所得出周期的真实性，文[14]给出了较好的判据：

(34)

其中，V2m为归一化的值。在归一化的V2m-Pn图中，若V2m=1.0，则有f=0，即在样本数据中没有表现出周期性；若V2m=0，则有f=1，这时样本中存在的最大周期能够从图中识别。在通常情况下，当f≥ 0.5时，样本数据表现出非常强的周期性；当f≤ 0.25时，则表现为不含有周期性。进一步分析V2m值，选择曲线较平滑部分的最小深度和噪声做实验，如果平滑部分的相对最小爆发值比平滑部分的V2m大10倍，则可以进一步讨论相应的短周期性。在计算中分组数m的大小十分重要，组数m分得越多，灵敏度越高，但每组数据中有少数的数据点在图中产生较大的噪声。分组数m较大时，还增加V2m的计算量。反之，有可能寻找不到周期。如图10，对类星体3C 279在B波段的光变周期分析中，很容易找到V2m及P值，并用文[14]的判据：f=(1-V2m)/V2m进行检验，结果如表1。

通过对表1分析发现，最小值V2m=0.352 1，f=1.840 1，类星体3C 279可能的光变周期为(2.81 ± 0.54)年。利用Jurkevich方法，分组数m=9，分析了类星体3C 279在B波段的光变数据，图10显示了B波段的分析结果。从表1可以看出，在满足判据f的条件下[14]，类星体3C 279在B波段的P2，P3与P1之间存在类星体简单的倍数关系：P2≈2P1，P3≈3P1，说明它们之间可能存在天文学倍频关系。

图9 类星体3C 279在B波段的光变周期

Fig.9 The light curve of quasar 3C 279 in B band

图10 Jurkevich方法分析类星体3C 279在B波段的光变周期

Fig.10 Jurkevich method analysis of the variability period of quasar 3C 279 in B band

3.2 类星体3C 454.3的数据点及光变周期

3C 454.3(PKS 2251+158, OY091)是一个比较亮、变化比较剧烈的类星体，并且在光学和射电波段存在比较明显的相关性[15]。本文的数据是从2008年6月23日到2017年7月30日，约765个数据点，获得了如图11的历史光变曲线。

利用Jurkevich方法，分组数m=9，分析了类星体3C 454.3在B波段的光变数据，图12显示了B波段的分析结果。从图12可以看出，类星体3C 454.3在B波段有3个明显的极小值，意味着类星体3C 454.3在B波段光变曲线中存在3个可能的周期，它们分别是：P1=457 d，P2=891 d和P3=1 320 d。根据文[16]的结果，要确定一个周期，数据样本的时间跨度要超过周期的6倍。文中数据样本的时间跨度大约2 500 d，小于P2，P3的6倍， 故P2，P3必须排除，并需要更多的观测数据确定它们。但是P2，P3与P1之间存在着简单的倍数关系：P2≈2P1，P3≈3P1，说明它们之间可能存在天文学倍频关系，意味着类星体3C 454.3在B波段的光变曲线中可能存在一个P1=457 d的真正的光变周期。

表1 类星体3C 279在B波段的周期分析表Table 1 The table of 3C 279 objects cycle analysis

图11 类星体3C 454.3在B波段的光变周期

Fig.11 The light curve of quasar 3C 454.3 in B band

图12 Jurkevich方法分析类星体3C 454.3在B波段的光变周期

Fig.12 Jurkevich method analysis of the variability period of quasar 3C 454.3 in B band

如图12，对类星体3C 454.3在B波段的光变周期分析中，很容易找到V2m及P值，并用文[14]的判据：f=(1-V2m)/V2m进行检验，结果如表2。

通过对表2分析发现，当V2m取最小值0.550 5，f=0.817，对应的光变周期可能为P=457 d，这既满足了图12的分析结果，也满足f判据[14]。

表2 类星体3C 454.3在B波段的周期分析表Table 2 The table of 3C 454.3 objects cycle analysis

4 讨论与结论

本文利用一个以2π为周期的正弦函数y=sinθ检验4种研究方法的可靠性。在实验中，使用的单位为弧度，选择0 rad为起点，步长为0.1 rad，不同的数据点个数共取15组，得到Jurkevich方法的分析结果为(6.3 ± 0.017) rad；时间补偿离散傅里叶变换分析方法的分析结果为(6.331 ± 0.130) rad；离散相关分析方法的分析结果为(6.287 ± 0.088) rad；功率谱密度方法的分析结果为(6.287 ± 0.045) rad。分析结果表明，Jurkevich方法、时间补偿离散傅里叶变换分析方法、离散相关分析方法和功率谱密度方法周期需要的数据点最低分别为50个、120个、100个和60个。获取最短的连续数据采样后，Jurkevich方法最有效。Jurkevich方法分析结果在4种方法中最精确可靠，且此计算方法简捷实用。

在实际应用中，发现Jurkevich方法十分依赖于分组数m，m越大分析结果越好，但会产生比较大的噪声；反之则有可能找不到周期[17]。如果数据分布不均匀，则会导致组间数据分布偏差很大，从而影响获得实际的周期。根据图6、图7和图8可以看出，对于不同的m值，V2m与Pn的关系曲线的倾斜程度不同，周期大的向下倾斜。结合f检验公式可知，V2m越小，f值越大，更容易满足f检验。由此可知f检验依赖于m值的大小。利用模拟数据寻找到Jurkevich方法的分组数为m=9时，分析结果最佳。最后用分组数m=9时的Jurkevich方法分析了类星体3C 279及3C 454.3的光变周期，得出类星体3C 279可能的光变周期为(2.81 ± 0.54)年，在文[18]中分析了类星体3C 279可能的光变周期为(130.6 ± 1.3) d，分析得到的可能周期大约是(130.6 ± 1.3) d的8倍。类星体3C 454.3可能的光变周期为457 d。在文[19]中分析了类星体3C 454.3可能的光变周期为12.39年，大约是本文得到周期的10倍。