2011数学建模A题-城市表层土壤重金属污染分析.doc

城市表层土壤重金属污染分析
摘要
对于问题一，首先采用克里格插值法根据附件中给的采样点的浓度数据对城区内每种重金属浓度值进行插值，进而绘制每种重金属的浓度的空间分布图，对此为了方便计算我们借助于Surfer软件进行绘图，然后我们采用单因子指数法和内梅罗综合指数法对各功能区的污染程度进行评价，但是这两种方法无法从自然异常中分离出人为异常，为了弥补其不足，采用地累积指数法[1]对土壤重金属污染做进一步评价，继而得到不同区域重金属的污染程度。评价结果如表9所示。
对于问题二采用主成分分析法，旨在利用降维的思想，把多指标转化为少数几个综合指标， 根据主成分分析法的一般步骤，首先对附件中给出的重金属浓度进行标准化处理，然后得到各金属之间的相关系数矩阵，求出相关系数矩阵的特征值和向量值，再得到因子的成分矩阵，确定出主成分的个数，计算出各因子的成分得分矩阵，最后通过算综合主成分中各个因子的权重系数得到污染性较大的因子，最后分析该污染主要原因。结果如表14和评价结果。
问题三我们用对流-扩散偏微分方程来进行描述，对流扩散方程是描述粘性流体运动的非线性偏微分方程模型，我们将对流扩散方程进行简化，即变为二阶椭圆形偏微分方程，利用有限插值数值法估计出污染物的浓度分布，并与实际检测值相比较，偏差较小的即为污染源的位置及源强，即将问题转化为非线性最优化问题，结果如表所示。
对于问题四在分析问题三模型的优缺点基础上，为更好的研究城市地质环境的演变模式，将问题三的模型中二维对流扩散偏微分方程扩展到三维 ，这样就可以全面考虑影响模型准确性的各相关参数，通过收集不同的地理，天气条件下地质元素的空间分布信息，进而根据这些信息，建立优化模型，即三维方程的参数，在此基础上通过模拟仿真进行分析。

关键字：单因子指数法 内梅罗综合指数法 地累积指数法 主成分分析法 偏微分模型

一 问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0-10厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二 问题分析

对于问题一，要根据问题附件中给出的采样点以及在采样点处8种主要重金属元素的信息，研究各金属元素的空间分布，并分析不同功能区的污染程度，附件中只给出了金属元素在采样点处的浓度，要给出在整个城区内的空间分布，需要通过插值得到更加密集的浓度分布值，由于城区内样本点间的重金属浓度存在空间相关性，因此可以通过克立格法进行插值。在这里我们借助Surfer8.0软件，也是通过克里格插值的方法生成重金属含量数据，并绘制等值线图。由于还需要分析不同区域重金属的污染程度，污染程度需要抽象为污染指数，通过比较污染指数与国家背景值，确定污染程度较大的金属元素，针对这个问题我们采用单因子指数法和内梅罗综合指数法[1]，但是这两种方法无法从自然异常中分离出人为异常，为了弥补其不足，采用地累积指数法[3]对土壤重金属污染做进一步评价，继而得到不同区域重金属的污染程度。
对于问题二要求通过数据说明重金属污染的主要原因，现有资料表明，某些重金属空间分布具有相关性，相关性较大的重金属可能在成因和来源上有一定的联系。在此我们选用主成分分析法来说明。主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指标。在实证问题研究中，为了全面、系统地分析问题，我们必须考虑众多影响因素。这些涉及的因素一般称为指标，在多元统计分析中也称为变量。因为每个变量都在不同程度上反映了所研究问题的某些信息，并且指标之间彼此有一定的相关性，因而所得的统计数据反映的信息在一定程度上有重叠。在用统计方法研究多变量问题时，变量太多会增加计算量和增加分析问题的复杂性，人们希望在进行定量分析的过程中，涉及的变量较少，得到的信息量较多。
根据主成分分析法的一般步骤，首先对附件中给出的重金属浓度进行标准化处理，然后得到各金属之间的相关系数矩阵，求出相关系数矩阵的特征值和向量值，再得到因子的成分矩阵，确定出主成分的个数，计算出各因子的成分得分矩阵，最后通过算综合主成分中各个因子的权重系数得到污染性较大的因子，最后分析该污染主要原因。
对于问题三要求分析重金属污染物的传播特征，建立模型并分析污染源的位置，对此我们用对流-扩散偏微分方程来进行描述，对流扩散方程是描述粘性流体运动的非线性偏微分方程模型，可以刻画很多自然现象，如污染物的扩散、降解、流体流动与传热等，我们将对流扩散方程进行简化，即变为二阶椭圆形偏微分方程，利用有限插值数值法估计出污染物的浓度分布，并与实际检测值相比较，偏差较小的即为污染源的位置及源强，即将问题转化为非线性最优化问题。

对于问题四在分析问题三模型的优缺点基础上，为更好的研究城市地质环境的演变模式，应该确定收集什么信息，有了这些信息如何建立模型来解决问题。首先将问题三的模型中二维对流扩散偏微分方程扩展到三维 ，这样就可以全面考虑影响模型准确性的各相关参数，为了得到参数的准确表达式，通过收集不同的地理，天气条件下地质元素的空间分布信息，进而根据这些信息，建立优化模型，即三维方程的参数，在此基础上通过模拟仿真进行分析。
三 问题假设
1.问题一中国家评价标准采用二级标准。
2.假设全部采样点在自然情况下产生，排除其他因素的干扰。
3.假设该城区没有受到灾难性的污染。
4.观测值能够反映一段时间内保持稳定的浓度分布。
5.假设污染源是均匀向外扩散的。
四 符号说明
污染物类别
土壤中污染物实测浓度
土壤中污染物的背景值
土壤单因子污染指数
为综合污染指数
元素污染指数平均值
元素污染指数中的最大值
地累积指数
各个区域元素的平均值

元素地球化学背景值
各地岩石差异导致的背景值变动系数
第个采样点的第种重金属的浓度值
第种重金属的均值
第种重金属的标准差
二元素相关系数
五 建模与求解
5.1问题一建模与求解
为了得到更加形象具体的重金属空间分布图，在这里我们借助Surfer8.0软件，也是通过克里格插值的方法生成重金属含量数据，并绘制等值线图。

图1 城区地形三维图

图2各功能区地形图

图3 城区浓度分布图

图4 城区浓度分布图

图5 城区浓度分布图

图6 城区浓度分布图

图7 城区浓度分布图

图8 城区浓度分布图

图9 城区浓度分布图

图10 城区浓度分布图
5.1.1单因子指数法
我们对该城区不同区域的污染程度进行评价，由于每个区域都存在八种重金属，所以我们对每个区域的八种重金属含量分别进行评价，评价模式分为单因子指数法和内梅罗综合指数法。首先建立单因子污染指数式为：
（1）
式中为土壤中污染物i的污染分指数，为土壤中污染物实测浓度，为土壤中污染物的背景值。然后通过求解上述单项污染指数式得到五个区域种污染元素的分指数，如表1（表2、表3、表4、表5）所示为生活区污染分指数。
表1 各区域重金属实测平均浓度

编号
As(μg/g)
Cd(ng/g)
Cr(μg/g)
Cu(μg/g)
Hg(ng/g)
Ni(μg/g)
Pb(μg/g)
Zn(μg/g)
1
6.27
289.96
69.02
49.40
93.04
18.34
69.11
237.01
2
7.25
393.11
53.41
127.54
642.36
19.81
93.04
277.93
3
4.04
152.32
38.96
17.32
40.96
15.45
36.56
73.29
4
5.71
360.01
58.05
62.21
446.82
17.62
63.53
242.85
5
6.26
280.54
43.64
30.19
114.99
15.29
60.71
154.24

表2 各区域单因子污染分指数
编号
As(μg/g)
Cd(ng/g)
Cr(μg/g)
Cu(μg/g)
Hg(ng/g)
Ni(μg/g)
Pb(μg/g)
Zn(μg/g)
1
2.01
2.23
2.01
3.74
2.66
1.49
2.23
3.43
2
2.01
3.02
1.72
9.66
18.35
1.61
3.00
4.03
3
1.12
1.17
1.26
1.31
1.17
1.26
1.18
1.06
4
1.59
2.77
1.87
4.71
12.77
1.43
2.05
3.52
5
1.74
2.16
1.41
2.29
3.29
1.24
1.96
2.24

表3 土壤单因子污染指数评价标准
分级指数
质量指数
分级指数
质量指数

清洁

轻污染

潜在污染

重污染

5.1.2内梅罗综合指数法
单项污染指数得到每个区域各种重金属污染程度，为了得到每个区域的综合污染程度，我们利用内梅罗综合指数法求解。内梅罗综合指数式为：
（2）
式中：为综合污染指数, 为各元素污染指数平均值, 为各元素污染指数中的最大值。
通过求解内梅罗综合指数式得到每个区域的综合污染指数（见表7，表8）
表4 各区域综合污染指数

I最大
I平均
综合污染指数P
生活区
1.66055
0.59623
0.882173
工业区
5.244954
1.434308
2,718,767

山区
0.188991
0.070955
0.100936
主干道区
2.248165
0.878292
1.206819
公园绿地
0.885529
0.377438
0.481306

表5 土壤内梅罗综合污染指数评价标准
等级
内梅罗综合污染指数
污染等级
Ⅰ
　≤0.7
清洁
Ⅱ
0.7＜≤1
警戒线
Ⅲ
1＜≤2
轻度污染
Ⅳ
2＜≤3
重度污染
Ⅴ
　3＜
重污染

5.1.3地累积指数法
虽然单因子指数法和内梅罗指数法均能对土壤重金属污染程度进行较为全面的评价，但无法从自然异常中分离人为异常，然而地累积指数法弥补了这项不足。地累积指数法表达式为：
（3）
式中：为地累积指数，是各个区域元素的平均值，是元素地球化学背景值，为各地岩石差异导致的背景值变动系数（一般取值为1.5）。
通过对地累积指数式的求解得到每个区域各种重金属的污染程度（见表9）
表6 每个区域各种重金属的污染程度

评价结果
经过内梅罗综合指数法和地累积指数法对该城不同区域重金属的污染程度分析，得到生活区中、、、、元素属于轻度污染，元素成为无污染元素，、污染为中度污染。综合污染等级为警戒级。工业区中、、元素属于轻度污染，、、、元素成为中度污染，则属于强度污染。综合污染等级为严重污染。山区、、、、、、、八种元素基本都对当地环境无污染，综合污染等级为清洁级。主干道路区、、、、、、 等元素基本为轻度污染，只有为强度污染，综合污染等级为轻度污染。公园绿地区大部分元素为轻度污染，有很少地方存在中度污染，综合污染等级为清洁级。
5.2问题二的建模与求解
根据问题的分析，某些重金属空间分布的含量也具有一定的相关性，相关性 较大的金属可能在成因和来源上有一定的联系。因此，需要对问题中所给出的8种重金属元素之间的关系进行分析。在这里我们选用主成分分析法来求解这个问题。

主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指标。在实证问题研究中，为了全面、系统地分析问题，我们必须考虑众多影响因素。这些涉及的因素一般称为指标，在多元统计分析中也称为变量。因为每个变量都在不同程度上反映了所研究问题的某些信息，并且指标之间彼此有一定的相关性，因而所得的统计数据反映的信息在一定程度上有重叠。在用统计方法研究多变量问题时，变量太多会增加计算量和增加分析问题的复杂性，人们希望在进行定量分析的过程中，涉及的变量较少，得到的信息量较多。
5.2.1对原始数据进行标准化
对原始数据标准化处理的一般公式为：
（4）
其中为第个采样点的第种重金属的浓度值，分别表示第种重金属的均值和标准差。
5.2.2求重金属含量相关系数矩阵
对于任意的两种重金属元素与元素在第个样本中的含量值分别为和（），则二元素相关系数的一般计算公式为
（5）
其中

式中分别为元素与元素样本点观测值的算术平均值。
首先我们通过将生活区的重金属元素原始数据导入软件进行主成分分析处理可以得到各重金属元素之间的相关系数如下表所示。（以生活区为例）

表17 各个因子间的相关系数矩阵

5.2.3求相关系数矩阵的特征值和特征向量
利用相关系数矩阵可求出相应的因子特征值和累计贡献率，用SPSS统计软件计算可得到生活区的各因子特征值和成分矩阵结果如下表所示。
表8 解释的总方差

表9 成分矩阵（载荷矩阵）

主成分个数提取原则为主成分对应的特征值大于1的前m个主成分，通过表21，可知提取3个主成分。用表16中数据除以主成分对应特征值便得到三个成分中每个指标对应的系数，即特征向量A1，A2，A3，将特征向量与标准化后数据相乘得到主成分表达式：F1=A1*ZX1; F2=A2*ZX2; F3=A3*ZX3。为了方便计算我们在这里用SPSS解决这个问题得到成分得分系数矩阵。
表10 各因子得分系数矩阵

以每个主成分所对应的特征值占所提取主成分总的特征值之和的比例作为权重计算主成分综合模型，即用第一主成分F1中每个指标所对应的系数乘上第一主成分F1所对应的贡献率再除以所提取的三个主成分贡献率之和，然后加上第二主成分F2对应系数乘以其贡献率百分比，再加上第三主成分F3对应系数乘以其贡献率百分比，就得到了综合主成分F中每个指标所对应的权重系数A。根据权重系数A值大小降序排列，系数大的对应的重金属元素就是主要污染的主要原因。这里只将生活区列举出来，其它区域分析方法与上述方法相同，就不做过多叙述，其他区域表格见附录。
表11 各区域因子的权重系数
生活区
权重系数
工业区
权重系数
山区
权重系数
Cr
0.34
Ni
0.34
Cu
0.35
Zn
0.34
Zn
0.33
As
0.23
Ni
0.27
As
0.31
Cr
0.23
Pb
0.22
Pb
0.31
Ni
0.22
Cd
0.21
Cd
0.29
Hg
0.21

Cu
0.2
Cr
0.29
Zn
0.17
Hg
0.11
Cu
0.23
Pb
0.06
As
0.1
Hg
0.22
Cd
0.04
主干道区
权重系数
公园绿地
权重系数

Pb
0.38
Pb
0.36

Cd
0.34
Cu
0.27

Cu
0.29
Hg
0.27

Ni
0.27
Zn
0.26

Zn
0.27
Cd
0.23

Cr
0.26
Cr
0.16

Hg
0.23
As
0.14

As
0.06
Ni
0.08

结果评价
对该城区不同区域采集的土壤、、、、、、、元素进行主成分分析，通过主成分分析法得到:
1> 生活区元素、在该城市土壤中的含量明显高于国家背景值。这是由于生活区存在大量的公路，公路两侧含铅汽油的燃烧和汽车轮胎的磨损的粉尘会增加土壤中、等元素的含量，而且由于生活区人类活动比较频繁，日常生活用品丢弃后成为垃圾也会导致这三种元素的增加。
2> 工业区主要污染元素为、、、、、。因为工业区“三废”排放，采矿和冶炼会增加这些元素在土壤中的含量。
3> 山区全部元素几乎都在国家背景值以下，只有元素为主要污染元素。这是由于山区远离城市和工业区的污染源，只有一部分金属矿山的开采会导致一些金属元素的外露。
4> 主干道路区、元素都高于背景值，其余元素都低于背景值。机动车尾气排放既是城市大气的主要污染源，也显著引起公路两侧土壤的重金属污染，汽车汽油、发动机、轮胎、润滑油和镀金部分都能燃烧或磨损而释放出

、元素。
5> 公园绿地区中为主要重金属元素，、、为次要重金属元素。公园绿地区农药化肥的使用会不同程度的影响重金属的污染，而且地下水灌溉和塑料薄膜会增加、、元素的土壤含量。
不难发现每个区域都受到不同程度的元素的污染，这是因为镉元素主要以硫化镉形式储存于锌矿、铅锌矿和铜铅锌矿中，土壤镉主要来源于锌矿、铅矿的冶炼、合金、电镀、化工厂等排放的废水，工业固废堆放，含有废旧电池的生活垃圾渗滤，污泥施肥以及过量或不恰当地使用化肥农药等。
总体来说：工业化程度越高的地区污染越严重，市区高于远郊和农村，地表高于地下，污染区污染时间越长重金属积累就越多，以大气传播媒介土壤重金属污染土壤具有很强的叠加性,熟化程度越高重金属含量就越高。
5.3 问题三模型建立与求解
通过对问题二的分析，我们知道了重金属的传播途径，比如重金属和的污染主要为汽车尾气排放所致，的污染为燃煤引起，它们都是经过空气来传播的，为此用对流-扩散偏微分方程来描述。主要是由于工业废水排放污染，经过地表径流的传播，则可以应用对流-扩散偏微分方程来进行描述。对流-扩散方程是描述粘性流体运动的非线性方程的线性化偏微分方程模型，能够刻画很多自然现象，如污染物的扩散，降解，流体流动与传热和电化学反应等。
首先建立对流-扩散偏微分方程传播模型。然后，依据问题所给的条件，可以知道这是对流-扩散偏微分方程的反问题。在求解时，将对流-扩散偏微分方程进行简化变为二阶椭圆型偏微分方程。利用椭圆型偏微分方程的有限差分数值解法正向估计出污染物的浓度分布，再代入不同的位置和强度参数，与实际检测值相比较，偏差最小的即为污染源的位置及源强，即将问题转化为非线性最优化问题，利用中工具箱求解椭圆型微分方程，则可得到各重金属元素的污染源位置及源强。
5.3.1 对流-扩散偏微分方程模型的建立

一般的对流-扩散偏微分方程模型为
（6）
其中为有界区域，为边界上的外法向量，为轴方向的流速，为轴方向的流速，为扩散系数，为污染物的降解率，为污染物的浓度分布，和分别表示多个污染点和污染强度，为狄拉克函数，即单位脉冲函数。
5.3.2对流-扩散偏微分方程反问题的求解
若已知的分布，并且已知和那么源项识别的反问题就是根据这些已知的分布确定源项，即确定污染源的位置和污染强度。
根据模型的假设，观测值为一段时间内稳定的浓度分布，即二阶稳态对流扩散方程。根据相关文献可知扩散系数=，污染物降解率=0.001，暂时不考虑横纵轴的流速，使轴方向的流速=0，
轴方向的流速=0.于是则可将上述的偏微分方程简化为
（7）
该方程为二阶椭圆型偏微分方程。要在污染区域中大概找出污染范围，首先选择污染源初始位置，以及污染强度。利用椭圆型偏微分方程的有限差分数值解正向估计出污染物的浓度分布，。

的分布与初始位置，以及污染强度的取值有关。用函数表示在初始位置，以及污染强度的条件下测得的污染物浓度分布，，，为第次估计值。将所得到的污染物浓度与实际检测值相比较，差距最小的即为最接近污染源的位置。即将问题转换为非线性最优化问题：
（8）
这里为选取的样点个数。
适时调整参数，，，使上面的问题达到最小值。使其达到最小值的，，即为污染源的位置坐标和污染强度。
5.3.3对流-扩散方程解的仿真模拟
对于简化的对流-扩散方程化为标准的二阶椭圆型偏微分方程：
（9）
用函数来近似表示单位脉冲函数，从而可得到偏微分方程为：
通过以上的分析，结合问题（1）的模型所得到的重金属污染分布结果可知，污染浓度最高的区域最有可能是污染源。因此搜索污染源时以差值所得到的浓度最高点作为搜索的初始值，并在该值附近进行搜索。对每个可能的污染点使用中的工具来求解椭圆型偏微分方程的解。比较不同污染点和源强下的污染分布和实际检测的污染浓度，选择两者误差最小的污染源位置和源强。用此方法可求得八种重金属的污染源位置和污染浓度如下：

表12 各种金属污染源坐标和浓度值
重金属
横坐标/
纵坐标/
浓度值(μg/g)

As
18134
10046
24.5449
13289
2515
19.27
Cd
23664
9790
725.9939
2124
2486
542.6413
Cr
3299
6018
320.3774
Cu
2383
3692
447.5421
Hg
2708
2295
1672.2
14100
2367
1508.1
15185
9121
1326.8
Ni
3299
6018
32.9542
Pb
4777
4897
165.78
Zn
13244
7056
1289.56
9132
4582
1156.7
5.4问题四模型建立与求解
在解决问题三的模型中，将复杂的对流-扩散方程简化得到污染物传播的基本模型，便于直观地从解空间中看到传播趋势，从而便于反问题求解污染源，但是没有考虑横纵轴的流速，对于随位置变化的扩散系数，以及降解率等以常数值来描述，具有局限性。无法满足对城市地质环境研究的需求。需要进一步改进。
5.4.1扩展二维对流扩散方程到三维
首先要研究哪些因素对地质环境的演变影响比较大，比如不仅受地域地形的影响，还受到天气(如风)、雨水冲刷和温度的影响，为了更准确的研究城市地质环境的演变，还需要收集发地区在不同天气条件下各种重金属的浓度分布。如在同一区域内雨天前的浓度分布与雨天后的浓度分布情况，在不同季节、不同温度条件下的浓度分布，以及风的因素随污染物浓度的影响。
因此，要实现这个目的，就需要将问题三中的二维对流-扩散模型扩展到三维的情形，即有三维污染物的扩散模型为
（10）

其中为污染物浓度；为轴方向的流速；为轴方向的流速；为轴方向的流速，表示扩散系数；为自净化系数，，表示沉降、稀释、化学及生物反应等，如当污染物为重金属时表示重金属形态转化系数，当污染物为有机物时表示有机物降解系数，亦可是各因素的综合；表示污染物浓度随时间的变化项。
在均匀稳定流中，三维对流-扩散模型的解为
（11）
其中为源强。
5.4.2 三维对流-扩散偏微分方程参数的确定
在这里主要是估计出更准确的分布参数，从而确定污染物的分布个数，并能够有效地模拟实际中地质模型的演变模式，问题就转换为准确估计三维对流-扩散方程的参数，进而用以模拟城市的地质演变模式。
在坡度和温度基本相同的条件下，考虑风的影响。某些重金属，如和主要来自于燃煤或工业废气引起的径流污染，这就使得风向和风速对这些重金属分布的影响很大，在大气中重金属会随着风的方向而漂移，从而改变了一定区域内的重金属含量浓度。
在不同的风向量条件下，收集区域土地的分布，这里为风向量，即有大小，又有方向。轴方向的流速为
（12）
将上式代入到（15.4）式中得到的数学模型为

当函数表示不同形式时，利用解析解得到不同的重金属浓度分布。将由不同函数所得到的污染物浓度物空间浓度分布与观测值进行比较，则可将问题转化为求解优化问题，即利用多项式逼近函数

，从而确定风对参数影响的近似函数表达式。
首先考虑风速对轴方向的影响，对模型进行仿真模拟得到的污染物浓度的空间分布趋势，说明风对污染物浓度的影响。
同理对于轴和轴方向，以及地形和温度等因素对三维对流-扩散方程的影响，类似的也可以用上面的模型进行估计修正，最后综合不同因素对参数的影响则可以确定参数形式的三维对流-扩散模型。从而可以从该模型出发，更好的研究城市地质演变的模式，有针对性的对城市的地质环境问题进行处理，这样会具有很高的实际应用价值。
六 模型评价
6.1评价
对于问题一,利用重金属浓度在空间的相关性，采样普通克里格插值法，绘出各种重金属的浓度空间分布，然后利用地质积累法和修正的内梅罗综合指数法对各功能区的污染程度进行分析，给出了各功能区被重金属污染的情况，同时给出了每种重金属的污染指数分布，以及在各功能区的污染指数直方图，多方面分析了重金属的污染程度，但克里格插值法有一定的局限性，在很大的程度上受局部采样点含量的影响较大，故对采样点所具有的代表性要求较高。
对于问题二，采用主成分分析法，用较少的有代表性的因子来说明众多因素变量所提取的主要信息，提取出了6个因子后计算其得分，这个模型在一定程度上降低了计算的复杂度，并用合理的方法对污染源进行了分析，其结果与实际情况基本相符，由于问题中所提供的数据是以1公里为网格进行取样的，没有给出具体的地理分布，所以在分析污染源的时是在插值的基础上进行分析的，有一定的模糊性，一些污染物的主要来源是与实际情况相结合确定出来的。
对于问题三，建立了二维对流—扩散偏微分方程来描述重金属污染的传播特征，这是比较合理的，在求解模型时将其简化为椭圆形偏微分方程模型，并将问题转化为非线性最优解，则较容易地得到问题的结果但在简化过程中忽略了一些变量因素对污染传播的影响，具有一定的局限性。

对于问题四，将二维对流—扩散方程扩展到三维，并对各种变化参数进行了讨论，使更充分地描述各参数对城市地质演变过程的影响，但没有给出更详细的在各种变化参数影响下的传播仿真模拟结果。

6.2模型的改进方向
在问题三的模型中，在应用对流-扩散方程的时候，相关参数的确定是根据专家咨询法给出的，有一定的主观性，在实际中为了更准确地模拟重金属的污染扩散，应该根据一定的科学方法适当地修正参数。对于偏微分方程的求解，可以应用更准确的算法进行求解。在实际问题中，对于偏微分方程的求解可以应用蒙特卡洛，等方法，在大量数据的基础上对偏微分方程的参数进行反演，这样可以提高精度，从而更好地应用于研究实际中城市地质环境模式。

七 参考文献

[1]何瑶，刘颖，陈玲，上海市表层土壤中金属元素含量及其分布特征 [J].中 国可持续发展论坛专刊，2008。
[2] 赵静、但琦，数学建模与数学实验[M],北京：高等教育出版，2008。
[3] 李尚志，数学建模竞赛教程[M]，南京：江苏教育出版社，1996。
[4] 王树禾，数学模型基础[M],台肥：中国科学技术大学出版社，1996。
[5] 姜启源，数学模型[M]，北京：高等教育出版社，1987。
[6] 薛定宇,高等应用数学问题的MATLAB求解[M]，北京：清华大学出版社，2004。
[7] 陈杰，MATLAB宝典[M]，北京：电子工业出版社，2007年。

八 附录

程序一：
load('chenshiturangjinshu.mat');
%[x,y]=meshgrid(x_zuobiao,y_zuobiao);

x=x_zuobiao;
y=y_zuobiao;
z=haiba;
[X,Y]=meshgrid(0:50:30000,0:50:20000);
Z=griddata(x,y,z,X,Y);
meshz(X,Y,Z);
shading interp;

程序二：
function pdemodel
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',1);
set(ax,'DataAspectRatio',[15000 10000 1]);
set(ax,'PlotBoxAspectRatio',[1 1 1]);
set(ax,'XLim',[0 30000]);
set(ax,'YLim',[0 20000]);
set(ax,'XTickMode','auto');
set(ax,'YTickMode','auto');
pdetool('gridon','on');

% Geometry description:
pderect([15161.538461538461 19776.923076923071 10207.692307692307 6515.3846153846134],'R1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R1')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(4,...
'dir',...

1,...
'1',...
'0')
pdesetbd(3,...
'dir',...
1,...
'1',...
'0')
pdesetbd(2,...
'dir',...
1,...
'1',...
'0')
pdesetbd(1,...
'dir',...
1,...
'1',...
'0')

% Mesh generation:
setappdata(pde_fig,'Hgrad',1.3);
setappdata(pde_fig,'refinemethod','regular');
setappdata(pde_fig,'jiggle',char('on','mean',''));
pdetool('initmesh')
pdetool('refine')
pdetool('refine')

% PDE coefficients:
pdeseteq(1,...

'1.0',...
'0.001',...
'30*exp(-((x-16000).^2+(y-8000).^2))',...
'1.0',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setappdata(pde_fig,'currparam',...
['1.0 ';...
'0.001 ';...
'30*exp(-((x-16000).^2+(y-8000).^2))';...
'1.0 '])

% Solve parameters:
setappdata(pde_fig,'solveparam',...
char('0','27210','10','pdeadworst',...
'0.5','longest','1','1e-4','','fixed','inf'))

% Plotflags and user data strings:
setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 1 0 1 1 0 1 0 0 1]);
setappdata(pde_fig,'colstring','');
setappdata(pde_fig,'arrowstring','');
setap