劉 奇，王 競，黃茂海，魏建彥

（1.中國科學院 國家天文臺 空間天文與技術重點實驗室，北京100101；2.中國科學院大學，北京100049；3.廣西大學 物理科學與技術學院 廣西相對論天體物理重點實驗室，廣西 南寧530004）

1 引 言

天文觀測圖像中的信息一般包含觀測源、背景、各種噪聲及宇宙射線等。其中，宇宙射線［1］是來自宇宙空間的各種高能粒子形成的射流，絕大多數是質子，還包含α粒子和少量其他原子核。在天文觀測圖像中宇宙線的灰度值通常高于周圍像素點的灰度值，并且分布規(guī)律具有隨機性。絕大多數的宇宙線在電荷耦合器件（Charge Cou?pled Device，CCD）圖像中呈現“點”狀或“線”狀分布，面積大的宇宙線也會和目標源形態(tài)輪廓相似。由于絕大多數天文圖像是地基觀測，受到地球大氣的影響，打在CCD靶面上的宇宙線是穿過大氣層，與大氣中的粒子發(fā)生相互作用后形成的次級粒子，常常被作為一種噪聲信號來處理。而在空間天文圖像中，由于不受地球大氣層的影響，CCD圖像中的宇宙線存在亮度高、同觀測目標源輪廓相似等特點，當圖像中出現非常多的宇宙線事件時對目標源的檢測會造成很大影響。

無論CCD圖像中的宇宙線呈現出哪種形態(tài)特征，后續(xù)處理分析的前提都是將它從圖像中提取出來。對于天文CCD圖像識別提取宇宙線，不同專家學者也提出了不同算法。基于多幀圖像，運用圖像對齊相減的方法確定宇宙線的位置。宇宙線打到CCD探測器上具有隨機性的特點，通過對相同的視場進行連續(xù)多幀拍攝，然后依據序列圖像信息確定宇宙線位置。Windhorst等［2］就是運用這種方法從哈勃空間望遠鏡寬視場相機圖像中定位、剔除了宇宙線事件。此類算法到目前為止已經相當成熟，但其缺點是需要對同一視場進行多次拍攝，降低了儀器的觀測時間和使用效率，并且無法通過單一圖像來定位宇宙線。

基于單幀圖像識別提取識別宇宙線的方法更具挑戰(zhàn)性。從數字圖像處理的角度出發(fā)，中值濾波可消除隨機噪聲，如果把宇宙線作為天文圖像中的隨機噪聲，對單幀圖像使用中值濾波可以很好地探測到宇宙線。但是由于天文CCD圖像作為高精度數字圖像，在中值濾波的同時會平滑掉恒星目標的邊緣，改變了其有用信息。因此，傳統(tǒng)的中值濾波并不適用于天文觀測圖像。Van等［3］提出了一種運用拉普拉斯邊緣檢測算子來檢測宇宙線的方法。此方法通過構建和設置兩個閾值來識別宇宙線。一個閾值量是反應每個像素點含有噪聲的噪聲比率，像素點噪聲越大，這個噪聲比率越高，大于一定閾值的噪聲比率認為是宇宙線候選體；第二個比率是原始圖像與拉普拉斯算子卷積后得到圖像與圖像的中頻信息的商，圖像的中頻信息是依據星像有對稱性而宇宙線沒有的特點，通過與不同中值濾波模板進行卷積相減得到的。通過設置這兩個閾值，算法不斷迭代，直到不能檢測到宇宙線為止。該方法在宇宙線檢測上非常有效，但算法復雜，處理速度較慢，而且對于形態(tài)輪廓與星像相似的宇宙線檢測準確率較低。相比于拉普拉斯邊緣檢測算法，Pych等［4］提出了一種簡單、直接的直方圖算法用來識別宇宙線。該算法直接將圖像分為若干子圖像，并分析其直方圖。由于絕大多數子圖像的直方圖分布比較緊湊，一旦有宇宙線事件，直方圖中在緊湊結構之外就會出現孤立的單點，從而被分辨出來。針對此方法，不同學者提出了不同的閾值設定定義，并且在算法實現過程中，其參數需要根據原始圖像的不同做出調整。相比拉普拉斯邊緣檢測算法，算法運行效率雖高，但其探測能力明顯低于前者。

通過統(tǒng)計分析月基光學望遠鏡（Lunar-based Ultraviolet Telescope，LUT）［5］數據，發(fā)現相比于地基望遠鏡，宇宙線打在LUT CCD靶面上的發(fā)生概率顯著提高，除地基望遠鏡圖像中出現的“點”狀或“線”狀特征外，還存在很多亮度高、影響面積大、形態(tài)輪廓與星像類似的宇宙線。對于此類形態(tài)輪廓與星像類似的宇宙線，傳統(tǒng)算法無法精確檢測處理?；诖藛栴}，本文提出了一種運用天文位置定標來識別提取宇宙線的方法。此方法與上述文獻中方法不同的是從天文圖像處理角度出發(fā)，不會引入其他噪聲影響，對于形態(tài)輪廓與星像相似的宇宙線識別探測效率更高。識別提取到宇宙線事件后，對其電子沉積數和入射角進行了分析，并將其形態(tài)進行提取、分類，制作出“CCD探測器宇宙線實體樣本庫”，為空間天文望遠鏡圖像的仿真提供經驗支持。

2 月基光學望遠鏡

LUT的結構與安裝示意圖如圖1所示。望遠鏡的主鏡直徑為150 mm，焦比為F/3.75，采用了里奇-克列基昂（R-C）系統(tǒng)設計。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觀測圖像統(tǒng)計發(fā)現，在大多數情況下約0.5 h內拍攝的圖像中雜散光的強度和大尺度結構基本不變。

本文根據LUT觀測的時間對圖片進行了分組，將連續(xù)觀測1 900 s以內的序列圖像分為一組。對某一幀圖像，將同組中除自身以外的所有圖像對齊后進行中值合并，宇宙線被剔除；同時，由于LUT步進指向跟蹤模式下，恒星在連續(xù)曝光的圖像上的位置總是有輕微的幾個像元的變化，所以通過中值合并，所有的天體目標也都被剔除了。因此，合并所得圖像僅保留了背景和CCD的“熱點”和“壞點”，得到的背景作為本幀圖像的背景模板。

最后，對每一幀圖像，減去圖像對應的背景模板，背景被扣除，同時CCD的“熱點”和“壞點”也被剔除，得到的圖像僅保留了恒星星像和宇宙線。

3.2 平場修正

平場修正是為了解決LUT CCD相機像面響應的非均勻性。LUT圖像的平場分為小尺度平場和大尺度平場兩部分［7］。小尺度平場反映了CCD像元間量子效率的不均勻性或像元有效面積的變化；大尺度不均勻性主要來自望遠鏡的光學系統(tǒng)，還有小部分來自于CCD鍍膜導致的大尺度效率變化。

對于小尺度平場，通過拍攝LUT內部平場燈圖像，即可測量像面上不同像元間響應的非均勻性。對于大尺度平場，是通過對標準星進行“抖動”觀測完成的。最終平場是大尺度平場與LED平面場的乘積。

3.3 天文位置定標

本文天文位置定標所用的星表為第谷星表（Tycho_2.0）［8］，主要是因為第谷星表的天文坐標精確高并且與LUT的觀測深度相近。但第谷星表采用雙色測光，王競等［9］利用恒星大氣模型計算出了LUT相應的近紫外AB星等。

天文位置定標的關鍵是將CCD圖像中的恒星星像位置和第谷星表中的天文坐標進行交叉認證。通過LUT的指向文件記錄獲得CCD圖像中心位置的J2000天文坐標，精度通常優(yōu)于1個角分。圍繞此中心位置提取交叉所用的第谷星表。

對于每幀圖像，這里使用SExtractor（Source-Extractor）軟件［10］提取出5～10顆明亮的恒星目標，并對它進行快速測光，獲得其近紫外AB星等。將恒星目標與對應的第谷星表進行位置匹配和亮度匹配。位置匹配交叉半徑設置為10角秒（約為LUT圖像的2個像元），亮度匹配閾值控制在2個星等以內（主要考慮了變星的影響）。

利用交叉匹配的結果生成CCD圖像上X-Y位置與天文坐標的位置對應關系，即世界坐標系統(tǒng)（World Coordinate System，WCS）的位置關系。通過WCS實現了對CCD圖像上每一個像元的對應天文位置。天文位置的定標精度優(yōu)于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等，由于其亮星飽和問題，星表在亮星部分不完備。所以將兩個星表綜合利用，保證了交叉星表的完備性。

為選取合適的交叉半徑，本文統(tǒng)計了6幀圖像中的候選體與兩個星表的最近距離。首先選取10角秒為半徑與兩個星表進行位置交叉，結果表明，相匹配的源中97.21%的距離都小于4角秒，小于LUT一個像元的尺寸4.76角秒。最終這里采用4角秒作為交叉半徑。

圖3為12幀圖像中宇宙線候選體與兩個星表交叉獲得的星像數量統(tǒng)計圖（彩圖見期刊電子版）。圖中，藍色表示與Tycho_2.0星表交叉獲得的恒星目標數，綠色表示與USNO_B1.0星表交叉獲得的恒星目標數，黃色為對兩個交叉結果取“并集”獲得的恒星目標數。很明顯，“并集”獲得到的恒星目標完備性更高。

圖3 星像數量統(tǒng)計圖Fig.3 Statistical chart of star quantity

圖4 是一幀LUT圖像宇宙線認證識別結果（彩圖見期刊電子版）。圖中，紅色圓圈標記宇宙線事件，藍色圓圈標記恒星目標。很明顯，通過以上處理流程能夠將較亮的宇宙線與恒星目標很好地識別區(qū)分開來。

圖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的統(tǒng)計結果如圖5所示。很明顯，宇宙線事件的FWHM明顯小于恒星目標，但仍有一部分相重疊。其中，90%的宇宙線樣本的FWHM小于1.6 pixel；95%的恒星目標的FWHM為1.5～2.2 pixel，而在這個區(qū)間內的宇宙線事件數量占宇宙線事件總數的8.6%。這8.6%的宇宙線事件在輪廓上與恒星目標完全類似，但本算法還是可以將它們識別出來，并和恒星目標區(qū)分開來。

圖5 半高全寬統(tǒng)計圖Fig.5 Histogram of FWHM

3.6.2 宇宙線識別完備性檢驗

對于天文CCD圖像，拉普拉斯算法［3］比直方圖算法［4］和萬能消除算法［6］更能準確地識別宇宙射線，且探測能力要明顯優(yōu)于兩者［12］。本文將基于天文位置定標的宇宙線識別算法和拉普拉斯算法的宇宙線識別結果進行了統(tǒng)計對比。

對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例樣本進行檢查，通過與恒星目標進行位置認證，發(fā)現其中的251例為恒星目標，即本算法沒有識別到的樣本中48.64%為恒星目標，是錯誤識別。

4 處理結果與分析

4.1 宇宙線事件電子沉積數分析

本文將單例宇宙線事件影響CCD像元上ADU（Analog-Digital Unit）值相加，根據CCD增益值GIAM=1.59（Electrons/ADU）計算出每例宇宙線事件的電子沉積數，取對數后做出了電子沉積數統(tǒng)計圖，并與HST探測器圖像中宇宙線事件的電子沉積數［13］進行對比，結果如圖6所示。

圖6 宇宙線樣本電子沉積統(tǒng)計圖Fig.6 Electron deposition by cosmic rays

通過與HST巡天相機（Advanced Camera for Surveys，ACS）上WFC與高分辨率探測器（High Resolution Channel，HRC）圖像中宇宙線電子沉積數分布進行對比，發(fā)現分布形態(tài)基本類似，但LUT圖像中宇宙線事件電子沉積數的峰值要略高于HST的電子沉積數峰值。最可能的原因是，宇宙射線打到CCD靶面上，其電子沉積數會隨著CCD探測器敏感層的厚度不同而發(fā)生變化［14］。

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

宇宙線打到CCD靶面時，留在CCD上的形態(tài)主要是由宇宙線粒子的能量和入射角度所決定的。為了量化定義宇宙線粒子對CCD帶來的影響，依據粒子射入CCD成像元件敏感層的瞬態(tài)效應，原子量大于1的粒子在器件中的軌跡可以認為是一條直線，粒子和CCD作用后產生的電荷信號會呈現集中的光點輪廓［15］。宇宙線與CCD作用的簡單形態(tài)模型如圖7所示。

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

宇宙線事件在CCD上留下的形態(tài)輪廓中，定義通過重心位置的最長線段為長軸，由重心指向長軸與宇宙線截取輪廓最遠交匯點的方向，就是宇宙線事件在CCD平面上的入射方向，如圖7所示的θ角。

圖8為一幀圖像上識別到的所有宇宙線事件入射角度示意圖。圖8中，箭頭方向為宇宙線事件在CCD平面上的入射方向。可以看出，一幀圖像中，CCD靶面上能夠接收到來自不同方向的宇宙線事件。

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

本文隨機選取12幀圖像，對識別到的1 280例宇宙線事件計算其入射角，并進行了統(tǒng)計分析。統(tǒng)計時平均分成了24個方向，每個方向包含15°區(qū)間，結果如圖9所示。在入射角為45°～60°，285°～300°上分別有167，211例事件，其他方向上平均為41例，這兩個方向上的宇宙線事件發(fā)生頻次明顯多于其他方向。

圖9 CCD靶面上宇宙線入射角統(tǒng)計圖Fig.9 Statistical chart of incidence angle of cosmic rays on CCD

在這些數據采集期間，宇宙線事件可認為是各向同性的，探測到的宇宙線事件集中在兩個方向上，應該歸咎于這兩個方向上的等效鋁厚度較薄。原則上依據輻射屏蔽模型［16］，如果已知嫦娥三號著陸器的物質結構模型，可以仿真計算各個方向上的等效鋁厚度。因暫時無法得到嫦娥三號著陸器的物質結構模型，本文在此不做定量分析。

4.3 建立“宇宙線實體樣本庫”

利用天文位置定標算法對230幀LUT圖像中的宇宙線事件進行識別提取，并依據其形態(tài)特征進行分類，制作出“宇宙線實體樣本庫”。此“宇宙線實體樣本庫”可作為輸入來源，為后續(xù)空間天文圖像中宇宙線的仿真研究提供支持。

對于230幀LUT圖像，首先依據每幀CCD圖像的背景噪聲值，將信號值小于5倍背景噪聲的像元賦值為0；然后對于識別出的恒星目標，以恒星目標的重心為中心，將恒星目標大于5倍背景噪聲的連通區(qū)域賦值為0；最后將連通像元數小于4倍、大于5倍背景噪聲的連通區(qū)域也賦值為0。經過此過程，CCD圖像中就只保留了宇宙線事件的信號值，其他像元值均為0。

以宇宙線的重心位置為中心，在CCD圖像中搜索非零值的連通區(qū)域，直至此連通區(qū)域周圍均為0值，定義此連通區(qū)域為一例宇宙線事件。從230幀圖像中共識別出29 494例宇宙線事件。

對于每一例宇宙線事件，本文截取出包含這一宇宙線事件完整信號的圖片，此圖片在宇宙線連通區(qū)域四邊各多出一個像元，因此圖片為長方形或方形。截取圖片樣例如圖10所示。

圖10 宇宙線事件樣例及對應的三維柱狀圖Fig.10 Cosmic ray samples and their three dimensional bar chart

對于29 494幅截取圖片，依據宇宙線事件信號值衰減的類型，將宇宙線事件分為兩類：衰減類型單一的單峰宇宙線事件，衰減類型非單一的多峰宇宙線事件。單峰宇宙線事件如圖10（a）～10（c）所示；多峰宇宙線事件如圖10（d）所示。統(tǒng)計得出，單峰宇宙線事件有26 459例，占事件總數的89.71%；多峰宇宙線事件僅有3 035例，占事件總數的10.29%。分類結果如表3所示。

表3 宇宙線事件的分類結果Tab.3 Classification results of cosmic rays

最后，將分類完成的宇宙線事件截取圖片編號后存儲，形成“宇宙線實體樣本庫”。

5 結 論

在空間天文CCD圖像中，宇宙線識別遇到的最大挑戰(zhàn)是將形態(tài)與星像類似的宇宙線事件與恒星目標區(qū)分開來。為了盡可能準確地識別宇宙線事件，本文基于嫦娥三號LUT獲得的CCD圖像，設計和實現了一套基于天文位置定標的宇宙線識別算法，該算法與常用的拉普拉斯算法相比，檢測出事件總數多出11.14%。對利用該算法識別出的宇宙線事件，進一步分析了電子沉積數及入射角的分布情況，發(fā)現宇宙線電子沉積分布與HST基本類似，入射角度集中在兩個特定方向上，初步判斷是由嫦娥三號著陸器內部屏蔽效應造成的。最后，制作出一個包含29 494例宇宙線事件截取圖片的空間天文CCD圖像“宇宙線實體樣本庫”，可為空間天文望遠鏡圖像仿真系統(tǒng)中宇宙線事件的仿真提供支持，具 有重要的應用價值。