城市地表重金属污染分布及污染源解析研究

城市地表重金属污染分布及污染源解析研究

摘要

本文研究的是城市中不同功能区重金属污染分布及重金属污染来源形式和污染源位置的相关问题。

针对第一问分析该城区内不同区域的重金属污染程度的问题，我们选用图示和指标衡量两步综合反映污染程度问题。先用matlab绘出该地区重金属总量以及每种重金属元素的空间分布图，后以采集样本中8种重金属浓度和背景值通过换算，得到地质累积指数，定量的分析污染程度。我们得到的结论为：工业区普遍受到较重的重金属污染，其次为主干道路区，再次是生活区和公园绿地区，山林区基本未受到重金属的污染。

第二问要求通过数据分析重金属污染的主要原因，我们发现各种重金属不是单一产生的，多种重金属在产生的过程中往往具有伴随性。针对这种现象，我们使用因子分析法来确定各种重金属元素之间的相关性，然后利用它们的相关性定量分析各因子对8种重金属元素产生的贡献值。最终通过已知数据可以得出： Ni、Cr、Cu的污染主要是由主干道路区汽车尾气产生， Hg的污染主要是由工业“三废”产生， As的污染主要是由生活区杀虫剂等药剂产生， Zn的污染主要是由生活区垃圾焚烧和主干道路区汽车轮胎磨损产生。

第三问要求建立模型找出污染源的位置。我们发现重金属传播必定是借助于介质进行传播，则找出这种介质并研究其运行规律对解决此问有很大帮助，最终我们确定重金属元素主要以大气颗粒物和大气水为主要载体向外传播，为此研究大气颗粒行为方式成为我们的解题关键。此问中我们运用高斯烟羽模型，推导出了由于大气颗粒物沉降带来的土壤重金属浓度分布服从高斯二维方程。最终我们通过读图确定8种重金属元素的污染源个数分别为7、14、8、5、6、5、5、7个，再利用模型求出了其具体的位置坐标。

针对问题三中的模型，我们进行了稳定性的检验。首先用matlab程序处理了所求的各种重金属污染源的位置坐标，求得了其各自相关点的距离，然后根据各个坐标点的间距，算得了其各自的污染源坐标点的相关系数，再与第二问中各种金属相关系数作比较，得出其误差平均为7.1%，误差的方差为0.005，这说明我们求得的污染源位置坐标具有较好的稳定性与可信度。

第四问我们主要考虑了将地形因素信息、大气因素信，在高斯烟羽模型的基础上，对模型进行了优化，并列出了一般扩散条件下、对流边界层扩散条件下、稳定边界层扩散条件下的模型。

关键词：地质累积指数因子分析法高斯烟羽模型稳定性分析

一、问题重述

对城市土壤地质环境异常的查证，以及如何应用查证获得的数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

现对某城市城区土壤地质环境进行调查。记生活区、工业区、山区、主干道路区及公园绿地区等，分别为1类区、2类区、……、5类区。将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。

现在已知采样点的位置、海拔高度及其所属功能区，8种主要重金属元素在采样点处的浓度，8种主要重金属元素的背景值等信息。利用数据，建立模型完成以下任务:

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

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

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

(4) 分析建立模型的优缺点，为更好地研究城市地质环境的演变模式，收集更多信息建立模型。

二、问题分析

问题一首先要确定衡量重金属污染程度的指标。通过查阅资料[1]选取用地质累积指数这一指标来衡量。根据附件二、三中的数据对每类区的8种重金属污染程度做指标分析。为给出该考察地区所采样本中所含8种主要元素的空间分布图。再对附件一中的数据进行整理，计算出该城区每个样本采集点8种重金属总量，并绘图表示，简单分析每个功能区重金属总体污染情况。再作8张子图，分别绘出每种重金属元素在该城区的污染情况，详细分析该城区内不同区域重金属污染程度。

问题二要求在数据分析的基础上确定重金属污染的主要原因。重金属污染的原因众多，且重金属污染往往产生的不是一种重金属元素，某些重金属元素之间是相互伴随产生的，即产生过程并非是独立的。所以要想分析出此城区重金属污染的主要原因，需要考虑各种重金属元素之间的相关性，然后利用它们的相关性定量分析各种主要污染原因对8种重金属元素产生的贡献，以此来确定每种重金属元素污染的主要原因。

问题三要求给定污染源位置模型。首先分析重金属污染物的传播特征，通过查找资料[4]得,在现代生活工业区中，重金属主要传播特征为以大气颗粒物和大气水为主要载体向外传播，而不是主要通过土壤或其它介质传播，这确定了我们研究大气的向外运动过程即为研究重金属传播的行为。本模型中我们选择高斯烟羽模型，并由大气颗粒物在空气中的传播规律，结合高斯烟羽模型，我们推导出了由于大气颗粒物沉降带来的土壤重金属浓度分布服从高斯二维方程，由此，我们通过已知数据，列出方程，求解得出污染源位置。

问题四要求分析模型优缺点，并收集所需更多信息条件下建立新模型。此问实质是

模型优化问题。对于我们之前建立的模型三是在一种理想状态下求解得的模型，而实际情况要远远复杂的多。因为重金属传播要不仅受到大气——即载体本身性质的影响，还会受到地形、地貌等一系列空间地理环境的影响。因此，要建立新的模型，我们需要收集与大气参数、地形参数有关的信息，在原有的模型上对有关参数进行修改。

三、 模型假设

1. 视大气为理想气体，遵守理想气体状态方程；

2. 污染源产生污染物的浓度是连续且均匀的；

3. 假设所有污染源离地面的高度为零；

4. 污染源是固定不动的，即其位置不会随时间发生改变；

5. 重金属由空气中沉降时为干沉降，即是由于重力作用或附着的颗粒相互碰撞而引起

的，不考虑风力等其他因素的作用； 6. 在水平方向，大气扩散系数呈各向同性。

四、 符号说明

ij I ：第i 类区第j 种重金属元素的地质累积指数；

ij C ：第i 类区所有第j 种重金属元素的样本中第j 种重金属的平均浓度，单位为

（/g g μ）；

j E ：j 种重金属元素的背景值，单位为（/g g μ）

； (,,)X x y z ：点（x,y,z ）处泄漏气体浓度，单位为kg/m 3

；

Q ：源强（即泄漏速度）

，单位为kg/s ； H ：污染源的有效高度，即污染源离地面的高度。

五、 模型建立

(一)

模型一

(1) 重金属元素在该城区的空间分布

为描述重金属在该城区的空间分布，我们首先利用题目附件一中的数据，用matlab 软件绘制出该地区的地貌及各个功能区分布图（如图1）（源程序见附录一）。

图1 城区地貌、区域分布模拟图（颜色表示不同功能区）

a)重金属元素总量的分布图

对题目附件二中的数据进行处理，统计出该城区每个采集样本点8种主要重金属元素总的含量，用matlab读取处理后数据，得到该城区各采集点中样本的重金属元素的

累积含量分布图（如图2）。

图2 重金属总量污染空间分布图

b)8种重金属元素的分布子图

根据附件一、二中的数据可分别整理出每种重金属元素在该城区的所有样本采集点

所对应的浓度值。并依据整理出的数据用matlab绘制出该地区每种重金属的空间浓度

分布图如下：

图3 As的空间浓度等值分布图

图4 Cd的空间浓度等值分布图

图5 Cr的空间浓度等值分布图

图6 Cu的空间浓度等值分布图

图7 Hg的空间浓度等值分布图

图8 Ni的空间浓度等值分布图

图9 Pb 的空间浓度等值分布图

图10 Zn 的空间浓度等值分布图

(2) 污染程度衡量指标及污染程度分析

为结合图1~10，并能定量分析出该城区不同区域重金属的污染程度，我们必须先确立一个合理的衡量污染程度的指标，通过查阅资料[1]并比较，我们选取地质累积指数这一指标来衡量。其定义式如下所示

2log (

)

ij ij j

C I kE =

——（1）

其中ij I 为第i 类区第j 种重金属元素的地质累积指数；ij C 为第i 类区所有第j 种重金属元素的样本中第j 种重金属的平均浓度，单位为（/g g μ）；k 为修正指数， 1.5k =，是考虑到由于成岩作用可能会引起背景值的变动；j E 第j 种重金属元素的背景值，单位为（/g g μ）。

地质累积指数分级标准与污染程度之间的相互关系见下表

表1：地质累积指数的划分级别

ij I

级别 污染程度 <0 1 无污染 0~1 2 轻度污染 1~2 3 中度污染 2~3 4 中度污染到强污染

3~4 5 强污染

4~5 6 强污染到极度污染

5~6

7

极度污染

(二)

模型二

通过查阅资料[2][3],我们采用多元统计数学方法之一的因子分析法，它根据多个实测变量（重金属元素）之间的相关性，运用数学变换，将多个变量转换为少数几个线性不相关的综合指标（不相关的几个主要原因），从而简化数据处理，其目的在于对大量观测数据，用较少的有代表性的因子（目标）来说明众多变量所提取的主要信息，提示出多个变量的因果关系。具体方法如下：

因子分析从变量的相关矩阵出发将一个m 维随机向量X （数据总库）分解成低于m 个且有代表性的公因子和一个特殊的m 维向量，是公因子数取得最佳的个数，从而使对m 维随机向量的研究转换成较小个数的公因子研究。

设一共有n （n=319）个样本，n 个指标构成样本空间X ： ()ij X x = n m ? 1,2...i n =；1,2...j m =。 之后我们对数据做如下处理步骤：

(1)原始数据的标准化，标准化公式为：

,()/ij ij j j X X X δ=-

其中ij X 为第i 个样本第j 个指标值，而j X 和j δ分别为j 指标的均值和标准差。标准化目的在于消除不同变量量纲的影响，而且标准化转化不会转变变量的相关系数。 （2）KMO 检验，判断其是否符合相关性分析条件。KMO 值定义为

222

ij i j

ij ij

i j

i j

r KMO r p ≠≠≠=+∑∑∑∑∑∑ 其判断标准是：

表2：KMO 检验法判断标准

0.9

若符合上表中的标准，则计算标准化数据的相关系数阵，相关性公式如下：

1

2

2

1

1

()()

()

()

N i i i XY N N i i i i X X Y Y r X X Y Y ===--=

--∑∑

∑

其中X,Y 表示不同的重金属元素。

（3）进行正交变换，使用方差最大法，其目的是使因子载荷两极分化，且旋转后因子仍正交。

（4）确定因子个数，计算因子得分，进行统计分析。 (三)

模型三

通过模型二结论和查找资料[4]可知，重金属污染物的传播特征主要是以大气颗粒与大气水为主要载体，后经过沉降，进入土壤等环境并以不同化学形态进入城市环境生态系统，其具有不可降解性，可长期存在于环境中。

依据假设1可知单点污染源所释放出的重金属在大气中水平方向、垂直方向的扩散，其浓度分布都可看做高斯分布，即重金属随着大气的扩散符合高斯烟羽模型[5]。则下风向任意一点X （x,y,z ）处气体浓度的函数为

22

(,,)()ay bz X x y z A x e

e

--= ——（2）

由概率统计理论可写出方差表达式为

220

2

20

y z y Xdy

Xdy

z Xdz

Xdz

δδ∞

∞

∞

∞

??=?????=??

??

?

?

——（3）

在假设2的条件下，可写出源强（污染物流出的流量）的积分公式

Q uXdydz +∞

+∞

-∞

-∞

=?

?

——（4）

式中y δ、z δ为泄漏气体在y 、z 方向分布的标准差，单位为m ；X （x,y,z ）为任意一处泄漏气体的浓度，单位为kg/m 3；u 为平均风速，单位为m/s ；Q 为源强（即泄漏速度），单位为kg/s 。

将（2）代入（3）式积分可得

221212y z a b δδ?=??

?

?=??

——（5） 将（2）、（5）代入（4）式中，积分可得

()2y z

Q

A x u πδδ=

——（6）

再将（5）、（6）式代入（2）式可得

22

22(,,)exp[()]222y z y z

Q y z X x y x u πδδδδ=-+

上式为无界空间，连续点源扩散的高斯模型公式，然而在实际中，由于地面的存在，

烟羽的扩散是有界的。可以把地面看做一面镜子，对泄漏气体起全反射作用。可把任意一点处的浓度看做两部分贡献之和：一部分是不存在地面时所造成的泄漏物浓度；另一部分是由于地面反射作用增加的泄漏物浓度。该处的泄漏物浓度即相当于不存在地面时由于（0,0，H ）的实源和位于（0,0，-H ）的像源在该点处所造成泄漏物浓度之和。

其中实源贡献为：

22

22

11()(,,)exp()exp()222a y z y z

Q y z H X x y z u πδδδδ-=-- 其中像源贡献为：

22

2

2

11()(,,)exp()exp()222b y z y z Q y z H X x y z u πδδδδ+=-- 则该处实际浓度为：

222222

(,,)(,,)(,,)

11()1()exp()[exp()exp()]2222a b y z y z z X x y z X x y z X x y z Q y z H z H u πδδδδδ=++-=-?-+- ——（7）

根据假设3，本题中H=0，则空气中任意一点的重金属污染浓度为

22

2

2

11(,,)exp()exp()22y z y z Q y z X x y z u πδδδδ=-- 其仍然服从高斯烟羽模型的形式，且根据假设4只要污染源的位置没有发生改变，

在污染源所在的空间四周，其大气中的重金属浓度分布就会一直呈高斯烟羽模型中的分布状态。而根据假设5，我们知道该空间内的重金属最终只会沉降在其下方的地表，由此，我们推断其沉降后的地表重金属的浓度分布也定会与其在空中的分布方式是一样的，所以，我们得出空气中的重金属沉降到土壤表层时，土壤表层中的重金属含量服从二维平面的高斯分布，即

22222

2

()()()()1

1

(,)exp{[2]}2(1)21x y y x xy xy x x y y

x y xy x y y x c x y μμμμρρσσσσπσσρ-----=

-+-- 由假设6可知，重金属在扩散为各项同性，所以上式中相关性系数0xy ρ=，则上式可写为：

2

222

()

()1

1(,)exp{[]}22y x x y x y

y x c x y μμπσσσσ--=-+ ——（8） 当某种重金属的污染源有多个时，则任意一点的重金属浓度即为每个污染源在该点

泄漏浓度的叠加。通过模型一中的重金属分布图，我们可以确定区域内每一种重金属污染源的个数你，然后我们将附录中给出的319个点的各种重金属浓度代入由n 个污染源叠加的方程组中从而求解出n 个污染源的空间坐标。 (四)

模型四

针对模型评价中模型三的缺点，为了更好地研究城市地质环境的演变模式，通过查资料[6]得，还需要收集以下信息来完善我们的模型：

(1) 气象参数：平均风速、水平向及垂直向湍流量脉动、温度梯度、位温、水平拉格

朗日时间尺度以及行星边界层参数等气象数据；

(2) 特定的重金属元素污染的特点，如存在形态、扩散的速度等；

(3) 大气扩散特征参数：下垫面类型、指数衰减、平均时间、地形高度等。 考虑到地形高度等地形条件信息，我们可以对模型三中（7）式进行优化得

00

(,,)(,,)(1)(,,)(,,)(,)(,)0.5(1)

(,,)(,,)T a H

a i

x y z f x y z f x y z Q

x y z p y x p z x u

f x y z dz x y z dz

z z z ρρρρρρ∞=?+-?=

??=?+ΦΦ==-??

式中(,,)T x y z ρ表示总浓度；(,,)a x y z ρ表示沿地形抬升的烟羽浓度；Φ表示烟羽质量与总烟羽质量的比值；Q 表示源的泄露速率；u 表示有效风速值；(,)p y x ，(,)p z x 分别表示水平方向、垂直方向浓度分布的概率密度函数；f 表示权函数；a z 表示有效高度；i z 表示该点地形的高度值。

在模型中加入气象参数信息，在模型三中公式（7）基础上可建立对流边界层扩散方程，如下：

2

22102

2

2

2(,,)(,,)(,,)(,,)

(,,)exp[]222(2)(2){exp[]exp[]}

22T d r p j d j m y y j

j j j j j

j

x y z x y z x y z x y z Q y x y z u z h mz z h mz ρρρρλρπσσσσ

σ

∞===++=?-??

--++-

+-

∑∑

式中(,,)T x y z ρ表示烟羽的总浓度；(,,)d x y z ρ表示污染源直接排放浓度；(,,)r x y z ρ表示虚拟源的排放浓度，其计算公式与(,,)T x y z ρ相似；(,,)pr x y z ρ表示夹卷源排放浓度，其

计算公式为简单的高斯扩散公式；j λ表示高斯分布的权系数，1λ为上升气流，2λ为下降气流；j h 为有效源高；j σ表示垂直扩散系数。

在模型中加大气扩散参数特征信息，在模型三中公式（7）基础上可建立稳定边界层扩散方程，如下：

2222

2

2(,,)(2)(2)1

{exp[]exp[]}2221

exp[]

22T z y p z p z z n z z z y y y

Q

x y z F F u

z h nh z h nh F y F ρσσπσσπσ∞=-∞

=

??-+++=?-+-=?-∑ 式中(,,)T x y z ρ表示烟羽的总浓度；z F 表示烟羽的稀释，使用边界层有限参数进行计算；

y F 表示烟羽的散布，使用边界层有限参数进行计算；p h 表示烟羽的高度；z h 表示垂直

混合层的极限高度；y σ、z σ表示烟羽在水平方向、垂直方向上的扩散参数。

六、 模型求解与分析

（一） 模型一求解与分析

通过对题目附件二、三中的数据处理得到第i 类区所有第j 种重金属元素的样本中

第j 种重金属的平均浓度ij C ，如下表所示

表3：各类区重金属的平均浓度

1类区 2类区 3类区 4类区 5类区 背景值Ej As （μg/g ） 6.27 7.25 4.04 5.71 6.26 1.8~5.4 Cd （ng/g ） 289.96 393.11 152.32 360.01 280.54 70~190 Cr （μg/g ） 69.02 53.41 38.96 58.05 43.64 13~49 Cu （μg/g ） 49.4 127.54 17.32 62.21 30.19 6.0~20.4 Hg （ng/g ） 93.04 642.36 40.96 446.82 114.99 19~51 Ni （μg/g ） 18.34 19.81 15.45 17.62 15.29 4.7~19.9 Pb （μg/g ） 69.11 93.04 36.56 63.53 60.71 19~43 Zn （μg/g ） 237.01 277.93 73.29 242.85 154.24 41~97 将上表中的有关数据依次代入公式（1）得到下表，其中背景值取题目附件中所给的平均值。

表4：每种重金属在五类区的平均地质累积指数 1类区 2类区 3类区 4类区 5类区

As 0.2155 0.4250 -0.4186 0.0805 0.2132 Cd 0.5724 1.0115 -0.3564 0.8846 0.5247 Cr

0.1715 0.1999 -0.2552 0.3201

-0.0916

Cu 1.3190 2.6874 -0.1931 1.6516 0.6086

Hg 0.8255 3.6130 -0.3581 3.0893 1.1311

Ni -0.0086272 0.1026 -0.2560 -0.0664 -0.2710

Pb 0.5717 1.0006 -0.3470 0.4502 0.3847

Zn 1.1953 1.4251 -0.4979 1.2304 0.5755

根据上表中的地质累积指数数据与表1的级别划分，再结合图2-图10具体分析该城区不同区域的重金属污染情况。

从图2中的样本采集点所在位置的重金属浓度空间分布看出，绝大多数山林区的重金属浓度要远远小于其它四类区重金属浓度。再结合表3中山林区（即3类区）的地质累积指标值可以看出，其8种重金属衡量指标值均小于0，根据表1中的划分类别属于无污染，即从总体情况来看，山林地区未受重金属污染。

图3为As的空间分布图，从图中的浓度等高线来看，该城区的绝大多数区域As的浓度较低，这一点从表3中各类区的As元素的平均地质累积指标大小也可看出，除了山林地区是无污染，其余也仅为轻度污染。但图2显示该城区任有几处区域As的含量较高，且主要是在工业区和主干道路区附近，说明这些地点重金属As污染较为严重。

图4为Cd的空间分布图，大多数区域的含量为中间水平，结合表3即为轻度污染水平，但在主干道路区和某些工业区附近Cd的含量远高出平均水平。

图5为Cr的空间分布图，结合表3可知大多数地区的含量较低，污染程度为无污染或轻度水平。但在某一个主干道路区附近其样本中重金属含量高达500多μg/g，则该地点的重金属Cr污染较为严重。

图6为Cu的空间分布图，结合表3可知除了在山林地区无污染，其它类区都有不同程度的污染，公园绿地有着轻度的Cu污染，其它类区均为中度污染。而图5显示在工业区与主干道路区内存在着Cu重金属含量异常高的地点，这表明该地点重金属Cu污染很严重。

图7为Hg的空间分布图，结合表3的平均地质累积指数可看出除了在山林区外重金属Hg在其它类区均有程度不同的污染，其在生活区为轻度污染，在公园绿地为中度污染，在工业区和主干道路区污染级别达到了强污染，在图6中同样可以看到此情况。

图8为Ni的空间分布图，结合表3可以看出除了在工业区有轻度的污染外，其余均为无污染。

图9为Pb的空间分布图，结合表3可知除了在山林区外，各类区均有不同程度的污染，从平均水平上来看，工业区为中度污染，其余的为轻度污染。

图10为Zn的分布图，结合表3可知除了在山林区外，重金属Zn在其余类区均有不同程度的污染，从平均水平上说，在公园绿地区为轻度污染，其余均为中度污染。而图9中也显示出在工业区的附近重金属Zn的含量要明显的高于其它区域，所以某些工业区附近Zn污染较为严重。

综上所述，工业区的重金属污染最为严重，已有一项重金属指标达到了强污染水平，一项达到中度与强污染水平之间，三项达到中度污染。其次为主干道路区，一项指标已经达到强度污染水平，其余在轻度污染或中度污染。就污染程度而言，再次是生活区和公园绿地区。而从总体情况来看，山林区基本未受到重金属的污染。

（二）模型二求解与分析

依据题中给出的数据，对此城区土壤单点样重金属元素含量的数据标准化处理。然后分步得出结论：

首先，算出KMO=0.778>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与其它元素的相关性很弱，剩下的元素之间的相关性基本上都一般。从成因来分析，相关性较好的元素可能在成因和来源上有一定的关联。

表5：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统计软件中计算出相应的因子的特征值和累值贡献率，得到表6。

表6：旋转前后因子的特征值和累值贡献率

因子

旋转前旋转后

总的特征

值

占总百分

量的百分

率

累计贡献

率

总的特征

值

占总百分

量的百分

率

累计贡献

率

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.75

6 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之间）为重要原因。

最后，因子分析的主要目的是将具有相近的因子载荷的各个变量置于一个公因子之下，正交方差最多旋转使每一个主因子只与最少个数的变量有关系，而使足够多的因子负荷均很小，以便对因子的意义作出合理的解释，输出见表7。

表7：各元素与各因子的关联度

元素因子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），查找资料[7]，可发现汽车燃烧的石化燃料所产生的尾气中含有Ni、Cr和Cu，那么因子1极有可能是汽车燃烧的石化燃料所产生的尾气。

元素Cd和Pb对于因子2的关联度分别达到0.877和0.832，是此因子产生的主要污染物，从相关性表里也可以看出Cd和Pb的相关性好（相关性系数达到0.660），查找资料[7]，可发现工业中电镀、采矿、冶炼、燃料、电池和化学工业中排放的废气、废水和废渣以及汽车排放的尾气中都含有比较多的Cd和Pb，由于因子1是汽车燃烧的石化燃料所产生的尾气，所以可以判断因子2主要是工业中电镀、采矿、冶炼、燃料、电池和化学工业中排放的废气、废水和废渣。

元素Hg对于因子3关联度达到0.953，是此因子产生的主要污染物，查找资料[7]，可发现工业中氯碱、塑料、电子等工业排放的汞废水中含有大量的Hg，所以可以判断因子3为工业中氯碱、塑料、电子等工业排放的汞废水。

元素As对于因子4的关联度达到0.970，查找资料[7]，可发现生活中使用的杀虫剂、杀菌剂、灭鼠剂和除草剂等用品产生的废气烟尘中含有As，而且其它元素对于此因子的关联度偏低，也正好验证了因子4很有可能是生活中使用的杀虫剂、杀菌剂、灭鼠剂和除草剂等用品产生的废气烟尘。

元素Zn的对于因子5的关联度达到0.904，其它的都比较低，查找资料[7]可发现汽车轮胎磨损产物以及生活中燃烧煤制品和垃圾焚烧产生的废气、烟尘，所以因子5为汽车轮胎磨损产物以及生活中燃烧煤制品和垃圾焚烧产生的废气、烟尘。

综上所述，通过分析数据可以得出此城区重金属污染的主要原因有以下五点：

a)汽车燃烧的石化燃料所产生的尾气，主要引起重金属Ni、Cr、Cu的污染；

b)工业中电镀、采矿、冶炼、燃料、电池和化学工业中排放的废气、废水和废渣，主

要引起重金属Cd、Pb的污染；

c)工业中氯碱、塑料、电子等工业排放的汞废水，主要引起重金属Hg的污染；

d)生活中使用的杀虫剂、杀菌剂、灭鼠剂和除草剂等用品产生的废气烟尘,主要引起重

金属As的污染；

e)汽车轮胎磨损产物以及生活中燃烧煤制品和垃圾焚烧产生的废气、烟尘，主要引起

重金属Zn的污染。

（三）模型三求解与分析

首先，确定每种重金属元素污染源的个数。观察图3至图10，可以确定重金属元素As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的污染源个数分别为7、14、8、5、6、5、5、7。

附件一、二已给出样点的坐标与相对应的8种重金属的浓度，公式（8）也以给出单污染源时据污染源一定距离的点处的重金属浓度关系式，用matlab读取附件一二中的数据，并设定其服从二维高斯分布，便可得到8种重金属元素每个污染源的坐标，如下表所示

表8：各重金属元素污染源坐标

重金属元素污染源坐标

As (18230,10060)、（12740,3168）、（27500,12110）、（4342，7640）、（6946，7268）、（1737,2795）、（18520,2050）

Cd （21460,11540）、（17660，3913）、（10130，1677）、（12450,2050）、（4342,2423）、（8104,2982）、（3184,5777）、（2026,2982）、（4683，5107）、（5812,8526）、（4342,10620）、（8904，16610）、（13370,4601）、（8237,4756）

Cr （3184,5777）、（4342，5845）、（10460,5353）、（3495,4383）、（3473,5777）、（10710,5591）、（4052,5218）、（3763,5404）

Cu （2315，3541）、（3184，5963）、（11000，16400）、（3473，5777）、（10710，8200）

Hg （2605，2609）、（2315，3541）、（15340,9193）、（13600，2423）、（8683，12110）、（6946，7454）

Ni （3184，5777）、（2315，3541）、（22290，12110）、（27500，12110）、（24310，12300）

Pb （4631，5032）、（2026，3354）、（3473，10250）、（4342，2423）、（14180，11930）

Zn （13600，9690）、（9551，4472）、（12740，2982）、（4920，7268）、（3184，5777）、（2315，3541）、（4342，2423）

七、模型稳定性分析

为了检验模型三的稳定性，我们采用matlab程序处理了所求的的各种重金属污染源的空间坐标，求得了其各自相关点的距离（源程序见附录二）。

然后根据各个坐标点的间距，我们算得了其各自的污染源坐标点的相关系数，入下表所示。

表9：各种重金属污染源坐标的相关系数

相关

As Cd Cr Cu Hg Ni Pb Zn 系数

As 1.000 0.235 0.161 0.200 0.429 0.314 0.229 0.367

Cd 0.235 1.000 0.446 0.357 0.250 0.286 0.371 0.347

Cr 0.161 0.446 1.000 0.468 0.229 0.775 0.800 0.804

Cu 0.200 0.357 0.468 1.000 0.400 0.880 0.440 0.600

Hg 0.429 0.250 0.229 0.400 1.000 0.267 0.500 0.450

Ni 0.314 0.286 0.775 0.880 0.267 1.000 0.440 0.400

Pb 0.229 0.371 0.800 0.440 0.500 0.440 1.000 0.429

Zn 0.367 0.347 0.804 0.600 0.450 0.400 0.429 1.000 然后我们在将上表中的数据与之前通过spss13.0处理所得的各种金属相关系数（见上文表5）作比较，得出我们用模型三求解得到的污染源坐标相关系数的误差，如表10所示。

表10：污染源坐标相关系数误差

相对

As Cd Cr Cu Hg Ni Pb Zn 误差

As 0 0.08 -0.02 -0.04 0.025 0.009 0.021 0.005

Cd 0.08 0 -0.270 0.101 0.057 0.013 0.043 0.019

Cr -0.02 -0.270 0 -0.120 0.630 0.082 0.089 0.096

Cu -0.04 0.101 -0.120 0 -0.040 0.077 -0.150 0.055

Hg 0.025 0.057 0.630 -0.040 0 -0.030 0.067 0.290

Ni 0.009 0.013 0.082 0.077 -0.030 0 -0.033 0.008

Pb 0.021 0.043 0.089 -0.150 0.067 0.033 0 -0.130

Zn 0.005 0.019 0.096 0.055 0.290 0.008 -0.130 0

上表中除了Hg与Cr污染源相关性系数的误差达到63%外，其余误差平均值为7.1%，再求出各个误差的方差为0.005。误差基本维持在7.1%左右，这表明模型求得的数据相关性与用spss13.0对原始数据处理所得的相关系数具有很好的一致性。而误差的方差为0.005很小，表明了我们的模型三的波动性比较小，稳定性强，可信度高。

八、模型推广

由于人类活动的日益频繁，其对自然界的影响也日益严重，含重金属的污染物通过各种途径进入土壤, 造成土壤中相应重金属元素的富集。食物链中重金属来源主要是植物通过被动从土壤中吸收, 对人类社会健康可持续发展严重危害。因此, 查明土壤中重金属污染物质来源, 并从源头上加以控制, 对实施污染治理具有十分重要的意义。

本论文中我们建立了四个模型，分别对重金属的污染的浓度分布、污染来源方式和污染源位置进行了探究。对于污染浓度分布模型，我们可以将其推广到地球地质物理学和环境科学的研究当中去，其对于我们研究不同地质的不同土壤分布有一定的借鉴性。对于模型二，其运用的因子分析法在统计学中寻找各个因素的内在联系时，具有很广的

参考面。而模型三、四，其可以推广到矿物勘探方面，具有很大的应用潜力，通过此模型，对模型上午更深一步处理，可推广到物探中的矿源探究上。

九、模型评价

优点：

1.通过matlab对大量样本数据进行处理，使其由319组数据变成形象的重金属浓度等

势图，这对于问题的说明很形象且具体；

2.模型二中应用因子分析法，它根据多个实测变量之间的相互关系，运用数学变换，

将多个变量转变为几个线性不相关的综合指标，从而简化数据处理；

3.模型三结合已经画出的浓度等势图先确定每一个重金属的污染源个数，然后再对已

知的多个污染源坐标进行求解，大大简化了模型求解的步骤，降低了模型求解过程中的难度。

缺点：

1.风可推动空气的高速流动，从而有助于大气中的重金属向外扩散，且会造成重金属

在大气中扩散各向异性；

2.模型三中在确立污染源时，只能确定点状污染源，而无法确定线状或带状污染源的

位置；

3.在考虑大气重金属污染颗粒分布特性时，没有将地形等因素考虑到，这使得模型在

地势不平的山区和有高达建筑物的市区时，有一定的偏离。

4.模型三中研究重金属传播特性时，我们没有考虑8种重金属每种的自身传播特点，

而是运用了同一种模型对其进行研究，有失偏颇。

十、参考文献

[1]张江华、赵阿宁、王仲复、柯海玲，陈清华，内梅罗指数和地质累积指数在土壤重

金属评价中的差异探讨，黄金，第8期/第31卷：43-46,2010年；

[2]王雄军、赖健清、鲁艳红、李德胜、周继华、王建武，基于因子分析法研究太原市

土壤重金属污染的主要来源，生态环境，第17卷第2期：671-676,2008年；

[3]史贵涛、陈振楼、许世远、王利、张菊、李海雯、李丽娜，上海城市公园土壤及灰

尘中重金属污染特征，环境科学，第28卷第2期：238-242，2007年；

[4]孔祥宇，成都理工大学，成都市中心地区大气重金属污染特征分析与区域扩散预测，

2010年硕士学位论文；

[5]高斯烟羽模型，https://www.360docs.net/doc/8f16976782.html,/view/632a23c56137ee06eff9182d.html，

2011/9/11；

[6]翟绍岩、赵敏、徐永清、刘明花，AERMOD模型原理及应用，https://www.360docs.net/doc/8f16976782.html,，

2011/9/10；

[7]百度百科“砷”https://www.360docs.net/doc/8f16976782.html,/view/34702.htm#14，

百度百科“镉污染”https://www.360docs.net/doc/8f16976782.html,/view/58884.htm#3，

百度百科“铬污染”https://www.360docs.net/doc/8f16976782.html,/view/1956687.htm#2，

百度百科“铜污染”https://www.360docs.net/doc/8f16976782.html,/view/584695.htm，

百度百科“汞污染”https://www.360docs.net/doc/8f16976782.html,/view/862254.htm#2，

百度百科“镍污染”https://www.360docs.net/doc/8f16976782.html,/view/1874982.htm，

百度百科“铅污染”https://www.360docs.net/doc/8f16976782.html,/view/616483.htm，

百度百科“锌污染”https://www.360docs.net/doc/8f16976782.html,/view/1694738.htm，

2011/9/10；

十一、附录

附录一城市地貌图和8种重金属空间分布图源程序

clc

clear

close all

D=load('总数据.dat');

M=load('JinShu.dat');

x=D(:,2);

y=D(:,3);

z=D(:,13);

c=D(:,5);

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

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

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

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

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

marker={'*','o','s','^','p'};

color={'k','r','y','c','g'};

mat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn','所有金属'};

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

for j=1:9

% 等高线图形

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]);

title([' ',mat{j}, ' '])

xlabel('X（m）')

ylabel('Y（m）')

colormap jet

colorbar

grid on

hold on

for i=1:5

loc=c==i;

plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',c olor{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

end

% 三维体图

figure

h=surf(xi,yi,zi);

set(h,'cdata',ci);

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')

附录二：求个污染源距离的matlab程序

clc;clear

ShapeX =[4631 2026 3473 4342 14180

13600 9551 12740 4920 3184 23154342];

ShapeY =[5032 3354 10250 2423 11930