2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

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

摘要

随着经济的发展，人类活动对城市地质环境的影响越来越大，因此对地质环境的研究分析非常重要。本文通过对某城区不同功能区域污染程度进行划分，分析了重金属污染的主要原因，并建立确定该城区不同重金属污染物来源的模型，最后给出为使研究更加深入而需要搜集的其他影响城市地质环境的信息，建立模型分析主要因素，以便更好研究城市地质环境演变模式。

针对问题一，题目要求给出8种重金属元素在该城区的空间分布。本文以采样点所在位置的坐标分别为x、y轴，以该采样点各种元素的浓度为z轴，利用matlab样条工具箱进行“v4”内插，作图得到8种重金属元素的水平面分布图，并将三维图像转化为等值线图，从图4-图11可以直观地看出各种重金属元素在空间的分布情况。结合图像，采用单因子指数法和内梅罗污染指数法对城区内不同区域重金属元素的污染程度进行评估，得到污染程度从高到低依次为：工业区，主干道路区，生活区，公园绿地区，山区。

针对问题二，题目要求分析说明重金属污染的主要原因。灰色关联分析法是分析系统原因的一个行之有效的方法。本文通过对各功能区中各重金属元素与内梅罗指数进行灰色关联分析，得到各重金属元素与重金属污染程度的灰色关联度，并对关联度大的重金属元素进行分析，分析这些元素的用途、来源、特性等，进而得到各区域中造成重金属污染的主要原因是：冶金业在该地区造成工业污染以及农业上滥用农药。

针对问题三，题目要求建立模型确定污染源的位置。首先从纵向，横向以及考虑自然因素的情况下分析重金属污染物的传播特征，详细情况如模型分析所示。针对这些传播特征，主要考虑对传播路径起主要影响的纵向和横向因素，运用类比的方法，把重金属污染物的传播方式等效为电磁场理论中电场能量的传递方式，即用通量的大小来体现某特定采样点是否为源的可能性。同时给出了线通量的定义及在本问题中的表达式，得到各元素线通量分布及等通量线，并得到各种元素可能的污染源坐标，如表5所示。 针对问题四，题目要求在收集到更多信息的前提下，建立模型以更好地研究城市地质环境的演变模式。本文认为，需要收集的信息还有：气温、降水量、蒸发量、水体污染情况、地面坡度、岩性特征、历史地震活动情况、工程地质。结合在问题一中已获得的土壤污染综合指数，采用主成分分析法建立模型，得到影响城市地质环境的主要因素，从而通过研究这些影响地质环境的主要因素获得城市地质环境的演变模式。

关键词：关键词：重金属元素污染 内梅罗指数 灰色关联分析 线通量模型 主成分分析

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

一、问题重述

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

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

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

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

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

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

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

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

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

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

二、符号说明

三、模型假设

1、假设所采用的采样方法科学合理，所得数据准确。

2、假设污染物的扩散只与浓度有关，且顺浓度梯度方向扩散，无外力参与。

3、假设不同功能区对应一种级别的污染程度，不存在具有相同污染程度的功能区。

4、假设数据是在污染物扩散稳定的情况下测得的。

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

四、问题分析

4.1 对重金属元素的空间分布和在不同城区的污染程度进行分析

8种重金属元素在空间的分布，用等高线图可以直观地表示出来。本文以采样点所在位置的坐标为x、y轴，用暖色系（灰度高）表示重金属含量高，以冷色系（灰度低）表示重金属含量低。由于题目所给采样点为离散的，故本文在画图前对数据进行了“v4”

[1,2]内插处理，利用matlab样条工具箱得到8种重金属元素在该城区的8个空间分布图。 8种重金属元素在该城区各区域的分布各不相同，要得到该城区内不同区域重金属的综合污染程度，还需要一个指标将8种不同的重金属元素的污染程度综合起来。内梅罗污染指数[3]是一种兼顾极值或称为突出最大值的计权型多因子环境质量指数，能够合理有效地评价环境的综合污染程度，内梅罗指数越大说明该地区污染越严重。

4.2 对重金属污染的主要原因进行分析

重金属污染是指由于人类活动，土壤中的微量有害元素在土壤中的含量超过背景值，过量沉积而引起的含量过高，统称为土壤重金属污染[4]。在不同的功能区，重金属的由来有所不同。本文采用灰色关联分析法[5]，通过对各功能区各重金属元素进行关联度分析,得到对每个功能区重金属污染贡献最大的几种元素，并针对功能区的不同特点分析贡献最大的几种重金属元素的由来，从而可以得到各功能区中重金属污染的主要原因。

4.3 对污染源位置确定进行分析

金属污染物的传播路径受到多种因素的影响，如风力，水流，海拔高度……由某一污染源散发出的重金属污染物，主要有以下传播特征：

（1）从纵向的角度上看，由于重力的作用，重金属污染物从高处向低处传播，在没有污染源进行补充的情况下，污染物浓度随海拔的减小而下降；

（2）从横向的角度上看，重金属污染物从污染源向四周呈放射性传播，在同一海拔高度，污染物浓度随与污染源的距离的增加而下降；

（3）考虑自然因素的影响，由于重金属污染物还存在于大气或者水流中，其传播方向与风向和河流等流向也有一定的相关关系。

在信息有限的情况下，我们只知道采样点所在位置的地理坐标以及该点处各种元素的浓度值，而不知道各点的气象情况和其他对浓度产生影响的因素。因此，想要通过已知数据分析出重金属元素的传播路线以确定源的所在位置，是十分困难的。

在问题一的分析中，我们可以得到重金属元素在该地区的分布等浓度线，然而，仅仅观察其分布图，认为浓度大的地方就是污染源是十分不科学的。因为某点的污染物浓度是重金属元素在一段时间内的积累效应，而对于一个特定的采样区域，污染物既有注入又有流失。例如：峡谷地带因为水的聚积作用，水流带来的污染物在此聚积，因此山谷的污染物浓度虽高，注入却远大于流失，显然不是源之所在。在此本文采用类比的方法，并忽略海拔的影响，此处的等浓度线可以等效为电磁场理论[6]中的电场等势面，重金属元素的传播路径等效为电场线。根据电磁场中的通量理论，污染源作为通量源，在取相同的包络面的情况下，源所在的地方其通量必定取得最大正值，且同时考虑了注入和流失两个方面的影响。按照以上原理，对该城区的各个网格点计算污染物的通量值，绘制相应的分布和等通量线，就能直观观察出源所在的大致位置。

4.4 对模型三以及影响地质环境的信息进行分析

对于模型三，本文通过建立通量模型确定了污染源的位置。该模型能够通过各采样点的重金属元素的线通量判断出该城区各种重金属元素的主要来源。模型类比电场通量模型，形象且容易理解。但由于题目所给信息不足，本文在建立模型时对数据进行了简

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

化处理，将离散无规律的采样点简化成规律的采样点，并且通过插值法对该城区中未采集到的点进行插值，由于该模型建立在等间距网格点上，故与实际采样点有一定的误差。

城市地质环境除了受土壤因素的影响外，还受到人类活动、气候、自然灾害等因素

[ 7]的影响，因此，为了更好地研究城市地质环境的演变模式，还应该收集的信息有：气

温、降水量、蒸发量、水体污染情况、地面坡度、岩性特征、历史地震活动情况、工程地质。有了以上的信息，通过对这些信息进行主成分分析，得到影响城市地质环境的主要因素，从而可以通过研究这些因素的情况而获得城市地质环境的演变模式。

五、模型建立

5.1 模型一的建立

该模型采用单因子指数法和内梅罗污染指数法对土壤重金属污染程度进行评价。

5.1.1 第i种重金属的污染指数Pi的计算公式 Pi=Xi （1） Si

，Si为第i种重金属元素的背景平均式中Xi为第i种重金属元素浓度实测值（µg/g）

值(µg/g)。

5.1.2 第j个采样点综合污染指数Cj的计算方法

Cj= （2）

式中Pjmean是第j个采样点所有重金属元素污染指数的平均值，Pjmax是第j个采样点所

有重金属元素污染指数的最大值。

5.1.3 第t类功能区的综合污染指数平均值Ct计算公式 ∑C(j)tn

Ct=j=1

nt （3）

式中Ct(j)为第t类功能区第j个采样点的内梅罗指数，nt为第t类功能区中的样本个

数。 5.1.4 对各区内梅罗指数平均值Ct进行排序。

从排序得到该城区内各类功能区的污染程度情况，再根据各功能区的各重金属元素的含量平均值和背景值进行比较，超出背景值范围则说明该功能区可能已经受到污染，从而确定各个功能区的污染程度。

5.2 模型二的建立

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

该模型采用灰色关联分析法对重金属污染的主要原因进行分析。

5.2.1 指定参考数据列

以内梅罗指数为参考数据列，得到第t类功能区的参考数据列为

x0t=(Ct(1),Ct(2),L,Ct(nt)) （4）

式中Ct(j)为第t类功能区第j个采样点的内梅罗指数（其中j=1,2,···，nt），nt为对第t类功能区的采样点个数。

5.2.2 确定比较数列

以8种重金属元素的浓度实测值为比较数列，如下所示

xit=(Xij(1),Xij(2),L,Xij(nt)) （5）

式中Xit(j)为第i种重金属元素在第t类功能区中的第j个采样点的浓度实测值

（µg/g），nt为第t类功能区的采样点个数。

5.2.3 对参考数列和比较数列进行初值化

为统一量纲并使所有数据列交于同一点，所有数据均用第一个数据除，得到一组新数列 y0t=(1,

Ct(2)C(n),L,ttCt(1)Ct(1)

X(2)X(n)yjt=(1,i,L,itXi(1)Xi(1) （6）

式中Ct(j)为第t类功能区第j个采样点的内梅罗指数（其中j=1,2,···，nt），Xit(j)为第i种重金属元素在第t类功能区中的第j个采样点的浓度实测值（µg/g），nt为第t类功能区中的采样点个数。

5.2.4 计算关联系数

第t类功能区中第i种元素的第k个样本点对参考数列的关联系数为 min( it(min))+ρmax( it(max))it ξit(k)=it （7） |y0t(k) yit(k)|+ρmax( it(max))it

其中

min( it(min))=min(min|y0t(k) yit(k)|) （8） ititk

max( it(max))=max(max|y0t(k) yit(k)|) （9） ititk

式中yit(k)为第t类功能区中第i种元素的第k个样本点元素浓度初值化后的值，ρ为

分辨系数，一般在0与1之间取，此处取ρ=0.5。

5.2.5 计算关联度

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

第t类功能区第i种金属元素的关联度定义为 1 rit=N∑ξ

k=1Nit(k) （10）

式中ξit(k)为第t类功能区中第i种元素的第k个样本点对参考数列的关联系数。

5.2.6 根据关联度大小进行排序

从排序得到各功能区各元素与重金属污染程度关联度的大小，从而得到主要重金属污染种类，对应于其在实际生产和生活中的使用，折射出重金属污染的主要原因所在。

5.3 模型三的建立

5.3.1矢量的通量

ur在矢量场A中，假设

urr dS=ndS （11）

ururururr为一有向曲面S的面元矢量，n为面S的外法线方向，则面元处的矢量A与面元矢量dS

urr点积称为A向n方向穿过dS的通量，记作

urur dΦ=A dS=AcosθdS （12）

ururur其中，θ为dS处的A与dS的夹角。

线

图1 面元矢量dS

及矢量的通量

ur将闭合曲面S的各面元上的dΦ相加，其和表示A向正侧穿过曲面S的通量，记作

urururrΦ= ∫SA dS= ∫SA ndS （13）

5.3.2通量值与正源、负源

ur从场通量与发出场通量的场源之间的关系来看，不论A代表何种量，只要式（13）

中的积分值Φ>0，就表明突出闭合面S的通量线多于穿入S的通量线，这时S内必有

ur发出通量线的源（称为A的正源）；而当Φ<0时，表明穿入S的通量线多于穿出S的通

ur量线，这时S内必有吸收能量的沟（或汇）（沟或汇称为A的负源）；而当Φ=0时，表

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

明穿出S的通量线等于穿入S的通量线，这时S内或者既无正源也无负源，或者两种源

ur的代数和为0，或者说S内没有矢量A的净源。

5.3.3线通量的定义

根据题目所给数据，可以得知采样点所对应的空间坐标，通过对海拔进行大小排序，发现海拔较大处均为山区，而其他功能分区的海拔较低，这也符合城市一般都位于地势起伏较小的平原地区的实际情况，而山区又不可能是重金属元素污染源的所在处，故可对通量模型进行简化，忽略海拔大小所造成的影响。

此处引入线通量的概念，如图2所示，

r图2 线元矢量dl及矢量的通量

一般情况下，穿过面S的通量为

Φ=∫SururA dS （14）

然而，当S的宽度趋近于零时，可以仿照面电流密度矢量的定义方式，规定线通量来描

r述通过线元dl的通量

urr Φl=∫A dl （15） l

5.3.4 通量模型的建立

正如解决问题一时所遇到的，已知采样点数据无法遍布整个平面，平面上有些区域是没有采样点的，而且采样点间的水平面坐标不是等间距的，在模型建立前同样运用“v4”插值法对数据进行处理，得到等间距网格所对应的不同元素浓度值。

如图3，设平面坐标系直角坐标系上的任意一非边缘点(x0,y0)，在其相邻位置上的点为(xi,yi)，其中i=1，2，3，4，分别表示上，下，左，右四个方向。

1,y1)

(x3,y)(x,y)

2

,y2)

图

3 网格点分布

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

记对应点第i种元素的浓度为aij，其中i=1,2,…,8，j=0,1,…,4。

由上文可知，相邻点是等间隔点，且间隔大小与人为设定的插值网格有关，记沿x轴正向的间隔为dx，沿y轴正向的间隔为dy。

点(x0,y0)所在位置的第i种元素浓度相对于第j个点所在位置相应元素浓度对距离的变化率为 kij=

其中 ai0 aijd1,i=1,2,...,8,j=1,2,3,4 （16）

dx,j=3,4 d1= （17） d,j=1,2 y

为了求出以点(x0,y0)作为源的情况下，矢量（在本问题中即浓度的变化）穿过闭合包络线（图中矩形）的线通量，必须知道矢量沿任意方向且与矩形相交所在的点的瞬时变化率，然而，浓度在任意点的变化无法通过已知数据准确反映出来，在所取的包络线度对于整个地区的线度来说足够小的情况下，认为浓度沿某个特定方向是均匀变化的，因此可以用两点间浓度的平均变化率来近似矢量沿两点连线方向且与矩形相交所在的点的瞬时变化率。

再者，如果要准确得到某点作为源时的线通量值，必须考虑矢量沿所有方向穿过包络的线元积分，这也是难以实现的，因为我们无法得到矢量的准确表达式。因此这里作以下理想化处理：矢量线只沿上，下，左，右四个方向辐射，且都为直线，即垂直穿过包络线。

综上所述，可以定义出通量模型中某点(x0,y0)作为重金属元素i的通量源，穿过矩形包络线的线通量为

Φli=∑kijd2,i=1,2,...,8,j=1,2,3,4 （18）

j=14

其中

dx,j=1,2 d2= （19） d,j=3,4 y

根据以上方法，可以求得网格线上各点的线通量值，通过观察线通量的分布和等线通量线的疏密程度，即可确定污染源的大致位置。

5.4 模型四的建立

该模型采用主成分分析法对影响城市地质环境的因素进行分析，主成分分析法的步骤如下所示。

5.4.1 新变量表示

用变量Q1、Q2···Qp表示气温、降水量、蒸发量、水体污染情况、地面坡度、岩

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

性特征、历史地震活动情况、工程地质以及土壤污染等综合指标。

5.4.2计算相关系数矩阵

z11 z21 Z= M zp1z12z22Mzp2Lz1p Lz2p (20) M Lzpp

其中zij（i，j=1，2，…，p）为原变量Qi与Qj的相关系数,且zij=zji(i≠j)，其计算公

式为

∑n

(Qki Qi)(Qkj Qj)

zij= （21）

5.4.3计算特征值与特征向量

a、解特征方程 λI Z=0 （22）

并按特征值的大小排列，使得

λ1≥λ2≥L≥λp≥0 （23）

b、分别求出对应于特征值λi的特征向量

β β11 β12 1p

ββ

21β 22

β2p

1= ,β2= ,...,

M Mβ

p= M

β

p1 βp2 βpp

5.4.4计算主成分贡献率及累计贡献率

第i个主成分的贡献率 fλi

i=2,L,p) （25） ∑p(i=1,

λk

k=1

第i个主成分的累计贡献率

（24）

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

∑λ

Fi=ik∑λ

k=1k=1p(i=1,2,L,p) （26） k

一般取累计贡献率达85%~95%的特征值λ1,λ2,L,λm所对应的第1、第2···第m个主成分，得到主成分表达式

Gi=β1iQ1+β2iQ2+···+βpiQp（i=1,2,...,m） （27）

5.4.5计算主成分载荷

Lij=ij (i=1,2,...,m;j=1,2,...,p) （28）

第i主成分在第j变量上的载荷越大，即Lij越大，说明该主成分对该变量的解释越

充分，用此法验证该模型所选取的主成分能否代表所有影响因素作为主要因素影响城市地质环境。

5.4.6 分析主成分与城市地质环境演变模型的关系

将所得主成分作为影响地质环境的主要因素，通过研究这些因素得到城市地质环境演变模式。对于此模型，由于没有实测数据，且题目只要求建立模型，故不涉及模型求解。

六、模型求解

6.1 对模型一进行求解

6.1.1利用matlab7.1运行附录中的程序p1.m得到8种重金属元素在该城区的空间分布

图

图4 As的浓度空间分布图 图5 Cd的浓度空间分布图

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

图6 Cr的浓度空间分布图 图7 Cu的浓度空间分布图

图8 Hg的浓度空间分布图 图9 Ni浓度的空间分布图

图10 Pb的浓度空间分布图 图11 Zn的浓度空间分布图

（1）As的浓度空间分布图如图4所示。由图可知，As在该城区的分布主要集中在坐标分别为（5000,8000），（13000,2500），（19000,10000），（27000,12000）的附近,其他区域虽然也有所分布,但浓度低。

（2）Cd的浓度空间分布图如图5所示。由图可知，Cd在该城区的分布比较分散，几乎遍布整个城区，且由颜色可知该城区Cd的污染严重。

（3）Cr的浓度空间分布图如图6所示。由图可知，Cr在该城区的分布主要集中在坐标（5000,6000）附近，在x=5000往右的区域Cr的含量很均匀，浓度偏低，在x=5000往左，Cr的含量逐渐减少，浓度以阶梯形式下降。

（4）Cu的浓度空间分布图如图7所示。由图可知，Cu集中分布在坐标（3000,3000）附近，除此之外的其他区域的Cu浓度都很低，基本可以忽略。

（5）Hg的浓度空间分布图如图8所示。由图可知，Hg在该城区的空间分布主要集中在坐标为（3000,3000），（14000,2000），（15000,9000）的点附近，除此之外的其他区域Hg分布均匀，且含量很低。

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

（6）Ni的浓度空间分布图如图9所示。由图可知，Ni在该城区的空间分布主要集中在坐标（4000,6000）附近。而在其他区域的分布相对分散，但浓度较低。

（7）Pb的浓度空间分布图如图10所示。由图可知，Pb主要分布在坐标为（2000,3000），（5000,5000），（4000,11000）附近,除此之外其他地区浓度很低。

（8）Zn的浓度空间分布图如图11所示。由图可知，Zn的分布主要集中在坐标为（10000,4500），（13000,10000）附近，这几处浓度最高,而在（4000,4000）和（13000,3000）附近浓度相对低些，除此之外其他区域的浓度很低。

6.1.2 利用matlab7.1运行附录中的程序p2.m，得到五类功能区的内梅罗平均指数：

[C1,C2,C3,C4,C5]=[4.5586,14.6816,1.6529,11.9054,3.7442] （29）

由式（29）可得：

C2>C4>C1>C5>C3 （30）

C值越大，表明受污染程度越严重，利用excel软件对各功能区各重金属元素浓度平均值进行处理，得到表1数据。 表1 各功能区各金属元素浓度平均值与背景值的比较

µg/g) ng/g) µg/g) µg/g) ng/g) µg/g) µg/g) µg/g) 干道路园绿地景值范由表1可以看到，除了山区的采样点外，其他所有采样点的值都超出了背景值的范围，这说明了除了山区，其他四个功能区都是重金属污染区。再由式（30），可对重金属污染区进行污染程度划分，得到各类功能区土壤中重金属元素的污染程度如表2所示。

表2 各功能区土壤中重金属污染程度

功能区 生活区 工业区 山区 公园绿地区 污染程度 中度污染 极重污染 无污染 重度污染 轻度污染

我们可以通过以下方法对表2的结果进行验证：以采样点所在位置为x、y轴，以每个样本点的内梅罗综合指数为z轴作图，得到图12如下。

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

图12 内梅罗综合指数在空间的分布与各功能区分布对比

由图可知，污染最严重的点附近主要是工业区，针对图中工业区的分布，发现有些地方虽然是工业区，但却不是污染最严重的点，这是因为这些工业区不是大型工业区，污染不是特别严重，因此可认为图中污染程度最严重的点分布着该城区的主要大型工业基地；主干道路区、生活区、公园绿地区、山区污染依次减小，从而验证了表2结论的正确性。

6.2 模型二的求解

利用matlab7.1运行程序附录中的p3.m，以附件1、2所给数据为参考数列和比较数列，求得关联度矩阵

0.9732 0.9363 0.9572 0.9160 0.9764 0.9675 0.9353 0.9030 0.9663 0.9447 0.9630 0.9322 0.9924 0.9629 0.9431 0.9275

r= 0.8693 0.7249 0.7572 0.7668 0.9236 0.7687 0.7062 0.7059 （31） 0.9668 0.9354 0.9619 0.9244 0.9833 0.9664 0.9479 0.9188 0.9437 0.8710 0.9342 0.9013 0.9735 0.9415 0.8895 0.8648

其中rij为第i类功能区第j种元素与该功能区的内梅罗指数的关联度。

由式（31）即可得到各功能区中各重金属元素与内梅罗指数C的关联度如表3所示。 表3 各功能区中各重金属元素与内梅罗指数的关联度

As Cd Cr Cu Hg Ni Pb Zn 生活0.9732 0.9363 0.9764 0.9675 区

工业0.9663 0.9447 0.9924 0.9629 0.954 区

山区 0.8693 0.7249 0.9236 0.7687 道路0.9668 0.9354 0.9833 0.9664 区

公园0.9437 0.8710 0.9735 0.9415 绿地

由表3可知，根据各功能区中各金属元素平均关联度可得污染程度：工业区>道路区>生活区>公园绿地>山区。从而验证了本文在第一个模型中所得的各功能区污染程度

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

的准确性。

各功能区中各重金属元素与内梅罗指数的关联度大小排序如表4所示。

表4 各功能区中各重金属元素与内梅罗指数的关联度大小排序

功能区 关联度排序

生活区 Hg>As>Ni>Cr>Cd>Pb>Cu>Zn

工业区 Hg>As>Cr>Ni>Cd>Pb>Cu>Zn

山区 Hg>As>Ni>Cu>Cr>Cd>Pb>Zn

道路区 Hg>As>Ni>Cr>Pb>Cd>Cu>Zn

公园绿地区

Hg>As>Ni>Cr>Cu>Pb>Cd>Zn

根据表4，作出以下分析：

该城区中Hg元素和As元素在五个功能区中所占关联度都分别为第一，第二位，说明Hg污染和As污染是该城区中的主要重金属污染物。Hg主要用于催化剂、锅炉、汞泵等，而As作为合金添加剂印刷合金、黄铜（冷凝器用）、蓄电池栅板、耐磨合金、高强结构钢及耐蚀钢等，说明该地区冶金工业发达，同时产生的工业污染也严重。并且，工厂没有做好工业污水、工业废弃物等的处理，从而导致土壤中Hg和As含量偏高。此外,Hg和As是农药的主要组成元素，该城区公园绿地区和山区的Hg和As浓度仍偏高，一方面是由于工业区污染的扩散，另一方面也是因为农药的过度使用。

Ni和Cr元素也是重金属污染的主要原因，它们都是不锈钢的主要组成成分。不锈钢是建筑用金属材料中强度最高的材料之一，具有良好的耐腐蚀性，能使结构部件永久地保持工程设计的完整性。含Cr不锈钢还集机械强度和高延伸性于一身，易于部件的加工制造，可满足建筑师和结构设计人员的需要。因此不锈钢应用非常广泛，在生活区、公园区、山区以及主干道路区都有所应用，因此这三个区域会出现Ni和Cr重金属污染物。同时，由于需求量大，工业区的不锈钢生产产业发达，造成Ni和Cr污染物的产生。

综上所述，该地区重金属污染的主要原因是冶金业在该地区造成工业污染以及农业上滥用农药。

6.3 模型三的求解

利用matlab7.1运行附录中的程序p4.m，得到如下图像。

图13 As元素线通量分布及等通量线 图14 Cd元素线通量分布及等通量线

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

图15 Cr元素线通量分布及等通量线 图16 Cu元素线通量分布及等通量线

图17 Hg元素线通量分布及等通量线 图18 Ni元素线通量分布及等通量线

图19 Pb元素线通量分布及等通量线 图20 Zn元素线通量分布及等通量线

观察以上8幅图像，对图中颜色越偏向暖色系的点，其线通量值越大，亦即相应点为污染源的可能性最大。通过程序运行所得到的数据，可以给出各种元素污染源所在水平面坐标如表5所示

表 5 各元素污染源坐标 元素种类

As

Cd

Cr

Cu

Hg

Ni

Pb

Zn

可能的污染源坐标 （27600，11800） （17700，3800） （3300，6000） （2400，3600） （2400，3600） （22200，12200） （4800，5000） （13800，9600） （18000，10000） （5700，8600） （4600，4500） （3300，6000） （15300，9200） （3000，6000） （2400，3800） （5100，7200） （4800，7400） （9300，4400）

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

七、模型优缺点

对于模型一，内梅罗指数法中最大污染指数对综合评价指数影响最大，突出了污染指数最大的污染物对环境的影响。从内梅罗指数法得出的效果来看，突出污染指数最大的污染物的影响让结果更合理，使各功能区的分布更清晰，减少了分析污染程度的难度。然而，由于内梅罗指数的综合评价指标并没有统一标准，因而对评价分级中人为因素的影响是不可忽略的。

对于模型二，灰色关联分析法并不依赖于数据量，也没有要求数据符合经典的分布规律。我们使用这个方法得出的结果与定性分析相吻合，符合一般的规律。然而，由于分辨系数的选取并不固定，得出的结果并不一定是最准确的。国内已经有学者指出[8]本文使用的关联度并不满足规范性[9]。

对于模型三，类比电磁场理论中的通量概念，并且仿照电流面密度矢量的定义方式，给出了线通量的概念，很好地避免了考虑影响不大的海拔因素，在建立计算线通量的表达式时，选择平面图形作为包络线，以及用平均浓度变化代替瞬时变化，都使得模型得到简化，计算简单又不失准确性。然而，在没有得到浓度变化率的准确表达式的情况下运用通量概念确定源位置，结果受采样点数据的影响，在采样点较稀疏的情况下，可能存在较大的误差，为简化计算所作出的理想化处理，也增加了结果的不稳定性。

对于模型四，使用主成分分析能从大量数据中提取出主成分，从而达到降维的目的。使用提取的主成分对地质环境的演变模式的进行分析可以避免数据自身的相关性，降低算法的复杂度。然而主成分分析只提取得分较前的几个元素，并不是完全包含原来数据的所有信息，这可能会对结果产生微小的影响。有时候由于主成份得分比较平均，要使用的主成份比较多，对分析过程的简化并没有想象中理想。

八、参考文献

[1] 薛定宇、陈阳泉，《高等应用数学问题的MATLAB求解》，北京：清华大学出版社，2008.10。

[2]楚天科技，《MATLAB R2008科学计算实例教程》，北京：化学工业出版社，2009.6。

[3]王春光、张思冲、辛蕊、罗娇赢、周晓聪，《哈尔滨市东郊菜地土壤重金属环境质量评价》，中国农学通报，2010.26（2）：262-266，2010。

[4]《土壤重金属污染》，/view/1735172.htm，2011.9.10。

[5]《灰色关联分析方法》，/p-63277922.html，2011.9.10。

[6]马冰然，《电磁场与电磁波》，广州：华南理工大学出版社，2007.8。

[7]周爱国、周建伟、梁合诚、唐朝晖、甘义群、孙自永等，《地质环境评价》，武汉：中国地质大学出版社，2008.5.1。

[8] 张绍良，张国良，《灰色关联度计算方法比较及存在问题分析》，系统工程，第14卷第3期：45-49,1996.5。

[9] 邓聚龙著，《灰色系统理论教程》，武汉：华中理工大学出版社，1990。

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

附录 附录

一、程序p1.m:p1.m:

clear;close all; %清空旧数据，关闭旧图像

data1=xlsread('F:\2011高教杯全国赛\cumcm2011Problems中文版\A\cumcm2011A附件_数据.xls','附件1'); %导入题目所给附件1

data2=xlsread('F:\2011高教杯全国赛\cumcm2011Problems中文版\A\cumcm2011A附件_数据.xls','附件2'); %导入题目所给附件2

x=data1(:,2); %储存水平面坐标的x值

y=data1(:,3); %储存水平面坐标的y值

z=data2(:,2:9); %储存8种不同元素的浓度值

for a=1:8

[x1,y1]=meshgrid(0:300:30000,0:200:20000); %建立网格线

z1=griddata(x,y,z(:,a),x1,y1,'v4'); %以“v4”法进行插值运算

figure;

contourf(x1,y1,z1,15); %作出8种不同元素的填充等浓度线 end

二、程序p2.m:p2.m:

clear;close all; %清空旧数据，关闭旧图像

data1=xlsread('F:\2011高教杯全国赛\cumcm2011Problems中文版\A\cumcm2011A附件_数据.xls','附件1'); %导入题目所给附件1

data2=xlsread('F:\2011高教杯全国赛\cumcm2011Problems中文版\A\cumcm2011A附件_数据.xls','附件2'); %导入题目所给附件2

x=data1(:,2); %储存水平面坐标的x值

y=data1(:,3); %储存水平面坐标的y值

s=[3.6 130/1000 31 13.2 35/1000 12.3 31 69]; %储存附件所给背景值

c=data2(:,2:9); %储存8种不同元素的浓度值

c(:,2)=c(:,2)/1000; %统一量纲，附件中Cd和Hg浓度的量纲与其他不同 c(:,5)=c(:,5)/1000;

for a=1:8

p(:,a)=c(:,a)/s(a); %污染物污染指数

end

for a=1:319

pmean(a,1)=sum(p(a,:))/8; %各采样点污染物污染指数的均值

end

for a=1:319

pmax(a,1)=max(p(a,:)); %各采样点污染物污染指数的最大值

end

C=sqrt((pmean.^2+pmax.^2)/2); %采用内梅罗指数法定义综合污染指数

[x1,y1]=meshgrid(0:300:30000,0:200:20000); %建立网格线

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

z1=griddata(x,y,C(:,1),x1,y1,'v4'); %以“v4”法进行插值运算

contourf(x1,y1,z1,10); %作出该地区综合污染指数等值线

for a=1:5

ans=find(data1(:,5)==a);

Cave(a)=sum(C(ans))/length(ans); %分别以各功能区为单位求污染物污染指数均值 end

hold on;

qu=data1(:,5);

new=[x,y,qu]; %导入水平面坐标与分区编号

b1(:,1)=find(new(:,3)==1);

b2(:,1)=find(new(:,3)==2);

b3(:,1)=find(new(:,3)==3);

b4(:,1)=find(new(:,3)==4);

b5(:,1)=find(new(:,3)==5);

plot(new(b1,1),new(b1,2),'+',new(b2,1),new(b2,2),'o',new(b3,1),new(b3,2),'*',new(b4,1),new(b4,2),'x',new(b5,1),new(b5,2),'d'); %作出该地区不同功能区的标示

三、程序p3.m: p3.m:

clear;close all; %清空旧数据，关闭旧图像

data=xlsread('F:\2011高教杯全国赛\cumcm2011Problems中文版\A\cumcm2011A附件_数据.xls','附件2'); %导入题目所给附件2

data1=xlsread('F:\2011高教杯全国赛\cumcm2011Problems中文版\A\cumcm2011A附件_数据.xls','附件1'); %导入题目所给附件1

data(:,1)=[]; %删掉编号所在列，方便数据处理

data(:,9)=[5.595065696872188;2.050686879701597;3.888298563664709;19.48276529716754;16.93280375665379;22.724400402542187;3.037349192403835;278.4314856449954;326.29517448187767;2.168613213605158;5.463305948625636;3.881272052582544;2.633239515018349;9.46334805067366;5.8640039047304;11.86692132078126;9.495010113989254;5.780916469665128;4.600638095010868;18.078460014979655;2.765714198188635;75.06194642916017;6.317656062533044;3.175982075050526;1.002105755941413;13.027235962605822;4.273270498136975;2.938675903149045;31.277250875574946;17.799553790448112;7.203477918187483;4.402250215125141;5.114637582756907;7.411040169305963;6.071798817309512;30.197995035158257;4.986040063219011;4.933408478760172;6.851727650580298;9.30454412611158;38.970506299568115;13.714554002857652;7.048582202760745;7.128746796901836;12.785511357350094;3.653523595443423;3.509614909556038;2.054432017198534;18.78500460165682;4.56868763930696;3.086422991477065;3.602471019982397;3.286155698475126;15.197347372190814;2.407994127442934;3.918089639224968;3.021220225505681;3.479181017287154;4.062489944841015;2.976220965345462;39.076789250677294;1.56838923950092;1.869842022929421;2.09594617952183;4.393470573406756;3.198853554555947;2.94663449909761;3.82610699445268;2.540367896060084;1.484074279882942;2.053723421933565;1.890947323462811;3.340249275117196;1.646648526651114;2.029100519052495;0.684083342613396;2.041502463345344;0.849988961415488;2.870055943613625;1.738665489003439;1.456952671918555;3.975748198588527;1.45770399240739;6.313275792533697;2.508112770368477;3.02604549399094;1.262863159123092;1.532824738176672;1.339330145544494;5.147561

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

479033372;1.442152129124555;1.728966302675296;35.02906634769579;1.222626258955965;8.991293152408028;3.116240408825262;1.285965859778781;2.392496703609948;8.394421365598285;1.909329903326118;1.700390094907993;1.46384563280339;2.091249430914044;1.989310778249561;1.820909288183112;0.721991608639142;3.556906244033831;1.50022577464769;4.430537788795548;2.021180805875407;1.107893919623452;1.513783484287548;1.499527256025515;2.10496561241004;1.132692260630582;1.164797235418949;1.055152011828799;1.388339303345372;1.164723447855671;2.104980456881918;1.155816659883849;1.433689945085654;1.420131846684555;0.848859610698792;1.507590772124207;0.893707751538136;1.424084754227707;2.809072379953484;2.217762634621828;1.760467032260697;0.868600910162271;0.785063959349192;1.115653623383892;0.973423789788422;4.683151313222338;0.956617527609808;0.680564271380159;2.091250806019084;3.360683902584098;2.823269963244356;4.591082146799722;3.216003158858953;14.932896440395785;3.38560415942609;27.454842601061785;1.54062408546135;2.492161589325188;2.068636869060615;1.559669888632168;5.348974599253076;3.8638929699764;5.183671248264863;3.86501621226541;7.110209527661803;2.012991019014493;1.61390149306727;11.396202111091059;1.39574520061232;2.94716250380622;2.12424866001441;6.801302747207917;4.764490882492501;7.556909863238746;3.106179652475045;4.604014655625449;

4.728045091860476;4.023840561468352;2.058523872346038;1.754391896416406;3.515992402566642;

1.8218092344078;1.181573234475986;4.032281433834817;6.072375228454682;6.111447343916569;2.642092177916958;2.859940944277829;17.26217648689107;5.064555238141181;3.046697499313291;3.669488093813069;303.59174924873116;2.550513604518666;1.782740505385051;1.770041262489848;

1.288107261720891;1.695395431355489;1.629784266635865;2.067574667971848;0.820072898168531;

4.34616170735235;4.456670678972494;2.048816916901536;1.317366468412547;1.979394676027896;1.964451108672633;1.854880051663096;1.638815068443622;1.519312170650063;2.721398466886224;1.104282138532609;1.285057712263679;1.742897069263565;1.47009434648751;1.597802319414337;1.670625620298231;1.097951746944716;2.525967930078434;1.69654927927184;0.857804999413691;1.302557912134571;1.268411868870746;0.961813721752549;1.620671172876886;2.024717153149426;3.46676234173421;2.544597956610798;1.897019235643718;2.935557000653834;1.673044970226424;7.038149538008591;4.964305528845143;5.881794510875182;7.489876278631503;11.762503936638574;5.581108391762609;3.463105883829255;4.023752563509703;1.045236086917905;2.714819672468936;5.788831243687037;36.89834816331109;2.457570213702902;16.557290350656668;3.986540321283792;4.260954110723686;1.97841652814125;2.539937843947859;4.149181176809315;8.541865093086031;2.727017730409122;3.550553333966364;2.65193029696082;7.419405357080878;2.252760548422841;4.361566021185795;2.908941131889256;1.547746029192779;2.200447613584008;2.097486629617464;1.674806247354249;5.80336820381127;8.235493594803586;1.365243527406383;5.466424278504436;1.746001650482619;281.1678260881151;2.027878277925154;1.603350507535081;1.425666817630259;4.410731736715362;1.629830593062367;4.944656546252916;2.397661503951931;1.8365046544483;4.558475475747136;2.083367003332384;1.505881303134054;3.805845787651712;2.876377612647909;2.694348598411806;3.720495304993889;9.019765426645925;4.127722734124237;8.539960986890623;1.826884025593284;1.853164031469232;2.668482120907493;2.489347930271103;2.924402194565862;1.793287839473828;1.298772726866198;1.945842155925883;1.741426780795412;1.969542680131738;2.481690124129614;1.687880546759312;1.243131257536124;1.736710369077962;1.153417193937324;1.795009845550388;1.637614702055958;2.033344927867501;2.073383273331175;1.187243663202835;1.465029466812806;1.527384345918749;1.111523107895161;1.747274548588028;1.129654592607302;1.241113417876147;1.020388252171966;1.115302011771711;1.694002865827037;1.604877964418506;1.886079647564806;1.831851534209567;2.042391605964359;3.203316339598711;6.829487504058084;1.980592

2011年全国大学生数学建模竞赛A题——城市表层土壤重金属污染分析：黄俊彬,李恪睿,陈泽君

93544035;1.463541689554989;5.161772125881541;1.491086375193536;2.080879390395252;1.490208609331122;1.333045146213944;1.740511663057965;2.25579729384667;];%由程序p2得到的各样本点综合污染指数

for a=1:9

data(:,a)=data(:,a)/data(1,a); %消除量纲并使数据交于同一点，用于灰色关联度分析 end

ans1=find(data1(:,5)==1);

newdata1=data(ans1,:); %找出1类区所属数据建立新矩阵

x01=newdata1(:,9); %以综合污染指数作为灰色关联度分析的初值序列 h1=size(newdata1,1);

for a=1:h1

d1(a,:)=abs(x01(a)-newdata1(a,1:8)); %求差序列

end

mind1=min(min(d1));

maxd1=max(max(d1));

rou=0.5; %分辨系数

for a=1:size(d1,1)

for b=1:size(d1,2)

ksi1(a,b)=(mind1+rou*maxd1)/(d1(a,b)+rou*maxd1); %计算关联系数序列

end

end

for a=1:size(ksi1,2)

r(1,a)=1/size(ksi1,1)*sum(ksi1(:,a)); %求关联度

end

ans2=find(data1(:,5)==2);

newdata2=data(ans2,:);

x02=newdata2(:,9);

h2=size(newdata2,1);

for a=1:h2

d2(a,:)=abs(x02(a)-newdata2(a,1:8));

end

mind2=min(min(d2));

maxd2=max(max(d2));

rou=0.5;

for a=1:size(d2,1)

for b=1:size(d2,2)

ksi2(a,b)=(mind2+rou*maxd2)/(d2(a,b)+rou*maxd2);

end

end

for a=1:size(ksi2,2)

r(2,a)=1/size(ksi2,1)*sum(ksi2(:,a));

end

ans3=find(data1(:,5)==3);

newdata3=data(ans3,:);

本文来源：https://www.bwwdw.com/article/0n3j.html