农产品产地环境安全评估：重金属转运模型研究

如上节所述，农作物吸收土壤重金属受多种因素影响，是多因子综合作用的结果，除土壤因素外，还受植物自身基因型及环境因素的作用，两者之间的关系非常复杂。长期以来，科研人员试图找到一条捷径或称桥梁来打通土壤-作物系统重金属转运特征预测或风险评估问题。近年来，虽多有突破，但距离建立稳定可靠的土壤-作物系统重金属转运量化模型还有很多工作要做。

土壤-作物系统重金属转运模型类型较多，大致可分为两类，一类是机理模型，二是经验模型。机理模型一般以重金属形态分析模型为主，通过解析土壤溶液不同形态重金属的组成及含量来计算植物可利用重金属含量进而建立土壤-植物系统重金属转运关系模型（田娜娜等，2017），如Mitserlich函数模型。机理模型往往比较复杂，所需参数多，且模型预测精度也不尽如人意，原因在于模型计算的重金属植物可利用态是否真的是植物可利用态还处于“灰箱”状态。第二类模型是经验模型，这类模型的形式多样，有简有繁。主要包括以下两种：

第一种是简单的系数模型。即通过简单的计算植物对重金属的吸收系数（吸收系数＝农作物重金属含量/土壤重金属含量）来反映植物对重金属的吸收性能。吸收系数也是美国国家环境保护局农用地风险评估模型中采用的土壤-作物重金属运移模型，但由于吸收系数的多变性及区域差异性，因此，吸收系数仅能用于风险评估第一层次即简单的风险评估中，对于更加准确的风险评估工作则需要建立更加精确的模型。

第二种是经验模型为回归模型。即采用线性或可转化为线性回归的非线性回归技术（如弗罗因德利希模型、平稳模型）构建土壤-作物重金属转运模型。模型构建常用的参数即自变量为影响农作物对重金属吸收的环境因子，如重金属全量、pH、SOM、CEC、Fed等，也有采用有效态重金属含量作为自变量的。下面就两种回归技术做一下简要使用说明。

一种回归技术为多元线性回归，模型形式为Y＝a0＋a1x1＋a2x2＋a3x3＋…＋ε，模型中a0为常数项，ax为影响农作物对重金属吸收的环境因子系数，xx为自变量，即影响农作物吸收重金属的环境因子，如pH、SOM、CEC、Fed等，因影响因子的复杂性，x可以是一次项，也可以是平方项或是经对数变化的数值，Y为农作物重金属含量，可以是未经变化的检测值，也可以是经对数变化的数值。这要看模型拟合的结果及原始数据的分布情况。采用逐步线性回归方法建立土壤-作物重金属转运模型，模型构建前需先对数据进行一定的分析，如数据的分布情况、异常值情况、数据之间的相关性等，以避免数据无法满足多元线性回归对参数正态性假设的要求，减少异常值的影响及避免数据的共线性问题。模型构建过程中，采用决定系数R2及方差分析方法检验构建模型的显著性，采用特征值和条件指数结合方法分析自变量的共线性，判断依据为自变量特征值越接近于0，共线性越明显；条件指数大于15说明可能存在共线性问题，当指数大于30时，说明有严重的共线性问题。采用中心化杠杆值分析统计量的异常值，采用Cook距离、DfFits和DfBetas来查找强影响点。一般认为，当杠杆值＜0.2时，对应个案的值是正常的，当0.2＜杠杆值＜0.5时，对应个案的值就比较危险了，而杠杆值＞0.5的情况是不允许出现的；一般可采用制作DfFits和DfBetas序列图的方式分析强影响点（李野等，2012）。近年来，采用多元线性回归技术构建土壤-作物重金属转运模型的研究也较为常见。如王梦梦等（2018）采用对数模型，将土壤全Cd含量、pH、SOM、CEC同时考虑对稻米Cd含量进行多元回归分析，预测当地水稻籽粒中Cd的含量。潘杨等（2015）采用多元回归分析的方法构建土壤-作物富集Cd的关系模型，模型引入了土壤Cd全量、土壤pH、土壤有机质含量3个参数。

采用多元线性回归构建土壤-稻米Cd转运模型时，由于自变量与因变量往往存在非线性关系，从而导致参数的引入与剔除存在误差，无法精确考虑在不同条件下不同土壤因素对稻米吸收Cd的影响能力的差别，最终影响模型拟合结果。在此，推荐一种混合线性回归的方法来构建土壤-作物重金属转运模型，模型的英文名称为Cubist。

Cubist模型是一种强有力的基于规则的预测模型，每个规则都关联一个多元线性回归子模型，规则与模型的匹配完善了单一模型带来的不足，从而提升了模型预测精度。近年来，Cubist模型广泛应用于遥感影像识别分类、数字土壤制图（Gray et al.，2015），将本方法应用于土壤-作物重金属转运模型构建的研究还较少。Cubist模型树由Rule Quest公司开发，是一种在叶子节点采用线性回归函数，表示一种分段式多元线性函数通过一系列的独立变量来预测一个变量的值。Cubist模型由分离数据的预测变量的一个或多个嵌套if-then语句组成。在这些分区中，使用线性回归模型来预测结果。Cubist方法非常适合于自变量对因变量的影响在不同条件下权重不同的情况，采用Cubist结合多自变量构建土壤-作物重金属转运模型，其精度较采用单一多元回归技术要提升几倍，更有甚者，采用多元回归无法构建的方法，在采用Cubist技术后会得到意想不到的结果。在R语言里有Cubist程序包可以免费使用。其使用命令格式为“mod1＜－cubist（x＝databasename[，－x]，y＝databasename＄x）”，模型预测结果展示命令为summary（mod1）。其预测命令为predict。例如，数据库为rice，因变量为数据库的第15列稻米Cd含量，列名称为R_Cd，由模型构建格式为“Rice_Cd＜－cubist（x＝rice[，－15]，y＝rice＄R_Cd）”，模型结果展示为summary（Rice_Cd）。表8-4和表8-5是采用Cubist方法构建的稻米富集Cd模型的形式（刘佳凤等，2018）。由图8-9可知，该模型预测的稻米Cd含量与实测稻米Cd含量相关性达到了0.87。

表8-4　Cubist预测稻米富集Cd模型

注：表中S和Cl的单位为mg/kg，Fe2O3和CaO的单位为％。

表8-5　不同因子在Cubist模型中的条件贡献及模型贡献（％）

图8-9　Cubist模型预测稻米Cd含量与实际检测值散点图

注：A图为培训数据，B图为验证数据。

表8-4表明，本研究构建的土壤-稻米Cd转运模型共有7个规则，其中规则1及规则3、4所含样品量最多，占全部培训数据集样品量的81.4％。在规则1中，以土壤pH为判断条件，当pH＞6.97时，即土壤处于碱性水平时，稻米Cd含量受pH及土壤Cd含量的影响。pH升高，稻米Cd含量下降；土壤Cd含量升高，稻米Cd含量随之升高。在规则1条件下，稻米Cd平均含量为0.04mg/kg。假定土壤pH固定为7，则稻米Cd含量达到食品安全标准0.2mg/kg，土壤Cd含量需达到12mg/kg，这在现实中基本不存在，可见在湖南碱性土壤种植的水稻，基本不存在超标情况。

在规则3中展示的是当土壤pH处于中性水平时，稻米Cd含量受pH及Fe含量影响的模型形式。根据规则3，稻米对Cd的吸收受土壤pH、Fe2O3含量的影响，两者影响力一负一正。对比规则1与规则3中pH的系数，可以看出，在中性土壤中，pH变化对稻米富集Cd的影响力是碱性土壤的10倍。

图8-10　土壤有机质与土壤S含量散点图

规则5、6、7含有18个稻米监测样品，虽然样本代表性不十分充分，但3种规则条件下稻米均处于严重超标水平，应引起重视。综合3个规则表明，在强酸性土壤（pH＜5.7）上稻米Cd污染严重，pH增加更会导致稻米Cd含量的上升，这在土壤pH控制时应引起注意。土壤低Fe水平（Fe2O3≤4.9％）增加Fe含量会抑制稻米对Cd的富集，但高Fe时（Fe2O3＞5.8％），会促进稻米对Cd的吸收，规则2也表明了这一点。(www.xing528.com)

综合本研究各规则中pH对稻米富集Cd的影响可见，其过程复杂，低pH及高pH均能抑制水稻对Cd的富集，但当土壤pH由强酸性向偏酸性上升时，稻米对Cd的富集增加，偏酸向中性过度时，稻米对Cd的富集显著下降，中性向碱性转化时稻米对Cd的富集也随之下降，但不如偏酸向中性过度明显。

作者收集了我国南方7个水稻种植区不同pH条件下稻米对Cd吸收系数的变化情况，（吸收系数＝稻米Cd含量/土壤Cd含量），结果如图8-11所示。由此可见，不同地区稻米对Cd的吸收系数虽有较大差异性，但pH对吸收系数变化的影响规律是一致的。这与pH影响土壤Cd存在的形态变化有直接关系，当pH变化时，土壤Cd与土壤有机质、黏土矿物及Fe/Mn氧化物结合的形态会产生显著的变化，从而影响稻米对Cd的富集。

图8-11　不同区域pH变化对稻米Cd吸收系数的影响

在规则7中，土壤Cl－也表现出对稻米富集Cd的影响，且这种影响以负作用为主，这与王芳等（2006）的研究结果一致。可能原因是Cl－能与土壤中的金属Cd2＋的络合作用抑制了Cd在作物体内的毒害效应。

除去以上两种经验模型外，还有一种模型称之为概率模型。在一个区域进行土壤-作物系统重金属污染一对一采样分析时，通常会发现农产品重金属含量的增加与土壤重金属含量之间并非线性关系，图8-12就展示了这种情况。就土壤Cd含量与稻米Cd含量之间的线性关系而言，无论是A种情况还是B种情况都不算是线性相关的情况，但这并不表示两者之间真的就没有关系，或者说无法建立两者之间的内在联系。由图8-12可以看出，2个图稻米重金属含量分布情况是不一样的，在土壤Cd含量在1mg/kg范围内时，A图稻米Cd含量一般不超过0.5mg/kg，而B图为不超过1.3mg/kg，在土壤Cd含量在1～3mg/kg时，两者之间也有一定的差别。我们称这种差异性为数据分布的差别，而这种差别可以用蒙特卡罗方法进行模拟，结果如图8-13～图8-16所示。由此可见，在不同重金属含量不同区域内，水稻对重金属Cd的吸收特征是不同的。

图8-12　不同区域土壤重金属Cd含量变化对应稻米重金属Cd含量变化关系图

图8-13　图8-12（A）中土壤Cd含量小于1mg/kg时稻米Cd含量拟合分布图

图8-14　图8-12（A）中土壤Cd含量大于1mg/kg时稻米Cd含量拟合分布图

图8-15　图8-12（B）中土壤Cd含量小于1mg/kg时稻米Cd含量拟合分布图

图8-16　图8-12（B）中土壤Cd含量大于1mg/kg时稻米Cd含量拟合分布图

蒙特卡罗法不同于确定性数值方法，它是用来解决数学和物理问题的非确定性的（概率统计的或随机的）数值方法。它是用一系列随机数来近似解决问题的一种方法，是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处理数学问题的一种手段。运用该近似方法所获得的问题的解更接近于实验结果，而不是经典数值计算结果。它起源于早期的用概率或近似概率的数学思想。

在应用蒙特卡罗方法分析污染物分布特征时，要求产生的随机数序列应符合各风险因子特定的概率分布，而产生各种特定的、不均匀的概率分布的随机数序列，通常采用的方法是先产生一种均匀分布的随机数序列，再设法转换成特定要求的概率分布的随机数序列，以此作为数字模拟试验的输入变量序列进行模拟求解。步骤为：①建立概率模型，即对各评价因子构造一个符合其特点的概率模型（随机事件、随机变量），包括对确定性问题，需把具体问题变为概率问题，建立概率模型。②产生随机数序列，作为健康风险评价的抽样输入进行数字模拟试验，得到大量的模拟试验值。③对结果进行统计处理（计算频率、均恒等特征值），得出结论和相应的准确度估计。

蒙特卡罗方法的应用有两种途径：仿真和取样。仿真是指提供实际随机现象的数学上的模仿的方法。取样是指通过研究少量的随机的子集来演绎大量元素的特性的方法。例如，通过采集少量的样本来演绎总体的分布趋势。因此，采用蒙特卡罗方法模拟稻米Cd分布中的不确定性因素，通过采集有限的样本来预测总体的情况，在理念上是完全可行的。

土壤-作物转运模型构建方法多种多样，在此再介绍两种方法，一种是决策树方法，另一种是随机森林回归方法，当然其他方法如人工神经网络方法等也可以应用，但在此不做介绍。