前言:
为了准备国赛,我开始了历年国赛论文的阅读。我把论文中的主要思想方法总结到下面,供日后查阅方便。文末会给出这道题相关的资源,包括题目、论文、数据以及求解使用的部分代码的链接。
一、2011年美赛A题的要求
相比美赛,国赛要中规中矩的多,不管是题目还是获奖论文。2011年美赛A题的要求如下【1】:
(1)给出8种主要重金属元素在该城区的空间分布,并分析该城区内不同区域重金属的污染程度。(按照功能划分,城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等)
(2) 通过数据分析,说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征,由此建立模型,确定污染源的位置。
(4) 分析你所建立模型的优缺点,为更好地研究城市地质环境的演变模式,还应收集什么信息?有了这些信息,如何建立模型解决问题?
二、求解思路概括
主要的其实就是前三问。拆分出来,任务的要求也很具体。下面结合论文【2】的内容简要介绍。
1、表示空间分布
作图题。然后根据图上的直观表现定性分析污染原因、污染程度、污染分布等。
2、分析原因
原因的话,其实就是一些功能区的特点造成的。所以文中用了单因素方差分析来说明不同功能区之间金属浓度的差异有显著性。为了说明不同区的污染程度,计算了尼梅罗污染指数。
3、污染源位置
简单粗暴直接假设污染源为重金属含量最高的地方。这一问变成了求重金属污染的完整的空间分布。用到的方法是半方差函数和克里格插值。最后还分析了一下各种元素以及海拔的相关性。
三、涉及到的知识技能的说明
1、数据预处理
数据超出 A ± 3 S A\pm3S A±3S 范围的为异常数据,当采用数据大于 A + 3 S A+3S A+3S时,以 A + 3 S A+3S A+3S代替;当采用数据小于 A − 3 S A-3S A−3S时,以 A − 3 S A-3S A−3S代替。
2、单因素方差分析【3】
要考察的因素称为试验指标(这里是某一种金属的含量),影响实验指标的因素称因素(这里是功能区),因素所处的不同状态称为水平(这里有5个水平)。考察单个因素对实验指标的影响情况,用单因素方差分析。
五个总体(五个功能区)的金属含量均值均值如果分别为
μ
1
μ_{1}
μ1、
μ
2
μ_{2}
μ2、
μ
3
μ_{3}
μ3、
μ
4
μ_{4}
μ4、
μ
5
μ_{5}
μ5
则有假设检验:
H
0
:
μ
1
=
μ
2
=
μ
3
=
μ
4
=
μ
5
H0:μ_{1}=μ_{2}=μ_{3}=μ_{4}=μ_{5}
H0:μ1=μ2=μ3=μ4=μ5
H
1
:
μ
1
、
μ
2
、
μ
3
、
μ
4
、
μ
5
不
全
相
等
H1:μ1、μ2、μ3、μ4、μ5不全相等
H1:μ1、μ2、μ3、μ4、μ5不全相等
于是可以构造相关检验统计量分析。
该检验统计量为
F
=
S
A
/
(
s
−
1
)
S
E
/
(
n
−
s
)
F=\frac {S_{A}/(s-1)}{S_{E}/(n-s)}
F=SE/(n−s)SA/(s−1)
其中
S
A
=
∑
j
=
1
s
∑
i
=
1
n
j
(
X
.
j
−
X
ˉ
)
2
S_{A}=\sum_{j=1}^s\sum_{i=1}^{n_{j}}(X_{.j}- \bar X)^2
SA=j=1∑si=1∑nj(X.j−Xˉ)2
S
E
=
∑
j
=
1
s
∑
i
=
1
n
j
(
X
i
j
−
X
.
j
ˉ
)
2
S_{E}=\sum_{j=1}^s\sum_{i=1}^{n_{j}}(X_{ij}- \bar {X_{.j}})^2
SE=j=1∑si=1∑nj(Xij−X.jˉ)2
S
T
=
S
A
+
S
E
S_T=S_A+S_E
ST=SA+SE
s
s
s是水平数,
n
n
n是样本数。
S
A
S_A
SA称效应平方和,衡量每个水平下样本与全体样本均值的差异,
S
E
S_E
SE称误差平方和,衡量同一水平下,每个样本与该水平样本均值的差异。
拒绝域为
F
=
S
A
/
(
s
−
1
)
S
E
/
(
n
−
s
)
≥
F
α
(
s
−
1
,
n
−
s
)
F=\frac {S_{A}/(s-1)}{S_{E}/(n-s)}\geq F_α(s-1,n-s)
F=SE/(n−s)SA/(s−1)≥Fα(s−1,n−s)
计算结果可以排成方差分析表。
方差来源 | 平方和 | 自由度 | 均方 | F比 |
---|---|---|---|---|
因素A | S A S_A SA | n − 1 n-1 n−1 | S ˉ A = S A s − 1 \bar S_{A}=\frac {S_{A}}{s-1} SˉA=s−1SA | F = S ˉ A S ˉ E F=\frac {\bar S_{A}}{\bar S_{E}} F=SˉESˉA |
误差 | S E S_E SE | n − s n-s n−s | S ˉ E = S E n − s \bar S_{E}=\frac {S_{E}}{n-s} SˉE=n−sSE | |
总和 | S T S_T ST | n − 1 n-1 n−1 |
3、半方差函数和克里格插值
在最初尝试复现半方差函数参数的求解时,使用了粒子群算法,参考了《基于粒子群算法与最小二乘拟合函数参数》。
但是由于参数的复杂性(有一些参数是用来给半方差函数分段的),结果很大程度取决于最初为待定系数设置的允许的取值范围,而且还没有解决克里格插值的问题。最终决定使用ArcMap软件复现这一步求解。使用方法参见我的另一篇文章《从0开始使用ArcGIS实现克里金(Kriging)插值(适用数学建模竞赛,详解)》。
这里对克里格插值的原理再做一个简要的介绍。
插值的原理如下述公式:
Z
V
∗
=
∑
i
=
1
n
λ
i
Z
(
x
i
)
Z_V^*=\sum_{i=1}^nλ_iZ(x_i)
ZV∗=i=1∑nλiZ(xi)
简单的说,就是:要想确定某个点的值,用周围的几个值
Z
(
x
i
)
Z(x_i)
Z(xi)与权重
λ
i
λ_i
λi相乘后相加,确定出该点的值。例如,最简单的权值确定办法可以用距离的倒数归一化后确定。
而克里格插值确定权重的方法要保证两点:一是要使估计无偏,二是要使估计值和实际值的差的平方和最小。即:
E
[
Z
V
∗
(
x
)
−
Z
V
(
x
)
]
=
0
E[Z_V^*(x)-Z_V(x)]=0
E[ZV∗(x)−ZV(x)]=0
V
a
r
[
Z
V
∗
(
x
)
−
Z
V
(
x
)
]
=
E
[
(
Z
V
∗
(
x
)
−
Z
V
(
x
)
)
2
]
→
m
i
n
Var[Z_V^*(x)-Z_V(x)]=E[(Z_V^*(x)-Z_V(x))^2]\rightarrow min
Var[ZV∗(x)−ZV(x)]=E[(ZV∗(x)−ZV(x))2]→min
λ
i
λ_i
λi的选取必须要保证上面两个条件,具体方法,需要用到拉格朗日乘数原理,这里不深入探讨。(可参考《克里格插值法》)
部分参考资料
(文中给出链接的,这里不再重复标注)
【1】2011高教社杯全国大学生数学建模竞赛题目:A题 城市表层土壤重金属污染分析
【2】王耀,李欣,王子嘉《城市表层土壤重金属污染分析》
【3】盛骤,谢式千,潘承毅《概率论与数理统计(第四版)》高等教育出版社