2011数学建模A题 城市表层土壤重金属污染分析

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

摘要

对于问题一，首先采用克里格插值法根据附件中给的采样点的浓度数据对城区内每种重金属浓度值进行插值，进而绘制每种重金属的浓度的空间分布图，对此为了方便计算我们借助于Surfer软件进行绘图，然后我们采用单因子指数法和内梅罗综合指数法对各功能区的污染程度进行评价，但是这两种方法无法从自然异常中分离出人为异常，为了弥补其不足，采用地累积指数法[1]对土壤重金属污染做进一步评价，继而得到不同区域重金属的污染程度。评价结果如表9所示。

对于问题二采用主成分分析法，旨在利用降维的思想，把多指标转化为少数几个综合指标，根据主成分分析法的一般步骤，首先对附件中给出的重金属浓度进行标准化处理，然后得到各金属之间的相关系数矩阵，求出相关系数矩阵的特征值和向量值，再得到因子的成分矩阵，确定出主成分的个数，计算出各因子的成分得分矩阵，最后通过算综合主成分中各个因子的权重系数得到污染性较大的因子，最后分析该污染主要原因。结果如表14和评价结果。

问题三我们用对流-扩散偏微分方程来进行描述，对流扩散方程是描述粘性流体运动的非线性偏微分方程模型，我们将对流扩散方程进行简化，即变为二阶椭圆形偏微分方程，利用有限插值数值法估计出污染物的浓度分布，并与实际检测值相比较，偏差较小的即为污染源的位置及源强，即将问题转化为非线性最优化问题，结果如表所示。

对于问题四在分析问题三模型的优缺点基础上，为更好的研究城市地质环境的演变模式，将问题三的模型中二维对流扩散偏微分方程扩展到三维，这样就

可以全面考虑影响模型准确性的各相关参数，通过收集不同的地理，天气条件下地质元素的空间分布信息，进而根据这些信息，建立优化模型，即三维方程的参数，在此基础上通过模拟仿真进行分析。

关键字：单因子指数法内梅罗综合指数法地累积指数法主成分分析法

偏微分模型

一问题重述

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

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

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

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

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

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

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

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

(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？

二问题分析

对于问题一，要根据问题附件中给出的采样点以及在采样点处8种主要重金属元素的信息，研究各金属元素的空间分布，并分析不同功能区的污染程度，附件中只给出了金属元素在采样点处的浓度，要给出在整个城区内的空间分布，需要通过插值得到更加密集的浓度分布值，由于城区内样本点间的重金属浓度存在空间相关性，因此可以通过克立格法进行插值。在这里我们借助Surfer8.0软件，也是通过克里格插值的方法生成重金属含量数据，并绘制等值线图。由于还需要分析不同区域重金属的污染程度，污染程度需要抽象为污染指数，通过比较污染指数与国家背景值，确定污染程度较大的金属元素，针对这个问题我们采用单因子指数法和内梅罗综合指数法[1]，但是这两种方法无法从自然异常中分离出人为异常，为了弥补其不足，采用地累积指数法[3]对土壤重金属污染做进一步评价，继而得到不同区域重金属的污染程度。

对于问题二要求通过数据说明重金属污染的主要原因，现有资料表明，某些重金属空间分布具有相关性，相关性较大的重金属可能在成因和来源上有一定的联系。在此我们选用主成分分析法来说明。主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指标。在实证问题研究中，为了全面、系统地分析问题，我们必须考虑众多影响因素。这些涉及的因素一般称为指标，在多元统计分析中也称为变量。因为每个变量都在不同程度上反映了所研究问题的某些信息，并且指标之间彼此有一定的相关性，因而所得的统计数据反映

的信息在一定程度上有重叠。在用统计方法研究多变量问题时，变量太多会增加计算量和增加分析问题的复杂性，人们希望在进行定量分析的过程中，涉及的变量较少，得到的信息量较多。

根据主成分分析法的一般步骤，首先对附件中给出的重金属浓度进行标准化处理，然后得到各金属之间的相关系数矩阵，求出相关系数矩阵的特征值和向量值，再得到因子的成分矩阵，确定出主成分的个数，计算出各因子的成分得分矩阵，最后通过算综合主成分中各个因子的权重系数得到污染性较大的因子，最后分析该污染主要原因。

对于问题三要求分析重金属污染物的传播特征，建立模型并分析污染源的位置，对此我们用对流-扩散偏微分方程来进行描述，对流扩散方程是描述粘性流体运动的非线性偏微分方程模型，可以刻画很多自然现象，如污染物的扩散、降解、流体流动与传热等，我们将对流扩散方程进行简化，即变为二阶椭圆形偏微分方程，利用有限插值数值法估计出污染物的浓度分布，并与实际检测值相比较，偏差较小的即为污染源的位置及源强，即将问题转化为非线性最优化问题。

对于问题四在分析问题三模型的优缺点基础上，为更好的研究城市地质环境的演变模式，应该确定收集什么信息，有了这些信息如何建立模型来解决问题。首先将问题三的模型中二维对流扩散偏微分方程扩展到三维，这样就可以全面考虑影响模型准确性的各相关参数，为了得到参数的准确表达式，通过收集不同的地理，天气条件下地质元素的空间分布信息，进而根据这些信息，建立优化模型，即三维方程的参数，在此基础上通过模拟仿真进行分析。

三问题假设

1.问题一中国家评价标准采用二级标准。

2.假设全部采样点在自然情况下产生，排除其他因素的干扰。

3.假设该城区没有受到灾难性的污染。

4.观测值能够反映一段时间内保持稳定的浓度分布。

5.假设污染源是均匀向外扩散的。

四 符号说明

i 污染物类别

i C 土壤中污染物实测浓度 i S 土壤中污染物的背景值 i P 土壤单因子污染指数

P 为综合污染指数 ave I 元素污染指数平均值

m ax I 元素污染指数中的最大值 Igeo 地累积指数

n C 各个区域元素n 的平均值 n B 元素n 地球化学背景值

k 各地岩石差异导致的背景值变动系数 ij C 第i 个采样点的第j 种重金属的浓度值

j C -

第j 种重金属的均值

j δ 第j 种重金属的标准差

xy r

二元素相关系数

五 建模与求解

5.1问题一建模与求解

为了得到更加形象具体的重金属空间分布图，在这里我们借助Surfer8.0软件，也是通过克里格插值的方法生成重金属含量数据，并绘制等值线图。

图1 城区地形三维图

图2各功能区地形图

图3 As城区浓度分布图

图4 Cd城区浓度分布图

图5 Cr城区浓度分布图

图6 Cu城区浓度分布图

图7 Hg城区浓度分布图

图8 Ni 城区浓度分布图

图9 Pb 城区浓度分布图

图10 Zn 城区浓度分布图 5.1.1单因子指数法

我们对该城区不同区域的污染程度进行评价，由于每个区域都存在八种重金属，所以我们对每个区域的八种重金属含量分别进行评价，评价模式分为单因子指数法和内梅罗综合指数法。首先建立单因子污染指数式为：

i

i i C P S

（1）

式中i P 为土壤中污染物i 的污染分指数，i C 为土壤中污染物实测浓度，i S 为土壤中污染物的背景值。然后通过求解上述单项污染指数式得到五个区域i 种污染元素的分指数，如表1（表2、表3、表4、表5）所示为生活区污染分指数。 表1 各区域重金属实测平均浓度

表2 各区域单因子污染分指数

表3 土壤单因子污染指数评价标准

5.1.2内梅罗综合指数法

单项污染指数得到每个区域各种重金属污染程度，为了得到每个区域的综合污染程度，我们利用内梅罗综合指数法求解。内梅罗综合指数式为：

2

)(2

max 2I I P ave +=

（2）

式中：P 为综合污染指数, ave I 为各元素污染指数平均值, m ax I 为各元素污染指数中的最大值。

通过求解内梅罗综合指数式得到每个区域的综合污染指数（见表7，表8）

表4 各区域综合污染指数

表5 土壤内梅罗综合污染指数评价标准

5.1.3地累积指数法

虽然单因子指数法和内梅罗指数法均能对土壤重金属污染程度进行较为全面的评价，但无法从自然异常中分离人为异常，然而地累积指数法弥补了这项不足。地累积指数法表达式为：

(

)2log n n C Igeo k B ??

=?????

（3）

式中：Igeo 为地累积指数，n C 是各个区域元素n 的平均值，n B 是元素n 地球化学背景值，k 为各地岩石差异导致的背景值变动系数（一般取值为1.5）。 通过对地累积指数式的求解得到每个区域各种重金属的污染程度（见表9）

表6 每个区域各种重金属的污染程度

评价结果

经过内梅罗综合指数法和地累积指数法对该城不同区域重金属的污染程度分析，得到生活区中As、Cd、Cr、Hg、Pb元素属于轻度污染，Ni元素成为无污染元素，Cu、Zn污染为中度污染。综合污染等级为警戒级。工业区中As、Ni、Cr元素属于轻度污染，Cu、Pb、Zn、Cd元素成为中度污染，Hg 则属于强度污染。综合污染等级为严重污染。山区As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种元素基本都对当地环境无污染，综合污染等级为清洁级。主干道路区As、Cd、Cr、Cu、Ni、Pb、Zn等元素基本为轻度污染，只有Hg为强度污染，综合污染等级为轻度污染。公园绿地区大部分元素为轻度污染，有很少地方存在中度污染，综合污染等级为清洁级。

5.2问题二的建模与求解

根据问题的分析，某些重金属空间分布的含量也具有一定的相关性，相关性较大的金属可能在成因和来源上有一定的联系。因此，需要对问题中所给出的8种重金属元素之间的关系进行分析。在这里我们选用主成分分析法来求解这个问

题。

主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指

标。在实证问题研究中，为了全面、系统地分析问题，我们必须考虑众多影响因素。这些涉及的因素一般称为指标，在多元统计分析中也称为变量。因为每个变量都在不同程度上反映了所研究问题的某些信息，并且指标之间彼此有一定的相关性，因而所得的统计数据反映的信息在一定程度上有重叠。在用统计方法研究多变量问题时，变量太多会增加计算量和增加分析问题的复杂性，人们希望在进行定量分析的过程中，涉及的变量较少，得到的信息量较多。

5.2.1对原始数据进行标准化

对原始数据标准化处理的一般公式为：

)8,...,2,1,319,...,2,1(/'

==??

?

??-=-

j i C C C j

j ij ij δ （4）

其中ij C 为第i 个采样点的第j 种重金属的浓度值，j j C δ,-

分别表示第j 种重金属的均值和标准差。

5.2.2求重金属含量相关系数矩阵

对于任意的两种重金属元素X 与元素Y 在第i 个样本中的含量值分别为i x 和

i y （n i ,...,2,1=），则二元素相关系数的一般计算公式为

yy

xx xy n

i n

i i

i

n

i i

i xy S S S y y x x y y x x r =

??

?

?

?

-???

??

-??? ??-??? ??-=

∑∑∑==-

-

=-

-1

1

2

21 （5）

其中 ∑∑∑∑===-=--=??? ??-??? ?

?-=n i i n i i n

i i i i n

i i xy y x n y x y y x x S 11111

2

111

22

1??? ??-=???

??-=∑∑∑===-

n i i n

i n

i i i xx x n x x x S

2

11

22

11??? ??-=???

??-=∑∑∑===-n i i n i i n

i i yy y n y y y S

式中-

-y x ,分别为元素X 与元素Y 样本点观测值的算术平均值。

首先我们通过将生活区的重金属元素原始数据导入SPSS 软件进行主成分分析处理可以得到各重金属元素之间的相关系数如下表所示。（以生活区为例） 表17 各个因子间的相关系数矩阵

5.2.3求相关系数矩阵的特征值和特征向量

利用相关系数矩阵可求出相应的因子特征值和累计贡献率，用SPSS 统计软件计算可得到生活区的各因子特征值和成分矩阵结果如下表所示。

表8 解释的总方差

表9 成分矩阵（载荷矩阵）

主成分个数提取原则为主成分对应的特征值大于1的前m个主成分，通过表21，可知提取3个主成分。用表16中数据除以主成分对应特征值便得到三个成分中每个指标对应的系数，即特征向量A1，A2，A3，将特征向量与标准化后数据相乘得到主成分表达式：F1=A1*ZX1; F2=A2*ZX2; F3=A3*ZX3。为了方便计算我们在这里用SPSS解决这个问题得到成分得分系数矩阵。

表10 各因子得分系数矩阵

以每个主成分所对应的特征值占所提取主成分总的特征值之和的比例作为权重计算主成分综合模型，即用第一主成分F1中每个指标所对应的系数乘上第一主成分F1所对应的贡献率再除以所提取的三个主成分贡献率之和，然后加上第二主成分F2对应系数乘以其贡献率百分比，再加上第三主成分F3对应系数乘以其

贡献率百分比，就得到了综合主成分F中每个指标所对应的权重系数A。根据权重系数A值大小降序排列，系数大的对应的重金属元素就是主要污染的主要原因。这里只将生活区列举出来，其它区域分析方法与上述方法相同，就不做过多叙述，其他区域表格见附录。

表11 各区域因子的权重系数

结果评价

对该城区不同区域采集的土壤As、Cd、Cr、Cu、Hg、Ni、Pb、Zn元素进行主成分分析，通过主成分分析法得到:

1> 生活区元素Zn、Cr在该城市土壤中的含量明显高于国家背景值。这是由于生活区存在大量的公路，公路两侧含铅汽油的燃烧和汽车轮胎的磨损的粉尘会增加土壤中Zn、Cr等元素的含量，而且由于生活区人类活动比较频繁，日常生活用品丢弃后成为垃圾也会导致这三种元素的增加。

2> 工业区主要污染元素为As、Cd、Cr、Ni、Pb、Zn。因为工业区“三废”排放，采矿和冶炼会增加这些元素在土壤中的含量。

3> 山区全部元素几乎都在国家背景值以下，只有Cu元素为主要污染元素。这是由于山区远离城市和工业区的污染源，只有一部分金属矿山的开采会导致一些金属元素的外露。

4> 主干道路区Pb、Cd元素都高于背景值，其余元素都低于背景值。机动车尾气排放既是城市大气的主要污染源，也显著引起公路两侧土壤的重金属污染，汽车汽油、发动机、轮胎、润滑油和镀金部分都能燃烧或磨损而释放出Hg、Zn元素。

5> 公园绿地区中Pb为主要重金属元素，Hg、Cu、Zn为次要重金属元素。公园绿地区农药化肥的使用会不同程度的影响重金属的污染，而且地下水灌溉和塑料薄膜会增加、Cu、Zn元素的土壤含量。

不难发现每个区域都受到不同程度的Cd元素的污染，这是因为镉元素主要以硫化镉形式储存于锌矿、铅锌矿和铜铅锌矿中，土壤镉主要来源于锌矿、铅矿的冶炼、合金、电镀、化工厂等排放的废水，工业固废堆放，含有废旧电池的生活垃圾渗滤，污泥施肥以及过量或不恰当地使用化肥农药等。

总体来说：工业化程度越高的地区污染越严重，市区高于远郊和农村，地表高于地下，污染区污染时间越长重金属积累就越多，以大气传播媒介土壤重金属污染土壤具有很强的叠加性,熟化程度越高重金属含量就越高。

5.3 问题三模型建立与求解

通过对问题二的分析，我们知道了重金属的传播途径，比如重金属b P 和d C 的污染主要为汽车尾气排放所致，g H 的污染为燃煤引起，它们都是经过空气来传播的，为此用对流-扩散偏微分方程来描述。n Z 主要是由于工业废水排放污染，经过地表径流的传播，则可以应用对流-扩散偏微分方程来进行描述。对流-扩散方程是描述粘性流体运动的非线性Burgers 方程的线性化偏微分方程模型，能够刻画很多自然现象，如污染物的扩散，降解，流体流动与传热和电化学反应等。 首先建立对流-扩散偏微分方程传播模型。然后，依据问题所给的条件，可以知道这是对流-扩散偏微分方程的反问题。在求解时，将对流-扩散偏微分方程进行简化变为二阶椭圆型偏微分方程。利用椭圆型偏微分方程的有限差分数值解法正向估计出污染物的浓度分布，再代入不同的位置和强度参数，与实际检测值相比较，偏差最小的即为污染源的位置及源强，即将问题转化为非线性最优化问题，利用matlab 中ool p det 工具箱求解椭圆型微分方程，则可得到各重金属元素的污染源位置及源强。

5.3.1 对流-扩散偏微分方程模型的建立 一般的对流-扩散偏微分方程模型为

??????

???Ω?∈=??>Ω?∈=Ω∈--+??????????+??????????=??+??+∑=.2111

),(),,(,0,),(),,,(),,(,),()

,(),(),(),(y x y x f n

c

t y x t y x f t y x c y x y y x x s y c y x e y x c y x e x y c v x c u c y x k q

i i i i δ（6） 其中Ω为有界区域，n 为边界上的外法向量，),(y x u 为x 轴方向的流速，)

,(y x v

为y 轴方向的流速，),(y x e 为扩散系数，),(y x k 为污染物的降解率，),(y x c 为污

染物的浓度分布，）（i i y x ,和),3,2,1(q i s i =分别表示多个污染点和污染强度，δ为

狄拉克函数，即单位脉冲函数。

5.3.2对流-扩散偏微分方程反问题的求解

若已知),(y x c 的分布，并且已知,,v u ),(y x k 和),(y x e 那么源项识别的反问题就是根据这些已知的分布确定源项∑=--q

i i i i y y x x s 1,)(δ，即确定污染源的位置

）（i i y x ,和污染强度i s 。

根据模型的假设，观测值为一段时间内稳定的浓度分布，即二阶稳态对流扩散方程。根据相关文献可知扩散系数),(y x e =-410）（s m /2，污染物降解率

),(y x k =0.001

）（1-s ，暂时不考虑横纵轴的流速，使x 轴方向的流速),(y x u =0， y 轴方向的流速),(y x v =0.于是则可将上述的偏微分方程简化为

∑=--=+???

?

????+??q

i i i i y y x x s c y c x c 1222),(001.0-δ

（7）

该方程为二阶椭圆型偏微分方程。要在污染区域中大概找出污染范围，首先选择

污染源初始位置）（∧

∧

00,y x ，以及污染强度∧

0s 。利用椭圆型偏微分方程的有限差分

数值解正向估计出污染物的浓度分布),(y x c ∧

，

Ω∈）（y x ,。 ),(y x c ∧的分布与初始位置）（∧∧00,y x ，以及污染强度∧

0s 的取值有关。用函数),,,,(000i i k k k y x s y x g ∧

∧

∧

∧

∧

表示在初始位置）（∧

∧

00,y x ，以及污染强度∧

0s 的条件下测得的

污染物浓度分布，0k x ∧，0k y ∧，0k s ∧

为第k 次估计值。将所得到的污染物浓度与实际检测值相比较，差距最小的即为最接近污染源的位置。即将问题转换为非线性最优化问题：