基于自適應M-H采樣的放射源定位算法

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

(1.西南科技大學 特殊環境機器人技術四川省重點實驗室，四川 綿陽 621010；2.西南科技大學 國防科技學院，四川 綿陽 621010)

核工業設備在我國得到了廣泛應用，據國家核安全局統計發布的最新數據，全國在用各類放射源 126991枚；各地方放射性廢物庫已收儲廢舊放射源 40341枚，國家放射源集中暫存庫以及由生產廠家回收的廢舊放射源共152109枚[1]。遺失的放射源會污染環境，嚴重威脅人們的身體健康甚至生命安全，快速準確的放射源定位與搜尋方法是核安全領域重要研究課題。

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

小型的個人劑量儀操作方便、使用靈活，是獲取輻射信號強度的常用傳感器；同時在放射源測量中存在著測量噪聲以及先驗信息不足等不確定因素，而概率統計是處理此類問題的有效方法[11]。基于此，建立放射源定位的貝葉斯推理模型，提出了改進的自適應M-H(Metropolis-Hastings)算法實現放射源的位置估計。使用Eu-152源設計實驗，對改進前后的M-H算法實驗結果進行了驗證。

1 輻射測量中的貝葉斯模型

1.1 輻射測量原理

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

(1)

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

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

在輻射測量中，可以探測放射源在單位時間內發射的放射性粒子數，該過程存在放射性計數統計漲落現象且服從泊松分布[13-15]。

(2)

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

1.2 輻射測量中貝葉斯推理模型

根據輻射測量原理，對環境中未知放射源的檢測，可用圖1所示的平面示意圖簡單說明。

圖1 放射源檢測示意圖

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

(3)

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

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

(4)

(5)

(6)

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

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

(7)

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

(8)

得到后驗分布p(θ|G,r)，估計放射源位置。以估計位置和真實位置的距離來評估放射源定位的精確度。

對于一般情形，點源的數量和位置都是未知的。除了定位放射源的位置，還要估計放射源數量。可采用漸進估計的方式，其后驗概率計算公式為

(9)

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

(10)

可得放射源數目r的MAP估計值為

(11)

放射源參數θ的MAP估計值為

(12)

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

2 自適應M-H采樣的放射源定位算法

2.1 M-H采樣算法估計放射源位置

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

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

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

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

① 構造提議函數。對于放射源參數θ，假設其Markov鏈中上一時刻狀態為θt-1。構造以θt-1為均值，標準差為σ的正態分布為提議函數q(θ*|θt-1)～N(θt-1,σ2)，其標準差σ的選取是關鍵。

③ 計算新狀態θ*的接受概率。將采樣得到的新狀態θ*通過式(8)計算p(θ*|G,r)。再根據式(13)，計算候選狀態θ*的接受概率α。

(13)

假設θ中各分量相互獨立,先驗分布p(θ|r)為均勻分布，提議函數為正態分布，式(13)可簡寫為

(14)

④ 更新粒子集。以概率α接受候選狀態θt=θ*，以概率1-α保持上一時刻狀態不變θt=θt-1。

(3) 估計放射源位置。

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

2.2 自適應M-H采樣算法估計放射源位置

采用隨機游走鏈的M-H算法，其初始點和提議函數標準差對整個算法的收斂、效率都有較大影響，它主要的不足是這兩者需要依靠人的經驗來選取調整。如果取值不合適，算法執行的最終結果可能收斂緩慢甚至不收斂。

通常，將后驗分布p(θ|G,r)附近的點為起點有利于在很短時間內就能得到有效的采樣，提高算法效率。對于提議函數標準差的選擇，也有不同的策略。一種策略是依據一條鏈的生成過程中對新狀態的接受率來選擇。根據文獻[16]～文獻[18]的研究，不同的接受率會使隨機游動采樣法的效率受到不同程度的影響。對于一維變量，最優接受率為0.44，將最大限度優化隨機游動采樣法的效率，大于五維以后會降為0.23左右，且與具體分布無關。

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

(2) 自適應調整提議函數標準差σ。

① 根據最優接受率0.23設置一個比較理想的范圍，如區間 [β-ε,β+ε] (β=0.23，ε為容忍度)。以第(1)階段得到的初始值θ0為采樣起始點，給定標準差初始值σ0，生成一條長為的L的Markov鏈(實驗部分L取值為3000)，即粒子集[θ0,θ1,θ2，…,θL]，記錄粒子集中接受候選狀態的次數為nnew。

② 估算接受比率Pn。

(15)

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

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

(16)

τ的計算依據是根據接受比率與理想值的接近程度進行調整。

在不斷調整標準差的過程中，還應當考慮到此類情況，在第n次經過隨機浮動后得到的標準差下，得到的接受率Pn比上一次得到的Pn-1(n>2)更遠離理想區間，即|Pn-β|>|Pn-1-β|，在這種情況，說明新得到的標準差σn不如上一次得到的σn-1理想，是應當舍棄的，同時接下來的調整應當在兩次中相對更好的條件下進行，也就是保留上次結果，置σn=σn-1，Pn=Pn-1，然后繼續下一輪調整，直到求得滿足接受率的標準差σn為止，保存標準差σ=σn。

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

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

3 實驗驗證及其分析

3.1 單放射源定位實驗

為驗證上述算法的可行性和有效性，設計如下實驗：在13 m×19 m的大廳，建立笛卡爾坐標系。采用1枚活度約為370 MBq的Eu-152放射源，以北京核儀廠的BH3084型個人劑量儀為探測器進行實驗。

將放射源完全暴露在空氣中，實驗選取9個隨機測量點，每個測量點測量時間為2 min，待讀數穩定，用同一個儀器測三次取平均值。放射源真實坐標為(4.8 m,4.8 m)，測量結果見表1。

表1 實驗測試點及劑量率

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

圖2 采樣軌跡圖

表2給出了實驗中使用的初始值θ0、提議函數初始標準差σ0和改進前后兩種算法的統計結果，可得M-H算法的估計誤差為0.26 m，自適應M-H算法的估計誤差為0.17 m。改進后的自適應M-H算法估計結果更為準確。

表2 估計結果統計

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

圖3 樣本自相關系數圖

3.2 雙放射源定位實驗

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

表3 實驗測試點及劑量率

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

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

圖4 樣本統計直方圖

4 結束語

定位點放射源是一個非常有挑戰的課題，而使用輕量級的傳感器結合有效的算法是一個重要的研究方向。 實驗結果表明：在室內，使用相對便宜的個人劑量儀，采用貝葉斯理論結合統計模擬方法能夠較好地定位點放射源。所提出的自適應M-H算法能對點放射源位置估計量的初始值和提議函數的標準差進行自適應調整，在效率和收斂速度方面都優于M-H采樣方法。下一步研究工作主要針對未知數量的多源、有障礙物情況下的點放射源定位。