城市表层土壤重金属污染分析

1、城市表层土壤重金属污染分析摘要本文以土壤重金属污染问题为背景，以数据为基础，研究某一城区的表层土壤重金属污染问题。建立了综合污染指数模型、AHP 法综合污染指数评价模型、灰色关联度模型和模型解决了题目中的四个问题。针对问题一: 首先建立 X、Y 和重金属浓度的三维坐标系，对给出的 319 组数据进行三次样条插值，得到八种主要重金属元素在城区内的空间分布图。然后利用综合污染指数和 AHP 法综合污染指数两种模型分别对该城区内五个区域的综合污染程度进行评估并进行对比验证。最终得到五个分区内八种元素的浓度分布和综合污染指数。针对问题二：通过建立灰色关联度模型，判断出八种重金属元素中哪几种之间关联度较

2、大，如 Ni 和 Cr 等存在明显的正相关，说明这些重金属很可能来自相同污染源。同时根据每种重金属元素在五个功能区浓度的分布差异和可能的污染来源，最终确定每种重金属污染的主要原因。针对问题三：首先对重金属污染物的途径进行分析，确定扩散方式。通过对条件的分析，最终利用连续点源扩散模型确定特征和浓度扩散特性，再由污染源周围浓度分布逆推出污染源具一定是高污染地区的浓度极值点的结论。置，并得到污染源位置不针对问题四：针对问题三模型的之处，增加风速或土壤含水率、ph 值等条件，提高模型与实际情况的贴近程度。最后选取增加风速一项信息，建立相应的 扩散模型，更好地研究城市地质环境的演变模式，尤其是污染物的扩

3、散模式。：型 重金属污染综合污染指数AHP 法综合污染指数灰色关联度模1问题重述1.1 问题背景随着城市经济的快速发展和城市的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。1.2 问题提出按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为 1 类区、2 类区、 5 类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境。为此，将所的城区划分为间距1 公里左右的网格子区域，按照每平方公里 1 个采样

4、点对表层土（0 仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照 2 公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤10 厘米深度）进行取样、景值。，并用 GPS采样点的位置。应用专门中元素的背附件 1 列出了采样点的位置、海拔高度及其所属功能区等信息，附件 2 列出了 8 种主要重金属元素在采样点处的浓度，附件 3 列出了 8 种主要重金属元素的背景值。现通过数学建模来完成以下任务：(1) 给出 8 种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2)(3)(4)通过数据分析，说明重金属污染的主要原因。分析重金属污

5、染物的特征，由此建立模型，确定污染源的位置。分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？一、问题分析2.1 问题（1）的分析第一问要求8种主要重金属元素在城区空间分布，在附件中给出319个采样散点坐标、浓度数据。鉴于三维图形可以直观、形象表现重金属元素在该城区的空间分布，应用插值法把一批离散样点转化为光滑曲面。常见插值方法很多，如多项式插值（lagrange 插值）、Hermite 插值、分段插位、Newton 插值等等。但是这些插值法的总体光滑度都不高，这对于地形图的应用是有缺陷的。考虑上述因作出浓度关于 X、素，选用整

6、体光滑度较高的三次样条插值法，通过Y 的三维图。第二问要分析该城区内不同区域重金属的污染程度。题目中已明确将城区分为五区，首先要求得各种元素在五区的分布图。然而，一种元素的污染指数无法反映该区综合污染程度，故应设法建立模型求各区综合污染指数。考虑到所给8种元素浓度差异较大，而浓度的大小与元素污染贡献相关，故应建立各元素污染2与综合污染的关系模型。最后分区求得综合污染指数，从而解决该城区不同区域重金属的污染程度问题。2.2 问题（2）的分析本问探求重金属污染的主要原因，即污染来源。在问题（1）的基础上，已可明确各元素分布规律和各区污染程度。知道，元素的化学性质对元素在土壤的相关性有一定影响，如

7、Cu、Cr 性质相似，关联度大，故在土壤中的迁移过程很相似，所以为探求来源可以建立模型研究土壤中重金属含量的相关性，若重金属含量有显著的相关性，说明可能来源相同。因此，建立相关性模型有助于分析重金属污染物的来源。2.3 问题（3）的分析首先通过对八种重金属元素海拔与浓度之间关系和表层土壤重金属元素的形成方式与存在机理的分析，确定了重金属污染物大气污染扩散模型。由于污染源是近似于连续均匀稳定的。在无 况下，污染物匀速在大气中向四周扩散，污染物的扩散服从扩散定律，即度成正比。在这些条件下，时间通过要研究法向面积的流量与它的浓度梯是连续点源污染扩散问题，由扩散定律、重金属元素物质质量守恒定律、气体泄

8、漏连续性定理出区域内的抛物线型偏微分方程，即列出连续点源扩散方程。与之相反，要确定污染源位置，可由污染源周围浓度分布逆推优化求解初污染源位置。具体数值实现方面，可以在区域极值周围进行优化求解，来确定污染来源。2.4 问题（4）的分析首先对第三问所建立的模型进行优缺点评估，然后为了更好地研究城市地质环境的演变模式，尤其是重金属污染物模式，应该让所建模型更加贴近实际情况并且增加信息量。比如收集城区尤其是污染源地区的风速分布，土壤含水率分布，地区 ph 值分布等信息。针对第三问的，在第四文中假设收集到的是城区内（污染源附近）的风速分布，建立相应的扩散模型解决问题。二、模型假设1、采样点测量准确，数据

9、真实反映污染程度；2、重金属污染物的扩散看作是空间某一连续点源向四周等强度地瞬时重金属污染物，重金属污染物在无穷空间扩散过程中不发生性质变化，且不计地形影响；3、重金属污染物扩散服从扩散定律，即的浓度梯度成正比；时间通过法向面积的流量与它4、重金属元素无法被降解，从源到采样点之间，污染物在途中保持质量守恒；5、受体污染物中某种元素是各污染源贡献的线性组合，各污染源之间互不相关；6、由各个污染源贡献的某元素的量有足够的差别，且采样点和分析期间变化不大；7、污染物浓度在 y、z 轴上的分布符合分布（正态分布）；8、在全部空间中风速是均匀的、稳定的。三、符号系统3、模型建立与求解5.1 问题一模型建

10、立与求解首先根据附件一中 X、Y 和各采样点重金属浓度数据，通过三次样条插值可得地形图(5.1)。4四Pi单项污染分指数P综综合污染指数Pi max所有元素污染指数最大值Piave第 i 样点的某种重金属的平均污染指数Ci样点此重金属的实测含量Si所在土壤适用的评价标准S一规则的球面面积Q1在(t,t t) 内通过 的流量Q2 内放射性物质的增量t任意扩散时刻C0污染源此重金属的实测含量（x0 , y0 , z0）污染源坐标（xi , yi , zi）样点坐标Q0污染源的重金属物质总量Q源强， kg / s x , y , z用浓度标准差表示的 x, y, z 轴上的扩散参数i强度的电源函数图

11、 5.1地形图根据附件中 X 与 Y 坐标和附件二对应的 8 种重金属元素浓度值，应MLAB，使用三次样条插值得到8 中主要重金属元素在城区的空间分布，如下图所示。图 5.2As 浓度分布图图 5.3Cd 浓度分布图5图 5.4Cr 浓度分布图图 5.5Cu 浓度分布图图 5.6Hg 浓度分布图6图 5.7Ni 浓度分布图图 5.8Pb 浓度分布图图 5.9Zn 浓度分布图7根据八种元素在城区的浓度分布图，需要建立模型求的城区的综合污染指数，这是因为单因子污染指数只能反映单个污染物的污染程度。当多种污染物在某一区域同时存在时，无法全面、综合地反映多种污染物的整体污染水平。因此，为了分析该城区不

12、同区域重金属的污染程度，需要计算各分区的综合污染指数P综 ，本文下述两种模型计算 P综 。5.1.1 模型一：综合污染指数模型首先由下式确定每种元素的单因子土壤污染指数2：Pi Ci / Si（1）若 Pi 1 则表示 i 样点已受到这种重金属的污染， Si 取为附件三中 i 样点背景值上限值。22 Pi m axPia v e2P根据下式计算 P ：（2）综综由上式得到每个测点的 Pi ，然后分别统计五个区的 Pi ，求各个区平均值，则到下表：5.1.1 重金属元素在各分区的综合污染柱状图综合污染指数模型兼顾单因子污染指数的平均值和最高值，能较全面地反映环境质量，而且可以突出污染较重的污染物

13、的作用，然而在计算过程中发现该模型只考虑了分指数的算数平均值和最高值，大大突出了污染最严重元素的贡献，未考虑其他重金属元素对综合污染的贡献。由此，引入模型二。5.1.2 模型二：综合污染指数评价模型为考虑各种元素对土壤污染的影响和作用，考虑采用权重系数来突出污染指8标在综合污染指数中的权重。计算公式如下：cijm1 n i1 p ppijpjp i jjimcj 1(3)i 0根据附件 3，每一种重金属元素其给出了背景值的平均值x 和标准偏差 ，规定重金属浓度在( x - 2 ，x 2 ) 范围内为类，浓度在( x 2 ，x 4 ) 范围内为类，浓度在( x 4 ， x 6 ) 范围内为类，浓

14、度在( x 6 ，+ ) 范围内为类。为评估不同浓度元素对土壤污染的作用，可以规定其系数如下：类为 0.2；类为 0.5；类为 1.0；类为 2.0。这些系数应该在现实中加以验证并做出改进。若重视某些点污染的影响可以加重其值，否则可以减小。由上述表述可知，每种元素污染指数 Pi 由平均值公式：5.2 八种重金属元素在五分区内浓度分布表由上表得到八种元素在五个区域分布的柱状图，如下:5.1.2As5.1.3Cd9分区As (g/g)Cd (ng/g)Cr (g/g)Cu (g/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)15.311404.90574.09779.827

15、146.626.79792.959384.16428.639655.54344.368243.0061250.45013.136155.031491.39131.85384.74521.20510.17122.9328.74819.43129.76643.732574.81950.444109.160856.4557.25786.145417.91454.874372.79819.77637.316195.2174.54670.416198.215求和24.4092092.809209.889479.4812471.67240.485423.9821521.4515.1.4Cr5.1.5Cu5

16、.1.6Hg5.1.7Ni5.1.8Pb5.1.9Zn针对城市表层土壤重金属污染的特点,参照土壤质量, 选择土壤的8 种采用主要重金属元素 As、Cd、Cr 、Cu、Hg、Ni、Pb、Zn 作为本次评价因子,AHP 法赋权值。依据层次分析法(AHP 法)确定权重的基本原理是构造判断矩阵，所谓判断矩阵就是征询们的意见，对诸因子的相互重要性给出评判，将判断结果用数值表示出来，统计汇总后写成矩阵形式(即判断矩阵)5。本文结合土壤污染的特征，通过相关资料的查找对每一个进行较客观的判断，经统计得到10之间的重要程度（即相对重要程度）的 88 A 判断矩阵。AHP判断用数量方式表达和处理, 实现了定性分析

17、和定量分析相结合,各个评判法把人的解决了评价模型中确定权重。5.1.10八种重金属元素权重表5.1.11A 判断矩阵表利用权重系数对五个区中的八种重金属元素进行综合评价，P综 表，如下表：到五个区的5.1.12五分区 P综 表利用 5.1.12 得下图;11背 景均 值As 3.6Cd 130Cr 31Cu 13.2Hg 35Ni 12.3Pb 31Zn 69综合11.4753.1152.3906.0484.1890.5532.9995.56859.68022.4005.0431.43118.41035.7271.0685.0017.122232.36630.5150.6520.6840.77

18、10.6550.7110.6270.4319.63641.0374.4221.6278.27024.4700.5902.7796.057149.67951.3542.8680.6382.8275.5780.3702.2712.87346.905AsCdCrCuHgNiPbZnAs12134112Cd1/211/21.521/21/21Cr12134112Cu1/32/31/314/31/31/32/3Hg1/41/21/43/411/41/41/2Ni12134112Pb12134112Zn1/211/21.521/21/21As (g/g)Cd (ng/g)Cr (g/g)Cu (g/g)H

19、g (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)均值5.67302.4453.5455.05300.1317.2561.74201.13背景均值3.61303113.23512.33169比值1.5742.3261.7274.1708.5751.4031.9922.915权重121341125.1.13重金属元素在各分区的平均综合污染柱状图由 5.1.1 和 5.1.13，两模型所求结果一致性很高，增强了结论的说服力。问题二模型建立与求解灰色关联度模型建立选取参考数列：x0 x0 (k) | k 1,2, n (假设有m 个比较数列xi xi (k) | k 1,2, n (则

20、称0 (n)其中k 表示时刻。i (n),i 1,2, mmin | x0 (t) xs (t) | max max | x0 (t) xs (t) |(k) stst(4)i| x (k) x (k) | maxmax | x (t) x (t) |0st0is为比较数列 xi 对参考数列 x0 在k 时刻的关联系数，其中 0,1 为分辨系数。in | x0 (t) xs (t) |tmaxmax | x0 (t) xs (t) |stms称式中大差。、分别为两级最小差及两级最一般来讲，分辨系数 越大，分辨率越大； 越小，分辨率越小。式定义的关联系数是描述比较数列与参考数列在某时刻关联程度的

21、一种指标，由于各个时刻都有一个关联度，因此信息显得过于分散，不便于比较，因此出定义：称为数列 xi 对参考数列 x0 的关联度给n1n i (k)ri k 1(5)由(5)易看出，关联度是把各个时刻的关联系数集中为一个平均值，亦即把12过于分散的信息集中处理。利用关联度4这个概念，可以对各种元素进行分析。在利用（4）和（5）计算关联度之前，需要对数据进行初始化处理，一般来讲，实际问题中的不同数列往往具有不同的量纲，而时，要求量纲要相同。因此，需首先对各种数据进行无量比较，要求所有数列有公共的交点。为了解决上述两个问题，行变换。在计算关联系数。另外，为了易于对给定数列进x )xn()x (, ,

22、 2( ), x(n) ，称x(x(给定数列数列。为原始数列的初始化5.2.2灰色关联度模型求解依照问题的要求，依次将每种重金属元素在319个点的浓度作为原始数列，并将其转变为相应的初始化数列，计算出八个数列两两之间的关联度如下表（这里 = 0.02）。5.2.1八种重金属元素的关联度5.2.2 重金属元素的相关性分析由上表易看出，与As元素相关性最大的是Ni元素，与Cd元素相关性最大的是Pb元素，与Cr元素相关性最大的是Ni元素，与Cu元素相关性最大的是Pb元素，与Hg元素相关性最大的是As 元素，与Ni元素相关性最大的是Cr元素，与Pb元素相关性最大的是Cd元素，与Zn元素相关性最大的是C

23、u元素。研究土壤中重金属全量的相关性可以推测重金属的来源是否相同，若重金属含量有显著的相关性，说明有相同来源的可能性较大，否则来源可能不同。对表层土壤各种重金属之间的相关性分析表明:在全部表层土壤测点中，Ni和Cr、As之间的相关系数分别为0.914和0.899，As和Cr、Hg之间的相关系数分别为0.868和0.826，均达到显著的相关水平，说明这些重金属来自相同污染源的概率较大。13元素As (g/g)Cd (ng/g)Cr (g/g)Cu (g/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)As (g/g)10.7140.8680.7240.8260.8990.7

24、510.688Cd (ng/g)10.7360.7810.6310.7270.8290.794Cr (g/g)10.7570.7400.9140.7840.725Cu (g/g)10.6360.7450.8030.801Hg (ng/g)10.6450.5010.452Ni (g/g)10.6240.568Pb (g/g)10.682Zn (g/g)1上表的结果还表明:Pb和Cu这两种元间存在着明显的正相关，因为这两种元属于亲硫元素，它们的构型和性质相近，在土壤中的地球化学行为十分相似。5.2.3分析城区重金属污染主要原因由重金属元素在各分区的平均综合污染柱状图，发现3区综合污染程度最低，5区

25、其次，2区污染程度最高，结果于情理之中。为更直观分析，作八种重金属元素在五分区的柱状图，如下：5.2.2八种元素在五分区柱状图由上图直出城区土壤重金属元素分布呈现一定的规律，表现为工业区、主干道路区污染很严重；生活区污染较严重；公园绿地区污染较轻，有六种金属元素超过背景值；山区受人为活动影响很少，污染最轻，各种元素均在背景值之下。整体而言，As,Ni元素的含量是最低的，其主要来自土壤本质；Hg、Cu、Cd的浓度值很大，通过上面八种单个元素在五分区的平均污染柱状图与综合污染柱状图对比分析，发现Cu、Hg和Cd三个柱状图与综合污染柱状图有很强的相似性，表明Cu、Hg和Cd对该城区土壤污染贡献最大，

26、是该城区重金属污染的首要元素。5.2.3八种重金属元素污染来源1主要污染物重金属元素（Hg）镉（Cd）铜（Cu）锌（Zn）铅（Pb）铬（Cr）镍（Ni）砷（As）制烧碱、 化物生产等工业废水和污泥，含、 蒸气冶炼、电镀、等工业废水、污泥和废气，肥料杂质冶炼、铜制品生产等废水、废渣和污泥，含铜冶炼、镀锌、纺织等工业废水和污泥、废渣，含锌，磷肥颜料、冶炼等工业废水，防爆燃烧排气，冶炼、电镀、制革、印染等工业废水和污泥冶炼、电镀、炼油、等工业废水和污泥、玻璃等工业废水、废气，硫酸、化肥、14经分析，推得：1区（生活区）含量最多污染元素为Hg、Zn、Cu，而且Cr元素在1区含量最高。 Cr元素区域明显

27、，主要来源于母质土壤。推测三种元素与生产生活废弃物堆放和生活区内车辆排出废气相关。2区（工业区）污染程度最重，重金属元素含量几乎均为最多。尤其以Hg元素为最高，该元素来源较单一，主要来自工业废气、废水的排放。如火力发电、钢铁冶炼、水泥制造等工厂都很可能是排放Hg的来源。3区（山区）重金属元素相差不大，均低于背景值。与其他区比较而言，Ni元素含量较高，此时Ni主要来自土壤本质，不为污染所致。4区（主干道路区）Hg的含量较高，一部分来源于工业区的污染，还与燃煤燃油相关；Pb主要来源于含铅的燃烧；Zn主要来源于汽车的轮胎磨损。5 区（公园绿地区）综合污染程度偏小，主要污染元素为 Hg、Zn、Cd 元

28、素，施用肥料如磷肥、含锌肥料等会带来 Cd、Zn 的污染，故推测这些元素主要来源于含这三种元素的化肥的使用。经整理、分析后得 5.2.4 表：5.2.4各分区污染来源分析表5.3 问题三模型建立与求解5.3.1模型思路引入由于重金属污染源为点源，而且无法向上扩散，故重金属污染物的可以看作模型。由于扩散相当于稀释，所以理论上哪里污染物浓度大就是污染源。但问题比较复杂，比如一个地方本身不是污染源，但周围到处是污染源，那周围的污染物都会“扩散”到这里并集中起来，这种地区土壤污染程度也很高，所以模型中，本地用地类型和周边用地类型都要考虑；地势高低也要考虑。初步估计此城市城区表层土壤重金属污染物可能的气

29、沉降、土壤迁移。途径大致有三种：地表水、大15AsCdCrCuHgNiPbZn生活区母质土壤堆肥、污水母质土壤含 Cu 化肥、受 2、4 影响母质土壤受 2、4 影响受 2、 4、5 影响工业区母质土壤冶炼、电镀等母质土壤冶炼厂、废水废渣工业三废，火力发电等母质土壤颜料、冶炼等废水废渣冶炼、镀锌山区母质土壤受周围影响 不大母质土壤受周围影响不大受周围影响不大母质土壤受周围影响不大受 周围 影响 不大主道路区母质土壤受 2、5 区影响母质土壤受 2 区影响汽车等燃烧母质土壤含铅 燃烧、交通运输汽 车轮 胎磨损公园绿地区含 As化肥绿化过程堆 肥、污染灌溉母质土壤含 Cu 化肥、含母质土壤受 2、

30、4 影响含 Zn、 磷肥地表水污染物的特点是污染浓度较大的地区一般属于附近区域内的地势较低点，而且重金属污染物属于守恒污染物质，大多数具有生物可给性，在水环境中不易为微生物所降解，而且能被水生生物积累富集。因此确定是否有重金属元素方式为地表水首先应该确定重金属元素浓度分布与高程的关系。以下给出每种重金属元素浓度最高的十个采样点在地形图上的位置：5.3.1每种重金属元素浓度最高的十个采样点在地形图上的位置由上图可知，每种重金属浓度较大的位置虽然地势一般较低但并不是周围一定区域内地势最低的地方，而且浓度较大的位置相对集中，不符合污染物由地表水的特点。同时，污染物在土体及和非饱和条件。污染物在土体及

31、水中迁移转化相当复杂，其水分状况可分为饱与水中迁移转化系统可分为系统内污染物迁移过程、系统内各相间交换过程(污染物固相液相浓度之间、污染物固相浓度之间、污染物液相浓度之间)、系统内污染物转化过程、系统内与系统外部污染物交换以及污染物在系统内发生物理、化学及生物作用的结果。而土壤迁移中横向影响的计算主要是考虑含水层内向扩散比，它是控制污染溶质扩散形态的重要。随着扩散比的减小，二维扩散流速形态亦随之成为细长形状如下图所示。因为所取土壤为每平方公里一个采样点的表层土（010 厘米深度），因此土壤迁移可以忽略不计。5.3.2 随扩散比的减小，二维扩散流速形态变化图16最后将重金属污染物的途径定为大气。

32、大气扩散的基本问题，是研究湍流传输与污染物浓度衰弱的关系问题。为了推断各种条件下污染物传输及浓度时空分布规律，必须将各种污染源、气象参数以及下垫面条件下的大气污染过程模式化。研究模式中的参数，以模式计算的形式给出空气污染浓度的时空分布规律。大气扩散模式很多。这里采取适用于开阔平坦地形条件下，续点源在正态分布假设下的扩散模式模型。5.3.2此模型建立模型首先是不考虑风速的，以重金属元素污染源正下方的地面为坐标0 0,原点，平均风向为 X 轴，铅直方向为 Z 轴，水平垂直于风向轴（X 轴）为Y 向，建立空间坐标系，则污染源距有效地面的高度为 H ，则金属元素污染源0,0）。位置坐标为时刻记作 t=

33、0，时刻 t 无穷空间中任意一点坐标将重金属元素从污染源为,zyz,t的浓度记为，若将重金属元素污染源的扩散看作是空间某污染物质，放射性物质在无穷空间扩散过程一连续点源向四周等强度地瞬时中不发生性质变化，且不计地形影响。梯度成正比，则：q i g adC时间通过法向面积的流量与浓度（1）i (i x, y, z) 是扩散系数， grad 表示梯度，负号表示由浓度高向浓度低的地方扩散。空间域 ，其体积为V ，包围 的曲面为S ，S 为一规则的球面, Sn (- x , - y ,1)zz(t,t t)内通过 的流量为：外法线向量为。则在t t nd dt（2）s 内重金属元素的增量为：Q2 C(

34、 x,Vy, z,t t)C(x,y,z,t) d V（3）从污染源的重金属污染物的总量为：t tQ p0dVdt（4）0t根据“质量守恒定律”和“气体泄漏连续性原理”，时间内通过所选曲面 S 的向外扩散的重金属污染物的量与S 曲面内重金属污染物增量之和，等于污染源在时间内向外的重金属污染物。则：171 Q2（5）即：t tVC(x, y, z, t t) C(x, y, z, t)dV q nd dt t t tp dVdt（6）0ts又根据曲面积分的 Gauss 公式： q nd divqdV(其中div 是散度记号)（7）sV t) C(x, y, z, t)C(x, y, z, t t

35、 tt tV tdV divqdVdt p dVdttt0tVt t lim C(x, y, z,t t) C(x, y, z,t) lim tkdiv(gradC)dt Ctttt0t0由以上两式得： C dV t divqdV t p dV t0tVV即为：C t dV divqdV p0dV（8）VVV由以上公式并利用积分中值定理得：C 2C 2C z2 ), t0,x, y, z 2Cidiv(gradC) i ( x2（9）tyy2z这是区域的抛物线型偏微分方程，根据假设扩散过程中浓度在 y 、 z 轴上的变化分布是分布，初始条件为作用在坐标原点的电源函数，记作C(x, y, z,

36、0) Q0 (x, y, z)（10）方程（9）满足方程（10）的解为：222x y e 4itzQC( x, y, z, t) 0, i( x, y, z)（11）( 4 t3) / 2i设检测时，t 已趋于无穷大，即污染已相对稳定，取 t=1000000s。且设污染源的坐 x y zi 标为（xy z0）,浓度为C0 ,假定此时的 0.0001 。对于任意一测量点，均有, 浓度为Ci ,则22004it2z0Q0exyz 222 z 2( x 2 )( y ) ( z z 2 )( t2Ci 00 0 i 0 i 0 i 0 4it4it4it e eOiz 2Ci Qi4it0 e( t

37、2i18则污染源的浓度可以表示为：( xi x0 )2 ( yi y0 )2 ( zi z0 )24itC C e（12）0i那么，对于任意假设的污染源，分别求出由在大区域内的各测量点而得到的浓度， C0 C01C02 C0m , m 1,2,其中 m 为大区域内测量点的个数。假设该区域内的极值点即为污染源，此时C0极 C01C02 C0m , m 1,2, ，对于污染源，真实的C0 将大于等于C极 ，且对于真实污染源的C0真里各数据的标准差将小于等于极值点做污染源的 C0极 里各数据的标准差。故在极值点周围做100m100m 的网格，对于有限范围的网格节点进行穷举求解，若满足ma x(C j

38、) C ，且jC ，其中 j 为网格点序号，j=1,2,，则第 j 个网格（C ）O极0极点将比极值点更接近污染源。这样依次进行比较，最后确定污染源的坐标为，浓度为C0 。此时，将污染源作为已知点再次进行三次样条插值，绘制地形，浓度三维分布图。求得每种元素的前十个最有可能的污染源的数据制成表格放到附录二。As 的 10 个可能污染源19Cd 的 10 个可能污染源Cr 的 10 个可能污染源Cu 的 10 个可能污染源20Hg 的 10 个可能污染源Ni 的 10 个可能污染源Pb 的 10 个可能污染源21Zn 的 10 个可能污染源的数值模拟，效果和建立的模型基本吻合，污染源出现在极值点附

39、近，但并是通过不是全部是在极值点上。将每种元素的可能污染源画在同一幅图上，进行比较，如图：5.4 问题四解答5.4.1 所建模型优点：（1）前提假设比较符合实际，即污染物的浓度符合正态分布模型经大量小尺度扩散实验证明粗略近似。（2）模式物理概念反映湍流扩散的随机性，其数算比较简单。（3）图表。模型具有坚实的实验基础，扩散参数的求取有比较可靠的经验公式或（4）对基本扩散模式作一些地形修正，便可以直接将其处理一些特殊条件下的大气扩散问题。（5）5.4.2 所建1、第三问所建扩散模式具有形式，因此其数学计算简单，计算量相对较少。模型缺点：模型浓度公式是建立在无风条件下，没有考虑到风场的改变。在风场改

40、变情况下无法使用上述模型。2、模型适用于不太复杂的地形条件，有局限性。5.4.3 还需收集的信息为更好地研究城市地质环境的演变模式，还应收集城区尤其是污染源地区的风速分布，土壤含水率分布，地区 ph 值分布等信息。这里假设收集到了城区的22X(*100m)Y(*100m)Z(*100m)162241061721619153731风速分布，建立模型解决问题。收集风速的信息，更符合实际情况，也可以更好地研究城市地质环境的演变模式，尤其是重金属污染物模式。5.4.4 在考虑风速下的模型建立当风速为u m/s 时，利用连续点源扩散模型分析重金属污染源周边污染物浓度的变化情况。此污染点源是有边界点源，污

41、染点源的实际高度为 H 。以污染点源在地面的投影点为坐标原点以风向方向为 x 轴，铅直方向为 z 轴，与x 轴水平面垂直方向为 y 轴建立三维坐标系，由于扩散过程中浓度在 y、z 轴上的变化分布符合分布，所以下风向的任一点C(x, y, z) 的浓度函数为：22C(x, y, z) A(x)eayebz（12）根据概率统计可以得出方差的表达式为：22y C(x, y, z)dyz C(x, y, z)dz y z 2 200（13）00C(x, y, z)dyC(x, y, z)dz进而源强的积分公式可以根据假设得出：Q uC(x, y, z)dydz把（12）式代入（13）积分可以得出：（1

42、4）a 12 y 212 z 2（15）b 将（12）式和（15）式代入（14）式可以得出：Q2 u y zA(x) （16）最后再将把（15）、（16）式代入（12）式可以得出：y2z2QC(x, y, z) exp 2 u y2 z （17）2y z由上，则得到在考虑风速下的重金属污染物模型，而确定污染源位置，即由污染源周围浓度分布逆推初污染源位置的推导过程同问题三。六、 模型分析23模型评价模型优点综合污染指数：突出污染较重的污染物的作用，考虑了分指数的算数平均值和突出了污染最严重元素的贡献率，克服了平均值法各个污染物分担的缺陷。综合污染指数：用 AHP 法对各分指数对综合污染的影响贡献

43、赋权AHP 法值，能更加全面的反映环境质量。不仅克服了模型一缺陷，同时又通过权重使主要污染因子对环境的突出影响得以体现。灰色关联度：可求出八种重金属元素两两之间哪些关联度最大，在通过数据得到每种重金属元素污染程度严重的地区在哪个功能区后，确定元性便于探讨污染来源。间的相关模型：物理概念清晰，具有很好的可植性，具有很高的计算效率和空间分辨率。既可以得到连续点源出污染源位置。6.1.2 模型缺点扩散方程，又可以由污染源周围浓度分布逆推综合污染指数：未考虑各分指数对综合污染的影响贡献的差异。这种通过直接平均计算得出的综合值只有在所有评价元素超低或比较一致的情况下还不至于偏差太大，如果在不同评价元素污

44、染指数分散的情况下，就会得出错误信息，掩盖真实污染情况。AHP 法综合污染指数：只是从每种重金属含量的角度加以评价，均化了和忽略了重金属毒性差异，实际生态效应。另外，权重的确定问题一直是学者的热点，权重选择根据具体情况差别较大，不具有一定的普遍性。灰色关联度：当在评价指标较多的情况下数据统计量大，易出现很多指标间相关度相似的情况，不利于问题分析。模型：适用于开阔平坦地形不复杂条件下，设下的扩散模式，有局限性。6.2 模型改进续点源在正态分布假1、鉴于综合污染指数和 AHP 法综合污染指数两缺点，在设置不同金属元素权重时，可以通过浓度和毒性相结合来寻找各元素的最佳权重。这样更真实的反映不同金属的

45、污染状况。2、城市土壤重金属污染实际上是一个全面而又复杂的过程，来源可分为两类：一类是土壤本身性质决定的重金属富集，一类是由人类活动造成。故还需对多种(如土壤样品、大气干湿沉降样品、水样、植物样、有机样品等)进行综合分析研究，这样可以更全面反映土壤重金属污染情况。七、 结论1、对城区内的 319 个土壤采样点按照不同功能区对八种重金属元素含量进行统计结果表明，城区的五个功能区中工业区和主干道路区重金属元素相对含量较高。这可能由两方面的原因造成:一是重金属元素的浓度与地形有关，高浓度经常出现于地势较低处，这两个功能区地势相对较低。另一方面，与人类活动有关，这两个地区中工业区排放的废水、废气较多，

46、而主干道路区中车辆尾气排放较多 2、从污染程度上来看，城区内表层土壤的重金属污染顺序为: HgCu Zn Cd PbCrAs;从污染区域来看，污染程度顺序为:工业区主干道路区生活区24公园绿地区山区。3、本文中运用了灰色关联度分析并发现了在城区内表层土壤中，除了 Hg 与 Zn两种元素间关联度较低外，八种重金属元素间均存在较为显著或极显著正相关关系。另外，As 与 Cr、Hg、Ni 有极显著的相关性;Cd 与 Pb 有显著关系;Ni 与 As、 Cr 之间均存在极显著的关系。Pb 和 Cu、Zn 和 Cu 呈显著的相关性。4、对重金属污染物的途径进行分析，经过排除确定为大气扩散。通过对条件的分

47、析，最终利用连续点源扩散模型确定特征和浓度扩散特性，再由污染源周围浓度分布逆推出污染源位置，并且污染源位置不一定是高污染地区的浓度极值点。5、增加风速或土壤含水率、ph 值等条件，都可以提高模型与实际情况的贴近程度，更好地研究城市地质环境的演变模式，尤其是污染物的扩散模式。参考文献12新，环境污染化学，：化学工业，2003，吉林市城市表层土壤中重金属污染状况研究，吉林大学学位论文，20093，崇明岛土壤重金属污染的空间分析，华东师范大学，2009及污染评价4社，2005 5，环境系统模型及数值模拟，：中国环境科学，土壤重金属污染的非线性可拓综合评价，土壤， 2004,36(2):151-156

48、,20046才，大气污染扩散的烟羽模型及其 GIS 集成研究J，环境监测管理与技术，2008，20（5）：17-1925附录一问题 3 程序%yuan_as.m求极值周围的污染源load load loadloadx.txt y.txt z.txtas.txt%读入数据，若算其他重金属，替换就好me=as;x=x./100;y=y./100;z=z./100;% 将坐标缩小100倍，使后面的计算（标准差，穷举）方便 m_me l=max(me);x0=x(l);y0=y(l);z0=z(l);me0=m_me;me00=m_mexl=x0-20 x0-10 x0+10 x0+20;% 以最大值

49、为中心分别向外延伸2000米，为大区域，求得大区域内的测量点yl=y0-20 y0-10 y0+10 y0+20;% 以极大值为中心的2000米以内的正方形为小区域，对其海拔和浓度进行三次样条插值，以便寻找污染源n=zeros(1,319); for i=1:319if(x(i)xl(2)if(x(i)yl(2)if(y(i)0)no(j)=i;% j=j+1;endendfor i=1:n_no2=no;yyy=no;zzz=no;meee=no;j=1;测量点的序号(i)=x(no(i);yyy(i)=y(no(i);zzz(i)=z(no(i);meee(i)=me(no(i);% 将小

50、区域测量点的x,y,z,浓度重新组end新的矩阵n=zeros(1,319); for i=1:319if(x(i)xl(1)26if(x(i)yl(1)if(y(i)0)no(j)=i;%j=j+1;测量点的序号endendfor i=1:n_nos2(i)=(x(no(i)-x0)2+(y(no(i)-y0)2+(z(no(i)-z0)2;xx(i)=x(no(i);yy(i)=y(no(i);zz(i)=z(no(i);mee(i)=me(no(i);% 将大区域的测量点的x,y,z,浓度重新组end新的矩阵xxi=xl(2):1:xl(3);yyi=yl(2):1:yl(3);% 将小

51、区域分划网格 zzi=griddata(xx,yy,zz,xxi,yyi,V4);% 在小区域内进行关于海拔的三次样条插值 c0=mee.*exp(s2/400);max_c0=max(c0);std_c0=std(c0);yuan=x0 y0 z0 me0;nn=11;mm=11;f=1; %定义判断指针，判断污染源是否为已知极值点for i=1:10%向外分层循环，确定每一层的个点的相对坐标 n1=ones(1,i*4)*i;m1=ones(1,i*4)*i;for j=1:2*in1(j)=i+1-j;m1(j+2*i)=i+1-j;endn2=-1*n1;m2=-1*m1; n=n1

52、n2;m=m1 m2;yuan_x=x0+n;yuan_y=y0+m;yuan_z=zzi(11+n,11+m); %将相对坐标转化为真实坐标 for k=1:8*i%对第i层的格点穷举，找出最优点for ii=1:n_nos2(ii)=(x(no(ii)-yuan_x(k)2+(y(no(ii)-yuan_y(k)2+(z(no(ii)-yuan_z(k)2;%距离平方endc01=mee.*exp(s2./400);max_c1=max(c01);std_c1=std(c01);%求c0，最大值，标准差27if(std_c1=me0)x0=yuan_x(k);y0=yuan_y(k);z0=yuan_z(k);me0=max_c1;nn=n(k)+11;mm=m(k)+11;%判断最优点 f=0;%若更改则污染源不是极值点endendend end if(f)fpr elsefpr xx=endf(原极值点)f(非原极值点)0;yy=yy y0;zz=zz z0;mee=mee me0;zzi=griddata