基于空间聚类与随机森林的稻米富集镉影响因素筛选研究

2019-08-26郭新蕾赵玉杰刘潇威周其文王夏晖李志涛朱智伟张铁亮王祖光

农业环境科学学报 2019年8期

郭新蕾，赵玉杰，刘潇威，周其文，王夏晖，李志涛，朱智伟，张铁亮，王祖光，张璠，孙扬

（1.农业农村部农产品质量安全环境因子控制重点实验室，天津 300191；2.农业农村部环境保护科研监测所，天津 300191；3.生态环境部环境规划院，北京 100012；4.中国水稻研究所，杭州 310006）

全国土壤污染普查结果表明，我国土壤重金属污染问题突出，尤其Cd污染在局部地区较为严重[1]。土壤Cd会通过植物吸收进入食物链对人体健康造成威胁[2]。稻米是我国的主粮，相关研究表明，水稻相对小麦、玉米等农作物而言更易富集Cd[3]。因此，防控稻米Cd污染是保障我国农产品质量安全的关键环节。稻米对Cd的富集受多种因素影响，主要包括：土壤pH、土壤有机质、铁锰氧化物、土壤Ca含量，以及S、P、Cl等阴离子[4]。大田条件下稻米对Cd的吸收受多种因素的综合影响，但总体而言各因子的影响能力不同，因此锁定某区域影响稻米Cd富集的主控因素是该区域稻米Cd污染治理的关键。探查影响稻米富集Cd主控因素的方法较多，目前采用的方法主要有相关性分析、多元线性回归、偏最小二乘回归[5]、聚类分析[6]、主成分分析[7]、Cubist混合线性回归[8]等。

影响稻米富集Cd的土壤环境因子具有区域相似性及区域变异性，因此稻米对Cd的富集能力也存在明显的区域差别，而以往的研究方法很少考虑稻米Cd富集的区域差异性。且仅基于点数据进行的回归预测研究通常存在异常值影响大、数据变异代表性不足、区域差异性不能很好地体现、模型解释影响因素之间非线性关系的程度有限等问题，并且以往影响稻米富集Cd的研究中环境因子检测参数少，没有全盘考虑土壤环境因子的影响，难以筛选出影响稻米富集Cd的主控因子[9]。近年来，随着空间自相关技术、随机模拟技术等数据挖掘技术的发展，新技术在环境主控因子筛查中的优势突显。与传统相关性分析技术相比，空间自相关尤其是双变量空间自相关将两个变量的关联性与空间位置耦合，可以在空间上将性质相近或差异明显的区域明确区分出来，便于我们更加集中地发现区域之间的关联性及差异性[10]。随机森林回归逐渐被用于土壤污染研究，随机森林算法训练速度快，能够处理高维度数据，是确定影响因素重要度排名的重要工具之一[11]。将两种技术结合并分别发挥各自在热点区筛选、主控因子筛查方面独特优势的相关研究还鲜见报道，因此本研究将双变量空间聚类与随机森林回归技术相结合用于长株潭地区影响稻米富集Cd土壤主控因子的筛选，从而为稻米Cd污染的修复治理及污染源解析热点筛查提供基础支撑。

1 材料与方法

1.1 样品采集

本研究样品采集自我国湖南省存在稻米Cd污染风险的某县主要稻米种植乡镇，采样时间为2016年10月下旬，采集的水稻类型为晚稻，采用格网方式布设采样点，每个采样点代表面积约为10 hm2。采集根际土和对应的水稻籽粒样品共2184对。拔出每株水稻后，从其根部收集粘连的土壤作为根际土，每个土壤样点由5个水稻根际土子样混合而成，土壤量为5 kg，采集深度为0～20 cm，同时对应采集水稻样品，每个水稻样品由5～8株水稻组成，采样时记录采样点经纬度及周边环境状况。样品采集、处理均按照《农田土壤环境质量监测技术规范》（NY/T 395—2012）的要求。样品带回实验室风干后，去除杂质，研磨过10目筛；将全部土壤样品混合均匀后采用四分法取土壤200 g，过20目筛；采用四分法取其中50 g用于土壤pH的检测，其余150 g采用四分法缩分至50 g，过100目筛，用于土壤多参数检测。稻米样品采集后用自来水冲洗干净，再用去离子水冲洗两遍，晒干后用砻谷机去壳，然后再粉碎制成样品，四分法过40目筛用于后续检测。

1.2 土壤和稻米理化性质分析

本研究除测定土壤常规理化参数外，根据研究目的还对土壤中金属及非金属元素含量进行测定，合计测定参数16项，包括pH、SOM（有机质）、DOC（溶解性有机碳）、Cd（镉）、S（硫）、Cl（氯）、P（磷）、Mn（锰）、Cu（铜）、Zn（锌）、Si（硅）、K（钾）、Mg（镁）、Ca（钙）、Na（钠）、Fe（铁）。其中pH、SOM的测定分别采用电位法、重铬酸钾滴定法，方法参见相关国家或行业标准[12-13]。DOC利用TOC分析仪测定[14]。土壤中重金属Cd全量采用三酸消解法（浓HNO3、HF和HClO4体积比为10∶4∶1），采用Agilent公司的7700x型号ICPMS测定[15]。其他元素含量采用PANalytical公司X射线能谱仪测定，仪器型号为Epsilon5，测试过程简述如下：取4 g过100目筛的土壤均匀平铺于压片机腔体内，再在土壤上部均匀平铺放置约13 g硼酸填满腔体，放置好装置后在10 MPa压力下进行压片20 s，保证压片厚度均一，取出压片使用X射线能谱仪进行测定[16]。稻米样品Cd经双酸（浓HNO3和HClO4体积比为10∶1）消解，采用Agilent公司的7700x型号ICP-MS测定[17]。抽取10%土壤和稻米样品进行双样检测，使用国家有证标准质控样品进行质量控制（GBW 07408、GBW 07423、GBW 07447），检测结果均在标准物质定值不确定度区间内。

1.3 分析方法简介

1.3.1 双变量空间自相关

Anselin等[10]提出双变量局部空间自相关分析，该方法与一般相关性分析不同，它在计算两个变量相关性的同时，还考虑了观测值的地理关系。因此，用双变量局部空间自相关分析来评估一个地理单元的某一变量与其他变量的相关关系，可以更贴切地反映变量的空间关联。对某个空间单元i的双变量局部空间自相关指数（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和周围邻居的观测值都相对较高，属于高高集聚；若iikl显著为正且zik小于0，则表明空间单元i和周围邻居的观测值都相对较低，属低低集聚；若iikl显著为负且zik大于0，则表明周围邻居的观测值远低于空间单元i上的值，属高低集聚；若iikl显著为负且zik小于0，则表明周围邻居的观测值远高于空间单元i上的值，属低高集聚[18]。

1.3.2 随机森林

随机森林是以决策树为基学习器的集成学习算法，通过筛选特征集和集成思想的引入，随机森林表现出比决策树更高的性能。为了计算一个特征的重要度，需要在其他特征不发生变化的前提下，对样本在该特征上的取值分布进行重新排列，这意味着样本在该特征上与预测变量的关系被破坏。然后，利用模型对“新”样本进行预测。两次预测产生的误差相减会出现3种情况，正值、负值和0，分别表示该特征对回归预测有正影响、负影响、没有影响。具体算法参见文献[19]。

1.4 数据分析

本文以GeoDa 1.12.1.59读入空间数据后，先构造空间权重矩阵，然后选择稻米Cd、土壤Cd作为研究变量进行双变量空间自相关分析；以R 3.5.2所载程序包Random Forest进行随机森林模型构建，并基于所得模型在袋外样本上计算各个变量的%IncMSE重要度；SPSS 22.0进行数据描述性统计分析。

2 结果与讨论

2.1 土壤Cd含量与稻米Cd含量空间自相关分析

本文对研究区的土壤Cd含量和稻米Cd含量两个变量进行双变量局部空间自相关分析，探讨土壤Cd含量与稻米Cd含量的空间聚集特征与差异规律，具体结果如图1和表1所示。聚类结果显示，研究区域可划分为5种类型（划定双变量空间自相关的z检验值P≤0.05），分别是土壤Cd高周边稻米Cd高（高高聚类区）、土壤Cd高周边稻米Cd低（高低聚类区）、土壤Cd低周边稻米Cd高（低高聚类区）、土壤Cd低周边稻米Cd低（低低聚类区）以及两者没有显著空间自相关关系的区域（不显著区）。高低聚类区、低低聚类区主要分布在研究区的北部，这些区域无论土壤Cd含量高或者低，稻米Cd含量均较低；高高聚类区、低高聚类区主要分布在研究区的中部、南部，这些区域无论土壤Cd含量高或者低，稻米Cd含量均较高；不显著区在研究区的大部分范围内呈不规律性分布。高高、高低聚类区虽然都为土壤Cd含量较高区域，但在土壤Cd含量平均值相近情况下，两者稻米Cd富集系数差异达9倍；低低、低高聚类区虽然都为土壤Cd含量较低区域，但在土壤Cd含量平均值相近情况下，两者稻米Cd富集系数差异达16倍。高高、低低、低高、高低、不显著等不同聚类区土壤与对应稻米Cd含量的相关系数分别为 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），大部分未达到显著水平，可见，即使是两者在空间上有显著聚集效应的区域，土壤Cd含量与稻米Cd含量的相关系数仍然较低，这说明在不同的区域还有其他更关键的因素控制着稻米对Cd的吸收。因此，本研究从5个区域的2184对样本中分别筛选出一定量的样本，共130个样本进行土壤理化参数检测，以探究其他环境因子对稻米Cd含量的影响作用。

图1 研究区空间自相关分析图Figure 1 Spatial autocorrelation in the study area

2.2 筛选样品的土壤理化参数描述性统计分析

从5个区域中分别筛选相同数量的样本，从而更能全面地体现研究区的状况。表2展示了影响稻米Cd富集的16个土壤理化参数的描述性统计结果。研究区土壤pH的范围较宽，偏度值为负及平均值小于7表明虽然大部分土壤偏酸性，但也有部分采样区土壤pH呈现碱性。可见经过空间聚类后筛选的样品有较好的区域代表性，更能说明区域的真实状况。检测数据统计结果表明，除土壤Cd外，其他重金属如Cu、Zn、As均未超标或仅有个别点位超标，基本不会影响水稻的生长发育和稻米对Cd的吸收，也基本不会影响稻米富集Cd主控因子筛查的研究。

表1 5个区域土壤Cd和稻米Cd含量的平均值和标准差Table 1 The average value and standard deviation of soil Cd and rice Cd content in five regions

表2 土壤理化参数分析结果Table 2 Soil physical and chemical parameter analysis

2.3 影响稻米富集Cd重要因素的筛选

利用空间聚类筛选的样品数据，采用随机森林回归预测稻米Cd含量，模型可解释76%的变量变异性，进而得出稻米Cd含量影响因素重要度排名图（图2）。本研究区影响稻米富集Cd最主要的因素为Ca、pH、Mn，其次是Fe、Si，再次是Zn、DOC、Cl、K、P、Mg、S、Cd、Na、Cu，影响程度逐渐减少，最后是SOM影响程度最小，在16个因素中土壤Cd对稻米富集Cd的影响排名在倒数第四，这与前文的两者相关系数分析结果一致。

土壤Ca含量（以CaO计）对稻米富集Cd的影响在所有的研究因素中排名首位，可见其对稻米富集Cd有重要影响，其作用机理一是Ca与Cd有拮抗作用[20]，二是Ca可以调控pH，进而调控稻米对Cd的吸收。图3为5个区域（HH：高高，LL：低低，LH：低高，HL：高低，N：不显著）土壤Ca含量（x）与土壤pH含量（y）平均值的相关关系图，关系模型为y=8.024 41-5.165 99×0.168 66（xR2=0.896 7）。这与Wang等[21]研究结果一致，说明在本研究区域内，通过调控Ca含量可显著影响土壤pH。通过模型计算表明，按照Ca含量每增加0.05%，土壤pH变化率不高于1.5%（pH仅变化0.1）计算，此时Ca含量为0.85%，土壤pH约为6.9，可见在土壤Ca含量高于0.85%时，土壤增加Ca虽然也会使土壤pH升高，但效果不明显。土壤Ca与土壤pH及其他各因素之间的相关关系如图4所示，图4表明土壤pH与土壤Ca显著正相关，且相关系数高于其他因素。pH的改变会引起土壤一系列理化性质的变化，进而改变重金属Cd的活性，影响稻米对Cd的吸收，主要作用机制为：（1）土壤pH增加会使土壤表面电荷增加，从而增加其吸附能力；（2）pH上升时重金属在氧化物表面的被动吸附大多转变为专性吸附；（3）pH升高使土壤中的有机质-金属络合物更稳定；（4）pH升高更易于形成重金属氢氧化物沉淀，导致重金属活性下降[22]。

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

图3 土壤Ca与稻米Cd和土壤pH相关关系图Figure 3 Correlation between soil Ca and rice Cd，soil pH

图4土壤理化参数和稻米Cd含量相关性分析图Figure 4 Correlation analysis of soil physical and chemical parameters and rice cadmium content

图5 为土壤Fe、Mn与稻米Cd含量的相关关系图。由图可见，土壤Fe、Mn增加与稻米Cd含量下降呈指数关系，低浓度的土壤Fe、Mn含量的增加，会使稻米Cd快速下降，高浓度效果则不显著，说明土壤Fe、Mn也是影响稻米富集Cd的关键因素，其原因在于土壤中Fe、Mn氧化物决定了水稻根表Fe、Mn氧化物胶膜量的多少和性质[23]，水稻根表形成的Fe、Mn胶膜对土壤中的Cd起到阻隔作用，阻挡了土壤中的Cd向水稻植物组织内部迁移[24]。土壤溶液中的铁锰氧化物也会与Cd离子形成胶体从而减少Cd的活性[25]。

随机森林回归结果表明研究区土壤Si含量也是影响稻米富集Cd的重要因素。由图4相关性分析图可知，Si与其他元素均呈负相关，与Fe、pH的负相关性最为明显。由图6可知，随着土壤Si含量的升高稻米Cd含量上升，可能原因在于土壤Si含量反映了土壤质地的变化，高Si含量的土壤沙性强、保肥能力弱、土壤Cd活性高，导致稻米Cd含量升高[26]。

SOM被认为是影响土壤重金属有效性的主要因素之一。Tang等[27]的研究表明，稻米Cd含量与SOM呈显著负相关。但也有研究认为[28]有机质对稻米吸收土壤Cd没有太大影响。本研究随机森林回归结果表明SOM对稻米富集Cd影响贡献率最小，图7也表明土壤SOM与稻米Cd的决定系数仅为0.052。但这不能证明土壤SOM对稻米Cd含量没有影响，因为本研究测定的是土壤有机质全量，没有对有机质的不同组分进行研究。何雨帆[29]的研究表明，有机质中含有的富里酸和胡敏酸对土壤Cd活性具有不同的影响作用，而测定有机质总量不能反映这种差异性，因此会出现本文的研究结果。

图5 土壤Fe、Mn与稻米Cd相关关系图Figure 5 Correlation between soil Fe,soil Mn and rice Cd

图6 土壤Si与稻米Cd相关关系图Figure 6 Correlation between soil Si and rice Cd

图7 土壤SOM与稻米Cd相关关系图Figure 7 Correlation between soil SOM and rice Cd

空间聚类后各个类别平均值的散点图相对于原始数据散点图，稻米Cd含量随土壤环境因子升高或降低的趋势更加明显，说明空间聚类方法将空间关联性相似的数据聚成一类后使得数据变化趋势更明显。将全体数据集划分为训练集和测试集，用训练集构建随机森林模型，用测试集进行检验得到预测值与真实值的R2为0.662，说明经验证模型的预测效果较好，由随机森林模型得到的土壤环境因子重要度排名较可靠。