1、2011高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置
2、报名号的话）： 13212002 所属学校（请填写完整的全名）： 湖北师范学院 参赛队员 (打印并签名) ：1. 易汀涓 2. 刘慧娜 3. 李端君 指导教师或指导教师组负责人 (打印并签名)： 陈琴 日期： 2011 年 9 月 12 日赛区评阅编号（由赛区组委会评阅前进行编号）：2011高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：城市表层土壤重金属污染分析摘要本文就重金属在某城市五类区域的浓度分布，利用Su
3、rfer软件绘出了八种重金属在该城区的分布图，我们保持单一变量法，对每种元素在整个区域的分布情况，利用会等高线思想汇出图1至图8。问题一根据单因子指数和内梅罗指数（和分别表示污染物的测量值和背景值）确定五个城区的金属污染程度。并利用（表示每种污染物的相对重要性比值，为各评价指标权重）对内梅罗指数进行了加权平均修正，得出：区域名称修正综合污染指数污染等级生活区2.569中度污染工业区8.374重度污染山区1.214轻度污染交通区5.784重度污染公园绿地区2.227中度污染在问题二中，本文先进行横纵比较，初步得工业区Hg的污染程度最严重，然后建立主成分分析模型，利用模型一中的单因子指数，建立矩阵
4、，利用matlab进行处理后，得到相关系数矩阵，计算的八个特征值分别为：6.3757，0.9437，0.4392，0.2415，0，0，0，0。据累计贡献率得出三个主成分，由此确定工业区Hg、Cu，交通区As、Hg，生活区Cu为重金属污染的主要原因。问题三从三个侧面进行数据分析，判断重金属污染物的三个传播特征：以污染源为圆心逐步向外扩散;以工业区为主要传播重心，随以重金属富集明显的路径进行传播;某些重金属不会扩散，其在土壤的稳定性很强，在一定时期内不会传播。参考地震源的位置确定采用Geiger定位方法，考虑单一污染物的浓度分布情况，利用Inglada方法进行初值的计算得（和分别为污染源、观测点
5、坐标，分别为观测点、污染源浓度）。利用Geiger方法迭代校正，确定Cu污染源位置在,的范围。在问题四中，我们充分考虑了模型三的静态特征，还应该收集土壤的固相、液相、气相和PH值以及重金属元素的毒理性的相关信息。此外，我们还化静态为动态加入时间参变量，提高了模型三的准确度。整合以上信息，可以利用差分方程建立模型求解。关键词：单因子指数 内梅罗指数 主成分分析法 Geiger定位法 Inglada方法一、问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下
6、城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（010 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等
7、信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求你们通过数学建模来完成以下任务：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、问题分析本问题是关于环境生物和地址学相结合的模型问题，附件给出了区域、采样点的坐标位置、海拔以及污染物浓度一系列的数据。四个问题从不同的
8、侧面要求我们讨论污染程度、主要原因、传播特征以及污染源的位置确定，最终是希望我们客观全面的研究地质环境的演变模式，从而有效地解决环境污染问题。问题一：基于附件1中采样点的位置和附件2中8中重金属元素在各个区域内的分布情况，运用surfer软件绘出下列8个等浓度图。考虑重金属在该城区的空间分布情况时，我们针对每种元素对不同的5个区域进行了绘制，在分析该城区重金属的污染程度时，我们针对每类区域的8种重金属的浓度，采用单因子指数法和内梅罗综合污染指数法进行污染程度的分析。并且在此基础上，我们对内梅罗指数进行了修正，得出比较客观的修正值，让结果变得更加准确。问题二：首先进行横向比较和纵向比较，初步进行
9、直观判断，然后进行主成分分析，得出三个主成分，和初步分析的结果基本吻合。问题三：分别从不同侧面进行分析，即：纵横坐标、海拔、区域以及重金属污染物固有的一些特性，然后跟确定地震源的确定方法进行比较和借鉴，我们采用了经典的Geiger定位法，并且还运用Inglada方法进行初值的计算。然后再用Geiger进行迭代校正。这样是结果更加精准。问题四：这一问属于开放性问题，所以在模型三的基础上，我们只要加上相关的参数，并且将模型三有静态改为动态模型，这样我们加上时间参变量，对其进行模型改造。利用差分方程进行有效的解决。三、模型假设1、假设每类区表层土壤中每种重金属含量连续变化；2、在对采样点分析时，忽略
10、土壤PH值、温度、气候等因素对测量带来的影响；3、在对每类区域的污染程度进行分析的时候，忽略考察的8种元素以外的重金属污染；4、忽略重金属由于土壤特性和氧化还原等一系列的物理和化学作用，给污染程度评价带来的误差。5、由于不好判断五类区的地理分布情况，我们就将问题理想化，考虑每类区的8种重金属污染程度；6、假设附件给出的数据对于建立模型已经足够充分；7、假设每个采样点的采集时间相同，重金属的含量不会因为采集时间的差异出现变化；8、忽略污染化合物之间的相关影响；9、忽略重金属污染物之间的交叉污染；10、在第三问中，我们将土壤理想化成成分均匀的介质。四、符号说明:土壤重金属污染物的环境质量指数:土壤
11、重金属污染物采样点的个数:土壤重金属污染物的实际测量值:土壤重金属污染物的背景值:本文重金属污染物种类总数：综合环境质量指数：相对重要性比值：各评价指标的权重：综合环境质量指数修正值：污染源坐标：污染源坐标校正量：浓度校正量：测定浓度的偏差五、模型的建立与求解 5.1问题一5.1.1 八种重金属元素的空间分布在考虑8种重金属在该城区的空间分布时，本文抓住问题的核心：单一元素在五类区域的空间分布。这样我们利用附件一和附件二使用matlab软件绘出8种重金属在整个区域的空间分布等高线（见图1-8）。图1 重金属元素As的空间分布图图2 重金属元素Cd的空间分布图图3 重金属元素Cr的空间分布图图4
12、 重金属元素Cu的空间分布图图5 重金属元素Hg的空间分布图图6 重金属元素Ni的空间分布图 图7 重金属元素Pb的空间分布图图8 重金属元素Zn的空间分布图5.1.2不同地区重金属的污染程度在分析该城区内不同区域重金属的污染程度时，本文使用单因子指数法和综合因子指数法。1、评价方法首先采用单因子污染指数法对土壤重金属污染进行评价，计算公式为：然后采用综合因子指数法对土壤重金属污染进行评价，计算公式为：式中，为土壤污染物的环境质量指数；为土壤污染物的实际测量值；为土壤污染物的背景值；是土壤污染物采样点的个数；为重金属种类总数。2、评价标准采用国家土壤重金属评价标准(GB15618-1995)1
13、进行评价，评价标准见表1、2：表1 土壤各元素污染程度分级标准等级划分1234单因子污染指数污染程度非污染轻度污染中度污染重度污染表2 土壤综合污染程度分级标准等级划分12345综合污染指数污染程度安全警戒线轻度污染中度污染重度污染3、评价结果使用EXCEL进行数据分析与处理，得到如下评价结果:表3. 不同区域重金属的污染指数区域单因子污染指数 综合污染指数 AsCdCrCuHgNiPbZn11.742.232.233.742.661.492.233.433.1722.013.021.729.6618.351.613.004.0313.5331.121.171.261.311.171.261.
14、181.061.2541.592.771.874.7112.771.432.053.529.4351.742.161.412.293.291.241.962.24 2.73单因子污染指数评价结果：生活区中，Cu和Zn属于重度污染；工业区中，Hg和Cu的污染程度及其严重；山区各重金属污染程度均较低；交通区中也是Hg的污染相当严重；公园绿地区各重金属污染程度相对还好。 5.1.3修正内梅罗指数法显然，用此方法的计算结果会突出最大污染物产生的影响，但却没有考虑污染因子的危害性差异 2 。会引入很大的误差，因此，将此式修改如下：将污染指数平均值改为加权平均，根据评价指标对环境及人体的危害性来确定各评价
15、指标的权重，具体做法如下：首先将综合评价指标按由小到大的顺序排列为，将最大的一个记为，令表示第种评价指标的相对重要性比值，则为各评价指标的权重。下面运用Excel进行简单计算，结果如下（表4）：表4 修正后的综合污染指数区域名综合指数相对重要性比值各评价指标的权重修正值生活区3.170.2210450.0999192.569工业区13.5310.452038.374山区1.250.0925960.0418561.214交通区9.430.6965490.3148615.784公园绿地区 2.730.2020540.0913352.227修正后得到不同城区的综合污染程度（表5）。表5 不同城区综合
16、污染程度区域名称综合污染指数污染等级生活区2.569中度污染工业区8.374重度污染山区1.214轻度污染交通区5.784重度污染公园绿地区2.227中度污染综合污染指数评价结果：工业区污染指数为8.374，属于重度污染，其次是交通区，污染指数为5.784，生活区和公园绿地区属于中度污染，污染指数分别为2.569和2.227，山区属于轻度污染，污染指数为1.214。5.2问题二5.2.1初步分析通过问题一中对于重金属污染的数据表3和表5，我们分别从横向和纵向进行分析，找出重金属污染的主要原因。I.横向比较横向比较五类综合污染指数，可以清晰地看出，工业区所占比例明显高出其他四类区域，交通区次之，
17、生活区和公园绿化区持平居后，山区的污染程度最轻。这说明工业污染导致重金属污染的成分最重，交通区主要以排放的污染气体为主，被污染的大气长期滞留，这样就会污染到土壤表层。对于生活区和公园绿化区，分析得知，他们的日常生活基本一致，然而公园绿化区的污染程度要比生活区稍轻一些，这根实际情况是相吻合的。下图给出了五类区的污染指数所占比例：图9 五类区域综合污染指数所占比例横向分析可以看出，工业区的综合指数最高，这与实际规律相符，说明工业区周围的治理污染工作还不到位，横向主要原因为工业污染。II.纵向比较纵向比较污染程度最严重的区域（工业区），分析8种重金属的单因子指数，这样可以明显找出导致污染的主要因素。
18、图10 工业区单因子污染指数由图10可以看出Hg的污染指数最高，工业区土壤中Hg明显富集，这与工业区的Hg元素及含Hg的化合物渗透土壤的污染物比较多，下面我们采用主成分分析法进行污染物的主要原因的确定。5.2.2主成分分析先利用问题一处理的结果，建立五个区域和八种元素的五行八列的单因子矩阵然后利用matlab计算X的相关系数方阵R其中为R中行列的元素，其计算公式为由此公式可以知道，然后计算上述方阵的特征值，令，由matlab计算结果为：然后进行排序，这样得到特征向量，得到贡献率：，累计贡献率：，利用matlab计算结果见下表：(表6)表6 特征值、贡献率、累计贡献率特征值6.37570.943
19、70.43920.24150.00000.00000.00000.0000贡献率0.79700.11800.05490.03020.00000.00000.00000.0000累计贡献率0.79700.91490.96981.00001.00001.00001.00001.0000表7 特征向量0.15190.35590.4057-0.2072-0.2221-0.68610.01210.35000.4885-0.1233-0.18680.44260.5757-0.1959-0.06810.3755-0.28440.13620.27510.29910.12340.18040.79140.2476
20、-0.0309-0.53680.51470.2521-0.31490.2269-0.31070.3678-0.36430.50610.1305-0.10970.38370.3307-0.46190.33490.52090.2952-0.2610-0.1244-0.42830.47790.11820.3635-0.4955-0.0429-0.59940.2553-0.3274-0.2628-0.06800.3839-0.0829-0.4541-0.1248-0.71800.25890.00530.20110.3851取累计贡献率达85%95%的特征值所对应的第一、第二第的主成分，那么就得到主成分
21、有三个，分别为第一主成分第二主成分第三主成分: 由此确定污染主要原因为：工业区Hg、Cu，交通区As、Hg，生活区Cu，分别为重金属污染的主要原因。5.3问题三5.3.1重金属污染物的传播特征对于重金属污染物的空间分布，结合数据的变化情况，我们首先从几个侧面分析污染物的传播特征，然后根据最速下降法进行方向追踪污染源的位置。我们整合了题目给的附表一和附表二的数据，建立了表一（见附录），有数据直观分析出，从平面角度：以污染物浓度较高的地理位置一圆心扩散，从区域划分角度：工业区交通区生活区公园绿化区山区。然而，值得注意的是，从元素的种类上分析时，我们发现，对于某些重金属不会扩散，也就是说，其在土壤的
22、稳定性很强，在一定时期内不会传播，那么我们在确定传染源的时候，就进行模糊聚类，将一些传播明显的元素进行梯度分析，进而得到最速下降最快的方向，也就得到了污染源的位置。5.3.2土壤重金属污染来源介绍土壤中重金属元素主要有自然来源和人为干扰输入两种途径。3在自然因素中，成土母质和成土过程对土壤重金属含量的影响很大。在各种人为因素中，则主要包括工业和交通等来源引起的土壤重金属污染。以下主要就受人为作用影响的土壤重金属污染来源进行介绍。1、不同工矿企业对重金属积累的影响工业过程中广泛使用重金属元素，工矿企业将未经严格处理的废水直接排放，使得它们周围的土壤容易富集高含量的有毒重金属。企业排放的烟尘、废气
23、中也含有重金属，并最终通过自然沉降和雨淋沉降进入土壤。矿业和工业固体废弃物在堆放或处理过程中，由于日晒、雨淋、水洗等，重金属极易移动，以辐射状、漏斗状向周围土壤扩散，固体废弃物也可以通过风的传播而使污染范围扩大。2、交通运输对土壤重金属污染的影响道路两侧土壤中的污染物主要来自汽车尾气排放及汽车轮胎磨损产生的大量含重金属的有害气体和粉尘的沉降，而污染元素则主要为Pb、Cu、Zn等元素。它们一般以道路为中心成条带状分布，强度因距离公路、铁路、城市以及交通量的大小有明显的差异。5.3.3定污染源的位置首先，根据附件，使用MATLB软件得到城区区域分布图（图11）图11 城区区域分布图为简化模型，本文
24、从Cu的污染分布着手，讨论污染源位置确定的模型，其他元素依次代入模型即可。污染源使用INGLADA方法进行初值计算4设观测点坐标为由图1图8初步估计污染源多于2个污染源到观测点的距离为对于污染源， （1）对于污染源， （2）将（1）（2）式展开后相减得到线性方程： （3）其中由于土壤为均匀介质，传播速度恒定，得等式 （4）将上述表格的数据每三个一组代入（3）式得到一系列的，从而反演得源位置。由附件所得320组数据，每四个一组，共分成80个组，代入（4）式中，得到80组污染源位置，进行统计后，的波动范围为02454，的波动范围为9871570，从而确定平面污染源的区域范围为。由题目中的采样点是在
25、间距1公里左右的网格子区域内，远小于和，（3）中的系数矩阵会有奇异性，解发散，这样更加体现此模型适宜这个问题。1.2用Geiger方法对初值进行迭代校正基于Geiger法的基本思想，建立空间直角坐标系： （5）式中分别为污染源坐标校正量，为浓度校正量，为空间偏导数，为测定浓度的偏差，可写成矢量形式：，（其中，） （6）（6）式中系数矩阵为矩阵： （7）将Inglada方法计算的初值代入（5）式求解，再用求解结果校正，如此反复迭代校正，直至满足一下两个条件。，为任一小量第次迭代，为任一小量（6）的求解方法有很多，比如牛顿-高斯算法5，但对于此问题过于复杂，又比如为提高稳定性，可采用中心化、定标化
26、、阻尼最小二乘法，但对于本题求解时，注意到遍历的数据比较多，矩阵较大，所以本文采用共轭梯度法求解。为标准评价土壤环境的污染程度，依据国家环保局颁布的土壤环境质量标准，土壤中各元素环境质量评价分级含量范围采用表8分级值。6表8 土壤金属污染评价标准（标准对比法）评价标准 As Cd Cr Cu Hg Ni Pb Zn 一类 15 0.2 90 35 0.15 40 35 100 二类 15-25 0.2-0.6 90-250 35-100 0.15-1.0 40-60 35-350 100-300 三类 25-40 0.6-1.0 250-300 100-400 1.0-1.5 60-200 3
27、50-500 300-500 超三类 40 1.0 300 400 1.5 200 500 500 5.4问题四5.4.1分析优缺点模型三中，我们大胆引进地震源的位置确定，采用的Geiger定位法是目前比较准确先进的定位方法。但是计算比较复杂。5.4.2收集信息基于以上问题的分析，除了题目给出的四个变量：样方横坐标，样方纵坐标，海拔和采样点的浓度值，我们将条件理想化，建立了以上三个模型，为了更好的研究城市地质环境的演变模式，还要丰富很多信息和相关专业知识。我们知道，更加全面客观的评价某地区的土壤污染程度，还需要掌握该区域土壤特质即土壤的物理结构和化学结构，包括土壤由固相、液相、气相和PH值以及
28、重金属元素的毒理性。在模型二中我们采用主成分分析法对主要原因的分析，还应该收集该区域的地形地貌，每个区域的地理分布状况，土壤特性以及气候影响。当然也不能忽视每个区域日常活动状况，其中包括动物植物微生物。在确定污染源的位置时，除了模型三所考虑的因素以外，重金属元素与化合物的相关性也应该纳入模型的计算，此外还应该考虑污染物与时间的关系，随着时间的推移，污染源的位置将随之变化，这样，静态模型就变为动态模型。5.4.3建立模型我们会基于上述城市信息，采用城市的地质环境的数据，运用时间序列分析中的自回归方法建立城市地质综合环境质量演变差分趋势预测模型。当预测对象随着时间呈现某种水平、上升或下降的趋势，并
29、且无明显的季节波动，又能找到一条合适的曲线反映这种变化趋势时，就可用时间为自变量，时序数值y为因变量，建立趋势模型：；正确识别并选择趋势模型是首要工作，基本方法有两种：图形识别法、差分计算7。在这里我们使用差分法，所谓差分，是变量的微小变化，根据历史数据计算差分把数据修匀，将非平稳时间序列转换为平稳序列，通过该平稳序列的表现，来发现该套用那个模型。假设时间序列为一阶向后差分为：二阶向后差分为：三阶向后差分为：多阶以此类推再根据差分特性确定使用的模型。表（9）表9 差分特性对应的模型差分特性使用模型一阶差分相等或大致相等一次线性模型二阶差分相等或大致相等二次线性模型三阶差分相等或大致相等三次线性
30、模型一阶差分比率相等或大致相等指数曲线模型一阶差分的一阶比率相等或大致相等修正指数曲线模型其中，多次线性模型一般形式为：，预测目标变量，为时间，为待定系数。当时，为线性模型；，为二次抛物线模型；，为三次多项式趋势模型。估计标准误差SE:，其中，为时间序列的长度，为限制条件，即利用最小二乘法时的所列方程数，可以简单理解为直线模型为2，抛物线模型为3，三次曲线为4，最高次数加1。指数曲线预测模型一般形式为：，为时间序列的初期水平，必须大于0。为时间序列的平均发展速度，可以为正，负值。修正指数曲线预测模型后，得，其中为时间，一阶差分的环比为常数，待定系数的求法采用“三段和值法”进行。六、模型评价与推
31、广6.1模型的优点在问题一中，我们采用目前学术界最常用的单因子指数法可以判断出环境中的主要污染因子，是综合评价指数的基础。为了更加全面合理地进行污染程度分析，我们采用尼梅罗综合指数法进行综合评价。它不仅考虑了各重金属对土壤作用的平均水平，更重要的是突出了较高浓度污染物对土壤环境质量的影响，因此，高浓度污染物对评价结果具有较大影响。这样突出了污染指数最大的污染物对环境质量的影响和作用，而且此方法操作简单，容易实现，给问题分析和结论提供了很好的依据。模型二中我们利用主成分分析法 进行精确计算，得出主成分。模型三中我们大胆引进地震源的位置确定，采用的Geiger定位法是目前比较准确先进的定位方法。模
32、型四中我们发散思维，综合考虑了不同侧面，对问题有一个全面客观的预测。6.2模型的缺点单因子指数法于尼梅罗指数法只能反映污染的程度而难于反映污染的质变特征8。模型二中的主成分分析法沿用了模型一中的结论，这样对模型一的依赖性就增强了，这会对模型造成使用的局限。参考文献1 于小俸 唐嶷 胡玉，壶瓶山自然保护区土壤重金属含量的空间分布及污染评价，中国农学通报，26(11)：358-361,2010。2 闫欣荣，修正的内梅罗指数法及其在城市地下饮用水源地水质评价中的应用，地下水，32(1):6-7,2010。3 邵学新 吴 明 蒋科毅，土壤重金属污染来源及其解析研究进展，4 牟磊育 赵仲和 张伟，用INGLADA与GEIGER方法实现近震精定位，中国地震，22(3)：294-302,2006。5 林峰 李庶林 薛云亮 徐宏斌，基于不同初值的微震源定位方法，岩石力学与工程学报，29(5)：996-1002,2010。6 杨贤智，李景锟，廖延梅，环境管理学，北京：高等教育出版社，1990。7 姜启源 谢金星 叶俊，数学模型，北京：高等教育出版社，2003。8 徐燕 李淑芹 郭书海 李凤梅 刘婉婷，土壤重金属污染评价方法的比较，安徽农业科学，36(11)：4615-4617,2008。19
