基于X射線熒光光譜與多特征串聯策略的土壤重金屬含量預測

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

(1.三峽大學 計算機與信息學院，湖北 宜昌 443002；2.三峽大學 湖北省農田環境監測工程技術研究中心，湖北 宜昌 443002；3.農產品質量安全追溯技術及應用國家工程實驗室，北京 100097)

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

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

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

1 原理與算法

1.1 試驗區間組合優化算法

區間組合優化(ICO)是基于模型集群分析策略(MPA)框架下的化學建模算法[8]，MPA強調從數據集中統計分析時，應通過隨機采樣，分析大量隨機產生的子集模型信息，從不同角度考察數據集的內在性質，獲得數據集的內在結構。ICO算法流程如圖1所示。

1.2 競爭性自適應重加權采樣算法

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

1.3 連續投影算法

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

1.4 多特征串聯算法

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

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

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

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

2 實驗部分

2.1 材 料

在周邊1 km范圍內無污染源的農田采集制樣土壤，為防止實驗器具影響和干擾樣本重金屬的濃度，樣本制作過程采用陶瓷用具。將土壤烘干研磨過0.45 mm孔篩后，參考《土壤環境質量標準》[17]中對Ⅰ、Ⅱ、Ⅲ類土壤中各元素含量的規定并結合農田土壤實際情況制作實驗樣本。每個樣本制作30 g土壤，共91個樣本。經初步測量，所采集農田土壤中As、Cu、Cr、Pb、Zn的含量依次為3.0、15、30、18、39 mg/kg。綜合國家土壤等級劃分標準和所采集土壤的具體情況，設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。按土壤需求濃度，將重金屬標準溶液用丙酮稀釋后混入土壤中，混合均勻后置于通風櫥中自然揮發，待土壤完全風干后，再研磨混勻，按指定標號放入樣品盒中，密封保存。

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

2.2 光譜數據采集

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

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

2.3 建模方法與評價

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

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

3 結果與討論

3.1 ICO算法選擇特征波長子區間

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

表2 不同子區間數PLS建模對比Table 2 Comparison of PLS modeling for different subinterval numbers

(續表2)

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

圖4 重金屬Cr進行ICO變量選擇(區間數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(迭代過程中各區間的采樣權重變化圖),B.RMSECV value of the sub-model extracted in each iteration of the ICO(ICO 每次迭代中提取的子模型的RMSECV值),C.ICO algorithm for wavelength range selection(ICO算法進行波長區間選擇)

圖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(波長回歸系數的變化趨勢)

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

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

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

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

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

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

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

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

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

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

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

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

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

3.5 模型廣適性比較研究

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

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

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

4 結 論

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