龍澤宇，劉永闊，楊立群，陳志濤

哈爾濱工程大學(xué)核安全與仿真技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江哈爾濱150001

隨著核技術(shù)應(yīng)用的日益廣泛，各類放射性源的安全問題顯得尤為重要。由于放射源放出的射線人眼無法觀測到，因此工作人員無法直接確定放射源的位置，放射源一旦“失控”，可能會(huì)嚴(yán)重威脅人們的身體健康甚至生命安全。如何快速準(zhǔn)確地對放射源進(jìn)行定位是核安全領(lǐng)域的重要研究課題。從查閱的國內(nèi)外相關(guān)文獻(xiàn)來看，針對放射源定位的研究通常分為2類：一是放射源搜尋定位，其中應(yīng)用比較多的是基于方向探測器的搜尋定位方法，國內(nèi)外研究人員根據(jù)放射源的輻射特性設(shè)計(jì)出各種各樣的探測器，由人工手持或是車載，根據(jù)探測器給出的放射源方向信息，在移動(dòng)過程中不斷靠近放射源的位置[1-2]，應(yīng)用較少的則是γ 相機(jī)這類對輻射區(qū)域直接進(jìn)行“熱點(diǎn)”成像從而進(jìn)行尋源的方法[3]；二是放射源估計(jì)定位，與搜尋定位不同，此類方法并不依靠儀器對放射源進(jìn)行尋找，而是利用輻射環(huán)境中部分探測點(diǎn)的劑量信息對放射源的位置直接進(jìn)行大致估計(jì)[4-5]。然而，目前的這2類方法又都有著各自所面臨的困境。第一類方法雖然定位精度高且定位結(jié)果可靠，但依靠人工直接在輻射環(huán)境下進(jìn)行搜尋工作有悖于輻射防護(hù)中合理可行盡量低的原則(as low as reasonable practical,ALARP)，若輻射場強(qiáng)度較高且搜尋時(shí)間過長，則可能會(huì)給工作人員的身體帶來不可逆的損傷，而車載類探測器和γ相機(jī)則受制于高昂的成本；第二類方法的主要問題則是所需要的探測點(diǎn)數(shù)量通常非常多，或是無法應(yīng)用于有障礙物遮擋的情形下，這使得此類方法的適用性較差。

針對目前放射源估計(jì)定位中存在的問題，本文提出了一種基于點(diǎn)核積分的放射源定位方法，該方法可以直接利用輻射環(huán)境中探測點(diǎn)處的劑量率值對放射源進(jìn)行定位分析，需要的探測點(diǎn)數(shù)目較少，且由于點(diǎn)核積分方法的引入，解決了環(huán)境中障礙物遮擋的問題，極大地提高了方法的適用性。同時(shí)，本文利用蒙特卡羅程序MCNP5設(shè)計(jì)了虛擬仿真實(shí)驗(yàn)，驗(yàn)證了該方法的有效性。

1 點(diǎn)核積分

目前在輻射防護(hù)領(lǐng)域中，用于輻射劑量計(jì)算的方法主要包括蒙特卡羅方法和點(diǎn)核積分方法2種[6]。蒙特卡羅方法又稱隨機(jī)抽樣方法，通過對大量粒子進(jìn)行輸運(yùn)計(jì)算來模擬真實(shí)的物理過程，使得計(jì)算結(jié)果通常都能很好地符合實(shí)際情況，因此被廣泛應(yīng)用于對輻射劑量精確度要求較高的情況下。點(diǎn)核積分方法則是一種半經(jīng)驗(yàn)性的計(jì)算方法，通過引入累積因子來考慮散射光子對輻射劑量的貢獻(xiàn)，與蒙特卡羅方法相比，雖然計(jì)算精確度有所下降，但計(jì)算效率得到了極大提升，這也是能夠?qū)Ψ派湓催M(jìn)行快速定位的關(guān)鍵原因。

1.1 計(jì)算公式

點(diǎn)核積分的核心思想是將放射源按照幾何形狀離散為若干各向同性點(diǎn)源的集合，單個(gè)點(diǎn)源在探測點(diǎn)處的劑量率為[7]

式中：En為放射源的光子能量，MeV；rp與rd則分別為點(diǎn)源與探測點(diǎn)的位置；C(En)為光子通量到劑量率之間的轉(zhuǎn)換系數(shù)，(mSv·h-1)/(1·cm-2·s)；B為累積因子；S(En)為En能光子對應(yīng)的源強(qiáng)；t(En)為平均自由程數(shù)，計(jì)算公式為

式中：i為點(diǎn)源和探測點(diǎn)之間的屏蔽層序號(hào)；ui(En)為En能光子在第i層屏蔽中的線減弱系數(shù)，1/m；di則為光子在第i層屏蔽中的穿行距離，m。將式(1)按放射源空間體積進(jìn)行積分，即可得到整個(gè)放射源在探測點(diǎn)處的劑量率：

1.2 累積因子

探測器探測到的γ光子包括放射源發(fā)射并直接到達(dá)探測點(diǎn)處的直穿光子，也包括出射光子與路徑上的屏蔽物發(fā)生相互作用后產(chǎn)生的散射光子。因此，探測點(diǎn)處的劑量由計(jì)算簡單的直穿光子項(xiàng)和計(jì)算困難的散射光子項(xiàng)這2部分組成[8]。點(diǎn)核積分通過引入累積因子B來對散射光子項(xiàng)進(jìn)行考慮，其物理意義為在考察點(diǎn)r處，某一輻射量的總值與直穿光子造成的同一輻射量的比值：

式中：x為光子平均自由程數(shù)；a、b、c、d、xk為不同光子能量、不同材料下的參數(shù)。

2 放射源定位

2.1 問題描述

在實(shí)際情況中，放射源都有著一定的空間體積，在進(jìn)行計(jì)算的時(shí)候應(yīng)考慮離散。但當(dāng)探測點(diǎn)和放射源之間的距離是源最大線度的10倍以上時(shí)，就完全可以將放射源視為點(diǎn)源處理。即使只有5～7倍，劑量計(jì)算的誤差也不會(huì)高于5%[10]。在放射源定位的問題上，放射源體積與整個(gè)輻射區(qū)域相比通?？梢院雎?，因此在本研究中，將所有的放射源均視為點(diǎn)源來處理。整個(gè)問題可由圖1所示的平面示意圖進(jìn)行簡單描述。

圖1 放射源定位問題示意

其中探測點(diǎn)可由參數(shù)向量D n=[xn yn dn](n=1,2,…)進(jìn)行表示，其中(xn,yn)為第n號(hào)探測點(diǎn)的平面坐標(biāo)，dn為該探測點(diǎn)處的劑量率實(shí)測值；放射源可由S m=[xmymEm](m=1,2,…)進(jìn)行表示，其中(xm,ym)為第m號(hào)放射源的平面坐標(biāo)，Em為該放射源的γ光子能量。

2.2 區(qū)域離散

從理論上來說，在進(jìn)行放射源定位之前，整個(gè)輻射區(qū)域內(nèi)任意點(diǎn)都有可能是該定位問題的解，因此為了減小解集的規(guī)模，同時(shí)也為了便于對放射源的定位結(jié)果進(jìn)行描述，本文將輻射區(qū)域進(jìn)行柵格化處理[11]。使用邊長為l的正方形柵格來劃分輻射區(qū)域，如圖2所示，以最終解得的柵格的中心坐標(biāo)來表示放射源的定位結(jié)果。

圖2 輻射區(qū)域柵格化示意

對于離散化的輻射區(qū)域，如果某柵格被障礙物占據(jù)，則可以直接認(rèn)為該柵格為非解，在計(jì)算過程中可以直接跳過。另外，柵格的邊長l也會(huì)對定位結(jié)果產(chǎn)生很大的影響，如果柵格過大，柵格的密度就會(huì)很小，最后的定位精確度就會(huì)變差；如果柵格過小，柵格的密度就會(huì)很大，會(huì)使計(jì)算效率明顯降低。通常情況下，柵格的邊長l取1 m 即可。

2.3 放射源源強(qiáng)反解

假設(shè)輻射環(huán)境中總共存在有N個(gè)探測點(diǎn)以及M個(gè)放射源，按照點(diǎn)核積分理論，每個(gè)探測點(diǎn)處的劑量率應(yīng)為各放射源貢獻(xiàn)之和，即

1)列主元消去算法。該算法是為減小方程組求解過程中的舍入誤差而提出的一種算法，與普通消去方法不同，其在每一次消元前都要進(jìn)行一次行替換，以避免絕對值很小的元素直接作為除數(shù)的情況出現(xiàn)。計(jì)算過程可大致分為選主元、行替換、消元、回代求解4個(gè)步驟。

2)NNLS算法。針對方程組Ax=b，其非負(fù)線性最小二乘問題在數(shù)學(xué)上一般表示為

2.4 可信度計(jì)算

需要指出的是，由于點(diǎn)核積分計(jì)算結(jié)果與實(shí)際劑量率總會(huì)存在一定的差異，因此此處的可信度并非絕對可信度，而是相對可信度。當(dāng)遍歷完所有的柵格后，即可選取其中可信度最高的柵格作為最終的放射源定位結(jié)果。

2.5 并行加速

盡管已經(jīng)將輻射區(qū)域進(jìn)行了離散，但當(dāng)輻射區(qū)域面積較大時(shí)，劃分的柵格數(shù)目仍然會(huì)很多，從而需要進(jìn)行大量的復(fù)雜方程組求解。因此為了提高計(jì)算效率，本文引入OpenMP[14]并行機(jī)制來加快計(jì)算速度。OpenMP 是一種適用于多核CPU的并行編程接口，對于需要進(jìn)行循環(huán)的結(jié)構(gòu)塊，OpenMP可以通過派生分支線程的方式來使其同時(shí)進(jìn)行計(jì)算，且當(dāng)所有的分支線程執(zhí)行完成后，才會(huì)執(zhí)行下一語句，不同線程里的結(jié)構(gòu)塊并不會(huì)發(fā)生沖突，其加速效率取決與CPU 核心的數(shù)目。引入OpenMP并行機(jī)制后的完整放射源定位方法流程如圖3所示。

圖3 放射源定位方法流程

3 實(shí)驗(yàn)驗(yàn)證及其分析

3.1 單放射源定位實(shí)驗(yàn)

為驗(yàn)證本文所提方法的可行性及有效性，通過蒙特卡羅程序MCNP5設(shè)計(jì)了如下虛擬仿真實(shí)驗(yàn)：在20 m×20 m 的輻射環(huán)境中，建立笛卡爾坐標(biāo)系，環(huán)境中心處為(0 m,0 m)，其中有一放射源S1，實(shí)際位置為(4.5 m,4.5 m)，其核素為Co-60，產(chǎn)生2 種能量的γ射線，分別為1.173、1.332 MeV，場景中的障礙物為普通混凝土墻，厚度均為15 cm，密度取2.56 g/cm3，并隨機(jī)選取8個(gè)探測點(diǎn)，如圖4所示。

圖4 單放射源定位實(shí)驗(yàn)

在隨機(jī)設(shè)定好放射源的源強(qiáng)后，即可直接通過MCNP5程序?qū)Ω魈綔y點(diǎn)劑量值進(jìn)行計(jì)算，因?yàn)閷?shí)際探測器會(huì)有一定的不確定度，因此將計(jì)算結(jié)果上下浮動(dòng)10%，從而更好地接近真實(shí)情況，最終結(jié)果見表1。

表1 探測點(diǎn)及劑量率

利用C++編程語言實(shí)現(xiàn)本文算法，并根據(jù)上述設(shè)計(jì)好的案例進(jìn)行實(shí)驗(yàn)。柵格劃分的邊長設(shè)為1 m，計(jì)算結(jié)果經(jīng)處理后如圖5所示，計(jì)算部分耗時(shí)約為0.15 s。仿真實(shí)驗(yàn)是在Windows10、64 位操作系統(tǒng)環(huán)境中進(jìn)行的，處理器為Intel Core i5-9300H 2.40 GHz。

圖5 單放射源定位實(shí)驗(yàn)結(jié)果

從圖5中可以看出，當(dāng)柵格邊長為1 m 時(shí)，最終確定的放射源位置為(4 m，4 m)，非常逼近放射源的實(shí)際位置，且可信度較高的其他柵格也均處于實(shí)際位置的周圍，證明了本文所提方法對單放射源定位有著較高的可行性。

3.2 雙放射源定位實(shí)驗(yàn)

與單放射源定位實(shí)驗(yàn)相比，其他條件不變，放射源S1的位置仍為(4.5 m,4.5 m)，另外在區(qū)域左下角(-2.6 m,-3.8 m)處增加一個(gè)放射源S2，核素種類為Cs-137，其γ光子能量為0.661 MeV，如圖6所示。

圖6 雙放射源定位實(shí)驗(yàn)

通過MCNP5程序?qū)Ω魈綔y點(diǎn)進(jìn)行計(jì)算并隨機(jī)上下浮動(dòng)10%后的劑量率見表2。

表2 探測點(diǎn)及劑量率

柵格劃分的邊長仍然為1 m，計(jì)算部分耗時(shí)約為25 s，其中可信度最大的前8種位置如表3所示，表中總偏差為各放射源定位結(jié)果和實(shí)際位置的偏差之和。

表3 雙放射源定位實(shí)驗(yàn)結(jié)果

從表3中可以看出，最終確定的放射源S1、S2的位置分別為(4.0 m,4.0 m)和(-3.0 m,-4.0 m)，與放射源的實(shí)際位置都比較貼近，且可信度較高的其他幾種位置也均處于實(shí)際位置的周圍，證明了本文方法對多放射源定位同樣可行、有效。

4 結(jié)論

本文將點(diǎn)核積分引入放射源的定位之中，提出了一種新的放射源定位方法，并通過MCNP5程序進(jìn)行了虛擬仿真實(shí)驗(yàn)，通過實(shí)驗(yàn)結(jié)果分析可以得出以下結(jié)論：

1)該定位方法需要的探測點(diǎn)數(shù)量較少，且能在有障礙物遮擋情況下對放射源進(jìn)行定位，比已有的放射源估計(jì)定位方法有著更高的適用性；

2)該定位方法對單放射源和多放射源都有著較高的可行性，但隨著放射源數(shù)量的增加，所需要的計(jì)算時(shí)間也會(huì)快速增長，因此在后續(xù)的工作中需要對該算法繼續(xù)進(jìn)行優(yōu)化，以進(jìn)一步降低計(jì)算所需要的時(shí)間。