γ能谱的数据处理.doc

精品文档能谱的数据处理由多道脉冲分析器获取的谱数据需要以一定的数学方法进行处理才能得到实验要求的最终结果。能谱的数据处理大致可以分为两个步骤。首先进行峰分析，即由能谱数据中找到全部有意义的峰，并计算出扣除本底之后每个峰的净面积。第二步是放射性核素的活度或样品中元素浓度的计算，即由峰位所对应的能量识别出被测样品中含有哪些放射性核素或被激发的元素，并且由峰的净面积计算出放射性核素的活度或元素在样品中的浓度。采用不同的物理实验方法，使用不同的探测器时，能谱的数据处理方法也有所不同。在本章中首先讨论在各种能谱数据处理中经常用到的峰分析方法，包括谱数据的平滑处理，本底扣除、寻峰、峰净面积计算和函数拟合法解谱。然后，以谱分析为例，讨论基于计算机的多道脉冲分析系统中的谱自动分析软件的工作原理。第一节 常用的峰分析方法一、谱数据的平滑处理由于射线和探测器中固有的统计涨落、电子学系统的噪声的影响，谱数据有很大的统计涨落。在每道计数较少时，相对统计涨落更大。谱数据的涨落将会使谱数据处理产生误差。其主要表现为在寻峰过程中丢失弱峰或出现假峰、峰净面积计算的误差加大等等。谱数据的平滑就是以一定的数学方法对谱数据进行处理，减少谱数据中的统计涨落，但平滑之后的谱曲线应尽可能地保留平滑前谱曲线中有意义的特征，峰的形状和峰的净面积不应产生很大的变化。对谱数据进行平滑处理通常使用数字滤波器。由信号分析理论的观点出发，我们可以把原始谱数据看成是噪声（即谱数据中的统计涨落）和信号（即峰函数和本底函数）的叠加。经过数字滤波器的处理可以提高信号噪声比。如图5-1-1所示，令第x道的原始谱数据为y(x)，经过数字滤波之后的谱数据为(5.1.1)其中，为数字滤波器的单位冲击响应函数，并有(5.1.2)y()Y(x)Y(x)图5-1-1 用数字滤波器进行谱的平滑处理由于谱数据是离散量，公式（5.1.1）、（5.1.2）的离散量表达形式为(5.1.3)(5.1.3)只要选择恰当的数字滤波器响应函数，就能够使平滑后的谱既保留了原始谱中的峰和本底的形状和大小，又得到最佳的信号噪声比。由频域的观点分析，谱中的统计涨落，即噪声的频谱分布在整个频率范围内，而峰函数和本底函数的频谱主要集中在低频范围。因此，使用一个低通滤波器进行滤波，可以使峰和本底信息都通过滤波器到达输出器，而噪声中的高频成分被滤波器抑制，从而提高了平滑后谱中的信号噪声比，减小了谱数据的统计涨落。图5-1-2中画出了二种常用的数字滤波器分别在道域和频域中的响应函数的图形。1. 最小二乘移动平滑方法1964年A. Savitzky和J.Egolay1提出了一个用于谱数据平滑处理的滤波器响应函数。其基本思想是，当求平滑之后谱的第m点数据时，先在原始谱数据第m点的左、右各取K个数据点，形成一个共有2K+1个数据点的窗口。在这个窗口中用多项式拟合原始谱数据，则拟合多项式在m点的值就是平滑后的谱在m点的值。当m值沿谱数据移动时，就可以得到整个平滑后的谱数据。这种方法称为最小二乘移动平滑法，或最小平方曲线拟合平滑法。原始谱数据为，平滑后谱数据为，在平滑窗口内，用q价多项式逼近原始谱数据时，平滑后谱第m点的值为(5.1.5)图5-1-2 SAVITZKY滤波器的响应函数(a) SAVITZKY滤波器在道域中的响应函数；(b) SAVITZKY滤波器在频域中的响应函数；1. SAVITZKY5点平滑滤波器；2. SAVITZKY11点平滑滤波器同时还可以把S（x）在m点的各阶导数值作为平滑后的谱在m点的各阶导数值。平滑后的谱在m点的各阶导数值。平滑后的谱在m点的p阶导数值为(5.1.6)根据上述原理，用最小二乘法函数拟合可以导出计算平滑后的谱数据和其各阶导数值的具体计算公式(5.1.7)规范化常数NK和权因子CK,j的值列在表5.1中。由表中查出NK和CK,j的值就可以写出平滑谱的计算公式。例如当平滑窗口选为5点时（K=2），5点平滑公式为(5.1.8)在文献【2】中，由数字滤波器理论可以推导出最小二乘移动平滑公式（5.1.7）中Ck,j/NK的一般计算公式。当平滑窗口为W=2K+1时，(5.1.9)KjK用这个公式计算得出的值，与表5.1中列出的数值相吻合。表5.1 最小二乘移动平滑计算公式中的NK值和CK,j值2K+1CK,j8-21-7-6-78-67-13-11-518420-36-4278799-21-334122164414-2-2391472169393-3-1421622484546120431672589597171421622484546122391472169393-3334122164414-24278799-21518420-3667-13-117-6-788-21NK32311051434292312135表5.2 最小二乘移动平滑法计算平滑谱的一阶导数公式中的NK值和CK,j值2K+1CK,j8748-7-9812922-6-643-41211133-5-930-14150-660300-4-1002-18334-1578-29486-3-902-17842-1796-532-14222-2-673-13843-1489-503-193-671-1-358-7506-832-296-126-58-80000000013587506832296126588267313843148950319367-13902178421796532142-2241002183341578294-86593014150660-30066434121-1133798-129228748NK-23256334152240245148118825212表5.3 最小二乘移动平滑计算平滑谱的二阶导数公式中的NK值和CK,j值2K+1CK,j840-72591-6125222-51191115-4-8-82628-3-15-29-5-175-2-20-48-10-6-802-1-23-53-13-9-17-3-10-24-56-14-10-20-4-21-23-53-13-9-17-3-12-20-48-10-6-8023-15-29-5-1754-8-8262851191115612522272591840NK397661881001429462427由公式（5.1.7）也可以计算平滑谱的各阶导数值，只不过权因子CK,j和规范化常数NK的值各不相同。表5.2中列出了采用不同的平滑窗口、用公式（5.1.7）计算平滑谱的一阶导数时的CK,j与NK的值。表5.3中列出了采用不同的平滑窗口，计算平滑谱的二阶导数时的CK,j与NK的值。根据表5.2平滑窗口为5点（K=2）时，平滑谱的一阶导数计算公式为(5.1.10)根据表5.3，平滑窗口为5点时，平滑谱在m点的二阶导数值为：(5.1.11)前面已经指出，平滑的本质是对谱曲线进行低通滤波，去掉高频成分，保留有用的低频信息。滤波的效果取决于低通滤波器的频谱特性。当式（5.1.7）中的权因子不同时滤波器的频谱特性不同，滤波的效果也不同。在某些实际应用的平滑程序中使用了不同于式（5.1.9）的权因子，例如在SPECTRAN-F3程序中使用的平滑公式如下：三点平滑公式为(5.1.12)五点平滑公式为(5.1.13)七点平滑公式为(5.1.14)这几个平滑公式的优点是权因子都是正数，平滑之后的谱数据不可能出现负值，从而提高了平滑之后的谱数据的可靠性。这在原始谱数据中本底很小、峰很高、而且峰的宽度很窄时是非常重要的。如果平滑之后的谱数据出现了负值（这显然是不合理的），可能使后续的计算程序在运行时产生错误。 2. 采用高斯滤波器的平滑方法如果把谱数据中的统计涨落看成是“白噪声”，当使用匹配滤波器进行滤波时，可以得到最佳的信号噪声比。所谓匹配滤波器，就是该滤波器在道域中冲击响应函数与峰函数互为镜象。一般情况下，谱中的峰函数可以近似为高斯函数。由于高斯函数是偶函数，所以匹配滤波器在道域中的响应函数也应该是高斯函数。实践表明，仅在平滑窗口比较大的情况下，使用高斯滤波器进行平滑处理，才能得到比最小二乘移动平滑方法稍好的平滑效果。因此，目前应用最广泛的还是最小二乘移动平滑方法。3. 谱平滑中的几个具体问题对谱进行平滑处理可以减少谱数据的统计涨落，从而减少了寻峰过程中假峰出现的几率，也可以减小峰净面积的计算误差。但是当滤波器的参数选择不当或平滑次数过多时也会产生某些缺点。例如，在寻峰时可能漏失弱峰，不能分辨距离很近的重峰等等。因此，如何选择滤波器的参数和平滑的重复次数是很重要的。(1) 平滑窗口的选择由公式（5.1.3）可以看出，离散量的卷积运算实际上是加权求和。当计算平滑后的谱的第m的数据时，需要在原始谱中第m点两边各取K个点（共2K+1个点）进行运算。我们把2K+1叫做平滑窗口。改变平滑窗口的大小对平滑效果有很大的影响。图5-1-3中 谱平滑处理的效果与平滑窗口大小的关系曲线图5-1-3中画出了经过平滑处理之后谱中统计涨落的改善与平滑窗口大小的关系曲线4。图中横坐标为最小二乘移动平滑方法的平滑窗口，纵坐标是标志着平滑效果的统计涨落改善因子。这个改善因子代表了平滑后与平滑前的涨落的比较，即峰的之内的面积与本底噪声的比值或峰高与本底噪声的比值的变化。由图5-1-3看出，当平滑窗口比较小时，随着平滑窗口增大，平滑的效果增大。在某一个平滑窗口时，改善因子达到最大值。但是当平滑窗口继续增大时，改善因子反而下降。这是因为在平滑窗口较小时，随着窗口的增大，谱中的统计涨落减小很快，但谱中的峰高和峰的形状变化不大。当平滑窗口超过一定数值时，经过滤波之后峰高将急剧下降，因而平滑效果反而降低。改善因子达到最大值时的平滑窗口称为最佳平滑窗口。最佳平滑窗口的大小与谱曲线中峰的宽度有关。当峰的半高宽FWHM比较大时，最佳平滑窗口也较大。因此，为了达到好的平滑效果，对于宽度不同的峰需要选取不同的平滑窗口。选取平滑窗口大小时应该考虑的另一个因素是平滑窗口大小对谱曲线形状的影响。当平滑窗口比峰的FWHM大很多时，平滑之后的谱中的峰将显著地变宽。这将使谱中原来相互靠得比较近的峰重叠得更加严重，从而使寻峰和峰净而积的计算更加困难。综合上面二个因素，在平滑处理时平滑窗口的大小要根据谱中峰的宽度来选择。一般的作法是选用的平滑窗口近似等于峰的半高度FWHM（以道为单位）。例如，当FWHM7时，取2K+1=5；7FWHM9时，取2K+1=7；9FWHM11时，取2K+1=9，等等。能谱曲线中峰的宽度随道址的增加而加大。我们可以把整个谱分成若干段，每段采用不同的平滑窗口。(2) 平滑重复次数在使用较小的平滑窗口时，对谱数据多次重复地进行平滑处理，可以更有效地减小谱数据中的统计涨落。一个均值为常数、服从正态分布的伪随机数系列，其平滑效果与平滑次数的关系曲线如图5-1-4所示。横坐标MS为平滑的重复次数，纵坐标是平滑之后该数列分布的标准偏差的相对值。由图中可以看出，随着平滑次数的增加，谱数据的统计涨落逐渐减小。但是在平滑次数大于3时，曲线下降得很平缓，再增加平滑次数对统计涨落的改善并不显著。多次平滑会使谱的形状产生畸变。这种畸变主要表现为峰高降低，峰宽增大，峰谷被填平。作为一个例子，图5-1-5和表5.4中给出了平滑次数对谱中峰形状的影响。图5-1-5是用Si（Li）探测器测得的55Fe的X射线谱。表5.4列出了使用五点最小二乘移动平滑方法时，峰位（mp）、峰高（h1）、峰谷（h2）和峰的关高宽（FWHM）与平滑次数MS的关系。由这个例子我们可以看出，当平滑次数增加时，峰位基本上维持不变，但是峰的半高宽增加，峰高降低，峰谷抬高。对于一个单峰来说，虽然峰高下降，但半高宽增加，因而峰的面积变化不大，不会使谱的定量分析产生很大的误差。在重叠峰的情况下，由于峰的展宽，可以会淹没位于某一强峰附近的弱峰。在复杂谱的分析中，这将会造成弱成分的漏失。图5-1-4 平滑效果与平滑次数的关系曲线 图5-1-5 55Fe的X射线谱表5.4 峰形状与平滑次数的关系MS（次数）mp（道址）h1（计数）h2（计数） FWHM（道数）01017.31192935024.511017.41182335024.551017.61160435224.9101017.61134735625.4301017.61047737527.5501017.6978541529.41001017.6855151833.4总之，需要对谱数据平滑多少次应考虑到改善谱的统计涨落、减少谱形畸变两个因素，根据谱数据的具体情况决定。在谱数据中各道计数较低，统计涨落较大的情况下，平滑次数可以多些。在谱数据中各道计数较大，或者谱形比较复杂的情况下，为了减少谱形的畸变，节省计算时间，平滑次数应当少一些，一般不多于3次。下面给出用FORTRAN语言编写的最小二乘移动平滑法进行谱数据平滑处理的程序。程序中变量的意义如下：SP：实数组，存放平滑前后的谱数据。A：实数组，工作单元暂存临近谱区左右边界的谱数据。B：实数组，存放平滑公式中的权因子。MB：整变量，谱区的左边界道址。MF：整变量，谱区的右边界道址。IP（1）l：IP数组的第一个元素，存放平滑次数。FW：实变量，谱区内的平均FWHM。SUBROUTINE SMTH （MB，MF）DIMENSION A（5），B（6）COMMON/FPDATA/IP（15），RP（10）COMMON/CURRP/FW，SG，CDL，CDR，STEP，NOSP，NPK，NPS，NPF，NN，NMCOMMON/ARRAY/SP（2100），PKL（160），T（128）IF（IP（1）LE） RETURNM11DO 1 J1，5J1MBJ61A（J）SP（J1）IF（FWLT4）GOTO 3IF（FWLT7） GOTO5IF（FWLT9） GOTO7IF（FWLT11） GOTO9B（1）333333B（2）27972B（3）13986B（4）2331B（5）14895B（6）41958J5GOTO 103B（1）5B（2）25J1GOTO 105B（1）375B（2）25B（3）625J2GOTO 107B（1）3125B（2）234375B（3）9375B（4）15625J3GOTO 109B（1）417249B（2）314685B（3） 6993 B（4）1282 5B（5）34965J410DO 40 IMB，MFSB（1）*SP（I）D0 30 L1，JL1L1L2ILL3IL3 S=S（SP（L2）SP（L3）*B（L1）SP（I5）=S40CONTINUEIMF50SP（I）SP（I5）II1IF（IGEMB） GOTO 50DO 60 K=1，5J1MBK660SP（J1）=A（K）IF（M1GEIP（1）RETURNM1M1GOTO 10END二、寻峰在谱数据中精确地计算出各个峰的峰位是能谱分析中的最关键的问题。在谱的定性分析中，只有正确地找到谱中全部峰的位置，才能根据主峰和各验证峰的能量来决定在被测样品中是否存在某种核素。在谱的定量分析中，尤其是用最小二乘法函数拟合进行重峰分析时，一般使用迭代法，峰位作为迭代参数的初值，如果峰位的误差很大，或混入了假峰，漏失了真峰，则会造成迭代次数增多，甚至不收敛，使迭代失败。由于谱结构的复杂和统计涨落的影响，从谱中正确地找到全部存在的峰是比较困难的。尤其是找到位于很高本底上的弱峰，分辨出相互靠得很近的重峰更为困难。谱分析对寻峰方法的基本要求如下：(1) 比较高的重峰分辨能力。能确定相互距离很近的峰的峰位。(2) 能识别弱峰，特别是位于高本底上的弱峰。(3) 假峰出现的几率要小。(4) 不仅能计算出峰位的整数道址，还能计算出峰位的精确值，某些情况下要求峰位的误差小于0.2道。很多作者对寻峰方法进行了研究，提出了很多有效的寻峰方法。有的方法分辨重峰的能力较强；有的方法适于在很高的本底上寻找弱峰；有的方法计算简单，寻峰速度快。但是各种不同的寻峰方法都可以概括为谱变换和峰判定二个步骤。谱变换一般采用线性滤波技术，其目的是减少谱数据的统计涨落，消除本底的影响，突出峰位信息。峰判定主要是根据预先设置的条件来识别真峰和剔除假峰。谱变换采用线性滤波技术，即对谱数据进行如下的线性变换：（m1，2n）(5.1.15)其中，是谱第m道的数据，是变换之后谱的第m道数据，Cj是滤波器单位冲击的响应函数决定的因子，n为谱中的数据点数，2K1为变换窗口。寻峰中使用的数字滤波器平平滑时使用的滤波器不同。平滑处理时使用低通滤波器，而寻峰中一般使用带通滤波器。谱中的本底成份随道址是缓慢变化的，在频域中属于低频成份，不能通过带通滤波器，因而在变换之后的谱中消除了本底成份的影响。此外，由于带通滤波器既能抑制高频噪声又能抑制低频噪声，从而能更有效地减少统计涨落，突出峰信息。图5-1-6中画出了变换前后谱的形状。由曲线（b）可以看出，变换之后的谱数据中的统计涨落并不能完全清除。此外，在康普顿边缘处也可能出现局部极大值。为了避免把较大的统计涨落和康普顿边缘误认为一个峰，剔除不真实的假峰，需要对变换之后的谱中的峰信息进行峰判定。只有峰信息满足判定条件时才确认存在一个真峰。峰判定包括峰高的统计性判定、峰净面积的统计性判定和峰形判定。几种不同的判定条件可以用在同一个寻峰方法中，也可以使用其中的一部分。采用不同的谱变换方法和不同的峰判定条件就形成了不同的寻峰方法。下面讨论几种主要寻峰方法的基本原理，并对它们的性能进行比较。1. 几种常用寻峰方法的基本原理(1) 对称零面积对合法：进行谱变换时使用的数字滤波器的冲击函数Cj满足下列关系式 （j=K，K）(5.1.16)即冲击函数围绕j0点左右对称，且对横轴所包围的面积为零，因而称为对称零面积对合法。由式（5.1.15），（5.1.16）可以看出，如果在峰区范围内本底谱是常数或按直线分布，则在变换之后的谱中本底的贡献为零。变换之后谱的形状完全反映了峰形的变化。下一步使用峰高统计判定条件来确定峰位。沿变换之后的谱进行检索，找出局部极大点，若极大值超过其均方根误差若干倍时，则认为找到了一个峰，该局部极大点所对应的道址就是峰位。沿道址在谱中进行检索，找出一系列的局部极大值，并进行判定就可以找到一系列的峰。上述的峰高判定条件可以写为：(5.1.17)其中TRH为寻峰阈值，是预先给定的常数。当寻峰阈值取得大时，能有效地剔除由统计涨落造成的假峰。但也可能同时漏掉了弱的真峰。当设定较小的寻峰阈值时，能找到全部的强峰和弱峰，但也可能会把由统计涨落造成的假峰当成真峰。因此，寻峰阈值的大小应当根据谱数据的具体情况和物理实验任务来确定。在一般情况下，寻峰阈值可选在25之间。利用匹配滤波器可以得到最佳信号噪声比，导出匹配滤波器寻峰方法。由于谱数据中的峰函数接近于高斯函数，因而把冲击函数为高斯函数的滤波器作为匹配滤波器。Cj的表达式如下：(5.1.18)其中，为高斯峰函数的标准偏差，2K1为变换窗口。实际上，与道址有关。可以利用谱仪系统的FWHM刻度公式求出第m道的FWHM值，由FWHM/2.355计算出第i道的值，在变换窗口内可以认为值是一个常数。图5-1-6 变换前后的谱曲线(a) 谱曲线(b) 经匹配滤波器变换之后的谱曲线图5-1-6中曲线（a）画出了谱曲线中的一个峰，曲线（b）是经匹配滤波器变换后，按式（5.1.17）求出的Rm值与道址m的关系曲线。在峰位处曲线出现局部极大值，而在峰的两侧出现局部极小值。两个极小值之间的距离可以作为峰的宽度信息，大约等于峰的十分之一高全宽度FWTM。下面是用FORTRAN语言编写的匹配滤波器法寻峰子程序，程序中使用的变量名和子程序名的意义如下：SP：实数组，存放谱数据；PKL：实数组，存放由寻峰程序找到的峰位的精确道址；C：实数组，存放匹配滤波器的冲击响应函数值；MB：整变量，寻峰谱区左边界在SP中的下标；MF：整变量，寻峰谱区右边界在SP中的下标；TRH1：实变量，寻峰阈值；FW：峰的关高宽NP3：在谱区内找到的峰的个数；NPS：谱区左边界道址；GETCRP：子程序名，由峰形刻度参数计算道址为X时的FW值；FDPS：子程序名，在某一整道址IPK附近计算出精确峰位PS。SUBROUTINE PKFD3（TRH1，TRH2，TRH3，MB，MF，NP3）DIMENSION C（5）COMMON/CURRP/FW，SG，CDL，CDR，STEP，NOSP，NPK，NPS，NPF，NN，NMCOMMON/FPDATA/IP（15），RP（1）COMMON/PSCB/NCSH，CHPK（15），CSIGMA（15），CCDL（15），CCDR（15），CSTEP（15）COMMON/ARRAY/SP（21），PKL（16），T（128）EQUIVALENCE（T（1），C（1）UMX.5*（MBMF）NPSCALL GETCRP（X，IERRi）FW2.355*SGN1=IFIX（1.1*FW）D65153*FWD2D*DD22*N11DO1J1，N2J1JN11C（J）EXP（J1*J1/D2）SUMSUMC（J）1CONTINUES1FLOAT（N2）SUMSUM/S1DO5J1，N25C（J）C（J）SUMW2.DO 2 IMB，MFDV.DVV.M1IN1M2IN1DO 8JM1，M2J1JM11C1SP（J）*C（J1）DVDVC18DVVDVVC1*C（J1）DVVSQRT（DVV）SDV/DVVIF（W2.LE.S）GOTO 15IF（W2.LE.W1）GOTO 15IF（W2.LE.TRH1）GOTO 15NP3NP31IPKI1CALL FDPS（W1，W2，S，IPK，PS）PKL（NP3）PS15W1W2W2S20CONTINUERETURNENDSUBROUTINE FDPS（A，B，C，K，PS）W12*BAIF（ABS（W1）.LT.1.E1）GOTO 2PSK.5*（CA）/W1GOTO 32PSK3RETURNEND在PKFD3子程序中，首先根据谱区内平均的FWHM值计算匹配滤波器的冲击响应函数值C（J），然后逐道计算变换谱，并找出变换谱中出局部极大值的整数道址。如果在该点由式（5.1.17）计算的Rm值超过了寻峰阈值TRH1，则找到了一个峰，峰位的整数道址为IPK。在IPK附近，由变换谱用二阶差值多项式求出精确峰位PS，并存入PKL数组中。在匹配滤波器法寻峰的计算过程中，需要重复地进行指数运算，因而运算速度较慢。为了简化运算提出了以冲击函数为矩形的滤波器代替匹配滤波器，这就是矩形滤波器法5。图5-1-7 矩形滤波器的冲击响应函数矩形滤波器的冲击函数如图5-1-7所示，Cj的表达式为（5.1.19）其中N为正奇数，等于变换窗口的三分之一。一般N值取近似于峰的FWHM。为了进一步提高寻找弱峰的能力，也可以把Cj取为（5.1.20）和匹配滤波器法相同，用矩形滤波器法寻峰时也是按道址顺序地进行谱交换，在变换后的谱中找出一系列的局部极大值。然后用公式（5.1.17）进行峰高判定。还可以利用变换谱中局部极大值两侧的局部极小点之间的距离进行峰宽判定。当满足判定条件时，确认该峰是一个真峰。该局部极大值所对应的道址就是峰位。矩形滤波器法寻峰的最大优点是计算方法简单，运算速度快，特别适用于运算速度较慢的微型计算机系统。这种方法也能在统计涨落较大的高本底上寻找弱峰。和匹配滤波器法一样，矩形滤波器寻峰方法的重峰分辨能力也比较差。下面是用FORTRAN语言编写的矩形滤波器法寻峰子程序PKFD5。程序中使用的变量名和子程序名与前面PKFD3中使用的变量和子程序名相同。SUBROUTINE PKFD5（TRH1，TRH2，TRH3，MB，MF，NP3）COM/MON/CURRP/FW，SG，CDL，CDR，STEP，NOSP，NPK，NPS，NPF，NN，NMCOMMON/ARRAY/SP（2100），PKL（160），T（128）W20.0DDV0.0XINPSCALL GETCRP（X，IERR1）FW2.355*SGN1IFIX（1.1*FW）MT（N12）/2M1I3*MT2M2I3*MT2DO 2 JM1，M2IF（J.GT.IMT）GOTO 105DVDVSP（J）DDVDDVSP（J）GOTO 201IF（J.GT.IMT1）GOTO 5DVDV2.0*SP（J）DDVDDV4.0*SP（J）2CONTINUEDDVSQRT（DDV）SDV/DVVIF（W2.LE.S）GOTO 50IF（W2.LE.W1） GOTO 50IF（W2.LE.TRH1） GOTO 50NP3NP31IPKI1CALL FDPS（W1，W2，S，IPK，PS）PKL（NP3）PS5W1W2W2S1CONTINUERETURNEND(2) 简单比较法寻峰6；简单比较法寻峰是最直观而又快速的一种求峰方法。在谱数据中，某一道的数据比其邻近的几道大很多时，则认为该道存在一个峰。在计算中边疆检索平滑后的谱数据，如在第m道满足(5.1.21)则第m道附近有一个峰。式中TRH是寻峰阈值。在第m道附近的谱数据中用二阶差值多项式计算出精确峰位。简单比较法适于寻找强单峰，在高本底上寻找弱峰和分辨重峰的能力都比较差。但算法简单，程序运行速度较快。下面给出用FORTRAN语言编写的简单比较法寻峰子程序PKFD4。其中调用了利用二阶差值多项式计算精确峰位的子程序FDPS。SUBROUTINE PKFD4（TRH1，TRH2，TRH3，MB，MF，NP3）COMMON/ARRAY/SP（2100），PKL（16），T（128）MDO 2 IMB，MFIF（SP（I）TRH1*SQRT（SP（I）IF（SP（I+2）.GT.W1）GOTO 20IF（SP（I2）.GT.W1） GOTO 20N1IS1SP（I）DO 10 J1，3，2N2N12JIF（SP(N2).LE.S1）GOTO 10N1N12JS1SP（N1）10CONTINUEIF（N1.EQ.M）GOTO 20MN1NP3NP31CALL FDPS（SP(N11)，SP(N1+1)，N1，PS）PKL（NP3）PS20CONTINUERETURNENDSUBROUTINE FDPS（A，B，C，K，PS）W12*BACIF（ABS(W1)）.LT1.0E10）GOTO 20PSK0.5*(CA)/W1GOTO 3020PSK30RETURNEND(3) 导数法寻峰：如果我们把谱看成是一个连续曲线，谱曲线和它的一、二、三阶导数的图形如图5-1-8所示。由图中可见，在峰顶处一阶导数值由正到负过零点。二阶导数在峰位处出现负的局部极小值。三阶导数在峰位附近由负到正过零点。因而利用谱曲线在峰位附近形状的特征，由谱曲线的斜率，曲率的变化可以准确地定出峰位。值得指出的是，利用谱曲线的曲率变化可以分辨重峰，确定重叠中各个组分峰的峰位。由图5-1-9中可以看出，由于两个组分峰相距很近，叠加而成的重峰只有一个极值。用二阶和三阶导数的变化特征仍然可以准确地定出每个组分峰的峰位。在二个组分峰的峰位处，出现负的极小值，由负到正过零点。因此，用二阶导数法和三阶导数法寻峰能够分辨重叠峰中各个组分峰的峰位。此外，利用峰和康普顿边缘形状的不同，也能很好地把它们区分开。由于上述特点，以及计算方法比较简单，导数法是目前使用最广泛的一种寻峰方法。图5-1-8 谱曲线及其各阶导数的图形 图5-1-9 重叠峰及其各阶导数的图形 一阶导数法导峰7：一阶导数法寻峰的计算步骤如下。首先用公式（5.1.10）计算出各道的平滑后的一阶导数值。按道址顺序检索出一阶导数值并求出由正到负的零点。如果第m道的一阶导数值为正，第m1道的一阶导数值为负，则有一个峰的峰位在m，m1道之间。用一阶导数值的线性插值求出峰位的精确数值。当找到一个峰位时，在此峰位的左侧求出一阶导数值的局部正极值的道址m1，在峰位的右侧求出一阶导数负极小值对应的道址m2。令W2m2m1，W2代表了峰的宽度。只有满足峰宽判定条件。0.8 FWHMW23FWHM(5.1.22)时，才认为该峰可能是一个真峰。W2太小时，该峰可能是由统计涨落造成的假峰。W2太大时，该峰可能是一个康普顿边缘。在这两种情况下，找到的峰都不是真峰，必需剔除。在实际应用中公式（5.1.22）中的系数3和0.8需要适当地进行调整，以适合该测量系统的具体情况。第二个判别条件是峰高判定。峰高判定条件为(5.1.23)式中，mp是一阶导数过零点的道址；ymp是峰位处的谱数据；Bmp是峰位处的本底值；TRH为预先设定的寻峰阈值，一般取24之间。在实际计算中为了避免本底计算的误差，对公式（5.1.23）进行如下的简化：设谱曲线为连续曲线，峰函数为高斯函数，连续谱函数为(5.1.24)式中x为道址，X为峰位，A为峰的净高度，B为本底函数。由于本底函数的斜率很小，谱的一阶导数近似等于而(5.1.25)由(5.1.23)、(5.1.24)、(5.1.25)式可以把峰高判定条件改写为(5.1.26)这个公式就是在一阶导数法寻峰程序中实际应用的峰高判定条件。图5-1-10 一阶导数法峰判定公式中的几个参数的图示在上面讨论的一阶导数法寻峰中，使用了式（5.1.22）作为峰宽判定条件。实际应用中发现这个判定条件在由于测量条件的变化而导致峰的形状改变时常常不易满足，从而造成假峰或漏失真峰。为此文献3中采用了一种先由谱中找出峰区的左、右边界，然后再在该峰区内用一阶导数法导峰的方法。其具体步骤如下：首先确定峰区的左、右边界首址。沿道址增加方向检索谱数据及其一阶导数值，当谱的一阶导数值为正，其值大于某个数值，且随着道址增加，谱的一阶导数值也增大时，该道址mL就是峰区的左边界道址。在mL道应满足下述条件(5.1.27)(5.1.28)(5.1.29)中期，ymL是mL道的谱数据。是mL道平滑谱的一阶导数值（参见式5-1-10和表5.2）。K为常数，是对应于某一置信度时正态分布的横坐标。K值与置信度的关系如表5.5所示。用于寻峰程序中的置信度称为峰置信度，可以预先指定，一般在0.600至0.999之间。当峰置信度选为0.90时，能够剔除假峰，而且能够找到大多数真峰。中mL道一阶导数值的标准偏差。当变换窗口为五点时(5.1.30)表5.5 K值与置信度的关系置信度K0.9993.0900.9952.5760.9902.3260.9751.9600.9501.6450.9001.2820.8000.8420.7500.6750.7000.5250.6000.2540.5000在谱区中沿道址增加的方向检索，在某一道址中，公式（5.1.27）、（5.1.28）、（5.1.29）满足时，此道址就是峰区的左边界道址mL，当道址继续增大，在某一道址中满足公式（5.1.31）、（5.1.32）、（5.1.33）时，认为该道址是峰区的右边界道址mR(5.1.27)(5.1.28)(5.1.29)在峰区的左、右边界确定之后，需要检验峰区的宽度。当下述公式(5.1.34)满足时，认为是一个正常峰区，否则认为峰区太窄。而式（5.1.34）中W是预先估计的最小峰区宽度。峰区确定之后，在mL至mR之间对谱数据的平滑一阶导数值进行检索。当一阶导数值由正变负时，则在该道附近存在一个峰。在该道附近对一阶导数值进行线性内插，求出一阶导数值过零点的精确位置就可以计算出精确峰位。下面给出在峰区内用一阶导数法寻峰的子程序PKFD1。在程序中，实数组SP中存放原始谱数据，实数组DIF1中存放谱数据各点的平滑一阶导数值。为了节省内存，另外设置了一个长度为256的数组SPA。首先调用子程序GETSP把谱数据由SP分段地传送到SPA中。再调用子程序DIFER1，计算这一段谱的平滑一阶导数值，计算结果存放在DEF1中。由SPA和DIF1计算峰区的左、右边界（IPKSTA、IPKEND），在峰区内判定一阶导数值由正到负的零点，并用线性内插法计算精确峰位，精确峰位的计算结果存放在PKL中。程序中使用了两个寻峰阈值：TRH1用于决定峰区的峰置信度，由操作员根据谱中峰的形状指定。SUBROUTINE PKFD1（TRH1，TRH2，TRH3，MB，MF，NP3）LOGICAL IWTESTDIMENSION SPA（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)，CSTEP(15)COMMON/ARRAY/SP（2100），PKL（160），T（128）FSIG（M1）1*SQPT(4*(SPA(M1+2)SPA(M1-2)SPA(M1+1)SPA(M1-1)IGMBNSP25610IF(IF#GE#MF) GOTO 700CALL GETSP（SPA，IG，NSP，MF）CALL DIFER1（SPA，DIF1，NSP）NSP1NSP-3IWTRH2I3IWTEST1W.TL.1050DO 300L1I,NSP1IF(DIF1(L1).LE.0) GOTO 300IF(SPA(L1).LE.0.001) GOTO 300IF(SPA(L1+1).LT.SPA(L1) GOTO 300IF(DIF1(L1+1).LE.DIF1(L1)-FSIG(L1)GOTO 300DTTRH1*SQRT(SPA(L1)IF(DIF1(L1).LE.DT) GOTO 300DT1TRH1*SQRT(SPA(L1+1)IF(DIEQ(L1+1).LE.DT1) GOTO 300GOTO 310300CONTINUEGOTO 31010 IPKSTAL1ISTOPNSPIF(.NOT.IWTEST)ISTOPNSP-1320IF(L1+2.GT.ISTOP) GOTO 600M1L1+2330IF(SPA(M1).LE.0.001) GOTO 200IF(DIF1(M1).GE.0.0) GOTO 390SIGF0.0IF(M1.LT.NSP-2)SIGFFSIG(M1)IF(DIF1 (M11).GT.DIF1 (M1)SIGF) GOTO 390IF(IWTEST)GOTO 340IF(DIF1(M1)-SIGF.LE.SQRT(SPR(M1) GOTO 390340IF(-DIF1(M1)-SIGF.LE.SQRT(SPR(M1)GOTO 200390IF(M1.EQ.ISTOP) GOTO 600M1M11GOTO 330200IPKENDM1WIPKEND-IPKSTAIF(W.LE.TRH2) GOTO 230JIPKSTAGOTO 250230IIPKENDGOTO 50250JJ1IF(J.GE.IPKEND) GOTO 230IF(DIF1 (J).GE.0.0) GOTO 250NP3NP31IF(NP3.GT.158) GOTO 700PKL(NP3)J2IGDIF1(J1)/(DIF1(J1)-DIF1(J)IF(PKL(NP3).GT.MF) GOTO 700400JJ1IF(J.GE.IPKEND) GOTO 230IF(DIF1 (J).kE.0.0) GOTO 400GOTO 250500LGNSPIG-10GOTO 10600IF(IPKSTA.LE.100) GOTO 650IGIPKSTALG10，GOTO 10650IGIG100GOTO 10700 RETURNENDSUBROUTINE DIFER1(SPA,DIF1,NSP)DIMENSION SPA(1)，DIF1(1)COMMON/CURRP/FW,SG,CDL,CDR,STEP,NOSP,NPK,NPS,NPF,NN,NMDO 10 I1,NSP10DIF1(I)0.0DIF1(1)SPA(2)SPA(1)DIF1(2)0.5*(SPA(3)SPA(1)DIF1(NSP)SPA(NSP)SPA(NSP1)DIF1(NSP1)0.5*(SPA(NSP)SPA(NSP2)NSP2NSP2IF(FW.LT.4) TOGO 30DO 20 I3,NSP220DIF1(I)0.1*(2.0*SPA(I2)SPA(I1)SPA(I1)SPA(I1)2.0*SPA(I2)RETURN30DO 40 I3,NSP240DIF1(I)0.5(SPA(I1)SPA(I1)RETURNEND一阶导数寻峰方法是最常用的一种寻峰方法。用这种方法能找出大部分单峰。计算方法较简单，运算速度比较快。其主要缺点是不能分辨相距很近的重叠峰。由图5-1-9可以看出，当两个峰相距很近时，谱曲线只有一个极值，一阶导数中只出现一个由正到负的过零点。在这种情况下，一阶导数寻峰方法把由两个组分峰重叠而成的重峰当成了一个单峰，从而漏失了一个有意义的峰。另一方面，当谱数据中的统计涨落比较大时，使用一阶导数寻峰方法可能找到很多个由统计涨落造成的假峰。在用较小的寻峰阈值寻找弱峰时，出现假峰的现象更为严重。对找到的峰进行净面积判定是降低假峰出现几率的有效方法。当峰的净面积比