用于能譜本底處理的階數(shù)自適應(yīng)型正交多項(xiàng)式模型法

楊萬(wàn)里,劉 沖,胡 銳

(1.南華大學(xué) 電氣工程學(xué)院,湖南 衡陽(yáng) 421001;2.中廣核久源(成都)科技有限公司,四川 成都 610059)

0 引 言

在能譜數(shù)據(jù)處理中一個(gè)重要的數(shù)據(jù)源處理過(guò)程是從復(fù)雜譜中扣除本底，本底扣除的準(zhǔn)確性及穩(wěn)定性直接影響凈峰面積計(jì)算的精確度。一個(gè)效果良好的本底扣除法[1]，對(duì)于低本底上的強(qiáng)峰、高本底上的弱峰、不同曲率本底上的峰等不同情況，處理后都能獲得準(zhǔn)確可靠的凈峰面積。除此之外，還要求盡可能減少需要用戶自己設(shè)置的參數(shù)，目的是方便大量的譜數(shù)據(jù)處理[2]。本底處理方法分為兩種思路[3-4]：本底估計(jì)(先扣除本底再計(jì)算凈峰面積)和本底模型(用數(shù)學(xué)函數(shù)描述本底，再計(jì)算凈峰面積時(shí)同時(shí)扣除)。本底估計(jì)包括：SNIP算法(scale normalization for image pyramids,SNIP)[5-7]、傅里葉變換法[8-9]、正交多項(xiàng)式擬合法[10]等。這類方法應(yīng)用的譜區(qū)間寬甚至全譜范圍、計(jì)算比較簡(jiǎn)單、速度快，但精度不高、參數(shù)較多且設(shè)置麻煩。本底模型包括：線性多項(xiàng)式模型[11]、正交多項(xiàng)式模型[3]等不同函數(shù)描述本底。它們根據(jù)高斯函數(shù)或改進(jìn)型高斯函數(shù)描述特征峰，采用最小二乘擬合法同時(shí)確定本底函數(shù)、特征峰函數(shù)的相關(guān)參數(shù)，即可求出凈峰面積。這類方法本底處理結(jié)果比較準(zhǔn)確，但譜區(qū)間較窄、計(jì)算復(fù)雜、速度不快、需要用戶根據(jù)譜形特點(diǎn)選擇本底函數(shù)及多項(xiàng)式函數(shù)的階數(shù)。本文介紹的階數(shù)自適應(yīng)型正交多項(xiàng)式模型法可以有效地將本底估計(jì)法和本底模型法結(jié)合起來(lái)，充分利用譜數(shù)據(jù)，建立起更加完善的數(shù)學(xué)模型。

1 算法理論分析

J.B.Olsen[10]提出用一系列正交多項(xiàng)式的線性組合擬合能譜，以迭代方式調(diào)整擬合權(quán)重(權(quán)重主要用于區(qū)分峰區(qū)和本底區(qū))，逐漸“剝除”譜峰，最終擬合曲線中只包含本底。由于采用正交多項(xiàng)式擬合，計(jì)算過(guò)程大大簡(jiǎn)化，且數(shù)值上更加穩(wěn)定，階數(shù)可以達(dá)到30，可用于本底曲率較大、復(fù)雜本底、擬合區(qū)間寬等情況。正交多項(xiàng)式的階數(shù)可由用戶設(shè)定。該方法的主要公式見(jiàn)式(1)～式(6)。

(1)

(2)

(3)

(4)

p0(i)=1p1(i)=i-α0

pj(i)=(i-αj-1)·pj-1(i)-βj-1·pj-2(i)

(5)

(6)

其中：i—道址，從1～N；

j—階數(shù)，從0～n；

bi—道址為i時(shí)的背景；

cj—合成譜階數(shù)為j時(shí)的系數(shù)；

r—自適應(yīng)參數(shù)；

pj(i)—道址為i時(shí)合成譜的階數(shù)；

yi—道址為i時(shí)的能譜；

wi—道址為i時(shí)的權(quán)重背景；

αj,βj和γj都是常數(shù)。

B.Vekemans等[3]提出用高斯函數(shù)表示能譜的特征峰，本底用正交多項(xiàng)式的線性組合函數(shù)表示，那么式(10)就可以表示能譜譜線。先用正交多項(xiàng)式擬合法估算本底，用式(1)表示，然后采用非線性最小二乘擬合法求解式(10)，調(diào)整正交多項(xiàng)式的系數(shù)Cj、計(jì)算出高斯函數(shù)的相關(guān)參數(shù)。該方法大大提高了譜線本底處理的精確度，但需要用戶設(shè)置正交多項(xiàng)式的階數(shù)。

正交多項(xiàng)式擬合法中正交多項(xiàng)式的階數(shù)一般由用戶根據(jù)本底的復(fù)雜程度、擬合區(qū)間設(shè)定。階數(shù)小了，不能夠準(zhǔn)確地描述本底；階數(shù)大了可能使擬合本底失真，也造成參數(shù)過(guò)多，顯然有一個(gè)最佳的階數(shù)。P.V.Espen[11]從統(tǒng)計(jì)學(xué)角度出發(fā)，研究發(fā)現(xiàn)如果第n+1階正交多項(xiàng)式的系數(shù)cn+1與它的標(biāo)準(zhǔn)差Scn+1相當(dāng)，說(shuō)明cn+1可能使擬合結(jié)果失真，如果滿足式(9)，那么階數(shù)n為最佳階數(shù)。

本文綜合上述方法提出階數(shù)自適應(yīng)型正交多項(xiàng)式模型法：從1階開(kāi)始，每次用正交多項(xiàng)式擬合法估算本底，計(jì)算最高階系數(shù)Cn的標(biāo)準(zhǔn)差Scn，如果不滿足式(9)，階數(shù)增加1；如此循環(huán)，直到第n+1階滿足式(9)，則最佳階數(shù)為n；把具有最佳階數(shù)式(1)代入式(10)，用L-M(Levenberg-Marquardt)算法求解，計(jì)算流程如圖1示。該方法不再需要用戶設(shè)定多項(xiàng)式階數(shù)，能自動(dòng)適應(yīng)各種復(fù)雜本底，也減少了多項(xiàng)式中不必要的參數(shù)個(gè)數(shù)，降低了式(10)的自由度，有利于L-M算法的收斂。

圖1 自適應(yīng)正交多項(xiàng)式模型方法流程圖

(7)

(8)

|cn|<2Scn

(9)

(10)

2 仿真驗(yàn)證

采用合成譜代替實(shí)測(cè)譜用于實(shí)驗(yàn)，并假設(shè)譜線已經(jīng)光滑處理，統(tǒng)計(jì)漲落較小，不以考慮。

為了證明階數(shù)自適應(yīng)型正交多項(xiàng)式模型法對(duì)本底處理的效果，用式(11)～式(13)中的p1～p3和多項(xiàng)式b4構(gòu)造合成譜區(qū)，其中p1～p3是高斯函數(shù)，用于模擬譜線特征峰；b4是4階多項(xiàng)式，用于模擬復(fù)雜本底。合成譜區(qū)見(jiàn)圖2，其中246～372道是感興趣區(qū)(regions of interest,ROI)。

(11)

(12)

(13)

b0=10.44

(14)

b1=0.02j+b0

(15)

b2=2.78·10-3j2+0.75j+b0

(16)

b3=4.3·10-5j3+2.09·10-2j2+2.52j+b0

(17)

b4=-8.57·10-7j4+4.46·10-4j3+

7.33·10-2j2+4.02j+b0

(18)

式中：

i—道址，從241到501；

j—變量，其值為i-240；

p—高斯函數(shù)；

b—多項(xiàng)式。

為了比較不同方法對(duì)各種曲率本底的處理效果，用p1～p3分別與b1、b2、b3、b4構(gòu)造出不同復(fù)雜程度的本底的合成譜區(qū)，分別用SNIP算法、傅里葉變換法、正交多項(xiàng)式擬合法、線性多項(xiàng)式本底模型法、正交多項(xiàng)式模型法和階數(shù)自適應(yīng)型正交多項(xiàng)式模型法處理這些合成譜的ROI的本底。

(19)

3 結(jié)果與討論

用階數(shù)自適應(yīng)型正交多項(xiàng)式模型法對(duì)圖2 ROI的本底進(jìn)行處理，處理效果如圖3所示。計(jì)算出正交多項(xiàng)式的最佳階數(shù)為4階，計(jì)算本底與真實(shí)本底十分接近，誤差在±0.5以內(nèi)，x2小于0.1。不同方法處理的結(jié)果見(jiàn)表1。

圖2 合成光譜和感興趣區(qū)

圖3 階數(shù)自適應(yīng)正交多項(xiàng)式模型方法處理本底的結(jié)果

從表1可知，對(duì)于不同曲率的本底被不同方法處理后的效果差別較大；本底模型法的處理結(jié)果整體優(yōu)于本底估計(jì)法。論文介紹的階數(shù)自適應(yīng)型正交多項(xiàng)式模型法不需用戶根據(jù)本底設(shè)置參數(shù)，且對(duì)不同曲率本底處理的效果相當(dāng)；其他方法或多或少需要根據(jù)本底的復(fù)雜程度調(diào)節(jié)一些參數(shù)，且在復(fù)雜本底時(shí)處理結(jié)果并不理想。

表1 采用不同方法進(jìn)行處理的不同曲率本底

SNIP算法、傅里葉變換法、正交多項(xiàng)式擬合都把本底當(dāng)做平緩變化部分(低曲率)、特征峰當(dāng)作快變化部分(高曲率)處理，并以此為切入點(diǎn)采用各自不同的方式，估計(jì)出本底。在伽馬能譜測(cè)量中，由于環(huán)境散射在低能端會(huì)形成一個(gè)變化較快的本底區(qū)域，該區(qū)域的本底曲率與環(huán)境中的散射材料密切相關(guān)[12-13]。另外，在含重疊峰的ROI，原本曲率較高的特征峰由于疊加效果，曲率降低。特別是多個(gè)特征峰相鄰，且峰面積相差較大時(shí)，面積小的特征峰相交區(qū)域曲率會(huì)嚴(yán)重降低。這樣不管是本底曲率高，還是由于重疊效應(yīng)使特征峰的曲率降低都會(huì)造成本底與特征峰的變化率有較大且不可忽略的重疊部分，上述算法很難區(qū)分它們。曲率重疊部分越大，本底估計(jì)的偏差越大，表1中曲率越高的本底，處理結(jié)果偏差越大，就是這個(gè)原因。而且這些方法需要根據(jù)不同曲率的本底調(diào)節(jié)參數(shù)，才能獲取較好的結(jié)果。

SNIP算法需要調(diào)節(jié)迭代次數(shù)，傅里葉變化法需要調(diào)節(jié)截至頻率、步進(jìn)頻率，正交多項(xiàng)式擬合法需要調(diào)節(jié)階數(shù)。這些參數(shù)都比較抽象，需要一定的應(yīng)用經(jīng)驗(yàn)，而且本底是隨著樣品成分、測(cè)量條件的變化而變化，不利于能譜的批量處理。正交多項(xiàng)式擬合法雖然在本底處理理論上進(jìn)步不大，但在數(shù)值算法上更簡(jiǎn)單、快捷、容易收斂，而且本底處理結(jié)果可以用正交多項(xiàng)式的線性組合表示，為相關(guān)的本底模型法提供了有利條件。

線性多項(xiàng)式模型法和正交多項(xiàng)式模型法不再只關(guān)注本底與特征峰變化率的不同，而需要建立更加具有物理意義的數(shù)學(xué)模型，特征峰用高斯函數(shù)表示，本底用數(shù)學(xué)函數(shù)表示。該方法充分利用了特征峰的數(shù)據(jù)，表1中計(jì)算結(jié)果也更加準(zhǔn)確，但依然需要根據(jù)本底的復(fù)雜程度設(shè)置多項(xiàng)式的階數(shù)。從表1可知，不同階數(shù)的擬合結(jié)果相差較大。而且線性多項(xiàng)式模型法階數(shù)不能大于4階[3]，大于4階后擬合結(jié)果可能出現(xiàn)振蕩，表1中線性多項(xiàng)式為5階時(shí)，結(jié)果出現(xiàn)振蕩，與真實(shí)本底相差很大，因此不能用于譜區(qū)間較寬、本底較復(fù)雜的情況。

自適應(yīng)型正交多項(xiàng)式模型法不僅充分利用特征峰的數(shù)據(jù)，也充分利用了本底數(shù)據(jù)，數(shù)學(xué)模型更加完善。在擬合本底時(shí)，根據(jù)本底的復(fù)雜程度自動(dòng)選擇合適的階數(shù)，無(wú)需人工設(shè)置。而高斯函數(shù)所需參數(shù)(峰位、標(biāo)準(zhǔn)差)，可以通過(guò)其它算法自動(dòng)求得近似值；峰面積是線性參數(shù)，初始值并不重要，都能滿足L-M算法。因此該方法不需要人工設(shè)置參數(shù)。表1中，該方法對(duì)不同曲率本底處理結(jié)果相當(dāng)，且都比較理想。

4 結(jié) 論

階數(shù)自適應(yīng)型正交多項(xiàng)式模型法建立了更優(yōu)良的數(shù)學(xué)模型，充分利用譜數(shù)據(jù)，能根據(jù)能譜本底的特點(diǎn)(復(fù)雜度、寬度)自動(dòng)選取最佳的正交多項(xiàng)式階數(shù)，不管是本底曲率低或高，經(jīng)過(guò)處理后都能獲得滿意的結(jié)果。該模型法不需要用戶設(shè)置任何參數(shù)，為快速、批量處理能譜數(shù)據(jù)提供了有力的支持，具有實(shí)用價(jià)值。對(duì)于某些能譜段的特征峰可能不太符合高斯響應(yīng)，要用修正的高斯函數(shù)描述，需要進(jìn)一步研究。