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

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

表1

Table 1

表1(Table 1)

表1 标准样品1中元素含量 Table 1 The content of elements in standard sample 1 元素 含量/% 元素 含量/% 元素 含量/% 元素 含量/% 元素 含量/% Na 0.35 Cl 0.003 7 Mn 0.017 As 0.000 14 Pb 0.000 87 MgO 2.01 K 4.16 Fe 7.6 Sr 0.009 Nb 0.001 43 SiO2 59.23 Ca 0.6 Co 0.002 1 Zr 0.009 6 Th 0.001 28 Al 18.82 Ti 0.4 Ni 0.003 7 Sn 0.000 2 Rb 0.020 5 P 0.069 V 0.008 8 Cu 0.004 2 Ba 0.045 Ga 0.002 6 S 0.006 Cr 0.009 9 Zn 0.005 5 W 0.000 079

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

图1

Fig.1

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

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

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图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

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

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

Wf(u, s)=f* Ψ¯s(u)=f* sndnθ¯s(t)dtn=sndndun(f* θ¯s)(u)(1)

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

图3

Fig.3

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

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

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

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图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

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

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

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

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

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

表2

Table 2

表2(Table 2)

表2 不同分解尺度的小波变换系数谱的寻峰结果 Table 2 Peak-finding results for wavelet transform coefficient spectrum with different decomposition scale 分解 尺度 峰数量 分解 尺度 峰数量 分解 尺度 峰数量 分解 尺度 峰数量 1 22 9 30 17 24 25 19 2 36 10 32 18 29 26 20 3 38 11 28 19 26 27 28 4 45 12 28 20 28 28 22 5 37 13 27 21 30 29 17 6 33 14 28 22 25 30 23 7 32 15 29 23 22 31 24 8 32 16 28 24 33 32 20

表2 不同分解尺度的小波变换系数谱的寻峰结果 Table 2 Peak-finding results for wavelet transform coefficient spectrum with different decomposition scale

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图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

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

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

E=0.007599N+0.134714(keV)(2)

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

表3

Table 3

表3(Table 3)

表3 寻峰结果 Table 3 Results of peak-finding 道址 对应 元素 Scale 道址 对应 元素 Scale 道址 对应 元素 Scale 142 Mg 2 569 Ba 1 1 090 W 3 179 Al 1 575 Ti 1 1 117 Zn 2 210 Si 1 632 V 2 1 202 Ga 21 247 P 4 692 Cr 2 1 382 As 30 299 S 4 757 Mn 2 1 387 Pb 30 328 Cl 8 823 Fe 1 1 660 Th 6 419 K 1 902 Co 1 1 744 Rb 10 452 Sn 1 963 Ni 3 1 911 Sr 10 466 Ca 1 1 047 Cu 9

表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射线的K_a1值分别为15.774(keV)和16.614(keV), 超出了X荧光分析仪元素分析范围。 可以看出, 含量比例高的元素在小波分解尺度为1和2时即可识别。 而含量比例低的元素, 有些在小波分解尺度取较大值时才能识别, 例如Ga, As和Pb元素。 因为一阶导数法对孤立的峰具有较高的灵敏度, 含量比例极低的元素W在小波分解尺度为3时就被识别。