核能譜單能峰快速高精度高斯函數(shù)擬合算法研究

馬英杰， 周 靖， 洪 旭， 周建斌， 王 敏， 萬(wàn)文杰

成都理工大學(xué)核技術(shù)與自動(dòng)化工程學(xué)院， 四川 成都 610059

核能譜單能峰快速高精度高斯函數(shù)擬合算法研究

馬英杰， 周 靖*， 洪 旭， 周建斌， 王 敏， 萬(wàn)文杰

成都理工大學(xué)核技術(shù)與自動(dòng)化工程學(xué)院， 四川 成都 610059

在核能譜分析中， 高斯函數(shù)最小二乘擬合法是計(jì)算單能峰凈峰面積常用的方法， 該方法精度較高， 但噪聲敏感性較強(qiáng)， 導(dǎo)致擬合出的高斯函數(shù)在峰位附近的殘差向量較大。 針對(duì)該問(wèn)題， 對(duì)高斯函數(shù)最小二乘擬合法進(jìn)行了詳細(xì)推導(dǎo)， 分析了峰位附近殘差向量較大的原因， 提出了一種基于高斯函數(shù)最小二乘擬合法的高斯函數(shù)加權(quán)最小二乘擬合法， 即在高斯函數(shù)最小二乘擬合法的基礎(chǔ)上， 引入了權(quán)重因子。 該權(quán)重因子與取對(duì)數(shù)后數(shù)據(jù)權(quán)重削弱趨勢(shì)相反或與數(shù)據(jù)本身趨勢(shì)相符， 以減小噪聲敏感性。 由于在求解高斯函數(shù)參數(shù)的過(guò)程中涉及到求逆矩陣運(yùn)算， 計(jì)算量較大， 耗時(shí)較長(zhǎng)， 為了提高實(shí)時(shí)性， 將求逆矩陣的運(yùn)算過(guò)程轉(zhuǎn)換為了簡(jiǎn)單的方程組運(yùn)算， 并給出了高斯函數(shù)的幅值、 中心及方差參數(shù)的快速求解公式。 將這兩種方法用于55Fe的特征X射線單能峰的實(shí)際擬合中， 結(jié)果表明， 高斯函數(shù)加權(quán)最小二乘擬合法效果均較好， 這說(shuō)明該方法降低了噪聲敏感性， 減小了高斯函數(shù)在峰位附近的殘差向量， 進(jìn)一步提高了擬合精度。 另外， 使用快速求解公式， 也減小了運(yùn)算量， 增強(qiáng)了實(shí)時(shí)性， 為在便攜式設(shè)備中的有效使用提供了可能。

核能譜； 單能峰； 權(quán)重因子； 加權(quán)最小二乘擬合

引 言

在核能譜分析中， 確定凈峰面積的方法基本上可分為兩類： (1)計(jì)數(shù)相加法， 即把峰內(nèi)的各道計(jì)數(shù)按照一定的公式直接相加。 這種方法比較簡(jiǎn)單， 只適用于確定單能峰面積； (2)高斯函數(shù)最小二乘擬合法， 即將峰內(nèi)各道計(jì)數(shù)擬合為一個(gè)高斯函數(shù)， 然后對(duì)這個(gè)函數(shù)進(jìn)行積分， 從而得到凈峰面積， 這種方法較計(jì)數(shù)相加法， 結(jié)果更準(zhǔn)確， 也適用于重疊峰[1-4]。

在進(jìn)行高斯函數(shù)最小二乘擬合時(shí)， 高斯函數(shù)的幅值、 中心及方差參數(shù)求解的準(zhǔn)確性決定了函數(shù)擬合的好壞， 是影響解譜效果的關(guān)鍵因素。 高斯函數(shù)參數(shù)求解的方法有： (1)文獻(xiàn)[5-9]使用的LM算法(levenberg-marquardt算法)具有速度快， 不易發(fā)散的特點(diǎn)， 但涉及到函數(shù)初始值的選取； (2)采用矩陣， 對(duì)參數(shù)進(jìn)行求解， 但需要進(jìn)行矩陣求逆運(yùn)算， 計(jì)算量較大。 此外， 采用高斯函數(shù)最小二乘擬合法得到的幅值、 中心及方差對(duì)噪聲有很強(qiáng)的敏感性[10]。 文獻(xiàn)[10]通過(guò)高斯函數(shù)的第二特征函數(shù)得到高斯函數(shù)的幅值、 中心及方差， 該方法具有求解速度快， 而且對(duì)噪聲不敏感的特點(diǎn)。

對(duì)單能峰的高斯函數(shù)最小二乘擬合法進(jìn)行研究， 旨在求解高斯擬合參數(shù)、 提高擬合精度、 降低噪聲敏感性、 增強(qiáng)擬合實(shí)時(shí)性等方面進(jìn)行改進(jìn)與提高。

1 高斯函數(shù)擬合原理

射線與物質(zhì)作用產(chǎn)生電離和激發(fā)的統(tǒng)計(jì)漲落以及譜儀固有的漲落特性， 使得探測(cè)器的輸出脈沖幅度圍繞某均值小幅漲落[11-12]。 若能較好的扣除本底， 則能譜單能峰的凈計(jì)數(shù)近似服從高斯分布， 即在譜峰區(qū)內(nèi)各道計(jì)數(shù)y(x)與道址x的關(guān)系為

(1)

式中， 變量x代表道址或能量，μ是峰高對(duì)應(yīng)的道址或能量，σ是表征峰形寬度的特征量， 它與半高寬度(FWHM)的關(guān)系為： FWHM=2.354 28σ。 若參數(shù)y0和σ已知， 則譜峰面積A可以通過(guò)積分計(jì)算得到

(2)

高斯函數(shù)最小二乘擬合的目的就是求出參數(shù)y0，μ和σ， 下面對(duì)求解過(guò)程進(jìn)行詳細(xì)推導(dǎo)。 對(duì)式(1)兩邊同時(shí)取對(duì)數(shù)得

(3)

將式(3)展開， 得

(4)

令

則式(4)可表示為

(5)

通過(guò)上面的變換， 將一個(gè)超越函數(shù)變?yōu)榱撕?jiǎn)單的二次函數(shù)， 從而將高斯函數(shù)擬合轉(zhuǎn)換為了二次多項(xiàng)式函數(shù)的最小二乘擬合。 若能解出系數(shù)c1，c2和c3， 便可根據(jù)式(6)求出高斯函數(shù)的參數(shù)。

(6)

下面對(duì)系數(shù)c1，c2和c3的求解過(guò)程進(jìn)行推導(dǎo)。 式(5)的實(shí)質(zhì)是一個(gè)超定方程組

(7)

將其簡(jiǎn)寫為

Ac=b

其中，

(8)

由于是超定方程組， 所以采用最小二乘法來(lái)求解。 設(shè)殘差向量表示為

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)

根據(jù)式(8)， 式(13)可等價(jià)表示為

(15)

根據(jù)式(14)和式(15)， 求解三元一次方程組， 即可求得系數(shù)c1，c2和c3， 再由式(6)解得σ2，μ，y0， 便可快速得到高斯函數(shù)的參數(shù)。 這樣的運(yùn)算過(guò)程避免了矩陣求逆的大量運(yùn)算， 簡(jiǎn)化了運(yùn)算過(guò)程， 提高了擬合速度， 算法也更為簡(jiǎn)潔， 更有利于在便攜式設(shè)備中實(shí)現(xiàn)。

采用上述步驟對(duì)高斯函數(shù)Y：y=400exp(-(x-50)2/(2×102))進(jìn)行最小二乘擬合， 結(jié)果如圖1所示(為了便于比較， 后面模擬使用的高斯函數(shù)均為Y， 數(shù)據(jù)用折線畫出)。

圖1 最小二乘擬合結(jié)果

將最小二乘擬合結(jié)果作線性分析， 結(jié)果如圖2所示。

圖2 最小二乘擬合結(jié)果線性分析

由圖2得， 最小二乘擬合結(jié)果與原高斯函數(shù)值之間的線性相關(guān)系數(shù)r為1， 說(shuō)明按照上面的方法進(jìn)行高斯函數(shù)的最小二乘擬合的結(jié)果是可信的。

實(shí)際的能譜中包含有噪聲， 這里以高斯函數(shù)Y與隨機(jī)噪聲疊加后的結(jié)果近似表示一個(gè)能譜單能峰， 按照上面的擬合步驟， 擬合結(jié)果如圖3中曲線b1、 圖4中曲線b1所示。

本研究使用構(gòu)式搭配分析軟件Coll, analysis 3.2a在R語(yǔ)言環(huán)境下進(jìn)行運(yùn)算，之后進(jìn)行Fisher精確檢驗(yàn)(Fisher exact test)，統(tǒng)計(jì)出槽位中的動(dòng)詞與目標(biāo)構(gòu)式的關(guān)聯(lián)強(qiáng)度。再使用WordNet 2.1進(jìn)行語(yǔ)義分析。R語(yǔ)言是成熟的編程語(yǔ)言，具有許多實(shí)用的程序包，如amap，cluster，fpc等，既能夠靈活地進(jìn)行數(shù)據(jù)檢索分析，又能提供個(gè)性化的數(shù)據(jù)統(tǒng)計(jì)，此外還具有繪圖功能[37]。

圖3 單能峰最小二乘擬合結(jié)果一

圖4 單能峰最小二乘擬合結(jié)果二

綜上， 最小二乘擬合結(jié)果對(duì)噪聲有很強(qiáng)的敏感性， 而實(shí)際的能譜中往往包含有噪聲， 因此減小噪聲對(duì)高斯擬合的影響， 具有實(shí)際意義。

2 高斯函數(shù)最小二乘擬合的改進(jìn)

為解決上述問(wèn)題， 在擬合公式中為殘差向量r添加一個(gè)權(quán)重因子W。 這里的權(quán)重因子應(yīng)該與取對(duì)數(shù)運(yùn)算后數(shù)據(jù)權(quán)重削弱趨勢(shì)相反或者與數(shù)據(jù)本身趨勢(shì)相符合。 這里， 取數(shù)據(jù)本身(取對(duì)數(shù)前的值)的歸一化值為權(quán)重因子， 即

(16)

結(jié)合式(9)， 設(shè)

s=W(b-Ac)=Wr

(17)

那么，

bTWTWb-bTWTWAc-cTATWTWb+cTATWAc

(18)

將ρ對(duì)向量c求導(dǎo)得

(19)

ATWTWAc=ATWTWb

(20)

同樣， 為了避免矩陣求逆運(yùn)算， 令

ATWTWA，Q=ATWTWb

(21)

則式(20)可表示為

Pc=Q

(22)

結(jié)合式(8)， 式(16)， 式(21)可等價(jià)表示為

(23)

根據(jù)式(22)和式(23)求解三元一次方程組， 即可求得系數(shù)c1， c2和c3， 再由式(6)解得σ2， μ， y0， 便可快速得到加權(quán)的高斯擬合函數(shù)。 對(duì)圖3、 圖4中的單能峰分別進(jìn)行加權(quán)最小二乘擬合， 擬合結(jié)果如圖3中曲線c1、 圖4中曲線c2所示。

對(duì)比圖3和圖4中擬合結(jié)果， 可以看出， 采用加權(quán)最小二乘擬合后， 在峰位附近的殘差向量減小， 擬合效果得到明顯改善， 說(shuō)明引入權(quán)重因子后， 降低了噪聲敏感性。

3 實(shí)驗(yàn)分析

實(shí)測(cè)55Fe標(biāo)樣的X射線能譜， 采用5點(diǎn)平滑進(jìn)行譜光滑， 使用SNIP方法扣除本底后， 分別采用最小二乘擬合法以及加權(quán)最小二乘擬合法的高斯函數(shù)擬合結(jié)果， 分別如圖5中曲線b、 曲線c所示。

圖5 55Fe全能峰擬合結(jié)果

將圖5中道址826—836的擬合結(jié)果進(jìn)行局部放大， 結(jié)果如圖6所示。

圖6 擬合結(jié)果局部放大

由圖6可知， 采用加權(quán)最小二乘擬合后， 在峰位附近的殘差向量較采用最小二乘擬合的殘差向量小。 采用優(yōu)良指數(shù)AIFOM(analytic improved figure of merit)對(duì)擬合效果的優(yōu)劣(擬合優(yōu)度)按照式(24)進(jìn)行評(píng)價(jià)：

(24)

式中np為擬合峰所包含的道址；n為本底和峰的道數(shù)之和；Ap為全能峰凈計(jì)數(shù)； Δy為第i道的實(shí)測(cè)值與擬合值之差。 判斷擬合效果好壞的標(biāo)準(zhǔn)是： AIFOM<0.01%， 擬合效果較好； 0.010.05%， 擬合效果差。

對(duì)55Fe標(biāo)樣的X射線能譜進(jìn)行多次實(shí)測(cè)， 對(duì)全能峰進(jìn)行擬合， 擬合優(yōu)度以及凈峰面積如表1所示(圖5中的55Fe實(shí)測(cè)譜的擬合優(yōu)度及凈峰面積見實(shí)測(cè)譜線a)。

表1 擬合優(yōu)度、 凈峰面積對(duì)比表

從表1可以看出， 5次測(cè)量中， 采用最小二乘擬合時(shí)， 擬合優(yōu)度有四次都大于0.01%， 擬合效果較差； 而采用加權(quán)最小二乘擬合后， 擬合優(yōu)度均小于0.01%， 擬合效果較好。 對(duì)同一個(gè)全能峰來(lái)說(shuō)， 采用加權(quán)最小二乘擬合后， 擬合優(yōu)度有所降低， 擬合效果明顯改善， 說(shuō)明擬合值越接近實(shí)際值。 比較加權(quán)前后擬合的凈峰面積， 變化率最高達(dá)2.76%， 最小為1.01%， 這說(shuō)明凈峰面積前后的變化還是比較大的。 在高放射性場(chǎng)合中， 凈峰面積的變化可能會(huì)更大。

綜上， 與最小二乘擬合法相比， 采用加權(quán)最小二乘擬合后， 降低了噪聲敏感性， 提高了高斯函數(shù)的擬合精度， 從而凈峰面積也更為準(zhǔn)確。

4 結(jié)束語(yǔ)

對(duì)單能峰的高斯函數(shù)最小二乘擬合的過(guò)程進(jìn)行了詳細(xì)推導(dǎo)， 在求解高斯函數(shù)的幅值、 中心及方差參數(shù)時(shí)， 將矩陣求逆運(yùn)算轉(zhuǎn)換為了簡(jiǎn)單的方程組運(yùn)算， 并給出了參數(shù)快速求解的具體公式， 大大簡(jiǎn)化了運(yùn)算過(guò)程， 提高了擬合速度。 分析了噪聲對(duì)高斯函數(shù)最小二乘擬合法的影響， 并引入了權(quán)重因子對(duì)該算法加以改進(jìn)。 實(shí)例證明， 改進(jìn)后的算法(加權(quán)最小二乘擬合法)減小了單能峰峰位附近的殘差向量， 降低了噪聲敏感性， 提高了高斯函數(shù)擬合精度。 結(jié)合快速求解公式， 加權(quán)最小二乘擬合法是一種擬合速度快且性能較優(yōu)良的高斯峰面積計(jì)算方法。

[1] PANG Ju-feng(龐巨豐). γ Spectrum Data Analysis(γ能譜數(shù)據(jù)分析). Xi’an： Shaanxi Science and Technology Press(西安： 陜西科學(xué)技術(shù)出版社), 1990, 681.

[2] QI Rong, MAO Yong, CHEN Xi-meng(齊 榮, 毛 永, 陳熙萌). Nuclear Techniques(核技術(shù)), 2008, 31(5): 330.

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

[4] Fudan University, Tsinghua University, Peking University(復(fù)旦大學(xué), 清華大學(xué), 北京大學(xué)). The Experimental Method of Nuclear Physics, Volume 1(原子核物理實(shí)驗(yàn)方法, 上冊(cè)). 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(同濟(jì)大學(xué)計(jì)算數(shù)學(xué)教研室). Modern Numerical Mathematics and Computing(現(xiàn)代數(shù)值數(shù)學(xué)和計(jì)算). Shanghai: Tongji University Press(上海: 同濟(jì)大學(xué)出版社), 2004. 56.

[7] XIE Zheng, LI Jian-hua, TANG Ze-ying(謝 政, 李建華, 湯澤瀠). Nonlinear Optimization(非線性最優(yōu)化). 2nd ed(第2版). Changsha: National University of Defence Technology Press(長(zhǎng)沙: 國(guó)防科技大學(xué)出版社), 2003. 1.

[8] YANG Zhi-hui, LI Yue-zhong, YU Qian, et al(楊志輝, 李躍忠, 余 倩， 等). Nuclear Electronics & Detection Technology(核電子學(xué)與探測(cè)技術(shù)), 2013, 33(10): 1271.

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

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

[11] WANG Shao-shun(王韶舜). Nuclear and Particle Physics Experimental Methods(核與粒子物理實(shí)驗(yàn)方法). Beijing: Atomic Press(北京: 原子能出版社), 1989. 208.

[12] WU Zhi-hua(吳治華). Nuclear Physics Experimental Methods(原子核物理實(shí)驗(yàn)方法). 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

國(guó)家自然科學(xué)基金項(xiàng)目(11475036， 41404108)資助

馬英杰， 女， 1970年生， 成都理工大學(xué)核技術(shù)與自動(dòng)化工程學(xué)院副教授 e-mail： ma.yingjie@qq.com *通訊聯(lián)系人 e-mail： zjsin@sina.cn

TL84

A

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