消除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.