流星频谱提取方法的研究与实现

朱正平，曹澍宇，李国主

（1 中南民族大学 电子信息工程学院，武汉 430074；2 中国科学院地质与地球物理研究所 中国科学院地球与行星物理重点实验室，北京 100029）

流星体是星际空间中的尘埃、颗粒和天体残骸，受地球引力影响进入地球大气层，当流星体进入高度为130～60 km 的大气层时，流星体与稠密大气分子相互作用，温度升高导致其表层物质熔融、蒸发，流星体蒸汽中的原子和离子被电离发光（流星现象）［1］，离化的原子和离子会在电离层E 层中长期存在，从而形成一个持久层，即Es 层［2-3］，另外，流星体在大气中的行进路径上也会产生等离子体尾迹［4-5］，这些等离子体尾迹也会显著地扰乱电离层.可以利用光学设备获得流星体的频谱信息，通过分析设备记录的流星频谱推断流星体的化学组成以及流星体表层物质蒸发的温度，这对研究流星体对空间环境的影响以及流星体的起源有重要意义［6-7］.

目前，国外多地的流星光学观测系统使用成像光谱仪直接获取流星频谱信息，尽管该方式可以获取高精度的流星频谱，但是由于流星出现时间和空间的不确定性导致光谱仪只能采集到可预测的流星雨事件的频谱信息，并且获取到的频谱数量较少［8-10］.为了获取大量的流星频谱数据，可以使用带有低照度CCD（Charge-coupled Device，电荷耦合元件）图像传感器的摄像机和光栅组合不间断的监控夜晚天空，利用光栅的衍射效应获取流星频谱，通过摄像机将流星频谱记录为视频并保存.但是利用光栅-摄像机组合获取的流星频谱视频会出现大量噪声以及频谱图像弯曲和非线性色散等问题.为了从光栅记录的频谱视频中提取出流星频谱信息，本文提出了一种从光栅获得的频谱视频中提取流星频谱的方法，通过对流星频谱图像进行去背景场、图像变换、倾斜校正、波长定标等处理，消除图像噪声与频谱图像畸变问题，使得频谱在图像上线性分布并且具有一个恒定的色散常数，最后为了提高方法的应用性，将方法的各个处理步骤集成到一款带有图形交互界面的程序中，该程序可以应用在流星观测系统的数据处理部分，能够实时处理海量的流星频谱视频数据，最后给出了使用该程序从流星频谱视频中提取出的频谱强度曲线并做了简要分析.

1 流星频谱观测方法

为了获取大量的流星频谱数据，使用图1 所示的流星观测系统，通过在摄像机镜头前安装光栅来获取流星频谱.流星体发出的光经过密集的光栅刻线被分离为不同波长的光并形成条状的频谱，使用CCD 摄像机将流星频谱记录为视频，最后经过UFOCapture 软件对采集到的频谱视频进行检测［11］，保留包含流星频谱信息的视频片段.

图1 流星观测系统Fig.1 Meteor observation system

光栅-摄像机组合获取的流星频谱视频存在的最大问题就是频谱图像的畸变以及频谱在图像上的非线性分布，这两个问题主要是流星体发出的光经光栅衍射和镜头折射产生的，因此流星观测系统的光路结构对于校正频谱图像畸变十分重要.

图2 是流星观测系统所用的光栅衍射原理图，流星发出的光以向量f(A，B，C)表示，经光栅衍射后 的 出 射 光 以 向 量 f′(A′，B′，C′)表 示 ，光 栅 方 程如下［12］：

图2 光栅衍射系统Fig.2 The theory of grating diffraction

λ表示光波长，m为光栅级数，级数为0时，频谱不会发生衍射，G为光栅常数.

图3 是成像原理图，光轴与Z 轴重合，光栅垂直于光轴安装.流星体发出的光(A，B，C)经过光栅衍射，不同波长的光被分离并且以不同的角度ρ出射，以(A′，B′，C′)代替其中波长为λ的衍射光，经过镜头折射最终落在图像传感器上点 P(x，y)，r 是 P 到光轴的距离，r可以使用衍射角ρ表示，如式（4），

图3 成像原理Fig.3 Imaging theory

f 为镜头焦距，g(ρ)是与镜头投影属性相关的函数［13］.

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

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

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

式（10），（11）定义了流星体发出的光经光栅衍射，不同波长的光被分离并映射到图像平面P 点的过程，P 点的位置只与衍射角和衍射光的波长以及镜头属性相关，适用于任何镜头，可以用于后续推导解决流星频谱图像畸变的方法.

2 流星频谱提取方法与实现

流星频谱图像弯曲和流星频谱的非线性分布是对流星频谱提取影响最大的问题.为了解决该问题，根据式（10），（11），通过变换频谱图像的投影方式g(ρ)来消除频谱图像弯曲，使频谱在图像上线性分布，获取频谱的线色散值.在进行图像变换之前，需要对频谱图像进行去背景场处理，消除背景噪声对后期图像变换以及频谱提取的影响.在图像变换后，频谱在整个图像平面上线性分布并且色散系数是一个恒定值，将频谱图像对应的像素灰度值相加，获得频谱强度信息，使用光谱灯实验确定的恒定的线色散值disp（单位为nm/pixel，即单位像素上波长的变化率）对频谱进行波长分配，获得频谱强度曲线.

2.1 去背景场

去背景场的目的是消除观测系统在记录流星事件时，由月光、人造光源或背景恒星带来的噪声.由于这些噪声光源在短时间内可以看作是恒定不变的，因此可以使用不包含流星频谱信息的图像作为背景噪声图像.在使用UFOCapture 截取的流星频谱观测视频中，前25 帧是不包含流星频谱的图像，将前25帧图像的灰度值相加并取平均，得到背景噪声数据，使用后续的流星频谱图像数据减去背景噪声数据就可消除背景光源带来的影响.

2.2 图像变换

流星光源发出的光被光栅衍射分离为不同波长的光，这些不同波长的光以不同的衍射角ρ出射，并经过镜头折射最终映射到CCD 上，理想镜头的g(ρ) = tanρ，代入式（10），（11），得到：

从上式看出，由于频谱中不同谱线的衍射角不同，导致频谱在图像平面上的分布不是线性的，并且实际使用的镜头都有畸变，不是完美的正切投影，这就使得频谱图像弯曲问题变得更加复杂.对于这些问题，我们通过改变投影函数g(ρ)来消除频谱的畸变.

假设投影方式为正弦投影，即g(ρ) = sinρ，由式（10），（11）得到：

通过将投影方式变换为正弦投影可以消除P点对衍射角ρ 的依赖，P 点在图像平面上横坐标x 与波长λ为线性关系，使得频谱在图像平面上为线性分布.

图像投影方式变换步骤如下：

如图4，经过正弦投影的像点到光轴中心的距离为 r′，

图4 正切投影与正弦投影Fig.4 Tangential projection and orthographic projection

利用式（4）的反函数求得ρ代入上式（17）得到：

求出上式（18）的反函数：

对于每个不同的镜头，由于自身的畸变导致各自的投影函数不同，为了完成图像变换需要得到每个光学系统的投影函数g(ρ)，因为g(ρ)和sin(ρ)都是单调递增的奇函数，因此可以做麦克劳林展开：

式（21）是理想镜头的投影变换公式，对于式（20）中实际镜头的图像变换参数的确定，可以通过在式（21）中引入镜头畸变系数p 来代替实际镜头的畸变，如式（22），利用光谱灯实验数据，使用式（22）拟合获取畸变参数和图像变换系数.

2.3 倾斜校正

经过频谱图像变换，消除了频谱图像弯曲的现象.为了提取出一维频谱，需要将频谱图像中频谱对应的那几行像素的灰度值相加，但是在相机拍摄曝光期间，流星体是快速移动的，因此记录下的频谱会在水平和垂直方向出现倾斜.为了保证一维频谱提取的正确性，需要对倾斜频谱进行校正，使频谱平行与图像平面X，Y 轴（本文中图像坐标系的X轴为图像水平方向，Y 轴为图像纵向，坐标原点为图像左上角）.我们可以旋转频谱图像以校正频谱倾斜，对X轴方向倾斜的校正：

rows为沿Y轴移动的行数，a为设置的旋转系数，x为选中区域的各个像素点的横坐标，x0是选中区域中心点的坐标.

相对的，对Y轴方向的倾斜校正为：

cols 为沿 X 轴移动的列数，b 为设置的旋转系数，y 为选中区域的各个像素点的纵坐标，y0是选中区域中心点的纵坐标.

2.4 波长定标

经过上述多个步骤提取出的流星频谱信息是每个像素点上的信号强度，波长定标的目的就是将像素点位置转换为对应的波长.由于在图像变换之后，频谱在图像平面上线性分布，即波长与对应像素点的位置的关系是线性的.因此，对于频谱的波长定标，可以通过选取零阶谱线或已知波长的谱线，如Na，Mg（零阶谱线没有出现在图像内）作为参考谱线，利用式（25）对各个谱线进行波长分配.

disp 为线色散值，x0，λ0分别为参考谱线的所在像素位置以及参考谱线的波长，x，λ 是目标谱线的像素位置和波长.

2.5 软件设计

在使用该方法处理流星观测系统记录的实例流星频谱视频数据时，由于处理步骤较多，使得该方法的应用性不强.针对该问题，该方法对流星频谱视频处理的各个步骤被集成到一款带有图像交互界面的程序中，图5 显示了该程序的执行流程以及几个关键步骤的实现方法.

图5 流星频谱视频处理程序流程图Fig.5 Flow chart of Meteor spectrum video processing program

3 实验结果与分析

为了确定式（20）图像变换公式以及频谱图像处理过程中所需要的各项参数，在观测设备安装后，需要对每个摄像机-光栅组合进行光谱灯实验.

本次实验设备：光谱灯为OceanOptics HG-1，光栅参数为600 L/mm（50 mm×50 mm），CCD 参数为最低照度0.0001 Lux，分辨率为576×720 pixel，镜头为Hikvision（f/0.95）.

3.1 参数确定

根据式（22），令r = x - x0，r′= x′- x0，可得

其中x 为图像中频谱实际像素位置，为拟合的频谱像素位置，x0为图像实际中心像素位置，p 为引入的镜头畸变系数.2

其中 λ 是 x 对应像素位置的频谱波长，λ0是 x0处对应的频谱波长，disp0为每次实验的频谱线色散值.

我们使用Hg-Ar 光谱灯进行实验，图6 是多次实验拍摄的光谱灯频谱组合图像，从中可以获取多组[λ，x，y]，即每次实验拍摄的频谱中各谱线的波长以及对应的横坐标和纵坐标，根据式（27），对获取的实验数据进行最小二乘拟合，得到多组谱线所在的纵坐标的平均值.

通过拟合的多组数据取平均得到：x0=360.83 pixel，c2= 4.6812e-7.

将纵向像素位置-y0看成disp0的函数，求出图7所示拟合曲线，如图7所示，其中蓝色数据点为上述拟合得到的多组(-y0，disp0)，红色数据点为拟合曲线极值点，对应的就是频谱的分辨率即线色散值disp =1.3470 nm/pixel以及图像中心y0= 290.20 pixel.

图7 纵向像素值与线色散值拟合曲线Fig.7 The fitting profiles of longitudinal pixel value and dispersion value

将上述得到的畸变参数c2代入式（22）可以得到：

3.2 频谱提取参数验证

图8（a）是光谱灯实验图像，从图中可以看出，在图像顶部和底部频谱出现了明显的弯曲现象，使用上述求得的图像中心坐标为（360.83，290.20），利用式（29）对光谱灯实验图像进行投影变换，图8（b）是变换后的图像，从图8（b）可以看出顶部，底部的频谱弯曲明显被消除.

图8 频谱图像变换Fig.8 Spectral image transformation

图9 是Sirius 标准频谱与测量频谱的对比图，由于没有校正CCD 每个像元对不同波长信号的响应度，所以与标准频谱对比，在一些频段差异比较大，但在图中标注的几条Sirius 的特征吸收谱线上，测量频谱与标准频谱的谱线重叠，由于该方法对频谱的波长定标是使用光谱灯实验确定的线色散值disp 结合式（25）给频谱中的各谱线分配波长，通过几条重叠的特征吸收谱线验证了谱线在图像上的分布是线性的，初步验证了该方法的有效性.

图9 标准频谱与测量频谱对比图Fig.9 Comparison of standard spectrum and measured spectrum

为进一步定量分析参数的准确性以及该方法的提取效果，假设处理后的频谱是线性分布的，并且线色散值disp为1.3470 nm/pixel，提取三次光谱灯实验图像中的谱线进行验证，如图10 所示，横坐标为波长与频谱阶数的乘积，纵坐标为处理后的提取谱线与已知的光谱灯实验谱线的误差，从图10中可以看出，最大提取误差在2.8 nm 左右，提取误差的均方根为0.8965 nm，结合线色散值1.3740 nm/pixel，像素误差约为0.66 pixel，这个误差是可以接受的，结合实验过程以及谱线提取方法，判断误差应该是来自光源的离散性.

图10 谱线提取误差Fig.10 Spectral line extraction error

3.3 实例分析

本次使用的是 2019.08.06 21∶44∶48 UT 北京空间环境国家野外站海南乐东站(18.4°N，109°E)观测到的英仙座流星雨事件，使用本文设计的程序进行实例分析.本次观测视频帧率25 fps，光栅刻线600 L/mm.

图11 分别是去背景场之前的图像以及去背景场后的流星频谱图像，在图11（a）中明显可以看到人造光源产生的大面积的白色噪点以及图中标注的背景恒星，这些噪点会对后续的频谱提取造成严重影响，经过去背景场后，从图11（b）中可以看出，由人造光源以及背景恒星产生的噪声被消除.

图11 去背景场Fig.11 Subtraction of background

图12 （a）是倾斜校准之前的频谱图像，从图中可以看出频谱中各谱线在图像纵向明显倾斜，从图12（b）中可以看出，经过对频谱图像的倾斜校正，明显改善了谱线在图像纵向的倾斜现象.

图12 频谱倾斜校正Fig.12 Tilt spectrum correction

通过上述几个步骤，可以将图像中频谱对应几行像素的灰度值相加，提取出像素位置和信号强度的关系，在波长校准中，选择零阶谱线或者已知波长的谱线作为参考谱线，使用频谱灯校准过程中得到的线色散值disp对频谱进行波长分配，绘制出频谱图.

图13 是利用本文所开发的程序从实例数据中提取出的频谱图，频谱线色散值disp为1.8204 nm/pixel.波长为0的谱线是零阶谱线也就是通过零级光栅没有发生衍射的光线，零阶谱线左侧是通过-1 级光栅的-1 阶谱线，由于本观测系统只能响应390 nm～800 nm范围的频谱，因此本文主要分析该范围内的频谱.根据Borovicka 在1994 年已验证的谱线信息［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 是条高温谱线，在 它 左 侧 的 是 Ca-I-2（422.6 nm）以 及 Mg-II-4（447.9 nm），由于没有校准设备的频谱响应，所以无法得到流星体中元素的相对丰度.但在400～600 nm范围内的谱线与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）谱线得到了验证.

图13 流星频谱曲线Fig.13 Meteor spectrum profiles

4 结语

本文研究了流星频谱的提取方法并开发了一款带有交互式的图形界面的程序用于处理流星频谱视频.通过去背景场消除噪声，通过图像变换消除了频谱图像弯曲，使频谱在图像上线性分布，通过倾斜校正消除了频谱倾斜，提取出了一维频谱，最后利用光谱灯实验确定的线色散值对频谱进行波长分配，得到了频谱强度曲线.通过将所提取的测量频谱和标准频谱对比，验证了方法的有效性和准确性.使用该方法以及在此基础上开发的程序可以从观测数据中获取大量的频谱信息，这对于进一步分析流星体的组成，获取流星体的起源、衰变等信息，研究流星体对地球空间环境的影响有重要意义.作为下一步的工作，我们将利用得到的频谱强度曲线结合设备的频谱响应曲线，获取流星体化学成分的相对丰度，分析流星体对空间环境的影响以及太阳系的起源和活动信息.

致谢：感谢北京空间环境国家野外观测站提供的观测数据.