融合光譜和空間特征的土壤重金屬含量極端隨機樹估算

于海洋 謝賽飛 郭靈輝 劉 鵬 張 平

(1.河南理工大學測繪與國土信息工程學院， 焦作 454003；2.河南理工大學自然資源部礦山時空信息與生態(tài)修復重點實驗室， 焦作 454003;3.河南省自然資源科學研究院河南省國土資源動態(tài)監(jiān)測重點實驗室, 鄭州 450053;4.河南省航空物探遙感中心遙感技術研究院, 鄭州 450053)

0 引言

礦業(yè)活動、冶金以及工業(yè)生產中產生的無機污染物通過大氣降塵和污水排放等途徑進入土壤，不斷積聚造成較為嚴重的土壤重金屬污染，這些污染物滲透進入土壤后移動性較差，殘留時間長，且不易被微生物降解，嚴重威脅生態(tài)安全。據統(tǒng)計，我國耕地土壤污染超標率高達19.4%，其中以重金屬污染最為嚴峻[1]。因此，高效、快速獲取土壤重金屬含量及空間分布，對于重金屬污染防治、農業(yè)生產和生態(tài)安全具有重要意義。

傳統(tǒng)土壤重金屬污染調查采用實地采集土壤樣品和化學分析的方法，需耗費大量的人力物力資源。遙感光譜反演方法為快速、高效的土壤重金屬污染信息獲取提供了可選方案，已有文獻對采用高光譜遙感技術監(jiān)測土壤重金屬污染進行了有益的嘗試[2]。主要針對光譜信息增強變換與特征選取、反演算法[3]等進行了分析，使用的光譜信息增強方法包括光譜微分(一階、二階等)[4]、連續(xù)統(tǒng)去除、高斯卷積平滑、多元散射校正[5]等，常用反演模型包括多元線性回歸(Multivariable linear regression，MLR)、偏最小二乘回歸[6-8]、支持向量機(Support vector machine，SVM)[9-10]、極限學習機[11]、人工神經網絡、隨機森林(Random forest，RF)[12]等。已有研究進行反演建模時僅考慮了土壤光譜特征，由于土壤光譜影響因素眾多，重金屬元素含量低，光譜信息弱，導致反演模型魯棒性差，泛化能力不強[13]。

土壤重金屬污染物主要來源于礦業(yè)活動、冶金以及工業(yè)生產中形成的大氣降塵、污水排灌等，潛在污染源長期存在，一般較為明確。因此，考慮構建污染源-匯空間特征量化污染物擴散與匯集空間影響因子，融合光譜特征建立土壤重金屬含量估算模型。近年來，極端隨機樹(Extremely randomized trees，ERT)集成學習方法在其他領域機器學習建模研究中表現優(yōu)異，具有魯棒性高、泛化能力強的特點[14]，本文將該算法引入土壤重金屬含量反演，以期進一步提升模型預測精度和泛化能力。

1 研究區(qū)概況

研究區(qū)位于河南省濟源市市區(qū)周邊的山前平原區(qū)(圖1)，分布于112.46°～112.68°E，35.01°～35.19°N之間，東、南分別以二廣高速、菏寶高速為界，西、北以山麓地形為界，面積263.5 km2。濟源市區(qū)周邊分布的部分金屬冶煉廠導致其耕地土壤形成了多種重金屬元素濃度異常分布區(qū)，其主要污染元素為鉛(Pb)、鉻(Cd)等，表層土壤鉛異常濃集中心位于鉛鋅冶煉廠附近，高值區(qū)Pb平均含量(質量比)為611 mg/kg，達到三級污染標準值的1.22倍。

2 實驗數據獲取與預處理

2.1 土壤光譜數據測定及預處理

野外土壤樣品采集時間為2019年10—11月，按照放射狀為主、環(huán)形補充的網狀方法布設采樣點(圖1)，采用GNSS RTK定位后采集表層土壤(0～20 cm)作為樣本，采集數量為249個。

圖1 研究區(qū)位置和土壤樣本分布Fig.1 Study area and distribution of soil samples

土壤光譜數據采用美國ASD公司的ASD Fieldspec3型光譜儀進行測量，光譜波長范圍為350～2 500 nm。在實驗室內將土樣干燥、研磨后過100目篩，使用高密度反射探頭測量土壤樣本光譜反射率，每個樣本連續(xù)測定10次，取其平均值作為最終的光譜反射率(圖2)。

圖2 土壤光譜曲線Fig.2 Soil spectral curves

因受外界測量條件和傳感器本身的影響，在土壤樣品光譜采集過程中存在不同程度的噪聲，故有必要對原始光譜進行降噪處理。采用Savizky-Golay(SG)卷積平滑法對光譜曲線進行平滑去噪。原始光譜反射率測量值與土壤重金屬含量之間的相關性較低，對濾波后的光譜反射率進行光譜變換增強光譜信息。研究選取的光譜變換方法包括：一階微分(First order differential reflectance，FD)、二階微分(Second order differential reflectance，SD)、多元散射校正(Multiplicative scatter correction，MSC)、標準正態(tài)變量變換(Standard normal variate，SNV)、連續(xù)統(tǒng)去除(Continuum removal，CR)、倒數一階微分(Reciprocal first derivative, RFD)、倒數二階微分(Reciprocal second derivative, RSD)。其中多元散射校正可以有效消除樣品散射的影響[15]；標準正態(tài)變量變換通過加權平均化消除固體的顆粒不均對光譜的影響；光譜微分可以有效消除基線和減弱背景干擾，分辨混合光譜，增強光譜特征[16]；連續(xù)統(tǒng)去除可以有效突出光譜曲線的吸收和反射特征[17]。

光譜特征中含有較多的冗余和共線性變量，在進行定量分析前，需篩選重金屬元素的特征波段，提高建模效率。將相關性分析中滿足P=0.01假設性檢驗的波段集合作為特征波段區(qū)域，然后采用連續(xù)投影算法[18]從特征波段區(qū)域內提取重金屬元素的特征波段。

2.2 潛在污染源空間特征量化

礦山開采、金屬冶煉過程中產生的無機污染物擴散以點源污染為主，通常以開采、冶煉廠為中心向四周擴散，同時受到風向、風速、地形以及降水等自然因子的影響，這些污染物擴散空間影響因子對于土壤污染濃度的分布產生直接影響，因此考慮引入適當的污染源-匯空間特征對上述污染擴散模型進行量化。

研究區(qū)內主要分布有較大的鉛鋅、鋼鐵等金屬冶煉廠5家，位置分布如圖1所示，這些冶煉廠是造成研究區(qū)土壤重金屬污染的潛在污染源。污染源污染物向四周擴散，樣本點與這些潛在污染源的距離以及方位關系是影響樣本點位置污染物累積量的重要因子，因此，主要選取了污染源與采樣點的空間距離和方位角作為空間特征因子。針對每個污染源分別計算以下2個空間特征：

(1)污染源與采樣點的空間距離

空間上距離污染源越近，一般污染物累積越多，因此引入空間距離特征描述距離因子對于污染源污染物向四周擴散的影響。根據采樣點與每個污染源的平面坐標采用歐氏距離公式進行計算，距離污染源越遠，受到污染源的影響越小，因此將該特征量化為距離的倒數進行建模。

(2)污染源與采樣點連線的方位角

在偏離每個污染源的不同方位，由于風向、風速的差異，污染物擴散會出現明顯變化，如果采樣點位于污染源主導風向的下風向，其污染物積聚濃度更高，因此，采樣點與污染源的方位角信息可以模擬在污染源不同方向上大氣擴散條件的差異。方位角空間特征是指以污染源O為起點、樣本點A為終點的連線與正北方向的夾角，描述了樣本采樣點與污染源的空間方位關系。方位角α具體計算步驟如下：首先計算潛在污染源O(x0,y0)與樣本M(xi,yi)之間的象限角β，計算公式為

(1)

其中

Δx=xi-x0Δy=yi-y0

如圖3所示，根據計算得到Δx、Δy來判斷象限角β位于第幾象限，并以此來計算方位角：①當Δx>0，Δy>0時，角β位于第Ⅰ象限，方位角α=β。②當Δx<0，Δy>0時，角β位于第Ⅱ象限，方位角α=180°-β。③當Δx<0，Δy<0時，角β位于第Ⅲ象限，方位角α=180°+β。④當Δx>0，Δy<0時，角β位于第Ⅳ象限，方位角α=360°-β。

圖3 方位角計算Fig.3 Azimuth calculation

在具體計算時，針對5個潛在污染源分別計算距離倒數因子和方位角因子2個空間特征，因此，每個樣本計算10個污染源空間特征。

2.3 污染擴散地形影響因子選取

地形地勢決定了水流流向、流速，對局部氣流風向、風速等也具有一定控制作用，從而對重金屬污染物的擴散和匯集產生影響，引入高程(Elevation)、坡度(Slope)、坡向(Aspect)、坡長因子(LS factor，LSF)[19]、形態(tài)特征(Morphometric features，MF)[20]、廣義表面指數(Generalized surface index，GSI)[20]、風效指數(Wind exposition，WE)[20]、地形濕度指數(Topographic wetness index，TWI)[21]，用于分析地形因子對于土壤重金屬污染物濃度的影響。以上8個地形因子均基于DEM數據計算，并通過插值提取相應樣本點位置特征值，DEM數據采用SRTM數據，網格尺寸為30 m。

2.4 土壤樣本重金屬元素含量測定

研究區(qū)土壤重金屬污染以Pb和Cd污染為主，重點針對這2種元素進行分析。土壤重金屬Pb、Cd含量采用XSERIES-2型電感耦合等離子體質譜儀測定，樣本的統(tǒng)計信息如表1所示。Cd含量與Pb含量之間相關系數達到0.825，兩者具有較高的同源性，發(fā)生協同作用的可能性較大。

表1 研究區(qū)重金屬Pb、Cd含量數據基本統(tǒng)計信息Tab.1 Basic content statistics of heavy metal Pb and Cd of study area

統(tǒng)計數據顯示重金屬Pb達到重污染的樣本占全部樣本的19.3%，重金屬Cd達到重污染的樣本占全部樣本的75.5%，研究區(qū)內存在較為嚴重的Pb、Cd污染，需全面加強對該地區(qū)土壤重金屬污染的監(jiān)測。

將采集的249個樣本按照約3∶1的比例隨機劃分為建模集和測試集，其中建模集186個樣本，測試集63個樣本。

3 研究方法

3.1 極端隨機樹

極端隨機樹類似于隨機森林方法，是一種由多棵決策樹構成的集成學習方法。隨機森林采用隨機采樣來選擇樣本集作為每個決策樹的訓練集[22]，該方法不能保證所有樣本能被充分利用，并且各決策樹之間可能存在相似性?；谝陨峡紤]，GEURTS等[23]提出極端隨機樹模型。

在極端隨機樹中，每棵決策樹均采用全部訓練集，訓練樣本的利用率高，能在一定程度上減少最終預測偏差；為了保證每棵決策樹間的結構差異，極端隨機樹在節(jié)點拆分時引入了更大的隨機性：從子數據集中隨機選取每個特征的判斷閾值，并選擇拆分效果最好的特征作為最優(yōu)判斷屬性。由于節(jié)點拆分判斷閾值的隨機性，極端隨機樹的泛化能力一般會優(yōu)于隨機森林方法。

一般以節(jié)點的不純度作為最優(yōu)判斷屬性的選取依據[24]，回歸類問題衡量節(jié)點不純度的函數一般選擇均方誤差(MSE)或平均絕對誤差(MAE)。選用MSE作為函數的節(jié)點不純度(G)，計算公式為

G(ui,vij)=

(2)

式中ui——某一個節(jié)點判斷屬性

vij——判斷屬性的取值

NS——當前節(jié)點所有訓練樣本個數

Xleft、Xright——左、右子節(jié)點的訓練樣本集合，yi∈Xleft，yj∈Xright

經過k輪訓練得到k棵結構不同的決策樹，最后通過投票或取平均的方式集合不同決策樹的預測結果hi(x)，得到模型的最終結果H(x)。在回歸類問題中，常采用平均的方式計算模型的最終結果，即

(3)

3.2 置換重要性指數

隨機森林和極端隨機樹算法可以基于不純度測量每個特征對模型預測的相對重要性，這種基于不純度計算特征重要性傾向于夸大連續(xù)特征或高基數屬性特征的重要性，另一種特征重要性計算方法是置換重要性(Permutation importance，PI)[22]，該指數是通過觀察每個預測屬性的隨機重排對模型預測精度的影響來直接衡量特征的重要性。

該方法計算過程為：首先訓練基線模型，并通過驗證集記錄R2得分為基準評分S。然后選定數據集中的一個特征要素Fj，打亂順序重新排列其屬性值為Fm,j(m表示M次打亂數據中某一次)，利用修改后數據集重新建立預測模型，通過驗證集計算R2得分Sm,j。特征重要性PFj是基準評分S與屬性值重新排列后數據集構建模型的評分之間的差異。公式為

(4)

PI指數計算方便，特征重要性評價準確，可解釋性較好。

3.3 反演模型精度評價

模型的精度評價采用驗證集的決定系數(R2)、均方根誤差(Root mean square error，RMSE)、相對分析誤差(Relative percent deviation，RPD)以及訓練集交叉驗證得分(Cross validate score，CVS)等。

相對分析誤差ERPD的計算公式為[25]

(5)

式中σ——驗證集樣本的標準差

e——均方根誤差

一般當ERPD<1.4時，模型無法對樣品進行預測；當1.4≤ERPD<2.0時，模型精度一般，具有粗略評估樣品的能力；當ERPD≥2.0時，模型具有較好的預測能力。

4 實驗與結果分析

針對光譜變換特征、污染源空間特征以及地形特征進行組合實驗，分別采用光譜、光譜與地形、光譜與空間以及光譜、空間與地形的特征組合進行實驗，分析不同建模特征的置換重要性，評價光譜特征以及污染擴散空間影響因子選取的有效性。反演模型同時選取了多元線性回歸、支持向量機、隨機森林、梯度提升決策樹(Gradient boosting decision tree，GBDT)等回歸模型作為參考，評估極端隨機樹估算模型的有效性和先進性。

4.1 特征重要性分析

4.1.1光譜特征分析

原始光譜反射率測量值與土壤重金屬含量相關系數較低，計算后的變換特征相關性顯著提高，其中CR1765(CR表示連續(xù)統(tǒng)去除光譜變換，1765表示波長位置為1 765 nm，其他編號含義相同)特征與土壤Pb含量最大相關系數達到-0.70。分別計算光譜變換特征與土壤Pb、Cd含量的相關性，對相關系數進行P=0.01水平上的假設性檢驗，將通過假設性檢驗的波段集合作為特征波段區(qū)域。基于SPA算法在特征波段區(qū)域內分別篩選不同元素共線性最小的有效特征波段組合，其中與Pb相關的篩選光譜特征為72個，與Cd相關的篩選光譜特征為65個。

圖4為單獨使用光譜變換特征構建Pb、Cd元素ERT估算模型時分析得到的置換重要性計算結果。每個特征重新排列計算10次，然后對計算所得PI值的均值和方差進行排序，在此展示的是PI均值最高的15個特征及其統(tǒng)計值。其中Pb元素重要性評價最高的特征為CR2262，其次為CR2174和SD1345，這3個特征PI值明顯高于其他特征。Cd元素重要性評價最高的特征為FD1802，該特征PI值明顯高于其他特征，剩余特征差異較小，其中PI值較高的包括MSC1799和SNV1729等。由統(tǒng)計結果可以看到，PI值較高的波段多為近紅外波段，與已有研究成果[2]較為吻合。

圖4 建模特征為光譜時ERT模型計算特征PI統(tǒng)計結果Fig.4 Statistical results of features PI calculated by ERT model when modeling features were spectrum

表2統(tǒng)計了不同回歸模型和建模特征土壤重金屬Pb、Cd含量反演測試集精度，其中僅使用光譜特征時Pb元素的ERT模型R2可達0.861，RPD為2.686，具有較高的定量反演精度，Cd元素的ERT模型R2可達0.736，RPD為1.945，具有粗略的預測能力，說明對于土壤重金屬元素含量反演，光譜特征雖然為弱信息，但由于污染物擴散中形成的落塵等改變了土壤組分和性狀，從而在土壤光譜特征中表現出來，變換后的土壤光譜特征能夠在一定程度上反映這種污染程度。

表2 不同回歸模型和建模特征土壤重金屬Pb、Cd含量反演精度對比Tab.2 Comparison of precision of soil heavy metal Pb and Cd inversion with different regression models and modeling features

4.1.2污染擴散地形影響因子分析

地形特征對于污染物的擴散和匯聚能夠產生一定的影響，從表2可以看到，當建模特征中加入8個地形特征后，Pb和Cd的建模精度均有明顯提升，Pb的ERT模型R2由0.861提升至0.912，Cd的ERT模型R2由0.736提升至0.800，其他統(tǒng)計值也得到了有效提升，說明了地形特征的有效性。

建模特征為光譜和地形組合時ERT模型計算特征置換重要性統(tǒng)計結果顯示，除光譜特征外地形廣義表面指數(DTM_GS)、高程(DTM_E)、坡度(DTM_S)和風效指數(DTM_WE)因子具有較高的PI值，說明這些地形特征能夠較好地反映地形對污染物的擴散和匯聚產生的影響。其中廣義表面指數、高程、坡度等因子與水流流向、流速以及土壤侵蝕與堆積等相關性較強，影響土壤污染物的運移與擴散。風效指數能夠在一定程度上反映地形對大氣污染物擴散產生的影響。

4.1.3污染源空間特征分析

污染源空間特征能夠較好地表征污染物擴散濃度分布，從表2可以看出，當建模特征中加入10個污染源空間特征后，Pb和Cd的建模精度均有極大提升，Pb的ERT模型R2由0.861提升至0.963， RMSE由43.185 mg/kg下降到22.301 mg/kg，下降了48.36%。Cd的ERT模型R2由0.736提升至0.914， RMSE由0.738 mg/kg下降到0.371 mg/kg，下降了49.73%，其他統(tǒng)計值也提升明顯，充分說明了引入污染源空間特征的有效性。

建模特征為光譜和污染源空間特征組合時ERT模型計算特征置換重要性統(tǒng)計結果顯示，10個污染源空間特征均在PI值最高的前15個特征中，且Pb模型的前7個特征、Cd模型的前9個特征均為污染源空間特征，從側面說明引入空間特征的有效性。該PI值也在一定程度上反映了不同污染源對于土壤重金屬污染的貢獻程度，研究區(qū)SB和WY 2個潛在污染源(圖5中SB、WY表示某冶煉廠，DIST表示空間距離倒數特征，ANGLE表示空間方位角特征)的貢獻度較高。

圖5 建模特征為光譜、污染源空間特征和地形特征組合時ERT模型計算特征PI統(tǒng)計結果Fig.5 Statistical results of features PI calculated by ERT model when modeling features were combination of spectrum, spatial features of pollution sources and topographic features

4.1.4多特征組合分析

當建模特征同時加入光譜、地形和污染源空間特征時，表2統(tǒng)計結果表明，建模精度與使用光譜和污染源空間特征組合時基本相同，變化較小，說明污染源空間特征優(yōu)勢更為明顯，地形因子與污染源空間因子有一定重疊。圖5多特征建模置換重要性統(tǒng)計結果也顯示出同樣的特點，整體PI值最高的特征為污染源空間特征，光譜特征次之，地形特征整體PI值最低，重要性最弱。

利用光譜、空間和地形特征組合建模時篩選出的PI值大于0.002的特征進行建模實驗，其中Pb選取特征15個，Cd選取特征14個。表2統(tǒng)計結果表明，篩選特征建模精度與使用全部特征時精度極為接近，Pb的ERT模型R2均為0.964，Cd的ERT模型R2分別為0.923和0.928，說明利用置換重要性進行特征篩選是有效的。

4.2 極端隨機樹模型估算精度分析

為評價土壤重金屬污染極端隨機樹ERT估算模型的先進性，選取MLR、SVM、RF、GBDT等回歸模型作為對比，表2測試集反演模型精度評價統(tǒng)計結果顯示，ERT模型在不同重金屬元素、不同特征集的反演建模中均取得了最優(yōu)精度，Pb的ERT模型的測試集R2達0.964，Cd的ERT模型R2為0.923，模型穩(wěn)定性最佳。整體上MLR模型反演精度最低，增加污染擴散空間特征時，Pb的MLR模型精度提升較明顯，Cd元素的MLR模型精度提升不大，MLR模型魯棒性較差。SVM模型與MLR模型相反，Pb的MLR模型精度提升不大，Cd的MLR模型精度提升明顯，但是模型穩(wěn)定性弱。RF和GBDT反演模型精度接近，優(yōu)于MLR和SVM模型，增加地形、空間特征時，反演模型精度均得到較大提升，模型魯棒性較高。

訓練集交叉驗證得分CVS對模型的泛化能力有較好的評價，在表2 CVS統(tǒng)計中ERT模型的優(yōu)勢也較為明顯，均優(yōu)于RF和GBDT模型。圖6給出了測試集實測真實值與ERT模型預測值的序列分布情況，結果表明重金屬Pb和Cd含量預測值與真實值的吻合度較高。

圖6 測試集實測真實值與ERT模型預測值序列分布Fig.6 Sequence distributions of measured true value of test set and predicted value of ERT model

5 結論

(1)僅使用光譜特征構建的Pb、Cd ERT估算模型具有較高的R2和RPD，說明變換后的土壤光譜特征能夠在一定程度上反映這種污染程度。其中Pb置換重要性評價最高的特征為2 262 nm連續(xù)統(tǒng)去除光譜變換特征，Cd重要性評價最高的特征為1 802 nm一階微分光譜變換特征。

(2)當建模特征中加入地形特征后，Pb和Cd的建模精度均有明顯提升，置換重要性統(tǒng)計結果顯示地形廣義表面指數、高程、坡度和風效指數等特征具有較高的PI值，說明這些地形特征能夠較好地反映地形對污染物的擴散和累積產生的影響。

(3)當建模特征中加入污染源空間特征后，Pb和Cd的建模精度均有極大提升，各項統(tǒng)計值改善明顯，充分說明了所提出構建污染擴散影響因子的有效性。污染源空間特征重要性分析也可以在一定程度上反映不同污染源對于土壤重金屬污染的貢獻度。

(4)光譜、地形和污染源空間特征組合建模結果表明PI值最高的特征為污染源空間特征，光譜特征次之，地形特征整體PI值最低。使用置換重要性指數優(yōu)選特征建立的估測模型與使用全部特征時建模最優(yōu)精度極為接近，說明了置換重要性指數特征篩選方法的有效性。

(5)與MLR、SVM、RF、GBDT等回歸模型對比，ERT估算模型在各項指標評價中優(yōu)勢明顯，其中Pb的ERT模型的測試集R2達0.964，Cd的ERT模型R2為0.923，ERT土壤重金屬估算模型估算精度較高，表明該方法反演土壤重金屬含量具有較高的可行性。

(6)提出構建潛在污染源空間特征量化污染物擴散空間影響因子，該方法適用于污染物來源較為明確的點源、線源污染類型，一般土壤重金屬污染物來源于礦業(yè)活動、冶金以及工業(yè)生產中形成的大氣降塵、污水排灌等，污染源明確，因此本文提出方法具有較高的推廣應用潛力。