秦始皇陵地宮宇宙射線繆子吸收成像模擬研究*

蘇寧 劉圓圓? 王力 程建平

1) (北京師范大學核科學與技術學院,射線束技術教育部重點實驗室,北京 100875)

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

1 引言

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

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

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

2 μ子吸收成像算法

2.1 μ子與物質相互作用原理

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

圖1 (a) 實驗測量得到的不同方向上的海平面μ子通量[20];(b)探測器探測到的μ子的方位角φ 和天頂角θ,其中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.

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

其中E代表μ子能量,為μ子穿過物質的密度對穿透長度的積分,a對應電離導致的能量損失率,bE對應其他作用導致的能量損失率.a,b與μ子能量和所穿透物質的種類有關.

2.2 圖像重建算法

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

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

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

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

其中S為探測器面積,Ω為立體角,t為測量時間.

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

第三步,重構ROI 的三維圖像.探測器在單個測量點得到的數據僅能給出ROI 相對該測量點的二維角坐標信息,為了獲得ROI 的三維信息,需要有兩個或兩個以上不同位置的測量點的數據.以兩個測量點為例,當ROI 為簡單的幾何體(如長方體)時,可以很方便地將兩個測量點觀測到的ROI 邊界角坐標對應起來.如圖2,在測量點1 觀測到ROI 的面abcd的角坐標為 (θ1,φ1),在測量點2 觀測到面abcd的角坐標為 (θ1′,φ′1),則可結合兩個測量點的位置信息,根據幾何關系計算出面abcd的三維坐標,以此類推,最終重建出ROI 的幾何大小和三維位置信息.

圖2 測量點與ROI 之間的幾何關系示意圖Fig.2.Geometric relationship between viewpoint and ROI.

3 模擬系統建立

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

3.1 秦始皇陵地宮模型

本文根據使用綜合地球物理方法得到的秦始皇陵考古數據建立了兩個秦始皇陵地宮簡化模型[22].其中模型1 為待測物體模型,包括了封土堆、土地和地宮,其中地宮由細夯土墻、宮墻、回填夯土和墓室組成;模型2 為參考模型,僅包括封土堆和土地.模型幾何及關鍵結構的尺寸如圖3 所示.表1 列出了地宮模型中不同區域的材質和密度,其中空氣定義為70%的氮氣和30%的氧氣,黃土的組分定義參考了文獻[23].

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

3.2 μ子探測器模型

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

圖3 秦始皇陵模型示意圖 (a) 模型1 內部結構示意圖;(a1) 模型1 俯視圖;(a2) 模型1 正視圖;(a3) 模型1 剖面3 示意圖;(a4) 模型1 剖面1 示意圖;(b) 模型2 示意圖(無內部結構);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 μ子源產生

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

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

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

圖4 根據Reyna 公式抽樣產生的1000 萬個μ子的動量和天頂角分布 (a) μ子數量隨μ子動量變化分布;(b) μ子數量隨μ子速度方向的天頂角變化分布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 模擬結果

4.1 墓室二維成像結果

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

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

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

圖5 兩個測量點得到的 f (θ,φ) 的二維投影圖 (a)測量點1 的 f (θ,φ) 投影圖,其中,t anθx=tanθcosφ ,tanθy=tanθsinφ ;(b)測量點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 墓室三維重建結果

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

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

5 總結

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