一种利用贝叶斯分析搜寻本底中弱放射源的方法与流程

1.本发明涉及辐射测定领域，尤其是弱信号识别领域，具体为一种利用贝叶斯分析搜寻本底中弱放射源的方法。


背景技术：

2.在当前国际形势下，核扩散风险依然存在，而全球面对气候变化的挑战也使得核能在能源配置中不可或缺。同时，核燃料交易使得核扩散风险加大。另外，放射源丢失或被盗事件也时有发生。例如，山东省海阳市某水泥有限公司于1998年4月2日，在停产半年后开工检修设备时发现，γ射线料位计的装有
60
co的铅罐被盗，该源活度为3.7
×
108bq；该源丢失9天后，在废品收购点找回。2004年2月6日，西北电力建设第一公司在陕西蒲城丢失一枚
137
cs，该源1969年购于英国，购入时活度为5ci。
3.搜寻放射源或核材料一般可分为两种方式：一种是固定探测器，探测通过物体；另一种是移动探测器，扫描某个给定区域。两种处理方式均涉及到弱信号探测难题，而移动探测器还涉及到本底快速变化的难题。因此，如何实现弱放射源的分析搜索，成为迫切需要解决的难题。
4.为此，本申请提供一种利用贝叶斯分析搜寻本底中弱放射源的方法。


技术实现要素：

5.本发明的发明目的在于：针对弱信号源分析搜索的难题，提供一种利用贝叶斯分析搜寻本底中弱放射源的方法。发明人基于改进后的方法，利用贝叶斯分析对采集的数据进行分析，进而推断源的位置和与探测器运动直线之间的距离，从而给出源是否存在、源的坐标、源强等关键信息。本申请能够用于弱信号的分析搜索，有利于提高搜索效率，得到弱信号的关键信息，具有重要的应用价值。
6.为了实现上述目的，本发明采用如下技术方案：
7.一种利用贝叶斯分析搜寻本底中弱放射源的方法，包括如下步骤：
8.步骤1，建立贝叶斯分析模型，根据数据样本与待推断参数的关系，建立似然函数f(x|θ
i
)；
9.步骤2，提供待推断参数的先验分布π(θ
i
)和数据；
10.步骤3，通过蒙特卡罗推断，获得参数的推断分布；基于推断分布，得到本底中弱放射源存在的概率。
11.所述步骤1中，根据如下公式：
[0012][0013]
建立全部待推断参数的似然函数；
[0014]
式中，b为本底计数；i0为探测器与放射源距离最近时的探测器探测到放射源的计数；x0为源与探测器运动直线最近的点；d为源与x0的距离。
[0015]
所述源与x0的距离d＝m
×
l，表示距离d是步长l的m倍。
[0016]
所述全部待推断参数的似然函数为：
[0017][0018]
其中，f(x|θ
i
)为θ
i
的似然函数，是θ
i
与样本x之间的函数关系；θ
i
为多维参数，在反演过程中均被视为随机变量；π(θ
i
)是随机变量θ
i
的概率分布，即先验分布。
[0019]
所述步骤2中，源与探测器运动直线的距离d为г分布，将d的分布替换为若干个由小到大的固定数值，并分别推断其它3个参数的分布；然后，根据不同d值下，其它三个参数的变化规律，推断出更具有高的置信度的d值，得到待推断参数的先验分布π(θ
i
)和数据。
[0020]
所述步骤3中，利用计算机语言编写蒙特卡罗推断过程，获得参数的推断分布。
[0021]
采用本申请，分析了在平行线扫描搜源过程中，探测器沿直线按一定步长采集的放射性信号。分析结果表明，本申请能识别出本底中极弱放射源信号。在测试时，当探测器探测到的源信号平均值为本底平均值30％时，可以准备识别出包含在其中的源信号；当探测器探测到的源信号平均值为本底平均值50％以上时，不仅可以给出探测到放射源信号的结果，还能在高置信度下给出几乎所有放射源的参数信息。
附图说明
[0022]
本发明将通过例子并参照附图的方式说明，其中：
[0023]
图1为探测器搜源几何关系图。
[0024]
图2为本底中放射源信号随探测器位置的变化图(源位置在序号150附近)。
[0025]
图3为典型的参数分布推断结果。
[0026]
图4为不同强度的源信号与本底叠加后的测量值随探测器位置的变化图(本底平均值均为12，图中源信号为本底平均值的100％)。
[0027]
图5为不同强度的源信号与本底叠加后的测量值随探测器位置的变化图(本底平均值均为12，图中源信号为本底平均值的50％)。
[0028]
图6为不同强度的源信号与本底叠加后的测量值随探测器位置的变化图(本底平均值均为12，图中源信号为本底平均值的30％)。
[0029]
图7为探测器探测到的源信号平均值为本底平均值100％时贝叶斯推断给出的源位置x0。
[0030]
图8为探测器探测到的源信号平均值为本底平均值50％时贝叶斯推断给出的源位置x0。
[0031]
图9为探测器探测到的源信号平均值为本底平均值30％时贝叶斯推断给出的源位置x0。
[0032]
图10为m为5时，x0、b、i0的贝叶斯推断结果图。
具体实施方式
[0033]
本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。
[0034]
本说明书中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的
替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。
[0035]
实施例1
[0036]
1、贝叶斯模型建立
[0037]
探测器测量的总计数t是源信号计数s与本底计数b的和，如下式(1)所示：
[0038]
t＝b+s
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)。
[0039]
给定一探测器，在源探测状态ξ
i
下探测器中源信号s的计数率期望值μ由源强i和距离r决定，令ξ
i
＝i/r2，得到下式(2)：
[0040][0041]
式中，a为探测器截面积；η
e
为探测器本征探测效率；d为给定探测器的本征探测特性参数，依赖于探测截面积和本征探测效率。
[0042]
源信号计数s服从概率密度分布p。当源信号较弱并且本底计数较少时，由于统计涨落的影响探测器单次测量的实际源信号s服从泊松分布p(j；μ
·
t)，是在源探测状态ξ
i
下的p(j；μ
·
t)分布，得到下式(3)：
[0043][0044]
式中，t为测量时间。
[0045]
同理，本底计数b同样服从泊松分布，如下式(4)所示：
[0046][0047]
式中，ω为本底计数率平均值，t为测量时间。
[0048]
因此，探测器总计数t服从如下分布p
t
，如下式(5)所示：
[0049][0050]
根据泊松分布的特性，探测器总计数t分布p
t
仍然是泊松分布，泊松参数为(d
·
ξ
i
+ω)
·
t，得到下式(6)：
[0051][0052]
采用平行线扫描方式进行搜源。当探测器沿直线作步进或连续运动，步进序列为x
i
，步进间隔为l或δl，得到图1的探测器搜源几何关系图。图1中，x是探测器位置，x0是源的真实位置，d为源到探测器运动之线的距离，l为探测器运动步长。
[0053]
每一个测量点对应一个测量状态ξ
i
，任意ξ
i
下探测器总计数t，如下式(7)所示：
[0054][0055]
式中，b为本底计数；i0为探测器与放射源距离最近时的探测器探测到放射源的计数；x0为源与探测器运动直线最近的点；d为源与x0的距离。
[0056]
贝叶斯方法是将先验信息与样本信息结合起来对参数的分布进行推断。在参数空间θ上，π(θ
i
)是随机变量θ
i
的概率分布，即先验分布。在获得样本x后，在x＝x条件下离散θ
i
的后验分布为π(θ
i
|x)，得到下式(8)：
[0057][0058]
求和项为x的边缘分布，用于概率的归一化处理。多维参数θ
i
在反演过程中均被视为随机变量。f(x|θ
i
)为θ
i
的似然函数，是θ
i
与样本之间的函数关系。在已知θ
i
的先验分布π(θ
i
)后，通过建立θ
i
与样本之间的关系模型，即可求出在已知样本信息下的各个参数θ
i
的贝叶斯估计。
[0059]
公式(7)即为各参数贝叶斯推断的似然函数。探测器在进入搜寻区域后，最先获得的数据认为只来自于天然本底。根据公式(8)可知，随着探测器越来越接近x0，源信号越来越强。越过x0后，信号呈对称下降趋势。探测器沿直线前进过程中的泊松分布计数在坐标轴上的分布见图2。图2中，x是探测器位置，b是探测器总计数。
[0060]
通过贝叶斯分析，可以得到b、i0、x0以及d等信息，从而实现本底中弱信号的识别与放射源初步定位。
[0061]
2、结果与分析
[0062]
(1)贝叶斯推断
[0063]
搜源的目的是给出测量数据中放射源的存在概率。当源信号较弱时，公式(4)中的泊松分布的本底b是测量数据的主要成分。如果能给出与放射源相关的3个参数i0、x0和d的概率分布，就可以对放射源存在的概率进行评估。由于探测器计数是服从泊松分布的随机分布，可以利用蒙特卡罗抽样方法进行分析。
[0064]
本申请的推断主要包括3个步骤。
[0065]
第一步，建立贝叶斯分析模型，根据数据样本与待推断参数的关系建立似然函数f(x|θ
i
)。根据公式(8)可知，在源参数的贝叶斯推断前，应给出待推断的参数的似然函数f(x|θ
i
)和参数的先验分布π(θ
i
)。公式(8)包含了全部需要推断的参数与探测器测量样本数据之间的关系，因此，公式(8)可作为全部待推断参数的似然函数。
[0066]
第二步，提供待推断参数的先验分布π(θ
i
)和数据。根据先验信息相容性原理，随着统计数据的无限增多，其先验信息的影响逐渐减弱。选择一个适当的先验分布可以用更少的样本更快速地逼近估计。本底b和源强i0先验分布为泊松分布。源位置x0的先验分布为所有测量点之间的平均分布，因为探测器在运动过程中的任意位置均有可能探测到源信号，它们是等概率的。源与探测器运动直线的距离d为г分布。由于距离平方反比率使得探测器探测到源信号的能力随距离的增加急剧下降，因此更多的应该是关注离探测器相对较近的位置，而距离非常远的源其概率相对小。在实际蒙特卡罗抽样过程中，距离d的分布很难脱离先验分布收敛于真实的d值，说明参数的似然函数中的参数之间的关系相对比较复杂。鉴于此，将d的分布替换为若干个由小到大的固定数值。然后，分别推断其它3个参数的分布，再根据在不同d值下其它三个参数的变化规律，推断哪个d值更具有高的置信度。
[0067]
第三步，利用计算机语言(如：c、python等)编写蒙特卡罗推断过程，获得各参数的推断分布。典型的参数推断结果见图3。图3中，x_1为源坐标，counts_b为本底计数，count_s1为信号计数，frequency为频率(就是概率分布)，sample value为采样值。
[0068]
(2)源是否存在及定位推断
[0069]
为了分析贝叶斯推断在弱放射源识别过程中的灵敏度，将本底平均值控制在12。在此基础上，调节源与探测器的距离，使得源信号在探测器中的最大值分别为本底的
100％、50％和30％。每种源信号强度的测量数据均为200个，见图4～图6所示(图4～图6中，横坐标为探测器位置点的序号，纵坐标为探测器计数)。图4～图6给出了不同强度的源信号与本底叠加后的测量值随探测器位置的变化；本底平均值均为12，图4中源信号为本底平均值的100％，图5中源信号为本底平均值的50％，图6中源信号为本底平均值的30％。
[0070]
为了使源与运动线之间的距离更直观，将距离d替换为m
×
l，即距离d是步长l的m倍。将公式(8)变为如下公式(10)：
[0071][0072]
在不同的距离(即m值)下，源位置x0的推断结果见图7-图9，图中标出了标准不确定度(sd)。
[0073]
由图7可以看出，探测器探测到的源信号平均值为本底平均值100％时，在不同距离m下贝叶斯推断给出的源位置x0在66.3附近。实验中真实的x0＝63。这个参数的贝叶斯推断结果均与实验值较好的吻合。在m由1增大到20的过程中，x0的分布展宽由大变小在变大。在m＝5时，x0的不确定度最小为1.1，表明在m在5附近具有最大概率。也就是说，源与探测器运动直线的距离d为探测器前进步长的5倍左右。在实际实验中，本底平均值为12counts/s，源发射率为12counts/s，探测器前进步长为0.1m，距离d为0.5m，m值为5。m值的推断结果也与实验值相吻合。这时，贝叶斯推断出源信号的发射率i0的分布在6.5counts/s至14.5counts/s之间，平均值为10.5counts/s(见图10)，与实验值12counts/s一致。
[0074]
在图8中，探测器探测到的源信号平均值为本底平均值50％时，贝叶斯推断给出的源位置x0的推断。结果表明，x0的分布在48-63之间，平均值56.1与实验值63较为接近，但是其不确定度比100％时大。m在3、6、和15时的不确定度最小。图8中，通过不确定度变化趋势给出的m值的置信度偏低。
[0075]
图9中，x0的分布在范围扩大至15-65，其离散度偏大，只能给出探测器探测到源信号的存在，但是给不出较为准确的源位置信息。源位置的不确定度变化趋势不明显，因此，无法给出有效的源距离探测器运动直线的距离信息。
[0076]
3、结论
[0077]
本申请中，分析了在平行线扫描搜源过程中，探测器沿直线按一定步长采集的放射性信号。分析结果表明，该方法能识别出本底中极弱放射源信号。当本底平均值为12cts/s时，可以识别出包含在其中的4cts/s源信号。当源信号更强时，不仅可以给出探测到放射源信号的结果，还能在高置信度下给出几乎所有放射源的参数信息，比如：位置x0、距离m和源发射率i0。
[0078]
本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。