2011年数学建模A题 城市表层土壤重金属污染分析(附所有图件的Matlab源代码)

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

摘要

通过建立地区三维地形图与各种元素丰度分布等值线图，结合地形、元素传播特征、风向，利用数学方法综合确定污染源的位置，并对各类污染作简要的分析。

对于问题一，利用Matlab，做出地形图与各元素等值线图，确定各分布。在通过几何平局法，确定总体重金属污染情况。

对于问题二，基于问题一，可知汽车尾气、工业“三废”、生活垃圾共同造成重金属污染。

对于问题三，由于各元素传播特征不同，传播能力弱的，其元素丰度高处就为污染源。例如Cr Cd。传播能力强者，由有相似分布者，可联合分析，认为其为同一污染源排出。如As、Hg、Zn。

关键词：污染物分布 Matlab 等值线

一. 问题重述

通过GPS记录了某城市城区的空间坐标，给出了不同功能区土壤表层0—10cm的8种主要重金属的含量，并给出了此城区的不同重金属的背景值，均值，标准偏差。

现通过数学建模来完成以下任务：

(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。

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

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

(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？

二 .问题分析

（1）问题一属于统计归纳的数学问题。通过数学方法进行图像拟合，从中发现规律，找到图中污染集中分布区域。之后，再把八种元素综合考虑，试图确定不同地区的污染程度。（2）对于问题二，可在第一问的基础上，结合各种重金属污染的传播特点，综合确定污染原因。

（3）问题三，对于不易传播的元素，可认为其元素丰度最高值处为污染源，对于其他元素，应通过其污染物散落范围综合求解得出。

（4）问题四，可在获知风向、水流方向的因素的情况下，综合利用微分方程求解污染源。

三、模型假设

（1）假设As、Hg、Zn由同一污染源排放。

（2）对于第一问，忽略高程的影响。

（3）假设重金属传播主要靠风和水流。

（4）假设风向和水流方向不随时间变化。

四、模型的建立与求解

1、问题一

对于第一小问，首先引入化学元素的丰度的概念。

元素丰度=所测元素浓度值/浓度背景值

分别求出八种元素的丰度，用Matlab绘制元素丰度等值线，并与功能区散点图结合。

上图为地形等高线与功能区结合图

还可利用Matlab绘制次地区三维地形图

由地形图可知，调查区东北角地形较高，西南地形较低。中部地区有山谷。但总体地形坡度不大（图中三个坐标的比例尺不同）。由于所给数据大都沿西南—东北走向分布，因此图中西北、东南脚地势低是由于kriging差值外插所致。

下一页将给出八种元素丰度等值线

As元素丰度分布情况

Cd元素丰度分布情况

Cr元素丰度分布情况

Cu元素丰度分布情况

Ni元素丰度分布情况

Hg元素丰度分布情况

Pb元素丰度分布情况

Zn元素丰度分布情况

对八种元素的丰度分布情况进行分析，可知：

As元素主要分布在城区的西部工业区、南部。

Cd元素主要分布在城区西部，尤其是西南角工业区。此外，位于（20000，,12000）的山谷交通区也有集中分布。

Cr元素集中分布于西南角。污染源极有可能在该点。

Cu元素集中分布于西南角（8号测点附近），其余地区基本无污染。

Ni元素主要在22测点（交通区）超标严重。

Hg分布于三个地区。西南角的工业区附近。城市南部的工业区附近。城市中部的山谷。

Pb主要分布在西南角，为工业区和生活区。

Zn分布有五个峰值，主要为西南角工业区，中部工业区以及山谷。

对于第二问，可将八种元素的污染影响平均化，即将八种元素在个点丰度值取均值，的平均丰度。

用Matlab作图。

分析图件，可得出以下结论：

①该城区的西南角落（8号测点附近）重金属污染较多，污染严重。

②该城区的南部（185号测点附近）的重金属污染也较严重。

③城区中部地势低洼的山谷部位，由于其四周有工业区，污染也较严重

④城市中心（大面积的交通区），污染并不大。

2、问题二

由第一问可知，该城市主要有三个重金属污染集中区域。分别在图中（2000，4000），(14000，2000)，(16000，90000)附近。前者由于靠近工厂，排放大量的废水、废气、废渣，且其地形处于下坡，上坡带还有居民区，城市垃圾也会顺坡而下。后两者属于山谷地区，相对于周边较开阔，风速和水流速度会下降，沉积物较多，易造成污染。

城市中土壤的重金属含量明显高于山区。城市中的工业区、交通线两侧污染严重，绿地污染较少。

公路两侧的污染以Pb、Zn、Cd、Cr、Cu为主，主要来自于汽油的燃烧，轮胎磨损产生的含锌粉尘。

Hg可以由于被土壤中的粘土矿物和有机物吸附，而进入土壤表层，产生污染。

总之，工业化程度越高，污染时间越长，污染越严重。

3、问题三

对于Cu元素，分布极为集中，Cu元素的污染源就在附近。即大约在（2500,4000）的工业区中。

Cd元素不易在土壤中传播，主要集中在污染源附近。主要为道路两侧，由汽车引起，其次为西南角的工业区，由于工厂排放产生。还有一个交通区和居民区所夹的山脚地带，Cd含量也很高，说明Cd也来源于生活垃圾。

Cr也不易在土壤传播，主要污染源在西南角的工业区。

Pb分布于西南角工业区，为污染源，且处于低地势区。另外，由于汽车尾气排放，交通感到两侧也有分布。

Ni只在西南角的一交通区（22测点）集中分布，可能由于公路引起。也可能是由于此处地势低，其西北方向的工厂排放的Ni汇集而成。

As、Hg、Zn的分布比较类似。可能由同一污染源排放。

经过对此三种元素丰度的加权计算，公式为：

加权丰度=（As丰度/As在此地区的平均丰度+Hg丰度/Hg在此地区的平均丰度+Zn丰度/Zn 在此地区的平均丰度）/3

用Matlab作图

通过二位等值线图与三维等值线图，

As、Hg、Zn综合图

As、Hg、Zn元素丰度分布三维图

可得三个元素丰度峰值点。即（2383,3994），（14000,2200），（15500,9200）

假设它们由同一污染源排放。确定污染源位置。

用三点过圆心法求得即求

(x-2383)^2+(y-3994)^2=r^2

(x-14000)^2+(y-2200)^2=r^2

(x-15500)^2+(y-9200)^2=r^2

此三元二次方程组的解，用Matlab解得

x=7068.0

y=8790.7

将此点上图

由图可知，污染源位于西北部工业区。

4、问题四

对于第一问，所建模型对风向、水流、地形的考虑不够，不能有效的确定污染物传播路径，因而不好确定污染源。

如果知道此地的主导风向和水流情况，再知道此地的地质情况，就可以了解元素迁移路径以及及此地微生物分解情况，从而建立更加精确的模型，进行求解

对于第二问，知道上述情况之后，通过流体力学知识，建立微分方程，计算多条传播路径，对污染源进行精确定位。

五、模型的评价与推广

1、模型的优点

①通过等值线图来分析元素富集情况，简单直观。

②忽略了一些复杂的实际情况，易于分析

2、模型的缺点

①在确定污染源时对风向、水流考虑不足。

②求元素在城区空间分布时，忽略了地形影响。

3、模型的推广

可推广至评价其他环境污染问题。如氯气泄露，求其影响范围等等。

附录

Matlab源码

求三元二次方程组

[x,y,r]=solve('(x-2383)^2+(y-3994)^2=r^2','(x-14000)^2+(y-2200)^2=r^2','(x-15500)^2+(y-9 200)^2=r^2')

绘图

clear

close all

%三维地形图

D=xlsread('PLOT1.xls');

M=xlsread('PLOT2.xls');

x=D(:,1);

y=D(:,2);z=D(:,3);c=D(:,4);

xi=linspace(min(x),max(x));

yi=linspace(min(y),max(y));

Zi=griddata(x,y,z,xi,yi','v4');

h=surf(xi,yi,Zi);

set(h,'EdgeAlpha',0.3)

colormap;

axis tight

xlabel('X');

ylabel('Y');

zlabel('Z');

title('地形图')

colorbar

mat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn','综合','As Hg Zn 综合'};

marker={'*','o','s','^','p'};

color={'k','r','m','c','b'};

str={'等高线','生活区','工业区','山林区','交通区','绿地区'};

str1={'等值线','生活区','工业区','山林区','交通区','绿地区'};

%地形等高线图

figure

contour(xi,yi,Zi,-150:20:300,'LineWidth',1.5);

xlabel('X');

ylabel('Y');

title(['地形等高线'])

axis tight

axis image

colorbar

grid on

hold on

for i=1:5

loc=c==i;

plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i}) end

legend(str,'location','SouthEast')

%各元素丰度等值线图

for i=1:10

Zi=griddata(x,y,M(:,i),xi,yi','v4');

figure

contour(xi,yi,Zi,'LineWidth',1.5);

xlabel('X');

ylabel('Y');

title([mat{i},'元素丰度等值线'])

axis tight

axis image

colorbar

grid on

hold on

for i=1:5

loc=c==i;

plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i});

end

legend(str1,'location','SouthEast')

end

%As Hg Zn 三维分布图

figure

contour3(xi,yi,Zi,10)

h=surface(xi,yi,Zi,'EdgeColor',[.8 .8 .8],'FaceColor','none');

set(h,'EdgeAlpha',0.5)

grid off

xlabel('X');

ylabel('Y');

zlabel('Z');

%As Hg Zn 污染源确定图

figure

contour(xi,yi,Zi,'LineWidth',1.5);

xlabel('X');

ylabel('Y');

title(['As Hg Zn 污染源确定图'])

axis tight

axis image

colorbar

grid on

hold on

for i=1:5

loc=c==i;

plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i});

end

plot(7068,8791,'p','markerfacecolor','g','MarkerEdgeColor','g');

legend([str1,'污染源'],'location','SouthEast')

绘图所需数据表格