2011高教社杯全国大学生数学建模竞赛_重金属污染分析

承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。

如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）：
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：
参赛队员(打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人(打印并签名)：
日期： 2011 年 9 月 12 日赛区评阅编号（由赛区组委会评阅前进行编号）：
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
城市表层土壤重金属污染分析
摘要
随着城市规模的不断扩大和工业化进程的加快,土壤污染问题日益严重。

本文通过对城市表层土壤重金属污染的情况进行分析，借助于单向污染指数法和内梅罗综合污染指数法，确定污染程度。

我们根据重金属元素的传播特征，建立重金属污染物的指数扩散模型，利用指数拟合的方法，确定污染源的位置。

对于问题一，我们根据所给数据，利用Matlab软件做出了各个重金属元素在该城区的空间分布图。

由图可以看出该城区内不同区域重金属的污染程度。

再通过建立单项污染指数模型和综合污染指数模型，求得不同重金属元素在不同功能区的污染程度，如
对于问题二，经过计算分析可得：Cd、Cu、Hg、Pb、Zn主要由2类区和4类区产生，As主要由2类区和1类区产生，Cr、Ni主要由4类区和1类区产生。

对于问题三，经过均衡性检验，海拔对结果影响较小，忽略不计。

重金属污染物以污染源为圆心，以递减的速率向周边传播，符合指数扩散特征。

因此，我们建立指数扩散模型来模拟重金属的扩散。

各污染源位置坐标如下（单位：m）：
As：(3913,7475)、(12660,3051)、(18938,10357)
Cd：(3788,2740)、(9315,6403)
Cr：(3688,5180)
Cu：(3908,2528)
Hg：(2365,2787)、(13450,2441)
Ni：(3395,5886)
Pb：(4707,5096)、(1936,2888)
Zn：(9168,4816)、(13768,9659)
对于问题四，可看成是对城区土壤污染程度的预测，应收集当地的河流路径数据、常年盛行风的方向及风速、通过植物每年平均流入生物界的数量这些信息，建立多元素影响下的重金属元素传播模型。

可借鉴高斯烟羽模型，预测该城市地质环境的演变。

最后，对建立的模型进行误差分析，所得数据较为精确，基本符合实际情况。

关键词：单向污染指数内梅罗综合污染指数局部指数扩散模型均衡性检验高斯烟羽模型
一、问题重述
由于重金属化学性质和生态效应的复杂性，土壤重金属污染一直是环境科学的重点研究课题。

我国对土壤重金属的研究大多数集中于农业土壤的研究。

随着城市化进程的加快，作为受人为活动剧烈影响的城市土壤也成为研究的热点领域。

对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。

在此基础上，对某城市城区土壤地质环境进行调查，采样记录每个样本所含的多种化学元素的浓度数据，并在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。

利用这些数据，需要解决的问题有：
1、确定8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。

2、通过数据分析，说明重金属污染的主要原因。

3、分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。

4、分析所建立模型的优缺点，并确定还应收集什么信息才能更好地研究城市地质环境的演变模式，有了这些信息后，建立模型解决问题。

二、问题分析
2.1问题1的分析
根据对所给数据的分析，在五类区域中，各个重金属的污染程度不同。

在研究8种主要重金属元素的分布的过程中，若只画出重金属浓度的分布图，只能大概看出各个区域的污染程度，不能清晰地看出区域和污染程度的具体对应关系。

因此，还需要对具体数据进行定量分析，选择单项污染指数法和综合污染指数法来评定污染的程度。

为了能够更清楚的表现重金属污染的情况，我们将分布图和数据图表结合起来。

在此过程中，由于数据数量太大，而且分析的是污染程度的大小，所以，我们可以将每种重金属浓度降序排序，选取浓度最高的若干个点，标记在图上，由此分析各个区域的不同重金属污染程度。

2.2问题2的分析
要找到重金属污染的主要原因，需要找到每种重金属元素污染的范围内各个功能区所占的比例，以此来确定污染的主要原因。

通过查找相关资料，可以得知重金属污染的原因有：工业污染、生活垃圾污染、交通污染、农业污染（农药、化肥和塑料薄膜使用、污水灌溉、污泥施肥）、工矿污染等。

本题中，需要对数据进行一定的处理，找到重金属污染的主要原因。

首先，对数据进行分类，我们按照功能区的不同将数据分为5组，可以得到每个功能区内各个重金属元素浓度的平均值，再利用单项污染指数法把数据进行标准化，将浓度按比例转化成单项污染指数，就可以看出重金属元素在各个功能区的分布。

由于平均单项污染指数并不能反映各个功能区内到底有多少个被污染，因此，我们可以筛选出达到每个重金属污染浓度的数据，进而按照每种重金属元素浓度超过标准值的比例来判定重金属污染的主要区域。

2.3问题3的分析
重金属污染物传播特点是由污染源向周边逐渐扩散。

要想确定污染源的位置，必须先确定重金属污染区域。

由问题1的图像可大致确定污染范围，再根据此区域的数据，建立污染物指数传播模型，可以求得重金属浓度与位置的函数关系，要想求得浓度最大值，可对所求表达式求偏导，偏导数为零的点即为浓度最大值点，并且需要确定污染源的范围，进行误差分析。

2.4问题4的分析
此问即是对之前所建立模型的改进，通过分析查找更多的影响重金属元素传播的因素，建立多元素影响下的重金属元素传播模型。

三、 模型假设
1、假设污染源全部在此城区内，不考虑城区外污染源的影响。

2、假设重金属传播过程中，只受周围功能区的影响，不会转移到生物界或大气中；
3、假设每种重金属都是可以在外界因素下可以自由传播的，不存在重金属元素在一个地区堆积，只有传播能力的不同；
4、假设在研究重金属元素传播的过程中，之前无重金属元素的积累，都是由污染源向周围传播，且分布是均匀连续的；
四、定义与符号说明
i P ：区域重金属i 的单项污染指数；
i C ：重金属i 含量的实测值(/)/(/)g kg ng kg μ； i S ：重金属i 含量的评价起始值(/)/(/)g kg ng kg μ；
P 综
：综合污染指数（综合反映各污染物对区域土壤的不同作用）； av P ：所有单项污染指数的平均值；
max P ：土壤环境中各单项污染指数中的最大值； 0Q ：进行指数拟合时，选用的浓度的初始值；
s :距离污染源的距离；
(,,)Q x y z ：点(,,)x y z 处的污染物浓度；
a ：指数扩散模型的系数； v ：风速；
x δ：正风向扩散系数；
y δ：偏风向扩散系数； z δ：垂直风向扩散系数。

五、模型的建立与求解
5.1准备工作
数据的处理：由于数据量比较大，需要进行初步处理。

分别按照不同功能区、各种重金属浓度大小分类，便于之后的问题求解。

5.2问题1模型建立与求解
首先，将8种主要重金属元素在该城区的空间分布用Matlab 程序（详见附录1）画出图形。

由于数据数量太大，我们将每种重金属浓度降序排序，选取浓度最高的50个点，标记在图上。

由于数据是按照每平方公里1个采样点进行采样的，因此可以确定此城区的大致轮廓，我们只考虑此城区内的污染源，由图像所看出的其他不属于本城区的污染源暂不考虑。

参考城区轮廓如下图1所示：
图1 城区大致轮廓
重金属元素在该城区的空间分布图如下： 注：图中颜色越深表示污染越严重。

+12(3(4(5( ：类区(生活区)
：类区工业区)：类区山区)：类区主干道路区)☆：类区公园绿地区)
图2 金属元素As的空间分布
图3 金属元素Cd的空间分布
图4 金属元素Cr的空间分布
图5 金属元素Cu的空间分布
图6 金属元素Hg的空间分布
图7 金属元素Ni的空间分布
图8 金属元素Pb的空间分布
图9 金属元素Zn的空间分布
由上面的图可知，As 元素有三个污染源，污染是由主干道路区和工业区造成的；Cd 元素有两个污染源，污染是由主干道路区造成的；Cr 元素有一个污染源，污染是由主干道路区造成的；Cu 元素有一个污染源，污染是由工业区造成的；Hg 元素有三个污染源，污染是由主干道路区造成的；Ni 元素有一个污染源，污染是由主干道路区造成的；Pb 元素有两个污染源，污染是由生活区和工业区造成的；Zn 元素有两个污染源，污染是由主干道路区造成的。

由上面分析大致得知：各个功能区分布在不同的浓度范围。

工业区和主干道路区污染很严重。

As 污染严重的范围内，工业区和主干道路区较多；Cd 污染严重的范围内，工业区和主干道路区较多；Cr 污染严重的范围内，生活区和主干道路区较多；Cu 污染严重的范围内，工业区和主干道路区较多；Hg 污染严重的范围内，工业区和主干道路区较多；Ni 几乎无污染；Pb 污染严重的范围内，工业区和生活区较多；Zn 污染严重的范围内，工业区和主干道路区较多。

为了更加确定的描述该城区内不同区域重金属的污染程度，我们引用了两种评价标准。

目前对于土壤污染综合评价的方法较多，主要有污染指数叠加法、内梅罗综合污染指数法、权重法、最大值法等。

通过比较分析，我们选取了比较常用的单项污染指数法和内梅罗综合污染指数法对重金属浓度进行了相应的比对。

（1）、单项污染指数法。

单项污染指数法的计算公式为
i
i i
C P S =
若01i P <≤为无污染，12i P <≤为轻微污染，23i P <≤为轻度污染，35i P <≤为中度污染，5i P >为重度污染。

同时，我们取重金属i 含量的评价起始值i S 为背景值范围中最大值。

用excel 求得各个重金属元素的单项污染指数后，可求得每个功能区的所有单项污染指数的平均值，结果如下表1所示：
由表一我们可以明显得看出不同重金属在各个功能区的污染程度，Hg 和Cu 污染很严重，与图像上反映的结果接近，可知，两种方法确定的重金属元素在该城区的空间分布是正确的。

（2）、内梅罗综合污染指数法。

内梅罗综合污染指数法的计算公式为：
P 综
内梅罗综合污染指数的分级标准见下表2： 表2 内梅罗综合污染指数的分级标准
根据内梅罗综合污染指数法的计算公式，带入重金属的av P 和max P 的值，求得P 综
如下表3：
表3 P 综
的值 由表2和表3的对照可知：该城区内不同区域重金属的污染程度为2类区>4类区>1
类区>5类区>3类区。

5.3问题2模型建立与求解
首先，按照1i P 筛选出每种重金属元素及其对应的功能区，其中共有319个样本
经过excel 表格筛选，将各个功能区污染的个数求出（详见附录2），并计算出功能
区内受到污染的比例
功能区污染的个数
功能区总个数
，结果如表5：
表5 各功能区重金属污染比例
由上表可以看出：工业区污染最严重，生活区和交通区大部分也受到污染，土壤质量最好的为山区。

可以推测：As 、Pb 元素污染是由于工业污染造成的；Cd 、Cu 、Hg 元素污染是由
工业污染和交通污染造成的；Cr 元素污染是由生活污染造成的；工业废气、废水、废渣，交通尾气，生活垃圾中含有大量重金属污染物，是重金属污染的主要原因。

5.4问题3模型建立与求解
通过建立局部指数扩散模型，得到土壤中重金属传播是按照指数传播的。

局部指数扩散模型的基本形式为：
0()as Q s Q e -=
在上式中：
s =
通过均衡性检验，海拔对污染物传播的影响较小，忽略不计。

因此，公式转化为：
0(,)Q x y Q e
-=为了更加方便的利用Matlab 中lsqcurvefit 函数求解多元非线性拟合，我们把上式两边取对数，得如下公式：
0ln (,)ln()Q x y Q =-由各重金属元素的空间分布图可近似确定污染源的个数和位置，以此来划定范围。

在范围内筛选数据，选出满足要求的点，做局部拟合，利用Matlab 编程（详见附录3）可求得污染源的位置00(,)x y ，结果为：
As ：(3913,7475)、(12660,3051)、(18938,10357) Cd ：(3788,2740)、(9315,6403) Cr ：(3688,5180) Cu ：(3908,2528)
Hg ：(2365,2787)、(13450,2441) Ni ：(3395,5886)
Pb ：(4707,5096)、(1936,2888) Zn ：(9168,4816)、(13768,9659)
由此可以看出，某些重金属元素的污染源接近，如Cr(3688,5180)和Ni(3395,5886)。

可以推测，重金属污染具有一定的相关性。

据此，我们做了各个重金属元素的相关性检验，得出结果如表6：
由上表可明显看出重金属元素之间的相关性联系。

Cu 和Cr 之间， Ni 和Cu 之间，
Pb 和Cu 之间，Zn 和Pb 之间的关系比较密切。

5.5问题4模型建立与求解
根据题目中所给数据建立的模型，由于缺乏某些必要的信息数据，如风速、时间等，导致考虑重金属污染物传播特征时，只能兼顾地理位置对重金属元素扩散的影响，所建模型不是特别理想。

要想更好地研究城市地质环境的演变模式，首先应考虑土壤中重金属污染物来源。

土壤中重金属的来源是多途径的，首先是成土母质本身含有重金属，不同的母质、成土过程所形成的土壤含有重金属量差异很大。

此外，人类工农业生产活动，也造成重金属对土壤的污染。

大气中重金属沉降、农药、化肥和塑料薄膜使用、污水灌溉、污泥施肥、含重金属废弃物堆积、金属矿山酸性废水污染等都是影响重金属污染物浓度的原因。

其次，还应考虑重金属污染物可能传播的途径，影响因素如雨水、河流、植物界吸收、风的传送即盛行风传播方向、机动车密度等等。

综合上述各种因素，经过筛选，研究城市地质环境的演变模式模型，还应收集当地的河流路径数据、常年盛行风的方向及风速、通过植物每年平均流入生物界的数量这些 信息。

通过收集到的信息，模拟高斯烟羽模型，综合考虑影响因素，可以建立多元素影响下的重金属元素传播模型：
2
22222112220
(,,)4x x y z y e z
x y z
Q Q x y z e
e
v δδδπδδδ-
⎡⎤
⎢⎥-⎢⎥--
⎢⎥⎢⎥-⎢⎥⎣
⎦
=
由分析知，海拔越高，污染物扩散越慢，因此在高斯烟羽模型的基础上，我们添加了一项关于z 变量表达式，若得到相关数据，模型即可得到验证，使模型的精确性大大提高。

六、模型评价及推广
6.1、模型优点
1)、通过大量数据的分析，更具有说服力，并且选取了具有代表性的部分点来代替大量的数据，简化了计算过程。

2)、本文利用单项污染指数法和综合污染指数法分析，将图形和图表数据结合在一起，共同描述该城区内不同区域重金属的污染程度，精确度极高。

3)、求解污染源问题即为利用所有数据进行拟合，求得浓度最大值点，我们建立了指数模型，优于二次函数模型。

6.2、模型缺点
1)、只考虑了地理位置对污染物传播的影响，而事实上还有很多影响因素，由于忽略部分的外界因素，会导致模型与实际情况有些差距。

6.3模型推广
由于城市的生产活动和日常生活向环境中输入巨量三废, 严重污染了城市及其周围的土体、水和空气, 改变了城市的地球化学环境。

我们建议，利用所建模型，建立实时监测预报模型，实时监测预报方法是集数据自动采集、处理、建模、预测预报及信息的短程、中程及远程传送于一体的集成技术与方法, 可以对城市地质环境和环境监测预报作出贡献，并利用计算机技术、电子技术和人工智能技术作为研究的支柱。

通过研究出的模型，确定风向、水流等对重金属污染物传播的影响，为污染源的处理提供了参照，如工厂的选址、工厂废水废渣的处理方案、生活垃圾的处理等等。

七、参考文献
[1]姜芝萍，杨俊衡.城市重点污染区土壤重金属污染评价标准探讨.安全与环境工程，2010，第17卷第一期：57页-64页。

[2]柴世伟,温琰茂等.广州市郊区农业土壤重金属污染评价分析.环境科学研究，2006,第19卷第4期:138页-142页。

[3]谭永基，数学模型，上海：复旦大学出版社，1999。

[4]赵静，数学建模与数学实验,北京：高等教育出版社，2009。

[5]王长海.污染物扩散的一种计算模式.交通环保，1997年第1期，9页-13页。

[6]汪家权,刘万茹等.基于单因子污染指数地下水质量评价灰色模型.合肥工业大学学报( 自然科学版),2002年10月,第25卷第5期:697页-702页。

八、附录
1、画重金属元素空间分布图的Matlab 程序（M 文件）：
注：由于数据量太大，原始附件中的数据便不再列出，以下程序中的A 矩阵是按照功能区由小到大排序列出的，如下矩阵所示：
()()()(/)(/)...(/)...x m y m m As g mg Cd ng mg Zn g mg A μμ⎡⎤=⎢⎥⎣⎦ 编号海拔功能区
Matlab 程序（M 文件）： x=A(:,2); y=A(:,3); z=A(:,6);
[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x),200)',linspace(min(y),max(y),200),'v4');%插值
pcolor(X,Y,Z);shading interp %伪彩色图 contourf(X,Y,Z,6) hold on
colormap (flipud(gray)); hold on
A1=[4777 4897 8 1 11.45 9277 16148 18 1 10.97 11121 16432 23 1 9.81 3573 6213 5 1 9.39 4043 1895 14 1 9.17 5382 3012 50 1 8.90 12153 12336 16 1 8.67 4592 4603 6 1 8.50 5101 4080 13 1 8.23 14298 7418 36 1 8.06 6534 5641 6 1 7.84 5636 133 17 1 7.78 4742 7293 9 2 21.87 4948 7293 6 2 18.38 1647 2728 6 2 14.08 9212 11305 5 2 10.74 5567 6782 7 2 10.53 7776 10613 9 2 10.27 2383 3692 7 2 9.62 12641 9560 11 2 9.58 5438 3994 10 2 9.35 7106 9467 44 2 8.90 5006 8846 6 2 8.67 8446 11200 4 2 8.23
8866 13143 3 2 8.23 10856 14727 41 2 8.23 9036 17538 3 2 7.78 18738 10921 53 3 10.99 15007 5535 70 3 9.17 20983 15862 93 3 8.00 19968 12961 42 3 7.63 18134 10046 41 4 30.13 12696 3024 27 4 23.72 6869 7286 18 4 16.58 3299 6018 4 4 9.84 14065 10987 25 4 9.13 2883 3617 15 4 8.94 2933 1767 7 4 8.72 17005 7212 33 4 8.67 17087 11933 43 4 8.50 10547 9591 32 4 8.23 4020 2990 27 4 8.00 74 781 5 4 7.84
8403 1075 6 4 7.78 9095 16414 29 5 11.68 10352 17133 31 5 10.74 7653 1952 48 5 9.35 4153 2299 73 5 8.72 14318 13569 30 5 8.23 3762 2170 30 5 8.06 ];
x1=A1(1:12,1);
y1=A1(1:12,2);
plot(x1,y1,'k+')
hold on
x2=A1(13:27,1);
y2=A1(13:27,2);
plot(x2,y2,'w*')
hold on
x3=A1(28:31,1);
y3=A1(28:31,2);
plot(x3,y3,'gs')
hold on
x4=A1(32:44,1);
y4=A1(32:44,2);
plot(x4,y4,'r>')
x5=A1(45:50,1);
y5=A1(45:50,2);
plot(x5,y5,'yp')
2、污染的功能区个数分布
3、求污染源的Matlab的M文件
注：所求图像程序相同，数据不同，因此以As的污染源为例。

xdata=[2486 5481 3299 3573 4741 5567 4742 4948 5291 5635 5394 5375 5006 5734];
ydata=[5999 6004 6018 6213 6434 6782 7293 7293 7349 7965 8631 8643 8846 9659 ]; zdata=[5.51 7.41 9.84 9.39 3.3 10.53 21.87 18.38 3.69 6.41 5.31 4.09 8.67 5.41]; data=[xdata; ydata];
zdata=log(zdata);
format long;
a0=[4 0.00009 4000 7000];
[a, resnorm] = lsqcurvefit(@myfun,a0,data,zdata)
function f= myfun(a, data);
x = data(1,:);
y = data(2,:);
f=a(1)-a(2)*sqrt((x-a(3)).^2+(y-a(4)).^2);。