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

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

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

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

1 點核積分

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

1.1 計算公式

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

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

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

1.2 累積因子

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

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

2 放射源定位

2.1 問題描述

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

圖1 放射源定位問題示意

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

2.2 區(qū)域離散

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

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

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

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

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

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

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

2.4 可信度計算

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

2.5 并行加速

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

圖3 放射源定位方法流程

3 實驗驗證及其分析

3.1 單放射源定位實驗

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

圖4 單放射源定位實驗

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

表1 探測點及劑量率

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

圖5 單放射源定位實驗結(jié)果

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

3.2 雙放射源定位實驗

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

圖6 雙放射源定位實驗

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

表2 探測點及劑量率

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

表3 雙放射源定位實驗結(jié)果

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

4 結(jié)論

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

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

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