2011年数学建模A题答案

1 问题重述

随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。

现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。

附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。

现要求你们通过数学建模来完成以下任务：

(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。

(2) 通过数据分析，说明重金属污染的主要原因。

(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？

2 模型假设

（1）污染物的浓度在短期内不会有大的变动。

（2）取样点的8种重金属浓度数据能较好的反映该区的污染情况。 （3）给出的样本的采样点的位置、海拔高度测量准备。

3 符号说明

pi 土壤中污染物i的环境质量指数

ci 物i的实测浓度

1

si 污染物i的评价标准

xi 第i种指标变量

xij 第i个评价对象的第j个指标的取值

rij 第i个评价对象的第j个指标的相关系数

bj 主成分yj的信息贡献率 ap 主成分yp的累积贡献率

Z 海拔

s 污染物的扩散速度

t 任意扩散时刻

Q0 从污染源放出的污染物总量 Q1 在t,tt内通过某空间域的流量

Q2 在某空间域污染物的增量

p0 为单位体积内污染物排放的速率

4 问题分析

为了正确说明及评价城市表层土壤重金属污染程度，本文按照题目要求，通过四个问题来说明城市表层土壤重金属污染程度。

问题一，给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度，本部分的工作从以下几个方面进行：

1．采用地统计学软件ArcGIS，建立8种主要重金属元素在该城区的空间分布图。

2.从给出的8种主要重金属元素在该城区的空间分布图中，直观的分析城区内不同区域重金属的污染程度。

3. 为了更客观、准确地对不同区域重金属的污染程度分析，采用目前国内外普遍采用的单因子指数法和内梅罗综合指数法建立模型，对城区的不同区域重金属污染程度进行分析。

4.通过比较得出不同区的重金属不同的污染指数，与背景值作比较，得出生活区、工业区、山区、主干道路区及公园绿地区中8种主要重金属元素的金属污染程度比较。

第二问题，通过数据分析，说明重金属污染的主要原因。本部分的工作从以

2

下几个方面进行：

1. 用SPSS软件对8种主要重金属元素进行相关分析，得出哪些金属具有强的相关性，从而得出具有相关性的金属我们认为它们是来来源一样。

2. 用matlab软件对8种主要重金属元素进行主成分分析，主成分分析方法是一种将因子纳入同一系统进行定量化研究、理论成熟的多元统计分析方法。通过分析变量之间的相关性，使得所反映信息重叠的变量被某一主成分替代，减少变量的数目，从而降低系统评价的复杂性。再以方差贡献率作为每个主成份的权重，由每个主成份的得分加权确定重金属污染的综合评价。

问题三，分析重金属污染物的传播特征，由此建立模型，确定污染源的位置. 本部分的工作从以下几个方面进行：

1.根据题目和有关资料对重金属污染物的传播特征的分析。 2.建立克里格插值法模型来确立污染源的位置。

问题四，分析模型的优缺点，收集相应信息，建立扩散模型 1、分析问题三建立模型的优缺点，以此确定需要收集的信息。 2、考虑影响重金属的污染因素，分类讨论建立扩散模型。

5 模型建立与求解

5.1重金属元素的空间分布

本文采用地统计学软件GIS，建立8种主要重金属元素在该城区的空间分布图。其中图一为功能区空间分布图，图二至图九为8种主要重金属元素的空间分布图。根据8种主要重金属元素用GIS建立空间分布图, 可以对每种金属的污染程度进行可视化表达, 直观显示区域污染情况的分布变化。

图一 功能区空间分布图

1粉色生活区 2 红色工厂 3 橘黄山区 4 蓝色主干道路 5绿色公园绿地

3

图二 As空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为As空间分布图（图二），其颜色的深浅代表As浓度，颜色越深说明该区域的As浓度越高，我们可以直观地看出As浓度较高是工厂、生活区、主干道路。

图三 Cd元素空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为Cd空间分布图（图三），其颜色的深浅代表Cd浓度，颜色越深说明该区域的Cd浓度越高，我们可以直观地看出Cd浓度较高是工厂、公园绿地。

4

图四 Cr元素空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为Cr空间分布图（图四），其颜色的深浅代表Cr浓度，颜色越深说明该区域的Cr浓度越高，我们可以直观地看出Cr浓度较高是工厂、山区、公园绿地。

图五 Cu元素空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为Cu空间分布图（图五），其颜色的深浅代表Cu浓度，颜色越深说明该区域的Cu浓度越高，我们可以直观地看出Cu浓度较高是工厂、山区、公园绿地。

5

图六 Hg元素空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为Hg空间分布图（图六），其颜色的深浅代表Hg浓度，颜色越深说明该区域的Hg浓度越高，我们可以直观地看出Hg浓度较高是工厂、公园绿地。

图七 Ni元素空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为Ni空间分布图（图七），其颜色的深浅代表Ni浓度，颜色越深说明该区域的Ni浓度越高，我们可以直观地看出Ni浓度较高是工厂、山区、公园绿地。

6

图八 Pb元素空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为Pb空间分布图（图八），其颜色的深浅代表Pb浓度，颜色越深说明该区域的Pb浓度越高，我们可以直观地看出Pb浓度较高是工厂、山区、公园绿地。

图九 Zn元素空间分布图

1生活区 2 工厂 3 山区 4 主干道路 5公园绿地

上图为Zn空间分布图（图九），其颜色的深浅代表Zn浓度，颜色越深说明该区域的Zn浓度越高，我们可以直观地看出Zn浓度较高是生活区、工厂、公园绿地。

5.2该城区内不同区域重金属的污染程度的分析[1][2][3]

7

5.2.1土壤污染的评价方法[11]

本文的土壤污染评价标准分别以《国家土壤环境质量标准》二级标准[8]以及按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值为参照标准，见表1。

表1 8种主要重金属元素的背景值 元素 As (ug/g) 平均值 3.6 130 31 13.2 35 12.3 31 69 标准偏差 0.9 30 9 3.6 8 3.8 6 14 范围 1.8~5.4 70~190 13~49 6.0~20.4 19~51 4.7~19.9 19~43 41~97 Cd(ug/g) Cr (ug/g) Cu (ug/g) Hg (ug/g) Ni (ug/g) Pb(ug/g) Zn (ug/g)

土壤是人类赖以生存的最基本的自然资源之一，也是生物可利用重金属的一个重要蓄积库，其所含的重金属通过食物链被植物、动物数十倍的富集,通过多种途径直接或间接地威胁人类安全和健康[7][8]。随着工业、城市污染的加剧和农用化学物质种类、数量的增加,土壤重金属污染日益严重。这种形势下迫切需要对所处的土壤环境质量做出客观、切实的综合评价,以此反映经济、技术发展对土壤质量、农业生产、生态环境乃至人类健康的影响,并为土地的可持续利用提供理论依据。在我国当前大规模农业地质环境调查活动中,有关土壤重金属污染状况评价工作业已列入议事日程,相关学科学者正在积极探讨简便有效实用的评价方法。而对于土壤重金属污染评价的方法采用目前国内外普遍采用的单因子指数法和内梅罗综合指数法。 a·单因子指数法

pici/si其中

pi为土壤中污染物i的环境质量指数; ci为物

i的实测浓度(mg kg-1);

si为污染物i的评价标准(mg kg-1)

b·内梅罗综合指数法

p1其中nn1ni1pinmaxp2i22

i1pi为土壤中各污染指数平均值; maxpi为土壤中各污染指数最大

8

值。

根据pi值变幅，基于污染指数主要划分如表2所示土壤质量分级 表2 土壤质量等级划分 指数范围 pi1 1表3该城市城区土壤重金属的含量

金属 最小值 最大值 As Cr Cu Hg Ni 1.61 15.32 2.29 8.57 4.27 30.13 920.84 2528.48 16000.00 142.5 472.48 3760.82 最大值所在区域 平均值 标准差 交通区 交通区 交通区 工业区 交通区 交通区 生活区 交通区 5.68 53.51 55.02 17.26 61.74 3.01 69.78 162.41 9.91 49.90 302.40 224.28 变异系数 52.99% 74.17% 130.41% 295.18% 57.42% 80.82% 168.08% Cd 40.00 1619.8 299.71 1624.44 542.00% Pb 19.68 Zn 32.86 201.20 338.17 表3数据显示，研究区土壤样品的8种主要重金属的平均含量分别为，As 5.68(ug/g)，Cd 302.40(ug/g)，Cr 53.51(ug/g)，Cu 55.02(ug/g)，Hg 299.71(ug/g)，Ni 17.26(ug/g)，Pb 61.74(ug/g)，Zn 201.20(ug/g)。变异系数反映了总体样本中各采样点的平均变异程度，在该城市的区域汞的变异系数最大，达到542.00%；而As得变异系数最小，为52.99%。八种重金属的平均变异程度由大到小的顺序为Hg，Cu，Zn，Cr，Pb，Cd，Ni ，As。Pb主要分布在生活区，Cu元素主要分布在工业区，其他元素主要分布在交通区。 5.2.3重金属的污染指数值

以相关元素背景值为参照标准，计算了研究区319个样点表层土壤各种金属的单因子污染指数值，数据列入表4.研究区土壤重金属单因子平均指数值：As ，Cd ，Cr， Cu， Hg ，Ni ，Pb， Zn分别为1.58、2.33、1.73、4.17、8.56、1.40、1.99、2.92，其中As，Cr，Ni ，Pb大于1小于2，属于轻微污染，但是最大值都超过2，说明该研究区部分地区出现As，Cr，Ni ，Pb明显外源污染现象。Cd和Zn超过2小于3，属于中度污染，而Cu，和Hg都大于3，属于严重污染。Cu和Hg为主要污染因子，重金属单因子平均指数值从大到小排列顺序为：

Hg>Cu >Zn> Cd>Pb > Cr> As >Ni.。从土壤样品中不同重金属含量的变异系数来看，都大于50.%说明该城市的样点重金属指数的平均变异程度普遍很高。

表4该城市地区重金属单因子指数统计分析 金属 最小值 最大值 平均值 标准差 As Cd Cr

变异系数 53.16% 74.25% 130.06% 0.45 0.31 0.49 8.37 12.46 29.7 1.58 2.33 1.73 9

0.84 1.73 2.25

Cu Hg Ni 0.17 0.24 0.35 0.63 0.48 191.55 457.14 11.59 15.24 54.5 4.17 8.56 1.40 1.99 2.92 12.32 46.49 0.81 1.61 4.91 295.44% 543.11% 57.86% 80.90% 168.15% Pb Zn 根据所给的数据将已区分好的各个区域归在一起，求出各种重金属元素在该区域的平均值，建立综合污染指数评价法模型；分析各种重金属元素在各个区域的污染指数来分析该城区内不同区域重金属的污染程度。由excel可以计算得到各种重金属元素在各个区域的污染指数。综合各区域各元素的污染指数（见附录一），得到下表5

表5不同功能区的土壤重金属的污染指数 区域 生活区 As Cd Cr Cu Hg Ni Pb Zn P 1.74 2.23 2.23 3.74 2.66 1.49 2.23 3.43 3.17 工业区 2.01 3.02 1.72 9.66 18.35 1.61 3.00 4.03 13.53 山区 1.12 1.17 1.26 1.31 1.17 1.26 1.18 1.06 1.25 主干道路区 1.59 2.77 1.87 4.71 12.77 1.43 2.05 3.52 9.43 公园绿地区 1.74 2.16 1.41 2.29 3.29 1.24 1.96 2.24 2.73 通过以上表格分析可知，生活区为严重污染，工业区为严重污染，山区为轻度污染，主干道路区为严重污染，公园绿地区为中度污染。 5.3重金属污染的主要原因 5.3.1土壤重金属的相关分析[5]

表6计算了土壤中重金属间的相关系数，从中可以看出，海拔与As、Cd、Cr、Cu、Ni、Pb、Zn之间有较好的相关性，说明海拔可能会影响元素的分布。As、Cd、Cr、Cu、Hg、Ni、Pb、Zn之间有较好的相关性，表明这些元素可能来源于统一土母质。As与Cd、Cr、Cu、Ni、Pb、Zn，Cd与Zn、Cu、Ni、Pb、Cr，Cu 与Ni、Pb、Zn，Ni与Pb、Zn， Pb与Zn之间存在显著的相关管，由此可初步推断，As、Cd、Cr、Cu、Ni、Pb、Zn七种元素可能同一来源，而Hg为其他来源。

表6海拔和土壤重金属元素间的相关系数

10

Correlations海拔海拔Pearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)NPearson CorrelationSig. (2-tailed)N1319-.2**.000319-.248**.000319-.152**.007319-.138*.014319-.084.134319-.163**.004319-.235**.000319-.178**.001319As-.2**.0003191319.255**.000319.1**.001319.160**.004319.0.251319.317**.000319.290**.000319.247**.000319Cd-.248**.000319.255**.0003191319.352**.000319.397**.000319.265**.000319.329**.000319.660**.000319.431**.000319Cr-.152**.007319.1**.001319.352**.0003191319.532**.000319.103.066319.716**.000319.383**.000319.424**.000319Cu-.138*.014319.160**.004319.397**.000319.532**.0003191319.417**.000319.495**.000319.520**.000319.387**.000319Hg-.084.134319.0.251319.265**.000319.103.066319.417**.0003191319.103.066319.298**.000319.196**.000319Ni-.163**.004319.317**.000319.329**.000319.716**.000319.495**.000319.103.0663191319.307**.000319.436**.000319Pb-.235**.000319.290**.000319.660**.000319.383**.000319.520**.000319.298**.000319.307**.0003191319.494**.000319Zn-.178**.001319.247**.000319.431**.000319.424**.000319.387**.000319.196**.000319.436**.000319.494**.0003191319AsCdCrCuHgNiPbZn**. Correlation is significant at the 0.01 level (2-tailed).*. Correlation is significant at the 0.05 level (2-tailed). 5.3.2土壤重金属的主成分分析[5]

根据所求出每个区域的各种金属的平均值，采用主成分分析法对其对进行综合评价分析。

5.3.2.1 模型的建立

1）对原数据进行标准化处理 假设进行主成分分析的指标变量为：x1,x2,,x，第i个评价对象的第j个指标

的取值为xij。将各个指标xij转换成标准化指标xij， xijxijxjdj， (i1,2,,)

1n1n其中xjxij,djxijxj,j1,2,ni1n1i12,，即xj,dj为第j个指标的样

本均值和样本标准差。对应的，称xij量。

2）计算相关系数矩阵R 相关系数矩阵Rrij

xijxjdj， (i1,2,,)为标准化指标变

11

,i,j1,2,, n1式中rij1,rijrji,rij是第i个指标与第j个指标的相关系数。rijk1xnkixkj3）计算特征值和特征向量 计算相关系数矩阵R的特征值120，及对应的特征向量

T1,2,,，其中j1j,2j,,n,由特征向量组成m个新的指标变量

y111x121x2n1xny12x122x2n2xn 2

yxxx1m12m2nmnm式中y1是第一主成分，y2是第二主成分，，ym是第m主成分。

4）选择p(pm)个主成分，计算综合评价值 a、计算特征值jj1,2, bj,m的信息贡献率和累积贡献率。称

jk1mj1,2,k,m

为主成分yj的信息贡献率；

pk1k1mpk

k为主成分y1,y2,,yp的累积贡献率，当p接近于1（p0.85,0.90,0.95）时，

,yp作为p个主成分，代替原来m个指标变量，

则选择前p个指标变量y1,y2,从而可对p个主成分进行综合分析。 b、计算综合得分：

bjyj

j1p其中bj为第j个主成分的信息贡献率，根据综合得分就可进行评价。

12

5.3.2.2 模型的求解

经过对土壤单点样重金属元素含量的数据标准化处理后，利用matlab软件进行主成分分析，可得出如下结果；

表7 变量相关矩阵 指标 As Cd Cr Cu Hg Ni Pb Zn As 1.000 .739 -.296 .473 .319 .100 .2 .590 Cd .739 1.000 -.546 .731 .762 .225 .830 .6 Cr -.296 -.546 1.000 -.780 -.842 -.015 -.523 .242 Cu Hg Ni Pb Zn .473 .731 -.780 1.000 .943 .595 .809 .292 .319 .762 -.842 .943 1.000 .423 .656 .205 .100 .225 -.015 .595 .423 1.000 .442 .501 .2 .830 -.523 .809 .656 .442 1.000 .585 .590 .6 .242 .292 .205 .501 .585 1.000 再利用SPSS统计软件可得到前几个因子的累积贡献率和初始特征值，如表8

所示；

表8 解释的总方差 初始特征值 因子 合计 1 2 3 4 5 6 7 8 4.828 1.676 1.109 0.368 3.913E-8 2.3E-9 2.194E-9 6.468E-18 提取平方和载入 方差百分比% 累积贡献值% 合计 方差百分比% 累积贡献值% 60.355 20.956 13.862 4.828 4.1E-7 3.304E-8 2.743E-8 8.085E-17 60.355 81.310 95.172 100.000 100.000 100.000 100.000 100.000 4.828 1.676 1.109 60.355 20.956 13.862 60.355 81.310 95.172 由表8可见，前三个主成分分别累计提取了总方差的95.172%，表明所提取的主成分能够较好地代表源数据所蕴涵的信息。将表8的前三个主成分的主要元素因子和土壤中的重金属元素的载荷矩阵反映在表9中；

表9 载荷矩阵 元素 1 As 2 因子 3 .730 .368 -.508 Cd .909 .104 -.245 Cr -.660 .736 .146 Cu Hg .928 -.272 .245 .861 -.422 .187 13

Ni Pb Zn .495 .274 .806 .942 .176 -.144 .549 .794 .070 从表9中可见，Cd、Cu和Pb在主成分1中具有非常高的载荷，Hg和As在主成分1中的载荷也比较高。研究证实Pb和Cd在来源上关系较密切，在城市土壤中Pb常被作为机动车污染源的标识元素，Pb主要来源于汽车尾气排放而汽车轮胎磨损可以释放Cd，再加上Cu可以来源于汽油和车体的磨损，从而可以推断出因子1可能主要代表汽车交通运输活动的影响。密集的工业活动是Cd和Hg的重要来源之一,其次Hg污染不仅与燃煤污染有关，在很大的程度上还与汽车尾气的排放有关；As和S一样在空气中以AsO2为主，也可以是由汽车燃料燃烧所释放的，进一步的加大了推断的可能性。由表格可以看出，第一主成分解释的原变量方差最多（60.355%），而其后的主成分所解释的方差逐渐在减少，因此可以认为第一主成分是污染的主要来源，鉴于此，可认为主干道路区上的尾气排放和工业活动是重金属Cd、Cu、Pb、Hg和As污染的主要原因。

在表9中体现了，Cr和Zn在主成分2中具有较高的载荷。有学者研究表明，厂矿企业的三废排放中含有大量的Zn和Cr，因此可以断定因子2为工业活动的影响。与此我们可以认为工业区的工业废气排放是重金属Zn和Cr的主要原因。

在表9中也体现了，Ni单独在主成分3中具有较高的载荷，且与其他元素关系不密切。由资料显示，主要受地球化学成因影响。废水中含有较多的Ni，因此也可以认为生活区的生活废水、污水是重金属Ni的主要原因。

综合以上所述，可知Ni使生活区的污染变严重，Cr和Zn使工业区的污染变严重，Cd、Cu、Pb、Hg和As使主干道路区的污染变严重。与第一问的结果一致，验证了两个模型的准确性。 5.4 寻找污染源

5.4.1克里格法（Kriging）模型的简介[6][7][8]

克里格法（Kriging）是地统计学的主要内容之一，从统计意义上说，是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲是对空间分布的数据求线性最优、无偏内插估计一种方法。克里格法的适用条件是区域化变量存在空间相关性。 克里格法，基本包括普通克里格方法（对点估计的点克里格法和对块估计的块段克里格法）、泛克里格法、协同克里格法、对数正态克里格法、指示克里格法、折取克里格法等等。随着克里格法与其它学科的渗透，形成了一些边缘学科，发展了一些新的克里金方法。如与分形的结合，发展了分形克里金法；与三角函数的结合，发展了三角克里金法；与模糊理论的结合，发展了模糊克里金法等等。 应用克里格法首先要明确三个重要的概念。一是区域化变量；二是协方差函数，三是变异函数 一区域化变量

当一个变量呈空间分布时，就称之为区域化变量。这种变量反映了空间某种属性的分布特征。矿产、地质、海洋、土壤、气象、水文、生态、温度、浓度等领域都具有某种空间属性。区域化变量具有双重性，在观测前区域化变量Zx是一个随机场，观测后是一个确定的空间点函数值。

14

区域化变量具有两个重要的特征。一是区域化变量Zx是一个随机函数，它具有局部的、随机的、异常的特征；其次是区域化变量具有一般的或平均的结构性质，即变量在x点与偏离空间距离为h的点xh处的随机量Zx与

Zxh具有某种程度的自相关，而且这种自相关性依赖于两点间的距离h与变量特征。在某种意义上说这就是区域化变量的结构性特征。 二、协方差函数[9] 协方差又称半方差，是用来描述区域化随机变量之间的差异的参数。在概率理论中，随机向量X与Y的协方差被定义为：

Cov(X,Y)EXEXYEY

区域化变量ZxZxu,xv,xw在空间点 x和xh处的两个随机变量Zx和 Zxh 的二阶混合中心矩定义为Zx的自协方差函数，即

CovZx,ZxhEZxZxhEZxEZxh

三、变异函数[10]

变异函数又称变差函数、变异矩，是地统计分析所特有的基本工具。在一维条件下变异函数定义为，当空间点x在一维x 轴上变化时，区域化变量 Zx在点 x和xh 处的值Zx 与Zxh 差的方差的一半为区域化变量Zx在

x轴方

1向上的变异函数，记为

h，即

x,hVarZxZxh2211EZxZxhEZxEZxh222

在二阶平稳假设条件下，对任意的h有， EZxhEZx 因此上式可以改写为：

x,hEZxZxh 2从上式可知，变异函数依赖于两个自变量x和 h，当变异函数x,h仅仅依赖于距离h 而与位置x无关时， 可改写成h，即

12 15

hEZxZxh 2设 Zx是系统某属性 Z在空间位置x处的值， Zx为一区域化随机变量，并满足二阶平稳假设，h为两样本点空间分隔距离，Zxi和 Zxih分别是区域化变量 在空间位置xi 和xih 处的实测值i1,2,据上式的定义，变异函数h的离散公式为：

1NhhZxiZxih 2Nhi1212Nh ，那么根

变异函数揭示了在整个尺度上的空间变异格局，而且变异函数只有在最大间

1隔距离处才有意义。

25.4.2 基于ArcGIS系统的实证分析[6][8][10]

克里格插值法是对区域变量进行无偏的最佳估值的可靠方法，它假设数据是服从正态分布的。ArcGIS中普通克里格插值包括4个部分功能：创建预测图、创建分位数图、创建概率图、创建标准误差的预测图。我们创建预测图，直观看出8种重金属污染源的区域部分和污染源的个数，并分别求出其相应污染源的x, y，z,以及浓度。 5.4.2.1选择一个参考点

首先我们先任意选择一个参考点（图十），以这个点为参考点把该城区划分为4个区，再在附近的4个区中各选20个点，把这个参考点与周边的点的浓度作比较。我们从交叉检验对话框（如图十一），我们给出这个点模型精确度评价。对不同的参数得到模型的比较中，可参考Prediction Error中的几个指标，符合以下标准的模型是最优的：标准平均值（Mean Standardized）的值为-0.03285，其值近于0，均方根预测误差（Root-Mean-Square）的值为208.1最小，平均标准误差的值263.9（Average Mean Error）的值最接近于均方根预测误差（Root-Mean-Square）的值208.1，标准均方根（Root-Mean-Square Standardized）为0.924，很接近1.所以我们选择这个点是优的。

16

图十 协方差模型对话框

图十一交叉检验对话框

5.4.2.2 污染源的预测

As为克里格内插预测生成图（图 十二），我们直观看出污染源有三个，我们可以与第一问8种金属元素的空间分布图直观看出污染源的区域情况，AS元素的污染源分布在工厂,我们由ArcGIS得出污染物的源头一的位置为x=2561.459m,y=8433.325m,z=5.45m,此污染源的浓度为9.15668ug/g,源头二的位置为x=8635.3m ,y=147.1m ,z=16.6m,此污染源的浓度为8.884ug/g，源头三的位置为X=18333.920m,y=10066.719m,z=45m,此污染源的浓度为6.249ug/g。

Cd为克里格内插预测生成图（图十三），我们直观看出污染源有两个，Cd

17

元素的污染源分布在工厂、公园绿地,我们由ArcGIS得出污染物的源头一的位置为x=8584.599m,y=7004.105m,z=39m,此污染源的浓度为587.234ug/g，源头二的位置为x=1540.588m,y=3737.317m,z=6.3m,此污染源的浓度为556.802ug。 Cr为克里格内插预测生成图（图十四），我们直观看出污染源有一个，我们由ArcGIS得出污染物的源头一的位置为x=2816.677m,y=5370.711,z=4.2m, 此污染源的浓度为145.109ug/g。

Cu为克里格内插预测生成图（图十五），我们直观有两个污染源有两个，我们由ArcGIS得出污染物的源头一的位置为x=9146.078m,y=83832.281m,z=38.7m,此污染源的浓度为

105.14ug/g,源头二的位置为

x=2255.197m,y=4451.927m,z=4.3m,此污染源的浓度为152..215ug/g。 Hg为克里格内插预测生成图（图十六），我们直观看出污染源有一个，我们由ArcGIS得出污染物的源头一的位置为x=1183.283m,y=3175.838m,z=8.1m, 此污染源的浓度为817.425ug/g。

Ni为克里格内插预测生成图（图十七），我们直观看出有一个污染源，我们由ArcGIS得出污染物的源头一的位置为x=1999.980m,y=5625.929m,z=2.65m, 此污染源的浓度为27.4574ug/g。

Pb为克里格内插预测生成图（图十八），我们直观看出污染源有一个，我们由ArcGIS得出污染物的源头一的位置为x=4292.155m,y=5133.039m,z=6.8m, 此污染源的浓度为144.109ug/g。

Zn为克里格内插预测生成图（图十九），我们直观看出污染源有两个，我们由ArcGIS得出污染物的源头的位置为x=8110.851m,5113.039m,z=18.8m,此污染源的浓度为397.344ug/g,源头二的位置为X=2827.842m,y=53.792m,z=3.6m,此污染源的浓度为545.8ug/g。

18

图十二 AS克里格内插预测生成图 图

图十四 Cr克里格内插预测生成图 图

19

图十三 Cd克里格内插预测生成 Cu克里格内插预测生成

图十五

图图十六 Hg克里格内插预测生成图

图十七 Ni克里格内插预测生成图

20

图十八 Pb克里格内插预测生成图 图十九 Zn克里格内插预测生成图

5.5模型的改进

5.5.1模型的优缺点

本文采用的克里格法（Kriging）模型，是克里格插值法寻找污染源，插值精确度高，这种方法要比因子分子法、聚类分析法寻找污染源更加优越而且比较具体、明确的求出所求点，相对来说所求值比较精确、清晰。图形结合，简化思维，感官效果比较好。

但是考虑的因素不够全面，未考虑其他污染因素对模型的影响。因为还有诸多因素影响着土壤重金属的传播，那是一个动态的、灵活的改变。所以对于污染源的的精确寻找还是有一定的困难的。

为更好地研究城市地质环境的演变模式，我们还需收集其他信息.如：土壤的PH值，生物降解量，城区人们生活和工业的垃圾排放量，降雨量对重金属扩散的影响，风中的传播速度等来描述各区域重金属污染的情况及它的扩展发散的形式，以便更好的了解重金属的污染趋势，做好防护准备。

具体的可将8种重金属污染物的传播方式分为两类。一类是主要在大气中传播，一类是主要在土壤中传播。对于在大气中传播的重金属污染物，我们建立重金属污染物在气体中扩散模型来确立污染源的位置；而对于另一类，我们则建立重金属污染物在土壤中扩散模型来确立污染源的位置。 5.5.2建立新模型

1）重金属污染物在气体中的扩散模型[12][13][14]

记污染物开始扩散的时刻为t0，并且以污染源作为扩散中心，则可将t 时刻无穷空间中任一点(x,y,z)的污染物浓度记为C(x,y,z,t)。假设单位时间内通过单位法向量面积的流量与浓度梯度成正比，则有

qkigradC

其中kiix,y,z是扩散系数，grad表示浓度梯度，其中的负号代表重金属粉尘的浓度是由高到低的地方扩散。

21

假设空间域的体积为 V，包围空间域的曲面为一规则的球面，设其表

xy面面积为S表,外法线向量为 n,,1，

zz则在t,tt内流通过空间域的流量可表示为： Q1tttqnddt

s空间域所包围的区域内污染的重金属的增量可表示为：

Q2Cx,y,z,ttCx,y,z,tdV

V由污染源排放出的该重金属的总量为：Q0由质量守恒可得出：Q0Q1Q2,即，

tttpdVdt

0C化简得：dVdivqdVP0

tVV根据A. FIck扩散微分方程式中：

CCCC2V2V2VUXUYUZkXkYkZ 22tXYZXYZ2其中：C为气体浓度；t为时间；ux,uy,uy为x,y,z方向风速；kx,ky,kz为x,y,z 方向上的扩散系数。

由此得到地面土壤表层的浓度：

C'x,y,z,tvswdCx,y,z,t x2y2z2s2t2

2）重金属污染物在土壤中的扩散模型[15]

记污染物开始扩散的时刻为t0，并且以污染源作为扩散中心，则可将t 时刻无穷平面中任一点(x,y)的污染物浓度记为 C(x,y,z,t)。假设单位时间内通过单位路程的流量与浓度梯度成正比，则有

qkigradC

其中kiix,y是扩散系数，grad表示浓度梯度，其中的负号代表重金属污染物的浓度是由高到低的地方扩散。

22

x假设在平面S上，设其传播路径为l,法线向量为 n,1，

y则在t,tt内流通过平面S的流量可表示为： Q1tttqndsdt

l平面S所包围的区域面积内污染的重金属的增量可表示为： Q2Cx,y,ttCx,y,td

s在有水气共存的土壤中,即土壤非饱和带中,污染物随水运移时,水中的污染物会与气体及土壤颗粒发生物质交换。水气两相间的物质交换平衡可用亨利定律来表示: CaHCw (22) 式中，Ca为气体中的污染物浓度mgcm3,Cw为水中的污染物浓度mgcm3,

H为亨利系数。不同物质的H值可在文献中查到。对于水中污染物与固体颗粒间的物质交换,经常用到以下这一简单的线性方程:

qsKdCw 式中，qs为固体颗粒表面的污染物浓度mgg1,Kd为污染物在水与固体颗粒间的分配系数cm3g1。

对于土壤中的水-固物质交换，可以使用Freundlich方程。它是一个经常被使用的非线性吸附等温线方程,其表达式为：

1/n qsKFCw

KF为Freundlich分配系数mg11/ng1cm3/n，1/n表示水-固吸附非线性

指数。

土壤中污染物的总浓度：Q2CwCaqs

由污染源排放出的该重金属的总量为：Q0由质量守恒可得出：Q0Q1Q2,即，

Cx,y,ttCx,y,tdqndsdtpdsdt stlt0ltttttttpdsdt

0l根据曲面积分Gauss公式得：qndsldivqd

s 23

C化简得：ddivqdP0

tss

6 模型验证

本文得到的模型都是通过统计软件SPSS和ArcGIS挑选出的最优模型，经过

检验模型显著程度很高，所得结果有较高的可靠性。

7 模型评价

7.1模型优点

1）本文针对各个问题分别采用了单因子指数法、内梅罗综合指数法、主成分分析法、克罗格插值法等方法，建立数学模型清晰合理。

2）运用spss13.0软件对数据进行处理和运算，降低运算量，简单易行，有很大的可操作性，而且所得的数据较为合理可靠。

3）运用ArcGIS软件画出各种重金属含量的在各个功能区的空间分布图，直观明了。

4）建立的ArcGIS软件寻找污染源，还能找到其相应污染源的浓度。 7.2模型缺点及改进

任何模型、系统都受到实际生活中的各种，本文的模型也不例外，为了简化模型，基本假设很多都是理想状态。所以际运用中还应该考虑污染因素对模型的影响，在应用的过程中应该根据实际情况进行灵活改变。

参考文献

[1] 蔡立梅，马瑾，周永章，黄兰椿等东莞市农业土壤重金属的空间分布特征及来源解析[J],环境科学，2008.12. [2] 郑海龙，陈杰，邓文靖等,城市边缘带土壤重金属空间变异及其污染评价[J],土壤学报，2006.1.

[3] 刘庆，王静，张建国等,基于GIS的农田土壤重金属空间分布研究[J],安全与环境学报[J]，2007年第2期.

[4]张乃明,土壤-植物系统重金属污染研究现状与展望[J],环境科学进展,1999,7(4):30-33.

[5] 杨善朝，张军舰,SPSS统计软件应用基础[M] ,桂林：广西师范大学出版社,2000.2.

[6] 苏姝，林爱文，刘庆华，普通插值法在空间中的运用[J],江南大学学报，2004年2月，第三卷第一期. [7] 克里格插值法，http://wenku.baidu.com/view/c9a7e727a5e9856a5612609c.html.

[8] 唐翔宇,朱永官,土壤中重金属对人体生物有效性的体外试验评估[J],环境

24

与健康杂志,2004,21(3):183-185.

[9]崔德杰,张玉龙,土壤重金属污染现状与修复技术研究进展[J],土壤通报,2004,35(3):366-370. [10]仲嘉亮，基于RS和GIS技术的气温空间插值方法的探讨[J]，沙漠与绿洲气象，2007年8月。

[11]长江水质的评价和趋势分析模型，

http://wenku.baidu.com/veiw/c91b3342a56- -bec0975e3e6.html. [12] 姜启源，谢金星，叶俊，数学模型（第三版），高等教育出版社，2003,1-12. [13] 王沫然，MATLAB与科学计算（第二版），电子工业出版社，2004,-90. [14] 王冬琳，数学建模及实验，国防工业出版社，2004,197-199. [15] 放射气体模型的预估模型.

http://wenku.baidu.com/view/8f3aaed076eeaeaad -1f330cf.html

附录：一

生活区 金属 As Cd Cr Cu Hg Ni 最小值 最大值 平均值 标准差 变异系数 0.65 0.67 0.6 0.74 0.34 0.72 0.79 0.63 3.18 8.03 24.01 18.85 15.71 2.67 15.24 41.93 1.74 2.23 2.23 3.74 2.66 1.49 2.23 3.43 山区 金属 As Cd Cr Cu Hg Ni 最小值 0.49 0.31 0.52 0.17 0.28 0.45 0.63 0.48 最大值 3.05 3.14 5.59 5.23 5.91 6.02 3.67 3.33 平均值 标准差 变异系数 1.13 1.17 1.26 1.32 1.16 1.26 1.17 1.06 0.50 0.60 0.79 0.81 0.79 0.85 0.57 0.45 最小值 0.59 1.4 3.44 3.53 2.91 0.46 2.31 6.36 最小值 工业区 平均最大值 值 6.08 2.01 标准差 1.16 变异系数 57.71% 33.91% 0.45 62.78% 0.88 8.41 3.02 1.80 59.60% 154.26% 0.50 9.21 1.71 1.40 81.87% 94.39% 0.96 191.55 9.66 31.00 320.91% 109.40% 0.34 385.71 18.35 63.22 344.52% 30.87% 0.35 3.39 1.61 0.67 41.61% 103.59% 1.01 14.03 3.00 2.72 90.67% 185.42% 0.82 23.57 4.03 5.01 124.32% 交通区 平均最大值 值 标准差 变异系数 Pb Zn 44.25% 0.45 8.37 1.59 0.90 56.60% 51.28% 0.39 12.46 2.77 1.87 67.51% 62.70% 1.43 1.09 1.87 2.62 140.11% 61.36% 0.93 103.40 4.71 9.07 192.57% 68.10% 0.24 457.14 12.77 62.07 486.06% 67.46% 0.50 11.59 1.43 0.95 66.43% 48.72% 0.71 5.85 2.05 1.05 51.22% 42.45% 0.59 54.50 3.52 5.56 157.95% 变异系数 pb Zn 公园绿地区 金属 最小值 最大值 平均值 标准差 25

As Cd Cr Cu Hg Ni 0.65 0.67 0.6 0.74 0.34 0.72 0.79 0.63 3.18 8.03 24.01 18.85 15.71 2.67 15.24 41.93 1.74 2.23 2.23 3.74 2.66 1.49 2.23 3.43 0.59 1.4 3.44 3.53 2.91 0.46 2.31 6.36 33.91% 62.78% 154.26% 94.39% 109.40% 30.87% 103.59% 185.42% Pb Zn

26