2011全国大学生数学建模竞赛A题一等奖论文.docx

p城市表层土壤重金属污染分析摘要本文通过对城市表层土壤受重金属污染的过程、实际情况和相关数据的分析， 运用多种数学模型对问题进行求解。在求解第 1 问时运用 Matlab 软件编程绘出二维等高线图，建立图形模型。直 观地给出了各金属元素的空间分布， 见正文第 6 页图 1。根据内梅罗污染指数模型：p + p=平均 最大2对数据整合依次得到各功能区的指数值(见正文第 8 页)，再与背景值的指数值作 比较得到生活区、工业区和交通区的污染程度较大，山区和公园绿地区的污染程 度较小的结论，同时在第一问的基础下，生成了有关各元素浓度与功能区之间的 统计直方图，建立统计模型，通过分析找到了重金属污染的主要原因为：工业污 染、交通污染和居民生活污染。污染程度较大的几种元素是： Cr、 Cu、 Pb 和 Zn。求解第三问时，将重金属在土壤中的传播等效为一种物质的紊流扩散，建立 了菲克扩散模型：?c + ? (c . v ) + ? (c . v ) + ? (c . v ) ? (c . w) = D ( ? 2 c + ? 2 c + ? 2 c ) + F?t ?x x ?y y ?z z ?z m ?x2 ?y 2 ?z 2 c通过逐步降维和高斯分布得到反应传播特征函数：c(x, y) = Q e 5 z0 4(v)Ex(y2)2h 几vEx根据分布曲线特征并结合图 1 找到了污染比较严重的几种元素污染源的位置为：Cr (3000， 6000)， Cu (2700， 3500)， Pb (2300， 3500)， Zn (3000， 6000)、Zn (9500， 4500)、 Zn (13500， 9500)。由于建立的菲克扩散模型时忽略了很多外在因素，为了更好地研究城市地质环 境的演变模式，还应收集诸如当地的自然环境(降雨量，温度等)，地质情况，自 然危害(地震，泥石流等)，人类活动因素的叠加，包括对资源的开发，人类对环 境的破坏以及保护等信息。而对流-弥散方程模型对扩散问题的研究是一个考虑了 多方面因素逐步精确化的模型，所以在第四问中我们选择了建立对流-弥散方程模 型：?a C(x, y,t)?ta = v(?C(x, y,t)?x + ?C(x, y,t)?y) + D(?2 C(x, y,t)?x2 + ?2 C(x, y,t)?y2)得到以下结论：空间相关性体现于：同一时刻空间上其他点对某一点的浓度都有 影响， 这种非局域性是的扩散过程加快， 即随着空间分数阶阶数 的减小， 溶质扩 散速度越快。基于此可以更准确的反应出城市地质环境的演变模式。关键词： 内梅罗污染指数模型；菲克扩散模型；对流-弥散方程模型1一、问题重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量 的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海 量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模 式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园 绿地区等，分别记为 1 类区、 2 类区、 5 类区，不同的区域环境受人类活动 影响的程度不同。现对某城市城区土壤地质环境进行调查。 为此， 将所考察的城区划分为间距 1 公里左右的网格子区域，按照每平方公里 1 个采样点对表层土(010 厘米深度) 进行取样、编号，并用 GPS 记录采样点的位置。应用专门仪器测试分析，获得了 每个样本所含的多种化学元素的浓度数据。另一方面，按照 2 公里的间距在那些 远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件 1 列出了采样点的位置、海拔高度及其所属功能区等信息，附件 2 列出 了 8 种主要重金属元素在采样点处的浓度，附件 3 列出了 8 种主要重金属元素的 背景值。现要求你们通过数学建模来完成以下任务：(1) 给出 8 种主要重金属元素在该城区的空间分布，并分析该城区内不同区 域重金属的污染程度。(2) 通过数据分析，说明重金属污染的主要原因。(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式， 还应收集什么信息？有了这些信息，如何建立模型解决问题？二、模型假设2.1 图表模型及内梅罗模型假设1、重金属污染物在在传播时不考虑海拔的影响。故重金属物质的空间分布于 海拔无关。2 、假设各功能区之间是连续的，不存在城区与非城区的过渡。3、在计算内梅罗指数时仅考虑已知的八种元素， 忽略其它元素对结果的影响。4、以题中所给背景值的内梅罗指数为标准值，而忽略内梅罗指数本身所赋予 的评价环境污染程度的标准。2.2、统计模型的假设1 、在采样测算的过程中所得的各数据之间没有相关性。22、假设背景值所提供的各重金属元素的浓度值为标准值、是可供参考比较的 对象3、当某区中某种重金属的平均浓度高出背景值浓度的一倍以上认为是污染严 重的。4、假设所采集的 319 份样品可以涵盖到整个城区，所以统计的局部结果可以 用来描述整个市区的情况。2.3、菲克扩散模型的假设1 、假设重金属污染物在土壤中的传播为流体的扩散。2 、假设重金属污染物向各个方向扩散的速率恒定且相等。3、假设重金属污染物被植物吸收及沉入到地表 10cm 以下的速率恒定。4 、假设在传播过程中所受的影响可看成是一个可综合的单一因子。5 、假设重金属污染物浓度在传播过程中产生了浓度场且该场是连续的。6、污染源距地表的深度是定值 z ，因此不考虑海拔的影响。 07、重金属污染物在土壤中运输时距地表的高度为定值，可忽略竖直方向的扩 散。8、假设单位时间发生污染量为零即 F 0。C2.4、对流-弥散方程模型模型的假设1 、假设渗透速度很小且方向不变。2 、假设重金属物质传播过程中周围土壤及其他物质是均匀的。3 、假设重金属物质的对流和弥散是一种相对稳定的状态，自身可变性小。三、符号说明P: 内梅罗综合环境污染指数；P ：某单一元素的内梅罗指数的平均值；平均P ：某单一元素的内梅罗指数的最大值；最大c ：重金属物质的浓度；t :重金属污染物传播的时间；3v :重金属沿 x 方向扩散的速度；xv :重金属沿 y 方向扩散的速度；yv z ：重金属沿 z 方向扩散的速度；w ：重金属被植物吸收及沉入地表 10cm 以下的速度；0w ：重金属被植物吸收及沉入地表 10cm 以下的速度的理想化定值；: 重金属沿 x、 y、 z 方向扩散的理想化相等速度；vD ：紊流扩散中分子扩散系数；mF :单位时间发生的污染量；CE :污染物的传播系数；E ：污染物沿 x 方向的传播系数；xE ：污染物沿 y 方向的传播系数；yE ：污染物沿 z 方向的传播系数；z ：分数阶数 ;Q :污染源的污染量；z ：污染源的距地表的竖直深度定值；0h ：污染物传播过程中距地表的深度设为定值；c ：沿 x 方向浓度值；1c ：沿 y 方向浓度值；2c ：沿 z 方向浓度值；3x、 y ：二维空间步长；D ：散项系数；T ：时间微元；c 、 A、 s、 n ：常数。04四、问题的分析与模型的建立4.1 分析与建模本题通过科学的方法对某城市城区的土壤地质环境进行了考察，得到一系列 的数据。而题目要求我们以所给数据来表现出该城区各种金属元素的分布和重金 属的污染情况，这让人很容易想到利用统计学的相关知识来对海量数据整合分析， 从而得到所需答案。于是我们根据需要对附件中所给数据进行初步提取，以取样 点的 x 值为横轴、以 y 值为纵轴、以一种重金属的浓度为等高线，运用 Matlab 软 件编程绘出二维等高线图，清晰地反映出重金属元素在城区的分布情况，这也就 是我们在解决第(1)问中的第一小问时所建立的图表模型。在解决该城区内不同 区域重金属污染程度问题上，我们想到评价环境中某一物质指标的因子，于是我们建立了内梅罗环境污染指数模型，即： P = P平均 + P最大 ，通过此模型来分析该城2区内不同区域重金属污染程度，这样第(1)问就求解结束。4.2 分析与建模第(2)问要求我们通过数据分析说明重金属污染的主要原因，显然在此种情 况下运用数学表达式模型来直观反映污染原因是不太容易的，而运用数据统计图 来反映这一抽象问题就可以将问题简化从而达到直观明了的效果，为此，我们建 立了统计模型(见 5.2 模型求解)。4.3 分析与建模重金属在土壤中传播，通常情况下是随着雨水在土壤中扩散，所以我们把它 的传播理解为一种流体的扩散，由于题目附件中给了我们 x、 y 和 z (海拔)的相 关数据，所以我们以重金属沿 x、 y 和 z 三个方向的传播建立数学模型，同时我们 知道重金属污染物在传播的过程中随着时间的推移 ,大部分重金属都会存留在土 壤中，因此它的浓度是随着离污染源距离越远而逐渐减小的，但总体上它是遵循 物质守恒定律的，通过搜集资料发现，此问题就完全等效为了流体扩散的问题， 把重金属理想化为流体，利用流体的相关知识我们建立以下菲克扩散的偏微分方 程模型：?c + ? (c . v ) + ? (c . v ) + ? (c . v ) ? (c . w) = D ( ? 2 c + ? 2 c + ? 2 c ) + F 1?t ?x x ?y y ?z z ?z m ?x2 ?y2 ?z 2 c5其中 c = c(x, y, z, t)4.4、分析与建模重金属污染物在土壤中的传播只是城市地质环境演变中的很小一部分，城市 地质环境的演变中还会受到诸多因素的影响，如当地的自然环境(降雨量，温度 等)，地质情况，自然危害(地震，泥石流等) ，人类活动因素的叠加，包括对资 源的开发，人类对环境的破坏以及保护，这些都会影响地质环境的演变。除此， 科学的研究表明：物质在实际溶质的运移中客观存在着弥散的尺度效应，溶质粒 子呈现反常扩散的现象。反常扩散本质上是非马尔科夫非局域性的运动，须考虑 运动过程中的时间相关性和空间相关性。此时粒子的扩散运动不再是布朗运动， 其平均平方位移是运移时间的非线性函。所以在考虑了这些许的外界因素，为了 更好的研究城市地质环境的演变模式，我们对其建立对流-弥散方程模型：?a C(x, y,t)?ta = -v(?C(x, y,t)?x + ?C(x, y,t)?y) + D(?2 C(x, y,t)?x2 + ?2 C(x, y,t)?y2)其中 0 0; ?a C(x, y,t)?ta为Caputo分数阶导数：?a C(x, y, t) = T(11- a) j 0(t) (t -n)-a ?C(?(x)n(,)y,n) dn(0 a D m 。再由模型假设知：v = v = v = v 。此外，被x y z植物吸收及渗入地表 10cm 以下的速度恒定为 w 。故上述模型表达式可化为：0?c + v(?c + ?c + ?c) w ?c = E ( ? 2 c + ? 2 c + 1 ? 2 c )?t ?x ?y ?z 0 ?z ?x2 ?y2 10 ?z 2当土壤没有受到重金属污染时，它的重金属浓度值应该为背景值浓度，即浓度函 数 C 应该满足初始条件：c(x, y, z,0) = c (x, y, z)0重金属污染物在土壤中向 x,y,z 三个方向传播将其等效成平板三维挡板原理，即 假设有：c(x, y, z,0) = c c (x, t) c (y, t) c (z, t)0 1 2 3则模型方程可以分离为以下三个方程：?c1 + v ?c1 E ? 2 c1 = 0 (1)?t ?x ?x2(2)?c ?c ? 2 c2 + v 2 E 2 = 0?t ?y ?y2?c3 + (v w ) ?c3 E ? 2 c3 = 0?t 0 ?z 10 ?z 2(3)由(1) (2) (3)分别解得： c = c (x, t)1 1c = c (y, t)2 2c = c (z, t)3 3代入假设，由高斯分布得：c(x, y, z, t) = Q e h( 4几 Et )3由于污染源的污染物在其周围是一种从小区域到大区域连续性的污染，故其产生的浓度场可看成是一种随着时间推移的浓度场的叠加，令 s = x2 + y2 + z 2 ，t)w此式可近似的表现出污染物浓度与其空间位置的关系，又由于题中所给样品均是 在地表 0 到 10cm 之间所收集的，所以我们可以忽略竖直方向的扩散因此有：c(x, y) = j+wdc(x, y, z) = Q e一 5 z0 一 4(v)Ex(y2)一w 2h 爪 vEx这样就得到了浓度与 x、 y 之间的函数关系，这也就是重金属污染物在土壤中传播 的特征方程，可以看出它是从高浓度向低浓度、由点源及面的传播的。同时，浓 度越高说明该取样点理论上离污染源越近。由 C 的函数表达式我们可以看到它是 一个正态分布函数，而任何一个正态分布函数总可以化为一个标准的正态分布函 数， 在一个标准的正态分布函数图像中可以看到它有且仅有一个极大值点。 所以 8 种重金属元素必定有 8 个浓度极大值点， 即至少每种元素有一个污染源， 而第 (2) 问的模型求解过程我们知道污染比较严重的元素为： Cr、 Cu、 Pb 和 Zn 四种，因此 这里我们只确定这四种污染物的污染源，分别根据相应元素的相关数据绘出其正 态分布图，再结合模型一的图象我们可以粗略确定出四种元素的污染源位置为：Cr ( 3000， 6000 )， Cu ( 2700， 3500 )， Pb ( 2300， 3500 )，Zn (3000， 6000)、 (9500， 4500)、 (13500， 9500)。当然，用相同的方法我们也可以粗略确定出其它几种污染较轻的元素的污染源的 位置，此处不再给出。5.4、对流弥散方程模型的分析求解由上式有：c(x, y, z) = limjt dc(x, y, z, t) = Q e一 v(22E(s一x)一 y) 一 4(5) 0 4h爪 vs我们运用数值法对二维时间分数阶对流-弥散方程进行求解，记模型四中的两 个 方 程 依 次 为 ( 1 ) 、 ( 2 ) ， 然 后 建 立 差 分 格 式 。 令x = i编 (i = 0,1,2,.);t = nT (n = 0,1,2,.) 其中 编x 和 编y 为二位空间步长。i x对方程(1)有段项进行有限差分：11?C(x, y, t) = + O(编x)?x ( xi , , yj ,tn ) 编x(3)?y (xi , , yj ,tn ) 编y?C(x, y, t) = + O(编y)(4)? 2 C(x, y, t) = + O(编x)2 )?x 2 (xi , , yj ,tn ) (编x)2(5)? 2 C(x, y, t) = + O(编y)2 )?y 2 (xi , , yj ,tn ) (编y)2(6)左端项则由式(2)变形为：?ta (xi ,yj ,tn ) T(1 _ a ) 0 ?n?a C(x, y, t) = 1 j t (t _n)_a ?C(x , y ,t )ijndn(7)12进一步将式(7)右端离散得到：?a C(x, y, t )?ta (xi , yj tn )= 1T(1_a)jttk +1 (tn _ n)_a dn k =0 kT(1 _ a ) T t n+ jtk +1 (t _ n)_a dn (8)= T 1_a (n + 1 _ k )1_a _ (n _ k )1_a )T(2 _ a) Tk =0+ i j k +1 i j kT(2 _ a) TT 1_a C(x ,y , t ) _ C(x ,y , t )上式右端括号内第一项为常数，记为A：并，考虑二位多空均匀介质稳定流条件下，式(1)可写为：?aC(x, y, t) = _v ?C(x, y, t) _ v ?C(x, y, t) + D ? 2 C(x, y, t) + D ? 2 C(x, y, t) (9)?t 2 xx ?x yy ?y xx ?x2 yy ?y 2将式(3)(4)(5)(6)(8)带入方程(9)，得到显式差分格式：i j k +1 i j k + A =T(2 _ a) TT 1_a C(x ,y , t ) _ C(x ,y , t )C(x , y , t ) _ C(x , y , t ) C(x , y , t ) _ C(x , y , t )_ v i=1 j n i j n _ v i j+1 n i j n xx 编x yy 编y(10)+ Dxx i_1 j n (编ix)2j n i+1 j n C(x , y , t ) _ 2C(x , y , t ) + C(x , y , t )C(x , y , t ) _ 2C(x , y , t ) + C(x , y , t )+ Dyy i j _1 n (编iy)2j n i j+1 n 这里着重研究方程的分数阶项的所表征的反常扩散的性质，故而假设渗透速度很 小且方向不变，因此对对流项的处理采用比较简单的应分格式。分数阶对流-弥散方程考虑了溶质运移的时空相关性。从而从理论上解释了实 际中存在的尺度效应。空间相关性体现于：同一时刻空间上其他点对某一点的浓 度都有影响，这种非局域性是的扩散过程加快，即随着空间分数阶阶数 Y 的减小，溶质扩散速度越快。 4所以，本题中，通过污染源释放污染物，在土壤中随着土壤溶剂向外扩散，由对流-弥散原理，在一定程度上可以很精确的测出某一地理位置上某种物质的含量，将其应用到对土壤污染物质的传播特征和污染源的确定上，就可以根据传播 特征和污染源的位置情况，推测出城市地质环境的演变过程，从而可以解决城市 在保护地质环境上的问题。六、模型优缺点的分析优点：我们通过用 MATLAB 画出实际坐标与各种金属元素浓度的等高线图，以 及利用 EXCEL 表格做出各区各种金属浓度的直方图，可以清楚地反映出各浓度的 含量，以及各区之间的比较，从而对分析污染源的位置提供了良好的参考。在求 解第三问时我们引入菲克扩散模型，把重金属在土壤中的传播等效为紊流扩散， 建立了菲克扩散模型，得到一个函数表达式，把抽象的传播特征问题规律化、数 学化，同时为了更精确的研究城市地质环境的演变模式我们在第四问中建立了对 流-弥散方程模型，使结果更精确化。缺点：我们在建立模型寻找污染源的位置时只是考虑相对理想化的模型，并 且运用了图文结合的方式，并没有严格意义上的客观性的结果，同时我们大量运 用了图表，通过图表来反应实际问题，缺少很严密的数学表达式或方程，具有一 定的主观性。参考文献1窦国仁， 紊流力学，人民教育出版社， 1981 年出版。2周璐莹、吴吉春、夏源，二维分数阶对流 -弥散方程的数值解，高校地质学 报，第 15 卷第四期： 569-575 页， 2009 年 12 月。3郭印，重金属污染物在天然河道中迁移分布的近似计算，中国知识网， 2011 年 9 月 10 日。4 周璐莹、吴吉春、夏源，二维分数阶对流-弥散方程的数值解，高校地质学 报，第 15 卷第四期： 569-575 页， 2009 年 12 月。13附录:某种元素的浓度值与 x、 y 生成二维等高线图的程序代码Z=9.84 8.50 4.69 9.62 4.49 1.77 9.39 7.41 3.50 3.50 4.49 18.38 10.74 8.94 6.14 16.5811.45 6.05 7.12 5.51 7.41 6.69 5.51 6.69 6.14 5.72 4.90 2.00 5.62 5.10 8.72 7.41 8.009.81 5.93 5.51 2.34 21.87 6.47 6.56 5.00 11.68 4.69 4.90 14.08 4.29 8.72 4.79 2.3810.53 9.17 7.56 8.23 8.67 4.90