吳希軍，張 杰，肖春艷，趙學(xué)亮, ，李 康，龐麗麗，史彥新，李少華

1. 燕山大學(xué)電氣工程學(xué)院河北省測(cè)試計(jì)量技術(shù)及儀器重點(diǎn)實(shí)驗(yàn)室，河北 秦皇島 066004 2. 河南理工大學(xué)資源與環(huán)境學(xué)院，河南 焦作 454000 3. 中國(guó)地質(zhì)調(diào)查局水文地質(zhì)環(huán)境地質(zhì)調(diào)查中心，自然資源部地質(zhì)環(huán)境監(jiān)測(cè)工程技術(shù)創(chuàng)新中心，河北 保定 071051 4. 河北先河環(huán)?？萍脊煞萦邢薰?，河北 石家莊 050000

引 言

土壤中的重金屬會(huì)對(duì)人體健康產(chǎn)生嚴(yán)重影響，對(duì)于重金屬的檢測(cè)顯得尤為重要[1]?，F(xiàn)有的土壤重金屬檢測(cè)技術(shù)有原子吸收光譜法、電感耦合等離子體法、原子熒光光譜法等[2]，這些檢測(cè)技術(shù)存在檢測(cè)時(shí)間長(zhǎng)、前期處理復(fù)雜、可能造成二次污染和成本較高的缺點(diǎn)，而X射線熒光光譜法具有檢測(cè)時(shí)間短、無污染且對(duì)多元素同時(shí)檢測(cè)等特點(diǎn)，更加滿足對(duì)土壤重金屬的檢測(cè)要求。

在實(shí)際的檢測(cè)中，土壤含水量、溫度、土壤粒徑大小都會(huì)對(duì)光譜的采集產(chǎn)生影響，其中土壤含水量會(huì)對(duì)光譜產(chǎn)生嚴(yán)重的干擾使得直接對(duì)土壤重金屬進(jìn)行估計(jì)精度較低。現(xiàn)有很多學(xué)者對(duì)去除含水量影響進(jìn)行了研究，Ji等[3]使用直接標(biāo)準(zhǔn)化法(direct standardization，DS)對(duì)濕土光譜進(jìn)行校正，使用校正后的光譜預(yù)測(cè)土壤屬性； 安曉飛等提出使用水分吸收指數(shù)對(duì)不同含水量樣本進(jìn)行分類并給與修正系數(shù)，來消除含水量的干擾； 胡曉艷等[4]提出了水分修正系數(shù)法(moisture determination index，MDI)有效降低了含水量的干擾，提高了利用干土土樣有機(jī)質(zhì)定量預(yù)測(cè)濕土有機(jī)質(zhì)的精度； Minasny等[5]使用外部參數(shù)正交化法(external parameter orthogonalization，EPO)去除近紅外光譜中土壤含水量的影響，對(duì)土壤有機(jī)碳含量進(jìn)行校正。這些算法具有易造成過度校正、忽視待估土壤屬性濃度與受含水量影響光譜之間關(guān)系等不足，且是對(duì)于去除含水量對(duì)土壤屬性預(yù)測(cè)精度提升的研究，對(duì)去除含水量影響提高土壤重金屬含量預(yù)測(cè)精度研究較少。

本研究采用非負(fù)矩陣分解算法(nonnegative matrix factorization，NMF)，對(duì)原始能譜矩陣進(jìn)行分段式分解，在保留與土壤重金屬有關(guān)的信息下，濾除能譜矩陣中因含水量引起的噪聲偏移，處理后建立基于偏最小二乘法預(yù)測(cè)(partial least squares regression，PLSR)的土壤重金屬含量反演模型，并對(duì)其進(jìn)行定量評(píng)估。

1 實(shí)驗(yàn)部分

1.1 儀器

CIT-3000SYB能量色散X射線熒光分析儀, 用于光譜檢測(cè)； 瑞紳葆PrepP-01 100T XRF大噸位壓片機(jī)，用于制作土壤樣品壓片； 金壇大地電動(dòng)攪拌器，用于水土充分?jǐn)嚢瑁?DHG-9141A電熱恒溫干燥箱，用于烘干原土壤中水分； 艾澤拉小型磨粉機(jī)，用于打磨土壤等。

1.2 試劑

土壤成分分析標(biāo)準(zhǔn)物質(zhì)GSS-4(GBW07404)、GSS-5(GBW07405)、GSS-7(GBW07407)、GSS-8(GBW07408)、GSS-15(GBW07429)和GSS-28(GBW07457)，購(gòu)自地球物理地球化學(xué)勘查研究所； 硝酸鎘分析純、硝酸鉻分析純、重金屬標(biāo)準(zhǔn)溶液(Pb，Zn，Cr，Cd，As)，購(gòu)自北京世紀(jì)奧科生物技術(shù)有限公司等。

1.3 樣品制備及測(cè)試

土壤樣本采樣于河北省保定市滿城區(qū)北莊村附近，采集耕層(0～20 cm)的土壤樣本。清除土樣中的雜物帶到實(shí)驗(yàn)室進(jìn)行烘干、用瑪瑙研缽碾壓后過10目(1 500 μm)的篩子； 再將土壤倒入研磨機(jī)進(jìn)行研磨，之后過200目(75 μm)的篩子，能過篩子的土壤倒入保鮮袋中備用。稱量不同質(zhì)量分析純和重金屬標(biāo)準(zhǔn)溶液，配制濃度為100，200，400，600和800 mg·kg-1不同含水量的重金屬溶液。華北平原的土壤含水量在15%左右，在土壤中加入配制好的重金屬溶液，制作出含水量10%～25%之間等步長(zhǎng)十個(gè)含水梯度的土壤樣品，將麥拉膜與樣品杯嵌套固定，裝入處理后的土壤并壓緊。

測(cè)定前先啟動(dòng)能量色散X熒光分析儀預(yù)熱一段時(shí)間，將標(biāo)準(zhǔn)樣放入儀器檢測(cè)臺(tái)的樣品腔中對(duì)儀器進(jìn)行校準(zhǔn)，使儀器處于最佳的工作狀態(tài)； 再將分析儀的探測(cè)窗對(duì)準(zhǔn)土壤樣品杯進(jìn)行檢測(cè)，檢測(cè)時(shí)間為200 s。

2 結(jié)果與討論

能量色散X射線熒光光譜檢測(cè)法原理是建立重金屬濃度與熒光強(qiáng)度之間的關(guān)系。本研究以重金屬鉛為例，基于莫塞來(Moseley)定律和普朗克方程可以通過Pb的X射線激發(fā)能量找到Pb的特征峰位置，以特征峰的凈峰面積為熒光強(qiáng)度。原始光譜經(jīng)過平滑去噪、本底扣除、基于NMF算法去除含水量影響后，對(duì)特征峰求積分得到凈峰面積，再使用偏最小二乘法建立重金屬含量和凈峰面積之間的關(guān)系。

2.1 異常樣本數(shù)據(jù)剔除

在實(shí)驗(yàn)過程中存在著人為因素例如儀器操作不當(dāng)、土壤樣品攪拌程度不夠等因素，導(dǎo)致小部分實(shí)驗(yàn)數(shù)據(jù)不準(zhǔn)確，使用這些偏差較大的數(shù)據(jù)進(jìn)行分析時(shí)會(huì)使得分析結(jié)果精度低，所以應(yīng)予以剔除。原始矩陣圖是在二維空間的坐標(biāo)系中以樣本點(diǎn)的形式體現(xiàn)數(shù)據(jù)的差異，計(jì)算了馬氏距離和使用了NJW聚類分析的方式，使原始矩陣圖切割成幾個(gè)子圖，使得幾個(gè)子圖間相似度最弱而每個(gè)子圖里樣本數(shù)據(jù)相似度強(qiáng)。圖1是相同含水量下的原始數(shù)據(jù)矩陣經(jīng)過NJW聚類分析結(jié)果圖。

圖1 NJW聚類分析結(jié)果 (a)： 原始樣本點(diǎn)； (b)： 聚類后樣本點(diǎn)Fig.1 NJW cluster analysis result (a): Original sample point; (b): Sample point after clustering

從聚類切割圖結(jié)果中發(fā)現(xiàn)了600 mg·kg-1中有一個(gè)數(shù)據(jù)點(diǎn)被劃分到400 mg·kg-1中不具有重復(fù)性，予以剔除。

2.2 土壤含水量對(duì)測(cè)定結(jié)果的影響

土壤中的水分會(huì)影響所制土壤樣品的均一性，而且水分會(huì)對(duì)源初級(jí)X射線和特征X射線產(chǎn)生吸收作用，加大X射線的散射效應(yīng)。故土壤中的水分的存在，降低了儀器對(duì)土壤中重金屬的特征峰強(qiáng)度，從而增加了儀器的檢測(cè)誤差。

每一含水量的測(cè)定光譜都是由相同重金屬含量同一含水量的多次測(cè)試結(jié)果進(jìn)行主成分分析(principal component analysis，PCA)加權(quán)處理得到的，可以準(zhǔn)確地反映這一含水量的真實(shí)光譜。圖2表示在重金屬含量800 mg·kg-1下重金屬Pb不同含水量的測(cè)定光譜，清楚顯示了水分對(duì)測(cè)試結(jié)果的影響，在重金屬特征吸收峰附近不同含水量之間的光譜相似性差、重疊程度低，隨著土壤含水量的增加檢測(cè)精度明顯下降，重金屬Pb的相對(duì)偏差從7.41%升高到27.78%。

圖2 不同含水量下濃度800 mg·kg-1重金屬Pb測(cè)定光譜Fig.2 Spectral determination spectra of 800 mg·kg-1 heavy metal Pb under different water contents

2.3 測(cè)定光譜預(yù)處理

在對(duì)樣本進(jìn)行測(cè)定獲取光譜的過程中由于外界環(huán)境、測(cè)定條件的差異會(huì)導(dǎo)致原始光譜中存在噪聲信號(hào)和基線漂移問題。光譜預(yù)處理的方法比較多，針對(duì)產(chǎn)生的噪聲信號(hào)，本研究使用了5點(diǎn)2次Savitzky-Golay卷積平滑去噪法，能夠非常有效地提高譜圖的信噪比； 針對(duì)儀器本身帶來的影響本研究采用了線性本底法扣除本底。

圖3是原始光譜經(jīng)過預(yù)處理的結(jié)果圖，(a)表示實(shí)際重金屬含量為400和600 mg·kg-1光譜經(jīng)過平滑后的效果圖，圖中可以看出光譜區(qū)間的噪聲明顯降低并且光譜的形狀和寬度沒有發(fā)生變化，說明這種方法可以去除噪聲，提高光譜的信噪比； (b)表示使用線性本底法進(jìn)行本底扣除，確定特征峰左右邊界計(jì)算本底面積，扣除本底面積。

2.4 基于NMF算法土樣含水量因素的去除

2.4.1 NMF算法

預(yù)處理后光譜的主要干擾是土壤含水量這一因素，常見的去除含水量影響算法有直接標(biāo)準(zhǔn)化法、外部參數(shù)正交化法，直接標(biāo)準(zhǔn)化法是對(duì)全譜端進(jìn)行整體處理易過擬合，外部參數(shù)正交化法忽視了變量前后之間的相互影響，本研究采用非負(fù)矩陣分解算法用于去除含水量的影響。

非負(fù)矩陣分解是由Lee和Seung獨(dú)立提出后，在圖像處理、人臉識(shí)別中得到廣泛的使用，較少應(yīng)用于X射線光譜分解和融合中。NMF算法是把一個(gè)非負(fù)矩陣分解成兩個(gè)非負(fù)矩陣因子的乘積，是一種有效的光譜分解的方法。

圖3 光譜預(yù)處理結(jié)果 (a): S-G平滑去噪； (b): 線性本底扣除本底Fig.3 Spectral preprocessing result (a): S-G smooth denoising; (b): Linear background minus background

理論上在相同重金屬濃度下光譜強(qiáng)度應(yīng)是幾乎一致的，由于不同含水量的影響，導(dǎo)致光譜強(qiáng)度存在明顯的差異。光譜數(shù)據(jù)的灰度值為非負(fù)的，所以可以建立一個(gè)初始非負(fù)矩陣用于光譜分解。將處理后的相同重金屬濃度不同含水量X射線熒光矩陣看作n個(gè)樣本m個(gè)能量段的非負(fù)數(shù)矩陣An×m，非負(fù)矩陣的行表示樣本數(shù)據(jù)，列表示能量端。在非負(fù)矩陣分解時(shí)應(yīng)在保證有效信息最大化的前提下進(jìn)行的，利用歐幾里得距離作為成本函數(shù)，用于量化分解精度即

(1)

先預(yù)設(shè)置端元數(shù)目r，其值小于行值、列值，對(duì)非負(fù)數(shù)矩陣An×m按照分解的數(shù)學(xué)表達(dá)式進(jìn)行分解

(2)

其中，Wn×r是端元光譜矩陣，其每一列代表著一個(gè)端元光譜，Hr×m是端元豐度矩陣，其每一列代表一個(gè)樣本不同端元的豐度值。

非負(fù)矩陣分解(NMF)數(shù)據(jù)分析具體操作過程如下：

(1)隨機(jī)生成端元光譜矩陣Wn×r；

(4)設(shè)置迭代次數(shù)，重復(fù)(2)和(3)的步驟直到成本函數(shù)不變或變化很小，迭代算法的收斂性在文獻(xiàn)[7]中已經(jīng)給出證明。 2.4.2 非負(fù)矩陣端元數(shù)目的確定

如果r過小會(huì)使得部分真實(shí)信息丟失，過大會(huì)達(dá)不到去除含水量這一因素影響的效果[8]，對(duì)于r的確定采用峰值信噪比(peak signal-to-noise ratio，PSNR)作為評(píng)價(jià)標(biāo)準(zhǔn)，峰值信噪比使用均方誤差(MSE)來定義[9]

(3)

其中，N表示樣本個(gè)數(shù)，下標(biāo)(i,k)表示第i個(gè)樣本第k個(gè)能量段表示的像素點(diǎn)，maxk表示第k個(gè)能量段的最大像素值。

圖4 r對(duì)峰值信噪比的影響Fig.4 Effect of r on peak signal to noise ratio

PSNR的值越大表示端元光譜矩陣Wn×r和端元豐度矩陣Hr×m融合成新的光譜矩陣的準(zhǔn)確度越高，在非負(fù)矩陣分解時(shí)端元數(shù)目設(shè)置從2到30，其間隔為1，圖4中可以清楚看出當(dāng)端元數(shù)目至10時(shí)趨于穩(wěn)定，后面數(shù)據(jù)波動(dòng)很小，所以端元數(shù)目設(shè)置為10。 2.4.3 NMF模型構(gòu)建

NMF算法的目的是獲取有用的光譜，去除含水量影響的干擾光譜。由非負(fù)矩陣分解得到端元光譜矩陣Wn×r和端元豐度矩陣Hr×m融合成新的光譜矩陣，新的光譜矩陣與原始光譜矩陣的差異矩陣為光譜差異矩陣。圖5表示了同一重金屬含量不同含水量下NMF算法處理前后的光譜。

當(dāng)r等于5時(shí)，除了21.69%和25%含水量較大的樣本其他樣本間光譜計(jì)數(shù)率曲線變化一致，光譜間的差異較小； 當(dāng)r等于10時(shí)，不同含水量樣本光譜間的差異微小，呈現(xiàn)出較好的相似性，重金屬特征峰附近差異也明顯消失，有效地降低了含水量對(duì)光譜的影響。

計(jì)算了樣本處理前含水量梯度兩兩之間的光譜數(shù)據(jù)相關(guān)系數(shù)，之間的相關(guān)系數(shù)相差較大，相關(guān)系數(shù)隨著含水量的增加而不斷減小，其中含水量在10%和25%之間的相關(guān)系數(shù)最小為0.72； 使用NMF算法處理后，樣本間的相關(guān)系數(shù)變化較小且接近于1，故進(jìn)一步說明了NMF算法消除了含水量對(duì)光譜的影響。

圖5 NMF算法處理前后的光譜 (a)： NMF處理前的光譜圖； (b)： NMF在端元數(shù)目r=5處理后的光譜圖； (c)： NMF在端元數(shù)目r=10處理后的光譜圖Fig.5 Spectrum before and after NMF algorithm processing

(a): Spectrum before NMF processing; (b): The spectrum of NMF processed by the number of end elementsr=5; (c): The spectrum of NMF after the number of end elementsr=10

2.5 PLSR，EPO-PLSR和 NMF-PLSR模型的建立與對(duì)比

偏最小二乘法預(yù)測(cè)(PLSR)是一種多變量統(tǒng)計(jì)分析方法，同時(shí)考慮了分解自變量光譜矩陣和目標(biāo)因變量重金屬含量矩陣的影響，將數(shù)據(jù)壓縮與回歸擬合相結(jié)合，使建立的模型具有更好的穩(wěn)健性，并能有效地解決自變量之間的多重共線性。本研究中使用了留一法交叉驗(yàn)證防止過擬合，主成分?jǐn)?shù)是由交叉驗(yàn)證決定的n=2。

圖6 PLSR模型、EPO-PLSR模型和NMF-PLSR模型 實(shí)測(cè)值和預(yù)測(cè)值對(duì)比 (a)： PLSR模型； (b)： EPO-PLSR模型； (c)： NMF-PLSR模型

將計(jì)算得到的凈峰面積和土壤含水量作為模型輸入，土壤重金屬Pb的含量作為輸出。將得到的150組數(shù)據(jù)集使用留出法進(jìn)行7/3分樣分成兩個(gè)互斥的集合，70%的數(shù)據(jù)即105組數(shù)據(jù)作為訓(xùn)練集，30%的數(shù)據(jù)即45組數(shù)據(jù)作為測(cè)試集。并采用EPO算法進(jìn)行去除含水量，建立了EPO-PLSR模型用作對(duì)比。模型的精度檢驗(yàn)及評(píng)價(jià)指標(biāo)使用建模決定系數(shù)R2、交叉驗(yàn)證均方根誤差RMSECV、平均絕對(duì)誤差MAE和相對(duì)分析誤差RPD，其中RPD是預(yù)測(cè)樣本標(biāo)準(zhǔn)差和預(yù)測(cè)均方根誤差RMSEP的比值，當(dāng)RPD<1.4時(shí)認(rèn)為模型預(yù)測(cè)是不可用的； 當(dāng)1.4≤RPD<2時(shí)認(rèn)為模型預(yù)測(cè)比較粗糙，精度不夠； 當(dāng)RPD≥2時(shí)認(rèn)為模型的預(yù)測(cè)能力較好[10]。

由圖6(a)和(b)對(duì)比可以看出，與PLSR模型相比，EPO-PLSR模型的R2和RPD分別提高了0.009 8和0.921 1，RMSECV和MAE分別降低了2.212 2和1.083，說明使用EPO-PLSR預(yù)測(cè)土壤重金屬的結(jié)果明顯優(yōu)于PLSR，相對(duì)于PLSR有了明顯改善，且PLSR的RPD值為1.53預(yù)測(cè)不夠準(zhǔn)確，只能實(shí)現(xiàn)在不同含水量下對(duì)土壤重金屬含量的粗略估算。經(jīng)過NMF-PLSR模型處理后的數(shù)據(jù)相對(duì)于PLSR模型改善較大，R2和RPD分別提高了0.019 7和1.029 2，RMSECV和MAE分別降低了2.386 3和1.439 6，相對(duì)于EPO-PLSR模型的R2和RPD分別提高了0.009 9和0.108 1，RMSECV和MAE分別降低了0.244 7和0.356 6，NMF-PLSR最佳擬合線更加接近1∶1線，說明了NMF-PLSR模型比EPO-PLSR模型預(yù)測(cè)能力強(qiáng)、模型解釋能力較好； 觀察預(yù)測(cè)值，NMF-PLSR模型在相同重金屬含量下預(yù)測(cè)值相差較小，進(jìn)一步說明了去除含水量的能力及預(yù)測(cè)的穩(wěn)定性。

3 結(jié) 論

通過實(shí)驗(yàn)的方法獲取了在相同重金屬含量條件下不同含水量土樣的光譜，分析了土壤含水量對(duì)光譜的影響，然后采用了非負(fù)數(shù)分解算法去除含水量的影響，利用峰值信噪比評(píng)價(jià)標(biāo)準(zhǔn)確定端元數(shù)目，并建立了NMF-PLSR模型對(duì)土壤重金屬含量進(jìn)行反演，得到了以下結(jié)論：

(1)隨著含水量的增加，土壤光譜反射率的非線性降低，相對(duì)偏差從7.41%升高到27.78%，儀器的檢測(cè)精度下降明顯，含水量對(duì)基于能量色散X射線光譜法的土壤重金屬含量反演影響顯著。

(2)當(dāng)端元數(shù)目為10時(shí)峰值信噪比趨于穩(wěn)定，經(jīng)NMF算法處理后土壤重金屬特征峰重疊明顯，光譜間的相似系數(shù)更加接近1，光譜間的差異得到了明顯去除。

(3)建立的NMF-PLSR模型相對(duì)于直接建立PLSR模型效果有顯著改善，對(duì)比EPO-PLSR模型也有一定的改善，R2和RPD分別提高了0.009 9和0.108 1，RMSECV和MAE分別降低了0.244 7和0.356 6，更加準(zhǔn)確實(shí)現(xiàn)土壤重金屬含量的反演，從而為土壤重金屬含量檢測(cè)提供了一定的技術(shù)支持。