丁文婉，劉圓圓,*，王 力，蘇 寧，程建平，韓建武，甄 剛

(1.北京師范大學 核科學與技術學院，生態(tài)環(huán)境部 北京師范大學 錦屏極低輻射本底測量聯(lián)合實驗室，射線束技術教育部重點實驗室，北京 100875；2.北京師范大學 物理學系，北京 100875；3.陜西省文物保護研究院，陜西 西安 710075)

近年來，輻射成像技術越來越多地被應用在地質[1]、考古[2]等領域，成為傳統(tǒng)物探技術的有力補充。對于山體、建筑物等尺寸較大物體的輻射成像，所使用的粒子既要有高穿透性，又要與待測物體有足夠多的相互作用。傳統(tǒng)成像技術所采用的X射線或伽馬射線穿透力弱，無法穿過這些物體，而像中微子這種穿透力極強，幾乎不與待測物體發(fā)生相互作用的射線，也不能較好地實現(xiàn)成像目的[3]。宇宙射線μ子不僅穿透力強，而且能沿其徑跡與待測物體持續(xù)發(fā)生相互作用并不斷損失能量，因此μ子更適用于大尺寸物體的探測[4]。

基于μ子對大尺寸物體進行輻射成像一般使用μ子吸收成像技術。μ子吸收成像主要通過探測宇宙射線μ子穿過待測物體前后的通量變化，結合已知的地形信息，來推測待測物體的密度分布等內部信息。早在1955年，George就使用了μ子吸收成像技術對澳大利亞一處隧道上方的巖石厚度進行了測量[5]。此后，該技術在不同的領域逐步得到應用。在火山成像方面，Nagamine等通過探測近海平面上穿過筑波山的μ子，獲得了該山體的大體輪廓，對μ子吸收成像技術應用于火山學進行了探索[6]。Tanaka等也基于μ子開展了火山成像技術的研究[7]，并在后續(xù)的研究中利用μ子成像獲取了火山內部的物質密度隨時間的變化情況[8-9]，這對實現(xiàn)火山噴發(fā)的監(jiān)控具有重要意義。在考古領域，μ子吸收成像技術最早的應用是在1970年，Alvarez等利用該技術探測了Giza金字塔，驗證了該技術可以用于對金字塔的探測[10]。到了2017年，Morishima等利用多種探測器對胡夫金字塔進行了μ子成像研究，首次發(fā)現(xiàn)了金字塔內部大走廊上方存在一條約30米長的空洞[11]。此外，μ子吸收成像技術在探測地下空洞[12]、礦體[13]等領域也有了一定的研究成果。除了μ子吸收成像，也可以利用μ子與物質相互作用發(fā)生散射的散射角與物質的原子序數(shù)有關的性質，通過測量μ子穿過待測物體后的方向變化，得到物體的密度、結構等信息，這樣的技術被稱為μ子散射成像技術。但由于散射成像探測器的空間布局限制了被測物體的尺寸，散射成像通常用于對小尺寸高原子序數(shù)物體的成像。

在μ子吸收成像的實際應用中，由于天然μ子通量較低，為了獲得更高的探測精度，一般需要花費較長的探測時間，比如像胡夫金字塔這樣的大型待測物體，只有大約1%的μ子能夠穿透金字塔到達探測器，對它的成像需要積累幾個月的數(shù)據(jù)[11]。鑒于此，通過計算機模擬來進行輔助分析是實際測量中重要的一步。一方面，在進行實地測量前，可以對待測對象進行成像模擬，以確定最佳的測量方案。另一方面，通過成像模擬輔助數(shù)據(jù)分析，可以得到待測物體更可靠的結構信息[8]。除此之外，成像模擬也可以用于對尚未應用過μ子成像技術的領域進行可行性研究[14-15]。

在μ子吸收成像模擬中，通常使用蒙特卡羅方法模擬每個μ子在待測物體中的輸運過程，進而根據(jù)探測器收集到的μ子信息得到μ子穿過待測物體后的通量分布。蒙特卡羅模擬通過概率統(tǒng)計的方法進行大量的抽樣實驗，可以比較逼真地描述μ子吸收成像的過程，但是其存在耗時過長的問題，尤其對尺寸較大的物體如金字塔或山體等，模擬所需的時間較長，模擬效率較低。因此研究快速μ子吸收成像模擬，提高模擬效率是非常必要的。目前針對快速μ子吸收成像模擬國內外已有一些研究，如Lesparre等總結了μ子在巖石中衰減的規(guī)律，并導出了簡單的公式，用來快速計算μ子穿過特定巖石厚度后對應的μ子通量[16]。Daniele Carbone等擬合了μ子穿過物體的質量厚度與能量損失的關系式，從而可以計算得到沿每個μ子入射方向上的巖石平均密度的異常[17]。張榮慶等提出了一種基于體素化能損投影的快速成像正演算法，通過連續(xù)采樣每一體素上的能量損失解出每個μ子的能量損失[18]。像這種通過能損計算得到通量分布的成像模擬近期也有一些研究[19]，能夠實現(xiàn)對山體、大型建筑等大尺寸待測對象的快速成像模擬，為探測方案優(yōu)化以及探測結果的預判等提供參考。本文以胡夫金字塔為例，對基于能損計算的快速μ子吸收成像模擬進行了研究，并與蒙特卡羅模擬進行了比較。該模擬方法根據(jù)待測物體的已知結構與探測器的位置信息，通過計算最小可穿透μ子能量得到相對于探測器的不同方向上的μ子通量分布，進而得到成像結果。

1 快速模擬的基本原理

1.1 宇宙射線μ子

如圖1所示，來自外太空的高能宇宙射線與大氣中的原子、分子等發(fā)生碰撞，會產生包括π介子、K介子在內的大量高能粒子，它們再進一步衰變產生μ子等次級粒子。μ子的平均壽命為2.2 μs，在海平面的平均能量為3～4 GeV[20]。μ子通量與其能量以及方向有關，如圖2所示，μ子通量大小隨能量的增加而逐漸減小，隨天頂角θ的減小而逐漸增大。

圖1 宇宙射線μ子產生過程Fig.1 Production process of cosmic ray muon

圖2 海平面測量得到的不同天頂角方向上的μ子動量譜[21]Fig.2 Muon momentum spectra in different zenith angular directions measured at sea level[21]

μ子穿過待測物體時，會與物質發(fā)生相互作用而損失能量。其在物質中的能量損失過程可用式(1)表示，其中：E為μ子穿過待測物體之前的初始能量；E′為μ子穿過不透明度為X0的物質之后的剩余能量。

(1)

不透明度X的定義由式(2)表示，單位為g/cm2；ρ為待測物體密度；x為μ子在待測物體中的位置，dx為微分路徑，對其路徑進行積分，得到：

(2)

相應地，若剩余能量E′=0，此時E為μ子穿過該待測物體所需要的最小能量Emin：

(3)

對于組成成分為單質的介質，可通過一種簡單的形式表示μ子的微分能量損失率：

(4)

式(4)中，右側第1項到第4項分別代表電離、韌致輻射、電子對效應、光核作用造成的能量損失率。當介質為化合物或混合物時，可將該介質視為單質加權平均之后的材料，則其能量損失率可表示為式(5)所示各組分能量損失的加權平均值，其中權重Wi為第i個元素在材料中的質量權重。

(5)

μ子吸收成像中一般通過μ子的通量變化來確認μ子在待測物體中的衰減程度。初始μ子的通量可通過選擇合適的模型來進行描述。μ子通量模型一般可以通過兩種方法獲得，第1種方法是通過蒙特卡羅模擬產生廣延大氣簇射，根據(jù)μ子在大氣中的傳播和衰減，獲得確定海拔下的μ子通量，比較常用的為CRY[22]。第2種方法是通過擬合實驗測得的μ子通量數(shù)據(jù)得到經(jīng)驗公式。例如由Bugaev等提出，并被Gaisser推廣的Gaisser模型[23]；描述垂直方向上μ子通量的Bugaev模型[24]；以及基于Gaisser模型基礎上改進的Tang模型[25]和基于Bugaev模型改進的Reyna模型[26]。以上描述海平面宇宙射線μ子通量分布的模型中，Reyna模型可描述所有天頂角方向的μ子，且對動量范圍為1 GeV/c≤pμ≤2 000/cosθGeV/c的μ子有效。因此，本文采用Reyna模型來描述μ子通量的分布。

1.2 基于能損計算的快速μ子吸收成像模擬

本文基于1.1中的能損計算和通量模型實現(xiàn)快速μ子吸收成像模擬。如圖3所示，首先基于GEANT4建模獲取待測物的不透明度信息，然后利用MATLAB根據(jù)能損計算得到穿過相應不透明度所需的最小能量，最后對通量模型進行積分運算，得到μ子剩余通量分布信息并成像。

圖3 快速μ子吸收成像模擬流程圖Fig.3 Flow chart of fast simulation of muon absorption radiography

首先，基于待測物體的幾何結構或高程數(shù)據(jù)建立模型，并計算相對于測量點各方向上μ子穿過待測物體的徑跡長度，結合相應區(qū)域的材料密度，可以得到不同方向上的待測物體不透明度X的分布。

得到不透明度X的分布后，可計算得到μ子穿過一定不透明度所需要的最小能量Emin的分布。該過程的能損計算有多種方法，例如可以根據(jù)式(4)、(5)求出不同能量的μ子穿過某一介質時對應的能量損失率-dE/dX，對能量損失率進行積分運算，即可得到μ子穿過該物體的能量損失。也可以把能量損失率看成與μ子能量無關的固定值，來計算特定不透明度的能量損失。圖4為根據(jù)Groom等[27]提供的μ子在不同元素、不同材料中的能量損失情況，得到的μ子在碳酸鈣中能量E與能量損失率-dE/dX的關系。圖4中點為Groom等提供的數(shù)據(jù)，線為能量E與能量損失率-dE/dX的擬合曲線，根據(jù)該曲線可計算出μ子穿過一定不透明度時所需要的最小能量。

圖4 μ子在碳酸鈣中的能量損失率Fig.4 Energy loss rate of muon in calcium carbonate

根據(jù)最小能量Emin的分布，可通過能譜模型，計算得到μ子的通量分布。本文計算μ子通量使用的Reyna模型的表達式為[26]：

I(p,θ)=

cos3θc1(pcosθ)-(c2+c3lg(pcos θ)+c4lg2(pcos θ)+c5lg3(pcos θ))

(6)

其中：I(p,θ)為μ子微分通量，單位為cm-2·sr-1·s-1·(GeV/c)-1；θ為天頂角。各參數(shù)的值分別為c1=0.002 53，c2=0.245 5，c3=1.288，c4=-0.255 5，c5=0.020 9。p為μ子的動量，動量和能量的關系如下：

(7)

其中：Eμ為μ子的靜止能量，約0.106 GeV；c為光速，當動量單位以GeV/c表示時，式(7)中的c=1。

2 快速模擬結果的初步驗證

本文以胡夫金字塔為例來實現(xiàn)基于能損計算的快速μ子吸收成像模擬。首先通過GEANT4建模并獲取不同方向上胡夫金字塔的不透明度信息，然后采用上文所述兩種方法計算各方向對應可探測的μ子最小能量，進而根據(jù)Reyna通量模型得到各方向上測得的通量分布。最后將該模擬方法與傳統(tǒng)蒙特卡羅模擬結果進行了比較。

2.1 胡夫金字塔不透明度數(shù)據(jù)的獲取

圖5為基于GEANT4建立的胡夫金字塔模型：金字塔為四棱錐，底面的邊長為230 m，高為139 m，材質為碳酸鈣，密度為2.2 g/cm3。探測器的擺放位置參考了Morishima等[11]對該金字塔的探測方案，探測器置于胡夫金字塔中的王后墓室。圖6為得到的相對于探測器的金字塔不透明度等高線圖，其中不透明度單位為g/cm2，θ為天頂角，范圍為0°到55°，φ為方位角，范圍為0°～360°。

a——俯視圖；b——正視圖藍色區(qū)域為大走廊，綠色區(qū)域為國王墓室，紅色區(qū)域為王后墓室

圖6 胡夫金字塔不透明度X等高線圖Fig.6 Contour map of opacity Xof Khufu pyramid

2.2 μ子最小能量的計算

依據(jù)1.2節(jié)中μ子能損的計算方法，可得到穿過每一對應不透明度X需要的最小能量Emin，如圖7所示。其中，黑色線為使用圖4中擬合的μ子能量損失率曲線計算得到最小能量；藍色線采用當μ子穿過不透明度為1 g/cm2的物質時，取能量損失為1.69 MeV[11]計算μ子的能量損失；作為對比，紅色線是由蒙特卡羅模擬得到不透明度和能量的關系。3種計算方法中，使用恒定的能量損失率來計算最小能量Emin的方法，在穿過不透明度達104g/cm2以上時，得出的能損結果與蒙特卡羅模擬結果差異較大。因此下文中將只使用由μ子能量損失率曲線計算得到的能損結果來進行快速模擬，并與蒙特卡羅模擬結果作比較。

圖7 不同方法得到的μ子穿過不同不透明度所需的最小能量Fig.7 Minimum energy required for muon to pass through different opacities by different methods

2.3 μ子通量的計算

基于1.2節(jié)中的μ子通量模型，根據(jù)μ子穿過金字塔一定不透明度所需的最小能量進行積分運算可以求出各方向上宇宙射線μ子經(jīng)過金字塔后的剩余通量分布。使用Reyna模型計算得到μ子通量成像結果如圖8a所示，通量單位為cm-2·sr-1·d-1。從圖中可看出，金字塔的四條邊緣線組成一交叉圖案，大走廊和國王墓室區(qū)域在圖中也清晰可見，而圖像中左上方的區(qū)域比其他區(qū)域通量高，是探測器所在房間即王后墓室的影響。

2.4 快速成像模擬結果分析

本文以蒙特卡羅模擬的結果作為對比，分別對快速成像模擬得到的成像效果、μ子通量的精確度、模擬時長進行評價。蒙特卡羅模擬方法利用GEANT4實現(xiàn)且使用的金字塔模型、μ子通量模型與快速成像模擬一致。如圖8所示，兩種方法得到的胡夫金字塔內部結構中的大走廊、國王墓室位置以及輪廓的成像結果基本一致。

a——快速μ子吸收成像模擬；b——蒙特卡羅模擬

取圖8中tanθsinφ=0時的μ子通量數(shù)據(jù)，得到一維的通量對比如圖9所示?？梢钥闯?，相較于蒙特卡羅模擬，快速成像模擬得到的通量為一條比較平滑的曲線，且通量值比蒙特卡羅模擬略高。由于本文使用的快速成像模擬是對每個成像的像素點上的不透明度進行數(shù)值計算，得到的通量結果與不透明度有對應的函數(shù)關系，所以一維的通量會是一條平滑的曲線。而蒙特卡羅方法是模擬每個μ子的輸運過程，所以通量結果會有漲落。對整體的通量分布而言，快速成像模擬得到的μ子通量值與蒙特卡羅模擬平均相差小于5%。這是因為快速成像模擬認為μ子沿直線傳播，忽略了傳播過程中的散射。實際上μ子在物體中的徑跡并不是直線，由于散射的作用μ子在待測物體中的徑跡會比快速成像模擬認為的直線徑跡略長，所以實際的能損會偏大，相應得到的μ子通量就會偏小。

圖9 快速μ子吸收成像模擬與蒙特卡羅模擬μ子通量結果對比Fig.9 Comparison of muon flux results between fast simulation of muon absorption radiography and Monte-Carlo simulation

蒙特卡羅模擬需根據(jù)能譜模型進行抽樣得到大量μ子，對每個μ子在每步上的能損和徑跡進行跟蹤模擬，記錄進入探測器的μ子的速度方向。蒙特卡羅模擬要得到足夠清晰的圖像，需要模擬大量的粒子，該過程耗時長，而快速μ子吸收成像模擬只是對每個成像像素上的數(shù)據(jù)進行數(shù)值計算，用時極短。圖8是在同一臺計算機上分別運行快速成像模擬以及蒙特卡羅模擬的結果，其中蒙特卡羅模擬的μ子數(shù)參考了Morishima等[11]實際測量的數(shù)據(jù)量，模擬了等效約60 d的μ子數(shù)。蒙特卡羅模擬用時約20 h，而快速μ子吸收成像模擬需要的時間不到5 min。

3 結論

針對大尺寸物體μ子吸收成像模擬中蒙特卡羅方法耗時過長的問題，本文通過μ子能損計算獲取不透明度與μ子所需最小能量的關系，使用Reyna通量模型計算μ子剩余通量，實現(xiàn)快速μ子吸收成像模擬。以胡夫金字塔為例，使用GEANT4、MATLAB軟件進行算法實現(xiàn)，得到μ子通量成像結果。通過與蒙特卡羅模擬結果的比較驗證了快速μ子吸收成像模擬方法的準確性。初步模擬結果表明，針對類似金字塔的大型物體的成像模擬，快速μ子吸收成像模擬方法所需的模擬時長僅為蒙特卡羅模擬的1/240左右，同時能得到與蒙特卡羅模擬非常接近的模擬結果，是針對大型物體μ子吸收成像模擬的有力工具。在后續(xù)的研究中，我們將嘗試加入散射作用的影響，以獲得更加逼真的模擬結果。