2011高教社杯全国大学生数学建模竞赛(同名3885)

高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）：A 我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：同济大学参赛队员(打印并签名)：1.冯建设2.赵云波3.刘雄飞指导教师或指导教师组负责人(打印并签名)：日期：2011年9月11日赛区评阅编号（由赛区组委会评阅前进行编号）：2011高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）一问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求通过数学建模来完成以下任务：(1)给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2)通过数据分析，说明重金属污染的主要原因。(3)分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4)分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、问题分析本题中，城市表层土壤金属污染分析需要综合不同区域，不同金属的综合影响，根据随机的数据采样点，考虑通过统计与插值的分析方法进行处理。对于问题一，需要给出8种主要重金属的空间分布，并分析城区不同区域的污染程度。由于采样数据给出的采样点相对来说空间位置比较散乱，首先可以考虑对此进行网格化的数据处理，然后通过Kriging方法进行空间散乱点的插值处理。Kriging是一种距离加权的插值方法，从地质统计学中借鉴而来,在空间数据场满足给定的统计分布特征的前提假设下进行插值。通过对函数的插值结果的观察，与原始采样数据的空间分布相比较，进一步的确定该方法可以获得的插值效果是否有效。对城区内不同区域重金属污染程度的综合评价，要对各重金属评价指标分别加权。权重是衡量因子集中某一因子对土壤污染程度影响相对大小的量,权重系数越大,则该因子对土壤的影响程度越大，考虑利用熵权法来确定各重金属评价指标的权重系数，在信息论中信息熵表示系统的有序程度，一个系统的有序程度越高，则信息熵越大，反之，一个系统的无序程度越高，则信息熵越小。所以，可以根据各项指标的指标值的差异程度，利用信息熵的这个工具计算出各指标的权重。然后由综合权重进行线性加和，可以得到各个区域的综合评定指标，同时根据金属含量背景值进行等级标准的划分，从而确定不同区域的污染程度。对于问题二，可以结合问题一中得到的模型，同时对采样数据进行简单的分析，根据各区域的污染程度的不同，可以感性的得到重金属污染的主要原因必定和污染最严重的区域有直接关系。对于问题三，需要建立重金属污染的传播特征模型，并确定污染源的位置。可以考虑根据扩散定律建立模型，首先借鉴一维扩散方程建立模型，但是这样的话在数据处理上必须首先根据采样点浓度特征大致确定污染源的位置，然后建立方程，根据采样数据提取信息求解，由于采样点较多，本问题的处理可能会遗漏部分有用信息。如果与实际偏差较大，可以进一步寻求其他解决方法。对于问题四，可以基于问题一二三得到的模型进行综合考虑，从数据的处理方式，以及地质演变过程中时间变量的影响等因素入手，寻求更好的建立模型的方式。三模型假设金属污染浓度场按稳定场处理，即各坐标点浓度不随时间变化。假设各区域的土壤特性相同。重力对金属污染的影响忽略。四变量与符号说明变量符号符号说明Xi取样点横坐标yi取样点纵坐标y(x点x处半差函数yyλ各点对应权重系数C估计值点与i点之间的变差函数值Ci点与j点之间的变差函数值X某金属浓度值矩阵X某金属浓度值归一化矩阵P第i个因素下第j个评价值的比重第个因素的熵值W各金属权重矩阵dj第j区的加权综合浓度指标五模型建立与求解模型需要给出8种主要重金属的空间分布，并分析城区不同区域的污染程度。由于采样数据给出的采样点相对来说空间位置比较散乱，首先可以考虑对此进行网格化的数据处理，然后通过Kriging方法进行空间散乱点的插值处理。Kriging是一种距离加权的插值方法，从地质统计学中借鉴而来,在空间数据场满足给定的统计分布特征的前提假设下进行插值。利用Kriging方法对各浓度散乱点进行插值处理，建立起浓度分布的空间曲面（曲面的高度值即用来表征重金属浓度值），具体介绍如下：Kriging方法就是对空间数据进行加权插值的权值设计方法。Kriging方法通过引进以距离为自变量的变差函数来计算权值。由于变差函数既可以反映变量的空间结构特性，又可以反应变量的随机分布特性，所以利用Kriging方法进行空间数据插值往往可以取得理想的效果。另外，通过设计变差函数，Kriging方法很容易实现局部加权插值，这样就克服了一般距离加权插值方法插值结果的不稳定性。1随机场和区域化变量首先，重金属的浓度数据场可表示为分布于空间的单值函数S=f(x，y，z)，由题意知S为标量，则数据场为标量场。运用统计学的方法来研究该数据场，首先将f看成随机函数，记为Z，依赖于多个自变量的随机函数，称为随机场。以空间点x的直角坐标为自变量的随机场称为一个区域化变量。区域化变量在观测前，可以看作是随机场；观测后就得到随机场的一个实现(一般都记作Z(x)，写法上不加区别)。浓度区域化变量同时反映地质变量的结构性和随机性特征。从地质学的观点来看，区域化变量可反映地质变量的以下特征：局部性、连续性、异向性、可迁性。2变差函数的构造假设空间点x只在一维x轴上变化，我们把区域化变量Z(x)在x和x+h(h为与x具有相同维数的距离向量)2个点处的值之差的方差之半定义为Z(x)在x方向上的变差函数，记为y(x，h)，即y(x进一步的，由于点x和h是在二维(或三维)空间中变化的，所以要考虑二维(或三维)变差函数。但它们都是以一维变差函数为基础的，只不过要考虑各向同性或各向异性，还要考虑结构的套合。这里暂时先以各向同性进行数据分析处理。3平稳性假设和本征假设当区域化变量Z(x)满足下列条件时，则称Z(x)满足二阶平稳(或弱平稳)。(1)在整个研究区域内，区域化变量Z(x)的数学期望存在，且等于常数，即E[Z(x)]=m(常数)，Px；(2)在整个研究区域内，区域化变量Z(x)的协方差函数存在且相同(即只依赖于滞后h，而与x无关)，即Cov{Z(x)在实际工作中经常连二阶平稳假设也不能满足，故提出本征假设。当区域化变量Z(x)的增量[Z(x)-Z(x+h)]满足下列条件时，则称Z(x)满足本征假设，或说Z(x)是本征的。(1)在整个研究区域内有，区域化变量Z(x)的数学期望存在，且等于常数，即E[Z(x)-Z(x+h)]=0，Px，Ph；(2)在整个研究区域内，增量[Z(x)-Z(x+h)]的方差函数存在且平稳(不依赖于x)，即Var[Z(x)-Z(x+h)]==2y4实验变差函数实验变差函数就是根据观测数据构造变差函数y(h)的估计值y*h。有了二阶平稳假设或本征假设，重金属浓度区域化变量Z(x)的增量[Z(x)-Z(x+h)]只依赖于分隔它们的h，而不依赖于具体位置x。这样，被向量h分割的每一个数据对{Z(xi),Z(xi+h)}（i=1，2，⋯，N(h)）可以看成是{Z(x),Z(x+h)}一次不同的实现（N(h)是被向量h相隔的数据对的个数）。这样便可以用求所有这些观测值的算术平均的方法来计算yy这就是实验变差函数的基本计算公式。5变差函数的理论模型为了对区域化变量的未知值作出估计，需要将实验变差函数拟合成相应的理论变差函模型，这些模型将直接参与Kriging方法计算。变差函数模型可以分为有基台值和无基台值两大类，这里利用常用的有基台值模型进行计算。3种常用的有基台值模型如下：(1)球状模型(亦称马特隆模型，在原点处为线性型)。球状模型的几何解释，是因为它起源于2个半径为a且球心距为2h的球体重叠部分的体积计算公式。它的一般的表示为y其中，C0为块金常数，C0+C为基台值，C为拱高，a当C0=0，C=1时，称为标准球状模型。其中，C0为块金常数，C0+C为基台值，C为拱高(2)指数函数模型(在原点处为线性型)。它的一般公式为y(h)=此处a不是变程。因为当h=3a时，有1-e-3≈0.95≈1。所以y(h)≈C0+C，故其变程为3a。当C0=0，C=1(3)高斯模型(在原点处为抛物线型)。它的一般公式为y(h)=此处a亦不是变程。因为当h=3a时，有1-e-3≈0.95≈1，所以y(h)≈C0+C，故其变程为3a。当C0=0，C=16Kriging方法插值设Z(x)是点x承载的区域化变量，且是二阶平稳(或本征)的。Zi(i=1，2，⋯，n)是一组离散的信息样品数据，即重金属浓度，它们是定义在点承载xi(i=1，2，⋯，n)上的。现要对点x0承载处的区域化变量进行估计，所用的估计量为Z0*=i=1nλiZi，它是n个数值的加权线性组合。Kriging方法的原则(1)无偏性条件若要使Z0*为Z0的无偏估计量，即要求E[Z0*-(2)Kriging方程组区域化变量在满足二阶平稳的条件下推导，可以得到估计方差的计算公式σ估计方差σE2对λ∂σ在无偏性条件下，为了使估计方差最小，这是个求条件极值的问题，要用到拉格朗日乘子法。令F=这里F是n个权系数λi和μ的(n+1)元函数，求出F对λi(i=1，2，⋯，n)和μ的偏导数，并令其为零，便得到下列∂F整理得j=1如果区域化变量只满足本征假设，而不满足二阶平稳假设，则利用协方差函数和变差函数的关系C(h)=C(0)-y(h)，可得用变差函数表示的Kriging方程组j=1这里，yi，j=y(xi，xj)=y(xi-xj)。具体算法实现过程通过数学软件matlab编程实现，（matlab代码见附件）。(3)kriging模型的具体实现随机取100个观测点作为原始数据，其分布如下图所示：值得指出的是，在初步的数据处理中，可以发现该城区的地形并不是规则的长方形，但是为了方便数据处理，这里在一个30*20的区域内随机取点，刚好可以完全涵盖该城区，但是这并不影响最终的结果，因为可以再得到模型之后再将不属于城区的区域挖去，这在matlab软件中是容易实现的。根据实验变差函数结果，作数据拟合，选择适当的变差函数模型进行插值，确定变差函数计算公式。这里有3种常用变差函数模型:球状模型、指数模型和高斯模型。但是这3种模型有3个关键参数需要确定，变程a、拱高C和块金常数C0。本题中块金常数C0=0，拱高C和变程a由实验变差函数图人工确定。对于第一种重金属As，由实验变差函数100个观测点的观测值得到的变差函数，可以计算得到y*可以求得近似值：C=20和a=3。然后确定变差函数模型，选定变差函数模型后，变差函数的计算公式可以显式写出。这样Kriging方程组的系数矩阵与增广矩阵都已经确定。然后进行Kriging插值，也就是求解Kriging方程组，确定加权系数，然后进行线性加权，即可得到该模型下的Kriging插值结果。（具体实现见附件matlab代码）经画图比较，在本题中高斯变差函数模型的符合程度较好，建立高斯变差函数模型之后得到如下图所示的函数图像：1号金属As分布原始数据结果1号金属As分布kriging插值结果同样的方式，对第二种重金属Cd浓度做变差处理，值得指出的是，为方便比较，这里不再进行随机取100个观测点的工作，而是直接利用在计算第一种金属As变差函数时候的随机点数据，在此基础上得到变差函数观测数据y*由实验变差函数观测图可以近似求得C=126000和a=4。然后确定变差函数模型，选定变差函数模型后，变差函数的计算公式可以显式写出。这样Kriging方程组的系数矩阵与增广矩阵都已经确定。然后进行Kriging插值，也就是求解Kriging方程组，确定加权系数，然后进行线性加权，即可得到该模型下的Kriging插值结果。经画图比较，在本题中高斯变差函数模型的符合程度较好，建立高斯变差函数模型之后得到如下图所示的图像：2号金属分布原始数据2号金属kriging插值数据相应的，利用kriging方法也可以进一步得到其六种元素的空间分布模型，作图如下：7确定各重金属评价指标的权重系数在综合评定该城区内不同区域重金属的污染程度时，还要对各重金属评价指标分别加权。权重是衡量因子集中某一因子对土壤污染程度影响相对大小的量,权重系数越大,则该因子对土壤的影响程度越大，在这里利用熵权法来确定各重金属评价指标的权重系数，具体原理及操作如下：方法优点：客观赋权法。背景：设有m个待评方案，n项评价指标，形成原始指标数据矩阵（这里m=319，n=8），对于某项指标，指标值的值越大，则该指标在综合评价中所起的作用。在信息论中信息熵表示系统的有序程度，一个系统的有序程度越高，则信息熵越大，反之，一个系统的无序程度越高，则信息熵越小。所以，可以根据各项指标的指标值的差异程度，利用信息熵的这个工具计算出各指标的权重。数据处理(i=1,2,…,m;j=1,2,…,n,)，这里取n=8，m=319。由于参与评价的各项指标有越大越优型、越小越优型，故需对矩阵中的特征值进行归一化处理，方法如下：X据此得到归一化矩阵X’:(i=1,2,…,m;j=1,2,…,n,)这里以越大越优型进行计算求解。2.计算第个因素下第个评价值的比重3.计算第个因素的熵值4.计算第个因素的差异系数对于给定的越大，因素评价值的差异性越小，则因素在综合评价中所起的作用越小。定义差异系数，则当因素越大时，因素越重要。5.定义权数，则就是熵权法确定的权重。记As、Cd、Cr、Cu、Hg、Ni、Pb、Zn的权重矩阵为W=[w1,w2,w3,w4,w5,w6,w7,w8]’计算得到：W=[0.0233350.0475360.0651130.1831340.490280.0222410.0445820.123778]’得到各取样点加权综合浓度指标D=XW=[d1,d2,…,d319]’；进一步得到i类区加权综合浓度指标，其中，dj为所有属于i类区的d，ni为其采样点个数。6.得到评价结果得到五个区的加权综合浓度指标如下：区12345加权综合浓度指标值0.0021670.0059420.0009530.0040960.001787对各重金属浓度背景值及标准差做相同归一化处理，得到归一化之后的背景平均值及标准差，然后对8种金属归一化之后的背景平均值和标准差进行加权求和，求得和各点加权综合浓度指标di平均值对应指标标准差对应指标0.0007990.000195在此基础上得到对污染程度等级进行分类的值如下：设轻度污染标准=平均值对应指标+3*标准差对应指标；中度污染标准=平均值对应指标+8*标准差对应指标；重度污染标准=平均值对应指标+15*标准差对应指标；即：污染程度无轻度污染中度污染重度污染得分<0.0013850.001385~0.0023610.002361~0.003727>0.003727故可得到各区的污染等级如下表：区1生活区2工业区3山区4主干道5公园绿地污染等级轻度污染重度污染无污染重度污染轻度污染由题设条件可知，按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同，通过以上的计算，得出的结论为：工业区与主干道严重污染，生活区与公园绿地轻度污染，山区几乎不受污染，这与实际情况符合的较好。根据该金属元素空间分布模型，可以确定工业区污染较为严重，故可以初步确定重金属污染的主演原因在于工业区活动对环境的影响。得到污染源分布如下：8传播特征第一问中已经确定了污染物的浓度场，记为u(x,y)，本例中，按稳定场处理。由扩散定律，得到扩散方程为：其中，F(x,y)为污染源强度，K(x,y)为（x，y）处土壤传播污染物的特性，f(x,y,u)为耗散强度。建立离散模型：设污染源为点源F(xi,yj)只有在源(xi,yj)处等于污染物的排放强度，在源处，恒为零。假设各分区内的土壤特性相同。则在非源处，则各区的其中，kj为j类区的平均土壤特性，Ki为i采样点处的土壤特性。因为城区内土壤特性相同，故f(x,y,u)可以简写为f(u)。取浓度场上的节点，建立离散方程组，通过求解方程组，可以分别得到F(x,y)，f(u)的一组数值，进行拟合，得出解析式，将解析式带回扩散方程，即可得出扩散规律。六、模型评价与改进1模型评价1．为了给出8种主要重金属的空间分布，并分析城区不同区域的污染程度。首先对空间位置比较散乱的采样点此进行网格化的数据处理，然后通过Kriging方法进行空间散乱点的插值处理，在空间数据场满足给定的统计分布特征的前提假设下进行插值。在求得相关参数后建立显示的高斯模型，通过对函数的插值结果的观察，与原始采样数据的空间分布相比较，可以发现二者具有较好的吻合度，结果令人满意。2.对城区内不同区域重金属污染程度的综合评价，要对各重金属评价指标分别加权。利用熵权法来确定各重金属评价指标的权重系数，熵权法的实际意义在这里体现得尤为明显，根据熵权法得到的相关系数均为正值，这一点也验证了熵权法在寻找个金属污染物权重时的正确性，然后由综合权重进行线性加和，得到各个区域的综合评定指标，同时根据金属含量背景值进行等级标准的划分，从而确定不同区域的污染程度，结果与实际完全符合，说明熵权法的运用是正确的，从而便于找到重金属污染的主要原因。2改进方向1.在建立重金属污染的传播特征模型，先假设了污染源的位置，然后考虑根据扩散定律建立模型，根据一维扩散方程建立模型，但是这样的话在数据处理上必须首先根据采样点浓度特征大致确定污染源的位置，然后建立方程，根据采样数据提取信息求解，由于采样点较多，本问题的处理可能会遗漏部分有用信息。如果与实际偏差较大，没有得到有效结果，可以进一步寻求其他解决方法。2.对于问题四，可以基于问题一二三得到的模型进行综合考虑，从数据的处理方式，以及地质演变过程中时间变量的影响等因素入手，寻求更好的建立模型的方式。若引入时间变量，将扩散模型修改为：如此一来，还需要知道浓度的时间变化率，即不同时刻采样点的浓度数值。然后要考虑二维(或三维)变差函数。但它们都是以一维函数为基础的，只不过要考虑各向同性或各向异性，以及结构的套合。或者可以进一步的加强考虑海拔Z的影响，得到四元的方程模型，然后建立求解。参考文献[1]周晓云,朱心雄.散乱数据点三角剖分方法综述[J].工程图学学报,1993,(01).[2]王靖波,潘懋,张绪定.基于Kriging方法的空间散乱点插值[J].计算机辅助设计与图形学学报,1999,(06).[3]乔家君.改进的熵值法在河南省可持续发展能力评估中的应用[J].资源科学,2004,(01).附件附件一重金属空间分布模型的建立与求解的matlab实现注：为减少篇幅，此处代码仅仅实现了对第一种金属元素含量的空间分布的求解，其余个元素的求解方式类似可以得到。%坐标网格化%zuobiao为采样点坐标矩阵%xx将坐标单位化为公里x=zuobiao(:,1)';y=zuobiao(:,2)';z=zuobiao(:,3)';xx=x/1000;yy=y/1000;[X,Y]=meshgrid(0:0.5:29,0:0.5:19);Z=griddata(xx,yy,z,X,Y,'v4');surf(X,Y,Z);shadinginterp%金属浓度数据网格化jinshu_1=jinshu(:,1)';J_1=griddata(xx,yy,jinshu_1,X,Y,'v4');jinshu_2=jinshu(:,2)';J_2=griddata(xx,yy,jinshu_2,X,Y,'v4');jinshu_3=jinshu(:,3)';J_3=griddata(xx,yy,jinshu_3,X,Y,'v4');jinshu_4=jinshu(:,4)';J_4=griddata(xx,yy,jinshu_4,X,Y,'v4');jinshu_5=jinshu(:,5)';J_5=griddata(xx,yy,jinshu_5,X,Y,'v4');jinshu_6=jinshu(:,6)';J_6=griddata(xx,yy,jinshu_6,X,Y,'v4');jinshu_7=jinshu(:,7)';J_7=griddata(xx,yy,jinshu_7,X,Y,'v4');jinshu_8=jinshu(:,8)';J_8=griddata(xx,yy,jinshu_8,X,Y,'v4');%注：产生100个随机坐标的代码只在对第一号金属进行变差函数观测的时候运行，在对其%他金属进行变差函数观测求解的时候只需利用当前的100个观测点即可，该代码不再运行%产生100个随机坐标%s_suiji即为随机观测点得坐标%fori=1:100%x_suiji(i,1:2)=[ceil(rand*30),ceil(rand*20)];%end%求变差函数的原始数据%bianchahanshu为变差函数矩阵bianchahanshu=[0,0];forj=1:100fork=1:100bianchahanshu=[bianchahanshu;...[sqrt((x_suiji(j,1)-x_suiji(k,1)).^2+(x_suiji(j,2)-x_suiji(k,2)).^2),...(J_1(x_suiji(j,2),x_suiji(j,1))-J_1(x_suiji(k,2),x_suiji(k,1))).^2]];endend%对变差函数矩阵进行处理，得到最终的变差函数矩阵form=1:600flag=find(bianchahanshu(:,1)==bianchahanshu(m,1));bianchahanshu(m,:)=sum(bianchahanshu(flag(:,1),:))/sum(flag==flag);bianchahanshu(flag(2:end,1),:)=[];mendbianchahanshu(:,2)=bianchahanshu(:,2)/2;%做出变差函数观测点的图像scatter(bianchahanshu(:,1),bianchahanshu(:,2))%根据变差函数，求每一个估计点权系数的系数矩阵以及增广矩阵%C为系数矩阵%C_plus为增广矩阵C=zeros(100,100);forp=1:100forq=1:100h=sqrt((x_suiji(p,1)-x_suiji(q,1)).^2+(x_suiji(p,2)-x_suiji(q,2)).^2);C(p,q)=20*(1-exp(-h/3));endendC_plus=zeros(101,101);C_plus(1:100,1:100)=C;C_plus(101,1:100)=1;C_plus(1:100,101)=1;C_plus(101,101)=0;%求解每一个估计点的权重系数矩阵，得到空间模型%suij