基于地统计学的土壤重金属分布与污染风险评价

杨奇勇 谢运球 罗为群 谷佳慧，2 曾红春，2

(1.中国地质科学院岩溶地质研究所岩溶生态系统与石漠化治理重点实验室， 桂林 541004；2.中国地质大学(北京)地球科学与资源学院， 北京 100083)

引言

土壤是农业生产的基础，也是生态环境的重要组成部分[1]。它是环境中重金属迁移、转化的重要媒介，既是重金属聚集的汇，也是向其他系统迁移的源[2]。土壤对污染物有一定的自净能力，但是这种自净能力是有一定限度的。随着工农业发展、矿山开采等，人类向土壤排放的污染物越来越多，当进入土壤中的污染物超过土壤承受能力时，土壤就会向外界输出污染物，从而危害人类健康[3-4]。因此，研究土壤重金属元素含量的空间分布并对其污染风险进行定量评估有着重要的意义。

地统计学方法在环境污染物的空间分布和风险评估研究方面得到了越来越广泛的应用[5-7]。普通克里格法(Ordinary Kriging, OK)能够对未采样区域实现最佳无偏性估计，是空间预测中最常用的方法。但是，在污染评估与风险管理的实际应用中，人们感兴趣的可能并不是空间某一位置污染物含量的具体数值，而是其超出或低于某一阈值的风险概率分布。指示克里格法(Indicator Kriging，IK)能给出单一变量在一定阈值限制下的风险概率，多变量指示克里格法(Multiple variable indicator Kriging，MVIK)能够将给定阈值的多个变量综合成一个变量，并绘出综合风险概率分布图[8-11]，因而在环境风险评估中有较好的应用前景。如GOOVAERTS等[8]应用IK 法绘制了Swiss Jura地区土壤Cd、Cu和Pb等重金属元素的概率分布图，对该区域的土壤污染进行了风险评估; CHU等[9]利用台湾中部的1082个土壤样品，采用MVIK法对重金属Cr、Cu、Ni和Zn进行土壤污染综合风险评价，指出在给定的重金属阈值下，调查区域土壤具有高的污染风险概率，风险概率高达0.82；LIN等[12]应用IK法对台湾Chunghua县的Cd、Cu、Ni、Zn等土壤重金属元素进行单一重金属污染风险评价，并用MVIK法对土壤重金属污染进行了综合风险评价。

本文选取云南省广南县广南图幅(1∶50 000)的野外调查数据，利用OK法进行土壤As、Cd、Cr、Cu、Pb、Hg、Ni、Zn等重金属含量的空间变异研究，采用GB15618—1995《土壤环境质量标准》中的二级评价标准，利用IK法和MVIK法对各项土壤重金属引起的土壤污染风险和污染综合风险进行评估，以期为当地矿山开采、农业生产和环境保护等提供科学依据。

1 材料与方法

1.1 采样区域概况

采样区位于云南省广南县内，总面积约为402.91 km2，其中岩溶面积约为197.52 km2，占采样区域总面积的49.02%(图1)。采样区属于亚热带季风气候，年平均气温为16.7℃，年平均降水量为1 042.1 mm。采样区西南地势低，为广南盆地，其余地区为山地、丘陵。海拔处在971～1 701 m之间。采样区内土壤类型主要是酸性母岩红壤、酸性母岩黄壤、棕泥土和黄紫泥土。主要土地利用类型有林地、水田、旱地。林地主要分布在北部山地和东部的岩溶丘陵区，水田主要分布在西部，旱地主要分布在东部岩溶洼地。岩溶区内存在石漠化现象，且东部石漠化较严重。区域内有矿石开采与采煤等活动。

图1 研究区位置与采样点分布Fig.1 Location of study area and distribution of sample sites in study area

1.2 样品采集与处理

按照DZ/T 0258—2014《多目标地球化学调查规范(1∶250 000)》要求进行间距为1 km的网格划分，网格内预布点充分考虑土地利用的代表性，预布点布置在土壤易于汇集的平缓地、洼地、山间平坝地等部位，相邻样点距离大于500 m。采样时间为2016年,采集0～20 cm深度的土柱，共获取土壤样品410件(图1)。每个样点由1个主样点和2个分样点组成，分样点与主样点距离超过25 m。利用手持GPS重新记录主样点位置作为采样点坐标。

样品加工、处理严格执行DZ/T 0258—2014《多目标地球化学调查规范(1∶250 000)》。土壤样品风干后，过20目筛。按照4 km21个点进行样品组合分析，共得到组合样品102个，对组合样品进行土壤重金属含量测试，采用密码插入法，按照分析样品数插入一定数量的准确度和精密度一级标准控制物监控测试，确保样品分析的准确度和精密度。重金属元素和pH值的分析方法及检出限见表1。

1.3 数据处理与分析

相关分析、经典统计分析在SPSS 16.0中进行。重金属含量和重金属指示变换值的半方差函数模拟在GS+7.0中完成，专题图的制作在ArcGIS 9.2中完成。地统计分析的方法、原理详见文献[13]。

土壤污染风险指示克里格分析方法的主要步骤如下[14]：①根据问题的要求确定各项土壤重金属污染评价的阈值。土壤重金属污染阈值采用GB15618—1995《土壤环境质量标准》中的二级评价标准(表2)对采样区8种重金属进行污染评价。②确定土壤重金属污染指标的指示函数。根据土壤重金属污染阈值，利用指示函数对相应的采样点数据进行二态指示变换，得到各样点的指示变换值Zij(1或0)，用来评价相应采样点上的土壤重金属污染状况。③在GS+7.0 中利用土壤重金属指示变换值进行空间变异分析，拟合得到最佳变异函数模型。④将变异函数模型参数输入到ArcGIS 9.2 中进行OK插值分析，得到满足相应阈值的土壤重金属元素污染风险的空间概率分布图，计算平均概率。⑤将多个土壤重金属指标合并为一个综合指标

表1 样品分析方法Tab.1 Testing methods for soil samples

(1)

式中k——重金属指标个数

i——采样点编号

⑥利用综合指数重复步骤③、④，获取土壤重金属污染综合指数空间分布图，计算综合指数污染风险平均概率。

表2 土壤重金属污染阈值Tab.2 Threshold of pollution for soil heavy metal mg/kg

1.4 空间评价精度验证

由于组合样点数量有限，研究采用交叉验证的方法对空间评价精度进行验证[15]。具体的做法是：①从所有样点中抽取一个样点作为验证样点，利用剩余样点进行空间预测。②将抽取的样点依次替换剩余样品中的其他样点作为验证样点，每次替换后都用剩余的样点进行一次空间预测，直到所有的样点都替换完毕。③选取平均绝对误差(MME)和均方根误差(RRMSE)[16]对空间预测精度进行评价

(2)

(3)

式中O——验证数据P——预测数据

n——验证数据的数目

2 结果与分析

2.1 土壤重金属的描述性统计分析

进行调查区土壤重金属含量统计分析，发现岩溶地貌区土壤重金属含量与非岩溶地貌区土壤重金属含量存在很大的差异性(表3)。从表3中可以看出，岩溶区土壤重金属含量最大值(除Ni元素外)、平均值要显著高于非岩溶地貌区，岩溶区pH值平均值也高于非岩溶地貌区。这与已有的研究结果一致[17-18]。参考表2中土壤重金属污染阈值标准，从平均值来看非岩溶地貌区只有Cd、Cu元素存在污染，但是岩溶地貌区除了As、Cr、Pb元素外其他重金属元素都存在污染。从偏度系数来看，所有采样点的重金属元素含量都不呈正态分布，但经过对数转换后，各项重金属元素基本呈正态分布，适合利用OK法进行空间预测。

2.2 土壤重金属元素间的相关性

土壤重金属元素与pH值相关分析结果(表4)表明8种重金元素两两之间都存在显著的正相关性(p<0.01)，相关系数在0.457～0.885之间，其中Zn元素和Cd元素的相关系数最高，As元素和Cd元素的相关系数最低。Zn元素与Cd、Cr、Pb、Hg、Ni，Cr元素与Pb、Ni，Pb元素与Hg的相关系数都在0.8以上，说明它们在来源上具有很大相似性，As与Cd、Cu的相关系数在0.5以下，说明三者的来源具有复合关系[19-20]。pH值与各项重金属元素之间相关系数在0.365～0.673之间，相关性显著(p<0.01)。

2.3 土壤重金属元素含量的空间分布

表5为土壤各项重金属元素含量对数转换后的半方差函数理论模型和参数。从表5可以看出，采样区土壤重金属元素含量主要选用球状、高斯和指数3种模型进行半方差函数拟合，各模型决定系数都大于0.8，处在0.896～0.979之间，说明所选取的模型能较好地反映各项土壤重金属元素含量的空间结构特[21]。块金值与基台值之比C0/(C0+C)通常作为衡量变量空间相关程度的尺度，若其值小于0.25，则为空间强相关；处在0.25～0.75为空间中等强度相关；大于0.75为弱相关[13]。从表5可以看出，除Cu元素含量(块基比为0.322)呈现中等偏强的空间相关性外，其他重金属元素含量都具有强烈的空间相关性。说明采样区域土壤重金属含量的空间异质性主要是由岩性、成土母质等内在因素引起的，而人类活动等随机性因素影响较小。利用获得的模型及其参数(表5)，在ArcGIS 9.2中分别对8种重金属元素含量进行空间插值，得到它们的空间分布图(图2)。

表3 研究区域不同地貌类型土壤重金属含量统计特征Tab.3 Characteristics of soil heavy metals contents for different landforms in study area

注：下标1、2、3分别表示非岩溶区、岩溶区和全部样点。

表4 土壤重金属元素间的相关系数Tab.4 Correlation coefficients of soil heavy metals

注：所有的相关性在p<0.01水平显著(双尾)。

表5 土壤重金属含量半方差函数的理论模型和参数Tab.5 Theoretical models and parameters of soil heavymetals contents

从图2可以看出，8种重金属元素含量总体上都是东南部、北部岩溶区含量高，西南部非岩溶地貌区的盆地中重金属元素含量低；但在细节上，8种重金属含量在西北部、北部空间分布有显著差异。As、Pb和Hg含量空间分布相似，Cd和Cr含量空间分布相似，Cu、Ni和Zn含量空间分布相似等。重金属含量空间分布的相似程度与它们在统计上的相关性强弱是一致的。

2.4 土壤重金属污染风险评价

利用表2中土壤重金属污染阈值对各项土壤重金属元素进行指示变换，发现采样区的Pb元素含量的指示变换值都为0，说明采样区土壤没有重金属Pb元素污染；而Cr元素的指示变换值也只有1个为1，其余都为0，说明采样区基本可以忽略土壤重金属Cr元素的污染；As、Cd、Cu、Hg、Ni和Zn等重金属元素指示变换值为1的样点数量依次为33、84、17、19、27、27个，也就是说这6种重金属元素分别有32.35%、82.35%、16.67%、18.63%、26.47%和26.47%的样点土壤有污染。因此，采样区的土壤重金属污染评价只考虑这6种重金属元素。利用这6种重金属元素的指示变换值分别进行OK空间插值，得到各项重金属元素土壤污染概率空间分布图(图3)。

图2 重金属含量空间分布Fig.2 Spatial distribution maps for soil heavy metal contents

图3 土壤重金属污染概率空间分布Fig.3 Probability pollution distribution maps of each soil heavy metal

从图3可以看出，除Cu元素外，其他重金属元素高概率污染区域都集中在岩溶区，其中As和Hg污染的高概率区域基本与岩溶地貌区分布一致(图1和图3)。将Cu元素的概率分布图与Cu元素含量分布图进行比较，发现Cu元素含量高值分布的东南部并没有出现大面积的土壤污染高概率区域。这一现象与采样区域东南部的高pH值分布有关，与朱礼学[22]的研究具有相似的结论。Cd元素污染概率分布主要集中在0.7～1.0的高概率区间，而其余5种重金属污染区域主要集中在0.1～0.3的低概率区域。统计分析结果表明As、Cd、Cu、Hg、Ni和Zn重金属污染的平均概率依次为0.326、0.805、0.185、0.192、0.267、0.270。因此，在本研究采用的重金属污染阈值下，调查区域土壤具有较高的Cd和As重金属污染风险。

利用综合指数变换值Zij进行空间插值，得到土壤重金属污染综合风险概率分布图(图4)。

图4 土壤重金属综合污染概率空间分布图Fig.4 Probability pollution distribution map of soil heavy metals

采样区域重金属污染综合风险概率处在0.09～0.78之间(图4)，其中0.09～0.3之间的低概率污染区域面积约为193.36 km2，占整个采样区的47.99%，主要分布在采样区域的西南部盆地到东北部丘陵一带。重金属污染综合风险在0.4～0.78之间的高概率污染区域面积约为69.95 km2，占整个采样区域的17.36%，主要分布在采样区域的东南部岩溶区。采样区域重金属污染综合风险平均概率为0.335，在本研究采用的重金属污染阈值下，具有较高的土壤污染风险。大量的研究表明岩溶区土壤重金属具有较高的背景值[17-18]，但是由于岩溶区富钙偏碱的环境，却很少有出现重金属诱发的地方病报道[22]。因此，针对岩溶区高背景的重金属值和富钙偏碱的岩溶环境，研究确定岩溶地貌区土壤重金属污染阈值有着实际的意义。

2.5 预测精度评估

经过交叉验证法检验，各项重金属元素含量的OK法预测、土壤污染风险概率的IK法及土壤重金属污染综合风险概率的MVIK法预测都具有较低的MME值和RRMSE值(表6)，表明各项预测都有较高的精度，因而评价结果是可信的。

表6 预测结果精度验证Tab.6 Precision evaluation for spatial prediction

3 结论

(1)土壤重金属之间均呈现极其显著的相关关系，其中Zn元素与Cd、Cr、Pb、Ni，Cr元素与Pb、Ni，Pb元素与Hg的相关系数都在0.8以上，具有较高的同源性，而As与Cd、Cu的相关系数在0.5以下，说明三者在来源上具有复合性。

(2)除Cu元素具有中等偏强的空间结构特征外，其余各项重金属元素都具有强烈的空间结构特征，其空间异质性主要是由岩性、成土母质等内在因素引起的。重金属含量高值主要分布在东南部和北部的岩溶地貌区，土壤重金属含量低值主要分布在西南部的非岩溶地貌区。

(3)指示克里格分析表明，采样区域不存在Pb和Cr重金属污染，而As、Cd、Cu、Hg、Ni和Zn等重金属元素存在不同程度的污染风险，其中Cd污染风险最高，达到0.805。土壤重金属综合污染风险为0.335，高风险区域主要在东南部的岩溶地貌区，低风险区域主要分布在西南部的非岩溶地貌区。

1 郑袁明，陈同斌，陈煌，等．北京市不同土地利用方式下土壤铅的积累[J]. 地理学报，2005,60(5)：791-797.

ZHENG Y M, CHEN T B, CHEN H, et al. Lead accumulation in soils under different land use types in Beijing City[J]. Acta Geographica Sinica, 2005, 60(5): 791-797. (in Chinese)

2 周萍，文安邦，史忠林，等. 三峡库区不同土地利用土壤重金属分布特征与污染评价[J/OL]. 农业机械学报, 2017, 48(7)：207-213. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20170726&journal_id =jcsam. DOI:10.6041/j.issn.1000-1298.2017.07.026.

ZHOU P, WEN A B, SHI Z L, et al. Distribution characteristics and pollution evaluation of soil heavy metals of different land use types in Three Gorges Reservior Region[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(7)：207-213. (in Chinese)

3 廖海军．北京市密云水库上游土壤重金属污染调查评价[J]. 中国环保产业，2007, 2(8): 23-26.

LIAO H J. Pollution evaluation of soil heavy metals in the upper Miyun reservoir of Beijing City[J]. Environmental Protection Industry of Protection Industry of China, 2007, 2(8): 23-26. (in Chinese)

4 WU G, KANG H, ZHANG X, et al. A critical review on the bio-removal of hazardous heavy metals from contaminated soils: issues, progress, eco-environmental concerns and opportunities[J]. Journal of Hazardous Materials, 2010, 174(1-3):1-8.

5 司涵，张展羽，吕梦醒，等. 小流域土壤氮磷空间变异特征分析[J/OL]. 农业机械学报, 2014, 45(3):90-96. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140316&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2014.03.016.

SI H, ZHANG Z Y, LÜ M X, et al. Spatial variability of soil nitrogen and phosphorus in small water shed[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(3):90-96. (in Chinese)

6 GHADERMAZI J, SAYYAD G, MOHAMMADI J, et al. Spatial prediction of nitrate concentration in drinking water using pH as auxiliary Co-kriging variable[J]. Procedia Environmental Sciences, 2011(3):130-135.

7 JANG C S, LIU C W, LU K L, et al. Delimitation of arsenic-contaminated groundwater using risk-based indicator approaches around blackfoot disease hyperendemic areas of southern Taiwan[J]. Environmental Monitoring and Assessment, 2007, 134(1-3): 293-304.

8 GOOVAERTS P，WEBSTER R，DUBOIS J P. Assessing the risk of soil contamination in the Swiss Jura using indicator geostatistics[J]. Enviromental and Ecological Statistics，1997，4(1): 31-48.

9 CHU H J, LIN Y P, JANG C S, et al. Delineating the hazard zone of multiple soil pollutants by multivariate indicator Kriging and conditioned Latin hypercube sampling[J]. Geoderma, 2010, 158(3-4): 242-251.

10 JUANG K W, LEE D Y. Comparison of three non-parametric Kriging methods for delineating heavy-metal contaminated soils[J]. Journal of Environmental Quality, 2000, 29(1): 197-205.

11 BELKHIRI L, MOUNI L, NARANY T S, et al. Evaluation of potential health risk of heavy metals in groundwater using the integration of indicator Kriging and multivariate statistical methods[J]. Groundwater for Sustainable Development, 2017, 4:12-22.

12 LIN Y P, CHENG B Y, GUEYSHIN S, et al. Combining a finite mixture distribution model with indicator Kriging to delineate and map the spatial patterns of soil heavy metal pollution in Chunghua County, central Taiwan[J]. Environmental Pollution, 2010, 158(1):235-244.

13 王政权. 地统计学及在生态学中的应用[M]. 北京: 科学出版社，1999: 65-132.

14 杨奇勇，杨劲松，余世鹏. 禹城市耕地土壤盐分与有机质的指示克里格分析[J]. 生态学报，2011，31(8): 2196-2202.

YANG Q Y, YANG J S, YU S P. Evaluation on spatial distribution of soil salinity and soil organic matter by indicator Kriging in Yucheng City[J]. Acta Ecologica Sinica，2011，31(8): 2196-2202. (in Chinese)

15 MYERS J C. Geostatistical error management[M]. New York: Van Nostrand Reinhold, 1997.

16 WU C F, WU J P, LUO Y M, et al. Spatial estimation of soil total nitrogen using Cokriging with predicted soil organic matter content[J]. Soilence Society of America Journal, 2009, 73(73):1676-1681.

17 覃朝科，易鹞，刘静静，等. 广西某铅锌矿区废水汇集洼地土壤重金属污染调查与评价[J]. 中国岩溶，2013, 32(3):318-324.

QIN C K, YI Y, LIU J J, et al. Investigation and evaluation on heavy metal pollution in the waste water collecting depression in a lead-zinc mine of karst area[J]. Carsologica Sinica, 2013, 32(3):318-324. (in Chinese)

18 邓琴，吴迪，秦樊鑫，等. 岩溶铅锌矿区土壤重金属污染特征[J]. 中国岩溶, 2017, 36(2):248-254.

DENG Q, WU D, QIN F X, et al. Pollution characteristics of heavy metals in soil of lead-zinc mining in karst areas[J]. Carsologica Sinica, 2017, 36(2):248-254. (in Chinese)

19 CHEN T, YAN B. Fixation and partitioning of heavy metals in slag after incineration of sewage sludge[J]. Waste Management, 2012, 32(5):957-964.

20 QUAN S X, YAN B, LEI Q, et al. Distribution of heavy metal pollution in sediments from an acid leaching site of e-waste[J]. Science of the Total Environment, 2014, 499(1): 349-355.

21 DUFFERA M, WHITE J G, WEISZ R. Spatial variability of Southeastern U.S. Coastal Plain soil physical properties: implications for site-specific management[J]. Geoderma, 2007, 137(3-4):327-339.

22 朱礼学. 土壤pH值在CaCO3在多目标地球化学调查中的研究意义[J]. 物探化探计算技术，2001, 21(2): 140-143.

ZHU L X. A study to the multi-aim geochemistry agrology survey[J]. Geophysical and Geochemical Exploration，2001, 21(2): 140-143. (in Chinese)