基于X射線熒光光譜與多特征串聯(lián)策略的土壤重金屬含量預(yù)測

任 順，張 雄，任 東，楊信廷,3，張 力

(1.三峽大學(xué) 計算機與信息學(xué)院，湖北 宜昌 443002；2.三峽大學(xué) 湖北省農(nóng)田環(huán)境監(jiān)測工程技術(shù)研究中心，湖北 宜昌 443002；3.農(nóng)產(chǎn)品質(zhì)量安全追溯技術(shù)及應(yīng)用國家工程實驗室，北京 100097)

近年來，有毒重金屬對我國土壤環(huán)境質(zhì)量危害日益嚴重，耕地、農(nóng)產(chǎn)品等污染問題亟待解決,因此，重金屬監(jiān)測已成為農(nóng)產(chǎn)品保護和安全生產(chǎn)的重要任務(wù)[1]。目前，土壤重金屬污染檢測方法有傳統(tǒng)實驗室檢測法與速測法，其中傳統(tǒng)實驗室檢測方法主要有原子吸收與原子熒光光譜法、電感耦合等離子質(zhì)譜法等，該類方法雖然準確度和精確度高，但儀器運行條件要求高且不適用于現(xiàn)場快速判斷；速測法有X射線熒光光譜法(XRF)、生物傳感器技術(shù)和激光誘導(dǎo)擊穿光譜法[2-3]等。其中，XRF具有檢測速度快、成本低且可同時檢測多種元素等優(yōu)點，已被廣泛用于土壤重金屬污染和農(nóng)產(chǎn)品檢測。相比實驗室測定方法，XRF法前處理簡單，對主要的重金屬可有效監(jiān)測和快速篩查，可為農(nóng)田和農(nóng)作物重金屬污染的監(jiān)測和防控、制定合理的農(nóng)業(yè)發(fā)展規(guī)劃提供科學(xué)依據(jù)[4]。XRF法快速檢測重金屬時，針對不同類型和性質(zhì)的土壤，需建立重金屬濃度預(yù)測模型，再根據(jù)實際檢測結(jié)果對預(yù)測模型優(yōu)化。由于土壤成分復(fù)雜，采用XRF儀獲取的光譜數(shù)據(jù)具有很高的空間復(fù)雜程度。此外，所得到的光譜數(shù)據(jù)往往由數(shù)以千計的波長點組成，而樣本數(shù)量常受制于實驗條件和成本，在進行化學(xué)建模算法時，常會出現(xiàn)統(tǒng)計學(xué)稱為維數(shù)災(zāi)禍的多項式復(fù)雜程度的非確定性(NP)難題。模型集群分析策略強調(diào)，要最大限度分析已有樣本集信息，通過隨機采樣利用大量子集模型信息，獲得數(shù)據(jù)集的內(nèi)在結(jié)構(gòu)，而不僅僅依靠單一模型信息，可更好地避免模型對樣本的依賴[5-7]。

本研究以土壤重金屬的X射線熒光光譜為研究對象，將0～26 keV范圍內(nèi)的光譜數(shù)據(jù)與濃度梯度法配制的銅(Cu)、鋅(Zn)、砷(As)、鉛(Pb)、鉻(Cr)5種重金屬污染土壤樣本含量值進行關(guān)聯(lián)，將多特征串聯(lián)策略區(qū)間組合優(yōu)化-競爭適應(yīng)性重加權(quán)采樣-連續(xù)投影(ICO-CARS-SPA)算法所提取出的特征變量作為輸入，建立偏最小二乘(PLS)模型預(yù)測土壤樣本中5種重金屬含量，并與一些常見變量選擇方法及其組合對比，結(jié)合光譜特征變量的優(yōu)選和模型的預(yù)測性能，實現(xiàn)了對農(nóng)田土壤重金屬含量高效、精準的定量檢測，可為土壤重金屬治理和污染防控以及制定科學(xué)的農(nóng)業(yè)發(fā)展規(guī)劃提供依據(jù)。

圖1 ICO算法流程示意圖Fig.1 Flowchart of ICO algorithm

1 原理與算法

1.1 試驗區(qū)間組合優(yōu)化算法

區(qū)間組合優(yōu)化(ICO)是基于模型集群分析策略(MPA)框架下的化學(xué)建模算法[8]，MPA強調(diào)從數(shù)據(jù)集中統(tǒng)計分析時，應(yīng)通過隨機采樣，分析大量隨機產(chǎn)生的子集模型信息，從不同角度考察數(shù)據(jù)集的內(nèi)在性質(zhì)，獲得數(shù)據(jù)集的內(nèi)在結(jié)構(gòu)。ICO算法流程如圖1所示。

1.2 競爭性自適應(yīng)重加權(quán)采樣算法

競爭性自適應(yīng)重加權(quán)采樣(CARS)算法選擇波數(shù)的方法是基于PLS模型回歸系數(shù)，衡量波長選擇的重要依據(jù)是PLS回歸系數(shù)絕對值的大小，通過挑選出交叉驗證均方根誤差(RMSECV)值最小的子集來對應(yīng)獲得與目標值最優(yōu)的波長組合。計算則采用指數(shù)衰減函數(shù)，逐次保留絕對值大的回歸系數(shù)對應(yīng)的波數(shù)點[9-10]。

1.3 連續(xù)投影算法

連續(xù)投影算法(SPA)是一種前向變量循環(huán)特征波長選取方法，該算法利用向量的投影分析，能找到含有最低限度冗余信息的變量組。SPA一般能夠較為有效地消除波長之間共線性，同時提高模型精確度[11-12]。

1.4 多特征串聯(lián)算法

X射線熒光光譜的高維度無法規(guī)避地含有較多噪聲和冗余。ICO采用波長區(qū)間進行變量優(yōu)選，能排除大量無用信息和噪聲波長點，且在每次優(yōu)化后的聯(lián)合區(qū)間上進行局部搜索策略，更好地體現(xiàn)了柔性收縮策略優(yōu)勢，減少偶然誤差的發(fā)生，最終可得到一組位置、組合、寬度等都經(jīng)過優(yōu)化的有效波長區(qū)間。但ICO算法選擇的波長數(shù)量一般較多，區(qū)間內(nèi)部不可避免地存在共線性和冗余問題，最終選中的變量區(qū)間可采用單一變量選擇算法進行精簡[13]。

CARS算法作為一種高效的波長選擇算法，通過選擇少量波長變量子集，可以給出更佳或相當?shù)念A(yù)測效果。由于該算法在迭代過程中同時引入了蒙特卡羅采樣(MCS)和自適應(yīng)加權(quán)采樣(ARS)兩個隨機因素，其單獨使用的穩(wěn)定性難以令人滿意[14-15]。

SPA算法生成的每一個波長組合中，新入選的波長點與上一個波長點的相關(guān)性均為最低，這也體現(xiàn)了SPA算法能夠有效消除波長之間的共線性[16]。但對于X射線熒光光譜而言，有效變量間的投影距離并不一定最大，通過SPA算法挑選出變量子集很可能存在部分無關(guān)信息變量甚至是干擾變量。

將ICO選擇的優(yōu)化組合區(qū)間結(jié)合與競爭適應(yīng)性重加權(quán)采樣法(CARS)和連續(xù)投影算法(SPA)串聯(lián)(ICO-CARS-SPA)進行特征波長提取，既能有效地鎖定特征波段區(qū)間，同時優(yōu)選出的最優(yōu)波長子集又能夠充分減少最優(yōu)區(qū)間內(nèi)的共線性和冗余，形成優(yōu)勢互補。

2 實驗部分

2.1 材 料

在周邊1 km范圍內(nèi)無污染源的農(nóng)田采集制樣土壤，為防止實驗器具影響和干擾樣本重金屬的濃度，樣本制作過程采用陶瓷用具。將土壤烘干研磨過0.45 mm孔篩后，參考《土壤環(huán)境質(zhì)量標準》[17]中對Ⅰ、Ⅱ、Ⅲ類土壤中各元素含量的規(guī)定并結(jié)合農(nóng)田土壤實際情況制作實驗樣本。每個樣本制作30 g土壤，共91個樣本。經(jīng)初步測量，所采集農(nóng)田土壤中As、Cu、Cr、Pb、Zn的含量依次為3.0、15、30、18、39 mg/kg。綜合國家土壤等級劃分標準和所采集土壤的具體情況，設(shè)I為各個元素初始濃度梯度，其中1～40號樣本的濃度梯度為I；41～50號樣本的濃度梯度為2*I；51～60號樣本濃度梯度為3*I；61～70號樣本濃度梯度為4*I；71～91號樣本濃度梯度為J*I；其中J=1,2,3,4,……,n，具體實施方案參照表1。按土壤需求濃度，將重金屬標準溶液用丙酮稀釋后混入土壤中，混合均勻后置于通風(fēng)櫥中自然揮發(fā)，待土壤完全風(fēng)干后，再研磨混勻，按指定標號放入樣品盒中，密封保存。

表1 土壤樣本配制方案Table 1 Soil sample allocation scheme

2.2 光譜數(shù)據(jù)采集

每次取一個樣本置于便攜式X射線土壤重金屬檢測儀(由三峽大學(xué)協(xié)同創(chuàng)新中心與北京農(nóng)業(yè)質(zhì)量標準與檢測技術(shù)研究中心聯(lián)合研制)上，獲取樣本在0～26 keV范圍內(nèi)共4 096個通道內(nèi)的光譜信息，每個樣本經(jīng)過轉(zhuǎn)動方向測量3次并求取其平均光譜。由于36和37號樣本的光譜數(shù)據(jù)遺失，同時根據(jù)配制的樣品濃度繪制91個樣本的實際濃度曲線，發(fā)現(xiàn)60號樣品中的Zn濃度和79號樣品中Pb濃度分布與配制濃度差異較大，分析原因可能是溶液配制時出現(xiàn)了問題，故剔除這4個樣本。剩余87個樣本的平均光譜圖見圖2，其中配置Ag靶微型X光管、電子冷卻Si-PIN探測器，設(shè)置電壓為30 kV，電流為40 μA，積分時間為300 ms。

圖2 87個樣本的平均光譜圖Fig.2 Average spectrogram of 87 samples

2.3 建模方法與評價

2.3.1 光譜數(shù)據(jù)集的劃分該方法以樣本集被測指標的理化參照值作為劃分標準，將剔除4個異常樣本后的87個不同濃度重金屬的土壤樣本通過濃度梯度法進行校正集和驗證集的劃分，將樣本按照2∶1分成2組，得校正集58個，驗證集29個，樣本分布如圖3。

圖3 校正集(A)和驗證集(B)樣本分布圖Fig.3 Distribution figure of calibration set(A) and prediction set(B)

3 結(jié)果與討論

3.1 ICO算法選擇特征波長子區(qū)間

采用ICO算法進行特征變量選擇，并建立PLS模型預(yù)測重金屬含量。不同區(qū)間劃分下試驗結(jié)果如表2所示，當子區(qū)間劃分數(shù)為20時，5種重金屬均獲得最小的均方根誤差值。

表2 不同子區(qū)間數(shù)PLS建模對比Table 2 Comparison of PLS modeling for different subinterval numbers

(續(xù)表2)

以生物毒性顯著的重金屬元素Cr為特征提取對象的結(jié)果如下：采用ICO算法，進行最優(yōu)聯(lián)合區(qū)間的選取，利用5折交叉驗證法建立PLS模型選擇特征變量。選擇最佳波長區(qū)間劃分數(shù)量20，PLS模型中最大潛變量數(shù)為10，所選子模型的比例為0.05，加權(quán)自舉采樣次數(shù)為1 000。在ICO算法迭代過程中，隨著迭代次數(shù)的增加，每個波長區(qū)間的采樣權(quán)重變化情況見圖4A，采樣權(quán)重值越接近1時，顏色越接近深紅色；采樣權(quán)重值越接近0時，顏色越接近深藍色；顏色介于深藍色和深紅色之間表明采樣權(quán)重值處于0和1之間。圖中還可見，第6個波長區(qū)間的采樣權(quán)重值在第2次迭代過程中已經(jīng)約為1，其采樣權(quán)重值在第3次迭代過程中依然有機會變得小于1，最終該區(qū)間由于權(quán)重過小被剔除；第9個波長區(qū)間在第1～5次迭代過程中顏色偏綠藍，權(quán)重系數(shù)約為0.3～0.5，但在第6次迭代過程中權(quán)重系數(shù)仍有機會變大，最終被選中。由ICO每次迭代中提取的子模型的RMSECV圖可見，經(jīng)過11次迭代后均方根誤差趨于穩(wěn)定(圖4B)，此時其值為22.624 5，最終挑選的波長區(qū)間為[4，9，13，20](圖4C)。在ICO最終選中的聯(lián)合區(qū)間內(nèi)引入了局部搜索策略，進行寬度的自動優(yōu)化，最終挑選了805個特征波長。

圖4 重金屬Cr進行ICO變量選擇(區(qū)間數(shù)20)Fig.4 ICO variable selection for heavy metal Cr(interval number 20)A.sampling weight change graph of each interval during the iteration process(迭代過程中各區(qū)間的采樣權(quán)重變化圖),B.RMSECV value of the sub-model extracted in each iteration of the ICO(ICO 每次迭代中提取的子模型的RMSECV值),C.ICO algorithm for wavelength range selection(ICO算法進行波長區(qū)間選擇)

圖5 CARS運算提取變量原理圖Fig.5 Schematic diagram of CARS operation extraction variablesA.variation trend of the wavelength variable(波長變量的變化趨勢),B.variation trend of the RMSEcv(RMSEcv的變化趨勢),C.trend of wavelength regression coefficient(波長回歸系數(shù)的變化趨勢)

3.2 ICO-CARS算法選擇特征波長

將ICO算法挑選出的805個特征波長進一步使用CARS剔除區(qū)間波長中的無關(guān)變量。由于CARS方法在迭代過程引入蒙特卡羅采樣(MCS)和自適應(yīng)加權(quán)采樣(ARS)兩個隨機因素，造成每次運行挑選的波長結(jié)果不盡相同，為增加波長選擇的穩(wěn)定性和可靠性，對CARS進行100次重復(fù)計算，得到其RMSECV平均值為8.642 0，標準偏差(STD)為0.686 3。圖5為基于CARS算法的波長變量篩選過程圖，其MCS抽樣運行次數(shù)N為50，PLS主成分數(shù)為10，采用5折交叉驗證。

由CARS算法篩選有效波長變量的變化趨勢可看到變量數(shù)呈指數(shù)函數(shù)下降(圖5A)，選擇變量的個數(shù)也從急劇減少到緩慢遞減，最終趨于穩(wěn)定。波長選擇過程采用5折交叉驗證得到的RMSECV變化趨勢，若RMSECV值減小，表明剔除了無關(guān)信息變量；若RMSECV值增大則表明剔除了有效信息變量。由圖5B可見，在1～20次的采樣過程中，RMSECV呈遞減趨勢，第20次時的RMSECV值最低，此后開始增加，表明采樣20次時有效剔除了ICO提取光譜區(qū)間中的無關(guān)信息。圖5C表示各波長回歸系數(shù)隨著采樣次數(shù)增加而變化，圖中“*”為RMSECV值最低點。由圖可見，當運行20次時，RMSECV值最低，此時保留的變量為79個。

3.3 ICO-CARS-SPA算法選擇特征波長

經(jīng)過ICO-CARS變量篩選后，再通過SPA算法進行波長精簡(圖6)。通過最小誤差均方根值(RMSE)來確定最終所選特征波長個數(shù)，RMSE值越小，表明模型穩(wěn)定性越好、精度越高。

圖6 ICO-CARS-SPA變量選擇Fig.6 ICO-CARS-SPA variable selectionA.RMSE of SPA model varies the number of variables(SPA模型RMSE隨變量個數(shù)的變化),B.optimal characteristic wavelength selected by SPA(SPA模型選擇的最優(yōu)特征波長)

結(jié)果顯示，RMSE隨著波長個數(shù)的增加呈逐漸減小的趨勢，當波長個數(shù)大于33時，RMSE值變化不再顯著，此時RMSE值為15.439 9(圖6A)。說明最終優(yōu)選的33個特征波長保留了更多有效信息，作為最優(yōu)的波長點個數(shù)，選取的波長點在實驗光譜中的索引如圖6B所示，所選波段占原始光譜信息的0.81%。

采用ICO-CARS-SPA最終優(yōu)選出的33個敏感波長為輸入變量，利用PLS對重金屬元素Cr含量建模，并按照相同特征提取步驟對Cu、Zn、As、Pb進行特征優(yōu)選和PLS建模。結(jié)果表明，采用多特征串聯(lián)策略ICO-CARS-SPA-PLS建模后5種重金屬的R2、RMSE和MRB均達到了滿意的效果(表3)。

3.4 不同特征波長選擇方法比較

本文建立的ICO-CARS-SPA方法與ICO、最小絕對值收斂和選擇算法(LASSO)、CARS、ICO-CARS和ICO-SPA等[18-19]其他波長選擇算法的對比結(jié)果見表3。其中LASSO是一種系數(shù)壓縮變量選擇方法，即通過引入懲罰項使影響較小或者無影響的自變量系數(shù)趨近于零，進而可實現(xiàn)只保留與響應(yīng)變量最相關(guān)的解釋變量，將大量無關(guān)變量剔除，但高維數(shù)據(jù)中常會面臨統(tǒng)計學(xué)稱為維數(shù)災(zāi)禍的NP難題，因此不可避免地造成LASSO特征波長選擇個數(shù)將小于樣本個數(shù)，無法保證數(shù)據(jù)集中所有重要信息被完全保留。由于土壤中重金屬的X射線熒光光譜數(shù)據(jù)樣本數(shù)據(jù)集的光譜波長數(shù)遠大于樣本個數(shù)，LASSO算法不適合直接應(yīng)用于波長選擇問題。

表3 不同波長選擇算法PLS建模性能比較Table 3 Comparison of PLS modeling performance of different wavelength selection

由表3可知，盡管ICO算法可將局部搜索策略應(yīng)用于最終挑選出的最優(yōu)聯(lián)合區(qū)間進行寬度的自動優(yōu)化，但區(qū)間挑選算法選擇的波長數(shù)量一般較多，不利于快速預(yù)測，且忽視了聯(lián)合區(qū)間內(nèi)部光譜間存在的相關(guān)性和冗余問題，因此可對優(yōu)選出的聯(lián)合區(qū)間進行進一步波長篩選。而ICO-CARS建模效果優(yōu)于ICO-SPA是因為對于X射線熒光光譜而言，有效變量間的投影距離并不一定最大，在初選算法波長數(shù)量較多的情況下，SPA篩選出的變量子集中可能包含一些無關(guān)信息甚至是干擾變量。

通過對ICO、LASSO、CARS、ICO-CARS、ICO-SPA和ICO-CARS-SPA波長優(yōu)選對比分析可知，采用波長優(yōu)選算法建模后5種重金屬的R2、RMSE和MRB相對于單一PLS建模大多有不同程度地提升，其中以ICO-CARS-SPA-PLS提升最明顯。表明ICO-CARS-SPA-PLS的多特征串聯(lián)提升策略從整體上比其他模型優(yōu)化方法建模效果更好，更能兼顧到其他元素，能夠增強上一輪變量選擇的過程與下一輪變量選擇的關(guān)聯(lián)性，根據(jù)每輪變量選擇的好壞進行動態(tài)調(diào)整，在一定程度上解決了變量區(qū)間選擇的“筑巢效應(yīng)”。

3.5 模型廣適性比較研究

為驗證ICO-CARS-SPA變量篩選模型具有較好的廣適性，加入100次隨機分組進行實驗，5種重金屬元素的評價參數(shù)及標準偏差結(jié)果見表4。

表4 5種重金屬100次隨機分組ICO-CARS-SPA-PLS定量模型預(yù)測結(jié)果Table 4 Prediction results of ICO-CARS-SPA-PLS quantitative model of 100 random groups of five heavy metals

由表3～4可見，由于Zn元素在土壤本底中的濃度較大，導(dǎo)致在制備Zn樣本時，其配制濃度跨度較大，建模結(jié)果顯示驗證集上Zn的決定系數(shù)提高，但RMSEp效果較差，這可能是由于出現(xiàn)了過擬合現(xiàn)象。因此在后續(xù)實驗中通過擴大樣本數(shù)量，對樣本中的重金屬濃度劃分更細致，對Zn元素做進一步改進，可提高其預(yù)測效果。由表4可知，ICO-CARS-SPA-PLS模型具有較好的廣適性，能夠?qū)崿F(xiàn)對樣本中5種重金屬的定量檢測。

4 結(jié) 論

通過多特征串聯(lián)策略的ICO-CARS-SPA算法，建立了X射線熒光光譜重金屬含量的定量檢測模型。實驗制備了91個土壤樣本，剔除4個問題樣本后，選取了87個樣本，并對比分析了幾種常見的光譜特征變量選擇方法及其組合。結(jié)果表明，以ICO-CARS-SPA算法提取的特征波長建立的PLS重金屬含量預(yù)測模型精度高、誤差小、模型推廣性能好?；赬射線熒光光譜結(jié)合多特征串聯(lián)策略算法ICO-CARS-SPA是有效的特征提取方法，降低了模型復(fù)雜度，有利于推動X射線熒光光譜技術(shù)在農(nóng)業(yè)信息領(lǐng)域上的進一步應(yīng)用。