2011全国赛论文 - 图文

2011高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 20111017 所属学校（请填写完整的全名）： 南京师范大学中北学院 参赛队员 (打印并签名) ：1.

陈锦 丁子荷

2. 3. 陶卫 指导教师或指导教师组负责人 (打印并签名)： 黄伯强

日期： 2011 年 9 月 12 日

赛区评阅编号（由赛区组委会评阅前进行编号）：

2011高教社杯全国大学生数学建模竞赛

编 号 专 用 页

评 阅 人 评 分 备 注 赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

城市表层土壤重金属污染分析

摘要

本文通过对某城市土壤的重金属污染分析，探讨人类活动对城市环境质量的影响，并研究在人类活动影响下城市地质环境的演变模式。

对于问题一，运用MATLAB软件及地理信息系统Arcgis软件，分别绘制出8种重金属元素在空间污染程度分布图。在分析不同区域重金属污染程度时，首先计算出某区域某一重金属的单项污染指数，接着求出该区域各重金属综合污染指数，对照国家环境质量标准，分析出不同区域重金属污染程度，得出结论：工业区、交通区中度污染，山区无污染。城区的综合污染指数由大到小依次为：交通区：1.1938；工业区：1.2912；生活区：1.0852；公园绿地区：0.9787；山区：0.8049。

在问题二中，对八种重金属元素进行相关性分析，将重金属元素分为四类。由问题二中各重金属污染指数分析得出:Pb、Cu、Zn主要来源于交通区，Cr、Ni主要来源于工业区，部分Cu、Zn来源于生活区。根据相关性较高的重金属元素在交通区、工业区、生活区、公园绿地区的污染分布图，验证了所得结论。

在问题三中，基于重金属污染物在土壤中以干湿沉降为主要传播方式，建立了重金属污染物在土壤中传播的数学模型，采用数值拟合的方法，将各重金属元素所有污染源对点的污染程度求和，即F(xi,yi,zi,?)，最后利用最小二乘法求出参数?i，根据以上模型，运用遍历搜索的方法找出污染源。计算得出全城区重金属元素As污染源坐标为（18134，,10046，41），（12696,3024,27），（4742,7293,9）（其余见表3-1）。 最后对问题三的模型进行优缺点分析，考虑到工业化进程加快、绿色植被增多等方面因素导致的地质环境的演变，收集绿色植被，人类活动，大气沉降等信息，建立t时间段总金属浓度Mt与三种信息量的数学模型，再将Mt合理加入到问题三中，使模型变得更为准确。

关键词：城市表层土壤；重金属污染；相关性分析

1.问题的重述

1

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

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

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

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

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

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

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

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

2.模型的基本假设

1.假设该城市各功能区内土壤性质无差异；

2.假设重金属污染的传播不受另一重金属污染物的影响，传播过程互不影响； 3.假设工业区中不考虑主干道路及其产生的污染；

4.污染物扩散过程中污染物浓度从圆心向四周递减变化； 5.采样点可以代表1平方公里的区域特征。

3.符号说明

2

Pji：采样点为j区i污染物的污染指数； Cji：j区i污染物的实测值； Si：i污染物的评价标准； P综：采样点的综合污染指数；

：采样点所有单项污染指数的平均值； P平均：采样点所有单项污染数的最大值； P最大f(x,y,z,?)：空间任意一点受某一重金属元素?污染的综合值的非线性拟合值； ai：空间一点到污染源的距离；

Fi'：第i点综合污染程度的实际值； Fi：第i点综合污染程度的理论值； ?i：变位参数；

S：理论值与实际值的差值的平方和； M：总金属浓度；

MAT：区域回归的大气沉降疏通量；

Ihum:人类活动引起的重金属增加的速率； Iplant：绿色植物吸收重金属的速率；

fitai：与f(xi,?)，f(xi,yi,?)或f(xi,yi,zi,?)（1?i?n）距离最近的值。

4.问题的分析

4.1 问题一的分析

先对附件二的数据筛选，剔除一些奇异值，以保证数据来源的合理性。接着运用MATLAB软件对研究区内实例数据进行插值，绘制8种重金属在城区的空间分布图。

采用国家土壤环境质量标准来评价各区域污染等级，先算出j区i污染物的污染指数Pji?Cji/Si污染起始值，接着以不同范围，求出所有单项污染指数的平均值，最后根据公式得到采样点的综合污染指数P综。由以上数据绘出不同重金属元素城区分布图和综合污染指数表，分析不同区域金属污染程度。 4.2 问题二的分析

推断重金属污染的主要原因复杂，考虑用因子分析和聚类分析两种分析方法来分析问题一中所得数据。由于变量的个数较多，因此之间可能存在关联，于是在因子分析时，做出8种金属重金属元素之间的相关性表Correlations。聚类分析是把那些污染程度接近的归为一类，运用MATLAB描绘出8种重金属元素地区污染分析图。最后由以上两种方法得到的图表和问题一中数据进行污染分析。 4.3 问题三的分析

重金属在土壤中的迁移方式可大致分为物理迁移、化学迁移和生活迁移，但考虑到模型的规律性，基于干湿沉降这种主要传播方式建立合理的重金属迁移模型。某金属元素可能会有多个污染源影响，每个污染源向外传播污染程度不同，因此，需要先对单个污染源研究分析。

污染源可能为中心污染区域内一点，污染传播以该点为中心，以半径递增污

3

染浓度递减方式向四周扩散。由此可得出空间中任意一点受到某元素单个污染源对它的污染，再以单个污染源为基础建立起所有可能的P个污染源对一点综合污染模型。用数值拟合方法得到某点的综合污染值，建立空间任意一点的综合污染的数学模型，非线性最小二乘法确定出参考的最优解。在此基础上，运用遍历搜索找出全城区各元素污染的污染严重值，这些严重值点即为若干污染源点。 4.4 问题四的分析

影响模型三的因素多种多样，如土壤含水率，绿色植物吸收金属速率等，收集绿色植被，农业活动，区域回归的大气沉降的信息，建立t时间段总金属浓 度Mt与三种信息量的数学模型，再将Mt合理加入到问题三中，使模型变得更为准确。

5.模型的建立与求解

5.1 问题一的求解

5.1.1 重金属元素的空间分布

根据题目中所给的附件一，运用MATLAB软件，进行二维插值，画出了该地 区的地貌图（程序见附录一）：

从图形中，我们可以初步判断出该地区地理特征，工业区集中在城市的低海拔地区，交通主干道密集， 生活区靠近工业区，山区人类活动相对较少，从而，我们大致可以估计出城市的重金属污染主要在工业区的工业污染，道路交通的废弃污染，人类生活污染。

进一步，由附件二，根据每个样本点重金属的观测数据，同样先利用Matlab软件进行插值拟合，分别描绘出重金属元素As,Cd,Cr,Cu,Hg,Ni,Pb和Zn在该城区的空间分布情况（以重金属As为例，全部图形见后面附录二）。

4

21.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104重金属As污染分布浓度线生活区工业区山区交通区绿地区2520151050-5-10Y

上面的平面图中，横纵轴的刻度是对点的定位，颜色的深浅表示污染程度的重与轻，每个图右方的图例精确地标出了重金属元素的数值（μg/g），方便了我们后面对不同区域重金属的污染程度的研究和对污染源的定位。

为了更好地反映出重金属污染物的空间分布情况，便于观察、分析、验证，借助与地理信息系统中的Arcgis软件，绘制出重金属空间分布图形。（以As为例，全部图形见后面附录三）

0.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934Xx 10

5

采用交叉验证法[1]（cross-validation）来验证插值的效果。即首先假定每一样本点的重金属浓度未知，都用周围样本点的值来估算，然后计算所有样本点的实际测量值与估算值的误差，以此来评判估值方法的优劣。采用平均误差（mean error,ME），平均绝对误差（mean absolute error,MAE），插值平均误差平方的平方根（Root mean squared interpolation error,RMSIE）作为评估插值方法的插值效果的标准。ME总体反映估计误差的大小，MAE可以估量估计值可能的误差范围，RMSIE可以反映利用样点数据的固执灵敏度和极值效应。RMSIE表达式：

RMSIE??(Zi?1na,i?Ze,i)2 n式中：Za,i为第i个样本点的实际测量值，Ze,i为估计值。N为用于检验的样本点数目。由此方法对上述插值方法进行检验，误差在允许范围内（如1%）。 5.1.2 不同区域重金属的污染程度

第一步，数据的筛选。剔除特异值时，先将区间平均值加减3倍标准值，从而得到一个范围，范围外则为特异值。

第二步，标准的制定。根据国家环境保护局科技标准司出台的土壤环境质量标准，得到如下标准值：

表1-1 土壤环境质量标准值[2] 级别 一级 二级 三级 土壤ph As 水田 As 旱地 Cd Cr 水田 Cr 旱地 Cu 农田 Cu 果园 Hg Ni Pb Zn 评价标准——土壤重金属评价分别采用研究区土壤重金属含量的自然背景值和国家一级标准值（为保护区域自然生态，维持自然背景的土壤环境质量的限制值）为参考标准。当评价标准为背景值时，相应的评价结果称之为重金属元素富集。当以国家一级标准为评价标准时，称之为重金属元素污染。

结论：经交叉验证法，检验若干数据后，与原数据对比可发现误差均小于1%。因此，我们认为，二维插值法得到的重金属元素污染分布图是有效的。

6

自然背景 15 15 0.20 90 90 35 -- 0.15 40 35 100 ?6.5 30 40 0.30 250 150 50 150 0.30 40 250 200 6.5~7.5 25 30 0.60 300 200 100 200 0.50 50 300 250 ?7.5 20 25 1.0 350 250 100 200 1.0 60 350 300 ?6.5 30 40 -- 400 300 400 400 1.5 200 500 500

表1-2 标准 污染级 As(μg/g) Cd(ng/g) Cr(μg/g) Cu(μg/g) Hg(ng/g) Ni(μg/g) Pb(μg/g) Zn(μg/g) 一级 15 200 90 35 150 40 35 100 二级 30 600 200 200 500 50 300 250 三级 40 1000 300 400 1500 200 500 500 第三步，方法的确定。评价方法[3]：将国家土壤环境质量标准中的一、二、三级标准分别作为土壤重金属污染评价的污染积累起始值、中度污染起始值和重度污染起始值，把土壤重金属污染分为4级（见表1-3）。

其中单项污染指数计算：

①若Cji?Si污染起始值，则Pji?Cji/Si污染起始值； ②若Si污染起始值?Cji?Si中度污染起始值，则

Pji?1?(Cji?Si中度污染起始值)/(Si中度污染起始值?Si污染起始值)；

③若Si中度污染起始值?Cji?Si重度污染起始值，则

Pji?2?(Cji?Si中度污染起始值)/(Si重度污染起始值?Si中度污染起始值)；

④若Cji?Si重度污染起始值，则

Pji?3?(Cji?Si重度污染起始值)/(Si重度污染起始值?Si中度污染起始值)。

22?P?P）/2?于是，综合污染指数计算：P综?（平均最大??。

1/2式中，Pji为采样点为j区i污染物的污染指数，Cji为j区i污染物的实测值，

Si为i污染物的评价标准，P综为采样点的综合污染指数，P平均为采样点所有单项污染指数的平均值，P为采样点所有单项污染数的最大值。 最大

表1-3 土壤重金属污染分级标准 污染等级 污染水平 Pji P综 一级 二级 三级 四级 Pji?1 1?Pji?2 2?Pji?3 P综?0.85 清洁 土壤轻度污染 0.85?P综?1.71 土壤中度污染 1.71?P?2.56 综Pji?3 7

P综?2.56 土壤污染相当严重

表1-4为不同的区域不同重金属元素的污染指数，表1-5为城区综合污染指数：

表1-4 各重金属元素污染指数 城区 生活区 As(μg/g) Cd(ng/g) Cr(μg/g) Cu(μg/g) Hg(ng/g) Ni(μg/g) Pb(μg/g) Zn(μg/g) 0.4107 工业区 0.4427 山区 0.2720 交通区 0.3587 公园绿地区 0.4127 1.1527 1.3608 0.7575 1.3341 1.1466 1.0609 0.5798 1.0890 1.2589 0.4903 0.2729 1.0955 0.8537 0.8280 0.7639

表1-5 城区综合污染指数 工业区 山区 0.5580 0.4958 0.4099 0.5772 0.4738 0.4438 0.4610 0.3638 0.4118 0.3793 1.0748 1.1427 1.0048 1.0993 1.0785 1.2959 1.5444 0.7089 1.5849 1.1194 城区 生活区 交通区 公园绿地区 P综 1.0852 1.2912 0.8049 1.1938 0.9787 第四步，重金属污染原因分析： 1.在生活区表层土壤综合污染指数达到轻度污染水平，以Cd和Cu富集程度相对较高，Pb和As的含量稍高于全国重金属元素背景值；

2.工业区污染指数接近中度污染水平，工业制造及燃煤产生的Hg、Zn等重金属元素对土地造成很大污染，在城区的综合污染指数1.2912是五个功能区中污染最大的地方；

3.由于山区无人为活动，城区综合污染指数0.8049，属于清洁区，不会对该地区造成污染；

4.在主干道路区车辆行驶频繁，轮胎与地面摩擦产生含Zn的较高颗粒物，以至于Zn在交通区的污染指数达到1.58，Zn在重金属元素均值图的交通区所占比例超过50%，分析得出Zn属于轻度污染；

5.因公园绿地区处于市区中心地带，水体、空气、地面垃圾对此影响较大，加上人群活动频繁，造成该区轻度污染。 5.2 问题二的求解

5.2.1 重金属污染的聚类分析和相关分析

在分析重金属污染的主要原因中考虑到不同金属重金属元素之间存在一定的相关性， 在SPSS软件下，得到相关性表Correlations。

8

表2-1 金属重金属元素相关系数

CorrelationsAs (μg/g)Pearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNPearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNPearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNPearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNPearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNPearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNPearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNPearson CorrelationSig. (2-tailed)Sum of Squares andCross-productsCovarianceNAs (μg/g)12908.5379.146319.255**.00055112.085Cd (ng/g)Cr (μg/g)Cu (μg/g)Hg (ng/g).255**.189**.160**.064.000.001.004.25155112.085173.30831912E+00712721.01340.003319.352**.00017649815550.254319115582804900.251319.532**.00019277476062.097319.103.066374342725025.48578.696319.397**.0004623870100922.4317.366319.265**.0003E+007Ni (μg/g)Pb (μg/g)Zn (μg/g).317**.290**.247**.000.000.0003026.6059.518319.329**.000234256.9736.657319.716**.000158404.3498.127319.495**.000254735.9801.056319.103.066530293.21667.589319131428.54198.832319.307**.00048549.701152.672319.436**.000467980.11471.63531913956.46943.888319.660**.000236494780550.481253.303319.431**.0001E+007Cd (ng/g)Cr (μg/g)173.30850619.421319319.189**.352**.001.00012721.01340.003319.160**.00425025.48578.696319.064.251100922.417649815550.254319.397**.000462387014540.471319.265**.0003E+00714540.47197039.041319319.532**.103.000.06619277476062.0973191844014226541.328319.417**.0004E+007110619.5319.495**.000254735.9801.056319.520**.0001348582374342711771.784319.417**.0004E+007110619.531918E+0082655400319.103.066530293.21667.589319.298**.00077331597436.94032908.914319319.383**.424**.000.000426597.43204043Cu (μg/g)1341.50110075.607319319.520**.387**.000.00013485826807305Hg (ng/g)4240.82421406.619319319.298**.196**.000.000773315924318.109319.307**.00048549.701152.6723191796837.82505.779319.494**.00026657488382.8553193E+007108254.4319.436**.000467980.11471.635319.494**.00026657488382.85531914E+007115078.7319Ni (μg/g)317.36697039.04111771.784319319319.317**.329**.716**.000.000.0003026.6059.518319.290**.00013956.46943.888319.247**.00080550.481253.303319234256.9736.657319.660**.00023649477436.940319.431**.0001E+00732908.914319158404.3498.127319.383**.000426597.41341.501319.424**.000320404310075.607319Pb (μg/g)Zn (μg/g)4240.82424318.109319319.387**.196**.000.000680730521406.6193193E+007108254.4319 由重金属元素的相关关系，把金属重金属元素分成四组，由问题一可知，山区污染几乎为0，因此不把它考虑在内。由表Correlations，按重金属元素相关程度密切原则，将金属重金属元素分成四组：1、Pb和Cd,Cu,Zn；2、Cr和Cu,Ni；3、 Ni 和Cr,Cu,Zn；4 、As。

聚类分析使用一定的相似距离测度方法不断把性质最接近的两类合并成更大的类，基本原则是使同一类中个体差异最小，不同类别之间差异最大。

运用聚类分析，描绘重金属聚类连线图，图中的1~8分别代表重金属元素As,Cd,Cr,Cu,Hg,Ni,Pb和Zn，连线代表了它们之间的相关性。从图中表现出的重金属元素的相关关系与相关分析得出的结论是吻合的。

**. Correlation is significant at the 0.01 level (2-tailed). 9

0.70.650.60.550.50.450.40.350.336427851

图2-1 重金属元素聚类图

5.2.2 重金属元素的污染原因分析

由问题一中对各区污染严重的重金属元素及问题二中四组相关性重金属元素数据，分别在各区绘制污染分布图（见图2-2,2-3,2-4,2-5）。

城市表层土壤重金属含量的高低不仅决定于本地区的地质背景，而且也受到城市工业活动、交通运输等因素的强烈影响。由问题一知山区综合污染指数0.8049。Pb的分布主要沿交通干道两侧及工业区呈现出较严重的带状污染（见图1），距离交通干道和工业区距离不同，重金属存在差异，如Pb在主干道区的单项污染指数1.09，工业区Pb污染均值为72.81，达到全城最高水平。经数据分析Cu重金属元素和很多重金属元素有很好的相关性，Zn重金属元素被广泛用作氧化剂和清洁剂，轮胎与地面的磨损可以产生含Zn的颗粒物，生活和工业污水中排放的Cu离子，因此工业区、主干道区、生活区的铜、锌污染极为突出。在燃煤为主要能源的工矿企业周围，Cr的平均污染指数接近0.5；但在生活区Cr污染严重，不妨推断工业废水和废气在人为或自然因素下进入生活区，人们不合理的生活方式可能也是一大原因。燃煤产生的Cd重金属元素在工业区产生了较大影响。Hg主要用于工业区制造科学设备（如气压计），相对其他区含量极多，Hg在工业区造成的污染是其他地区的2-5倍。全球约66%的精炼Ni用于制造不锈钢产品，As与其化合物被运用在农药、除草剂、与许多种的合金中，两者颇有联系，在生活、工业、交通区产生的影响大致一样。

相关性较高的重金属元素在交通区、工业区、生活区、公园绿地区的污染分布图：

10

21.81.61.41.210.80.60.40.2000.511.522.5x 1034x 1042交通区Pb元素Cu元素Zn元素x 1041.81.61.41.210.80.60.40.2000.511.52工业区Cr元素Ni元素2.5x 1034

图2-2 交通区分布 图2-3 工业区分布

21.81.61.41.210.80.60.40.2000.511.522.5x 1034x 104生活区Cu元素Zn元素21.81.61.41.210.80.60.40.20x 104公园绿地区Cd元素Hg元素00.511.522.5x 104

图2-4 生活区分布 图2-5 公园绿地区分布

由上图可以验证上述分析得出的结论：Pb在交通区污染比较严重，Cr在工业区污染比较严重，Cu在生活区污染比较严重。 5.3 问题三的求解

重金属在土壤中的迁移方式多种多样，人为因素而导致土壤中的金壤迁移的情况不具有规律性，因此不作讨论。而自然力量主导的重金属迁移又以干湿沉降、泥沙传播为主要传播途径。重金属迁移特征用此方法研究分析，建立模型。

图1-1到图1-8可知分布图，在某些区域重金属污染极其严重，不妨猜想污染源可能在其中心区域内一点，并以特定传播方式向空间传播，而空间中任意一点受到的某种重金属元素污染极可能是多个污染源对该点的综合污染，据此考虑，尝试得到某重金属元素单个污染源对一点的污染，再以单个污染源为基础建立起所有可能的p个污染源对一点综合污染模型。

fitai值的选定：ai是与f(xi,?)，f(xi,yi,?)或f(xi,yi,zi,?)（1?i?n）距离最近的值。金属污染等值线ai与等值线上金属污染浓度的对应关系：

11

rellrant?a1????m1?rellrant?m2?a2????空间任意一点受某一重金属元素?污染的综合值用非线性拟......?rellrant?a????mn?n合值f(x,y,z,?)表示。

污染程度Mn半径An金属元素可能污染源半径a1浓度等高线m1

图3-1单个污染源扩散示意图

（1）所有可能污染源连线都平行于x轴或y轴时，建立如下模型

f(xi,?)?fit?a1,a2???an??xixdx?rellrant?m1,m2???mn?

f(x,?)??1f(x1,?)??2f(x2,?)??????pf(xp,?)

（2）p个所有可能污染源只在x、y平面坐标对城区某点产生污染

12

f(xi,yi,?)??xixdx?dyyyifit?a1,a2???an??rellrant?m1,m2???mn?

f(x,y,?)??1f(x1,y1,?)??2f(x2,y2,?)??????pf(xp,yp,?)

（3）p个所有可能污染源在x、y、z空间中对某点产生污染

?f(x,y,z,?)?iiixixdx?dy?dzyzyizifit?a1,a2???an??rellrant?m1,m2???mn?

f(x,y,z,?)??1f(x1,y1,z1,?)??2f(x2,y2,z2,?)??????pf(xp,yp,zp,?)

用非线性最小二乘法确定参数[5]

非线性最小二乘法基本思想：根据某重金属元素可能污染源F(xi,yi,zi,?)的

坐标值及污染值求得空间若干点的某重金属元素综合污染值，通过与附件的若干点处的测量值进行比较，在0??i?10这个范围内以步长为0.01用遍历搜索法求出?1,?2????n的局部最优解，使得理论值与实际值的差值的平方和S?(F1'?F1)2?(F2'?F2)2?????(Fn'?Fn)2，当S取得最小值，此时?i(1?i?n)为所求的最佳值。即如下最小二乘拟合模型

minS(?i,?)??(Fi'?Fi)2，

i?1n最后得到变位参数的局部最优解。 用Matlab软件求出各参数（附录6），拟合出多项式： As:F(x,y)?7.9473?0.0003x?0.0002y Cd:F(x,y)?438.0403?0.0196x?0.0238y Cr:F(x,y)?102.7611?0.0098x?0.008y Hg:F(x)?1.1811?0.0001x Cu:F(x,y)?172.9358?0.222x?0.015y Ni:F(x,y)?21.7307?0.001x?0.0011y Zn:F(x,y)?396.9287?0.032x?0.0344y Pb:F(x,y)?122.4738?0.0095x?0.0058y 污染源位置的确定：

由上述模型，可以计算出全城区各点的某金属重金属元素浓度（污染量），则以步长为10，运用遍历搜索找出全城区各重金属元素污染的污染严重值（2-4

13

个为标准），这些严重值即为若干污染源。

经模型得到污染源的位置：

表3-1：污染源位置 污染源 坐标（m） 浓度值 As重金属元素的污染源 （18134,10046,41） 30.13(μg/g) （12696,3024,27） 23.72(μg/g) （4742,7293,9） 21.87(μg/g) Cd重金属元素的污染源 （21439,11383,450） 1619.8(ng/g) （1647,2728,6） 1092.9(ng/g) （4777,4897,80） 1044.5(ng/g) Cr重金属元素的污染源 （3299,6018,4） 920.84(μg/g) （3526,4357,70） 258.15(μg/g) Cu重金属元素的污染源 （1647,2728,6） 308.61(μg/g) Hg重金属元素的污染源 （2708,2295,22） 16000(ng/g) （8629,12086,1） 1801(ng/g) （10685,5528,34） 900(ng/g) Ni重金属元素的污染源 （3299,6018,4） 142.5(μg/g) Pb重金属元素的污染源 （1647，2728,6） 434.8(μg/g) （4777,4897,8） 472.48(μg/g) Zn重金属元素的污染源 （13797,9621,18） 3760.82(μg/g) （4948,7293,6） 1626.02(μg/g) 根据某重金属元素污染源的位置和污染浓度数据，不难发现，假设的存在实际污染源与猜想值有一定偏差；但若干可能污染源和实际测量污染源位置几乎重合，污染浓度数据值也很靠近，这在一定程度上验证了猜想。 5.4 问题四的求解

题三模型的优点：1）形象直观地反映出各地金属污染的主要来源；

2）通过模型可以从污染区域中找出污染源。

问题三模型的缺点：1）只对金属重金属元素的主要污染物进行传播分析，省去了若干离散的污染程度低的污染物研究；

2）在金属污染物的传播方式中，我们将各地区干湿沉

降、泥沙转移的传播方式标准化、统一化，没有考虑到土壤水分配系数对流速度等对金属转移的影响；此外，其他自然力（如水流、风力）对金属重金属元素转移没有进行分析。

信息收集及改进模型： 1.信息收集

由于工业化进程，导致部分工厂搬到主干路（交通方便），造成该地区地质As,Hg等重金属元素上升，生活区绿色植被显著提高，植物的根被对金属Pb、Cr等的吸收，使金属含量减少。工业区排放的废气中含有不少金属重金属元素，在风力（大气）作用下污染全城区。

信息数据：绿色植物吸收重金属的速率；农业活动引起的重金属增加的速率；区域回归的大气沉降疏通量。

2.改进

某重金属元素所有污染源对空间一点污染传播所有时间t?t2?t1

14

F(x,y,z,?)??1f(x1,y1,z1,?)??2f(x2,y2,z2,?)??????pf(xp,yp,zp,?)

F(x,y,z,?)??1f(x1,y1,z1,?)??2f(x2,y2,z2,?)??????pf(xp,yp,zp,?)?Mt在时间t内，Mt??(Ihum?Iplant)dt?MAt

t1t2

式中，Mt：总金属浓度；MAt：区域回归的大气沉降疏通量；Ihum：人类活 动引起的重金属增加的速率；Iplant：绿色植物吸收重金属的速率。

6.模型的评价

6.1 模型的优点 （1）利用Matlab软件插值方法绘图画出各重金属元素在城区的空间分布图，形象，清晰明了；

（2）对插值检验利用交叉验证法，很好地验证了插值方法的效果；

（3）利用重金属污染分级标准很容易对城区的重金属污染程度做出判断； （4）聚类分析，相关性分析很好地分析出各重金属污染之间的关系，说明了污染的来源；

（5）通过拟合，最小二乘法等相关分析，能够准确的找出污染源的位置。 6.2 模型的缺点

（1）样本点的取样限制会造成污染分析的误差；

（2）由于自然环境的不同，污染源的扩散方式会不同，会影响各样本的污染程度；

（3）影响污染源迁移的因素有很多，难以准确的描绘出污染的传播特征。

7.模型的改进和推广

7.1 模型的改进

我们可以在本文的基础上增加一些自然因素，地质特征对模型的改进，准确的找出污染源，对污染源的传播特性准确的描绘。

7.2 模型的推广

通过对重金属元素污染传播的模式的研究，我们建立的模型同样可以应用在油田污染，酸雨污染等环境保护工程。

15

8.参考文献

[1] Holdaway M R. Spatial modeling and interpolation of monthly temperature using Kriging. Clim. Res.,1996,24:1835-1845； [2]国家环境保护局科技标准司；

[3]李保杰等，基于地统计的矿业城市土壤重金属污染研究，江苏农业科学，524-558页，2011年第39卷第3期； [4]王学松 秦勇，徐州城市表层土壤中重金属环境风险测度与源解析，地球化学，92-93页，2006年第35卷第1期；

[5]赵静等，数学建模与数学实验，高等教育出版社，313-316页，2007年。

16

9.附录

附录一 城市三维立体图形程序 clc clear close all

D=load('zuobiao.txt'); x=D(:,2); y=D(:,3); z=D(:,4); c=D(:,5);

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'};

str={'等高线','生活区','工业区','山区','交通区','绿地区'}; % 三维体图 figure

h=surf(xi,yi,zi); set(h,'cdata',ci); colormap hsv

title('三维立体图') xlabel('X') ylabel('Y') zlabel('Z') 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

legend(str,'location','best') colormap summer

附录二 基于matlab的城市重金属污染空间分布图形

17

21.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104重金属As污染分布浓度线生活区工业区山区交通区绿地区2520151050-5-10Y0.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934Xx 1021.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104

重金属Cd污染分布浓度线生活区工业区山区交通区绿地区1000500y0-500-1000-15000.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934xx 10

18

21.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104重金属Cr污染分布浓度线生活区工业区山区交通区绿地区6004002000-200-400-600-800y0.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934xx 1021.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104

重金属Cu污染分布浓度线生活区工业区山区交通区绿地区2000150010005000-500-1000y0.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934xx 10

19

21.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104重金属Hg污染分布浓度线生活区工业区山区交通区绿地区14000120001000080006000400020000-2000-4000y0.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934xx 10-6000

21.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104重金属Ni污染分布浓度线生活区工业区山区交通区绿地区10050y0-50-1000.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934xx 10

20

21.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104重金属Pb污染分布浓度线生活区工业区山区交通区绿地区450400350300250200150100500-50y0.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934xx 1021.91.81.71.61.51.41.31.21.110.90.80.70.60.50.40.30.20.1x 104

重金属Zn污染分布浓度线生活区工业区山区交通区绿地区300025002000150010005000-500-1000-1500y0.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.934xx 10

21

附录三 基于地理信息系统Arcgis软件的重金属空间分布图形

22

23

24

25

附录四 表格

重金属元素 As (μg/g) Cd (ng/g) Cr (μg/g) Cu (μg/g) Hg (ng/g) Ni (μg/g) Pb (μg/g) Zn (μg/g) 表一 城市土壤重金属元素限值表 平均值 标准偏差 范围 3.6 0.9 1.8~5.4 130 30 70~190 31 9 13~49 13.2 3.6 6~20.4 35 8 19~51 12.3 3.8 4.7~19.9 31 6 19~43 69 14 41~97 标准值 2 8 38 290 2 153 80 2292

表二 平均值 城区 生活区 工业区 山区 交通区 As(μg/g) Cd(ng/g) Cr(μg/g) 6.16 6.64 4.08 5.38 261.06 50.22 344.33 44.62 151.50 36.89 333.65 51.95 258.65 42.64 269.84 45.27 5956.61 36.68 Cu(μg/g) Hg(ng/g) Ni(μg/g) 45.05 49.68 17.16 50.76 28.98 38.33 216.08 86.97 240.60 40.94 128.05 114.59 122.23 5488.61 17.75 18.44 14.55 16.47 15.17 16.48 2.72 Pb(μg/g) 54.82 72.81 36.26 61.32 55.81 56.20 175.40 Zn(μg/g) 144.38 181.66 70.89 187.73 117.91 140.51 2323.97 公园绿地区 6.19 平均值 方差 5.69 1.01 相关系数 As(μg/g) As(μg/g) 1.000 Cd(ng/g) 0.255 0.189 Cr(μg/g) Cu(μg/g) 0.160 Hg(ng/g) 0.064 Ni(μg/g) 0.317 0.290 Pb(μg/g) Zn(μg/g) 0.247 附录五 程序1

Cd(ng/g) 0.255 1.000 0.352 0.397 0.265 0.329 0.660 0.431 表三 相关系数 Cr(μg/g) Cu(μg/g) Hg(ng/g) Ni(μg/g) Pb(μg/g) Zn(μg/g) 0.189 0.160 0.064 0.317 0.290 0.247 0.352 0.397 0.265 0.329 0.660 0.431 1.000 0.532 0.103 0.716 0.383 0.424 0.532 1.000 0.417 0.495 0.520 0.387 0.103 0.417 1.000 0.103 0.298 0.196 0.716 0.495 0.103 1.000 0.307 0.436 0.383 0.520 0.298 0.307 1.000 0.494 0.424 0.387 0.196 0.436 0.494 1.000 重金属元素空间分布：

figure(1)

x=xlsread('D:\\sj.xls','附件1','B4:B318'); y=xlsread('D:\\sj.xls','附件1','C4:C318');

26

z=xlsread('D:\\sj.xls','附件1','F4:F318'); cx=[100:100:30000]; cy=[100:100:20000];

cz=griddata(x,y,z,cx,cy','v4'); contourf(cx,cy,cz); hold on;

x=xlsread('D:\\sj.xls','附件1','B4：B47'); y=xlsread('D:\\sj.xls','附件1','C4：C47'); scatter(x,y) hold on;

x=xlsread('D:\\sj.xls','附件1','B48：B82'); y=xlsread('D:\\sj.xls','附件1','C48：C82'); scatter(x,y) hold on;

x=xlsread('D:\\sj.xls','附件1','B83：B148'); y=xlsread('D:\\sj.xls','附件1','C83：C148'); scatter(x,y) hold on;

x=xlsread('D:\\sj.xls','附件1','B149：B283'); y=xlsread('D:\\sj.xls','附件1','C149：C283'); scatter(x,y) hold on;

x=xlsread('D:\\sj.xls','附件1','B283：B319'); y=xlsread('D:\\sj.xls','附件1','C283：C319'); scatter(x,y)；

set(gca,'XTick',[0:1000:30000]); set(gca,'YTick',[0:1000:20000]); title('重金属As污染分布')

附录六 程序2

分区重金属元素平均浓度： a=load('1qu.txt');

a1=a(:,5);a2=a(:,6);a3=a(:,7);a4=a(:,8);a5=a(:,9);a6=a(:,10);a7=a(:,11);a8=a(:,12);

ave1=mean(a1);ave2=mean(a2);ave3=mean(a3);ave4=mean(a4);ave5=mean(a5);ave6=mean(a6);ave7=mean(a7);ave8=mean(a8);

ave(1,:)=[ave1 ave2 ave3 ave4 ave5 ave6 ave7 ave8]; a=load('2qu.txt');

a1=a(:,5);a2=a(:,6);a3=a(:,7);a4=a(:,8);a5=a(:,9);a6=a(:,10);a7=a(:,11);a8=a(:,12);

ave1=mean(a1);ave2=mean(a2);ave3=mean(a3);ave4=mean(a4);ave5=mean(a5);ave6=mean(a6);ave7=mean(a7);ave8=mean(a8);

ave(2,:)=[ave1 ave2 ave3 ave4 ave5 ave6 ave7 ave8]; a=load('3qu.txt');

a1=a(:,5);a2=a(:,6);a3=a(:,7);a4=a(:,8);a5=a(:,9);a6=a(:,10);a7=a(:,11);a8=a(:,12);

ave1=mean(a1);ave2=mean(a2);ave3=mean(a3);ave4=mean(a4);ave5=mean(a5);ave6=mean(a6);a

27

ve7=mean(a7);ave8=mean(a8);

ave(3,:)=[ave1 ave2 ave3 ave4 ave5 ave6 ave7 ave8]; a=load('4qu.txt');

a1=a(:,5);a2=a(:,6);a3=a(:,7);a4=a(:,8);a5=a(:,9);a6=a(:,10);a7=a(:,11);a8=a(:,12);

ave1=mean(a1);ave2=mean(a2);ave3=mean(a3);ave4=mean(a4);ave5=mean(a5);ave6=mean(a6);ave7=mean(a7);ave8=mean(a8);

ave(4,:)=[ave1 ave2 ave3 ave4 ave5 ave6 ave7 ave8]; a=load('5qu.txt');

a1=a(:,5);a2=a(:,6);a3=a(:,7);a4=a(:,8);a5=a(:,9);a6=a(:,10);a7=a(:,11);a8=a(:,12);

ave1=mean(a1);ave2=mean(a2);ave3=mean(a3);ave4=mean(a4);ave5=mean(a5);ave6=mean(a6);ave7=mean(a7);ave8=mean(a8);

ave(5,:)=[ave1 ave2 ave3 ave4 ave5 ave6 ave7 ave8]; ave

附录七 程序3

污染指数： clear

s=[15 200 90 35 150 40 35 100 30 600 200 200 500 50 300 250 40 1000 300 400 1500 200 500 500];

c=[6.16 261.06 50.22 45.05 86.97 17.75 54.82 144.38 6.64 344.33 44.62 49.68 240.60 18.44 72.81 181.66 4.08 151.50 36.89 17.16 40.94 14.55 36.26 70.89 5.38 333.65 51.95 50.76 128.05 16.47 61.32 187.73 6.19 258.65 42.64 28.98 114.59 15.17 55.81 117.91 ];

for i=1:5 for j=1:8

if c(i,j)>s(3,j)

p(i,j)=3+(c(i,j)-s(3,j))/(s(3,j)-s(2,j)); end

if s(2,j)

p(i,j)=2+(c(i,j)-s(2,j))/(s(3,j)-s(2,j)); end

if s(1,j)

p(i,j)=1+(c(i,j)-s(1,j))/(s(2,j)-s(1,j)); end

if c(i,j)<=s(1,j)

p(i,j)=c(i,j)/s(1,j); end end end p

p=[0.4107 1.1527 0.5580 1.0609 0.5798 0.4438 1.0748 1.2959

28

0.4427 1.3608 0.4958 1.0890 1.2589 0.4610 0.2720 0.7575 0.4099 0.4903 0.2729 0.3638 0.3587 1.3341 0.5772 1.0955 0.8537 0.4118 0.4127 1.1466 0.4738 0.8280 0.7639 0.3793 ];

for i=1:5

pz(i)=sqrt((mean(p(i,:))^2+max(p(i,:))^2)/2); end pz

1.1427 1.0048 1.0993 1.0785 1.5444 0.7089 1.5849 1.1194

附录八 程序4

聚类分析;

load sj2.txt %把原始数据保存在纯文本文件sj2.txt 中 r=corrcoef(sj2) %计算相关系数矩阵

d=tril(r); %取出相关系数矩阵的下三角重金属元素 for i=1:10 %对角线重金属元素化成零 d(i,i)=0; end d=d(:);

d=nonzeros(d); %取出非零重金属元素 d=d';d=1-d; z=linkage(d); dendrogram(z)

附录九 程序5

重金属污染来源分析： plot(x,y,'o') hold on

plot(x1,y1,'*') hold on

plot(x2,y2,'*') hold on

plot(x3,y3,'x')

legend('交通区','Pb重金属元素','Cu重金属元素','Zn重金属元素') plot(x,y,'o') hold on

plot(x1,y1,'+') hold on

plot(x2,y2,'*')

legend('工业区','Cr重金属元素','Ni重金属元素') plot(x,y,'o') hold on

plot(x1,y1,'+')

29

hold on

plot(x2,y2,'*')

legend('生活区','Cu重金属元素','Zn重金属元素') plot(x,y,'o') hold on

plot(x1,y1,'+') hold on

plot(x2,y2,'*')

legend('公园绿地区','Cd重金属元素','Hg重金属元素')

附录十 程序6

拟合：

x = data(:,1:2); X=[ones(319,1) x]; B=zeros(3,1); for i=4:11

Y = data(:,i);

[b,bint,r,rint,stats]=regress(Y,X); B=[B b]; bint; stats; end B

30

本文来源：https://www.bwwdw.com/article/5w1r.html