羅 鋒，劉 超，趙永恒

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

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

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

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

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

1 方法介紹

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

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

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

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

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

1.2 篩選參考點(diǎn)

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

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

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

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

(1)丟棄局部極小值點(diǎn)

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

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

(1)

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

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

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

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

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

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

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

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

(2)

(3)

(4)

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

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

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

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

(5)

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

(6)

同理，對(duì)于流量值[2]：

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

(7)

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

(8)

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

(9)

(10)

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

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

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

1.3 平滑度計(jì)算

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

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

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

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

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

s=c*s0.

(11)

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

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ù)參考點(diǎn)集中數(shù)據(jù)點(diǎn)數(shù)n選擇平滑度比例系數(shù)c

Table3Selectionofproportionalitycoefficientcofsmoothingfactoraccordingtothenumberofdatapointsninthereferencepointset

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

2 實(shí)驗(yàn)結(jié)果

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

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

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

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

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

2.2 誤差分析和適用范圍

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

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

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

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

圖10 不同波長(zhǎng)覆蓋范圍的歸一化結(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ù)譜，在國(guó)內(nèi)外研究的基礎(chǔ)上提出了基于固定窗口劃分的樣條函數(shù)擬合方法，通過(guò)大量實(shí)驗(yàn)修正其中某些參數(shù)的經(jīng)驗(yàn)取值。隨后，利用更多的數(shù)據(jù)進(jìn)行驗(yàn)證，結(jié)果表明，本方法在拓寬適用范圍的前提下實(shí)現(xiàn)了自動(dòng)化處理，與其他方法相比具有更好的普適性，在需要大規(guī)模自動(dòng)化處理的場(chǎng)合有獨(dú)特的優(yōu)越性。

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

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