管 弦 魏 星 李子錕,2 樊海軍 張濟(jì)鵬 孫 濤
1(國民核生化災(zāi)害防護(hù)國家重點實驗室 北京 102205)
2(成都理工大學(xué) 地學(xué)核技術(shù)四川省重點實驗室 成都 610059)
隨著放射源在能源、工業(yè)、農(nóng)業(yè)、醫(yī)療等領(lǐng)域的廣泛應(yīng)用,輻射安全要求也日益嚴(yán)格,截至2021年,國內(nèi)生產(chǎn)、銷售、使用放射性同位素和射線裝置的單位有92 453家,目前,在用放射源共156 539枚,當(dāng)年共發(fā)生5 起輻射事故,其中放射源丟失事故就占3起[1]。遺失放射源對社會與公眾的生命財產(chǎn)安全造成嚴(yán)重威脅,而常見放射源通常較小,發(fā)生丟失事故后只能使用輻射探測儀定位。
傳統(tǒng)放射源定位方法是通過固定式監(jiān)測節(jié)點或者便攜式輻射測量儀遍歷待測區(qū)域,耗時長且效率低,需要專業(yè)人員進(jìn)入未知輻射水平區(qū)域內(nèi)操作,難以保障工作人員的生命健康安全。目前已有較多放射源快速定位相關(guān)研究,張振朝等[2]設(shè)計了一種由NaI、CsI、鍺酸鉍(Bi4Ge3O12,BGO)3種晶體與鉛耦合組成的γ 射線方向探測器,利用各晶體性能差異特點,快速確定γ射線入射方向。周楊等[3]設(shè)計了一套基于3 種NaI(Tl)晶體的γ 射線源定位陣列系統(tǒng),通過建立各探測器全能峰占比與入射光子角度間關(guān)系函數(shù)實現(xiàn)對放射源快速定位。楊靜遠(yuǎn)等[4]設(shè)計了多探測器放射源定位系統(tǒng),基于核電廠環(huán)境監(jiān)測數(shù)據(jù),建立指紋定位方法。Willis 等[5]設(shè)計一套四陣列式NaI(Tl)放射源定位裝置,估計入射射線角度。Durbin 等[6]使用機(jī)器學(xué)習(xí)算法GLAM,基于“十字形”陣列NaI(Tl)探測器對放射源進(jìn)行快速定位。Fragkos等[7]利用小型CZT探測網(wǎng)絡(luò),實現(xiàn)了對無屏蔽放射源的精準(zhǔn)定位。眾多研究成果表明,利用陣列探測器可實現(xiàn)放射源的快速、精準(zhǔn)定位。
本文針對快速尋源定位問題,設(shè)計了一套NaI(Tl)探測陣列放射源定位裝置,利用蒙特卡羅方法,模擬4 個探測器計數(shù)率隨入射射線角度的變化曲線,基于支持向量機(jī)SVM(Support Vector Machine)回歸擬合獲得計數(shù)率-角度響應(yīng)模型,使用模擬數(shù)據(jù)初步驗證模型的有效性,并開展了放射源角度定位與位置定位物理實驗,進(jìn)一步驗證了探測裝置和定位模型可精準(zhǔn)、快速定位放射源位置。
探測陣列放射源定位裝置由4個10 cm×10 cm×10 cm 的NaI(Tl)閃爍探測器組成,裝置結(jié)構(gòu)如圖1所示。4個探測器呈“田”字型緊靠排列擺放在水平面上,探測器之間。在放射源定位過程中,源與探測裝置的距離通常較遠(yuǎn),遠(yuǎn)大于10 倍探測裝置尺寸(2 m),此時放射源對探測裝置所展開的空間角較小,因此,可近似認(rèn)為到達(dá)探測裝置表面的源γ射線為平行束[8-9]。方向為點源到探測裝置中心的連線向量,此時各探測器對源射線的響應(yīng)只與入射角度和能量相關(guān)。
圖1 NaI(Tl)探測陣列放射源定位裝置結(jié)構(gòu)示意圖Fig.1 Structure diagram of array NaI (Tl) detector device for radiation source location
為探究入射角度和探測器響應(yīng)之間的關(guān)系,利用蒙特卡羅模擬軟件SuperMC,開展平行束入射數(shù)值模擬,模擬模型如圖2 所示。以陣列探測器中心為原點建立直角坐標(biāo)系,以象限分別對探測器進(jìn)行編號。為提升采樣效率,降低模擬方差,源項設(shè)置為半徑0.5 m的圓形面源,面源中心與原點間的連線為源平面中心軸法線,連線長度為50 cm,取法線反方向與y軸正方向的夾角為入射射線參考角度θ,初級射線為發(fā)射角度(θ)的γ 平行束,模擬能量為0.661 MeV,發(fā)射粒子總數(shù)為109個,面源與探測器間為標(biāo)準(zhǔn)大氣壓空氣。圖3分別展示了θ分別為60°和210°時,γ射線在各探測器內(nèi)的碰撞分布。
圖2 平行束入射蒙卡模擬模型示意圖Fig.2 Diagram of Monte Caro simulation model with parallel beam γ-ray
圖3 γ射線在陣列探測器內(nèi)的碰撞分布 (a) θ=60°,(b) θ=210°Fig.3 Collision distribution of particles in array detector (a) θ = 60°,(b) θ = 210°
從圖3 中可得出,平行束入射時,4 個探測器之間互相屏蔽,當(dāng)初級γ 射線入射角度發(fā)生變化后,γ射線在各探測器內(nèi)發(fā)生碰撞概率存在明顯變化??梢?,4 個探測器所測得能譜計數(shù)與入射射線角度存在一定關(guān)系,因此,模擬入射角度θ范圍[10°~360°],每隔10°取一發(fā)射角,4個探測器的全譜總計數(shù)隨角度分布如圖4所示。
圖4 陣列探測器模擬總計數(shù)隨入射射線角度分布曲線Fig.4 Distribution curves of MC simulation gross counts with γ-ray incidence angles
圖4縱坐標(biāo)為模擬各探測器發(fā)生能量沉積的光子數(shù),從中可得出,對所有探測器,其在固定角度范圍內(nèi)呈較為平滑的變化趨勢,4 個探測器的總計數(shù)與入射射線角度存在唯一對應(yīng)關(guān)系,且總計數(shù)與入射射線角度θ呈一定函數(shù)關(guān)系。但當(dāng)θ變化至90°、180°、270°、360°處時,探測器總計數(shù)變化趨勢發(fā)生明顯變化,這是由陰影效應(yīng)所導(dǎo)致,當(dāng)θ處于[90°,180°]范圍內(nèi),進(jìn)入1號探測器的部分初級γ射線被4號探測器所屏蔽,進(jìn)入4 號探測器的初級射線則沒有被屏蔽,而當(dāng)θ處于[0°,90°]范圍內(nèi)時則相反。致使在二維空間內(nèi)難以通過單個函數(shù)描述總計數(shù)和入射角度對應(yīng)關(guān)系,為解決這一問題,本文使用支持向量線性回歸方法對入射射線角度進(jìn)行預(yù)測。
支持向量機(jī)(Support Vector Machine,SVM)是一種基于統(tǒng)計學(xué)習(xí)理論的前饋網(wǎng)絡(luò)學(xué)習(xí)方法,其通過核函數(shù),將數(shù)據(jù)從原始低維空間映射至高維歐式空間,建立超平面帶,在高維空間內(nèi)實現(xiàn)對數(shù)據(jù)線性回歸,從而解決非線性回歸問題[10]。
假定現(xiàn)有一組測量數(shù)據(jù),包含各探測器總計數(shù)c與對應(yīng)入射射線角度θ,以改組數(shù)據(jù)為訓(xùn)練集,其如式(1)所示。
式中:R為實數(shù)空間;R4為四維實數(shù)空間,對應(yīng)4個探測器計數(shù);l為測量點數(shù)量。
從圖4 中可看出,此數(shù)據(jù)集無法使用線性函數(shù)描述輸入和輸出的對應(yīng)關(guān)系,因此,將原始數(shù)據(jù)映射到高維空間,具體如式(2)所示[11]。
經(jīng)過映射Φ(xi)后,數(shù)據(jù)集TΦ在高維空間H內(nèi)為近似線性數(shù)據(jù)。
給定高維空間的H一個超平面,其表達(dá)式為[11-12]:
式中:ω為權(quán)重向量;b為偏置。將該超平面沿Y軸方向上下平移ε形成一個超平面帶,使得所有的數(shù)據(jù)點都被包含在內(nèi),將非線性回歸問題轉(zhuǎn)化為凸優(yōu)化問題,如式(4)所示[11-12]。
式中:加和為ε-不敏感訓(xùn)練誤差的總和;C為懲罰因子,用于權(quán)衡訓(xùn)練誤差和懲罰項‖ω‖2;ξi和ξi' 為兩個非負(fù)松弛變量,其代表允許數(shù)據(jù)點越過超平面帶的程度[12-13]。
為求解該最優(yōu)化問題,引入拉格朗日函數(shù),如式(5)所示[12-13]。
式(5)中αi、a'i、ξi、ξ'i為拉格朗日乘子,將上式轉(zhuǎn)化為對偶問題,如式(6)所示[13]。
在[0,C]范圍內(nèi)尋得最優(yōu)解,則凸二次規(guī)劃問題的最優(yōu)解如式(7)所示[14]。
角度回歸預(yù)測函數(shù)如式(8)所示。
為降低運算維度,使用高斯核函數(shù)代替映射內(nèi)積,則角度回歸預(yù)測函數(shù)如下所示[14-15]。
使用入射角度在[10°,180°]的模擬各探測器總計數(shù)作為訓(xùn)練集,采用高斯核函數(shù),使用留一驗證法獲得最佳核函數(shù)參數(shù),建立線性回歸模型。
采用各探測器能譜計數(shù)占比作為預(yù)測模型輸入,真實角度作為輸出,如式(10)所示:
式(10)中,cji為第i個測量點的j號探測器能譜計數(shù)。由于陣列探測器的尺寸相對于離源距離較小,可認(rèn)為在源與探測器的距離遠(yuǎn)大于探測器尺寸時,各探測器對放射源的張角以及距離一致,此時能譜計數(shù)占比只與入射射線角度相關(guān)。則計數(shù)占比誤差為:
式中:ΔN為能譜計數(shù)標(biāo)準(zhǔn)誤差。
將計數(shù)占比與真實角度作為訓(xùn)練集,其表達(dá)式如式(12):
式中:n為訓(xùn)練集內(nèi)的測量點數(shù)量;Tx為訓(xùn)練集輸入,其是由4個探測器能譜計數(shù)之比構(gòu)建的矩陣;Ty為訓(xùn)練集輸出,其是由各測量點對應(yīng)真實角度組成的列向量。
建立模型過程中,從訓(xùn)練集內(nèi)隨機(jī)抽取一個測量點的計數(shù)占比pi和真實角度θi為測試集,用于測試預(yù)測模型的輸出與真實輸出Ty間的誤差,多次迭代直至達(dá)到誤差最小,模型公式如式(13)所示:
式中:Tx,i為第i個測量點4個探測器能譜計數(shù)之比組成的列向量;θ為預(yù)測模型輸出角度;μ為高斯核函數(shù)參考向量;σ為高斯核函數(shù)參數(shù);αˉ為最優(yōu)權(quán)重向量;bˉ為最優(yōu)偏置。μ、σ、αˉ、bˉ均通過求解凸優(yōu)化問題得到[9]。通過此回歸模型即可對射線入射角度進(jìn)行預(yù)測。
使用蒙特卡羅模擬[10°~180°]范圍內(nèi)的模擬各探測器總計數(shù)能譜計數(shù)之比作為訓(xùn)練集獲得模型最優(yōu)參數(shù),在[190°~360°]范圍內(nèi)的模擬各探測器能譜計數(shù)之比作為測試集,測試入射角度預(yù)測模型,測試結(jié)果如表1所示。
表1 MC模擬放射源角度預(yù)測結(jié)果(190°~360°)Table 1 Angle predictions of MC simulation radioactive source
從表1 中可得出,基于SVM 角度回歸預(yù)測模型,最大角度絕對偏差為2.53°,平均絕對偏差為0.52°,對于模擬數(shù)據(jù)可精確估計入射射線角度。
為驗證角度回歸預(yù)測模型定位效果,開展放射源定位物理實驗,探測裝置由4 個10 cm×10 cm×10 cm的NaI(Tl)探測器陣列組合而成。探測器安裝于一旋轉(zhuǎn)平臺之上,平臺可在水平方向360°旋轉(zhuǎn),放射源固定在一線性軌道之上,軌道長10 m,可遠(yuǎn)程操控放射源至軌道任意位置。實驗裝置實物如圖5所示。
圖5 探測器裝置實物圖(a) NaI(Tl)陣列探測器,(b) 陣列旋轉(zhuǎn)臺Fig.5 Physical photos of array NaI(Tl) detector radiation source location device(a) Array NaI(Tl) detector, (b) Array rotation platform
圖6 線性軌道實物圖(a) 軌道整體,(b) 放射源托架平臺Fig.6 Physical diagram of linear orbits(a) Orbits, (b) Radioactive source bracket platform
為獲取陣列探測器對入射射線角度的響應(yīng)關(guān)系,開展放射源角度響應(yīng)實驗,源位于軌道中心5 m處,探測器與放射源的連線與軌道垂直,為確保入射射線為平行束,兩者直線距離為4 m,調(diào)整探測器與源在同一水平面高度,通過旋轉(zhuǎn)平臺改變陣列探測器與放射源的相對角度。實驗裝置擺放以及測量過程如圖7所示。
實驗過程中,記錄4個探測器的總計數(shù),單次測量時間為3 min,每次測量完畢后將平臺沿逆時針方向旋轉(zhuǎn)10°,圖7 展示了θ=0°時各探測器位置分布。實驗結(jié)果如圖8所示。
圖8 放射源角度響應(yīng)實驗結(jié)果Fig.8 Experiment results of radioactive angle response
從圖8 中可得出,實驗測得各探測器總計數(shù)隨入射射線角度的變化趨勢與模擬結(jié)果基本一致,然而,由于實驗測量總計數(shù)受放射源活度和放射源離探測器距離影響,難以直接通過探測器能譜總計數(shù)建立線性回歸預(yù)測模型,考慮到輻射事件的獨立性,對于固定能量的γ 平行束,各探測器總計數(shù)與所有探測器凈計數(shù)和之比與射入探測器的γ射線注量無關(guān),因此實驗采用的訓(xùn)練集n如式(14):
式中:Ni為第i號探測器的能譜總計數(shù),則ni的誤差可表示為:
訓(xùn)練集n隨入射射線角度分布如圖9所示。
圖9 角度響應(yīng)實驗?zāi)茏V計數(shù)占比Fig.9 Proportions of gross count from energy spectra of radioactive angle response experiment
同樣將[10°~180°]范圍內(nèi)的實驗測得各探測器能譜計數(shù)之比和真實角度分別作為訓(xùn)練集輸入和輸出,獲得模型最優(yōu)參數(shù),并建立角度預(yù)測回歸模型,在[190°~360°]范圍內(nèi)的實驗測得各探測器能譜計數(shù)之比作為測試集作為測試集輸入獲得預(yù)測角度,結(jié)果如表2所示。
表2 角度響應(yīng)實驗角度預(yù)測結(jié)果Table 2 Prediction of angle response experiment results
從表2 中可知,角度預(yù)測結(jié)果最大偏差出現(xiàn)在360°處,為6.5°,最小偏差出現(xiàn)在300°處,為0.18°,平均偏差為3.35°,表明在源與探測器距離固定情況下,預(yù)測模型能較為準(zhǔn)確預(yù)測入射射線角度。
為驗證角度預(yù)測回歸模型是否適用于不同距離下的場景,開展放射源位置響應(yīng)實驗,如圖10所示,實驗過程中,以軌道最左端為零點,遠(yuǎn)程操控137Cs點源從軌道左端起點(20 cm)處,移動至軌道右端終點(940 cm)處至軌道中心處,每隔40 cm 取一坐標(biāo)點進(jìn)行測量,單次測量時間為3 min,探測器與放射源平行,記錄4個探測器的全譜計數(shù),實驗過程中分別改變陣列探測器離軌道距離S,實驗結(jié)果如圖11所示。
圖10 放射源位置響應(yīng)實驗示意圖Fig.10 Diagram of radioactive location response experiment
將位置響應(yīng)實驗結(jié)果依照式(14)構(gòu)建測試集輸入,通過模型獲得預(yù)測角度,各距離實驗的角度預(yù)測結(jié)果如圖12。
圖12 位置響應(yīng)實驗角度預(yù)測平均偏差Fig.12 Predicted mean angle deviation of location response experiment
從圖12 中可得出,隨著探測器離軌道距離增大,角度預(yù)測平均偏差也隨之降低,這是由于當(dāng)源距離探測器較近時,點源對探測器的張角較大,此時入射射線不可視為平行束,因此誤差較大,當(dāng)源與探測器距離相對較遠(yuǎn)時(≥3 m),角度預(yù)測平均誤差低于3°,表明平行束入射時,模型能精確預(yù)測入射射線角度。
在實際放射源定位場景中,放射源與定位探測器往往不處于同一水平面,為評估探測器與放射源相對高度差對定位效果的影響,開展探測陣列高度響應(yīng)實驗,實驗過程中,探測陣列與軌道的直線距離調(diào)整為4 m,探測陣列中軸線與軌道平行,放射源移動過程與位置響應(yīng)實驗一致,單次測量時間為3 min,通過操縱陣列旋轉(zhuǎn)臺,將探測陣列相對軌道的高度分別調(diào)至0.3 m、0.4 m、0.5 m、0.6 m,記錄各探測器的全譜計數(shù),實驗結(jié)果如圖13所示。
圖13 高度響應(yīng)實驗結(jié)果 (a) 相對高度0.3 m,(b) 相對高度0.4 m,(c) 相對高度0.5 m,(d) 相對高度0.6 mFig.13 Results of relative height response experiment(a) Relative height 0.3 m, (b) Relative height 0.4 m, (c) Relative height 0.5 m, (d) Relative height 0.6 m
從圖13 和表3 中可得出,不同高度下的探測器全譜計數(shù)響應(yīng)曲線未發(fā)生明顯變化,隨著探測器高度提升,雖然探測器與點源距離增大,但各探測器能譜計數(shù)反而提升,這可能由于γ 射線在垂直方向上不同角度入射方形探測器時,探測器對射線束的探測截面積也隨之增大,接收到初級射線的數(shù)量發(fā)生改變,致使對點源探測效率發(fā)生改變,導(dǎo)致隨探測器高度增高,能譜計數(shù)增大。使用角度預(yù)測模型對實驗數(shù)據(jù)進(jìn)行處理,定位平均偏差如圖14所示。
表3 高度響應(yīng)實驗結(jié)果Table 3 Results of relative height response experiment
圖14 高度響應(yīng)實驗角度預(yù)測平均偏差Fig.14 Predicted mean angle deviation of relative height response experiment
從圖14中可得出,探測陣列與放射源相對高度為0 m 時定位偏差最小,隨著相對高度增加偏差逐漸增大,但變化幅度較小,偏差均小于5°,表明在一定相對高度范圍內(nèi),角度預(yù)測模型能夠較為精確地對放射源進(jìn)行定位。
在實際定位場景中,考慮到探測器不同活度和不同能量的放射源的探測能力存在一定差異,定位能力必然存在一定不同,為探究對不同能量射線源的定位能力,開展蒙特卡羅模擬,模擬模型與圖1中的模型基本相同,分別設(shè)置射線能量為:50 keV、100 keV、200 keV、300 keV、500 keV、750 keV、1 000 eV、1 500 keV、2 000 keV、3 000 keV,模擬粒子數(shù)為109個,得到各探測器的歸一化全譜計數(shù)。使用角度預(yù)測模型定位放射源,結(jié)果如圖15所示。
圖15 各能量射線入射角度預(yù)測平均偏差Fig.15 Predicted mean angle deviation with various energy γ-ray
從圖15 中可得出,低能射線源(50~1 000 keV)的定位角度偏差相對較低,最大值為0.41°,高能射線源(1 500~3 000 keV)的定位角度偏差相對較高,最大值為0.77°,角度偏差整體隨能量而增大,這是由于兩個方面的原因所導(dǎo)致:一是當(dāng)射線能量越高,探測器的探測效率越低;二是當(dāng)射線能量高于1.022 MeV 后,發(fā)生電子對效應(yīng)產(chǎn)生的逃逸光子有可能進(jìn)入相鄰探測并發(fā)生能量沉積,致使響應(yīng)規(guī)律發(fā)生變化,但由于本文設(shè)計的定位裝置晶體體積較大,因此這兩個因素對定位結(jié)果影響較小,整體定位角度偏差較小。
在以上模擬的基礎(chǔ)上,探究探測陣列與定位模型可定位不同能量射線源的最小活度。物理實驗過程中探測器的本底計數(shù)率約為1 600,重構(gòu)后不同活度下探測器全譜計數(shù)如式(16)所示:
式中:N為探測器全譜計數(shù);I為放射源射線強(qiáng)度,為模擬歸一化全譜計數(shù);Poisson(1 600)為服從泊松分布的本底計數(shù)。設(shè)置射線強(qiáng)度范圍為[104,1010],使用不同強(qiáng)度下的探測器全譜計數(shù)作為輸入,獲得定位偏差變化規(guī)律。圖16 為能量為50 keV 時定位偏差隨射線強(qiáng)度的變化曲線。
圖16 定位平均偏差隨γ射線源強(qiáng)度變化Fig.16 Varition of mean localization deviation with the primary γ-ray intensity
從圖16 中可得出,隨著射線強(qiáng)度的增大,定位偏差也逐漸降低,定義定位偏差為5°時的射線強(qiáng)度為最小可定位強(qiáng)度,不同能量的最小可定位強(qiáng)度如圖17所示。
圖17 放射源最小可定位強(qiáng)度Fig.17 Minimum localizable intensity of radiation source
在探測陣列距離放射源4 m 時,源強(qiáng)度達(dá)到最小可定位強(qiáng)度即可確保射入探測陣列的光子所產(chǎn)生的全譜計數(shù)能夠?qū)崿F(xiàn)精準(zhǔn)定位,而探測陣列表面的射線強(qiáng)度與源離探測陣列的距離的平方成反比。因此,基于圖17中最小可定位強(qiáng)度可獲得不同距離下的最小可定位強(qiáng)度。
針對丟失放射源快速定位問題,提出一種基于SVM的放射源角度回歸預(yù)測模型,該模型通過的各探測器能譜計數(shù)之比,快速預(yù)測放射源與探測器間的相對角度。放射源角度響應(yīng)實驗角度預(yù)測結(jié)果顯示,模型預(yù)測角度最大偏差為6.5°,最小偏差為0.18°,平均偏差為3.35°,表明在放射源與探測器距離固定的場景下,回歸模型能夠準(zhǔn)確預(yù)測相對角度;放射源位置響應(yīng)實驗角度預(yù)測結(jié)果顯示,隨著探測器與放射源的距離增大,平均偏差逐漸降低,最大偏差為9.21°,最小偏差為1.77°,說明當(dāng)入射射線越趨近于平行束時,回歸模型預(yù)測精度越高;通過相對高度測試實驗,證明定位模型在一定相對高度差內(nèi)均可實現(xiàn)對放射源的精準(zhǔn)定位;基于不同能量射線入射蒙卡模擬結(jié)果,獲得了不同能量下的最小可定位放射源強(qiáng)度,該方法可快速、精確預(yù)測放射源與陣列探測器的相對角度。
作者貢獻(xiàn)聲明管弦負(fù)責(zé)算法設(shè)計、實驗數(shù)據(jù)收集與處理、文章的起草與修訂;孫濤負(fù)責(zé)研究的提出與設(shè)計和文章的最終修訂;魏星負(fù)責(zé)實驗方案設(shè)計;李子錕負(fù)責(zé)開展實驗與蒙特卡羅模擬;樊海軍負(fù)責(zé)資料調(diào)研;張濟(jì)鵬負(fù)責(zé)算法實現(xiàn)與測試。