全国大学生数模竞赛论文城市表层土壤重金属污染分析

城市表层土壤重金属污染分析摘 要随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。本文针对某城市地表土壤重金属测量数据，开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式。问题1：对8种重金属元素进行空间散乱点插值和Kriging插值，利用ArcGIS地统计软件得到8种重金属的空间分布图，运用单项污染指标法和内梅罗综合污染指标法分别得出每一种重金属在不同功能区的污染等级：AsCd CrCuHgNiPbZn一区轻污染中污染中污染重污染中污染轻污染中污染重污染二区中污染重污染轻污染重污染重污染轻污染重污染重污染三区轻污染轻污染轻污染轻污染轻污染轻污染轻污染轻污染四区轻污染轻污染轻污染重污染重污染轻污染中污染中污染五区轻污染轻污染轻污染中污染重污染轻污染轻污染轻污染综合污染等级如下：功能区污染指标污染等级生活区3.4349重污染工业区4.0279重污染山区1.0622轻度污染交通区3.5196重污染公园绿地区2.2354中度污染问题2：对8种土壤重金属元素含量运用相关性分析，主成分分析和聚类分析发现，该市土壤中重金属污染主要受铁路交通，农业农药塑料薄膜使用，污灌，生活垃圾堆放，采矿这六个方面影响问题3：重金属的传播特征一般是离污染源越远，污染程度越小。将问题1中空间分布图对应在该区域的所有插值点作为样本点，建立层次规划模型，目标层为整个区域的某元素加权平均浓度，准则层为该元素污染源，决策层为该元素污染源内样本点的浓度，建立非线性规划模型，得到两组最优权重表示污染源内的点在该污染源中所占比例，表示污染源内的点对样本点加权浓度的影响度。根据越小，污染源对加权平均浓度的影响力越大，对As进行模拟，得到，可以判断显然为污染源，是污染源的可能性很小。越大，其对应点的As浓度越大，越接近污染源的中心。问题4：将单一的测定土壤中的重金属浓度，扩展为测定空气，水域，植物和土层（050cm）的各重要元素浓度，并将取样间距缩小,定期对其取样，将m种元素n个时期的元素浓度均值作为研究对象，以每次检测时各元素污染指数均值作为参考数列，进行灰色关联度分析，得到各个时期主要污染元素的演变；并对每一种元素浓度污染指标进行灰色预测，得到各元素的进入土壤和在土壤中的传播规律及演变模式，结合地理统计软件ArcGIS，画出各元素在不同介质中的空间分布图，对各图进行比较，分析，全面地分析城市地质环境的演变模式。 关键词：地统计学，内梅罗综合污染评价，多元分析，灰色预测 一、问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（010 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求通过数学建模来完成以下任务：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、模型假设1. 假设该城市土壤，地质环境短时间内不发生变化。2. 假设采集回的样本能很好的反应采样点的实际情况。3. 假设题目所给数据真实准确。4. 假设在土壤采样，化验，分析过程中未出现任何操作错误或条件改变。5. 假设采样过程中未出现特大严重重金属污染情况三、问题分析根据题目所给信息，对城市表层土壤重金属污染进行分析讨论后，解决了如下四个问题： （1）土壤重金属空间分布是应用地统计学方法对变量空间结构的表达与模拟，为了更好的说明土壤中重金属的分布情况，首先要对数据进行预处理，由于在土壤采样，化验，分析过程中会出现操作错误或条件改变，这将会造成变量Kriging插值连续表面中断，使得实验半方差函数1发生畸变，甚至会掩盖变量固有的空间结构特征，因此首先要剔除异常值，然后进行变量正态分布检验，然后计算半方差函数获取土壤变量空间自相关性，最后用Kriging插值法2获取土壤变量的空间分布图；应用综合污染指数法，能全面综合地反映土壤的污染程度，所以要分析城区内不同区域重金属的污染程度，需采用兼顾单元素污染指数平均值和最大值的尼梅罗综合污染指数法3，计算公式如下： （2）为了说明重金属污染的主要原因，首先要对这8种重金属元素进行相关性分析，以确定它们的污染来源是否相同，或具有关联出现，然后再对该8种重金属元素进行主成分分析，考虑占主导地位的成分的影响，使具有强相关性的元素归并为一类。最后在前面的多元分析下说明重金属污染的主要原因。 （3）重金属的传播特征一般是离污染源越远，污染程度越小。由于同一区域土壤中重金属污染的来源可以是单一的，也可以是多途径的，因此由第（1）问各元素在不同功能区浓度的空间分布与第（2）问中对各元素污染原因分析，可以确定各污染源的大概位置以及来自于同一污染源的元素种类(即如果污染点和污染源均相同，则可以得出这些元素来自同一污染源)，得出8种元素的可能污染源位置和来自同一可能污染源的元素种类，在个可能污染源周围划分一个面积为的区域，将问题（1）中空间分布图对应在该区域的所有插值点作为样本点（以下称样本点），建立以整个区域的某元素加权平均浓度为目标层，每种元素污染源为准则层，每种元素污染源内样本点的浓度为决策层的非线性规划模型。（4）研究城市地质环境的演变模式是个复杂的系统过程，仅凭土壤中的重金属浓度进行分析远远不够，为此，我们将单一对土壤中重金属浓度的测定，扩展为测定空气，水域，植物和土层（050cm）的各重要元素浓度，并将取样间距缩小，定期对其取样，记录。运用地理统计软件ArcGIS，画出各元素在不同介质中的空间分布图，对各图进行比较，分析；将m种元素n个时期的元素浓度均值作为研究对象，以每次检测时各元素浓度的均值作为参考数列，进行灰色关联度分析；对每一种元素浓度污染指标进行灰色预测，分析结果，得到各元素的进入和传播规律及这种方式转移情况，并得出各元素土层分布规律，从而全面地分析了城市地质环境的演变模式。四、符号说明定义符号说明半方差函数样本间距间距为h的“样本对”对数样本点的空间坐标样本点出的元素浓度值插值点处的估计值污染物的污染指数污染物的实测值的平均值污染物的评价标准单项污染指数法评价标准重金属元素j的污染指数单元素污染指数最大值单元素污染指数平均值普通Kriging系数样本数样本点的加权平均浓度元素第i个污染源t污染源个数第i个污染源中样本点的个数I整个区域内插值点的个数污染源内的点对样本点加权浓度的影响度污染源内的点在该污染源中所占比例名词解释：1.城市土壤：指由于人为的、非农业作用形成的，并且由于土地的混合、填埋或污染而形成的厚度大于或等于50cm的城区或郊区土壤2.地质环境：指地壳上部包括岩石指地壳上部包括岩石、水、气和生物在内的互相关联的系统五、模型建立与求解5.1问题（1）模型的建立与求解5.1.1重金属元素空间分布模型的建立我们根据题意在解决重金属元素空间分布时建立了两个模型：模型一：根据给出数据，进行散乱点插值，利用matlab做出三维坐标图（程序见附录一）。一般情况下，x，y，z是维数相同的矩阵，x，y是网格坐标矩阵，z是网格点上的高度矩阵，c用于指定在不同高度下的颜色范围。c省略时，Matlab认为c=z，也即颜色的设定是正比于图形的高度的，但是此题图像中颜色的深浅用于标定金属元素的浓度。（红黄蓝绿紫），图（1）到图（8）分别表示As,Cd ,Cr,Cu,Hg,Ni,Pb,Zn在空间的浓度大小。 图（1） 图（2） 图（3） 图（4） 图（5） 图（6） 图（7） 图（8）模型二：（1）为使Kriging插值表面连续，要剔除异常值，利用matlab（程序见附录二），首先要对数据经行预处理（见附件1）。（2）正态分布检验 利用spss数据统计分析功能，通过P-P图检验数据是否符合指定的分布，当数据符合指定的分布时，P-P图中各点近似呈现一条直线，并且我们对于P-P图中各点不成直线但有一定规律的，对其数据进行转换，使转换后数据更接近指定分布。对各重金属含量概率服从做如下表（1）所列假设： 表（1）假设8种重金属元素服从对数正态分布重金属元素AsCdCrCuHgNiPbZn概率分布对数对数对数对数对数对数对数对数得到各重金属元素P-P图如下所示：图（9）8种重金属元素的P-P图因此题目所给数据符合指定的对数正态分布（3）半方差函数分析首先建立如下所示的半方差函数： (1)其中是样本间距，又称为位差；是间距为的“样本对”对数，为空间坐标，为处的重金属浓度。然后我们又利用了Kriging插值，插值函数如下： (2)其中为处得估计值，为处得观测值，为普通Kriging系数，为观测点个数。模型一是基于三维非网格形式的插值拟合，对于土说明壤中重金属的空间分布误差较大，所以我们选择模型二，Kriging插值，较能说明空间分布。5.1.2重金属元素空间分布模型的求解结合多元统计学和GIS技术利用ArcGis地统计学软件对土壤的属性数据进行数学分析和空间制图可以很直观的了解研究城市土壤重金属元素的来源和空间分布特征如下图所示（以As元素为例，其它元素分布图见附录三）： 图（10） As元素空间分布图 5.1.3重金属污染程度模型的建立与求解污染的评价方法：（1） 单项污染指数法 (3) 其中为污染物的污染指数，为污染物的实测值的平均值，即对本题所给附表二中给出的数据进行处理，为污染物的评价标准单项污染指数法评价标准，即本题中附表三给出的背景值：非污染，轻污染，中污染，重污染。（程序见附录四）表（2）各区域8种重金属污染程度 元素 区域As,Cd ,Cr,Cu,Hg,Ni,Pb,Zn一区轻污染中污染中污染重污染中污染轻污染中污染重污染二区中污染重污染轻污染重污染重污染轻污染重污染重污染三区轻污染轻污染轻污染轻污染轻污染轻污染轻污染轻污染四区轻污染轻污染轻污染重污染重污染轻污染中污染中污染五区轻污染轻污染轻污染中污染重污染轻污染轻污染轻污染数据分析：将每种金属元素在五个区中污染最严重的用红色标出，可以看出二区红色最多，因为二区是工业区，这符合实际。而三区没有一个红色，因为三区是山区，离污染源最远，这也符合实际。（2）内梅罗综合污染指标法单因子污染指数法只能分别反应各个污染物的污染程度，不能全面，综合地反应土壤的污染程度时，需将单因子污染指数按一定方法综合起来进行评价，即应用综合污染指数法评价。重金属元素中和污染评价采用兼顾单元污染指数平均值和最大值的内梅罗综合污染指数法。计算公式如下： ， (4) (5)式中： 为监测点单项污染指标平均值中的最大值，为检测点单项污染指标平均值（是基于单项污染指数法中的八种元素各自的平均值在对其求和除以八得出的结果）， 为单项指标的背景值，为监测点的综合污染指标，为j监测点所有污染物单项污染指数中的最大值与背景值的商；为j检测监测点所有污染物单项污染指数的平均值与背景值的商。表（3） 土壤内梅罗指数评价标准 等级 内梅罗污染指数 污染等级 一 清洁（安全） 二 尚清洁（警戒线） 三 轻度污染 四 中度污染 五 重污染根据公式，计算五个区的的内梅罗指数分别为（程序见附录五）： 生活区： 3.4349(重污染) 工业区： 4.0279(重污染) 山区： 1.0622(轻度污染) 交通区： 3.5196(重污染) 公园绿地区：2.2354(中度污染) 由数据分析可以看出对于污染指标工业区指标值最大，而山区的污染指标值最小，这种情况也符合实际，各区各重金属元素污染的分布图（见附录五）5.2问题（2）模型的建立与分析 为说明土壤中重金属污染的主要原因，先要对数据的统计学特征进行初步了解，通过描述统计分析和相关性分析，能进一步弄清土壤重金属含量数据的表观特征，为其污染的主要原因的说明提供一定的依据和理论支持。各重金属元素的背景值如下表（4）所示： 表（4）8种重金属元素的背景值元素平均值标准偏差范围As (g/g)5.4Cd (ng/g)1303070190Cr (g/g)3191349Cu (g/g)20.4Hg (ng/g)3581951Ni (g/g)19.9Pb (g/g)3161943Zn (g/g)691441975.2.1土壤重金属含量的描述统计分析特征描述统计分析主要研究如何将实验或调查得到的大量数据简缩成有代表性的数据，使其能客观，全面地反映这组数据的全貌，可以将其所提供的信息充分显现出来，为下一步统计分析提供依据，对该市319个城市土壤样品8种重金属含量进行统计分析，结果见表（5）表（5）城市土壤重金属含量统计描述元素样本数最小值最大值均值标准差偏度峰度变异系数（%）As3191.6130.135.67653.02433.32419.69653.28Cd319401619.8302.3962224.98762.0245.64174.4Cr31915.32920.8453.509770.00189.445104.166130.82Cu3192.292528.4855.0167162.915112.755180.685296.12Hg3198.5716000299.71131629.53988.59674.148540.7Ni3194.27142.517.26189.94147.10981.84457.59Pb31919.68472.4861.740950.05784.40628.30881.08Zn31932.863760.82201.2026339.23256.4953.591168.6 与该市土壤环境背景值进行比较，可以发现，8种重金属元素在该地区的平均含量分别是其自然环境背景值的1.58,2.33,1.73,4.17,37.46,4.54,10.29,14.07倍，都具有明显的累积现象；变异系数（C-V）是表征样本间变异程度的重要参数，C-V30%为强变异。该市城市土壤样品中8种重金属含量差异很大（见表（5），变异系数排序依次为：AsNiCdPbCrZnCuHg，这些重金属的变异系数均高于30%，变异系数越大说明土壤重金属污染受人为活动干扰越强烈，或者可以理解为污染程度越大，从变异系数中可以看出，Hg含量变异最为明显，其次为Cu，Zn。5.2.2土壤重金属含量之间的统计分析 研究土壤重金属含量之间的相关性，在一定程度上可反映这些元素污染程度的相似性或污染元素有相似的来源，利用spss统计分析软件包，对该市城市土壤和重金属含量做pearson相关分析，主成分分析及聚类分析。（1）相关性分析：一般情况下，如果元素元素含量间显著相关，说明他们来源相同的可能性较大，这一来源既有可能出自天然，也有可能来自人类活动造成的符合污染，如表（6）所示，由表可知，城市土壤重金属除As-Hg外其余其余都具有较高的相关性，表明它们之间的同源性很高，而As-Hg的相关性很弱。表（6）重金属含量的pearson相关矩阵（n=319）元素AsCdCrCuHgNiPbZnAs10.255*0.189*0.160*0.0640.317*0.290*0.247*Cd0.255*10.352*0.397*0.265*0.329*0.660*0.431*Cr0.189*0.352*10.532*0.103*0.716*0.383*0.424*Cu0.160*0.397*0.532*10.417*0.495*0.520*0.387*Hg0.0640.265*0.103*0.417*10.103*0.298*0.196*Ni0.317*0.329*0.716*0.495*0.103*10.307*0.436*Pb0.290*0.660*0.383*0.520*0.298*0.307*10.494*Zn0.247*0.431*0.424*0.387*0.196*0.436*0.494*1注：*P0.05水平下，显著相关；*P0.01水平上，极显著相关（2）主成分分析：对城市土壤中这8种重金属As，Cd,Cr,Cu,Hg,Ni,Pb和Zn进行主成分分析，当把因子固定为6个时，通过spss进行主成分分析时，这6个因子所对应的累积贡献率为93.156%（见表（7），因而可以判定这6个因子都为它的主成分因子，从主成分分析的旋转成分矩阵（见表（8）上可以看出 表（7）6个因子所对应的累积贡献率解释的总方差成份初始特征值提取平方和载入旋转平方和载入合计方差的 %累积 %合计方差的 %累积 %合计方差的 %累积 %13.56044.50044.5003.56044.50044.5001.82722.83422.83421.15014.37758.8771.15014.37758.8771.52519.06341.8973.96512.06370.941.96512.06370.9411.04513.06854.9654.7689.59680.537.7689.59680.5371.03212.89867.8635.5787.22087.756.5787.22087.7561.02012.74480.6076.4325.39993.156.4325.39993.1561.00412.54993.156表（8）6个主成分因子的得分系数矩阵旋转成份矩阵a成份123456As0.1350.1300.0160.9740.0430.084Cd0.2230.9180.1480.0860.0310.126Cr0.8590.187-0.0160.0040.2450.141Cu0.3950.1770.2730.0310.8100.110Hg0.0170.1370.9670.0160.1730.069Ni0.8900.1010.0510.1960.1030.170Pb0.0730.7180.0610.1620.4960.275Zn0.2610.2350.0820.0970.1220.916提取方法 :主成分分析法。 旋转法 :具有 Kaiser 标准化的正交旋转法。 第一主成分中主要是Ni和Cr，第二主成分中主要是Cd和Pb，第三主成分中Hg占主导作用；第四主成分中As占主导作用；第五主成分中Cu占主导作用，第六主成分中Zn占主导作用。结合相关性分析和主成分分析结果我们可以得出：在主成分因子1中，Ni，Gr的因子负荷量较高，而且他们之间存在显著相关性，这说明它们污染来源可能相同，可能主要来源于工业污染。在主成分因子2中，Cd和Pb的因子负荷量较高，可能来源于矿物的开采含铅汽油的排放。在主成分因子3中，Hg的因子负荷量较高，占主导作用，除来源于土壤母质外，还包括人为因素诸如：工业生产废料，城市生活垃圾堆放，农业施肥以及污灌等。在主成分因子4中，As的因子负荷量最高，占主导作用，除来源于地质岩石外，还有农业化肥，矿山开采，矿石冶炼等人为因素影响，且可能由于点污染。在主成分因子5中，Cu的因子负荷量较高，占主导作用，可能由于铜矿开发产生的尾石，矿石，三废的排放，有机肥，含铜农业化学物质等，含铅汽油的排放。在主成分因子6中，Zn占主导作用，与其它几种重金属含量关系较小可能与其独特的积累特征如点源污染等有关。综上所述，可以推断出该市城市土壤重金属污染与工业和城市生活垃圾污染物的积累，农业，城市交通有关。（三）聚类分析：然后采用R型聚类分析对八种元素进行相关性处理：使用的方法是最短距离法 它的直观意义为两个类中最近两点间的距离，R型聚类树形图图见下： 图（11）重金属元素聚类树形图注：树状图的的横坐标1，2, 3，4,5,6,7，8分别代表As,Cd ,Cr,Cu,Hg,Ni,Pb,Zn这8种元素。得出的结论与主成分分析类似，因此我们这样总结：该市城市土壤重金属污染与工业和城市生活垃圾污染物的积累，农业，城市交通有关。5.3问题（3）模型的建立与求解我们建立了两个模型5.3.1模型一的建立与求解由分析根据问题（1），（2）的结果可初步判断：Cd有三个污染源，Pb有两个污染源，其中两个污染源是共同的；Cr和Ni是同一污染源；As有两个污染源；Cu有一个污染源；Hg有三个污染源；Zn有两个污染源。设为某元素第i个污染源，设为第i个污染源的面积，有第一问可知内插值点的个数，内第j点的坐标记为，元素浓度记为，以内的插值点作为样本点，对于同一污染源的N种元素其内样本点的坐标取值相同，绘制如下层次分析流程图：目标层准则层决策层 加权平均浓度污染源污染源污染源样本点样本点点样本点样本点样本点样本点 因此我们建立了非线性规划模型：目标函数： (6)约束条件： (7)又上可得到最优权值和，为可能污染源对加权平均浓度的影响力，越小，污染源对加权平均浓度的影响力越大，我们可以设定一个阀值，若0.1，则可判断是污染源；表示第i个可能污染源中第j个样本点在该污染源中所占的比重，显然，越大，该点越接近于污染源中心。利用Lingo软件对模型进行求解（程序见附录六）：对As进行模拟，得到，可以判断显然为污染源，是污染源的可能性很小，区域浓度差异可能是As的沉积引起的。观察所有，可知越大，其对应点的As浓度越大，也可理解为越接近传染源的中心，显然，这是与As的空间分布相符的。同理可得到其它各重金属元素的污染源的位置。5.3.2模型二的建立与求解（1）.区域网格搜索法：假设在指定的区域内只由一个污染源头造成的，然后对于指定的区域内在坐标下求出金属浓度的最大值，作为最大污染指标（），然后根据坐标做成合适的方阵，将此方阵作为指标方阵。接着将城区插值形成若干个与指标方阵大小相同的方阵，把金属浓度作为Z轴。将指标方阵中的污染指数与城区的污染指数做方差求和，式子如下：在对城区中的每个方阵区域的进行比较，找出最小的即为金属浓度与指标最为接近的。因为指标是污染最为严重的，所以最小的的区域认为是污染源头。 第二问中通过聚类分析已经将是Ni和Cr分成一类， Cd和Pb分成一类， Hg，As，Cu，Zn分别为一类。基于MATLAB求出金属经过分类的和所对应的区域：金属Ni和CrCd和PbHg AsCuZn1.90E+063.51E+052.17E+098.05E+031.67E+083.12E+08对应区域(上下代表X，Y轴)（3000,4000）（5000,6000）（11000,12000）（4000,5000）（2000,3000）（2000,3000）（3000,4000）（8000,9000）（3000,4000）（2000,3000）（2000,3000）（2000,3000）拿出比较有代表的As进行研究，根据第一问已经得出有金属浓度的因素决定颜色的深浅，在第一问中的三维图进行变换视角，和As得出数据进行验证：从图像中可以看出在此处的得出的数据符合第一问中给出的值，所以模型正确。5.4问题（4）模型的分析，建立与求解5.4.1 首先对所建立的模型进行分析。模型的优点：运用多种插值方法进行比较，选择较优的插值方法，减小了由于方法不当引起的系统误差。能比较直观地分析土壤，重金属的浓度变化。运用Kriging插值法，误差小，与实际吻合度高，且外推能力强。比较同一重金属元素在不同区域的浓度分布，能直观地反映某一重金属元素在不同区域的污染程度。用内梅罗综合污染指标法综合地反应土壤的污染程度，解决了单因子污染指数法的不全面性。模型的缺点：相邻采样点距离较大，运用Kriging插值可能存在误差。运用Kriging插值时未考虑海拔影响。只考虑了重金属在土壤中的传播，忽略了河流，大气，其它元素对地质环境的影响。5.4.2需要收集的信息：将城区划分为间距为1公里左右的网格子区域，按照每平方公里随机均匀的采取5个样本点，再对每个样本点的表层土（0-50cm每隔10cm取一个样品）进行取样，编号，并用GPS记录采样点位置，记录所含元素浓度数据。对该城市各城区主要河流进行采样，编号，并记录位置和元素含量数据。对该城市各个功能区的空气进行各元素的浓度测量，记录采样点位置，元素含量值，并编号。对该城市各区的植物进行采样，编号，并用GPS记录采样点位置，并记录植物中元素含量值。由于城市土壤的岩石中所含元素是经年累月的结果，需考虑在采样期间是否进行过大面积的土地翻新对同一采样点利用同一方法定期对上面5个采样种类进行采样，记录浓度，时间。5.4.3模型的建立：1.运用Kriging插值法分别作出每一种元素在地表，水域，大气的空间分布图，由图像得出各元素在地表，水域，大气中的浓度分布情况，得出各统计量。2.对该城市城区采集的植物中各元素浓度进行Kriging插值，并绘制其空间分布图3.对该市空气中元素浓度进行Kriging插值，并绘制空间分布图，得出统计量。4.绘制个各元素在土壤中的垂直方向上的变化图，分析各元素在各土层的浓度变化。5.每次定期取样均进行以上4个过程。以定期取样所得统计量的均值作为原始数据，假设做了N次定期检测，第j种元素第i次的浓度均值为。（1） 关联度分析设影响地质环境演变的元素有m种，将每次检测的各元素浓度均值与背景值P（其余元素的背景值的确定方法与题給雷同）的比值作为参考数列其中k表示时刻。假设有m个比较数列则比较数列对参考数列在k时刻的关联系数为： （12）其中为分辨系数比较数列对参考数列关联度为： （13）由上可得出m种元素的关联度，关联度越大，说明在污染物中所占的比重就越大。（二）灰色预测（1） 数据检验与处理对每一种元素进行灰色预测，设参考数据为（表示第k次检测土壤中某元素浓度的平均值与背景值的比值），计算数列的级比： (14)如果所有的级比都落在可容覆盖内，则数列可以作为模型GM（1,1）的数据进行灰色预测。否则，需要对数列作必要的变换处理，使其落入可容覆盖内。即取适当的常数c，作平移变换： (15)则使数列的级比（2） 建立模型建立模型GM（1,1），对已知参考数列做一次累加生成数列： (16)记取的加权均值，则， 为确定的参数，记 的白化微分方程模型为： (17)其中a是发展灰度，b是内生控制灰度。由于，取为灰导数，为背景值，则方程（17）相应的灰微分方程为或 (18)即矩阵形式为，其中用最小二乘法球的参数估计值为 于是方程（18）有响应（特解）4则可以得到预测值 (19)（3） 检验预测值残差检验：令残差为，计算 (20)如果0.2,则可认为达到一般要求；如果0.1,则认为达到较高要求。级比偏差值检验首先由参考数据，计算出级比，再用发展系数求出相应的级比偏差 (21)如果0.2,则可认为达到一般要求；如果0.1，则认为达到较高要求。（3） 预测预报由GM（1,1）所得到的指定时区内的预测值，实际问题的需要，给出相应的预测预报。参考文献1秦耀东，土壤空间变异研究中的半方差问题/view/cf87a4ceda38376baf1faeef.html，2011/9/102/view/c9a7e727a5e9856a5612609c.html3/view/02201b8a680203d8ce2f24cf.html4孙洪泉，地质统计学及其应用，徐州；中国矿业大学出版社，1990.91935雷志栋，杨诗秀，田间土壤水分渗透的空间分布，水力学报，1987（3）；196刘会灯，朱飞，MATLAB编程基础与典型应用，北京：人民邮电出版社，2008.77姜启源，谢金星，叶俊，数学模型，北京：高等教育出版社，2003.8附录附录一：load F:/2.txtload F:/1.txtnum=X2;x1=num(:,1);x2=num(:,2);x3=num(:,3);x4=num(:,4);x5=num(:,5);x6=num(:,6);x7=num(:,7);x8=num(:,8);y1=X1(:,1);y2=X1(:,2);y3=X1(:,3);y4=X1(:,4);z1=0:300:29000;z2=0:300:19000;z1,z2=meshgrid(z1,z2);z3=griddata(y1,y2,y3,z1,z2,v4);z4=griddata(y1,y2,x1,z1,z2,v4);mesh(z1,z2,z3,z4);z4=griddata(y1,y2,x2,z1,z2,v4);mesh(z1,z2,z3,z4); z4=griddata(y1,y2,x3,z1,z2,v4);mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x4,z1,z2,v4);mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x5,z1,z2,v4);mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x6,z1,z2,v4);mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x7,z1,z2,v4);mesh(z1,z2,z3,z4)z4=griddata(y1,y2,x8,z1,z2,v4);mesh(z1,z2,z3,z4)1.7417932.2304722.22643.7426652.6583051.4912422.2292383.4349082.0142753.0239321.7228769.66180618.353011.6107053.0013174.0279351.1233591.171691.2567641.3119151.1701731.2564061.1792231.0622351.5855682.7693421.8727074.71325212.766361.4322852.049493.5196361.7399212.1580221.4076132.2872513.2854781.2430661.9583412.235395附录二：x1=5.6765 0.3023 53.509 55.106 0.29971 17.2618 61.740 201.20;x2=3 0.2246 69.9 162.7 1.627 9.9 50 338.7;x3=14.67650.9761263.2090543.2060 5.1807100000000046.9618211.74001217.3000x4=3.3235 0.3715 156.1910 432.9940 4.5813 12.4382 88.2600 814.9000load F:/2.txtnum=X2;for i=1:319if num(i,1)x4(1)num(i,1)=num(i,1);else if num(i,1)2*x1(1) num(i,1)=0;else if num(i,1)(1/3)*x1(1) num(i,1)=0; else if num(i,1) (1/3)*x1(1) num(i,1)=x1(1)-3*x3(1); else if num(i,1)x3(1)&num(i,1)2*x1(1) num(i,1)=x1(1)+3*x3(1); end end end endendend附录三： Cd元素空间分布图 Cr元素空间分布图 Cu元素空间分布图 Hg元素空间分布图 Ni元素空间分布图 Pb元素空间分布图 Zn元素空间分布图附录四：num=xlsread(F:/1.xls);for i=1:319num(i,6)=num(i,6)./1000;num(i,9)=num(i,9)./1000;endx=num(1:44,5:12);x1=num(45:80,5:12);x2=num(81:146,5:12);x3=num(147:284,5:12);x4=num(285:319,5:12);x5=3.6000 0.130 31.0000 13.2000 0.035 12.3000 31.0000 69.0000;for i=1:8t1(i)=mean(x(:,i)/x5(i);t2(i)=mean(x1(:,i)/x5(i);t3(i)=mean(x2(:,i)/x5(i);t4(i)=mean(x3(:,i)/x5(i);t5(i)=mean(x4(:,i)/x5(i);endw1=sqrt(max(t1(i)2+mean(t1(i)2)/2);w2=sqrt(max(t2(i)2+mean(t2(i)2)/2);w3=sqrt(max(t3(i)2+mean(t3(i)2)/2);w4=sqrt(max(t4(i)2+mean(t4(i)2)/2);w5=sqrt(max(t5(i)2+mean(t5(i)2)/2);附录五： As元素在各区的分布 Cd元素在各区的分布 Cr元素在各区的分布 Cu元素在各区的分布 Hg元素在各区的分布 Ni元素在各区的分布 Pb元素在各区的分布 Zn元素在各区的分布附录六：model:sets:x/1.10/;y/1.10/;com(x,y):w1,w2,z1,z2,z3;endsetsmin=sum(com(i,j):(z0-z3(i,j)2);z0=w11*sum(com(i,j):w1(i,j)*z1(i,j)+w22*sum(com(i,j):w2(i,j)*z2(i,j);w11+w22=1;sum(com(i,j):w1(i,j)=1;sum(com(i,j):w2(i,j)=1;data:z1=0.2269920.6174920.9915621.3441321.6693621.9603922.2089922.4052422.5370522.5897220.4393320.8006221.1417421.4570721.7400821.9831622.17733