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

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

1. 燕山大學電氣工程學院河北省測試計量技術及儀器重點實驗室，河北 秦皇島 066004 2. 河南理工大學資源與環境學院，河南 焦作 454000 3. 中國地質調查局水文地質環境地質調查中心，自然資源部地質環境監測工程技術創新中心，河北 保定 071051 4. 河北先河環?？萍脊煞萦邢薰?，河北 石家莊 050000

引 言

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

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

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

1 實驗部分

1.1 儀器

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

1.2 試劑

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

1.3 樣品制備及測試

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

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

2 結果與討論

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

2.1 異常樣本數據剔除

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

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

從聚類切割圖結果中發現了600 mg·kg-1中有一個數據點被劃分到400 mg·kg-1中不具有重復性，予以剔除。

2.2 土壤含水量對測定結果的影響

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

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

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

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

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

2.4.1 NMF算法

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

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

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

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

(1)

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

(2)

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

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

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

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

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

(3)

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

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

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

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

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

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

圖5 NMF算法處理前后的光譜 (a)： NMF處理前的光譜圖； (b)： NMF在端元數目r=5處理后的光譜圖； (c)： NMF在端元數目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)是一種多變量統計分析方法，同時考慮了分解自變量光譜矩陣和目標因變量重金屬含量矩陣的影響，將數據壓縮與回歸擬合相結合，使建立的模型具有更好的穩健性，并能有效地解決自變量之間的多重共線性。本研究中使用了留一法交叉驗證防止過擬合，主成分數是由交叉驗證決定的n=2。

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

將計算得到的凈峰面積和土壤含水量作為模型輸入，土壤重金屬Pb的含量作為輸出。將得到的150組數據集使用留出法進行7/3分樣分成兩個互斥的集合，70%的數據即105組數據作為訓練集，30%的數據即45組數據作為測試集。并采用EPO算法進行去除含水量，建立了EPO-PLSR模型用作對比。模型的精度檢驗及評價指標使用建模決定系數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預測土壤重金屬的結果明顯優于PLSR，相對于PLSR有了明顯改善，且PLSR的RPD值為1.53預測不夠準確，只能實現在不同含水量下對土壤重金屬含量的粗略估算。經過NMF-PLSR模型處理后的數據相對于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模型在相同重金屬含量下預測值相差較小，進一步說明了去除含水量的能力及預測的穩定性。

3 結 論

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

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

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

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