基于RF-XGBoost的土壤鎘污染影響因子及空間分布研究

馮鋒，王育紅，左雨芳

（江蘇師范大學地理測繪與城鄉規劃學院，江蘇 徐州 221000）

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

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

盡管上述方法相較于傳統的預測方法極大提升了準確度，但由于不同研究區內的重要影響因子具有較大差異，在構建預測模型時將影響力極小的影響因子輸入模型，將導致預測精度降低?；诖?，本研究面向貴州省貴陽、遵義和畢節三市的土壤Cd污染，初選20 個原始影響因子，創新性地利用隨機森林（Random Forest，RF）對原始影響因子進行篩選，再利用極端梯度提升（eXtreme Gradient Boosting，XGBoost）算法對篩選后的影響因子進行訓練，提高預測精度，以充分了解三市土壤Cd污染空間分布狀況。

1 材料與方法

1.1 研究區概況

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

1.2 研究方法

1.2.1 地累積指數法

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

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

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

1.2.2 RF算法

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

除用于預測外，RF也被廣泛用于評估影響因子的特征重要性，即貢獻率。特征重要性求解主要基于決策樹構建時選擇節點的指標（Gini指數），先將計算得到的每個特征在每一棵決策樹上進行節點分割時的Gini指數差值匯總取平均值后，再計算其占所有特征Gini指數的總變化值的百分比，該值即為特征重要性[11]。

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

則第i個特征在節點n上分裂前后的Gini 指數差值為：

式中：VIMin表示特征i在節點n上的Gini 指數變化量；GIn表示節點n的Gini 指數；GIl表示節點n分裂后的左節點Gini 指數；GIr表示節點n分裂后的右節點Gini 指數。

假設特征i在決策樹t中出現的節點集合為N，則可求得特征i在決策樹t中的Gini指數變化量（VIMti）：

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

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

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

1.2.3 XGBoost算法

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

式中：fk為第k個基學習器；xi為第i個樣本的特征向量為第i個樣本的預測值。

為提高預測精度，XGBoost 算法構建損失函數L代表模型的偏差，同時與梯度提升樹不同的是，XGBoost 算法引入正則項Ω 以抑制模型復雜度，故最終得到的目標函數（Obj）為：

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

假設XGBoost 訓練完成后共生成了K棵決策樹，則將每個樣本結構映射到每棵決策樹的葉節點上的分數相加即可得到預測值，公式如下：

式中：F為預測值；f（x）為某一棵決策樹的模型；ωq（x）為決策樹q的所有葉節點分數組成的集合；T為決策樹q的葉節點數量。

1.3 數據描述

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

圖1 樣本點分布圖Figure 1 Sample points distribution map

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

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

表2 樣本點數據基本結構與形式Table 2 Basic structure and form of sample points data

2 結果與討論

2.1 描述性統計分析

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

2.2 基于RF的影響因子貢獻率分析

對收集的170 個樣本點計算地累積指數，分級結果如表3 所示。研究區內77.65%的樣本點無污染，輕微污染、中度污染和中強污染樣點分別占14.71%、5.88%、1.76%。

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

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

圖3 各影響因子貢獻率Figure 3 Contribution rate of each influencing factor

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

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

（3）其余12 項影響因子貢獻率均低于0.05，包括植被、土壤質地、河流等，說明在研究區內，這些因子與土壤Cd污染的空間分布關聯性較低。

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

2.3 影響因子篩選

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

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

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

表4 實驗過程記錄Table 4 Record of experimental process

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

2.4 XGBoost結果分析

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

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

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

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

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

3 結論

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

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

（3）RF-XGBoost 模型的精度和穩定性顯著高于RF、XGBoost 模型，結合地統計學可準確客觀地預測土壤Cd 污染空間分布特征，該模型在土壤重金屬污染空間分布研究中具有推廣性，可為宏觀掌握大尺度范圍內土壤重金屬污染分布提供參考。