馬英杰， 周 靖， 洪 旭， 周建斌， 王 敏， 萬文杰

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

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

馬英杰， 周 靖*， 洪 旭， 周建斌， 王 敏， 萬文杰

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

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

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

引 言

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

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

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

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

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

(1)

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

(2)

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

(3)

將式(3)展開， 得

(4)

令

則式(4)可表示為

(5)

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

(6)

下面對系數(shù)c1，c2和c3的求解過程進行推導。 式(5)的實質(zhì)是一個超定方程組

(7)

將其簡寫為

Ac=b

其中，

(8)

由于是超定方程組， 所以采用最小二乘法來求解。 設(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)可等價表示為

(15)

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

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

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

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

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

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

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

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

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

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

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

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

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

(16)

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

s=W(b-Ac)=Wr

(17)

那么，

bTWTWb-bTWTWAc-cTATWTWb+cTATWAc

(18)

將ρ對向量c求導得

(19)

ATWTWAc=ATWTWb

(20)

同樣， 為了避免矩陣求逆運算， 令

ATWTWA，Q=ATWTWb

(21)

則式(20)可表示為

Pc=Q

(22)

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

(23)

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

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

3 實驗分析

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

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

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

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

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

(24)

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

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

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

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

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

4 結(jié)束語

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

[1] PANG Ju-feng(龐巨豐). γ Spectrum Data Analysis(γ能譜數(shù)據(jù)分析). Xi’an： Shaanxi Science and Technology Press(西安： 陜西科學技術(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(復旦大學, 清華大學, 北京大學). 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(同濟大學計算數(shù)學教研室). Modern Numerical Mathematics and Computing(現(xiàn)代數(shù)值數(shù)學和計算). Shanghai: Tongji University Press(上海: 同濟大學出版社), 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(長沙: 國防科技大學出版社), 2003. 1.

[8] YANG Zhi-hui, LI Yue-zhong, YU Qian, et al(楊志輝, 李躍忠, 余 倩， 等). Nuclear Electronics & Detection Technology(核電子學與探測技術(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(核與粒子物理實驗方法). 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年生， 成都理工大學核技術(shù)與自動化工程學院副教授 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