賈文寶陳奕澤黑大千趙 冬程 璨,3

1（南京航空航天大學(xué) 南京 210000）

2（蘭州大學(xué) 蘭州 730000）

3（核技術(shù)應(yīng)用教育部工程研究中心東華理工大學(xué) 南昌 330013）

瞬發(fā)γ射線成像（Prompt Gamma-ray Activation Image，PGAI）技術(shù)是一種新型元素成像方法，該技術(shù)基于瞬發(fā)γ射線中子活化分析（Prompt Gamma Neutron Activation Analysis，PGNAA）并結(jié)合準(zhǔn)直測(cè)量手段，可以實(shí)現(xiàn)對(duì)樣品內(nèi)部元素空間分布分析。由于中子具有良好的穿透能力，這一優(yōu)勢(shì)使該技術(shù)十分適用于大塊樣品的成分和結(jié)構(gòu)分析。近年來(lái)，該技術(shù)已被應(yīng)用于文物、天文和安全檢測(cè)等眾多領(lǐng)域的分析工作。PGAI技術(shù)的基礎(chǔ)是PGNAA技術(shù)，其利用中子與核素發(fā)生非彈性散射反應(yīng)或輻射俘獲反應(yīng)，在極短的時(shí)間內(nèi)產(chǎn)生能夠表征核素信息的特征γ射線，通過(guò)采集特征γ射線能譜并對(duì)樣品內(nèi)核素進(jìn)行分析。

由于PGAI技術(shù)在多個(gè)領(lǐng)域具有重要應(yīng)用前景，包括我國(guó)在內(nèi)的多個(gè)國(guó)家均開(kāi)展了大量的研究。目前PGAI測(cè)量平臺(tái)主要集中在反應(yīng)堆中子源或大型加速器中子源上搭建，目前較為成熟的測(cè)量平臺(tái)位于匈牙利布達(dá)佩斯的中子中心（Budapest Neutron Center，BNC）的反應(yīng)堆和德國(guó)慕尼黑工業(yè)大學(xué)FRMⅡ反應(yīng)堆上［1-4］。在這兩個(gè)PGAI測(cè)量平臺(tái)上已經(jīng)完成大量實(shí)驗(yàn)，包括對(duì)古文物檢測(cè)、核燃料位置識(shí)別、化石生物結(jié)構(gòu)分析等，樣品處的熱中子通量為107～109s-1·cm-2，可以實(shí)現(xiàn)較高分辨率（2～5 mm量級(jí)）的元素成像［5-7］。隨著新型的探測(cè)器和探測(cè)技術(shù)發(fā)展，美國(guó)國(guó)家標(biāo)準(zhǔn)與技術(shù)研究院（National Institute of Standards and Technology，NIST）反應(yīng)堆上開(kāi)展了基于康普頓相機(jī)的PGAI研究［8］，但是對(duì)γ射線的能量分辨率以及對(duì)高能γ射線的探測(cè)效率受到限制，仍有待探索。韓國(guó)原子能科學(xué)院也開(kāi)展了基于符合測(cè)量的PGAI研究，對(duì)相關(guān)探測(cè)器系統(tǒng)進(jìn)行了模擬探索研究［9］。國(guó)內(nèi)方面，在中國(guó)先進(jìn)研究堆（China Advanced Research Reactor，CARR）以及中國(guó)綿陽(yáng)研究堆（China Mianyang Research Reactor，CMRR）上也開(kāi)展了PGAI技術(shù)的研究工作，但目前處于驗(yàn)證性研究的階段［10］。

然而，基于反應(yīng)堆中子源的PGAI技術(shù)無(wú)法實(shí)現(xiàn)現(xiàn)場(chǎng)應(yīng)用，極大地限制了該技術(shù)的應(yīng)用場(chǎng)景?；谕凰刂凶釉春椭凶影l(fā)生器的PGAI技術(shù)在現(xiàn)場(chǎng)應(yīng)用中更具優(yōu)勢(shì)，在先前的研究中，利用DT中子發(fā)生器搭建了PGAI測(cè)量平臺(tái)，對(duì)B樣品、金屬樣品進(jìn)行了元素成像，空間分辨率可達(dá)到約4 cm［11-12］。但小型中子源存在中子通量低的缺點(diǎn)，為保證重建圖像質(zhì)量，要想在合理的測(cè)量時(shí)間內(nèi)進(jìn)行測(cè)量，則需要增大測(cè)量體積，使得測(cè)量結(jié)果的空間分辨率下降。

編碼成像技術(shù)為解決空間分辨率與測(cè)量時(shí)間之間的矛盾提供了重要的思路，由于編碼開(kāi)孔數(shù)量多，探測(cè)系統(tǒng)的探測(cè)效率高，同時(shí)單個(gè)編碼孔尺寸小，保證了重建圖像的分辨率。目前編碼成像技術(shù)主要應(yīng)用于編碼孔徑成像、時(shí)序編碼成像、熒光光譜分析（X-ray Fluorescence，XRF）等技術(shù)中［13-15］。在編碼孔徑成像中，往往采用陣列探測(cè)器對(duì)信號(hào)進(jìn)行收集，但陣列探測(cè)器無(wú)法實(shí)現(xiàn)特征峰的能量測(cè)量，不適合PGAI測(cè)量。而在時(shí)序編碼成像中利用時(shí)間調(diào)制信號(hào)強(qiáng)度實(shí)現(xiàn)單個(gè)探測(cè)器的位置甄別，代替了陣列探測(cè)器。在XRF成像技術(shù)中，研究者近年同樣將目光聚焦于能夠提供更高X射線熒光傳輸效率的編碼孔準(zhǔn)直器上，以彌補(bǔ)該技術(shù)中常采用的單孔準(zhǔn)直器傳輸效率低的缺陷。

本文通過(guò)MCNP，對(duì)基于多編碼板準(zhǔn)直的PGAI技術(shù)進(jìn)行模擬研究。以氯（Cl）作為目標(biāo)元素進(jìn)行二維元素成像，采用高純鍺探測(cè)器對(duì)經(jīng)過(guò)不同編碼板準(zhǔn)直器準(zhǔn)直的γ射線進(jìn)行收集，利用MLEM迭代算法對(duì)多組γ射線信號(hào)進(jìn)行分析，擬實(shí)現(xiàn)對(duì)1 cm×1 cm像素大小的元素空間分布圖像。

1 方法與模擬設(shè)計(jì)

1.1 多編碼準(zhǔn)直器PGAI測(cè)量方法

基于多編碼準(zhǔn)直器的PGAI利用中子束對(duì)整塊樣品進(jìn)行輻照，并用多塊編碼板準(zhǔn)直器對(duì)產(chǎn)生的γ射線進(jìn)行準(zhǔn)直，被準(zhǔn)直的γ射線被后方的高純鍺探測(cè)器接收，基本原理如圖1所示。采集到的特征γ射線計(jì)數(shù)為：

圖1 基于多編碼板準(zhǔn)直的PGAI技術(shù)原理示意圖Fig.1 Schematic diagram of PGAI measurement based on multi-coded collimators

式中：Q為探測(cè)器所采集的待測(cè)量核素產(chǎn)生的特征γ射線凈計(jì)數(shù)；NA為阿伏伽德羅常數(shù)；m為該核素的質(zhì)量；M為該核素的原子量；η為該核素的天然同位素豐度；φ為中子通量；σ為中子與該核素的反應(yīng)截面；λ為特征γ射線分支比；ε為探測(cè)效率；t為測(cè)量時(shí)間。

將單個(gè)編碼孔大小作為樣品的最小像素單元，探測(cè)器收集的特征γ射線計(jì)數(shù)由編碼準(zhǔn)直板開(kāi)孔位置對(duì)應(yīng)的所有像素所激發(fā)出來(lái)的特征γ射線乘上相應(yīng)探測(cè)效率累積而來(lái)。此時(shí)，式（1）可以表示為：

致病菌的篩選和培養(yǎng)步驟完成后，需要對(duì)其進(jìn)行生化實(shí)驗(yàn)鑒定。在鑒定過(guò)程中發(fā)現(xiàn)致病菌能發(fā)酵葡萄糖、乳糖、甘露醇、麥芽糖，且能產(chǎn)生酸氣，但是不能發(fā)酵蔗糖。而通過(guò)M.R試驗(yàn)、V.P試驗(yàn)和H2S試驗(yàn)后發(fā)現(xiàn)試驗(yàn)結(jié)果均為陰性，有利于接下來(lái)采取菌種DNA測(cè)序診斷措施。對(duì)菌種進(jìn)行DNA測(cè)序時(shí)，檢測(cè)人員需要取分離純化菌株接種于營(yíng)養(yǎng)瓊脂培養(yǎng)基上，然后利用相關(guān)儀器及技術(shù)進(jìn)行檢測(cè)，根據(jù)同源性和相似性分析最終得出所鑒定的菌株為大腸桿菌。

式中：φn表示像素n位置處的中子通量；εn表示像素n產(chǎn)生的特征γ射線被探測(cè)器采集的探測(cè)效率；mn為像素n處待測(cè)核素的質(zhì)量；K為式（1）中其他常數(shù)參數(shù)的乘積。

采用i塊不同開(kāi)孔方式的編碼準(zhǔn)直板進(jìn)行多次測(cè)量后，將得到如下方程：

式（3）中，[Q]可由采集的γ射線能譜獲取，[ε]中每個(gè)像素點(diǎn)對(duì)不同編碼板時(shí)，探測(cè)效率都是不同的，需要通過(guò)模擬獲取，而[φ]會(huì)受到樣品的中子自屏效應(yīng)影響，需要模擬獲取中子自屏修正因子。當(dāng)上述參數(shù)都獲得后，便可以通過(guò)最大似然期望最大化（Maximum Likelihood Expectation Maximization，MLEM）算法，計(jì)算出[m]，即重建元素分布圖像［16］。MLEM是建立在泊松模型上的統(tǒng)計(jì)迭代算法，中子與物質(zhì)發(fā)生輻射俘獲反應(yīng)或非彈性散射激發(fā)出特征γ射線的過(guò)程同樣服從泊松分布，同時(shí)統(tǒng)計(jì)迭代算法能較方便在圖像邊界、噪聲位置添加約束條件，獲得質(zhì)量較高的重建圖像［17］。迭代公式為：

式中：f k+1j表示像素j迭代了k+1次后的重建灰度值；f kj表示像素j迭代了k次后的重建灰度值；Hij表示第i塊編碼板準(zhǔn)直時(shí)，像素j的響應(yīng)函數(shù)，在本研究中表示中子通量與探測(cè)效率的乘積。

1.2 MCNP模型

本研究使用蒙特卡羅模擬軟件MCNP5對(duì)設(shè)定樣品中的Cl元素進(jìn)行檢測(cè)仿真，采用的核數(shù)據(jù)庫(kù)為ENDF/B-VI。利用MCNP5建立的幾何模型如圖2所示，測(cè)量裝置由中子源、樣品、編碼準(zhǔn)直器、γ射線探測(cè)器組成，其幾何中心均在y軸正方向上。

圖2 元素成像測(cè)量裝置MCNP模型Fig.2 MCNP model of element imaging measuring device

中子源設(shè)置為單向面源，截面積為6 cm×6 cm，以y軸正方向發(fā)射能量為0.025 eV的熱中子，粒子數(shù)設(shè)置為十億，以確保模擬計(jì)數(shù)的統(tǒng)計(jì)性。

樣品尺寸設(shè)置為6 cm×6 cm×1 cm（長(zhǎng)×寬×厚），像素點(diǎn)大小為1 cm×1 cm（長(zhǎng)×寬），每個(gè)像素點(diǎn)以一定的空間分布填充氯化鈉（NaCl）和碳酸鈉（Na2CO3），其中NaCl又設(shè)置了兩種不同的密度，樣品成分分布情況以及像素編號(hào)如圖3所示。

圖3 樣品成分分布情況Fig.3 Composition distribution of sample

準(zhǔn)直器選擇Fe作為材料，橫截面為8 cm×8 cm，其中編碼開(kāi)孔區(qū)域是樣品面積對(duì)應(yīng)的6 cm×6 cm范圍內(nèi)，同時(shí)向四周覆蓋有1 cm Fe以屏蔽從周圍漏射的γ射線，如圖2所示。準(zhǔn)直長(zhǎng)度（y軸方向）為10 cm。編碼方式為經(jīng)過(guò)挑選的隨機(jī)開(kāi)孔方式，開(kāi)孔率為50%。一共設(shè)計(jì)了36塊編碼板準(zhǔn)直器，設(shè)計(jì)的原則是盡量保證每個(gè)像素點(diǎn)對(duì)應(yīng)的開(kāi)孔次數(shù)在18次左右。圖4給出了設(shè)計(jì)的部分編碼圖案。

圖4 幾種典型的編碼方式Fig.4 Several typical pattern of coded-aperture collimators

γ射線探測(cè)器選擇能量分辨率高的高純鍺探測(cè)器，模型根據(jù)型號(hào)為GMX30P4-70的高純鍺探測(cè)器進(jìn)行構(gòu)建［18］，其鍺晶體直徑為66.8 mm，長(zhǎng)70.7 mm。

本模擬以熱中子與核素發(fā)生輻射俘獲反應(yīng)為基礎(chǔ)。樣品中Cl元素與熱中子反應(yīng)發(fā)射的能量在2 MeV以下的部分特征γ射線能量由表1給出［18］，最終選擇能量為1.164 MeV的特征γ射線對(duì)Cl元素分布圖像進(jìn)行重建。

表1 Cl的特征γ射線能量Table 1 Characteristic gamma ray energy of Cl

本研究需要對(duì)中子通量和探測(cè)效率進(jìn)行模擬。MCNP模擬軟件中，可以用F4計(jì)數(shù)卡對(duì)樣品柵元內(nèi)的平均中子通量進(jìn)行模擬，同時(shí)也可以用來(lái)對(duì)中子場(chǎng)進(jìn)行修正；可用F8計(jì)數(shù)卡模擬探測(cè)器中的能量沉積以及探測(cè)效率，同時(shí)通過(guò)與E8計(jì)數(shù)能量卡結(jié)合使用，可以模擬探測(cè)器能譜，從而獲得Cl元素特征γ射線計(jì)數(shù)。

2 結(jié)果和討論

2.1 γ射線能譜

用MCNP中F8計(jì)數(shù)卡記錄高純鍺探測(cè)器中Ge晶體內(nèi)的γ射線能量沉積情況，模擬獲取的γ射線能譜圖如圖5（a）所示，表1中部分能量的特征峰在圖5（a）中清晰可見(jiàn)。選擇能量為1.164 MeV的特征γ射線計(jì)數(shù)進(jìn)行Cl元素圖像重建，計(jì)算該能量的特征峰峰面積，在不同編碼板準(zhǔn)直器準(zhǔn)直后，采集到的Cl元素特征γ射線計(jì)數(shù)會(huì)產(chǎn)生變化，這種變化是隨編碼方式的變化而變化的，這也是能夠通過(guò)MLEM算法重建圖像的基礎(chǔ)，部分編碼準(zhǔn)直后，1.164 MeV特征峰如圖5（b）所示，圖中也展現(xiàn)了峰面積的差異。

圖5 γ射線能譜圖(a)，不同編碼準(zhǔn)直板準(zhǔn)直下的Cl元素1.164 MeV特征峰(b)Fig.5 γ ray spectrogram(a)and Cl element 1.164 MeV characteristic peak under different coded collimators(b)

2.2 探測(cè)效率模擬結(jié)果

探測(cè)效率與特征γ射線能量、探測(cè)器靈敏體積和形狀、探測(cè)器與源所張立體角等因素有關(guān)。中子活化樣品中Cl元素后，產(chǎn)生特征γ射線，此時(shí)可以將每一個(gè)像素點(diǎn)都看成一個(gè)能量為1.164 MeV的γ射線源，因此模擬獲取探測(cè)效率矩陣時(shí)，直接將該能量的γ射線源依次設(shè)置在每個(gè)像素點(diǎn)內(nèi)來(lái)獲取每個(gè)像素點(diǎn)對(duì)應(yīng)某一開(kāi)孔方式的編碼板準(zhǔn)直器的探測(cè)效率，并用F8計(jì)數(shù)卡記錄，這樣不僅可以加快模擬的速度，同時(shí)也不用再考慮γ射線自吸收的問(wèn)題。

對(duì)于編碼開(kāi)孔的位置，γ射線能透過(guò)準(zhǔn)直板到達(dá)探測(cè)器位置，而未開(kāi)孔的位置，γ射線被屏蔽，因此探測(cè)效率矩陣會(huì)與編碼開(kāi)孔方式所對(duì)應(yīng)，通過(guò)模擬也印證了這一點(diǎn)，但產(chǎn)生的特征γ射線是各向同性發(fā)射的，某一像素點(diǎn)四周的開(kāi)孔情況也會(huì)影響到探測(cè)效率，因此探測(cè)效率矩陣和開(kāi)孔方式有一定差異，如圖6所示，為圖4給出的編碼方式下獲得的探測(cè)效率灰度圖。

圖6 部分探測(cè)效率灰度圖Fig.6 Gray scale image of partial detection efficiency

2.3 Cl元素空間分布重建

模擬獲取多編碼板準(zhǔn)直后探測(cè)器采集的γ射線能譜圖，得到特征γ射線計(jì)數(shù)。對(duì)樣品中Cl元素的位置和相對(duì)含量進(jìn)行分析，采用最大似然期望最大化迭代算法進(jìn)行圖像重建，引入相對(duì)偏差df［19］以及結(jié)構(gòu)相似性（Structural Similarity，SSIM）［20］對(duì)重建圖像的質(zhì)量進(jìn)行評(píng)價(jià)。

相對(duì)偏差定義為：

式中：fj為重建圖像中像素j的灰度值；f0j為原圖像中像素j的灰度值。該參數(shù)能在已知原圖像的情況下，更加直觀地評(píng)價(jià)重建圖像的質(zhì)量。

結(jié)構(gòu)相似性定義為：

式中：μx為圖像x的均值；σx為圖像x的標(biāo)準(zhǔn)差；y同理；σxy為協(xié)方差；c1、c2為常數(shù)項(xiàng)，避免分母為0；SSIM的范圍為-1～1，兩個(gè)圖像越相似，這個(gè)值越接近1。該參數(shù)從亮度、對(duì)比度和結(jié)構(gòu)三個(gè)方面來(lái)對(duì)比兩張圖像。

對(duì)模擬樣品中Cl元素空間分布情況成像，結(jié)果如圖7所示。重建圖像與原圖像的相對(duì)偏差為0.065 8，SSIM為0.952 1，可以看出，對(duì)Cl元素位置分析十分準(zhǔn)確，密度相同的像素灰度值相近，但不同密度像素點(diǎn)的灰度值之間差異與實(shí)際質(zhì)量之間的差異不對(duì)應(yīng)，造成這一現(xiàn)象的原因是中子自屏效應(yīng)，因此利用MCNP對(duì)中子自屏效應(yīng)進(jìn)行修正。

圖7 未考慮中子自屏效應(yīng)的Cl元素重建圖像Fig.7 Cl element reconstruction without considering neutron self shielding effect

均勻中子面源發(fā)射的中子束在到達(dá)樣品表面時(shí)，中子通量是均勻的，但中子穿透樣品后，樣品內(nèi)部核素種類以及分布不均勻，會(huì)使得樣品內(nèi)部的中子場(chǎng)不均勻，這被稱為中子自屏效應(yīng)［21］。中子自屏效應(yīng)會(huì)影響核素的定量分析，因此需要通過(guò)一些手段來(lái)對(duì)其進(jìn)行修正。

在已知模擬樣品內(nèi)部的核素分布情況、密度的情況下，采用MCNP中的F4計(jì)數(shù)卡對(duì)該預(yù)設(shè)樣品各個(gè)像素點(diǎn)處的中子自屏修正因子進(jìn)行計(jì)算。將像素點(diǎn)柵元內(nèi)的F4計(jì)數(shù)與無(wú)限稀釋情況下該像素點(diǎn)柵元內(nèi)的F4計(jì)數(shù)之比作為該像素點(diǎn)的中子自屏修正因子fn［21］，可表示為：

由于NaCl密度不同，不同像素的中子自屏因子不同，密度大的像素中子自屏效應(yīng)越明顯，fn越小，反之亦然。對(duì)模擬Cl元素分布重建圖像進(jìn)行自屏修正后，如圖8（a）所示，對(duì)Cl元素的定量分析更加準(zhǔn)確，相對(duì)偏差減小為0.002 3，SSIM提高到0.998 4。含Cl元素的像素點(diǎn)的測(cè)量結(jié)果與真值的偏差均在10%以內(nèi)，如圖8（b）所示。在實(shí)際未知樣品的測(cè)量中，中子自屏修正是一個(gè)無(wú)法忽視的問(wèn)題，常用的中子自屏修正方法除了本文提到的模擬修正外，還可以結(jié)合熱中子照相技術(shù)或利用活化箔監(jiān)測(cè)等手段進(jìn)行實(shí)驗(yàn)修正。

圖8 中子自屏修正后的Cl元素重建圖像(a)，含Cl的11個(gè)像素模擬結(jié)果與真值的偏差(b)Fig.8 Reconstructed image of Cl element after neutron self screen correction(a),deviation between simulation results of 11 pixels with Cl and true value(b)

從模擬結(jié)果來(lái)看，將編碼孔應(yīng)用于PGAI技術(shù)是可行的，若要將該方法實(shí)施到實(shí)際實(shí)驗(yàn)中，需要注意以下幾個(gè)問(wèn)題：

1）對(duì)探測(cè)效率的刻度要準(zhǔn)確。由于每個(gè)像素對(duì)應(yīng)不同編碼準(zhǔn)直板時(shí)的探測(cè)效率都是不同的，通過(guò)實(shí)驗(yàn)一個(gè)一個(gè)測(cè)量并不合適，可以通過(guò)本文提到的用MCNP模擬獲取探測(cè)效率的方式，這要求已知探測(cè)器內(nèi)部結(jié)構(gòu)的尺寸、材料等信息。

2）中子自屏修正。針對(duì)實(shí)際樣品，采用合適的中子自屏修正方法，中子自屏效應(yīng)會(huì)影響到元素的定量分析。

3）可以通過(guò)同時(shí)添加中子束編碼準(zhǔn)直器，實(shí)現(xiàn)對(duì)樣品的三維成像。

3 結(jié)語(yǔ)

小型中子源PGAI技術(shù)具有很好的應(yīng)用前景，但受到小型中子源中子通量低的影響，成像分辨率較低。本文介紹了一種將編碼板準(zhǔn)直器結(jié)合到PGAI技術(shù)中用于元素成像的手段，旨在解決小型中子源成像分辨率低的問(wèn)題，并通過(guò)蒙特卡羅模擬手段驗(yàn)證了該方法的可行性。從模擬結(jié)果上來(lái)看，成功做到了1 cm×1 cm像素大小的元素空間分布成像，重建圖像與原圖像相似度高，相對(duì)偏差小。因此，從拓展PGAI技術(shù)應(yīng)用范圍的角度來(lái)看，基于編碼板準(zhǔn)直器的PGAI技術(shù)具有一定的研究?jī)r(jià)值，值得后續(xù)開(kāi)展相關(guān)實(shí)驗(yàn)。

作者貢獻(xiàn)聲明賈文寶：負(fù)責(zé)文章框架和最終版本的修訂；陳奕澤：負(fù)責(zé)模擬仿真及數(shù)據(jù)分析和文章起草撰寫(xiě)；黑大千：負(fù)責(zé)文章資料的查閱及整理；趙冬：負(fù)責(zé)文章資料的查閱；程璨：負(fù)責(zé)文章的修改和整體把握。