2012数学建模国赛一等奖论文.doc

2011高教社杯全国大学生数学建模竞赛 承 诺 书 我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。 我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。 我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。 我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）： 2011高教社杯全国大学生数学建模竞赛 编 号 专 用 页 赛区评阅编号（由赛区组委会评阅前进行编号）： 赛区评阅记录（可供赛区评阅时使用）： 评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：城市表层土壤重金属污染分析摘要本题主要通过对重金属污染浓度的分析，画出浓度的空间分布图，建立相关的模型，得出该地区的重金属污染程度、污染原因以及污染源的所在位置。针对不同重金属具有不同的污染源，对不同区域重金属元素浓度值进行分析。运用克拉格差值法，画出各重金属的浓度等高线图。从水平和垂直方向考虑浓度分布，在水平面与海拔上对金属浓度分别分析，得出在三维空间内浓度的分布状况。考虑到区域内各种重金属元素污染程度，及不同区域内元素的综合污染状况不同，先运用地积累法，对各区域内不同元素分别考虑其污染程度，再对每个区域用内梅罗综合污染指数法得出该地区综合污染程度。根据给出的样本数据对各种重金属元素之间的相关性进行分析，得出相关性矩阵，对于相关性比较明显的重金属，我们可以认为其来源大致相同，即有相同的污染源。然后采用KMO法和Bartlett法对原始数据进行主成分分析适宜性检验，在此基础上进一步对样本数据进行主成分分析，确定影响比较显著的几种重金属元素，从而简化模型，方便分析重金属污染的主要原因。在分析重金属污染传播的特征时，可以结合大气扩散模型和水体对流扩散模型的偏微分方程建立重金属元素的对流扩散模型。为求具体污染源的位置，在对扩散模型简化的基础上，利用反演扩散模型方程来求污染源。反演方程可以用遗传算法来实现，把问题转化为建立一个非线性最优化问题。比较用差分方法反演扩散模型与遗传算法反演扩散模型的优缺点，得出差分方法是建立在把偏微分方程转化为差分方程进而求解各个点处的浓度随时间变化差异性确定污染源，针对样本数据分布不规律，需要进行插值处理，会给模型的精度带来较大影响。而遗传算法具有极好的全局优化性能,减少了陷入局部极值的风险, 特别对非线性优化问题具有很强的适应性，所以采用遗传算法反演扩散方程，具有较好的合理性。但由于未考虑模型里的对流项以及海拔的影响，导致其误差较大，如能提供该地区具体的工业排放区及交通要道等的具体密集情况，进一步减小污染源位置的范围和个数的初始值，可以提高遗传算法得出的污染源的具体位置的准确性。关键词：等高线图 相关性 主成分分析 遗传算法 反演 扩散方程 目录1、 问题重述.32、 问题分析.32.1、对问题一的分析32.2、对问题二的分析32.3、对问题三的分析42.4、对问题四的分析43、 模型假设.44、 符号说明.45、 模型的建立与求解.55.1、问题一的模型建立与求解5 5.1.1、重金属空间分布.5 5.1.2、重金属的污染程度的分析.85.2、问题二的模型建立和求解11 5.2.1、重金属相关性分析.11 5.2.2、土壤重金属主成分分析13 5.3、问题三的模型建立和求解15 5.3.1 污染物传播特征的分析155.3.2模型建立和求解155.3.3模型可行性分析185.4问题四的分析和求解195.4.1优缺点分析195.4.2模型优化196、参考文献 207、附录.211、问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（010 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求你们通过数学建模来完成以下任务：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？2、问题分析2.1 对问题一的分析题目要求统计出8重金属元素在城区的空间分布，需要考虑重金属元素在X，Y，Z三维空间的浓度分布特征。由附件给出的数据分析得，当X，Y的坐标确定之后，该点的海拔高度Z相应确定下来。所以在统计重金属的空间分布时，可以先考虑元素浓度和平面位置的关系，再通过海拔分布图分析元素分布与海拔之间的关系。对这八种重金属元素分别单独考虑，利用附件中各元素在城区的样本浓度，采用克立格插值法，应用SUPFER 8.0软件对各元素含量分布进行空间分析，可分别绘制出8种重金属元素的空间等高线分布图。利用附件给出了的各重金属元素在当地的背景值，可以运用地积累指数法来分析该城区不同区域8种重金属元素各自的污染程度。同时，也可以运用内梅罗指数法综合评判该城区各区域内的重金属污染程度。2.2 对问题二的分析 根据附件给出的数据，我们分区域考虑，确定每个区域里8种重金属元素相互之间的相关性，相关性较显著的重金属元素可以认为来自于同一个污染源，由此可以大致确定每个区域污染源的种类。因为每个区域中重金属的种类较多，相互之间又具有比较明显的相关性，所以我们可以在采用KMO法和Bartlett法对原始数据进行主成分分析适宜性检验的基础上，进一步考虑使用主成分分析法，由此分析出每个区域影响比较显著的若干种重金属，从而确定出重金属污染的主要原因。2.3 对问题三的分析重金属元素传播途径主要有水流传播、土壤渗透、大气扩散等。在城区我们主要考虑以大气扩散方式传播。由此，我们通过对流扩散方程建立模型，从中体现重金属传播的主要特性。在随时间的传播过程中，污染源所在地的浓度随时间的变化率较大。由此可以通过对模型的反演确定污染源的大概位置。2.4 对问题四的分析前三个问题的求解模型虽然可以比较合理地解决问题，但是也存在一定的局限性。由此我们可以进一步收集该城区土壤成分的信息，建立更优化的模型。3、模型假设1、固定时间内污染元素浓度保持稳定。2、元素污染浓度在水平面与海拔高度的分布成独立的相关性。3、污染元素的扩散只受水平面的扩散影响。4、污染元素在水平面的扩散在各方向的扩散速度相同。5、污染源的位置大致固定，不随时间变化。4、符号说明Ci： 重金属污染物i的实测浓度Bi： 重金属元素n的背景值Igeo：地积累指数PI： 重金属污染物i的污染指数Si： 重金属污染物i的国家标准值Pmax：某区域中重金属单项污染指数最大值Pave：某区域所有重金属单项污染指数的平均值IPI： 某区域土壤中重金属综合污染指数： x方向上的网格步长Ex ： x方向的扩散系数Ey ： x方向的扩散系数Ez ： x方向的扩散系数5、模型的建立和求解5.1问题一的模型建立和求解5.1.1 重金属空间分布的求解克里格插值法（Kriging）是用协方差函数和变异函数来确定高程变量随空间距离而变化的规律,以距离为自变量的变异函数,计算相邻高程值关系权值,在有限区域内对区域化变量进行无偏最优估计的一种方法，是地统计学的主要方法之一。根据采样点各重金属元素的含量分布，采用克立格插值法，应用SUPFER 8.0软件对海拔高度和8种重金属的含量分布进行空间分析，可绘出以下空间分布图。 图5-1-1 海拔分布图图5-1-2 As空间分布图由图5-1-2分析得，As元素浓度较高部分在左图呈明显黑色的三个区域中，说明这三个区域As元素污染较严重。同时考虑浓度与海拔的关系，右图呈现各浓度主要分布在海拔高度为0100m之间的范围内，且与海拔无特别关系。图5-1-3 Cd空间分布图由图5-1-3分析得，Cd元素浓度较高部分在左图中的几个颜色较深的区域，且该元素在城区内的污染都比较严重。由右图得在海拔高度为0100m的范围内均有浓度分布，且浓度较高部分分布在海拔较低区域内。图5-1-4 Cr空间分布图由图5-1-4分析得，Cr元素在城区内总体分布较少，但有一块很明显的黑色区域，表明污染主要集中在这一块。同时，右图的海拔图显示污染主要分布于0-100m海拔之间，且浓度量大致在0-100 (g/g)之间。 图5-1-5 Cu空间分布图 由图5-1-5分析得，Cu元素在城区内的总体分布较少，高浓度主要集中于两块明显黑色的区域，表明这些区域污染较严重。右图的海拔图显示污染主要分布于0-100m海拔之间，且与海拔无特别关系。图5-1-6 Hg空间分布由图5-1-6分析得，Hg元素在城区内的总体分布较少，高浓度主要集中于三块明显黑色的区域，表明这些区域污染较严重。右图的海拔图显示污染主要分布于0-100m海拔之间，且与海拔无特别关系。图5-1-7 Ni空间分布图由图5-1-7分析得，Ni元素在城区内的总体分布较少，高浓度主要集中于二块明显黑色的区域，表明这些区域污染较严重。右图的海拔图显示污染主要分布于0-100m海拔之间， 且与海拔无特殊关系。 图5-1-8 Pb空间分布图 由图5-1-8分析得，Pb元素浓度较高部分在左图中的两个颜色较深的区域，说明这两个区域As元素污染较严重。由右图得在海拔高度为0100m的范围内均有浓度分布，且浓度较高部分大致分布在海拔较低区域内。 图5-1-9 Zn空间分布图由图5-1-9分析得，Zn元素浓度较高部分在左图中的几个颜色较深的区域，说明这几个区域Zn元素污染较严重。由右图得在海拔高度为0100m的范围内均有浓度分布，且与海拔无特殊关系。综上所述，因为各种重金属有不同的来源，比如交通污染，工业污染，生活垃圾污染等，所以呈现出不同的区域分布。但从海拔方面来分析，其主要分布在0100m的低海拔区域。5.1.2、重金属的污染程度的分析 1）模型建立与求解 我们使用地积累指数法建立模型。Muller曾提出地积累指数定量评价沉积物中的重金属污染程度，并规定了相应的污染程度级别划分标准（表5-1-1），这种评价方法也可以用来评价土壤中的重金属的污染程度及其分级情况。其计算公式为： 式中Ci为重金属元素i在样本中的浓度，Bi为重金属元素i的背景值，K是为考虑各地岩石差异可能会引起背景值的变动而取的系数（一般取值K=1.5）。表5-1-1 Muller 地积累指数分级表根据附件表格给出的各重金属元素在多个采样点的浓度，同时我们也知道各个采样点对应的区域，所以我们可以用Matlab计算出每种重金属在每个区域的平均浓度Ci。Bi取各元素背景值的平均值代入计算。则可以算出各元素在各个区域的地积累指数Igeo。表5-1-2生活区重金属污染程度表5-1-3工业区重金属污染程度表5-1-4山区重金属污染程度表5-1-5主干道路区重金属污染程度表5-1-6公园绿地区重金属污染程度2）为了全面反映各污染物对土壤的不同作用，突出高浓度污染物对环境质量的影响，采用内梅罗综合污染指数法再次建立一个模型。内梅罗指数不仅考虑到各种污染物的平均污染状况，并且强调最严重的污染物对污染指数的影响。因此，可采用单因子污染指数( PI) 和内梅罗综合污染指数( IPI) 对城区土壤重金属的污染程度进行综合评价。单因子污染指数法：式中，PI 为重金属污染物i的污染指数;Ci为重金属污染物i的实测浓度，Si为重金属污染物i的国家标准值。内梅罗综合污染指数法:式中，Pmax为一个区域采样点中重金属单项污染指数最大值；Pave为一个区域所有重金属单项污染指数的平均值。IPI为一个区域土壤中重金属综合污染指数。评判标准如下：IPI1，表示无污染；1 IPI2，表示轻度污染；2 3，表示严重污染。由此可以进一步计算出内梅罗综合污染指数( IPI)。我们把城区分为5个区域分别计算。计算过程中Ci取一个区域中重金属i的采样值的平均浓度，Si取为重金属i的国家标准值。通过Matlab编程计算出的结果如下表5-1-7所示。表5-1-7各区域综合污染程度5.2、 问题二的模型建立和求解5.2.1、重金属元素相关性分析城市土壤重金属来源于成土母质和人类活动，同一来源的重金属之间存在着相关性，根据相关性可以判断土壤重金属污染来源是否相同。如果重金属之间存在显著的正相关，则其来源可能相同，否则来源可能不止一个。应用Matlab编程，可以计算出各个区域内8种重金属相互之间的相关性，并对其相关性做出分析。表5-2-1生活区重金属元素相关性在生活区采集的土壤样品中，As和Ni，Cd和Pb呈现比较明显的相关性，相关系数分别是0.605和0.802.，表明As和Ni可能具有相同的污染源，Cd和Pb可能具有相同的污染源。但其他元素之间的相关性都不明显，所以生活区的重金属来源相对复杂，很难于某种污染源进行确认。表5-2-2工业区重金属元素相关性在工业区采集的土壤样品中，Cd、Cr、Cu、Hg、Pb、Zn两两之间都有比较明显的相关性，表明这六种重金属元素可能具有相同的污染源，As和Ni具有明显的相关性，可能具有相同的污染源。表5-2-3山区重金属元素相关性在山区采集的土壤样品中，Cd、Pb、Zn两两之间呈现明显的相关性，所以此三者可能有相同的污染源。As和Cu、Cu和Hg之间有明显的相关性Cr和Ni、Ni和Zn也有比较明显的相关性。表5-2-4交通区重金属元素相关性在主干道路区采集的土壤样品中，Cd、Cr、Cu、Ni、Pb、Zn两两之间具有明显的相关性，所以可能有相同的污染源。As、Hg与其自身外七种重金属的相关性都不明显，两者与其他重金属的来源可能不同。表5-2-5公园区重金属元素相关性在公园绿地区采集的土壤样品中，As、Cr、Ni两两之间的相关性比较明显，可能有相同的污染源。Cu、Pb、Zn两两之间的相关性比较明显，可能有相同的污染源。5.2.2、土壤重金属主成分分析为了进一步明确该城区不同功能区土壤中重金属的污染来源, 利用Matlab软件对上述不同功能区土壤重金属的8项指标进行主成分分析。在利用主成分分析进行源解析前,采用KMO法和Bartlett法对原始数据集进行主成分分析适宜性检验。经计算, 其KMO值为0.73（大于0. 5的最小值要求），Bartlett检验统计的显著水平小于0. 01,因此可认为相关系数矩阵与单位阵有显著差异,以上检验结果均表明原始数据集适合进行主成分分析。数据处理及分析结果见下表。表5-2-6生活区主成分分析 表5-2-7工业区主成分分析生活区的8种重金属的全部信息可由3个主成分表示。即对前3个主成分进行分析已经能够反映全部数据的大部分信息。结合表中数据可以看出:三个主成分中第一主成分的贡献率为45.2%，特点表现为因子变量在元素Pb、Cd上有较高的载荷。研究结果表明汽车尾气颗粒物中铅含量较高,在怠速状态下向大气的排放量分别为01571 mg/min。又知Pb和Cd在该区相关性明显，所以大致可认为交通污染是生活区的一个主要污染原因。第二主成分的贡献率为14.2%，特点表现为因子变量在元素As、Zn上有较高的载荷，因为生活垃圾中会富含比较多的As和Zn，所以生活区生活垃圾的堆积可能是重金属污染的主要原因。第三主成分的贡献率为13.4%，特点表现为因子变量在元素Ni上有较高的载荷，普遍认为，Ni和Zn主要来源于轮胎的磨损，所以也表明交通污染可能是生活区的污染来源之一。工业区的8种重金属的全部信息可由2个主成分表示。第一主成分的贡献率为65.7%，特点表现为因子变量在元素Cr、Cu、Hg、Pb、Zn上有较高的载荷,且这些重金属在该区相关性显著，可以认为来源相同，工业污染排放是工业区重金属污染的主要原因。第二主成分的贡献率为15.8%，特点表现为因子变量在元素As、Hg上有较高的载荷, 这两种元素通常被视为燃煤和垃圾燃烧的标志元素，所以燃煤污染也是工业区重金属污染的一个主要来源。 表5-2-8山区主成分分析 表5-2-9主干道路区主成分分析山区的8种重金属的全部信息可以由3个主成分表示。第一主成分的贡献率为38.0%，特点表现为因子变量在元素Zn上有较高的载荷。第二主成分的贡献率为25.5%，特点表现为因子变量在元素Cd、Pb上有较高的载荷，第三主成分的贡献率为19.4%，特点表现为因子变量在元素Hg上有较高的载荷。综合考虑，因为山区的污染程度较轻，所以重金属污染的主要原因可能是山区富含较多的矿石。主干道路区的8种重金属的全部信息可以由3个主成分表示。第一主成分的贡献率为46.9%，特点表现为因子变量在元素Pb、Ni、Cu、Cd上有较高的载荷，所以交通污染中的汽车尾气排放和轮胎磨损可能是重金属污染的主要原因。第二主成分中Hg具有较高的载荷，所以也可能受到工业污染如仪表、电气工业等的影响。第三主成分中As具有较高的载荷，所以可能受到工业上排放As的工厂如化工、冶金业等的影响。表5-2-10公园绿地区主成分公园绿地区的8种重金属的全部信息可以由3个主成分表示。第一主成分的贡献率为48.9%，特点表现为因子变量在元素Cd、Cr、Pb、Zn上有较高的载荷，这些重金属元素主要来源可以认为是交通污染。 第三主成分中Hg具有较高的载荷，有研究发现由于大量的荧光灯、汞灯等气体放电灯的广泛使用,使得绿地表层土壤中汞的含量远远超过背景值。所以荧光灯等的大量使用可能是Hg污染的主要原因。5.3 问题三模型的建立和求解5.3.1 污染物传播特征的分析重金属元素传播途径主要有水流传播、土壤渗透、大气扩散等。在城区我们主要考虑以大气扩散方式传播。由此，我们通过扩散方程建立模型，从中体现重金属传播的主要特性。5.3.2模型建立和求解由问题分析得，我们主要采用对流扩散方程建立模型。上式中C为重金属浓度，Ex，Ey，Ez为在x，y，z三个方向的扩散系数，ux，uy，uz为在x，y，z三个方向的对流系数。模型一：基于污染源的（5-3-1）扩散模型，建立的逆求污染源模型的基本思想是通过扩散方程离散化，即把连续的定解区域用有限个离散点构成的网格来代替， 这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。针对本问题由于污染情况都大致在0100海拔高度以及忽略对流影响，以及考虑在x、y方向上的扩散系数相同，得到二维方程： 5-3-2把上式离散化得到方程： 由于污染源点处实际上具有浓度随时间变化相对较大的情况，问题转化为求浓度随时间的变化率最大问题，进一步转化为在等步长得时间内求浓度向方向的扩散率。这样模型可简化为只对相邻网格求浓度差，通过比较浓度差的值，来相应确定。对模型一的求解，由于本题给出的数据具有差异性，这里先利用Matlab对数据插值拟合为后续的x，y取网格化数据提供基础。这里插值对本题模型精度造成极大影响。由此我们提出了利用遗传算法反演扩散特征进而反求污染源的模型二。模型二， 在重金属传播特征的基础上引入污染源的确定模型如下： 5-3-3上式中C为重金属浓度，Ex，Ey，Ez为在x，y，z三个方向的扩散系数，ux，uy，uz为在x，y，z三个方向的对流系数。Mi为污染源的排放强度。 为污染源处对浓度变化率的影响。进而可以通过确定该项的值来确定污染源的大致位置。 记，0|r|(1/4) + (rand(2/4) + (rand(3/4) + (rand0) + 1; %5 %Te = (rand(1/3) + (rand(2/3) + (rand0) + 1; %4 %Te = (rand(1/2) + (rand0) + 1; %3 %Te = (rand0) + 1; %2 N(i,:)=(randperm(siz)=Te); endigBest = 1; %历史最佳个体索引NBest = N(1,:); %个体最佳值F = zeros(1,Num);data=xlsread(a.xls);V=data(1:siz,2);X=data(1:siz,11);X=(X+1000)./(max(X)+1000);%D = xlsread(C:Documents and SettingsAdministrator桌面数模2011最短距离D.xls);%CAS = xlsread(C:Documents and SettingsAdministrator桌面数模2011A区发案率.xls);%计算适应度for k = 1:Num F(k) = fun(N(k,:),V,X);endit = 1;while it Iter %选择两个个体交叉，变异 ind = randperm(Num); for k = 1:2:Num %交叉操作 n1 = N(ind(k),:); n2 = N(ind(k+1),:); n1 n2 = intercross(n1,n2,pinter); %变异操作 n1 = mutate(n1,pmut); n2 = mutate(n2,pmut); F1 = fun(n1,V,X); F2