汪雪元，何劍鋒*，劉 琳，聶逢君

1. 東華理工大學(xué)江西省放射性地學(xué)大數(shù)據(jù)技術(shù)工程實驗室，江西 南昌 330013 2. 東華理工大學(xué)江西省核地學(xué)數(shù)據(jù)科學(xué)與系統(tǒng)工程技術(shù)研究中心，江西 南昌 330013 3. 東華理工大學(xué)軟件學(xué)院，江西 南昌 330013

引 言

X射線熒光光譜(XRF)分析中，找到光譜的所有特征峰和準確地確定光譜中各個峰的峰位是XRF定性分析的關(guān)鍵問題。 常用的尋峰方法有簡單比較法、導(dǎo)數(shù)法[1-2]、協(xié)方差法、高斯乘積函數(shù)法、對稱零面積變換法、高斯擬合法、線性擬合法和傅里葉自去卷積法等，以及新引入的連續(xù)小波變換法[3-6]、Boosted-Gold反卷積法[7-8]、模擬退火法、蟻群算法、遺傳算法和人工神經(jīng)網(wǎng)絡(luò)法等新的尋峰方法[9]。 這些算法各有優(yōu)缺點，光譜分析過程中尋峰方法的選擇需要根據(jù)光譜的具體情況和具體分析要求而定。

導(dǎo)數(shù)法尋峰原理簡單，算法速度快。 一階導(dǎo)數(shù)對孤立的弱峰靈敏度高，但無法分辨重峰[10]，且導(dǎo)數(shù)法對XRF信號的信噪比有較高要求。 連續(xù)小波變換法尋峰，基本不受信號中噪聲的影響，并且有較強的重峰分辨能力[3, 5]。 但和一階導(dǎo)數(shù)法相比，連續(xù)小波變換法的弱峰識別能力不夠。 針對導(dǎo)數(shù)法和連續(xù)小波變換法尋峰的優(yōu)缺點，本文提出了一種綜合利用連續(xù)小波變換和一階導(dǎo)數(shù)聯(lián)合進行尋峰的自適應(yīng)算法，該算法具備一階導(dǎo)數(shù)法的弱峰識別能力和連續(xù)小波變換法的重峰分辨能力。

1 導(dǎo)數(shù)法尋峰

1.1 導(dǎo)數(shù)法簡介

導(dǎo)數(shù)法是一種常用的光譜尋峰方法。 將譜線看成一條連續(xù)的曲線，利用導(dǎo)數(shù)求得曲線的極值點、拐點，進而可以確定光譜特征峰的峰址、峰寬。 XRF中元素的特征峰峰型是一個類似高斯分布的函數(shù)。 極大值點的一階導(dǎo)數(shù)為向下的過零點。 一階導(dǎo)數(shù)法對孤立的峰具有較高的靈敏度，并且峰址定位精確[10]。 但其不能分辨重峰，且對噪聲敏感。

1.2 導(dǎo)數(shù)法尋峰分析

表1為某石油錄井巖性分析標(biāo)準樣品1中元素含量，為各元素的質(zhì)量分數(shù)。 標(biāo)準樣品1的EDXRF光譜數(shù)據(jù)由Si-PIN探測器測得，實測譜線如圖1所示。

表1 標(biāo)準樣品1中元素含量Table 1 The content of elements in standard sample 1

圖1 標(biāo)準樣品1實測譜線Fig.1 Measured spectrum of standard sample 1

對標(biāo)準樣品1實測譜線在未平滑的情況下用一階導(dǎo)數(shù)法尋峰，識別出139個峰，如圖2(a)所示。 對標(biāo)準樣品1實測譜線用均值平滑算法進行平滑并扣除本底后尋峰，當(dāng)均值平滑窗口大小為5時，識別出43個峰，如圖2(b)所示。 但從圖2(b)可以觀察到，峰9和峰10包含的重疊峰并沒有識別出來。 可知： 一階導(dǎo)數(shù)法具有較強的弱峰識別能力，但重疊峰分辨能力較差。

當(dāng)均值平滑窗口大小為3，7，9，11，13，15，17和19時，識別出的峰個數(shù)分別為70，28，20，19，18，17，15和15。 可知，當(dāng)均值平滑窗口逐漸增大時，識別出的峰逐漸減少。 且可看出，當(dāng)均值平滑窗口較小時，平滑窗口大小變動對譜峰識別能力影響較大； 當(dāng)均值平滑窗口較大時，平滑窗口大小變動對譜峰識別能力影響接近最大。

圖2 一階導(dǎo)數(shù)法尋峰效果(a)： 未平滑的尋峰結(jié)果； (b)：均值平滑后尋峰結(jié)果-窗口大小為5Fig.2 Peak-finding results of first-order derivative method(a)： Peak-finding results without smoothing； (b)： Peak-finding results after mean smoothing-window size is 5

2 小波變換法尋峰

2.1 連續(xù)小波變換法簡介

由小波變換的性質(zhì)可以將f的小波變換Wf(u,s)寫成式(1)的形式。

(1)

式(1)中，Ψ為小波函數(shù)；θ為一使小波函數(shù)Ψ具有n階消失矩的快速衰減函數(shù)。 從式(1)可以看出，信號光滑后的n階導(dǎo)數(shù)可以通過小波變換來實現(xiàn)。 可見，小波變換具有自動去噪功能。 由于XRF的特征峰近似于高斯函數(shù)，高斯系列小波gausn常用于XRF的尋峰。 其中，gaus2，gaus4，gaus6和gaus8的作用類似于導(dǎo)數(shù)法尋峰中的偶數(shù)階導(dǎo)數(shù)，常用于重峰分解，函數(shù)如圖3。

由圖3可知，隨階數(shù)的增加，高斯小波的形狀變得更窄，即重峰分解能力更強。 但小波函數(shù)兩側(cè)的波瓣也變大，導(dǎo)致變換后的光譜中會出現(xiàn)干擾峰，不利于光譜中的弱峰識別。 小波函數(shù)gaus2兩側(cè)沒有值大于0的波瓣，選用gaus2對光譜進行處理，不會增加直接正向的干擾峰。

2.2 連續(xù)小波變換法尋峰過程

連續(xù)小波變換尋峰過程： (1)對XRF進行連續(xù)小波變換，得到不同分解尺度的小波變換系數(shù)； (2)選擇合適分解尺度的小波變換系數(shù)[11]，例如可以選擇峰幅值最接近原始譜線峰幅值的小波變換系數(shù)譜； (3)從選定小波變換系數(shù)譜中尋峰。 用gaus2對標(biāo)準樣品1實測譜線進行連續(xù)小波變換，得到不同分解尺度的小波變換系數(shù)。 圖4顯示了分解尺度1～8的小波變換系數(shù)。

圖3 偶數(shù)階高斯小波函數(shù)Fig.3 Even-order Gauss wavelet functions

圖4 分解尺度1～8時的小波變換系數(shù)Fig.4 Wavelet transform coefficients when the decomposition scale changes from 1 to 8

通過對不同分解尺度小波變換系數(shù)的分析，可知： (1)分解尺度較小時，小波變換系數(shù)含有大量噪聲，但對重峰分解能力最強； (2)隨著分解尺度的增加，系數(shù)逐漸變得光滑，分辨率逐漸變大，弱峰識別能力增強，但重峰識別能力變?nèi)酰?(3)當(dāng)分解尺度變得更大時，部分弱峰被平滑掉，逐漸變得不可識別。 從峰高可以看出圖4中標(biāo)準樣品1實測譜線的小波變換系數(shù)d6最接近原始譜線。 對小波變換系數(shù)d6進行尋峰，結(jié)果如圖5所示。 尋到的峰個數(shù)為33個。 對這些峰進行識別，識別出20種元素： Mg，Al，Si，P，S，Cl，K，Sn，Ca，Ti，V，Cr，F(xiàn)e，Co，Ni，Cu，W，Zn，Ga和Th。 而其余9種元素并未被識別出來。

圖5 連續(xù)小波變換法尋峰Fig.5 Peak-finding results by continuous wavelet transform method

從上述分析可知，選擇最優(yōu)分解尺度小波變換系數(shù)的常規(guī)小波變換尋峰法無法同時滿足較強的重疊峰分辨能力和弱峰識別能力。

3 連續(xù)小波變換與導(dǎo)數(shù)法聯(lián)合尋峰

3.1 連續(xù)小波變換與導(dǎo)數(shù)法聯(lián)合尋峰過程

本文提出的連續(xù)小波變換與導(dǎo)數(shù)法聯(lián)合尋峰的自適應(yīng)算法步驟如下： (1)對原始光譜用小波變換迭代法扣除本底[6]； (2)對扣除本底后光譜用高斯小波gaus2進行變換，得到不同分解尺度的小波變換系數(shù)CFi(1≤i≤Scale，Scale為一最大分解尺度，例如32)； (3)對不同分解尺度的小波變換系數(shù)譜用一階導(dǎo)數(shù)法尋峰，得到峰集合PKi(1≤i≤Scale)； (4)合并集合PKi(1≤i≤Scale)，得到UPK。 在合并過程中，以PK1為UPK初始值，然后依次將PK2，PK3，…，PKScale中的首次出現(xiàn)的峰加入UPK。 (5)由能量E與道址N之間的關(guān)系表達式計算出X熒光分析儀測量元素范圍內(nèi)輕元素的K系和重元素的L系特征X射線對應(yīng)的道址； (6)根據(jù)集合UPK和元素特征X射線對應(yīng)的道址，確定原始光譜所含元素并確定道址。

3.2 連續(xù)小波變換與導(dǎo)數(shù)法聯(lián)合尋峰實例分析

按照上述算法步驟對標(biāo)準樣品1的實測譜線扣除本底，連續(xù)小波變換后得到不同分解尺度的小波變換系數(shù)CFi(1≤i≤Scale，Scale=32)。 然后對小波變換系數(shù)譜尋峰，步驟如下：

(1)小波變換系數(shù)譜尋峰。 對小波變換系數(shù)譜CFi尋峰，得到不同分解尺度峰集合PKi(1≤i≤Scale)，尋峰結(jié)果如表2所示。 分解尺度為1的小波變換系數(shù)譜尋峰結(jié)果如圖6(a)所示，識別的峰個數(shù)為22。 從圖6(a)可以看出，原始光譜中能夠直觀辨別的重疊峰均已被識別出。

表2 不同分解尺度的小波變換系數(shù)譜的尋峰結(jié)果

(2)合并峰集合。 合并不同分解尺度峰集合PKi(1≤i≤Scale)，得到的UPK含有107個峰，如圖6(b)所示。 有些元素的特征X射線在UPK中有多條，它們是不同分解尺度下尋峰的結(jié)果，在后續(xù)步驟中會確定哪條該保留。

(3)確定X熒光分析儀的能量刻度。 由不同的標(biāo)準樣品得到能量E和道址N的關(guān)系式如式(2)所示。 計算出輕元素(原子序數(shù)11≤Z≤45的元素)的K系和重元素(原子序數(shù)45

E=0.007 599N+0.134 714 (keV)

(2)

(4)確定原始光譜所含元素。 根據(jù)集合UPK和元素特征X射線標(biāo)準刻度，進行元素識別，得到結(jié)果如表3所示。 表中Scale值表示對應(yīng)元素是在相應(yīng)小波分解尺度下被識別的。

圖6 自適應(yīng)小波變換導(dǎo)數(shù)法尋峰(a)： 分解尺度為1的尋峰結(jié)果； (b)： 合并尋峰結(jié)果Fig.6 Peak-finding by adaptive method based on derivative and wavelet transform(a)： Peak-finding results for scale 1； (b)： Union of peak-finding results

表3 尋峰結(jié)果Table 3 Results of peak-finding

(5)尋峰結(jié)果分析。 由表1和表3比較后可知： 除元素Na，Zr和Nb外的26種元素都被確定。 X熒光分析儀在低能端的探測效率較低，原子序數(shù)靠前的Na元素?zé)o法被探測到。 當(dāng)?shù)乐種為2 048時，由式(2)可以計算出該道址對應(yīng)的能量為15.697 5(keV)，而Zr和Nb的特征X射線的Ka1值分別為15.774(keV)和16.614(keV)，超出了X熒光分析儀元素分析范圍。 可以看出，含量比例高的元素在小波分解尺度為1和2時即可識別。 而含量比例低的元素，有些在小波分解尺度取較大值時才能識別，例如Ga，As和Pb元素。 因為一階導(dǎo)數(shù)法對孤立的峰具有較高的靈敏度，含量比例極低的元素W在小波分解尺度為3時就被識別。

4 結(jié) 論

分析比較了導(dǎo)數(shù)法和連續(xù)小波變換法的XRF尋峰效果，總結(jié)了導(dǎo)數(shù)法和連續(xù)小波變換法在對XRF進行尋峰時的優(yōu)缺點，并提出了一個綜合利用高斯小波gaus2和一階導(dǎo)數(shù)對XRF進行尋峰的自適應(yīng)算法。 該算法對不同分解尺度的小波變換系數(shù)譜進行一階導(dǎo)數(shù)尋峰，然后合并尋峰結(jié)果，最后進行元素識別。 本文提出的自適應(yīng)尋峰算法具備一階導(dǎo)數(shù)法的弱峰識別能力，小波變換尋峰法的重疊峰分辨能力和自動去噪能力。 新方法對所有分解尺度的小波變換系數(shù)譜進行尋峰并合并尋峰結(jié)果，優(yōu)于只對某一分解尺度小波變換系數(shù)譜進行尋峰的常規(guī)小波變換尋峰法。 按照新算法得到的特征峰集合可以作為后續(xù)精確尋峰、峰面積計算工作的初始值。