基于改進粒子濾波的未知放射源定位方法

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

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

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

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

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

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

1.1 放射源定位

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

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

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

(1)

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

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

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

(2)

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

(3)

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

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

(4)

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

2 基于改進粒子濾波的放射源定位算法

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

2.1 粒子初始化

(5)

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

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

(6)

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

(7)

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

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

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

(8)

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

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

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

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

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

(9)

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

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

(10)

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

(11)

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

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

3 仿真實驗與分析

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

3.1 重采樣改進后的效果分析

設(shè)置1個放射性點源的位置坐標(biāo)為(xs=85 m,ys=100 m)、放射性活度As=294 MBq。仿真實驗參數(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%，本底輻射計數(shù)率nb=10 s-1，粒子個數(shù)N=5 000，每個觀測點處輻射計數(shù)器的測量時間相等且滿足τ=10 s。在搜尋放射源的過程中，每記錄1個觀測值都會通過濾波過程更新粒子集，并估計出放射源的位置和放射性活度，然后移動機器人向估計的放射源方向移動15 m到達新的觀測位置，進行新一輪的測量和估計。

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

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

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

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

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

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

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

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

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

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

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

表1 150 m×150 m搜尋區(qū)域內(nèi)無屏蔽物干擾時未知放射性點源的定位結(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)存在屏蔽物干擾時未知放射性點源的定位結(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算法提高了放射源的估計精度；2) 基于A-SR算法的粒子濾波器估計放射源時，在無屏蔽物干擾的環(huán)境中，位置誤差達0.04 m，活度誤差達0.97 MBq，在有屏蔽物干擾環(huán)境中，位置誤差達1.58 m，活度誤差達4.99 MBq，屏蔽物對放射源估計的結(jié)果影響很大；3) 有屏蔽物干擾情況下結(jié)束循序估計時的觀測次數(shù)多于無屏蔽物干擾情況下的觀測次數(shù)，同時，由于A-SR算法增加了粒子多樣性，因此基于A-SR算法的粒子濾波器滿足結(jié)束循序估計條件時的觀測數(shù)較基于SR算法的粒子濾波器的觀測次數(shù)多。

4 結(jié)論

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

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