放射源定位范例-基于贝叶斯估计和马尔可夫链蒙特卡洛方法的单放射源

作者: 蒜薹 | 来源:发表于2019-04-14 10:29 被阅读0次

放射源定位范例-基于贝叶斯估计和马尔可夫链蒙特卡洛方法的单放射源定位研究

简介

本研究以放射源坐标和强度为未知量，也就是待估算变量，其中，坐标考虑二维平面。放射源数量依然为一个，因此，本研究有三个待估计变量。

在构造概率模型方面，已经没有任何问题了，本研究有能力考虑任意数量且具有任意强度的放射源，任意材料组成和任意分布的空间介质，任意数量、坐标的探测器，并且具有任意有效探测面积和固有特虐效率。但是，在求取最优解的方法上，以前研究所使用的简单蒙特卡罗方法只在待估计变量数量小的时候，具有实际应用价值，随着待估计变量数量的增加，计算量是指数增加的，因此，在面对多变量的概率模型时，需要考虑其他方法。

本研究以马尔可夫链蒙特卡洛(Markov Chain Monte Carlo, MCMC)方法为概率模型的求解方法。这是一种统计学方法，将蒙特卡洛随机抽样的思想结合到马尔可夫链中，构造了一种可以按照已知分布进行随机抽样的方法。MCMC的特点是，所提供的已知分布任意，可以非常复杂，可以离散，也可以连续。在本研究中，这种已知分布就是探测器观测概率模型的似然函数，或者后验函数。

Hite(2019,a)使用似然函数作为已知分布，对给定数量放射源的坐标和强度进行估计，使用高斯分布作为马尔可夫链状态转移矩阵，使用Metropolis-Hastings(M-H)采样提升采样效率。

统计学模型构造

泊松分布和光子的观测概率模型

探测器对单个光子的响应过程遵循二项分布，而在观测过程中，探测器会对数量极大的光子群进行响应，并且在本研究中认为光子之间互相独立，则整个观测过程是遵循泊松分布的。

第Di个探测器对第Si个放射源的响应数量，也就是观测到的 $\gamma$ 光子数量为 $I_{Di,Si,Observed}$ ，而理论上应该观测到的数量为 $I_{Di,Si,Expected}$ ，这种事件的发生概率可以使用泊松分布进行计算。
$P(\boldsymbol r_{Si},I_{Si}|Obs_{Di}) \equiv Poisson(I_{Di,Si,Observed},I_{Di,Si,Expected})$

$P(\boldsymbol r_{Si},I_{Si}|Obs_{Di}) = \frac{{I_{Di,Si,Expected}}^{I_{Di,Si,Observed}}}{I_{Di,Si,Observed}!} \cdot e^{-I_{Di,Si,Expected}}$

单放射源多探测器观测的概率模型

按照本研究的设定，环境中只有一个放射源，探测器的数量不止一个，并且分布于环境中的不同位置。假设经过一次观测，则每一台探测器都有一个属于自己的观测数值，也就是对 $\gamma$ 光子的响应次数。设定，放射源编号S0；探测器编号从D0开始，至DN-1；编号为Di的探测器，所观测到的响应次数为 $I_{Di,S0,Observed}$ ，理论上这台探测器的响应次数为 $I_{Di,S0,Expected}$ 。对于这种设定下的观测结果，可以将每一个探测器的观测事件发生概率进行累积，并且以乘法的方式累积，用以表示本次观测的总概率 $P(\boldsymbol r_{S0},I_{S0}|Obs_{All})$ ，其计算方法如下。
$P(\boldsymbol r_{S0},I_{S0}|Obs_{All}) = \prod_{Di=D0}^{DN-1} {P(\boldsymbol r_{S0},I_{S0}|Obs_{Di})}$

$P(\boldsymbol r_{S0},I_{S0}|Obs_{All}) = \prod_{Di=D0}^{DN-1} {\frac{{I_{Di,S0,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,Expected}}}$

似然函数的构造

似然函数指以观测数值为条件，放射源的第Plani种分布情况为结果的概率，以 ${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0})$ 表示，可以发现，似然函数就是上一节所描述的概率模型， $P(\boldsymbol r_{S0,Plaini},I_{S0}|Obs_{All})$ 。
${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) \equiv P(\boldsymbol r_{S0,Plaini},I_{S0}|Obs_{All})$

${\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) = \prod_{Di=D0}^{DN-1} {\frac{{I_{Di,S0,Expected}}^{I_{Di,S0,Observed}}}{I_{Di,S0,Observed}!} \cdot e^{-I_{Di,S0,Expected}}}$

与最基本的极大似然估计不同，MCMC方法中，不能对似然函数求取对数，也不能省略似然函数中不变的部分，也就是分母部分，而需要全部使用。

马尔可夫链蒙特卡洛方法的实现

MCMC计算流程

马尔可夫链蒙特卡洛(Markov Chain Monte Carlo, MCMC)方法可以分为几个步骤，首先制作带求解变量的初始值；其次进入循环，并随机产生带求解变量；接着，在循环内，对初始值的方案计算似然函数的函数值，对随机产生的方案计算似然函数的函数值；接着，在循环内，对上一步两个似然函数值求比值，并根据比值判断当前循环内随机产生的方案符合似然函数的分布；接着，在循环内，进入下一次循环；最后，在样本数量获取满足要求时，结束循环，并从这些样本中，找出最大似然函数值对应的样本，即最优方案。

值得注意的是，在循环中，仅仅使用似然函数，并未出现后验概率函数，这是Hite(2019,a)的工作内容，并不表示MCMC在做放射源定位时，只能使用似然函数。

initialize
    Ns = 1e8  // desired number of samples from posterior
  xs = []   // lists storing X-Axis coordinates of samples 
  ys = []   // lists storing Y-Axis coordinates of samples 
  Is = []   // lists storing Intensities of samples 
  x0 = 1.   // The starting value for X-Axis coordinates of the initializing sample
  y0 = 2.   // The starting value for X-Axis coordinates of the initializing sample
  I0 = 1e8  // The starting value for X-Axis coordinates of the initializing sample
  xs.APPEND(x0)  // Set the initializting sample as the first sample in the list
  ys.APPEND(y0)
  Is.APPEND(z0)
done

while LENGTH(xs)<Ns do
  // prepare referenced plan and current plan
    sizeOfSamples = xs.SIZE()
    xReferenced = xs[sizeOfSamples-1]
    yReferenced = ys[sizeOfSamples-1]
    IReferenced = Is[sizeOfSamples-1]
    
    xCurrent = NormalDistribution(xReferenced, xSigma)
    yCurrent = NormalDistribution(yReferenced, ySigma)
    ICurrent = NormalDistribution(IReferenced, ISigma)
    
    // get alpha
    alpha = LikeliHood(PlanCurrent)/LikeliHood(PlanReferenced)
    
    // get gamma
    gamma = UniformDistribution(0,1)
    
    // compare
    if gamma<=alpha then
        xs.APPEND(xCurrent)
        ys.APPEND(yCurrent)
        Is.APPEND(ICurrent)
    end if
end while

RETURN xs, ys, Is

MCMC循环方式

循环是MCMC的核心步骤，主要任务是产生方案和筛选方案，方案使用随机数生成器产生，使用似然函数等进行筛选。在本研究中，似然函数为泊松分布的累乘，而在编程计算中发现，泊松分布的计算结果往往是无限大或者无限小。其原因是泊松分布中包含了指数、阶乘，虽然理论计算结果在0至1的区间内，但是在实际计算中，需要分别计算泊松分布的各个部分，恰好就是这几个部分，其独立运算的结果会超过量程。如泊松分布的一个部分如下，当前方案中，若探测器响应次数的期望值为10000，观测值为12000，那么计算结果就达到 $10^{48000}$ ，这已经远远超过了常用计算机硬件配置下C/C++编译器对double型变量的量程。
${I_{Di,S0,Expected}}^{I_{Di,S0,Observed}} = 10000^{12000} = 10^{48000}$
巧合的是，MCMC的Metropolis方法提供了一种数学结构，使用户可以通过求取对数的方法来控制计算数值的范围，以解决这个问题。而Metropolis方法的产生动机，是为了解决MCMC方法效率低下问题。

Metropolis方法改进了MCMC循环部分的方案筛选方法，设定某一次循环内，用于产生随机方案的方案为参考方案，以符号 $Referenced$ 表示，参考方案的似然函数为 ${\cal{L}}_{Obs,PlanReferenced}({\boldsymbol r}_{S0},I_{S0})$ ；被参考方案产生的随机方案以符号 $Current$ 表示，代表当前循环的随机方案，其似然函数为 ${\cal{L}}_{Obs,PlanCurrent}({\boldsymbol r}_{S0},I_{S0})$ 。
$min\{ 1, \frac {{\cal{L}}_{Obs,PlanCurrent}({\boldsymbol r}_{S0},I_{S0})} {{\cal{L}}_{Obs,PlanReferenced}({\boldsymbol r}_{S0},I_{S0})} \}$

Metropolis筛选方法的核心如上，构造了两个似然函数的比例关系。我们可以通过取对数的方式，缩小计算结果的数值范围，使C/C++变量的量程足够使用。
$min\{ ln(1), ln( \frac {{\cal{L}}_{Obs,PlanCurrent}({\boldsymbol r}_{S0},I_{S0})} {{\cal{L}}_{Obs,PlanReferenced}({\boldsymbol r}_{S0},I_{S0})} ) \}$
考虑方案Plani的似然函数的对数形式，
$ln{\cal{L}}_{Obs,Plani}({\boldsymbol r}_{S0},I_{S0}) = \sum_{Di=D0}^{DN-1} (I_{Di,S0,Observed} \cdot ln(I_{Di,S0,Plani,Expected})-I_{Di,S0,Plani,Expected}) -DN \cdot ln(I_{Di,S0,Observed}! )$
当前方案与参考方案的比例的对数形式如下，为了简化表达，省略一些共同符号。
${\cal{L}}_{Obs,PlanCurrent}({\boldsymbol r}_{S0},I_{S0}) \equiv {\cal{L}}_{PlanCurrent}({\boldsymbol r},I)$

${\cal{L}}_{Obs,PlanReferenced}({\boldsymbol r}_{S0},I_{S0}) \equiv {\cal{L}}_{PlanReferenced}({\boldsymbol r},I)$

$I_{Di,S0,Observed} \equiv I_{Di,Observed}$

$I_{Di,S0,Plani,Expected} \equiv I_{Di,Plani}$

${\cal{L}}_{PlanCurrent}({\boldsymbol r},I) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ln(I_{Di,PlanCurrent})-I_{Di,PlanCurrent}) -DN \cdot \sum_{Di=D0}^{DN-1} lnI_{Di,Observed}$

${\cal{L}}_{PlanReferenced}({\boldsymbol r},I) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ln(I_{Di,PlanReferenced})-I_{Di,PlanReferenced}) -DN \cdot \sum_{Di=D0}^{DN-1} lnI_{Di,Observed}$

综合上述简化修改，对数形式如下。
$ln( \frac {{\cal{L}}_{PlanCurrent}({\boldsymbol r},I)} {{\cal{L}}_{PlanReferenced}({\boldsymbol r},I)} ) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ln(I_{Di,PlanCurrent})-I_{Di,PlanCurrent}) -DN \cdot \sum_{Di=D0}^{DN-1} lnI_{Di,Observed} \\ - ( \sum_{Di=D0}^{DN-1} (I_{Di,S0,Observed} \cdot ln(I_{Di,S0,Plani,Expected})-I_{Di,S0,Plani,Expected}) -DN \cdot \sum_{Di=D0}^{DN-1} lnI_{Di,Observed} ) )$
化简得到
$ln( \frac {{\cal{L}}_{PlanCurrent}({\boldsymbol r},I)} {{\cal{L}}_{PlanReferenced}({\boldsymbol r},I)} ) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ln(I_{Di,PlanCurrent})-I_{Di,PlanCurrent}) \\ - (\sum_{Di=D0}^{DN-1} (I_{Di,Observed} \cdot ln(I_{Di,PlanReferenced})-I_{Di,PlanReferenced}) )$

$ln( \frac {{\cal{L}}_{PlanCurrent}({\boldsymbol r},I)} {{\cal{L}}_{PlanReferenced}({\boldsymbol r},I)} ) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ln(I_{Di,PlanCurrent}) - I_{Di,Observed} \cdot ln(I_{Di,PlanReferenced}) ) \\ - \sum_{Di=D0}^{DN-1} ( I_{Di,PlanCurrent} - I_{Di,PlanReferenced} )$

$ln( \frac {{\cal{L}}_{PlanCurrent}({\boldsymbol r},I)} {{\cal{L}}_{PlanReferenced}({\boldsymbol r},I)} ) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ( lnI_{Di,PlanCurrent} - lnI_{Di,PlanReferenced} ) ) \\ - \sum_{Di=D0}^{DN-1} ( I_{Di,PlanCurrent} - I_{Di,PlanReferenced} )$

$ln( \frac {{\cal{L}}_{PlanCurrent}({\boldsymbol r},I)} {{\cal{L}}_{PlanReferenced}({\boldsymbol r},I)} ) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ln \frac{I_{Di,PlanCurrent}} {I_{Di,PlanReferenced}} ) - \sum_{Di=D0}^{DN-1} ( I_{Di,PlanCurrent} - I_{Di,PlanReferenced} )$

$ln( \frac {{\cal{L}}_{PlanCurrent}({\boldsymbol r},I)} {{\cal{L}}_{PlanReferenced}({\boldsymbol r},I)} ) = \sum_{Di=D0}^{DN-1} ( I_{Di,Observed} \cdot ln \frac{I_{Di,PlanCurrent}} {I_{Di,PlanReferenced}} - ( I_{Di,PlanCurrent} - I_{Di,PlanReferenced} ) )$

MCMC估算结果和讨论

MCMC估算结果

MCMC采样结果从第6个样本开始就已经符合真实情况了，放射源坐标均在 $(1m,2m)$ 附近，放射源强度均在 $10^8$ 左右。

仅对坐标制作直方图，由于计算时间短，样本数量少，但是可以看出采样结果均围绕在真实坐标附近。

![img] (
https://github.com/PascalXie/DataForMCMC/blob/master/temp34.png?raw=true)

结果讨论

MCMC是一个有效的算法，但是抽样和判断复杂，导致计算时间长，难以达到实时定位的效果。

缩短计算时间的方法包括采用多线程抽样。

参考文献

[Hite, et al., 2019, a] Jason Hite, John Mattingly, Bayesian Metropolis methods for source localization in an urban environment, Radiation Physics and Chemistry, Volume 155, 2019, Pages 271-274, ISSN 0969-806X, https://doi.org/10.1016/j.radphyschem.2018.06.024. (http://www.sciencedirect.com/science/article/pii/S0969806X17307867)

[Hite, et al., 2019, b] Jason Hite, John Mattingly, Dan Archer, Michael Willis, Andrew Rowe, Kayleigh Bray, Jake Carter, James Ghawaly, Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Volume 915, 2019, Pages 82-93, ISSN 0168-9002, https://doi.org/10.1016/j.nima.2018.09.032. (http://www.sciencedirect.com/science/article/pii/S016890021831177X)

放射源定位范例-基于贝叶斯估计和马尔可夫链蒙特卡洛方法的单放射源

放射源定位范例-基于贝叶斯估计和马尔可夫链蒙特卡洛方法的单放射源定位研究

简介

统计学模型构造

泊松分布和光子的观测概率模型

单放射源多探测器观测的概率模型

似然函数的构造

马尔可夫链蒙特卡洛方法的实现

MCMC计算流程

MCMC循环方式

MCMC估算结果和讨论

MCMC估算结果

结果讨论

参考文献

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读