重金属污染源的数学建模

《重金属污染源的数学建模》由会员分享，可在线阅读，更多相关《重金属污染源的数学建模（18页珍藏版）》请在装配图网上搜索。

1、2011高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包扌占电话、电子邮 件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问 题。我们知道，抄袭别人的成果是违反竞赛规则的，如果引用别人的成果或其他 公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正 文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反 竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）:A我们的参赛报名号为（如果赛区设

2、置报名号的话）：20001032所屈学校（请填写完整的全名）：中国人民解放军国防科学技术大学参赛队员（打卬并签名）：1李腾烦2. 袁雪强3. 赵泉朴指导教师或指导教师组负责人（打印并签名）：杨文强FI期：2011 年 9 月 9_|赛区评阅编号（由赛区组委会评阅前进行编号）:2011高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）:赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）:基于MCMC的城市表层土壤重金属污染分析摘要：本文针对某城市城区表层土壤重金属污染分析的问题，

3、从各分区重金属元 素污染的空间分布、污染程度、污染原因及传播途径等多方而进行了全面的统计 与研究，同时采用贝叶斯推理和马尔科夫链蒙特卡洛法(MCMC),利用MATLAB 进行编程，对所设计的模型进行求解，得到了各重金属污染源的地理位置坐标， 最后对城市地质环境演变模式的问题进行了探究，得出了初步的结论。首先，根据所给数据信息，利用MATLAB及EXCEL作图，得到了该地区的地势 图及各重金属污染程度的分布图。对于问题1,通过统计各重金属的重污染点在各区域的分布和计算不同区域各中 金属污染程度之间的相关系数，得到了各重金屈的大致空间分布，并采用改进的 内梅罗指数法，得出不同区域重金属的污染程度。

4、对于问题2,在问题1所得数据及结论的基础上，结合实际情况，对各重金屈污 染的主要原因进行了较为全面的分析。对于问题3,通过对污染的传播特征的研究，基于城市的地质特性，得到空气传 播为污染主要传播方式的结论，并由此建立关于气体传播的概率模型，采用贝叶 斯推理和马尔科夫链蒙特卡洛法，利用MATLAB进行50000次抽样，得到收敛 的坐标值，从而计算出各重金属污染源的地理位置坐标分别为：As(2742,7294)、 Cd (2382,3693)、Cr (3470,2309)、Cu (2707,2295)、Hg(2708,2294). Ni (1647,2729)、Pb (2882,3617)、Zn

5、(4152,2299)。对于问题4,以题中城市土壤重金属污染的演变规律为例进行了一定深度的探究, 并对所需信息及模型设计方案进行了初步的讨论。关键词：内梅罗指数法相关系数贝叶斯推理MCMC法MH算法一、问题的提岀与分析(一) 背景介绍随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影 响口显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数 据资料开展城市环境质星评价，研究人类活动影响下城市地质环境的演变模式， 口益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地 区等，分别记为1类区、2类区5类区，不同的区域环境受人

6、类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1 公里左右的网格子区域，按照每平方公里1个采样点对表层土(0-10厘米深度) 进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得 了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那 些远离人群及工业活动的门然区取样，将其作为该城区表层土壤中元素的背景值。（二）要解决的问题现给出采样点的位置、海拔高度及其所属功能区等信息，8种主要重金属元素在 采样点处的浓度及8种主要重金属元素的背景值，解答以下问题：1. 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区

7、域重 金属的污染程度。2. 通过数据分析，说明重金属污染的主要原因。3. 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。4. 分析所建立模型的优缺点，并说明为更好地研丸城市地质环境的演变模式， 还应收集什么信息，有了这些信息，如何建立模型解决问题。（三）问题分析整个问题可以分为三类分别进行求解：问题1和问题2可将其归结为对所给数据的统计与分析，可以通过作图或列表等直观的形式针对 不同金屈和区域进行相关统计，再对统计的结果加以比对与分析，即可得出结论。 问题3确定污染源的位置实质就是要找到使周围污染达到最大的点，这样可以转化为概 率模型，在统计数据的基础上，运用贝叶斯推理马尔科夫链

8、蒙特卡洛法可加以计 算和确定。问题4开放类问题，可根据现有资料进行回答。二、模型假设1题目中所给出的样本均具有普遍性及代表性，而非所测地区中的特殊情况。2. 该地区的空气流动没有显著的方向性特征。3. 污染物浓度在収样前后较长一段时间内不会有显著的变化。三、符号说明P :土壤污染元素综合污染指数C,： 土壤中污染元素i的实测值土壤中污染元素i的评价标准值s：任一闭合曲面m3：s内指定所围的区域M1 :通过S流入Q的质星m2：d内减少的质暈m3：q内由F浓度的变化引起的质量增加kl衰减系数8:贝叶斯公式模型参数3加:模型参数的个数y贝叶斯公式观测数据T（X兀）：马尔可夫蒙特链卡洛法中的分布接受概

9、率乃（文）：目标概率密度函数S：距离污染源的距离p:沉降系数h:海拔高度%:临界海拔1：重金属元素吸附系数四、模型建立与求解（-）空间分布的特征根据所给数据，对采样点按区域进行重新的排用，由此作出各采样点的地势分布 图（图一），乂由各种金屈在各区域的浓度值，作出金属污染程度在各区域的分 布图（以Pb （铅）为例，见图二，其它金屈元素分布图见附件），由此分析可 得到各金属元素在该城区的空间分布。图一采样点地势分布图图二Pb污染程度分布图以Pb为例进行分析，经统计,Pb的重度污染区域主要集中于图的左下方， 且Pb的重度污染点在工业区占62.8%,另有20.9%重度污染点则分布在左下方 的交通区中，

10、由此可得Pb的污染区域主要分布于工业区和左下方的部分交通区, 对于其它金属元素，可以用相同的方法，还可以通过计算不同区域各金属元素污 染程度的相关系数，找出各金属元素污染之间的联系，这样可以有效提高分析效 率和准确性。例如对表一中工业区的各金属元素污染程度相关系数（其它区域见 附件）的比较可以得出，Cr与Cu、Cr与Hg、Cu与Hg、Cd与Pb、Cd与Zn、 Pb与Zn相关系数较大，可以认为它们具有相同或相近的空间分布。表一工业区各金属元素污染程度相关系数CdCrCuHgNiPbZnAs0.3286010.3796280.1529320.1812950.6897250.3946630.5177

11、31Cd0.5409850.5664580.5331710.4886660.8292160.753606Cr0.9197260.9021390.698250.6753990.69507Cu0.983480.502780.669660.621675Hg0.4791390.6124650.590405Ni0.577620.634098Pb-0.738756棍据以上数据分析，可以得到齐甫金属元素大致的空间分布，分别为: l.Hg、Cd、Pb、Zn、Cr均主要分布于工业区和左下方的部分交通区；2. As Cu主要分布于工业区；3. Ni的分布比较分散，没有明显的富集现象。（二）重金属污染程度的划分为了

12、全面、综合地反映污染物的整体污染水平，需要一种同时考虑污染物综 合污染水平的多因子评价方法，这里采用内梅罗指数法（内梅罗抬数评价法的 修正及其应用），其计算公式如下：式中，P为土壤污染元素综合污染指数：C,为土壤中污染元素，的实测值：为 土壤中污染元素i的评价标准值。内梅罗指数土壤污染评价标准如表二。表二内梅罗指数土壌污染评价标准等级内梅罗污染指数污染等级IP0.8清洁（安全）II0.7P2.50尚清洁（警戒线）III2.50 P4.25轻度污染IV4.25 P7.20重污染但由丁内梅罗指数法过分突出极大值对水质污染的彫响，会出现评价项目中 只有一项指标值偏高，而其他指标值均较低也会使综合评分

13、值偏高的情况，故对 内梅罗法进行改进如下：其中Cmaxn工qi=lill式中：q为权重值前项的平均值,n值根据评价数据确定：C,为前n项的值；m为前门项中h = %的项数。对于以上8种重金属，取权重值如表三:表三各重金属权重值AsCdCrCuHgNiPbZn权重值0.0700.1400.0140.0200.6980.0300.0700.020计算得到域内不同重金属的污染程度评价结果如表四:表四各区域内重金属污染程度AsCdCrCuHgNiPbZn生活区IIIIIIIIIIIIIIIIIIIIVTsIPIXIIIIIIIIIVVIIIIIIV山区IIIIII11IIIII11II交通区IIIII

14、IIVVVIIIIIIIV公园绿地区IIIIIIIIIIIVIIIIIIV（三）重金属污染原因分析根据（一）、（二）屮所得数据及结论，将各重金属污染的原因做如下分析： 1对于除Ni以外的其他重金属元素，均有60%以上的重污染点分布在工业区， 由此可认为重金属污染主要源白工业区，这主要是由于工业区中的工厂在生产过 程中产生了含该种金属的废弃物，在排放时也没有进行必要的净化处理而直接排 入水体或大气中，从而造成了污染的形成与扩散，2. Hg、Cd、Pb、Zn四种金属元素的均有约20%重污染点分布在交通区，这些污 染主要来1于公路上行驶车辆燃烧汽汕或柴汕所产生的尾气沉降及轮胎的磨损。 3金属元素Zn

15、在生活区和公园绿地区也有一定的分布，这可能是由于该区域所 使用的肥料中含有Zn元素。4. 金属元素Ni的污染点较为分散，可认为该污染与人为因素关系不大，主要來 自于土壤母质本身。（四）重金属污染物的传播特征重金属一般不会直接通过土壤进行传播，它一般通过如下方式进行传播：1. 空气传播工业区中的工厂会排放出大量的工业废气，这些废气中可能含有大量的重金 属颗粒，这些颗粒会随空气的流动而分散到城市各处，从而造成污染，另外，废 气也可能会随降雨渗入土壤。2. 水源传播工厂排出的工业废水中含有大量的重金属离子及沉淀物，这些物质均可随水 渗入土壤，造成污染。3. 工业废渣和生活垃圾废渣及生活垃圾的堆积会造

16、成重金属元素的大量积聚。另外，重金属的传播受地形因素的影响较为明显，从数据来看，该城市整体 海拔较低，地势起伏不大。相对1000米的取样间隔，300米的高度差不会因高 空吸附或雨水冲刷造成污染物沉积异常。但山势对风向的影响势必导致污染物浓 度区域性变化。图三城区地势图如上图的标注，山脚、较易形成污染物富集；由于I、II山谷的 存在，经过山谷的大风在A、B区域减速后极易产生污染物富集：此外，迎向污 染源的山坡受到的污染会明显大于背向污染源的山坡。从Hg、Cd的污染分布情 况可以验证这儿点（如图四）：图四Cd污染程度分布图(五) 建立模型确定污染源在进行污染源的确定之前先进行筛选，如果污染是由点污

17、染源引起的，那么 污染源就一定处于污染程度最大，即重金属浓度最高的采样点的附近。由于城市中建筑物较多，人类活动频繁，导致局部气流复杂，因此通过空气 传播占主要地位，乂由于城市缺少地表水系，雨水及其它污水会通过下水管排出 城外，故影响相对较小，因此只讨论通过空气传播的模型：设u(x,y,z)是t时刻点(x,y)处某种重金属污染物的浓度。任取一个闭合曲而S,它所围的区域是Q,由于气体扩散，由t到At时刻的这段吋间内，通过 S流入Q的质量为：由高斯公式得:仆 T fff / r d, 52uM产 J页+ b 话+ c 左)dxdydzdt其中a2fb2fc2分别是沿xyz方向的扩散系数。由丁衰减(例

18、如吸收、代谢等)，Q内的质量减少为:t+AtM3 = | JJJk2udxdydzdtt Q其中k?为衰减系数。由于重金属污染物在土壤中很难降解，因此衰减系数可以取为0。由质量守怛定律，在Q内由于扩散和衰减的合作用，枳# J-0内物质的总量为Mi-M2从另外的一个角度來看，Q内由于浓度的变化引起的质量增加为:M3 = jjju( x, y, z, t + At) - u( x, y, z, t) dxdydzQ显然M3 = MMy即52udxdydzdt由的任意性得:站丟+b-芳+內k-u = a设扩散源在点(,yo)处，则此扩散问题满足以下方程组：由假设可知，如果认为经过相当长的时间后，扩散

19、已经停止，则物质处于平衡状 态,有:dll a=0即：芈+b?芈+ C 次 x dy由丁污染物在空中传播时，地而上污染物的沉积暈与z方向上的传播无关，故上 述方程可以改写成：7 dd t r Aa - + b ku = 0d2x d2y因此问题转化成二维含源反演问题，可以用贝叶斯推理來估计参数a?和彷2。贝叶斯推理的基础是贝叶斯定理，它可以表述如下：p(e)p(yl&)p(y)P(e)P(yl&)其中,0为模型参数，y为观测数据，p(&)为参数的先验概率密度数，p(y|&)为 似然函数，p(&|y)为参数的后验概率密度函数。在贝叶斯推理中，p(&)表示在 未获得测量数据之前，对模型参数分布的认

20、识，主要来源于以往数据、经验和主 观判断等p(y|&)代表模型参数拟合测量数据的程度，越大表示拟和效果越好， 反之则差。p(O|y)表示了获得测皐数据以后模型参数的分布规律，即为在统计 反演意义下的反问题的解。 其步骤通常分为以下三步：1. 确定先验概率密度数先验概率为均匀分布：p(0n假设参数之间相互独立则:一般情况下，对于参数仿，如果仅知道它分布的区间珂飞,那么可以假设bx-ax0其他mp(o)= f!p(q)其中加为模型参数的个数。2. 确定似然函数似然函数的构造对于贝叶斯推理的结果具有很大影响，一般可以认为测量误 差服从正态分布sN（“ 亍），此时似然函数可以表示为：p(y|) =re

21、xP(2心P（f（8）y）T（f（&）y）Icr其中n为测暈个数，f（e）可以用解析解表示，一般情况下，f为数值解。3. 对后验概率密度函数抽样，进而获得参数的估计值选取马尔科夫链蒙特卡洛法（MarkovchainMont Carlo,简记为MCMCJo算法 为MH算法，表示如下：设定参数初始值马e 0,12000并设定k = l；（2）proposal分布丁（兀）中生成召6（0,12000,计算接受概率：a(x,) = inin 1,（幼T （兀I文）龙（海）T（刃忑）生成0,1之间均匀分布的随机数仁（4）如果设置+1 = X ；否则+1 = o当k = K（迭代次数）时，停止；否则k =

22、k + l,继续第（2）步。其中，龙（刃为目标概率密度函数，在贝叶斯推理中为似然函数。（六）模型的求解通过MATLAB编程，使用MH算法进行50000次抽样，使H标坐标值趋于 半稳，最终得到各重金属元素污染源的地理位置坐标如表五。（图五为Hg元素 进行抽样后的迭代曲线图，可以看到图线最终趋于一条半稳的直线。）图五Hg元素抽样后的迭代曲线图表五各重金属元素污染源地理位置坐标AsCdCrCuHgNiPbZnxfil27422382347027072708164728824152Y值72943693230922952294272936172299(七) 对城市地质环境的演变模式的探究地质环境的演变模

23、式是一个地域在一段较长时期内所表现出來的地质变化 规律与趋势，因此，要对城市地质环境的演变模式进行研究，就必须举握所研究 城市在长期的地质变化过程中形成的规律，以研究题中城市土壤重金属污染的演 变规律为例，就需要收集该城市每年生活、匸业等重要污染源的废弃物排放量， 以及每年的生物降解量，另外，该城市的水文特征及气候特征也会对重金属元素 扩散产生重要的影响.因此还应该考虑降雨量及气候条件的信息。拥有这些数据 后，可以建立因子分析法，回归分析，曲线拟合等模型进行分析，进一步还可使 用马尔科夫链等方法对该城市未来的污染状况进行预测与评估。五、模型评价优点：1. 通过改进后的内梅罗指数法引入了权值，避

24、免了极大值对结果的较大影响, 提高了数据的准确性和可信度。2. 通过马尔科夫链蒙特卡洛法(MCMC)简化了微分方程的求解，降低了计 算时的复杂度。3. 通过对地形的综合分析，剔除了地形影响较大的区域，简化了模型。 缺点：1. 没有兼顾降雨、降解等其他对污染产生作用的因素的影响。2. 对污染的传播规律还有进一步探索的余地。六、参考文献1 朱嵩、刘国华、毛根海、程伟平、黄跃飞，利用贝叶斯推理估计二维含源 对流扩散方程参数，四川大学学报，第40卷第2期，2008年5月。2 李尚志、程继新，数学建模教程，2011年3 杨忠平、卢文喜、辛欣，长春市城区表层土壤重金属污染來源解析， 2011 年 9 月 104 赵彦峰、郭恒亮、孙志英，基于土壤学知识的主成分分析判断土壤重金属 来源，地理科学，第28卷第1期，2008年2月5 李亚松，张兆吉，费宇虹，王昭，内梅罗指数评价法的修正及其应用， 水资源保护，第25卷第并期，2009年11月