区域农田土壤重金属污染源解析方法研究现状及展望

Research Status and Prospects of Apportionment Methods for Heavy Metal Pollution Sources of Farmland Soil at Regional Scale

摘要: 重金属污染来源精细解析是区域农田土壤开展污染消减与阻控措施的首要措施。农田土壤重金属来源复杂多样，且污染物的迁移与累积过程受到地形、气象、植被等多种因素的影响，导致污染源对区域不同位置农田产生差别化污染贡献。为精细解析该差别化污染贡献，需使用多个类别型和连续型代理变量表征各污染源空间输出强度及影响因素，并建模代理变量与农田土壤重金属含量之间的关系。而在区域尺度这一关系往往存在空间异质性，现有方法能否有效解决该问题。对此，本文在分析农田土壤重金属污染来源及其空间异质性污染贡献的基础上，梳理了区域农田土壤重金属污染源解析方法的研究现状及其发展趋势，并展望了新的源解析方法，为我国各级政府和相关组织制定既经济又环保的农田土壤重金属污染风险防控措施提供技术支撑。

Abstract: Fine apportionment of the heavy metal pollution sources is the prerequisite for pollution reduction and prevention for farmland soil. The sources of heavy metals in farmland soil are complex and diverse, and the migration and accumulation process of pollutants are affected by multiple factors such as topography, meteorology, vegetation, etc., resulting in different pollution contributions from pollution sources to farmland soil in different places in the region. In order to analyze the different pollution contribution accurately, it is necessary to use multiple categorical and continuous proxy variables to characterize the emission intensity of each pollution source and corresponding influencing factors, and to model the relationship between proxy variables and heavy metal content in farmland soil. The relationships between the heavy metal contents and proxy variables are spatially heterogeneous, thus whether existing methods can solve this problem deserves investigation. To answer this question, we first analyze the sources of heavy metal pollution in farmland soil and their spatially heterogeneous pollution contribution, and then summarize the research status and prospects of apportionment methods of heavy metal pollution sources of farmland soil at a regional scale. This study aims to provide technical support for the government and related organizations to take measures to prevent and control heavy metal pollution risks in farmland soil that are both economical and environmentally friendly.

1. 引言

农田土壤是农业生产的重要基础。然而，我国农田土壤重金属污染形势严峻 [1] ，由于土壤重金属具有隐蔽性、滞后性、累积性以及不可逆性等污染特点，农田土壤一旦被重金属污染便难以修复，会长期影响农产品数量和质量安全等，最终危害人类健康和社会经济发展。对此，国务院于2016年印发《土壤污染防治行动计划》，制定了全国土壤污染防治工作的行动纲领，重点部署了截止2050年的土壤污染防治工作任务和目标，要求严控新增污染。而农田土壤重金属风险防控的关键在于源头阻控 [2] 。因此，为了制定既经济又环保的污染消减与阻控措施，首先需要精细解析污染来源及其对区域不同位置农田的差别化污染贡献。

农田土壤重金属源于多种自然源与人为源的叠加，同时各类污染物的迁移与累积受到地形、气象、植被等多种因素的影响。为了精细解析污染源对区域不同位置农田的差别化污染贡献，需要使用多个类别型和连续型代理变量(本文称为“多维混合类型代理变量”)表征各污染源空间输出强度及影响因素，并建模代理变量与农田土壤重金属含量之间的关系。由于污染源的空间位置和尺度的差异，以及污染物迁移与累积特征的空间相关性和异质性，导致这种关系在局部地理空间范围内存在相关性，但是在区域不同部分间经常存在显著差异 [3] [4] ，即存在空间异质性关系。定量解析这种空间异质性关系不仅能够获得污染源对区域不同位置农田的差别化污染贡献，还能提高对区域总体污染贡献解析的准确性，而忽略这种关系会导致全局建模的失败 [5] 。

通常，土壤重金属污染源解析方法按照分析对象可分为污染源分析、受体分析、污染源与受体联合分析三类，能够对污染来源进行定性或定量研究。但上述被广泛应用的方法无法有效处理空间异质性关系。对此，污染源与土壤重金属之间的空间位置关系已被越来越多的学者认识到其对污染源解析的重要意义。

空间异质性关系建模是空间统计领域的基本问题。目前针对该问题的研究已经取得了部分进展，如地理加权回归与空间多水平建模。但是这些方法无法有效处理多维混合类型代理变量。而机器学习领域的基于结构化算子的基函数扩展模型(如提升树、SVM、神经网络等)在处理多维混合类型代理变量方面具有显著的优势，但不能建模空间异质性关系。因此，需要结合两种方法的优势，研发一种能够处理多维混合类型代理变量的空间异质性关系模型，用于精细解析具体污染源对区域不同位置农田的差别化污染贡献。

本文首先分析了农田土壤重金属污染来源及其空间异质性污染贡献，在此基础上梳理了区域农田土壤重金属污染源解析方法的研究现状，并探讨了其发展趋势，最后展望了新的源解析方法。

2. 农田土壤重金属污染来源及其空间异质性污染贡献

土壤中重金属元素的来源分为自然源和人为源两大类 [6] [7] [8] [9] [10] 。自然源包括成土母质、岩石风化、火山喷发等，它们在气候、生物、地形等因素作用下，形成了土壤重金属背景值 [11] [12] 。人为源主要包括工业生产、农业生产、交通运输和居民生活等，通过大气沉降、灌溉与径流、固废与堆肥向农田输入重金属，其输入和累积过程受到气象、地形和植被等因素的影响 [13] [14] [15] 。污染源的尺度差异及影响因素的空间相关性和分异性，导致污染源对土壤的污染贡献同时呈现空间相关性和异质性。对于自然源贡献，除了全国尺度上的地带性分布与东西共轭关系外，在市县以及更小的尺度上，由于局部地形和生物等因素的影响，形成了垂向分异性规律和斑块式的空间异质性分布格局 [16] 。人为源中，点源污染的污染物以企业或生产场地为中心沿风向和水向扩散，污染贡献在空间上呈现非线性渐变趋势 [4] ，而非点源污染则存在明显的斑块特征(如污灌、农药化肥等)或条带化分布(如交通) [17] [18] 。因此，农田土壤重金属污染源多样且影响因素复杂，且空间尺度差异较大，污染源对区域农田土壤重金属的污染贡献既存在局部的空间相关性又存在全局上的空间异质性 [19] [20] 。

传统基于产排污系数的污染源强计算方法只能得到排放总量，无法反映污染源对区域不同位置的污染输出量。因此一般使用代理变量反映其输出强度的空间分布，如采用土壤母质、土壤类型代理自然源，使用与污染源间距离的倒数反映点源污染强度，使用污灌分布、土地利用类型、种植模式、距离道路距离或道路密度等反映非点源污染强度 [17] [21] [22] 。在受体分析方法中，还会使用主成分得分表征对应类型污染源的强度 [23] 。虽然代理变量能够反映污染源输出的空间差异，但它并不是直接的污染输出量，还需要结合污染迁移与累积影响因素，并与土壤重金属污染含量进行建模校正，才能定量计算各自的污染贡献 [24] 。而在区域尺度上，代理变量与真实污染输出分布之间难以存在全局一致的线性比例关系，而且未知历史污染源的贡献以及未能被考虑的影响因素均存在空间异质性，导致各类代理变量与土壤重金属含量之间的关系在区域不同部分存在差异，即存在空间异质性关系(统计学上称为空间非平稳关系或随空间变化的关系) [5] [24] 。

因此，各类污染源的污染贡献存在局部的空间相关性与全局的空间异质性，且空间变异尺度差别较大。可以使用多维混合类型代理变量表征污染源在区域不同空间位置的污染输出量及其影响因素变化，但是代理变量与土壤重金属含量之间存在空间异质性关系。在区域农田土壤重金属污染源解析模型构建中，需要妥善处理这种空间异质性关系 [25] 。

3. 农田土壤重金属污染源解析方法及其发展趋势

污染源解析是对污染来源进行定性或定量研究的方法。土壤重金属源解析方法按照分析对象可分为污染源分析、受体分析、污染源与受体联合分析三类(见表1)。

Table 1. Apportionment methods for heavy metal pollution sources of farmland soil

表1. 农田土壤重金属污染源解析方法

污染源分析方法主要是源排放清单法。该方法在对各类污染源状况调查和统计的基础上，依据污染源的活动水平与排放因子模型估算污染源排放量 [26] 。详细的污染源清单建立难度较大，且污染源对区域土壤的污染贡献存在空间异质性，不同来源的重金属在土壤中的累积能力也有差异，因此该方法不能估算污染源对区域不同位置农田的差别化污染贡献，难以用于指导精准防控。

受体分析法以污染受体(土壤)作为主要研究对象，目前以传统多元统计方法为主。主要包括相关性分析(Correlation Analysis, CoA)、主成分分析(Principal Component Analysis, PCA)、因子分析(Factor Analysis, FA)、聚类分析(Cluster Analysis, CA)等降维方法，PMF和UNMIX等分解方法，以及主成分分析–多元线性回归(Principal Component Analysis-Multiple Linear Regression, PCA-MLR)、绝对主成分分析–多元线性回归(Absolute Principal Component Analysis-Multiple Linear Regression, APCA-MLR)等组合方法 [27] [28] [29] 。该类方法基于土壤中各类重金属之间的线性关系对土壤重金属进行归类，然后依据先验知识判断每一类土壤重金属可能的污染源，不能精确匹配具体的污染源。分解方法和组合方法虽然能给出每类污染源在样点处的贡献，但是在解析过程中采用全局线性模型，没有考虑污染源与土壤重金属之间的空间异质性关系。

污染源与受体联合分析法通过建立具体污染源与土壤重金属含量之间的关系模型，定量解析各污染源对土壤重金属污染的贡献，主要包括化学质量平衡法、空间模式匹配法、空间回归法和机器学习法。化学质量平衡法通过对应污染源排放成分谱与受体中各污染物的浓度，建立线性比例关系解析污染源贡献，如元素比值法、同位素比值法和形态分析法 [26] 。该方法能够定量解析具体污染源对样点处的局部贡献，但是假设各污染源排放物的化学组成稳定且存在明显差异，不适用于复杂污染源的情况，同时还需要较为完整的污染源指纹图谱，代价较高。空间模式匹配法则只需要污染源位置、类型等简单信息，通过匹配由克里金(Kriging)、序贯高斯模拟(Sequential Gaussian Simulation, SGCS)等空间插值方法获得的土壤重金属空间分布模式，或匹配由多变量因子克里金(Multivariate Factorial Kriging, MFK)得出的各尺度因子空间分布特征，定性识别各尺度污染源 [30] [31] 。

近些年，越来越多的学者认识到污染源与土壤重金属之间的空间位置关系对污染源解析的重要意义，除了将空间模式匹配方法与传统多元统计方法相结合，定性或定量解析污染源 [32] - [37] ，还采用距离线性回归模型(Multiple Linear Regression with Distance Model, MLRD)和地理加权回归模型(Geographically Weighted Regression, GWR)定量解析污染源与土壤重金属含量之间的空间异质性关系 [11] [25] [38] 。但是MLRD只能以距离倒数的方式建模点源污染贡献，忽略了其他类型源的污染贡献，容易过高估计点源的污染贡献。GWR能够定量建模污染源与土壤重金属含量之间的空间异质性关系，但是难以处理表征众多污染源及影响因素的多维混合类型代理变量。机器学习方法在多维混合类型变量处理方面具有明显的优势，逐渐被引入土壤重金属污染源解析领域，如条件推断树(conditional inference tree, CIT) [39] 、Cubist模型 [40] 、随机梯度提升树(Stochastic Gradient Boosting Tree, SGBT) [24] 、随机森林(Random Forest, RF) [41] 。但是这些机器学习方法不能处理空间异质性关系，且不能解析空间差别化污染贡献。最常用的RF在处理多维混合类型变量时得出的全局重要性严重有偏，会误导对污染贡献的判断 [42] [43] 。

因此，传统的污染源分析方法、受体分析方法以及化学质量平衡法，不能匹配具体的污染源或者解析其对区域不同位置农田的差别化污染贡献。利用农田污染大数据以及污染源和受体之间的空间位置关系进行源解析，已成为新的发展趋势。但目前采用的空间模式匹配法、空间回归法和机器学习法，由于不能建模多维混合类型代理变量与土壤重金属含量间的空间异质性关系，难以提高解析精度。

4. 新方法展望

地理学第一定律和第二定律不仅阐述了单个目标变量在地理空间上的相关性与异质性，还包括变量之间关系在地理空间上的相关性与异质性 [44] [45] 。除依照先验知识分区外，监督学习模型是目前最主要的空间异质性关系建模工具。它将土壤重金属污染源及影响因素代理变量作为预测变量，将土壤重金属含量作为响应变量，通过对样点数据的学习构建预测变量与响应变量之间的关系模型，解析土壤重金属污染源贡献 [24] [46] 。目前，主要的空间异质性关系监督学习模型都是从非空间的监督学习模型经过空间化改进而来，按照模型的构建方式可以分为空间全局模型、空间多水平模型和空间局部模型三类(如表2)。

Table 2. Supervised learning models for spatial heterogeneous relationships

表2. 空间异质性关系监督学习模型

空间全局模型是从非空间全局线性模型及其扩展模型中通过两种空间化方式改进而来。第一种方式直接将地理空间坐标作为独立的预测变量加入非空间模型，而不改变方法本身。常见模型包括多元线性回归、广义相加模型以及基于结构化算子的基函数扩展模型。其中基于结构化算子的基函数扩展模型，如梯度提升树(Gradient Boosting Decision Tree, GDBT)、支持向量机(Support Vector Machine, SVM)，采用结构化基函数(或希尔伯特空间的核函数)并使用正则化处理方法，能够快速有效地处理多维混合类型变量并防止过拟合 [46] 。但是这种方法只能够获得响应变量与预测变量之间的全局关系及响应变量的空间依赖，不能给出除坐标外其他预测变量与响应变量之间的空间局部关系。其中GDBT在提升学习过程中，虽然会重点拟合局部样点，但其损失函数依然是全局的，是在全局损失函数控制下的局部重点拟合，因此同样不能得出预测变量与响应变量之间的空间局部关系 [47] 。第二种方式将基函数扩展模型中的参数扩展为空间坐标的函数，得到空间扩展模型(Spatial Expansion Method, SEM) [48] 。该类模型能够解析各预测变量与响应变量之间的空间异质性关系，但是所解析的非平稳关系只能是某种预先设定的变化模式，而且这种变化模式必须是全局稳定的，因此在实际中难以应用。

空间多水平模型通过将非空间多水平模型的属性水平转化为地理空间水平，将水平内的属性分组转化为空间分区，并基于不同空间水平与分区的方差分解将空间异质性关系分解为多层线性关系的组合 [49] [50] ，主要包括随机效应模型、随机系数模型。它们能够考虑预测变量的空间尺度差异，建模多水平异质性关系。但是该类方法需要预先确定关系的空间水平和同水平内关系差异的空间分区，而关系差异的空间分区本身就是一个空间异质性关系探索的过程 [51] 。对于政治经济现象，由于其存在行政区内的一致性和行政区之间的差异性，可以按照行政区划进行分区。但是对于土壤重金属污染这种自然与人为混合作用的现象，难以在建模前完成空间关系差异的分区。同时，这类方法难以处理土壤重金属污染源解析所涉及的多维混合类型代理变量 [52] 。

空间局部模型通过将非空间局部线性模型中的多维属性空间邻近性转化为二维(或三维)地理空间邻近性，避免了由于变量数目增加而带来的“维度灾难”，并利用空间局部相关性关系解析全局空间异质性关系，如贝叶斯空间变系数模型(Bayesian Spatially Varying Coefficient, BSVC) [53] [54] 和GWR模型 [48] 。BSVC基于地理空间邻近关系所产生的自相关性构建空间局部条件概率分布函数，获得地理空间局部位置上预测变量系数的概率分布。但是该方法需要预先设定各参数的先验分布，且难以处理多维混合类型代理变量。GWR基于地理空间距离核函数与加权最小二乘法，通过为每个空间位置建立局部线性模型，进行地理空间异质性关系的解析 [55] 。但是由于它在局部模型构建时采用传统的加权最小二乘法，容易产生共线性问题；同时在处理类别型变量时，会引入过多的哑变量，加剧共线性问题，严重影响模型的准确性与稳定性 [56] 。为了缓解GWR的共线性问题，Wheeler [57] [58] 结合多元线性回归的正则化技术，分别提出了地理加权岭回归(Geographically Weighted ridge Regression, GWRR)与地理加权Lasso回归(Geographically Weighted Lasso, GWL)。它们能够在一定程度上解决GWR的共线性问题，但是由于采用线性基函数，会增加拟合参数进而增大模型方差，而且难以处理类别型变量，依然不能较好地处理农田土壤重金属污染源解析中的多维混合类型代理变量。另外，空间局部模型通过在区域内每个点上单独构建局部相关性模型的方式分解区域的空间异质性关系，不仅带来大量的冗余计算，而且无法处理不同污染源及影响因素的空间变异尺度差别较大的问题。

总之，现有的空间异质性关系监督学习模型，虽然在处理多维混合类型变量和空间异质性关系方面各有优势，但是均难以同时解决。可结合地理空间局部加权与多水平分解，将全局异质性关系转化为能适配污染源及影响因素空间变异尺度的局部相关性关系，同时结合基于结构化算子的基函数扩展模型，增强多维混合变量处理能力，研发能够精细解析污染源对区域不同位置农田差别化污染贡献的模型工具。

5. 结论

农田土壤重金属污染源多样且影响因素复杂，污染源对区域不同位置农田土壤重金属存在空间异质性污染贡献。作为区域农田土壤开展污染消减与阻控的首要措施，精细解析重金属污染源需使用多维混合类型代理变量，并建模代理变量与农田土壤重金属含量之间的空间异质性关系。当前常用的区域农田土壤重金属污染源解析方法无法有效处理这一空间异质性关系，因此，利用农田污染大数据以及污染源和受体之间的空间位置关系进行源解析，已成为新的发展趋势。为实现污染源对区域不同位置农田差别化污染贡献的精细解析，可结合地理空间局部加权与多水平分解，将全局异质性关系转化为能适配污染源及影响因素空间变异尺度的局部相关性关系，同时结合基于结构化算子的基函数扩展模型，研发新的源解析方法工具，为我国各级政府和相关组织制定既经济又环保的农田土壤重金属污染风险防控措施提供技术支撑。

基金项目

国家自然科学基金项目(42271428)。

参考文献

NOTES

^*第一作者：丁海元，对外经济贸易大学统计学院在职人员高级课程研修班学员。

^#通讯作者：高秉博，中国农业大学土地与科学技术学院副教授。