吕柏楠, 王 超, 师华定, 李明秋

1.河南理工大学, 河南 焦作 454003 2.生态环境部土壤与农业农村生态环境监管技术中心, 北京 100012

土壤是自然生态系统的重要组成部分，是人类赖以生存和农业生产的重要物质基础[1]. 近年来，我国环境问题日益严重，工矿、农业土壤环境状况不容乐观[2]. 2014年全国土壤污染状况调查显示，我国耕地土壤污染点位超标率达19.4%，主要污染物为重金属，其中镉污染点位超标率达到7.0%，其次分别为镍、砷、铜、汞、铅、铬等[3]；据2016年原农业部监测系统调查数据，我国每年因为土壤重金属污染而使粮食产量降低1.0×107t以上[4]. 重金属污染物在土壤中不断富集，不仅会造成耕地质量下降、农产品产量降低等负面影响，还会通过食物链进入人体，威胁人类健康[5]. 因此，明晰重金属污染物的来源，是进行土壤重金属污染防治的重要前提.

空间插值技术可以通过合理布设土壤点位，实现研究区土壤污染数据的全覆盖. 因此，如何进行空间模拟预测及方法的选用是准确预测区域土壤属性空间分布特征的关键[6]. 克里格(Kriging)插值为土壤属性的空间预测提供了一种无偏最优估值方法，是目前地统计学中应用最广的方法之一[7-8]. 作为克里格的衍生方法，指示克里格是一种不依赖于正态分布假设条件的非参数估值方法，可以处理高程度的空间变异现象而不必削减重要的高值数据，适合用于解决存在多特异值的土壤污染数据[9-10]. 段艺芳等[11]采用指示克里格法，发现陕北地区8种土壤重金属生态风险在空间上均表现出由南向北、由东南向西北逐渐增加的趋势；杨奇勇等[12]采用指示克里格法，发现云南省广南县除Cu元素外，其他重金属元素高概率污染区域都集中在岩溶区.

污染源解析技术是针对污染物的来源进行定性(源识别)或定量(源解析)判断的方法[13]. 目前，受体模型法是当前土壤污染物源解析研究中最主要的技术方法. 受体模型种类较多，常用的有化学质量平衡(CMB)模型[14]、主成分分析/因子分析-多元线性回归(PCA/FA-MLR)模型[15-16]、正定矩阵因子分解(PMF)模型[17]和UNMIX模型[18]等. 其中，主成分分析/因子分析-多元线性回归(PCA/FA-MLR)模型在土壤重金属源解析领域应用较早，它不需要准确的源成分谱数据，只需要知道受体信息及排放源组成，使用较为广泛[19]. Li等[20]采用FA-MLR模型，揭示了重庆4个典型农业县的工业排放、燃料消耗和当地农业对8种重金属的主要人为贡献为64%～90%；Yang等[21]利用PCA-MLR模型，揭示了湖北省青山区东部地区2010—2014年期间增加的重金属浓度有82.5%来自工业/农业活动、17.5%来自交通源.

云南是矿产资源大省，素有“有色金属王国”之称，矿产资源矿种齐全、储量丰富，金子河流域铜矿、铅锌矿、金矿等矿产资源丰富，开发历史悠久，黑色金属矿和有色金属矿采选企业数量较多. 因此，该研究以位于金子河流域中游的铜厂乡耕地为研究区域，运用地统计学方法和主成分分析-多元线性回归(PCA-MLR)模型对研究区土壤中的Cd、Hg、As、Pb、Cr 5种重金属的空间分布及来源进行解析，探寻金子河两岸农用地土壤重金属污染成因，以期为该区域土壤污染防控提供科学依据.

1 材料与方法

1.1 研究区概况

研究区位于云南省红河哈尼族彝族自治州金平县西部(22°43′N～22°52′N、102°55′E～103°03′E)，如图1所示，地势南高北低，最高海拔 1 191 m，最低海拔413 m，年平均气温13.6 ℃，年均降雨量800～1 000 mm，主导风向为西南风；土壤类型以黄棕壤为主，成土母岩以灰岩和玄武岩为主；研究区西岸区域以旱地为主，东岸区域以水田为主；主要重金属污染企业有黑色金属(铁、铜)矿采选厂、有色金属(金)矿采选厂以及有色金属矿山尾矿库，主要分布在研究区的北部和西部区域.

图1 金子河流域矿区分布及研究区采样点分布Fig.1 Distribution of mining areas in the Jinzi River Basin and distribution of sampling points in the study area

1.2 样品的采集与分析

土壤采样点采用网格布点法进行布设[22]，网格大小为250 m×250 m，采样时用GPS进行精确定位，采集表层(0～20 cm)土壤样品76个，并选取其中8个点位采集深层(100～150 cm)土壤样品. 采用梅花形布点法采集4～6个子样，每个子样均使用竹铲垂直采集，去除杂草、砾石、肥料块等杂质，等份组合成一个混合样. 室内将样品进一步用球磨机碾碎过200目尼龙筛，土壤重金属全量采用HNO3-HCl-HF-HClO4溶样，采用电感耦合等离子体原子发射光谱法(ICP-AES)测定土壤中Cd、Cr、Pb的含量，采用原子荧光法测定土壤中的As、Hg的含量[23].

1.3 数据处理

研究区土壤中5种重金属元素的描述性统计分析与样品污染源解析均借助SPSS 19.0软件完成；半方差函数借助GS+9.0软件完成；指示克里格插值借助ArcGIS 10.2软件完成.

1.4 研究方法

1.4.1内梅罗综合污染指数法

内梅罗综合污染指数(Nemerow integrated pollution index，)由美国学者Nemerow[24]提出，是在综合各因子平均污染状况的基础上，兼顾极值和均值，有效评价研究区各采样点重金属的综合污染状况，是目前进行综合污染指数计算最常用的方法之一[25]，计算公式：

(1)

式中，PN为内梅罗综合污染指数，Pimax为采样点重金属i污染物单项污染指数中的最大值，Piave为采样点重金属污染物单项污染指数的平均值.

内梅罗指数评价法分级标准见表1.

表1 内梅罗综合污染指数评价法分级标准

1.4.2指示克里格法

指示克里格(indicator kriging)方法是由Journel提出的，被广泛应用于具有极端值和高偏离度特性的土壤污染研究[26-29]. 指示克里格预测值的大小代表大于或小于特定阈值的概率，即从指示数据得出的非抽样地点的期望值等于变量的累计分布[30]. 指示克里格具体计算方法[31]：①根据指标的阈值确定指示器；②计算最佳指示半方差函数；③利用半方差函数模型参数在ArcGIS 10.2软件中进行指示克里格插值，得到土壤重金属污染概率的空间分布图.

1.4.3PCA-MLR模型

主成分分析-多元线性回归法(PCA-MLR)首先利用主成分分析法定性识别污染源，在获得污染因子数量后，对各示踪元素的浓度进行多元线性回归(MLR)，就可得到重金属浓度与各示踪元素浓度的回归式，回归系数用于计算各源的相对贡献率[32-33].

2 结果与讨论

2.1 土壤重金属含量描述

由土壤重金属含量的描述性统计分析结果(见表2)可知，研究区土壤中Cd、Hg、As、Pb、Cr的平均含量分别为3.17、0.16、44.34、779.57和358.85 mg/kg，均高于云南省土壤背景值. 其中，Pb的平均含量(779.57 mg/kg)远高于云南省土壤背景值(40.60 mg/kg)，前者约是后者的19.20倍；其次，Cd的平均含量(3.17 mg/kg)约是云南土壤背景值(0.22 mg/kg)的14.55倍，表明这2种重金属在土壤中的富集程度较高，可能受到的人为活动较剧烈. 依据GB 15618—2018《土壤环境质量 农用地土壤污染风险管控标准(试行)》可知，研究区土壤中Cd、Hg、As、Pb、Cr含量低于筛选值的点位比例分别为26.32%、100%、56.58%、39.47%、17.11%，超过管制值的点位比例分别为40.79%、0、5.26%、40.79%、2.63%，说明研究区Cd、Pb污染较为严重. 从变异系数看，Pb、Cd的变异系数最大，分别为148%、139%，属于高程度变异；其次为As(95%)、Hg(76%)、Cr(55%)，为中等程度变异. 其中，Pb、Cd在研究区分布极不平衡，表明土壤中Pb、Cd的含量受到了较为明显的人为影响.

表2 样品描述性统计

2.2 内梅罗综合污染指数评价分析

内梅罗综合污染指数反映的研究区土壤点位重金属的综合污染情况如表3所示. 研究区综合污染指数范围为1.73～131.69，平均值为18.28，处于重度污染程度的土壤点位比例为90.79%，说明研究区土壤整体处于重度污染水平. 其中，土壤中Cd、Hg、As、Pb、Cr的单因子污染指数范围分别为0.14～96.75、0.62～13.50、0.07～7.07、0.23～178.77、1.45～19.59，而土壤点位中Cd、Pb、Cr属于重度污染的比例分别为63.16%、59.21%和78.95%，表明研究区Cd、Pb和Cr污染比较严重.

表3 土壤重金属内梅罗综合污染指数

2.3 空间分布特征

选择GB 15618—2018的筛选值作为5种重金属元素的阈值，研究区土壤样点中的Hg含量均低于筛选值，对农产品质量安全、农作物生长或土壤生态环境的风险低，因此不进行空间插值分析.

2.3.1半方差函数分析

对各元素变量进行半方差函数模型拟合，分析局部地区变量的空间变异及相关程度. 拟合结果显示：Cd、As、Pb这3种重金属的决定系数均大于0.9，所选拟合模型均符合地统计要求，而Cr的决定系数仅为0.09，说明Cr元素在研究区范围内不具有空间相关性，不能进行半方差函数分析，不适宜空间局部插值. 块基比能够反映随机因素导致的空间变异占总体变异的比例[35]，当块基比值小于0.25时表示变量具有显著的空间自相关，块基比值范围在0.25～0.75时表示变量具有中等自相关，块基比值大于0.75时表示变量自相关程度不强[36]. 变程大小反映了当地重金属空间异质性的尺度或空间自相关尺度[37].

基于GS+ 9.0软件，获得研究区Cd、As、Pb的最佳指示半方差函数模型(见图2)及参数(见表4). Cd、As的半方差函数符合球形模型，Pb符合高斯模型. Cd、As、Pb的块基比分别为0.192、0.113、0.066，均在0.25以下，都表现为显著的空间自相关性，说明在当前阈值条件下，Cd、As、Pb条件概率的空间分布是由结构性因素(如土壤母质、地形、土壤类型等)共同作用引起的. As的变程为 2 523.88 m，Cd的变程为 1 723.00 m，Pb的变程为 1 190.04 m. 从变程数据来看，研究区域土壤中Cd、As空间分布的均一性相比Pb要强，而Pb在土壤中的分布存在小范围的变异现象[38].

表4 土壤重金属含量的理论模型及相关参数

图2 指示半方差函数模型Fig.2 Indicating semivariance function model

2.3.2指示克里格分析

根据Cd、As、Pb的指示半方差函数模型，利用ArcGIS 10.2软件进行克里格插值，得到研究区耕地土壤Cd含量大于0.3 mg/kg、As含量大于40 mg/kg、Pb含量大于90 mg/kg的污染概率分布结果(见图3).

从图3可以看出，研究区西岸土壤中Cd、As、Pb的污染风险较高，特别是Cd和Pb，在绝大部分区域的污染风险概率高于0.7，且Cd和Pb的污染风险概率分布相似，可能存在相同的污染源；而As仅在北部区域的污染风险概率低于0.3，在中南部绝大部分区域的污染风险概率高于0.7，而且与Cd和Pb的污染风险概率分布相似，三者可能存在相同的污染源. 而西岸地区属于海拔为500～1 000 m的山地区域，富含矿产资源，且附近存在多家金属矿选企业. 因此，推断西岸重金属的主要污染来自矿产资源的开发及金属冶炼所产生的废水、废气、废渣[39-41].

图3 表层土壤Cd、As和Pb污染概率分布Fig.3 Pollution probability distribution of Cd, As and Pb in top layer soils

研究区的东岸地区，土壤Cd、As、Pb的污染风险较低，其中，Cd的污染风险概率仅在北部区域高于0.7，中南部绝大部分区域低于0.3；As的污染风险概率在东岸绝大部分区域低于0.3；Pb仅在北部小部分区域存在概率高于0.7的污染风险，其余大部分区域的污染风险概率低于0.3. 而东岸附近无企业，因此，推断东岸重金属的主要污染来自土壤母质或土壤类型等自然源.

总体而言，Cd、As、Pb污染的高概率区域主要分布在研究区的西部与西南部地区，北部地区是Cd、Pb污染的高概率区域，而Cd、As、Pb污染的低概率区域主要分布在研究区的东部以及东南部地区，Cd、As、Pb污染风险分布具有显著的地域性差异.

2.4 土壤重金属的来源解析

根据研究区重金属元素的空间分布特征，将研究区的东西两岸(金子河为界)及研究区整体分别利用PCA-MLR模型进行来源解析，探究局部(东西两岸)与整体(研究区)污染源的关系.

2.4.1源类型识别结果

由表5可见，对于研究区整体而言，PCA-MLR源1主要载荷元素为Pb、Cd、As和Hg. 金属矿山的开采冶炼、工业废渣和矿渣堆放等是土壤中Pb、Cd污染的主要来源[42-43]；As常见于多种金属矿床中，尤其多伴生于金矿床[44]；金属冶炼工艺是环境中汞元素的主要来源[45]. 而依据实地调查可知，研究区的北部地区存在有色金属矿(金矿)采选厂及尾矿库，西部地区存在黑色金属矿(铁矿、铜矿)采选厂和有色金属(金矿、铅锌矿)采选厂. 因此，将因子1识别为以选冶为主的工业源. 对于东岸和西岸而言，PCA-MLR源1主要载荷元素均为Cd、As和Pb. 其中，东岸的源成分谱中As占比较大，因此，将东岸因子1识别为以金矿选冶为主的工业源；西岸的源成分谱中As所占比重较小，因此，将西岸因子1识别为以黑色金属矿(铁矿、铜矿)和铅锌矿选冶为主的工业源.

表5 土壤样品的主成分因子载荷

对于研究区整体、东岸及西岸而言，PCA-MLR源2主要载荷元素均为Cr、Hg. 该研究区Cr的变异系数较小，受人为活动影响较小. 由于深层土壤中的重金属含量高于表层土壤，可能受土壤成土母质的影响[46]. 由图4可见，研究区大部分表层土壤中Cr的含量低于深层土壤，因此，将因子2识别为自然源，受成土母质影响.

图4 研究区Cr元素表层和深层土浓度对比Fig.4 Contrast of Cr element surface and deep soil concentration in the study area

2.4.2源贡献率结果

由图5可见，研究区整体、东岸和西岸的土壤重金属来源均为2种，分别是工业源和自然源. 对于研究区整体而言，工业源的相对贡献率最大，为87.21%；对于东岸而言，自然源的相对贡献率最大，为93.46%；对于西岸而言，工业源的相对贡献率最大，为91.02%. 因此，成土母质等自然因素是东岸土壤重金属含量的主要来源；工业矿选冶炼是西岸土壤重金属含量的主要来源，也是整个研究区土壤重金属含量的主要来源.

图5 研究区整体、东岸和西岸表土中 重金属各类源的相对贡献率Fig.5 The relative contribution rate of various sources of heavy metals in the surface soils of the whole study area, the east bank and the west bank

造成东西岸区域土壤重金属污染主要来源不同的原因，主要是西岸区域周围较多的矿产资源开采冶炼过程中产生的废水、烟粉尘及冶炼渣、尾矿库尾砂中的残余重金属淋滤迁移，以及关闭或搬迁企业历史遗留重金属废渣常年随雨水冲刷、淋溶，进入周边土壤. Xu等[47]研究发现，矿业对Cd、Pb、As的主要影响发生在周边地区5 km范围内，其中，西岸区域距离矿产开采冶炼场所基本在7 km范围内，而东岸区域距离矿产开采冶炼场所基本超过7 km. 尾矿库作为矿山开采主要污染区域，尾矿库中的尾矿砂随风力迁移，及废弃矿山构筑物淋滤水力作用驱动，而影响周围环境[48]. 西岸区域海拔高且落差大，尾矿砂中的重金属受降水等水力作用影响范围较广，而东岸距离尾矿库较远，地势较为平缓，受到来自尾矿库的影响较小. 因此，要控制研究区耕地周围重点污染企业“三废”的排放，加强对关停或搬迁的遗留工业场地周边污染的监测，才能有效防范工业产生的重金属污染对周边耕地的影响.

3 结论

a) 金子河流域土壤中Cd、Hg、As、Pb、Cr含量低于GB 15618—2018《土壤环境质量 农用地土壤污染风险管控标准(试行)》的筛选值的点位比例分别为26.32%、100%、56.58%、39.47%、17.11%，含量超过管制值的点位比例分别为40.79%、0%、5.26%、40.79%、2.63%；研究区域内重金属综合污染指数较高，处于重度污染程度的土壤点位比例为90.79%，且Cd、Pb和Cr的污染比较严重.

b) 从空间分布特征可以看出，元素Cd、As、Pb污染的高概率区域主要分布在金子河流域西部与西南部地区，Cd、Pb污染的高概率区域主要分布在研究区的北部地区，Cd、As、Pb污染风险分布具有显著的地域性差异.

c) 金子河流域东西两岸表土中重金属污染的主要来源不同. 成土母质等自然因素是东岸土壤重金属含量的主要来源；工业矿选冶炼是西岸土壤重金属含量的主要来源，也是整个研究区土壤重金属含量的主要来源.