馮鋒，王育紅，左雨芳

（江蘇師范大學(xué)地理測繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇 徐州 221000）

鎘（Cadmium，Cd）是一種重金屬元素，具有非生物降解性，在土壤中富集不但改變土壤理化性質(zhì)，直接降低農(nóng)作物產(chǎn)量，而且造成農(nóng)作物Cd 含量超標(biāo)。Cd 在經(jīng)食物鏈進(jìn)入人體后，會導(dǎo)致人體骨骼嚴(yán)重軟化，并對肝臟、腎臟造成損害[1-2]?！度珖寥牢廴緺顩r調(diào)查公報(bào)》顯示，我國約有16.7%的農(nóng)田土壤遭到重金屬污染，其中40%以上為Cd 污染。因此為有效防治土壤重金屬污染，我國相繼制定出臺了《土壤污染防治行動計(jì)劃》《中華人民共和國土壤污染防治法》等一系列政策、計(jì)劃和法律法規(guī)，以期“守護(hù)土壤健康，助力高質(zhì)發(fā)展”。

土壤Cd 污染具有極強(qiáng)的空間變異性，即使是相鄰區(qū)域，由于各種因素影響，污染也不盡相同，因此當(dāng)前對于土壤Cd 污染的研究多集中于小尺度的研究區(qū)，并應(yīng)用經(jīng)過實(shí)地布點(diǎn)采樣、實(shí)驗(yàn)室分析獲取的不同點(diǎn)位的Cd 含量進(jìn)行污染特征及分布評價(jià)、污染來源解析以及影響因子分析的模式，而對于中大尺度研究區(qū)，由于區(qū)域面積過大，受限于人力、物力、財(cái)力，常難以實(shí)現(xiàn)[3]。但大尺度的研究分析通常能揭示土壤Cd 污染宏觀上的空間分布規(guī)律，同時(shí)在大尺度上進(jìn)行污染影響因子分析對污染的防治具有重要的決策支撐作用。傳統(tǒng)的土壤重金屬污染空間分布預(yù)測方法包括反距離加權(quán)、徑向基函數(shù)、克里格插值等，但上述方法僅根據(jù)未知點(diǎn)到樣本點(diǎn)的距離估算污染狀況，而重金屬在土壤內(nèi)的富集受多種因子影響，包括眾多自然因子（成土母質(zhì)、成土過程、土壤質(zhì)地等）和人為因子（交通污染、工業(yè)污染、農(nóng)業(yè)施肥灌溉等），因此地統(tǒng)計(jì)學(xué)等插值方法預(yù)測結(jié)果通常精度較低[4-5]。近年來，隨著機(jī)器學(xué)習(xí)的成熟和普及，支持向量機(jī)、隨機(jī)森林、神經(jīng)網(wǎng)絡(luò)等算法被應(yīng)用于土壤重金屬污染的空間分布研究中，例如：胡夢珺等[6]采用隨機(jī)森林模型預(yù)測了蘭州市主城區(qū)校園地表灰塵中重金屬污染程度，并與傳統(tǒng)的空間插值結(jié)果進(jìn)行了對比；范俊楠等[7]利用BP 神經(jīng)網(wǎng)絡(luò)對湖北省重點(diǎn)區(qū)域行業(yè)企業(yè)周邊的土壤重金屬污染進(jìn)行了預(yù)測；Omondi 等[8]使用隨機(jī)森林模型對內(nèi)羅畢和特里爾卡河匯流區(qū)的土壤重金屬進(jìn)行了預(yù)測，并計(jì)算了不同影響因子的貢獻(xiàn)率。

盡管上述方法相較于傳統(tǒng)的預(yù)測方法極大提升了準(zhǔn)確度，但由于不同研究區(qū)內(nèi)的重要影響因子具有較大差異，在構(gòu)建預(yù)測模型時(shí)將影響力極小的影響因子輸入模型，將導(dǎo)致預(yù)測精度降低。基于此，本研究面向貴州省貴陽、遵義和畢節(jié)三市的土壤Cd污染，初選20 個原始影響因子，創(chuàng)新性地利用隨機(jī)森林（Random Forest，RF）對原始影響因子進(jìn)行篩選，再利用極端梯度提升（eXtreme Gradient Boosting，XGBoost）算法對篩選后的影響因子進(jìn)行訓(xùn)練，提高預(yù)測精度，以充分了解三市土壤Cd污染空間分布狀況。

1 材料與方法

1.1 研究區(qū)概況

貴陽市、遵義市和畢節(jié)市位于貴州西部，地勢南高北低，喀斯特地貌面積占一半以上。三市礦產(chǎn)資源豐富，煤、鐵、磷等的儲量在全國名列前茅。長期以來，三市的經(jīng)濟(jì)發(fā)展依托于礦產(chǎn)資源的開采，因此三市是受到土壤Cd污染的典型研究區(qū)。

1.2 研究方法

1.2.1 地累積指數(shù)法

地累積指數(shù)又名Muller 指數(shù)，于1969 年提出，此后被廣泛用于水、土壤、沉積物等環(huán)境中的重金屬污染評價(jià)[9]。其計(jì)算公式為：

式中：Igeo為地累積指數(shù)；Ci為重金屬i的實(shí)測值，mg·kg-1；Bi為重金屬i的地球化學(xué)背景值，mg·kg-1；1.5 為考慮環(huán)境差異造成的背景值浮動而加入的修正常數(shù)。根據(jù)Igeo得分將土壤重金屬污染劃分為7個等級，具體標(biāo)準(zhǔn)如表1所示。

表1 地累積指數(shù)與污染程度等級劃分Table 1 The geo-accumulation index and classification of pollution degree

1.2.2 RF算法

RF 是2001 年由Breiman[10]在貝爾實(shí)驗(yàn)室創(chuàng)立的隨機(jī)決策森林方法基礎(chǔ)上提出的一種集成學(xué)習(xí)算法。該算法通過自助法（bootstrap）重采樣，從原始的訓(xùn)練數(shù)據(jù)內(nèi)抽取若干樣本構(gòu)建出多個弱學(xué)習(xí)器——分類回歸樹（Classification and Regression Tree，CART），將其組合匯總后生成強(qiáng)學(xué)習(xí)器以解決分類和回歸兩類預(yù)測問題。

除用于預(yù)測外，RF也被廣泛用于評估影響因子的特征重要性，即貢獻(xiàn)率。特征重要性求解主要基于決策樹構(gòu)建時(shí)選擇節(jié)點(diǎn)的指標(biāo)（Gini指數(shù)），先將計(jì)算得到的每個特征在每一棵決策樹上進(jìn)行節(jié)點(diǎn)分割時(shí)的Gini指數(shù)差值匯總?cè)∑骄岛螅儆?jì)算其占所有特征Gini指數(shù)的總變化值的百分比，該值即為特征重要性[11]。

假設(shè)一個訓(xùn)練好的RF 模型由T棵CART 組成，RF 所用訓(xùn)練樣本共包含F(xiàn)個特征自變量（X1,X2,X3,…,XF），其因變量共有K個類別取值。對于RF中任意一棵編號為t（1≤t≤T）的CART，假設(shè)其共包含N個節(jié)點(diǎn)，則該樹任意節(jié)點(diǎn)n（1≤n≤N）的Gini 指數(shù)（GItn）可用式（2）計(jì)算。式中，pnk表示節(jié)點(diǎn)n上樣本屬于第k類（1≤k≤K）的經(jīng)驗(yàn)概率值。

則第i個特征在節(jié)點(diǎn)n上分裂前后的Gini 指數(shù)差值為：

式中：VIMin表示特征i在節(jié)點(diǎn)n上的Gini 指數(shù)變化量；GIn表示節(jié)點(diǎn)n的Gini 指數(shù)；GIl表示節(jié)點(diǎn)n分裂后的左節(jié)點(diǎn)Gini 指數(shù)；GIr表示節(jié)點(diǎn)n分裂后的右節(jié)點(diǎn)Gini 指數(shù)。

假設(shè)特征i在決策樹t中出現(xiàn)的節(jié)點(diǎn)集合為N，則可求得特征i在決策樹t中的Gini指數(shù)變化量（VIMti）：

同理可得在所有決策樹T中，特征i的Gini 指數(shù)變化量總和（VIMi）：

最后對所求得的第i個特征的Gini指數(shù)變化量進(jìn)行歸一化，得到其特征重要性得分（VIMi′）：

RF 采用有放回的采樣并按照預(yù)定數(shù)目隨機(jī)選擇參與決策樹的特征，因此其對異常值和噪聲具有較高的容忍度，不容易出現(xiàn)過擬合，整體上具有較強(qiáng)的泛化能力、數(shù)據(jù)挖掘能力和預(yù)測準(zhǔn)確率，曾被譽(yù)為代表集成學(xué)習(xí)技術(shù)水平的最好算法之一[12]。

1.2.3 XGBoost算法

XGBoost 是一種優(yōu)化的分布式梯度提升樹算法，其基本思想是通過不斷增加新的決策樹參與訓(xùn)練以擬合預(yù)測值與真實(shí)值之間的殘差，并利用集成思想獲得最終的預(yù)測值[13]。其計(jì)算公式如下：

式中：fk為第k個基學(xué)習(xí)器；xi為第i個樣本的特征向量為第i個樣本的預(yù)測值。

為提高預(yù)測精度，XGBoost 算法構(gòu)建損失函數(shù)L代表模型的偏差，同時(shí)與梯度提升樹不同的是，XGBoost 算法引入正則項(xiàng)Ω 以抑制模型復(fù)雜度，故最終得到的目標(biāo)函數(shù)（Obj）為：

式中：Ω(fk)為第k棵決策樹的正則項(xiàng)；l(yi,)為第i個樣本的預(yù)測誤差。

假設(shè)XGBoost 訓(xùn)練完成后共生成了K棵決策樹，則將每個樣本結(jié)構(gòu)映射到每棵決策樹的葉節(jié)點(diǎn)上的分?jǐn)?shù)相加即可得到預(yù)測值，公式如下：

式中：F為預(yù)測值；f（x）為某一棵決策樹的模型；ωq（x）為決策樹q的所有葉節(jié)點(diǎn)分?jǐn)?shù)組成的集合；T為決策樹q的葉節(jié)點(diǎn)數(shù)量。

1.3 數(shù)據(jù)描述

本研究土壤Cd 含量數(shù)據(jù)來源于文獻(xiàn)[14]提供的附件資料，該附件以經(jīng)緯度及數(shù)值方式整理了我國土壤As、Cu、Pb、Cr、Zn和Cd含量數(shù)據(jù)，經(jīng)篩選后共獲得我國2006—2016 年內(nèi)公開發(fā)表的2 450 篇文獻(xiàn)內(nèi)土壤重金屬含量數(shù)據(jù)，收集的文獻(xiàn)中的樣本點(diǎn)布設(shè)方法大多基于隨機(jī)采樣，土壤樣品數(shù)據(jù)采集深度為0～20 cm，土壤顆粒大小集中于0.5～4.0 mm，Cd 含量采用原子吸收分光光度法測定。提取附件內(nèi)貴陽市、遵義市和畢節(jié)市的土壤Cd 含量數(shù)據(jù)進(jìn)行預(yù)處理，根據(jù)Zscore 剔除離群值點(diǎn)，最終得到170 個Cd 樣本點(diǎn)數(shù)據(jù)，其分布如圖1所示。

圖1 樣本點(diǎn)分布圖Figure 1 Sample points distribution map

通過閱讀文獻(xiàn)[15-19]以及考慮到數(shù)據(jù)獲取的難易程度，本研究選擇土壤質(zhì)地、年平均氣溫、植被指數(shù)、高程、河流、土地利用類型、國民生產(chǎn)總值、人口密度、道路等影響因子，為便于數(shù)據(jù)的查找和使用，將數(shù)據(jù)組織成如下格式（表2）：OID 為每個樣本點(diǎn)的唯一ID編號；CON 為Cd含量（mg·kg-1）；F01～F13為重金屬含量自然影響因子；F14～F20 為人為影響因子。F01 為黏土含量；F02 為砂土含量；F03 為粉砂土含量；F04為土壤侵蝕度；F05 為年均氣溫；F06 為年均濕潤指數(shù)；F07 為年均干燥度指數(shù)；F08 為歸一化植被指數(shù)；F09 為高程；F10 為坡向；F11 為坡度；F12 為與一級河流距離；F13 為與二級河流距離；F14 為土地利用類型；F15 為人口密度；F16 為與主干道距離；F17 為與次干道距離；F18 為與高速路距離；F19 為與鐵路距離；F20 為人均GDP。各影響因子空間分布特征如圖2所示。

圖2 影響因子分布圖Figure 2 Influence factors distribution map

表2 樣本點(diǎn)數(shù)據(jù)基本結(jié)構(gòu)與形式Table 2 Basic structure and form of sample points data

2 結(jié)果與討論

2.1 描述性統(tǒng)計(jì)分析

研究區(qū)內(nèi)Cd 含量最小值為0.02 mg·kg-1，最大值為4.10 mg·kg-1，平均值為0.68mg·kg-1，標(biāo)準(zhǔn)差為0.86，變異系數(shù)為125.37%。貴州省的土壤Cd 背景值為0.66 mg·kg-1，貴陽、遵義和畢節(jié)三市的土壤Cd 含量平均值僅比背景值高出0.02 mg·kg-1，表明整體上研究區(qū)內(nèi)土壤Cd污染程度較低。變異系數(shù)反映樣本分布的空間均衡性，研究區(qū)內(nèi)Cd 的變異系數(shù)為125.37%，揭示研究區(qū)內(nèi)Cd 的分布極為不均衡，受到了較大的外源影響。

2.2 基于RF的影響因子貢獻(xiàn)率分析

對收集的170 個樣本點(diǎn)計(jì)算地累積指數(shù)，分級結(jié)果如表3 所示。研究區(qū)內(nèi)77.65%的樣本點(diǎn)無污染，輕微污染、中度污染和中強(qiáng)污染樣點(diǎn)分別占14.71%、5.88%、1.76%。

表3 研究區(qū)土壤Cd污染地累積指數(shù)分級Table 3 Classification of soil Cd based on geo-accumulation index in the study area

依據(jù)Igeo等級為樣本點(diǎn)設(shè)置標(biāo)簽，Ⅰ、Ⅱ、Ⅲ、Ⅴ級的標(biāo)簽分別為0、1、2、3。對170 個樣本點(diǎn)，選擇其中70%作為訓(xùn)練集，余下的30%作為測試集導(dǎo)入RF 模型內(nèi)，模型的主要參數(shù)如下：n_estimators=100，criterion=‘gini’，max_features=‘sqrt’，random_state=2022，n_jobs=-1，基于python 語言利用PyCharm 編譯器實(shí)現(xiàn)模型。得到的各影響因子的貢獻(xiàn)率如圖3所示，可知：

圖3 各影響因子貢獻(xiàn)率Figure 3 Contribution rate of each influencing factor

（1）貢獻(xiàn)率排名前三的影響因子分別為土壤侵蝕度（0.100）、高程（0.088）和年均氣溫（0.084）。土壤侵蝕度對土壤Cd 污染的貢獻(xiàn)率最高，表明其與研究區(qū)土壤Cd 污染的分布關(guān)系最為密切，原因在于研究區(qū)內(nèi)多為喀斯特地貌，該地貌具有土層淺薄、土被連續(xù)性差的特點(diǎn)，這極大降低了土壤對于Cd 元素的吸收和吸附能力，在雨水作用下，含Cd的顆粒懸浮物易隨泥沙和地表徑流進(jìn)行遷移；高程對土壤Cd 污染的貢獻(xiàn)率位居第二，主要原因在于研究區(qū)內(nèi)山地較多，地勢起伏變化大，而相對較高的地形對Cd 的遷移擴(kuò)散起到阻隔作用；年均氣溫對土壤Cd 污染的貢獻(xiàn)率位居第三，原因在于貴州境內(nèi)氣溫較高，地區(qū)降水豐富、地表徑流量大，進(jìn)而導(dǎo)致Cd擴(kuò)散。

（2）除上述3 項(xiàng)影響因子外，貢獻(xiàn)率超過0.05 的影響因子還有人均GDP（0.079）、人口密度（0.078）、與高速路距離（0.067）、與主干道距離（0.060）以及年均干燥度指數(shù)（0.052）。其中4項(xiàng)為人為影響因子，說明研究區(qū)內(nèi)土壤Cd 污染受到人為影響較為嚴(yán)重，印證了研究區(qū)內(nèi)Cd 含量變異系數(shù)大，是受到了較大的外源性影響，且交通為主要污染源之一。年均干燥度指數(shù)是蒸發(fā)量與降水量的比值，用于衡量地區(qū)的氣候干燥程度，環(huán)境越干燥，表明蒸發(fā)量越大，則土壤內(nèi)的水分越少，從而抑制Cd隨水分?jǐn)U散。

（3）其余12 項(xiàng)影響因子貢獻(xiàn)率均低于0.05，包括植被、土壤質(zhì)地、河流等，說明在研究區(qū)內(nèi)，這些因子與土壤Cd污染的空間分布關(guān)聯(lián)性較低。

（4）研究區(qū)內(nèi)對土壤Cd 污染貢獻(xiàn)率排名前三的均為自然因子，且貢獻(xiàn)率超過0.05 的影響因子中，自然影響因子（0.324）高于人為影響因子（0.284），說明在大尺度研究區(qū)內(nèi)，自然因子對土壤Cd 污染分布起決定性作用。

2.3 影響因子篩選

準(zhǔn)確率（Accuracy，ACC）是指測試樣本中模型預(yù)測“正確”的樣本的占比，其值在[0，1]范圍內(nèi)，值越大，模型的分類結(jié)果越準(zhǔn)確；Kappa系數(shù)用于判斷分類模型分類結(jié)果的好壞，是衡量模型分類精度的重要指標(biāo)，其值在[-1，1]范圍內(nèi)，值越接近1，分類結(jié)果越好。將兩者結(jié)合進(jìn)行影響因子的篩選有助于降低評估的誤差。其計(jì)算公式如下：

式中：TP、TN、FP、FN分別為真正、真負(fù)、假正、假負(fù)的樣本數(shù)。

將所有影響因子的貢獻(xiàn)率按從低到高進(jìn)行排序，從貢獻(xiàn)率最低的因子進(jìn)行刪除，每刪除一個影響因子則重新構(gòu)建一個新的RF 模型。以本研究170 個樣本點(diǎn)作為基礎(chǔ)，在保證訓(xùn)練集、測試集和模型各參數(shù)不變的前提下，迭代構(gòu)建RF模型，并記錄每一個模型的ACC和Kappa系數(shù)，探索最優(yōu)的影響因子集合。實(shí)驗(yàn)記錄的ACC和Kappa系數(shù)如表4所示。

表4 實(shí)驗(yàn)過程記錄Table 4 Record of experimental process

由表4 可知，從實(shí)驗(yàn)1 到實(shí)驗(yàn)8，隨著影響因子輸入的減少，ACC和Kappa系數(shù)呈現(xiàn)上下波動，并在實(shí)驗(yàn)4時(shí)，ACC和Kappa系數(shù)同時(shí)達(dá)到最大值，說明隨著影響因子的刪除，輸入模型內(nèi)的干擾信息逐步減少，輸入的數(shù)據(jù)集得到優(yōu)化，即刪除黏土含量、粉砂土含量、土地利用類型3 項(xiàng)影響因子后，模型的準(zhǔn)確度最高，一致性最好。

2.4 XGBoost結(jié)果分析

以篩選優(yōu)化后得到的17 個影響因子為基礎(chǔ)，引入XGBoost 模型進(jìn)一步預(yù)測貴陽、遵義和畢節(jié)三市土壤Cd污染分級。經(jīng)反復(fù)測試后模型的主要參數(shù)設(shè)置如下：booster=‘gbtree’，n_estimators=100，num_class=4，max_depth=5，subsample=0.6，min_child_weight=1，learning_rate=0.1，seed=12。為驗(yàn)證RF-XGBoost 模型的性能，除選取ACC、Kappa系數(shù)兩個指標(biāo)衡量模型精度外，再次引入F1_score指標(biāo)同經(jīng)影響因子篩選后的RF 以及未經(jīng)影響因子篩選的XGBoost 模型進(jìn)行穩(wěn)定性對比。F1_score的計(jì)算公式如下：

式中：precision為精確率；recall為召回率。3種模型的性能指標(biāo)對比結(jié)果如表5所示。

表5 3種模型的性能指標(biāo)對比Table 5 Comparison of performance indexes of the three models

由表4 可知，RF-XGBoost 模型的ACC、Kappa系數(shù)和F1_score均高于另外兩種模型，其中ACC均提升了0.039 3，Kappa系數(shù)分別提升了0.059 2、0.091 4，F(xiàn)1_score 分別提升了0.250 4、0.270 1，表明RF-XGBoost 模型在土壤Cd 污染分級的預(yù)測上具有較好的準(zhǔn)確性和穩(wěn)定性。

為充分了解研究區(qū)內(nèi)的Cd 污染空間分布，基于《土壤環(huán)境監(jiān)測技術(shù)規(guī)范》中的規(guī)定，以2.5 km的精度在研究區(qū)內(nèi)均勻布設(shè)41 950個預(yù)測點(diǎn)，利用訓(xùn)練好的RF-XGBoost 模型預(yù)測各點(diǎn)位的污染分級后，使用地統(tǒng)計(jì)學(xué)普通克里格插值獲取研究區(qū)內(nèi)的污染分布情況，由于普通克里格插值的估算數(shù)值為連續(xù)數(shù)值，故按照[0，0.5）、[0.5，1.5）、[1.5，2.5）、[2.5，3）的間隔使用重分類功能對土壤Cd 污染空間分布進(jìn)行調(diào)整，最終得到研究區(qū)內(nèi)Cd 污染的分布圖，如圖4 所示。由圖可知研究區(qū)內(nèi)的Cd 污染主要呈現(xiàn)斑塊狀，但總體污染程度較低，污染主要集中于畢節(jié)市及遵義市東北方向。其中畢節(jié)市出現(xiàn)較大范圍的中度-中強(qiáng)污染帶，原因在于畢節(jié)市礦產(chǎn)資源豐富，礦業(yè)經(jīng)濟(jì)為該市的支柱經(jīng)濟(jì)，并且畢節(jié)市的礦產(chǎn)開采歷史長，導(dǎo)致Cd長期富集，形成污染。此外，畢節(jié)市土壤侵蝕程度較高，Cd 易隨地表徑流進(jìn)行遷移擴(kuò)散，又因該市的地勢起伏較大，山脈較多，對Cd 遷移擴(kuò)散有阻隔作用，從而導(dǎo)致Cd富集，形成中度-中強(qiáng)污染帶。

3 結(jié)論

（1）研究區(qū)的土壤Cd 含量平均值略高于貴州省的背景值，整體污染程度較低，但土壤Cd污染分布極不均衡，受到較大的外源性影響。

（2）土壤侵蝕度、高程和年均氣溫3 項(xiàng)自然影響因子對研究區(qū)土壤Cd 污染貢獻(xiàn)率最高，揭示在較大尺度研究區(qū)內(nèi)自然環(huán)境是造成土壤Cd富集的主要影響因素；人均GDP、人口密度、與高速路距離、與主干道距離4 項(xiàng)人為影響因子對土壤Cd 污染貢獻(xiàn)率較高，印證研究區(qū)內(nèi)的土壤Cd 含量受到較大外源性影響，并提示在研究區(qū)內(nèi)，人類活動是造成土壤Cd污染的重要來源，其中交通污染源需重點(diǎn)防治。

（3）RF-XGBoost 模型的精度和穩(wěn)定性顯著高于RF、XGBoost 模型，結(jié)合地統(tǒng)計(jì)學(xué)可準(zhǔn)確客觀地預(yù)測土壤Cd 污染空間分布特征，該模型在土壤重金屬污染空間分布研究中具有推廣性，可為宏觀掌握大尺度范圍內(nèi)土壤重金屬污染分布提供參考。