数学建模论文-城市重金属污染及地质环境演变分析.doc

2011高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从a/b/c/d中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）：中国石油大学（北京） 参赛队员 (打印并签名) ：1.崔天赐 2. 李彦奇 3. 刘逸飞 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2010高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：城市重金属污染及地质环境演变分析摘要在问题1中，我们首先根据离散的300余个样本的坐标及其相应的重金属浓度，用matlab分别拟合绘制8种重金属元素浓度随坐标位置变化的浓度分布渐变图，用不同颜色代表不同的污染程度，定性直观地看出8种重金属元素在该城区的空间分布。同时，我们引入污染指数，利用地累积指数法对8种元素在5块区域的污染程度进行定量计算，并且对污染程度进行分级，不仅定量分析了不同区域的重金属污染情况，而且对第一部分的拟合图做了检验，两部分吻合很好。在问题2中，我们建立主成分分析模型，运用主成分分析法，用较少的新变量（主成分）代替原来较多的旧变量（8种重金属元素），既减少了变量个数，又保证了新变量能够覆盖大部分信息。这样便可以比较容易地通过主成分挖掘出各污染元素与区域功能之间的关系，进而确定重金属污染物的来源。 在问题3中，我们分析了重金属污染物在地层中的扩散属于“分子扩散”而无“对流扩散”，符合菲克扩散定律的传播特征。我们根据空间中任一体积元内金属污染物的物料守恒，结合反映金属元素传播特征的菲克扩散定律，建立金属元素在地层中的扩散模型，得到浓度关于坐标的拉普拉斯方程。我们采用有限差分法进行数值求解，在问题一浓度分布图的基础上，利用直观浓度最大点附近的若干已知样本数据作为边界条件，借助excel进行迭代运算，使所得到的数值解与真实情况最大误差为0.001。再用matlab筛选出浓度最大点的坐标即为要求的污染源位置。由此分别求得8种重金属的污染源。在问题4中，我们首先分析了研究地质环境转换的具体内容，也就是既要研究不同区域地质环境状况的“横向对比”，又要研究不同时间地质环境状况的“纵向对比”，对前三问建立的模型进行了优缺点分析。并且提出要从物质组成、地质结构、动力作用三个方面补充收集14类信息，建立层次分析模型，得到了一套地质环境状况评价指标，也就是各因素所占权重（权重的变化可以看出演化的过程）及综合评价指标。再将不同时期的评价体系进行对比，即可得出城市地质环境在的演变模式。关键词：地污染指数；主成分分析；菲克定律；层次分析一、 问题重述为了对某城市城区土壤地质环境进行调查，研究人类活动影响下城市地质环境的演变模式，现将城区按功能划分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，进行研究，不同的区域环境受人类活动影响的程度不同。将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土进行取样、编号，并用gps记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。这样便得到了采样点的位置、海拔高度及其所属功能区；8种主要重金属元素在采样点处的浓度；8种主要重金属元素的背景值。依据以上信息，通过数学建模解决以下问题：（1）给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。（2）通过数据分析，说明重金属污染的主要原因。（3）分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。（4）分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、 问题分析对于问题1，分析题目要求可以看出它有两部分构成：重金属元素在各区域的空间分布分析，以及不同区域重金属污染程度评价。要分析重金属元素在各区域的空间分布，结合题目所给出的离散的样品点，我们想到用matlab拟合出每种重金属元素的分布图，可以直观地看出重金属的分布特征。对于第二部分，要想对不同区域重金属的污染程度进行评价，只要引入一个污染评价指标便可以定量评价每个区域的污染程度。对于问题2，要分析重金属污染的原因，也就是各种金属的来源，需要找出区域的功能与不同污染元素之间的关系。但由于污染元素一共有8种，变量因子过多，且原始数据纷杂，不利于找出内在联系，于是想到用较少的新变量代替原来较多的旧变量，而且保证新变量尽可能多的保留原来变量所反映的信息。这样就可以通过分析少量新变量之间的关系，挖掘出内在联系。主成分分析法便可以实现这一目的。对于问题3，要得到污染源，必须要得到地层中重金属元素的浓度分布，以求浓度最大点。首先我们要明确重金属元素在地层中的扩散规律，然后根据这一规律建立扩散模型。通过查阅大量文献，我们总结出重金属元素在地层中的扩散只有“分子扩散”，而无“对流扩散”，满足菲克定律。这样我们根据特定地层区域内重金属“物料守恒”建立了扩散模型，求得了地层中的浓度分布，进而得到浓度最大点也就是污染源。对于问题4，要研究地质环境演变时，既要研究不同区域地质环境状况的“横向对比”，又要研究不同时间地质环境状况的“纵向对比”，并且要明确地质环境这一概念所包含的因素（重金属污染只是其中一种）。明确了以上两点也就明确了我们研究地质环境演变规律所需要的信息。而在“横向对比”过程中，我们必须要对地质环境状况有一套评价指标，也就是各因素所占权重（权重的变化可以看出演化的过程）及综合评价指标。“纵向对比”只需要对比不同时间同一区域的地质环境状况即可。三、 模型假设重金属污染物浓度在空间的分布连续；重金属污染物在地层中的扩散系数为常数；该城区的重金属污染物分布在短时间内不随时间变化，是稳定状态。四、 模型建立与求解4.1重金属元素分布与不同区域污染程度（问题1）4.1.1各重金属元素在不同区域的分布图要分析8种重金属元素在该城区的空间分布，我们想到了根据给出的样本测量数据，把8种元素分开讨论，剔除样本中的特异点后，对于每一种元素都用matlab软件的“griddata”语句绘制出该种元素在5个区域内的连续分布（在不同区域内的浓度大小）的网格图，进而直观地看出不同区域每一种重金属污染的程度。如下图，是通过matlab绘制的各种重金属元素的污染分布图（绘图源代码见 附录matlab程序第1问）。由于图片所占页面空间太大，这里只贴出as、cd元素的分布图作为代表，剩下6种元素的分布图见附件。对于8种重金属元素的分布图来说，红色表示浓度最大，紫色表示浓度最小。颜色由暖色至冷色表示浓度递减。图1：as元素在各区域的分布从图1中可以直观地看出as元素主要分布在工业区、生活区与绿地周围，而在交通区分布较少。在山区大部分地区分布较少，少部分地方有分布。图2：cd元素在各区域的分布 对于图2，可以直观地看出cd元素主要分布在工业区与交通区周围，次之分布在生活区与绿地，在山区分布较少。另外，由其它6种重金属元素的分布图（见附件图像第1问）可以看出cr元素主要分布在生活区与交通区，在工业区、绿地与山区分布较少。对于cu元素，发现工业区与交通区分布较多，其它地方分布较少。对于hg元素，发现只有少部分生活区与交通区分布较多，其它部分基本上没有受到hg元素的污染。对于ni元素，发现在山区、生活区与工业区分布较多，在交通区与绿地分布较少。对于pb元素，发现在交通区分布最多，在部分生活区与部分工业区也有所分布，而在其它区域都分布较少。对于zn元素，发现在生活区与工业区分布较多，其它区域分布较少。4.1.2 地累积指数法评价不同区域重金属污染程度地累积指数法是一种常用的研究沉积物中重金属污染的定量指标，其计算公式为： （1）式中：是实测元素在地层中的含量；为元素在地层中的背景值；是考虑了各地岩石差异可能会引起的变动而取得变动系数（一般情况下取值为1.5）。【1】得到的地积累指数可分为不同的等级，分级情况与污染程度的对应关系见下表：表1：重金属污染级别、地积累指数（）和分级比较项目污染程度分类清洁轻度污染偏中度污染中度污染偏重污染重污染严重污染5级别0123456对于本题，是题目中所给的各种元素的背景值。是通过题目所给的数据统计得到的各区域中各种重金属元素的平均含量，如下表：表2：各区域中各种重金属元素的平均含量元素平均含量ascd cr cu hg ni pb zn 生活区6.270289.969.0149.4093.0418.3469.10237.0工业区7.251393.153.40127.5642.319.8193.04277.9山区4.044152.338.917.3140.9515.4536.5573.2交通区5.703359.058.0262.02446.117.5963.39242.3绿地区6.104263.142.8330.04109.414.8860.19145.9把代入到式（1）后，得到5块区域内各种重金属元素的地积累指数值（），如下表：表3：各区域内各种重金属元素的地积累指数值ascdcr cuhgnipbzn生活区0.81.1571.1541.9041.410.581.1561.78污染级别12222122工业区1.011.5960.7843.2724.1970.691.5852.01污染级别22145123山区0.1670.2280.320.3910.2260.330.2370.087污染级别11111111交通区0.6631.4650.92.233.6720.521.0321.812污染级别12134122绿地区0.7611.010.461.1861.6450.280.9571.081污染级别12122112通过上表各区域的污染级别可以看出：对于生活区，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均属于轻度污染。另外，为了对比各区域的总污染程度，我们对各个区域内的所有重金属元素的地积累指数求和，得到各区域的总地积累指数（），结果如下表：表4：各区域的总地积累指数（）区域生活区工业区山区交通区绿地区1.21.90.241.50.9通过上表比较可以发现不同区域的污染程度为：工业区交通区生活区绿地区山区。通过对比4.1.1与4.1.2的结论，发现重金属元素的分布情况与各区域内的重金属污染程度相互呼应，吻合的很好。该城区的重金属污染特点，我们可以总结为：该城区的主要污染元素为hg、cu、zn；且污染主要集中在工业区与交通区。4.2主成分分析模型（问题2）要研究重金属污染的主要原因，就要发现不同区域内主要污染物与该区域功能的关联性，进而确定各区域污染原因之间的关联性，以确定整个城区污染的原因。但由于导致每个区域污染的重金属元素众多，没法直接通过重金属元素之间的关系确定各区域污染原因之间的关联性。于是想到用较少的新变量代替原来较多的旧变量，而且保证新变量尽可能多的保留原来变量所反映的信息。主成分分析法便可以实现，主成分分析法是把原来多个变量划为少数几个综合指标的一种统计分析方法。 4.2.1主成分分析模型的建立假定有个地理样本，每个样本共有个变量描述，这样就构成了一个阶的地理数据矩阵： （2）要想从这么多变量的数据中抓住地理事物的内在规律，自然要在维空间中加以考察，这是比较麻烦的。为了克服这一困难，就需要用较少的几个综合指标来代替原来较多的变量指标，这些综合指标（即新变量）选取最简单的形式就是取原来变量指标的线性组合，适当调整组合系数，使新的变量指标之间相互独立且代表性最好。若原来的变量指标为，它们的综合指标新变量指标为，则 （3）在（3）式中，系数由下列原则来决定：（）相互无关；是的一切线性组合中方差最大者；是与不相关的的所有线性组合中方差最大者；是与都不相关的的所有线性组合中方差最大者。这样决定的新变量指标分别称为原变量指标的第一，第二，第主成分。其中，在总方差中占的比例最大，的方差依次递减。在实际问题的分析中，常挑选前几个最大的主成分，这样既减少了变量的数目，又抓住了主要矛盾，简化了变量之间的关系。从以上分析可以看出，找主成分就是确定原来变量()在诸主成分（）上的载荷（），从数学上容易知道，它们分别是的相关矩阵个较大的特征值所对应的特征向量。模型计算步骤：计算相关系数矩阵： （4）为原变量与的相关系数，计算公式为： （5）计算特征值与特征向量。先解特征方程，常用雅克比法求出特征值，并使其按大小顺序排列：。然后，分别求出对于特征值的特征向量，要求，即，其中表示向量的第个分量。计算主成分贡献率及累计贡献率。贡献率： （6）累计贡献率： （7）一般累计贡献率达到80%以上的特征值所对应的为第1、第2、第个主成分。计算主成分载荷 （8）4.2.2模型求解计算相关系数矩阵对所给的数据进行分区域标准化处理，利用式（5）得到5个区域的相关系数矩阵如下：（具体计算源代码见 附录matlab相关系数）生活区：1.00 0.26 0.36 -0.43 -0.41 -0.28 -0.37 -0.50 0.26 1.00 0.29 0.08 -0.20 -0.07 0.05 -0.36 0.36 0.29 1.00 -0.14 -0.06 -0.16 -0.09 -0.10 -0.43 0.08 -0.14 1.00 0.38 0.24 0.53 0.29 -0.41 -0.20 -0.06 0.38 1.00 0.35 0.50 0.40 -0.28 -0.07 -0.16 0.24 0.35 1.00 0.38 0.15 -0.37 0.05 -0.09 0.53 0.50 0.38 1.00 0.20 -0.50 -0.36 -0.10 0.29 0.40 0.15 0.20 1.00 工业区：1.00 0.35 0.51 -0.56 -0.48 -0.44 -0.32 -0.34 0.35 1.00 0.07 -0.21 -0.37 -0.24 -0.29 -0.30 0.51 0.07 1.00 -0.26 -0.30 -0.19 -0.13 -0.15 -0.56 -0.21 -0.26 1.00 0.33 0.38 0.15 0.18 -0.48 -0.37 -0.30 0.33 1.00 0.54 0.57 0.53 -0.44 -0.24 -0.19 0.38 0.54 1.00 0.92 0.90 -0.32 -0.29 -0.13 0.15 0.57 0.92 1.00 0.98 -0.34 -0.30 -0.15 0.18 0.53 0.90 0.98 1.00 山区：1.00 0.03 0.03 -0.17 -0.07 -0.01 -0.11 0.12 0.03 1.00 0.61 -0.18 0.28 -0.23 -0.30 -0.09 0.03 0.61 1.00 -0.33 0.23 -0.18 -0.37 -0.05 -0.17 -0.18 -0.33 1.00 -0.29 0.11 0.53 0.08 -0.07 0.28 0.23 -0.29 1.00 0.07 0.09 0.25 -0.01 -0.23 -0.18 0.11 0.07 1.00 0.36 -0.01 -0.11 -0.30 -0.37 0.53 0.09 0.36 1.00 0.51 0.12 -0.09 -0.05 0.08 0.25 -0.01 0.51 1.00 交通区：1.00 0.51 0.38 -0.11 -0.24 -0.27 -0.28 -0.08 0.51 1.00 0.32 -0.09 -0.12 -0.10 -0.09 -0.12 0.38 0.32 1.00 -0.11 -0.13 -0.17 -0.16 -0.07 -0.11 -0.09 -0.11 1.00 0.12 0.14 0.09 0.00 -0.24 -0.12 -0.13 0.12 1.00 0.37 0.42 0.21 -0.27 -0.10 -0.17 0.14 0.37 1.00 0.89 0.01 -0.28 -0.09 -0.16 0.09 0.42 0.89 1.00 0.03 -0.08 -0.12 -0.07 0.00 0.21 0.01 0.03 1.00 绿地区：1.00 0.51 0.04 -0.40 -0.10 -0.26 -0.11 -0.33 0.51 1.00 0.03 0.12 0.05 0.32 0.15 -0.21 0.04 0.03 1.00 0.17 0.01 0.01 -0.07 -0.01 -0.40 0.12 0.17 1.00 0.20 0.64 0.10 0.12 -0.10 0.05 0.01 0.20 1.00 0.50 0.54 -0.01 -0.26 0.32 0.01 0.64 0.50 1.00 0.36 -0.03 -0.11 0.15 -0.07 0.10 0.54 0.36 1.00 0.13 -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 表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个主成分即可；对于工业区只需求出前2个主成分即可；对于山区，只需求出前3个主成分即可；对于交通区，只需求出前3个主成分即可；对于绿地区，只需求出前3个主成分即可。求主成分载荷对于5个区域，利用式（8），分别求出原变量在各个选定的较大主成分上的载荷。如下表：（具体算法见 附录matlab第2问；详细数据见 附录附表第2问）表10：各区域主成分载荷原变量主成分生活区工业区山区z1z2z3z1z2z1z2z3x1（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 原变量主成分交通区绿地区z1z2z3z1z2z3x1（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 （说明：由于计算载荷时我们只精确到小数点后一位，所以，有的载荷值经过“四舍五入”，为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重金属元素在地层中的扩散特点由于我们研究表层土壤中重金属元素的分布，所以其扩散规律不涉及到随流体流动而产生的“对流扩散”，也就是说其扩散情况只有“分子扩散”。分子扩散完全是因为组分分布不均匀，即在空间中存在浓度梯度，导致这些组分依靠分子热运动从高浓度带扩散到低浓度带，最后趋于平衡。根据描述物质“分子扩散”现象的规律菲克扩散定律，在单位时间内通过垂直于扩散方向的单位截面积的扩散物质流量（称为扩散通量）与浓度沿该截面外法向的方向导数成正比，也就是说，浓度沿截面法向单位距离的变化率越大，扩散通量越大。4.3.2扩散模型的建立在地层中，若各点金属元素浓度不同，就会有重金属元素由高浓度带向低浓度带扩散。确定污染源的位置，就是要找到重金属元素的浓度最高点，故问题归结为研究地层中重金属元素的浓度分布，即根据重金属元素在地层中的扩散特点建立浓度关于位置的函数关系。 假设地层中任意一点在时刻t的重金属污染物浓度为。假设在时间区间内，流过点处面积为的曲面重金属流量为，则由菲克定律的得： （9） 其中比例系数称为重金属的扩散系数（大于0），在一般情况下都把它作为常数处理。是曲面的外法线单位向量。负号表示流量方向与浓度梯度方向相反。在地层中任取一封闭曲面，记它所包围的区域为。于是，从时刻t1到时刻t2，通过曲面流出该区域的重金属污染物质量为： （10）若内有污染源不断产生污染物，设污染源密度为（在时刻t于任一点处单位时间、单位体积所产生的污染物质量），则在时刻t1到t2内产生的污染物质量为： （11）由于物质扩散，在时刻t1到t2内，中的污染物质量的改变量为： （12）由物料守恒定律有：，即 （13）式（13）就是我们建立的重金属污染物在地层中的扩散模型。根据这个方程便可以求得。由于我们研究的是该城区地层在这一确定时刻的污染情况，所以与时间无关，即的值为0。另外，在研究该问题的时候，我们在短时间内不考虑新的污染物的注入，所以。而又为常数，这样我们的模型可以简化为即： （14）由于方程（14）是偏微分方程,只有在问题具有高度对称的情况下,才能求出解析解,但我们研究的重金属元素在地层中的扩散受到地貌、浓度分布、以及其他复杂因素的影响，显然不能符合高度对称这一条件。对于这些情况,只能寻求它的数值解。【2】4.3.3模型求解（拉普拉斯方程的数值求解）有限差分法是求解拉普拉斯数值解的一种常用算法【3】，其思想是用差分代替导数，用网格将求解区域覆盖，对于平面拉普拉斯方程，第i行，第j列小格的值由拉普拉斯五点差分格式给出： （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种重金属元素的污染源位置坐标，如下表：表11：8种重元素的污染源位置坐标asx/my/mz/mg/ghgx/my/mz/mng/g118234.36 10063.09 16.05 30.02 16946.42 7267.79 25.06 1860.30 227496.26 12299.33 72.86 18.00 28683.03 12112.98 4.76 1756.89 312735.11 2981.66 22.45 23.69 322286.44 10435.80 39.60 1730.44 44052.08 7826.85 10.06 22.50 cdx/my/mz/mng/gcux/my/mz/mg/g121418.14 11367.57 39.89 1600.14 11736.61 2795.30 12.41 303.14 217655.49 3727.07 31.04 1232.87 210419.64 8013.20 52.52 286.39 pbx/my/mz/mg/gnix/my/mz/mg/g14777.00 4897.00 7.16 472.48 122286.44 12299.33 83.32 62.56 （说明：每种元素的含量均有多个峰值，也就是有多个污染源。序号1,2表示污染源的编号）4.4城市地质环境演变模式分析（问题4）4.4.1所建模型局限性在前3问建立的所有模型中，城区地层评价指标都针对于城区以及不同区域的重金属污染物而建立的，这些评价指标都只能评价重金属污染状况，而重金属污染状况只是地质环境的很小的一个方面，所以我们前面所建立的模型虽然能比较准确的评价地层重金属污染的情况，但却不足以全面地反应地层的地质环境状况，模型反应的信息太过于局限；要研究该城区地质环境的演化模式，必须要对该城区的地质状况进行纵向对比，也就是不同时间的地质环境状况对比。而我们在前3问建立的模型，都是反映该城区地层在这一特定时间的地质环境状况，不包含时间这一参数，也就不能对地质环境状况进行纵向对比。反映不出地质环境状况随时间的变化情况。4.4.2研究地质环境状况演变所需要的信息地质环境指由岩石圈、水圈和大气圈组成的环境系统。城市发展对环境的改造与扰动主要通过三个方面, 即: 大规模的地面和地下工程建设；