城表层土壤重金属污染

1、城市表层土壤重金属污染分析摘要随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。本文针对某城市八种重金属进行污染浓度分布建模分析，研究其不同地区受污染程度和污染原因，并确定污染源。针对问题一，首先对附件一中数据进行处理，然后利用MATLAB插值拟合曲线，得到八种重金属元素在该城市的空间分布图，运用单因子污染指数法和内梅罗综合污染指数法分别确定单种重金属元素的环境污染指数和元素的综合污染指数后进行污染等级评价，最总得出结论，在五大功能区中，工业区污染较为严重，山区污染污染程度不高，较清洁，且该城市受Zn、Hg、Cu元素污染较严重。针对问题二，重金属元素污染之间存在

2、必然的联系，同一区域的污染元素之间有较强的相关性，因此我们采用因子分析法来研究多个变量的相关性。通过对原始数据的处理，将成因进行归类，总结出几个比较客观的成因线索，最终导出成因结果。本文使用SPS漱件对八种重金属元素进行相关性分析，最终得出结论：该城市的重金属污染主要原因为以下方面：工业污染源、交通污染源、燃煤污染源、农药污染源。针对问题三，讨论重金属元素在土壤中的传播问题，我们建立抛物型的偏微分方程模型，解出这个偏微分方程的Cauchy问题的解。由于污染源污染范围有限，所以我们挑选出污染较为严重的一处，通过反演的方法拟合出污染源的位置(x0,y0,z0),本文仅针对Zn、Hg、Cu这类污染较

3、严重的元素，得到一处或多处污染源。针对问题四，我们首先提出模型的优缺点，然后根据优缺点提出模型的改进建议。考虑到各区域某一重金属元素污染下的地质累积污染程度强弱状况、地质演变因素，我们收集材料，在条件允许下，建立优化模型。例如，考虑地质累积指数法，建立更为精确的污染程度分析，使得问题一的结果更加可信；考虑大气干湿度对重金属土壤污染扩散的影响，建立相关数学模型，能够估计n年后该城市的土壤污染程度。关键词：单因子污染指数法、内梅罗综合污染指数法、因子相关性分析法、抛物型的偏微分方程模型一、问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异

4、常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土(010厘米深度)进行取样、编号，并用GPS己录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作

5、为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求你们通过数学建模来完成以下任务：(1)给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2)通过数据分析，说明重金属污染的主要原因。(3)分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4)分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二、模型假设1、假设该城市各个功能区相互独立，不交叉。2、假设该

6、城市仅受这八种重金属污染，无其他污染物干扰3、假设八种重金属污染源之间无相互干扰。4、假设八种金属污染物在分析期间分布浓度无明显变化。三、符号说明Pi=1,2,3,4,5,6,7,8环境质量指数Cii=1,2,3,4,5,6,7,8评价对象污染物i的实测质量分数Si=1,2,3,4,5,6,7,8污染物i的评价标准，一般取二类标准唉内梅罗综合污染指数Pmax元素环境质量指数最大值Pave兀素环境质量指数平均值u（x,y,z,t）t时刻（x,y,z）处的一种重金属元素浓度a2,b2,c2分别为沿x,y,z方向上的扩散系数k2衰减系数（x»另,4,5）/=1,2,n取样点的横坐标、纵坐标

7、、高度及一种元素的浓度四、问题一重金属污染分布及污染程度4.1 问题分析利用附件一给出的取样点位置以及所属功能区，附件二给出的八种重金属在对应位置的浓度。运用指数法计算计算出八种重金属元素在指数法条件下的各个指标值，从而确定区域污染程度。初步分析计算得到如下指标：表一：该城市八种重金属含量参数统计表73?r背景值最大值最小值平均值标准差As(Ng/g)3.630.131.615.683.02Cd(Ng/g)1301.620.040.30.22Cr(Rg/g)31920.8415.3253.5169.89Cu(Ng/g)13.22528.482.2955.02162.66Hg(邑g/g)3516

8、0.010.31.63Ni(/g)12.3142.54.2717.269.93Pb(Ng/g)31472.4819.6861.7449.98Zn(Ng/g)693760.8232.86201.20338.704.2 数据处理1)首先将中金属浓度单位统一为Ng/g。2)根据单因子污染指数法求出各重金属元素在不同浓度的环境质量指数R,根据标准得出各个重金属元素质量评价等级。3)利用R求得元素环境质量指数最大值Pmax和平均值Rave,根据综合污染指数法的计算公式，求得内梅罗综合污染指数，确定综合元素质量评价等级。4.3 模型建立4.3.1 单因子污染指数法单因子污染指数法是指在某一污染物影响下的环

9、境污染指数。其能够反映污染物在该环境中的污染程度。根据下式计算单因子污染程度等级，并进行分级。CR=，i=123,4,5,6,7,8(1)Si其中，国家土壤环境质量标准的二级标准经查阅为表二：国家二级标准单因子污染评价各污染标准(mg/kg)toKAsCdCrCuHgNiPbZnSi250.3200500.350250200表三：城市表层土壤单因子污染分级标准污染分指数R<11EPi<22<R<3R>3质量等级清洁潜在污染轻污染重污染4.3.2 内梅罗综合污染指数法前文对该城市单个元素的污染做了评价，但不能反映各元素共同作用对城市表层土壤的复合污染。鉴于此为突出环

10、境要素中含量最大的污染物随环境质量的影响，因此采用内梅罗综合污染指数法进行综合评价，具公式为：D_Pmax+Pave(P=q2(内梅罗综合污染指数分级标准见下表：表四：内梅罗综合污染指数分级标准等级1111111VV综合指数凄-0.70.7；P综.11:二嗪-22:二P综-3嗪3污染等级清洁警戒限轻度污染重度污染重污染4.4 模型求解根据附件一中给出的三维空间坐标进行插值处理，得到该城市功能区分布图。16000140001200014JQQD加0口60004Q002000X1。"图一该城市功能区分布图4.4.1 重金属元素的空间分布图二As的分布图图三Cd的分布图XW依据附件所给数据

11、，利用MATLABS?,得出八种重金属元素在该城市的空间分布等高线图，如下，图四Cr的分布图图五Cu的分布图图六Hg的分布图图七Ni的分布图图八Pb的分布图图九Zn的分布图4.4.2不同区域重金属污染程度分析依据单因子污染指数公式，得出各重金属元素污染指数，应用单因子污染指数法对这八种重金属元素进行污染程度评价，得到下表：表五：城市表面土壤重金属单因子污染指数统计表AsCdCrCuHgNiPbZn最大值1.215.404.6050.5753.332.851.8918.80最小值0.060.130.080.050.030.090.080.16平均值0.231.010.271.101.000.35

12、0.251.01污染程度清洁潜在污染清洁潜在污染清洁清洁清洁潜在污染表六：城市表面土壤重金属单因子污染指数统计表污染程度AsCdCrCuHgNiPbZn清洁3148331617226031488207潜在4113360355228134污染轻污染0116085160360重污染170280018从表七看出，该城市表层土壤污染并不严重，但有部分地区受Cd,Cu,Hg,Zn元素重污染，且Zn重污染区最多达到18个。卜面给出不同功能区的受污染程度分析如下，表七：生活区重金属污染程度表AsCdCrCuHgNiPbZn最大值0.763.067.442.272.212.822.8910.98最小值0.16

13、0.580.210.320.080.220.700.51平均值0.421.620.741.290.580.461.131.74污染清洁潜在清洁潜在清洁清洁潜在潜在程度污染污染污染污染表八：工业区重金属污染程度表AsCdCrCuHgNiPbZn最大值1.693.132.869.082.221.172.746.75最小值0.110.760.170.420.080.110.890.66平均值0.491.950.851.721.360.501.262.01污染程度清洁潜在污染清洁潜在污染潜在污染清洁潜在污染轻污染表九：山区重金属污染程度表AsCdCrCuHgNiPbZn最大值0.732.151.162

14、.051.380.521.892.10最小值0.120.270.180.080.060.090.130.39平均值0.271.000.430.580.270.350.280.84污染程度清洁清洁清洁清洁清洁清洁清洁清洁表十:交通区重金属污染程度表AsCdCrCuHgNiPbZn最大值3.033.899.215.7615.081.001.6813.87最小值0.110.330.170.410.060.150.630.48平均值0.391.800.631.440.930.431.111.83污染程度清洁潜在污染清洁潜在污染清洁清洁潜在污染潜在污染表十一：公园绿地区重金属污染程度表AsCdCrCuH

15、gNiPbZn最大值0.783.041.062.272.870.731.895.96最小值0.180.650.180.300.070.190.770.44平均值0.421.460.480.950.590.381.101.32污染清洁潜在清洁清洁清洁清洁潜在潜在程度污染污染污染依据内梅罗综合污染指数法的计算公式得到五个功能区的内梅罗综合污染指数P综（取平均值），建立表格：表十二：不同区域重金属污染程度功能区生活区工业区山区交通区公园绿地区嗪1.652.010.931.891.35污染程度轻度污染中度污染警戒限轻度污染轻度污染依据上面几个表格得出结论，工业区污染状况最为严重，其次是生活区和交通区，

16、然后是公园绿地区，而山区污染状况最好，污染程度都为清洁。五、问题二重金属污染主要原因5.1 问题分析问题二要求我们通过数据分析得出重金属污染的主要原因，由问题一得出八种元素在该城市空间分布以及不同功能区的污染程度。在不同功能区中，重金属的来源可能不同，但也可能来源于同一污染源，因此每处土壤中的污染元素存在一定的相关性，不互相独立，且相互影响。因此使用因子分析算法来解决这个问题，通过因子分析法并利用SPS欹件各个因子空间分布，确定各因子集区域，从而得出污染主要原因。5.2 模型建立及求解通过对因子分析算法的了解，我们发现该城市重金属元素含量的数据特征符合因子分析的要求。在该城市环境质量评价基础上

17、，对存在污染的土壤的污染来源及分布进行分析。运用SPSS22软件进行因子分析，得到如下结论：首先给出八种重金属元素的相关系数矩阵：表十三:八种重金属污染相关系数表AsCdCrCuHgNiPbZnAs10.2480.180.1510.060.3130.2840.24Cd0.24810.3130.3670.2530.3160.6360.392Cr0.180.31310.5090.0860.7160.340.386Cu0.1510.3670.50910.410.4870.4940.355Hg0.060.2530.0860.4110.0950.2870.181Ni0.3130.3160.7160.48

18、70.09510.2920.427Pb0.2840.6360.340.4940.2870.29210.454Zn0.240.3920.3860.3550.1810.4270.4541由表可见，Ni和Cr的相关性最好，相关系数最大为0.716。其次为Pb和Cd,相关系数为0.636。以下依次是Cu和Cr相关系数为0.509,而其它元素之间的相关性并不好。从成分上分析，相关信号的元素在成因和来源上有一定的联系。其原因分析：1）工业生产中，锲铭合金是常用合金。2）铅和镉元素在蓄电池的制造中会以气态形式散发。且他们都是工业排放的废物的主要元素。从问题一中各元素在不同功能区中的污染浓度看出，铅元素不仅在

19、工业区，交通区和生活区都占有较大比重，说明汽车尾气污染也很严重。3）绿地公园区的污染不大，其重金属遗留原因应是农药残留。因子分析的关键在于利用相关系数矩阵求出相应的因子特征值和累计贡献率，依然可用用SPSS22软件计算得出，结果如下：明的燮昊数女恩言十元件起始特徵值摄取平方和载入循璟平方和载入燮昊的累加*就十燮昊的累加燮昊的累加13.56044.50044.5003.56044.50044.5002.51831.47831.47821.15014.37758.8771.15014.37758.8772.19227.40058.8773.96512.06370.9414.7689.59680.5

20、375.5787.22087.7566.4325.39993.1567.3013.76996.9248.2463.076100.000撷取方法：主元件分析从表中看出在累计方差达87.756%勺前提下，分析得五个主因子，可以看出五个主因子为源资料提供了87.756%的信息，满足因子分析原则。因子分析的主要目的是将具有相近的因子荷载的各个变量置于一个公因子之下，正交方差最大旋转使每一个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小，以便对因子的意义作出更合理的解释，输出结果见表。旋傅元彳矩障a元件12As(gg/g).452.130Cd(ng/g).350.679Cr(gg/g).

21、846.149Cu(gg/g).488.591Hg(ng/g)-.135.776Ni(gg/g).883.088Pb(gg/g).369.739Zn(gg/g).551.431根据软件得出结果，基于旋转理论：变量与某一个因子的联系系数绝对值（载荷）越大，则该因子与变量关系越近。正交因子解说明：因子1为Cr、Ni和Cu的组合，这说明这几种土壤重金属污染物可能是同一来源或相似来源；因子2为Cd和Pb的组合表明两者可能有相似的来源。Hg在土壤的污染严重，污染区域大多处于某一工业区或某大型污染企业，其来源比较单一。表十四：中国土壤重金属污染来源表重金属CrHgAsPbNiCrHgAsPbCrCdNiP

22、bCuZnHgCdPbNiZnHgCrCd来源矿产开采、冶炼、加工排放的废气、废水和废渣煤和石油燃烧过程中排放的飘尘电镀工业废水塑料、电池、电子工业排放的废水工业排放的废水染料、化工制革工业排放的废水汽车尾气Pb农药、肥料AsCuCd最后给出结论：该城市重金属污染来源主要有四个方面：交通污染源，燃煤污染源，工业污染源，化肥农药污染源，其中Zn、Pb、Cu主要来自汽车尾气的排放，即交通污染源;Hg主要来自工业燃煤及居民采暖用煤;As污染主要来自化工、冶金、炼焦、火力发电、造纸、玻璃、毛革、电子工业等工业污染；CdCr、Ni污染主要来自金属冶炼。六、问题三重金属污染源确定6.1 问题分析问题三要求

23、我们确定重金属元素的污染源，因此我们要先了解污染物的传播特征，考虑到金属元素在土壤中并不固定，存在迁移和变化规律，建立迁移扩散模型。再此，我们建立非线性抛物型偏微方程作为模型来定性解决问题，考虑到扩散，采用衰减的扩散过程数学模型。随后利用反演的方法逐步找到污染源。6.2 模型建立6.2.1 有衰减的扩散过程数学模型考虑到元素扩散和衰减的合作用，有物质不灭定理得：.2-2.2(1)二u2;U2U21u2二a-2b2c-2-ku2222.t二x；ycz方程（1）是常系数线性抛物方程，它就是有衰减的扩散过程的数学模型（对于具体问题，尚需匹配求解）。在本题中就是要通过初始条件t=0时污染源的条件来进行

24、求解。6.2.2线性抛物方程的Cauchy问题求解设污染源在点（X0,y0,Z0）处，则此扩散模型满足Cauchy问题：一.2.2.2一u20u,2：u2;u.2(2)二a2b2c-2-ku二t二x二y一zu(x,y,z,0)=M、(xX0)、(y-y。)、(z-z0)0.x=0其中S(x)为Dirac函数，即6(x)=4,d(x)dx=1,二，x=0解上述方程得到,，MMu(x,y,z,t)=8二tabc*二t22expj(x-%)(yy0)(z-4)224a2t4b2t4c2t2-k2tj(3)6.2.3参数估计在解（3）中存在四个未知参数，a、b、c、k,他们分别是扩散与衰减过程中的扩散

25、系数和衰减系数的算术平方根。令取样时刻的t=1,对（3）式两边取对数得222M(x-x0)(y-y0)(z-z0),2X-、lnu(x,y,z,t)=ln2-2-+2-+-=0/-+k2t5(4)8二tabc*二t4a2t4b2t4c2t22x=(2SzxL”(XzMz(z-z0)21=，:二,M一,2；=In3-Inabc-k(2二)则（4）变为(5)W=Inu(x，y，z，1)-：X：YZ；利用已知数据（xi，yi，zi，ui），i=1,2，.，n,可求得（X，Y，Z，W），i=1,2,.,n的数据，用三元回归分析方法求出网P、K名的估值如下,AAAA.”w-gx+By+/Z)(6).in

26、ininin其中W二'Wk,X=_'Xk,Y=-Yk,Z=-Zknk4nk鼻nk鼻nk、隹?满足方程liia+1-28+113y=I-。«I21s+1226+123'=I20(:7)l3ia+|32+l33'=I30nnn110c(Xk-X)(Wk-W),120c(Yk-Y)(Wk-W),130c(Zk-Z)(Wk-W)k4kJknnn一2一2一2lii=£(Xk-X),l22=£(Yk-Y)33=£(Zk-Z)(8)kikikdnn112c(Xk-X)(Yk-Y),li3-"-Y)(Zk-Z),kzikWnl2

27、3=(Yk-Y)(Zk-Z),l12=l21,l31=l13,l32=123k1由ct、6、/可求得g2>b2、c2的估值，即利用a=-工,G=-4=-口,AAAabc又由于2M-k=w+lnabcln3(9)(2二)由(5)式可得1,再把a、bc带入(9)式得，2Mk=6+lnabclnr(10)(2二)32222至此得到参数a2、b2、c2、k2的估计值a、b、c、k,把他们带入(2)式分别代替a2、b2、c2、k2,则得不含未知数的解u(x,y,z,t)得近似表达式。6.3模型求解通过第一问求解，我们知道该城市主要受Zn、Cd、Hg,这里主要讨论Zn元素的某一污染区域，讨论用反演法

28、求污染源。从第一问中图九发现Zn主要污染区域在1(x,y)|x12000,15000,y8500,11000)。(12855,8945,18)、(13797,9621,18)、(14325,8666,23)、1)筛选后区域内的点有(12641,9560,11)、(14000,8970,14(14207,9980,14)、(14065,10987,25)、(12734,10344,32)2)初步选择(13797,9621,18)为污染源，将这个数据及t=1带入(3)式，这时u是关于x,y,z的方程且有参数a,b,c,k；3)计算得a=425.73,b=472.8,c=16.655,k=23.56

29、;4)换一个假设污染源，重复上述步骤得误差最小。最终得到污染源位置13522.3,9448.0,23.9得出结论见下表表十五:污染较严重重金属污染源表t11污染源坐标ZnCdHg(13522.3,9448.0,23.9)(945.2,4562.9,69.3)21456.4,11433.2,164.3)(2017.2,3112.3,12.3)(15387.1,9190.0,25.0)(18522.8,3563.9,25.1)(2501.1,3095.5,48.2)七、问题四模型优缺点7.1 模型优缺点7.1.1 优点：1）问题一用单因子污染指数法对八种重金属的污染指数进行分析，得出单个元素的污染

30、指数，并用内梅罗综合污染指数法求取了综合污染指数，对不同功能区的污染程度进行了总结，使得评价更客观，然后用MATLABJ线拟合方法画出了空间分布图，能更加直观看见金属分布状况。2）问题2中运用了因子分析算法，考虑到每个重金属元素污染成因并不是相互独立的，运用SPSS22软件建立了每个重金属元素的相关系数矩阵，通过主因子的属性确定重金属元素的污染原因。优点在于将每个重金属元素联系起来而不是单独分析每个元素的产生原因，这样更贴近实际。7.1.2 中的抛物型偏微分方程形式很好，每一个系数都有其代表的实际意义，符合重金属污染物在土壤中的迁移规律。7.1.3 缺点1）城市土壤污染分析复杂，本文仅对表层土

31、壤样品进行分析，其评价结果不能十分客观和完全真实地反映城市的重金属污染情况。2） 土壤环境质量评价应注意其历史变化过程，不仅要对土壤环境质量现状作出评价，也要对土壤环境质量回顾和土壤环境质量预测影响作出评价。本文仅对现状进行了分析。3）问题3中，利用反演的方法寻找污染源点肯定有所误差。4）没有考虑自然地质变化过程及人类活动影响。7.2收集资料地质环境演变1）重金属元素不仅来自人类活动的影响，同时大气干湿沉降是重金属的主要来源，同时我们假设重金属大气干湿沉降速率此后保持不变，在上述假设条件下，我们需要知道各元素的年平均大气干湿沉降速率，分别对重金属大气干湿沉降对表层土壤中重金属含量的累积影响进行

32、分析，利用下列公式计算：Ct=Co+v*t230其中Ct表示t年后表层重金属含量，Co表示表层土壤重金属含量现状，V表示重金属大气干湿沉降速率，t表示预测年限。利用该计算公式可粗略估计t年后表层土壤中重金属含量。2）若给出了城市主要工厂的具体位置或该工厂的污染排放是否达标等信息则可以更为准确地分析该城市的污染原因及污染源的具体位置。3）我们需要知道上述八种重金属元素在深层土壤中的浓度数据，再利用问题一所述模型对各功能区深层土壤的重金属污染状况作物污染评价。八、模型改进1）问题一中除了采用以上三种污染指数分析法外，还可采用潜在生态危害指数法，研究沉积物中重金属对环境的影响，该种方法不仅反映了某一

33、特定环境中的各种污染物的影响，而且也反映了多种污染物的综合影响，并且用定量的方法划分出潜在生态危害的程度。2）主因子分析的步骤：不仅需要对题目给出的所有数据进行主因子分析，最好还应有针对各个功能区的主因子分析。九、参考文献1姜启源，谢金星，叶俊，数学模型（第三版），北京：高等教育出版社，2003附录(1)城市取样点分布及等高线示意图clearjinshu=xlsread('wurannongdu.xls');%fe属浓度pos=xlsread('zuobiao.xls');%标x=pos(:,1);y=pos(:,2);k=pos(:,4);z=jinshu(:

34、,2);xx,yy=meshgrid(0:500:28654,0:500:188449);zz=griddata(x,y,z,xx,yy,'v4');figure;subplot(2,1,1);holdon;mesh(xx,yy,zz);holdon;plot3(x,y,z,'*');axis(0286540184490408);subplot(2,1,2);contour(xx,yy,zz,20);holdon;zx=zz(1:end-1,2:end)-zz(1:end-1,1:end-1);zy=zz(2:end,1:end-1)-zz(1:end-1,1:end-1);quiver(xx(1:end-1,1:end-1),yy(1:end-1,1:end-1),zx,zy);(2)功能区分布图v,c=voronoin(x,y);figure;holdon;axis(028654018449);fori=1:length(c)patch(v(ci,1),v(ci,2),k(i);plot(x(i),y(i);end(3)重金属污染分布示意图(以As为例)clearj