基于NMF-PLS對含水量影響下土壤重金屬含量反演模型研究

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

1. 燕山大學電氣工程學院河北省測試計量技術(shù)及儀器重點實驗室，河北 秦皇島 066004 2. 河南理工大學資源與環(huán)境學院，河南 焦作 454000 3. 中國地質(zhì)調(diào)查局水文地質(zhì)環(huán)境地質(zhì)調(diào)查中心，自然資源部地質(zhì)環(huán)境監(jiān)測工程技術(shù)創(chuàng)新中心，河北 保定 071051 4. 河北先河環(huán)保科技股份有限公司，河北 石家莊 050000

引 言

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

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

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

1 實驗部分

1.1 儀器

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

1.2 試劑

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

1.3 樣品制備及測試

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

測定前先啟動能量色散X熒光分析儀預熱一段時間，將標準樣放入儀器檢測臺的樣品腔中對儀器進行校準，使儀器處于最佳的工作狀態(tài)； 再將分析儀的探測窗對準土壤樣品杯進行檢測，檢測時間為200 s。

2 結(jié)果與討論

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

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

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

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

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

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

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

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

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

2.3 測定光譜預處理

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

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

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

2.4.1 NMF算法

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

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

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

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

(1)

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

(2)

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

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

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

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

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

(3)

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

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

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

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

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

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

圖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模型的建立與對比

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

圖6 PLSR模型、EPO-PLSR模型和NMF-PLSR模型 實測值和預測值對比 (a)： PLSR模型； (b)： EPO-PLSR模型； (c)： NMF-PLSR模型

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

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

3 結(jié) 論

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

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

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

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