数学建模A题城市表层土壤重金属污染分析附所有图件的Matlab源代码.docx
《数学建模A题城市表层土壤重金属污染分析附所有图件的Matlab源代码.docx》由会员分享，可在线阅读，更多相关《数学建模A题城市表层土壤重金属污染分析附所有图件的Matlab源代码.docx（34页珍藏版）》请在冰点文库上搜索。
数学建模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作图。
分析图件，可得出以下结论：
1该城区的西南角落（8号测点附近）重金属污染较多，污染严重。
2该城区的南部（185号测点附近）的重金属污染也较严重。
3城区中部地势低洼的山谷部位，由于其四周有工业区，污染也较严重
4城市中心（大面积的交通区），污染并不大。
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=
y=
将此点上图
由图可知，污染源位于西北部工业区。
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'）
绘图
clear
closeall
%三维地形图
D=xlsread（''）;
M=xlsread（''）;
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',
colormap;
axistight
xlabel（'X'）;
ylabel（'Y'）;
zlabel（'Z'）;
title（'地形图'）
colorbar
mat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn','综合','AsHgZn综合'};
marker={'*','o','s','^','p'};
color={'k','r','m','c','b'};
str={'等高线','生活区','工业区','山林区','交通区','绿地区'};
str1={'等值线','生活区','工业区','山林区','交通区','绿地区'};
%地形等高线图
figure
contour（xi,yi,Zi,-150:
20:
300,'LineWidth',;
xlabel（'X'）;
ylabel（'Y'）;
title（['地形等高线']）
axistight
axisimage
colorbar
gridon
holdon
fori=1:
5
loc=c==i;
plot（x（loc）,y（loc）,marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i}）
end
legend（str,'location','SouthEast'）
%各元素丰度等值线图
fori=1:
10
Zi=griddata（x,y,M（:
i）,xi,yi','v4'）;
figure
contour（xi,yi,Zi,'LineWidth',;
xlabel（'X'）;
ylabel（'Y'）;
title（[mat{i},'元素丰度等值线']）
axistight
axisimage
colorbar
gridon
holdon
fori=1:
5
loc=c==i;
plot（x（loc）,y（loc）,marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i}）;
end
legend（str1,'location','SouthEast'）
end
%AsHgZn三维分布图
figure
contour3（xi,yi,Zi,10）
h=surface（xi,yi,Zi,'EdgeColor',[.8.8.8],'FaceColor','none'）;
set（h,'EdgeAlpha',
gridoff
xlabel（'X'）;
ylabel（'Y'）;
zlabel（'Z'）;
%AsHgZn污染源确定图
figure
contour（xi,yi,Zi,'LineWidth',;
xlabel（'X'）;
ylabel（'Y'）;
title（['AsHgZn污染源确定图']）
axistight
axisimage
colorbar
gridon
holdon
fori=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'）
绘图所需数据表格
PLOT1
x（m）
y（m）
海拔（m）
功能区
74
781
5
4
1373
731
11
4
1321
1791
28
4
0
1787
4
2
1049
2127
12
4
1647
2728
6
2
2883
3617
15
4
2383
3692
7
2
2708
2295
22
4
2933
1767
7
4
4233
895
6
5
4043
1895
14
1
2427
3971
2
1
3526
4357
7
4
5062
4339
5
4
4777
4897
8
1
5868
4904
16
4
6534
5641
6
1
5481
6004
0
4
4592
4603
6
1
2486
5999
2
1
3299
6018
4
4
3573
6213
5
1
4741
6434
5
5
5375
8643
15
1
5635
7965
29
4
5394
8631
12
4
5291
7349
10
4
4742
7293
9
2
4948
7293
6
2
5567
6782
7
2
7004
6226
11
4
7304
5230
10
1
7048
4600
24
4
8180
4496
15
4
9328
4311
24
1
9090
5365
20
4
8049
5439
18
4
8077
6401
29
1
8017
7210
39
4
6869
7286
18
4
7056
8348
37
1
7747
8260
49
4
8457
8991
21
4
9460
8311
45
4
9062
7639
45
4
9319
6799
49
4
10631
6472
57
4
10685
5528
34
4
10643
4472
45
4
11702
4480
71
4
11730
5532
54
4
11482
6354
61
4
10700
8184
50
4
10630
8774
29
4
11678
8618
17
4
11902
7709
30
4
13244
7056
37
4
12746
8450
21
4
12855
8945
18
4
13797
9621
18
4
14325
8666
23
3
15467
8658
17
1
12442
4329
65
4
13093
4339
56
4
13920
5354
79
4
14844
5519
62
1
16569
6055
78
1
16387
6609
44
5
16061
7352
28
5
15658
7594
24
4
14298
7418
36
1
14177
6684
35
4
15092
6936
32
5
12778
5799
93
4
17044
10691
93
3
17087
11933
43
4
17075
12924
25
4
17962
12823
25
4
18413
11721
88
3
19007
11488
84
3
18738
10921
53
3
17814
10707
64
4
18134
10046
41
4
17198
9810
37
4
17144
9081
20
4
18393
9183
26
4
19767
8810
46
4
21006
8819
55
4
21091
9482
43
4
22846
9149
69
4
23664
9790
46
3
22304
10527
40
4
21418
10721
35
1
21439
11383
45
4
20554
11228
43
4
20101
10774
40
4
21072
10404
32
4
20215
9951
31
4
18993
12371
78
3
19968
12961
42
3
21766
12348
67
4
22674
12173
52
2
22535
11293
54
3
25221
5795
27
4
26453
5577
11
1
26416
6508
14
1
27816
5581
11
3
25361
6423
49
3
24065
7353
104
3
25998
7032
51
3
27177
7771
17
3
26424
8639
8
4
26073
8807
38
3
24631
9422
76
3
24702
9522
62
3
25461
9834
68
3
24813
10799
46
4
26086
11094
53
3
26015
12078
57
3
27700
11609
165
3
27696
11621
169
3
27346
13331
100
3
26591
13715
126
3
27823
14737
189
3
27232
14482
150
3
24580
13319
107
3
24153
12450
71
3
22965
13535
78
3
23198
13523
62
4
24685
14278
98
3
28654
8755
23
3
24003
15286
90
3
21684
13101
114
3
22193
12185
79
3
17079
5894
81
3
15255
5110
110
3
15007
5535
70
3
3518
2571
59
5
3469
2308
52
5
3762
2170
30
5
3927
2110
27
5
4153
2299
73
5
3267
793
0
5
4684
1364
37
5
5495
1205
9
5
5664
1653
13
5
5541
2093
26
5
5451
2757
92
5
4020
2990
27
4
4026
3913
13
4
5101
4080
13
1
5438
3994
10
2
5382
3012
50
1
5314
2060
40
4
5503
1127
6
1
5636
133
17
1
6605
374
6
1
7093
1381
45
4
7100
2449
89
4
6837
3490
28
4
7906
3978
22
4
8045
3052
39
4
8394
2035
27
4
8403
1075
6
4
8079
0
16
4
9663
1288
3
4
9469
2286
15
4
9178
3299
42
4
9095
3975
26
4
10225
3821
19
4
10210
2789
19
4
10340
1764
7
4
11557
1581
7
4
11415
2585
12
4
11649
3515
27
1
12734
4015
43
2
12696
3024
27
4
12400
2060
13
4
12591
1063
18
1
13765
1353
15
4
13694
2357
33
4
13855
3345
79
1
14862
2524
28
1
14896
1603
4
2
15387
729
8
1
15810
2307
8
1
16032
3061
35
4
15801
3966
115
3
15087
3512
69
5
16872
2798
10
4
17734
3629
14
4
16823
4207
67
5
17008
4775
82
3
17203
6218
40
1
17005
7212
33
4
16947
7487
41
2
16301
8299
24
1
17904
8287
25
1
18303
7385
39
5
18438
6539
22
4
18556
5588
15
5
18954
4874
4
4
18012
4414
20
4
19072
8519
36
4
20282
8590
57
4
21475
8540
85
3
21450
7555
58
4
20261
7586
29
3
19569
7348
70
3
19411
6934
28
3
19501
6091
9
4
20582
6548
13
5
19909
5300
3
4
21018
5764
9
4
22176
5492
25
4
23359
5325
28
3
23238
6502
169
3
22624
4818
27
3
21703
6591
65
3
5006
8846
6
2
5734
9659
3
4
6395
10443
4
2
7405
10981
6
2
8446
11200
4
2
7612
11938
2
2
7912
12840
1
4
8866
13143
3
2
9296
13102
9
4
9475
12000
9
2
9212
11305
5
2
8629
12086
1
2
7776
10613
9
2
8622
10638
4
4
9237
9872
28
4
8307
9726
14
4
7106
9467
44
2
6423
8831
40
2
7458
8920
36
2
8904
8868
24
4
10547
9591
32
4
10398
10360
0
4
10395
11203
8
1
11529
11243
16
4
11563
10298
12
4
11646
9381
14
升级会员 