土壤是多种重金属元素的归宿, 土壤多金属复合污染已经成为自然界普遍的污染形势之一^[1].而金属冶炼活动、工厂排放、化肥农药投入和污水灌溉等一系列人为活动是造成土壤多金属污染的主要原因^[2].Hu等^[3]的研究发现, 中国珠江三角洲地区由于受到改革开放以来的工业和农业的迅速发展, 已经形成了一定量的土壤多金属污染.Quinton等^[4]的研究表明, 长期的地表径流、土壤流失也会引起土壤多金属的复合污染, 且不同的土壤类型对于重金属的复合污染程度不同.所以, 土壤多金属共同污染已经成为生态系统最为普遍的问题之一.

而土壤重金属污染评价作为确定地区性土壤污染程度和制定污染控制策略的重要参考而被人们广泛使用, 评价方法的选择对于重金属污染的评价结果影响较大^[5].目前, 国内外有关土壤重金属污染的评价方法已有很多^[6], 其中主要方法有: 单因子指数评价法、内梅罗综合污染指数法^[7]、地质累积指数法^[8]和潜在生态风险评价法^[9]等.而地质累积指数法作为判断区域性土壤重金属累积最简洁的方法, 常规是通过一个地区土壤点位元素实测含量对比土壤元素背景值来判断累积程度.而有些学者选取不同的土壤背景值以及土壤重金属含量本身就会带来不确定性(由于采样点位多少导致平均含量不同), 从而对土壤累积性评价造成一定差异, 并且污染程度分级标准跨度较大的缺点也会使评价结果不够准确^{[10, 11]}.

近年来一些学者结合三角模糊数来降低传统地质累积指数法对于评价结果的不确定性^{[12, 13]}, 在实际应用中结合α-截集技术、区间数隶属度计算方法, 建立土壤重金属污染评价模型, 能够规避由于不确定性引起的评价不准确的问题, 这种方法特别适用于由于精确度不高或数据丰富度不足的情况^{[14, 15]}, 并成功地应用于土壤重金属污染风险评价与水环境评价中, 使评价结果更加真实可靠.但通过综合分析改进后的地质累积指数评价方法仍存在不足之处, 主要表现在: 通过结合三角模糊数来增加模糊区间的定义, 最终得出的污染等级仍为一个确定的等级, 仅能体现研究区整体的土壤污染水平, 这就导致局部土壤污染水平可能被放大或缩小, 不能体现一个研究区局部污染情况, 特别是当研究区点位重金属含量差异较大时, 污染程度会被扩大化.

鉴于此, 本文以地质累积指数法为研究基础, 结合三角模糊数将研究区土壤重金属实测含量与区域土壤背景值模糊化后, 加入拉丁超立方抽样法(LHS)进行随机模拟平行化分析, 使评价结果更加合理化.为探究不同含量区间此方法的适用性, 选择对郑州市重金属含量差异较大的5种土壤类型(潮土、褐土、风沙土、黄绵土和新积土)进行独立评价, 再对比研究区整体评价结果, 通过找出不同情况下方法应用的差异性, 以期为土壤重金属累积性评价提供新思路.

研究选择中国河南省中部地区郑州市作为研究对象, 地理坐标为东经112°42′~114°14′, 北纬34°16′~34°58′, 是中国中部地区重要的农业中心城市, 总面积7 446 km².郑州市土壤类型表现出以西为褐土(图 1), 以东为潮土的特征, 并伴有部分新积土、黄绵土和风沙土, 5种土壤类型面积比约为: 12.5∶8.7∶2.2∶1.8∶1.

图 1

Fig. 1

图 1 研究区位置、土壤类型和采样点位示意 Fig. 1 Location of study area, soil type, and sampling point sites

本研究样点于2016年5月采集, 采用网格布点法对采样区进行均匀布点, 为尽量避免人为因素对土壤重金属含量的影响, 用5点采样法在样点周围100 m范围内等距离采集4个样品, 与采样点混合成一个复合样品.本研究共采集213个表层样(0~20 cm).所有样品在晾样间风干一周, 再去除石块, 植物残根和其他杂质, 研磨碾碎后过0.149 mm尼龙筛并混合均匀.采用硝酸-氢氟酸-高氯酸消解体系, 结合石墨消解仪自动消解.利用电感耦合等离子体质谱仪测定5种重金属(Cd、Hg、As、Pb和Cr)含量.

在综合考虑采样点的平均值和标准差的基础上^[16], 将污染物含量和土壤背景值分别表示为: A₁=(a₁, b₁, c₁), A₂=(a₂, b₂, c₂), 将a取值为数据最小值和数据均值减去2倍标准差之间的较大值, b表示为数据均值, c取值为数据最大值和数据均值加上2倍标准差后的较小值, 其中a≤b≤c^[17].隶属度函数μ_Ã(x)表示为^{[18, 19]}:

(1)

计算过程中利用α-截集技术能够将数据转化为区间形式的可信度水平, 再经过隶属度函数计算地区土壤重金属累积等级, 能够简化计算过程, 方法详见文献[20, 21]. α=0.9是容易被人们接受的置信水平, 故本文所有的置信水平α均取0.9.

地质累积指数算法表示为:

(2)

式中, C_n为土壤点位平均含量, B_n为土壤背景值.依据地质累积指数值(I_geo)把土壤中重金属污染程度分为7个等级(见表 1).

表 1 (Table 1)

表 1 地质累积指数区间及污染程度分级 Table 1 Range of geological accumulation index and classification of pollution degree 地质累积指数 污染级别 污染程度 Igeo＜0 0 清洁 0≤Igeo＜1 1 轻度污染 1≤Igeo＜2 2 偏中度污染 2≤Igeo＜3 3 中度污染 3≤Igeo＜4 4 偏重度污染 4≤Igeo＜5 5 重度污染 Igeo＞5 6 严重污染

表 1 地质累积指数区间及污染程度分级 Table 1 Range of geological accumulation index and classification of pollution degree

利用式(1)中的隶属度函数形成的曲线与x轴围成的面积0.5(c-a)与该隶属度函数的比值, 再通过逆变换法^[23], 可以得到三角模糊数随机模拟的公式:

(3)

随机模拟就是通过计算机程序随机产生[0, 1]区间上均匀分布的随机数u₁, u₂, …, u_m, (m为随机模拟次数), 再通过公式(3)可以得到一系列x的模拟结果x₁, x₂, …, x_m, 这样就可以把三角模糊数之间的运算转化为实数之间的运算.

根据表 1的地质累积指数区间分级, 经推导基于随机模拟三角模糊数地质累积指数各污染程度的可信度P可表示为:

(4)

根据公式(4)得出综合地质累积指数各污染程度的可信度.

表 2为研究区表层土壤重金属含量的基本描述性统计.由于研究区不同土壤类型的重金属含量具有一定差异, 所以对不同土壤类型进行分别评价, 可避免不同土壤特征的干扰, 还可以对不同重金属含量的情况进行适用性分析.总体来说, 风沙土重金属含量较小, 黄绵土重金属含量相对较大. 5种重金属有一个共同点, 那就是在风沙土中的平均含量都是最低的.对于Cd来说, 风沙土和新积土的平均含量在总体水平之下, 黄绵土平均含量是最高的, 是风沙土平均含量的2.14倍.其他3种类型土壤Cd平均含量从大到小为: 褐土 > 新积土 > 潮土, 其中潮土标准差最大, 表明在潮土中不同点位Cd含量差别最大.对于重金属Hg, 除了风沙土中平均ω(Hg)为0.024 mg·kg^-1, 其余4种类型的土壤中平均ω(Hg)差别不大, 其中新积土最高为0.057mg·kg^-1, 其他3种重金属含量从大到小排序为: 褐土 > 潮土 > 黄绵土.与Cd一样, 平均ω(As)最高的土壤类型为黄绵土, 为12.09mg·kg^-1, 这要比平均ω(As)最小的风沙土高了1.24倍.同样, 在黄绵土中ω(Pb)也是最大的, 为26.01mg·kg^-1, 其他4种土壤中从大到小排序为: 新积土(24.87mg·kg^-1)、褐土(24.52mg·kg^-1)、潮土(21.58mg·kg^-1)和风沙土(19.22mg·kg^-1).对于Cr来说, 5种土壤类型中含量从大到小的排序与Pb相同, 在黄绵土、新积土、褐土、潮土和风沙土的含量分别为71.16、70.47、68.91、66.86和56.23mg·kg^-1.

表 2 (Table 2)

表 2 不同土壤类型点位基础性统计数据 Table 2 Basic statistical data of different soil types 重金属 项目 潮土 风沙土 褐土 黄绵土 新积土 总体 Cd 最小值/mg·kg-1 0.059 0.069 0.070 0.124 0.099 0.059 最大值/mg·kg-1 0.826 0.185 0.604 0.394 0.308 0.826 平均值/mg·kg-1 0.161 0.103 0.172 0.220 0.163 0.167 标准差 0.098 0.036 0.074 0.060 0.048 0.081 Hg 最小值/mg·kg-1 0.014 0.013 0.019 0.023 0.034 0.013 最大值/mg·kg-1 0.232 0.050 0.308 0.160 0.123 0.308 平均值/mg·kg-1 0.051 0.024 0.053 0.047 0.057 0.049 标准差 0.038 0.010 0.034 0.026 0.023 0.033 As 最小值/mg·kg-1 3.76 3.75 5.39 9.80 7.26 3.75 最大值/mg·kg-1 13.65 10.05 20.55 25.75 12.55 25.75 平均值/mg·kg-1 8.28 5.39 10.65 12.09 10.54 9.60 标准差 2.41 1.72 2.46 3.01 1.32 3.00 Pb 最小值/mg·kg-1 17.1 16.5 18 18.9 19.6 16.5 最大值/mg·kg-1 36.2 23.2 78.3 53.8 33.1 78.3 平均值/mg·kg-1 21.58 19.22 24.52 26.01 24.87 23.36 标准差 3.20 1.39 7.97 6.58 3.30 6.12 Cr 最小值/mg·kg-1 34.9 30.50 24.00 58.00 45.80 24.00 最大值/mg·kg-1 87.6 76.50 116 91.8 96.2 116 平均值/mg·kg-1 66.86 56.23 68.91 71.16 70.47 67.55 标准差 11.58 11.85 12.10 7.29 9.51 11.86 pH 最小值 7.37 7.42 5.88 7.84 7.32 5.88 最大值 8.54 8.47 8.54 8.40 8.35 8.54 平均值 8.18 8.19 7.92 8.09 7.99 8.05 标准差 0.18 0.25 0.48 0.13 0.22 0.35

表 2 不同土壤类型点位基础性统计数据 Table 2 Basic statistical data of different soil types

利用ArcScene软件将研究区DEM与5种重金属(Cd、Hg、As、Pb和Cr)以及pH空间插值结合成立体的空间插值图(图 2), 能够更加直观地了解研究区重金属含量与高程之间的联系.研究区pH低值区出现在研究区西部地势较高的部分, 而随着向东地势的逐渐降低, pH随之逐渐升高, 但pH最高值与最低值相差不大, 研究区整体pH偏弱碱性.重金属Cd含量在西北部达到最高, 继而向西南部随着DEM的升高Cd含量逐渐降低到中等水平, 而向东南部随着DEM的降低Cd含量逐渐降低到较低水平.Hg的含量态势则与Cd略有不同, 最高值区出现在北部地势较低地区, 低值区依旧是在东部地区, 而中部形成了一条蜿蜒曲折的Hg含量较低地区, 被南、西、北三面Hg含量较高地区包围.As和Pb的含量则与研究区DEM趋势相同, 在地势高的地区两种重金属含量较高, 随着地势的降低, 两种重金属含量随之降低.对于Cr而言, 除研究区东南部有部分Cr含量低值区, 其余部分含量均相对较高.

图 2

Fig. 2

图 2 研究区重金属立体空间分布 Fig. 2 Stereo spatial distribution of heavy metals in the study area

利用公式(1), 根据土壤点位的基础性统计数据中的最大值、最小值、均值和标准差, 得出地质累积指数中三角模糊数的最小可能值a₁, 最可能值b₁以及最大可能值c₁, 采用α-截集技术将三角模糊数模糊化处理成图 3的区间形式.

图 3

Fig. 3

图 3 经α-截集技术处理后的研究区土壤重金属实测含量区间 Fig. 3 Measured concentration range of heavy metals in soil treated by α-cutting technique in the study area

模糊化处理后的相对区间越大, 表明土壤重金属污染程度的不确定性越大.从图 3可以看出, 经α-截集技术处理后, 相比其他4种土壤, 风沙土在5种重金属中的区间含量是最小的, 在Cd、Hg、As和Pb中的含量区间绝对宽度都是最小的, 表明风沙土在这4种重金属中的点位含量离散程度较低, 数据不确定性较小; 而黄绵土在Cr中的含量区间绝对宽度是最小的, 其余4种重金属在Cr中含量区间绝对宽度从大到小依次为: 褐土＞风沙土＞潮土＞新积土.潮土中的Cd和Hg含量区间绝对宽度都是相对最大的, 而褐土中As和Pb的相对含量区间绝对宽度是最大的.

目前, 学者在区域土壤重金属累积性研究中的地球化学背景值选择上, 大多数选择的都是该省的土壤背景值^[24], 而这对于小区域的土壤重金属累积性评价来说是不够准确的, 因为会存在省内地域上重金属含量的差异, 从而导致评价结果不准确.也有些学者将地质背景值、普通页岩的平均值和砂岩中的地球化学背景值当作研究区背景值参与评价^[25], 这都会造成评价结果的差异性.因此, 考虑到不同背景值的影响, 本研究中采用河南省土壤重金属背景值^[22], 将10%作为最大可能值, -10%作为最小可能值, 并通过α-截集技术将背景值处理成区间形式, 结果见表 3.

表 3 (Table 3)

表 3 河南省土壤背景值三角模糊化结果1) Table 3 Triangulation of soil background values in Henan province 项目 Cd Hg As Pb Cr 背景值 0.074 0.03 11.4 19.6 63.8 三角模糊化 (0.067, 0.074, 0.081) (0.027, 0.03, 0.033) (10.26, 11.4, 12.54) (17.64, 19.6, 21.56) (57.42, 63.8, 70.18) α-截集 [0.073, 0.075] [0.029 7, 0.030 3] [11.286, 11.514] [19.404, 19.796] [63.162, 64.438] 1)利用表 3中的河南省土壤背景值当做公式(2)中的Bn

表 3 河南省土壤背景值三角模糊化结果1) Table 3 Triangulation of soil background values in Henan province

根据研究区土壤重金属实测含量区间以及河南省土壤背景值经α-截集技术计算出的区间值, 得出研究区土壤重金属地质累积指数区间(图 4).如果区间值介于表 1中两个地质累积等级之间, 则表明对于污染程度的判断具有较大模糊性.根据图 4结果可知, 研究区土壤重金属地质累积指数中除了黄绵土中的Cd在两个污染等级之间, 其余的所有地质累积指数区间均没有跨越两个污染等级, 表明各元素污染等级的不确定性程度不大.

图 4

Fig. 4

图 4 研究区土壤重金属地质累积指数区间 Fig. 4 Index range of heavy metal geological accumulation in the soil in the study area

由于受自然过程以及人类活动的影响, 污染物受土壤或水环境中物理、化学和生物过程等一系列作用后, 污染程度具有一定的不确定性^[26].所以为减少不确定度, 将研究区重金属实测数据与地球化学背景值三角模糊化后的评价结果分别进行5 000~40 000次的蒙特卡罗抽样(MCS)和拉丁超立方抽样(LHS)下的随机模拟.因模拟结果过于繁杂, 所以只列出Cd的随机模拟结果(图 5).可以看出, 利用MCS模拟在30 000次随机试验后, 结果收敛, 而利用LHS模拟在20 000次的随机模拟后的结果收敛, 说明LHS模拟效率更高.为比较两种方法的误差, 故本文对两种方法分别进行30 000次下, 置信水平为90%和95%的随机模拟绝对误差分析(表 4).可以看出, 利用LHS随机模拟结果绝对误差都小于0.000 5, 而利用MCS随机模拟方法除了Cr的绝对误差较小以外, 其余4种金属的绝对误差都在0.002以上.所以, 在最终的模拟中选取30 000次下的LHS方法来评价5种金属的累积性.

图 5

Fig. 5

图 5 模糊化参数的Cd随机模拟 Fig. 5 Random simulation of cadmium with fuzzy parameters

表 4 (Table 4)

表 4 30 000次模拟下的MCS方法和LHS方法的绝对误差分析 Table 4 Absolute error analysis of MCS method and LHS method under 30 000 simulations 重金属 MCS LHS 绝对误差 90% 95% 90% 95% MCS LHS Cd 0.467 7 0.464 0 0.463 5 0.463 3 0.003 7 0.000 2 Hg -0.106 1 -0.100 1 -0.098 6 -0.098 9 0.006 0 0.000 3 As -0.940 1 -0.937 0 -0.937 9 -0.938 3 0.003 1 0.000 4 Pb -0.371 0 -0.368 6 -0.365 5 -0.365 3 0.002 4 0.000 2 Cr -0.539 8 -0.540 2 -0.540 0 -0.540 0 0.000 4 0

表 4 30 000次模拟下的MCS方法和LHS方法的绝对误差分析 Table 4 Absolute error analysis of MCS method and LHS method under 30 000 simulations

表 5为三角模糊数下的研究区5种土壤类型与总体的地质累积指数评价结果, 而图 6为结合LHS下的研究区5种金属地质累积等级可信度.由于研究区5种金属地质累积指数均小于3, 所以只列出小于3的部分.可以看出, 前者得到的是确定的地质累积等级, 而后者得到的结果为各地质累积等级的可能性.根据结果对比显示, 研究区Hg存在较大不确定性, Cd有较小不确定性, 而As、Pb和Cr的污染等级基本无争议.仅用三角模糊数方法的评价结果中, 潮土、褐土、黄绵土和风沙土Cd和Hg的地质累积指数均为1级轻微污染水平, 而结合LHS评价方法的结果中, 潮土、褐土、黄绵土和风沙土的Cd有超过65%的可能性为1级轻微污染; 潮土Hg的地质累积指数有64%的可能为0级清洁水平; 褐土中Hg有48%的可能为1级轻微污染水平, 而有47%的可能为0级清洁水平, 两者比例基本持平; 黄绵土Hg的地质累积指数为0级的可能性(56%)要大于为1级的可能性(44%); 新积土的Hg为0级的概率为40%, 占有一定比例; 而总体Hg只有36%的可能为1级, 超过55%的可能性为0级.所以说, 基于三角模糊数地质累积指数法评价出的研究区潮土、黄绵土以及总体的Hg评价结果与结合LHS法的评价结果具有一定差异, 这是由于研究区土壤点位Hg的含量差别较大, 而决定评价结果的是研究区土壤点位Hg含量的平均值, 这就导致含量低值区平均含量被扩大化, 使评价结果准确性降低.而研究区土壤点位Cr含量都较低, 评价结果基本都小于0, 这就导致方法的改进对于评价结果没有改变.

表 5 (Table 5)

表 5 基于三角模糊数的土壤重金属地质累积指数评价等级1) Table 5 Evaluation grade of soil heavy metal geological accumulation index based on triangular fuzzy number 重金属 潮土 风沙土 褐土 黄绵土 新积土 总体 Cd 1 0 1 1 1 1 Hg 1 0 1 1 1 1 As 0 0 0 0 0 0 Pb 0 0 0 0 0 0 Cr 0 0 0 0 0 0 1)数值表示等级, 标准参照表 1

表 5 基于三角模糊数的土壤重金属地质累积指数评价等级1) Table 5 Evaluation grade of soil heavy metal geological accumulation index based on triangular fuzzy number

图 6

Fig. 6

(a)潮土, (b)风沙土, (e)褐土, (d)黄绵土, (e)新积土, (f)总体 图 6 结合拉丁超立方抽样下的三角模糊数评价结果 Fig. 6 Combined with the triangular fuzzy number evaluation results under Latin hypercube sampling

(1) 郑州市5种土壤类型重金属含量差异较大, 风沙土中5种重金属点位平均含量都是最小的, 而黄绵土除了Hg之外, 其余4种重金属平均含量均为最大.除了Cr之外, 研究区东北部其余4种重金属含量较高, 而Cr则是西部和东南部含量较高, 土壤pH整体呈现由西向东逐渐升高的态势.

(2) 不论是基于三角模糊数的地质累积指数法还是结合LHS方法, 对于重金属含量差异较小或所有点位地质累积指数评价结果均小于0的研究区, 方法上的改进不能改变评价结果.也就是说, 当评价结果由于差异较小都趋于统一地质累积区间时, 或地质累积指数评价结果已经达到最小, 像本文的As、Pb和Cr, 虽然结合LHS方法能够评价出研究区的污染程度还有较小的可能性为轻微污染, 但在88%的可信度以上为清洁水平.不过结合LHS方法能够指出虽然研究区整体处于清洁水平, 但还是有很小一部分受到了污染, 再结合GIS方法能够用于人们研究更小区域的污染集中情况.

(3) 基于三角模糊数的地质累积指数法评价结果为确定的污染等级, 而这一确定的污染等级是由研究区总体的重金属平均值决定的, 所以在一定程度上低值区的点位平均值被高值区点位平均值放大了, 这就导致部分地区的污染水平被放大, 如本研究的风沙土Cd和Hg的污染等级均为0, 但研究区总体的污染水平为1.而这一观点反过来亦成立, 某些区域的污染水平可能被缩小.这就是常规和基于三角模糊数地质累积指数的局限性, 而结合LHS方法使评价结果变成各个累积等级的可能性后, 就缩小了这个缺点, 能够评价出研究区局部地区的土壤污染程度.