基于EMD-AR譜的3C 273光變周期分析*

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

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

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

1 引言

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

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

2 方法原理

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

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

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

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

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

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

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

2.2 AR譜估計方法

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

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

j為虛數(shù)單位，ω為角變量. 以上(5)式和(6)式給出的模型被稱為自回歸模型， 簡稱AR模型， 獲得的功率譜估計稱為AR譜. 求得s(m)的前p+1個自相關(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)首先對觀測數(shù)據(jù)進行EMD分解， 得到各階模態(tài)函數(shù)“C1”-“Cn”;

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

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

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

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

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

圖2 類星體3C 273的觀測數(shù)據(jù): (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.

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

4 仿真分析

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

圖3給出了仿真信號EMD分解后的各IMF分量及其對應(yīng)的功率譜. 從圖3可以清楚地看到EMD分解在空間域表現(xiàn)為從小尺度到大尺度， 在頻域則表現(xiàn)為從高頻到低頻的分解過程. 分析圖3 (a)“C1”分量和“C2”分量的時域波形圖及其圖3 (b)對應(yīng)的頻譜， 30 Hz和70 Hz的2個頻率成分已被精確地分解出來. 計算“C1”-“R”分量與原始信號的相關(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ù)量級. 選擇相關(guān)系數(shù)大的“C1”-“C2”分量重構(gòu)數(shù)據(jù)， 并進行AR譜估計. 圖4所示為采用本文提出的EMD-AR譜方法得到的功率譜估計. 可以看出， 功率譜的精度很高，30 Hz和70 Hz的頻率成分已被精確地分離和表征出來.

圖3 仿真信號EMD分解后的各IMF分量及其對應(yīng)的功率譜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)所示， “星劃”線為原始的兩個正弦函數(shù)疊加的仿真信號， “點劃”線為加入的“Bumps”脈沖噪聲. 加噪后的仿真信號如圖5 (b)所示. 依照本文提出的方法， 先對信號進行EMD分解，得到的各級IMF分量如圖6(a)所示. 分別計算他們與原始仿真信號的相關(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ù)，并進一步作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光變曲線是以年為單位的數(shù)據(jù)樣本， 在1887-2016年的觀測數(shù)據(jù)中， 每年隨機保留1個數(shù)據(jù). 對提取出的數(shù)據(jù)進行三次樣條插值預(yù)處理(本文三次樣條插值采用的邊界條件是非結(jié)點邊界條件， 由于本數(shù)據(jù)曲線兩端數(shù)據(jù)差的絕對值也較小， 非結(jié)點邊界條件有利于對周期的估算， 下同.)， 然后對其進行EMD分解， 自適應(yīng)地得到6階模態(tài)函數(shù)“C1”-“R”， 如圖7 (a)所示. 計算他們與原始信號的相關(guān)系數(shù)分別為: 0.5455、0.5611、0.3500、0.2371、0.0232、0.2338. 去除相關(guān)度很小的“C5”分量，用其余分量重構(gòu)信號， 進一步估算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階模態(tài)函數(shù)“C1”-“R”， 如圖8 (a)所示. 分別計算他們與原始信號的相關(guān)系數(shù)為: 0.4979、0.3705、0.4546、0.4411、0.2935、0.0509、0.1208、0.1593. 去除相關(guān)度很小的“C6”分量， 用其余分量重構(gòu)信號， 再進一步估算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階模態(tài)函數(shù)“C1”-“R”， 如圖9 (a)所示. 計算他們與原始信號的相關(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)信號，再進一步估算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段數(shù)據(jù)為4 h內(nèi)連續(xù)記錄的192個觀測數(shù)據(jù)，采用三次樣條插值的方法補充完整數(shù)據(jù)為240采樣點， 然后對其進行EMD分解. 如圖10 (a)所示， 光變曲線被自適應(yīng)地分解， 得到7階模態(tài)函數(shù)“C1”-“R”，計算他們與原始信號的相關(guān)系數(shù)分別為: 0.5864、0.5097、0.4326、0.2836、0.3491、0.0515、0.3796.去除相關(guān)度很小的“C6”分量， 用其余分量重構(gòu)信號， 再進一步估算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段數(shù)據(jù)的處理， 對G2光變曲線進行EMD分解后得到8階模態(tài)函數(shù)“C1”-“R”，如圖11(a)所示. 分別計算他們與原始信號的相關(guān)系數(shù)為:0.7559、0.4842、0.2828， 0.2951、0.1245、0.0329、0.1922、0.1800. 去除相關(guān)度很小的“C6”分量，用其余分量重構(gòu)信號， 再進一步估算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段數(shù)據(jù)進行EMD分解后得到7階模態(tài)函數(shù)“C1”-“R”， 如圖12 (a)所示. 計算他們與原始信號的相關(guān)系數(shù)分別為: 0.6786、0.3648、0.3090、0.3922、0.1689、0.2518、0.1924.7個相關(guān)度系數(shù)都在同一個數(shù)量級， 所以保留所有分量重構(gòu)信號， 再進一步估算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 結(jié)果與討論

為便于分析與討論， 將3條長時標周期光變曲線估算的光變周期列表， 并與一些文獻中得到的類星體3C 273的光變周期作對比， 如表1所示(P1-P6表示6個估算的周期). 本文得到的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é)果， 本文得到的類星體3C 273的長時標周期為: 21.23、11.02、5.51、13.51、4.69、3.79、2.76 yr.值得一提的是，因為原始數(shù)據(jù)跨距了130 yr， 由于歷史原因， 早期數(shù)據(jù)的可靠性無法保證， 我們只用近半個世紀的數(shù)據(jù)重復(fù)做了分析， 得到的結(jié)果是一致的.

表1 不同方法估算的類星體3C 273的長時標周期Table 1 Long time scale period of quasar 3C 273 estimated by different methods

同樣的， 為便于分析， 將3段短時標周期數(shù)據(jù)估算的光變周期列表， 如表2所示(P1-P9表示9個估算的周期). G1和G2光變曲線計算所得的240 min周期及G3光變曲線所得的180 min周期， 分別與他們的數(shù)據(jù)等長， 不能作為周期判定. 3條光變曲線分別在30、15和7.5 min， 10和5 min以及6和3 min附近都出現(xiàn)了峰值， 且表現(xiàn)為倍周期關(guān)系. 綜合3段光變曲線的分析， 可預(yù)測類星體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 結(jié)論

本文利用類星體3C 273長達130 yr的觀測數(shù)據(jù)， 提取得到年、月和日的光變數(shù)據(jù)， 作為長時標周期估算的光變曲線; 同時選取3段連續(xù)性很好， 幾乎是等間隔的以“min”為單位的觀測數(shù)據(jù)， 作為對該星體進行短時標周期估算的數(shù)據(jù)樣本. 采用最新發(fā)展起來的經(jīng)驗?zāi)B(tài)分解與傳統(tǒng)的AR模型譜估計的組合方法， 將具有非線性、非平穩(wěn)特性的光變數(shù)據(jù)分解為具有單分量信號特征的模態(tài)函數(shù)， 再選取相關(guān)系數(shù)高的分量重構(gòu)信號后估算AR譜， 以此計算光變曲線的周期. 通過對類星體3C 273的光變資料分析， 我們得到類星體3C 273的長時標周期為:21.23、13.51、11.02、5.51、4.69、3.79、2.76 yr，與此前文獻[33-34， 38-40]計算的結(jié)果一致; 短時標周期為: (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譜的方法對計算Blazar天體的光變周期是可靠的， 為天體的周期光變研究提供了一條新思路.