湖南省郴州市苏仙区重点污染企业影响区的土壤重金属污染源解析

徐 源， 师华定,*， 王 超， 费 杨， 于靖靖， 舒 密

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

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

3.北京信息科技大学， 北京 100101

近年来，随着城市化进程的加快，重金属通过污水灌溉、大气干湿沉降、工业废渣等途径进入土壤[1-2]，重金属污染已成为当今土壤污染中污染面积最广、危害最大的环境问题之一[3-4]. Cd、Hg、As、Pb、Cr作为重金属“五毒”元素，在土壤中具有移动性大、毒性高、无法降解等特点，在生产活动中容易被作物吸收富集，不仅严重影响作物的产量和品质，还可以通过食物链在人体积累，危害人体健康[5-7]. 研究土壤“五毒”元素污染的来源以及污染源的主要影响区域，对维持区域环境健康具有重要意义[8-9].

目前，土壤重金属污染来源解析主要利用基于重金属总量的受体模型，常用的受体模型有化学质量平衡(CMB)模型、主成分分析/因子分析-多元线性回归(PCA/FA-MLR)模型、正定矩阵因子分解(PMF)模型、UNMIX模型等[10-11]. 其中，PMF模型具有适用性广、不需要测量源成分谱、能对因子贡献作非负约束等优点，得到广泛的应用. 例如，黄华斌等[12]利用PMF模型解析出研究区土壤重金属的来源；董騄睿等[13]通过PMF模型解析结果与Pb稳定同位素比值结果进行对比分析，指出PMF模型可以较好地应用于土壤重金属源解析研究.

然而，单一的PMF模型的解析结果往往比较笼统，对污染源的判别中主观性较强，且缺乏直观视觉效果，因此较多学者将其与地统计学方法结合运用，进一步得到污染源在空间上的分布状况[14-16]. 这类研究虽对解析结果有了很大改进，但在源成分谱未知的情况下，对于污染物的具体来源依然无法准确识别. 另外，关于PMF模型对土壤重金属的源解析研究中，目前大多仅停留在污染源大类上，如自然源、工业源、农业源等，其中对工业源的探讨并未涉及具体的重点污染行业. 双变量莫兰指数方法(bivariate Moran′s I)通过指数的计算对空间内两种要素作比较，表征两种要素间的空间相关性. 周俊驰[17]采用莫兰指数对耕地土壤重金属进行分析，识别出研究区内Cd、As、Hg这3种重金属在耕地土壤中的富集与工矿开采排放关系密切；于靖靖等[18]利用莫兰指数和广义加性模型探讨了重点污染企业对重金属Cd的影响，表明依据空间自相关理论可以对特定的污染源做很好的识别. 鉴于此，该研究将地统计学方法与PMF模型结合起来，在利用传统PMF模型进行污染源解析的基础上，利用莫兰指数表征不同行业与污染源的空间关系，明确工业源下的不同重点污染企业对土壤重金属的影响，以辅助解析验证PMF模型的有效性，以期为区域企业管理和土壤污染治理提供参考.

1 材料与方法

1.1 研究区概况

研究区域位于湖南省郴州市苏仙区(112°53′55″E～113°16′22″E、25°30′21″N～26°03′29″N)，海拔150～160 m，以山地为主，岗地、水面较少，年降雨量 1 452.1 mm，属中亚热带季风湿润气候，土地总面积为 1 342.27 km2. 该区域地处郴州市西河流域，西河是郴江的一级支流、湘江的三级支流. 苏仙区矿产资源丰富，素有“有色金属之乡”的美誉，有世界罕见的柿竹园多金属矿，国有硚口铅锌矿，玛瑙山锰矿，东波铅锌矿以及年产煤数万吨的许家洞、街洞和栖凤渡煤矿等. 发达的矿业生产给当地带来了一系列环境问题.

1.2 数据来源与预处理

利用ArcGIS 10.3在研究区范围内按照2 km×2 km布设302个网格，用GPS对302个采样点精确定位. 采样时采用梅花形布点法对土壤样品进行采集，用木铲采集每个网格表层(0～20 cm)的土壤样品至少6个，样品混合均匀后按四分法获取每个采样点采集的土壤不少于2 kg，装入样品袋中带回，并记录现场土样的相关信息及周边土地利用现状. 将采集的土壤样品弃去杂草、砾石、动植物残体等杂物，混匀风干后用研钵磨碎过100目(孔径相当于0.15 mm)筛后保存. 采用HCl-HNO3-HF微波密闭消解技术进行土壤样品消解，采用电感耦合等离子体原子发射光谱法(ICP-AES)测定土壤中w(Cd)、w(Cr)、w(Pb)，采用原子荧光法测定土壤中w(As)、w(Hg). 研究区采样点及主要潜在污染源分布见图1.

图1 苏仙区采样点及主要潜在污染源分布概况

1.3 数据处理

研究区土壤中5种重金属元素的描述性统计分析、正态分布检验和转换、Pearson相关性分析，借助SPSS 19.0软件完成；样品污染源解析借助PMF 5.0完成；空间自相关分析借助GeoDa软件完成；克里金插值借助ArcCIS 10.3完成.

1.4 研究方法

1.4.1正定矩阵因子分解(PMF)模型

PMF模型是由Paatero等[19]提出的一种因子分析受体模型. 该模型将采样数据矩阵(X)分解成因子贡献矩阵(G)、因子成分矩阵(F)以及残差矩阵(E)[20-21]：

(1)

式中，a为受体样品个数，b为所测的化学物质种类，p为主因子数(即主要源个数).

PMF模型基于加权最小二乘法进行限定和迭代计算，利用样品的重金属浓度和不确定度数据进行各样点的加权计算，使得目标函数Q最小化. 目标函数Q定义[20]如下：

式中，zcd为第c(c=1,2,…,a)个样品中第d(d=1,2,…,b)个元素的含量，gck为源k对第c个样品的相对贡献，ucd为第c个样品中第d个元素含量的不确定性大小，fkd为源k中第d个元素的含量，ecd为残差.

1.4.2莫兰指数

该研究利用莫兰指数作为测度指标，探讨属性值之间是否具有特殊的空间形态[22-23]，分为单变量莫兰指数和双变量莫兰指数. 其中，单变量莫兰指数可以指出区域同一属性值的分布是聚集、离散或者随机模式[24-25]，双变量莫兰指数揭示了空间中某一要素的一个指标与其相邻位置要素的另一个指标的依赖关系[18]. 二者计算方法[26-27]分别见式(3)(4).

(3)

(4)

式中，I为单变量莫兰指数，为双变量莫兰指数，xi、xj分别为要素i、j的属性值，为属性值的平均值，是第二要素的平均值，wi,j分别为要素i和j之间的空间权重，n为要素总数.

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

(5)

(6)

VarI=E(I2)-E2(I)

(7)

式中，E(I)为莫兰指数的期望值，E(I2)为莫兰指数方差的期望值. 当Z得分大于1.96或小于-1.96时，说明该要素在95%置信区间内呈现明显的聚集或离散特征；当Z得分介于-1.96～1.96之间时，则说明该要素在95%置信区间内呈随机分布.

2 结果与讨论

2.1 土壤重金属含量描述

研究区土壤中5种重金属元素含量的描述性统计结果(见表1)显示，w(Cd)、w(Hg)、w(As)、w(Pb)的平均值均超过了湖南省土壤环境背景值，分别是湖南省土壤环境背景值的15.0、3.2、4.6、7.3倍，表明这4种重金属元素含量的积累受人为活动的影响，其中，w(Cd)、w(As)、w(Pb)的平均值均超过GB 15618—2018《土壤环境质量 农用地土壤污染风险管控标准(试行)》农用地土壤污染风险筛选值，说明可能存在土壤污染风险；w(Cr)的平均值和中值均未超过湖南省土壤环境背景值，表明Cr未受到人为活动影响或影响较小. 该研究区5种重金属元素含量的变异系数差异较大，数据离散程度较高，表明5种重金属元素含量的空间差异比较明显.

表1 研究区土壤中各重金属元素含量

2.2 相关性分析

将5种重金属元素含量正态化处理后，采用SPSS 19.0对苏仙区土壤中5种重金属元素含量进行Pearson相关性分析，结果(见表2)显示，w(Cd)、w(As)、w(Pb)两两之间均呈显著正相关(P<0.001)，且Pearson相关系数均在0.7以上，表明这3种元素间同源性很强.w(Hg)、w(Pb)之间呈显著正相关(P<0.001)，但Pearson相关系数较小，说明这两种元素亦存在部分同源性.w(As)、w(Cr)之间呈显著负相关(P<0.05)，表明两种元素间不存在同源性.

表2 各重金属元素含量间的Pearson相关系数

2.3 土壤重金属的来源解析

利用PMF 5.0对302个土壤样品中5种重金属元素的来源进行解析，并给出各因子的贡献率. 初步将各重金属元素数据载入后信噪比(S/N)均大于3，定义为“strong”. 运行模型后，多次调试重金属元素的“strong”“weak”以及改变因子个数进行多次迭代计算，最终确定4个因子数，重金属元素含量实测值与模拟值的拟合系数均大于0.75，实测Q值与理论Q值的偏差小于10%，解析结果较好.

图2 研究区土壤重金属污染源贡献率

由图2可见，因子1对5种重金属元素的贡献率各不相同，其中对Cr的贡献率最高，达到88.2%. 该研究区w(Cr)平均值低于湖南省土壤环境背景值，受人为活动影响较小，因此，将因子1识别为自然源，受成土母质影响. 因子2对Hg的贡献率最高，为80.6%，其次是Pb，贡献率为24.2%，二者含量相关性较强，表明二者存在部分同源性. 有色金属冶炼、燃煤、金矿和汞矿活动是我国最主要的汞排放源[29]. Hg和Pb的环境迁徙性较强，多项研究指出大气沉降是土壤中Hg、Pb元素的主要来源[30-32]，因此，将因子2识别为工业源中的大气干湿沉降源. 因子3对Cd和Pb两种元素的贡献率较高，分别为63.2%和65.9%；Cd和Pb含量相关性较强，显示同源性较强. 金属矿山的开采、冶炼、重金属尾矿、冶炼废渣和矿渣堆放等是土壤中Cd、Pb污染的主要来源[33-35]. 因此，将因子3识别为工业生产过程中直接排放的工业废弃物. 因子4对As的贡献率最大，为75.1%. 人为As来源众多，且与Cd、Pb、Hg等元素来源相似，含砷矿石的开采、运输、加工等各环节都有损耗；另外，砷化合物作为原料的玻璃、颜料、原药、纸张的生产以及煤的燃烧等过程，都可产生含砷废水、废气和废渣[36-37]，因此，将因子4识别为工业混合源.

综上，该研究区表层土壤中5种重金属元素的来源为自然源和3类工业源，其中，自然源的综合贡献率为41.60%，因子2、3、4三类工业源的综合贡献率分别为22.05%、20.05%和16.30%.

图3 主要污染源贡献率的空间分布

2.4 污染源贡献率的空间分布

为探究各污染源的主要影响区域，将PMF模型解析得到的各源对单个样品的贡献率在ArcGIS 10.3下用普通克里金插值，得到各源贡献率的空间分布情况(见图3). 因子1代表自然源，该源贡献率的高值点主要分布在苏仙区北部，源贡献率均超过52.75%，低值点主要分布在苏仙区中部，说明苏仙区北部受人为活动影响较小，中部受人为活动影响较大，污染程度较高. 整体来看，因子1对苏仙区大部分区域的土壤样品都有着较大的贡献，贡献率大多在32.39%以上，说明因子1作为自然源对土壤样品的影响具有普遍性. 因子2代表工业源中的大气干湿沉降源，该源贡献率的高值点主要分布在苏仙区中东部，贡献率大多超过36.84%. 该高值区主要地形为山地，西部和北部邻近区域分布大量的矿山开采和重金属冶炼企业，地理位置特殊；整体来看，该污染源对大部分地区的影响较小，贡献率大多低于24.62%，说明因子2作为迁徙源具有区域聚集性. 因子3代表直接排放的工业废弃物，该源贡献率的高值点主要分布在苏仙区中部及南部少部分地区，说明这些区域受工矿企业影响较大，人为污染较为严重. 因子4贡献率的高值点主要分布在苏仙区中部及北部，整体来看，该源对大部分区域影响较小，贡献率较低，大多在26.43%以下.

图4 重点污染行业核密度分布

2.5 因子贡献率与企业密度的空间相关性

按照GB/T 4754—2017《国民经济行业分类》将研究区内342家污染企业分为四大类，其中采选业(08黑色金属矿采选业、09有色金属矿采选业)215个，冶炼和压延加工业(31黑色金属冶炼和压延加工业、32有色金属冶炼和压延加工业)105个，化学原料和化学制品制造业(26化学原料和化学制品制造业)12个，其他行业(废弃资源综合利用业、仓储业、生态保护和环境治理业等)10个. 利用ArcGIS 10.3绘制全部企业及不同行业的核密度图，分析各污染源与不同行业间的空间关系，明晰3类工业源的具体类别. 由图4可见，全部企业密度高值区主要分布在苏仙区中部及西南少部分地区，其中采选业主要聚集在苏仙区中部，冶炼和压延加工业主要聚集在苏仙区中部、东北部及西南少部分地区，化学原料和化学制品制造业主要聚集在中北部，其他行业主要聚集在苏仙区中部少部分地区.

运用GeoDa软件分别计算重点污染企业及各因子贡献率的单变量莫兰指数，结果(见表3)显示，在P<0.001显著性水平下，研究区重点污染企业及各因子贡献率的单变量莫兰指数均为正，且数值较大，Z得分均大于1.96，说明不同行业及各因子贡献率在区域内的分布具有聚集性特征. 单变量莫兰指数分析结果均通过显著性检验，具有统计学意义，可为后续双变量检验提供理论支撑.

为探讨各源贡献率与企业间的空间依赖关系，将各源贡献率的空间分布结果与各行业的核密度结果进行双变量莫兰指数分析，结果(见表4)显示，因子1在P<0.001显著性水平下，与全部企业、采选业、化学原料和化学制品制造业的双变量莫兰指数为负，表明三者对因子1的贡献存在负影响. 在重点污染企业分布密集区域，工业源贡献率较高，自然源贡献率较低. 因子1与冶炼和压延加工业、其他行业间未通过显著性水平(P>0.1)，因其各自单变量莫兰指数均通过显著性检验，说明要素本身在空间上呈现聚集特征，但因子1与二者间却不存在空间上的依赖关系，呈空间随机性. 所得结果进一步验证了因子1作为自然源的结论.

表3 研究区重点污染企业及污染源贡献率的单变量莫兰指数

表4 研究区重点污染企业与污染源贡献率间的双变量莫兰指数

因子2在P<0.001显著性水平下，与全部企业和4类重点污染行业间的双变量莫兰指数均为负，表明在4类企业分布较密集的区域内，因子2的贡献率多为低值，说明因子2为非直接工业源. 因子2对Hg元素的贡献最大，气态单质汞(GEM)、气态氧化汞(GOM)和颗粒态汞(PBM)可以在大气中传输，尤其气态单质汞可随大气运动进行跨界传输[38-39]. 宋正城[40]指出，有色金属冶炼活动对周边表层土壤中汞的影响最严重的区域并非位于冶炼厂区，而是位于冶炼厂数公里外，说明汞的积累具有空间迁移性. 所得结果进一步验证了因子2作为非直接工业源中的大气干湿沉降源的结论. 对于迁徙源，污染源与其影响区域不存在对应的空间位置关系，难以识别到具体行业中，其详细成因还需进一步确认.

因子3在P<0.001显著性水平下，与全部企业和采选业的双变量莫兰指数较大，表明因子3与重点污染企业中的采选业在空间上存在显著正相关，在采选业分布密集的区域内，因子3的贡献率亦多为高值；因子3与其余3类行业均不具有正相关关系或相关性不明显，所得结果进一步明确了因子3为采选业排放的工业废弃物，该污染源与其贡献率的高值区存在对应的空间位置关系，因此该类废弃物指在空间上未经远距离主动或被动运输转移，直接对排放口周围的土壤造成污染.

因子4在P<0.001显著性水平下，与4类重点污染行业间的双变量莫兰指数均为正，其中与化学原料和化学制品制造业的双变量莫兰指数最大，表明4类重点污染行业对因子4的贡献率均有正向影响，在重点污染企业分布密集的区域内，因子4的贡献率亦多为高值，进一步说明因子4为各行业的工业混合源，且以化学原料和化学制品制造业排放的工业废弃物为主.

该研究首先对研究区土壤中5种重金属元素进行污染源大类的定性识别，并计算出每个污染源的贡献率，然后依据污染源贡献率与企业密度间的空间关系，将工业源细分到具体的重点污染行业上，对于存在相似贡献的工业源，亦有很好地表征和识别. 秦治恒等[41]利用空间自相关方法，发现湘江流域中土壤Cd含量分布与污染企业的位置存在很强的相关性. 于靖靖等[18]基于莫兰指数和广义加性模型发现包含苏仙区在内的湘江子流域重点污染企业影响区中，采选业对土壤中Cd的影响最大. 这两项研究都利用了空间自相关理论对特定污染源的定性识别，但尚未对污染源的贡献量进行探讨，陈雅丽等[42]统计了近10年关于土壤重金属污染源解析的研究成果，统计结果指出，对湖南省而言，工业活动及其产生的大气沉降是土壤中重金属累积的主要来源. 其中，在矿产业发达的湖南省，Cd主要受采矿业的影响，Hg主要受工业活动产生的大气沉降影响，Pb一方面来自工业活动，另一方面是含铅灰尘的沉降进入表层土，As在工业发达地区受土壤母质和工业活动的影响，Cr主要受土壤母质影响. 笔者所得结果与上述研究结果基本一致. 对比可知，PMF模型与莫兰指数结合的方法，一方面可以定量解析出各污染源的影响程度及主要影响区域；另一方面可以达到明确相关重点污染企业对土壤中5种重金属元素的影响，从而有针对性地提出区域企业管理、重点污染源监管及土壤环境治理策略.

3 结论

a) 研究区5种重金属元素受4类源影响. 其中，自然源是Cr的主要来源，对Cr的贡献率为88.2%；大气干湿沉降是Hg的主要来源，对Hg的贡献率为80.6%；采选业排放的工业废弃物是Cd和Pb的主要来源，对Cd和Pb的贡献率分别为63.2%、65.9%；以化学原料和化学制品制造业排放的废弃物为主的工业混合源是As的主要来源，对As的贡献率为75.1%.

b) 研究区中部受重点污染企业影响严重，且主要以采选业排放的废弃物为主；此外，中部偏北地区主要以化学原料和化学制品制造业废弃物为主；中部偏东地区受大气干湿沉降影响严重，该地区连接3条支流上游，应警惕河流输送和污水灌溉等活动引起的面源污染；西北部地区受重点污染企业影响较小.

c) PMF模型与莫兰指数相结合的方法，使解析结果从工业源大类明晰到不同重点污染行业上，不仅能更好地识别有相似贡献的工业源，而且可以辅助优化验证PMF的有效性.