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

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

0 引言

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

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

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

1 研究區(qū)概況

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

2 實(shí)驗(yàn)數(shù)據(jù)獲取與預(yù)處理

2.1 土壤光譜數(shù)據(jù)測(cè)定及預(yù)處理

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

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

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

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

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

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

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

礦山開采、金屬冶煉過(guò)程中產(chǎn)生的無(wú)機(jī)污染物擴(kuò)散以點(diǎn)源污染為主，通常以開采、冶煉廠為中心向四周擴(kuò)散，同時(shí)受到風(fēng)向、風(fēng)速、地形以及降水等自然因子的影響，這些污染物擴(kuò)散空間影響因子對(duì)于土壤污染濃度的分布產(chǎn)生直接影響，因此考慮引入適當(dāng)?shù)奈廴驹?匯空間特征對(duì)上述污染擴(kuò)散模型進(jìn)行量化。

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

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

空間上距離污染源越近，一般污染物累積越多，因此引入空間距離特征描述距離因子對(duì)于污染源污染物向四周擴(kuò)散的影響。根據(jù)采樣點(diǎn)與每個(gè)污染源的平面坐標(biāo)采用歐氏距離公式進(jìn)行計(jì)算，距離污染源越遠(yuǎn)，受到污染源的影響越小，因此將該特征量化為距離的倒數(shù)進(jìn)行建模。

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

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

(1)

其中

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

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

圖3 方位角計(jì)算Fig.3 Azimuth calculation

在具體計(jì)算時(shí)，針對(duì)5個(gè)潛在污染源分別計(jì)算距離倒數(shù)因子和方位角因子2個(gè)空間特征，因此，每個(gè)樣本計(jì)算10個(gè)污染源空間特征。

2.3 污染擴(kuò)散地形影響因子選取

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

2.4 土壤樣本重金屬元素含量測(cè)定

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

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

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

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

3 研究方法

3.1 極端隨機(jī)樹

極端隨機(jī)樹類似于隨機(jī)森林方法，是一種由多棵決策樹構(gòu)成的集成學(xué)習(xí)方法。隨機(jī)森林采用隨機(jī)采樣來(lái)選擇樣本集作為每個(gè)決策樹的訓(xùn)練集[22]，該方法不能保證所有樣本能被充分利用，并且各決策樹之間可能存在相似性。基于以上考慮，GEURTS等[23]提出極端隨機(jī)樹模型。

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

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

G(ui,vij)=

(2)

式中ui——某一個(gè)節(jié)點(diǎn)判斷屬性

vij——判斷屬性的取值

NS——當(dāng)前節(jié)點(diǎn)所有訓(xùn)練樣本個(gè)數(shù)

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

經(jīng)過(guò)k輪訓(xùn)練得到k棵結(jié)構(gòu)不同的決策樹，最后通過(guò)投票或取平均的方式集合不同決策樹的預(yù)測(cè)結(jié)果hi(x)，得到模型的最終結(jié)果H(x)。在回歸類問(wèn)題中，常采用平均的方式計(jì)算模型的最終結(jié)果，即

(3)

3.2 置換重要性指數(shù)

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

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

(4)

PI指數(shù)計(jì)算方便，特征重要性評(píng)價(jià)準(zhǔn)確，可解釋性較好。

3.3 反演模型精度評(píng)價(jià)

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

相對(duì)分析誤差ERPD的計(jì)算公式為[25]

(5)

式中σ——驗(yàn)證集樣本的標(biāo)準(zhǔn)差

e——均方根誤差

一般當(dāng)ERPD<1.4時(shí)，模型無(wú)法對(duì)樣品進(jìn)行預(yù)測(cè)；當(dāng)1.4≤ERPD<2.0時(shí)，模型精度一般，具有粗略評(píng)估樣品的能力；當(dāng)ERPD≥2.0時(shí)，模型具有較好的預(yù)測(cè)能力。

4 實(shí)驗(yàn)與結(jié)果分析

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

4.1 特征重要性分析

4.1.1光譜特征分析

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

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

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

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

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

4.1.2污染擴(kuò)散地形影響因子分析

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

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

4.1.3污染源空間特征分析

污染源空間特征能夠較好地表征污染物擴(kuò)散濃度分布，從表2可以看出，當(dāng)建模特征中加入10個(gè)污染源空間特征后，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)計(jì)值也提升明顯，充分說(shuō)明了引入污染源空間特征的有效性。

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

圖5 建模特征為光譜、污染源空間特征和地形特征組合時(shí)ERT模型計(jì)算特征PI統(tǒng)計(jì)結(jié)果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多特征組合分析

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

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

4.2 極端隨機(jī)樹模型估算精度分析

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

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

圖6 測(cè)試集實(shí)測(cè)真實(shí)值與ERT模型預(yù)測(cè)值序列分布Fig.6 Sequence distributions of measured true value of test set and predicted value of ERT model

5 結(jié)論

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

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

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

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

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

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