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

城市表层土壤重金属污染分析摘要本文通过处理已给的青岛市各采样点的数据，对城市表层土壤重金属含量进行了有效分析。通过Kriging克里金插值拟合建立模型，得出八种重金属元素的空间分布特征、人类活动及工业发展给该地区增加的重金属量。通过各种重金属元素在不同功能区的高污染占比和相关性分析得出的6类相关元素来确定污染物。接着通过已确定的污染物的特征猜想污染源的位置，并通过拟合验证，然后作图得出具体污染源坐标。最后考虑五种不同的污染形式，并查阅资料确定它们对重金属污染的不同影响，以此来预测十年之后污染状况。对于问题（1），利用MATLAB的v4样条插值分析海拔数据确定出青岛市的海平面，确定市区的有效面积并进行区域划分。通过MATLAB的克里金工具箱对各重金属元素进行Kriging插值，得出由于人类及工业发展给该地区增加的重金属量和比率分布、8种重金属元素平均增加量比率在整个青岛的空间分布。并将污染物增加量比率与污染程度挂钩。我们可以发现污染较为集中的是西南部的交通区和工业区，污染最严重的西南交通区平均每种元素超标近160倍。对于问题（2），我们为确定污染物，分析各元素高污染在功能区占比。随后为进一步确定相关性，通过相关性分析，划分六类相关元素：As；Cu、Hg； Cu、Pb；Cr、Cu、Ni；Cr、Ni、Zn；Cd、Pb、Zn。说明造成重金属污染的污染物及各功能区的综合污染情况为：交通区工业区生活区公园绿地区山区。对于问题（3），我们就问题（2）中的污染物研究它们的传播方式和所在分布区域特征。并猜想污染状况最严重的点即为我们所要求的污染源。以As（砷）为例分析其浓度分布，并结合其他特征峰值点验证出猜想正确。元素分类进一步确定污染源坐标，对于问题（4），通过地质积累指数-内梅罗指数法、单因子综合指数-内梅罗指数法两种方式分别计算地区受污染情况，与问题（1）中数据进行比较验证。然后通过间接比较的方法结合青岛市的情况，预测十年后的污染状况，除了山区以外，剩下四个功能区的重度污染比率均超过30%，其中工业区和交通区的重度污染比率更是超过50%，环境形势变得十分严峻。关键词Kriging插值 相关性分析 回归分析 地质累积指数 内梅罗综合指数法1. 问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（010 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。题中附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现需要通过数学建模完成以下任务：1、给出8种主要重金属元素在该城区的空间分布特征，选择合适的插值方法，估算由于人类及工业发展给该地区增加的重金属量。2、通过题目提供数据分析，建模说明造成重金属污染的污染物3、根据污染物的传播特征和分布区域特征，建模确定污染源。4、进一步补充数据资料，验证你们在问题1中得到结果的合理性，由此建立模型，预测如果不加以控制十年后该地区的污染状况。2. 问题分析问题（一）：此题需要利用所给采样点海拔和所在功能区的数据可以构建起一个三维模型，运用v4差值求解得到城市区域网格分布。接下来题目的关键，是怎样将各采样点坐标与该位置重金属浓度相联系。我们在划分城市功能区域的基础上，利用各重金属浓度数据，并采用克里金工具箱的Kriging克里金插值拟合建立模型，我们建立起了这种联系。对采样点重金属含量进行插值拟合得到8种主要重金属元素的空间分布特征，之后通过结合各重金属浓度的标准值分析得出增加量。问题（二）：首先初步的将8种元素单独分析，利用每种高污染占比数值与所在功能区的对应关系推出污染物；再对各重金属元素进行相关性分析，将其分为几类组合从而确定出污染物，结合采样点所在功能区分析出污染源。问题（三）：题目的关键是要确定污染物的传播方式和分布区域特征对污染源的影响。我们以As为例来探讨了污染物的传播特征和区域分布特征，并由此判断污染源就是污染物增加量比率的峰值处，即污染程度最高的点。问题（四）：我们通过问题（二）中的结论，采用间接预测的方式，为每种污染方式找到相应的污染功能区对应，通过变化的参考量变化的比率来估算污染程度的变化，并由此预测十年后的重金属污染情况。3. 模型假设1、题目所给的数据合理正确，每个网格内的采样点具有随机代表性。2、测量的数据为同时测得，并在一定的时间范围内有效。3、不考虑地下水、河流及海洋对于污染物的影响。4、划分的网格区域较为稳定，在短时间内不会发生变化。5、城区周围数值连续变化，不存在突变情况。4. 变量与符号说明As： 砷元素的化学符号；Cd： 镉元素的化学符号；Cr： 铬元素的化学符号；Cu： 铜元素的化学符号；Hg： 汞元素的化学符号；Ni： 镍元素的化学符号；Pb： 铅元素的化学符号；Zn： 锌元素的化学符号；ij(i=1,2,3,8 j=1,2,3,319) 第i种元素在第j个采样点的重金属元素的增加量比率cij(i=1,2,3,8 j=1,2,3,319) 第i种元素在第j个采样点的重金属元素的现含量ci,sta(i=1,2,3,8) 各重金属元素含量的标准值ave 八种重金属元素平均增加量比率pi pi(i=1,2,3,319) 两种不同重金属元素在第i个采集点p p(i=1,2,3,319) 两种不同重金属元素在所有采样点浓度的平均值p 两种重金属元素的相关系数I 地质积累指数P 单因子污染指数5. 模型的建立与求解5.1. 问题（1）模型的建立与求解5.1.1. 城市区域分布我们选定横坐标（030000）、纵坐标（020000）为青岛市的城市区域。我们横坐标每隔300m取一点，纵坐标每隔200m取一点，共取到101*101的网格。由于青岛市是沿海城市，我们选定的区域也有着丰富的海区。我们利用MATLAB的v4插值求解，通过319个采样点的海拔，将海拔低于0m的区域定义为海洋，确定了青岛市的三维图像，与地图相吻合。同时对于采样点的功能区，我们认为插值点功能区与直线最近采样点功能区相同。为此我们划分了城市的功能区，如下图所示图1 城市网格区域分布图图2 城市网格区域分布图（平面，图中白点为采样点）图中可以看到山区主要位于城市海拔比较高的区域，工业区由于航运便捷、季风影响主要分布在西海岸沿海，交通区分布较为广泛，生活区、公园绿地区分布较为零散，主要夹杂在交通区内。这些都符合正常的城市规划。5.1.2. 各重金属元素的空间分布我们在划分区域的基础上，采用克里金工具箱的Kriging克里金插值法对采样点重金属含量进行插值拟合。Kriging插值法广泛地应用于地下水模拟、土壤制图等领域，是一种很有用的地质统计格网化方法。它首先考虑的是空间属性在空间位置上的变异分布，确定对一个待插点值有影响的距离范围，然后用此范围内的采样点来估计待插点的属性值。该方法在数学上可对所研究的对象提供一种最佳线性无偏估计(某点处的确定值)的方法。它是考虑了信息样品的形状、大小及与待估计块段相互间的空间位置等几何特征以及品位的空间结构之后，为达到线性、无偏和最小估计方差的估计，而对每一个样品赋与一定的系数，最后进行加权平均来估计块段品位的方法。但它仍是一种光滑的内插方法，在数据点多时，其内插的结果可信度较高。图3至图10为各元素的分布图。图3 As（砷）浓度分布图砷的分布比较分散。集中分布在东部山区，广泛分布在西南部交通区。有一处最高峰（小范围高浓度的分布），即污染最严重的区域，在东部山区；另外有三处污染较为严重，一处是中部交通区、一处在西南部沿海地带的交通区、最后一出事西部沿海的工业区；另外，西南部大面积都有较高值的砷浓度分布，说明整个西南部交通区均有较严重的砷污染。图4 Cd（镉）浓度分布图镉主要分布在西部沿海的交通区，部分集中分布在东部山区和西部沿海的生活区。污染最严重的是西南部沿海的交通区；其次是西南部沿海的生活区、东部山区以及西部沿海工业区；中南部的交通区均有大面积较为平缓的峰，说明在中南部的交通区均有大范围的轻度污染。图5 Cr（铬）浓度分布图铬的分布很集中。铬主要分布在西南部交通区和工业区、东部偏北的山区，部分分布于西南部沿海的生活区铬的浓度分布有两个明显的峰，分别是西偏南沿海的交通区和生活区，其次是东部山区；另外西部沿海及西南部沿海的工业区也有较严重的污染；西南部的交通区都有较大面积的轻微铬金属污染。图6 Cu（铜）浓度分布图铜的分布很集中。铜主要集中分布在小范围分布在西南角生活区和工业区，部分分布在西南部交通区。 有两个集中的峰值，分别是西南部沿海的生活区，工业区和交通区；其次西南部沿海的交通区受到轻微的污染。铜大部分分布在西南角交通区，部分分布在西南部生活区、西南部和中部工业区、中部的山区图7 Hg（汞）浓度分布图汞的分布很集中。主要分布在东南部的交通区，其次分布在中部工业区和生活区和南部的交通区。汞的浓度分布有五个主要的峰值，最高峰是西南部的交通区，其次是西南部沿海的工业区和中部地区小范围的工业区和山区；另外在南部沿海的交通区和西部沿海生活区都有一个较低的峰，说明有较为严重的汞污染；西南部和南部的交通区生活区都广泛受轻度汞金属污染；其他地区几乎不受污染图8 Ni（镍）浓度分布图镍主要小范围集中分布在沿海交通区、工业区和公园绿地区，另外部分零散分布在东部和中部山区，广泛分布在中部及西部交通区。汞浓度分布有两个主要峰值，分别是东部山区、西偏南部沿海的交通区和生活区；另外在东部和中部山区西南部公园区均有小范围的集中分布，在西部工业区有较为严重的汞污染，其他大部分交通区均有轻微汞污染。图9 Pb（铅）浓度分布图铅集中分布在西部及西南部沿海交通区和生活区以及部分工业区。有两个主要峰值，分别为西部沿海交通区和西部生活区，另外西部工业区也有较高的铅浓度值，有较严重的铅污染；西南部交通区大范围的受较严重的铅污染，中部交通区和东部和东南部山区均受轻微铅污染。图10 Zn（锌）浓度分布图锌的分布范围最广。锌集中分布在西南部及中部交通区和生活区，少部分分布在公园区和夹杂于交通区内的山区和工业区。没有明显的峰值，但是有西南部以及中部有大面积及严重的污染，包括西部和西南部的交通区、西南部额中部的生活区、中部的工业区和部分山区、西部沿海的公园区。另外西部山区、东北部沿海公园区有部分轻微污染。5.1.3. 各重金属元素增加量由于人类及工业的发展，城市的重金属元素含量在不断增加，已得到每种元素的平均增加量。结论如下表：表一：八种元素增加量元素AsCdCrCuHgNiPbZn平均增加量2.08172.4022.5141.82264.714.9630.74132.20从表格中可以看出，Cd，Hg，Zn，Cu污染程度高。因为单纯评判增加量没有意义，我们定义增加量比率为ij：ij=cij-ci.staci,sta100%表示第i种重金属元素在第j个取样点的增加量比率。其中cij为第i种重金属元素在第j个取样点 ci.sta为第i种重金属元素的标准值。（i=1,2,3,8 j=1,2,3,319）我们再定义八种元素平均增加量比率为ave，j:ave，j=18i=18ij由此得到每个网格的元素平均增加量比率（负数抹零），结果如图11所示：图11 元素平均增加比率上图体现的是8种重金属元素平均增加量比率在整个青岛的空间分布，我们可以发现污染较为集中的西南部的交通区和工业区。污染最严重的西南交通区平均每种元素超标近160倍。我们知道，重金属增加量越大，污染状况也会越严重。所以，对增加量的分析可化为对各功能区的污染状况分析和各种金属元素在不同区域高污染占比的比较分析。污染倍数与污染等级存在如下表的关系：表二：土壤污染分级标准平均增加量比率污染等级0无污染01轻度污染12中度污染23重度污染 将插值后各采样点的污染等级按照城市区域划分归类得到结果如下表所示：表三：区域平均污染等级污染等级海洋生活区工业区山区交通区公园区无污染28942761671442920182轻度污染015610936439571中度污染015811273246273中重度污染068559522547重度污染021016085641102总计289486860327182643475 上表就无污染比例： 工业区生活区交通区公园区山区。工业区的无污染比例最低，所受污染范围最大。就重度污染占统计总数的比率来看，污染程度为：山区公园区生活区工业区交通区。所以，我们可以得出结论：交通区和工业区的污染程度较重，山区和公园区污染程度较轻。但是可能会出现个别元素污染严重影响平均数值的情况，我们将各个元素单独分析其高污染比率，得出下表（单位为%）表四：八种元素高污染比率占比重污染比率AsCdCrCuHgNiPbZn平均海洋000000000生活区1.0411.064.4922.0042.483.007.3717.8613.66工业区4.4815.925.3123.8847.432.6512.4417.2516.17山区2.544.492.4719.0634.952.560.101.998.55交通区2.5016.085.0325.7738.781.559.6126.6417.00公园区2.0313.122.1120.0041.471.265.8917.8912.97加权平均1.767.922.7515.9727.691.614.1510.809.08从各重金属元素平均含量的地区分布来分析。竖向看每一种重金属元素在各功能区的高污染占比，As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的最高含量分别出现在工业区，交通区，工业区，交通区，工业区，工业区，工业区，交通区，从而可以发现工业区和交通区占了重金属污染的主导。再横向看各个城市功能区各种重金属元素的污染比率占比，工业区和交通区在多项元素上比率占比都很高。这个分析结果这与该城区各重金属元素单点最高值分布基本相符。从最右侧的“平均”一栏可以看出，污染程度为：山区公园区生活区工业区交通区。与平均污染等级结果基本一致。5.2. 问题（2）模型的建立与求解5.2.1. 通过重金属高污染占比来初步分析得出污染物我们再次分析表四。首先，我们纵向分析。As：砷元素在五个区中高污染占比最大的为工业区（4.48%），其次为山区（2.54%）和交通区（2.50%）。Cd：镉元素在五个区中高污染占比最大的为交通区（16.08%），其次为工业区（15.92%）。Cr：铬元素在五个区中高污染占比最大的为工业区（5.31%），其次为交通区（5.03%）。Cu：铜元素在五个区中高污染占比最大的为交通区（25.77%），其次为工业区（23.88%）。Hg：汞元素在五个区中高污染占比最大的为工业区（47.43%），其次为生活区（42.48%）和公园区（41.47%）。Ni：镍元素在五个区中高污染占比最大的为生活区（3.00%），其次为工业区（2.65%）和山区（2.56%）。Pb：铅元素在五个区中高污染占比最大的为工业区（12.44%），其次为交通区（9.61%）。Zn：锌元素在五个区中高污染占比最大的为交通区（26.64%），其次为公园区（17.89%）和生活区（17.86%）。通过查阅资料，可知：砷的主要污染来源是含砷金属采矿，冶金，化学制药，玻璃工业中的脱色剂，各种杀虫剂，硝酸盐药物，化肥，皮革，农药等等。镉的主要污染来源是电镀，采矿，有色金属冶炼，燃料，电池和化学工业等排放的废渣，污水。铬的主要污染来源是劣质化妆品原料，皮革制剂，金属部件镀铬部分等。铜的主要污染来源是铜锌矿的开采和冶炼、金属加工、机械制造、钢铁生产等。含铜废水灌溉农田，农药的施用。汞的主要污染来源是原仪表厂，食盐电解，贵金属冶炼，化妆品，照明用灯，齿科材料，燃煤，水生生物等。镍的主要污染来源是岩石风化，河流输送,开采，矿物燃烧，工业生产。铅的主要污染来源是冶炼金属产生的废水废渣，汽车排出的含铅废气。锌的主要污染来源是矿开采、冶炼加工、机械制造以及镀锌、仪器仪表、有机合成和造纸等工业的排放。对比表四的结论，我们可以归纳出青岛市8种重金属污染得污染物，结论如下：砷：主要原因是工业区的冶金、化学制药等工业三废的排放，其次是山区的砷金属采矿，以及交通区的道路表面沥青磨损导致砷进入到环境中；镉：主要原因是交通区冬季防滑剂的使用；次要原因是因为工业区电镀，采矿，冶炼，燃料，电池和化学工业等排放的污水。铬：主要原因是工业区劣质化妆品原料，皮革制剂，金属部件镀铬部分，工业颜料以及鞣革，橡胶和陶瓷原料等；次要原因是公路区汽车燃油尾气排放。铜：主要原因是公路区车辆制动衬面磨损；次要原因是工业区铜锌矿的开采和冶炼、金属加工、机械制造、钢铁生产。汞：主要原因是工业区煤和石油燃烧的污染，废水污染；次要原因是次要原因是公园绿地区的有机肥、化肥和农药的大量使用、生活区的生活垃圾如废旧电器，电路板，光管，电池等的不合理堆放以及生活中化工产品在建筑施工和其他服务业的大量使用；部分原因是交通区燃油尾气排放。镍：主要原因是生活区的生活垃圾的不合理堆放和各种含重金属的化工产品在生活中的使用；次要原因是山区岩石风化和工业区的镍矿的冶炼和镀镍、机器制造和金属加工等产业。铅：主要原因是工业区冶炼，五金，机械，电镀，化妆品，染发剂，釉彩碗碟，燃煤，次要原因是公路区燃油尾气的排放。锌：主要原因是公路区的燃油尾气排放和冬季防冻剂的使用，次要原因是工业区的矿开采、冶炼加工、机械制造以及镀锌、仪器仪表、有机会合成和造纸等工业的排放和公园绿地区的有机肥，化肥和农药的大量使用和生活区含重金属的生活垃圾不合理堆放。5.2.2. 重金属相关性分析进一步确定污染物计算由于污染物可能会产生多种重金属元素，考虑到一些重金属元素之间可能存在一定的联系，因此我们对各重金属元素进行了聚类分析，拟对造成污染的重金属元素作相关性分析。聚类分析方法是基于数值分类的思想建立起来的，它是研究样品或指标分类的一种多元统计分析方法。相关性考虑的是重金属元素两两之间的关系。如果不同的元素之间存在很高的相关性，我们有理由相信它们可能是同一污染物在污染源产生的。相关性的系数定义如下：设pi、pi分别为两种不同金属在第i个采样点的浓度，p、p分别为两种金属所有采样点浓度的平均值，则相关性系数为p=(pi-p)(pi-p)(pi-p)2(pi-p)2我们计算了两两不同种类的元素之间的相关性，结果如下表所示：表五：八种元素两两相关性元素AsCdCrCuHgNiPbZnAs0.25470.1890.15970.06440.31660.28990.2469Cd0.25470.35240.39670.26470.32940.66030.4312Cr0.1890.35240.53160.10320.71580.38280.4243Cu0.15970.39670.53160.41670.49460.520.3873Hg0.06440.26470.10320.41670.10290.29810.1958Ni0.31660.32940.71580.49460.10290.30680.4364Pb0.28990.66030.38280.520.29810.30680.4937Zn0.24690.43120.42430.38730.19580.43640.4937我们认为相关性大于0.3视为有弱相关性，大于0.4视为有较强相关性。研究发现，上三角矩阵中，As元素仅与Ni有弱相关性；Cd与Cr、Cu、Ni有弱相关性，与Pb、Zn有较强的相关性；Cr与Pb有弱相关性，与Cu、Ni、Zn有较强相关性；Cu与Zn有弱相关性，与Hg、Ni、Pb有较强相关性；Hg与其它元素没有相关性；Ni与Pb有弱相关性，与Zn有较强相关性；Pb与Zn有较强相关性。5.2.3. 重金属污染的污染物分析将有较强相关性的元素之间连线，我们可以得出以下六类元素组合：根据相关性分析，并比照5.2.1的结论，总结出以下结论：第一类，As；主要是由于采矿，冶金，化学制药以及农业生产中的化肥、农药。第二类，Cu、Hg；主要是由于金属冶炼加工、机械制造、钢铁生产和公路汽车燃油尾气排放。第三类，Cu、Pb；主要是由于矿石开采、金属冶炼加工、机械制造。第四类，Cr、Cu、Ni；主要是由于工业区金属加工，公路汽车燃油尾气排放。第五类，Cr、Ni、Zn；主要是由于公路上汽车尾气排放，道路表面磨损。第六类，Cd、Pb、Zn；主要是由于工业区的矿开采、冶炼加工、机械制造以及农业生产中有机肥、化肥和农药的大量使用。5.2.4. 模型结论综合上述两种方法的出的结论，我们可以得出问题（2）的结论，即各功能区污染程度和污染物的总结。结论如下：显然交通区的污染程度最大，分析其来源为工业生产三废的排放、汽车尾气排放生的大量含重金属的有害气体和粉尘经雨水、重力等作用富集于土壤中。需要说明的是工业废气和汽车尾气中主要含Cr、Ni、Cu、Zn、Cd，这些元素在交通区的高污染占比都排在前几位。工业区的污染程度也很大，在某些重金属元素上甚至超过了交通区的污染程度，并且其总污染指数于交通区的总污染指数相差不大。分析其污染的主要来源有重金属矿冶炼以及相关的制造产业、矿渣堆放等，并且其扩散和传播污染的途径也是比较广泛的。生活区的污染程度居中，其重金属的来源主要有生活垃圾、生活污水、废旧电子产品和废旧电池等。Cd、Cu、Hg、Pb四种元素在生活区的高污染占比很高，这些元素正是生活垃圾、生活污水和废旧电子产品中含量较高的。公园绿地区的污染程度较轻，其重金属的主要来源有农业农药、化肥和经处理后的工业废水、生活废水的灌溉等。另外，公园绿地区的植物对金属元素有一定的吸收作用，可以降低土壤中重金属元素的含量。山区的污染程度最低，其重金属污染来源于重金属矿石的开采和岩石风化。山区的污染程度最低主要是因为山区不是人类活动的集中区，尤其东北部大片山地又远离污染严重的交通区和工业区。另外与其海拔也有一定的关系，山区中植物和一些微生物也有净化的作用。5.3. 问题（3）模型的建立与求解5.3.1. 污染物传播特征与分布区域特征分析我们认为污染物是以污染源为中心，以不同的传播方式向四周扩散，而造成不同污染物浓度分布。由问题（2）得出的污染物，我们总结出主要有以下几种传播方式：1. 大气传播的重金属主要来源于工业生产、汽车尾气排放及汽车轮胎磨损产生的大量含重金属的有害气体和粉尘等，它们主要分布在工况的周围和公路的两侧，大气中的大多重金属是经自然沉降和雨淋沉降进入土壤的。2. 经过自然沉降与雨淋进入土壤的重金属污染，主要以工矿烟囱，废物堆和公路为中心，向四周及两侧扩散：由城市-郊区-农区，随距城市的距离加大而降低污染浓度。3. 在公路，铁路区域重金属污染源于含铅汽油的燃烧，汽车轮胎磨损产生的粉尘等。重金属污染元素主要为Pb，Zn，Cr，Cd。污染成条状分布，以其为轴向两侧重金属污染强度逐渐减弱。并且随时间的推移，公路，铁路上的重金属污染有很强的叠加性。4. 河流的重金属污染往往受污染源矿山控制，也可能因为地表，地下的工业废水汇集而加重。河流在单一污染源的上游到下游，由于金属元素迁移能力减弱和水体自身净化能力，化学污染强度逐渐降低。从城区重金属的分布区域特征看出，它们都是从浓度峰值出发以不规则图形向四周扩散，我们不妨猜想该中心就是污染物浓度峰值。分析浓度分布强度，我们以As元素为例，研究污染物与距离之间的关系。我们选取As在中部的一个高浓度峰值点，以该峰值点为圆心，考察在峰值点半径为1.5km内的格点污染浓度变化，用cftool拟合作出图12所示。图12 As污染源与半径关系（1）图12采用y=aebx拟合，拟合的函数为y=48.18*e-0.001594x其中r2=0.7063，相关程度不是很大。我们将图中的点按照所属类别分开，得到图13图13 As污染源与半径关系（2）该污染源周围仅存在生活区、工业区和山区，图中可以很明显地看出山区的污染情况明显好于工业区。由于生活区距离较远且含量极低，可不作拟合。分别将工业区和山区作拟合。结果如下y=73.63*e-0.001903x（工业区，r2=0.8531）y=45.3*e-0.001587x（山区，r2=0.7555）将格点如此归类后相关性有了明显的增加，且从拟合函数看来工业区明显浮于山区之上，由此可见工业区相对于山区污染物浓度更高，污染物浓度随距离增加而减少的更为缓慢。从传播特征分析，工业区内重金属污染物随大气传播时，污染物没有被净化，浓度依然保持较高水平。而山区的重金属污染随河流自上而下的流淌而净化，因此浓度减少更为迅速，并且易于呈现出椭圆形的分布。从图3的污染物分布区域特征也可以看出，污染物浓度的分布特征为不规则图形，大部分呈现为椭圆形或者短带状，这也是因为污染物自圆心心向周围扩散的传播方式不同而造成的。我们再选取几个特征明显的污染源进行更细致的拟合分析，在不同功能区块内，就传播特征来看，污染物从源头扩散向四周，依据传播方式不同，随距离增加而减少的程度不同，导致污染物分布区域图形成不规则形状。综合以上分析，将污染物的高峰值点作为污染物的源头的是十分合理的。因此我们可以判断污染状况最严重的点即为污染源。5.3.2. 污染源的位置坐标求解按照问题（2）中得出的六类元素组合分别分析其可能的污染源坐标，计算方法如下：第一步，计算出每组各元素在城区网格内的元素平均增加量比率。第二步，作出平面图直接观察有几个峰值及峰值位于的大区域。第三步，在大区域内直接查找出峰值及峰值位于的坐标。5.3.2.1. As组合As分布如图14所示：图14 As组合分布图As组合分布较为零散，但是大致有四个污染源（从强到弱，下同）：西北山区、西南交通区、西海岸沿海工业区、中部交通区，通过对MATLAB图像观察得出四个污染源坐标分别为（27900，12400）、（1500，200）、（6300，8400）、（18900，10200）。5.3.2.2. Cu、Hg组合Cu、Hg组合分布如图15所示：图15 Cu、Hg组合分布图Cu、Hg组合大致有五个污染源：西南交通区、西南工业区、西海岸沿海工业区、中部交通区、东岸交通区，通过对MATLAB图像观察得出五个污染源坐标分别为（2100，200）、（2400，3400）、（3600，5400）、（15300，9400）、（13800，2400）。5.3.2.3. Cu、Pb组合Cu、Pb组合分布如图16所示：图16 Cu、Pb组合分布图Cu、Pb组合大致有两个污染源：西南工业区、西海岸沿海工业区，通过对MATLAB图像观察得出两个污染源坐标分别为（2100，3600）、（3600，5600）。5.3.2.4. Cr、Cu、Ni组合Cr、Cu、Ni组合分布如图17所示：图17 Cr、Cu、Ni组合分布图Cu、Pb组合大致有两个污染源：西海岸沿海工业区、西南工业区，通过对MATLAB图像观察得出两个污染源坐标分别为（3600，5600）、（2100，3400）。5.3.2.5. Cr、Ni、Zn组合Cr、Ni、Zn组合分布如图18所示：图18 Cr、Ni、Zn组合分布图Cr、Ni、Zn组合分布较为零散，大致有八个污染源（从强到弱）：西部交通区、中南部交通生活区、西海岸沿海工业区、中部交通区（两处）、中部交通工业区（三处），通过对MATLAB图像观察得出八个污染源坐标分别为（3600，5800）、（9600，4800）、（6000，8400）、（7500，7200）、（10500，7400）、（13200，11400）、（14100，9600）、（14400，8200）。5.3.2.6. Cd、Pb、Zn组合Cd、Pb、Zn组合分布如图19所示图19 Cd、Pb、Zn组合分布图Cd、Pb、Zn组合大致有5个污染源：西海岸沿海工业区、西部交通区、中部交通区（两处）、中南部交通生活区，通过对MATLAB图像观察得出五个污染源坐标分别为（6300，8400）、（4800，5600）、（7500，7000）、（8700，8400）、（9600，4600）。5.3.3. 污染源分析通过计算及比对，我们大致知道了各个污染元素的分类及对应坐标，我们在地图上标注出各个污染源，与划分区域分布图对比，得出下图：图20 污染源分布对比图图20体现了污染源分布的区域特征。可以发现重金属元素污染源有许多位置重叠。并且对照图11发现，组合分析出来的污染源重叠部分与8种元素平均污染状况最严重的区域几乎重合，这更加说明我们的猜想的结论较合理。5.4. 问题（4）模型的建立与求解5.4.1. 验证问题（1）中结果为了验证问题（1）中的增加量比率与污染程度的关系，通过查阅相关资料，我们确定了采用地质累积指数-内梅罗综合指数法和单因子综合污染指数-内梅罗综合指数法两种方法来验证增加比率与污染程度关系。地质累积指数-内梅罗综合指数法侧重对于采样点综合情况的把握，单因子综合污染指数-内梅罗综合指数则侧重与对各个元素综合情况的把握，各有特点。5.4.1.1. 地质累积指数-内梅罗综合指数地质累积指数，又称地质累积指数，是Muller于1949年提出的针对沉积物种的重金属污染程度的评价方法，它不仅考虑了自然地质过程造成的背景值的影响，而且也充分注意了人为活动对重金属污染的影响，因此，该指数不仅反映了重金属分布的自然变化特征，而且可以判别人为活动对环境的影响，是区分人为活动影响的重要参数。地质累积指数定义为Im,n=log2um,nkcn,sta其中k为修正数值，通常取1.5。Im,n为在第m个采样点第n种元素的地质累积指数。um,n为在第m个采样点的第n种元素的浓度。为了评价各功能区重金属的污染程度，需要综合考虑各种元素对环境的影响，即需将各单个指标进行合成。常用的环境质量综合评价指标是内梅罗指数法，这种方法特别考虑了污染最严重的因子。记采样点的所有金属元素地质累积指数的最大值和平均值分别为Imax和Iave，则此处的内梅罗指数的定义为：Im=Imax2+Iave22地质积累指数污染划分表如下：表六：地质积累指数污染等级表地质积累指数分级污染程度0I11无污染，轻度污染1I22中度污染2I33中重度污染345极重度污染利用地质累积指数-内梅罗综合指数评价系统污染程度的具体步骤为：第一步，计算各个采样点每种元素的地质积累指数。第二步，利用内梅罗指数法计算每个采样点的综合地质累计指数。第三步，用kriging插值法对指数插值。第四步，统计每个功能区的地质积累指数，对照地质积累指数污染分级表，划分区域和分级。第五步，对照表三，作出形式类似的区域污染等级表，并作对比。此种方法强调了对与功能分区的重视，而内梅罗综合指数法更是侧重了某种元素含量的单一格点。由上表我们通过MATLAB编程计算了各个功能区的污染情况，如表七所示：表七：区域平均污染分级污染分级海洋生活区工业区山区交通区公园区12894447262178714062772027615886057810030597153303464038358135135048771022139总计289486860327182643475依照高污染所占比例，可以得出污染程度：山区公园区生活区交通区工业区，虽然结果中交通区与工业区与此前结果相反，但是此两区相差不大，污染指数均明显高于其它功能区，仍然在可接受范围内。此外，从表五与表二的对比不难看出，表五高污染（分级4、5）的比例明显比表二的小，原因可能是评价标准不同，即所谓的“重污染”的标准不同所致。5.4.1.2. 单因子综合污染指数-内梅罗综合指数评价方法采用指数法。定义单因子综合污染指数为：Pm,n=um,ncn,sta其中Pm,n表示第m个采样点第n种元素的污染指数。um,n为在第m个采样点的第n种元素的浓度。为了侧重极端污染情况，我们仍然需要引入内梅罗指数，记某一金属元素在所有采样点污染指数的最大值和平均值分别为Pmax和Pave，此处的内梅罗指数的定义为：Pn=Pmax2+Pave22内梅罗指数的污染情况等级表如下：表八：内梅罗指数污染等级表内梅罗指数分级污染程度P0.70清洁0.7P11尚清洁1P22轻度污染2P33中度污染P34重度污染利用单因子综合污染指数-内梅罗综合指数评价系统污染程度的具体步骤为：第一步，计算各个采样点每种元素的单因子污染指数。第二步，利用内梅罗指数法计算每个功能区各单项元素的污染指数。第三步，对照表四，划分内梅罗综合指数分等级。此方法重在考察对于单一元素对于整体污染情况的影响，对于环境变量显得更加敏感。结果如下表所示：表九：八种元素内梅罗指数重污染比率AsCdCrCuHgNiPbZn海洋00000000生活区1.473.489.126.506.691.126.8918.49工业区2.593.733.5464.86161.981.386.4210.51山区1.311.402.171.822.532.261.731.54交通区3.445.3311.2535.00191.834.282.7724.01公园区1.493.421.293.7516.110.983.408.91从表中可以明显看出交通区和工业区在评价指标中几乎所有的统计项都超出了重度污染的标度，这可能是因为几乎所有的元素的极大值都处在工业区或者生活区，大大影响了此二区的污染程度。5.4.2. 十年后城市污染情况预测由于重金属元素分布、扩散的以及城市发展的复杂性，十年后的污染情况变得极难预测，我们只能粗略地通过间接方式计算出总体污染情况。下表为根据问题（2）中的分析对于各个功能区的污染程度和对应的参照数值来源。表十：污染区主要污染原因污染形式主要影响范围预测参考年比值（%）尾气排放交通区机动车保有量114.41工业污染工业区工业产值（第三产业）111.80水源污染公园绿地区、生活区、山区污染水源排放量102.38生活污染生活区人口101.52农业污染公园绿地区化肥使用量102.64通过查阅大量的资料，设预测参考值在当年均为1，尾气污染、工业污染、水源污染、生活污染、农业污染十年后的参考值分别为3.8427、3.0508、1.2652、1.1628、1.2977，将主要影响范围内的影响区间按照功能区相乘可得综合污染参考系数，交通区、工业区、公园绿地区、生活区、山区分别为3.8427、3.0508、1.2652*1.2977=1.6419、1.2652*1.1628=1.4712、1.2652。将此污染指数按照功能区与图11相乘，可得预测的十年后污染指数分布图：图21 预测10年后污染图可以发现，极值高达500倍，与图11的极值160倍相比增加了两倍多。表十：区域平均污染等级污染等级海洋生活区工业区山区交通区公园区无污染2894208621305416133轻度污染01414436818667中度污染015210448934970中重度污染0926341822546重度污染02753301381467159总计289486860327182643475可以发现，除了山区以外，剩下四个功能区的重度污染比率均超过30%，其中工业区和交通区的重度污染比率更是超过50%，环境形势变得十分严峻。6. 模型分析与评价我们对于问题中所采用的数学模型，一一进行分析评价。6.1. 插值模型采用v4插值和kriging插值两种插值方法，可以相对准确客观地反映出地质水平及元素浓度分布。6.2. 增加量比率模型能够相对合理地应用于对于土壤因为人类活动而造成了增加量变化，从而体现在污染程度上面。6.3. 相关性模型对于土壤元素按照两两之间的相关性聚类分析，为分析污染物及污染源提供了途径及依据。6.4. 回归模型对污染源附近作回归模型分析，有利于分析土壤重金属在各个功能区的分布和传播方式。6.5. 地质累计指数-内梅罗综合指数模型一种能比较客观地反映土壤污染情况的模型，侧重对于采样点综合情况的把握和分析。6.6. 单因子综合污染指数-内梅罗综合指数模型一种能比较客观地反映土壤污染程度的模型，侧重与对各个元素综合情况的把握和分析。7. 模型总结模型可能存在如下几个优缺点。优点：1、 结合青岛市地图，合理分析其功能区分布。2、 采用较为常用合理的kriging插值模型，成功给出重金属各项空间分布。3、 利用相关性算法，成功找到聚类分析的金属，得出的结果也相对准确。4、 采用两种方式验证问题（1）中的结果，各有侧重，斧正结果的准确性。5、 考虑青岛市的发展，成功预测十年后污染情况。缺点：1、 未考虑海拔、海洋对于土壤变化的影响。2、 在问题一中未引入内梅罗指数法，仅仅取平均数，导致可能结果对于某种由单一金属造成的污染未加以侧重。3、 预测时考虑面相对单一，考虑因素少。8. 参考文献1韩中庚.数学建模方法及其应用.北京:高等教育出版社,2005.6.2司守奎,孙玺菁.数学建模算法与应用.北京:国防工业出版社,2011.83张德丰.MAT