版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
优质文本优质文本#/84间还要多。下面给出用峰区内协方差法编制的寻峰子程序pkfd7。SUBROUTINEPKFD7(TRH”TRH2,TRH3,MB,MF,NP3)LOGICALIWTESTDIMENSIONSPA(256),DIF1(256)COMMON/CURRP/FW,SG,CDL,CDR,STEP,NOSP,NPK,NPS,NPF,NN,NMCOMMON/FPDATA/IP(15),RP(10)COMMON/PSCB/NCSH,CHPK(15),CSIGMA(15),CCDL(15),CCDR(15)1CSTEP(15)COMMON/ARRAY/SP(21OO),PKL(160),T(128)FSIG(M1)=.1*SQRT(4.*(SPA(M1+2)+SPA(M1_2))+SPA(M1+1)+SPA(M1_1))IG=MBNSP=25610IF(IG.GE.MF)GOTO700CALLGETCRP(SPA,IGNSP,MF)CALLDIFER1(SPA,DIF1,NSP)NSP]=NSP_3iw=trh2I=3IWTEST=IW.LT.1050DO300L]=I,NPS]IF(DI1(L1).LE.0)GOTO300IF(SPA(L1).LE.0.001)GOTO300IF(SPA(L1+]).LT.SPA(L1))GOTO300IF(DIF1(L11).LE.DIF1(L1)_FSIC(L1))GOTO300DT=TRH1*SQRT(SPA(I1))IF(DIF1(L1).LE.DT)GOTO300DT^TRHJSRTCSPAg丿)IF(DIF1(L1]).LE.DT1)GOTO300GOTO310300CONTINUEGOTO500310IPKSTA=L]ISTOP=NSPIF(.NOT.IWTEST)ISTOP=NSP]320IF(L1+2.GT.ISTOP)GOTO600
330IF(SPA(叫).LE.0.001)GOTO200IF(DIF])(M1).GE.0.0)GOTO390SIGF=0.0IF(M].LT.NSP_2)SIGF=FSIG(M1)IF(DIF1(M^1).GT.DIF1(M1)+SIGF)GOTO390IF(IWTEST)GOTO340IF(DIF1(M1).GT.DIF1(MU1)+SIGF)GOTO390340IF(—DIF1(M1)—SIGF.LE.SQRT(SPA(M1))GOTO200390IF(M1.EQ.ISTOP)GOTO600M=M11+1GOTO330200IPKEND=M]W=IPKEND—IPKSTAIF(W.LE.TRH2)GOTO230IPKST]=IPKSTA+IGIPKEN1=IPKEND+IGGOTO250230I=IPKENDGOTO50250CALLPKFD2(TRH3,TRH2,TRH3,IPKST1,IPKEN1,NP3)GOTO230500IG=500IG=NSP+IG—10GOTO10600IF(IPKSTA.LE.100)GOTO650IG=IPKSTA+IG—10GOTO10650IG=IG+100GOTO10700RETURNEND(5)线性拟合寻峰方法[11]前面讨论的匹配滤波器寻峰方法由于有很强的抑制统计涨落的能力,能够在很高的本底上探测弱峰而假峰出现的几率小。但是匹配滤波器寻峰方法分辨重峰的能力较差。吸取匹配滤波器方法的优点,同时用一阶导数法和线性拟合双重峰的技术来提高分辨重峰的能力,形成了一种新的寻峰方法,称为线性拟合寻峰方法。首先分析匹配滤波器法分辨重峰能力较差的原因。图5-1-12画出了谱曲线和由式(5.1.17)计算得出的Rm与m的关系曲线。谱曲线中有五个峰,峰1是个单峰,峰2、峰3和峰4、峰5是两个重峰。由图5-1-12看出,Rm的局部极值只有三个超过了阈值TRH,因而用匹配滤波器法
只能找到三个峰,即mp1、mp2、mp4。错误地把第四和第五两个峰当成了一个峰,漏掉了第三个峰。由Rm的图形可以看到,在mp3的位置Rm仍然有一个局部极值。不过此局部极值淹没在第二个峰产生的负尾巴中,mp3处的Rm值没能超过寻峰阈值。但是mp3处的5Rm超过3阈值TRH。对于由第四、第五两个峰叠加而成的重峰,其特点是峰的宽度较大。为了能够找到全部的峰位,可以米用下述方法。第一步,对谱数据用匹配滤波器进行变换,并求出每点的Rm值。第二步,用一阶导数法在Rm中找出所有的局部极大点。若极大值超过阈值,则该点处至少存在一个峰。第三步,对Rm的局部极大值没有超过阈值,而其净高度5Rm超过阈值的点进行峰净面积判定。如果在峰宽度之内,扣除本底之后的净面积大于总面积误差的若干倍,则在此Rm局部极大值所对应的道址存在一个真峰。第四步,对每个已经找到的峰进行宽度检验。在Rm曲线中,每个峰位两侧负的极小点间的距离W反映出峰的宽度的大小。当W大于某一数值时,认为该峰是一个重峰,否则该峰是一个单峰,这个临界宽度选为2.3FWHM比较合适。判定重峰区之后,下一步的问题是如何确定重峰区中各个组分峰的峰位。文献[12]中介绍了一种线性拟合重叠高斯峰的方法。这种拟合方法不需要迭代初值,即不需要预先确定峰位,而拟合的结果中包含有各个组分峰的准确峰位。用这种方法可以确定重峰区中相互靠得很近的组分峰的峰位。在线性拟合寻峰方法中,使用匹配滤波器进行谱变换,可以提高探测弱峰的能力,降低假峰出现的几率;引入线性拟合重叠高斯峰技术,可以提高分辨重峰的能力,克服匹配滤波器寻峰方法的弱点。所以线性拟合寻峰方法是一种比较好的寻峰方法。图5-1-13是线性拟合寻峰方法的程序框图。图5-1-12采用匹配滤波器进行谱变换之后Rm值与道址关系曲线
图5-1-13线性拟合寻峰程序框图2.几种常用寻峰方法的比较前面叙述了几种常用的寻峰方法的基本原理,下面对这几种寻峰方法的性能和优缺点进行简单的讨论。(1)寻找弱峰的能力由于谱数据中存在很大的统计涨落,在寻峰过程中必然导致两类误差。第一类误差是由于统计涨落而造成的假峰。设寻峰过程中出现假峰的概率为P0。第二类误差是由于很弱的真峰被统计涨落淹没,使寻峰程序漏失了真峰。设P1为寻峰程序能探测出弱峰的概率。P0、P1与寻峰方
法、寻峰阈值、峰高与本底的比值等因素有关。对于好的寻峰方法,应用较小的P0和较大的P1,即能探测出较多的弱峰而假峰出现的几率较小。同时也要求在改变寻峰阈值时P0、P1的变化比较缓慢,以使这种寻峰方法有更大的适用范围和稳定性。文献[13]中从理论上计算出假峰出现的概率P0。假定谱数据是只包括常数本底的白噪声时Po=Po=Jfds卜drEXP[—2(S2—2tsr+r2)]—gf2(5.1.45)其中,f=trh®[-912何1)2]-2,tViqq=B艺C2,Q=q=B艺CC,C.为滤波器的冲击响应函数,TRH为寻峰11j1221jj—1Jj=—Kj=—K阈值,B为本底谱的高度。在匹配滤波器、广义二阶差分滤波器等三种不同寻峰方法的情况下,由式(5.1.45)计算出的P0与TRH的关系如图5-1-14所示。由图可以看出,提高寻峰阈值可以有效地降低假峰出现的概率。在相同的寻峰阈值的情况下,使用匹配滤波器寻峰方法假峰出现的概率最小。虽然提高寻峰阈值能够降低假峰出现的概率,但是必然会使真峰漏失的概率增大,即探测出弱峰的概率P1降低。为了评价不同寻峰方法探测弱峰的本领,以实验的方法做出P1与P0的关系曲线。实验是用计算机生成的模拟实验谱进行的。模拟谱中各道计数=B+XHEXP[-j==B+XHEXP[-j=1(m—m)2pi—2q2]+R(m)<B+XPHEXP[1-2q2(m-m)2pj——]>2q2(5.1.46)其中ym是m道的谱数据,B=2500为各道的本底计数,是一个常数,q=4为高斯峰函数的标准偏差,mpj为峰位,Np为峰的个数,H为峰高,R(m)是均值为零、方差为1的符合正态分布的伪随机数系列。图5-1-14假峰出现的几率与寻峰阈值TRH的关系曲线匹配滤波器寻峰方法;广义二阶差分寻峰方法;矩形滤波器寻峰方法。
以不同的寻峰方法编制成FORTRAN程序,在不同的相对峰高a=H、帀的情况下对模拟谱寻峰,计算出P0和P],其曲线如图5-1-15所示。wnIL—亠_wnIL—亠_]fLu-11图5-1-15图5-1-15不同寻峰方法P1与P0关系曲线图5-1-16计算机生成的二重峰模拟谱曲线由图可见,匹配滤波器法具有最好的探测弱峰的能力,曲线比较平缓。在假峰出现的概率较小的P1较大,仍然能够找到大部分弱的真峰。一阶导数寻峰方法寻找弱峰的能力最差。(2)分辨重峰的能力能否分辨出相距很近的峰是衡量寻峰方法好坏的另一个重要指标。以计算机生成的模拟谱对各种寻峰方法检验表明,各种寻峰方法分辨重峰的能力有很大的差别。仍然用公式(5.1.46)生成模拟谱。设置不同的峰位mp和峰高H可以得到几组模拟谱数据。对于二重峰(图5-1-16)峰高比分成1:1、1:0.5、1:0.25、1:0.1(各组分峰的峰高为20000、10000、5000、2000)四组,每组谱数据中峰位距离又分成0.5q、Q、1.5Q、2q、2.5Q、3Q、3.5Q、4q、4.5Q、5q(Q=4)共十种情况。三重峰(图5-1-17)中,组分峰的峰高比分成1:0.5:1、1:0.25:1(组分峰的峰高分别为20000、10000、50000)两组。每组谱数据中峰距又分成q、1.5q、2q、2.5q、3q、3.5q、4q、4.5q、5q(q=4)共九种情况。用模拟谱对几种不同的寻峰方法进行检验,各种不同寻峰方法能够分辨出的最小峰位距离列于表5.6中。由表中的数据可以看出,线性拟合寻峰方法有图5-1-17计算机生成的三重峰模拟谱曲线由上面的讨论可知,不同的寻峰方法各有各自的特点。匹配滤波器法寻找弱峰的能力最强,
假峰出现的概率小,适于弱放射性核素的分析。而线性拟合法具有最好的的分辨重峰的能力,适用于通用的Y谱分析系统中。矩形滤波器法和一阶导数法虽然寻峰能力较差,但计算方法简单、运算速度快,适用于运算速度比较慢的微型计算机系统。表5.6几种常用的寻峰方法分辨重峰能力的比较能分辨的最、小峰距方法、二重峰一重峰峰高比峰高比1:11:0.51:0.251:0.11:0.5:11:0.25:1线性拟合法O1.5Q1.5Q3.5Q3Q3.5Q协方差法1.5Q2.5Q2.5Q3.5Q3.5Q4Q匹配滤波器法2.5Q2.5Q3.5Q4Q4Q矩形滤波器法2.5Q3Q3.5Q3.5Q4.5Q一阶导数法4Q4.5Q4.5Q45Q三、本底扣除在能谱曲线中,我们感兴趣的峰常常是叠加在很高的本底谱上。在谱定量分析的过程中,为了计算某种能量射线的强度需要求出峰的净面积为了计算峰的净面积,必需扣除峰区内的本底。这里我们讨论的本底不是指被测样品之外的其它放射源的干扰,而是指由被测射线与探测器(或其周围介质)通过不同的物理过程产生的或被测样品在射线的作用下通过不同的激发过程而造成的干扰。例如在Y能谱分析中,一般利用光电峰的面积来计算核素的活度。但是光电峰常常叠加在更高能量的Y射线的康普顿电子谱上。各种更高能量的Y射线在探测器中产生的康普顿电子谱叠加在一起构成了本底谱。又如在带电粒子激发的X射线测量中,我们测量的是入射粒子在样品中激发出的特征X射线的能量和强度,但是入射粒子产生的轫致辐射、样品中二次电子的轫致辐射、被测样品中由核激发态产生的Y射线引起的康普顿散射等等都能产生很强的本底。扣除本底的方法可分为二种:一种是设法求出本底谱在整个测量能区中的分布模式,从整个谱数据中逐道减去本底在该道的计数。得到不包含本底的谱数据,这种方法称为全谱本底扣除法。另一种方法是在包含有我们感兴趣的峰的一个很窄的谱段中,认为本底谱是按直线或多项式规律分布,在这个谱段中从谱数据中逐道减去由直线或多项式分布模拟的本底数据。这种方法称为峰区本底扣除法。1.全谱本底扣除法使用全谱本底扣除法扣除本底必须首先求出在整个谱的范围内各道的本底值,然后逐道由谱数据中减去各道的本底计数,才能得到扣除本底之后各道谱数据。在某些物理测量中,可以写出本底谱分布的解析表达式。例如在电子激发的X射线分析中,入射电子在靶中与原子核的库伦场作用时产生的轫致辐射是造成本底的主要因素。轫致辐射谱是随能量变化的连续谱,可以写出轫致辐射谱的解析表达式。但是在很多物理测量中要找到一个能够在全谱范围内描述本底分布的
解析表达式是很困难的。在这种情况下,可以设法从所测得的能谱数据中找出本底分布的规律。其基本思想是在谱数据中找出位于各个峰区之外并且能够反映出本底大小的足够多的数据点。以一定的方法把这些数据点连接起来形成本底谱。然后再由谱数据中减去本底谱得到扣除本底之后的谱数据。下面以用Ge(Li)探测器测得的Y谱为例来讨论这种本底扣除方法[14]用Ge(Li)探测器测量Y能谱时,形成本底谱的因素有两个:第一是高能Y光子在探测器中的康普顿散射产生的康普顿谱边缘和平缓的康普顿坪。另一个因素是Y光子的能量在探测器灵敏感区中全部被吸收之前,由于小角散射产生能量损失,致使全能峰在低能部分产生一个拖长的尾部。我们把这部分能量吸收不完全的光子造成的低能尾部也看成是峰下本底的一部分,在本底扣除的过程中对它进行校正。在一个峰区内低能尾部事件在每道计数中所占的比例是该峰内比这一道高的全部峰净计数的函数。因此,低能尾部事件造成的本底的分布函数应是峰净计数分布(峰区内谱曲线减去康普顿电子谱)的积分函数。峰区内的本底分布可以看成是连续的康普顿电子谱(可以近似地用直线来模拟)和峰净计数的积分函数两部分叠加而成的连续曲线。图5-1-19用插入新的极小值的方法修正本底谱图5-1-19用插入新的极小值的方法修正本底谱在整个谱区内本底分布的具体计算步骤如下。①对谱数据进行平滑处理。②把整个谱分割成若干个相邻的谱段,每个谱段的宽度取为整个谱中峰的平均FWHM的三倍或五倍。在每个谱段中找出谱数据的最小值,并把最小值所对应的道址记录下来。为了不丢失能代表本底谱的数据点,这个步骤可以用不同的谱段宽度重复地进行二次。第一次谱段宽度为3FWHM,第二次为5FWHM。③剔除重峰区中两个峰之间的谷点。当很多个峰重叠在一起时,尤其是重峰区较宽时,可能会把重峰区中峰谷误认为是本底点(图5-1-18)。因此,在第二步中找出的各最小点的横坐标间的距离小于2FWHM时,把相邻的两个点中计数高的那一点作为重峰区中的峰谷予以以剔除。④把找到的最小点用直线联成本底曲线。如果在任何一个0.5FWHM的区间内,直线段上的各个点都高于谱数据,则认为这个线段不符合真正的本底谱,必需加以修正。如图5-1-19所示,A、B、C三个点是第二步中找到的最小点,虚线是由直线连接而成的本底谱曲线。在BC区间内(BC>0.5FWHM),由直线连接而成的本底谱中的各点都比谱数据大,则在本底谱与谱数据之间相差最大的点D插入一个新的最小值。用直线连成修正之后的本底谱如实线所示。进行上述修正可以使由直线连成的本底谱更接近于真实的本底谱。⑤在每个区间内由谱数据和直线本底之差计算积分函数,把积分函数规一化,使它和直线本底叠加之后,在区间的两个边界上(即上两步骤中找出的最小点)本底数值维持不变。然后再一次采用第四步骤中的方法对本底谱进行修正。⑥最后对由上述步骤计算出的本底谱数据进
行平滑处理,就得到了比较符合实际情况的本底数据。上面讨论了在整个谱区的范围内计算本底分布的方法。求得每道的本底数据之后,逐道地由谱数据减去本底数据,就得到了扣除本底之后的谱数据。文献[14]中指出,采用上述方法扣除本底之后再进行谱的定量分析,可以大大地提高分析精度,还可以更多地分析出谱中存在的弱成分。2.在峰区内扣除本底的方法当探测器的能量分辨率比较高,能谱曲线中峰的宽度比较窄时,最常用的方法是峰区内扣除本底的方法。首先在峰位的两侧找出峰区的左、右边界道址mL和mRO由于峰区的宽度很窄,可以用通过mL道和mR道两点谱数据的一条直线来模拟峰区内本底的分布(图5-1-20中的虚线)描述本底分布的另一种方法是在峰区外每侧各取若干点,用多项式对这些点进行函数拟合,用拟合所得的多项式来模拟峰区内本底的分布由谱数据逐道地减去用直线或多项式模拟的本底数据就得到扣除本底之后的谱数据。易而易见,在探测器能量分辨率很差、峰很宽时(例如用NaI探测器测得的丫谱),用这种方法将会产生很大的误差。使用峰区内扣除本底的方法扣除本底时,若要获得较好的精度,关键在于峰区左、右边界道址的选择。若峰区选得过窄,则用来模拟本底的直线偏高,扣除本底之后的谱数据偏小。当峰区选得过宽时,峰区的左右边界道可能落在邻近的另一个峰的尾部上,也会造成本底扣除的误差。图5-1-20用直线来模拟峰区的本底分布下面我们讨论确定峰区左、右边界道的方法。对于一个单峰区,当峰形在峰位两侧比较对称时,可以由峰的FWHM计算峰区的左、右边界道址。峰区的宽度取为3FWHNM,FWHNM的值可以根据峰位mp由测量系统的FWHNM刻度公式计算。由于峰形对称,左、右边界道和峰位的距离都是1.5FWHNM。m=INT(m—1.5FWHM+0.5)(5.1.47)LPmr=INT(m+1.5FWHM+0.5)(5.1.48)式中mp是峰位,符号INT的含义是取整数。当峰形函数在峰位两侧不对称时,峰位不在峰区的中点,需要根据峰形函数分别计算出左、右边界道址。例如,对于存在有低能尾部的峰,其峰形函数可用式(5.1.490)来描述(参见图5-1-21)。(5.1.49)y=HEXP[—(m—m)2/2q2],m三m—Jm(5.1.49)y=HEXP[J(2m—2m+J)/2q2],m
式中H为峰高,mp为峰位,b是高斯函数的标准偏差,J为接点的道址和峰位之间的距离。在峰位的左侧,有一个接点,其道址为mp—J。在接点的右侧,峰函数是高斯函数。在接点的左侧,峰函数用指数曲线来描述。这时峰区的左、右边界道址为m=INT(m—1.12FWHM2/J—0.5J+0.5)(5.1.50)LPm=INT(m+1.5FWHM+0.5)(5.1.51)Rp图5-1-21带有低能尾部的峰函数的图形另外一种确定峰区边界的方法是用峰区外两侧谱数据中的平稳点。这是一种根据谱的实际形状来确定峰区边界的方法。其基本思想是由峰位向两侧外推,当谱形变得趋于平缓时,即认为找到了相应的峰区边界道址。具体的计算公式见式(5.1.27)—(5.1.33)。计算峰区边界点的第三种方法是用数字滤波器处理谱数据,由滤波之后的谱数据中的特征点来决定峰区的边界道址。如图5-1-6所示。经过匹配滤波器滤波之后的谱数据在峰位两侧有两个负的极小值。可以把这两个极小值所对应的道址作为峰区的边界道址。必需指出,上面讨论的确定峰区边界道址的方法都只适用于单峰区,确定重峰区边界道址的方法将在以后进行讨论。在峰区内扣除本底方法中,产生误差的另一个原因是由于在边界道中计数较低、统计涨落较大,从而造成了通过峰区边界点的本底直线的不确定性。为了减小这种误差,左右边界道中的谱数据常常采用相邻几道谱数据的平均值,以减小统计涨落的影响。上面讨论了扣除本底的二种方法:全谱本底扣除法和峰区内本底扣除法。一般说来,这两种方法都只适合于由Ge(Li)、Si(Li)等能量分辨率较高的探测得的能谱。在谱的定量分析中,首先用全谱本底扣除法扣除本底,然后再进行寻峰、重峰分析等操作。这样可以提高谱定量分析的精度,提高测量系统对弱成分的灵敏感度。但是这种方法的计算工作量比较大。在探测器的能量分辨率比较高,峰区的宽度不大时,使用峰区内本底扣除法也能得到比较好的精度。因而峰区内本底扣除法是目前Ge(Li)、Si(Li)探测器测谱定量分析中最广泛使用的一种方法。四、直接法计算峰的净面积峰净面积的计算是能谱定量分析的主要内容之一。在峰区内谱曲线和横轴之间所包围的面积称为峰的总面积(At);本底曲线之下所包围的面积称为本底面积(AB);峰的总面积减去本底面积就是峰的净面积(An)o峰的净面积和探测器所接收的某种能量的射线的强度成正比。
计算峰净面积的方法有二种:第一种方法称为直接法。在直接法中,首先确定峰区的左、右边界道址mL、mR,在mL、mR范围内求出谱数据的各道计数和,得到峰的总面积。根据本底分布曲线求出峰区内各道本底的计数和,得到本底面积。由峰总面积减去本底的计数和,得到本底面积。由峰总面积减去本底面积就得到峰的净面积。计算峰净面积的第二种方法是函数拟合法。在拟合法中,通常用一个角解析表达式即谱函数来描述峰区内谱的分布。谱函数可以看成是若干个峰函数和本底函数的叠加。在峰区范围内,用最小二乘法以谱函数拟合谱数据,可以求出峰函数中的各个参数。在峰区范围内计算峰函数的定积分可以求出峰的净面积。图5-1-22峰的净面积和本底面积直接法比较简单,运算速度快。在峰的高度比较大,能量分辨率比较高,峰区比较窄的情况下能够得到比较满意的计算精度。但是一般说来,直接法只适于单峰的净面积计算。在重峰区内不能正确地计算出各个组分峰的净面积。函数拟合法的计算过程比较复杂,其优点是计算精度高。在峰高比较小,本底比较大的情况下也能得到比较满意的计算结果。函数拟合法不仅适合于单峰区,也适合于双峰,特别是在重峰区内能够计算出几个相互重叠的组分峰各自的净面积。文献[15]中介绍了很多种计算峰净面积的方法。几种方法的主要区别在于如何选择边界道址、如何扣除本底、各道计数相加时赋予各道计数相等的权重还是按计数大小而赋以不同的权重。这几种方法各有特点,下面进行简单的介绍。1.全峰面积法(TPA法)在全峰面积法中,峰区边界道址选在峰本身的贡献很小而边界道中的计数基本上反映本底计数的地方。当谱的能量分辨率较高、峰区较窄时,一般以通过谱数据中峰区边界点的直线来描述本底的分布。由峰区内谱数据各道的计数和减去梯形本底面积得到峰的净面积片一一厶图5-1-23峰外本底扣除法计算净面积
AN=ATAN=AT-AB-(ymL+y)(m-m+1)/2mRRL(5.1.52)式中,mL、mR是峰区的左、右边界道址,y、y是mL、mR道的谱数据。由式(5.1.52)可LRmmLR以看出,对于一个叠加在高本底上的弱峰,净面积AN的相对统计误差很大。为了减少峰区边界道计数y、y的统计涨落对峰净面积计算的影响,可以把峰边界道外侧若干道的计数平均起LmR来作为边界道计数来计算本底。如图5-1-23所示,设在峰区左侧取nL点,峰区右侧取nR点。峰左侧本底各道的计数和为ABL(5.1.53)ABL(5.1.53)m=+任峰右侧各道的计数和为ABRm=mABRm=m+1Ry(5.1.54)则峰的净面积为A=(2则峰的净面积为A=(2Ry)—(-bl+-BR)(m—m+1)/2NmnnRLm=LR=(送y)—(红+仏)ynn2m=LR(5.1.55)若取n=n=n,则有LRA=(送y)—N(A+A)Nm2nBLBRm=根据误差传递公式可以计算出峰净面积的统计方差(5.1.56)Q2AN—m、‘AL)2(Q2AN—m、‘AL)2(%+2n2LnR(5.1.57)当nL=nR=n'N=m一m+1时,
RLQ2AN=(送y)—(弓)2(A2BLm=mL+ABR)(5.1.58)由于全峰面积法计算区内整个峰的面积,因此当因测量系统的计数率改变而引起峰的形状的改变时或者在整个谱区内各个峰的宽度不同时,都能够得到可靠的计算结果。但是当峰区很宽时,峰区内的本底分布不能用直线模拟,否则将产生比较大的误差。对于Ge(Li)、Si(Li)等探测器来说,由于能量分辨率较高,一般峰区很窄。因此,全峰面积法是Ge(Li)、Si(Li)谱仪中使用最广泛的计算峰净面积的方法。2.Covel1法计算峰的净面积在Covell方法中,峰区的边界道不是选在峰位两侧谱曲线变化的平稳点上,而是选在峰位两侧谱曲线斜率较大的点上(图5-1-24)。以连接两个边界点的直线描述本底的分布,由峰区的总面积减去本底面积得出峰净面积。设在峰位mp两侧各取N点计算出峰区边界道址m=m—NLPm=m+NRp(5.1.59)(5.1.60)则净面积为AL=(mRym)-(N+新mT+ymR)m=mL式中,y、y分别为峰边界道址mT、LmRL净面积的统计方差为q2=A+(N+1)(N-i)(yANN'22mR处的谱数据mL+ymR)(5.1.61)(5.1.62)图5-1-24Covell方法计算峰的净面积用Covell方法求得的峰净面积比用全峰面积法计算的峰净面积小。显然,用Covell法计算出的Y谱光电峰的净面积不等于在探测器内产生光电作用的光子数,但是这个净面积与产生光电作用的光子数成正比。因此,如果谱仪系统在进行效率刻度的时和谱定量分析都同样地采用Covell法计算峰的净面积,仍然可以正确地计算出Y射线的强度。Covell方法的优点有二个:第一,由于峰区比较窄,本底分布可以更好地用直线模拟。特别是当峰座落在康谱顿边缘等非直线本底上的情况下,在较窄的峰区宽度内用直线方程模拟本底的分布也能得到较好的计算精度。因此,Covell法不仅适用于由Ge(Li)、Si(Li)等探测器获得的谱中峰净面积的计算,而且也适用于由NaI探测器获得的Y谱中峰净面积的计算。第二,由于边界道址中的计数y、y较大,其相对统计涨落较小,因而提高了本底扣除的精度,减小了LmR峰净面积计算的误差。在Covell方法中,峰边界道的选择对峰净面积计算精度有很大的影响。从理论上说,存在着一个最佳的N值使得峰净面积的相对统计误差Qan/an最小。当N取得比较大时,Covell法接近全峰面积法。当N取得过小时,计算出的净面积太小,其误差也大。另一方面当多道分析器系统的计数率变化时,常常引起峰的形状发生变化。一般当计数率增高时,能量分辨率将显著变坏,峰展宽,峰高降低。在Covell法中,由于边界道选在峰两侧斜率较大的点上,峰形状的改变将引起净面积数值的变化,使得峰净面积的计算结果与系统的计数率有关。在计数率变化时,峰净面积的计算产生附加的误差。在短半衰期核素的能谱测量中,注意到这一点是很重要的。总之,Covell法与全峰面积法相比各有优缺点。前者对本底的形状不敏感,适用于不同探测器测得的谱,但易受峰形变化的影响。全峰面积法正好相反,对本底的形状比较敏感,但受峰形变化的影响较小。全峰面积法多用的能量分辨率较高的谱的定量分析中。3.Wason方法计算峰的净面积Wason方法把全峰面积法和Covell方法结合起来。和全峰面积法一样,以峰两侧的平稳点之间联成一条直线来模拟本底的分布,而峰区的左、右边界道址的选择方法和Covell法相同(图5-1-25)。设决定本底分布的峰两侧平稳点为mB、mB,峰区边界点为mL、mR。图5-1-25中阴blbrLR影部分所表示的峰净面积为An=(送y)-(N+2)(y+y)(5.1.63)Nm2mBLmBRm=式中,y、y分别为本底边界点mB、mB道中的谱计数。净面积的统计方差为BLmBRBLBRQ2=(送y)-(N+£)2(A+A)(5.1.64)ANm^^BLBR''m=在Wason方法中,峰区内道数较少,因而本底分布的形状对净面积计算的影响较小。又由于本底取得较低,提高了峰净面积和本底面积的比值。和Covell方法一样,在Wason方法中由于峰边界道选在峰两侧斜率较大的点上,峰净面积对峰形的变化很敏感。当由于测量系统的计数率的变化或其它原因引起峰形状有较大的变化时,峰净面积的计算将会产生较大的附加误差。图5-1-25Wason法计算净面积4.Quittner方法计算峰的净面积在前面讨论的三种计算峰净面积的方法中都是用直线来模拟峰区内本底的分布。在计算叠加在高本底上的弱峰时,这种模拟将产生很大的误差。在Quittner方法中用三次多项式曲线来描述峰区内本底的分布。首先在峰外两侧选取两个平稳的谱段进行二次多项式拟合(图5-1-26中峰两侧的实线)。在这二个似合谱段中选取两个参考点mB、mB,其高度分别为y、y,LRmBLmBR拟合的二次曲线在参考点处的斜率分别是qL和qR。假定峰区内的本底分布符合三次多项式,经过简单的推导可以得出第m道的本底TOC\o"1-5"\h\zd(、「―q,—2q3(y—y)*、B=y+q(m—m)+[RL+IbrIbl](m—m)2(5.1.65)mmBLLBLl+l(l+l)2BL(5.1.65)LRvLR7q+q3(y—y)+[-^L乩+IbrmB^](m—m)3(l+l)2(l+l)3BL'lr,'lr7式中B是第m道的本底值,lL、lR分别是峰两侧参考点和峰位间的距离。mlR峰的净面积为AN=(5r(y+B)(5.1.66)Nm=在峰的净面积和总面积相比很小的情况下,用Quittner方法计算峰净面积要比Covell法精确。此外,由于用曲线来模拟峰区内本底的分布,可以允许峰区宽度比较大。因而这种方法也可以用
于NaI谱中峰净面积的计算。图5-1-26Quitter法计算峰的净面积Quittner方法的主要限制是在峰的两侧必需有较宽的平缓谱段,也就是该峰与邻近峰的距离要比较远。在几个峰之间的距离比较近时,用这种方法计算峰净面积将产生较大的误差。五、复杂谱的解析前面介绍了计算单峰净面积的几种方法。当分析由多种放射性核素组成的混合样品谱时,能谱曲线的形状比较复杂。特别是相距很近的峰叠加在一起,形成了重峰。用计算单峰净面积的方法不能计算出重峰中各个组分峰的净面积。这时必须对复杂谱进行解析才能求出各个组分峰的净面积,从而计算出各种射线的能量和强度,确定混合样品中各种核素的含量。复杂谱中的各道计数可以看成是各种能量的射线在该道产生的计数的线性叠加。由这个原理出发,形成了几种不同的解谱方法。解谱的方法可以分为两类:一类使用标准谱,包括剥谱法、逆矩阵法和最小二乘法解谱等;另一类不使用标准谱,常用的是函数拟合法。用标准谱的方法解谱时,首先必须知道被测样品中存在着哪几种核素,分别单独测出每种核素的标准谱。可以认为被测样品谱是各标准谱的线性叠加。计算出被测混合样品中各种核素强度的比例关系。为了保证线性叠加的假定成立,标准谱和被测混合样品谱必须在同样的条件下测量。谱仪系统的能量分辨率、探测效率、能量刻度等必须保持不变,尤其是不应随计数率的不同产生显著的变化。由于这些因素的限制,标准谱方法解谱大多用于不太复杂的谱的解析。使用标准谱解谱的方法读者可参阅有关文献[15]由Ge(Li)、Si(Li)等能量分辨率比较高的探测器测得的、X射线能谱的解析通常使用函数拟合法。在这种方法中,首先把谱划分成若干个谱区间,每个谱区间包含若干个有意义的峰。在每个谱区间中写出表征谱形的谱函数解板表达式。由测得的谱数据计算出谱函数中的各个参数。由这些参数可以计算出这个谱区间内每个组分峰的净面积,从而求出样品中各核素的强度。函数拟合法解谱的最大优点是不必预先知道样品中含有的核素的种类。此外,在测量系统的计数率等测量条件变化使谱形有较大的改变的情况下,这种方法仍然能够达到比较满意的精度。既不需要标准谱,也提高了谱数据处理的自动化程度。这些优点使函数拟合法解谱得到了广泛的应用。
最小二乘法函数拟合法解谱的基本原理在谱曲线的每个谱区间中包含有若干个叠加在本底上的峰。在不考虑谱数据的统计涨落时,在某一个谱区间中谱曲线可以用函数yi=F(i,P),(i=0,1,……,n-1)(5.1.67)来描述。n为该谱区间内谱数据点数,yi为不考虑统计涨落时谱区间内第i点的谱数据,P=(P1,P2,…,PK)为待定参数向量,K个参数决定了谱曲线的形状。这些参数包括本底函数中的本底参数和峰函数中的峰形参数。例如,我们假定本底函数用二次多项式描述,峰函数用高斯函数描述,则谱函数可以写为(5.1.68)=F(i,P)=P1+P2i+P3i2(5.1.68)叽EXP—(i•-P3j+2)2/2+3]j=1(i=0,1,,n—1)式中j为该谱区间内峰的序号:Np为谱区间内峰的个数;P1,P2,P3是描述本底曲线形状的本底参数;P3+1为第j个峰的峰高;P3j+2为第j个峰在该区间内的位置;P3j+3为第j个峰高斯函数的标准差,代表了第j个峰的宽度。写出谱函数表达式之后,把谱函数和测得的谱数据进行拟合,以求出K个待定参数。在最小二乘法意义上的函数拟合就是求解待定参数P,使满足各点谱数据yi和谱函数的残差平方和为最小,即:Q=y.-F(i,P)]2(5.1.69)iii=0为最小。Q称为目标函数:Wi为第i点的权重因子,在初次拟合时可取Wi=1/yi;y.为测得的谱在谱区间中的第i个数据。为了使Q达到极小,各参数应满足如下的方程组:电丽1(5.1.70)<(5.1.70)2即残差平方和对各参数的偏导数都是零,这个方程组称为正规方程。求解正规方程,可以求出K个待定参数,这样就求出了谱区间内各个峰的峰位、峰高的峰宽度。对每个峰的峰函数积分,可以计算出谱区间内各个峰的净面积。求解正规方程的数值计算方法依正规方程的性质不同而异。当正规方程是一个线性方程组时(例如峰函数是高斯函数,而峰位和峰宽为已知量不参加拟合时),可以采用消去法或分解法等求解线性方程组的计算方法。当正规方程是非线性方程组时,通常采用迭代法求解。应用最多的是高斯-牛顿法的改进算法:麦夸脱(Marquardt)法和变尺度法。迭代法是一种试探性的数值求解方法。我们以高斯-牛顿法来说明求解正规方程的步骤:先给P设定一组初始近似值P(o),(j=1,2,…,K),称为迭代初值,并记真值与初值之差为A.,则jP=P(0)+Aj,jJ(j=l,2,于是确定件的问题化为确定修正值Aj。对谱函数F(',P)在p(0)附近作台劳展开,并略去高次项可得dFF(i,P)-F+-°0dF1其中F=F(i,P(0),P(0),…P(0))0\12K7人dFA+oidF2(5.1.71)dFdF(i,P,P,…P)0=^^2K—dF.dP.JJP=P(0)11P=P(0)11(5.1.72)P=P(0)KK将式(5.1.71)代入式(5.1.69)可得dF人dF+亠A+.…0-dP1dF1K)]2(5.1.73)dQdPi將=2W事寻+…ji=01JSwdFKi=0idFKj如果设-Ewi討i-F0)]i=0j(5.1.74)ajiajyn-1dFdFi=0idFdFiiW些(y-F)i=0idFi0j(l,j=1,2,…,K)(5.1.75)可以得到下述方程组aA+aA++aA=a1111221kk1yaA+aA++aA=a211222kk2y(5.1.76)Wk1A1+ak2A2+…+akkAk=aky求解这个K阶线性方程组可以求出A,进而计算出P,于是P,=P(o)+A.(j=1,2,-K)(5.1.77)由于在台劳展开的过程中略去了高次项,因而由式(5.1.77)算出的q并不是真值。为了更进一步逼近真值,可以用求出的P.代替原来的初值Pj®,重复上述运算而得到新的Aj,这个过程反复迭代直到IAJ的值小于预先给定的误差为止。'J综上所述,用高斯-牛顿法求解非线性最小二乘问题,可归纳为如下的步骤:人为地指定迭代初值P(0),(j=1,2,-K),并记P=P(o)+A.,把求解pj化求解A.的jjjjJj问题。用台劳展开式把求解Aj的问题线性化,变成求解一个线性方程组。以Aj修正P.作为新的初值,重复①、②的步骤直至满足预先给定的误差条件为止。在迭代过程中,如果AJ逐次减小,逐渐逼近P.的真解,则迭代过程是收敛的。但如果初值选得与真值相差较大,将使台劳展开式严重失真,使迭代所得的P.更远离真解,这时迭代过程不收敛,出现发散现象。在能谱数据似合过程中,作为初值的峰位、峰高、峰宽和本底参数很难预先选得很精确。在迭代过程中可能会出现收敛速度很慢,甚至出现发散现象。为了放宽对迭代初值的限制,麦夸脱提出了改进方案。麦夸脱法与高斯-牛顿法基本相似,其差别仅在于把迭代时求解A.的线性方程组改写jS1+d)A1+ai2A2+…+aiKAK=气(5.1.78)ai2A1+(a22+d)A2+…+a2KAK="2y(5.1.78)SK1A1+aK2A2+…+(aKK+d)AK=°K为与方程式(5.1.76)不同,这个方程在系数矩阵A的对角元中加入了一个常数d,称为阻尼因子,加入阻尼因子之后,求解式(5.1.78)实际上是求解A.使jRF©FQ鼻£W」yi-(F+詁A1+…+詁Ak)]2+d送A2「0暫1叭K冃j=Q+d送A2jj=1达到极小值。在迭代过程中不断改变阻尼因子。开始d选取比较大的数值,使得麦夸脱法接近梯度法。理论上可以证明,只要d值选得足够大,总会使残差平方和Q值随着迭代过程逐渐减少。但阻尼因子越大,步长A越小,收敛速度越慢。为了加快迭代的速度,在迭代过程中按一定的比例逐渐地减小阻尼因子d。在d接近为零时,麦夸脱法退化为高斯-牛顿法。阻尼因子较大时,A很小,在这种情况下,以误差条件来判定迭代结果的精度可能出现假收敛的现象。因此,只有在d逐渐减小到一定的数值以下时,才能进行是否满足收敛条件的判定。实际计算的经验表明,麦夸脱法可以允许初值有更大的选取范围。在初值偏离真值20—30%的情况下,仍然能够得到比较好的计算结果。这在进行函数拟合时是很有利的。有关函数拟合计算方法更详细的讨论参见有关文献[16]。在谱曲线的函数拟合中,谱的拟合区间中峰的个数决定了待定的拟合参数的个数。在谱区间内最多允许存在几个峰取决于下面两个因素:首先,当拟合参数过多时,系数矩阵A可能呈现严重“病态”造成求解的困难。在实际应用中,当数据点300个、待定参数达到45个时,麦夸脱法仍然有效。另一个需要考虑的因素是待定参数多时将占用更大的计算机内存空间,运算时间也要增多,这就需要使用更大、更快的计算机。在用最小二乘法函数拟合法解谱时,为了能得到较好的精度,必须处理好以下几个问题:首先是如何合理地划分谱拟合区间,使得在谱区间中包含有应该分解的重峰,但谱区间又不能太大,以节省计算机的内存和运算时间。第二是正确地选择峰函数和本底函数的数学表达式,使得在不同的物理测量条件下,测量系统的计数率不同时都能够正确地反映出谱曲线的真实形状,以提高函数拟合的精度。第三是如何正确地选择迭代初值,使初值尽可能地接近解的真值,以加快迭代
收敛的速度,避免出现发散现象。第四是合理地选用迭代收敛判定条件,当拟合结果不理想时,如何对拟合进行改善。最后必须指出,在拟合区间内找出全部有意义的组分峰的峰位是保证函数拟合结果的可靠性的先决条件。由式(5.1.68)可以看出,谱函数的解析表达式是建立在各个组分峰的基础之上的。如果寻峰程序不能分辨出重峰中组分峰的个数,漏失了有意义的峰或出现了假峰,则谱函数表达式就不正确。在此基础上进行函数拟合就不会导致正确的计算结果。下面我们对上述几个实际问题进行讨论,最后再简单地介绍适用于微型和小型计算机的非迭代函数拟合方法。拟合区间的划分通常,谱数据的函数拟合不是在整个能谱范围内进行,而是把整个谱划分成若干个小的谱区间(在每个谱区间内包含有诺干个有意义的峰),这些谱区间称为拟合区间。划分谱区间的目的是为了提高拟合精度和节省计算机内存。由式(5.1.68)可以看出,谱函数包括两个部分:本底函数和峰函数。在某些情况下(例如在电子激发的X射线能谱分析中),可以近似地写出本底函数的理论模型。但是,在很多情况下,在整个谱范围内写出本底函数的解析表达式是比较困难的。而在比较小的谱区间内,本底函数可以近似地用多项式描述。显然,在谱区间很大时,用多项式描述本底将会造成较大的误差。此外,谱拟合程序在计算机内运行时,程序占用的内存容量与参加拟合的数据点数和拟合区间内峰的个数有关。在较小的拟合区间内进行拟合时,占用的内存容量较小,这在使用小型或微型计算机时尤为重要。划分拟合区间在寻峰之后进行,根据谱中各个峰位之间距离划分拟合区间。拟合区间的划分原则如下:相互重叠的峰应处于同一个拟合区间之内。拟合区间的边界应位于本底点上。孤立的单峰自己占有一个拟合区间。拟合区间内包含的数据点数和峰的个数不应超出拟合程序所允许的最大数值,否则必须进行特殊的处理。根据上述原则,在不同的函数拟合程序中,可以采用不同的划分拟合区间的方法。例如在SAMPO程序[17]中,如果两个峰位之间的距离小于6FWHM,则认为这两个峰是相互重叠的峰,应放在同一拟合区间内拟合。如果第三个峰的峰位与第二个峰的峰位之间的距离仍然小于6FWHM,则这三个峰的同一拟合区间内进行拟合,依此类推。以重峰中最左边的一个峰的峰位向左再推移3FWHM作为拟合区间的左边界道址;从最右边的峰位向右推移3FWHM作为拟合区间的右边界道址。在谱数据涨落比较大或探测器的峰函数拖有较长的尾巴时,为了更精确地反映本底的形状,可以把拟合区间再向外扩展若干道,但是不能超过谱曲线的局部最小点,进入另一个峰的下降段。又如在小型计算机上运行的SPECTRAN-F程序⑶采用如下的方法划分拟合区间:对于任何两个相邻的峰,如果峰位距离小于3FWHM,则认为是重峰,应放在同一个拟合敬意中。在一个拟合区间内,最多可以包含七个相互重叠的峰。如果拟合区间内只有一个孤立的单峰,由峰位向两侧各扩散1.5FWHM得到拟合区间的左右边界。和SAMPO程序相比,Spectran-F程序中的拟合区间要小得多,适合在内存容量不大的小型计算机中运行。3•谱函数的选择函数拟合法解谱的主要思想是通过调整谱函数中的参数逼近原始谱数据,从而求解出一组最
优化的参数。因此,预先假定的谱函数F(i,P)是不能够反映谱曲线的真实形状是非常重要的。如果谱函数曲线的形状与实测的谱曲线之间相差较大,将会产生较大的似合误差,特别是在计算重峰区中的弱组分峰时误码差更大。由射线和探测器相互作用的基本理论和电子学系统的响应函数出发,写出一个简单而又实用的谱函数的解析表达式是比较困难的。目前大都是使用经验公式描述谱函数。在通用的Y射线或X射线谱仪中使用的谱函数应满足下述要求。在拟合区间内,谱函数应能正确地描述不同种类[如平面型Ge(Li)、同轴型Ge(Li)、高纯Get和Si(Li)探测器]和不同性能探测器测出的能谱曲线的形状。图5-1-27中画出了用同轴型Ge(Li)探测器测得的能量为1332keV的丫射线能谱曲线的一部分。由此图可以看出,在峰顶附近,谱曲线近似地符合高斯曲线,由峰两侧的形状可以看出在峰下有一个台阶状的本底;在峰的低能一侧有一个由于射线能量在探测器的灵敏区内不完全吸收造成的拖长的尾部。不同质量的探测器测出的谱的形状有明显的差别。在由高质量探测器测出的谱曲线中,峰的宽度比较窄,低能尾部比较小。当探测器的质量比较差时,峰的宽度比较宽、低能尾部也较高。对于不同类型和不同质量的探测器,峰低能侧尾部的开始点出现在0.01—0.5倍峰高的范围内。值得指出,既使对于同一个探测器,在使用了较长的时间之后,其能量分辨率也会变坏,低能尾部将会抬高。谱函数能否正确地描述低能尾部的形状对解谱的精度会产生很大的影响。例如,在一个很强的峰的低能侧存在一个弱峰时,这个弱峰很可能叠加在强峰的低能尾部上。如果谱函数不能正确地描述探测器响应的低能尾部的形状,就不可能正确地计算出弱峰的大小。谱函数应能适应测量系统计数率的变化。当测量系统的平均计数率增高到每秒1万次以上时,由于电子学系统中的脉冲堆积效应,在峰的高能一侧将产生明显的拖长的尾部。采用脉冲堆积排除电路、极-零补偿、基线恢复器等电子学技术可以减少脉冲堆积效应的影响,但不能完全消除。此外,电子学系统的漂移和能谱稳定器的使用也会使峰高能侧的尾部加大。因此,为了避免计数率变化带来的误差,谱函数应能适应由于计数率不同和电子学系统的漂移造成的峰形状的改变。这在短寿命核素的测量中尤为重要。上面所述的两个要求使得在解谱过程中对弱峰能够达到比较高的灵敏度,对重峰有比较强的分辨能力。在两个峰相距很近(例如峰位距离为FWHM)而峰高比又很大时(如两峰的高度比为100)能够正确地计算出弱峰的强度。谱函数表达式应尽可能地简单,待定参数抢占数不能太多,以减小计算工作量,提高程序运行的稳定性和可靠性。这在使用小型或微型计算机解谱时更为重要。谱函数表达式中的某些参数在整个测量能区内随能量的变化应是平缓的或者这些参数的变化对拟合结果的影响不大。这样,在函数拟合之前我们就能够比较准确地估计出这些参数的值,或者将这些参数取为固定的常数,以减少计算工作量。在拟合区间之内,谱函数可以看成为峰函数和本底函数的叠加。把谱曲线看成是连续曲线时F(x)=B(x)+豹S(x)(5.1.79)=1其中,在谱函数的连续函数表达式中x是自变量(道址),S(x)代表峰函数,B(x)是本底函数,Np为拟合区间内峰的个数。下面首先讨论本底函数。在Y谱中形成本底的来源有三个。首先是其它辐射源产生的辐射本底,其次是更高能量的Y射线的康普顿散射。第三种来自被测的这种
能量Y射线本身,由于丫射线能量在探测器的灵敏区中没有完全被吸收,造成某些计数分布在比全能峰能量更低的各道内,形成了台阶状的本底分布。虽然很多文献都把台阶状的分布看成是本底的一部分,但是这部分计数实际上是该种能量Y射线强度贡献的一部分,因此台阶的高度与峰的净高度成正比。在Y谱中,本底函数B(x)是由上述三种原因产生的本底的线性叠加。一般说来,辐射本底是一个常量,与道址无关;高能T射线的康普顿电子谱可以近似地用多项式分布描述;关于台阶状的本底,文献中提出了很多种经验公式。文献[18]中综合了近年来发表的描述台阶状本底的多种经验公式,如表5.7所示。表中列出的函数是规一化之后的台阶状函数。其中x是道址;X是峰位的道址;yi是第i道中谱的净计数;W是峰的FWHM;b是高斯分布的标准偏差;erf是误差函数\-2nJx\-2nJxEXP(-t2)dt02erf(x)=2-erf(x)。表5.7台阶状本底函数的经验表达式*序号发表日期台阶状本底函数表达式BS1969阶跃函数fc;:BS21969V1—EXP[P(x一X)/2b2],x<X01,x>XBS319701厶X—Pb—x;+arctg-2兀'Pb1BS41970V1—1/[1+(x—X-)2],x<Xb0,x>XV,BS51977V1—1/[1+(—1)2],x<Xb0,x>XBo19731/[1+EXP[(x—X)/0.75W]S6BS71974erfc[—(x—X)/P]]BS819762erfc[(x—X)/j2b]BS91975V2[2—EXP2(x—X)/b]x<X2EXP2[—2(x—X)/b]x>X、2BS1019791—为y./为y.iii=xi=0*表中的P],P2是待定的拟合参数。
这十种台阶函数按其形状可以分为三类,图5-1-28中画出了以B、B和B为代表的三S1S4S8类台阶状本底函数的图形。函数B的形状与B相似,而B、B、X、BB与B相S2S4S3S5S6S7S9S8似。图5-1-28三类台阶状本底函数曲线在早期的一些丫谱定量分析程序(如GAMMA-S、SPECTRAN-F、SAMPO等)中,没有考虑台阶状本底的存在,直接用多项式描述整个本底函数。文献[18]指出,引入台阶状的本底函数可以显著地提高拟合质量。采用同样的峰函数,以同样的拟合方法对由不同的探测器测量的光电峰进行拟合表明,在其它条件相同的情况下,引入台阶状本底函数可以使单位自由度的残差平方和XR由2.97下降到1.92。但是采用几种不同的台阶状本底函数相比较,拟合质量相差不大。值R得指出的是,当一个弱峰位于一个强峰的低能侧并和强峰重叠时,引入台阶状本底函数可以大大地提高对弱峰的探测灵敏度和分析精度。综上所述,本底函数B(x)可以用于个多项式(通常用一次或二次多项式)和台阶状本底函数之和描述B(x)=P+Px+Px2+B(x)(5.1.80)1x3s引入台阶状本底函数可以提高拟合质量,但不同种类的台阶函数对拟合质量的影响不大。因此,从简化拟合程序来说,应尽量采用简单的台阶状函数表达式。峰函数是谱函数表达式中的重要组成部分,所选用的峰函数是否能够反映峰的真实形状对拟合质量有很大影响。前面已经指出,一个单峰的形状不仅和探测器的类型和质量有关,测量系统的平均计数率、电子仪器的噪声、成形滤波网络的选择、极-零补偿、基线恢复和稳谱器的使用等因素对峰的形状也都有很大的影响。一般情况下,单峰的峰顶附近可以用高斯函数描述。在峰的两侧,则偏离高斯函数,出现拖长的尾部。在峰位两侧高斯函数和尾部函数的连接点分别称为左、右接点。左接点的位置和低能尾部的高度和探测器的类型和质量有关。分辨率高、质量好的探测器接点距离峰位比较远,接点道址内的计数与峰高相比较小,低能尾部所占的面积对峰的净面积贡献不大。能量分辨低的探测器或探测器的质量变坏时,接点道址距峰中心较近,接点道内的计数较高,低能尾部的面积在峰的净面积中占有较大的比例。右接点的位置和高能尾部的大小与谱仪系统的计数率和脉冲堆积的情况有关,当测量系统内脉冲堆积严重、主放大器中的成形滤波网络选用不当、基线恢复器性能不良时,峰高能侧的尾部将变大。在系统中使用稳谱器时,也可能造成峰高能侧的尾部变大。
近十年来,对峰函数的近似解析表达式进行了大量的研究工作。一类研究工作是寻找能精确地表达峰形状的经验表达式,尽可能地提高拟合质量。这些央函数一般比较复杂,每个峰函数中包含有八个甚至十几个待定的拟合参数。另一类研究工作是寻找一些比较简单的峰函数(有两个到四个拟合参数),在一定条件下能够达到比较满意的拟合质量。这些比较简单的峰函数适于在小型或微型计算机上运行的拟合程序中采用。文献[18]中列举了近年来发表的和一些著名的丫谱分析程序中使用的15种峰形函数,并对这些峰形函数的性能进行了实验比较。这些峰函数大致可分为几类:①标准高斯函数②在峰区内各段使用不同的函数,这几个函数在接点处平滑地连接起来。例如,用两个b不同的高斯函数分别描述峰两侧的形状,这两个高斯函数在峰顶处平滑地联接起来。也可以在峰顶附近采用高斯函数,左、右接点之外用指数函数描述尾部的形状,这三个函数在左、右接点处平滑地联接起来。③高斯函数与不同数量、不同形式的尾部函数叠加。④采用变形的高斯函数和其它函数。表5.8列出了几种典型的规一化之后的峰函数表达式。表5.8几种典型的峰函数表达式序号PpS2PS3发表日期1969SPECTRANF(1978)SAMPO(1969序号PpS2PS3发表日期1969SPECTRANF(1978)SAMPO(1969)GAMANL台阶状本底函数表达式EXP[—(x—X)2/2b-2]EXP[—(x—X)2/2b2,EXP[J(2x—2X+J)/2b2,EXP[—(x—X)2/2b2,EXP[J(2x—2X+J)/2b2,EXP[J'(2X—2x+J')/2b2,EXP[—(x—X)2/2b2],C24WMC(1972)C24WMC(1972)EXP[—(x—X)2/2b2]+P^EXP[-P2(x—X)]{1—EXP£(x—X)2/2b2},x>XPS5(1973)PSPS5(1973)PS6HYPERMET(1976)EXP[—(x—X)2/2bf,EXP[-(x-X)2/2bg,e—(x—X)2/2b2+丄p2e(x-x)/P2x—X2bP51erfc(4—+5)+P32bP52P223e(x—X)/PerfC(+'畫卩。)2bp^2p4EXP[—(x—X)2/2b2]+片EXP[-(x-P2)2/2P32]PS7(1979PS7(1979)EXP[—(x—X)2/2b2]+片EXP[-(x-P2)/什]x>P2-警2P4为了比较不同峰函数的拟合质量,文献[18]中对由不同质量的探测器测得的60Co的能量为1332keV的单峰,以相同的拟合方法、不同的峰函数对数据进行拟合,计算出的单位自由度残差平方和xR列于表5.9中。D]、D2、D3是三个质量不同的探测器,D1是灵每区体积为65cm3的同轴Ge(Li)探测器;已经使用了九年,峰的FWHM为2keV,在峰的低能侧有比较大的尾部;
D2是新的同轴Ge(Li)探测器,体积为120cm3,峰的FWHM为2keV,在峰低能侧的尾部也比较大;D3是13cm3的平面型Ge(Li)探测器,峰的FWHM为1.7keV,用它测得的谱中峰的低能尾部最小。由表5.9中可以看出:①采用高斯函数作为峰函数时拟合质量最差。当使用较好的探测器时,拟合质量得到改善。由此说明在峰函数中引入一些参数来描述峰形状中的拖长的尾部对改善拟合质量是非常重要的。在探测器的质量比较差时,这一点尤为重要。②在表5.8中列举的峰函数中,当采用拟合参数的个数多的峰函数时,能得到更好的拟合质量。此外,需要指出,在计算表5.9中的结果时,谱数据是在低计数率(约2000计数/秒)的情况下测得的,测量系统中使用了极-零补偿电路和基线恢复器。在获取的谱中,峰的高能侧基本上没有拖长的尾部。因此,表5.9的数据中没有包括由于计数率的变化所引起的高能尾部产生的影响。表5.9不同峰函数拟合结果比较峰函数序号待定参数个数x2RD1D2D3PS1360.049.913.2PS2440.126.513.5PS468.153.153.04PS776.572.522.76PS81.541.772.36上面讨论了用不同的峰函数拟合一个单峰时拟合方差的大小。在实际应用中需要考虑的另个重要的问题是在拟合一个重峰区时叠加在强峰两侧的弱峰的分析灵敏度和精度。设想当一个很弱的峰叠加在一个强峰的尾部上时,如果选用的峰函数不能准确地描述尾部的变化形状,就不能正确地计算出弱峰的净面积。在这方面,很多文献的作者作了不少的工作。文献[19]中,采用误差函数描述台阶状本底的分布,以高斯函数加上指数尾部作为峰形函数,对65Zn和46SC的能量为111.56keV和1120.3keV的两个相互重叠的峰(峰距为1.8FWHM)进行分析,在两个峰的强度比为0.01时,计算每个峰的强度,系统误差小于2%。文献[20]中,采用不对称高斯函数作为峰形函数,同时使用阶跃函数或在低能部分按指数变化的台阶状本底函数。当两个峰的距离为3FWHM、强度比为0.001时,射线强度计算结果的系统误差小于2%。甚至在分析由116mIn和59Fe产生的确良1093keV和1095keV的两个峰时(峰距为FWHM),在两个峰的强度相差两个数量级时,仍能得到满意的计算结果。在选用峰函数时,不仅要考虑到拟合质量,还需要考虑一些其它的实际因素,其中最主要的是在拟合过程中需要求解的参数的个数。在参数较少时,可以使用内存容量较小、运算速度较慢的小型或微型计算机,而且函数拟合的过程也比较稳定。因此,对于通用的谱分析程序,一般选用含有3—5个拟合参数的峰函数。在计数率不高,测量系统中包含有脉冲堆积排除电路、基线恢复器和准确地调整极-零补偿电路的情况下,可以选用高斯函数附加低能尾部函数构成峰函数。在计数率较高或系统中包含有能谱稳定器的情况下,需要在峰函数中引入描述高能尾部的函数。迭代初值与迭代收敛条件的选取
在函数拟合法解谱过程中,求解非线性最小二乘问题通常使用迭代法。在迭代过程开始之前,首先要选取各个待定拟合参数的初值P=(P(0),P(0),,P(0))。还要规定迭代的收敛条件和迭代12K终止条件。当收敛条件满足时,就求得了各个拟合参数的估值。当迭代收敛的速度很慢或迭代发散时,由迭代终止条件来中止迭代过程。图5-1-29迭代初值选取方法示意图在选取代初值时,应使各参数的初值尽可能地接近这些参数的真值。这样可以加快迭代收敛的速度,如果初值偏离真值太远,可能造成迭代发散,不能求出正确的结果。在函数拟合中也可以使用一些不需初值的计算方法估算初值,例如可以采用随机尝试法删取初值。但是在Z、X射线能谱分析中,最常用的方法最藉助于寻峰和系统刻度的结果,以一些简单的方法近似地求出各个参数初值。下面我们以一个具体的例子来说明初值的选取方法。图5-1-29中画出了一个拟合区间中的谱曲线。假定峰函数选用高斯函数,本底函数选用二次多项式。在峰区内由寻峰程序找到了两个峰,其峰位的道址分别是mp1、mp2,峰区内谱函数的离散量形式的表达式为y=P+Pi+Pi2+PEXP[—(i—P)2/2P2]i1234L\5丿6」(5.1.81)+PEXP[—(i—P)2/2P2]789(i=0,1,2,・・・n—1)其中,yi为拟合峰区内第i点的谱数据。式(5.1.81)中共有9个特定的拟合参数。在选取拟合参数的初值时,首先用前面已讨论过的确定峰区边界的方法在谱数据中找出峰区的两个边界点(y、mL)、(y、mR),以通过这二点的直线近似地代替本底函数。用寻峰程序找到的峰位TOC\o"1-5"\h\zLmR•L<R.mp1、mp2计算出峰位参数P5、P8的初值;用H]、H2作为峰高参数P4、P7的初值;P6和P9代表高斯函数的标准偏差,其初值可以根据峰位的道址由系统的峰形刻度参数计算出来。由此得出P(0)=ymLP(o)=(y—y)/(mR—mL)mmRLRLP(0)=03P(o)=y一[P(o)+P(0)(m一m)]mp112p1LP(o)=m—mp1LP(0)=0根据系统的峰形刻度曲线计算6P(0)=y—[P(0)+P(0)(m一m)]7mp212p2L
P(°)=m-mp2LP(0)根据系统的峰形刻度曲线计算上面以高斯峰函数为例讨论了如何选取拟合参数9的初值。当采用其它形式的谱函数时,也可以使用类似的方法选取拟合参数的初值。下面我们讨论如何确定迭代程序中的收敛条件。在函数拟合的迭代过程中,在两种情况下终止迭代。一种情况是满足了迭代收敛条件,这时对拟合参数的估值达到了预先给定的精度,函数拟合取得了满意的结果。另外一种情况是当迭代达到了预先设定的次数或连续几次迭代都不能使迭代质量得到明显的改善时,也要使迭代过程中止。在后一种情况下,函数拟合过程即使再进行下去也不能取得满意的结果。出现这种情况一般是由于迭代初值选取不当,谱函数的表达式与真实的谱形相差太大,寻峰程序漏失了真峰或混入假峰等因素造成的。在某些拟合程序中,以拟合参数的相对误差小于某一数值作为迭代收敛条件。例如在SAMPO程序中,迭代的收敛条件为:在相邻两次迭代中,求出的各个拟合参数相对误差小于10-8。此外,以最大迭代次数为100或连续四次迭代中单位自由度的残差平方和X2没有显著的改善作为迭代的终止条件。在SAMPO程R序中,收敛条件是比较苛刻的,达到这个收敛条件需要花费较多的运算时间。但是SAMPO程序是在比较大的计算机上运行的,这个问题并不大。对于在运算速度比较慢的小型或微型计算机上运行的拟合程序,为了减少运算时间,收敛条件必需放宽。实际上,我们可以由各个拟合参数的物理意义出发,合理地分别规定每个参数的必需达到的精度。这样,不仅迭代收敛条件容易满足,而且也能使拟合结果达到满意的精度。例如,在SPECTRAN-F程序中,迭代的收敛条件规定如下:由相邻的两次迭代求出的峰位值差小于0.01道。即18XjlVO.Ol(j=1,2…,竹),Np为该拟合区间内峰的个数。相邻的两次迭代求得的峰宽值相差小于千分之一,即|SqI/g<0.005。b为描述峰曲线峰顶部分的高斯函数的标准偏差。相邻的两次迭代求出的各个峰的峰高的变化小于峰高的标准误差的一半。当峰高的标准误差的一半小于2时,取为2。相邻两次迭代中,尾部接点J的变化小于Q/<2,即lSJI<qI•込。当上述四个条件同时满足时,则认为迭代过程收敛,各个参数的估值达到了满意的精度。需要指出,在SPECTRAN-F中使用麦夸脱法求解非线性最小二乘问题时,为了避免出现由于阻尼因子d太大而造成假收敛,在迭代过程中只是在阻尼因子下降到小于一定的数值之后方使用收敛条件判定迭代过程是否应该结束。如果经过多次迭代,收敛条件仍不能满足,在下述二代迭代终止条件之一来结束迭代过程。一个是迭代次数达到了15次。另外一个终止条件是当经过了四次迭代之后,继续进行迭代时目标函数Q下降不显著。即当lQ(h+1)-Q(h)|/Q(h)<Q-3时,迭代过程终止。其中Q(h+1)为在第k+1次迭代之后由式(5.1.69)计算的残差平方和,Q(h)为第k次迭代过程之后求出的残差平方和。在这两种情况下,收敛条件没有满足,对谱数据进行的函数拟合没有取得满意的结果。函数拟合法解谱的误差在函数拟合法解谱过程中,通常以XR检验和衡量拟合质量。RXR=」F[y.-F(i,P)]2/y.(5.1.82)Rn-J11
XR为单位自由度的残差平方和。n为拟合区间内的数据点数,K是特定拟合参数的个数,n-KR称为自由度,yi是由拟合区间左边界算起的第i道的谱数据。根据一般误差理论,XR也是一个随机变数,有其自己的分布规律。如果函数拟合的结果F(i,P)与yi的真实分布吻合,XR值应在一个合理的数值范围内。如果XR值超出了这个合理的数值范围,则说明F(i,P)不能反R映谱的真实形状,拟合过程中计算出的参数P的误差比较大或者是一组完全不可相信的结果(例如有漏失成分等)。实际上,XR的大小与自由度的大小和谱数据的统计涨落的大小有关。一些文献的作者根据对谱数据拟合的实际经验提出了XR的合理的数值范围。R.L.Heath[2i]等认为xR应在2—4之间。B.R.Kowalskis]认为x2可以接受的数值应在1—25之间。P.A.AarniQ23]指出,根RTOC\o"1-5"\h\z据实际经验XR可以接受的数值应小于25。因此,可以认为当XR在1—5之间时,拟合质量比较RR好,而当XR>25时,拟合结果是不可靠的。在SPECTRAN-F程序中,当XR>4时,程序将打RR印出信息,提示用户需要检查在拟合过程中是否存在比较大的系统误差,以便进一步改善拟合质量。由XR的定义中可以看出,XR包括了两类误差,一个是各道谱数据本身的统计涨落造成的偶RR然误差;另一部分是拟合中存在的系统误差。确切地说,只有第二部分误差才反映了拟合质量的优劣。为了确切地反映拟合中产生的系统误差的大小,有些作者提出了一些新的误差检验因子[23-25]。在对谱函数的优劣进行比较的试验性计算中,或在系统的峰形刻度等程序中,使用这些新的检验因子能够更灵敏地测试出由于谱函数选用不当,峰形参数不精确等因素造成的系统误差。由函数拟合求解出的参数P的方差可以由下式估算q2Pj=(A)-1”XR(5.1.83)TOC\o"1-5"\h\z式中Q2P/是由拟合求出的第j个参数的方差,A是当迭代终止时由式(5.1.76)表示的线性议程组的系数矩阵,(A)-1..是系数矩阵的逆矩阵中的第j个对角元素。由式(5.1.83)表示的P.的方程ijJ仍然是一个近似估计,没有考虑到各个参数之间的协方差。另外Q2P.中也没有包括由系统刻度引起的误差。J各拟合参数的方差求出之后,根据所采用的峰净面积计算公式可计算出各组分峰的净面积的方差。例如在SPECTRAN-F程序中,由下式计算各个组分峰的净面积A=HA/豹H(5.1.84)JJtk'丿k=1其中,Aj为重峰区中第j个组分峰的净面积,H.为拟合中求出的第j个组分峰的峰高;Np为重峰区内组分峰的个数;At为重峰区的总净面积。At是把重峰区看成是一个单峰区采用前面讨论过的计算单峰净面积的方法计算出来的。由误差传递公式计算出A.的方差QA.JAj°AjAtHk一H°AjAtHk一Hj1\2宠Hk丿k=1丿2H°2H+j—j宠Hk-k=1-2°At(5.1.85)AHkAk=1式中°H2式中°H2是第j个组分峰峰高的方差;°At是重峰区总净面积At的方差。上面讨论了采用函数拟合法解谱求出的各个拟合参数和净面积的误差。下面我们讨论影响函数拟合法解谱的精度和可靠性的一些实际因素,主要是重峰区内各组分峰的峰位之间的距离和各组分峰的峰高之间的比例对拟合精度和拟合结果的可靠性的影响。在一定程度上,这个因素决定了函数拟合法分解重峰的能力。前已讨论过,正确地求解重峰和弱峰首先要求谱函数能够正确地反映测量系统的响应函数,也就是谱函数能正确地描述谱曲线中峰和本底的大小和形状。但是,即使是采用了比较准确的谱函数,在组分峰相距很近,峰高比很大时,函数拟合法解谱并不能得到满意的结果。由理论上准确地分析这个问题是比较困难的,在一些文献中7.17]采用实验的方法分析函数拟合法解谱的误差与组分峰的距离和峰高比的关系,得出了有意义的结果。文献[17]对SAMPO程序进行了实验研究。SAMPO程序采用表5.8中序号这PS3的峰函数,本底函数采用二次多项式。首先由寻峰结果和峰形刻度参数进行一次线性拟合求出各拟合参数的初值,然后用可变尺度法进行的函数拟合,求出峰高、峰位等参数,用峰函数的数值积分计算组分峰的净面积。用一个实验谱测试SAMPO程序分解重峰的能力。在实验谱的一个拟合区间中包含一个二重峰,这个重峰位于每道计数为1800的康普顿坪上。重峰中能量较高的组分峰的峰高是固定的,峰高与康普顿坪高度的比值是1.6。在组分峰位距离(以FWHM为单位)不同、峰高比不同的情况下,进行函数拟合,计算出两个组分峰的强度,求出峰强度比值。在表5.10中列出了计算结果。表5.10SAMPO程序求解重峰的精度拟合
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 《妈妈的爱》教案设计
- 人教版五年级数学上册【全册教案】
- 河流清淤疏浚服务合同
- 人教版二年级语文下册教学计划教案
- 物流配送中心管理准则
- 墙面施工合同:美术馆内部装修
- 网络综艺策划副导演聘用协议
- 宾馆水暖设施更新工程协议
- 工业用水施工合同
- 大同市工业档案管理准则
- 肺结节介绍课件
- 山西陆合集团恒泰南庄煤业有限公司矿山矿产资源开发、地质环境保护与土地复垦方案
- 酒店账单-水单-住宿
- 2023年山东春季高考数学试题word版(含答案解析)
- 我的连衣裙【经典绘本】
- 农村公路畅通工程质量检测方案第三方检测及交工验收
- 中国石油化工集团公司职工违纪违规行为处分规定
- 深圳市某河道排涝工程监理规划
- 机械加工工时定额标准计算手册
- 课堂教学评价标准
- 2021年中国环卫行业及环卫设备(环卫装备)行业现状及趋势分析
评论
0/150
提交评论