冯海霞，陈建军，邓建榕，赵永恒

(1. 中国科学院国家天文台，北京 100101；2. 中国科学院大学，北京 100049；3. 中国科学院光学天文重点实验室 (国家天文台)，北京 100101)

宇宙线是来自外太空的高能粒子，1912年，奥地利地理学家赫斯在电离室测量电流时第一次发现了宇宙线。宇宙线中质子和氦核占98%，重核1%，电子1%，微量其他粒子[1]。当初级宇宙线进入地球大气层时，质子和氦核等高能粒子与大气中的原子核相互作用产生级联的广延大气簇射，生成大量的次级粒子。在使用电荷耦合器件观测目标源时，这些粒子在CCD上留下不同的径迹，其中包括电子、α粒子、X射线以及μ子等。靠近大气底部的宇宙线主要是簇射产生的次级μ子。在CCD图像的宇宙线中，占据主要成分的是μ子。μ子和电子等带电粒子在CCD上留下明显的特征径迹。电子大部分是由环境中放射性同位素产生的，并非来源于宇宙线[2]。α粒子和X射线等其他粒子在CCD图像中呈点状，不易与CCD芯片本身产生的电子噪声分开。因此，研究CCD图像中的宇宙线，需要对宇宙线中的μ子进行甄选。

当前探测宇宙线的装置很多，例如西藏羊八井国际宇宙线观测站、印度乌蒂的GRAPES-3、阿根廷的皮埃尔·俄歇观测站(Pierre Auger Observatory)以及四川的高海拔宇宙线观测站等。这些装置采用一系列探测器阵列和望远镜来探测宇宙线粒子。本文利用大天区面积多目标光纤光谱天文望远镜(the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST，又叫郭守敬望远镜)探测μ子。文[3]系统介绍了郭守敬望远镜的结构、原理及应用。郭守敬望远镜配置有16个光谱仪和32台4 k × 4 k的CCD相机。CCD芯片大小为49.2 mm × 49.2 mm，像素大小为12 μm × 12 μm，工作温度在-100 ℃左右。从2011年先导巡天至今，已经积累了一百多万张原始图像，每幅图像都有很多宇宙线事件，可以用来研究宇宙线的性质和分布。图像中的宇宙线个数与图像的曝光时间有关，将曝光时间归一化到分钟。在单位时间内，图像中的宇宙线事件个数为5～60。

CCD的工作原理是利用光电效应将光子转化为电子，电容器吸收释放的电子，然后将收集的电子转化为电压，最后以数字信号的形式输出。μ子有带负电和带正电两种，电荷数为1。它的质量约为106 MeV，比电子重207倍。当它穿过CCD时，可以电离径迹上的硅原子，受到库伦散射和轫致辐射的影响比电子小，在CCD上留下可识别的直线径迹。

文[2]介绍CCD图像中的宇宙线事件关于质心呈现不同的分布，结合宇宙线事件的形态特征甄选μ子。本文基于郭守敬望远镜的观测数据，对图像进行分析处理，根据μ子在图像中的形态特征对其进行甄选。首先从图像中提取宇宙线像素，采用拉普拉斯边缘检测法从原始图像中提取宇宙线候选像素列表。在提取的宇宙线候选像素列表中包含一些坏像素和噪点，需要对这些坏像素和噪点进行处理。数据处理过程中，为了得到独立完整的宇宙线事件，采用凝聚层次聚类算法将宇宙线像素聚类成宇宙线事件。提取宇宙线事件的特征，结合μ子在图像中的形态特征甄选出宇宙线μ子。最后对甄选技术进行评估。

1 拉普拉斯边缘检测法提取宇宙线候选像素列表

在郭守敬望远镜的观测数据中，单次曝光的CCD图像中含有上百条流量值很强的光纤光谱。这些图像中，不仅包含所要观测的目标信息，同时还包括一些噪声。在对图像中的μ子进行甄选时，首先从图像中获取宇宙线候选像素列表。检测图像宇宙线的方法有很多，例如经典的中值滤波法[4]、拉普拉斯边缘检测法[5]、基于直方图的快速算法[6]以及万能噪声消除算法[7]等。本文采用改进的拉普拉斯边缘检测法[8]提取图像中的宇宙线候选像素列表。相比于其他方法，该方法在图像背景流量值很强的情况下，可以有效地从单次曝光的多光纤光谱中检测宇宙线。它是在文[5]算法的基础上进行改进,使得当宇宙线落在光纤光谱轮廓内时可以被识别出来。

改进的拉普拉斯边缘检测算法的流程：

(1)对原图像进行预处理，然后与尺寸大小为3 × 3的拉普拉斯模板进行卷积运算得到图像I；结合泊松噪声和读出噪声，在原图像上构造剔除宇宙线的噪声模型N；再根据I和N的比率关系，生成一个原始的宇宙线候选像素列表；

(2)在原图像的基础上去除步骤(1)产生的宇宙线候选像素列表，生成一个新的图像。对图像中的每一条光纤进行处理，将距离光纤轨迹中心8个像素以内的像素点组成一个模块。结合光谱在空间方向和色散方向的分布情况，拟合这些模块的图像轮廓(这些拟合的图像轮廓不包括宇宙线)，得到拟合后的图像I1；

(3)根据泊松噪声和读出噪声，对拟合后的图像I1构造剔除宇宙线的噪声模型N1；最后将原图像和拟合后图像的残差与噪声模型N1进行比较，生成宇宙线候选像素列表。

采用拉普拉斯边缘检测法提取宇宙线候选像素列表的结果见图1。图1(a)是从CCD图像中截取的一部分原始图像；图1(b)对应于图1(a)采用拉普拉斯边缘检测法提取后的宇宙线像素，图中的背景为黑色，对应的流量值为0。从图1可以直观地看出，该算法可以有效地从图像中提取宇宙线像素。

2 坏像素和噪点的处理

在使用拉普拉斯边缘检测法提取的宇宙线候选像素列表中，包括宇宙线、CCD上的坏像素以及算法导致的噪点。CCD上的坏像素是本身的坏点，基本在同一批拍摄的图像中的同一位置出现。可以利用同一目标的多次曝光图像标识这些坏像素，然后将它们去除。在使用拉普拉斯边缘检测法提取宇宙线候选像素列表时，如果图像中的光纤轨迹中心发生偏移，会在拟合图像轮廓时产生较大的误差，最终在生成宇宙线列表的基础上产生拟合的偏差，这些偏差就是算法导致的噪点。这些噪点在图像中呈现几列虚线，虚线处的像素点个数远超过宇宙线在图像中每一列可能的像素点个数(虚线处的像素点个数为150～500，而图像中每一列的宇宙线像素点个数为8～50)。利用这一特点，在图像中找出异常的列来标注算法导致的噪点，然后将它们从图像中去除。

图1 拉普拉斯边缘检测法提取宇宙线

Fig.1 Extraction of cosmic rays with Laplacian edge detection algorithm

3 凝聚层次聚类算法聚类宇宙线事件

从图像中提取宇宙线像素后，采用凝聚层次聚类算法将宇宙线像素聚类成一个个独立的宇宙线事件。本文采用改进的经典凝聚层次聚类算法。

层次聚类试图在不同层次对数据集进行划分，从而形成树形的聚类结构。数据集的划分可采用自底向上的聚合策略，也可采用自顶向下的分拆策略[9]。凝聚层次聚类是层次聚类中的一种，它对数据集的划分采用自底向上的聚合策略。主要思路是将数据集中的每一个样本看成单个的聚类簇，然后计算每个聚类簇之间的距离，找出其中距离最近的两个聚类簇将它们合并，重复上述距离计算和合并的过程，直至达到预设的聚类簇个数。

在经典的凝聚层次聚类算法中，算法的关键是计算类间距的方法和确定聚类簇的个数。每一个聚类簇是一个样本集合，可以采用计算集合的某种距离来计算类间距，常见的距离有最小距离、最大距离和平均距离。将宇宙线像素聚类成宇宙线事件，使得宇宙线事件之间相互独立，采用最小距离计算类间距。最小距离计算公式：

(1)

其中，Ci，Cj为给定的聚类簇；{ma}a=1,2, ...,k为Ci的样本点集；{mb}b=1,2, ...,l为Cj的样本点集；dist为计算ma,mb样本点的距离函数。

在处理图像时，无法预知图像中宇宙线事件的个数，所以无法预设聚类簇的个数。需要设定类间距dlim作为终止运算的条件。具体算法如下：

(1)假设数据集为D，数据集中的数据是图像中宇宙线所对应的像素，将每个像素看作一个初始的聚类簇分别对应C1,C2, ...,Cn。

D={C1,C2, ...,Cn} .

(2)

(2)计算所有聚类簇之间的类间距，如果两个聚类簇各自只含有一个样本点时，类间距是这两个样本点的欧氏距离；反之，类间距是两个聚类簇之间的最小距离。如果两个聚类簇的类间距不大于设定的阈值dlim，将这两个聚类簇合并成一个聚类簇。采用欧氏距离和最小距离计算类间距：

(3)

其中，Ci,Cj为数据集中的任意两个聚类簇；{ma}a=1,2, ...,k为Ci的样本点集；{mb}b=1,2, ...,l为Cj的样本点集；xa和ya对应ma在图像上的位置信息；xb和yb对应mb在图像上的位置信息。

(3)更新数据集。

(4)重复步骤(2)、(3)，直到数据集中所有聚类簇的类间距大于设定的阈值，结束运算。

当图像中出现大量宇宙线事件时，需要逐个计算聚类簇之间的距离，计算量剧增，运行速率减慢。结合宇宙线事件是有计数的连续像素这一特点，通过判断与聚类簇最外层元素距离为dlim的范围内是否存在其他的聚类簇优化算法，如果存在其他聚类簇，将其他聚类簇与该聚类簇合并；如果不存在，将数据集中未处理的聚类簇重复上述判断过程，直到所有聚类簇的类间距大于dlim，结束运算。使用该方法将宇宙线像素聚类成宇宙线事件时，如果两个宇宙线事件的最小距离大于dlim，说明两个宇宙线事件相互独立；如果最小距离不大于dlim，说明是同一个宇宙线事件。在提取宇宙线像素的过程中，当宇宙线流量值和背景值比较接近时，会出现提取的宇宙线像素在图像中不连续。尽管这种情况很少出现，为减小宇宙线像素缺失带来的误差，结合宇宙线在图像中呈现随机稀疏分布的特性，设定dlim为3对提取后的宇宙线像素聚类，可以有效地聚类宇宙线事件。

4 特征提取

将图像中的宇宙线像素聚类成宇宙线事件后，结合μ子在图像中的形态特征，对宇宙线事件进行特征提取。在特征提取的过程中，主要使用主成分分析和皮尔森(Pearson)相关系数。通过主成分分析获取宇宙线事件的特征值，通过皮尔森相关系数判断宇宙线事件是否呈线性，以此甄选宇宙线事件中的μ子。

4.1 主成分分析

主成分分析也称为Karhunen-Loeve扩展，是一种经典的特征提取和数据表示技术[10]。文[11]对主成分分析进行了总的概括。主成分分析也可以用于对数据进行降维处理，目的是找一组最优的标准正交向量基，使得数据投影到这组最优标准正交向量基上的方差最大。宇宙线事件是二维数据，假设在一个宇宙线事件中含有n个像素，其样本集M={m1,m2, ...,mn}；将所有样本进行中心化：

(4)

定义样本集的散度矩阵为St：

(5)

(6)

假设从高维空间投影到低维空间的投影矩阵为w，使投影后样本点方差最大化的优化目标为

Stw=λw.

(7)

对St做特征值分解得到特征值λ，将特征值λ进行降序排序，选取最大的特征值，最大的特征值所对应的特征向量就是主成分的解，从而找到最优的标准正交向量基，得到宇宙线事件的特征向量。将宇宙线事件在特征向量方向上的投影作为宇宙线事件的特征。两个特征向量对应于两个特征值。其中一个特征值对应宇宙线的长轴，即在CCD上宇宙线直线径迹所对应的轴；另一个特征值对应宇宙线的短轴，即宇宙线在图像上的宽度。μ子在图像上的径迹是直线型，且径迹长度无法确定，但宽度有限，是个小量。选取标记宽度的特征作为μ子的特征，并将其记为投影宽度。

4.2 皮尔森相关系数

相关系数是用来衡量变量之间线性相关程度的统计量。定义相关系数的方式有很多，最常用的是皮尔森相关系数，皮尔森相关系数用于描述两个变量之间是否存在线性关系。文[12]对皮尔森相关系数进行了详细介绍。假设两个变量分别是x和y，它们的皮尔森相关系数用r表示，则有

(8)

其中，μx，μy分别为变量x，y的均值。r的取值范围为[-1, 1]，|r|越接近1时，两个变量之间的线性相关性越显著。当r> 0时，两个变量之间呈正相关；当r< 0时，两个变量之间呈负相关；当r=0时，两个变量之间彼此独立。

宇宙线μ子在CCD上呈直线轨迹，则μ子在图像中的横、纵坐标存在线性关系，利用皮尔森相关系数标记μ子的这一特征。由于μ子在图像中的位置有时呈现正相关，有时呈现负相关，为了显示线性相关性质的显著程度，在传统的皮尔森相关系数的基础上取绝对值。下文提到的r为皮尔森相关系数的绝对值。

5 实验结果及分析

利用宇宙线μ子与其他粒子在图像中的形态特征差异甄选μ子。在图像中最常见的粒子是μ子和 “worms”(“worms” 是多级散射的电子，详细介绍见文[2])。α粒子在图像上极少出现。CCD图像中还存在一些形状奇特的宇宙线事件，因无法确定其形态特征，难以判断属于何种宇宙线粒子，本文未做讨论。图2是从郭守敬望远镜获得的图像，图2(a)是宇宙线中的μ子，图2(b)是宇宙线中的 “worm”， “worm” 形状呈现一定的弧度。

图2 郭守敬望远镜CCD图像的μ子和 “worm”。(a) 宇宙线μ子；(b) 宇宙线 “worms”

Fig.2 Cosmic-rays muon and “worm” in the CCD image of LAMOST. (a) Cosmic-ray muon; (b) Cosmic-ray “worms”

利用上文介绍的处理过程，从2015年郭守敬望远镜观测的图像中随机选取200多张，得到3万多个宇宙线事件，将选取的宇宙线事件的位置索引信息作为样本集C1。宇宙线事件的位置索引信息可以反映宇宙线事件在图像上呈现的形态。从样本集中通过人眼识别的方式，选取1 000个μ子和1 000个 “worms” 作为数据样本1。统计这些宇宙线事件的投影宽度和皮尔森相关系数，结果见表1。表1中宇宙线事件投影宽度的均值记为μwp，标准差记为σwp；皮尔森相关系数的均值记为μr，标准差记为σr。

表1 关于投影宽度和皮尔森相关系数，宇宙线μ子、“worms” 的特征对比

Table 1 Feature comparison between muons and “worms” in cosmic rays of the “Projection width” and the Pearson correlation coefficient

Type of cosmic raysThe projection width (μwp±σwp)Pearson correlation coefficient (μr±σr)muons1.43 ± 0.500.93 ± 0.06“worms”5.10 ± 2.030.73 ± 0.15

根据表1可以直观地看出，μ子和 “worms” 的区别。在投影宽度中，μ子的特征值小，而 “worms” 的特征值大；在皮尔森相关系数中，由于μ子在CCD中呈直线径迹，对应的横纵坐标的线性相关程度更显著，对应的皮尔森相关系数值更接近1，而 “worms” 由于在CCD中发生多次散射，径迹有一定的弧度，对应的皮尔森相关系数值较μ子的更小。

通过人眼识别的方式从样本集C1中随机选取2 500个μ子和2 500个非μ子的其他宇宙线事件(包含 “worms”、X射线、可能的α粒子等其他粒子)作为数据样本2。这5 000个宇宙线事件的投影宽度和皮尔森相关系数的分布结果如图3(a)。在图3(a)中，利用投影宽度和皮尔森相关系数两个特征可以将绝大多数的μ子和其他宇宙线区分开。在图中，皮尔森相关系数大于0.8和投影宽度小于1的范围内，μ子和其他宇宙线事件不能很好地区分。观察这部分与μ子特征值相近的宇宙线事件，发现这些宇宙线事件的像素个数很少。分析皮尔森相关系数和宇宙线事件像素个数的关联，二者的分布结果见图3(b)。图3(a)中不易区分的μ子和其他宇宙线事件，在图3(b)中可以通过宇宙线事件的像素个数有效地区分。

图3 不同特征值的分布结果。(a) 投影宽度和皮尔森相关系数的分布结果；(b) 皮尔森相关系数和像素点个数的分布

Fig.3 The distribution of different characteristic values. (a) Distribution of projection width vs pearson correlation coefficient;(b) Distribution of pearson correlation coefficient vs the number of pixel

对于像素个数少的宇宙线事件，只结合投影宽度和皮尔森相关系数甄选μ子会有一定的误差。当甄选μ子时，需要结合宇宙线事件的像素个数。图4显示了数据样本2中宇宙线事件像素个数的统计结果。横坐标表示宇宙线事件的像素点个数，纵坐标表示像素点个数对应的宇宙线事件的个数。图4(a)是整个数据集的统计结果。图中宇宙线事件的像素点个数越多，对应的宇宙线事件个数越少。基于图4(a)，选取像素点个数不大于30的宇宙线事件进行统计，结果见图4(b)。图中，宇宙线μ子的像素点个数大于8，其他宇宙线事件的像素点个数大于5。在非μ子的其他宇宙线中，将近一半的宇宙线的像素点个数小于8。在像素点个数大于8的情况下，μ子的个数比其他宇宙线的个数多。因此，将像素点个数设置为8，对宇宙线中μ子和其他宇宙线事件进行分类。

图4 宇宙线事件的像素点个数的统计结果

根据宇宙线事件特征值的不同分布，设定不同的特征值参数，对μ子进行甄选。从样本集C1中，随机选取2 500个μ子和2 500个非μ子的宇宙线在CCD上的位置索引信息作为测试数据集C2，设定不同的特征值参数对宇宙线中的μ子进行甄选，最后的甄选结果如表2。表2显示了甄选结果的准确性，表中第1列是特征参数的设定，其中，μmw为μ子投影宽度的均值，σmw为μ子投影宽度的标准差，μmr为μ子皮尔森相关系数的均值，σmr为μ子皮尔森相关系数的标准差；第2列是正确甄选宇宙线中μ子的个数；第3列是被错分成μ子的宇宙线事件个数；第4列是甄选结果的正确性；最后一列是甄选结果的灵敏度，用于衡量算法对μ子的识别能力。假设在整个数据集中，μ子的个数为M，非μ子的宇宙线事件个数为F。在甄选结果中，被正确甄选为μ子的宇宙线事件个数为TM，被正确甄选为非μ子的宇宙线事件个数为TF，则甄选结果的准确率公式：

(9)

甄选结果的灵敏度公式：

(10)

从表2可以看出，当对投影宽度、皮尔森相关系数和宇宙线事件的像素个数设定不同阈值时，得到不同的准确率；当特征值设定合适的阈值时，可以有效地从宇宙线事件中甄选μ子。表2特征参数的选择需要结合表1中μ子不同特征的均值和标准差。通过统计宇宙线μ子的投影宽度和皮尔森相关系数，得到描述不同特征的具体数值参数以及参数的波动情况。结合不同特征的参数波动，对测试数据集中的μ子进行甄选，得到最终的结果。

当对宇宙线特征设置合适的阈值时，得到的准确率约为98%，说明该方法可以有效地甄选宇宙线中的μ子。该方法对μ子的识别度很高，但是会将与μ子特征相似、不是μ子的宇宙线事件错分成μ子。

尽管该算法在数据处理过程中存在误差，但是通过分析μ子的本征特征，结合数据集中μ子和其他宇宙线事件的形态差异，依然可以有效地甄选宇宙线中的μ子。该方法的核心是特征参数的设定。如果特征值的参数设定恰当，最后的甄选结果更可观，也可以通过训练大量数据集的方式调整特征参数，优化方法。根据表2中最优的甄选条件(表2第2行对应的特征参数)，对郭守敬望远镜6号红端的CCD相机在2017年11月拍摄的图像进行处理。统计这一个月曝光时间为600.0 s的图像中μ子的个数。统计结果见图5。图中横坐标为拍摄CCD图像的修正儒略分钟数，纵坐标为图像中μ子的个数。图中蓝色竖线对应于第二横坐标拍摄CCD图像的日期。图5直观地显示出宇宙线μ子在这一个月的涨落情况。

表2 设定不同特征参数时μ子的甄选结果

图5 2017年11月郭守敬望远镜6号红端CCD图像中宇宙线μ子个数的统计结果

6 总 结

本文探讨了一种简单有效地对单次曝光CCD图像中的μ子进行甄选的方法，结合宇宙线在图像上的形态特征，达到高效处理数据的效果。该方法使用拉普拉斯边缘检测法提取宇宙线像素，采用聚类算法将图像中的宇宙线像素聚类成宇宙线事件，聚类得到宇宙线事件后，通过特征提取和宇宙线事件之间的形态特征差异对μ子进行甄选。在这个过程中，甄选准确率依赖于特征参数的设定，可以通过不断的统计和测试，设定更优的参数以及增加提取的特征，提高甄选结果的准确性。

在获取训练和测试样本集的过程中，通过人眼识别的方式区分宇宙线μ子和其他宇宙线事件，会导致选择的数据样本集有一定的偏差，影响之后的统计和甄选结果。

致谢：感谢中国科学院高能物理研究所胡红波老师的有益讨论。