劉婷婷，彭青玉，3

(1.暨南大學(xué)計(jì)算機(jī)科學(xué)系，廣東 廣州 510632；2.廣東省高等學(xué)校光電信息與傳感技術(shù)重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510632；3.中國科學(xué)院光學(xué)天文聯(lián)合開放實(shí)驗(yàn)室，云南 昆明 650011)

宇宙射線是來自宇宙空間的各種高能粒子形成的射流，主要包括質(zhì)子、粒子和少量其他原子核。它們可以看做是天文圖像中的鹽噪聲，明顯高于周圍像素點(diǎn)的灰度值，并且具有隨機(jī)分布的規(guī)律。部分宇宙射線很有可能位于待測星的星像中，改變了原有像素點(diǎn)的信息，造成星像中心位置的偏移和光度測量值的偏差。因此，消除宇宙射線噪聲對于天文圖像的分析是非常重要的。

早期最簡單直接的方法就是對同一視場拍攝多幅圖像，然后將那些好的圖像的像素點(diǎn)替代被宇宙射線污染的像素點(diǎn)。目前這一類的算法已發(fā)展到很成熟的階段，其中具有代表性的有Windhorst等人的方法[1]。然而這些方法都依賴于對同一視場拍攝多幅圖像。

眾所周知，中值濾波方法[2]對去除隨機(jī)噪聲有良好的效果。然而該方法在對CCD圖像所有像素點(diǎn)進(jìn)行濾波的同時(shí)，也平滑了星像的邊緣，改變了星像中的有用信息?？梢娭兄禐V波方法不適用于高精度測量的CCD圖像處理。

因此，基于單幅圖像來消除宇宙射線更具有挑戰(zhàn)性。目前有基于單幅圖像消除宇宙射線方法[3]及線性濾波方法[4]。針對這一問題本文比較了3種較新的消除宇宙射線的算法：Laplacian邊緣檢測算法[5](簡稱Laplacian算法)，基于直方圖的快速算法[6](簡稱直方圖算法)以及萬能噪聲消除算法[7](簡稱萬能算法)。采用云南天文臺1m望遠(yuǎn)鏡拍攝的多幅圖像，研究星像和星系中人工及IRAF[8]軟件添加宇宙射線及其識別和剔除，通過編程實(shí)現(xiàn)3種算法來處理這些圖像，并分析比較每種算法的優(yōu)缺點(diǎn)。Laplacian算法被廣泛用于剔除宇宙射線，Van[5]指出該方法非常有效，但也認(rèn)為這種算法比較復(fù)雜，處理速度較慢。而直方圖方法由于算法簡單因此處理圖像的速度比較快。萬能方法聲稱可以處理任何噪聲，將該算法用于CCD圖像宇宙射線的消除，并與前述兩種算法進(jìn)行分析比較。為了進(jìn)一步探討如何準(zhǔn)確地替代宇宙射線像素點(diǎn)，將這3種算法分別用于處理實(shí)拍圖像，并采用兩類方法來替代宇宙射線像素點(diǎn)的灰度值：中值濾波方法、曲面擬合方法。

本文第1節(jié)將介紹上述3種算法的原理；第2節(jié)給出程序的設(shè)計(jì)與實(shí)現(xiàn)；第3節(jié)將通過對不同的CCD圖像進(jìn)行處理，給出實(shí)現(xiàn)過程和結(jié)果，比較3種算法的優(yōu)缺點(diǎn)，同時(shí)探討宇宙射線像素點(diǎn)的替代方法；第4節(jié)將對各種算法進(jìn)行綜合比較并得出結(jié)論。

1 算法介紹

1.1 Laplacian算法

Laplacian算法是通過構(gòu)造和設(shè)置兩個(gè)閾值來識別宇宙射線。它首先通過原始圖像I進(jìn)行子采樣放大(子采樣因子為fs)，放大后的圖像為I(2)。然后與Laplacian模板做卷積，如式(1)。▽2f為Laplacian算子，該算子是線性二次微分算子，用來獲取有大的灰度變化的圖像邊緣信息，當(dāng)然包括了宇宙射線影響的部分。之后將圖像L恢復(fù)到原始圖像分辨率得到圖像L+。

(1)

另一方面，對原始圖像I用來構(gòu)造剔除宇宙射線的噪聲模型，如式(2)。該模型僅由泊松噪聲和讀出噪聲組成，不包括宇宙射線噪聲。其中g(shù)是以光電子/ADU為單位的增益因子，σrn是以電子為單位的讀出噪聲，M5是5×5的中值濾波模板。卷積運(yùn)算M5*I得到的是有用信號的泊松噪聲平方的估計(jì)。

(2)

至此可以得到噪聲的比率S，如式(3)，其中fs是子采樣因子。這個(gè)比率S反映了圖像中每個(gè)像素點(diǎn)含有噪聲的比率，像素點(diǎn)噪聲越大則這個(gè)比值越大。然后設(shè)置閾值σlim，認(rèn)為S值大于σlim的那些像素點(diǎn)為候選的宇宙射線。

(3)

根據(jù)星像有對稱性而宇宙射線卻沒有的性質(zhì)，可以進(jìn)一步區(qū)別宇宙射線與星像，從而構(gòu)造一個(gè)精細(xì)結(jié)構(gòu)如式(4)。其中M3表示3×3的中值濾波模板，M7表示7×7的中值濾波模板。卷積運(yùn)算M3*I得到圖像的中、低頻信息，而(M3*I)*M7得到圖像的低頻信息，兩者相減得到圖像的中頻信息。容易推斷星像所在區(qū)域的F值會(huì)較大，而宇宙射線所在區(qū)域算出的該值會(huì)很小。

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

(4)

因此構(gòu)造第二個(gè)比率，如式(5)。通過設(shè)置閾值flim，認(rèn)為T大于flim的像素點(diǎn)為候選的宇宙射線，

(5)

最后同時(shí)設(shè)置兩個(gè)條件S>σlim和T>flim來探測宇宙射線。

該方法僅僅考慮用探測點(diǎn)附近的區(qū)域內(nèi)的中值濾波來替代宇宙射線點(diǎn)的原始灰度，是否有更好的替代算法沒有深入討論。算法在實(shí)現(xiàn)過程中可以迭代執(zhí)行，直到不能探測到新的宇宙射線才終止。

1.2 直方圖算法

直方圖方法是一種簡單且能直接用于識別宇宙射線的算法。該算法首先將圖像分成若干個(gè)子圖，然后分析這些子圖的直方圖(絕大多數(shù)直方圖分布是非常緊湊的，宇宙射線則會(huì)以孤立的單個(gè)點(diǎn)分布在這些緊湊范圍之外)，最后根據(jù)直方圖分布的特性設(shè)置閾值從而探測宇宙射線。對于閾值的設(shè)定范圍，Pych[6]建議為th*σ，th一般約等于3.0。每個(gè)子圖的σ計(jì)算公式，如式(6)，其中ci為子圖中的像素灰度值,n為子圖中的像素個(gè)數(shù)。

(6)

宇宙射線有時(shí)由多個(gè)像素點(diǎn)組成，該算法利用增長半徑的方法進(jìn)一步探測宇宙射線。它是以探測到的宇宙射線點(diǎn)為中心，在一定的半徑范圍內(nèi)繼續(xù)探測是否有其它的宇宙射線點(diǎn)。對探測到的宇宙射線點(diǎn)，通過設(shè)定一個(gè)范圍，然后在這個(gè)范圍內(nèi)求出未受宇宙射線影響的所有像素點(diǎn)的平均值，再用該平均值去替代宇宙射線點(diǎn)。

該算法在編程實(shí)現(xiàn)過程中需要輸入若干個(gè)參數(shù)，這些參數(shù)可以根據(jù)圖像的不同做調(diào)整，最終找到合適的值：

(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)算法。該算法的思想是通過排序統(tǒng)計(jì)像素點(diǎn)與周圍像素點(diǎn)之間的絕對值差，從而設(shè)置合適的ROAD閾值來探測識別噪聲點(diǎn)。為了更加清楚地理解ROAD算法，引用萬能算法[7]中的圖來具體分析，如圖1。

圖1 ROAD值的計(jì)算Fig.1 Calculation of the values of ROAD

2 程序的設(shè)計(jì)與實(shí)現(xiàn)

采用Visual C++ 6.0在Windows環(huán)境中編程[9]實(shí)現(xiàn)以上所述的3種算法。其中Laplacian算法和直方圖算法具體實(shí)現(xiàn)過程見圖2。其中圖2(a)和圖2(b)分別表示這兩種算法實(shí)現(xiàn)的流程圖。萬能算法由于算法本身比較簡單文中省去其流程圖。

3 實(shí)驗(yàn)與結(jié)果分析

以云南天文臺1m望遠(yuǎn)鏡拍攝的圖像進(jìn)行資料分析，實(shí)驗(yàn)分別采用Laplacian算法、直方圖算法和萬能算法來消除宇宙射線，重點(diǎn)研究星像和星系上人工及利用IRAF軟件添加宇宙射線的圖像。

根據(jù)上文提到的宇宙射線的性質(zhì)，在圖像的星像和星系上人工添加宇宙射線點(diǎn)，該宇宙射線點(diǎn)的灰度值要大于m+3N，其中m為該像素點(diǎn)的原始灰度值，N為公式(2)中的噪聲。

3.1 Laplacian算法處理星像和星系上的宇宙射線

Laplacian算法分別用于處理星像和星系上的宇宙射線，如圖3、4。圖3(a)表示在原始星像上坐標(biāo)為(1005,143)的位置，人工添加灰度值為1426的宇宙射線(并用黑色點(diǎn)表示以示區(qū)別)。星像周圍的亮點(diǎn)為原圖中帶有的宇宙射線點(diǎn)。通過設(shè)置兩個(gè)閾值分別為σlim=0.7、flim=2且經(jīng)過2次迭代運(yùn)算后，可以準(zhǔn)確識別到星像上和星像周圍的宇宙射線，如圖3(b)，并沒有產(chǎn)生假的探測。圖3(c)是用3×3中值濾波的方法剔除噪聲后的圖像，顯然，剔除宇宙射線噪聲后的視覺效果良好。

(a)Laplacian算法實(shí)現(xiàn)的流程圖 (b)直方圖算法實(shí)現(xiàn)的流程圖圖2 兩種算法實(shí)現(xiàn)的流程圖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)為在原始圖像的星系上添加宇宙射線點(diǎn)，該點(diǎn)的坐標(biāo)為(1031，1185)，灰度值為1424。圖4(b)顯示了經(jīng)過該算法處理后，探測到星系上的噪聲點(diǎn)。圖4(c)為消除宇宙射線后的圖像，圖像的視覺效果也很好。針對這幅圖像該算法設(shè)置的兩個(gè)閾值分別為σ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)清楚看出星像周圍留有未消除的宇宙射線。采用的參數(shù)分別為：子圖長寬都為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相比較，這兩種方法都能精確探測到星系上的宇宙射線。這是由于星系上像素點(diǎn)的灰度值變化比星像上灰度值變化要緩慢，因此星系上的宇宙射線比較容易探測到。該圖像采用的參數(shù)分別為：子圖長寬為8，th取值為3.2，增長半徑為1，消除宇宙射線的最大、最小半徑分別為2、1。

通過實(shí)驗(yàn)比較得出，Laplacian算法比直方圖算法更能準(zhǔn)確地識別和探測宇宙射線，且探測能力要明顯優(yōu)于直方圖算法。但從處理速度方面考慮，直方圖算法要快很多。

3.3 萬能算法處理星像和星系上的宇宙射線

實(shí)驗(yàn)發(fā)現(xiàn)對于在星像上添加宇宙射線的天文圖像，由于宇宙射線的ROAD值可以很大也可以很小，而星像上未受宇宙射線影響的像素點(diǎn)的ROAD值同樣具有此特點(diǎn)。如圖7，兩者ROAD值非常相似。因此萬能算法僅僅通過設(shè)定ROAD閾值，是較難區(qū)分星像和宇宙射線的。

萬能算法用于處理星系上的宇宙射線，如圖8，該方法通過設(shè)置ROAD閾值，可以精確剔除星系上的宇宙射線。

圖7 萬能算法用于星像上有宇宙射線的情況 ：B表示未受宇宙射線影響的背景區(qū)域，X表示未受宇宙射線影響的星像區(qū)域，A1，A2，A3，A4，A5表示背景上的宇宙射線，Y表示星像上的宇宙射線。圖中的表給出了這些區(qū)域上像素點(diǎn)的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望遠(yuǎn)鏡所拍攝。如圖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方法識別和剔除宇宙射線的效果較好。

為了進(jìn)一步研究宇宙射線對星像位置和光度的影響，重點(diǎn)討論利用實(shí)拍圖像中受宇宙射線影響的星像作為研究對象，對其位置和光度進(jìn)行分析。

(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 剔除實(shí)拍圖像中星像上的宇宙射線

將上述3種方法分別用于處理實(shí)際拍攝中星像受到宇宙射線影響的圖像。該圖像是2009年2月15日晚由云南天文臺1m望遠(yuǎn)鏡所拍攝。如圖13、14、15，分別為Laplacian算法、直方圖算法、萬能消除算法處理該圖，圖13(a)中的星像由于受到宇宙射線的影響，在測量該星像位置時(shí)，測量結(jié)果有嚴(yán)重的偏差。圖13(b)、14(b)、15(b)分別為3種算法識別的宇宙射線，其中區(qū)域A、B為星像上的宇宙射線，區(qū)域C為背景上的宇宙射線。圖13(c)、14(c)、15(c)為用下述替代方法替代宇宙射線后的圖像。

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

(a) (b) (c)

(a) (b) (c)圖15 萬能算法剔除實(shí)拍圖像中星像上的宇宙射線Fig.15 Rejection of real cosmic rays in an actual image with the universal noise removal algorithm

在消除宇宙射線噪聲之后，重新測量該星像的位置時(shí)，發(fā)現(xiàn)測量結(jié)果有明顯的改進(jìn)。改進(jìn)的好壞取決于如何替代宇宙射線像素點(diǎn)。采用的3種替代宇宙射線的算法分別為：3×3像素范圍的中值濾波法、5×5像素范圍的曲面擬合法和7×7像素范圍的曲面擬合法。曲面擬合公式如式(7)所示，其中x、y為像素點(diǎn)的坐標(biāo)位置，參數(shù)a、b、c、d、e、f為待求的擬合參數(shù)。

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

(7)

在連續(xù)曝光的一系列圖像，給定的未受宇宙射線影響的兩顆星之間的距離以像素為單位與光度差隨時(shí)間的變化過程是緩變的。如果其中有一顆星受到噪聲的影響，該緩變過程將會(huì)有突變。其中星像位置的測量采用高斯擬合[10]的方法，星像光度的測量采用孔徑測光的方法。

圖16 星像位置隨時(shí)間(北京時(shí))的變化Fig.16 Different positions of a star at different moments(in the Beijing time)

下面討論宇宙射線對星像位置及光度的影響，圖16為圖13(a)、13(c)、14(c)、15(c)中星像與其相鄰的一顆未受宇宙射線影響的亮星(圖中未給出)之間的距離隨時(shí)間變化的情況。可以看出，3種剔除宇宙射線的方法在20h35m42s(北京時(shí))這個(gè)時(shí)刻符合較好，且比受到宇宙射線影響的圖像的位置有較大改善(約0.16pixel)。圖17為圖13(a)、13(c)、14(c)、15(c)中星像與其相鄰的一顆未受宇宙射線影響的亮星(圖中未給出)之間的光度差隨時(shí)間變化的情況?？梢钥闯?，3種剔除宇宙射線的方法在20h35m42s(北京時(shí))這個(gè)時(shí)刻符合較好，且比受到宇宙射線影響的圖像的光度有較大提高(約0.04mag)，基本上達(dá)到預(yù)期目標(biāo)。

圖17 星像光度隨時(shí)間(北京時(shí))的變化Fig.17 Different measured magnitudes of a star at different moments(in the Beijing time)

4 結(jié) 論

對于星系上受宇宙射線影響的CCD圖像，Laplacian算法、直方圖算法、萬能算法都能很好地用于識別宇宙射線，對于星像上受宇宙射線影響的CCD圖像，Laplacian算法較其它兩種算法更能精準(zhǔn)地識別噪聲。從算法處理速度考慮，直方圖方法和萬能方法由于算法簡單因此處理速度更快。通過實(shí)驗(yàn)發(fā)現(xiàn)，宇宙射線的替代算法是值得深入研究的。

致謝:感謝暨南大學(xué)計(jì)算機(jī)科學(xué)系張慶豐老師、孟小華老師、李展老師為本文提出建設(shè)性建議。

[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, 著，阮秋琦，阮宇智 等譯.數(shù)字圖像處理[M]. 北京:電子工業(yè)出版社，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++圖像處理程序設(shè)計(jì)[M].北京:清華大學(xué)出版社, 2007.

[10]李展, 彭青玉, 韓國強(qiáng). CCD圖像數(shù)字定心算法的比較[J]. 天文學(xué)報(bào),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.