数学建模A题 城市表层土壤重金属污染分析之欧阳法创编

2011高教社杯全国大学生数
学建模竞赛
承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。

如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一
项填写）：A
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：重庆交通大学
参赛队员 (打印并签名) ：1. 陈训教
2. 范雷
3. 陈芮
指导教师或指导教师组负责人(打印并签名)：胡小虎
日期：2011 年9 月12日赛区评阅编号（由赛区组委会评阅前进行编号）：
2011高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：
全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：
城市表层土壤重金属污染分析
摘要
本文针对城市表层土壤重金属污染做出了详细的分析，对于本题中所提出的问题一，我们利用MATLAB软件对所给的数值进行空间作图，然后分别作出了八种重金属元素的空间分布特征，然后，我们利用综合指数（内梅罗指数）评价的方法，对五个区域进行了综合评价，得出结果令人满意。

对于问题二，我们根据第一问和题目所给的数据进行综合分析，得出了重金属污染的主要原因来自于交通区含铅为主的大量排放，和工业区污水的大量排放等等。

对于问题三，我们通过对问题一中的八张重金属元素空间分布的图可以看出，发现大多数金属都呈中心发散性传播，同时经过分析，我们发现，如果考虑大气传播和固态传播，很难得出结论，在交通区，由于是汽车尾气造成的传播，发现重金属的传播无规律可循等，所以，我们考虑液态形式的传播，以针对地表水污染物的物理运动过程，以偏微分方程为建模基础，
通过和假设和模型参数的估计，得出了可能污染源位置，最后，我们对模型进行了稳定性检验即灵敏性分析和拟合检验，发现在参数变化在10%左右，模型的稳定性良好。

最后我们全面分析了模型的优缺点，，最后可以用MATLAB软件得出相应的结果。

为更好地研究城市地质环境的演变模式，测定污染源范围还应收集该地区的每年生活、工业等重要污染源的垃圾排放量，地下水流动方向以及每年的生物降解量，降雨量对重金属元素扩散的影响。

一但有污染证据，我们可以在该污染源附近沿地下水流动方向设定更多采样点，由此，我们可以构造一个三维公式来计算污染物质浓度的浮动就可以模拟三维空间内的重金属分布影响。

关键字：表层土壤重金属污染 MATLAB 内梅罗指数偏微分方程
稳定性检验灵敏性分析地质演变生物降解量
一、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。

对城
市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点，对城区重金属污染分析以及传播越来越有其必要性。

对于本题中所提出的问题一，我们利用MATLAB软件对所给的数值进行空间作图，分别作出了八种重金属元素的空间分布特征图，然后，我们运用综合指数（内梅罗指数）评价的方法，对五个区域进行了综合评价，得出结果令人满意。

对于问题二，我们根据第一问和题目所给的数据进行综合分析，得出了重金属污染的主要原因来自于交通区大量排放的含铅为主污染物，和工业区污水的大量排放等等。

对问题三，通过分析第一问中所给定各种元素空间分布规律，和查阅大量资料，我们发现，如果考虑大气传播和固态传播，很难得出结论，在交通区，由于是汽车尾气造成的传播，发现重金属的传播无规律可循等，所以，我们考虑液态形式的传播，以针对地表水污染物的物理运动过程，以偏微分方程为建模基础，通过和假设和模型参数的估计，得出了可能污染源位置，最后，我们对模型进行了稳定性检验即灵敏
性分析，发现在参数变化在10%左右，模型的稳定性良好。

为更好地研究城市地质环境的演变模式，测定污染源范围还应收集该地区的每年生活、工业等重要污染源的垃圾排放量，地下水流动方向以及每年的生物降解量，降雨量对重金属元素扩散的影响。

一旦有污染证据，我们可以在该污染源附近沿地下水流动方向设定更多采样点，由此，我们可以构造一个三维公式来计算污染物质浓度的浮动就可以模拟三维空间内的重金属分布影响。

二、问题分析
针对本题所提出的问题一，我们利用题给的数据运用MATLAB软件作出相应的八中重金属元素空间分布图，同时考虑到分析和评价城区累不同金属的污染度，我们根据数据应用内梅罗指数综合评价城区内不同区域的重金属污染度。

得到的结果令人满意。

针对问题二，通过对第一问和原始数据的分析，我们得出了重金属污染的主要原因还是交通区汽车含铅气体的大量排放，工业区不达标的污水大量排放等等。

对问题三，通过分析第一问中所给定各种元素空
间分布规律，和查阅大量资料，我们发现，如果考虑大气传播和固态传播，很难得出结论，在交通区，由于是汽车尾气造成的传播，发现重金属的传播无规律可循等，所以，我们考虑液态形式的传播，以针对地表水污染物的物理运动过程，以偏微分方程为建模基础，通过和假设和模型参数的估计，得出了可能污染源位置，最后，我们对模型进行了稳定性检验即灵敏性分析，发现在参数变化在10%左右，模型的稳定性良好。

证明了我们模型的正确性。

对于问题四，
三、模型假设
假设一：问题中附件给出的原始数据真实，有效。

假设二：（0～10cm）的地表是由土和沙砾所构成。

假设三：所有的土和蓄水层的性质在浸透到水层和未浸透水层都是均匀的。

假设四：稳定、均匀内的水流只发生在整个未浸透水层的垂直方向，以及在浸透水层水平（纵
向）平面中。

沿地表水流动方向。

假设五：物理过程起着重要作用，在此模型中只考虑物理过程（扩散），不考虑在其过程中的化
学反应。

假设六：所有的污染源均为点源。

四、符号说明
：单项内梅罗指标数
：的平均值
：的最大值
：内梅罗综合指数
：污染物实测值
：污染金属背景值的最小值：污染金属背景值的平均值
：污染金属背景值最大值
：水平扩散系数（）
：垂直扩散系数（）
：污染物质的浓度（）
：背景浓度（）
:污染物源处的浓度（）
：渗透系数（）
：水力传导性（）
：水力梯度
：污染物质的排放率（）
：有效多孔性
:污染物质的排放率（）
：阻滞因子
:复合参数
：污染开始的时间（）
：地下水的流速（）
：函数
：污染源的坐标
五、模型的建立与求解
5.1.1
问题一的模型建立
从附件-1中得到的数据，运用MATLAB软件，绘制了下面的三维地表图形，图形中的黑点就是取样的地点。

（z代表的是海拔）。

下面建立该城区内不同区域重金属的污染程度，为了更好地评价该城区不同区域的重金属污染程度，我们引入目前比较广泛应用的评价指标——内梅罗指数，进行综合评定，单项污染指数法能够比较直观地反映环境中各项污染指标的情况；内梅罗综合指数法不仅考虑到了所有评价因子单项污染程度的平均水
平，而且还考虑到了最大污染指数，因此能够更为科学、综合的反映评价区域内总体环境质量状况。

评价公式如下：
单项指标数：
（1）
其中：为污染物实测值，是依据附件给出的平均值，偏差，以及范围所定出的界限值。

综合污染指数（）采用内梅罗综合指数法：
（
）
2
其中：为单项指标数的最大值，为单项指标数的平均值。

给定污染综合指数等级的划分评定表，见表一：
表一：土壤污染评价分级标准
等级划分污染等级污染水平
10.7安全清洁
2警戒线尚清洁
3轻度污染土壤轻度污染
4中度污染土壤中度污染
5重度污染土壤重度污染5.1.2
问题一的模型求解
依据原始数据和公式（1）（2），利用MATLAB编
程求得八种重金属元素的不同区域的单项内梅罗指
数，见表二：
表二：八种重金属元素的不同区域的单项内梅罗指数
As (μg/g)Cd (ng/g)Cr (μg/g)Cu (μg/g)Hg (ng/g)Ni (μg/g)Pb (μg/g)Zn (μ生活区0.4180.9330670.766889 1.970.6202670.4585 1.158651 2.12工业区0.483333 2.1330140.593444 2.221543 2.28530.495250.269953 2.25山区0.269333 1.0154670.4328890.5773330.2730670.38625 1.0072560.86交通区0.380667 2.0857290.645 2.034886 2.122350.4405 1.132698 2.14公园绿
地区
0.417333 1.8702670.484889 1.00960.76660.38225 1.11958 1.60
为了能更好评出等级，我们求出各区域的内梅罗综合
指数，也即求每个区域的八种重金属元素的值的平
均值，最终评价出污染等级，其果见表二，最后得出
五大区域的污染等级见（表三）
表三：五大地区的污染等级
城区不同区
域区域内梅罗综合指
数
污染等级污染水平
生活区 1.718093轻度污染
工业区 1.9205133轻度污染
山区0.8350922警戒线
交通区 1.79966263轻度污染
公园绿地区 1.4854143轻度污染5.2
分析说明重金属污染的主要原因
5.2．1 该市表层土壤重金属含量基本状况分析
将实际测得的该市五大区域土壤样品重金属含量与自然区土壤重金属含量背景值进行比较发现，五大
区8大重金属含量平均值均高于背景值（表4）并且8大金属含量大部分样点超过自然区表层土壤均值，占总样点数的百分比最小为67.71%，最大为88.09%。

这说明表土层8大重金属均有外源物质的进入，并有了一定的积累。

表四
5.2．2 表层土壤重金属污染总体评价
根据五大区域8种重金属内梅罗指数（表二）求出各元素的单项内梅罗平均指数（见表五），再对照土壤的等级评价标准（表一）对表层土壤污染进行评价，
评价得出Cd、Cu、Hg、Zn这四种重金属指数达到了污染等级指数，污染级污染样点分别占79.94%、88.08%、67.71%、78.68%，均属于大面积轻度污染，其中Pb临近警戒线。

这就可以看出，Cd、Cu、Hg、Zn是构成污染的的主要原因，需要加强控制并采取相应措施进行降排，而Pb需要预防性控制，其余的As、Cr、Ni均较为安全。

表五：各元素的单项内梅罗平均指数
5.2．3地区对四种轻度污染元素的影响
经过对样点金属含量的详细分析，可以看出各种重金属在不同的区域的分不存在着很大的差异性。

于是针对达到污染等级的四种重金属在各个区域的单项内梅罗指数作出其分布图（图9）。

通过对各金属污染指数的分析得到这四中重金属元素的主要污染来源区为工业区和交通区。

工业区内和交通区内的采样点的单项污染指数明显高于其他地区。

工业排放和增加了重金属的含量，通过图9知道在工业区内的四种污染等级的重金属元素污染指数均超过了2且为两种以上的重金属复合污染，其四种重金属元素含量均达到了较
高的水平。

交通区域污染指数也较高，也为多种重金属复合污染，这说明该市的交通表较发达，车辆排污水较高，大量排放的尾气，和车胎摩擦所带来的重金属是造成交通区污染的主要原因。

而造成生活区污染的主要是Zn、Cu两种重金属元素，主要是生活产生的废弃物品，相比之下其Cd、Hg的含量较低。

Cd、Zn是造成公园绿地区污染的主要两种元素。

图9：四种重金属在各个区域的单项内梅罗指数图5.2．4 结果与讨论
综上分析得到该市表层土壤重金属含量的特征表现为：该是大部分地区表层土壤中的As、Cd、Cr、Cu、Hg、Ni、Pb和Zn8大重金属均有外源物质的进入。

测试样点单项污染指数Cd、Cu、Hg、Zn达到污染等级，污染级污染样点分别占79.94%、88.08%、67.71%、78.68%，均属于大面积轻度污染。

工业污染和交通污染为主要的原因。

其中工业污染造成的影响最大，是最重要的原因，其污染指数达到1.92。

工业区内大部分的土壤重金属单项污染指数和综合指数叫其他区的污染明显要高。

工厂排放的重金属污染物通过雨水河流、空气扩散到周围使得离工厂
较近的周边土壤受到污染。

该市由于交通发达，汽车尾气的排放量大，将尾气中含有的大量重金属元素直接排放到空气中，由于空气的传播速度较快这也造成了交通繁忙区周边受到了较严重大较大面积的污染。

对于生活区污染级元素为Cu、Zn两种，其根据上面分析主要来自生活废弃重金属的危险废物，比如废旧电器、电路板、光管、电池等，这些废弃物的污染是生活区的主要污染原因，应对生活垃圾进行分类和回收。

生活污水同样也是造成生活区污染的重要原因之一。

而公园绿地区由于土地裸露面积比较大引用被污染的河水对绿地进行浇灌，以及农药的喷洒极易对其土壤表层造成污染。

5.3 问题三的模型建立与求解
5.3.1问题三的模型建立
在建立模型之前，我们先说明一下，我们在这里只建立污水的物理运动过程，也就是说只建立关于液体传播的模型，因为根据，几种重金属污染物的空间分布，发现，呈气态是不大可能的，无法确定污染源，拿交通区来说，通过资料查核和数据显示，交通区所产的污染是由于汽车含铅尾气的排放造成的，这
样，就无法去确定污染源具体位置，所以，经分析和研究，我们选定污染物是以液态形式，也即随水流传播发撒。

建立了一偏微分方程为基础，通过合理的假设和参数估计，估计出污染源的可能位置。

针对污水的物理运动过程，即对流。

扩散和阻滞，我们用给出的各个采样点的位置坐标和采样点的海拔高度的数据利用对海拔的线性插值，画出该城区的基本地形，考虑到污水在水中的运动状态，水是以空间平动形式流动的，也就是说，不考虑垂直扩散，假设污水是以均匀流动且是以为污染源为坐标原点，以水流方向为X轴；设污染地点在原点的连续性污染过程，选用了液体中污染物质的二维扩散偏微分方程：
(3)
初始条件：
； (4)
边界条件：
(5)
(6)
找出在出时点源的解析解：
(7)
其中:
； (8)
是hantush函数，且有
(9)
(10)
其中常量参数：，，都可以从资料中查得到。

在计算前，我们首先按前面的假设对所有用到的参数分下类：
数据处理中污染源的坐标和时间是未知量，从而的值也就是位置的，因此，，，和都是变量。

由于我们只考虑水的二维运动，所以以每天每平方英尺加仑来衡量，据资料查得=265，地表水的速度，按达西定律由
，其中表示水力梯度，我们假定地表水的流动是一维的，这里取。

是有效的多孔性，由于题目是讨论地下10厘米的水位，所以，我们估计天。

弥散系数，该系数融合了两种形式的扩散，横向扩
散和纵向扩散，查得资料其值为25。

阻滞因子是基于污染物的特征和地表结构，由于根据资料表明他对污染物的扩散不是很大，这里我们就取为1，
我们用如下步骤来估计污染源的位置和和转移坐标：————直污染源为新的坐标原点，
————置新的轴和地下水量方向平行，————置新的轴垂直于新的轴，
我们构造一个方程来计算污染物随流体的运动，我们计算在每个采样点的浓度改变，并与数据集中的变化作比较，反复地修正污染源，的值直到满意为止，收敛准则是数据和预测值间的残差的平方和，要求极小的目标函数是：
(11)
其中是第个采样点污染源的浓度数据实测值，是第个采样点的预测值，是背景浓度值。

5.3.2问题三的模型求解
对问题三的模型的求解，我们首先考虑从每一个元素进行分析，利用该模型分别求出八种重金属元素的污染源位置，然后我们为了得到综合污染源，也就是说综合污染源是排出多种重金属的污染源，根据题
意，种种采样点不在山区的是以每平方公里进行采样
的，如果单元素污染点之间距离相差在2000m一类
的，我们就进行单元素污染点之间的合并。

下面我们就以重金属Hg元素作为分析：
通过对该城区地形图和Hg的空间分布图来综合
分析，，可以看到，重金属元素在一些区域浓度普遍
高于周围其他区域。

结合大气沉降和地下水渗透以及
流动等自然模型，推断得知，污染源大致位于重金属
浓度较高的区域或其周围。

然后我们找出分布图中亮
点区域对应的地形图区域，确定出污水的一维流动方
下图所示：
确定了X轴和大致区域然后，我们利用节点搜索
的算法（程序见附录），反复调整，直到搜索
出函数
（12）
的最小值，即可得到污染源的最终坐标位置。

Hg
元素的污染源坐标为：（2509,2993），
（13987,2875），（14974,8905）。

反复运用此模型最
终解得八种单元素污染源坐标结果如下表：
Zn Pb Ni Hg Cu Cr Cd As
140802010330025092400394531331805
90335210139875520
4507875028757295
3300468514974
623047508905
为了得到多元素污染源，依据上表，将单元素污
染点之间距离相差在2000m一类的，我们就进行单
元素污染点之间的合并得到最终结果，分别是
（19987,15389），（15789,10027），
（12584,8257），（9375,6831），（5922,3781），
（4723,2573），（2699,5892）七个多元素污染源点。

5.4.1应增加搜集的信息
为更好地研究城市地质环境的演变模式，测定污染源
范围还应收集该地区的每年生活、工业等重要污染源
的垃圾排放量，地下水流动方向以及每年的生物降解
量，降雨量对重金属元素扩散的影响。

一但有污染证
据，我们可以在该污染源附近沿地下水流动方向设定
更多采样点，由此，我们可以构造一个三维公式来计
算污染物质浓度的浮动，如下：
（13）
因为渗透是一个连续过程，我们规定污染源的作用如
同一个阶梯函数（连续地）并满足下列边界条件：
（14）
对于在时刻的瞬时点源，该方程有一个形为
（15）
的解析式，其中
（16）
最后我们再画出一个中心在近似点源半径为100米得圆内，从地表取了一些土样并分析其重金属成分求得最大之，因此，我们就能精确得识别污染源位置。

六、模型三的检验
为了验证我们模型的正确性和稳定性，我们做出了模型的敏感性分析，在任然以每平方公里为取样单位的话，只要位置波动值能控制在一千以内，那就说明我们的模型的稳定性是比较好的。

我们分别改变模型中常数，，，和的值，以模型中的值为基准上下波动10%，并计算相应的污染源，（程序代码）和模型中一样，见下表：
变动的参数参数较小10%情
况下的位置变化参数不变的位置
变化
参数增加10%情
况下的位置变化
300～6300400～700
200～5300180～460
400～8000350～460
200～6000330～650从上表中可以看出，数值波动都在1000一下，也也就是说，当参数变化是，对于单元素污染源的位置波动比较大，但是动控制在1000以内，也就是说，没有超过其极限值。

所以验证了我们都模型的稳定性。

我们通过对单元素污染源点坐标周围的点用模型公式中进行计算，得出相应的模拟浓度，并进行Excel软件进行拟合，相应的点在表中对应相应的浓度，得到如下图形：
通过图可以看出，峰值相差比较大，拟合效果不是很好，但他任然展示了相似的趋势，并且有很好的相似性。

证明了我们的模型正确性。

七、问题四模型优缺点的分析及优化
模型优缺点的分析
内梅罗综合指数模型
对用于污染程度评价的模型一所使用的内梅罗指数法模型进行分析
优点：数学过程简便。

物理概念清晰，评价方式简单便于决策。

对数据的处理考虑到了各个
散点数据间的联系，加入权重进行综合
内梅罗指数排名。

缺点：其描述的环境质量是非连续的，分级标准建立在二值逻辑基础上，它的截然性和非连续性造
成了相差很小的污染指数强度间可能会出于两
种不同的等级。

污染源定位模型
优点：
●模型有很好的实践性，而所给的算法几乎没
有时间的复杂性。

对于所给问题的数据规
模，我们采用格点搜索法求最优解。

●模型得出了数据与计算值鱼很好的一致性，
它是快速、有效和稳定的。

●至于对数据的简化计算，准确性并没有降
低。

缺点：
●如果考虑的区域比较大，就会有一些误差。

●为了降低计算的复杂性，我们间滑落影响污
染扩散的纵向地下水流这会影响到结果的精
确度，还有在地形方面只考虑了污染源的水
平定位没有考虑其海拔高度的定位。

八、模型的推广与改进
8.1.1模型的推广：
由于我们的模型具有一定实用性和稳定性，所以可以将该模型推广到具有稳定地形的关于水污染的模型中，我们充分考虑到在模型中遇到的一些问题，设定参数在可变的情况下具有一定的稳定性。

所以也可以应用到一些地形较稳定的地形中，也能达到良好的效果。

8.1.2模型的改进：
我们所建立的模型是将两种弥散系数a近似估计为一实数25ft，这是在外界条件相同的情况下的一种假定，而实际中不同地方土壤的PH值是不同的，特别是污染区。

水-土壤系统pH升高能明显地降低各元素在土壤中的吸附,促进其在土壤中的迁移,且吸附常数(Kf)与土壤有机质含量、粘土含量呈正相关,而与土壤pH呈负相关.pH值高金属元素在土壤中淋溶贡献较大,且淋溶量随雨量的增大而增大.同时在土壤中的淋溶与土壤性质密切相关,有机质含量和粘粒含量较高的土壤对其的持留能力较强，由此我们可对模型进行优化。

PH对传播的影响主要可以体现在阻滞力系数以及渗透力系数上，对此我们假设
（17）
并且令
(18)
带入模型得到优化后的模型如下：
(19)
在针对模型进行取样采点进行统计计算，即可得到，的值，获得后就是完善的模型.
九、参考文献
【1】孙树瑜，曾爱武．王树楹。

等．规整填料塔中精馏过程的三维模拟lⅡ)：模型的验证及液相
分布和混台行为对精馏过程的影响J] 化工学报
1998，49(5)：503—565
【2】夏星辉,陈静生.土壤重金属污染治理方法研究进展[J].环境科学,1997，18(3)72-76
【3】郑州市郊区农业区划办公室．郑州市郊区土壤普查办公室．郑州市郊区壤R1．1986，39．【4】魏秀琴．郑州市东北郊污水灌区重金属元素对。