基于地统计学的土壤污染源解析模型差异对比

2022-11-04王彬姜坤师华定徐嘉礼吴海波罗春晖

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

王彬，姜坤，师华定，徐嘉礼，吴海波，罗春晖

（1.浙江益壤环保科技有限公司，浙江绍兴 312000；2.绍兴文理学院，浙江绍兴 312000；3.生态环境部土壤与农业农村生态环境监管技术中心，北京 100012；4.中国环境科学研究院土壤与固体废物环境研究所，北京 100012）

土壤是地球系统的组成部分，更是人类及动植物赖以生存且不可替代的自然资源[1-2]。随着工业化、农业化的不断推进，我国土壤重金属污染形势逐渐加重，各个区域的污染程度不均，污染成因也不尽相同[3-8]，污染源解析显得尤为重要。近年来，国内外学者在土壤重金属源解析方面进行了大量研究[9-11]，主要采用源排放清单法、化学质量平衡模型、多元统计模型、先进数学算法以及空间分析法。其中，源排放清单法是基于污染源重金属投入通量直接建立污染源清单数据库，例如SHI等[12]汇总了我国农田土壤重金属的排放清单；化学质量平衡模型是基于污染源和土壤重金属的元素组成，根据质量守恒定律建立平衡模型，例如LIU等[13]、CHEN等[14]分别利用元素比值法、同位素比值法对研究区进行了污染源解析；多元统计模型是基于土壤重金属元素组成的源解析方法，例如杜展鹏等[15]采用PCAAPCS-MLR模型得到了污染源对滇池草海和外海水质的贡献程度，柴磊等[16]基于PMF模型对兰州耕地土壤重金属污染源进行了解析，刘玲玲等[17]采用UNMIX模型对北京城区公园土壤重金属污染源进行了解析；先进数学算法是基于土壤重金属元素组成和其他变量的解析方法，如HU等[18]利用条件推断树和有限混合模型识别了珠江三角地区土壤重金属污染来源和作用机制，宋志廷等[19]利用随机森林模型构建了天津武清地区表层土壤重金属源汇的量化关系；空间分析法是基于土壤重金属元素组成和空间关系的解析方法，如任加国等[20]利用空间关联指数法对沱江上源支流土壤重金属的空间变异及其相关性进行了解析。

以上关于源解析的研究方法较为单一，且局限性较大，如源排放清单法需要大量的数据清单，同位素比值法只能针对特定污染物（Cd、Hg、Pb）进行溯源，先进数学算法无法判断解析结果对数据空间变异性的解释程度等，且目前在流域尺度上基于地统计学方法将多种源解析模型进行对比分析的研究极少，而松阴溪流域是典型的灌溉水型河流，其周边土地利用方式包括了耕地、园地以及成林，重工业也较为发达，适用于土壤污染源解析研究，因此本研究以该流域为研究对象，将多元统计模型和空间分析方法相结合，对研究区土壤重金属污染类型及其贡献率进行深入剖析，并将不同模型的解析结果进行对比分析，旨在为区域性土壤污染源解析工作提供方法借鉴和参考。

1 材料与方法

1.1 研究区概况

松阴溪是瓯江上游主要支流，位于浙江省境西南部，源出遂昌县垵口乡北园附近，先自南而北，后折向东南，流经松阳县，在丽水市大港头入瓯江中游大溪。松阴溪全长114 km，河宽约100 m，流域面积2 055 km2，灌溉着浙南的松古盆地。流域内有耕地约10 666.7 hm2，人口约22万人，占松阳县总人口的90%以上。该流域周边道路纵横交错，分布着大量的企业，尤其是制造业和采矿业居多，且附近居民长期使用流域内的水资源进行农田灌溉，对园地进行大量的农药喷洒，其周边土壤环境污染问题日趋严重。

1.2 样品的采集、制备与测定

利用ArcGIS10.6软件对研究区表层土壤点位进行布设，所有点位均布设在农用地（耕地、林地、园地），原则上按照每2 500 m×2 500 m网格布设一个点位，以网格中心点作为采样点，根据周边具体情况适当调整。共布设62个表层土壤点位，对所有点位做2 km缓冲区后划定影响范围，并进行适当调整，如图1所示。根据布设点位的经纬度，利用GPS对每个计划采样点进行精确定位，误差不超过30 m，采样人员到达计划采样点位后，需观察其是否符合土壤采样的代表性要求，在允许范围内优选采样点，位移距离不超过100 m，且低洼地、陡坡地、住宅、道路、沟渠附近等不布设点位。

图1 研究区采样点位、企业分布及土地利用情况Figure 1 Sampling points,enterprise distribution and land use in the study area

样品采集以确定点位为中心划定采样区域，一般为50 m×50 m，采用蛇形法进行混合取样，采样前先用铁铲切割一个大于取土量的25 cm深的土方，然后用木铲去除铁铲接触面后装入样品袋，垂直挖土，保持采样量上下一致，之后将采集完的土壤样品置于阴凉处保存。样品制备时，首先取适量土壤样品平铺在玻璃板上，避免阳光直射且室内温度不超过40℃，自然风干，去除杂质，用木棒将土块粉碎后过60目尼龙筛，再过100目尼龙筛，混合均匀后待测。采用原子荧光法和等离子体质谱法（ICP-MS）测定重金属含量，分析方法精密度和准确度采用国家土壤标准物质GSS-15和室内平行样品进行质量控制。

1.3 研究方法

1.3.1 PCA-APCS-MLR模型

PCA-APCS-MLR（主成分-绝对主成分-多元线性回归）模型主要是通过降维充分解读样本数据所携带信息的差别，最后根据多元线性回归系数来计算各类污染源的贡献率[21-22]。首先需要计算样本浓度的主成分得分，其公式为：

式中：(Az)k为主成分的得分值，i为样本元素所在列号，wij为第j主成分的因子矩阵系数，zk为k观测点污染物浓度标准化后的值，ck为k观测点污染物浓度，c为污染物浓度的平均值，σ为污染物浓度的标准差，由于(Az)k为标准化的值，不能直接用于计算主成分的原始贡献，必须将标准化的因子得分转化为非标准化的绝对因子得分，才能用于主成分对污染指标的贡献分析，其计算公式为：

式中：APCSjk为绝对主成分得分值，(Az)jk为第j主成分的得分值，(A0)j为0值下主成分得分值，i为样本元素所在列号，Sij为第j主成分的因子矩阵系数，(Z0)i为观测点污染物浓度0时标准化后的值，ci为污染物浓度的平均值，σi为污染物浓度的标准差。

以污染物浓度作为因变量，以绝对主成分得分为自变量，建立多元线性回归方程，并利用回归系数计算污染源的贡献率，其计算公式为：

式中：m代表某种污染源类型，aim为回归系数，bi为常数项，PCim为针对元素i、污染源m的贡献率，aim×为元素i的所有样本绝对主成分因子得分平均值。

该模型可通过多元线性回归方程对污染物含量进行反预测，从而得到每个样点对应的预测值。一般采用标准误差和决定系数来分析模型精度，计算公式为：

式中，SE为标准误差，R为决定系数，n为样本数量，xi和yi分别为污染物含量的实测值和预测值。SE越小，R越大，代表模型精度越高，反之越低。

1.3.2 PMF模型

PMF（正定矩阵因子分解）模型法是对因子分析模型的一种优化，该模型已被成功运用于土壤污染源解析方面，它不仅克服了因子分析中负贡献的情况，且模型运行不需要属于源成分谱，通过样本不确定度来预估模型误差，故其解析结果更有意义[23-26]。PMF模型将原始数据分解成因子成分矩阵（F）和因子贡献矩阵（G）以及一个残差矩阵（E），之后根据不确定性计算各组分的误差，并通过最小二乘法计算排放源的贡献度。PMF模型计算公式为：

PMF模型为得到最优的因子成分矩阵（F）和因子贡献矩阵（G），从而定义了一个目标函数Q，使其趋于自由度最小。其目标函数计算公式为：

式中：Xij为第i个样品中第j个元素的浓度，p为污染物来源数量，gik为源k对第i个样品的贡献，fkj为源k中第j个元素的浓度，eij为残差矩阵即PMF模型中未能解释样品浓度矩阵Xij的部分，uij为第i个样品中第j个元素浓度的不确定性大小。不确定度可以根据元素的实测浓度与相应的方法检出限（MDL）的关系确定，如果元素浓度的实测值没有达到或者刚好达到方法检出限时，不确定度计算公式为：

当元素远远超出了相应的检出限时，不确定度的计算公式为：

式中：σ为标准差，c为元素浓度。

该模型运行结果对各污染物的标准误差和决定系数均进行了统计，并绘制了污染物含量预测值与实测值的受体拟合图，用以分析该模型的精度。

1.3.3 UNMIX模型

UNMIX模型是一种适用于解决混合问题的数学模型，受体点被认为是未知源组分的线性组合，这些来源对每个样本都贡献了一个未知的数量[27-29]。该模型假定来源的构成和贡献都是正的。利用特定物种选择的浓度数据估计每个样本的源个数、源组成以及源对该样本的贡献，即在m个源的N个土壤样品中，n个被分析的物种中的某一物种j的浓度可以用公式表示：

式中：Cij代表第i个样品中第j个物种的浓度；Ujk代表第j个物种在源k中的质量分数，表示源组成；Dik代表源k在第i个样品中的总量，视为源贡献率；S代表各个源组成的标准偏差。

该模型运行结果同PMF模型类似，可直接观察到污染物含量预测值与实际值的拟合情况，从而进一步判断其精度。

2 结果与讨论

2.1 土壤重金属含量统计

研究区土壤重金属含量统计情况见表1。对62个土壤样品的统计结果进行分析可知，土壤pH的平均值和中位值分别为4.93和4.95，说明研究区土壤整体上呈弱酸性。Cd、Hg、As、Pb、Cr、Cu、Zn、Ni 8种重金属含量的平均值和中位值均未超过《土壤环境质量农用地土壤污染风险管控标准（试行）》（GB 15618—2018）风险筛选值，说明研究区土壤重金属含量大部分处于较低水平。仅Cd、As两种重金属含量的最大值超过了风险筛选值但低于风险管制值，说明这两种重金属对研究区造成了一定污染，其他重金属污染较轻。从偏度来看，数据分布均为右偏度，但经过Log变换后均符合正态分布；从峰度来看，数据均比正态分布高耸狭窄，说明其数据更集中在平均值附近。

表1 松阴溪流域土壤重金属含量描述性统计Table 1 Descriptive statistics of soil heavy metals content in Songyin creek

2.2 不同源解析模型结果

2.2.1 PCA-APCS-MLR模型结果

利用SPSS26.0软件对重金属进行因子分析前，需要先对样本数据进行检验，其中KMO检验结果为0.68，满足最低0.5的要求；Bartlett检验观测值为311.42，显著性为0.00，拒绝零假设，说明该样本数据总体适合因子分析。对研究区样本数据进行主成分提取，且选择特征值大于1的因子作为主因子，最终得到3个主因子。为了更加突出主因子的载荷能力，使结果更易于解读，采用最大方差法对主因子进行旋转，其结果见表2。从表2可以看出，3个主因子特征值变化幅度较小，累计方差贡献率为79.94%，反映了样本携带的绝大部分的变异，未旋转前3个主因子可解释方差分别为40.83%、24.93%、14.18%，旋转后变为36.05%、29.42%、14.47%，模型得到优化。从表2可以看出，旋转后的成分矩阵可以更好地识别因子特征，第一主因子的特征元素为Cd、Pb、As、Cu、Zn，第二主因子的特征元素为Cr、Ni，第三主因子的特征元素为Hg。通过PCA-APCS-MLR模型进一步计算各污染因子对各个元素所有样点的贡献率，负贡献率改为0，其计算结果见图2。可以看出，大部分贡献率较高的重金属为主因子对应的特征元素，而第二主成分对Cu的贡献率也较高。

图2 PCA-APCS-MLR模型对各重金属污染源的贡献率Figure 2 Contribution rate of heavy metal pollution sources by PCA-APCS-MLR

表2 PCA-APCS-MLR模型主因子解及成分矩阵Table 2 Main factors solution and component matrix of PCA-APCS-MLR

2.2.2 PMF模型结果

利用EPA PMF5.0软件将样本数据及其不确定度（根据方法检出限求得）输入模型，观察8种重金属数据的信噪比，将比值＜0.5设置为“bad”，0.5～1.0设置为“weak”，比值＞1.0设置为“strong”。利用旋转工具观察不同旋转因子Fpeak下目标函数Q值的大小，多次迭代运行，发现旋转因子Fpeak=-0.5时Q值最小，最终模型结果见图3。从样点来看，3种污染源的贡献变化幅度均较大，说明污染源的区域性差异较为显著；从元素来看，污染源1对Hg的贡献较大，污染源2对Cd、As的贡献较大，对Pb、Cu、Zn也有一定的贡献，污染源3对Cr、Ni的贡献较大。

图3 PMF模型对各重金属污染源的贡献率Figure 3 Contribution rate of heavy metal pollution sources by PMF

2.2.3 UNMIX模型结果

利用EPA UNMIX6.0软件对样本数据进行模拟，首先对数据进行处理，利用Suggest Exclusion工具自动将浓度方差大于50%的数据进行筛选，结果显示8种重金属中仅Hg的浓度方差略高于50%，由于该元素为源解析重要因素，因此也可纳入模型。UNMIX模型是利用几何的“边”对样本数据进行分析，以各重金属的浓度和作为基底，单个重金属元素作为因变量，从而分析各个组分的边缘诊断结果，通过运行该模型可以发现，除Cd、Hg外都具有明显的上边界或下边界，通常认为大部分元素具有边界时，模型基本是可靠的，且Cd、Hg两种元素均无缺失值，各组分信噪比均符合要求，因此保留所有组分。

将8种重金属元素的浓度均作为初始物种，将浓度和（total）作为模型输入组分的总物种和标准物种，计算出符合要求的3源方案，结果显示MinRsq=0.86，MinSig/Noise=2.83，均符合模型要求，说明该模型可以解释86%的物种方差，源解析方案结果可靠。利用Excel进一步计算源贡献率，如图4所示。可以看出，污染源1对Cr、Ni的贡献较大，其他元素相对较小；污染源2对Hg的贡献较大，对Pb、Zn也有一定的贡献；污染源3对Cd、As的贡献较大，对Pb、Cu、Zn也有一定的贡献。

图4 UNMIX模型对各重金属污染源的贡献率Figure 4 Contribution rate of heavy metal pollution sources by UNMIX

2.3 污染源空间识别

为进一步对污染源进行识别，利用ArcGIS10.6对前文3种模型的样点贡献值进行地统计插值分析，如图5所示。对于PCA-APCS-MLR模型来说，污染源1贡献值较高的大部分区域分布的企业较少，道路情况一般，周边无地表水，且无其他明显特征，因此判别为“自然源”；污染源2贡献值较高的区域主要分布着大量企业，主要包括采矿业和制造业，相关研究表明[30-32]，企业在冶炼加工、机械制造、废气排放时均会造成重金属污染，因此判别为“工业源”；污染源3贡献值较高的区域主要分布着复杂的道路（省道、县道以及城镇街道），且建有大型的停车场，由上文可知，该污染源的特征元素是Hg，根据以往研究表明[33-35]，交通因素对Hg的累积具有较大的作用，因此判别为“交通源”。对于PMF来说，污染源1贡献值较高的区域与PCA-APCS-MLR模型的污染源3类似，且部分区域属于农田和果园，当地村民利用附近未处理过的地表水进行灌溉，且相关研究表明[36-37]，果园采用大量的化学药剂也会导致土壤重金属的累积，因此判别为“交通-农药-污灌源”；污染源2贡献值较高的区域无明显特征，因此判别为“自然源”；污染源3与PCAAPCS-MLR模型的污染源2较为类似，因此判别为“工业源”。对于UNMIX模型来说，污染源1贡献值较高的区域分布着大量的企业，因此判定为“工业源”；污染源2贡献值较高的区域与PMF模型的污染源1类似，因此判定为“交通-农药-污灌源”；污染源3贡献值较高的区域无明显特征，因此判定为“自然源”。

图5 各源解析模型中污染源贡献率的空间分布Figure 5 Spatial distribution of pollution source contribution rate in each source analysis model

2.4 源解析模型结果对比

为了更好地对比源解析模型结果，首先需要对其模型精度进行对比。其中PCA-APCS-MLR模型通过多元线性回归计算各重金属元素的标准误差和决定系数，而PMF模型和UNMIX模型在运行结果中可直接观察预测值和实测值的拟合曲线，并给出了标准误差和决定系数（相关系数的平方），统计结果见表3。从元素上来讲，除Pb、Zn外，其他元素的标准误差均较低，且大部分元素的决定系数在0.7以上，说明各模型整体拟合效果较好；从模型上来讲，各模型的标准误差和决定系数较为相似，其中PMF模型的标准误差相对较小，可能是由于该模型在分析时添加了样本的不确定性数据，从而降低了整体的预测误差。

表3 各源解析模型对重金属元素的拟合精度Table 3 Fitting accuracy of each source analytical model to heavy metal elements

结合地统计学方法对各模型的污染源类型及其贡献进行分析，可以发现3种模型之间存在一定的差异性和相似性，从污染源类型来讲，通过以上3种源解析模型都能识别出工业源与自然源，但利用PCAAPCS-MLR模型仅识别出交通源，而利用PMF模型、UNMIX模型却识别出交通-农药-污灌源，这可能是由于PCA-APCS-MLR模型缺乏敏感性，对一些贡献度较小的污染源不能甄别处理，只对主要的源成分进行筛选和识别，而PMF模型、UNMIX模型对源成分比较敏感，且对于相似的源成分进行归类，因此两种模型没有将交通源、农药源以及污灌源进行区分。从污染源贡献来讲，PCA-APCS-MLR模型的源总体贡献率差异较大，而PMF模型、UNMIX模型差异较小，可能是由于通过PCA-APCS-MLR模型分析未识别出农药源和污灌源，且容易出现负贡献，使其他源对As、Zn、Cr、Ni、Hg的贡献度特别高，而PMF模型、UNMIX模型事先假设了受体点没有负贡献，且对于源识别较为全面，因此其源总体贡献率的计算较为准确合理。总体来说，PCA-APCS-MLR模型在判别主要污染源类型时具有较大优势，但容易出现负贡献，贡献率容易出现偏高或偏低的情况；而PMF模型、UNMIX模型难以区分成分相似的源类型，但在计算源贡献率时具有一定的优势，不会出现负值或偏高或偏低的情况。因此，在土壤污染源解析中，有时需要将不同的受体模型结合起来，更加精确地对源类型及其贡献率进行综合判别和计算，对研究区土壤污染源类型和总贡献率进行重新整合和计算，其中源类型遵循细分原则，源贡献以PMF模型、UNMIX模型为主，得到最终结果为自然源30.1%、工业源37.4%、交通源18.0%、农药-污灌源14.5%。