城市表层土壤重金属污染分析
摘要
本文以城市表层土壤重金属污染为研究对象，分析了8种主要重金属元素在该城区的空间分布、污染程度和原因，建立模型确定污染源的位置，并对模型改进使之能更好的研究城市地质环境的演变模式。
对于问题一，首先对数据预处理剔除异常点并进行无量纲化处理，绘出8种重金属的海拔图和等高线图，得出空间分布为：重金属污染集中在海拔较低处，污染程度随海拔的升高而降低，铬、铜、汞、铅、锌的分布集中在一起。镉分布较广，各功能区分布不一等结论。然后建立尼梅罗模型，由第m 区i 组j 种重金属实际值c ij ，j 种重金属背景值的平均值s j ，得到单因子污染指数F ij =c ij /s j ，求得尼梅罗污染指数为：
P mj =
根据此式得到不同区域各个重金属的尼梅罗污染指数，对每个区域的8种重金属的尼梅罗污染指数求和，得到结果如下： 生活区1.2143，工业区1.4885，山区0.4007，主干道路区0.6715，公园绿地区0.4370。 得出污染的程度为：工业区>生活区>主干道路区>公园绿地区>山区
对于问题二，建立主成分分析模型，对5个区域分别进行主成分分析，运用Spss 得到相关系数矩阵，可知8个变量适合做主成分分析，然后对8个原指标进行解释总方差计算，提取四个主成分，计算成分矩阵，得出结论：
生活区，Zn , Ni , Hg 为主要污染元素。在工业区，Cu , As , Hg , Pb 为主要污染元素。在山区，Zn , As 为主要污染元素。在主干道路区，Zn , As , Hg , Cd 为主要污染元素，在公园绿地区，Zn , Cu , Hg 为主要污染元素。
对于问题三，采用菲克定律，得到重金属速率随时间关系：
C -C i C -C i C -C i dm
=-D ⋅4π⋅r 1r 00≈-D ⋅4π⋅r 020=-D ⋅4π⋅r 020
dt r 0-r i r i r i
将8种重金属的测量值带入，最后得到了如下可能是污染源的坐标：As （17087,11933），
Cr （4592,4603）Cu （3299,6018）（18134,10046）Cd （7405,10981），（20951,10046），
Hg （2383,2692），（13694,2357）Ni （4592,4603）Pb （1647,2728）Zn （2383,3692），（14173,11941）
对于问题四，结合实际情况，分析城市地质环境的演变模式。本文将能影响城市的因素先分为两大类：自然因素和人为因素。再将其各细分为三个次级因素。通过网络调查可以得出各因素在不同时期的实际参数，再通过专家评估和自己的判断，得出每种因素对地质环境演变的影响程度，并对其赋权。最后建立模型得到其演变模式。
关键词 重金属污染 尼梅罗指数 主成分分析 菲克定律
一、问题重述
1.1 背景分析
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。重金属的来源非常广泛，传统上可以分为工业来源和农业来源。随着我国城市化进程的加快，一些有往于城市所特有的污染也随之产生。另外近几年来突发性环境事故频发，重金属污染事件也层出不穷，环境事故也成为重金属污染的重要来源之一。城市工业“三废”排放，金属采矿和冶炼，家庭燃煤，生活垃圾，汽车尾气排放都增加了城市土壤重金属的负荷。重金属污染环境的主要有汞、铅、铬、锌镉、铜等。其中汞的毒性最大，铬、铅、锌等也有相当大毒性。此外还有砷，砷虽不属于金属. 但它的毒性与重金属相似，因此归于重金属一类阐述，称为类金属。目前对我国土壤污染比较普遍的重金属有汞、铬、砷。
对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
1.2 需要解决的问题
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、模型假设
结合本题的实际，为了确保模型求解的准确性和合理性，我们排除了一些位置因素的干扰，提出以下几点假设：
1、假设重金属的物理化学性质稳定。 2、题目所给的数据准确无误，客观真实。
3、假设在采样前该城区没有发生重金属泄露事故。 4、城市地表土壤中的重金属污染物含量成正态分布。
三、符号说明
为了便于问题的求解，我们给出以下符号说明：（其他未说明的符号在文中第一次
四、问题分析
4.1 问题的重要性
土壤重金属污染在一定时期内不表现出对环境的危害性，当其含量超过土壤承受力或限度，或土壤环境变化时，重金属有可能突然活化，引起严重的生态危害。重金属对土壤中的农作物有很强的毒害作用，除此之外，土壤中的重金属还会经由雨水淋滤及地表径流作用转移进入地表水系统，对饮用水安全构成威胁。对城市土壤地质环境异常的查证，成为人们关注的问题。 4.2 问题的思路分析
4.2.1 问题一的分析
对于问题一，本文首先作出了319个采样点的三维地势图，随后为了准确分析第一、二、三、四、五区的重金属污染程度，剔除了20个异常点，随后建立了尼梅罗指数法的模型，得出了5个区环境污染程度的级别。
4.2.2 问题二的分析
问题二要求通过数据分析，说明重金属污染的主要原因。本文再次建立了主成份分析法的数学模型，对5个区内的8种重金属对环境影响的孰强孰弱进行了研究，得到了5个区对环境影响最大的重金属元素。再结合相关文献资料，得出5个区内重金属污染的原因。
4.2.3 问题三的分析
问题三中需要分析重金属污染的传播特征，而在实际生活中，金属的应用主要是合金和离子溶液。往往是通过化学反应来实现其应用，故在问题三中，本文将重金属污染的传播理想化为化学扩散，运用菲克定律，分析不同浓度的地区之间重金属的扩散，追本溯源，从而分析出污染源。
4.2.4 问题四的分析
问题四需要先对问题三中的模型进行评析，得出其结构简单，便于理解，却又过于理想化的特点。然后结合实际情况，分析城市地质环境的演变模式。本文将能影响城市
的因素先分为两大类：自然因素和人为因素。再将其各细分为三个次级因素。通过网络调查可以得出各因素在不同时期的实际参数，再通过专家评估和自己的判断，得出每种因素对地质环境演变的影响程度，并对其赋权。最后建立模型得到其演变模式。
五、模型的建立与求解
经过以上的分析和准备，我们将逐步建立以下数学模型，进一步阐述模型的实际建立过程。
5.1 问题一的求解
5.1.1 8种重金属异常点的处理
对数据预处理，根据8种重金属的浓度随采样点的分布，画出数据图。由于砷、铬、铜、镍、铅、锌这6种重金属的浓度单位为ug /g ，镉、、汞的单位为ng /g ，考虑到不同单位的值若反映在同一图内会有较大误差，故分别作出如下两图：
图1 砷等6金属浓度随采样点的分布 图2 镉、汞浓度随采样点的分布 从两张图可以看出，砷、铬、铜、镍、铅、锌的浓度不服从正态分布，且有多个异常点，为了问题一计算污染程度的准确性，本文将浓度大于500ug /g 的点剔除掉；同理，铬、汞的分布也不服从正态分布并存在多个异常点，为了问题一计算污染程度的准确性，将浓度大于4000ng /g 的采样点剔除。 5.1.2 元素浓度的无量纲化处理
在利用Spss 统计软件数据进行聚类分析的时候，因为单位不统一需要进行无量纲化处理，我们采用均值化方法，即每一个变量除以该变量的平均值，即
x x i '=i
标准化以后各变量的平均值都为1，标准差为原始变量的变异系数。该方法在消除量纲和数量级影响的同时，保留了各变量取值差异程度上的信息，差异程度越大的变量对综合分析的影响也越大。
5.1.3 8种重金属在该城区的空间分布 对数据预处理，进行异常值的过滤，用Matlab 绘出对8种重金属元素的取点的坐标和与其相对应的海拔值的3D 图，以便定性分析重金属分布与海拔高度的关系。
图3 重金属分布图
从图中可知，重金属污染集中在海拔较低处，污染程度随海拔的升高而降低。
根据8种重金属元素浓度值与坐标关系，画出各个重金属元素的二维等高线图，以As 元素为例，其余7种重金属见附录。
y
1
由各个重金属的分布可以得出：
1、铬、铜、汞、铅、锌的分布集中在一起。 2、镉分布较广，各功能区分布不一。 3、工业区呈现混合污染的态势。
5.1.4 不同区域重金属的污染程度
土壤重金属评价方法包括生态危害指数法、地积累指数法、单因子尼梅罗指数法等。由于运用尼梅罗指数的方法，数学过程简洁，加权过程中避免了权系数中主观因素的影响，故采用此方法来评价土壤重金属的污染程度。
0.5
1.5
2
2.5
x 10
4
图4 As 的空间分布
5.1.5 单因子污染指数
首先用给出的数据，计算出单因子污染指数，确定单个金属对各城区土壤的污染，其公式如下:
F ij =c ij /s j
其中定义c ij 第m 区i 组j 种重金属实际值，s j 重金属背景值的平均值，F ij 为第m 区土壤j 种重金属单项污染指数。
5.1.6 尼梅罗综合污染指数
继而求出尼梅罗综合污染指数如下：
P mj =
其中P mj 为在m 区j 种重金属污染物的尼梅罗指数, F max =max(F ij ) 为监测点土壤重
金属最大单项污染指数, F =
为监测点所有重金属单项污染指数平均值。 n
求得不同区域重金属的尼梅罗污染指数结果如下：
∑F
i =1
n
ij
通过对同一个区域的8种重金属元素尼梅罗污染指数求和∑P mj
j =1
(m =1,2,3,4,5) ，
得到结果如下：
未受污染。
故得到结果如下：
污染的程度为：工业区>生活区>主干道路区>公园绿地区>山区, 生活区、工业区的污染较严重，主干道路区已受到污染，山区、公园绿地区未收到污染。
5.2 问题二主成分分析模型的建立与求解
5.2.1 模型的准备
主成分分析是把原来多个变量化为少数几个综合指标的一种统计分析方法，从数学角度来看，这是一种降维处理技术。假定有n 个样本，每个样本共有p 个变量，构成一个n ⋅p 阶的数据矩阵。
⎡x 11x 12 x 1p ⎤⎢x ⎥x x 21222p ⎥ X =⎢⎢ ⎥⎢⎥x x x ⎢np ⎥⎣n 1n 2⎦
记x 1, x 2, x p 为原变量指标，z 1, z 2, z m m ≤p 为新变量指标，得到新指标为：
⎧z 1=l 11x 1+l 12x 2+ +l 1p x p ⎪z =l x +l x + +l x ⎪22112222p p
⎨
⎪
⎪⎩z m =l m 1x 1+l m 2x 2+ +l mp x p
从以上的分析可以看出，主成分分析的实质就是确定原来变量x j （j =1,2, , p ）在主）成分z （上的荷载 l ij （i =1,2, , m ; j =1,2, , p ）。 i i =1,2, , m
系数l ij 的确定原则：
1、z i 与z j (i ≠j ; i , j =1,2, , m ) 相互确定无关，
2、z 1是x 1, x 2, , x p 的一切线性组合中方差最大者，z 2是与z 1不相关的x 1, x 2, , x p 的所有线性组合中方差最大者„z m 是与z 1, z 2, , z m -1都不相关的x 1, x 2, , x p 的所有线性组合中方差最大者。
则新变量指标z 1, z 2, , z m 分别称为原变量指标x 1, x 2, , x p 的第一，第二，„，第m 主成分。
从数学上容易知道，从数学上可以证明，它们分别是的相关矩阵的m 个较大的特征值所对应的特征向量。
5.2.2 主成分分析中相关量的计算步骤 (1)相关系数矩阵的计算
为原变量x i 与x j 的相关系数，r ij =r ji ，其计算公式为：
r ij （i , j =1,2, , p ）
r ij =
∑(x
n
ki
-i )(x kj -j )
(2)特征值的计算
解特征方程λI -R =0，常用雅可比法（Jacobi ）求出特征值，并使其按大小顺序排列λ1≥λ2≥ , ≥λp ≥0； (3)计算主成分贡献率及累计贡献率
贡献率:
λi
∑λ
k =1
p
(i =1, 2, , p )
k
累计贡献率：
∑λ∑λ
k =1k =1p
i
k
(i =1, 2, , p )
k
一般取累计贡献率达85%~95%的特征值所对应的第一、第二、„、第m （m ≤p ）个主成分。
(4)计算主成分载荷
l ij =p (z i , x j ) =ij (i , j =1,2, , p ) (5)计算各主成分的得分：
⎡z 11z 12 z 1m ⎤⎢z ⎥z z 21222m ⎥ Z =⎢⎢ ⎥⎢⎥z z z nm ⎦⎣n 1n 2
运用主成分分析的方法，对5个区域分别进行主成分分析，本文以第一区（生活区）为例，对8种重金属进行主成分分析，具体结果如下：
5.2.2.1 运用Spss 得到相关系数矩阵
在两种金属元素的相关系数大于0.5时，表示这两种金属的相关性很强。由相关系数矩阵可知，变量之间存在较强的相关关系，适合做主成分分析。 5.2.2.2 得到8个原指标的解释总方差
合用此4个做主成分分析。 5.2.2.3 得到成分矩阵
上表中每一列表示一个主成分作为原来变量线性组合的系数，也就是主成分分析模型中的系数l ij 。
比如第一个主成分所在列的第一个系数0.655，表示第一个主成分和原来的第一个变量As 之间的线性相关系数，这个系数越大表明主成分对该变量的代表性越大。
根据主成分分析模型和因子载荷，可以得到4个主成份与原来8个变量之间的线性组合表达式如下：
⎧Y 1=0.655x 1+0.634x 2+0.732x 3+0.699x 4+0.513x 5+0.599x 6+0.743x 7+0.860x 8⎪Y =0.491x -0.513x +0.359x +0.131x -0.383x +0.651x -0.453x -0.242x ⎪212345678
⎨
Y =0.302x -0.096x -0.224x -0.368x +0.695x +0.2x -0.076x -0.144x 12345678⎪3⎪⎩Y 4=0.334x 1+0.124x 2-0.477x 3+0.444x 4-0.066x 5-0.067x 6+0.082x 7-0.285x 8在主成分表达式中，由于Zn 在Y 1中比例最大，Ni 在Y 2中比例最大，Hg 在Y 3中比例最大且三者相关系数都超过了0.5，又Y 4中没有相关系数超过0.5的变量，故得到以下结论：
在生活区，Zn , Ni , Hg 为主要污染元素。
同理，对其它4个区进行Spss 分析，得到如下结论：
在工业区，Cu , As , Hg , Pb 为主要污染元素。在山区，Zn , As 为主要污染元素。在主干道路区，Zn , As , Hg , Cd 为主要污染元素，在公园绿地区，Zn , Cu , Hg 为主要污染元素。 造成这一污染的主要原因如下：
重金属具有富集性，很难在环境中降解且Zn , Hg , Ni , Cu , Pb , Cr , Cd 具有很强的生物毒性。由于工业生产产生大量废物，汽车尾气排放，汽车轮胎磨损产生大量含重金属的有害气体和粉尘，Cu , As , Hg , Pb 主要在工业区富集，随时间的推移，工业区的土壤重金属具有很强的叠加性。Zn , Ni , Hg 主要在生活区富集，生活区内由于含重金属的废弃物堆积，施用含Zn , Ni , Hg 等的农药，不合理的施用化肥以及大量使用农用塑料薄膜都会造成重金属污染，综合风向等多种因素的影响可能会发生迁移且富集，从而造成公园绿地区的重金属污染。Zn , As , Hg , Cd 主要在主干道路区富集, 交通区如公路、铁路等交通设施的两侧汽车尾气排放及汽车轮胎磨损产生的大量含重金属的有害气体和粉尘的沉降，而污染元素则主要为Zn , As , Hg , Cd 。
综合考虑全城区的重金属污染状况，可以发现八种重金属的浓度含量与远离人群和工业活动地区所取的背景值相比较均有较大幅度的提升，其中以铜、汞和锌三种重金属
浓度含量的提升幅度最为明显，表明该城区内可能有铜锌矿的开采和冶炼，金属加工，机械制造和钢铁生产等有关工业活动，另外从汞的浓度变化幅度值我们也可以推测该区的工业处理水很大一部分被用在了灌溉农作物和养殖水生生物。
另外从每种重金属元素在该城市的不同区的分布情况来看，砷元素、镉元素、镍元素和铅元素在每个区的浓度分布和离散性大致是均匀的，表明这些元素可能通过大气金属沉降的的方式来对该城市的五个区进行污染。我们再看铬元素、铜元素、汞元素和锌元素的浓度分布情况，它们各自在五个区内的浓度分布情况不同，这表明工业区中有多家工厂，它们对工业‘三废’的处理并不到位，以至于工厂周边坏境较严重地受到了重金属污染。
从重金属污染物的来源考虑：工业过程中广泛使用重金属元素，工矿企业将未经严格处理的废水直接排放，使得它们周围的土壤容易富集高含量的有毒重金属。其次，来自于市区铜矿工业厂、硫酸厂、油漆厂、采矿和化学工业产生大量废气及粉尘微粒，由于风的输送，这些细微颗粒的Pb ，从工业废物堆扩散至周围地区，最终通过自然沉降和雨淋沉降进入土壤。另外，矿业和工业固体废弃物在堆放或处理过程中，由于日晒、雨淋、水洗等，重金属极易移动，以辐射状、漏斗状向周围土壤扩散。农业生产，尤其是近代农业生产过程中含重金属的化肥、有机肥、城市废弃物和农药的不合理施用以及污水灌溉等，都可以导致土壤中重金属的污染。化肥中品位较差的过磷酸钙和磷矿粉中含有微量的As , Cd 重金属元素。与传统的有机肥肥源相比，当前有机肥肥源大多来源于集约化的养殖场，大多使用饲料添加剂。目前的饲料添加剂中常含有高含量的Cu 和Zn ，这使得有机肥料中的Cu 、Zn 含量也明显增加并随着肥料施入农田。许多农用化学品如Cu 制剂，含Hg 、As 的制剂使用后也会使土壤遭受污染。利用污水灌溉已成为农业灌溉用水的重要组成部分。农用塑料薄膜生产应用的热稳定剂中含有Cd 、Pb ，在大量使用塑料大棚和地膜过程中都可以造成土壤重金属的污染。可以得出农业生产生活上化肥、农药农用化学品的施用是土壤重金属污染的主要途径之一。最后就是交通运输对土壤重金属污染的影响。道路两侧土壤中的污染物主要来自汽车尾气排放及汽车轮胎磨损产生的大量含重金属的有害气体和粉尘的沉降，而污染元素则主要为Pb 、Cu 、Zn 等元素 。它们一般以道路为中心成条带状分布，强度因距离公路、铁路、城市以及交通量的大小有明显的差异。
5.3 问题三的求解
在实际生活应用中，人们主要利用的金属并不是单质，而是各种合金和离子溶液。而污染环境的各种金属物质在各种介质中不是固定不变的，它们的形态，位置，浓度往往会随着时间的变化而发生改变，从而使污染扩散。在第三问中，本文结合菲克定律对重金属污染的传播进行分析。
菲克第一定律指出，在稳态扩散系的条件下，单位时间内通过垂直于扩散方向的单位面积的扩散物质量（通称扩散通量）与该截面处的浓度梯度成正比。为简便起见，仅考虑单向扩散问题。设扩散沿x 轴方向进行，菲克第一定律的表达式为：
dm ∂C
=-D () （*） Adt ∂x
∂C
J =-D
∂x
公式（*）表示在∆t 时间，沿x 方向通过x 处截面所迁移的物质的量∆m 与x 处，浓度梯度成正比。
∂C
式中：J 为扩散通量(atoms /(m 2⋅s )) 或kg /(m 2⋅s )) ；D 为扩散系数(m 2⋅s ) , 为浓
∂x
度梯度(atoms /(m 3⋅m ) 或kg /(m 3⋅m )) ；“-”号表示扩散方向为浓度梯度的反方向，即扩散由高浓度向低浓度区进行。此方程又称为扩散第一方程。
实际上，大多数重要的扩散是非稳态的，在扩散过程中扩散物质的浓度随时间而变化，即dc /dx ≠0。为了研究这种情况，根据扩散物质的质量平衡，在菲克第一定律的基础上推导出了菲克第二定律，用以分析非稳态扩散。在球面坐标系中，菲克第二定律的表达式为
∂C D ∂2∂C
=2(r ) ∂t r ∂r ∂r
在第三问中考虑到重金属污染在空气中传播量极少，故在本文中忽略了空气传播的影响，将整个城市理想化为一个平面图，并建立极坐标系。原直角坐标系的原点即为极坐标系的原点。假设污染源的位置为O 0(r 0, θ0) ，受污染点的位置为O i (r i , θi ) ，结合上面
⎧⎪C r =r i =C i
公式，两点的重金属浓度分别为C 0, C 1， 且边界条件⎨
C =C 0⎪⎩r =r 0
∂C
=cos θ=a 则有 r 2
∂r
根据边界条件，确定待定常数a , b
r r (C -C i ) a =0i 0
r 1-r 0
C r -C 1r i b =00
r 0-r i a
解得 C =-+b
r
求得浓度分布
rr (C -C i ) C 0r 0-C i r i
C (r ) =-i 00 （1） +
r (r 0-r i ) r 0-r i
dm ∂C
=a 的关在实际中，往往需要求出单位时间内通过路程的扩散量，并利用r 2
∂r dt
系
C 0-C i dm dC
（2） =JA =-D ⋅4πr 2=4πDa =4πDrr i 0
dt dr r 0-r i
一般说来，这种相变速度较慢，而且涉及的范围较广，因此可将重金属的扩撒过程当作准稳态扩散处理。由球对称稳态扩散的分析结果式（2），并利用C 0 >>C 1，即污染
dm
源很小、扩散范围很大的条件。由扩散通过路程的流量，其负值即为重金属的扩散
dt
速率。
C -C i C -C i C -C i dm
=-D ⋅4π⋅r 1r 00≈-D ⋅4π⋅r 020=-D ⋅4π⋅r 020
dt r 0-r i r i r i
将个别污染异常点剔除后将各点的坐标转换为极坐标，将各点的测量值带入方程。当t →0 时，即初状态时，即可得出污染源。不过由于在剔除了部分异常点，故需进行多次计算，最终可能得到一个或数个污染源。将8种重金属的测量值带入，最后得到了如下可能是污染源的坐标：
As （17087,11933），（18134,10046） Cd （7405,10981），（20951,10046） Cr （4592,4603）， Cu （3299,6018）
Hg （2383,2692），（13694,2357） Ni （4592,4603） Pb （1647,2728）
Zn （2383,3692），（14173,11941）
5.4 问题四的求解
缺点：问题三中的模型将重金属污染的扩散理想化为均匀介质中的扩散，忽略了实际情况中介质的多样化和叠加情况。忽略了海拔因素的影响，在实际情况中污染往往从地势较高的地区向地势较低的地区扩散。
优点：运用科学的模型将重金属扩散简化为相似的化学扩散形式，使问题简化，方便对污染源的确定。将坐标系转化为极坐标系，使污染源和被污染的地点关系更加直观。
为了研究城市地质环境的演变模式，必须同时考虑自然因素和社会因素。自然因素包括了：地壳运动，河流因素，气候变化等。人为因素包括了农业因素，生活垃圾，工业排污等。
图5 地质因素
本题已知的重金属污染在上题的分析中可以得出，在工业区的污染最为严重，故可将本题中的重金属污染视作工业排污。
还需要搜集的信息主要是当地的土壤样本，河流的布局及流量，当地气候情况，垃圾处理方式，农业规划力度。
通过网上调查可以得到系列相关参数如下表所示
再通过专家评估和自己判断，对每个参数进行赋权分别为w 1, w 2, w 3, w 4, w 5, w 6, 。 最后根据实际测得的数据进行求和得：
G i =∑w j y ij
j =1
6
由此便可得出不同时期各种因素对该城市的地质环境的影响，从而得出其演变模式。
六、模型的检验
在对问题三污染源的位置进行判定时，当时采用的是菲克定律，现在以金属As 为例，用“二次筛选法”判断造成As 污染的污染源位，检验问题三结果的正确性。
首先用Excel 做出As 污染与海拔的散点图，通过对散点图的分析得知：As 正常浓度在1.8~5.4，由上图，海拔在100米之上时，As 的浓度几乎都在正常值，可以认为受As 污染的地方就是在海拔100米之下，可以认为As 由低处向高出传播是比较难的。主要是由于空气传播或者人为污染，而只有人为污染和其他特殊情况，才会使重金属向海拔高的地方传播。
两次筛选数据，确定污染源
（1）第一次筛选数据，由图8分析及重金属传播特性可知，污染源在海拔100米以下的地区，因此本文只截取海拔在100米以下的数据来确定污染源。
（2）第二次筛选数据，用Matlab 画出海拔100米以下数据As 污染浓度的等高线图如下
图6 As 污染浓度的等高线图
由图9可以看出，污染源最有可能在：
(3000~4000,6000~9000),(17000~19000,9000~11000),
因为在这几个地方，As 的浓度最高，并且离这几个点越远，As 的浓度越小。
在筛选后的区域(3000~4000,6000~9000),(17000~19000,9000~11000),分别取这样的数据：坐标(x , y ), 海拔和As 的浓度。本文将选出的数据列表如下：
表6 坐标(x , y ), 海拔和As 的浓度
x 10
4
然后我们做出这12个的采样点的浓度变化曲线如下
图7 12个的采样点的浓度变化曲线
从图10可以看出6，9，12可能是污染源，但是再次考虑到污染物的传播特征，低洼处的污染可能是由于海拔较周围地区低，周围地区普遍受污染且二次污染严重，污染物大都传播到这里无法排出而在这大量积累，造成这个地方As 的浓度也会非常高，所以某些As 浓度高处也可能不是污染源。又因为编号12周围大部分地区都比这个地区海拔高，所以不能确定它是不是污染源，而编号6，9处却比周围大多地区海拔高，所以可以说6，9处有As 的污染源。即(17087,11933),(18134,10046)是污染源。
与问题三的结论一致，从而检验了问题三结果的合理性。其他7重金属同理。
七、模型的评价与改进
7.1模型的优点
1、在对问题一进行污染程度的求解前，对附表给出的319组数据中的20个异常点进行剔除，使得得出的5个区的污染程度较精确。
2、运用内梅罗指数的方法，数学过程简洁，运算方便。物理概念清晰，对于一个评价区，我们只计算出它的综合指数，再按照相应的分级标准，便可知道该评价区某环境要素的综合环境质量状况，便于决策者做出综合决策。
3、问题二采用了主成分分析的数学模型，较为合理的得出了影响生活区，工业区，山区，主干道路区，公园绿地区的重金属元素。
4、合理恰当地使用了表格和图形来代替文字说明，使数据的体现和意思的表达更加清晰。
5、由于本题的数据量较大，所以利用Matlab 来处理数据，方便准确。
7.2模型的缺点
1、本题中有些重金属在水中还是有微弱的反应的，故模型的建立和求解出来的答案与实际有一定的差距。
2、由于题中没有说明各个区域中的主要工商业，故在确定污染源的来源时不能准确定位。
3、内梅罗指数的方法过分突出极大值对污染的影响，未考虑次大值对污染的影响，未考虑危害性最大的污染因子对水质的影响。
7.3 模型的改进
对于问题一，本文采用的是传统的尼梅罗综合污染指数法, 而传统尼梅罗综合指数过分突出污染指数最大的重金属污染物对环境质量的影响和作用，在评价时可能会人为地夸大或缩小一些因子的影响作用，使其对环境质量评价的灵敏性不够高，在某些情况下，它的计算结果难以区分土壤环境污染程度的差别。
所以现在对其进行改进，改进如下：
F max +Fw '
计算公式为P mj =F ' max =，Fw 为权重值最大的污染因
2子的c i 与s ij 比值，其他符号意义同第一问。
改进的尼梅罗指数法具有如下优点：
1、以国家最新土壤环境质量标准划分污染等级。污染等级的划分随标准改变而改变，具有划时代性。
2、待测土壤尼梅罗指数的计算方法，与划分土壤污染等级时各类标准土壤尼梅罗指数的计算方法相同，所采用的基准也一样，具有同一性。 3、考虑了土壤中危害性最大的污染因子对土壤的影响。
八、模型的推广
本模型的建立对于研究减少城市污染问题和保护环境具有重要意义，尤其在当今以高能耗高污染的生产模式为主的工业时代，该模型的建立对于研究城市规划，和工厂位置选择以及交通建设时具有重要的参考价值，同时利用该模型也可以研究物质扩散现象的规律。
九、参考文献
[1] 杨启航, 数学建模, 北京:高等教育出版社,2003；
[2] 韩中庚，数学建模方法及其应用，北京高等教育出版社，2005； [3] 姜启源，谢金星，叶俊，数学模型，北京，2003；
[4] 韩中庚，数学建模方法及其应用，北京高等教育出版社，2005。
十、附录
y
00.511.522.5
x 10
4
图8 cd
y
2.5
x 10
4
0.511.52
图9 cr
y
0.5
11.522.5
x 10
4
图
10 cu
y
00.511.522.5
x 10
4
图11 hg
y
0.5
11.522.5
x 10
4
图12 Ni
y
0.5
1
1.52
2.5
x 10
4
图13 Pb
y
1
1、去除20个异常点的程序 clc
clear all
X=xlsread('data1.xls'); %导入数据
%%%%%取出单位是ug/g的六个污染元素的浓度 A1=X(:,1); A3=X(:,3); A4=X(:,4); A6=X(:,6); A7=X(:,7); A8=X(:,8);
%%%%%%处理, 找出浓度大于500ug/g的异常值并去除 k1=find(A1>500) A1(k1)=0;
k3=find(A3>500) A3(k3)=0; length(A3)
k4=find(A4>500) A4(k4)=0;
k6=find(A6>500) A6(k6)=0;
k7=find(A7>500) A7(k7)=0;
k8=find(A8>500) A8(k8)=0;
%%%%%取出单位是ng/g的六个污染元素的浓度 A2=X(:,2); A5=X(:,5);
0.51.522.5
x 10
4
图14 Zn
%%%%%%处理, 找出浓度大于4000ng/g的异常值并去除 k2=find(A2>4000)
A2(k2)=0;
k5=find(A5>4000)
A5(k5)=0;
X=[A1,A2,A3,A4,A5,A6,A7,A8]
xlswrite('chuli1.xls',X) %导出数据
2、319个选样点的三维地势图程序
clc
clear all
close all
d=xlsread ('data.xls');%导入数据
x=d (:,2);%采样点坐标x 值
y=d (:,3);%采样点坐标y 值
z=d (:,4);%采样点坐标z(海拔) 值
nx=linspace(min(x),max(x),100);
ny=linspace(min(y),max(y),100);
[xx,yy]=meshgrid (nx,ny);
zz=griddata(x,y,z,xx,yy,'v4');
surfl (xx,yy,zz);
shading interp
set (gcf,'color','w');
colormap(cool);
hold on
for i=1:44
plot3 (d(i,2),d(i,3),d(i,4),'g*')
end
hold on
for i=45:80
plot3(d(i,2),d(i,3),d (i,4),'r+')
end
hold on
for i=81:146
plot3(d(i,2),d(i,3),d(i,4),'w+');
end
hold on
for i=147:284
plot3(d(i,2),d(i,3),d(i,4),'y+');
end
hold on
for i=285:319
plot3(d(i,2),d(i,3),d(i,4),'k*');
end
title('选样点的三维地势图')
20
3、问题二中八种金属等高线图
clc
clear all
close all
A=xlsread('data2');
x=A(:,1);y=A(:,2);
[x,y]=meshgrid(0:1000:29000,100:300:19000); z3=griddata(A(:,1),A(:,2),A(:,3),x,y,'v4');
figure(1)
[c,h] =contour(x,y,z3);
clabel(c,h)
ylabel('y');
xlabel('图1 As的空间分布');
z4=griddata(A(:,1),A(:,2),A(:,4),x,y,'v4');
figure(2)
[c,h] =contour(x,y,z4);
clabel(c,h)
ylabel('y');
xlabel('图2 Cd的空间分布')
z5=griddata(A(:,1),A(:,2),A(:,5),x,y,'v4');
figure(3)
[c,h] =contour(x,y,z5);
clabel(c,h)
ylabel('y');
xlabel('图3 Cr的空间分布')
z6=griddata(A(:,1),A(:,2),A(:,6),x,y,'v4');
figure(4)
[c,h] =contour(x,y,z6);
clabel(c,h)
ylabel('y');
xlabel('图4 Cu的空间分布')
figure(5)
z7=griddata(A(:,1),A(:,2),A(:,7),x,y,'v4');
figure(5)
[c,h] =contour(x,y,z7);
clabel(c,h)
ylabel('y');
xlabel('图5 Hg的空间分布')
z8=griddata(A(:,1),A(:,2),A(:,8),x,y,'v4');
figure(6)
[c,h] =contour(x,y,z8);
clabel(c,h)
ylabel('y');
xlabel('图6 Ni的空间分布')
21
z9=griddata(A(:,1),A(:,2),A(:,9),x,y,'v4');
figure(7)
[c,h] =contour(x,y,z9);
clabel(c,h)
ylabel('y');
xlabel('图7 Pb的空间分布')
z10=griddata(A(:,1),A(:,2),A(:,10),x,y,'v4');
figure(8)
[c,h] =contour(x,y,z10);
clabel(c,h)
ylabel('y');
xlabel('图8 Zn的空间分布')
4、模型检验中海拔小于100米As 等高线图 clc
clear all
close all
A=xlsread('data2');
xx=A(:,1);yy=A(:,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% T1=A(:,1);T2=A(:,2);T3=A(:,3);
k=find(T3
U1=T1(k);U2=T2(k);U3=T3(k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[xx,yy]=meshgrid(0:1000:29000,100:300:19000); %%%%第一种金属污染源
figure(1)
z3=griddata(T1(k),T2(k),T3(k),xx,yy,'v4');
[c,h] =contour(xx,yy,z3);
clabel(c,h)
22
城市表层土壤重金属污染分析
摘要
本文以城市表层土壤重金属污染为研究对象，分析了8种主要重金属元素在该城区的空间分布、污染程度和原因，建立模型确定污染源的位置，并对模型改进使之能更好的研究城市地质环境的演变模式。
对于问题一，首先对数据预处理剔除异常点并进行无量纲化处理，绘出8种重金属的海拔图和等高线图，得出空间分布为：重金属污染集中在海拔较低处，污染程度随海拔的升高而降低，铬、铜、汞、铅、锌的分布集中在一起。镉分布较广，各功能区分布不一等结论。然后建立尼梅罗模型，由第m 区i 组j 种重金属实际值c ij ，j 种重金属背景值的平均值s j ，得到单因子污染指数F ij =c ij /s j ，求得尼梅罗污染指数为：
P mj =
根据此式得到不同区域各个重金属的尼梅罗污染指数，对每个区域的8种重金属的尼梅罗污染指数求和，得到结果如下： 生活区1.2143，工业区1.4885，山区0.4007，主干道路区0.6715，公园绿地区0.4370。 得出污染的程度为：工业区>生活区>主干道路区>公园绿地区>山区
对于问题二，建立主成分分析模型，对5个区域分别进行主成分分析，运用Spss 得到相关系数矩阵，可知8个变量适合做主成分分析，然后对8个原指标进行解释总方差计算，提取四个主成分，计算成分矩阵，得出结论：
生活区，Zn , Ni , Hg 为主要污染元素。在工业区，Cu , As , Hg , Pb 为主要污染元素。在山区，Zn , As 为主要污染元素。在主干道路区，Zn , As , Hg , Cd 为主要污染元素，在公园绿地区，Zn , Cu , Hg 为主要污染元素。
对于问题三，采用菲克定律，得到重金属速率随时间关系：
C -C i C -C i C -C i dm
=-D ⋅4π⋅r 1r 00≈-D ⋅4π⋅r 020=-D ⋅4π⋅r 020
dt r 0-r i r i r i
将8种重金属的测量值带入，最后得到了如下可能是污染源的坐标：As （17087,11933），
Cr （4592,4603）Cu （3299,6018）（18134,10046）Cd （7405,10981），（20951,10046），
Hg （2383,2692），（13694,2357）Ni （4592,4603）Pb （1647,2728）Zn （2383,3692），（14173,11941）
对于问题四，结合实际情况，分析城市地质环境的演变模式。本文将能影响城市的因素先分为两大类：自然因素和人为因素。再将其各细分为三个次级因素。通过网络调查可以得出各因素在不同时期的实际参数，再通过专家评估和自己的判断，得出每种因素对地质环境演变的影响程度，并对其赋权。最后建立模型得到其演变模式。
关键词 重金属污染 尼梅罗指数 主成分分析 菲克定律
一、问题重述
1.1 背景分析
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。重金属的来源非常广泛，传统上可以分为工业来源和农业来源。随着我国城市化进程的加快，一些有往于城市所特有的污染也随之产生。另外近几年来突发性环境事故频发，重金属污染事件也层出不穷，环境事故也成为重金属污染的重要来源之一。城市工业“三废”排放，金属采矿和冶炼，家庭燃煤，生活垃圾，汽车尾气排放都增加了城市土壤重金属的负荷。重金属污染环境的主要有汞、铅、铬、锌镉、铜等。其中汞的毒性最大，铬、铅、锌等也有相当大毒性。此外还有砷，砷虽不属于金属. 但它的毒性与重金属相似，因此归于重金属一类阐述，称为类金属。目前对我国土壤污染比较普遍的重金属有汞、铬、砷。
对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
1.2 需要解决的问题
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、模型假设
结合本题的实际，为了确保模型求解的准确性和合理性，我们排除了一些位置因素的干扰，提出以下几点假设：
1、假设重金属的物理化学性质稳定。 2、题目所给的数据准确无误，客观真实。
3、假设在采样前该城区没有发生重金属泄露事故。 4、城市地表土壤中的重金属污染物含量成正态分布。
三、符号说明
为了便于问题的求解，我们给出以下符号说明：（其他未说明的符号在文中第一次
四、问题分析
4.1 问题的重要性
土壤重金属污染在一定时期内不表现出对环境的危害性，当其含量超过土壤承受力或限度，或土壤环境变化时，重金属有可能突然活化，引起严重的生态危害。重金属对土壤中的农作物有很强的毒害作用，除此之外，土壤中的重金属还会经由雨水淋滤及地表径流作用转移进入地表水系统，对饮用水安全构成威胁。对城市土壤地质环境异常的查证，成为人们关注的问题。 4.2 问题的思路分析
4.2.1 问题一的分析
对于问题一，本文首先作出了319个采样点的三维地势图，随后为了准确分析第一、二、三、四、五区的重金属污染程度，剔除了20个异常点，随后建立了尼梅罗指数法的模型，得出了5个区环境污染程度的级别。
4.2.2 问题二的分析
问题二要求通过数据分析，说明重金属污染的主要原因。本文再次建立了主成份分析法的数学模型，对5个区内的8种重金属对环境影响的孰强孰弱进行了研究，得到了5个区对环境影响最大的重金属元素。再结合相关文献资料，得出5个区内重金属污染的原因。
4.2.3 问题三的分析
问题三中需要分析重金属污染的传播特征，而在实际生活中，金属的应用主要是合金和离子溶液。往往是通过化学反应来实现其应用，故在问题三中，本文将重金属污染的传播理想化为化学扩散，运用菲克定律，分析不同浓度的地区之间重金属的扩散，追本溯源，从而分析出污染源。
4.2.4 问题四的分析
问题四需要先对问题三中的模型进行评析，得出其结构简单，便于理解，却又过于理想化的特点。然后结合实际情况，分析城市地质环境的演变模式。本文将能影响城市
的因素先分为两大类：自然因素和人为因素。再将其各细分为三个次级因素。通过网络调查可以得出各因素在不同时期的实际参数，再通过专家评估和自己的判断，得出每种因素对地质环境演变的影响程度，并对其赋权。最后建立模型得到其演变模式。
五、模型的建立与求解
经过以上的分析和准备，我们将逐步建立以下数学模型，进一步阐述模型的实际建立过程。
5.1 问题一的求解
5.1.1 8种重金属异常点的处理
对数据预处理，根据8种重金属的浓度随采样点的分布，画出数据图。由于砷、铬、铜、镍、铅、锌这6种重金属的浓度单位为ug /g ，镉、、汞的单位为ng /g ，考虑到不同单位的值若反映在同一图内会有较大误差，故分别作出如下两图：
图1 砷等6金属浓度随采样点的分布 图2 镉、汞浓度随采样点的分布 从两张图可以看出，砷、铬、铜、镍、铅、锌的浓度不服从正态分布，且有多个异常点，为了问题一计算污染程度的准确性，本文将浓度大于500ug /g 的点剔除掉；同理，铬、汞的分布也不服从正态分布并存在多个异常点，为了问题一计算污染程度的准确性，将浓度大于4000ng /g 的采样点剔除。 5.1.2 元素浓度的无量纲化处理
在利用Spss 统计软件数据进行聚类分析的时候，因为单位不统一需要进行无量纲化处理，我们采用均值化方法，即每一个变量除以该变量的平均值，即
x x i '=i
标准化以后各变量的平均值都为1，标准差为原始变量的变异系数。该方法在消除量纲和数量级影响的同时，保留了各变量取值差异程度上的信息，差异程度越大的变量对综合分析的影响也越大。
5.1.3 8种重金属在该城区的空间分布 对数据预处理，进行异常值的过滤，用Matlab 绘出对8种重金属元素的取点的坐标和与其相对应的海拔值的3D 图，以便定性分析重金属分布与海拔高度的关系。
图3 重金属分布图
从图中可知，重金属污染集中在海拔较低处，污染程度随海拔的升高而降低。
根据8种重金属元素浓度值与坐标关系，画出各个重金属元素的二维等高线图，以As 元素为例，其余7种重金属见附录。
y
1
由各个重金属的分布可以得出：
1、铬、铜、汞、铅、锌的分布集中在一起。 2、镉分布较广，各功能区分布不一。 3、工业区呈现混合污染的态势。
5.1.4 不同区域重金属的污染程度
土壤重金属评价方法包括生态危害指数法、地积累指数法、单因子尼梅罗指数法等。由于运用尼梅罗指数的方法，数学过程简洁，加权过程中避免了权系数中主观因素的影响，故采用此方法来评价土壤重金属的污染程度。
0.5
1.5
2
2.5
x 10
4
图4 As 的空间分布
5.1.5 单因子污染指数
首先用给出的数据，计算出单因子污染指数，确定单个金属对各城区土壤的污染，其公式如下:
F ij =c ij /s j
其中定义c ij 第m 区i 组j 种重金属实际值，s j 重金属背景值的平均值，F ij 为第m 区土壤j 种重金属单项污染指数。
5.1.6 尼梅罗综合污染指数
继而求出尼梅罗综合污染指数如下：
P mj =
其中P mj 为在m 区j 种重金属污染物的尼梅罗指数, F max =max(F ij ) 为监测点土壤重
金属最大单项污染指数, F =
为监测点所有重金属单项污染指数平均值。 n
求得不同区域重金属的尼梅罗污染指数结果如下：
∑F
i =1
n
ij
通过对同一个区域的8种重金属元素尼梅罗污染指数求和∑P mj
j =1
(m =1,2,3,4,5) ，
得到结果如下：
未受污染。
故得到结果如下：
污染的程度为：工业区>生活区>主干道路区>公园绿地区>山区, 生活区、工业区的污染较严重，主干道路区已受到污染，山区、公园绿地区未收到污染。
5.2 问题二主成分分析模型的建立与求解
5.2.1 模型的准备
主成分分析是把原来多个变量化为少数几个综合指标的一种统计分析方法，从数学角度来看，这是一种降维处理技术。假定有n 个样本，每个样本共有p 个变量，构成一个n ⋅p 阶的数据矩阵。
⎡x 11x 12 x 1p ⎤⎢x ⎥x x 21222p ⎥ X =⎢⎢ ⎥⎢⎥x x x ⎢np ⎥⎣n 1n 2⎦
记x 1, x 2, x p 为原变量指标，z 1, z 2, z m m ≤p 为新变量指标，得到新指标为：
⎧z 1=l 11x 1+l 12x 2+ +l 1p x p ⎪z =l x +l x + +l x ⎪22112222p p
⎨
⎪
⎪⎩z m =l m 1x 1+l m 2x 2+ +l mp x p
从以上的分析可以看出，主成分分析的实质就是确定原来变量x j （j =1,2, , p ）在主）成分z （上的荷载 l ij （i =1,2, , m ; j =1,2, , p ）。 i i =1,2, , m
系数l ij 的确定原则：
1、z i 与z j (i ≠j ; i , j =1,2, , m ) 相互确定无关，
2、z 1是x 1, x 2, , x p 的一切线性组合中方差最大者，z 2是与z 1不相关的x 1, x 2, , x p 的所有线性组合中方差最大者„z m 是与z 1, z 2, , z m -1都不相关的x 1, x 2, , x p 的所有线性组合中方差最大者。
则新变量指标z 1, z 2, , z m 分别称为原变量指标x 1, x 2, , x p 的第一，第二，„，第m 主成分。
从数学上容易知道，从数学上可以证明，它们分别是的相关矩阵的m 个较大的特征值所对应的特征向量。
5.2.2 主成分分析中相关量的计算步骤 (1)相关系数矩阵的计算
为原变量x i 与x j 的相关系数，r ij =r ji ，其计算公式为：
r ij （i , j =1,2, , p ）
r ij =
∑(x
n
ki
-i )(x kj -j )
(2)特征值的计算
解特征方程λI -R =0，常用雅可比法（Jacobi ）求出特征值，并使其按大小顺序排列λ1≥λ2≥ , ≥λp ≥0； (3)计算主成分贡献率及累计贡献率
贡献率:
λi
∑λ
k =1
p
(i =1, 2, , p )
k
累计贡献率：
∑λ∑λ
k =1k =1p
i
k
(i =1, 2, , p )
k
一般取累计贡献率达85%~95%的特征值所对应的第一、第二、„、第m （m ≤p ）个主成分。
(4)计算主成分载荷
l ij =p (z i , x j ) =ij (i , j =1,2, , p ) (5)计算各主成分的得分：
⎡z 11z 12 z 1m ⎤⎢z ⎥z z 21222m ⎥ Z =⎢⎢ ⎥⎢⎥z z z nm ⎦⎣n 1n 2
运用主成分分析的方法，对5个区域分别进行主成分分析，本文以第一区（生活区）为例，对8种重金属进行主成分分析，具体结果如下：
5.2.2.1 运用Spss 得到相关系数矩阵
在两种金属元素的相关系数大于0.5时，表示这两种金属的相关性很强。由相关系数矩阵可知，变量之间存在较强的相关关系，适合做主成分分析。 5.2.2.2 得到8个原指标的解释总方差
合用此4个做主成分分析。 5.2.2.3 得到成分矩阵
上表中每一列表示一个主成分作为原来变量线性组合的系数，也就是主成分分析模型中的系数l ij 。
比如第一个主成分所在列的第一个系数0.655，表示第一个主成分和原来的第一个变量As 之间的线性相关系数，这个系数越大表明主成分对该变量的代表性越大。
根据主成分分析模型和因子载荷，可以得到4个主成份与原来8个变量之间的线性组合表达式如下：
⎧Y 1=0.655x 1+0.634x 2+0.732x 3+0.699x 4+0.513x 5+0.599x 6+0.743x 7+0.860x 8⎪Y =0.491x -0.513x +0.359x +0.131x -0.383x +0.651x -0.453x -0.242x ⎪212345678
⎨
Y =0.302x -0.096x -0.224x -0.368x +0.695x +0.2x -0.076x -0.144x 12345678⎪3⎪⎩Y 4=0.334x 1+0.124x 2-0.477x 3+0.444x 4-0.066x 5-0.067x 6+0.082x 7-0.285x 8在主成分表达式中，由于Zn 在Y 1中比例最大，Ni 在Y 2中比例最大，Hg 在Y 3中比例最大且三者相关系数都超过了0.5，又Y 4中没有相关系数超过0.5的变量，故得到以下结论：
在生活区，Zn , Ni , Hg 为主要污染元素。
同理，对其它4个区进行Spss 分析，得到如下结论：
在工业区，Cu , As , Hg , Pb 为主要污染元素。在山区，Zn , As 为主要污染元素。在主干道路区，Zn , As , Hg , Cd 为主要污染元素，在公园绿地区，Zn , Cu , Hg 为主要污染元素。 造成这一污染的主要原因如下：
重金属具有富集性，很难在环境中降解且Zn , Hg , Ni , Cu , Pb , Cr , Cd 具有很强的生物毒性。由于工业生产产生大量废物，汽车尾气排放，汽车轮胎磨损产生大量含重金属的有害气体和粉尘，Cu , As , Hg , Pb 主要在工业区富集，随时间的推移，工业区的土壤重金属具有很强的叠加性。Zn , Ni , Hg 主要在生活区富集，生活区内由于含重金属的废弃物堆积，施用含Zn , Ni , Hg 等的农药，不合理的施用化肥以及大量使用农用塑料薄膜都会造成重金属污染，综合风向等多种因素的影响可能会发生迁移且富集，从而造成公园绿地区的重金属污染。Zn , As , Hg , Cd 主要在主干道路区富集, 交通区如公路、铁路等交通设施的两侧汽车尾气排放及汽车轮胎磨损产生的大量含重金属的有害气体和粉尘的沉降，而污染元素则主要为Zn , As , Hg , Cd 。
综合考虑全城区的重金属污染状况，可以发现八种重金属的浓度含量与远离人群和工业活动地区所取的背景值相比较均有较大幅度的提升，其中以铜、汞和锌三种重金属
浓度含量的提升幅度最为明显，表明该城区内可能有铜锌矿的开采和冶炼，金属加工，机械制造和钢铁生产等有关工业活动，另外从汞的浓度变化幅度值我们也可以推测该区的工业处理水很大一部分被用在了灌溉农作物和养殖水生生物。
另外从每种重金属元素在该城市的不同区的分布情况来看，砷元素、镉元素、镍元素和铅元素在每个区的浓度分布和离散性大致是均匀的，表明这些元素可能通过大气金属沉降的的方式来对该城市的五个区进行污染。我们再看铬元素、铜元素、汞元素和锌元素的浓度分布情况，它们各自在五个区内的浓度分布情况不同，这表明工业区中有多家工厂，它们对工业‘三废’的处理并不到位，以至于工厂周边坏境较严重地受到了重金属污染。
从重金属污染物的来源考虑：工业过程中广泛使用重金属元素，工矿企业将未经严格处理的废水直接排放，使得它们周围的土壤容易富集高含量的有毒重金属。其次，来自于市区铜矿工业厂、硫酸厂、油漆厂、采矿和化学工业产生大量废气及粉尘微粒，由于风的输送，这些细微颗粒的Pb ，从工业废物堆扩散至周围地区，最终通过自然沉降和雨淋沉降进入土壤。另外，矿业和工业固体废弃物在堆放或处理过程中，由于日晒、雨淋、水洗等，重金属极易移动，以辐射状、漏斗状向周围土壤扩散。农业生产，尤其是近代农业生产过程中含重金属的化肥、有机肥、城市废弃物和农药的不合理施用以及污水灌溉等，都可以导致土壤中重金属的污染。化肥中品位较差的过磷酸钙和磷矿粉中含有微量的As , Cd 重金属元素。与传统的有机肥肥源相比，当前有机肥肥源大多来源于集约化的养殖场，大多使用饲料添加剂。目前的饲料添加剂中常含有高含量的Cu 和Zn ，这使得有机肥料中的Cu 、Zn 含量也明显增加并随着肥料施入农田。许多农用化学品如Cu 制剂，含Hg 、As 的制剂使用后也会使土壤遭受污染。利用污水灌溉已成为农业灌溉用水的重要组成部分。农用塑料薄膜生产应用的热稳定剂中含有Cd 、Pb ，在大量使用塑料大棚和地膜过程中都可以造成土壤重金属的污染。可以得出农业生产生活上化肥、农药农用化学品的施用是土壤重金属污染的主要途径之一。最后就是交通运输对土壤重金属污染的影响。道路两侧土壤中的污染物主要来自汽车尾气排放及汽车轮胎磨损产生的大量含重金属的有害气体和粉尘的沉降，而污染元素则主要为Pb 、Cu 、Zn 等元素 。它们一般以道路为中心成条带状分布，强度因距离公路、铁路、城市以及交通量的大小有明显的差异。
5.3 问题三的求解
在实际生活应用中，人们主要利用的金属并不是单质，而是各种合金和离子溶液。而污染环境的各种金属物质在各种介质中不是固定不变的，它们的形态，位置，浓度往往会随着时间的变化而发生改变，从而使污染扩散。在第三问中，本文结合菲克定律对重金属污染的传播进行分析。
菲克第一定律指出，在稳态扩散系的条件下，单位时间内通过垂直于扩散方向的单位面积的扩散物质量（通称扩散通量）与该截面处的浓度梯度成正比。为简便起见，仅考虑单向扩散问题。设扩散沿x 轴方向进行，菲克第一定律的表达式为：
dm ∂C
=-D () （*） Adt ∂x
∂C
J =-D
∂x
公式（*）表示在∆t 时间，沿x 方向通过x 处截面所迁移的物质的量∆m 与x 处，浓度梯度成正比。
∂C
式中：J 为扩散通量(atoms /(m 2⋅s )) 或kg /(m 2⋅s )) ；D 为扩散系数(m 2⋅s ) , 为浓
∂x
度梯度(atoms /(m 3⋅m ) 或kg /(m 3⋅m )) ；“-”号表示扩散方向为浓度梯度的反方向，即扩散由高浓度向低浓度区进行。此方程又称为扩散第一方程。
实际上，大多数重要的扩散是非稳态的，在扩散过程中扩散物质的浓度随时间而变化，即dc /dx ≠0。为了研究这种情况，根据扩散物质的质量平衡，在菲克第一定律的基础上推导出了菲克第二定律，用以分析非稳态扩散。在球面坐标系中，菲克第二定律的表达式为
∂C D ∂2∂C
=2(r ) ∂t r ∂r ∂r
在第三问中考虑到重金属污染在空气中传播量极少，故在本文中忽略了空气传播的影响，将整个城市理想化为一个平面图，并建立极坐标系。原直角坐标系的原点即为极坐标系的原点。假设污染源的位置为O 0(r 0, θ0) ，受污染点的位置为O i (r i , θi ) ，结合上面
⎧⎪C r =r i =C i
公式，两点的重金属浓度分别为C 0, C 1， 且边界条件⎨
C =C 0⎪⎩r =r 0
∂C
=cos θ=a 则有 r 2
∂r
根据边界条件，确定待定常数a , b
r r (C -C i ) a =0i 0
r 1-r 0
C r -C 1r i b =00
r 0-r i a
解得 C =-+b
r
求得浓度分布
rr (C -C i ) C 0r 0-C i r i
C (r ) =-i 00 （1） +
r (r 0-r i ) r 0-r i
dm ∂C
=a 的关在实际中，往往需要求出单位时间内通过路程的扩散量，并利用r 2
∂r dt
系
C 0-C i dm dC
（2） =JA =-D ⋅4πr 2=4πDa =4πDrr i 0
dt dr r 0-r i
一般说来，这种相变速度较慢，而且涉及的范围较广，因此可将重金属的扩撒过程当作准稳态扩散处理。由球对称稳态扩散的分析结果式（2），并利用C 0 >>C 1，即污染
dm
源很小、扩散范围很大的条件。由扩散通过路程的流量，其负值即为重金属的扩散
dt
速率。
C -C i C -C i C -C i dm
=-D ⋅4π⋅r 1r 00≈-D ⋅4π⋅r 020=-D ⋅4π⋅r 020
dt r 0-r i r i r i
将个别污染异常点剔除后将各点的坐标转换为极坐标，将各点的测量值带入方程。当t →0 时，即初状态时，即可得出污染源。不过由于在剔除了部分异常点，故需进行多次计算，最终可能得到一个或数个污染源。将8种重金属的测量值带入，最后得到了如下可能是污染源的坐标：
As （17087,11933），（18134,10046） Cd （7405,10981），（20951,10046） Cr （4592,4603）， Cu （3299,6018）
Hg （2383,2692），（13694,2357） Ni （4592,4603） Pb （1647,2728）
Zn （2383,3692），（14173,11941）
5.4 问题四的求解
缺点：问题三中的模型将重金属污染的扩散理想化为均匀介质中的扩散，忽略了实际情况中介质的多样化和叠加情况。忽略了海拔因素的影响，在实际情况中污染往往从地势较高的地区向地势较低的地区扩散。
优点：运用科学的模型将重金属扩散简化为相似的化学扩散形式，使问题简化，方便对污染源的确定。将坐标系转化为极坐标系，使污染源和被污染的地点关系更加直观。
为了研究城市地质环境的演变模式，必须同时考虑自然因素和社会因素。自然因素包括了：地壳运动，河流因素，气候变化等。人为因素包括了农业因素，生活垃圾，工业排污等。
图5 地质因素
本题已知的重金属污染在上题的分析中可以得出，在工业区的污染最为严重，故可将本题中的重金属污染视作工业排污。
还需要搜集的信息主要是当地的土壤样本，河流的布局及流量，当地气候情况，垃圾处理方式，农业规划力度。
通过网上调查可以得到系列相关参数如下表所示
再通过专家评估和自己判断，对每个参数进行赋权分别为w 1, w 2, w 3, w 4, w 5, w 6, 。 最后根据实际测得的数据进行求和得：
G i =∑w j y ij
j =1
6
由此便可得出不同时期各种因素对该城市的地质环境的影响，从而得出其演变模式。
六、模型的检验
在对问题三污染源的位置进行判定时，当时采用的是菲克定律，现在以金属As 为例，用“二次筛选法”判断造成As 污染的污染源位，检验问题三结果的正确性。
首先用Excel 做出As 污染与海拔的散点图，通过对散点图的分析得知：As 正常浓度在1.8~5.4，由上图，海拔在100米之上时，As 的浓度几乎都在正常值，可以认为受As 污染的地方就是在海拔100米之下，可以认为As 由低处向高出传播是比较难的。主要是由于空气传播或者人为污染，而只有人为污染和其他特殊情况，才会使重金属向海拔高的地方传播。
两次筛选数据，确定污染源
（1）第一次筛选数据，由图8分析及重金属传播特性可知，污染源在海拔100米以下的地区，因此本文只截取海拔在100米以下的数据来确定污染源。
（2）第二次筛选数据，用Matlab 画出海拔100米以下数据As 污染浓度的等高线图如下
图6 As 污染浓度的等高线图
由图9可以看出，污染源最有可能在：
(3000~4000,6000~9000),(17000~19000,9000~11000),
因为在这几个地方，As 的浓度最高，并且离这几个点越远，As 的浓度越小。
在筛选后的区域(3000~4000,6000~9000),(17000~19000,9000~11000),分别取这样的数据：坐标(x , y ), 海拔和As 的浓度。本文将选出的数据列表如下：
表6 坐标(x , y ), 海拔和As 的浓度
x 10
4
然后我们做出这12个的采样点的浓度变化曲线如下
图7 12个的采样点的浓度变化曲线
从图10可以看出6，9，12可能是污染源，但是再次考虑到污染物的传播特征，低洼处的污染可能是由于海拔较周围地区低，周围地区普遍受污染且二次污染严重，污染物大都传播到这里无法排出而在这大量积累，造成这个地方As 的浓度也会非常高，所以某些As 浓度高处也可能不是污染源。又因为编号12周围大部分地区都比这个地区海拔高，所以不能确定它是不是污染源，而编号6，9处却比周围大多地区海拔高，所以可以说6，9处有As 的污染源。即(17087,11933),(18134,10046)是污染源。
与问题三的结论一致，从而检验了问题三结果的合理性。其他7重金属同理。
七、模型的评价与改进
7.1模型的优点
1、在对问题一进行污染程度的求解前，对附表给出的319组数据中的20个异常点进行剔除，使得得出的5个区的污染程度较精确。
2、运用内梅罗指数的方法，数学过程简洁，运算方便。物理概念清晰，对于一个评价区，我们只计算出它的综合指数，再按照相应的分级标准，便可知道该评价区某环境要素的综合环境质量状况，便于决策者做出综合决策。
3、问题二采用了主成分分析的数学模型，较为合理的得出了影响生活区，工业区，山区，主干道路区，公园绿地区的重金属元素。
4、合理恰当地使用了表格和图形来代替文字说明，使数据的体现和意思的表达更加清晰。
5、由于本题的数据量较大，所以利用Matlab 来处理数据，方便准确。
7.2模型的缺点
1、本题中有些重金属在水中还是有微弱的反应的，故模型的建立和求解出来的答案与实际有一定的差距。
2、由于题中没有说明各个区域中的主要工商业，故在确定污染源的来源时不能准确定位。
3、内梅罗指数的方法过分突出极大值对污染的影响，未考虑次大值对污染的影响，未考虑危害性最大的污染因子对水质的影响。
7.3 模型的改进
对于问题一，本文采用的是传统的尼梅罗综合污染指数法, 而传统尼梅罗综合指数过分突出污染指数最大的重金属污染物对环境质量的影响和作用，在评价时可能会人为地夸大或缩小一些因子的影响作用，使其对环境质量评价的灵敏性不够高，在某些情况下，它的计算结果难以区分土壤环境污染程度的差别。
所以现在对其进行改进，改进如下：
F max +Fw '
计算公式为P mj =F ' max =，Fw 为权重值最大的污染因
2子的c i 与s ij 比值，其他符号意义同第一问。
改进的尼梅罗指数法具有如下优点：
1、以国家最新土壤环境质量标准划分污染等级。污染等级的划分随标准改变而改变，具有划时代性。
2、待测土壤尼梅罗指数的计算方法，与划分土壤污染等级时各类标准土壤尼梅罗指数的计算方法相同，所采用的基准也一样，具有同一性。 3、考虑了土壤中危害性最大的污染因子对土壤的影响。
八、模型的推广
本模型的建立对于研究减少城市污染问题和保护环境具有重要意义，尤其在当今以高能耗高污染的生产模式为主的工业时代，该模型的建立对于研究城市规划，和工厂位置选择以及交通建设时具有重要的参考价值，同时利用该模型也可以研究物质扩散现象的规律。
九、参考文献
[1] 杨启航, 数学建模, 北京:高等教育出版社,2003；
[2] 韩中庚，数学建模方法及其应用，北京高等教育出版社，2005； [3] 姜启源，谢金星，叶俊，数学模型，北京，2003；
[4] 韩中庚，数学建模方法及其应用，北京高等教育出版社，2005。
十、附录
y
00.511.522.5
x 10
4
图8 cd
y
2.5
x 10
4
0.511.52
图9 cr
y
0.5
11.522.5
x 10
4
图
10 cu
y
00.511.522.5
x 10
4
图11 hg
y
0.5
11.522.5
x 10
4
图12 Ni
y
0.5
1
1.52
2.5
x 10
4
图13 Pb
y
1
1、去除20个异常点的程序 clc
clear all
X=xlsread('data1.xls'); %导入数据
%%%%%取出单位是ug/g的六个污染元素的浓度 A1=X(:,1); A3=X(:,3); A4=X(:,4); A6=X(:,6); A7=X(:,7); A8=X(:,8);
%%%%%%处理, 找出浓度大于500ug/g的异常值并去除 k1=find(A1>500) A1(k1)=0;
k3=find(A3>500) A3(k3)=0; length(A3)
k4=find(A4>500) A4(k4)=0;
k6=find(A6>500) A6(k6)=0;
k7=find(A7>500) A7(k7)=0;
k8=find(A8>500) A8(k8)=0;
%%%%%取出单位是ng/g的六个污染元素的浓度 A2=X(:,2); A5=X(:,5);
0.51.522.5
x 10
4
图14 Zn
%%%%%%处理, 找出浓度大于4000ng/g的异常值并去除 k2=find(A2>4000)
A2(k2)=0;
k5=find(A5>4000)
A5(k5)=0;
X=[A1,A2,A3,A4,A5,A6,A7,A8]
xlswrite('chuli1.xls',X) %导出数据
2、319个选样点的三维地势图程序
clc
clear all
close all
d=xlsread ('data.xls');%导入数据
x=d (:,2);%采样点坐标x 值
y=d (:,3);%采样点坐标y 值
z=d (:,4);%采样点坐标z(海拔) 值
nx=linspace(min(x),max(x),100);
ny=linspace(min(y),max(y),100);
[xx,yy]=meshgrid (nx,ny);
zz=griddata(x,y,z,xx,yy,'v4');
surfl (xx,yy,zz);
shading interp
set (gcf,'color','w');
colormap(cool);
hold on
for i=1:44
plot3 (d(i,2),d(i,3),d(i,4),'g*')
end
hold on
for i=45:80
plot3(d(i,2),d(i,3),d (i,4),'r+')
end
hold on
for i=81:146
plot3(d(i,2),d(i,3),d(i,4),'w+');
end
hold on
for i=147:284
plot3(d(i,2),d(i,3),d(i,4),'y+');
end
hold on
for i=285:319
plot3(d(i,2),d(i,3),d(i,4),'k*');
end
title('选样点的三维地势图')
20
3、问题二中八种金属等高线图
clc
clear all
close all
A=xlsread('data2');
x=A(:,1);y=A(:,2);
[x,y]=meshgrid(0:1000:29000,100:300:19000); z3=griddata(A(:,1),A(:,2),A(:,3),x,y,'v4');
figure(1)
[c,h] =contour(x,y,z3);
clabel(c,h)
ylabel('y');
xlabel('图1 As的空间分布');
z4=griddata(A(:,1),A(:,2),A(:,4),x,y,'v4');
figure(2)
[c,h] =contour(x,y,z4);
clabel(c,h)
ylabel('y');
xlabel('图2 Cd的空间分布')
z5=griddata(A(:,1),A(:,2),A(:,5),x,y,'v4');
figure(3)
[c,h] =contour(x,y,z5);
clabel(c,h)
ylabel('y');
xlabel('图3 Cr的空间分布')
z6=griddata(A(:,1),A(:,2),A(:,6),x,y,'v4');
figure(4)
[c,h] =contour(x,y,z6);
clabel(c,h)
ylabel('y');
xlabel('图4 Cu的空间分布')
figure(5)
z7=griddata(A(:,1),A(:,2),A(:,7),x,y,'v4');
figure(5)
[c,h] =contour(x,y,z7);
clabel(c,h)
ylabel('y');
xlabel('图5 Hg的空间分布')
z8=griddata(A(:,1),A(:,2),A(:,8),x,y,'v4');
figure(6)
[c,h] =contour(x,y,z8);
clabel(c,h)
ylabel('y');
xlabel('图6 Ni的空间分布')
21
z9=griddata(A(:,1),A(:,2),A(:,9),x,y,'v4');
figure(7)
[c,h] =contour(x,y,z9);
clabel(c,h)
ylabel('y');
xlabel('图7 Pb的空间分布')
z10=griddata(A(:,1),A(:,2),A(:,10),x,y,'v4');
figure(8)
[c,h] =contour(x,y,z10);
clabel(c,h)
ylabel('y');
xlabel('图8 Zn的空间分布')
4、模型检验中海拔小于100米As 等高线图 clc
clear all
close all
A=xlsread('data2');
xx=A(:,1);yy=A(:,2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% T1=A(:,1);T2=A(:,2);T3=A(:,3);
k=find(T3
U1=T1(k);U2=T2(k);U3=T3(k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[xx,yy]=meshgrid(0:1000:29000,100:300:19000); %%%%第一种金属污染源
figure(1)
z3=griddata(T1(k),T2(k),T3(k),xx,yy,'v4');
[c,h] =contour(xx,yy,z3);
clabel(c,h)
22