湘江子流域重点污染企业影响区土壤重金属镉污染源识别

2020-05-02于靖靖师华定王明浩刘孝阳

环境科学研究 2020年4期

于靖靖，师华定,2*，王明浩,3，刘孝阳,2

1.中国环境科学研究院土壤与固废环境研究所，北京 100012 2.生态环境部土壤与农业农村生态环境监管技术中心，北京 100012 3.清华大学环境学院，北京 100084

土壤作为植物生长和生物生产的媒介，是人类农业生产的基础，维持着地球生态系统的平衡.研究表明，工业化和城市化导致了土壤中重金属污染程度的加剧[1]，土壤中重金属含量的升高会对人类和动物造成威胁[2-4].重金属元素中Cd(镉)的生物毒性最强，其在环境中生物有效性高、迁移能力强，有毒有害，容易通过直接吸入灰尘或植物富集进入食物链[5-9].因此，土壤重金属污染问题尤其是土壤Cd污染日益受到国内外研究学者的广泛关注[10-14].根据原环境保护部和原国土资源部于2014年联合发布的《全国土壤污染状况调查公报》[15]，我国土壤环境状况总体不容乐观，部分地区土壤污染较重，土壤总点位超标率为16.1%，其中Cd的点位超标率为7.0%.有研究[16-18]表明，在工业和城市快速发展地区的人为影响因素中，汽车尾气排放、工业废物、大气灰尘和气溶胶的沉积等是重金属污染的重要来源.因此，研究土壤中重金属的空间分布和主要来源，对于识别污染热点区域和评估潜在污染源具有重要意义.

土壤重金属污染源识别需要借助多元统计分析方法，如主成分分析(FCA)法、因子分析(FA)法、聚类分析、化学质量平衡(CMB)法、富集因子分析法、相关分析法、计算机成图法、同位素示踪技术等[19-21].Reza等[22]利用多元统计分析和相关分析法识别印度东北部矿区附近土壤重金属污染源的结果表明，Cd和Pb来源于人为因素，特别是煤矿开采活动，而Ni和Cr来源于自然因素或人为因素；ZHANG等[23]运用主成分分析法、皮尔逊相关系数法研究黄河灌区农业生产中重金属同源性的结果表明，交通运输、农村居民点和水体是土壤中Cd、Hg、Pb、Cr、Cu、Zn、Ni污染的主要来源，其中交通运输是最重要的因素；Cicchella等[24]利用GeoDA和ArcGIS软件绘制各元素和关联因子的分布图，利用污染因子和污染程度对土壤污染进行分析，结果表明，土壤污染具有不同的人为来源，Cd、Hg、Pb、Cr、Cu、Zn等元素污染主要出现在城市和工业地区，并在交通水平较高的道路附近出现，而其他元素(Co、Ni、Se、Tl、V)的含量与背景值有较大的一致性；谢伟城等[25]采用铅同位素地球化学示踪法对湘江长潭株段河床沉积物重金属进行了分析，发现长沙、湘潭段沉积物重金属的人为源比例小于67%，株洲段沉积物重金属的人为源比例达70%.然而，现有研究大多只是对整个研究区域内土壤重金属的潜在污染源进行定性研究，缺乏在显著影响区域更精细地把不同污染源的贡献率进行定量化.因此，了解掌握污染分布、识别与定量污染来源是污染防控工作的前提.

双变量局部莫兰指数方法(bivariate local Moran′s I)通过指数的计算对空间内两种要素作比较，并通过聚类分布图直观展示两种要素间的空间相关性.周俊驰[26]采用局部莫兰指数对研究区耕地土壤重金属进行分析，识别出研究区内铅锌矿及钨矿主要产区——高龙镇是重金属Cd、As、Hg的高值聚集区，说明这3种重金属在耕地土壤中的富集与工矿开采排放关系密切.贺祥等[27]基于广义加性模型(GAM)分析南京市PM2.5浓度影响因素时发现，SO2、CO、NO2等影响因素的模型拟合度较优，模型解释度较高.这表明双变量局部莫兰指数方法与GAM结合的方法可以很好地应用于土壤重金属污染来源的解析工作.

该研究利用双变量局部莫兰指数方法，首先确定研究区土壤中w(Cd)与不同行业类型密度的热点区，再利用GAM在热点区内定量分析各行业类型对土壤w(Cd)的贡献率，以期为区域企业管理、重点污染源监管提供参考.

1 材料与方法

1.1 研究区概况

湖南省是享誉全国的有色金属之乡，大多数有色金属和稀有金属矿藏的开采和冶炼主要分布在湘江流域内，湘江水系土壤重金属污染区域主要集中在郴州、衡阳和株洲3个江段[28-30].该研究选取位于湖南省南部郴州市的湘江流域内头山、大河滩两个子流域为研究区，包括桂东县、汝城县、资兴市、永兴县、苏仙区、永兴县、北湖区，东界江西赣州市，南邻广东韶关市，西接湖南永州市，北连湖南衡阳市、株洲市，素称湖南省的“南大门”.研究区地理位置为112°37′E～114°15′E、25°16′N～26°21′N，地处南岭山脉与罗霄山脉交错处，长江水系与珠江水系分流的地带，境内以山丘为主，地貌复杂多样，地势整体上呈东南高、西北低的特征.郴州市属于湘江支流(耒水、春陵水等)的上游，矿产资源丰富，被誉为“世界有色金属博物馆”.《2018年湖南省统计年鉴》数据显示，郴州市拥有工业企业1 076个，其中重工业713个，排全省第三.《2018年国民经济和社会发展统计公报》显示，湖南省10种有色金属产量为55.6×104t，比2017年下降5.6%；铁矿石原矿量为56.0×104t，比2017年增长21%.湖南省仅县属以上21个矿山的废水排放量就高达 2 375×104m3/a，废水排放量较大，且郴州境内河流中大多数金属元素含量均超过GB 8978—1996《中华人民共和国污水综合排放标准》的地面水或工业废水排放标准，如w(Cd)超过排放标准1.2～9倍[31-33]，工矿企业废水排放是导致土壤重金属污染的主要原因.

1.2 数据来源

通过查阅《2017年中国环境统计年鉴》《2017年湖南省统计年鉴》以及实地信息采集等方式核实并收集研究区相关污染企业数据.2018年5—9月，在对研究区内自然环境与工业污染源、污染物排放、人口状况等社会环境调研的基础上，采用系统布点法，借助2 km×2 km的网格在研究区内均匀布点共计2 638个(见图1).在每个采样点用梅花形采样法对土壤样品进行采集，在1.07×104km区域内，用不锈钢铲取5～6个表层(0～20 cm)土壤的子样品，充分混合，四分法后剩余约2 kg封装带回.将采集的土壤样品自然风干后碾碎、过筛，采用HCl-HNO3-HF微波密闭消解技术进行土壤样品的消解，用电感耦合等离子体原子发射光谱法(ICP-AES)测定土壤中w(Cd).在采集土壤样品时，注意记录每个土壤样品的编号、采样地点、采样经纬度以及周边污染企业等信息.

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

1.3 研究方法

1.3.1莫兰指数

空间自相关是邻近要素在空间上存在相近或相似的规律，主要通过空间自相关指标来测度.空间自相关的测度指标有很多，通常以莫兰指数(Moran′s I)的应用最为广泛[34-35]，它分为全局莫兰指数和局部莫兰指数两种，全局莫兰指数可以指出区域属性值的分布是聚集、离散或者随机分布模式；局部莫兰指数高值表明有相似变量值的面积单元在空间上聚集，低值表明不相似变量的面积单元在空间上聚集.具体计算方法如下：

(1)

(2)

(3)

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

(4)

(5)

VarI=E(I2)-E2(I)

(6)

式中，E(I)为莫兰指数的期望值，E(I2)为莫兰指数方差的期望值.全局莫兰指数值介于-1～1之间，大于0为正相关，小于0为负相关，绝对值越大表明空间自相关性越强，空间聚集越明显；相反，则空间自相关性越弱，空间分布越随机.局部莫兰指数所得到的空间聚类包括高-高、低-低、高-低、低-高4种类型.在土壤重金属污染中，高-高空间聚类可视为区域热点.

1.3.2广义加性模型

广义加性模型(GAM)是广义线性模型的推广[36-37]，该模型不需预先设定模型的参数，可在原始数据驱动下自动选择合适的多项式，对响应变量和多个解释变量进行非线性拟合，适用于各种分布函数(如高斯分布、二项式分布、泊松分布、卡方分布等)[37].其表达式一般可表示为

(7)

式中：g(μ)为响应变量的连接函数；E[Y]为响应变量Y的数学期望；a为常数截距；f1，f2，…，fm为连接解释变量的平滑函数，一般采用样条平滑法、局部加权回归散点平滑法、薄板平滑样条法；xi为解释变量.分析结果中用edf、P、方差解释率参数来表征模型的拟合优度，其中edf为估计自由度，用来判断响应变量与各解释变量是否属于线性关系；P为分析结果的显著性水平，评估各解释变量对响应变量影响的相关性；方差解释率表征模型对因变量总体变化的解释率.

2 结果与讨论

2.1 土壤中w(Cd)的分布

分析样点数据可知，研究区土壤中w(Cd)为0.032～14.01 mg/kg，平均值为0.63 mg/kg，中位值为0.31 mg/kg，二者均高于湖南省土壤背景值(0.08 mg/kg)[38]，表明土壤中w(Cd)可能受到人为活动的影响.峰度可表征数据集的平坦度，研究区土壤中w(Cd)的峰度(273.43)较高，表明数据分布较集中.偏度可表征数据的正态性，研究区土壤中w(Cd)的偏度(12.78)大于0，表明数据是右偏分布，即存在较多极端高值.

对土壤中w(Cd)进行对数变换，并用Surfer 11制作w(Cd)的等值线图(见图2)，结果显示，土壤中w(Cd)高值区主要分布在该流域西部的苏仙区、北湖区、永兴县，低值区主要分布在桂东县和汝城县东部.

图2 研究区土壤中w(Cd)的等值面分布Fig.2 Contour map of soil Cd content in the research area

2.2 土壤中w(Cd)与企业分布的相关性

以研究区内478家污染企业分布与 2 638 个土壤重金属样点数据为基础，采用空间分析方法对污染企业密度与土壤重金属污染的网格数据进行空间相关分析，并基于不同区域的空间关系类型信息进行污染源空间分布特征的判断.获取基础数据后，在ArcGIS 10.5中进行企业空间分布的核密度分析，将核密度值和采样点数据转换成覆盖整个研究区的1 km×1 km格网数据，在GeoDa软件中进行双变量空间相关性分析，从而判断污染源与土壤w(Cd)在空间分布上的协同特征.

图3 研究区重点污染企业分布Fig.3 Distribution map of key polluting enterprises in the research area

按照GB/T 4754—2017《国民经济行业分类》将研究区内478家(见图3)污染企业分为4种重点行业，其中采选业(黑色金属矿采选业、有色金属矿采选业)269个，冶炼和压延加工业(黑色金属冶炼和压延加工业、有色金属冶炼和压延加工业)132个，化学品制造业(化学原料和化学制品制造业)34个，其他行业(废弃资源综合利用业、仓储业、生态保护和环境治理业等)43个.通过GeoDa软件分别计算全部企业及4种重点行业的全局莫兰指数(见表1).由表1可见，研究区全部企业及4种重点行业的全局莫兰指数均为正，其中全部企业的空间相关性在P<0.05 水平下显著，其他4种重点行业的空间相关性均在P<0.001 水平下显著.不同行业类别的空间相关性存在一定差异，但研究区4种重点行业企业分布的全局莫兰指数均较高，说明不同行业企业在区域内的分布集聚性较强.研究区全部企业的全局莫兰指数比单个行业的全局莫兰指数低，说明研究区内478家企业分布较离散，总体来看没有形成显著的集聚区；但区域内单个行业的企业呈现显著的集聚分布特征，尤其是冶炼与压延加工业，其集聚性最为明显，这主要受资源禀赋、地形条件、交通等的影响.

表1 研究区不同行业类型的全局莫兰指数Table 1 Global Moran Index for different industry types in the research area

为进一步分析重点污染行业对土壤中w(Cd)的影响，采用Bivariate Local Moran′s I方法进行局部莫兰指数分析，分别绘制土壤w(Cd)与行业总体及4种重点污染行业的空间相关性聚类分布结果(见图4、5)，共分为5类，分别是土壤w(Cd)高值与企业密度高值聚类(H-H)、土壤w(Cd)低值与企业密度低值聚类(L-L)、土壤w(Cd)低值与企业密度高值聚类(L-H)、土壤w(Cd)高值与企业密度低值聚类(H-L)和无聚类.分析结果表明，研究区内行业总体及4种重点污染行业均与土壤w(Cd)呈显著相关.

土壤w(Cd)与行业总体及4种重点行业之间均表现出较强的H-H聚类趋势，主要分布在研究区中西部，说明该区域重点污染企业密度较高，环境污染比较严重，可能与该区域相关污染企业的不规范生产与未经处理的工业废水排放有关，同时该区域位于河流走向的下游，土壤重金属Cd随着河流的迁移作用而在此处富集.宋杰[39]研究表明，湘江流域的长期工业活动，使上游河岸沉积物各深度剖面中Cd、Pb、As、Cr等元素的含量都低于下游剖面.聚类结果为L-L类型的区域主要集中在研究区东北部，说明该区域污染程度低，企业密度低，总体上产业发展从产生环境污染到出现环境问题，会滞后很长时间，存在较为显著的滞后性.

聚类结果为L-H类型同样只在采选业、冶炼和压延加工业中存在，在采选业中主要分布在研究区东南部，在冶炼和压延加工业中主要分布在研究区北部和东南部，且在H-H聚类结果的外围，该区域环境污染不明显，但随着产业聚集度的升高，有可能是以后污染的高速增长区，而且该区域位于湘江流域上游，即使存在较高密度的企业，但该企业违规排放的污染物并没有对土壤重金属Cd产生明显影响，或者污染物随着河流的迁移作用并未在该区域积累.聚类结果为H-L类型只在采选业、冶炼和压延加工业中存在，主要集中在研究区西北部和中部，该区域重金属污染严重，但企业密度较低，可能是由于农业生产、土壤母质等其他环境因素造成的污染，也有可能是河流上游污染企业违规排放的污染物在此处汇集[30]，如曾清如等[40]研究发现，湖南省东河与西河流域非污染土壤中w(Cd)、w(As)远高于湖南省同类土壤的背景值，属于高背景区，其成因还需进一步确认.

注：“无”表示土壤w(Cd)与企业密度之间无聚类；“H-H”表示土壤w(Cd)高值与企业密度高值聚类；“L-L”表示土壤w(Cd)低值与企业密度低值聚类；“L-H”表示土壤w(Cd)低值与企业密度高值聚类；“H-L”表示土壤w(Cd)高值与企业密度低值聚类.下同.图4 研究区土壤中w(Cd)与行业总体密度聚类分析结果Fig.4 Cluster analysis results of soil Cd content and total firm density in the research area

图5 研究区土壤中w(Cd)与各类行业密度聚类分析结果Fig.5 Cluster analysis results of soil Cd content and various industry densities in the research area

2.3 土壤重金属Cd污染源识别

以土壤w(Cd)作为响应变量，以不同行业的企业密度作为解释变量，建立GAM相应关系模型，采用R 3.5.3软件中mgcv软件包对数据进行分析.

根据局部莫兰指数相关性分析可知，聚类结果为H-H类型的区域是企业污染显著影响区.为进一步明确4种重点行业与土壤Cd污染的相关性，选取4种重点行业聚类结果都是H-H类型的区域，每种行业选定126个网格数据.将采选业行业密度、冶炼和压延加工业行业密度、化学品制造业行业密度、其他行业密度4个影响因素作为解释变量，w(Cd)作为响应变量，采用样条平滑函数对单个解释变量构建GAM，并分析每个解释变量对响应变量的显著性影响及方差解释率(见表2).结果表明，其他行业密度与土壤w(Cd)呈显著相关(P<0.01)，而其他3个变量与土壤w(Cd)均呈极显著相关(P<0.001)，表明各影响因素单独作为w(Cd)解释变量时均具有统计学意义，这与采用局部莫兰指数进行分析的结果相同；采选业行业密度、冶炼和压延加工业行业密度、化学制品制造业行业密度对w(Cd)影响的模型解释率(46.2%～51.2%)较高，表明这3个影响因素与w(Cd)构建模型的拟合度最优；其他行业密度对w(Cd)的解释率(19.3%)较低，表明其他行业密度单独对w(Cd)构建模型的拟合优度较差，模型对w(Cd)的解释能力也较差，但它通过了显著性检验，同样具有统计学意义.

表2 土壤中w(Cd)与各影响因素的GAM分析结果Table 2 GAM model analysis results of soil Cd content and influencing factors

注：*** 表示在0.001水平下变量是显著的；** 表示在0.01水平下变量是显著的.

edf=1表明解释变量与响应变量w(Cd)呈线性关系，edf>1表明解释变量与响应变量w(Cd)间呈非线性关系，其值越大，意味着非线性影响能力越强.该研究结果(见表2)表明，4种解释变量均与土壤w(Cd)具有非常显著的非线性关系，所以土壤w(Cd) 受这4种解释变量的复杂非线性综合影响.

GAM可以较好地分析受企业影响区域内不同行业污染企业对w(Cd)的影响程度.分析结果显示，在企业污染显著影响区内4种重点行业企业密度对土壤w(Cd)空间分布影响排序为采选业>化学品制造业>冶炼和压延加工业>其他行业，这种排序反映了不同行业类型的企业对土壤w(Cd)的影响程度，对于受企业影响区内的土壤Cd污染，采选业是其主要污染来源.众所周知，在金属采选过程中产生的废水、废渣、废气及大量堆积的尾矿库，重金属会在生物地球化学作用下被释放和迁移，从而对周边土壤、河流以及大气环境等造成污染[41-42].LIU等[43]对郴州市铅锌矿尾矿库溃坝区域的研究发现，在进行紧急处理后，部分土壤仍受到Cd、As、Pb、Cu等重金属的严重污染.魏本杰等[44]研究表明，湘江流域某冶炼厂周边土壤已受到重金属污染，其中Cd污染程度最高，Zn、Hg、Pb也存在一定污染.湖南省有色金属矿藏丰富，且主要分布在湘江中下游的长沙、株洲、湘潭、衡阳与郴州五城市及周边地区，随着金属采矿、冶炼活动影响的进一步加剧，湘江中下游己成为湖南省河流沿岸土壤重金属污染问题最突出的地区，这与笔者所得土壤Cd污染与采矿、冶炼活动关系密切的研究结果相类似.

该研究首先对整个研究区域内土壤Cd的潜在污染源进行定性研究，然后在污染热点区域不同污染源的贡献率进行定量分析.土壤中Cd污染成因复杂，陈秀端等[45]通过绝对主成分分数/多元线性回归与地统计相结合的方法定量解析土壤中多种重金属的主要来源为自然源、交通源、化石燃料燃烧源.秦治恒等[46]基于空间位置利用空间自相关方法，发现湘江流域土壤中w(Cd)分布与相关污染企业的位置存在很高的相关性.对比可知，Bivariate local Moran′s I与GAM结合的方法可以达到明确土壤中w(Cd)高值分布区域，并进一步定量解析相关污染重点企业对土壤中w(Cd)的影响，从而有针对性地提出区域企业管理、重点污染源监管及土壤环境治理策略.