趙奉奎，徐曉美，呂立亞

(1.南京林業(yè)大學(xué) 汽車與交通工程學(xué)院，江蘇 南京 210037；2.東南大學(xué) 儀器科學(xué)與工程學(xué)院，江蘇 南京 210096)

光譜本底常使得凈峰面積的估算結(jié)果過(guò)大和峰位估算結(jié)果偏移[1]。如能量色散型X射線熒光光譜(EDXRF)分析作為一種多元素分析技術(shù)，能夠同時(shí)提供元素的定性和定量分析結(jié)果，但由于其特征峰疊加在本底之上，因此，分析時(shí)為獲取特征峰信息，必須扣除本底的影響[2]。

為了扣除本底，研究人員已經(jīng)提出了多種算法，一種方法是從硬件的角度出發(fā)，通過(guò)改進(jìn)光路，如設(shè)計(jì)三角形光路結(jié)構(gòu)抑制本底的產(chǎn)生[3]，更多的研究是從軟件的角度對(duì)已生成的光譜進(jìn)行本底扣除，常采用多項(xiàng)式擬合[4-5]、傅立葉變換[6]、削峰法[2,7]和小波變換(WT)[8-11]，以及近幾年研究較多的神經(jīng)網(wǎng)絡(luò)等[12]。如劉察等[13-14]將信號(hào)變換到小波空間，以減輕譜信號(hào)中基線的變化。大多數(shù)方法依賴于特征峰和本底在頻域分布在不同區(qū)域的特性，但當(dāng)兩者在頻域有重合時(shí)，一般方法不易將其分開。小波變換具有多分辨率分析的特點(diǎn)，能夠?qū)π盘?hào)進(jìn)行多尺度分析，非常適合區(qū)分特征峰和本底，然而直接利用小波變換進(jìn)行本底扣除容易導(dǎo)致光譜波形畸變。Galloway等[8]提出了一種迭代小波變換(IWT)方法，相比于傳統(tǒng)小波變換能夠?qū)π盘?hào)進(jìn)行更加細(xì)化的分解，可根據(jù)信號(hào)形態(tài)區(qū)分頻帶重復(fù)的信號(hào)，但未研究終值迭代條件。筆者利用信息熵作為迭代小波終值條件，改進(jìn)了迭代小波變換[15]方法，但計(jì)算稍復(fù)雜。

基于此，本文提出了一種更加簡(jiǎn)單的迭代終值準(zhǔn)則計(jì)算方法，改進(jìn)了基于迭代小波變換的本底扣除方法，旨在精確扣除本底，同時(shí)避免小波變換導(dǎo)致的波形畸變。

1 迭代小波變換理論

1.1 小波變換理論

小波變換能夠在時(shí)域和頻域同時(shí)進(jìn)行信號(hào)分析。在小波變換過(guò)程中，利用選定的小波基，即母小波函數(shù)ψ，對(duì)信號(hào)進(jìn)行分解。在實(shí)際運(yùn)用時(shí)，更多的是利用離散小波變換(DWT)。為了實(shí)現(xiàn)多分辨率分析，還需要構(gòu)造尺度函數(shù)φ。經(jīng)壓縮因子2j和平移因子2jn壓縮平移的小波函數(shù)和尺度函數(shù)如式(1)～(2)所示：

(1)

(2)

根據(jù)Mallat快速小波變換，信號(hào)f分解的逼近系數(shù)aj[n]和細(xì)節(jié)系數(shù)dj[n]可根據(jù)式(3)獲得:

aj[n]=〈f,φj,n〉

dj[n]=〈f,ψj,n〉

(3)

通過(guò)與濾波器h[n]與g[n]的卷積和下采樣，根據(jù)式(4)計(jì)算逼近系數(shù)aj+1[n]和細(xì)節(jié)系數(shù)dj+1[n]：

(4)

圖1 離散小波分解過(guò)程Fig.1 Process of discrete wavelet decomposition

1.2 迭代小波變換理論

根據(jù)收斂定義，對(duì)正整數(shù)n，存在正整數(shù)N，當(dāng)n>N時(shí)，對(duì)任意指定的ε>0，都有|fn+1-fn|max<ε，則認(rèn)為f收斂。對(duì)于光譜分析，如果連續(xù)3次迭代過(guò)程中fn和fn+1之間的最大誤差均小于指定ε，則認(rèn)為f收斂，具體計(jì)算過(guò)程將在下述步驟中詳述。

圖2 在不同分解層次進(jìn)行小波變換所得低頻逼近圖Fig.2 Approximations obtained through wavelet transform at different levels

圖3 基于改進(jìn)迭代小波變換的本底扣除流程圖Fig.3 Program flow chart of background removing method based on the improved IWT

進(jìn)行迭代小波變換之前，需先確定分解層數(shù)。隨著迭代小波變換過(guò)程的進(jìn)行，光譜中的特征峰成分會(huì)逐步消減。為盡量減少人工干預(yù)，所有小波變換均在相同的分解層次上進(jìn)行。最優(yōu)分解層次的確定通過(guò)在不同分解層次上進(jìn)行小波變換，比較其低頻逼近的方法實(shí)現(xiàn)。最優(yōu)的分解層次需使信號(hào)經(jīng)小波分解后，所有的本底成分均包含在低頻逼近中，盡管其中可能仍含有一些特征峰成分。經(jīng)過(guò)試驗(yàn)比較，選取DB4小波基進(jìn)行小波分解。信號(hào)在不同分解層次進(jìn)行小波變換所得低頻逼近如圖2所示。圖中A7、A8、A9分別表示第7層、第8層、第9層分解的低頻逼近。當(dāng)對(duì)信號(hào)進(jìn)行較低分解層次的小波變換時(shí)，如第7層，所得低頻本底振蕩較為劇烈，與實(shí)際本底相差較大。當(dāng)對(duì)光譜信號(hào)進(jìn)行較高分解層次的小波變換時(shí)，如第8層和第9層，所得本底的變化趨勢(shì)與實(shí)際本底比較一致，但第9層所得低頻逼近過(guò)于平坦，所以選擇第8層作為最優(yōu)分解層次。在后續(xù)的迭代小波變換過(guò)程中，均在第8層對(duì)信號(hào)進(jìn)行小波分解。

基于迭代小波變換的本底扣除過(guò)程如圖3所示，具體流程如下：

①fm[i]用來(lái)表示經(jīng)m-1次迭代小波變換后的光譜，其中m≥1為整數(shù)，表示迭代次數(shù)，i表示光譜通道，原始光譜標(biāo)記為f1[i]。

②在優(yōu)化分解層次上對(duì)fm[i]進(jìn)行小波分解，將得到的逼近系數(shù)am[i]作為第m次估計(jì)的本底。

③計(jì)算|am[i]-am-1[i]|max=errm，a0與f1[i]相等，表示原始光譜。將errm和用戶指定的誤差ε進(jìn)行比較，比較次數(shù)由times進(jìn)行索引，times初始值為0。當(dāng)對(duì)所有的通道均存在errm<ε時(shí)，表示相鄰兩次所估計(jì)的本底足夠一致，將times=times+1，跳至步驟④，否則，重置times=0，并轉(zhuǎn)至步驟⑤。

④如果times≥3，則連續(xù)4次小波變換所得估計(jì)本底足夠一致，認(rèn)為所得估計(jì)本底已經(jīng)收斂，最后一次小波變換所得低頻逼近am[i]作為最終估計(jì)的本底，從原光譜中減去am[i]即可實(shí)現(xiàn)本底扣除。

⑤比較所有通道值的fm[i]和am[i]，將各通道最小值替換fm[i]，即：

(5)

通過(guò)削去較大值的方法，迭代地去除低頻逼近中的特征峰成分，令m自增，轉(zhuǎn)至步驟②。

2 實(shí)驗(yàn)結(jié)果與討論

2.1 仿真光譜的本底扣除

為驗(yàn)證算法的有效性，首先利用仿真光譜對(duì)算法進(jìn)行分析。已有研究表明，EDXRF光譜中的特征峰可由高斯峰表征，本底可由解析函數(shù)表征[1]。因此，本研究采用高斯峰為特征峰，多項(xiàng)式函數(shù)為本底，構(gòu)建仿真光譜。

為驗(yàn)證算法是否對(duì)噪聲敏感，在光譜上增加了信噪比為38 dB的噪聲，仿真光譜如圖4A所示。

根據(jù)“1.2”所述小波分解層次選擇方法，選取第8層為最優(yōu)分解層。相鄰2次小波分解所得errm見表1。結(jié)果顯示，隨著迭代過(guò)程的進(jìn)行，errm明顯衰減。由于最大峰高7 000，選取ε=10可以認(rèn)為2次所得本底足夠一致。10次迭代后，連續(xù)3次errm<ε，表明所得估計(jì)本底已近似收斂。

Iteration number(m)12345678910errm6 12060018067301711864

圖5 迭代小波變換估計(jì)的光譜與傳統(tǒng)小波變換估計(jì)的光譜對(duì)比Fig.5 Comparison of the estimated spectra with IWT and non-iterative WT

迭代過(guò)程見圖4B，圖中第1個(gè)估計(jì)本底(1st Estimated bakground)表示第1次小波變換后得到的本底。由圖可知，第一次小波變換得到的本底與實(shí)際信號(hào)本底相差很大。隨著迭代小波變換的進(jìn)行，所得估計(jì)本底與仿真光譜的本底逐漸靠近。7次迭代后，所得估計(jì)本底收斂到仿真光譜本底。從原始仿真光譜中減去估計(jì)的本底即可實(shí)現(xiàn)本底扣除，去除本底后光譜見圖5(IWT spectrum)。同時(shí)也可以看到，在存在噪聲的前提下，該算法仍能夠準(zhǔn)確估算本底。

作為對(duì)比，對(duì)該仿真光譜利用多項(xiàng)式擬合法和傳統(tǒng)小波分解進(jìn)行本底估計(jì)。利用DB4小波對(duì)仿真光譜在第8層進(jìn)行分解，不進(jìn)行迭代，小波逼近系數(shù)作為估計(jì)的本底。實(shí)驗(yàn)發(fā)現(xiàn)，從原始光譜中去除小波變換估計(jì)的本底使光譜產(chǎn)生了較大畸變(圖5 WT spectrum)。相比之下，迭代小波變換能準(zhǔn)確扣除本底，且未引入波形畸變。

2.2 實(shí)驗(yàn)光譜的本底扣除

為進(jìn)一步驗(yàn)證算法的有效性，利用迭代小波變換對(duì)EDXRF光譜儀所測(cè)光譜進(jìn)行實(shí)驗(yàn)。光譜儀采用Si-Pin探測(cè)器，高壓為45 kV。通過(guò)對(duì)光譜的觀察和實(shí)驗(yàn)，選擇DB6在第8層對(duì)光譜進(jìn)行分析。光譜迭代分解過(guò)程見圖6A。由圖可知，經(jīng)一次小波分解后，所得估計(jì)本底與光譜實(shí)際本底相差較大，而經(jīng)小波分解迭代6次后，所得估計(jì)本底趨于收斂。

迭代小波變換過(guò)程中的errm列于表2中，考慮到特征峰最高值為8 000，ε=10時(shí)認(rèn)為所得估計(jì)本底足夠接近。經(jīng)7次迭代后，err7<ε，9次迭代后，可認(rèn)為迭代過(guò)程收斂。表中數(shù)據(jù)顯示的本底變化趨勢(shì)與圖6A中一致。

表2 實(shí)驗(yàn)光譜相鄰兩次估計(jì)本底誤差Table 2 Errs of two consecutive estimated backgrounds of the experimental spectrum

圖6 實(shí)驗(yàn)光譜迭代小波變換過(guò)程(A)和原始實(shí)驗(yàn)光譜與扣除本底光譜對(duì)比圖(B)
Fig.6 Process of iteration for the experimental spectrum(A) and comparison of the original spectrum and background-subtracted spectrum(B)

從原始光譜中減去第9次小波變換所得估計(jì)本底，即可去除本底。原始實(shí)驗(yàn)光譜與本底扣除光譜對(duì)比如圖6B所示，由圖可知，本底已完全扣除，去除本底后的光譜基線與水平坐標(biāo)非常接近，且波形未發(fā)生畸變，峰位未偏移，即使弱峰也保持完整，因此該算法也為低含量元素的分析奠定了基礎(chǔ)。

3 結(jié) 論

本文基于迭代小波變換的本底扣除改進(jìn)算法提出了一種迭代停止標(biāo)準(zhǔn)。該算法通過(guò)對(duì)仿真光譜和實(shí)驗(yàn)光譜的分析，準(zhǔn)確扣除了光譜本底，證明了算法的有效性，為光譜凈峰面積的計(jì)算和峰位的計(jì)算奠定了基礎(chǔ)。與傳統(tǒng)小波變換本底扣除法相比，該算法能夠?qū)π盘?hào)進(jìn)行更加精細(xì)的分解，從而分離特征峰和本底重復(fù)的部分，且不會(huì)導(dǎo)致波形畸變。與多項(xiàng)式擬合法相比，本底扣除結(jié)果更加準(zhǔn)確。本研究迭代的思想也可由傅立葉或多項(xiàng)式擬合來(lái)實(shí)現(xiàn)，為光譜分析提供了有力的工具。