基于自適應(yīng)M-H采樣的放射源定位算法

王明生，肖宇峰，劉 冉，劉 成2，楊 川

(1.西南科技大學(xué) 特殊環(huán)境機(jī)器人技術(shù)四川省重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621010；2.西南科技大學(xué) 國(guó)防科技學(xué)院，四川 綿陽(yáng) 621010)

核工業(yè)設(shè)備在我國(guó)得到了廣泛應(yīng)用，據(jù)國(guó)家核安全局統(tǒng)計(jì)發(fā)布的最新數(shù)據(jù)，全國(guó)在用各類放射源 126991枚；各地方放射性廢物庫(kù)已收儲(chǔ)廢舊放射源 40341枚，國(guó)家放射源集中暫存庫(kù)以及由生產(chǎn)廠家回收的廢舊放射源共152109枚[1]。遺失的放射源會(huì)污染環(huán)境，嚴(yán)重威脅人們的身體健康甚至生命安全，快速準(zhǔn)確的放射源定位與搜尋方法是核安全領(lǐng)域重要研究課題。

關(guān)于放射源定位方法，國(guó)內(nèi)外都取得了不少成果。Sun等人[2]介紹了一種使用編碼孔徑相機(jī)的方法估計(jì)放射源位置，但其設(shè)備較重，使用范圍有限。張譯文等人[3]建立了基于多探測(cè)器單元的方向信息和計(jì)數(shù)信息的放射源定位和活度計(jì)算方法，用于區(qū)域內(nèi)的放射源定位和活度監(jiān)測(cè)。左國(guó)平等人[4]設(shè)計(jì)了一種三角圓筒鉛屏蔽的NaI定位探測(cè)器，鉛屏蔽層厚度影響角度分辨能力。Miller等人[5]介紹了一種定位方法，用于機(jī)載平臺(tái)從空中大范圍的放射源定位。Sharma等人[6]介紹了一種基于統(tǒng)計(jì)的網(wǎng)格細(xì)化方法定位放射源；Pahlajani等人[7]通過(guò)部署傳感器網(wǎng)絡(luò)確定放射源，兩者適用于固定場(chǎng)景。Bai等人[8]采用最大似然估計(jì)方法定位放射源，沒(méi)有利用放射源的先驗(yàn)信息。另外一些學(xué)者通過(guò)信號(hào)強(qiáng)度對(duì)射頻標(biāo)簽進(jìn)行定位[9-10]。

小型的個(gè)人劑量?jī)x操作方便、使用靈活，是獲取輻射信號(hào)強(qiáng)度的常用傳感器；同時(shí)在放射源測(cè)量中存在著測(cè)量噪聲以及先驗(yàn)信息不足等不確定因素，而概率統(tǒng)計(jì)是處理此類問(wèn)題的有效方法[11]?；诖?，建立放射源定位的貝葉斯推理模型，提出了改進(jìn)的自適應(yīng)M-H(Metropolis-Hastings)算法實(shí)現(xiàn)放射源的位置估計(jì)。使用Eu-152源設(shè)計(jì)實(shí)驗(yàn)，對(duì)改進(jìn)前后的M-H算法實(shí)驗(yàn)結(jié)果進(jìn)行了驗(yàn)證。

1 輻射測(cè)量中的貝葉斯模型

1.1 輻射測(cè)量原理

在均勻空氣中，點(diǎn)放射源d處的劑量當(dāng)量率H(μSv/h)根據(jù)文獻(xiàn)[12]可由式(1)近似計(jì)算。

(1)

式中,A為放射源活度(MBq)；E為衰變時(shí)釋放的伽馬射線能量(MeV)；d為距離(m)。

式(1)可改寫為H=I/d2，其中I=AE/6，可以看作在離源1 m處的劑量當(dāng)量率。

在輻射測(cè)量中，可以探測(cè)放射源在單位時(shí)間內(nèi)發(fā)射的放射性粒子數(shù)，該過(guò)程存在放射性計(jì)數(shù)統(tǒng)計(jì)漲落現(xiàn)象且服從泊松分布[13-15]。

(2)

式中，z∈N+為實(shí)測(cè)的每分鐘計(jì)數(shù)值(min-1)，λ=Mη，M為理論計(jì)數(shù)值，η為探測(cè)器的效率，λ為期望計(jì)數(shù)值。z可表示為包含一個(gè)能量響應(yīng)常數(shù)Energynumber的算式[12]，即z=H×Energynumber。

1.2 輻射測(cè)量中貝葉斯推理模型

根據(jù)輻射測(cè)量原理，對(duì)環(huán)境中未知放射源的檢測(cè)，可用圖1所示的平面示意圖簡(jiǎn)單說(shuō)明。

圖1 放射源檢測(cè)示意圖

如圖1所示，在輻射場(chǎng)內(nèi)存在r≥1個(gè)未知的點(diǎn)源，源s(s=1,2,…,r)用參數(shù)向量θs=[xs,ys,Is]T∈R2×R+表示，(xs,ys)為放射源的平面坐標(biāo)，Is為放射源1 m處輻照強(qiáng)度計(jì)數(shù)值，所有點(diǎn)源集合表示為θ=[θ1,θ2,…,θr]。gi=[xi,yi,zi]T∈R2×R+表示在該平面某位置(xi,yi)的一個(gè)測(cè)量點(diǎn)，zi為測(cè)量點(diǎn)處的計(jì)數(shù)值。i=1,2，…，m，表示不同的測(cè)量點(diǎn)。假設(shè)各測(cè)量點(diǎn)之間相互獨(dú)立，觀測(cè)向量集G=[g1,g2,…,gi]，則觀測(cè)向量G的似然函數(shù)為

(3)

式中,p為式(2)定義的泊松分布。

p(gi|θ,r)=P(zi;λi(θ))

(4)

(5)

(6)

式(5)中，忽略空氣介質(zhì)的衰減，μb為背景輻射值。

假設(shè)放射源參數(shù)向量θ的先驗(yàn)分布為p(θ|r)，根據(jù)貝葉斯理論可得θ的后驗(yàn)分布為

(7)

式中，p(G|r)可視為常數(shù)，因此聯(lián)合式(3)，式(7)可寫為

(8)

得到后驗(yàn)分布p(θ|G,r)，估計(jì)放射源位置。以估計(jì)位置和真實(shí)位置的距離來(lái)評(píng)估放射源定位的精確度。

對(duì)于一般情形，點(diǎn)源的數(shù)量和位置都是未知的。除了定位放射源的位置，還要估計(jì)放射源數(shù)量?？刹捎脻u進(jìn)估計(jì)的方式，其后驗(yàn)概率計(jì)算公式為

(9)

在p(θ,r|G)基礎(chǔ)上，可得到r個(gè)放射源概率的一般意義計(jì)算公式：

(10)

可得放射源數(shù)目r的MAP估計(jì)值為

(11)

放射源參數(shù)θ的MAP估計(jì)值為

(12)

本文主要研究r已知情況下的放射源定位。

2 自適應(yīng)M-H采樣的放射源定位算法

2.1 M-H采樣算法估計(jì)放射源位置

馬氏鏈蒙特卡羅(MCMC)模擬方法是從某一目標(biāo)分布f(x)中抽取樣本的一類方法。常用采樣方法有Gibbs采樣和M-H采樣。Gibbs采樣應(yīng)用于多維分布的采樣，需要滿足的條件是：多維度分布中任意一個(gè)子分量關(guān)于其他子分量的條件分布已知，以此條件分布作為提議函數(shù)。在放射源位置估計(jì)問(wèn)題中，條件分布p(xs|ys,Is,r)、p(ys|xs,Is,r)、p(Is|xs,ys,r)未知，無(wú)法使用Gibbs采樣。而使用隨機(jī)游走鏈的M-H采樣，通常構(gòu)造提議函數(shù)為正態(tài)分布函數(shù)，正態(tài)分布的性質(zhì)由均值和標(biāo)準(zhǔn)差決定，更易于實(shí)現(xiàn)。

本文采用M-H算法，對(duì)放射源參數(shù)向量θ的后驗(yàn)分布p(θ|G,r)采樣，然后對(duì)樣本進(jìn)行統(tǒng)計(jì)分析，間接得到后驗(yàn)分布p(θ|G,r)的統(tǒng)計(jì)性質(zhì)，完成放射源位置估計(jì)，主要包括以下幾個(gè)步驟。

(1) 初始化采樣點(diǎn)，令t=0。

(2) 采樣新粒子與更新粒子集，當(dāng)t=1,2,…,n時(shí)。

① 構(gòu)造提議函數(shù)。對(duì)于放射源參數(shù)θ，假設(shè)其Markov鏈中上一時(shí)刻狀態(tài)為θt-1。構(gòu)造以θt-1為均值，標(biāo)準(zhǔn)差為σ的正態(tài)分布為提議函數(shù)q(θ*|θt-1)～N(θt-1,σ2)，其標(biāo)準(zhǔn)差σ的選取是關(guān)鍵。

③ 計(jì)算新?tīng)顟B(tài)θ*的接受概率。將采樣得到的新?tīng)顟B(tài)θ*通過(guò)式(8)計(jì)算p(θ*|G,r)。再根據(jù)式(13)，計(jì)算候選狀態(tài)θ*的接受概率α。

(13)

假設(shè)θ中各分量相互獨(dú)立,先驗(yàn)分布p(θ|r)為均勻分布，提議函數(shù)為正態(tài)分布，式(13)可簡(jiǎn)寫為

(14)

④ 更新粒子集。以概率α接受候選狀態(tài)θt=θ*，以概率1-α保持上一時(shí)刻狀態(tài)不變?chǔ)萾=θt-1。

(3) 估計(jì)放射源位置。

對(duì)采樣得到的粒子集θ=[θ0,θ1,θ2，…，θn]進(jìn)行統(tǒng)計(jì)模擬分析，得到θ各分量粒子集的期望E(xs)和E(ys)，以期望值估計(jì)放射源位置。

2.2 自適應(yīng)M-H采樣算法估計(jì)放射源位置

采用隨機(jī)游走鏈的M-H算法，其初始點(diǎn)和提議函數(shù)標(biāo)準(zhǔn)差對(duì)整個(gè)算法的收斂、效率都有較大影響，它主要的不足是這兩者需要依靠人的經(jīng)驗(yàn)來(lái)選取調(diào)整。如果取值不合適，算法執(zhí)行的最終結(jié)果可能收斂緩慢甚至不收斂。

通常，將后驗(yàn)分布p(θ|G,r)附近的點(diǎn)為起點(diǎn)有利于在很短時(shí)間內(nèi)就能得到有效的采樣，提高算法效率。對(duì)于提議函數(shù)標(biāo)準(zhǔn)差的選擇，也有不同的策略。一種策略是依據(jù)一條鏈的生成過(guò)程中對(duì)新?tīng)顟B(tài)的接受率來(lái)選擇。根據(jù)文獻(xiàn)[16]～文獻(xiàn)[18]的研究，不同的接受率會(huì)使隨機(jī)游動(dòng)采樣法的效率受到不同程度的影響。對(duì)于一維變量，最優(yōu)接受率為0.44，將最大限度優(yōu)化隨機(jī)游動(dòng)采樣法的效率，大于五維以后會(huì)降為0.23左右，且與具體分布無(wú)關(guān)。

該階段主要工作是得到一個(gè)采樣起始點(diǎn)θ0，使其后驗(yàn)概率值p(θ0|G,r)大于閾值ρ(ρ≥0)，以此讓采樣的起始點(diǎn)靠近后驗(yàn)分布p(θ|G,r)附近，具體過(guò)程如下。

(2) 自適應(yīng)調(diào)整提議函數(shù)標(biāo)準(zhǔn)差σ。

① 根據(jù)最優(yōu)接受率0.23設(shè)置一個(gè)比較理想的范圍，如區(qū)間 [β-ε,β+ε] (β=0.23，ε為容忍度)。以第(1)階段得到的初始值θ0為采樣起始點(diǎn)，給定標(biāo)準(zhǔn)差初始值σ0，生成一條長(zhǎng)為的L的Markov鏈(實(shí)驗(yàn)部分L取值為3000)，即粒子集[θ0,θ1,θ2，…,θL]，記錄粒子集中接受候選狀態(tài)的次數(shù)為nnew。

② 估算接受比率Pn。

(15)

若P0∈[β-ε,β+ε]，則得到滿足要求的標(biāo)準(zhǔn)差σ=σ0；反之，則以σ0為中心，τ為偏差做一次隨機(jī)游動(dòng)，得到一個(gè)新的標(biāo)準(zhǔn)差σ1，即σ1～N(σ0,τ2)。然后以同樣的方法得到P1，以此類推計(jì)算P2,…,Pn。τ的計(jì)算為

τ=σn×(|Pn-β|÷β)

(16)

τ的計(jì)算依據(jù)是根據(jù)接受比率與理想值的接近程度進(jìn)行調(diào)整。

在不斷調(diào)整標(biāo)準(zhǔn)差的過(guò)程中，還應(yīng)當(dāng)考慮到此類情況，在第n次經(jīng)過(guò)隨機(jī)浮動(dòng)后得到的標(biāo)準(zhǔn)差下，得到的接受率Pn比上一次得到的Pn-1(n>2)更遠(yuǎn)離理想?yún)^(qū)間，即|Pn-β|>|Pn-1-β|，在這種情況，說(shuō)明新得到的標(biāo)準(zhǔn)差σn不如上一次得到的σn-1理想，是應(yīng)當(dāng)舍棄的，同時(shí)接下來(lái)的調(diào)整應(yīng)當(dāng)在兩次中相對(duì)更好的條件下進(jìn)行，也就是保留上次結(jié)果，置σn=σn-1，Pn=Pn-1，然后繼續(xù)下一輪調(diào)整，直到求得滿足接受率的標(biāo)準(zhǔn)差σn為止，保存標(biāo)準(zhǔn)差σ=σn。

(3) 執(zhí)行M-H算法估計(jì)放射源位置。

以第(1)階段中得到的θ0為初始值和第(2)階段中得到的σ為標(biāo)準(zhǔn)差，執(zhí)行M-H算法對(duì)放射源參數(shù)θ的后驗(yàn)分布p(θ|G,r)實(shí)施采樣，對(duì)采樣的粒子集θ=[θ0,θ1,θ2，…，θn]進(jìn)行統(tǒng)計(jì)模擬分析，得到θ各分量粒子集的期望E(xs)和E(ys)，以期望值估計(jì)放射源位置。

3 實(shí)驗(yàn)驗(yàn)證及其分析

3.1 單放射源定位實(shí)驗(yàn)

為驗(yàn)證上述算法的可行性和有效性，設(shè)計(jì)如下實(shí)驗(yàn)：在13 m×19 m的大廳，建立笛卡爾坐標(biāo)系。采用1枚活度約為370 MBq的Eu-152放射源，以北京核儀廠的BH3084型個(gè)人劑量?jī)x為探測(cè)器進(jìn)行實(shí)驗(yàn)。

將放射源完全暴露在空氣中，實(shí)驗(yàn)選取9個(gè)隨機(jī)測(cè)量點(diǎn)，每個(gè)測(cè)量點(diǎn)測(cè)量時(shí)間為2 min，待讀數(shù)穩(wěn)定，用同一個(gè)儀器測(cè)三次取平均值。放射源真實(shí)坐標(biāo)為(4.8 m,4.8 m)，測(cè)量結(jié)果見(jiàn)表1。

表1 實(shí)驗(yàn)測(cè)試點(diǎn)及劑量率

使用Matlab編程實(shí)現(xiàn)M-H算法和自適應(yīng)M-H算法，對(duì)實(shí)驗(yàn)采集的數(shù)據(jù)進(jìn)行處理?？稍O(shè)放射源位置為x∈[1,5]、y∈[1,5]上的均勻分布，迭代4500次。如圖2所示，給出了M-H算法和自適應(yīng)M-H算法的采樣軌跡圖，從圖2(a)中可以看出M-H算法在迭代約500次后進(jìn)入平穩(wěn)狀態(tài)，能估計(jì)出放射源位置。但樣本路徑混合性能較差，表現(xiàn)為在真值附近波動(dòng)小，沒(méi)有對(duì)待估變量的支撐域進(jìn)行充分采樣。在相同初始條件下，依據(jù)2.2節(jié)的描述，設(shè)置自適應(yīng)M-H算法中接受率容忍度ε=0.05。ρ通常是一個(gè)很小的值，取較大值表明采樣初始值靠近后驗(yàn)分布附近，但計(jì)算量大，本文取ρ=2×10-10。w的值根據(jù)采樣空間大小決定，目的是消除采樣隨機(jī)性，減少極端采樣值的影響，文中取w=10。從圖2(b)中可以看出自適應(yīng)M-H算法約迭代100次后進(jìn)入平穩(wěn)狀態(tài)，其收斂速度約為M-H算法的5倍。同時(shí)采樣的樣本具有良好的混合性能，表現(xiàn)為在支撐域附近劇烈的震動(dòng)，在真值附近進(jìn)行充分的采樣。

圖2 采樣軌跡圖

表2給出了實(shí)驗(yàn)中使用的初始值θ0、提議函數(shù)初始標(biāo)準(zhǔn)差σ0和改進(jìn)前后兩種算法的統(tǒng)計(jì)結(jié)果，可得M-H算法的估計(jì)誤差為0.26 m，自適應(yīng)M-H算法的估計(jì)誤差為0.17 m。改進(jìn)后的自適應(yīng)M-H算法估計(jì)結(jié)果更為準(zhǔn)確。

表2 估計(jì)結(jié)果統(tǒng)計(jì)

此外，還希望得到的樣本是獨(dú)立的，為檢驗(yàn)樣本間的自相關(guān)性，如圖3所示，給出了兩種算法得到的樣本自相關(guān)系數(shù)圖。在圖3中，M-H算法得到的樣本自相關(guān)曲線，下降較為平緩，說(shuō)明樣本間有較大的自相關(guān)性，得到的樣本集中獨(dú)立的樣本數(shù)較少，算法效率不高。也說(shuō)明在軌跡圖(見(jiàn)圖2(a))中，樣本曲線在真值附近波動(dòng)小的原因，采樣不夠充分。自適應(yīng)M-H算法得到的樣本自相關(guān)曲線具有很陡的下降曲線，說(shuō)明得到的樣本集中樣本間自相關(guān)性小，采樣更充分，表明與M-H算法取得相同長(zhǎng)度的Markov鏈中包含有更多的獨(dú)立樣本，算法效率更高。

圖3 樣本自相關(guān)系數(shù)圖

3.2 雙放射源定位實(shí)驗(yàn)

同單放射源定位實(shí)驗(yàn)的場(chǎng)地和條件，使用2枚活度約為370 MBq的Eu-152放射源，將其中一枚放入鉛罐，分別放置在s1(4m,4.8m)和s2(4.8m,4.8m)進(jìn)行實(shí)驗(yàn)，測(cè)量結(jié)果見(jiàn)表3。

表3 實(shí)驗(yàn)測(cè)試點(diǎn)及劑量率

單放射源定位實(shí)驗(yàn)已經(jīng)證明了自適應(yīng)M-H算法具有更好的定位效果，所以直接使用自適應(yīng)M-H算法對(duì)雙放射源定位實(shí)驗(yàn)采集的數(shù)據(jù)進(jìn)行處理。具體的分析過(guò)程同單放射源實(shí)驗(yàn)，限于篇幅，不在贅述，

樣本直方圖可以直觀表現(xiàn)出對(duì)各空間位置的估計(jì)，圖4給出了迭代6000次的采樣結(jié)果統(tǒng)計(jì)直方圖。圖4(a)中擬合的正態(tài)曲線均值，即統(tǒng)計(jì)結(jié)果為(3.7 m,4.9 m)、圖4(b) 中擬合的正態(tài)曲線均值，即統(tǒng)計(jì)結(jié)果為(5.0 m,4.9 m)，可分別對(duì)應(yīng)于放射源s1和s2，能夠完成雙放射源的定位，證明了定位方法對(duì)多源的有效性。

圖4 樣本統(tǒng)計(jì)直方圖

4 結(jié)束語(yǔ)

定位點(diǎn)放射源是一個(gè)非常有挑戰(zhàn)的課題，而使用輕量級(jí)的傳感器結(jié)合有效的算法是一個(gè)重要的研究方向。 實(shí)驗(yàn)結(jié)果表明：在室內(nèi)，使用相對(duì)便宜的個(gè)人劑量?jī)x，采用貝葉斯理論結(jié)合統(tǒng)計(jì)模擬方法能夠較好地定位點(diǎn)放射源。所提出的自適應(yīng)M-H算法能對(duì)點(diǎn)放射源位置估計(jì)量的初始值和提議函數(shù)的標(biāo)準(zhǔn)差進(jìn)行自適應(yīng)調(diào)整，在效率和收斂速度方面都優(yōu)于M-H采樣方法。下一步研究工作主要針對(duì)未知數(shù)量的多源、有障礙物情況下的點(diǎn)放射源定位。