2011年国赛高教杯数学建模
A题 城市表层土壤重金属污染分析
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
整体求解过程概述(摘要)
本文通过对附件中所给海量数据进行处理分析，找出主要污染原因及污染物传播特性，进而来研究人类活动影响下城市地质环境的变化模式，对地质环境的保护起到积极作用。
对于问题一，先用matlab软件对数据插值拟合，作出各元素在该城区的空间分布图。从图中可直接地观察各元素对不同地点的污染程度，再将附表中所给数据标准化。本文采用目前学术界应用最普遍的单因子浸染指数法和内梅罗综合浸染数法对不同功能区进行逐一分析。得出各个区域重金属污染程度，其中生活区、工业区与主干道路区-重度污染、山区-轻度污染、公园绿地区-中度污染，且它们之间综合污染指数排序为工业区>主干道路区>生活区>公园绿地区>山区。
对于问题二，根据问题一中计算出的单因子浸染指数，并作出各功能区不同元素的单因子浸染指数柱状图。通过比较分析发现，某些元素在各功能区单因子指数相近（如As、Ni、Cr），说明同种元素在整个城区的污染原因相同。通过matlab软件分别算出在整个城区和各功能区各元素之间的相关度，可以看出有些元素呈现一定相关性，若相关显著可知其间污染原因可能相同。为了进一步明确不同功能区土壤中重金属的污染来源，再以公园绿地区为例利用主成分分析方法，得出公园绿地区八个变量的全部信息可由四个主成分表示，并分析污染物的来源。根据上述模型进行分析，可知As的主要污染原因为：在某些有色金属的开发和冶炼中，砷化物的广泛利用，煤的燃烧。Ni和Cr相关度较高，主要污染原因有：地球化学成因影响；整个城区的多种人为活动产生的综合污染。除这些元素外，再对其他元素的污染进行分区考虑，根据其它元素之间的相关度可知：生活区Cu、Cd、Pb主要污染原因为生活废弃物；Hg、Zn和其他元素均不相关，说明其主要污染原因是重金属的运移。工业区污染原因主要为工业生产污染。山区轻度污染，主要原因为重金属的运移。主干道路区污染比较复杂，其主要污染原因为交通运输、道路建设、人为活动、金属运移等。公园绿地区处于中度污染，主要污染原因为人类活动和重金属的运移。
对于问题三，综合考虑对流、弥散、扩散吸附和微生物降解等作用, 根据质量守恒原理, 建立了污染物在饱和及非饱和土壤中运移的对流扩散数学模型:
a_1 (∂^2 C)/(∂x^2 )+a_2 (∂^2 C)/(∂y^2 )+a_3 (∂^2 C)/(∂z^2 )+a_4 ∂C/∂x+a_5 ∂C/∂y+a_6 ∂C/∂z+a_7 x+a_8 y+a_9 z+a_10+a_11 C=∂C/∂t根据附表数据拟合出扩散系数a_i（i =1…10），再利用初始浓度及长时间的稳定态作为初值条件构成微分方程的定解问题，最后采用数值模拟方法预测污染物的迁移过程、迁移范围及浓度分布特征。并由此来分析污染物的传播特征并确定污染源的位置（见表5），结合问题二中的主成分分析的结果，确定最终污染源有7个。
对问题四，问题三模型建立后，发现如果能够收集到城市详细的卫星图、城市历史地质资料、同一位置不同时间的重金属浓度、城市土壤更多的测定值等信息，能够更好的建立污染物迁移速度的演变模式。可以建立通过聚类分析法、人工神经网络的SOFM网络得出城市不同时期的综合评价结果，从而分析出城市地质环境的演变模式。
模型假设：
1、假设附件所给数据准确无误；
2、假设污染源排放污染物稳定、连续；
3、假设污染物在自然环境下自由传播；
4、假设污染物在传播过程中质量守恒；
5、假设污染物传播只与空间有关，与时间无关；
问题分析：
本题主要是研究表层土壤重金属污染的问题，主要目的是通过对附表中所给八种主要重金属的采样数据进行分析，从而得到这八种元素的空间分布，并确定各功能区的污染程度；再根据数据分析，先找出各金属间的相关度，确定其是否来自同一污染源，再根据相关分类标准，找出重金属污染的主要原因；然后根据数据建立适当函数，找出污染物传播特征，并建立适当模型，最终确定污染源的物质。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：
%syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
a=[-0.000109287652109 0.000635260069389 0.516315789473684 0.025660770031218 -0.021669595782074 12.330000000000000 8017 7210 39 1
-0.000757120898578 -0.000024470387428 -0.666666666666666 -0.189154518950437 -1.158571428571430 -16.220000000000000 4948 7293 6 1
-0.000096945211774 0.000090358584587 -0.015898771367521 -0.015780845556514 -0.056978193146417 -0.351730769230769 13855 3345 79 1
-0.000457307232306 0.004842383083857 1.479861111111110 -0.072228260869565 -0.055145228215768 -53.160000000000000 2383 3692 7 1
-0.000308833734969 0.000211980435812 4.667500000000000 0.193295454545455 -0.060533807829182 -25.530000000000000 8622 10638 4 1
0.000062607723798 -0.000109998421723 0.889880952380952 0.008987435328899 -0.012854122621565 1.737142857142860 12153 12336 16 1
-0.000052640377967 0.000012890933239 -0.008995925263367 0.000516320474777 -0.001855010660981 0.008055555555556 26015 12078 57 1
-0.005422959353530 0.000760551532441 5.000000000000010 -1.230736842105260 0.628602150537634 116.920000000000000 4777 4897 8 1
0.002510316980250 0.038800479841974 0.472566280566281 -0.054024691358025 -0.158550724637681 -1.182702702702700 3518 2571 59 1
0.000042033121506 -0.000011309266827 0.046103040086091 -0.059379157427938 0.024301270417423 -0.425079365079365 18413 11721 88 1
];
%X=[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10]';
b=[124.24
130.36
54.98
381.64
94.93
53.42
28.30
472.48
128.60
38.47
];
X=a\b
function readdata()
data=xlsread('data','sheet1','B4:E322');
save data data %坐标信息
data1=xlsread('data','sheet2','B4:I322');
save data1 data1 %浓度信息
beij=xlsread('data','sheet3','B4:C11');
save beij beij %背景信息
%************************************************
%程序功能：画Cu污染分布等值线图
%建立日期：2011/9/9
%程序员：Hec1990
%************************************************
function drawcu()
clc,clear
load data
load data1
Cu=data1(:,4);
X=data(:,1);
Y=data(:,2);
Cu_XY=[X Y Cu];
[x y]=meshgrid(0:100:3e4,0:100:2e4);
z=griddata(Cu_XY(:,1),Cu_XY(:,2),Cu_XY(:,3),x,y,'v4');
mesh(x,y,z)
%ezsurfc(z,x,y)
figure
[c h]=contour(x,y,z);
clabel(c,h)
hold on
[c1 d1]=find(data(:,4)==1);
x1=data(c1,1);
y1=data(c1,2);
plot(x1,y1,'r*');
[c2 d2]=find(data(:,4)==2);
x2=data(c2,1);
y2=data(c2,2);
hold on
plot(x2,y2,'k*')
[c3 d3]=find(data(:,4)==3);
x3=data(c3,1);
y3=data(c3,2);
plot(x3,y3,'c*')
[c4 d4]=find(data(:,4)==4);
x4=data(c4,1);
y4=data(c4,2);
plot(x4,y4,'y*')
[c5 d5]=find(data(:,4)==5);
x5=data(c5,1);
y5=data(c5,2);
plot(x5,y5,'g*')
grid on
legend('等值线','生活区','工业区','山区','交通区','绿地区');
cu=reshape(z,301*201,1);
a=reshape(x,301*201,1);
b=reshape(y,301*201,1);
save cu a b cu %训练神经网络数据集
save dist c1 c2 c3 c4 c5 %不同区域
%************************************************
%************************************************
function comp()
load data1
load beij
for i=1:8
a(:,i)=(data1(:,i)-beij(i,1))/beij(i,2);
end
As=a(:,1);Cd=a(:,2);Cr=a(:,3);Cu=a(:,4);Hg=a(:,5);Ni=a(:,6);Pb=a(:,7);Zn=a(:,8)
;
load dist
As1=As(c1,1);Cd1=Cd(c1,1);Cr1=Cr(c1,1);Cu1=Cu(c1,1);Hg1=Hg(c1,1);Ni1=Ni(c1,1);P
b1=Pb(c1,1);Zn1=Zn(c1,1);
As2=As(c2,1);Cd2=Cd(c2,1);Cr2=Cr(c2,1);Cu2=Cu(c2,1);Hg2=Hg(c2,1);Ni2=Ni(c2,1);P
b2=Pb(c2,1);Zn2=Zn(c2,1);
As3=As(c3,1);Cd3=Cd(c3,1);Cr3=Cr(c3,1);Cu3=Cu(c3,1);Hg3=Hg(c3,1);Ni3=Ni(c3,1);P
b3=Pb(c3,1);Zn3=Zn(c3,1);
As4=As(c4,1);Cd4=Cd(c4,1);Cr4=Cr(c4,1);Cu4=Cu(c4,1);Hg4=Hg(c4,1);Ni4=Ni(c4,1);P
b4=Pb(c4,1);Zn4=Zn(c4,1);
As5=As(c5,1);Cd5=Cd(c5,1);Cr5=Cr(c5,1);Cu5=Cu(c5,1);Hg5=Hg(c5,1);Ni5=Ni(c5,1);P
b5=Pb(c5,1);Zn5=Zn(c5,1);
aa=[10 30 2 5 40 5 5 1];
w=[aa(1)/sum(aa) aa(2)/sum(aa) aa(3)/sum(aa) aa(4)/sum(aa) aa(5)/sum(aa)
aa(6)/sum(aa) aa(7)/sum(aa) aa(8)/sum(aa)];
y1=[mean(As1) mean(Cd1) mean(Cr1) mean(Cu1) mean(Hg1) mean(Ni1) mean(Pb1)
mean(Zn1)];
y2=[mean(As2) mean(Cd2) mean(Cr2) mean(Cu2) mean(Hg2) mean(Ni2) mean(Pb2)
mean(Zn2)];
y3=[mean(As3) mean(Cd3) mean(Cr3) mean(Cu3) mean(Hg3) mean(Ni3) mean(Pb3)
mean(Zn3)];
y4=[mean(As4) mean(Cd4) mean(Cr4) mean(Cu4) mean(Hg4) mean(Ni4) mean(Pb4)
mean(Zn4)];
y5=[mean(As5) mean(Cd5) mean(Cr5) mean(Cu5) mean(Hg5) mean(Ni5) mean(Pb5)
mean(Zn5)];
qu1=sum(w.*y1);qu2=sum(w.*y2);qu3=sum(w.*y3);qu4=sum(w.*y4);qu5=sum(w.*y5);
