基于PMF模型的湖北随县三里岗地区土壤重金属源解析

董健彪，谢淑云，田 欢，杨文兵，李 华，周徐维，Diego Armando Pinzon Nunez

(1.中国地质大学(武汉)地球科学学院，湖北 武汉 430074； 2.中国地质大学(武汉)材料与化学学院，湖北 武汉 430074； 3.湖北省地质局第八地质大队，湖北 襄阳 441002； 4.青海省地质调查局，青海 西宁 810000)

1 引言

近年来，大量重金属元素通过采矿、工业废料排放、农药、化肥等途径进入土壤系统，同时由于其在环境中的持久性、生物累积性、难降解性及其与有益元素结合造成新的风险等特点，受到人们的广泛关注(杨希 等，2021；Zhang et al.，2021)。2014年《全国土壤污染状况调查公报》显示，我国土壤重金属污染情况研究也受到广泛关注，点位超标率达到了16.1%，其中As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种重金属是最主要的无机污染物。

土壤重金属源解析工作是污染评价及治理工作的前提，具有重要的现实意义(Hu et al.，2018)。学界流行的源解析方法很多，如主成分分析法(PCA)、正定矩阵因子分解法(Positive Matrix Factorization)等(董騄睿 等，2015)。PMF模型是美国环保署推荐的源解析方法，起初多用于研究大气以及水污染，该方法在识别污染源贡献比例上表现良好，且在与GIS技术结合时可以更有效地识别污染源(Dong et al.，2019；Guan et al.，2018；Liang et al.，2017)，因此近年来越来越多的学者利用PMF模型识别土壤重金属污染来源(Hu et al.，2018)。

三里岗镇地处湖北省随州市，是湖北省香菇主产区之一，香菇产品远销海外，近年出口香菇中出现重金属超标产品，对于研究区农产品销量造成一定影响(杨定国 等，2014)，因此有必要对研究区展开土壤地球化学调查，探讨土壤中重金属元素的含量。本文运用PMF模型结合GIS方法，以三里岗镇为研究区，对土壤样品中的As、Cd、Cr、Cu、Ni、Pb、Zn和Hg元素的含量和来源进行评价和解析，并结合GIS方法分析不同污染源的源贡献率分布情况，以期为土壤重金属污染的治理和管控提供理论依据和技术支持。

2 材料与方法

2.1 区域地质概况

图1 三里岗镇地质图(改自陈婕 等(2018))Fig.1 Geological map of Sanligang town(modified after Chen Jie et al.2018)1—第四纪；2—白垩纪；3—泥盆纪；4—志留纪；5—奥陶纪-志留纪；6—寒武纪；7— 震旦纪；8—元古代；9—一级构造单元界限；10—二级构造单元界限；11—三级构造单元界限

2.2 数据来源与处理

本文数据来源于湖北省“金土地”工程-高标准基本农田地球化学调查(项目编号HBJTD20140108)，该项目由湖北省地质局第八地质大队主持完成，完成了对随县三里岗镇的1∶50000表层土壤样测量，平均采样密度为 5.19 个点/km2。共采集表层土壤1696件，其中重复样34件。

本文使用Excel 2019和SPSS 21进行数据整理、计算和分析，利用GS+(Version7)进行半变异函数分析，通过EPA PMF 5.0进行重金属源解析，利用ArcGIS 10.2进行空间分析和图件制作。

2.3 污染评价方法

本文利用单因子指数法和内梅罗指数法对研究区重金属污染情况进行评价(杨涛毅 等，2011；Meng et al.，2011)。单因子指数法是将土壤中某污染元素的含量值Ci除以评价标准值Si，根据结果Pi对土壤中单一污染元素进行污染评价。公式如下：

Pi=Ci/Si

(1)

内梅罗综合污染指数法利用多种元素的单因子污染指数的平均值(Pivag)以及最大值(Pimax)进行计算，兼顾了多元素的综合特征以及单一因素的极端影响。公式如下：

(2)

2.4 半变异函数

半变异函数是地质统计学中用来描述区域化变量性质的一种常用方法，通过半变异函数，可以对区域化变量的随机性特征以及结构性特征进行分析，得到研究变量的空间分布规律(耿治鹏 等，2022)。半变异函数公式为

(3)

式中，y(h)为半变异函数值，N(h)表示距离向量为h时的采样点数目，h表示空间距离，Z(x)为区域化变量，一般表示为

Z(x)=(xw，xv，xw)

(4)

2.5 PMF模型

PMF全称正定因子矩阵分解法(Positive Matrix Factorization)，是由Paatero et al.(1994)基于因子分析原理所提出的一种分析模型，利用样品的浓度和不确定度，通过最小二乘法计算因子数(徐源，2021)。

PMF模型将重金属元素含量视作i×j阶的矩阵，将该矩阵拆分为三个因子矩阵，分别是因子贡献矩阵G(i×k)，因子成分矩阵F(k×j)以及残差矩阵E(i×j)，公式表示为

(5)

其中Xij为原始矩阵第i个样品第j个化学元素，Gik为第k个因子对于第i个样品的因子贡献度，Fjk是第k个因子对第j个元素的贡献浓度，p则表示因子数(薛建龙，2014)。

PMF模型定义了一个目标函数Q：

(6)

式中，c表示样品序数，d表示元素序数，k表示因子的序数，因此zcd表示第c个样品第d个元素的含量值，gck表示第k个因子对第c个样品的贡献度，dk表示第k个因子对第d个元素的贡献浓度，ucd表示第c个样品第d个元素的不确定度，ecd为残差。

不确定度计算公式为

(7)

当元素浓度小于或等于检出限(MDL)时，的不确定度计算公式为

(8)

3 结果与讨论

3.1 土壤重金属元素参数特征

表1 随县土壤重金属元素描述性统计表①Table 1 Descriptive statistical table of heavy metal elements in soil of Sui county

依据GB15618-2018《土壤环境质量 农用地土壤污染风险管控标准(试行)》(以下简称《标准》)中的风险筛选值为评价标准值对研究区水田和其他用地样品分别统计分析，结果见表2。表2显示，两类用地中，Cd元素污染样品最多，水田中超标样品数为32，占比为5.6%，其他用地中超标样品数为379，占比为38.9%，其次是As元素，在两类用地中超标样品数均超过了5%，而Cr、Ni、Zn元素超标样品数均在2%左右， Cu、Pb、Hg元素无超标样品。整体来看，研究区Cd元素超标率远远超过《全国污染调查公报》中的全国点位平均超标率7%，存在一定的风险(吕悦风 等，2019)。

基于单因子指数法与内梅罗指数法对研究区水田与其他用地土壤中重金属进行评价，评价标准值表2已给出，评价结果见表3，污染评价标准见表4(邓通德 等，2022)。从表3中可以看出研究区8种重金属元素的单因子指数均值皆小于1，研究区8种重金属在全区的平均污染程度不高，除Cu、Pb、Hg外其余元素均有污染样品，但只有Cd元素污染样品数超过了20%，形成一定规模。综合污染指数P综均值为0.73，根据表4中评价标准，研究区整体污染程度处于尚清洁等级，有253件P综>1的污染样品，占比达到16.39%，说明研究区污染情况值得重视。

表2 重金属超标样品数目统计②Table 2 Statistics of the number of heavy metal samples exceeding the standard

表3 重金属污染指数描述统计量③Table 3 Descriptive statistics of heavy metal pollution index

续表

表4 污染评价标准Table 4 Standards for pollution assessment

3.2 土壤半变异函数分析

利用GS+软件对数据进行半变异函数分析，并根据决定系数R2越趋于1、残差平方和RSS越趋于0则模型拟合度越好的原则选取合适的模型，半变异函数分析要求数据符合正态分布或近正态分布(吴敏 等，2016；蔡大为 等，2020)，而8种元素原数据都不符合正态分布特征，进行半变异函数分析前对数据进行转换，转换前后的偏峰系数对比见表5。半变异函数分析结果见表6，8种元素的决定系数R2均大于0.7，残差平方和都比较小，说明模型拟合较好。

表5 偏峰系数检验Table 5 Test of skewness coefficient and kurtosis coefficient

表6 半变异函数分析Table 6 Semi-variogram analysis

3.3 土壤重金属PMF源解析

利用EPA PMF 5.0对8种重金属元素进行PMF源解析，根据Qrob/Qexp比值快速下降处因子数最佳(李娇 等，2019)将因子数确定为4。表7列出了4个因子对8种元素的源成分谱和源贡献率。图2是各因子对各元素的源贡献率，利用ArcGIS对PMF源贡献率进行克里格插值，绘制出4个因子的源贡献率分布图(图3)。

表7 PMF源解析出的各源成分谱及其贡献率Table 7 Component spectrum and contribution rate of each source resolved by PMF source

图2 各元素因子源贡献率Fig.2 Factor source contribution rate of each element

图3 因子源贡献率分布图Fig.3 Distribution of factor source contribution rates

根据表7及图2，因子1主要代表As元素的来源，对As元素的源贡献率高达71.78%，对Pb元素也有一定的源贡献率(22.45%)。由于我国煤资源中As元素含量很高(白向飞，2003)，高砷煤燃烧过程中As便会进入环境中(Finkelman et al.，2018)，但煤燃烧同样会导致Hg污染，而因子1对于Hg的源贡献率却很低，因此因子1应代表其他来源。通过观察图1和图3中因子1源贡献率分布图可以看出，因子1高源贡献率与震旦系(Z)以及白垩系(K)地层重合度很高，因此推测因子1可能与地层有较大关联。表8统计了1544组样品、不同年代地层中8种重金属元素的含量，有35组样品属于岩浆岩，由于占比少，未参与统计。从表6中可以看出，震旦系以及白垩系地层中的As含量(15.13 mg/kg和14.79 mg/kg)和Pb(25.70 mg/kg和26.66 mg/kg)含量相比其他地层更高，同时也高于研究区含量均值(As：11.52 mg/kg；Pb：24.42 mg/kg)。因此因子1代表的来源主要是震旦系和白垩系地层母岩。

表8 不同年代地层8种重金属元素含量④Table 8 Contents of eight heavy metals in strata of different ages

从表7及图2中看出，因子2对Hg、Pb两种元素的源贡献率高达76.96%、63.28%，对As、Cd、Zn等元素也表现出了一定的影响，源贡献率分别为19.00%、16.24%、25.80%。有研究表明，交通排放是一个重要的重金属(包括Hg、Pb、Cd等元素)来源(Liu et al.，2012；Yang et al.，2017)，而三里岗乡村道路分布全区，从图3也可以看出，因子2源贡献率超过24.20%的区域在全区分布广泛，因此因子2应该受到了交通排放所带来的影响，但因为研究区为乡镇，车辆流通量小，所以交通排放到土壤中的重金属含量并不高，所以并未有样品的Hg、Pb含量超过风险筛选值。其次，图3中因子2高源贡献率区域(>34.34%)分布在北部、中部以及南部，与图1元古界(Pt)地层重合程度较高，说明因子2可能也受到了来自元古界地层母岩的影响。半变异函数分析表明Pb受到自然因素和人为因素共同影响，因此因子2解释为元古界地层母岩与道路交通排放的共同影响较为合理。

从表7以及图2中看出，因子4对于Cr、Cu、Ni的影响最高，分别达到了88.40%、82.14%和88.77%，同时对Zn元素的贡献率也达到了42.46%。从图1及图3中因子4源贡献率分布可以看出，因子4高源贡献率(>39.87%)区域与奥陶系到志留系(O-S)地层高度重合，从表8中也可以看出，Cr、Cu、Ni在寒武系到志留系地层中的含量远高于其余地层中的含量，Zn元素含量也高于研究区背景值，所以因子4可以解释为奥陶系到志留系地层母岩的影响。

4 结论

本文利用半变异函数及PMF模型对三里岗镇8种重金属元素的含量及来源进行了分析，结果表明：

(1)8种重金属元素除Hg、Pb外的均值都高于全国背景值，依据《标准》中所给出的风险筛选值进行评价，Cd元素超标率超过20%，内梅罗综合污染指数法显示三里岗镇污染样品占比为16.39%，平均污染指数为0.73，处于尚清洁水平。

(2)半变异函数分析结果表明，8种重金属元素的块金效应均<50%，表明研究区重金属元素含量主要受自然因素影响。

(3)PMF源解析得到了4个来源，依据不同地层中重金属元素的含量以及因子源贡献率分布情况将因子1解释为震旦系和白垩系地层母岩的影响、因子2解释为元古界地层母岩与乡村道路交通排放的混合源、因子3解释为震旦系和寒武系地层母岩的影响、因子4解释为奥陶系到志留系地层母岩的影响。

致谢：本研究数据主要来自湖北省地质局第八地质大队和“金土地”工程-高标准基本农田地球化学调查课题(项目编号HBJTD20140108)，受中国地质大学(武汉)中央高校基金生态地球化学团队项目(CUG170104)联合资助。