2011 高教社杯全国大学生数学建模竞赛
承 诺 书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网
上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的
资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参
考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规
则的行为，我们将受到严肃处理。
我们参赛选择的题号是（从 A/B/C/D 中选择一项填写）： A
我们的参赛报名号为（如果赛区设置报名号的话）： 3979
所属学校（请填写完整的全名）： 广东金融学院
参赛队员 (打印并签名) ：1. 蔡宗奇
2. 曾晓骏
3. 陈友辉
指导教师或指导教师组负责人 (打印并签名)： 周雪刚
日期： 2011 年 9 月 12 日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2011 高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
评
阅
人
评
分
备
注
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
城市表层土壤重金属污染分析
摘要
本文利用某城区土壤样本数据对土壤的重金属污染情况进行分析，得到 8 种重金属
元素在该城区的空间分布以及不同区域内重金属的污染程度，然后根据数据分析的结果
推断出重金属污染的主要原因，并建立模型确定该城区污染源的大概位置，最后在信息
的收集足够的前提下，建立出新的模型，对城市地下水污染进行了分析与预测。
对于问题一，首先，我们对数据进行分析及合理假设，利用三次插值法进行插值，
绘制出 8 种重金属元素在该城区的等浓度线图。然后我们结合单项污染指数法和内梅罗
综合污染指数法对该城区重金属元素的分布进行评价，得出了该城区的土壤中重金属元
素浓度均超标且受到污染的样点数达到总样点数的 %的结论，其中受到污染范围最
大的为交通区，其次分别为生活区、工业区、山区和公园绿地区。
对于问题二，我们运用了因子分析法对剔除特大值和特小值后的数据进行分析。首
先我们对数据进行标准化处理，计算标准化后的相关系数阵、特征值和特征向量。然后
对数据进行正交变换，使用方差最大法，确定了 6 个主要因子后计算因子得分，统计分
析找出了 6 个因子的联系，得出了 Cr 和 Ni 污染主要是由工厂排污排气造成，Cd 和 Pb
污染主要是交通尾气所造成，Hg 污染可能为交通工具中的汞蒸汽所致，As 污染可能与
农业生产中过度使用农药所致，Cu 污染主要是居民生活与交通工具所造成，Zn 污染多
为厂矿企业的三废排放以及交通工具尾气排放所造成。
对于问题三，首先我们分析了重金属污染物的传播特征，并建立了污染源扩散模型，
通过观察某个静态点不同时段的图片，得出了污染源扩散浓度在二维空间上符合正态分
布的结论，在此基础上，我们通过合理假设与推断，建立了高斯扩散模型
最后通过最小二乘法对曲线进行拟合，估算出参数的值，得出了各重金属元素污染源的
坐标。
对于问题四，考虑到问题三中我们建立的模型的缺点主要是数据不够多，导致确认
到的污染源只存在于理想的情况下。我们认为应该采集更多的样点数，以减少误差，使
计算得到的污染源位置更精确。为了更好地研究城市地质环境的演变模式，我们选择了
作为水圈中重要组成部分的地下水作为研究的一个侧面，此外还应收集该城区的历年来
的垃圾掩埋量、生活废水排放量、工业废水排放量、工业废气排放量、地下水总硬度超
标率以及各种重金属离子含量浓度超标率和浓度，运用灰色关联度分析法并建立了
GM(1,1)模型来对地下水水质进行模拟。
关键词：单项污染指数 内梅罗综合污染指数 因子分析法 高斯模型 GM(1,1)模型
2 2
0 0
2 2
( ) ( )1
ln ln ln ln ( )
22 y x y x
y y x xq
c
u
一、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响
日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开
展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关
注的焦点。
按照一般城市的功能划分，城区一般可分为生活区、工业区、山区、主干道路区及
公园绿地区等，分别记为 1 类区、2 类区、……、5 类区，不同的区域环境受人类活动
影响的程度不同。
对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距 1 公里左
右的网格子区域，按照每平方公里 1 个采样点对表层土（0~10 厘米深度）进行取样、
编号，并用 GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的
多种化学元素的浓度数据。另一方面，按照 2 公里的间距在那些远离人群及工业活动的
自然区取样，将其作为该城区表层土壤中元素的背景值。
现有采样点的位置、海拔高度及其所属功能区等信息以及 8 种主要重金属元素在采
样点处的浓度和 8 种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出 8 种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金
属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收
集什么信息？有了这些信息，如何建立模型解决问题？
二、问题分析
对于问题一的分析
问题一要求我们给出 8 种主要重金属元素在该城区的空间分布，这要求我们对题中
所给数据进行分析后画图。首先要对表中数据进行分析，检查数据是否存在异常的情况。
然后运用 Matlab 软件对数据进行插值，参考地理上用的等高线图，画出 8 种元素的等
浓度线图，就可以看出 8 种重金属元素在城区的分布。分析该城区内不同区域重金属的
污染程度，我们选择可以计算某一评价指标的单因子污染指数法和兼顾极值或称突出最
大值的计权型多因子环境质量指数内梅罗综合污染指数法，通过计算，得出该城区内不
同区域重金属的污染程度。
对于问题二的分析
问题二要求我们通过数据分析来说明重金属污染的主要原因，这就要求我们对题目
中所给的数据进行深层的分析。运用因子分析法，首先，我们对该地区的重金属元素含
量浓度的样本数据进行剔除特大值和特小值处理，再进行标准化处理，以消除不同变量
的量纲的影响。然后，计算标准化后的数据的相关矩阵，求出相关系数的特征值和特征
向量，以得到各个元素之间不同的相关性。使用因子分析法的关键就是利用相关系数矩
阵求出相关因子的特征值和累计贡献率。为了使因子载荷两级分化，我们需要对原始数
据矩阵进行正交变换，使用方差最大法。最后，通过确认因子个数，计算各因子得分进
行统计分析，结合实际情况，就能够说明重金属污染的主要原因。
对于问题三的分析
根据重金属污染物的传播特征，我们可以得知 8 种重金属污染物的传播特征，选择
其中一种扩散模型进行分析，建立污染源扩散模型，分析重金属污染物传播的特点与形
态，建立相应的模型来确定污染源的位置。
对于问题四的分析
由于问题三所建立模型在缺失数据的情况只考虑到一种传播情况，所得污染源位置
是建立在最理想化的基础上的，因此必须获得必要的数据如当地水流、风速等情况才能
进一步分析与检验该污染源的位置。在获得所需要的数据后，我们可以建立新的模型来
研究该地区的地质环境变化。
三、模型假设
1.假设题目中所给的数据都真实可靠；
2.假设该地区在近期以及未来一段时间都没有发生任何可能影响数据大小的灾害或人
为的改造，如地震或建造大型公共设施等；
3.把传染源传播污染物质的过程看作是在污染源所在的位置向四周等强度瞬时释放金
属污染物质，不计风力和大地的影响；
4.传播途径中，污染物本身是被动的、保守的，即污染物和空气无相对运动；
5.污染物浓度在各风向上分布为正态分布；
6.全部高度风速均匀稳定；
7.源强是连续均匀稳定的；
8.扩散中污染物的质量是守恒的（不考虑转化）。
四、符号说明与名词解释
符号说明
符号 说明
i 角标，i=1，…，n
j 角标，i=1，…，n
区域重金属 i 的单项污染指数
第 i 种重金属浓度的实测值（μg/g 或者 ng/g）
第 i 种重金属浓度的评价起始值（μg/g 或者 ng/g）
综合污染指数（综合反映各污染物对区域土壤的不同作用）
所有单项污染指数的平均值
土壤环境中各单项污染指数中的最大值
第 i 个样本的第 j 个指标值
第 j 个指标的均值
第 j 个指标的标准差
c 空间点（x，y，z）的污染物浓度（ng/g）
q 源强，单位时间污染物排放量(ng/s)
u 平均风速(m/s)
烟气扩散系数，与大气稳定度和水平距离有关，并随 x 的增大而增加
iP
iC
iS
P
综
P
平均
maxP
ijx
jx
j
y
名词解释
重金属污染[1]：重金属系指密度 以上约 60 种元素或密度在 以上的 45 种元
素。砷、硒是非金属，但是它的毒性及某些性质与重金属相似，所以将砷、硒列入重金
属污染物范围内。环境污染方面所指的重金属主要是指生物毒性显著的汞、镉、铅、铬
以及类金属砷，还包括具有毒性的重金属锌、铜、钴、镍、锡、钒等污染物。
五、模型的建立与求解
问题一的求解
对数据源进行处理
对题目所给数据进行分析，发现数据中有个别几个数据与其他数据相比显得特别巨
大，比如 Hg（汞）元素在编号为 9 的点的数值为 16000，明显比所测数据都大，但是通
过分析我们认为，这几个点有可能为污染源或者垃圾掩埋点，所以并不把它们认定为异
常数据。通过对整体数据分析后，我们发现除了上述几个特别的样点外，其他数据在正
常范围之内。
8 种主要重金属元素在该城区的空间分布
由上述，我们排除了数据源异常的情况，根据题目所给出的数据，我们利用 Matlab
软件（程序代码见附录一）画出该城区的各功能区域的散点图（见图 1）。
图 1 该城区各功能区散点分布
接着，根据题目所给数据，画出 8 种重金属元素在该城区的等浓度分布图（程序代
码见附录二）
0 1 2 3
x 10
4
0
1
2
x 10
4 o-生 活 区 *-工 业 区 + 山 区 ◇ -交 通 区 星 形 -绿 地 区
图 2 As 元素等浓度分布图
图 3 Cd 元素等浓度分布图
As等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
5
10
15
20
25
Cd等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
200
400
600
800
1000
1200
1400
1600
图 4 Cr 元素等浓度分布图
图 5 Cu 元素等浓度分布图
Cr等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
100
200
300
400
500
600
700
800
900
Cu等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
0
500
1000
1500
2000
2500
图 6 Hg 元素等浓度分布图
图 7 Ni 元素等浓度分布图
Hg等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
0
2000
4000
6000
8000
10000
12000
14000
16000
Ni等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
20
40
60
80
100
120
140
图 8 Pb 元素等浓度分布图
图 9 Zn 元素等浓度分布图
由上述八种重金属元素的等浓度图，我们可以很直观的看出这八种重金属元素在各
Pb等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
50
100
150
200
250
300
350
400
450
Zn等 浓 度 分 布 图
0 1 2
x 10
4
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
0
500
1000
1500
2000
2500
3000
3500
个区域的浓度高低的区间分布。
该城区各功能区的重金属污染程度分析
由于该城区的重金属污染包括 8 种重金属元素污染，因此我们考虑利用单项污染指
数法来对各种元素污染程度进行评价，然后再利用内梅罗综合污染指数法综合评价该城
区总体重金属的污染程度。
（1）单项污染指数法
为了分析该城区各种重金属的污染程度，我们利用如下单项污染指数法模型：
，
其中 为区域第 i 种重金属的单项污染指数； 为第 i 种重金属元素含量的实测值；
为第 i 种重金属含量的评价起始值。若 >1,则表示该区域受到污染。
区域土壤背景值代表了未曾受到或者相对受到“三废”污染较少的土壤中有毒物质的
平均值，具有显著的区域特点[4]。土壤总量污染起始值评价是以区域土壤中某污染物浓
度的平均值加 2 倍标准差所得结果作为区域土壤污染的起始值[5]，并将其作为评价标准
对各采集样点分析结果进行评价。其具体步骤是，利用题目给出的 8 种主要重金属元素
的背景值数据（见表 1）。
表 1 重金属总量的土壤污染起始值计算结果
重金属
元素
含量平均值
（μg/g 其中 Hg 为 ng/g）
标准差 含量起始值
（μg/g 其中 Hg 为 ng/g）
As
Cd 130 30 190
Cr 31 9 49
Cu
Hg 35 8 51
Ni
Pb 31 6 43
Zn 69 14 97
以平均值加 2 倍标准差作为 ，利用单项污染指数法模型计算得到该城区重金属元
素总量的单项污染指数（见表 2）。
表 2 重金属元素总量的单项污染指数值统计结果
单项污染指数重金属
元素
样品数量/个
最大值 最小值 平均值
超标率%
As 319 49
Cd 319 65
Cr 319 34
Cu 319 70
Hg 319 48
Ni 319 22
Pb 319 55
Zn 319 57
由表 2 可知，该城区的土壤重金属元素都有超标，其中以 Cu 的含量超标最为严重，
为 70%,其他七种元素也超标，部分较为严重。
根据以上数据计算得到各个功能区重金属元素的单项污染指数（见表 3）。
i
i
i S
C
p
iP iC iS
iP
iS
表 3 各区单项重金属的单项污染指数值统计结果
超标率%功能区
As Cd Cr Cu Hg Ni Pb Zn
样本个数 319 319 319 319 319 319 319 319
生活区
工业区
山区
交通区
公园绿地区
由上表可得，各区的各种重金属元素均超标，其超标率最高的元素统计如表 4
表 4 各个功能区超标率最高的元素
功能区 生活区 工业区 山区 交通区 公园绿地区
元素 Cu % Cu % Cd/Cu
%
Cu % As %
根据表 4，得知该地区的 Cu 元素超标十分严重，严重影响到了当地的土壤。
（2）内梅罗综合污染指数法。
为了综合评价该城区总体重金属的污染程度，我们选择了内梅罗综合污染指数法[2]
进行分析。 内梅罗综合污染指数法的计算公式为：
其中 为综合污染指数（综合反映各污染物对区域土壤的不同作用）； 为所有单项
污染指数的平均值； 为土壤环境中各单项污染指数中的最大值。内梅罗指数[2]是一
种兼顾极值或称突出最大值的计权型多因子环境质量指数，内梅罗综合指数的分级标准
[3]如表 5：
表 5 内梅罗综合污染指数的分级标准
污染等级 污染程度
1 安全 清洁
2 警戒线 尚清洁
3 轻污染 超标
4 中污染 土壤、作物受中度污染
5 重污染 土壤、作物受严重污染
2
1
2
max
2
2
PP
P 平均
综
综
P
平均
P
maxP
综
P
综
P
综
P
综
P
综
P
综
P
表 6 为采用内梅罗综合污染指数法对区域土壤中重金属总量污染起始值统计计算后
的结果。
表 6 内罗梅综合污染指数法统计结果
评价分级
评价标准 总样品
数/个
清洁
/个
尚清
洁/个
轻度污
染/个
中度污
染/个
重度污
染/个
清洁、
尚清
洁所
占比
例%
轻度、
中度、
重度污
染所占
比例%
重金属总
量污染起
始值评价
319 16 46 122 53 82
由表 6 可知，采用重金属污染起始值评价时，轻度、中度和重度污染的样点数达到
了 %，而清洁、尚清洁的样点的约占 %。
表 7 受轻度、中度和重度污染区域在各个功能区的分布统计
功能区 1 2 3 4 5
所占比例%
由表 7 可见，在 5 个功能区中，受到污染最多的为交通区，其次为生活区、工业区、
山区和公园绿地区。
综上，我们得知，该地区的污染问题十分严重，所测的 8 种重金属元素均超标，且
超标率最高达到了 70%，而 %的污染地区位于该地区的交通地带，其次是生活区、
工业区、山区和公园绿地区。
问题二的求解
按照问题分析中对问题二的分析思路，借助 SAS 统计软件的 insight 模块运用因子
分析法[6]对数据进行分析。在因子分析数学模型中，通过正交的方差最大旋转法使每一
个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小。变量或因子
的重要程度都是以其方差大小来衡量。因子旋转后每个变量因子负荷代表着在系统中的
作用或重要性程度，以各个变量目标因子载荷平方与因子方差贡献率乘积作为变量的权
重，构成一个判别污染来源的综合指标。
数据标准化处理
对该地区重金属元素含量浓度的样本数据做剔除特大值与特小值后的处理，为了消
除不同变量的量纲的影响，我们对数据进行了标准化处理，而且标准化转化不会改变变
量的相关系数。
标准化公式为：
其中 为第 i 个样本的第 j 个指标值，而 和 分别为 j 指标的均值和标准差。
计算标准化数据的相关系数阵，求出相关系数的特征值和特征向量
首先，给出该城市表层土壤 As、Cd、Cr、Cu、Hg、Ni、Pb 和 Zn 八种重金属元素
原始含量浓度数据的相关系数矩阵。如图 10 所示：
j
jij
ij
xx
x
ijx jx j
图 10 变量相关矩阵
可见，Cr 和 Ni 的相关性最好，相关系数最大，为 ，其次为 Pb 和 Cd，相关
系数为 ，以下依次是 Cu 和 Cr，Pb 和 Cu 的相关性较好，相关系数分别是
和 ，Ni 和 Cu，Zn 和 Pb 的相关系数分别为 和 。其他元素之间的相
关性并不是很好。从成因上来分析，相关性较好的元素可能在成因和来源上有一定关联。
因子分析的关键就是利用相关系数矩阵求出相关的因子的特征值和累计贡献率，
用 SAS 统计软件计算可得出，见图 11：
图 11 特征值和累计贡献率
在累积方差为 %（90%）的前提下，分析得到 6 个主因子，可以看出 6 个主因
子提供了源资料 %的信息，满足因子分析的原则。从上图还可以看出因子 1 可能
为该城市土壤重金属污染的最重要的污染源，对该城市重金属污染的贡献最大，因子 2、
因子 3、因子 4、因子 5、因子 6 对该城市重金属污染有重要作用。
进行正交变换，使用方差最大化法
其目的是使因子载荷两级分化，而且旋转后的因子仍然正交。
图 12 旋转前因子载荷矩阵
图 13 方差极大正交旋转后因子载荷矩阵
图 14 旋转后的因子载荷矩阵
确定因子个数，计算因子得分，进行统计分析
由以上图 12、图 13 和图 14 可见，旋转前后因子载荷的变量结果基本一直。变量与
某一个因子的联系系数绝对值（载荷）越大，则该因子与变量关系越近。正交因子解说
明：因子 1 为 Cr 和 Ni 的组合，因子 2 为 Cd 和 Pb 的组合，因子 3 为 Hg,因子 4 为 As,
因子 5 为 Cu，因子 6 为 Zn。Cr 和 Ni 、Cd 和 Pb 可能是同一个来源，而且这两组元素
正好是相关性最好的两组元素。
结论
因子 1 为元素 Cr 和 Ni 的组合，结合第一题的重金属元素等浓度线图，我们可以得
出以下结论，该城市表层土壤 Cr 和 Ni 基本没污染，只有个别点富集程度较高。这个可
能是由于该富集中心位于工厂排污口或排气口附近。
因子 2 为元素 Cd 和 Pb 的组合，对两种元素的重金属等高浓度线图进行对比我们可
以看出两种元素在空间分布上为集中于一个带状的污染带。且两种元素的富集中心主要
坐落于交通区。这个主要因为交通区汽车等交通工具所排放的尾气污染物中均含有铅元
素和铬元素等。所以我们可以认为因子 2 的两种重金属污染是有交通工具的尾气排放造
成的。
因子 3 为元素 Hg，从该元素的等浓度线图上分析，该元素的富集中心主要集中于交
通区，这也是因为交通工具的尾气排放中的汞蒸汽所导致的。
因子 4 为元素 As，从该元素的等浓度线图上分析，该元素的富集区中心为于山区，
而砷是农作物农药当中的主要元素。所以我们可以认为砷污染的主要原因是由于山区进
行农作物种植时农药的过度使用。
因子 5 为元素 Cu,从该元素的等浓度线图上看，该元素的富集中心位于生活区与交
通区附近，这可能是因为生活区的商业活动、居民生活和居民交通工具的使用所导致的
土壤中的铜污染。
因子 6 为元素 Zn,从该元素的等浓度线图上看，锌元素的富集中心位于工业区于交
通区相间的区域，这可能是因为厂矿企业的三废排放和交通工具尾气的排放所导致的锌
污染。
针对问题三的建模与求解
重金属污染物的传播特征分析
根据重金属污染物的传播特征，可以初步将所有传播特征分成五类：自然沉降和雨
淋沉降（中心向四周扩散型）、含重金属元素物质的燃烧（带状分布有叠加性）、干湿沉
降（吸附或固定使富集）、人为的不合理使用（过量形成富集）、其他，而结合题中所给
的 8 种重金属污染物的特点可知，根据以上分类原则，一种重金属可能同时属于几种不
同的传播特征类，因此我们先仅按其主要传播特征将其归类，考虑前四大类，可以得到
下表 8：
表 8 8 种主要金属污染物的传播特征
种类 传播特征
一 自然沉降和雨淋沉降（中心向四周扩散型）
二 含重金属元素物质的燃烧（带状分布有叠加性）
三 干湿沉降（吸附或固定使富集）
四 人为的不合理使用（过量形成富集）
因此根据上表的各类重金属污染物的传播特点，由第一类污染物是从污染源中心向
四周扩散传播的，可以知道该类最适用于求解污染源，即只要找到污染物浓度最高的地
方就是所要寻找的污染源。
污染源扩散模型
扩散的数学模型基于一个基本假设[9]，即：在溶液中，溶质的扩散遵循 Fick 第一定
律，用数学表达式表示为：
， ()
式中， 为溶质在溶液中的扩散通量， 表示溶质分子在溶液中的扩散系数， 表
示的是浓度梯度，负号表示溶质分子沿着浓度降低的方向扩散。
在带有重金属元素的大气中，重金属元素的扩散同样遵循 Fick 第一定律，即：
， ()
式中，F 表示的是元素在大气中的扩散通量，D 表示的是元素在大气中的扩散系数。
根据质量守恒定律，可导出物质扩散的基本微分方程：
， ()
式中， ， ， 为重金属元素在 x，y，z 方向的扩散通量。
当 D 为常数时，将 ， ， 代入式()，得到污染源扩散的基本微分方程：
， ()
由于土壤溶质的扩散模式复杂多样，为了讨论的方便，我们不对其它因素加以考虑，
例如溶质自重、边界条件、空间形状等，只对最基本的扩散情况（即瞬时点污染源）进
行讨论。其他的情况大多可通过对其在时间和空间上的积分来求得。
设 t =0 时，在(x，y，z)处污染物质量为 M，则其扩散方程为：
， ()
若坐标原点在污染源发生处时，扩散方程的解析解为：
， ()
根据上式，当 时，扩散方程的解析解为：
()
从式()中可看出，污染源扩散浓度与扩散系数 D、坐标位置(x,y,z)以及污染物质
量 M、时间 t 密切相关，根据研究的需要，我们主要关注污染源浓度与空间位置之间的
关系，对其他因素进行简化或采取赋值的方式进行处理，以突出研究的重点。
重金属污染物的传播特征
为了进一步研究重金属污染物的扩散规律，给出时间 t、扩散系数 D、污染物质量 M
分别取某个静态点的图片：
X
C
DF
00
0F 0D X
C
X
C
DF
0
z
F
y
F
x
F
t
C zyx
xF yF zF
xF yF zF
)(
2
2
2
2
2
2
z
C
y
C
x
C
D
t
C
2
2
2
2
2
2
z
C
D
y
C
D
x
C
D
t
C
zyx
)
2
1
)(
2
1
)(
2
1
(),,,( 444
222
tD
z
z
tD
y
y
tD
x
x
zyx e
tD
e
tD
e
tD
Mtzyxc
zyx DDD
)
4
exp(
)(8
),,,(
222
2
3 Dt
zyx
Dt
M
tzyxc
图 15 D=1，M=1，t=1 浓度分布模拟图 图 16 D=1，M=1，t=5 浓度分布模拟图
根据图 15 和图 16 可看出：在污染物质量和扩散系数相同的情况下，瞬时点污染源
的扩散范围随时间的增加而不断扩大，从图中色轴及标注的数值可以看出其浓度也在不
断降低。
图 17 D=1，M=100，t=5 浓度分布模拟图 图 18 D=5，M=100，t=5 浓度分布模拟图
从图 17 和图 18 了解到：在扩散质量和时间相同的情况下，扩散系数越大，扩散范
围越大，污染源浓度也越低。
由上述分析可以得出，污染源扩散浓度在二维空间上符合正态分布，其二维正态分
布的联合概率密度函数为：
建立模型
在这里，我们选择建立高斯扩散模型[10]（即中心向四周扩散模型）。
1.高斯模式的有关假定
设定坐标系。坐标系取排放点（无界源、地面源或高架源排放点）在地面的投影点
为原点，主风向为 x 轴，y 轴在水平面内垂直于 x 轴，正方向在 x 轴的左侧，z 轴垂直
于水平面，向上为正，即右手坐标系。作如下假设：
1）把传染源传播污染物质的过程看作是在污染源所在的位置向四周等强度瞬时释
放重金属污染物质，不计风力和大地的影响；
2）传播途径中，污染物本身是被动的、保守的，即污染物和空气无相对运动；
3）污染物浓度在各风向上分布为正态分布；
4）全部高度风速均匀稳定；
5）源强是连续均匀稳定的；
6）扩散中污染物的质量是守恒的（不考虑转化）。
2
2
2
2
21
21
2
1
2
1
22
21
)())((2)(
)1(2
1
exp·
12
1
),(
yyxx
yxf
2.无界空间连续点源扩散模式
由正态分布假定，得下风向任一点的浓度分布
，
方差的表达式
， ，
由假定条件，得出源强积分式
（单位时间物料守恒）
其中浓度 c 为未知数， z=0 时, a,b (= )为待定系数,分别与 x,y 有关函数
积分后，可以解出四个未知数：得到高斯模型
，
其中，c 表示空间点（x，y，0）的污染物浓度，ng/g;q 表示源强，单位时间污染物排放
量，ng/s;u 是平均风速，m/s； 、 为烟气扩散系数，与大气稳定度和水平距离有关。
结合本题所给数据，我们对高斯模型两边取以 e 为底的对数，得：
，
进而，我们可以推出关于污染源原点 的模型：
模型求解
由于该模型中需要确认的参数有 ， ， ， ， ，故采用曲线拟合来确
定这些参数，使得拟合值与真实值之间的误差达到最小。而常用的判断曲线拟合效果好
坏的办法就是最小二乘法，因此我们运用最小二乘法来对曲线进行拟合。
根据各类重金属污染物的传播特点，只要找到污染物浓度最高的地方就是所要寻找
的污染源。因此选择样本点时，应该选择在平均浓度较高，也就是在元素等浓度线图中
浓度较高且集中的区域进行选点。
运用 软件对所求参数进行估算（运行程序见附录三），从而得到传染源的
位置如表 9
表 9 各金属传染源坐标
重金属元素 As Cd/Pb Cr/Ni Cu
污染源点坐标 （12715，3042） (21349,25055) (3712,4996) (11357,39560)
污染源浓度 1400ng/g
重金属元素 Hg Zn
污染源点坐标 （7167，2054） (12701,2843)
污染源浓度 3700ug/g
2 2
( , , ) ( ) ay bxc x y z A z e e
2
2 0
0
y
y cdy
cdy
2
2 0
0
x
x cdx
cdx
q ucdydx
2
1
2
2 2
2 2( , ,0) exp[ ( )]2 22 y xy x
q y x
c x y
u
y x
2 2
2 2
1
ln ( , ,0) ln ln ln ( )
22 y x y x
q y x
c x y
u
0 0 0( , , )x y z
2 2
0 0
2 2
( ) ( )1
ln ln ln ln ( )
22 y x y x
y y x xq
c
u
ln
2
q
u y
z 0y 0x
针对问题四的建模与求解
对于解决问题三时所建立的模型，该模型的优点是所需要的数据量相对较少，而缺
点是该模型建立的前提是忽略了天气、地貌和水流等影响。所求的的结果是在相对理想
化情况才存在的，在实际中借鉴意义不大。
地质环境是自然环境的一种，是指由岩石圈、水圈和大气圈组成的环境系统。在长
期的地质历史演化的过程中，岩石圈和水圈、岩石圈和大气圈之间、大气圈和水圈之间
进行物质迁移和能量转换，组成了一个相对平衡的开放系统。人类和其他生物依赖地质
环境生存发展，同时，人类和其他生物又不断改变着地质环境。而地下水[7]作为水圈中
极其重要的组成部分，由于水量稳定，水质好，一直是农业灌溉，工矿和城市的重要水
源之一。随着经济的发展，人类生产活动中排出的有毒有害物质使得地下水质持续恶化。
所以我们选择地下水污染作为研究城市地质环境演化模式的一个侧面，对地下水进行水
质模拟，建立模型以便对地下水域内的地下水域状况进行了解。
为了对地下水水质进行模拟，我们需要采集该城市历年来的垃圾掩埋量，生活废水
排放量，工业废水排放量，工业废气排放量，地下水总硬度超标率，各种重金属离子含
量浓度超标率和浓度，
地下水污染成因分析
由于地下水水质受到气候，地下水位，污染源的存在以及互相影响非线性的关系。
所以这里我们采用灰色关联度分析法来分析产生地下水污染的成因。
(1)建立参考数列和比较数列
设历年的地下水总硬度数据为参考数列 ，设垃圾掩埋量数据 ，生活废水排放量
数据 ，工业废水排放量数据 和工业废气排放量数据 作为比较数列。
(2)对各类列进行初值化处理
(y=1,2,3,4)
(3)计算参考数列在 k 点的灰色关联系数：
(i=1,2,3,4)，
为分辨系数，具体取值规则如下[11]：
记 为所有差值绝对值的均值，即
记 ， 则 的 取 值 为 ： 当 >3 时 ， ；
。
该模型侧重于总体分析， 不仅可以调节 的大小，而且可以调节它的变化区
间， 的下界值随 增大而增大，下界值增大说明区间变小，分辨率变低，分辨效
果则不明显。根据因素见的关联分析可选择不同的分辨系数，在计算时，一般取 =
便可得到满意的分辨率。当 则无法进行关联分析，此时所有 ，即关联
系数蜕变为一个点[8]。
(4)计算灰色关联度
oy 1y
2y 3y 4y
)(max
)(
ky
ky
y
j
j
j
)()(maxmax*)()(
)()(maxmax*)()(minmin
)(
00
00
kykykyky
kykykyky
k
ikii
ikiiki
i
v
m
i
n
k
iv kykynm 1 1
0 )()(*
1
max
v max v
时，v
)(ki
)(ki
1)( ki
(i=1,2,3,4)
(5)排关联序
在得出各个序列的灰色关联度排序后我们就可以比较出城市地下水总硬度这一指标
与垃圾掩埋量， 生活废水排放量，工业废水排放量和工业废气排放量中那一个指标关
联系数较高。从而得出地下水总硬度超标率的主要原因。依照上述(1)-(5)的做法我们同
样可以用灰色关联度分析法对造成城市地下水各种重金属离子超标的成因做出判断。
地下水污染预测
地下水中的污染是一个包含了物理作用，化学作用和生物作用的过程。由于作用复
杂，我们可以利用灰色系统理论，把地下水环境视为灰色系统，他污染物看作是在一定
范围内变化的灰色变量。从而我们可以建立起 GM（1，1）模型[12]。
（1）GM（1，1）模型建立原理。
设原始数列为 为原始非负时间序列，预测的
结果就是想知道
为累加声称序列，即 =
具体地说，就是：
GM(1,1)模型的白化微分方程为：
a 为待辨识参数；u 为待辨识内生变量；设待辨识变量 ，按最小二乘法得：
，
其中 ，y= 。
于是我们可以得到灰色预测模型为：
(2) 残差检验和后残差检验：
模型建立后，我们还应该对模型进行残差检验和后验差检验。残差检验有两种，绝
对误差和相对误差检验：
绝对误差
，
相对误差
n
k
ii kN
r
1
)(
1
)}(,),2(),1({ )0()0()0()0( nQQQQ
),2(),1( )0()0( nQnQ
)()0( tQ )()1( tQ niiQ
t
m
3,2,)1(
1
)0(
niiQiQiQ
QQ
i ,,2),1()()(
)1()1(
)()0()1(
)0()1(
uaQ
dt
dQ
)1(
)1(
u
a
a
yBBBa TT 1)(
1))()1((
2
1
-
1
1))3()2((
2
1
1))2()1((
2
1
-
)1()1(
)1()1(
)1()1(
nQnQ
QQ
QQ
B
)(
)3(
)2(
)0(
)0(
)0(
nQ
Q
Q
nt
a
u
e
z
u
QtQ at 3,2,1,))1(()1( )0()1(
),2,1(),()()(
)0(
)0()0( niiQiQi
，
其中 。
后残差检验：
首先计算原始数列 ，其定义为：
，
，
，
然后计算残差数列 ,其定义为：
，
，
，
由此计算方差比 C= ，最后计算小误差概率 ，
。
根据下面的预测精度等级划分表，如表10所示
表10 预测精度等级划分表
小误差概率p值 方差比c值 预测精度等级
> < 好
> < 合格
> < 勉强合格
不合格
最后我们将各种重金属离子含量浓度的超标率代入模型中，对将来该城市地下水的
污染进行预测。
六、模型评价
本文通过建立单项污染指 数法模型和内梅罗综合指数评价模型，对该地区土壤中
的重金属元素进行评价分析，从而客观地了解到土壤中的 8 种重金属元素的分布、超标
率以及各重金属原子分布来说明不同区域的重金属污染程度。
由于题目所给数据较少，因此拟合出的曲线和污染源点的位置并不是其准确位置，
而在实际中，突发事件也是必须考虑到的，情况将更加复杂，而该模型并不能很好地解
决突发事件对该地区影响，因而这也是该模型需要改进的方向。
),2,1(%,100*
)(
)(
)(
)0(
)0(
)0( ni
iQ
i
i
),2,1(),1()1(),1()()(
)1()0()1()1(
)0( nIQQiQiQiQ
0
)0( SX 的均方差
1
2
0
0
n
S
S
2
1
)0()0(2
0 )(
n
i
QiQS
n
i
iQ
n
Q
1
)0()0( )(
1
1
)0( S的均方差
1
2
1
1
n
S
S
2
1
)0()0(2
1 )(
n
i
iS
n
i
i
n 1
)0()0(
)(
1
0
1
S
S
0)0()0( S
七、参考文献
[1]佚名，重金属污染来源、分布、治理方法，第一页，出版年不详
, 2011-9-11.
[2]内梅罗综合指数分析法， 百度百科， ,
2011-9-11．
[3]HJ/T 166-2004,土壤环境监测技术规范.
[4]姜芝萍，杨俊衡，城市重点污染区土壤重金属污染评价标准探讨——以衡阳市某区
为例，安全与环境工程，第 17 卷第 1 期:第 59-60 页.
[5]蔡信德，郭杨，不同标准对城市土壤重金属质量分数的评价,环境科学研究，2005，
26（5）:第 191-197 页.
[6]王雄军，基于因子分析法研究太原市土壤重金属污染的主要来源，生态环境，2008，
17（2）:第 671-676 页.
[7]地下水，百度百科，
[8]吕锋 灰色系统关联度之分辨系数的研究,系统工程理论与实践，1997.
[9]李明，吴超，粉尘点污染扩散模型的可视化研究，中南大学资源与安全工程学院.
[10]黄金杰，杨桂花，马骏驰，基于高斯大气污染评价模型，计算机仿真，2011，28
（2）:
第 101-106 页.
[11]吕锋,灰色系统关联度之分辨系数的研究系统工程理论与实践，1997，第 51 页.
[12]温丽华，灰色系统理论及其应用，计算机应用技术，2003,第 15-17 页.
附录
附录一
qqu1=[4043 1895 2427 3971 4777 4897 6534 5641 4592 4603 2486 5999 3573
6213 5375 8643 7304 5230 9328 4311 8077 6401 7056 8348 15467 8658 14844
5519 16569 6055 14298 7418 21418 10721 26453 5577 26416 6508 5101 4080 5382
3012 5503 1127 5636 133 6605 374 11649 3515 12591 1063 13855 3345 14862
2524 15387 729 15810 2307 17203 6218 16301 8299 17904 8287 10395 11203 15467
12080 16428 9069 16289 10072 12153 12336 11958 13313 10800 13282
9277 16148 11121 16432 12625 16259 17981 18449];
qu2=[0 1787 1647 2728 2383 3692 4742 7293 4948 7293 5567 6782
22674 12173 5438 3994 12734 4015 14896 1603 16947 7487 5006 8846
6395 10443 7405 10981 8446 11200 7612 11938 8866 13143 9475 12000
9212 11305 8629 12086 7776 10613 7106 9467 6423 8831 7458 8920
11646 9381 12641 9560 14000 8970 14207 9980 15140 11101 16440 13232
10022 12204 9333 14631 10856 14727 12644 14943 9036 17538
10599 17980
];
qu3=[14325 8666 17044 10691 18413 11721 19007 11488 18738 10921
23664 9790 18993 12371 19968 12961 22535 11293 27816 5581
25361 6423 24065 7353 25998 7032 27177 7771 26073 8807 24631 9422
24702 9522 25461 9834 26086 11094 26015 12078 27700 11609
27696 11621 27346 13331 26591 13715 27823 14737 27232 14482
24580 13319 24153 12450 22965 13535 24685 14278 28654 8755
24003 15286 21684 13101 22193 12185 17079 5894 15255 5110
15007 5535 15801 3966 7008 4775 1475 8540 20261 7586 19569 7348
19411 6934 23359 5325 3238 6502 22624 4818 21703 6591 12734 10344
14405 18032 14074 16516 14262 15129 20591 13549 20983 15862
20177 17642 18906 16346 18467 17001 16607 17365 15952 18397
22605 14301 23146 15382 22046 17634 23785 17643 25981 18051
27380 18202 23325 16701 26852 16114
];
qu4=[74 781 1373 731 1321 1791 1049 2127 2883 3617 2708 2295 2933 1767
3526 4357 5062 4339 5868 4904 5481 6004 3299 6018 5635 7965 5394 8631
5291 7349 7004 6226 7048 4600 8180 4496 9090 5365 8049 5439 8017 7210
6869 7286 7747 8260 8457 8991 9460 8311 9062 7639 9319 6799 10631 6472
10685 5528 10643 4472 11702 4480 11730 5532 11482 6354 10700 8184
10630 8774 11678 8618 11902 7709 13244 7056
12746 8450 12855 8945 13797 9621 12442 4329 13093 4339 13920 5354
15658 7594 14177 6684 12778 5799 17087 11933 17075 12924 17962 12823
17814 10707 18134 10046 17198 9810 17144 9081 18393 9183 19767 8810
21006 8819 21091 9482 22846 9149 22304 10527 21439 11383 20554 11228
20101 10774 21072 10404 20215 9951 21766 12348 25221 5795 26424 8639
24813 10799 23198 13523 4020 2990 4026 3913 5314 2060 7093 1381
7100 2449 6837 3490 7906 3978 8045 3052 8394 2035 8403 1075 8079 0
9663 1288 9469 2286 9178 3299 9095 3975 10225 3821 10210 2789
10340 1764 11557 1581 11415 2585 12696 3024 12400 2060 13765 1353
13694 2357 16032 3061 16872 2798 17734 3629 17005 7212 18438 6539
18954 4874 18012 4414 19072 8519 20282 8590 21450 7555 19501 6091
19909 5300 21018 5764 22176 5492 5734 9659 7912 12840 9296 13102
8622 10638 9237 9872 8307 9726 8904 8868 10547 9591 10398 10360
11529 11243 11563 10298 14065 10987 12727 7691 15198 10100
15248 9106 16267 11058 16440 12068 15412 12982 14269 12877
13277 13204 13175 12238 12632 17949 14624 14004 16629 14481
18470 14411 19041 15769 17414 15476 15748 15728 25021 16290
5985 2567
];
qu5=[4233 895 4741 6434 16387 6609 16061 7352 15092 6936 3518 2571
3469 2308 3762 2170 3927 2110 4153 2299 3267 793 4684 1364
5495 1205 5664 1653 5541 2093 5451 2757 15087 3512 16823 4207
18303 7385 18556 5588 20582 6548 14173 11941 15517 17034 14482 12692
14318 13569 10352 17133 9095 16414 10510 15314 13954 5615
10142 1662 17765 3561 6924 5696 4678 3765 6182 2005 7653 1952
];
qq1=qqu1(:,1);
qq2=qqu1(:,2);
qq3=qu2(:,1);
qq4=qu2(:,2);
qq5=qu3(:,1);
qq6=qu3(:,2);
qq7=qu4(:,1);
qq8=qu4(:,2);
qq9=qu5(:,1);
qq10=qu5(:,2);
plot(qq1,qq2,'o',qq3,qq4,'*',qq5,qq6,'+',qq7,qq8,'s',qq9,qq10,'p')
title('o-生活区 *-工业区 + 山区 ◇-交通区 星形-绿地区');
附录二
yuansu=[
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
];
x=d(:,1);
y=d(:,2);
z1=yuansu(:,2);
z2=yuansu(:,3);
z3=yuansu(:,4);
z4=yuansu(:,5);
z5=yuansu(:,6);
z6=yuansu(:,7);
z7=yuansu(:,8);
z8=yuansu(:,9);
nx=linspace(min(x),max(x),200);
ny=linspace(min(y),max(y),200);
figure
[xx,yy]=meshgrid(nx,ny);
zz1=griddata(x,y,z1,xx,yy,'cubic');
contour(xx,yy,zz1,20);
title('As 等浓度分布图');
figure
[xx,yy]=meshgrid(nx,ny);
zz2=griddata(x,y,z2,xx,yy,'cubic');
contour(xx,yy,zz2,40);
title('Cd 等浓度分布图');
figure
zz3=griddata(x,y,z3,xx,yy,'cubic');
contour(xx,yy,zz3,40);
title('Cr 等浓度分布图');
figure
zz4=griddata(x,y,z4,xx,yy,'cubic');
contour(xx,yy,zz4,40);
title('Cu 等浓度分布图');
figure
zz5=griddata(x,y,z5,xx,yy,'cubic');
contour(xx,yy,zz5,40);
title('Hg 等浓度分布图');
figure
[xx,yy]=meshgrid(nx,ny);
zz6=griddata(x,y,z6,xx,yy,'cubic');
contour(xx,yy,zz6,25);
title('Ni 等浓度分布图');
figure
[xx,yy]=meshgrid(nx,ny);
zz7=griddata(x,y,z7,xx,yy,'cubic');
contour(xx,yy,zz7,20);
title('Pb 等浓度分布图');
figure
[xx,yy]=meshgrid(nx,ny);
zz8=griddata(x,y,z8,xx,yy,'cubic');
contour(xx,yy,zz8,20);
title('Zn 等浓度分布图');
附录三
function f=qq(a,x)
f=a(1)-log(a(2))-log(a(3))-1/2.*(x(:,1)-a(4)).^2./a(2).^2-1/2.*(x(:,2)-a(5)).^2./a(3).^2;
%AS例子
xdata=as(1:length(as),2:3);
ydata=as(1:length(as),1);
x0=[1 1 1 1 1]; %³õÖµ
option=optimset('MaxFunEvals',10000);
[x,fval]=lsqcurvefit(@qq,x0,xdata,ydata)
qq(x,x(4:5))
