对城市表层土壤重金属污染的探讨(1)

1、2011高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名

2、号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：2011高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：对城市表层土壤重金属污染的探讨摘要本文主要针对城市表层土壤重金属污染的问题进行了研究，首先通过MATLAB进行插值，得到了每种元素的浓度的分布图，之后采用三种不同的

3、判断方法进行分析，并得到如下结论：1.用内梅罗单因子指数法分析得到，工业区和交通区的Cd为轻度污染，工业区的Cu达到重度污染的水平，交通区的Cu和Hg达到轻度污染的水平，同时，生活区、工业区、交通区的Zn达到了轻度污染的水平。Cd，Cu，Zn的污染程度总体已经达到轻度污染水平。2.用内梅罗综合指数法分析得到工业区、交通区污染指数达到轻度污染水平3.用负荷指数法分析可知，五个不同城市区域的污染程度均已达到轻度污染，工业区的污染水平已经达到了中度污染。以上三种分析方法从不同侧重点分析污染特征，各有利弊，综合使用可以更全面分析，结果虽有一定差距，但可以接受。第二问是针对重金属污染的主要原因进行说明，

4、本文中选取相关矩阵和变异系数两个量来进行衡量，通过相关矩阵探究重金属污染元素两两之间的相关性，相关性较强的元素可以综合分析，由变异系数确定重金属污染分布的差异性，为了说明问题更加准确，本文对五种不同城市区域分别研究，最后研究城市总体的情况，并由此得到结论。第三问是要确定污染源所在位置，由模型假设可得，重金属污染传播特征符合高斯扩散模型。依据推导出的方程，对等式两边取对数，得到二次二项式之后对其进行多元非线性拟合，求出污染源的坐标位置。在该过程的MATLAB模拟时，采用了“环带辐射搜索法”对所需的点进行搜索并拟合曲线，得到8种元素共19个污染点（有的污染点之间距离很近，可以合为一个小区域），总体

5、来看，可分为六个大污染区域和四个小污染点，经过拟合优度检验，结果较为精确。第四问对第三问从原理角度进一步分析，对模型进行完善，根据土壤动力学理论的相关知识，列出重金属离子扩散的偏微分方程。针对偏微分方程的求解，使用交替方向隐式格式（ADI）的差分方法并设计算法对解决。从而建立研究城市地质环境的演变规律的模型。关键词：内梅罗指数法 负荷指数法 相关矩阵 变异系数 高斯扩散模型 交替隐式差分方法一、问题重述人类活动对于城市土壤环境的影响越来越受到人们的关注，按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度

6、不同。现对某城市城区土壤地质环境进行调查。根据附录中所给出的信息对以下问题进行分析：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、问题分析对于问题一，采用插值的方法在三围和二维图中直接将重金属浓度的变化趋势变现出来即可，制图软件使用MATLAB。对于土壤重金属污染程度分析，我们对不同区域的总和污染程度进行探讨，

7、在这里本文采用了单因子指数法、改进的内梅罗指数法和污染负荷指数法三种方法综合判定，得到各个区的污染程度。三、模型假设1.土壤各项同性均匀2.同种金属离子在不同方向上扩散系数相同3.扩散（传播）过程中忽略重力的影响4.忽略高度对于扩散（传播）的影响5.忽略机械弥散对于扩散（传播）的影响四、定义与符号说明内梅罗单因子指数 内梅罗综合指数 负荷指数 变异系数 流量 表示污染源释放重金属离子的总量 相关指数 土壤容重 空间域 体积 相关指数 土壤溶质对流通量 土壤水分通量 土壤溶质的浓度 土壤平均空隙流速 土壤体积含水率 溶质分子的扩散通量 土壤平均空隙流速 扩散系数 五、模型的建立与求解5.1模型一

8、的建立与求解5.1.1模型一的建立第一问第一部分设计使用MATLAB对于图像的处理，对数据的水平分布散点进行插值，绘制出不同元素的浓度等高线图和三维立体图。土壤重金属污染的主要污染元素因为城市区域职能的不同而不同，污染浓度也必定不同，因此由内梅罗单因子指数法就此问题展开讨论，从而得到每一不同的城市区域的主要污染元素是非常必要的。内梅罗单因子指数法在此，本文首先利用内梅罗单因子指数法，通过单因子评价可以确定每种重金属污染物的危害程度，以污染指数表示，通过每种重金属含量实测值和评价标准相比计算。污染指数用表示，。但是，由于重金属离子本身在自然界中就会存在，因此要对该公式进行修正，即分子分母同时减去

9、自然背景值，修正之后公式为，据此式得出五种不同城市区域内的八种金属分别的污染程度。改进的内梅罗综合指数法但是，仅仅由单因子指数法判定不能反映土壤污染的整体情况，需要用一个方式建立一个合适的评价系数，通过这个值得高低综合判断污染程度的大小，本文在此采用了改进的内梅罗综合指数法，设立，即采样点的综合污染指数，令，为采样点重金属污染指数中的最大值，为单因子指数平均值，在传统的内梅罗指数法中，直接取算术平均值，但是本文在这里对此进行修正，本文认为，如果单纯考虑算术平均数，那么就默认所有单位计量的重金属对于动植物和人类的危害是相同的，但这个假设是不成立的，重金属的生理毒性有高低，本文按照其对于环境的影响

10、程度，对此采用加权平均。即本文中，其中的的取值如下：HgPbCdAsZnCuCrNi类别权重33332222将相关数据代入即可得到结果。由上述方法可得综合污染指数，但是该方法仍然存在缺陷，即该方法污染过分突出了含量最大的污染物的作用，在评价时，人为地夸大或缩小了一些因子的影响，使其灵敏度不够高，为了避免这个问题，本文引入了另一种常见对重金属污染的评价方法负荷指数法。负荷指数法该种方法是从事重金属污染水平的分级研究中提出来的一种评价方法，该指数由评价区域所包含的多种重金属成分共同构成，并使用求积的统计法，通过这种方法能对整个区域各个点位各个重金属进行定量评价。并对各点的污染程度进行分析，并对各点

11、的污染程度进行评级，以及重金属在时间空间上的变化趋势，该种方法能够避免污染指数叠加造成的人为夸大或缩小一部分因子对其影响的情况。设元素的最高污染指数为，为元素的实测含量，为元素的背景值。某一点的污染指数()为：，为评价点个数，代入数据求解即可得到结论。5.1.2模型一的求解如下图所示图像为八种重金属污染物的分布，三维图Z轴表示浓度，二维图中采用等高线的方法，由不同的颜色代表不同的浓度，暖色调表示浓度高，冷色调表示浓度低。之所以采用两种示意图，是因为考虑到两种示意图各自的优缺点，三维图显示较为直观，但是该方法不易读出具体数据，二维图虽没有三维图直观，但是却能具体反映出浓度的变化情况。（二维图详见

12、附录）如下图所示为As元素的浓度分布 Cd元素分布如下图Cr元素分布如下图 Cu元素分布如下图Hg元素分布如下图 Ni元素分布如下图Pb元素分布如下图 Zn元素分布如下图内梅罗单因子指数法结果由MATLAB分析单因子指数法得到的结果如下：As (g/g)Cd (ug/g)Cr (g/g)Cu (g/g)Hg (ug/g)Ni (g/g)Pb (g/g)Zn (g/g)1区0.23430.940950.17360.9837820.2190210.2181330.1740021.2825092区0.32031.547710.102333.1069522.2919070.2711790.283291

13、1.5948663区0.03900.131290.036350.1118820.0224760.1138550.0253690.032784区0.18491.353030.123531.3319271.5540470.1919530.1485581.3271365区0.23370.885550.05770.4617310.3018560.1079320.1356560.650704总0.18251.014120.102781.13630.9989060.1791280.140371.00918该方法的判断标准如下表：土壤综合污染程度分析标准值污染水平非污染轻度污染中度污染重度污染由该标准可以判

14、断，2区（工业区）和4区（交通区）的Cd超标，为轻度污染，2区（工业区）的Cu超标严重，达到重度污染的水平，同样，4区(交通区)的Cu也超过标准，达到轻度污染的水平，4区（交通区）的Hg也超过标准，达到轻度污染的水平，同时，1区（生活区）、2区（工业区）、4区（交通区）的锌都超标，达到了轻度污染的水平。Cd，Cu，Zn的污染程度总体已经达到轻度污染水平，Hg的污染水平已经非常接近轻度污染，因此该城市重金属污染状况应该引起重视。内梅罗综合指数法结果由内梅罗综合指数法，得到的结果如下：五区污染程度内梅罗综合指数12345总0.8756561.4975410.2599631.0411480.7127

15、740.862572内梅罗综合指数的参考标准如下：土壤综合污染程度分级标准综合污染指数污染程度安全警戒线轻度污染中度污染重度污染由此标准，可知工业区和交通区的污染程度已经达到轻度污染程度，居民区和公园区的污染指数也已经达到了警戒线，只有山区的污染指数较低，在安全范围之内，总体上看也是安全的，但是也需要引起重视，防止污染进一步恶化。负荷指数法结果背景值选取题目所给背景值，得到结果如下：污染指数法（题目背景值）12345整体1.8335612.3484481.0602211.9209211.577991.689651土壤综合污染程度分级标准污染指数污染程度安全轻度污染中度污染重度污染得到的结果显示

16、，各个区域都是轻度污染级别，工业区的污染情况已经达到中度污染，这与用内梅罗综合指数法得到的结论有一定差距。但是由于评价侧重点不同，本文认为，这种差异是可以接受的。结论分析 三种表达方式，第一种具体表现了每一种污染物的状况，第二种和第三种都是衡量一个地区污染指数的表达方式，两种方法各有侧重，第一种方法突出了影响最重要的一个因素，但是此种方法弊端也同样在此，这样可能会把影响因子人为的夸大或缩小，使得结果可能相对不够准确，负荷指数法的优点在于避免指数加和关系对于结果的歪曲，并且能够对任意给定的区域作出定量判断，问题在于没有重视不同污染物的背景差别，权重认为一致，这是该方法的问题所在。5.2模型二的建

17、立和求解5.2.1模型二的建立要考虑重金属污染的原因，那就不能孤立看待每一种重金属元素，本文认为，要通过相关系数求出重金属元素之间的关系，相关系数越高的两种重金属，说明他们关联越大，因此本问先求出相关系数矩阵，再求出变异系数，变异系数计算公式为，目的是观察观测值的变异程度，据此判断人类活动对于重金属元素在土壤中含量的影响。5.2.2 模型二的求解1.求变异系数并逐元素逐区分析求出变异系数，由此来判断人类活动对于重金属含量的影响，使用SPSS的描述统计量，得到了如下结果：描述统计量金属N极小值极大值均值标准差As3191.6130.135.67653.02429Cd31940.001619.80

18、302.3962224.98760Cr31915.32920.8453.509770.00179Cu3192.292528.4855.0167162.91510Hg3198.5716000.00299.71131629.53978Ni3194.27142.5017.26189.94142Pb31919.68472.4861.740950.05776Zn31932.863760.82201.2026339.23254由此我们得到变异系数如下表：AsCdCrCuHgNiPbZn0.5340.7441.3082.9615.4370.5760.8111.686此数据与上海师范大学硕士生尹骏的硕士毕业论

19、文上海市城郊梯度土壤重金属空间分布特征和评价研究一文中得到的数据相去甚远，该文章中得到的上海市闵行区土壤重金属浓度变异系数最大的元素锌，变异系数也仅仅0.59，此文中的数据与之相比，差异性非常大，因此本文认为有必要对每个元素分城市区域详细的进行讨论，本文首先将对砷元素展开讨论。如下图所示（描述统计量见附录）。As的变异系数（分区）As1As2As3As4As50.340.5850.4450.5670.323该数据组为砷元素的分布情况，由此可以看出，砷分布只在居民区和公园绿地区分布较为平均，在其他地区，分布差异均较大，本文认为，这种现象的发生与砷元素本身的来源有关。除了土壤中本身含有的砷元素外，

20、砷元素的来源主要是造纸工厂，工厂的废气排放和含砷废水排放，另外矿山采矿（主要是锌矿）也会产生大量砷元素，这些最后都会进入土壤，这也就是山区和工厂区变异系数大的原因，即这类工厂周围砷元素明显富集，其他地区影响相对较低。Cd的变异系数（分区）Cd1Cd2Cd3Cd4Cd50.6330.6040.5150.6760.841 镉的整体分布都较为平均，因为其在工业中用途相当广泛，因此，工业区污染虽重但变异系数却较低。生活区和交通区由于汽车尾气的影响使得其分布呈现不均匀的态势，即车辆通行多的地区、堵车严重的地区镉含量要高于其他地区。Cr变异系数（分区）Cr1Cr2Cr3Cr4Cr51.5630.8240.

21、6311.4060.340铬元素的分布存在较大差异性，在生活区和交通区差异尤其大，最大可达60多倍，原因在于，铬主要来自工厂排放废水废气、石油化工产品燃烧，汽车尾气中有铬的存在，而在车流量大易堵车之处，土壤中铬含量将会明显的高，这个事经过其他实验数据验证的，因此可以解释生活区和交通区铬元素分布不均的现象。而工厂区金属加工、电镀、制革等行业中为了防止工业生产过程中循环水对设备的腐蚀常加入铬酸盐，铬酸盐的广泛的采用使得工厂区的铬污染虽然较其他区域重，但是其用量的原因使得该地区铬变异系数却并不大。Cu的变异系数（分区）Cu1Cu2Cu3Cu4Cu50.9543.2530.6201.9320.751C

22、u的来源也同样分为两部分，工业生产废渣废气和汽车尾气，少数工厂会产生大量含铜的工业三废，使得这部分地区土壤中Cu含量大大高于其他地区，这也就解释了工厂区的Cu变异系数大的原因，在交通区，由于汽车尾气中含有一定量的铜随汽车尾气的排放而进入大气，随后以降水的形式进入土壤中，使得这一地区的土壤中Cu含量明显偏高，这也与数据中显示的情况吻合，本文认为，在交通要道和工业生产中，大量使用Cu元素的工厂附近，其土壤中Cu离子浓度都要明显高于其他地区。Hg的变异系数（分区）Hg1Hg2Hg3Hg4Hg50.9543.4940.6804.8791.950据资料表明，汞的产生主要来自化石燃料，特别是煤炭的燃烧，这

23、是汞元素最主要的排放途径，尤其是燃煤电厂，据统计其排放总量占到总量的25%以上。排放到大气中的汞元素最终大部分进入土壤沉积，而且汞元素的特点是极易造成跨区域污染，污染与排污点距离呈现反相关，所以这也就解释了各个城市区域汞元素变异系数均较大的原因。并且汞污染在工业及交通区污染较为严重。土壤中镍元素来源也较为广泛，但是从该地区情况可以看出，镍元素污染并不严重，在此略去不作分析。土壤中铅含量与汽车尾气，工业三废的使用均有关，情况较为复杂，但本城市铅污染并不严重，可能由于含铅工业没有在该地大量发展，城市交通区铅污染指数也并不高，因此对其变异系数也不做分析。Zn的变异系数（分区）Zn1Zn2Zn3Zn4

24、Zn51.8711.2620.4221.5841.497该城市的Zn污染较为严重，生活区，工业区，交通区污染指数均超过标准达到轻度污染级别，并且同一区域内差别也很大，在生活区和交通区，本文认为产生这种现象的原因是汽车尾气的排放和大气中含锌降雨降落的结果，因此Zn污染严重程度与距离污染源距离反相关，这就解释了Zn的变异系数除了山区都很大的原因，并且Zn污染在该市较为严重，应该对污染较重的生活区采取一定的措施进行保护。2.分析土壤污染各个元素间相关系数相关矩阵 AsCdCrCuHgNiPbZnAs1.000.255.189.160.064.317.290.247Cd.2551.000.352.39

25、7.265.329.660.431Cr.189.3521.000.532.103.716.383.424Cu.160.397.5321.000.417.495.520.387Hg.064.265.103.4171.000.103.298.196Ni.317.329.716.495.1031.000.307.436Pb.290.660.383.520.298.3071.000.494Zn.247.431.424.387.196.436.4941.000从总体上看，Ni和Cr相关性最高，达到0.716，Pb和Cd相关性次之，0.660，Pb和Cu相关性也在0.5以上，说明在总体上这几个元素之间有很

26、大关联。主要原因在于工厂中两者的联合使用和汽车尾气的影响。下面本文将对每一部分具体分析。对于1区（生活区）的相关系数矩阵如下：相关矩阵As1Cd1Cr1Cu1Hg1Ni1Pb1Zn1As110.3810.2380.5310.2930.6050.45-0.017Cd10.38110.3490.4990.3970.2830.8020.346Cr10.2380.34910.3760.150.5270.4160.412Cu10.5310.4990.37610.1980.4340.5020.238Hg10.2930.3970.150.19810.2110.340.242Ni10.6050.2830.52

27、70.4340.21110.30.334Pb10.450.8020.4160.5020.340.310.328Zn1-0.0170.3460.4120.2380.2420.3340.3281此部分Pb与Cd相关度最高，与总表吻合，Ni和Cr的相关度与总体相比明显偏低，Ni和As相关度在该地区明显高于总表。Pb和Cu的相关度水平与总表接近，本文认为Pb与Cd的来源在该区域均为主要来自汽车尾气，其中一部分直接沉积到土壤，另一部分被大气吸收后又随降雨等降回地面被土壤吸收，因此两者相关度最高。Ni与Cr在生活中的相关度低本文认为其主要与大气和地下水有关，即生活区的由人为排放的Ni和Cr很少，主要是由大

28、气运动造成的，即将工业区和交通区排放的这两种元素随降水等形式降落到土壤内，并且地下水流动，将这些重金属离子带入到1区土壤中。Pb和Cu的相关度与总表接近，也主要是汽车尾气的影响，至于Ni和As的相关度提高很大，可能由于工业区之间相关系数很大造成。对于2区（工业区）的相关系数矩阵如下：As2Cd2Cr2Cu2Hg2Ni2Pb2Zn2As21.000.394.407.177.204.701.440.543Cd2.3941.000.554.558.529.519.839.764Cr2.407.5541.000.916.899.693.685.704Cu2.177.558.9161.000.984.4

29、93.670.627Hg2.204.529.899.9841.000.471.616.597Ni2.701.519.693.493.4711.000.584.629Pb2.440.839.685.670.616.5841.000.756Zn2.543.764.704.627.597.629.7561.000工业区各元素之间的相关系数比其他地区明显增大，这与该地区的工业生产常常会同时排放多种元素，造成这一地区元素间联系较为紧密。对于3区（山区）的相关系数矩阵如下：相关矩阵 As3Cd3Cr3Cu3Hg3Ni3Pb3Zn3As31.000-.425-.013.463.057-.033-.259-.

30、355Cd3-.4251.000-.041-.083.131-.079.757.591Cr3-.013-.0411.000.197-.079.934.070.582Cu3.463-.083.1971.000.552.180.057.039Hg3.057.131-.079.5521.000-.144.117.104Ni3-.033-.079.934.180-.1441.000-.043.531Pb3-.259.757.070.057.117-.0431.000.568Zn3-.355.591.582.039.104.531.5681.000由于山区受人类活动影响较小，因此该地区的土壤元素分布较多

31、的呈现自然状态，当然受到大气环流和降水以及地下水等影响，大气中相关系数较高的元素也会同时降落到地面，被土壤吸收之后导致该部分区域几种元素之间相关系数增大，Ni和Cr呈现非常明显的相关性，并高出其他地区，说明该地区自然分布和认为因素都起到重要作用。对于4区（交通区）的相关系数矩阵如下：相关矩阵 As4Cd4Cr4Cu4Hg4Ni4Pb4Zn4As41.000.117.138.090-.005.226.057.186Cd4.1171.000.372.422.211.348.613.292Cr4.138.3721.000.894.011.869.427.395Cu4.090.422.8941.000

32、.032.886.505.431Hg4-.005.211.011.0321.000.039.265.118Ni4.226.348.869.886.0391.000.394.502Pb4.057.613.427.505.265.3941.000.481Zn4.186.292.395.431.118.502.4811.000该部分污染物的产生原因相对简单，主要为煤以及汽油天然气等燃烧，其中主要的是汽车尾气，由图中可以看出汽车尾气中富含的元素之间彼此均呈现一定程度的相关性，Ni，Cr，Cu彼此之间相关性很高。Pb与上述元素之间相关性也不可忽略。对于5区（公园绿地区）的相关系数矩阵如下：相关矩阵 As

33、5Cd5Cr5Cu5Hg5Ni5Pb5Zn5As51.000.138.410.030-.046.371.077.174Cd5.1381.000.223.311.078.079.415.180Cr5.410.2231.000.295-.002.612.265.180Cu5.030.311.2951.000.106.177.664.351Hg5-.046.078-.002.1061.000.042.146.047Ni5.371.079.612.177.0421.000.062.420Pb5.077.415.265.664.146.0621.000.362Zn5.174.180.180.351.04

34、7.420.3621.000Ni和Cr，Pb和Cu呈现一定的线性关系，这一部分来源本文认为也是以汽车尾气为主，这与4区情况吻合。3.结论通过对相关性和变异系数的分析，本文得到如下结论，即土壤重金属污染主要是由人为因素而非自然因素造成的，每一部分的相关性和总体的相关性呈现较大差异，工业区的金属元素相关性明显高于其他地区，生活区，交通区的汽车尾气排放对于该地区金属污染结果影响较大。工业区主要来自三废排放，公园区主要来自大气环流和地下水的影响，山区主要受自然因素、大气环流和地下水影响。4.建议生活区的有些地方，某些污染指数超标，比如生活区有的地方As元素超标，甚至超过了大部分的工业区的指数，这可能与

35、该地区的大气环流特性、地下水特性、地形特性均有关，应加以妥善治理，如不能解决，则要尽快远离此地。对于交通区的治理是一个相对麻烦的问题，因为世界各国各城市交通区污染都较为严重，纽约、伦敦、东京也是如此，越大的城市治理起来就越麻烦，除采用无铅汽油等防止恶化的方法之外，种一些可降解重金属的植物，对一些经常堵车的交通区采取一定措施保证畅通，都是现阶段可行的办法。对于工业区，由于这是污染最主要的地区，因此限制不合格的三废排放是最重要的，对不合乎规定的要及时停产，另外，工厂一定不能修在城市的上风区和河流上游，带来的危害特别大，对其他城市区域土壤均带来严重影响。工业区与生活区尽量不要相连，就算相连，生活区也

36、不能建在工业区的下风区和下游。5.3模型三的建立和求解1模型三的建立（1）传播特征分析本文对污染模型进行探究后发现，在理想状态下，金属离子的传播服从扩散定律，即单位时间通过单位法向面积的流量与它的浓度梯度成正比，为了简化模型，本文在此忽略重力对于传播的影响，并且认为离子的扩散系数是处处相同且恒定不变的，并且土壤是均匀的。设开始扩散的时刻，时刻无穷空间中任意一点坐标为的浓度记为,根据假设2，单位时间通过单位法向面积的流量与浓度梯度成正比，则可得到：是扩散系数，表示梯度，负号表示由浓度高向浓度低的地方扩散。考察空间域，其体积为，包围的曲面为，为一规则的球面（计算时由于球面比半球面简单并且对于最后结

37、果没有影响，因此本文计算方式上还是按照球面来计算，但只能取下半球面，因为从地表只能向下扩散）外法线向量为。则在内通过的流量为：内重金属离子的增量为：从污染源释放的重金属离子物质的总量为：根据质量守恒定律，单位时间内通过所选曲面的向外扩散的重金属元素与曲面内重金属元素物质增量之和，等于污染源在单位时间内向外放出的污染物。则：即：又根据曲面积分的Gauss公式：因此：由以上公式并利用积分中值定理得： 这是无界区域的抛物线型偏微分方程，根据假设，初始条件为作用在坐标的函数，记作 表示污染源释放重金属离子的总量方程的解为： （2）模型的数学变换由于根据第一问的数据分析可知在海拔较高处污染很小，而在海拔

38、低处污染的程度与海拔关系亦不明显，因此本文中忽略海拔对于结果的影响，将方程简化为二维方程，即对该方程两边取对数，得到如下结果： 令，为未知常量，可得如下方程：中，可以直接由数据得到，在此本文中对曲线进行拟合，通过拟合得到的曲线求出其中的即可。（3）算法的实现此算法首先要解决的是筛选污染源周围主要受到该污染点影响的点，在此，本文认为，极值点的大体范围已经可以通过第一问的图像和数据得到，然后采用“环带辐射搜索法”，即以假定的极值点周围以一定的半径开始向外搜索，搜索时，以环带形式向外辐射，比较最外环带内所有点与次外环带内的重金属污染浓度的最大值之间的关系，如果最外环带中的点存在比次外环带的点的浓度的

39、最大值大，搜索即可停止，取最外环带内的圆的点，如果最外环带内所有的点都比次外环带低，那么就再向外以上述方式进行搜索，直到发现这样的点或者将搜索范围扩展到边界为止。之后通过MATLAB拟合曲线即可得到结论。2.模型的求解按上述方法，本文从初始点进行搜索，每次向外拓展250米，即最大圆的半径比次大圆大250米，得到每种污染源的位置如下表所示：求出污染源位置搜索初始点金属xyx(m)y(m)海拔(m)1As18041.999789.7911813410046411As11761.533071.307126963024271As5044.8796760.2194742729392Cd21441.951

40、0778.212143911383452Cd6532.48710366.7863951044342Cd2020.1982735.7911647272863Cr3852.0355438.7283299601844Cu2220.5233234.3282383369275Hg1973.6953148.31227082295225Hg12925.52167.742136942357335Hg15327.878298.871152489106166Ni3959.7775624.9273299601846Ni21415.8211831.762219312185797Pb4356.1325467.37547

41、77489787Pb1928.3323424.0291647272868Zn14338.989098.664137979621188Zn8883.9154873.06693284311248Zn3376.7595748.2633299601848Zn12508.232902.484126963024273.对结果的分析通过图示比较局部搜索浓度最高值和理论污染源的位置，对其误差情况进行粗略分析。可见理论值与实际求得的值很接近，有明显的6个部分呈现集中分布的区域和4个分散的区域。工厂区的污染源最多，山区和生活区也有污染源分布，在此本文认为污染源主要来自工厂生产，这与实际情况相吻合。对误差的分析拟合

42、优度检验拟合的优劣可以用误差的平方和来衡量。假设实际测得的值为 ，其平均值为，根据拟合曲线所求得的理论值为 ，则可得其误差平方和为，均方差为，若误差平方和与均方差的比值越小，则说明,实际观察值与估计值越接近，拟合的越好。为此可定义相关指数，其计算公式为： ,对每种元素任取一个点进行分析，结果如下： 元素编号R2（相关指数）Ni1350.9755Cd950.9734Cr220.946Cu80.9317Hg90.8747Ni220.9562Pb160.9754Zn610.9916各组的值都接近于1，据此本文认为经过该方法检验，拟合效果较好，确定污染源位置较为精确。5.4模型四的建立和求解1.模型三

43、的优缺点评价 模型三建立的较为简洁且相对准确，使用拟合的方法求出了坐标，但是模型三中的问题也很突出，没有对土壤中重金属离子传播的机理进行深入的分析，建立模型时过于理想化，如认为传播之后的分布曲线是高斯分布，忽略海拔对于结果的影响，对于降雨等因素都予以忽略。2.模型四的建立 考虑到模型三的缺点，模型四的建立基于土壤动力学原理，对该过程中的土壤浓度的变化从机理上进行分析，从而在3问的基础上继续予以深化。土壤溶质（包括污染物）在土壤中随着土壤水分的运移而运动，同时土壤溶质在运移的过程中还会受到土壤基质、土壤溶质本身的物理化学性质以及土壤中其他化学成分的影响。土壤溶质运移的物理机理主要包括对流、扩散和

44、机械弥散。其中，为土壤溶质对流通量，表示土壤水分通量，表示土壤溶质的浓度。，是土壤平均空隙流速，是土壤体积含水率。由以上两式可知，。 扩散的表达式为，其中表示溶质分子的扩散通量，表示扩散系数。由于机械弥散的问题对结果影响较小，在此忽略。列出方程如下：在二维状态下，可得方程如下：（）（为土壤容重）3.模型四所需变量土壤容重，土壤体积含水率，土壤平均空隙流速，扩散系数。 4.模型四的求解二元偏微分方程正常情况下无解析解，因此采用有限差分的方法求其数值解，这里采用P-R格式的差分方法，即实质是在时间的第层和第层之间插入一个第层作为过渡层，即一个时间段分两步进行，在前的步长内对方向采用隐式的差分格式，

45、另一个方向采用显式的差分格式，在后的时间内正好相反，即对方向采取显式差分格式，对方向采用隐式的差分格式。用和分别和方向上的空间步长，表示时间步长，按照交替隐式的差分方法，对方程展开并建立差分方程，经整理后可得两式： 该计算方式以由上述两式组成的方程组即是原式的交替方向隐式格式的差分方程，其计算是根据初始条件即的时刻，可由前一式计算出1/2时刻的浓度值，因为这时前一式的右边各项均是已知的，所以对左边各项可采取固定 的方法，对方向按照一维格式的隐式方法求解，即追赶法，因为在这一步计算中仅计算方向也就是方向，所以又称为向方向追赶。然后由后一式计算出一个时间步长的浓度，再把代入前式计算出的浓度值，以此

46、类推，反复迭代即可得到各个时刻的浓度。在收集到土壤容重，土壤体积含水率，土壤平均空隙流速，扩散系数后，各项的系数便为常数，上述方程可改进为下式：其中和表示为的系数，均为常数，为以为自变量的函数，在已知时各点的浓度可通过以上模型使用MATLAB进行以时间为步长的仿真，从而求得时刻各点的浓度。由此方法更好地研究城市地质环境的演变模式。六、模型的评价及优化对本文中出现的模型，其实有很多可以改进的地方，模型二的分析已经相对到位，对数据具体的分析已经十分详细，但是这里没有对8个因子进行因子分析，没有对影响最大的元素究竟是哪种进行探究，这是二问模型的不足之处，对二问的优化就是对八种元素进行因子分析，求出影

47、响最大的几种主成分元素。对于模型三，曲线拟合的方法在此是最适合的，由于数据量的缺乏，因此想要建立模型求出污染点，只能在第一问中求得的解基础上想办法，对于数据搜寻的过程，我们独创了“环形辐射搜索法”，简便且准确，这都是第三问中的优点，不足之处就在于，第三问中的假设确实对结果影响较小，但是如果想再优化时，可以考虑把三问中的二维模型改为三维，即考虑高度的影响。对于模型四，我们只是给出了分析方法，但根据现有数据，无法进行计算，因此我们只是对于计算方法给出了详细的解释。有具体数据时可以进行计算。该方法的优点是直接深入问题本质，缺点是二阶偏微分方程无解析解，只能给出数值解，计算相对较为麻烦。而对模型再进行简化，又会带来计算误差。七、参考文献 1 尹骏.上海市城郊梯度土壤重金属空间分布特征和评价研究.上海师范大学硕士研究生论文.2010.42 成伟.数据挖掘技术支持下的土壤重金属污染评价系统的研究.浙江大学博士研究生论文.200