(完整word版)多道γ能谱分析软件中寻峰算法比较总结
黑糊糊-中学生500字作文
自动寻峰
由于谱结构的复杂和统计涨落的影响,从谱中正确地找到全部存在
的峰是比较困难的。
尤其是找到位于很高本底上的弱峰,分辨出相互靠得很近的重峰更为困难。
谱分析对寻峰方法的基本要求如下:
(1)
比较高的重峰分辨能力。能确定相互距离很近的峰的峰位。
(2)
能识别弱峰,特别是位于高本底上的弱峰。
(3) 假峰出现的几率要小。
(4)
不仅能计算出峰位的整数道址,还能计算出峰位的精确值,某些情况下要求峰位的误
差小于0.2道。
很多作者对寻峰方法进行了研究,提出了很多有效的寻峰方法。
目的 :
判断有没有峰存在
确定峰位(高斯分布的数学期望),以便把峰位对应的道址,转换成能量
确定峰边界——为计算峰面积服务(峰边界道的确定,直接影响峰面积的计算)
分为两个步骤:谱变换和峰判定
要求:支持手动自动寻峰,参数输入,同时计算并显示峰半高
宽、精确峰位、峰宽等信
息,能够区分康普顿边沿和假峰
感兴区内寻峰
人工设置感兴趣大小,然后在感兴区内采用简单方法寻峰
重点研究:对感兴区内的弱峰寻峰、重峰的分解
对于一个单峰区,当峰形在峰位两侧比较对称
时,可以由峰的FWHM计算峰区的左、右
边界道址。峰区的宽度取为3FWHM,FWHM的值可以根
据峰位m
p
由测量系统的FWHM
刻度公式计算。由于峰形对称,左、右边界道和峰位的距离都是1.5FWHNM。
m
L
INT(m
p
1.5FWHM0.5)
m
R
I
NT(m
p
1.5FWHM0.5)
式中m
p
是峰位,INT的含义是取整数。
对于存在有低能尾部的峰,其峰形函数描述(参见图)。
y
m
HEXP[
(mm
p
)
2
2
2
]
,m≥mp-J
,m≤mp-J
y
m
HEXP[J(2m2m
p
J)2
2
]
式中H为峰高,mp为峰位,
是高斯函数的标准偏差
,J为接点的道址和峰位之间的距
离。在峰位的左侧,有一个接点,其道址为mp-J。在接点的右侧,
峰函数是高斯函数。
在接点的左侧,峰函数用指数曲线来描述。这时峰区的左、右边界道址为
m
L
INT(m
p
1.12FWHM
2
J0.5J
0.5)
m
R
INT(m
p
1.5FWHM0.5)
带有低能尾部的峰函数的图形
全谱自动寻峰
基于核素库法:能量刻度完成后,根据
核素库中的能量计算对应的道址,在各个道址附
近(左右10道附近)采用简单的寻峰方法(导数法)
方法:
根据仪器选择开发
IF函数法简单比较法(适于寻找强单峰,速度快)
满足条件:
data
im
data
i
kdata
i
data
im
可认为有峰存在
然后在data i-m至data i+m中找最大值,对应的道值即为峰位
k:找峰阈值,根据高斯分布,一般k取值1—1.5
常用5点、7点极大值法(m取2,3)
判定峰是否有意义
一般,用R=N0
Nb≥ R0确定峰是否有意义
R为峰谷比, R0为设定值 (经验值)
N0为净峰幅度与基底之和
Nb为基底计数
int CMmcaView::SearPeakCompare(int Beginch,
int Endch, int m, float k)
高斯乘积函数找峰法(可靠性差,不建议采用)
描述谱峰形状的函数主要是高斯函数
G(i)
A
exp(ii
0
)
2
2
2
则由相邻的数据点定义
2
一个新的函数(第一高斯乘积
函数,只与
FWHM2.3556
有关):
P
m
(i
)
G(i)G(im1)11.092m
exp()m2
2G(i2)G(im)H
m是步长(用道表示),是高斯乘积函数的阶数,则Pm(i)称为第
m阶高斯乘积函数。找
峰的灵敏度与m有关,随m的增加灵敏度提高。
为避免基线参数的影响,最好扣除本底后,再应用高斯乘积函数找峰。
考虑统计涨落的影响,把判断无峰存在的1变为一个“单位带”。即峰的判断为:
1ky
i
P
m
(i)
(1ky
i
)
无峰
有峰
k3
峰位的确
定:由Pm(i)过1的两点求平均来确定;峰边界的确定:“单位带”下限的两个最
端点;半高宽的确
定:函数Pm(i)在“1”上的截距;组合峰的确定:在乘积函数的两个峰之
间没有处于“带内”的乘积函数值
导数法(一阶、二阶、三阶)
1
y
N
m
'
i
jm
C
m
jy
ij
Nm为规范化常数,Cj平滑的变换系数。
3次多项式5点光滑一阶导数公式:(可以采用)
y
i
'
1
(y
i2
8y
i1
8y
i1
y<
br>i2
)
峰位确定:一阶导数值由正变负=0处;峰边界确定:一阶导数
12<
br>由负变正=0处
CalculateDifferential(0, size, m,
differ);
for (int j = m; j <= size-m; j++)
{
for(int i=1;i<=m;i++)
{
if(differ[j-i])>0&&differ[j-i]>maxtemp)
{maxtemp=differ[j-i]; nmax=j-i;}
if(differ[j+i])<0&&differ[j+i]
if
((nmin-nmax)>0.8*fwhm && (nmin-nmax)<3*fwhm)
FWHM参数根据仪器能量分辨率可人工确定,fwhm~20
peakposition[p++]=j+0.5;保持峰位对应的道址
}
5点光滑二阶导数公式(软件中推荐采用)
1
y
i
''
(2y
i2
y
i1
2y
i
y
i12y
i2
)
7
7点二阶导数
(5*(coun
tsdata[j-3]+countsdata[j+3])-3*(countsdata[j-1]+co
untsdata[j+1])-4*countsdata[j])42;
y
i
'
1
(22.0y
i3
67.0y
i2
58.0y
i1
58.0y
i1
67.0y
i2
22.0y
i3
)
252.0
软件中推荐采用11点以上的公式
峰位确定:二阶导数最小值对应的道址;峰边界确定:二阶导数正极大值点
for (int
j = m; j <= size-m; j++)m~30
{
int
maxtemp=-0.5,mintemp=-0.5;
If(differ[j]<-0.05)
for(int i=1;i<=m;i++)
{
if(differ[j-i]>maxtemp)
{maxtemp=differ[j-i]; nmax=j-i;}
if(differ[j+i]>mintemp) {mintemp=differ[j+i];
nmin=j+i;}
}
if ((nmin-nmax)>0.8*fwhm &&
(nmin-nmax)<3*fwhm)
FWHM参数根据仪器能量分辨率可人工确定,fwhm~20
peakposition[p++]=j+0.5;保持峰位对应的道址
}
D<
br>i
n
i1
2n
i
n
i
jk
cn
j1
1
0
jij
12
k
2
D
i
c
j
n
ij
;
c
j
nij
jk
jk
Significan
ce of 2
nd
Derivative:
SD
i
p
2
j
2
2
1
2
c<
br>j
100
expjp
2
p
2
where p is the assumed peak width.
k: Go upto c
j
10
-6
Peak
‘found’
when S > Threshold
试验:系
列1为处理后的原始能谱,系列2为5点一阶导数,系列3为5点二阶导数,系列4为对称零面积法寻峰
只要选择好合适的寻峰阈值,足以满足准确寻找到全能峰,并剔除假峰(如康普顿边沿,反散射峰)
5点光滑三阶导数公式判定各感兴区是单峰还是重峰
1
y
i
'''
(y
i2
2y
i1
2y
i1
y<
br>i2
)
2
峰位确定:三阶导数由负变正=0处;峰边界确定:三阶导数由正变负=0处
判定峰是否有意义0.8FWHM ≤ N ≤ 3FWHM
峰高判定条件
|
max
TRHym
p
e0.5
|y
m
这个公式就是在一阶导数法寻峰程序中实际应用的峰高判定条件。
CalculateDifferential(Beginch, Endch, m,
differ);
int
CMmcaView::SearPeakDifferential(int Beginch, int
Endch, int fwhm, int differ[], int m)
{
int n1=0, differ[Endch-Beginch+1], nmax=0, nmin=0,
maxtemp, mintemp,temp;
maxtemp=differ[0]; mintemp=differ[0];
for (int j = 1; j <= Endch-Beginch; j++)
{
temp=differ[j-1];
if(_copysign(temp,differ[j])!=differ[j-1] &&
differ[j]<0) n1=j+Beginch;
if(differ[j]
}
if ((nmin-nmax)>0.8*fwhm && (nmin-
nmax)<3*fwhm) return n1;
else return
(0);
}
对称零面积法(推荐自动寻峰中采用,可探测弱峰和重峰)
面积为零
的“窗”函数与实验谱数据进行褶积变换,且要求“窗”函数为对称函数。对
线性基底的褶积变换将为零
,只有存在峰的地方不为零。
y
'
i
jm
C
m
j
y
ij
j
2
jm
C
m
j
0C
j
C
j
m
1k
2
exp[
2
]
匹配滤波器法(类峰形
函数)
C
j
exp[
2
]
2
<
br>2m1
km
2
~
y
i
峰
判定准则
R
i
~
y
i
jm
Cdata
j
m
ij
1
2
f
m
2
C
j
data
ij
jm
2m+1为变换宽度,
FWHM2.3556
为峰宽参数,若变换后的y'和其均方根误差的比值超
过预先给定
的寻峰阈值(f),则认为找到了一个峰。
峰位的确定:Ri的正极值对应的道址;峰边界的确定:R
i的正峰两边相邻的两个极小值
之间的距离可以作为峰的宽度信息;半宽度:两过零截距。
CalculateArea(0, size, m, fwhm, area, R);
for (int j = m; j <= size-m; j++)
{if(area[j]>0&&R[j]>fh)
for(int
i=1;i<=m;i++)
{
if(area[j-i])>0&&area[j-i]
if
((nmin2-nmin1)>0.6*fwhm && (nmin2-nmin1)<=2*fwhm)
peakposition[p++]=j+0.5;保持峰位对应的道址
}
协方
差法(曲线拟合寻峰,计算机寻峰中采用,可分辨重峰,比较好的寻峰方法,但计
算较为复杂,运算速度
较慢)
1975年等提出了一种新的寻峰方法,称为协方差法。用一个峰形函数与实验
谱数据
逐段拟合(一个高斯形函数与实验谱yi的协方差)
y
ij
y
i
'
C
j
b
i
,Cj为峰形高斯函数
C
j
EXP[2.773(jH
2
)]
H为峰FWHM,y'i为拟合峰高,bi<
br>为本底常数(在峰区内假定不变)
mjm
y
i
'
jm
g
gC
jj
jm
mm
j
jm
mm
j
y
ij
j
2
j
jm
gC
g<
br>jj
jm
m
jm
mm
j
y
ij<
br>
g
gC
jm
m
(
<
br>g
j
C
j
)
2
mm
g
j
g
j
C
j
y
ij
g
j
C
j
g
j
y
ij
y<
br>i
'
jmjm
用
R
i
'
jmjm
f
(f判峰阈值)判定是否存在峰
12
y
i
mm
m
m
22
gggC(gC)
j
j
jj
jj
jm
jm
jmjm
2
4ln2()
Cj通常为纯峰形函数
高斯函数:C
j
exp<
br>
,H为峰的FWHM
H
m
j
g
j为各道计数的权重因子
g
j
1
y
ij
4exp[2(jH)]
或g
j
y
ij
参数选择:H的取值最好与实验谱峰的半宽度接近,2m+1一般取2H左右最好,f一般取
2~5
峰位确定:当Ri为极大值对应的道址;峰边界确定: Ri为负极大值处对应的道址
为了更
好地分辨出落在一个强峰‘肩部’上的弱峰,可以在一个峰的左半部分和右半部
分别计算R
i<
br>值,寻找相互靠得很近的组分峰。
线性拟合寻峰方法(适合于在峰区内分辨重峰)
吸
取匹配滤波器方法的优点,同时用一阶导数法和线性拟合双重峰的技术来提高分辨重
峰的能力,形成了一
种新的寻峰方法,称为线性拟合寻峰方法。
Deconvolution method
First the background is removed (if desired),
then Markov spectrum is calculated (if desired),
then the response function is generated
according to given sigma and deconvolution is
carried
out.
可以提供多种算法,方便自行选择
总结
1.对于弱峰,数据光滑前,高斯乘积函数法和协方差法不能
使用,若先光滑再找峰,又
容易影响重叠峰的分辨;而导数法和对称零面积变换法,无论峰的统计质量如
何,均可
使用。
2.从统计假峰及高基底的抑制能力及重峰的分辨能力来看,一、三阶导数法
和对称零面
积变换法是较好的。对于一、三阶导数法,可先用适当多数据点的一阶导数法找峰,选
取适当的灵敏度常数,以抑制假峰;然后用少点的三阶导数法(或用一阶导数法重复三
次)检查是否有
漏峰和重峰。对称零面积变换法同理。
3.从高基底的抑制能力和弱峰识别的准确度来看,对称零面积
变换法最好。(在计算机
自动找峰程序中,最好采用对称零面积变换法。)
参考资料 4.对找到的峰进行净面积判定是降低假峰出现几率的有效方法。当峰的净面积比峰的总面
积(峰的
净面积和本底面积之和)的标准偏差大若干倍时,才确认该峰是一个真峰,否
则认为它是假峰,予以剔除
。
峰的判弃主要是利用峰面积来进行判定真假峰。对于给定的灵敏因子S,若峰的净面积为NetAR
EA,
峰的宽度为Width。这些参数满足下式认为峰有意义,应保留,否则将找到的此峰丢弃。此式
为:
NetAREAWidth
S
(AREAWidth)
12
S越大灵敏度越高,一般情况下S=3。
参考文献
SPECTRAN-F Version 2 Listings Volume.4 Common
Subroutines and Data structures, CANBERRA
Industries, Inc. 1980
能谱的数据处理(原文).doc
otti, Nucl. Instr. Meth. 50 (1967) 309
Mariscotti Algorithm modified by Routti and
Prussin:
and n, Nucl. Instr. Meth. 72 (1969)
125