蘇寧 劉圓圓? 王力 程建平

1) (北京師范大學(xué)核科學(xué)與技術(shù)學(xué)院,射線束技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,北京 100875)

2) (北京師范大學(xué)物理學(xué)系,北京 100875)

1 引言

宇宙射線是來自于太空的高能射線.對(duì)宇宙射線的觀測(cè)研究[1,2]不僅推動(dòng)了天體物理、粒子物理等基礎(chǔ)研究的發(fā)展,同時(shí)研究人員提出可以利用宇宙射線中的繆子(μ子)能量高、穿透能力強(qiáng)的特點(diǎn),通過探測(cè)穿透待測(cè)物體的宇宙射線μ子的通量、角分布等信息實(shí)現(xiàn)對(duì)待測(cè)物體的成像,這一技術(shù)被稱為μ子成像技術(shù).μ子成像技術(shù)是一種輻射成像技術(shù),根據(jù)所利用的μ子與物質(zhì)相互作用性質(zhì)的不同,μ子成像技術(shù)可分為兩類:散射成像和吸收成像.散射成像利用宇宙射線μ子穿透物質(zhì)時(shí)發(fā)生的多重庫(kù)侖散射的散射角大小與材料的原子序數(shù)有關(guān)的特性,適用于對(duì)小型、高密度、高原子序數(shù)的物體的成像[3].吸收成像利用宇宙射線μ子穿透物質(zhì)時(shí)的能量損失率與穿透路徑上物質(zhì)的密度、厚度有關(guān)的特性來對(duì)物體的密度結(jié)構(gòu)成像.由于宇宙射線μ子能量很高、穿透能力強(qiáng),μ子吸收成像技術(shù)可以對(duì)尺度達(dá)千米量級(jí)的物體成像[4].

具體來說,對(duì)μ子吸收成像的研究最早可追溯至1955 年,George[5]使用μ子探測(cè)器對(duì)澳大利亞一處地下工廠設(shè)施上方的巖層厚度進(jìn)行測(cè)量.20 世紀(jì)60 年代末期,Alvarez 領(lǐng)導(dǎo)的實(shí)驗(yàn)組[6]使用μ子探測(cè)器對(duì)金字塔中可能存在的墓室進(jìn)行搜尋,這是此技術(shù)在考古領(lǐng)域的首次應(yīng)用.此后近三十年的時(shí)間里,對(duì)μ子吸收成像技術(shù)的研究進(jìn)展緩慢,主要是對(duì)這一技術(shù)應(yīng)用于火山[7]、地下洞穴[8]、礦體探測(cè)[9]的模擬分析和初步的實(shí)驗(yàn)應(yīng)用.進(jìn)入21 世紀(jì)后,得益于探測(cè)器技術(shù)的進(jìn)步,μ子吸收成像技術(shù)得到了快速發(fā)展.日本[4]、意大利[10]等國(guó)家的科學(xué)家利用μ子吸收成像技術(shù)對(duì)多座火山的內(nèi)部結(jié)構(gòu)進(jìn)行了較好的成像,并開展了對(duì)火山內(nèi)部密度結(jié)構(gòu)動(dòng)態(tài)監(jiān)測(cè)的研究[11],嘗試結(jié)合重力數(shù)據(jù)進(jìn)行反投影成像以提高成像精度[12],成像不再局限于二維,而是發(fā)展到三維的密度結(jié)構(gòu)成像[13].在對(duì)金字塔[14]、地下洞穴[15]、礦體[13]等不同類型觀測(cè)目標(biāo)的成像方面也有了一些重要的研究成果.特別是2017 年,Morishima 等科學(xué)家[14]分別將不同的μ子探測(cè)器先后放置在胡夫金字塔內(nèi)部和外側(cè)進(jìn)行探測(cè),發(fā)現(xiàn)并驗(yàn)證了大走廊上方存在一個(gè)至少30 m 長(zhǎng)的未知墓室.這一研究結(jié)果表明μ子吸收成像技術(shù)在對(duì)大型文物的無損探測(cè)方面具備極強(qiáng)的應(yīng)用潛力.而我國(guó)作為一個(gè)歷史悠久的文明古國(guó),現(xiàn)存大量文物遺跡亟待考古研究.針對(duì)帝王陵墓等一些大型文物內(nèi)部結(jié)構(gòu)的無損探測(cè),目前考古學(xué)中傳統(tǒng)的地球物理方法存在著一定的局限性,例如電法、磁法易受地面環(huán)境的干擾;地震波法成本較高[16];探地雷達(dá)法的探測(cè)深度范圍在幾米到十幾米之間浮動(dòng)[17],可探測(cè)深度相對(duì)較淺.且每一種方法僅對(duì)特定的一兩種物性敏感,存在多解性的問題[18].將μ子吸收成像技術(shù)應(yīng)用于考古領(lǐng)域,可以成為對(duì)傳統(tǒng)地球物理方法的重要補(bǔ)充.

目前國(guó)內(nèi)對(duì)μ子吸收成像的研究較少,且集中在地質(zhì)探測(cè)方面,在μ子吸收成像應(yīng)用于大型文物成像方面的研究幾乎空白.本文以秦始皇陵地宮為研究對(duì)象,利用蒙特卡羅模擬的方法研究μ子吸收成像應(yīng)用于帝王陵墓無損探測(cè)的可行性.

2 μ子吸收成像算法

2.1 μ子與物質(zhì)相互作用原理

宇宙射線μ子是一種輕子,它主要由原初宇宙射線與大氣層中的原子相互作用產(chǎn)生的π 介子和K 介子的衰變產(chǎn)生,能量很高,在相對(duì)論效應(yīng)的作用下能夠在發(fā)生衰變之前到達(dá)海平面.海平面觀測(cè)到的宇宙射線μ子的平均能量在3—4 GeV,平均通量約為1 cm—2·min—1,是到達(dá)海平面的宇宙射線的主要成分[19].到達(dá)地面的μ子的天頂角θ可以是0°—90°間的任意值,μ子通量隨θ的增大而減小.圖1(a)為實(shí)驗(yàn)中測(cè)得的不同天頂角方向入射的海平面宇宙射線μ子動(dòng)量分布,圖1(b)為天頂角θ、方位角φ示意圖.

圖1 (a) 實(shí)驗(yàn)測(cè)量得到的不同方向上的海平面μ子通量[20];(b)探測(cè)器探測(cè)到的μ子的方位角φ 和天頂角θ,其中xOy為水平面Fig.1.(a) Sea-level muon flux at different zenith angles measured in experiment[20];(b) zenith angle θ and azimuth angle φ of the muon detected by a detector.The xOy plane represents for horizontal plane.

μ子穿過物質(zhì)會(huì)通過電離、軔致輻射、電子對(duì)產(chǎn)生、核相互作用等方式損失能量,其平均能量損失率可一般地表達(dá)為[4]

其中E代表μ子能量,為μ子穿過物質(zhì)的密度對(duì)穿透長(zhǎng)度的積分,a對(duì)應(yīng)電離導(dǎo)致的能量損失率,bE對(duì)應(yīng)其他作用導(dǎo)致的能量損失率.a,b與μ子能量和所穿透物質(zhì)的種類有關(guān).

2.2 圖像重建算法

由(1)式可知,μ子穿透一定量的物質(zhì)后平均損失的能量 ΔE可表示為

其中X0為μ子穿透物質(zhì)的密度對(duì)穿透路徑的積分.X0越大,則μ子穿透物質(zhì)平均損失的能量越多,能穿透物質(zhì)的μ子越少,故可根據(jù)探測(cè)器探測(cè)到的穿透待測(cè)物體的μ子數(shù)量的多少推斷μ子穿透物質(zhì)的X的大小.輔助一些已知信息,即可計(jì)算未知量.如已知穿透路徑上的平均密度時(shí),可計(jì)算出路徑長(zhǎng)度;已知穿透路徑長(zhǎng)度時(shí),可計(jì)算出穿透路徑上的平均密度.

本文基于對(duì)待測(cè)物體掌握的先驗(yàn)知識(shí)(待測(cè)物體幾何尺寸、部分區(qū)域的密度分布),假設(shè)待測(cè)物體內(nèi)部存在密度分布不同于先驗(yàn)知識(shí)的感興趣區(qū)域(region of interest,ROI),所使用的μ子吸收成像的圖像重建算法可概括為以下3 步:

第一步,獲取μ子通量.假設(shè)Nμ(θ,φ) 表示探測(cè)器在 (θ,φ) 方向上探測(cè)到的μ子數(shù),則探測(cè)器探測(cè)到的μ子通量Φ0(θ,φ) 為

其中S為探測(cè)器面積,Ω為立體角,t為測(cè)量時(shí)間.

第二步,與先驗(yàn)知識(shí)對(duì)比,獲取ROI 邊界的二維角坐標(biāo)信息.利用已知的待測(cè)物體的幾何尺寸和密度信息建立參考模型,模擬得到宇宙射線μ子穿透參考模型后剩余的μ子通量Φs(θ,φ),定義μ子通量差f(θ,φ)=Φ0(θ,φ)-Φs(θ,φ) .當(dāng)f(θ,φ)/=0 時(shí),代表 (θ,φ) 方向上的密度分布與利用已知信息構(gòu)建的模型的密度分布不同,(θ,φ) 位于ROI 內(nèi).若f(θ,φ)＜0,則代表 (θ,φ) 方向上存在密度偏高的結(jié)構(gòu)(如高密度礦物),若f(θ,φ)＞0 ,則代表(θ,φ)方向上存在密度偏低的結(jié)構(gòu)(如空穴).可根據(jù)f(θ,φ) 發(fā)生突變處對(duì)應(yīng)的 (θ,φ) 確定ROI 邊界的角坐標(biāo).

第三步,重構(gòu)ROI 的三維圖像.探測(cè)器在單個(gè)測(cè)量點(diǎn)得到的數(shù)據(jù)僅能給出ROI 相對(duì)該測(cè)量點(diǎn)的二維角坐標(biāo)信息,為了獲得ROI 的三維信息,需要有兩個(gè)或兩個(gè)以上不同位置的測(cè)量點(diǎn)的數(shù)據(jù).以兩個(gè)測(cè)量點(diǎn)為例,當(dāng)ROI 為簡(jiǎn)單的幾何體(如長(zhǎng)方體)時(shí),可以很方便地將兩個(gè)測(cè)量點(diǎn)觀測(cè)到的ROI 邊界角坐標(biāo)對(duì)應(yīng)起來.如圖2,在測(cè)量點(diǎn)1 觀測(cè)到ROI 的面abcd的角坐標(biāo)為 (θ1,φ1),在測(cè)量點(diǎn)2 觀測(cè)到面abcd的角坐標(biāo)為 (θ1′,φ′1),則可結(jié)合兩個(gè)測(cè)量點(diǎn)的位置信息,根據(jù)幾何關(guān)系計(jì)算出面abcd的三維坐標(biāo),以此類推,最終重建出ROI 的幾何大小和三維位置信息.

圖2 測(cè)量點(diǎn)與ROI 之間的幾何關(guān)系示意圖Fig.2.Geometric relationship between viewpoint and ROI.

3 模擬系統(tǒng)建立

本文基于Geant4 模擬宇宙射線μ子在秦始皇陵中的輸運(yùn)過程.Geant4 是歐洲核子中心使用C++語(yǔ)言開發(fā)的一款模擬粒子輸運(yùn)過程的蒙特卡羅軟件包,能提供構(gòu)建幾何模型、跟蹤粒子徑跡、模擬粒子與物質(zhì)的相互作用等功能[21].本文模擬系統(tǒng)的建立主要包括3 個(gè)方面:秦始皇陵地宮模型、μ子探測(cè)器、海平面宇宙射線μ子源.

3.1 秦始皇陵地宮模型

本文根據(jù)使用綜合地球物理方法得到的秦始皇陵考古數(shù)據(jù)建立了兩個(gè)秦始皇陵地宮簡(jiǎn)化模型[22].其中模型1 為待測(cè)物體模型,包括了封土堆、土地和地宮,其中地宮由細(xì)夯土墻、宮墻、回填夯土和墓室組成;模型2 為參考模型,僅包括封土堆和土地.模型幾何及關(guān)鍵結(jié)構(gòu)的尺寸如圖3 所示.表1 列出了地宮模型中不同區(qū)域的材質(zhì)和密度,其中空氣定義為70%的氮?dú)夂?0%的氧氣,黃土的組分定義參考了文獻(xiàn)[23].

表1 秦始皇陵地宮模型材質(zhì)及密度定義表[22]Table 1.Material and density definition table of the Qinshihuang Mausoleum model[22].

3.2 μ子探測(cè)器模型

本文采用了理想μ子探測(cè)器模型,參考常見的μ子探測(cè)器幾何結(jié)構(gòu),將探測(cè)器模型表面積設(shè)為1 m × 1 m,其作用為提取到達(dá)探測(cè)器區(qū)域的μ子的速度方向信息.選取探測(cè)器位置時(shí)考慮了以下4 個(gè)因素:1)探測(cè)器位置的海拔高度需要低于地宮區(qū)域,且避開已知的文物埋藏區(qū)域;2)由于μ子通量隨天頂角增大而顯著減少,為保證μ子的計(jì)數(shù)率,應(yīng)避免探測(cè)器對(duì)待測(cè)區(qū)域探測(cè)方向的天頂角過大;3)目標(biāo)探測(cè)區(qū)域要盡可能離探測(cè)器近以提高計(jì)數(shù)率;4)目標(biāo)探測(cè)區(qū)域內(nèi)不同密度的待探測(cè)結(jié)構(gòu)在探測(cè)方向上盡量沒有重疊.為了給出成像目標(biāo)的三維信息,選定兩個(gè)測(cè)量點(diǎn),如圖3(a4)所示.測(cè)量點(diǎn)1 選取在μ子計(jì)數(shù)率最高的地宮正下方,探測(cè)器水平放置,埋深80 m;測(cè)量點(diǎn)2 選取在非垂直方向的地宮邊緣下方,位于地宮中心正南方80 m 處,探測(cè)器埋深89.5 m,探測(cè)器表面與水平面成40°夾角,使探測(cè)平面法線方向指向墓室中心,從而增大此方向上的μ子計(jì)數(shù)率.

圖3 秦始皇陵模型示意圖 (a) 模型1 內(nèi)部結(jié)構(gòu)示意圖;(a1) 模型1 俯視圖;(a2) 模型1 正視圖;(a3) 模型1 剖面3 示意圖;(a4) 模型1 剖面1 示意圖;(b) 模型2 示意圖(無內(nèi)部結(jié)構(gòu));Fig.3.Model of Qinshihuang Mausoleum:(a) Inner structure of Model 1;(a1) top view of Model 1;(a2) front view of Model 1;(a3) profile 3 of Model 1;(a4) profile 1 of Model 1;(b) Model 2 (no inner structure).

3.3 μ子源產(chǎn)生

模擬中μ子源根據(jù)海平面μ子能譜經(jīng)驗(yàn)公式抽樣產(chǎn)生.關(guān)于海平面μ子能譜分布的模型有很多種,例如Gaisser[24],Reyna[25],以及Smith 和Duller[26]等提出的模型.比較這些模型發(fā)現(xiàn),Reyna 總結(jié)的海平面μ子能譜公式能在1 GeV/c ≤p≤2000/cosθGeV/c的范圍內(nèi)擬合全部天頂角范圍內(nèi)的μ子能譜分布,適合于μ子吸收成像中μ子源的抽樣[27].其表達(dá)式為

I(p,θ) 為海平面μ子微分通量,單位為(GeV/c)-1·cm-2·sr-1·s-1;p為μ子動(dòng)量,單位為GeV/c;θ為μ子速度方向天頂角大小;c1=0.00253 ,c2=0.2455,c3=1.288,c4=-0.2555 ,c5=0.0209 .

模擬中使用的μ子源通過對(duì)(4)式在p=30—1000 GeV/c,θ=0°—70°的范圍內(nèi)抽樣產(chǎn)生,以減少對(duì)無法穿透土地到達(dá)探測(cè)器的低能μ子的模擬,所使用的μ子源的動(dòng)量和天頂角分布如圖4.

圖4 根據(jù)Reyna 公式抽樣產(chǎn)生的1000 萬(wàn)個(gè)μ子的動(dòng)量和天頂角分布 (a) μ子數(shù)量隨μ子動(dòng)量變化分布;(b) μ子數(shù)量隨μ子速度方向的天頂角變化分布Fig.4.Momentum spectrum and zenith angle distribution of the 10 million muons sampled by Reyna formula:(a) Momentum spectrum of the sampled muons;(b) zenith angle distribution of the sampled muons.

4 模擬結(jié)果

4.1 墓室二維成像結(jié)果

對(duì)模型1 和模型2 均模擬了等效于實(shí)際測(cè)量一年的μ子量,對(duì)兩個(gè)測(cè)量點(diǎn)在模型1 和模型2 測(cè)得的μ子通量作差得到f(θ,φ) 的分布,如圖5 所示,圖中灰度值表示f(θ,φ)(單位:cm-2·sr-1·d-1) .

圖5(a)中間的白色矩形區(qū)域f(θ,φ) ＞ 0,說明此方向上存在密度偏低的結(jié)構(gòu),即墓室,根據(jù)白色矩形關(guān)于圖像中心的對(duì)稱性可知測(cè)量點(diǎn)1 位于墓室正下方;周圍的深色矩形區(qū)域f(θ,φ) ＜ 0,說明此方向上存在密度偏高的結(jié)構(gòu),即墓室周圍的墻體; t anθy=0 附近深色矩形區(qū)域邊緣處向內(nèi)凹陷的部分對(duì)應(yīng)于墻體的不連續(xù)部分.圖5(b)中的白色扇形區(qū)域f(θ,φ) ＞ 0,對(duì)應(yīng)于墓室;周圍的深色弧形區(qū)域f(θ,φ) ＜ 0,對(duì)應(yīng)于墻體;弧形區(qū)域與白色扇形區(qū)域之間的區(qū)域f(θ,φ) ＜ 0,對(duì)應(yīng)于地宮開挖范圍內(nèi)回填的夯土.對(duì)比圖5(a)和圖5(b)可知,相對(duì)于測(cè)量點(diǎn)1,測(cè)量點(diǎn)2 除可以顯示墻體、墓室外,還可以看出地宮開挖范圍回填夯土所在區(qū)域,這是因?yàn)樵跍y(cè)量點(diǎn)1 的探測(cè)方向上墻體總是與回填夯土區(qū)域有重疊,導(dǎo)致兩者的位置在圖5(a)所示的角度投影圖上也是相互重疊的.

為了便于墓室三維重建,選擇圖3(a1)所示剖面1、剖面2 處的f(θ,φ) 進(jìn)行分析,從而定位墓室墻壁在這兩個(gè)剖面處的角坐標(biāo).如圖5 所示,虛線表示剖面方向,圓點(diǎn)表示剖面處的墻壁的角坐標(biāo)位置,其中點(diǎn)B(23.75°,270°)和E(38.8°,0°)對(duì)應(yīng)于剖面1 處的南墻,點(diǎn)A(25.64°,90°)和F(58.0°,0°)對(duì)應(yīng)于剖面1 處的北墻,點(diǎn)C(35.75°,180°)對(duì)應(yīng)于剖面2 處的西墻,D(37.23°,0°)對(duì)應(yīng)于剖面2處的東墻.

圖5 兩個(gè)測(cè)量點(diǎn)得到的 f (θ,φ) 的二維投影圖 (a)測(cè)量點(diǎn)1 的 f (θ,φ) 投影圖,其中,t anθx=tanθcosφ ,tanθy=tanθsinφ ;(b)測(cè)量點(diǎn)2 的 f (θ,φ) 投影圖Fig.5.Two-dimensional projection of f (θ,φ) obtained at viewpoint 1 and 2:(a) Distribution of f (θ,φ) obtained at viewpoint 1,where the t anθx=tanθcosφ ,t anθy=tanθsinφ ;(b) distribution of f (θ,φ) obtained at viewpoint 2.

4.2 墓室三維重建結(jié)果

墓室的大小和位置是秦始皇陵考古研究最關(guān)心的問題之一,利用4.1 節(jié)中得到的墓室墻壁的角坐標(biāo),結(jié)合兩個(gè)測(cè)量點(diǎn)的三維坐標(biāo)可以重建墓室的三維位置.如圖6(a)所示,由幾何關(guān)系可計(jì)算出墓室南北方向長(zhǎng)53.35 m,墓室中心埋深22 m.再由圖5(a)中墓室東墻、西墻的角坐標(biāo)結(jié)合墓室埋深可計(jì)算出墓室東西方向長(zhǎng)度為85.82 m,如圖6(b)所示.重建得到的墓室墻壁邊長(zhǎng)相較于理論值的差異約為7%,埋深差異約為6%.該重建結(jié)果相對(duì)于真值的差異與模擬中探測(cè)器位置、統(tǒng)計(jì)誤差、算法等有關(guān).在未來的實(shí)際測(cè)量中,還需要考慮探測(cè)器的角分辨率、幾何接受度、探測(cè)效率,以及本底和噪聲等對(duì)徑跡重建和圖像重建的影響.

圖6 墓室三維重建結(jié)果 (a) 剖面1 處重建結(jié)果;(b)剖面2 處重建結(jié)果Fig.6.Three-dimensional reconstruction results of the chamber:(a) Reconstruction result at Profile 1;(b) reconstruction result at Profile 2.

5 總結(jié)

本文使用GEANT4,基于現(xiàn)有的秦始皇陵考古數(shù)據(jù)以及理想μ子探測(cè)器,對(duì)秦始皇陵地宮μ子吸收成像進(jìn)行了模擬研究.成像結(jié)果初步驗(yàn)證了秦始皇陵地宮μ子吸收成像的可行性,單視角獲得的μ子通量投影數(shù)據(jù)可以清晰地反映出地宮內(nèi)部的不同結(jié)構(gòu),并能給出墓室墻壁的二維角坐標(biāo)信息,利用兩個(gè)視角下獲得的通量數(shù)據(jù)可以重建出墓室大小和三維位置,重建得到的墓室邊長(zhǎng)和墓室中心埋深相對(duì)于理論值的差異在7%左右.本文僅是對(duì)μ子吸收成像應(yīng)用于秦始皇陵地宮無損探測(cè)的初步研究,后續(xù)研究中將進(jìn)一步細(xì)化秦始皇陵模型和μ子探測(cè)器模型,調(diào)整測(cè)量點(diǎn)的位置,增加測(cè)量點(diǎn)數(shù)量,優(yōu)化多視角三維圖像重建算法,深入分析影響圖像重建結(jié)果的各種因素.