楊萬里,劉 沖,胡 銳

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

0 引 言

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

1 算法理論分析

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

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

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

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

yi—道址為i時的能譜；

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

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

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

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

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

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

(7)

(8)

|cn|<2Scn

(9)

(10)

2 仿真驗證

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

為了證明階數(shù)自適應(yīng)型正交多項式模型法對本底處理的效果，用式(11)～式(13)中的p1～p3和多項式b4構(gòu)造合成譜區(qū)，其中p1～p3是高斯函數(shù)，用于模擬譜線特征峰；b4是4階多項式，用于模擬復(fù)雜本底。合成譜區(qū)見圖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—多項式。

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

(19)

3 結(jié)果與討論

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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