礦山土壤重金屬含量估算研究

饒振興，馬超振，邵鵬遠，翟會君，顧展飛

(1.河南省地質礦產勘查開發局 第四地質勘查院，河南 鄭州 450001； 2.河南省地質礦產勘查開發局 第二地質環境調查院，河南 鄭州 450001； 3.鄭州航空工業管理學院，河南 鄭州 450001)

土壤資源關于農業生產安全，并直接影響人類的身體健康。近年來，因煤礦開采而導致的土壤污染問題逐步受到人們的關注。煤礦開采雖能促進社會和經濟的發展[1]，但不合理的開發將會造成一系列的環境問題，包括水資源污染、地面塌陷、溫室氣體排放污染以及噪聲污染等[2-5]，嚴重影響生態環境。其中，煤矸石等廢料的長期堆放會導致土壤中的重金屬異常富集[6]，再經由農作物的吸收間接進入人體內，危及人類健康[7]。因此，高效地估算土壤中的重金屬含量，尤其是在礦山環境條件下，意義重大。

傳統的估算方法一般是通過實驗室的精密分析，但成本高昂且操作效率低[8-9]。隨著科技水平的發展，遙感技術在高光譜領域發展迅速，為土壤污染物的高效估算提供了可能[10-12]，但基于高光譜的重金屬含量反演推斷技術仍大多以實驗室的光譜數據分析為基礎，一定程度上限制了其使用效率。為此，已有不少專家針對重金屬含量和野外土壤高光譜的關系，嘗試建立直接的映射函數，并取得了部分成果[13-14]，如陶超等[15]采用支持向量機(SVM)作為分類器，將一定數量的室內樣本數據用于分類器訓練，同時對野外樣本數據進行分類測試，通過光譜變換后進行了重金屬污染定性分類模型遷移試驗，結果表明不同試驗區域的砷(As)、鉛(Pb)和鋅(Zn)等重金屬污染狀況可直接利用室外光譜采樣進行研究。鄒濱等[16]依托湖南衡陽市某礦區，針對土壤鎘(Cd)含量測定，開展了實驗室與野外原位DS(direct standardization，直接校正法)光譜聯合反演模型的有效性驗證。試驗結果表明，將野外原位光譜DS轉換校正結合融合實驗室光譜先驗知識的模型精度最高，其決定系數R2達到了0.72，對野外原位光譜模型的可靠性提升具有重要作用。張霞等[17]針對野外高光譜實時監測土壤重金屬含量過程中存在的光譜數據冗余問題，基于雄安新區的野外農田土壤樣本，分別從反演機理和波段選擇算法2個角度出發，提出了相應的反演方法，實現了土壤樣本中鉛含量估算精度的有效提升。綜上所述，在野外土壤高光譜與重金屬含量估算的模型關系方面，較多研究集中于可行性或有效性的驗證以及存在問題的完善方面，缺乏基于偏最小二乘模型的高光譜估算模型研究及其具體應用方法的介紹。因此，本文以某城市煤礦環境下的周邊土地為研究基礎，為將土壤中隱蔽性較強的重金屬污染物含量進行高效估算與表征，建立了偏最小二乘模型的高光譜估算模型，其中偏最小二乘模型的建立基礎為土壤反射光譜數據與光譜變換指標，從而突破了上述基于實驗室數據進行重金屬含量估算的局限性。

1 礦區概況

1.1 基本情況

某城市北部礦區內土壤成土類型較多，土壤pH值在6.5～7.5。第四系松散沉積物在該城市北部礦區內具有普遍的分布，厚度達50 m左右，部分區域超過100 m。砂壤土和砂土為區域內主要的表層土質，具有良好的耕性。土壤中鉀、磷、氮等含量較為豐富，整體評價為肥沃。該礦區經過長期的煤礦開采，對周邊土地環境影響巨大，原有土地利用結構遭到破壞，可耕地面積減少，同時大量的采空區形成改變了原有地貌，使農田出現下沉。與此同時，該礦區的周邊土壤不可避免地遭到了不同程度的污染，土壤重金屬富集污染較嚴重，造成周邊居民身體健康的潛在隱患。

1.2 主要治理對策及原理

為更有針對性地采取土壤污染治理對策，實現可持續發展，需對土壤重金屬含量進行估算和評價。為此，提出并綜合應用了一種根據現場高光譜數據實現反演模型建立的方案，從而準確地估算出土壤中的重金屬含量。其基本原理為：土壤中的重金屬含量較低，其對于光譜的響應來說是非常微弱的，即無法通過光譜表面直接得出土壤重金屬的光譜反應。但重金屬與土壤中的黏土或鐵氧化物等較易聚集在一起，且光譜響應對此理化性質比較敏感，故重金屬含量可先與土壤中的鐵氧化物建立一定的函數關系，再通過上述理化性質建立高光譜模型，從而間接進行重金屬含量的估算[18-20]。同時，在研究方法方面，為解決高光譜波段數據量較大、不同波段間相關性強等問題，擯棄傳統的逐步回歸和主成分方法，采用了偏最小二乘方法，進一步完善了高光譜估算模型。

2 數據采集及預處理

2.1 土壤樣本制備

基于格網GIS技術，結合礦區實地情況，按1.5 km2范圍采集1個樣點的原則進行控制，根據布點方案采集土壤樣本36個，如圖1所示。樣本采集時，盡量選擇土壤粒度及含水量相近的區域，以避免土壤環境對光譜的影響。

圖1 礦區土壤采樣點示意Fig.1 Schematic diagram of soil sampling points in mining area

2.2 野外光譜測定

采用手持地物光譜儀進行本次野外光譜測定。測定時，選取太陽輻射強烈、晴朗無云的中午時段，風力不大于3級且空氣能見度在10 km以上。操作人員宜穿深色衣服，操作時保持與測量的樣本相距約500 mm，同時，光譜儀跟水平面的法線的夾角需控制在一定范圍內，該值一般為10°。采集時，每10 min定標1次。為便于后期取平均數，在每一個采樣點處，均需要進行10條光譜數據的測量。

2.3 土壤理化性質測試

采集的土壤樣本作通風陰干處理，經壓碎研磨后，過篩備用。在試驗過程中，擬通過電感耦合等離子體原子發射光譜法對土壤的理化性質進行檢測，主要測定土壤中Cr、Mn、Cu、Fe、Ni、Cd、Zn、As、Pb等重金屬的含量。該方法在測定時，先進行微波消解，冷卻后加超純水搖勻，再按提前確定好的參數檢測樣本里元素的發射強度，并基于校準函數計算出重金屬元素的含量。該方法在操作時，不僅快速精準，而且可同時對多種元素進行測定，另外還具有線性范圍寬的優勢，應用廣泛。

2.4 其他相關數據采集

在土壤的野外高光譜數據采集前，采用影像融合方法，將2.5 m全色波段數據與Spot-5衛星中10 m多光譜數據進行融合，由此獲得礦區2.5 m高清彩色影像。同時，由于人為耕種與施肥等因素對土壤的理化性質存在較大影響。因此，搜集了礦區不同采樣點的土地利用類型。另外，此次試驗還收集了該地區土壤類別空間分布數據，以更好地考慮重金屬元素臨界含量與容量在不同土壤類別的差異。

2.5 高光譜數據預處理

土壤光譜數據由于是野外采集的，故需要在模型擬合以及光譜指標變換前，對相關數據進行平滑降噪、剔除不合格光譜等預處理，并進行均值化。

3 土壤重金屬污染估算

3.1 高光譜估算模型建立方案

在模型建立時引入偏最小二乘回歸模型，相對于傳統的多元線性等回歸模型，其不僅吸收了主成分分析方法中提取主成分的優勢，在進行回歸建模時綜合了自變量信息，而且解決了變量中多重相關的問題。在完成建模后，為更好地檢驗實測值與預測值的關系，評估方法的準確性，提出了分析方法中常用的相關系數、均方根誤差等指標參數。

3.2 提取高光譜指標信息

在獲取土壤光譜后，其一般含有繁雜的各種信息，難以分清光譜波段的相關關系。為更高效地實現有效數據獲取，對土壤高光譜特征進行更準確的表征，基于如圖2所示的高光譜反射率趨勢，提取出4個關鍵的指標，具體涵蓋有去除包絡線、反射率倒數對數以及一階、二階微分指標。

圖2 土壤的反射率光譜走勢示意Fig.2 Schematic diagram of reflectance spectrum trend of soil

(1)去除包絡線。每個土樣的高光譜曲線整體發展趨勢類似，為有效體現土樣波譜的反射峰和吸收峰，使之歸一化到統一的背景值中(0～1)，放大微小的波譜變化，可進行去除包絡線處理，如圖3所示。

圖3 土壤反射率光譜去除包絡線曲線 Fig.3 Removal envelope curve of soil reflectance spectrum

(2)為增強光譜曲線在坡度上的響應，消除空氣中氣體對光譜的影響，并在一定程度上突出對可見光波段的光譜波段的差異性，可進行反射光譜的倒數對數變換，如圖4所示。

圖4 土壤反射率光譜倒數對數曲線Fig.4 Reciprocal logarithm curve of soil reflectance spectrum

(3)受儀器自身的影響，光譜儀在使用中所采集的光譜基線會隨著時間而出現漂移，不利于后期偏最小二乘法的建模。為此，可采取光譜微分處理，不僅可以去除基線漂移現象，還能提取到光譜曲線中較為隱蔽的信息。同時，更高的波譜分辨率可以在一階和二階微分處理后的反射光譜中實現，如圖5所示。

圖5 土壤反射率光譜一階微分和二階微分曲線Fig.5 Soil reflectance spectrum first and second order differential curves

3.3 對重金屬的相關性分析

根據高光譜估算理論原理，為滿足高光譜建模條件，需對土壤中各重金屬元素含量與Fe、Mn含量的相關性進行判定。對于重金屬元素的相關系數，一般來說，當其值大于0.5，同時呈現為顯著相關的，則表明兩者呈強顯著相關性，據此可以實現高光譜的建模分析。本次土壤中幾種重金屬元素含量的相關系數矩陣見表1，可知，Fe、Mn的相關系數為0.94，體現了兩者的共生關系。Ni、As、Cr、Pb、Zn等與Fe的相關系數均在0.5以上，說明均呈強顯著相關。而Cu與Fe、Mn的雖呈顯著相關，但相關系數未超過0.5；Cd與Fe既未呈顯著相關，相關系數也較小。因此，在上述重金屬元素中，不滿足高光譜建模要求的只有Cu與Cd元素。

表1 重金屬相關系數矩陣Tab.1 Correlation coefficient matrix of heavy metal

3.4 提取光譜特征波段

野外現場測定試驗中，光譜分辨率設定為1 nm，通過測定獲得的光譜波長區間在350～1 050 nm。高分辨率在提高測量數據精度的同時，也帶來了非常多的無效信息，導致后續建模困難。因此，可將光譜分辨率增大十倍，即從原先的1 nm變更至10 nm，達到簡化模型、提升工效的目的。在這之后，接著計算樣本重金屬含量與上述4個光譜指標的相關函數，并按照所得關系，優先出較為敏感的波段，以此建立相關性函數模型。其中可將波段作為高光譜建模的相關性標準如下：不論是正相關還是負相關，均需滿足0.05的顯著性水平，同時相關系數需大于0.5。

3.5 偏最小二乘回歸模型的建立

根據重金屬含量與不同光譜指標的相關性分析結果，結合特征波段和范圍，將36個土壤樣本分為驗證樣本與建模樣本兩部分，前者設定為10個，后者建模分析用取26個。在建模完成后，對于所得結果，分別采用決定系數(R2)和均方根誤差(RMSE)進行精度對比，見表2(考慮篇幅，僅列出元素Ni的建模結果)。其中，為防止出現比較離譜的數據結論，采用了決定系數的指標，一般其值大于0.3時，表示具有一定的參考意義。而對于均方根誤差來說，其值與精度呈正相關，與誤差為負相關。在表2所示Ni元素的回歸模型中，建模精度最高的是二階微分，在建模樣本中，決定系數與均方根誤差分別為0.869 3、16.980 0；在驗證樣本中，決定系數與均方根誤差分別為0.798 3、17.990 0，對于土壤中重金屬含量的估算要求基本是可以滿足的。而對于二階微分外的其余光譜指標，精度均不是很高，不適用于重金屬含量估算。

表2 重金屬元素Ni的建模結果Tab.2 Modeling results of Ni heavy metal elements

3.6 估算結果精度評價

為更加好地比較各光譜指標的估算精度，直觀展示建模樣本與驗證樣本的實測值、預測值，并對建模結果進行評價，制作了相關數據的散點圖。本文僅以元素Ni和Mn為例，如圖6所示。

圖6 元素Ni、Mn的數據散點Fig.6 Data scatter diagram of elements Ni and Mn

在圖6中，三角形為建模樣本，圓形為驗證樣本，斜線為y=x的趨勢線，模型擬合效果主要由散點與趨勢線的距離體現，距離越近，則表示精度越高。由圖6(a)可知，除少數樣本距離趨勢線較遠外，大部分的散點與趨勢線均較近。因此，Ni元素的二階微分回歸模型擬合效果較好。由圖6(b)可看出，其大部分散點亦分布在趨勢線周圍，但上下位置分別有部分建模樣本偏離趨勢線較遠，故Mn元素一階微分回歸模型的精度低于Ni。據此可同樣分析出其余重金屬元素的模型精度，本文不再贅述。

由模型精度分析可知，基于偏最小二乘法的應用模型能有效估算重金屬含量，但最優估算模型隨重金屬元素種類的不同而變化，對于重金屬元素Cr、Zn、As、Pb、Mn、Ni，其效果依次逐漸變好。其中，Ni元素對應估算模型最適合的光譜指標為二階微分，Mn元素對應為一階微分，Pb元素對應為反射率，As元素對應為反射率，Zn元素對應為二階微分，Cr元素對應為去除包絡線。由此可知，為提高建模精度，不同的重金屬元素應嘗試運用不同的光譜指標；另外，元素Ni、Mn、Zn最優模型對應的光譜指標均為微分變換值，說明在提高建模精度的措施方面，利用微分方法處理光譜反射率是一種有效的途徑。

4 結語

(1)土壤中的重金屬較易與黏土、鐵氧化物以及有機質等聚集，土壤的這些理化性質在高光譜數據中均有明顯的光譜響應，這為利用高光譜模型間接估算重金屬含量提供了基礎原理。

(2)土壤的性質與光譜的關系為非線性，因此應用非線性的偏最小二乘等方法，可極大提高預測的精確性；采取去除包絡線、反射光譜倒數對數、微分數據處理等方法，既可消除高光譜噪聲，也能提高建模精度。

(3)土壤重金屬估算的建模結果顯示，所述估算模型基本能滿足估算要求，同時也有良好的精度保證。另外，對光譜反射率進行微分處理是提高建模精度的有效方法。