承 诺 书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 淮北师范大学 参赛队员 (打印并签名) ：1. 孙丰凯
指导教师或指导教师组负责人 (打印并签名) ： 李昌文
日期： 年 月 日
赛区评阅编号（由赛区组委会评阅前进行编号）：
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
城市表层土壤重金属污染分析及污染源位置的确定
摘要
土壤重金属污染问题是当今环境科学研究的重要内容, 特别是城市土壤重金属污染受到更多重视。城市人口集中、工业发达、交通拥挤等因素导致土壤中重金属不仅含量高而且种类较多。一旦城市土壤被重金属污染, 城市居民的健康将受到危害, 因此, 世界上许多国家, 尤其是欧美一些发达国家十分关注城市土壤重金属的污染状况。
针对问题一：模型一：我们根据数据用MATLAB 软件编程绘制出各功能区的空间分布图（用不同的颜色表示不同的区域）和各种重金属元素的浓度与空间分布的结合图，并将这两个图进行对比大致判断不同区域中各种重金出属在土壤表层的污染程度。我们又对模型一进行改进建立模型二更加准确的确定不同功能区的污染程度，对此我们用SPSS 软件计算出不同功能区的8种重金属元素的各类统计量，再采用MULLER 地积累指数法计算出地积累指数，根据其分级标准进一步确定其污染程度。
针对问题二：我们首先利用SPSS 软件对附件2的8种元素的浓度数据进行降维因子分析，得出各元素之间的相关矩阵，通过分析各个权重大致得出元素之间的关系，然后采用相关系数质心聚类法对原始数据进行聚类分析，当距离为15~20时主要分为3类，及重金属污染的主要原因有3种。
针对问题三：首先我们作海拔与各种重金属污染程度相关关系的显著性检验，排除了海拔对重金属浓度的影响，然后分析重金属一般的物理迁移、生物迁移等传播特征，由此我们建立二维数据差值与拟合模型，用IDW 算法求出足够精细分割的节点上的重金属浓度值，采用计算机绘图，以曲代直拟合模拟确定每种金属污染源的位置，最后结合问题一中功能图和问题二中污染的主要原因检验模型的合理信。
针对问题四：我们通过分析问题三中模型的优缺点，发现问题三的模型只能确定某一时刻的重金属的污染源位置和浓度，并不能确定重金属污染的长期发展 趋势，因此我们考虑收集各重金属在不同时刻的浓度，进而建立时间序列模型解决该问分析出的问题，并利用模型更好的研究城市地质环境的演变模式。
关键字： 降维因子分析 聚类分析 显著性检验 二维拟合 IDW算法 时间序列模型
一、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、模型假设
1、假设题目所给数据都真实有效；
2、假设该城区远离人群及工业活动的自然区还没被人类活动所影响； 3、假设我们不考虑数据中的异常值，并将其剔除； 4、假设该城区在二维上为规则的矩形区域； 5、假设不考虑计算机二维拟合所带来的误差。
三、符号说明
Q L ：下四分位数 Q U ：上四分位数 I g e :o ；地积累指数
不能同区浓度的均 Cn :各类重金属元素在功 值
K :考虑各地岩石差异可能会引起背景值的变动而取的系数(我们取K =1. 5);
元平素均值 Bn :附件3里不同重金属污染的 ；
z i ：每种金属的第 值i 个测定的浓度 B ：延迟算子 ∇n x t ：x t 的n 阶差分
Φ（B ） ：自回归系数多项式
Θ（B ）：移动平均系数多项式 ∇n x t ：x t 的n 阶差分
四、问题分析及模型的建立与求解
4.1 问题1的分析及模型的建立与求解
4.1.1数据的处理
根据非参数统计分析中的异常值检测方法，称Q L -3⨯IQL 和Q U +3⨯IQL 为两个外篱笆（Outer Fences），其中Q L 、Q U 分别为下四分位数和上四分位数，
IQL =Q L -Q U , 收集的数据中处于外篱笆之外的数称为强异常值，我们需要将此类数据剔除。
在该题目中我们对不同功能区不同重金属分别考虑强异常值，通过计算求得而附录2数据全为正值，故只需考虑剔除大于Q U +3⨯IQL Q L -3⨯IQL 全为负值，
的异常值后的数据，通过计算得到不同功能区不同重金属浓度的强异常值（见表1），然后得到所有剔除异常值后的数据（见附录5），然后我们根据这个筛选过后的数据进行问题分析与解决。
表1
4.1.2 问题1的分析
本题要求绘出8种主要重金属元素在该城区的空间分布和分析不同区域内重金属的污染程度，我们可以根据附件1用计算机仿真技术得出该城区的地形图及其MATLAB 程序（见附录1）。事实上，如果我们将8种主要重金属元素一起考虑并与地形图相结合绘制出其空间分布图的话，我们可以发现该空间图是相当混乱的，根本没办法进行分析。所以，我们将各种重金属分开考虑，然后绘制出各种重金属元素的空间分布图，并根据其分布图可以大致得出不同区域的重金属污染程度。为了使我们的结果更加可信，我们又对上述模型进行改进得出确切的各功能区的污染程度。然而，我们在绘图是发现有些重金属浓度有明显的异常值，破坏了其真实的空间分布图，所以我们在绘图分析时舍弃了这些异常值。 4.1.3 问题1模型一的建立与求解
我们根据对题目附件1所给数据进行分析，用MATLAB 编程（见附录2）绘出各个功能区在该城区内的分布（见图1，用不同的颜色表示不同的功能区），使各功能区的位置更加直观。我们又对剔除异常值后的数据进行分析用MATLAB 编程（见附录3）绘制出8种重金属元素分别在该城区的不同功能区的空间分布图（颜色越趋向与红色表示浓度越深），各种重金属元素对应一个图，本文仅以重金属Cd 为例（见图2），其余各图可由附录2的编程类似可得。然后，将两图进行对比分析，该城区内不同区域重金属的污染程度。
图1
图
2
我们根据图1可得各个功能区的大致范围，然后根据图2的颜色参照图1大致可判断生活区、工业区、山区、交通区和公园绿地区Cd 的污染程度为：中、中、轻、中和轻2。各区域中其余各重金属污染元素的污染程度可类似分析得到，我们在此不再赘述。
4.1.4 问题1模型二的建立与求解
为了增加我们分析结果的可信度，我们又用SPSS 软件计算出不同功能区的8种重金属元素的各类统计量，再采用MULLER 地积累指数法计算出地积累指数，根据其分级标准进一步确定其污染程度。
我们根据附件2数据通过SPSS 软件求得5个功能区的8种重金属污染元素
的浓度平均值（见表2):
接着我们采用MULLER 地积累指数法来研究重金属污染的程度，其表达式为：Igeo =log 2[Cn
/(K ⨯Bn )]，式中Cn 是上表中元素n 在各功能区的浓度的平均值；Bn 是附件3所给的各种元素的背景值；K 为考虑各地岩石差异可能会引起背景值的变动而取的参数（本论文中我们取K 值为1）。
然后根据上表和地积累指数公式计算得8种重金属污染元素的地积累指数
Igeo （见表3）：
根据MULLER 地积累指数分级标准即：
当Igeo ≤0时，为污染； 当02时，为重度污染。
结合表2，对比上述分级标准我们得出了该城区内不同区域的重金属的污染程度（见表4）：
4.2 问题2的分析及模型的建立与求解 4.2.1 问题2的分析
通过我们对数据的初步分析并结合实际，我们发现各金属之间存在着一定的相关关系，我们考虑对8种金属首先进行因子分析，初步确定各金属之间的相关系数，然后通过聚类分析将金属分为几类，进而得出重金属污染的原因。 4.2.1 问题2的模型建立与求解
下面我们利用SPSS 软件对筛选过后的数据（见附录5）进行降维因子分析得到各金属之间的相关系数，其相关矩阵如下表（见表5）：
从表中我们可以看出Zn 和Cd 、Cu 、Pb 的相关系数分别为0.670、0.751、0.720，这些数据相对1来讲都比较高，说明Zn 和Cd 、Cu 、Pb 有较强的相关性，这表明这4种重金属的污染来自于同一原因。而Ni 与As 、Cr 的相关系数分别为0.615、0.706，这说明Ni 与As 、Cr 有较强的相关性，这表明这3种重金属的污染来自于同一原因。对于Hg 我们再通过进一步的分析寻找这种重金属污染的原因。
下面我们继续利用SPSS 软件对筛选过的数据（见附录4）进行聚类分析，我们采用相关系数质心聚类法进行聚类分析，得到如下的树状图（见图3）。
图3
由图3可以得出Zn 、Cd 、Cu 、Pb 是一类的，Ni 、As 、Cr 是另一类的，这与我们在上面的降维因子分析的结果是一致的，这说明我们的分析是正确的，而聚类分析将Hg 划分为Ni 、As 、Cr 的一类结合实际是可信的。有研究表明，Pb 主要来源于汽车尾气排放，Cu 、Zn 、Cd 主要来源于汽油燃烧、车体的磨损等，交通活动产生的灰尘经外力扬起沉积在周围土壤中，导致土壤中Pb 、Zn 、Cu 、Cd 含量升高，故可断定Zn 、Cd 、Cu 、Pb 可能主要来源于交通活动，而Ni 、As 、Cr 、Hg 多是工厂废水废气中的化合物的元素，故可能来自于工业活动。
因此通过数据分析，特别上面我们进行的降维因子分析和聚类分析，我们可以断定Zn 、Cd 、Cu 、Pb4种重金属污染的主要原因是交通活动，而Ni 、As 、Cr 、Hg 四种重金属污染的主要原因是工业活动。 4.3 问题3的分析及模型的建立与求解 4.3.1 问题3的分析
由以上两问我们知道了每种重金属的污染程度，又用聚类分析说明了重金属的主要原因，又知重金属一般只能发生各种形态的相互转化和分散、富集过程（即迁移），土壤表层的重金属有可能发生水平迁移，如生物迁移等，题目给出了该地区319组空间三维地点、不同重金属的含量（采样），首先我们对海拔做与各种重金属污染关系的显著性检验，从而排除海拔这一维的影响，由此我们可以建立二维数据差值与拟合模型，求出足够精细分割的节点上的重金属浓度值，用MATLAB 编程，采用计算机绘图、以曲代直、拟合找出每种金属污染源的位置，最后需要对结果进行必要的检验。
4.3.2 问题3的模型建立与求解
对于海拔与各种重金属污染之间关系的显著性检验：我们通过分析数据发现海拔只是处在0-350左右，对于各重金属元素的影响较小，所以我们据此对海拔
与各重金属的关系进行显著性检验；首先我们做出原假设：H 0为海拔与重金属污染程度没有关系，备择假设H 1为海拔与重金属污染程度没有关系。
据此我们使用SPSS 软件的配对样本t 检验（置信区间为95%）我们得到了海拔与各类重金属污染之间的相关系数（见表6）。
表6：成对样本相关系数
由上表中的相关系数可以看出，海拔与各重金属元素的浓度的相关系数在显著性水平水平小于5%（表中Siq 均小于0.05）时均为负值，故海拔与各重金属的污染程度没有相关关系，因此我们不能拒绝原假设，即海拔对各重金属的污染程度几乎没有影响，可以忽略，至此我们排出了海拔这个可能影响各种金属污染程度的因子。
根据实际情况我们认为重金属浓度曲面是光滑的，曲面的一阶、二阶导数是连续的，因为在现实中重金属污染源产生后对周围区域污染在不发生异常情况下不存在陡然影响的情况。另外，重金属污染浓度是一个按区域来划分的变量，在某个位置与其周围区域的浓度是相互依赖的，但这种依赖的作用 随距离的增大而减小，即每一个已知的数据点（采样点）影响周围的每一个未知点，一个给定的数据点（采样点）已未知点越近作用就越大。我们根据已知散乱的数据点简单易于计算的函数z =f (x , y ) （已知曲面）(x i , y i , z i )构造一个具有一定光滑性、
作为z =g (x , y ) （未知曲面）的近似，使曲面z =f (x , y ) 通过已知数据点，即
f (x i , y i ) =z i ，然后计算点(x , y ) 的函数值f (x , y ) 作为g (x , y ) 的近似值。
]⨯[0, 18000]是x ，y 平面上的一个规则的矩我们假定该城区为R :[0, 21000
形区域，然后在x 轴和y 轴上取定分割
∆x :0=x 0
∆y :0=y 0
R 上的一个矩形网络分割
∆=∆x ⨯∆y ，将R 分成m 个子矩形R ij :[x i -1, x i ]⨯[y i -1, y i ]，直线x =x i (i =0, 1, , n )
和直线y =y j (j =0, 1, , m ) 称为∆分割的两簇网格线，共有(m +1)(n +1) 个交点（节点），考虑到该城区区域x 、y 的大小，我们规定在x 、y 方向分别剖分290
2100018000
=100（m ）=100（m ）和190个区间，区间间隔dx =，dy =。 210180
采用地球科学上的反距离权重法IDW ，已知点(x i , y i ) 的值为
已知点⎧z i ,
⎪319
，推测点(x , y ) 的值为z ，有f (x , y ) =⎨w z , 未知点 z i (i =1, 2, , 31) 9
∑i i ⎪⎩i =1
其中权重为w i =
1/d i
n i =1
p p i
∑1/d
所有未知点(x , y ) 处的值由各, d i =(x -x i ) 2+(y -y i ) 2，
已知点加权得到，权重为(x , y ) 到各点距离的p 次方反比，p 作为距离(x , y ) 近的反映了未知点位于已知点的偏倚情况：p 越大，(x i , y i ) 作用的相对大于的影响值，
当(x i , y i ) 距离(x , y ) 越近其相对作用越大，越远相对作用越小。本题考虑到分析重金属污染源的位置，我们取p=2。另外，我们知道第二问中当地积累指数的值
大于2时，该金属为重度污染，因此我们结合附件3中的背景值以及地积累指数法的公式Igeo =log 2[Cn /(K ⨯Bn )]（K=1），我们可以得出各种金属达到重污染时的浓度的临界值（见表6）。二维拟合出一些大于这个临界值的区域，即为污染源的位置。
表6
我们用附件所给的不同地点(x , y ) Cd 重金属含量数据，然后编写了一个
matlab 程序，我们得到Cd 浓度的二维曲面拟合图及污染源位置图如下。（程序见附录4）
上图蓝色区域即为所求Cd 重金属的污染源位置（共30个），根据其坐标系我们确定其中两个Cd 污染源大概位置为[14000]⨯[12000]，, 14250, 13000。 [2000, 2900]⨯[2000, 3000]（单位：m ）
按照此模型可以确定余下七种污染源的位置，之后就可以做好对污染源的重点监测与控制，减少重金属对该城区的影响。 4.3.3 模型结果的检验
我们根据模型确定了Cd 重金属污染源的确切位置，以下我们检验此结果
的合理性。第一问中根据附件1所给的319个不同点所属的五个功能区，用MATLAB 拟合了该城区五个功能区的大致区域，我们发现由模型确定的Cd 重金属污染源区域基本都在交通区，在问题二中分析认为Cd 重金属主要来自交通区，主要来源为汽油、车体的磨损等，交通活动产生的灰尘经外力扬起沉积在周围土壤中，导致土壤中Cd 含量升高，故模型确定的污染源位置可信。
4.4 问题4的分析及模型的建立与求解 4.4.1 问题4的分析
首先我们分析模型3的优缺点，在改善模型3的基础上，为了更好地研究城
市地质环境的演变模式，我们考虑收集各种金属浓度在不同时间的浓度值，建立新的时间序列模型，使模型更能反映该城区的不同区域的不同重金属的污染程度和污染源的位置。另外，由化学知识我们也可考虑重金属元素所在的化合物，因为不同的重金属化合物有着不同的毒性（如甲基汞比无机汞污染严重、6价铬比4价铬污染严重）。在此我们只考虑前者建立模型。
4.4.2 问题4的模型建立与求解
模型3的优点：操作简单，在充分利用已知点的信息的情况下，能够准确确定不同重金属的污染源的位置（区域）； 模型3的缺点：只能确定采样那一时刻的污染源的位置，不能准确反映该地区连续的污染源的位置，随机性较大，不能预测重金属污染的发展，不能去找到合适的切入点措施去治理。
下面我们说明收集各种重金属不同时间时的浓度（如2010年9月-2011年9月）后建立关于重金属污染的时间序列模型。由于城区不同区域受到了人类活动所造成的污染，故各重金属污染浓度的时间序列一定是非平稳时间序列。 我们首先要利用应用时间序列所学知识将其化为平稳时间序列，具体方法是对样本时间序列进行随机时序分析，采用差分运算的方法充分提取序列中蕴藏的确定性信息，找出序列发展的长期趋势，进而确定污染治理的方案。
首先我们假设所收集的某重金属（如As ）随时间变化的浓度序列为，根据Cramer 分解定理，该非平稳序列可分解为： X （t x 1, x 2, x 3 ）
j
x t =∑βj ⋅t +ψ(B ) a t ， （1）
j =0d
式中，d
在此分解定理下，对{xt }的d 阶差分就可以将该浓度的非平稳时间序列中蕴含的确定性信息和发展的长期趋势充分提取。其中展开1阶差分，有
∇x t =x t -x t -1等价于x t =x t -1+∇x t ，显然这就是一个自回归过程。同样展开二阶差分有∇2x t =∇(∇x t ) =∇(x t -x t -1) =2x t -1-x t -2+∇2x t 也是一个自回归过程。故展开任意一个关于浓度的非平稳时间序列{xt }的d 阶差分，都可得到一个d 阶的自
d
i
回归过程： x t =∑(-1) i +1C d x t -i +∇d x t ， （2）
i =1
i 等价于 ∇x t =∑(-1) i C d x t -i ， （3）
d
i =0
d
利用这种自回归方式我们就将各重金属的浓度非平稳序列转化为平稳时间序列，称其为差分平稳时间序列，在此我们将差分平稳序列记为Y 。 （）t y 1, y 2, y 3, 下面我们通过ARIMA 模型对上面已化为差分平稳时间序列的重金属浓度序列进行预测。通过（3）式我们可以看到差分后的序列{∇d x t }等于原序列的加权和，故又可以拟合自回归移动平均（ARMA ）模型，我们称这种为求和自回归移动平均（ARIMA ）模型。我们引入后移算子B 后，我们得到各重金属浓度的ARIMA 模型为： Φ（B ） (1-B ) d y t =Θ(B ) εt ， （4）也可利用历史观察值的线性函数表示，在此不再赘述（（4）式中的Φ（B ）为d 阶自回归系数多项式）。（4）式中ψ1，ψ2，ψ3， 的值由如下等式确定： Φ（B ） (1-B ) d ψ(B ) =Θ(B ) ， （5）
2p 其中Φ（B ） 称为p 阶自回归系数多项式，Θ（B ） =1-φ（B ）-φB - φB 12p
=1-θ1B -θ2B 2- -θq B q 称为q 阶移动平均系数多项式。我们如果类似于AMRA 模
*
（B ）型将Φ记为广义自相关函数，则有：
*d
Φ （B ）=Φ（B ）（1-B ）=1-φ1B -φ2B 2- （6）
其中（6）式中的B d 为延迟d 期后的转换系数，φn 为系数。容易验证
ψ1，ψ2，ψ3， 满足如下递推公式：
⎧ψ1=φ1-θ1
⎪ψ=φψ+φ-θ⎪21122
⎨， （7）
⎪⎪⎩ψj =φ1ψj -1+ +φp +d ψj -p -d -θj
⎧0，j
其中ψj =⎨
⎩1, j =0,
θj =0, j >q
。
那么y t +s 的真实值为：
y t +s =(εt +s +ψ1εt +s -1+ +ψs -1εt +1) +(ψs εt +ψs +1εs -1+ ) ， 由于εt +s , εt +s -1, 的不可获得性，所以
∧
y t +的估计值（即预报值）只能为：
**
y t (s ) =ψ0εt +ψ1*εt -1+ψ2εt -2+ ，
真实值与预报值的均方误差为：
E [y t +s -y t (s )]=(1+ψ+ +ψ) σε+
∧
2
2
1
2s -1
2
∑(ψ
j =0
∞
*22-ψ （8） s +j j ) σε，
要使均值误差即上述（8）式最小，当且仅当ψj *=ψs +j ，所以在均方误差最小原则下，s 期的预测值为：
y t (s ) =ψs εt +ψs +1εt -1+ψs +2εt -2+ ， （9） 此时， s期的预报误差为：
e t (s ) =εt +s +ψ1εt +s -1+ +ψs -1εt +1 ， （10） 此时我们就得到了我们的预测真实值，它等于预测值加上预报误差，即（9）+
（10）：
∧
y t +s =(εt +s +ψ1εt +s -1+ +ψs -1εt +1) +(ψs εt +ψs +1εs -1+ 。
综上，通过对所预测的各重金属不同时间的浓度变化我们或者分析求得的数据，或者进行作图，都可以找到各重金属元素的污染程度的变化，进而找到污染的原因和污染的长期发展趋势，从而找到切入点，进而更好地研究城市地质环境的演变模式，减弱人类活动对环境的重金属元素污染程度影响。
五、模型的评价与推广
5.1 模型的评价
优点：
(1) 本文针对不同问题插入了一些图表增加论文的可读性，使我们的模型更加直观易懂；
(2) 对于模型的重金属的空间分布图，我们用MATLEB 绘制其浓度等值线与空间分布结合图然后和功能区分布图对照使问题更直观，由图便可以直接得出各重金属在不同功能区的污染程度；
(3) 我们首先对原始数据进行了处理，舍弃了一些明显的异常值，使我们的模型更加符合实际情况，增加了模型的可信度；
(4) 充分利用MATLAB ，SPSS 等软件进行编程求解，所得误差较小，数据准确合理；
(5) 该模型实用性强，对现实有较强的指导意义；
(6) 模型中应用了降维因子分析、聚类分析、二维拟合和时间序列分析等数学方法解决问题，得到了较好的结果，是我们的结论具有了较强的理论依据； 缺点：
(1) 我们由原始数据进行建模得出结果，但没有搜集实际数据进行检验； (2) 我们忽略了海拔的因素对浓度的影响，而且在二维拟合计算机模拟中，模拟结果会存在误差 。 5.2 模型的推广
此模型虽然建立在分析重金属污染的情况下，但它同样也适用于研究其他污染物，比如大气污染，另外，问题三建立的二维拟合模型，在现实世界中只要对二维变量分析区域情况都能适用，比如已知海域某点的深度确定整个海底的空间三维图，对近海航行安全起到一定的作用。
六、参考文献
[1] 王燕. 应用时间序列分析(第二版) [M],北京:中国人民大学出版社.2008.12 [2] 肖华勇. 实用数学建模与软件应用[M],西安:西北工业大学出版社.2008.11 [3] 王静龙, 梁小筠. 非参数统计分析[M],北京:高等教育出版社.2006.4 [4] 刘会灯, 朱飞.MATLAB 编程基础与典型应用[M],北京:人民邮电出版社.2008.7
[5] 薛薇. 统计分析与SPSS 的应用(第三版)[M],北京:中国人民大学出版社.2011.1
[6]百度百科. 重金属污染http://baike.baidu.com/view/282970.htm,2011.9.9 [7]百度百科. 各种化学元素http://baike.baidu.com/view/30675.htm,2011.9.9
七、附录
附录1：该城区的地形图
其程序为：
clc,clear;
A=[];%输入附件1中的第2、3、4、5列数据 x=A(:,1);y=A(:,2);z=A(:,3);
scatter(x,y,5,z,0:10:30000)%散点图 figure
[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x),200)',linspace(min(y),max(y),200),'v4' ); %插值
pcolor(X,Y,Z);shading interp %伪彩色图 figure,contourf(X,Y,Z) %等高线图 figure,surf(X,Y,Z)%三维曲面
附录2：各功能区的分布图的MATLAB 编写程序
clc,clear ;
D=[];%输入附件1中的x(m),y(m),海拔(m)和功能区数据 x=D(:,1); y=D(:,2); z=D(:,3); c=D(:,4);
xi=linspace(min(x),max(x),100);%差值 yi=linspace(min(y),max(y),100); [xi,yi]=meshgrid(xi,yi); zi=griddata(x,y,z,xi,yi); ci=griddata(x,y,c,xi,yi); marker={'*', 'o' , 's' , '^', 'p' }; color={'k' , 'r' , 'y' , 'c' , 'b' }; figure
h=surf(xi,yi,zi); set(h,'cdata' ,ci);
colormap hsv
title(' 功能区分布图' ) xlabel('X' ) ylabel('Y' ) colorbar hidden off hold on for i=1:5
loc=c==i;
plot3(x(loc),y(loc),z(loc),marker{i},'markerfacecolor' ,color{i}); %对不同功能区绘不同的颜色 end
附录3：各重金属污染元素的浓度等值线图与功能区结合的MATLAB 编写程序 clc,clear;
A=[];%输入某一重金属污染元素的x 、y 和相应的浓度值； x=A(:,1);y=A(:,2);z=A(:,3);
xi=linspace(min(x),max(x),100); yi=linspace(min(y),max(y),100); [xi,yi]=meshgrid(xi,yi); zi=griddata(x,y,z,xi,yi);
figure,contourf(xi,yi,zi) %等高线图 hold on;
A=[] %输入功能1的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'p' ); grid on ; hold on ;
A=[] %输入功能2的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'ro' ); hold on ;
A=[] %输入功能3的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'k*'); hold on ;
A=[] %输入功能4的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'gd' ); hold on ;
A=[] %输入功能5的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'m
附录4：Cd 浓度的二维曲面拟合图及污染源位置图的MATLEB 编写程序
clc,clear;
data=[] %输入某一重金属污染元素的x 、y 和相应的浓度值； plot(data(:,1),data(:,2),'*'); xx=[0,0;
29000,0; 29000,19000; 0,19000; 0,0]; hold on ;
plot(xx(:,1),xx(:,2),'r' ); pause
d=zeros(319,1); w=zeros(319,1);
m=290; %x的区间数 n=190; %y的区间数
dx=29000/m; %x的间隔 dy=19000/n; %y的间隔
point=(m+1)*(n+1); %总点数
T=zeros(point,3); %储存各点数据 X=zeros(n+1,m+1); Y=zeros(n+1,m+1); Z=zeros(n+1,m+1); p=2; k=0; tp=0;
for ix=1:m+1 for iy=1:n+1
x=(ix-1)*dx; %x的坐标 y=(iy-1)*dy; %y的坐标
for i=1:319 %（x,y ）格点到各已知各点的距离
d(i)=sqrt((x-data(i,1))^2+(y-data(i,2))^2); w(i)=1/d(i)^p; %各点的权重 end ;
s=sum(w);
w=w/s; %权值归一化
z=sum(data(:,3).*w); %求得网格点 X(iy,ix)=x; Y(iy,ix)=y; Z(iy,ix)=z;
k=k+1; %储存各点(x,y,z)的坐标 T(k,1)=x;
T(k,2)=y; T(k,3)=z; if z>520 tp=tp+1; D(tp,1)=x;
D(tp,2)=y; %获得污染源的位置 end end ; end ;
subplot(2,1,1);
surf(X,Y,Z) %作曲面图 subplot(2,1,2);
contour(X,Y,Z); %作等值线图 hold on
plot(D(:,1),D(:,2),'*') %作出污染源位置的点 title(' 某金属浓度等值线和污染源位置图' ); grid on hold off
承 诺 书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 淮北师范大学 参赛队员 (打印并签名) ：1. 孙丰凯
指导教师或指导教师组负责人 (打印并签名) ： 李昌文
日期： 年 月 日
赛区评阅编号（由赛区组委会评阅前进行编号）：
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
城市表层土壤重金属污染分析及污染源位置的确定
摘要
土壤重金属污染问题是当今环境科学研究的重要内容, 特别是城市土壤重金属污染受到更多重视。城市人口集中、工业发达、交通拥挤等因素导致土壤中重金属不仅含量高而且种类较多。一旦城市土壤被重金属污染, 城市居民的健康将受到危害, 因此, 世界上许多国家, 尤其是欧美一些发达国家十分关注城市土壤重金属的污染状况。
针对问题一：模型一：我们根据数据用MATLAB 软件编程绘制出各功能区的空间分布图（用不同的颜色表示不同的区域）和各种重金属元素的浓度与空间分布的结合图，并将这两个图进行对比大致判断不同区域中各种重金出属在土壤表层的污染程度。我们又对模型一进行改进建立模型二更加准确的确定不同功能区的污染程度，对此我们用SPSS 软件计算出不同功能区的8种重金属元素的各类统计量，再采用MULLER 地积累指数法计算出地积累指数，根据其分级标准进一步确定其污染程度。
针对问题二：我们首先利用SPSS 软件对附件2的8种元素的浓度数据进行降维因子分析，得出各元素之间的相关矩阵，通过分析各个权重大致得出元素之间的关系，然后采用相关系数质心聚类法对原始数据进行聚类分析，当距离为15~20时主要分为3类，及重金属污染的主要原因有3种。
针对问题三：首先我们作海拔与各种重金属污染程度相关关系的显著性检验，排除了海拔对重金属浓度的影响，然后分析重金属一般的物理迁移、生物迁移等传播特征，由此我们建立二维数据差值与拟合模型，用IDW 算法求出足够精细分割的节点上的重金属浓度值，采用计算机绘图，以曲代直拟合模拟确定每种金属污染源的位置，最后结合问题一中功能图和问题二中污染的主要原因检验模型的合理信。
针对问题四：我们通过分析问题三中模型的优缺点，发现问题三的模型只能确定某一时刻的重金属的污染源位置和浓度，并不能确定重金属污染的长期发展 趋势，因此我们考虑收集各重金属在不同时刻的浓度，进而建立时间序列模型解决该问分析出的问题，并利用模型更好的研究城市地质环境的演变模式。
关键字： 降维因子分析 聚类分析 显著性检验 二维拟合 IDW算法 时间序列模型
一、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、模型假设
1、假设题目所给数据都真实有效；
2、假设该城区远离人群及工业活动的自然区还没被人类活动所影响； 3、假设我们不考虑数据中的异常值，并将其剔除； 4、假设该城区在二维上为规则的矩形区域； 5、假设不考虑计算机二维拟合所带来的误差。
三、符号说明
Q L ：下四分位数 Q U ：上四分位数 I g e :o ；地积累指数
不能同区浓度的均 Cn :各类重金属元素在功 值
K :考虑各地岩石差异可能会引起背景值的变动而取的系数(我们取K =1. 5);
元平素均值 Bn :附件3里不同重金属污染的 ；
z i ：每种金属的第 值i 个测定的浓度 B ：延迟算子 ∇n x t ：x t 的n 阶差分
Φ（B ） ：自回归系数多项式
Θ（B ）：移动平均系数多项式 ∇n x t ：x t 的n 阶差分
四、问题分析及模型的建立与求解
4.1 问题1的分析及模型的建立与求解
4.1.1数据的处理
根据非参数统计分析中的异常值检测方法，称Q L -3⨯IQL 和Q U +3⨯IQL 为两个外篱笆（Outer Fences），其中Q L 、Q U 分别为下四分位数和上四分位数，
IQL =Q L -Q U , 收集的数据中处于外篱笆之外的数称为强异常值，我们需要将此类数据剔除。
在该题目中我们对不同功能区不同重金属分别考虑强异常值，通过计算求得而附录2数据全为正值，故只需考虑剔除大于Q U +3⨯IQL Q L -3⨯IQL 全为负值，
的异常值后的数据，通过计算得到不同功能区不同重金属浓度的强异常值（见表1），然后得到所有剔除异常值后的数据（见附录5），然后我们根据这个筛选过后的数据进行问题分析与解决。
表1
4.1.2 问题1的分析
本题要求绘出8种主要重金属元素在该城区的空间分布和分析不同区域内重金属的污染程度，我们可以根据附件1用计算机仿真技术得出该城区的地形图及其MATLAB 程序（见附录1）。事实上，如果我们将8种主要重金属元素一起考虑并与地形图相结合绘制出其空间分布图的话，我们可以发现该空间图是相当混乱的，根本没办法进行分析。所以，我们将各种重金属分开考虑，然后绘制出各种重金属元素的空间分布图，并根据其分布图可以大致得出不同区域的重金属污染程度。为了使我们的结果更加可信，我们又对上述模型进行改进得出确切的各功能区的污染程度。然而，我们在绘图是发现有些重金属浓度有明显的异常值，破坏了其真实的空间分布图，所以我们在绘图分析时舍弃了这些异常值。 4.1.3 问题1模型一的建立与求解
我们根据对题目附件1所给数据进行分析，用MATLAB 编程（见附录2）绘出各个功能区在该城区内的分布（见图1，用不同的颜色表示不同的功能区），使各功能区的位置更加直观。我们又对剔除异常值后的数据进行分析用MATLAB 编程（见附录3）绘制出8种重金属元素分别在该城区的不同功能区的空间分布图（颜色越趋向与红色表示浓度越深），各种重金属元素对应一个图，本文仅以重金属Cd 为例（见图2），其余各图可由附录2的编程类似可得。然后，将两图进行对比分析，该城区内不同区域重金属的污染程度。
图1
图
2
我们根据图1可得各个功能区的大致范围，然后根据图2的颜色参照图1大致可判断生活区、工业区、山区、交通区和公园绿地区Cd 的污染程度为：中、中、轻、中和轻2。各区域中其余各重金属污染元素的污染程度可类似分析得到，我们在此不再赘述。
4.1.4 问题1模型二的建立与求解
为了增加我们分析结果的可信度，我们又用SPSS 软件计算出不同功能区的8种重金属元素的各类统计量，再采用MULLER 地积累指数法计算出地积累指数，根据其分级标准进一步确定其污染程度。
我们根据附件2数据通过SPSS 软件求得5个功能区的8种重金属污染元素
的浓度平均值（见表2):
接着我们采用MULLER 地积累指数法来研究重金属污染的程度，其表达式为：Igeo =log 2[Cn
/(K ⨯Bn )]，式中Cn 是上表中元素n 在各功能区的浓度的平均值；Bn 是附件3所给的各种元素的背景值；K 为考虑各地岩石差异可能会引起背景值的变动而取的参数（本论文中我们取K 值为1）。
然后根据上表和地积累指数公式计算得8种重金属污染元素的地积累指数
Igeo （见表3）：
根据MULLER 地积累指数分级标准即：
当Igeo ≤0时，为污染； 当02时，为重度污染。
结合表2，对比上述分级标准我们得出了该城区内不同区域的重金属的污染程度（见表4）：
4.2 问题2的分析及模型的建立与求解 4.2.1 问题2的分析
通过我们对数据的初步分析并结合实际，我们发现各金属之间存在着一定的相关关系，我们考虑对8种金属首先进行因子分析，初步确定各金属之间的相关系数，然后通过聚类分析将金属分为几类，进而得出重金属污染的原因。 4.2.1 问题2的模型建立与求解
下面我们利用SPSS 软件对筛选过后的数据（见附录5）进行降维因子分析得到各金属之间的相关系数，其相关矩阵如下表（见表5）：
从表中我们可以看出Zn 和Cd 、Cu 、Pb 的相关系数分别为0.670、0.751、0.720，这些数据相对1来讲都比较高，说明Zn 和Cd 、Cu 、Pb 有较强的相关性，这表明这4种重金属的污染来自于同一原因。而Ni 与As 、Cr 的相关系数分别为0.615、0.706，这说明Ni 与As 、Cr 有较强的相关性，这表明这3种重金属的污染来自于同一原因。对于Hg 我们再通过进一步的分析寻找这种重金属污染的原因。
下面我们继续利用SPSS 软件对筛选过的数据（见附录4）进行聚类分析，我们采用相关系数质心聚类法进行聚类分析，得到如下的树状图（见图3）。
图3
由图3可以得出Zn 、Cd 、Cu 、Pb 是一类的，Ni 、As 、Cr 是另一类的，这与我们在上面的降维因子分析的结果是一致的，这说明我们的分析是正确的，而聚类分析将Hg 划分为Ni 、As 、Cr 的一类结合实际是可信的。有研究表明，Pb 主要来源于汽车尾气排放，Cu 、Zn 、Cd 主要来源于汽油燃烧、车体的磨损等，交通活动产生的灰尘经外力扬起沉积在周围土壤中，导致土壤中Pb 、Zn 、Cu 、Cd 含量升高，故可断定Zn 、Cd 、Cu 、Pb 可能主要来源于交通活动，而Ni 、As 、Cr 、Hg 多是工厂废水废气中的化合物的元素，故可能来自于工业活动。
因此通过数据分析，特别上面我们进行的降维因子分析和聚类分析，我们可以断定Zn 、Cd 、Cu 、Pb4种重金属污染的主要原因是交通活动，而Ni 、As 、Cr 、Hg 四种重金属污染的主要原因是工业活动。 4.3 问题3的分析及模型的建立与求解 4.3.1 问题3的分析
由以上两问我们知道了每种重金属的污染程度，又用聚类分析说明了重金属的主要原因，又知重金属一般只能发生各种形态的相互转化和分散、富集过程（即迁移），土壤表层的重金属有可能发生水平迁移，如生物迁移等，题目给出了该地区319组空间三维地点、不同重金属的含量（采样），首先我们对海拔做与各种重金属污染关系的显著性检验，从而排除海拔这一维的影响，由此我们可以建立二维数据差值与拟合模型，求出足够精细分割的节点上的重金属浓度值，用MATLAB 编程，采用计算机绘图、以曲代直、拟合找出每种金属污染源的位置，最后需要对结果进行必要的检验。
4.3.2 问题3的模型建立与求解
对于海拔与各种重金属污染之间关系的显著性检验：我们通过分析数据发现海拔只是处在0-350左右，对于各重金属元素的影响较小，所以我们据此对海拔
与各重金属的关系进行显著性检验；首先我们做出原假设：H 0为海拔与重金属污染程度没有关系，备择假设H 1为海拔与重金属污染程度没有关系。
据此我们使用SPSS 软件的配对样本t 检验（置信区间为95%）我们得到了海拔与各类重金属污染之间的相关系数（见表6）。
表6：成对样本相关系数
由上表中的相关系数可以看出，海拔与各重金属元素的浓度的相关系数在显著性水平水平小于5%（表中Siq 均小于0.05）时均为负值，故海拔与各重金属的污染程度没有相关关系，因此我们不能拒绝原假设，即海拔对各重金属的污染程度几乎没有影响，可以忽略，至此我们排出了海拔这个可能影响各种金属污染程度的因子。
根据实际情况我们认为重金属浓度曲面是光滑的，曲面的一阶、二阶导数是连续的，因为在现实中重金属污染源产生后对周围区域污染在不发生异常情况下不存在陡然影响的情况。另外，重金属污染浓度是一个按区域来划分的变量，在某个位置与其周围区域的浓度是相互依赖的，但这种依赖的作用 随距离的增大而减小，即每一个已知的数据点（采样点）影响周围的每一个未知点，一个给定的数据点（采样点）已未知点越近作用就越大。我们根据已知散乱的数据点简单易于计算的函数z =f (x , y ) （已知曲面）(x i , y i , z i )构造一个具有一定光滑性、
作为z =g (x , y ) （未知曲面）的近似，使曲面z =f (x , y ) 通过已知数据点，即
f (x i , y i ) =z i ，然后计算点(x , y ) 的函数值f (x , y ) 作为g (x , y ) 的近似值。
]⨯[0, 18000]是x ，y 平面上的一个规则的矩我们假定该城区为R :[0, 21000
形区域，然后在x 轴和y 轴上取定分割
∆x :0=x 0
∆y :0=y 0
R 上的一个矩形网络分割
∆=∆x ⨯∆y ，将R 分成m 个子矩形R ij :[x i -1, x i ]⨯[y i -1, y i ]，直线x =x i (i =0, 1, , n )
和直线y =y j (j =0, 1, , m ) 称为∆分割的两簇网格线，共有(m +1)(n +1) 个交点（节点），考虑到该城区区域x 、y 的大小，我们规定在x 、y 方向分别剖分290
2100018000
=100（m ）=100（m ）和190个区间，区间间隔dx =，dy =。 210180
采用地球科学上的反距离权重法IDW ，已知点(x i , y i ) 的值为
已知点⎧z i ,
⎪319
，推测点(x , y ) 的值为z ，有f (x , y ) =⎨w z , 未知点 z i (i =1, 2, , 31) 9
∑i i ⎪⎩i =1
其中权重为w i =
1/d i
n i =1
p p i
∑1/d
所有未知点(x , y ) 处的值由各, d i =(x -x i ) 2+(y -y i ) 2，
已知点加权得到，权重为(x , y ) 到各点距离的p 次方反比，p 作为距离(x , y ) 近的反映了未知点位于已知点的偏倚情况：p 越大，(x i , y i ) 作用的相对大于的影响值，
当(x i , y i ) 距离(x , y ) 越近其相对作用越大，越远相对作用越小。本题考虑到分析重金属污染源的位置，我们取p=2。另外，我们知道第二问中当地积累指数的值
大于2时，该金属为重度污染，因此我们结合附件3中的背景值以及地积累指数法的公式Igeo =log 2[Cn /(K ⨯Bn )]（K=1），我们可以得出各种金属达到重污染时的浓度的临界值（见表6）。二维拟合出一些大于这个临界值的区域，即为污染源的位置。
表6
我们用附件所给的不同地点(x , y ) Cd 重金属含量数据，然后编写了一个
matlab 程序，我们得到Cd 浓度的二维曲面拟合图及污染源位置图如下。（程序见附录4）
上图蓝色区域即为所求Cd 重金属的污染源位置（共30个），根据其坐标系我们确定其中两个Cd 污染源大概位置为[14000]⨯[12000]，, 14250, 13000。 [2000, 2900]⨯[2000, 3000]（单位：m ）
按照此模型可以确定余下七种污染源的位置，之后就可以做好对污染源的重点监测与控制，减少重金属对该城区的影响。 4.3.3 模型结果的检验
我们根据模型确定了Cd 重金属污染源的确切位置，以下我们检验此结果
的合理性。第一问中根据附件1所给的319个不同点所属的五个功能区，用MATLAB 拟合了该城区五个功能区的大致区域，我们发现由模型确定的Cd 重金属污染源区域基本都在交通区，在问题二中分析认为Cd 重金属主要来自交通区，主要来源为汽油、车体的磨损等，交通活动产生的灰尘经外力扬起沉积在周围土壤中，导致土壤中Cd 含量升高，故模型确定的污染源位置可信。
4.4 问题4的分析及模型的建立与求解 4.4.1 问题4的分析
首先我们分析模型3的优缺点，在改善模型3的基础上，为了更好地研究城
市地质环境的演变模式，我们考虑收集各种金属浓度在不同时间的浓度值，建立新的时间序列模型，使模型更能反映该城区的不同区域的不同重金属的污染程度和污染源的位置。另外，由化学知识我们也可考虑重金属元素所在的化合物，因为不同的重金属化合物有着不同的毒性（如甲基汞比无机汞污染严重、6价铬比4价铬污染严重）。在此我们只考虑前者建立模型。
4.4.2 问题4的模型建立与求解
模型3的优点：操作简单，在充分利用已知点的信息的情况下，能够准确确定不同重金属的污染源的位置（区域）； 模型3的缺点：只能确定采样那一时刻的污染源的位置，不能准确反映该地区连续的污染源的位置，随机性较大，不能预测重金属污染的发展，不能去找到合适的切入点措施去治理。
下面我们说明收集各种重金属不同时间时的浓度（如2010年9月-2011年9月）后建立关于重金属污染的时间序列模型。由于城区不同区域受到了人类活动所造成的污染，故各重金属污染浓度的时间序列一定是非平稳时间序列。 我们首先要利用应用时间序列所学知识将其化为平稳时间序列，具体方法是对样本时间序列进行随机时序分析，采用差分运算的方法充分提取序列中蕴藏的确定性信息，找出序列发展的长期趋势，进而确定污染治理的方案。
首先我们假设所收集的某重金属（如As ）随时间变化的浓度序列为，根据Cramer 分解定理，该非平稳序列可分解为： X （t x 1, x 2, x 3 ）
j
x t =∑βj ⋅t +ψ(B ) a t ， （1）
j =0d
式中，d
在此分解定理下，对{xt }的d 阶差分就可以将该浓度的非平稳时间序列中蕴含的确定性信息和发展的长期趋势充分提取。其中展开1阶差分，有
∇x t =x t -x t -1等价于x t =x t -1+∇x t ，显然这就是一个自回归过程。同样展开二阶差分有∇2x t =∇(∇x t ) =∇(x t -x t -1) =2x t -1-x t -2+∇2x t 也是一个自回归过程。故展开任意一个关于浓度的非平稳时间序列{xt }的d 阶差分，都可得到一个d 阶的自
d
i
回归过程： x t =∑(-1) i +1C d x t -i +∇d x t ， （2）
i =1
i 等价于 ∇x t =∑(-1) i C d x t -i ， （3）
d
i =0
d
利用这种自回归方式我们就将各重金属的浓度非平稳序列转化为平稳时间序列，称其为差分平稳时间序列，在此我们将差分平稳序列记为Y 。 （）t y 1, y 2, y 3, 下面我们通过ARIMA 模型对上面已化为差分平稳时间序列的重金属浓度序列进行预测。通过（3）式我们可以看到差分后的序列{∇d x t }等于原序列的加权和，故又可以拟合自回归移动平均（ARMA ）模型，我们称这种为求和自回归移动平均（ARIMA ）模型。我们引入后移算子B 后，我们得到各重金属浓度的ARIMA 模型为： Φ（B ） (1-B ) d y t =Θ(B ) εt ， （4）也可利用历史观察值的线性函数表示，在此不再赘述（（4）式中的Φ（B ）为d 阶自回归系数多项式）。（4）式中ψ1，ψ2，ψ3， 的值由如下等式确定： Φ（B ） (1-B ) d ψ(B ) =Θ(B ) ， （5）
2p 其中Φ（B ） 称为p 阶自回归系数多项式，Θ（B ） =1-φ（B ）-φB - φB 12p
=1-θ1B -θ2B 2- -θq B q 称为q 阶移动平均系数多项式。我们如果类似于AMRA 模
*
（B ）型将Φ记为广义自相关函数，则有：
*d
Φ （B ）=Φ（B ）（1-B ）=1-φ1B -φ2B 2- （6）
其中（6）式中的B d 为延迟d 期后的转换系数，φn 为系数。容易验证
ψ1，ψ2，ψ3， 满足如下递推公式：
⎧ψ1=φ1-θ1
⎪ψ=φψ+φ-θ⎪21122
⎨， （7）
⎪⎪⎩ψj =φ1ψj -1+ +φp +d ψj -p -d -θj
⎧0，j
其中ψj =⎨
⎩1, j =0,
θj =0, j >q
。
那么y t +s 的真实值为：
y t +s =(εt +s +ψ1εt +s -1+ +ψs -1εt +1) +(ψs εt +ψs +1εs -1+ ) ， 由于εt +s , εt +s -1, 的不可获得性，所以
∧
y t +的估计值（即预报值）只能为：
**
y t (s ) =ψ0εt +ψ1*εt -1+ψ2εt -2+ ，
真实值与预报值的均方误差为：
E [y t +s -y t (s )]=(1+ψ+ +ψ) σε+
∧
2
2
1
2s -1
2
∑(ψ
j =0
∞
*22-ψ （8） s +j j ) σε，
要使均值误差即上述（8）式最小，当且仅当ψj *=ψs +j ，所以在均方误差最小原则下，s 期的预测值为：
y t (s ) =ψs εt +ψs +1εt -1+ψs +2εt -2+ ， （9） 此时， s期的预报误差为：
e t (s ) =εt +s +ψ1εt +s -1+ +ψs -1εt +1 ， （10） 此时我们就得到了我们的预测真实值，它等于预测值加上预报误差，即（9）+
（10）：
∧
y t +s =(εt +s +ψ1εt +s -1+ +ψs -1εt +1) +(ψs εt +ψs +1εs -1+ 。
综上，通过对所预测的各重金属不同时间的浓度变化我们或者分析求得的数据，或者进行作图，都可以找到各重金属元素的污染程度的变化，进而找到污染的原因和污染的长期发展趋势，从而找到切入点，进而更好地研究城市地质环境的演变模式，减弱人类活动对环境的重金属元素污染程度影响。
五、模型的评价与推广
5.1 模型的评价
优点：
(1) 本文针对不同问题插入了一些图表增加论文的可读性，使我们的模型更加直观易懂；
(2) 对于模型的重金属的空间分布图，我们用MATLEB 绘制其浓度等值线与空间分布结合图然后和功能区分布图对照使问题更直观，由图便可以直接得出各重金属在不同功能区的污染程度；
(3) 我们首先对原始数据进行了处理，舍弃了一些明显的异常值，使我们的模型更加符合实际情况，增加了模型的可信度；
(4) 充分利用MATLAB ，SPSS 等软件进行编程求解，所得误差较小，数据准确合理；
(5) 该模型实用性强，对现实有较强的指导意义；
(6) 模型中应用了降维因子分析、聚类分析、二维拟合和时间序列分析等数学方法解决问题，得到了较好的结果，是我们的结论具有了较强的理论依据； 缺点：
(1) 我们由原始数据进行建模得出结果，但没有搜集实际数据进行检验； (2) 我们忽略了海拔的因素对浓度的影响，而且在二维拟合计算机模拟中，模拟结果会存在误差 。 5.2 模型的推广
此模型虽然建立在分析重金属污染的情况下，但它同样也适用于研究其他污染物，比如大气污染，另外，问题三建立的二维拟合模型，在现实世界中只要对二维变量分析区域情况都能适用，比如已知海域某点的深度确定整个海底的空间三维图，对近海航行安全起到一定的作用。
六、参考文献
[1] 王燕. 应用时间序列分析(第二版) [M],北京:中国人民大学出版社.2008.12 [2] 肖华勇. 实用数学建模与软件应用[M],西安:西北工业大学出版社.2008.11 [3] 王静龙, 梁小筠. 非参数统计分析[M],北京:高等教育出版社.2006.4 [4] 刘会灯, 朱飞.MATLAB 编程基础与典型应用[M],北京:人民邮电出版社.2008.7
[5] 薛薇. 统计分析与SPSS 的应用(第三版)[M],北京:中国人民大学出版社.2011.1
[6]百度百科. 重金属污染http://baike.baidu.com/view/282970.htm,2011.9.9 [7]百度百科. 各种化学元素http://baike.baidu.com/view/30675.htm,2011.9.9
七、附录
附录1：该城区的地形图
其程序为：
clc,clear;
A=[];%输入附件1中的第2、3、4、5列数据 x=A(:,1);y=A(:,2);z=A(:,3);
scatter(x,y,5,z,0:10:30000)%散点图 figure
[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x),200)',linspace(min(y),max(y),200),'v4' ); %插值
pcolor(X,Y,Z);shading interp %伪彩色图 figure,contourf(X,Y,Z) %等高线图 figure,surf(X,Y,Z)%三维曲面
附录2：各功能区的分布图的MATLAB 编写程序
clc,clear ;
D=[];%输入附件1中的x(m),y(m),海拔(m)和功能区数据 x=D(:,1); y=D(:,2); z=D(:,3); c=D(:,4);
xi=linspace(min(x),max(x),100);%差值 yi=linspace(min(y),max(y),100); [xi,yi]=meshgrid(xi,yi); zi=griddata(x,y,z,xi,yi); ci=griddata(x,y,c,xi,yi); marker={'*', 'o' , 's' , '^', 'p' }; color={'k' , 'r' , 'y' , 'c' , 'b' }; figure
h=surf(xi,yi,zi); set(h,'cdata' ,ci);
colormap hsv
title(' 功能区分布图' ) xlabel('X' ) ylabel('Y' ) colorbar hidden off hold on for i=1:5
loc=c==i;
plot3(x(loc),y(loc),z(loc),marker{i},'markerfacecolor' ,color{i}); %对不同功能区绘不同的颜色 end
附录3：各重金属污染元素的浓度等值线图与功能区结合的MATLAB 编写程序 clc,clear;
A=[];%输入某一重金属污染元素的x 、y 和相应的浓度值； x=A(:,1);y=A(:,2);z=A(:,3);
xi=linspace(min(x),max(x),100); yi=linspace(min(y),max(y),100); [xi,yi]=meshgrid(xi,yi); zi=griddata(x,y,z,xi,yi);
figure,contourf(xi,yi,zi) %等高线图 hold on;
A=[] %输入功能1的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'p' ); grid on ; hold on ;
A=[] %输入功能2的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'ro' ); hold on ;
A=[] %输入功能3的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'k*'); hold on ;
A=[] %输入功能4的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'gd' ); hold on ;
A=[] %输入功能5的x(m),y(m)值； x=[:,1];y=[:,2]; plot(x,y,'m
附录4：Cd 浓度的二维曲面拟合图及污染源位置图的MATLEB 编写程序
clc,clear;
data=[] %输入某一重金属污染元素的x 、y 和相应的浓度值； plot(data(:,1),data(:,2),'*'); xx=[0,0;
29000,0; 29000,19000; 0,19000; 0,0]; hold on ;
plot(xx(:,1),xx(:,2),'r' ); pause
d=zeros(319,1); w=zeros(319,1);
m=290; %x的区间数 n=190; %y的区间数
dx=29000/m; %x的间隔 dy=19000/n; %y的间隔
point=(m+1)*(n+1); %总点数
T=zeros(point,3); %储存各点数据 X=zeros(n+1,m+1); Y=zeros(n+1,m+1); Z=zeros(n+1,m+1); p=2; k=0; tp=0;
for ix=1:m+1 for iy=1:n+1
x=(ix-1)*dx; %x的坐标 y=(iy-1)*dy; %y的坐标
for i=1:319 %（x,y ）格点到各已知各点的距离
d(i)=sqrt((x-data(i,1))^2+(y-data(i,2))^2); w(i)=1/d(i)^p; %各点的权重 end ;
s=sum(w);
w=w/s; %权值归一化
z=sum(data(:,3).*w); %求得网格点 X(iy,ix)=x; Y(iy,ix)=y; Z(iy,ix)=z;
k=k+1; %储存各点(x,y,z)的坐标 T(k,1)=x;
T(k,2)=y; T(k,3)=z; if z>520 tp=tp+1; D(tp,1)=x;
D(tp,2)=y; %获得污染源的位置 end end ; end ;
subplot(2,1,1);
surf(X,Y,Z) %作曲面图 subplot(2,1,2);
contour(X,Y,Z); %作等值线图 hold on
plot(D(:,1),D(:,2),'*') %作出污染源位置的点 title(' 某金属浓度等值线和污染源位置图' ); grid on hold off