城市表层土壤重金属污染分析
城市表层土壤重金属污染分析
摘要
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。生活废水、汽车尾气、工业排放等都会引起城市表层土壤的重金属污染, 对人体造成很大的危害。因此城市表层土壤重金属研究是非常有意义的。本文根据所给数据，运用了插值拟合结合内梅罗综合污染指数评价法，建立BP 神经网络模型，利用遗传算法，研究人类活动影响下城市地质环境的演变模式。
对于问题一，先用MA TLAB 软件对所给数据进行处理，插值拟合得出8种主要重金属元素在该城区的空间分布图，再用内梅罗指数法建立模型进行求解。首先用EXCEL 对数据进行分析，得出各区的8种重金属的平均浓度。然后结合MA TLAB 软件求出各区的单项污染指数和综合污染指数，进而得出各区的综合污染等级。生活区、工业区、主干道路区为重污染区，公园绿地区为中度污染，山区为轻度污染。
对于问题二，用内梅罗指数法建立模型进行求解，首先用EXCEL 对数据进行分析，得出每种金属在各区的单项污染指数，用EXCEL 处理，进而分析得出重金属污染主要来源于工业废水、生活垃圾、农药等。Cu 和Hg 对各区的污染都很大，As 和Cr 污染主要来源于生活垃圾，Pb 污染主要来源于汽车尾气的排放。
对于问题三，为了确定污染源的位置，我们建立了BP 神经网络模型，该模型能准确反映重金属污染物的传播特征并得出重金属浓度的空间分布，但难以准确地找到重金属污染物的传染源。由于传染源位于最大浓度坐标点，我们利用遗传算法可以搜索到污染源的位置，并绘制了八种重金属的污染源汇总图。
对于问题四，神经网络模型的优点是能学习和存贮大量的输入-输出模式映射关系，而无需事前揭示描述这种映射关系，并具有较强的自学习能力，充分利用了海拔高度的信息；缺点是训练要求样本点容量较大，网络的收敛速度慢，需要较长的训练时间。可以通过搜集前几年该城区八种重金属浓度的采样数据，岩石、土壤、大气、水和生物等因素的相关数据，进一步拓展神经网络模型，得到该城市地质环境的演变模式。
关键字： 插值拟合 内梅罗指数法 神经网络 遗传算法
一、 问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、„„、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
我们通过数学建模解决了以下问题： 问题一：给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
问题二：通过数据分析，说明重金属污染的主要原因。
问题三：分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
问题四：分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，收集信息并建立模型解决问题。
二、 问题分析
问题一属于空间分布和综合评价问题。重金属的传播过程是一个扩散的过程，通常物质扩散模型中物质从高浓度向低浓度扩散且其浓度的分布是连续的，据此我们可以用附表中所给的采样点污染数据为基础借助MA TLAB 软件[1]进行插值拟合得出8种主要重金属污染物在整个城区的空间分布图。对于该城区内不同区域重金属的污染程度的研究可以借助我国《土壤监测技术规范》(HJ/T 166-2004)[2]中推荐的内梅罗指数法进行评价，求出不同区域重金属的污染等级。
针对问题二，根据内梅罗指数评价法，得出每种金属在各区的单项污染指数，进而分析各区受到重金属污染的原因。
针对问题三，将重金属污染传播过程看作网络输入与输出的一类非线性映射, 但是仅通过网络的输入输出数据难以准确寻找重金属污染物浓度的极值（即污染源），而遗传算法具有全局的非线性寻优能力。所以建立神经网络模型，用遗传算法确定污染源的位置来寻找污染源。
针对问题四，首先对神经网络模型进行优缺点分析，然后根据影响城市演化模型的因素，分析还应搜集的数据以及模型如何建立的问题。
三、模型假设
1. 假设题目所给数据真实有效。
2. 假设只考虑题目中所给的8中重金属，不考虑其它重金属。 3. 假设重金属浓度随空间是连续变化的。
4. 假设重金属污染物的传播只受到污染物在液体中的扩散和其他物质对污染物的吸附作用的影响。
四、符号说明
符号设定
P ij C ij S j P N P j,ave P j,max
符号说明
区域i 中第j 个重金属的污染分指数 区域i 中第j 个重金属的实测浓度
第j 元素的评价标准
综合污染指数 平均单项污染指数 最大单项污染指数
五、模型建立及求解
5.1 问题1的求解
5.1.1 重金属的空间分布图
本节主要表述了该城市的地形以及功能区，给出8种主要重金属元素在该城区的空间分布，运用MA TLAB 软件对所给数据进行处理，插值拟合得出该城区的地形图以及8种主要重金属元素在该城区的空间分布图，再用MA TLAB 软件对所给数据进行分析得出功能区散点图：
4
x 10
4
图1：该城区的地形图 图2：功能区散点图
x 10
4
图3：As 浓度(μg/g)的空间分布
x 10
4
图5：Cr 浓度(μg/g)的空间分布
x 10
4
图7：Hg 浓度 (ng/g)的空间分布
x 104
图4：Cd 浓度 (ng/g)的空间分布
x 10
4
图6：Cu 浓度(μg/g)的空间分布
x 10
4
图8：Ni 浓度(μg/g)的空间分布
x 10
4
x 10
4
图9：Pb 浓度(μg/g)的空间分布 图10：Zn 浓度(μg/g)的空间分布
说明：
图1的Z 轴为海拔高度，X 、Y 轴为地理坐标值（单位：m ）。 图2-图10 的X 、Y 轴为地理坐标值（单位：m ）。
5.1.2 模型建立
对于该城区内不同区域重金属的污染程度的研究可以借助我国《土壤监测技术规范》(HJ/T 166-2004)[2]中推荐的内梅罗指数法进行评价，求出不同区域重金属的污染等级。 内梅罗指数法是当前国内外进行综合污染指数计算的最常用的方法之一。该方法先求出各因子的分指数(超标倍数) ，然后求出个分指数的平均值，取最大分指数和平均值计算。 (1) 单项污染指数
通过单因子评价，可以确定主要的重金属污染物及其危害程度。一般以污染指数来表示，以重金属含量实测值和评价标准相比除去量纲来计算污染指数： P ij
C ij
, （1）
S j
注：（1）式中P ij 为区域i 中第j 个重金属的污染分指数；C ij 为区域i 中第j 个重金属的实测浓度；S j 为第j 元素的评价标准。当Pij ≤1时，表示土壤未受该因子污染，当Pij >1时，表示土壤受该因子污染。 （2）综合污染指数
单项污染指数只能反映各个重金属元素的污染程度，不能全面地反映土壤的污染状况，而综合污染指数兼顾了单项污染指数平均值和最高值，可以突出污染较重的重金属污染物的作用。综合污染指数计算方法如下： P N
（2）
注：（2）式中P N 是采样点的综合污染指数；P j,ave 为平均单项污染指数；P j,max 为最大单项
污染指数。内梅罗综合污染指数的分级标准见表5-1
表5-1 内梅罗综合污染指数的分级标准
5.1.3 模型求解
本文以背景值作为评价标准进行求解，用EXCEL 对文中所给数据进行分类，把数据分入1类区、2类区、3类区、4类区、5类区。然后得出各区重金属含量的平均值，如下表：
表5-2 各区重金属含量平均值
然后根据公式（1）、（2）结合MA TLAB 软件算得各区重金属单项污染指数和综合污染指数，如下表：
表5-3 各区重金属单项污染指数和综合污染指数
再由内梅罗综合污染指数的分级标准得出各区的综合污染等级，如下表：
表5-4 各区的综合污染等级 于中等污染区，山区属于轻度污染区。
5.2 问题2的求解
根据问题一得出的各区重金属单项污染指数和综合污染指数，用EXCEL 制成表格和图标进行直观分析。
5.2.1 各区重金属污染原因分析
由表5-3中的单向污染指数绘制各区重金属单项污染指数柱状图。分析表5-4和图11，根据实际情况分析重金属污染的主要原因。
图11：各区重金属单项污染指数表
由表5-4可以看出，5类区的污染程度排序为：工业区>主干道路区>生活区>公园绿地区>山区
工业区污染最大，其中Cu 和Hg 对工业区的污染最大。分析其污染的主要来源有金属矿山的开采、冶炼、重金属尾矿、冶炼废渣和矿渣堆放等。
主干道路区的污染也很大，其污染主要体现为大气污染，来自工业生产排放的废气、汽车尾气排放生的大量含重金属的有害气体和粉尘。
生活区的污染主要来源有生活垃圾、污水、医疗垃圾和细菌等。
公园绿地区重金属污染程度较轻，其主要来源有农业农药和化肥。但是，植物有吸收金属矿物的作用，相对减轻了重金属的危害。
山区重金属污染最轻，表明不仅重金属的来源相对较少，而且，重金属污染与重力和海拔也有一定关系，山区中植物和一些微生物都有净化的作用。
5.2.2 重金属污染来源分析
由表5-3中的单向污染指数绘制重金属在各区的单项污染指数柱状图，并分析重金属污染的主要原因。
图12：重金属在各区的单项污染指数
从图中可以看出，Cu 和Hg 对各区的污染最大。其他重金属也对各区有影响。其主要来源如下表所示：
表5-5 重金属污染主要来源
5.3 问题3的求解
由于定性分析重金属污染物的传播特征需要较强的专业知识背景，而所给数据的变量不
多（变量仅坐标、高度与功能区），统计学方法又不足以完全揭示变量与污染物传播特征的关系。一方面，人工神经网络具有较强的自组织、自适应与自学习能力，能够在未完全了解重金属污染物传播机理的情况下，完成自变量、变量间与重金属污染物浓度之间的非线性映射；另一方面，将重金属污染传播过程看作网络输入与输出的一类非线性映射, 但是仅通过网络的输入输出数据难以准确寻找重金属污染物浓度的极值（即污染源），而遗传算法具有全局的非线性寻优能力。综合上述考虑，我们建立神经网络结合遗传算法求解重金属污染源的数学模型。
5.3.1 建立BP 神经网络模型 （一）BP 神经网络原理
BP 神经网络是一种多层前馈神经网络, 该网络的主要特点是信号前向传递, 误差反向传播。在前向传递中, 输入信号从输入层经隐含层逐层处理, 直至输出层。每一层的神经元状态只影响下一层神经元状态。如果输出层得不到期望输出, 则转入反向传播, 根据预测误差调整网络权值和阈值, 从而使BP 神经网络预测输出不断逼近期望输出[3]。 （二）建立BP 神经网络步骤
首先运用MA TLAB 中函数“mapminmax ”实现数据归一化处理，考虑到地形对重金属污染物传播的影响，构建BP 神经网络的隐藏层，经过多次测试确定隐藏层节点数N=100，学习速度为0.01，训练次数为300。matlab 神经网络工具箱中函数“traingdm ”实现了BP 神经网络训练学习。为了得到真实的数据值，运用MA TLAB 中函数“mapminmax ”对输出数据进行反归一化。最后分析预测输出与期望输出之间的误差。 （三）结果分析
以As 为例，训练结束后的神经网络性能图为：
10
M e a n S q u a r e d E r r o r (m s e )
10
10
10
10
131 Epochs
图13：As 的神经网络性能图
如图13，训练在第125次迭代过程达到均方误差最小，MSE=0.0015604时训练结束。 用训练好的神经网络预测重金属As 污染，见图14-图15。
神经网络预测误差
输出
样本
图14：As 的预测与期望 图15：As 的预测误差
由图14-图15知，预测输出和期望输出相差不大，虽然个别点上的误差较大FF0C 由于污染物的浓度在空间上的分布是渐变的，不会影响污染源的寻找，所以基本符合预处理的要求。
用训练好的神经网络模拟出了该城区8种重金属的污染传播特征图，如下：
图16：该城区8种重金属的污染传播特征图
5.3.2 建立遗传算法求解模型 (一) 遗传算法原理
BP 神经网络模型模型能准确反映重金属污染物的传播特征，并得出重金属浓度的空 间分布，但是难以准确地找到重金属污染物的传染源。而遗传算法具有全局的非线性寻优能力。综合上述考虑，我们建立遗传算法求解重金属污染源的数学模型。算法流程图如下：
图17：遗传算法流程
（二）结果分析
以As 为例，As 的污染源的性能图像如下：
遗传算法查找As 污染性能图
浓度
进化代数
图18：As 污染性能图
As 的最大浓度为13.2624μg/g；坐标x=1857.3;坐标y=3208.2。所以我们可以确定As 的污染源在（1857.3，3208.2）附近，其浓度为13.2624μg/g。由此能说明该模型能较为准确地反映污染源位置。
根据遗传算法得出8种重金属污染源坐标以及浓度，见下表：
表5-6 8种重金属污染源
观察样本点中最大值点的位置和浓度，比较发现8种重金属通过模型得到的污染源的位置基本处于样本点中最大值点的附近，反映该模型对污染源位置的寻找具有较好的准确性；同时，计算得到的Cu 和Pb 在污染源位置的浓度高于样本点浓度的最大值，效果较佳；但是As 在污染源位置的浓度低于样本点浓度的最大值，效果略差。
5.4 问题4的求解
5.4.1 模型的优缺点
在解决问题一和问题二时，我们运用了内梅罗指数法，它的优点是数学过程简便，物理概念清晰，评价方式简单便于决策，对数据的处理考虑到了各个散点数据间的联系，加入权重进行综合内梅罗指数排名；缺点是内梅罗综合指数过分突出污染指数最大的重金属污染物对环境质量的影响和和作用，在评价时可能会人为地夸大或缩小一些因子的影响作用，使其对环境质量评价的灵敏性不够高，在某些情况下，它的计算结果难以区分土壤环境污染程度的差别。
由于浓度分布是多种因素共同影响的结果，很难准确地给出各种因素共同作用的公式。在解决问题三时，我们引用了BP 神经网络模型，它很好地解决了这个问题。神经网络具有较强的自组织、自适应与自学能力，能够在未完全了解重金属污染物传播机理的情况下，完成自变量、因变量之间的非线性映射。但是由于神经网络构建所需数据量较大，而题目所给原始数据只有300 多个，此数据量较小，所以不能更加精确地拟合出重金属的空间分布情况。
5.4.2 可能需要的信息
六、模型推广
本文考虑的只是重金属对土壤的污染问题，我们可以把它推广到重金属对植物和动物的影响，从而有利于对农作物的培育和动物的养殖，甚至可以确定对人体带来的危害，也可以应用到其它金属元素对土壤的污染和影响，从而研制促进农作物生长的化肥，有利于农业的发展。
参考文献
[1] 张志涌，《精通MATLAB 6.5版》[M].北京：北京航天航空大学出版社，2003，234-302。 [2] HJ/T166-2004，《土壤环境监测技术规范》[S]. 北京：中国标准出版社，2004。 [3] matlab 中文论坛编注，matlab 神经网络30 个案例分析，北京：北京航天航空大学出版社，2010.
附录
附录一：城区地形分布图的MATLAB 程序：
clc,clear load data
Z=data(:,3); %竖直坐标 X=data(:,1); %横坐标 Y=data(:,2); %纵坐标
[x y]=meshgrid(0:100:3e4,0:100:2e4); %根据坐标范围划分，作为插值使用 z=griddata(X,Y,Z,x,y,'v4');%该命令实现插值，v4是插值方法 mesh(x,y,z) %绘制空间图
附录二：功能区分布散点图的MATLAB 程序：
clear load data load data1 figure
[c1 d1]=find(data(:,4)==1); %即找到功能区1, 返回的c1是每个数据对应的位置编号，即行标号；d1就是1
x1=data(c1,1); %功能区1的横坐标 y1=data(c1,2); %功能区1的纵坐标 plot(x1,y1,'rsquare');
[c2 d2]=find(data(:,4)==2); %同上 x2=data(c2,1); y2=data(c2,2); hold on
plot(x2,y2,'ksquare')
[c3 d3]=find(data(:,4)==3); x3=data(c3,1); y3=data(c3,2);
plot(x3,y3,'msquare')
[c4 d4]=find(data(:,4)==4); x4=data(c4,1); y4=data(c4,2);
plot(x4,y4,'ysquare')
[c5 d5]=find(data(:,4)==5); x5=data(c5,1); y5=data(c5,2);
plot(x5,y5,'gsquare')
附录三：重金属在该城区空间分布图的MATLAB 程序：
load data load data1 As=data1(:,1); Cd=data1(:,2); Cr=data1(:,3); Cu=data1(:,4); Hg=data1(:,5); Ni=data1(:,6); Pi=data1(:,7); Zn=data1(:,8); X=data(:,1); Y=data(:,2);
[x y]=meshgrid(0:100:3e4,0:100:2e4); z1=griddata(X,Y,As,x,y,'v4'); z2=griddata(X,Y,Cd,x,y,'v4'); z3=griddata(X,Y,Cr,x,y,'v4'); z4=griddata(X,Y,Cu,x,y,'v4'); z5=griddata(X,Y,Hg,x,y,'v4'); z6=griddata(X,Y,Ni,x,y,'v4'); z7=griddata(X,Y,Pi,x,y,'v4'); z8=griddata(X,Y,Zn,x,y,'v4'); mesh(x,y,z1) mesh(x,y,z2) mesh(x,y,z3) mesh(x,y,z4) mesh(x,y,z5) mesh(x,y,z6) mesh(x,y,z7) mesh(x,y,z8) contour(x,y,z1) contour(x,y,z2) contour(x,y,z3) contour(x,y,z4) contour(x,y,z5)
contour(x,y,z6) contour(x,y,z7) contour(x,y,z8)
附录四：单项污染指数和综合污染指数求解的MATLAB 程序：
a=[6.27 289.96 69.02 49.4 93.04 18.34 69.11 237.01;7.25 393.11 53.41 127.54 642.36 19.81 93.04 277.93;4.04 152.32 38.96 17.32 40.96 15.45 36.56 73.29;5.71 360.01 58.05 62.21 446.82 17.62 63.53 242.85;6.26 280.54 43.64 30.19 114.99 15.29 60.71 154.24]; b=[3.6,130,31,13.2,35,12.3,31,69]; b=b'; for i=1:5 for j=1:8;
c(i,j)=a(i,j)/b(j,1); end end >> c c =
1.7417 2.2305 2.2265 3.7424 2.6583 1.4911 2.0139 3.0239 1.7229 9.6621 18.3531 1.6106 1.1222 1.1717 1.2568 1.3121 1.1703 1.2561 1.5861 2.7693 1.8726 4.7129 12.7663 1.4325 1.7389 2.1580 1.4077 2.2871 3.2854 1.2431 c=c';
>> mean(c)
ans =
2.4693 5.4270 1.1913 3.8386 2.0393 Pj,max 3.7424 18.3531 1.3121 z=[];
>> x=[2.4693 5.4270 1.1913 3.8386 2.0393]; y=[3.7424 18.3531 1.3121 12.7663 3.2854]; for i=1:5;
z(i)=sqrt((x(i)^2+y(i)^2)/2); end >> z z =
3.1704 13.5331 1.2532 9.4264 2.7343
2.2294 3.0013 1.1794 2.0494 1.9584 12.7663 3.4349 4.0280 1.0622 3.5196
2.2354 3.2854