城市表层土壤重金属污染的逆扩散模型.docx

此文档收集于网络，如有侵权，请联系网站删除城市表层土壤重金属污染的逆扩散模型及污染源确定摘要本文主要依据某城市城区土壤表层重金属污染的采样信息，建立了土壤重金属污染的逆扩散模型，并确定了污染源的分布。首先，本文通过查阅相关资料，对土壤重金属污染的相关问题进行了解释。根据所给319个采样点的GPS信息和功能区数据，使用MATLAB得到城区的平面栅格图。再对采样点的坐标数据及8种重金属的浓度信息，通过平面二维插值和曲线拟合，得到了不同重金属污染物的浓度等值线分布图。将所得浓度等值线分布图与城区的平面栅格图进行合成，从而得到每种重金属污染物功能区分布图，通过比较发现，重金属污染物最严重的为工业区，其次是交通区和生活区。然后，由国家环境质量标准评算公式，计算出5个功能区的综合污染指数以及每种重金属元素的综合污染指数，并给出了重金属污染的主要原因。在分析重金属污染物的传播特征时，以污染物的浓度和高程为主要的影响因素。重金属污染物在传播过程主要以扩散为主，从高浓度向低浓度进行。取扩散的逆过程，污染物浓度由低向高转移，从而确定出污染源的位置。同时取采样点的反向高程，结合重金属浓度逆扩散过程，运用非饱和土壤水分运移模型中的BuckinghamDarcy定律和土壤中溶质对流弥散吸附方程，得到重力势与浓度的影响方程，再结合控制变量法得到经验方程，建立污染物的逆扩散模型。通过模型的求解，确定出城区的重金属污染源位置、污染程度及污染源周围采样点的点号。最后，从建立模型所用思想方法的科学性，对模型的优缺点进行了讨论，并给出了模型改进方法，得到了优化的土壤污染物的逆扩散模型。关键词 重金属 二维插值 污染源 浓度 逆扩散模型一.问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（010厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现要求通过数学建模来完成以下任务：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，确定还应收集的信息。通过这些信息，建立模型进一步解决问题。二.问题分析本题是一个城区重金属污染分布、传播和污染源确定的问题。重金属污染，是指由重金属或其化合物造成的环境污染。主要由采矿、废气排放、污水灌溉和使用重金属制品等人为因素所致1。本题的关键是根据采集的数据得到8种重金属元素在该城区的分布，以及不同区域的重金属污染程度，建立模型得出污染物的传播特征，确定污染源的位置。针对问题（1），根据采样点GPS数据，先将采集点分布的区域进行栅格化，得到城区功能区的平面栅格图。然后从不同重金属元素的平面坐标考虑，将采集的数据去除采样点的高程数据，导入MATLAB中进行二维插值拟合，得到平面坐标下的浓度等值线分布图，进而确定8种重金属元素的城区空间分布。再将栅格图与浓度等值线图合成，确定出城区内不同区域重金属的污染程度。针对问题（2），根据国家环境质量标准，需要对每种重金属元素的单因子污染指数和每个采样点元素的多因子污染指数进行考虑。分析采样数据，分别给出单因子的污染指数和多因子综合污染指数，进而比较数据得出主要污染重金属和功能区污染状况。对问题（3），要分析重金属污染物的传播特征，需要确定影响污染物传播的主要因素。通过相关资料，污染物在土壤中的传播主要是扩散作用的影响，即从高浓度向低浓度扩散，同时重力作用过程也参与其中。因此，需要对采样点的重金属浓度和高程进行考虑。污染源为浓度较高的区域。正常情况下浓度从高向低扩散，但寻找污染源需要由低浓度到高浓度进行逆向分析，为此综合考虑高程和浓度因素。以反向高程（用某一固定高程值与现有高程的差值）作为采样点的高程值，从而建立重金属污染物逆扩散模型，确定污染源所在的栅格区域，得到污染源周围采样点的点号。根据建立的模型，对模型的优缺点进行评价。考虑到实际中污染物的传播不仅与高程和浓度有关，还与采样点的土壤物理化学特性和水特性有关等有关，综合考虑各种因素，减小扩散模型误差，对问题（3）中的模型进行扩充完善，从而得到较为客观完整的逆扩散模型。三.模型假设1.假设城区和自然区所有土壤的物理化学特性都相同，自然区采样的背景值为金属浓度的标准值2.所采集的数据均有效和代表性，且均在所要求的限差内；3.采样点污染物的浓度值即为所在区域1平方公里内的污染物浓度值；4.城区1平方公里中存在两个或两个以上的采样点时，同种重金属污染物浓度为区域内所有采样值的平均值。四.符号说明i：采样点标号。j:1，2，,8，依次代表8种重金属元素As,Cd,Cr,Cu,Hg,Ni,Pb,Zn。Pij:采样点i的第j种重金属元素的污染指数。 Cij:采样点i的第j种重金属元素的实测值。Sj:8种重金属元素的背景值。P: 采样点所有单项污染指数的平均值。Pmax：采样点所有单项污染指数的最大值。Pz: 样点的综合污染指数。Cl：运移区域浓度。Ch：被运移区域浓度。C:单次逆扩散的浓度变化值。E(i,j,k):表示8种元素在不同点位的污染程度的三维数组。五.模型建立求解1.问题（1）首先对整个城区的水平坐标按步长1km进行归化处理,将整个区域划分为1km1km的格网，并用不同颜色分别代表采样点所处功能区，得到城区的栅格图，见图1,其中x30,yCuZnCrPbCdNiAs。查阅资料，8种重金属产生的原因：As：制造农药、防腐剂、染料和医药等，工业生产使用较多；Cd：工业生产，有电镀、采矿、冶炼、燃料、电池和化学工业等排放的废水；Cr：电镀、化妆品原料、皮革制剂、工业颜料以及鞣革、橡胶和陶瓷原料；Cu：输电、工业生产如；Hg：工业生产，如仪表厂、实验电解、贵金属冶炼化妆品等；Pb：电池，含铅汽油，工业生产油漆、蓄电池、冶炼，机械等；Ni：电镀、有磁性物质生产制造；Zn：工业原料、电池。3.问题（3）在研究重金属污染物在土壤中运移特性时，重金属污染物应满足扩散定律与质量守恒定律。根据扩散定律，即扩散物质在单位时间内沿法线方向n流过单位面积的曲面的质量与物质浓度Cx,y,z,t沿法线方向n的方向导数Cn成正比。由此可得扩散物质在单位时间内沿着法线方向n流过的面积为dS的曲面的质量dm为：dm=-D(x,y,z)CndSdt 其中D(x,y,z)为扩散系数，其值与土壤的含水率、土壤密度、物理化学特性、重力、土水势等有关。规定重金属污染物的总质量保持不变，即在运移过程中不产生新的污染物，也不会有污染物被清除。单位时间通过扩散曲面的流量区域内重金属污染物增量相等：sdmnd=vdivdmdv 综上所述得到扩散方程:Cn=CxDCx+CyDCy+CzDCz 在实际应用时，由于已知条件的限制，为了得到扩散系数D(x,y,z)，本文进行了一定的假设与经验数据，并且根据运移结果推算污染源。由此可借鉴简化模型得到推算土壤污染源的数学模型。重金属污染物实际中是在三维空间运移，而题目中没有体现垂直浓度差异，所以本文将重金属污染物运移简化到二维网格模型中，高度将作为计算重力势的属性条件，同时由于题目中没有给出土壤属性，以及考虑到建模可行性与实际性，本文将只考虑重力影响和浓度影响。重金属污染物在土壤中从高浓度区域向低浓度区域运移，在重力作用下像低洼地势运移。重金属污染物胶体由于收到地球重力影响，产生重力势，将迫使重金属污染物向地势较低的方向汇集。从热力学上看，溶质浓度较高的区域，溶质的吉布斯自由能较高，由于系统是始终向低能量方向运动或转化，最终达到能力平衡的状态，溶液中溶质就得从高浓度向低溶液扩散3。从动力学上看，在溶质浓度高的地方，水分子的间隙里含有大量的溶质分子或离子，这样很难再继续溶入，而溶质浓度低的地方间隙较大，对溶质分子或离子有吸引的作用，溶液中溶质就得从高浓度向低溶液扩散。、逆扩散模型对比非饱和土壤水分运移模型中的BuckinghamDarcy定律与土壤中溶质对流弥散吸附方程，为了得到重力势与浓度影响方程，可以使用控制变量法得到经验方程。因为要根据扩散的结果确定污染源的位置，在这里需要考虑扩散的逆过程，即逆扩散模型。、高程相同，浓度不同计算浓度运移经验公式。根据菲克定律，本题假定一理想扩散模型，通过检测中心点浓度变化拟合图像，其图像是一类抛物线曲线，先期传播较快，随着时间积累传播速率减慢。曲线无限接近于x轴但不接近，符合客观实际。根据曲线所反映的情况，建立逆扩散经验公式，用以确定污染源：C1=-Cl1-e-CCh （6）其中，Cl、Ch为运移与被运移区域浓度，较低的为Cl，C=Ch-Cl。考虑低浓度的重金属污染物，在推算污染源时为一极小量，为了简化算法循环，优化算法结构，将浓度低于正常范围最大值的区域予以舍弃，即为零。实际中，结合野外采集的重金属分布及浓度的数据不难发现，自然界土壤中本身存在一定量的重金属物质，寻找污染源时，这类低浓度区域不用列入推算范围。计算模拟中也发现，低浓度地区的重金属演变速率明显慢于高浓度地区。由于采样点是离散的，为了提高精度，对经计算得到的新重金属污染物浓度地域分布矩阵进行二维内插。高程分布是连续的，也要进行二维内插。由于山地区高程变化较为剧烈，在对正常的高程进行作图时会发现大片的非山地采样区域的高程变化不能得到很好的体现。在为了方便分析污染传播特性和理解推算逆扩散的算法，需要对采样点高程进行处理得到反向高程，取高于采样点最大值高程的某一固定高程值与采样点高程的差值作为采样点的反向高程。如图4所示。图中粗线条显示的即为实际中低洼的地区。图4 反向等高线在进行二维内插后浓度进矩阵的尺寸被拉大（即分辨率增加），所以对内插后的矩阵进行处理，使内插后的矩阵与原矩阵尺寸相同，方便比对和迭代计算。依次选取浓度非零点进行循环计算，将点周围8个点浓度进行比对，带入公式（6），将四周转移浓度值计入单次逆扩散矩阵。考虑到对角线上的4各点对于网格模型来说传播时间较长，所以要在上式的基础上将结果除以2(在图5B)，即得到所需的单位时间的浓度变化量。由于重金属污染物的运移具有整体性和同时型，考虑到算法计算的依次性，与实际情况不符，故将运移的重金属污染物存入单次逆扩散矩阵，一次循环后整体运移。 A B图5 金属砷逆扩散示意图6.对25步进行多次循环计算可得到高程相同浓度不同的情况下逆扩散模型。、高差影响因素在考虑高差因素的影响时，既要考虑需要分析垂直渗透因素，也要考虑平面逆扩散因素，但所给出的数据量不足以进行以上分析。为了简化模型的计算，对高程的单个因子的逆扩散推算可以类比于单因子（浓度差）的逆扩散模型，同时从高差对浓度影响的特性考虑，最后得到关于高程的逆扩散方程：C2=-Chl1-e-h （7）其中，Chl为相对低洼区的浓度, Chh为相对高处的浓度，h=Chh-Chl。为常系数，调整他能够得到不同的逆扩散推算的结果。最后可求得单次、单向的逆扩散的浓度：C(i,j)=C1+C2 （8）在其后需要进行逐点的循环运算，可得到单次逆扩散的浓度变化矩阵。图6 单次逆扩散矩阵将逆扩散矩阵加到原浓度矩阵上即完成了一次推算，要确定污染源则需要将对上述推算进行迭代。本文依次对所给出的八种重金属进行迭代运算，大多收敛特性较好，其中元素铬的推算所需的迭代次数较多，在这里对其进行演示(如图7)，其他元素的推算结果见附录3。 A1次 B100次 C1000次 D3000次图7 金属Cr逆扩散迭代计算结果分布图最后采用三维数组E(i,j,k)对八种元素的污染程度进行综合分析, 其中一组i、j、k可以表示一种元素在污染源的二维坐标及其污染程度。进一步需要知道各元素污染程度的分级标准以及他们的权值，在这里参照表7中的权值以及分级标准4。元素权重分级标准(g/g)轻度污染中度污染重度污染As0.21015.415.1215.1223.52.5Cd0.13000.190.380.380.540.54Cr0.074249115.4115.4300300Cu0.037020.443.7143.71100100Hg0.30390.0510.0260.0261.001.00Ni0.100019.969.569.65175.0175Pb0.13004370.070.0350.0350Zn0.014997115.30115.3300.0300表7 八种元素的污染分级标准及其权重值对八种元素的污染情况进行综合加权计算，即得到综合污染源的位置，图8是用1公里为单位的网格表示的该城区的污染污染源位置及其污染程度。 图8污染源分布根据采样信息可以得到污染源附近的采样点号。轻度污染源中度污染源重度污染源103 135 250 25925 27 42 43 49 53 65 81 153 154 221 254 2637 8 13 1719 29 30表8 污染源附近的采样点号4.问题（4）本文结合污染物扩散特性，建立了重金属污染的逆扩散模型，主要优点：算法简单，计算快捷，避免了复杂的积分、梯度计算，在数据量庞大时，对数据处理依然较快精度较好，结果符合实际客观规律，基本上确定了污染源。考虑到样本点采集的随机性，分布不均匀，通过适当插值以及综合评价，减小了样本误差产生的影响。模拟了重金属污染物反扩散过程，能够直观的查看重金属污染物分布随时间的变化情况。模型的不足：将三维空间重金属污染物运移活动简化为二维平面考虑，忽略了重金属污染物的垂直运移。本文规定土壤为无运动的稳态环境，忽略了土壤环境变化对重金属污染物运移的影响。没有考虑重金属污染物挥发、富集、移除或增加。对于问题（3）所建立的重金属污染物浓度的逆扩散模型，只考虑了土壤中单一浓度扩散对污染物浓度的主要影响作用，而水特征曲线是研究土壤中污染物扩散的重要手段，主要是通过对对流、扩散的影响来影响重金属污染物的扩散。还应收集该城区的土壤水特征曲线。常用的土壤水特征曲线经验公式较多，本文选用Van Genuchten提出的经验公式：-rs-r=11+hn1-1n 式中，为土壤进气吸力的相关系数，n是与形状相关的系数，为土壤含水率，S为饱和土壤的含水率，r为土壤剩余含水率，h表示土壤基质吸力。当=S时，即土壤饱和状态时，土壤基质吸力为零，而土壤基质吸力大小将影响重金属污染物的运移能力，即确定扩散系数D5。此外，采样点土壤本身的物理化学特性对污染物扩散转移也具有较大影响，包括土壤的物理化学特性，如pH值，土壤密度，通透性，土壤含水率，土壤的紧实程度等。对所建立的重金属污染浓度逆扩散模型中，对不同的采样点加入物理化学特性内容以及水特征曲线，确定相关的逆扩散浓度矩阵，确定综合的污染物浓度扩散模型，研究该城市地质环境的演变模式。参考文献1 互动在线（北京）科技有限公司，互动百科重金属污染/wiki/%E9%87%8D%E9%87%91%E5%B1%9E%E6%B1%A1%E6%9F%93，2011.9.10。2 环境保护部环境标准研究所，中国环境标准网，，2011.9.9。3 雷志栋，杨诗秀等，土壤水动力学M，北京：清华大学出版社，1988.4 李雪梅，邓小文等，污染因子权重及区域环境质量综合评价分级标准的确定，干旱区资源与环境，1003-7578，20105 Van Genuchten M Th. A closed form equation for predicting the hydraulic conductivity of unsaturated soil J，U.S:Soil Science Society of America Journal,1980,33:892898.6 马莉，数学实验与建模，北京：清华大学出版社，2010。7 姜启源，谢金星等，数学模型，北京：高等教育出版社，2003。8 李学垣等，土壤化学M，北京：高等教育出版社,2005。附录1. Cd,Cr,Cu,Hg,Ni,Pb,Zn各元素的在浓度等值线图2. Cd,Cr,Cu,Hg,Ni,Pb,Zn各元素的在浓度区域分布图3. 逆扩散推算结果4. 逆扩散模型主程序clearclcp=load(D:MATLAB_7.0work建模2011高教社杯2011Q1.txt);E=load(D:MATLAB_7.0work建模2011高教社杯2011Q2.txt);B=load(D:MATLAB_7.0work建模2011高教社杯2011Q3.txt);P=zeros(319,4);for i=1:319 m=p(i,1)/1000; n=p(i,2)/1000; P(i,1)=floor(m)+1; P(i,2)=floor(n)+1; P(i,3)=p(i,3); P(i,4)=p(i,4);endfor p=1:8 E1=E(:,p); for i=1:318 for j=i:318 if P(i,1)=P(j+1,1)&P(i,2)=P(j+1,2) E1(j+1,1)=(E1(i,1)+E1(j+1,1)/2; else continue end end High(P(i,1),P(i,2)=P(i,3); Zone(P(i,1),P(i,2)=P(i,4); Elem_1(P(i,1),P(i,2)=E1(i,1); end A_High=interp2(High,1); A_Zone=interp2(Zone,1); A_Elem_1=interp2(Elem_1,1)-B(p,1)-2*B(p,2); for i=1:57 for j=1:37 if A_Elem_1(i,j)0 A_Elem_1(i,j)=0; end end end for i=1:28 for j=1:18 G_Elem_1(i,j)=(A_Elem_1(i*2-1,j*2-1)+A_Elem_1(i*2,j*2-1)+A_Elem_1(i*2-1,j*2)+A_Elem_1(i*2,j*2)/4; G_High(i,j)=(A_High(i*2-1,j*2-1)+A_High(i*2,j*2-1)+A_High(i*2-1,j*2)+A_High(i*2,j*2)/4; end end G_Elem_1(i+1,j+1)=A_Elem_1(i*2+1,j*2+1); G_High(i+1,j+1)=A_High(i*2+1,j*2+1); for k=1:20 C_kuo=zeros(29,19); for i=1:27 for j=1:17 if G_Elem_1(i+1,j+1)=0 continue else C_z=0; for m=0:2 for n=0:2 if m=1&n=1 continue elseif G_Elem_1(i+m,j+n)=0 continue elseif G_Elem_1(i+m,j+n)G_Elem_1(i+1,j+1) G_Elem_1s=G_Elem_1(i+m,j+n); G_Elem_1B=G_Elem_1(i+1,j+1); else G_Elem_1s=-G_Elem_1(i+1,j+1); G_Elem_1B=-G_Elem_1(i+m,j+n); end if G_High(i+m,j+n)G_High(i+1,j+1) G_High_s=G_High(i+m,j+n); G_High_B=G