一种模拟能谱曲线仿真方法与流程

本发明涉及一种核能领域，尤其涉及一种模拟能谱曲线仿真方法。

背景技术：

核衰变过程在时间上是随机发生的，衰变过程释放射线(能量)也是随机的，但对其发生的时间间隔及能量值作统计分析，可得知核衰变过程有着以下特性：在发生核衰变的时间间隔上近似服从指数分布；核衰变过程对外释放的能量(即能谱)近似服从高斯分布。

基于核衰变过程存在以上特性，传统方式的仿核信号发生器是以服从不同分布的随机数来模拟核信号特征，即以服从指数分布随机数模拟核信号在时间间隔统计特性；以高斯分布随机数模拟核信号在幅度上的统计特性。但仅以服从高斯分布的随机数来模拟核信号在幅度上的统计特性是不确切的，模拟得到的高期分布曲线与实际能谱曲线存在着较大的误差，不能准确地反映核信号特性；与此同时，各核素的能谱统计特性不尽相同，因此对不同种类核素的幅度特性模拟则需要产生不同参数的高斯分布随机数与之匹配，这是不实际的，在际操作过程中难以实现。

鉴于传统意义上的仿核信号发生器存在的种种弊端，本文提出一种全新方法用以解决以上问题。

技术实现要素：

本发明的目的是提供一种模拟能谱曲线仿真方法，以实现对核能谱线进行仿真。

为了解决上述技术问题，本发明提供了一种模拟能谱曲线仿真方法，其特征在于，包括如下步骤：步骤S1，获取实际核能谱曲线图；以及步骤S2，对实 际核能谱曲线图进行处理以获得模拟能谱曲线。

进一步，所述模拟能谱曲线仿真方法还包括：步骤S3，通过反演比较模拟能谱曲线与实际能谱曲线，以获得模拟能谱曲线与实际能谱曲线之间的误差。

进一步，步骤S2中对实际核能谱曲线图进行处理以获得模拟能谱曲线的方法包括：步骤S21，对实际核能谱曲线图进行曲线识别并将能谱曲线数值化，以得到能谱曲线各点的数值；步骤S22，通过蒙特卡罗方法随机直接抽样这一组数值得到关于各核能级的随机数，以模拟核衰变过程的随机性；以及步骤S23，对所述随机数进行统计处理得到所述模拟能谱曲线。

进一步，所述步骤S21中对实际核能谱曲线图进行曲线识别并将能谱曲线数值化以得到能谱曲线各点的数值的方法包括：将实际核能各谱曲线图的图像经滤波去噪处理后再将实际能谱曲线图进行显示，并根据显示的实际能谱曲线图临摹该核能谱曲线的各关键点获得能谱曲线数据，以建立能谱曲线数据库；或将实际核能各谱曲线图的图像经滤波、降噪预处理、曲线识别、曲线特征提取，以及进行插值处理以完善及修复缺失能谱曲线的各点数据，以建立能谱曲线数据库。

进一步，步骤S22中通过蒙特卡罗方法随机直接抽样这一组数值得到关于各核能级的随机数以模拟核衰变过程的随机性的方法包括：核信号时间统计特性模拟；以及核信号幅度统计特性模拟。

进一步，所述核信号时间统计特性模拟的方法包括：通过服从指数分布的随机数实现核信号时间统计特性模拟，其中

指数分布的随机数由(0，1]均匀分布的随机数通过反函数法变换得到，且(0，1]均匀分布随机数适于通过线性同余法求得。

进一步，所述核信号幅度统计特性模拟的方法包括：

通过对实际核能谱曲线识别并数字化得到各能级幅值及计数率，再通过蒙特卡罗方法直接抽样并输出所述随机数；其中

对实际核能谱曲线的识别并数字化的过程包括：

步骤S221，对实际能谱曲线图进行滤波、降噪；

步骤S222，通过最大类间分割法求出阈值，并将能谱曲线图进行二值化处理，再由像素点扫描法提取出能谱曲线上各点的数值即坐标；

步骤S223，对能谱曲线进行修补及数值化。

进一步，所述通过蒙特卡罗方法直接抽样并输出所述随机数的方法，即

通过蒙特卡罗方法直接抽样能谱曲线及曲线上的各点数值，以获得一系列随机数，从而以模拟核衰变过程的随机性。

进一步，所述步骤S221中对实际能谱曲线图进行滤波的方法，即对实际能谱曲线图进行维纳滤波处理，以滤除能谱曲线图中的高斯噪声。

进一步，所述步骤S223中对能谱曲线进行修补及数值化的方法包括：通过三次样条插值法来填补在能谱曲线特征提取的过程中缺失的数据点，并通过坐标的比例扩伸以获得能谱曲线图上各点的数值。

本发明的有益效果是，本发明的模拟能谱曲线仿真方法通过对实际核能谱曲线图进行曲线识别并将能谱曲线数值化，从而得到能谱曲线各点的数值(即核能谱的能级和各能级的计数率)，再通过蒙特卡罗方法随机直接抽样这一组数值以得到关于各核能级的随机数，从而模拟核衰变过程的随机性，最终再对所述随机数进行统计处理得到模拟能谱曲线，通过反演对比模拟能谱曲线与实际能谱曲线以确定仿核信号发生器的可靠性与精确性。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1是本发明的模拟能谱曲线仿真方法的原理框图；

图2是本发明的模拟能谱曲线仿真方法流程图；

图3是本发明的所述步骤S2中对实际核能谱曲线图进行处理以获得模拟能谱曲线的方法流程图；

图4是本发明所涉及的生成n＝10000个(0,1)均匀分布随机数分布图；

图5是本发明的指数分布随机数分布图；

图6是本发明的对以上指数分布随机数取值范围均匀划分1000个组矩并进行统计的统计图；

图7是本发明的提取到的能谱曲线特征图；

图8是本发明的能谱曲线初步模拟效果图；

图9是本发明的三次样条插值后的效果图；

图10是本发明的最终得到该能谱曲线模拟效果图；

图11示出了模拟核信号随机发生过程的效果图；

图12示出了采用蒙特卡罗方法直接抽样最终效果图。

具体实施方式

现在结合附图对本发明作进一步详细的说明。这些附图均为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。

如图1所示，本发明的模拟能谱曲线仿真方法通过对实际核能谱曲线图进行曲线识别并将能谱曲线数值化，从而得到能谱曲线各点的数值(即核能谱的能级和各能级的计数率)，再通过蒙特卡罗方法随机直接抽样这一组数值以得到关于各核能级的随机数，从而模拟核衰变过程的随机性，最终再对所述随机数进行统计处理得到模拟能谱曲线，通过反演对比模拟能谱曲线与实际能谱曲线以确定仿核信号发生器的可靠性与精确性。

本发明的具体实施方式如以下实施例所示。

如图2所示，本发明的一种模拟能谱曲线仿真方法，包括如下步骤：

步骤S1，获取实际核能谱曲线图；

步骤S2，对实际核能谱曲线图进行处理以获得模拟能谱曲线。

可选的，所述模拟能谱曲线仿真方法还包括：

步骤S3，通过反演比较模拟能谱曲线与实际能谱曲线，以获得模拟能谱曲线与实际能谱曲线之间的误差。

进一步，如图3所示，所述步骤S2中对实际核能谱曲线图进行处理以获得模拟能谱曲线的方法包括：

步骤S21，对实际核能谱曲线图进行曲线识别并将能谱曲线数值化，以得到能谱曲线各点的数值；步骤S22，通过蒙特卡罗方法随机直接抽样这一组数值得到关于各核能级的随机数，以模拟核衰变过程的随机性；以及步骤S23，对所述随机数进行统计处理得到所述模拟能谱曲线。

具体的，所述步骤S21中对实际核能谱曲线图进行曲线识别并将能谱曲线数值化以得到能谱曲线各点的数值的方法包括：

将实际核能各谱曲线图的图像经滤波去噪处理后再将实际能谱曲线图进行显示，并根据显示的实际能谱曲线图临摹该核能谱曲线的各关键点获得能谱曲线数据，以建立能谱曲线数据库；或将实际核能各谱曲线图的图像经滤波、降噪预处理、曲线识别、曲线特征提取，以及进行插值处理以完善及修复缺失能谱曲线的各点数据，以建立能谱曲线数据库。

其中，所述步骤S22中通过蒙特卡罗方法随机直接抽样这一组数值得到关于各核能级的随机数以模拟核衰变过程的随机性的方法包括：核信号时间统计特性模拟和核信号幅度统计特性模拟。

所述核信号时间统计特性模拟的方法包括：通过服从指数分布的随机数实现核信号时间统计特性模拟，其中指数分布的随机数由(0，1]均匀分布的随机数通过反函数法变换得到，且(0，1]均匀分布随机数适于通过线性同余法求得。

具体的，通过线性同余法求得(0，1]均匀分布随机数的方法如下：

线性同余法的递推公式如下：

x_i+1≡λx_i+c(mod M) (1)

其中λ、c为常数。选取的初始x₁称为种子，对随机数的生成质量有一定影响，其值分别在1～2¹⁶＝65535之间选取。为了便于在计算机上使用，通常取

M＝2^S，其中S为计算机中二进制的最大可能有效位数。

图4是取10000个(0，1]随机数分布情况

指数分布随机数的产生方法，即指数分布随机数可由反函数法实现，具体过程如下：

设随机变量X的分布函数服从指数分布：

F(x)＝1-e^-ax，x≥0 (3)

其中，a是一个时间常数，e是自然底数。

由上式可以，F(x)∈[0，1)，且在定义域内单调递减，因此函数F(x)在0～+∞必有反函数，求其反函数：

由于0＜1-F(x)≤1，因此上式可以简化为

由式(5)可得知由符合(0，1]均匀分布的随机数抽样得到服从指数分布随机数x。

取由以上单位均分布随机数通过反函数法产生的指数分布随机数分布图如图5。对以上指数分布随机数取值范围均匀划分1000个组矩并进行统计，最终统计图如图6所示。

所述核信号幅度统计特性模拟的方法包括：通过对实际核能谱曲线识别并数字化得到各能级幅值及计数率，再通过蒙特卡罗方法直接抽样并输出所述随机数；其中对实际核能谱曲线的识别并数字化的过程包括：

步骤S221，对实际能谱曲线图进行滤波、降噪；步骤S222，通过最大类间分割法求出阈值，并将能谱曲线图进行二值化处理，再由像素点扫描法提取出能谱曲线上各点的数值即坐标；步骤S223，对能谱曲线进行修补及数值化。

具体的，所述通过蒙特卡罗方法直接抽样并输出所述随机数的方法，即通过蒙特卡罗方法直接抽样能谱曲线及曲线上的各点数值，以获得一系列随机的随机数，从而以模拟核衰变过程的随机性。

所述步骤S221中对实际能谱曲线图进行滤波的方法，即对实际能谱曲线图进行维纳滤波处理，以滤除能谱曲线图中的高斯噪声，从而尽可能的减小噪声带来的干扰。

所述核信号幅度统计特性模拟的方法的具体实施过程如下：

所述步骤S221中对实际能谱曲线图进行滤波、降噪的具体实施步骤如下：

通过维纳滤波对实际能谱曲线图进行滤波、降噪处理，即所述维纳滤波器为一种线性滤波器，并且还是一种基于最小均方误差准则、对平稳过程的最优 估计器。

假设维纳滤波器输入信号为s(t)，叠加噪声n(t)。输出信号x(t)通过滤波器g(t)由下面卷积运算得到：

x(t)＝g(t)*(s(t)+n(t)) (6)

对于估计出的信号x(t)，期望等同于s(t)。

其误差为：e(t)＝s(t+d)-x(t) (7)

方差为：e²(t)＝s²(t+d)-2s(t+d)x(t)+x²(t) (8)

其中s(t+d)是所期望的滤波器输出。

将x(t)写成卷积积分，即

即可计算平方误差为：

其中R_s是s(t)的自相关函数，R_x是x(t)的自相关函数，R_xs是x(t)和s(t)的自相关函数。维纳滤波的最终目的就是求最优的g(t)，使得E(e²)最小。

所述步骤S222中通过最大类间分割法求出阈值，并将能谱曲线图进行二值化处理，再由像素点扫描法提取出能谱曲线上各点的数值即坐标；

最大类间方差法的具体算法过程如下：

设一幅图像的灰度值为1～m，其中灰度值为i的像素点个数为n_i，N表示图像像素点总数，那么灰度值为i出现的概率为：

令灰度值大于阈值k为c₁组，即c₁＝{1～k}，灰度值大于阈值k的则为c₂组， c₂＝{k+1～m}，那么C₁和C₂出现的概率分别为：

计算得到C₁和C₂的灰度均值为：

其中，那么可得：

μ_r＝ω₁·μ₁+ω₂·μ₂ (16)

由此可计算两组间的方差σ²为：σ²(k)＝ω₁(μ₁-μ_r)²+ω₂(μ₂-μ_r)² (17)

将式(16)代入式(17)可得：σ²(k)＝ω₂ω₂(μ₂-μ₁)²

那么最佳阈值T^*＝Arg max{σ²(k)}，0≤k＜m-1 (18)

求得分割阈值T^*＝0.6353。

所述步骤S223中对能谱曲线进行修补及数值化的具体步骤如下：

实际核能谱曲线图经滤波去噪、二值化后为提取出核能谱曲线上各点的数值即坐标，需提取能谱曲线特征，并将曲线数值化。具体过程如下：

首先，直线识别，即通过扫描核能谱曲线二值图的行和列，识别出核能谱 图中的直线；

其次，定点，由识别的直线判断出能谱曲线所在坐标系的横、纵坐标，并定位原点，一般由上至下，从左往右扫描，识别出的第一条直线就为横、纵坐标；

第三，能谱曲线特征提取。为减少图像中边框及坐标对曲线的影响，需将边框滤除。滤除边框后再由像素点扫描方法逐行或逐列扫描像素点为0的点(二值图像中黑色为0，白色为1)。

最后，曲线数值化。提取出曲线后，通过计算扫描到的能谱曲线有效点到扫描原点的横行与纵行距离确定该像素点在图中的位置，最后通过乘以扩大坐标的比例因子得到该像素点的坐标值。

最终提取能谱曲线特征效果如图7所示。

以及能谱曲线初步模拟效果如图8所示。

进一步，从图7和图8可看到，得到的模拟的能谱曲线图较之原能谱曲线在某些点的数据发生缺失。为尽可能真实地反映实际能谱曲线特性，需要对缺失的数据进行填补修复。

具体的，所述通过三次样条插值法来填补在能谱曲线特征提取的过程中缺失的数据点，并通过坐标的比例扩伸以获得能谱曲线图上各点的数值，以实现对缺失的数据有效地填补与修复。

所述三次样条插值法来填补在能谱曲线特征提取的过程中缺失的数据点的具体算法如下：

定义区间[a，b]上的分段函数S(x)，若满足：

①S(x)在每个子区间[x_i，x_i+1]上是一个三次多项式函数；

②S(x)在整个区间[a，b]上有连续的二阶导数。

则称S(x)为区间[a，b]上关于a＝x₀＜x₁＜…＜x_n＝b的一个三次样条函数。从而三次样条插值问题为：给定函数g(x)的n+1个节点x₀，x₁，…，x_n得函数y₀，y₁，…，y_n，求一个三次样条函数S(x)，使其满足：

S(x_j)＝y_j，j＝0，1，...，n (19)

其中，函数S(x)称为g(x)的三次样条插值函数。

如果S(x)是f(x)的三次样样条插值函数，那么必须满足以下条件：

①插值条件，即

S(x_j)＝y_j，j＝0，1，...，n-1

②连续性条件，即

③一阶导数连续条件，即

④二阶导数连续条件，即

通过三次样条插值后的效果图，如图9所示，由其局部放大图可以看到，三次样条插值后的数据点较为平滑，较为逼近实际值。

实际核能谱曲线模拟效果，即对实际核能谱曲线图经上述图像处理后，最 终得到该能谱曲线模拟效果图如图10所示。

具体的，通过蒙特卡罗方法直接抽样能谱曲线及曲线上的各点数值，以获得一系列随机的随机数，从而以模拟核衰变过程的随机性。

图11示出了模拟核信号随机发生过程的效果图；

图12示出了采用蒙特卡罗方法直接抽样最终效果图(本图由实际能谱曲线图数值化后得到能级和计数率这一数组后，随机抽样过程并统计得到的。本图适于证明由蒙特卡罗抽样的合理性和准确性)。

即通过以上数字图像处理过程得到了模拟能谱曲线及曲线上各点的数值(横坐标为道址Channel，纵坐标为计数率Count)，再以蒙特卡罗方法直接抽样这一组数据就能得到一系列随机的能级(能级例如但不限于通过多道分析器量化得到的，所述道址是指核衰变过程释放的能量经多道分析器量化后得到的)随机数，从而以模拟核衰变过程的随机性。最后再对所述随机数进行统计，可以得到模拟能谱曲线图，这样一方面可以验证系统的可靠性与精确性，另一方面也可以反演于多道分析器，用以标定多道分析器的精确性。

采用蒙特卡罗方法，对某事件A发生的概率P(A)＝p(未知),进行模拟计算，具体计算方法：

(1)进行N次重复独立抽样试验，计算事件A发生次数为n_A。

引入随机变量X_i，表示第i次试验中事件A发生次数，令

则有

(2)计算事件A在N次重复独立抽样试验中的发生频率f_N，为

(3)当N充分大时，以概率f_N作为概率P(A)＝p的估计值为

(4)要求估计值为概率P(A)＝p的无偏估计，即

以及直接抽样，即对核信号在时间和幅度上的特性是以两组服从不同分布的随机数来模拟的，而随机数是离散的，不连续的。对于离散随机序列的抽样，直接抽样法非常理想的。

离散型分布直接抽样法具体抽样过程如下：

设离散随机变量X的取值范围为X_i(i＝0，1，2，3…)，其概率分布为

P(X＝X_i)＝P_i(i＝0,1,2,3……)。其中P_i≥0，

(1)产生(0，1)区间上均匀分布的随机数r；

(2)求得正整数n＝0，1，2...，使得r满足

(3)抽取离散随机变量X的抽样值为X＝X_n。而当0<r≤P₀时X＝X₀；

(4)重复步骤(1)、(2)、(3)直到抽取n个样本值。

由于产生(0,1)均匀分布的随机数r若在区间的概率为

即事件出现的概率等同于事件X＝X_n发生的概率。

又因为随机数r服从(0，1)上的均匀分布，其概率密度函数为

其分布函数如下：

故产生的随机数r抽中样本值为X＝x_n的概率为

由此可知，由直接抽样法抽取到(X＝X_n)的概率等价于随机数X_n在随机数序列X₁，X₂，...X_n出现的频率。

对于直接抽样法可靠性可由以下证明：

设X为离散型随机变量，其概率分布为P_i＝P{X＝X_i}，其中i＝1，2，…。X分别以P_i取得X_i，则事件|X-E(X)|≥ε表示随机变量X取得所有满足不等式|X_i-E(X)|≥ε的可能值X_i，则

∵

∴

∴

由于事件X＝X_i(i＝0，1，2，…N)发生的概率为p_i(0<p_i<1)，则X≠X_i的概率则为1-p_i，而每次X＝X_i发生的概率是不变，且每次抽样结果同其它各次抽取结果无关。因此X＝X_i单个事件是一次贝努利试验，那么抽样n次，则为n重伯努利试验。若令事件A(X＝X_i)发生的次数为n_A，即n_A～B(n，p)。由于X₁，X₂，…，X_n是n个相互独立且服从参数为p的0-1分布的随机变量，且

有D(Xi)＝p(1-p)，i＝1，2，...n。给定任意ε＞0则有

由(4.31)式可推得

而

因此可推得

化简得

即当抽取的次数n越大，抽样后事件A出现的次数与抽样总数的频率比越接近于事件A发生的概率。

由直接抽样抽取随机数其误差为：

令因此

即是p的无偏估计，

即抽样n的次数越大，估计值越接近理论值p。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。