2011高教社杯全国大学生数学建模竞赛
承 诺 书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 河南科技大学 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 2011 年 9 月 12 日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2011高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）： 赛区评阅记录（可供赛区评阅时使用）： 评 阅 人 评 分 备 注 全国统一编号（由赛区组委会送交全国前编号）： 全国评阅编号（由全国组委会评阅前进行编号）：
城市表层土壤重金属污染分析
摘要
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。城市工业、经济的发展，污水排放和汽车尾气排放等均能引起城市表层土壤重金属污染。而重金属污染对城市环境和人类健康造成了严重的威胁，因此对城市表层土壤重金属污染的研究具有重大意义。
对于问题1，先用MATLAB软件对所给数据进行处理，插值拟合得出8种主要重金属元素在该城区的空间分布图；再用内梅罗综合污染指数评价法建立模型进行求解。首先用EXCEL对数据进行分析，得出各区的8种重金属的平均浓度；然后结合MATLAB软件求出各区的单项污染指数和综合污染指数，进而得出各区的综合污染等级，如下表： 区域 生活区 工业区 山区 主干道路区 公园绿地区 污染等级 重污染 重污染 轻度污染 重污染 中度污染 对于问题2，先借助SPSS软件对各种重金属元素的浓度和海拔做相关性分析，得出各种元素之间及其与海拔之间的相关系数矩阵和相关度；然后结合第一问给出的空间分布图和区域散点图，参照主要重金属含量土壤单项污染的指数，分析得出各重金属污染的主要原因主要来自工业区、主干道路区和生活区。
对于问题3，由上述问题的分析可以认为重金属的分布是连续的，物质的扩散从高浓度向低浓度进行。在模型一数据处理基础上建立遍历搜索模型，结合MATLAB软件求出重金属空间分布中的极值点即可能的污染源，得出极值点后再结合《国家土壤环境质量标准》通过MATLAB软件对极值点进行筛选，得出8种重金属元素的主要污染源。
对于问题4，对所建立的模型进行分析，找出了各个模型的优缺点。然后分析影响城市地质演化模型的因素，为更好地研究城市地质环境的演变模式，从动态和多元的角度出发，还应搜集采样点的长期动态数据和岩石、土壤、大气、水和生物等因素的相关信息，分别建立动态动态传播模型和城市地质环境的综合评价预测模型。
关键词：梅罗综合污染指数评价法 污染等级 相关矩阵 遍历搜索模型 污染源
一 、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0-10 厘米深度）进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自
然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息有了这些信息，如何建立模型解决问题
二、模型假设
1、假设题目中所给数据可靠无误。
2、假设问题一中各区平均的污染程度可以看做该区的污染程度。 3、假设问题二中只考虑题目中所给的8中重金属，不考虑其它重金属。 4、假设重金属传播特征不受风向等因素影响。
三、问题分析
(一) 问题1的分析：
问题1属于空间分布和综合评价问题，重金属的传播过程是一个扩散的过程，通常物质扩散模型中物质从高浓度向低浓度扩散且其浓度的分布是连续的，据此我们可以用附表中所给的采样点污染数据为基础借助MATLAB软件[1]进行插值拟合得出8种主要重金属污染物在整个城区的空间分布图。对于该城区内不同区域重金属的污染程度的研究可以借助我国《土壤监测技术规范》(HJ/T 166-2004)[2]中推荐的内梅罗综合污染指数法进行评价，求出不同区域重金属的污染等级。 (二) 问题2的分析：
问题2要求通过数据分析来说明重金属污染的主要原因。首先可以对重金属和海拔进行相关性分析，得出相关矩阵和相关度，再结合问题一求出的结论分析出重金属可能的主要来源和重金属污染的主要原因。 (三) 问题3的分析：
由问题一的分析我们得知重金属的分布是连续的，同时我们还可以知道物质的扩散是从高浓度向低浓度进行的，在扩散模型中某区域浓度最高的点可能就是扩散源，所以重金属空间分布中的极值点就可能是重金属的传播模型中污染源。因此问题三的求解就转化为在模型一所拟合出的重金属空间分布曲面上搜索极值的问题。搜索极值的现代算法有模拟退火，遗传算法，鱼群算法等多种。考虑的模型中所搜索的域有限，且目标解数目不确定，遍历搜索是较好的方法。得出极值点后再结合国家土壤环境质量标准筛选出污染源。
(四) 问题4的分析：
首先应对问题一，二，三所建立的模型进行优缺点分析然后根据影响城市演化模型的因素，分析还应搜集的数据以及模型如何建立的问题。
四、符号说明
符号设定 Pij Cj Sj PN Pj,ave Pj,max z 符号说明 区域i中第j个重金属的污染分指数 第j个重金属的实测浓度 第j元素的评价标准 综合污染指数 平均单项污染指数 最大单项污染指数 浓度分布矩阵 注：在此没有设定的符号在下文中会具体说明。 五、模型的建立及求解
一、问题一的求解：
1.1 用MATLAB软件对所给数据进行插值拟合得出调查区的地形图和8种主要重金属元素在该城区的空间分布图，再用MATLAB软件对所给数据进行分析得出功能区散点图：
图1：调查地区的地形图 图2：功能区散点图
图3：砷和镉在该城区的空间分布图 图4：铬和铜在该城区的空间分布图 图5：汞和镍在该城区的空间分布图 图6：铅和锌在该城区的空间分布图
说明：图1的Z轴为海拔高度，X、Y轴为地理坐标值（单位：m）。
图2 的X、Y轴为地理坐标值（单位：m）。
图3-图6的Z轴为重金属元素的浓度（单位：μg/g），X、Y轴为地理坐标值（单位：m）。
1.2 模型建立：
土壤环境质量单项污染指数主要用来评价某一污染物的污染程度，指数小污染轻，指数大污染则重。但区域内土壤环境质量作为一个整体和外区域进行比较时除用单项污染指数外，还常用综合污染指数。综合污染指数可以综合判断某土壤多种污染物的联合污染效应。
目前土壤环境质量评价方法有很多，各有优点和缺点。本文根据我国《土壤监测技术规范》(HJ/T 166-2004)[2]中推荐的内梅罗综合污染指数法进行评价。在计算某个区域某种重金属单项污染指数(分指数)的基础上，再计算该区域多种重金属的综合污染指数。单项污染指数和综合污染指数的计算公式如下：
Pij=Cj/Sj （1）
PN=(P2j,ave+P2j,max)/2 （2）
当Pij≤1时，表示土壤未受该因子污染，当Pij >1时，表示土壤受该因子污染。内梅罗综合污染指数反映了各污染物对土壤的作用，同时突出了高浓度污染物对土壤环境质量的影响。根据HJ/T 166-2004,内梅罗综合污染指数的分级标准（见表1），得出各个区域的污染等级。
表1：内梅罗综合污染指数的分级标准 等级 1 2 3 4 内梅罗污染指数 PN≤0.7 0.7＜PN≤1.0 1.0＜PN≤2.0 2.0＜PN≤3.0 污染等级 清洁（安全） 尚清洁（警戒线） 轻度污染 中度污染 5 PN＞3.0 重污染 1.3 模型求解： 本文以背景值作为评价标准进行求解，用EXCEL对文中所给数据进行分类，把数据分入1类区、2类区、3类区、4类区、5类区。然后得出各个区里面主要重金属含量的平均值，可看作各区中主要重金属含量值。如下表：
表2：各区重金属含量的平均值 As Cd Cr Cu Hg Ni Pb Zn 区域 (μg/g) (ng/g) (μg/g) (μg/g) (ng/g) (μg/g) (μg/g) (μg/g) 1 2 3 4 6.27 289.96 69.02 49.4 93.04 18.34 69.11 237.01 93.04 277.93 36.56 73.29 7.25 393.11 53.41 127.54 642.36 19.81 4.04 152.32 38.96 5.71 360.01 58.05 17.32 40.96 15.45 62.21 446.82 17.62 63.53 242.85 6.26 280.54 43.64 30.19 114.99 15.29 60.71 154.24 5 然后根据公式（1）、（2）结合MATLAB软件算得各区重金属单项污染指数和综合污染指数，如下表：
表3：各区重金属单项污染指数和综合污染指数 综合污染区域 单 项 污 染 指 数 指数 1 2 3 4 5 As Cd Cr Cu Hg 2.6583 18.3531 1.1703 12.7663 3.2854 Ni Pb Zn 4.028 3.1704 13.5331 1.2532 9.4264 2.7343 1.7417 2.2305 2.2265 3.7424 2.0139 3.0239 1.7229 9.6621 1.1222 1.1717 1.2568 1.3121 1.5861 2.7693 1.8726 4.7129 1.7389 2.158 1.4077 2.2871 1.4911 2.2294 3.4349 1.6106 3.0013 1.2561 1.1794 1.0622 1.4325 2.0494 3.5196 1.2431 1.9584 2.2354 再由内梅罗综合污染指数的分级标准得出各区的综合污染等级，如下表： 表4：各区综合污染等级 区域 污染等级 生活区 5 重污染 工业区 5 重污染 山区 3 轻等污染 主干道路区 5 重污染 公园绿地区 4 中等污染 从表中可以看出，该城区内生活区、工业区、主干道路区属于重污染区，公园绿地区属于中等污染区，山区属于轻度污染区。 二、问题二的求解： 2.1 模型建立：
用SPSS11．0统计软件对各种重金属元素浓度和海拔做相关性分析，得出各种元素与元素之间和元素与海拔之间的相关系数矩阵及其相关性，结合第一问得出的空间分布图和区域散点图，参照主要重金属含量土壤单项污染的指数，分析得出各重金属污染的主要原因。
2.2 模型求解：
2.2.1 根据题中所给数据，以As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种重金属元素浓度和海拔作相关性分析，经SPSS11．0统计软件进行相关性分析，得出该市表层土壤As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种重金属原始含量数据和海拔的相关系数矩阵，如图所示
图7: 重金属原始含量数据和海拔的相关系数矩阵
可见各重金属浓度均和海拔成负相关，即海拔越高，其含各种重金属浓度越低；Cr和Ni的相关性最好，相关系数最大，为0.716，其次为Pb和Cd，相关系数为0．660，以下是Cr和Cu的相关性较好，相关系数是0.532，其它元素之间的相关性并不是很好。从成因上来分析，相关性较好的元素可能在成因和来源上有一定的关联。结合第一问中8种主要重金属元素在该城区的空间分布可以看出，Cr和Ni、Pb和Cd可能是来自同一来源。 2.2.2 根据空间分布图、区域散点图和主要重金属含量土壤单项污染的指数进行分析：
对于Cr和Ni，在来源上关联较密切，该市表层土壤Cr和Ni基本未污染，只有个别点富集程度较高，污染达到中度污染，该富集中心的位置主要分布在生活区周边和主干道区周边，这可能是由于生活废水的排放和交通源汽车尾气的排放等原因造成的。
对于Pb和Cd，在来源上关联较密切， Pb和Cd的高含量点主要分布在交通繁忙的主干道路区周边和工业区周边,这可能是因为Pb和Cd来自该市中心交通源汽车尾气的排放、汽车轮胎的磨损和冶炼厂的废水、尘埃和废渣,以及电镀、电池、颜料、塑料稳定剂、涂料工业的废水等。所以可以说Pb和Cd的污染主要是由于主干道污染和工业污染。
对于Cu,该市表层土壤Cu基本未污染，只有个别点富集程度较高，污染达到中度污染，该富集中心的位置主要分布在生活区周边，这可能是由城市商业活动、城市居民生活累加到土壤中的Cu。
对于Hg, 其高含量点主要分布在交通繁忙的主干道路区周边和工业区周边, Hg污染
的一个主要原因是由于燃煤造成的，无论是工业用煤还是居民用煤，而且燃烧方式落后。工业排放也是表层土壤Hg污染的另一个重要来源，主要在大面积污染的几个工业浓集中心。
对于Zn,其高含量点也主要分布在交通繁忙的主干道路区周边和工业区周边, 这主要是由于汽车尾气的排放和厂矿企业的三废排放。
对于As，该市表层土壤As基本都是轻度或中度污染，只有个别点富集程度较高，该富集中心的位置主要分布在工业区周边，主要来源可能是工厂的废水排放。
综上所述，可以认为工业区、主干道路区和生活区的活动是造成该城区表层土壤重金属污染的主要原因。 三、问题三的求解： 3.1 模型建立：
依据问题一得出的各重金属元素在该城区的空间分布，得到浓度分布矩阵Z（Z是100×100的矩阵），进而结合MATLAB软件建立搜索模型。
Z是100×100的矩阵，借鉴元胞的思想建立一个100×100规模的二维网格，将元素浓度分布矩阵对应放入，其中每一个元素占据其中一个格子。根据问题分析可知：污染源存在于二维网格中的某些格子中。并且污染源所在格子元素浓度大于周围格子的元素浓度。二维元胞自动机(规则四方网格划分)的邻居通常有以几种形式如图2所示：黑色元胞为中心元胞，灰色元胞为该元胞的邻居。（参考文献[4]）
图8：元胞邻居模型
分析三种邻居模型发现第二种模型最适合。第二种邻居模型中污染源存在的格子z(i,j)应满足：
对于与外界不相邻的格子
z(i,j)>z(i-1,j) （3）
z(i,j)>z(i+1,j) （4）
z(i,j)>z(i,j-1) （5） z(i,j)>z(i,j+1) （6） z(i,j)>z(i-1,j-1) （7） z(i,j)>z(i-1,j+1) （8） z(i,j)>z(i+1,j-1) （9） z(i,j)>z(i+1,j+1) （10）
对于边界处的格子理论上应满足 以左边界为例
z(i,1)>z(i-1,1) （11） z(i,1)>z(i-1,2) （12）
z(i,1)>z(i,2) （13） z(i,1)>z(i+1,2) （14） z(i,1)>z(i+1,1) （15）
对于顶角处的格子理论上应满足： 以左边界为例
z(1,1)>z(1,2) （16） z(i,j)>z(2,2) （17） z(i,j)>z(2,1) （18）
为了简化模型在此不予考虑，即认为对于边界和顶角处不存在污染源。 通过搜索模型可以求出重金属空间分布中的极值点即可能的污染源，再结合国家土壤环境质量标准[3]（下表5）通过MATLAB软件对极值点进行筛选出，求出重金属的主要污染源。
表5：国家土壤环境质量标准 级别 As (μg/g) Cd (ng/g) Cr (μg/g) Cu (μg/g) Hg (ng/g) Ni (μg/g) Pb (μg/g) Zn (μg/g) 一级 二级 三级 15 25 30 200 300 1000 90 300 400 35 100 400 150 500 1500 40 50 200 35 300 500 100 250 500 3.2 模型求解： 3.2.1 根据问题一中得出的砷元素在城区的空间分布（见图3）。得到浓度分布矩阵Z（矩阵较大未附出），结合MATLAB软件建立搜索模型进行搜索得出砷元素在空间分布极大值61个。用同样方法得出其他7种重金属在空间分布极大值个数（见表6）。
表6：八种金属元素空间分布极大值个数 As Cd Cr Cu Hg Ni Pb Zn 元素 61 60 57 62 63 60 53 58 个数 3.2.2 运用scatter函数画出各重金属元素空间分布极大值点的散点图使数据可视化。得到各种重金属元素空间分布极大值点的散点图如下：
图9：砷、镉、铬、铜的空间分布极大值散点图 图10：汞、镍、铅、锌的空间分布极大值散点图
3.2.3 结合国家土壤环境三个等级的质量标准通过MATLAB软件对极值点进行分级筛选：首先用国家土壤环境一级质量标准进行筛选，得出筛选结果，再用国家土壤环境二级质量标准对一级指标得出的点进行筛选，依次类推，最终得到筛选结果如表7
表7：不同国标等级下的极大值个数 元素 As Cd Cr Cu Hg Ni Pb Zn 一级个数 6 57 13 45 39 6 53 57 二级个数 1 52 3 20 17 3 3 25 三级个数 1 11 2 5 11 0 1 18 从筛选的结果中选出适当的点作为重金属的主要污染源，所选点个数和点坐标如下列各表：
表8：重金属主要污染源个数 元素 As Cd Cr Cu Hg Ni Pb Zn 个数 6 11 3 5 11 6 3 18 表9：砷污染源二维坐标及其浓度值 As (μg/g) 15.061 23.641 16.121 23.175 30.032 18.971 X/m 18900 12900 7200 4500 18300 27600 Y/m 2200 3200 7400 7800 10200 12200 表10：镉污染源二维坐标及其浓度值 Cd (ng/g) 1068.8 1458.6 1401.9 1321.9 1121.4 1054.9 X/m 4500 2400 2400 17700 17700 5100 Y/m 2600 3400 3600 4000 4200 5200 Cd (ng/g) 1054.9 1264.4 1024 1267.8 1263.8 1578.6 X/m 5100 3600 6000 4800 4800 21600 Y/m 5200 6000 8600 11200 11400 11600 表11：铬污染源二维坐标及其浓度值 Cr (μg/g) 747.81 304.81 976.76 X/m 4800 10800 3600 Y/m 4800 5600 6000 表12：铜污染源二维坐标及其浓度值 Cu (μg/g) 2759.4 2609.8 2622.3 2565.2 1391.9 X/m 2400 2700 2400 2700 3600 Y/m 3600 3600 3800 3800 6000 表13：汞污染源二维坐标及其浓度值 Hg (ng/g) 16385 14487 15460 15427 1839 2333 13434 13411 11432 1692 1723 X/m 3000 13800 2700 2700 7200 3300 15300 15300 15600 22500 8700 Y/m 2600 2600 3400 3600 7400 8200 9200 9400 9400 10600 12200 表14：镍污染源二维坐标及其浓度值 Ni (μg/g) 146.08 70.587 69.355 X/m 3600 22200 27600 Y/m 6000 12200 12200 表15：铅污染源二维坐标及其浓度值 Pb (μg/g) 527.92 485.07 354.06 X/m 2100 5100 3600 Y/m 3400 5200 10600 表16：锌污染源二维坐标及其浓度值 Zn (μg/g) 1485.6 550.9 1631.5 1457 1749.2 3092.1 2801.4 1965.3 1961.9 X/m 4500 8100 12900 12900 2400 9600 9600 3600 3600 Y/m 2600 3200 3200 3400 3600 4600 4800 5800 6000 Zn (μg/g) 1111.9 1064.7 5320 5859 5265 5379 3664.5 2985.5 552.1 X/m 5400 5400 12900 8100 6000 9600 13800 13800 6000 Y/m 7200 7400 7800 8400 8600 8600 9800 10000 11000 运用scatter函数画出各重金属元素主要污染源的散点图使数据可视化。各种重金属元素主要污染源的散点图如下：
图11:砷、镉主要污染源的散点分布图 图12:铬、铜主要污染源的散点分布图 图13:汞、镍主要污染源的散点分布图 图14:铅、锌主要污染源的散点分布图
四、 问题四的求解：
优点：解决问题一的第一小问时，我们用MATLAB对原始数据进行差值拟合。由于所给数据采样点的不规则性，首先使用griddata函数对所给数据进行插值规整得出一个X,Y分别等步长的某种元素的浓度分布矩阵。在规整的浓度分布矩阵基础上分别使用pcolor，contourf，contour，surf等函数绘出了各种重金属元素在城区的空间分布。通过综合比较之后选定三维surf曲面建立重金属元素的空间分布模型，直观明了。解决问题二时，我们用SPSS对各种重金属元素浓度和海拔做因子分析，得出各种元素浓度和海拔相关性，各元素浓度和海拔呈现负相关，正好验证了第一问中求得的山区各重金属浓度最低，污染程度最轻这一结果。模型三中依据模型一中建立的浓度分布矩阵建立了遍历搜索模型。该模型能够有效且快速的找出空间极大值，即可能的污染源。然后结合国家土壤环境质量标准对污染源进行筛选，能方便的求出各种重金属元素的主要污染源。模型三的另一优点是可以根据筛选标准的高低，方便的区分不同污染源的污染程度的高低，有利于相关人员根据污染程度的高低采取不同的治理措施。
缺点：解决问题一第二小问时，我们把各区内采样点重金属浓度实测值的平均值用作各区重金属浓度的实测值，经过内梅罗综合污染指数评价法进行求解得出的各区污染等级只能反映各区的平均污染等级，不能反映各个采样点各自的污染等级。解决问题二时，我们忽略了该市风向、天气等因素对重金属污染的影响。
附表所给数据是静态的，无法根据所给数据建立城区污染的动态演化过程。为此我们还可以在原有采样点进行定期采样，获得重金属元素的动态传播模型。城市地质环境是一个涉及到地球岩石圈表层的岩石、土壤、大气、水和生物的复杂系统。为建立城市的地质演进模型，还应搜集岩石、土壤、大气、水和生物等因素的相关信息，进而建立城市地质环境的综合评价预测模型（参考文献[5]）。
六、模型的改进与推广
对于问题一所建模型，我们在求各区污染程度的时候，仅仅拿各区重金属的平均浓度进行分析，得出的结果只能反映各区的平均污染等级，不能反映各区在不同位置的污染等级。所以要想得出各区在不同位置的污染等级，需进一步求出各种重金属的空间分布函数。
对于问题二所建模型，我们在分析重金属污染的主要原因时，仅考虑了城市内各区造成的污染，忽略了该城市周边农田中农药的使用等因素造成的污染。要更好的分析出重金属污染的主要原因，我们还需对该城市周边农田中农药的使用等因素造成的污染进行调查分析。
对于问题三所建模型，我们在研究城市地质环境的演变模式时，我们仅对城市海拔进行了分析。而地质环境是一个涉及到地球岩石圈表层的岩石、土壤、大气、水和生物的复杂系统，为建立城市的地质演进模型，还应搜集岩石、土壤、大气、水和生物等因素的相关信息，进而建立城市地质环境的综合评价预测模型。
七、参考文献
[1] 张志涌，《精通MATLAB 6.5版》[M].北京：北京航天航空大学出版社， 234-302，2003。
[2] HJ/T166-2004，《土壤环境监测技术规范》[S]. 北京：中国标准出版社，2004。 [3] GB15618-1995，《土壤环境质量标准》[S].北京：中国标准出版社，1995。 [4] 祝红芳 王从庆，《机器人路径规划的元胞自动机算法》[J].江西：江西科学，第27卷第1期，36-40，2009.2。
[5] 周涛发 岳书仓 柏林，《城市地质环境及其评价与保护》[J].安徽：合肥工业大学学报，第20卷第3期，22-27，1997。
八、附录
附录一：
1.1 城区地形分布图的MATLAB程序：
A=xlsread('F:\\A\\cumcm2011A附件_数据.xls',1,'A4:E322'); x=A(:,2);y=A(:,3);z=A(:,4); scatter(x,y,5,z)%散点图 figure
[X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值 pcolor(X,Y,Z);shading interp%伪彩色图 title('功能区')
figure,contourf(X,Y,Z) %等高线图 figure,contour(X,Y,Z) title('功能区')
figure,surf(X,Y,Z)%三维曲面
1.2 功能区分布散点图的MATLAB程序：
A=xlsread('F:\\A\\cumcm2011A附件_数据.xls',1,'A4:E322'); x=A(:,2);y=A(:,3); z=A(:,5);
x1=find(z==1); x=x(x1(:)); y=y(x1(:));
scatter(x,y,20,'d') hold on;
x=A(:,2);y=A(:,3); x2=find(z==2); x=x(x2(:)); y=y(x2(:));
scatter(x,y,20,'h') hold on;
x=A(:,2);y=A(:,3); x3=find(z==3); x=x(x3(:)); y=y(x3(:));
scatter(x,y,20,'s') hold on;
x=A(:,2);y=A(:,3); x4=find(z==4); x=x(x4(:)); y=y(x4(:));
scatter(x,y,20,'p') hold on;
x=A(:,2);y=A(:,3); x5=find(z==5); x=x(x5(:)); y=y(x5(:));
scatter(x,y,20,'x') title('功能区分布')
legend('生活区','工业区','山区','主干道区','公园绿地区')
1.3 重金属在该城区空间分布图的MATLAB程序：
A=xlsread('F:\\A\\cumcm2011A附件_数据.xls',1,'A4:E322'); B=xlsread('F:\\A\\cumcm2011A附件_数据.xls',2,'B4:I322'); x=A(:,2);y=A(:,3); for k=1:8 z=B(:,k);
scatter(x,y,5,z)%散点图 figure
[X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值 pcolor(X,Y,Z);shading interp%伪彩色图 title('功能区')
figure,contourf(X,Y,Z) %等高线图 figure,contour(X,Y,Z) title('功能区')
figure,surf(X,Y,Z)%三维曲面 end
附录二：单项污染指数求解的MATLAB程序：
a=[6.27 289.96 69.02 49.4 93.04 18.34 69.11 237.01 7.25 393.11 53.41 127.54 642.36 19.81 93.04 277.93 4.04 152.32 38.96 17.32 40.96 15.45 36.56 73.29 5.71 360.01 58.05 62.21 446.82 17.62 63.53 242.85 6.26 280.54 43.64 30.19 114.99 15.29 60.71 154.24 ] function f=fun(a) m=size(a,1); n=size(a,2); c=[] b=[3.6
130 31 13.2 35 12.3 31 69 ]; b=b'; for i=1:5 for j=1:n;
c(i,j)=a(i,j)/b(j) end end
附录三：由各区的平均单项污染指数Pj,ave和最大单项污染指数Pj,max求各区的综合污染指数的MATLAB程序：
Pj,ave 2.46935 5.426975 1.19135 3.838588 2.03925 Pj,max 3.7424 18.3531 1.3121 12.7663 3.2854 function f=fun4(x)
a=((x(1)^2+x(2)^2)/2)^(1/2) 附录四：
4.1 重金属元素砷、镉、铬、铜的污染源分布图的MATLAB程序： A=xlsread('F:\\A\\cumcm2011A附件_数据.xls',1,'A4:E322'); B=xlsread('F:\\A\\cumcm2011A附件_数据.xls',2,'B4:I322'); ss={'As ','Cd','Cr ','Cu','Hg','Ni','Pb','Zn'}; x=A(:,2);y=A(:,3); for k=1:4 z=B(:,k);
[X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值 z=Z';
for i=2:99 for j=2:99 if
(z(i,j)>z(i-1,j))&&(z(i,j)>z(i+1,j))&&(z(i,j)>z(i,j+1))&&(z(i,j)>z(i,j-1))&&(z(i,j)>z(i-1,j-1))&&(z(i,j)>z(i-1,j+1))&&(z(i,j)>z(i+1,j-1))&&(z(i,j)>z(i+1,j+1));
z(i,j)=1000; end; end; end;
[ii,jj]=find(z==1000);
disp(ii');disp(jj');
subplot(2,2,k),scatter(ii,jj,'*'),title(ss{k}) end
4.2 重金属元素汞、镍、铅、锌的污染源分布图的MATLAB程序： A=xlsread('F:\\A\\cumcm2011A附件_数据.xls',1,'A4:E322'); B=xlsread('F:\\A\\cumcm2011A附件_数据.xls',2,'B4:I322'); ss={'As ','Cd','Cr ','Cu','Hg','Ni','Pb','Zn'}; x=A(:,2);y=A(:,3); for k=5:8 z=B(:,k);
[X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值 z=Z';
for i=2:99 for j=2:99 if
(z(i,j)>z(i-1,j))&&(z(i,j)>z(i+1,j))&&(z(i,j)>z(i,j+1))&&(z(i,j)>z(i,j-1))&&(z(i,j)>z(i-1,j-1))&&(z(i,j)>z(i-1,j+1))&&(z(i,j)>z(i+1,j-1))&&(z(i,j)>z(i+1,j+1));
z(i,j)=1000; end; end; end;
[ii,jj]=find(z==1000); disp(ii');disp(jj');
subplot(2,2,k-4),scatter(ii,jj,'*'),title(ss{k}) end
附录五：
筛选污染源的MATLAB程序（以砷为例，其它元素的程序只需调整参数即可）： A=xlsread('F:\\A\\cumcm2011A附件_数据.xls',1,'A4:E322'); B=xlsread('F:\\A\\cumcm2011A附件_数据.xls',2,'B4:I322'); ss={'As ','Cd','Cr ','Cu','Hg','Ni','Pb','Zn'}; x=A(:,2);y=A(:,3); z=B(:,1);
[X,Y,Z]=griddata(x,y,z,linspace(0,30000)',linspace(0,20000),'v4');%插值 z=Z'; k=z;
for i=2:99 for j=2:99 if
(z(i,j)>z(i-1,j))&&(z(i,j)>z(i+1,j))&&(z(i,j)>z(i,j+1))&&(z(i,j)>z(i,j-1))&&(z(i,j)>z(i-1,j-1))&&(z(i,j)>z(i-1,j+1))&&(z(i,j)>z(i+1,j-1))&&(z(i,j)>z(i+1,j+1));
z(i,j)=1000; end;
end; end;
[ii,jj]=find(z==1000); c=[ii,jj]; n=size(ii,1); for i=1:n
kk(i)=k(ii(i),jj(i)); end pp=kk';
kkk=find(kk>500); kkk=kkk'; n=size(kkk,1) for i=1:n
rr(i)=pp(kkk(i)); end rr
for i=1:n
cc(i,1)=c(kkk(i),1); cc(i,2)=c(kkk(i),2); end
scatter(cc(:,1),cc(:,2),'*') axis([0,100,0,100]) title(ss{1})
数学建模A题优秀论文
