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

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

摘 要

本文针对城市表层土壤受重金属污染问题进行综合分析。首先运用Matlab 软件求解8种主要重金属元素在城区空间分布，土壤受污染主要原因，传播特征以及为今后如何更好研究地质演变问题分别建立了相应的数学模型，并对其求解结果作出了分析。

针对问题1，根据各种污染物浓度在不同区域内分布的随机性，利用空间内插值法，以城市位置为平面“横纵向”，污染物浓度为“竖向”建立三维空间模型，得到这8种重金属元素的空间分布图，直观反映出污染物浓度的等值线，得到污染物浓度的分布规律和各种重金属元素的污染浓度范围区块。然后通过建立污染负荷指数法模型算出各区的)(PLI 值，生活区)(PLI 为1.8336，工业区)(PLI 为2.1573，山区)(PLI 为1.0602，交通区)(PLI 为1.9209，公园绿地区)(PLI 为1.5780；结果表明工业区的污染程度最高为2级强污染，其他区为1级中等污染。

针对问题2，根据问题1中不同功能区重金属的污染程度与附件2的数据，利用Pearson 相关性分析和主成分分析法，建立了相关性模型与主成分模型，对应找出各区重金属污染的主要因子，得到该城区的重金属污染主要元素为：u g Cd C H Pb 、、、，污染主要原因为冶炼、化物生产等工业废水、污泥，汽车尾气排放。

针对问题3，根据重金属污染物的圆型传播特征，利用物理热力学第二定理，分别建立了一维土壤迁移模型、沉降模型和最优解模型，确定了生活区有3个污染源：分别为样本点20、306、259；工业区的污染源为样本点261 ；山区的污染源为样本点62 ；交通区有2个污染源：分别为样本点245、292 ；公园绿地区的污染源为样本点315 。其结果与重金属元素在该城区的空间分布图作对比，吻合得很好。

针对问题4，通过对问题1,2,3所建模型优缺点的分析，为更好地研究城市地质环境的演变模式，还应收集城区每年生活、工业等重要污染源的垃圾排放放量、每年的生物降解量、排污企业个数的增减值、PH 值、河流所经区域等信息。根据这些信息建立时间序列模型，用此模型预测未来时间段该城区的土壤中重金属含量变化情况，结合问题1,2,3所建模型可以更好地为今后如何使用专业降浓剂保证土壤良好性的用剂量提供参考依据。

关键词：重金属 空间分布 污染负荷指数 主成份分析 最优解

一、问题重述

随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。

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

附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。

通过数学建模来完成以下任务：

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

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

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

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

二、问题分析

2.1 问题一的分析

通过对附件1的数据进行处理求出该城区的城区状貌，并对该城区利用三维插值法求出8种主要重金属元素在该城区的空间分布。并对8种金属分别进行查找出各自范围区块。

利用最高污染系数)

(CF分别求出8种重金属的各自污染系数，然后利用污染负荷指数法)

(PLI求出319个地点的重金属污染程度。根据这些地点的污染程度对5个不同区域进去综合评价得出不同区域的污染程度。

2.2 问题二的分析

研究重金属污染的主要原因既为研究重金属污染的主要来源，而同一活动可以产生许多污染物，同一种污染物的来源也不尽相同。因此利用Pearson相关分析和主成份分析方法解析重金属污染的元素和分析土壤重金属污染物的来源和类别。通过对附件2、3的数据进行分析综合运用这些统计分析方法来分析5个不同区域的8种重金属污染物的来源，并对5个不同区域查找相关性比较强的区域。最后综合5个不同区域的来源分析出该城区的重金属污染的主要原因。

2.3 问题三的分析

对于重金属的传播特征形式，利用生物学基础对各种金属在该城区的传播形式进行假设分析出三种传播特征：圆型传播、条型传播、散点传播。结合前两问求出的结果和8种主要金属元素浓度含量的空间分布图分别求出属于圆型传播、条型传播、散点传播的金属元素，由此可对这8种重金属元素综合归类为那种传播特征。 由于污染源在不同区域的位置对各个区域的影响有一个直接影响值，在污染源数目较多的场合下，要确定最优解的模型污染源是很困难的，而且在处理数据时还有考虑该地区的地势情况中的地下金属传播途径的迁移形式和大气中大气降沉形式，则综合考虑该城区的污染源就更困难。因此通过对模型的研究，建立迁移模型、沉降模型对土壤金属流动进行分析，再结合最优解模型找出该城区5个不同区域的污染源。 2.4 问题四的分析

对于前三问建立的模型和结果进行分析查找出所建立的模型对该城区和不同区域的优缺点。

由于该城区的城市布局和降水量，PH 值等不知道，因此需要去完善模型并预测地质的演变，就得收集城区每年生活工业等重要污染源的垃圾排放放量以及每年的生物降解量，降雨量，PH 值，河流所经区域等。通过这些数据建立起各重金属元素与时间有关的时间预测模型，由此可预测出未来时间段内的土壤重金属浓度，即知道该城区地质环境的演变模式。且根据不同功能区的不同安全浓度指标，以及考虑大气压降和地势海拔等因素，制定出未来时间段内如何使用专业降浓剂保证土壤良好性的用剂量，这样，既达保证了经济效益又保护了土壤环境。

三、问题假设

1.该城区城市布局合理；

2.该城区土壤重金属元素来源不受其他城区影响，数据采样都是可靠准确；

3.假设每个样本点都能很好的反映该平方公里的实际情况；

4.污染物的排放瞬时完成，且排放速率恒定；

5.气体的传播服从扩散定律，非稳态；

6.污染物的沉降速率恒定；

四、符号说明

i CF ：元素的最高污染系数； i C ：元素i 的实测含量；

i C 0：元素i 的评价标准，既背景值。 P L I

：污染负荷指数; c :土壤重金属在液相的浓度； s :土壤重金属在固相中的浓度；

ij M ：在区域点,i j 重金属浓度；

M ：总重金属的浓度；

t Q ：土壤重某金属在t 年后的含量； t ：预测年限；

n ：评价元素的个数；

z o n e

P L I 为区域污染负荷指数； m 为评价点的个数（采样点的个数）。

五、问题一的解答

5.1空间分布

对附件1的319个随机地点中的）海拔（m m y m x ),(),(做一个三维曲线拟合得该城区的城区状貌。如图1.

图1该城区状貌

由图1可以得到该城区的城区状貌，然后利用附件2的8种金属元素浓度代替海拔()m 就可以得出8种金属元素在该地区的空间分布图如图2所示。

3

4

0.5

1

1.5

2

2.5x 10

4

-15

-10-5051015202530x 10

4

-1500

-1000

-500

500

1000

1500

图2 8种主要金属元素浓度的空间分布

由图2的8种主要金属元素浓度空间分布图可对比出每种元素的浓度高涵盖点区，再结合附件1、2对数据进行处理，则可以得出该种元素的污染浓度范围区域，其具体情况如表1所示。

x 10

4

-1000

-800

-600-400-2000200400600800

x 10

4

-1000

-50005001000150020002500

x 10

4

-5000

5000

10000

15000

0.5

1

1.5

2

2.5x 10

4

-100

-50

50

100

0.5

1 1.5

2 2.5x 10

4

-1500

-1000-5000500100015002000250030003500

0.5

1 1.5

2 2.5x 10

4

100

200

300

400

500

表1 8种元素的污染浓度范围

5.2污染程度

污染负荷指数法模型

污染负荷指数法是Tomlinson 等人从事重金属污染水平的分级研究中提出来的一种评价方法。该指数又评价区域所包含的多种重金属成分共同构成，并使用了求积的统计法，通过这种方法能对整个区域各个点位各种重金属进行定量评价，并对各点的污染程度进行分级，能直观反映对环境污染最严重的元素和各种元素对环境污染的贡献程度，以及重金属在时间、空间上的变化趋势，应用比较方便]1[。能避免污染指数加和关系造成的对评价结果歪曲的现象，并能对任意给定的区域进行定量的判断。他也采用研究区土壤背景值为评价标准，因此它能很好的判断土壤的综合人为污染情况。

首先根据某一点的测量金属含量，进行最高污染系数（CF ）的计算]2[:

i

i i C C

CF 0= （1）

某一点的污染负荷指数（PLI ）为：

n n CF CF CF CF PLI 321??= （2）

某一区域的污染负荷指数（zone PLI ）为：

m m

z o n e P L I P L I P L I P L I P L I 321??= （3）

表2 污染负荷指数与污染程度之间的关系

通过Matlab 程序可以计算出各个地点的最高污染系数和污染负荷指数（程序见附录1），并结合这些点的数据综合评价得出如表3所示的五个区域的污染负荷指数（zone PLI ）。

表3 五个区域的污染负荷指数（zone PLI ）

由上表可知该城区内不同区域重金属的污染程度为中等污染的区域有生活区、山区、交通区、公园绿地区。强污染的区域有工业区。

六、问题二的解答

研究土壤重金属污染的主要原因需要调查污染土壤中各种污染物的影响因素。土壤重金属污染影响因素主要指产生重金属污染物的设备、位置、场所等污染源]3[。这些污染源产生的重金属污染有的直接进入土壤，有的则是通过大气和水体再沉降进入土壤，造成城区土壤的重金属污染。按重金属污染物产生的部门，土壤重金属污染源主要分为工业污染源、交通污染源、生活污染源、农药和肥料等。具体来源参见表4中国土壤重金属的主要来源]4[。

表4 中国土壤重金属的主要来源

同一区域可以产生许多污染物，同一种污染物的来源也不尽相同。相关分析和主成分分析等常规统计方法常用于污染物源解析中。为了更好地分析土壤重金属污染物的来源和类别，综合运用这些统计分析方法来分析5个不同区域的8种重金属污染物的来源。最后综合5个不同区域的来源分析出该城区的重金属污染的主要因素。 6.1相关分析模型

相关分析就是研究两个或两个以上变量之间相互关系的统计分析方法。在环境研究中常用来定性的分析研究各环境变量之间的相互关系，而变量之间的相互关系用相关系数和相关系数的显著水平来表示]6[。相关分析法有多种，环境研究中最常用的是Pearson 相关分析法。

相关系数矩阵 p p ij r R ?=)(

),2,1,(,1

1

p j i n x

x

r n

k kj

ki ij =-=

∑= （4） 式中：ij ji ij ij r r r r ,,1==是第i 指标与第j 指标的相关系数。

在重金属污染源解析中，当两种重金属含量绝对值最大，表明这两种重金属很可能来自同种污染源。

1).生活区的相关矩阵结果见表表5所示。

表5 生活区Pearson 相关矩阵

从相关性分析结果可以发现，土壤中Cd 与Pb 显著正相关，且相关性较强；其次As 与Ni 、Cu 也到达了显著的正相关。以此类似可作出其他4个区域的Pearson 相关矩阵（见附录2）。 6.2主成份分析模型

多变量研究中，由于变量的个数很多，并且彼此往往存在一定的相关性，因此使观察的数据反映的信息在一定程度上重叠。因子分析则是通过一种降维方法进行简化得到综合指标。综合指标之间既互不相关，又能反映原来的观察指标的信息]5[。

1.计算特征值和特征向量

计算相关系数矩阵，及对应的特征向量120p L λλλ≥≥≥，及对应的特征向量

12,,,p u u u L ，其中12(,,,)T j j j np u u u u L =，由特征向量组成m 个新的指标变量

11112121

21212222

1122n n n n m m m

nm n y u x u x u x y u x u x u x y u x u x u x L L L L L L L L L L L L L L =+++??=+++????=+++? （5）

式中1y 是第一主成分，2y 是第二主成分，L ，m y 是第m 主成分。 2.选择()p p m ≤个主成分，计算综合评价值

A.计算特征值()1,2,,j j m L λ=的信息贡献率和累积贡献率。称

1

(1,2,,)

j

j m

k

k b j m L λλ

==

=∑ （6） 为主成分j y 的信息贡献率；

1

1p

k

k p m

k

k a λλ

===

∑∑ （7）

为主成分12,,,p y y y L 的累积贡献率，当p a 接近于1(0.85,0.90,0.95)p a =）时，

则选择前p 个指标变量12,,,p y y y L 作为p 个主成份，代替原来m 个指标变量，从而可对p 个主成分进行综合分析。

B.计算综合得分：∑==p

j j j y b z 1 （8）

其中：j b 为第j 个主成分的信息贡献率，根据综合得分就可进行评价。由此可分别计算得8个不同区域的特征值和主成份因子；1).生活区的结果如表6所示。

表6 生活区解释的总方差

对生活区域的土壤样品中8中重金属进行主成份分析，提取出3个因子解释了总方差的72.797%，其中因子1解释了总方差的45.199%，因子2解释了总方差的14.165%，因子3解释了总方差的13.432%。表7给出了生活区的成份得分系数矩阵（因子载荷矩阵）

表7 生活区的成份得分系数矩阵（因子载荷矩阵）

量与某个因子的联系系数绝对值（荷载）越大，则该因子与变量关系越近。根据因子载荷矩阵，主成份1为Cd 、Hg 、Pb 的组合，这说明这几种土壤重金属污染物可能是同一来源或相似来源。主成份2为Cr 、Zn 。主成份3为As 、Cu 、Ni 。主成份1的方差累积贡献率达到45.199%，这可以解释Cd 、Hg 、Pb 为该生活区土壤重金属的主要污染元素。

通过因子分析结果参照土壤重金属的主要来源，认为该生活区土壤重金属污染物可能来源于生活化物排污，废渣，废水以及汽车尾气。Pearson相关分析和主成份分析结果表明该生活区的Cd、Hg、Pb可能都主要来自同一污染源。

同理结合Pearson相关分析和主成份分析可得出其他另外4个区域的污染源。（表见附录2），表8列出了五个不同区污染源。

表8 五个不同区污染源

源和污染元素来分析发现工业区与交通区的相关性比较大，则两者可以规划为同一污染区，污染源主要是废气、废水；生活区和山区、公园绿地区的相关性比较大，则三者可规划为同一污染区，污染源主要是化物、农药。

综合分析得出该城区的重金属污染主要元素：Cu、Cd、Hg、Pb。重金属污染物的主要原因是：冶炼、化物生产等工业废水、污泥和废气、肥料杂质，汽车防爆燃料尾气，含Cu、Hg农药。

七、问题三的解答

7.1传播特征分析

对于重金属的传播特征形式，利用生物学基础对各种金属在该城区的传播形式进行假设分析出三种传播特征：圆型传播、条型传播、散点传播。现结合求出的问题一、二结果进行分析得出传播特征情况，特别是结合8种主要金属元素浓度含量的空间分布图可以看出属于圆型传播的金属元素为：As、Cd、Cu、Hg、Pb、Zn；条型传播的金

属元素为：Cr、Ni、As；散点传播的金属元素为：As、Cd；由此可对这8种重金属元素综合归类为圆型传播特征。

7.2 模型研究

由于污染源在不同区域的位置对各个区域的影响有一个直接影响值，在污染源数目较多的场合下，要确定最优解的模型污染源是很困难的，而且在处理数据时还有考虑该地区的地势情况中的地下金属传播途径的迁移形式，而在大气中也通过大气降沉形式给土壤输送金属元素。则综合考虑该城区的污染源就更困难。因此通过对模型的研究，建立迁移模型、沉降模型对土壤金属流动进行分析，再结合最优解模型找出该城区的污染

源。

7.3 模型的建立 1.迁移模型

迁移模型描述由对流和弥散引起的土壤溶质迁移。对于重金属来说，在对流弥散过程中常伴随有较为强烈的吸附或分解过程，且存在作物根系吸收和排出。则此类土壤溶液流态可分为稳态流和非稳态流或瞬态流。由于该城区的流态并未知道，因此我们假设该城区为非稳态流。非稳态迁移模型为：

()()()t z qC C

q D S C z z t t ,,ψθθρθ-?

?????-????=??+?? （9）

式中 C:土壤重金属在液相的浓度，mg/L;θ:土壤含水率，3m /3m ;ρ:土壤干容重，kg/L ；S:土壤重金属在固相中的浓度，g g /μ；()q D ,θ：弥散系数，s m /2；q ：水流流速，m/s ；

t ：时间 ，s ；z ：土壤深度坐标，m ；()t z ,ψ：由植物根系引起的溶质吸收或排放率，s L mg ?/

()()t z R K z z t ,+?

??

???????=??φθθ （10） 式中 ()θK ：土壤水力传导率，m/s ；φ：土壤水势，m ；()t z R ,：植物根系吸水函数，s m m ?33/

2.沉降模型

沉降模型描述在风力或降水作用引起的土壤溶质累计。重金属含量的变化方程为

L ij Ar Am ij

Q I I t

M -+=??, （11）

式中 ij M ：在区域点i,j 重金属浓度，2/hm g ；t ：时间，s ；AM I ：区域回归的大气沉降重金属输入通量，()s hm g ?2/；ij Agr I ,：

由于农业活动引起的重金属输入通量，()

s hm g ?2/；L Q ：重金属渗流通量，()

s hm g ?2/。

方程（11）进一步拓展为连续变化的函数

n L m

c A g r Am M k M k I I dt

dM 1

--+= (12)

式中 M ：总重金属的浓度，2/hm g ；c k ：重金属作物吸收率系数，1/a ；L k ：重金属

渗流率系数，1/a ；m 、n ：常数。方程（12）中

z C zC M s t θρ+= （13）

()

k

b k

c b z Y k k 01ρ= （14） n

F w L zK q k 1

1????

?

?=ρ （15）

式中 ： t C ：土壤溶液的重金属浓度，kg mg /；s C ：被土壤吸附的重金属浓度，L mg /；k Y ：作物产量，()s hm kg ?2/；k b 0、k b 1：回归系数；w q ：达西水流速率，()

s m L ?2/；F K ：

Freundlich 参数，kg L /。

在应用上述方程进行计算时，为了估计方程输出的随机性，把输入参数作为随机变量处理，看作为正态或对数正态分布，给出均值和标准差，例如（12）式Am I 取值为9.021±

3.最优解模型

确定污染源最佳位置的问题，就是要在定的研究区域内寻找一个点位来确定污染源，使得污染源在这一点位上能够保证对于所研究区域内的一组污染源最优。

在一区域内找最优点，即在二维凸边形可行区域内求一点),(y x ，使得在区域内该点到给定的一组点),(i i i b a p 的最小加权距离达到最大，即

1m a x ((,))m i n (,);(,)i n

f x y Wind x y x y s ≤≤=∈ （16）

式中 ：),(y x f 为系统的影响程度，s 是所研究范围内的凸可行区域，n 为s 中现有受影

响对象的数量，),(y x d i 为位于点),(y x 上的污染源与位于),(i i i b a p 上的第i 个受影响对象的欧氏距离，即2

1

22

))()((),(i i i b y a x y x d -+-=。，i w 代表第i 个受影响的对象与污染源之间的相对不适度的大于零的加权因子。

由(16)式给出的表达式是不易求解的， 但是由于该模型的总体性质， 通过分离两项优化操作就可将方程写成下列的标准数学规化式：

（17） st n i b y a x w Z i i i 3,2,1,)()(22=-+-≤ （18） （19）

且s b a p i i i ∈),(，s 为凸边形。

这里权重是有限数， 可行空间由在正棱形范围内的锥面z 方向上形成边界 对于标准的数学规划式(17)、(18)，(19)， 应用规化论的算法，根据不同的精度要求，可以确定寻查步长，用试探法进行求解。

由以上三个模型并结合空间分布图进行分块求解，综合附件1可大致确定该城区生活区、工业区、山区、交通区、公园绿地区的污染源分别为下表所示：

表9: 生活区污染源位置

表10：工业区污染源位置

max Z (,)x y s ∈

表11：山区污染源位置

表12：交通区污染源位置

表:13：公园绿地区污染源位置

由上表可知生活区有3个污染源：分别为样本点20、306、259；工业区的污染源为样本点261 ；山区的污染源为样本点62 ；主干道路区有2个污染源：分别为样本点245、292 ；公园绿地区的污染源为样本点315 。

八、问题四的解答

8.1模型的评价

模型优点：

1. 模型是由简单到复杂一步步建立的，增强可读性、理解性，逻辑性强；

2. 在建模过程中将一些给定的数据参数化，使得模型更具一般性，可应用范围广；

3. 通过计算机模拟的方法较好地对近似计算所得结果做出了验证；

4. 问题2中的总体各种金属相关性累计量比较少，分别对5个不同区进行各自相关性分析和主成份分析，到得的效果比较好；

5. 问题3中能比较充分考虑土壤含重金属的迁移、沉降地质变化情况。

模型缺点：

本文在计算过程中只使用了对原始数据进行分析，而原始数据之间的距离比较大，不能很好的拟合出该城区的状貌和计算浓度间的关系。

问题1中的模型适合用于其它方面的等级程度划分，将模型拟合出的浓度函数关系式通过一定的变换，找到更多的数据点，近似优化函数关系，这样会得到的结果比较精确；

问题2中分别对5个不同区进行各自相关性分析和主成份分析，其中得到的交通区结果比较差，因此需要对数据再进行处理才能用。

问题3中的模型适合用于最优选址问题，也可以利用其它更为简便的参数估计方法进行求解，如粒子群算法、遗传算法等。

由于收集到的数据比较单一，没有其他与重金属传播有关的数据，因此在建立模型求解问题时准确性不高且不具有预测的功能。 8.2收集信息

为了更好的研究城市地质环境的演变模式，应收集城区每年生活工业等重要污染源的垃圾排放放量以及每年的生物降解量，排污企业的个数，降雨量，PH 值，河流所经区域等影响因子。 8.3信息分析模型 时间预测模型

由以上收集的信息，要解决分析出城区地质环境的演变模式，即是要了解未来一段时间内土壤中各种重金属污染物的含量及变化趋势，由此便于提前采取措施防止或减轻各种重金属污染所造成的负面影响。于是对城区内土壤中各种重金属污染物的主要来源进行分析，分析出影响这些重金属含量时间变化的主要影响因子，在此基础上建立各种土壤重金属含量的时间预测模型，最后调查分析预测模型中各变量的最可能的变化速率和取值范围并利用这些模型对土壤中各重金属含量进行时间预测。 根据土壤环境容量模型推导出时间预测模型]5[：

()

()K K QK K Q Q t t t --+=1/10 （20）

式中：t Q 为土壤重某金属在t 年后的含量(/)ug g ；0Q 为土壤中某重金属的起始含量(/)ug g ;Q 为每年外界重金属进入土壤的量(/)ug g ；K 为土壤重金属的年残留率；

t 为预测年限。

该模型的特点是不区分各个具体输出项(如重金属随作物带出土壤、重金属在土壤中向下层迁移等)与土壤中重金属含量的数学关系，仅用年残留率K 来计算土壤中重金属的残留与降解。

通过这些数据建立起各重金属元素与时间有关的时间序列模型，预测出未来时间段内的土壤重金属浓度，且根据不同功能区的不同安全浓度指标，以及考虑大气压降和地势海拔等因素，制定出未来时间段内如何使用专业降浓剂保证土壤良好性的用剂量，这样，既达保证了经济效益又保护了土壤环境。

参考文献

[1]. 中国环境检测总站主编, 中国土壤元素背景值[M], 北京: 中国环境科学出版社, 1990;

[2]. 徐争启, 倪师军, 应用污染负荷指数法评价攀枝花地区金沙江水系沉积物中的重金属[J], 四川环境, 2004,23(3):64-67;

[3]. 崔航. 长春市土壤重金属分布规律及土壤环境质量评估[J]. 吉林大学. 2005; [4]. 陈怀满. 中化学物质的行为与环境质量[M],京: 科技出版社, 2002;

[5]. 林艳. 基于地统计学与GIS 的土壤重金属污染评价与预测[J]. 中南大学, 2009; [6]. 吴春发. 复合污染土壤环境安全预测预警研究[D], 浙江大学, 2008;

[7]. 隋红建. 土壤重金属迁移模拟研究的现状与展望, 农业工程学报, 2006;

[8].庄世坚. 污染源最佳位置的确定. 2011.9.11.

附录：

附录1：

1.空间分布程序

clear all

clc

A=xlsread('Book2011');

x=A(:,1);

y=A(:,2);

z=A(:,12);

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

figure

[X,Y,Z]=griddata(x,y,z,linspace(0,28654)',linspace(0,18449),'v4');%插值pcolor(X,Y,Z);

shading interp%伪彩色图

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

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

2.不同区域污染程度程序

clear all

clc

C1=xlsread('Book20111生活区');

C2=xlsread('Book20111工业区');

C3=xlsread('Book20111山区区');

C4=xlsread('Book20111交通区');

C5=xlsread('Book20111公园绿地区');

A2=xlsread('Book20112')';

Iplzone_1=1;

for i=1:size(C1,1)

F_1=C1(i,:)./A2;

Ipl=(F_1(1)*F_1(2)*F_1(3)*F_1(4)...

*F_1(5)*F_1(6)*F_1(7)*F_1(8))^(1/8);

Iplzone_1=Iplzone_1*Ipl;

end

Iplzone_1=Iplzone_1^(1/size(C1,1)) %生活区污染负荷指数

Iplzone_2=1;

for i=1:size(C2,1)

F_2=C2(i,:)./A2;

Ipl=(F_2(1)*F_2(2)*F_1(3)*F_2(4)...

*F_2(5)*F_2(6)*F_1(7)*F_2(8))^(1/8);

Iplzone_2=Iplzone_2*Ipl;

end

Iplzone_2=Iplzone_2^(1/size(C2,1)) %工业区污染负荷指数

Iplzone_3=1;

for i=1:size(C3,1)

F_3=C3(i,:)./A2;

Ipl=(F_3(1)*F_3(2)*F_3(3)*F_3(4)...

*F_3(5)*F_3(6)*F_3(7)*F_3(8))^(1/8);

Iplzone_3=Iplzone_3*Ipl;

end

Iplzone_3=Iplzone_3^(1/size(C3,1)) %山区区污染负荷指数Iplzone_4=1;

for i=1:size(C4,1)

F_4=C4(i,:)./A2;

Ipl=(F_4(1)*F_4(2)*F_4(3)*F_4(4)...

*F_4(5)*F_4(6)*F_4(7)*F_4(8))^(1/8);

Iplzone_4=Iplzone_4*Ipl;

end

Iplzone_4=Iplzone_4^(1/size(C4,1)) %交通区污染负荷指数Iplzone_5=1;

for i=1:size(C5,1)

F_5=C5(i,:)./A2;

Ipl=(F_5(1)*F_5(2)*F_5(3)*F_5(4)...

*F_5(5)*F_5(6)*F_5(7)*F_5(8))^(1/8);

Iplzone_5=Iplzone_5*Ipl;

end

Iplzone_5=Iplzone_5^(1/size(C5,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)

成份 1 -0.105 0.122 0.218

0.277

0.272 0.028 0.142 0.105 2

0.505

0.145

0.001 -0.154

-0.155

0.343

0.14

0.221

山区：

山区解释的总方差

山区的成份得分系数矩阵（因子载荷矩阵） As (μg/g) Cd (ng/g) Cr (μg/g) Cu (μg/g) Hg (ng/g) Ni (μg/g) Pb (μg/g)

Zn (μg/g) 成份 1 0.022 -0.042 0.404 0.072 -0.138 0.415 -0.039 0.24 2 0.202 -0.39 0.046 0.004 -0.166 0.063 -0.376 -0.24

3

0.395

0.006

-0.01

0.489

0.425

-0.032 0.033

-0.03

交通区解释的总方差

交通区的成份得分系数矩阵（因子载荷矩阵）

As (μg/g)

Cd

(ng/g)

Cr

(μg/g)

Cu

(μg/g)

Hg

(ng/g)

Ni

(μg/g)

Pb

(μg/g)

Zn

(μg/g)

成份1 -0.076 0.246 0.066 0.217 0.146 -0.012 0.259 0.245

2 0.41 -0.107 0.284 0.051 0.134 0.378 -0.065 -0.005

公园绿地区解释的总方差

公园绿地区的成份得分系数矩阵（因子载荷矩阵）

As (μg/g)

Cd

(ng/g)

Cr

(μg/g)

Cu

(μg/g)

Hg

(ng/g)

Ni

(μg/g)

Pb

(μg/g)

Zn

(μg/g)

成份1 -0.075 0.252 0.066 0.304 -0.016 -0.023 0.307 0.3

2 -0.396 -0.084 -0.325 0.072 -0.015 -0.378 0.054 -0.007

3 0.19 -0.15

4 -0.078 -0.004 0.862 -0.11

5 0.234 -0.12