繆子對小尺寸中低原子序數(shù)物質(zhì)三維成像技術(shù)的模擬研究

季選韜，羅思遠(yuǎn)，朱 坤，彭肖宇，祝 錦，肖 敏，羅鳳嬌，王曉冬

(南華大學(xué) 核科學(xué)技術(shù)學(xué)院，湖南 衡陽 421001)

宇宙射線繆子是初級宇宙射線與大氣層通過級聯(lián)簇射作用產(chǎn)生的高能帶電粒子[1]，是大自然的本底射線源，質(zhì)量約為電子的207倍，到達(dá)海平面的平均能量為3～4 GeV(最高能量可達(dá)TeV)，服從與cos2θ呈正比的分布[2]，通量約1 cm-2·min-1。因?yàn)榭娮哟┩感詮?qiáng)，所以繆子成像技術(shù)被廣泛應(yīng)用于各類交叉學(xué)科研究中。當(dāng)前宇宙射線繆子成像技術(shù)主要可分為兩大類。第1類是應(yīng)用在考古、地質(zhì)勘察領(lǐng)域的繆子透射成像技術(shù)，1970年，Alvarez利用繆子穿過物質(zhì)前后通量的改變，發(fā)現(xiàn)了金字塔內(nèi)的暗室[3]。隨后該技術(shù)被廣泛應(yīng)用于活火山監(jiān)測[4-6]、煤礦探索[7]、冰川測量[8]、月球陰影調(diào)查[9]、胡夫金字塔內(nèi)未知的密室成像[10]。2020年，國內(nèi)研究團(tuán)隊(duì)通過測量常熟地下隧道內(nèi)外的繆子通量分布，首次成功重建了隧道的圖像[11]。第2類是應(yīng)用于核材料檢測領(lǐng)域的繆子散射成像技術(shù)，2003年美國洛斯阿拉莫斯國家實(shí)驗(yàn)室(LANL)首次提出了基于繆子穿過物質(zhì)前后的散射角信息和POCA(point of closest approach)算法，在短時(shí)間內(nèi)實(shí)現(xiàn)了對高原子序數(shù)Z物質(zhì)的成像[12]。隨后該技術(shù)應(yīng)用在繆子對鈾、钚等重核材料的監(jiān)測[13-14]和集裝箱檢測[15]等核安保方面的相關(guān)研究。為實(shí)現(xiàn)好的成像效果，更多的先進(jìn)探測器被研制[16-19]。福島核電站事故后，美國研究團(tuán)隊(duì)利用散射成像的原理重建了反應(yīng)堆內(nèi)散落的核物質(zhì)圖像[20]，成為散射成像技術(shù)的一次重要應(yīng)用。然而，當(dāng)前的繆子透射成像和散射成像技術(shù)均僅使用了繆子本身的數(shù)量和其徑跡信息，且主要應(yīng)用于高Z或大尺寸的物質(zhì)。實(shí)際上，繆致次級粒子也攜帶了豐富的待測物體信息。Mrdja等提出了基于繆子與繆致次級粒子符合探測技術(shù)的繆子成像技術(shù)[21-23]，在實(shí)驗(yàn)上對牛股骨重建了三維圖像，首次實(shí)現(xiàn)了繆子對中低Z物質(zhì)的成像，拓寬了應(yīng)用領(lǐng)域。

為進(jìn)一步研究基于繆子與繆致次級粒子符合探測技術(shù)的繆子對中低Z物質(zhì)的成像技術(shù)，本文模擬繆致次級粒子的產(chǎn)生和輸運(yùn)行為，研究基于繆子與繆致次級粒子符合探測技術(shù)的符合繆子徑跡篩選，并結(jié)合有限角度成像算法ASD-POCS(adaptive steepest descent projection on convex sets)，實(shí)現(xiàn)繆子對中低Z物質(zhì)的三維成像。

1 繆致次級粒子的研究

1.1 Geant4模擬參數(shù)

1) 探測器系統(tǒng)的幾何建模

本文利用Geant4軟件建立的成像系統(tǒng)的幾何模型如圖1所示。整體由兩部分組成，第1部分是繆子徑跡探測器模塊，其作用是獲取入射繆子的徑跡。該模塊由系統(tǒng)上方3層位置靈敏探測器構(gòu)成，如圖1中灰色部分所示，相鄰兩個(gè)位置靈敏探測器距離5 cm，單個(gè)探測器是尺寸為40 cm×40 cm×1 cm的理想探測器(探測效率100%)，均設(shè)置了100 μm的位置分辨率。第2部分是次級粒子探測器模塊，用于獲取次級粒子的各類信息。該模塊由4個(gè)閃爍體探測器按照四面環(huán)繞的方式排布。單個(gè)探測器的尺寸是50 cm×50 cm×5 cm，并依據(jù)EJ200型號將其定義為C9H10的塑料閃爍體探測器，如圖1黃色部分所示。白色立方體是待測物體，被放置在中央，其尺寸為3 cm×3 cm×3 cm。在模擬中，當(dāng)繆子進(jìn)入系統(tǒng)并與待測物體發(fā)生作用時(shí)，依據(jù)作用截面產(chǎn)生繆致次級粒子并從待測物體中射出，入射繆子的位置信息和繆致次級粒子信息分別被兩個(gè)模塊在1個(gè)符合時(shí)間窗內(nèi)被記錄。

圖1 探測器成像系統(tǒng)的幾何模型示意圖Fig.1 Schematic diagram of geometric model of imaging system

2) 入射粒子及物理過程

Geant4模擬中，粒子源是CRY天然繆子源[24]。模擬時(shí)物理過程被限定為G4Decay-Physics、G4RadioactiveDecayPhysics、G4Em-standardPhysics、G4OpticalPhysics 4種模塊。

1.2 繆致次級粒子產(chǎn)生過程

繆子作為高能的帶電粒子，主要通過4種相互作用損失能量：繆子的電離、繆子的軔致輻射、繆子的電子對效應(yīng)和繆子的光核反應(yīng)。當(dāng)繆子穿過鐵時(shí)，4種相互作用的能量損失隨繆子能量的變化曲線如圖2所示。從圖2可知，對于平均能量為3～4 GeV的天然繆子源，主要是通過電離作用損失能量。由于繆子本身的能量非常高，通過電離作用產(chǎn)生的電子大部分是能量較高的δ電子。因此繆致次級粒子主要包括兩種粒子：δ電子經(jīng)過電磁級聯(lián)簇射作用產(chǎn)生的次級電子和次級γ。

圖2 不同電磁相互作用下繆子的能量損失隨繆子能量的變化曲線[25]Fig.2 Variation of energy loss with muon energy for different electromagnetic interactions[25]

1.3 繆致次級粒子數(shù)量隨繆子穿過的待測物體厚度的變化趨勢

圖3示出了不同材料中原初電離電子、次級γ和次級電子的數(shù)量與繆子穿過待測物體厚度的關(guān)系。其中，原初電離電子是通過Geant4軟件限定了物理過程后獲取的，是一種理想情況下繆致原初電離電子數(shù)據(jù)，次級電子指繆子電離產(chǎn)生的電子，包括從待測物體出射的原初電離電子和原初電離電子繼續(xù)電離產(chǎn)生的電子。為獲得全部的次級粒子，模擬中利用離開待測物體的粒子數(shù)據(jù)取代了閃爍體探測器獲取的，但在圖像重建時(shí)利用的仍是閃爍體探測器獲取的數(shù)據(jù)。圖3a、b、c分別展示了在水、鐵、鉛 3種材料中原初電離電子數(shù)量、次級粒子總數(shù)、次級γ數(shù)量、次級電子數(shù)量隨繆子穿過的待測物體厚度的變化。從圖中可知，在繆子穿過的厚度不超過10 cm時(shí)，原初電離電子、次級電子、次級γ的數(shù)量與繆子穿過的厚度近似呈線性關(guān)系。這是因?yàn)椋蹼婋x電子的產(chǎn)額正比于繆子的沉積能量。繆子能量與物質(zhì)相互作用主要由電離損失能量，因此繆子的能量損失率用-dE/dx表示，其中dE為繆子的微分能量損失，dx為繆子的單位路徑。因?yàn)樵赿x路徑上繆子的能量變化非常小，即dE幾乎保持不變，繆子的能量損失僅與穿過的總厚度相關(guān)，所以繆致的原初電離電子隨繆子穿過的厚度的增加而線性增加。

a——水中原初電離電子、次級γ和次級電子的數(shù)量與繆子穿過待測物體厚度的關(guān)系；b——鐵中原初電離電子、次級γ和次級電子的數(shù)量與繆子穿過待測物體厚度的關(guān)系；c——鉛中原初電離電子、次級γ和次級電子的數(shù)量與繆子穿過待測物體厚度的關(guān)系；d——不同材料中原初電離電子數(shù)量對比；e——不同材料中次級γ數(shù)量對比；f——不同材料中次級電子數(shù)量對比圖3 不同材料中原初電離電子、次級γ和次級電子的數(shù)量與繆子穿過待測物體厚度的關(guān)系 Fig.3 Variation of muonic primary ionized electrons, secondary gamma and secondary electrons counts with thickness of muon passes through in different materials

次級電子數(shù)量遠(yuǎn)小于原初電離電子的數(shù)量，這是因?yàn)榇蟛糠值牡湍茉蹼婋x電子因?yàn)樽晕招?yīng)而被待測物體吸收。從待測物體出射的次級電子可看作由兩部分構(gòu)成，第1部分是產(chǎn)生于待測物體表層(小于0.4 mm)的次級電子，對于水，這部分占比小于25%；第2部分是產(chǎn)生于待測物體內(nèi)部的，來源于高能的原初電離電子。由于產(chǎn)生于表層的次級電子數(shù)量與待測物體的形狀、大小無關(guān)，是大致不變的值。而高能的原初電離電子同樣與待測物體厚度近似呈線性關(guān)系，所以次級電子也呈現(xiàn)類似線性關(guān)系。對于次級γ，其主要來源于原初電離電子的電磁級聯(lián)簇射中的軔致輻射。軔致輻射的能量損失率正比于原初電離電子能量，即正比于繆子的沉積能量，因?yàn)榭娮拥某练e能量僅與穿過的總厚度相關(guān)，從而次級γ數(shù)量也與繆子穿過的厚度近似存在線性關(guān)系。次級γ也有可能來源于繆子直接軔致輻射產(chǎn)生，由于軔致輻射的截面與粒子的質(zhì)量m2呈反比，而繆子的質(zhì)量約為電子的207倍，所以，繆子直接軔致輻射產(chǎn)生的次級γ的數(shù)量相較于原初電離電子的可忽略不計(jì)。考慮到實(shí)驗(yàn)上難僅通過閃爍體探測器區(qū)分次級電子與次級γ，計(jì)算兩種次級粒子的總和與實(shí)際的實(shí)驗(yàn)條件更相符。

圖3d、e、f對比了不同材料產(chǎn)生的原初電離電子、次級γ和次級電子數(shù)量。對于原初電離電子，由于繆子的電離能損與待測物體的Z呈正比，所以水產(chǎn)生原初電離電子數(shù)量少于鐵和鉛；對于次級γ，因?yàn)樵蹼婋x電子的電磁級聯(lián)中的軔致輻射截面與Z2呈正比，從而水產(chǎn)生的數(shù)量遠(yuǎn)小于鐵和鉛；對于次級電子，水的線性程度較好，這是因?yàn)榇渭夒娮又饕a(chǎn)生于水的內(nèi)部，來源于高能的原初電離電子。鐵和鉛的線性程度較差，其次級電子主要產(chǎn)生于待測物體表層。

1.4 繆致次級粒子能譜

圖4示出了水、鐵、鉛的次級粒子能譜。圖4a為次級γ能譜，從圖可知，隨Z的增加，能譜的峰值逐漸右移，這是因?yàn)榭娮拥碾婋x損失截面和軔致輻射的截面均與待測物體的Z呈正比。圖4b是次級電子能譜，從圖可知，鐵和鉛產(chǎn)生的繆致次級電子能譜趨勢變化不大，僅是計(jì)數(shù)的變化，這是因?yàn)殍F和鉛產(chǎn)生的次級電子主要來源于待測物體表層。水的次級電子能譜的峰值低一些，這是因?yàn)樗拇渭夒娮又饕獊碓从趦?nèi)部產(chǎn)生的高能原初電離電子，經(jīng)過自吸收后其能量減小。3種材料的繆致次級電子能譜的峰值均在0.1 MeV以上，在實(shí)驗(yàn)上可與本底環(huán)境噪聲區(qū)分開。水產(chǎn)生的次級γ的能量較低，在實(shí)驗(yàn)上難以從環(huán)境噪聲中區(qū)分，但由于水的繆致次級粒子是以次級電子為主，所以水仍可通過次級粒子的信息成像。

2 成像原理

2.1 符合繆子徑跡

圖5a示出了基于符合繆子徑跡對物體進(jìn)行成像的原理圖。在1個(gè)很短的時(shí)間窗內(nèi)，若繆子徑跡探測模塊與次級粒子探測模塊同時(shí)響應(yīng)，則該入射的繆子徑跡被定義為符合繆子徑跡。這個(gè)很短的時(shí)間窗取決于繆致次級粒子從產(chǎn)生到打到閃爍體探測器上的時(shí)間。為更好探測次級粒子，本研究使用了具有快時(shí)間響應(yīng)的EJ200塑料閃爍體探測器，其上升時(shí)間0.9 ns，衰減時(shí)間2.1 ns，脈沖寬度FWHM=2.5 ns。但實(shí)驗(yàn)上會存在兩種情況導(dǎo)致無法獲取全部的次級粒子，第1種是如果1個(gè)入射繆子產(chǎn)生了多個(gè)次級粒子，且打在了同一片閃爍體探測器上，那么該片閃爍體探測器的時(shí)間分辨難以區(qū)分多個(gè)次級粒子；第2種是若多個(gè)次級粒子打在不同的探測器上，且都產(chǎn)生了可觀測的信號，由于是在1個(gè)時(shí)間窗內(nèi)，所以在實(shí)驗(yàn)上也僅被認(rèn)為是1個(gè)信號。因此在篩選符合繆子徑跡的過程中，如果1個(gè)入射繆子產(chǎn)生了多個(gè)次級粒子，符合繆子徑跡數(shù)量也僅增加1個(gè)。在整個(gè)過程中，假定入射繆子徑跡不發(fā)生改變，即符合繆子徑跡是一直線。

a——繆致次級γ能譜；b——繆致次級電子能譜圖4 繆致次級粒子能譜Fig.4 Energy spectrum of muonic secondary particles

a——符合繆子徑跡定義的示意圖；b——沿某一方向下符合繆子徑跡的數(shù)量隨沿該方向下待測物體厚度的變化圖5 符合繆子徑跡的定義與特性Fig.5 Definition and characteristic of coincide muon trajectory

結(jié)合圖3的結(jié)論，在繆子穿過的厚度不超過10 cm的范圍內(nèi)，次級γ和次級電子的數(shù)量之和隨繆子穿過的厚度呈線性增長，可推出，當(dāng)沿某一方向下的待測物體厚度增加時(shí)，沿該方向入射的繆子則更易于產(chǎn)生繆致次級粒子，從而沿該方向的符合繆子徑跡數(shù)量更多。利用Geant4軟件模擬了沿某一方向下符合繆子徑跡的數(shù)量與沿該方向下待測物體厚度的關(guān)系，結(jié)果如圖5b所示。水、鐵和鉛3種材料的某一方向下符合繆子徑跡數(shù)量均隨沿該方向下待測物體厚度的增加而近似線性增加。但符合繆子徑跡的線性關(guān)系相較于次級粒子數(shù)量的差一些，特別是鉛的線性關(guān)系。這是因?yàn)樵诋a(chǎn)生多個(gè)次級粒子的符合事件中，符合繆子徑跡數(shù)量僅增加1個(gè)，這使得次級粒子數(shù)量的線性關(guān)系沒有被充分利用，這種情況在次級粒子數(shù)量最多的鉛中最明顯。所以當(dāng)獲取了沿某一方向下符合繆子徑跡的數(shù)量，則可反推出沿該方向下待測物體的厚度，如果獲取的徑跡足夠多，則可對待測物體重建三維圖像，此原理與醫(yī)學(xué)成像的斷層成像原理類似。但由于繆子的天頂角服從cos2θ呈正比的分布，即沿水平或接近水平入射方向的天然繆子較少，從而符合繆子徑跡被限定在一有限角度的天頂角范圍內(nèi)，因此需結(jié)合有限角度的醫(yī)學(xué)成像ASD-POCS算法重建圖像。

2.2 ASD-POCS算法

ASD-POCS算法[26]是一種將聯(lián)合代數(shù)重建算法[27](SART)與總變分方法[28](TV)相結(jié)合的算法，其主要應(yīng)用于醫(yī)學(xué)成像中稀疏角度投影問題或有限角度投影問題，與天然繆子角度不完善的情況符合。該算法通過受約束的TV最小化得到離散圖像f，如下式所示：

f*=argmin‖f‖TVsubject to

Af=p,f≥0

其中：f為參與迭代的離散待重建圖像；A為系統(tǒng)矩陣；p為投影數(shù)據(jù)。ASD-POCS算法的成像過程可分為以下3個(gè)步驟。

第1步，創(chuàng)建投影數(shù)據(jù)p，和迭代停止條件ε，如下式所示：

投影數(shù)據(jù)p相當(dāng)于密度函數(shù)μ(x,y,z)沿某射線做線積分，應(yīng)用到繆子成像領(lǐng)域時(shí)，可用沿某方向的符合繆子徑跡密度作為投影數(shù)據(jù)p，而迭代停止條件ε的取值取決于投影數(shù)據(jù)質(zhì)量以及配置參數(shù)。

第2步，通過SART將所有投影數(shù)據(jù)進(jìn)行迭代，不斷減少重建圖像的估計(jì)投影數(shù)據(jù)與真實(shí)投影數(shù)據(jù)的殘差，如下式所示：

當(dāng)獲得180°的投影數(shù)據(jù)時(shí)，可直接通過迭代將圖像完整重建出來。但由于探測器的面積有限，而符合繆子徑跡相當(dāng)于一同時(shí)穿過繆子徑跡探測器和待測物體的直線，所以能探測到的符合繆子徑跡會被限定在小于180°的有限角度內(nèi)，如圖6所示，所以僅通過第2步無法獲得完整的圖像。

圖6 在有限角度下收集符合繆子徑跡的示意圖Fig.6 Schematic of collection of coincide muon trajectories at limited angle range

第3步，采用最速下降算法(SD)使TV最小化，如下式所示：

min TV(f) subject toAf=p,f≥0

完成第3步后，再回到第2步，這樣交替進(jìn)行第2步和第3步，直到滿足迭代停止條件，最終獲得一離散圖像f的最優(yōu)解。

3 繆致次級粒子的研究

3.1 不同材料的成像效果對比

利用ASD-POCS算法和符合繆子徑跡信息重建的三維圖像如圖7所示，從左到右3個(gè)立方體填充的材料為鉛、鐵、水，體積均為3 cm×3 cm×3 cm。入射的繆子數(shù)量為5×107個(gè)，換算為天然宇宙射線繆子的通量約為20 d的實(shí)際成像時(shí)間。由于鐵和鉛產(chǎn)生的符合繆子徑跡數(shù)量遠(yuǎn)大于水，對3種物質(zhì)同時(shí)成像可能導(dǎo)致重建圖像中的水接近透明。為更好比較不同材料的成像效果，Geant4模擬時(shí)是在相同條件下分別進(jìn)行模擬再匯于同一三維矩陣中。當(dāng)待測物體為水時(shí)，次級電子對成像的貢獻(xiàn)達(dá)到90%。當(dāng)待測物體為鐵和鉛時(shí)，次級γ對成像的貢獻(xiàn)占主要，分別為72.5%和82.1%。從圖7可知以下結(jié)論。

1) 鐵的重建圖像內(nèi)部無空洞，噪聲點(diǎn)最少，成像效果最好。

2) 水的內(nèi)部存在較少中空，噪聲點(diǎn)較多，這是由兩個(gè)原因?qū)е拢皇且驗(yàn)榉峡娮訌桔E數(shù)量較少，從而SART迭代的驅(qū)動矩陣的賦值略小，再經(jīng)過ASD-POCS算法迭代后仍存在大量離散的噪聲點(diǎn)，不過該問題可通過增加成像時(shí)間來解決；二是因?yàn)椴糠执渭夒娮蛹挟a(chǎn)生于待測物體表層(小于0.4 mm)，導(dǎo)致邊緣一層像素的數(shù)據(jù)量大，從而導(dǎo)致內(nèi)部略微中空。模擬結(jié)果可發(fā)現(xiàn)，水的成像效果仍好于鉛。

3) 鉛的內(nèi)部中空的部分相較另兩種材料較大，噪聲點(diǎn)也較多，這是因?yàn)樵阢U中心產(chǎn)生的次級粒子更難逃逸出鉛塊，所以穿過鉛中心部分的符合繆子徑跡較少，導(dǎo)致成像的中空和大量的噪聲點(diǎn)。因此，該技術(shù)適合于對低、中序數(shù)物質(zhì)成像，對中Z物質(zhì)的成像效果最好，高Z物質(zhì)的重建圖像中間可能會有中空。

3.2 復(fù)雜待測物體的成像效果

本文構(gòu)建了“USC”的幾何模型，其高度均為2 cm，字母“U”和“C”的凹槽寬度為2 cm，字母“S”的凹槽寬度為0.8 cm，3個(gè)字母全部填充材料為鐵，入射繆子數(shù)量同樣是5×107，重建的圖像如圖8所示。從圖中可知，重建的圖像可清晰區(qū)分“USC”，尺寸與位置與幾何建模匹配，且可區(qū)分“S”中0.8 cm的凹槽。但成像的不足也很明顯，3個(gè)字母均有微小的傾斜，且字母周圍的噪聲點(diǎn)較多，這是物理模型與算法的數(shù)學(xué)模型不匹配導(dǎo)致的，算法中假定沿某一方向的符合繆子徑跡的數(shù)量與沿該方向下的待測物體厚度呈線性關(guān)系，但實(shí)際的過程中僅是近似呈線性關(guān)系，這導(dǎo)致ASD-POCS算法最后的收斂結(jié)果噪聲點(diǎn)較多，甚至導(dǎo)致了重建圖像的輕微傾斜。可通過優(yōu)化算法的數(shù)學(xué)模型改善噪聲點(diǎn)，實(shí)現(xiàn)更好成像精度。

a——鉛、鐵、水的三維成像的俯視圖；b——鉛、鐵、水的三維成像的三維圖圖7 不同材料的三維成像結(jié)果對比圖Fig.7 Comparison chart about 3D imaging result of different materials

a——“USC”重建圖像的俯視圖；b——“USC”重建圖像的三維圖圖8 幾何模型“USC”的三維成像結(jié)果Fig.8 3D imaging result of "USC" geometric model

4 小結(jié)

本文研究了基于繆子與繆致次級粒子符合探測技術(shù)與有限角度成像算法的繆子對小尺寸中低Z的三維成像技術(shù)，包括模擬次級粒子的特性、研究基于符合繆子徑跡的成像原理、編寫ASD-POCS有限角度的三維成像算法。結(jié)果表明，次級粒子由次級γ和次級電子組成；其能譜均大部分大于0.1 MeV，可與本底區(qū)分；兩種次級粒子數(shù)量之和與繆子穿過的待測物體厚度近似呈線性關(guān)系，所以基于繆子及繆致次級粒子的符合探測技術(shù)建立的某一方向下的符合繆子徑跡數(shù)量也與該方向下的待測物體厚度近似呈線性關(guān)系，以此為成像原理再結(jié)合ASD-POCS成像算法進(jìn)行三維成像。重建的三維圖像表明，該技術(shù)更適合于對小尺寸中、低Z物質(zhì)成像，對小尺寸中Z物質(zhì)的成像效果最好，高Z的重建圖像內(nèi)部會有中空。圖像中可區(qū)分字母“S”里0.8 cm的凹槽，通過進(jìn)一步的算法優(yōu)化可實(shí)現(xiàn)更好成像精度。