1、12011 高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料） ，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从 A/B/C/D 中选择一项填写）： A 我们的参赛报名号为（如
2、果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 北方民族大学 参赛队员 (打印并签名) ：1. 林秀萍 2. 陆彬 3. 任柳 指导教师或指导教师组负责人 (打印并签名)： 日期：2011 年 9 月 12 日赛区评阅编号（由赛区组委会评阅前进行编号）：22011 高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：3城市表层土壤重金属污染分析摘要城市土壤是城市环境的重要组成部分，对其的研究已日益成为人们关注
3、的焦点，因此建立相应的模型进行描述。针对具体问题，给出相应的解法和结果。对于问题 1，首先是定性分析，通过绘制曲面图，等高线图，箱式图，折线图，伪彩色图和散点图，直观的给出采样地点浓度在空间的分布，并预言镍元素在山区有较多的矿藏。其次定量分析，通过对每种元素平均值，方差和背景的对比，描述不同区域内重金属元素的统计特征。由于每种元素的毒性不同，采用加权平均数的方法引入污染程度的指标，得到山区的污染程度最低，和背景相差甚微。而其他的区域则有不同程度的污染。其中工业区污染最为严重，而交通区各重金属元素浓度的值都很突出，这是值得我们深思的。对于问题 2，首先总体概括土壤污染的原因，然后对数据进行相关性
4、分析，最后对具体的元素分别阐述。从元素的性质得到可能的污染源，从应用范围得到浓度如此分布的原因。针对问题 3 从热力学和动力学的角度分析，根据质量作用定律，建立带有统计意义的微分方程模型。这个模型是一个定态模型（与时间无关的平衡） 。用拟合的方法找到污染源。模型引入传播系数，它是描述污染源污染能力的重要指标，描述了在相同距离内受污染源影响的大小，同时引入半衰半径，描述污染浓度降低到一半时候的距离，也能够描述突发事件后人们的安全区的距离。并通过铜的计算得到与实际符合很好的结果。针对问题 4 由于没有具体数据，因此建立抽象函数的模型。根据地质学的原理，引入时间项，时间与地层的深度有关系，建立采样深
5、度和距离的二元函数。然后在加上突发事件的影响进行修正。这个抽象函数的功能在理论上能确定污染源的位置和发展情况，污染源周围的影响随时间变化。针对每一问都给出了相应的 matlab 运算程序，使得模型的实用性增加关键词：重金属浓度 微分方程 统计 半衰半径4摘要 3一、问题重述： 5二、模型假设： 5三、符号说明 5四、问题求解 5（1）问题一求解： 5定性分析 5定量分析计算 9（2）问题二求解： 12总体描述 12具体分析 12（3）问题三求解 13模型的建立 13模型求解 14（4）问题四求解 14模型优缺点 14模型推广 14算法 15关键程序指令 15参考文献： .165一、问题重述：按
6、照功能划分，城市可分为生活区、工业区、山区、主干道路区及公园绿化区等，分为记为 1 类区、2 类区、5 类区，不同的区域环境受人类活动影响的程度不同。现要求你们通过数学建模来完成以下任务：(1) 给出 8 种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、模型假设：1）不考虑大规模的地质运动；2）不考虑河流的影响；3）假设重金属
7、的传播途径是向四周扩散；三、符号说明：权重系数；iw：对应重金属元素在每个功能区的平均含量；ix：污染源浓度；0c：第 点重金属的浓度；i：传播系数；k：第 点离污染源的距离；ir：半衰半径.2/1：是突发事件带来的浓度变化j：是深度s四、问题求解（1）问题一求解：定性分析根据题目给出的不同地理坐标处的重金属浓度数据，要分析 8 种元素的空间分布，最直观的方式就是用图形。用 MATLAB 绘制不同种类的图形，从不同的角度来表现各种金属的浓度分布。1）各元素在全体平面的分布曲面图6（运行附件 9 程序 zong.m）图 1Hg 元素的全体分布 图 2 As 元素的全体分布图 3 Cd 在全体的分
8、布 图 4Cr 在全体的分布图 5 Zn 在全体的分布 图 6 Cu 在全体的分布7图 7 Ni 在全体的分布 图 8 Pb 在全体的分布图 18，是 8 种元素在全体采样平面的浓度分布。 （全部 8 张图见附件 1） 。图中底平面代表地理坐标平面，高度代表了浓度。很直观的得到土壤中不同重金属在同一地域内的空间分布在峰的数量，位置，高度方面都有很大差异。并能初步判断该金属的污染程度和位置。这是因为土壤成因和类型不同以及人类活动影响必然导致土壤重金属分布差异。2）等高线图曲面图由于遮挡，会有些部分看不见，因此作各元素在全体分布的等高线图，来进行信息的补充。 （由附件 9 程序 zong.m 同时
9、得到）图 9 Hg 在全体的等高线图 图 10As 在全体的等高线图以图 9,10 为代表说明等高线图补充了曲面图由于视角带来的信息缺失。按照红橙黄绿蓝的顺序浓度逐渐减小，增加了各种元素的浓度情况的直观性。 （全8 张图见附件 2）为了方便的显示不同元素在不同的功能区的差异，按元素和地区不同画出曲面图和等高线图（运行附件 9 程序 fen.m）8图 11 Cd 在工业区的浓度曲面 图 12 Cd 在工业区浓度等高线图以图 11-12 为代表 得出的结论和第一张图相差不多，不予赘述， （全部582=80 张图见附件 3、4） 3）箱式图由于不同功能区中人类活动不同，不同元素的浓度也会表现不同。因
10、此按区域类型对每种重金属做箱式图。 （运行附件 9 程序 tu.m）用 Cd-区域箱式图作代表分析从整体情况上看，镉元素在不同功能区浓度的均值，数据分布范围，集中程度，和个别突出点的数量、大小都有差异。同时图中标注了背景线，很直观的看到生活区，工业区，交通区和公园绿地区都有较高的分布浓度。交通区和工业区的平均分布相差不多，但是在交通区有很高的点，这可能与交通污染使之浓度升高有关。而山地的平均浓度较低，和背景平均相差甚微，从最值上看也比其他地区小很多，说明污染比较小。 （全部 8 张图见附件 5）图 11 Cd 在不同区域的箱式图 图 12 不同金属在不同区域的浓度对比4）折线图为了大致反映不同
11、元素在不同功能区的分布差异因此画出以元素为区分的，不同功能区的浓度折线图，同一颜色的线，代表同一功能区。横坐标表示不同元素编号。交点代表对应功能区域在特定元素浓度。另外，Cd 元素和 Hg 元素很小是因为原始数据量纲的差异。 （运行附件 9 程序 tu.m 同时得到）95）伪彩色图伪彩色图反映了连续的浓度变化，补充了等高线图的缺失。(运行附件 9 程序 fen.m 同时得到)（全部图片见附件 6、7）图 13 Pb 的浓度在 xy 平面伪彩色图 图 14 工业区 As 的浓度随海拔分布6）散点图浓度分布随高度的散点图，形式相差不大，以工业区 As 元素浓度和和海拔关系。( 运行附件 9 程序
12、fen.m 同时得到)（全部图片见附件 8）得到结论：浓度大致随海拔递减。而且数据密度也变得稀少，说明海拔越高，人类活动越弱，该地区受影响越小。不同功能区在平面的分布散点图反映了采样点在平面内的分布状态：均匀程度、密集程度。图中左下角区域点比较密集，左上角比较稀疏，从颜色上区分了不同区域，右上角为山地，左下角为工业区、生活区。图 14 不同功能区在平面的分布图定量分析计算1）首先统计各种数据用 MATLAB 对数据进行处理（运行附件 9 程序 tu.m 同时得到） ，不同区域的重金属平均含量如表 110表 1 不同区域的重金属平均含量As (g/g)Cd (ng/g)Cr (g/g)Cu (g
13、/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)生活区 6.2705 289.96 69.01849.403 93.041 18.34269.106 237.01工业区 7.2514 393.11 53.409127.54 642.36 19.81293.041 277.93山地区 4.0441 152.32 38.96 17.317 40.956 15.45436.556 73.294交通区 5.708 360.01 58.05462.215 446.82 17.61763.534 242.85公园绿地区6.2637 280.54 43.63630.192 114.9
14、9 15.29 60.709 154.24背景值 3.6 130 31 13.2 35 12.3 31 69由表 1 的平均值可以得到不同的功能区的不同金属元素的平均浓度，反应了这类型区域被不同金属的污染水平。其中同一种元素在工业区的浓度普遍偏高，山地区的普遍偏低表 2 不同地区金属含量的标准方差As (g/g)Cd (ng/g)Cr (g/g)Cu (g/g)Hg (ng/g)Ni (g/g)Pb (g/g)Zn (g/g)生活区 2.1502 183.68 107.8947.163 102.9 5.662372.325 443.64工业区 4.2443 237.58 44.002414.9
15、4 2244.1 8.370285.367 350.83山地区 1.7993 78.38 24.59510.732 27.854 10.42817.733 30.943交通区 3.2383 243.39 81.606120.22 2180.3 11.78732.527 384.78公园绿地区2.0235 235.84 14.84 22.682 224.28 4.974345.843 230.92背景值 0.9 30 9 3.6 8 3.8 6 14标准差值可以得到数据集中程度，也就是反映了有没有突出的点，方差数值越大意味着越分散。也就是说这个区域的对应元素有超乎寻常的浓度，要么非常高，要么非常
16、低。其中工业区的数值普遍偏高，结合平均值，说明有很多的浓度非常高的地区。山地区的波动不算大，也就是说明山地区的情况差不多。浓度变化不大。112）突出点的处理表 3 不同功能区平均值最值情况As Cr Cr Cu Hg Ni Pb Zn功能区平均值最大值 2.01433.02392.22649.661818.3531.61073.00134.0279平均值最大值出现的功能区编号2 2 1 2 2 2 2 2同一元素标准差最大值 4.71598.113111.988115.26280.513.101814.22831.688标准差最大值出现的功能区编号2 4 1 2 2 4 2 1分析数据可得，除
17、 Cr 以外，其它重金属的最大值出现的区域都在工业区，说明工业区的污染最为严重。这个现象完全符合实际。工业区属于制造业，在产品生产过程中一定会往外排放三废，废水经过灌溉会把重金属污染物落入土壤，废气通过扩散也会沉降到突然中，使得重金属含量升高。标准偏差的最大值多出现在工业区，说明了工业区的成分复杂。3）综合评价加权平均数评价方法，做出污染程度的定量描述。由于不同元素的化学性质不同，相同含量的重金属对环境所造成的污染是有高低之分。查阅生活饮用水卫生标准和环境质量标准污染物排放标准 1，如表 4表 4 污染物排放标准砷 mg/L 镉 mg/L 铬 mg/L 铅 mg/L 汞 mg/L 硒 mg/L
18、 铜 mg/L 锌 mg/L0.01 0.005 0.05 0.01 0.001 0.01 0.05 0.05标准的高低描述了每种金属的毒性大小。标准数值越高，毒性越小，所以在评价各区域的污染指标时就引入了一个权重系数。为消除量纲的差异，先对标准作无量纲化处理，取各元素的浓度与砷元素浓度的比。由于毒性与标准负相关，因此用标准的倒数作为权重系数。用 线 性 加 权 函 数 作 为 综 合 评 价 模 型 ， 对 每 个 功 能 区 污 染 程ijxwy81度 进 行 综 合 评 价 ， 其 中 是 权 重 系 数 ， 是 对 应 重 金 属 元 素 在 每 个 功 能 区i i的 平 均 含 量
19、 。 计 算 得生活区污染严重程度 y=166.3155工业区污染严重程度 y=219.0882山地区污染严重程度 y=82.68223交通区污染严重程度 y=164.6723公园绿地区污染严重程度 y= 129.587背景的污染程度 y=70.15数值越高，污染越严重。8 种重金属在不同的功能区具有不同程度的污染，综合污染程度工业区生活区 交通区 公园绿地区 山地区。而山地区的指标非常接近背景程度值。说明由于山地由于交通不便，人类活动少，不易产生大的生活聚居区，不易进行工业生产。12（2）问题二求解：总体描述城市土壤重金属的主要来源有 4 种 2，分别为:(1) 交通污染：汽车尾气排放、汽车
20、汽油、发动机、轮胎、润滑油和镀金部分都能燃烧或磨损而释放出 Pb、Cd、Cu 、Zn 等重金属进入周围土壤中；(2) 工矿业污染：矿厂冶炼加工、电镀、塑料、电池和化工都导致Zn、Pb 的急剧增加；(3) 燃料；(4) 生活垃圾。具体分析针对不同的元素的用途不同功能区污染原因。1）砷的化合物还用于制造农药、防腐剂、染料和医药等。砷作合金添加剂生产铅制弹丸、印刷合金、黄铜（冷凝器用） 、蓄电池栅板、耐磨合金、高强结构钢及耐蚀钢等. 3因此砷元素在不同功能区箱式图才会出现在工业区和交通区平均值和突出点比较高，生活区的平均含量不低， ，而且方差不大，说明砷元素在在工业区有专门的需求。2）汞元素用于汞化
21、合物的合成，或作为催化剂、颜料、涂料等；有的还作为药物，口服、过量吸入其粉尘及皮肤涂布时均可引起中毒。此外，雷汞(Hg(ONC)2.1/2H2O)用于制造雷管等 .4因此汞的箱式图表现出生活区和公园比较低。工业区有非常高的点说明可能有工厂生产汞或者以汞为原料。3）世界上的锌全部消费中大约有一半用于镀锌，约 10%用于黄铜和青铜，不到 10%用于锌基合金，约 7.5%用于化学制品，约 13%用于制造干电池，以锌饼、锌板形式出现。 5因此生活区才有很高的点出现，可能是电池等的回收点，而工业区有高值则是因为消耗对应的合金造成扩散所致。4）铬元素主要用于合金(如供铬钢用的铬铁合金或如镍铬电阻合金)及电
22、镀(如汽车减震器和内部装璜以及切削工具)。 6铬在工业区只有一个非常高的点，说明可能在工业区有一个专门消耗它的点。5）镉元素作为合金元素能配成很多合金，镍镉和银镉电池具有体积小、容量大等优点。镉的氧化电位高，广泛用于电镀，铁、钢、铜的保护。 7因此才表现出生活区的平均值虽然低，但是仍然有少数点的浓度几乎和工业区相同。因为镉电池消耗比较大。而作合金时就会大规模应用在工业生产和生活消耗。6）铜，电气，轻工，机械制造建筑，国防。在电气，电子工业中应用最广，用量最大，占总消费的一半以上。如电缆，电机等。 8因此，铜的箱式图也表现得生活区和工业区的平均值，和方差，以及个别突出点都相差不大，而且工业区有一
23、个非常突出的点，应该是铜的冶炼或者是加工工厂。7）铅，重要用途是蓄电池。常用作颜料，铅白保险丝。 9因此铅在除了山地区的浓度都不低，是因为颜料，电池等无论在生活区，工业区，交通区都普遍使用。8）镍主要用来制造不锈钢，和其他抗腐蚀合金。也用作加氢催化剂。 10因此镍的箱式图在大部分区域突出的点很少，因为常常和钢伴生，山地的平均13值和背景一样，但是却有非常突出的点，因此推断这个地区可能有很大的镍矿。总结上述 8 种元素的图，发现交通区总有非常高的值，说明这个区域可能是这些重金属产品的集散地或者是中转站。工业区有不同的表现说明有很多类型的工厂来生产或者消耗对应元素。而生活区的部分则说明该元素工业产
24、品的需求程度。表 5 重金属在总区域两两泊松相关系数值As Cd Cr Cu Hg Ni Pb Zn As 1 .Cd .255 1 .Cr .189 .352 1 .Cu .160 .397 .532 1 .Hg 0.064 .265 0.103 .417 1 .Ni .317 .329 .716 .495 0.103 1 .Pb .290 .660 .383 .520 .298 .307 1 .Zn .247 .431 .424 .387 .196 .436 .494 1表 5 使用 SPSS16.0 中的 Pearson 相关系数分析法，分析元素相关性，有几种元素，比如 Ni 和 Cr
25、的相关性很高，说明他们的伴生关系。因为都大量用于钢材的合金和电镀。同时也说明了上述污染部分原因。(具体区域相关性表见附件 10)（3）问题三求解模型的建立金属废弃物的扩散以废弃堆为中心向四周扩散，我们找出在一定范围内重金属浓度最高的点，污染源就应该在附近。根据热力学第二定律，浓度梯度越大，传播动力越大，传播越快。浓度越大，传播越远，这些特点质量作用定律11描述相似。而相对论中时间空间平权，因此空间 等价时间 。建立质量作用rt定律的微分方程模型。观察等高线图，各个线间距大致不等，因此不是线性衰减。最简单情况一级动力方程：（1）iikcdr对（1）式两边积分化简得：（2）ii0ln将两坐标的距离
26、 代入（2）式得：)(ir（3）20200 )()(lnyxkciii 式（3）中的 、 、 、 可以通过数据拟合得到0xy对于式子的分析， 代表了污染源浓度，污染源浓度越大，对外影响越大。c值得大小则反映了衰减特性。也就是元素在传播时候的阻力和动力效果相抵k14消之后的效果。特殊情况时候， ，即半衰半径，也就是在理论上，沿污染源附2/0cikr/2ln/1近的径向，每经过一个半衰半径，浓度变成一半。这个值也可以描述污染源的对外扩张能力。 越大说明越不容易衰减，污染影响越广，要经过很远，才能不应向人们的正常生产生活。因此这个指标也可以用来描述管理者对这个污染源的管理好坏。模型求解运行附件 9
27、程序 fitfunc.m 同时得到得到 Cu 元素的污染源大约在（2210.8 ，3146.1） ，中心点的浓度理论估计是 2.7672107 微克每克。理论浓度分布图为图 15 ，和 Cu 真实采样浓度等高线图 16 相比较，峰的尖锐程度和位置吻合很好。 图 15 Cu 元素理论浓度分布 图 16 Cu 在 xy 平面的浓度等高线得到衰减半径 42.663 米，说明这个污染源衰减比较快，也就是这个源的管理者控制管理的很好。（4）问题四求解模型优缺点由于模型是使用质量作用定律，这个定律是基于大量物质的统计结论，因此，本模型也是大量的统计规律，虽然不能准确描述浓度与坐标的关系，但是能描述该地点最
28、可能的状态。上一模型不仅给出了浓度与污染源距离的关系，还给出了该污染源污染能力的大小的衡量指标和现有区域内的浓度最高的污染源的污染能力，方便了描述和监控。同时给出了半衰半径，如果遇到土壤严重污染的地区，经过一个半衰半径，浓度衰减到原来的一半，大致判断要撤离多远，污染源的严重性才不会影响正常的生活该模型通过调节区域大小，来确定当前区域浓度最高的污染源和其特征，因此可以多次调节区域来对比不同的污染源的污染能力。该模型是定态的平衡模型，假定了扩散平衡，也就是时间无限长的情况，由于生产力的发展，人类改变环境的能力越来越强，等到达到平衡是不现实的，因此还要加入时间项，来描述随时间的变迁，以便做出更好的预
29、测模型推广由于在地质学上深度表现为时间的累积，同一地点的不同深度的浓度，则15表明了该地点污染物浓度随时间变化的情况，并记录了相应的历史变迁。由于不同的时间下可能有突发事件，比如天灾人祸，迁徙动土等对浓度产生影响，因此加入突发事件项。 )(jckrs方程分析：当浓度不随深度变化，也就是不随时间变化，方程退化成第三问的形式。当浓度不随距离变化，这样就能预言了在短时间内的污染物浓度变化。由于没有数据，只能给出相应的算法和相应的部分 matlab 指令 12。算法找到同一地点的深度和浓度的数据画出特定地点（即同一横坐标和纵坐标下的浓度和深度）图形。观察不同坐标图形是否相似，进行深度的变形，使之按照习
30、惯描述浓度和时间关系用函数式拟合关系式组合上一问的式子，得到浓度在空间和时间上的分布关键程序指令a) 按照横坐标和纵坐标排序：aa=sortrows(a,x,y)b) 提取第一关键字的坐标：b=a(:,x)c) 查看有几种该坐标；c=unique(b)d) 以上一步的元素为索引，提取特定坐标的数据行：cc=find(aa(x)=c(i)e) 提取上述特定坐标的对应的所有值:h=aa(cc,:)f) 提取特定横坐标下，第二关键字的坐标：b=h(:,y)g) 查看有几种该坐标；c1=unique(b)h) 以上一步的元素为索引，提取特定坐标的数据行：cc1=find(h(x)=c1(i)i) 提取
31、上述特定坐标的对应的所有值:d=h(cc1,:)j) 深度变换，由于深度越大，时间越久远，因此将深度变换成高度。即以最深的地方作为 0 基准面，相上的高度，gaodu=max(aa(:,shendu)- d(:,shendu)k) 为了直观，画出单一坐标下的浓度和深度的图下形：scatter(gaodu，d(:,nongdu)拟合浓度和深度的关系：比如用多项式 polyfit(gaodu, d(:,nongdu)16参考文献：1百度百科, 生活饮用水卫生标准, http:/ 年 9 月 11 日；2肖锦华,中国城市土壤重金属污染研究进展及治理对策,环境科学与管理,第 34 卷第 4 期:25-
32、28,2009 年 4 月;3百度名片, 砷, http:/ 年 9 月 11 日；4百度名片, 汞, http:/ 年 9 月 11 日；5百度名片, 锌, http:/ 年 9 月 11 日；6百度词典, 铬, http:/ 年 9 月11 日；7百度名片, 镉, http:/ 年 9 月 11 日；8百度名片, 铜, http:/ 年 9 月 11 日；9百度百科, 铅, http:/ 年 9 月 11 日；10百度名片 ,镍, http:/ 2011 年 9 月 11 日.11胡忠鲠等 .现代化学基础（第二版）.北京:高等教育出版社,2005.12登薇.matlab 函数速查手册（修订版） 。北京:人民邮电出版社， 2010
