2011数学建模国家一等奖论文A题

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

摘要

随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出，本文就八种重金属元素在城市不同位置（区域）的污染浓度数据等进行建模分析，讨论这八种重金属元素对城市不同功能区的污染程度和污染产生的原因，从而确定污染源，更好地解决问题。

问题一，我们首先对附件中的数据进行处理，应用单因子污染指数法、尼梅罗综合污染指数法及地质累积指数法从不同角度对城市土壤重金属污染程度进行了不同指标下的评价。在得出各取样点的重金属污染指数后，即能确定这八种重金属在该城区的空间分布，同时由所给的各取样点所处的区域，得出各区域重金属的污染程度。总的来说，工业区污染较为严重，山区重金属污染程度较低，且该城市受到Hg、Zn、Cu的污染较为严重，受到Ni、As的污染程度较弱。

问题二，每个土壤样品的重金属元素组成并不是各自、互不相干的，来自同一种源的那些成分之间存在较强的相关性，因此我们采用因子分析法来研究多个变量的相关性。通过将庞杂的原始数据按照成因上的联系进行归类，由果及因地理出几条比较客观地成因线索，以导出正确的成因结论。本文在土壤环境质量评价结果基础上将城市土壤八种重金属元素的数据系列使用分析软件SPSS进行因子分析，进而判别城市重金属的污染来源及其分布规律。最后得出结论：该城市重金属污染来源主要有四个方面：交通污染源，燃煤污染源，工业污染源，冶金农药污染源。

问题三，针对重金属元素在土壤中的扩散问题，我们采用抛物型的偏微分方程模型，我们可以解出这个偏微分方程的Cauchy问题的解。由于污染源的污染范围是有限的，我们利用前两问找出某种元素污染较为严重的某一区域，利用反演的方法拟合出污染源的位置(x0,y0,z0)，并用此方法找出其他元素的污染源位置。我们主要对污染较为严重的Zn、Cd、Hg元素进行研究，得出了重金属元素Zn、Cd、Hg的一个或几个污染源位置。

问题四，我们首先提出了模型的优缺点，并有针对性的提出模型的改进意见。考虑到城市环境地质演变模式的相关因素，我们提出了一些需要收集的信息，在信息量相对充足的情况下，建立相应的数学模型。例如我们考虑到大气干湿沉降对重金属含量的影响，建立了与大气干湿沉降速率有关的数学模型,可以利用该模型粗略预测t年后在重金属污染下的城市环境地质状况。

关键词：单因子污染指数法、尼梅罗综合污染指数法、地质累积指数法、因子相关性分析法、抛物型的偏微分方程模型、反演、城市环境地质演变模式

一、问题重述

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

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

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

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

现要求你们通过数学建模来完成以下任务：

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

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

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

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

二、模型假设

1、综合评价区域内重金属污染状况时，假设污染程度只受题目中8个污染物的影响。

2、假设附件中所给数据都是真实的、可靠的。

3、为了模型的准确性，假设5个功能区相互之间无交叉。

4、为了保证对重金属污染来源测定的准确性，我们假设源到采样点之间，污染物在途中保持质量守恒。

5、考虑到使用因子分析法的基础，我们假设受体污染物中某种元素是各污染源贡献的线性组合，各污染源是之间互不相关的。

6、假设各个污染源贡献的某元素的量（称因子载荷）有足够的差别，而且它在采样点和分析期间变化不大。

7、假设土壤的密度是均匀的。

三、符号说明

——单项污染分指数                    

——评价对象中污染物的实测含量值     

——重金属累积起始值

——重金属污染起始值

——重金属重污染起始值

——尼梅罗综合污染指数

——所有元素污染指数最大值

——所有元素污染指数平均值

——地质累积指数

——沉积岩中的地球化学背景值

——t时刻点处的一种元素的浓度

——分别是沿x,y,z方向上的扩散系数

——衰减系数

——取样点的横坐标、纵坐标、高度及一种元素的浓度

四、问题一 重金属污染分布

4.1问题分析

问题1要求我们给出八种重金属元素在该城区的空间分布以及不同区域内重金属的污染程度，附件1给出了取样点的位置以及所属的功能区，附件2给出了八种重金属在相应位置的浓度。

目前国内外污染评价的方法较多，最具代表性的有：指数法、模糊综合评价法、层次分析法等，其中指数法包括：单因子污染指数法、综合污染指数法，地质累积指数法，在这里，我们主要采用指数法。

由各指数法的污染程度计算方法，通过计算确定各位置八种重金属元素在各指数法下的指标值，与各指数法的污染分级指标相比较，从而确定各位置八种重金属元素的污染程度。

经过初步的分析与计算，我们得到该城市整体重金属含量特征参数统计表格：

表1  该城市整体重金属含量特征参数统计表

元素	采样点 （个）	背景值 (μg/g)	最大值 (μg/g)	最小值 (μg/g)	平均值 (μg/g)	标准差
As	319	3.6	30.13	1.61	5.68	3.02
Cd	319	0.13	1.61	0.044	302.40	0.22
Cr	319	31	920.84	15.32	53.51	69.
Cu	319	13.2	2528.48	2.29	55.02	162.66
Hg	319	0.035	16.00	0.008	0.30	1.63
Ni	319	12.3	142.50	4.27	17.26	9.93
Pb	319	31	472.48	19.68	61.74	49.98
Zn	319	69	3760.82	32.86	201.20	338.70

4.2数据的处理

①首先将重金属元素浓度单位统一为μg/g。

    ②根据单因子污染指数法的计算方法，求出各重金属元素在不同浓度范围内的污染指数，根据标准得出各区域在八种重金属单项污染的土壤质量等级。

    ③由上述的，得出所有元素污染指数最大值以及所有元素污染指数平均值，根据综合污染指数法的计算公式，求出尼梅罗综合污染指数，从而确定各区域在八种重金属污染下的综合污染土壤质量等级。

    ④根据地质累积指数法的计算公式以及附件3中各重金属元素的背景值，求出各重金属元素地质累积指数，从而确定各区域在某一种重金属元素污染下的地质累积污染程度强弱状况。

    ⑤根据附件一中所给取样点三维空间坐标进行插值处理，得出该城区各功能区的分布示意图（具体程序见附录2）：

图1  各功能区分布示意图

4.3模型的建立

   4.3.1单因子污染指数法

   单因子污染指数是指在某一污染物影响下的环境污染指数可以分别反映各污染物对环境的污染程度。本文根据下式计算各重金属元素的单因子污染指数，按单因子污染指数分级标准对重金属元素污染程度进行分级：

    

式中：—单项污染分指数；—评价对象中污染物的实测含量值；—重金属累积起始值；——重金属污染起始值；——重金属重污染起始值

而在城市表层土壤污染评价中，本文以国家土壤环境质量标准的一级标准作为重金属元素累积起始值（），国家土壤环境质量标准的二级标准作为污染起始值（），国家土壤环境质量标准的三级标准作为重污染起始值（）

表2：国家标准单因子污染指数评价各污染等级起始值（mg/kg）

元素	累积起始值()	污染起始值()	重污染起始值()
Cd	0.15	0.30	1.0
Hg	0.15	0.30	1.5
As	15	25	30
Pb	35	250	500
Cu	30	50	400
Cr	90	200	300
Zn	85	200	500
Ni	40	50	200

表3：城市表层土壤单因子污染指数分级标准

污染分指数	<1	1<2	2<3	3
质量等级	清洁	潜在污染	轻污染	重污染

4.3.2尼梅罗综合污染指数法

前文虽对单个元素对城市表层土壤的污染程度做了评价分析，但土壤中各种重金属元素常常相伴而生，共同对土壤造成污染，因此前文针对单个元素的污染评价并不能反映各元素的共同作用对城市表层土壤环境造成的复合污染。鉴于此为突出环境要素中含量最大的污染物对环境质量的影响，采用尼梅罗综合污染指数法对研究区城市近地表灰尘及表层土壤重金属污染程度进行综合评价，计算公式为：

                     

式中：——尼梅罗综合污染指数；——所有元素污染指数最大值；——所有元素污染指数平均值。尼梅罗综合污染指数既反映了各种重金属对城市表层土壤的作用，同时又突出了高含量重金属对城市表层土壤环境质量的影响。尼梅罗综合污染指数分级标准见表4：

表4：尼梅罗综合污染指数分级标准

等级	Ⅰ	Ⅱ	Ⅲ	Ⅳ	Ⅴ
综合指数	0.7	0.7<1.0	1.0<2.0	2.0<3.0	3.0
污染等级	清洁	警戒限	轻度污染	中度污染	重污染

4.3.3 地质累积指数法

前面虽然考虑了人为污染因素对城市表层土壤环境造成的破坏，但却没有考虑由于自然成因作用可能引起背景值变动的因素，而从环境地球化学的角度出发的地质累积指数法注意到了这一点，因此，本文采用地质累积指数法对城市土壤重金属污染程度进行了评价，以揭示自然地质过程和人为活动共同对重金属污染的影响。其计算公式为：

                        

其中—地质累积指数，—元素在灰尘及土壤中的实测含量值，—沉积岩中的地球化学背景值，1.5为修正指数，是考虑到由于成岩作用可能会引起背景值的变动，通常用来表征沉积特征、岩石地质等其它影响。

按受污染程度强弱，将地质累积指数分为7个级别（表4），0-6级表示污染程度由无到极强，最高一级（7级）的元素含量可能达背景值的几百倍,见表5：

表5：地质累积污染指数分级标准

	级别	污染程度
<0	0	无污染
0-1	1	无污染到中度污染
1-2	2	中度污染
2-3	3	中度污染到强污染
3-4	4	强污染
4-5	5	强污染到极强污染
>5	6	极强污染

4.4模型的求解

4.4.1重金属元素的空间分布

依据附件中所给数据，利用MATLAB编程求解（具体程序见附录2），得出这8种重金属元素在该城区的空间分布状况的等值线图，见如下8个图解：

图2  As的分布状况                          图3  Cd的分布状况

                             

图4  Cr的分布状况                     图5  Cu的分布状况

图6   Hg的分布状况                      图7  Ni的分布状况

图8  Pb的分布状况                         图9  Zn的分布状况 

4.4.2不同区域重金属污染程度分析

①依据单因子污染指数法的计算公式，我们得出各元素的单因子污染指数，应用单因子污染指数法对城市土壤中Hg、Cd、Pb、Cu、Cr、Zn、As、Ni八种重金属进行污染程度评价，具体评价时我们采用单因子污染指数的平均值进行评价，具体见下表：

表6：城市表层土壤重金属单因子污染指数统计结果I

	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
最大值	3.0260	3.8854	9.2084	9.0814	15.0833	1.6167	2.89	13.8694
最小值	0.1073	0.2667	0.1702	0.0763	0.0571	0.1068	0.5623	0.3866
平均值	0.3838	1.51	0.5835	1.2246	0.7571	0.4250	1.0827	1.5570
污染 程度	清洁	潜在 污染	清洁	潜在 污染	清洁	清洁	潜在 污染	潜在 污染

表7：城市表层土壤重金属单因子污染指数统计结果Ⅱ

污染程度	各污染等级采样单元数(个)
	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
清洁	314	83	316	172	260	314	88	107
潜在 污染	4	113	3	60	35	5	228	134
轻污染	0	116	0	85	16	0	3	60
重污染	1	7	0	2	8	0	0	18

从表7可见，从各元素的平均单因子污染指数来看，城市表层土壤受重金属污染并不十分严重，仅存在受Cu、Zn和Cd的潜在污染。但平均值仅能反映全区的平均状况，掩盖了局部区域受污染的程度。从单个采样点的评价结果统计（表6）来看，部分采样单元的城市表层土壤受Cd、Cu、Hg、Zn的污染程度仍然较高。在所有采样单元中，受轻污染以上程度污染的采样单元：Cd有7个，Hg有8个，Cu有3个，Zn有18个

以上我们给出了整个城市表层土壤重金属单因子污染指数统计结果，下面我们将给出5个功能区城市表层土壤重金属单因子污染指数统计结果：（由于数据较多，不便在正文中给出，具体表格见附录1）：

由这个5个表可得出：工业区重金属污染状况最为严重，且在该区内，Zn元素相对其他几个元素而言污染较为明显；其次是生活区和交通区，有四个元素污染程度处于潜在污染的情形。而山区的重金属污染状况最好，各个元素的污染程度都处于清洁情形。

②依据尼梅罗综合污染指数法的计算公式，我们得出5个功能区的尼梅罗综合污染指数（为了比较简便的说明污染程度，我们对各区域内所有采样点的求出平均值），按照尼梅罗综合污染指数分级标准确定不同区域重金属的污染程度，见下表:

表8不同区域重金属的污染程度

功能区	生活区	工业区	山区	交通区	公园绿地区
	1.6538	2.0117	0.9268	1.8876	1.3478
污染 程度	轻度污染	中度污染	警戒限	轻度污染	轻度污染

由表8可看出工业区污染较为严重，山区重金属污染程度较低，同时由城区污染示意图(图10)可得出>2即中度污染及其以上的取样点集中在交通区和工业区附近， 2即中度污染以下的取样点集中在山区和公园绿化区附近，这可能是由于工业区污染物：废气、废水等的排放较多，导致附近地区的土壤中重金属污染程度较高：而山区重金属污染程度较低，原因可能是:一是由于山区受人类活动的影响小，二是山区在工业区的上风向且远离工业区，不易受到工业污染物的影响 。

图10  尼梅罗综合污染指数法下的城区污染示意图

③依据地质累积指数法的计算公式，我们得出每个重金属元素在各取样点的，在这里给出319个取样点在地质累积指数法下的污染状况，见表9：

表9：城市表层土壤重金属地质累积污染指数评价结果

污染程度	取样点数（个）
	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
无污染	161	119	198	88	167	211	161	151
无污染到中度污染	148	126	104	119	68	103	108	92
中度污染	7	59	12	73	46	3	45	52
中度污染到强污染	3	14	3	29	21	2	2	13
强污染	0	1	0	8	5	0	3	8
强污染到极强污染	0	0	2	0	5	0	0	2
极强污染	0	0	0	2	7	0	0	1

由上表可得出:取样点中受到Hg的污染最为严重，有7个取样点受到Hg的极强污染；其次是Zn的污染较为严重，中度污染及其以上的取样点较多；相对而言，取样点中Ni的污染程度最低，有211个取样点处于无污染的程度，而强污染、强污染到极强污染、极强污染的个数都为0；其次是As的污染程度较低，有161个取样点处于无污染的程度，而强污染、强污染到极强污染、极强污染的个数也都为0。由于这些取样点是根据一定标准所选取的，根据这些点污染状况由部分到整体的思想，可得出该城市在地质累积指数法下的污染状况，即该城市受到Hg、Zn、Cu的污染较为严重，受到Ni、As的污染程度较弱。

五、问题二 重金属污染的主要原因

5.1问题分析

问题2要求我们通过数据分析，说明重金属污染的主要原因，由问题1得出了八种重金属元素在该城区的空间分布以及不同区域的重金属污染程度，在不同区域内，重金属污染的原因不同，即重金属的来源不同，但是也有可能不同的重金属来自同一污染源，从而每个土壤样品的重金属元素组成并不是各自、互不相干的，来自同一种源的那些成分之间存在较强的相关性。因此我们想到了使用因子分析算法来解决这个问题，通过因子分析法并采用SPSS软件得出各个因子空间分布的等值线图，将此等值线与城市各功能区空间分布图，确定各因子集中的区域，即可确定该因子的属性，从而得出重金属污染原因。

5.2模型建立和求解

通过对因子分析算法的了解，我们发现该城市重金属元素含量的数据特征符合因子分析的要求。在对该城市土壤环境质量评价基础上，对存在污染的319个土壤采样点其污染来源及分布规律进行分析。以Hg、Cd、Pb、As、Cu、Cr、Zn、Ni八种重金属元素为指标作因子分析，从而为长春市表层土壤重金属污染成因的解释提供一定的理论依据。运用SPSS13.0统计软件进行因子分析，得出以下结论：首先给出该城市表层319个土壤采样点Hg、Cd、Pb、As、Cu、Cr、Zn、Ni八种重金属原始数据的相关系数矩阵：

表10：变量相关系数矩阵

指标	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
As	1.000	0.255	0.1	0.160	0.0	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.1	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.0	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

由表可见，Ni和Cr的相关性最好，相关系数最大为0.716。其次为Pb和Cd，相关系数为0.660，以下依次是Cu和Cr相关系数为0.532，其它元素之间的相关性并不好。从成因上来分析，相关性较好的元素可能在成因上和来源上有一定的关联。

因子分析的关键就是利用相关系数矩阵求出相应的因子特征值和累计贡献率，用SPSS13.0统计软件计算可得出，结果见表11：

表11特征值和累计贡献率

因子	旋转前			旋转后
	特征值	占总变量得百分率（%）	累计贡献率（%）	特征值	占总变量得百分率（%）	累计贡献率（%）
F1	3.560	44.500	44.500	2.097	26.212	26.212
F2	1.150	14.377	58.877	2.067	25.835	52.047
F3	0.965	12.063	70.941	1.258	15.721	67.768
F4	0.768	9.596	80.537	1.021	12.769	80.537

在累积方差为80.537%的前提下，分析得到4个主因子，而且从表11可以看出旋转前后总累计贡献率没有变化，即总信息量没有损失。从表11还可以得出，旋转之后，主因子F1、主因子F2以及主因子F3的方差累计贡献率就占70.941%，这可以解释为因子F1、因子F2、因子F3可能为该城市土壤重金属主要污染因子，而因子F4对该城市重金属污染相对而言作用较小。

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

表12：旋转前因子载荷矩阵

元素	因子
	F1	F2	F3	F4
As (μg/g)	0.426	-0.200	0.681	0.551
Cd (ng/g)	0.711	0.281	0.282	-0.322
Cr (μg/g)	0.735	-0.444	-0.303	-0.046
Cu (μg/g)	0.756	0.125	-0.365	0.137
Ni (μg/g)	0.723	-0.515	-0.190	0.137
Pb (μg/g)	0.7	0.314	0.237	-0.248
Zn (μg/g)	0.699	-0.037	0.123	-0.241
Hg (ng/g)	0.408	0.673	-0.297	0.449

表13：方差极大正交旋转后因子载荷矩阵

元素	因子
	F1	F2	F3	F4
As (μg/g)	0.134	0.173	0.026	0.969
Cd (ng/g)	0.131	0.851	0.131	0.096
Cr (μg/g)	0.878	0.242	0.034	0.009
Cu (μg/g)	0.577	0.330	0.546	-0.025
Ni (μg/g)	0.877	0.154	0.036	0.218
Pb (μg/g)	0.170	0.840	0.228	0.117
Zn (μg/g)	0.399	0.628	0.020	0.096
Hg (ng/g)	0.003	0.146	0.942	0.035

变量与某一个因子的联系系数绝对值（荷载）越大，则该因子与变量关系越近。表12和表13的正交因子解说明：因子F1为Cr、Ni、Cu的组合，因子F2为Cd、Pb、Zn的组合，因子F3为Hg，因子F4为As。

为了更好的进行分析、评价，利用因子得分所得到的4个因子经过方差极大正交旋转后的表层土壤重金属污染严重的319个单点样本在4个主因子上的得分表（由于数据量较大，论文中不便给出），可作出各个因子空间分布的等值线图，更能直观地说明各个元素在空间平面上的分布特征：

  图11 因子F1等值线                          图12  因子F2等值线 

图13  因子F3等值线                           图14  因子F4等值线

结合上面的4个主因子空间分布的等值线图、图1（城市各功能区分布示意图）以及各重金属元素在该城区的空间分布状况的等值线图，可得出以下结论：

主因子F1以Cr、Ni、Cu元素的载荷为主，这与前面重金属之间的相关分析结果相一致，且与城市功能区分布示意图相比较得出主因子F1主要分布在工业区与交通区相交区域，且Cr、Ni这两种重金属元素主要来自工业上冶金所产生的污染，Cu来自汽车尾气和轮胎磨损产生的粉尘。即可得出一个结论：主因子F1主要来自工业与交通的污染。

主因子F2以Zn、Pb、Cd元素的载荷为主，Pb和Zn常被作为机动车污染源的标识元素，所以该组元素主要来自于交通运输过程，被称为“交通因子”。从因子F2图可以看出，交通因子在空间分布上近似可认为是一个带状的污染源，分布比较广泛，这主要是因为Pb、Zn主要来自汽车尾气和轮胎磨损产生的粉尘。尤其是在宽城区一带，该处交通发达，大部分一级市内交通线路都在这附近，与城市功能区分布示意图相比较，在生活区与交通区相邻位置，交通因子会达到一个高值点，说明该位置机动车流量大，汽车排放尾气量较大，造成重金属Zn、Pb等的污染较为严重。

主因子F3为Hg，Hg是该城区污染最严重的重金属之一，这在我国城市中非常常见，尤其是在我国北方城市中表现得尤其突出。其中最大的一个原因是燃煤污染，Hg污染的高峰集中在生活区和工业区，这说明居民采暖用煤和工业用煤成为表层土壤汞污染的主要原因。另外除了燃煤来源外，工业三废排放也是城市表层土壤Hg污染的另一个重要来源。

主因子F4为As，As主要来源于工业污染，工业上排放砷的主要有化工、冶金、炼焦、火力发电、造纸、玻璃、毛革、电子工业等，同时农药的使用也是造成As污染的一方面原因。由As的污染分布状况、城市功能区分布图以及根据附件中的数据得出的城市等高线图，可以看出As污染主要分布在工业区域生活区的交界处，同时在海拔相对周围较低的地区，这就可以解释As的污染分布状况是很明显的富集点，总是在海拔较低的地区达到一个高值点。

综上，得出以下结论：该城市重金属污染来源主要有四个方面：交通污染源，燃煤污染源，工业污染源，冶金农药污染源，其中Zn、Pb、Cu主要来自汽车尾气的排放，即交通污染源；Hg主要来自工业燃煤及居民采暖用煤；As污染主要来自化工、冶金、炼焦、火力发电、造纸、玻璃、毛革、电子工业等工业污染；Cd、Cr、Ni污染主要来自金属冶炼。

六．问题三污染源位置确定

6.1问题分析

问题3要求我们分析重金属污染物的传播特征，从而确定污染源的位置，由于污染物在土壤中并不是固定的，它们可以通过挥发、水土流失等形式进入大气和水体中，所以应考虑重金属元素在土壤中迁移、转化的规律，建立相应的重金属迁移、扩散模型，通过污染点来确定污染源位置。由于大多数与反应扩散有关的现象，都能由现象或非线性抛物型偏微分方程作为数学模型来定性或定量地加以解决，考虑到污染物在扩散过程中会有所挥发或流失，故我们采用由衰减的扩散过程的数学模型，且污染源对周边环境的影响范围是有限的。随后，我们利用反演的方法，逐步找出污染源的位置。

6.2模型建立

6.2.1有衰减的扩散过程数学模型

考虑重金属元素在土壤中的扩散和衰减合作用，由物质不灭定律得

                 (1)

    方程(1)是常系数线性抛物方程，它就是有衰减的扩散过程的数学模型，对于具体问题，尚需与相应的定解条件（初始条件与边界条件等）匹配才能求得确定情况下的解。在本体中就是要通过初始t=0时污染源的条件来进行求解。

6.2.2线性抛物方程的Cauchy问题的解

设污染源在点处，则此扩散问题满足Cauchy问题

                  (2)

其中为Dirac函数，即   

解上述方程得到

        (3)

6.2.3参数估计

在Cauchy问题(2)的解(3)中，有四个未知参数，它们分别是扩散与衰减过程中的扩散系数与衰减系数的算术平方根。令取样时刻的t=1（取样时间的选择不同只会造成参数值的不同，对问题没有实际影响）。对(3)式两边去对数，得

    (4)

令  

则(4)式变为

                              (5)

利用已知数据，可以求出的数据，用三元回归分析方法求出的估计值如下：

                                     (6)

其中                   (7)

满足方程组

                        (8)

                      (9)

由可求得的估计值，即  

又由于

                                  (10)

由（5）式可得，再把带入（9）式，得

                                  (11)

至此得到参数的估计值，把它们代入（2）式分别代替，则得不含未知参数的解得近似表达式。

6.3模型的求解：用反演的方法求污染源的位置

   1）通过第一问的分析，我们发现这一地区造成污染的重金属元素主要有Zn，Cd，Hg，这里主要讨论Zn元素的某一污染区域，讨论用反演的方法求污染源的位置。在图9中，我们发现Zn元素的污染位置主要集中在两块，其中一块的区域范围为。

①筛选在这个区域内的点有：(12855,45,18)、(13797,9621,18)、(14325,8666,23)、(121,9560,11)、(14000,70,14)、(14207,9980,14)、(14065,10987,25)、(12734,10344,32)。

②我们初步选取一个污染源位置(13797，9621，18）,将这个数据及t=1带入(3)式中，这时u是关于的方程且有参数；

③利用筛选过的观测点数据（留两个检验拟合度用）和上面讨论的参数估计法，计算出的值。；

④用那两个数据检验上一步得出的系数，发现误差较大，重新选取污染源位置重复上面的步骤，直至最后最小二乘误差小于1；

最后，得出近似污染源的位置为(13522.3,9448.0,23.9),此时。

2）类似的可求出Zn另一污染区域的污染源及Cd和Hg元素污染的污染源，如下表所示（具体程序见附录2）：

表14   Zn、Cd、Hg三种元素的污染源位置

元素	x(m)	y(m)	z(m)
Zn	13522.3	9448.0	23.9
	9458.2	4562.9	69.3
Cd	21456.4	11433.2	1.3
	2017.2	3112.3	12.3
Hg	15387.1	9190.0	25.0
	18522.8	3563.9	25.1
	2501.1	3095.5	48.2

通过与实际数据对比发现，这些污染源位置点一般不是一定区域内某重金属元素含量的极大值点，但是这些点距处在元素含量的较大值处，这与实际情况是相符的，该模型是可行的。

七、问题四 模型的优缺点

7.1问题分析

问题4要求我们分析上面所建立的模型的优缺点，同时为了更好地研究城市地质环境的演变模式，收集一些题目中未给出的信息，在这些信息下重新改进模型，完善模型，以便于更好地应用于其他城市地质模式。

7.2模型的优缺点

7.2.1优点：

①问题1中运用了单因子指数分析法、尼梅罗综合污染指数法、地质累积指数法对八种重金属元素在该城市的空间分布进行了分析，并给出了空间分布状况的等值线图，更可直观地看出分布状况。

②问题1中我们不仅考虑了单个因子对土壤的污染状况作出评价，同时对土壤重金属污染程度进行了综合评价；另外，我们不仅考虑了人为污染因素对城市表层土壤环境造成的破坏，同时运用地质累积指数法表明自然地质过程对重金属污染的影响。这样建立的模型较为完整地揭示了自然地质过程和人为活动共同对重金属污染的影响。

③问题2中运用了因子分析算法，考虑到每个重金属元素污染成因并不是相互的，运用SPSS软件建立每个重金属元素间的相关系数矩阵，利用相关系数矩阵求出了相应的因子特征值和累计贡献率，通过主因子的属性确定重金属元素的污染原因。优点在于将每个重金属元素联系起来而不是单独分析每个元素的产生原因，这样更贴近实际，得出的结果也相对较为准确。

③问题3中，得到的抛物型偏微分方程形式很好，每一个系数都有其代表的实际意义，符合重金属污染物在土壤中的迁移规律。

 7.2.2缺点

①本论文使用的所有数据仅限于表层土壤样品重金属元素的测试结果，而城市重金属污染评价实际上是一个全面而又复杂的过程，需要对多种样品中重金属含量进行综合研究，仅对表层土壤样品的数据进行分析评价，其评价结果不能十分客观和完全真实地反映城市的重金属污染情况。

②本论文所采用的模型仅对土壤环境质量现状作出评价，然而土壤环境质量的研究不单就某个时间，某个整体而言，应与其发展阶段和区域土壤类型相结合，也就是说土壤环境质量评价应注意其历史变化过程，不仅要对土壤环境质量现状作出评价，也要对土壤环境质量回顾和土壤环境质量预测影响作出评价。

③问题3中，利用反演的方法寻找污染源点肯定有所误差，且划定某一元素污染源的影响范围是往往只剩几个观测点，拟合数据较少也会造成误差较大。而且偏微分方程的解要求有初值，这对我们求解污染源点位置增加了难度。

7.3收集信息研究城市环境地质演变模式

地质环境主要指的是自地表面下的坚硬壳层，即岩石圈。地质环境是地球演化的产物。岩石在太阳能作用下的风化过程，使固结的物质出来，参加到地理环境中去，参加到地质循环以至星际物质大循环中去。

1）重金属元素不仅来自人类活动的影响，同时大气干湿沉降是重金属的主要来源，同时我们假设重金属大气干湿沉降速率此后保持不变，在上述假设条件下，我们需要知道各元素的年平均大气干湿沉降速率，分别对重金属大气干湿沉降对表层土壤中重金属含量的累积影响进行分析，利用以下的计算公式：

                  

其中——t年后表层土壤中重金属含量，——表层土壤中重金属现状含量；V——重金属大气干湿沉降速率；t——预测年限。利用该计算公式可得出t年后表层土壤中重金属含量，对土壤的重金属污染状况进行预测预警，即可粗略预测t年后在重金属污染下的城市环境地质状况。

2）以上我们研究的都是城市表层土壤的重金属的污染状况，但城市环境地质演变还应考虑到深层土壤中的重金属污染状况，所以我们需要知道上述八种重金属元素在深层土壤中的浓度数据，再利用问题一所述模型对各功能区深层土壤的重金属污染状况作物污染评价，即可得到深层土壤在各个指标下的污染状况，结合对比表层土壤的重金属污染状况可以更为完整地对城市环境地质演变模式作出分析。

3）若给出了城市主要工厂的具体位置或该工厂的污染排放是否达标等信息，则可以更为准确地分析该城市的污染原因及污染源的具体位置。

八、模型的改进

本次建模由于时间、精力等原因的有限，有很多可以改进却未来得及改进的地方：

1、问题一中对每种金属污染，在各个功能区做出其均值及方差，测试到的污染物的浓度普遍呈正态分布，比较同一种重金属污染物在不同功能区的分布，则看出每个功能区的主要污染物。

2、问题一中除了采用以上三种污染指数分析法外，还可采用潜在生态危害指数法，研究沉积物中重金属对环境的影响，该种方法不仅反映了某一特定环境中的各种污染物的影响，而且也反映了多种污染物的综合影响，并且用定量的方法划分出潜在生态危害的程度。

3、主因子分析的步骤：不仅需要对题目给出的所有数据进行主因子分析，最好还应有针对各个功能区的主因子分析（附录3给出了各个功能区的八个元素的相关性系数矩阵），根据各个功能区的主因子可以降低远距离地区因子的影响，提取出的主因子更能准确的反映各地的（污染）状况。

4、问题二的因子分析结束后，增加一个步骤，称为“聚类”：因子分析结束后每个采样点成为了一个“向量”(x,y,z,f1,f2…fk)，其中（x,y,z）为该采样点的地理坐标，（f1,f2,…fk）为该采样点的各主因子得分，对任意两个“向量”做出欧式空间距离，将这种距离相近的采样点归为同一类，最终将所有的点都归为几类，这样,对归结完的每一类，都是具有相近的地理位置（x,y,z）或相近的主因子（F1，F2,…Fk），从而在整体上能更好的解释污染现象

5、“聚类”操作完成后，考虑到聚类的标准（欧式空间距离大小），不同类具有不同的主要污染物，在此基础上建立不同的污染物传播模型。

九．参考文献

[1]杨忠平，长春市城市重金属污染的生态化学特征及其来源解析：159-163，2008

[2]陈高武，重庆都市圈土壤重金属元素迁移富集及生态效应研究： 73-77，2008

[3]翟航，长春市土壤重金属分布规律及土壤环境质量评价研究：42-43，48-50，59-69，2007

[4]王树禾，扩散问题的偏微分方程模型：69-75，1990

[5]纪书华，多孔介质中重金属反应性运移的数值模拟研究：30-34，2009

[6]卢纹岱，《SPSS统计分析》，北京：电子工业出版社，2003

[7]苏金明，阮沈勇，《MATLAB实用教程》，北京：电子工业出版社，2005

[8]姜启源，谢金星，叶俊，《数学模型（第三版）》，北京：高等教育出版社，2003

附录1：（城市各功能区重金属污染状况）

元素	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
最大值	0.7633	3.0636	7.4446	2.2681	2.2083	0.8200	2.89	10.9782
最小值	0.1560	0.5787	0.2051	0.3243	0.0800	0.2223	0.6980	0.5102
平均值	0.4180	1.6260	0.7418	1.2850	0.5817	0.4586	1.1273	1.73
标准差	0.1433	0.6088	1.0733	0.6832	0.5437	0.1416	0.3401	1.6387
污染程度	清洁	潜在污染	清洁	潜在污染	清洁	清洁	潜在污染	潜在污染

生活区重金属污染状况

工业区重金属污染状况

元素	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
最大值	1.6870	3.1327	2.8558	9.0814	2.2175	1.1700	2.7392	6.7534
最小值	0.1073	0.7633	0.1711	0.4233	0.0786	0.1068	0.26	0.6627
平均值	0.4929	1.9521	0.5842	1.7214	1.3617	0.4988	1.2590	2.0140
标准差	0.3152	0.5676	0.4412	1.3973	2.1949	0.2196	0.3796	1.3185
污染程度	清洁	潜在 污染	清洁	潜在 污染	潜在 污染	清洁	潜在 污染	轻污染

山区重金属污染状况

元素	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
最大值	0.7327	2.1537	1.1576	2.0545	1.3786	0.5200	1.3667	2.0993
最小值	0.1180	0.2667	0.1800	0.0763	0.03	0.1318	0.5623	0.3866
平均值	0.2696	0.9950	0.4303	0.5761	0.2730	0.3821	0.53	0.8393
标准差	0.1199	0.49	0.2593	0.3469	0.1857	0.2369	0.1702	0.2950
污染程度	清洁	清洁	清洁	清洁	清洁	清洁	清洁	清洁

交通区重金属污染状况

元素	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
最大值	3.0260	3.8854	9.2084	5.7567	15.0833	1.005	1.6813	13.8694
最小值	0.1073	0.3340	0.1702	0.4113	0.0571	0.1548	0.62	0.4814
平均值	0.3904	1.8002	0.6313	1.4450	0.9284	0.42	1.1070	1.8290
标准差	0.2907	0.6591	0.8136	0.7559	2.0858	0.1600	0.1902	1.4365
污染程度	清洁	潜在 污染	清洁	潜在 污染	清洁	清洁	潜在 污染	潜在 污染

公园绿地区重金属污染状况

元素	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
最大值	0.7787	3.0356	1.0571	2.2666	2.8661	0.7275	1.49	5.96
最小值	0.1847	0.80	0.1812	0.3013	0..667	0.1900	0.7683	0.4369
平均值	0.4176	1.4571	0.4845	0.9523	0.5931	0.3822	1.1028	1.3247
标准差	0.1349	0.6812	0.1635	0.4353	0.6116	0.1243	0.2285	0.9699
污染程度	清洁	潜在 污染	清洁	清洁	清洁	清洁	潜在 污染	潜在 污染

附录2

1)取样点在城市的分布图及城市等高线示意图程序:

clear;

WS=[7.84    153.80  44.31   20.56   266.00  18.20   35.38   72.35 

5.93    146.20  45.05   22.51   86.00   17.20   36.18   94.59 

4.90    439.20  29.07   .56   109.00  10.60   74.32   218.37 

6.56    223.90  40.08   25.17   950.00  15.40   32.28   117.35 

6.35    525.20  59.35   117.53  800.00  20.20   169.96  726.02 

14.08   1092.90 67.96   308.61  1040.00  28.20 434.80  966.73 

8.94    269.80  95.83   44.81   121.00  17.80   62.91   166.73 

9.62    1066.20  285.58  2528.48 13500.00 41.70 381. 1417.86 

7.41    1123.90   88.17   151.  16000.00  25.80   172.36  926.84 

8.72    267.10  65.56   29.65   63.00   21.70   36.94   100.41 ];

    （ws为八个重金属的浓度矩阵，由于数据量较大只给出一部分数据）

Points=[ 74      781     5   4

1373    731     11  4

1321    1791    28  4

0      1787    4   2

1049    2127    12  4

17    2728    6   2

2883    3617    15  4

2383    3692    7   2

2708    2295    22  4

2933    1767    7   4];

(Points为取样点x,y,海拔和功能区矩阵，由于数据量较大只给出一部分数据)

X=Points(:,1);  Y=Points(:,2); Kind=Points(:,4); 

Z=WS(:,2);   

 [XS,YS]=meshgrid(0:500:28654,0:500:18449);

ZS=griddata(X,Y,Z,XS,YS, 'v4');

figure; 

subplot(2,1,1); hold on;

mesh(XS,YS,ZS); hold on; plot3(X,Y,Z,'*');

axis([0 28654 0   18449 0 408]);

subplot(2,1,2); 

contour(XS,YS,ZS,20); hold on;

Zx=ZS(1:end-1,2:end)-ZS(1:end-1,1:end-1);

Zy=ZS(2:end,1:end-1)-ZS(1:end-1,1:end-1);

quiver(XS(1:end-1,1:end-1),YS(1:end-1,1:end-1),Zx,Zy);

2）城市功能区分布图程序：

[v,c]=voronoin([X,Y]);

figure; hold on;

axis([0 28654 0   18449]);

for i = 1:length(c) 

    patch(v(c{i},1),v(c{i},2),Kind(i));

    plot(X(i),Y(i));

end

3)重金属污染分布状况示意图程序（以下只给出Zn元素的一部分数据）：

WS=[7.84    153.80  44.31   20.56   266.00   18.20   35.38   72.35 

5.93    146.20  45.05   22.51   86.00    17.20   36.18   94.59 

4.90    439.20  29.07   .56   109.00   10.60   74.32   218.37 

6.56    223.90  40.08   25.17   950.00   15.40   32.28   117.35 

6.35    525.20  59.35   117.53  800.00   20.20   169.96  726.02 

14.08   1092.90 67.96   308.61  1040.00  28.20   434.80  966.73 

8.94    269.80  95.83   44.81   121.00   17.80   62.91   166.73 

9.62    1066.20 285.58  2528.48 13500.00 41.70   381.  1417.86 

7.41    1123.90 8.17    151.  16000.00 25.80   172.36  926.84 

8.72    267.10  65.56   29.65   63.00    21.70   36.94   100.41 ];

Points=[74     781     5   4

1373   731     11  4

1321   1791    28  4

0      1787    4   2

1049   2127    12  4

17   2728    6   2

2883   3617    15  4

2383   3692    7   2

2708   2295    22  4

2933   1767    7   4

4233   5 6   5];

X=Points(:,1); Y=Points(:,2);

Z=WS(:,8);

[XS,YS]=meshgrid(0:500:28654,0:500:18449);

ZS=griddata(X,Y,Z,XS,YS,'v4');

figure;

contour(XS,YS,ZS,15);

hold on;

plot (9328,4311,'*');

hold on;

plot (13797,9621,'+');

4）主因子空间分布等值线图程序(同样由于数据量较大，只给出一部分数据)：

A=[

 74     781     4   -0.03716    0.02207    -0.58109

1373    731     4   0.02152    -0.35653    -0.4182

1321    1791    4   -0.55108    0.25968     0.36732

0       1787    2   -0.2273    -0.17722    -0.056

1049    2127    4   -0.03695    1.39972     0.90152

17    2728    2   -0.54295    5.56314     2.06958

2883    3617    4   0.20973     0.5588     -0.54

2383    3692    2   4.40099     0.02876     11.022

2708    2295    4   -1.13543    1.02304     7.14967

2933    1767    4   0.21997     0.359      -0.73475];

x=A(:,1);y=A(:,2);z=A(:,6);

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

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

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

Zi = griddata(x,y,z,Xi,Yi);     

figure;

contour(Xi,Yi,Zi,15)

5)问题3中利用反演法求参数a,b,c,k的程序：

clear;

data=[

12855   45    18  102.86 

13797   9621    18  3760.82 

14325   8666    23  85.61 

121   9560    11  90.07 

14000   70    14  94.34 

14207   9980    14  117.87 

    ];

x0=13797;       

y0=9621;

z0=18;

M=1.60*1.0e+08;

X=((data(:,1)-x0).^2)/4; x=mean(X);

Y=((data(:,2)-y0).^2)/4; y=mean(Y);

Z=((data(:,3)-z0).^2)/4; z=mean(Z);

W=log(data(:,4));        w=mean(W);

l(1,1)=sum((X-x).^2);

l(2,2)=sum((Y-y).^2);

l(3,3)=sum((Z-z).^2);

l(1,2)=sum((X-x).*(Y-y));

l(1,3)=sum((X-x).*(Z-z));

l(2,3)=sum((Y-y).*(Z-z));

l(2,1)=l(1,2);

l(3,1)=l(1,3);

l(3,2)=l(2,3);

m(1)=sum((X-x).*(W-w));

m(2)=sum((Y-y).*(W-w));

m(3)=sum((Z-z).*(W-w));

zm=inv(l)*m';

a=(-1/zm(1)).^0.5;

b=(-1/zm(2)).^0.5;

c=(-1/zm(3)).^0.5;

ib=w-zm'*[x y z]';

k=(ib+log(a*b*c)-log(M/(2*pi^0.5)^3))^0.5

a

b

c

k

6）检验拟合的参数a,b,c,k的准确性

function f=nongdu(x,y,z)

a=534.59;

b=324.72;

c=2.;

k=2.21;

M=1.60*1.0e+08;

x0=13797;       

y0=9621;

z0=18;

f=M/(8*pi*a*b*c*(pi)^0.5)*exp(-(x-x0)^2/(4*a*a)-(y-y0)^2/(4*b*b)

-(z-z0)^2/(4*c*c)-k*k);

附录3：各个功能区八个元素的相关性系数矩阵

生活区

指标	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
As	1.000	0.381	0.238	0.531	0.293	0.605	0.450	0.017
Cd	0.381	1.000	0.349	0.499	0.397	0.283	0.802	0.346
Cr	0.238	0.349	1.000	0.376	0.150	0.527	0.416	0.412
Cu	0.531	0.499	0.376	1.000	0.198	0.434	0.502	0.238
Hg	0.293	0.397	0.150	0.198	1.000	0.211	0.340	0.242
Ni	0.605	0.283	0.527	0.434	0.211	1.000	0.300	0.334
Pb	0.450	0.802	0.416	0.502	0.340	0.300	1.000	0.328
Zn	0.017	0.346	0.412	0.238	0.242	0.334	0.328	1.000

工业区

指标	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
As	1.000	0.329	0.380	0.153	0.181	0.690	0.395	0.518
Cd	0.329	1.000	0.541	0.566	0.533	0.4	0.829	0.754
Cr	0.380	0.541	1.000	0.920	0.902	0.698	0.675	0.695
Cu	0.153	0.566	0.920	1.000	0.983	0.503	0.670	0.622
Hg	0.181	0.533	0.902	0.983	1.000	0.479	0.612	0.590
Ni	0.690	0.4	0.698	0.503	0.479	1.000	0.578	0.634
Pb	0.395	0.829	0.675	0.670	0.612	0.578	1.000	0.739
Zn	0.518	0.754	0.695	0.622	0.590	0.634	0.739	1.000

山区

指标	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
As	1.000	0.291	0.113	0.527	0.075	0.078	0.205	0.176
Cd	0.291	1.000	0.066	0.090	0.246	0.049	0.766	0.606
Cr	0.113	0.066	1.000	0.3	0.006	0.945	0.107	0.627
Cu	0.527	0.090	0.3	1.000	0.505	0.358	0.122	0.252
Hg	0.075	0.246	0.006	0.505	1.000	0.045	0.226	0.170
Ni	0.078	0.049	0.945	0.358	0.045	1.000	0.028	0.629
Pb	0.205	0.766	0.107	0.122	0.226	0.028	1.000	0.590
Zn	0.176	0.606	0.627	0.252	0.170	0.629	0.590	1.000

交通区

指标	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
As	1.000	0.121	0.139	0.092	0.004	0.228	0.060	0.188
Cd	0.121	1.000	0.373	0.424	0.211	0.351	0.615	0.294
Cr	0.139	0.373	1.000	0.4	0.012	0.869	0.428	0.395
Cu	0.092	0.424	0.4	1.000	0.032	0.886	0.506	0.432
Hg	0.004	0.211	0.012	0.032	1.000	0.040	0.266	0.118
Ni	0.228	0.351	0.869	0.886	0.040	1.000	0.396	0.503
Pb	0.060	0.615	0.428	0.506	0.266	0.396	1.000	0.482
Zn	0.188	0.294	0.395	0.432	0.118	0.503	0.482	1.000

公园绿地区

指标	As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
As	1.000	0.358	0.6	0.107	0.176	0.691	0.265	0.285
Cd	0.358	1.000	0.5	0.500	0.054	0.433	0.598	0.712
Cr	0.6	0.5	1.000	0.357	0.023	0.739	0.397	0.509
Cu	0.107	0.500	0.357	1.000	0.136	0.267	0.756	0.521
Hg	0.176	0.054	0.023	0.136	1.000	0.048	0.3	0.063
Ni	0.691	0.433	0.739	0.267	0.048	1.000	0.168	0.298
Pb	0.265	0.598	0.397	0.756	0.3	0.168	1.000	0.748
Zn	0.285	0.712	0.509	0.521	0.063	0.298	0.748	1.000