马英杰， 周 靖， 洪 旭， 周建斌， 王 敏， 万文杰

成都理工大学核技术与自动化工程学院， 四川 成都 610059

核能谱单能峰快速高精度高斯函数拟合算法研究

马英杰， 周 靖*， 洪 旭， 周建斌， 王 敏， 万文杰

成都理工大学核技术与自动化工程学院， 四川 成都 610059

在核能谱分析中， 高斯函数最小二乘拟合法是计算单能峰净峰面积常用的方法， 该方法精度较高， 但噪声敏感性较强， 导致拟合出的高斯函数在峰位附近的残差向量较大。 针对该问题， 对高斯函数最小二乘拟合法进行了详细推导， 分析了峰位附近残差向量较大的原因， 提出了一种基于高斯函数最小二乘拟合法的高斯函数加权最小二乘拟合法， 即在高斯函数最小二乘拟合法的基础上， 引入了权重因子。 该权重因子与取对数后数据权重削弱趋势相反或与数据本身趋势相符， 以减小噪声敏感性。 由于在求解高斯函数参数的过程中涉及到求逆矩阵运算， 计算量较大， 耗时较长， 为了提高实时性， 将求逆矩阵的运算过程转换为了简单的方程组运算， 并给出了高斯函数的幅值、 中心及方差参数的快速求解公式。 将这两种方法用于55Fe的特征X射线单能峰的实际拟合中， 结果表明， 高斯函数加权最小二乘拟合法效果均较好， 这说明该方法降低了噪声敏感性， 减小了高斯函数在峰位附近的残差向量， 进一步提高了拟合精度。 另外， 使用快速求解公式， 也减小了运算量， 增强了实时性， 为在便携式设备中的有效使用提供了可能。

核能谱； 单能峰； 权重因子； 加权最小二乘拟合

引 言

在核能谱分析中， 确定净峰面积的方法基本上可分为两类： (1)计数相加法， 即把峰内的各道计数按照一定的公式直接相加。 这种方法比较简单， 只适用于确定单能峰面积； (2)高斯函数最小二乘拟合法， 即将峰内各道计数拟合为一个高斯函数， 然后对这个函数进行积分， 从而得到净峰面积， 这种方法较计数相加法， 结果更准确， 也适用于重叠峰[1-4]。

在进行高斯函数最小二乘拟合时， 高斯函数的幅值、 中心及方差参数求解的准确性决定了函数拟合的好坏， 是影响解谱效果的关键因素。 高斯函数参数求解的方法有： (1)文献[5-9]使用的LM算法(levenberg-marquardt算法)具有速度快， 不易发散的特点， 但涉及到函数初始值的选取； (2)采用矩阵， 对参数进行求解， 但需要进行矩阵求逆运算， 计算量较大。 此外， 采用高斯函数最小二乘拟合法得到的幅值、 中心及方差对噪声有很强的敏感性[10]。 文献[10]通过高斯函数的第二特征函数得到高斯函数的幅值、 中心及方差， 该方法具有求解速度快， 而且对噪声不敏感的特点。

对单能峰的高斯函数最小二乘拟合法进行研究， 旨在求解高斯拟合参数、 提高拟合精度、 降低噪声敏感性、 增强拟合实时性等方面进行改进与提高。

1 高斯函数拟合原理

射线与物质作用产生电离和激发的统计涨落以及谱仪固有的涨落特性， 使得探测器的输出脉冲幅度围绕某均值小幅涨落[11-12]。 若能较好的扣除本底， 则能谱单能峰的净计数近似服从高斯分布， 即在谱峰区内各道计数y(x)与道址x的关系为

(1)

式中， 变量x代表道址或能量，μ是峰高对应的道址或能量，σ是表征峰形宽度的特征量， 它与半高宽度(FWHM)的关系为： FWHM=2.354 28σ。 若参数y0和σ已知， 则谱峰面积A可以通过积分计算得到

(2)

高斯函数最小二乘拟合的目的就是求出参数y0，μ和σ， 下面对求解过程进行详细推导。 对式(1)两边同时取对数得

(3)

将式(3)展开， 得

(4)

令

则式(4)可表示为

(5)

通过上面的变换， 将一个超越函数变为了简单的二次函数， 从而将高斯函数拟合转换为了二次多项式函数的最小二乘拟合。 若能解出系数c1，c2和c3， 便可根据式(6)求出高斯函数的参数。

(6)

下面对系数c1，c2和c3的求解过程进行推导。 式(5)的实质是一个超定方程组

(7)

将其简写为

Ac=b

其中，

(8)

由于是超定方程组， 所以采用最小二乘法来求解。 设残差向量表示为

r=b-Ac

(9)

bTb-bTAc-cTATb+cTATAc=

bTb-2bTAc+cTATAc

(10)

(11)

要使ρ的最小值存在， 必须有-2ATb+2ATAc=0， 即

ATAc=ATb

(12)

M=ATA，N=ATb

(13)

则式(12)可表示为

Mc=N

(14)

根据式(8)， 式(13)可等价表示为

(15)

根据式(14)和式(15)， 求解三元一次方程组， 即可求得系数c1，c2和c3， 再由式(6)解得σ2，μ，y0， 便可快速得到高斯函数的参数。 这样的运算过程避免了矩阵求逆的大量运算， 简化了运算过程， 提高了拟合速度， 算法也更为简洁， 更有利于在便携式设备中实现。

采用上述步骤对高斯函数Y：y=400exp(-(x-50)2/(2×102))进行最小二乘拟合， 结果如图1所示(为了便于比较， 后面模拟使用的高斯函数均为Y， 数据用折线画出)。

图1 最小二乘拟合结果

将最小二乘拟合结果作线性分析， 结果如图2所示。

图2 最小二乘拟合结果线性分析

由图2得， 最小二乘拟合结果与原高斯函数值之间的线性相关系数r为1， 说明按照上面的方法进行高斯函数的最小二乘拟合的结果是可信的。

实际的能谱中包含有噪声， 这里以高斯函数Y与随机噪声叠加后的结果近似表示一个能谱单能峰， 按照上面的拟合步骤， 拟合结果如图3中曲线b1、 图4中曲线b1所示。

本研究使用构式搭配分析软件Coll, analysis 3.2a在R语言环境下进行运算，之后进行Fisher精确检验(Fisher exact test)，统计出槽位中的动词与目标构式的关联强度。再使用WordNet 2.1进行语义分析。R语言是成熟的编程语言，具有许多实用的程序包，如amap，cluster，fpc等，既能够灵活地进行数据检索分析，又能提供个性化的数据统计，此外还具有绘图功能[37]。

图3 单能峰最小二乘拟合结果一

图4 单能峰最小二乘拟合结果二

综上， 最小二乘拟合结果对噪声有很强的敏感性， 而实际的能谱中往往包含有噪声， 因此减小噪声对高斯拟合的影响， 具有实际意义。

2 高斯函数最小二乘拟合的改进

为解决上述问题， 在拟合公式中为残差向量r添加一个权重因子W。 这里的权重因子应该与取对数运算后数据权重削弱趋势相反或者与数据本身趋势相符合。 这里， 取数据本身(取对数前的值)的归一化值为权重因子， 即

(16)

结合式(9)， 设

s=W(b-Ac)=Wr

(17)

那么，

bTWTWb-bTWTWAc-cTATWTWb+cTATWAc

(18)

将ρ对向量c求导得

(19)

ATWTWAc=ATWTWb

(20)

同样， 为了避免矩阵求逆运算， 令

ATWTWA，Q=ATWTWb

(21)

则式(20)可表示为

Pc=Q

(22)

结合式(8)， 式(16)， 式(21)可等价表示为

(23)

根据式(22)和式(23)求解三元一次方程组， 即可求得系数c1， c2和c3， 再由式(6)解得σ2， μ， y0， 便可快速得到加权的高斯拟合函数。 对图3、 图4中的单能峰分别进行加权最小二乘拟合， 拟合结果如图3中曲线c1、 图4中曲线c2所示。

对比图3和图4中拟合结果， 可以看出， 采用加权最小二乘拟合后， 在峰位附近的残差向量减小， 拟合效果得到明显改善， 说明引入权重因子后， 降低了噪声敏感性。

3 实验分析

实测55Fe标样的X射线能谱， 采用5点平滑进行谱光滑， 使用SNIP方法扣除本底后， 分别采用最小二乘拟合法以及加权最小二乘拟合法的高斯函数拟合结果， 分别如图5中曲线b、 曲线c所示。

图5 55Fe全能峰拟合结果

将图5中道址826—836的拟合结果进行局部放大， 结果如图6所示。

图6 拟合结果局部放大

由图6可知， 采用加权最小二乘拟合后， 在峰位附近的残差向量较采用最小二乘拟合的残差向量小。 采用优良指数AIFOM(analytic improved figure of merit)对拟合效果的优劣(拟合优度)按照式(24)进行评价：

(24)

式中np为拟合峰所包含的道址；n为本底和峰的道数之和；Ap为全能峰净计数； Δy为第i道的实测值与拟合值之差。 判断拟合效果好坏的标准是： AIFOM<0.01%， 拟合效果较好； 0.010.05%， 拟合效果差。

对55Fe标样的X射线能谱进行多次实测， 对全能峰进行拟合， 拟合优度以及净峰面积如表1所示(图5中的55Fe实测谱的拟合优度及净峰面积见实测谱线a)。

表1 拟合优度、 净峰面积对比表

从表1可以看出， 5次测量中， 采用最小二乘拟合时， 拟合优度有四次都大于0.01%， 拟合效果较差； 而采用加权最小二乘拟合后， 拟合优度均小于0.01%， 拟合效果较好。 对同一个全能峰来说， 采用加权最小二乘拟合后， 拟合优度有所降低， 拟合效果明显改善， 说明拟合值越接近实际值。 比较加权前后拟合的净峰面积， 变化率最高达2.76%， 最小为1.01%， 这说明净峰面积前后的变化还是比较大的。 在高放射性场合中， 净峰面积的变化可能会更大。

综上， 与最小二乘拟合法相比， 采用加权最小二乘拟合后， 降低了噪声敏感性， 提高了高斯函数的拟合精度， 从而净峰面积也更为准确。

4 结束语

对单能峰的高斯函数最小二乘拟合的过程进行了详细推导， 在求解高斯函数的幅值、 中心及方差参数时， 将矩阵求逆运算转换为了简单的方程组运算， 并给出了参数快速求解的具体公式， 大大简化了运算过程， 提高了拟合速度。 分析了噪声对高斯函数最小二乘拟合法的影响， 并引入了权重因子对该算法加以改进。 实例证明， 改进后的算法(加权最小二乘拟合法)减小了单能峰峰位附近的残差向量， 降低了噪声敏感性， 提高了高斯函数拟合精度。 结合快速求解公式， 加权最小二乘拟合法是一种拟合速度快且性能较优良的高斯峰面积计算方法。

[1] PANG Ju-feng(庞巨丰). γ Spectrum Data Analysis(γ能谱数据分析). Xi’an： Shaanxi Science and Technology Press(西安： 陕西科学技术出版社), 1990, 681.

[2] QI Rong, MAO Yong, CHEN Xi-meng(齐 荣, 毛 永, 陈熙萌). Nuclear Techniques(核技术), 2008, 31(5): 330.

[3] Fu Chen, Wang Nanping. Nuclear Science and Techniques, 2010. 21(4): 214.

[4] Fudan University, Tsinghua University, Peking University(复旦大学, 清华大学, 北京大学). The Experimental Method of Nuclear Physics, Volume 1(原子核物理实验方法, 上册). Beijing: Atomic Press(北京: 原子能出版社), 1985. 374.

[5] Marquardt D W. Journal of Society Industry Apply Mathematics, 1963, 11(6): 431.

[6] Department of Computation Mathematics of Tongji University(同济大学计算数学教研室). Modern Numerical Mathematics and Computing(现代数值数学和计算). Shanghai: Tongji University Press(上海: 同济大学出版社), 2004. 56.

[7] XIE Zheng, LI Jian-hua, TANG Ze-ying(谢 政, 李建华, 汤泽潆). Nonlinear Optimization(非线性最优化). 2nd ed(第2版). Changsha: National University of Defence Technology Press(长沙: 国防科技大学出版社), 2003. 1.

[8] YANG Zhi-hui, LI Yue-zhong, YU Qian, et al(杨志辉, 李跃忠, 余 倩， 等). Nuclear Electronics & Detection Technology(核电子学与探测技术), 2013, 33(10): 1271.

[9] Morhac M, Matousek V. Applied Spectroscopy, 2008, 62(1): 91.

[10] GU Min, GE Liang-quan(顾 民, 葛良全). Nuclear Techniques(核技术)， 2009， 32(11): 864.

[11] WANG Shao-shun(王韶舜). Nuclear and Particle Physics Experimental Methods(核与粒子物理实验方法). Beijing: Atomic Press(北京: 原子能出版社), 1989. 208.

[12] WU Zhi-hua(吴治华). Nuclear Physics Experimental Methods(原子核物理实验方法). Beijing: Atomic Press(北京: 原子能出版社), 1997.

*Corresponding author

Study on the High Speed and Precision Gaussian Function Fitting Algorithm for Nuclear Single Spectral Peak

MA Ying-jie, ZHOU Jing*, HONG Xu, ZHOU Jian-bin, WANG Min, WAN Wen-jie

The College of Nuclear Technology and Automation Engineering, Chengdu University of Technology, Chengdu 610059, China

In nuclear spectrum, Gaussian function least square fitting is a commonly used method. Usually the method has high precision, but it is very much sensitive to noise, which causes that the residual vector is larger near the peak in the Gaussian function. To solve the problem, Gaussian function least square fitting was deduced particularly, and the causes are analyzed. As a result, Gaussian function weighted least square fitting is proposed, i.e., a weight factor, which had an opposite tendency to the data weight reduction tendency after taking logarithm, or it had the same tendency to the origin data. This was introduced based on Gaussian function least square fitting to reduce noise sensitivity. In the process of solving Gaussian parameter, to improve the real-time performance, the solution process of inverse matrix was transferred to the solution process of simple equations because the computation of inverse matrix was time consuming. Gaussian function parameter, amplitude, center value and variance, were given with the fast calculation formulas. By applying these two methods to the practical fitting of55Fe characteristic X-ray single spectrum peak, respectively, the results show that Gaussian function weighted least square fitting is more satisfactory. It indicates the proposed method can decrease the noise sensitivity and reduce the residual vector near the peak; in addition, the fitting precision is also improved. What's more, the real-time performance is improved by applying fast calculation formulas, which makes it possible to apply the proposed method to portable equipment efficiently.

Nuclear spectrum; Single spectrum peak; Weight factor; Weighted least square fitting

Dec. 20, 2015; accepted Apr. 15, 2016)

2015-12-20，

2016-04-15

国家自然科学基金项目(11475036， 41404108)资助

马英杰， 女， 1970年生， 成都理工大学核技术与自动化工程学院副教授 e-mail： ma.yingjie@qq.com *通讯联系人 e-mail： zjsin@sina.cn

TL84

A

10.3964/j.issn.1000-0593(2016)08-2373-05