城市市区表层土壤重金属污染的分析模型 - 图文

城市市区表层土壤重金属污染的分析模型

摘要：

本文围绕某城市区域表层土壤重金属污染问题展开讨论，对该城区内土壤重金属的空间分布污染程度、污染的主要原因、污染源的确定分别建立了模型，并对其结果做出了分析。

问题1，依据附件1和附件2所给的采样点信息，应用GIS软件绘制出8种主要重金属元素在该区域的空间分布，并利用单项污染指数和综合污染指数模型，分析得到该城区不同区域内重金属的污染程度：在生活区,Cr、Zn为重度污染，Cd、Cu、Pb为中度污染，Hg为轻度污染；在工业区，Cd、Cu、Hg、Zn为重度污染，Pb为中度污染，As、Cr、Ni为轻度污染；在山区，Cd 、Cu、Cr、Ni、Pb为轻度污染；在主干道路区，Cd、Cu、Hg、Cr 、Zn为重度污染，As、Ni为中度污染，Pb为轻度污染；在公园绿地区，Zn为重度污染，Cd、Hg为中度污染，Cu、Pb为轻度污染。

问题2，应用数理统计分析方法，对8种土壤中重金属元素之间的相关性和主成分进行分析，利用MATLAB编程的计算孙结果可知：Cr和Ni之间关系最为密切，含量受彼此影响较大，来源大致相同。Pb和Cd、Cr和Cu及Cu和Pb相关性水平较高，说明来源极为相近。重金属As、Cd、Cr、Cu、Ni、Pb、Zn在主成分1上类聚，重金属Hg在主成分2上类聚。主成分1、2不是仅仅来源于自然界，除受成土母质作用外，人为因素的输入也对其产生了明显的影响。

问题3，在土壤重金属含量空间结构特征的理论基础上，根据重金属污染物的传播特征，建立半方差函数理论模型。首先，采用域法识别特异值，对原始数据进行处理，使其符合正态分布；其次，借助GS+9.0软件系统，完成半方差函数理论模型间的自动拟合、最优选择、参数计算及半方差图，统计得到各项参数；然后，对重金属含量、重金属空间结构特征进行分析，绘制出土壤中8种主要重金属含量的Krining插值图；最后，将插值图叠加，可以得出污染最严重的区域的范围是：以采样点（2883,3617,15）、（4043,1895,14）、（3520,4357,7）、（5062,4339,5）围成的区域（即为污染源的位置）。

问题4，在问题3的基础上，对半方差函数理论模型进行分析评价。为了更好地研究城市地质环境的演变模式，分析重金属污染物的传播特征，现考虑采样时间对元素污染程度的影响，应在不同时间段进行样品采样，建立数理模型，得到改进后的重金属污染程度与时间的函数关系。

最后，对文中所建模型进行合理评价，并提出改进方向，Krining是一种优化的插值方法，模型改进后可以广泛的应用到地质、土壤、水文等领域。

关键词: Krining插值 主成分分析 GIS GS

1

1.问题重述

1.1问题的实际意义

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

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

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

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

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

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

(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？

2.问题分析

2.1关于问题1的分析

该问题要求给出8种主要金属元素在该城区的空间分布，并分析该城区不同区域重金属的污染程度。题目中给出了该城市5个区域219个采样点位置、海拔高度及其8种主要元素的浓度等信息。要得到各金属元素在该城区的空间分布，需要把附件1采样点的信息和附件2中8种重金属元素在采样点的浓度在图中反映出来。对此在GIS环境下利用空间数据插值方法研究8种重金属元素的空间分布特征，得到城区重金属元素的污染程度。

2

为了进一步分析该城区不同区域重金属元素的污染程度。采用土壤重金属污染评价方法，依照国家环境质量一级、二级、三级标准值，利用单项污染指数和综合污染指数模型，在EXCEL中对附件1和附件2所给数据进行处理，得出不同区域不同元素的污染等级。 2.2关于问题2的分析

问题2要求通过数据分析，说明重金属污染的主要原因。文中应用统计分析方法分析环境中的污染原因，主要是对8种重金属元素进行相关性及主成分分析，从整体角度揭示该城市表层土壤重金属污染的主要原因。

在分析过程中，以主成分分析法为主要的研究方法研究该城区土壤中重金属的污染原因。首先对采样数据进行标准化处理，即通过数学处理使各变量的均值为0，均方差为1，以消除由于量纲和单位所造成的差异，保证所得的数据能够真正体现原始的特性；然后将标准化的数据经过MATLAB中编程及SPSS软件，可以得到变量的特征值及因子载荷矩阵、因子得分；最后综合土壤中重金属之间的相关性和主成分因子载荷得到重金属污染的主要原因。 2.3关于问题3的分析

该问题要求分析重金属污染物的传播特征，建立模型，确定污染源的位置。 分析重金属的传播特征，可以对其空间结构特征进行研究。重金属含量的空间结构特征就是指土壤中重金属含量的空间相关性及其变化规律。空间结构特征分析同样是通过确定半方差函数最佳拟合模型，并利用最佳拟合模型及其参数（块金值C0、基台值C0?C、块金效应C0/C0?C和变程A）分析城市表层土壤中各种金属元素的空间相关性及其变异规律、变异程度及产生变异的原因。利用变程的大小分析变量的可逆性，即空间依赖程度，从而从理论上认识表层土壤中重金属含量在二维平面空间上的空间变化特征及其影响因素。

利用半方差函数理论模型，按照半方差函数理论模型的要求，首先对原始数据采用域法识别特异值的方法进行处理，即样本平均值加3倍标准差，在此区间以外的数据均定为特异值，然后分别用正常的最大和最小值代替特异值，在半方差分析中将这8种元素的浓度进行对数转换，连同采样点的地理坐标输入软件GS+9.0，拟合半方差函数，并选择最佳拟合模型及其参数。对重金属含量、重金属空间结构特征进行分析，绘制出土壤中8种主要重金属含量的Krining插值图；

3

最后，将插值图叠加，得出污染最严重的区域，即为污染源的位置。 2.4关于问题4的分析

在问题3的基础上，对半方差函数理论模型进行分析评价。该模型在建立的过程中，未考虑到采样时间对重金属污染物的传播影响，为更好地研究城市地质环境的演变模式，还应当了解该城市不同时间的土壤重金属含量值，进一步建立模型。

3.符号说明

符号 涵 义 Pi Ci Si 采样点i污染物的污染指数 i污染物的实测值 i污染的评价标准 采样点的综合污染指数 P综P平均 采样点所有单项污染指数的平均值 采样点所有单项污染指数的最大值 P最大G?r? 累计贡献率 4.模型假设

⑴ 题目中给的采样点信息量足够大，能客观的反映现实情况； ⑵ 忽略地形、植物、空气对重金属污染的影响；

⑶ 在短时期内，城区各采样点的土壤重金属含量保持稳定； ⑷ 不考虑PH值对土壤金属污染含量的影响。

5.模型建立与求解

5.1各重金属元素在城区的空间分布及城区不同区域重金属污染程度模型（第1问）

5.1.1主要重金属元素在该城区的空间分布模型I

由于所给数据的样点是随机分布，且数量比较大，因此在作空间分布图时，采用Krining空间数据插值法，利用GIS软件平台ArcGIS9.3中普通Krining法对区域内变量进行空间拟合，得出所给8种主要重金属元素在该城区的空间分布（图5—1）。

4

5

6

7

图5—1 8种主要重金属元素在该城区的空间分布

注：图5—1中离散点为采样点，线条为等高线，填充颜色由浅到深表示浓度依次增加。

8

5.1.2 城区内不同区域重金属的污染等级模型Ⅱ

采用国家土壤环境质量标准（GB15618——1995）中的一、二、三级标准（表5—1）分别作为土壤重金属污染评价的污染积累起始值、中度污染起始值和重度污染起始值，把土壤重金属污染分为4个等级（表5—2）。建立单项污染指数和综合污染指数模型，评价不同城区8种元素的污染等级。具体方法如下： ⑴ 评价标准

表5—1 国家环境质量各级标准值（μg/g）表[2]

项目 As Cd 一级标准 15 0.2 二级标准 25 0.3 三级标准 30 1 ⑵ 评价方法 单项污染指数计算：

Cr 90 300 400 Cu 35 100 400 Hg 0.15 0.5 1.5 Ni 40 50 200 Pb 35 300 500 Zn 100 250 500 则P若Ci?Si污染起始值，i?Ci/Si污染起始值；若Si污染起始值?Ci?Si中度污染起始值，

则Pi?1??Ci?Si污染起始值?/?Si中度污染起始值?Si污染起始值?；

若Si中度污染起始值?Ci?Si重度污染起始值，

则Pi?2??Ci?Si中度污染起始值?/?Si重度污染起始值?Si中度污染起始值?；

若Ci?Si重度污染起始值，则Pi?3??Ci?Si重度污染起始值?/?Si重度污染起始值?Si中度污染起始值?

综合污染指数计算：P综=???P2平均+P2最大?/2??12

表5—2 土壤重金属污染级别分类表[1]

污染等级 1 2 3 4 Pi Pi?1 P综 P综?0.85污染水平 清洁 土壤轻度污染 1?Pi?2 0.85?P综?1.71 土壤中度污染 2?Pi?3 1.71?P综?2.56Pi?3 P综?2.56土壤重度污染 9

在EXCEL中运用上述算法处理附件1和附件2的数据，对模型进行求解，得到不同区域8种重金属元素的污染等级图表（表5—3.图5—2.）。

表5—3 不同区域重金属污染等级 As Cd Cr Cu Hg Ni Pb Zn 生活区 0.615 2.143 4.584 1.914 1.498 0.664 2.169 8.965 污染等级 1 3 4 3 2 1 3 4 工业区 1.243 2.581 1.421 7.214 10.649 0.899 2.079 5.471 污染等级 2 4 2 4 4 2 3 4 山区 0.552 1.626 1.032 1.130 0.844 1.553 1.113 1.414 污染等级 1 2 2 2 1 2 2 2 主干道路区 2.157 2.999 5.820 4.471 12.392 1.875 1.338 11.408 污染等级 3 4 4 4 4 3 2 4 公园绿地区 0.625 2.330 0.805 1.615 2.046 0.581 1.440 4.714 污染等级 1 3 1 2 3 1 2 4 14 12 10 8 6 4 2 0 生活区 工业区 山区 主干道路区 公园绿地区 As Cd Cr Cu Hg Ni Pb Zn 图5—2 不同区域重金属污染等级

由表5—3和图5—2可知城市不同区域金属污染含量的污染程度： 生活区：Zn>Cr>Pb>Cd>Cu>Hg>Ni>As 工业区：Hg>Cu>Zn>Cd>Pb>Cr>As>Ni 山区：Cd>Ni>Zn>Cu>Pb>Cr>Hg>As 主干道路区：Hg>Zn>Cr>Cu>Cd>As>Ni>Pb 公园绿地区：Zn>Cd>Hg>Cu>Pb>Cr>As>Ni

综合各区域的重金属污染程度，显然该城区中Zn的污染程度最严重，Hg、Cr、Cd的污染程度较为严重。

5.2分析重金属污染的主要原因模型Ⅲ（第2问）

10

5.2.1 土壤中重金属之间的相关性分析

在EXCEL中对数据处理得到城区土壤重金属含量的相关系数见表5—4.

表5—4 城区土壤重金属含量的相关系数表

元素 As Cd Cr Cu Hg Ni Pb Zn As 1 0.255 0.189 0.16 0.064 0.317 0.29 0.247 Cd 0.255 1 0.352 0.397 0.265 0.329 0.66 0.431 Cr 0.189 0.352 1 0.532 0.103 0.716 0.383 0.424 Cu 0.16 0.397 0.532 1 0.417 0.495 0.52 0.387 Hg 0.064 0.265 0.103 0.417 1 0.103 0.298 0.196 Ni 0.317 0.329 0.716 0.495 0.103 1 0.307 0.436 Pb 0.29 0.66 0.383 0.52 0.298 0.307 1 0.494 Zn 0.247 0.431 0.424 0.387 0.196 0.436 0.494 1 从重金属元素之间的相关系数可以看出，Cr和Ni相关系数最大，达到了0.715，说明它们之间关系最为密切，含量受彼此影响较大，来源大致相同，Pb和Cd、Cr和Cu及Cu和Pb相关性也较密切，说明来源极为相近。 5.2.2 土壤中重金属的主成分分析[4]

利用主成分分析的方法分析土壤中重金属污染的来源，需要4个步骤。 步骤⑴： 计算相关系数矩阵

?r11?r21R??????rp1?r12r22?rp2?r1p??r2p???????rpp??

式中rij?i,j?1,2,…,p?为原变量的xi与xj之间的相关系数，其计算公式为

nrij?

?(xk?1ki?xi)(xkj?xj)2?(xk?1nki?xi)?(xk?1nkj?xj)2

因为是R实对称矩阵（即rij?rji），所以只需计算上三角元素或下三角元素即可。

步骤⑵：计算特征值与特征向量

11

首先解特征方程?I?R?0，通常用雅可比法（Jacobi）求出特征值

?i(i?1,2,…,p)，并使其按大小顺序排列，即?1??2?……??p?0；然后分别求

2出对应于特征值?i的特征向量ei(i?1,2,…,p)。这里要求ei=1，即?eij?1，其

j?1p中eij表示向量ei的第j个分量。

步骤⑶：计算主成分贡献率及累计贡献率 主成分的贡献率为

zi??i??k?1p?1,2,…,p?

k累计贡献率为

G?r??????k?1k?1pik(i?1,2,…,p)

k一般取累计贡献率达85—95%的特征值?1,?2,…,?m所对应的第1、第2，……，第m?m?p?个主成分。 步骤⑷: 计算主成分载荷

其计算公式为

lij?p(zi,xj)??ieij(i,j?1,2,…,p)

得到各主成分的载荷以后，还可以进一步计算，得到各主成分的得分

?z11?zZ??21????zn1

z12z22?zn212

?z1m??z2m???????znm?

利用MATLAB编程（见附录1）可计算出变量的特征值、贡献率，如表5—5.图5—3.

表5—5 特征值、贡献率

数值 1 2 3 4 5 6 7 8

特征值4321012345678特征值 3.5602 1.1505 0.9653 0.7673 0.5778 0.4315 0.3015 0.2459 贡献率/(﹪) 44.5025 14.3813 12.0663 9.5913 7.2225 5.3938 3.7688 3.0708 累计贡献率/（﹪） 44.5025 58.88375 70.95 80.54125 87.76375 93.1575 96.92625 100 特征值

图5—3 特征值曲线

主成分分析一般要求累计的贡献率G?r??85%，对总方差贡献率小于5%的主成分可以不考虑。从表5—5、图5—3看出：共有2种主成分影响着具体沉积物的污染状况，即可以用2个化学过程来说明该区域的污染成因。

表5—6 因子载荷

主成分载荷 元素 As Cd Cr Cu Hg Ni Pb

主成分1 0.426 0.711 0.735 0.756 0.408 0.723 0.764 13

主成分2 -0.200 0.281 -0.444 0.125 0.673 -0.515 0.314 Zn 0.699 -0.037 由表5—5表明，重金属As、Cd、Cr、Cu、Ni、Pb、Zn在主成分1上类聚，重金属Hg在主成分2上类聚。As、Cd、Cr、Cu、Ni、Pb、Zn在主成分1上有较高的载荷，As、Cd、Cr、Cu、Ni、Pb、Zn平均值均高于土壤背景值，说明成土母质对土壤中As、Cd、Cr、Cu、Ni、Pb、Zn不起控制作用，因此认为主成分1主要来源于自然界以及人为因素，而不是仅仅来源于自然界。主成分2主要包括Hg元素，而Hg元素的平均值也高于土壤背景值，说明除受成土母质作用外，人为因素的输入也对其产生了明显的影响。

Cd、Zn在主干道路区达到最高值，说明Cd、Zn主要来源于汽车尾气排放、轮胎的老化磨损，Cu和Zn还来源于汽油、车体的磨损等。交通活动产生的灰尘经外力扬起沉积在周围土壤中，导致土壤中Zn和Cu含量升高，说明Cd、Zn、Cr可能主要受交通运输活动的影响。初步断定Zn、Cd主要来源于交通活动，Cu可能来源于工业活动。因此重金属的主要来源于交通及工业活动，属于人为来源。 5.3确定重金属污染源的位置模型Ⅳ（第3问）

首先，采用域法识别特异值，对原始数据进行处理，得到符合正态分布的合理数据（见附录2）。

然后，借助GS+9.0软件系统完成实验半方差函数和理论模型间的自动拟合、最优选择、参数计算及半方差图绘制，对符合正态分布的8种重金属的半变异拟合模型如图5—4.

Zn As

Cd Cr

Cu Hg

14

Ni Pb

图5—4 8种重金属的变异函数模型

再根据拟合的变异函数模型，可以统计得到各项参数，对重金属含量、重金属空间结构特征进行分析。利用GS+9.0软件自动拟合模型，得到表层土壤重金属含量的半方差函数理论模型相关参数见表5—7.

表5—7土壤重金属含量的半方差函数理论模型相关参数 元素 预测模型 C0 C0?C C0/C0?C有效变程A Rss 决定系数 As 线性 15640.3012 0.465 0.917 Cd 线性 15640.3012 0.000178 0.66 Cr 线性 15640.3012 164719 0.271 Cu 线性 15640.3012 333684 0.738 Hg 线性 15640.3012 0.000244 0 Ni 线性 15640.3012 42.2 0.762 Pb 线性 15640.3012 181096 0.723 Zn 线性 15640.3012 1.28E+08 0.565 进而，利用GS+9.0软件绘制出土壤中8种主要重金属含量的Krining插值图。 As

4.128 6.106 0.342 0.034 0.047 0.282 717.176 933.090 0.231 1365.84 2212.359 0.383 0.086 0.086 0 35.584 45.762 0.222 786.981 1387.157 0.433 23305.55 34566.939 0.326 15

Cd

Cr

Cu

16

Hg

Ni

Pb

17

Zn

图5—5 土壤中8种主要重金属含量的Krining插值图

注：以上图中不同颜色表示污染程度不同，白色区域表示污染最为严重的区域。

表5—8 各重金属污染程度较高的采样点 元素 编号 As 154、237、7、250、12、138 Cd 33、246、253 Cr 32 Cu 163、44、14 Hg 157、225 Ni 189、153、36、150、228、250、35、80 Pb 223、224、15、239、40 Zn 43 将以上8种重金属的Krining插值图叠加，可以得出污染最严重的区域的范围是：以采样点（2883,3617,15）、（4043,1895,14）、（3520,4357,7）、（5062,4339,5）

18

围成的区域，即为污染源的位置。

5.4为更好地研究城市地质环境的演变模式，进一步建立模型Ⅴ （第4问） 在问题3的基础上，对半方差函数理论模型进行分析评价。为了更好地研究城市地质环境的演变模式，分析重金属污染物的传播特征，现考虑采样时间对元素污染程度的影响，应在不同时间段进行样品采样，建立数理模型，得到改进后的重金属污染程度与时间的函数关系。

6.模型的优缺点

6.1模型的优点：

⑴ 兼顾了单因子污染指数的平均值和最高值，突出污染较重的污染物的作用，给予其较大的权值，能较为全面的反映沉积物的总体质量；

⑵ 采用域法识别特异值的方法处理数据，排除特异值，减小模型误差；

⑶ Krining是一种优化的插值方法，可以广泛的应用到地质、土壤、水文等领域；

⑷ 用重金属含量的空间结构特征分析重金属的传播特征，方法合理。 6.2模型的缺点：

⑴ 未能考虑重金属元素的毒性差别及其和不同母岩中重金属元素的浓度差别，评价时未能识别其污染机理；

⑵ 忽略了时间、空气等因素对土壤中重金属污染的影响。

7.模型的改进方向

在问题2中，应用理统计方法，分析土壤中重金属污染的主要原因时，利用重金属元素之间相关性分析和主成分分析法的同时，可以进行土壤重金属的聚类分析，进一步验证上述方法，从而保证模型的合理准确性。

参考文献

[1] 尹君.基于GIS绿色食品基地土壤环境质量评价方法研究[J].农业环境保护，

2001,20(6).

[2] 黄国峰，吴启堂.绿色食品产地土壤环境质量现状评价标准的修正[J].农业

环境保护，2000，19(2).

[3] 李保杰,于法展，纪亚洲.徐州市九里矿区土壤重金属插值分析及污染评价

[A].测绘科学，2010.

[4] 王海东，方凤满，谢宏芳，王翔，黄春海.芜湖市区土壤重金属污染评价及

来源分析[A].城市环境与城市生态，2010.

[5] 韩忠庚.数学建模竞赛：获奖论文评选与点评/韩忠庚主编.北京：科学出版

社，2007.

19

2的MATLAB程序3的数据处理

20

附录

附录一：问题 附录二：问题

本文来源：https://www.bwwdw.com/article/yqzr.html