城区表层土壤重金属污染分析模型南开大学数学建模

2024年1月9日发(作者：初一半期数学试卷)

城区表层土壤重金属污染分析模型

摘要本文针对城市表层土壤重金属污染的分析与评价，在综合考虑了各种不定因素的影响条件下，综合运用模糊数学对各种重金属的污染程度给出判定。利用基础几何与回归分析结合的方法建立了确定污染源的数学模型。并用数学软件Matlab对数据进行很好的处理，得到较为准确的数字信息，解决了所需解决的问题。

模型一中，对于重金属元素的空间分布及该城区内不同区重金属的污染程度的问题，文中根据单个元素在不同地域的污染程度，准确画出各种重金属元素的空间分布图。

依据土壤重金属污染程度模糊性，采用模糊隶属函数模型推导得出数据，辅以加权平均模型，得出各个区的污染级别。生活区、工业区、山区、主干道路区及公园绿地区的污染级别分别时三级、三级、二级、三级、三级，从而断言其污染程度。并通过Pearson相关性系数分析法来判断元素之间的相关关系，并结合重金属元素的空间分布图说明重金属污染的主要原因。

模型二中，对于重金属元素的传播特征的理解是它的分布特点和动态传播规律。元素分布的特点是通过模型一中的重金属元素的空间分布结合它们之间的相关性来描述的。元素的动态变化趋势函数是通过基于基本几何和区域散点线性关系得出的，即：

F(xi,i)＋x12tans{21＋x1[(x1x2)tan＋x2tan]2

1 ＋......＋xi[(xixj)tan＋xjtan]}2通过取As重金属污染浓度空间分布图中的一块图形进行定性的分析得出潜意识矢量方程y2212.81.0457x。通过计算得到污染源的浓度为27.1039ug/g。再通过将计算得出的浓度带入动态变化趋势函数得到坐标(4587.9,7010.4)。这个点的坐标就是污染源的坐标。

模型三，模型三基于物力、物理化学、流体力学等基础学科的知识建立了非保守性污染物质在土壤中迁移转化一维控制方程。通过公式的推倒，得到方程：

1xv\'txv\'tC(x,t)erfc()1N()

22D\'t2D\'t

得到的结论是:污染物浓度在时空上变化规律的近似解析解在形式上与正态概率分布函数相同的。

关键字模糊数学 Pearson相关性系数分析法隶属函数一维弥漫方程

渗透方程

1．问题重述

1．1 问题的提出

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

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

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

(2) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置，并分析所建立模型的优缺点。

(3)为了更好地研究城市地质环境的演变模式，还应收集什么信息？以及有了这些信息，如何建立模型解决问题？

1．2 问题的分析

针对问题一不同的功能区对土壤的影响大小是不同的。要评定8种重金属在该城区的空间分布和对该城区不同区域的污染程度，仅用数据描述过于抽象。因此采用matlab 成图技术直接观察这8种主要重金属元素在该城区的空间分布，这样比仅用数据描述更为形象。

在评定重金属元素污染程度的时候，由于考虑到元素毒性问题，不能仅凭元素所占含量来划分它们的污染程度，因此必须对元素进行权重的衡量。

要说明重金属元素污染的主要原因，需要结合它们的空间分布图来说明。不同元素的污染源不同，但是它们之间存在着一定的相关关系，从这种关系里找出它们污染的主要原因。

针对问题二要得到传播方程，必须从污染物浓度随传播路径的变化规律中找到一般的方程。它们的关系应类似于Vi(Xi,Yj,Zk,,Ld)的关系。得到传播方程后，对给定的区域取一些离散的点，并做出回归方程。由方程进一步推倒出污染源的位置。

针对问题三由于问题二中所建立的方程过于简陋，在现实生活中运用不大，只能作为一个理想的模型。问题三所要解决的就是，实现方程与实际生活关系。能运用方程解决现与实相关的问题。

2．模型的假设与符号的约定

2．1 模型的假设

1 每个监测点一种重金属污染浓度平均值就是该区污染的评价值

2 背景值作为该区的污染评价标准

3 元素与元素之间不发生化学反应，它们之间的性质在短时期内不会改变

2．2 符号的约定

Ci 第i种污染物的检测浓度

Si 第i种污染物评价标准浓度

xi 该重金属含量的实测值

wi 第i个因子的权重(i1,2,3,4,5)

(xi) 某种重金属元素对一级重金属环境质量的隶属度函数

(xi) 某种重金属元素对二级重金属环境质量的隶属度函数

(xi) 某种重金属元素对三级重金属环境质量的隶属度函数

ai 该重金属对应于一级土壤重金属环境质量状况的标准值

bi 该重金属对应于二级土壤重金属环境质量状况的标准值

ci 该重金属对应于三级土壤重金属环境质量状况的标准值

U 污染物评价因素集合

V 评价等级集合

R{UV}组成一个模糊矩阵,即关系模糊矩阵

S1,S2,S3该指标对应的各土壤重金属环境质量级别的标准值

3．模型的建立与求解

3．1 建模过程

模型1.1

对于重金属元素在该城区的空间分布，本文采用matlab 绘制出实验的319个样点的GPS坐标，生成相应的散点图，录入各点重金属浓度，并以此浓度作图生成等值线图，最终得到重金属空间分布格局图。具体见附录。

8种主要重金属元素在该城区的空间分布图

图一土壤重金属As含量（μg/g）空间分布图

图二土壤重金属Cd含量（ng/g）空间分布图

图三土壤重金属Cr含量（μg/g）空间分布图

图四土壤重金属Cu含量（μg/g）空间分布图

图五土壤重金属Hg含量（ng/g）空间分布图

图六土壤重金属Ni含量（μg/g）空间分布图

图七土壤重金属Pb含量（μg/g）空间分布图

图八土壤重金属Zn含量（μg/g）空间分布图

注：点（.）、小圆（o）、叉子符（x）、加号（+）、星号（*），分别为生活区、工业区、山区、交通区、公园绿地区区的标志。

模型1-2

对题目所给的附录2表格中，筛选出各金属元素的区域分布情况，并把单个区域所存在的元素浓度取平均值得到各检测点浓度表。同时，把附录3中的背景值浓度做5个划分，分别为无污染、尚污染、轻污染、中污染和重污染。

表一城区各功能区监测点评价浓度表

As Cd Cr Cu Hg Ni Pb Zn

(μg/g)

(ug/g)

(μg/g) (μg/g)

(ug/g)

(μg/g) (μg/g) (μg/g)

一区 6.2705 0.28996 69.018 49.403 0.093041 18.342 69.106 237.01

二区 7.2514 0.39311 53.409 127.54 0.64236 19.812 93.041 277.93

三区 4.0441 0.15232 38.96 17.317 0.040956 15.454 36.556 73.294

四区 5.708 0.36001 58.054 62.215 0.44682 17.617 63.534 242.85

五区 6.2637 0.28054 43.636 30.192 0.11499 15.29 60.709 154.24

表二该地区土地污染等级评价标准

一级（背景二级

值）

无污染尚污染轻污染中污染

1.8 2.7 3.6 4.5

0．07 0．1 0．13 0．16

13 22 31 40

6 9.6 13.2 16.8

0．019 0．027 0．035 0．043

4.7 8.5 12.3 16.1

19 25 31 37

41 55 69 83

评价等级

污染程度

As (μg/g)

Cd (ug/g)

Cr (μg/g)

Cu (μg/g)

Hg (ug/g)

Ni (μg/g)

Pb (μg/g)

Zn (μg/g)

三级

重污染

5.4

0．19

20.4

0．051

19.9

由于考虑到元素的毒性，必须对各个元素进行权重的确定，确定权重的公式这里采用反映土壤各种重金属元素相对含量大小的加权法该计算权重的方法在一定程度上可以反映污染超标的轻重对因子权重的影响该方法计算权重的一般公式为：

WiCi/SiCi1n

i/Si并用matlab编写权重矩阵。

表三: 由matlab计算得到元素毒性权重表

(μg/g)

0.08817

(ug/g)

0.11291

(μg/g)

(ug/g)

(μg/g)

0.1127 0.18945 0.13456 0.07548 0.11284 0.17388

0.0371 0.06912 0.09277 0.04639 0.06964 0.03968 0.22255 0.42272

0.11785 0.12292 0.13185 0.13763 0.12277 0.13181 0.12372 0.11144

0.05163 0.09018 0.06098 0.15348 0.41573 0.04664 0.06674 0.11461

0.10664 0.13227 0.08627 0.14019 0.20137 0.07619 0.12003 0.13701

为了进行模糊运算需要确定隶属度函数并以隶属度来描述土壤污染状况的模糊界线使用表2的值作为各个隶属度函数的拐点,用分段的直线函数模拟,可以得到不同PH 值下各种重金属污染状况的隶属度函数. 一般来说 ,各种土壤重金属环境质量状况的隶属度函数可以用下面三个分段函数表示某种重金属元素对一级重金属环境质量的隶属度函数：

1xiai(xi)(bixi)（/biai）aixibi

0xibi某种重金属元素对二级重金属环境质量的隶属度函数：

0(xi)(xiai)/(biai)(cx)/(cb)iiiixiai,xiciaixibi

bixici某种重金属元素对三级重金属环境质量的隶属度函数：

0xibi(xi)(xibi)（/cibi）bixici

1xici加权平均模型

bjwirij,j1,，m

i1n式中bj为最终评价结果对应于第j个等级的隶属度，wi为对应的权重，rij为模糊关系矩阵R中的对应元素，n为参评因子个数，m为所划分的等级数，类似的计算得到评价向量B(b1,b2,...bn)由于该模型计箅结果已经自动归一化，所以以该集合中最大值所对应的级别作为最终评价结果[3]。

用Matlab计算模糊矩阵M：

对于矩阵R,R{UV}，有5个矩阵，分别有

01000000.75330000.628301000000.55780100001000，R，R00.3782R1230000.625010000.2050.79500.01160.988400.5850000.537010001000000.8466010100000101000.2980.70201001001，R

R4500010100.30040.699600.60660.39340001010010010.24670.37170.44220.5718

0.3750.4150.4630.1534模糊矩阵复合运算模型的评价：

评价向量的计算：计算各采样点的关系模糊矩阵和权重模糊矩阵，分别用单因素决定模型和加权平均模型进行评价。对评价结果进行归一化处理，得到各采样点的评价向量

评价分值的计算：根据土壤重金属环境质量各级别的控制意义，取土壤重金属环境质量一级水平的分值为100,二级水平的分值为80三级水平的分值为60,按照公式(4) ,取k 值为2,计算在两种评价模型下各采样点的评价分值。从评价分值可以直观的看出各采样点土壤重金属环境质量的优劣，可比性大大增强。

表四各区污染程度分析表

评价模型

加权平均模型

评价结果一区二区三区四区

一级隶属 0 0 0 0

二级隶属

0.0154 0.00043036 0.60623 0.014011

三级隶属

0.9845

评价等级三级

污染程度严重污染

评价分值 60.31

0.99957

三级

严重污染

60.01

五区

0.071929

0.38688 0.98599 0.92807

二级三级三级

中度污染重污染重污染

71.71 60.28 61.44

3．2．1模型解释

依据表四的结果，按照最大隶属的原则，可以评价出各个区的土壤重金属环

境质量级别。但需要注意的是，当某个区具有两个最大隶属度时，需要考虑次大值更靠近哪个最大值，以该最大值对应的级别为最终评价结果。

模型优缺点：加权评价模型主要适用于各个参评因素超标情况接近，即不存在单因素否决，且评价出发点为希望评价体现不同参评因子对土壤重金属环境质量的综合影响的情况。而单因素决定模型适用于个别参评重金属元素超标过大，严重影响土壤环境质量，评价出发点为希望评价体现单因素否决的情况。根据实际情况，采用加权评价模型能够更准确的评价环境污染程度，但是对于山区的评价有待商榷。因为环境没有过多的人类活动的时刻影响，所以外来的重金属沉降积累很容易占较大的权重。由此，评价结果缺乏可信度。

3．3重金属污染的主要原因的分析

使用Matlab中的Pearson相关性系数分析法，对8种重金属元素，按功能区分类，进行两两之间的相关性指数分析[2]。

表五土壤中8种元素的相关性

As Cd Cr Cu

1 0.255** 0.189** 0.160**

1 0.532**

Hg Ni Pb Zn

0.064 0.317** 0.290** 0.247**

0.103 0.716** 0.383** 0.424**

0.103 0.298** 0.196**

1 0.307** 0.436**

1 0.494**

1 0.352** 0.397** 0.265** 0.329** 0.660** 0.431**

1 0.417** 0.495** 0.520** 0.387**

Zn 1

**. 在 .01 水平（双侧）上显著相关。

观察发现，7种元素两两之间的相关性十分明显。除了之间的相关性为显著相关外，其余的重金属之间都达到了极显著相关关系。

将重金属元素两两为一组，按照相关系数R值的由高到低来排列，得出前7位的顺序：Ni-Cr>Pb-Cd>Cu-Cr>Pb-Cu>Zn-Pb，它们之间的相关系数分别为：0.716，0.660，0.520，0.494，0.436。其余重金属间的相关系数都小于04。

由此表明：Ni-Cr Pb-Cd Cu-Cr Pb-Cu Zn-Pb。这7组重金属之间的相关性实验中，对土壤重金属分布起到主导作用。

相关性分析可以帮助揭示重金属元素的污染来源，按土壤所在的不同功能区分类，将该城区中的319个样点的重金属数据进行两两之间的相关性指数分析。

表六生活区中8种元素的相关性

As Cd Cr Cu Hg Ni Pb Zn

1 0.381* 0.238 0.531** 0.293 0.605** 0.450** -0.017

1 0.349* 0.499** 0.397** 0.283 0.802** 0.346*

1 0.376* 0.150 0.527** 0.416** 0.412**

1 0.198 0.434** 0.502** 0.238

1 0.211 0.340* 0.242

1 0.300* 0.334*

1 0.328*

*. 在 0.05 水平（双侧）上显著相关。

**. 在 .01 水平（双侧）上显著相关。

观测发现，生活区土壤金属元素中与As含量与Ni、Pb含量，Cd含量与Cu、Pb含量，Cr含量与Ni、Pb、Zn含量，Cu含量与Ni、Pb含量, 达到了O．01水平的极显著正相关。相关系数分别是0.605，0.531,0.450，0.499，0.802，0.527，0.416，0.412，0.434，0.502。As含量与Cd含量，Cd含量与Cr、Zn含量，Cr含量与含量含量，Hg含量与Pb含量, Ni含量与Pb、Zn含量, Pb含量与Zn含量, 达到了O.05水平的显著正相关。

由于Cd和Pb主要来自交通污染源，而Cr、Zn来自自非交通污染源。这说明了在生活区，以工业污染源为主的重金属元素之间有着明显的交互作用，交通污染源为主的重金属元素之间也有着明显的交互作用，但这两种污染源之间的交互作用并不明显。

表七工业区中8种元素的相关性

As Cd Cr Cu Hg Ni Pb Zn

1 0.329 0.380* 0.153 0.181 0.690** 0.395* 0.518**

1 0.541** 0.566** 0.533** 0.489** 0.829** 0.754**

1 0.920** 0.902** 0.698** 0.675** 0.695**

1 .983** 0.503** 0.670** 0.622**

1 0.479** 0.612** 0.590**

1 0.578** 0.634**

1 0.739**

*. 在 0.05 水平（双侧）上显著相关。

**. 在 .01 水平（双侧）上显著相关。

观测发现，工业区土壤金属元素中，As含量与与Ni含量，Cd含量与Pb、Zn含量，Cr含量与Cu、Hg、Ni、Pb、Zn含量，Cu含量与Hg、Pb、Zn含量，Hg含量与Pb含量，Ni含量与Zn含量，Pb含量与Zn含量，达到了O．01水平的极显

著正相关。相关系数分别是：0.690，0.829，0.754，0.920，0.902，0.698，0.675，0.695，0.983，0.670，0.622，0.612，0.634，0.739。As含量与Cr、Pb含量，达到0.05水平的显著正相关。相关系数是：0.380，0.395。

这说明了，在工业区，无论是交通污染源内部和工业污染源内部，还是交通污染源与工业污染源之间，都存在着明显的交互作用；体现出在这一区域范围内，2种污染源对土壤重金属的共同作用。

表八山区中8种元素的相关性

As Cd Cr Cu Hg Ni Pb Zn

1 -0.291* 0.113 0.527** 0.075 0.078 -0.205 -0.176

1 0.066 0.090 0.246* 0.049 0.766** 0.606**

1 0.364** 0-.006 0.945** 0.107 0.627**

1 0.505** 0.358** 0.122 0.252*

1 -.045 0.226 0.170

1 0.028 0.629**

1 0.590**

*. 在 0.05 水平（双侧）上显著相关。

**. 在 .01 水平（双侧）上显著相关。

观测发现，山区土壤重金属元素中，As含量与Cu含量，Cd含量与Pb、Zn含量，Cr含量与Ni、Zn含量，Cu含量与Hg含量，Ni含量与Zn含量，Pb含量与Zn含量，达到了O．01水平的极显著正相关。相关系数是：0.527，0.766，0.606，0.945，0.627，0.505，0.629，0.590。Cd元素与Hg元素，Cu元素与Zn元素，达到了0.05水平的显著正相关。相关系数是：0.246，0.252。

由于Cr与Ni、Zn主要来自于汽车运输，而其他元素没有极显著相关或着存在负相关关系。且有六种元素的实际值在背景值的区间里。故而可显示出，山区土壤污染源主要是交通运输，而受工业污染源的影响较小。

表九交通区中8种元素的相关性

As Cd Cr Cu Hg Ni Pb Zn

1 0.121 0.139 0.092 -0.004 0.228** 0.060 0.188*

1 0.373** 0.424** 0.211* 0.351** 0.615** 0.294**

1 0.894** 0.012 0.869** 0.428** 0.395**

1 .032 0.886** 0.506** 0.432**

1 0.040 0.266** 0.118

1 0.396** 0.503**

1 0.482**

**. 在 .01 水平（双侧）上显著相关。

*. 在 0.05 水平（双侧）上显著相关。

观测发现，交通区土壤重金属元素中，Cd含量与Cu、Pb含量，Cr含量与Cu、

Ni、Pb含量，Cu含量与Ni、Pb、Zn含量，Ni含量与Zn含量，Pb含量与Zn含量，达到了O．01水平的极显著正相关。相关系数分别是：0.424，0.615，0.894，0.869，0.428，0.886，0.506，0.432，0.503，0.482。As含量与Zn含量，Cd含量与Hg含量，达到了0.05水平的显著正相关。相关系数是：0.188，0.211。

由于Cd、Pb、Cr、Ni、Zn元素主要来自于交通污染源，而且Cd与Pb显著相关较大，据此推测可能来自于工业污染源。Hg与其他的元素出现负相关关系，表明Hg存在着不同的污染源。

表十公园绿化区中8种元素的相关性

As Cd Cr Cu Hg Ni Pb Zn

1 0.358* 0.689** 0.107 0.176 0.691** 0.265 0.285

1 0.564** 0.500** 0.054 0.433** 0.598** 0.712**

1 0.357* 0.023 0.739** 0.397* 0.509**

1 0.136 0.267 0.756** 0.521**

1 -0.048 0.389* 0.063

1 0.168 0.298

1 0.748**

*. 在 0.05 水平（双侧）上显著相关。

**. 在 .01 水平（双侧）上显著相关。

观测发现，公园绿化区土壤重金属元素中，As含量与Cr、Ni含量，Cd含量与Zn含量，Cr含量与Ni含量，Cu含量与Pb含量，Pb含量与Zn含量，达到了O．01水平的极显著正相关。相关系数分别是：0.689，0.691，0.712，0.739，0.756，0.748。As含量与Cd含量，Cr含量与Cu、Pb含量，Hg含量与Pb含量，达到了O．05水平的显著正相关。相关系数分别是：0.358，0.357，0.397，0.389。

由于Cr、Ni、Cd、Zn、Pb来自于交通污染，且有以交通污染源为主的Cu、Pb、Zn之间存在着明显的交互作用，显示出居民区的土壤重金属受到交通污染源的控制，受工业污染源的影响较小。Ni与其他元素出现负相关关系，说明Hg还存在着其他污染源。Hg与其他元素没有相关性，其含量与公园绿化区的背景值最为接近，

说明Hg它很有可能来源于自然环境中的岩石风化。

3．3．1 土壤重金属污染评价

综合以上信息分析得出：这八种土壤重金属含量分布比较均匀。将全部重金属元素与背景值相比较知道均比背景值大出许多，此地存在中污染区和重污染区。

在生活区，重金属元素比重差异不是太大(由权重表可知)，但含量却是不小，均大于背景值，但Zn比重最大，占了17.388％。在工业区，重金属元素比重差异很大，尤是Hg元素的比重最大，占了42.272％。通过资料查询知晓，Hg主要来自于工业废水，故而可知工业区的废水排放很不符合环境标准。观察浓度表，所有元素含量都大大超过背景值。在山区，各元素权重最为均匀，可是依然处于污染状态。在交通区，Hg元素的权重亦是最大，是41.573％，来自于交通污染源的各元素交互作用明显，含量也是不可小觑，其值都超过了背景值。在公园绿地区来自于交通污染源的重金属元素和来自于非交通污染源的重金属元素，或者各

自内部之间都存在明显的交互作用，体现出不同污染源对土壤重会属分布的共同作用。由此对比该地区土地污染等级评价标准，可得生活区、工业区、交通区、公园绿地区都达到了重污染水平，而山区也是中污染水平。

综合八种重金属元素含量的大小，不难得出五个区的污染程度排序：工业区 > 生活区 > 交通区 > 公园绿地区 > 山区。反映了重金属污染从地区内部到远郊递减的规律。

3．4 模型二

3．4．1 传播特征的分析

本文对题中传播特征的解释是在这片地域的重金属污染物的分布特点和动态变化趋势。

1.分布特点：

对于As这种元素，它的分布主要集中于工厂较多的工业园区，和道路交通较为发达的地区。

对于Cd这种元素，它的分布主要集中于交通要道地区附近，并且有扩散的潜能，如在公园河山区土壤中也出现了这种元素，在从元素的空间分布图来看，工业区主要建在山区附件。因此，山区也有也存在了这种元素。类似，公园绿地和居民区都存在较多的重金属元素Cd。

对于Cr这种元素，它的分布主要集中于交通地区，和工业园区。Cr主要用于合金(如供铬钢用的铬铁合金或如镍铬电阻合金)及电镀(如汽车减震器和内部装璜以及切削工具)，电镀和制造特种钢。而Cr的来源主要以铬铁矿或铬铅矿形式存在。因此，工业地区可能较多。

对于Cu这种元素，它的分布主要集中于工业区，由于Cu的用途之广，所以该城区的污染浓度非常大

对于Hg这种元素，它的分布主要集中于工业地区，以及对于汞矿的开采与汞的冶炼。

对于Ni这种元素，它的分布主要集中在交通，工业园区，分布面积广泛。

对于Pb这种元素，它的分布同Cd，是重工业地区污染排放物之一

对于Zn这种元素，它的分多存在于交通地区。

2．一般动态变化趋势

同一污染源，不同元素的迁移变化趋势是不同的。它是一个类似于

Vi(Xi,Yj,Zk,,Ld)的函数，即三维、浓度、距离关于平面弥漫速度的关系。在不考虑特殊的情况下，重金属浓度越高，它的扩散速度越大，扩散的距离就越远。同时考虑到地形影响Vi的情况，在海拔高的地方，污染物迁移速度变慢，污染较轻（此处无污染源）。海拔低的地方，污染物迁移速度快。

3．4．2 模型假设

1．城区的重金属污染是一种连续的稳定场

2．污染源以点污染源为主

3．未知的污染源的扩散方式是同心圆模式

4．污染物的浓度变化方式不受高度影响，并且以相同的梯度减小

5 对于单个污染物浓度来说，它的最大点就是该地区的污染源

6 不考虑高度对污染物扩散方式的影响

3．4．3符号的约定

xn 横坐标

yn 纵坐标

 各重金属的浓度

 能量极大点

 行走极限点

Li 行走极限点到能量极大点的距离

 浓度变化

li 颇高

3．4．4 模型建立

对第一问中，取As重金属污染浓度空间分布图中的一块图形进行定性和定量的分析。

由于污染源的不确定性，假设该地区的污染源为(x0,y0)。在这一扩散过程中，选取浓度为背景值的点(x1,x2),再选取周围大于背景值的点(x2,y2),(x3,y3)……(xn1,yn1),(x0,y0),其中,(x0,y0),(xn1,yn1)为我们选取点中浓度比较大的一个或某几个点。分别作出它们的线型关系。

由于污染源是以等浓度梯度/s扩散的一种源，且不考虑高度的影响，因此当距离为Li时，这点的浓度就是si/s。

如下示意图所示：

上图为扩散的理想模型，背景值点是某一元素达到一级标准的点，我们把这个点称为行走极限点，记为。当这个元素到达这个点的时候，它就不在向外扩散了。而位置浓度扩散点就是我们要找的污染源

现假设污染源的扩散能力极强，足以穿透10cm的地面，我们把这些点称为能量极大点，记为。随着污染物的迁移，R逐渐的增大，浓度逐渐变小，这是一个及其漫长的过程，我们把这种过程成为弥漫时间。

现在，由已知的一系列点(x2,y2),(x3,y3)……对其做一维线性回归方程y01x,称这个方程为的潜意识一维矢量方程，的过错称为记忆过程。我们就知道了这个方程的任意的坐标。

：

在问题一中的元素空间分布图中取任意元素的局部等浓梯度变化图，找出该元素的背景值浓度，以该点为点向比它浓度更大的地方做的潜意识一维矢量方程y01x，其中x,y是有界域的值。

任意带入该界域的x值，称这对点为的原始最终点。当该点以最接近点的方式进行移动时。

污染物浓度扩散的分布是均等分布的，因此它的浓度和上图所切割的理想面1积成正比，点形成的浓度梯度角为tan（L以米计算）。

10Lx12tan当取x1时，l1x1tan,面积s1；

2当取x2时，l2x2tan,面积s2……

1x1[(x1x2)tan＋x2tan]

21xi[(xixj)tan＋xjtan]

2当取xi时，lixitan ，面积si面积和浓度的关系

/s为点的浓度

那么把L分为n个等分的x段，那么每段的浓度为/s0,

某一点的浓度为ss

ii1n到此我们确定了简单的一维弥漫方程：

F(x,i)ss

ii1nF(xi,i)x12tans{21 ......＋xi[(xixj)tan＋xjtan]}21＋x1[(x1x2)tan＋x2tan]2

做法如下图：

当确定了的潜意识一维矢量方程时，在界域内给定x点，就可由简单的一维弥漫方程得到它的浓度i，和相比：

如果i0时,我们则认为过程不够彻底，并继续朝矢量方程的方向前进。直到i时，我们则认为的记忆过程基本完成。若并确认这点为所要找的污染源，这种通过随机取点来找污染源的方式所找的污染源是不准确的点污染源，而已有一定半径的面污染源。

模型运用

现以As的部分图来说明：

作出上图的一维矢量方程y2212.81.0457x

任意给出多个点的坐标a(5200,7650.44),b(5000,7441.3),c(4800,7232.16),

d(4700,7127.59),(4742,7293),带入一维弥漫方程

F(xi,i)x12tans{21＋x1[(x1x2)tan＋x2tan]2

1 ......＋xi[(xixj)tan＋xjtan]}23tan25tan2(2n1)tan2xx......x)s20L222

x35(2n1)2(x2x2......x)s20L20L20L20L(x

As的背景值点为5.4ug/g，由浓度对比图看出该的横坐标为5200，并从5200开始，以直线y2212.81.0457x向左前进，通过ab到达c得到x578.76，通过a到达得到L580.9,最后计算出si28.8279平方米，s29.0484平方米，21.87ug/g,5.4ug/g,所以最后得到的浓度为F(xi,i)＝27.1039ug/g。本文认为这点就为基本的污染源点。最后通过这点反求出x的值为，带入上面的矢量方程得到y，所以污染源的坐标是

(4587.9,7010.4)。

考虑到污染源的模糊性，以c为半径左r84.1m的模糊污染源。并在这个理想的模糊污源里用以上的方法继续寻找最终的理想点污染源。

3．5 模型三

城市地质环境演变模式：城市地质环境同自然地质环境的相同点在于它们拥有相同的水圈、大气圈和岩层圈，由于重金属元素对大气圈和岩层圈的影响几乎可以忽略，我们就认为城市地质环境的演变中，重金属元素只影响水圈。

模型二中没有考虑到重金属元素的对于大地的渗透作用，实际上，重金属元素的渗透也是不可忽视的。如果他们过多被排放，就会透过渗透作用到地下水层，严重的影响了地下水层的水质。因此我们要控制排放量，通过减少排放量从而减少它对地下水层的污染。

模型三基于物力、物理化学、流体力学等基础学科的知识建立了非保守性污染物质在土壤中迁移转化一维控制方程。

还需要收集的数据还要：

Ss为比贮水量

污染组分i在地下水中的浓度

介质体积含水率，对于饱和土壤，则是介质有效孔隙度

Dij

弥散系数张量

方向的实际渗流速度或平均孔隙流速，m/s；

Rk,i n个不同反应中第k个反应的溶解污染组分i的产率

笛卡儿坐标

模型建立

在稳定均匀流及x轴与流速向量方向一致的情况下，如果不考虑给排水和生物降解作用，并忽略弥散系数与空间和时间的关系，可由

nCiCiCi\'W*

(Dij)(Civi)Rk,i

txixjxik1

得到非保守性污染物质（与土壤介质发生吸附、沉淀或生化反应等作用）在饱和土壤中迁移转化的一维控制方程：

C2CCnD2vRk

txk1x对于有机物污染问题，吸附作用是影响污染物迁移转化的主要因素。如果把吸附和解吸完全看作一个可逆过程，则：

Rkk1nbS

et

如线性平衡吸附方程CsKdC能够很好地描述介质的吸附特性，那么上式中的Rk可表示为：

k1nRkk1nbCsC

bKdtt

代入计算整理后可得：

bC2CC[1Kd]D2v

txx令Rd1(b/)Kd，D\'D/Rd，v\'v/Rd，其中Rd为阻滞系数，表示多孔介质对于非保守性物质的阻滞能力，则上式可简化为：

C2CCD\'2v\'

txx在一维、稳定流场和连续注入定污染物浓度（C0）时，方程的定解条件如下：

C(x,0)0 0＜x＜

C(0,t)C0 0＜t＜

C(,t)0 0＜t＜

借助于Laplace变换，并代入上述边界条件，方程的解析解为：

xv\'txv\'tC(x,t)1v\'xC(x,t)erfcexp

erfc\'\'C02D\'2Dt2Dt在xv\'/D\'>10的情况下，式（2-19）中右端第二项比第一项小的多，可以忽略不计，这时污染物浓度在时空上的变化规律可以用以下近似解表示：

C(x,t)C1xv\'terfc()

C022D\'t经过适当变量替换后，式（2-20）可以表示为：

C(x,t)1xv\'txv\'terfc()1N()

22D\'t2D\'t

式中erfc(x)和N(x)分别为余误差函数和正态分布函数，即：

erfc(xv\'t2D\'t)2xv\'texp(y2)dy

2D\'t2y2D\'texp(N())dy

22D\'t2xv\'t1xv\'t由上述求解可以看出，污染物浓度在时空上变化规律的近似解析解在形式上与正态概率分布函数相同，所以目前分析室内一维、稳定流场和连续注入定示踪剂浓度的实验数据，以确定多孔介质弥散系数的方法主要有正态分布函数法和反函数法。

但是正态分布函数法存在明显缺点，而反函数法在考虑示踪剂和多孔介质之间的吸附与解吸作用时不能够应用。为此，本研究借用反函数变换思路，在利用保守性物质确定土壤介质的有效孔隙度基础上，通过相关计算确定非保守性物质[6]BTEX在介质中的弥散系数D和阻滞系数Rd。

模型的评价与改进

模型一的优点：在元素浓度分布图中，能直观的反映单个元素在该区的污染程度；在重金属污染的评价中，考虑到元素毒性的特点，建立判定元素毒性的权重矩阵，清晰的判断了该区的污染程度；在重金属污染程度的影响中，利用了Person相关性系数分析法，合理的说明了元素之间的相关性，并结合元素浓度分布图说明了重金属污染程度的影响。

模型一的缺点：模糊数学在判定元素毒性方面不够准确

模型二的优点：简单的运用数学中的几何规律推到出了一维弥漫方程。

模型二的缺点：假设过多，没有考虑到实际因素带来的结果

模型的改进：对于模型一来说，在利用加权平局来评价元素模糊性的同时，可再利用单因素评价法来进行对比。

参考文献

[1] 姜起源．数学模型（第三版）[M]．北京：高等教育出版社，2009．

[2] 范拴喜甘卓亭李美娟张掌权周旗《土壤重金属污染评价方法进展》中国农学通报 26（17） 310－315 2010

[3] 朱青周生路孙兆金王国梁《两种模糊数学模型在土壤重金属综合污染评价中的应用与比较》环境保护科学 30（3） 53－57 2004

[4] 邓蔚．MATLAB函数速查手册（修订版）[M]。北京：人民邮电出版社，2010.

[5]陈崇希，李国敏，地下水溶质运移理论及模型，武汉：中国地质大学出版社，1996

[6]郭建青，钱会，分析一维砂柱弥散试验数据的反函数法，水利学报，1999，2：43~48

附录

附录一

MATLAB运算程序如下：

xi=[6.2705 0.28996 69.018 49.403 0.093041 18.342 69.106 237.01

7.2514 0.39311 53.409 127.54 0.64236 19.812 93.041 277.93

4.0441 0.15232 38.96 17.317 0.040956 15.454 36.556 73.294

5.708 0.36001 58.054 62.215 0.44682 17.617 63.534 242.85

6.2637 0.28054 43.636 30.192 0.11499 15.29 60.709 154.24];

ai=[1.8 0.07 13 6 0.019 4.7 19 41];

bi=[3.6 0.13 31 13.2

ci=[5.4 0.19 49 20.4

r1=[0 0 1

0 0 1

0 0.205 0.795

0 0 1

0 0 1];

r2=[0 0 1

0 0 1

0 0.0116 0.9884

0 0 1

0 0 1];

r3=[0 0.7533 0.2467

0 0.6283 0.3717

0 0.5578 0.4422

0 0.3782 0.5718

0 0.625 0.375

0 0.585 0.415

0 0.537 0.463

0 0.8466 0.1534];

r4=[0 0 1

0 0 1

0 0.3004 0.6996

0 0 1

0 0 1];

r5=[0 0 1

0.035 12.3

0.051 19.9

31 69];

43 97];

0 0 1

0 0.298 0.702

0 0 1

0 0.6066 0.3934

0 0 1

0 0 1];

for i=1:5

temp=xi(i,:);

for j=1:8

w(i,j)=(temp(:,j)./bi(:,j))./(sum(xi(i,:)./bi(1,:)));

end

R1=w(1,:)*r1;

R2=w(2,:)*r2;

R3=w(3,:)*r3;

R4=w(4,:)*r4;

R5=w(5,:)*r5;

R=[R1;R2;R3;R4;R5]\';%模糊矩阵R

Rf1=R(1,:).*100;

Rf2=R(2,:).*80;

Rf3=R(3,:).*60;

Rf=[Rf1+Rf2+Rf3];%最终得分

z=[7.84 5.93 4.9 6.56 6.35 14.08 8.94 9.62 7.41 8.72 5.93 9.17 5.72 4.49 5.51 11.45 6.14 7.84 7.41 8.5 5.51 9.84 9.39 3.3 4.09 6.14 5.31 3.69 21.87 18.38 10.53 3.5 6.35 5.51 4.49 3.5 5.51 4.29 4.29 6.56 16.58 7.41 5.93 4.69 4.9 5.31 4.29 5.51 4.69 7.2 5.31 4.9 4.9 3.89 3.69 3.11 3.89 3.89 2.91 3.3 4.9 4.09 2.91 2.72 3.11 3.3 3.3 6.14 3.69 4.49 3.11 8.06 3.69 3.69 3.5 2.72 8.5 1.77 2.53 3.69 6.14 10.99 6.35 30.13 3.89 2.91 1.96 6.98 2.91 5.93 5.93 7.41 4.29 4.69 5.72 6.77 4.49 3.11 2.91 4.29 7.63 5.93 2.34 2.91 5.72 2.34 6.56 4.69 6.35 5.1 4.69 3.5 4.69 4.49 3.3 2.91 4.09 5.72 4.9 4.9 4.09 5.93 2.91 2.72 2.34 2.53 3.89 2.34 2.34 1.96 2.72 3.11 3.5 2.72 1.77 3.5 2.53 9.17 7.41 5.72 8.06 5.72 8.72 6.77 6.56 6.77 6.98 6.14 6.69 8 6.69 8.23 9.35 8.9 3.77 5.41 7.78 5.62 5.41 4.58 6.91 5 5.62 6.91 7.78 6.26 7.56 4.79 2.77 6.26 7.34 4.17 5 7.56 5.62 6.05 4.79 23.72 6.47 4.17 5 4.17 6.26 5 1.61 4.58 5.41 3.57 2 3.37 6.47 3.37 6.91 2.38 7.56 8.67 4.79 5.41 5.83 5.62 4.38 2.77 5.41 5.62 6.69 7.56 2.77

1.61 5.62 2.97 4.58 5.2 3.17 5 7.56 5.2 6.05 2.77 2.38 5.2 8.67 5.41 7.12 4.58 8.23 4.58 6.05 8.23 4.17 3.97

6.69 4.17

9.13 7.34

4.38 5.62

3.77 5.41

5.62 3.77

3.37 3.17

3.17 2.97

8.23 10.74

4.79 7.56

x=[74 1373 1321

4043 2427

2486 3299

4948 5567

8077 8017

10631 10685

11902 13244

13920 14844

12778 17044

18134 17198

22304 21418

21766 22674

25998 27177

26015 27700

22965 23198

15007 3518

5664 5541

5503 5636

8403 8079

11557 11415

13855 14862

17734 16823

18438 18556

19569 19411

22624 21703

8866 9296

7106 6423

11646 12641

15140 15198

14269 13277

11121 10856

14262 14624

10.74 3.77 10.27 5.2 6.47 6.47 8.9 3.378.23 4.38 5.2 5.2 5.41 5 9.58 3.17 3.775.62 4.79 5.2 3.37 7.34 3.57 4.38 5.416.05 6.26 5.2 6.91 4.58 8.67 6.47 7.1210.97 9.81 8.23 5.41 2.77 7.78 6.473.57 4.38 6.91 7.12 4.38 3.97 8 3.572.97 7.12 5.62 5.83 3.97 3.17 2.773.17 1.8 3.17 3.77 4.79 3.57 6.47 7.3411.68 7.34 6.05 5.41 6.26 6.47 6.479.35];

0 1049 1647 2883 2383 2708 2933 42333526 5062 4777 5868 6534 5481 45923573 4741 5375 5635 5394 5291 47427004 7304 7048 8180 9328 9090 80496869 7056 7747 8457 9460 9062 931910643 11702 11730 11482 10700 10630 1167812746 12855 13797 14325 15467 12442 1309316569 16387 16061 15658 14298 14177 1509217087 17075 17962 18413 19007 18738 1781417144 18393 19767 21006 21091 22846 2366421439 20554 20101 21072 20215 18993 1996822535 25221 26453 26416 27816 25361 2406526424 26073 24631 24702 25461 24813 2608627696 27346 26591 27823 27232 24580 2415324685 28654 24003 21684 22193 17079 152553469 3762 3927 4153 3267 4684 54955451 4020 4026 5101 5438 5382 53146605 7093 7100 6837 7906 8045 83949663 9469 9178 9095 10225 10210 1034011649 12734 12696 12400 12591 13765 1369414896 15387 15810 16032 15801 15087 1687217008 17203 17005 16947 16301 17904 1830318954 18012 19072 20282 21475 21450 2026119501 20582 19909 21018 22176 23359 232385006 5734 6395 7405 8446 7612 79129475 9212 8629 7776 8622 9237 83077458 8904 10547 10398 10395 11529 1156314000 14207 14065 12734 12727 14173 1546715248 16428 16289 16267 16440 16440 1541213175 12153 11958 10800 10022 9333 927712644 12625 9036 10599 12632 14405 1407416629 18470 20591 20983 20177 19041 18906

18467 17414 15748 15517 16607 15952 22605 23146 22046 23785 25981 27380 25021 23325 26852 17981 14482 14318 10352 9095 10510 13954 10142 17765 6924 4678 6182 5985 7653];

y=[781 731 1791 1787 2127 2728 3617 3692 2295 1767 895 1895 3971 4357 4339 4897 4904 5641 6004 4603 5999 6018 6213 6434 8643 7965 8631 7349 7293 7293 6782 6226 5230 4600 4496 4311 5365 5439 6401 7210 7286 8348 8260 8991 8311 7639 6799 6472 5528 4472 4480 5532 6354 8184 8774 8618 7709 7056 8450 8945 9621 8666 8658 4329 4339 5354 5519 6055 6609 7352 7594 7418 6684 6936 5799 10691 11933 12924 12823 11721 11488 10921 10707 10046 9810 9081 9183 8810 8819 9482 9149 9790 10527 10721 11383 11228 10774 10404 9951 12371 12961 12348 12173 11293 5795 5577 6508 5581 6423 7353 7032 7771 8639 8807 9422 9522 9834 10799 11094 12078 11609 11621 13331 13715 14737 14482 13319 12450 13535 13523 14278 8755 15286 13101 12185 5894 5110 5535 2571 2308 2170 2110 2299 793 1364 1205 1653 2093 2757 2990 3913 4080 3994 3012 2060 1127 133 374 1381 2449 3490 3978 3052 2035 1075 0 1288 2286 3299 3975 3821 2789 1764 1581 2585 3515 4015 3024 2060 1063 1353 2357 3345 2524 1603 729 2307 3061 3966 3512 2798 3629 4207 4775 6218 7212 7487 8299 8287 7385 6539 5588 4874 4414 8519 8590 8540 7555 7586 7348 6934 6091 6548 5300 5764 5492 5325 6502 4818 6591 8846 9659 10443 10981 11200 11938 12840 13143 13102 12000 11305 12086 10613 10638 9872 9726 9467 8831 8920 8868 9591 10360 11203 11243 10298 9381 9560 8970 9980 10987 10344 7691 11941 12080 11101 10100 9106 9069 10072 11058 12068 13232 12982 12877 13204 12238 12336 13313 13282 12204 14631 16148 16432 14727 14943 16259 17538 17980 17949 18032 16516 15129 14004 14481 14411 13549 15862 17642 15769 16346 17001 15476 15728 17034 17365 18397 14301 15382 17634 17643 18051 18202 16290 16701 16114 18449 12692 13569 17133 16414 15314 5615 1662 3561 5696 3765 2005 2567 1952];

x1=[4043 2427 4777 6534 4592 2486 3573 5375 7304 9328 8077 7056 15467 14844 16569 14298 21418 26453 26416 5101 5382 5503 5636 6605 11649 12591 13855 14862 15387 15810 17203 16301 17904 10395 15467 16428 16289 12153 11958 10800 9277 11121 12625 17981];

y1=[1895 3971 4897 5641 4603 5999 6213 8643 5230 4311 6401 8348 8658 5519 6055 7418 10721 5577 6508 4080 3012 1127 133 374 3515 1063 3345 2524 729 2307 6218 8299 8287 11203 12080 9069 10072 12336 13313 13282 16148 16432 16259 18449];

x2=[0 1647 2383 4742 4948 5567 22674 5438 12734 14896 16947 5006 6395 7405 8446 7612 8866 9475 9212 8629 7776 7106 6423 7458 11646 12641 14000 14207 15140 16440 10022 9333 10856 12644 9036 10599];

y2=[1787 2728 3692 7293 7293 6782 12173 3994 4015 1603 7487 8846 10443 10981 11200 11938 13143 12000 11305 12086 10613 9467 8831 8920 9381 9560 8970 9980 11101 13232 12204 14631 14727 14943 17538 17980];

x3=[14325 17044 18413 19007 18738 23664 18993 19968 22535 27816 25361 24065 25998 27177 26073 24631 24702 25461 26086 26015 27700 27696 27346 26591 27823 27232 24580 24153 22965 24685 28654 24003 21684 22193 17079 15255 15007 15801 17008 21475 20261 19569 19411 23359 23238 22624 21703 12734 14405 14074 14262 20591 20983 20177 18906 18467 16607 15952 22605 23146 22046 23785 25981 27380 23325 26852];

y3=[8666 10691 11721 11488 10921 9790 12371 12961 11293 5581 6423 7353 7032 7771 8807 9422 9522 9834 11094 12078 11609 11621 13331 13715 14737 14482 13319 12450 13535 14278 8755 15286 13101 12185 5894 5110 5535 3966 4775 8540 7586 7348 6934 5325 6502 4818 6591 10344 18032 16516 15129 13549 15862 17642 16346 17001 17365 18397 14301 15382 17634 17643 18051 18202 16701 16114];

x4=[74 1373 1321 1049 2883 2708 2933 3526 5062 5868 5481 3299 5635 5394 5291 7004 7048 8180 9090 8049 8017 6869 7747 8457 9460 9062 9319 10631 10685 10643 11702 11730 11482 10700 10630 11678 11902 13244 12746 12855 13797 12442 13093 13920 15658 14177 12778 17087 17075 17962 17814 18134 17198 17144 18393 19767 21006 21091 22846 22304 21439 20554 20101 21072 20215 21766 25221 26424 24813 23198 4020 4026 5314 7093 7100 6837 7906 8045 8394 8403 8079 9663 9469 9178 9095 10225 10210 10340 11557 11415 12696 12400 13765 13694 16032 16872 17734 17005 18438 18954 18012 19072 20282 21450 19501 19909 21018 22176 5734 7912 9296 8622 9237 8307 8904 10547 10398 11529 11563 14065 12727 15198 15248 16267 16440 15412 14269 13277 13175 12632 14624 16629 18470 19041 17414 15748

25021 5985];

y4=[781 731 1791 2127 3617 2295 1767 4357 4339 4904 6004 6018 7965 8631 7349 6226 4600 4496 5365 5439 7210 7286 8260 8991 8311 7639 6799 6472 5528 4472 4480 5532 6354 8184 8774 8618 7709 7056 8450 8945 9621 4329 4339 5354 7594 6684 5799 11933 12924 12823 10707 10046 9810 9081 9183 8810 8819 9482 9149 10527 11383 11228 10774 10404 9951 12348 5795 8639 10799 13523 2990 3913 2060 1381 2449 3490 3978 3052 2035 1075 0 1288 2286 3299 3975 3821 2789 1764 1581 2585 3024 2060 1353 2357 3061 2798 3629 7212 6539 4874 4414 8519 8590 7555 6091 5300 5764 5492 9659 12840 13102 10638 9872 9726 8868 9591 10360 11243 10298 10987 7691 10100 9106 11058 12068 12982 12877 13204 12238 17949 14004 14481 14411 15769 15476 15728 16290 2567];

x5=[4233 4741 16387 16061 15092 3518 3469 3762 3927 4153 3267 4684 5495 5664 5541 5451 15087 16823 18303 18556 20582 14173 15517 14482 14318 10352 9095 10510 13954 10142 17765 6924 4678 6182 7653];

y5=[895 6434 6609 7352 6936 2571 2308 2170 2110 2299 793 1364 1205 1653 2093 2757 3512 4207 7385 5588 6548 11941 17034 12692 13569 17133 16414 15314 5615 1662 3561 5696 3765 2005 1952];

xi=linspace(min(x),max(x),319);

yi=linspace(min(y),max(y),319);

[xi,yi]=meshgrid(xi,yi);

zi=griddata(x,y,z,xi,yi,\'v4\');

contour(xi,yi,zi,30)

contourf(xi,yi,zi,30)

colorbar

hold on

scatter(x1,y1,\'.\')%生活区坐标

hold on

scatter(x2,y2,\'o\')%工业区坐标

hold on

scatter(x3,y3,\'x\')%山区坐标

hold on

scatter(x4,y4,\'+\')%交通区坐标

hold on

scatter(x5,y5,\'*\')%公园绿地区

xlabel(\'x(m)\');

ylabel(\'y(m)\');

更多推荐

污染,重金属,土壤

城区表层土壤重金属污染分析模型南开大学数学建模

发布评论取消回复

最近发表

热门文章

标签列表

城区表层土壤重金属污染分析模型南开大学数学建模

相关文章

发布评论取消回复

最近发表

热门文章

标签列表