[理学]城市重金属污染.doc

城市表层土壤重金属污染分析摘要对于问题一，根据题中的数据，用Matlab软件绘出每个元素在该城市的空间分布及等值线图（见图2），对此图分析可知该城市表层土壤中重金属的含量分布：As含量平均值的高低顺序为工业区生活区公园绿地区交通区山区；Cd含量平均值的高低顺序为工业区交通区生活区公园绿地区山区；Cr含量平均值的高低顺序为生活区交通区工业区公园绿地区山区；Cu含量平均值的高低顺序为工业区交通区生活区公园绿地区山区；Hg含量平均值的高低顺序为生活区工业区交通区公园绿地区山区；Ni含量平均值的高低顺序为工业区生活区交通区山区公园绿地区；Pb含量平均值的高低顺序为工业区生活区公园绿地区交通区山区；Zn含量平均值的高低顺序为工业区交通区生活区公园绿地区山区。由单因子污染指数分析可知，该城市各功能区表层土壤均存在一定程度的As、Cd、Cr、Cu、Hg、Ni、Pb、Zn污染。其中，Hg的污染指数最高，其次为Cu，As的污染指数最低。由综合污染指数分析可知，该城市的生活区、工业区、交通区均属于重度污染；公园绿地区属于中度污染；山区属于轻度污染。对于问题二，以As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种重金属元素单因子污染指数均值作为因子分析，经SPSS16.0统计软件进行因子分析求其相关矩阵、特征值、累积贡献率。提取特征根较大的前三个因子，累积贡献率为97.317%（90%），满足因子分析的的原则。经过分析，该城区土壤重金属污染的主要来源为：工矿企业三废排放、交通区汽车尾气排放和居民生活中的污染物。对于问题三，基于不同区的采样点的各种重金属元素含量不同，再结合第一问中的单因子污染指数均值，确定每个区每个采样点浓度的权值。在此基础上建立数学规划模型进行求解，最终确定的生活区的污染源的坐标为（9328,4311,24），工业区污染源坐标为（2383,3692,7），山区污染源坐标为（22193,12185,79），交通区污染源坐标为（3299,6018,4），公园绿地区污染源坐标为（4153,2299,73）。对于问题四，为了弥补模型三的不足，建立对流弥散模型，运用Laplace变换求出重金属污染的滞留浓度。由菲克定律可知，重金属元素会在土壤表层不同方向进行扩散，考虑到地质元素的演变模式，在风力、水力和植物根系作用下，重金属元素会在环境三因子之间转移，基于此搜集有关信息，来确定土壤中重金属的滞留浓度、液相浓度、吸附浓度之间的关系以及植物对重金属的吸收率，从而重新确定319个采样点的重金属元素的浓度，然后依照问题三的权值确定方法和模型找出污染源。关键词：重金属污染；空间分布；因子分析法；欧氏距离；对流弥散32一、问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（010 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求你们通过数学建模来完成以下任务：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、问题分析针对问题一，首先利用题目附件1中的采样点的位置数据、附件2中的主要重金属元素在采样点处的浓度，利用MATLAB软件画出它们的空间分布图并分析；然后，以该城市自然区8种主要重金属元素的背景值为评价标准，采用单因子污染指数法和内梅罗综合污染指数评价法，计算土壤重金属单项污染指数和综合污染指数，并分析该城区内不同区域重金属的污染程度。针对问题二，通过该城市表层土壤中8种重金属的含量数据，应用统计手段及SPSS16.0统计软件，采用因子分析法对该城市表层土壤重金属采样数据进行详细分析研究，再结合问题一中各种重金属元素的空间分布，便可确定重金属污染的主要原因。针对问题三，视土壤中重金属元素浓度为稳态，根据污染源在不同区域的位置对于该区内的采样点的影响有轻有重，并且采样点位于五个区，考虑到重金属元素的浓度在每个区的差别较大，因此我们分五个区来讨论，对于每个区分别找出污染源。对不同区域分别考虑8种重金属元素浓度对采样点重金属污染影响的权重，用试探法在空间内确定一点使得该点到所有采样点的加权距离达到最大，则该点即为污染源。然而在污染源数目较多的情况下，要确定最优解的模型往往是困难的，我们从数学规划模型出发建立了确定污染源位置的模型。 针对问题四，在第三问中把土壤中的重金属元素的含量视为稳态来建模，而实际上，重金属元素在土壤中含量是非稳态的，含重金属浓度较高的污染表土容易在风力、水力和作物根系的作用下进入大气和水体中，根据菲克定律，其浓度会在各个方向进行同同程度的扩散，这里只考虑水平面上和竖直方向的扩散，根据以上分析建立一维反应溶质运移的对流弥散模型，在初始浓度为零时，半无限一维空间内第三类边界条件下利用Laplace变换进行求解。根据同种元素不同方向扩散浓度之间的关系，确定具体每个采样点的新浓度，然后对得到的新浓度按照模型三的模型进行求解，最终确定污染源的位置。三、模型假设1.假设采样点的设置全面，测量数据准确；2.假设该城市城区按照功能划分的5类区，短时间内各个城区的人口、环境没有较大的变动；3.假设重金属传播扩散只是在水平面上和竖直方向上；4.假设重金属污染主要与浓度大小有关；5.假设每个区受污染程度互不相关。四、符号说明表层土壤中重金属元素的污染指数元素的实测值（ug/g）元素的评价标准（ug/g）样点表层土壤中重金属的综合污染指数样点表层土壤中重金属元素所有单因子污染指数的最大值样点表层土壤中重金属元素所有单因子污染指数的平均值采样点的横坐标（m）采样点的纵坐标(m)采样点的海拔（m）系统的影响程度所研究范围内的凸可行域中现有受影响对象的数量第个受影响的对象与污染源之间的相对不适度的大于零的加权因子五、模型的建立与求解5.1问题一模型的建立与求解5.1.1表层土壤重金属的空间分布根据题目附件1、附件2中所给的数据，先进行数据的相关处理，作出不同区域表层土壤重金属平均值的柱状分布图（见图1），再利用MATLAB作出8种主要重金属元素在该城区的空间分布图（见图2）。图1不同区域表层土壤重金属的柱状分布图(a)表层土壤重金属As的空间分布图及等值线图(b)表层土壤重金属Cd的空间分布图及等值线图(c)表层土壤重金属Cr的空间分布图及等值线图(d)表层土壤重金属Cu的空间分布图及等值线图(e)表层土壤重金属Hg的空间分布图及等值线图(f)表层土壤重金属Ni的空间分布图及等值线图(g)表层土壤重金属Pb的空间分布图及等值线图(h)表层土壤重金属Zn的空间分布图及等值线图图2城市表层土壤重金属的空间分布图及等值线图该城市的8种主要重金属元素的空间分布如图2所示。其中，图2(a)为土壤重金属As的空间分布图及等值线图，其中以工业区含量为最高，生活区和公园绿地区次之；图2（b）为土壤重金属Cd的空间分布图及等值线图，可以看出工业区、交通区含量最高，山区含量最低；图2（c）为土壤重金属Cr的空间分布图及等值线图，可知生活区含量最高，工业区含量次之；图2（d）为土壤重金属Cu的空间分布图及等值线图，在工业区的含量最高，其均值高达127.54ug/g，山区最低，其均值为17.32ug/g；图2（e）为土壤重金属Hg的空间分布图及等值线图，可见山区含量最低，生活区含量最高，其均值是山区的23.25倍；图2（f）为土壤重金属Ni的空间分布图及等值线图，工业区、生活区含量最高，其它区域含量相差不大；图2（g）为土壤重金属Pb的空间分布图及等值线图,工业区含量最高；图2（h）为土壤重金属Zn的空间分布图及等值线图,工业区含量最大，生活区、交通区次之并两区含量相差很小。5.1.2表层土壤重金属污染程度评价（1）评价标准采用该城市自然区土壤背景值作为此次表层土壤重金属As、Cd、Cr、Cu、Hg、Ni、Pb、Zn污染程度的评价标准。表1 As、Cd、Cr、Cu、Hg、Ni、Pb、Zn污染程度的评价标准（ug/g） 土壤背景值AsCdCrCuHgNiPbZn3.60.133113.20.3512.33169（2）评价方法运用单因子指数法和内梅罗指数综合评价法对该城市表层土壤重金属As、Cd、Cr、Cu、Hg、Ni、Pb、Zn污染程度进行现状评价。a.单因子指数法b.内梅罗综合污染指数评价法内梅罗指数反映了各污染元素对土壤的作用，同时突出了高浓度污染元素对土壤环境质量的影响。(3)污染指数分级标准根据计算得出的各区域中As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的单因子污染指数和内梅罗综合污染指数，参照绿色食品产地环境质量现状评价纲要（试行）（中国绿色食品发展中心，1994），将该城区表层土壤的单因子污染指数和内梅罗综合污染指数进行了分级（表2、表3）。表2单因子指数污染分级标准等级划分1234单因子污染指数污染程度非污染轻度污染中等污染重度污染表3内梅罗指数综合污染分级标准等级划分12345综合污染指数污染程度安全警戒线轻度污染中等污染重度污染5.1.3表层土壤重金属污染程度评价结果与讨论（1）不同功能区表层土壤重金属的含量水平表4不同功能区表层土壤重金属含量（单位：ug/g；n：该区采样点）元素生活区（n=44）工业区（n=36）最大值最小值平均值标准偏差最大值最小值平均值标准偏差As11.452.346.272.1321.871.617.254.18Cd1.040.090.290.181.090.110.390.23Cr744.4618.4669.02106.66285.5815.4053.4043.39Cu248.859.7349.4046.622528.4812.70127.54409.14Hg0.550.010.930.1013.560.010.642.21Ni32.808.8118.345.6041.704.2719.818.25Pb472.4824.4369.1171.50434.8031.2493.6484.14Zn2893.474337237.01438.571626.0256.33277.93345.92元素山区（n=66）交通区（n=138）最大值最小值平均值标准偏差最大值最小值平均值标准偏差As10.991.774.041.7930.101.615.713.22Cd0.410.040.150.081.600.050.360.24Cr173.3416.2039.9624.41920.8015.3258.0581.31Cu69.062.2917.3210.651364.9012.3462.21119.79Hg0.210.010.040.0316.000.010.4421.72Ni74.035.5115.4510.35142.506.1917.6211.74Pb113.8419.6836.5617.60181.5022.0163.5332.41Zn229.8032.8673.3030.713760.8040.92242.85383.39元素公园绿地区（n=35）最大值最小值平均值标准偏差As11.682.776.262.00Cd1.030.100.280.23Cr96.2816.3143.6414.63Cu143.319.0430.1922.36Hg1.340.100.110.21Ni29.107.6015.294.90Pb227.4026.8966.7145.18Zn1389.3937.14154.24227.60生活区、工业区、山区、交通区和公园绿地区表层土壤中的重金属As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的平均含量均高于该城市自然区土壤背景值。其中，生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属As的平均含量分别是该城市土壤背景值的1.74倍、2.01倍、1.12倍、1.59倍和1.74倍；生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属Cd的平均含量分别是该城市土壤背景值的2.23倍、3.00倍、1.15倍、2.77倍和 2.15倍；生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属Cr的平均含量分别是该城市土壤背景值的2.23倍、1.72倍、1.29倍、1.87和1.41倍；生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属Cu的平均含量分别是该城市土壤背景值的3.74倍、9.66倍、1.31倍、4.71和2.29倍；生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属Hg的平均含量分别是该城市土壤背景值的26.57倍、18.29倍、1.14倍、12.57和3.14倍；生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属Ni的平均含量分别是该城市土壤背景值的1.49倍、1.61倍、1.26倍、1.43和1.24倍；生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属Pb的平均含量分别是该城市土壤背景值的2.23倍、3.02倍、1.18倍、2.05和2.15倍；生活区、工业区、山区、交通区和公园绿地区表层土壤中重金属Zn的平均含量分别是该城市土壤背景值的3.43倍、4.03倍、1.06倍、3.52和2.24倍(表4)。（2）不同功能区表层土壤重金属污染程度评价本文采用单因子污染指数法和内梅罗综合污染指数评价法对该城市的生活区、工业区、山区、交通区、公园绿地区表层土壤中重金属As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的污染程度进行了评价，结果见表5所示。表5城市不同功能区表层土壤重金属污染指数区域单因子污染指数均值内梅罗综合污染指数生活区1.742.232.333.7426.571.492.233.4319.18工业区2.013.001.729.6618.291.613.024.0313.49山区1.121.151.291.311.141.261.181.061.25交通区1.592.771.874.7112.571.432.053.529.29公园绿地区1.742.151.412.293.141.242.152.242.65图3不同功能区表层土壤重金属污染程度比较由单因子污染指数分析可知，该城市各功能区表层土壤中As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的单因子污染指数均大于1，说明各功能区表层土壤均存在一定程度的As、Cd、Cr、Cu、Hg、Ni、Pb、Zn污染。As、Cd、Cr、Cu、Hg、Ni、Pb、Zn在生活区、工业区、山区、交通区和公园绿地区的单因子污染指数分别为1.122.01、1.153.00、1.292.33、1.319.66、1.1426.57、1.241.61、1.183.02、1.064.03，Hg的污染指数最高，其次为Cu，As的污染指数最低。在该城市的五个区域中Hg污染最严重，在山区属于轻度污染，其余四个区域都属于重度污染；其次为Cu，在山区属于轻度污染，在公园绿地区属于中度污染，其余三区都属于重度污染；As的污染较轻，除了工业区为中度污染外，其余各城区均属于轻度污染。由综合污染指数分析可知，该城市的生活区、工业区、交通区表层土壤As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的综合污染指数均大于3，综合污染指数分别为19.18、13.49、9.29，均属于重度污染；公园绿地区的表层土壤重金属综合污染指数为2.65，属于中度污染；山区的表层土壤重金属综合污染指数为1.25，属于轻度污染。由上述分析和图3可见，该城市的生活区、工业区表层土壤重金属污染最为严重，其次为交通区，山区污染较轻。该城市各功能区表层土壤重金属的污染状况基本上反映了该市土地利用格局，与实际情况相符合。5.2问题二模型的建立与求解5.2.1模型的准备（1）研究方法因子分析从变量的相关矩阵出发将一个维的随机向量分解成低于个且有代表性的公因子和一个特殊的维向量，使其公因子取得最佳的个数，从而使对维的随机向量的研究转化成对较少个数的公因子的研究。设有个样本个指标构成样本空间，因子分析过程一般经过以下步骤：(1)原始数据的标准化，标准化的公式为，其中为第个样本的第个指标值，而和分别为第指标的均值和和标准差。标准化的目的在于消除不同变量的量纲的影响，而且标准化转化不会改变变量的量纲的影响，而且标准化转化不会改变变量的相关系数。(2)计算标准化数据的相关系数阵，求出相关系数矩阵的特征值和特征向量。(3)进行正交变换，使用方差最大法。其目的是使因子载荷两极分化，而且旋转后的因子仍然正交。(4)确定因子个数，计算因子得分，进行统计分析。因子分析只强调变量的离差（变量化），而不强调变量在样品中的比重（百分含量）。因子分析的数学模型中，通过正交的方差最大旋转法使每一个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小。变量或因子的重要程度都是以其方差大小来衡量的。载荷每个变量因子负荷代表着在系统中的作用或重要性程度，以各个变量目标因子载荷平方与因子方差贡献率乘积作为变量的权重，构成一个判别污染主要原因的综合指标，而且因子分析是一个客观计算同主观思维相结合的过程，其最终结果需要结合实际与已知信息进行主观评价。（2）方差最大正交旋转法的基本原理因子旋转方法有正交旋转和斜交旋转两类，对公共因子作正交旋转就是对载荷矩阵作一正交变换，右乘正交矩阵，使得有更鲜明的实际意义。旋转以后的公共因子向量为，它的各个分量为也是互不相关的公共因子。根据正交矩阵的不同选取方式，将构造出不同的正交旋转的方法。实践中，我们常用的方法是最大误差旋转法。令，则的第j列元素平方的相对误差可定义为，j=1,2,r用除以是为了消除各个原始数据变量对公共因子依赖程度不同的影响，选择除数是因为第行平方和=取是为了消除符号不同的影响。全部公共因子各自载荷之间的总方差为现在就是要寻找一个正交矩阵，经过对已知的载荷矩阵的正交变换后，新的因子载荷矩阵中的元素能使最大。5.2.2模型的建立与求解对某城市城区土壤地质环境重金属含量统计结果见表6。表6 土壤含两种金属分析表元素平均值标准差变异系数最小值最大值背景值超标率As5.683.020.531.805.403.600.58Cd0.300.220.740.070.190.131.00Cr53.5169.891.3113.0049.0031.000.73Cu55.02162.662.966.0020.4013.201.00Hg0.301.635.430.020.050.041.00Ni17.269.930.584.7019.9012.300.40Pb61.7449.980.8119.0043.0031.001.00Zn201.21338.701.6841.0097.0069.001.00 由表可以看出，Cd、Hg、Pb、Zn的含量都超出了当地金属含量的背景值。变异系数反映了样品变异程度的统计量。Hg的变异系数最大，达到5.43。对比平均值发现，Zn和Cr的含量很高，从大到小超标率为Hg Cu Zn Cd Pb Cr As Ni。以上分析表明，该城区土壤已存在一定程度的重金属积累现象。由此可知，该城市各个功能区表层土壤重金属元素含量的数据特征完全符合因子分析的要求，本题以As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种主要重金属元素的污染指数作因子分析，找出污染的主要原因。以下通过问题一模型的单因子污染指数均值（见表5），应用SPSS16.0统计软件进行因子分析，可以得出以下的结果。（1）表7给出了参与因子分析的样本观测数据（即各个区域重金属污染指数）各个变量的均值、标准差和参与计算的样品数N。表7描述统计量均值标准差NAs1.64000.327805Cd2.26000.716385Cr1.70400.374945Cu4.34203.247415Hg12.342010.583345Ni1.40600.156625Pb2.12600.654095Zn2.85601.199355（2）给出该城市表层土壤As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种重金属单因子污染指数均值的相关系数矩阵，如表8所示。表8相关系数矩阵AsCdCrCuHgNiPbZnAs1.0000.8640.4800.7660.6140.6910.9790.846Cd0.8641.0000.5330.8310.5780.7600.8890.946Cr0.4800.5331.0000.3470.9430.6840.4500.754Cu0.7660.8310.3471.0000.5530.8980.8810.825Hg0.6140.5780.9430.5531.0000.8410.6200.809Ni0.6910.7600.6840.8980.8411.0000.7890.889Pb0.9790.8890.4500.8810.6200.7891.0000.874Zn0.8460.9460.7540.8250.8090.8890.8741.000可见，Pb和As的相关性最好，相关系数最大，为0.979，以下依次是Zn和Cd，Hg和Cr的相关性较好，相关系数分别为0.946和0.943，Cr和As、Cu和Cr、Pb和Cr的相关性不是很好，相关系数分别为0.480、0.347、0.450，其它元素之间的相关性都很好。从成因上来分析，相关性较好的元素可能在成因上有一定的关联。（3）图4显示的是因子分析结果的碎石图。碎石图的纵轴为特征根，横轴为特征根序号，特征根按大小顺序进行排序。该图有一个明显的拐点，该点之前是与大因子连接的陡峭的折线，之后是与小因子连接的平缓折线。从碎石图中可以看出，前三个公因子对应的特征值较大，之后比较平缓，从而可以看出应提取的公共因子个数。图4碎石图（4）表9为各公共因子的特征根和贡献率。公因子的特征根，从大到小排列；贡献率为各公因子所解释的方差占总方差的百分比，即各公因子的特征根占总特征根的百分比；累积贡献率为各公因子方差占总方差的累计百分比，前三个公因子的累计百分比为97.317%，表明这三个因子可解释重金属元素的绝大多数信息，满足因子分析的原则。表9提取3个公因子的特征根和贡献率因子数初始情形提取后的情形因子旋转后的情形特征根贡献率累积贡献率特征根贡献率累积贡献率特征根贡献率累积贡献率16.27478.42478.4246.27478.42478.4243.40542.56442.56421.11313.91092.3351.11313.91092.3352.59732.46775.03030.3994.98297.3170.3994.98297.3171.78322.28797.31740.2152.683100.00054.36E-165.458E-15100.00068.64E-171.081E-15100.00075.38E-176.733E-16100.0008-1.53E-16-1.914E-15100.000提取方法：主成分分析从表9还可看出旋转前后总的累积贡献率没有发生变化，即总的信息量没有损失。旋转之后，主因子1的方差贡献率最大，为42.564%，主因子2和主因子3的方差贡献率都大于20%。这可以解释为因子1可能为该城市的主要污染原因，对该市的重金属污染的贡献率最大，因子2和因子3对该市重金属污染有重要作用。（5）因子分析的主要目的是将具有相近的因子载荷的各个变量置于一个公因子之下，正交方差最大旋转使每一个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小，以便对因子的意义作出更合理的解释，输出结果见表10和表11。表10 未旋转的因子载荷矩阵因子123Zn0.9830.0420.055Ni0.9270.100-0.361Pb0.925-0.3230.114Cd0.913-0.2390.144Ar0.889-0.2740.310Cu0.873-0.347-0.342Hg0.8310.532-0.035Cr0.7190.6790.132表11 旋转后的因子载荷矩阵因子123As0.9140.2790.220Pb0.8630.2390.413Cd0.8280.3050.365Zn0.6890.5690.416Cr0.2220.9680.092Hg0.2890.8880.318Cu0.6050.1450.783Ni0.4060.5460.732由表10和表11可见，旋转前后因子载荷的变量结果基本一致。变量与某一个因子的联系系数绝对值（载荷）越大，则该因子与变量关系越近。正交因子解说明：因子1为As、Pb和Cd的组合，因子2为Cr和Hg的组合，因子3为Cu和Ni的组合。5.2.3结果分析结合8种主要重金属在该城区的空间分布及对照前面的统计分析结果，可以得出以下结论： 因子1为元素As、Pb和元素Cd的组合，主要来源为印染厂的废水排放及交通源汽车尾气的排放。 因子2为元素Cr和元素Hg的组合，污染较为严重，主要来源就是厂矿企业的污水排放及汽车尾气的排放。因子3为元素Cu和元素Ni的组合，主要来源于化工行业、塑料、橡胶和印染行业的三废排放，以及城市商业活动、城市居民生活累加到土壤中的Cu，以及交通源。 综上所述，结合元素之间的相关性，可得出该城区土壤重金属污染的主要原因是：工矿企业三废排放、居民燃煤污染、交通源汽车尾气排放和居民生活污染。5.3问题三模型的建立与求解5.3.1重金属污染物的传播特征分析重金属污染物并不是静止不变的存在于土壤中，而是会在三种环境因子间进行传播。城市重金属污染的传播受自然因素和人为因素双重作用的影响。自然因素主要指风力、重力和降水等作用，人为因素主要指交通流和人流作用。传播模式见图5。蒸发重力沉降降水风力作用重力沉降降水风力作用人为作用交通流作用渗漏、降水灌溉大气土壤水体图5重金属在城市各环境因子的传播途径由于人类的生产和生活将大量的重金属污染物排入土壤中，使用农药和用污染水灌溉。各环境因子中的污染物通过自然和人为原因进行转化。灰尘中的重金属在交通流和人流作用下被扬起，也可以进入水体成为水体沉淀物，也可能成为土壤的一部分；大气中的被吸附在颗粒物表面通过重力沉降也回归到土壤中；而土壤中的重金属由于降水和渗漏进入水体。城市土壤是重金属的载体，它既是空气中可悬浮颗粒、灰尘的归宿地，也是水体污染通过灌溉、渗漏、循环而最终沉积的场所。 在建模的过程中，我们主要考虑了重金属元素的浓度和传播距离对其传播的影响。距离污染源越近，传播距离越小，污染越严重，反之污染越轻。5.3.2模型的准备确定污染源位置的问题，就是要在一定的研究区域内寻找一个点，即污染源位置，使得污染源设置在这一点位上能够保证对于所研究区域的一组受污染对象的总体污染程度达到最小，我们采用最佳定位的数学模型，建立此模型的前提：(1)已经设定了坐标系。(2)每一分布点的坐标是已知的，即分布点的位置是确定的。(3)每一分布点不是等价的，各有自己的权重因子，且这些权重值都已确定。(4)从污染源位置到各分布点的距离是可以计算的。在满足这些前提下，上述污染源位置的确定问题在数学上可以叙述如下：在三维凸边体可行区域内求一点,使得在区域内该点到给定的一组点的加权距离达到最大，即 （1）式中，位于点上的污染源与位于上的第个受影响对象的欧式距离，即。5.3.3权值的确定根据题目中附件1给出了319个分布点的位置坐标，首先把5个区的对应的所有点分别找出来，对每个区的每个分布点计算它在这个区所占的权重值。我们根据第一问单因子污染指数均值和不同元素在该点的浓度确定各个元素在各个点的比重后作归一化处理得到每个点在其所在区的权重值。对于生活区，有24个编号分别为12、13、16、18、20、21、23、25、33、36、39、42、63、67、68、72、94、106、107、152、154、156、157、158、176、180、183、184、186、187、195、198、199、243、254、258、259、267、268、269、272、273、276、306的采样点，我们对每个采样点计算浓度影响对其权值。 以编号12的采样点为例，我们定义该点为生活区的第一个点，由此对生活区的所有采样点重新编号为1到24。重新编号后，编号为1的采样点的相关数值见表12。表12 生活区编号为1的采样点的相关数值生活区12345678AsCdCrCuHgNiPbZn各种元素浓度9.170.28743.9445.770.16819.762.74223.16单因子均值1.742.232.233.4726.571.492.233.43为了计算权重，我们作如下定义由此，对剩下的23个元素分别做相同处理，求得24个值，则可得权重的计算公式为依照上述采样点的权重计算方法分别求出五个区每个采样点的权重。权重因子越大表明该点受污染源影响程度越高，反之则越低。5.3.4模型的建立与求解在满足模型所有前提的情况下，由（1）式给出的表达式是不易求解的，但是由于该模型的总体特性，通过分离两项优化操作就可将其转化成如下标准优化数学模型,则污染源最佳位置的选定问题的数学模型为 （2） （3）对每个区进行建模，可见式（2）、（3）为标准数学规划模型，应用规划论的算法，根据不同的精度要求，用试探法进行求解，程序详见附录1。由此可得到每个区污染源的相对位置（见表13），再由图6、图7可以直观地看出每个区域污染源的位置。表13每个区的污染源的相对位置坐标生活区工业区山区交通区公园绿地区X(m)932823832219332994153Y(m)431136921218560182299海拔(m)24779473图6采样点的空间位置及不同功能区的污染源位置图图7采样点的空间位置及不同功能区的污染源位置的局部放大图5.4问题四的求解5.4.1问题三模型的优缺点（1）模型的优点污染源位置的确定方法，则推广到二维坐标系上为超市、市场的最佳位置的选择提供了一种简单实用的方法，可以确定出影响某一区域人口的最佳地点。 此法也可以用来解决相反的问题，即已知污染源确定可影响的范围。权重因子的值可由污染源的源强给定。这一问题的解决可以从空间分布上确定住房选址、社交活动、交通运行在影响区以外或是以内，对进行区域的选择具有很好的很好的参考价值。（2）模型的缺点此模型确定权重时，只考虑了采样点土壤中每个点的重金属的具体浓度对其位置的影响，忽略了土壤中的水分对其浓度的影响、重金属自身的吸附降解和作物根系对重金属的吸入和排出，可能使计算结果存在一定误差。5.4.2需要另外采集的信息（1）土壤中水分的信息土壤受到重金属污染后，含重金属浓度较高的污染表土容易在风力和水力的作用下分别进入水体中，导致大气污染、地表水污染等问题。重金属元素在土壤中传播时，常发生以下反应，即不同重金属元素之间的相互作用、同一种重金属元素在环境中三因子中的转换等。如沉淀与溶解、吸附与解吸、通过植物进行交换等。（2）土壤层采样点的深度本题主要研究城市表层重金属含量，根据环境土壤学表层土壤可分为两层：耕作层和犁底层。耕作层：受耕作、施肥、灌溉影响最强烈的土壤层，厚度一般约为20cm左右。此层易受生产活动和地表生物、气候条件的影响，一般疏松多孔，干湿交替频繁、温度变化大，物质转化快，含有效养分多。根系主要分布与这一层，一般约占全部根系总量的60%。犁底层：位于耕作层之下，厚约6-8cm。典型的犁地层很紧实，孔隙度小，非毛管孔隙（大孔隙）少，毛管空隙（小空隙）多，所以通气性差，透水性不良，结构常呈片状，甚至有明显可见的水平层理。这是经常受耕畜和犁的压力以及通过降水、灌溉使黏力沉积而形成的。（3）植物体系植物体系通过生物代谢对污染物进行吸附、交换、沉淀或降解作用，污染物通过根毛细胞的作用积累与植物的茎、叶和果实部分，由此可以得出在植物比较密集的地区，植物体系对重金属的浓度影响较大。所以植物体系可作为研究地质环境演变模式的信息。5.4.3模型准备（1）菲克第二定律菲克第二定律是在第一定律的基础上推导出来的，菲克第二定律指出，在非稳态扩散过程中，在距离x处，浓度随时间的变化率等于该处的扩散通量随距离变化率的负值。所以可得如果扩散系数T与浓度无关，则上式可以写成上式中，C为扩散物质的体积浓度（kg/m）,t为扩散时间，x为距离。实际上，固溶体中溶质原子的扩散系数T是随浓度变化的，为了使求解扩散方程简单些，往往近似把T看做恒量处理。（2）水动力弥散系数的测定水动力弥散系数是表征在一定流速下，多孔介质对某种污染物质弥散能力的影响。水动力弥散系数可通过室内弥散实验测定，其实验原理如下：以人工配置的均质各向岩样，进行示踪剂注入实验。其计算公式为： （4）式中：时刻计算点的浓度；观测定的浓度峰值；计算点的坐标； 时间；弥散系数地下水的实际流速常系数，；峰值到达时间，。根据（4）式绘制理论曲线，进而可以求得水动力弥散系数。（3）土壤体积含水量和土壤中水的平均空隙流速 土壤体积含水量可通过专用仪器测得，如Sentek公司的Diviner2000便携式土壤水分监测系统，能够帮助用户实现在一定区域范围内的多点多深度土壤水分测量。土壤中水的平均空隙流速也可通过专用测量仪器测量。5.4.4模型建立在土壤中，重金属局部的传播过程，即吸附、降解等，常常与土壤中有机质的数量和类型以及土壤中微生物活性有关，而土壤中有机质的含量及微生物活性的大小又与土壤深度有关。通常有机质的浓度和微生物活性在表层土壤内最大，而后随着深度的增加而减小。通过实验数据表明，土壤中有机质的含量及微生物活性的大小与土壤深度的关系呈S型曲线。MFlury et al(1998)将溶质的降解和吸附过程与土壤深度的关系用如下函数进行了概化，并且用实验数据进行了验证。式中：、随深度变化的降解和吸附系数；、常数，代表和的最大值；深度；经验参数，其大小决定了降解和吸附过程随深度变化的快慢程度。人们在研究时通常将降解和吸附系数作为常数处理，即只考虑某一范围内的平均情况。重金属的传播非稳态，存在吸附和分解，并且有作物根系吸出和排出，由此对耕作层建立一维对流弥散模型式中：土壤体积含水量；土壤容重；土壤中水的平均孔隙流速；土壤中重金属的滞留浓度；土壤中重金属吸附浓度；由植物根系引起的重金属吸收或排除率。 对流弥散模型描述由对流和弥散引起的土壤溶质迁移。对于重金属来说，在对流和弥散过程中常伴随有较为强烈的的吸附或分解过程。我们考虑初始浓度为零时半无限一维空间的对流弥散模型运移情况，初边值条件如下式中：土壤中重金属在液相的浓度。此模型可用Laplace变换进行求解。 根据表层土壤的分类及其特点，以及考虑到重金属元素的沉积，同时采样点的浓度会向水平方向和竖直方向扩散，直至滞留浓度最小为止。对于重金属耕作层和犁地层浓度，后者中重金属浓度较大。又由于犁底层根系很少，所以对于犁底层可不考虑作物根系的作用，故此层重金属浓度主要是滞留浓度，而耕作层主要是可在环境三因子中相互转化的重金属浓度。 把每个采样点的浓度分为三部分来考虑，即滞留浓度、液相浓度和作物根系对浓度的影响来考虑，其中作物根系浓度可有实验测得，根据这三个量与深度和高度的关系，确定采样点的新的浓度大小，之后再与单因子均值作乘积，归一化处理后得每个采样点浓度的权值，然后利用问题三的模型对欧式距离加权求得最大值，从而确定每个区新的污染源。六、参考文献1 宋来忠，王志明，数学建模与实验M，北京：科学出版社，20052 杨桂元，李天胜，徐军，数学建模应用实例M，合肥：合肥工业大学出版社，20073 佟春生，畅建霞，王义民，系统工程的理论与方法概述M，北京：国防工业出版社，20054 王雄军，赖健清，鲁艳红，基于因子分析法研究太原市土壤重金属污染的主要来源J，生态环境学报，17（2）：671-676，20085 马旺海，曹斌，杨进峰，城市重金属污染特征J，中央民族大学学报，17（1）：66-73，20086 谢小进，康建成，李卫江，上海市不同工业类型土壤重金属含量差异性分析J，生态环境学报，18（5）：1777-1782，20097 张德生，沈冰，土壤中反应溶质运移的对流弥散模型及其解析解J，西安理工大学学报，17（2）：122-126，20018 隋红建，吴璇，崔岩山，土壤重金属迁移模拟研究的现状与展望J，农业工程学报，22（6）：197-200，20069 王力宾，顾光同，多元统计分析：模型、案例及SPSS应用M