基于UNMIX模型和莫兰指数的湖南省汝城县土壤重金属源解析

刘昭玥， 费 杨， 师华定,*， 莫 莉， 漆静娴， 王 超

1.中国环境科学研究院土壤与固体废物环境研究所， 北京 100012

2.生态环境部土壤与农业农村生态环境监管技术中心， 北京 100012

土壤环境质量安全对推动农业、社会经济快速发展和保障人体健康至关重要[1-2]. 在我国经济快速发展的同时，环境成本也在不断增加，重金属通过大气沉降、污水灌溉、工业排放等途径进入土壤[3-4]，使得我国土壤污染面积逐年扩大，污染程度不断加剧. 由于重金属通常可在低浓度下表现出高毒性的特点，且停留时间长、生物有效性持久，因此对生态系统和人体健康产生了极大威胁[5-6]. 另外，土壤重金属污染受到成土母质和工农业等人类活动的多重影响，成因复杂，且具有高度的空间异质性[7-9]. 因此，在当前土壤环境问题日趋严峻的背景下，土壤环境的保护和治理迫切需要准确掌握区域土壤重金属污染的来源，这是实现土壤污染源头管控的前提和关键.

在土壤污染源解析方法中，受体模型法由于可避开各污染源的组成信息，具有直接测定受体环境的优势，是当前污染物源解析研究中最主要也最常用的技术手段[10]. 目前应用较多的受体模型主要包括多元统计分析法[11-12]、化学质量平衡法[13-14](chemical mass balance，CMB)、UNMIX模型[15-16]、正定矩阵因子分析[17-18](positive matrix factorization，PMF)等. 其中，UNMIX模型是一种较新型的多元受体模型，其通过对受体中各物种含量的数据空间进行降维来识别污染源，无需事先掌握污染源个数、成分谱和数据的不确定度等[19]，克服了其他受体模型在结果中给出负贡献、人为添加干扰等缺点，使分析结果更准确、实用. UNMIX模型在土壤重金属源解析领域的应用已逐渐受到关注. 刘玲玲等[20]采用UNMIX模型探究北京城区公园5种土壤重金属的来源，发现Cr主要来自土壤母质和特殊工业源，贡献率为21.38%；Cd主要来自工农业混合源，贡献率为35.43%；Cu、Zn、Pb主要来自交通源，贡献率为43.19%. 张军等[21]使用UNMIX受体模型分析了宝鸡市区工业城市土壤重金属污染的来源途径，发现Cd、Cu和As污染主要由工农业活动所致，Zn和Ni污染由交通排放造成，Cr、Pb和Mn污染主要受“混合源”影响.

但是，由于土壤污染问题的复杂性和多变性，单一的源解析方法难以满足目前对污染源贡献定量解析的要求. 受体模型方法主要集中于数据间的数学关系，忽视了数据的空间属性，缺乏对污染来源的空间特征分析[22]，同时，对其解析结果的合理解释也是受体模型应用过程中的重点和难点. 双变量莫兰指数可通过指数计算来表征两种要素的空间相关性. 周俊驰[23]采用局部莫兰指数识别出研究区内铅锌矿及钨矿主要产区是重金属Cd、As、Hg的高值聚集区. 徐源等[24]使用莫兰指数识别出污染源与不同类别企业的空间关系，辅助说明了源解析的结果. 基于以上分析，该研究拟采用UNMIX模型与莫兰指数结合的方式，对湖南省汝城县土壤中Cd、Hg、As、Pb、Cr、Cu、Ni七种重金属的来源进行分析；在得到UNMIX模型解析结果的基础上，使用局部莫兰指数识别各源的高贡献值区域与潜在污染源的空间关系，一方面可从数据的数学关系与空间属性两个角度给出源解析结果，同时也可辅助说明并验证模型的解析结果，以期为土壤污染的源头管控工作提供参考.

1 材料与方法

1.1 研究区概况

汝城县(113°16′E～113°59′E、25°19′N～25°52′N)位于湖南省郴州市东南部，县域内水系分属湘江、珠江、赣江，有“鸡鸣三省，水注三江”之称. 汝城县位于南岭北麓多金属成矿带，矿藏资源分布较广，矿种较多，特别是有色金属矿十分丰富，矿床、矿点、工矿企业众多. 工矿企业生产排出的废水、废气是该地区土壤重金属污染的重要原因. 湖南省《“十一五”湘江流域水污染防治规划(待批稿)》中指出：随废水流入湘江的有害物质中以有色金属的“贡献”最大，主要有Hg、Cd、Pb、Zn、Cu、Cr等重金属. 郴州市作为湘江水系的土壤重金属污染集中江段之一，仅县属以上21个矿山的废水排放量就高达2 375×104m3/a，境内河流中大量重金属元素含量超过GB 8978—1996《污水综合排放标准》的地面水或工业废水排放标准[25-26]. 湘江流域水系沉积物的重金属含量统计结果显示，Hg、Pb含量的最大值均出现在郴州境内的东河[27]. 湖南省长株潭地区的大气沉降监测数据显示，Pb、As、Cd在大气沉降中的浓度较高[28]. 郴州地区4个不同方位的降尘监测点中，Pb和Cd的检出率分别达到100%和90%[29].

1.2 数据来源

土壤点位数据来自研究区随机布点采样，共计233个点位. 在考虑均衡的情况下进行随机布点，每个点位周围随机采集5个深度为0～20 cm的表层土壤样品，按四分法取1.0 kg为土壤样品，装入样品袋中带回，采样过程中用GPS记录点位位置. 所有土样带回实验室，经室内自然风干，剔除植物残体和石块，研钵磨碎后过100目(孔径相当于0.15 mm)筛. 用pH计测定土壤样品的pH. 土壤样品中的各重金属元素在经HCl-HNO3-HClO4-HF混合消解法进行电热板消解后，采用电感耦合等离子体光谱仪测定其含量. 研究区点位及主要潜在污染源分布见图1.

图1 土壤采样点位及主要潜在污染源分布

图 2 土壤重金属含量空间分布

通过查阅《2019年中国环境统计年鉴》《2019年湖南省统计年鉴》以及实地信息采集等方式核实并收集研究区相关污染企业数据，涉重金属企业共73家，包括56家采选业(36家有色金属矿采选业、20家黑色金属矿采选业)，11家冶炼和压延加工业(7家有色金属冶炼和压延加工业、4家黑色金属冶炼和压延加工业)，4家化学原料和化学制品制造业和2家其他行业(1家公共设施管理业、1家石油加工业). 涉重企业排放的废水、废气、废渣中包含Cd、Pb、Hg、As等多种重金属元素.

1.3 研究方法

1.3.1UNMIX模型

UNMIX模型是美国环境保护局开发的一种较新型的多元受体模型，该模型采用奇异值分解(SVD)方法对数据空间进行降维，可判断源数量、源组成和源对各样本的贡献率. UNMIX模型的原理可以用式(1)[30]表示：

(1)

式中：Cab为第a个样品中第b个物种的含量；Ubk为第b个物种在源k(k=1，…，m)中的含量，代表源的组成；Dak为源k在第a个样品中的总量，代表源的贡献率；S为各源组成的标准偏差.

UNMIX系统通常会建议用户选择最小信噪比大于2〔Min(S/N)>2〕且浓度拟合系数在0.8以上的运算结果.

1.3.2莫兰指数

在众多空间相关性的测度指标中，莫兰指数的应用最为广泛[31]，可分为全局莫兰指数和局部莫兰指数. 全局莫兰指数用于表示属性值在空间上是否聚集，局部莫兰指数可进一步指示属性值高值或低值的聚集区域. 具体计算方法[32]：

(2)

(3)

对计算所得到的莫兰指数，利用Z分布进行显著性检验，检验公式：

(4)

(5)

VarI=E(I2)-E2(I)

(6)

式中，E(I)为莫兰指数的期望值，E(I2)为莫兰指数方差的期望值. 当|Z|>1.96时，说明该要素在95%置信区间内呈明显的聚集或离散特征，否则为随机分布.

1.4 数据处理

土壤重金属含量的描述性统计在SPSS 25软件中完成；使用EPA UNMIX 6.0软件进行污染源解析；基于莫兰指数的空间相关性分析在GeoDa软件中进行.

2 结果与讨论

2.1 土壤重金属含量描述性统计结果

研究区土壤中7种重金属含量的描述性统计结果如表1所示. Cd、Hg、As、Pb、Cr、Cu、Ni含量的平均值分别为0.29、0.18、19.91、44.17、66.31、28.67、25.16 mg/kg，除Cr和Ni外，Cd、Hg、As、Pb、Cu含量的平均值均高于当地背景值，分别为背景值的2.95、1.96、1.42、1.49、1.10倍，说明该区域土壤重金属含量受到人为因素的影响. 根据变异系数(CV)划分标准[36]，Cd、Hg、Pb、Cr、Cu、Ni均表现为高变异(50%100%)，表明研究区各点位重金属含量的变化幅度较大，空间差异明显，局部可能存在重金属富集情况. 以GB 15618—2018《土壤环境质量 农用地土壤污染风险管控标准(试行)》中筛选值为标准，未发现Hg含量超标点位，Cd和As超标点位比例在10%以上，其他重金属超标点位较少，研究区土壤污染程度整体较轻.

表1 土壤重金属含量统计分析

2.2 土壤重金属含量空间分布

在ArcGIS 10.2中分别对7种重金属元素含量进行反距离权重插值，得到研究区土壤重金属含量的空间分布. 如图2(a)(b)所示，Cd和Hg含量的空间分布呈现出明显的面状污染特征，以污染源为高值中心，在附近一定距离处形成环状区域. 同时，Cd和Pb的高值区〔见图2(a)(d)〕并未与企业密集区重合，而是出现在距工业企业一定距离处，因此其可能随工业废气经过了远距离传输. 研究区中北部河流密集区Hg含量明显较高. As含量在研究区中部、南部以及东北部较高〔见图2(c)〕，这些地区均有大量矿产采选、加工企业. 由图2(e)(f)(g)可以看出，Cr、Cu和Ni含量的空间分布相似，即低值区范围较大，高值区主要集中在东北和西南方向. 研究区土壤重金属的分布情况在空间上直观反映了土壤重金属污染的可能来源，为模型解析结果的解释提供了支撑.

图 3 不同源对研究区各土壤点位的贡献值

2.3 基于UNMIX模型的土壤重金属源解析

使用UNMIX 6.0软件对研究区7种重金属的来源进行解析，得到4种源对各点位和各重金属元素的贡献(见图3、4). UNMIX模型结果的拟合系数Rsq.=0.82，代表82%的方差可由该模型解释，大于系统要求的最小值〔Min(Rsq.)>0.8〕；信噪比S/N=2.49，大于系统要求的最低值〔Min(S/N)>2〕. 因此，这4个源得出的解析结果是可信的.

由图4可见，源1对As的贡献率最高，达到93.44%. As常以伴随元素的方式存在于多种重金属矿中，其受热易挥发的特点导致As污染易扩散，有色金属冶炼、矿物资源开发和工业废物排放(如冶炼和化石燃料燃烧产物)是工矿企业周围土壤中砷的主要人为来源之一[37-38]. 结合图2(c)的土壤As含量空间分布可以看到，土壤As含量在企业密集区较高，且高值区面积较大. 同时，从图3可以看出，源1对69、86、114、229号点位的贡献值明显较高，而这些点位附近均存在较多的矿产采选、加工企业. 曾庆庆等[39]研究发现，As含量较高的区域均在矿产采选及冶炼企业附近，工业排放对As的贡献率达到67.4%. 故初步推断源1是以矿业为主的工业源.

图4 不同源对研究区土壤重金属的贡献率

源2对Hg的贡献高达61.43%，在131～136号点位出现了明显的高贡献值，该部分点位分布在研究区北部河流密集区域. Hg元素多伴生于有色金属的硫化物矿床中，研究区内工矿企业众多，在矿石开采、冶炼等矿业活动中，各种形态的Hg可随工业“三废”直接或间接地进入土壤. 该研究区内土壤Hg含量的空间分布结果显示，Hg含量在河流附近较高，可能由于水体受到工业排放的污染. Wang等[40]在对安徽省某流域的河流底泥进行重金属污染源解析时发现，工矿企业排放的废水是河流底泥中重金属的主要来源. 被废水污染河流中的重金属通过河流输送、农业灌溉、地下水渗透等途径进入土壤，在周边地区形成污染. 已有研究表明，灌溉水是土壤中Hg的重要来源之一[41]. 综上，初步得到源2为工业排放导致的污水源.

源3对Cd和Pb的贡献率分别高达79.38%和70.27%，其次是对Cu的贡献率，为38.92%. 从源3对各点位的贡献值来看，高贡献值点位较多，主要分布在研究区西南、东南及中部工业密集处等多个区域，故考虑其为影响范围较大的面源污染. 研究表明，Cd广泛存在于各类矿石中，Pb及其化合物是重要的工业原料[42]. 矿产开采、冶炼以及化石燃料燃烧等工业活动排放的“三废”是土壤中Cd、Pb的重要来源[43-44]. 同时，Cd是一种相对容易移动的元素，可通过风化作用进入环境. 目前我国铅锌冶炼工业原料以精矿为主，Cd等污染物主要存在于精矿中，在铅锌冶炼工业的备料、制酸等工序排放的废气中均有Pb、Cd的排出，并易沉降于企业周围土壤. 研究区土壤重金属含量空间分布结果也显示，Cd和Pb含量的高值区主要出现在距企业密集区一定距离处，并受到大气传输的影响. 关于Cu的来源，已有多项研究[45-46]表明其与有色金属的冶炼有关，受到工业源的影响. 综上，将源3识别为以大气沉降源为主、工业源为辅的混合源.

源4对多种重金属都有一定贡献，其中对Ni和Cr的贡献率最高，分别达到82.95%和82.02%. 从图3 来看，源4对各点位的贡献值大多保持在较低水平，且变化不大. 根据描述性统计结果，Ni和Cr的平均值均未超过当地背景值. 张晓文[35]在对湖南省某工业区土壤重金属污染状况进行研究时也发现，Cr和Ni的含量低于背景值，说明这两种元素受自然因素的影响较大；另外从元素的空间分布情况来看，二者的空间分布相似，低值区面积大、范围广. Xu等[47]对嘉兴市土壤重金属进行源解析时也发现，成土母质在较大程度上影响了土壤中Cr和Ni含量. 因此，将源4解释为自然源.

2.4 主导源空间分布

根据UNMIX模型给出的各源对各点位的贡献率，将对各点位贡献最高的污染源作为主导源进行筛选分类，可以得到4种源各自的高贡献区域(见图5)，有助于增强模型解析的可信度并进行下一步分析.

图5 各污染源主导点位的分布

如图5所示，主导源为源1的点位主要集中在企业密集区，结合土壤重金属含量的空间分布也可以看出，源1贡献最大的是As，其高值区与工矿企业的密集区分布一致. 李艳玲等[48]在使用主成分分析方法对工业城市的农田土壤重金属来源进行分析时也发现，工业污染源的得分高值区主要分布在工业发达区. 这进一步说明了源1为工业直接污染源，企业密集区的土壤重金属含量主要受到工矿企业排放的影响. 以源2为主导源的点位在研究区中部、南部出现沿河流分布的趋势，北部河流密集区点位明显增多，与土壤Hg含量的空间分布(见图2)类似，刘瑞平等[49]在对矿区河流中7种重金属含量进行分析时发现，河水中Hg含量超标最严重，最大超标倍数达到 3 099 倍.

综上，Hg污染在河流沿岸地区较为严重，河流附近土壤的重金属污染主要由受到工业排放影响的水源引起. 以源3为主导源的点位少量分布在工矿企业密集区的西北部，大量集中在研究区东南部，与该区域潜在污染源的分布存在较大差异，说明在受到工业企业影响的同时，大气迁移的作用更为显著. 另外，研究区位于季风气候区，东南风和西北风交替盛行的特点使得企业密集区西北和东南方向的土壤更易受到来自大气沉降的污染，与主导源为源3的点位分布特征一致，验证了源3为大气沉降和工业直接排放混合源的解释. 以源4为主导源的点位数量多，分布广泛而均匀，且主要分布在工矿企业较少的研究区外围. 结合图2中Cu、Cr、Ni的空间分布可以看出，源4 贡献最高的这3种重金属元素含量的高值也存在于远离企业密集的外围区域，而目前已有较多研究[50-52]表明，Cu、Cr、Ni主要受成土母质控制，其含量的高低与母质成分有很大关系，与该文中源4为自然源的结论相符. 主导源的分布情况直观反映了各类污染源的主要影响区域，也增强了UNMIX模型的解析结果的可信度.

2.5 因子贡献率与潜在污染源的空间相关性

在ArcGIS 10.2软件中提取研究区潜在污染源(企业、河流)的核密度数据，运用GeoDa软件计算模型各因子贡献率与潜在污染源核密度的单变量莫兰指数，结果如表2所示. 除源2在P<0.003水平下显著外，其余潜在污染源和因子均在P<0.001水平下显著，Z得分均大于1.96，莫兰指数均为正. 这说明各因子贡献率及潜在污染源的分布在研究区域内存在聚集性特征，单变量莫兰指数分析均通过显著性检验，具有统计学意义，为后续双变量莫兰指数分析提供了理论基础.

表2 研究区潜在污染源核密度及因子贡献率的单变量莫兰指数

为更加准确地反映各因子与研究区内潜在污染源的相关关系，对各源的贡献率与潜在污染源的核密度进行双变量局部莫兰指数分析，得到4个因子贡献率与潜在污染源核密度的空间相关性聚类分布结果(见图6). 聚类类型共分为5类，分别是因子贡献高值与潜在污染源核密度高值聚类(H-H)、因子贡献低值与潜在污染源核密度低值聚类(L-L)、因子贡献低值与潜在污染源核密度高值聚类(L-H)、因子贡献高值与潜在污染源核密度低值聚类(H-L)和无聚类.

注： “不显著”表示因子贡献值与潜在污染源核密度无聚类；“H-H”表示因子贡献高值与潜在污染源核密度高值聚类；“L-L”表示因子贡献低值与潜在污染源核密度低值聚类.

由图6可见，源1贡献值与企业核密度的聚类结果在研究区中部企业密集区几乎全部为H-H型，即企业数量最多的中部区域，源1的贡献几乎均为高值. 史广等[53]研究发现，冶炼厂附近区域土壤As含量较高，冶炼厂的生产过程对生产区域及周边局部地区产生了严重的As污染. 该研究中源1对As的贡献率高达93.4%，而土壤As含量的变异系数高达115%，属于极度变异，可能受点源污染的影响较大，符合工矿业直接污染的特征. 在局部莫兰指数分析中，源1的高贡献值与企业的高密度存在对应的空间关系，说明该类污染在空间上未经大气或河流的远距离运输，进一步证实了源1为工矿企业生产活动对周边土壤造成的直接污染.

源2贡献值与河流密度的聚类结果显示，显著聚类类型绝大部分为H-H型和L-L型，即源2的高贡献值几乎仅在河流附近出现，说明源2贡献值的高低与河流位置具有较强的相关性. 该研究区内工矿企业沿河分布的趋势明显，大量工矿企业建在河流附近，采矿废水和选矿废液的直接排放、采矿废石和尾矿等固体废弃物的堆积和淋滤，可导致重金属元素随工业废水进入附近河流[54]. 有研究[28]表明，汝城县所处的郴州市是湘江水系中土壤重金属污染集中江段之一. 当地表水达到污染水平时，受污染河流中重金属元素以灌溉用水、渗滤等形式进入周边土壤，造成土壤重金属污染[55]. 我国有26.4%的土壤污染是由污水灌溉造成的[56]. 陈志凡等[57]研究发现，黄河中下游某城郊农田土壤中As含量的高值区分布在当地河流附近，判断其主要受到污灌的影响. 该研究中，源2对Hg的贡献值最高，而图2显示土壤Hg含量在河流附近出现高值，综合局部莫兰指数分析结果，确定源2为污水源.

在将源3分布与企业核密度和河流密度做局部莫兰指数分析时发现，聚类结果中虽然存在部分H-H型区域，但并未与企业或河流密集区完全重合，还存在相当面积的L-H型聚类；另外，在其与企业、河流密度的聚类中，研究区东南部始终存在较大面积的H-L型聚类，且贡献值与河流密度的H-H型聚类仅在同时有大量企业分布的中部地区出现，北部河流附近并未出现H-H型聚类. 这说明企业分布与源3贡献率存在一定相关性，但并不是造成源3高贡献值的主要原因，而河流与源3的相关性不大，使其高贡献值区域偏离潜在污染源分布的原因更可能是大气的迁移作用. 已有研究[58-60]表明，大气干湿沉降是造成土壤重金属污染的重要原因之一，特别是Cd、Pb、Hg更易通过大气沉降进入土壤. 综合考虑该研究中源3对Cd、Pb、Cu的高贡献率与空间分析结果，证实了源3为大气沉降为主、工业排放为辅的混合源.

在源4与企业、河流密度的聚类结果中，源4高贡献值显著聚集的区域几乎全部为H-L型聚类，即源4的贡献高值集中出现在远离企业、河流的地区；而企业、河流附近的显著聚类类型也几乎全部为L-H型，即企业、河流附近未出现源4的高贡献值，这些区域受源4的影响较小. 这说明源4的贡献值几乎没有受到企业、河流的影响，其主要贡献区域在远离企业、河流的研究区外围，符合自然源的特征. 另外，源4对Cr和Ni的贡献率均在80%以上，在重金属含量的空间分布结果中，Cr、Ni在土壤中的分布规律十分相似，说明其受到同种因素影响的可能性很大. 因此，源4为自然源的解释是合理的.

该研究在使用UNMIX受体模型进行源解析的基础上，结合重金属含量空间插值和双变量莫兰指数分析，在空间关系上验证了模型的解析结果，补充说明了研究区土壤重金属的来源. 按照贡献率的高低，研究区土壤重金属的4种来源为自然源、大气沉降与工业直接排放混合源、污水源和工业直接排放源，其贡献率分别为41.87%、33.10%、13.27%和11.76%，其中人为来源的贡献率为58.13%，是土壤中重金属的主要来源. 值得注意的是，工业集聚区的工矿企业生产活动除直接污染土壤外，通过排放废水、废气污染空气和水环境，进而对土壤产生的间接影响也尤为突出. Liang等[61]在对湖南省涟源市土壤重金属来源进行研究时发现，大气沉降是仅次于成土母质的第二大来源. Wei等[62]发现，湖南省水口山矿业发达地区土壤Cd、Pb等重金属主要来自冶炼厂废气的排放. 该文的研究结果与上述结果基本一致. 因此，在矿业发达的湖南省，土壤污染治理也应同时关注工矿企业的废水、废气等间接影响因素，重视大气沉降、水源污染等对土壤造成的影响，从而有针对性地开展企业监管与土壤污染防控工作.

3 结论

a) 湖南省汝城县土壤中除Cr、Ni外，Cd、Hg、As、Pb、Cu含量的平均值均高于当地背景值. 7种重金属含量的变异性强，变化范围大. 部分点位的Cd、As、Pb、Cu、Cr、Ni含量超过GB 15618—2018《土壤环境质量 农用地土壤污染风险管控标准(试行)》的筛选值，Cd和As污染相对突出. 该区域土壤环境受人为活动影响较大，存在局部的重金属污染，污染程度整体较轻.

b) 研究区土壤中7种重金属受自然源、大气沉降与工业直接排放混合源、污水源和工业直接排放源的影响，其贡献率分别为41.87%、33.10%、13.27%和11.76%. 其中，自然源是Cr和Ni的主要来源，Cd和Pb主要受大气沉降和工业直接排放的影响，水源污染对Hg的贡献率最高，As主要来自直接工业源，Cu则受到自然源和工业源的综合作用.

c) 经过插值和莫兰指数等空间分析方法的验证，UNMIX模型解析效果良好. 受体模型和空间分析方法的联合应用可以更加全面地解析土壤重金属来源.