1、2011 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛 承承 诺诺 书书 我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从 A/B/C/D 中选择一项填
2、写）：A 我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：西安电子科技大学 参赛队员(打印并签名)：1.欧阳照玮 2.李娟 3.王小磊 指导教师或指导教师组负责人 (打印并签名)：日期：2012 年 8 月 9 日 赛区评阅编号（由赛区组委会评阅前进行编号）：2011 高教社杯全国大学生数学建模竞赛高教社杯全国大学生数学建模竞赛 编编 号号 专专 用用 页页 赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评 阅 人 评 分 备 注 全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：1 城市表层土
3、壤重金属污染分析城市表层土壤重金属污染分析 摘摘 要要 本文围绕城市表层土壤重金属污染问题，对土壤重金属的污染程度、污染的主要原因、污染源的确定等问题分别建立了模型，并对结果进行了详细的分析。针对问题一，运用地统计学中的 Kriging 插值法对 8 种重金属元素的分布进行模拟。对每一种重金属元素，分别建立指数、球形和高斯变异模型，选取三种函数中方均根预测值最小的作为变异函数。运用 Matlab 绘制各个元素的三维空间分布图。随后计算各功能区原始样本点的土壤内梅罗综合污染平均指数，参照等级划分表确定各功能区的污染程度。结果显示，工业区污染最为严重，其次是交通区和生活区，公园绿地区和山区尚为清洁
4、，土壤受污染程度最轻。针对问题二，本文建立因子分析改进模型分析土壤中重金属的来源。运用SPSS 软件对 8 种重金属进行主成分分析，得到 3 个主因子。累积变异量达到了70.941%。将 8 种金属分为三类：Pb,Cu,Cr,Ni,Cd,Zn 的污染受到人类的工业生产的影响非常明显；Pb 和 Hg 的分布主要是交通工具尾气排放的影响；As 的分布和自然界中成土作用有关，受人类影响较小。针对问题三，首先对土壤中重金属元素建立一维对流弥散方程，求解得到原方程参数值及浓度分布公式。取不同的扩散时间，分析各地区重金属元素含量变化，得到重金属污染物的传播特性。最后，寻找每种重金属元素的非正常极大值区域，
5、并采用局部搜索的方法，利用实际数据拟合浓度分布公式，选取四个方向拟合优度均符合条件的点，作为该区域内的污染源。最终得到重金属的 23 个点污染源。针对问题四，加入对污染通过大气传播的途径的考虑，建立高斯扩散模型。与问题三的模型相结合，与实际情况相结合考虑重金属污染的传播扩散受风速、风向、降水、水流等因素的影响。最后，对模型进行了评价，分析了模型的优缺点，针对如何更准确找到污染源给出了自己的结论。关键词：关键词：Kriging 插值 变异函数 内梅罗综合污染指数 因子分析 对流弥散方程 高斯扩散 2 一、问题重述一、问题重述 随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影
6、响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为 1 类区、2 类区、5 类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1 公里左右的网格子区域，按照每平方公里 1 个采样点对表层土（0-10 厘米深度）进行取样、编号，并用 GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照
7、2 公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件 1 列出了采样点的位置、海拔高度及其所属功能区等信息，附件 2 列出了 8 种主要重金属元素在采样点处的浓度，附件 3 列出了 8 种主要重金属元素的背景值。要求你们通过数学建模来完成以下任务：(1)给出 8 种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2)通过数据分析，说明重金属污染的主要原因。(3)分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4)分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建
8、立模型解决问题？二、问题的分析二、问题的分析 近几十年，由于城市化和工业化的高速发展，土壤重金属污染问题日益引发人们关注。本文针对土壤重金属污染问题，主要解决四个相关问题，分别为确定重金属元素的空间分布并分析其在不同区域的污染程度，说明重金属污染的主要原因，确定污染源的位置，对模型进行评价与扩展并分析该城市地质环境的演变模式。2.1 2.1 问题一：确定重金属元素空间分布并分析不同区域重金属的污染程度问题一：确定重金属元素空间分布并分析不同区域重金属的污染程度 问题一要求确定 8 种重金属的空间分布并分析不同区域重金属的污染程度。由于受到人为活动和气候因素的影响，重金属元素在土壤中的分布满足一
9、定的随机性和结构性。传统统计学方法不考虑测定参数和测定位置之间的关系，只根据空间不同位置的测定数据，无法进行区域内的优化插值计算，更无法表现出重金属元素的空间分布，所以本文采用地统计学方法。Kriging 插值法是地统计学中一种重要的数据处理方法。可用于探测研究对象的空间自相关结构(或空间变异结构),并估计模拟变量值。通过 Matlab 绘制出不同重金属元素在土壤空间中的分布图。分析其空间分布特征及其变异规律。把重金属空间分布图与功能区分布图进行比较，对 5 个区域的八种重金属元素浓度数据进行处理，分析不同地区重金属污染的分布特点。求出五个区域八种重金属浓度的内梅罗综合污染指标，通过与国家环境
10、二级标准和该地区背景值相3 比较，以此来刻画重金属的污染程度。2.2 2.2 问题二：说明重金属污染的主要原因问题二：说明重金属污染的主要原因 问题二要求说明该城区重金属污染的主要原因，土壤中重金属含量可能受到人为或自然原因的影响，在遵循合理性、代表性、系统性以及可获得性的原则下，本文采用因子分析法来判别土壤重金属元素的来源。因子分析法是将具有相关性的多个原始指标的评价问题转换为较少的、新的综合指标的评价问题。所以，利用 SPSS 先分析 8 种重金属元素的相关性，以 8种重金属元素浓度为单个因子，利用因子分析法，得到影响土壤整体变异性的 3个重金属污染主因子。根据主因子的变异特征向量大小，将
11、 8 种重金属元素进行分类。由于各类因子所受影响不同，结合前面的 8 种重金属空间分布图，在一定程度上能反映出人类活动和自然过程的影响程度，从而得到重金属污染的主要原因。2.3 2.3 问题三：分析重金属传播特性并确定污染源位问题三：分析重金属传播特性并确定污染源位置置 问题三要求分析重金属污染物的传播特性，确定污染源的位置。根据土壤中重金属浓度距离污染源越远浓度越低的特性，引用一维对流弥散改进模型模拟重金属元素的传播模式，从而分析重金属污染物的传播特性。利用问题一中的插值数据建立二维矩阵，在 Matlab 中求出矩阵极值，再结合附件 3 中土壤重金属背景值，找出城区内不满足正常重金属元素含量
12、范围的极大值点，这样便大概确定了污染源所在的范围。2.4 2.4 问题四：模型评价并研究城市地质环境的演变模式问题四：模型评价并研究城市地质环境的演变模式 问题四要求确定收集信息，并建立模型使更好地研究城市地质环境的演变模式。研究地质环境的演变模式可简化为分析预测地质中重金属元素含量的变化。而在之前的模型建立过程中，我们没有考虑过地形、高度及沉降的影响，这样会对重金属的含量预测造成一定的误差。故我们考虑收集风向、风速等数据，对大气中的重金属元素建立高斯扩散模型，从而更加精确的分析预测土壤中重金属元素含量的变化，为研究地质环境的演变模式提供了一定的依据。三、模型假设与符号说明三、模型假设与符号说
13、明 3.1 3.1 条件假设条件假设 1.1.假设重金属元素在土壤和水中化学反应均匀。2.2.假设各区域成土母质中含重金属的浓度是相同的。3.3.假设各区域重金属分布稳定，污染源排放量不变。4.4.假设采样数据合理，能够反映当地的土壤质量。5.5.假设污染物浓度保持不变，且其他地点的污染物含量为叠加形式，即在较短时间内，污染物不会发生降解等使污染物浓度降低的情况。3.2 3.2 符号说明符号说明 符号符号 符号解释符号解释 ()r h 变异函数 0C 块金常数 1C 拱高 4 RMSPE 均方根预测误差 iP 重金属i的污染指数 iC 污染物i的实测含量 iS 污染物二级参考标准值 P综 内梅
14、罗综合污染指数 iY 重金属浓度标准化后的可观测指标 iX 公共因子 ija 因子载荷 k 污染物浓度随时间变化量 1Q 空间域内污染物通过的流量 2Q 空间域内污染物的增量 0Q 污染源释放的重金属总量 R 扩散影响半径 注：其他未注明的符号具体在文章中说明注：其他未注明的符号具体在文章中说明 四、模型准备四、模型准备 4.1 4.1 功能区分布图功能区分布图 利用 Matlab 软件，根据附件 1 的采样地理位置，绘出二维等高线图，并标出功能区分布情况，如图 1 所示：图图 1 该城区该城区 5 个功能区分布图个功能区分布图 在图 1 中，“*o+x”分别代表生活区、工业区、公园绿地区、交
15、通区5 和山区。从图 1 可以看出，右上角地势高，左下角地势低，生活区均匀分布于图形主对角线左下方，工业区位于副对角线上方且布局沿对副角线方向，公园绿地区称斜“V”字型，分布于图形偏左侧，交通区密集分布与整个图形中，山区则位于地势较高处。五、模型的建立与求解五、模型的建立与求解 5.1 5.1 问题一：问题一：确定重金属元素空间分布并分析不同区域重金属的污染程度确定重金属元素空间分布并分析不同区域重金属的污染程度 5.1.1 5.1.1 模型模型 I:I:土壤重金属的空间变异性模型土壤重金属的空间变异性模型 模型的建立模型的建立 （1 1）Kriging 空间插值法是以变异函数为基础的，其参数
16、设置和变异函数模型的选择对内插值效果影响很大。变异函数的计算一般要求数据符合正态分布，所以本文先对 8 种重金属浓度数据进行正太分布检验。检验数据的正态分布性有多种方法：频率分布直方图法、卡法检验法、Q-Q 图、P-P 图等。本文利用 SPSS 统计软件，对土壤重金属元素浓度数据取对数后绘制出 Q-Q 图，检验其正态分布性。8 中重金属元素浓度数据的 Q-Q 图，如图 2 所示：图图 2 八种重金属元素的正态八种重金属元素的正态 Q-Q 图（对数处理后）图（对数处理后）由上图可见，8 种重金属元素浓度的正态分布 Q-Q 图大致是一条直线，所以可看作服从正态分布。所以可以建立变异函数模型，进行
17、Kriging 差值。（2 2）利用 Kriging 插值法进行空间插值时，均要求变量具有符合本征假设的规定。即增量 ()()Z xZ xh（1）的方差函数存在且平稳(不依赖于x)。式中()Z x和()Z xh表示空间区域内距离为h的两个位置的观测值。即当空间距离h较小时,估计点与样点(已知高程点)的相关性较高、变异性较小；反之,估计点与样点的相关性较小、变异性较大。随着h的增加,变异函数()r h呈缓慢增加或不再增加,这时的()r h 称为临界变异6 值。理论上当0h 时,()0r h,但有时在原点附近出现不连续现象,称为块金效应（Nugget Effect)。可见,变异函数能同时描述区域化
18、变量的随机性和结构性,从而在数学上对区域化变量进行严格分析,是分析空间变异规律和空间结构的有效工具。Kriging 法提供了不同的理论变异函数模型,通过其结构及各项参数,从不同角度反映了空间变异性结构；常用的模型有高斯模型、球面模型和指数模型。高斯模型：高斯模型：空间相关性随距离增大而逐渐变小,当距离趋于无穷大时,相关性趋于 0。公式为:2201()1 exp(/)r hCCha（2）该模型变程为3a。如果核方差相对于与空间变化有关的随机变化很小时,可选用较弯曲的高斯曲线。球面模型：球面模型：空间相关性随距离增大而逐渐变小,距离大于一定值后,空间相关性为 0。公式为：301301003()()
19、022hhhr hCChaaaCCha（3）上式中0C为块金常数,表示随机变异的量。0CC为基台值,表示变量空间变异的结构性方差。C为拱高。a为变程，即曲线到达基台值时所对应的分离距离。当存在明显的变程和梁，核方差也很重要但数值不太大时，可用球面模型进行半方差拟合。指数模型：指数模型：空间相关性随距离增大呈指数衰减,当距离趋于无穷大时,相关性趋于 0。公式为:01()1 exp(/)r hCCh a（4）其变程为 3a。如果存在明显的核方差和梁，而没有渐变的变程，则可用指数模型拟合。模型的求解模型的求解 由前面分析可知，八种重金属元素的浓度均成对数正态分布状态，就适合用Kriging 法插值。
20、变异函数模型的选取要比较三个变异函数模型的优劣，采用均方根预测误差（RMSPE）作为变异函数模型效果判别指标。公式如下：2*112()()inikikkiRMSPEZxZ xn（5）其中，*()ikZx和()ikZ x分别表示检测点的Z检测值和Z估计值。在计算均方根预测误差时，我们从采样点中随机预留出一定数量的点作为检测点，而采用其他的采样点作为插值的源数据，通过 Kriging 插值计算各检测点7 的估计值。利用 Matlab 作图，其结果表示如下：图图 3 三种变异函数效果判别指标图三种变异函数效果判别指标图 图 3 中，横坐标表示三种变异函数，1、2、3 分别为采用指数模型、球面模型、高
21、斯模型。纵坐标表示三种变异函数下的RMSPE值。由上图可知，对于 8 种重金属的高斯模型的RMSPE明显高于前两种变异函数，所以，高斯模型不适合作为 Kriging 插值的变异函数。比较前两种模型可知，对于 Cd 元素，指数模型的RMSPE大于球形模型，说明用球形模型差值效果较指数模型好，同理可知其余 7 种重金属，用指数模型差值效果比球形好。在 Matlab 中，插值后得到各重金属的浓度空间分布图如下：8 图图 4 重金属浓度空间分布图重金属浓度空间分布图 从图 4 可以得出 Pb、Zn、Ni、Cr、Cu、Cd 在该地区的西南部的浓度较高，该地区地势平坦，是工业区集中的地方。Hg 的分布较为
22、特殊，在交通区附近含量较高。As 的分布则与上述几种重金属的分布有明显不同，有待进一步分析。5.1.2 5.1.2 模型模型 IIII：不同区域重金属的污染指数模型：不同区域重金属的污染指数模型 本小节利用单因子和内梅罗综合污染指数法，得到各个功能区的污染指数。其中内梅罗指数特别考虑了污染最严重的因子，在加权过程中避免了权系数中主观因素的影响，是目前仍然应用较多的一种环境污染指数法。1 1 模型的建立模型的建立 国家环境保护局和国际技术监督局联合发布的 土壤环境质量标准 中明确给出了三级土壤环境质量标准，其中第二级适用于保护农业生产，维护人体健康。故下面将 8 中元素的浓度与自然背景值和二级标
23、准进行对比，具体图形如下：生活区生活区 结论：Cd 和 Zn 元素的污染程度生活区接近于二级标准，As 和 Ni元素对该区域几乎没有污染。因为该区域地势较为平坦，经常会有生活垃圾堆积，其中含有较多的 Cd 和Zn，有一定程度的污染。0100200300400500600背景值 测量值 二级标准 9 工业区工业区 结论：对工业区污染很严重的元素有 Cd、Cu、Hg 和 Zn，其他元素也有 污染但 不严重。因为该区域存在铜矿业工厂、硫酸厂、油漆厂、采矿和化学工业厂等，由于风、河流等作用，使重金属元素不断向周围扩散。山区山区 结论：整体来说，各元素对山区的污染都很小，接近于背景值。由于山区地势较高，
24、重金属很难从污染源扩散到此。交通区交通区 结论：由于汽车尾气、车辆零件之间的摩擦等原因，使得交通区重金属 Cd 严重超标，Hg和Zn元素污染也较为严重，但低于二级标准。公园绿地区公园绿地区 结论：该区域 Cd 元素污染较为严重，可能是由农药、化肥和塑料薄膜 过度使 用所致。图图 5 8 种重金属背景值、测量值和种重金属背景值、测量值和二级标准值对比图二级标准值对比图 根据各个采样点 8 种重金属单因子污染指数值，分别求出 5 个功能区 8 种重0100200300400500600700背景值 测量值 二级标准 0100200300400500600背景值 测量值 二级标准 010020030
25、0400500600背景值 测量值 二级标准 0100200300400500600背景值 测量值 二级标准 10 金属元素污染指数均值，利用 Matlab 绘图，结果如图 6 所示：图图 6 五个功能区八种重金属污染指标对比图五个功能区八种重金属污染指标对比图 结论：图 6 反映了五个区域中工业区污染程度最严重，其次交通区，山区受污染程度最轻。土壤污染程度的确定采用单因子污染指数法和内梅罗综合污染指数法。单因子污染指数法，表达式为：iiiCPS（6）式中，iP为土壤中重金属污染物i，（1,2,.,8i）的环境污染指数；iC为污染物i的实测含量(/g g)；iS为污染物i的二级参考标准（/g
26、g)。（详见附录）内梅罗综合污染指数法，表达式为：22max=2iiiiC SC SavP综（/）（/）（7）式中，P综为该区域的综合污染指数；2maxiiC S（/）为土壤污染物中污染指数最大值；2iiC Sav（/）为土壤污染指数平均值。单因子污染指数和内梅罗综合污染指数的土壤质量分级标准见表 1：表表 1 基于污染指数的土壤质量分级基于污染指数的土壤质量分级 等级 污染指数 污染程度 污染水平 I P 0.7 安全 清洁 II 0.7 P1 警戒级 尚清洁 III 1 P2 轻污染 土壤作物均受到轻污染 IV 2 P3 中污染 土壤作物均受到中污染 V 3 P 重污染 土壤作物均受到重污
27、染 2 2 模型的求解模型的求解 利用公式（6），求出各个采样点的单因子污染指数。将各个采样点单因子污00.20.40.60.811.21.4AsCdCrCuHgNiPbZn生活区 工业区 山区 交通区 公园绿地区 11 染指数最大值与均值代入公式（7），求出各个检测点的八种重金属综合污染指标。与国家环境质量分级标准比较，得到 5 个区域的土壤质量，如表 2 所示：表表 2 5 个区域土壤综合污染指数个区域土壤综合污染指数 功能区 内梅罗综合指数 污染等级 生活区 0.7862 II 工业区 1.1274 III 山区 0.4028 I 交通区 0.9794 II 公园绿地区 0.7337 I
28、 结论：由表 2 2 可知，工业区受污染最严重，土壤污染程度最为严重；其次为交通区和生活区，污染程度已达到等级 II；公园绿地区和山区土壤质量最好，污染程度最轻。5.2 5.2 问题二：重金属污染的主要原因问题二：重金属污染的主要原因 因子分析法其主要原理是降维的思想，通过研究指标体系的内在结构关系，把多指标转化为几个相互独立而且包含原有指标大部分信息的综合指标的多元统计方法。新的综合指标称为主成分或公因子。利用各主成分的因子得分计算出每个评价对象的综合值，并以此作为综合评价的依据。5.2.1 5.2.1 模型模型 I:I:因子分析模型因子分析模型 1 1 模型的建立模型的建立 本文将 8 种
29、重金属元素浓度作为因子分析的 8 个变量因子，建立指标体系，即因子分析模型如下：111112211221122222881182288mmmmmmYXXXYXXXYXXX（8）其中，128,Y YY是已标准化的可观测的指标。12,mX XX为公共因子，i是特殊因子，它与公共因子之间彼此独立。ij是指标iY在公共因子jX上的载荷，因子载荷ij的统计含义是指标iY在公共因子jX上的相关系数，表示iY与jX线性相关程度。12,iiim说明了指标iY依赖于各个公共因子的程度。12,jjmj说明了公共因子jX与各个指标的联系程度。2 2 模型的求解模型的求解 在 SPSS 中导入原始矩阵，将原始数据指标
30、标准化，再考虑收集到的原有变量之间是否存在一定的线性关系。我们利用 SPSS 软件，借助变量的相关系数矩阵、卡方检验和 KMO 检验方法进行分分析。其软件结果如下表：12 表表 3 因子分析检验结果因子分析检验结果 KMO 和和 Bartlett 的检验的检验 取样足够度的 Kaiser-Meyer-Olkin 度量。.778 Bartlett 的球形度检验 近似卡方 905.711 df 28 Sig.000 由表3 3可知，对数据进行KMO和Bartlett检验，KMO的度量为0.778.Bartlett的球形度检验的近似2为 905.711.自由度为 28.双侧检验值为 0.000。说明
31、各个元素浓度之间存在较强的相关性，说明可以进行因子分析。用 SPSS 做出 8 种重金属浓度的因子分析碎石图，如图？所示：图图 7 因子分析碎石图因子分析碎石图 由 SPSS 生成的因子分析碎石图可以看出前 3 个因子的特征值明显较大，固选取 3 个主因子较为合适。这里在 SPSS 中将提取公共因子的标准值设置为 0.8，做方差最大正交旋转，已提取前 3 个主因子，得到如下结果：表表 4 特征值和累计贡献率特征值和累计贡献率 成份 初始特征值 合计 方差的%累积%1 3.560 44.500 44.500 2 1.150 14.377 58.877 3.965 12.063 70.941 4.
32、768 9.596 80.537 5.578 7.220 87.756 6.432 5.399 93.156 7.301 3.769 96.924 8.246 3.076 100.000 提取平方和载入 旋转平方和载入 合计 方差的%累积%合计 方差的%累积%3.560 44.500 44.500 2.228 27.849 27.849 1.150 14.377 58.877 1.776 22.204 50.053 13 .965 12.063 70.941 1.671 20.888 70.941 上表中第二列是因子变量的方差贡献，又称特征值，它是衡量因子重要程度的指标，例如第一行的特征值为
33、3.560，后面描述因子的特征值递减。第三列是各因子变量的方差贡献率，表示该因子描述的方差占原有变量总方差的比例。第四列是因子变量的累计方差贡献率，表示前m个因子描述的总方差占原有变量的总方差的比例。第五列和第七列则是从初始解中按照一定标准（在前面的分析中是设定了提取因子的标准是特征值大于 0.8）提取了 3 个主因子后对原变量总体的描述情况。各列数据的含义和前面第二列到第四列相同。由表中的数据可以看出在提取了三个主因子后累计的方差贡献率超过了 70%。三个主因子的方差贡献分别为 44.5%、14.38%、12.06%。已经较好的体现整体数据的特征。下表为极大正交旋转后的因子载荷矩阵表：表表 5 极大正交旋转后因子载荷矩阵表极大正交旋转后因子载荷矩阵表 成份 1 2 3 Pb.764.314.237 Cu.756.125-.365 Cr.735-.444-.303 Ni.723-.515-
