城市表层土壤重金属污染分析本科生毕业设计（论城市表层土壤重金属污染分析二级学院：院专业：数学年级：学号：2作者姓名：指导教师：完成日期：2013年5月3日目录1 引言 (2)1.1 问题由来 (2)1.2 相关信息 (2)1.3 问题的提出 (2)2 基本假设 (2)3 符号说明 (3)4 问题分析 (4)5 模型的准备 (4)6 模型的建立与求解 (5)6.1 问题1的分析与求解 (5)6.1.1 模型的建立 (5)6.1.2 模型的结果 (5)6.1.3 污染程度的单项污染指数评价 (7)6.1.4 综合污染程度的评价模型 (9)6.2 问题2的分析与求解 (10)6.2.1 模型的建立 (10)6.2.2 模型的结果 (11)6.3 问题3的分析与求解 (12)6.3.1 土壤重金属污染物的来源 (12)6.3.2 重金属污染物的传播特征 (13)6.3.3基于matlab的数据处理与三维数据插值模型 (13)6.4 问题4的分析与求解 (19)6.4.1 模型的优点 (19)6.4.2 数据的收集 (20)6.4.3 模型的建立与求解 (20)7 模型的推广 (21)8 结语 (21)参考文献 (22)附录 (23)城市表层土壤重金属污染分析摘要：本文主要对某城市表层土壤地质环境重金属污染的情况进行分析.首先，利用MATLAB绘制出8种重金属元素在该城区的空间分布，通过内梅罗指数评价法来分析该城区内不同区域重金属的污染程度.接着，通过SPSS，用相关性分析和主成分分析，来分析重金属污染的主要原因.然后，建立三次多项式插值模型和高斯模型，确定污染源的位置.最后，运用灰色预测模型对土壤中的重金属的变化进行预测，较好地研究了城市地质环境的演变模式.关键词：重金属污染；内梅罗指数评价法；相关性分析；主成分分析；高斯扩散模型；灰色预测模型中图分类号：X53Analysis of Heavy Metal Pollution of Urban Topsoil Abstract: This paper analyses the heavy metal pollution situation of urban topsoil geological environment. First，using MATLAB to draw the spatial distribution of eight heavy metals in the city, and using nemerow pollution index method for the analysis of the heavy metal pollution of the city zone in different areas. Then，by using SPSS, pollution causes of heavy metals are analyzed with correlation analysis and principal component analysis. After that, establishing cubic polynomial interpolation model and Gaussian diffusion model to confirm the pollution sources. At last, grey prediction model is set up to forecast the changes of heavy metals in the soil, preferably forecasting the evolution of the geological environment model.Key words: Pollution of heavy metals; Nemerow pollution index method ; Correlation analysis; Principal component analysis; Gaussian diffusion model;Grey prediction model.1 引言1.1 问题由来随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出.对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点.1.2 相关信息按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同.现对某城市城区土壤地质环境进行调查.为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土壤（010 厘米深度）进行取样、编号，并用GPS记录采样点的位置.应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据.另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值.附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值.1.3问题的提出(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度.(2) 通过数据分析，说明重金属污染的主要原因.(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置.(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？2 基本假设(1)假设采样数据确定可信，可以真实反映出当地重金属的污染情况.(2)假设暂不考虑影响土壤金属分布的物理性质.(3)假设暂不考虑自然环境对金属浓度的影响.(4)假设重金属元素的背景值可以真实地反映出当地不受污染时的真实值.(5)污染物浓度在y、z轴上的分布符合高斯分布（正态分布）.(6)源强是连续均匀的.(7)在扩散过程中污染物质量是守恒的.3 符号说明4 问题分析这是一个关于分析城市土壤重金属污染分布和寻找污染源的问题.对于给出的海量数据，如何寻找数据间的关系和建立符合数据变化的函数模型是至关重要的，因此，对于问题一，先给出8种重金属元素在该城区的空间分布，并把不同的区域用不同的符号表示，不同浓度的区域用不同的颜色区分，更好地分析该城区内不同区域重金属的污染程度，然后再采用内梅罗指数评价法对该城区进行污染程度评价.对于问题二，通过SPSS进行相关性分析，如果相关性分析的结果大于0.7，就进行主成分分析，找出重金属污染的主要原因，否则，表示主成分分析效果不佳.对于问题三，首先分析重金属污染物的传播特征，由此建立三次多项式插值模型，确定污染源的位置.对于问题四，先分析所建模型的优缺点，再收集有关的信息，更好地研究城市地质环境的演变模式.5 模型的准备[1]根据附件1中的采样数据，用matlab可以画出当地的空间地形图(图1)，编程实现(附录1).图1 该地的空间地形图从以上图中可以看出该地东北方向的地势较高，西南方向地势较低，且是从西南向东北的带状区域.↑生活区①工业区□山区交通区公园绿图2 该地的等高线图用matlab可以画出该地的等高线图(图2)，编程实现(附录2).从图2中，可以看出各区域在该地的地形分布，通过此图较容易分析8种重金属在该地的空间分布，以及不同区域的污染程度，为问题1的求解做好准备.6 模型的建立与求解6.1 问题1的分析与求解6.1.1.模型的建立[1]根据问题分析及模型准备，为了比较容易的观察出8种重金属的在城区分布，分别把每种重金属在采样点的浓度与各区域的划分画在同一个图上. 编程实现(附录3).6.1.2.模型的结果8种重金属的空间分布等高线图如图3：AsCdCr Cu HgNiZnPb图3 8种重金属在该城区的分布图从图3中可以看出对该地污染较大的重金属有：砷、镉、铅、锌；对该地造成中等污染的重金属有：汞；对该地污染较轻的重金属有：铬、铜、镍.其中，在生活区，污染较严重的重金属是：镉、铬、铜、锌；污染较轻的重金属是：砷、汞、镍、铅.在工业区，污染较严重的重金属是：汞、铜、锌；污染较轻的重金属是：砷、镉、铅、铬、镍.在山区，各种金属对其的影响都很小，其中，相对大一点的是：镍.在交通区，污染较严重的重金属是：汞、铜、锌、镉，污染较轻的重金属是：砷、铬、镍、铅.在公园绿地区，污染较严重的重金属是：汞、镉、镍、锌，污染较轻的重金属是：砷、铬、铜、铅.上述结果只是从图中观察得出，可信度不算很高，下面用污染程度的单项污染指数评价法来评价各种金属在各区域的污染程度.6.1.3. 污染程度的单项污染指数评价法[2]要知道各种重金属在不同功能区中的污染程度，首先，要知道各种重金属在不同区域的采样点中的实测值.根据实测值与背景值，就可以用单因子指数法来评价各种重金属在不同区域中的污染程度，其计算公式：=CiPiSi，其中，Pi表示重金属的污染指数，Ci表示重金属的实测值，Si表示重金属的评价标准.对照《土壤环境质量标重准》，土壤重金属污染评价的分级标准为，Pi<1:未污染;1<2Pi≤:轻度污染;2<3Pi≤:中度污染;3Pi≥:严重污染.下面用Excel 统计得出不同功能区中重金属元素的含量平均值，并将得出的数据作为Ci，8种重金属的背景值分别为：3.6，0.13，31，13.2，0.035，12.3，31，69 (mg/kg)，把它们作为Si.得到各种金属在不同区域中的污染程度，如表2.表1 各种金属不同区域中的平均值(mg/kg) 元素As Cd Cr Cu Hg Ni Pb Zn 生活区 6.27 0.29 69.02 49.40 0.093 18.34 69.11 237.01 工业区7.25 0.39 53.41 127.54 0.64 19.70 93.04 277.93 山区 4.04 0.15 38.96 17.32 0.041 15.45 36.56 73.29 交通区 5.71 0.36 58.05 62.21 0.45 17.62 63.53 242.85公园绿地区6.25 0.28 43.64 30.19 0.115 15.29 60.71 154.24 表2 各种金属不同区域中的污染指数元素 As Cd Cr Cu Hg Ni Pb Zn 生活区 1.74 2.23 2.23 3.74 2.66 1.49 2.23 3.43 工业区 2.01 3.02 1.72 9.66 18.35 1.60 3.00 4.03 山区 1.12 1.17 1.26 1.31 1.17 1.26 1.18 1.06 交通区 1.59 2.77 1.87 4.71 12.77 1.43 2.05 3.52 公园绿地区1.742.162.162.293.291.241.962.24从表1中，可以得出，在生活区，Cu 、Zn 属于严重污染，Hg 、Cd 、Cr 、Pb 属于中度污染，As 、Ni 属于轻度污染；在工业区，Hg 、Cu 、Zn 、Cd 、Pb 属于严重污染，As 属于中度污染，Cr 、Ni 属于轻度污染；在山区，污染较轻，8种重金属都属于轻度污染；在交通区，Hg 、Cu 、Zn 属于严重污染，Cd 、Pb 属于中度污染，Cr 、As 、Ni 属于轻度污染；在公园绿地区，Hg 属于严重污染，Cu 、Zn 、Cd 属于中度污染，Pb 、As 、Cr 、Ni 属于轻度污染.6.1.4.综合污染程度的评价模型[3]由于单因子指数评价法只能估计每种重金属在某一区域的污染程度，但如果要评价某一区域的综合污染程度，就要把各种重金属在同一区域的污染指数结合起来，下面就用内梅罗指数评价法来评价五个功能区的污染程度. 综合污染评价不仅兼顾了单因子污染指数平均值和最高值，还可以突出污染较重的重金属污染物的作用.其计算公式为：,22max2i P P P +=）（综其中P 综是采样点的综合污染指数；max i P 为i 采样点重金属污染物单项污染指数中的最大值；∑==ni i P n P 11为单因子污染指数平均值.采用《土壤环境质量标准》中国家一级标准作为环境质量的标准，8种重金属的标准值分别为：15，0.2，90，35，0.15，40，35，100(mg/kg).用Excel 表格统计得到各种金属在不同区域中的单项污染程度，如表3.表3 各种金属不同区域中的单项污染指数元素 AsCdCrCuHgNiPbZn生活区0.4180 1.4500 0.7669 1.4114 0.6200 0.4585 1.9746 2.3701工业区 0.4833 1.95000.59343.64404.2667 0.4925 2.6583 2.7793山区 0.2693 0.7500 0.4329 0.4949 0.2733 0.3862 1.0446 0.7329交通区 0.3807 1.8000 0.64501.7774 3.0000 0.4405 1.81512.4285公园绿地区0.41671.4000 0.4849 0.8626 0.7667 0.3822 1.7346 1.5424得出上表后，利用Excel 表格统计得出单因子污染指数平均值P ，再通过内罗梅综合评价公式计算得到各种金属不同区域中的综合污染指数，如表4，编程实现(附录4). 内梅罗指数以1、2、3、4、5 为界，划分为清洁、轻污染、污染、重污染、严重污染 5 个等级.表4 各种金属不同区域中的综合污染指数功能区污染指数污染等级生活区 1.8733 轻污染工业区 3.3653 重污染山区0.8341 清洁交通区 2.3832 污染公园路地区 1.3980 轻污染从上表可以看出，工业区的污染最严重，污染等级为重污染，其次，交通区，污染等级为污染，生活区和公园绿地区的污染较轻，为轻污染，山区的污染最轻，清洁.6.2 问题2的分析与求解[4]由问题的求解中，知道山区的人为污染比较轻，而其他四个不同功能区都受到不同程度的人类活动污染.其中，生活区的重金属污染可能来源于污水的排放、煤烟的排放和垃圾的污染，工业区的重金属污染可能来源冶炼金属的废气、废水排放和废弃物中的金属污染，公园路地区的重金属污染可能来源于农药的使用.6.2.1.模型的建立由于不同的重金属可能来源于相同的污染源，下面先进性KMO检验.该检验从比较观测变量之间的简单相关系数和偏相关系数的相对大少出发，其值的变化从0到1，一般KMO大于0.9说明适合做因子分析，如果KMO大于0.7，说明还可以做因子分析，若KMO 过小，表明不适合做因子分析.相关性检验结果如表5，编程实现(附录5).从表5中可以看出KMO值为0.778，表明进行主成分分析的效果还可以，下面用SPSS 进行主成分分析，得到表6、7、8、图4，编程实现(附录6).6.2.2.模型的结果表6 公共因子方差表6显示了公共因子方差，初始公共因子方差表示因子提取前各个变量的全部公共因子的载荷系数平方和，提取公因子方差，是根据某种原则提取的公共因子，公共因子数小于等于变量数.图4 特征碎石图表7表示全部解释方差表，图4显示主成分分析中特征的变化情况.图中曲线存在一个明显的拐点，保留前3个主成分能够概括原始数据的绝大部分信息，再结合表7，前3个因子贡献占总体方差为70.941%.表8 主成分载荷矩阵主成分1 2 3As .426 -.200 .681Cd .711 .281 .282Cr .735 -.444 -.303Cu .756 .125 -.365Hg .408 .673 -.297Ni .723 -.515 -.190Pb .764 .314 .237Zn .699 -.037 .123表8表示主成分载荷矩阵，从表8可以看出，Pb 、Cu、Cr、Ni、Cd、Zn在第一主成分的载荷较大，即与第一主成分的相关系数较高；Hg在第二主成分上的载荷较大，即与第二主成分相关程度较高；As在第三主成分上的载荷较大，即相关程度较高.由主成分分析表和问题1的各种金属不同区域中的单项污染指数，可以得到下面的结论：Pb 、Cd、Zn、As主要来源于工业区，主要是源于颜料厂、冶金厂、电镀厂的废气、废水、废渣的排放，垃圾的焚烧、煤的燃烧和农药的使用. Cu、Hg的污染主要来源于工业区和交通区，主要是源于含Pb汽油的燃烧和汽车轮胎磨损产生的粉尘. Cr、Ni主要来源于工业区、生活区和山区，主要是源于电镀、皮革、颜料、劣质化妆品原料等铬化合物制造企业所排放的“三废”和生活所用的化妆品、照明用灯、燃煤以及山区不合理的采矿等.在公园路地区，Cu、Hg、Cr、Cd都有所污染，可能跟过往的车辆有关.6.3 问题3的分析与求解6.3.1 土壤重金属污染物的来源土壤重金属污染物主要来工业废水、废气的排放和生活污水的排放、农药的施用和生活垃圾堆放、焚烧，以及矿产资源的不合理开采和炼制等.6.3.2 重金属污染物的传播特征重金属污染物在土壤中传播属于土壤中溶质的传播.土壤中的重金属污染物与其它的污染物具有不一样的物理化学性质，因此也具有不一样的环境效应.远距离污染源对环境土壤的影响是一个十分复杂的问题，因为污染源周围土壤中污染物的含量受到许多因素的支配，如污染物的状态、排放浓度和总量、污染源的高度、气温降水、地貌特征、土壤的组成和性质等等都会产生影响.因此，影响重金属污染物的扩散有很多因素.但是，从问题1的图3来看，8种重金属在浓度较高的点的附近的分布多数都呈同心圆，而在浓度较低的点的附近没有明显的扩散趋势，说明相应扩散的主要原因不是浓度.因此，可以利用 MATLAB软件的griddata命令下的'v4'三次多项式插值法分别对8种重金属污染物的浓度进行数据插值，并观察得出图像的极大值，进而为确定污染源的位置提供依据.6.3.3基于matlab的数据处理与三维数据插值模型[1]1. 模型的建立利用MATLAB软件的griddata命令下的'v4'三次多项式插值法分别对8种重金属污染物的浓度进行数据插值，编程实现(附录7).画出图像后，观察图像的极值点，根据附件1、2的数据，利用MATLAB的搜索函数find，先从重金属的最高浓度开始查找，然后找到该点对应的x，y坐标，最后把找到的坐标跟图像对比，符合的就是图像的极值点，找到第一点后，再用次高浓度查找，直到找到所有的极值点为止.2. 模型的结果图5 As的浓度图从图5中，可以看出As的浓度图中分别有2个极值点，可以大概确定极值点的位置为：(18134，10046)、(12696，3024).图6 Cd的浓度图从图6中，可以看出Cd的浓度图中分别有1个极值点，可以大概确定极值点的位置为：(21439，11383).图7 Cr的浓度图从图7中，可以看出Cr的浓度图中分别有2个极值点，可以大概确定极值点的位置为：(3299，6018)、(4592，4603).图8 Cu的浓度图从图8中，可以看出Cu的浓度图中分别有2个极值点，可以大概确定极值点的位置为：(2383，3692)、(3299，6018).图9 Hg的浓度图从图9中，可以看出Hg的浓度图中分别有4个极值点，可以大概确定极值点的位置为：(2383，3692)、(2708，2295)、(13694，2357)、(15248，29106).图10 Ni的浓度图从图10中，可以看出Ni的浓度图中分别有1个极值点，可以大概确定极值点的位置为： (3299，6018).图11 Pb的浓度图从图11中，可以看出Pb的浓度图中分别有2个极值点，可以大概确定极值点的位置为： (1647，2728)、(4777，4897).图12 Zn的浓度图从图11中，可以看出Zn的浓度图中分别有5个极值点，可以大概确定极值点的位置为： (9328，4311)、(13797，9621)、(12696，3024)、(4948，7293)、(2383，3692).根据各种重金属浓度图中的极值点，可以确定各种重金属污染源的位置，如表9.表9 各种重金属污染源的位置元素污染源位置/mAs (18134，10046)、(12696，3024)Cd (21439，11383)Cr (3299，6018)、(4592，4603)Cu (2383，3692)、(3299，6018)(2383，3692)、(2708，2295)、(13694，2357)、(15248，Hg29106)Ni (3299，6018)Pb (9328，4311)、(13797，9621)Zn(9328，4311)、(13797，9621)、(12696，3024)、(4948，7293)、(2383，3692)由上表可以看出，As 和Zn 有共同的污染源，来自交通区的(12696，3024)；Cr 、Cu 、Ni 有共同的污染源，来自交通区的(3299，6018)；Cu 、Hg 、Zn 有共同的污染源，来自工业区的(2383，3692)；Pb 、Zn 有共同的污染源，来自生活区的(9328，4311)和交通区的(13797，9621).3.模型的改进[5]根据重金属在土壤中的扩散特点，该问题的解决还可以选用高斯扩散模型进行建模. 污染源选为坐标原点，污染处重金属元素的浓度记为(,,,0)C x y z ，则时刻s 处的无穷空间中的任一点重金属浓度为 (,,,)i i i i C x y z s ，单位距离内通过单位法向面积的流量与浓度梯度成正比，则i q gradC δ=-⋅()=1,2,3,i i δ⋅⋅⋅是扩散系数，grad 表示梯度，负号表示由浓度高向浓度低的地方扩散. 考察空间域Ω ，其体积为V ，包围Ω 的曲面为1S ， 1S 为一规则的球面， 1S 外法线向量为=-,-,1x y n z z ⎛⎫⎪⎝⎭. 则在(),s s s +∆内通过Ω的流量()11. 1s sss Q q n d ds σ+∆=⎰⎰⎰Ω内重金属元素的增量为()()()2,,,-,,, 2VQ C x y z s s C x y z s dV =+∆⎡⎤⎣⎦⎰⎰⎰又由于污染物释放的重金属的总量为()00 3s ssQ p dVdt +∆Ω=⎰⎰⎰⎰由质量守恒定律，有()012 4Q Q Q =+根据曲面面积的Gauss 积分公式有()1. 5s Vq n d divq dV σ=⎰⎰⎰⎰⎰其中div 表示散度记号，由()1式至()5式并利用积分中值定理整理得()()222222, 0,-,, 6i i y z C C C C div gradC s x y z s x y z δδδδ⎛⎫∂∂∂∂==++>∞<<∞ ⎪∂∂∂∂⎝⎭ 这是无界区域的偏微分方程，根据假设(5)，初始条件为作用在坐标原点的污染源函数，可记作()0(,,,0)(,,)7C x y z Q x y z δ=0Q 表示释放的重金属元素总量， (),,x y z δ是单位强度的扩散系数.方程()6满足方程()7的解为()()()2223/2,,, , ,,44i i Q x y z C x y z s e i x y z k s k s δπδ++== 其中k 为t 时间内传播的速度.6.4 问题4的分析与求解6.4.1模型的优缺点1.模型的优点(1)模型一8种重金属元素空间分布图比较清晰、直观，分析也比较详细；对污染程度，采用内梅罗单项指数评价法和综合指数评价法，从不同的角度分析，得到相同的结论，更准确的判断出各个区域的污染情况.(2)模型二在进行主成分分析前，先进行相关性分析，判断该模型是否适合进行主成分分析，使得到的数据更具有准确性.(3)模型三使用 MATLAB软件的'v4'三次多项式插值法对8种重金属进行插值，比较形象的观察出各种金属在该城区空间分布的极值点，从而得到污染源的位置.2.模型的缺点(1)模型一中，8种重金属元素空间分布图的点都是离散的点，利用MATLAB的多项式插值法画图，得出的数据有些误差.(2) 模型三的分布中，忽略了污染源对风向、水流的考虑，而且直接从图中观察得出污染源的位置，得出的数据不够准确.6.4.2数据的收集为了更好地研究城市地质环境的演变模式，还应收集不同时间、不同季节中重金属的污染浓度；更深层土壤中的污染浓度；风向、水流、降雨量对重金属污染的影响；生活区、工业区和交通区等重要污染区的污水、废气和废弃物的排放量以及土壤的降解速度.6.4.3模型的建立与求解[6]由于该题目中给出的信息是少量的、不完全，因此，在本问中，选用灰色预测模型，它是建立数学模型并做出预测的一种预测方法.模型原理如下：对于给定的原始数据序列：(0)(0)(0)(0){(1),(2),,()}x x x x N =，经一次累加得到一个新序列：(1)(1)(1)(1){(1),(2),,()}x x x x N =.设(1)x 满足一阶常微分方程(1)(1) 1dx ax u dt()()()()()()111=0.5+0.5-1z k x k x k()()()()01z k az k u +=其中a,u 为待估计参数，a 为 发展灰数，u 为内生控制灰数.设ˆa为待估计参数向量ˆ=a au ⎛⎫⎪⎝⎭. 按最小二乘法求解得：()()1ˆ=y 2T T aB B B - 其中 (0)(0)(0)T ((2),(3),,()),y x x x N =(1)(1)12(1)(1)12(1)(1)12[(1)(2)]1[(2)(3)]1[(1)()]1x x x x B x N x N ⎡⎤-+⎢⎥-+⎢⎥=⎢⎥⎢⎥--+⎢⎥⎣⎦. 把()2式代入1式，并解微分方程，得到灰色预测模型为：(1)(0)(1)[(1)].ai u u x i x e a a-+=-+通过灰色预测模型，可以得出以后不同时间各种重金属在土壤中的浓度，通过分析上述方法得到的数据，就可以预测城市未来的地质环境的演变模式. 编程实现(附录8).7 模型的推广本文中所建立的模型是金属物在土壤中扩散的的模型，该模型可以用在金属物在大气中扩散的模型、石油在海洋中泄露的模型、放射性气体在空间中扩散的模型以及炸弹爆炸的烟雾浓度的扩散模型等.8 结语随着我国经济的发展和城市人口的增加，土壤重金属污染问题越来越受到科学家的关心与研究.本文主要是通过建立数学模型，分析了城市的污染程度、确定了污染源，并建立模型预测城市未来的发展演变.通过本文的研究，发现土壤重金属污染已经比较严重，而且重金属的除去与清洁需要一个相当漫长的过程，因此，人们必须要提高对城市土壤的保护意识，并且做好对土壤重金属的清洁工作.参考文献：[1] 王沫然.MATLAB与科学计算[M].北京：电子工业出版社，2012，10.[2]王新智等.基于单项指数法和模糊综合评价法对松花江吉林市段水质的评价[J].环境科学与管理，2012，9，9期，37卷： 184-187.[3]寇文杰等.内梅罗指数法在水质评价中存在问题及修正[J].南水北调与水利科技，2012，8，4期，10卷：38-47.[4] 吕振通，张凌云.SPSS统计分析与应用[M]. 北京：机械工业出版社，2009，6.[5]石东伟，陈冬娜.高斯扩散模型在确定污染源位置中的应用[J].河南科技学院学报， 2012，4，2期，40卷：55-58.[6] 唐丽芳，贾冬青，孟庆鹏.用MATLAB实现灰色预测GM(1，1)模型[J].沧州师范专科学校学报，2008，6，2期，24卷：36-37.[7] 吴建国.数学建模案例精编[M].北京：中国水利水电出版社，2005，5.[8] 陈泽，占海明.MATLAB在科学计算中的应用[M].北京：电子工业出版社，2011，6.[9] Mark M.Meerschaert.数学建模方法与分析[M].北京：机械工业出版社，2005，6.[10] Frank R. Giordano，etc .First Course in Mathematical Modeling[M]. 北京：机械工业出版社，2009,10.[11] Schneider，Leo J.Mathematical Modelling[J]. The Mathematics Teacher，1996， 89（7）： 608.附录附录1把附件1中的取样点位置数据表命名为data.x=data(:,1); y=data (:,2); z=data (:,3);[X,Y,Z]=griddata(x,y,z,linspace(0,27823)',linspace(0,18449),'v4');figure,surf(X,Y,Z)附录2利用excel 把采样点位置数据按不同区域分开，分别存为data1，data2，data3，data4，data5.以下程序中x,y,z与附录1中表示的数据相同.x1= data1(:,1); y1=data1 (:,2); z1=data1 (:,3);[X,Y,Z]=griddata(x,y,z,linspace(0,27823)',linspace(0,18449),'v4');contourf(X,Y,Z) ,hold onscatter3(x1,y1,z1,'*'),hold onx2= data2(:,1); y2=data2 (:,2); z2=data2 (:,3);scatter3(x2,y2,z2,'ro'),hold onx3= data3(:,1); y3=data3 (:,2); z3=data3(:,3);scatter3(x3,y3,z3,'ks') ,hold onx4= data4(:,1); y4=data4 (:,2); z4=data4(:,3);scatter3(x4,y4,z4,'mp') ,hold onx5= data5(:,1); y5=data5 (:,2); z5=data5(:,3);scatter3(x5,y5,z5,'v')附录3把附件2的8种主要重金属的浓度表存为DATA.由于8种重金属的空间分布等高线图的程序都差不多，下面以As的空间分布等高线图为例.下面程序中的x，y与附录1中表示的数据相同.a=DATA(:,1)；[X,Y,A]=griddata(x,y,a,linspace(0,27823)',linspace(0,18449),'v4');contourf(X,Y,A); hold onx1= data1(:,1); y1=data1 (:,2); z1=data1 (:,3);scatter3(x1,y1,z1,'*'),hold onx2= data2(:,1); y2=data2 (:,2); z2=data2 (:,3);scatter3(x2,y2,z2,'ro'),hold onx3= data3(:,1); y3=data3 (:,2); z3=data3(:,3);scatter3(x3,y3,z3,'ks') ,hold onx4= data4(:,1); y4=data4 (:,2); z4=data4(:,3);scatter3(x4,y4,z4,'mp') ,hold onx5= data5(:,1); y5=data5 (:,2); z5=data5(:,3);scatter3(x5,y5,z5,'v')附录4矩阵a 表示下表中各种重金属在不同区域中的单项污染指数元素As Cd Cr Cu Hg Ni Pb Zn 生活区 0.4180 1.4500 0.7669 1.4114 0.6200 0.4585 1.9746 2.3701工业区 0.4833 1.9500 0.5934 3.6440 4.2667 0.4925 2.6583 2.7793山区 0.2693 0.7500 0.4329 0.4949 0.2733 0.3862 1.0446 0.7329交通区 0.3807 1.8000 0.6450 1.7774 3.0000 0.4405 1.8151 2.4285公园绿地区 0.4167 1.4000 0.4849 0.8626 0.7667 0.3822 1.7346 1.5424 生活区：P = mean(a(1,:)); max i P =2.3701; P 综=sqrt((P ^2+max i P ^2)/2) 工业区：P = mean(a(2,:)); max i P =4.2667; P 综=sqrt((P ^2+max i P ^2)/2) 山区：P = mean(a(3,:)); max i P =1.0446; P 综=sqrt((P ^2+max i P ^2)/2)交通区：P = mean(a(4,:)); max i P =3.0000; P 综=sqrt((P ^2+max i P ^2)/2) 公园路地区：P = mean(a(5,:)); max i P =1.7346; P 综=sqrt((P ^2+max i P ^2)/2) 附录5根据SPSS 统计分析与应用书籍中的步骤操作，得出以下程序.GET DATA/TYPE=XLS/FILE='D:\课件\毕业论文\A\浓度表.xls'/SHEET=name '附件2'/CELLRANGE=full/READNAMES=on/ASSUMEDSTRWIDTH=32767.DATASET NAME DataSet1 WINDOW=FRONT.FACTOR/VARIABLES V1 V2 V3 V4 V5 V6 V7 V8/MISSING LISTWISE/ANALYSIS V1 V2 V3 V4 V5 V6 V7 V8/PRINT CORRELATION KMO EXTRACTION ROTATION/FORMAT SORT/PLOT EIGEN/CRITERIA FACTORS(8) ITERATE(25)/EXTRACTION PC/CRITERIA ITERATE(25)/ROTATION VARIMAX/SAVE REG(ALL)/METHOD=CORRELATION.p{color:0;font-family:Monospaced;font-size:13pt;font-style:normal;font-weight:normal;text-decoration:none}附录6GET DATA/TYPE=XLS/FILE='D:\课件\毕业论文\A\浓度表.xls'/SHEET=name '附件2'/CELLRANGE=full/READNAMES=on/ASSUMEDSTRWIDTH=32767.DATASET NAME DataSet1 WINDOW=FRONT.FACTOR/VARIABLES V1 V2 V3 V4 V5 V6 V7 V8/MISSING LISTWISE/ANALYSIS V1 V2 V3 V4 V5 V6 V7 V8/PRINT INITIAL EXTRACTION/PLOT EIGEN/CRITERIA FACTORS(3) ITERATE(25)/EXTRACTION PC/ROTATION NOROTATE/SAVE REG(ALL)/METHOD=CORRELATION.p{color:0;font-family:Monospaced;font-size:13pt;font-style:normal;font-weight:normal;text-decoration:none}附录7由于8种重金属的空间分布图的程序都差不多，下面以As的空间分布图为例.x=data(:,1);y= data(:,2); a=DATA(:,1);[X,Y,A]=griddata(x,y,a,linspace(0,27823)',linspace(0,18449),'v4');mesh(X,Y,A);hold on;plot3(x,y,a,'o')附录8由于8种重金属的高斯扩散模型的程序都差不多，下面以As的高斯扩散模型为例. x0= DATA(:,1);n=length(x0);for i=1:nx1(i)=sum(x0(1:i));endfor k=2:nz(k)=0.5*x1(k)+0.5*x1(k-1);endfor i=1:n-1b(i,1)=-z(i+1); y(i)=x0(i+1);endb(:,2)=1; y=y'; au=b\y; ycl(1)=x0(1);for k=1:nc=x0(1)-au(2)/au(1);ycl(k+1)=c*exp(-au(1) *k)+au(2)/au(1);endyc0(1)=x0(1);for k=1:nyc0(k+1)=ycl(k+1)-ycl(k);enddisp(uint16(yc0(2:1:n+1)))；。
