劉浩杰，肖宇峰,*，張 華，田星皓，張 秤

(1.西南科技大學(xué) 特殊環(huán)境機(jī)器人技術(shù)四川省重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621000；2.中國原子能科學(xué)研究院 核技術(shù)應(yīng)用研究所，北京 102413)

核技術(shù)的發(fā)展已走過了100多年的歷程，在國防現(xiàn)代化建設(shè)、工業(yè)、農(nóng)業(yè)、生命科學(xué)、材料科學(xué)、信息科學(xué)、環(huán)境保護(hù)和人民健康等方面發(fā)揮著重要作用，但因?yàn)楸O(jiān)管不當(dāng)導(dǎo)致的放射源丟失也給人類生產(chǎn)和生活帶來巨大威脅。

放射性物質(zhì)的丟失將帶來大場景下搜尋未知放射源的困難，準(zhǔn)確、有效地確定其位置是開展應(yīng)急處置工作的首要條件。針對放射源的定位問題，國內(nèi)外都展開了相應(yīng)的研究，主要方法分為固定式探測和移動(dòng)式探測兩類。固定式探測又分為單點(diǎn)探測和多點(diǎn)探測。單點(diǎn)探測主要圍繞射線檢測研發(fā)各類探測儀器，已有大量研發(fā)成果和產(chǎn)品，適合區(qū)域較小、放射性分布較明確的場景，如三角圓筒鉛屏蔽的NaI探測器[1]、NaI(Tl)晶體探測器[2]、高純鍺探測器[3]、厚GEM探測器[4]等。多點(diǎn)探測通過部署傳感器網(wǎng)絡(luò)來發(fā)現(xiàn)放射源，大量的探測結(jié)點(diǎn)帶來高額的維護(hù)成本，且放射源的定位誤差也較大[5-7]。移動(dòng)式探測通過移動(dòng)平臺(tái)搭載輻射計(jì)數(shù)器來發(fā)現(xiàn)放射源，這類研究較新穎且實(shí)用價(jià)值較高，主要問題是需要人工遙控平臺(tái)探測、定位放射源，而且在大場景下探測未知源時(shí)效率低下[8-9]。綜合來看，固定式探測的準(zhǔn)確性受測量位置和測量區(qū)域的約束，不適合大場景下的未知源搜尋。相比之下，移動(dòng)式探測不受測量位置和區(qū)域的限制，具有很好的靈活性。

針對移動(dòng)式探測遇到的大場景下未知源探測效率低下問題，本文提出一種可應(yīng)用于自主移動(dòng)機(jī)器人上的放射源定位方法：在遞推貝葉斯估計(jì)模型上，通過改進(jìn)的粒子濾波器融合機(jī)器人自定位數(shù)據(jù)和輻射計(jì)數(shù)值，進(jìn)而估計(jì)放射源位置和放射性活度。

1 放射源定位參數(shù)估計(jì)模型

1.1 放射源定位

假設(shè)某區(qū)域內(nèi)存在1個(gè)位置和放射性活度均未知的放射性點(diǎn)源(忽略尺寸大小)，為簡單起見，只考慮二維平面區(qū)域而不計(jì)放射源的垂直高度。為估計(jì)這個(gè)未知放射源的位置和放射性活度，利用自主定位的移動(dòng)機(jī)器人搭載輻射計(jì)數(shù)器進(jìn)行移動(dòng)探測，如圖1所示，進(jìn)而根據(jù)探測信息實(shí)現(xiàn)放射源定位。向量gk=[xk,yk,zk]T表示在該平面某位置(xk,yk)的1個(gè)觀測點(diǎn)，zk為觀測點(diǎn)處的輻射計(jì)數(shù)。其中k∈N*(N*為正整數(shù)集)，表示不同的觀測點(diǎn)。用X=[xs,ys,As]T表示放射源參數(shù)向量，(xs,ys)為放射源的平面坐標(biāo)，As為放射源的活度。

圖1 放射源檢測示意圖Fig.1 Schematic diagram of radioactive source detection

已有研究[10-11]證明，核衰變過程中輻射計(jì)數(shù)器中的計(jì)數(shù)值服從泊松分布，假如某個(gè)位置處的輻射計(jì)數(shù)率為μ，則輻射計(jì)數(shù)器在τ時(shí)間內(nèi)檢測到計(jì)數(shù)值為z∈N(N為自然數(shù)集)的概率P(z;λ)為：

(1)

式中，λ=μτ為泊松分布的均值，且其方差σ2=λ。

如果在觀測位置(xk,yk)用輻射計(jì)數(shù)器測量一段時(shí)間(τk)得到計(jì)數(shù)值zk，則zk的似然函數(shù)為：

P(zk|X)=P(zk;λk(X))

(2)

式中：P(zk|X)為式(1)定義的泊松分布；λk(X)為平均輻射計(jì)數(shù)值，若放射源的活度為As，且每次衰變只發(fā)出1個(gè)光子，光子發(fā)射率為100%，則λk(X)可近似表示為：

(3)

1.2 遞推貝葉斯估計(jì)在放射源參數(shù)估計(jì)中的應(yīng)用

在信號(hào)處理領(lǐng)域，遞推貝葉斯估計(jì)方法可利用帶有噪聲的觀測信息來最優(yōu)地估計(jì)系統(tǒng)的狀態(tài)或參數(shù)，即每接收1個(gè)新的觀測信息，都會(huì)遞歸地更新系統(tǒng)狀態(tài)的后驗(yàn)概率密度函數(shù)，根據(jù)這個(gè)函數(shù)即可得出系統(tǒng)狀態(tài)的最優(yōu)估計(jì)。設(shè)放射源參數(shù)向量X=[xs,ys,As]T在經(jīng)過第k-1次測量值的處理后得到后驗(yàn)分布P(X|z1:k-1)，其中z1:k-1={z1,…,zk-1}為測量值的累積集，當(dāng)接收到第k次測量值zk時(shí)，后驗(yàn)分布更新為P(X|z1:k)，即：

(4)

式中：P(zk|X)為式(2)定義的似然函數(shù)；P(zk|z1:k-1)為歸一化常數(shù)。利用更新后的后驗(yàn)分布可估計(jì)放射源的位置和放射性活度。

2 基于改進(jìn)粒子濾波的放射源定位算法

在濾波過程中，為避免因遞推次數(shù)增加而造成樣本退化現(xiàn)象，一般需進(jìn)行重采樣，但由于傳統(tǒng)重采樣只復(fù)制少數(shù)高權(quán)值粒子而舍棄多數(shù)低權(quán)值粒子，易造成粒子多樣性缺失。本文提出改進(jìn)算法，即高權(quán)值粒子的放射源參數(shù)值由滿足某正態(tài)分布的隨機(jī)數(shù)代替，通過該算法可增加粒子的多樣性、提高定位精度?；诟倪M(jìn)粒子濾波算法估計(jì)放射源參數(shù)的主要過程包括粒子初始化、權(quán)值更新與歸一化、自優(yōu)化重采樣、放射源參數(shù)估計(jì)。

2.1 粒子初始化

(5)

式中：U為均勻分布；(xmin,xmax)、(ymin,ymax)分別為搜索區(qū)域橫向與縱向的搜索范圍；(Amin,Amax)為粒子放射性活度的范圍。

2.2 權(quán)值更新與歸一化

(6)

然后，將新的權(quán)值進(jìn)行歸一化：

(7)

2.3 自優(yōu)化重采樣

在粒子濾波過程中，不可避免地會(huì)出現(xiàn)樣本退化現(xiàn)象，即隨著遞推計(jì)算次數(shù)的增加，粒子權(quán)值的方差會(huì)增大，具體表現(xiàn)為極少數(shù)粒子的權(quán)值很大，而大多數(shù)粒子的權(quán)值幾乎為0，從而導(dǎo)致大量的計(jì)算浪費(fèi)在無用粒子的更新上。為解決樣本退化問題，目前最常用的方法是重采樣，其核心思想是保持粒子總數(shù)不變的情況下，淘汰權(quán)值小的粒子，保留并復(fù)制權(quán)值大的粒子。

在重采樣前，為確定粒子退化程度，引入有效樣本數(shù)Neff，Neff可用下式近似估計(jì)：

(8)

當(dāng)有效樣本數(shù)Neff小于閾值Nthr時(shí)進(jìn)行重采樣，否則保持粒子不變。文獻(xiàn)[12]實(shí)驗(yàn)證明，當(dāng)閾值Nthr設(shè)置為2N/3時(shí)(N為粒子總數(shù))，效果最佳。

在精確度方面，由于系統(tǒng)重采樣的估計(jì)誤差最小[13-15]，所以選擇系統(tǒng)重采樣(SR)算法解決樣本退化問題，其主要過程如下：

雖然重采樣算法對樣本退化問題有一定的抑制作用，但隨著重采樣次數(shù)的增加，由于只復(fù)制高權(quán)值粒子，將導(dǎo)致樣本多樣性的缺失。為克服樣本多樣性匱乏問題，采用優(yōu)化高權(quán)值粒子的改進(jìn)重采樣算法——自優(yōu)化重采樣(A-SR)算法，具體過程如下。

步驟1，首先根據(jù)粒子的權(quán)值判斷有效樣本Neff是否小于閾值Nthr，如果滿足條件，則利用SR算法對粒子進(jìn)行篩選和復(fù)制，否則保持原來的粒子集。

步驟2，考慮到步驟1中可能導(dǎo)致多樣性缺失，對復(fù)制的高權(quán)值粒子進(jìn)行以下優(yōu)化：

(9)

考慮到隨著迭代次數(shù)的增加，粒子集會(huì)越來越收斂，活動(dòng)范圍縮小，所以影響粒子波動(dòng)的方差應(yīng)隨觀測次數(shù)k(即迭代次數(shù))的增加而減小。合適的σp與σr能增加粒子的多樣性，使對放射源位置和活度估計(jì)的精確度得到一定的提高，因此它們的選值是改善粒子濾波算法結(jié)果的關(guān)鍵。

2.4 放射源參數(shù)估計(jì)

(10)

重復(fù)以上過程將使粒子不斷地收斂，放射源的估計(jì)值也會(huì)逐漸向?qū)嶋H值靠攏，當(dāng)粒子空間橫向長度與縱向長度中的最大值小于或等于某個(gè)閾值Dthr時(shí)，即：

(11)

2.5 放射源定位算法實(shí)現(xiàn)步驟

在粒子濾波框架中，引入放射源定位參數(shù)值估計(jì)和自優(yōu)化重采樣思想，完成基于改進(jìn)粒子濾波的放射源定位算法。算法步驟如下。

3 仿真實(shí)驗(yàn)與分析

為驗(yàn)證方法的準(zhǔn)確性和有效性，在matlab軟件下開展仿真實(shí)驗(yàn)，計(jì)算機(jī)采用intel(2.8 GHz)處理器，4 GB內(nèi)存。實(shí)驗(yàn)過程如下：1) 在150 m×150 m平面區(qū)域內(nèi)設(shè)置1個(gè)放射源；2) 假設(shè)移動(dòng)機(jī)器人攜帶輻射計(jì)數(shù)器進(jìn)入該區(qū)域，并探測所處位置點(diǎn)的輻射計(jì)數(shù)值；3) 執(zhí)行改進(jìn)粒子濾波算法，利用測量的位置坐標(biāo)和輻射計(jì)數(shù)值更新粒子集，進(jìn)而得到放射源位置和放射性活度的估計(jì)值；4) 對比分析放射源參數(shù)的設(shè)置值和估計(jì)值，評(píng)價(jià)估計(jì)方法的準(zhǔn)確性和有效性。

3.1 重采樣改進(jìn)后的效果分析

設(shè)置1個(gè)放射性點(diǎn)源的位置坐標(biāo)為(xs=85 m,ys=100 m)、放射性活度As=294 MBq。仿真實(shí)驗(yàn)參數(shù)設(shè)置如下：xmin=0 m、xmax=150 m、ymin=0 m、ymax=150 m、Amin=30 MBq、Amax=1 GBq，有效探測面積S=3.14×10-3m2，(本征)探測效率ρ=30%，本底輻射計(jì)數(shù)率nb=10 s-1，粒子個(gè)數(shù)N=5 000，每個(gè)觀測點(diǎn)處輻射計(jì)數(shù)器的測量時(shí)間相等且滿足τ=10 s。在搜尋放射源的過程中，每記錄1個(gè)觀測值都會(huì)通過濾波過程更新粒子集，并估計(jì)出放射源的位置和放射性活度，然后移動(dòng)機(jī)器人向估計(jì)的放射源方向移動(dòng)15 m到達(dá)新的觀測位置，進(jìn)行新一輪的測量和估計(jì)。

為突出重采樣算法改進(jìn)前后的差異，暫時(shí)不考慮屏蔽物的影響，分別采用SR算法和A-SR算法估計(jì)放射源，其粒子分布情況示于圖2，橫軸和縱軸分別表示二維平面中x方向和y方向上的位置坐標(biāo)。由圖2可看出，經(jīng)過多次迭代估計(jì)，由于SR算法只篩選出高權(quán)值粒子進(jìn)行復(fù)制，導(dǎo)致了多樣性匱乏，這會(huì)使粒子集信息量大幅降低而影響估計(jì)精度。A-SR算法增加粒子多樣性的效果顯著，使粒子集更加稠密。

統(tǒng)計(jì)重采樣改進(jìn)前后的位置誤差和放射性活度誤差，結(jié)果示于圖3。由圖3可知，A-SR算法較SR算法的估計(jì)誤差小，表明改進(jìn)算法提高精度的效果顯著。

圖2 SR算法和A-SR算法估計(jì)的放射源粒子分布Fig.2 Particle distribution by SR and A-SR method

圖3 重采樣改進(jìn)前后的位置誤差和放射性活度誤差分析Fig.3 Analysis of position error and radioactivity error before and after resampling improvement

通過以上對比實(shí)驗(yàn)分別得到了最佳σp和σr值。但當(dāng)σp和σr同時(shí)取最佳值參與計(jì)算時(shí)，會(huì)因?yàn)檫^多地增加粒子的多樣性而影響最后估計(jì)的精度，所以此時(shí)σp和σr的值應(yīng)適當(dāng)下調(diào)。經(jīng)過大量實(shí)驗(yàn)對比發(fā)現(xiàn)，當(dāng)σp=50且σr=45 MBq時(shí)，誤差最小，為驗(yàn)證該結(jié)果，將σp=60且σr=0 Bq、σp=0且σr=50 MBq、σp=50且σr=45 MBq時(shí)估計(jì)的各項(xiàng)誤差進(jìn)行比較，結(jié)果示于圖4e、f。由圖4e、f可知，當(dāng)σp=50且σr=45 MBq時(shí)，估計(jì)值與設(shè)置值的位置誤差以及放射性活度誤差均較小，且誤差波動(dòng)更平緩。以上仿真結(jié)果表明，當(dāng)選取最佳σp和σr參數(shù)時(shí)，提高精度的效果最為顯著。

a——σp影響的位置誤差；b——σp影響的放射性活度誤差；c——σr影響的位置誤差；d——σr影響的放射性活度誤差；e——3種最佳情況下的位置誤差對比；f——3種最佳情況下的放射性活度誤差對比圖4 σp和σr影響的各項(xiàng)誤差對比分析Fig.4 Comparative analysis of error affected by σp and σr

3.2 有無屏蔽物干擾下定位的對比分析

為模擬無屏蔽物干擾情況下的放射源定位，假設(shè)1個(gè)放射性點(diǎn)源位置設(shè)置為(xs=40 m，ys=105 m)，放射性活度As=294 MBq，如圖5所示，其余參數(shù)與3.1節(jié)相同。對于放射源的估計(jì)過程，主要分為2個(gè)階段。階段1：讓搭載輻射計(jì)數(shù)器的可自主定位移動(dòng)機(jī)器人從起點(diǎn)開始進(jìn)行測量和估計(jì)(為驗(yàn)證算法的普適性，設(shè)置2個(gè)不同起點(diǎn)：起點(diǎn)1(10 m，10 m)和起點(diǎn)2(120 m，20 m))，在完成放射源估計(jì)后，移動(dòng)機(jī)器人向估計(jì)源方向移動(dòng)20 m到達(dá)新的觀測位置，繼續(xù)測量和估計(jì)，循序進(jìn)行下去(圖5a)。階段2：當(dāng)機(jī)器人與估計(jì)放射源之間距離小于20 m時(shí)，采取環(huán)繞估計(jì)源移動(dòng)的方式獲取其周圍的測量數(shù)據(jù)，環(huán)繞的方向不定，每次移動(dòng)10 m，直到粒子集滿足條件：max_length≤Dthr=5 m，停止估計(jì)(圖5b)。

在實(shí)際環(huán)境中，搜索區(qū)域一般存在屏蔽物的干擾，如圖6所示，若其中黑色部分是由材質(zhì)相同的混凝土組成的屏蔽物，則根據(jù)文獻(xiàn)[16]可設(shè)其線衰減系數(shù)μ=0.05，這種情況下對放射源的估計(jì)與無屏蔽物干擾時(shí)的估計(jì)過程相同，只是增加了移動(dòng)探測的避障操作。對比圖5、6可發(fā)現(xiàn)，在2種不同環(huán)境中均能較準(zhǔn)確地估計(jì)放射源的位置，但當(dāng)存在屏蔽物干擾時(shí)放射源位置的估計(jì)誤差相對較大。

圖5、6是基于A-SR算法的粒子濾波器在2種不同環(huán)境中對單個(gè)未知放射源進(jìn)行估計(jì)的情況，為達(dá)到對比2種算法定位效果的目的，再利用基于SR算法的粒子濾波器同樣在2種環(huán)境中對未知放射源進(jìn)行估計(jì)，統(tǒng)計(jì)以 (10 m，10 m)為起點(diǎn)進(jìn)行估計(jì)時(shí)這2種算法在2種不同環(huán)境中的估計(jì)數(shù)據(jù)，包括估計(jì)平面位置及誤差、估計(jì)放射性活度及誤差，結(jié)果列于表1、2。

圖5 無屏蔽物干擾下的放射源定位Fig.5 Location of radioactive source without shielding interference

表1 150 m×150 m搜尋區(qū)域內(nèi)無屏蔽物干擾時(shí)未知放射性點(diǎn)源的定位結(jié)果Table 1 Location result of unknown radioactive point source in 150 m×150 m search area without shielding interference

表2 150 m×150 m搜尋區(qū)域內(nèi)存在屏蔽物干擾時(shí)未知放射性點(diǎn)源的定位結(jié)果Table 2 Location result of unknown radioactive point source in 150 m×150 m search area with shielding interference

由表1、2可知：1) 位置誤差和活度誤差表明A-SR算法提高了放射源的估計(jì)精度；2) 基于A-SR算法的粒子濾波器估計(jì)放射源時(shí)，在無屏蔽物干擾的環(huán)境中，位置誤差達(dá)0.04 m，活度誤差達(dá)0.97 MBq，在有屏蔽物干擾環(huán)境中，位置誤差達(dá)1.58 m，活度誤差達(dá)4.99 MBq，屏蔽物對放射源估計(jì)的結(jié)果影響很大；3) 有屏蔽物干擾情況下結(jié)束循序估計(jì)時(shí)的觀測次數(shù)多于無屏蔽物干擾情況下的觀測次數(shù)，同時(shí)，由于A-SR算法增加了粒子多樣性，因此基于A-SR算法的粒子濾波器滿足結(jié)束循序估計(jì)條件時(shí)的觀測數(shù)較基于SR算法的粒子濾波器的觀測次數(shù)多。

4 結(jié)論

針對在復(fù)雜環(huán)境中難以搜尋未知放射源的問題，提出一種基于改進(jìn)粒子濾波的未知放射源定位方法，該方法結(jié)合輻射計(jì)數(shù)器和移動(dòng)機(jī)器人提供的測量數(shù)據(jù)，在線估計(jì)單個(gè)未知放射源的參數(shù)，即位置和放射性活度。仿真實(shí)驗(yàn)結(jié)果表明此方法有效，在150 m×150 m搜尋區(qū)域中對未知放射性點(diǎn)源(活度設(shè)定為294 MBq)進(jìn)行定位，無屏蔽物情況下位置定位誤差約為0.04 m，放射性活度估計(jì)誤差約為0.97 MBq；有屏蔽物情況下位置定位誤差約為1.58 m，放射性活度估計(jì)誤差約4.99 MBq。

實(shí)際應(yīng)用中，放射源的數(shù)量和大小可能均未知，將來還需對本文方法做進(jìn)一步優(yōu)化。另外，可與輻射成像方法結(jié)合以提高定位精度，通過本文方法在較大范圍內(nèi)初步確定放射源，再利用輻射成像工具完成更精確的定位。