余 蓮，張 雄，張皓晶，李富婷，徐小林，任國(guó)偉，吳月承

(云南師范大學(xué)物理與電子信息學(xué)院，云南 昆明 650500)

對(duì)長(zhǎng)周期光變的研究是分析耀變體性質(zhì)的重要方法之一，它的確定關(guān)系到耀變體的結(jié)構(gòu)和輻射等相關(guān)的多個(gè)物理量的估算，比如輻射區(qū)域的半徑、噴流的多普勒因子和黑洞質(zhì)量等，長(zhǎng)周期光變?yōu)槔碚撃Ｐ偷慕⑻峁┮欢ǖ膮?shù)。光變可分為長(zhǎng)時(shí)標(biāo)光變、中等時(shí)標(biāo)光變、短時(shí)標(biāo)光變[1]。對(duì)于長(zhǎng)時(shí)標(biāo)光變和中等時(shí)標(biāo)光變來說，由于受到觀測(cè)儀器、月相和天氣等多種因素的影響，很難得到比較完整的光變的觀測(cè)數(shù)據(jù)序列[2]。在研究光學(xué)劇變類星體光變周期性時(shí)，通過研究長(zhǎng)周期光變分析方法，可以獲取周期光變分析的最佳參數(shù)，用有限的實(shí)際觀測(cè)數(shù)據(jù)序列獲得最佳周期估算值。

通常分析周期的方法有：功率譜法、Jurkevich方法、小波、結(jié)構(gòu)函數(shù)等，這些方法的特點(diǎn)是可以對(duì)高于奈奎斯特采樣定律均勻采樣的數(shù)據(jù)進(jìn)行可靠分析[3-5]。但是實(shí)際的數(shù)據(jù)分析，特別是針對(duì)天體觀測(cè)所獲得耀變體的長(zhǎng)時(shí)標(biāo)光變周期分析中，這些方法的使用受到許多條件的限制，比如說傅里葉分析[6]要求連續(xù)的等間隔采樣，缺失數(shù)據(jù)點(diǎn)的處理引進(jìn)了一些不真實(shí)的信息。所以這些方法應(yīng)用在天體的光變周期分析中，增大了確定周期的誤差。而文[7]提出的Jurkevich方法，是一種建立在期待值均方誤差基礎(chǔ)上的統(tǒng)計(jì)方法，對(duì)于處理不等間隔觀測(cè)數(shù)據(jù)有很大幫助[7]。本文利用4種方法對(duì)這一問題進(jìn)行討論，并提出相應(yīng)的改進(jìn)，使得在不等間隔時(shí)間序列的類星體光變曲線中，尋找周期性的分析計(jì)算變得更為簡(jiǎn)便準(zhǔn)確。

1 類星體長(zhǎng)周期光變的計(jì)算方法

1.1 Jurkevich方法

(7)

顯然，V2m

在周期分析中按照不同的試驗(yàn)周期P把數(shù)據(jù)折疊，然后把折疊好的資料按相位從小到大排列，分成m組。若給定一個(gè)觀測(cè)時(shí)間為t，試驗(yàn)周期為P，則其對(duì)應(yīng)的相位定義為[7]

(8)

其中，φ(t)滿足0≤φ(t)≤1；t0為時(shí)間原點(diǎn)。計(jì)算出每組數(shù)據(jù)的方差V2l和總方差V2m，如果所得的實(shí)驗(yàn)周期等于實(shí)際周期，則V2m的值將達(dá)到最小，由實(shí)驗(yàn)周期Pn與V2m的關(guān)系曲線中的V2m最小值對(duì)應(yīng)的時(shí)間即可得到樣本的周期。

1.2 時(shí)間補(bǔ)償離散傅里葉變換分析方法

時(shí)間補(bǔ)償離散傅里葉變換方法是計(jì)算光變周期最常用的方法之一，此方法最早由文[9]在1981年提出。在過去的研究中文[10]用該方法分析PKS 1510-089紅外光變周期。通過對(duì)1，sinωt，cosωt作Gram-Schmidt正交化，得到3個(gè)正交向量，將數(shù)據(jù)投影到3個(gè)正交向量上得到了頻譜。具體過程如下：

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)

小括號(hào)表示求兩個(gè)向量的內(nèi)積。根據(jù)以下關(guān)系確定h0，h1，h2：

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

(13)

在不均勻采樣的情況下，加權(quán)的時(shí)間補(bǔ)償離散傅里葉變換：

(14)

在許多類星體的觀測(cè)中，觀測(cè)數(shù)據(jù)f(ti)的精度各不相同，考慮到這個(gè)問題，引進(jìn)權(quán)重方程ωi=ω(ti)，重新定義內(nèi)積：

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

(15)

(16)

(17)

(18)

在內(nèi)積中引入權(quán)重后得到回歸系數(shù)：c0=0，

(19)

c1=a1∑ωif(ti)cosxi，

(20)

(21)

(22)

由線性回歸理論可知0≤I(ω)≤Q，式中：

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

(23)

利用這一性質(zhì)，引進(jìn)標(biāo)準(zhǔn)化因子：統(tǒng)計(jì)量S(ω)=I(ω)/Q，稱這個(gè)量為譜相關(guān)系數(shù)，對(duì)于所有頻率ω，0≤S(ω)≤1。

1.3 離散相關(guān)分析方法

離散相關(guān)分析方法是Edelson和Krolik用來分析兩組離散數(shù)據(jù)相關(guān)性的方法之一，文[11]用此方法分析PKS 2155-304的光變。該方法可以指出兩個(gè)變量時(shí)間序列的相關(guān)性和時(shí)間延遲，并可應(yīng)用于周期性分析[11]。具體步驟如下：

首先計(jì)算兩組數(shù)據(jù)的離散相關(guān)函數(shù)值，如數(shù)組ai和bi，則離散相關(guān)函數(shù)值為

(24)

其次，計(jì)算DCF(τ)值。通過時(shí)間延遲Δtij=ti-tj把兩組數(shù)聯(lián)系起來，假如時(shí)間延遲為τ，在區(qū)間τ± Δτ/2中有M個(gè)Δtij，則DCF(τ)值為

(25)

再次離散相關(guān)函數(shù)的誤差為

(26)

對(duì)于所得到的離散相關(guān)圖，如果峰值在0的右邊，表明數(shù)組ai早于數(shù)組bi的變化。反之，數(shù)組ai遲于數(shù)組bi的變化。

1.4 功率譜密度分析方法

功率譜密度的定義[12]是如果u(t)是一個(gè)可以進(jìn)行傅里葉變換的函數(shù)，則

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

(27)

因?yàn)閡(t)是實(shí)函數(shù)，u(ν)是一個(gè)復(fù)函數(shù)，它們之間滿足Parseval公式

(28)

若u(t)表示光譜，則等式左端表示u(t)在(-,)上的總能量，右端的函數(shù)稱為能譜密度，它是一個(gè)非負(fù)實(shí)數(shù)，表示單位頻率所具有的能量。從數(shù)學(xué)意義上講，大多數(shù)u(t)不能進(jìn)行傅里葉變換，如果使用截尾函數(shù)ur(t)截取u(t)，

(29)

那么對(duì)于持續(xù)時(shí)間有限的截尾函數(shù)ur(t)而言可以進(jìn)行傅里葉變換，變換式：

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

(30)

(31)

根據(jù)截尾函數(shù)的定義，將上式兩端除以2T并令T→，得到

(32)

與能譜密度的定義相對(duì)應(yīng)，上式右端函數(shù)稱為功率譜密度，記為

(33)

從表達(dá)式中可見功率譜密度是一個(gè)非負(fù)實(shí)數(shù)，從整個(gè)功率譜密度的推導(dǎo)可以看出，功率譜密度是一個(gè)頻域中的量，它直接反應(yīng)了在頻域中不同頻率對(duì)應(yīng)的值。

2 天文觀測(cè)數(shù)據(jù)中周期信號(hào)的模擬檢驗(yàn)

2.1 天文觀測(cè)周期信號(hào)的模擬檢驗(yàn)結(jié)果

為檢驗(yàn)上述4種研究方法的可靠性，用一個(gè)模擬的周期信號(hào)作為天文觀測(cè)數(shù)據(jù)分析上述4種方法。在這里以2π為周期的正弦函數(shù)y=sinθ，檢驗(yàn)分析4種研究方法的準(zhǔn)確性，使用的單位為弧度。在實(shí)驗(yàn)中選擇0 rad為起點(diǎn)，步長(zhǎng)為0.1 rad，不同的數(shù)據(jù)點(diǎn)個(gè)數(shù)共取15組，所有研究方法考慮噪聲等因素的影響。

圖1為檢驗(yàn)Jurkevich方法分析天體周期的可靠性，取80～360個(gè)數(shù)據(jù)點(diǎn)，每組增加20個(gè)數(shù)據(jù)點(diǎn)，通過Jurkevich方法分析后得到15組不同數(shù)據(jù)點(diǎn)的V2m-P曲線，進(jìn)而得到的分析結(jié)果為(6.30 ± 0.02) rad。

圖1 Jurkevich方法對(duì)sin函數(shù)不同數(shù)據(jù)點(diǎn)的周期性分析
Fig.1 Periodic analysis of sin functions′ different data points by using Jurkevich method

圖2為檢驗(yàn)時(shí)間補(bǔ)償離散傅里葉變換分析方法分析天體周期的可靠性，取80～360個(gè)數(shù)據(jù)點(diǎn)，每組增加20個(gè)數(shù)據(jù)點(diǎn)，用時(shí)間補(bǔ)償離散傅里葉變換分析方法分析后得到15組不同數(shù)據(jù)點(diǎn)的DCDFT-Frequency曲線，進(jìn)而得到的分析結(jié)果為(6.33 ± 0.13) rad。

圖2 時(shí)間補(bǔ)償離散傅里葉變換分析方法(DCDFT)對(duì)sin函數(shù)不同數(shù)據(jù)點(diǎn)的周期性分析
Fig.2 Periodic analysis of sin functions′ different data points by using DCDFT method

圖3為檢驗(yàn)離散相關(guān)分析方法分析天體周期的可靠性，取80～360個(gè)數(shù)據(jù)點(diǎn)，每組增加20個(gè)數(shù)據(jù)點(diǎn)，用離散相關(guān)分析方法分析后得到15組不同數(shù)據(jù)點(diǎn)的DCF-Delay曲線，進(jìn)而得到的分析結(jié)果為(6.287 ± 0.088) rad。

圖4為檢驗(yàn)功率譜密度分析方法分析天體周期的可靠性，取80～360個(gè)數(shù)據(jù)點(diǎn)，每組增加20個(gè)數(shù)據(jù)點(diǎn)，用功率譜密度分析方法分析后得到15組不同數(shù)據(jù)點(diǎn)的Power-Frequency曲線，進(jìn)而得到的分析結(jié)果為(6.287 ± 0.045) rad。

由圖5得到，Jurkevich方法、時(shí)間補(bǔ)償離散傅里葉變換分析方法、離散相關(guān)分析方法和功率譜密度分析方法周期需要的數(shù)據(jù)點(diǎn)最低分別為50個(gè)、120個(gè)、100個(gè)和60個(gè)。獲取最短的連續(xù)數(shù)據(jù)采樣后，Jurkevich方法最有效。

2.2 利用天文模擬數(shù)據(jù)尋找Jurkevich方法的最佳參數(shù)

圖3 離散相關(guān)分析方法(DCF)對(duì)sin函數(shù)不同數(shù)據(jù)點(diǎn)的周期性分析
Fig.3 Periodic analysis of sin functions′ different data points by using DCF method

如圖6，當(dāng)取360個(gè)數(shù)據(jù)點(diǎn)，分組數(shù)為1～6組時(shí)，要求不等間隔的數(shù)據(jù)樣本周期數(shù)不少于6個(gè)[7]，因此周期性可忽略；當(dāng)分組數(shù)為7～12組時(shí)，在第10組、第11組和第12組可能出現(xiàn)了倍周期，即可能剛好是第1個(gè)周期的重復(fù)出現(xiàn)。即使出現(xiàn)極小微差的偽周期，周期性最好的仍為第9組。

當(dāng)取720個(gè)數(shù)據(jù)點(diǎn)，如圖7，分組數(shù)為1～6組時(shí)，同樣沒有周期性；當(dāng)分組數(shù)為7～12組時(shí)，在第10組和第12組出現(xiàn)了倍周期，即可能剛好是第1個(gè)周期的重復(fù)出現(xiàn)；在第11組出現(xiàn)了明顯的偽周期，因此第9組周期性的分析結(jié)果最佳。

如圖8，當(dāng)取1 500個(gè)數(shù)據(jù)點(diǎn)，分組數(shù)為1～6組時(shí)，周期性不明顯，所以這幾組的周期可忽略。當(dāng)分組數(shù)為7～12組時(shí)，在第8組、第10組和第12組出現(xiàn)了倍周期，第11組的周期性顯然沒有第9組明顯?？偟膩碚f，周期性的最佳分析結(jié)果為第9組。下面用具有多個(gè)觀測(cè)數(shù)據(jù)的兩個(gè)源驗(yàn)證此結(jié)果。

圖4 功率譜密度分析方法(PSD)對(duì)sin函數(shù)不同數(shù)據(jù)點(diǎn)的周期性分析
Fig.4 Periodic analysis of sin functions′ different data points by using PSD method

圖5 sin函數(shù)用時(shí)間補(bǔ)償離散傅里葉變換分析方法 (DCDFT)、Jurkevich方法、離散相關(guān)分析方法 (DCF) 和功率譜密度分析方法 (PSD) 進(jìn)行周期分析時(shí)要求的最低數(shù)據(jù)點(diǎn) (右邊的兩個(gè)圖形從上向下分別為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函數(shù)的V2m-P，左圖分組數(shù)為1～6組，右圖分組數(shù)為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函數(shù)的V2m-P，左圖分組數(shù)為1～6組，右圖分組數(shù)為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函數(shù)的V2m-P，左圖分組數(shù)為1～6組，右圖分組數(shù)為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 計(jì)算并分析類星體3C 279和3C 454.3的光變周期

3.1 類星體3C 279的數(shù)據(jù)點(diǎn)及光變周期

本文研究的類星體3C 279在B波段和3C 454.3在B波段數(shù)據(jù)點(diǎn)主要是從網(wǎng)站(http://www.astro.yale.edu/smarts/glast/home.php)獲取。

從2008年2月5日到2017年7月19日，由以上網(wǎng)站搜尋到類星體3C 279在B波段的數(shù)據(jù)點(diǎn)有661個(gè)，B波段的光變曲線如圖9，可以看出有幾次大的爆發(fā)，最亮?xí)rB波段星等為14.9，光變曲線中有6個(gè)以上周期存在，這些均滿足Jurkevich方法中確認(rèn)長(zhǎng)周期存在的必要條件[13]。

利用Jurkevich方法分組數(shù)m=9，圖10可得到3C 279在B波段的光變周期(2.81 ± 0.54)年。周期Pn與V2m的關(guān)系曲線中有很多的極小值，為了有效地判斷周期Pn與V2m關(guān)系曲線中所得出周期的真實(shí)性，文[14]給出了較好的判據(jù)：

(34)

其中，V2m為歸一化的值。在歸一化的V2m-Pn圖中，若V2m=1.0，則有f=0，即在樣本數(shù)據(jù)中沒有表現(xiàn)出周期性；若V2m=0，則有f=1，這時(shí)樣本中存在的最大周期能夠從圖中識(shí)別。在通常情況下，當(dāng)f≥ 0.5時(shí)，樣本數(shù)據(jù)表現(xiàn)出非常強(qiáng)的周期性；當(dāng)f≤ 0.25時(shí)，則表現(xiàn)為不含有周期性。進(jìn)一步分析V2m值，選擇曲線較平滑部分的最小深度和噪聲做實(shí)驗(yàn)，如果平滑部分的相對(duì)最小爆發(fā)值比平滑部分的V2m大10倍，則可以進(jìn)一步討論相應(yīng)的短周期性。在計(jì)算中分組數(shù)m的大小十分重要，組數(shù)m分得越多，靈敏度越高，但每組數(shù)據(jù)中有少數(shù)的數(shù)據(jù)點(diǎn)在圖中產(chǎn)生較大的噪聲。分組數(shù)m較大時(shí)，還增加V2m的計(jì)算量。反之，有可能尋找不到周期。如圖10，對(duì)類星體3C 279在B波段的光變周期分析中，很容易找到V2m及P值，并用文[14]的判據(jù)：f=(1-V2m)/V2m進(jìn)行檢驗(yàn)，結(jié)果如表1。

通過對(duì)表1分析發(fā)現(xiàn)，最小值V2m=0.352 1，f=1.840 1，類星體3C 279可能的光變周期為(2.81 ± 0.54)年。利用Jurkevich方法，分組數(shù)m=9，分析了類星體3C 279在B波段的光變數(shù)據(jù)，圖10顯示了B波段的分析結(jié)果。從表1可以看出，在滿足判據(jù)f的條件下[14]，類星體3C 279在B波段的P2，P3與P1之間存在類星體簡(jiǎn)單的倍數(shù)關(guān)系：P2≈2P1，P3≈3P1，說明它們之間可能存在天文學(xué)倍頻關(guān)系。

圖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的數(shù)據(jù)點(diǎn)及光變周期

3C 454.3(PKS 2251+158, OY091)是一個(gè)比較亮、變化比較劇烈的類星體，并且在光學(xué)和射電波段存在比較明顯的相關(guān)性[15]。本文的數(shù)據(jù)是從2008年6月23日到2017年7月30日，約765個(gè)數(shù)據(jù)點(diǎn)，獲得了如圖11的歷史光變曲線。

利用Jurkevich方法，分組數(shù)m=9，分析了類星體3C 454.3在B波段的光變數(shù)據(jù)，圖12顯示了B波段的分析結(jié)果。從圖12可以看出，類星體3C 454.3在B波段有3個(gè)明顯的極小值，意味著類星體3C 454.3在B波段光變曲線中存在3個(gè)可能的周期，它們分別是：P1=457 d，P2=891 d和P3=1 320 d。根據(jù)文[16]的結(jié)果，要確定一個(gè)周期，數(shù)據(jù)樣本的時(shí)間跨度要超過周期的6倍。文中數(shù)據(jù)樣本的時(shí)間跨度大約2 500 d，小于P2，P3的6倍， 故P2，P3必須排除，并需要更多的觀測(cè)數(shù)據(jù)確定它們。但是P2，P3與P1之間存在著簡(jiǎn)單的倍數(shù)關(guān)系：P2≈2P1，P3≈3P1，說明它們之間可能存在天文學(xué)倍頻關(guān)系，意味著類星體3C 454.3在B波段的光變曲線中可能存在一個(gè)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，對(duì)類星體3C 454.3在B波段的光變周期分析中，很容易找到V2m及P值，并用文[14]的判據(jù)：f=(1-V2m)/V2m進(jìn)行檢驗(yàn)，結(jié)果如表2。

通過對(duì)表2分析發(fā)現(xiàn)，當(dāng)V2m取最小值0.550 5，f=0.817，對(duì)應(yīng)的光變周期可能為P=457 d，這既滿足了圖12的分析結(jié)果，也滿足f判據(jù)[14]。

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

4 討論與結(jié)論

本文利用一個(gè)以2π為周期的正弦函數(shù)y=sinθ檢驗(yàn)4種研究方法的可靠性。在實(shí)驗(yàn)中，使用的單位為弧度，選擇0 rad為起點(diǎn)，步長(zhǎng)為0.1 rad，不同的數(shù)據(jù)點(diǎn)個(gè)數(shù)共取15組，得到Jurkevich方法的分析結(jié)果為(6.3 ± 0.017) rad；時(shí)間補(bǔ)償離散傅里葉變換分析方法的分析結(jié)果為(6.331 ± 0.130) rad；離散相關(guān)分析方法的分析結(jié)果為(6.287 ± 0.088) rad；功率譜密度方法的分析結(jié)果為(6.287 ± 0.045) rad。分析結(jié)果表明，Jurkevich方法、時(shí)間補(bǔ)償離散傅里葉變換分析方法、離散相關(guān)分析方法和功率譜密度方法周期需要的數(shù)據(jù)點(diǎn)最低分別為50個(gè)、120個(gè)、100個(gè)和60個(gè)。獲取最短的連續(xù)數(shù)據(jù)采樣后，Jurkevich方法最有效。Jurkevich方法分析結(jié)果在4種方法中最精確可靠，且此計(jì)算方法簡(jiǎn)捷實(shí)用。

在實(shí)際應(yīng)用中，發(fā)現(xiàn)Jurkevich方法十分依賴于分組數(shù)m，m越大分析結(jié)果越好，但會(huì)產(chǎn)生比較大的噪聲；反之則有可能找不到周期[17]。如果數(shù)據(jù)分布不均勻，則會(huì)導(dǎo)致組間數(shù)據(jù)分布偏差很大，從而影響獲得實(shí)際的周期。根據(jù)圖6、圖7和圖8可以看出，對(duì)于不同的m值，V2m與Pn的關(guān)系曲線的傾斜程度不同，周期大的向下傾斜。結(jié)合f檢驗(yàn)公式可知，V2m越小，f值越大，更容易滿足f檢驗(yàn)。由此可知f檢驗(yàn)依賴于m值的大小。利用模擬數(shù)據(jù)尋找到Jurkevich方法的分組數(shù)為m=9時(shí)，分析結(jié)果最佳。最后用分組數(shù)m=9時(shí)的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倍。