流星頻譜提取方法的研究與實(shí)現(xiàn)

朱正平，曹澍宇，李國主

（1 中南民族大學(xué) 電子信息工程學(xué)院，武漢 430074；2 中國科學(xué)院地質(zhì)與地球物理研究所 中國科學(xué)院地球與行星物理重點(diǎn)實(shí)驗(yàn)室，北京 100029）

流星體是星際空間中的塵埃、顆粒和天體殘骸，受地球引力影響進(jìn)入地球大氣層，當(dāng)流星體進(jìn)入高度為130～60 km 的大氣層時，流星體與稠密大氣分子相互作用，溫度升高導(dǎo)致其表層物質(zhì)熔融、蒸發(fā)，流星體蒸汽中的原子和離子被電離發(fā)光（流星現(xiàn)象）［1］，離化的原子和離子會在電離層E 層中長期存在，從而形成一個持久層，即Es 層［2-3］，另外，流星體在大氣中的行進(jìn)路徑上也會產(chǎn)生等離子體尾跡［4-5］，這些等離子體尾跡也會顯著地?cái)_亂電離層.可以利用光學(xué)設(shè)備獲得流星體的頻譜信息，通過分析設(shè)備記錄的流星頻譜推斷流星體的化學(xué)組成以及流星體表層物質(zhì)蒸發(fā)的溫度，這對研究流星體對空間環(huán)境的影響以及流星體的起源有重要意義［6-7］.

目前，國外多地的流星光學(xué)觀測系統(tǒng)使用成像光譜儀直接獲取流星頻譜信息，盡管該方式可以獲取高精度的流星頻譜，但是由于流星出現(xiàn)時間和空間的不確定性導(dǎo)致光譜儀只能采集到可預(yù)測的流星雨事件的頻譜信息，并且獲取到的頻譜數(shù)量較少［8-10］.為了獲取大量的流星頻譜數(shù)據(jù)，可以使用帶有低照度CCD（Charge-coupled Device，電荷耦合元件）圖像傳感器的攝像機(jī)和光柵組合不間斷的監(jiān)控夜晚天空，利用光柵的衍射效應(yīng)獲取流星頻譜，通過攝像機(jī)將流星頻譜記錄為視頻并保存.但是利用光柵-攝像機(jī)組合獲取的流星頻譜視頻會出現(xiàn)大量噪聲以及頻譜圖像彎曲和非線性色散等問題.為了從光柵記錄的頻譜視頻中提取出流星頻譜信息，本文提出了一種從光柵獲得的頻譜視頻中提取流星頻譜的方法，通過對流星頻譜圖像進(jìn)行去背景場、圖像變換、傾斜校正、波長定標(biāo)等處理，消除圖像噪聲與頻譜圖像畸變問題，使得頻譜在圖像上線性分布并且具有一個恒定的色散常數(shù)，最后為了提高方法的應(yīng)用性，將方法的各個處理步驟集成到一款帶有圖形交互界面的程序中，該程序可以應(yīng)用在流星觀測系統(tǒng)的數(shù)據(jù)處理部分，能夠?qū)崟r處理海量的流星頻譜視頻數(shù)據(jù)，最后給出了使用該程序從流星頻譜視頻中提取出的頻譜強(qiáng)度曲線并做了簡要分析.

1 流星頻譜觀測方法

為了獲取大量的流星頻譜數(shù)據(jù)，使用圖1 所示的流星觀測系統(tǒng)，通過在攝像機(jī)鏡頭前安裝光柵來獲取流星頻譜.流星體發(fā)出的光經(jīng)過密集的光柵刻線被分離為不同波長的光并形成條狀的頻譜，使用CCD 攝像機(jī)將流星頻譜記錄為視頻，最后經(jīng)過UFOCapture 軟件對采集到的頻譜視頻進(jìn)行檢測［11］，保留包含流星頻譜信息的視頻片段.

圖1 流星觀測系統(tǒng)Fig.1 Meteor observation system

光柵-攝像機(jī)組合獲取的流星頻譜視頻存在的最大問題就是頻譜圖像的畸變以及頻譜在圖像上的非線性分布，這兩個問題主要是流星體發(fā)出的光經(jīng)光柵衍射和鏡頭折射產(chǎn)生的，因此流星觀測系統(tǒng)的光路結(jié)構(gòu)對于校正頻譜圖像畸變十分重要.

圖2 是流星觀測系統(tǒng)所用的光柵衍射原理圖，流星發(fā)出的光以向量f(A，B，C)表示，經(jīng)光柵衍射后 的 出 射 光 以 向 量 f′(A′，B′，C′)表 示 ，光 柵 方 程如下［12］：

圖2 光柵衍射系統(tǒng)Fig.2 The theory of grating diffraction

λ表示光波長，m為光柵級數(shù)，級數(shù)為0時，頻譜不會發(fā)生衍射，G為光柵常數(shù).

圖3 是成像原理圖，光軸與Z 軸重合，光柵垂直于光軸安裝.流星體發(fā)出的光(A，B，C)經(jīng)過光柵衍射，不同波長的光被分離并且以不同的角度ρ出射，以(A′，B′，C′)代替其中波長為λ的衍射光，經(jīng)過鏡頭折射最終落在圖像傳感器上點(diǎn) P(x，y)，r 是 P 到光軸的距離，r可以使用衍射角ρ表示，如式（4），

圖3 成像原理Fig.3 Imaging theory

f 為鏡頭焦距，g(ρ)是與鏡頭投影屬性相關(guān)的函數(shù)［13］.

φ 是 P 在圖像平面上的方位角，P(x，y)可以使用r和方位角φ表示為：

而方位角φ和衍射角ρ都可以使用光向量表示為：

將等式（8），（9）代入式（5），（6），（1），（2）使用向量分量和衍射角表示：

式（10），（11）定義了流星體發(fā)出的光經(jīng)光柵衍射，不同波長的光被分離并映射到圖像平面P 點(diǎn)的過程，P 點(diǎn)的位置只與衍射角和衍射光的波長以及鏡頭屬性相關(guān)，適用于任何鏡頭，可以用于后續(xù)推導(dǎo)解決流星頻譜圖像畸變的方法.

2 流星頻譜提取方法與實(shí)現(xiàn)

流星頻譜圖像彎曲和流星頻譜的非線性分布是對流星頻譜提取影響最大的問題.為了解決該問題，根據(jù)式（10），（11），通過變換頻譜圖像的投影方式g(ρ)來消除頻譜圖像彎曲，使頻譜在圖像上線性分布，獲取頻譜的線色散值.在進(jìn)行圖像變換之前，需要對頻譜圖像進(jìn)行去背景場處理，消除背景噪聲對后期圖像變換以及頻譜提取的影響.在圖像變換后，頻譜在整個圖像平面上線性分布并且色散系數(shù)是一個恒定值，將頻譜圖像對應(yīng)的像素灰度值相加，獲得頻譜強(qiáng)度信息，使用光譜燈實(shí)驗(yàn)確定的恒定的線色散值disp（單位為nm/pixel，即單位像素上波長的變化率）對頻譜進(jìn)行波長分配，獲得頻譜強(qiáng)度曲線.

2.1 去背景場

去背景場的目的是消除觀測系統(tǒng)在記錄流星事件時，由月光、人造光源或背景恒星帶來的噪聲.由于這些噪聲光源在短時間內(nèi)可以看作是恒定不變的，因此可以使用不包含流星頻譜信息的圖像作為背景噪聲圖像.在使用UFOCapture 截取的流星頻譜觀測視頻中，前25 幀是不包含流星頻譜的圖像，將前25幀圖像的灰度值相加并取平均，得到背景噪聲數(shù)據(jù)，使用后續(xù)的流星頻譜圖像數(shù)據(jù)減去背景噪聲數(shù)據(jù)就可消除背景光源帶來的影響.

2.2 圖像變換

流星光源發(fā)出的光被光柵衍射分離為不同波長的光，這些不同波長的光以不同的衍射角ρ出射，并經(jīng)過鏡頭折射最終映射到CCD 上，理想鏡頭的g(ρ) = tanρ，代入式（10），（11），得到：

從上式看出，由于頻譜中不同譜線的衍射角不同，導(dǎo)致頻譜在圖像平面上的分布不是線性的，并且實(shí)際使用的鏡頭都有畸變，不是完美的正切投影，這就使得頻譜圖像彎曲問題變得更加復(fù)雜.對于這些問題，我們通過改變投影函數(shù)g(ρ)來消除頻譜的畸變.

假設(shè)投影方式為正弦投影，即g(ρ) = sinρ，由式（10），（11）得到：

通過將投影方式變換為正弦投影可以消除P點(diǎn)對衍射角ρ 的依賴，P 點(diǎn)在圖像平面上橫坐標(biāo)x 與波長λ為線性關(guān)系，使得頻譜在圖像平面上為線性分布.

圖像投影方式變換步驟如下：

如圖4，經(jīng)過正弦投影的像點(diǎn)到光軸中心的距離為 r′，

圖4 正切投影與正弦投影Fig.4 Tangential projection and orthographic projection

利用式（4）的反函數(shù)求得ρ代入上式（17）得到：

求出上式（18）的反函數(shù)：

對于每個不同的鏡頭，由于自身的畸變導(dǎo)致各自的投影函數(shù)不同，為了完成圖像變換需要得到每個光學(xué)系統(tǒng)的投影函數(shù)g(ρ)，因?yàn)間(ρ)和sin(ρ)都是單調(diào)遞增的奇函數(shù)，因此可以做麥克勞林展開：

式（21）是理想鏡頭的投影變換公式，對于式（20）中實(shí)際鏡頭的圖像變換參數(shù)的確定，可以通過在式（21）中引入鏡頭畸變系數(shù)p 來代替實(shí)際鏡頭的畸變，如式（22），利用光譜燈實(shí)驗(yàn)數(shù)據(jù)，使用式（22）擬合獲取畸變參數(shù)和圖像變換系數(shù).

2.3 傾斜校正

經(jīng)過頻譜圖像變換，消除了頻譜圖像彎曲的現(xiàn)象.為了提取出一維頻譜，需要將頻譜圖像中頻譜對應(yīng)的那幾行像素的灰度值相加，但是在相機(jī)拍攝曝光期間，流星體是快速移動的，因此記錄下的頻譜會在水平和垂直方向出現(xiàn)傾斜.為了保證一維頻譜提取的正確性，需要對傾斜頻譜進(jìn)行校正，使頻譜平行與圖像平面X，Y 軸（本文中圖像坐標(biāo)系的X軸為圖像水平方向，Y 軸為圖像縱向，坐標(biāo)原點(diǎn)為圖像左上角）.我們可以旋轉(zhuǎn)頻譜圖像以校正頻譜傾斜，對X軸方向傾斜的校正：

rows為沿Y軸移動的行數(shù)，a為設(shè)置的旋轉(zhuǎn)系數(shù)，x為選中區(qū)域的各個像素點(diǎn)的橫坐標(biāo)，x0是選中區(qū)域中心點(diǎn)的坐標(biāo).

相對的，對Y軸方向的傾斜校正為：

cols 為沿 X 軸移動的列數(shù)，b 為設(shè)置的旋轉(zhuǎn)系數(shù)，y 為選中區(qū)域的各個像素點(diǎn)的縱坐標(biāo)，y0是選中區(qū)域中心點(diǎn)的縱坐標(biāo).

2.4 波長定標(biāo)

經(jīng)過上述多個步驟提取出的流星頻譜信息是每個像素點(diǎn)上的信號強(qiáng)度，波長定標(biāo)的目的就是將像素點(diǎn)位置轉(zhuǎn)換為對應(yīng)的波長.由于在圖像變換之后，頻譜在圖像平面上線性分布，即波長與對應(yīng)像素點(diǎn)的位置的關(guān)系是線性的.因此，對于頻譜的波長定標(biāo)，可以通過選取零階譜線或已知波長的譜線，如Na，Mg（零階譜線沒有出現(xiàn)在圖像內(nèi)）作為參考譜線，利用式（25）對各個譜線進(jìn)行波長分配.

disp 為線色散值，x0，λ0分別為參考譜線的所在像素位置以及參考譜線的波長，x，λ 是目標(biāo)譜線的像素位置和波長.

2.5 軟件設(shè)計(jì)

在使用該方法處理流星觀測系統(tǒng)記錄的實(shí)例流星頻譜視頻數(shù)據(jù)時，由于處理步驟較多，使得該方法的應(yīng)用性不強(qiáng).針對該問題，該方法對流星頻譜視頻處理的各個步驟被集成到一款帶有圖像交互界面的程序中，圖5 顯示了該程序的執(zhí)行流程以及幾個關(guān)鍵步驟的實(shí)現(xiàn)方法.

圖5 流星頻譜視頻處理程序流程圖Fig.5 Flow chart of Meteor spectrum video processing program

3 實(shí)驗(yàn)結(jié)果與分析

為了確定式（20）圖像變換公式以及頻譜圖像處理過程中所需要的各項(xiàng)參數(shù)，在觀測設(shè)備安裝后，需要對每個攝像機(jī)-光柵組合進(jìn)行光譜燈實(shí)驗(yàn).

本次實(shí)驗(yàn)設(shè)備：光譜燈為OceanOptics HG-1，光柵參數(shù)為600 L/mm（50 mm×50 mm），CCD 參數(shù)為最低照度0.0001 Lux，分辨率為576×720 pixel，鏡頭為Hikvision（f/0.95）.

3.1 參數(shù)確定

根據(jù)式（22），令r = x - x0，r′= x′- x0，可得

其中x 為圖像中頻譜實(shí)際像素位置，為擬合的頻譜像素位置，x0為圖像實(shí)際中心像素位置，p 為引入的鏡頭畸變系數(shù).2

其中 λ 是 x 對應(yīng)像素位置的頻譜波長，λ0是 x0處對應(yīng)的頻譜波長，disp0為每次實(shí)驗(yàn)的頻譜線色散值.

我們使用Hg-Ar 光譜燈進(jìn)行實(shí)驗(yàn)，圖6 是多次實(shí)驗(yàn)拍攝的光譜燈頻譜組合圖像，從中可以獲取多組[λ，x，y]，即每次實(shí)驗(yàn)拍攝的頻譜中各譜線的波長以及對應(yīng)的橫坐標(biāo)和縱坐標(biāo)，根據(jù)式（27），對獲取的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行最小二乘擬合，得到多組譜線所在的縱坐標(biāo)的平均值.

通過擬合的多組數(shù)據(jù)取平均得到：x0=360.83 pixel，c2= 4.6812e-7.

將縱向像素位置-y0看成disp0的函數(shù)，求出圖7所示擬合曲線，如圖7所示，其中藍(lán)色數(shù)據(jù)點(diǎn)為上述擬合得到的多組(-y0，disp0)，紅色數(shù)據(jù)點(diǎn)為擬合曲線極值點(diǎn)，對應(yīng)的就是頻譜的分辨率即線色散值disp =1.3470 nm/pixel以及圖像中心y0= 290.20 pixel.

圖7 縱向像素值與線色散值擬合曲線Fig.7 The fitting profiles of longitudinal pixel value and dispersion value

將上述得到的畸變參數(shù)c2代入式（22）可以得到：

3.2 頻譜提取參數(shù)驗(yàn)證

圖8（a）是光譜燈實(shí)驗(yàn)圖像，從圖中可以看出，在圖像頂部和底部頻譜出現(xiàn)了明顯的彎曲現(xiàn)象，使用上述求得的圖像中心坐標(biāo)為（360.83，290.20），利用式（29）對光譜燈實(shí)驗(yàn)圖像進(jìn)行投影變換，圖8（b）是變換后的圖像，從圖8（b）可以看出頂部，底部的頻譜彎曲明顯被消除.

圖8 頻譜圖像變換Fig.8 Spectral image transformation

圖9 是Sirius 標(biāo)準(zhǔn)頻譜與測量頻譜的對比圖，由于沒有校正CCD 每個像元對不同波長信號的響應(yīng)度，所以與標(biāo)準(zhǔn)頻譜對比，在一些頻段差異比較大，但在圖中標(biāo)注的幾條Sirius 的特征吸收譜線上，測量頻譜與標(biāo)準(zhǔn)頻譜的譜線重疊，由于該方法對頻譜的波長定標(biāo)是使用光譜燈實(shí)驗(yàn)確定的線色散值disp 結(jié)合式（25）給頻譜中的各譜線分配波長，通過幾條重疊的特征吸收譜線驗(yàn)證了譜線在圖像上的分布是線性的，初步驗(yàn)證了該方法的有效性.

圖9 標(biāo)準(zhǔn)頻譜與測量頻譜對比圖Fig.9 Comparison of standard spectrum and measured spectrum

為進(jìn)一步定量分析參數(shù)的準(zhǔn)確性以及該方法的提取效果，假設(shè)處理后的頻譜是線性分布的，并且線色散值disp為1.3470 nm/pixel，提取三次光譜燈實(shí)驗(yàn)圖像中的譜線進(jìn)行驗(yàn)證，如圖10 所示，橫坐標(biāo)為波長與頻譜階數(shù)的乘積，縱坐標(biāo)為處理后的提取譜線與已知的光譜燈實(shí)驗(yàn)譜線的誤差，從圖10中可以看出，最大提取誤差在2.8 nm 左右，提取誤差的均方根為0.8965 nm，結(jié)合線色散值1.3740 nm/pixel，像素誤差約為0.66 pixel，這個誤差是可以接受的，結(jié)合實(shí)驗(yàn)過程以及譜線提取方法，判斷誤差應(yīng)該是來自光源的離散性.

圖10 譜線提取誤差Fig.10 Spectral line extraction error

3.3 實(shí)例分析

本次使用的是 2019.08.06 21∶44∶48 UT 北京空間環(huán)境國家野外站海南樂東站(18.4°N，109°E)觀測到的英仙座流星雨事件，使用本文設(shè)計(jì)的程序進(jìn)行實(shí)例分析.本次觀測視頻幀率25 fps，光柵刻線600 L/mm.

圖11 分別是去背景場之前的圖像以及去背景場后的流星頻譜圖像，在圖11（a）中明顯可以看到人造光源產(chǎn)生的大面積的白色噪點(diǎn)以及圖中標(biāo)注的背景恒星，這些噪點(diǎn)會對后續(xù)的頻譜提取造成嚴(yán)重影響，經(jīng)過去背景場后，從圖11（b）中可以看出，由人造光源以及背景恒星產(chǎn)生的噪聲被消除.

圖11 去背景場Fig.11 Subtraction of background

圖12 （a）是傾斜校準(zhǔn)之前的頻譜圖像，從圖中可以看出頻譜中各譜線在圖像縱向明顯傾斜，從圖12（b）中可以看出，經(jīng)過對頻譜圖像的傾斜校正，明顯改善了譜線在圖像縱向的傾斜現(xiàn)象.

圖12 頻譜傾斜校正Fig.12 Tilt spectrum correction

通過上述幾個步驟，可以將圖像中頻譜對應(yīng)幾行像素的灰度值相加，提取出像素位置和信號強(qiáng)度的關(guān)系，在波長校準(zhǔn)中，選擇零階譜線或者已知波長的譜線作為參考譜線，使用頻譜燈校準(zhǔn)過程中得到的線色散值disp對頻譜進(jìn)行波長分配，繪制出頻譜圖.

圖13 是利用本文所開發(fā)的程序從實(shí)例數(shù)據(jù)中提取出的頻譜圖，頻譜線色散值disp為1.8204 nm/pixel.波長為0的譜線是零階譜線也就是通過零級光柵沒有發(fā)生衍射的光線，零階譜線左側(cè)是通過-1 級光柵的-1 階譜線，由于本觀測系統(tǒng)只能響應(yīng)390 nm～800 nm范圍的頻譜，因此本文主要分析該范圍內(nèi)的頻譜.根據(jù)Borovicka 在1994 年已驗(yàn)證的譜線信息［14］，從圖上可以看到如Mg-I-2（518.5 nm）、Na-I-1（589.5 nm）以及大氣譜線N-I-21（648.6 nm），N-I-3（746.8 nm），O-I-1（777.4 nm）.Mg-I-2 是條高溫譜線，在 它 左 側(cè) 的 是 Ca-I-2（422.6 nm）以 及 Mg-II-4（447.9 nm），由于沒有校準(zhǔn)設(shè)備的頻譜響應(yīng)，所以無法得到流星體中元素的相對豐度.但在400～600 nm范圍內(nèi)的譜線與Borovicka 和Betlem 在1997 年觀測到的英仙座流星雨的譜線類似［15］并且Ca-I-2（422.6 nm）、Mg-II-4（447.9 nm）、Mg-I-2（518.5 nm）、Na-I-1（589.5 nm）譜線得到了驗(yàn)證.

圖13 流星頻譜曲線Fig.13 Meteor spectrum profiles

4 結(jié)語

本文研究了流星頻譜的提取方法并開發(fā)了一款帶有交互式的圖形界面的程序用于處理流星頻譜視頻.通過去背景場消除噪聲，通過圖像變換消除了頻譜圖像彎曲，使頻譜在圖像上線性分布，通過傾斜校正消除了頻譜傾斜，提取出了一維頻譜，最后利用光譜燈實(shí)驗(yàn)確定的線色散值對頻譜進(jìn)行波長分配，得到了頻譜強(qiáng)度曲線.通過將所提取的測量頻譜和標(biāo)準(zhǔn)頻譜對比，驗(yàn)證了方法的有效性和準(zhǔn)確性.使用該方法以及在此基礎(chǔ)上開發(fā)的程序可以從觀測數(shù)據(jù)中獲取大量的頻譜信息，這對于進(jìn)一步分析流星體的組成，獲取流星體的起源、衰變等信息，研究流星體對地球空間環(huán)境的影響有重要意義.作為下一步的工作，我們將利用得到的頻譜強(qiáng)度曲線結(jié)合設(shè)備的頻譜響應(yīng)曲線，獲取流星體化學(xué)成分的相對豐度，分析流星體對空間環(huán)境的影響以及太陽系的起源和活動信息.

致謝：感謝北京空間環(huán)境國家野外觀測站提供的觀測數(shù)據(jù).