城市表层土壤重金属污染分析
摘要
问题一中本文利用Surfer 软件，根据提供的数据，绘制出8种重金属的空间立体含量图。再运用地质累积指数法，用I PL 以及I PLzone 表示出单一金属污染负荷指数和8种重金属元素的区域污染负荷指数。由3维立体图可看出As 砷元素的污染程度在1,2,5 区的含量较高，污染较为严重；Cd 镉元素在1,2,4,5区污染较为严重；Cr 铬元素在1,4区污染严重；Cu 铜元素在五个区分布普遍较高，尤其以2,4区污染最为严重；Hg 元素在2区污染较为严重。根据I PLzone 得知，工业区属强污染，其他区域均属于中等污染。
问题二将8种重金属含量进行了主成分分析。通过对比不同金属对环境污染贡献值的差异，分析出污染的主要可能原因，生活区：As ，Cu ，Ni 为主要污染物质，原因可能与居民生活垃圾或该区交通污染有关；工业区：Cr ，Cu ，Hg ，Pb ，Zn 为主要污染物质，导致原因可能是各类工厂废气、废水、废渣的排放。山区：Cr ，Ni 和Zn 为主要污染物质，原因可能是冶金，采矿等活动造成的。交通区：Cr ，Cu 和Ni 为主要污染成分，原因可能是汽车尾气，轮胎老化，车体磨损。公园绿地区：Cr ，Pb 和Zn 为主要污染物质，原因可能是由于汽车尾气排放和工业活动导致公园土壤重金属含量的增加。
问题三中通过查阅资料，分析出不同重金属在地表迁移传播的特点，并利用SPSS 计算出金属间的相关性系数，得出Pb 与Cd 和Hg,Cu 与Cr 具有显著相关性，判断它们可能出自相似污染源或同一污染源。对8种金属浓度进行曲面拟合，将拟合方程结合等值线图，建立了双向逼近模型，确定出主要污染源的位置是：（1）圆心（1862.46, 2827.52），半径1487.60；（2）圆心（3840.71, 5253.45），半径1190.08，单位为Km 。
问题四中将地表重金属含量作为预测地质环境的指标，通过对第三问条件的延伸，在收集了地表坡度、植被覆盖率和年降水量三种数据的基础上，建立F(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y+G(α, β, χ) ；其中G(α, β, χ) 为三个影响因素集成的函数。
关键词：重金属污染 地质累积指数 主成分分析 曲面拟合 双向逼近模型
1问题重述
随着人口急剧增长、工农业和科学技术水平的不断提高，人类正以前所未有的规模和强度开发自然资源和破坏生态环境。其中重金属污染是引起的主要问题之一。
现对某城市城区土壤地质环境进行调查。将该城市划分为五个区，分别是生活区、工业区、山区、主干道路区及公园绿地区，并标号为1、2、3、4、5。将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
本题目要求：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
2基本假设
1 不发生大规模泥石流，山体滑坡等地质灾害，不发生旱涝灾害等自然不可逆灾害现象。
2 假设该城区土层表面土质均一稳定。 3 取样期间，植被覆盖率没有明显变化；
4 假设取样区域内没有河流径流，不会对污染物造成不可预测的扩散。 5 假设所取样区域不发生大规模降水，且气候变化平缓，温度变化不大。
3主要变量符号说明
4问题的分析
此题目中，第一问要求求解8种金属在该区域中的空间分布，需要对数据种的横从坐标与8种不同的重金属进行三维效果处理，从而能够可以从空间上给出其分布情况。
而第二问，需要求解重金属污染的主要原因，这就要求通过对现有的数据进行处理，本文利用了主成分分析方法，对污染原因进行了说明。
第三问，提出根据不同金属在土壤中横向的扩散传播情况进行模型建立，我们通过相关系数的模型建立，分析处理，得出污染源的坐标位置。
第四问，在第三问的基础上，综合考虑三个方面的因素（地表坡度、植被覆盖率、年降水量）进行分析，增加G(α, β, χ) 这一函数来表征对重金属含量的影响。
5问题一的模型建立与求解
5.1 模型的建立
5.1.1各金属元素的空间分布
在利用Surfer 软件的条件下，进行x,y, 以及污染物浓度的绘制，形象直观的得到关于横纵坐标与8种污染物浓度的三维图像：
As
Cd
Cr
Cu
Hg
Ni
Pb
图1
Zn
5.1.2分析污染程度
通过建立污染负荷指数模型，对8种重金属进行对不同地区的污染程度大小的计算求解。污染指数负荷法是Tomlinson 等在从事重金属研究水平的分级研究中提出来的一种评价的方法，该方法被广泛的应用于土壤和河流沉积物重金属污染的评价[1 2]。某一点的污染负荷指数公式如下：
F i =C i /C o i I P L =
式中，Fi 为元素i 的污染系数；Ci 为元素i 的实测质量分数（mg/Kg）；C oi 为元素i 的评价标准，即背景值，一般选用全球页岩平均值作为重金属的评价标准（mg/Kg）；I PL 为某一点污染负荷指数；n 为评价元素的个数。
某一区域的污染负荷指数(I PLzone ) 为：
I PLzone =
式子中，I PLzone 为区域污染负荷指数，n 为评价点的个数（即采样点的个数）。
5.2模型求解
根据地积累指数法的公式，分别算出各区域各元素的I PL ，并统计结果如下：
根据上表并不能完全了解清楚每一个区域的不同程度的污染物的分布情况，通过计算I PLzone ，得出各区域的污染负荷指数。以下是五区的污染程度分级差异情况：
5.3模型补充
通过以上计算分析，得出的是用8种重金属物质的浓度来衡量不同功能区的污染程度，为了补充前一模型的不足，下面用地址累积指数法来进行单一重金属对于不同功能区的衡量：
地质累积指数不仅反映了重金属分布的自然变化特征，而且可以判别认为活动对人为环境的影响，是区分人为活动影响的重要参数[3-4]。其模型如下：
I geo
⎡C n ⎤=log 2⎢⎥
⎣1.5BE n ⎦
其中，C n 表示样品中元素n 的浓度，BE n 表示基线浓度，1.5为修正参数，通常用来表征沉积特征、岩石、地质及其他影响。
一般而言，地质累积指数有2种分级，Forstner 等将其分为7级[5]，而Anon 则分为5级[6]，不同的级别分别代表不同的重金属污染程度，本文采用Anon ，见表3。
表3 地质积累指数分级情况
Igeo 值（Anon ）
5
级别 1 2 3 4 5
污染程度 无污染 无污染到轻度污染
中度污染 中度污染到强度污染
极强污染
6问题二的数据分析与求解
6.1该市土壤重金属主成分分析
为了分析重金属污染的主要原因，运用SPSS 软件对上述不同功能区重金属的各项指标进行主成分分析，数据处理及分析结果见表4：
以上数据结果说明：利用SPSS 软件得出的第一组数据，即生活区的主成分提取出4个（λc >0.8），贡献率分别为45.199%，14.165%，13.432%，10.083%；对表格进行分析，可以得到第一主成分Cd 和Cu 、Pb 有较高的载荷，反映出了住宅区内居民产生的生活垃圾或与该区交通污染有关。第二主成分中Zn 具有较高的载荷，反映出了生活区可能存在如汽车轮胎老化等现象。其他两种主成分较为次要，这里不加以分析。
工业区提取出了2个主成分进行分析（λc >0.8），贡献率分别是65.670%，15.786%，对表格进行分析，可以得到第一主成分中Cr 和Cu 及Hg 、Pb 、Zn 有一定程度的高载荷，反映出这一工业区污染主要原因有电镀、采矿、冶炼、燃料、电池和化学工业等。第二主成分中As 有较高载荷，说明污染来自于各种化学制剂，包括农药，杀虫剂在内的物质。
山区提取出了3个主成分进行分析（λc >0.8），贡献率分别是38.022%，25.445%，19.362%，对表格进行分析，可以发现Cr 和Ni 及Zn 有较高载荷，说明污染原因有汽车轮胎磨损以及煤燃烧产生的粉尘、烟尘中均含有有毒重金属。第二主成分As 的载荷较高；第三主成分中Cu 和Hg 的载荷较高，说明污染与金属加工，日常用品制造等物质有密切关系。
交通区提取出了3个主成分进行分析（λc >0.8），贡献率分别是46.883%，16.084%，12.424%,对表格分析，得知第一主成分中Cr,Cu,Ni 具有较高载荷；第二主成分中Hg 具有较高载荷，说明其原因可能有金属冶炼和燃煤；第三主成分中As 具有很高载荷，达到0.939。
绿地区提取出了3个主成分进行分析（λc >0.8），贡献率分别是48.840%，20.182%，13.202%，对表格进行分析，可以得出第一主成分中Cr ，Pb,Zn 有较高的载荷，说明污染原因可能是工业产生的废物，废气，废渣；第三主成分中，Hg 具有较高的载荷，达到了0.889。
7问题三的模型建立与求解
在这一问中，通过查阅资料[7]了解到8种金属元素在降雨量不发生剧烈改变，气候变化平缓的情况下，大多数不会产生较为明显的横向迁移传播，本文用
SPSS 软件，对所给数据进行了一定量的处理，得到8种金属之间相关性系数，从中得出了某几种金属之间的相关性结论，相关性系数表格见附录，这里不予以赘述。
7.1. 各金属的相关性及浓度的曲线拟合
生活区：As 与Cu; As与Ni; Cr与Ni; Cd与Pb; 具有较强的相关性，说明其可能来自于同一污染源。
工业区： Cu 与Cr ；Cu 与Hg ；Pb 与Zn; As与Ni ；Cr 与Hg; Cr 与Ni; Cr与Pb; Cu与Zn 之间存在较强的相关性，说明这几组重金属物质之间存在一定相似的污染源。
山区：Cr 与Ni ；Zn 与Cd ；Ni 与Pb ；具有较强的相关性，说明其污染源可能源于同一个范围。
交通区：Cr 与Cu; Cr 与Ni; Cu 与Ni; Cd 与Zn; 具有较强的相关性，说明其污染源相似。
公园绿地区：Cr 与Ni; Cu与Zn; 具有较强的相关性，说明其污染源相似。 本文利用Matlab 软件进行重金属浓度曲面拟合，拟合得出的曲面如下图所示(这里只表示As 元素的曲面拟合，其他元素见附录2) ：
As
拟合出的函数的方程式如下： As:
f(x,y) = p00 + p10*x + p01*y
p00 = 7.777 (7.663, 7.891)
p10 = -0.0001683 (-0.0001748, -0.0001617) p01 = 2.846e-005 (1.905e-005, 3.786e-005) Cd:
f(x,y) = p00 + p10*x + p01*y
p00 = 488.8 (480.4, 497.2)
p10 = -0.01223 (-0.01271, -0.01175) p01 = -0.00194 (-0.00263, -0.00125) Cr:
f(x,y) = p00 + p10*x + p01*y
p00 = 85.53 (82.74, 88.31)
p10 = -0.001988 (-0.002148, -0.001828) p01 = -0.0005593 (-0.0007881, -0.0003304) Cu:
f(x,y) = p00 + p10*x + p01*y
p00 = 113.6 (109, 118.1)
p10 = -0.00405 (-0.004311, -0.003789) p01 = -0.0006835 (-0.001057, -0.0003101) Hg:
f(x,y) = p00 + p10*x + p01*y
p00 = 732.3 (669.5, 795)
p10 = -0.0185 (-0.0221, -0.01489) p01 = -0.02318 (-0.02833, -0.01802) Ni:
f(x,y) = p00 + p10*x + p01*y
p00 = 22.07 (21.71, 22.44)
p10 = -0.0003019 (-0.000323, -0.0002809) p01 = -5.326e-005 (-8.337e-005, -2.315e-005) Pb ：
f(x,y) = p00 + p10*x + p01*y
p00 = 111 (109.1, 112.8)
p10 = -0.003123 (-0.003231, -0.003016) p01 = -0.0004823 (-0.0006362, -0.0003284) Zn:
f(x,y) = p00 + p10*x + p01*y
p00 = 398 (387, 409.1)
p10 = -0.01267 (-0.0133, -0.01203) p01 = -0.003615 (-0.004525, -0.002706) 7.2各种金属浓度的等值线
利用Surfer 画出不同金属元素的等值线图，建立等值线与MATLAB 双向逼近模型，从而确定所得的污染源的位置坐标。
Surfer 软件等值线如下所示：
15010050
As
Cd
15010050
Cr Cu
[***********]005000
5000
10000
15000
20000
25000
Hg
Ni
Pb Zn
通过计算，可以确定8种重金属污染源的主要范围为： 1、污染源（1862.46, 2827.52），半径1487.60； 2、污染源（3840.71, 5253.45），半径1190.08； 3、污染源（9397.57, 4475.32），半径732.35； 4、污染源（12731.69, 2919.06），半径1281.62； 5、污染源（13620.79, 9693.37），半径640.81； 6、污染源（18021.83, 9968.01），半径732.35； 7、污染源（21489.31, 11478.50），半径686.58。(单位：Km)
8问题四的模型建立与求解
8.1模型优缺点：
优点：在问题三中，通过Surfer 软件绘制出的重金属等值线与MATLAB 软件绘制出的曲面拟合函数曲面逼近模型，能够较为精确地反映出污染源的坐标；
缺点：根据现有条件建立的基于两款强大的应用型软件的双向逼近模型，并没有考虑因温度，地质灾害，气流扰动，海拔高度变化等因素而造成的影响，无法根据实时的条件进行调试，还需要进行大量数据统计之后，才能进行数据等值线的绘制和数据函数的曲面函数图像拟合。 8.2需收集信息
为更好地研究城市地质环境的演变模式， 本文认为还应当收集地质灾害，水土流失情况，人类活动强度，植被覆盖率，地标坡度等因素的信息。 8.3模型的建立
基于收集到的信息，本文只提取三个因素予以考虑，将α, β, χ分别表示为地表坡度，植被覆盖率和年降水量， 在第三问的基础上，将拟合函数中x 的幂
指数升级为2次，以更好的表示出数据的特点。用G(α, β, χ) 这一综合指标函数考虑三个因素对于地表重金属含量浓度的影响。现建立函数方程如下： F(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y+G(α, β, χ) ；
式中：F(x,y)代表重金属含量浓度，p00， p10， p20，p11，分别代表x,y 的系数。
参考文献
[1]L.Anecological risk index for aquatic pollution control:a sedimen-tological approach[J].Water Research,1980,14(8):975-1001.
[2]何云峰，朱广伟，陈英旭，等，运河杭州段沉积物中重金属的潜在生态风险研究[J],浙江大学学报：农业与生命科学版，2002，28（6）；669-674；
[3]縢彦国，虞先国，倪师军，等，应用地质累积指数评价沉积物中重金属污染选择地球化学背景的影响[J],环境科学与技术，2002，25（2）：7-9
[4]縢彦国，虞先国，倪师军，等，应用地质累积指数评价攀枝花地区土壤重金属污染[J],重庆环境科学，2002，24（4）；25-31。
[5]FORSTNER U.AHLF W.CALMANO W, er al.sediment criteria development-contributions from environmental geochemistry to water quality management[C]//HELING D,ROTHE P,FORSTNER U,et al,Sediments and environmental geochenmistry selected aspects and case histories. Springer-Verlag, Berlin Heideberg,1990,12(6);311-338;
[6]ANON D L ,Ein fliessge was seroko system[J],Regierung sprasidium Giessen, Niedernhausen, 1994,5(2):87。
[7]房存金，土壤中主要重金属污染物的迁移转化及治理，2010年8月，39卷第4期，当代化工。
附录
附录1：各区域各点的I PL
附录2：各元素的曲面拟合图像
Cd
Cr
Cu
Hg
Ni
Pb
Zn
附录3：MATLAB 拟合曲面函数时源程序代码 clc clear close all
D=ZuoBiao; M=JinShu;
As=M(:,2);Cd=M(:,3);Cr=M(:,4);Cu=M(:,5);Hg=M(:,6);Ni=M(:,7); Pb=M(:,8);Zn=M(:,9); 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);
%Asi=linspace(min(As),min(As),100); [xi,yi]=meshgrid(xi,yi); zi=griddata(x,y,z,xi,yi); ci=griddata(x,y,c,xi,yi);
Asi=griddata(x,y,As,xi,yi);Cdi=griddata(x,y,Cd,xi,yi); Cri=griddata(x,y,Cr,xi,yi); Cui=griddata(x,y,Cu,xi,yi); Hgi=griddata(x,y,Hg,xi,yi);
Nii=griddata(x,y,Ni,xi,yi);Pbi=griddata(x,y,Pb,xi,yi); Zni=griddata(x,y,Zn,xi,yi);
marker={'*','o','s','^','p'}; color={'k','r','y','c','b'};
mat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn'};
str={'等高线',' 生活区',' 工业区',' 山林区',' 交通区',' 绿地区'};
for j=1:8
% 等高线图形 figure
contourf(xi,yi,zi,0:10:500);
% set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2) % clabel(C,h,[0:10:50,50:50:300,300:100:500]); title(['金属 ',mat{j}, ' 二维等高线分布图']) xlabel('X') ylabel('Y')
colormap summer
colorbar grid on
hold on for i=1:5
loc=c==i;
plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i}); end
legend(str,'location','best') for k=1:length(x)
text(x(k)-200,y(k)+200,num2str(M(k,j+1)),'fontsize',8); end end
% 三维体图 figure
h=surf(xi,yi,zi); set(h,'cdata',ci); colormap hsv
title('三维图立体（颜色条表示分类）') xlabel('X') ylabel('Y') 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
str{1}='三维图';
legend(str,'location','best')
附录5：各区域各金属元素的相关系数
生活区
工业区
山区
交通区
公园绿地区
城市表层土壤重金属污染分析
摘要
问题一中本文利用Surfer 软件，根据提供的数据，绘制出8种重金属的空间立体含量图。再运用地质累积指数法，用I PL 以及I PLzone 表示出单一金属污染负荷指数和8种重金属元素的区域污染负荷指数。由3维立体图可看出As 砷元素的污染程度在1,2,5 区的含量较高，污染较为严重；Cd 镉元素在1,2,4,5区污染较为严重；Cr 铬元素在1,4区污染严重；Cu 铜元素在五个区分布普遍较高，尤其以2,4区污染最为严重；Hg 元素在2区污染较为严重。根据I PLzone 得知，工业区属强污染，其他区域均属于中等污染。
问题二将8种重金属含量进行了主成分分析。通过对比不同金属对环境污染贡献值的差异，分析出污染的主要可能原因，生活区：As ，Cu ，Ni 为主要污染物质，原因可能与居民生活垃圾或该区交通污染有关；工业区：Cr ，Cu ，Hg ，Pb ，Zn 为主要污染物质，导致原因可能是各类工厂废气、废水、废渣的排放。山区：Cr ，Ni 和Zn 为主要污染物质，原因可能是冶金，采矿等活动造成的。交通区：Cr ，Cu 和Ni 为主要污染成分，原因可能是汽车尾气，轮胎老化，车体磨损。公园绿地区：Cr ，Pb 和Zn 为主要污染物质，原因可能是由于汽车尾气排放和工业活动导致公园土壤重金属含量的增加。
问题三中通过查阅资料，分析出不同重金属在地表迁移传播的特点，并利用SPSS 计算出金属间的相关性系数，得出Pb 与Cd 和Hg,Cu 与Cr 具有显著相关性，判断它们可能出自相似污染源或同一污染源。对8种金属浓度进行曲面拟合，将拟合方程结合等值线图，建立了双向逼近模型，确定出主要污染源的位置是：（1）圆心（1862.46, 2827.52），半径1487.60；（2）圆心（3840.71, 5253.45），半径1190.08，单位为Km 。
问题四中将地表重金属含量作为预测地质环境的指标，通过对第三问条件的延伸，在收集了地表坡度、植被覆盖率和年降水量三种数据的基础上，建立F(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y+G(α, β, χ) ；其中G(α, β, χ) 为三个影响因素集成的函数。
关键词：重金属污染 地质累积指数 主成分分析 曲面拟合 双向逼近模型
1问题重述
随着人口急剧增长、工农业和科学技术水平的不断提高，人类正以前所未有的规模和强度开发自然资源和破坏生态环境。其中重金属污染是引起的主要问题之一。
现对某城市城区土壤地质环境进行调查。将该城市划分为五个区，分别是生活区、工业区、山区、主干道路区及公园绿地区，并标号为1、2、3、4、5。将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
本题目要求：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
2基本假设
1 不发生大规模泥石流，山体滑坡等地质灾害，不发生旱涝灾害等自然不可逆灾害现象。
2 假设该城区土层表面土质均一稳定。 3 取样期间，植被覆盖率没有明显变化；
4 假设取样区域内没有河流径流，不会对污染物造成不可预测的扩散。 5 假设所取样区域不发生大规模降水，且气候变化平缓，温度变化不大。
3主要变量符号说明
4问题的分析
此题目中，第一问要求求解8种金属在该区域中的空间分布，需要对数据种的横从坐标与8种不同的重金属进行三维效果处理，从而能够可以从空间上给出其分布情况。
而第二问，需要求解重金属污染的主要原因，这就要求通过对现有的数据进行处理，本文利用了主成分分析方法，对污染原因进行了说明。
第三问，提出根据不同金属在土壤中横向的扩散传播情况进行模型建立，我们通过相关系数的模型建立，分析处理，得出污染源的坐标位置。
第四问，在第三问的基础上，综合考虑三个方面的因素（地表坡度、植被覆盖率、年降水量）进行分析，增加G(α, β, χ) 这一函数来表征对重金属含量的影响。
5问题一的模型建立与求解
5.1 模型的建立
5.1.1各金属元素的空间分布
在利用Surfer 软件的条件下，进行x,y, 以及污染物浓度的绘制，形象直观的得到关于横纵坐标与8种污染物浓度的三维图像：
As
Cd
Cr
Cu
Hg
Ni
Pb
图1
Zn
5.1.2分析污染程度
通过建立污染负荷指数模型，对8种重金属进行对不同地区的污染程度大小的计算求解。污染指数负荷法是Tomlinson 等在从事重金属研究水平的分级研究中提出来的一种评价的方法，该方法被广泛的应用于土壤和河流沉积物重金属污染的评价[1 2]。某一点的污染负荷指数公式如下：
F i =C i /C o i I P L =
式中，Fi 为元素i 的污染系数；Ci 为元素i 的实测质量分数（mg/Kg）；C oi 为元素i 的评价标准，即背景值，一般选用全球页岩平均值作为重金属的评价标准（mg/Kg）；I PL 为某一点污染负荷指数；n 为评价元素的个数。
某一区域的污染负荷指数(I PLzone ) 为：
I PLzone =
式子中，I PLzone 为区域污染负荷指数，n 为评价点的个数（即采样点的个数）。
5.2模型求解
根据地积累指数法的公式，分别算出各区域各元素的I PL ，并统计结果如下：
根据上表并不能完全了解清楚每一个区域的不同程度的污染物的分布情况，通过计算I PLzone ，得出各区域的污染负荷指数。以下是五区的污染程度分级差异情况：
5.3模型补充
通过以上计算分析，得出的是用8种重金属物质的浓度来衡量不同功能区的污染程度，为了补充前一模型的不足，下面用地址累积指数法来进行单一重金属对于不同功能区的衡量：
地质累积指数不仅反映了重金属分布的自然变化特征，而且可以判别认为活动对人为环境的影响，是区分人为活动影响的重要参数[3-4]。其模型如下：
I geo
⎡C n ⎤=log 2⎢⎥
⎣1.5BE n ⎦
其中，C n 表示样品中元素n 的浓度，BE n 表示基线浓度，1.5为修正参数，通常用来表征沉积特征、岩石、地质及其他影响。
一般而言，地质累积指数有2种分级，Forstner 等将其分为7级[5]，而Anon 则分为5级[6]，不同的级别分别代表不同的重金属污染程度，本文采用Anon ，见表3。
表3 地质积累指数分级情况
Igeo 值（Anon ）
5
级别 1 2 3 4 5
污染程度 无污染 无污染到轻度污染
中度污染 中度污染到强度污染
极强污染
6问题二的数据分析与求解
6.1该市土壤重金属主成分分析
为了分析重金属污染的主要原因，运用SPSS 软件对上述不同功能区重金属的各项指标进行主成分分析，数据处理及分析结果见表4：
以上数据结果说明：利用SPSS 软件得出的第一组数据，即生活区的主成分提取出4个（λc >0.8），贡献率分别为45.199%，14.165%，13.432%，10.083%；对表格进行分析，可以得到第一主成分Cd 和Cu 、Pb 有较高的载荷，反映出了住宅区内居民产生的生活垃圾或与该区交通污染有关。第二主成分中Zn 具有较高的载荷，反映出了生活区可能存在如汽车轮胎老化等现象。其他两种主成分较为次要，这里不加以分析。
工业区提取出了2个主成分进行分析（λc >0.8），贡献率分别是65.670%，15.786%，对表格进行分析，可以得到第一主成分中Cr 和Cu 及Hg 、Pb 、Zn 有一定程度的高载荷，反映出这一工业区污染主要原因有电镀、采矿、冶炼、燃料、电池和化学工业等。第二主成分中As 有较高载荷，说明污染来自于各种化学制剂，包括农药，杀虫剂在内的物质。
山区提取出了3个主成分进行分析（λc >0.8），贡献率分别是38.022%，25.445%，19.362%，对表格进行分析，可以发现Cr 和Ni 及Zn 有较高载荷，说明污染原因有汽车轮胎磨损以及煤燃烧产生的粉尘、烟尘中均含有有毒重金属。第二主成分As 的载荷较高；第三主成分中Cu 和Hg 的载荷较高，说明污染与金属加工，日常用品制造等物质有密切关系。
交通区提取出了3个主成分进行分析（λc >0.8），贡献率分别是46.883%，16.084%，12.424%,对表格分析，得知第一主成分中Cr,Cu,Ni 具有较高载荷；第二主成分中Hg 具有较高载荷，说明其原因可能有金属冶炼和燃煤；第三主成分中As 具有很高载荷，达到0.939。
绿地区提取出了3个主成分进行分析（λc >0.8），贡献率分别是48.840%，20.182%，13.202%，对表格进行分析，可以得出第一主成分中Cr ，Pb,Zn 有较高的载荷，说明污染原因可能是工业产生的废物，废气，废渣；第三主成分中，Hg 具有较高的载荷，达到了0.889。
7问题三的模型建立与求解
在这一问中，通过查阅资料[7]了解到8种金属元素在降雨量不发生剧烈改变，气候变化平缓的情况下，大多数不会产生较为明显的横向迁移传播，本文用
SPSS 软件，对所给数据进行了一定量的处理，得到8种金属之间相关性系数，从中得出了某几种金属之间的相关性结论，相关性系数表格见附录，这里不予以赘述。
7.1. 各金属的相关性及浓度的曲线拟合
生活区：As 与Cu; As与Ni; Cr与Ni; Cd与Pb; 具有较强的相关性，说明其可能来自于同一污染源。
工业区： Cu 与Cr ；Cu 与Hg ；Pb 与Zn; As与Ni ；Cr 与Hg; Cr 与Ni; Cr与Pb; Cu与Zn 之间存在较强的相关性，说明这几组重金属物质之间存在一定相似的污染源。
山区：Cr 与Ni ；Zn 与Cd ；Ni 与Pb ；具有较强的相关性，说明其污染源可能源于同一个范围。
交通区：Cr 与Cu; Cr 与Ni; Cu 与Ni; Cd 与Zn; 具有较强的相关性，说明其污染源相似。
公园绿地区：Cr 与Ni; Cu与Zn; 具有较强的相关性，说明其污染源相似。 本文利用Matlab 软件进行重金属浓度曲面拟合，拟合得出的曲面如下图所示(这里只表示As 元素的曲面拟合，其他元素见附录2) ：
As
拟合出的函数的方程式如下： As:
f(x,y) = p00 + p10*x + p01*y
p00 = 7.777 (7.663, 7.891)
p10 = -0.0001683 (-0.0001748, -0.0001617) p01 = 2.846e-005 (1.905e-005, 3.786e-005) Cd:
f(x,y) = p00 + p10*x + p01*y
p00 = 488.8 (480.4, 497.2)
p10 = -0.01223 (-0.01271, -0.01175) p01 = -0.00194 (-0.00263, -0.00125) Cr:
f(x,y) = p00 + p10*x + p01*y
p00 = 85.53 (82.74, 88.31)
p10 = -0.001988 (-0.002148, -0.001828) p01 = -0.0005593 (-0.0007881, -0.0003304) Cu:
f(x,y) = p00 + p10*x + p01*y
p00 = 113.6 (109, 118.1)
p10 = -0.00405 (-0.004311, -0.003789) p01 = -0.0006835 (-0.001057, -0.0003101) Hg:
f(x,y) = p00 + p10*x + p01*y
p00 = 732.3 (669.5, 795)
p10 = -0.0185 (-0.0221, -0.01489) p01 = -0.02318 (-0.02833, -0.01802) Ni:
f(x,y) = p00 + p10*x + p01*y
p00 = 22.07 (21.71, 22.44)
p10 = -0.0003019 (-0.000323, -0.0002809) p01 = -5.326e-005 (-8.337e-005, -2.315e-005) Pb ：
f(x,y) = p00 + p10*x + p01*y
p00 = 111 (109.1, 112.8)
p10 = -0.003123 (-0.003231, -0.003016) p01 = -0.0004823 (-0.0006362, -0.0003284) Zn:
f(x,y) = p00 + p10*x + p01*y
p00 = 398 (387, 409.1)
p10 = -0.01267 (-0.0133, -0.01203) p01 = -0.003615 (-0.004525, -0.002706) 7.2各种金属浓度的等值线
利用Surfer 画出不同金属元素的等值线图，建立等值线与MATLAB 双向逼近模型，从而确定所得的污染源的位置坐标。
Surfer 软件等值线如下所示：
15010050
As
Cd
15010050
Cr Cu
[***********]005000
5000
10000
15000
20000
25000
Hg
Ni
Pb Zn
通过计算，可以确定8种重金属污染源的主要范围为： 1、污染源（1862.46, 2827.52），半径1487.60； 2、污染源（3840.71, 5253.45），半径1190.08； 3、污染源（9397.57, 4475.32），半径732.35； 4、污染源（12731.69, 2919.06），半径1281.62； 5、污染源（13620.79, 9693.37），半径640.81； 6、污染源（18021.83, 9968.01），半径732.35； 7、污染源（21489.31, 11478.50），半径686.58。(单位：Km)
8问题四的模型建立与求解
8.1模型优缺点：
优点：在问题三中，通过Surfer 软件绘制出的重金属等值线与MATLAB 软件绘制出的曲面拟合函数曲面逼近模型，能够较为精确地反映出污染源的坐标；
缺点：根据现有条件建立的基于两款强大的应用型软件的双向逼近模型，并没有考虑因温度，地质灾害，气流扰动，海拔高度变化等因素而造成的影响，无法根据实时的条件进行调试，还需要进行大量数据统计之后，才能进行数据等值线的绘制和数据函数的曲面函数图像拟合。 8.2需收集信息
为更好地研究城市地质环境的演变模式， 本文认为还应当收集地质灾害，水土流失情况，人类活动强度，植被覆盖率，地标坡度等因素的信息。 8.3模型的建立
基于收集到的信息，本文只提取三个因素予以考虑，将α, β, χ分别表示为地表坡度，植被覆盖率和年降水量， 在第三问的基础上，将拟合函数中x 的幂
指数升级为2次，以更好的表示出数据的特点。用G(α, β, χ) 这一综合指标函数考虑三个因素对于地表重金属含量浓度的影响。现建立函数方程如下： F(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y+G(α, β, χ) ；
式中：F(x,y)代表重金属含量浓度，p00， p10， p20，p11，分别代表x,y 的系数。
参考文献
[1]L.Anecological risk index for aquatic pollution control:a sedimen-tological approach[J].Water Research,1980,14(8):975-1001.
[2]何云峰，朱广伟，陈英旭，等，运河杭州段沉积物中重金属的潜在生态风险研究[J],浙江大学学报：农业与生命科学版，2002，28（6）；669-674；
[3]縢彦国，虞先国，倪师军，等，应用地质累积指数评价沉积物中重金属污染选择地球化学背景的影响[J],环境科学与技术，2002，25（2）：7-9
[4]縢彦国，虞先国，倪师军，等，应用地质累积指数评价攀枝花地区土壤重金属污染[J],重庆环境科学，2002，24（4）；25-31。
[5]FORSTNER U.AHLF W.CALMANO W, er al.sediment criteria development-contributions from environmental geochemistry to water quality management[C]//HELING D,ROTHE P,FORSTNER U,et al,Sediments and environmental geochenmistry selected aspects and case histories. Springer-Verlag, Berlin Heideberg,1990,12(6);311-338;
[6]ANON D L ,Ein fliessge was seroko system[J],Regierung sprasidium Giessen, Niedernhausen, 1994,5(2):87。
[7]房存金，土壤中主要重金属污染物的迁移转化及治理，2010年8月，39卷第4期，当代化工。
附录
附录1：各区域各点的I PL
附录2：各元素的曲面拟合图像
Cd
Cr
Cu
Hg
Ni
Pb
Zn
附录3：MATLAB 拟合曲面函数时源程序代码 clc clear close all
D=ZuoBiao; M=JinShu;
As=M(:,2);Cd=M(:,3);Cr=M(:,4);Cu=M(:,5);Hg=M(:,6);Ni=M(:,7); Pb=M(:,8);Zn=M(:,9); 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);
%Asi=linspace(min(As),min(As),100); [xi,yi]=meshgrid(xi,yi); zi=griddata(x,y,z,xi,yi); ci=griddata(x,y,c,xi,yi);
Asi=griddata(x,y,As,xi,yi);Cdi=griddata(x,y,Cd,xi,yi); Cri=griddata(x,y,Cr,xi,yi); Cui=griddata(x,y,Cu,xi,yi); Hgi=griddata(x,y,Hg,xi,yi);
Nii=griddata(x,y,Ni,xi,yi);Pbi=griddata(x,y,Pb,xi,yi); Zni=griddata(x,y,Zn,xi,yi);
marker={'*','o','s','^','p'}; color={'k','r','y','c','b'};
mat={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn'};
str={'等高线',' 生活区',' 工业区',' 山林区',' 交通区',' 绿地区'};
for j=1:8
% 等高线图形 figure
contourf(xi,yi,zi,0:10:500);
% set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2) % clabel(C,h,[0:10:50,50:50:300,300:100:500]); title(['金属 ',mat{j}, ' 二维等高线分布图']) xlabel('X') ylabel('Y')
colormap summer
colorbar grid on
hold on for i=1:5
loc=c==i;
plot(x(loc),y(loc),marker{i},'markerfacecolor',color{i},'MarkerEdgeColor',color{i}); end
legend(str,'location','best') for k=1:length(x)
text(x(k)-200,y(k)+200,num2str(M(k,j+1)),'fontsize',8); end end
% 三维体图 figure
h=surf(xi,yi,zi); set(h,'cdata',ci); colormap hsv
title('三维图立体（颜色条表示分类）') xlabel('X') ylabel('Y') 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
str{1}='三维图';
legend(str,'location','best')
附录5：各区域各金属元素的相关系数
生活区
工业区
山区
交通区
公园绿地区