建模A题城市表层土壤重金属污染分析竞赛论文

2021高教社杯全国大学生数学建模竞赛 承 诺 书 我们仔细阅读了中国大学生数学建模竞赛的竞赛规那么. 我们完全明白，在竞赛开始后参赛队员不能以任何方式〔包括 、电子邮件、网上咨询等〕与队外的任何人〔包括指导教师〕研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规那么的, 如果引用别人的成果或其他公开的资料〔包括网上查到的资料〕，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规那么，以保证竞赛的公正、公平性。如有违反竞赛规那么的行为，我们将受到严肃处理。

我们参赛选择的题号是〔从A/B/C/D中选择一项填写〕：
我们的参赛报名号为〔如果赛区设置报名号的话〕：
所属学校〔请填写完整的全名〕：
参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)：
日期：
年 月 日 赛区评阅编号〔由赛区组委会评阅前进行编号〕：
2021高教社杯全国大学生数学建模竞赛 编 号 专 用 页 赛区评阅编号〔由赛区组委会评阅前进行编号〕：
赛区评阅记录〔可供赛区评阅时使用〕：
评 阅 人 评 分 备 注 全国统一编号〔由赛区组委会送交全国前编号〕：
全国评阅编号〔由全国组委会评阅前进行编号〕：
城市表层土壤重金属污染分析 摘要 研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。本文就城市表层土壤重金属污染问题使用Matlab、Excel、Surfer等软件进行了多层次的分析，通过图像再现以及通过相关专业知识建立模型来解决问题。

首先根据表格中各采样点数据和各元素的背景值，计算出各个采样点的各种元素与相对应元素的背景值的比值。考虑到单个采样点上元素分布的随机性，我们规定当元素与背景值的比值小于1时全部取为1， 创造性地提出非污染法采样法来进一步处理数据。并用最终得到的有效数据通过计算机模拟重现8种主要重金属元素的空间分布，构建重金属污染综合指数来评估该城区不同区域重金属污染程度。

从相关专业知识出发建立大气降尘重金属积累模型来研究重金属污染物的传播特征好，求解污染源位置。我们找出各个功能区中污染比拟严重的集中元素，通过生活常识和查阅相关资料分析这些重金属元素的来源。我们选取该城区污染程度最高Hg作为代表元素分析其主要来源。查阅资料得出重金属主要通过大气沉降进入土壤。所以本文在建立Hg元素传播特征时我们主要考虑的是大气沉降。考虑到Hg在土壤中具有累积效应，且每年累积量得多少与采样点和污染源之间的距离有关。根据累积效应公式构造出每个采样点Hg含量与该采样点到污染源距离之间的关系模型，通过数据计算和分析最终确定了污染源的位置。

最后通过将数据所反映的现实与模型比拟来优化模型，以求更好地研究城市地质环境的演变模式。取得了较好的效果，除了误差小之外，该模型的最大优点在于可以对将来重金属污染趋势进行预测。

关键词：计算机模拟、污染指数 、大气沉降、 累积效应 一、问题的提出 随着城市经济的快速开展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。

现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土〔0~10 厘米深度〕进行取样、编号，并用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。

附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。

现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。

(2) 通过数据分析，说明重金属污染的主要原因。

(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。

(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？ 二、问题的分析 〔1〕重金属元素在该城区的空间分布即为其浓度的空间分布，根据不同采样点的浓度分布，利用Matlab软件画出其等高线，通过图像就可以清晰的看出每一种元素的空间分布图。为了得到不同区域内重金属的污染程度，我们可以建立污染指数的计算公式，得到每种元素在每个区内的污染指数，通过比拟就可以正确的分析不同区域内重金属的污染程度了。通过数据我们可以得到不同区域重金属的污染指数，通过对不同种类重金属元素赋合理的权重可以表示不同种类重金属元素对环境影响的程度的大小，两者结合可以得出每个区受重金属污染的综合指数。

〔2〕 我们用Excel将数据进行处理，把数据按功能区进行排序，这样就得到了各个功能区元素的分布情况，通过生活常识和查阅相关资料得出各个功能区重金属污染的主要原因。

〔3〕重金属主要通过大气沉降、随地表水流动、灌溉、雨水淋溶等四种方式进入土壤，而在这四种方式中，通过大气沉降为最主要的方式。通过前两问的分析，我们看出该城区有五片地区污染比拟严重，且各个区在空间上是相互独立的。我们猜测他们有各自的污染源。我们选取污染最严重的一个地区〔〕进行分析，且考虑到该区域点海拔相差不大，予以忽略。且该区域Cu、Pb等污染也比拟严重，他们之间可能具有一定的相关性。重金属在环境中具有累积效应，且累积量多少与采样点与污染源之间的距离有关，通过这一关系我们构造了Hg的这样累积模型结合Hg在该区域的分布模型和表格中的数据综合分析和计算最终确定污染源的位置。

〔4〕所有模型都有优点缺点，主要是从和实际情况相比拟而言的。选取模型时我们尽量先利用数学或相关专业知识进行分析，尽量用和实际情况相符的模型。当然模型不可防止的要进行很多简化或理想化处理而这些往往是产生误差的所在。对于收集信息来更好地研究城市地质环境的演变模式，我们觉得需收集大局部和城市地质环境的演变有关的量的信息，对城市地质环境的演变影响大的应该优先收集，其他稍微次要的信息也要收集。通过收集到的更完全的数据可以对已有模型进行合理优化，就可以更好的解决问题。

三、模型假设 1．假设题目所给的数据是完全正确，没有问题的；

2．假设重金属元素背景值对环境的污染根本没有，其临界值可以做为判定是否受到污染的标准；
不同区域上只有采样点某种重金属浓度超过临界值，才判定此区受这种重金属的污染；

3．假设八种主要重金属元素在传播上具有一定的相关性；

4．假设污染源产生的污染物重金属粒子粒径几乎都大于1微米，在重力和惯性力的作用下很容易沉降，进而污染周边地域。

四、符号说明 ---环境污染指数；

--第区采样点数目；

--第区的临界值；

--对应各区的重金属浓度在临界值以上的数目；

--Hg每年在土壤中的累积量；

--Hg在土壤中的残留率；

--区域土壤背景值；

--Hg的年输入量；

--采样点N年的Hg元素总累计量；

--采样点离污染源的距离；

--采样点的横纵坐标；

--污染源的坐标；

--采样点Hg年输入量常数。

五、模型的建立与求解 5.1问题一 为了得到8种主要重金属元素在该城区的空间分布，我们利用给出的8种主要重金属元素在采样点处的浓度及相应采样点的地理坐标来建立模型，画出8种主要重金属元素在该地区的空间分布图。

以As元素为例，各个取样点的X、Y轴坐标，以相应取样点的As元素浓度为Z轴值。利用Matlab [1]画出As元素在该地区的空间分布三维图。对As元素的浓度进行合理的设定，就可以画出As元素浓度的等高线，投影到XoY平面上，最终得到As元素污染浓度的二维空间分布图。如图1所示：
图1 As元素在该地区的空间分布 从图1中可以明显地看出As元素浓度的几个高点，也就是As元素的重污染区。第一个重污染区位于工业区的第317〔X=6182，Y=2005〕，第318〔X=5985，Y=2567〕取样点周围，As元素浓度分别到达21.87(μg/g)和23.72(μg/g)；
第二个位于交通区的第178〔X=12696，Y=3024〕取样点周围，As元素浓度到达23，72(μg/g)；
第三个位于交通区的第84〔X=17814，Y=10707〕取样点周围，As元素浓度到达30.13(μg/g)；
第四个位于山区的第100〔X=18993，Y=12371〕取样点周围，As元素浓度到达20(μg/g)。

同样可以相同的方法做出其他七种主要重金属元素在该地区的空间分布图。详见图2至图8。

图2　Pb元素在该地区的空间分布 图 3　　Ni元素在该地区的空间分布 图 4　　Hg元素在该地区的空间分布 图 5　　Cu元素在该地区的空间分布 图 6　　Cr元素在该地区的空间分布 图 7 Cd元素在该地区的空间分布 图 8 Zn元素在该地区的空间分布 评估该城区内不同区域重金属的污染程度时，我们以该地区这8种主要重金属的背景值为参考标准，计算出5个区域内8种主要重金属的污染程度指标。对于不同区域重金属的污染程度，常用的均值法是将采样点重金属浓度求均值，通过在背景值根底上产生的等级划分得到其污染程度，而事实上，临界值以下的点是非污染点，对区域并未造成污染，而由于受采样点中非污染点的影响，计算结果并不科学，有人用去点法解决这个问题，这样解决了非污染采样点的问题，但是对于污染较小的区域，大量的去点可能会造成更大的误差。因此，并不可行。为了解决这个问题，我们在此根底上模拟均值，建立了一个新的环境污染指数来表征重金属在不同区域的污染程度，其计算公式如下：
公式〔1〕 其中为第区采样点数目；

为第区得临界值；

为对应各区的重金属浓度在临界值以上的数目。

其特点是将低于临界值的点对临界值比值定义为标准值1，优点在于解除了临界值以下的采样点对模型造成的影响，缺点在于可能由于大量的取1而使高浓度点过于均值化，降低了不同的值间的差距，但事实上，并不影响各种元素在不同区域上污染程度的比拟与分析，因此此算法是可行的。结果如表格1所示：
表格 1 主要重金属污染程度指数 区 元 素 生活区 工业区 山区 交通区 公园绿地区 As 1.2419 1.4371 1.0596 1.1866 1.2369 Cd 1.6211 2.0980 1.0854 1.9647 1.6115 Cr 1.5348 1.2715 1.0994 1.3280 1.0695 Cu 2.4665 6.2650 1.1222 3.0841 1.4800 Hg 2.0214 12.7097 1.1114 8.9541 2.4690 Ni 1.0745 1.1692 1.1059 1.0789 1.0367 Pb 1.6943 2.1926 1.0940 1.5482 1.5116 Zn 2.5354 2.8931 1.0451 2.5715 1.7371 考虑到不同种类的重金属元素对环境污染的程度不同，例如很少量的Cd〔镉〕、Hg〔汞〕等就能对环境产生极大的污染，而Zn〔锌〕等要到达相比照拟大的浓度才能对环境产生显著影响。所以必须对表格1数据进行进一步处理。

我们知道不同的重金属元素对环境危害程度具有相对的差异性，我们参考了中华人民共和国环境保护部给出的土壤环境质量标准GB15618-1995[3]与Hakanson指数法[4]，参考之前对这方面研究的并结合生活常识对不同种类的主要重金属元素赋合理的权重来表示不同种类的主要重金属元素对环境影响的程度的大小。这8种主要重金属元素的权重如表格2所示：
表格 2 各金属元素的权重 元素 As Cd Cr Cu Hg Ni Pb Zn 权重 0.114 0.340 0.02 0.057 0.341 0.057 0.057 0.013 将表格1中各个区域的各种主要重金属元素污染程度指数乘以相应的权重，就可以得到对每个区域主要重金属污染的综合评价，最后得到的结果如表格3所示：
表格 3 各区域重金属污染的综合评价 区域 生活区 工业区 山区 交通区 公园绿地区 综合污染指数 1．74 5．820 1.095 4.241 1．805 从表格3中可以看出工业区和交通区的主要重金属污染程度最为严重，这可能和工业生产过程中排放的大量废水、废气、废渣以及交通区的汽车排放的大量尾气和车辆磨损有关；
生活区污染程度较高，可能是有城市的大量生活垃圾导致的；
公园绿地区污染程度只比生活区稍微高一些，公园绿地区的环境质量并不像人们所想象的那样是城市中的最好的地方，可能是因为公园绿地区植被丰富，土壤有机物含量高更容易吸附重金属，再和生活区受到大致相同程度的人类活动影响时，重金属污染反而高一些；
山区受到人类活动影响的程度最低，污染最为轻微的地方也是山区。

5.2问题二 在城市开展过程中，重金属或含重金属的污染物会通过各种途径进入土壤，造成城市范围内土壤的严重污染。而距离城市越近，其所受到的污染就越严重，正如该城市中的情况，由表1可看出，该城市五个区域除了山区由于距离城市较远污染不大外，其余四个区域都受到了一定程度的污染，特别是工业区和交通区，污染很严重。

图9 该地区的地形等高线图 A 根据数据我们计算得到表3，即综合污染指数在不同区域的分布为工业区>交通区>公园绿地区>生活区>山区。Hg污染程度在工业区远远高于其他区，而通过查阅资料我们可知，城市土壤表层中Hg的主要来源于煤燃烧和垃圾燃烧，而城市中煤燃烧主要在工业区。两者相吻合，所以我们可以做出判断工业区Hg污染主要来源于工业区大量的煤燃烧。我们还发现交通区中A点〔如下列图所示〕的Hg含量也很高，我们画出了该地区的地形图，由图我们可以看出A点地势较低而周围地势较高，Hg可以随着流水传播，所以我们猜测该点Hg污染严重的原因是由于地势过低。

对于Zn、Cu、Pb三种重金属，可以发现主要集中在工业区、交通区和生活区，且工业区最大，交通区次之，而这Pb元素是汽车尾气产生的重要重金属元素，而且Cu和Zn还是工业生产中产生的废水内常常含有的元素，Pb也是轮胎磨损产生的粉层中的污染元素。因此，我们推断这三种元素污染的主要原因为汽车尾气排放，工业生产废料废渣及轮胎磨损等产生的，后经自然降雨或人为因素才进入土壤中的。

查阅资料我们可知一般情况下土壤表层中Cr主要来源于工业废水废气废渣及工厂电镀排放得废液，受人为影响很小，表格中所给数据该城区的五个功能区的Cr含量相差并不大，这与我们得到的结论相符。同样，Ni与Cd元素，在采样点污染指数上，工业区也都大于其它区，且根据五个区的分布，靠近工业区，其值整体呈上升趋势，因此断定其污染依然是工业生产产生的废气废水废渣。

对于其他原因，如污泥排放，农药、化肥使用，污水灌溉等等也会对重金属的传播与污染造成一定的作用，比方As、Cu、Cd就常见于化肥农药中[2]。

综上所述，各种重金属污染的主要原因如下表4所示 污染物 污染原因 Cr、Hg、As、Pb 、Ni 矿产开采、工业排放的废水、废气、废渣 Cr、Hg、As、Pb 煤和石油天然气燃烧、垃圾燃烧 Cr、Cd、Ni、Pb、Cu、Zn 工业电化学排放的废水 Hg、Cd、Pb、Ni、Zn 塑料、电池、电子工业排放的废水 Cr、Cd 燃烧、化工制革工业排放的废水 Pb 汽车尾气 As、Cu、Cd 农药、化肥 5.3问题三 由于我们发现Hg对该城区污染尤为严重，且可以看出Hg污染严重的地区在工业区。所以说Hg污染来源于工业生产活动，且通过观察各种重金属元素在该城区分布我们可以看出在工业区Hg污染比拟严重的地方还是Cu、Pb等元素浓度较高的地方。所以可以认为这些元素之间具有一定的相性，所以在该问的解答中确定Hg污染源位置是极其重要的。

现在我们考虑Hg的传播特征，Hg通常造成土壤污染的途径有以下几种途径1、污染物随大气传输而迁移扩散。2、污染物虽地表水流动、补给、渗入而迁移。3、污染物通过灌溉在土壤中积累。4、固体废弃物受自然降水淋溶作用，转移或渗入土壤。5、固体废弃物受风力作用产生转移。Kloke.A.D.RSavcrbcck. [5]等调查研究发现，对工业兴旺地区大气沉降对土壤系统中重金属累积奉献率在各种外源输入因子中排在首位。大气汞通过沉降进入土壤后，被土壤中的粘土矿物和有机物的吸附或固定，富集于土壤表层，或为植物吸收而转入土壤，所以在一般情况下很难发生转移[6]。所以我们考虑Hg的传播特征时主要考虑它的大气沉降特征。

Hg元素在土壤中具有累积效应用公式表示,式中 --Hg每年在土壤中的累积量；

--Hg在土壤中的残留率；

--区域土壤背景值；

--Hg的年输入量。

年后，Hg在该地区的土壤中的累积量为；

；

公式〔2〕 公式〔3〕 在上式中每年的输入量既包括大气中的沉降量也包括自然淋溶迁移量。而在土壤中Hg背景值稳定的情况下Hg背景值每年的损失量等于自然淋溶量即能到达动态平衡。此时考虑土壤中Hg每年的输入量就是大气中Hg的沉降量。所以对上式进行修正得到下式 公式〔4〕 据研究，一般重金属在土壤中不易被自然淋溶迁移，残留率一般在左右，所以我们取=0.9. Hg随大气沉降，所以在污染源周围区域单位面积上每年Hg输入量与它们离污染源的距离成反比。我们假设，于是 。

公式〔5〕 取这一重金属污染区进行研究。因为研究区域采样点海拔相差不大，我们忽略了海拔的影响取为做外表距离，设污染源处坐标为，所以 。因为是关于时间的一个常数，也是一个与时间相关的常数，所以我们令。对于一确定的年份是确定的，于是公式〔5〕可表达为 公式(6) 由重金属空间分布图，我们先粗略的取为〔3000,3000〕，计算出一下表格：
表格 5 C值确实定 Hg含量〔ng/g〕 〔ng/g〕 800.00 749.00 1049 2127 950.00 899.00 0 1787 1040.00 989.00 1647 2728 1339.29 1288.29 4684 1364 13500.00 13449.00 2383 3692 16000.00 15949.00 2708 2295 我们取,当求得值之后我们再对修正得到的坐标为〔2260，3813〕。几次对进行修正的值如下表 表格 6 C值和污染源位置的修正 序号 1 〔2260,3813〕 2 〔2134，3762〕 3 〔2030，3713〕 由修正值我们可以看出越来越精确，所以我们确定污染源位置为〔2030,3713〕。我们观察重金属元素在该区的分布发现在点〔2030,3713〕附近元素的浓度等高线分布最密，这也证明了污染源点确定是比拟精确的，模型到达了一定的精度。同理采用相同模型我们可以确定出其他几个重金属污染区污染源的位置 5.4问题四 〔1〕分析8种主要重金属元素在该城区的空间分布时，我们采取用Surfer软件进行绘图，同时显示出各种主要重金属元素在该城区的空间分布的外表图和等值线图。能够非常清楚直观表现出各种重金属元素在该城区的独特的空间分布，便于开掘更多隐含信息。

分析该城区内不同区域重金属的污染程度时，我们充分地参考了当地的背景值和国家土壤标准，相应采用了非污染点采样法和Hakanson指数法，利用综合污染指数能非常系统的反映各个区域内整体受污染情况。缺点是不能很好的全面反映各个区域内单一的一种重金属的污染情况。

〔2〕对于第三问我们利用对数据进行的分析得到的图像、表格等并查阅相关资料得到的了以大气沉降为主的传播特性。并以此建立大气降尘重金属积累模型来确定污染源的位置。

该模型的主要优点是在长时间观测的模式下，重金属累积量随时间的不断加大，因为偶然因素而使每年的观测值受到的影响相对而言就会越来越小，精确度非常的高。第二个最为显著的优势在于此模型可以对未来相当长时间〔理论上可达10年〕的该城区的重金属污染变化趋势做出准确的预测。另外一个优势是可以通过预测来研究城市地质环境的演变模式。主要缺点是在对某一年的数据进行处理时得到的误差烧大。

〔3〕对于更好地研究城市地质环境的演变模式，只需有足够数量的、高质量的上述大气降尘重金属积累模型所需的数据即可。首先这就需要长期的观测记录大气颗粒物沉降中干沉降和湿沉降所占比例。其次还需要长期检测记录重金属在土壤中的年积累量、区域土壤背景值、重金属污染物的年输入量、重金属污染物在土壤中的残留率。再次还得记录当地工矿企业、交通运输、养殖种植等产业在该城市的分布变化以及技术进步或采用绿色能源带来的环境红利。最后还得记录当地水源分布，气候特征的状况。获得到这些信息后，只需要依据信息对该模型进行必要的修正就能很好地解决完问题。

六、参考文献 [1]张智星，MATLAB程序设计与应用，北京：清华大学出版社，2002. [2]杨维，高雅玲等，毗邻铁矿的千山景区土壤重金属污染分析与评价，第26卷：151-153，2021. [3]土壤环境质量标准GB15618-1995, ，2009年9月10日. [4]HAKANSONL.An ecological risk index for aquatic pollution control a sediment ecological approach [J],Water Research,1980,14(8):975-1001. [5]Kloke.A.D,R.Saverbeek.ChangingMetalCyclesandHumanHealth[M] .Bcrlin:springer-verlag,1984:113-141. [6]赵柯，大气降尘对土壤重金属累积量估算方法探讨，环境科学与管理，第32卷第11期：55-57,2007. 附录：
根底程序 x1=[4043;2427;4777;6534;4592;2486;3573;5375;7304;9328;8077;7056;15467;14844;16569;14298;21418;26453;26416;5101;5382;5503;5636;6605;11649;12591;13855;14862;15387;15810;17203;16301;17904;10395;15467;16428;16289;12153;11958;10800;9277;11121;12625;17981;0;1647;2383;4742;4948;5567;22674;5438;12734;14896;16947;5006;6395;7405;8446;7612;8866;9475;9212;8629;7776;7106;6423;7458;11646;12641;14000;14207;15140;16440;10022;9333;10856;12644;9036;10599;14325;17044;18413;19007;18738;23664;18993;19968;22535;27816;25361;24065;25998;27177;26073;24631;24702;25461;26086;26015;27700;27696;27346;26591;27823;27232;24580;24153;22965;24685;28654;24003;21684;22193;17079;15255;15007;15801;17008;21475;20261;19569;19411;23359;23238;22624;21703;12734;14405;14074;14262;20591;20983;20217;18906;18467;16607;15952;22605;23146;22046;23785;25981;27380;23325;26852;74;1373;1321;1049;2883;2708;2933;3526;5062;5868;5481;3299;5635;5394;5291;7004;7048;8180;9090;8049;8017;6869;7747;8457;9460;9062;9319;10631;10685;10643;11702;11730;11482;10700;10630;11678;11902;13244;12746;12855;13797;12442;13093;13920;15658;14177;12778;17087;17075;17962;17814;18134;17198;17144;18393;19767;21006;21091;22846;22304;21439;20554;20211;21072;20215;21766;25221;26424;24813;23198;4020;4026;5314;7093;7100;6837;7906;8045;8394;8403;8079;9663;9469;9178;9095;10225;10210;10340;11557;11415;12696;12400;13765;13694;16032;16872;17734;17005;18438;18954;18012;19072;20282;21450;19501;19909;21018;22176;5734;7912;9296;8622;9237;8307;8904;10547;10398;11529;11563;14065;12727;15198;15248;16267;16440;15412;14269;13277;13175;12632;14624;16629;18470;19041;17414;15748;25021;5985;4233;4741;16387;16061;15092;3518;3469;3762;3927;4153;3267;4684;5495;5664;5541;5451;15087;16823;18303;18556;20582;14173;15517;14482;14318;10352;9095;10510;13954;10142;17765;6924;4678;6182;7653]; y1=[1895;3971;4897;5641;4603;5999;6213;8643;5230;4311;6401;8348;8658;5519;6055;7418;10721;5577;6508;4080;3012;1127;133;374;3515;1063;3345;2524;729;2307;6218;8299;8287;11203;12080;9069;10072;12336;13313;13282;16148;16432;16259;18449;1787;2728;3692;7293;7293;6782;12173;3994;4015;1603;7487;8846;10443;10981;11200;11938;13143;12000;11305;12086;10613;9467;8831;8920;9381;9560;8970;9980;11101;13232;12204;14631;14727;14943;17538;17980;8666;10691;11721;11488;10921;9790;12371;12961;11293;5581;6423;7353;7032;7771;8807;9422;9522;9834;11094;12078;11609;11621;13331;13715;14737;14482;13319;12450;13535;14278;8755;15286;13101;12185;5894;5110;5535;3966;4775;8540;7586;7348;6934;5325;6502;4818;6591;10344;18032;16516;15129;13549;15862;17642;16346;17001;17365;18397;14301;15382;17634;17643;18051;18202;16701;16114;781;731;1791;2127;3617;2295;1767;4357;4339;4904;6004;6018;7965;8631;7349;6226;4600;4496;5365;5439;7210;7286;8260;8991;8311;7639;6799;6472;5528;4472;4480;5532;6354;8184;8774;8618;7709;7056;8450;8945;9621;4329;4339;5354;7594;6684;5799;11933;12924;12823;10707;10046;9810;9081;9183;8810;8819;9482;9149;10527;11383;11228;10774;10404;9951;12348;5795;8639;10799;13523;2990;3913;2060;1381;2449;3490;3978;3052;2035;1075;0;1288;2286;3299;3975;3821;2789;1764;1581;2585;3024;2060;1353;2357;3061;2798;3629;7212;6539;4874;4414;8519;8590;7555;6091;5300;5764;5492;9659;12840;13102;10638;9872;9726;8868;9591;10360;11243;10298;10987;7691;10100;9106;11058;12068;12982;12877;13204;12238;17949;14004;14481;14411;15769;15476;15728;16290;2567;895;6434;6609;7352;6936;2571;2308;2170;2110;2299;793;1364;1205;1653;2093;2757;3512;4207;7385;5588;6548;11941;17034;12692;13569;17133;16414;15314;5615;1662;3561;5696;3765;2005;1952]; z1=[43.9400;80.3500;94.7800;57.6500;744.4600;54.6400......;57.5100;138.3700;67.2200;77.2700;35.6600......36.19000;35.7600;57.3600]; nx1=linspace(min(x1),max(x1),44); ny1=linspace(min(y1),max(y1),44); [x1x1,y1y1]=meshgrid(nx1,ny1); zz1=griddata(x1,y1,z1,x1x1,y1y1,＇v4＇); contour(x1x1,y1y1,zz1); hold on; x=[4043;2427;4777;6534;4592;2486;3573;5375;7304;9328;8077;7056;15467;14844;16569;14298;21418;26453;26416;5101;5382;5503;5636;6605;11649;12591;13855;14862;15387;15810;17203;16301;17904;10395;15467;16428;16289;12153;11958;10800;9277;11121;12625;17981]; y=[1895;3971;4897;5641;4603;5999;6213;8643;5230;4311;6401;8348;8658;5519;6055;7418;10721;5577;6508;4080;3012;1127;133;374;3515;1063;3345;2524;729;2307;6218;8299;8287;11203;12080;9069;10072;12336;13313;13282;16148;16432;16259;18449]; x5=[0;1647;2383;4742;4948;5567;22674;5438;12734;14896;16947;5006;6395;7405;8446;7612;8866;9475;9212;8629;7776;7106;6423;7458;11646;12641;14000;14207;15140;16440;10022;9333;10856;12644;9036;10599]; y5=[1787;2728;3692;7293;7293;6782;12173;3994;4015;1603;7487;8846;10443;10981;11200;11938;13143;12000;11305;12086;10613;9467;8831;8920;9381;9560;8970;9980;11101;13232;12204;14631;14727;14943;17538;17980]; x2=[14325;17044;18413;19007;18738;23664;18993;19968;22535;27816;25361;24065;25998;27177;26073;24631;24702;25461;26086;26015;27700;27696;27346;26591;27823;27232;24580;24153;22965;24685;28654;24003;21684;22193;17079;15255;15007;15801;17008;21475;20261;19569;19411;23359;23238;22624;21703;12734;14405;14074;14262;20591;20983;20217;18906;18467;16607;15952;22605;23146;22046;23785;25981;27380;23325;26852]; y2=[8666;10691;11721;11488;10921;9790;12371;12961;11293;5581;6423;7353;7032;7771;8807;9422;9522;9834;11094;12078;11609;11621;13331;13715;14737;14482;13319;12450;13535;14278;8755;15286;13101;12185;5894;5110;5535;3966;4775;8540;7586;7348;6934;5325;6502;4818;6591;10344;18032;16516;15129;13549;15862;17642;16346;17001;17365;18397;14301;15382;17634;17643;18051;18202;16701;16114]; x3=[74;1373;1321;1049;2883;2708;2933;3526;5062;5868;5481;3299;5635;5394;5291;7004;7048;8180;9090;8049;8017;6869;7747;8457;9460;9062;9319;10631;10685;10643;11702;11730;11482;10700;10630;11678;11902;13244;12746;12855;13797;12442;13093;13920;15658;14177;12778;17087;17075;17962;17814;18134;17198;17144;18393;19767;21006;21091;22846;22304;21439;20554;20211;21072;20215;21766;25221;26424;24813;23198;4020;4026;5314;7093;7100;6837;7906;8045;8394;8403;8079;9663;9469;9178;9095;10225;10210;10340;11557;11415;12696;12400;13765;13694;16032;16872;17734;17005;18438;18954;18012;19072;20282;21450;19501;19909;21018;22176;5734;7912;9296;8622;9237;8307;8904;10547;10398;11529;11563;14065;12727;15198;15248;16267;16440;15412;14269;13277;13175;12632;14624;16629;18470;19041;17414;15748;25021;5985]; y3=[781;731;1791;2127;3617;2295;1767;4357;4339;4904;6004;6018;7965;8631;7349;6226;4600;4496;5365;5439;7210;7286;8260;8991;8311;7639;6799;6472;5528;4472;4480;5532;6354;8184;8774;8618;7709;7056;8450;8945;9621;4329;4339;5354;7594;6684;5799;11933;12924;12823;10707;10046;9810;9081;9183;8810;8819;9482;9149;10527;11383;11228;10774;10404;9951;12348;5795;8639;10799;13523;2990;3913;2060;1381;2449;3490;3978;3052;2035;1075;0;1288;2286;3299;3975;3821;2789;1764;1581;2585;3024;2060;1353;2357;3061;2798;3629;7212;6539;4874;4414;8519;8590;7555;6091;5300;5764;5492;9659;12840;13102;10638;9872;9726;8868;9591;10360;11243;10298;10987;7691;10100;9106;11058;12068;12982;12877;13204;12238;17949;14004;14481;14411;15769;15476;15728;16290;2567]; x4=[4233;4741;16387;16061;15092;3518;3469;3762;3927;4153;3267;4684;5495;5664;5541;5451;15087;16823;18303;18556;20582;14173;15517;14482;14318;10352;9095;10510;13954;10142;17765;6924;4678;6182;7653]; y4=[895;6434;6609;7352;6936;2571;2308;2170;2110;2299;793;1364;1205;1653;2093;2757;3512;4207;7385;5588;6548;11941;17034;12692;13569;17133;16414;15314;5615;1662;3561;5696;3765;2005;1952]; plot(x,y,＇*＇,x5,y5,＇^＇,x2,y2,＇+＇,x3,y3,＇p＇,x4,y4,＇x＇); xlabel(＇x-axis＇),ylabel(＇y-axis＇),zlabel(＇z-axis＇);