BFOSC長縫光譜自動處理軟件

宋得陽，張君波，王 靚*

(1. 中國科學(xué)院南京天文光學(xué)技術(shù)研究所，江蘇 南京 210042；2. 中國科學(xué)院天文光學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室(南京天文光學(xué)技術(shù)研究所)，江蘇 南京 210042；3. 中國科學(xué)院大學(xué)，北京 100049；4. 中國科學(xué)院國家天文臺，北京 100101；5. 中國科學(xué)院光學(xué)天文重點(diǎn)實(shí)驗(yàn)室 (國家天文臺)，北京 100101)

我國自主研制的第1臺2 m級望遠(yuǎn)鏡——國家天文臺興隆基地的2.16 m望遠(yuǎn)鏡，配備了北京暗天體攝譜儀[1]。這臺儀器具備5種不同的觀測模式，其中長縫光譜是應(yīng)用最為廣泛的一種模式，可以在不同的縫寬、光柵和濾光片組合之間進(jìn)行切換。2016年更換CCD相機(jī)后最新的色散、波長范圍等參數(shù)如表1[1-3]。據(jù)公開資料統(tǒng)計(jì)，2021年興隆觀測基地為2.16 m望遠(yuǎn)鏡用戶分配的252個(gè)觀測夜中，使用BFOSC的有165夜，占比超過65%。

表1 BFOSC不同光柵及基本性能參數(shù)Table 1 Grating parameters of BFOSC

BFOSC測光模式下的數(shù)據(jù)可以使用文[4]開發(fā)的程序處理。但長縫光譜部分尚沒有專門的數(shù)據(jù)處理軟件，觀測者一般自行使用IRAF[5-6]等軟件利用交互方式對原始光譜圖像進(jìn)行處理，獲得一維譜。但這類軟件并不是針對BFOSC特點(diǎn)設(shè)計(jì)的[7]，沒有針對長縫光譜的特點(diǎn)進(jìn)行優(yōu)化，且美國國立光學(xué)天文臺(National Optical Astronomy Observatory, NOAO)已經(jīng)宣布自2013年起不再對IRAF進(jìn)行官方支持。

本文針對BFOSC長縫光譜，開發(fā)了一套光譜處理軟件，實(shí)現(xiàn)了自動化處理光譜數(shù)據(jù)。該程序可以讀取用戶指定的光譜儀原始數(shù)據(jù)，自動進(jìn)行圖像合并、波長定標(biāo)、目標(biāo)源光譜抽取等步驟，在關(guān)鍵步驟中保存中間結(jié)果并進(jìn)行可視化。

程序引入了優(yōu)化抽譜方法，自動屏蔽了宇宙線，提高了最終數(shù)據(jù)產(chǎn)品的精度。用戶可以使用該程序快速實(shí)現(xiàn)批量數(shù)據(jù)處理，大大減輕了用戶的工作量，有效提升了工作效率。

1 系統(tǒng)整體框架

1.1 環(huán)境配置

本程序基于Python3編程語言編寫，依賴一些開源和免費(fèi)的第3方軟件包，包括Astropy(用于FITS文件的讀取與存儲)，Numpy，Scipy(用于科學(xué)計(jì)算)和Matplotlib(可視化)。本軟件在www.zenodo.org/record/7866030#.ZEi9uM5ByUk以開源形式發(fā)布。

1.2 程序結(jié)構(gòu)與流程

程序主要包括數(shù)據(jù)讀取、數(shù)據(jù)合并和光譜抽取3部分。其中數(shù)據(jù)讀取包括掃描生成觀測日志、讀取各個(gè)步驟數(shù)據(jù)。數(shù)據(jù)合并包括本底合并、平場合并與改正以及定標(biāo)燈紅藍(lán)端合并。光譜抽取包括本底改正、平場改正、目標(biāo)星定位、圖像畸變改正、背景改正、波長定標(biāo)和一維譜抽取等步驟。具體流程如圖1。用戶只需要輸入原始數(shù)據(jù)文件所在的路徑，程序自動執(zhí)行必要的步驟并存儲各步驟產(chǎn)生的中間結(jié)果，以便后續(xù)使用，對關(guān)鍵步驟進(jìn)行可視化并得到最終數(shù)據(jù)產(chǎn)品。

圖1 程序流程圖Fig.1 Pipeline flow chart

其中本底(bias)是指沒有光子輸入時(shí)，CCD像素在通電時(shí)(零秒曝)的讀出值。平場(flat)是指輸入光源為均勻光的情況下CCD獲得的圖像，反映了CCD的各個(gè)像素對光子的響應(yīng)差異。定標(biāo)燈是指使用空心陰極射線燈的光譜，包含波長已知的發(fā)射譜線，用于標(biāo)定CCD圖像上波長與像素之間的對應(yīng)關(guān)系?？茖W(xué)目標(biāo)是光譜儀的使用者拍攝的天體目標(biāo)，例如彗星、恒星和星系等。

1.3 程序使用

本文介紹的程序是針對BFOSC的長縫光譜設(shè)計(jì)，用于將光譜儀獲取的原始觀測數(shù)據(jù)(通常為FITS格式的光譜圖像)轉(zhuǎn)換為可供科學(xué)分析的一維光譜，即目標(biāo)天體的波長與流量的對應(yīng)關(guān)系，最終輸出數(shù)據(jù)以FITS格式保存。用戶需要在操作系統(tǒng)的終端運(yùn)行本程序，并輸入原始數(shù)據(jù)文件和觀測日志所在的路徑，程序自動執(zhí)行必要的步驟，存儲中間結(jié)果并進(jìn)行可視化，得到最終數(shù)據(jù)產(chǎn)品，即目標(biāo)天體的一維譜。這些結(jié)果和中間過程文件都存儲在指定的計(jì)算機(jī)目錄中。此外，用戶也可以更改運(yùn)行模式，手動對目標(biāo)星進(jìn)行標(biāo)記，程序會根據(jù)標(biāo)記位置尋找目標(biāo)星準(zhǔn)確輪廓。

2 實(shí)際處理問題與解決方法

2.1 掃描生成觀測日志

2.16 m望遠(yuǎn)鏡的觀測者在觀測同時(shí)經(jīng)常會記錄電子版的觀測日志，但是為了避免格式不一致，程序會掃描原始數(shù)據(jù)所在的計(jì)算機(jī)目錄，讀取其中FITS格式的文件，根據(jù)每個(gè)文件的頭部(Header)信息，結(jié)合用戶記錄的觀測日志，重新生成可以機(jī)讀的格式化觀測日志，包含文件名、圖像類型、觀測目標(biāo)源的名稱、曝光時(shí)間、觀測開始時(shí)間等信息。

2.2 本底改正

通常情況下，觀測者在觀測開始前或結(jié)束后在CCD相機(jī)快門關(guān)閉的情況下進(jìn)行多次零秒曝光，獲得一組本底圖像，通過圖像合并獲得本底的數(shù)值并抑制CCD讀出噪聲帶來的誤差。這個(gè)步驟提供了兩種不同的合并模式，即均值模式和中值模式，通常情況下默認(rèn)為均值模式。本底合并后生成本底圖像，并將合并后的本底圖像以FITS格式單獨(dú)存儲為一個(gè)文件。之后對需要操作的數(shù)據(jù)直接扣除保存的本底合并數(shù)據(jù)即完成本底改正。

2.3 目標(biāo)源定位

本文介紹的程序采用以下方法確定目標(biāo)天體的光譜在CCD圖像上的位置：從光譜圖像的左側(cè)開始，沿x軸方向每隔50個(gè)像素取圖像的截面最高處作為初始位置，取上下各20像素范圍內(nèi)的數(shù)據(jù)，用高斯函數(shù)進(jìn)行擬合，獲取縱向(y方向)的中心位置與半高全寬，并將中心位置當(dāng)作新的初始位置。重復(fù)上述過程，直到擬合參數(shù)收斂。接下來對所有y方向的中心位置與橫坐標(biāo)(x值)的關(guān)系進(jìn)行三階多項(xiàng)式函數(shù)擬合，作為目標(biāo)源光譜的中心位置曲線。光譜圖像如圖2。此外程序還為用戶提供了手動選取目標(biāo)源的模式，在該模式下用戶可以在彈出的圖形界面上用鼠標(biāo)標(biāo)注目標(biāo)源的大致位置，程序在其標(biāo)記的位置附近按上述流程尋找目標(biāo)天體的光譜。

圖2 (a)中心位置與半高全寬圖；(b)平均輪廓圖；(c)真實(shí)目標(biāo)源位置圖；(d)輪廓中心擬合圖Fig.2 (a)The central positions and full-width half-maximum (FWHM)of the spectral profiles，(b)the average profile is plotted in the lower left panel with black curve. (c)the target image，and (d)the positions and FWHMs of the center fittings

2.4 平場改正

由于CCD的各個(gè)像素對光子的響應(yīng)略有差異，即便輸入光源為均勻光，CCD獲得的圖像也是不均勻的。在拍攝過程中一般進(jìn)行多次曝光，將獲得的多幅平場圖像進(jìn)行疊加，以獲得較高的信噪比并抑制讀出噪聲。然后使用三次樣條插值補(bǔ)齊數(shù)據(jù)中的壞像元數(shù)據(jù)，之后使用Savitzky-Golay濾波過程[8]獲得主平場圖像。每一幅目標(biāo)天體的光譜圖像需要除以主平場圖像進(jìn)行平場改正。

2.5 視場光學(xué)畸變改正

由于CCD相機(jī)的光學(xué)畸變，CCD圖像上狹縫所成的像存在一定程度的彎曲。在波長定標(biāo)燈的光譜圖像上選取沿著y軸中間位置的一維譜作為基準(zhǔn)光譜，在y軸的上下兩個(gè)方向，每間隔200個(gè)像素進(jìn)行一次同樣的一維譜提取操作，將得到的光譜與基準(zhǔn)光譜計(jì)算交叉相關(guān)函數(shù)

(1)

其中，x和y分別為待改正的光譜與基準(zhǔn)光譜；n為光譜取點(diǎn)個(gè)數(shù)；k為相對像素差值；SDx和SDy分別為兩光譜的標(biāo)準(zhǔn)差。交叉相關(guān)函數(shù)最大值對應(yīng)的像素之差即為實(shí)際的像素偏移量。根據(jù)各行的像素差獲得改正函數(shù)，并生成改正圖像，如圖3，不同顏色的線代表y方向不同位置的光譜。經(jīng)過視場畸變改正后能有效提高一維光譜的波長定標(biāo)精度，避免系統(tǒng)誤差。

圖3 (a)畸變改正后的FeAr定標(biāo)燈譜的一維譜圖像；(b)畸變改正前的一維譜圖像Fig.3 (a)The FeAr spectra at different positions along the slit direction after the distortion correction；(b)before the distortion correction

2.6 背景改正

實(shí)際觀測過程中，我們獲得的光譜不僅來自于目標(biāo)源(圖4(a))，也有天光背景以及儀器內(nèi)部的雜散光。目標(biāo)星的光譜也包含這些成分，因此為了獲取目標(biāo)源的真實(shí)流量，我們需要對圖像數(shù)據(jù)進(jìn)行背景改正。

圖4 (a)目標(biāo)源天體的原始光譜圖像；(b)本例中測定天空背景采用的像素區(qū)域Fig.4 (a)The original 2D image of the target source；(b)the CCD image areas that are used to determine the background level

程序一般會控制目標(biāo)星使其落在沿y軸中間附近的位置。程序自動選取CCD圖像上距離目標(biāo)星足夠遠(yuǎn)、CCD像素一致性比較好的行(例如圖4(b)中第340～500行以及1 600～1 820行)作為測定背景光大小的區(qū)域。在選取過程中，程序會計(jì)算該區(qū)域每一行沿著x軸(即色散方向)的平均值和標(biāo)準(zhǔn)差，對于強(qiáng)度超過均值以上3倍標(biāo)準(zhǔn)差的行，可能包含目標(biāo)源附近其他天體的光譜，這種情況在密集星場中尤為常見。程序在計(jì)算背景值時(shí)對這些行進(jìn)行了屏蔽，以免獲得比真實(shí)值顯著偏高的背景流量。目標(biāo)源光譜附近的像素值減去背景光譜，即為目標(biāo)天體的光譜。

2.7 波長定標(biāo)

波長定標(biāo)的目的是在CCD圖像上建立像素位置與波長的映射關(guān)系。本文介紹的程序包含自動對波長進(jìn)行標(biāo)定，具體處理流程如下：(1)讀取該模式下提供的已經(jīng)定標(biāo)完成的模板光譜數(shù)據(jù)、模板特征發(fā)射線波長與像素位置對應(yīng)數(shù)據(jù)；其次對要處理的光譜數(shù)據(jù)與模板光譜數(shù)據(jù)進(jìn)行交叉相關(guān)函數(shù)計(jì)算，獲得兩個(gè)光譜數(shù)據(jù)位置的偏移量，將讀入的模板像素位置數(shù)據(jù)增加偏移量獲得新的像素位置數(shù)據(jù)。(2)依次取新的像素位置左右各10個(gè)像素的一維光譜，用廣義高斯函數(shù)Aexp[-(a-c)/α]β+B進(jìn)行擬合，獲得峰值位置對應(yīng)的像素?cái)?shù)值。結(jié)果如圖5和圖6，其中每幅圖中的圓點(diǎn)代表一維光譜數(shù)據(jù)，實(shí)線代表擬合曲線，虛線代表擬合得到的譜線中心位置。(3)用五階多項(xiàng)式函數(shù)進(jìn)行擬合，獲得波長與像素位置的對應(yīng)關(guān)系函數(shù)，并保存全部波長與像素位置數(shù)據(jù)組，畫出波長定標(biāo)函數(shù)圖像，如圖7。

圖5 用于波長證認(rèn)的FeAr譜線峰值位置擬合Fig.5 Peak positions of the identified FeAr lines

圖6 用于波長證認(rèn)的FeAr譜線峰值位置擬合Fig.6 Peak positions of the identified FeAr lines

圖7 (a)多項(xiàng)式擬合得到的波長(y)與像素(x)的關(guān)系，藍(lán)色點(diǎn)代表定標(biāo)譜線；(b)擬合殘差Fig.7 (a)The relation of fitted wavelength v.s. pixel numbers，where blue dots represent the fitting emission lines；(b)the residuals

此外，由于BFOSC采用的定標(biāo)燈本身的特性，其紅端存在大量明亮的Ar元素的發(fā)射線，而藍(lán)端以Fe線為主，它們的強(qiáng)度差異較大，如果曝光時(shí)間過長會導(dǎo)致紅端飽和，而曝光時(shí)間過短會導(dǎo)致藍(lán)端譜線沒有足夠的強(qiáng)度。因此，我們通常拍攝一次長曝光和一次短曝光，分別作為藍(lán)端和紅端的波長定標(biāo)譜。在本文的例子中，紅藍(lán)兩端的波長定標(biāo)圖像曝光時(shí)間分別為30 s和300 s。因此，在進(jìn)行定標(biāo)之前本方法對兩組光譜數(shù)據(jù)進(jìn)行截取之后合并處理生成一個(gè)新的定標(biāo)光譜圖像，處理過程中計(jì)算交叉相關(guān)函數(shù)時(shí)對紅端和藍(lán)端分別求解計(jì)算。

2.8 優(yōu)化抽譜

優(yōu)化抽譜方法采用加權(quán)求和，保證具有足夠流量的情況下包含盡量少的噪聲像素，對暗弱源尤其有效，能顯著提高信噪比[9]。優(yōu)化抽譜過程主要分為兩步：(1)獲得目標(biāo)星的精確空間輪廓，確定每個(gè)像素的加權(quán)值；(2)根據(jù)空間輪廓的形狀，排除明顯高于輪廓的像素點(diǎn)，即扣除宇宙線。在我們開發(fā)的光譜數(shù)據(jù)處理程序中，第1步首先默認(rèn)空間輪廓寬度為36個(gè)像素，輪廓方向每個(gè)像素具有相同的權(quán)重值，然后計(jì)算光譜數(shù)據(jù)方差，根據(jù)方差改正加權(quán)值與空間輪廓，重復(fù)此過程迭代2～3次，最終獲得精確的空間輪廓。第2步，先根據(jù)上一步得到的輪廓提取科學(xué)目標(biāo)的一維譜，并與普通抽譜方法得到的一維譜計(jì)算差值，對其中大于5倍標(biāo)準(zhǔn)偏差的像素點(diǎn)進(jìn)行屏蔽，我們認(rèn)為這些像素是受到宇宙線照射的像素。重復(fù)上述步驟直至收斂，至此我們獲得了屏蔽宇宙線后的科學(xué)目標(biāo)一維譜。

3 結(jié)果比較

為了驗(yàn)證程序的可靠性，我們使用IRAF對超新星SN2023eoc的同一組光譜數(shù)據(jù)進(jìn)行了處理并繪制了一維譜圖像，如圖8，其中，黃色線是本軟件自動化處理獲得的結(jié)果，藍(lán)色線是IRAF處理結(jié)果。由圖8可以看出，兩種處理方法結(jié)果基本一致，并且本軟件有效自動屏蔽了宇宙線，在流量比較低的藍(lán)端，信噪比明顯優(yōu)于IRAF的處理結(jié)果。IRAF需要在處理過程中對每條光譜手動屏蔽宇宙線，相比之下本文方法節(jié)省了處理時(shí)間且精度更高。

圖8 結(jié)果比較。藍(lán)色曲線為IRAF處理結(jié)果，黃色曲線為本軟件處理結(jié)果Fig.8 Comparison of results. The blue curve is the IRAF result，and the yellow curve is the result of this pipeline

4 總 結(jié)

本文系統(tǒng)介紹了一款針對2.16 m望遠(yuǎn)鏡BFOSC光譜儀的長縫光譜數(shù)據(jù)處理軟件，主要流程包括生成日志、本底和平場的合并與改正、目標(biāo)星定位、視場畸變改正、背景改正、波長定標(biāo)、優(yōu)化抽譜。本文的數(shù)據(jù)處理軟件使用方便，相比IRAF需要手動處理，本軟件只需要用戶輸入原始數(shù)據(jù)所在的路徑即可，處理過程透明可視，每個(gè)步驟完成都會生成相應(yīng)的圖片并保存數(shù)據(jù)以便重復(fù)使用，工作效率高，能批量化自動處理光譜數(shù)據(jù)且自動扣除宇宙線數(shù)據(jù)。與IRAF處理結(jié)果相比效果更好，具有較高的實(shí)用價(jià)值。

致謝：感謝興隆 2.16 m望遠(yuǎn)鏡全體工作人員的支持。