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

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

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

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

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

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

1 快速模擬的基本原理

1.1 宇宙射線μ子

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

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

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

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

(1)

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

(2)

相應(yīng)地，若剩余能量E′=0，此時(shí)E為μ子穿過該待測(cè)物體所需要的最小能量Emin：

(3)

對(duì)于組成成分為單質(zhì)的介質(zhì)，可通過一種簡(jiǎn)單的形式表示μ子的微分能量損失率：

(4)

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

(5)

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

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

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

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

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

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

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

根據(jù)最小能量Emin的分布，可通過能譜模型，計(jì)算得到μ子的通量分布。本文計(jì)算μ子通量使用的Reyna模型的表達(dá)式為[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為μ子的動(dòng)量，動(dòng)量和能量的關(guān)系如下：

(7)

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

2 快速模擬結(jié)果的初步驗(yàn)證

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

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

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

a——俯視圖；b——正視圖藍(lán)色區(qū)域?yàn)榇笞呃龋G色區(qū)域?yàn)閲?guó)王墓室，紅色區(qū)域?yàn)橥鹾竽故?/p>

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

2.2 μ子最小能量的計(jì)算

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

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

2.3 μ子通量的計(jì)算

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

2.4 快速成像模擬結(jié)果分析

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

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

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

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

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

3 結(jié)論

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