張建鳴 李志偉 劉芳? 李景太 冒鑫 程雅蘋 龐捷 馮鑫茁 倪四道 歐陽曉平 韓然

1) (華北電力大學(xué)核科學(xué)與工程學(xué)院,北京 102206)

2) (北京衛(wèi)星環(huán)境工程研究所,可靠性與環(huán)境工程重點(diǎn)實驗室,北京 100094)

3) (中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院,大地測量與地球動力學(xué)國家重點(diǎn)實驗室,武漢 430077)

4) (西北核技術(shù)研究所,西安 710024)

繆子透射成像方法是一種基于宇宙線繆子穿過目標(biāo)物體前后的通量變化,進(jìn)而獲得其內(nèi)部密度結(jié)構(gòu)的無損探測成像方法.繆子透射成像方法假設(shè)繆子在低原子序數(shù)物質(zhì)中沿直線運(yùn)動,但實際上多重庫侖散射作用會使繆子路徑一定程度上偏離直線,有可能對成像精度造成影響.為此,本文使用Geant4 軟件包開展了繆子透射成像蒙特卡羅模擬,針對數(shù)米尺度多種密度結(jié)構(gòu)的模型,定量分析了打開和關(guān)閉多重庫侖散射物理過程時對成像精度的影響.結(jié)果表明: 對于數(shù)米尺度標(biāo)準(zhǔn)巖石物質(zhì),繆子透射成像方法能夠很好地恢復(fù)內(nèi)部密度異常幾何特征；但多重庫侖散射作用對目標(biāo)物體內(nèi)部區(qū)域近垂向繆子通量造成的偏差可達(dá)5%,而在目標(biāo)物體邊界區(qū)域的偏差可達(dá)13%.因此,需要宇宙線繆子透射成像中考慮多重庫侖散射作用的影響,以獲得數(shù)米尺度目標(biāo)物體更準(zhǔn)確的絕對密度值成像結(jié)果.

1 引言

宇宙線繆子來自于初級宇宙射線與大氣相互作用產(chǎn)生的次級粒子的衰變,具有衰變時間長、能量高、穿透性強(qiáng)(穿透數(shù)百米厚地層)等特點(diǎn)[1].繆子與物質(zhì)主要發(fā)生庫侖散射和電離能損兩種相互作用,由此發(fā)展出兩類成像方法: 散射成像和透射成像[2?5].散射成像技術(shù)利用繆子多重庫侖散射角大小與物質(zhì)的原子序數(shù)相關(guān)的性質(zhì),適用于區(qū)分幾厘米至幾十厘米、密度相似、原子序數(shù)有較大差異的不同材料.2003 年,美國洛斯·阿拉莫斯國家實驗室[6]首次實現(xiàn)了繆子多重庫侖散射成像,可區(qū)分圓柱體鎢塊和支撐用的兩條鋼軌支架.隨后該實驗室先后發(fā)展了多種繆子散射的圖像重建算法,如最鄰近點(diǎn) (point of closest approach,PCA)、最大似然估計 (maximum likehood scattering and displacement,MLSD)、最大后驗估計 (maximum a posteriori,MAP)等算法[7?11],可有效區(qū)分鎢W、鐵Fe、鋁Al 等材料.繆子透射成像方法是基于繆子穿過物體時產(chǎn)生的能量損失與穿過路徑的密度長度相關(guān)的原理,通過分析穿過物體前后的繆子通量分布,即可實現(xiàn)對數(shù)米級至數(shù)百米尺度物體內(nèi)部密度結(jié)構(gòu)成像.1955 年,George[12]使用探測器測量了地表和礦井隧道內(nèi)的繆子通量,根據(jù)繆子通量的衰減測定了隧道上方的巖層厚度.20 世紀(jì)60 年代末,Alvarez[13]首次將繆子應(yīng)用到考古領(lǐng)域,測量了金字塔的內(nèi)部結(jié)構(gòu).2017 年,Morishima 等[14]基于繆子成像技術(shù)發(fā)現(xiàn)并驗證了金字塔內(nèi)存在一個至少30 m 長的巨大密室.繆子透射成像技術(shù)還被廣泛應(yīng)用于火山監(jiān)測,隧道探測等領(lǐng)域.中國空間技術(shù)研究院于2019 年完成了隧道上覆層結(jié)構(gòu)的繆子成像,并分析了觀測時長及數(shù)據(jù)誤差對繆子透射成像的影響[15,16],繼而完成了黑龍江五大連池火山的繆子成像觀測,獲得了火山內(nèi)部高分辨的密度結(jié)構(gòu)圖像[17],顯示出繆子成像在地球科學(xué)研究中的良好潛力.北京師范大學(xué)蘇寧等[18]使用宇宙線繆子對秦始皇陵地宮透射成像開展了模擬研究,利用兩個視角的投影重建了模型的大小和三維位置,成像結(jié)果顯示繆子成像可獲得皇陵內(nèi)部的主要結(jié)構(gòu)特征.

鑒于對物體內(nèi)部密度結(jié)構(gòu)的敏感性,繆子透射成像技術(shù)也被應(yīng)用于小尺度物體內(nèi)部密度結(jié)構(gòu)成像.單一繆子散射成像在垂直方向易造成圖像邊界模糊和偽影,而利用繆子吸收信息能限制目標(biāo)物體的厚度,從而對圖像重建加以修正.中國科學(xué)技術(shù)大學(xué)何偉波[19]使用繆子多模態(tài)成像方法,同時聯(lián)合利用繆子的散射和透射信息,可在較短時間內(nèi)得到質(zhì)量更好的圖像重建結(jié)果.而繆子透射模式和散射模式均可完成不同工況下的反應(yīng)堆堆芯成像[20?22].

針對繆子在介質(zhì)中直線傳播假設(shè)的合理性及多重庫侖散射作用對于繆子透射成像精度的影響,本文開展了數(shù)米級小尺度物體的透射成像模擬,定量分析了多重庫侖散射對小尺度物體透射成像精度的影響,可以為小尺度物體繆子透射成像工作提供理論參考.

2 繆子成像方法

帶電繆子穿過物體時,繆子受到物體內(nèi)原子核的庫侖電場作用偏轉(zhuǎn),偏轉(zhuǎn)的方向與入射方向的夾角稱為散射角.多次庫侖散射平面角分布近似于均值為零的高斯分布,其平面角的均方根 (root mean square,RMS)可表示為

其中,θ為散射角,p,β分別為繆子的動量和相對速度,c為光速,Q為電荷數(shù),X0為物質(zhì)的輻射長度,

其中,ρ為物體密度,Z為物質(zhì)的原子序數(shù),A為相對原子質(zhì)量.在10–3

繆子穿過物體時,同時伴隨著電離作用能量損失過程,平均能量損失率可以描述為

a項與電離損失相關(guān),b項為韌致輻射、μ核反應(yīng)以及電子對產(chǎn)生等輻射能量損失項,?dE為繆子的能量損失,x為質(zhì)量厚度,是密度分布ρ(ε)在繆子路徑L上的線積分:

從(3)式可以看出,繆子的能量損失與穿過物體的密度和路徑有關(guān).繆子穿過一定質(zhì)量厚度的物體后,其剩余能量存在一個最小值Emin,即最小穿透能量,當(dāng)繆子入射能量Eμ

通過分析繆子穿過物體前后通量的變化,可以計算繆子穿透物質(zhì)的質(zhì)量厚度x,再根據(jù)先驗信息,得到物體的長度或密度分布.

對小尺度物體進(jìn)行透射成像時,將待測物體放在探測器中間,待測物體上方探測器用于記錄入射繆子位置信息,下方探測器用于記錄繆子出射位置信息.當(dāng)探測器間不放置任何模型時,通過上層探測器挑選近似垂直入射的繆子,在下層探測器探測到這些近似豎直入射的繆子的通量分布,記為N(o;i,j);探測器間放置模型時,用同樣的方法得到下層探測器的繆子分布,記為N(a;i,j).定義各像素繆子計數(shù)的比值 Ratio(i,j)為

Ratio 越小,說明該像素對應(yīng)豎直方向上的繆子被吸收的越多,該像素對應(yīng)方向上的密度與厚度之積也越大.因此,當(dāng)待測物體厚度一定時,Ratio 的分布可反映物體內(nèi)部的密度結(jié)構(gòu)分布差異.當(dāng)利用繆子進(jìn)行小尺度物體的透射成像時,繆子的傳播路徑示意圖如圖1 所示: 一部分繆子被吸收,一部分繆子穿透物體,穿透繆子由于散射作用,出射方向與位置可能會發(fā)生一定程度的偏移.

圖1 繆子穿透物體的路徑Fig.1.Path of the muon through the object.

本文首先通過模擬不同能量的繆子穿過不同厚度巖石的出射分布,驗證多重庫侖散射模擬過程的正確性,并對不同厚度的巖石-空洞模型進(jìn)行透射成像,對比開關(guān)散射物理過程的成像結(jié)果,從而定量分析多重庫侖散射作用對透射成像精度的影響,完成小尺度物體繆子透射成像和多重庫侖散射作用對透射成像精度影響的研究.

3 模擬環(huán)境搭建和可靠性驗證

本文基于Geant4 (Geometry and Tracking)軟件包,針對以標(biāo)準(zhǔn)巖石為主體的多種結(jié)構(gòu)模型進(jìn)行了宇宙線繆子成像模擬,包括Geant4 散射模塊的可靠性驗證、單能繆子穿過不同厚度巖石的散射影響分析、開關(guān)多重庫侖散射過程對透射成像精度的影響分析等.模擬計算使用的宇宙線繆子源由宇宙射線生成程序(Cosmic-ray shower generator,CRY)產(chǎn)生,通過設(shè)置不同時間、不同緯度、高度繆子的能量、位置、方向等信息,產(chǎn)生相應(yīng)分布的宇宙線繆子源[25],模擬產(chǎn)生的107個繆子能量分布見圖2.

圖2 CRY 隨機(jī)抽樣產(chǎn)生繆子的能量分布Fig.2.Energy distribution of muons generated by random sampling of CRY.

3.1 多重庫侖散射模塊驗證

圖3 為模擬計算模型,仿真空間和探測器材料設(shè)置為真空,待測物體尺寸為3.0 m×3.0 m×10 cm.在Geant4 程序,設(shè)置繆子為單能點(diǎn)源,入射能量為3 GeV,垂直入射到待測物質(zhì),繆子數(shù)量為106個,物理過程包括繆子散射物理過程.計算統(tǒng)計繆子穿過物體后的散射角分布,并將模擬結(jié)果與(1)式計算的散射角RMS 值進(jìn)行對比.

圖3 散射模塊驗證模型(繆子為單能點(diǎn)源,垂直穿過待測物體)Fig.3.Verification model of scattering module (A single energy muon at a point source,passing vertically through the object to be measured).

對比結(jié)果見表1,模擬使用的繆子數(shù)目為106個,統(tǒng)計誤差為取散射角中心98%分布計算RMS 值,得到材料U 的散射角分布RMS 值為31.7 mrad,與(1)式計算值29.30 mrad 相對偏差約8.20%.同時本文多次改變平板材料并計算得到繆子穿過Cu,Al,U,Fe 等物質(zhì)后的散射角分布.鑒于(1)式計算值的誤差約為11%,模擬值均在經(jīng)驗公式計算值的誤差范圍內(nèi),驗證了散射模擬所建立的物理模型及物理過程的可靠性.

表1 3 GeV 的繆子穿過10 cm 厚不同材料的散射角Table 1.Multiple scattering for 3 GeV muons passing through 10 cm of various materials.

3.2 繆子最小穿透能量

對于一定密度長度的物質(zhì),穿透該物質(zhì)需要的繆子最小能量Emin,可通過繆子的連續(xù)阻止本領(lǐng)(continuous slowing down ability,CSDA)[26]表格給出,CSDA 結(jié)合繆子能損公式和實際觀測數(shù)據(jù)給出了繆子能量和不同材料的平均穿透厚度的關(guān)系.本文使用Geant4 模擬不同能量的繆子能穿透的標(biāo)準(zhǔn)巖石的厚度,重新研究了繆子能量和標(biāo)準(zhǔn)巖石的最小穿透厚度的關(guān)系.模擬設(shè)置單能點(diǎn)源繆子在巖石上表面中心位置垂直入射,繆子能量范圍設(shè)置為0.4—20.0 GeV,對于每個單能粒子模型入射106次,并記錄繆子能量沉積為零時豎直方向的位置坐標(biāo).結(jié)果顯示: 500 MeV 的繆子就可穿透1.5—1.8 m厚的標(biāo)準(zhǔn)巖石(圖4(a)),10 GeV 的繆子穿透的標(biāo)準(zhǔn)巖石厚度在15—22 m 間(圖4(c)),同能量的繆子所能穿透的厚度是一個區(qū)間范圍.標(biāo)準(zhǔn)巖石的厚度與最小穿透能量Emin基本呈線性關(guān)系(圖5).

圖4 單能繆子豎直穿透的標(biāo)準(zhǔn)巖石厚度分布 (a) 0.5 GeV;(b) 1.0 GeV;(c) 10.0 GeVFig.4.Standard rock thickness distribution for vertical penetration of muons with different energies: (a) 0.5 GeV;(b) 1.0 GeV;(c) 10.0 GeV.

圖5 不同能量繆子能豎直穿透的標(biāo)準(zhǔn)巖石的平均厚度(誤差棒給出1σ 誤差)Fig.5.Average thickness of standard rocks penetrated vertically by muons of different energies (Error bar is 1σ).

3.3 庫侖散射特征分析和繆子透射成像模型

為了分析低能繆子穿透物體后的多重庫侖散射影響,實驗?zāi)M了不同能量的繆子穿透不同厚度的巖石的出射點(diǎn)分布.計算模型見圖6,探測器中間放置100 m×100 m×h的標(biāo)準(zhǔn)巖石,其厚度h為1,3,5,10,30 m.待測物體上方放置探測器1 和2,待測物體下方緊貼著放置探測器3,用于記錄繆子的出射位置,并結(jié)合探測器4 記錄到的繆子位置信息計算繆子的出射角度,探測器1 和2,探測器3 和4 間距均為1 m.繆子在探測器1 中心位置上方2 cm 位置點(diǎn)垂直入射,繆子入射能量根據(jù)待測物體的厚度h,選取入射能量大于Emin的繆子.

圖6 多重散射特征分析模型 (100 個1.5 GeV 的繆子,垂直穿過3 m 厚度的標(biāo)準(zhǔn)巖石(紅色線條),藍(lán)色為繆子徑跡)Fig.6.Multiple scattering feature analysis model (100 muons of 1.5 GeV vertically passes through a standard rock with a thickness of 3 m (red line),the blue is the muon track).

為了定量研究多重庫侖散射對透射成像的影響,實驗設(shè)置了使用繆子對不同厚度模型透射成像模擬.模型示意圖如圖7 所示,其中待測物體尺寸為1.8 m×1.8 m×0.8h的巖石體,內(nèi)部包括一個尺寸為0.9 m×0.9 m×0.6h的空氣腔體,巖石和空氣腔體中心重合,h=1,3,5 m.巖石模型放置在探測空間坐標(biāo)原點(diǎn),模型上下方各是兩層間距為1 m,尺寸為3 m×3 m×2 cm 的探測器,繆子在探測器1上表面入射.為了提高模擬效率,入射繆子的天頂角限定在3°以內(nèi).模擬粒子物理過程分別包括開啟和關(guān)閉庫侖散射過程的模擬.

圖7 巖石-空腔模型透射成像模擬Fig.7.Rock-cavity model transmission imaging simulation.

4 模擬結(jié)果與討論

4.1 單能繆子穿透巖石的多重庫侖散射影響

能量分別為0.7,1.0,3.0 GeV 的繆子垂直穿透1 m 厚的標(biāo)準(zhǔn)巖石的出射位置分布,如圖8(a)—(c)所示.從圖8(a)可明顯看出,700 MeV 的繆子穿透1 m 厚度的標(biāo)準(zhǔn)巖石,多重庫侖散射使得部分繆子偏離入射方向,粒子出射位置呈現(xiàn)中間密集,越往邊緣越稀疏的分布,多重庫侖散射引起的部分繆子路徑偏移距離最大可達(dá)十幾厘米.而3.0 GeV 的繆子穿透1 m 厚度的標(biāo)準(zhǔn)巖石,多重庫侖散射引起的部分繆子路徑偏移距離在1 cm 內(nèi),散射引起的徑跡偏移很小.通過不同能量點(diǎn)源繆子垂直穿過不同厚度巖石的路徑偏移量(圖9(a))和散射角圖(圖9(b))可以看出,能量較低的繆子散射角和累計路徑偏移距離較大,累計路徑散射偏移距離與累計散射角隨入射繆子能量的增加而迅速減小.如圖9(a)所示,3 GeV 的繆子穿透5 m 厚度的巖石引起的路徑平均偏移距離為18 cm,占穿過路徑長度的3%左右,平均立體散射角接近0.13 rad.總體來說,對于一定厚度的物體,能量接近Emin的繆子(即剩余能量比較小的繆子)的散射作用最為明顯.

圖8 點(diǎn)源繆子豎直穿過1 m 厚標(biāo)準(zhǔn)巖石的出射位置分布 (a) 0.7 GeV;(b) 1.0 GeV;(c) 3.0 GeVFig.8.Distribution of emission positions of point source muons vertically passing through 1 m thick standard rock: (a) 0.7 GeV;(b) 1.0 GeV;(c) 3.0 GeV.

圖9 單能點(diǎn)源繆子穿過不同厚度的標(biāo)準(zhǔn)巖石的多重庫侖散射影響 (a) 路徑平面平均偏移長度(誤差棒給出0.5σ 誤差);(b) 累計散射角Fig.9.Multiple Coulomb scattering effects of single-energy point source muons passing through standard rocks with different thicknesses: (a) Average offset length of the path plane (the error bar is 0.5σ);(b) cumulative scattering angle.

4.2 散射作用對小尺度物體透射成像的影響

以不同厚度模型為例,定量分析了多重庫侖散射作用對繆子透射成像精度影響.圖10 自上而下為巖石厚度為0.8,2.4,4.0 m 的模型,使用的繆子為經(jīng)過篩選的天頂角小于3°的宇宙線繆子,開啟與關(guān)閉多重庫侖散射物理過程的透射成像結(jié)果以及各像素上兩種結(jié)果的差值.從不同厚度的模型成像結(jié)果(圖10(a)—(c)和圖10(d)—(f))可以看出,不論有無散射過程,巖石與內(nèi)部空腔的邊界能夠清晰地看到,繆子透射成像方法可以很好地恢復(fù)待測物體密度異常的空間分布和幾何形態(tài).但多重庫侖散射作用使得穿過高密度區(qū)域邊緣的繆子更容易散射到低密度區(qū)域,增大高密度區(qū)域與低密度區(qū)域的通量衰減差異(圖10(g)—(i)),高密度區(qū)域繆子減少得更多,低密度區(qū)域繆子減少得更少,同時造成密度異常邊界附近出現(xiàn)一些通量增加的假象(圖10(a)—(c)).

圖10 繆子穿過不同厚度模型前后各像素上計數(shù)的比值(俯視圖) (a),(d),(g) 0.8 m;(b),(e),(h) 2.4 m;(c),(f),(i) 4.0 m (黑色實線為巖石模型的外邊界,黑色虛線為巖石與空腔的分界)Fig.10.The ratio of the counts on each pixel before and after the muons pass through the models with different thicknesses(top view): (a),(d),(g) 0.8 m;(b),(e),(h) 2.4 m;(c),(f),(i) 4.0 m (Solid black line is the outer boundary of the rock model,the dashed black line is the boundary between the rock and the cavity).

圖11(a)—(c)和圖11(d)—(f)給出不同厚度下A-A′(圖10(a)—(c)黃色虛線)和B-B′(圖10(a)—(c)紅色虛線)切線經(jīng)過像素上繆子穿過模型前后計數(shù)的比值分布,A-A′切線在巖石和外邊界附近,B-B′切線橫穿了巖石-空腔模型.從圖11(a)—(c)可以看出,不開啟庫侖散射過程,A-A′切線上各像素位置的繆子穿過物體前后的通量比值幾乎為1.0,而開啟庫侖散射過程,A-A′切線上的各像素位置的繆子穿過物體前后的通量比值普遍大于1.0,多重庫侖散射作用使得繆子散射到高密度邊界區(qū)域外.由圖11(d)—(f)可以看出,開啟散射作用,巖石區(qū)域的繆子計數(shù)減少比值比不開散射過程的多,而空腔區(qū)域的繆子計數(shù)減少比值比不開散射過程的整體偏多,散射作用增大了不同密度物質(zhì)的繆子的通量衰減差異,同時使繆子的出射位置分布更加離散.例如,對于穿過2.4 m 厚度的標(biāo)準(zhǔn)巖石,近垂直方向的繆子因多重庫侖散射作用造成的通量衰減占總通量衰減的5%左右,邊界附近的占比達(dá)到了13%,引起的路徑平均偏移量約為6 cm,占穿透路徑長度的2%.

圖11 不同厚度模型切線經(jīng)過像素的計數(shù)比值分布 (a),(d) 0.8 m;(b),(e) 2.4 m;(c),(f) 4.0 m (A-A′ 線在巖石和外邊界附近,B-B′ 切線橫穿了巖石-空腔模型 (圖10 (a)—(c)))Fig.11.Distribution of count ratios of the tangent lines passing through pixels for different thickness models: (a),(d) 0.8 m;(b),(e) 2.4 m;(c),(f) 4.0 m (A-A′ line is near the rock and the outer boundary,and B-B′ tangent traverses the rock-cavity model(Fig.10 (a)–(c))).

4.3 繆子散射角度挑選后的透射成像

增加探測器可以準(zhǔn)確地測出繆子穿透物體前后的方向變化,方向變化較小的繆子直線性更好,有利于提高成像的位置精度.本文使用探測器1 和2 測出繆子的入射天頂角,探測器3 和4 測出繆子的出射天頂角,只使用兩者天頂角相差小于1°的繆子用于成像,篩選掉散射角較大的繆子.對于不同厚度模型的巖石-空氣模型,成像結(jié)果如圖12 所示,可以看出經(jīng)過方向挑選后的繆子成像結(jié)果并沒有出現(xiàn)巖石外邊緣周邊通量增加的假象,空氣與巖石的邊界更加清晰,模型的異常的位置結(jié)構(gòu)更接近于不開啟散射作用的圖像(圖10(d)—(f)),但巖石區(qū)域的繆子通量變化更加明顯,與不做散射角度篩選的成像結(jié)果對比,繆子的通量減少得更多.這是因為在沒有放置模型的情況下,繆子穿過的是空氣,繆子散射角普遍小于1°,能統(tǒng)計到的繆子幾乎沒有變化,但在探測區(qū)域放置巖石-空腔模型后,一部分繆子穿過模型后的散射角超過了1°,沒有被統(tǒng)計到,導(dǎo)致該像素上的N(a;i,j)減小,通量的比值相較于未經(jīng)過散射角篩選的成像結(jié)果減小.

圖12 經(jīng)過繆子散射角篩選后不同厚度模型繆子透射成像結(jié)果(俯視圖) (a) 0.8 m;(b) 2.4 m;(c) 4.0 mFig.12.Transmission imaging results of the model of different thicknesses after the screening of the muon scattering angle(top view): (a) 0.8 m;(b) 2.4 m;(c) 4.0 m.

5 結(jié)論

本文基于Geant4 軟件包,開展了多重庫侖散射作用對小尺度物體透射成像精度影響的研究.研究結(jié)果顯示繆子透射成像方法能夠很好地恢復(fù)待測物體密度異常的空間分布和幾何形態(tài),初步驗證了使用近垂直繆子對小尺度物體密度結(jié)構(gòu)進(jìn)行透射成像的可行性,表明繆子透射成像方法在小尺度物體內(nèi)部密度結(jié)構(gòu)探測中也有很好的成像效果.但是,多重庫侖散射作用一定程度上改變了繆子的運(yùn)動方向和徑跡,影響了高密度區(qū)域與低密度區(qū)域的通量衰減差異,造成在密度差異較大物質(zhì)邊界附近出現(xiàn)通量增加的假象,增大了界面附近的通量衰減差異.通過挑選散射角小于1°的繆子透射成像的方法,可以改善密度差異較大物質(zhì)邊界附近通量增加的假象,但繆子的通量變化更加明顯,對于由通量衰減比值恢復(fù)密度絕對值的準(zhǔn)確性也會產(chǎn)生一定影響.因此,在基于宇宙線繆子進(jìn)行小尺度物體透射成像時,使用通量衰減比值來恢復(fù)物體絕對密度值要注意多重庫侖散射引起的誤差.