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

A题城市表层土壤重金属污染分析摘要通过建立地区三维地形图与各种元素丰度分布等值线图，结合地形、元素传播特征、风向，利用数学方法综合确定污染源的位置，并对各类污染作简要的分析。对于问题一，利用Matlab，做出地形图与各元素等值线图，确定各分布。在通过几何平局法，确定总体重金属污染情况。对于问题二，基于问题一，可知汽车尾气、工业“三废”、生活垃圾共同造成重金属污染。对于问题三，由于各元素传播特征不同，传播能力弱的，其元素丰度高处就为污染源。例如CrCd。传播能力强者，由有相似分布者，可联合分析，认为其为同一污染源排出。如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.0y=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-9200)^2=r^2')绘图clearcloseall%三维地形图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;axistightxlabel('X');ylabel('Y');zlabel('Z');title('地形图')colorbarmat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn','综合','AsHgZn综合'};marker={'*','o','s','^','p'};color={'k','r','m','c','b'};str={'等高线','生活区','工业区','山林区','交通区','绿地区'};str1={'等值线','生活区','工业区','山林区','交通区','绿地区'};%地形等高线图figurecontour(xi,yi,Zi,-150:20:300,'LineWidth',1.5);xlabel('X');ylabel('Y');title(['地形等高线'])axistightaxisimagecolorbargridonholdonfori=1:5loc=c==i;plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i})endlegend(str,'location','SouthEast')%各元素丰度等值线图fori=1:10Zi=griddata(x,y,M(:,i),xi,yi','v4');figurecontour(xi,yi,Zi,'LineWidth',1.5);xlabel('X');ylabel('Y');title([mat{i},'元素丰度等值线'])axistightaxisimagecolorbargridonholdonfori=1:5loc=c==i;plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i});endlegend(str1,'location','SouthEast')end%AsHgZn三维分布图figurecontour3(xi,yi,Zi,10)h=surface(xi,yi,Zi,'EdgeColor',[.8.8.8],'FaceColor','none');set(h,'EdgeAlpha',0.5)gridoffxlabel('X');ylabel('Y');zlabel('Z');%AsHgZn污染源确定图figurecontour(xi,yi,Zi,'LineWidth',1.5);xlabel('X');ylabel('Y');title(['AsHgZn污染源确定图'])axistightaxisimagecolorbargridonholdonfori=1:5loc=c==i;plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i});endplot(7068,8791,'p','markerfacecolor','g','MarkerEdgeColor','g');legend([str1,'污染源'],'location','SouthEast')绘图所需数据表格PLOT1x(m)y(m)海拔(m)功能区7478154137373111413211791284017874210492127124164727286228833617154238336927227082295224293317677442338956540431895141242739712135264357745062433954477748978158684904164653456416154816004044592460361248659992132996018443573621351474164345553758643151563579652945394863112452917349104474272939249487293625567678272700