陳 偉 周建斌 方 方 洪 旭 萬文杰,3 喻 杰,3 王懷平

1（成都工貿(mào)職業(yè)技術(shù)學(xué)院信息工程學(xué)院 成都 611731）

2（成都理工大學(xué)核技術(shù)與自動化工程學(xué)院 成都 610059）

3（四川新先達測控技術(shù)有限公司 成都 610052）

X射線熒光光譜（X-Ray Fluorescence，XRF）分析作為一種無損檢測技術(shù)，技術(shù)特點上適合于樣品中主、次、痕量多元素的同時測定，方便用于現(xiàn)場、原位和活體分析，已廣泛應(yīng)用于地質(zhì)、礦石、冶金、考古、生物、生命、環(huán)境、材料、毒性物品等領(lǐng)域［1］。在測定痕量成分時，由于背景值的影響，特征峰容易部分或全部被掩蓋，給核素信息提取帶來困難。XRF背景值產(chǎn)生的原因如下：其一是由樣品引起，原級X射線譜在樣品中產(chǎn)生康普頓散射和瑞利散射，其強度隨樣品成分而發(fā)生變化；其二是樣品產(chǎn)生的射線與儀器相互作用引起；宇宙射線、樣品中的放射性核素與電子線路也可產(chǎn)生背景值［2?5］。

通過改進前端探測方法和后期譜線處理技術(shù)可以提高XRF對痕量信息的處理能力。文獻［6］應(yīng)用靶材的康普頓散射線作內(nèi)標(biāo)，用能量色散X熒光分析（Energy-Dispersive X-ray Fluorescence，EDXRF）譜儀測定了從鈦到鐵和從鈀到鋇的痕量核素。文獻［7］通過對配有三維偏振光學(xué)的EDXRF光譜儀的測量條件進行優(yōu)化，完成了大豆樣品中Mg、P、Cl、K、Cu等微量元素的測定。文獻［8］利用XRF分析方法結(jié)合干灰化樣品前處理方法研究了不同品牌大米中微量元素的含量，同時用毛細管X光透鏡聚焦的微束X射線熒光掃描了單粒大米中多元素的分布。文獻［9］采用粉末壓片制樣，用波長色散X射線熒光光譜法對大米中的微量元素含量進行測定，討論了各目標(biāo)元素測量條件，并運用經(jīng)驗系數(shù)法和散射射線作內(nèi)標(biāo)對基體效應(yīng)進行校正。文獻［10］使用新的聚合物樹脂凝膠從樣品溶液中消除鈾基質(zhì)，通過全反射X射線熒光（Total reflection X-ray Fluorescence，TXRF）光譜測定法測定鈾樣品中的痕量元素。

在XRF分析中，譜線平滑可減小噪聲，有利于核素定性定量信息的提取。常用能譜平滑方法有算術(shù)滑動平均法、重心法、快速傅里葉變換法、高斯法、小波分析、B樣條法等?，F(xiàn)有能譜平滑方法各有優(yōu)勢，但不是在任何情況下都有效?；趯?shù)字Sallen-Key（S-K）濾波器的研究［11?14］，本文提出一種多重Sallen-Key（Multiple Sallen-Key，MSK）平滑技術(shù)，并將其用于EDXRF的實測譜線中痕量核素信息的提取。

1 多重數(shù)字S-K平滑原理

1.1 數(shù)字S-K濾波原理

Sallen-Key濾波器［15］在1955年由Sallen和Key提出，該濾波電路釆用部分正反饋控制，具有較高的品質(zhì)因素。在核電子學(xué)中，低通S-K濾波器廣泛用于核脈沖信號的濾波成形，可以在較少的級數(shù)下得到準(zhǔn)高斯波形的輸出。

文獻［11］對圖1中的低通S-K濾波器應(yīng)用基爾霍夫電流定律，建立了式（1）中的數(shù)學(xué)模型：

式中：y對應(yīng)圖1中的Vo；x對應(yīng)Vin；R1=R2=R；C1=C2=C；a=（R3+R4）/R3。式（1）為連續(xù)信號通過S-K濾波器時，輸出與輸入信號的關(guān)系。

圖1 低通S-K濾波電路原理圖Fig.1 Block diagram of S-K low-passfilter circuit

在數(shù)字信號處理時，需對連續(xù)信號做數(shù)字離散化。在模數(shù)轉(zhuǎn)換時，遵從奈奎斯特（Nyquist）采樣定理，采樣周期記為Ts，則第n個采樣點對應(yīng)的采樣時間為nTs，結(jié)合數(shù)值微分原理：

將式（2）和式（3）代入式（1）得：

式（4）是以固定的采樣頻率對連續(xù)信號做數(shù)字 離散后表達式，合并整理可得：

令：k=（RC）/Ts，上式變?yōu)椋?/p>

整理可得：

式 中：k=RC/Ts；a=（R3+R4）/R3；x（nTs）為輸 入 數(shù)據(jù)；y（nTs）為濾波后數(shù)據(jù)。k為濾波頻率因子，可改善輸出信號的平滑度，調(diào)節(jié)信號寬度；a為濾波幅值因子，用于調(diào)節(jié)輸出信號的幅值。

S-K濾波器具有較好的濾波效果［11?14］，濾波效果可通過k、a進行調(diào)節(jié)?，F(xiàn)有對數(shù)字S-K濾波的研究與應(yīng)用主要在信號處理方面，本文將數(shù)字S-K濾波應(yīng)用于核能譜平滑處理和痕量核素分析中。

1.2 MSK平滑原理

受S-K濾波器中RC電路等延遲因素影響，經(jīng)數(shù)字S-K濾波器處理后的數(shù)據(jù)會產(chǎn)生時間延遲，經(jīng)其處理后的核能譜數(shù)據(jù)會發(fā)生譜峰漂移（圖3）。發(fā)生譜峰漂移的譜線不利于在已經(jīng)刻度好的能譜工具中應(yīng)用；在高背景低含量核素的能譜分析中，一次譜線平滑無法達到有效特征峰信息提取的目的，為此，提出了多重（Multiple Sallen-Key，MSK）平滑（以下稱MSK平滑）。

MSK平滑的工作原理框圖如圖2所示。參數(shù)寄存器用于存儲平滑參數(shù)：平滑次數(shù)Times、k、a，第1緩存器，用于在平滑開始時暫存輸入的譜數(shù)據(jù)流，在平滑過程中暫存經(jīng)反向數(shù)字S-K濾波器后的數(shù)據(jù)流，并在平滑結(jié)束后輸出數(shù)據(jù)流；第2緩存器，用于暫存經(jīng)過正向數(shù)字S-K濾波器后的數(shù)據(jù)流；數(shù)據(jù)處理流程如下：平滑開始時，實測譜數(shù)據(jù)流送入第1緩存器；首先做正向平滑，第1緩存器中的譜數(shù)據(jù)流按道址由低向高順序通過數(shù)字S-K濾波器；然后再做反向平滑，方向轉(zhuǎn)換功能將第2緩存器中的譜數(shù)據(jù)流按道址由高向低送入反向數(shù)字S-K濾波器，處理后數(shù)據(jù)送入第1緩存器，如此便完成1次MSK平滑循環(huán)；MSK根據(jù)Times中的平滑次數(shù)，決定是否執(zhí)行下一次平滑循環(huán)。運行MSK平滑后的核能譜數(shù)據(jù)，譜峰位置不變，方便用原能譜的能量刻度信息進行數(shù)據(jù)分析。

圖2 MSK原理框圖Fig.2 Block diagram of MSK principle

取某大米樣品的X熒光實測能譜，設(shè)平滑參數(shù)為：k=1.02，a=1，做正向數(shù)字S-K平滑，平滑前后的譜數(shù)據(jù)如圖3所示。圖3中v1為實測譜，v2為正向SK平滑后譜線。分析圖3中的Ni Kα1與Ni Kβ1特征峰可知，執(zhí)行正向S-K平滑后，譜峰峰位向右漂移了2道。

圖3 譜峰位置對比圖Fig.3 Contrast chart of spectral peak position

取平滑參數(shù)：k=1.02，a=1，Times=1，對上述實測能譜數(shù)據(jù)做MSK平滑，平滑前后的譜線如圖4所示。圖4中v1為實測譜數(shù)據(jù)（以下簡稱原譜），v2為MSK平滑后譜數(shù)據(jù)（以下簡稱平滑譜）。由圖4可以看出，Ni Kα1與Ni Kβ1特征峰可知，執(zhí)行MSK平滑后，特征峰峰位保持不變。

圖4 MSK峰位移動對比圖Fig.4 Contrast chart of MSK spectral peak position movement

在現(xiàn)場數(shù)據(jù)分析中，對已經(jīng)刻度好的能譜分析工具，譜峰漂移不利于數(shù)據(jù)分析。尤其在分析高背景低含量核素時，常規(guī)譜線平滑技術(shù)已無法有效處理核能譜的統(tǒng)計漲落，MSK可以有效地解決這一問題。

2 實驗與分析

鎘是一種對人體有極大危害的重金屬，主要通過食物鏈在體內(nèi)富集，GB2762-2017《食品安全國家標(biāo)準(zhǔn)食品中污染物限量》中規(guī)定了食品中的鎘限量指標(biāo)，其中大米中的鎘含量≤0.2 mg·kg?1。大米中元素的分析手段有原子吸收光譜技術(shù)、電感耦合等離子體質(zhì)譜技術(shù)、原子熒光光譜技術(shù)、中子活化分析技術(shù)、氣相色譜－質(zhì)譜技術(shù)、X熒光分析法等，其中X熒光分析法因?qū)嶒灄l件要求和測量方法簡單快速，適于對大米中重金屬鎘的快速篩查。

樣品制備方法：將干燥的大米樣品用家用研磨機磨碎，在30 MPa的壓力下制成直徑為25 mm、高度為20 mm的圓柱形餅樣，用食品保鮮膜包裹并標(biāo)記。以下實驗所用大米標(biāo)準(zhǔn)樣品中的鎘含量，是依據(jù)GB5009.15-2014，采用石墨爐原子吸收光譜法測定的；實驗中的X射線能譜采用自主研發(fā)的EDXRF型X熒光分析儀測量，儀器的檢出限為0.04 mg·kg?1。譜數(shù)據(jù)平滑處理在MATLAB 2016b中完成。

2.1 MSK平滑譜峰分析

熒光儀的X-tube與樣品間采用Ni濾片，取含量為0.19 mg·kg?1的鎘大米標(biāo)準(zhǔn)樣品，測量時間為1 800 s。取平滑參數(shù)k=1.02，a=1，Times值為依次取1、2、3、5、7、9，對上述實測譜線做MSK平滑。鎘核素的感興趣區(qū)（Region of Interest，ROI）取1 652~1 722道，以下所有ROI值與此相同。對MSK平滑前后的ROI數(shù)據(jù)做瀑布圖，結(jié)果如圖5所示。

圖5 不同次數(shù)的MSK平滑處理結(jié)果Fig.5 Resultsof MSK smoothing with different times

由圖5可知，Times參數(shù)對含背景單能高斯峰的峰形影響較大。當(dāng)Times值較小時，平滑譜中可見因背景干擾引起的弱峰；隨著Times值的增大，平滑譜的峰形雖有一定的壓縮與展寬，但峰位維持不變。譜峰形狀及光滑度變好。與原譜相比，經(jīng)MSK平滑后，譜線中的噪聲得到較好地抑制。

2.2 多次高斯平滑與多重高斯平滑對比分析

高斯平滑也稱高斯濾波，是一種按照高斯函數(shù)的特性來計算權(quán)值的線性平滑濾波器，其實質(zhì)是高斯濾波模板與測量信號的卷積。高斯平滑的主要參數(shù)有高斯模板窗寬和標(biāo)準(zhǔn)差。高斯平滑對抑制服從正態(tài)分布的噪聲非常有效，廣泛應(yīng)用于圖形圖像和能譜數(shù)據(jù)輸出處理。

對前述實驗中0.19 mg·kg?1的鎘大米標(biāo)準(zhǔn)樣實測譜線，先做MSK平滑，平滑參數(shù)為：k=1.02，a=1，Times=9；再對原譜做多次高斯平滑（Multi-scale Gaussian Smoothing，MGS）［16］，平滑參數(shù)為：高斯模板窗寬為7，標(biāo)準(zhǔn)差為1.22，平滑次數(shù)為50。選取鎘元素ROI區(qū)域作圖（圖6）。圖6中，I區(qū)為鎘核素的ROI譜線放大圖，v1為實測譜，v2為多次高斯平滑后的平滑譜，v3為MSK平滑譜。

圖6 鎘大米實測能譜及其平滑對比圖Fig.6 Measured spectrum of cadmium riceand itssmoothing contrast chart

由圖6可以看出，兩種平滑方法所得的譜線高度符合。取經(jīng)兩種平滑方法處理后的ROI區(qū)能譜數(shù)據(jù)做線性擬合分析，結(jié)果如圖7所示?？梢钥闯?，兩種方法得到結(jié)果的擬合優(yōu)度R2為0.999 99。說明只要平滑參數(shù)選取適當(dāng)，同一能譜分別經(jīng)MSK平滑和多次高斯平滑處理后，可達到高度相近的處理效果。

圖7 高斯平滑與MSK平滑在ROI區(qū)的線性符合分析Fig.7 Liner analysisof MSK smoothing and Gaussian smoothing in ROI

在MATLAB中分別記錄MSK平滑與高斯平滑的運算時間。得到MSK平滑運算時間1.5 ms，高斯平滑運算時間11.8 ms。結(jié)合兩種方法平滑效果和高斯平滑的原理可知，在相同平滑效果下，MSK平滑次數(shù)更少，代碼執(zhí)行效率更高。同時，由于是迭代式算法，MSK更容易在嵌入式系統(tǒng)中實現(xiàn)［17］。

2.3 標(biāo)準(zhǔn)樣線性刻度與分析

取§2.2實驗中MSK平滑前后的譜線，分別用統(tǒng)計敏感的非線性迭代剝峰算法（Statistics-sensitive Nonlinear Iterative Peak-clipping，SNIP）扣除背景（圖8）。圖8中v1為扣SNIP背景后的實測譜，v2為先經(jīng)MSK處理再扣除SNIP背景后的平滑譜；圖8中I區(qū)為鎘核素的ROI譜線放大圖，可以直觀看出，對于峰背比較低的特征峰，利用MSK平滑方法處理后的能譜，可以有效抑制統(tǒng)計漲落，特征峰更有利于核素的定性定量信息提取。

圖8 扣除SNIP背景后的0.19 mg·kg?1能譜Fig.8 The0.19 mg·kg?1 spectrum after eliminating SNIP background

依據(jù)GB5009.15-2014中的石墨爐原子吸收光譜測定方法得到標(biāo)準(zhǔn)樣大米中鎘的含量值，然后取0.110 mg·kg?1、0.260 mg·kg?1、0.446 mg·kg?1、0.640 mg·kg?1、0.768 mg·kg?1的標(biāo)準(zhǔn)樣，采用EDXRF X熒光分析儀做現(xiàn)場能譜測量，測量時間為1 800 s。設(shè)平滑參數(shù)：k=1.02，a=1，Times=9，對上述實測譜線做MSK平滑。選取上述各標(biāo)準(zhǔn)樣的實測譜和MSK平滑譜，扣除SNIP背景后，計算各譜線中的鎘核素ROI凈峰面積，統(tǒng)計結(jié)果列于表1。

表1 實測標(biāo)準(zhǔn)樣的ROI峰面積Table1 ROIpeak area of measured referencematerial

取表1中的標(biāo)準(zhǔn)樣含量值為y坐標(biāo)，原譜凈峰面積、MSK處理譜凈峰面積分別為x坐標(biāo)，對其進行線性擬合分析，可得工作曲線：

式中：yi為被測核素含量；b為工作曲線截距；xi為被測元素ROI峰面積。線性擬合結(jié)果見圖9。分析可知：標(biāo)樣含量值與原譜凈峰面積之間的關(guān)系為yi=0.230 66xi?0.330 88，R2=0.982 77，擬合曲線見v1；標(biāo)樣含量值與平滑譜凈峰面積間的工作曲線為yi=0.226 01xi?0.013 487，R2=0.998 24，擬合曲線見v2；分析擬合優(yōu)度R2指標(biāo)可知，MSK處理提高了刻度曲線的線性。

圖9 標(biāo)準(zhǔn)樣含量值與ROI峰面積間的線性擬合結(jié)果Fig.9 Linear fitting result between referencematerials content valueand ROIpeak area

取采用石墨爐原子吸收光譜法測量所得含量值為0.100 mg·kg?1、0.130 mg·kg?1、0.160 mg·kg?1、0.190 mg·kg?1、0.360 mg·kg?1、0.519 mg·kg?1、0.563 mg·kg?1的鎘大米標(biāo)準(zhǔn)樣（以下稱標(biāo)準(zhǔn)樣含量值為樣品標(biāo)稱值）做X熒光現(xiàn)場測量，測量后譜線做MSK平滑，平滑參數(shù)與刻度實驗時相同，處理后的ROI凈峰面積代入式（8）中對應(yīng)參數(shù)的工作曲線，求得的被測樣的回歸分析含量值，以及各含量值與樣品標(biāo)稱值之間的相對誤差見表2。分析表2可知，經(jīng)MSK平滑處理后的含量，與標(biāo)準(zhǔn)含量值間的相對誤差減小，說明MSK處理可以提高樣品含量值的測量精度。

表2 標(biāo)樣含量的回歸分析及其相對誤差Table 2 Regression analysis value and relative error of reference materials contents(mg·kg?1)

綜上可得，核能譜在MSK平滑處理后，其痕量核素的ROI峰面積與含量值之間仍能維持良好的線性。MSK平滑處理可以減小被測樣含量提取中的誤差，提高測量精度。本文提出的基于MSK平滑的痕量核素定量信息提取方法具有很好的應(yīng)用價值。

3 結(jié)語

本文通過對模擬S-K電路的數(shù)字化分析，提出MSK譜線處理方法，并將之用于鎘大米含量分析。MSK方法的參數(shù)設(shè)置方便，可有效抑制X熒光測量中高背景對痕量核素特征峰定量信息提取時的影響。實驗表明，MSK能夠以更簡潔的遞推運算和更短的代碼執(zhí)行時間，實現(xiàn)與多次高斯平滑方法相同的平滑效果。

取鎘大米標(biāo)準(zhǔn)樣的XRF實測譜線，經(jīng)MSK方法處理后，可以提高工作曲線的線性。取另一批大米標(biāo)準(zhǔn)樣做XRF分析，把處理后的鎘ROI峰面積代入上述工作曲線中求XRF分析含量值，將該值與采用GB5009.15-2014方法所得的標(biāo)準(zhǔn)含量值做誤差分析，可得MSK處理方法能減小的測量誤差，提高樣品的測量精度。

作者貢獻聲明陳偉負責(zé)算法、文章的起草和最終版本的修訂；周建斌、方方負責(zé)研究的提出及設(shè)計；洪旭負責(zé)算法程序的設(shè)計；萬文杰、喻杰負責(zé)圖表的設(shè)計與優(yōu)化；王懷平負責(zé)資料的收集及整理。