2011数学建模A题全国一等奖论文(城市重金属污染及地质演变).pdf

2011 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛 承承 诺诺 书书 我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网 上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。 我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的 资料（包括网上查到的资料） ，必须按照规定的参考文献的表述方式在正文引用处和参 考文献中明确列出。 我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规 则的行为，我们将受到严肃处理。 我们参赛选择的题号是（从 A/B/C/D 中选择一项填写） ： 我们的参赛报名号为（如果赛区设置报名号的话） ： 所属学校（请填写完整的全名） ：中国石油大学（北京） 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日 赛区评阅编号（由赛区组委会评阅前进行编号）： 2010 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛 编编 号号 专专 用用 页页 赛区评阅编号（由赛区组委会评阅前进行编号）： 赛区评阅记录（可供赛区评阅时使用）： 评 阅 人 评 分 备 注 全国统一编号（由赛区组委会送交全国前编号）： 全国评阅编号（由全国组委会评阅前进行编号）： 1 城市重金属污染及地质环境演变分析城市重金属污染及地质环境演变分析 摘要 在问题 1 中，我们首先根据离散的 300 余个样本的坐标及其相应的重金属浓度，用 matlab 分别拟合绘制 8 种重金属元素浓度随坐标位置变化的浓度分布渐变图， 用不同颜 色代表不同的污染程度，定性直观地看出 8 种重金属元素在该城区的空间分布。同时， 我们引入污染指数I，利用地累积指数法对 8 种元素在 5 块区域的污染程度进行定量计 算，并且对污染程度进行分级，不仅定量分析了不同区域的重金属污染情况，而且对第 一部分的拟合图做了检验，两部分吻合很好。 在问题 2 中，我们建立主成分分析模型，运用主成分分析法，用较少的新变量（主 成分）代替原来较多的旧变量（8 种重金属元素） ，既减少了变量个数，又保证了新变量 能够覆盖大部分信息。 这样便可以比较容易地通过主成分挖掘出各污染元素与区域功能 之间的关系，进而确定重金属污染物的来源。 在问题 3 中，我们分析了重金属污染物在地层中的扩散属于“分子扩散”而无“对 流扩散” ，符合菲克扩散定律的传播特征。我们根据空间中任一体积元内金属污染物的 物料守恒，结合反映金属元素传播特征的菲克扩散定律，建立金属元素在地层中的扩散 模型，得到浓度关于坐标的拉普拉斯方程。我们采用有限差分法进行数值求解，在问题 一浓度分布图的基础上，利用直观浓度最大点附近的若干已知样本数据作为边界条件， 借助 Excel 进行迭代运算，使所得到的数值解与真实情况最大误差为 0.001。再用 matlab 筛选出浓度最大点的坐标即为要求的污染源位置。由此分别求得 8 种重金属的污染源。 在问题 4 中，我们首先分析了研究地质环境转换的具体内容，也就是既要研究不同 区域地质环境状况的“横向对比” ，又要研究不同时间地质环境状况的“纵向对比” ，对 前三问建立的模型进行了优缺点分析。并且提出要从物质组成、地质结构、动力作用三 个方面补充收集 14 类信息，建立层次分析模型，得到了一套地质环境状况评价指标， 也就是各因素所占权重（权重的变化可以看出演化的过程）及综合评价指标。再将不同 时期的评价体系进行对比，即可得出城市地质环境在的演变模式。 关键词：地污染指数；主成分分析；菲克定律；层次分析 2 一、 问题重述 为了对某城市城区土壤地质环境进行调查， 研究人类活动影响下城市地质环境的演 变模式，现将城区按功能划分为生活区、工业区、山区、主干道路区及公园绿地区等， 分别记为 1 类区、2 类区、5 类区，进行研究，不同的区域环境受人类活动影响 的程度不同。 将所考察的城区划分为间距 1 公里左右的网格子区域， 按照每平方公里 1 个采样点 对表层土进行取样、编号，并用 GPS 记录采样点的位置。应用专门仪器测试分析，获得 了每个样本所含的多种化学元素的浓度数据。另一方面，按照 2 公里的间距在那些远离 人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。 这样便得到了采样点的位置、海拔高度及其所属功能区；8 种主要重金属元素在采 样点处的浓度；8 种主要重金属元素的背景值。 依据以上信息，通过数学建模解决以下问题： （1）给出 8 种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重 金属的污染程度。 （2）通过数据分析，说明重金属污染的主要原因。 （3）分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 （4）分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收 集什么信息？有了这些信息，如何建立模型解决问题？ 二、 问题分析 对于问题 1，分析题目要求可以看出它有两部分构成：重金属元素在各区域的空间 分布分析， 以及不同区域重金属污染程度评价。 要分析重金属元素在各区域的空间分布， 结合题目所给出的离散的样品点， 我们想到用matlab拟合出每种重金属元素的分布图， 可以直观地看出重金属的分布特征。对于第二部分，要想对不同区域重金属的污染程度 进行评价，只要引入一个污染评价指标便可以定量评价每个区域的污染程度。 对于问题 2，要分析重金属污染的原因，也就是各种金属的来源，需要找出区域的 功能与不同污染元素之间的关系。但由于污染元素一共有 8 种，变量因子过多，且原始 数据纷杂，不利于找出内在联系，于是想到用较少的新变量代替原来较多的旧变量，而 且保证新变量尽可能多的保留原来变量所反映的信息。 这样就可以通过分析少量新变量 之间的关系，挖掘出内在联系。主成分分析法便可以实现这一目的。 对于问题 3，要得到污染源，必须要得到地层中重金属元素的浓度分布，以求浓度 最大点。首先我们要明确重金属元素在地层中的扩散规律，然后根据这一规律建立扩散 模型。通过查阅大量文献，我们总结出重金属元素在地层中的扩散只有“分子扩散” ， 而无“对流扩散” ，满足菲克定律。这样我们根据特定地层区域内重金属“物料守恒” 建立了扩散模型，求得了地层中的浓度分布，进而得到浓度最大点也就是污染源。 对于问题 4，要研究地质环境演变时，既要研究不同区域地质环境状况的“横向对 比” ，又要研究不同时间地质环境状况的“纵向对比” ，并且要明确地质环境这一概念所 包含的因素（重金属污染只是其中一种） 。明确了以上两点也就明确了我们研究地质环 境演变规律所需要的信息。而在“横向对比”过程中，我们必须要对地质环境状况有一 3 套评价指标，也就是各因素所占权重（权重的变化可以看出演化的过程）及综合评价指 标。 “纵向对比”只需要对比不同时间同一区域的地质环境状况即可。 三、 模型假设 重金属污染物浓度在空间的分布连续； 重金属污染物在地层中的扩散系数为常数； 该城区的重金属污染物分布在短时间内不随时间变化，是稳定状态。 四、 模型建立与求解 4.1 重金属元素分布与不同区域污染程度（问题 1） 4.1.1 各重金属元素在不同区域的分布图 要分析8种重金属元素在该城区的空间分布， 我们想到了根据给出的样本测量数据， 把 8 种元素分开讨论，剔除样本中的特异点后，对于每一种元素都用 matlab 软件的 “griddata” 语句绘制出该种元素在 5 个区域内的连续分布 （在不同区域内的浓度大小） 的网格图，进而直观地看出不同区域每一种重金属污染的程度。如下图，是通过 matlab 绘制的各种重金属元素的污染分布图（绘图源代码见 附录Matlab 程序第 1 问） 。 由于图片所占页面空间太大，这里只贴出 As、Cd 元素的分布图作为代表，剩下 6 种元素的分布图见附件。 对于 8 种重金属元素的分布图来说，红色表示浓度最大，紫色表示浓度最小。颜色 由暖色至冷色表示浓度递减。 图 1：As 元素在各区域的分布 4 从图 1 中可以直观地看出 As 元素主要分布在工业区、生活区与绿地周围，而在交 通区分布较少。在山区大部分地区分布较少，少部分地方有分布。 图 2：Cd 元素在各区域的分布 对于图 2， 可以直观地看出 Cd 元素主要分布在工业区与交通区周围， 次之分布在生 活区与绿地，在山区分布较少。 另外，由其它 6 种重金属元素的分布图（见附件图像第 1 问）可以看出 Cr 元 素主要分布在生活区与交通区，在工业区、绿地与山区分布较少。 对于 Cu 元素，发现工业区与交通区分布较多，其它地方分布较少。 对于 Hg 元素，发现只有少部分生活区与交通区分布较多，其它部分基本上没有受 到 Hg 元素的污染。 对于 Ni 元素， 发现在山区、 生活区与工业区分布较多， 在交通区与绿地分布较少。 对于 Pb 元素，发现在交通区分布最多，在部分生活区与部分工业区也有所分布， 而在其它区域都分布较少。 对于 Zn 元素，发现在生活区与工业区分布较多，其它区域分布较少。 4.1.2 地累积指数法评价不同区域重金属污染程度 地累积指数法是一种常用的研究沉积物中重金属污染的定量指标，其计算公式为： 2 log n geo n C I KB （1） 式中： n C是实测元素n在地层中的含量； n B为元素n在地层中的背景值；K是考虑 了各地岩石差异可能会引起的变动而取得变动系数（一般情况下取值为 1.5） 。 【1】 得到的地积累指数可分为不同的等级，分级情况与污染程度的对应关系见下表： 5 表 1：重金属污染级别、地积累指数（ geo I）和分级比较 项目 污染程度分类 清洁 轻度污染 偏中度污染 中度污染 偏重污染 重污染 严重污染 geo I 5 级别 0 1 2 3 4 5 6 对于本题， n B是题目中所给的各种元素的背景值。 n C是通过题目所给的数据统计 得到的各区域中各种重金属元素的平均含量，如下表： 表 2：各区域中各种重金属元素的平均含量 元素平 均含量 As Cd Cr Cu Hg Ni Pb Zn 生活区 6.270 289.9 69.01 49.40 93.04 18.34 69.10 237.0 工业区 7.251 393.1 53.40 127.5 642.3 19.81 93.04 277.9 山区 4.044 152.3 38.9 17.31 40.95 15.45 36.55 73.2 交通区 5.703 359.0 58.02 62.02 446.1 17.59 63.39 242.3 绿地区 6.104 263.1 42.83 30.04 109.4 14.88 60.19 145.9 把, nn B C代入到式 （1） 后， 得到 5 块区域内各种重金属元素的地积累指数值 （ geo I） ， 如下表： 表 3：各区域内各种重金属元素的地积累指数值 As Cd Cr Cu Hg Ni Pb Zn 生活区 geo I 0.8 1.157 1.154 1.904 1.41 0.58 1.156 1.78 污染级别 1 2 2 2 2 1 2 2 工业区 geo I 1.01 1.596 0.784 3.272 4.197 0.69 1.585 2.01 污染级别 2 2 1 4 5 1 2 3 山区 geo I 0.167 0.228 0.32 0.391 0.226 0.33 0.237 0.087 污染级别 1 1 1 1 1 1 1 1 交通区 geo I 0.663 1.465 0.9 2.23 3.672 0.52 1.032 1.812 污染级别 1 2 1 3 4 1 2 2 绿地区 geo I 0.761 1.01 0.46 1.186 1.645 0.28 0.957 1.081 污染级别 1 2 1 2 2 1 1 2 通过上表各区域的污染级别可以看出： 对于生活区，Cd、Cr、Cu、Hg、Pb、Zn 这几种重金属元素属于偏中度污染；As、Ni 两种重金属元素属于轻度污染。 对于工业区，Hg 是重污染；Cu 是偏重污染；Zn 是中度污染；As、Cd、Pb 属于偏中 度污染；Cr、Ni 属于轻度污染。 对于山区，所有的重金属元素都是轻度污染。 对于交通区，Hg 是偏重污染；Cu 是中度污染；Cd、Pb、Zn 均属于偏中度污染；As、 Ni 属于轻度污染。 对于绿地区，Cd、Cu、Hg、Zn 属于偏中度污染；As、Cr、Ni、Pb 均属于轻度污染。 另外，为了对比各区域的总污染程度，我们对各个区域内的所有重金属元素的地积 6 累指数求和，得到各区域的总地积累指数（ tol I） ，结果如下表： 表 4：各区域的总地积累指数（ tol I） 区域 生活区 工业区 山区 交通区 绿地区 tol I 1.2 1.9 0.24 1.5 0.9 通过上表比较可以发现不同区域的污染程度为： 工业区交通区生活区绿地区山 区。 通过对比 4.1.1 与 4.1.2 的结论， 发现重金属元素的分布情况与各区域内的重金属 污染程度相互呼应，吻合的很好。该城区的重金属污染特点，我们可以总结为：该城区 的主要污染元素为 Hg、Cu、Zn；且污染主要集中在工业区与交通区。 4.2 主成分分析模型（问题 2） 要研究重金属污染的主要原因， 就要发现不同区域内主要污染物与该区域功能的关 联性，进而确定各区域污染原因之间的关联性，以确定整个城区污染的原因。 但由于导致每个区域污染的重金属元素众多， 没法直接通过重金属元素之间的关系 确定各区域污染原因之间的关联性。于是想到用较少的新变量代替原来较多的旧变量， 而且保证新变量尽可能多的保留原来变量所反映的信息。主成分分析法便可以实现，主 成分分析法是把原来多个变量划为少数几个综合指标的一种统计分析方法。 4.2.1 主成分分析模型的建立 假定有n个地理样本，每个样本共有p个变量描述，这样就构成了一个np阶的地 理数据矩阵： 11121 21222 12 p p nnnp xxx xxx X xxx （2） 要想从这么多变量的数据中抓住地理事物的内在规律，自然要在p维空间中加以考 察，这是比较麻烦的。为了克服这一困难，就需要用较少的几个综合指标来代替原来较 多的变量指标，这些综合指标（即新变量）选取最简单的形式就是取原来变量指标的线 性组合，适当调整组合系数，使新的变量指标之间相互独立且代表性最好。 若原来的变量指标为 12 ,., p x xx， 它们的综合指标新变量指标为 12 ,.,() m z zzmp， 则 111 11221 221 12222 1 122 , , . , pp pp mmmmpp zl xl xl x zl xl xl x zl xlxlx （3） 在（3）式中，系数 ij l由下列原则来决定： , ij z z（; ,1,2,.,ij i jm）相互无关； 7 1 z是 12 ,. p x xx的一切线性组合中方差最大者； 2 z是与 1 z不相关的 12 ,., p x xx的所 有线性组合中方差最大者； m z是与 121 ,., m z zz 都不相关的 12 ,., p x xx的所有线性 组合中方差最大者。 这样决定的新变量指标 12 ,., m z zz分别称为原变量指标 12 ,., p x xx的第一， 第二， ， 第m主成分。其中， 1 z在总方差中占的比例最大， 23 ,., m z zz的方差依次递减。在实际 问题的分析中，常挑选前几个最大的主成分，这样既减少了变量的数目，又抓住了主要 矛盾，简化了变量之间的关系。 从以上分析可以看出，找主成分就是确定原来变量 j x(1,2,.,jp)在诸主成分 i z （1,2,.,im）上的载荷 ij l（1,2,., ;1,2,.,im jp） ，从数学上容易知道，它们分别 是 12 ,., p x xx的相关矩阵m个较大的特征值所对应的特征向量。 模型计算步骤： 计算相关系数矩阵： 11121 21222 12 . . . p p pppp rrr rrr R rrr （4） ( ,1,2,., ) ij r i jp为原变量 i x与 j x的相关系数， ijji rr，计算公式为： 1 22 11 ()() ()() n kiikjj k ij nn kiikjj kk xxxx r xxxx （5） 计算特征值与特征向量。 先解特征方程0RE，常用雅克比法求出特征值，并使其按大小顺序排列： 12 .0 p 。 然后， 分别求出对于特征值 i 的特征向量(1,2,., ) i e ip， 要求1 i e ， 即 2 1 1 p ij j e ， 其中 ij e表示向量 i e的第j个分量。 计算主成分贡献率及累计贡献率。 贡献率： 1 i p k k （6） 累计贡献率： 1 1 i k k p k k （7） 一般累计贡献率达到 80%以上的特征值 12 ,., m 所对应的为第 1、第 2、第 ()m mp个主成分。 计算主成分载荷 ( ,)( ,1,2,., ) ijijiij lp z xe i jp （8） 8 4.2.2 模型求解 计算相关系数矩阵 对所给的数据进行分区域标准化处理，利用式（5）得到 5 个区域的相关系数矩阵 如下： （具体计算源代码见 附录matlab相关系数） 生活区： ij r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 r 1 r 1.00 0.26 0.36 -0.43 -0.41 -0.28 -0.37 -0.50 2 r 0.26 1.00 0.29 0.08 -0.20 -0.07 0.05 -0.36 3 r 0.36 0.29 1.00 -0.14 -0.06 -0.16 -0.09 -0.10 4 r -0.43 0.08 -0.14 1.00 0.38 0.24 0.53 0.29 5 r -0.41 -0.20 -0.06 0.38 1.00 0.35 0.50 0.40 6 r -0.28 -0.07 -0.16 0.24 0.35 1.00 0.38 0.15 7 r -0.37 0.05 -0.09 0.53 0.50 0.38 1.00 0.20 8 r -0.50 -0.36 -0.10 0.29 0.40 0.15 0.20 1.00 工业区： ij r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 r 1 r 1.00 0.35 0.51 -0.56 -0.48 -0.44 -0.32 -0.34 2 r 0.35 1.00 0.07 -0.21 -0.37 -0.24 -0.29 -0.30 3 r 0.51 0.07 1.00 -0.26 -0.30 -0.19 -0.13 -0.15 4 r -0.56 -0.21 -0.26 1.00 0.33 0.38 0.15 0.18 5 r -0.48 -0.37 -0.30 0.33 1.00 0.54 0.57 0.53 6 r -0.44 -0.24 -0.19 0.38 0.54 1.00 0.92 0.90 7 r -0.32 -0.29 -0.13 0.15 0.57 0.92 1.00 0.98 8 r -0.34 -0.30 -0.15 0.18 0.53 0.90 0.98 1.00 山区： ij r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 r 1 r 1.00 0.03 0.03 -0.17 -0.07 -0.01 -0.11 0.12 2 r 0.03 1.00 0.61 -0.18 0.28 -0.23 -0.30 -0.09 3 r 0.03 0.61 1.00 -0.33 0.23 -0.18 -0.37 -0.05 4 r -0.17 -0.18 -0.33 1.00 -0.29 0.11 0.53 0.08 5 r -0.07 0.28 0.23 -0.29 1.00 0.07 0.09 0.25 6 r -0.01 -0.23 -0.18 0.11 0.07 1.00 0.36 -0.01 7 r -0.11 -0.30 -0.37 0.53 0.09 0.36 1.00 0.51 8 r 0.12 -0.09 -0.05 0.08 0.25 -0.01 0.51 1.00 交通区： ij r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 r 1 r 1.00 0.51 0.38 -0.11 -0.24 -0.27 -0.28 -0.08 9 2 r 0.51 1.00 0.32 -0.09 -0.12 -0.10 -0.09 -0.12 3 r 0.38 0.32 1.00 -0.11 -0.13 -0.17 -0.16 -0.07 4 r -0.11 -0.09 -0.11 1.00 0.12 0.14 0.09 0.00 5 r -0.24 -0.12 -0.13 0.12 1.00 0.37 0.42 0.21 6 r -0.27 -0.10 -0.17 0.14 0.37 1.00 0.89 0.01 7 r -0.28 -0.09 -0.16 0.09 0.42 0.89 1.00 0.03 8 r -0.08 -0.12 -0.07 0.00 0.21 0.01 0.03 1.00 绿地区： ij r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 r 1 r 1.00 0.51 0.04 -0.40 -0.10 -0.26 -0.11 -0.33 2 r 0.51 1.00 0.03 0.12 0.05 0.32 0.15 -0.21 3 r 0.04 0.03 1.00 0.17 0.01 0.01 -0.07 -0.01 4 r -0.40 0.12 0.17 1.00 0.20 0.64 0.10 0.12 5 r -0.10 0.05 0.01 0.20 1.00 0.50 0.54 -0.01 6 r -0.26 0.32 0.01 0.64 0.50 1.00 0.36 -0.03 7 r -0.11 0.15 -0.07 0.10 0.54 0.36 1.00 0.13 8 r -0.33 -0.21 -0.01 0.12 -0.01 -0.03 0.13 1.00 计算特征值、特征向量，求贡献率与累计贡献率 根据式（6） （7）可以求得各个区域的主成分的贡献率与累计贡献率，见下表： （具 体计算源代码见 附录matlab 程序第 2 问） 表 5：生活区特征值及贡献率 表 6：工业区特征值及贡献率 主成分 特征值 贡献率 累计贡献率 主成分 特征值 贡献率 累计贡献率 1 2.97 0.37 0.37 1 3.96 0.49 0.37 2 1.42 0.18 0.55 2 1.51 0.19 0.68 3 0.97 0.12 0.67 3 0.91 0.11 0.80 4 0.83 0.10 0.77 4 0.73 0.09 0.89 5 0.58 0.07 0.85 5 0.49 0.06 0.95 6 0.45 0.06 0.90 6 0.33 0.04 0.99 7 0.41 0.05 0.95 7 0.07 0.01 1.00 8 0.38 0.05 1.00 8 0.01 0.00 1.00 表 7：山区特征值及贡献率 表 8：交通区特征值及贡献率 主成分 特征值 贡献率 累计贡献率 主成分 特征值 贡献率 累计贡献率 1 3.96 0.31 0.37 1 2.62 0.33 0.37 2 1.51 0.19 0.50 2 1.47 0.18 0.51 3 0.91 0.14 0.64 3 1.06 0.13 0.64 4 0.73 0.13 0.76 4 0.96 0.12 0.76 5 0.49 0.10 0.87 5 0.70 0.09 0.85 6 0.33 0.06 0.93 6 0.63 0.08 0.93 7 0.07 0.04 0.97 7 0.46 0.06 0.99 8 0.01 0.03 1.00 8 0.10 0.01 1.00 10 表 9：绿地区特征值及贡献率 主成分 特征值 贡献率 累计贡献率 1 2.38 0.30 0.37 2 1.73 0.22 0.51 3 1.21 0.15 0.67 4 0.93 0.12 0.78 5 0.84 0.10 0.89 6 0.44 0.06 0.94 7 0.27 0.03 0.98 8 0.20 0.02 1.00 一般情况下认为前几个较大的主成分累计贡献率达到 65%便可以代替原来的多个指 标，能基本覆盖其中的信息。 所以，根据以上 5 个区域的贡献率表格可以看出，对于生活区只需求出前 3 个主成 分 123 ,z z z即可；对于工业区只需求出前 2 个主成分 12 ,z z即可；对于山区，只需求出前 3 个主成分 123 ,z z z即可； 对于交通区， 只需求出前 3 个主成分 123 ,z z z即可； 对于绿地区， 只需求出前 3 个主成分 123 ,z z z即可。 求主成分载荷 对于 5 个区域，利用式（8） ，分别求出原变量 128 ,.,x xx在各个选定的较大主成分 上的载荷。如下表： （具体算法见 附录matlab第 2 问；详细数据见 附录附表 第 2 问） 表 10：各区域主成分载荷 原变量 主成分 生活区 工业区 山区 Z1 Z2 Z3 Z1 Z2 Z1 Z2 Z3 X1（As） -0.7 -0.9 -0.8 0.0 0.1 -0.1 -0.2 0.4 X2（Cd） -0.3 -0.2 0.6 0.0 -0.1 0.0 0.8 0.3 X3（Cr） 0.3 0.3 -0.1 0.0 0.0 -0.1 -0.6 -0.6 X4（Cu） 0.3 -0.3 -0.5 0.1 -0.2 0.4 -0.3 0.3 X5（Hg） 0.1 -0.4 0.4 0.0 0.0 0.1 -0.3 0.5 X6（Ni） 0.0 0.0 -0.1 -0.1 0.4 0.2 0.1 -0.2 X7（Pb） -0.3 0.3 -0.1 0.2 -0.1 -0.4 0.0 0.0 X8（Zn） -0.3 -0.1 -0.1 -0.1 -0.1 0.2 0.1 -0.2 原变量 主成分 交通区 绿地区 Z1 Z2 Z3 Z1 Z2 Z3 X1（As） 0.0 1.2 -0.1 -0.6 0.8 0.5 X2（Cd） 0.0 -0.8 0.1 0.6 -0.5 -0.1 X3（Cr） 0.0 -0.1 -0.1 -0.1 -0.2 -0.1 X4（Cu） 0.0 0.0 -0.1 0.2 0.6 0.0 X5（Hg） 0.0 0.1 0.7 0.3 0.0 0.6 X6（Ni） -0.6 0.1 -0.2 -0.4 -0.2 0.1 X7（Pb） 0.5 0.0 -0.1 0.0 0.1 -0.3 X8（Zn） 0.0 0.0 -0.1 0.0 0.0 0.1 11 （说明：由于计算载荷时我们只精确到小数点后一位，所以，有的载荷值经过“四 舍五入” ，为 0.0） 结果分析： 对于生活区，第一主成分的特点表现在因子变量在 Cr、Cu、Hg 的浓度上有较高的 正载荷，反映了 Cr、Cu、Hg 的富集程度；第二主成分的特点表现在因子变量 As 的浓度 上有很大的正载荷，反映了 As 的富集程度；第三主成分的特点表现在因子变量 Hg 上有 很大的正载荷，反映了 Hg 的富集程度。 对于工业区，第一主成分的特点表现在因子变量 Cu、Pb 上有较高的正载荷，反映 了 Cu、Pb 的富集程度；第二主成分的特点表现在因子变量 As、Ni 上有较高的正载荷， 反应了 As、Ni 的富集程度。 对于山区，第一主成分的特点表现在因子变量 Cu、Ni、Zn 上有较高的正载荷，反 应了 Cu、Ni、Zn 的富集程度；第二主成分的特点表现在因子变量 Cd 上具有很高的正载 荷，反应了 Cd 的富集程度；第三主成分的特点表现在因子变量 As、Cd、Cu、Hg 上具有 很高的正载荷，反映了 As、Cd、Cu、Hg。 对于交通区，第一主成分的特点表现在因子变量 Pb 上具有较高的正载荷，反映了 Pb 的富集程度；第二主成分的特点表现在因子变量 As 上具有很高的正载荷，反映了 As 的富集程度；第三主成分的特点表现在因子变量 Cd、Hg 上具有较高的正载荷，反映了 Cd、Hg 的富集程度。 对于绿地区，第一主成分的特点表现在因子变量 Cd、Hg 上具有较高的正载荷，反 映了 Cd、Hg 的富集程度；第二主成分的特点表现在因子变量 As、Cu 上具有较高的正载 荷，反映了 As、Cd 的富集程度；第三主成分的特点表现在因子变量 As、Hg 上具有较高 的正载荷，反映了 As、Hg 的富集程度。 总结上面的分析可以发现：对于 As，在 5 个区域的第一主因素上都没有反应，在生 活区、工业区、交通区以及绿地区的第二主因素上有所反应，说明 As 元素的污染主要 来源于人类生活、工业排放、汽车尾气排放以及大气沉降等；对于 Cd，只在绿地区的第 一主因素上，可能来源与汽车尾气或大气沉降；对于 Cr，只在生活区的第一主元素上有 所反应，表明起污染主要来源于人类生活；对于 Cu，在生活区、工业区以及山区的第一 主元素上均有所反应，说明其污染来源于人类生活、工业排放和土壤本身；对于 Hg，在 生活区与绿地区的第一主因素上均有所反应，说明其主要来源于人类生活和大气沉降； 对于 Ni，只在山区的第一主因素上有所反应，说明其来源于土壤本身；对于 Pb，发现 其在交通区与工业区的第一主因素上有所反应， 说明其来源于汽车尾气排放与工业排放； 对于 Zn，发现其只在山区的第一主因素上有所反应，说明其污染主要来源与土壤本身。 4.3 重金属在地层中的扩散模型（问题 3） 4.3.1 重金属元素在地层中的扩散特点 由于我们研究表层土壤中重金属元素的分布， 所以其扩散规律不涉及到随流体流动 而产生的“对流扩散” ，也就是说其扩散情况只有“分子扩散” 。分子扩散完全是因为组 分分布不均匀，即在空间中存在浓度梯度，导致这些组分依靠分子热运动从高浓度带扩 散到低浓度带，最后趋于平衡。根据描述物质“分子扩散”现象的规律菲克扩散定 律，在单位时间内通过垂直于扩散方向的单位截面积的扩散物质流量（称为扩散通量） 与浓度沿该截面外法向的方向导数成正比，也就是说，浓度沿截面法向单位距离的变化 12 率越大，扩散通量越大。 4.3.2 扩散模型的建立 在地层中，若各点金属元素浓度不同，就会有重金属元素由高浓度带向低浓度带扩 散。确定污染源的位置，就是要找到重金属元素的浓度最高点，故问题归结为研究地层 中重金属元素的浓度分布， 即根据重金属元素在地层中的扩散特点建立浓度关于位置的 函数关系。 假设地层中任意一点( , , )x y z在时刻t的重金属污染物浓度为( , , )u x y z。假设在时间 区间, t tdt内，流过点( , , )M x y z处面积为ds的曲面重金属流量为dQ，则由菲克定律 的得： ( , , )( , , ) u dQk x y zdSdtk x y z gradu n dSdt n （9） 其中比例系数( , , )k x y z称为重金属的扩散系数（大于0），在一般情况下都把它作 为常数处理。n是曲面的外法线单位向量。负号表示流量方向与浓度梯度方向相反。 在地层中任取一封闭曲面S，记它所包围的区域为V。于是，从时刻t1到时刻t2， 通过曲面S流出该区域V的重金属污染物质量为： 222 111 1 ()()() ttt ttt SVV Qkgradu nds dtdiv k gradu dv dtkudv dt （10） 若V内有污染源不断产生污染物， 设污染源密度为( , , , )F x y z t（在时刻t于任一点处 单位时间、单位体积所产生的污染物质量），则在时刻t1到t2内产生的污染物质量为： 2 1 2 ( , , , ) t t V QF x y z t dv dt （11） 由于物质扩散，在时刻t1到t2内，V中的污染物质量的改变量为： 22 11 321 ( , , , )( , , , ) tt tt VVV uu Qu x y z tu x y z tdvdt dvdv dt tt （12） 由物料守恒定律有： 321 QQQ，即 ( , , , ) u k uF x y z t t （13） 式（13）就是我们建立的重金属污染物在地层中的扩散模型。根据这个方程便可以 求得( , , , )uu x y z t。 由于我们研究的是该城区地层在这一确定时刻的污染情况，所以u与时间无关，即 u t 的值为0。另外，在研究该问题的时候，我们在短时间内不考虑新的污染物的注入， 所以( , , , )0F x y z t 。而k又为常数，这样我们的模型可以简化为0u 即： 222 222 0 uuu xyz （14） 由于方程 （14） 是偏微分方程,只有在问题具有高度对称的情况下,才能求出解析解, 但我们研究的重金属元素在地层中的扩散受到地貌、浓度分布、以及其他复杂因素的影 响，显然不能符合高度对称这一条件。对于这些情况,只能寻求它的数值解。 【2】 13 4.3.3 模型求解（拉普拉斯方程的数值求解） 有限差分法是求解拉普拉斯数值解的一种常用算法 【3】，其思想是用差分 (, )( ,) , u xx yu x yy xy 代替导数, uu xy ，用网格将求解区域覆盖，对于平面拉普拉 斯方程，第 i 行，第 j 列小格的u值由拉普拉斯五点差分格式给出： 1,1,1,1 1/4() ijijiji ji j uuuuu （15） 我们根据这个算法来求上述扩散模型的拉普拉斯方程的数值解，以确定 As 元素的 一个污染源为例，具体操作步骤如下： 在由问题 1 得到的 As 元素浓度分布图上，大致确定浓度最高点的位置，取以 （12380，10043） ， （12380,3532） ， （15107，10043） ， （15107，3532）为顶点的矩形框 将其及附近的 8 个已知浓度点圈在内； 将矩形框划分为 125136 个边长为 20m 的小正方形，得到 126137 个交点，将 8 个已知浓度点的位置分别取为离其最近的交点处； 用 Excel 中的单元格代表交点， 单元格中的数值代表 As 元素浓度， 建立 126137 的矩阵， 在相应单元格中输入 8 个已知的浓度值， 在其他未知浓度值的单元格中输入 0， 将 Excel 表初始化； 如图 3 所示，利用 Excel 作为计算工具利用公式（15）对各浓度值进行 32767 次 迭代运算，使其最大误差为 0.001。 根据以上步骤，我们求得了矩形区域内 126137 个点处的 As 元素浓度值，即求得 了该范围内满足 8 个已知点浓度条件的拉普拉斯方程的数值解。利用 matlab 编程找到 其中的浓度最大点，即为该区域 As 元素污染源。根据以上步骤，我们在 Excel 中输入 相关数据以及公式，进行迭代计算，得到了 8 种元素的浓度最高的坐标，也就是污染源 的位置。下图是 As 元素一个污染源的计算截图。 图 3：Excel 执行有限差分法截图 这样，我们便得到了 8 种重金属元素的污染源位置坐标，如下表： 14 表 11：8 种重元素的污染源位置坐标 As X/m Y/m Z/m g/g Hg X/m Y/m Z/m ng/g 1 18234.36 10063.09 16.05 30.02 1 6946.42 7267.79 25.06 1860.30 2 27496.26 12299.33 72.86 18.00 2 8683.03 12112.98 4.76 1756.89 3 12735.11 2981.66 22.45 23.69 3 22286.44 10435.80 39.60 1730.44 4 4052.08 7826.85 10.06 22.50 Cd X/m Y/m Z/m ng/g Cu X/m Y/m Z/m g/g 1 21418.14 11367.57 39.89 1600.14 1 1736.61 2795.30 12.41 303.14 2 17655.49 3727.07 31.04 1232.87 2 10419.64 8013.20 52.52 286.39 Pb X/m Y/m Z/m g/g Ni X/m Y/m Z/m g/g 1 4777.00 4897.00 7.16 472.48 1 22286.44 12299.33 83.32 62.56 （说明：每种元素的含量均有多个峰值，也就是有多个污染源。序号 1,2表示污 染源的编号） 4.4 城市地质环境演变模式分析（问题 4） 4.4.1 所建模型局限性 在前 3 问建立的所有模型中， 城区地层评价指标都针对于城区以及不同区域的重 金属污染物而建立的，这些评价指标都只能评价重金属污染状况，而重金属污染状况只 是地质环境的很小的一个方面， 所以我们前面所建立的模型虽然能比较准确的评价地层 重金属污染的情况，但却不足以全面地反应地层的地质环境状况，模型反应的信息太过 于局限； 要研究该城区地质环境的演化模式，必须要对该城区的地质状况进行纵向对比， 也就是不同时间的地质环境状况对比。而我们在前 3 问建立的模型，都是反映该城区地 层在这一特定时间的地质环境状况，不包含时间t这一参数，也就不能对地质环境状况 进行纵向对比。反映不出地质环境状况随时间的变化情况。 4.4.2 研究地质环境状况演变所需要的信息 地