城市表层土壤重金属污染分析数学建模国家二等奖

2011高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：参赛队员(打印并签名)：1.2.3.指导教师或指导教师组负责人(打印并签名)：日期：年月日赛区评阅编号（由赛区组委会评阅前进行编号）：2011高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：城市表层土壤重金属污染分析摘要本文以某城区为例，对其土壤地质环境重金属含量进行取样分析，找出取样点被污染的主要原因并寻找到取样点的污染源。为此，对于第（1）和第（2）小题我们对取样点的数据统计分析，用MATLAB软件画出8张元素的分布图，各代表同一重金属元素在不同功能区不同含量范围内的分布点，统计出点的数目，对均值和方差进行对比得出该城区内不同区域重金属的污染程度和主要污染原因。针对第（3）小题，我们建立了2个模型来研究，干尘扩散模型——模型一，高斯水流扩散模型——模型二。模型一是为了求解重金属元素附着在干尘粒上的扩散方式和范围，运用烟雾扩散模型，结论是根据取样点的某金属含量程度来判断该点距污染源的距离，则污染源在以该距离为半径的圆的圆周上，在根据取样点的金属含量高低分布情况，找到污染源。模型二是研究重金属在土壤中的扩散情况，在相间物质交换为平衡的条件下,可用阻滞系数来表示其影响，可得到一个微分方程带入参数，可以找到某金属的污染源坐标。关键词：烟雾扩散模型图例统计污染物运移物质交换MATLAB1问题的提出随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。本文根据某个城市的5个区域的采样数据分析出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。通过数据分析，说明重金属污染的主要原因。分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求通过数学建模来完成以下任务：(1)给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2)通过数据分析，说明重金属污染的主要原因。(3)分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4)分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？2问题的分析对于问题和问题，主要采用图形对照法，根据附件中的数据使用MATLAB图形处理软件图形化，分析出同一重金属元素在不同功能区及不同金属元素在同一功能区的含量分布情况，找出污染的主要原因；而对于问题，为了更好的研究重金属传播特性和确定污染源位置，我们建立了两种模型：模拟重金属在空气中的传播，建立干尘扩散模型——模型一；考虑到金属在土壤里的传播，建立土壤扩散模型——模型二。在问题中，对模型的优缺点进行了评价，提出了其它在模型中没有考虑到的一些会影响重金属污染程度的因素，并在模型一和模型二基础上进一步优化。问题一：金属元素符号说明：砷镉铬铜汞镍铅锌重金属含量分段说明：、、、、、、、含量安全范围分别为：、、、、、、、，我们分段时从均值开始到安全范围上限，以后的范围断是人为分段（由于相关数据难以查到），具体分布可以参照表1经过参考附录中的坐标点及海拔我们用软件实现了该城区的三维模拟再根据已知数据用画出俯瞰的功能区分布图如图2示，经过图形处理，得到该市大概的地区分布图形如示：由砷的相关数据画出各个区砷含量分布：由和知砷在部分工业区含量在以上，属于重度污染，还有部分工业区砷含量在，交通区砷含量普遍在内，属于轻度污染；由镉的相关数据得出该城市镉的含量分布用比，可以得出该城市在山区镉含量为，属于轻度污染；在公园和生活区污染相对要严重，含镉量在；在交通区含镉量在；而工业区含镉量则在之间，污染相当严重。由铬元素的数据可以得到铬元素的分布：通过和对比，该城市普遍铬元素含量在，在生活区含量在，而工业区和交通区污染则更严重，含量在之间；由附件中铜元素的数据可以画出铜元素在该城市的分布：通过和的对比，发现在工业区和交通区铜元素含量在之间，其它功能区含量普遍低于；根据附件中汞元素的分布数据得出：与对比后可以看出工业区﹑交通区以及地处工业区附近的公园绿地区含量偏高，在之间，生活区和山区汞含量则低于；从附件中镍元素的数据得到了镍元素在该城区的分布：通过和比较，工业区和交通区以及生活区镍含量偏高，在之间，而其它功能区则普遍低于；同样的，由附件中的数据可以得到铅和锌元素在城市的分布和：通过和图3对比铅以及锌在交通区﹑业区以及生活区含量高，铅含量在80μg/g～480μg/g之间，锌含量则是500μg/g～3800μg/g的范围内；上述图对比为宏观上的判断，针对此，我们列出了同一重金属元素在不同含量段在整个城区的分布，并统计出每个含量段的取样点数，如表1As(μg/g)Cd(ng/g)Cr(μg/g)Cu(μg/g)含量范围样点数含量范围样点数含量范围样点数含量范围样点数3.6-5.488130-1906431-4914813.2-20.4585.4-1513190-50015949-1009920.4-15021115-30.135500-100039100-2005150-500101000-20007200-5003500-10000500-100021000-26002Hg(ng/g)Ni(μg/g)Pb(μg/g)Zn(μg/g)含量范围样点数含量范围样点数含量范围样点数含量范围样点数35-516112.3-19.916831-438669-977151-1509719.9-507043-8010997-16068150-5004150-100180-15053160-50094500-10007100-1501150-2509500-1500131000-16003250-48031500-38005表1通过表1的分析，、、、、在该城区属于重度污染；、达到污染；分析，、、、、在该城区属于重度污染；、达到污染；属于轻度污染综上，得出重金属分布情况：八种重金属元素在工业区、交通区含量较高，污染较严重；在生活区和公园绿地区含量相对要少，属于轻度污染；而山区各种金属元素含量除了接近工业区和交通区的区域，重金属含量都在安全范围内。问题2：主要思路：分别对五个功能区重金属元素的分布情况进行分析，找出该区某重金属含量超过安全范围的主要原因。通过均值和偏差的对比得讨论如下：功能区均值砷（As）镉（Cd）铬（Cr）铜（Cu）汞（Hg）镍（Ni）铅（Pb）锌（Zn）生活区7.48502.97118.16174.8953.6623.37116.22501.39工业区4.21323.8350.7950.4778.9814.6960.96264.86山区5.11224.9539.5925.0890.5214.9444.81104.51交通区5.83286.7243.9641.61312.8117.7157.92162.75公园绿地区5.39236.0538.9118.4147.4614.8741.0592.31标准偏差砷（As）镉（Cd）铬（Cr）铜（Cu）汞（Hg）镍（Ni）铅（Pb）锌（Zn）生活区3.93281.54166.97415.523081.3519.4996.14547.46工业区1.48162.9745.2452.79146.886.1225.11607.79山区3.71249.7422.1321.75261.599.7933.31165.03交通区2.46173.1614.7829.891719.685.3531.21159.96公园绿地区2.39189.9016.338.0448.645.8812.4267.64对照图和表所列数据，得以下的结论：1、生活区铅和锌的含量相对其他重金属含量较高，考虑到人们的生活环境，如人们所用到的家具和房间的涂料，以及家用汽车的增加导致尾气的排放、生活垃圾等都含有大量的铅；含锌化合物主要用作油漆、油墨、塑料、广告颜料、化妆品、蓄电池的原料或填充剂。所以可以看出该市的人们在生活区使用的涂料、化妆品、电池和相关商品较多并且处理不的不完善，导致含锌量较高。2、工业区所有重金属的含量都很高。长期以来我国由于科技和经济欠发达，居民的环境保护意识薄弱，工业生产过程中产生的废水、废气、废渣未经处理直接排放。其中含有大量的重金属，经过自然的沉降雨水的淋溶等途径进入土壤导致Hg、Pb、Cu、Cr、As、Ni、Zn、Cd含量逐年增加导致工业区所有重金属含量都很高。3、山区各重金属含量大部分都在正常范围内，部分地区属于轻度污染。主要原因是受到工业区和生活区的影响。工业区通过空气传播和水流传播以及生活区的生活垃圾污染导致部分山区的轻度污染。4、交通区铅和镉、铜、锌的含量较高。交通区汽车燃油、润滑油和镀金部分的燃烧或磨损释放出大量重金属元素。一般以道路为中心成带状分布，导致路边土壤中这几种重金属元素明显高于非交通线路两侧的土壤。5、公园绿地区汞、锌、铅的含量相对较高。首先考虑到该区所处的地理位置,地处生活区、工业区、交通区附近。所以公园较严重的重金属污染可能来源于交通、工业污染和居民生活废弃物等污染。问题3：一模型的假设模型一干尘扩散模型1、假设无风情况下，重金属看作在地面某一点向四周等强度的在无穷空间扩散，不计地形影响。2、重金属物质附着在干沉上的扩散服从扩散定律，即单位时间通过单位法平面面积的流量与它的浓度梯度成正比。3、重金属物质在穿过降雨时，其强度由于雨水的吸收而减少，减少比率为常数。4、风的影响只考虑重金属污染源的某一风向、某一风速的影响。风的影响范围，随着时间而减小。模型二土壤水流扩散模型1、重金属在土壤多孔介质中有三种存在形式，假设相间物质交换为平衡2、延滞系数唯一。二、模型的建立1.模型概述重金属扩散的范围和强度主要受到污染源的自然扩散、风以及降雨的影响。自然扩散符合烟雾扩散模型，在给定污染源总含量等数据时，可以求出静态扩散范围。风力具有加大扩散面积的作用，而降雨具有吸收稀释重金属元素含量的作用。2.参数示意污染源浓度污染源重金属扩散的速度表示某点的浓度污染源重金属扩散的总量扩散的时间降雨对重金属的吸收率有风时候的速度有风时候的方向风力影响随时间减小量泄漏点距有效地面的高度(m)分别为用浓度标准差表示的轴上的扩散参数3.模型的求解求解主要思想：在无降雨情况下污染源重金属总量等于扩散后的总量，所以有，当有降雨的情况下，扩散后金属含量有变化，吸收后剩余量为.3.1模型一建立过程建立空间坐标系：将重金属开始扩散时设为，以污染源为坐标原点，建立以正东为x正方向，正北为y正方向，建立三维直角坐标系，如下：3.1.1推导重金属随干沉扩散范围半径（无风无降雨的扩散模型）时刻t无穷空间中任意一点的重金属含量记为P。根据假设2，单位时间通过单位法相面的流量为：·····················s是扩散系数，grad表示梯度，负号表示由浓度高向浓度低的地方扩散。考察空间域Ω，Ω的体积为V，包围Ω的曲面为S，S的外法线向量为n,则在[t,t+Δt]内通过Ω的流量为：··················而Ω内重金属物质含量的增量为：·····················由扩散前后金属总量不变得：································根据曲面面积分的奥氏公式：·······················其中是散度记号。由可得：无界区域的抛物型偏微分方程······················根据假设1，初始条件为作用在坐标原点的点源函数，可以记作：···························得：·······这个结果表明。对于任意时刻t，重金属物质浓度C的等值面是球面x2+y2+z2=R2,并且随着球面半径R的增加C的只是连续减少的；当R—>∞或t—>∞时，C—>0通过将x2+y2+z2=R2代入方程（8）可以解出R：··············R为在给定扩散系数s，影响时间t，重金属含量为C，以及污染源总含量W的情况下的受自然扩散影响范围的半径，即在不考虑风力、降雨情况下的影响范围。通过matlab，可以实现该模型的图像生成，如图一。在此假定：s=0.00001(扩散系数)，W=3000（μg/g）（污染源总量，假定），C=103.20（重金属在检测点物质浓度，这里取铅），t=1000000(秒)（扩散时间）r=0.5(假设雨水对于重金属吸收率为50%)ws=0.04(公里每小时)（假设作用在扩散源的风速为40米/小时）wa=0.785（假设风向与x0y面平行，且与x轴成正45度角）wd=0.3（假设风力影响范围从R开始，每天减弱300米）我们调用函数ZJ_Shu(k,Q,C,t,r,wv,we,wj)（函数本身可参考附录），加入参数为ZJ_Shu(0.0001,3000,103.2,1000000,0.5,0.04,0.785,0.3)，符号意思看参数说明表，结果为：经过去不同重金属的均值，我们得到元素取值（平均值）半径RAs3.624.9564×10^3Cd13015.3267×10^3Cr3118.1306×10^3Cu13.220.4128×10^3Hg3517.8477×10^3Ni12.320.6261×10^3Pb3118.1306×10^3Zn6916.4310×10^33.1.2考虑风和降雨影响的模型建立降雨因素：假设降雨对于放射性物质的吸收率为r,则方程（9）变为：风力因素：只考虑产生在辐射原点某一角度，某一强度的风力影响。风力的影响范围从R开始每天减少，并且对于影响范围内每一点的作用相同。3.2模型结果分析模型的分析结果说明在给定的条件下，扩散的最远距离为24.95公里。由于现实中风向的不确定性，出于安全考虑，在污染源周边24.95公里范围的地域内都有污染，而根据距离，可以判断出污染源在以距离为半径的圆周上，寻找到工业区或者是生活区。模型二：重金属元素在土壤中的传播模型分析：重金属元素在土壤多孔介质中有三种可能的存在形式：:溶于水中、挥发为气体及吸附于固体颗粒。挥发性有机污染物在水、气、固体颗粒三相间的物质交换与分配是决定其运移的重要因素,在相间物质交换为平衡的条件下,可用阻滞系数来表示其影响。污染物与土壤固体颗粒间的非平衡吸附解吸是相间交换中影响污染物运移的最重要的机制,由于天然土壤具有固有的不均质性,必须用多个反应系数才能准确描述污染物与固体颗粒间的非平衡吸附解吸。在土壤环境中,一系列的机制控制着污染物的运移:一、地下水流决定了污染物的运动方向和速率;二、扩散使污染物产生纵向及横向的转移;三、污染物与土壤颗粒中有机质及矿物质之间的吸附解吸、污染物在土壤包气带中的水气界面处的物质交换使污染物的运移受到阻滞作用;四、由于具有挥发性,污染物还随气体迁移、扩散;五、土壤中的生物与化学作用使污染物降解,生成无害物质或其他有害物质。要预测污染物的运移和其归宿,必须对土壤—污染物—空气—水这一复杂的相互作用的系统及污染物在其中的诸多迁移机制有充分理解,并把该系统模型化从而得到污染源的位置。模型的建立符号说明气体中的污染物浓度亨利系数水中的污染物浓度颗粒表面的污染物浓度污染物在水与固体颗粒间的分配系数分配系数土壤含水率土壤含气率土壤干容重孔隙水流速度水平距离扩散系数比例系数距离传播速度时间某已知点的横坐标某已知点的纵坐标污染源的横坐标污染源的纵坐标1气—水—固三相间的物质交换在有水气共存的土壤中,即土壤非饱和带中,污染物随水运移时,水中的污染物会与气体及土壤颗粒发生物质交换。水气两相间的物质交换平衡可用亨利定律来表示式中,为气体中的污染物浓度,为水中的污染物浓度,H为亨利系数。不同物质的H值可在文献中查到。对于水中污染物与固体颗粒间的物质交换,经常用到以下这一简单的线性方程:=式中,为固体颗粒表面的污染物浓度,为污染物在水与固体颗粒间的分配系数。有机物质在水与固体颗粒间的分配与土壤颗粒中的有机质含量、土壤颗粒的大小与表面积等因素直接相关,而在天然土壤中,土壤的这些性质往往不是均一的,用线性的方程来描述水—固间的物质分配与实测数据往往不很吻合,而非线性方程则更好地表示了其间的物质交换。方程是一个经常被使用的非线性吸附等温线方程,其表达式为:式中,为分配系数,表示水—固吸附非线性指数。方程、都属经验公式,其中的参数值可通过专门的实验测出。污染物由气相吸附在固体颗粒表面的量一般很小,可以忽略不计。假设气相为静止,且不考虑污染物在气相中的扩散,污染物在非饱和带中随水运移可由下列一维对流—扩散方程来表示:式中,，分别为土壤含水率与含气率,为土壤干容重,为孔隙水流速度,为水平距离,为扩散系数。方程左边第二、第三项分别表示污染物在固相与气相中的分配。在水—气、水—固之间的物质分配为平衡的条件下,可把方程、(或)代入,且方程两边同除以得:为延滞系数,延滞系数Rt代表污染物运移时由于发生了水—气及水—固体颗粒之间的物质交换而产生的延滞作用。在不同的土壤含水率条件及不同的污染物运移承载相(气或水)情况下,延滞系数的表达式是不同的,如考虑污染物随气相在非饱和带中运移,则可由污染物通过土柱的时间与不发生水—固或气—固吸附时非挥发性试剂通过土柱的时间的对比实验来求得。近来有实验表明,实验测得的污染物在非饱和土柱中随气相运移时的延滞系数大于由求得的预测值。这部分额外的延滞效应则应归于污染物在水—气界面上的积累。污染物在水—气界面上的积累不仅使污染物运移受阻,而且其影响还相当大,如等在实验中发现(三氯乙烯)的滞后是由于在水—气界面处的积累造成的。因此,挥发性有机污染物在水—气—固三相的分配可能包括以下几个方面:污染物在水—气两相间的分配、污染物在水—固体颗粒两相间的分配及污染物在水—气界面处的累积。同时考虑这三个方面后的污染物随气体运移的延滞系数为:式中,为污染物在水气界面处累积的吸附系数,为有效水气界面面积。由所求出的延滞系数可以看出和即浓度有关，我们可以假设延滞系数即扩展系数，显然延滞系数和传播速度成正比而从以污染源为起点延滞系数将逐渐降低，设污染源的位置为任意一地点（采样的点）位置为。单位时间传播的距离的积分就是污染源到这点的距离从而由得根据和已知数据便可以求出污染源的位置。用MATLAB求解并取代表性坐标如元素坐标As(μg/g)（5013，7202）Cd(ng/g)（2724，2576）（4034，6005）Cr(μg/g)（4632，4315）（3504，6013）Cu(μg/g)（2367，3443）（3845，6012）Hg(ng/g)（4812，1845）（4865，7011）Ni(μg/g)（4055，6089）Pb(μg/g)（2002，2951）（4990，5310）Zn(μg/g)（14024，9529）（3970，5978）（5372，7537）经过数据分析，主要污染源在该城区的西南方向，集中了工业区，交通区，生活区，海拔在4——8之间。问题44模型的评价对于干尘扩散模型——模型一优点:考虑了无风无雨时候的理想情况和有风有雨的情况，把现实生活中的大部分情况都包括在内，并且得出了不同情况时候的传播关系。缺点：考虑的因素较少，忽略了人为，湿度，自然环境的吸收能力等因素。得到的结果的准确性自然较低。对于土壤扩散模型——模型二优点：通过文献，得到了直观且简单的计算方法和关系式。缺点：也只简单的求出延滞系数，通过所求的式子得到各个量之间的关系，得到浓度和位置的一个关系式。对不同的重金属的传播途径的分析和利用不够透彻。为更好地研究城市地质环境的演变模式，还应收集的信息有:本地区不同地方的降水量、地基的复杂程度、人们的活动范围以及将来社会的发展趋向（如哪里新建住房等）等。有了以上这些信息。我们认为应综合的根据以下各方面的原因分析。（1）降雨——降雨的多少必定影响这个地区的土壤和污染等；（2）地基的复杂程度——地基的不同也对这个地区的承受能力和变化有影响；（3）人们的活动范围——人们的活动对这个地方的演变起到重要的作用，人多的地方演变的必定更快。（4）将来社会的发展趋向——社会的发展趋向对环境的影响更是至关重要，直接影响人们活动的地域，从而影响环境。5模型的推广1、本模型可以扩展到其他污染物的污染，只需检测出相应相应的污染物在各个区域的含量分布，就可以找出污染源；2、本模型对扩散型污染物都有一定的适用性，可以参照此解决一些扩散性污染的问题；3、本模型具有较强的应用价值。环境的污染容易破坏生态链，从而影响人类的生活，重金属的污染，甚至可以直接危害人类的健康。该模型可以通过控制一些因素找出被污染城市的污染源，从而及时制定应急措施，减小损失。参考文献杨桂元，数学模型应用实例，合肥：工业大学出版社，2007年6月第一版汪国强，数学建模优秀案例选编，华南：华南理工大学出版社，2001年8月第一版周博谢东来，MATLAB科学计算，三河市：机械工业出版社，2010年5月第一版附录：图形实现主代码A=[......];%此矩阵数据太多，在此省略bianhao=A(:,1);%编号As=[A(:,2)];Cd=[A(:,3)];Cr=[A(:,4)];Cu=[A(:,5)];Hg=[A(:,6)];Ni=[A(:,7)];Pb=[A(:,8)];Zn=[A(:,9)];x=[A(:,10)];y=[A(:,11)];z=[A(:,12)];gnqu=[A(:,13)];X_hang=length(gnqu);%scatter(x,y,5,z)%散点图[X,Y,Z]=griddata(x,y,z,linspace(0,28654)',linspace(0,18449),'v4');%插值%pcolor(X,Y,Z);%shadinginterp%伪彩色图%所有监测点网络图fig1=figure(1);set(fig1,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonplot(x,y,'-*')title('检测点网络图');xlabel('x');ylabel('y');fig2=figure(2);set(fig2,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdon%Ask=0;fori=1:X_hangif3.6<=As(i)&As(i)<=5.4plot(x(i),y(i)+250,'^');holdonelseif5.4<As(i)&A(i)<=15plot(x(i)+250,y(i),'s','markerfacecolor','g');holdonelseifAs(i)>15plot(x(i),y(i)-250,'d','markerfacecolor','r');text(x(i)+200,y(i)+200,sprintf('(%0.1f)',As(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendtitle('As的分布','backgroundcolor','y','color','b','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');text(1000,17900,'

-重度污染As(μg/g)>15','backgroundcolor','k','color','w');text(1000,16900,'△-轻度污染3.6<=As(μg/g)<=5.4','backgroundcolor','k','color','w');text(1000,15900,'□-污染5.4<As(μg/g)<=15','backgroundcolor','k','color','w');%Cdfig3=figure(3);set(fig3,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonfori=1:X_hangifCd(i)<130&Cd(i)<190plot(x(i),y(i),'v','markerfacecolor','k');%text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','k','fontsize',10,'fontweight','bold');holdonelseif190<=Cd(i)&Cd(i)<=500plot(x(i),y(i),'p','markerfacecolor','g','markersize',10);%text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','k','fontsize',10,'fontweight','bold');holdonelseif500<Cd(i)&Cd(i)<=1000plot(x(i),y(i),'o','markerfacecolor','r');%text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','r','fontsize',10,'fontweight','bold');holdonelseif1000<Cd(i)&Cd(i)<2000plot(x(i),y(i),'d','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cd(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendtitle('Cd的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');text(100,17900,'▲-130<Cd(ng/g)<190','color','w','backgroundcolor','k');text(100,16900,'★-190<=Cd(ng/g)<=500','color','w','backgroundcolor','k');text(100,15900,'●-500<Cd(ng/g)<=1000','color','w','backgroundcolor','k');text(100,14900,'◆-1000<Cd(ng/g)<2000','color','w','backgroundcolor','k');%Crfig4=figure(4);set(fig4,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonfori=1:X_hangif31<=Cr(i)&Cr(i)<=49plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif49<Cr(i)&Cr(i)<=100plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif100<Cr(i)&Cr(i)<=200plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif200<Cr(i)&Cr(i)<=500plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif500<Cr(i)&Cr(i)<=1000plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cr(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-32<Cr(μg/g)<49','color','w','backgroundcolor','k');text(100,16900,'★-49<=Cr(μg/g)<=100','color','w','backgroundcolor','k');text(100,15900,'●-100<Cr(μg/g)<=200','color','w','backgroundcolor','k');text(100,14900,'◆-200<Cr(μg/g)<500','color','w','backgroundcolor','k');text(100,13900,'■-500<Cr(μg/g)<1000','color','w','backgroundcolor','k');title('Cr的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Cufig5=figure(5);set(fig5,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonfori=1:X_hangif13.2<=Cu(i)&Cu(i)<=20.4plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif20.4<Cu(i)&Cu(i)<=150plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif150<Cu(i)&Cu(i)<=500plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif500<Cu(i)&Cu(i)<=1000plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif1000<Cu(i)&Cu(i)<=2600plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Cu(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-13.2<Cu(μg/g)<20.4','color','w','backgroundcolor','k');text(100,16900,'★-20.4<=Cu(μg/g)<=150','color','w','backgroundcolor','k');text(100,15900,'●-150<Cu(μg/g)<=500','color','w','backgroundcolor','k');text(100,14900,'◆-500<Cu(μg/g)<1000','color','w','backgroundcolor','k');text(100,13900,'■-1000<Cu(μg/g)<2600','color','w','backgroundcolor','k');title('Cu的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Hgfig6=figure(6);set(fig6,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonfori=1:X_hangif35<=Hg(i)&Hg(i)<=51plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif51<Hg(i)&Hg(i)<=150plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif150<Hg(i)&Hg(i)<=500plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif500<Hg(i)&Hg(i)<=1000plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif1000<Hg(i)&Hg(i)<=1600plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Hg(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-35<Hg(ng/g)<51','color','w','backgroundcolor','k');text(100,16900,'★-51<=Hg(ng/g)<=150','color','w','backgroundcolor','k');text(100,15900,'●-150<Hg(ng/g)<=500','color','w','backgroundcolor','k');text(100,14900,'◆-500<Hg(ng/g)<1000','color','w','backgroundcolor','k');text(100,13900,'■-1000<Hg(ng/g)<=1600','color','w','backgroundcolor','k');title('Hg的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Nifig7=figure(7);set(fig7,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonfori=1:X_hangif12.3<=Ni(i)&Ni(i)<=19.9plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif19.9<Ni(i)&Ni(i)<=50plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif50<Ni(i)&Ni(i)<=100plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif100<Ni(i)&Ni(i)<=150plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Ni(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendtext(100,17900,'▲-12.3<Ni(μg/g)<19.9','color','w','backgroundcolor','k');text(100,16900,'★-19.9<=Ni(μg/g)<=50','color','w','backgroundcolor','k');text(100,15900,'●-50<Ni(μg/g)<=100','color','w','backgroundcolor','k');text(100,14900,'◆-100<Ni(μg/g)<150','color','w','backgroundcolor','k');title('Ni的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Pbfig8=figure(8);set(fig8,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonfori=1:X_hangif31<=Pb(i)&Pb(i)<=43plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif43<Pb(i)&Pb(i)<=80plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif80<Pb(i)&Pb(i)<=150plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif150<Pb(i)&Pb(i)<=250plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif250<Pb(i)&Pb(i)<=480plot(x(i),y(i),'s','markerfacecolor','r','markersize',12);text(x(i)+200,y(i)+200,sprintf('(%0.1f)',Pb(i)),'color','r','fontsize',10,'fontweight','bold');holdonendendendendendendtext(100,17900,'▲-31<Pb(μg/g)<43','color','w','backgroundcolor','k');text(100,16900,'★-43<=Pb(μg/g)<=80','color','w','backgroundcolor','k');text(100,15900,'●-80<Pb(μg/g)<=150','color','w','backgroundcolor','k');text(100,14900,'◆-150<Pb(μg/g)<250','color','w','backgroundcolor','k');text(100,13900,'■-250<Pb(μg/g)<=480','color','w','backgroundcolor','k');title('Pb的分布','color','b','backgroundcolor','y','fontsize',12,'fontweight','bold');xlabel('x');ylabel('y');%Znfig9=figure(9);set(fig9,'menubar','none','position',[2551001000600]);contourf(X,Y,Z);%等高线图holdonfori=1:X_hangif69<=Zn(i)&Zn(i)<=97plot(x(i),y(i),'v','markerfacecolor','k');holdonelseif97<Zn(i)&Zn(i)<=160plot(x(i),y(i),'p','markerfacecolor','g');holdonelseif160<Zn(i)&Zn(i)<=500plot(x(i),y(i),'o','markerfacecolor','b');holdonelseif500<Zn(i)&Zn(i)<=1500plot(x(i),y(i),'d','markerfacecolor','r','markersize',10);holdonelseif1