嫦娥三号月基光学望远镜图像中的宇宙线识别

刘 奇，王 竞，黄茂海，魏建彦

（1.中国科学院 国家天文台 空间天文与技术重点实验室，北京100101；2.中国科学院大学，北京100049；3.广西大学 物理科学与技术学院 广西相对论天体物理重点实验室，广西 南宁530004）

1 引 言

天文观测图像中的信息一般包含观测源、背景、各种噪声及宇宙射线等。其中，宇宙射线［1］是来自宇宙空间的各种高能粒子形成的射流，绝大多数是质子，还包含α粒子和少量其他原子核。在天文观测图像中宇宙线的灰度值通常高于周围像素点的灰度值，并且分布规律具有随机性。绝大多数的宇宙线在电荷耦合器件（Charge Cou⁃pled Device，CCD）图像中呈现“点”状或“线”状分布，面积大的宇宙线也会和目标源形态轮廓相似。由于绝大多数天文图像是地基观测，受到地球大气的影响，打在CCD靶面上的宇宙线是穿过大气层，与大气中的粒子发生相互作用后形成的次级粒子，常常被作为一种噪声信号来处理。而在空间天文图像中，由于不受地球大气层的影响，CCD图像中的宇宙线存在亮度高、同观测目标源轮廓相似等特点，当图像中出现非常多的宇宙线事件时对目标源的检测会造成很大影响。

无论CCD图像中的宇宙线呈现出哪种形态特征，后续处理分析的前提都是将它从图像中提取出来。对于天文CCD图像识别提取宇宙线，不同专家学者也提出了不同算法。基于多帧图像，运用图像对齐相减的方法确定宇宙线的位置。宇宙线打到CCD探测器上具有随机性的特点，通过对相同的视场进行连续多帧拍摄，然后依据序列图像信息确定宇宙线位置。Windhorst等［2］就是运用这种方法从哈勃空间望远镜宽视场相机图像中定位、剔除了宇宙线事件。此类算法到目前为止已经相当成熟，但其缺点是需要对同一视场进行多次拍摄，降低了仪器的观测时间和使用效率，并且无法通过单一图像来定位宇宙线。

基于单帧图像识别提取识别宇宙线的方法更具挑战性。从数字图像处理的角度出发，中值滤波可消除随机噪声，如果把宇宙线作为天文图像中的随机噪声，对单帧图像使用中值滤波可以很好地探测到宇宙线。但是由于天文CCD图像作为高精度数字图像，在中值滤波的同时会平滑掉恒星目标的边缘，改变了其有用信息。因此，传统的中值滤波并不适用于天文观测图像。Van等［3］提出了一种运用拉普拉斯边缘检测算子来检测宇宙线的方法。此方法通过构建和设置两个阈值来识别宇宙线。一个阈值量是反应每个像素点含有噪声的噪声比率，像素点噪声越大，这个噪声比率越高，大于一定阈值的噪声比率认为是宇宙线候选体；第二个比率是原始图像与拉普拉斯算子卷积后得到图像与图像的中频信息的商，图像的中频信息是依据星像有对称性而宇宙线没有的特点，通过与不同中值滤波模板进行卷积相减得到的。通过设置这两个阈值，算法不断迭代，直到不能检测到宇宙线为止。该方法在宇宙线检测上非常有效，但算法复杂，处理速度较慢，而且对于形态轮廓与星像相似的宇宙线检测准确率较低。相比于拉普拉斯边缘检测算法，Pych等［4］提出了一种简单、直接的直方图算法用来识别宇宙线。该算法直接将图像分为若干子图像，并分析其直方图。由于绝大多数子图像的直方图分布比较紧凑，一旦有宇宙线事件，直方图中在紧凑结构之外就会出现孤立的单点，从而被分辨出来。针对此方法，不同学者提出了不同的阈值设定定义，并且在算法实现过程中，其参数需要根据原始图像的不同做出调整。相比拉普拉斯边缘检测算法，算法运行效率虽高，但其探测能力明显低于前者。

通过统计分析月基光学望远镜（Lunar-based Ultraviolet Telescope，LUT）［5］数据，发现相比于地基望远镜，宇宙线打在LUT CCD靶面上的发生概率显著提高，除地基望远镜图像中出现的“点”状或“线”状特征外，还存在很多亮度高、影响面积大、形态轮廓与星像类似的宇宙线。对于此类形态轮廓与星像类似的宇宙线，传统算法无法精确检测处理。基于此问题，本文提出了一种运用天文位置定标来识别提取宇宙线的方法。此方法与上述文献中方法不同的是从天文图像处理角度出发，不会引入其他噪声影响，对于形态轮廓与星像相似的宇宙线识别探测效率更高。识别提取到宇宙线事件后，对其电子沉积数和入射角进行了分析，并将其形态进行提取、分类，制作出“CCD探测器宇宙线实体样本库”，为空间天文望远镜图像的仿真提供经验支持。

2 月基光学望远镜

LUT的结构与安装示意图如图1所示。望远镜的主镜直径为150 mm，焦比为F/3.75，采用了里奇-克列基昂（R-C）系统设计。LUT指向反光镜放在一个二维转动平台上，方位方向的指向范围是-28°～+13°，俯仰方向的指向范围是+20°～+38°，指向精度为0.05°。望远镜的视场（Field of View，FOV）是1.36°×1.36°，像元分辨率为4.76″。表1所示为LUT的主要技术指标，更详细的技术指标请参考文献［5］。

表1 月基光学望远镜的主要技术指标Tab.1 Main characteristics of LUT

图1 月基光学望远镜结构与安装示意图Fig.1 Schematic diagram of structure and assembly of Lunar-based Ultraviolet Telescope（LUT）

3 数据处理流程

通过分析经典的中值滤波法、拉普拉斯边缘检测法［3］、基于直方图的快速算法［4］以及万能噪声消除算法［6］等，结合LUT天文观测图像的特点，本文提出了基于天文位置定标的宇宙线识别提取算法，其数据处理流程如图2所示。主要步骤如下：（1）对原始天文图像进行预处理操作，包括对原始图像进行背景扣除和平场修正；（2）对图像进行天文位置定标；（3）对高于5倍背景噪声的信号进行提取，获得宇宙线候选体；（4）对提取出的信号与星表进行位置交叉证认，剔除与星表中天文位置匹配的恒星目标，剩下的即是要提取和识别的宇宙线事件。

图2 宇宙线提取数据处理流程Fig.2 Data processing outline of cosmic ray extraction

3.1 背景扣除

背景扣除的主要目的是为了剔除杂散光背景、CCD自身“热点”和“坏点”的影响，获得保留恒星目标和宇宙线信息的图像。LUT由于是在月昼期间观测，因此受到很强的杂散光影响。LUT的杂散光主要来自于太阳光的散射。杂散光的强度和大尺度结构随着LUT反光镜和太阳相对位置的变化而变化。孟宪民等［7］对LUT观测图像统计发现，在大多数情况下约0.5 h内拍摄的图像中杂散光的强度和大尺度结构基本不变。

本文根据LUT观测的时间对图片进行了分组，将连续观测1 900 s以内的序列图像分为一组。对某一帧图像，将同组中除自身以外的所有图像对齐后进行中值合并，宇宙线被剔除；同时，由于LUT步进指向跟踪模式下，恒星在连续曝光的图像上的位置总是有轻微的几个像元的变化，所以通过中值合并，所有的天体目标也都被剔除了。因此，合并所得图像仅保留了背景和CCD的“热点”和“坏点”，得到的背景作为本帧图像的背景模板。

最后，对每一帧图像，减去图像对应的背景模板，背景被扣除，同时CCD的“热点”和“坏点”也被剔除，得到的图像仅保留了恒星星像和宇宙线。

3.2 平场修正

平场修正是为了解决LUT CCD相机像面响应的非均匀性。LUT图像的平场分为小尺度平场和大尺度平场两部分［7］。小尺度平场反映了CCD像元间量子效率的不均匀性或像元有效面积的变化；大尺度不均匀性主要来自望远镜的光学系统，还有小部分来自于CCD镀膜导致的大尺度效率变化。

对于小尺度平场，通过拍摄LUT内部平场灯图像，即可测量像面上不同像元间响应的非均匀性。对于大尺度平场，是通过对标准星进行“抖动”观测完成的。最终平场是大尺度平场与LED平面场的乘积。

3.3 天文位置定标

本文天文位置定标所用的星表为第谷星表（Tycho_2.0）［8］，主要是因为第谷星表的天文坐标精确高并且与LUT的观测深度相近。但第谷星表采用双色测光，王竞等［9］利用恒星大气模型计算出了LUT相应的近紫外AB星等。

天文位置定标的关键是将CCD图像中的恒星星像位置和第谷星表中的天文坐标进行交叉认证。通过LUT的指向文件记录获得CCD图像中心位置的J2000天文坐标，精度通常优于1个角分。围绕此中心位置提取交叉所用的第谷星表。

对于每帧图像，这里使用SExtractor（Source-Extractor）软件［10］提取出5～10颗明亮的恒星目标，并对它进行快速测光，获得其近紫外AB星等。将恒星目标与对应的第谷星表进行位置匹配和亮度匹配。位置匹配交叉半径设置为10角秒（约为LUT图像的2个像元），亮度匹配阈值控制在2个星等以内（主要考虑了变星的影响）。

利用交叉匹配的结果生成CCD图像上X-Y位置与天文坐标的位置对应关系，即世界坐标系统（World Coordinate System，WCS）的位置关系。通过WCS实现了对CCD图像上每一个像元的对应天文位置。天文位置的定标精度优于1角秒。

3.4 提取宇宙线候选体

本文利用SEXtractor软件对LUT图像提取宇宙线候选体，输出的参数有X-Y位置坐标、J2000天文坐标、流量计数和椭率等。

提取候选体的判据有：（1）单像元计数大于5倍的背景噪声；（2）连通像元数大于等于4，这是由恒星星像决定的；（3）距离CCD四边大于20个像元。

3.5 认证识别宇宙线事件

宇宙线候选体样本中包含了一些恒星目标，需要剔除这些恒星才能准确识别出宇宙线事件。将宇宙线候选体样本与Tycho_2.0星表和US⁃NO_B1.0星表进行位置交叉认证，剔除样本中的恒星目标，获得宇宙线样本。

本文选取两个星表的原因主要考虑了以下两点：（1）Tycho_2.0星表的极限星等为11.5，无法涵盖LUT图像中的全部星像；（2）虽然US⁃NO_B1.0极限星等为21等，由于其亮星饱和问题，星表在亮星部分不完备。所以将两个星表综合利用，保证了交叉星表的完备性。

为选取合适的交叉半径，本文统计了6帧图像中的候选体与两个星表的最近距离。首先选取10角秒为半径与两个星表进行位置交叉，结果表明，相匹配的源中97.21%的距离都小于4角秒，小于LUT一个像元的尺寸4.76角秒。最终这里采用4角秒作为交叉半径。

图3为12帧图像中宇宙线候选体与两个星表交叉获得的星像数量统计图（彩图见期刊电子版）。图中，蓝色表示与Tycho_2.0星表交叉获得的恒星目标数，绿色表示与USNO_B1.0星表交叉获得的恒星目标数，黄色为对两个交叉结果取“并集”获得的恒星目标数。很明显，“并集”获得到的恒星目标完备性更高。

图3 星像数量统计图Fig.3 Statistical chart of star quantity

图4 是一帧LUT图像宇宙线认证识别结果（彩图见期刊电子版）。图中，红色圆圈标记宇宙线事件，蓝色圆圈标记恒星目标。很明显，通过以上处理流程能够将较亮的宇宙线与恒星目标很好地识别区分开来。

图4 宇宙线识别结果示意图Fig.4 Schematic diagram of cosmic ray recognition re⁃sults

3.6 算法先进性验证

3.6.1 与恒星星像类似的宇宙线识别效果检验

在天文图像宇宙线识别中，最关键的是剔除与恒星目标类似的宇宙线事件。

针对从所有LUT图像中提取出的宇宙线事件和恒星目标，使用IRAF.PSF_measure软件做了图像轮廓拟合，计算出每个样本的半高全宽（Full Width at Half Maximum，FWHM）。

点扩散函数（Point Spread Function，PSF）拟合中，本文采用普遍使用的Moffat函数表示的PSF轮廓强度分布公式［11］，即：

式中：I0为中心亮度值，B为背景值，α值与FWHM相关，β的值决定曲线轮廓。Moffat函数的FWHM为：

PSF拟合所得FWHM的统计结果如图5所示。很明显，宇宙线事件的FWHM明显小于恒星目标，但仍有一部分相重叠。其中，90%的宇宙线样本的FWHM小于1.6 pixel；95%的恒星目标的FWHM为1.5～2.2 pixel，而在这个区间内的宇宙线事件数量占宇宙线事件总数的8.6%。这8.6%的宇宙线事件在轮廓上与恒星目标完全类似，但本算法还是可以将它们识别出来，并和恒星目标区分开来。

图5 半高全宽统计图Fig.5 Histogram of FWHM

3.6.2 宇宙线识别完备性检验

对于天文CCD图像，拉普拉斯算法［3］比直方图算法［4］和万能消除算法［6］更能准确地识别宇宙射线，且探测能力要明显优于两者［12］。本文将基于天文位置定标的宇宙线识别算法和拉普拉斯算法的宇宙线识别结果进行了统计对比。

对230帧LUT序列图像进行数据分析，本文提出的算法共得到29 731例宇宙线事件；拉普拉斯算法使用IRAF.LA.Cosmic软件包通过设置sigclip=4.5，sigfrac=0.5，经过4次迭代运算后，选取相同大于5倍背景噪声亮度阈值，共得到26 750例宇宙线事件。本文基于天文位置定标的识别算法多识别出2 981例宇宙线事件，占拉普拉斯算法识别结果的11.14%。

对两种方法识别到的宇宙线事件进行了位置交叉证认，两种算法识别到共同宇宙线样本26 234例。结果如表2所示，在拉普拉斯算法检测到的样本中，98.07%可被本算法检测到，具有很高重合率。

表2 本文算法与拉普拉斯算法的结果比较Tab.2 Comparison of proposed algorithm and Laplace al⁃gorithm

对本文算法没有识别到而拉普拉斯算法识别到的516例样本进行检查，通过与恒星目标进行位置认证，发现其中的251例为恒星目标，即本算法没有识别到的样本中48.64%为恒星目标，是错误识别。

4 处理结果与分析

4.1 宇宙线事件电子沉积数分析

本文将单例宇宙线事件影响CCD像元上ADU（Analog-Digital Unit）值相加，根据CCD增益值GIAM=1.59（Electrons/ADU）计算出每例宇宙线事件的电子沉积数，取对数后做出了电子沉积数统计图，并与HST探测器图像中宇宙线事件的电子沉积数［13］进行对比，结果如图6所示。

图6 宇宙线样本电子沉积统计图Fig.6 Electron deposition by cosmic rays

通过与HST巡天相机（Advanced Camera for Surveys，ACS）上WFC与高分辨率探测器（High Resolution Channel，HRC）图像中宇宙线电子沉积数分布进行对比，发现分布形态基本类似，但LUT图像中宇宙线事件电子沉积数的峰值要略高于HST的电子沉积数峰值。最可能的原因是，宇宙射线打到CCD靶面上，其电子沉积数会随着CCD探测器敏感层的厚度不同而发生变化［14］。

4.2 CCD靶面上宇宙线事件入射角度分析

宇宙线打到CCD靶面时，留在CCD上的形态主要是由宇宙线粒子的能量和入射角度所决定的。为了量化定义宇宙线粒子对CCD带来的影响，依据粒子射入CCD成像元件敏感层的瞬态效应，原子量大于1的粒子在器件中的轨迹可以认为是一条直线，粒子和CCD作用后产生的电荷信号会呈现集中的光点轮廓［15］。宇宙线与CCD作用的简单形态模型如图7所示。

图7 宇宙线事件和CCD作用示意图Fig.7 Model of cosmic rays hit on CCD

宇宙线事件在CCD上留下的形态轮廓中，定义通过重心位置的最长线段为长轴，由重心指向长轴与宇宙线截取轮廓最远交汇点的方向，就是宇宙线事件在CCD平面上的入射方向，如图7所示的θ角。

图8为一帧图像上识别到的所有宇宙线事件入射角度示意图。图8中，箭头方向为宇宙线事件在CCD平面上的入射方向。可以看出，一帧图像中，CCD靶面上能够接收到来自不同方向的宇宙线事件。

图8 CCD靶面上宇宙线入射方向示意图Fig.8 Incident direction of cosmic rays on CCD

本文随机选取12帧图像，对识别到的1 280例宇宙线事件计算其入射角，并进行了统计分析。统计时平均分成了24个方向，每个方向包含15°区间，结果如图9所示。在入射角为45°～60°，285°～300°上分别有167，211例事件，其他方向上平均为41例，这两个方向上的宇宙线事件发生频次明显多于其他方向。

图9 CCD靶面上宇宙线入射角统计图Fig.9 Statistical chart of incidence angle of cosmic rays on CCD

在这些数据采集期间，宇宙线事件可认为是各向同性的，探测到的宇宙线事件集中在两个方向上，应该归咎于这两个方向上的等效铝厚度较薄。原则上依据辐射屏蔽模型［16］，如果已知嫦娥三号着陆器的物质结构模型，可以仿真计算各个方向上的等效铝厚度。因暂时无法得到嫦娥三号着陆器的物质结构模型，本文在此不做定量分析。

4.3 建立“宇宙线实体样本库”

利用天文位置定标算法对230帧LUT图像中的宇宙线事件进行识别提取，并依据其形态特征进行分类，制作出“宇宙线实体样本库”。此“宇宙线实体样本库”可作为输入来源，为后续空间天文图像中宇宙线的仿真研究提供支持。

对于230帧LUT图像，首先依据每帧CCD图像的背景噪声值，将信号值小于5倍背景噪声的像元赋值为0；然后对于识别出的恒星目标，以恒星目标的重心为中心，将恒星目标大于5倍背景噪声的连通区域赋值为0；最后将连通像元数小于4倍、大于5倍背景噪声的连通区域也赋值为0。经过此过程，CCD图像中就只保留了宇宙线事件的信号值，其他像元值均为0。

以宇宙线的重心位置为中心，在CCD图像中搜索非零值的连通区域，直至此连通区域周围均为0值，定义此连通区域为一例宇宙线事件。从230帧图像中共识别出29 494例宇宙线事件。

对于每一例宇宙线事件，本文截取出包含这一宇宙线事件完整信号的图片，此图片在宇宙线连通区域四边各多出一个像元，因此图片为长方形或方形。截取图片样例如图10所示。

图10 宇宙线事件样例及对应的三维柱状图Fig.10 Cosmic ray samples and their three dimensional bar chart

对于29 494幅截取图片，依据宇宙线事件信号值衰减的类型，将宇宙线事件分为两类：衰减类型单一的单峰宇宙线事件，衰减类型非单一的多峰宇宙线事件。单峰宇宙线事件如图10（a）～10（c）所示；多峰宇宙线事件如图10（d）所示。统计得出，单峰宇宙线事件有26 459例，占事件总数的89.71%；多峰宇宙线事件仅有3 035例，占事件总数的10.29%。分类结果如表3所示。

表3 宇宙线事件的分类结果Tab.3 Classification results of cosmic rays

最后，将分类完成的宇宙线事件截取图片编号后存储，形成“宇宙线实体样本库”。

5 结 论

在空间天文CCD图像中，宇宙线识别遇到的最大挑战是将形态与星像类似的宇宙线事件与恒星目标区分开来。为了尽可能准确地识别宇宙线事件，本文基于嫦娥三号LUT获得的CCD图像，设计和实现了一套基于天文位置定标的宇宙线识别算法，该算法与常用的拉普拉斯算法相比，检测出事件总数多出11.14%。对利用该算法识别出的宇宙线事件，进一步分析了电子沉积数及入射角的分布情况，发现宇宙线电子沉积分布与HST基本类似，入射角度集中在两个特定方向上，初步判断是由嫦娥三号着陆器内部屏蔽效应造成的。最后，制作出一个包含29 494例宇宙线事件截取图片的空间天文CCD图像“宇宙线实体样本库”，可为空间天文望远镜图像仿真系统中宇宙线事件的仿真提供支持，具 有重要的应用价值。