土壤重金属传播数学模型建立

城市土壤重金属污染分布模型建立分析及优化

论文摘要

本题目研究的对象是重金属污染在城市的土壤中的扩散与传播情况，亟待解决城市重金属污染分布问题。本组认为，需要建立的模型从总体上说，可以依据物理数学方程中给出的二阶扩散方程，通过随机模拟方程参数数值，对扩散方程中的未知参数进行拟合求解，从而得到契合度最高的扩散方程以及定性判定污染源的种类，与定量确定污染源位置。

建模过程中的核心，是我们对扩散方程的拟合进行各方面的优化（包括参数的值、个数、方程的形式等等）从而得出比较契合原始数据的重金属分布的扩散方程。通过对模型方程的分析与描述来总结出重金属在城市中传播的特点，最后，我们需要更为精确的模型，所以引入海拔Z这个参量建立优化模型，并且定性地提出其他影响因素。

首先，本组对于每一种重金属污染物，以该城市地平面为x-y平面，利用MatLab 对数据进行粗略的分析绘图（包括地理图像），从直观上对污染情况有个了解，从而定性地回答问题1和2。在初步的分析中，本组发现功能区分布相互交叉，并非严格区分区域，而且污染源不只一点。这样，扩散方程所要考虑的叠加和参数选择就比较复杂。本组开始试图引入海拔这个变量时，发现虽然海拔与污染分布也存在着一定的联系（地势、距离），但是求解过程中也就是尝试求出问题三的模型方程时，直接将海拔等因素考虑进去时会使求解方程的复杂程度大大增加。系数如果为海拔Z的函数，势必会使方程的求解非常复杂。但是如果不考虑，则数据与真实误差会增加。这是本组模型建立中的难点。

所以我们希望在建立模型中，首先尽可能少地排除次要因素来初步求出拟合方程（考虑二阶扩散），再将各种因素逐一考虑进去，建立起一个更加严谨的优化模型。通过对方程进行不断加入变量优化拟合，从而得到形式比较精炼又相对精度较高的优化方程。运用这个方程，我们很快就能对污染源位置进行确定，并且进一步的通过定量的分析加深了对该城市的金属污染的认识，得到了一些新的结论。这部分也是本模型的创新点。

建立模型后，对于数据的敏感性分析、误差分析，也是本组工作中重要的一部分。对于数据的处理，我们运用到了MatLab、Excel等软件。最后本组得到了一个相对理想的模型分析结果，以验证本模型建立的合理性。也对最终优化的模型提出了自己的解释。希望能对研究该城市的地质变化与人类活动的关系提供更有利的依据。所有计算建模工作以及软件使用，都是由本组成员共同学习完成的。

关键词：重金属，污染分布，扩散方程，数据拟合，系数优化，敏感度分析。

注意方程的解释、组合方程的解释、取点的解释

一、问题重述

通过对于某城市土壤地质环境调查所得到的数据，包括采样数据以及化学重金属元素浓度的数据，分析如下几个问题：

1.1通过8种重金属元素在该城区的空间分布，分析不同区域的污染程度。

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

1.3建立数学模型，分析传播特征。

1.4分析模型，提出优化。

二、问题分析与初步思路

2.1题目给出的数据分为三部分

●采样点的地理数据（已将该城市的地理条件进行格点划分，数据包括X轴，

Y轴以及海拔z轴数据）；

●每一个采样点所获得的重金属含量数据（As，Cd，Cr，Cu，Hg，Ni，Pb，

Zn的具体含量）；

●主要重金属元素的背景值，用于污染程度比对。

2.2问题具体分析及难点

本问题的实质，是重金属元素的一个扩散分布问题。所以我们尝试用点源持续输入一定浓度的物质的扩散模型来解释重金属污染。而如何对于某一种元素的具体分布，拟合出一个若干输入点的函数，可以更好地契合原始数据，是我们亟待解决的问题。其中包括函数方程形式的确定（2阶扩散、1阶衰减(优化模型时引入)），衰减参数以及每一个方向上扩散系数的确定，涉及到偏微分方程的解决以及分子动力学的知识。问题的解决要点在于：

●本问题与单纯的点源一次输入有一些不同，因为假设污染源给土壤在单位

时间内给予恒量的污染程度，所以需要考虑点源持续输入。

●城市里划分的区域是交叉的，以至于可能在非工业区集中的区域也会出现

少量的工业区样本点。所以当我们考虑某一重金属的扩散的时候，其扩散

原点不会仅有一个。

●2阶扩散、1阶衰减的扩散方程的系数拟合，需要一定的样本点支持。而在

同一片区的不同点源扩散区域重叠时，我们需要考虑扩散方程的叠加。

●不同区域、不同污染源、不同重金属元素所对应的扩散方程的系数都是不

相同的，所以我们至少需要拟定8种扩散方程。

●在优化模型的时候，一旦加入海拔（Z）对于扩散系数的影响（高度对于土

壤中重金属分布的影响），这样函数的建立会非常困难。

本题目给出的几个亟待解决的问题中，问题1、2可以用图像结合数据进行定

性分析，而本组侧重于建立一个与原始数据契合度相对较高的数学模型，更进一步的创新点是，我们希望能够结合海拔的因素，进行考虑地势因素的浓度分布分析。

2.3 初步思路

对于整个城市区域，首先建立一个扩散系数为常值情况下的扩散模型。以对于一种重金属为例，首先通过Matlab 对原始数据进行整理筛选。以背景值为衡量，将超出该金属背景值浓度n 倍的样本点（n 在5.2中进行了解释）为点源输入的初始点。建立方程形式如下：

2222(,,)x y C C C

E E k x y z t x y

δ???=?+?+???? 对于方程中的未知系数，本组进行格点搜索取值（取一定的步长）赋给未知量，

再利用此函数得出原始数据坐标点上的金属浓度，与原始重金属数据进行契合度的分析，横向比较所有格点搜索得到的函数形式，取到一个较为合理的拟合方程，来代表这个城市区域内的该重金属的分布情况。

最后，我们加入对于海拔Z 的分析。海拔高度的变化直接影响水流以及其他可能影响重金属传输渗透的因素，所以当考虑海拔这个变量的时候，扩散系数不再是一个常量，方程也更加复杂，但是无疑精确度将会更加好。可以优化我们前面给出的模型。

三、 模型假设

3.1 环境假设：

● 各个重金属污染的扩散过程相互独立，即一种重金属的污染浓度不会影响

另外重金属。

● 不考虑重金属除土壤内扩散的其余所有过程，如人为搬运、空气传播等。 ● 土壤结构均匀恒定，没有土壤固体结构吸附重金属，忽略气候变化，视为

理想环境。

● 海拔对于重金属扩散的作用小至忽略，即参数设为常值。（在之后的优化

模型过程中，此条不作为假设之一，将把海拔作为影响因素。）

3.2 函数假设

● 污染为点源持续输入，其扩散方程也仅考虑两阶扩散，衰减仅考虑在土壤

内的一届衰减。

● 水平方向上的每一个方向重金属扩散的系数视为相等。

● 点源输入的输入点选择时，忽略一些浓度不高的样本点，将浓度显著较高

的样本点作为输入点的位置建立函数。

● 点源输入的持续性体现在单位时间施加恒量的污染物。

四、 符号说明

x E ………… 给定坐标系X 方向扩散系数； y E

………… 给定坐标系Y 方向扩散系数； k

………… 点源单位时间持续输入浓度； t

………… 时间； (,)x y δ

…………δ函数； C

………… 浓度； ,x y

………… 坐标表示； ,s s x y

………… 点源输入点坐标；

五、 模型初步建立与部分问题解答

5.1 扩散模型的初步建立

模型建立主要依靠二阶扩散方程，目标是拟合一个与原始污染情况契合度最好的扩散方程。

土壤重金属扩散，近似地可以用二阶扩散方程（Laplace 方程）进行解释，本组也依靠此方程表示来建立优化前的数学模型。模型的建立推导过程见（6.1模型建立与分析），下面给出简单步骤： 1）以初始点源为（0,0）建立方程如下：

2222(,,)x y C C C

E E k x y z t x y

δ???=?+?+???? 其中(,)x y δ是δ函数，定义见（6.1）

2）解出通解，并如果定义(,)s s x y 为我们所找到的点源输入初始点，则得到该点的浓度表达式：

()(

)22ln +(,)s s x y x x y y k E E C x y ??--???

'=3）对于某给定的重金属元素，对x E ，y E ，k 进行格点搜索跨步取值，带入原始坐标得到一组新浓度数据，与原始数据进行契合度检验，求出契合度最好的那一组x E ，y E ，k 。

4）通过此拟合函数找出污染源、污染分布。

5.2 部分数据图像化

● 该城市地理形态图像（X-Y-海拔）

● 该城市不同采样点不同金属元素的浓度（与地理形态方向一致）

参考国家土壤重金属污染评价方法：/i i i P C S =。i S 为背景值。 其中i P 的值小于1为没有污染，小于等于2且大于1为轻度污染，小于等于3且大于2为中度污染，大于3则为重度污染。

通过Matlab 编程，借助每一个元素的背景值，筛选出每一种元素在该城市浓度超出背景值n （n 一般取1、3）倍以上的情况，并予以标示。

极端数据即峰值如何取：将n 取到相对比较大的情况下，使含量满足这个条件的区域小于等于5个。（整个过程用MatLab 进行模拟划分）

FANGCHENG DE FU BIAO

如果考虑高度与重金属扩散有关，则需要对海拔进行分区，然后在每一个分区中进行讨论计算。

>1 >3 >5

1.2.3. >1 >3 >8

1.2.3. >1 >3 >8

1.2.3. >1 >3 >50 1.2.3. >1 >3 >15 1.2.3.As

Cd

Cr

H

Cu

>1 >3 >

1.2.3. >1 >3 >20 1.2.3. >1 >3 >10 1.2.3.定性分析：

重金属污染主要集中在该城市的工业区以及交通区。这两个区域内的重金属含量经常出现极端值，由此可以断定，该两个区域中的污染源应当为主要污染源。

六、 扩散模型建立与结果分析

6.1扩散模型推导过程

设点源坐标为（0,0）根据涨落理论中的物质二阶扩散方程（Laplace ）列出方

程如下：

Ni

Zn

Pb

2222(,,)x y C C C

E E k x y z t x y

δ???=?+?+???? 其中(,)x y δ函数定义如下，满足：

0(00)1.(,)(0)2.(,)1x y x y x y x y dt δδ+∞

-∞

?≠≠?=???∞==??

?=???或、 对于齐次式：

22220x y C C C

E E t x y

???-?-?=??? 可以解得：

22141x y t

l C e

t

+-=21l C =

进而考虑非齐次：

2222(,)x y C C C

E E k x y t x y

δ???-?-?=??? （1） 若有：

2222(,)(,)x y C x y C C

E E k x y t x y

δ'???-?-?=??? （2） 则C '满足（1），则只需要解（2），接下来一系列的变量代换解决：

*,C C u v k

'

=

==

有：2*2*

22

)C C u v

δ??+=-?? 进一步解积分方程；

222

222

22*1

)(,)u v a u v a p q C dS

dudv

p q dpdq δδ+≤+≤+≤??=-

==

??

??

??

而经过计算得：

222

**

2()

ln u v a

C r C dS r π+≤???=??

所以：

*C =所以：

C '=所以得出了目标函数：

()(

)22ln +(,)s s x y x x y y k E E C x y ??--???

'=这就是我们需要拟合的一个方程部分。

如何解释(,)C x y '在一定情况下为负值的事实？ 考虑其次项

2244110

()x y E E x y e l C t t ??-- ?+

?

++?

???

=→→∞ ?

??

2

1l

C =

考虑最终的稳定态的时候，2l C 可以不考虑。故得到方程的解为：

()(

)222ln +(,)s s x y l x x y y k E E C x y pC C p p ??--???

''+=+=+ 这也解释了为什么(,)C x y '可能为负。至此，我们明确了模型需要的框架与方程形式。

6.2模型拟合方法、结果与分析 拟合方法（结合程序）：

建模拟合结果以及分析（二维浓度分布图）

拟合时取20个点源峰值作为输入点。原因：如果只取极端数据样本点作为点源输入初始点，则峰值过少，本组认为这样拟合会将其余峰值的样本点忽略。所以考虑到计算量与精确度，本组选择20个点源拟合，每个点源拟合5个变量。结果证明这样是合理的。

每一个元素的扩散方程是许多点源所确定的扩散方程的组合形式，所以每一个方程涉及到上百个变量赋值。故不在论文中给出具体方程，方程形式见（6.1）

Cr原始数据分布/拟合后分布Cu原始数据分布/拟合后分布

Hg原始数据分布/拟合后分布Ni原始数据分布/拟合后分布

Pb原始数据分布/拟合后分布Zn原始数据分布/拟合后分布

拟合分析

上述八个对比图表明，我们组所做的最优预测整体趋势和题目所给出的数据还是基本吻合的，尤其是预测组中浓度较高点的位置和个数与实际采点数据中浓度较高点的位置和个数是一致的，因为考虑到模型方程的连续性和插值画法的误差，所以在浓度高的点到浓度低的点之间存在一定的误差是可以接受的。进一步观测我们发现当实际采点数据中高浓度点相对集中，并且在非高浓度区域浓度波动不大（除了个别峰值外其他区域浓度基本不高且相对平稳）时，如最优预测图拟合情况较为理想，比如Cr、Cu的分布对比图。因为当高浓度点很多很杂乱时，该污染源的传播扩散就会变得相对复杂许多，污染源之间会交叉影响，使得用基本函数描述就会变得困难很多，再加上，如果污染源越多，那么受到不同环境的因素的影响的情况就越多，在各种情况的影响下我们所得出的最优预测就很难面面俱到，就像As的分布对比图。可以说如果能在总体趋势上与实际数据保持一致，那在相当大的程度上说明我们所得的最优传播扩散方程是有较好的可靠性的。

模型缺点

由于加入海拔Z这个参数后方程将会变得非常复杂，所以在模型建立过程中，我们仅认为海拔对于扩散系数的影响比较小。而这个假设与实际情况其实相差还是有一定距离的。所以，

七、误差分析以及敏感度分析

7.1 契合度分析的式子给出以及误差分析：

设120......a a 为某金属元素在区域中的20个点源，在i a 点(1,...,20)i =有扩散方程的5个变量,,,,x y s s E E k x y 需要拟合。所以一共有100个变量值需要拟合，每一个i a 成立一个最优拟合方程：

()(

)222ln +(,)s s x y l x x y y k E E C x y pC C p p ??--???

''+=+=+

最后整个区域的分布是满足20个点的方程组合F 。

则对于F ，带入原始数据的坐标，我们可以得到一组新浓度：

120

,......,C C ''（原始浓度为：120,......,C C ） 契合度与浓度分布中最大的值得商给出：

1

20441()

i i i i C C Fit Maximum C =??????'- ???

???

?=∑

结论：结合实际图像与上表，我们可以得知，单就峰值分布来看，我们的模型与原始数据的契

合度还是在合理范围之内的。从另外一面说，也就是模型建立是合理的。

7.2 数据敏感度分析

对于每一个元素来说，取20个峰值为点源，每个点源拟合出5个最优系数。现在给这100个系数以不同方向的10%与5%的数据扰动。然后通过相同的过程得出契合度，与最优解本身的契合度进行敏感度偏差百分比分析，得出下表：

八、模型优化及评价

本题目给出的数据只有关于空间和浓度的，这在相当大的程度上已经简化了模型的建立过程。我们知道影响重金属传播的因素可以分成三类：传播来源、传播途径还有衰变因素。

在来源方面我们通过查阅资料发现金属污染的来源相当广泛，包括大气中重金属沉降,农药、化肥和塑料薄膜使用, 污水灌溉, 污泥施肥，含重金属废弃物堆积，金属矿山酸性废水污染等。这些污染途径涉及了包括大气，水源，土壤等各方面的环境。各方面的环境变化势必会在一定程度上影响土壤中重金属传播的速率。

在传播途径中土壤本身的各种性质是决定性的。比如土壤孔性、土壤结构性、土壤耕性、土壤胶体、吸收性、酸碱性等等这些都会影响这金属传播的不同。另外在对本题目的分析后我们发现其实海拔对金属的传播也是有着紧密的联系的。最后不得不提的是水资源的分布也很大程度上影响着传播速率，因为在土壤中的金属基本是以离子的形式存在着的，这些离子有很多是易溶于水的，那么当水循环充分且迅速时传播的速率势必大大增加，反之类似。

最后对于金属的衰变情况，我们知道金属在土壤中是较难分解的，虽然微生物还有植被对部分金属有吸收和分解的功能但是效果不是很明显，所以很多情况下对于土壤的改善还是用了化学和物理的方法。

综合以上分析我们觉得在考虑实际建模时至少加入以下几个影响因子：1、应该按照土壤本身的性质对区域划分，分别测定不同性质的土壤中各种金属传播速率因子的大小2、为将城市的水资源分布设定等级，根据等级不同来设定扩散方程中系数的不同3、加入海拔这个变化因子从立体的角度来分析金属在土壤中的传播针对以上四点，本组认为：对于因子1、2也就是土壤本身性质和水资源分布来说，他们的影响具有普遍性。换句话说，不同区域的扩散方程应该根据上述两点因素的不同乘以一个变化量或者加上一个变化量，以表明不同土壤和水资源性质对相应区域的整体上的影响。对于因子3——海拔的变化，本组认为海拔的影响主要是体现在如何计算距离上，本质上距离上的计算方法并不影响扩散方程本身的形式，只是影响优化拟合时的标准。我们之前是将直线距离转换成折线距离来计算，但是还是与实际的地形情况有出入，所以在实际检测中如果能用卫星扫描来精确确定地区各点距离的话势必会大大增加优化拟合的准确性。

九、参考书目

其中书籍的表述方式为：

[编号] 作者，书名，出版地：出版社，出版年。

参考文献中期刊杂志论文的表述方式为：

[编号] 作者，论文名，杂志名，卷期号：起止页码，出版年。

参考文献中网上资源的表述方式为：

[编号] 作者，资源标题，网址，访问时间（年月日）。

[1]D. Bleecker, David 《基础偏微分方程》.

[2]复旦大学数学系《数学物理方程》人民教育出版社 1979.5

[3]Ockendon.John《应用偏微分方程》

[4]卢纹岱《SPSS for Windows 统计分析》电子工业出版社

[5]胡晓冬，董辰辉《Matlab从入门到精通》人民邮电出版社

[6]苏汝铿高等教育出版社《统计物理学》

土壤中重金属污染现状与防治方法

郑喜珅1，鲁安怀1，高翔1，赵谨2，郑德圣2（1：北京大学地质学系，北京100871；2：中国地质大学材料学院，北京100083）