基于样条函数的恒星光谱自动归一化方法
罗锋1,2, 刘超1,2, 赵永恒1,2     
1. 中国科学院国家天文台, 北京 100101;
2. 中国科学院大学, 北京 100049
摘要: 恒星的观测谱一般由连续谱、谱线和噪声组成,其中连续谱是黑体辐射导致的辐射流量随波长变化的光滑连续光谱。光谱分类及恒星物理参数估计等研究依赖于连续谱及谱线信息的准确提取。因此光谱数据处理的工作主要是拟合连续谱,并通过对光谱进行归一化来提取谱线特征。连续谱拟合的方法主要有多项式拟合、中值滤波、小波滤波等。已有的方法在低信噪比、宇宙线信号干扰、存在发射线等情况下,有不同程度的局限性,体现在鲁棒性和准确度上。目前,针对郭守敬望远镜的107条光谱没有自动化方法应用到归一化上的问题,研究并开发一种适用于不同的温度、信噪比及波长覆盖范围,并能够自动化处理的恒星光谱归一化方法,显得十分迫切。在仔细分析不同类型光谱的基础上,提出了一种基于固定窗口划分的连续谱拟合方法。该方法对光谱中能够体现连续谱特征的数据点进行筛选提取,通过细微地控制样条函数平滑度产生更加准确的连续谱。使用郭守敬望远镜中不同光谱型、温度范围、波长覆盖范围的光谱进行实验,结果表明,该方法具有良好的精度和普适性。
关键词: 连续谱归一化    郭守敬望远镜    恒星光谱    样条函数    
Automatic Normalization Method of Stellar Spectrum Based on Spline Function
Luo Feng1,2, Liu Chao1,2, Zhao Yongheng1,2     
1. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: The observed spectrum of stars is generally composed of continuous spectrum, spectral line and noise. The continuous spectrum is the smooth continuous spectrum of radiation flux varying with wavelength caused by blackbody radiation. Spectral classification and stellar physical parameter estimation depend on the accurate extraction of continuous spectrum and spectral line information. Therefore, the main work of spectral data processing is to fit continuous spectrum and extract spectral line features by normalization. At present, the methods of continuous spectrum fitting mainly include polynomial fitting, median filtering and wavelet filtering. Existing methods have limitations to varying degrees in the case of low SNR, interference of cosmic ray signals and existence of emission lines, which are mainly reflected in robustness and accuracy. For the time being, there is no automated method applying to the normalization of the 107 spectra from LAMOST. In the period of avalanche of astronomical data, it is very urgent to research and develop a spectral normalization algorithm of stars with better universality and automatic processing that can be applied to a wider range of temperature, SNR and wavelength coverage. On the basis of careful analysis of different types of spectra, a continuous spectrum fitting method based on fixed window dividing is proposed. This method can filter and extract the data points in the spectrum which can reflect the characteristics of continuous spectrum, and produce more accurate continuous spectrum by fine controlling the smoothness of spline function. Experiments are carried out using spectra of different spectral types, temperature ranges and wavelength coverage ranges in LAMOST, and the results show that the proposed algorithm has good accuracy and universality.
Key words: Continuous spectrum normalization    LAMOST    Stellar spectrum    Spline function    

作为国家重大科学工程,大天区面积多目标光纤光谱天文望远镜(the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST, 又叫郭守敬望远镜)突破了天文望远镜大视场与大口径难以兼得的难题[1],成为目前国际上口径最大的大视场望远镜,是我国光学望远镜研制的又一里程碑。截至2018年6月,郭守敬望远镜完成了7年的巡天观测,获取的光谱数首次超过千万量级,成为世界上第一个获取光谱数超千万的巡天项目。遗憾的是,对于这107条光谱,目前没有自动化方法应用到归一化上,本文主要对这方面进行研究。

每条观测谱都由连续谱、谱线和噪声组成[2]。通过谱线可以对恒星的化学组成进行分析,并且由于多普勒效应,谱线还携带着视向速度的信息[3]。因此,从原始光谱中准确地提取谱线是恒星光谱研究及光谱数据处理工作的重要步骤,对后续的研究有重要意义。

要消除连续谱首先要对连续谱进行拟合,目前连续谱拟合方法主要有多项式拟合、中值滤波、小波滤波、形态滤波器等[4]。此外,利用天文软件[5]进行半自动化处理的方法是人工筛选提取认为合适的数据点[6]并选择函数形式,然后由软件进行多项式拟合得到连续谱。这种方法依赖人工参与,无法满足目前巡天项目中海量光谱的实时处理需求。国外的斯隆数字巡天(Sloan Digital Sky Survey, SDSS)中的数据处理程序SSPP(Segue Stellar Parameter Pipeline)采用分段多项式拟合的方法[7]:首先将光谱分为红端和蓝端两部分,去除强巴尔末线系后,分别对蓝端使用9阶、红端使用4阶多项式分别拟合,丢弃3个标准差以外的数据点;然后拼接两段连续谱再进行一次9阶多项式拟合,得到最终的连续谱。这种方法稍显繁琐,且计算量大。

本文提出的方法是将光谱划分为数个固定的(包含相同的像素数量)窗口,窗口内选中值点,通过筛选策略丢弃部分参考点,这样做能够有效地避开发射线、吸收线及宇宙线的干扰。通过对参考点之间“夹角”的控制,使用非常小的平滑度构造样条曲线,实现对连续谱的精确拟合。

1 方法介绍

本文提出的连续谱归一化方法包含3个步骤:(1)将一条待处理光谱根据数据点总数划分为数个窗口,每个窗口选出流量中值点形成参考点集;(2)制定一些筛选策略,将参考点集中有可能影响拟合结果的数据点丢弃;(3)根据参考点流量最大值所在波长位置选取对应的样条函数平滑度经验值进行拟合,得到连续谱。

1.1 划分窗口构造参考点集

为了有足够的数据点作为拟合函数的控制点,根据输入的待处理光谱数据点总数,采用如表 1的策略划分窗口。设数据点总数为n,每个窗口包含像素数,即窗宽为w

表 1 窗宽与数据点总数之间的对应关系 Table 1 The relationship between windows width and the total number of data points
数据点总数n 0~1 000 1 000~1 500 1 500~2 500 2 500~3 500 3 500~4 500 4 500~n
窗宽w 13 15 17 19 21 23

根据表 1确定窗宽后对光谱进行窗口划分,遍历所有窗口,每个窗口内选取流量值等于窗口中所有点中值的像素作为参考点,并将其加入参考点集Srp

1.2 筛选参考点

1.1节中构造的参考点集Srp只体现了光谱的大致形状,其中有些点可能位于吸收线或者发射线上,如图 1图 2

图 1 有些参考点位于吸收线上 Fig. 1 Some reference points are on the absorption line
图 2 有些参考点位于发射线上 Fig. 2 Some reference points are on the emission line

理论上,连续谱应非常接近黑体辐射谱,而黑体谱可以划分为两个单调区间,峰值左边是单调增区间,右边是单调减区间[8],希望用来构造拟合曲线的参考点符合并体现这一特征。为了实现这一目标,采取两个步骤对参考点集再次筛选。

(1) 丢弃局部极小值点

Srp中一个参考点的索引为i,流量为f(i),若i满足:

$ f\left( {i - 1} \right) > f\left( i \right)\;并且\;f\left( i \right) < f\left( {i + 1} \right),i \ge 1, $ (1)

则将参考点iSrp中删去。每次从i=1开始到遍历完Srp中所有参考点结束为一趟,在一趟丢弃结束时,有些点可能成为新的局部极小值点,此时需要再丢弃一趟。设丢弃的趟数为t,需要保证以下两点:(1)对连续谱正常的光谱,t次丢弃之内Srp中元素数量不再减少;(2)对连续谱异常[9]的光谱,t次丢弃后Srp中有足够的参考点用来构造拟合曲线。实验表明,合适的t取值为t=len(Srp)//30 + 2。

经这一步骤处理后,图 1图 2中两条光谱的参考点集如图 3图 4

图 3 图 1中光谱经t次丢弃局部极小值后的参考点集 Fig. 3 The reference point set after the local minimum is discarded t times in the spectrum in Fig. 1
图 4 图 2中光谱经t 次丢弃局部极小值后的参考点集 Fig. 4 The reference point set after the local minimum is discarded t times in the spectrum in Fig. 2

(2) 丢弃“夹角”较大的点

为了获得更加精确的连续谱拟合曲线,需要较小的平滑度构造样条函数。如果有的参考点和其左右两点连线构成的夹角过大,则会给拟合曲线带来较大的偏离,这样的参考点应从Srp中删去,如图 5

图 5 位于发射线上的数据点应从参考点集中删去 Fig. 5 The data point on the emission line should be removed from the reference point set

图 6,在参考点集中取一点i,设其对应的波长为w(i),流量为f(i),与它左右两点间连线的夹角为θ(i),组成的三角形边长分别为a, b, c。在几何上,有:

$ a = \sqrt {{{\left[ {w\left( i \right) - w\left( {i - 1} \right)} \right]}^2} + {{\left[ {f\left( i \right) - f\left( {i - 1} \right)} \right]}^2}} , $ (2)
$ b = \sqrt {{{\left[ {w\left( {i + 1} \right) - w\left( i \right)} \right]}^2} + {{\left[ {f\left( {i + 1} \right) - f\left( i \right)} \right]}^2}} , $ (3)
$ c = \sqrt {{{\left[ {w\left( {i + 1} \right) - w\left( {i - 1} \right)} \right]}^2} + {{\left[ {f\left( {i + 1} \right) - f\left( {i - 1} \right)} \right]}^2}} . $ (4)
图 6 使用余弦定理计算“夹角” Fig. 6 Calculate the angle by use of Cosine theorem

根据余弦定理计算θ(i)的值。实际上,流量值并不对应长度单位,wf的单位很难统一。由于光谱数据中流量是未定标的相对流量,直接用余弦定理计算的夹角随着光谱流量值所在的区间不同及光谱数据处理中必要的乘除运算而变化。为此做如下处理:

对于Srp中的一个参考点i,设其与左边数据点波长数值的差为wd-,右边为wd+,则:

$ w{d^ - } = w\left( i \right) - w\left( {i - 1} \right), $ (5)
$ w{d^ + } = w\left( {i + 1} \right) - w\left( i \right). $ (6)

同理,对于流量值[2]

$ f{d^ - } = f\left( i \right) - f\left( {i - 1} \right), $ (7)
$ f{d^ + } = f\left( {i + 1} \right) - f\left( i \right). $ (8)

引入比例因子rw, rf,令:

$ {r_w} = \frac{{w{d^ - }}}{{w{d^ - } + w{d^ + }}}, $ (9)
$ {r_{\rm{f}}} = \frac{{f{d^ - }}}{{f{d^ - } + f{d^ + }}}. $ (10)

计算判据r=|rw-rf|的值并引入阈值r0。若r > r0,则认为点i属于夹角较大的点,会给精确拟合带来干扰,应从Srp中除去。通过实验发现,效果较好的r0取值为0.5。

通过比值消去量纲,则不论wf取何种单位,rw, rf都是不变的。使用rw, rf作为判断依据寻找这样的参考点更为可靠。图 5中光谱经此步骤处理后的结果如图 7

图 7 位于发射线上的参考点被成功移除 Fig. 7 The reference point on the emission line was successfully removed
1.3 平滑度计算

样条函数[10]拟合效果的好坏,除了与参考点的选取有关,还取决于平滑度的选择。对于同一条光谱,同样的参考点集,不同的平滑度会有不同的效果,如图 8

图 8 平滑度对拟合曲线的影响 Fig. 8 The influence of smoothing factor on fitting curve

较小的平滑度可以更好地照顾到光谱图像上起伏的细节,但过小容易产生振动频繁的现象,使拟合曲线无法体现连续谱;较大的平滑度不会有剧烈的振动,但是过大会使拟合曲线过于平缓,也无法很好地体现连续谱。应该指出的是,对于同一个平滑度,在不同的流量区间,样条曲线会有不同的表现。在拟合之前,将光谱数据进行最大值标准化,即flux=flux/max(flux)。

经实验分析,温度越高的光谱需要越小的平滑度照顾光谱在蓝端较急剧的转折,温度越低的光谱需要越大的平滑度弱化参考点上下起伏的趋势。而一条光谱的温度在连续谱上体现为流量峰值所在的位置。采用在参考点集中选取流量最大值,用其所在波长位置估计该条光谱的温度区间,选择对应的平滑度s0。为了能够适用于局部归一化的场合,对于参考点较少的情况,需要更小的平滑度。为此,引入默认值为1的比例系数c,用于在参考点较少时将之与s0相乘的结果作为最终的平滑度取值。具体数据列在表 2表 3中。

表 2 根据参考点集中流量最大值点所在波长位置选择平滑度 Table 2 Selection of smoothing factor according to the wavelength position of the maximum flux point in the reference point set
wave/nm 370~476 476~688 688~794 794~900
s0 0.01 0.02 0.03 0.05
表 3 根据参考点集中数据点数n选择平滑度比例系数c Table 3 Selection of proportionality coefficient c of smoothing factor according to the number of data points n in the reference point set
n 1~29 30~59 60~89 90~119 120~149 150~180
c 0.2 0.35 0.5 0.65 0.8 1

在构造样条曲线时使用的平滑度为:

$ s = c * {s_0}. $ (11)
2 实验结果 2.1 不同情况下的处理结果

为检验所提出方法的处理效果,从郭守敬望远镜中随机抽取了不同类型的光谱10 000条进行实验。从中选出不同温度的光谱进行拟合,如图 9。图中左侧蓝色实线为原始光谱,红色实线为拟合曲线;右侧绿色实线为对应的归一化谱。由图 9可见,本方法在不同温度时做到了较好的拟合。

图 9 不同温度下的归一化结果 Fig. 9 Normalization results at different temperatures

图 10,左侧桃红色实线是原始光谱,绿色实线为样条函数拟合曲线;右侧靛蓝色实线为对应的归一化谱。选取不同的波长覆盖范围后,本方法能够自动应用于需要局部归一化的场合,不需要人工干预。

图 10 不同波长覆盖范围的归一化结果 Fig. 10 Normalization results of different wavelength coverage range

图 11,左侧是原始光谱及拟合曲线,在不同的信噪比情况下,本方法拟合效果受到的影响不大。与文[11]中专为低质量光谱开发的归一化方法处理效果非常接近。

图 11 不同信噪比的归一化结果 Fig. 11 Normalization results under different SNR

图 12,在光谱中含有发射线及宇宙线的情况下,本方法仍能有效识别并进行准确的拟合。

图 12 含有发射线或宇宙线情况下的归一化结果 Fig. 12 Circumstances containing emission lines or cosmic rays
2.2 误差分析和适用范围

由于仪器原因,光谱中常见的仪器形变模式有3种,如图 13 (a)。为了测试本文提出的方法是否能够有效地处理这些仪器导致的光谱形变,在有效温度4 000~8 000 K的范围内随机抽取9条光谱,每条光谱归一化后的结果都乘以这3种形变曲线,并再次归一化。图 13 (b)的直方图中显示了两次归一化结果残差的分布,可以看出本文算法的误差平均值在10-3量级,误差弥散在10-2量级。

图 13 郭守敬望远镜的形变模式及归一化误差分析 Fig. 13 Deformation mode in LAMOST and error analysis of normalization

经实验分析,本方法适用条件如下:信噪比范围5~600,波长覆盖范围370~900 nm,有效温度范围3 000~50 000 K。其中,信噪比低于5的光谱处理结果不稳定,总体精度有随着信噪比下降而降低的趋势。有效温度小于3 000 K的光谱处理结果同样不稳定,大量光谱总的结果变得不可信。除此之外,信噪比高于600,波长小于370 nm或大于900 nm,有效温度高于50 000 K的光谱未做实验。另外,本方法适用于中低色散的光谱,高色散光谱需要对参数进行微调。

需要指出,对于本方法涉及的参数w, t, r0, s0, c,它们的选择是为了适应仪器特性,对于其他来源的光谱,该方法仍然可用但参数初值需要略微调整。

3 结论与展望

归一化是光谱数据处理的重要环节,这一步骤处理结果的准确度直接影响后续数据分析结果的正确性以及精确程度。为了获得更加准确的连续谱,在国内外研究的基础上提出了基于固定窗口划分的样条函数拟合方法,通过大量实验修正其中某些参数的经验取值。随后,利用更多的数据进行验证,结果表明,本方法在拓宽适用范围的前提下实现了自动化处理,与其他方法相比具有更好的普适性,在需要大规模自动化处理的场合有独特的优越性。

恒星光谱的归一化方法有很多种,在不同的场合各有优劣。本文的方法致力于实现郭守敬望远镜海量光谱实时处理的自动化及普适性拓展。涉及到的参数值是在大量实验的基础上,由经验总结出来的。另外,在开发算法的过程中,有一个重要问题未得到很好的解决,即缺乏一个严格的最优解评价标准来衡量归一化结果的好坏,使得算法仅能对一些明显的错误(如连续谱出现负值)做出基本的容错处理。如果能够研究出可量化并且方便计算机实现的归一化结果评价标准,则可以对处理结果进行评价,通过参数自动修改后再次处理直至获得最优解。

随着天文数据的急剧增长,还需要更加高效和智能的归一化方法及数据处理程序。下一步将进行低温星、特殊星(如碳星等)的光谱归一化研究。

参考文献
[1] 赵永恒. LAMOST天体光谱巡天[J]. 物理, 2015, 44(4): 205–212 DOI: 10.7693/wl20150401
[2] ZHAO J K, ZHAO G, CHEN Y Q, et al. Automatic normalization and equivalent-width measurement of high-resolution stellar spectra[J]. Chinese Journal of Astronomy and Astrophysics, 2006, 6(6): 689–696. DOI: 10.1088/1009-9271/6/6/07
[3] 赵颖玥, 罗阿理. DA白矮星视向速度测量[J]. 天文研究与技术, 2019, 16(1): 1–7 DOI: 10.3969/j.issn.1672-7673.2019.01.001
[4] STARCK J L, SIEBENMORGEN R, GREDEL R. Spectral analysis using the wavelet transform[J]. The Astrophysical Journal, 1997, 482(2): 1011–1020. DOI: 10.1086/apj.1997.482.issue-2
[5] 赵景昆, 赵刚, 陈玉琴. 高分辨率恒星光谱处理与元素丰度分析软件平台[J]. 天文研究与技术——国家天文台台刊, 2007, 4(2): 153–158
[6] 赵瑞珍, 罗阿理. 天体光谱信号的连续谱归一化新方法[J]. 光谱学与光谱分析, 2006, 26(3): 587–590 DOI: 10.3321/j.issn:1000-0593.2006.03.052
[7] LEE Y S, BEERS T C, SIVARANI T, et al. The SEGUE stellar parameter pipeline. Ⅰ. description and comparison of individual methods[J]. The Astronomical Journal, 2008, 136(5): 2022–2049. DOI: 10.1088/0004-6256/136/5/2022
[8] 曾谨言. 量子力学:卷Ⅰ[M]. 北京: 科学出版社, 2013: 3-4.
[9] 于敬敬, 潘景昌, 孟凡龙, 等. 基于距离度量的LAMOST光谱中连续谱异常的自动检测[J]. 光谱学与光谱分析, 2017, 37(7): 2246–2249
[10] SCHUMAKER L L. Spline functions: basic theory[M].[S.l.]: Cambridge University Press, 2007.
[11] 吴明磊, 潘景昌, 衣振萍, 等. 恒星低质量光谱的连续谱拟合方法[J]. 光谱学与光谱分析, 2018, 38(3): 963–967
由中国科学院国家天文台主办。
0

文章信息

罗锋, 刘超, 赵永恒
Luo Feng, Liu Chao, Zhao Yongheng
基于样条函数的恒星光谱自动归一化方法
Automatic Normalization Method of Stellar Spectrum Based on Spline Function
天文研究与技术, 2019, 16(3): 300-311.
Astronomical Research and Technology, 2019, 16(3): 300-311.
收稿日期: 2018-12-18
修订日期: 2019-03-01

工作空间