数学建模论文城市表层表层土壤重金属污染分析.doc

《数学建模论文城市表层表层土壤重金属污染分析.doc》由会员分享，可在线阅读，更多相关《数学建模论文城市表层表层土壤重金属污染分析.doc（31页珍藏版）》请在三一办公上搜索。

1、城市表层表层土壤重金属污染分析摘要 本文采用离散点数据插值法拟合原始数据，确定重金属元素空间分布情况，确定污染区域；通过内梅罗综合指数法对各区域的污染程度进行评定，确定各区域污染程度。通过数据分析，明确污染主要来自主干道路区和工业区，工业活动是污染主要来源，人为活动是主要影响方式。污染传播用扩散方程来描述，通过对不同元素的扩散方程解进行拟合，并对其进行相关系数分析，据此来获得污染传播情况；建模结果表明扩散微分方程能很好的描述污染传播情况；忽略海拔对污染扩散的影响，根据扩散方程的解取最大值获得某种重金属元素污染源的位置不止一个，分别为：As(18313,10106)、（12687,3034），C

2、d（21489,11424）、（4255,2349）等。关键字：内梅罗综合指数 离散点数据插值 扩散自由定律 二维扩散方程 目录问题一 重金属元素空间分布及污染程度21.1 重金属元素空间分布与分析21.2 重金属污染程度判别71.2.1单因子污染指数法71.2.2内梅罗综合污染指数7问题二 重金属污染分析102.1 重金属污染来源与成因分析102.1.1 重金属之间相关性分析102.1.2 各功能区重金属元素圆饼图分析112.2 重金属污染程度132.2.1从整体的角度分析13问题三 重金属的传播特征和污染源位置153.1特征分析153.2 模型假设153.3 符号说明153.4传播模型建立

3、与求解163.5求污染源位置坐标173.6模型验证193.7 确定重金属污染源的位置21问题四 模型的评价与改进224.1模型的优点224.2缺点和改进224.2.1模型缺点224.2.2模型改进224.3 模型的推广23参考文献24附录25问题重述近年来，城市污染已成为人们关注的焦点。现通过对某城区进行调查，将该城区分为一公里左右的网格，对网格内的表层土壤进行采样调查，得到采样点的位置，海拔，所属功能区，及采样点重金属浓度。通过对所得样本值进行分析，处理；数据分析后建立模型，利用给出的主要重金属元素背景值，依据所得数据完成以下任务：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城

4、区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析所建立的模型的优缺点，确定还应收集的信息；依据这些信息，用建立的模型解决问题，从而对污染进行预测和研究城市地质环境的演变模式。问题一 重金属元素空间分布及污染程度1.1 重金属元素空间分布与分析由于各种重金属在各个区域的浓度都不同，所以对不同区域的影响程度也不同，根据题中所给采样点的坐标和各重金属的浓度，用MATLAB软件绘制附录一出该地区的三维图和重金属在该区的三维分布图及8种重金属的浓度为等值线的三维地形图附录二。图该地区三维图图的等

5、浓度三维分布图 图的等浓度三维分布图图的等浓度三维分布图图的等浓度三维分布图图的等浓度三维分布图图的等浓度三维分布图图的等浓度三维分布图图的等浓度三维分布图从上面的八个污染元素的浓度立体图和二维平面图形可以看出来的主要污染来源是主干道区和工业区，的主要污染来源于工业区、主干道区、生活区，的主要污染来源是山区、主干道区区及公园绿地区，的主要污染来源是生活区、工业区、山区及主干道路区，的主要污染是来源是工业区、生活区及山区，的主要污染来源是山林区和工业区，的主要污染来源于工业区和生活区，的主要污染来源于生活区和主干道区。1.2 重金属污染程度判别土壤污染评价方法采用单因子污染指数法和内梅罗综合污染

6、指数法1.2.1单因子污染指数法计算公式为： 式中， 为土壤中污染物i的环境质量指数；为污染物i的实测质量分数（mgkg-1）；为污染物i的评价标准(mgkg-1)1，一般取二类标准。 1.2.2内梅罗综合污染指数计算公式为：=式中， 为某地区的综合污染指数；为土壤污染物中污染指数最大值；为土壤污染物中污染指数平均值。经计算得出各地区单因子污染指数如下表表1 各地区单因子污染指数表地区污染指数重金属1生活区2工业区3山区4交通区5公园绿地去As0.59530.76320.54840.84660.599Cd1.05451.11280.68311.28451.043Cr1.16440.75160.

7、59491.27730.4829Cu1.22123.64420.65712.67120.9314Hg0.80193.76060.49774.05551.2059Ni0.71510.78430.9461.26540.6663Pb0.95010.93790.50070.6390.693Zn2.50221.95210.77862.82971.7571根据农田土壤环境质量监测技术规范NYT3952000中的土壤污染分级标准进行评价表2 等级评价综合污染指数污染等级污染水平P0.7安全清洁0.7P1.0警戒限尚清洁1.0P2.0轻污染土壤污染物超标2.0P3.0中污染土壤受中度污染P3.0重污染土壤受污

8、染已相当严重可知各地区各种重金属污染水平如下表表 3 重金属污染水平 功能区地区污染指数重元素金属1生活区2工业区3山区4交通区5公园绿地区As清洁尚清洁清洁尚清洁清洁Cd超标超标清洁超标超标Cr超标尚清洁清洁超标清洁Cu超标重污染清洁中污染尚清洁Hg尚清洁重污染清洁重污染超标Ni尚清洁尚清洁尚清洁超标清洁Pb尚清洁尚清洁清洁清洁清洁Zn中污染超标尚清洁中污染超标 由上表可知，除山区外，其余各区均受到重金属不同程度的污染，其中工业区污染最严重，其次为交通区.在城市中各区与自然区相比，Cd,Cr除山区外均超标，Cu和Hg在工业区和交通区严重污染，各区的污染程度为：工业区主干道生活区公园绿地山区。

9、说明该地区工业区人为污染程度最大，山区最小，但各地区均受到人为污染影响。 问题二 重金属污染分析2.1 重金属污染来源与成因分析对于问题二：则要求通过数据分析说明重金属污染的主要原因，即要分析重金属污染的主要来源和成因。在上题的结论基础上，依据各种重金属在不同区域的平均浓度大小，通过分析比较，并考虑各区域的特点，利用MATLAB做出各个功能区的重金属污染所占比例的圆饼图，如下所示 表4 各功能区的重金属平均浓度重金属功能区Cd (ng/g)Cr (g/g)Cu (g/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)Cd (ng/g)289.9669.0249.4093.0

10、418.3469.11237.01289.96393.1153.41127.54642.3619.8193.04277.93393.11152.3238.9617.3240.9615.4536.5673.29152.32360.0158.0562.22446.8217.6263.53242.85360.01280.5443.6430.19114.9915.2960.71154.24280.542.1.1 重金属之间相关性分析在进行重金属的相关性分析：用SPSS 对重金属污染数据进行因子分析，分析8 种金属污染情况的相关性，相关系数表如下表所示。表5 各重金属相关矩阵AsCdCrCuHgNiPb

11、ZnAs1.0000.2550.1890.1600.0640.3170.2900.247Cd0.2551.0000.3520.3970.2650.3290.6600.431Cr0.1890.3521.0000.5320.1030.7160.3830.424Cu0.1600.3970.5321.0000.4170.4950.5200.387Hg0.0640.2650.1030.4171.0000.1030.2980.196Ni0.3170.3290.7160.4950.1031.0000.3070.436Pb0.2900.6600.3830.5200.2980.3071.0000.494Zn0.

12、2470.4310.4240.3870.1960.4360.4941.000从表中可以看出，Cu 与Cr、Hg、Ni、Pb 呈显著正相关;Cd 与Pb 呈显著正相关，Zn 与Cd 呈显著正相关; Ni 与Cr 呈高度显著正相关，Cr 与Zn 呈显著正相关;Ni 与Zn 呈显著正相关; Pb 和Zn 呈显著正相关。土壤中Cr、Ni、Zn 表现出了较高的正载荷，说明这些金属有相似的物源，以及相似的迁移和累积特征。这些金属污染主要源于工业三废、生活污水、交通废气排放。2.1.2 各功能区重金属元素圆饼图分析根据题所给数据进行分析处理，由问题一得出来的平均值进行分析，利用MATLAB编程画出各个功能区

13、的重金属所占比例的圆饼图附录三，如下所示： 图工业区 图生活区 图山区 图公园绿地区 图交通区从圆饼图可以清楚的看出工业区和主干道区的各种金属所占比例基本一致，可以推算出工业区和主干道区重金属污染物来源相同，主要来源是工业制造和交通废气排放。生活区、山区及公园绿地区各种金属所占比例基本一致，可以推算出生活区、山区及公园绿地区重金属污染物来源相同，主要来源于生活污水和重金属废物垃圾的堆积。 则通过各种金属相关性和功能区圆饼图分析更能个较为准确判定出重金属的污染程度和污染源。2.2 重金属污染程度 根据上述两个方法的结合与对比，我们可以得出下面的污染程度关系式：工业区交通区生活区公园绿地区山区。与

14、实际情况相符。2.2.1从整体的角度分析 工业区污染最大，占整体区域的36%。分析其污染的主要来源有金属矿山的开采、冶炼、重金属尾矿、冶炼废渣和矿渣堆放等。而且其扩散和传播污染的途径也是非常广泛的。 交通区重金属的污染程度相对也较大，其污染主要体现为大气污染，来自工业生产排放的废气、汽车尾气排放生的大量含重金属的有害气体和粉尘。 生活区重金属污染居中，重金属的主要来源有生活垃圾、污水、医疗垃圾和细菌等。 公园绿地区重金属污染程度较轻，其主要来源有农业农药和化肥。但是，植物有吸收金属矿物的作用，相对减轻了重金属的危害。 山区重金属污染最轻，只占到整体区域的8%，表明不仅重金属的来源相对较少，而且

15、，重金属污染与重力和海拔也有一定关系，山区中植物和一些微生物都有净化的作用。2.2.2 从局部的角度分析 从表（3）可以看出八种重金属的主要分布区域，再根据这些区域的区别可以分析出各金属污染物的成因。其主要来源如下表所示： 图5 污染物来源污染物来源砷As矿山开采、冶炼、加工排放的三废、飘尘、农药等镉Cd电镀、塑料电池电子、染料化工工业废水铬Cr电镀、染料化工工业废水、飘尘、矿山开采、冶炼、主要来源是汽车燃料和轮胎磨损。铜Cu电镀工业废水、农药汞Hg塑料电池电子工业废水、矿山开采、冶炼、加工排放的三废镍Ni电镀工业废水、矿山开采、冶炼、加工排放的三废铅Pb电镀工业废水、矿山开采、冶炼、加工排放

16、的三废、汽车尾气锌Zn电镀工业废水问题三 重金属的传播特征和污染源位置3.1特征分析重金属的传播可认为是由一个圆形小面元的重金属污染源不断的输出重金属造成的，重金属在出圆形小面元后可认为满足自由扩散定律在三维空间自有扩散，由于海拔的变化与空间距离相比在尺度上可认为不变，所以扩散速度不会因为海拔变化而不同，在时间足够长以后重金属分布稳定。基于上述假设，任意位置重金属的浓度只取决源的重金属浓度和该点与源的距离。3.2 模型假设(1) 该城区土壤是均匀的；(2) 采样在所属功能区是随机，且均匀的，样本具有代表性；(3) 重金属在在该城市传播的时间足够长；(4) 重金属在该城市的传播趋于稳定；(5)

17、重金属的传播没有方向性，不考虑重金属元素的重力；(6) 对数据处理时先不考虑人为客观影响。(7) 忽略自然界对重金属的净化作用。3.3 符号说明 :重金属元素的浓度（As,Cd,Cr,Cu,Hg,Ni,Pb,Zn），本文中对八种重金属浓度进行单独分析，例如对重金属元素As分析时用表示其浓度，对其他金属元素浓度表示时也用表示。 :重金属元素的背景浓度（As,Cd,Cr,Cu,Hg,Ni,Pb,Zn），本文中对八种重金属浓度进行单独分析，例如对重金属元素As分析时用表示其浓度，对其他金属元素浓度表示时也用表示。 :X方向上第i点的坐标值(i=1,2,.,N)，例如对样本点进行计算时X方向上第100

18、个点用X100来表示。 : Y方向上第j点的坐标值(j=1,2,.,N)，例如对样本点进行计算时Y方向上第100个点用Y100来表示。 N : 总共采样点数。 : 点（，）的坐标表示。: 平面两点的距离。 P : 土壤质量综合指数，值可以为小数。 : 某地第i个评价项目的连续的评价分类。 U: 重金属浓度: 极角: 以源点为圆心的微元半径: 以源点为圆心的微元半径重金属浓度: 以样点为终点的向量: 以源点点为终点的向量3.4传播模型建立与求解 重金属污染源的圆形小面元在面元边界上保持浓度不变，即： 重金属在面元外自由扩散满足自由扩散方程：在面元外通解为代入边界条件傅立叶展开解得上述方程为重金属

19、传播方程，当时，此时我们确定该点为污染源。3.5求污染源位置坐标 污染源可已通过将两个采样点带入重金属传播方程求得一个可选污染源点，为了达到更高的精度，先求得多个可选源点，然后将采样点与可选污染源距离的倒数作为权重求取最佳可选污染源。当中心不在原点时传播方程变为求源点即为求出重金属传播方程，需要确定的坐标值和的值，所以取峰值点附近的数据Bi，设其中重金属浓度最大，次大和第三大的点为B0，B1，B2；将B0，B1，B2代入 由方程（1）,（2）,（3）可确定的第一源点可取值R0()第二步：将B1，B2，B3代入 由方程（4）,（5）,（6）可确定的第二源点可取值R1(),通过这种方法可以确定一系

20、列可选污染源源点Rn()（n=0，1，2，.）直到将所有数据Bi代入，然后进行加权求平均确定最佳可选源点。最佳可选源点表达为其中为了求取，先令，求均值获得均值RU即：所以求得重金属传播方程为 (即为)3.6模型验证通过上面算法对区域内的Hg求源点的=(13699.72,2354.01) 将该值代入附近的Hg浓度值和相应的位置坐标值代入得RU值如下： 表6 RU值X12734 12696 12400 12591 13765 13694 Y4015 3024 2060 1063 1353 2357 RU109496 192991 127858 98674 162514 86887 X14862 1

21、4896 15387 15810 16032 15801 Y2524 1603 729 2307 3061 3966 RU34085 26849 89037 86572 102386 113905 X1484413920124421309313855 15087 Y55195354432943393345 3512 RU1581901383731030079339485272 79535 求得重金属的传播方程为，在上表中RU标准差=31050。对和进行线性拟合得到相关系数=0.9368；由此可以证明本模型在重金属Hg在上的传播可以用重金属传播方程描述。其他点也有类似性质。图15 传播源的坐标图

22、 3.7 确定重金属污染源的位置为了求解方便，取离源最近的5个点即浓度最大的五个点对该曲面进行拟合附录四，运用上述方法求得八种重金属的污染源如下表：表7 八种重金属的的污染源位置元素AsAsCdCdCdCrCrCrx(m)18313 12687 21489 4255 21163 3279 4497 10695 y(m)10106 3034 11424 2349 9652 5917 4501 5485 元素CuCuHgHgHgNiNiNix(m)2408 3238 15206 13672 2518 3269 22065 24207 y(m)3871 5930 9151 2249 2363 588

23、1 12072 12369 元素PbPbPbZnZnZnx(m)1624 4765 14161 3160 9426 13702 y(m)2853 5065 11908 6000 4491 9667 由上表可得：一种重金属元素在该城区是分布的远不止一个，比较符合问题一中重金属元素的空间分布，说明模型建立比较合理。 问题四 模型的评价与改进4.1模型的优点(1)由于本模型是基于二维扩散方程推导出来的，所以能预测源点的位置和任意一点的浓度值。(2)本模型原理简单，假设条件通俗易懂，模型推导结果简明，易于推广。4.2缺点和改进4.2.1模型缺点(1)本模型计算步骤较多，计算较为麻烦。(2)本模型所求的

24、解为稳态解，不能描述重金属传播的动态过程。(3)本模型没有考虑水流动的作用和海拔对污染传播的影响。4.2.2模型改进由于需要获得较高计算精确度，所以计算较为复杂，所以需要计算机辅助完成繁杂的计算，为了使计算精简需对，对计算步骤进行模块化。为了描述重金属传播的动态过程，需要在模型引入外加条件的影响F(r,t)，改进的模型表达为，其中 (为土壤密度)由于模型重金属传播还受地表水和地下水流动的影响，所以还需要提供水的流动信息，为了引入水流动的影响我们引入影响函数，引入该因子后本模型变为：对于此模型我们取边界条件： 采用分离变量法求通解，然后用上述距离倒数作权重加权的方法拟合曲面，并确定拟合效果。4.

25、3 模型的推广本模型具有极大的准确性，适合于环保部门的环境监测，用于查找污染企业的位置，由于只需几个点检测即可确定污染企业位置，为环保部门减轻压力。也可以作为矿业部门的重金属矿源定位，由于重金属矿石具有重金属污染源的类似性质，重金属矿石中的重金属元素满足自由扩散定律向外扩散，满足重金属传播方程。所以通过外围点的检测，可以准确定位矿源的位置，给探矿部门提供极大方便。本模型可应用于农业生产，化肥在平坦田里的施肥效果类似与重金属自由扩散，所以可以通过本模型的拟合来研究施肥效果，由于施肥时无法保证绝对均匀，所以需要研究不均匀施肥效果，本模型可以通过源处肥料的扩散来研究周围肥料浓度。 参考文献1 杨启航

26、,数学建模,北京:高等教育出版社,2003年。2孙贤斌，重金属污染对生物群落结构及空间分布的影响N，应用生物学报，2007 年。3李建，郑春江，环境背景数据手册M，北京：中国坏境科学出版社，1989年。4 韩中庚，数学建模方法及其应用，北京：高等教育出版社，.2005。5 苏金明，阮沈勇，Matlab实用教程（第版），北京：电子工业出版社,2008 年。6谢小进，上海城郊地区城市化进程与农用土壤重金属污染的关系研究N,资源科学，2009年。7 王雄军等，基于因子分析法研究太原市土壤重金属污染的主要来源(J)，生态环境，17(2)：671-676,2008年。8 赵明登, 徐林春, 严佳庆,等.

27、 求解二维扩散方程的有限近似解法J，水动力学研究与进展, 2005, 9: 629-635。附录附录一x=A(:,1);y=A(:,2);z=A(:,3);X,Y,Z=griddata(x,y,z,linspace(min(x),max(x),200),linspace(min(y),max(y),200),v4);%插值figure,surf(X,Y,Z)%三维曲面附录二clcclearclose all D=xlsread(D:c.xls,1,A4:E322);M=xlsread(D:c.xls,2,A4:I322);x=D(:,2);y=D(:,3);for j=1:8 z=M(:,j+

28、1) c=D(:,5); xi=linspace(min(x),max(x),100); yi=linspace(min(y),max(y),100); xi,yi=meshgrid(xi,yi); zi=griddata(x,y,z,xi,yi); ci=griddata(x,y,c,xi,yi); marker=p,s,d,*,o; color=b,m,r,c,y; mat=As,Cd,Cr,Cu,Hg,Ni,Pb,Zn; str=等浓度线,生活区,工业区,山林区,交通区,绿地区; % 等高线图形 figure contourf(xi,yi,zi,0:1:20); title(重金属 ,m

29、atj, 二维等浓度线分布图) xlabel(x(m) ylabel(y(m) colormap summer colorbar grid on hold on for i=1:5 loc=c=i; plot(x(loc),y(loc),markeri,markerfacecolor,colori,MarkerEdgeColor,colori); end legend(str,location,best) for k=1:length(x) text(x(k)-200,y(k)+200,num2str(M(k,j+1),fontsize,10); endend附录三clear allclaa=

30、D(:,5);a1=find(a=1);as=M(:,2);cd=M(:,3);cr=M(:,4);cu=M(:,5);hg=M(:,6);ni=M(:,7);pb=M(:,8);zn=M(:,9);a11=as(a1);a12=cd(a1);a13=cr(a1);a14=cu(a1);a15=hg(a1);a16=ni(a1);a17=pb(a1);a18=zn(a1);a111=sum(a11);a112=sum(a12);a113=sum(a13);a114=sum(a14);a115=sum(a15);a116=sum(a16);a117=sum(a17);a118=sum(a18);

31、b1=size(a11);b2=size(a12);b3=size(a13);b4=size(a14);b5=size(a15);b6=size(a16);b7=size(a17);b8=size(a18);c1=a111/b1(1,1);c2=a112/b2(1,1);c3=a113/b3(1,1);c4=a114/b4(1,1);c5=a115/b5(1,1);c6=a116/b6(1,1);c7=a117/b7(1,1);c8=a118/b8(1,1);d=c1, c2,c3,c4,c5,c6,c7,c8;pie(d,1 0 1 0 1 0 1 0)labels=As 1%,Cd 35%

32、,Cr 8%,Cu 6%,Hg 11%,Ni 2%,Pb 8%,Zn 28%附录四clcclear allM=xlsread(D:c.xls,2,A4:I322);D=xlsread(D:c.xls,1,A4:E322); a=M(:,9); b=find(111=a); for i=1:length(b) xyz(i,:)=D(b(i),1),D(b(i),2),D(b(i),3),D(b(i),4),M(b(i),9) ;end length(b) xlswrite(h.xls,xyz,A1:E153)clc clear alld=xlsread(D:h.xls);a=d(:,5);x,y

33、=sort(a);x0=0;x1=0;y1=0;y0=0;for i=1:3:153n=i;a1=d(y(n),2);a2=d(y(n+1),2);a3=d(y(n+2),2);b1=d(y(n),3);b2=d(y(n+1),3);b3=d(y(n+2),3);h1=b22+a22-a12-b12;h2=2*b2-2*b1;h3=b32+a32-a12-b12;h4=2*b3-2*b1;w1=h1 h2;h3 h4;f1=2*a2-2*a1;f2=2*b2-2*b1;f3=2*a3-2*a1;f4=2*b3-2*b1;w2=f1 f2;f3 f4;w11=det(w1);w22=det(w2);s1=2*a2-2*a1;s2=b22+a22-a12-b12;s3=2*a3-2*a1;s4=b32+a32-a12-b12;k1=s1 s2;s3 s4;g1=2*(a2-a1);g2=2*(b2-b1);g3=2*(a3-a1);g4=2*(b3-b1);k2=g1 g2;g3 g4;k11=det(k1);k22=det(k2);x0=w11/w22;y0=k11/k22;x1=x0;x0=x0+x1y1=y0;y0=y0+y1end