2011高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的
表
关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf
述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：
参赛队员(打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人 (打印并签名)：
日期：2013年 8 月26 日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2011高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
评
阅
人
评
分
备
注
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
城市表层土壤重金属污染分析
摘 要
本文主要是对重金属污染浓度的分析，画出浓度的空间分布图，画出浓度的空间分布图，建立相关的模型，得出该地区的重金属污染程度、污染原因以及污染源的所在位置。
问题一要求我们对重金属污染浓度进行具体分析。首先采用Cubic插值法，利用Matlab软件画出五个功能区在城区的位置分布图及8种重金属元素空间分布图，以颜色深浅对重金属污染程度进行直观表现，颜色浅则表示浓度大并将两种分布图进行对比，得到8种重金属的具体分布区域。其次建立双权重因子模糊综合评判模型，考虑浓度和毒性双重因子作为权重判定各区域内不同元素的污染程度，最终得出：工业区普遍受到较重的重金属污染，其次为主干道路区，再次是生活区和公园绿地区，山林区基本未受到重金属的污染。
问题二要求通过数据分析重金属污染的主要原因，我们发现各种重金属不是单一产生的，多种重金属在产生的过程中往往具有伴随性。针对这种现象，我们使用因子分析法来确定各种重金属元素之间的相关性，然后利用它们的相关性定量分析各因子对8种重金属元素产生的贡献值。最终通过已知数据可以得出： Ni、Cr、Cu的污染主要是由主干道路区汽车尾气产生， Hg的污染主要是由工业“三废”产生， As的污染主要是由生活区杀虫剂等药剂产生， Zn的污染主要是由生活区垃圾焚烧和主干道路区汽车轮胎磨损产生。
问题三在分析重金属污染的传播特征时，可以结合大气扩散模型和水体对流扩散模型的偏微分方程建立重金属元素的对流扩散模型。为求具体污染源的位置，在对扩散模型简化的基础上，利用反演扩散模型方程来求污染源。反演方程可以用遗传算法来实现，把问题转化为建立一个非线性最优化问题。以为例，其污染源坐标为
问题四要求我们分析已用模型的优缺点，并改进模型以研究更好的城市地质环境的演变模式。通过搜集不同历史时间不同位置的数据，利用时间序列分段线性拟合各个指标值，较准确地恢复任一历史时间的土壤质量指标，根据重金属元素浓度变化的大致趋势，采用三次指数平滑法对未来各指标变化进行预测。
关键词：双权重因子模糊综合评判模型 因子分析 遗传算法 反演
扩散方程
一、 问题提出
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、交通区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0-10 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、 问题分析
2.1对问题一的分析
题目要求统计出八种重金属元素在城区的空间分布，需要考虑重金属元素在X,Y,Z三维空间的浓度分布特征。由附件给出的数据分析得，当X,Y的坐标确定之后，该点的海拔高度Z相应确定下来。所以在统计重金属的空间分布时，可以先考虑元素浓度和平面位置的关系，再通过海拔分布图分析元素分布与海拔之间的关系。
对这八种重金属元素分布单独考虑，利用附件中各元素在城区的样本浓度，采用Cubic插值法采用Matlab软件对各元素含量分布进行空间分析，可以分别绘制出8种重金属元素的空间等高线分布图。
当考虑不同区域内各元素的污染程度时，现行重金属污染评价方法一般采用污染物浓度超标赋权法。其中，地累积指数评价模型法[1]( Forstner，1989,Forstner et al.1993 )是Muller 提出的一种能够较好地表征土壤中重金属富集程度的定量指标，但其评价法只考虑了浓度而未考虑到重金属元素本身的生物毒性。对于不同重金属，因污染物个体的毒性级别不同，污染物浓度超标赋权法有可能掩盖某些低浓度有机组分的毒性作用。因此本文采取双权重因子模糊综合评判数学模型，将污染物浓度和毒性级别指数加权叠加，并作归一化处理，得到某污染组分的权重公式考虑引入模糊综合评价来进行基于各元素浓度、毒性响应的多因素综合评价问题。结合环境土壤法给出的标准，对各区的综合污染进行进一步精确的分析。
2.2对问题二的分析
对于问题（2），分析重金属污染的主要原因，我们主要要依据问题（1）中得出的各个功能区的不同重金属元素的污染状况，得出某一重金属元素主要是由哪一种污染物所带来的。首先进行每个区域内八种重金属元素的相关性分析，相关性比较显著的重金属元素可以认为来自同一个污染源，由此可以大致确认每个区域污染源的种类。因为每个区域中重金属的种类众多，相互之间又具有比较明显的相关性，所以我们可以采用KMO法和Bartlett法对原始数据进行主成分分析适宜性检验的基础上，进一步考虑使用主成分分析法，由此分析出每个区域营销比较显著的若干种金属，从而确认出重金属元素的主要原因。
2.3对问题三的分析
基于一个基本假设，我们发现重金属污染物的扩散符合对流扩散模型，则通过求解扩散微分方程，可以得到一个重金属污染物浓度的函数表达式，我们发现，对于特定的重金属污染物，如果假设出几个污染源的位置以及相应的排放强度，即可求出整个区域的污染物浓度。根据附件已经给出的某一区域的污染物浓度实测值，我们可以检验假设的污染源位置是否正确，通过不断验证，一定可以找到几个污染源的位置，它们可以使得求出的这种重金属污染物的浓度分布与实际的分布最为接近，那么这几个污染源的位置即就是我们需要的污染源的位置。那么原问题便转化为了一个非线性最优化的问题，在matlab软件中，我们借助遗传算法，可以通过不断迭代，最终找出这样几个位置，满足要求，进而可以得出不同重金属的污染源位置。
2.4对问题四的分析
前三个问题的求解模型可以比较合理地解决问题，但是也在一定程度上存在局限性。由此我们可以进一步收集该城区突然惩罚的信息，建立更优化的模型。
三、 模型的假设
（1）假设所有数据都是真实可靠的；
（2）假设污染源为重金属含量最高的地方；
（3）假设污染源的位置大概固定，不随时间变化。
（4）假设污染元素在水平面的扩散在各个方向的扩散速度相同。
（5）重金属污染物的扩散方式为简单的对流扩散，符合对流扩散方程；
（6）污染物的浓度在海拔高度上大致都分布在0~100以内，我们只考虑
方向上的扩散；
（7）由于扩散到一定程度，扩散的影响因素将会增多，我们假设存在理想最大扩散距离
；
四、 符号说明
表1
元素
As
Cd
Cr
Cu
Hg
Ni
Pb
Zn
标号/j
1
2
3
4
5
6
7
8
表2
污染因子实测值
污染因子对应于各级土壤重金属环境质量状况的标准值
因子毒性指数
因子权重值
的相关系数
重金属浓度在
三个方向的扩散系数
重金属浓度在
三个方向的对流系数
污染源的排放强度
污染源的理想最大扩散距离
函数
与解原微分方程等效的新函数
五、 模型建立与求解
5.1问题一的模型建立与求解
5.1.1重金属空间分布的求解
5.1.1.1模型建立与求解
根据采样点各重金属元素的含量分布，采用Cubic插值法，应用Matlab软件对海拔高度和8种重金属的含量分布进行空间分析，可汇出以下空间分布图。
图1：海拔分布图
图2：As空间分布图
图3：Cd空间分布图
图4：Cr空间分布图
图5：Cu空间分布图
图6：Hg空间分布图
图7：Ni空间分布图
图8：Zn空间分布图
通过比较8种重金属元素分布图与功能分布图，可以得出8种重金属元素的分布情况，如下表3：
表3: 8种重金属元素的分布
元素
主要分布区域
工业区、西南方的公园绿地区、西南方和东南方的交通区
工业区、西南方的交通区、公园绿地区、西南方和西北方的山区
城区西南的工业区、少量分布在城区东北的工业区
城区南部以及中部的工业区
城区西南部的工业区和生活区的交接处
城区西南部的交通区和工业区以及山区密集的地方
城区南部的生活区以及工业区
城区南部的生活区、交通区和工业区
5.1.1.2模型评价
运用Cubic插值法，画图精准，有助于空间浓度分析。
5.1.2不同区域内重金属的污染程度的求解
双权重因子模糊综合评判模型是基于各元素浓度、毒性响应的双因素综合评价问题[2]。将污染物浓度和毒性级别指数加权叠加，并作归一化处理，得到某污染组分的权重公式考虑结合环境土壤法给出的标准，对各区的综合污染进行进一步精确的分析。
5.1.2.1建立模型
双权重因子模糊综合评判决策的方法与步骤如下：
（1）设评判因子
，代表8种金属，评价集
，代表国家标准中5个土壤等级，分别为清洁、尚清洁、轻污染、中污染和重污染。
（2）各因素的隶属度确定采用呈梯形分布的隶属函数。评价单项污染因子
为隶属度函数的拐点，用分段直线函数模拟，。其中
为污染因子对应于各级土壤重金属环境质量状况的标准值，
表示污染因子实测值。
当
时，
当
时，
公式（1）
当
时，
（3）得到由评价
（各金属）在评价集
上的隶属度所组成的8*5的关系模糊矩阵
。
（4）
为评价集
对因子的权重分配构成的权重向量，
，
为第
种污染因子所对应的权重。确定权重时，以
，
表示因子权重值，
为因子毒性指数，将污染物浓度和毒性指数加权叠加，作归一化处理，得权重公式：
公式（2）
公式（3）
(5)据矩阵运算法则，
确定了一个模糊映射，把评价因子集
上的权重向量
映射到评价集
，得模糊子集
，即评价等级
，据最大隶属度原则确定评价等级。
5.1.2.2模型求解
（1） 先对数据进行预处理。将公式（1）中的
用每个区8种重金属元素含的均值替代，对五个区8种重金属元素分别进行处理，结果如下表4：
表4：各区各个重金属含量的均值
区域
生活区
工业区
山区
主干道路区
公园绿地区
As (μg/g)
6.270455
7.251389
4.044091
5.708043478
6.263714286
Cd (ng/g)
0.289961
0.393111
0.15232
0.360014
0.280543
Cr (μg/g)
69.01841
53.40917
38.9597
58.05391304
43.636
Cu (μg/g)
49.40318
127.5358
17.31727
62.21492754
30.19171429
Hg (ng/g)
0.093041
0.642355
0.040956
0.446823
0.114992
Ni (μg/g)
18.34227
19.81167
15.45379
17.61710145
15.28971429
Pb (μg/g)
69.10636
93.04083
36.55591
63.5342029
60.70857143
Zn (μg/g)
237.0086
277.9275
73.29424
242.8548551
154.2422857
（2）搜集资料，确定等级标准。利用国家《土壤环境质量标准》(GB 15618—1995)确定的区域土壤评价，选取甘肃省土壤重金属元素背景值和临界含量作为
的标准参照；考虑重金属生物体的毒性作用时选取Hakanson的重金属毒性响应系数所确定的毒性指数作为
的标准参照，最终确立土壤重金属污染程度分级标准及生物毒性指数，如下表5：
表5：土壤重金属污染程度分级标准及生物毒性指数(mg·kg ) [3]
元素/等级
1清洁
2尚清洁
3轻污染
4中污染
5重污染
毒性指数
As (μg/g)
82.56
95.43
124.26
164.28
201.38
3
Cd (ng/g)
0.1204
0. 2523
0.6
1.4
2
2
Cr (μg/g)
74.88
99.54
150
350
500
5
Cu (μg/g)
28.37
40.63
120
280
400
4
Hg (ng/g)
0.092
0.2592
0.45
1.05
1.5
1
Ni (μg/g)
77.51
82.56
94.21
102.56
154.25
3
Pb (μg/g)
23.35
36.09
150
350
500
4
Zn (μg/g)
83.68
116.75
240
560
800
6
（3）根据公式（1）运用Matlab软件计算各元素隶属值，得到各区域的关系模糊矩阵
,以生活区为例，其他区域详细见附件。
=
（4）根据公式（2）、（3）运用Matlab软件，将污染物浓度和毒性级别指数加权叠加并作归一化处理，得到某污染组分的权重如下表6：
表:6：该城市各类区各参评因子的权重值
区域
生活区
工业区
山区
主干道路区
公园绿地区
As (μg/g)
0.0129
0.0145
0.0107
0.0112
0.0125
Cd (ng/g)
0.3694
0.2485
0.1993
0.1286
0.4113
Cr (μg/g)
0.0902
0.0922
0.0823
0.141
0.0652
Cu (μg/g)
0.053
0.0592
0.0659
0.0741
0.0603
Hg (ng/g)
0.3717
0.4801
0.5438
0.5371
0.3701
Ni (μg/g)
0.0543
0.0684
0.0435
0.0459
0.0436
Pb (μg/g)
0.0378
0.0415
0.0452
0.0392
0.0362
Zn (μg/g)
0.0777
0.0785
0.0635
0.0801
0.0568
（5）模糊矩阵复合运算，继续用Matlab软件处理得到各评价样品对评价等级的隶属度,又根据最大隶属度原则，确定各样品的的污染程度，此即为土壤环境质量分级，最终结果该城区内不同区域8种主要重金属元素的污染程度分析结果如下表：
表7：模糊综合评价结果
区域
隶属度
污染等级
生活区
0.526786
0.407573
0.132641
0
0
1（清洁）
工业区
0.1751
0.168614
0.573177
0.166009
0
3（污染）
山区
0.96077
0.093246
0.000185
0
0
1（清洁）
交通区
0.1981
0.181409
0.676976
0.000715
0
3（污染）
公园绿地区
0.491847
0.505642
0.058511
0
0
2（轻污染）
（5）综合分析：根据表7可以看出，生活区和山区污染度最低，属于标准里的清洁范畴；公园绿地区轻微污染，而工业区和交通区污染最严重，其中又以交通区的污染程度为最；而纵观表6各区的隶属度值的大小，我们可以得出以下结论：
对生活区的分析：
污染排序
,其中
元素尚清洁
对工业区的分析：
污染排序
,其中
元素尚清洁
对山区的分析：
山区没有重金属污染，尚清洁
对交通区的分析：
污染排序
,其中
元素尚清洁
对公园绿地区的分析：
污染排序
,其中
元素尚清洁。
5.1.2.3模型评价
优点：双权重因子模糊综合评判模型考虑浓度和毒性双重因子作为权重，查找资料引入毒性这一影响污染程度的指标，使不同评价对象间可比性增强，采用评价分值，即利用评价向量的分量形成权重，通过确定各等级对应的分值，对各等级的分值进行加权平均，得到评价分值。从评价分值可以直观的看出各采样点土壤重金属环境质量的优劣。
缺点：双权重因子模糊综合评判模型计算冗长，耗时较长，所以对于采样样本数目较多的情况并不可行，并且没有考虑人为因素的影响。
适用：建议对于少量样本测量时进行使用。
5.2问题二的模型建立与求解
城市土壤重金属来源于成土母质和人类活动，同一来源的重金属之间存在着
相关性，根据相关性可以判断土壤重金属污染来源是否相同。如果重金属之间存
在显著的正相关，则其来源可能相同，否则来源可能不止一个。
5.2.1建立模型
通过查阅资料[4][5],我们采用多元统计数学方法之一的因子分析法，它根据多个实测变量（重金属元素）之间的相关性，运用数学变换，将多个变量转换为少数几个线性不相关的综合指标（不相关的几个主要原因），从而简化数据处理，其目的在于对大量观测数据，用较少的有代表性的因子（目标）来说明众多变量所提取的主要信息，提示出多个变量的因果关系。具体方法如下：
因子分析从变量的相关矩阵出发将一个m维随机向量X（数据总库）分解成低于m个且有代表性的公因子和一个特殊的m维向量，是公因子数取得最佳的个数，从而使对m维随机向量的研究转换成较小个数的公因子研究。
设一共有n（n=319）个样本，n个指标构成样本空间X：
；
。
之后我们对数据做如下处理步骤：(1)原始数据的标准化，标准化公式为：
其中
为第i个样本第j个指标值，而
和
分别为j指标的均值和标准差。标准化目的在于消除不同变量量纲的影响，而且标准化转化不会转变变量的相关系数。
（2）KMO检验，判断其是否符合相关性分析条件。KMO值定义为
其判断标准是：
表8：KMO检验法判断标准
0.9
0.7（符合因子分析法的条件），再给出8种重金属元素原始含量数据的相关系数矩阵。各元素相关系数矩阵如表5所示，可见As和Ni、Pb的相关性分别为0.317、0.290（一般），而As和Cd、Zn的相关性更低，As和Hg几乎无相关性，同样可以看出其他两两之间的相关性，Cd和Pb的相关性较好达到0.660，Cr和Ni、Cu和Pb、Cu和Cr相关性分别达到0.716、0.520和0.532，而Hg只与Cu的相关性为0.417（一般），Hg与其它元素的相关性很弱，剩下的元素之间的相关性基本上都一般。从成因来分析，相关性较好的元素可能在成因和来源上有一定的关联。
表9：8种重金属元素相关性系数矩阵
相关系数
As
Cd
Cr
Cu
Hg
Ni
Pb
Zn
As
1.000
0.255
0.189
0.160
0.064
0.317
0.290
0.247
Cd
0.255
1.000
0.352
0.397
0.265
0.329
0.660
0.431
Cr
0.189
0.352
1.000
0.532
0.103
0.716
0.383
0.424
Cu
0.160
0.397
0.532
1.000
0.417
0.495
0.520
0.387
Hg
0.064
0.265
0.103
0.417
1.000
0.103
0.298
0.196
Ni
0.317
0.329
0.716
0.495
0.103
1.000
0.307
0.436
Pb
0.290
0.660
0.383
0.520
0.298
0.307
1.000
0.494
Zn
0.247
0.431
0.424
0.387
0.196
0.436
0.494
1.000
然后，利用相关系数矩阵在SPSS13.0统计软件中计算出相应的因子的特征值和累值贡献率，得到表10，如下：
表10：旋转前后因子的特征值和累值贡献率
因子
旋转前
旋转后
总的特征值
占总百分量的百分率
累计贡献率
总的特征值
占总百分量的百分率
累计贡献率
1
3.560
44.500
44.500
2.059
25.733
25.733
2
1.150
14.377
58.877
1.767
22.087
47.820
3
0.960
12.063
70.941
1.217
15.215
63.035
4
0.768
9.596
80.537
1.026
12.824
75.858
5
0.578
7.220
87.756
0.952
11.898
87.756
在积累方差为87.756%（>80%）符合因子分析法的前提下，分析得到5个主因子，这5个因子共占了源资源的87.756%，这5个因子也就是我们要分析得到的重金属污染的五个主要原因。从表中可以看出主因子1和2占的贡献率达到25.733%和22.087%，显然是重金属污染的最主要的原因，因子3、4、5（贡献率在11至16之间）为重要原因。
最后，因子分析的主要目的是将具有相近的因子载荷的各个变量置于一个公因子之下，正交方差最多旋转使每一个主因子只与最少个数的变量有关系，而使足够多的因子负荷均很小，以便对因子的意义作出合理的解释，输出见表11。
表11：各元素与各因子的关联度
元素
因子1
因子2
因子3
因子4
因子5
As
0.131
0.154
0.023
0.970
0.076
Cd
0.170
0.877
0.082
0.102
0.132
Cr
0.882
0.209
0.002
0.014
0.146
Cu
0.614
0.362
0.505
-0.022
0.029
Hg
0.013
0.134
0.953
0..030
0.083
Ni
0.864
0.089
0.019
0.222
0.196
Pb
0.195
0.832
0.191
0.121
0.208
Zn
0.271
0.285
0.101
0.091
0.904
从表7数据可得，元素Ni、Cr和Cu的对因子1的关联度分别达到0.882、0.864和0.614，是此因子产生的主要污染物，从相关性表里也可以看出Ni和Cr的相关性好（相关性系数达到0.716），查找资料，可发现汽车燃烧的石化燃料所产生的尾气中含有Ni、Cr和Cu，那么因子1极有可能是汽车燃烧的石化燃料所产生的尾气。
元素Cd和Pb对于因子2的关联度分别达到0.877和0.832，是此因子产生的主要污染物，从相关性表里也可以看出Cd和Pb的相关性好（相关性系数达到0.660），查找资料[6],可发现工业中电镀、采矿、冶炼、燃料、电池和化学工业中排放的废气、废水和废渣以及汽车排放的尾气中都含有比较多的Cd和Pb，由于因子1是汽车燃烧的石化燃料所产生的尾气，所以可以判断因子2主要是工业中电镀、采矿、冶炼、燃料、电池和化学工业中排放的废气、废水和废渣。
元素Hg对于因子3关联度达到0.953，是此因子产生的主要污染物，查找资料[6],可发现工业中氯碱、塑料、电子等工业排放的汞废水中含有大量的Hg，所以可以判断因子3为工业中氯碱、塑料、电子等工业排放的汞废水。
元素As对于因子4的关联度达到0.970，查找资料[6],可发现生活中使用的杀虫剂、杀菌剂、灭鼠剂和除草剂等用品产生的废气烟尘中含有As，而且其它元素对于此因子的关联度偏低，也正好验证了因子4很有可能是生活中使用的杀虫剂、杀菌剂、灭鼠剂和除草剂等用品产生的废气烟尘。
元素Zn的对于因子5的关联度达到0.904，其它的都比较低，查找资料[6],可发现汽车轮胎磨损产物以及生活中燃烧煤制品和垃圾焚烧产生的废气、烟尘，所以因子5为汽车轮胎磨损产物以及生活中燃烧煤制品和垃圾焚烧产生的废气、烟尘。
综上所述，通过分析数据可以得出此城区重金属污染的主要原因有以下五点：
a) 汽车燃烧的石化燃料所产生的尾气，主要引起重金属Ni、Cr、Cu的污染；
b) 工业中电镀、采矿、冶炼、燃料、电池和化学工业中排放的废气、废水和废渣，主要引起重金属Cd、Pb的污染；
c) 工业中氯碱、塑料、电子等工业排放的汞废水，主要引起重金属Hg的污染；
d) 生活中使用的杀虫剂、杀菌剂、灭鼠剂和除草剂等用品产生的废气烟尘,主要引起重金属As的污染；
e) 汽车轮胎磨损产物以及生活中燃烧煤制品和垃圾焚烧产生的废气、烟尘，主要引起重金属Zn的污染。
5.3问题三的模型建立与求解
5.3.1.污染物传播特征的前提假设
重金属元素传播途径主要有水流传播、土壤渗透、大气扩散等。在城区中我们主要考虑污染以大气扩散。为了研究方便，作如下假设（重新申明）：
（1）重金属污染物的扩散方式为简单的对流扩散，符合对流扩散方程；
（2）污染物的浓度在海拔高度上大致都分布在0~100以内，我们只考虑
方向上的扩散；
（3）由于扩散到一定程度，扩散的影响因素将会增多，我们假设存在理想最大扩散距离
；
5.3.2.模型建立
我们根据对流扩散方程建立模型[7]：
其中
为重金属浓度，
在
三个方向的扩散系数,
为对流系数。
模型一：由于附件中给的数据都是离散化的点样本，我们将空间的平面看成有限个点的网格组成，于是我们将原来的方程作离散化处理，用差商代替原方程的微商，通过求解离散点上的解来近似原方程的解，转化如下：
根据对流扩散假设，污染源中心处重金属浓度会随着时间慢慢减小，而且其浓度随时间的变化率与周围相比比较大，那么原问题就转化为求等步长时间内浓度
向
方向的变化率问题。
5.3.3.模型求解
上式c为重金属浓度，Ex，Ey，Ez为三个方向上的扩散系数，ux，uy，uz为对流系数，Mi为污染源的排放强度。R=f（x，y，z），L为污染源的理想最大扩散距离。
简化上式可得
求解上式微分方程，对解析解进一步变换，结果如下：
则上式求c的问题可以转化为求v函数的问题
对结果再一步变化可得：
由边界值可解出：
其中：
最后求得v函数的解析解为：
我们已经有了各个样本点的浓度数据
，由上式可知，
可由
，
……
，我们先研究一种特定的金属污染物，则其
为定值，只要我们假设出它的污染源的位置，即
，那么我们就可以得到所有样本点的浓度数据，如果这些样本点的浓度数据和实际的符合度很好，那么就说明我们的假设是正确的，那么原问题就转化成了下面的非线性最优化问题：
则满足上式的
即就是我们需要的污染源的位置，再转化为原始三维坐标即可通过matlab的工具箱，我们采用遗传算法可以大大提升算法效率和准确度，得到不同金属污染的污染源位置坐标如下表：
表12：各重金属污染的污染源位置坐标表
坐标
区域
编号
浓度
7056
8348
37
1
42
7.41
12644
14943
43
2
275
5.41
12734
10344
32
3
251
7.34
23664
9790
46
3
92
114.3
20261
7586
29
3
209
188
14405
18032
152
3
280
398.4
26416
6508
14
1
107
36.73
18413
11721
88
3
80
59.61
20177
17642
276
3
288
23.15
5375
8643
15
1
25
11.82
19569
7348
70
3
210
9.95
22624
4818
27
3
219
28.97
9328
4311
24
1
36
170
18738
10921
53
3
82
58.93
14262
15129
66
3
282
19
9212
11305
5
2
231
27.91
15007
5535
70
3
138
31.5
25981
18051
173
3
301
9.7
13855
3345
79
1
183
54.98
17008
4775
82
3
194
38.35
20983
15862
93
3
287
27.54
12625
16259
66
1
276
149.52
16947
7487
41
2
197
104.9
22624
4818
27
3
219
105.31
5.4问题四的解答
5.4.1优缺点分析
本文比较有特点的几个模型中，优点明显的有：
（1） 双权重因子模糊综合评判模型考虑了环境质量界限模糊性的特性，打破用一个确定性指标来评价环境质量的缺陷，具有一定的合理性和优越性，并且结合了不同元素的毒性响应系数，对污染程度的判断更合理。
（2） 反演扩散模型比较合理地模拟重金属在该城区内的扩散。通过对模型的简化，比较简单的将模型转换为二维偏微分方程。由于遗传算法有很好的的全局优化性能，减少了陷入局部极值的风险，特别是对非线性优化问题有很强的适应性。传统的反演方法在源强和浓度的关系上，只能给出一种虚拟的联系，而遗传算法可以给出较为直接的关系，从而为污染预报和污染的控制提供了基础。用遗传算法来反演，方法简单且通用性强。
而缺点则有：
由于不考虑模型里的对流项以及海拔上面的影响，导致其自身的确定误差较大，其次遗产算法本身就是一个枚举算法，算法在运行中具有随机性，导致其求出的污染点的分布不是完全符合，但具有稳定的收敛性，是可以接受的。
5.4.2城市地质环境的演变模式模型
5.4.2.1 所需要搜集的数据
查阅相关的数据[8]，得到了有关地层的几个重要性指标：土壤重金属元素含量，pH值，有机质及养分含量，不同种类无机矿物的含量，N，P等营养元素的含量，农药、有机污染物残留含量等多方面影响土壤地质环境的因素。要能研究地质演化的信息，必须搜集不同历史时间的数据。对于数据的要求为，将所要研究的区域数据的取样点越密集，研究结果越准确。
5.2.2.2 时间序列分析
根据已有研究[8]，在城市发展对土壤质量影响的分析、规律统计和建模方面，主要有三种模型：基于基本统计理论的回归分析和自回归趋势分析，灰色关联分析和灰色建模、基于GIS的图形空间分析模型。我们可以利用时间序列分段线性拟合各个指标值，较准确地恢复任一历史时间的土壤质量指标，根据重金属元素浓度变化的大致趋势，采用三次指数平滑法对未来各指标变化进行预测。
六、 参考文献
[1].张江华等，《内梅罗指数和地质累积指数在土壤重金属评价中的差异探讨》，《环保与分析》，第31卷第8期：P43-44，2010
[2]. 张志斌，王轲，蒲瑞丰，《基于双权重因子模糊综合模型的金昌市土壤重金属污染评价》，《西北师范大学学报》，2009
[3]. 国家《土壤环境质量标准》(GB 15618—1995)
[4].王雄军、赖健清、鲁艳红、李德胜、周继华、王建武，基于因子分析法研究太原市土壤重金属污染的主要来源，生态环境，第17卷第2期：671-676,2008年；
[5].史贵涛、陈振楼、许世远、王利、张菊、李海雯、李丽娜，上海城市公园土壤及灰尘中重金属污染特征，环境科学，第28卷第2期：238-242，2007年；
[6].百度百科“砷” #14，
百度百科“镉污染”#3，
百度百科“铬污染”#2，
百度百科“铜污染” ，
百度百科“汞污染” #2，
百度百科“镍污染” ，
百度百科“铅污染” ，
百度百科“锌污染” ，
2011/9/10；
[7].闵涛等，《二维稳态对流---扩散方程反演的迭代算法》，《水动力学研究与进展》2007.11
[8].吴新民等，《长江三角洲土壤质量演变趋势分析》，《地理与地理信息科学》，第51 卷第2 期：P88-89，2008
七、 附件
7.1作图程序
clc
clear
close all
x=xlsread('2011a.xls','B4:B322');
y=xlsread('2011a.xls','C4:C322');
z=xlsread('2011a.xls','D4:D322');
c=xlsread('2011a.xls','E4:E322');
M=xlsread('2011a.xls',2,'B4:I322');
xi=linspace(min(x),max(x),150);
yi=linspace(min(y),max(y),150);
[xi,yi]=meshgrid(xi,yi);
for i=1:8
Mi(:,:,i)=griddata(x,y,M(:,i),xi,yi);
end
ci=griddata(x,y,c,xi,yi);
zi=griddata(x,y,z,xi,yi,'cubic');
marker={'*','o','s','^','p'};
color={'k','r','y','c','b'};
mat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn'};
str={'等高线','生活区','工业区','山林区','交通区','绿地区'};
for j=1:8
%等高线图形
a=M(:,j);
figure
contourf(xi,yi,Mi(:,:,j));
% set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
% clabel(C,h,[0:10:50,50:50:300,300:100:500]);
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')
end
%三维体图
figure
surf(xi,yi,zi);
colormap hsv
title('三维图立体（颜色条表示分类）')
xlabel('X')
ylabel('Y')
colorbar
hidden off
hold on
for i=1:5
loc=c==i;
plot3(x(loc),y(loc),z(loc),marker{i},'markerfacecolor',color{i});
end
str{1}='三维图';
legend(str,'location','best')
7.2模糊关系矩阵程序
clc
clear
x=xlsread('x.xlsx');
s=xlsread('s.xlsx');
for dtt=1:5
for i=1:8
for stn=1:5
if stn==1
if x(i,dtt)<=s(i,stn)
A(i,stn,dtt)=1;
elseif x(i,dtt)>=s(i,stn)&&x(i,dtt)<=s(i,stn+1)
A(i,stn,dtt)=(s(i,2)-x(i,dtt))/(s(i,2)-s(i,1));
elseif x(i,dtt)>=s(i,stn+1)
A(i,stn,dtt)=0;
end
elseif stn==5
if x(i,dtt)<=s(i,stn-1)
A(i,stn,dtt)=0;
elseif x(i,dtt)<=s(i,stn)&&x(i,dtt)>=s(i,stn-1)
A(i,stn,dtt)=(x(i,dtt)-s(i,stn))/(s(i,stn)-s(i,stn-1));
elseif x(i,dtt)>=s(i,stn)
A(i,stn,dtt)=1;
end
elseif stn==2
if x(i,dtt)>=s(i,stn+1)||x(i,dtt)<=s(i,stn-1)
A(i,stn,dtt)=0;
elseif x(i,dtt)>s(i,stn-1)&&x(i,dtt)=s(i,stn)&&x(i,dtt)<=s(i,stn+1)
A(i,stn,dtt)=(s(i,stn+1)-x(i,dtt))/(s(i,stn+1)-s(i,stn));
elseif x(i,dtt)==s(i,stn)
A(i,stn,dtt)=1;
end
elseif stn==3
if x(i,dtt)>=s(i,stn+1)||x(i,dtt)<=s(i,stn-1)
A(i,stn,dtt)=0;
elseif x(i,dtt)>s(i,stn-1)&&x(i,dtt)=s(i,stn)&&x(i,dtt)<=s(i,stn+1)
A(i,stn,dtt)=(s(i,stn+1)-x(i,dtt))/(s(i,stn+1)-s(i,stn));
elseif x(i,dtt)==s(i,stn)
A(i,stn,dtt)=1;
end
elseif stn==4
if x(i,dtt)>=s(i,stn+1)||x(i,dtt)<=s(i,stn-1)
A(i,stn,dtt)=0;
elseif x(i,dtt)>s(i,stn-1)&&x(i,dtt)=s(i,stn)&&x(i,dtt)<=s(i,stn+1)
A(i,stn,dtt)=(s(i,stn+1)-x(i,dtt))/(s(i,stn+1)-s(i,stn));
elseif x(i,dtt)==s(i,stn)
A(i,stn,dtt)=1;
end
end
end
end
end
7.3遗传算法主程序
clear all;
%遗传算法主函数
um=100;
Iter=10;%迭代次数
pinter=0.9;%交叉概率
pmut=0.5;%变异概率
%产生初始种群
siz= 161;
for i=1:um
Temp=rand;
Te=3;
N(i,:)=(randperm(siz)<=Te);
end
igBest=1;
NBest=N(1,:);
F=zeros(1,um);
data=xlsread('a.xlsx');
V=data(1:siz,9);
X=data(1:siz,11);
[V,index]=sort(V,'descend');
X=X(index);
%计算适应度
for k=1:um
F(k)=Vfunction(N(k,:),V,X);
end
it=1;
while itrand
ind1=find(s);
ind2=find(s~=1);
if length(ind1)~=0& length(ind2)~=0
iRa1=round(rand*(length(ind1)-1))+1;
iRa2=round(rand*(length(ind2)-1))+1;
s(ind1(iRa1))=1-s(ind1(iRa1));
s(ind2(iRa2))=1-s(ind2(iRa2));
end
else
s=s;
end
%交叉
function[s1,s2]=intercross(s1,s2,pinter)
iTemp=0;
if pinter>rand
ind1=find(s1);
ind2=find(s2);
inter=intersect(ind1,ind2);
ind1=setdiff(ind1,inter);
ind2=setdiff(ind2,inter);
if length(ind1)~=0&length(ind2)~=0
iRa=round(rand*(length(ind1)-1)+1);
iRa2=round(rand*(length(ind2)-1)+1);
s1(ind1(iRa))=1-s1(ind1(iRa));
s1(ind2(iRa2))=1-s1(ind2(iRa2));
s2(ind2(iRa2))=1-s1(ind2(iRa2));
s2(ind1(iRa))=1-s2(ind1(iRa));
end
else
s1=s1;
s2=s2;
end
7.4.问题一中各区域关系模糊矩阵
（以表格形式给出），如下表13-16：
表13：工业区关系模糊矩阵
1
0
0
0
0
0
0.595021
0.404979
0
0
1
0
0
0
0
0
0
0.952901
0.047099
0
0
0
0.679408
0.320592
0
1
0
0
0
0
0
0.500037
0.499963
0
0
0
0
0.881477
0.118523
0
表14：山区关系模糊矩阵
1
0
0
0
0
0.758001
0.241999
0
0
0
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
0
0.99591
0.00409
0
0
表15：主干道路区关系模糊矩阵
1
0
0
0
0
0
0.690209
0.309791
0
0
1
0
0
0
0
0
0.728047
0.271953
0
0
0
0.016653
0.983347
0
0
1
0
0
0
0
0
0.759071
0.240929
0
0
0
0
0.991079
0.008921
0
表16：公园绿地区关系模糊矩阵
1
0
0
0
0
0
0.918772
0.081228
0
0
1
0
0
0
0
0.85141
0.14859
0
0
0
0.86249
0.13751
0
0
0
1
0
0
0
0
0
0.783877
0.216123
0
0
0
0.695803
0.304197
0
0
浙公网安备 33021202002483