多道γ能谱分析软件中寻峰算法比较总结.docx

1、多道能谱分析软件中寻峰算法比较总结自动寻峰由于谱结构的复杂和统计涨落的影响，从谱中正确地找到全部存在的峰是比较困难的。尤其是找到位于很高本底上的弱峰，分辨出相互靠得很近的重峰更为困难。谱分析对寻峰方法的基本要求如下：(1)比较高的重峰分辨能力。能确定相互距离很近的峰的峰位。(2)能识别弱峰，特别是位于高本底上的弱峰。(3)假峰出现的几率要小。(4)不仅能计算出峰位的整数道址，还能计算出峰位的精确值，某些情况下要求峰位的误差小于 0.2 道。很多作者对寻峰方法进行了研究，提出了很多有效的寻峰方法。目的 ：判断有没有峰存在确定峰位（高斯分布的数学期望） ，以便把峰位对应的道址，转换成能量确定峰边界

2、为计算峰面积服务（峰边界道的确定，直接影响峰面积的计算）分为两个步骤：谱变换和峰判定要求：支持手动 /自动寻峰，参数输入，同时计算并显示峰半高宽、精确峰位、峰宽等信息，能够区分康普顿边沿和假峰感兴区内寻峰人工设置感兴趣大小，然后在感兴区内采用简单方法寻峰重点研究：对感兴区内的弱峰寻峰、重峰的分解对于一个单峰区，当峰形在峰位两侧比较对称时，可以由峰的 FWHM 计算峰区的左、右边界道址。峰区的宽度取为 3FWHM ，FWHM 的值可以根据峰位 mp 由测量系统的 FWHM刻度公式 计算。由于峰形对称，左、右边界道和峰位的距离都是 1.5FWHNM 。mLINT (mp1.5FWHM0.5)mRI

3、NT (mp1.5FWHM0.5)式中 mp 是峰位， INT 的含义是取整数。对于存在有低能尾部的峰，其峰形函数描述（参见图） 。y mHEXP (mmp ) 2 / 22 ，mmpJy mHEXP J (2m2mpJ ) / 2 2 ，mmpJ式中 H 为峰高， mp 为峰位， 是高斯函数的标准偏差， J 为接点的道址和峰位之间的距离。在峰位的左侧，有一个接点，其道址为 mpJ。在接点的右侧，峰函数是高斯函数。在接点的左侧，峰函数用指数曲线来描述。这时峰区的左、右边界道址为mINT (mp1.12FWHM 2/ J 0.5J 0.5)LmRINT (mp1.5FWHM0.5)带有低能尾部的

4、峰函数的图形全谱自动寻峰基于核素库法：能量刻度完成后，根据核素库中的能量计算对应的道址，在各个道址附近（左右 10 道附近）采用简单的寻峰方法（导数法）方法：根据仪器选择开发IF 函数法 /简单比较法 (适于寻找强单峰，速度快 )满足条件： datai m datai k datai datai m 可认为有峰存在然后在 data i-m 至 data i+m 中找最大值，对应的道值即为峰位k：找峰阈值，根据高斯分布，一般 k 取值 11.5 常用 5 点、 7 点极大值法（ m 取 2，3）判定峰是否有意义一般，用 R=N0 / Nb R0 确定峰是否有意义R 为峰谷比， R0 为设定值 （

5、经验值）N0 为净峰幅度与基底之和Nb 为基底计数int CMmcaView:SearPeakCompare(int Beginch, int Endch, int m, float k)高斯乘积函数找峰法 （可靠性差，不建议采用）描述谱峰形状的函数主要是高斯函数G (i )Aexp(ii 0 ) 2 / 2 2 则由相邻的数据点定义2一个新的函数（第一高斯乘积函数，只与FWHM2.3556有关）：Pm (i )G(i )G (i m 1)exp( 11.092m ) m2G (i 2)G (i m)H 2m 是步长（用道表示），是高斯乘积函数的阶数，则Pm(i)称为第 m 阶高斯乘积函数。找

6、峰的灵敏度与 m 有关，随 m 的增加灵敏度提高。为避免基线参数的影响，最好扣除本底后，再应用高斯乘积函数找峰。考虑统计涨落的影响，把判断无峰存在的1 变为一个 “单位带 ”。即峰的判断为：1k /yi无峰Pm (i )k /yik 3(1) 有峰峰位的确定：由Pm(i)过 1 的两点求平均来确定；峰边界的确定： “单位带 ”下限的两个最端点；半高宽的确定：函数 Pm(i)在“1上”的截距；组合峰的确定：在乘积函数的两个峰之间没有处于 “带内 ”的乘积函数值导数法（一阶、二阶、三阶）yi 1mC j yi jN m jmNm 为规范化常数， Cj 平滑的变换系数。3 次多项式 5 点光滑一阶导

7、数公式：（可以采用）yi 1 ( yi 2 8 yi 1 8yi 1 yi 2 ) 峰位确定：一阶导数值由正变负 =0 处；峰边界确定：一阶导数12由负变正 =0 处CalculateDifferential(0, size, m, differ);for (int j = m; j = size-m; j+)for(int i=1;i0&differj-imaxtemp) maxtemp=differj-i; nmax=j-i;if(differj+i)0&differj+i0.8* fwhm & (nmin-nmax)3* fwhm )/FWHM 参数根据仪器能量分辨率可人工确定， fwh

8、m20peakpositionp+=j+0.5;/ 保持峰位对应的道址5点光滑二阶导数公式（软件中推荐采用） 1yi (2 yi 2 yi 1 2 yi yi 1 2 yi 2 )/7 点二阶导数(5*(countsdataj-3+countsdataj+3)-3*(countsdataj-1+countsdataj+1)-4*countsdataj)/42;yi1(22.0 yi 3 67.0 yi 2 58.0 yi 1 58.0 yi 1 67.0 yi 2 22.0 yi 3 )252.0软件中推荐采用 11 点以上的公式峰位确定： 二阶导数最小值对应的道址；峰边界确定：二阶导数正极大

9、值点for (int j = m; j = size-m; j+)/m30int maxtemp=-0.5,mintemp=-0.5;If(differj -0.05)for(int i=1;imaxtemp) maxtemp=differj-i; nmax=j-i;if(differj+imintemp) mintemp=differj+i; nmin=j+i;if (nmin-nmax)0.8* fwhm & (nmin-nmax) Threshold试验：系列1 为处理后的原始能谱，系列2 为 5 点一阶导数，系列3 为 5 点二阶导数，系列4 为对称零面积法寻峰只要选择好合适的寻峰阈值

10、，足以满足准确寻找到全能峰，并剔除假峰（如康普顿边沿，反散射峰）5 点光滑三阶导数公式判定各感兴区是单峰还是重峰yi 1 ( yi 2 2 yi 1 2 yi 1 yi 2 )2峰位确定：三阶导数由负变正 =0 处；峰边界确定：三阶导数由正变负=0 处判定峰是否有意义 0.8FWHM N 3FWHM峰高判定条件| ym |max TRH ymp e 0.5 /这个公式就是在一阶导数法寻峰程序中实际应用的峰高判定条件。CalculateDifferential(Beginch, Endch, m, differ);int CMmcaView:SearPeakDifferential(int Be

11、ginch, int Endch, int fwhm, int differ, int m)int n1=0, differEndch-Beginch+1, nmax=0, nmin=0, maxtemp, mintemp,temp; maxtemp=differ0; mintemp=differ0;for (int j = 1; j = Endch-Beginch; j+)temp=differj-1;if(_copysign(temp,differj)!=differj-1 & differj0) n1=j+Beginch ; if(differjmaxtemp) maxtemp=diff

12、erj; nmax=j+Beginch;if (nmin-nmax)0.8*fwhm & (nmin-nmax)3*fwhm)else return (0);return n1;对称零面积法（推荐自动寻峰中采用，可探测弱峰和重峰）面积为零的“窗”函数与实验谱数据进行褶积变换，且要求“窗”函数为对称函数。对线性基底的褶积变换将为零，只有存在峰的地方不为零。yimmC j yi jC j0 C jC jjmjm匹配滤波器法 （类峰形函数） C jj21mk22 exp2 exp2m221 k mmC j datai jyij m1f峰判定准则 Riyim2C 2jdata i jjm2m+1 为变

13、换宽度，FWHM 2.355 6 为峰宽参数， 若变换后的 y和其均方根误差的比值超过预先给定的寻峰阈值（ f ），则认为找到了一个峰。峰位的确定： Ri 的正极值对应的道址；峰边界的确定： Ri 的正峰两边相邻的两个极小值之间的距离可以作为峰的宽度信息；半宽度：两过零截距。CalculateArea(0, size, m, fwhm, area, R);for (int j = m; j 0&Rjfh)for(int i=1;i0&areaj-i0&areaj+i0.6*fwhm & (nmin2-nmin1)=2*fwhm)peakpositionp+=j+0.5;/ 保持峰位对应的道址协

14、方差法（曲线拟合寻峰，计算机寻峰中采用，可分辨重峰，比较好的寻峰方法，但计算较为复杂，运算速度较慢）1975 年 H.P.BLOK 等提出了一种新的寻峰方法，称为协方差法。用一个峰形函数与实验谱数据逐段拟合（一个高斯形函数与实验谱 yi 的协方差）yi j yiC j bi ，Cj 为峰形 /高斯函数 Cj EXP 2.773( j / H 2 ) H 为峰 FWHM ，yi 为拟合峰高， bi为本底常数（在峰区内假定不变） m j mmmmmg jg j C j yijg j C jg j yi jyij mj mj mj mmmmg jg j C j2(g jC j )2jmjmjmmmm

15、myig jg j C j yi jg j C jg j yi jf （f 判峰阈值）判定是否存在峰用 Rijmj mjmj myimmmm1 2g jg jg j C j2(g j C j )2jmjmj mjmCj 通常为纯峰形函数 高斯函数： C jexp4 ln 2 (j )2，H 为峰的 FWHMHgj 为各道计数的权重因子 g j1或g jexp2( jH )4yi jyi j参数选择： H 的取值最好与实验谱峰的半宽度接近，2m+1 一般取 2H 左右最好 ，f 一般取25峰位确定：当 Ri 为极大值对应的道址；峰边界确定： Ri 为负极大值处对应的道址为了更好地分辨出落在一个强

16、峰肩部上的弱峰，可以在一个峰的左半部分和右半部分别计算 Ri 值，寻找相互靠得很近的组分峰。线性拟合寻峰方法（适合于在峰区内分辨重峰）吸取匹配滤波器方法的优点，同时用一阶导数法和线性拟合双重峰的技术来提高分辨重峰的能力，形成了一种新的寻峰方法，称为线性拟合寻峰方法。Deconvolution methodFirst the background is removed (if desired), then Markov spectrum is calculated (if desired), then the response function is generated according to

17、 given sigma and deconvolution is carried out.可以提供多种算法，方便自行选择总结1对于弱峰，数据光滑前，高斯乘积函数法和协方差法不能使用，若先光滑再找峰，又容易影响重叠峰的分辨；而导数法和对称零面积变换法，无论峰的统计质量如何，均可使用。2从统计假峰及高基底的抑制能力及重峰的分辨能力来看，一、三阶导数法和对称零面积变换法是较好的。对于一、三阶导数法，可先用适当多数据点的一阶导数法找峰，选取适当的灵敏度常数，以抑制假峰；然后用少点的三阶导数法（或用一阶导数法重复三次）检查是否有漏峰和重峰。对称零面积变换法同理。3从高基底的抑制能力和弱峰识别的准确度来

18、看，对称零面积变换法最好。 （在计算机自动找峰程序中，最好采用对称零面积变换法。 ）参考资料4.对找到的峰进行净面积判定是降低假峰出现几率的有效方法。 当峰的净面积比峰的总面积（峰的净面积和本底面积之和）的标准偏差大若干倍时，才确认该峰是一个真峰，否则认为它是假峰，予以剔除。峰的判弃 主要是利用峰面积来进行判定真假峰。 对于给定的灵敏因子 S，若峰的净面积为 NetAREA，峰的宽度为 Width。这些参数满足下式认为峰有意义，应保留，否则将找到的此峰丢弃。此式为：NetAREA WidthS ( AREA / Width )1 2S 越大灵敏度越高，一般情况下 S=3。参考文献SPECTRAN-F Version 2 Listings V olume.4 Common Subroutines and Data structures, CANBERRA Industries, Inc. 1980 能谱的数据处理（原文） .docM.A.Mariscotti, Nucl. Instr. Meth. 50 (1967) 309 Mariscotti Algorithm modified by Routti and Prussin: J.T.Routti and S.G.Prussin, Nucl. Instr. Meth. 72 (1969) 125