夏 可，張世文，沈 強，楊邵文，夏沙沙，胡青青

(1. 安徽理工大學 測繪學院，安徽 淮南 232001； 2． 安徽理工大學 地球與環(huán)境學院，安徽 淮南 232001)

1 引 言

中國作為世界人口大國之一，同時也是一個能源消耗大國。據(jù)相關資料顯示，目前我國礦山約有103 187個，大中型礦山9 399個，小型礦山和個體礦93 788個。我國礦產(chǎn)資源豐富，許多礦山已有近百年的開采歷史，為我國的經(jīng)濟建設和地區(qū)發(fā)展做出了極大的貢獻，與此同時，隨著礦產(chǎn)資源的開采，對礦山周邊的環(huán)境造成了嚴重的污染和破壞。習近平總書記在黨的十九大報告中指出，堅持人與自然和諧共生，樹立金山銀山就是綠水青山的理念，因此礦山生態(tài)修復以及土地復墾顯得愈發(fā)重要。截至2017年底，我國新增礦山恢復面積約4.43萬hm2，治理下的礦山約10 032個[1]。

土地復墾與生態(tài)修復是治理礦山的有力手段，礦區(qū)土地與環(huán)境監(jiān)測則是土地復墾與生態(tài)修復的前提和基本要求。在礦山開采過程中造成的礦業(yè)廢棄地污染以重金屬為主，重金屬污染對人體造成的危害巨大，因此土地復墾中對重金屬的檢測尤為重要。重金屬檢測可以為土地復墾提供方案選擇、復墾效果評價，實現(xiàn)長期動態(tài)監(jiān)測。傳統(tǒng)的重金屬檢測以實驗室測定為主，例如原子熒光法(AFS)、電感耦合等離子體法(ICP)、電感耦合等離子質譜法(ICP-MS)、原子吸收法(AAS)等，其效率低、成本高，無法實現(xiàn)大面積檢測，且長期動態(tài)監(jiān)測比較困難。近年來，高光譜遙感技術以其獨特的優(yōu)勢，在很多領域研究中得到了良好的發(fā)展與應用，這也為土壤重金屬檢測提供了一個全新的途徑。國內外眾多專家與學者利用高光譜遙感對土壤中各成分都進行過大量的研究，王維等發(fā)現(xiàn)一階微分高光譜反演模型對土壤重金屬Cu具有較好的估算潛力[2]；Srivastava等利用偏最小二乘回歸模型對印度旁遮普印度恒河平原的土壤有機碳含量進行了快速的預測[3]；Babaeian等利用不同光譜分辨率下的光譜特征提取土壤水力參數(shù)，并用HYPRES和Rosetta pedotransfer函數(shù)進行評價，結果表明估計的性能取決于水力參數(shù)的類型以及輸入信號的頻譜分辨率[4]；沈強等通過一階微分、倒數(shù)對數(shù)、連續(xù)統(tǒng)去除法結合多元逐步回歸模型，發(fā)現(xiàn)礦業(yè)廢棄地重構土壤中As的最佳反演模型為連續(xù)統(tǒng)去除法逐步回歸模型[5]；李晉華等通過多光譜數(shù)據(jù)進行多元散射校正預處理，并結合偏最小二乘模型對玉米成分進行預測，取得了較好的效果[6]；李旭青等利用光譜數(shù)據(jù)對水稻冠層氮素含量建立反演模型，結果表明隨機森林算法可有效地對水稻冠層氮素含量進行解釋、且所需樣本少、不易擬合[7]。

前人在土壤成分含量估測方面做了大量的研究并取得了較好的結果，但多集中在單一預測模型，單一預測模型往往不能全面反映事物的信息，信息的缺失又將會導致預測的偏差[8]，因此為了減弱單一預測模型中隨機因素的影響，提高模型預測精度，本文以四川省古藺縣石屏硫廠廢棄地復墾土壤為研究對象，基于Matlab和Python編程平臺以多種光譜數(shù)據(jù)預處理結合線性模型與非線性模型進行對比，在此基礎上尋求較優(yōu)單一模型，探討較優(yōu)模型組合方式在土壤重金屬含量高光譜估算中應用的可行性，以期為基于高光譜遙感技術下的礦業(yè)廢棄地復墾土壤重金屬快速檢測提供理論和方法支持。

2 實 驗

2.1 研究區(qū)概況

研究區(qū)位于四川省古藺縣石屏鄉(xiāng)(東經(jīng)105°59′54″～106°01′57″，北緯28°01′3″～28°02′51″之間)，地處四川盆地與貴州高原的過渡地帶，具有四川盆地氣候和高原氣候特征，四季分明、雨熱同期，年降水量偏少，溫度適中，光照條件較充足，年平均氣溫在13.8～18.6 ℃之間。研究區(qū)總面積約為4 451.89畝，從上世紀50年代始建，主要從事硫磺開采和冶煉，隨著時間的推移，礦區(qū)廢棄物堆積如山，對周邊土壤和環(huán)境造成嚴重破壞。研究區(qū)于2013年進行復墾，復墾方向主要為林地、耕地和草地，復墾措施主要為覆土、平整和土壤改良，并于2014年底完成復墾。

2.2 數(shù)據(jù)獲取

數(shù)據(jù)主要來源于實驗室測定和光譜測量。在研究區(qū)內采用網(wǎng)格隨機布點，共采集65個土壤樣本，土壤采樣深度為0～20 cm，主要涉及的復墾地塊類型有礦渣堆、旱地、林地、水田。將采集的土樣經(jīng)自然風干、研磨、過100目篩，各取100 g，分成兩份，一份用于實驗室測定土壤重金屬含量，一份用于光譜測量。實驗室土壤重金屬的測定采用電感耦合等離子法，光譜數(shù)據(jù)利用美國ASD(Analytical Spectral Devices)公司生產(chǎn)的Fieldspec 4便攜式地物光譜儀測定。在進行光譜測量時，選擇光線較暗的實驗室，提前預熱機器10 min，調整好探頭與光源位置，將樣本放在培養(yǎng)皿中，下面用黑布平鋪，實驗之前進行白板校正，實驗過程中每10個樣本進行一次白板校正。每個樣本測定10次，取其平均值作為最終的光譜數(shù)據(jù)，利用RS3軟件進行光譜數(shù)據(jù)收集，ViewSpecPro軟件進行數(shù)據(jù)后處理。利用K-S(Kennard-Stone)算法劃分建模樣本和驗證樣本，K-S算法能夠根據(jù)光譜變量間的歐式距離，在特征空間中均勻地選取光譜差異較大的樣本[9]，較傳統(tǒng)隨機選樣或根據(jù)距離和含量選樣更為合理。

2.3 數(shù)據(jù)預處理

光譜測量過程中由于儀器的原因、人為因素的影響以及外界環(huán)境的干擾都會對光譜數(shù)據(jù)造成影響。而這些噪聲對特征波段的選取和模型的建立都會造成影響，為此我們首先對光譜數(shù)據(jù)進行S-G平滑處理以減弱噪聲[10]，其次再對數(shù)據(jù)進行一階微分[11](First order differential，F(xiàn)DR,公式(1))、離差標準化[12](Deviation standardization, DS,公式(2))、多元散射校正[13-14](Multiplicative scatter correction,MSC,公式(3)～(5))和連續(xù)統(tǒng)去除法[15-16](Continuum removal,CR)等預處理變換。

(1)

(2)

其中，λi為第i個光譜波長，R(λi)為對應光譜波長的反射率，Δλ為λi+1與λi之間的差值，R(λmax)、R(λmin)分別為某一樣本的最大反射率和最小反射率。

(3)

(4)

(5)

2.4 建模方法

2.4.1 偏最小二乘

偏最小二乘(Partial least squares，PLS)的基本原理是通過最小誤差的平方和找尋一組數(shù)據(jù)的最佳匹配函數(shù)，它是傳統(tǒng)的多元線性回歸、典型相關分析和主成分分析的集合體[17-18]。較傳統(tǒng)多元線性回歸而言，偏最小二乘能夠允許自變量直接存在多重相關性，可以在樣點個數(shù)低于變量個數(shù)的情況下實行建模。建模原理為：建立m×n的光譜矩陣X，n×l的重金屬含量檢測矩陣Y，其中m為光譜波段數(shù)，n為樣品個數(shù)，l為重金屬種類。

將X、Y進行分解，公式如下：

X=TPT+E，

(6)

Y=UQT+F，

(7)

其中，U、T為得分矩陣，P、Q為載荷，E、F為PLSR模擬時的殘差矩陣。對U、T做線性回歸，B為關聯(lián)系數(shù)矩陣，有U=TB，重金屬含量預測公式為：

Y預測=T計算BQ=X測量PTBQ.

(8)

2.4.2 人工神經(jīng)網(wǎng)絡

人工神經(jīng)網(wǎng)絡(Artificial neural network，ANN)是一種仿生物神經(jīng)網(wǎng)絡技術，由大量神經(jīng)元連接而成，主要分為輸入層、隱藏層和輸出層。人工神經(jīng)網(wǎng)絡具有自適應和自主學習的能力，因而在處理大量隨機性數(shù)據(jù)和非線性數(shù)據(jù)方面具有明顯優(yōu)勢[19]。

2.4.3 隨機森林

隨機森林(Random forest，RF)算法最早由Leo Breiman和Adele Cutler所提出，能夠處理高維特征且不易產(chǎn)生過擬合，對于大數(shù)據(jù)而言，模型訓練速度較快，對數(shù)據(jù)的適用能力較強[20-21]?；驹頌椋涸O訓練樣本個數(shù)為M，特征數(shù)目為N，從M個訓練集中以bootstrap取樣方式取樣M次，對未抽取樣本進行預測和評估誤差，針對每一節(jié)點隨機選取n個特征，根據(jù)這n個特征選擇最佳分裂方式，其次在每個訓練子集上構建決策樹，最后根據(jù)每棵決策樹的輸出取平均或投票作為最終結果。

2.5 熵值權系數(shù)組合模型

組合模型是將多個模型綜合在一起，賦予不同的權重進而得到預測模型。組合模型預測的核心就是如何確定模型的權重，由于預測模型中預測值與實測值之間的殘差具有不確定因素，而熵值法在組合模型中有較好的適用性[22]。其計算公式為：

(9)

(10)

(11)

di=1-hi,

(12)

(13)

3 結果與討論

3.1 數(shù)據(jù)預處理結果

原始光譜數(shù)據(jù)受到噪聲、樣本背景和其他無關成分等干擾，選用合適的預處理方法能夠消除噪聲，提升模型預測能力。因此在建立重金屬定量反演模型中，光譜數(shù)據(jù)預處理十分重要。對原始光譜進行S-G平滑處理，光譜波段范圍為350～2 500 nm，由于光譜在采樣過程中光譜波段兩端產(chǎn)生較大的噪聲，為此，剔除了350～449 nm和2 451～2 500 nm范圍內的波段數(shù)據(jù)，共采集2 000組波段數(shù)據(jù)，如圖1所示。圖2為不同預處理結果，依次為FDR、DS、MSC和CR。

圖2 不同預處理光譜反射曲線

3.2 土壤重金屬含量分析

為了獲悉pH和不同重金屬含量之間的差異，利用SPSS軟件做描述性統(tǒng)計，如表1所示。

由表1可知，研究區(qū)的pH范圍在2.90～8.28之間，均值為6.38，土壤呈酸性居多，較少部分為弱堿性，主要是由于該研究區(qū)長期進行采礦活動所導致。從含量統(tǒng)計分析表中可以看出，5種重金屬之間的含量差異比較明顯，其中Cd和Hg的含量偏低，Cr和Ni的含量較高；從標準差來看，Cr的離散程度最大，Hg的離散程度最小；從變異系數(shù)來看，5種重金屬變異系數(shù)均在30%以上，其中Cd的變異系數(shù)最大，達到114.44%，這可能是由于先前的采礦活動對周邊環(huán)境影響較大，導致不確定因素增加，引起Cd分布發(fā)生異變；土壤中過量的Cd會抑制植物的生長，亦能通過食物鏈嚴重危害人體的健康，因此Cd異變應引起注意。

3.3 相關性分析

重金屬在不同的波段對光譜的吸收強度也不同，因此不同波段所顯示的相關性也不同。將經(jīng)預處理后的光譜數(shù)據(jù)與土壤重金屬含量做相關性分析，找尋顯著相關的波段，以期為后期建模打下基礎。相關性大小用皮爾森系數(shù)R表示，圖3分別為一階微分(a)、離差標準化(b)、多元散射校正(c)和連續(xù)統(tǒng)去除法(d)的相關系數(shù)曲線圖。

圖3 土壤重金屬含量與不同變換形式下光譜反射率的相關系數(shù)

從圖3中可以看出，經(jīng)不同的預處理，不同重金屬相關性幅度變化較大，但都有不同程度的提升。經(jīng)一階微分變換，Cd、Cr、Ni、As、Hg相關性較原始光譜相關性均有所提高，對Cd和Cr相關性提高最為明顯，在1 910，2 206 nm處絕對值分別達到0.475和0.532，Hg相關性在1 910 nm處達到極大值，為-0.706；5種重金屬與光譜曲線的相關性系數(shù)基本都低于0.6，但大部分波段達到了水平為0.05以上的顯著相關，部分波段達到水平為0.01以上的極顯著相關。經(jīng)離差標準化變換，5種重金屬相關系數(shù)沒有一階微分提升明顯，但多種重金屬相關性整體提高，以Cr和Ni最為明顯，Cr波段相關性絕對值整體位于0.4以上，Ni為0.3以上；對于Cd、As和Hg而言效果不為明顯，其中Cd效果最不理想。經(jīng)連續(xù)統(tǒng)去除法變換后，Cd、Cr、Ni、As和Hg在1 394，2 208，2 212，1 418，1 415 nm處達到極值，分別為-0.398，-0.484，-0.481，-0.57，-0.513，相關性雖說有所提升，但并不是最佳。經(jīng)多元散射校正變換后，Ni和As相關性極好，相關系數(shù)分別在1 929 nm和568 nm處達到極值，為-0.679和0.715，整體波段也表現(xiàn)較高的相關性；Hg相關性在1 600 nm處達到極大值0.636，未有一階微分達到的極大值高，但整體波段達到了很好的相關性，有一半波段相關性絕對值位于0.5以上；Cd在1 279 nm處達到極大值0.498，較比其他3種預處理，多元散射校正效果達到最佳；Cr在2 316 nm處達到極值-0.423，提升效果不為明顯?？傮w來說，經(jīng)多元散射校正變換后，多種重金屬達到極高的相關性，表明多元散射校正能夠有效去除噪聲以及由散射帶來的基線漂移等干擾，增強與樣品成品相關的光譜信息[23]。

3.4 單一模型建立與分析

利用K-S算法將樣本劃分為建模集和驗證集(建模樣本50個，驗證樣本15個)，對樣本進行交叉驗證。通過對光譜數(shù)據(jù)不同的預處理做相關性分析，選擇相關性較大和顯著性波段分別用于PLS、ANN和RF模型的建模。采用決定系數(shù)(Determination coefficient,R2)、均方根誤差(Root mean square error,RMSE)和相對分析誤差(Relative percent deviation,RPD)指標對模型進行精度評定。其中決定系數(shù)作為數(shù)據(jù)擬合程度的一個評定，越接近于1，效果越好；均方根誤差反映預測數(shù)據(jù)的精密度，用來衡量預測值與真值之間的偏差，RMSE越小表明預測精度越高；相對分析誤差指預測偏差，它是標準差與均方根誤差的比值，RPD的高低反映模型的預測能力，當RPD>2.0時，模型的預測能力很好，當2.0>RPD>1.5時，模型的預測能力一般，當RPD<1.5時，表明建模失敗[24]。表2～表4依次為PLR、ANN和RF建模結果。

表2 偏最小二乘模型(PLS)

表3 人工神經(jīng)網(wǎng)絡模型(ANN)

表4 隨機森林模型(RF)

從表2～表4中可以看出，不同的預處理對3種建模方法來說差異較大，從驗證集的R2、RMSE和RPD來看，多元散射校正結合隨機森林建模，R2與RPD普遍達到了0.80和2.0以上。偏最小二乘模型預測效果一般，其中As經(jīng)多元散射校正變換后效果最好，R2與RPD分別為0.82和2.14。離差標準化結合神經(jīng)網(wǎng)絡對于含量相對較高的Cr、Ni預測效果較好，R2和RPD最高達到0.83和2.8，表明在進行重金屬高光譜含量估算中需要考慮重金屬含量對建模反演效果的影響。較比前兩種建模方式，隨機森林表現(xiàn)了優(yōu)異的估測能力，經(jīng)多元散射校正變換多種重金屬的R2和RPD有明顯提升，其中As效果最好，R2、RMSE和RPD分別達到0.89，2.00，2.77，其次是Hg、Ni和Cd，R2均達到0.80以上，Cr較比其他3種預處理變換效果無明顯差異。

根據(jù)以上分析可知，結合MSC和RF的優(yōu)勢，相比其他預處理，MSC-RF模型整體來說要略勝一籌，其次DS-ANN對重金屬含量相對偏高的元素也表現(xiàn)出了較好的預測能力?？傮w來說，采用多元散射校正結合隨機森林建立的重金屬反演模型效果最好。

3.5 組合模型與驗證結果

單一模型自身存在一定的局限性，而組合模型能夠“取長補短”，發(fā)揮多種模型的優(yōu)勢[25]；組合模型的關鍵問題在于單一預測模型的選取以及相應權重系數(shù)的確定，本文針對不同重金屬，分別選取兩種較優(yōu)單一模型，利用熵值法確定模型權重系數(shù)，結果如表5和圖4所示。

表5 組合模型參數(shù)

Notes：Yrepresentative predictive value，Capital letters represent models，The subscript represents the preprocessing transformation.

圖4 土壤重金屬含量較優(yōu)模型與組合模型預測散點圖比較

熵值法能夠根據(jù)單種預測模型預測誤差序列的變異程度來確定組合模型的權重系數(shù)而且計算簡單[26]，提供了一種客觀賦權的方法。通過利用熵值法進行模型組合，從表5中土壤重金屬含量的預測結果可以看出，驗證集精度相比傳統(tǒng)單一模型驗證集精度有了顯著的提高，其中As的R2值達到最高，相比其最優(yōu)模型R2由0.89提高至0.91，RMSE由2.00降低至1.85；對于Cd、Cr和As線性模型和非線性的組合表現(xiàn)出了優(yōu)異的估測能力，尤其對于Cd，R2和RPD分別由最優(yōu)單一預測模型的0.83和1.76提升至0.85和2.39，RMSE由0.18減少至0.16；表明利用熵值法確定的組合模型估算土壤重金屬含量是可行的。圖4為重金屬含量實測值與預測值散點圖，散點越接近于1∶1對角線，效果越好，從圖中可以看出，組合模型散點較多集中在對角線附近，且散點趨勢線與對角線之間的角度差很小，表明結合多種模型能夠對重金屬含量的預測產(chǎn)生一個很好的效果。由于土壤為礦物質、有機質和水分等物質組成的復雜有機整體，各成分之間相互影響、相互作用，土壤光譜易受區(qū)域性和地域性影響，因此不同區(qū)域內的土壤光譜有所差異。本文以四川古藺礦業(yè)廢棄地復墾土壤重金屬為研究對象，所建立的重金屬含量反演模型在其他區(qū)域的適用性還有待于進一步研究。

4 結 論

針對高光譜反射率反演土壤重金屬含量問題，本文采取4種預處理變換并結合PLS、ANN和RF 3種建模方法建立礦業(yè)廢棄地復墾土壤反射率與土壤重金屬的單一定量模型，并采用熵值法進行較優(yōu)模型組合，得到如下主要結論：

(1)原始土壤光譜反射率和重金屬之間的相關性較弱，經(jīng)不同預處理變換后，部分波段相關性有了明顯的提升；其中MSC總體效果較好，對As和Ni最為明顯，相關性系數(shù)極值達到0.715和-0.679，且整體波段都取得了較好的相關性；其次FDR對Hg較為明顯，相關性極值達到0.706；DS對重金屬含量相對偏高的Cr的相關性起到了很好的提升效果，絕大部分波段相關性系數(shù)位于0.4以上，最大極值達到-0.530。

(2)預處理變換和模型的選取不同對預測結果影響較大，利用3種精度指標對模型進行評價，其中MSC-RF反演效果最好，驗證集As、Hg、Ni和Cd的R2、RMSE和RPD分別達到：0.89，2.00，2.77；0.84，0.03，2.52；0.83，4.78，2.72；0.83，0.18，1.76；DS-ANN對Cr效果較好，3種精度指標分別為：0.83，15.66，2.19；利用MSC和DS結合非線性模型效果要優(yōu)于線性模型，由于這兩種變換能夠顯著提升整體波段的相關性，為RF和ANN模型提供了大量的輸入?yún)?shù)；結果表明，針對不同的建模應選取適宜的預處理變換方法。

(3)從單一模型和組合模型預測結果來看，單一模型雖具有運算速度快、操作簡便等優(yōu)勢，但組合模型能綜合利用多種模型信息，減弱單一模型中隨機因素帶來的影響，應提高組合模型在土壤光譜預測方面的應用性。熵值法根據(jù)預測模型的變異程度確定權重，避免了權重確定的主觀性，客觀地反映了單一預測模型在組合模型中的重要性，利用熵值法進行兩種較優(yōu)模型組合，較傳統(tǒng)單一模型，精度指標均有所提高，表明采用熵值組合模型對土壤重金屬含量進行預測效果更好。