2011年全国数学建模大赛A题获奖论文

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

摘要

本文旨在对城市土壤地质环境的重金属污染状况进行分析，建立模型对金属污染物的分布特点、污染程度、传播特征以及污染源的确定进行有效的描述、评价和定位。

对于重金属空间分布问题，首先基于克里金插值法，应用Surfer 8软件对各数据点的分布情况进行模拟，得到了直观的重金属污染空间分布图形；随后，分别用内梅罗综合污染指数以及模糊评价标准和模型对城区内不同区域重金属的污染程度进行了评判。

对于金属污染的主要原因分析问题，基于因子分析法、问题一的结果和对各个金属污染物的来源分析等因素，判断出金属污染的主要原因有：工业生产、汽车尾气排放、石油加工并推测该区域是镍矿富集区。随后讨论了污染源之间的相互关系和不同金属的污染贡献率。

针对污染源位置确定问题，我们建立了两个模型：模型一以流程图的形式出现，基于污染传播的一般规律建立模型，求取污染源范围，模型作用更倾向于确定污染源的位置；模型二基于最小二乘法原理，建立了拟合二次曲面方程，在有效确定污染源的同时也反映了其传播特征，模型更加清楚，理论性也更强。

在研究城市地质环境的演变模式问题中，我们对针对污染源位置确定问题所建模型的优缺点进行了评价，同时建立了考虑了时间，地域环境和传播媒介的污染物传播模型，从而反映了地质的演变。

综上所述，本文模型的特点是从简单的模型建立起，强更准确的数学模型发展，逐步达到目标期望。

关键词：重金属污染，克里金插值最小二乘法因子分析流程图

一、问题重述

1.1问题背景

随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。评价和研究城市土壤重金属污染程度，讨论土壤中重金属的空间分布，研究城市土壤重金属污染特征、污染来源以及在环境中迁移、转化机理，并对城市环境污染治理和城市进一步的发展规划提出科学建议，不仅有利于城市生态环境良性发展，有利于人类与自然和谐，也有利于人类社会

健康和城市可持续发展[1]

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

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

1.2 目标任务

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

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

(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，分析还应收集的信息，并进一步探索怎样利用收集的信息建立模型及解决问题。

二、 模型假设

1）忽略地下矿源对污染物浓度的影响；

2）认为海拔对污染物的分布较小，故只在少数模型中讨论其作用； 3）认为题目中的采样方式是科学的，能够客观反映污染源的分布。

三、 符号说明

3.1第一问中的符号说明

i p ——污染物i 的环境污染指数 i C ——污染物i 的实测值 i S ——污染物i 的背景值

m ax (/)i i C S ——土壤污染指数的最大值 (/)i i avg

C S ——土壤污染指数的平均值

i x ——

i 区域土壤的平均值

m ax()i x ——i 区域土壤评价因素上限 min()i x ——i 区域土壤评价因素下限

3.2 问题三中模型一的符号说明

,,x y z ——采样点空间坐标值 q ——所在功能区

L

——污染源 C ——浓度

R ——基准半径

rr ——污染源范围半径

3.3 问题三中模型二的符号说明

,x y ——采样点坐标值

Z

——浓度

C ——模型系数 Q ——总误差

3.4 问题四符号说明

t ——城市历史时间统计值

()i

t Y ——第i 种元素在t 时刻的污染浓度；

（i=1,2,…,8,依次代表As 、Cd 、Cr 、Cu 、Hg 、Ni 、Pb 、Zn ）

il

V

→

——第i 种元素在液态载体中的迁移矢量；

ig

V

→

——第i 种元素在气态载体中的迁移矢量；

is

V

→

——第i 种元素在固态载体中的迁移矢量；

()ig

t ω、()il t ω、()is

t ω—— 第i 种元素在气、液、固三种状态下在t 时刻的

加权值 ()N t ——其他因素

四、 模型的建立与求解

4.1问题一

4.1.1 金属元素在该城区的空间分布

针对题目中对于地理坐标的要求，我们选用Surfer 8软件对各数据点的分布情况进行直观的图形求解。首先将数据采集进该软件，接下来对数据点进行网格化处理。在这里我们采用的是克里金（Kriging ）插值法。该插值法是一种地质统计网格化方法。它首先考虑的是空间属性在空间位置上的变异分布，确定对一个待插值点有影响的距离范围，然后用此范围内的采样点来估计待插点的属性值。由于它考虑了采样点的形状、大小及与待估计地段相互间的空间位置等集合特征以及品位的空间结构之后，为达到线性、无偏和最小估计方差的估计，而对每个采样点赋予一定的系数，最后加权平均来估计块段品位的方法。[2]同时由于克里金插值法是一种光滑插值法，相对于多元回归插值法以及最小曲率插值法等其他插值法在解决像本题这一类数据点较多的问题内插的结果可信度较高，且在空间问题上具有合理的解决方式，所以采用该插值法。最后进行制图和统计，具体方法如下：（1）利用该软件的[地图]中的[等值线图]和[线框图]功能制作出金属浓度地理分布图；（2）综合等值线图及线框图对重金属高浓度区域进行标注，分别求出其二维坐标范围；（3）从统计数据中分理处该范围内的取样点数据，得到该城区8种主要重金属元素在该城市的空间分布。

图1 各取样点的As 含量分布图

500010000150002000025000

x(m)

各取样点中As 含量的分布图

05000

10000

15000

y (m )

以上为各取样点的As含量分布图（其余元素分布图见附件1），通过分析可得该元素主要分布在如下三个区域：（1）3703

以上组合图形，一方面弥补了文中略去的统计数据，另一方面也更为直观的显示出了重金属元素高浓度范围区域。为接下来高浓度划分区域提供了直接的依据。

表1 8种重金属元素高浓度区域

8种重金属元素高浓度区域划分

x1 x2 y1 y2

As 1698 8234 2464 8973 11263 14292 1762 4602

16437 20351 8430 12004

Cd 813 18513 918 6093 16811 23347 8679 13445

4693 8982 8679 12424

7621 11501 15487 17938

Cr 2345 5272 3437 7386

Cu 0 6089 1633 10688

Hg 915 3911 1225 4425 12284 15007 1157 3471

13673 16520 7915 10314

Ni 1135 5770 1989 8970 17005 27149 7007 14969

Pb 507 6293 782 6977 12284 15620 10381 13445

Zn 541 6123 918 8884 7961 14837 1667 11062

表2 土壤8种重金属元素统计值

元素含量最小值含量最大均值标准差变异系数

As 1.61 30.13 5.68 3.02 1.88

Cd 0.04 1.62 0.3024 0.22 1.37

Cr 15.32 920.84 53.51 70 0.76

Cu 2.29 2528.48 55.02 162.92 0.34

Hg 8.57 16000 299.71 1629.54 0.18

Ni 4.27 142.5 17.26 9.94 1.74

Pb 19.68 472.48 61.74 50.06 1.23

Zn 32.86 3760.82 201.2 339.23 0.59

表2中，从标准差和变异系数可见,，各重金属含量的离散程度较大，其中

As 、Cd 、Ni 、Pb 的变异系数均较大，变异系数超过100 % ；Cr 达76% ，均达强变异程度，说明土壤中Zn 、Pb 、Cd 、Hg 受外界干扰比较明显，空间分异较大；其余元素的变异系数为27%-58% 之间，变异强度属中等，说明它们受外源影响相对轻微或更普遍更均匀。 4．1．2 污染指数模型

这里我们首先采用土壤重金属的单项污染指数评价方法，在此基础之上采用内梅罗综合污染指数法并参照土壤污染水平分极标准国家二级标准（GB15618-1995）。其中i p 为污染物i 的环境污染指数，i C 为污染物i 的实测值，

i S 为污染物i 的标准值，m ax (/)i i C S 为土壤污染指数的最大值，(/)i i avg

C S 为土壤

污染指数的平均值。

单项污染指数模型：

/i i i p C S =

内梅罗综合污染指数：

2

2

max (/)(/)2

i i i i avg

i C S C S P +=

将各个功能区的内梅罗综合污染指数计算出来的值见表3

表3 各功能区综合污染指数

区域 P(As)

P(Cd)

P(Cr) P(Cu) P(Hg)

P(Ni)

P(Pb)

P(Zn) 1 2.5642 5.8962 17.0539 13.5907 11.2695 2.1604 10.8919 29.7514 2 4.5256 6.3175 6.627 135.6196 273.0498 2.6541 10.1423

16.905 3 2.3002 2.3668 4.0525 3.814 4.2589

4.3476

2.7273

2.4719

4 6.0234 9.025

5 21.04

6 73.1892 323.3748 8.2545 4.3859 38.6209 5 2.6032 5.7798

2.4112

7.8454

27.1573 1.8898

5.3686 14.3258

根据土壤污染水平分级标准采用国家土壤环境二级标准（GB15618- 1995）。土壤污染综合污染指数可分为以下等级：（1）综合污染指数 >3 为重污染；（2）综合污染指数 2-3 为中污染；（3）综合污染指数 1-2 为轻污染；（4）综合污染指数 0.7-1 为警戒级；（5）综合污染指数 ≤0.7 为安全级。

从表中我们可以直观的看到，五个功能区的环境污染都十分的严重，最轻的也为轻度污染，交通区和工业区的污染尤其严重——所有元素都为，而居民区的铅和锌污染很严重，山区的污染程度最轻。但是这组数据评判标准有一些问题：某些污染指数如Hg 超标十分严重，超出标准几百倍且在山区也为重污染，这不禁使我们感到这套评价体系很可能并不使用于该地区的土壤环境评估。

对此，我们引入模糊评价标准和模型，并以Hg 为例说明其使数据的评价作用更符合实际。

4.1.3 土壤环境标准隶属函数模型 土壤清洁度的隶属函数:

1

m ax()()m ax()m in()0

i i

i i i x x u x x x ??

-?=?-?

?? m i n ()m i n ()m a x (

)

m a x ()i i i

i i i i

x x x x x x x ≤≤≤≥ 其中i x 为i 区域土壤的平均值，m ax()i x 为i 区域土壤评价因素上限，min()i x 为i 区域土壤评价因素下限。其中m ax()i x 和 min()i x 可由国家评价标准结合本题的背景数据共同完成。以Hg 为例，其国家标准见表4

表4土壤环境质量标准（GB15618－1995）

级别元素 一级 二级 三级 汞（mg/g ） ≤

0.15

0.30

0.50

1.0

1.5

上表中的一级标准即对应的是背景值，这与题目数据中所给的背景值相差较大，故将其更改为我们的背景值，同时，为了更加清晰的看到各区的污染程度，我们采用五级标准制，见表5

表5结合本题制定的Hg 土壤质量标准值

Hg 土壤环境质量标准值

污染等级 一级 二级 三级 四级 五级 污染严重程度 清洁 轻度污染 中度污染 严重污染 极严重污染 Hg 标准浓度ng/g 35 300 500 1000 1500 现在我们由Hg 的模糊隶属函数即可求出其污染等级，见表6 Hg 土壤环境标准隶属函数：

11,

35;300(),35300;

2650,

300.H g i x x

u x x x ≤??

-?=≤

≥??

235

,265500(),2000,

H g i x x

u x -???

-?=??

???

35300;300500;500,35.x x x x ≤<<<≥≤

()

3300

,2001000,5000,

H g i x x u x -???

-?=?????

300500;5001000;300,1000.

x x x x <≤<<<≥

4

500

,

500

1500 (),

500

0,

H g i x

x

u x

-

?

?

?

-

?

=?

?

?

?

?

5001000;

10001500;

500,1500.

x

x

x x

<≤

<<

<≥

5

1;

1000 ();

500

0;

H g i x

u x

?

?

-

?

=?

?

??

1500;

5001000;

500.

x

x

x

≥

<≤

<

表6 Hg浓度在五个区域中的模糊聚类评价等级

Hg浓度在五个区域中的模糊聚类评价等级

区域 1 2 3 4 5

污染等级一级三级一级三级一级由此表，我们看出交通区和工业区的Hg污染等级为中度污染，其他几个区的污染等级为清洁，比较符合我们的实际情况。

用同样的方法，我们求的了其他七种元素的污染等级，见表7

表7其他元素浓度在五个区域中的模糊聚类评价等级

其他元素浓度在五个区域中的模糊聚类评价等级

1 2 3 4 5

As 一级二级一级一级一级

Cd 三级三级一级三级三级

Cr 二级二级一级二级一级

Cu 二级四级一级二级二级

Ni 二级五级三级五级二级

Pb 二级三级二级二级二级

Zn 三级四级一级三级二级

与此同时，我们在附录中还提供了一种单项污染指数模型的修正模型，希望对该模型的建立有所帮助。

4.2问题二

对于问题二，意在说明重金属污染的主要原因。

4.2.1基于因子分子法判断重金属污染的主要原因

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

因子分析从变量的相关矩阵出发将一个m 维的随机向量X分解成低于m个且有代表性的公因子和一个特殊的m维向量，使其公因子数取得最佳的个数，从而使对m维随机向量的研究转化成对较少个数的公因子的研究。[3]题中土壤单点样重金属元素含量的数据特征完全符合因子分析的要求，在这

里以Hg 、Cd 、Pb 、As 、Cu 、Cr 、Ni 、Zn 八种重金属元素指标作因子分析，这样在解释各指标变化异常时可以着重讨论综合指标因子，同时为题中重金属污染成因的解释提供一定的理论依据。下面对各功能区土壤采点样重金属元素含量的数据标准化处理后，经SPSS13.0 for windows 统计软件进行因子分析，可得出以下结果。

1）给出题中表层土壤Hg 、Cd 、Pb 、As 、Cu 、Cr 、Ni 、Zn 八种重金属原始含量数据的相关系数矩阵。

表8 相关系数矩阵

Correlation Matrix

1.000.255.189.160.064.317.290.247.255 1.000.35

2.397.265.329.660.431.189.352 1.000.532.10

3.716.383.42

4.160.397.532 1.000.417.49

5.520.387.064.265.103.417 1.000.103.298.19

6.31

7.329.716.495.103 1.000.307.436.290.660.383.520.29

8.307 1.000.494.247

.431

.424

.387

.196

.436

.494

1.000

As Cd Cr Cu Hg Ni Pb Zn

Correlation

As

Cd Cr Cu Hg Ni Pb Zn

表中可见，Ni 和Cr 的相关系数最大，为0.716，相关性最好，从成因上来分析，相关性较好的元素可能在成因和来源上有一定的关联。其原因肯能是： 1.工业生产的相互关联性，镍和铬合金是工业上常用的合金。

2.考虑到表7中山区也含有很高浓度的镍，可以推断该城区富含镍矿，而对镍矿的开采和镍铬合金的冶炼是其污染的主要原因。

其次为Pb 和Cd ，相关系数为0．660，这确实出乎我们的意料，因为一般镍镉合金出现在一起比较多，但了解了Pb 和Cd 的物理性质后，我们不难发现这两种金属的一些共同点：

1.它们都广泛用于蓄电池的制造且易于以气态的方式进行传播

2.它们都是工业废水的主要组成部分[4] 。所以推断其污染原因是工业废水的排放。另外，我们从表7中可以看到Pb 的污染等级很平均，那是因为像铅在石油中含量很大，故在工业区中其污染原因是石油能源的燃烧与化工产品的生产；而汽车也以石油为动力源，故在交通区，其也有广泛的分布；而化工产品的使用，生活区的铅含量有时会很大，这在附件数据中也有反映，铅的最大浓度出现在生活区。

以下依次是Cr 和Cu ，Pb 和Cu ，相关系数分别为0.532，0.520，其它元素之间的相关性并不是很好。

2）利用相关系数矩阵求出相应的因子的特征值和累计贡献率，见表9

表9 特征值和累计贡献率

Total Variance Explained

3.5604

4.50044.500 3.56044.50044.500 2.0592

5.73325.7331.15014.37758.877 1.15014.37758.877 1.76722.08747.820.96512.06370.941.96512.06370.941 1.21715.21563.035.7689.59680.537.7689.59680.537 1.02612.82475.858.5787.22087.75

6.578

7.220

87.756

.952

11.898

87.756

.432 5.39993.156.301 3.76996.924.246

3.076

100.000

Component 12345678

Total % of Variance

Cumulative %

Total % of Variance

Cumulative %

Total % of Variance

Cumulative %

Initial Ei g envalues

Extraction Sums of Squared Loadings Rotation Sums of Squared Loadings

Extraction Method: Principal Component Analysis.

从表中可见，在累积方差为87.756％ (>85％ )的前提下，分析得5个主因子，可以看到5个主因子提供了源资料87.756％的信息，满足因子分析的原则。从上表可以看出旋转前后总的累计贡献率有发生变化，即总的信息量没有损失。而且，表中显示旋转之后，主因子1和主因子2的方差贡献率均为22％左右，主因子2到主因子5的方差贡献率的范围为11.898％到15.215％之间。这可以解释为因子1和因子2可能为本市土壤重金属污染的最重要的污染源，从而推测出金属污染的原因主要是本市重金属污染的贡献最大，因子3、因子4、因子5对本市重金属污染有重要作用。而从上面分析和下面的因子组成可知因子1和因子2主要是工业污染和交通废气污染，故更加确定了污染原因。 3）因子的相关分析

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

表10正交方差最大旋转后的矩阵

Rotated Com ponent Matrix

a

.882.209.002.014.146.864.089.019.222.196.614.362.505-.022.029.170.877.082.102.132.195.832.191.121.208.013.134.953.030.083.131.154.023.970.076.271

.285

.101

.091

.904

Cr Ni Cu Cd Pb Hg As Zn

1

2345Component

E xtraction Method: P rincipal Component Analysis. Rotation Method: Varimax w ith Kaiser Normalization.

Rotation converged in 5 iterations.

a.

基于旋转理论：变量与某一个因子的联系系数绝对值(载荷)越大，则该因子与变量关系越近。正交因子解说明：因子1为Cr 、Ni 和Cu 的组合，这说明这几种土壤重金属污染物可能是同一来源或相似来源；因子2为Cd 和Pb 的组合表明两者可能有相似的来源；因子3为Hg ；因子4为As ；因子5为Zn 。

这里，我们对Hg 的污染原因作出分析。Hg 在土壤的污染严重，污染区域大多处于某一工业区或某大型污染企业，其来源比较单一，废气和废水是其污染的主要来源。[5]

表11中国土壤重金属污染物来源[6]

来源

重金属

矿产开采、冶炼、加工排放的废气、废水和废渣

Cr Hg As Pb Ni Mo

煤和石油燃烧过程中排放的飘尘

Cr Hg As Pb 电镀工业废水

Cr Cd Ni Pb Cu Zn 塑料、电池、电子工业排放的废水

Hg Cd Pb Ni Zn

Hg 工业排放的废水

Hg 染料、化工制革工业排放的废水

Cr Cd 汽车尾气 Pb 农药、化肥

As Cu Cd

4.3问题三

4.3.1 模型综述

对于问题的第三问我们主要采取了两个模型：

一是以全体坐标为研究对象，通过假设污染源来进行比较和检索，从而确定污染源的范围。

二是以第一问中的图像和污染源范围为基础，利用已经筛选了的数据，基于最小二乘法的原理，作出其污染源的目标函数，直接求出污染源的空间坐标。 4.3.2 模型一

根据第一问所得出的8种主要重金属元素在该城区的空间分布图，以及统计得出的八种元素的浓度排序，首先在元素浓度集中分布区选择一个浓度较高点(,,,)i i i i L x y z q ，其中,,x y z 为其空间坐标，q 为其所在功能区，由于题目中所提到的采样特点，我们可以找到其相邻的两个采样点1111(,,,)i i i i L x y z q ----和

1111(,,,)i i i i L x y z q ++++，我们根据污染物扩散的基本原理，认为一般情况下：污染物

浓度与距污染源的空间距离成反比关系，反映到所选的浓度最高点上，即：

12

2

2

111()()()

i

i i i i i i nC C x x y y z z ±±±±=

-+-+-

我们在这里设这个最初距离为1i R ±。

所以我们要做的第一步是：若现实情况下1i i C C ±<，我们则保留

(,,,)i i i i L x y z q ，若不然，我们就用1i C ±中较大的点代替(,,,)

i i i i L x y z q ，并重复上

面的比较。

接着我们进行第二步，将1i R ±中的较大值赋值给R 作为基准值，其编号赋值给m 随后顺藤摸瓜，找到另一个与1111(,,,)i i i i L x y z q ±±±±相邻的采样点

1111

(,,,)m m m m L x y z p ±±±±，计算它与(,,,)i i i i L x y z q 的距离r 。

在第三步中，我们对R 与r 进行比较，若R r <：就对m C 与1m C ±做比较，若

m

C >1m C ±，则将1m ±赋值给m ，r 的值赋值给R ，继续顺藤摸瓜。若m C <1m C ±，

则将1m C ±再与i C 做比较，若i C >1m C ±，我们则考虑这可能是个意外因素，比如1m C ±处的地势较低或是顺风区等等，仍将1m ±赋值给m ，r 的值赋值给R ，继续顺藤摸瓜；若i C <1m C ±，则重复第一步中有关置换浓度最高点的方法。若R r >，则直接对1m C ±和i C 进行比较，若i C >1m C ±，则则将1m ±赋值给m ，r 的值赋值给R ，继续顺藤摸瓜，若不然则重复第一步中有关置换浓度最高点的方法。

第四步中，我们要对这个循环算法选一个终止条件，在参照了第一问图中污染物聚集区的大体范围后，我们以3i m ->，即以一个在其相邻上下三个采样点中浓度保持递减的采样点为污染源范围基准点。同时，我们以rr 为半径来划定污染源的范围。

11

2

i i R R rr +-+=

下面，我们用一个流程图来形象的描述我们的求解过程： 图3 流程图（见下页）

开始

开始区域浓度最大值Ci

区域浓度最大值Ci 1

i i c c ±>是

用其中较大者代替

用其中较大者代替否

i c 计算 和计算 和1i R +1

i R -选取R 值较大者作为新的比较基准其角标赋值于m

选取R 值较大者作为新的比较基准其角标赋值于m Yes

R

计算r 的值m next

C C >R=r

R=r NO

i next

C C >R=r R=r Yes

NO

|i-m|>4

|i-m|>4Yes

No

Yes

2

f r R R -=

输出

输出

由于要输投入的数据量较大，在这里，我们仅以As 为例作出其污染源范围，相关程序附于附件上。所得到的污染源范围是以（74，781，5）为原点，以1042m 为半径的区域范围。 4.3.3 模型二

我们建立这个模型的基本思路是：依照各元素的浓度与坐标点的关系，利用最小二乘法拟合曲面方程。同时，通过对所求的曲面方程求极值，得到的极值点一方面可以在空间上反映重金属元素的分布，另一方面可结合极值点附近的元素浓度分布，来判断污染源的位置。 1.曲面方程的规范化形式及其系数

本题中，我们采用二次曲面方程表示各浓度与坐标点的关系。二次曲面方程可表

示为[7]

201102100100,0

22

2

()i j i i

j i

j i j

Z C x C xy C y C x C y C C x y

--===+++++=

∑∑，

2.曲面方程的拟合

对于给定的一组数据点k k ,y ,z k （x ），k=1,2,…,N ，求作m 次(m=N)多项式

,0

2

()i j i

i j i

j i j

Z C

x y

--===

∑∑

使总误差

,2

00

1

()i j

i k k N

i j i

k k j m

i j

Q z C

x y ---

===?

?=

??

?

?∑

∑∑

为最小。

这里Q 可看作是关于Cij (i,j 0,1,,m )==?的函数，共有(m+1)(m+2)/2个未知量ij C ,所以上述拟合多项式的构造问题可以归结为求多元函数的极值问题，即i,j-i C 应满足

,0i j i

Q C -?=?，因而得

,00

1

()0i j i

k k

N

i j i

i j i k k k k j m

i j

z C

x y x y ----

===???

???

=???????

???

∑

∑∑

化简并交换求和符号后得，

10

0i ,j-i 1

C N k N

I J I

J j I i I i k k

k

k k j m i j k x y z

x y =??-+--+??????

===??

=??∑??

∑∑∑[8]

本题中以Cu 元素为例，在其空间分布图的基础上，缩小范围拟合其浓度与坐标的关系。

表 12 Cu元素集中分布数组

编号x(m) y(m) 海拔(m) Cu (μg/g)

13 2427 3971 2 26.57

21 2486 5999 2 29.01

7 2883 3617 15 44.81

150 **** **** 27 46.01

23 3573 6213 5 104.89

14 3526 4357 7 123.27

6 164

7 272

8 6 308.61

22 3299 6018 4 1364.85

8 2383 3692 7 2528.48

利用最小二乘法求二次曲面方程的过程如下（X 、Y单位取km）：

（1）求出坐标x、y立方、四次方对应的值（见附录）；

（2）计算z 与坐标x、y平方的乘积（见附录）；

（3）列出包含6个未知系数的方程组，利用matlab求解出的系数为：

C=1.0e+003 *

-0.7776 0.2040 0.1115 3.5721 -1.8003 0.0071

即所拟合的二次曲面方程为：

22

=+-+++-+

Z x xy y x y

1.0e003*(0.77760.20400.1115 3.5721 1.80030.0071)对于污染源，我们可以考虑对曲面方程求极大值，相当于对

22

' 1.0e003*(0.77760.20400.1115 3.5721 1.80030.0071) =+---+-

Z x xy y x y

求极小值。

也即22

=---+

min''0.77760.20400.1115 3.5721 1.8003

Z x xy y x y

（忽略了共同比例和常数项）

4.4问题四

4.4.1模型的优缺点

在此我们针对第三问中的模型进行评判。

对于模型一，其优点是：

1)层次分明，易于理解；

2）覆盖的数据全面；

3）考虑因素全面，比如考虑进了海拔的影响。

其缺点是：

1）数据量大，输入繁琐枯燥；

2）污染源表示为以采样点为中心的圆形区域，无法精确到点。

对于模型二，其优点是：

1）数学模型清晰；

2）可以将污染源确定到点，进一步缩小了污染源范围；

3）数学原理明确；

4）需要处理的数据量小。

其缺点是：

1）所取的拟合点有限，且模型不可能通过所有采样点； 2）忽略了海拔对采样点的影响。

对于这两个模型，由于受资料的限制，无法讨论污染源扩散与时间，天气，温度的关系。

4.4.2 地质环境演变模型

研究城市地质环境的演变环境首先要对该城市在不同时期内的水陆分布有初步的坐标定位以便作出比较。从土壤分类上可分为饱和土壤和非饱和土壤。饱和土壤的意思是土壤中所有孔隙全部充满水。一般情况下，土壤都是出于水分不饱和状态的。只有在淹渍条件下，才会出于水分饱和状态。所以假设城市土壤绝大部分为非饱和土壤。将土壤视为非饱和多孔介质, 描述非饱和多孔介质中多种迁移场量耦合的热质迁移机制, 与污染物在土壤中对流- 弥散迁移规律相结合, 并且考虑土壤固体骨架对重金属污染物吸附的影响因素, 建立相关的数学模型.得出有益的结论.总体上，可以将土壤中重金属迁移方式分为以下几种：1）土壤中迁移 2）水流中迁移 3）大气中沉降及降雨。以上三种形式中重金属分别是夹杂在气液固三种形式的水载体上被土壤吸附。 第一步，根据性质建立方程： （1）连续性方程

非饱和土壤中, 存在液相的蒸发或水蒸汽的冷凝, 因而在连续性方程中出

现蒸发率或冷凝率m ?

；混合气体中包含空气和水蒸汽; 水蒸汽的迁移是在整体气体迁移

g

V

→

上，叠加了蒸汽的扩散运动

V

V

→

，从而呈现

g

V

V

V

→

→

+

液相：

()()

l

l

N

l

l l

m V ρερετ

?

→+=-??

； （1）

混合气体：

()()

g

g

N

g g g

m V ρε

ρ

ετ

?

→

+=-?? ； （2）

水蒸气：

()g

V

g V

g V N m V V ρε

ρε

τ

?

?

???→→

++=??

????

?

??

； （3）

（2）动量方程

在非饱和土壤中, 液相流动受到毛细抽吸力、Darcy 阻力和重力的作用; 将孔隙中的空气和蒸汽看作为理想混合气, 不仅考虑混合气的整体迁移, 而且也考虑蒸汽相对于空气的扩散运动, 混合气体除了受到气相压差、黏性力、重力作用外, 还要克服固体骨架对它的达西阻力。

1）液相：

2

._g l

l l

l

l

l l

l

l

l g l

l

l

g

l

l

m

g g g D V g

N V N V V V

V V N V

K

K K

εεενρτ

ε?

→

→

??

→→→

→

→→→

??+-=--

-

+-

? ?????

??

-————（4） 2）气相：

2

._g

g g g

g g

g V l

g l g

g

g

g

g

g

m

g g V P

g

N

V V N V V V V

V V N V

K K

εερρ

τ

ε

?

→

→

????

??

→→→→→→→→

+++=--

-

+-

? ? ?????

??

??

————（5） （3）蒸汽扩散方程

在蒸汽随空气迁移的同时, 由于温度梯度和湿分的变化, 引起蒸汽相对于空气的扩散，因此蒸汽迁移的速度应为二者的迭加值, 即:

,,V

g

V d

TV

LV

l

V d

V

V

V

N N

T D D V

ε=+→

=--

（6）

（4） 能量方程

固、气、液三相处于局部热力学平衡态, 即有: T = T l = T g = Ts , 建立能量平衡方程如下:

()

()

()

()

,,,,()

P l l

P g g

l

g

m

l

l

l P l g P g m

l

g

l

g

c c N c N T T T V V c N c N N

N m T

T

T

V V ρρρε

ε

τ

ρρε

ελ?

→→

????

??+?

+?+?????

?????

??→

→

????

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

?

?

??

??

————（7）

在以上方程中, 下标l 、g 、v 、d 、m 分别表示液体、气体、蒸汽、空气和表观量; T 、ε、P 、V →

分别表示温度、相含量、压力和速度矢量; 量、压力和速度矢量; ,V d

V

为蒸汽相对于空气的扩散速度; 这里气体压力采用分压不分容

的原理, p c 、ρ、v 分别表示定压比热容、密度和运动黏度; τ、g 、g →

分别表示时间、重力加速度和重力矢量; v D 、l D 、l K 、g K 、m λ分别表示水蒸汽向空气的质扩散率、水液向多孔骨架的扩散率、非饱和导水率、非饱和导气率和表观导热系数; T V D 、LV D 分别为温度梯度和湿分梯度引起的蒸汽扩散运动的扩散系数。

（5）重金属污染物在土壤中的对流- 弥散-吸附方程

重金属污染物在土壤中的迁移受土壤中液态流体流动即对流、污染物在土壤中的扩散, 以及土壤中固体骨架对污染物吸附的影响. 综合以上因素, 建立以

下对流- 弥散- 吸附方程:

2

2

2

2

111s

C X

C Y

s

q

C C C C C u vc D

D

y

x x

y

ε

ετ

τ

++=

+

-

-?????????? （8）

土壤中重金属污染物在土壤固体骨架和土壤溶液中的吸附平衡方程:

1s

q

bc q

bc

=

+ （9） 将吸附平衡方程( 9)结合于方程( 8)中, 得出:

()

()()()

()2

2

2

2

////21/11C X

C Y

l

l

s s s y

y

C x

C

u C

x

v C

D

D

C q b

bC τ

εε+

--=

??

???????

?

?

???

????

??????+-+ ??????

? （10）

在以上方程中, C 为重金属污染物在土壤溶液中的浓度; C X

D

、C Y D 分别

是溶液中的重金属污染物在水平方向和竖直方向的扩散系数;

q 、s q 、s ε分别

表示重金属污染物在土壤固体骨架上的吸附量、最大吸附量和土壤固体骨架的相对含量; b 为吸附常数;

l

u

、l v 分别是土壤中液相流体在水平方向和竖直方向

的速度; x 、y 分别为水平和竖直方向的坐标。

本题中可以将土壤视为非饱和土壤, 以Philip 和deV ries 、Whitaker 、Lu ikov 3种理论为基础, 结合了污染物在土壤中对流- 弥散-吸附特性, 建立了描述重金属污染土壤的热质迁移数学方程. 方程组( 1)~ ( 7), 以及方程( 10)共有8个方程, 含有8个物理场量（T 、P 、l ε、l V 、g V 、v V 、m ?

、C ），故方程组封闭, 具有可解性。[9]

第二步，建立污染特性方程 参数说明：

t:城市历史时间统计值

()i

t Y ：第i 种元素在t 时刻的污染浓度；

（i=1,2,…,8,依次代表As 、Cd 、Cr 、Cu 、Hg 、Ni 、Pb 、Zn ） il

V

→

：第i 种元素在液态载体中的迁移矢量；

ig

V

→

：第i 种元素在气态载体中的迁移矢量；

is

V

→

：第i 种元素在固态载体中的迁移矢量；

il

V

→

、

ig

V

→

、

is

V

→

与取样点和该元素污染源的三维坐标有关；

()ig

t ω、()il t ω、()is

t ω依次为第i 种元素在在气、液、固三种状态下在t 时刻的加权值，

与该段时间城市的地表水陆范围有关。

N （t ）：其他因素，例如：重金属矿由地下开采地上导致的地表重金属含量变化；商业的发展导致工业区及其他区域的迁移。

()()()()()1

1

1

i

ig

il

is

t t t t N t ig

il

is

Y V

V

V

ωω

ω

=?

+?

+?

+→

→

→

通过均匀取样以及前三问的算法确定()i

t Y

，进而求出()ig t ω、()il t ω、()is t ω的数值，

说明重金属污染物的三种传播途径在不同时刻所占的比重的差异，最终分析出该地区随时间变化产生的地质演化。

五、参考文献

[1] Kim K Myung J H ，Ahn J S ，et a1．Heavy metal contamination in dusts and stream sediment in the Taejon area[J]．Korea J Geochemical Exploration,1998，64：409-41 9．

[2] 闫昊明. Surfer8绘图指南. 2006年-06?第二版

[3]王雄军，赖健清，基于因子分析法研究太原市土壤重金属污染的主要来源.生态环境 2008,17(2):671-676.

[4]Miguel D E, J imenezD GM, Llamas J F, et al1 The overlooked contribution of compost app lication to the trace element load in the urban soil ofMadrid ( Spain) [ J ]1 The Science of the Total Environment, 1998, 215 (1 - 2) : 113 - 1221。

[5] 王应刚, 辛晓云, 郭翠花, 太原市土壤中汞污染及成因研究[ J ]. 生态学杂志, 2003, 22 (5) ∶40—42

[6]陈怀满，土壤中化学物质的行为与环境质量[M].北京：科学出版社，2002. [7][8]霍晓程，李小平，用最小二乘法拟合曲面方程.赤峰学院学报(自然科学版)第25卷第6期，2000年6月.

[9]陈威, 杨亦霖, 张爱国, 高泽世. 非饱和土壤中重金属污染物迁移机理分析安徽大学学报2020年9月底34卷第5期 2010年

[10]宋丽叶.非饱和土壤水流问题的特征差分方法及数值模拟http://202.204.72.14/kns50/detail.aspx?QueryID=17&CurRec=60

2005-10-17

，