基于机器学习识别中微量元素与我国七大片区产地稻米镉、砷的富集规律

2023-11-09牟力言刘春湘陈敏秦莉林大松BatsaikhanBayartungalag

农业环境科学学报 2023年10期

牟力言，刘春湘，陈敏，秦莉*，林大松，Batsaikhan Bayartungalag

（1.农业农村部环境保护科研监测所，天津 300191；2.农业农村部环境保护科研监测所湘潭综合实验站，湖南湘潭 411100；3.蒙古科学院地理与生态地质研究所，乌兰巴托 1568683）

土壤重金属污染是影响当下我国经济发展和生态安全的首要问题之一[1-2]。2014 年《全国土壤调查公报》指出，国内土壤污染超标率为16.1%，无机型污染类别占比82.8%；其中，Cd、As超标率远高于其他无机污染物，分别占比7%、2.7%，并在农业型土壤中更高。由土壤-稻米作物系统呈现的超标Cd、As迁移转化规律及其影响推动人们对污染效应与健康关系问题的认识[3-4]，尤其是食物链传递形式转移带来的Cd、As 污染问题广泛存在[5-6]。大量的研究表明，土壤理化性质[7]、土壤类型[8]、土地利用方式[9]、作物种类及其农艺措施[10]、土壤重金属形态及含量[11]等都是影响稻米Cd、As等重金属污染效应的关键因子。因此，不得不考虑以主控因子搭建模型来预测Cd、As 在作物中的传递，如将土壤pH 和SOM 作为胡萝卜（Daucus carotaL.）Cd 吸收主控因子构建的土壤Cd 预测模型[12]。中微量元素（Ca、Mg、Cu、Fe、Mn 和Zn 等）由于其作为水稻生长过程中必需的有益营养元素而凸显其重要性。这些元素的缺乏和过量都有可能对稻米Cd、As 的富集产生促进或抑制，包括Ca、Zn、Si、Se 等元素带来的影响[13-17]。然而，这些文献并未提供有关中微量元素与稻米对Cd、As 富集关联的生物有效性或毒性的充足证据。因此，利用中微量元素来探索与稻米Cd、As富集的浓度分配关系是必要的，尤其对在大尺度地域水平上稻米Cd、As 富集差异之间的联系具有现实意义。

机器学习（ML）算法，如决策树算法（CART）[18]、随机森林（RF）[19]和人工神经网络（ANN）[20]等，已成为揭示多元因子间隐藏关系的强大工具。它们在污染物毒性风险预测[21]、污染修复预测[22]和材料合成设计[23]等问题上引起了广泛关注并得到认可。该算法可以通过中微量元素输入学习和模拟来预测稻米Cd、As 富集是否超标，进而探讨稻米Cd、As 污染与中微量元素贡献之间的空间错配关系。同时，利用该模型可以预测土壤中的重金属吸附，并绘制基于土壤重金属吸附能力的全球分布[24]；也可通过植被、海拔、土壤质地和气候等表观数据计算出生态区域重金属浓度，用以评估区域内重金属的生物积累[25]。尽管它们被认为是黑箱模型，但该模型的可解释性使我们能够同时分析多种因素对于目标的贡献，并探究其联系[22]。决策树和随机森林均属于基于信息论的分类器，十分适合离散特征的处理。而对于离散特征，其他很多模型都需要对其进行编码，得到很稀疏的编码向量再进行模型拟合。并且，决策树和随机森林实现起来较为简单，易于理解和可视化，规则也易于表达，在当前的相关研究中应用广泛，方法成熟。考虑到环境的多元复杂性，基于机器学习可提取土壤-稻米体系Cd、As污染空间关系的特定识别规则并绘制结果，可以为重金属污染的准确风险评估和有针对性的防治措施提供有价值的参考。

目前基于机器学习手段的中微量元素识别稻米Cd、As 富集规律研究较少。为了探究全国不同区域以中微量元素预测的稻米Cd、As富集及其规律，做出以下假设：中微量元素（Ca、Mg、Cu、Fe、Mn、Zn）与稻米Cd、稻米无机As 的富集存在关联，并对其积累过程产生正向贡献。为验证该假设，本研究基于全国点对点监测数据和多采样点数据，通过建立和优化CART、RF和土壤-稻米生物有效性模型，以期为区域稻米Cd、As 含量预测及土壤Cd、As 的生态安全阈值的国家标准制定提供可靠的数据参考。

1 材料与方法

1.1 数据来源与预处理

本研究数据来源于农业农村部环境监测总站对全国25个省份的土壤和稻米的2020—2021年间的例行监测数据，以及湖南湘潭、云南禄丰、江西新余、浙江大田试验监测站的采样数据。数据范围涵盖全国七大片区（华东、华北、华中、华南、西南、西北和东北片区）的土壤和稻米Cd、As 数据；被考察的中微量元素有土壤交换性Ca、土壤交换性镁（Mg）、土壤有效铜（Cu）、土壤有效Fe、土壤有效锰（Mn）、土壤有效Zn；土壤性质参数有：pH、SOM、阳离子交换量（CEC）。共9 515组数据，用于CART的预测模型训练、RF对稻米Cd、As 主控因子筛选、以及土壤-稻米体系有效性模型分析，具体的数据筛选原则如下：①数据必须来自于土壤和稻米的点对点协同监测；②数据必须有本研究考察的中微量元素指标和土壤性质参数；③选择的数据在每个分区覆盖的样本量足够，具体为：华东2 650 组、华北1 462 组、华中1 310 组、华南1 365 组、西南1 420组、西北488组、东北820组。

1.2 指标测定

本研究中所涉及的土壤理化性质和中微量元素，以及土壤稻米Cd、As 的测定方法参照国家现行标准执行。指标及其对应的标准编号如下：

土壤pH：NY/T 1121.2—2006；SOM：NY/T 1121.6—2006；CEC：LY/T 1243—1999；土壤总Cd：HJ 766—2015；稻米Cd：GB 5009.268—2016；稻米As：GB 5009.268—2016；土壤总As：GB/T 22105.2—2008；土壤交换性Ca：NY/T 1121.13—2006；土壤交换性Mg：NY/T 1121.13—2006；土壤有效铜、铁、锰、锌：HJ 804—2016。

1.3 模型设计与优化

1.3.1 土壤-稻米体系中稻米Cd、As 迁移的影响因素确定

先前的研究中，关于土壤-稻米体系中Cd、As 迁移的影响因素在土壤理化性质、金属形态及含量、稻米基因型以及其他农艺措施等方面进行了更多的探讨[26-28]。另外，土壤pH、SOM 和CEC 在产地稻米Cd、As 富集过程中的重要性影响在我们先前的研究中也进行了研究[29]。目前，中微量元素在土壤和作物中的分布迁移被广泛关注[30-31]，尤其是其他中微量元素对Cd、As 富集可能带来的影响[32]。通过文献调查，我们做出进一步假设：中微量元素（Ca、Mg、Cu、Fe、Mn、Zn）与稻米Cd、As 的富集存在一定的关联，并对其积累过程产生正向贡献。因此，有必要构建模型来预测和验证中微量元素在识别稻米Cd、As 超标富集规律中的重要性。

1.3.2 决策树模型

CART是一种分类与回归模型[33]。它是一棵二叉树，即每个节点下面包含两个子节点。在CART 的生成过程中，根节点包含所有样本，然后按照分裂准则，根节点又被分为两个子节点，重复执行这个过程，直至节点不可再分为止。即基于特征对实例的分类，通过样本学习，最终达到关联关系预测的效果[18]。根据先前研究中CART构建步骤，其通常包括3个步骤：特征选择、决策树的生成和决策树的修剪。

1.3.3 随机森林

RF 是通过引入随机属性加入到CART 算法，从而构建bagging 与CART 相结合的模型算法[34]，并通过多颗决策树中每颗树的投票结果获取得到最优结果[35]。为了确定全国七大片区稻米产地土壤Cd、As污染的主控因子，我们建立了RF模型。具体地，通过调用RF 算法，设置最大特征数、最大深度、CART 树数量等参数；然后输入最优的mtry 和ntree，将中微量元素因子指标（Ca、Mg、Cu、Fe、Mn、Zn）与土壤理化指标（pH、CEC、SOM）作为输入变量；将重金属Cd、As的BCF（生物富集系数）作为输出变量；建立回归模型并对各因子实现重要性分析。具体的计算公式如下：

式中：m，n，t分别是基础指标因子总数、分类树和单颗树节点数；DGkhj为第k个指标因子在第h颗树的第j个节点上基尼指数减少值；Pk为第k个基础指标因子的重要性。

RMSE（Root Mean Square Error）均方根误差和回归系数（R2）用以评价模型准确性的指标，即以预测值与真实值偏差的平方与观测次数n比值的平方根。本研究以此衡量预测值与真实值之间的偏差，并对异常值作出敏感性判断[35]。

式中：n表示观测次数，Yi表示预测值，f（xi）表示真实值。RMSE取值[0，+∞），当预测值与真实值的误差为0 时，RMSE为0，即完美模型；误差越大，RMSE越大。实际过程中，RMSE值在0.2～0.5 之间，即可认为预测效果较好。

R2表征依变数Y的变异中有多少百分比可由控制的自变数X来解释[35]。其意义在于判断自变量对因变量的解释程度，即拟合优度，其拟合优度越大，自变量对因变量的解释程度就越高，自变量引起的变动占总变动的百分比就越高。其取值范围通常为[0，1]。

1.3.4 生物有效性模型

BCF 通常以污染物在生物体内浓度与其生存土壤环境中浓度的比值计算得到[29]，公式如下：

式中：Crice为稻米中Cd/As 含量（mg·kg-1）；Csoil为土壤中的总Cd/As含量（mg·kg-1）。

通过多元线性回归分析，构建生物有效性模型，其基本形式如下：

式中：BCF为生物富集系数；pH为土壤pH；SOM为土壤有机质，g·kg-1；CEC为土壤阳离子交换量cmol·kg-1；Ca、Mg、Cu、Fe、Mn和Zn分别表示土壤交换性Ca、土壤交换性Mg、土壤有效Cu、土壤有效Fe、土壤有效Mn和土壤有效Zn含量，mg·kg-1；a～i均为无量纲参数，表示土壤理化性质对生物富集系数的影响程度；k为方程的截距，表示稻米富集重金属Cd的固有敏感性。

1.4 统计分析

本研究中采用Excel 对数据进行整理分析；采用SPSS 22.0 进行中微量元素、土壤理化性质与稻米Cd/As 的BCF进行多元回归分析；使用Python 3.7 进行CART 与RF 模型的归一化、标准化处理，并实现模型的构建与优化。

2 结果与讨论

2.1 模型的优化与评估

我们延伸大量文献集中pH、SOM 和CEC 等土壤理化性质对稻米Cd、As超标影响贡献的研究思路，尝试建立了少有研究的中微量元素对稻米Cd、As 超标的影响预测模型。中微量元素驱动的稻米Cd、As 的吸附、溶解和共沉淀等过程可能是稻米Cd、As加速富集的原因之一[32]。参照CART原理[33，36]，将6种中微量元素（Ca、Mg、Cu、Fe、Mn 和Zn 元素丰缺指数结果）为特征输入（X），Cd 和As 的超标情况为结果输出（Y），带入到CART算法中进行训练。随机选择70%的原始数据作为训练样本集，随机选取30%的原始数据作为测试样本集，进行100次交叉循环验证，利用训练样本集训练CART模型，当测试样本集的准确率最高时，输出该CART 树。CART 树构建完成后以预测样本遍历所构建的CART 模型，将模型预测结果与样本实际结果进行比较，从而得到决策树预测精度。

依据《土壤环境质量农用地土壤污染风险管控标准（试行）》GB 15618—2018 和《食品中污染物限量》GB 2762—2022，并结合实际情况以样本Cd 含量超过0.2 mg·kg-1、As 含量超过0.35 mg·kg-1被认定为超标。通过训练CART 模型，Cd、As 超标模型的预测准确率最低分别为82.13%、85.57%，最高分别为95.55%、97.55%；其平均准确率分别为88.10%、90.34%。足够高的预测准确率表明中微量元素是影响稻米Cd、As超标情况的重要判别指标，也表明土壤中微量元素（Ca、Mg、Cu、Fe、Mn 和Zn）与稻米Cd、As富集有着显著的相关性。

通过进一步量化，我们评估了基于RF 的全国七大片区稻米Cd、As污染主控因子分析模型的准确性。如表1 所示，不同片区的Cd 和As 模型对应的RMSE分别在0.11～0.44和0.18～0.48之间；不同片区的Cd和As 模型对应的R2分别在0.60～0.86 和0.54～0.83 之间。表明全国七大片区稻米Cd、As 污染主控因子的分析模型具有较高的精确度。这同时支持了我们对全国七大片区稻米Cd、As 超标富集主控因子贡献率的分析结果。

表1 全国七大片区稻米产地土壤重金属Cd、As污染主控因子分析模型评估指标Table 1 Assessment parameters of the main control factor analysis model for heavy metal Cd and As pollution in the soils of rice production areas in seven major regions in China

2.2 稻米镉的主控因子分析

为了探究全国七大片区Cd富集究竟由哪种环境因子主导，我们进行了不同区域和不同主控因子对Cd富集的贡献率比对分析（图1）。整体上，不同区域主控的富集驱动因子存在明显差异，以及同一主控因子在不同区域提供的稻米Cd 富集之间存在明显差异。由单一因子主要驱动的Cd富集在不同区域上具体表现为：华东片区pH 的贡献占主导、华南片区的Ca 元素和东北片区的SOM 分别占主导贡献；以及由两个或三个因子主要驱动的稻米Cd富集表现在不同区域上呈现为：华北片区的Fe 和pH 贡献占主导，西北片区的Mn 和CEC 贡献占主导，华中片区的Fe、Zn和CEC 贡献占主导；对于西南片区，其不同因子之间的贡献水平表现出较小差别的梯度减小贡献，依次为Zn＞CEC＞Mn＞pH＞Fe＞SOM＞Ca＞Cu＞Mg，但均表现出对稻米Cd 富集影响的相当贡献（图1a）。同样地，同一主控因子在不同区域的贡献率也表现出明显差异（图1b）。这种由地域差异引起的环境因子驱动的Cd富集结果可见，Ca元素在华南片区、Fe元素在华北和华南片区、SOM 在东北片区以及pH 在华东和华北片区的贡献相比其他区域明显较高。

图1 全国七大地理片区Cd主控因子贡献率Figure 1 Contribution rate of Cd main control factors in China′s seven major geographical regions

区域水平上考察不同环境因子对稻米Cd富集是识别主控因子准确性的重要支持。先前的研究表明，土壤Cd 含量与土壤pH、土壤交换性Ca、土壤有效Zn呈现出显著的相关性关系，并且对耕作手段、地域特征等因素带来的土壤理化性质、中微量元素对稻米Cd 富集可能的影响做出推断，如CEC、Zn 和Ca 等[37]。和君强等[38]通过模型构建了湖南长沙地区稻米Cd 污染土壤中pH和SOM与Cd富集之间的显著关联；基于随机森林回归的主控因子识别模型，其Ca、pH、Mn 等也是影响湖南某县稻米Cd 富集的主要影响因素[39]。东北片区特有的土壤类型中丰富的SOM 对稻米Cd富集的贡献显而易见[40]。受土壤质地和类型的影响，华北和华中片区中土壤Fe 可能通过其有效性形态的变化来改变水稻根际Fe 膜对Cd 的吸附固定量，进而影响对稻米Cd的富集[41]。同时，土壤中Ca对土壤pH的调控作用，可能是华南片区稻米Cd 过度富集的原因之一。总之，比较分析发现全国七大片区稻米Cd富集的主控因子存在区域差异，也存在主控因子贡献率上的差异。这种差异基于区域元素化学循环，也可能与中微量元素和土壤中Cd 之间的拮抗作用有关，也可以是基于对土壤理化性质如pH的调控进而对稻米Cd的富集产生调控[39，42-43]。

综上，9个环境因子在七大地理片区中对稻米Cd的富集贡献和影响存在明显的地域差异，这也说明影响稻米Cd富集的土壤环境因子的地域差异性。中微量元素指标因子对各大地理片区稻米Cd的富集影响显著，除了与地域母质，区域施肥、耕作制度差异等相关外，可能还与不同地理片区的特异的大气沉降条件等有关[42]。因此，在华东片区，应重点监控pH、SOM等土壤理化指标的变化情况；华北和华中片区应重点监控Fe 指标与土壤pH 的变化情况；华南片区应重点监控Ca 指标、pH、CEC 等土壤理化指标；而西南、西北、东北等片区应重点监控Zn、Mn 元素指标、CEC、SOM等土壤理化指标。

2.3 稻米砷的主控因子分析

与Cd 的分析过程类似。通过比较发现，全国七大片区稻米As主控富集驱动因子在区域尺度上同样存在一定差异，同区域不同因子的贡献率明显不同（图2）。由单一因子主要驱动的As 富集在不同区域上具体表现为：Fe 元素在华东、华南和西南片区的贡献相比于其他因子尤为突出，Zn 在华北片区的因子贡献中占据主导地位，包括西北片区的Mn 和东北片区的Mg；另外，华中片区的Cu、Zn和Fe同时具有较高的贡献。与Cd 不同的是，在因子水平上Fe 表现出对稻米As 富集明显的特异性贡献，除了华北片区。而其他因子在不同区域也表现出相当的贡献水平。

图2 全国七大地理片区As主控因子贡献率Figure 2 Contribution rate of As main control factors in China′s seven major geographical regions

土壤Fe 含量和形态依赖于稻米土壤淹水和排水条件而变化[44]，进而对土壤有效态As 的含量和形态产生调节作用。主要表现为Fe氧化物和As（Ⅴ）的协同还原过程，加速溶解产生的As（Ⅲ）造成稻米As 的吸收风险[44-46]。这解释了多片区Fe 元素对稻米As 富集的主要贡献，尤其是土壤Fe 含量较高、雨水较广的华东、华南和西南片区等地。另外，Zn 元素可促进水稻土壤中As 的甲基化和脱甲基化过程，并影响水稻根系对As的吸收和运输；Mg元素与土壤中As（Ⅲ）或有机甲基As 的协同作用可加速稻米As 富集并影响其在水稻不同部位的分配；而Mn 氧化物在土壤中的吸附与As 有关[47-49]。总之，不同片区稻米As 富集可能依赖于这些中微量元素在土壤中的形态过程，从而表现出明显差异。这可能是华北、西北和东北地区的Zn、Mn、Mg主控稻米As富集的潜在原因。

不同区域的土壤中As主控因子贡献率存在明显差异，其中Fe、Zn、Mn、Mg 等元素对稻米中As 的富集影响最为显著。Fe 元素在华东、华南、西南、东北片区对稻米中As 的富集影响显著，推断Fe 元素是影响我国水稻产区As 富集的重要因素，这与前人的研究结论一致。As 污染的水稻土壤中Fe 与As 是普遍存在的，并且Fe 的一系列氧化还原反应是导致As 迁移的主要原因[50]。Zhao 等[6]建议旱作、干湿交替可减少水稻As 吸收和在籽粒中的积累，但是会增加重金属Cd 的吸收。这进一步说明土壤质地、金属形态和管理模式等对稻米富集无机金属As的影响。有效的人为监测管控措施可能极大减小稻米对As的富集。例如，对于华东、华南和西南片区，应加强监测土壤Fe含量和形态，及时排水以减少Fe 氧化物还原释放五价无机As；对于华北和华中片区，应加强监测土壤Zn含量，并施用适量Zn 肥以保证足够Zn 供给；同时应施用适量生物质肥料以增加土壤微生物活性并提高土壤pH值。

2.4 土壤-稻米体系生物有效性模型

基于随机森林算法分析得到全国七大地理片区Cd、As 主控因子贡献率的统计结果，进一步提取了Cd、As污染基于上述指标空间关系的特定识别规则，并构建了Cd、As 的土壤-稻米体系生物有效性模型（表2、表3）。量化结果显示，9 个因子与稻米Cd、As的BCF呈显著相关（P＜0.05）。我们逐步引入上述因子考察基于土壤-稻米体系生物有效性模型的解释能力。在不考虑区域差异与重金属本身特性的情况下，引入基于pH、SOM 和CEC 的三因子量化模型，不同片区Cd 对应的决定系数在0.303～0.400 之间（表2，P＜0.05），不同片区As对应的决定系数在0.300～0.340之间（表3，P＜0.05）。在考虑区域差异与重金属本身特性的情况下，对所考察因子进行重要性排序，构建基于Top5 的五因子量化模型，不同片区Cd 对应的决定系数在0.450～0.583 之间（表2，P＜0.05），不同片区As 对应的决定系数在0.450～0.634 之间（表3，P＜0.01），表明五因子量化的结果能很好地解释其对稻米Cd、As富集的影响。

表3 七大片区As土壤-稻米体系联合量化预测模型Table 3 Joint quantitative prediction equation of As soil-crop system in seven large areas

为了进一步得出更优结果，我们将以上因子全部引入模型，结果不同片区Cd 和As 对应的决定系数最大可达0.680和0.664（表2，表3，P＜0.05），较高的确定系数表明九因子模型能显著解释pH、SOM、CEC、Ca、Mg、Cu、Fe、Mn 和Zn 对稻米Cd、As 富集的影响。然而，环境差异会引起土壤理化性质、微生物群落以及稻米自身的一系列变化，进而会导致模型决定系数的降低[29]，具体表现在本研究上的区域空间跨度大、土壤类型和水稻品种的显著差异，以及田间环境中其他土壤理化因子对稻米富集Cd、As的影响，如氧化还原电位、电导率、土壤机械组成等。在实际生产过程中，因子数量带来的资源成本问题，同时可以选择基于Top5因子模型结果来做一些综合评估。

受地域差异、土壤背景、农业过程等综合因素的影响，九因子对稻米中重金属Cd、As富集的贡献率呈现出差异，这种差异反映了由环境复杂性导致的土壤和稻米中重金属污染的空间分布规律。基于土壤-稻米的生物有效性模型，进一步确立了本模型对大尺度范围筛查土壤和水稻中Cd、As 超标识别的重要性地位，以及中微量元素对土壤和稻米中Cd、As富集的影响。尤其是如Fe 对As 富集在华东、华南和西南等片区的区域型特异性影响，应该进一步被广泛关注，并采取有效的管控措施。