1、数学建模A题 城市表层土壤重金属污染城市表层土壤重金属污染分析摘要城市化进程的加快伴随着城市环境的污染加重，重金属污染便是一个重要的原因。针对不同的问题建立了不同的数学模型，结合所给的海量数据，求解模型，对土壤的重金属污染进行分析。针对问题一，首先根据所给的319个样本点的数据，采用双线性插值算法建立了空间分布模型，分别作出了8种重金属在该城区的空间分布图(见附录3)，清晰地表现出8种重金属元素在该城区的浓度分布；然后利用潜在生态危害指数法建立了重金属元素污染程度模型，最后得到8种重金属元素分别在5个区的污染程度以及综合污染程度，具体分析结果为：生活区为中等污染，工业区为很强污染，山区为轻微污
2、染，主干道区为很强污染，公园绿地区为中等污染。针对问题二，建立了主成分分析模型，借助SPSS软件对数据进行处理，得到影响5个区的主要因子。影响生活区的主要重金属元素为Pb,Zn,Ni，影响工业区的主要重金属元素为Cr，As，影响山区的主要重金属元素为Zn,Cd,Hg，影响主干道区的主要重金属元素为Cu，Hg，As，影响公园绿地区的主要重金属元素为Cd，Hg，然后对各个元素进行了原因分析。针对问题三，通过分析重金属元素的传播特征，确定主要传播特征为扩散，并据此建立了一维有源扩散模型，利用最小二乘法对模型进行求解，分别确定了8种重金属元素在该城区的污染源位置坐标，共有31个污染源(具体数据见文中表
3、5.3.1)，进一步对污染源进行精确定位，发现实际上存在一个污染源向外扩散多种重金属元素，最终确定污染源的实际总数量为18个。针对问题四，首先对前三问的模型进行了优缺点的评价；然后在扩散的基础上，结合了以下3个影响重金属元素传播的因素：(1)、重力、环境等引起的扩散系数各向异性；(2)、重金属污染物传播中的对流特征；(3)、污染源强度的可变性。据此建立了三维对流扩散模型，利用分裂浅水波方程的方法，对模型进行简化，得到模型求解方法。关键词：重金属污染；空间分布；双线性插值算法；潜在生态危害指数；主成分 分析；扩散模型一 问题的重述随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量
4、的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土(010 厘米深度)进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些
5、远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。现通过数学建模来完成以下任务：(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度；(2) 通过数据分析，说明重金属污染的主要原因；(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置；(4) 分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？二 问题分析2.1 问题一的
6、思路分析问题一需要根据所提供的数据来建立模型，绘出8种重金属元素在该城区内的空间分布，数据提供的信息包括空间三维坐标及浓度值,可以采用双线性插值算法，建立空间分布模型，运用Matlab软件便可以作出8个重金属元素的空间分布图。在分析污染程度的时候，需要对五个区域分别进行分析，采用潜在生态危害指数法，便能够得到8种重金属元素分别在5个区的污染程度以及综合污染程度。2.2 问题二的思路分析问题二要求分析污染的主要原因，因为题目所给变量较多，可以建立一个主成分分析模型，用较少的变量代替较多的变量，然后用SPSS软件进行分区求解，再结合问题一所得到的8种元素的空间分布图，便能够确定引起污染的主要因素及
7、这些元素污染的主要原因。2.3 问题三的思路分析问题三中需要先分析污染物的传播特征，忽略风速、海拔等因素对传播的影响，依此建立有源扩散模型，并且可以将问题简化为一维扩散，得到的便是一维泊松方程，采用最小二乘法来确定污染源的位置坐标。2.4 问题四的思路分析优缺点评价需要客观地分析模型的优势及缺陷。结合问题三，综合考虑重金属元素的传播特性，如扩散系数的各向异性、对流的影响及污染源强度的不稳定性等，建立出一个三维对流扩散方程，求解方法要求易于实现，分裂浅水波方程的方法可以大大简化求解的难度。三 模型假设1、假设该城区不受周围其他城区的影响；2、假设采样后各个重金属元素的浓度参数不再发生变化，即已经
8、达到稳定状态；3、假设每个污染源重金属元素发生扩散时互不影响；4、假设在问题三中，空气对流、重力等因素对扩散没有影响；5、假设在问题三中，污染源的强度及扩散系数恒定；四 符号说明符号符号所表示的意义C重金属浓度值Cfi第i种重金属的污染系数Ci第i种重金属的浓度实测值Cmi第i种重金属的浓度计算值Cni第i种重金属的浓度背景值Eri第i个重金属元素的潜在生态危害系数Tri第i种重金属的毒性系数D扩散系数Dx、Dy、Dz三个方向上的扩散系数S污染源强度r距离x、y、h三维坐标m1、m2积分产生的待定系数计算值与实测值之间的差u、v、w三个方向上的对流速度t时间五 模型的建立与求解5.1 问题一模
9、型的建立与求解5.1.1 模型的建立(1) 空间分布模型双线性插值算法1，又叫双线性内插，其核心思想是在两个方向分别进行一次线性插值。图5.1.1 双线性插值示意图如图5.1.1所示，若求P(x,y)点的浓度C(x,y)，已知四个临近点Q11(x1,y1)，Q12(x1,y2)，Q21(x2,y1)，Q22(x2,y2)的浓度，对于标准的双线性插值算法，按照下面的过程进行插值。X方向的线性插值： (5.1.1) (5.1.2)Y方向的线性插值： (5.1.3)基于以上的模型，利用Matlab软件，只需要调用griddata函数，即使二维空间数据不是标准的网格，该函数也能自动实现双线性插值，得到
10、精度很高的空间分布，利用该模型也可以得到海拔的三维图。(2) 污染程度模型重金属是潜在危害的重要污染物，因此本文采用潜在生态危害指数法来分析不同区域重金属的污染程度2,3。单个重金属元素的潜在生态危害系数Eri： (5.1.4)综合潜在生态危害指数RI： (5.1.5)式中，Eri为第i种重金属的潜在生态危害系数；Tri为第i种重金属的毒性系数；Cfi为第i种重金属的污染系数；Ci为第i种重金属的浓度实测值；Cni为背景值。5.1.2 模型的求解(1) 空间分布模型的求解按照以上模型，利用Matlab作出该城区的三维海拔图，从图中可以看出，该城区地势不平坦，生活区及主干道路区主要位于低海拔区。
11、图5.1.2 该城区的三维海拔图利用同样的模型，分别导入8种重金属元素的数据，作出如图5.1.35.1.10的浓度空间分布图，图中各个采样点数字代表海拔，颜色代表浓度，大图见附录。 图5.1.3 金属元素As在该城区的空间分布 图5.1.4 金属元素Cd在该城区的空间分布 图5.1.5 金属元素Cr在该城区的空间分布 图5.1.6 金属元素Cu在该城区的空间分布 图5.1.7 金属元素Hg在该城区的空间分布 图5.1.8 金属元素Ni在该城区的空间分布 图5.1.9 金属元素Pb在该城区的空间分布 图5.1.10 金属元素Zn在该城区的空间分布(2) 污染程度模型的求解从参考资料中能够查到8种
12、重金属的毒性系数及污染程度分级表3，如表5.1.1与5.1.2。表5.1.1 重金属的毒性系数重金属As Cd Cr Cu Hg Ni Pb Zn 毒性系数10302540551表5.1.2 潜在生态危害分级表污染程度轻微中等强很强极强Eri320RI1200利用EXCEL将五个区域内的采样点筛选出来，分别用上述模型求解，得到各个区的Eri和RI，如表5.1.3与5.1.4所示。表5.1.3 各个区的Eri及对应污染程度区域AsCdCrCuHgNiPbZn117.466.94.518.7106.37.511.13.4轻微中等轻微轻微强轻微轻微轻微220.190.73.448.3734.18.1
13、15.04.0轻微强轻微中等极强轻微轻微轻微311.235.22.66.646.86.35.91.1轻微轻微轻微轻微中等轻微轻微轻微415.983.13.723.6510.77.210.23.5轻微强轻微轻微极强轻微轻微轻微517.464.72.811.4131.46.29.82.2轻微中等轻微轻微强轻微轻微轻微表5.1.4 各个区的RI及对应污染程度区域12345RI235.87923.82115.50657.83246.05污染程度中等很强轻微很强中等从上面两个表格便能够得到该城区内不同区域重金属的综合污染程度及不同重金属在各个区域内的污染程度。可以得到以下结论：a、As、Cr、Ni、Pb
14、、Zn、Cu六种重金属对各个区域的污染都很小；b、Hg的污染很强；c、Cd对不同区域的污染强度不同，对工业区和主干道区污染最强；d、工业区和主干道区污染的综合污染最强。5.2 问题二模型的建立与求解5.2.1 模型的建立利用主成分分析法4来确定8种重金属中占主要原因的因素，其核心思想是降维，即用较少的变量替代较多的变量，用这些较少的变量来反映较多的变量反映的情况，将问题简单化。具体的步骤如下：(1) 标准差标准化标准差标准化是将某变量中的观察值减去该变量的平均数，然后除以该变量的标准差，即： (5.2.1)经过标准差标准化后，各变量将有约一半观察值的数值小于0，变量的平均数为0，标准差为1。经
15、标准化的数据都是没有单位的纯数量。对变量进行的标准差标准化可以消除量纲(单位)影响和变量自身变异的影响。(2) 相关系数矩阵 (5.2.2)式中，为原变量的与之间的相关系数，其计算公式为： (5.2.3)因为R是实对称矩阵，所以只需计算上三角元素或下三角元素即可。(3) 特征值与特征向量首先解特征方程，通常用雅可比法(Jacobi)求出特征值，并使其按大小顺序排列，即；然后分别求出对应于特征值的特征向量。这里要求，其中，表示向量的第j个分量。(4) 主成分贡献率及累计贡献率主成分zi的贡献率为： (5.2.4)累计贡献率为： (5.2.5)一般把累计贡献率达8595%的特征值所对应的第1、2，
16、m(mp)个成分规定为主成分。(5) 主成分载荷计算公式为: (5.2.6)根据主成分在原变量中所占的载荷，计算出对应的变量。5.2.2 模型的求解本题要求分析数据，计算污染的主要原因，就是对这8种污染源进行分析，得出污染土壤的主要影响因素，按照上面的模型，借助SPSS软件，得出结果：表5.2.1 生活区主成分的分析表解释的总方差主成分成份初始特征值提取平方和载入合计方差%累积%合计方差%累积%123As3.61645.19945.1993.61645.19945.1990.669-0.646-0.010Cd1.13314.16559.3651.13314.16559.3650.7840.17
17、1-0.417Cr1.07513.43272.7971.07513.43272.7970.6430.2340.493Cu0.80710.08382.8800.80710.08382.8800.729-0.2460.024Hg0.5246.54489.4240.5246.54489.4240.4920.130-0.437Ni0.4455.55694.9800.4455.55694.9800.686-0.2530.523Pb0.2362.95497.9340.2362.95497.9340.8030.112-0.348Zn0.1652.066100.0000.1652.066100.0000.50
18、10.6910.267从表中看出，主成分1与Pb呈现较强的正相关，主成分2与Zn呈现较强的正相关，主成分3与Ni呈现较强的正相关，综合可得，Pb，Zn，Ni就可以在很大程度上反映原变量的情况，可以确定为主成分。其中Pb的主要污染源是生活区中交通尾气的排放，Zn的主要来源是工业污染，从第一问作的空间分布图中可以看出，生活区比较接近工业区，工业区污染是造成生活区收到污染的原因之一。表5.2.2 工业区主成分的分析表解释的总方差主成分成份初始特征值提取平方和载入合计方差%累积%合计方差%累积%12As5.25465.67065.6705.25465.67065.6700.5180.758Cd1.26
19、315.78681.4551.26315.78681.4550.7860.074Cr0.7819.76391.2180.7819.76391.2180.916-0.206Cu0.2673.33694.5540.2673.33694.5540.868-0.463Hg0.2272.83497.3880.2272.83497.3880.845-0.459Ni0.1491.86899.2560.1491.86899.2560.7670.421Pb0.0510.63499.8900.0510.63499.8900.8580.049Zn0.0090.110100.0000.0090.110100.0000
20、.8590.188从表中可以看出，主成分1与Cr呈现出较强的正相关，主成分2与As呈现较强的正相关，综合可得，Cr，As可以在很大程度上反映原变量的情况，可以确定为主成分。由于该区是工业区，工厂会比较多，Cr是一种重要的环境污染元素，主要来源是工厂排放的“三废”，As也是一种工业废水中排放的重金属元素，所以该区的土壤的重金属污染主要是由工业污染引起的。表5.2.3 山区主成分的分析表解释的总方差主成分成份初始特征值提取平方和载入合计方差%累积%合计方差%累积%123As3.04238.02238.0223.04238.02238.022-0.0090.6680.478Cd2.03625.445
21、63.4672.03625.44563.4670.601-0.6780.141Cr1.54919.36282.8291.54919.36282.8290.7610.480-0.365Cu0.6898.61691.4450.6898.61691.4450.5170.4690.615Hg0.2493.11694.5620.2493.11694.5620.324-0.0740.747Ni0.2312.88397.4450.2312.88397.4450.7370.499-0.409Pb0.1612.01799.4630.1612.01799.4630.605-0.6300.173Zn0.0430.5
22、37100.0000.0430.537100.0000.905-0.173-0.183从表中可以看出，主成分1与Zn呈现出较强的正相关，主成分2与Cd呈现出较强的负相关，主成分3与Hg呈现出较强的正相关，综合可得，Zn，Cd，Hg可以在很大程度上反映原变量的情况，可以确定为主成分。主要是由于锌矿的开采和运输等造成的。Cd主要存在于含Zn的矿石中，而Hg主要是由工业生产引起的，也有可能是因为工业区扩散。表5.2.4主干道路区主成分的分析表解释的总方差主成分成份初始特征值提取平方和载入合计方差%累积%合计方差%累积%123As3.75146.88346.8833.75146.88346.8830.
23、235-0.1310.939Cd1.28716.08462.9671.28716.08462.9670.6210.4580.011Cr0.99412.42475.3910.99412.42475.3910.874-0.319-0.151Cu0.7499.35984.7510.7499.35984.7510.906-0.250-0.196Hg0.6858.56093.3110.6858.56093.3110.1700.7430.022Ni0.3364.19597.5050.3364.19597.5050.888-0.320-0.023Pb0.1161.45098.9550.1161.45098.
24、9550.7030.477-0.057Zn0.0841.045100.0000.0841.045100.0000.6430.1150.217由该表可以看出：主成分1与Cu呈现出较强的正相关，主成分2与Hg呈现出较强的正相关，主成分3与As呈现出较强的正相关，综合可得，Cu，Hg，As可以在很大程度上反应原变量的情况，可以确定为主成分。化石燃料的燃烧，城市垃圾污染是造成Cu的主要污染源，矿物燃料的燃烧会造成Hg的污染。As一般是由工厂三废引起的。表5.2.5公园绿地区主成分的分析表解释的总方差主成分成份初始特征值提取平方和载入合计方差%累积%合计方差%累积%123As3.90748.84048.
25、8403.90748.84048.8400.636-0.5690.349Cd1.61520.18269.0231.61520.18269.0230.8110.103-0.232Cr1.05613.20282.2251.05613.20282.2250.809-0.426-0.003Cu0.5466.82989.0540.5466.82989.0540.6790.457-0.151Hg0.3364.20393.2570.3364.20393.2570.2030.3540.889Ni0.2423.03096.2880.2423.03096.2880.663-0.6220.012Pb0.2132.6
26、5998.9470.2132.65998.9470.7820.5440.097Zn0.0841.053100.0000.0841.053100.0000.7980.286-0.242由表可以看出：三个主成分与Cd，Hg呈现出较强的相关性，可以确定为主成分。工业污染和交通运输是造成Cd污染的主要原因，公园绿地一般在交通区附近，所以交通运输和尾气排放等会引起绿地区的污染。矿物燃料的燃烧会造成Hg的污染。综上所述，该地区的污染原因主要是工业污染、尾气排放和矿石开采引起土壤的重金属污染。5.3 问题三模型的建立与求解5.3.1 模型的建立重金属污染物主要以泥沙颗粒为载体迁移转化5，在自然情况下，假设污
27、染源的强度是一定的，并且忽略风向、地壳运动等不稳定因素的影响，因此，对于重金属污染物的传播特性，可以运用粒子源在自由空间内的扩散模型来确定。由粒子的扩散特性，先建立一维有源扩散模型。已经假设污染源的强度是一定的，且各向扩散系数相同，则扩散满足一维泊松方程，如下： (5.3.1)式中，负号表示扩散方向总是沿高浓度指向低浓度。对(5.3.1)式积分两次，便能够得到一维有源的扩散模型，如下： (5.3.2)式中，m1、m2为积分产生的待定常数，r为污染点到污染源的直线距离，即：从第一问中的空间分布图可以得到污染源的大概位置及个数，取每个污染源周围大概3公里范围内采样点来求解确定该污染源的位置，如图5.3.1所示。图5.3.1 采样点选取示意图利用最小二乘法便可以确
