基于樣條函數(shù)的恒星光譜自動歸一化方法*

羅 鋒，劉 超，趙永恒

(1. 中國科學(xué)院國家天文臺，北京 100101；2. 中國科學(xué)院大學(xué)，北京 100049)

作為國家重大科學(xué)工程，大天區(qū)面積多目標(biāo)光纖光譜天文望遠(yuǎn)鏡(the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST, 又叫郭守敬望遠(yuǎn)鏡)突破了天文望遠(yuǎn)鏡大視場與大口徑難以兼得的難題[1]，成為目前國際上口徑最大的大視場望遠(yuǎn)鏡，是我國光學(xué)望遠(yuǎn)鏡研制的又一里程碑。截至2018年6月，郭守敬望遠(yuǎn)鏡完成了7年的巡天觀測，獲取的光譜數(shù)首次超過千萬量級，成為世界上第一個獲取光譜數(shù)超千萬的巡天項目。遺憾的是，對于這107條光譜，目前沒有自動化方法應(yīng)用到歸一化上，本文主要對這方面進(jìn)行研究。

每條觀測譜都由連續(xù)譜、譜線和噪聲組成[2]。通過譜線可以對恒星的化學(xué)組成進(jìn)行分析，并且由于多普勒效應(yīng)，譜線還攜帶著視向速度的信息[3]。因此，從原始光譜中準(zhǔn)確地提取譜線是恒星光譜研究及光譜數(shù)據(jù)處理工作的重要步驟，對后續(xù)的研究有重要意義。

要消除連續(xù)譜首先要對連續(xù)譜進(jìn)行擬合，目前連續(xù)譜擬合方法主要有多項式擬合、中值濾波、小波濾波、形態(tài)濾波器等[4]。此外，利用天文軟件[5]進(jìn)行半自動化處理的方法是人工篩選提取認(rèn)為合適的數(shù)據(jù)點[6]并選擇函數(shù)形式，然后由軟件進(jìn)行多項式擬合得到連續(xù)譜。這種方法依賴人工參與，無法滿足目前巡天項目中海量光譜的實時處理需求。國外的斯隆數(shù)字巡天(Sloan Digital Sky Survey, SDSS)中的數(shù)據(jù)處理程序SSPP(Segue Stellar Parameter Pipeline)采用分段多項式擬合的方法[7]：首先將光譜分為紅端和藍(lán)端兩部分，去除強(qiáng)巴爾末線系后，分別對藍(lán)端使用9階、紅端使用4階多項式分別擬合，丟棄3個標(biāo)準(zhǔn)差以外的數(shù)據(jù)點；然后拼接兩段連續(xù)譜再進(jìn)行一次9階多項式擬合，得到最終的連續(xù)譜。這種方法稍顯繁瑣，且計算量大。

本文提出的方法是將光譜劃分為數(shù)個固定的(包含相同的像素數(shù)量)窗口，窗口內(nèi)選中值點，通過篩選策略丟棄部分參考點，這樣做能夠有效地避開發(fā)射線、吸收線及宇宙線的干擾。通過對參考點之間 “夾角” 的控制，使用非常小的平滑度構(gòu)造樣條曲線，實現(xiàn)對連續(xù)譜的精確擬合。

1 方法介紹

本文提出的連續(xù)譜歸一化方法包含3個步驟：(1)將一條待處理光譜根據(jù)數(shù)據(jù)點總數(shù)劃分為數(shù)個窗口，每個窗口選出流量中值點形成參考點集；(2)制定一些篩選策略，將參考點集中有可能影響擬合結(jié)果的數(shù)據(jù)點丟棄；(3)根據(jù)參考點流量最大值所在波長位置選取對應(yīng)的樣條函數(shù)平滑度經(jīng)驗值進(jìn)行擬合，得到連續(xù)譜。

1.1 劃分窗口構(gòu)造參考點集

為了有足夠的數(shù)據(jù)點作為擬合函數(shù)的控制點，根據(jù)輸入的待處理光譜數(shù)據(jù)點總數(shù)，采用如表1的策略劃分窗口。設(shè)數(shù)據(jù)點總數(shù)為n，每個窗口包含像素數(shù)，即窗寬為w。

表1 窗寬與數(shù)據(jù)點總數(shù)之間的對應(yīng)關(guān)系Table 1 The relationship between windows width and the total number of data points

根據(jù)表1確定窗寬后對光譜進(jìn)行窗口劃分，遍歷所有窗口，每個窗口內(nèi)選取流量值等于窗口中所有點中值的像素作為參考點，并將其加入?yún)⒖键c集Srp。

1.2 篩選參考點

1.1節(jié)中構(gòu)造的參考點集Srp只體現(xiàn)了光譜的大致形狀，其中有些點可能位于吸收線或者發(fā)射線上，如圖1、圖2。

圖1 有些參考點位于吸收線上
Fig.1 Some reference points are on the absorption line

圖2 有些參考點位于發(fā)射線上
Fig.2 Some reference points are on the emission line

理論上，連續(xù)譜應(yīng)非常接近黑體輻射譜，而黑體譜可以劃分為兩個單調(diào)區(qū)間，峰值左邊是單調(diào)增區(qū)間，右邊是單調(diào)減區(qū)間[8]，希望用來構(gòu)造擬合曲線的參考點符合并體現(xiàn)這一特征。為了實現(xiàn)這一目標(biāo)，采取兩個步驟對參考點集再次篩選。

(1)丟棄局部極小值點

設(shè)Srp中一個參考點的索引為i，流量為f(i)，若i滿足：

f(i-1)>f(i)并且f(i)

(1)

則將參考點i從Srp中刪去。每次從i=1開始到遍歷完Srp中所有參考點結(jié)束為一趟，在一趟丟棄結(jié)束時，有些點可能成為新的局部極小值點，此時需要再丟棄一趟。設(shè)丟棄的趟數(shù)為t，需要保證以下兩點：(1)對連續(xù)譜正常的光譜，t次丟棄之內(nèi)Srp中元素數(shù)量不再減少；(2)對連續(xù)譜異常[9]的光譜，t次丟棄后Srp中有足夠的參考點用來構(gòu)造擬合曲線。實驗表明，合適的t取值為t=len(Srp)//30 + 2。

經(jīng)這一步驟處理后，圖1和圖2中兩條光譜的參考點集如圖3、圖4。

圖3 圖1中光譜經(jīng)t次丟棄局部極小值后的參考點集
Fig.3 The reference point set after the local minimum is discardedttimes in the spectrum in Fig.1

圖4 圖2中光譜經(jīng)t次丟棄局部極小值后的參考點集
Fig.4 The reference point set after the local minimum is discardedttimes in the spectrum in Fig.2

(2)丟棄 “夾角” 較大的點

為了獲得更加精確的連續(xù)譜擬合曲線，需要較小的平滑度構(gòu)造樣條函數(shù)。如果有的參考點和其左右兩點連線構(gòu)成的夾角過大，則會給擬合曲線帶來較大的偏離，這樣的參考點應(yīng)從Srp中刪去，如圖5。

圖5 位于發(fā)射線上的數(shù)據(jù)點應(yīng)從參考點集中刪去
Fig.5 The data point on the emission line should be removed from the reference point set

如圖6，在參考點集中取一點i，設(shè)其對應(yīng)的波長為w(i)，流量為f(i)，與它左右兩點間連線的夾角為θ(i)，組成的三角形邊長分別為a,b,c。在幾何上，有：

(2)

(3)

(4)

圖6 使用余弦定理計算 “夾角”
Fig.6 Calculate the angle by use of Cosine theorem

根據(jù)余弦定理計算θ(i)的值。實際上，流量值并不對應(yīng)長度單位，w和f的單位很難統(tǒng)一。由于光譜數(shù)據(jù)中流量是未定標(biāo)的相對流量，直接用余弦定理計算的夾角隨著光譜流量值所在的區(qū)間不同及光譜數(shù)據(jù)處理中必要的乘除運(yùn)算而變化。為此做如下處理：

對于Srp中的一個參考點i，設(shè)其與左邊數(shù)據(jù)點波長數(shù)值的差為wd-，右邊為wd+，則：

wd-=w(i)-w(i-1)，

(5)

wd+=w(i+1)-w(i) .

(6)

同理，對于流量值[2]：

fd-=f(i)-f(i-1),

(7)

fd+=f(i+1)-f(i) .

(8)

引入比例因子rw,rf，令：

(9)

(10)

計算判據(jù)r=|rw-rf|的值并引入閾值r0。若r>r0，則認(rèn)為點i屬于夾角較大的點，會給精確擬合帶來干擾，應(yīng)從Srp中除去。通過實驗發(fā)現(xiàn)，效果較好的r0取值為0.5。

通過比值消去量綱，則不論w與f取何種單位，rw,rf都是不變的。使用rw,rf作為判斷依據(jù)尋找這樣的參考點更為可靠。圖5中光譜經(jīng)此步驟處理后的結(jié)果如圖7。

圖7 位于發(fā)射線上的參考點被成功移除
Fig.7 The reference point on the emission line was successfully removed

1.3 平滑度計算

樣條函數(shù)[10]擬合效果的好壞，除了與參考點的選取有關(guān)，還取決于平滑度的選擇。對于同一條光譜，同樣的參考點集，不同的平滑度會有不同的效果，如圖8。

圖8 平滑度對擬合曲線的影響
Fig.8 The influence of smoothing factor on fitting curve

較小的平滑度可以更好地照顧到光譜圖像上起伏的細(xì)節(jié)，但過小容易產(chǎn)生振動頻繁的現(xiàn)象，使擬合曲線無法體現(xiàn)連續(xù)譜；較大的平滑度不會有劇烈的振動，但是過大會使擬合曲線過于平緩，也無法很好地體現(xiàn)連續(xù)譜。應(yīng)該指出的是，對于同一個平滑度，在不同的流量區(qū)間，樣條曲線會有不同的表現(xiàn)。在擬合之前，將光譜數(shù)據(jù)進(jìn)行最大值標(biāo)準(zhǔn)化，即flux=flux/max(flux)。

經(jīng)實驗分析，溫度越高的光譜需要越小的平滑度照顧光譜在藍(lán)端較急劇的轉(zhuǎn)折，溫度越低的光譜需要越大的平滑度弱化參考點上下起伏的趨勢。而一條光譜的溫度在連續(xù)譜上體現(xiàn)為流量峰值所在的位置。采用在參考點集中選取流量最大值，用其所在波長位置估計該條光譜的溫度區(qū)間，選擇對應(yīng)的平滑度s0。為了能夠適用于局部歸一化的場合，對于參考點較少的情況，需要更小的平滑度。為此，引入默認(rèn)值為1的比例系數(shù)c，用于在參考點較少時將之與s0相乘的結(jié)果作為最終的平滑度取值。具體數(shù)據(jù)列在表2和表3中。

在構(gòu)造樣條曲線時使用的平滑度為：

s=c*s0.

(11)

表2 根據(jù)參考點集中流量最大值點所在波長位置選擇平滑度

Table 2 Selection of smoothing factor according to the wavelength positionof the maximum flux point in the reference point set

wave/nm370～476476～688688～794794～900s00.010.020.030.05

表3 根據(jù)參考點集中數(shù)據(jù)點數(shù)n選擇平滑度比例系數(shù)c

Table3Selectionofproportionalitycoefficientcofsmoothingfactoraccordingtothenumberofdatapointsninthereferencepointset

n1～2930～5960～8990～119120～149150～180c0.20.350.50.650.81

2 實驗結(jié)果

2.1 不同情況下的處理結(jié)果

為檢驗所提出方法的處理效果，從郭守敬望遠(yuǎn)鏡中隨機(jī)抽取了不同類型的光譜10 000條進(jìn)行實驗。從中選出不同溫度的光譜進(jìn)行擬合，如圖9。圖中左側(cè)藍(lán)色實線為原始光譜，紅色實線為擬合曲線；右側(cè)綠色實線為對應(yīng)的歸一化譜。由圖9可見，本方法在不同溫度時做到了較好的擬合。

如圖10，左側(cè)桃紅色實線是原始光譜，綠色實線為樣條函數(shù)擬合曲線；右側(cè)靛藍(lán)色實線為對應(yīng)的歸一化譜。選取不同的波長覆蓋范圍后，本方法能夠自動應(yīng)用于需要局部歸一化的場合，不需要人工干預(yù)。

如圖11，左側(cè)是原始光譜及擬合曲線，在不同的信噪比情況下，本方法擬合效果受到的影響不大。與文[11]中專為低質(zhì)量光譜開發(fā)的歸一化方法處理效果非常接近。

如圖12，在光譜中含有發(fā)射線及宇宙線的情況下，本方法仍能有效識別并進(jìn)行準(zhǔn)確的擬合。

2.2 誤差分析和適用范圍

由于儀器原因，光譜中常見的儀器形變模式有3種，如圖13(a)。為了測試本文提出的方法是否能夠有效地處理這些儀器導(dǎo)致的光譜形變，在有效溫度4 000～8 000 K的范圍內(nèi)隨機(jī)抽取9條光譜，每條光譜歸一化后的結(jié)果都乘以這3種形變曲線，并再次歸一化。圖13(b)的直方圖中顯示了兩次歸一化結(jié)果殘差的分布，可以看出本文算法的誤差平均值在10-3量級，誤差彌散在10-2量級。

經(jīng)實驗分析，本方法適用條件如下：信噪比范圍5～600，波長覆蓋范圍370～900 nm，有效溫度范圍3 000～50 000 K。其中，信噪比低于5的光譜處理結(jié)果不穩(wěn)定，總體精度有隨著信噪比下降而降低的趨勢。有效溫度小于3 000 K的光譜處理結(jié)果同樣不穩(wěn)定，大量光譜總的結(jié)果變得不可信。除此之外，信噪比高于600，波長小于370 nm或大于900 nm，有效溫度高于50 000 K的光譜未做實驗。另外，本方法適用于中低色散的光譜，高色散光譜需要對參數(shù)進(jìn)行微調(diào)。

需要指出，對于本方法涉及的參數(shù)w,t,r0,s0,c，它們的選擇是為了適應(yīng)儀器特性，對于其他來源的光譜，該方法仍然可用但參數(shù)初值需要略微調(diào)整。

圖9 不同溫度下的歸一化結(jié)果
Fig.9 Normalization results at different temperatures

圖10 不同波長覆蓋范圍的歸一化結(jié)果
Fig.10 Normalization results of different wavelength coverage range

圖11 不同信噪比的歸一化結(jié)果
Fig.11 Normalization results under different SNR

圖12 含有發(fā)射線或宇宙線情況下的歸一化結(jié)果
Fig.12 Circumstances containing emission lines or cosmic rays

圖13 郭守敬望遠(yuǎn)鏡的形變模式及歸一化誤差分析
Fig.13 Deformation mode in LAMOST and error analysis of normalization

3 結(jié)論與展望

歸一化是光譜數(shù)據(jù)處理的重要環(huán)節(jié)，這一步驟處理結(jié)果的準(zhǔn)確度直接影響后續(xù)數(shù)據(jù)分析結(jié)果的正確性以及精確程度。為了獲得更加準(zhǔn)確的連續(xù)譜，在國內(nèi)外研究的基礎(chǔ)上提出了基于固定窗口劃分的樣條函數(shù)擬合方法，通過大量實驗修正其中某些參數(shù)的經(jīng)驗取值。隨后，利用更多的數(shù)據(jù)進(jìn)行驗證，結(jié)果表明，本方法在拓寬適用范圍的前提下實現(xiàn)了自動化處理，與其他方法相比具有更好的普適性，在需要大規(guī)模自動化處理的場合有獨特的優(yōu)越性。

恒星光譜的歸一化方法有很多種，在不同的場合各有優(yōu)劣。本文的方法致力于實現(xiàn)郭守敬望遠(yuǎn)鏡海量光譜實時處理的自動化及普適性拓展。涉及到的參數(shù)值是在大量實驗的基礎(chǔ)上，由經(jīng)驗總結(jié)出來的。另外，在開發(fā)算法的過程中，有一個重要問題未得到很好的解決，即缺乏一個嚴(yán)格的最優(yōu)解評價標(biāo)準(zhǔn)來衡量歸一化結(jié)果的好壞，使得算法僅能對一些明顯的錯誤(如連續(xù)譜出現(xiàn)負(fù)值)做出基本的容錯處理。如果能夠研究出可量化并且方便計算機(jī)實現(xiàn)的歸一化結(jié)果評價標(biāo)準(zhǔn)，則可以對處理結(jié)果進(jìn)行評價，通過參數(shù)自動修改后再次處理直至獲得最優(yōu)解。

隨著天文數(shù)據(jù)的急劇增長，還需要更加高效和智能的歸一化方法及數(shù)據(jù)處理程序。下一步將進(jìn)行低溫星、特殊星(如碳星等)的光譜歸一化研究。