最小二乘拟合解谱方法的研究

《最小二乘拟合解谱方法的研究》由会员分享，可在线阅读，更多相关《最小二乘拟合解谱方法的研究（49页珍藏版）》请在装配图网上搜索。

1、成都理工大学毕业设计（论文）最小二乘拟合解谱方法的研究摘 要随着探测器和电子技术、计算机技术的高速发展，核辐射测量仪器越来越多的采用能谱分析技术，它可以一次测量获得更多的有用信息，是快速、可靠地确定待测样品中放出各种能量的放射性核素或元素的性质及其强度的重要手段。几十年来，它在核物理研究、原子能各种领域、核防护、环境监测、资源勘查等方面都发挥了巨大的作用。最小二乘拟合解谱的应用有其独到之处，特别是在重峰分析方面，可以用此方法获得我们肉眼分辨不开的重叠峰，得到重峰的相关信息。本文主要介绍用高斯峰形函数（高斯峰）加本底函数拟合实测谱线，用非线性最小二乘拟合法解出拟合参数的方法，并将此方法形成了程序

2、。程序主要包括单峰、多峰分析和重峰分析。程序分析实际的谱线得到的拟合结果可用图像直接绘出以便直观的看出与原始谱线的符合程度，分析所得的峰高，峰位，半高宽，面积等数据信息可保存在文件中，以便于进行后续的分析和研究。关键词：非线性最小二乘；能谱分析；高斯拟合The Research of Nonlinear Approximation Analysis of SpectrumAbstract：With the rapid development of computer and detector and electronic technology, more and more nuclear rad

3、iation instruments adopt spectrum analyzing technology, which can get more useful information in one time, is a fast and reliable method of determining elements or elemental properties and its intensity, these elements can emit several energy and form particular spectrum. In the past tens of years,

4、it plays a big role in nuclear physics, each field of atomic energy, nuclear protection, environmental monitoring and resource exploration.The application of least square method to analysis spectrum has its original points, especially in multiple peaks analysis, using this method can get weak peaks

5、which were invisible and submerged in other peaks. We can obtain relevant information of submerged peaks by analyzing and fitting peak positions. In the passage, we introduce a method of fitting Gaussian peaks and background and then deriving fitting parameters by least square method, besides, writi

6、ng a program for this method. The program mainly includes Gaussian peaks (single peak and multi-peak) fitting analysis and multiple peaks analysis. The fitting results derived from analyzing actual spectrum line by the program can be drawn directly in the form of pictures so as to detect coincidence

7、 extent between actual spectrum-line and the original spectrum-line, following data messages such as peak height, peak position, FWHM, area from analysis results can be reserved in file so as to Follow-up analysis and research. Key words: nonlinear least squares; spectrum analysis; Gaussian fitting目

8、录第1章前言11.1研究意义11.2国内外研究现状11.3主要工作内容4第2章 谱线的特点52.1核辐射测量的特点52.1.1核辐射是核衰变的产物52.1.2核辐射的能量具有特征性52.1.3核素的含量与特征辐射的强度存在正比关系52.2仪器谱62.2.1 射线和物质的相互作用62.2.2 统计涨落的影响82.2.3 其它作用的影响8第3章 解谱的基本理论93.1 谱数据的获取原理93.2 高斯拟合法的原理103.2.1扣除基底求“净”全能峰103.2.2峰形函数加基线函数同时拟合的峰面积法123.2.3单峰曲线拟合分析方法123.2.4重峰的分析方法133.3 算法介绍133.3.1非线性模

9、型133.3.2梯度和海赛矩阵的计算153.3.3勒温伯格马阔特方法16第4章 分析软件的设计与实现204.1分析软件的功能介绍204.2分析软件界面214.3数据处理流程图及其内容介绍224.3.1 数据处理流程图224.3.2解谱算法的主要函数代码244.3.3主要程序代码28第5章 拟合结果的分析305.1常数基底的分析325.1.1单峰325.1.2多峰（宽谱段）335.1.3重叠峰345.2 一次基底函数的分析355.2.1 单峰365.2.2 多峰（宽谱段）375.2.3重叠峰385.3对比分析及讨论405.3.1常基底与一次本底函数的对比405.3.2拟合法与累计面积法的面积值对

10、比41结论42致谢43参考文献4445第1章 前言1.1研究意义随着探测器和电子技术、计算机技术的高速发展，核辐射测量仪器越来越多的采用能谱分析技术，进行谱数据的采集和分析。能谱分析技术，可以一次测量获得更多的有用信息，是快速、可靠地确定待测样品中放出各种能量辐射的放射性核素或元素的性质及其强度的重要手段，是一种比较直观的仪器分析技术。几十年来，它在核物理研究、原子能各种领域、核防护、环境监测、资源勘查等方面都发挥了巨大的作用。这种分析技术主要包括能谱的测量技术以及数据的收集和分析技术。前者主要是样品的制备和测量方法，否则即为数据的获取和分析方法。然而，由于在测量过程中多种因素的影响，实际测得

11、的谱线是一个复杂的、相互干扰的、连续谱，要从这复杂的谱线中得到更多、更可靠的信息，就需要对谱的分析方法进行研究。本次实习就是研究如何利用非线性最小二乘拟合法，对核辐射测量谱线进行解析，建立合适的拟合函数，将拟合得到的谱线结果与实验谱对比，如果拟合效果良好可利用拟合的数据得到我们需要的峰面积，峰位以及能量等解谱相关的量，并在此基础上研究利用此方法对重峰进行分析求解。该工作获得的成果可直接应用到实际解谱中去，同时也为解谱方法的教学提供一个软件工具。1.2国内外研究现状自从1948年Hofstodter开始使用NaI（TI）晶体探测伽玛射线1，1950年以Wilkinson的模拟-数字转换方法为基础

12、的多道分析器出现以后，由探测器和多道分析器的组合，构成了伽玛能谱数据的获取系统。随着探测技术的不断发展，于六十年代初，又出现了高分辨率的Ge（Li）探测器，这就要求多道分析器增加道数和提高工作速度。这样，由Ge（Li）探侧器和高道数多道分析器构成的数据获取系统，使数据量大大增加，此外，由于放射性测量技术的多样化，一个实验所得到的数据量就相当庞大，使得用以往“硬件”型多道幅度分析器来获取数据受到限制。于是从六十年代开始，使用计算机作为射线能谱的数据获取装置。计算机具有优良的数据处理能力，能够把获取到的大量信息减化为少量必要的信息。它既具有“硬件”型多道幅度分析器的特长，也具有计算机“软件”的主要

13、功能，既是数据获取装置又是数据处理设备，越来越显示出它的生命力。所谓数据处理1，就是要从测量数据中，确定物质的含量。对于能谱分析来说，从能谱测量系统获取的数据中，经过一定的数据处理和分析方法，求出待测试样的定性、定量结果。如果试样的成分已知，就是要求出其中各成分的量。对于未知组成的试样，数据分析包括定性分析和定量分析两步。第一步，利用各种核（元）素的射线能量的不同，以及辐射多种能量射线的核（元）素，射线相对强度比的差异，从预先备有的核（元）素特征中，细致筛选出可能的核（元）素，确定被测试样中所含有的成分。第二步，在定性分析的基础上，根据能谱的特征、复杂程度、标准的放射源谱的具备情况以及数据处理

14、手段等条件，选择恰当的解析方法，或者几种方法结合使用，对能谱进行定量分析。能谱的定量分析方法1，从国外来看，于五十年代末期开始经历这样一个发展过程，峰面积法逐次差引法（剥带法）逆矩阵法（解联立方程组）逐道最小二乘法峰面积法。对于NaITI）谱，在电子计算机没有得到广泛应用以前，人们对于伽玛能谱数据的分析，只能手工计算，或者借助于台式计算器处理。这样，就不可能进行复杂的数学分析，只能进行各道计数累加的简单峰面积法。然而，由于NaI（TI）探测器的分辨率有限，较为复杂的混合放射性核素的能谱，简单的峰面积法就无能为力，促使人们进一步发展更加精确而较复杂的解析方法。随着电子计算机技木的发展和应用，为人

15、们用更复杂的数学分析作能谱更精确的定量分析提供有利条件。因此，逐道最小二乘法就应运而生，而且在精确地解析NaI（TI）谱的发展中，不断地补充和完善起来。半导体探测器的出现，使得能谱探测技术发生了革命性的变化。由于半导体探测器具有很高的分辨率，即使混合核素较为复杂的能谱，大部分的射线谱峰均能孤立分开。这样，又使得能谱的定量很析变得简明，峰面积法又重现了新的生命力。累加计数的峰面积法，固然在许多场合仍然十分有用。但是计算机的数据处理能力，使得峰面积法又有新的发展，出现了更加确切反映峰的特征、复杂分析函数拟合的峰面积法。这种函数拟合法能更准确地求得真正面积。目前NaI（TI）探测器和半导体探测器，由

16、于它们各具优点，根据不同的探测对象择而使用，因而均有广泛应用。九十年代，数据的处理手段主要是两个方面，一方面是应用小型计算机作射线能谱数据的获取，同时利用计算机“软件”的功能作数据的输入、处理和输出，即“在线”运行。并用直观的图形显示法将结果显示出来，可以一边观察呈现在显像管上的图形结果，一边通过计算机的控制机构或者用光笔通过显像管向计算机传达新的运算指令，实现人和机器的直接“对话”。另一方面是应用大型的通用计算机作复杂的能谱数据处理，即将多道分析器测得的能谱数据，通过磁盘、磁带或穿孔纸带送到计算机进行处理，即“离线”运行，对混合的能谱作定性和定量的解析工作，从而给出试样中所含的各种放射性核素

17、的性质和活度。目前能够看到的谱数据处理方面的文章不多，相应的程序就更少了。找到的相近的文章有国外的T.Vidmar7等人写的一篇介绍用单放射性核素谱线性合成为可以尽量匹配测量谱的研究的文章；成都理工大学的王玉双3专门针对野外X荧光仪的分析程序，主要为了满足野外硬件受限的条件下而采用简单的数学拟合的方法来分析谱，以便在野外现场的得到的更多的元素信息。文中用到了非线性最小二乘法求高斯拟合参数的方法，但毕竟受到野外仪器硬件方面的限制，其用法也受到了些影响，难以达到人们的满意程度。另外，从中获知国际原子能机构提供的免费插件AXIL，成了当今X荧光分析的一个重要软件。AXIL现在包括了DOS和Windo

18、ws两个版本的软件，它具有自动解谱、能量标定、含量分析等功能。该软件虽然分析精度高、速度快，但是他是基于操作系统的英文界面软件，使用起来不灵活，操作起来不方便。国内很多大型微机软件都是在这个软件的基础上进行的改进。国外也有多家研究机构开发出基于Windows操作系统的荧光分析软件，如winAxIL、GUPIx等，但是价格昂贵，多与硬件一起搭配使用。所以需要发展适合本国习惯的、算法优良、操作简便的分析软件。在伽玛能谱方面，中国地质大学的刘煌根4针对Na（TI）伽玛能谱仪测量的全谱数据分辨率不高，低能特征峰都重叠在一起，但又需要得到这些重叠的特征峰的信息，论文中用构造三高斯函数为基函数的径向基函数

19、神经网络（RBFNN）来分析能谱重叠峰。在对获得的全能谱原始数据进行适当的预处理后，用带约束调节的三高斯峰函数拟合感兴趣的低能重叠峰，从而获得三高斯峰的九个参数，并用MATLAB编制程序绘制各分高斯峰的图形，检验参数并确定是最好的三高斯模型参数。这是一种新的分析思路有一定的借鉴价值。1.3主要工作内容实习期间，针对研究任务主要进行了如下工作：l 通过前期调研，了解当前最小二乘拟合解谱方法的应用现状及现有软件等，在此基础上完成了系统的设计和工作的安排计划；l 对VC+语言的学习准备，对C语言数值算法的理解；l 对谱分析方法的学习。包括数据的光滑，高斯拟合的理论基础与各种拟合函数和实际符合情况的学

20、习；l 根据系统的设计和工作计划，应用VC+设计了软件框架，编写核数据处理软件显示（谱线高低调整，扩展、收缩），高斯拟合等部分；l 在完成软件开发的基础上进行了初步应用试验，并进行了初步分析；l 完成了论文的编写。第2章 谱线的特点2.1核辐射测量的特点核辐射测量具有以下三个特点。2.1.1核辐射是核衰变的产物原子核自发的发生核结构的变化，由一个元素的原子核转变为另一个元素的原子核，同时伴随放射出粒子或电磁辐射的现象，称为放射性衰变。放出的粒子有：粒子、粒子、辐射、碎片、特征X射线等。2.1.2核辐射的能量具有特征性某核素放出的射线（粒子、辐射、特征X射线）一定具有某种特定能量。所以如果探测到

21、某种能量的辐射，则一定存在某种核素。即：放射性核素与辐射的能量间存在一一对应关系。而射线的能量和形成的脉冲信号的最大幅度成正比，且脉冲的峰值幅度经AD转换，转换成一个与峰值幅度成正比的数字值，即道址ch，所以谱线中的道址对应于射线的能量，也对应于放射性核素。这就是进行谱的定性分析的基础。2.1.3核素的含量与特征辐射的强度存在正比关系从核衰变的规律我们可以知道，核素的量越多，在一定的时间内，发生衰变放出的射线就越多，即核素的含量与特征辐射的强度存在正比关系。而谱数据采集所得的每道的计数，就是记录不同能量射线的个数，即统计各特征能量射线的辐射强度，所以每道的计数对应着核素的量，这就是进行谱的定量

22、分析的基础。由以上特性可知，核辐射所形成的原始特征谱应该是线谱，即只有在有关射线能量对应的道址上才有计数，其它的道址计数为零。如图2-1为137Cs和60Co所形成的辐射特征谱，谱中只有在661keV、1173keV和1332keV三个能量处有计数，其它能量处为零，所以谱图中就是三条线，即线谱。137Cs-661keV计数60Co-1173keV, 1332keV能量（道址)图2-1 137Cs和60Co所形成的辐射特征谱2.2仪器谱然而，经过探测器的探测、仪器的采集、测量，实际谱数据采集所获取的谱线并非如此，而是一条复杂的、具有高斯分布特征的连续谱，我们把这种谱叫仪器谱。这也是我们实际要分析

23、的谱，这种谱线主要由射线与物质的相互作用、统计涨落、以及其它效应等作用的结果。下面以NaI 谱仪为例进行说明。2.2.1 射线和物质的相互作用射线的能量范围一般在0.01MeV10MeV之间，在这范围内，射线与物质的相互作用有三种类型：光子与核外电子的作用，包括以下三种形式：l 瑞利散射弹性散射。l 康普顿散射非弹性散射。l 光电效应束缚电子吸收全部光子能量，并逃离原子形成光电子。光子与核或核外电子库仑场的作用，包括以下三种形式：l 电子对效应形成正负电子对。l 德尔布茹克散射或核的势散射弹性散射。光子与原子核（或单个核子）的作用，包括以下三种形式：l 汤姆逊散射弹性散射。l 共振散射非弹性散

24、射。l 光致分裂或核的光效应光子能量被吸收引起光核反应。在以上各种相互作用形式中，主要是光电效应、康普顿散射和电子对生成三种效应，其它的作用不超过总几率的1%，所以我们主要讨论光电效应、康普顿散射和电子对效应。1.光电效应当光子与物质原子的束缚电子相互作用时，光子把全部能量转交给某个束缚电子，使之发射出去，并具有一定的动能，而光子本身消失，这种过程叫光电效应，发射出的电子叫光电子。假设光子的能量为Er=hv0，电子的结合能为，则光电子的动能Eph为： （2-1）一般光子的能量要比电子的结合能大得多，因此，近似认为光电子的动能等于光子的能量。而光电子的能量全部沉积在闪烁体中，形成脉冲信号。2.康

25、普顿散射在发生这种作用过程中，光子与原子的核外电子发生非弹性散射，入射光子只是损失掉一部分能量，损失掉的能量交给电子，产生反冲电子逸出原子，同时光子本身被散射，运动方向和能量发生变化。散射光子有两种可能性：或逃逸出晶体，或继续在晶体中产生次级效应（光电效应或康普顿效应等）。若是逃逸出晶体，则只有反冲电子在闪烁体中沉积能量；若留在晶体中，则继续发生作用直到能量全部沉积或射线逸出晶体。此时沉积在晶体中的能量等于所有过程产生光电子的能量之和。若所有后面效应产生的电子能量与首次康普顿效应产生的反冲电子能量之和正好等于原始射线能量，则和光电效应产生的脉冲幅度是一样的，相当于光子的全部能量。多次康普顿效应

26、，使光子能量不断降低，形成能量从0- hv0的连续分布，所以由康普顿散射所沉积在晶体中的能量分布在整个小于等于入射光子的能量0- hv0范围内。这也是连续谱产生的根本原因。3.电子对效应当光子能量大于1.02MeV时，就存在形成电子对效应的几率，即光子完全被吸收而产生一正、负电子对。当能量大于1.02MeV的光子通过原子序数较大的物质时，将容易出现这种现象。电子对的动能为：。它们在NaI（Tl）晶体中有三种可能性：湮灭光子能量全部消耗在晶体中，它们的总能量是1.02MeV加到上面讲的电子对的动能中去，等于入射光子的能量。一个湮灭光子其中有一个逃逸出闪烁体，于是沉淀在晶体中的能量比入射光子的能量

27、少0.51MeV，从而形成单逃逸峰。两个湮灭光子全部逃逸出晶体，对应的能量为。2.2.2 统计涨落的影响由于在对射线的探测与采集过程中存在着统计涨落，所以谱线具有高斯分布特征11的形状，即：输出脉冲幅度围绕一平均幅度涨落。因此，单一能量入射射线，通过探测器后，得到一个呈高斯分布的、有一定宽度的谱线，通常我们用“能量分辨率”表示高斯分布曲线的宽度。2.2.3 其它作用的影响另外，还有以下的作用，使所测得的谱线更加复杂。探测器外面的物质，如铅屏蔽、样品材料等物质的康普顿散射、韧致辐射对能谱的贡献（消耗在原子或分子的激发上）。l 能谱中的特征X射线峰；l X射线逃逸峰；l 两个光子能量的“加和”；l

28、 大样品对能谱的影响；l 自然本底对样品能谱的贡献；l 其他。第3章 解谱的基本理论3.1 谱数据的获取原理在进行能谱测量时，首先要使射线的能量转变为电脉冲，并使脉冲的幅度与入射射线的能量成正比。多道ADC（模数转换）把转换成的电脉冲转换成数字量，称为道址ch。但是谱数据获取过程中的ADC和普通用途的ADC有差别。普通的AD转换，所转换的模拟信号都是较宽的信号，和此相比ADC的转换时间相对很短，对于一个信号可以进行多次采样转换，所以一般都是每隔一定的时间对模拟信号进行一次转换，把当前的模拟信号幅度转换成数字量，然后以时间t为横坐标轴，把每个转换的数据按时间先后顺序和值的大小画点，就可得到模拟信

29、号的形态。而核信号都是很窄的信号，一般几个s，和ADC的转换时间相当，甚至比ADC的转换时间还要窄，所以不能像上面那样获取整个信号的形态。但是我们知道，一个脉冲信号对应一个射线，而脉冲信号的峰值幅度和射线的能量成正比，所以每个脉冲信号的峰值幅度才是我们所需要的。因而，我们只需要在信号峰值到达时，对峰值幅度进行A/D转换，从而得到峰值的数字量，该值也对应于射线的能量高低。在核能谱测量中，我们就是把所测能量范围（转换为脉冲信号的幅度）由ADC平分成若干份，每个能量都对应一个存储空间，存放该能量射线的个数。然后开始测量，每来一个射线，即在该射线能量对应的存储单元中取出原计数，加1再放回原单元保存，这

30、样统计每个能量的射线有多少个，即可得到每个能量的射线的计数，该计数对应于产生该能量射线的核素（或元素）的量的多少。测量结束后，以射线的能量（脉冲幅度，AD转换的值-即道址）为横坐标，计数值为纵坐标，把所统计的各能量射线的计数，按能量由低到高依次画点，即可得到一条谱线。地址单元存储器中各单元的数据变化 216 0-1 25 0-1 18 0-117 0-116 0-1-2-3-4-5 9 0-18 0-1图3-1 核能谱测量的ADC数据采集过程随着电子技术的发展，目前已有很高速的AD转换芯片，转换时间可达ns量级，可以把核信号像普通ADC数据采集那样定时进行采集，由于转换时间快，则一个脉冲上可以

31、进行多次采样转换，这样即可把整个脉冲信号的形态都采集下来，然后在计算机中用软件方法获取每个信号的最大幅度值，即对应于射线的能量，再统计不同幅度值的脉冲个数，把每个幅度值（能量）都对应一个存储空间，存放该幅度的脉冲统计个数，即对应射线能量的射线个数。由上可知，核能谱测量所获取的谱数据实际上就是一组按能量统计出来的射线个数的计数值，也就是一个数组，通常用无符号长整形数组来描述（例如：unsigned long data）。3.2 高斯拟合法的原理3.2.1扣除基底求“净”全能峰所谓基底就是叠加于峰区的自然本底、该射线的康普顿部分和更高能量辐射的康普顿连续谱或其它弱辐射峰等计数的总和。 指数基底扣除

32、（在半对数坐标纸上为直线）也是常用的一种基底扣除法。我们用指数曲线去逼进峰区的基底计数分布，以数学的形式1表示为 （3-1）其中为左边界道上光滑后的计数，bi为第i道上的基底计数。当i=时，基底计数与右边界道上光滑后的计数相等，代入（3-1）式得： 则： 代入（3-1） 式得： （3-2）故“净”全能峰第i道的计数为： （3-3）按（3-3）即可求出“净”全能峰。指数基底扣除，对于在半对数坐标上为谷型或峰形的基底，基底计数小于待求峰计数10%时，对面积计算结果的影响还不太显著，随着基底计数的增加，对面积计算结果的影响越来越大，结果很不可靠。因此，如果不加甄别地均用指数函数基底扣除，就会导致错误

33、的结果。如何判别叠加于待求峰上的基线形状，然后采用什么样的函数基底扣除法，还值得很好地去研究。 用二次曲线拟合代替高斯函数的做法太粗糙，尽管很简单。我们必须从真正的高斯函数出发才有可能得到较好的拟合结果。假定谱线在“净”全能峰的形状可以用纯高斯函数来描述为： （3-4）其中为峰址，Cm为峰道址上的计数（峰高），H为峰的半宽度。由（3-4）式求积分10得： （3-5）用方程组求出的参数和代入（3-5）式即得全能峰的面积3.2.2峰形函数加基线函数同时拟合的峰面积法上面所述的峰面积法，确定拟合峰区以后先扣除基底然后对“净”峰进行拟合。这里再介绍确定峰区以后将基线函数与峰形函数结合于一起，和原始峰区

34、数据进行拟合，最后由峰形函数求积分以获得峰面积的方法。峰形函数仍然用纯高斯函数、基线函数还是指数函数，那么描述全能峰的函数1为： （3-6）其中H为峰的半宽度，bl为左边界道上的基底计数，其它符号均与上述公式相同。3.2.3单峰曲线拟合分析方法在上面已经叙述了谱峰函数拟合的峰面积法：一种是先扣除指数基底，然后用纯高斯函数（或改型双曲正割函数）拟合计算峰面积，另一种是将指数基线加纯高斯函数一起进行拟合，然后由高斯函数的参数计算峰面积。实验表明，用简单的高斯函数描述峰形状，有的简单的谱线可以拟合得很好。但对于其他的一些峰，用简单的高斯函数来描述峰的形状，拟合就不理想，尤其是在峰的低能侧，拟合得很不

35、好。当然对于距离较远很孤立的单峰来说，用（3-1）和（3-2）的方法或者用线性基线或低次多项式基线加简单高斯函数进行拟合，仍然可以获得较精确的结果。只要效率刻度谱峰采用相同的拟合方法，一样可以获得比较准确的拟合效果。 然而，在重峰分析中，描述峰的数学模型不恰当，将会影响峰面积的准确计算，从而影响射线辐射强度测定的可靠性。特别是在一个多道谱中，既有单峰也有重叠峰的情况下，更需保证单峰和重叠峰分析精度的一致性。为了获得满意的描述全能峰的数学模型，还得从单峰的曲线拟合分析做起。3.2.4重峰的分析方法 由单峰曲线拟合的分析方法，获得正确的峰形函数和阶梯函数后，即可建立正确的单峰拟合函数1： （3-7

36、）其中基线函数B（i）是由峰区外两边的连续谱确定。对于重叠峰组来说，基线函数同样可由重叠峰区外两边的连续谱确定，通常它取为一个L次多项式： （3-8）其中（k=1.2.L+1）为待定参数。这样，描述重叠峰组的拟合函数为： （3-9）其中，m为重叠峰组内峰的数目，j为组内峰的序号。重叠峰区的拟合函数为（不考虑阶梯函数时） （3-10）3.3 算法介绍3.3.1非线性模型我们现在考虑，模型和m个未知参数集 （kl，2M）之间是非线性关系的拟合。我们定义优值函数通过求它的最小值确定最佳拟合参数。但是在非线性关系时，求最小值必须迭代进行。对参数给定一个初始值，我们设计个改善这个初始实验解的过程。这个过

37、程接着不断重复进行，直到值不再增长（或者不再明显增长）为止。这个问题和我们已经解决的求一般非线性函数极小化的问题有何不同呢？表面上有不同，但不是完全不同：当非常接近极小值时，我们希望函数趋近于二次型2 5，我们可以将它写为： （3-11）其中，d是长度为M的向量，而D是M×M的矩陈。如果这种近似是个很好的近似，则我们知道怎样从当前的实验参数值一步跳到极小值的也就是： （3-12）另一方面，式（3-11）对于在处进行极小比的函数形状是一个比较差的局域近似。在这种情况下，我们首先要做的是，沿负梯度的方向走步，如同在最速下降法中那样，换言之，即： （3-13）其中常数要足够小。以保证不偏离

38、下降方向。 为了利用式（3-12）或式（3-13），我们必须能够计算在任意参数集a上函数的梯度。为了利用式（3-12），我们也必须计算矩阵D，它是在a取任意值时，优值函数的二阶导数矩阵（海赛矩阵）。在这里我们没有直接计算海赛矩阵的方法。我们只能估算需要极小比的函数以及（在某些情况下）它的梯度。因此我们必须求助于迭代方法，不仅仅是因为我们的函数是非线性，而且是为了我们要构造海赛矩阵的内容。在这里，情况要简单得多。我们精确地知道的形式，因为它建立在我们自己已经规定的模型函数的基础上因此海赛（Hessian）矩阵对我们来说是已知的。只要我们愿意，我们便可无约束地采用等式（3-12）。使用式（3-13

39、）的唯一原因在于，式（3-12）不能成功地改善拟合，并标志着式（3-11）作为一个局域近似的不可行。3.3.2梯度和海赛矩阵的计算需要拟合的模型2 5是： （3-14） 优值函数是： （3-15）对参数a的梯度，在取极小值时为0，它的各个分量是：k=1, .2, M （3-16）再求一次导数得到： （3-17）可以通过如下定义而去掉因子2： （3-18）和式（3-12）中的D的关系是，利用式（3-18）可以把式（3-13）重写成一组线性方程组： （3-19）这个方程组解出增量和当前的近似相加给出下一步的近似。在最小二乘方问题的文献中，短阵a通常被称为曲率矩阵，等于海赛矩阵的二分之一。等式（3-

40、13），即最速下降等式，可写为： （3-20）注意，海赛矩阵（3-17）的分量 依赖于基函数对它们参数的一阶导数以及二阶导数：某些处理过程常常不加评注地忽略了二阶导数。我们也忽略二阶导数，但是稍加评注后才这样做。 当式（3-16）中的梯度依赖于导数时，二阶导数便发生了、因为对梯度的微分必定包括像这里的项。当二阶微分项是0时自然可以去掉（例如在（3-18）的线性情况下），或者和一阶导数项相比足够小以致可以忽略。在实际应用中还有另外一个因为小而可忽略的可能情况：在方程（3-17）中。乘以二阶导数的项是。对于一个成功的模型来说，这项应当是每一点的随机测量误差。这个误差可正可负，并且一般应该和模型无关

41、。因此，当对i累加求和时，二阶导数项就相互抵消了。在拟合模型比较差，或者在模型数据中，不存在通过反号补偿而抵消的数据点时若把二阶导数项包括进来，则会导致不稳定性。从现在起。我们将一直用下式作为的定义： （3-21）3.3.3勒温伯格马阔特方法1.两个重要的洞察 在勒温伯格（Levenberg）的一个早期建议的启示下，马闹特（Marquardt）提出了一个非常巧妙的方法2 5 6，它是在逆海赛矩阵方法（3-19）的极限和最速下降法（3-20）之间进行平滑调和。后一方法在远离最小值时运用，当接近最小值时它逐渐切换到前方法，这种勒温伯格，马阔特方法（也称马阔特方法）在实际应用中非常奏效，并且已经成为

42、非线性最小二乘方问题的标准。 这种方法基于两个基本的但是非常重要的洞察，考虑方程（3-20）中的“常数”，它应该是什么？大小数量级有多大？由什么来确定它的尺度大小？在梯度中没有关于这些答案的任何信息。梯度只说明了一个坡度而没有说明这个坡度延伸到多远。马阔特首先洞察到。海赛矩阵的各个分量，尽管在一些精确场合是无用处的，但却对问题的尺度大小及数量级提供了一些信息。 量是无量纲的，也就是说它是一个纯粹的数；这可以从它的定义式（3-20）中明显看出。另一方面，是由量度的，它很可能有量钢，例如为，或者千瓦小时。或者其它形式的单位（实际上，的各个分量可能有不同的量纲）。因此在和之间的比例常数必须具有量纲。

43、检查一下的各个分量，会发现只有一个量明显地具有这种量纲，那就是即对角元的倒数。因此这些对角元确定了“常数”的尺度。但是这种尺度可能本身就很大，因此我们引入个人为的因子，而将它除以常数无量钢），并尽可能置以减少步长。换言之，用下式代替方程（3-20）： （3-22）上式要求是正的，这可以由定义式（3-21）保证一这也是采用式（3-21）的另一个原因。马阔特的第二个洞察在于等式（3-22）和（3-19）可以结合起来，如果我们按照如下规定定义一个新矩阵： （3-23）那么式（3-22）和（3-19）均可用下式代替： （3-24）当非常大时，矩阵强制成以对角元为主的，此时等式（3-24）和（3-22）

44、几乎等价另一方面，当趋近于零时，等式（3-24）便回到了式（3-22）。2.算法的具体过程 对拟合参数a给定初始值后，我们建议采用如下的马阔特步骤： 计算。 取一个适中的值，比如 。解线性方程（3.25）得到的，并且计算 如果，则将扩大10倍（或者其它的具有意义的倍数），再回到（）。如果，将减小10倍，将代替a，再回到（）。 我们还必须知道终止迭代的条件。收敛后继续进行选代是不必要的，也是浪费的，因为求最小值充其量是对参数n的统计估计（收敛的程度决定于机器的精度以及舍入误差的极限）。当的变化量1时，在统计的角度上参数的变化是毫无意义的。 进一步，我们经常发现极小值位于一个复杂的拓扑形状的平谷处

45、，参数围绕此极小值徘徊。其原因在于马阔待方法是对正规方程的推广，因此它也遇到和正规方程样的问题，即在极小值处的退化。由于零主元也有可能导致的直接失败，但可能性不大。更经常遇到的情况是，一个小的主元将会产生一个大的补值而被排斥，的值随之增加。对于足够大的，矩阵完全是正定的，不可能有很小的主元。因此这种方法的确能倾向于远离零主元，但付出的代价是，在一个下降坡度不大的谷地处做最速下降，从而使趋近于最小值比较慢。 从这些考虑我们建议，在实际中，在减少很小量的第一次或第二次时停止迭代，例如绝对值的减小量小于0.1，或者（以防舍入误差阻止获得这个精度）某些分数值如0.001。不要在增加的那一步停止迭代；这

46、只说明还没有调整得非常合适。一旦可以接受的极小值找到了，就可以令，并计算矩阵： （3-25）如同前而所讨论的一样，上式可用来计算拟合参数a的协方差矩阵（参看下节）。 下面的一对函数程序准确地体现了非线性参数估计的马阔持方法。特别要注意的是，在数组中必须输入分量值为1和0，它们分别对应于数组中相应的参数值是需要被拟合的，还是需要保持固定在它们的输入值。 程序mrqmin执行马阔特方法的一个迭代。在第一次调用它时令alamda<0,它标志着程序的开始。alamda首先设定后，在这个设定值下的调用将作为下次迭代的值；作为迄今找到的最佳估计参数以及对应的返回。当你对收敛程度满意时在进行最后一次调

47、用前令alamda为0，短阵alpha和covar（在所有先前的调用中它们被用作工作空间），就是收敛参数值的曲率矩阵和协方差。参量alpha、a和chisq在两次调用之间不能变化。 alamde除了在最后一次调用时设为零外。任何其它两次调用之间也保持不变。当遇到使增加的一个步骤时，chisq和a返回它们的输入值（最佳值）。但是alamda将设置为一个增大的值返回。 程序mrqmin要调用程序mrqcof用来计算矩阵 （参见（3-21）式，和向量（3-16）和等式（3-18）。在 mrqcof 中又要调用用户提供的程序funcs（x，a，y，dyda）。它用来计算输入值：及时的模型函数值以及微分

48、向量。第4章 分析软件的设计与实现4.1分析软件的功能介绍此软件实现的基本功能如下： 1）读入测量仪器获得的原始数据文件，并分上下窗口显示其直观的谱数据，下窗口显示所有1024道址的数据，上窗口则显示经过处理后的谱线。 2）随光标移动自动计算和显示对应的道址、计数。 3）窗口中的谱线可以放大缩小的显示、当谱线放缩时是以光标作为基准以便显示适合观察的谱线。 4）上窗口中根据不同数据自动生成道址计数刻度线，随着数据的不同显示刻度线跟随变化。 5）上窗口中的数据显示，包括光滑后的谱线以绿色显示，以蓝色显示拟合得到的数据。 6）对原始数据的谱光滑处理，由于本设计主要研究高斯拟合解谱所以借用了其他人做的

49、最小二乘拟合或滑动褶积变换法作为本设计的数据光滑为数据拟合做初始变换工作。 7）按使用的意图分析拟合指定的道址，并提供了不同的拟合公式形式，可由使用自己设定拟合的初始参数以便更好的得到拟合结果，也可直接确定使用程序默认提供的初始参数。 8）针对光滑后的数据或未光滑的数据进行解谱运算，如若成功解谱，则在上窗口显示高斯拟合数据、本底数据。 9）软件具有分解重峰的功能，可根据使用者的判断选定道址区域进行重峰分析，若拟合成功显示拟合谱线。10）生成解谱结果14 15：显示谱线中所包含峰位道址、峰位计数、峰位面积，并且可以以文本.txt格式保存。本核数据处理软件到总体功能结构图如4-1图所示：核分析软件

50、谱线光滑高斯拟合分析高斯重峰分析寻峰谱线的显示图4-1 分析软件功能结构图4.2分析软件界面程序的初始启动界面8 10 12 13如下图4-2所示，程序菜单栏分为:文件、显示、谱分析、查看、窗口、帮助。其中，文件菜单内的子菜单分为：打开、关闭、退出。显示菜单内的子菜单分为：显示谱段的选择、谱线显示的调整、光标的位置的控制、谱线的扩展、收缩。谱分析菜单内的子菜单分为：能量刻度、谱光滑、寻峰、计算峰面积、高斯拟合分析、高斯重峰分析。帮助菜单内的子菜单分为：关于。工具栏很简单，只有新建、打开数据、左右移动光标按钮、上下放大缩小谱线按钮及帮助按钮。主界面的上窗口显示分析过后的谱线，下窗口显示原始谱线。

51、状态栏显示在上窗口中光标移动处的道址、计数。图4-2 软件运行的初始界面4.3数据处理流程图及其内容介绍4.3.1 数据处理流程图1.高斯拟合分析流程高斯拟合分析所实现的数据处理流程图如图4-3所示，其中矩形表示流程过程：菱形表示在某一过程中需要判断；每一流程有相应的箭头指示。此数据处理流程图的内容介绍如下：1）读取需要分析的核数据，放入指定的数组为分析处理做准备。2）对数据进行光滑处理消除统计涨落对分析的影响。根据实际分析的需要也可以选择直接对数据进行高斯拟合处理得到我们想要的结果。3）如果我们选择了高斯拟合分析，则直接进入如图4-3所示的流程中，我们只需要按对话框提示输入相应的信息，完成分

52、析过程。4）最后以图形的方式直观到显示拟合到结果，同时拟合出来的参数可以按自己意愿保存起来以供后续分析。分析数据到导入输入分析到道址自动获得拟合的峰数输入拟合初始参数参数求解高斯拟合结果显示是否保存结果返回保存结果YYNN图4-3 高斯拟合流程图2.高斯重峰分析流程高斯重峰分析所实现的数据处理流程图如图4-4所示，其中有很多地方与高斯拟合分析是相似的，主要的区别就在于我们要手动输入需要用多少个峰来拟合我们分析的谱线和需要手动输入确定我们预先知道的峰位信息，以便我们得到淹没峰的拟合。分析数据到导入输入分析的道址输入拟合的峰个数输入拟合初始参数参数求解高斯重峰分析结果显示是否保存结果返回保存结果Y

53、YNN图4-4 高斯重峰分析流程4.3.2解谱算法的主要函数代码1. mrqmin函数（马阔特步骤）void CSmcaView:mrqmin（float *x, float *y, float *sig, int ndata, float *a, int *ia, int ma, float covar120120, float alpha120120, float *chisq, float *alamda） int j,k,l,m; static int mfit; /*mfit 是一个静态临时变量*/ static float ochisq,*atry,*beta,*da,*oneda

54、; /*中间过度变量*/ if （*alamda < 0.0） atry=new floatma+1; beta=new floatma+1; da=new floatma+1; for （mfit=0,j=1;j<=ma;j+） if （iaj） mfit+;/*oneda是一个临时数组*/ oneda=new float*mfit+1; /为指针数组的每个元素分配一个数组 for （int i = 1; i <=mfit; i+） onedai = new float2; *alamda=（float）0.001; mrqcof（x,y,sig,ndata,a,ia,ma

55、,alpha,beta,chisq）;/*回给得到alpha和beta*/ ochisq=（*chisq）; /*chisq X平方的值*/ for （j=1;j<=ma;j+） atryj=aj; /*将aj给过渡变量*/ for （j=0,l=1;l<=ma;l+） /*通过增加对角元，改变线性化的拟合矩阵*/ if （ial） for （j+,k=0,m=1;m<=ma;m+） if （iam） k+; covarjk=alphajk; /*将alpha变换后给covar 再用它调mrqcof*/ covarjj=alphajj*（1.0f+（*alamda）; one

56、daj1=betaj; gaussj（covar,mfit,oneda,1）; /*矩阵求解 解等式得到增量da 与当前近似相加给出下步迭代*/ for （j=1;j<=mfit;j+） daj=onedaj1; if （*alamda = 0.0） /*收敛*/ covsrt（covar,ma,ia,mfit）/*将协方差矩阵排成拟合系数的真实顺序*/return; for （j=0,l=1;l<=ma;l+） /*试验成功了吗？*/ if （ial） atryl=al+da+j; /*相当于把A和da相加 给A进行变换*/ mrqcof（x,y,sig,ndata,atry,i

57、a,ma,covar,da,chisq）; /*这次调用 没回来 错误*/ /*是不是声明的问题？*/ if （*chisq < ochisq） /*成功了 则接收新的解*/ *alamda *=（float）0.1; ochisq=（*chisq）; for （j=0,l=1;l<=ma;l+） if （ial） for （j+,k=0,m=1;m<=ma;m+） if （iam） k+; alphajk=covarjk; betaj=daj; /*daj来自于betaj，回赋给betaj*/ al=atryl; /*将调整了的A值回赋给A*/ else /*失败了 则增加*alamda和返回*/ *alamda *=（float） 10.0; *chisq=ochisq; 在以上的步骤实现的过程中又调用了其他重要的函数如下所述。2. mrqcof函数此函数被mrqmin 用来计算线性化了的拟合矩阵alpha 和向量beta 以及。void CSmcaView:mrqcof（float *x, float *y, float *sig, int ndata, float *a, int *ia, int ma, float alpha120120, float *beta, float *chis