单项污染指数 内梅罗综合污染指数 因子分析法 高斯模型 GM(1,1)

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

承 诺 书

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

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

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

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

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 3979 所属学校（请填写完整的全名）： 广东金融学院

参赛队员 (打印并签名) ：1. 蔡宗奇 2. 曾晓骏 3. 陈友辉 指导教师或指导教师组负责人 (打印并签名)： 周雪刚

日期： 2011 年 9 月 12 日

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

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

编 号 专 用 页

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

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

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

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

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

摘要

本文利用某城区土壤样本数据对土壤的重金属污染情况进行分析，得到8种重金属元素在该城区的空间分布以及不同区域内重金属的污染程度，然后根据数据分析的结果推断出重金属污染的主要原因，并建立模型确定该城区污染源的大概位置，最后在信息的收集足够的前提下，建立出新的模型，对城市地下水污染进行了分析与预测。

对于问题一，首先，我们对数据进行分析及合理假设，利用三次插值法进行插值，绘制出8种重金属元素在该城区的等浓度线图。然后我们结合单项污染指数法和内梅罗综合污染指数法对该城区重金属元素的分布进行评价，得出了该城区的土壤中重金属元素浓度均超标且受到污染的样点数达到总样点数的80.6%的结论，其中受到污染范围最大的为交通区，其次分别为生活区、工业区、山区和公园绿地区。

对于问题二，我们运用了因子分析法对剔除特大值和特小值后的数据进行分析。首先我们对数据进行标准化处理，计算标准化后的相关系数阵、特征值和特征向量。然后对数据进行正交变换，使用方差最大法，确定了6个主要因子后计算因子得分，统计分析找出了6个因子的联系，得出了Cr和Ni污染主要是由工厂排污排气造成，Cd和Pb污染主要是交通尾气所造成，Hg污染可能为交通工具中的汞蒸汽所致，As污染可能与农业生产中过度使用农药所致，Cu污染主要是居民生活与交通工具所造成，Zn污染多为厂矿企业的三废排放以及交通工具尾气排放所造成。 对于问题三，首先我们分析了重金属污染物的传播特征，并建立了污染源扩散模型， 通过观察某个静态点不同时段的图片，得出了污染源扩散浓度在二维空间上符合正态分布的结论，在此基础上，我们通过合理假设与推断，建立了高斯扩散模型

q1(y?y0)2(x?x0)2lnc?ln?ln?y?ln?x?(?) 222?y?x2?u最后通过最小二乘法对曲线进行拟合，估算出参数的值，得出了各重金属元素污染源的

坐标。

对于问题四，考虑到问题三中我们建立的模型的缺点主要是数据不够多，导致确认到的污染源只存在于理想的情况下。我们认为应该采集更多的样点数，以减少误差，使计算得到的污染源位置更精确。为了更好地研究城市地质环境的演变模式，我们选择了作为水圈中重要组成部分的地下水作为研究的一个侧面，此外还应收集该城区的历年来的垃圾掩埋量、生活废水排放量、工业废水排放量、工业废气排放量、地下水总硬度超标率以及各种重金属离子含量浓度超标率和浓度，运用灰色关联度分析法并建立了GM(1,1)模型来对地下水水质进行模拟。

关键词：单项污染指数 内梅罗综合污染指数 因子分析法 高斯模型 GM(1,1)模型

1

一、问题重述

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

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

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

现有采样点的位置、海拔高度及其所属功能区等信息以及8种主要重金属元素在采样点处的浓度和8种主要重金属元素的背景值。

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

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

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

(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？

二、问题分析

2.1对于问题一的分析 问题一要求我们给出8种主要重金属元素在该城区的空间分布，这要求我们对题中所给数据进行分析后画图。首先要对表中数据进行分析，检查数据是否存在异常的情况。然后运用Matlab软件对数据进行插值，参考地理上用的等高线图，画出8种元素的等浓度线图，就可以看出8种重金属元素在城区的分布。分析该城区内不同区域重金属的污染程度，我们选择可以计算某一评价指标的单因子污染指数法和兼顾极值或称突出最大值的计权型多因子环境质量指数内梅罗综合污染指数法，通过计算，得出该城区内不同区域重金属的污染程度。 2.2 对于问题二的分析

问题二要求我们通过数据分析来说明重金属污染的主要原因，这就要求我们对题目中所给的数据进行深层的分析。运用因子分析法，首先，我们对该地区的重金属元素含量浓度的样本数据进行剔除特大值和特小值处理，再进行标准化处理，以消除不同变量的量纲的影响。然后，计算标准化后的数据的相关矩阵，求出相关系数的特征值和特征向量，以得到各个元素之间不同的相关性。使用因子分析法的关键就是利用相关系数矩阵求出相关因子的特征值和累计贡献率。为了使因子载荷两级分化，我们需要对原始数据矩阵进行正交变换，使用方差最大法。最后，通过确认因子个数，计算各因子得分进

2

行统计分析，结合实际情况，就能够说明重金属污染的主要原因。 2.3对于问题三的分析

根据重金属污染物的传播特征，我们可以得知8种重金属污染物的传播特征，选择其中一种扩散模型进行分析，建立污染源扩散模型，分析重金属污染物传播的特点与形态，建立相应的模型来确定污染源的位置。 2.4对于问题四的分析

由于问题三所建立模型在缺失数据的情况只考虑到一种传播情况，所得污染源位置是建立在最理想化的基础上的，因此必须获得必要的数据如当地水流、风速等情况才能进一步分析与检验该污染源的位置。在获得所需要的数据后，我们可以建立新的模型来 研究该地区的地质环境变化。

三、模型假设

1.假设题目中所给的数据都真实可靠；

2.假设该地区在近期以及未来一段时间都没有发生任何可能影响数据大小的灾害或人 为的改造，如地震或建造大型公共设施等；

3.把传染源传播污染物质的过程看作是在污染源所在的位置向四周等强度瞬时释放金 属污染物质，不计风力和大地的影响；

4.传播途径中，污染物本身是被动的、保守的，即污染物和空气无相对运动； 5.污染物浓度在各风向上分布为正态分布； 6.全部高度风速均匀稳定； 7.源强是连续均匀稳定的；

8.扩散中污染物的质量是守恒的（不考虑转化）。

四、符号说明与名词解释

4.1 符号说明 符号 i j Pi Ci Si P综 P平均 说明 角标，i=1，?，n 角标，i=1，?，n 区域重金属i的单项污染指数 第i种重金属浓度的实测值（μg/g或者ng/g） 第i种重金属浓度的评价起始值（μg/g或者ng/g） 综合污染指数（综合反映各污染物对区域土壤的不同作用） 所有单项污染指数的平均值 土壤环境中各单项污染指数中的最大值 第i个样本的第j个指标值 第j个指标的均值 第j个指标的标准差 空间点（x，y，z）的污染物浓度（ng/g） 源强，单位时间污染物排放量(ng/s) 平均风速(m/s) 烟气扩散系数，与大气稳定度和水平距离有关，并随x的增大而增加 Pmax xij xj ?j c q u ?y

3

4.2 名词解释

重金属污染[1]：重金属系指密度4.0以上约60种元素或密度在5.0以上的45种元素。砷、硒是非金属，但是它的毒性及某些性质与重金属相似，所以将砷、硒列入重金属污染物范围内。环境污染方面所指的重金属主要是指生物毒性显著的汞、镉、铅、铬以及类金属砷，还包括具有毒性的重金属锌、铜、钴、镍、锡、钒等污染物。

五、模型的建立与求解

5.1 问题一的求解

5.1.1 对数据源进行处理

对题目所给数据进行分析，发现数据中有个别几个数据与其他数据相比显得特别巨大，比如Hg（汞）元素在编号为9的点的数值为16000，明显比所测数据都大，但是通过分析我们认为，这几个点有可能为污染源或者垃圾掩埋点，所以并不把它们认定为异常数据。通过对整体数据分析后，我们发现除了上述几个特别的样点外，其他数据在正常范围之内。

5.1.2 8种主要重金属元素在该城区的空间分布

由上述，我们排除了数据源异常的情况，根据题目所给出的数据，我们利用Matlab 7.0 软件（程序代码见附录一）画出该城区的各功能区域的散点图（见图1）。

4x 10o-生活区 *-工业区 + 山区 ◇-交通区 星形-绿地区21.81.61.41.210.80.60.40.2000.511.522.5x 1034图1 该城区各功能区散点分布

接着，根据题目所给数据，画出8种重金属元素在该城区的等浓度分布图（程序代码见附录二）

4

18000160001400012000100001580006000400052000010202500.511.522.5x 104As等浓度分布图图2 As元素等浓度分布图

18000160001400012000100008008000600040002000060040020000.511.5Cd等浓度分布图22.5x 104

1600140012001000图3 Cd元素等浓度分布图

5

18000900160008001400070012000100008000600040002000060050040030020010000.511.5Cr等浓度分布图22.5x 104图4 Cr元素等浓度分布图

1800016000140001200015001000080006000400020000000.511.5Cu等浓度分布图22.5x 104

250020001000500图5 Cu元素等浓度分布图

6

180001600014000120001000080006000400020000160001400012000100008000600040002000000.511.5Hg等浓度分布图22.5x 104图6 Hg元素等浓度分布图

1800016000140001001200010000800060004000200004080140

120602000.511.5Ni等浓度分布图22.5x 104图7 Ni元素等浓度分布图

7

1800016000140001200010000800060004000450400350300250200150100200050000.511.5Pb等浓度分布图22.5x 104图8 Pb元素等浓度分布图

18000160003000140001200010000800060001000400050020000000.511.5Zn等浓度分布图22.5x 104

3500250020001500图9 Zn元素等浓度分布图

由上述八种重金属元素的等浓度图，我们可以很直观的看出这八种重金属元素在各

8

个区域的浓度高低的区间分布。

5.1.3 该城区各功能区的重金属污染程度分析

由于该城区的重金属污染包括8种重金属元素污染，因此我们考虑利用单项污染指数法来对各种元素污染程度进行评价，然后再利用内梅罗综合污染指数法综合评价该城区总体重金属的污染程度。 （1）单项污染指数法

为了分析该城区各种重金属的污染程度，我们利用如下单项污染指数法模型：

Cpi?i，

Si其中Pi为区域第i种重金属的单项污染指数；Ci为第i种重金属元素含量的实测值；Si为第i种重金属含量的评价起始值。若Pi>1,则表示该区域受到污染。

区域土壤背景值代表了未曾受到或者相对受到“三废”污染较少的土壤中有毒物质的平均值，具有显著的区域特点[4]。土壤总量污染起始值评价是以区域土壤中某污染物浓度的平均值加2倍标准差所得结果作为区域土壤污染的起始值[5]，并将其作为评价标准对各采集样点分析结果进行评价。其具体步骤是，利用题目给出的8种主要重金属元素的背景值数据（见表1）。

表1 重金属总量的土壤污染起始值计算结果 重金属含量平均值 标准差 含量起始值 元素 （μg/g其中Hg为ng/g） （μg/g其中Hg为ng/g） As 3.6 0.9 5.4 Cd 130 30 190 Cr 31 9 49 Cu 13.2 3.6 20.4 Hg 35 8 51 Ni 12.3 3.8 19.9 Pb 31 6 43 Zn 69 14 97 以平均值加2倍标准差作为Si，利用单项污染指数法模型计算得到该城区重金属元素总量的单项污染指数（见表2）。

表2 重金属元素总量的单项污染指数值统计结果 重金属样品数量/个 单项污染指数 超标率% 元素 最大值 最小值 平均值 As 319 5.57963 0.298148 1.05 49 Cd 319 8.525263 0.210526 1.59 65 Cr 319 18.79265 0.312653 1.09 34 Cu 319 123.9451 0.112255 2.70 70 Hg 319 313.7255 0.168039 5.88 48 Ni 319 7.160804 0.214573 0.87 22 Pb 319 10.98791 0.457674 1.44 55 Zn 319 38.77134 0.338763 2.07 57 由表2可知，该城区的土壤重金属元素都有超标，其中以Cu的含量超标最为严重，为70%,其他七种元素也超标，部分较为严重。

9

根据以上数据计算得到各个功能区重金属元素的单项污染指数（见表3）。

表3 各区单项重金属的单项污染指数值统计结果 超标率% 功能区 As Cd Cr Cu Hg Ni Pb Zn 319 319 319 319 319 319 319 319 样本个数 68.18 72.73 45.45 77.27 56.82 27.27 59.10 61.36 生活区 61.1 86.11 41.67 94.44 69.44 27.27 80.56 77.782 工业区 16.67 25.76 22.73 25.76 22.73 19.70 19.70 13.64 山区 49.28 76.81 36.23 82.61 52.17 19.57 66.67 71.40 交通区 公园绿地区 74.29 54.29 25.71 68.57 51.43 14.29 54.29 48.57 由上表可得，各区的各种重金属元素均超标，其超标率最高的元素统计如表4

表4 各个功能区超标率最高的元素 功能区 生活区 工业区 山区 交通区 公园绿地区 元素 Cu 77.27% Cu 94.44% Cd/Cu 25.76% Cu 82.61% As 74.29% 根据表4，得知该地区的Cu元素超标十分严重，严重影响到了当地的土壤。

（2）内梅罗综合污染指数法。

为了综合评价该城区总体重金属的污染程度，我们选择了内梅罗综合污染指数法[2]

进行分析。 内梅罗综合污染指数法的计算公式为：

?P平均?Pmax?? P综????2??其中P综为综合污染指数（综合反映各污染物对区域土壤的不同作用）；P平均为所有单项

2212污染指数的平均值；Pmax为土壤环境中各单项污染指数中的最大值。内梅罗指数[2]是一种兼顾极值或称突出最大值的计权型多因子环境质量指数，内梅罗综合指数的分级标准[3]

如表5：

表5 内梅罗综合污染指数的分级标准 P综 污染等级 污染程度 1 2 3 4 5

P综?0.7 0.7?P综?1.0 安全 警戒线 轻污染 中污染 重污染 清洁 尚清洁 超标 土壤、作物受中度污染 土壤、作物受严重污染 1.0?P综?2.0 2.0?P综?3.0 P综?3.0 10

表6为采用内梅罗综合污染指数法对区域土壤中重金属总量污染起始值统计计算后的结果。

表6 内罗梅综合污染指数法统计结果 评价分级 清洁、轻度、 尚清中度、评价标准 总样品清洁 尚清轻度污中度污重度污洁所重度污 数/个 /个 洁/个 染/个 染/个 染/个 占比染所占例% 比例% 重金属总 量污染起319 16 46 122 53 82 19.4 80.6 始值评价

由表6可知，采用重金属污染起始值评价时，轻度、中度和重度污染的样点数达到了80.6%，而清洁、尚清洁的样点的约占19.4%。

表7 受轻度、中度和重度污染区域在各个功能区的分布统计 功能区 1 2 3 4 5 所占比例% 14.4 14.0 13.6 47.5 10.5 由表7可见，在5个功能区中，受到污染最多的为交通区，其次为生活区、工业区、山区和公园绿地区。

综上，我们得知，该地区的污染问题十分严重，所测的8种重金属元素均超标，且超标率最高达到了70%，而47.5%的污染地区位于该地区的交通地带，其次是生活区、工业区、山区和公园绿地区。

5.2 问题二的求解 按照问题分析中对问题二的分析思路，借助SAS统计软件的insight模块运用因子

[6]

分析法对数据进行分析。在因子分析数学模型中，通过正交的方差最大旋转法使每一个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小。变量或因子的重要程度都是以其方差大小来衡量。因子旋转后每个变量因子负荷代表着在系统中的作用或重要性程度，以各个变量目标因子载荷平方与因子方差贡献率乘积作为变量的权重，构成一个判别污染来源的综合指标。 5.2.1数据标准化处理

对该地区重金属元素含量浓度的样本数据做剔除特大值与特小值后的处理，为了消除不同变量的量纲的影响，我们对数据进行了标准化处理，而且标准化转化不会改变变量的相关系数。 标准化公式为：

xij?xjxij?

?j其中xij为第i个样本的第j个指标值，而xj和?j分别为j指标的均值和标准差。 5.2.2 计算标准化数据的相关系数阵，求出相关系数的特征值和特征向量

首先，给出该城市表层土壤As、Cd、Cr、Cu、Hg、Ni、Pb和Zn八种重金属元素原始含量浓度数据的相关系数矩阵。如图10所示：

11

图10 变量相关矩阵

可见，Cr和Ni的相关性最好，相关系数最大，为0.7158，其次为Pb和Cd，相关系数为0.6603，以下依次是Cu和Cr，Pb和Cu的相关性较好，相关系数分别是0.5315和0.5200，Ni和Cu，Zn和Pb的相关系数分别为0.4946和0.4937。其他元素之间的相关性并不是很好。从成因上来分析，相关性较好的元素可能在成因和来源上有一定关联。 因子分析的关键就是利用相关系数矩阵求出相关的因子的特征值和累计贡献率，用SAS统计软件计算可得出，见图11：

图11 特征值和累计贡献率

在累积方差为93.16%（90%）的前提下，分析得到6个主因子，可以看出6个主因子提供了源资料93.16%的信息，满足因子分析的原则。从上图还可以看出因子1可能为该城市土壤重金属污染的最重要的污染源，对该城市重金属污染的贡献最大，因子2、因子3、因子4、因子5、因子6对该城市重金属污染有重要作用。 5.2.3 进行正交变换，使用方差最大化法 其目的是使因子载荷两级分化，而且旋转后的因子仍然正交。

图12 旋转前因子载荷矩阵

图13 方差极大正交旋转后因子载荷矩阵

图14 旋转后的因子载荷矩阵

12

5.2.4 确定因子个数，计算因子得分，进行统计分析

由以上图12、图13和图14可见，旋转前后因子载荷的变量结果基本一直。变量与某一个因子的联系系数绝对值（载荷）越大，则该因子与变量关系越近。正交因子解说明：因子1为Cr和Ni的组合，因子2为Cd和Pb的组合，因子3为Hg,因子4为As,因子5为Cu，因子6为Zn。Cr和Ni 、Cd和Pb可能是同一个来源，而且这两组元素正好是相关性最好的两组元素。 5.2.5结论

因子1为元素Cr和Ni的组合，结合第一题的重金属元素等浓度线图，我们可以得出以下结论，该城市表层土壤Cr和Ni基本没污染，只有个别点富集程度较高。这个可能是由于该富集中心位于工厂排污口或排气口附近。

因子2为元素Cd和Pb的组合，对两种元素的重金属等高浓度线图进行对比我们可以看出两种元素在空间分布上为集中于一个带状的污染带。且两种元素的富集中心主要坐落于交通区。这个主要因为交通区汽车等交通工具所排放的尾气污染物中均含有铅元素和铬元素等。所以我们可以认为因子2的两种重金属污染是有交通工具的尾气排放造成的。

因子3为元素Hg，从该元素的等浓度线图上分析，该元素的富集中心主要集中于交通区，这也是因为交通工具的尾气排放中的汞蒸汽所导致的。

因子4为元素As，从该元素的等浓度线图上分析，该元素的富集区中心为于山区，而砷是农作物农药当中的主要元素。所以我们可以认为砷污染的主要原因是由于山区进行农作物种植时农药的过度使用。

因子5为元素Cu,从该元素的等浓度线图上看，该元素的富集中心位于生活区与交通区附近，这可能是因为生活区的商业活动、居民生活和居民交通工具的使用所导致的土壤中的铜污染。

因子6为元素Zn,从该元素的等浓度线图上看，锌元素的富集中心位于工业区于交通区相间的区域，这可能是因为厂矿企业的三废排放和交通工具尾气的排放所导致的锌污染。

5.3 针对问题三的建模与求解

5.3.1 重金属污染物的传播特征分析

根据重金属污染物的传播特征，可以初步将所有传播特征分成五类：自然沉降和雨淋沉降（中心向四周扩散型）、含重金属元素物质的燃烧（带状分布有叠加性）、干湿沉降（吸附或固定使富集）、人为的不合理使用（过量形成富集）、其他，而结合题中所给的8种重金属污染物的特点可知，根据以上分类原则，一种重金属可能同时属于几种不同的传播特征类，因此我们先仅按其主要传播特征将其归类，考虑前四大类，可以得到下表8：

表8 8种主要金属污染物的传播特征 种类 传播特征 一 自然沉降和雨淋沉降（中心向四周扩散型） 二 含重金属元素物质的燃烧（带状分布有叠加性） 三 干湿沉降（吸附或固定使富集） 四 人为的不合理使用（过量形成富集） 因此根据上表的各类重金属污染物的传播特点，由第一类污染物是从污染源中心向四周扩散传播的，可以知道该类最适用于求解污染源，即只要找到污染物浓度最高的地方就是所要寻找的污染源。

5.3.1.1污染源扩散模型

13

扩散的数学模型基于一个基本假设[9]，即：在溶液中，溶质的扩散遵循Fick第一定律，用数学表达式表示为：

?C ， (5.1) F0??D0?X?C式中，F0为溶质在溶液中的扩散通量，D0表示溶质分子在溶液中的扩散系数，表示

?X的是浓度梯度，负号表示溶质分子沿着浓度降低的方向扩散。

在带有重金属元素的大气中，重金属元素的扩散同样遵循Fick第一定律，即：

?C ， (5.2) F??D?X式中，F表示的是元素在大气中的扩散通量，D表示的是元素在大气中的扩散系数。

根据质量守恒定律，可导出物质扩散的基本微分方程：

?C?Fx?Fy?Fz????0 ， (5.3) ?t?x?y?z式中，Fx，Fy，Fz为重金属元素在x，y，z方向的扩散通量。

当D为常数时，将Fx，Fy，Fz代入式(5.2)，得到污染源扩散的基本微分方程：

?C?2C?2C?2C?D(2?2?2) ， (5.4) ?t?x?y?z由于土壤溶质的扩散模式复杂多样，为了讨论的方便，我们不对其它因素加以考虑，例如溶质自重、边界条件、空间形状等，只对最基本的扩散情况（即瞬时点污染源）进行讨论。其他的情况大多可通过对其在时间和空间上的积分来求得。

设t =0时，在(x，y，z)处污染物质量为M，则其扩散方程为：

?C?2C?2C?2C?Dx2?Dy?Dz2 ， (5.5) 2?t?x?y?z若坐标原点在污染源发生处时，扩散方程的解析解为：

c(x,y,z,t)?M(12?Dxte?x24Dxt)(12?Dyt?ey24Dyt)(12?Dzte?z24Dzt) ， (5.6)

根据上式，当Dx?Dy?Dz时，扩散方程的解析解为：

c(x,y,z,t)?M328(?Dt)从式(5.7)中可看出，污染源扩散浓度与扩散系数D、坐标位置(x,y,z)以及污染物质量M、时间t密切相关，根据研究的需要，我们主要关注污染源浓度与空间位置之间的关系，对其他因素进行简化或采取赋值的方式进行处理，以突出研究的重点。

5.3.1.2重金属污染物的传播特征

为了进一步研究重金属污染物的扩散规律，给出时间t、扩散系数D、污染物质量M分别取某个静态点的图片：

x2?y2?z2exp?() (5.7)

4Dt 14

图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了解到：在扩散质量和时间相同的情况下，扩散系数越大，扩散范围越大，污染源浓度也越低。

由上述分析可以得出，污染源扩散浓度在二维空间上符合正态分布，其二维正态分布的联合概率密度函数为：

2??(x??12)2?(x??1)(y??2)(y??2)??11??f(x,y)?·exp?????? 2?222??2(1??)???2??1?21??1212?????

5.3.2 建立模型

在这里，我们选择建立高斯扩散模型[10]（即中心向四周扩散模型）。 1.高斯模式的有关假定

设定坐标系。坐标系取排放点（无界源、地面源或高架源排放点）在地面的投影点为原点，主风向为x轴，y轴在水平面内垂直于x轴，正方向在x轴的左侧，z轴垂直于水平面，向上为正，即右手坐标系。作如下假设：

1）把传染源传播污染物质的过程看作是在污染源所在的位置向四周等强度瞬时释放重金属污染物质，不计风力和大地的影响；

2）传播途径中，污染物本身是被动的、保守的，即污染物和空气无相对运动； 3）污染物浓度在各风向上分布为正态分布； 4）全部高度风速均匀稳定； 5）源强是连续均匀稳定的；

6）扩散中污染物的质量是守恒的（不考虑转化）。

15

2.无界空间连续点源扩散模式

由正态分布假定，得下风向任一点的浓度分布

c(x,y,z)?A(z)e?aye?bx，

22方差的表达式

2?y???0y2cdy??0cdy2，?x???0x2cdx?，

?0cdx 由假定条件，得出源强积分式

（单位时间物料守恒）q?????????ucdydx

其中浓度c为未知数， z=0时, a,b (=

)为待定系数,分别与x,y有关函数 2?2积分后，可以解出四个未知数：得到高斯模型

qy2x2c(x,y,0)?exp[?(2?)]， 22?y2?x2?u?y?x1其中，c表示空间点（x，y，0）的污染物浓度，ng/g;q表示源强，单位时间污染物排

放量，ng/s;u是平均风速，m/s；?y、?x为烟气扩散系数，与大气稳定度和水平距离有关。

结合本题所给数据，我们对高斯模型两边取以e为底的对数，得：

q1y2x2lnc(x,y,0)?ln?ln?y?ln?x?(2?2)，

2?y?x2?u

进而，我们可以推出关于污染源原点(x0,y0,z0)的模型：

1(y?y0)2(x?x0)2lnc?ln?ln?y?ln?x?(?) 222?y?x2?uq5.3.3 模型求解

2?u定这些参数，使得拟合值与真实值之间的误差达到最小。而常用的判断曲线拟合效果好坏的办法就是最小二乘法，因此我们运用最小二乘法来对曲线进行拟合。

根据各类重金属污染物的传播特点，只要找到污染物浓度最高的地方就是所要寻找的污染源。因此选择样本点时，应该选择在平均浓度较高，也就是在元素等浓度线图中浓度较高且集中的区域进行选点。

运用matlab7.0软件对所求参数进行估算（运行程序见附录三），从而得到传染源的位置如表9

表9各金属传染源坐标 重金属元素 As Cd/Pb Cr/Ni Cu 污染源点坐标 （12715，3042） (21349,25055) (3712,4996) (11357,39560) 污染源浓度 23.7299ng/g 1400ng/g 464.5288ng/g 748.9036ng/g 重金属元素 Hg Zn 污染源点坐标 （7167，2054） (12701,2843) 污染源浓度 3700ug/g 405.42ng/g

16

由于该模型中需要确认的参数有lnq，?y，?z，y0，x0，故采用曲线拟合来确

5.4 针对问题四的建模与求解

对于解决问题三时所建立的模型，该模型的优点是所需要的数据量相对较少，而缺点是该模型建立的前提是忽略了天气、地貌和水流等影响。所求的的结果是在相对理想化情况才存在的，在实际中借鉴意义不大。

地质环境是自然环境的一种，是指由岩石圈、水圈和大气圈组成的环境系统。在长期的地质历史演化的过程中，岩石圈和水圈、岩石圈和大气圈之间、大气圈和水圈之间进行物质迁移和能量转换，组成了一个相对平衡的开放系统。人类和其他生物依赖地质环境生存发展，同时，人类和其他生物又不断改变着地质环境。而地下水[7]作为水圈中极其重要的组成部分，由于水量稳定，水质好，一直是农业灌溉，工矿和城市的重要水源之一。随着经济的发展，人类生产活动中排出的有毒有害物质使得地下水质持续恶化。所以我们选择地下水污染作为研究城市地质环境演化模式的一个侧面，对地下水进行水质模拟，建立模型以便对地下水域内的地下水域状况进行了解。

为了对地下水水质进行模拟，我们需要采集该城市历年来的垃圾掩埋量，生活废水排放量，工业废水排放量，工业废气排放量，地下水总硬度超标率，各种重金属离子含量浓度超标率和浓度，

5.4.1 地下水污染成因分析

由于地下水水质受到气候，地下水位，污染源的存在以及互相影响非线性的关系。所以这里我们采用灰色关联度分析法来分析产生地下水污染的成因。 (1)建立参考数列和比较数列

设历年的地下水总硬度数据为参考数列yo，设垃圾掩埋量数据y1，生活废水排放量数据y2，工业废水排放量数据y3和工业废气排放量数据y4作为比较数列。 (2)对各类列进行初值化处理

yj?yj(k)maxyj(k)(y=1,2,3,4)

(3)计算参考数列在k点的灰色关联系数：

minminy0(k)?yi(k)??*maxmaxy0(k)?yi(k)ik(i=1,2,3,4)， ?i(k)?iky0(k)?yi(k)??*maxmaxy0(k)?yi(k)ik?为分辨系数，具体取值规则如下[11]：

记?v为所有差值绝对值的均值，即

1mn?v???y0(k)?yi(k)

m*ni?1k?1

记????max?v，则?的取值为：当?m?max?3?v时，1.5?????2??。

ax>3?v时，?????1.5??；

该模型侧重于总体分析，?不仅可以调节?i(k)的大小，而且可以调节它的变化区间，?i(k)的下界值随?增大而增大，下界值增大说明区间变小，分辨率变低，分辨效果则不明显。根据因素见的关联分析可选择不同的分辨系数，在计算时，一般取?=0.5便可得到满意的分辨率。当???则无法进行关联分析，此时所有?i(k)?1，即关联系数蜕变为一个点[8]。 (4)计算灰色关联度

17

1ri?N??(k)(i=1,2,3,4)

ik?1n (5)排关联序

在得出各个序列的灰色关联度排序后我们就可以比较出城市地下水总硬度这一指标与垃圾掩埋量， 生活废水排放量，工业废水排放量和工业废气排放量中那一个指标关联系数较高。从而得出地下水总硬度超标率的主要原因。依照上述(1)-(5)的做法我们同样可以用灰色关联度分析法对造成城市地下水各种重金属离子超标的成因做出判断。 5.4.2 地下水污染预测

地下水中的污染是一个包含了物理作用，化学作用和生物作用的过程。由于作用复杂，我们可以利用灰色系统理论，把地下水环境视为灰色系统，他污染物看作是在一定范围内变化的灰色变量。从而我们可以建立起GM（1，1）模型[12]。 （1）GM（1，1）模型建立原理。

设原始数列为Q(0)?{Q(0)(1),Q(0)(2),????,Q(0)(n)}为原始非负时间序列，预测的结果就是想知道Q(0)(n?1),Q(0)(n?2),????

Q(t)为累加声称序列，即Q(t)=?Q(0)(i?1),i?2,3????n

(0)(1)tm?1具体地说，就是：

(1)(0)??Q(1)?Q(1) ?(1)(0)(i)??Q(i)?Q(i)?Q(i?1),i?2,????,nGM(1,1)模型的白化微分方程为：

dQ(1)?aQ(1)?u dt??a?a为待辨识参数；u为待辨识内生变量；设待辨识变量a???u??，按最小二乘法得：

??a?(BTB)?1BTy，

??1(1)?(1)-(Q(1)?Q(2))1?Q(0)(2)??2??(0)??1?(1)(1)Q(3)?其中B???2(Q(2)?Q(3))1?，y=?。

????????????1????(0)1?-(Q(1)(n?1)?Q(1)(n))1???Q(n)?????2? 于是我们可以得到灰色预测模型为：

uuQ(1)(t?1)?(Q(0)(1)?)e?at?,t?1,2,3????n

za (2) 残差检验和后残差检验：

模型建立后，我们还应该对模型进行残差检验和后验差检验。残差检验有两种，绝对误差和相对误差检验： 绝对误差

? 相对误差

(0)(i)?Q(i)?Q(i),(i?1,2,????n)，

(0)?(0) 18

41 42 43 44 45 46 377.14

16.58 459.00 94.79 47.17 1900.00 19.90 71.32 215.10 7.41 337.30 77.27 248.85 90.00 20.10 99.58 210.00 5.93 568.10 75.14 118.16 135.00 23.80 111.54 572.96

4.69 599.00 69.05 122.18 121.00 19.80 102.72 427.04

4.90 635.50 68.42 227.76 176.00 19.50 96.33 538.98 5.31 600.70 44.65 45.10 51.00 15.50 65.87 186.33 47 4.29 567.60 48 5.51 228.50 49 4.69 568.60 50 7.20 214.70 51 5.31 151.90 52 4.90 343.30 53 4.90 293.90 54 3.89 312.90 55 3.69 315.90 56 3.11 416.30 57 3.89 374.00 58 3.89 344.30 59 2.91 252.90 60 3.30 503.40 61 4.90 303.80 62 4.09 127.00 63 2.91 265.00 64 2.72 278.90 65 3.11 751.20 66 3.30 361.30 67 3.30 488.00 68 6.14 227.00 69 3.69 347.40 70 4.49 136.00 71 3.11 327.10 72 8.06 113.10 73 3.69 270.50 74 3.69 160.30 75 3.50 305.50 76 2.72 70.90 19.45 9.12 77 8.50 261.30 78 1.77 119.80 79 2.53 468.80 80 3.69 150.70 81

6.14 100.30

60.25 48.67 46.00 16.10 63.74 208.06 49.27 30.85 62.00 22.90 45.93 102.04

306.02 70.41 900.00 16.80 79.67 196.73 50.33 40.16 156.00 20.80 47.76 403.98 47.24 24.44 140.00 17.30 37.49 92.55 42.01 58.81 80.00 13.80 79.07 275.82 60.29 51.03 53.00 12.60 75.93 278.37 33.79 277.82 55.00 14.00 68.24 295.61 45.43 34.05 55.00 12.60 62.84 196.33

57.88 47.64 167.00 11.90 116.19 242.04 45.17 50.19 35.00 15.00 58.11 157.35

35.29 47.87 100.00 15.10 133.72 141.02 45.98 71.54 32.14 14.40 42.99 146.22 38.74 30.46 36.43 7.20 53.73 102.86 56.02 65.86 63.21 40.05 90.69 3760.82 27.58 23.99 30.00 11.93 57.47 85.61 35.66 29.39 24.64 9.23 60.54 122.96 43.43 32.61 64.29 9.90 53.40 135.71 53.11 53.80 27.86 10.46 60.27 155.00 47.54 52.28 25.71 9.11 113.46 218.27 51.18 34.55 37.50 10.80 54.62 125.92 42.15 67.04 49.29 16.31 34.28 82.96 37.76 19.97 26.79 10.01 54.41 221.22 36.56 23.07 21.43 14.96 34.19 78.98 25.98 23.73 25.71 9.79 63.81 138.06 52.40 20.81 65.36 19.69 29.56 62.24 33.12 57.85 25.71 13.50 62.04 118.16 38.29 26.08 25.71 14.29 40.13 82.86 39.50 30.86 31.07 14.74 61.89 148.88 15.00 7.09 22.73 32.86 45.97 24.83 51.43 18.00 43.39 84.49 15.32 13.34 8.57 6.19 26.31 47.76 37.04 32.03 45.00 12.15 65.25 178.98 59.61 19.00 34.29 24.98 38.47 89.08 37.49 20.23 34.29 14.85 29.29 61.94

24

82 83 84 85 86 87 88 89 90 10.99 109.80 56.07 69.06 58.93 20.70 38.87 63.27 6.35 91.80 36.12 16.91 36.43 12.49 27.01 47.76

30.13 743.90 49.03 26.18 27.86 17.66 72.76 182.04 3.89 416.80 37.04 23.78 22.50 11.48 54.45 105.00 2.91 369.80 36.34 52.48 22.50 7.99 42.02 84.08 1.96 194.00 18.08 16.17 26.79 6.98 40.27 94.69 6.98 50.10 41.02 14.25 17.14 13.39 26.57 40.92 2.91 198.80 28.21 19.24 13.93 9.56 47.81 94.80 5.93 886.60 42.69 28.12 43.93 21.15 94.64 163.27 91 5.93 128.90 92 7.41 114.30 93 4.29 232.90 94 4.69 132.80 95 5.72 1619.80 43.48 96 6.77 282.50 97 4.49 180.60 98 3.11 386.60 99 2.91 345.00 100 4.29 95.60 22.49 17.15 85.71 10.13 27.97 67.24 101 7.63 87.10 45.83 14.83 30.00 14.63 29.25 48.78 102 5.93 203.00 103 2.34 353.00 104 2.91 233.20 105 5.72 174.30 106 2.34 87.60 18.46 9.73 107 6.56 245.60 108 4.69 167.90 109 6.35 111.40 110 5.10 94.60 77.92 20.34 23.57 28.69 25.92 58.47 111 4.69 111.60 112 3.50 85.50 26.33 24.88 36.43 11.36 26.77 64.08 113 4.69 169.50 114 4.49 138.00 115 3.30 131.40 116 2.91 41.00 41.77 12.50 17.14 17.10 19.68 78.37 117 4.09 129.70 118 5.72 148.00 119 4.90 108.10 120 4.90 132.50 121 4.09 162.90 122 5.93 88.80 52.41 15.30 10.71 19.91 28.62 63.88 123 2.91 206.70 124 2.72 121.80 125

2.34 231.10

47.52 16.31 12.86 17.66 33.51 91.73 48.34 21.45 35.36 16.54 35.83 63.88 29.17 40.02 1714.29 9.79 38.65 95.41 36.11 17.28 20.36 15.53 37.03 82.76 15.50 20.36 15.41 30.99 57.55 41.97 52.80 27.86 18.34 49.10 104.90 37.23 18.70 27.86 11.93 36.45 63.98 35.93 26.38 24.64 12.26 60.00 157.24 40.46 152.21 23.57 15.53 58.05 170.71 35.97 16.88 15.00 14.51 45.83 89.49 24.53 12.70 11.79 9.00 58.80 89.08 24.92 21.62 85.71 8.33 45.20 100.10 33.83 29.45 20.36 13.05 42.10 71.43 13.93 8.89 24.43 43.37

36.73 61.30 55.71 14.18 47.24 114.29 33.15 18.96 60.00 15.98 33.46 55.71 28.82 59.17 206.79 11.70 28.02 61.53 24.57 12.09 31.07 8.55 27.12 43.88 39.11 22.51 25.71 15.98 39.53 82.96 34.52 35.54 50.36 12.60 25.45 52.45 35.97 11.29 43.93 11.03 30.06 61.94 26.83 10.12 40.71 7.76 28.84 68.16 36.73 14.21 52.50 12.60 27.15 57.86 22.72 9.40 35.36 8.89 32.39 69.39 79.52 18.67 42.86 27.68 28.30 92.24 33.79 12.96 30.00 12.94 29.17 74.90 17.46 12.02 31.07 5.51 53.79 79.18 19.98 7.71 24.64 8.10 27.30 57.96 36.48 12.80 37.50 17.66 29.02 96.53

25

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

2.53 128.60 21.75 9.10 30.00 9.56 22.77 58.06 3.89 126.80 50.01 12.74 53.57 18.68 29.40 69.90 2.34 79.70 74.36 27.49 9.64 45.45 21.22 87.65 2.34 294.10 42.68 12.27 79.29 15.53 83.44 99.49 1.96 238.30 19.71 25.22 19.29 6.64 66.24 142.86 2.72 109.30 18.52 12.55 19.29 9.79 26.06 61.43 3.11 69.40 22.93 7.56 26.79 7.65 22.09 53.06 3.50 155.70 33.30 12.08 42.86 11.36 28.05 51.53 2.72 109.30 19.83 14.54 19.29 9.45 34.56 47.55

1.77 205.50 173.34 27.81 41.79 74.03 55.71 229.80 3.50 53.20 35.15 7.91 20.36 11.14 24.90 34.39 2.53 40.00 19.08 2.29 15.00 7.88 24.15 54.29 9.17 82.80 63.88 24.85 11.79 31.50 21.68 76.53

7.41 222.40 31.10 38.50 142.50 8.66 128.60 133.98 5.72 164.80 49.44 27.98 126.43 13.61 37.64 76.12 8.06 194.00 41.79 31.43 213.21 13.61 45.52 127.35 5.72 310.10 56.40 37.11 138.21 16.43 59.63 154.49

8.72 1024.90 77.61 71.24 128.57 19.35 227.40 1389.39 6.77 492.80 37.86 35.61 142.50 13.39 114.13 118.98 6.56 202.30 33.28 37.67 1339.29 11.36 140.62 110.92 6.77 97.20 38.67 18.12 24.64 14.74 28.77 49.90 6.98 208.60 41.00 28.55 67.50 15.19 93.83 84.59 6.14 143.40 41.67 22.57 88.93 12.49 38.14 65.61 6.69 158.10 36.61 16.22 36.00 12.91 31.67 63.57

8.00 332.60 77.06 46.01 240.00 24.47 68.10 259.29 6.69 401.80 84.94 60.11 138.00 19.16 96.76 227.55 8.23 756.40 42.73 87.52 63.00 19.26 88.74 184.69 9.35 407.50 55.54 61.83 112.00 24.05 66.82 208.27 8.90 307.30 54.39 57.21 326.00 25.72 131.93 256.94 3.77 242.10 30.93 32.13 28.00 11.56 50.60 144.69 5.41 178.90 29.54 23.73 52.00 9.89 49.84 118.88 7.78 315.50 49.76 28.03 550.00 18.95 45.73 109.29 5.62 134.60 25.33 19.10 45.00 11.66 40.50 87.14 5.41 235.60 36.88 48.80 43.00 14.06 53.61 213.47 4.58 203.80 39.03 24.18 87.00 16.66 53.09 138.88

6.91 568.50 54.59 113.46 264.00 23.22 82.40 399.90 5.00 506.50 59.45 70.71 202.00 26.13 78.01 334.39 5.62 880.00 78.29 121.12 293.00 25.61 171.14 540.00

6.91 250.30 39.09 50.46 81.00 20.41 59.17 172.65 7.78 161.90 42.08 77.92 186.00 18.22 62.81 171.84 6.26 249.40 54.70 81.74 75.00 25.51 57.92 171.12 7.56 248.20 42.23 69.39 57.00 16.56 44.54 165.10 4.79 156.30 41.16 32.40 83.00 14.58 35.50 103.88

26

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

2.77 120.30 49.85 18.38 43.00 25.51 37.93 115.41 6.26 429.20 54.47 56.60 45.00 19.78 87.50 230.41 7.34 205.10 43.63 23.51 47.00 19.78 41.33 104.29 4.17 113.80 37.27 17.06 34.00 14.47 26.74 64.49 5.00 221.30 30.75 70.32 47.00 16.14 47.20 185.41 7.56 580.10 39.78 85.46 270.00 17.49 66.69 315.92 5.62 171.00 31.23 25.73 292.00 15.72 31.09 110.61 6.05 365.00 35.92 30.91 110.00 17.91 44.75 147.35 4.79 289.70 42.12 41.10 57.00 17.49 76.80 237.45

23.72 452.10 37.22 43.61 160.00 16.35 73.27 1635.92 6.47 847.60 53.17 59.00 96.00 19.05 68.82 186.53 4.17 310.70 40.70 42.64 58.00 14.99 115.59 177.76 5.00 317.50 42.46 30.96 162.00 17.81 59.94 221.94 4.17 357.50 38.78 43.89 14901.00 17.70 98.22 349.80 6.26 387.20 38.03 30.06 85.00 21.97 54.98 142.55 5.00 196.50 50.03 18.56 29.00 25.82 33.90 84.80 1.61 295.80 15.40 18.35 19.00 4.27 40.42 106.53 4.58 129.00 31.09 18.93 38.00 15.10 29.76 69.80 5.41 204.90 40.16 25.86 41.00 16.24 33.28 80.00 3.57 174.80 31.13 26.21 26.00 11.14 31.77 94.08 2.00 197.00 78.36 22.98 19.00 24.05 43.73 104.18 3.37 106.70 16.31 9.04 17.00 7.60 26.89 37.14

6.47 738.00 47.35 24.07 135.00 19.58 34.09 101.73 3.37 788.70 28.90 17.62 34.00 15.51 26.04 91.02 6.91 133.10 40.58 32.64 27.00 22.39 33.97 115.20 2.38 149.50 27.18 11.19 54.00 14.47 38.35 59.49 7.56 206.60 55.79 24.44 22.00 28.63 29.73 81.73 8.67 78.90 58.87 21.65 26.00 23.43 26.52 70.51

4.79 178.60 32.31 29.01 45.00 14.26 47.98 104.90 5.41 250.10 34.32 16.25 15.00 18.53 41.29 90.20 5.83 89.20 54.90 23.28 19.00 21.76 24.95 53.57 5.62 453.50 37.22 21.92 40.00 16.56 41.51 87.35 4.38 92.40 38.70 15.03 21.00 13.85 24.22 47.24 2.77 198.70 18.11 15.33 19.00 9.68 34.15 95.92 5.41 260.10 43.81 26.46 33.00 13.22 39.49 87.04 5.62 101.10 49.01 19.96 23.00 20.82 26.00 71.53 6.69 146.50 43.39 19.89 36.00 17.49 33.24 60.00 7.56 66.50 49.24 14.67 19.00 17.08 25.48 41.12 2.77 119.00 27.49 10.62 44.00 10.52 31.94 53.27 1.61 434.80 21.29 12.84 18.00 8.85 56.14 106.43 5.62 188.00 30.52 26.03 58.00 14.26 31.96 77.45 2.97 77.50 16.91 9.95 34.00 7.50 28.15 41.63 4.58 130.80 35.61 19.59 23.00 16.14 29.87 61.12 5.20 134.50 29.37 16.33 32.00 13.54 29.83 59.90

27

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

3.17 101.80 32.23 10.68 12.00 12.81 33.20 51.12 5.00 232.00 32.43 24.23 37.00 15.10 36.08 101.84 7.56 121.10 62.91 31.03 44.00 27.28 28.36 88.16 5.20 238.70 47.18 59.66 24.00 19.68 62.57 108.47 6.05 151.20 92.76 31.43 79.00 32.07 31.57 130.61 2.77 202.60 56.54 18.09 49.00 28.11 26.14 82.14

2.38 351.60 64.86 28.97 65.00 17.28 112.62 105.31 5.20 173.80 58.28 19.99 51.00 22.49 31.61 69.80 8.67 423.30 61.41 117.83 216.00 30.30 178.88 293.16

5.41 645.30 50.27 80.60 113.00 23.32 99.42 358.27 7.12 967.70 36.41 70.71 66.00 16.87 119.35 457.96 4.58 445.70 28.32 132.05 139.00 12.91 87.44 303.06 8.23 420.60 35.81 73.86 561.00 22.49 60.91 241.63 4.58 257.40 42.25 44.84 260.00 19.58 46.90 192.65 6.05 248.70 44.36 42.79 150.00 19.89 72.20 170.10 8.23 350.30 59.07 36.58 177.00 24.57 52.56 170.00 4.17 109.90 29.09 13.71 31.00 11.35 25.08 51.63

3.97 306.20 38.30 36.22 115.00 13.54 60.63 145.51 10.74 306.50 53.98 97.05 255.00 27.91 69.40 177.86 3.77 515.60 28.74 77.47 1801.00 15.10 70.18 178.57 10.27 175.90 57.88 37.50 65.00 26.97 45.34 100.41 5.20 417.80 35.71 47.58 801.00 16.03 94.93 204.08 6.47 368.20 44.22 65.70 91.00 18.33 76.88 247.04 6.47 393.00 35.82 70.17 135.00 18.53 72.96 254.80 8.90 117.50 32.54 20.51 41.00 14.47 41.45 56.33 3.37 379.10 32.13 37.81 49.00 13.74 70.37 208.37 6.69 320.30 87.90 66.99 78.00 21.90 123.60 230.11 4.17 613.40 50.34 149.26 133.00 16.60 110.20 447.83

8.23 220.50 61.48 41.96 84.00 22.40 58.37 153.27 4.38 229.30 32.69 60.40 64.00 13.70 62.62 204.81 5.20 334.30 47.05 40.37 76.00 18.40 61.03 211.47

5.20 284.30 40.47 44.73 349.00 16.00 119.46 184.61 5.41 230.10 41.65 35.71 57.00 16.90 43.91 118.70 5.00 551.10 71.75 71.23 74.00 22.30 111.62 209.50 9.58 201.00 39.11 25.26 125.00 10.70 99.35 90.07 3.17 203.80 26.57 24.29 23.00 11.90 41.69 94.34 3.77 355.20 26.25 32.87 14.00 9.60 59.22 117.87 9.13 133.70 52.99 22.10 26.00 24.80 37.73 95.48 7.34 95.60 47.85 19.51 14.00 20.80 28.52 57.68

5.62 352.60 44.57 58.98 51.00 13.30 69.95 531.45

4.79 556.20 50.87 143.31 92.00 19.10 180.05 388.69 5.20 113.70 41.85 20.73 12.00 19.00 31.87 57.27

28

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

3.37 240.50 28.04 22.63 261.00 11.70 35.74 92.46 7.34 120.30 54.19 21.90 27.00 23.00 29.63 81.01

3.57 514.10 42.34 47.67 13800.00 17.30 69.96 269.89 4.38 246.30 29.90 21.84 84.00 14.10 59.00 95.17 5.41 158.00 46.86 24.02 31.00 19.30 36.27 79.24 4.38 211.30 27.79 19.01 34.00 12.00 38.38 81.42

5.62 236.50 35.95 66.52 199.00 13.90 40.98 193.67 6.05 193.00 40.60 24.88 27.00 14.40 33.53 84.86 6.26 169.70 44.26 88.12 46.00 17.20 42.71 97.67 5.20 320.10 35.92 36.86 68.00 16.50 58.46 162.85 6.91 180.20 54.08 27.01 37.00 18.40 44.13 118.91 4.58 351.80 55.39 78.07 87.00 16.90 69.55 188.88 8.67 245.70 47.79 27.55 35.00 18.40 53.42 98.81 6.47 86.80 41.12 15.46 23.00 15.90 37.53 70.18

7.12 367.80 92.02 49.80 97.00 16.30 41.26 321.12 3.77 467.10 49.03 34.44 45.00 15.40 60.83 132.86 5.41 364.70 40.34 40.93 79.00 18.70 83.32 175.34 10.97 248.50 40.61 61.52 81.00 17.20 76.19 168.05 9.81 171.80 75.38 163.20 30.00 26.30 45.27 125.16 8.23 409.90 44.67 66.92 80.00 36.00 96.85 197.63 5.41 302.50 34.22 27.60 408.00 14.80 68.70 218.24 2.77 236.20 42.67 16.35 62.00 9.40 41.88 149.52 7.78 114.50 56.38 26.96 36.00 22.40 31.24 75.91 6.47 165.20 73.40 42.73 40.00 19.70 84.13 95.69 5.62 380.40 46.63 28.31 48.00 14.60 83.82 155.98 3.77 398.40 29.57 18.64 60.00 10.50 113.84 172.53 3.57 268.60 28.11 23.20 64.00 12.20 54.52 101.00 4.38 126.50 28.57 20.57 19.00 12.10 25.17 53.10 6.91 290.30 47.87 28.90 34.00 14.80 44.26 94.96 7.12 228.40 40.29 25.15 37.00 15.30 40.10 83.19 4.38 305.50 52.44 22.92 13.00 19.10 45.21 109.33 3.97 407.60 35.65 22.33 11.00 18.60 60.36 121.62 8.00 96.60 19.42 11.26 12.00 7.50 27.54 47.38 3.57 185.50 23.15 13.42 34.00 9.50 29.49 92.36 3.37 288.70 26.12 15.10 18.00 10.30 30.14 63.83 3.17 90.40 16.20 8.30 32.00 7.00 44.31 44.67 2.97 285.40 26.86 15.00 65.00 13.30 38.63 68.20 7.12 100.00 46.15 19.43 14.00 23.00 22.01 65.70 5.62 306.90 42.02 25.21 40.00 19.20 36.99 141.50 5.83 319.50 43.32 25.89 54.00 15.80 40.98 83.40 3.97 100.10 21.69 11.96 40.00 7.90 42.79 69.24 3.17 218.60 39.51 15.26 34.00 10.50 50.98 84.24 2.77 239.80 26.06 15.82 35.00 10.20 41.43 80.90 3.17 156.80 19.73 8.34 42.00 7.60 39.21 71.01

29

299 2.97 281.10 28.56 11.42 48.00 12.60 37.95 81.01 300 3.17 142.50 36.75 9.93 43.00 13.30 32.61 61.64 301 1.80 195.50 28.53 7.32 36.00 9.70 37.41 50.92 302 3.17 153.90 20.90 8.21 37.00 7.60 31.37 38.21 303 3.77 104.20 30.34 12.34 24.00 11.80 39.31 57.16 304 4.79 72.10 65.54 11.55 35.00 19.80 26.04 47.58 305 3.57 190.80 31.33 10.67 65.00 15.70 51.56 94.02 306 6.47 282.90 52.68 20.34 25.00 22.90 32.53 103.50 307 7.34 149.00 44.22 20.14 33.00 16.00 35.43 147.75 308 8.23 121.30 43.29 31.63 86.00 11.40 33.21 46.86 309 10.74 479.20 96.28 29.23 98.00 25.30 80.36 112.35 310 11.68 870.50 70.84 35.17 302.00 29.10 78.15 435.44 311 7.34 279.00 51.25 27.95 44.00 22.50 51.20 117.66 312 6.05 162.00 36.22 17.91 35.00 14.20 36.41 61.02 313 5.41 907.00 43.08 36.48 10.00 14.50 41.02 121.20 314 6.26 132.90 42.59 16.58 27.00 16.20 35.52 63.31 315 6.47 197.00 38.18 21.09 64.00 18.60 40.18 168.05 316 6.47 100.70 36.19 13.31 42.00 11.50 34.34 56.23 317 4.79 119.10 35.76 19.71 44.00 9.90 39.66 67.06 318 7.56 63.50 33.65 21.90 60.00 12.50 41.29 60.50 319 9.35 156.00 57.36 31.06 59.00 25.80 51.03 95.90 ];

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);

30

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]; %3??μ

option=optimset('MaxFunEvals',10000); [x,fval]=lsqcurvefit(@qq,x0,xdata,ydata) qq(x,x(4:5))

31

32

本文来源：https://www.bwwdw.com/article/u957.html