2011A数学建模城市污染物.doc

2011高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 （打印并签名）：1. 2. 3. 指导教师或指导教师组负责人（打印并签名）： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2011高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：城市表层土壤重金属污染分析摘 要随着工农业的不断发展，我国土壤污染问题日益突出，进入土壤的持久性有机污染物、重金属污染等与日俱增，土壤环境面临着前所未有的压力。因此，了解土壤污染物来源是切实有效地控制土壤污染，保障环境安全和农业可持续发展的重要前提。本文针对城市表层土壤重金属污染问题，在合理的假设下，先绘制出重金属元素空间分布图；再根据单因子指数与内梅罗综合污染指数分析方法对各功能区污染程度分析，说明重金属污染的主要原因；接着，在此基础上建立二维污染物空间扩散模型，采用搜索算法在污染较为严重的区域搜索出适量的样本点，再采用多元非线性回归方法对二维扩散模型进行拟合，求解相应的参数；最后通过对第三问模型优缺点的分析，建立时空扩散模型，更加精确反映城市地质环境的演变情况。对于问题1，根据附件数据，调用Matlab软件中的surf函数绘制出重金属元素空间分布图，采用三维图像表示城区地貌特征，用灰度描绘该地受污染程度，清晰反映出在不同地貌，不同功能区的重金属污染物空间分布特点。选用单因子指数和内梅罗指数模型评价污染状况，运用Visual C+ 6.0编程计算出不同功能区单因子指数与内梅罗综合指数，由数据反映出工业区和交通区Hg和Cu严重超标，各功能区染程度由重到轻依次为：工业区、交通区、生活区、公园绿地区、山区。对于问题2，对各功能区的污染物含量进行分析，分析推测该城市存在规模较大或者数量较多的Hg,Cu,Zn严重超标的工厂。为进一步揭示表层土壤重金属污染来源，对各个功能区的不同重金属污染元素进行相关性分析，运用spss软件计算出各相关系数，揭示了不同功能区内土壤表层重金属污染的不同来源，分析可得，交通区污染主要来源于汽车尾气排排放，工业区污染主要来源于工厂的超标排放。对于问题3，根据污染物扩散的特征，建立一维污染物扩散模型：，以此推广建立污染物的二维空间扩散模型：；再运用搜索算法在局部极值点附近搜索适量样本点；采用多元非线性拟合方法对二维空间扩散模型进行拟合，得到区域的最大值，即为污染源所在地,以Hg的分布为例说明。对于问题4，通过对上述模型优缺点的分析，进一步搜集相关数据，依据重金属传播所具有的富集性和不易分解性等特征，在空间分布模型的基础上引进时间变量，突出时间对于污染物扩散的影响，建立时空扩散模型，以更好地研究地质环境的演化模式。本文描述了城市土壤重金属在城市不同功能区的含量分布、不同功能区各元素间分布格局的差异性、给出定量描述城区不同区域重金属的污染程度评价模型，通过对数据的分析以及数据间相关性的分析推测不同污染物的来源，并根据重金属的污染物的传播特征，建立确定污染源位置的模型，为合理规划和利用城市土壤、改造和提高城市环境质量、保障人类健康等提供重要依据。关键词：城市土壤 重金属污染 单因子指数 内梅罗综合指数 扩散模型 非线性拟合一、 问题重述11背景资料背景值：未受人类活动明显影响的土壤本身的化学元素组成和含量称为土壤背景值，或土壤环境背景值。影响土壤背景值的因素很复杂。包括数万年以来人类活动的综合影响，风化、淋溶、淀积等地球化学作用的影响，生物小循环的影响，母质成因、质地与有机物含量的影响等等。因此，土壤背景值是个范围值，而不是个确定值。1.2问题描述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（010 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。1.3需要解决的问题：问题（1）：给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。问题（2）：通过数据分析，说明重金属污染的主要原因。问题（3）：分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。问题（4）：分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、 问题分析2.1 问题1： 先对采集点做分布进行分析，若分布较不均匀，则需做空间插值处理，随后根据附录数据并运用matlab的surf函数进行绘图，得到用三维坐标反映四维信息图像，其中以白色到黑色的渐变过程描绘污染程度的由轻到重的变化。再采用最为常用的土地污染评价方法单因子指数法和内梅罗综合污染指数法对五个不同功能区的污染状况进行评价。2.2 问题2： 由问题1得到各功能区不同种类的重金属含量的单因子指数，并由此分析其重金属污染物的主要类型，推测污染物来源的性质；因为对不同污染物的之间的相关性分析可以帮助分析污染源的性质，通过选择污染较为严重的几种元素进行相关性分析，放映了不同功能区污染来源的差异，并推测各功能区重金属污染的主要原因。2.3 问题3：通过污染物的传播特点，建立二维污染物空间分布模型；考虑到污染源的位于污染严重的附近区域，仅在污染严重的区域搜索样本点，极大地简化算法，且能准确反映污染源的位置；采用多元非线性拟合方法求解二维高斯分布函数，确定局部极值所对应的坐标轴，即污染源所在地；对相关参数进行解释，说明重金属污染物的传播特点。2.4问题4： 对模型一和模型三存在的优缺点进行评价，通过对模型优缺点的分析，了解模型存在的问题，收集更多的信息，借用这些信息对模型进行改善，更好地研究城市地质环境的演变模式。思路的流程图如下：三、 符号说明 ：表示第i重金属的单因子指数 ：表示第i功能区的内梅罗综合指数：i从1到n分别表示n个一维观测点 ：表示污染源中心：表示污染源排污总量：表示第t年污染源的排污总量：表示观测点位置为x的污染物含量：表示二维污染源分布的中心位置：i从1到n分别表示n个二维观测点4、 模型的假设假设1:样本数据真实可靠，可以客观评价城市重金属的分布状况。假设2:各重金属污染物之间没有相互反应，传输中没有引起某种物质的消除和形成。假设3:各重金属污染物在传播的过程中总含量没有衰减。假设4:该城市重金属均由本城市市民生产和生活活动产生，无外来重金属的影响。假设5:污染严重的局部地区仅存在唯一的污染源。假设6:污染源的的传播作用范围有限。5、 模型的建立与求解数据的采集：对采集点做散点图分布，观测采集点的地区分布，如下所示：图1采集点在地区空间上的分布相对均匀，基本覆盖城市的生活区，工业区，山区，交通区和公园绿地区，获得的信息能够较为全面的反映该城市各不同功能区的土壤污染状况。5.1问题1模型的建立与求解为了分析8种重金属污染的分布情况，首先做出这8种重金属污染的分布图，以对重金属污染分布有总体认识；然后根据所给出的数据，用内梅罗综合污染指数法，定量的分析8种重金属的污染分布。考虑到海拔因素对人类活动存在重大影响，为直观表现在不同地貌下，人类活动对城市不同功能区土壤重金属的影响。以三维图像表示城市的地貌特征，用灰度描绘各种重金属污染物的污染严重程度（其中纯白代表无污染）。8种主要重金属元素在该城区的空间分布如下图所示（源代码见附录1）：图1-1 As污染分布图图1-2 Cd污染分布图图1-3 Cr污染分布图图1-4 Cu污染分布图图1-5 Hg污染分布图图1-6 Ni污染分布图图1-7 Pb污染分布图图1-8 Zn污染分布图从各重金属空间分布图可大致看出：1、 山区和公园绿地区重金属污染程度较轻，生活区、交通区和工业区污染较为严重。2、 Hg分布较广，各功能区污染程度不一；各区As、Cd、Ni含量普遍较低。3、 重金属污染集中分布在海拔较低处，污染程度随着海拔的升高而减轻。4、 工业区呈现混合污染态势。1）污染物评价模型的建立土壤重金属评价方法包括生态危害指数法、地积累指数法、单因子指数法和内梅罗综合污染指数法、模糊贴近度方法、潜在生态危害指数法等。因为内梅罗污染指数反映了各污染物对土壤的作用，并且特别考虑了污染最严重的因子，突出了高浓度污染物对土壤环境质量的影响，其环境质量指数在加权过程中避免了权系数中主观因素的影响。故本次研究所采用目前在土壤重金属研究领域被广泛运用的单因子指数法和内梅罗综合指数法来进行重金属污染等级划分。1）单因子指数法： (1.1)为土壤中i种重金属的单项污染指数为土壤中i种重金属的实测值为i种重金属的土壤环境环境背景值2）内梅罗综合污染指数： (1.2)j是城市功能区编号i为重金属元素种类n为重金属污染物种类是j功能区单项污染因子2）污染物评价模型的求解及分析 对数据预处理：将附录1各观测点功能区分类B,C,D,E列复制黏贴到附录2中第J列，将数据复制黏贴到文件名为a.txt的txt文档中，保存路径如：“C:/a.txt”。利用Visual C+ 6.0编程实现各功能区单因子指数和内梅罗综合指数的计算（源代码见附录2）。 单因子指数AsCd Cr Cu Hg Ni Pb Zn 功能区11.74182.23052.22643.74272.65831.49122.22923.4349功能区22.01433.02391.72299.661818.3531.61073.00134.0279功能区31.12341.17171.25681.31191.17021.25641.17921.0622功能区41.58562.76931.87274.713312.76641.43232.04953.5196功能区51.73992.1581.40762.28733.28551.24311.95832.2354表1-1各功能区不同重金属的单因子指数等级划分1234单因子指数11.0-2.02.0-3.03.0污染等级非污染轻污染中污染重污染表1-2单因子指数污染程度分级注：来源于中国绿色食品发展中心绿色食品产地环境质量现状评价导则2000 年由数据发现，在工业区和交通区的重金属Hg含量严重超标，可推测该城市存在大型化工厂或电池制造厂等高排放量工厂。内梅罗综合指数功能1功能2功能3功能4功能5Min2.09049.27690.82426.46251.8689Average3.170613.5331.25319.42642.7344Max6.675724.96532.826317.36675.0243表1-3内梅罗综合指数等级划分12345综合污染指数0.70.7-1.01.0-2.02.0-3.03.0污染等级安全警戒级轻污染中污染重污染表1-4内梅罗综合指数污染程度分级（注：来源于中国绿色食品发展中心绿色食品产地环境质量现状评价导则2000 年） 图1-9 各不同功能区的内梅罗综合系数5.2问题2分析1）数据分析：内梅罗综合污染指数受污染最严重的因子影响较大，一方面突出表现了高浓度的重金属元素Hg对工业区和交通区土壤环境质量的污染，另一方面也掩盖了其它重金属元素对该地区的影响。因此，我们对剔除Hg后的数据再次进行计算，将两次结果进行对比，更加全面地反映出不同功能区的受污染情况。内梅罗综合指数功能区1功能区2功能区3功能区4功能区5Min2.08144.73440.82392.47651.4318Average3.16017.28591.25463.79372.0852Max6.66215.78052.83318.06014.2902表2-1排除Hg后的内梅罗综合指数剔除Hg影响后，发现该城市工业区，交通区，生活区受其它重金属元素污染严重；公园绿地区污染情况也不容忽视，远离城市的山区受污染相对较轻。2）通过数据特征，定性描述重金属污染原因背景资料：城市土壤重金属的来源存在两种途径：一种是来自于成土母质；另一种是外源输入，重金属分布与城市化过程中人类的活动密切相关。大量的城市土壤研究表明，城市重金属污染主要涉及As、Cd、Cr、Cu、Ni、Hg、Pb、Zn。外源输入的重金属主要来源于： 工业污染源。主要是基础工业特别是采矿和冶炼业是向环境中释放重金属的主要污染源。 生活污染源。生活垃圾中含有的重金属能渗到土壤中;生活废水中的有机物可以吸收重金属,用这些水灌溉也可污染土壤。 交通污染源。机车的燃料中所含的重金属如对土壤的污染。工业污染源采矿冶金等行业三废排放：Hg、Cr、Cd、As、Pb、Zn交通污染源汽车尾气排放：Pb、Zn、Cd、Cr、Cu生活污染源生活废弃物：Cr、Cu、Zn、Pb 表2-2 各行业主要污染物类型注:中国环保信息网，城市土壤的重金属污染及生物治理，2010年 7对照表一，Cr,Ni,Cd,As含量较为接近相对背景值，由背景资料知道这几种元素受母岩控制较大，而人类活动对其影响较小。而Hg,Cu,Zn,Pb等重金属含量在人类活动频繁的地势较为平坦的区域内含量明显增加，其中Hg，Cu含量明显幅度较大，工业区Hg平均值超标18.4倍,Cu超标9.7倍，交通区Hg平均值超标12.8倍，Cu超标4.7倍，说明该城市存在Hg，Cu排放严重超标的企业或工厂，并且其排放的污染物已经对该城市造成严重影响。各重金属间的相关性分析：各土壤重金属之间的依存关系，能帮助揭示重金属元素的污染来源，选用受人类活动影响较大的Hg,Cu,Pb和Zn为例。分析元素之间的相关性，分析污染的来源。由于Hg,Cu,Pb,Zn其数量值表现出等级特征，故采用spearman分析Hg,Cu,Pb,Zn之间的等级相关性。使用spss中Correlate 命令项spearman分别对功能区一，功能区二和功能区四的Hg，Cu，Pb,Zn做相关性分析，得到下表：相关系数RHg-CuHg-PbHg-ZnCu-PbCu-ZnPb-Zn功能区一0.6450.7360.7740.7440.7730.877功能区二0.5930.4110.5730.8030.860.778功能区四0.6230.6120.6790.7820.8840.891表2-3由上图可知：第一：Cu，Zn，Pb间的相关性较强，其中Zn,Pb自相关性较强，有资料表明受交通污染源控制时，交通道路断面土壤重金属污染程度取决于其距离交通道路的远近。离公路越远，受交通污染源的影响越小，土壤Pb、Zn污染越轻，递减趋势明显。土壤中Zn和Pb含量随道路远近的变化趋势具有较强的一致性，说明土壤Pb和Zn源于相同的污染源为交通污染源。第二：不同功能区土壤重金属间依存关系各不相同，反映各地区重金属来源的途径存在差异。在工业区，部分土壤重金属间依附关系不明显，而部分金属间的依附关系显著，这反映该功能区受到复合型污染源（即同时排放不同重金属元素的污染源）影响，而公园绿地区、山区重金属污染相对小。综上所述:该城市工业区，交通区，生活区污染均较为严重，其中重金属元素Hg含量最为突出，可推测该地区存在大型或数量较多的Hg严重超标的企业；Cu，Zn，Pb等重金属元素浓度也严重超标，由于Cu,Zn,Pb在各地区含量均存在很强的相关性，故推测Cu,Zn,Pb等重金属元素来源于超标排放的企业或工厂。不同功能区，反映出污染源的性质不同，如工业区反映重金属污染源的多样性，生活区和交通区反映污染源单一，分别为生活废弃物和汽车尾气排放。5.3问题3模型的建立由于重金属是以一定比例存在于水体、悬浮物和底泥之中，因此，重金属在水体中的迁移将取决于它在这三相中的运动规律。考虑到重金属在水体中溶解态、悬浮态和底泥态的物理运动物质交换，吸附和解吸物质交换，以及沉降和再悬浮物质交换过程复杂，且对城市的河流，悬浮物，底泥等量值不明确，无法给出准确信息建模。故从总体角度将重金属输送迁移基本方程简化为平面扩散模型。因为重金属的传播以污染源为中心，向四周散开，浓度最大的地方就是污染源，将问题转化为求局部最大值。假设在局部足够小的区域内对于某类重金属存在唯一的重金属污染源， 经过足够长得时间后，此污染源对周围地区照成了重金属污染，在该地区寻找n个观测点，其重金属污染物含量为，对该过程的模拟图如下所示： 图3-11）一维污染源传播模型：对于这一特定的重金属污染源,假设污染源的总排放量为，且污染物在空间的分布形式如下图所示： 图3-2对污染物的直线扩散方式进行分析，记坐标轴为的观测点的污染物含量为，若受体距离污染源越远，则污染源对其作用效果越不明显，即传播系数为一随距离递减的曲线，则在的观测点，其污染物的一维分布形式满足：其中 为污染物的传播系数，为污染物的传播范围。因为重金属污染在传播过程中，性质稳定，且具有富集的特点，空间传播速度较慢，故认为为一常数，由污染物的性质决定。当有：分离变量求解得到：简化模型： （3.1）发现模型的分布呈倒扣钟形，且近似于高斯分布曲线，其中当时，函数取得最大值，即为污染源位置所在地， 表示污染物的空间分布状态。2）二维的污染物扩散模型重金属污染物分布在平面上，依据上述分析，对一维污染物模型进行推广，忽略空间的各向异性，可以得知重金属在空间上的分布为二维高斯分布形式，其形状为一倒扣钟形，其扩散模型为： （3.2）搜索样本点的算法：取n个样本点， ， 对对函数模型进行拟合。因为局部地区的污染源位置大致位于极大值的所处的周围区域，故选取极大值附近的点作为拟合函数模型的样本点具有实际意义。以Hg为例说明：图3-3 如图3-3所示：存在三个等浓度线密集的区域，表明空间分布明显存在三个污染源地区，分别位于，和编号Hg污染物含量x轴y轴81350023833692131112427397171212883361761040164727281477352643579160002700229513914235182571151138402639131401263469230015024040202990三个区域内，故在该范围内搜索极值点，并选取周围附近10个点作为样本点，算法实现见附录2，在VC6.0运行通过程序，搜索结果分别如下表3-1,3-2,3-3所示：编号Hg污染物含量x轴y轴257138001524891066324.641546786582562715198101006230143258666258841642890692482314000897019815163018299249141420799802593116289100726163.21137979621表3-1 表3-2 编号Hg污染物含量x轴y轴12667471828613737314950017873109132117915800104921276104016472728144142.5326779310632933176791600027082295140126.4334692308表3-3参数拟合求解：对模型两边取对数，得到。对该模型进行拟合，程序代码见附录3。拟合的模型如下污染源1： 污染源2：污染源3：得到三个污染源分别为：污染源1：，污染源2：和污染源3：。经查表分析知道污染源附近均存在污染特别严重的区域，例如：污染源1附近存在点A（15248,9106）污染源2附近存在B(2700,2295)以及C(2383, 3692)，而污染源3附近存在D(2708, 2295),故模型较为可信。 5.4问题4的模型建立问题3中的重金属污染源传播模型以及污染源位置的确定模型：优点：通过污染物的传播特点，建立二维污染物空间分布模型，准确反映了污染物的平面分布状态，考虑到污染源的位于污染严重的附近区域，仅在污染严重的区域搜索样本点，极大地简化算法，且能准确反映污染源的位置。采用多元非线性拟合二维高斯分布函数，确定局部极值所对应的坐标轴，确定污染源所在地。缺点：二维的污染物扩散模型仅反映了空间的污染物积累对于周围地区污染扩散分布的影响，却不能准确反映重金属污染传播的时空分布关系，而污染物扩散模型是对污染物扩散的模拟，是对土壤重金属污染物随时间和空间迁移转化规律的描述，即用于描述土壤中重金属污染物质与时间、空间定量关系的数学方程。1）一维时空扩散模型由重金属污染源的空间分布，确定污染源的大致区域，然后收集该地区工业，交通业发展的状况，引进时间因素t，对二维的污染物扩散模型进行修正，分析过程如图所示： 图4-1如图所示：以一维污染源模型为例说明，其中u为污染源位置，x轴为距离各观测点的位置，y轴表示各位置污染物含量的占总污染物含量的比例。污染物排放的扩散与在空间的分布形式与污染物排放的时间存在相互制约的关系，当污染物排放时间较短，分布形式类似Q1所示：污染物集中分布在u位置及其周围地区；当污染物排放时间较长时，分布形式类似Q3所示：污染物相对均匀分布在u位置及其周围地区。模型（3.2）是对污染物在空间的富集状态的描述，对于某一年份污染物排量在空间上的分布也服从模型（3.2），为了体现不同年份的影响效果，由高斯分布模型知道钟形的形状即污染物的空间分布状况由参数决定。引入时间的影响，即第t年份污染物排放总量对空间分布的影响模型如下： （4.1）已知道随t的增大而增大，且当时，有，可考虑线性模型或指数模型： （4.2）采用上述搜索算法，在局部区域内搜索n个观测点，因为重金属性质稳定，不易发生化学变化等改变物质属性的，且传播过程易出现富集，人工治理难度大，所以各观测点的污染物的含量为污染源多年排放量的线性累积和，假设该地区污染企业存在m年，故得到： （4.3） 2）二维扩散模型假设局部地区的污染源坐标为,对上述一维模型推广应用得到重金属的二维空间分布模型为： （4.4）搜索n个样本点，则这n个点满足上述模型：其中或，选取模型以实际为参照标准。模型综合考虑空间和时间对于污染物扩散的影响，更加精确描述重金属污染物的传播方式，为研究城市地质环境的演变模式提供理论基础。六、 模型的评价优点：1、 问题1采用Matlab绘制三维立体图，并引入第四维颜色量表示污染物浓度，清晰反映出在不同地貌，不同功能区的重金属污染物分布状态。采用内梅罗指数法对城区污染程度进行评价，数学过程简洁，运算方便。2、 问题2从数据独立性和相关性两方面分析污染源性质，清晰反映出主要污染物和污染物源头的性质。3、 问题3采用高斯分布模型模拟重金属污染物的传播特性，并假设局部地区存在单一污染源，极大简化问题和搜索算法，准确确定污染源位置。4、 考虑到污染物扩散的时间积累效果，使结果更加精确，更加贴近实际，得出的结果对重金属污染的防治具有现实指导意义。缺点：1、 重金属的传播模型没有体现基于功能区土壤，气候，水流等对金属污染物传播方式的影响。2、 模型忽视了污染源分布的随机性。改进以及推广 搜集该城市各功能区的土壤，气候，水流等信息，分析传播媒介的主要性质，确定传播系数等数值，建立基于土壤，气候，水流的信息分析的对流扩散模型。以此作为重金属污染的传播模型，能更加准确研究城市地质环境的演变模式。.5.6.参考文献：1 黄芳，水土污染空间分析及源辨析， 2006年。2 李亚松、张兆吉、费宇红、王昭，内梅罗指数评价法的修正及其应用，中国地质科学院水文地质环境地质研究所，2009年。3 郑喜坤、鲁安怀、高翔、赵谨、郑德圣，污染物指标，北京大学地质系、中国地质大学材料学院，2008年。4 韩中庚，数学建模方法及其应用，高等教育出版社，2008年。 5 姜启源，数学模型，高等教育出版社，2003年。6 赵静，但琦，数学建模与数学实验，高等教育出版社，2008年。7 韩中庚，数学建模竞赛，科学出版社，2006年8 中国绿色食品发展中心，绿色食品产地环境质量现状评价技术导则,2001年。9 中国环保信息网，城市土壤的重金属污染及生物治理，/news/gongyegufei/news_47340.html， 2011年9月11日附录1：Matlab 四维作图代码：%As图像%其它图像作图程序类似x=xlsread(cumcm2011A附件_数据.xls,1,B4:B322);y=xlsread(cumcm2011A附件_数据.xls,1,C4:C322);h=xlsread(cumcm2011A附件_数据.xls,1,D4:D322);x1=0:300:30000;y1=0:200:20000;X,Y=meshgrid(x1,y1);H=griddata(x,y,h,X,Y,V4);z1=xlsread(cumcm2011A附件_数据.xls,2,B4:B322);Z1=griddata(x,y,z1,X,Y,V4);col=zeros(101,101,3);for i=1:101 for j=1:101 col(i,j,1)=Z1(i,j)*(-0.0351)+0.9435; col(i,j,2)=Z1(i,j)*(-0.0351)+0.9435; col(i,j,3)=Z1(i,j)*(-0.0351)+0.9435; endendsurf(X,Y,H,col);shading interp附录2：/程序计算各功能区各类金属的单因子指数以为各功能区的内梅罗综合系数/搜索样本点/程序在VC6.0运行通过 #include#include#include#include#include#includeusing namespace std;struct Area /各污染物种类见命名int ID;double As;double Cd;double Cr;double Cu;double Hg;double Ni;double Pb;double Zn;double x,y;double high;int kind;Area area320;void ReadDate()string str;str=C:a.txt; /数据用txt格式保存，路径C:a.txtfstream infile(str.c_str();if(infile.is_open() /成功打开txt文件for(int i=1;iareai.ID;infileareai.Asareai.Cdareai.Crareai.Cuareai.Hgareai.Niareai.Pbareai.Zn;infileareai.xareai.y;infileareai.high;infileareai.kind;area0.ID=0; /area0存入各土壤污染物的背景值，实例为背景值的平均值area0.As=3.6; area0.Cd=130;area0.Cr=31;area0.Cu=13.2;area0.Hg=35;area0.Ni=12.3;area0.Pb=31;area0.Zn=69;elsecoutCan not open the file.;exit(1);double P69;/Pji表示功能区j的第i种金属污染物的平均单因子指标void GetPji() int count6; /countj用于记录j功能区计数memset(P,0,sizeof(P);/数组初始化memset(count,0,sizeof(count);for(int i=1;i320;i+)switch(areai.kind)/对数据进行功能区分类处理case 1:P11+=(areai.As/area0.As);P12+=(areai.Cd/area0.Cd); P13+=(areai.Cr/area0.Cr);P14+=(areai.Cu/area0.Cu); P15+=(areai.Hg/area0.Hg);P16+=(areai.Ni/area0.Ni); P17+=(areai.Pb/area0.Pb);P18+=(areai.Zn/area0.Zn); count1+;break;case 2:P21+=(areai.As/area0.As);P22+=(areai.Cd/area0.Cd); P23+=(areai.Cr/area0.Cr);P24+=(areai.Cu/area0.Cu); P25+=(areai.Hg/area0.Hg);P26+=(areai.Ni/area0.Ni); P27+=(areai.Pb/area0.Pb);P28+=(areai.Zn/area0.Zn); count2+;break;case 3:P31+=(areai.As/area0.As);P32+=(areai.Cd/area0.Cd); P33+=(areai.Cr/area0.Cr);P34+=(areai.Cu/area0.Cu); P35+=(areai.Hg/area0.Hg);P36+=(areai.Ni/area0.Ni); P37+=(areai.Pb/area0.Pb);P38+=(areai.Zn/area0.Zn); count3+;break;case 4:P41+=(areai.As/area0.As);P42+=(areai.Cd/area0.Cd); P43+=(areai.Cr/area0.Cr);P44+=(areai.Cu/area0.Cu); P45+=(areai.Hg/area0.Hg);P46+=(areai.Ni/area0.Ni); P47+=(areai.Pb/area0.Pb);P48+=(areai.Zn/area0.Zn); count4+;break;case 5:P51+=(areai.As/area0.As);P52+=(areai.Cd/area0.Cd); P53+=(areai.Cr/area0.Cr);P54+=(areai.Cu/area0.Cu); P55+=(areai.Hg/area0.Hg);P56+=(areai.Ni/area0.Ni); P57+=(areai.Pb/area0.Pb);P58+=(areai.Zn/area0.Zn); count5+;break;for(int k=1;k6;k+)for(int j=1;j9;j+)Pkj=Pkj/countk; /cout 输出cout 各功能区的各种污染物的单因子系数：endl;cout As Cd Cr Cu Hg Ni Pb Zn endl;for(int l=1;l6;l+)cout第l个功能区: ;for(int j=1;j9;j+)coutfixedsetw(7)setprecision(4)Plj ;coutendl;coutendl;double p6; /pj表示j功能区的内梅罗综合系数void Getpj() double maxvaule;double ave;for(int i=1;i6;i+)ave=Pi1;maxvaule=Pi1;for(int j=2;j9;j+)if(maxvaulePij) /得到j功能区的最大单因子数maxvaule=Pij;ave+=Pij;pi=sqrt(maxvaule*maxvaule)+(1.0/8.0*ave)*(1.0/8.0*ave)/2.0);/cout 输出Pcout排除Hg的影响下各功能区的内梅罗综合系数：endl;for(int k=1;k6;k+)cout功能区k ;coutpkendl;coutendl;void GetpjDeletHg() /求解内梅罗系数（排除Hg含量影响 ）double maxvaule;double ave;for(int i=1;i6;i+)ave=Pi1;maxvaule=Pi1;for(int j=2;j9;j+)if(j=5)continue;elseif(maxvau