蔡劍華 劉備 楊江河， 庹滿先 聶建軍 汪勝輝 樊軍輝

(1 湖南文理學(xué)院數(shù)理學(xué)院 常德 415000)

(2 廣州大學(xué)天體物理中心 廣州 510006)

1 引言

耀變體是活動(dòng)星系核中觀測(cè)性質(zhì)極其特殊的子類(lèi)， 而光變是它們最突出的觀測(cè)特征， 大多數(shù)表現(xiàn)為長(zhǎng)期的無(wú)規(guī)則快速光變， 有的光變表現(xiàn)出周期性[1-7]. Blazar的光變表現(xiàn)出光譜隨亮度的變化， 光變時(shí)標(biāo)不但帶來(lái)了輻射區(qū)大小信息， 甚至帶來(lái)了Blazar噴流強(qiáng)弱的信息[8-16]. 周期性光變可能是Blazar內(nèi)部結(jié)構(gòu)的反映[8]， 但是由于觀測(cè)數(shù)據(jù)分布的不均勻性， 使得從光變曲線中分析周期性會(huì)有一定的不確定性. 在長(zhǎng)期的研究中， 人們尋找天體光變周期的常用方法有:周期圖譜法[17-18]、相位分析法[19]、相關(guān)函數(shù)方法[20-21]、Jurkevich方法[22]、模型譜估計(jì)法[23-25]、小波分析法[26-28]和Hilbert-Huang變換[29]等. 后來(lái)， 為了更好地處理天體觀測(cè)中的非等間隔數(shù)據(jù)， 人們又提出了CLEAN方法[30]、CLEANEST方法[31-32]和基于小波的向量投影方法[32]. 這些方法在天文觀測(cè)和理論研究中得到了廣泛應(yīng)用， 并有作者不斷將最新發(fā)展起來(lái)的現(xiàn)代信號(hào)處理方法應(yīng)用到天體物理領(lǐng)域. 3C 273是最早發(fā)現(xiàn)的類(lèi)星體， 是研究最多的天體之一. 對(duì)類(lèi)星體3C 273光變周期的分析， 已有學(xué)者開(kāi)展了相關(guān)的工作. Smith等[33]研究表明3C 273在1887-1963年的光變曲線中有5-6個(gè)不同周期的波動(dòng)， 并稱(chēng)其有12.7-15.2 yr的周期(P). Babadzhanyants等[34]分析了3C 273在1887-1991年的B波段光變曲線， 發(fā)現(xiàn)P= 13.4 yr; 2001年樊軍輝等[35-36]研究了3C 273的110 yr光變數(shù)據(jù)，發(fā)現(xiàn)了其P=2.0、(13.65±0.20)和(22.5±0.2)yr的周期.2007年Fan等[37]還研究了它的射電光變周期性， 發(fā)現(xiàn)其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]根據(jù)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]綜合自己以前的工作和對(duì)1998-2008年觀測(cè)數(shù)據(jù)的分析， 進(jìn)一步得出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]基于密歇根大學(xué)射電天文臺(tái)數(shù)據(jù)庫(kù)(Michigan Radio Astronomy Observatory， UMRAO)從1965年到2012年所收集的3C 273天體47 yr的8.0 GHz射電波段觀測(cè)數(shù)據(jù)， 得到3C 273的平均周期是16.9、7.5和2.5 yr.

本文試圖在前人工作的基礎(chǔ)上， 結(jié)合近期發(fā)展起來(lái)的非平穩(wěn)信號(hào)處理方法-經(jīng)驗(yàn)?zāi)B(tài)分解和傳統(tǒng)自回歸(Auto-regressive，AR)模型譜估計(jì)的優(yōu)點(diǎn)，提出一種組合的功率譜估計(jì)方法， 并將他們應(yīng)用于類(lèi)星體3C 273的光變周期估算， 旨在為天體的周期光變研究探索一條新途徑.

2 方法原理

2.1 經(jīng)驗(yàn)?zāi)B(tài)分解

經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition， EMD)， 又叫Huang變換， 它是由美籍華人科學(xué)家N. E. Huang于1998年提出的一種新型信號(hào)分解方法. 與被廣泛使用的小波分解不同， 經(jīng)驗(yàn)?zāi)B(tài)分解無(wú)須選擇基函數(shù)， 特別適合于對(duì)非平穩(wěn)、非線性信號(hào)的處理. 經(jīng)驗(yàn)?zāi)B(tài)分解采用“篩”選的方法， 將信號(hào)x(t)分解為一組固有模態(tài)函數(shù)(Intrinsic Mode Function， IMF)的和[41-42]，t表示時(shí)間.

其中n表示分解得到的最大階數(shù)，i表示取其中的各階，Ci(t)即為各階模態(tài)函數(shù)(為作圖和描述方便常簡(jiǎn)為Ci)，R(t)為剩余分量， 具體分解過(guò)程見(jiàn)文獻(xiàn)[41]. 固有模態(tài)函數(shù)滿足信號(hào)關(guān)于時(shí)間軸局部對(duì)稱(chēng)和零點(diǎn)數(shù)與極點(diǎn)數(shù)相等或至多相差1的兩個(gè)條件， 包含了原信號(hào)不同時(shí)間尺度的局部特征信息.EMD分解實(shí)質(zhì)上是對(duì)非平穩(wěn)數(shù)據(jù)進(jìn)行平穩(wěn)化處理的一個(gè)過(guò)程， 在空間域表現(xiàn)為從小尺度到大尺度的層層濾波， 在頻域則表現(xiàn)為從高頻到低頻的層層濾波. 同時(shí)， Huang等[41-42]也證明了經(jīng)驗(yàn)?zāi)B(tài)分解的完備性， 即分解得到的各固有模態(tài)以及剩余分量相加能夠恢復(fù)原始信號(hào). 這些特征為信號(hào)的相關(guān)性分析和構(gòu)造濾波器抑制噪聲提供了條件[43-44].

若要消除高頻部分的噪聲， 一個(gè)包含多個(gè)IMF的低通時(shí)空濾波器可構(gòu)造為

k為要保留的IMF的最小階數(shù). 若噪聲僅存在低頻部分， 高通時(shí)空濾波器為

h為要保留的IMF的最大階數(shù). 若噪聲在高低頻都存在， 則帶通時(shí)空濾波器為

本文就是在EMD分解后得到各階分量， 并計(jì)算其與原始光變數(shù)據(jù)的相關(guān)度， 以此構(gòu)造濾波器，取相關(guān)系數(shù)大的系列IMF分量來(lái)進(jìn)行進(jìn)一步譜估計(jì)的.

2.2 AR譜估計(jì)方法

隨機(jī)過(guò)程s(m)的AR模型可以表示為[45-46]

m為離散數(shù)據(jù)點(diǎn)，al(l=1，2，···，p)為AR模型的系數(shù)，l為各系數(shù)序列標(biāo)號(hào)，p為系數(shù)序列標(biāo)號(hào)的最大值.u(m)是一個(gè)方差為δ2的白噪聲序列. 求解AR模型系數(shù)常用的方法有Yule-Walker法、Burg法和最小二乘法， 本文采用的是Yule-Walker法. AR模型常用的階數(shù)判定準(zhǔn)則有3種: FPE (Final Predietion Error)準(zhǔn)則、AIC (Akaike Information Criterion)準(zhǔn)則和BIC (Bayesian Information Criterion)準(zhǔn)則， 本文采用的是AIC準(zhǔn)則. 隨機(jī)過(guò)程的功率譜估計(jì)定義為[47-48]

j為虛數(shù)單位，ω為角變量. 以上(5)式和(6)式給出的模型被稱(chēng)為自回歸模型， 簡(jiǎn)稱(chēng)AR模型， 獲得的功率譜估計(jì)稱(chēng)為AR譜. 求得s(m)的前p+1個(gè)自相關(guān)函數(shù)， 代入公式即可求得s(m)的功率譜.

2.3 方法步驟

基于EMD-AR譜的天體光變周期分析流程如圖1所示， 描述如下:

圖1 基于EMD-AR譜的天體光變周期分析流程Fig.1 Analysis flow of variability periodicity based on EMD-AR spectrum

(1)首先對(duì)觀測(cè)數(shù)據(jù)進(jìn)行EMD分解， 得到各階模態(tài)函數(shù)“C1”-“Cn”;

(2)計(jì)算各階IMF分量與原始觀測(cè)數(shù)據(jù)的相關(guān)度， 去除相關(guān)系數(shù)小的分量， 用余下的IMF分量“Ch”-“Ck”重構(gòu)數(shù)據(jù)，k和h的取值由(2)-(4)式確定;

(3)最后估算其AR譜， 進(jìn)一步分析天體的光變周期.

3 數(shù)據(jù)樣本

本文分析的數(shù)據(jù)來(lái)自廣州大學(xué)樊軍輝教授課題組在其系列論文[16， 18， 35-37， 39， 49]中收集和觀測(cè)的類(lèi)星體3C 273光變數(shù)據(jù)， 時(shí)間跨度為1887-2016年， 長(zhǎng)達(dá)130 yr， 共有4342個(gè)觀測(cè)數(shù)據(jù). 由于歷史原因和天體觀測(cè)的特點(diǎn)， 有些月份甚至年份的觀測(cè)數(shù)據(jù)缺失. 但在后期， 尤其是在2000年后， 數(shù)據(jù)觀測(cè)記錄的頻次很高， 出現(xiàn)了多段連續(xù)幾個(gè)小時(shí)觀測(cè)， 且以1 min為時(shí)間間隔記錄的觀測(cè)數(shù)據(jù)， 這些數(shù)據(jù)為3C 273短時(shí)標(biāo)周期的分析提供了條件.

用于長(zhǎng)時(shí)標(biāo)周期分析的數(shù)據(jù)如下. 在已有的1887-2016年的共4342個(gè)觀測(cè)數(shù)據(jù)中， 根據(jù)時(shí)標(biāo)每年隨機(jī)提取出一個(gè)數(shù)據(jù)， 共112個(gè)觀測(cè)數(shù)據(jù)， 為便于分析我們將此設(shè)為L(zhǎng)1光變曲線. 在4342個(gè)觀測(cè)數(shù)據(jù)中， 按時(shí)標(biāo)每月隨機(jī)提取出一個(gè)數(shù)據(jù)， 共469個(gè)觀測(cè)數(shù)據(jù)， 設(shè)為L(zhǎng)2光變曲線. 同上， 在觀測(cè)數(shù)據(jù)中，按時(shí)標(biāo)隨機(jī)每天提取出一個(gè)數(shù)據(jù)， 共1254個(gè)觀測(cè)數(shù)據(jù)， 設(shè)為L(zhǎng)3光變曲線. L1-L3光變曲線的時(shí)域波形如圖2 (a)-(c)所示.

圖2 類(lèi)星體3C 273的觀測(cè)數(shù)據(jù): (a)-(c)用于長(zhǎng)時(shí)標(biāo)周期計(jì)算; (d)-(f)用于短時(shí)標(biāo)周期計(jì)算.Fig.2 Observational data of quasar 3C 273: (a)-(c) for long time scale period calculation， (d)-(f) for short time scale period calculation.

用于短時(shí)標(biāo)周期分析的數(shù)據(jù)如下. 分析4342個(gè)觀測(cè)數(shù)據(jù)， 發(fā)現(xiàn)有很多段較連續(xù)的以“min”為時(shí)間間隔的觀測(cè)數(shù)據(jù)， 為我們進(jìn)行天內(nèi)光變的研究提供了條件. 選取較長(zhǎng)的3段: Julian day (JD)2454459.4701-JD2454459.6516數(shù)據(jù)段中， 在連續(xù)的4 h內(nèi)有192個(gè)觀測(cè)數(shù)據(jù)， 我們將此設(shè)為G1光變曲線; JD2454460.4704-JD2454460.6459數(shù)據(jù)段中， 在連續(xù)的4 h內(nèi)有198個(gè)觀測(cè)數(shù)據(jù)， 設(shè)為G2光變曲線; JD2454461.4713-JD2454461.6165數(shù)據(jù)段中，在連續(xù)的3 h內(nèi)有165個(gè)觀測(cè)數(shù)據(jù)， 設(shè)為G3光變曲線. G1-G3光變曲線的時(shí)域波形如圖2 (d)-(f)所示. 這3段觀測(cè)數(shù)據(jù)連續(xù)性很好， 幾乎是等間隔的以“min”為單位的采樣觀測(cè)， 分別對(duì)它們進(jìn)行周期分析， 嘗試得到以“min”為單位的3C 273短時(shí)標(biāo)周期.

4 仿真分析

為驗(yàn)證采用本文方法進(jìn)行功率譜估計(jì)的有效性， 作者首先對(duì)一組仿真信號(hào)進(jìn)行了處理. 設(shè)計(jì)兩個(gè)頻率為30 Hz和70 Hz正弦分量的疊加信號(hào)并對(duì)其采樣， 采樣頻率為1000 Hz、采樣點(diǎn)數(shù)為1459以構(gòu)成仿真信號(hào). 采用本文提出的EMD-AR譜估計(jì)方法， 對(duì)其進(jìn)行功率譜的估計(jì).

圖3給出了仿真信號(hào)EMD分解后的各IMF分量及其對(duì)應(yīng)的功率譜. 從圖3可以清楚地看到EMD分解在空間域表現(xiàn)為從小尺度到大尺度， 在頻域則表現(xiàn)為從高頻到低頻的分解過(guò)程. 分析圖3 (a)“C1”分量和“C2”分量的時(shí)域波形圖及其圖3 (b)對(duì)應(yīng)的頻譜， 30 Hz和70 Hz的2個(gè)頻率成分已被精確地分解出來(lái). 計(jì)算“C1”-“R”分量與原始信號(hào)的相關(guān)度分別為: 0.7548、0.7075、0.0006、0.0025、0.0023、0.0016、0.0011. 顯然， “C1”-“C2”分量與原始數(shù)據(jù)的相關(guān)系數(shù)比“C3”-“R”分量的相關(guān)系數(shù)大得多， 他們不在同一數(shù)量級(jí). 選擇相關(guān)系數(shù)大的“C1”-“C2”分量重構(gòu)數(shù)據(jù)， 并進(jìn)行AR譜估計(jì). 圖4所示為采用本文提出的EMD-AR譜方法得到的功率譜估計(jì). 可以看出， 功率譜的精度很高，30 Hz和70 Hz的頻率成分已被精確地分離和表征出來(lái).

圖3 仿真信號(hào)EMD分解后的各IMF分量及其對(duì)應(yīng)的功率譜Fig.3 IMF components and its power spectrum of simulation signal after EMD decomposition

圖4 EMD-AR譜方法估算的仿真信號(hào)的功率譜Fig.4 Power spectrum of simulation signal estimated by EMD-AR method

為進(jìn)一步驗(yàn)證本文方法可有效過(guò)濾非周期性信號(hào)且不會(huì)得到虛假的周期信號(hào)， 作者在仿真上加入了一些強(qiáng)度隨機(jī)、寬度隨機(jī)的“Bumps”脈沖噪聲信號(hào)， 再進(jìn)一步作譜估計(jì). 如圖5 (a)所示， “星劃”線為原始的兩個(gè)正弦函數(shù)疊加的仿真信號(hào)， “點(diǎn)劃”線為加入的“Bumps”脈沖噪聲. 加噪后的仿真信號(hào)如圖5 (b)所示. 依照本文提出的方法， 先對(duì)信號(hào)進(jìn)行EMD分解，得到的各級(jí)IMF分量如圖6(a)所示. 分別計(jì)算他們與原始仿真信號(hào)的相關(guān)系數(shù)為:0.7255、0.6946、0.0006、0.0006、0.0001、0.0119、0.0039. 取相關(guān)系數(shù)大的“C1”-“C2”分量重構(gòu)數(shù)據(jù)，并進(jìn)一步作AR譜估計(jì)， 得到的功率譜如圖6 (b)所示. 顯然， 在圖6 (b)中30 Hz和70 Hz的頻率成分也非常清晰， 本文方法能正確地將周期信息從含噪信號(hào)中辨識(shí)、提取出來(lái)， 可有效過(guò)濾這些非周期性信號(hào)且不會(huì)得到虛假的周期信號(hào)， 為正確估算光變周期提供了條件.

圖5 加噪前后的仿真信號(hào): (a)原始仿真信號(hào)和“Bumps”脈沖噪聲; (b)加入“Bumps”脈沖噪聲后的仿真信號(hào).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 加噪仿真信號(hào)的EMD分解(a)和本文方法估算的加噪仿真信號(hào)的功率譜(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 長(zhǎng)時(shí)標(biāo)周期計(jì)算

L1光變曲線是以年為單位的數(shù)據(jù)樣本， 在1887-2016年的觀測(cè)數(shù)據(jù)中， 每年隨機(jī)保留1個(gè)數(shù)據(jù). 對(duì)提取出的數(shù)據(jù)進(jìn)行三次樣條插值預(yù)處理(本文三次樣條插值采用的邊界條件是非結(jié)點(diǎn)邊界條件， 由于本數(shù)據(jù)曲線兩端數(shù)據(jù)差的絕對(duì)值也較小， 非結(jié)點(diǎn)邊界條件有利于對(duì)周期的估算， 下同.)， 然后對(duì)其進(jìn)行EMD分解， 自適應(yīng)地得到6階模態(tài)函數(shù)“C1”-“R”， 如圖7 (a)所示. 計(jì)算他們與原始信號(hào)的相關(guān)系數(shù)分別為: 0.5455、0.5611、0.3500、0.2371、0.0232、0.2338. 去除相關(guān)度很小的“C5”分量，用其余分量重構(gòu)信號(hào)， 進(jìn)一步估算AR譜如圖7 (b)所示.從圖7的功率譜中， 可以計(jì)算出L1光變曲線的光變周期為11.02和4.69 yr.

圖7 L1光變曲線的EMD分解(a)及其光變周期估計(jì)(b)Fig.7 EMD decomposition of L1 light curve (a) and estimation of variability periodicity (b)

類(lèi)比L1的處理， 對(duì)L2光變曲線進(jìn)行EMD分解后得到8階模態(tài)函數(shù)“C1”-“R”， 如圖8 (a)所示. 分別計(jì)算他們與原始信號(hào)的相關(guān)系數(shù)為: 0.4979、0.3705、0.4546、0.4411、0.2935、0.0509、0.1208、0.1593. 去除相關(guān)度很小的“C6”分量， 用其余分量重構(gòu)信號(hào)， 再進(jìn)一步估算AR譜如圖8 (b)所示. 從圖8的功率譜中， 計(jì)算得到L2光變曲線的光變周期為13.51、2.76和1.46 yr.

圖8 L2光變曲線的EMD分解(a)及其光變周期估計(jì)(b)Fig.8 EMD decomposition of L2 light curve (a) and estimation of variability periodicity (b)

同上處理方法， 對(duì)L3光變曲線進(jìn)行EMD分解后得到11階模態(tài)函數(shù)“C1”-“R”， 如圖9 (a)所示. 計(jì)算他們與原始信號(hào)的相關(guān)系數(shù)為: 0.0404、0.0493、0.0822、0.1269、0.1341、0.3268、0.4112、0.4100、0.2474、0.0939、0.4152. 去除相關(guān)度很小的“C1”、“C2”、“C3”和“C10”分量， 用其余分量重構(gòu)信號(hào)，再進(jìn)一步估算AR譜如圖9 (b)所示. 從圖9的功率譜中， 可以計(jì)算出L3光變曲線的光變周期為: 21.23、11.02、5.51、4.69、3.79、2.04 yr.

圖9 L3光變曲線的EMD分解(a)及其光變周期估計(jì)(b)Fig.9 EMD decomposition of L3 light curve (a) and estimation of variability periodicity (b)

5.2 短時(shí)標(biāo)周期計(jì)算

G1段數(shù)據(jù)為4 h內(nèi)連續(xù)記錄的192個(gè)觀測(cè)數(shù)據(jù)，采用三次樣條插值的方法補(bǔ)充完整數(shù)據(jù)為240采樣點(diǎn)， 然后對(duì)其進(jìn)行EMD分解. 如圖10 (a)所示， 光變曲線被自適應(yīng)地分解， 得到7階模態(tài)函數(shù)“C1”-“R”，計(jì)算他們與原始信號(hào)的相關(guān)系數(shù)分別為: 0.5864、0.5097、0.4326、0.2836、0.3491、0.0515、0.3796.去除相關(guān)度很小的“C6”分量， 用其余分量重構(gòu)信號(hào)， 再進(jìn)一步估算AR譜如圖10 (b)所示. 從圖10的功率譜中， 可以計(jì)算出G1光變曲線的周期為: 240、30、15.3、10、7.3、6.4、4.2 min.

圖10 G1光變曲線的EMD分解(a)及其光變周期估計(jì)(b)Fig.10 EMD decomposition of G1 light curve (a) and estimation of variability periodicity (b)

類(lèi)比G1段數(shù)據(jù)的處理， 對(duì)G2光變曲線進(jìn)行EMD分解后得到8階模態(tài)函數(shù)“C1”-“R”，如圖11(a)所示. 分別計(jì)算他們與原始信號(hào)的相關(guān)系數(shù)為:0.7559、0.4842、0.2828， 0.2951、0.1245、0.0329、0.1922、0.1800. 去除相關(guān)度很小的“C6”分量，用其余分量重構(gòu)信號(hào)， 再進(jìn)一步估算AR譜如圖11 (b)所示. 從圖11的功率譜中， 可以計(jì)算出G2光變曲線的周期為: 240、30、14.9、7.5、6.2、4.4、3.5、3.02 min.

圖11 G2光變曲線的EMD分解(a)及其光變周期估計(jì)(b)Fig.11 EMD decomposition of G2 light curve (a) and estimation of variability periodicity (b)

同上處理方法， 將G3段數(shù)據(jù)進(jìn)行EMD分解后得到7階模態(tài)函數(shù)“C1”-“R”， 如圖12 (a)所示. 計(jì)算他們與原始信號(hào)的相關(guān)系數(shù)分別為: 0.6786、0.3648、0.3090、0.3922、0.1689、0.2518、0.1924.7個(gè)相關(guān)度系數(shù)都在同一個(gè)數(shù)量級(jí)， 所以保留所有分量重構(gòu)信號(hào)， 再進(jìn)一步估算AR譜如圖12 (b)所示. 從圖12的功率譜中， 可計(jì)算出G3光變曲線的周期為: 180、29.1、15.1、10、7.2、5、4.4、3.5、3.03 min.

圖12 G3光變曲線的EMD分解(a)及其光變周期估計(jì)(b)Fig.12 EMD decomposition of G3 light curve (a) and estimation of variability periodicity (b)

5.3 結(jié)果與討論

為便于分析與討論， 將3條長(zhǎng)時(shí)標(biāo)周期光變曲線估算的光變周期列表， 并與一些文獻(xiàn)中得到的類(lèi)星體3C 273的光變周期作對(duì)比， 如表1所示(P1-P6表示6個(gè)估算的周期). 本文得到的21.23 yr的周期與Fan等[39]中給出(21.10±0.14) yr的結(jié)論一致. 本文得到的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的周期呈較好的倍周期關(guān)系. 本文得到的4.69 yr、2.76 yr的周期與Vol’vach等[38]給出的(4.9±0.3) yr和(2.8±0.3)yr的周期一致. 綜合上述分析結(jié)果， 本文得到的類(lèi)星體3C 273的長(zhǎng)時(shí)標(biāo)周期為: 21.23、11.02、5.51、13.51、4.69、3.79、2.76 yr.值得一提的是，因?yàn)樵紨?shù)據(jù)跨距了130 yr， 由于歷史原因， 早期數(shù)據(jù)的可靠性無(wú)法保證， 我們只用近半個(gè)世紀(jì)的數(shù)據(jù)重復(fù)做了分析， 得到的結(jié)果是一致的.

表1 不同方法估算的類(lèi)星體3C 273的長(zhǎng)時(shí)標(biāo)周期Table 1 Long time scale period of quasar 3C 273 estimated by different methods

同樣的， 為便于分析， 將3段短時(shí)標(biāo)周期數(shù)據(jù)估算的光變周期列表， 如表2所示(P1-P9表示9個(gè)估算的周期). G1和G2光變曲線計(jì)算所得的240 min周期及G3光變曲線所得的180 min周期， 分別與他們的數(shù)據(jù)等長(zhǎng)， 不能作為周期判定. 3條光變曲線分別在30、15和7.5 min， 10和5 min以及6和3 min附近都出現(xiàn)了峰值， 且表現(xiàn)為倍周期關(guān)系. 綜合3段光變曲線的分析， 可預(yù)測(cè)類(lèi)星體3C 273的短時(shí)標(biāo)周期為: (30±1)、(15±0.3)、(7.5±0.2)、(10±0.1)、(5±0.6)和(6±0.4)、(3±0.5) min.

表2 不同光變曲線估算的類(lèi)星體3C 273的短時(shí)標(biāo)周期Table 2 Short time period of quasar 3C 273 estimated by different light curves

6 結(jié)論

本文利用類(lèi)星體3C 273長(zhǎng)達(dá)130 yr的觀測(cè)數(shù)據(jù)， 提取得到年、月和日的光變數(shù)據(jù)， 作為長(zhǎng)時(shí)標(biāo)周期估算的光變曲線; 同時(shí)選取3段連續(xù)性很好， 幾乎是等間隔的以“min”為單位的觀測(cè)數(shù)據(jù)， 作為對(duì)該星體進(jìn)行短時(shí)標(biāo)周期估算的數(shù)據(jù)樣本. 采用最新發(fā)展起來(lái)的經(jīng)驗(yàn)?zāi)B(tài)分解與傳統(tǒng)的AR模型譜估計(jì)的組合方法， 將具有非線性、非平穩(wěn)特性的光變數(shù)據(jù)分解為具有單分量信號(hào)特征的模態(tài)函數(shù)， 再選取相關(guān)系數(shù)高的分量重構(gòu)信號(hào)后估算AR譜， 以此計(jì)算光變曲線的周期. 通過(guò)對(duì)類(lèi)星體3C 273的光變資料分析， 我們得到類(lèi)星體3C 273的長(zhǎng)時(shí)標(biāo)周期為:21.23、13.51、11.02、5.51、4.69、3.79、2.76 yr，與此前文獻(xiàn)[33-34， 38-40]計(jì)算的結(jié)果一致; 短時(shí)標(biāo)周期為: (30±1)、(15±0.3)、(7.5±0.2)、(10±0.1)、(5±0.6)和(6±0.4)、(3±0.5) min， 這些可能的周期還存在3組大致倍周期的關(guān)系. 分析結(jié)果也表明， 基于EMD-AR譜的方法對(duì)計(jì)算Blazar天體的光變周期是可靠的， 為天體的周期光變研究提供了一條新思路.