克里金插值法测土壤污染源

1、-1-A 题 城市表层土壤重金属污染分析摘要本文研究的是某城区土壤重金属污染分析问题，建立了重金属污染物传播的模型， 并依据污染物传播模型确定了污染源位置。为了对研究区内各个重金属元素的污染程度， 本文采用内梅罗单因子指数法对研究 区土壤重金属污染进行评价，其中 Hg、Cu、Zn 是研究区最主要的污染元素，且 Hg、 Cu 生态危害较大，主要分布在工业区；为了全面地反映土壤污染情况，本文采用内梅 罗综合指数法对研究区重金属污染进行评价，其中工业区达到轻度污染程度，交通区达 到了污染的警戒线；由于要充分考虑人为活动对重金属污染的影响，本文最后采用地累 积指数法对重金属污染进行再度评价，并使用 G

2、S+9.0 软件对地累积指数法的评价结果 进行 kriging 插值得到污染分布图，这些图直观的表明部分区域尤其是工业区Hg、Cu、Cr 达到中污染程度， Ni 、Pb、 Zn 污染程度在中度以下。为建立重金属污染传播模型，本文对研究区内各个重金属元素进行半方差分析，根 据半方差分析结果选择合适的模型进行克里金插值。 由于克里金插值法对极值位置的预 测并不准确，因此采取对采样点加入人为的扰动方法进行多次克里金插值，将统计结果 较为稳定的极值点作为预测的污染源。关键词：土壤重金属，内梅罗指数法，地累积指数法， kriging 插值-1-一、问题重述改革开放以来，随着经济的高速发展，环境污染问题也

3、随之而来，其中水污染和大 气污染是重中之重。由于土壤污染具有滞后性、隐蔽性和长期性，土壤污染所带来的食 物安全问题和生态安全问题往往不能引起人们足够的重视。随着近代工农业的发展，重 金属己成为当今世界倍受关注的一类公害。随着工业、城市污染的加剧和农用化学物质 种类、数量的增加，土壤重金属污染日益严重，污染程度加剧，面积逐年扩大。由于重 金属污染物在土壤中移动性差、滞留时间长、不能被微生物降解，并可经水、植物等介 质最终影响人类健康。因此，土壤重金属污染问题己经成为当今环境科学研究的重要内 容。现对某城市城区土壤地质进行调查。附件 1 列出了采样点的位置、海拔高度及其所属功能区等信息，附件 2列

4、出了 8种 主要重金属元素在采样点处的浓度，附件 3列出了 8种主要重金属元素的背景值。问题一，弄清 8 种主要重金属元素在城区的空间分布，并能对该城区内不同区域重 金属的污染程度进行分析对有效控制定污染治理措施有重要意义。问题二，若能通过分析数据，说明重金属污染的重要原因就能够有效防范重金属污 染。问题三，由于重金属的污染具有较大的危害所以能够分析重金属污染的传播特征并 由此建立模型确定污染源的位置就更加急迫了。问题四，若能收集更多的信息改进所建立的模型，就能够更好的研究城市地质的演 变模式，这对于更高效更准确地控制与防范重金属污染有很大帮助。二、问题的分析对于问题一，已知地区空间区域内部分

5、点（x,y,z）的重金属浓度，探求一个完整的、 连续的空间区域浓度模型。本文采用 MATLAB 绘制出空间 3 维地质图，并结合内梅罗 指数法、地累积指数评价法，对该城区内不同区域重金属的污染程度进行分析和评价。对于问题二，结合问题一的污染浓度结果，运用结合各个区域的功能作用，以及各 个工农业发展的实际情况，做出客观分析，得到重金属污染的主要原因。对于问题三，问题要求分析重金属污染物的传播特征，由此建立模型，确定污染源 的位置。首先我们做出假设，每种重金属都是可以在外界因素下可以自由传播的，不存 在重金属元素在一个地区堆积，只有传播能力的不同。由于海拔因素对其影响不大，故 忽略海拔因素对其传播

6、的影响。首先，利用 matlab作图，画出以坐标x、y及重金属浓 度大小的差值等高线图。这样，我们就可以看到在整个城市中某一重金属元素的大致分 布。其次，由此图及第一问计算出各点的污染程度可以看出重金属的传播特征。其传播 大致符合自然界规律，高浓度向低浓度传播，但由于其重金属元素的特殊性，其传播速 度较慢，且不同重金属的传播速度不同。由此就可以得到污染源的特征就是其污染元素 的含量与周围含量相比相对较大，即区域极值。最后，根据所有数据的特征，建立如下 模型来寻找其污染源。由上可知，将其寻找污染源的问题转化为寻找图像空间极值的问 题。空间极值点的检测：本文对研究区内各个重金属元素进行半方差分析，

7、根据半方差 分析结果选择合适的模型进行克里金插值。 由于克里金插值法对极值位置的预测并不准 确，因此采取对采样点多次加入人为的随机扰动后再进行克里金插值，将统计结果较为 稳定的极值点作为预测的污染源。对于问题四，为更好地研究城市地质环境的演变模式，应该结合当地的天然地貌， 查明其分布位置、形态分布特征、组合特征、过度关系与相对时代。还应该适当的结合 热工地貌：主要包括露天采矿场、人工边坡、水库与大坝等。结合气象与水文调查，调 查当地水文地质，以及植被情况调查等。-2-三、模型的假设1假设题中污染源均为点状污染源。2、假设重金属污染物传播大致符合自然界规律，高浓度向低浓度传播。3、假设污染源的特

8、征就是其污染元素的含量与周围含量相比相对较大，即区域极 值。4、假设各金属元素之间具有一定的相关性，相关性较强的金属元素可认为来源于 同一污染源。四、模型的建立与求解（一）问题一土壤重金属污染评价可以为土地的合理利用、环境管理、环境规划和土壤污染的综 合防治提供科学依据。目前存在多种评价方法，每一种方法都有它的适用范围，很多新 方法的提出是为了弥补前面某种方法的不足和缺陷，因而每种方法都有自己的优点，可 以很好地反映土壤环境质量的某个或多个方面。本文采用内梅罗指数法和地累积指数法 对土壤重金属污染进行评价通过matlab自带的V4插值法做出的重金属元素空间分布图形如下:-3-AsCdCrCuH

9、gD> QNi-4-5-#-PbZn图1重金属元素空间分布根据以上分布图形可以直观的看出重金属的分布情况，图上颜色由深到浅代表重金属元素含量分布的递减情况，由图可以明显看出重金属呈集中的点源分布。且工业区分 布相对密集，交通去和生活区次之，公园绿地区和山区分布较少。（1）内梅罗指数法a）单因子指数法（1）Si为土壤环境质量通过单因子指数法可以计算各样点的超标率，通过单因子评价，可以确定主要的重 金属污染物及其危害程度，一般以污染指数来表示，以重金属含量实测值和评价标准相 比出去量纲来技术污染指数；式中：Pi为i重金属元素的污染指数；Ci为重金属含量实测值; 标准值国家二级标准值（国家二级

10、标准）。但是上式中，S为土壤环境质量标准值，没有扣除自然背景值，因此不能区分外来因素的作用。为了突出外来因素的作用，对单因子指数用土壤环境背景值进行修正，修 正公式如下：各散点不同金属单因子污染指数见附录（）单因子指数法污染程度评判标准如下：表1单因子污染程度评定标准Pi值<11 vR 兰22cR 兰3>3污染等级0123污染程度无污染中等污染强污染极强污染b）综合指数法：单因子指数只能反映各个重金属元素的污染程度，不能全面的反映土壤的污染状 况，而综合污染指数兼顾了单因子污染指数平均值和最高值，可以突出污染较重的重金 属污染物的作用。综合污染指数计算方法如下：P综 S（p2 Pi

11、max2）/2（ 3）但是由于不同重金属对土壤环境、生态环境的影响不同，采用加权计算法来求平均 值比较合适，改进公式如下：对于权重 W的确立，Swaine按照重金属对环境的影响程 度，将环境研究中人们都比较关注的微量元素分成了三类，因一类、二类、三类微量元 素环境重要性逐渐下降，分别赋值为 3、2、1作为权重。本研究涉及的几种重金属其类别和权重分配如下表:表2综合指数法各金属元素权重值HgPbCdAsZnCuCrNi类别IIIInnnn权重33332222综合污染指数分级如下表:表3综合污染指数分级土壤综合污染等级土壤综合污染指数污染程度污染水平1P综w 0.7安全清洁20.7 v P综V 1

12、.0警戒线尚清洁31.0 v P综V 2.0轻污染污染物超过起始污染值，作物开始污染42.0 v P综w 3.0中污染土壤和作物污染明显5P综3.0重污染土壤和作物污染严重从以上计算公式可以看出，内梅罗综合指数过分突出污染指数最大的重金属污染物 对环境质量的影响和作用，在评价时可能会人为的夸大或缩小一些因子的影响作用，使 其对环境质量评价的灵敏性不够高，在某些情况下，它的计算结果难以区分土壤环境污 染程度的差别。（2）地累积指数法地累积指数是德国学者 Muller于1969年提出，广泛用于研究沉积物及其它物质中 重金属污染程度的定量指标，它不仅考虑了沉积成岩作用等自然地质过程造成的背景值 的影

13、响，同时充分注意了人为活动对重金属污染的影响。因此，该指数不仅反映了重金 属分布的自然变化特征，而且可以判别认为活动对环境的影响，是区分人为活动影响的 重要参数，今年来它被国内外的学者专家广泛应用于人为活动产生的重金属对土壤污染 的评价。其计算公式如下：ige。=log2【cn/（k Bn）（4）根据I值将污染等级分为6级，为了区别以国家二级标准作基准的污染评价，本文 采取累积污染的提法，对应污染程度为无污染至极强污染。表4地累积法污染分级Iv 00112233445> 5级数0123456累积污染程度无无中中中强强强极 强极强本文研究区各金属元素地累积法平均指数见下表：表5本研究去个金

14、属元素平均指数元糸HgPbCdAsZnCuCrNi评价结果-0.07762170.30568131-0.08180370.702897490.27371226-0.22634690.150746080.32683593积累污染程度无-中无-中无-中无-中无 -中无无 -中无-中由于表格数据不够直观，所以通过以上数据利用GS+9.0绘制污染程度平面分布图如下：AsCd«M41»»ltn2KMnBMI他也2M54CrCu9551191032&fr54D«5IIfliaJ2B«54NiHg«5511P1032SS54i!tiin2

15、軸&4PbZn图2各金属元素指数分布图由图可以直观的看出生活区主要存在 Cd、Cr、Hg、Cu、Ni、Zn重金属污染，其中 Cd和Cr及Zn尤为突出。工业区的重金属污染较为严重，Cu、Hg尤其严重，Cr、Zn、Pb也存在。交通区除Hg污染相当严重外，Zn、Cd、Cr、Ni的污染也令人担忧。山区 和公园绿地整体污染较轻但山区人存在 Zn污染，公园绿地存在 Cd和Zn污染。（二）问题二我们知道按重金属污染物产生的部门，土壤重金属污染源主要分为工业污染源、交 通污染源生活污染源农药和肥料等。工业污染按工业类型可分为冶金工业、化学工业、 造纸工业、制革工业、纺织工业、印染工业、动力工业等。重金

16、属污染随工业废气废水 及废物排放，经过大气和水体最终进入土壤，使土壤中这些重金属含量明显增加或引起 污染。除此之外，交通运输也是土壤重金属的重要来源。机动车尾气作为环境污染源， 是土壤化学污染的一个重要来源。随着我国运输业的发展，机动车拥有量迅速增加，致 使交通干线旁土致使交通干线旁土壤重金属污染日趋严重。车辆排放的污染物可以通过多种途径被传播或者沉降到环境中如大气扩散，干湿沉降等方式，通过路面径流排放到 路沟或者路旁农田。而城镇居民聚集处也产生大量生活废水、废物，导致重金属最终沉 积在土壤中，造成污染。水系作为污染传播的一个载体，其分布广、范围大，对区域内 重金属的污染影响不可忽视。我国工矿

17、生产、农业活动、污灌等人为活动都是造成土壤 重金属污染的重要来源。表6中国土壤重金属的主要来源来源重金属矿场开采、冶炼、加工排放的废气、废水和废渣Cr、Hg、As、Pb、Ni、Mo煤和石油燃烧过程中排放的飘尘Cr、Hg、As、Pb电镀工业废水Cr、Cd、Ni、Pb、Cu、Zn塑料、电池、电子工业排放的废水Hg、Cd、P b、Ni、Zn :Hg工业排放的废水Hg染料、化工制革工业排放的废水Cr、Cd汽车尾气Pb农药、化肥As、Cu、Cd为了更好地分析城区表层土壤中重金属污染物的来源和类别，本文通过相关分析研究区各种重金属污染物的来源。相关分析就是研究两个或两个以上变量之间相互关系的 统计分析方法

18、。在环境研究中常用来定性的分析研究各环境变量之间的相互关系。环境 研究中最常用的是pearson相关分析法，它所计算出相关系数取值于一之间。当相关系 数为正，表示两变量为正相关，取值越大相关性越强而当相关系数为负，表示两变量负 相关，值越小相关性越强；当pearson相关系数为0或接近0，表示这两个变量不相关或几乎不相关。为了分析研究区各种土壤重金属的来源，我们采用了相关分析对研究区 重金属含量数据进行了相关分析，我们采用了pearson相关分析对这8种重金属含量数据进行了相关性分析，结果见表六。通常情况下通过以下取值范围判断变量的相关强度： 相关系数 0.8-1.0 极强相关0.6-0.8

19、强相关0.4-0.6中等程度相关0.2-0.4 弱相关0.0-0.2极弱相关或无相关表7研究区各金属元素相关系数AsCdCrCuHgNiPbZnAs-0.25470.18900.15970.06440.31660.28990.2469Cd-0.35240.39670.26470.32940.66030.4312Cr-0.53160.10320.71580.38280.4243Cu0.41670.49460.52000.3873Hg-0.10290.29810.1958Ni-0.30680.4364Pb-0.4937Zn-由上表分析得，强相关：Cr、Ni ； Cd、Pb ;中等程度相关：Cr、C

20、u ; Cu、Hg ; Cu、Ni ; Cu、Pb ; Cd、Zn ; Cr、Zn ;Ni、Zn ； Pb、Zn ；弱相关：As、Cd ; Cd、Cr ; Cd、Cu ; Cd、Hg ; As、N ; Cd、Ni ; As、Pb ; Cr、Pb ; Hg、Pb ; Ni、Pb ; As、Zn ; Cu、Zn ;土壤中Cr、Ni; Cd、Pb;显著正相关，且相关性较强；其次为 Cu、Cr; Pb、Cu;影响因素分析结果表8各功能区内梅罗指数分布功能区AsCdCrCuHgNiPbZnP综合1r 0.10.34 10.140.42 10.120.160.140.930.68 120.140.560.

21、081.321.310.20.231.151.033r 0.020.050.030.05 j0.010.080.020.020.06 :40.080.490.10.560.890.140.120.960.745r 0.10.310.050.190.180.080.110.460.35结合相关程度分析和影响因素分析结果，得出在工业区排放的Cu、Hg、Zn的污染浓度较高，其主要来源为工业成产中的废物废水的排放。而在交通区，Hg、Zn的浓度较高，其主要来源于汽车的尾气排放。而在生活区，随着人为活动的增强或者减弱，土 壤中各种重元素浓度也随之变化， 其主要污染以Zn为主。而在山区及公园绿化区等地， 由

22、于自然环境的优势，重金属污染较少。（三）问题三（1）克里格理论基础：a）对于任何一种估计方法，都不能要求计算值和它的实际值完全一样，也就是说， 偏差不可避免。但可以要求一种估计方法满足这样的条件明实际值与估计值的偏差的平 均为，或者说估计误差的期望应该等于零，即无偏。E（Zv-Z；）=0（ 5）b）估计值与实际值之间的单个偏差尽可能的小，即误差平方的期望值（估计方差） 应最小。2 * *-e 二 Var（Zv -Zv）二 E（Zv -Zv）二 Min（6）最合理的估计方法应当提供一个无偏估计且估计方差为最小的估计量。最常用的方法是用样品的加权平均求估计值，即：nZv 八'aZa（ 7）

23、ad最常式中入是加权因子，是各样品在估计时的影响大小，估计方法的好坏就取决于 如何计算或选择加权因子。以最小的估计方差给出待估单元平均值的无偏线性估计 量的方法，即所谓克里格法，它是一种最佳的局部估计方法。与直接使用普通克里格插 值法相比，复杂半变异函数模型以及趋势分析模型的盲目使用容易导致错误的数值估 计。因此，本文的克里格插值以空间分布特征研究为基础，直接利用分析所得的半变异 函数的理论模型进行普通克里格插值。克里格法实质上是利用区域化变量的原始数据和半方差函数的结构特点，对未采样点的区域化变量的取值进行线性无偏最优估计的一种方法。具体的讲，它是根据待估样 点有限邻域内若千己测定的样点数据

24、，在认真考虑了样点的形状、大小和空间相互位置 关系，它们与待估样点相互空间位置关系，以及半方差函数提供的结构信息之后，对该 待估样点进行的一种线性无偏最优估计。它既考虑了采样表面的总趋势变化统计特征， 又考虑采样表面特征的相关变化和随机变化。克里格法将上述三项分别称为采样表面的结构项、相关项和随机噪音。克里格法和 以往各种内插法的不同还在于它最大限度地利用了空间采样所提供的各种信息。在估计未知样点数值时，它不仅考虑了落在该样点的数据，而且还考虑了邻近样点的数据不仅 考虑了待估样点与邻近已知样点的空间位置，而且还考虑了各邻近样点彼此之间的位置关系。此外克里格法还利用已有观测空间分布的结构特征，这

25、种估计比其它传统方法更 精确、更符合实际。但克里格内插法是根据无偏估计和方差最小两项要求来确定加权系 数，其中关键的函数是半方差函数，如果半方差函数和相关分析结果表明区域化变量的 空间相关性并不存在，则此方法不适用。下面详细介绍一下半方差函数。半方差函数半方差函数定义为：区域化变量z(x)和z(Xi +h)的增量平方的数学期望 即区域化变量增量的方差。半方差函数既是距离 h的函数，又是方向a的函数。其计算 公式如下：1 N(h) 2r(h)z(xJz(Xj h)( 9)N(h)式中：r(h)为半方差函数值，半方差函数曲线图是半方差函数 r(h)对距离h的坐标 图形。N(h)是被分隔的数据对的数

26、量，z°)和z(Xi h)分别是点在X和Xh处样本的测 量值，h是两分隔样点的距离。半方差值与对应距离共同组合半方差图。 图中Co是块金方差，反映了最小取样距离 内由于样点变异性和测量误差所引起的方差，C0 Ci是基台值，是半方差随距离增加到一定程度后出现的平稳值 C0是结构方差，是基台值和块金系数之间的差值,a是变程， 当h_a时，可认为样点完全独立，不受空间影响，当h<a时，则要考虑一定的空间相关关系。半方差图揭示了空间样本变异与各个样本分离的偏离距离的关系，根据这种关 系可以选择使得样本方差和样本数目最佳的样本之间的距离。因此研究其空间变异的半方差函数模型由于土壤采样点各

27、属性值取值受误差影响, 应是有基台值的函数模型。主要有以下模型Qh = 0QU4 人 3(10)(11)球状模型 r(h)二 C0 C(一一 ),0 : h 乞 aa 2aCo C,h aQ,h =0咼斯模型r(h)h2c2C0 C(1 _ea ),h00,h 二 0(12)指数模型r(h)二h(Co+C(2),hcO为了使理论模型能最充分地描述所研究的某一区域化变量的变化规律，在建立模型的过程中要对模型进行最优拟合。通常情况下，根据半方差函数分布图初步选择几种合 适模型进行拟合，通过比较平均误差、均方根预测误差、平均标准差、标准化误差和平 均标准化误差等参数和预测误差图来选择最优的模型。由于

28、特异值的存在会对变异函数具有显著的影响，因此数据处理前剔除这些特异值是十分有必要的。本文采用域法识别特异值，即样本平均值加减倍标准差，在此区间以 外的数据均定为特异值，将其剔除。而数据的非正态分布会使得半变异函数产生比例效 应，进而造成实验半变异函数产生畸变，抬高基台值和块金值，增大估计误差变异函数 点的波动大，甚至会掩盖其固有的结构，因此应该消除比例效应。这些数据经常用对数 变换后参加计算，从而消除比例效应。本研究中，由于重金属 Hg、Cu具有强相关性，重金属 Cd Zn、Pb具有较强的相 关性，因此可认为分别来自于两类污染源。由于重金属As的含量较低，不形成污染，因此不讨论。重金属 Cr、

29、Ni的含量虽然不低，但生物毒性很小，因此也不考虑。重金 属含量数据的统计结果见表9。表9半方差函数 儿糸块金值(Co)基台值(Co+C)变程(m)块金值/基台值Co/(C+Co)理论模型残差平方和RSSHg、Cu9.119639.11963156401球面0.02Cd、Zn、Pb0.141.94515200.071979球面0.1731L5 T口口口口口口口口口V&rii&gram本文采用GS+9.0软件进行数据统计分析，半变异函数的计算及其拟合。步长的选 择采用距离组方法，即用一定距离范围内所有样点距离的平均值作为步长分组，本文步 长划分以为间距。本文选取拟合度最好的模型类型最适合模型决定系数较大，残差平方 和较小，计算了研究区土壤中种重金属含量的实验变异函数，并根据实验变异函数特点 进行了理论变异函数的拟合，拟合计算采用