2011年全国大学生数学建模竞赛A题论文

城市表层土壤重金属污染分析

摘要

土壤作为城市环境的重要组成部分，不仅提供人类生存所需的各种营养物质，而且接受来自工业和生活废水、固体废物、农药化肥、及大气降尘等物质的污染．很容易导致金属元素的蓄积，从而造成土壤重金属的污染．本文讨论了某城市表层土壤重金属污染的空间分布分布状况，地区污染程度，以及污染传播特征，污染源等，建立了相应的几何与数学模型或算法，得到了较好的结果，为防治城市表层重金属污染，保护和提高土壤资源和生态环境，提供参考．

对于问题一：通过给定数据的相关分析，不考虑地形高低对污染浓度变化的影响，用Matlab软件编程绘制个重金属元素污染浓度空间分布三维网格图和二维等高线图，综合研究该城市各功能区的空间分布以及污染程度分布．建立了Muller地积累指数模分析模型：，确定污染程度水平分级标准，通过统计计算，分析了各重金属在不同功能区的污染状况及程度．结论是：主干道路区和工业区的重金属元素的污染最严重，其他次之．

对于问题二，为说明重金属元素污染的主要原因，采用单因子指数模型和内梅罗综合指数模型进行综合指标评价分析，结合问题一中统计数据进行综合分析，得到个重金属元素在各功能区及城区的综合污染程度指标．污染最严重的功能区是主干道路区，其次按照污染程度从大小的顺序依次为：工业区、生活区、公园绿地区、山区．主干道路区土壤表层重金属元素含量很高，且种类多．根据地区的差异性和元素的特殊性，分析出重金属污染Hg和Cu污染是最严重的污染源，且污染最严重的地区在主干道路区和工业区．这些污染主要由于含铅汽油的燃烧、汽车轮胎磨损产生的含锌粉尘、工业污水的排放、生活废水的排放、化肥农药的多度使用、金属矿山的开采．详细情况见正文．

对于问题三，为了找出该城区的污染源，在分析出重金属元素的主要传播特征之后，考虑大气空间传播情况，建立了微分方程模型，通过模型求解分析，用其等效的向内（向污染浓度较高的方向搜索）搜索算法，计算确定了重金属元素主要污染源的位置，其中As较严重的中心污染源坐标分别为：（5291，7349，10）、（12696,3024，27）、（18134、10046、41）、（17814,10707,64）、（27700,11609,165）．这五个污染源主要分布在主干道路区．（5291，5739，10），（12696,3024，27），（17814,10707，64）分布在工业区，其它两种污染源分布在生活区．其余元素的中心污染源见正文．

对于问题四，需对前面所建立的模型进行分析与评价并进行模型的优化，在详细分析了前三个问题的求解模型及过程之后，评价出所建立模型的优缺点．在问题三中，重金属元素除了在大气中传播以外，还通过水土流动传播．另外，前几个模型都是静态的，但污染物传播的过程与时间有关，是一个动态的过程．最后建立了一个扩算方程模型进行优化，能够为更好的研究城市地理环境的演变模式做贡献．

关键词：重金属污染 地积累指数模型 单因子指数模型 内梅罗综合污染指数 微分方程模型

一、问题重述

1．1 问题背景

随着工业发展和城市化进程的加剧，通过交通运输、工业排放、市政建设和大气沉降等造成城市重金属污染越来越严重．对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究城市不同功能区表层土壤重金属污染特征和污染空间分布性，以便更好的研究人类活动影响下城市地质环境的演变模式．

本文 就如何应用查证获得的海量数据资料展开城市环境质量评价，研究地质环境的演变模式建立数学模型．附录1列出了采样点的位置、海拔高度及其所属功能区等信息，附录2列出了8种列出了8种主要重金属元素在采样点处的浓度，附录3列出了8种主要重金属元素的背景值．

1．2 需要解决的问题有

(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度．

(2) 通过数据分析，说明重金属污染的主要原因．

(3) 分析重金属污染物的传播特征，由此建立数学模型，确定污染源的位置．

(4) 分析所建立模型的优缺点，为更好的研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？

二、问题分析

该题目一方面通过GPS记录了该城市大量样本点的位置以及所属功能区，再应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据，通过这两个表的数据 就大致可以提取出一些对于解决问题的重要信息，另一方面，题目给出了自然区各样本点的重金属元素的背景值，作为重金属污染情况的指标．对于 分析研究各个样本点的污染程度至关重要． 利用Matlab软件进行三维网格图和等高线图的制作并结合相关的数据统计分析， 可以分析该城区不同区域重金属的污染程度．后面 利用地积累指数法和内梅尔综合评价指数对数据进行处理，分析污染严重的功能区和重金属．结合图形的分析以及模型的建立 综合分析重金属污染物的传播特征．接着对模型进行一定的优化处理，使得处理的结果更准确．

三、模型假设

1、假设题目所给的数据合理正确．

2、该区域的划分是稳定的，不会出现大的变动．

3、不考虑观测误差、随机误差和其他外在因素所产生的误差．

4、重金属在大气中无穷空间扩散，不受风的影响，其扩散服从热传导定律．

5、重金属污染程度连续变化，大气中重金属元素浓度连续变化．

6、界限不明显区域有扩大、缩小、消失的过程，穿过大气进入仪器的重金属含量只有浓度大小之分，浓度大小由仪器灵敏度确定．

四、变量与符号说明

五、模型建立与求解

针对问题一，首先想到的是用Matlab软件编程，进行三维网格图、三维曲面图、等高线图和散点图的制作．

5．1 问题(1)的分析、模型建立与求解：

5．1．1 问题（1）的分析

对于问题一，首先来分析一下， 要给出8种主要重金属在该城区的空间分布， 就必须确定每个重金属元素与他们所对应的地区之间的联系．刚好题目给出了每个样本点的各元素浓度，那么 是不是可以将每种重金属元素含量浓度含量与该目标点所在的功能区建立联系？由此 想到利用Matlab软件画出每种元素在该城区的三维曲面和空间曲面图．同时 在分析不同区域重金属的污染程度时，考虑到这个污染程度是否可以量化．并且是否能够建立一种模型将这种指标量化．这道问题还要求考虑每个功能区的污染程度， 知道每个功能区的每种重金属污染程度又是不一样的．那 通过什么指标来判断每个功能区的污染程度大小，这也是 为什么用权重作为评价每个功能区的污染程度的指标．

5．1．2 问题（1）的模型建立

该城区受这八种重金属元素As、Cd、Cr、Cu、Hg、Ni、Pb、Zn污染程度不一样．题目提供了每种重金属元素的背景值，那么 怎么样利用这些背景值和每种元素相关的浓度确定不同区域重金属的污染程度？所以 需要找出一种方法来准确的分析该城区内不同区域重金属的污染程度，并且最好能够量化．建立8种主要重金属元素在该城区的空间分布模型如下：引入了一种用于研究沉积物及其他物质中重金属污染程度的该区内不同地域重金属的污染程度的定量指标——地积累指数又称Muller指数法，Muller指数法表达式为:

式中Fn表示污染物重金属元素n的浓度；Bn 表示第n种重金属元素的背景值上限，C为考虑各地岩石差异可能会引起背景值的变动而取得一系列系数(一般取值为1．5)，用来表征沉积特征、岩石地质及其他影响．Muller地积累指数评价和分级标准分级标准具体详见表1

表1：地积累指数分级标准

该方法指标主要是通过每种重金属元素测得的实际浓度以及他们的相关背景值，计算出每种元素的地积累指数．然后根据上面这张表 就可以判断出每种元素的污染级别，这样就可以对每种元素的污染情况进行分析．然后 再利用Matlab软件对题目所给数据进行处理，画出相应的网格曲面图和等高线曲线图．这里需要对Matlab进行编程，首先利用每个样本点的横坐标、纵坐标、海拔高度建立等高线图，程序语句见附录一．通过该图，可以直观的看到该城区各功能区的空间分布．但是这张图不能反映出8种主要元素在城区的污染情况， 需要借助于各种主要元素的浓度．所以 需要再建立一张等高曲线图以及相应的网格曲面图，将主要元素的浓度作为第三坐标，命令语句见附录一．

5．1．3 问题（1）的求解过程

首先通过Matlab软件，调用每个样本点的位置相关数据．就是以海拔为第三坐标，并且对每个功能区进行颜色区分，画出该城区每个功能区的二维等高线图．最后把每个样本点显示在图上．得到如下这张图：

图一：重金属As空间二维等高线分布图

这张图只反映出了该城区各功能区的空间分布，还不能看出每种重金属污染的情况． 将每种重金属元素的浓度在图上反应出来，做出该城区重金属污染的二维等高线图．具体程序语句见附录二，得到如下这张图：

图二：重金属As分布平面图

同时为了对应这张As含量分布平面图， 也画出了三维网格曲面图（图三）．

图三：重金属As含量分布的空间三维图

从空间三维图三中可以看到，有一处的波峰很高说明该处污染情况很严重，有二处处于波峰说明污染情况比较严重的主要有二处，还有一处面积比较广且所处高度稍微低一点这表明该处所受污染情况相对严重且污染的范围较广；同样分析二维等高线图二，图中有一处等高线之间的间距越来越密集且颜色很深表明该处受污染情况很严重，有二处等高线比较密集颜色相对较深表明这二处的污染情况相对严重，还有一处等高线间的距离较密集但是所包围的面积较广说明该处的污染也较严重且污染的面积很广．再结合前面的数据他们中心污染源的坐标分别为：（5291，5739），（12696,3024），（17814,10707）．都是分布在工业区，还有一处污染级别不是特别严重,但是在该处存在着污染源，此处刚好是山林密集区．通过观察图三，会发现刚好有三个点处于波峰，还有个点波峰稍微偏低，但还是能很直观的看出来．

再来看一下，Cd这种重金属的城区各功能区的二维等高线图，分布平面图，空间分布图（图四、图五）：

图四：重金属Cd空间二维等高线分布图

图五：重金属Cd含量空间分布平面图

以及相应的三维网格曲面图（图六）：

图六：重金属Cd含量空间分布图

从空间分布图六中可以看到，污染情况比较严重且面积比较广的主要有一处，还有五处污染也相对严重．以及几处小的污染；同样从二维分布图五可以看出等高曲线所谓面积有一处颜色很深，说明该区域污染情况很严重，同时也观察到又五处等高曲线所围的面积颜色比较深，这说明了这五处区域污染情况相对严重，很明显的是有一处等高曲线所围成的面积比较广且颜色较深，表明了该区域有一处污染情况较严重且污染面积比较广，由此可见不管是从二维还是三维图形进行分析的结果是相吻合的．再结合前面的数据它金属Cd中心污染源的坐标为：（22304,10527）．分布在主干道路区，还有一处污染级别不是特别严重．再观察图三，会发现刚好有三个点处于波峰．

如此，通过同样的方法，都能够得到对其它六种种重金属在该城区的空间分布以及污染情况的了解（参见附录三）通过观察每种元素的三维曲面图以及等高曲线图．很容易观察到，每种重金属对该城区都存在或大或小的污染．其中有些地区是存在多种重金属污染，并且污染很严重，通过观察这8张图 会发现这六种元素Cd,Cr,Cu,Hg,Ni,Pb在横坐标在[3000,4000],纵坐标在[3000,6000]这个区域内含量都非常高， 大致可以判定这段区域属于重度污染区．

下面将题目中所给的数据用excel进行分类处理，得到样本点的地积累指数．然后运用数学统计法得到各种元素污染程度数据分布表，通过这些表 就可以确定该城区内不同区域重金属的污染程度．统计该表时，是通过统计每个功能区的总样本点个数，然后通过地积累指数法分别计算出每种样本点的地积累指数，并判断他们的所在的污染级别．然后统计每种污染级别下，各功能区的污染点数占总点数的百分比也就是说的权重，通过该权重 就能够分析出每种重金属元素的污染程度大小，以及污染所波及的范围．从而得到每种重金属元素污染最严重的地区．通过Excel对数据运算，得到重金属元素As污染情况分布表：

表二：As污染程度分布数据表

下面通过同样的数据处理，得到Cd污染程度数据分布表：

表三：Cd污染程度数据分布表

其它六种元素的污染程度数据分布表见附录三．

表中数值0的意义是在该污染级别下不存在观测的样本点．这是个大样本事件，可以认为该级别污染很轻微，甚至不存在这种级别的污染．而百分比越大，就说明在该污染级别下涉及的样本点比较多，污染波及范围较广．

5．1．4 问题（1）的结果分析

5．1．4．1 As这种重金属污染情况分析

由该表 可以看出各个区域受As的污染程度，其中一类区即是生活区31．82%无污染，63．64%轻度—中度污染，4．55%为中等污染，无强污染和及严重污染的情况；二类区即是工业区38．89%不受重金属污染，52．78%受轻度—中度污染，5．56%受中等污染，2．78%受中等—强污染；三类区即是山区大多数不受污染，只有15．15%受轻度—中度污染，1．51%受中等污染；四类区即是主干路区47．83%不受污染，50．00%受轻度—中度污染，0．72%受中等污染，1．45%受中等—强污染；五类区即是公园绿地区大多数受轻度—中度污染，25．71%不受污染，2．86%受中等污染． 再结合相应的几何图形，会发现在四区存在三个很明显的污染源，在污染源附近 会看到，有很多二区的样本点．有个别一区的点，说明这种元素对一区的影响相对来说轻点．所以 由分析可知工业区受污染最严重，污染面积达到了61．11%，其次是生活区、主干道路区，生活区污染面积都达到了50%以上，也就是说这三个区有至少一半的土壤受到该元素的不同程度的污染．其余功能区受污染程度就次之．

5．1．4．2 Cd这种重金属污染情况分析

由该表可以看出各个区域受Cd的污染程度，其中一类区即是生活区29．55%无污染，54．55%轻度—中度污染，13．64%为中等污染，无强污染和及严重污染的情况；二类区即是工业区16．77%不受重金属污染，44．44%受轻度—中度污染，30．56%受中等污染，8．33%受中等—强污染；三类区即是山区大多数不受污染，只有75．76%受轻度—中度污染，21．21%受中等污染；四类区即是主干路区23．91%不受污染，44．2%受轻度—中度污染，26．09%受中等污染，5．07%受中等—强污染；五类区即是公园绿地区大多数受轻度—中度污染，48．57%不受污染，31．43%受轻度-重度污染，11．43%受中等污染，8．57%受中等-强污染．再结合相应的几何图形，会发现在四区存在三个很明显的污染源，在污染源附近 会看到，有很多二区的样本点．有个别一区的点，说明这种元素对一区的影响相对来说轻点．所以 由分析可知工业区受污染最严重，污染面积达到了61．11%，其次是生活区、主干道路区，生活区污染面积都达到了50%以上，也就是说这三个区有至少一半的土壤受到该元素的不同程度的污染．其余功能区受污染程度就次之．

5．1 这六种重金属Cr、Cu、Hg、Ni、Pb、Zn污染情况分析

由于重金属含量越多，说明该地区的重金属污染程度越严重．

Cr污染最严重的有一处，该中心污染源的坐标为：（3299，6018），所在地区为主干道路区，一定程度上波及到了生活区．一区和四区存在强-及严重污染，一区波及面积达到了52．27%，四区波及面积达到了41．3%，该元素污染最严重的就是生活区．

Cu污染最严重的有一处，该中心污染源的坐标为：（2427，3971），所在地区为生活区，一定程度上波及到了工业区和主干道路区．一区和四区存在及严重污染，一区污染波及范围达到了84．09%，四区污染波及范围达到了84．06%，该元素污染最严重的就是生活区和主干道路区．

Hg污染最严重的有一处，中心污染源的坐标为：（3299，6018），所在地区为主干道路区，一定程度上波及到了生活区．一区和四区存在及严重污染，一区污染波及范围达到了54．55%，四区污染波及范围达到了50．74%，该元素污染最严重的就是主干道路区．

Ni污染最严重的有一处，中心污染源的坐标为：（3299，6018），所在地区为主干道路区，一定程度上波及到了生活区．一区、二区和四区存在及严重污染，一区污染波及范围达到了90．91%，二区污染波及范围达到了94．44%， 四区污染波及范围达到了93．48%，该元素污染最严重的就是主干道路区和生活区．

Pb污染最严重的有二处，中心污染源的坐标为：（2383，3692）、（5062，4339），所在地区为生活区和主干道路区，一定程度上波及到了工业区．一区和四区存在及严重污染，一区污染波及范围达到了52．73%，四区污染波及范围达到了80．87%，该元素污染最严重的就是主干道路区．

Zn污染最严重的有一处，中心污染源的坐标为：（14065，10987），所在地区为主干道路区，一定程度上波及到了工业区．四区存在及严重污染，四区污染波及范围达到了67．39%，该元素污染最严重的就是主干道路区．所以，该城区不同区域重金属污染最严重的区域是主干道路区和工业区，其次是生活区、公园绿地区、山区．

5．2 问题（2）的求解：

5．2．1 问题（2）的分析

通过问题一的分析，可粗劣的判断哪几种元素污染比较大，哪个功能区污染比较严重，但是 怎么样才能具体到哪个功能区污染最严重，被污染的功能区的土壤哪种重金属污染最严重？所以，针对问题二给出的数据分析，不能简单的进行数据处理．为了使得所寻找出来的原因更有说服力，用两种方法分别进行说明和验证，还要进行综合指标评价．最后确定了最严重的污染地区以及污染最严重的相关元素，根据地区的差异性和元素的特殊性，才能说明重金属污染的主要原因．

5．2．2 数据的统计分析

首先 通过数据的处理，建立每个功能区各重金属元素的污染程度样本所占的百分比表．一功能区的相关百分比数据如下：

表四：一功能区各重金属污染程度所占百分比

在此功能区从总体来看，重金属污染程度处于中等-强污染，其中主要污染来自重金属元素Ni，另外在该区域有少数地方Cu污染及严重．

表五： 二功能区各重金属污染程度所占百分比

在该功能区重金属Hg和重金属Ni的污染极为严重，尤其是在该区域的某些地方．由此可见，在此功能区照成重金属污染的罪魁祸首为重金属元素Hg和重金属元素Ni．

通过这两张表， 会发现有些地区之所以污染严重，主要是因为个别元素污染所导致的．所以 要分析重金属污染的原因，就得分析该重金属在该功能区为什么会产生污染．其它三个功能区各重金属污染程度百分比见附录三．通过该附录表 可以看到在该功能区里，重金属污染程度较轻，污染等级集中在轻度污染及以下． 再观察功能区四，重金属污染十分严重，大多数重金属污染元素都集中在在各个功能区，但是在这个功能区，Pb污染级别比较轻，没有中度甚至以上级别的污染． 再看功能区五，从总体上分析，该地区重金属污染中等、强污染几乎没有，正因如此造成重金属污染的少数种类重金属元素就凸显出来了——Ni元素和Hg元素．纵观整体，分析所有的功能区， 很容易发现造成重金属污染的主要重金属元素，他们就是Ni元素和Hg元素．

知道前面的数据分析理由不充分，只是一个粗劣的判断．为了综合前面处理的数据，准确找出各个功能区污染的主要元素． 需要利用单因子指数法和内梅罗综合污染指数法进行综合评价．

5．2．3 单因子指数法和内梅罗综合污染指数法的建立与求解

单因子指数法是目前国内土壤重金属的单项污染指数评价方法之一，其计算公式为：，式中为重金属污染物n的环境质量指数；为重金属污染物i的实测值；为重金属污染物的评价标准．﹥1表示污染； =1或﹤1表示无污染；且值越大，则污染物越严重．为了更全面的反应各重金属对土壤的不同作用．突出高浓度重金属对环境质量的影响， 采用内梅罗综合污染指数法．其计算公式为:,式中表示重金属污染物种污染指数的最大值；表示重金属污染物中污染指数的平均值．土壤污染水平分级标准采用国家土壤环境二级标准．土壤污染综合污染指数分级标准为综合污染指数>3为重污染，2~3为中污染，1~2为轻污染，0．7~1为警戒级，≤ 0．7为安全级．

下面为了找到每种元素在该城区的综合污染指数，借助于Matlab循环计算．编写如下系列命令见附录七．运行程序结果为As综合污染指数：p=4．0093，分别运行另外几种程序，得到每种重金属元素的综合评价指标，简单结果如下表：

表六：每种重金属元素的综合污染指标

这样处理后 很容易看到哪种重金属元素污染比较严重，或者对哪个区的污染比较严重，但是 并不能从整个城区的角度上分析哪个重金属元素的污染程度大小．为此， 综合处理以上的有关综合污染指数，同样的 运用Mallab循环计算，用同样的方法编写程序命令， 得到各个区域的每种重金属元素的综合评价指标以及相应的简单结果如下：

表七：每种重金属元素的综合污染指标

5．2．4 城区重金属污染情况综合分析

知道综合污染指数越大，即P值越大，说明该重金属污染越严重．对比上面二个表 知道，Hg元素污染最严重，其次按照污染程度从大到小排序依次为：Cu Zn, Cr, Pb, Cd, Ni, As．综合前面的数据统计分析结果， 可以断定Hg元素的污染是引起重金属污染的主要原因，而Hg元素的浓度之所以变得那么大，与该元素所分布的功能区密切相关．通过观察， 发现Hg这种元素污染最严重的地区是在主干道路区．Cu这种元素污染最严重的地区是工业区．Zn这种元素污染最严重的地区是在主干道路区．Cr这种元素污染最严重的地区是在主干道路区．Pb这种元素污染最严重的地区是在生活区和工业区．Cd这种元素污染最严重的地区是在主干道路区．Ni这种元素污染最严重的地区是在主干道路区．As这种元素污染最严重的地区是在主干道路区．

通过分析每个功能区的加权平均数，发现污染最严重的功能区是主干道路区．其次按照污染程度从大到小的顺序依次排序为：工业区、生活区、公园绿地区、山区．之所以主干道路区污染最为严重就是因为该区土壤表层重金属元素含量很高，且种类多．

5．2．5 重金属污染的主要原因

根据上面的分析结果， 知道Hg这种重金属元素污染最为严重的原因就是主干道路区的含量太多．公路、铁路两侧土壤中的重金属污染，主要是Pb、Zn、Cd、Cr、Co、Cu的污染为主．它们来自于含铅汽油的燃烧，汽车轮胎磨损产生的含锌粉尘等．随着时间的推移，公路、铁路土壤重金属污染具有很强的叠加性．这也是为什么主干道路区各种重金属污染都很严重．

其次，工业区的重金属污染也很严重，主要是因为经过自然沉降和雨淋沉降进入土壤的重金属污染，主要以工矿烟囱、废物堆和公路为中心，向四周及两侧扩散；由城市—郊区—农区，随距城市的距离加大而降低，特别是城市的郊区污染较为严重．此外，大气汞的干湿沉降，也可以引起土壤中汞的含量增高．大气汞通过干湿沉降进入土壤后，被土壤中的粘土矿物和有机物的吸附或固定，富集于土壤表层，或为植物吸收而转入土壤，造成土壤汞的浓度的升高．还由于部分矿产开发中选矿、冶炼工艺水平落后，个别矿区没有环保治理设备，废水、废气排放而带来的大量废弃物的产生未经处理直接投放环境，而其中的重金属也随着自然的沉降、雨水的淋溶等途径进入土壤，进入正常循环的生态系统，造成重金属污染严重危害人们的生产生活．

生活区的重金属污染主要来自生活废水的排放，化肥农药的多度使用，化肥中品位较差的过磷酸钙和磷矿粉中含有微量的As、Cd重金属元素．含铅及有机汞的农药发挥作用的同时也为土壤重金属污染埋下了祸根，造成土壤的胶质结构改变，营养流失，对农作物的产量及品质都造成极大的不良影响．目前的饲料添加剂中也常含有高含量的Cu和Zn，这使得有机肥料中的Cu、Zn含量也明显增加并随着肥料施入农田．

公园绿地区的重金属污染主要来自城市绿地区，施用肥料如磷肥、含锌肥料等也会带来Hg的严重污染．

山区的重金属污染主要是来自金属矿山的开采、冶炼、重金属尾矿、冶炼废渣和矿渣堆放等，可以被酸溶出含重金属离子的矿山酸性废水，随着矿山排水和降雨使之带入水环境（如河流等）或直接进入土壤，都可以间接或直接地造成土壤重金属污染．

总的来说，随着城区工业化的进程越来越快，各种重金属的排放会越来越多，严重影响了人们的生活健康和环境．

5．3 问题（3）的模型建立与求解

题目要求 先根据金属污染物的传播特征，然后由此建立模型，通过模型的求解，确定污染源的位置．

5．3．1 重金属污染物的传播特征

根据问题一中建立的重金属元素分布等高线图， 能够了解城区重金属的分布规律．它们是由污染源向四周逐渐扩散的，它们的扩散方式主要有：

1、大气中的重金属主要来源于工业生产、汽车尾气排放及汽车轮胎磨损产生的大量含重金属的有害气体和粉尘等．它们主要分布在工矿的周围和公路、铁路的两侧．大气中的大多数重金属是经自然沉降和雨淋沉降进入土壤的．

2、公路、铁路两侧土壤中的重金属污染，主要是Pb、Zn、Cd、Cr、Co、Cu的污染为主．它们来自于含铅汽油的燃烧，汽车轮胎磨损产生的含锌粉尘等．它们成条带状分布，以公路、铁路为轴向两侧重金属污染强度逐渐减弱；随着时间的推移，公路、铁路土壤重金属污染具有很强的叠加性．

3、经过自然沉降和雨淋沉降进入土壤的重金属污染，主要以工矿烟囱、废物堆和公路为中心，向四周及两侧扩散；由城市—郊区—农区，随距城市的距离加大而降低，特别是城市的郊区污染较为严重．

4、大气汞的干湿沉降也可以引起土壤中汞的含量增高．大气汞通过干湿沉降进入土壤后，被土壤中的粘土矿物和有机物的吸附或固定，富集于土壤表层，或为植物吸收而转入土壤，造成土壤汞的浓度的升高．5、在河流中不同河段的重金属污染往往受污染源（矿山）控制，河流同一污染源的下段自上游到下游,由于金属元素迁移能力减弱和水体自净化能力的适度恢复，金属化学污染强度逐渐降低．

5．3．2 微分方程模型的建立与求解分析

热传导定律：单位时间通过单位法向量面积的流量与污染程度梯度成正比

1） 重金属污染浓度的变化规律：在通过的流量

内重金属元素的改变量：

曲面积分的奥氏公式

，

，

2）重金属元素污染程度的变化规律：

初始条件为．

Q：表示重金属元素的排放总量，单位浓度的点源函数

对任意的等值面是球面当增大时，C减小．

仅当，对任意点,．

3) 空气中的重金属元素含量在大气中传播的变化规律

沿沿方向的大气中重金属元素总量

沿方向的重金属元素污染程度

大气中重金属元素的总量与重金属元素污染程度成正比

，

，未进入大气中的重金属元素含量的重金属元素含量为

4） 仪器灵敏度与大气中重金属元素浓度界限

仪器灵敏度 ，当, 观测结果为浓度小，设污染源在，仪器在，观测到的浓度界限为 ~浓度小（仪器无法观测）区域边界 ．

5） 重金属元素浓度小区域边界的变化规律

，

，

，

，

对任意t，重金属元素浓度小区域边界是圆周，

重金属元素浓度小区域半径．

6） 结果分析

，（最大值），．

观测不到重金属元素边界达到最大的时刻，可以遇到重金属污染程度为零污染的时刻．

注：当，当．

5.2.3 微分方程模型的求解

从上述微分方程模型，和求解过程看某些参数 无法直接确定，但这个算法可以考虑从污染浓度等高线（污染浓度作立标）零污染边界向内收寻污染源，简化为从等高线节点向内搜索污染源，逐步减少搜索半径，最终确定污染源．

利用Matlab求解该微分方程,考虑如下算法：

1. k=0, 任取初值点（包括二维坐标、浓度）；

2. 搜索初值点附近浓度相同的样本点至少两个, %寻找等高线；

3. 拟合前述所得样本点的等高线；并向内搜索样本点，如果该样本点浓度增高，则以该样本点为初始点，重复第2步；否则

4. 继续向内搜索样本点，重复第3步；

5. 如果无法在初始样本点附近找到新的样本点，使其满足浓度增高，则该样本点即为污染源，输出该样本点；

6. 循环以上步骤，继续寻找下一个污染源．

根据该算法确定出每种重金属元素的污染源，所在位置如下：

针对As元素：污染情况比较严重区域的主要有五处，它们的编号分别为：28、178、84、83、121．再结合前面的数据他们中心污染源的坐标分别为：（5291，7349，10）（12696,3024，27）、（18134、10046、41）、（17814,10707,64）、（27700,11609,165）．

针对Cd元素：重金属Cd的污染程度严重的地区有第一区1个，第二区2个，第四区3个，第五区1个，它们的编号为16,6,8,95，22,143，且坐标分别如下：（4777,4897,8），（1647,2728,6），（2383,3692,7）（22304,10527,40），（3299,6018,40）,（2708，2295,22），（4156,2299,73）；

针对Cr元素：重金属Cr的污染程度严重的地区有第四区1个，它的编号为：22，以及它的地理位置为（3299，6018，4）；

针对Cu元素：重金属Cu的污染程度最严重的地区为第二区，第四区．它们编号为：8,22，以及它们的地理位置为（2383,3692,7），（3299,6018,4）；

针对Hg元素：重金属Hg的污染程度最严重的地区为第二区1个，，第四区为3个，它们编号为：8,9，，182，257以及它们的地理位置为（2383,3692,7），（2708,2295,22），（13694,2357,3374），（15248,9106,16）；

针对Ni元素：重金属Ni的污染程度最严重地区为第四区，第三区，它们编号为：22,135以及它们的地理位置为（3299,6018,4），（22193,12185,79）；

针对Pb元素：重金属Pb的污染程度最严重次序为第一区＞第四区＞第二区＞第四区＞第五区，它们编号为：16,6,8,20,143以及它们的地理位置为（4777,4897,8），（1321,1791,28）；

针对Zn元素：重金属Zn的污染程度最严重的地区为第四区1个．它的编号为：250以及它们的地理位置为（14065，10987，25）；

六、模型的分析与评价

6．1 针对问题四的求解

我们分别分析前面的模型优缺点，结果如下：

1 对于空间分布几何模型：该模型能够直观显示各种重金属的空间分布以及各种重金属的污染程度，及地区分布，方法简单，易操作，但该模型无法从量化的角度分析该污染情况，该几何模型熟悉相关软件．

2． 对于地累积指标模型：的模型较简单，计算方便，且无量纲化，可以用于各种元素的污染程度比较，但数据计算量比较大，而且污染程度指标的分级是否完全合理有待于 进一步实践检验．

3．对于采用单因子指数模型： ，反映了各种重金属对土壤的不同污染程度，但是单因子指数法针对的是静态的重金属元素，但事实上，重金属元素会随某些原因而发生改变，从而产生误差．

对于采用内梅罗综合污染指数模型：是综合反映8种重金属元素的总的污染程度的大小的综合指标，或者同种元素在某个地区污染程度的大小的综合指标，项对指标值影响较大，从公式看若某地区有重大污染源，那么该指标值较大，有时甚至相当大；事实上， 计算出的指标有超过200的情况．反之，如果该指标值较大和很大，若是计算的地区综合指标，说明该地区有严重的污染源；若是计算的某观测点8种元素的综合指标，说明某种元素该观察点有强的、甚至是非常强的污染源，不能具体得出是哪种重金属元素产生的污染．其缺点是不能反映地区污染分布状况． 需要地累积指数模型与该模型结合使用．

4． 对于加权平均模型：，它把综合污染指标量纲化，且消去了偶然的重大突变源，便于量化各功能区的污染程度，该方法简单，易于统计，方便 对最初模型进行验证，该方法具有较强的实施性．

5． 对于微分方程模型：虽能判别出重金属污染浓度的变化规律，但是该方法只能判别大气中重金属的传播特征，具有一定的局限．且计算量比较大．

七、模型的改进

7．1 问题四中信息的收集

在整个问题的求解过程中， 运用了微分方程模型来判别重金属污染浓度的变化规律，虽然能初步判断出重金属在大气中的传播特征 但该问题是要讨论重金属的传播特征，所以 还应考虑到重金属在地下的传播特征．为了更好地研究城市地质环境的演变模式，还应该收集的信息有：

分析土壤中重金属污染物的传播特征， 选用BCR连续提取法分析土壤中重金属的形态分析，然后判断重金属在土壤中的传播特征．BCR连续提取分析法的原理是指用各种提取剂对土壤重金属的各个形态进行提取，进而采取一定的方法测量其各个形态含量．土壤中重金属形态分析方法中共有的或比较重要的形态定义如下：

可交换态重金属：是指吸附在土、腐殖质及其他成分上的金属，对环境变化敏感，已于迁移转化，能被植物吸收；

碳酸盐结合态重金属：是指土壤中的重金属元素在碳酸盐矿物上形成的共沉淀结合态，对环境条件特别是pH值最敏感：当pH下降时，易重新释放出来而进入环境；当pH升高时，有利于碳酸盐的形成．

铁锰氧化物结合态重金属：一般是以矿物的外囊物和细分散颗粒存在，活性的铁锰氧化物比表面积大，吸附或共沉淀阴离子而成．当pH值和氧化还原电位较高时，有利于铁锰氧化物的形成，铁锰氧化物的结合态反应了人文活动对环境的污染．

有机结合态重金属:使土壤中各种有机物与土壤中的金属螯合而成的，反应水生生物活动及人类排放富含有机物的污水结果．

残渣态重金属：一般存在硅酸盐、原生和次生矿物等土壤的晶格中，是自然地质风化的结果，在自然条件下不易释放，能长期稳定在沉积物中，不易为植物吸收．主要受矿物成分及岩石风化和土壤侵蚀影响．

假如 用BCR连续提取法分析出该功能区土壤中有可交换重金属、有机结合态重金属等，那么 就可以说该功能地区有的重金属易于迁移转化，并且是通过植物进行传播的，有的是与土壤中各种有机物结合进行传播的．同理对不同的区域采用BCR连续分析出该区域的重金属形态，然后分析出该区域的重金属传播特征．

7．2 偏微分方程的建立

大气中重金属沉降、污水(泥)灌溉、废弃物堆积、金属矿山污染金属会在污染源附近形成土壤污染沉降（扩散），土壤污染在土壤中的扩散计算通常偏微分方程模型．重金属在土壤当中会因为分子扩撒、地表径流、地下径流、生物扩散等原因来传播，它可以模拟成一个扩散方程或者运用拉普拉斯模型等来解决．

假如所给的数据认为是在平衡态测得，重金属污染物从高浓度处向低浓度出扩散，与扩散同时，一部份重金属污染物由于外界因素进入地下被土壤吸收，扩散系数与吸收系数都是常数 ，但是经过一段时间，进入土壤的重金属污染物达到了平稳状态．在这种假设之下，设是t时刻在点处重金属污染物的浓度，是t时刻在点处被土壤吸收的重金属污染物的浓度．又设重金属污染物在各个方向的传播系数k,吸收系数h,于是有

（1）

又，即是重金属在该区域被土壤吸收的系数与重金属在该区域的浓度成正比，带入上式得

（2）

对于上式关于t从0到+∞积分得到

由于最后无重金属，故，又开始时（）被吸收的元素，故平衡状态在t=+∞时达到，这 ，于是由（2）得

（3）

其中是开始时的浓度分布，近似于污染源的脉冲函数．把此污染源取为坐标原点（0,0,0），则，L时污染源处含量，于是（记）

（4）

作傅里叶变换得

(5)

,

再作反变换得

（6）

其中C是可计算常数．则

=

，

于是解得

（7）

（7）

中的D为可计算常数．

用这种模型就可以进行参数估计，达到模型的优化，精确的计算出污染源的位置．

八、模型的应用与推广

微分方程模型是一类常用的数学模型，常用于描述对象的特性随时间（或空间）演变的过程，预测未来动态，研究控制手段．本文提出的关于制动试验系统的数学模型对于实际的试验是有参考作用的． 设计的控制方法简单易行，控制精度高，非常适于工程应用．

(1)对于现实世界的变化，人们关注的往往是变量之间的变化率，或变化速度、加速以及所处的位置随时间的发展规律之中的规律一般可以写成一个微分方程或方程．所以实际问题中，有大批的问题可以用微分方程来建立数学模型，涉及的领域包括物理学、化学、天文学、生物学、人口、资源等等．

参考文献：

[1] 姜启源．数学模型（M）．第三版．北京：高等教育出版社，2003．8

[2] 张德丰．MATLAB数值计算方法（M）．北京：机械工业出版社，2010．1

[3] 汪晓银，周保平．数学建模与数学实验（M）．北京：科学出版社，2010．2

[4] 刘绍贵．南昌市城郊表层土壤重金属污染特征研究（EB/OL）． http://wenku．baidu．com/view/e8e14d10f18583d0496459dd．html, 2010-04-02．

[5] 伊俊，柳云龙．上海市城郊土壤重金属空间分布及其污染评价（EB/OL）．

http://wenku．baidu．com/view/bd85958b84868762cbaed505．html，2010．

附录

附录一： 程序

1． 城区各功能区的地形二维等高线污染浓度散点图分布图程序

clc

clear

close all

D=load('坐标数据．dat');

M=load('浓度含量shuju．dat');

x=D(:,2);

y=D(:,3);

z=D(:,4);

c=D(:,5);

xi=linspace(min(x),max(x),300);

yi=linspace(min(y),max(y),300);

[xi,yi]=meshgrid(xi,yi);

zi=griddata(x,y,z,xi,yi);

ci=griddata(x,y,c,xi,yi);

marker={'d','^','s','o','V'};

color={'b','y','c','r','m'};

mat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn'};

str={'等高线','生活区','工业区','山林区','交通区','绿地区'};

for j=1:8

% 等高线图形

figure

contourf(xi,yi,zi,0:10:500);

% set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)

% clabel(C,h,[0:10:50,50:50:300,300:100:500]);

legend( '生活区','工业区','山林区','交通区','绿地区');

%legend( {'k','r','y','c','b'});

title(['重金属 ',mat{j}, ' 空间二维等高线分布图'])

xlabel('X---轴')

ylabel('Y---轴')

colormap summer

colorbar

grid on

hold on

for i=1:5

loc=c==i;

plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i});

end

% legend(str,'location','best')

for k=1:length(x)

text(x(k)-200,y(k)+200,num2str(M(k,j+1)),'fontsize',8);

end

legend( '功能一区——生活区','功能二区——工业区','功能三区——山林区','功能四区——交通区','功能五区——绿地区');

end

2． 重金属 ‘**’含量程序

clc; clear;

load shuju1．txt

A=shuju1;

x=A(:,1);y=A(:,2);z=A(:,3);

scatter(x,y,5,z)%散点图

figure(1)

[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x),200)',linspace(min(y),max(y),200),'v4');%插值

pcolor(X,Y,Z);

shading interp%伪彩色图

figure(2),contourf(X,Y,Z) %等高线图

legend('污染程度分布');

title('重金属As含量空间分布平面图')

xlabel('x—轴'),ylabel('y—轴');%zlabel('z—轴')

figure(3),surf(X,Y,Z)%三维曲面

legend('污染程度分布');

title('重金属As含量空间分布平面图')

xlabel('x—轴'),ylabel('y—轴');zlabel('z—轴')

3 重金属综合污染指数程序

%Cd综合污染指数

Grandvalue=190;

maxvalue=0;

sum=0;

C=[

352．60

556．20

113．70

240．50

120．30

514．10

246．30

158．00

211．30

236．50

]; % 可以输入任意长度的污染浓度数据

C=C'

n=length(C);

for i=1:n;

if C(i)>maxvalue

maxvalue= C(i);

else

sum=sum+ C(i)/Grandvalue;

end

end

p=sqrt(((sum/n)^2+(maxvalue/Grandvalue)^2)/2)

附录二： 污染状况空间分布图

附录三：污染程度数据分布表

附表1 Cd污染程度数据分布表

附表2 Cr污染程度数据分布表

附表3 Cu污染程度数据分布表

附表4 Hg污染程度数据分布表

附表5 Ni污染程度数据分布表

附表6 Pb污染程度数据分布表

附表7 Zn污染程度数据分布表

附表8 三功能区各重金属污染程度所占百分比

附表9 四功能区各重金属污染程度所占百分比

附表10 五功能区各重金属污染程度所占百分比