数模全国赛论文.docx
《数模全国赛论文.docx》由会员分享，可在线阅读，更多相关《数模全国赛论文.docx（24页珍藏版）》请在冰豆网上搜索。
数模全国赛论文
城市表层土壤重金属污染分析
摘要：
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。
对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
土壤重金属污染日益严重，重金属可通过吸入，附入和食物的富集对人体健康和环境产生较大影响。
对重金属的污染研究对于综合治理环境，城市规划等有重要的指导意义在本文中将研究的城市分为五大功能区分别为生活区，工业区，山区，交通区五大功能区。
利用Matlable插值计算五大功能区的横纵坐标海拔和浓度的数据，得到区域平面图和8种主要重金属元素在该城区的空间分布。
运用单因子评价中的分指数法逐一计算每种重金属元素的污染分指数，然后运用内梅罗（N.L.Nemerow）综合污染指数法得到总体的污染指数，参照农业部行业标准《中华人民共和国国家标准GB15618-1995土壤环境质量标准》的分级法得到污染程度。
通过对比数据，利用相关矩阵得到8大金属元素之间的相关程度，从而得到污染的来源主要是工矿企业和汽车尾气。
运用特异函数得到X方向和Y方向上的金属元素浓度变化趋势。
使用二维差值得到污染源的大体位置，即图中等浓度线最密集区域。
污染源位置的确定对于环境整治有重要指导作用。
最后，本文讨论了本模型的改进方向，和改进模型。
关键词：
土壤重金属；污染；内梅罗综合污染指数；二维差值
一、问题的背景和重述
1.1问题背景
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。
对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。
为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10厘米深度）进行取样、编号，并用GPS记录采样点的位置。
应用专门仪器测试分析，获得了每个样本所含的8种化学元素的浓度数据。
另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中8种元素的背景值。
1.2问题重述
（1）给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
（2）通过数据分析，说明重金属污染的主要原因。
（3）分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
（4）分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？
有了这些信息，如何建立模型解决问题？
二、模型假设
考虑到现实中重金属元素的扩散收到多方面因素的影响，为使计算简化，先对模型进行简化。
（1）每一采样点获得的数据都可靠，不需剔除个别数据。
（2）当地的水文气象稳定，无重大地质变化。
三、符号说明
土壤中污染物的污染分指数；
内梅罗（N.L.Nemerow）综合污染指数；
土壤中污染物实测含量；
污染物的评价标准；
污染物的背景值；
污染物的标准差；
污染物中污染指数最大值；
四、模型分析
4.1.1八种主要重金属元素在该城区的空间分布
五大功能区的地理坐标数据相对散乱，数据间距互不相同，但与重金属含量数据一一对应。
可以对地理坐标数据和土壤重金属含量数据进行二维插值，以得到区域平面图和8种主要重金属元素在该城区的空间分布。
4.1.2该城区内不同区域重金属的污染程度
各采样点所获得的重金属数据，基本能够反映出城区的重金属污染情况，可以运用单因子评价中的分指数法逐一计算每种重金属元素的污染分指数。
由于单因子污染指数只能分别反映单个污染物的污染程度，不能全面、综合地反映多种污染物的整体污染水平，故可采用内梅罗（N.L.Nemerow）综合污染指数法对各区域的污染程度进行综合评价。
4.2重金属污染主要原因
一般地认为，我国的污染主要来自工矿企业、交通运输业，居民造成的生活垃圾污染影响较小。
而工厂所产生的重金属污染受其原料、产品、工艺决定，往往集中在几种互相关联的元素。
同样地，交通运输业所产生的重金属污染决定于成品油。
所以可以用因子分析法寻找重金属污染的原因。
4.3重金属污染物的传播特征并确定污染源
地质统计学中的变异函数能够较好地描述区域化变量的空间分布结构性与随机性特征，从而使针对区域化变量的空间变异性分析得以实现。
重金属浓度数据相对与地理数据一一对应，可在地图上绘出重金属污染物等浓度线，并与变异函数结合，以判断污染源的位置。
五、模型的建立与求解
5.1.1八种主要重金属元素在该城区的空间分布
已获得的离散地理坐标数据和对应的重金属浓度数据，与整体城区面积相比数据量太小，故可以通过插值的方法，去寻取近似函数。
运用计算机绘出城区的5种功能区区域分布图（见图1）和8种重金属污染物浓度的区域分布图（见图2）
图1功能区分布图
图2不同金属元素浓度平面分布图
5.1.2该城区内不同区域重金属的污染程度
土壤环境质量现状评价方法与大气、水质的现状评价方法基本相同，均通过计算污染指数进行。
土壤环境质量单因子评价通常采用分指数法，即通过逐一计算土壤中各主要污染元素的污染分指数，已确定污染程度。
土壤污染分指数计算公式为：
（1）
其中
。
在实际情况中，多种污染物往往同时污染某一区域或土壤，单因子评价难以表示它们的整体污染水平，因此需要一种同时考虑多种污染物综合污染水平的多因子评价方法。
即将单因子污染指数按一定方法综合，较全面反映污染状况。
由于污染指数
消除了量纲，为综合指数的获取提供了方便。
内梅罗（N.L.Nemerow）综合污染指数法是国内普遍采用的综合评价方法之一。
可它以全面反映各污染物对土壤的不同作用，并突出高浓度污染物对环境质量的影响。
[1]
内梅罗（N.L.Nemerow）综合污染指数法公式为：
（2）
通过公式
（1）
（2）得到数据（见表1）。
土壤质量分级是土壤质量评价的基本内容。
以单项污染指数和综合污染指数，按照农业部行业标准《中华人民共和国国家标准GB15618-1995土壤环境质量标准》的分级法将土壤污染划分为安全级、警戒级、轻污染、中污染、重污染五个等级（见表2）。
表1各功能区的污染指数
污染指标
区域
As
Cd
Cr
Cu
Hg
Ni
Pb
Zn
内梅罗综合
污染指数P
1生活区
1.16
1.53
1.41
2.42
1.82
0.92
1.61
2.44
2.09
2工业区
1.34
2.07
1.09
6.25
12.59
0.99
2.16
2.87
9.28
3山区
0.75
0.80
0.79
0.85
0.80
0.78
0.85
0.76
0.82
4交通区
1.06
1.89
1.18
3.05
8.76
0.89
1.48
2.50
6.46
5公园
绿地区
1.16
1.48
0.89
1.48
2.25
0.77
1.41
1.59
1.87
表2土壤污染分级
等级
内梅罗综合污染指数
污染等级
1
P≤0.7
清洁（安全）
2
0.7≤P＜1.0
尚清洁（警戒）
3
1.0≤P＜2.0
轻污染
4
2.0≤P＜3.0
中污染
5
P≥3.0
重污染
由表1，表2可见，污染程度由大到小顺序为工业区>交通区>生活区>公园绿地区>山区。
生活区属于中污染程度，工业区属于极重污染，山区属于尚清洁警戒状态，交通区属于较重污染，公园绿地区属于轻污染程度。
5.2重金属污染主要原因
这里我们运用因子分析法得出重金属污染的主要原因。
因子分析法从变量的相关矩阵出发将一个m维的随机向量分解成低于m个且有代表性的公因子和一个特殊的m维向量，使其公因子数取得最佳的个数，从而使对m维随机向量的研究转化成对较少个数的公因子的研究。
设有n个样本，n个指标构成样本空间X，
，
原始数据的标准化，标准化的公式为，
为第i个样本的第j个指标值，而
和
分别为j指标的均值和标准差。
标准化的目的在于消除不同变量的量纲的影响，而且标准化转化不会改变变量的相关系数。
首先得到八种重金属含量数据的相关系数矩阵，如表3所示。
可见，Cr和Ni的相关性最好，相关系数最大，为0.7158，其次为Pb和Cd，相关系数为0.6603，以下依次是Cr和Cu,Pb和Zn的相关性较好，相关系数分别为0.5316和0.4937，Ni和Zn的相关系数为0.4364，其它元素之间的相关性并不是很好。
从成因上来分析，相关性较好的元素可能在成因和来源上有一定的关联。
表38种金属元素的相关矩阵
As
Cd
Cr
Cu
Hg
Ni
Pb
Zn
As
1.0000
Cd
0.2547
1.0000
Cr
0.1890
0.3524
1.0000
Cu
0.1597
0.3967
0.5316
1.0000
Hg
0.0644
0.2647
0.1032
0.4167
1.0000
Ni
0.3166
0.3294
0.7158
0.4946
0.1029
1.0000
Pb
0.2899
0.6603
0.3828
0.5200
0.2981
0.3068
1.0000
Zn
0.2469
0.4312
0.4243
0.3873
0.1958
0.4364
0.4937
1.0000
从表1中可以看出Cr和Ni的浓度分布基本一致，这与理论得出的一样，说明Cr和Ni的污染来源确实有很大关联。
Cr和Ni的分布成面积形，集中位于城市的西南部的小块区域，这块区域主要是交通区，生活区和工业区，而生活区最严重，交通区次之，说明Cr和Ni污染的原因是含Ni、Cr元素的以及生活废水的排放。
Cd和Pb在来源上关联较密切，在空间分布上近似可认为是一个带状的污染源，呈带状分布。
这主要因为Pb大部分来自市交通区汽车尾气的排放，在主干道路区和在公园绿地区与工业区交界区有明显的富集，形成一个高值区，说明该处交通发达。
根据文献，城市中的铅一般主要来源的厂矿企业部门为颜料厂，冶炼等行业的废水，橡胶和农药厂，而镉主要来源于颜料行业和石油化厂。
可以得出一个基本结论：
市内交通尾气的排放和汽车轮胎的磨损很可能是土壤铅污染的基本来源，表层土壤中的Cd含量主要富集在工业区说明市镉主要来源可能是是颜料行业和石油化厂。
[2]
As的污染较轻，各功能区的污染程度差不多。
富集区主要在工业区，所以As污染的主要来源可能是工业废水和原料的污染。
高含量Cu主要集中在市西南部工业区和交通区，为小面积型污染，主要来源可能是化工行业、塑料、橡胶和印染行业的废水排放，以及城市商业活动和交通累加到土壤中的铜。
Hg是该市污染最严重的重金属元素之一，这在我国城市中非常常见，尤其是在我国北方城市中表现得尤其突出。
其中一个最大的原因是燃煤污染。
图2中可出，汞污染也属于面积型污染，主要集中在三块小面积区域，工业占了相当大比重，且工业区Hg污染极为严重，交通区也十分严重。
这有最可能原因是该市的燃煤工业集中在这三个区域，比如冶炼企业。
还有化工企业的废渣废气，废水。
另外，Hg污染与汽车尾气的排放也有很大关系。
Zn在该市的污染大体是中度污染。
生活区，工业区和交通区的富集程度差不多。
说明Zn主要来源于工业三废，工矿企业，居民生活垃圾和交通尾气和轮胎的磨损。
[3]
5.3重金属污染物的传播特征并确定污染源
通过对地理坐标数据和污染物浓度数据进行二维线性插值，获得等浓度分布图。
地质统计学中的变异函数能够较好地描述区域化变量的空间分布结构性与随机性特征，从而使针对区域化变量的空间变异性分析得以实现。
本文即运用变异函数对污染源进行确定。
变异函数被定义为：
（3）
式中，h为距离滞后，或称步长，E表示数学期望，
为在位置x处的变量值，
为在与位置x偏离h处的变量值。
随着距离段的变化，可计算出一系列的变异函数值。
以h为横坐标，
为纵坐标作图，使得到了变异函数图。
从计算公式可见，变异函数实际上是一个协方差函数，是同一个变量在一定相隔距离上差值平方的期望值。
差值越小，说明在此距离段上该变量值的相关性越好；反之亦然。
[4]
（3）式是较为严格的数学定义，同时适用于空间上连续分布的变量。
但在实际工作中，采样点常常是离散的，为此将上式改写为：
（4）
式中，
为距离等于h的点对数，
为处于点
处变量的实测值，
为与点
偏离h处变量的实测值。
实际上，式（4）将式（3）中的期望符号改为求算术均值，便可实际应用了。
此变异函数被称为实验变异函数。
出于地球化学数据的波动性常常比较大，同一个距离段内所有样点对时间差值的平方可能较为离散，常常不会服从正态分布，那么，上式中采用求算术均值的做法就会有偏差。
计算出的实验变异函数值，只是不同方向上的一些不连续的点，还必须经过理论模型拟合得出理论模型，才能较好地分析区域化变量的空间变异性。
[5]
有了变异函数，就可以应用地质统计学的理论和方法对参数的空间分布进行研究。
研究主要有两个方面的内容：
一方面是应用变异函数对参数的空间分布进行结构分析和变异性分析：
另一方面是应用结构分析的结果，结合二维线性插值绘制的等值线图，分析重金属元素的传播特性和污染源的位置。
[6]
如图3是8种元素的等浓度分布图，从这个图中我们便可以很好的观察出各种重金属元素的大体分布形状和区域，对比不同功能区的分布特点得出污染物的大体来源。
Zn
Pb
Ni
Hg
Cu
Cr
As
Cd
图38种元素的等浓度分布图
通过图3可以很明显得看出8种重金属元素污染源位于等浓度线最密集的区域分布情况和传播趋势。
Zn富集区域主要成小面积状，污染源很有可能位于富集最高的地方有两处（3700,9700），（9100,4500）。
Pb主要成环状分布,有一块富集程度相当高的区域，污染源位于点（18400,10100）附近。
Ni主要成环状分布,有一块富集程度相当高的区域，污染源位于点（21600,11400）附近。
Hg主要集中在一小块区域，相对集中，污染源在点（3900,5100）附近。
Cu就是富集在一块区域，污染源在（2400,3500）附近。
Cr富集在三块小面积区域，污染源在（15200,9300），（13700,2200），（1800,2900）附近。
Cd主要成环状分布,有一块富集程度相当高的区域，污染源在点（3300,5900）附近。
As主要成环状分布,有两块富集程度相当高的区域，污染源在点（1900,2200）（4700,4900）附近。
5.4模型的评价及改进
5.4.1模型的优缺点
（1）本模型在处理数据时使用了二维插值，处理数据简单明了，具有可信性，实用性。
（2）考虑到各重金属元素分布之间的联系，并且与各功能区结合，具有一定的现实意义。
（3）使用变异函数，能够较好地描述区域化变量的空间分布结构性与随机性特征。
（4）本模型的主要缺点就是，在插值运算时对误差考虑地不详细，误差分析不足。
5.4.2模型的改进
在研究重金属元素空间分布过程中，数据采集没考虑到时间的影响，数据的可信性较低，对数据应进行不间断的周期性重复测量；在研究传播特征时，缺少对时间因素的考虑。
从我国的看各种重金属污染主要是由近30年的人类活动所引起。
在改革开放后，工业和农业飞速发展，城市人口增多，工业三废和生活废物、废水急剧增加，环境污染也越来越严重。
综合现有的研究与前文的分析得知研究区土壤污染的主要来源是这些工业、城镇与交通污染排放物，因而得知土壤污染与工业、城镇等污染排放的增长趋势一致，将有工业活动的人类活动对土壤的影响追溯至20世纪80年代。
这样，就可以把土壤零污染年设为1981年。
截止2011年，土壤中重金属污染积累历史约为30年。
因此我们假定30年前为研究区的“零污染”时间，还没有这些污染源时，研究区各处土壤中各种重金属含量等于该种重金属在当地背景值的平均值，且假定30年来土壤污染以加速度累积增长，从而计算得到壤重金属积累的速度。
当某点处的含量值低于或等于背景值时，我们把它看成30年来都未受人为污染的点，假设其即不净化也不积累。
从而可以计算出30年来土壤污染的加速度与任意一年的污染速度。
[7]
（1）单位质量土壤重金属积累
式中，S为研究区30年以来单位质量土壤重金属积累量；
为2011年（采样时问）土壤中该种重金属污染物的含量，
为该种土壤重金属在该区域的背景值。
（2）土壤重金属积累加速度
式中，A为研究区30年以来的土壤重金属积累加速度；T为积累时间，即30年。
（3）土壤重金属积累速度
式中，V为2011年土壤重金属积累速度。
根据前面的土壤重金属累积预测公式推导出了乐观情景、无突变情景和悲观情景下时间t对各处土壤中各种重金属污染物含量的计算公式：
（1）乐观情景下土壤重金属累积预测公式
（2）无突变情景卜土壤重金属累积预测公式
（3）悲观情景下土壤重金属累积预测公式
[8]
以此模型充实原模型，收集各时间段内的数据，使模型同时可以充分体现出时间与空间关系。
六、参考文献
[1]张丛，环境评价教程[M]，北京：
中国环境科学出版社，2002。
[2]王雄军等，基于因子分析法研究太原市土壤重金属污染的主要来源[J]，生态环境，200817
（2），671-676。
[3]顾继光，周启星，王新，土壤重金属污染的治理途径及其研究进展[J]，应用基础与工程科学学报，2003，1l
（2）：
143～151。
[4]黄镇国，张伟强，珠江河口磨刀门的整治与地貌演变[J]，地理与地理信息科学，2005，2l（6）：
61～73。
[5]储国正，铜陵狮子山铜金矿田成矿系统及其找矿意义，中国地质大学（北京）博士学位论文，2003。
[6]祁志宏，土壤重金属污染的空间分析及污染评价，合肥工业大学（合肥）硕士学位论文，2005。
[7]侯景儒，黄竞先,地质统计学的理论与方法,北京：
地质出版社，1990。
[8]林艳，统计学与GIS的土壤重金属污染评价与预测，中南大学硕士学位论文，2009。
主要程序段
（1）
clc,clear
A=xlsread（'cumcm2011A附件_数据.xls','¸'附件1','A4:
E322'）;
C=xlsread（'cumcm2011A附件_数据.xls','附件2','A4:
I322'）;
x=A（:
2）;y=A（:
3）;z=A（:
4）;qy=A（:
5）;
As=C（:
2）;Cd=C（:
3）;Cr=C（:
4）;Cu=C（:
5）;
Hg=C（:
6）;Ni=C（:
7）;Pb=C（:
8）;Zn=C（:
9）;
[X,Y,Z]=griddata（x,y,z,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,Z）
title（'采样地形图'）;
[X,Y,AS]=griddata（x,y,As,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,AS）
title（'As浓度随平面地形分布图'）;
[X,Y,CD]=griddata（x,y,Cd,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,CD）
title（'Cd浓度随平面地形分布图'）;
[X,Y,CR]=griddata（x,y,Cr,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,CR）
title（'Cr浓度随平面地形分布图'）;
[X,Y,CU]=griddata（x,y,Cu,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,CU）
title（'Cu浓度随平面地形分布图'）;
[X,Y,NI]=griddata（x,y,Ni,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,NI）
title（'Ni浓度随平面地形分布图'）;
[X,Y,HG]=griddata（x,y,Hg,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,HG）
title（'Hg浓度随平面地形分布图'）;
[X,Y,PB]=griddata（x,y,Pb,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,PB）
title（'Pb浓度随平面地形分布图'）;
[X,Y,ZN]=griddata（x,y,Zn,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,ZN）
title（'Zn浓度随平面地形分布图'）;
[X,Y,QY]=griddata（x,y,qy,linspace（0,28654）',linspace（0,18449）,'v4'）;
figure,contourf（X,Y,QY）
title（'区域分布图'）;
（2）
functionR=byhs（w,v）;
R
（1）=0;
fori=1:
318
N=fix（log（abs（v（i）-v（i+1））））
forj=1:
N-1
ZZ=w（j）-w（j+1）;
pf=ZZ*ZZ;
s=sum（pf）;
end
r=s/（2*N）;
R（i+1）=r;
end
clc,clear%将附件2中重金属元素浓度粘贴到紧邻附件1的右侧
B=xlsread（'cumcm2011A附件_数据','附件1','A4:
M322'）;
A=sortrows（B,2）;
x=A（:
2）;y=A（:
3）;
As=A（:
6）;Cd=A（:
7）;Cr=A（:
8）;Cu=A（:
9）;
Hg=A（:
10）;Ni=A（:
11）;Pb=A（:
12）;Zn=A（:
13）;
Asx=byhs（As,x）;
Cdx=byhs（Cd,x）;
Crx=byhs（Cr,x）;
Cux=byhs（Cu,x）;
Hgx=byhs（Hg,x）;
Nix=byhs（Ni,x）;
Pbx=byhs（Pb,x）;
Znx=byhs（Zn,
升级会员 