基于點核積分的放射源定位方法

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

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

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

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

1 點核積分

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

1.1 計算公式

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

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

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

1.2 累積因子

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

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

2 放射源定位

2.1 問題描述

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

圖1 放射源定位問題示意

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

2.2 區域離散

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

圖2 輻射區域柵格化示意

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

2.3 放射源源強反解

假設輻射環境中總共存在有N個探測點以及M個放射源，按照點核積分理論，每個探測點處的劑量率應為各放射源貢獻之和，即

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

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

2.4 可信度計算

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

2.5 并行加速

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

圖3 放射源定位方法流程

3 實驗驗證及其分析

3.1 單放射源定位實驗

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

圖4 單放射源定位實驗

在隨機設定好放射源的源強后，即可直接通過MCNP5程序對各探測點劑量值進行計算，因為實際探測器會有一定的不確定度，因此將計算結果上下浮動10%，從而更好地接近真實情況，最終結果見表1。

表1 探測點及劑量率

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

圖5 單放射源定位實驗結果

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

3.2 雙放射源定位實驗

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

圖6 雙放射源定位實驗

通過MCNP5程序對各探測點進行計算并隨機上下浮動10%后的劑量率見表2。

表2 探測點及劑量率

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

表3 雙放射源定位實驗結果

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

4 結論

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

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

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