消除CCD图像中宇宙射线的算法的比较

刘婷婷，彭青玉，3

(1.暨南大学计算机科学系，广东 广州 510632；2.广东省高等学校光电信息与传感技术重点实验室，广东 广州 510632；3.中国科学院光学天文联合开放实验室，云南 昆明 650011)

宇宙射线是来自宇宙空间的各种高能粒子形成的射流，主要包括质子、粒子和少量其他原子核。它们可以看做是天文图像中的盐噪声，明显高于周围像素点的灰度值，并且具有随机分布的规律。部分宇宙射线很有可能位于待测星的星像中，改变了原有像素点的信息，造成星像中心位置的偏移和光度测量值的偏差。因此，消除宇宙射线噪声对于天文图像的分析是非常重要的。

早期最简单直接的方法就是对同一视场拍摄多幅图像，然后将那些好的图像的像素点替代被宇宙射线污染的像素点。目前这一类的算法已发展到很成熟的阶段，其中具有代表性的有Windhorst等人的方法[1]。然而这些方法都依赖于对同一视场拍摄多幅图像。

众所周知，中值滤波方法[2]对去除随机噪声有良好的效果。然而该方法在对CCD图像所有像素点进行滤波的同时，也平滑了星像的边缘，改变了星像中的有用信息。可见中值滤波方法不适用于高精度测量的CCD图像处理。

因此，基于单幅图像来消除宇宙射线更具有挑战性。目前有基于单幅图像消除宇宙射线方法[3]及线性滤波方法[4]。针对这一问题本文比较了3种较新的消除宇宙射线的算法：Laplacian边缘检测算法[5](简称Laplacian算法)，基于直方图的快速算法[6](简称直方图算法)以及万能噪声消除算法[7](简称万能算法)。采用云南天文台1m望远镜拍摄的多幅图像，研究星像和星系中人工及IRAF[8]软件添加宇宙射线及其识别和剔除，通过编程实现3种算法来处理这些图像，并分析比较每种算法的优缺点。Laplacian算法被广泛用于剔除宇宙射线，Van[5]指出该方法非常有效，但也认为这种算法比较复杂，处理速度较慢。而直方图方法由于算法简单因此处理图像的速度比较快。万能方法声称可以处理任何噪声，将该算法用于CCD图像宇宙射线的消除，并与前述两种算法进行分析比较。为了进一步探讨如何准确地替代宇宙射线像素点，将这3种算法分别用于处理实拍图像，并采用两类方法来替代宇宙射线像素点的灰度值：中值滤波方法、曲面拟合方法。

本文第1节将介绍上述3种算法的原理；第2节给出程序的设计与实现；第3节将通过对不同的CCD图像进行处理，给出实现过程和结果，比较3种算法的优缺点，同时探讨宇宙射线像素点的替代方法；第4节将对各种算法进行综合比较并得出结论。

1 算法介绍

1.1 Laplacian算法

Laplacian算法是通过构造和设置两个阈值来识别宇宙射线。它首先通过原始图像I进行子采样放大(子采样因子为fs)，放大后的图像为I(2)。然后与Laplacian模板做卷积，如式(1)。▽2f为Laplacian算子，该算子是线性二次微分算子，用来获取有大的灰度变化的图像边缘信息，当然包括了宇宙射线影响的部分。之后将图像L恢复到原始图像分辨率得到图像L+。

(1)

另一方面，对原始图像I用来构造剔除宇宙射线的噪声模型，如式(2)。该模型仅由泊松噪声和读出噪声组成，不包括宇宙射线噪声。其中g是以光电子/ADU为单位的增益因子，σrn是以电子为单位的读出噪声，M5是5×5的中值滤波模板。卷积运算M5*I得到的是有用信号的泊松噪声平方的估计。

(2)

至此可以得到噪声的比率S，如式(3)，其中fs是子采样因子。这个比率S反映了图像中每个像素点含有噪声的比率，像素点噪声越大则这个比值越大。然后设置阈值σlim，认为S值大于σlim的那些像素点为候选的宇宙射线。

(3)

根据星像有对称性而宇宙射线却没有的性质，可以进一步区别宇宙射线与星像，从而构造一个精细结构如式(4)。其中M3表示3×3的中值滤波模板，M7表示7×7的中值滤波模板。卷积运算M3*I得到图像的中、低频信息，而(M3*I)*M7得到图像的低频信息，两者相减得到图像的中频信息。容易推断星像所在区域的F值会较大，而宇宙射线所在区域算出的该值会很小。

F=(M3*I)-[(M3*I)*M7].

(4)

因此构造第二个比率，如式(5)。通过设置阈值flim，认为T大于flim的像素点为候选的宇宙射线，

(5)

最后同时设置两个条件S>σlim和T>flim来探测宇宙射线。

该方法仅仅考虑用探测点附近的区域内的中值滤波来替代宇宙射线点的原始灰度，是否有更好的替代算法没有深入讨论。算法在实现过程中可以迭代执行，直到不能探测到新的宇宙射线才终止。

1.2 直方图算法

直方图方法是一种简单且能直接用于识别宇宙射线的算法。该算法首先将图像分成若干个子图，然后分析这些子图的直方图(绝大多数直方图分布是非常紧凑的，宇宙射线则会以孤立的单个点分布在这些紧凑范围之外)，最后根据直方图分布的特性设置阈值从而探测宇宙射线。对于阈值的设定范围，Pych[6]建议为th*σ，th一般约等于3.0。每个子图的σ计算公式，如式(6)，其中ci为子图中的像素灰度值,n为子图中的像素个数。

(6)

宇宙射线有时由多个像素点组成，该算法利用增长半径的方法进一步探测宇宙射线。它是以探测到的宇宙射线点为中心，在一定的半径范围内继续探测是否有其它的宇宙射线点。对探测到的宇宙射线点，通过设定一个范围，然后在这个范围内求出未受宇宙射线影响的所有像素点的平均值，再用该平均值去替代宇宙射线点。

该算法在编程实现过程中需要输入若干个参数，这些参数可以根据图像的不同做调整，最终找到合适的值：

(1)输入子图的长度和宽度x,y(文中的图像选用的x,y都为8);

(2)输入th(文中的值分别为3.2、3.3);

(3)输入探测宇宙射线的增长半径radius(采用radius=1);

(4)输入消除宇宙射线的最大半径和最小半径(分别为2和1)。

1.3 万能算法

万能算法的核心是ROAD(Rank-Ordered Absolute Differences)算法。该算法的思想是通过排序统计像素点与周围像素点之间的绝对值差，从而设置合适的ROAD阈值来探测识别噪声点。为了更加清楚地理解ROAD算法，引用万能算法[7]中的图来具体分析，如图1。

图1 ROAD值的计算Fig.1 Calculation of the values of ROAD

2 程序的设计与实现

采用Visual C++ 6.0在Windows环境中编程[9]实现以上所述的3种算法。其中Laplacian算法和直方图算法具体实现过程见图2。其中图2(a)和图2(b)分别表示这两种算法实现的流程图。万能算法由于算法本身比较简单文中省去其流程图。

3 实验与结果分析

以云南天文台1m望远镜拍摄的图像进行资料分析，实验分别采用Laplacian算法、直方图算法和万能算法来消除宇宙射线，重点研究星像和星系上人工及利用IRAF软件添加宇宙射线的图像。

根据上文提到的宇宙射线的性质，在图像的星像和星系上人工添加宇宙射线点，该宇宙射线点的灰度值要大于m+3N，其中m为该像素点的原始灰度值，N为公式(2)中的噪声。

3.1 Laplacian算法处理星像和星系上的宇宙射线

Laplacian算法分别用于处理星像和星系上的宇宙射线，如图3、4。图3(a)表示在原始星像上坐标为(1005,143)的位置，人工添加灰度值为1426的宇宙射线(并用黑色点表示以示区别)。星像周围的亮点为原图中带有的宇宙射线点。通过设置两个阈值分别为σlim=0.7、flim=2且经过2次迭代运算后，可以准确识别到星像上和星像周围的宇宙射线，如图3(b)，并没有产生假的探测。图3(c)是用3×3中值滤波的方法剔除噪声后的图像，显然，剔除宇宙射线噪声后的视觉效果良好。

(a)Laplacian算法实现的流程图 (b)直方图算法实现的流程图图2 两种算法实现的流程图Fig.2 Flowcharts of the two algorithms

(a) (b) (c)图3 Laplacian方法剔除星像中的宇宙射线Fig.3 Rejection of cosmic rays in an image of a star with the Laplacian algorithm

图4(a)为在原始图像的星系上添加宇宙射线点，该点的坐标为(1031，1185)，灰度值为1424。图4(b)显示了经过该算法处理后，探测到星系上的噪声点。图4(c)为消除宇宙射线后的图像，图像的视觉效果也很好。针对这幅图像该算法设置的两个阈值分别为σlim=0.5、flim=1.3。

(a) (b) (c)图4 Laplacian方法剔除星系上宇宙射线Fig.4 Rejection of cosmic rays in an image of a galaxy with the Laplacian algorithm

3.2 直方图算法处理星像和星系上的宇宙射线

用该方法分别处理星像和星系上的宇宙射线，如图5、6。其中图5(b)的星像上的宇宙射线能被探测到，但星像周围的噪声并没有被识别。由图5(c)清楚看出星像周围留有未消除的宇宙射线。采用的参数分别为：子图长宽都为8，th取值为3.3，增长半径为1，消除宇宙射线的最大、最小半径分别为2、1。

(a) (b) (c)图5 直方图方法剔除星像中的宇宙射线Fig.5 Rejection of cosmic rays in an image of a star with the histogram algorithm

(a) (b) (c)图6 直方图的方法剔除星系上宇宙射线Fig.6 Rejection of cosmic rays in an image of a star with the histogram algorithm

使用该方法处理图4(a)，剔除噪声的过程，如图6。与图4相比较，这两种方法都能精确探测到星系上的宇宙射线。这是由于星系上像素点的灰度值变化比星像上灰度值变化要缓慢，因此星系上的宇宙射线比较容易探测到。该图像采用的参数分别为：子图长宽为8，th取值为3.2，增长半径为1，消除宇宙射线的最大、最小半径分别为2、1。

通过实验比较得出，Laplacian算法比直方图算法更能准确地识别和探测宇宙射线，且探测能力要明显优于直方图算法。但从处理速度方面考虑，直方图算法要快很多。

3.3 万能算法处理星像和星系上的宇宙射线

实验发现对于在星像上添加宇宙射线的天文图像，由于宇宙射线的ROAD值可以很大也可以很小，而星像上未受宇宙射线影响的像素点的ROAD值同样具有此特点。如图7，两者ROAD值非常相似。因此万能算法仅仅通过设定ROAD阈值，是较难区分星像和宇宙射线的。

万能算法用于处理星系上的宇宙射线，如图8，该方法通过设置ROAD阈值，可以精确剔除星系上的宇宙射线。

图7 万能算法用于星像上有宇宙射线的情况 ：B表示未受宇宙射线影响的背景区域，X表示未受宇宙射线影响的星像区域，A1，A2，A3，A4，A5表示背景上的宇宙射线，Y表示星像上的宇宙射线。图中的表给出了这些区域上像素点的ROAD 值

(a) (b) (c)图8 万能消除方法剔除星系上宇宙射线Fig.8 Rejection of cosmic rays in an image of a galaxy with the universal noise removal algorithm

3.4 IRAF软件添加宇宙射线

将IRAF[8]软件用于图像中添加宇宙射线，该图像是2009年2月25日晚由云南天文台1m望远镜所拍摄。如图9(a)为无宇宙射线的原始图像，图9(b)为IRAF软件添加的宇宙射线，图9(c)为添加完噪声后的图像。图10、11、12分别为用Laplacian方法、直方图方法、万能消除方法剔除宇宙射线。其中图10(b)、11(b)、12(b)分别为3种方法识别的宇宙射线，与图9(b)比较得出，Laplacian方法几乎识别到所有的宇宙射线，万能消除方法看上去也都识别到所有的宇宙射线，但还是有部分噪声没有探测到，且存在小部分假的探测，而直方图方法则有更多的宇宙射线未识别。图10(c)、11(c)、12(c)分别为3种方法消除宇宙射线后的图像，且与图9(a)比较，其中图10(c)最接近图9(a)。相比较这3种方法，Laplacian方法识别和剔除宇宙射线的效果较好。

为了进一步研究宇宙射线对星像位置和光度的影响，重点讨论利用实拍图像中受宇宙射线影响的星像作为研究对象，对其位置和光度进行分析。

(a)无宇宙射线的原始图像 (b)使用IRAF添加的宇宙射线 (c)添加宇宙射线后的图像图9 IRAF软件添加宇宙射线Fig.9 Simulation of cosmic-ray hits with the IRAF

(a)添加宇宙射线后的图像 (b)Laplacian方法识别宇宙射线 (c)消除宇宙射线后的图像图10 Laplacian方法剔除宇宙射线Fig.10 Rejection of cosmic rays with the Laplacian algorithm

(a)添加宇宙射线后的图像 (b)直方图方法识别宇宙射线 (c)消除宇宙射线后的图像图11 直方图方法剔除宇宙射线Fig.11 Rejection of cosmic rays with the histogram algorithm

(a) 添加宇宙射线后的图像 (b)万能消除方法识别宇宙射线 (c)消除宇宙射线后的图像图12 万能消除方法剔除宇宙射线Fig.12 Rejection of cosmic rays with the universal noise removal algorithm

3.5 剔除实拍图像中星像上的宇宙射线

将上述3种方法分别用于处理实际拍摄中星像受到宇宙射线影响的图像。该图像是2009年2月15日晚由云南天文台1m望远镜所拍摄。如图13、14、15，分别为Laplacian算法、直方图算法、万能消除算法处理该图，图13(a)中的星像由于受到宇宙射线的影响，在测量该星像位置时，测量结果有严重的偏差。图13(b)、14(b)、15(b)分别为3种算法识别的宇宙射线，其中区域A、B为星像上的宇宙射线，区域C为背景上的宇宙射线。图13(c)、14(c)、15(c)为用下述替代方法替代宇宙射线后的图像。

(a) (b) (c)图13 Laplacian算法剔除实拍图像中星像上的宇宙射线Fig.13 Rejection of real cosmic rays in an actual image with the Laplacian algorithm

(a) (b) (c)

(a) (b) (c)图15 万能算法剔除实拍图像中星像上的宇宙射线Fig.15 Rejection of real cosmic rays in an actual image with the universal noise removal algorithm

在消除宇宙射线噪声之后，重新测量该星像的位置时，发现测量结果有明显的改进。改进的好坏取决于如何替代宇宙射线像素点。采用的3种替代宇宙射线的算法分别为：3×3像素范围的中值滤波法、5×5像素范围的曲面拟合法和7×7像素范围的曲面拟合法。曲面拟合公式如式(7)所示，其中x、y为像素点的坐标位置，参数a、b、c、d、e、f为待求的拟合参数。

f(x,y)=ax2+bxy+cy2+dx+ey+f

(7)

在连续曝光的一系列图像，给定的未受宇宙射线影响的两颗星之间的距离以像素为单位与光度差随时间的变化过程是缓变的。如果其中有一颗星受到噪声的影响，该缓变过程将会有突变。其中星像位置的测量采用高斯拟合[10]的方法，星像光度的测量采用孔径测光的方法。

图16 星像位置随时间(北京时)的变化Fig.16 Different positions of a star at different moments(in the Beijing time)

下面讨论宇宙射线对星像位置及光度的影响，图16为图13(a)、13(c)、14(c)、15(c)中星像与其相邻的一颗未受宇宙射线影响的亮星(图中未给出)之间的距离随时间变化的情况。可以看出，3种剔除宇宙射线的方法在20h35m42s(北京时)这个时刻符合较好，且比受到宇宙射线影响的图像的位置有较大改善(约0.16pixel)。图17为图13(a)、13(c)、14(c)、15(c)中星像与其相邻的一颗未受宇宙射线影响的亮星(图中未给出)之间的光度差随时间变化的情况。可以看出，3种剔除宇宙射线的方法在20h35m42s(北京时)这个时刻符合较好，且比受到宇宙射线影响的图像的光度有较大提高(约0.04mag)，基本上达到预期目标。

图17 星像光度随时间(北京时)的变化Fig.17 Different measured magnitudes of a star at different moments(in the Beijing time)

4 结 论

对于星系上受宇宙射线影响的CCD图像，Laplacian算法、直方图算法、万能算法都能很好地用于识别宇宙射线，对于星像上受宇宙射线影响的CCD图像，Laplacian算法较其它两种算法更能精准地识别噪声。从算法处理速度考虑，直方图方法和万能方法由于算法简单因此处理速度更快。通过实验发现，宇宙射线的替代算法是值得深入研究的。

致谢:感谢暨南大学计算机科学系张庆丰老师、孟小华老师、李展老师为本文提出建设性建议。

[1]Windhorst R A , Franklin B E, Neuschaefer L W. Cosmic rays in multi-orbit images with the HST wide field planetary camera 2[J]. PASP, 1994, 106：798.

[2]Gonzalez R C, Woods R E, 著，阮秋琦，阮宇智 等译.数字图像处理[M]. 北京:电子工业出版社，2007.

[3]Zhu Z Q, Ye Z F. Detection of cosmic-ray hits for single spectroscopic CCD images[J]. PASP, 2008, 120(869):814-820.

[4]Rhoads J E. Cosmic-ray rejection by linear filtering of single images[J]. PASP, 2000, 112(771): 703-710.

[5]Van Dokkum P G. Cosmic-ray rejection by Laplacian edge detection[J]. PASP, 2001, 113(789):1420-1429.

[6]Pych W. A fast algorithm for cosmic-ray removal from single images[J]. PASP, 2004, 116(816):148-153.

[7]Garnett R , Huegerich T , Chui C, et al. A universal noise removal algorithm with an impulse detector[J].IEEE T Image Processing, 2005,14(11):1747-1754.

[8]Wells L A , Bell D J. Cleaning images of bad pixels and cosmic rays using IRAF.1994.

[9]杨淑莹. VC++图像处理程序设计[M].北京:清华大学出版社, 2007.

[10]李展, 彭青玉, 韩国强. CCD图像数字定心算法的比较[J]. 天文学报,2009,50:340-348.

LI Zhan,PEIVG Qing-yu,HAN Guo-qiang.Comparison of Digital Centering Algorithms Based on CCD Images[J].AAS,2009,50(3):340-348.