核数据处理课程设计报告正式稿

核数据处理课程设计题目：γ能谱谱数据分解方法研究

英文题目：γ

Spectroscopy spectral data

Decomposition Method

学院：核工程与地球物理学院学生姓名：刘煜锴专业年级：核工程与核技术0 9级学号：0 9 0 2 1 1 3 5 指导教师：吴和喜老师

二零一二年六月

γ能谱谱数据分解方法研究

摘要

由多道脉冲分析器获取的谱数据需要以一定的数学方法进行处理才能得到实验要求的最终结果。能谱的数据处理大致可以分为两个步骤。首先进行峰分析，即由能谱数据中找到全部有意义的峰，并计算出扣除本底之后每个峰的净面积。第二步是放射性核素的活度或样品中元素浓度的计算，即由峰位所对应的能量识别出被测样品中含有哪些放射性核素或被激发的元素，并且由峰的净面积计算出放射性核素的活度或元素在样品中的浓度。

采用不同的物理实验方法，使用不同的探测器时，能谱的数据处理方法也有所不同。在本课程设计报告中采用能谱数据处理中经常用到的峰分析方法，包括谱数据的平滑处理，本底扣除、寻峰、峰净面积计算和函数拟合法解谱。然后，以γ谱分析为例，讨论基于计算机的多道脉冲分析系统中的谱在matlab软件平台的数据分解。

γ能谱，解谱，平滑处理，寻峰寻峰

关键词：

γ

Spectroscopy spectral data Decomposition

Method

Abstract

Spectral data obtained by multi-channel pulse analyzer need to be the end result of certain mathematical methods to deal with to get the experimental requirements. The energy spectrum of data processing can be divided into two steps. First peak analysis, all the significant peaks from the spectral data, and calculate the net area of each peak after deducting the background. The second step is Radionuclides or the calculation of element concentrations in the samples, energy identified by the peak position corresponding to the tested samples containing radionuclides or excited elements and calculate the net area of the peak radionuclides or elements in the sample concentration.

Different physical experiments, using a different detector, the energy spectrum of data processing methods are different. This course is designed report using the peak analysis method often used in the spectral data processing, including spectral data smoothing, background subtraction, peak search, peak net area calculations and functions fitting solution spectrum. Then, spectral analysis, for example, the discussion of the data decomposition in matlab software platform based on the spectrum in the multi-channel analyzer computer.

Key Words：γSpectroscopy, deconvolution, smoothing, peak searching to find the peak

目录

目录 (4)

1前言 (5)

1.1 γ能谱分析的目的、意义 (5)

1.2 γ能谱分析的主要难题 (5)

1.3 γ能谱解谱的研究简介 (6)

1.4 本次课程设计γ能谱分析的研究内容 ............................................................................................. 6 2 γ

谱数据的平滑 ....................................................................................................................................... 8 2.1 重心法 (8)

2.2 多项式最小二乘移动平滑法 (8)

2.3傅里叶变换法 (9)

2.4 小波变换法 (10)

3 峰位与峰区的确定 (10)

4 峰面积方法简述 (11)

5 峰相关参数简介 (12)

6 能量刻度 (13)

7 程序的测试 (14)

7.1实验结果 (14)

7.2数据分析 (17)

7.3说明 (17)

8 心得体会 (18)

1前言

1.1 γ能谱分析的目的、意义

随着社会的发展，人民生活水平的提高，核能应用，环境辐射污染等方面越来越引起人们的重视，因此核辐射监测技术在经济发展，国防建设等方面起着越来越重要的作用，人们更加关注放射性污染对人体所造成的危害。

γ能谱测量作为一种重要的核地球物理方法，是解决地球科学、环境科学等有关

γ能谱仪作为核辐射测量中最常用的数据采集，分析仪器，其通过问题的主要手段之一。

γ射线能量的测量可识别发射γ射线源的核素种类，对辐射强度测量能够获得本γ射对

γ辐射总量。γ能谱仪常被用十家居及环境的辐射监测，铀矿线的核素含量或活度，及

勘探，地质填图，油气勘测，寻找各种金属和非金属矿产，食品卫生检验等多个领域，因

γ能谱测量技术、分析技术有着广阔的发展前景。

此

γ能谱分析任务主要包括γ能谱的测量和数据的采集，谱数据的解谱分析等。主

γ谱数据平滑、谱寻峰、核素识别、峰区边界确定、本底扣除、全能要解谱分析内容包括

峰函数拟合、净峰面积求解、核素含量计算及其活度分析等。

1.2 γ能谱分析的主要难题

γ能谱仪。其中目前用十放射性核素测量的Y能谱仪主要有两种:NaI(Tl)和HPGe

γ能谱仪的能量分辨率高(对于60CO的1332keV全能峰可以达到2. OkeV左右)，能HPGe

γ能谱仪，其解谱技术相对比较简单，在复杂核素组成的量线性好，被认为是当今最好的

样品中也可以得到较好的分析结果，但是由于HPGe探头需要在液氮冷却的环境下进行测

γ能量，目_价格也比较昂贵，使用维护不便，其使用的环境受到一定的限制。而NaI(Tl)

谱仪由十使用维护方便，探测效率高(约是HPGe的几十倍)，仪器成本低等多方面占有优势，因此在放射监测中仍被广泛的应用，但Na工(T1>能谱仪能量分辨率相对较低(对于137CS 的661keV一般为7%左右)，对于核素种类较多，谱线较复杂的样品进行解谱分析时却无法胜任，因此必需对其Y能谱进行较为复杂的算法分析才能得到比较准确的结果。

γ能谱仪进行放射性γ射线谱测量时，因其能量分辨率较低，如果所测物 NaI(T1)

γ谱质中所含核素种类多、含量低、介质复杂，在高本底等复杂背景的测量环境下，所测

线会相当复杂，如果解谱算法精度不高，会影响最终分析结果的不确定性，主要的影响因

γ能谱仪仪器性能的影响，γ能谱仪的测量统计涨落干扰，γ能谱仪的素有:NaI(T1)

本底干扰。

1.3 γ能谱解谱的研究简介

γ能谱分析任务主要包括γ谱数据平滑、寻峰、核素识别、峰区边界确定、本底扣全能峰拟合，净峰面积求解和核素含量计算及其活度分析等。

γ谱数据平滑方法主要有:傅里叶变换、几点平滑法、指数平滑法、马尔可夫平滑法、

γ能谱展开，同时利用能域和频域的信息进行平小波变换(即在联合能量一频率空间内把

滑处理，被认为是一种非常理想的能谱平滑处理方法)、数字滤波器，重心法等。目前实

际应用中常用几点平滑法，马尔可夫平滑法。尽管上述平滑方法各有优势，但是大多都需

要谱分析者根据实际经验设定平滑参数或选择平滑函数等，才能得到较为准确的平滑效

果，且并不是任何情况下总是有效的，特别是在能谱数据中含大量密集干扰背景的环境下，

平滑效果并不理想，可能会出现假峰或遗漏真峰等。

γ能谱寻峰方法主要有:简单比较寻峰法、高斯乘积函数找峰法、五点极值法、道数寻峰法、导数法(二阶导数寻峰精度高，一阶导数确定峰区边界准确)、对称零面积法、协方差法(精度高，速度慢)，斜宽寻峰法(速度快)等，各有优缺点。但这些寻峰方法却没有

γ能谱的寻峰问题一直是一种方法能够同时达到寻峰速度快，寻峰精度又高的要求。复杂

γ能个比较难以解决的问题，因此应充分考虑各种寻峰法的优点，联合使用几种寻峰法对

谱进行联合寻峰是一种比较理想的解决办法。

γ能谱本底扣除方法主要有:直线本底、阶梯/阶跃本底，抛物线本底等二种传统本底扣除方法，目前大部分学者认为阶跃本底是一种比较准确的扣除本底方法。近年来也有部分学者利用傅立叶变换，小波变换等方法在扣除本底方面上进行相关的研究，并已取得一定的研究成果。国外部分学者在阶跃本底的基础上采用SNIP (Sensitive

Nonlinear:Iterative Peak Algorithms)本底扣除方法、，被认为是目前比较理想的本底扣除方法。但不管采用哪种本底扣除，只有充分利用全能峰尽可能多的有效计数，才能认为是比较理想的本底扣除方式，才可能得到比较好的本底扣除效果。在本底扣除时，全能峰峰区边界的选择也会影响到本底扣除的质量，因此精确定位峰区边界也成为一个必须考虑的因素。

1.4 本次课程设计γ能谱分析的研究内容

γ能谱测量、γ能谱的物理形成原理、目前常用解谱算法分析本次课程设计通过对

γ能谱仪的计数统计涨落问题，重心法、多项式最小二乘移动平滑做了较好的探究，针对

傅里叶变换法、小波变换法来对谱数据进行平滑处理;然后采用弱峰检测能力较好、重叠

峰分辨能力好的对称零面积寻峰法进行寻峰;并通过了二阶插值多项式拟合来准确的计算峰位，还确定了全能峰峰区边界道址知道了全能峰的峰位的道址及其对应的全能峰的能量后，求出了能量刻度方程曲线二次多项式y=kx+b。然后利用瓦森峰面积法来计算各个寻得的峰的峰面积。

最后对所求的全能峰净计数，全能峰净面积求解等几个关键问题进行了统一分析验证，以便能达到较为可靠准确的解谱质量，并与标准参考数据进行对比验证。

γ能谱解谱方法研究的主要内容框架（程序流程图）如图一所示：本课程设计的

γ能谱解谱方法研究内容

图一

2 γ谱数据的平滑

由于核衰变和γ能谱仪固有的统计涨落、电子电路系统的噪声，γ射线散射的影响

γ能谱中，每个道址上的计数可能与理论期望值有着很大的偏差，故会形成一等，在实测

γ谱形。这种涨落不仅会影响全能峰峰位的确定和净峰面积的计种带有计数统计涨落的

算，而且也会掩盖掉弱峰，漏判可探测的核素，或是出现假峰，反而会识别出并不存在的

γ能谱的定性定量分析产生不可避免的分析误差。为了减少γ谱数据统计涨落核素，对

γ谱峰全部的重要特征，以便对γ能谱进行更为可靠的分析，的影响，而又尽可能保留

γ能谱原始数据进行平滑或去噪处理，这就是γ能谱平滑处理技术。

必须先对所测

γ谱数据是按整数道址离散存储的，由于相邻各道的计数存在着一定的数学关联，所以利用数学统计方法逐道对谱数据进行修正，是可以消除大部分统计涨落的影响，且仍会保留原始谱的重要特征。一些经典的谱平滑算法，例如，傅里叶变换法、离散函数褶积滑动变换法、多项式最小二乘移动平滑法（也叫几点平滑法）、马尔可夫平滑法、小波变换

γ能谱展开，利用能谱在能域和频域的信息进行平滑处法(即在联合能量一频率空间内把

理，部分人认为这是一种非常理想的能谱平滑处理方法)、数字滤波器，重心法等。本课程设计采用了重心法，多项式最小二乘移动平滑法、傅里叶变换法、小波变换法这四种平滑方法进行平滑。

2.1 重心法

所谓重心法就是假设光滑后的数据是原来数据的重心

重心法二点光滑公式:

重心法五点光滑公式:

2.2 多项式最小二乘移动平滑法

1964年A. Savitzky和J.Egolay提出了一个用于谱数据平滑处理的滤波器响应函数。其基本思想是，当求平滑之后谱的第m点数据时，先在原始谱数据第m点的左、右各取K

个数据点，形成一个共有2K+1个数据点的窗口。在这个窗口中用多项式拟合原始谱数据，则拟合多项式在m 点的值就是平滑后的谱在m 点的值。当m 值沿谱数据移动时，就可以得到整个平滑后的谱数据。这种方法称为最小二乘移动平滑法，或最小平方曲线拟合平滑法。 滤波之后的谱数据为：

∑+-=+=K K

j j m j m y g y

∑+-==K

K j j g 1

当平滑窗口选为5点时（K=2），5点平滑公式为：

)31217123(35

12112++---+++-=m m m m m m y y y y y y 2.3傅里叶变换法

该方法首先进行傅里叶变换，将一个时域(能域)函数变换到频域，然后在频域中滤波，最后再进行傅里叶逆变换，将频域函数变换到时域(能域)中。

令y （E)表示原始谱，~

y (E)表示光滑谱，s(E)表示噪音。即:

使用大写字母表示傅里叶变换后的函数，能量频率分布函数为:

离散形式为:

,其中K=0,1,......N-1, N 为总

道数。

频率谱分布大致与能域谱的相关长度成反比，~

y (E)一般有多道相关性，s(E)的相关长度为一道。所以，~Y (E)的频率范围大概为1-2弧度/道，S(ω)大概为π弧度/道。原始谱

是噪声谱和光滑谱的叠加，光滑谱主要集中于低频范围，可以采用一个低通滤波函数F(ω)，即：

但是它具有边缘效应。采用高斯滤波，不产生边缘效应，即：

傅里叶逆变换后，为：

其中σ为高斯宽度，σ=355.2H

,H 为半宽度。

2.4 小波变换法

令函数()n s 表示γ

能谱信号()n f 被噪声污染后的能谱分布，其中变量n 为道址。

则基本噪声模型为：

其中()n e 为信号噪声，σ为信号强度，则γ能谱的小波降噪方法为：

1)选定一种小波，对原始能谱()n s 连续做需要层次的小波分解;

2)在小波域中对分解得到的系数进行阀值处理，包括软阈值法和硬阈值法。

3)将处理之后的系数进行逆离散小波变换，重构出降噪后的能谱()n f 。

3 峰位与峰区的确定

本课程设计中采用的寻峰方法为对称零面积寻峰法。

对称零面积变换法其实就是用面积为零的窗函数与实验能谱数据进行褶积变换，并且要求窗函数是对称函数。对线性基底的褶积变换将为零，只有存在峰的地方不为零。

进行谱变换时使用的数字滤波器的冲击函数(零面积函数)Tj 满足下列关系式：

?????==∑-=-K K j j j j 0

T T T （j=－K ，……，＋K ）

即冲击函数围绕j ＝0点左右对称，且对横轴所包围的面积为零，因而称为对称零面积对合法。如果在峰区范围内本底谱是常数或按直线分布，则在变换之后的谱中本底的贡献为零。变换之后谱的形状完全反映了峰形的变化。下一步使用峰高统计判定条件来确定峰位。沿变换之后的谱进行检索，找出局部极大点，若极大值超过其均方根误差若干倍时，则认为找到了一个峰，该局部极大点所对应的道址就是峰位。沿道址在谱中进行检索，找出一系列的局部极大值，并进行判定就可以找到一系列的峰。上述的峰高判定条件可以写为：

TRH

21K K j j m 2j K K j j m j *m *m m ≥???

? ??=?=∑∑-=+-=+y T y T y y R 其中TRH 为寻峰阈值，是预先给定的常数。当寻峰阈值取得大时，能有效地剔除由统计涨落造成的假峰。但也可能同时漏掉了弱的真峰。当设定较小的寻峰阈值时，能找到全部的强峰和弱峰，但也可能会把由统计涨落造成的假峰当成真峰。因此，寻峰阈值的大小应当根据谱数据的具体情况和物理实验任务来确定。

当然对称零面积函数还可以选择矩形函数。本课程设计的程序中就是采用的矩形函数。它的思想与上面相同。

4 峰面积方法简述

本课程设计中，选用了Wasson峰面积法进行峰面积计算。因为结合了上述两种方法的优点：取基线时采用总峰面积法，确定峰区时采用科沃尔法，这样，可以减少科沃尔法基线过高的影响，同时也避免了总峰面积法确定的峰区边界对误差贡献过高的影响。当然为了便于比较，还采用了与总峰面积法（TPA）进行峰面积计算。

5 峰相关参数简介

，实际应用中，峰址：通过寻峰算法，对各道的计数进行处理，找出可能的峰道址i

0,

一般要给出峰位置处的道数，相应的能量值。寻峰算法一般采用高斯对称零面积卷积法当

变换谱与它的标准偏差SS(i）出现的正值大于给定的灵敏因子时，才可能存在峰，峰址为

平滑后峰区内最大计数道址处。

峰左右边界：是确定峰面积的重要参数之一。由于全能峰大概成高斯分布，在峰的左右边界处，除去本底的计数很少，而相对统计涨落又很大。有些全能峰的边界拉的很长，正确截取边界值直接影响后期活度计算的结果精度。采用微分法可以确定峰左右边界，这里我使用对称零面积法，最终确定峰的左右边界值 L 和R。

净峰面积：为累积峰面积计数的参数之一，通过净峰面积可以计算出放射性同位素的活度。净峰面积就是去除本底计数后峰区计数总和。在计算净峰面积时，本底去除方法至关重要，最简单的是线性本底去除，令 s (i )为光滑后的能谱图的第i道计数值，则线性本底公式为：

6 能量刻度

能量刻度就是假定能量与道址成线性关系，用一系列已知能量的峰及其相应的道址，用线性拟合的方法得到一个数学表达式。

（1）能量刻度方程

设其表达式为：E=ai+b；

式中：E：能量值；i：道址； a和b为未知系数；

求出a、b后，代入上式就能得出能量和道址的函数关系式，这样我们就能知道每一道对应的射线的能量值，通过相应的能量峰的峰位前后位置比较，可以判断谱线是否产生了漂移，这关系到后面的谱线分析。

（2）能量刻度需要注意的问题

为了准确地进行能量刻度和确定未知能量，需注意以下几个问题：

①要求放射源的能量是精确的已知值；

②要求峰位确定的准确；

③要求考虑非线性；

④要求谱仪的稳定性好。

（3）本课程设计的能量刻度方法

至少确定两个精确峰位的能量和其对应的道址，并输入到matlab的命令窗中。为了保证刻度的精确度，所选的峰位的能量值在50~4000KeV之间。输入道址与对应能量后，在本课程设计中，采用的方法是最小二乘法线性拟合。

7 程序的测试

7.1实验结果

本次课程设计采用的谱数据为iaea-1995文件夹下iaearfnw TSTSPEC 里面的数据。首先来看看CALIB.ASC。READ_ME.TXT中说明了这个谱数据包含的部分峰得峰位与对应能量如下：

Channel Energy (keV)

301 122.06

1281 511.00

1661 661.66

2097 834.84

2951 1173.24

3207 1274.54

3353 1332.50

运行程序，其中参数选择为：选择傅里叶变换法平滑，高斯滤波器A=1,FWHM=4，对称零面积法的参数是K=2，H=3,b=1得出的结果为：

各峰的道址如下：

301 1281 1661 2098 2951 3207 3353

峰所在道址对应的计数:

1.154515e+004 8.699714e+003 6.831708e+003

2.917624e+003

3.426126e+003

2.729146e+003 2.872052e+003

各峰的准确道址如下：

3.005700e+002 1.280555e+003 1.660610e+003 2.097531e+003 2.951262e+003 3.206988e+003 3.353291e+003

各峰的左边界道址如下：

296 1274 1656 2093 2946 3202 3348

各峰的右边界道址如下：

305 1287 1665 2102 2956 3212 3358

能量刻度曲线系数为：0.396595 2.787151

即该曲线方程为：E=0.396474*i+3.124935

采用瓦森峰面积法

各个峰的面积为:

3.768820e+004 3.198056e+004 2.349197e+004 9.690379e+003 1.237195e+004

9.915058e+003 1.036585e+004

如果采用总峰面积法的话

下面采用总峰面积法

各个峰的面积为:

4.395662e+004 4.981342e+004 2.813315e+004 1.175224e+004 1.561551e+004

1.255346e+004 1.320438e+004

输出的平滑前与平滑后的谱的对比如下图7.1：

010002000300040005000

6000700080009000

4道址能量平滑前与平滑后的谱

图7.1平滑前与平滑后的谱

刻度好的能量刻度曲线如下图7.2所示：

010002000300040005000

6000700080009000

0500

1000

1500

20002500

3000

3500

道址能量k e v 刻度好的能量刻度方程曲线

图7.2 刻度好的能量刻度曲线

寻峰之后的道址与计数的关系图如图7.3所示：

010002000300040005000

6000700080009000

-20000

2000

4000

6000

8000

10000

12000

道址计数

寻峰之后的道址与计数的关系

图

图7.3寻峰之后的道址与计数的关系图

7.2数据分析

通过实验数据中的峰位与对应能量与READ_ME.TXT中说明的部分峰峰位与对应能量对比可以发现，数据吻合的非常好。只有第四的峰的道址相差1道。而且经能量刻度后，把峰的道址代入能量刻度曲线，算得的对应得能量也吻合的很好。

对于峰的面积计算，通过瓦森峰面积法与总峰面积法的比较可以得出：选用了Wasson 峰面积法进行峰面积计算。因为结合了总峰面积法与科沃尔峰面积法两种方法的优点：取基线时采用总峰面积法，确定峰区时采用科沃尔法，这样，可以减少科沃尔法基线过高的影响，同时也避免了总峰面积法确定的峰区边界对误差贡献过高的影响。虽然总峰面积法对一给定的全能峰可以获得最大的总计数，但是由于边界道对峰面积贡献不大，但却使得误差显著增加。

7.3说明

由于READ_ME.TXT中说明了这个谱数据包含的部分峰得峰位与对应能量只是针对比较强峰而言，还可以使用重心法平滑、选择多项式最小二乘移动平滑法、小波变换平滑谱数据，改变参数就可以寻得相对比较弱的峰。例如：选择最小二乘移动平滑法，窗口宽度为5，或选择小波变换法平滑，对称零面积法的参数是K=2，H=3,b=1。得出的各峰的道址如下：

181 301 337 1281 1661 2098 2951 3207 3353

这样的话就比以上一次的参数运行得出的峰的个数增加了两个（相对而言的弱峰），由于不同的参数，不同的谱数据文件得出的数据，图形各不相同，在这里就不一一列举了。

8 心得体会

本次的课程设计虽然结束了，但是努力学习并没有结束。这是对我所学的课本知识的测验，也是对是对编程思想的强化。

课程设计是对所学知识的综合理解与应用，它不仅要求我们对matlab软件特别熟悉，还让我懂得发现问题，分析问题，并解决问题可以提高自己的能力。同时通过分析解决问题，加深对所学知识的理解与掌握。通过这次课程设计我学到了很多的东西，在课设期间我几乎每天都认真思考更好的办法，将自己的想法通过编程实现。很多时候我坚持认为自己的想法没有问题，但就是怎么也不能达到预期的结果。当经过很多次的修改、重新运行，终于成功的那一刻，内心的成就感和自豪感是那么的强烈。

经过这次课程设计，使我懂得了：采用不同的物理实验方法，使用不同的探测器时，能谱的数据处理方法也有所不同。在本课程设计报告中采用能谱数据处理中经常用到的峰分析方法，包括谱数据的平滑处理，本底扣除、寻峰、峰净面积计算解谱。然后，以γ谱分析为例，讨论基于计算机的多道脉冲分析系统中的谱在matlab软件平台的数据分解。

γ能谱谱数据分解方法。

深化理解了