基于空間聚類與隨機(jī)森林的稻米富集鎘影響因素篩選研究

郭新蕾，趙玉杰，劉瀟威，周其文，王夏暉，李志濤，朱智偉，張鐵亮，王祖光，張 璠，孫 揚(yáng)

（1.農(nóng)業(yè)農(nóng)村部農(nóng)產(chǎn)品質(zhì)量安全環(huán)境因子控制重點(diǎn)實(shí)驗(yàn)室，天津 300191；2.農(nóng)業(yè)農(nóng)村部環(huán)境保護(hù)科研監(jiān)測(cè)所，天津 300191；3.生態(tài)環(huán)境部環(huán)境規(guī)劃院，北京 100012；4.中國(guó)水稻研究所，杭州 310006）

全國(guó)土壤污染普查結(jié)果表明，我國(guó)土壤重金屬污染問(wèn)題突出，尤其Cd污染在局部地區(qū)較為嚴(yán)重[1]。土壤Cd會(huì)通過(guò)植物吸收進(jìn)入食物鏈對(duì)人體健康造成威脅[2]。稻米是我國(guó)的主糧，相關(guān)研究表明，水稻相對(duì)小麥、玉米等農(nóng)作物而言更易富集Cd[3]。因此，防控稻米Cd污染是保障我國(guó)農(nóng)產(chǎn)品質(zhì)量安全的關(guān)鍵環(huán)節(jié)。稻米對(duì)Cd的富集受多種因素影響，主要包括：土壤pH、土壤有機(jī)質(zhì)、鐵錳氧化物、土壤Ca含量，以及S、P、Cl等陰離子[4]。大田條件下稻米對(duì)Cd的吸收受多種因素的綜合影響，但總體而言各因子的影響能力不同，因此鎖定某區(qū)域影響稻米Cd富集的主控因素是該區(qū)域稻米Cd污染治理的關(guān)鍵。探查影響稻米富集Cd主控因素的方法較多，目前采用的方法主要有相關(guān)性分析、多元線性回歸、偏最小二乘回歸[5]、聚類分析[6]、主成分分析[7]、Cubist混合線性回歸[8]等。

影響稻米富集Cd的土壤環(huán)境因子具有區(qū)域相似性及區(qū)域變異性，因此稻米對(duì)Cd的富集能力也存在明顯的區(qū)域差別，而以往的研究方法很少考慮稻米Cd富集的區(qū)域差異性。且僅基于點(diǎn)數(shù)據(jù)進(jìn)行的回歸預(yù)測(cè)研究通常存在異常值影響大、數(shù)據(jù)變異代表性不足、區(qū)域差異性不能很好地體現(xiàn)、模型解釋影響因素之間非線性關(guān)系的程度有限等問(wèn)題，并且以往影響稻米富集Cd的研究中環(huán)境因子檢測(cè)參數(shù)少，沒(méi)有全盤考慮土壤環(huán)境因子的影響，難以篩選出影響稻米富集Cd的主控因子[9]。近年來(lái)，隨著空間自相關(guān)技術(shù)、隨機(jī)模擬技術(shù)等數(shù)據(jù)挖掘技術(shù)的發(fā)展，新技術(shù)在環(huán)境主控因子篩查中的優(yōu)勢(shì)突顯。與傳統(tǒng)相關(guān)性分析技術(shù)相比，空間自相關(guān)尤其是雙變量空間自相關(guān)將兩個(gè)變量的關(guān)聯(lián)性與空間位置耦合，可以在空間上將性質(zhì)相近或差異明顯的區(qū)域明確區(qū)分出來(lái)，便于我們更加集中地發(fā)現(xiàn)區(qū)域之間的關(guān)聯(lián)性及差異性[10]。隨機(jī)森林回歸逐漸被用于土壤污染研究，隨機(jī)森林算法訓(xùn)練速度快，能夠處理高維度數(shù)據(jù)，是確定影響因素重要度排名的重要工具之一[11]。將兩種技術(shù)結(jié)合并分別發(fā)揮各自在熱點(diǎn)區(qū)篩選、主控因子篩查方面獨(dú)特優(yōu)勢(shì)的相關(guān)研究還鮮見報(bào)道，因此本研究將雙變量空間聚類與隨機(jī)森林回歸技術(shù)相結(jié)合用于長(zhǎng)株潭地區(qū)影響稻米富集Cd土壤主控因子的篩選，從而為稻米Cd污染的修復(fù)治理及污染源解析熱點(diǎn)篩查提供基礎(chǔ)支撐。

1 材料與方法

1.1 樣品采集

本研究樣品采集自我國(guó)湖南省存在稻米Cd污染風(fēng)險(xiǎn)的某縣主要稻米種植鄉(xiāng)鎮(zhèn)，采樣時(shí)間為2016年10月下旬，采集的水稻類型為晚稻，采用格網(wǎng)方式布設(shè)采樣點(diǎn)，每個(gè)采樣點(diǎn)代表面積約為10 hm2。采集根際土和對(duì)應(yīng)的水稻籽粒樣品共2184對(duì)。拔出每株水稻后，從其根部收集粘連的土壤作為根際土，每個(gè)土壤樣點(diǎn)由5個(gè)水稻根際土子樣混合而成，土壤量為5 kg，采集深度為0～20 cm，同時(shí)對(duì)應(yīng)采集水稻樣品，每個(gè)水稻樣品由5～8株水稻組成，采樣時(shí)記錄采樣點(diǎn)經(jīng)緯度及周邊環(huán)境狀況。樣品采集、處理均按照《農(nóng)田土壤環(huán)境質(zhì)量監(jiān)測(cè)技術(shù)規(guī)范》（NY/T 395—2012）的要求。樣品帶回實(shí)驗(yàn)室風(fēng)干后，去除雜質(zhì)，研磨過(guò)10目篩；將全部土壤樣品混合均勻后采用四分法取土壤200 g，過(guò)20目篩；采用四分法取其中50 g用于土壤pH的檢測(cè)，其余150 g采用四分法縮分至50 g，過(guò)100目篩，用于土壤多參數(shù)檢測(cè)。稻米樣品采集后用自來(lái)水沖洗干凈，再用去離子水沖洗兩遍，曬干后用礱谷機(jī)去殼，然后再粉碎制成樣品，四分法過(guò)40目篩用于后續(xù)檢測(cè)。

1.2 土壤和稻米理化性質(zhì)分析

本研究除測(cè)定土壤常規(guī)理化參數(shù)外，根據(jù)研究目的還對(duì)土壤中金屬及非金屬元素含量進(jìn)行測(cè)定，合計(jì)測(cè)定參數(shù)16項(xiàng)，包括pH、SOM（有機(jī)質(zhì)）、DOC（溶解性有機(jī)碳）、Cd（鎘）、S（硫）、Cl（氯）、P（磷）、Mn（錳）、Cu（銅）、Zn（鋅）、Si（硅）、K（鉀）、Mg（鎂）、Ca（鈣）、Na（鈉）、Fe（鐵）。其中pH、SOM的測(cè)定分別采用電位法、重鉻酸鉀滴定法，方法參見相關(guān)國(guó)家或行業(yè)標(biāo)準(zhǔn)[12-13]。DOC利用TOC分析儀測(cè)定[14]。土壤中重金屬Cd全量采用三酸消解法（濃HNO3、HF和HClO4體積比為10∶4∶1），采用Agilent公司的7700x型號(hào)ICPMS測(cè)定[15]。其他元素含量采用PANalytical公司X射線能譜儀測(cè)定，儀器型號(hào)為Epsilon5，測(cè)試過(guò)程簡(jiǎn)述如下：取4 g過(guò)100目篩的土壤均勻平鋪于壓片機(jī)腔體內(nèi)，再在土壤上部均勻平鋪放置約13 g硼酸填滿腔體，放置好裝置后在10 MPa壓力下進(jìn)行壓片20 s，保證壓片厚度均一，取出壓片使用X射線能譜儀進(jìn)行測(cè)定[16]。稻米樣品Cd經(jīng)雙酸（濃HNO3和HClO4體積比為10∶1）消解，采用Agilent公司的7700x型號(hào)ICP-MS測(cè)定[17]。抽取10%土壤和稻米樣品進(jìn)行雙樣檢測(cè)，使用國(guó)家有證標(biāo)準(zhǔn)質(zhì)控樣品進(jìn)行質(zhì)量控制（GBW 07408、GBW 07423、GBW 07447），檢測(cè)結(jié)果均在標(biāo)準(zhǔn)物質(zhì)定值不確定度區(qū)間內(nèi)。

1.3 分析方法簡(jiǎn)介

1.3.1 雙變量空間自相關(guān)

Anselin等[10]提出雙變量局部空間自相關(guān)分析，該方法與一般相關(guān)性分析不同，它在計(jì)算兩個(gè)變量相關(guān)性的同時(shí)，還考慮了觀測(cè)值的地理關(guān)系。因此，用雙變量局部空間自相關(guān)分析來(lái)評(píng)估一個(gè)地理單元的某一變量與其他變量的相關(guān)關(guān)系，可以更貼切地反映變量的空間關(guān)聯(lián)。對(duì)某個(gè)空間單元i的雙變量局部空間自相關(guān)指數(shù)（Bivariate Local Moran′s i）定義如下：

式中：ωij為空間單元i、j之間的空間連接矩陣；Xik是空間單元i屬性k的值；Xjl是空間單元j屬性l的值；k、l是屬性 k、l的平均值；σk、σl是屬性 k、l的方差。若iikl顯著為正且zik大于0，則表明空間單元i和周圍鄰居的觀測(cè)值都相對(duì)較高，屬于高高集聚；若iikl顯著為正且zik小于0，則表明空間單元i和周圍鄰居的觀測(cè)值都相對(duì)較低，屬低低集聚；若iikl顯著為負(fù)且zik大于0，則表明周圍鄰居的觀測(cè)值遠(yuǎn)低于空間單元i上的值，屬高低集聚；若iikl顯著為負(fù)且zik小于0，則表明周圍鄰居的觀測(cè)值遠(yuǎn)高于空間單元i上的值，屬低高集聚[18]。

1.3.2 隨機(jī)森林

隨機(jī)森林是以決策樹為基學(xué)習(xí)器的集成學(xué)習(xí)算法，通過(guò)篩選特征集和集成思想的引入，隨機(jī)森林表現(xiàn)出比決策樹更高的性能。為了計(jì)算一個(gè)特征的重要度，需要在其他特征不發(fā)生變化的前提下，對(duì)樣本在該特征上的取值分布進(jìn)行重新排列，這意味著樣本在該特征上與預(yù)測(cè)變量的關(guān)系被破壞。然后，利用模型對(duì)“新”樣本進(jìn)行預(yù)測(cè)。兩次預(yù)測(cè)產(chǎn)生的誤差相減會(huì)出現(xiàn)3種情況，正值、負(fù)值和0，分別表示該特征對(duì)回歸預(yù)測(cè)有正影響、負(fù)影響、沒(méi)有影響。具體算法參見文獻(xiàn)[19]。

1.4 數(shù)據(jù)分析

本文以GeoDa 1.12.1.59讀入空間數(shù)據(jù)后，先構(gòu)造空間權(quán)重矩陣，然后選擇稻米Cd、土壤Cd作為研究變量進(jìn)行雙變量空間自相關(guān)分析；以R 3.5.2所載程序包Random Forest進(jìn)行隨機(jī)森林模型構(gòu)建，并基于所得模型在袋外樣本上計(jì)算各個(gè)變量的%IncMSE重要度；SPSS 22.0進(jìn)行數(shù)據(jù)描述性統(tǒng)計(jì)分析。

2 結(jié)果與討論

2.1 土壤Cd含量與稻米Cd含量空間自相關(guān)分析

本文對(duì)研究區(qū)的土壤Cd含量和稻米Cd含量?jī)蓚€(gè)變量進(jìn)行雙變量局部空間自相關(guān)分析，探討土壤Cd含量與稻米Cd含量的空間聚集特征與差異規(guī)律，具體結(jié)果如圖1和表1所示。聚類結(jié)果顯示，研究區(qū)域可劃分為5種類型（劃定雙變量空間自相關(guān)的z檢驗(yàn)值P≤0.05），分別是土壤Cd高周邊稻米Cd高（高高聚類區(qū)）、土壤Cd高周邊稻米Cd低（高低聚類區(qū)）、土壤Cd低周邊稻米Cd高（低高聚類區(qū)）、土壤Cd低周邊稻米Cd低（低低聚類區(qū)）以及兩者沒(méi)有顯著空間自相關(guān)關(guān)系的區(qū)域（不顯著區(qū)）。高低聚類區(qū)、低低聚類區(qū)主要分布在研究區(qū)的北部，這些區(qū)域無(wú)論土壤Cd含量高或者低，稻米Cd含量均較低；高高聚類區(qū)、低高聚類區(qū)主要分布在研究區(qū)的中部、南部，這些區(qū)域無(wú)論土壤Cd含量高或者低，稻米Cd含量均較高；不顯著區(qū)在研究區(qū)的大部分范圍內(nèi)呈不規(guī)律性分布。高高、高低聚類區(qū)雖然都為土壤Cd含量較高區(qū)域，但在土壤Cd含量平均值相近情況下，兩者稻米Cd富集系數(shù)差異達(dá)9倍；低低、低高聚類區(qū)雖然都為土壤Cd含量較低區(qū)域，但在土壤Cd含量平均值相近情況下，兩者稻米Cd富集系數(shù)差異達(dá)16倍。高高、低低、低高、高低、不顯著等不同聚類區(qū)土壤與對(duì)應(yīng)稻米Cd含量的相關(guān)系數(shù)分別為 0.153（P=0.008）、-0.064（P=0.238）、-0.151（P=0.002）、0.008（P=0.939）、0.034（P=0.281），大部分未達(dá)到顯著水平，可見，即使是兩者在空間上有顯著聚集效應(yīng)的區(qū)域，土壤Cd含量與稻米Cd含量的相關(guān)系數(shù)仍然較低，這說(shuō)明在不同的區(qū)域還有其他更關(guān)鍵的因素控制著稻米對(duì)Cd的吸收。因此，本研究從5個(gè)區(qū)域的2184對(duì)樣本中分別篩選出一定量的樣本，共130個(gè)樣本進(jìn)行土壤理化參數(shù)檢測(cè)，以探究其他環(huán)境因子對(duì)稻米Cd含量的影響作用。

圖1 研究區(qū)空間自相關(guān)分析圖Figure 1 Spatial autocorrelation in the study area

2.2 篩選樣品的土壤理化參數(shù)描述性統(tǒng)計(jì)分析

從5個(gè)區(qū)域中分別篩選相同數(shù)量的樣本，從而更能全面地體現(xiàn)研究區(qū)的狀況。表2展示了影響稻米Cd富集的16個(gè)土壤理化參數(shù)的描述性統(tǒng)計(jì)結(jié)果。研究區(qū)土壤pH的范圍較寬，偏度值為負(fù)及平均值小于7表明雖然大部分土壤偏酸性，但也有部分采樣區(qū)土壤pH呈現(xiàn)堿性?？梢娊?jīng)過(guò)空間聚類后篩選的樣品有較好的區(qū)域代表性，更能說(shuō)明區(qū)域的真實(shí)狀況。檢測(cè)數(shù)據(jù)統(tǒng)計(jì)結(jié)果表明，除土壤Cd外，其他重金屬如Cu、Zn、As均未超標(biāo)或僅有個(gè)別點(diǎn)位超標(biāo)，基本不會(huì)影響水稻的生長(zhǎng)發(fā)育和稻米對(duì)Cd的吸收，也基本不會(huì)影響稻米富集Cd主控因子篩查的研究。

表1 5個(gè)區(qū)域土壤Cd和稻米Cd含量的平均值和標(biāo)準(zhǔn)差Table 1 The average value and standard deviation of soil Cd and rice Cd content in five regions

表2 土壤理化參數(shù)分析結(jié)果Table 2 Soil physical and chemical parameter analysis

2.3 影響稻米富集Cd重要因素的篩選

利用空間聚類篩選的樣品數(shù)據(jù)，采用隨機(jī)森林回歸預(yù)測(cè)稻米Cd含量，模型可解釋76%的變量變異性，進(jìn)而得出稻米Cd含量影響因素重要度排名圖（圖2）。本研究區(qū)影響稻米富集Cd最主要的因素為Ca、pH、Mn，其次是Fe、Si，再次是Zn、DOC、Cl、K、P、Mg、S、Cd、Na、Cu，影響程度逐漸減少，最后是SOM影響程度最小，在16個(gè)因素中土壤Cd對(duì)稻米富集Cd的影響排名在倒數(shù)第四，這與前文的兩者相關(guān)系數(shù)分析結(jié)果一致。

土壤Ca含量（以CaO計(jì)）對(duì)稻米富集Cd的影響在所有的研究因素中排名首位，可見其對(duì)稻米富集Cd有重要影響，其作用機(jī)理一是Ca與Cd有拮抗作用[20]，二是Ca可以調(diào)控pH，進(jìn)而調(diào)控稻米對(duì)Cd的吸收。圖3為5個(gè)區(qū)域（HH：高高，LL：低低，LH：低高，HL：高低，N：不顯著）土壤Ca含量（x）與土壤pH含量（y）平均值的相關(guān)關(guān)系圖，關(guān)系模型為y=8.024 41-5.165 99×0.168 66（xR2=0.896 7）。這與Wang等[21]研究結(jié)果一致，說(shuō)明在本研究區(qū)域內(nèi)，通過(guò)調(diào)控Ca含量可顯著影響土壤pH。通過(guò)模型計(jì)算表明，按照Ca含量每增加0.05%，土壤pH變化率不高于1.5%（pH僅變化0.1）計(jì)算，此時(shí)Ca含量為0.85%，土壤pH約為6.9，可見在土壤Ca含量高于0.85%時(shí)，土壤增加Ca雖然也會(huì)使土壤pH升高，但效果不明顯。土壤Ca與土壤pH及其他各因素之間的相關(guān)關(guān)系如圖4所示，圖4表明土壤pH與土壤Ca顯著正相關(guān)，且相關(guān)系數(shù)高于其他因素。pH的改變會(huì)引起土壤一系列理化性質(zhì)的變化，進(jìn)而改變重金屬Cd的活性，影響稻米對(duì)Cd的吸收，主要作用機(jī)制為：（1）土壤pH增加會(huì)使土壤表面電荷增加，從而增加其吸附能力；（2）pH上升時(shí)重金屬在氧化物表面的被動(dòng)吸附大多轉(zhuǎn)變?yōu)閷Ｐ晕剑唬?）pH升高使土壤中的有機(jī)質(zhì)-金屬絡(luò)合物更穩(wěn)定；（4）pH升高更易于形成重金屬氫氧化物沉淀，導(dǎo)致重金屬活性下降[22]。

圖2 變量重要度（%IncMSE）排名Figure 2 Ranking of variable importance（%IncMSE）

圖3 土壤Ca與稻米Cd和土壤pH相關(guān)關(guān)系圖Figure 3 Correlation between soil Ca and rice Cd，soil pH

圖4土壤理化參數(shù)和稻米Cd含量相關(guān)性分析圖Figure 4 Correlation analysis of soil physical and chemical parameters and rice cadmium content

圖5 為土壤Fe、Mn與稻米Cd含量的相關(guān)關(guān)系圖。由圖可見，土壤Fe、Mn增加與稻米Cd含量下降呈指數(shù)關(guān)系，低濃度的土壤Fe、Mn含量的增加，會(huì)使稻米Cd快速下降，高濃度效果則不顯著，說(shuō)明土壤Fe、Mn也是影響稻米富集Cd的關(guān)鍵因素，其原因在于土壤中Fe、Mn氧化物決定了水稻根表Fe、Mn氧化物膠膜量的多少和性質(zhì)[23]，水稻根表形成的Fe、Mn膠膜對(duì)土壤中的Cd起到阻隔作用，阻擋了土壤中的Cd向水稻植物組織內(nèi)部遷移[24]。土壤溶液中的鐵錳氧化物也會(huì)與Cd離子形成膠體從而減少Cd的活性[25]。

隨機(jī)森林回歸結(jié)果表明研究區(qū)土壤Si含量也是影響稻米富集Cd的重要因素。由圖4相關(guān)性分析圖可知，Si與其他元素均呈負(fù)相關(guān)，與Fe、pH的負(fù)相關(guān)性最為明顯。由圖6可知，隨著土壤Si含量的升高稻米Cd含量上升，可能原因在于土壤Si含量反映了土壤質(zhì)地的變化，高Si含量的土壤沙性強(qiáng)、保肥能力弱、土壤Cd活性高，導(dǎo)致稻米Cd含量升高[26]。

SOM被認(rèn)為是影響土壤重金屬有效性的主要因素之一。Tang等[27]的研究表明，稻米Cd含量與SOM呈顯著負(fù)相關(guān)。但也有研究認(rèn)為[28]有機(jī)質(zhì)對(duì)稻米吸收土壤Cd沒(méi)有太大影響。本研究隨機(jī)森林回歸結(jié)果表明SOM對(duì)稻米富集Cd影響貢獻(xiàn)率最小，圖7也表明土壤SOM與稻米Cd的決定系數(shù)僅為0.052。但這不能證明土壤SOM對(duì)稻米Cd含量沒(méi)有影響，因?yàn)楸狙芯繙y(cè)定的是土壤有機(jī)質(zhì)全量，沒(méi)有對(duì)有機(jī)質(zhì)的不同組分進(jìn)行研究。何雨帆[29]的研究表明，有機(jī)質(zhì)中含有的富里酸和胡敏酸對(duì)土壤Cd活性具有不同的影響作用，而測(cè)定有機(jī)質(zhì)總量不能反映這種差異性，因此會(huì)出現(xiàn)本文的研究結(jié)果。

圖5 土壤Fe、Mn與稻米Cd相關(guān)關(guān)系圖Figure 5 Correlation between soil Fe,soil Mn and rice Cd

圖6 土壤Si與稻米Cd相關(guān)關(guān)系圖Figure 6 Correlation between soil Si and rice Cd

圖7 土壤SOM與稻米Cd相關(guān)關(guān)系圖Figure 7 Correlation between soil SOM and rice Cd

空間聚類后各個(gè)類別平均值的散點(diǎn)圖相對(duì)于原始數(shù)據(jù)散點(diǎn)圖，稻米Cd含量隨土壤環(huán)境因子升高或降低的趨勢(shì)更加明顯，說(shuō)明空間聚類方法將空間關(guān)聯(lián)性相似的數(shù)據(jù)聚成一類后使得數(shù)據(jù)變化趨勢(shì)更明顯。將全體數(shù)據(jù)集劃分為訓(xùn)練集和測(cè)試集，用訓(xùn)練集構(gòu)建隨機(jī)森林模型，用測(cè)試集進(jìn)行檢驗(yàn)得到預(yù)測(cè)值與真實(shí)值的R2為0.662，說(shuō)明經(jīng)驗(yàn)證模型的預(yù)測(cè)效果較好，由隨機(jī)森林模型得到的土壤環(huán)境因子重要度排名較可靠。

3 結(jié)論

（1）空間聚類與隨機(jī)森林回歸相結(jié)合是篩選影響稻米富集Cd土壤主控因子的有效手段，空間聚類分析為區(qū)域代表性樣點(diǎn)的采集提供了基礎(chǔ)，隨機(jī)森林回歸是鎖定影響稻米富集Cd主控因子的有效手段。

（2）空間聚類結(jié)果表明，研究區(qū)高高、高低聚類區(qū)稻米Cd富集系數(shù)差異達(dá)9倍，低低、低高聚類區(qū)稻米Cd富集系數(shù)差異達(dá)16倍。

（3）本研究區(qū)影響稻米富集Cd最主要的因素為Ca、pH、Mn，其次是Fe、Si，再次是Zn、DOC、Cl、K、P、Mg、S、Cd、Na、Cu，影響程度逐漸減小，最后是SOM影響程度最小。

（4）通過(guò)調(diào)控Ca含量可顯著影響土壤pH，進(jìn)而影響稻米對(duì)Cd的吸收。隨著土壤Fe、Mn的升高稻米Cd含量呈指數(shù)下降。土壤Si含量與土壤Fe、pH等均呈負(fù)相關(guān)，Si含量增加，稻米Cd含量相應(yīng)升高。有機(jī)質(zhì)對(duì)稻米富集Cd的影響不顯著。