2. 粒子技术与辐射成像教育部 重点实验室 (清华大学), 北京 100084
2. Key Laboratory of Particle and Radiation Imaging (Tsinghua University) Ministry of Education, Beijing 100084, China
天体发射源周围存在反射介质,因此观测得到的信号中会有反射成分。利用相应的方法,分析反射成分可以获得超新星遗迹的前身星、活动星系核(Active Galactic Nuclei,AGN)核心黑洞的质量等重要信息。
反射信号与发射信号之间有一定的时间延迟,利用时延可以研究反射介质的相对位置、速度等参数及相关信息。例如,可以利用回光现象确定超新星前身星。超新星爆发时周围若存在大量的尘埃,其爆发时发出的光线抵达尘埃所在的地方,尘埃会散射光线,同时尘埃会吸收超新星的光线,并以不同的波长辐射,即是超新星的回光现象[1]。回光到达地球的时间较爆发的时间有延迟,但可以保存爆发时与光线有关的信息,利用回光分析可确定超新星遗迹周围星际介质的分布、超新星遗迹的年龄及距离等信息,并且可以获得超新星爆发时的光谱以确定超新星的类型[2]。通过观测回光,文[3]证实了Cassiopeia A Supernova是来自红超巨星氦核坍缩的IIb型超新星。观测到第谷超新星的436年后,Krause等人观测到其回光,并获得回光图像与光谱,确定了第谷超新星为正常Ia型超新星[4]。由于回光现象较超新星爆发信号有延迟,因此利用回光现象可以进行距离的测定及光谱分析,对超新星的诸多细节进行研究,例如确定前身星、测量尘埃环的半径等[5]。
信号处理中,利用相关函数可以有效地检测带有噪声的信号,以及信号的周期性,测量信号的时延长度。在提取反射信号与发射信号间的延迟信息时,相关函数是常用的方法[6]。
天体物理中,在测量活动星系核核心黑洞质量时所用到的反响映射(Reverberation Mapping)原理,就是利用信号的相关函数检测时延长度,通过时延长度确定宽线区(Broad-Line Region,BLR)与核心的距离,利用这一距离可以确定活动星系核核心黑洞的质量[7]。
在活动星系核统一模型中,核心有由吸积驱动的超大质量黑洞,吸积盘的外侧是宽线区,核心黑洞的连续辐射传播到宽线区,由于光致电离效应,此处的尘埃会受激产生发射线,同时,宽线区的尘埃会绕黑洞旋转,因此发射线会有多普勒展宽,利用维里定理和宽线区的谱线展宽,可以计算核心黑洞的维里质量[8, 9, 10]。由于光由核心黑洞传播到宽线区有时间延迟,同时光致电离的特征时间较传播时延很短,因此宽线区的半径可以由这一时延与光速的乘积描述[11]。通过计算连续光变曲线和发射线光变曲线的互相关函数(Cross-Correlation Function,CCF),可以确定平均时延长度与相关系数的关系,可以将使得相关系数值最大时的时延作为线谱与连续谱的时延,实际求解时一般通过求解互相关函数的质量中心(Centroid)确定平均时延长度[12]。反响映射的主要工作就是利用转移方程及观测量,确定平均时延长度。具体的求解过程利用了相关原理,通过相关函数可以确定发射线与连续谱的时延。
发射天体与反射天体之间的距离过近,同时都离地球极其遥远时,将不能区分发射信号与反射信号,此时接收信号是发射信号和反射信号的叠加。当反射体与发射源相对静止时,发射信号与反射信号之间的延迟固定不变,可以直接求解接收信号的自相关,且会出现两个相关峰值,峰值的间隔即给出了延迟大小。而当反射体与发射源之间有相对运动时,反射信号和发射信号的延迟将成为一个随时间变化的函数,此时对接收的信号求解自相关将得不到第2个相关峰,即不能直接给出延迟信息。本文提出了一种基于对接收信号在时域做变换的方法进行求解,可以有效地在有相对运动时求解延迟信息,具体的做法在基本原理部分进行说明。
双星系统是由两个围绕共同质心运动的恒星组成的恒星系统,较亮的一颗恒星被称作主星,另一颗被称作伴星。当双星系统中包含了致密天体,如白矮星、中子星或黑洞时,来自另一颗恒星(donor)的气体会在致密天体周围吸积(accretion),产生大量的辐射。致密天体是白矮星时,被称作激变变星(cataclysmic variable)[13]。致密天体是中子星或者黑洞时,即是X射线双星。X射线双星会有更高的亮度,本文提出的方法可以在X射线双星系统得到更好的应用。
本文包含3部分内容,第1部分介绍了利用相关原理提取延迟信息的方法与缺陷,并提出了一种可以在发射天体与反射介质间有相对运动时提取延迟信息的方法,在此基础上给出了对双星系统轨道的估算方法;第2部分利用该方法,进行了MATLAB的模拟,验证了这种方法的可行性;第3部分是对该方法以及模拟结果的总结。
1 基本原理利用相关原理可以提取时延信息,对于发射源与反射介质没有相对运动的情形,可以直接利用接收信号的自相关函数提取反射信号的时间延迟,从而给出反射介质的位置信息。当发射源与反射介质之间存在相对运动时,反射信号的时间延迟是时间的函数,此时已经不能直接用相关原理处理。双星系统沿椭圆轨道运动,因此不能直接使用自相关函数计算时延。利用参数空间搜索的方法,实现信号在时域的变换,可以计算反射信号与发射信号的时延,进行双星系统轨道的估算。
1.1 相关原理这里研究的发射信号为有限带宽随机信号,其自相关信号是有一定展宽的脉冲函数。反射体与发射源相对静止,即信号传到接收者的时间差为常数。还需要做如下几点具体的说明。
首先,该发射信号具有一定的相干长度。由于处理信号的主要手段是进行相关运算,得到较好的相关结果的基本要求即是待处理信号具有一定的相干长度。因此只进行对具有良好的相干长度的信号进行处理。
其次,发射信号要足够强。由于被观测的双星系统距离地球很远,导致信号衰减幅度很大,同时信号反射时会大幅衰减,若发射信号不够强,可能导致接收的信号中反射信号成分过于微弱,无法提取所需的信息。
同时,由于被观测的双星系统与地球相距甚远,且伴星和主星间的距离较小,故可以认为接收的发射信号和反射信号是平行光。
记接收的信号为R(t),发射信号为S(t),在双星系统中,反射介质是伴星,记伴星的反射系数为κ,反射信号与发射信号的时间延迟为d,
反射信号可以表示为
SRef(t)=κS(t-d),
(1)
接收信号可以表示为
R(t)=S(t)+κS(t-d),
(2)
发射信号的自相关函数为
CR(τ)=(1+κ2)CS(τ)+κCS(τ+d)+κCS(τ-d).
(5)
|
| 图 1 发射信号自相关函数。发射信号是模拟产生的有限带宽伪随机信号,信号长度为10 s,相干长度为200 ms,采样间隔为1 ms Fig. 1 The Auto-Correlation Function (ACF) of a sequence of transmitted signals. Here,the transmitted signals are simulated pseudo random signals of a limited bandwidth. The signal sequence has a total length of 10s,a coherence length of 200ms,and a sampling interval of 1ms |
如图 2,时间延迟为d=1 500 s。可以看出除τ=0 s处有明显峰值外,在τ=1 500 s处也出现了明显的峰值,此处峰值与τ=0处主峰值之比的理论值为
,κ为伴星的反射系数。在数值模拟中,由于信号是有限长度的,因此实际值会小于
。
|
| 图 2 接收信号自相关函数。发射信号是模拟的有限带宽伪随机信号,信号长度为4 000 s,相干长度为100 ms,采样间隔为10 ms,反射信号时间延迟为1 500 s,反射系数为0.5,由于τ<0没有物理含义,只画出τ>0的部分 Fig. 2 The Auto-Correlation Function (ACF) of a sequence of received signals. Here,the received signals are from transmitted signals simulated to be pseudo random signals of a limited bandwidth. The transmitted-signal sequence has a total length of 4000s,a coherence length of 100ms,and a sampling interval of 10ms. The time delay of the reflected signals is 1500s,and the reflection coefficient is 0.5. We only show the ACF for τ>0 because τ<0 has no physical meaning |
存在相对运动时,反射信号的时间延迟是时间的函数,记作d(t),则反射信号为
SRef(t)=κS[t-d(t)],
(6)
R(t)=S(t)+κS[t-d(t)],
(7)
|
| 图 3 接收信号自相关函数。发射信号是模拟的有限带宽伪随机信号,长度为2 000 s,相干长度为100 ms,采样间隔为50 ms,反射体作初始距离为600 ls,速度为0.06 c的匀速直线运动,反射系数为0.5,由于τ<0没有物理含义,只画出τ>0的部分 Fig. 3 The Auto-Correlation Function (ACF) of a sequence of received signals. Here,the received signals are from transmitted signals simulated to be pseudo random signals of a limited bandwidth. The transmitted-signal sequence has a total length of 2000s,a coherence length of 100ms,and a sampling interval of 50ms. A reflector is assumed here to be moving in a constant velocity of a speed 0.06c. The initial distance of the reflector from the signal emitter is 600ls. The reflector has a reflection coefficient 0.5. We only show the ACF for τ>0 because τ<0 has no physical meaning |
记y=f(t)=t-d(t),接收信号可改写为R(t)=S(t)+κS[f(t)],其中d(t)由双星的轨道及运动信息确定,可以给出待定参数的具体表达式,这些参数是不随时间变化的,可以完全描述双星的轨道和运动。
f(t)是非线性的,可以对其局部线性化,并把斜率归一化,也就是说可将y=f(t)映射成为yC=t+t0,其中t0是与时间t无关的常数。这里给出一种具体的做法。如图 4,将f(t)细分成若干段,每一段可以认为是线性的,每段端点的坐标可以容易给出,由此可以计算每一段的两点式,对(ti,yi)点与(ti+1,yi+1)点之间的线段两点式为
,令
,代入两点式,即可得yC=y,对每一小段进行相同的操作即可得到需要的映射结果,这样映射得到的直线没有常数项t0,方便后续的相关计算。
|
| 图 4 时域映射示意图。将y=f(t)分割成若干可以看作线性的小段,将(ti,yi)点与(ti+1,yi+1)点之间的线段局部放大画在图中。写出其两点式,即可以利用给出的映射方式得到映射后的直线yC=t Fig. 4 Illustration of the map of time to time delay. The map divides time into small intervals in each of which the map function y=f(t) is approximately linear. The map within an arbitrary time interval is shown in the plot after being magnified,with the end points of the interval and corresponding map values marked as (ti,yi) and (ti+1,yi+1). The line yC=t can be determined using the linear-approximation formula of the mapping method |
在对R(t)的时域做了如上映射后,可以得到新的信号:
RC(t)=S[D(t)]+S(t+t0),
(9)
接收信号R(t)与映射后信号RC(t)的互相关函数为
CRRC(τ)=κCS(τ-t0).
(12)
可以看出,互相关函数会在τ=t0处出现峰值,使用本文的修正方法,t0=0,因此在τ=0处出现峰值。
做上述时域映射需要给出f(t)=t-d(t)的除时间t外不含其它变量的表达式。在双星系统中,d(t)由轨道参数与运动参数唯一确定,但这些参数的具体值需要求解。可以给出这些参数的大致取值范围,遍历此范围内的值,每一个值对应的f(t)进行时域映射操作后计算接收信号R(t)与映射后信号RC(t)的互相关,当互相关有明显的峰值,并且峰值达到最大,即相关性最好时,可以认为此组参数是所要求解的轨道参数与运动参数。
因此进行参数空间搜索,对每一个参数组做时域映射并计算接收信号与映射后信号的互相关,找出互相关峰值最大时对应的轨道参数和运动参数,即是需要估算的双星系统的轨道参数与运动参数。
2 双星系统轨道估算的数值模拟利用时域映射的方法,用模拟产生的信号进行数值模拟,搜索模拟接收数据中的反射信号,验证该方法。
模拟时采用的发射信号为有限带宽伪随机信号,在生成接收信号时,由于t时刻的反射信号是由该时刻之前对应的发射信号给出的,因此需要截取生成的发射信号中间的一段信号进行模拟。
针对X射线双星系统进行模拟,考虑反射体与发射天体、接收者在同一平面内的二维情形,且伴星围绕主星以椭圆轨道运行,该椭圆轨道的长轴和发射天体与观测者的连线垂直。建立以主星所在的椭圆焦点为极点的极坐标系,极径记作r,极角记作θ。这种情形可以看做是一个Edge On的双星系统。
描述一个三维空间双星系统的轨道参数(Orbital Elements)有:倾角(Orbital inclination)、升交点黄经(Ecliptic longitude of the ascending node)、近星点经度(Ecliptic longitude of the periastron)、半长轴(Semimajor axis)、离心率(Eccentricity)、近星点时刻(Time at the periastron)[14]。在本文中,考虑二维平面内的情形,可以使用更少的轨道参数描述双星运动。
模拟伴星围绕主星做椭圆轨道运动的情形,选择利用离心率、半长轴、起始观测时的极角作为描述伴星轨道形状的轨道参数,伴星单位时间内扫过的面积是常数,可以用这一个参数作为描述伴星的运动的参数。所选择的4个参数可以完全描述伴星的轨道和运动。
记该椭圆模型的离心率为e;半长轴为a;起始观测时极角记作θ0;伴星单位时间内扫过的面积是
;观测时间长度为t;反射系数为κ。
模拟时,如下设定该X射线双星系统的参数。伴星的反射系数为0.1,其运动和轨道采取如下的参数:e=0.01,a=600 ls,θ0=0 rad,
=0.04 ls2/s。
生成长度为40 000 s的序列作为发射信号,由于存在时延,计算t时刻反射信号时需要之前的发射信号,因此计算接收信号时,只将后20 000 s的序列作为被接收的发射信号,同时被接收的反射信号可以利用整个40 000 s的序列通过计算时延给出。
参数e的搜索范围为e∈[0.001,0.02],搜索步长为0.001;参数a的搜索范围为e∈[590 ls,610 ls],搜索步长为0.5 ls;参数θ0的搜索范围为θ0∈[-1 rad,1 rad],搜索步长为0.01 rad;参数
的搜索范围为
∈[-0.04 ls2/s,1.2 ls2/s],搜索步长为0.000 1 ls2/s。
用上述参数在全空间进行搜索,计算每组参数对应的接收信号R(t)和映射后信号RC(t)的相关函数,记录各组参数对应的相关函数的峰值,所有峰值的最大值对应的一组参数为:e=0.01,a=600 ls,θ0=0 rad,
=0.04 ls2/s。与模拟给定的参数吻合得很好。
利用搜索的结果,下面给出相关函数峰值随各个参数变化的曲线。
将其他参数固定在搜索求解出的数值,在给出的搜索范围改变θ0的值,给出相关函数峰值随参数θ0变化的曲线,在θ0=0处出现了明显的峰值。θ0表示观测初始时伴星的位置,图 5说明,通过搜索,很好地确定了伴星初始位置。
|
| 图 5 相关峰值随θ0的变化。将相关峰值归一化,使得相关峰值的最大值为1 Fig. 5 The variations of normalized peak values with θ0 values |
将其他参数固定在搜索求解出的数值,在给出的搜索范围改变
的值,给出相关函数峰值随参
数
变化的曲线,在
=0.04 ls2/s处出现了明显的峰值。
可以表征运动的速度,图 6说明,通过搜索,很好地确定了伴星的运动速度。
|
图 6 相关峰值随 的变化。将相关峰值归一化,使得相关峰值的最大值为1
Fig. 6 The variations of normalized peak values with values
|
将其他参数固定在搜索求解出的数值,同时在给出的搜索范围改变参数e、a,画出相关函数峰值随参数e、a变化,在e=0.01、a=600 ls处出现了明显的峰值。参数e、a可以确定伴星的椭圆轨道形状,图 7说明,通过搜索,很好地确定了伴星轨道。
|
| 图 7 相关峰值随e、a的变化。将相关峰值归一化,使得相关峰值的最大值为1 Fig. 7 The variations of normalized peak values with e and a values |
通过搜索,可以确定参数e、a、θ0、
,这4个参数可以完全确定伴星的运动、位置和轨道,因此,该方法可以实现双星系统的轨道估算,并且可以给出更丰富的信息。
提取反射信号的时延,可以用于确定双星系统的轨道和运动信息。反射信号混叠在发射信号中,同时由于伴星与主星之间的相对运动而相位发生了很大的变化,直接计算相关函数无法实现时延的提取。通过对接收信号进行时域映射的方法,可以有效地在有相对运动时提取反射信号,将信号在时域进行映射处理,消除了由相对运动带来的影响,建立椭圆轨道的运动模型,使用e、a、θ0、
4个参数进行描述,利用相关原理,可以准确地搜索反射信号,之后,通过数值模拟的方法验证了本方法的可行性,在具体情形下准确地找到了对应的反射信号,实现了双星系统的轨道估算。
在实际的X射线双星系统中,伴星反射信号的强度与两星的距离D及伴星的半径RDS相关,即反射系数
。如果κ过小,反射信号将过于微弱,很不利于轨道的估算。在模拟的过程中发现,如果κ小于某一极限值,程序不能有效地估算轨道,这一极限值与采样的间隔和观测时间的长度有关,采样间隔越短,观测时间越长,对应的极限值越小,在本文模拟的情形下,κ<0.01时将难以进行估算。因此,在实际的X射线双星系统估算中,大的κ有利于轨道的估算。X射线双星系统吸积方式主要有两种,一种是星风吸积(Wind-fed accretion),另一种是洛希瓣溢出(Roche-lobe overflow)。在大质量X射线双星系统(High-Mass X-ray Binary,HMXB)中,伴星质量大,星风强烈,星风被致密星体吸积,星风吸积占主导,此时的伴星半径小,距离远,反射系数κ小。在低质量X射线双星系统(Low-Mass X-ray Binary,LMXB)中,伴星质量小,星风弱,伴星离开主序向巨星和超巨型演化时会充满洛希瓣(Roche Lobe),吸积的主导方式为洛希瓣溢出,此时的伴星半径大,距离近,反射系数κ会很大。对于洛希瓣吸积的X射线双星系统GRS1915+105[15],距离D=107.7 ± 9.9R⊙,RDS=RRoche=24 ± 4R⊙,对应的κ∝0.05。因此,在洛希瓣吸积时的轨道估算确实可能会更为容易。
| [1] | 田文武, 杨雪娟. 当前超新星遗迹研究中的若干热点问题[J]. 天文学进展, 2010, 28(2): 97-111. Tian Wenwu, Yang Xuejuan. Light echoes and remnants of historically galactic supernovae[J]. Progress in Astronomy, 2010, 28(2): 97-111. |
| [2] | Rest A, Suntzeff N B, Olsen K, et al. Light echoes from ancient supernovae in the Large Magellanic Cloud[J]. Nature, 2005, 438(7071): 1132-1134. |
| [3] | Krause O, Birkmann S M, Usuda T, et al. The Cassiopeia A supernova was of type IIb[J]. Science, 2008, 320(5880): 1195-1197. |
| [4] | Krause O, Tanaka M, Usuda T, et al. Tycho Brahe's 1572 supernova as a standard Type Ia as revealed by its light-echo spectrum[J]. Nature, 2008, 456(7222): 617-619. |
| [5] | Gouiffes C, Rosa M, Melnick J, et al. Light echoes from SN 1987 A[J]. Astronomy and Astrophysics, 1988, 198(1-2): L9-L12. |
| [6] | 葛新成, 罗大成, 曹勇. 相关函数在数字信号处理中的应用[J]. 电光与控制, 2006, 13(6): 78-80. Ge Xincheng, Luo Dacheng, Cao Yong, Application of correlation function in digital signal processing[J]. Electronics Optics & Control, 2006, 13(6): 78-80. |
| [7] | Peterson B M. Reverberation mapping of active galactic nuclei[J]. Publications of the Astronomical Society of the Pacific, 1993, 105(685): 247-268. |
| [8] | Netzer H, Peterson B M. Reverberation mapping and the physics of active galactic nuclei[M]//Maoz D, Sternberg A, Leibowitz E M. Astronomical Time Series: Astrophysics and Space Science Library. 1997: 85-108. |
| [9] | Peterson B M, Ferrarese L, Gilbert K M, et al. Central masses and broad-line region sizes of active galactic nuclei. II. A homogeneous analysis of a large reverberation-mapping database[J]. The Astrophysical Journal, 2004, 613(2): 682-699. |
| [10] | Peterson B M. The central black hole and relationships with the host galaxy[J]. New Astronomy Reviews, 2008, 52(6): 240-252. |
| [11] | Peterson B M, Horne K. Reverberation mapping of active galactic nuclei[C]//Planets to Cosmology: Essential Science in the Final Years of the Hubble Space Telescope. 2006: 89-98. |
| [12] | Peterson B M. Variability of active galactic nuclei[C]//Advanced Lectures on the Starburst-AGN Connection, Proceedings of a conference held in Tonantzintla. 2001: 3-67. |
| [13] | Smith R C. Cataclysmic variables[J]. Contemporary Physics, 2006, 47(6): 363-386. |
| [14] | Benacquista M. An introduction to the evolution of single and binary stars[M]. New York: Springer, 2012: 13-28. |
| [15] | Greiner J. The binary components of GRS 1915+105[C]//X-ray Binaries to Gamma-Ray Bursts: Jan van Paradijs Memorial Symposium. 2003: 111-120. |


