小波变换导数法X射线荧光光谱自适应寻峰研究

汪雪元，何剑锋*，刘 琳，聂逢君

1. 东华理工大学江西省放射性地学大数据技术工程实验室，江西 南昌 330013 2. 东华理工大学江西省核地学数据科学与系统工程技术研究中心，江西 南昌 330013 3. 东华理工大学软件学院，江西 南昌 330013

引 言

X射线荧光光谱(XRF)分析中，找到光谱的所有特征峰和准确地确定光谱中各个峰的峰位是XRF定性分析的关键问题。 常用的寻峰方法有简单比较法、导数法[1-2]、协方差法、高斯乘积函数法、对称零面积变换法、高斯拟合法、线性拟合法和傅里叶自去卷积法等，以及新引入的连续小波变换法[3-6]、Boosted-Gold反卷积法[7-8]、模拟退火法、蚁群算法、遗传算法和人工神经网络法等新的寻峰方法[9]。 这些算法各有优缺点，光谱分析过程中寻峰方法的选择需要根据光谱的具体情况和具体分析要求而定。

导数法寻峰原理简单，算法速度快。 一阶导数对孤立的弱峰灵敏度高，但无法分辨重峰[10]，且导数法对XRF信号的信噪比有较高要求。 连续小波变换法寻峰，基本不受信号中噪声的影响，并且有较强的重峰分辨能力[3, 5]。 但和一阶导数法相比，连续小波变换法的弱峰识别能力不够。 针对导数法和连续小波变换法寻峰的优缺点，本文提出了一种综合利用连续小波变换和一阶导数联合进行寻峰的自适应算法，该算法具备一阶导数法的弱峰识别能力和连续小波变换法的重峰分辨能力。

1 导数法寻峰

1.1 导数法简介

导数法是一种常用的光谱寻峰方法。 将谱线看成一条连续的曲线，利用导数求得曲线的极值点、拐点，进而可以确定光谱特征峰的峰址、峰宽。 XRF中元素的特征峰峰型是一个类似高斯分布的函数。 极大值点的一阶导数为向下的过零点。 一阶导数法对孤立的峰具有较高的灵敏度，并且峰址定位精确[10]。 但其不能分辨重峰，且对噪声敏感。

1.2 导数法寻峰分析

表1为某石油录井岩性分析标准样品1中元素含量，为各元素的质量分数。 标准样品1的EDXRF光谱数据由Si-PIN探测器测得，实测谱线如图1所示。

表1 标准样品1中元素含量Table 1 The content of elements in standard sample 1

图1 标准样品1实测谱线Fig.1 Measured spectrum of standard sample 1

对标准样品1实测谱线在未平滑的情况下用一阶导数法寻峰，识别出139个峰，如图2(a)所示。 对标准样品1实测谱线用均值平滑算法进行平滑并扣除本底后寻峰，当均值平滑窗口大小为5时，识别出43个峰，如图2(b)所示。 但从图2(b)可以观察到，峰9和峰10包含的重叠峰并没有识别出来。 可知： 一阶导数法具有较强的弱峰识别能力，但重叠峰分辨能力较差。

当均值平滑窗口大小为3，7，9，11，13，15，17和19时，识别出的峰个数分别为70，28，20，19，18，17，15和15。 可知，当均值平滑窗口逐渐增大时，识别出的峰逐渐减少。 且可看出，当均值平滑窗口较小时，平滑窗口大小变动对谱峰识别能力影响较大； 当均值平滑窗口较大时，平滑窗口大小变动对谱峰识别能力影响接近最大。

图2 一阶导数法寻峰效果(a)： 未平滑的寻峰结果； (b)：均值平滑后寻峰结果-窗口大小为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 连续小波变换法简介

由小波变换的性质可以将f的小波变换Wf(u,s)写成式(1)的形式。

(1)

式(1)中，Ψ为小波函数；θ为一使小波函数Ψ具有n阶消失矩的快速衰减函数。 从式(1)可以看出，信号光滑后的n阶导数可以通过小波变换来实现。 可见，小波变换具有自动去噪功能。 由于XRF的特征峰近似于高斯函数，高斯系列小波gausn常用于XRF的寻峰。 其中，gaus2，gaus4，gaus6和gaus8的作用类似于导数法寻峰中的偶数阶导数，常用于重峰分解，函数如图3。

由图3可知，随阶数的增加，高斯小波的形状变得更窄，即重峰分解能力更强。 但小波函数两侧的波瓣也变大，导致变换后的光谱中会出现干扰峰，不利于光谱中的弱峰识别。 小波函数gaus2两侧没有值大于0的波瓣，选用gaus2对光谱进行处理，不会增加直接正向的干扰峰。

2.2 连续小波变换法寻峰过程

连续小波变换寻峰过程： (1)对XRF进行连续小波变换，得到不同分解尺度的小波变换系数； (2)选择合适分解尺度的小波变换系数[11]，例如可以选择峰幅值最接近原始谱线峰幅值的小波变换系数谱； (3)从选定小波变换系数谱中寻峰。 用gaus2对标准样品1实测谱线进行连续小波变换，得到不同分解尺度的小波变换系数。 图4显示了分解尺度1～8的小波变换系数。

图3 偶数阶高斯小波函数Fig.3 Even-order Gauss wavelet functions

图4 分解尺度1～8时的小波变换系数Fig.4 Wavelet transform coefficients when the decomposition scale changes from 1 to 8

通过对不同分解尺度小波变换系数的分析，可知： (1)分解尺度较小时，小波变换系数含有大量噪声，但对重峰分解能力最强； (2)随着分解尺度的增加，系数逐渐变得光滑，分辨率逐渐变大，弱峰识别能力增强，但重峰识别能力变弱； (3)当分解尺度变得更大时，部分弱峰被平滑掉，逐渐变得不可识别。 从峰高可以看出图4中标准样品1实测谱线的小波变换系数d6最接近原始谱线。 对小波变换系数d6进行寻峰，结果如图5所示。 寻到的峰个数为33个。 对这些峰进行识别，识别出20种元素： Mg，Al，Si，P，S，Cl，K，Sn，Ca，Ti，V，Cr，Fe，Co，Ni，Cu，W，Zn，Ga和Th。 而其余9种元素并未被识别出来。

图5 连续小波变换法寻峰Fig.5 Peak-finding results by continuous wavelet transform method

从上述分析可知，选择最优分解尺度小波变换系数的常规小波变换寻峰法无法同时满足较强的重叠峰分辨能力和弱峰识别能力。

3 连续小波变换与导数法联合寻峰

3.1 连续小波变换与导数法联合寻峰过程

本文提出的连续小波变换与导数法联合寻峰的自适应算法步骤如下： (1)对原始光谱用小波变换迭代法扣除本底[6]； (2)对扣除本底后光谱用高斯小波gaus2进行变换，得到不同分解尺度的小波变换系数CFi(1≤i≤Scale，Scale为一最大分解尺度，例如32)； (3)对不同分解尺度的小波变换系数谱用一阶导数法寻峰，得到峰集合PKi(1≤i≤Scale)； (4)合并集合PKi(1≤i≤Scale)，得到UPK。 在合并过程中，以PK1为UPK初始值，然后依次将PK2，PK3，…，PKScale中的首次出现的峰加入UPK。 (5)由能量E与道址N之间的关系表达式计算出X荧光分析仪测量元素范围内轻元素的K系和重元素的L系特征X射线对应的道址； (6)根据集合UPK和元素特征X射线对应的道址，确定原始光谱所含元素并确定道址。

3.2 连续小波变换与导数法联合寻峰实例分析

按照上述算法步骤对标准样品1的实测谱线扣除本底，连续小波变换后得到不同分解尺度的小波变换系数CFi(1≤i≤Scale，Scale=32)。 然后对小波变换系数谱寻峰，步骤如下：

(1)小波变换系数谱寻峰。 对小波变换系数谱CFi寻峰，得到不同分解尺度峰集合PKi(1≤i≤Scale)，寻峰结果如表2所示。 分解尺度为1的小波变换系数谱寻峰结果如图6(a)所示，识别的峰个数为22。 从图6(a)可以看出，原始光谱中能够直观辨别的重叠峰均已被识别出。

表2 不同分解尺度的小波变换系数谱的寻峰结果

(2)合并峰集合。 合并不同分解尺度峰集合PKi(1≤i≤Scale)，得到的UPK含有107个峰，如图6(b)所示。 有些元素的特征X射线在UPK中有多条，它们是不同分解尺度下寻峰的结果，在后续步骤中会确定哪条该保留。

(3)确定X荧光分析仪的能量刻度。 由不同的标准样品得到能量E和道址N的关系式如式(2)所示。 计算出轻元素(原子序数11≤Z≤45的元素)的K系和重元素(原子序数45

E=0.007 599N+0.134 714 (keV)

(2)

(4)确定原始光谱所含元素。 根据集合UPK和元素特征X射线标准刻度，进行元素识别，得到结果如表3所示。 表中Scale值表示对应元素是在相应小波分解尺度下被识别的。

图6 自适应小波变换导数法寻峰(a)： 分解尺度为1的寻峰结果； (b)： 合并寻峰结果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 寻峰结果Table 3 Results of peak-finding

(5)寻峰结果分析。 由表1和表3比较后可知： 除元素Na，Zr和Nb外的26种元素都被确定。 X荧光分析仪在低能端的探测效率较低，原子序数靠前的Na元素无法被探测到。 当道址N为2 048时，由式(2)可以计算出该道址对应的能量为15.697 5(keV)，而Zr和Nb的特征X射线的Ka1值分别为15.774(keV)和16.614(keV)，超出了X荧光分析仪元素分析范围。 可以看出，含量比例高的元素在小波分解尺度为1和2时即可识别。 而含量比例低的元素，有些在小波分解尺度取较大值时才能识别，例如Ga，As和Pb元素。 因为一阶导数法对孤立的峰具有较高的灵敏度，含量比例极低的元素W在小波分解尺度为3时就被识别。

4 结 论

分析比较了导数法和连续小波变换法的XRF寻峰效果，总结了导数法和连续小波变换法在对XRF进行寻峰时的优缺点，并提出了一个综合利用高斯小波gaus2和一阶导数对XRF进行寻峰的自适应算法。 该算法对不同分解尺度的小波变换系数谱进行一阶导数寻峰，然后合并寻峰结果，最后进行元素识别。 本文提出的自适应寻峰算法具备一阶导数法的弱峰识别能力，小波变换寻峰法的重叠峰分辨能力和自动去噪能力。 新方法对所有分解尺度的小波变换系数谱进行寻峰并合并寻峰结果，优于只对某一分解尺度小波变换系数谱进行寻峰的常规小波变换寻峰法。 按照新算法得到的特征峰集合可以作为后续精确寻峰、峰面积计算工作的初始值。