2016, Vol. 37
2.农业部环境保护科研监测所, 天津 300191
2.Agro-Environmental Protection Institute, Ministry of Agriculture, Tianjin 300191, China
近年来，随着产地环境调查工作的开展，我国基本农田重金属污染物空间分布情况已基本掌握[1~4]，而重金属污染来源问题成为防控其污染必须解决的关键问题之一. 重金属来源复杂多样，不但有母质贡献源，而且人为贡献也越发突出，如企业点源、 农业投入品、 大气降尘面源等，而且区域重金属污染还可通过社会经济发展指标差异性间接反映出来. 正是来源的复杂性及空间变异不确定性，造成了重金属源解析的困难性. 目前，科研人员对重金属来源解析技术进行的研究主要集中在对汇变异性的分析及对源特征的专家判断来建立源-汇之间的关系，采用的方法主要有主成分分析、 正定主成分分析、 地质统计学、 混合分布模型等，主成分、 正定主成分分析是一种数据降维分析方法，可以有效将有共源特征的重金属进行归类进而根据归类特征探寻污染源[5~8]，地质统计学方法是建立在空间相关理论基础上以变异函数为基本工具以揭示污染物空间分布特征而用于判定污染源的统计技术[9~11]，混合分布模型则是将数据分布进行分开拟合而确定人为源及自然源的一种有效工具[12, 13]. 近年来，随着随机模拟技术的发展，随机森林回归技术逐步被应用于土壤污染来源解析中来，随机森林回归克服了传统技术对数据的苛刻要求，不必担心过度拟合问题，是解决源污染分担计算的重要工具之一[14]. 因此发挥各自特长，将多种技术融合应用于产地重金属污染源解析中是污染源解析技术发展的必然选择，而相关研究还明显不足.
本研究选择在天津西北部武清区，武清是天津市主要集约化农区，同时也是北京排污河、 北运河及永定河等北京、 河北、 天津排污河流的主要通道，有50余年的污灌历史. 关于本区域土壤及水体有机污染物研究较多[15~18]，而很少关于本区域重金属污染状况及污染来源的研究. 本研究在武清区采用了578个土壤样品，采用主成分分析、 地质统计学及随机森林回归相结合的方法探寻农田重金属污染源并估计各源相对贡献，以期为本区域重金属污染治理及环境管理提供有效支持.
1 材料与方法 1.1 研究区概况，样品采集及分析天津市武清区属大陆型半干旱季风气候区,一般年份降水量在500~600 mm之间,7~8月降水占全年降水量的80%左右. 汛期洪水来水量年际变化很大，枯水季节河道一般无天然径流,只有少量污水. 北京排污河、 北运河及永定河的污水大部分在本区域内被用来灌溉. 武清区污灌自从1982年开始，到2000年以后，由于水资源越来越缺乏,污灌问题更加突出[19]. 研究区种植的农作物以小麦、 玉米为主，土壤类型主要由潮土、 壤质潮土、 沙质潮土、 黏化潮土、 盐化潮土、 湿潮土、 石灰岩类淋溶潮土等组成[20].
本研究在武清区29个乡镇共布设578个土壤表层(0~20 cm)采样点，见图 1所示. 每个样品代表面积约为100 hm2，每个样品由4~6个子样混合组成，采样时记录采样点经纬度坐标. 样品采集后自然风干捡去杂草、 石粒后过2 mm筛备用. 采用四分法取100 g，用玛瑙球磨机研磨至全部通过孔径0.149 mm尼龙筛，混匀后用于测定重金属含量. 准确称取0.5000 g准备好的土壤样品，采用硝酸-氢氟酸-高氯酸三酸混合液微波消解后采用ICP-MS(安捷伦7500) 检测重金属Cu、 Zn、 Pb、 Cr、 Cd、 Ni含量[21]，另准确称取0.5000 g土壤采用50%王水微波消解后采用原子荧光法检测As、 Hg含量，仪器为吉天AFS-9130[22, 23]. 采用复样、 盲样及穿插标样检测等方法控制检测精度及准确度.
|
图 1 研究区采样点位分布示意 Fig. 1 Distribution of sampling sites in the study region |
采用随机森林回归计算各源对土壤重金属的贡献分担，首先需获取研究区影响土壤重金属含量及空间变异性的基础数据，本研究收集的数据主要有以下5类. ①基础地理数据，主要包括可能引起重金属变异的污灌河流、 主要道路、 居民地位置等，并以此为依据计算各监测点到最近河流、 道路及居民地的距离； ②土壤基础数据，主要包括土壤类型，各乡(镇)耕地面积； ③投入品数据： 主要为各乡(镇)猪、 牛、 鸡养殖量，并以此为基础，结合各乡(镇)耕地面积计算单位面积耕地猪、 牛、 鸡等承压压强； ④污染源数据，主要为排放重金属的企业分布情况； ⑤社会经济发展数据，包括各乡镇面积、 人口数量、 农业人口数量、 各乡镇工业生产总值、 各乡镇粮食平均单产等，并以此为基础结合各乡镇面积，计算单位面积人口、 农业人口压力及工业生产压强. 基础地理数据依据全国1∶25万地形图获取. 土壤基础数据根据天津土种志[20]及武清区统计年鉴获得. 投入品，社会经济发展数据均来自于武清区统计年鉴[24]. 污染源数据原自于Google Earth软件查询结果.
1.3 高斯条件模拟值法高斯条件模拟值技术是地统计模拟技术中最常用的方法，高斯条件模拟技术首先要创建一个基于标准正态分布(平均值=0，方差=1) 绘制的随机分配值的格网. 然后，将由简单克里格法创建的半变异函数模型应用到栅格. 这样可以确保栅格值遵循输入数据集中的空间结构. 生成的栅格构成一个非条件实现，而且通过每次使用不同的包含正态分布值的栅格可生成更多的非条件实现[25, 26]. 如果已选择条件模拟，则会通过克里格法将非条件栅格条件化. 该过程利用每个位置处的克里格估计值(预测值)确保模拟值遵循输入数据值，且平均来看，克里格法预测得到复制. 高斯地统计模拟的一般工作流程包括： 一是数据准备如降聚，正态化转化，趋势移除等； 二是创建简单克里格模型，如计算并拟合实验半变异函数，确定领域搜索策略等，模型的选择、 拟合优度的评估及搜索领域的确定采用交叉验证法，一般采用标准化平均误差及均方根标准化误差两项指标来评估结果的无偏性及最优性，标准化平均误差应尽可能地为 0，而均方根标准化误差应尽可能地接近于 1； 三是创造实现，如运行简单克里格模型，建立高斯条件模拟基础面，对数据进行反变换处理，增加移除趋势等； 四是进行高斯条件模拟计算处理，如选择模拟次数，确定模拟边界，确定计算统计量等[27].
1.4 随机森林回归随机森林回归算法是由 Breiman提出的基于决策树分类器的融合算法，其基本思想是基于统计学理论，利用 bootstrap 重抽样方法从原始样本中抽取多个样本，对每个 bootstrap 样本构建决策树，然后将所有决策树预测平均值作为最终预测结果. 随机森林回归可以看成是由很多弱预测器(决策树)集成的强预测器. 具有： ① 较少的参数调整； ②不必担心过度拟合； ③适用于数据集中存在大量未知特征； ④能够估计哪个特征在分类中更重要以及⑤数据集中存在大量噪声时同样可以取得很好的预测性能等优点. 具体计算步骤参见文献[28].
1.5 数据分析采用的软件数据描述性统计、 主成分分析采用的软件为IBM SPSS V20，数据条件模拟值及其图件制作采用软件为ArcGIS Desktop 10.3,随机森林回归采用的软件为R V3.1.2.
2 结果与讨论 2.1 重金属含量描述性统计表 1 给出了武清区土壤8项重金属Cd、 Pb、 Cu、 Zn、 Cr、 Ni、 As、 Hg描述性统计结果，天津市背景值及武清区背景值参见文献[29, 30]. 海河流域平原区土壤重金属平均含量参见文献[2]. 可见，研究区土壤重金属平均含量与海河流域平原区基本相当，Cd、 Cr、 As略低，Pb、 Cu、 Zn、 Ni、 Hg略高，与天津市背景值相比，Cd、 Pb、 As、 Hg有所增加，而Cu、 Zn、 Cr、 Ni依然低于背景值，这与张文具等人研究结果一致[30]. 但与本区域背景值相比，除Cr减少外，其它7项重金属含量均有增加. 结合Cr的中值及峰度、 偏度统计值可知，Cr基本处于背景水平，而其它7项重金属均受外界环境的影响出现增加态势，其中Cd、 Pb、 Hg的峰度、 偏度值较高，偏态分布明显，且其平均值也有增加，尤其是Hg增加突出，可知Cd、 Pb、 Hg受到外界污染，且存在点源性污染高值.
|
表 1 研究区表层土壤重金属含量描述性统计结果 Table 1 Summary statistics of metal concentrations in topsoils in the study region |
2.2 主成分分析
主成分分析是一种对高维数据进行降维从而更好了解重金属之间关系的方法，本研究对8项重金属进行主成分分析，同时为了更好展示结果差异性，采用标准化最大方差法对主成分矩阵进行旋转，结果见表 2所示. 基于特征值大于1的限制，本研究共获得3个主成分，可解释数据68.27%的变异性. 主成分F1代表总变异的35.24%，主要有 Cu、 Zn、 Cr、 Ni、 As等5项重金属组成，描述性统计可知，研究区这5项重金属除Cr外虽有一定累积，但总体增加量不高，监测值基本服从正态分布，异常高值少，因此可初步判断，这5项重金属变异主要来源于土壤类型即背景的差异性. Pb、 Cd属于第二组分F2，代表总变异的19.18%，Pb、 Cd两种元素都在研究区内有一定的累积，而且从峰度、 偏度值可知，两者都有大值存在，说明受到点源的污染. Google Earth查询结果可知，在Pb、 Cd高值区共有23家重金属排污企业，主要企业类型为电镀、 电子、 化工类及畜禽养殖. Hg属于第三组分F3，代表总变异的13.85%，本区域Hg背景值为0.04mg·kg-1，而现在平均值达到了0.08mg·kg-1,中值也达到了0.07 mg·kg-1，平均值增加了2倍，可见，Hg存在较严重的外源污染情况，1982年天津市农科院中心实验室对北京排污河沿岸水体、 地下水及土壤监测结果表明，河水、 地下水及土壤都受到Hg的污染，水体污染尤其严重[31]，可见，Hg明显受到污灌的影响.
|
|
表 2 标准化最大方差法旋转后的重金属主成分矩阵承载 Table 2 Loadings of heavy metals on VARIMAX-rotated factors of different datasets |
2.3 重金属空间分布
数据经正态变换后采用高斯条件模拟值法绘制研究区8项重金属空间分布插值图，见图 2所示. 高斯条件模拟所依据的简单克里格拟合模型及模型参数见表 3所示. 块金值与基台值比值[C0/(C0+C)]一方面反映了样品分析误差及采样误差在微域内存在变异性的大小，另一方面则反映了人为因素对重金属空间变异性的影响[32]，块金比越高，反映人为因素对重金属空间变异性影响越强，中等强度的变异性(块金比在0.25~0.75之间)则反映了重金属空间变异受到人为及自然因素的双重作用，而弱变异则反映了重金属空间变异的本征特征[10]. 在8项重金属中，土壤As块金值比最小，表明As空间变异性受外界因子影响最小，其变异主要是受土壤类型变异影响，这与主成分分析结果一致，王馨慧等对北京排污河上游凉水河沉积物As含量的监测结果也表明，本区域As并非来源于污灌[33]. 而Hg、 Pb、 Cu、 Zn、 Cd、 Ni、 Cr块金比均属于中等强度变异性，说明还有外界因素影响其分布，而Hg受到的影响最强. 这一特点由重金属的变程大小也可以看出，As、 Ni、 Cr的变程均大于2 km，而受外界影响大的Hg、 Pb、 Cu、 Zn、 Cd变程均在1.7 km左右，由于8项重金属所采用的模型相同，变程的小大直接反映了空间自相关距离的大小，受外界人为影响大，其空间自相关自然下降，变程则相应减少.
|
图 2 武清区重金属(Zn、 Cu、 Ni、 Cr、 As、 Pb、 Cd、 Hg)含量分布 Fig. 2 Spatial distribution maps of heavy metals (Zn，Cu,Ni,Cr,As,Pb,Cd and Hg) in topsoils in the Wuqing district in Tianjin,China |
|
|
表 3 重金属空间变异半变异函数拟合模型及参数 Table 3 Best-fitted semivariogram models and parameters of the soil heavy metals |
图 2表明，8项重金属的空间分布有差异性也有相同性，高值区主要集中在武清区县城西南地区以及北京排污河一线，还有部分重金属在研究区东南部靠近大毕庄附近也有高值区. 具体而言，重金属Zn高值区主要沿北京排污河及永定河一线，距离河流越近，Zn含量越高，可见，Zn的变异必然是受到了污灌的影响. Cu的分布基本与Zn相同，主要差别在于Cu在东南部高值区与黏化潮土分布更加一致，而沿北京排污河的高值区更窄，说明Cu是受土壤类型及污灌的双重影响. Ni高值区分布与黏化潮土区基本一致，说明土壤类型是其主要影响因素. Cr、 As的变异由图 2可知与土壤类型有更直接的联系，黏质土壤Cr、 As含量明显高于其它类型土壤. Pb的污染点状比较突出，说明有排污企业的污染，如靠近北辰区大毕庄的高值区，在大孟庄附近的高值区. Hg的高值区主要沿北京排污河分布，说明主要是由于污灌引起的. Cd的分布与其它重金属有明显的差异，呈现四周高，中间低的趋势，显然这种变异性不是由土壤类型的变化造成的，但似乎也与污灌关系不密切.
2.4 重金属来源的随机森林回归分析研究区8项重金属描述性统计分析、 主成分分析及地质统计学分析逐步明析了重金属含量增加及空间变异性产生的原因. 但以上分析基本是对重金属汇展开的探讨，土壤重金属含量增加及空间变异的原因分析在很大程度上是基于对已有知识的掌握及专家判断，并没有在源-汇之间建立起真实的量化联系. 且在某些情况下专家判断存在较大的不确定性，更有时无法给出明显的答案，如对重金属Cd空间变异原因的探讨. 正是基于上以问题，本研究采用了随机森林回归技术构建重金属源-汇量化关系. 结果见图 3所示. 所构建的模型对重金属Cu、 Zn、 Pb、 Cd、 Hg、 As、 Ni、 Cr空间变异可解释性分别为94.43%、 95.64%、 90.85%、 91.41%、 91.54%、 90.11%、 96.39%、 96.11%，均达到极显著水平. 图 3可知，Cd的空间变异第一影响因素为粮食单产，也就是说乡(镇)粮食平均单产越高，土壤Cd的污染越重，高投入与高产出的关联性及Cd污染与土壤类型的不一致性说明农业投入品如矿物肥对产地Cd的长期累积性应引起重视,这与Chen 等[34]的研究结果一致，再结合描述性统计及主成分分析结果可知，Cd的变异性是以面源污染为主，点源污染为辅的. 土壤Pb变异主要贡献源为粮食单产、 鸡粪承载及距离居民地距离远近，可见，Pb的来源比较复杂，即有面源污染也有点源污染. Hg最主要变异源为距离河流及道路的远近，由于本研究主要道路及河流具有一定的重叠性，可见，Hg主要是来源于污灌. Cu、 Zn的变异源与上文分析的一致，Cu的变异更多来自于土壤类型的变化，而Zn的变异则主要由土壤类型及污水灌溉造成的. Ni、 As、 Cr的变异也主要来源于土壤类型的变化，虽然土壤类型的变化承担率有所不一.
|
图 3 不同贡献源对不同重金属变异贡献比例 Fig. 3 Sources importance scores for spatial variation of the heavy metals in the topsoils in the study region |
(1) 描述性统计结果表明，除Cr以外，天津武清区7项重金属平均含量均有所增加，其中Hg平均含量增加了2倍，Cd、 Pb增加次之，Hg、 Cd、 Pb均存在点源污染引起的高值点.
(2) 主成分分析结果表明，8项重金属可提取3个主成分，并基本判断Cu、 Zn、 As、 Ni、 Cr的变异来源主要是土壤类型的差异性，Pb、 Cd变异与点源污染有关，Hg主要来源于污灌.
(3) 地质统计学半变异函数拟合结果表明，除As外，其余7项重金属变异性明显受到人为活动的影响，Hg、 Pb影响最为强烈. 重金属空间分布表明，Cu、 Zn的变异与污灌有必然的联系，Ni、 Cr、 As的变异与土壤类型差异性有更直接的联系.
(4) 随机森林回归在进一步量化源-汇关系基础上，印证了以上判断，并进一步锁定土壤Cd的变异性是受到农业投入品及点源污染的双重影响.
| [1] | 环境保护部, 国土资源部. 全国土壤污染状况调查公报[EB/OL]. http://www.zhb.gov.cn/gkml/hbb/qt/201404/t20140417_270670.htm, 2014-04-17. |
| [2] | 中国地质调查局. 中华人民共和国多目标区域地球化学调查——海河流域平原区[M]. 北京: 中国地质出版社, 2014 : 16 -173. |
| [3] | Yang Z F, Yu T, Hou Q Y, et al. Geochemical evaluation of land quality in China and its applications[J]. Journal of Geochemical Exploration,2014,139 : 122–135 . |
| [4] | Teng Y G, Wu J, Lu S J, et al. Soil and soil environmental quality monitoring in China:a review[J]. Environment International,2014,69 : 177–199 . |
| [5] | Liu P, Zhao H J, Wang L L, et al. Analysis of heavy metal sources for vegetable soils from Shandong province, China[J]. Agricultural Sciences in China,2011,10 (1) : 109–119 . |
| [6] | Micó C, Recatalá L, Peris M, et al. Assessing heavy metal sources in agricultural soils of an European Mediterranean area by multivariate analysis[J]. Chemosphere,2006,65 (5) : 863–872 . |
| [7] | Esmaeili A, Moore F, Keshavarzi B, et al. A geochemical survey of heavy metals in agricultural and background soils of the Isfahan industrial zone, Iran[J]. CATENA,2014,121 : 88–98 . |
| [8] | Qu M K, Li W D, Zhang C R, et al. Source apportionment of heavy metals in soils using multivariate statistics and geostatistics[J]. Pedosphere,2013,23 (4) : 437–444 . |
| [9] | Wei C Y, Wang C, Yang L S. Characterizing spatial distribution and sources of heavy metals in the soils from mining-smelting activities in Shuikoushan, Hunan Province, China[J]. Journal of Environmental Sciences,2009,21 (9) : 1230–1236 . |
| [10] | Gu Y G, Li Q S, Fang J H, et al. Identification of heavy metal sources in the reclaimed farmland soils of the pearl river estuary in China using a multivariate geostatistical approach[J]. Ecotoxicology and Environmental Safety,2014,105 : 7–12 . |
| [11] | Lu A X, Wang J H, Qin X Y, et al. Multivariate and geostatistical analyses of the spatial distribution and origin of heavy metals in the agricultural soils in Shunyi, Beijing, China[J]. Science of the Total Environment,2012,425 : 66–74 . |
| [12] | Zhong B Q, Liang T, Wang L Q, et al. Applications of stochastic models and geostatistical analyses to study sources and spatial patterns of soil heavy metals in a metalliferous industrial district of China[J]. Science of the Total Environment,2014,490 : 422–434 . |
| [13] | Yang S Y, Chang W L. Use of finite mixture distribution theory to determine the criteria of cadmium concentrations in Taiwan farmland soils[J]. Soil Science,2005,170 (1) : 55–62 . |
| [14] | Hu Y N, Cheng H F. Application of stochastic models in identification and apportionment of heavy metal pollution sources in the surface soils of a large-scale region[J]. Environmental Science & Technology,2013,47 (8) : 3752–3760 . |
| [15] | 张枝焕, 陶澍, 沈伟然, 等. 天津地区主要河流沉积物中多环芳烃化合物的组成与分布特征[J]. 环境科学学报,2005,25 (11) : 1507–1516. |
| [16] | 张枝焕, 陶澍, 叶必雄, 等. 天津地区表层土壤和河流沉积物中甾烷与五环三萜烷的分布特征及污染源分析[J]. 生态环境,2005,14 (2) : 157–164. |
| [17] | 张翔宇, 张枝焕, 陶澍, 等. 天津地区表层土壤中菲系列化合物的分布特征和污染源分析[J]. 农业环境科学学报,2005,24 (1) : 74–78. |
| [18] | 刘瑞民, 王学军, 张巍. 天津地区表层土壤PAHs多尺度空间分布特征[J]. 农业环境科学学报,2008,27 (3) : 844–849. |
| [19] | 于卉, 郭勇, 刘德文, 等. 武清县污水灌溉情况及影响分析[J]. 水资源保护,2000 (4) : 3–6. |
| [20] | 蒋德勤. 天津土种志[M]. 天津: 天津科学技术出版社, 1990 : 69 . |
| [21] | Method 3052, Microwave assisted acid digestion of siliceous and organically based matrices[S]. |
| [22] | GB/T 22105.2-2008, 土壤质量总汞、总砷、总铅的测定原子荧光法第2部分:土壤中总砷的测定[S]. |
| [23] | GB/T 22105.1-2008, 土壤质量总汞、总砷、总铅的测定原子荧光法第1部分:土壤中总汞的测定[S]. |
| [24] | 武清区统计. 2006年武清统计年鉴[EB/OL]. http://www.tjwq.gov.cn/tongjj/tjzl/list.shtml, 2015-02-09. |
| [25] | Dietrich C R, Newsam G N. A fast and exact method for multidimensional Gaussian stochastic simulations:extension to realizations conditioned on direct and indirect measurements[J]. Water Resources Research,1996,32 (6) : 1643–1652 . |
| [26] | Dietrich C R, Newsam G N. A fast and exact method for multidimensional gaussian stochastic simulations[J]. Water Resources Research,1993,29 (8) : 2861–2869 . |
| [27] | Chiles J P, Delfiner P. Geostatistics:modeling spatial uncertainty[M]. New York: John Wiley & Sons, 1999 : 449 -471. |
| [28] | 崔东文, 金波. 基于随机森林回归算法的水生态文明综合评价[J]. 水利水电科技进展,2014,34 (5) : 56–60. |
| [29] | 张文具, 范华义. 天津市土壤中Cd、Hg、As、Pb、Cu、Zn、Ni环境标准制订[J]. 城市环境与城市生态,2002,15 (3) : 47–48. |
| [30] | 范华义, 张文具. 天津市土壤元素背景值的地域差异及成因分析[J]. 城市环境与城市生态,2002,15 (1) : 13–14. |
| [31] | 沈碧贞. 北京排污河污水灌溉对武清县的土壤、作物与地下水的影响[J]. 天津农业科学,1985 (1) : 33–37. |
| [32] | 戴明新, 师荣光, 赵玉杰, 等. 四川泸县农业土壤Cd含量空间变异性研究[J]. 农业环境科学学报,2007,26 (3) : 1093–1099. |
| [33] | 王馨慧, 单保庆, 唐文忠, 等. 北京市凉水河表层沉积物中砷含量及其赋存形态[J]. 环境科学,2016,37 (1) : 180–186. |
| [34] | Chen W P, Chang A C, Wu L S. Assessing long-term environmental risks of trace elements in phosphate fertilizers[J]. Ecotoxicology and Environmental Safety,2007,67 (1) : 48–58 . |