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

羅 鋒，劉 超，趙永恒

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

作為國家重大科學工程，大天區面積多目標光纖光譜天文望遠鏡(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

根據表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(i-1)>f(i)并且f(i)

(1)

則將參考點i從Srp中刪去。每次從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 discardedttimes in the spectrum in Fig.1

圖4 圖2中光譜經t次丟棄局部極小值后的參考點集
Fig.4 The reference point set after the local minimum is discardedttimes 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。在幾何上，有：

(2)

(3)

(4)

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

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

對于Srp中的一個參考點i，設其與左邊數據點波長數值的差為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)

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

通過比值消去量綱，則不論w與f取何種單位，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中。

在構造樣條曲線時使用的平滑度為：

s=c*s0.

(11)

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

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 根據參考點集中數據點數n選擇平滑度比例系數c

Table3Selectionofproportionalitycoefficientcofsmoothingfactoraccordingtothenumberofdatapointsninthereferencepointset

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

2 實驗結果

2.1 不同情況下的處理結果

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

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

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

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

2.2 誤差分析和適用范圍

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

經實驗分析，本方法適用條件如下：信噪比范圍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，它們的選擇是為了適應儀器特性，對于其他來源的光譜，該方法仍然可用但參數初值需要略微調整。

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

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

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

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

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

3 結論與展望

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

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

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