2011全国一等奖论文

2011高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

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

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 20111431 所属学校（请填写完整的全名）： 江苏大学 参赛队员 (打印并签名) ：1. 武尧晰 2. 陈明明 3. 姚佳丽 指导教师或指导教师组负责人 (打印并签名)： 教练组

日期： 2011 年 9 月 12日

赛区评阅编号（由赛区组委会评阅前进行编号）：

2011高教社杯全国大学生数学建模竞赛

编 号 专 用 页

评 阅 人 评 分 备 注 赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

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

摘要

本文讨论了城市表层土壤重金属污染的问题，通过运用空间插值法、二维稳态污染物质浓度场基本模型，对重金属元素对城区土壤的污染进行了分析和评价。

在问题一中，首先利用反距离加权空间插值法（IDW）对采样点进行插值，得到了8种重金属元素在整个城区内的空间分布图。然后采用了Muller提出的地积累指数分级来定量评估区域内金属的污染程度，据此评价方法分析了在每一个区域内8种元素的污染程度，其中生活区内的元素As、Cr、Hg、Pb属于轻度污染，元素Cd、Zn属于轻度—中等污染，元素Cu属于中等污染。同时，比较五个区域的综合地积累指数，得出污染程度从高到低的区域依次是：工业区、交通区、生活区、公园绿地区、山区。

在问题二中，通过对各采样点8种元素的地积累指数进行统计分析，得到每种元素污染程度最大的区域，然后在该区域内分析该元素产生污染的主要原因。结论为：元素As、Hg、Ni、Pb污染的主要原因来自工业区，元素Cd、Cu、Zn来自工业区与交通区，元素Cr来自生活区内。

在问题三中，通过对重金属元素传播特性的分析，得出二维稳态污染物浓度扩散模型。而求取污染源位置是一个典型的反问题，即要通过对已知采样点数据进行拟合确定出污染源的位置。由于一种元素在整个城区内的污染源可能不止一个，而我们建立了一个适用于单个污染源的浓度扩散模型，因此要先确定几个污染源候选区域，然后用这些区域内采样点的值来拟合污染物的传播方程，从而反推出污染源的位置，其中As元素的污染源位置有：（4717,6788）、（12630,2938）、（18510,10200）、（9781,16480），（单位：米）。

在问题四中，首先评价了问题三中模型的优缺点，针对模型的缺点，为更好地研究城市地质环境的演变模式，还需要在土壤的不同深度进行采样，了解土壤中重金属的浓度随时间的变化规律，以及收集土壤对重金属元素的降解能力等信息。 关键字：城市土壤、重金属污染分析、污染物的传播、反问题分析

1

1、问题重述

1.1问题背景

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

1.2对于城市表层土壤的采样及归类

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

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

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

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

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

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

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

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

2

2、模型的假设

（1）假设各采样点的重金属污染情况可以代表其所在网格的污染情况； （2）假定题目附录中给定的数据真实可靠，具有较好的代表性； （3）假设重金属污染源释放污染物的量是稳定的； （4）假设土壤对重金属无分解能力，且不考虑土壤的迁移； （5）假设重金属元素在土壤表层各个方向扩散速度均相等。

3、符号说明

Igeo：表示地积累指数；

Cn ：表示某元素实测质量分数； Bn ：表示某元素背景质量分数；

Z*(x0)：表示x0处的估计值；

i： 表示第i个样点对估计点的贡献权重； C：表示污染物在环境介质中的浓度； Cil：第i种元素在第j个采样点的浓度； f：表示污染物的推流迁移通量； u：表示环境介质的流速分量； Im：表示污染物分子的扩散通量；

Dm：表示分子扩散系数，是指单位浓度梯度下物质在单位时间内通过单位面积的量；

pi：表示某个功能区内，第i种重金属元素地积累指数大于1的采样点数占本区域内所

有采样点数的频率；

i：表示第i种重金属元素背景值的标准偏差； Ui：表示第i种重金属元素背景值的平均值；

ai：表示第i种重金属元素背景值的所取范围的下限； bi：表示第i种重金属元素背景值的所取范围的上限。

4、问题分析

人类活动对城市地质环境的影响越来越大，而土壤中重金属的污染与人类的生活息息相关，所以了解重金属的在城市中的分布，找到污染原因和污染源，对重金属污染的治理至关重要。 4.1问题一的分析

问题一要求给出重金属元素在该城区的空间分布，一个城区数百个采集点显然不能准确得到各元素的空间分布，因此我们考虑用反距离加权法[2]进行空间插值，将离散的测量数据转换为连续的分布曲面，从而得到各个重金属元素在该城区内的分布图；评价不同区域重金属的污染程度时，需要给出一个评价指标，我们采用Muller提出的地积累指数[3]作为指标来定量评价土壤中重金属的污染程度。

3

4.2问题二的分析

在问题一中数据处理的基础上，通过对各采样点的地积累指数进行统计分析，对于某一种重金属元素，可以得到其在每个区域内地积累指数大于1的采样点数占所在区域总采样点数的比例，比较所得的这五个比例，值越大则说明该区域对这种元素的污染贡献越大，从而得出重金属元素污染的主要原因来自哪个区域，进而根据该区域的特点来分析重金属污染的具体原因。 4.3问题三的分析

本问题要求根据各采样点的污染物浓度，综合分析采集到的数据及重金属元素的传播特性来反推出污染源的位置。这是一个典型的反问题。而反问题则是根据事物的演化结果，由可观测的现象来探求事物的内部规律或所受的外部影响，由表及里，索隐探秘，起着倒果求因的作用。反问题区别于传统正问题的主要特点是反问题属于病态问题，即无法同时满足求解过程中解的存在性、单一性、以及数值稳定性等三个要素。因此在这种情况下，将无法得到污染源的精确解。于是先对污染源产生的重金属元素的传播特性进行分析，得出重金属元素污染物的传播特性方程，再通过对已知采样点数据来拟合从而反推出污染源的位置。

5、模型的建立与求解

5.1 问题一的模型建立与求解

5.1.1分析8种重金属元素在该城区的空间分布

反距离加权法[2]是最常用的加权平均法之一，它考虑了空间数据值与距离反相关的空间特性，其公式为：

Z(x0)iZ(xi) （1）

*i1N其中Z*(x0)为x0处的估计值，N为用于差值的临近点的个数，Z*(xi)为样点xi处的实测值，i为第i各样点对估计点的贡献权重，其计算公式为：

dpiNi0 （2）

(di0p)i1其中

i1Ni1，di0 为待估点与样点间的距离，p为距离的幂，以最小平均误差最为选

择标准，它控制随着距离增大权系数的变化，此外p也决定了内插结果的平滑效果。 在MATLAB中运用反距离加权法进行空间插值运算[4]，作出8种元素在城区内的分布图，并用散点表示出采样点所在的功能区，8种元素的空间分布图如下图所示：

4

图1 元素As的空间分布图

图2 元素Cd的空间分布图

图3 元素Cr的空间分布图

5

图4 元素Cu的空间分布图

图5 元素Hg的空间分布图

图6 元素Ni的空间分布图

6

图7 元素Pb的空间分布图

图8 元素Zn的空间分布图

在图1-8中，横纵坐标表示采样点的位置，颜色的深浅代表重金属元素污染程度的高低，颜色深的区域说明该区域内重金属元素的污染程度高，并结合散点图来表示不同功能区，其中不同的散点符号代表不同的功能区内的采样点，散点符号所代表的具体功能区见图上的标注。

综合比较8种元素的空间分布图，可以大致看出在整个城区内：元素As、元素Cr的深色面积比其他元素大得多，说明这两种元素在整个城区内污染比较严重；而元素Cu、元素Hg、元素Pb、元素Zn次之；元素Ni最小。在此基础上可以定性分析出整个城区内的各元素的污染程度，不同功能区的重金属污染程度定量分析将在5.1.2中详细讨论。

5.1.2分析不同功能区重金属的污染程度

Muller曾提出用地积累指数[3]定量评价沉积物中的重金属污染程度，并规定了相应的污染程度级别划分标准（见表1）。这种评价方法也可用来评价土壤中重金属的污染程度及其分级情况。

7

地积累指数的计算公式为：

Igeolog2(Cn/1.5Bn) (3) 其中Cn为元素实测质量分数，mgkg1，Bn为该元素背景质量分数，mgkg1

表1 Muller地积累指数分级

地积累指数(Igeo) 分级 污染程度

5Igeo10 6 极严重污染

4Igeo5 5 强—极严重污染

3Igeo4 4 强污染 2Igeo3 3 中等—强污染

1Igeo2 2 中等污染 0Igeo1 1 轻度—中等污染 Igeo0 0 无污染

采用上述评价方法来评价该城区不同功能区重金属的污染程度，并从两个方面来分别考虑：一是分别评价每个功能区内8种重金属元素各自的污染程度；二是在一个功能区内综合考虑8种元素的污染程度得出某个功能区整体重金属的综合污染程度，然后在五个功能区之间比较综合污染程度。

根据假设各采样点的污染情况可以代表其所在网格的污染情况，那么就可以利用Muller地积累指数级数对单个采样点进行评价，通过对某功能区内所有采样点综合评价就可以得到该功能区的重金属污染程度。

（1） 评价各个功能区内各重金属元素的污染程度

由于Muller公式：Igeolog2(Cn/1.5Bn)中要用到元素的背景值，而本题中却给定了元素的背景值的均值、标准偏差、范围，而这些背景值的区间范围是由大量的测量结果得出的，因此我们在这里假设背景值服从正太分布。其概率公式为：

21e2i （4） 2i式中：Pi(x)为第i个元素的背景值取x时的概率，x(ai,bi)。

针对第i种元素在第j个采样点都有一个地积累指数，而地积累指数即为该采样点的浓度与每一个可能的背景值的地积累指数的加权平均值(权重为该背景值出现的概率)。

即，此问题的模型：

biPi(x)Igeoijlog2(Ci,j/(1.5*x))*dx （5） aiq(xUi)2Pi(x)

其中qPi(x)dx是归一化因子，即Pi(x)/qdx1。

aiaibibi而要计算某元素的地积累指数，即为每个元素采样点的浓度代入公式（5）计算出

城区内319个采样点8种重金属元素的地积累指数。在不同功能区内，我们发现大多数采样点的地积累指数分布在小于0（对应的指数级数为0，无污染程度）、0~1（对应的

8

级数为1，轻度—中等污染）的范围内，较少数的采样点指数大于1(对应级数大于等于2，中等及以上程度污染)，所以我们在Muller地积累指数分级的基础上进行简化，如下表所示：

表2 简化后的Muller地积累指数分级

地积累指数(Igeo) 分级 污染程度

1Igeo 2 中等及以上污染 0Igeo1 1 轻度—中等污染

Igeo0 0 无污染

在本模型中采用上表作为评价标准标准。

用Matlab按照简化后的Muller地积累指数级数以上得到的地积累指数按照不同级数进行采样点个数的统计，并作出柱状图[1]如下所示

图9 生活区8种元素的地积累指数分布

图10 工业区8种元素的地积累指数分布

图11 山区8种元素的地积累指数分布

图12 交通区8种元素的地积累指数分布

9

图13 公园绿地区8种元素的地积累指数分布

图9-13中，纵轴表示落在一个地积累指数段的采样点个数，横轴表示8种重金属元素。

根据图9-13的柱状图所示就可以得出每个功能区中不同元素的污染程度，具体分析为：从图9中生活区内44个采样点8种元素的地积累指数可以看出，对于元素As：有14个采样点无污染，占总采样点一定比例；28个采样点属于轻度—中等污染，占总采样点的比例较大；只有2个采样点属于中等及以上污染程度，占总采样点个数的比例很小，综合考虑生活区中元素As应该属于轻度污染。对于元素Cd：有12个采样点属于无污染程度，24个采样点属于轻度—中等污染程度，8个采样点属于中等及以上污染程度，轻度—中等污染的采样点数比例较大，并且，无污染与中度及以上污染综合考虑对于元素Cd生活区应属于轻度—中等污染程度；同理可知生活区内对于元素Cr轻度污染、元素Cu中度污染、元素Hg轻度污染、元素Ni轻度污染、元素Pb轻度污染、元素Zn轻度—中等污染。

用相同的方法分析其余四个区域内的污染情况，得到的结论如表3所示：

生活区 工业区 山区 交通区 公园绿地区 As 轻度污染 轻度污染 无污染 轻度污染 表3 各区各元素的污染情况 Cd Cr 轻度—中等污染 中度污染 无污染 轻度—中等污染 轻度—中等污染 轻度污染 轻度污染 轻微污染 轻度污染 轻度污染 Cu 中等污染 中度污染 轻微污染 中等及以上程度 轻度—中等污染 轻度—中等污染 生活区 工业区 山区 交通区 公园绿地区

Hg 轻度污染 中度污染 轻微污染 轻度—中等污染 轻度—中等污染 Ni 轻度污染 轻度污染 无污染 轻度污染 轻度污染 10

Pb 轻度污染 轻度—中等污染 无污染 轻度—中等污染 轻度污染 Zn 轻度—中等污染 轻度—中等污染 无污染 轻度污染 轻度污染

（2）在五个功能区之间比较其整体的污染程度

在每一个功能区内，计算出第i种元素的平均地积累指数Igeoi，其值如下表所示：

表4 在各区域内各金属元素的地积累指数平均值 生活区 工业区 山区 交通区 绿化公园区 As 0.1645 0.2500 -0.4956 -0.0185 0.1650 Cd 0.3697 0.8291 -0.4950 0.6275 0.2019 Cr 0.1987 0.0211 -0.3945 0.0783 -0.1084 Cu 0.9215 1.4219 -0.3502 1.1194 0.4161 Hg 0.2296 1.3468 -0.5681 0.4933 0.2145 Ni -0.0194 0.0270 -0.3991 -0.1260 -0.2656 Pb 0.2602 0.6964 -0.4341 0.3038 0.1405 Zn 0.4922 0.8813 -0.5686 0.6763 0.0728

取8种元素的平均地积累指数Igeoi的平均值就可以得出得出第j类功能区的平均地积累指数：Ijgeo18Igeoi，其值如下表所示： 8i1

生活区 0.327132 表5 各功能区的平均地积累指数 工业区 山区 交通区 0.684191 -0.46315 0.394281 绿地公园区 0.104595 从表5的数据可以分析出，在该城区内污染程度由高到低的区域依次是：工业区、交通区、生活区、绿地公园区、山区。 5.2问题二的模型建立与求解

本模型要分析重金属污染的主要原因，而不同重金属产生污染的原因是不相同的，所以在模型中，我们对8种重金属分别进行了讨论。经过分析可知，地积累指数Igeo越大的采样点对重金属污染的贡献就越大，根据Muller地积累指数分级，可知当某元素的地积累指数Igeo1时，污染程度为中等及其中等以上水平，所以可以认为大于1的采样点都有可能是污染产生的主要原因所在地。在某个功能区内，对于第i个元素，我们把它的地积累指数大于1的采样点数占本区域内所有采样点的概率pi作为污染影响程度的评价指标，并把pi最大的功能区作为元素污染源所在地，然后在此功能区内分析产生污染的原因。

8个元素在各个功能区内地积累指数大于1的采样点占本区域内所有采样点的概率pi如下表所示

11

表6 各区内地积累指数大于1的的采样点占本区域内所有采样点的概率 生活区 工业区 山区 交通区 绿化公园区 As Cd Cr 0.1136 0.0556 0.0303 0.0725 0.0286 Cu 0.4091 0.5278 0.0303 0.5145 0.1143 Hg 0.2727 0.5000 0.0152 0.3188 0.2857 Ni 0 0.0556 0.0303 0.0145 0 Pb 0.1591 0.2778 0.0303 0.1884 0.1714 Zn 0.2955 0.3889 0.0152 0.3333 0.1143 0.0455 0.1818 0.1389 0.3889 0.0152 0.0303 0.0217 0.3261 0.0571 0.2000 分析表4中的数据可以得出，在工业区中元素As、元素Hg、元素Ni、元素Pb的地积累指数大于1的采样点占总采样点的概率比其他区域的该项比例值要大，所以这三种重金属污染的主要原因应该在工业区内进行分析；在工业区与交通区元素Cd、元素Cu、元素Zn的地积累指数大于1的采样点占总采样点的概率值相当，并且比其他区域的该项比例值要大，所以这三种重金属污染的主要原因要在这两个区域内进行分析；在生活区内元素Cr的地积累指数大于1的采样点占总采样点的概率比其他区域的该项比例值要大，所以此元素的污染的主要原因应该在生活区内进行分析。

工业区产生的重金属污染可能来自于工厂生产产生的废弃、废渣、废水的乱排乱放，以及矿山的开采和金属冶炼等工业生产活动；交通区产生的重金属污染可能来自含铅汽油的燃烧，汽车轮胎磨损产生的含锌粉尘等；生活区的重金属污染主要来自生活垃圾的随意丢弃、堆放，以及燃煤不完全燃烧。 5.3问题三的模型建立与求解

5.3.1重金属污染物的传播特征与模型建立

（1）概念引进[5]

①介质的迁移运动——推流迁移

推流迁移是指污染物在气流或水流作用下产生的空间位置上的转移。单纯的推移作用不能降低污染物的质量和浓度，数学抽象用迁移通量（单位时间内通过单位面积的量）来描述污染物的推流迁移：

fxuxCfyuyC （6） fuCzz式中：fx、fy、fz——在x、y、z方向上的污染物推流迁移通量，ML2T1； ux、uy、uz——环境介质在x、y、z方向上的流速分量，LT1；

C——污染物在环境介质中的浓度，ML3。

②分子的扩散

分子扩散是指由分子的随机热运动引起的质点分散现象。它存在于污染物的所有运动过程中。分子扩散过程符合菲克（Fick）第一定律，其质量通量与扩散物质的浓度梯度成正比，可由下式表示：

12

CIDmm,xxC （7） Im,yDmyIDCm,zmz式中：Im,x、Im,y、Im,z——x、y、z方向上的污染物分子扩散通量，ML2T1；

Dm——分子扩散系数，是指单位浓度梯度下物质在单位时间内通过单位面积的量，L2T1。

由于扩散都是由浓度高处向浓度低处的质点迁移，故都应取负梯度方向。 （2）建立二维稳态模型

假设重金属污染物的污染源在原点，且不考虑由土壤的自我降解能力等因素而产生的污染物的减少，

由多重积分的意义可知，通过闭合曲面S从时刻t到时刻tt流入的污染物质量为：

M1tttS(DxCCCcosDycosDzcos)dSdt (8) xyz其中cos,cos,cos为S的外法向余弦。

由高斯定理可知：

tt2C2C2CM1(Dx2Dy2Dz2)dxdydzdt （9)

txyz由于重力等因素的作用，重金属元素会有一定的移动速度，污染物从时刻t到时刻tt流出S的污染物质量为：

M2ttt[uC(x,y,z,t)cosuC(x,y,z,t)cosuC(x,y,z,t)cos]dSdt (10)

Sxyz同理，有高斯定理可知：

M2ttt(uxCCCuyuz)dxdydzdt （11） xyz由质量守恒定理得：

t时间段内滞留在土壤中的重金属污染物的质量为：MM1M2所以重金属污染的四维数学模型为：

C2C2C2CCCC （12） Dx2Dy2Dz2uxuyuztxyzxyz假设重金属污染源在短时间内释放污染物的量是稳定的，且当这个过程时间足够长时，土壤中重金属的浓度可以看做将不随时间的变化而变化，只随着位置的不同而发生

C0，据此我们可以建立一个二维稳态模型，其表达式改变。在点源稳定的基础上，t为：

2C2C2CCCCDx2Dy2Dz2uxuyuz0 （13） xyzxyz对于z轴方向，其由重力等因素引起的移动远远大于其扩散作用，而且重力等因素主要集中z轴方向，所以其在x轴、y轴方向的移动可以忽略，即Dzuxuy0。故（13）可以简化为：

13

2C2CCDx2Dy2uz0 （14） xyz在污染源处，污染物的质量为M0，即初始条件为：

C(x,y)M0(x)(y) （15）

重金属元素污染物在土壤中的二维扩散是无边界的，即：x,y 边界条件为：

limC(x,y)0,limC(x,y)0xy(x)

1,x0,0,x0(y)1,y00,y0 （16）

求解微分方程（14）并代入初始条件和边界条件，得到解析解为：

ux2y2M0 （17） C(x,y,z)expz4*zDxDy4zDxDy由于只在土壤表层取样，并且只有横向扩散信息，无法了解到纵向的扩散情况。所以在模型中不考虑重金属污染物的纵向扩散。即式（17）中各点的浓度与z无关，，所以可以认为 uz、z的值对C的值无影响，可以认为是常数。假设重金属污染物在扩散平面内各个方向的速度均相等，即有：DxDy。对（17）式中常数项进行简化，令：

4zDx4zDy，D'所以，式（17）可以化简为： auu4*zDxDyzzM0x2y2C(x,y)a*expD'D' （18）

而以上的讨论是在污染源在坐标原点的前提下进行的。只要通过坐标转换，将污染源的坐标平移到(X,Y)。即得到本问题的模型：

(xX)2(yY)2C(x,y)a*exp （19） 'D'D便可由已知采样点的浓度反推出污染源的位置。而要求取污染源位置是一个典型的反问题，即要通过对已知采样点数据进行拟合确定出污染源的位置。由于一种元素在整个城区内的污染源可能不止一个，而我们建立了一个适用于单个污染源的浓度扩散模型，因此要先确定几个污染源候选区域，然后用这些区域内采样点的值来拟合污染物的传播方程，从而反推出污染源的位置。

5.3.2 模型求解

由于不同种元素的污染可能是由不同的功能区引起的，将8中元素放在一起考虑情况过于复杂，因此将这8种元素分开考虑。

以As元素为例，以其在整个城区内所有采样点的位置为x轴、y轴，以其浓度值为z轴，利用反距离加权法[2]进行空间插值拟合出这种元素三维曲面图及二维等高线分布图，如下图所示：

14

图14 As元素三维曲面城区分布图 图15 As元素二维等高线城区分布图

在分布图中，以颜色深浅度表示As元素对其所在区域的污染程度，颜色越深即表明污染越严重，并且我们规定黑色标志污染程度最严重。因此在图中可以根据颜色最深区域大致估计出局部污染严重点的坐标范围，也即污染源点的可能所在域。通过MATLAB编程在上述坐标范围内进行搜索寻找最大值，即局部污染最严重的平面坐标点，而由于图中黑色区域可能不止一处，即污染源可能不止一个，针对每一可能的污染源，结合分布图上污染较严重的区域范围给定一个搜索半径，寻找落在此区域的采样点的值。

由于在作图的数据是经过插值后的数据，在此中寻找到得局部最严重的污染点的坐标可能与在此小区域的采样点的局部最大值有差异，这样便会由于插值而带来一定得误差。又由于拟合后的局部最大值与原数据的最大值空间距离相差不会太远。因此只要控制好上一步的搜索不长，变可以得到与数据的最大值。再以此最大值出发，给出一个搜索半径，寻找落在真实最大值附近采样点。

由这些离散点坐标位置处的As元素浓度，结合重金属污染物的传播特征模型在MATLAB中通过sftool工具箱拟合出它的传播特征表达式、三维曲面图像及污染源的地理坐标。其拟合效果图如下图所示：

图16 污染源拟合图

图中各个离散点表示落在污染源候选点周围的采样点，曲面中的最高点表示污染源位置，而曲面分布表示此污染源的传播特性图。针对每一种元素，重复以上的算法，得出各元素污染源的位置如下表所示：

15

元素 表7各元素污染源的位置 污染源位置 元素 X（m） Y（m） 4717 12630 18510 9781 3695 4601 8381 17570 19000 9781 3695 4610 6788 2938 10200 16480 7760 4792 1154 3151 11750 16480 7760 4796 Cu Hg Ni Pb Zn 污染源位置 X（m） 2086 4500 1965 12640 18510 3809 18570 2099 4780 3471 10750 12600 Y（m） 2845 7000 2887 3001 10210 7762 10120 2838 4770 6000 4143 9952 As Cd Cr

5.4问题四的模型建立 5.4.1模型的评价

在问题三中运用了二维稳态污染物质浓度场基本模型，结合反距离空间插值法有效地拟合出了8种重金属污染物的浓度分布图。但是没有考虑土壤的自我降解能力，而且由于采样点只在土壤表层采样，因此无法考虑重金属污染元素在纵向的扩散情况。而在实际情况下，这些因素对研究土壤的污染程度有一定的影响。 5.4.2需采集的信息与模型的建立

在第三问模型中，由于所给信息量不足，建模过程中我们根据三维扩散模型的理论建立了二维稳态污染物扩散模型，未考虑重金属污染物的纵向扩散以及假设不存在土壤的降解系数，这显然误差较大，实际效用不是很大，对该城区的地质环境的演变预测不能提供更好的研究，为此要更好的研究城市地质环境演变模式，我们还需要以下采集信息：在土壤的不同深度进行采样；在三年内各个季度中旬进行一次取样，并统计本季度天气情况以及自然环境的变化；采集土壤对重金属元素的降解能力。

由于土壤中微生物作用而使土壤具有一定得降解物质的能力，通过统计各区域有机质的类型，根据土壤有机质分解动力学可计算得出土壤的降解能力K。在问题三中，二维稳态污染物质浓度场基本模型的建立过程中若加上土壤的降解能力K这个信息，就可以知道从t时刻到tt时刻水域内的污染物减少量为：

M3ttt[KC(x,y,z,t)]dxdydzdt （20）

那么从t时刻到tt时刻水域内滞留的污染物的质量

MM1M2M3 （21）

即得：

C2C2C2CCCCDx2Dy2Dz2uxuyuzKCtxyzxyz （22）

16

式（22）与问题三中的模型有类似的初始条件和边界值，因此可以计算出上式的解析表达式：

2M0(xuxt)2(yuyt)(zuzt)2C(x,y,z,t)exp[]exp(Kt)1.54Dxt4Dyt4Dzt8(t)DxDyDz （23）

本模型也是已知采样点处污染物的浓度要放过来求污染源的位置，而这又是一个关于时

间的函数。因此很自然的求解思路是将式(23)中的时间取反，即让时间步长为负数，从而再现污染物的传播历史过程，直至找到污染源。然而，简单地把时间步长取为负数会带来数值稳定性问题。

M0(xux(t))2C(x,y,z,)exp[08((t))1.5DDD4D(t)xxyzt(yuy(t))24Dy(t)(zuz(t))2]exp(K(t))d4Dz(t) （24）

其中为小于t的某个时刻。

6、参考文献

[1]白世强，城市土壤重金属污染监测与磁响应和生态风险研究，浙江大学，博士学位论文，2007.04.01。

[2]李湘凌，合肥义城地区土壤重金属污染评价中典型插值方法的对比，吉林大学学报，41：222—227，2011。

[3] Muller G．Index of geoaccumulafion in sediments of the Rhine river[J]．Geojoumal，2：108，1969。

[4]宋叶志，MATLAB数值分析与应用，北京：机械工业出版社，2009。

[5]卓金武，MATLAB在数学建模中的应用，北京：北京航空航天大学出版社，2009。

17