龔春慧 曾國(guó)強(qiáng) 葛良全 湯曉斌 譚承君

1（成都理工大學(xué)核技術(shù)與自動(dòng)化工程學(xué)院 成都 610059）

2（南京航空航天大學(xué)核科學(xué)與工程系 南京 210016）

航空γ能譜測(cè)量在地質(zhì)找礦領(lǐng)域有著廣泛的應(yīng)用，它是將航空γ能譜儀系統(tǒng)安裝在飛機(jī)上，在飛機(jī)飛行過(guò)程中測(cè)量介質(zhì)放出的γ射線(xiàn)，通過(guò)分析獲得的γ能譜數(shù)據(jù)而獲得地表介質(zhì)鈾、釷、鉀含量[1]，因此對(duì)低能射線(xiàn)的探測(cè)并不關(guān)注。而在核事故中泄露放射性核素發(fā)出的γ射線(xiàn)能量普遍偏低，如最可能出現(xiàn)的137Cs的特征γ射線(xiàn)能量為661.6keV，131I為364.5keV。

使用航測(cè)系統(tǒng)用于核應(yīng)急監(jiān)測(cè)時(shí)，由于飛機(jī)外殼相當(dāng)于數(shù)厘米厚的鋁殼，對(duì)低能γ射線(xiàn)的衰減很?chē)?yán)重，采用探測(cè)器內(nèi)置的測(cè)量模式的探測(cè)效率很低。其次，由于空氣對(duì)低能γ射線(xiàn)的衰減也很?chē)?yán)重，為減少空氣層對(duì)γ射線(xiàn)的衰減，應(yīng)當(dāng)在盡可能低的高度進(jìn)行飛行測(cè)量，而固定翼飛機(jī)的最低飛行高度僅為約100m[2]。因此在使用該類(lèi)航測(cè)系統(tǒng)對(duì)環(huán)境中的低能射線(xiàn)進(jìn)行探測(cè)時(shí)，存在探測(cè)效率低等問(wèn)題，本文對(duì)GRS-16航測(cè)系統(tǒng)在地面輻射監(jiān)測(cè)中的應(yīng)用進(jìn)行模擬仿真，對(duì)以地面輻射為監(jiān)測(cè)目標(biāo)的航空γ能譜測(cè)量的方法研究做一些初步探索和研究。

1 探測(cè)器模型

Na(Tl)晶體密度大、探測(cè)效率高、發(fā)光效率好，且能通過(guò)鍛壓等技術(shù)制作體積較大的晶體，在航空 γ能譜測(cè)量領(lǐng)域得到廣泛應(yīng)用[3]。我國(guó)航空 γ能譜測(cè)量主要由石家莊核工業(yè)航測(cè)遙感中心承擔(dān)，使用的是 GRS-16航空多道 γ譜儀系統(tǒng)。GRS-16 由 15 條 10.16cm×10.16cm×40.64cm NaI(Tl)晶體組成，分裝成3箱，每箱包括五條晶體，成“品”字列排列，每條晶體外表包裹一層2mm的坡莫合金，箱體為厚度為2mm的鋁殼，如圖1。模型物質(zhì)成分如表1。

圖1 探測(cè)器幾何模型Fig.1 Detector geometric model.

表1 探測(cè)器模型的物質(zhì)成分Table 1 Ingredient of the detector model.

2 理論基礎(chǔ)

2.1 Monte Carlo方法

MCNP程序(Monte Carlo Neutron and Photo Transport Code)是由美國(guó)洛斯阿拉莫斯國(guó)家實(shí)驗(yàn)室的Monte Carlo小組（稱(chēng)為X-6組）在一系列程序工作基礎(chǔ)上，集中編制的一個(gè)大型多功能 Monte Carlo計(jì)算程序，可用于計(jì)算中子、光子和電子的聯(lián)合輸運(yùn)問(wèn)題以及臨界問(wèn)題，中子能量范圍從10?11–20MeV ， 光 子 和 電 子 的 能 量 范 圍 從1keV–1000MeV[4?5]。

在 MCNP輸入卡的描述中除使用了柵元、曲面、材質(zhì)和源的描述外，最關(guān)鍵的是聯(lián)合使用 F8和E8兩個(gè)計(jì)數(shù)卡。F8卡可以用于計(jì)算電子脈沖能量幅度分布，通過(guò)高斯展寬值設(shè)定模擬不同探測(cè)效率的探測(cè)器。E8卡稱(chēng)作記數(shù)能量卡，運(yùn)用它可劃分一組能量沉積箱[6]。在航測(cè)中由于有各種干擾粒子進(jìn)入探測(cè)器，為記錄從地面發(fā)射而進(jìn)入探測(cè)器的有效粒子就要?jiǎng)澏ㄒ粋€(gè)記數(shù)的有效能量范圍，E8就可以用于劃分一組有效能量沉積箱[7]。

MCNP程序中的重復(fù)結(jié)構(gòu)卡對(duì)于有相同結(jié)構(gòu)柵元的描述非常方便。重復(fù)結(jié)構(gòu)能力導(dǎo)致柵元概念的擴(kuò)充Universe。在一個(gè)U單元是一個(gè)柵元或一批柵元的集合，F(xiàn)ill卡能將U單元按指定的順序?qū)λ枋龅膸缀文Ｐ瓦M(jìn)行填充。TRCL卡使得只需一次掃描界定在形狀和尺寸是等同的只是在幾何位置上不同的幾個(gè)柵元的曲面。在本文模擬中，晶體分裝為三箱，且每箱內(nèi)部結(jié)構(gòu)相同，所以編寫(xiě)MCNP5輸入卡文件時(shí)，使用重復(fù)結(jié)構(gòu)卡以減少用戶(hù)必須提供的輸入數(shù)據(jù)量以及有大量重復(fù)結(jié)構(gòu)問(wèn)題所需的機(jī)器內(nèi)存。

2.2 等效質(zhì)量厚度的原理

等效質(zhì)量厚度是基于虛擬源原理：通過(guò)改變空氣密度來(lái)模擬不同的探測(cè)高度，對(duì)于輻射源經(jīng)過(guò)厚度為h1密度為ρ1的介質(zhì)1（航測(cè)系統(tǒng)空中測(cè)量）和厚度為h2密度為ρ2的介質(zhì)2（模擬實(shí)驗(yàn)）衰減后，在空中一點(diǎn)的未散射光子通量相同，對(duì)航測(cè)系統(tǒng)空中測(cè)量與模擬實(shí)驗(yàn)這兩種情況，到達(dá)飛機(jī)外殼任何一點(diǎn)的未散射光子通量相同，即作用于航測(cè)系統(tǒng)的輻射場(chǎng)相同，則航測(cè)系統(tǒng)的響應(yīng)也就相同，如圖2所示。

圖2 地面輻射監(jiān)測(cè)Fig.2 Ground radiation monitoring.

則：

同理，

由式(2)、(3)和(6)，有Φ1=Φ2，即按照式(4)選擇相同的質(zhì)量厚度，對(duì)同樣的輻射源來(lái)說(shuō)，同一能量γ射線(xiàn)經(jīng)過(guò)厚度為h1密度為ρ1的介質(zhì)1散射后到達(dá)飛機(jī)底部相同一點(diǎn)的未散射光子注量，與經(jīng)過(guò)厚度為h2密度為ρ2的介質(zhì)2散射后的相同，所以采用等效質(zhì)量厚度方法是可行的。

2.3 最小可探測(cè)活度MDA

最小可探測(cè)活度(Minimum Detectable Activity,MDA)是探測(cè)系統(tǒng)的重要性能指標(biāo)之一，如果某核素的活度大于其MDA值，則其可以被探測(cè)系統(tǒng)較可靠地測(cè)量到。航空γ能譜測(cè)量系統(tǒng)MDA[8?9]可用下式表示：

式中，TL為測(cè)量活時(shí)間，航測(cè)系統(tǒng)的探測(cè)周期一般為1s或數(shù)秒，TL一般為1s；p為γ射線(xiàn)的發(fā)射概率；εm為有效探測(cè)效率；LD為最小可探測(cè)限。

式中，ε為特征能量的光峰效率；S為源面積。

由于航測(cè)是在空中測(cè)量，探測(cè)范圍大，在模擬過(guò)程中，地面以有限面源代替無(wú)限面源。探測(cè)器與源的最大距離可以通過(guò)衰減系數(shù)來(lái)計(jì)算，取u(E)為能量為E的γ射線(xiàn)在空氣中的線(xiàn)性衰減系數(shù)，此處E為1.33MeV，從而確定源面積：

如式(10)所示，取航測(cè)系統(tǒng)與源的斜距為1km，足以用于模擬 γ射線(xiàn)能量小于1.33MeV的真實(shí)環(huán)境。源面積為：

式中，k為單邊置信度因子；σ0為分布的標(biāo)準(zhǔn)偏差(σ02=2B)；B為峰區(qū)本底。

式中，B1和B2分別為能譜感興趣區(qū)域(Region of interest, ROI)左邊m道和右邊m道的總計(jì)數(shù)；N為ROI的計(jì)數(shù)之和。本文中，m值取為1。

2.4 全能峰高斯展寬系數(shù)

采用實(shí)際航測(cè)的 γ能譜，譜儀的半高寬(Full Width at Half Maximum, FWHM)隨γ射線(xiàn)能量的關(guān)系如表2。

表2 FWHM隨γ射線(xiàn)能量的關(guān)系Table 2 Relationship between the FWHM and γ-ray energy.

半高寬FWHM與能量E的關(guān)系可用下式表示：

式中，a、b、c為系數(shù)，F(xiàn)WHM、E、a的單位是MeV，b的單位是的單位為 MeV?1。在MCNP輸入文件中用式(14)來(lái)保證輸出能譜的峰展寬嚴(yán)格和實(shí)驗(yàn)譜一致，即模擬能量沉積譜在探測(cè)器里面的分辨函數(shù)和實(shí)際探測(cè)器的分布函數(shù)相一致。

根據(jù)式(14)和表2中能量及對(duì)應(yīng)半高寬的數(shù)據(jù)，在Matlab中擬合得到MCNP中展寬a、b、c的值如下：a= 0.0282，b= 0.0135，c=?0.3186。用該值作為GEB卡輸入實(shí)現(xiàn)高斯能量展寬。

3 結(jié)果分析

在實(shí)際航測(cè)過(guò)程中，航測(cè)譜儀系統(tǒng)通常安裝在飛機(jī)上，由于飛機(jī)外殼對(duì)低能γ射線(xiàn)衰減很?chē)?yán)重，為了對(duì)探測(cè)器內(nèi)置和探測(cè)器外掛進(jìn)行比較，本文分別對(duì)兩種測(cè)量模式進(jìn)行模擬。

航測(cè)測(cè)量高度主要在 80–120m 范圍內(nèi)[10]，根據(jù)§2.2提出的等效質(zhì)量厚度方法，在模擬中將空氣密度擴(kuò)大 10倍，則飛機(jī)與源距離相應(yīng)縮短為8–12m。同時(shí)，該模型假設(shè)源粒子發(fā)射方向各向同性，為了避免計(jì)算效率太低，本文采用二次源方法進(jìn)行減小方差模擬計(jì)算，具體實(shí)施方法為：在探測(cè)器下方1cm處設(shè)置垂直距離為8–12m的平面記錄初始源的放射性粒子進(jìn)入該平面的粒子數(shù)及能譜，將此作為二次源模擬計(jì)算探測(cè)器的脈沖幅度分布。根據(jù)式(11)，可求得當(dāng)高度為 100m時(shí)源半徑設(shè)置為965m。

3.1 探測(cè)器外置

為了簡(jiǎn)化計(jì)算過(guò)程，加速計(jì)算結(jié)果的收斂速度，使用Monte Carlo模擬軟件MCNP5建立了一個(gè)相對(duì)簡(jiǎn)化的三維蒙特卡羅模型，如圖2，其中NaI(Tl)探測(cè)器幾何模型如圖1。圖3分別為MCNP5模擬探測(cè)器外置測(cè)量模式下，航測(cè)系統(tǒng)在不同高度時(shí)對(duì)137Cs和131I探測(cè)所得到的結(jié)果。圖3(a)中能量為0.662MeV的峰是137Cs的全能峰，在0.112MeV附近有一個(gè)峰，該峰相當(dāng)于一次散射射線(xiàn)。模擬中未考慮核素特征 X射線(xiàn)的發(fā)射。圖3(b)中能量為0.364MeV的峰是131I的全能峰，0.104MeV所在的峰為散射峰。

圖3 探測(cè)器外置時(shí)不同高度下137Cs (a)和131I (b)探測(cè)結(jié)果Fig.3 Pulse height distributions of external detectors for137Cs (a) and 131I (b) at different heights.

以上模擬出來(lái)的結(jié)果乘以源活度即可轉(zhuǎn)換為計(jì)數(shù)率，鑒于固定翼飛機(jī)的最低飛行高度僅為約100m，分別取137Cs與131I的沉積活度為37kBq·m?2、80kBq·m?2，在 100m 高度對(duì)其進(jìn)行探測(cè)，得到的計(jì)數(shù)率譜線(xiàn)如圖4所示。其中，航測(cè)譜儀距離地面100m處時(shí)，源面積：

對(duì)沉積活度為 37 kBq·m?2的137Cs，放射性活度為：

對(duì)沉積活度為 80 kBq·m?2的131I，放射性活度：

對(duì)137Cs而言，

對(duì)131I而言，

3.2 探測(cè)器內(nèi)置

編寫(xiě)MCNP5輸入卡文件時(shí)，假設(shè)飛機(jī)機(jī)艙為10cm厚的7075型鋁合金，密度為0.281g·cm?3，元素含量百分比分別為 0.4%Si、0.5%Fe、1.8%Cu、0.3%Mn、2.5%Mg、0.23%Cr、5.6%Zn、0.2%Ti、88.47%Al。圖5為MCNP5模擬不同高度下航測(cè)系統(tǒng)對(duì)137Cs和131I的探測(cè)得到的結(jié)果。同樣地，圖5(b)中能量為0.364MeV的峰是131I的全能峰。圖5(a)、(b)均可看出，探測(cè)效率隨著航測(cè)高度的增加在減小。

圖5 探測(cè)器內(nèi)置時(shí)不同高度下137Cs (a)和131I (b)探測(cè)結(jié)果Fig.5 Pulse height distributions of internal detectors for137Cs (a) and 131I (b) at different heights.

同樣地，以上模擬出來(lái)的結(jié)果乘以源活度即可轉(zhuǎn)換為計(jì)數(shù)率，鑒于固定翼飛機(jī)的最低飛行高度僅為約100m，分別在100m高度對(duì)137Cs和131I進(jìn)行探測(cè)，分別取137Cs與131I的沉積活度為37kBq·m?2、80kBq·m?2，得到的計(jì)數(shù)率譜線(xiàn)如圖6所示。

圖6 探測(cè)器內(nèi)置100m處對(duì)137Cs、131I計(jì)數(shù)率譜線(xiàn)Fig.6 Spectra of internal detector in the present works for137Cs and 131I at the height of 100m.

對(duì)137Cs而言，

綜上所述，使用蒙特卡羅方法計(jì)算探測(cè)器外置與探測(cè)器內(nèi)置兩種測(cè)量模式下，航測(cè)系統(tǒng)對(duì)137Cs和131I的最小可探測(cè)活度如表3所示。

表3 航測(cè)系統(tǒng)對(duì)137Cs和131I的MDA模擬計(jì)算Table 3 MDAs for 137Cs and 131I in simulation, ascalculated from their most intense γ-rays, for the external detectors and internal detectors.

從圖3(a)和 5(a)可以看到，探測(cè)器外置時(shí)，在高度為80m、90m、100m、110m、120m全能峰峰位處對(duì)應(yīng)的相對(duì)計(jì)數(shù)率分別為1.097×10?8cm?2、9.563×10?9cm?2、8.279×10?9cm?2、7.233×10?9cm?2、6.346×10?9cm?2，相應(yīng)地，探測(cè)器內(nèi)置時(shí)對(duì)應(yīng)的相對(duì)計(jì)數(shù)率分別為 9.478×10?9cm?2、8.263×10?9cm?2、7.164×10?9cm?2、6.250×10?9cm?2、5.485×10?9cm?2；從圖3(b)和5(b)可以看到，探測(cè)器外置時(shí)，在高度為80m、90m、100m、110m、120m對(duì)應(yīng)的相對(duì)計(jì)數(shù)率分別為 1.265×10?8cm?2、1.079×10?8cm?2、9.134×10?9cm?2、7.778×10?9cm?2、6.699×10?9cm?2，相應(yīng)地，探測(cè)器內(nèi)置時(shí)對(duì)應(yīng)的相對(duì)計(jì)數(shù)率分別為1.054×10?8cm?2、8.959×10?9cm?2、7.608×10?9cm?2、6.497×10?9cm?2、5.598×10?9cm?2。綜上可知，在不同的探測(cè)高度，相同能量對(duì)應(yīng)的計(jì)數(shù)隨高度的增加而減小，說(shuō)明探測(cè)器距離源越遠(yuǎn)，其探測(cè)效率越低。同時(shí)，結(jié)合圖4和圖6，在100m處探測(cè)器外置、探測(cè)器內(nèi)置時(shí)，對(duì)137Cs探測(cè)時(shí)其全能峰峰位處對(duì)應(yīng)的相對(duì)計(jì)數(shù)率分別為 953.941s?1、821.788 s?1，對(duì)131I探測(cè)時(shí)其全能峰峰位處對(duì)應(yīng)的相對(duì)計(jì)數(shù)率分別為 2275.643 s?1、1894.897 s?1，可知相同的能量對(duì)應(yīng)的計(jì)數(shù)探測(cè)器外置比探測(cè)器內(nèi)置大。同時(shí)，根據(jù)表3的探測(cè)器外置和探測(cè)器內(nèi)置MDA值對(duì)比，可以看出在相同探測(cè)條件下，探測(cè)器外置優(yōu)于探測(cè)器內(nèi)置。這是由于飛機(jī)外殼對(duì)γ射線(xiàn)有衰減，模擬結(jié)果與理論相符。為了避免飛機(jī)外殼對(duì)γ射線(xiàn)的衰減，實(shí)際航測(cè)時(shí)，建議使用探測(cè)器外置方式進(jìn)行測(cè)量。

圖7為該航測(cè)系統(tǒng)空中飛行時(shí)在天然本底地區(qū)的測(cè)量能譜，計(jì)算其在137Cs和131I特征γ射線(xiàn)能區(qū)的本底，然后根據(jù)式(7)計(jì)算其對(duì)137Cs和131I的最小可探測(cè)活度分別為 244 Bq·m?2和 386 Bq·m?2（測(cè)量周期為1s，飛行高度為100m）。

圖7 空中飛行時(shí)實(shí)測(cè)譜線(xiàn)Fig.7 Measured spectrum in airborne detection.

由式(7)可知，本底越小，探測(cè)限越低，MDA越小。圖7所示這條譜線(xiàn)為天然本底地區(qū)的測(cè)量能譜，其放射性核素137Cs和131I的含量均很低。

4 結(jié)語(yǔ)

通過(guò)Monte Carlo程序模擬計(jì)算，建立了航空測(cè)量能譜的模擬方法，通過(guò)實(shí)驗(yàn)結(jié)果計(jì)算并得到全能峰高斯展寬系數(shù)a、b、c；同時(shí)得到探測(cè)器外置以及探測(cè)器內(nèi)置時(shí)不同高度對(duì)137Cs和131I的探測(cè)結(jié)果，并對(duì)MDA進(jìn)行計(jì)算。在相同測(cè)量條件下，探測(cè)器外置測(cè)量模式優(yōu)于探測(cè)器內(nèi)置測(cè)量模式。

實(shí)際航測(cè)時(shí)，在完全相同情況下，采用鋁合金包裹NaI晶體對(duì)于137Cs的探測(cè)效率相對(duì)不銹鋼提高 8.7%，鋁合金的原子序數(shù)以及密度都低于不銹鋼。本文模擬坡莫合金包裹NaI晶體，考慮到探測(cè)效率對(duì)MDA的影響，可使用碳纖維包裹NaI晶體以提高探測(cè)效率。同時(shí)模擬仿真的條件相對(duì)實(shí)際情況比較理想，譬如對(duì)飛機(jī)的模擬僅采用10cm厚的7075型鋁合金做簡(jiǎn)單替代，考慮的問(wèn)題也相對(duì)簡(jiǎn)單，在實(shí)際情況下，航測(cè)譜儀在進(jìn)行測(cè)量時(shí)，在測(cè)量地面放射性的同時(shí)，空氣中放射性也有涉及。在下一步工作中，同時(shí)對(duì)地面輻射、空氣輻射進(jìn)行模擬仿真研究，并對(duì)比分析在不同探測(cè)手段下的MDA，探索不同放射性環(huán)境下最為合適的用于航空環(huán)境放射性調(diào)查的探測(cè)器及其安裝方法，為下一階段航空放射性測(cè)量?jī)x器的研發(fā)做一定的指導(dǎo)工作。

1 谷懿. 航空伽瑪能譜測(cè)量大氣氡校正方法研究[D]. 成都理工大學(xué), 2010 GU Yi. Study on radon influence estimation in airborne gamma-ray spectrometry[D]. Chengdu University of Technology, 2010

2 李惠彬, 李君利. 航空放射性測(cè)量技術(shù)對(duì)核事故泄露核素的探測(cè)能力評(píng)估[J]. 核電子學(xué)與探測(cè)技術(shù), 2013,33(5): 531–534 LI Huibin, LI Junli. Evaluation of airborne radiation survey used for monitoring leakage of radionuclides in nuclear accident[J]. Nuclear Electronics and Detection Technology, 2013, 33(5): 531–534

3 林帆, 張慶賢, 葛良全, 等. 多晶體 γ探測(cè)器相互影響的蒙特卡洛數(shù)值模擬[J]. 核技術(shù), 2010, 33(9): 675–678 LIN Fan, ZHANG Qingxian, GE Liangquan,et al. A Geant4-based M-C simulation on the interaction of multi-crystal gamma spectrometer[J]. Nuclear Techniques,2010, 33(9): 675–678

4 X-5 Monte Carlo Team. MCNP-A General Monte Carlo N-Particle Transport Code[R]. Version 5, Los Alamos National Laboratory, Report LA-UR-03-1987, 2003

5 許淑艷. 蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用（修訂版）[M]. 北京: 原子能出版社, 2006 XU Shuyan. Monte Carlo methods in experimental nuclear physics (revised edition)[M]. Beijing: Published by Atomic Energy Press, 2006

6 GRS-16/GRS-10/GRS-2 Gamma Spectrometer Operation Manual[R]. Calibration Program PEICalib, Version 5.0,2006: 4–5

7 孫健, 來(lái)永芳, 張力軍, 等. 無(wú)人直升機(jī)航測(cè)系統(tǒng)探測(cè)效率的M-C計(jì)算[C]. 第十四屆全國(guó)核電子學(xué)與核探測(cè)技術(shù)學(xué)術(shù)年會(huì)論文集(下冊(cè)), 2008: 877–811 SUN Jian, LAI Yongfang, ZHANG Lijun,et al.Calculation of detecting efficiency for unmanned-airplane airborne radioactivity measuring system by Monte Carlo[C]. 14thsession of the National Nuclear Science and Nuclear Detection Technology Annual Meeting Proceedings(b), 2008: 877–811

8 師全林, 常永福, 張佳媚. 最小可探測(cè)活度 MDA的統(tǒng)計(jì)學(xué)解釋[C]. 全國(guó)第四屆核儀器及其應(yīng)用學(xué)術(shù)會(huì)議論文集, 2003: 188–195 SHI Quanlin, CHANG Yongfu, ZHANG Jiamei. A statistical interpretation to MDA[C]. Proceedings of the 4thNational Conference on Nuclear Instrument & Its Application, 2003: 188–195

9 Bagatelas C, Tsabaris C, Kokkoris M,et al. Determination of marine gamma activity and study of the minimum detectable activity (MDA) in 4pi geometry based on Monte Carlo simulation[J]. Environmental Monitoring and Assessment, 2010, 165: 159–168

10 賈明雁, 黃雄亮, 蘇旭軍, 等. 基于虛擬源原理的航空放射性測(cè)量系統(tǒng)刻度裝置[P]. CN 201732166 U, 2011 JIA Mingyan, HUANG Xiongliang, SU Xujun,et al.Scale device of aerial radioactivity measuring system based on virtual source principle[P]. CN 201732166 U,2011