对长周期光变的研究是分析耀变体性质的重要方法之一,它的确定关系到耀变体的结构和辐射等相关的多个物理量的估算,比如辐射区域的半径、喷流的多普勒因子和黑洞质量等,长周期光变为理论模型的建立提供一定的参数。光变可分为长时标光变、中等时标光变、短时标光变[1]。对于长时标光变和中等时标光变来说,由于受到观测仪器、月相和天气等多种因素的影响,很难得到比较完整的光变的观测数据序列[2]。在研究光学剧变类星体光变周期性时,通过研究长周期光变分析方法,可以获取周期光变分析的最佳参数,用有限的实际观测数据序列获得最佳周期估算值。
通常分析周期的方法有:功率谱法、Jurkevich方法、小波、结构函数等,这些方法的特点是可以对高于奈奎斯特采样定律均匀采样的数据进行可靠分析[3-5]。但是实际的数据分析,特别是针对天体观测所获得耀变体的长时标光变周期分析中,这些方法的使用受到许多条件的限制,比如说傅里叶分析[6]要求连续的等间隔采样,缺失数据点的处理引进了一些不真实的信息。所以这些方法应用在天体的光变周期分析中,增大了确定周期的误差。而文[7]提出的Jurkevich方法,是一种建立在期待值均方误差基础上的统计方法,对于处理不等间隔观测数据有很大帮助[7]。本文利用4种方法对这一问题进行讨论,并提出相应的改进,使得在不等间隔时间序列的类星体光变曲线中,寻找周期性的分析计算变得更为简便准确。
1 类星体长周期光变的计算方法 1.1 Jurkevich方法文[7]基于天文测量中的非均匀测量问题提出了一种统计方法,文[8]曾用该方法分析蝎虎天体PKS 0735+178的光变特性。假设观测样本数据为N个,Xi为单次测量值,X为所有测量值的平均值,V2为测量数据样本方差,S为样本标准偏差,则有
$ \bar X = \frac{1}{N}\sum\limits_{i = 1}^N {{x_i}} , $ | (1) |
$ {V^2} = \sum\limits_{i = 1}^N {X_i^2 - N\overline {{X^2}} } , $ | (2) |
$ S = \sqrt {{V^2}/\left( {N - 1} \right)} . $ | (3) |
若样本划分为m组,则第l组的统计参数应为
$ \overline {{X_l}} = \frac{1}{{{m_l}}}\sum\limits_{i = 1}^{{m_l}} {{X_i}} , $ | (4) |
$ V_l^2 = \sum\limits_{i = 1}^{{m_l}} {X_i^2} - {m_l}\overline {X_l^2} , $ | (5) |
$ S_l^2 = V_l^2/\left( {{m_l} - 1} \right). $ | (6) |
对应m组的总方差为
$ V_m^2 = \sum\limits_{i = 1}^m {V_l^2} . $ | (7) |
显然,Vm2<V2。上述公式表明,当自变量相同的n个样本计算后则变为一个样本。如果分组的平均值相等,那么期待值方差等于分组的平均方差;如果分组的平均值不同,那么期待值方差超过分组的平均方差[7]。
在周期分析中按照不同的试验周期P把数据折叠,然后把折叠好的资料按相位从小到大排列,分成m组。若给定一个观测时间为t,试验周期为P,则其对应的相位定义为[7]
$ \phi \left( t \right) = \bmod \left( {\frac{{t - {t_0}}}{P},1} \right), $ | (8) |
其中,ϕ(t)满足0≤ϕ(t)≤1;t0为时间原点。计算出每组数据的方差Vl2和总方差Vm2,如果所得的实验周期等于实际周期,则Vm2的值将达到最小,由实验周期Pn与Vm2的关系曲线中的Vm2最小值对应的时间即可得到样本的周期。
1.2 时间补偿离散傅里叶变换分析方法时间补偿离散傅里叶变换方法是计算光变周期最常用的方法之一,此方法最早由文[9]在1981年提出。在过去的研究中文[10]用该方法分析PKS 1510-089红外光变周期。通过对1,sinωt,cosωt作Gram-Schmidt正交化,得到3个正交向量,将数据投影到3个正交向量上得到了频谱。具体过程如下:
$ {H_0} = 1;{H_1} = \cos \omega t;{H_2} = \sin \omega t. $ | (9) |
正交化后:
$ {h_0} = {a_0}{H_0}, $ | (10) |
$ {h_1} = {a_1}{H_1} - {a_1}{a_0}\left( {{h_0},{H_1}} \right), $ | (11) |
$ {h_2} = {a_2}{H_2} - {a_2}{a_0}\left( {{h_0},{H_2}} \right) - {a_2}{a_1}\left( {{h_1},{H_2}} \right), $ | (12) |
小括号表示求两个向量的内积。根据以下关系确定h0,h1,h2:
$ \left( {{h_0},{h_0}} \right) = \left( {{h_1},{h_1}} \right) = \left( {{h_2},{h_2}} \right) = 1. $ | (13) |
在不均匀采样的情况下,加权的时间补偿离散傅里叶变换:
$ F\left( \omega \right) = \left( {f,{h_1} + i{h_2}} \right)/{a_0}\sqrt 2 , $ | (14) |
在许多类星体的观测中,观测数据f(ti)的精度各不相同,考虑到这个问题,引进权重方程ωi=ω(ti),重新定义内积:
$ \left( {{g_1},{g_2}} \right) = \sum {{\omega _i}{g_1}\left( {{t_i}} \right){g_2}\left( {{t_i}} \right)} . $ | (15) |
在内积中引入权重后得到:
$ a_0^{ - 2} = \sum {{\omega _i}} , $ | (16) |
$ a_1^{ - 2} = \sum {{\omega _i}{{\left( {\cos {x_i}} \right)}^2}} - a_0^2{\left( {\sum {{\omega _i}\cos {x_i}} } \right)^2}, $ | (17) |
$ a_2^{ - 2} = \sum {{\omega _i}{{\left( {\sin {x_i}} \right)}^2}} - a_0^2{\left( {\sum {{\omega _i}\sin {x_i}} } \right)^2} - a_1^2{\left[ {\sum {{\omega _i}\sin {x_i}\cos {x_i}} - a_0^2\left( {\sum {{\omega _i}\sin {x_i}} } \right)\\ \left( {\sum {{\omega _i}\cos {x_i}} } \right)} \right]^2} $ | (18) |
在内积中引入权重后得到回归系数:
$ {c_0} = 0, $ | (19) |
$ {c_1} = {a_1}\sum {{\omega _i}f\left( {{t_i}} \right)\cos {x_i}} , $ | (20) |
$ {c_2} = {a_2}\sum {{\omega _i}f\left( {{t_i}} \right)\sin {x_i}} - {a_1}{a_2}{a_0}\left[ {\sum {{\omega _i}\sin {x_i}\cos {x_i}} - a_0^2\left( {\sum {{\omega _i}\sin {x_i}} } \right)\left( {\sum {{\omega _i}\cos {x_i}} } \right)} \right]. $ | (21) |
频率ω处的强度由下式给出:
$ I\left( \omega \right) = c_1^2 + c_2^2. $ | (22) |
由线性回归理论可知0≤I(ω)≤Q,式中:
$ Q = \left( {f,f} \right) = \sum {{\omega _i}f{{\left( {{t_i}} \right)}^2}} , $ | (23) |
利用这一性质,引进标准化因子:统计量S(ω)=I(ω)/Q,称这个量为谱相关系数,对于所有频率ω,0≤S(ω)≤1。
1.3 离散相关分析方法离散相关分析方法是Edelson和Krolik用来分析两组离散数据相关性的方法之一,文[11]用此方法分析PKS 2155-304的光变。该方法可以指出两个变量时间序列的相关性和时间延迟,并可应用于周期性分析[11]。具体步骤如下:
首先计算两组数据的离散相关函数值,如数组ai和bi,则离散相关函数值为
$ UDC{F_{ij}} = \frac{{\left( {{a_i} - \bar a} \right)\left( {{b_j} - \bar b} \right)}}{{\sqrt {\sigma _{\rm{a}}^2 \times \sigma _{\rm{b}}^2} }}, $ | (24) |
其中,a和b分别为两组数据的平均值;σa和σb分别为对应的标准偏差。
其次,计算DCF(τ)值。通过时间延迟Δtij=ti-tj把两组数联系起来,假如时间延迟为τ,在区间τ ± Δτ/2中有M个Δtij,则DCF(τ)值为
$ DCF\left( \tau \right) = \frac{1}{M}\sum {UDC{F_{ij}}\left( \tau \right)} . $ | (25) |
再次离散相关函数的误差为
$ \sigma \left( \tau \right) = \frac{1}{{M - 1}}{\left\{ {\sum {{{\left[ {UDC{F_{ij}} - DCF\left( \tau \right)} \right]}^2}} } \right\}^{\frac{1}{2}}}. $ | (26) |
对于所得到的离散相关图,如果峰值在0的右边,表明数组ai早于数组bi的变化。反之,数组ai迟于数组bi的变化。
1.4 功率谱密度分析方法功率谱密度的定义[12]是如果u(t)是一个可以进行傅里叶变换的函数,则
$ u\left( \nu \right) = \sum {u\left( t \right){{\rm{e}}^{ - 2i{c_v}t}}} , $ | (27) |
因为u(t)是实函数,u(ν)是一个复函数,它们之间满足Parseval公式
$ \sum {{u^t}\left( t \right)} = \sum {{{\left| {u\left( \nu \right)} \right|}^2}} , $ | (28) |
若u(t)表示光谱,则等式左端表示u(t)在(-∞, ∞)上的总能量,右端的函数称为能谱密度,它是一个非负实数,表示单位频率所具有的能量。从数学意义上讲,大多数u(t)不能进行傅里叶变换,如果使用截尾函数ur(t)截取u(t),
$ ur\left( t \right) = \left\{ {\begin{array}{*{20}{c}} {u\left( t \right),}&{\left| t \right| \le T}\\ {0,}&{\left| t \right| > T} \end{array}} \right., $ | (29) |
那么对于持续时间有限的截尾函数ur(t)而言可以进行傅里叶变换,变换式:
$ {u_T}\left( \nu \right) = \sum {{u_T}\left( t \right){{\rm{e}}^{2i{c_\nu }t}}} = \sum {u\left( t \right){{\rm{e}}^{2i{c_\nu }t}}} , $ | (30) |
同样它也满足Parseval公式:
$ \sum {u_T^2\left( t \right)} = \sum {{{\left| {{u_T}\left( v \right)} \right|}^2}} . $ | (31) |
根据截尾函数的定义,将上式两端除以2T并令T→∞,得到
$ \frac{1}{{2T}}\sum {{u^2}\left( t \right)} = \sum {\frac{1}{{2T}}{{\left| {ur\left( v \right)} \right|}^2}} , $ | (32) |
与能谱密度的定义相对应,上式右端函数称为功率谱密度,记为
$ P\left( \nu \right) = \mathop {\lim }\limits_{T \to \infty } \frac{{{{\left| {{u_T}\left( \nu \right)} \right|}^2}}}{{2T}}. $ | (33) |
从表达式中可见功率谱密度是一个非负实数,从整个功率谱密度的推导可以看出,功率谱密度是一个频域中的量,它直接反应了在频域中不同频率对应的值。
2 天文观测数据中周期信号的模拟检验 2.1 天文观测周期信号的模拟检验结果为检验上述4种研究方法的可靠性,用一个模拟的周期信号作为天文观测数据分析上述4种方法。在这里以2π为周期的正弦函数y=sinθ,检验分析4种研究方法的准确性,使用的单位为弧度。在实验中选择0 rad为起点,步长为0.1 rad,不同的数据点个数共取15组,所有研究方法考虑噪声等因素的影响。
图 1为检验Jurkevich方法分析天体周期的可靠性,取80~360个数据点,每组增加20个数据点,通过Jurkevich方法分析后得到15组不同数据点的Vm2-P曲线,进而得到的分析结果为(6.30 ± 0.02) rad。
图 2为检验时间补偿离散傅里叶变换分析方法分析天体周期的可靠性,取80~360个数据点,每组增加20个数据点,用时间补偿离散傅里叶变换分析方法分析后得到15组不同数据点的DCDFT-Frequency曲线,进而得到的分析结果为(6.33 ± 0.13) rad。
图 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方法的最佳参数利用Jurkevich方法分析周期时,分别取360、720和1 500个数据点,同时取分组数1~12组、试验周期P=2π,令步长为0.1 rad,计算得到不同x(x=P)所对应的y(y=Vm2)值。
如图 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组。下面用具有多个观测数据的两个源验证此结果。
3 计算并分析类星体3C 279和3C 454.3的光变周期 3.1 类星体3C 279的数据点及光变周期本文研究的类星体3C 279在B波段和3C 454.3在B波段数据点主要是从网站(http://www.astro.yale.<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与Vm2的关系曲线中有很多的极小值,为了有效地判断周期Pn与Vm2关系曲线中所得出周期的真实性,文[14]给出了较好的判据:
$ f = \frac{{1 - V_m^2}}{{V_m^2}}, $ | (34) |
其中,Vm2为归一化的值。在归一化的Vm2-Pn图中,若Vm2=1.0,则有f=0,即在样本数据中没有表现出周期性;若Vm2=0,则有f=1,这时样本中存在的最大周期能够从图中识别。在通常情况下,当f ≥ 0.5时,样本数据表现出非常强的周期性;当f ≤ 0.25时,则表现为不含有周期性。进一步分析Vm2值,选择曲线较平滑部分的最小深度和噪声做实验,如果平滑部分的相对最小爆发值比平滑部分的Vm2大10倍,则可以进一步讨论相应的短周期性。在计算中分组数m的大小十分重要,组数m分得越多,灵敏度越高,但每组数据中有少数的数据点在图中产生较大的噪声。分组数m较大时,还增加Vm2的计算量。反之,有可能寻找不到周期。如图 10,对类星体3C 279在B波段的光变周期分析中,很容易找到Vm2及P值,并用文[14]的判据:f=(1-Vm2)/Vm2进行检验,结果如表 1。
P | Vm2 | f=(1-Vm2)/Vm2 | T/a |
413 | 0.700 6 | 0.427 3 | 1.13 |
829 | 0.529 | 0.860 4 | 2.27 |
1 222 | 0.352 1 | 1.840 1 | 3.35 |
通过对表 1分析发现,最小值Vm2=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,说明它们之间可能存在天文学倍频关系。
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的真正的光变周期。
如图 12,对类星体3C 454.3在B波段的光变周期分析中,很容易找到Vm2及P值,并用文[14]的判据:f=(1-Vm2)/Vm2进行检验,结果如表 2。
通过对表 2分析发现,当Vm2取最小值0.550 5,f=0.817,对应的光变周期可能为P=457 d,这既满足了图 12的分析结果,也满足f判据[14]。
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值,Vm2与Pn的关系曲线的倾斜程度不同,周期大的向下倾斜。结合f检验公式可知,Vm2越小,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倍。
[1] | 谢照华.耀变体的吸积盘与喷流关联及粘滞系数的研究[D].昆明: 云南大学, 2012. http://cdmd.cnki.com.cn/Article/CDMD-10673-1013134637.htm |
[2] | 王文广.耀变体的光变特性及色指数分析[D].昆明: 云南师范大学, 2017. http://cdmd.cnki.com.cn/Article/CDMD-10681-1017730514.htm |
[3] | FAN J H, ADAM G, XIE G Z, et al. Correlation between the gamma-ray and the radio emissions[J]. Astronomy & Astrophysics, 1998, 338(1): 27–30. |
[4] | FAN J H, KURTANIDZE O, LIU Y, et al. Optical monitoring of two brightest nearby quasars, PHL 1811 and 3C 273[J]. The Astrophysical Journal Supplement, 2014, 213(2): article id. 26. DOI: 10.1088/0067-0049/213/2/26 |
[5] | RIEGER F M, MANNHEIM K. Implications of a possible 23 day periodicity for binary black hole models in Mkn~501[J]. Astronomy & Astrophysics, 2012, 359(3): 948–952. |
[6] | WEBB J R, SMITH A G, LEACOCK R J, et al. Optical observations of 22 violently variable extragalactic sources-1968-1986[J]. The Astronomical Journal, 1988, 95(2): 374–397. |
[7] | JURKEVICH I. A method of computing periods of cyclic phenomena[J]. Astrophysics and Space Science, 1971, 13(1): 154–167. DOI: 10.1007/BF00656321 |
[8] | 余莲, 张雄, 王文广, 等. 蝎虎天体PKS 0735+178的光变特性分析[J]. 天文研究与技术, 2018, 15(1): 10–16 |
[9] | FERRAZ-MELLO S. Estimation of periods from unequally spaced observations[J]. The Astronomical Journal, 1981, 86(4): 619. |
[10] | 罗双玲, 张雄, 王文广. 时间补偿离散傅里叶变换分析PKS 1510-089红外光变周期[J]. 云南师范大学学报(自然科学版), 2016, 36(5): 1–4 DOI: 10.7699/j.ynnu.ns-2016-057 |
[11] | FAN J H, LIN R G. The variability analysis of PKS 2155-304[J]. Astronomy & Astrophysics, 2000, 355(3): 880–884. |
[12] | 张蓉竹, 蔡邦维, 杨春林, 等. 功率谱密度的数值计算方法[J]. 强激光与粒子束, 2000(6): 661–664 |
[13] | 张雄, 谢光中, 白金明. 用Jurkevich方法分析计算BL Lac天体的光变周期[J]. 天体物理学报, 1998(3): 33–41 |
[14] | KIDGER M R, TAKALO L, SILLANP A. A new analysis of the 11-year period in OJ287:confirmation of its existence[J]. Astronomy and Astrophysics, 1992, 264: 32–361. |
[15] | PYATUNINA T B, KUDRYAVTSEVA N A, GABUZDA D C, et al. Frequency-dependent time delays for strong outbursts in selected blazars from the Metsahovi and UMRAO monitoring data bases-Ⅱ[J]. Monthly Notices of the Royal Astronomical Society, 2007, 381(2): 797–808. DOI: 10.1111/j.1365-2966.2007.12281.x |
[16] | HEIDT J, WAGNER S J. Statistics of optical intraday variability in a complete sample of radio-selected BL Lac objects[J]. Astronomy and Astrophysics, 1996, 305: 42–52. |
[17] | 刘凯博, 杨东末, 侯德东, 等. Jurkevich方法在3C 120天体光变周期分析中的应用[J]. 云南师范大学学报(自然科学版), 2009, 29(5): 1–5, 16 DOI: 10.3969/j.issn.1007-9793.2009.05.001 |
[18] | 李怀珍. Blazar天体的光变及能谱分布的研究[D].昆明: 云南大学, 2011. http://cdmd.cnki.com.cn/article/cdmd-10673-1011212105.htm |
[19] | 苏成悦. 类星体PKS 2251+158的长期光学光变周期分析[J]. 天体物理学报, 2000, 20(1): 11–16 |