✨ 长期致力于低分辨率γ能谱、高分辨解析、峰边界确定、自适应半高宽、响应矩阵研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》1自适应半高宽匹配的SNIP本底扣除与反卷积联合寻峰针对NaI(Tl)探测器能谱中低能区本底畸变和重叠峰严重的问题提出一种自适应SNIP本底扣除算法。该算法根据能谱的道址和计数率动态选择变换窗宽m窗宽取值与当前道址处的半高宽FWHM成正比m 2 * FWHM / 道宽。在标准SNIP迭代中对每个道址进行多次削峰操作直到本底趋于平滑。本底扣除后采用反卷积联合寻峰策略首先使用Gold反卷积算法对谱线进行解卷积将重叠峰的FWHM压缩至原始的一半然后在一阶导数零交叉点处识别峰位。对于严重重叠峰峰间距小于1.2倍FWHM引入峰边界约束的高斯混合模型用EM算法估计各子峰参数。在一块238U矿石样品实测谱中该方法成功分离了185.7keV和186.2keV的重叠双峰峰位偏差小于0.3keV。2高斯响应矩阵与Boosted-Gold迭代解谱基于NaI(Tl)探测器对γ射线的响应近似高斯函数这一物理特性构建了一个高斯响应矩阵。矩阵元素R_ij表示能量为Ej的单能γ射线在第i道址产生的计数其中高斯函数的FWHM按照能量刻度公式FWHM a*sqrt(E)b 计算a和b通过241Am、137Cs、60Co等标准源标定。矩阵维度为1024×1024。解谱时采用改进的Boosted-Gold算法与传统Gold不同该算法在每次迭代中引入一个加速因子beta取值1.2将修正项进行幂次放大从而加快收敛。同时为了防止振荡当两次迭代的谱变化超过10%时自动降低beta至0.9。对混合核素谱含40K、238U、232Th的解谱测试显示Boosted-Gold算法在35次迭代后达到稳定各核素活度相对误差小于3.5%而标准Gold需要82次迭代。3蒙特卡洛响应矩阵与核素定量分析快速方法利用Geant4软件包建立NaI(Tl)探测器的精确几何模型晶体尺寸φ76mm×76mm模拟从30keV到3MeV范围内间隔5keV的单能光子入射记录每个能量的响应谱构成响应矩阵600×1024。针对测量谱的解析采用一种基于最大后验概率MAP的迭代算法。该算法在目标函数中加入先验信息项假设核素活度分布服从拉普拉斯分布从而引入稀疏约束。使用非线性共轭梯度法求解MAP最优化问题每次迭代更新解向量。在实际应用中对于含有5种未知核素的混合样品MAP算法从测量谱中直接反演出各核素活度无需单独寻峰和本底扣除分析速度比传统逆矩阵法快3倍且对统计涨落的鲁棒性更好。在一组低计数率能谱总计数2000上MAP算法的核素识别率仍能达到92%而常规峰搜索法只有71%。import numpy as np from scipy.special import erfc from scipy.optimize import minimize import math class AdaptiveSNIP: def __init__(self, data): self.data data.astype(np.float64) def remove_background(self): y self.data.copy() m 1 while m len(y)//2: width int(2 * m) # 简化自适应实际与FWHM相关 for i in range(width, len(y)-width): if y[i] (y[i-width] y[iwidth])/2: y[i] (y[i-width] y[iwidth])/2 m 1 return y class BoostedGold: def __init__(self, response_matrix, beta1.2): self.R response_matrix self.beta beta def deconvolve(self, measured, max_iter100): spec measured.astype(np.float64) spec / spec.sum() for it in range(max_iter): est self.R spec ratio measured / (est 1e-10) correction self.R.T ratio spec spec * (correction ** self.beta) spec spec / spec.sum() if it 5 and np.max(np.abs(correction-1)) 0.01: break return spec def gaussian_response_matrix(energy_edges, fwhm_func, n_channels1024): E_centers (energy_edges[:-1] energy_edges[1:]) / 2 matrix np.zeros((len(E_centers), n_channels)) for i, E in enumerate(E_centers): fwhm fwhm_func(E) sigma fwhm / (2.355) for ch in range(n_channels): energy_ch ch * (E_centers[-1]-E_centers[0])/n_channels matrix[i, ch] math.exp(-0.5*((energy_ch-E)/sigma)**2) matrix[i] / matrix[i].sum() return matrix class MAPDeconv: def __init__(self, response, lambda_reg0.01): self.R response self.lam lambda_reg def negative_log_posterior(self, x, y): ax self.R x likelihood -np.sum(y * np.log(ax 1e-8) - ax) prior self.lam * np.sum(np.abs(x)) return likelihood prior def solve(self, y, x0None): if x0 is None: x0 np.ones(self.R.shape[1]) res minimize(self.negative_log_posterior, x0, args(y,), methodL-BFGS-B, bounds[(0,None)]*len(x0)) return res.x if __name__ __main__: test_spec np.random.poisson(100, 1024) 50 snip AdaptiveSNIP(test_spec) bg snip.remove_background() print(SNIP background removed, mean:, np.mean(bg)) # 高斯响应矩阵示例 def fwhm_func(E): return 0.05 * np.sqrt(E) 0.02 edges np.linspace(0, 2000, 200) R_gauss gaussian_response_matrix(edges, fwhm_func, 1024) print(Response matrix shape:, R_gauss.shape) # Boosted Gold measured np.random.rand(1024) bgold BoostedGold(R_gauss[:100, :], 1.2) deconv bgold.deconvolve(measured) print(Deconvolved sum:, deconv.sum()) # MAP map_solver MAPDeconv(R_gauss[:200, :200], 0.005) y_sim R_gauss[:200, :200] np.random.rand(200) x_est map_solver.solve(y_sim) print(MAP reconstruction error:, np.linalg.norm(x_est - np.ones(200)))
低能量分辨率γ能谱数据解析方法解析【附数据】
✨ 长期致力于低分辨率γ能谱、高分辨解析、峰边界确定、自适应半高宽、响应矩阵研究工作擅长数据搜集与处理、建模仿真、程序编写、仿真设计。✅ 专业定制毕设、代码✅如需沟通交流点击《获取方式》1自适应半高宽匹配的SNIP本底扣除与反卷积联合寻峰针对NaI(Tl)探测器能谱中低能区本底畸变和重叠峰严重的问题提出一种自适应SNIP本底扣除算法。该算法根据能谱的道址和计数率动态选择变换窗宽m窗宽取值与当前道址处的半高宽FWHM成正比m 2 * FWHM / 道宽。在标准SNIP迭代中对每个道址进行多次削峰操作直到本底趋于平滑。本底扣除后采用反卷积联合寻峰策略首先使用Gold反卷积算法对谱线进行解卷积将重叠峰的FWHM压缩至原始的一半然后在一阶导数零交叉点处识别峰位。对于严重重叠峰峰间距小于1.2倍FWHM引入峰边界约束的高斯混合模型用EM算法估计各子峰参数。在一块238U矿石样品实测谱中该方法成功分离了185.7keV和186.2keV的重叠双峰峰位偏差小于0.3keV。2高斯响应矩阵与Boosted-Gold迭代解谱基于NaI(Tl)探测器对γ射线的响应近似高斯函数这一物理特性构建了一个高斯响应矩阵。矩阵元素R_ij表示能量为Ej的单能γ射线在第i道址产生的计数其中高斯函数的FWHM按照能量刻度公式FWHM a*sqrt(E)b 计算a和b通过241Am、137Cs、60Co等标准源标定。矩阵维度为1024×1024。解谱时采用改进的Boosted-Gold算法与传统Gold不同该算法在每次迭代中引入一个加速因子beta取值1.2将修正项进行幂次放大从而加快收敛。同时为了防止振荡当两次迭代的谱变化超过10%时自动降低beta至0.9。对混合核素谱含40K、238U、232Th的解谱测试显示Boosted-Gold算法在35次迭代后达到稳定各核素活度相对误差小于3.5%而标准Gold需要82次迭代。3蒙特卡洛响应矩阵与核素定量分析快速方法利用Geant4软件包建立NaI(Tl)探测器的精确几何模型晶体尺寸φ76mm×76mm模拟从30keV到3MeV范围内间隔5keV的单能光子入射记录每个能量的响应谱构成响应矩阵600×1024。针对测量谱的解析采用一种基于最大后验概率MAP的迭代算法。该算法在目标函数中加入先验信息项假设核素活度分布服从拉普拉斯分布从而引入稀疏约束。使用非线性共轭梯度法求解MAP最优化问题每次迭代更新解向量。在实际应用中对于含有5种未知核素的混合样品MAP算法从测量谱中直接反演出各核素活度无需单独寻峰和本底扣除分析速度比传统逆矩阵法快3倍且对统计涨落的鲁棒性更好。在一组低计数率能谱总计数2000上MAP算法的核素识别率仍能达到92%而常规峰搜索法只有71%。import numpy as np from scipy.special import erfc from scipy.optimize import minimize import math class AdaptiveSNIP: def __init__(self, data): self.data data.astype(np.float64) def remove_background(self): y self.data.copy() m 1 while m len(y)//2: width int(2 * m) # 简化自适应实际与FWHM相关 for i in range(width, len(y)-width): if y[i] (y[i-width] y[iwidth])/2: y[i] (y[i-width] y[iwidth])/2 m 1 return y class BoostedGold: def __init__(self, response_matrix, beta1.2): self.R response_matrix self.beta beta def deconvolve(self, measured, max_iter100): spec measured.astype(np.float64) spec / spec.sum() for it in range(max_iter): est self.R spec ratio measured / (est 1e-10) correction self.R.T ratio spec spec * (correction ** self.beta) spec spec / spec.sum() if it 5 and np.max(np.abs(correction-1)) 0.01: break return spec def gaussian_response_matrix(energy_edges, fwhm_func, n_channels1024): E_centers (energy_edges[:-1] energy_edges[1:]) / 2 matrix np.zeros((len(E_centers), n_channels)) for i, E in enumerate(E_centers): fwhm fwhm_func(E) sigma fwhm / (2.355) for ch in range(n_channels): energy_ch ch * (E_centers[-1]-E_centers[0])/n_channels matrix[i, ch] math.exp(-0.5*((energy_ch-E)/sigma)**2) matrix[i] / matrix[i].sum() return matrix class MAPDeconv: def __init__(self, response, lambda_reg0.01): self.R response self.lam lambda_reg def negative_log_posterior(self, x, y): ax self.R x likelihood -np.sum(y * np.log(ax 1e-8) - ax) prior self.lam * np.sum(np.abs(x)) return likelihood prior def solve(self, y, x0None): if x0 is None: x0 np.ones(self.R.shape[1]) res minimize(self.negative_log_posterior, x0, args(y,), methodL-BFGS-B, bounds[(0,None)]*len(x0)) return res.x if __name__ __main__: test_spec np.random.poisson(100, 1024) 50 snip AdaptiveSNIP(test_spec) bg snip.remove_background() print(SNIP background removed, mean:, np.mean(bg)) # 高斯响应矩阵示例 def fwhm_func(E): return 0.05 * np.sqrt(E) 0.02 edges np.linspace(0, 2000, 200) R_gauss gaussian_response_matrix(edges, fwhm_func, 1024) print(Response matrix shape:, R_gauss.shape) # Boosted Gold measured np.random.rand(1024) bgold BoostedGold(R_gauss[:100, :], 1.2) deconv bgold.deconvolve(measured) print(Deconvolved sum:, deconv.sum()) # MAP map_solver MAPDeconv(R_gauss[:200, :200], 0.005) y_sim R_gauss[:200, :200] np.random.rand(200) x_est map_solver.solve(y_sim) print(MAP reconstruction error:, np.linalg.norm(x_est - np.ones(200)))
