龙口市污灌区重金属污染影响因子分析*

刘建华 顾振飞 曹见飞，2 徐明雪 王翠秀# 吴泉源

(1.山东师范大学地理与环境学院，山东 济南 250358；2.伟志股份公司，福建 泉州 362200)

面对全球水资源匮乏的现状，污水灌溉已经成为许多地区农田灌溉的主要方式，然而污水中含有的重金属元素将不断在土壤中富集，并会通过农作物进入食物链影响人类健康，由此产生的风险引起学者们的广泛关注[1-2]。对污灌区土壤重金属进行影响因子分析有助于重金属污染的有效治理，对预防重金属污染产生的潜在危害也具有重要的意义。污灌区土壤重金属的含量同时受到多个因子的影响，传统的影响因子研究方法主要有Pearson相关系数法、主成分分析法、地理加权回归等方法[3-5]。这类回归统计方法虽然可以测量和估计各影响因子的影响程度，但这些分析方法都是在假设预测变量独立的情况下进行的，忽略了驱动因素和重金属浓度之间的空间关系，并且受影响因子间多重共线性的阻碍，从而产生较大的方差，降低了模型的推理能力[6]。2010年，WANG等[7]研发了地理探测器模型并用于地理环境因子分析，目前已有多名学者将地理探测器模型引入土壤重金属污染的影响因子解析中，并且取得理想的效果[8-11]。为此，本研究采用地理探测器模型综合分析污灌区重金属污染的影响因子，解决以往研究中影响因子解释度低的问题。

龙口市是山东省重要的老工业基地，同时也是胶东农产品生产基地，该地污水灌溉历史悠久，污水中的有害物质在土壤中不断积累，使得农田生态安全面临严峻的挑战，农田污染现状调查及农田生态影响分析等都是亟待解决的问题。本研究将重金属富集因子特征分析与地理探测器模型分析相结合，评估龙口市污灌区重金属的污染水平与空间分布，分析各影响因子对重金属污染的贡献，探索不同影响因子之间的协同或拮抗作用，在此基础上，比较了地理探测器模型与Pearson相关系数法的优缺点，以期为污灌区重金属空间分布特征与影响因子分析提供一种高速有效的方法，为农田规划与环境治理提供借鉴。

1 材料与方法

1.1 研究区概况

研究区位于山东省龙口市北部平原地区(37°34′35″N～37°44′12″N， 120°20′56″E～120°40′47″E)，区域总面积409 km2。研究区内成土母质以冲积和海积物为主，土壤类型主要包括褐土、棕壤和风砂，pH为6.5～8.0，成土条件基本一致。研究区距海较近，属温带季风性气候，年平均气温12 ℃左右，年均降水量600 mm左右，年降水分配不均，主要集中在5—9月，具有春秋季旱、夏季涝的特点[12]。研究区内大部分区域为农田，主要种植小麦、玉米等作物，核心建成区主要分布着芦头工业区、住宅区、商业区等。研究区内煤炭、金、高岭土、砂矿等矿产资源丰富，煤炭、采矿、铸造和电镀等工业发达，其中工业残渣的堆放、废液排放都会对周围的环境产生一定的影响。由于当地降水难以满足农业需求，经污水处理厂处理后的生活废水、工业污水经由黄水河及永汶河后进入各支渠用于农田灌溉，加剧农田土壤中各类重金属元素的富集，区域内主要污染重金属为Cr、Ni、Pb、Zn、Cu、Co、As、Cd。

1.2 样品采集和理化分析

于2020年5、11月采集土壤样品，采样期间降水稀少，且农作物处于成长的关键期，一般会采用河水灌溉。根据研究区各种土地利用类型占比及主要污染源(大型工矿业)分布情况布设采样点，样品采集深度为0～20 cm，以30 m×30 m为一个样方，清理土壤表面后用小木铲进行采集，四分法混合均匀后剩余约1 kg装入聚乙烯自封袋中低温封装带回。本研究共设70个采样点，皆采用全球定位系统(GPS)定位，采样点分布见图1。首先挑除过滤土壤样品中的粗大杂质，然后使用玛瑙研钵磨细直至通过200目孔径筛。对预处理后的土壤样品进行重金属检测分析，采用石墨炉原子吸收光谱法(GF-AAS)测定Cd含量；用电感耦合等离子体原子发射光谱法(ICP-OES)测定Cu、Pb、Zn、Co、Ni含量；采用氢化物发生-原子荧光光谱法(HG-AFS)测定As含量；采用火焰原子吸收分光光度法测定Cr、Mn含量。每个样品重复测量3次取平均值，用去离子水调节土壤样品土水比为1 g∶5 mL后测定pH。最后，根据国家标准土壤样品(GBW7401)进行准确度和精密度的控制，测量结果均符合质量控制要求。

图1 研究区采样点分布Fig.1 Distribution of sampling points in the study area

1.3 数据来源与统计分析

研究区位于龙口市污灌区，灌溉系统和农业投入品对当地土壤重金属污染均存在一定影响。但研究区范围较小，农产品相对单一，各采样点灌溉习惯和农业投入情况差异较小，因此不将灌溉系统和农业投入品作为影响因子进行研究。本研究主要考察7个影响因子对重金属污染的影响，其中包括3个人为因素(工业密度(DI)、人口密度(DP)和道路密度(DR))和4个自然因素(植被覆盖度(NDVI)、土壤母质(PMT)、土壤有机质(SOM)和高程(DEM))。其中，NDVI、PMT、DEM数据来源于2014年山东省土壤类型分布图；DP数据来源于龙口市统计年鉴；使用ArcGIS 10.0软件采样点为中心建立2 000 m的缓冲区，计算每个缓冲区内的DI和DR数据；SOM为点数据，与70个土壤采样点一一对应，实测获得。采用地理探测器模型进行分析时需要先将数据离散化，因此采用ArcGIS 10.0软件中的自然断裂分类方法将连续数据转换为离散数据，以满足地理探测器模型的数据格式要求。为了评价地理探测器模型的效果，将分析结果与Pearson相关系数法的分析结果进行比较。

1.4 研究方法

1.4.1 富集系数

富集系数可以用来粗略评价环境中重金属的污染程度(富集程度)并分析重金属来源，富集系数计算见式(1):

(1)

式中:EFi为重金属i的富集系数；Ci为土壤样品中重金属i的质量浓度，mg/kg；Cn为土壤样品中参考元素的质量浓度，mg/kg；Bi为重金属i的背景值，mg/kg；Bn为参考元素的背景值，mg/kg。

对研究区各重金属元素的含量特征进行分析发现，Mn的变异性较小，受人为活动影响较小[13]，且含量水平基数与其他元素相比较大，因此本研究中选择Mn作为参考元素。在重金属来源分析时，认为富集系数≤1.5时重金属主要受自然因素影响，当富集系数>1.5时重金属主要受人为因素影响[14]。

1.4.2 地理探测器模型

本研究采用地理探测器模型分析影响因子对研究区重金属污染的影响。地理探测器模型是一种分析地理事物空间分异及其驱动力的空间统计方法，可用于分析因子的交互作用对因变量的影响。地理探测器模型由风险探测器、因子探测器、生态探测器和交互探测器4大模块组成，通过计算单个或多个因子与因变量间的解释力(q)分析因子对因变量的影响程度。q的取值范围为0～1，q=0时表示因子与因变量间无关联，q=1时表示因变量可全部由该因子解释，q越大表明因子对因变量的影响越大。在因子交互探测中，交互探测模块通过比较2个单因子a1、a2交互后的解释力(即q(a1∩a2))与a1、a2单独作用时的解释力(即q(a1)、q(a2))大小，判别两者间是否存在交互作用及其作用强弱，交互探测表达式见表1[15]。

表1 交互探测表达式Table 1 Expression of interactive detection

2 结果与讨论

2.1 污灌区重金属污染的描述性统计与评价

研究区土壤样品pH为7.00～7.85，平均值为7.43，属于中性偏碱性土壤，土壤样品中重金属质量浓度与富集系数统计见表2。Cu、Cr、Ni、Pb、Zn、As、Cd、Co的平均值均高于龙口市土壤背景值，最大值分别为背景值的9.87、1.92、3.17、7.70、3.35、3.25、8.06、1.79倍，说明研究区土壤中重金属受人为干扰较大，集聚趋势比较明显。本研究选取Mn作为参考元素，计算得到研究区Cu、Cr、Ni、Pb、Zn、As、Cd、Co的富集系数平均值分别1.96、1.32、1.49、1.75、1.69、1.56、3.76、1.28。其中，Cu、Pb、Zn、As、Cd的富集系数大于1.5，可以被认为是受到人为因素影响的主要重金属元素，因此在研究区中富集现象明显；Co、Cr、Ni富集系数低于1.5，且在研究区范围内标准偏差均小于等于0.2，说明这3种重金属的富集主要受自然因素的影响较大，受到人为因素影响较小。

表2 研究区重金属质量浓度和富集因子统计Table 2 Statistical of mass concentration of heavy metals and enrichment factor in study area

2.2 重金属富集系数的空间分布

图2为研究区70个采样点8种重金属富集系数的空间分布情况。8种重金属的富集系数分布趋势总体相似，基本表现为西北部富集系数相对偏高，这是因为该地区位于东海工业园区附近(发电厂东部)；此外，富集系数的高值点位于研究区中东部，靠近徐福镇的煤炭开采地区域。8种重金属中，Cd的富集系数普遍偏高，表明该元素污染较为严重，可能是农资产品、地膜以及除草剂在农田中普遍使用造成的[16-17]；Co、Cr的富集系数在全区范围内均偏小，富集程度较轻。总体看来，研究区重金属富集系数的分布情况与工业区的分布具有相似性，主要受到工业园区、煤矿区、农作种植等的影响。

2.3 地理探测器模型分析

2.3.1 因子分析

利用地理探测器模型中的因子探测器模块评估7个影响因子对8种重金属污染的解释力，结果见表3。总体看来，7种影响因子中DI对8种重金属的平均解释力最大,其次为SOM和PMT，说明工业是重金属污染的主要来源。实地调研发现，龙口市许多工厂、煤矿区过去直接向河流系统排放工业废水，以往流域内的发电厂、电镀、冶炼等企业也存在将未经处理或处理不彻底的污水排放至河流及农田的情况，因此研究区内工业分布是造成重金属污染的主要原因[18]。8种重金属元素中，DI对Cd、Pb的解释力较大，说明工业是影响土壤样品Cd、Pb污染的主要影响因子。研究区内煤矿区分布较多，长期以来的煤炭开采、煤矸石的堆放使污染物随雨水进入河流，煤炭废渣以及汽车尾气排放的Cd、Pb经过大气的干湿沉降迁入土壤，这些都是造成附近土壤中Pb、Cd富集的主要原因[19]；此外，农业活动中为提高产量而使用的尿素、磷肥也会使得土壤中Pb、Cd含量增加[20-21]。SOM对Cu和Zn的解释力较大，主要是因为SOM对重金属有较强的络合作用。有研究表明，猪饲料中添加了大量的Cu、Zn来防治疾病、加快生长，但是90%(质量分数)以上会以猪粪的形式排出，且猪粪中有机质丰富，络合作用较强[22-23]，因此SOM是土壤中Cu、Zn富集的主要影响因子。Cr受到PMT、DP、SOM及DI共同作用，说明Cr污染的影响因子较为复杂；Co、As和Ni等主要受PMT的影响，人类活动对这类重金属在土壤中的富集影响较小，这与前人的研究结果总体一致[24]。本研究中DR对各重金属富集的影响处于较低水平，这可能与大多数工业区位于郊区，DR通常较低，这在某种程度上可能掩盖了DR对重金属污染的影响。

图2 研究区8种重金属富集系数的空间分布Fig.2 Spatial distribution of 8 heavy metals enrichment factor in study area

2.3.2 交互作用分析

从图2(g)可以看出，Pb分布的空间异质性较大，在影响因子分析中更具有代表性，因此以Pb为例考察影响因子间对重金属污染的交互作用。表4为影响因子间交互作用对Pb污染的解释力，通过比较单一因子解释力和交互作用解释力的大小关系，判断两个影响因子的交互作用效果。经比较，NDVI、SOM、DI、DR、DP、DEM、PMT间的交互作用均为非线性增强，也即影响因子间单因子解释力合计小于影响因子交互作用的解释力，Pb分布受到各影响因子间交互作用的影响。其中，DI和DP交互作用对Pb的解释力达0.921，是影响Pb分布的主要因子，其他影响因子也会加强工业以及人口因素对Pb分布的影响。

影响因子间交互作用非线性增强，显示了人为因素和自然因素在控制污灌区重金属污染上存在相互作用而非相互独立，影响因子间可能通过复杂的联动机制增强对污灌区重金属浓度的影响。对其他重金属的影响因子分析时也发现了相同的规律，可见研究区重金属污染是多个影响因子共同作用的结果。

2.4 地理探测器模型与Pearson相关系数的比较

Pearson相关系数法是影响因子相关性分析的常用方法，因此采用Pearson相关系数法来验证地理探测器模型分析结果的可靠性，并以Pb影响因子分析为例进行对比。由图3可见，两种方法在确定最主要影响因子时结果一致，均认为DI是Pb污染的主要影响因子，其次为NDVI。然而，NDVI、DR、DP、PMT的Pearson相关系数为负，但地理探测器模型得到的解释力为正。这是因为地理探测器模型对影响因子相互关系的探测既包括线性关系也存在非线性关系，而Pearson相关系数探讨的是影响因子与重金属间的线性关系，分为正相关和负相关，如果影响因子与重金属之间出现非线性关系时，则会出现错误的判断。总体而言，Pearson相关系数法验证了地理探测器模型对污灌区中Pb污染主要影响因子的识别，但在较弱影响因子评估中存在一定差异。在对Zn、Pb、Cu、Cr、Ni、As的影响因子分析中，地理探测器模型和Pearson相关系数法的评估结果几乎相同，但对Co、Cd的影响因子分析中有所偏差。与Pearson相关系数法相比，地理探测器模型在检测重金属污染和潜在影响因子间的非线性关系中具有一定优势。

表3 各影响因子对8种重金属污染的解释力Table 3 The explanatory power of different influence factors for 8 heavy metals pollution

表4 影响因子间交互作用对Pb污染的解释力Table 4 The explanatory power of interactions between influence factors for Pb pollution

图3 地理探测器模型与Pearson相关系数法的比较Fig.3 Comparison of geographical detector model and Pearson correlation index method

3 结 论

(1) 研究区内土壤重金属呈现不同程度的富集，从富集系数平均值来看，Cu、Pb、Zn、As、Cd的富集程度较大，主要受到人为因素的影响；Co、Cr、Ni富集程度较低，主要受自然因素的影响较大，受人为因素影响较小。从富集系数空间分布来看，8种重金属富集系数分布情况相似，主要受到工业园区、煤矿区、农作种植等的影响。

(2) 地理探测器因子探测结果表明，7种影响因子中DI对8种重金属的平均解释力最大,其次为SOM和PMT。不同重金属的主要影响因子也存在不同，Cd、Pb主要归因于DI，Zn、Cu主要归因于SOM，As、Co、Ni主要归因于PMT，Cr分布受到PMT、DP、SOM、DI的共同作用影响。

(3) 基于地理探测器模型交互作用分析，不同影响因子之间的交互作用强弱不同但均为非线性增强，污灌区重金属污染是多种影响因子共同作用的结果。

(4) 与Pearson相关系数法比，地理探测器模型可以有效分析因子间的交互作用对重金属污染的影响，在检测潜在影响因子间的非线性关系中具有一定优势。