多晶體航空γ能譜儀的地面輻射監(jiān)測模擬仿真

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

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

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

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

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

1 探測器模型

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

圖1 探測器幾何模型Fig.1 Detector geometric model.

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

2 理論基礎(chǔ)

2.1 Monte Carlo方法

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

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

MCNP程序中的重復(fù)結(jié)構(gòu)卡對于有相同結(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)相同，所以編寫MCNP5輸入卡文件時(shí)，使用重復(fù)結(jié)構(gòu)卡以減少用戶必須提供的輸入數(shù)據(jù)量以及有大量重復(fù)結(jié)構(gòu)問題所需的機(jī)器內(nèi)存。

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

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

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

則：

同理，

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

2.3 最小可探測活度MDA

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

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

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

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

如式(10)所示，取航測系統(tǒng)與源的斜距為1km，足以用于模擬 γ射線能量小于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í)際航測的 γ能譜，譜儀的半高寬(Full Width at Half Maximum, FWHM)隨γ射線能量的關(guān)系如表2。

表2 FWHM隨γ射線能量的關(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)來保證輸出能譜的峰展寬嚴(yán)格和實(shí)驗(yàn)譜一致，即模擬能量沉積譜在探測器里面的分辨函數(shù)和實(shí)際探測器的分布函數(shù)相一致。

根據(jù)式(14)和表2中能量及對應(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í)際航測過程中，航測譜儀系統(tǒng)通常安裝在飛機(jī)上，由于飛機(jī)外殼對低能γ射線衰減很嚴(yán)重，為了對探測器內(nèi)置和探測器外掛進(jìn)行比較，本文分別對兩種測量模式進(jìn)行模擬。

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

3.1 探測器外置

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

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

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

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

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

對137Cs而言，

對131I而言，

3.2 探測器內(nèi)置

編寫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模擬不同高度下航測系統(tǒng)對137Cs和131I的探測得到的結(jié)果。同樣地，圖5(b)中能量為0.364MeV的峰是131I的全能峰。圖5(a)、(b)均可看出，探測效率隨著航測高度的增加在減小。

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

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

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

對137Cs而言，

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

表3 航測系統(tǒng)對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)可以看到，探測器外置時(shí)，在高度為80m、90m、100m、110m、120m全能峰峰位處對應(yīng)的相對計(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)地，探測器內(nèi)置時(shí)對應(yīng)的相對計(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)可以看到，探測器外置時(shí)，在高度為80m、90m、100m、110m、120m對應(yīng)的相對計(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)地，探測器內(nèi)置時(shí)對應(yīng)的相對計(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。綜上可知，在不同的探測高度，相同能量對應(yīng)的計(jì)數(shù)隨高度的增加而減小，說明探測器距離源越遠(yuǎn)，其探測效率越低。同時(shí)，結(jié)合圖4和圖6，在100m處探測器外置、探測器內(nèi)置時(shí)，對137Cs探測時(shí)其全能峰峰位處對應(yīng)的相對計(jì)數(shù)率分別為 953.941s?1、821.788 s?1，對131I探測時(shí)其全能峰峰位處對應(yīng)的相對計(jì)數(shù)率分別為 2275.643 s?1、1894.897 s?1，可知相同的能量對應(yīng)的計(jì)數(shù)探測器外置比探測器內(nèi)置大。同時(shí)，根據(jù)表3的探測器外置和探測器內(nèi)置MDA值對比，可以看出在相同探測條件下，探測器外置優(yōu)于探測器內(nèi)置。這是由于飛機(jī)外殼對γ射線有衰減，模擬結(jié)果與理論相符。為了避免飛機(jī)外殼對γ射線的衰減，實(shí)際航測時(shí)，建議使用探測器外置方式進(jìn)行測量。

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

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

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

4 結(jié)語

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

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

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

2 李惠彬, 李君利. 航空放射性測量技術(shù)對核事故泄露核素的探測能力評(píng)估[J]. 核電子學(xué)與探測技術(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 林帆, 張慶賢, 葛良全, 等. 多晶體 γ探測器相互影響的蒙特卡洛數(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 孫健, 來永芳, 張力軍, 等. 無人直升機(jī)航測系統(tǒng)探測效率的M-C計(jì)算[C]. 第十四屆全國核電子學(xué)與核探測技術(shù)學(xué)術(shù)年會(huì)論文集(下冊), 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 師全林, 常永福, 張佳媚. 最小可探測活度 MDA的統(tǒng)計(jì)學(xué)解釋[C]. 全國第四屆核儀器及其應(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 賈明雁, 黃雄亮, 蘇旭軍, 等. 基于虛擬源原理的航空放射性測量系統(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