馮海霞，陳建軍，鄧建榕，趙永恒

(1. 中國科學院國家天文臺，北京 100101；2. 中國科學院大學，北京 100049；3. 中國科學院光學天文重點實驗室 (國家天文臺)，北京 100101)

宇宙線是來自外太空的高能粒子，1912年，奧地利地理學家赫斯在電離室測量電流時第一次發(fā)現(xiàn)了宇宙線。宇宙線中質子和氦核占98%，重核1%，電子1%，微量其他粒子[1]。當初級宇宙線進入地球大氣層時，質子和氦核等高能粒子與大氣中的原子核相互作用產(chǎn)生級聯(lián)的廣延大氣簇射，生成大量的次級粒子。在使用電荷耦合器件觀測目標源時，這些粒子在CCD上留下不同的徑跡，其中包括電子、α粒子、X射線以及μ子等?？拷髿獾撞康挠钪婢€主要是簇射產(chǎn)生的次級μ子。在CCD圖像的宇宙線中，占據(jù)主要成分的是μ子。μ子和電子等帶電粒子在CCD上留下明顯的特征徑跡。電子大部分是由環(huán)境中放射性同位素產(chǎn)生的，并非來源于宇宙線[2]。α粒子和X射線等其他粒子在CCD圖像中呈點狀，不易與CCD芯片本身產(chǎn)生的電子噪聲分開。因此，研究CCD圖像中的宇宙線，需要對宇宙線中的μ子進行甄選。

當前探測宇宙線的裝置很多，例如西藏羊八井國際宇宙線觀測站、印度烏蒂的GRAPES-3、阿根廷的皮埃爾·俄歇觀測站(Pierre Auger Observatory)以及四川的高海拔宇宙線觀測站等。這些裝置采用一系列探測器陣列和望遠鏡來探測宇宙線粒子。本文利用大天區(qū)面積多目標光纖光譜天文望遠鏡(the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST，又叫郭守敬望遠鏡)探測μ子。文[3]系統(tǒng)介紹了郭守敬望遠鏡的結構、原理及應用。郭守敬望遠鏡配置有16個光譜儀和32臺4 k × 4 k的CCD相機。CCD芯片大小為49.2 mm × 49.2 mm，像素大小為12 μm × 12 μm，工作溫度在-100 ℃左右。從2011年先導巡天至今，已經(jīng)積累了一百多萬張原始圖像，每幅圖像都有很多宇宙線事件，可以用來研究宇宙線的性質和分布。圖像中的宇宙線個數(shù)與圖像的曝光時間有關，將曝光時間歸一化到分鐘。在單位時間內，圖像中的宇宙線事件個數(shù)為5～60。

CCD的工作原理是利用光電效應將光子轉化為電子，電容器吸收釋放的電子，然后將收集的電子轉化為電壓，最后以數(shù)字信號的形式輸出。μ子有帶負電和帶正電兩種，電荷數(shù)為1。它的質量約為106 MeV，比電子重207倍。當它穿過CCD時，可以電離徑跡上的硅原子，受到庫倫散射和軔致輻射的影響比電子小，在CCD上留下可識別的直線徑跡。

文[2]介紹CCD圖像中的宇宙線事件關于質心呈現(xiàn)不同的分布，結合宇宙線事件的形態(tài)特征甄選μ子。本文基于郭守敬望遠鏡的觀測數(shù)據(jù)，對圖像進行分析處理，根據(jù)μ子在圖像中的形態(tài)特征對其進行甄選。首先從圖像中提取宇宙線像素，采用拉普拉斯邊緣檢測法從原始圖像中提取宇宙線候選像素列表。在提取的宇宙線候選像素列表中包含一些壞像素和噪點，需要對這些壞像素和噪點進行處理。數(shù)據(jù)處理過程中，為了得到獨立完整的宇宙線事件，采用凝聚層次聚類算法將宇宙線像素聚類成宇宙線事件。提取宇宙線事件的特征，結合μ子在圖像中的形態(tài)特征甄選出宇宙線μ子。最后對甄選技術進行評估。

1 拉普拉斯邊緣檢測法提取宇宙線候選像素列表

在郭守敬望遠鏡的觀測數(shù)據(jù)中，單次曝光的CCD圖像中含有上百條流量值很強的光纖光譜。這些圖像中，不僅包含所要觀測的目標信息，同時還包括一些噪聲。在對圖像中的μ子進行甄選時，首先從圖像中獲取宇宙線候選像素列表。檢測圖像宇宙線的方法有很多，例如經(jīng)典的中值濾波法[4]、拉普拉斯邊緣檢測法[5]、基于直方圖的快速算法[6]以及萬能噪聲消除算法[7]等。本文采用改進的拉普拉斯邊緣檢測法[8]提取圖像中的宇宙線候選像素列表。相比于其他方法，該方法在圖像背景流量值很強的情況下，可以有效地從單次曝光的多光纖光譜中檢測宇宙線。它是在文[5]算法的基礎上進行改進,使得當宇宙線落在光纖光譜輪廓內時可以被識別出來。

改進的拉普拉斯邊緣檢測算法的流程：

(1)對原圖像進行預處理，然后與尺寸大小為3 × 3的拉普拉斯模板進行卷積運算得到圖像I；結合泊松噪聲和讀出噪聲，在原圖像上構造剔除宇宙線的噪聲模型N；再根據(jù)I和N的比率關系，生成一個原始的宇宙線候選像素列表；

(2)在原圖像的基礎上去除步驟(1)產(chǎn)生的宇宙線候選像素列表，生成一個新的圖像。對圖像中的每一條光纖進行處理，將距離光纖軌跡中心8個像素以內的像素點組成一個模塊。結合光譜在空間方向和色散方向的分布情況，擬合這些模塊的圖像輪廓(這些擬合的圖像輪廓不包括宇宙線)，得到擬合后的圖像I1；

(3)根據(jù)泊松噪聲和讀出噪聲，對擬合后的圖像I1構造剔除宇宙線的噪聲模型N1；最后將原圖像和擬合后圖像的殘差與噪聲模型N1進行比較，生成宇宙線候選像素列表。

采用拉普拉斯邊緣檢測法提取宇宙線候選像素列表的結果見圖1。圖1(a)是從CCD圖像中截取的一部分原始圖像；圖1(b)對應于圖1(a)采用拉普拉斯邊緣檢測法提取后的宇宙線像素，圖中的背景為黑色，對應的流量值為0。從圖1可以直觀地看出，該算法可以有效地從圖像中提取宇宙線像素。

2 壞像素和噪點的處理

在使用拉普拉斯邊緣檢測法提取的宇宙線候選像素列表中，包括宇宙線、CCD上的壞像素以及算法導致的噪點。CCD上的壞像素是本身的壞點，基本在同一批拍攝的圖像中的同一位置出現(xiàn)?？梢岳猛荒繕说亩啻纹毓鈭D像標識這些壞像素，然后將它們去除。在使用拉普拉斯邊緣檢測法提取宇宙線候選像素列表時，如果圖像中的光纖軌跡中心發(fā)生偏移，會在擬合圖像輪廓時產(chǎn)生較大的誤差，最終在生成宇宙線列表的基礎上產(chǎn)生擬合的偏差，這些偏差就是算法導致的噪點。這些噪點在圖像中呈現(xiàn)幾列虛線，虛線處的像素點個數(shù)遠超過宇宙線在圖像中每一列可能的像素點個數(shù)(虛線處的像素點個數(shù)為150～500，而圖像中每一列的宇宙線像素點個數(shù)為8～50)。利用這一特點，在圖像中找出異常的列來標注算法導致的噪點，然后將它們從圖像中去除。

圖1 拉普拉斯邊緣檢測法提取宇宙線

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

3 凝聚層次聚類算法聚類宇宙線事件

從圖像中提取宇宙線像素后，采用凝聚層次聚類算法將宇宙線像素聚類成一個個獨立的宇宙線事件。本文采用改進的經(jīng)典凝聚層次聚類算法。

層次聚類試圖在不同層次對數(shù)據(jù)集進行劃分，從而形成樹形的聚類結構。數(shù)據(jù)集的劃分可采用自底向上的聚合策略，也可采用自頂向下的分拆策略[9]。凝聚層次聚類是層次聚類中的一種，它對數(shù)據(jù)集的劃分采用自底向上的聚合策略。主要思路是將數(shù)據(jù)集中的每一個樣本看成單個的聚類簇，然后計算每個聚類簇之間的距離，找出其中距離最近的兩個聚類簇將它們合并，重復上述距離計算和合并的過程，直至達到預設的聚類簇個數(shù)。

在經(jīng)典的凝聚層次聚類算法中，算法的關鍵是計算類間距的方法和確定聚類簇的個數(shù)。每一個聚類簇是一個樣本集合，可以采用計算集合的某種距離來計算類間距，常見的距離有最小距離、最大距離和平均距離。將宇宙線像素聚類成宇宙線事件，使得宇宙線事件之間相互獨立，采用最小距離計算類間距。最小距離計算公式：

(1)

其中，Ci，Cj為給定的聚類簇；{ma}a=1,2, ...,k為Ci的樣本點集；{mb}b=1,2, ...,l為Cj的樣本點集；dist為計算ma,mb樣本點的距離函數(shù)。

在處理圖像時，無法預知圖像中宇宙線事件的個數(shù)，所以無法預設聚類簇的個數(shù)。需要設定類間距dlim作為終止運算的條件。具體算法如下：

(1)假設數(shù)據(jù)集為D，數(shù)據(jù)集中的數(shù)據(jù)是圖像中宇宙線所對應的像素，將每個像素看作一個初始的聚類簇分別對應C1,C2, ...,Cn。

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

(2)

(2)計算所有聚類簇之間的類間距，如果兩個聚類簇各自只含有一個樣本點時，類間距是這兩個樣本點的歐氏距離；反之，類間距是兩個聚類簇之間的最小距離。如果兩個聚類簇的類間距不大于設定的閾值dlim，將這兩個聚類簇合并成一個聚類簇。采用歐氏距離和最小距離計算類間距：

(3)

其中，Ci,Cj為數(shù)據(jù)集中的任意兩個聚類簇；{ma}a=1,2, ...,k為Ci的樣本點集；{mb}b=1,2, ...,l為Cj的樣本點集；xa和ya對應ma在圖像上的位置信息；xb和yb對應mb在圖像上的位置信息。

(3)更新數(shù)據(jù)集。

(4)重復步驟(2)、(3)，直到數(shù)據(jù)集中所有聚類簇的類間距大于設定的閾值，結束運算。

當圖像中出現(xiàn)大量宇宙線事件時，需要逐個計算聚類簇之間的距離，計算量劇增，運行速率減慢。結合宇宙線事件是有計數(shù)的連續(xù)像素這一特點，通過判斷與聚類簇最外層元素距離為dlim的范圍內是否存在其他的聚類簇優(yōu)化算法，如果存在其他聚類簇，將其他聚類簇與該聚類簇合并；如果不存在，將數(shù)據(jù)集中未處理的聚類簇重復上述判斷過程，直到所有聚類簇的類間距大于dlim，結束運算。使用該方法將宇宙線像素聚類成宇宙線事件時，如果兩個宇宙線事件的最小距離大于dlim，說明兩個宇宙線事件相互獨立；如果最小距離不大于dlim，說明是同一個宇宙線事件。在提取宇宙線像素的過程中，當宇宙線流量值和背景值比較接近時，會出現(xiàn)提取的宇宙線像素在圖像中不連續(xù)。盡管這種情況很少出現(xiàn)，為減小宇宙線像素缺失帶來的誤差，結合宇宙線在圖像中呈現(xiàn)隨機稀疏分布的特性，設定dlim為3對提取后的宇宙線像素聚類，可以有效地聚類宇宙線事件。

4 特征提取

將圖像中的宇宙線像素聚類成宇宙線事件后，結合μ子在圖像中的形態(tài)特征，對宇宙線事件進行特征提取。在特征提取的過程中，主要使用主成分分析和皮爾森(Pearson)相關系數(shù)。通過主成分分析獲取宇宙線事件的特征值，通過皮爾森相關系數(shù)判斷宇宙線事件是否呈線性，以此甄選宇宙線事件中的μ子。

4.1 主成分分析

主成分分析也稱為Karhunen-Loeve擴展，是一種經(jīng)典的特征提取和數(shù)據(jù)表示技術[10]。文[11]對主成分分析進行了總的概括。主成分分析也可以用于對數(shù)據(jù)進行降維處理，目的是找一組最優(yōu)的標準正交向量基，使得數(shù)據(jù)投影到這組最優(yōu)標準正交向量基上的方差最大。宇宙線事件是二維數(shù)據(jù)，假設在一個宇宙線事件中含有n個像素，其樣本集M={m1,m2, ...,mn}；將所有樣本進行中心化：

(4)

定義樣本集的散度矩陣為St：

(5)

(6)

假設從高維空間投影到低維空間的投影矩陣為w，使投影后樣本點方差最大化的優(yōu)化目標為

Stw=λw.

(7)

對St做特征值分解得到特征值λ，將特征值λ進行降序排序，選取最大的特征值，最大的特征值所對應的特征向量就是主成分的解，從而找到最優(yōu)的標準正交向量基，得到宇宙線事件的特征向量。將宇宙線事件在特征向量方向上的投影作為宇宙線事件的特征。兩個特征向量對應于兩個特征值。其中一個特征值對應宇宙線的長軸，即在CCD上宇宙線直線徑跡所對應的軸；另一個特征值對應宇宙線的短軸，即宇宙線在圖像上的寬度。μ子在圖像上的徑跡是直線型，且徑跡長度無法確定，但寬度有限，是個小量。選取標記寬度的特征作為μ子的特征，并將其記為投影寬度。

4.2 皮爾森相關系數(shù)

相關系數(shù)是用來衡量變量之間線性相關程度的統(tǒng)計量。定義相關系數(shù)的方式有很多，最常用的是皮爾森相關系數(shù)，皮爾森相關系數(shù)用于描述兩個變量之間是否存在線性關系。文[12]對皮爾森相關系數(shù)進行了詳細介紹。假設兩個變量分別是x和y，它們的皮爾森相關系數(shù)用r表示，則有

(8)

其中，μx，μy分別為變量x，y的均值。r的取值范圍為[-1, 1]，|r|越接近1時，兩個變量之間的線性相關性越顯著。當r> 0時，兩個變量之間呈正相關；當r< 0時，兩個變量之間呈負相關；當r=0時，兩個變量之間彼此獨立。

宇宙線μ子在CCD上呈直線軌跡，則μ子在圖像中的橫、縱坐標存在線性關系，利用皮爾森相關系數(shù)標記μ子的這一特征。由于μ子在圖像中的位置有時呈現(xiàn)正相關，有時呈現(xiàn)負相關，為了顯示線性相關性質的顯著程度，在傳統(tǒng)的皮爾森相關系數(shù)的基礎上取絕對值。下文提到的r為皮爾森相關系數(shù)的絕對值。

5 實驗結果及分析

利用宇宙線μ子與其他粒子在圖像中的形態(tài)特征差異甄選μ子。在圖像中最常見的粒子是μ子和 “worms”(“worms” 是多級散射的電子，詳細介紹見文[2])。α粒子在圖像上極少出現(xiàn)。CCD圖像中還存在一些形狀奇特的宇宙線事件，因無法確定其形態(tài)特征，難以判斷屬于何種宇宙線粒子，本文未做討論。圖2是從郭守敬望遠鏡獲得的圖像，圖2(a)是宇宙線中的μ子，圖2(b)是宇宙線中的 “worm”， “worm” 形狀呈現(xiàn)一定的弧度。

圖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。宇宙線事件的位置索引信息可以反映宇宙線事件在圖像上呈現(xiàn)的形態(tài)。從樣本集中通過人眼識別的方式，選取1 000個μ子和1 000個 “worms” 作為數(shù)據(jù)樣本1。統(tǒng)計這些宇宙線事件的投影寬度和皮爾森相關系數(shù)，結果見表1。表1中宇宙線事件投影寬度的均值記為μwp，標準差記為σwp；皮爾森相關系數(shù)的均值記為μr，標準差記為σr。

表1 關于投影寬度和皮爾森相關系數(shù)，宇宙線μ子、“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

根據(jù)表1可以直觀地看出，μ子和 “worms” 的區(qū)別。在投影寬度中，μ子的特征值小，而 “worms” 的特征值大；在皮爾森相關系數(shù)中，由于μ子在CCD中呈直線徑跡，對應的橫縱坐標的線性相關程度更顯著，對應的皮爾森相關系數(shù)值更接近1，而 “worms” 由于在CCD中發(fā)生多次散射，徑跡有一定的弧度，對應的皮爾森相關系數(shù)值較μ子的更小。

通過人眼識別的方式從樣本集C1中隨機選取2 500個μ子和2 500個非μ子的其他宇宙線事件(包含 “worms”、X射線、可能的α粒子等其他粒子)作為數(shù)據(jù)樣本2。這5 000個宇宙線事件的投影寬度和皮爾森相關系數(shù)的分布結果如圖3(a)。在圖3(a)中，利用投影寬度和皮爾森相關系數(shù)兩個特征可以將絕大多數(shù)的μ子和其他宇宙線區(qū)分開。在圖中，皮爾森相關系數(shù)大于0.8和投影寬度小于1的范圍內，μ子和其他宇宙線事件不能很好地區(qū)分。觀察這部分與μ子特征值相近的宇宙線事件，發(fā)現(xiàn)這些宇宙線事件的像素個數(shù)很少。分析皮爾森相關系數(shù)和宇宙線事件像素個數(shù)的關聯(lián)，二者的分布結果見圖3(b)。圖3(a)中不易區(qū)分的μ子和其他宇宙線事件，在圖3(b)中可以通過宇宙線事件的像素個數(shù)有效地區(qū)分。

圖3 不同特征值的分布結果。(a) 投影寬度和皮爾森相關系數(shù)的分布結果；(b) 皮爾森相關系數(shù)和像素點個數(shù)的分布

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

對于像素個數(shù)少的宇宙線事件，只結合投影寬度和皮爾森相關系數(shù)甄選μ子會有一定的誤差。當甄選μ子時，需要結合宇宙線事件的像素個數(shù)。圖4顯示了數(shù)據(jù)樣本2中宇宙線事件像素個數(shù)的統(tǒng)計結果。橫坐標表示宇宙線事件的像素點個數(shù)，縱坐標表示像素點個數(shù)對應的宇宙線事件的個數(shù)。圖4(a)是整個數(shù)據(jù)集的統(tǒng)計結果。圖中宇宙線事件的像素點個數(shù)越多，對應的宇宙線事件個數(shù)越少?；趫D4(a)，選取像素點個數(shù)不大于30的宇宙線事件進行統(tǒng)計，結果見圖4(b)。圖中，宇宙線μ子的像素點個數(shù)大于8，其他宇宙線事件的像素點個數(shù)大于5。在非μ子的其他宇宙線中，將近一半的宇宙線的像素點個數(shù)小于8。在像素點個數(shù)大于8的情況下，μ子的個數(shù)比其他宇宙線的個數(shù)多。因此，將像素點個數(shù)設置為8，對宇宙線中μ子和其他宇宙線事件進行分類。

圖4 宇宙線事件的像素點個數(shù)的統(tǒng)計結果

根據(jù)宇宙線事件特征值的不同分布，設定不同的特征值參數(shù)，對μ子進行甄選。從樣本集C1中，隨機選取2 500個μ子和2 500個非μ子的宇宙線在CCD上的位置索引信息作為測試數(shù)據(jù)集C2，設定不同的特征值參數(shù)對宇宙線中的μ子進行甄選，最后的甄選結果如表2。表2顯示了甄選結果的準確性，表中第1列是特征參數(shù)的設定，其中，μmw為μ子投影寬度的均值，σmw為μ子投影寬度的標準差，μmr為μ子皮爾森相關系數(shù)的均值，σmr為μ子皮爾森相關系數(shù)的標準差；第2列是正確甄選宇宙線中μ子的個數(shù)；第3列是被錯分成μ子的宇宙線事件個數(shù)；第4列是甄選結果的正確性；最后一列是甄選結果的靈敏度，用于衡量算法對μ子的識別能力。假設在整個數(shù)據(jù)集中，μ子的個數(shù)為M，非μ子的宇宙線事件個數(shù)為F。在甄選結果中，被正確甄選為μ子的宇宙線事件個數(shù)為TM，被正確甄選為非μ子的宇宙線事件個數(shù)為TF，則甄選結果的準確率公式：

(9)

甄選結果的靈敏度公式：

(10)

從表2可以看出，當對投影寬度、皮爾森相關系數(shù)和宇宙線事件的像素個數(shù)設定不同閾值時，得到不同的準確率；當特征值設定合適的閾值時，可以有效地從宇宙線事件中甄選μ子。表2特征參數(shù)的選擇需要結合表1中μ子不同特征的均值和標準差。通過統(tǒng)計宇宙線μ子的投影寬度和皮爾森相關系數(shù)，得到描述不同特征的具體數(shù)值參數(shù)以及參數(shù)的波動情況。結合不同特征的參數(shù)波動，對測試數(shù)據(jù)集中的μ子進行甄選，得到最終的結果。

當對宇宙線特征設置合適的閾值時，得到的準確率約為98%，說明該方法可以有效地甄選宇宙線中的μ子。該方法對μ子的識別度很高，但是會將與μ子特征相似、不是μ子的宇宙線事件錯分成μ子。

盡管該算法在數(shù)據(jù)處理過程中存在誤差，但是通過分析μ子的本征特征，結合數(shù)據(jù)集中μ子和其他宇宙線事件的形態(tài)差異，依然可以有效地甄選宇宙線中的μ子。該方法的核心是特征參數(shù)的設定。如果特征值的參數(shù)設定恰當，最后的甄選結果更可觀，也可以通過訓練大量數(shù)據(jù)集的方式調整特征參數(shù)，優(yōu)化方法。根據(jù)表2中最優(yōu)的甄選條件(表2第2行對應的特征參數(shù))，對郭守敬望遠鏡6號紅端的CCD相機在2017年11月拍攝的圖像進行處理。統(tǒng)計這一個月曝光時間為600.0 s的圖像中μ子的個數(shù)。統(tǒng)計結果見圖5。圖中橫坐標為拍攝CCD圖像的修正儒略分鐘數(shù)，縱坐標為圖像中μ子的個數(shù)。圖中藍色豎線對應于第二橫坐標拍攝CCD圖像的日期。圖5直觀地顯示出宇宙線μ子在這一個月的漲落情況。

表2 設定不同特征參數(shù)時μ子的甄選結果

圖5 2017年11月郭守敬望遠鏡6號紅端CCD圖像中宇宙線μ子個數(shù)的統(tǒng)計結果

6 總 結

本文探討了一種簡單有效地對單次曝光CCD圖像中的μ子進行甄選的方法，結合宇宙線在圖像上的形態(tài)特征，達到高效處理數(shù)據(jù)的效果。該方法使用拉普拉斯邊緣檢測法提取宇宙線像素，采用聚類算法將圖像中的宇宙線像素聚類成宇宙線事件，聚類得到宇宙線事件后，通過特征提取和宇宙線事件之間的形態(tài)特征差異對μ子進行甄選。在這個過程中，甄選準確率依賴于特征參數(shù)的設定，可以通過不斷的統(tǒng)計和測試，設定更優(yōu)的參數(shù)以及增加提取的特征，提高甄選結果的準確性。

在獲取訓練和測試樣本集的過程中，通過人眼識別的方式區(qū)分宇宙線μ子和其他宇宙線事件，會導致選擇的數(shù)據(jù)樣本集有一定的偏差，影響之后的統(tǒng)計和甄選結果。

致謝：感謝中國科學院高能物理研究所胡紅波老師的有益討論。