景区灭火的数学模型(共27页)

1、精选优质文档-倾情为你奉上西南交通大学峨眉校区2017年全国大学生数学建模竞赛第二次预选赛试题 A题：景区灭火的数学模型姓名刘琦王杨学号专业交运城轨电气工程联系电话QQ摘要本文以景区防火问题为背景，在通过多种方式插值方法修补残缺等高线的基础上，结合等高线高程信息建立景区三维地形模型并巧妙求出其地表总面积，最后通过图论模型制定最优灭火路线，达到及时控制和消灭火情的目的。问题一中，面对有5处需要修补的残缺等高线图，我们先根据各处残缺是否垂直于X轴的不同特点，通过人机交互分式将5处残缺线分离出来，并对原图进行去噪处理，进而将问题转换为一元插值问题。根据问题特点，我们选择了三次样条插值和三次多项式插值

2、，并对二者插值结果进行比较，发现三次多项式插值效果更好，因此我们选择三次多项式插值结果作为最终修补结果。最后我们将分区域修补好的等高线与去噪后的等高线重叠并保存，便可得到修补完全的等高线图。问题二中，我们在问题一得到完整等高线图的基础上，首先将8种高程信息不同的等高线从完整等高线种分离出来，并通过坐标变换，将地图上已知等高线的各点坐标转换为实际点坐标。这里考虑DTM（Digital Terrain Model）数字地面模型，由于高程值的集合是已知的，每一条等高线对应一个已知的高程值，可以得出其等高线模型，再使用MATLAB工具箱对多种不同插值方法得到的三维图进行比较，最后通过立方插值和最近点插

3、值的混合插值，得到了更多地表曲面的点再进行拟合处理，进而得到了较为理想的景区三维地形模型。对于景区地表面积的求解，我们可以通过对景区地表曲面进行曲面积分进行求解，但这种方式需要景区地表曲面表达式，在多项式拟合效果不理想的情况下，我们采用了三角网模型（TIN），以直代取，巧妙合理而又不失准确性的求出了景区地表面积。问题三我们将确定最优的灭火路线问题转换为求解景区地表曲面指定两点之间的最短曲线段问题。我们引入图论中赋权图的概念，将曲面上任意两点间曲线段的长度作为权值，赋给其对应投影所对应的网格的边，进而将求三维立体空间内最短曲线段的问题转换为求平面赋权图网络中指定顶点间具有最小权的路问题，采用迪克

4、斯特拉（Dijkstra）算法，按距A点从近到远的顺序，依次求得A点到图中各顶点的最短路和距离，直至B点。但计算量是远超计算机所能承受的，为了减少计算机计算量，我们采用了逐步优化的方法，对512*512个网格，先以8个网格为边长，64个小网格为一个大网格，进行求解，进而对每个大网格的路线进行二次优化，以减少误差。首先利用Excel处理数据，通过C+编写程序，最后得到了平面赋权图中A点到B点具有最小权的路，并利用MATLAB将其对应的曲线在三维地表曲面中表示了出来。 关键词： 等高线 DTM 图论 Dijkstra算法 三角网 C+ MATLAB专心-专注-专业一、问题提出1. 知识背景随着时代

5、和经济的发展，森林公园等自然景区逐渐成为了人们出游时一个不可或缺的选择，防火等安全问题也显得日益重要。本题基于景区防火的问题，要求在补全等高线的基础上建立其三维地形图，并且通过构建数学模型制定最优灭火路线，达到及时控制和消灭火情的目的。2. 需要解决的问题现旅游管理部门想在景区发生火灾时能及时控制和消灭火情，请你利用附件提供的数据通过建立数学模型解决下面三个问题：（1） 本题中给出的等高图存在局部破损，利用数学模型修补好该等高图。（2） 在修补好等高图的基础上，通过数学模型构建出景区的三维地形图，同时计算其地表面积。（3） 等高图中B点发生火灾，从固定的消防站A出发去灭火，通过建立模型确定最优

6、的灭火路线。附件：景区地形等高图说明：1.该图水平及竖直方向以10m每像素为单位，山高以50m为单位。2.实际图形见附件，为512×512像素。 二.基本假设1. 不考虑曲面坡度等其他因素对于救火员速度的影响，救火员在景区地表行走时速度恒定，即在空间内任意两点内救火员速度均保持不变。2.救火员接收到报警电话后就立即出发，景区表面每一处均不会阻碍救火员行进。3.计算面积时以50米海拔处为景区地平线。3. 符号说明符号意义单位备注网格的x坐标无全局变量网格的y坐标无全局变量网格的z坐标无全局变量S景区地表总面积全局变量G赋权图无全局变量V顶点集无全局变量E边的集合无全局变量W邻接矩阵无全

7、局变量到这条边是否在最短路径上无0-1矩阵结点到对应的空间曲线长度m全局变量结点到对应的空间曲线d对应的权值无全局变量D两个顶点之间具有的权无全局变量4 问题分析4.1问题背景随着我国的旅游行业进一步发展，选择观赏自然景区的人越来越多，在对一些旅游景区的消防安全以及灭火救援的要求上也有所提高。制定合理的灭火路线，及时控制和消灭火情，从而达到保障游客人身安全和保护自然景观的目的，直接影响了我国社会的安定和谐，同时对旅游经济的快速发展也具有重大意义。4.2分析问题在问题一中，将补全等高线转化为一元插值问题，利用插值法进行曲线修补，最终得到完整的等高图。 在问题二中，需要在第一问得到的完整等高线基础

8、上得到景区地形三维模型，我们可以采用二维插值的方式进行求解，对于表面积的求解，我们可以采用对曲面进行曲面积分的方法或采用三角网模型，以曲代直，进而得到景区地表表面积。在问题三中，需要我们制定最优的救火路线，我们可以转换为求解地表曲面上A点到B点的之间最短曲线段的问题，进而可以将曲面进行投影，引入赋权图概念，将最短曲线段的求解进一步转换为投影平面内求解对应两点的最小权路的问题。一元插值缺失的等高图二维插值问题一补全的等高图问题二景区三维地形图最短路径问题问题二景区地表表面积问题三的最优救火路径4.3数据的前期处理 由于5处缺口分别分布在4条等高线上，所以要对图像进行预处理，即将五处缺口分离为4张

9、新的图像。对于图像分割，我们可以采取灰度阈值分割法、基于区域的分割方法以及基于边缘的分割方法等等，但本题给出的图像较为简单，同时考虑到效率问题，所以采用人机交互的方法，利用画图板的橡皮擦手动分割图像。分割后可以得到以下4张图片。图4.1 图4.2 图4.3 图4.45 模型的建立与求解5.1问题一模型建立与求解5.1.1问题一的分析 由数据的前期处理，我们将五处缺口分别分散在四张图片中，图中只保留了缺口部分的黑点。然后利用MATLAB中的imread函数可以读出其灰度矩阵，则对等高线的补全就转化为了一元插值问题，利用插值法可以分别补全每个缺口，最后将得到的四张修补后的图片与最初的缺失等高图进行

10、图层重叠，重叠后可得到完整的等高图。5.1.2问题一模型的建立 由问题知我们需要读出每张图片的灰度矩阵，然后选择当中部分点进行插值法，首先采用人为选点的方法，但在同一图片选取不同组的点进行插值，对结果进行比较后发现，人为选点的方法具有很大的主观性，会产生较大的误差，所以该方法不可取。于是利用图片中保留的所有点，即插值的对象为整个矩阵。我们分别采用对线性插值、三次多项式插值、三次样条插值等多种方法进行插值运算，对曲线进行修补，比较其插值结果，发现三次多项式插值的结果更加合理，所以本问题全部采用三次多项式插值的方法。 5.13问题一的模型求解 对等高线的补全可以转化为一元插值问题，利用MATLAB

11、中imread函数读出每张图片的灰度矩阵，选取合适的插值方法即可解决。1. 计算思路：经过插值结果的比较，发现三次多项式插值更加适合本题，用MATLAB中的imread函数读出每张图片的灰度矩阵，然后运用三次多项式插值，可以补全图片中的缺失部分，其余剩下的图片采取同样的方法可以补全缺失部分，最后将处理后的图片与原等高图在同一窗口中显示，即得补全的等高图。2. 计算步骤：Step1：利用imread函数读出图1的灰度矩阵；Step2：图4.1缺失部分可以由灰度矩阵得到的坐标点数据进行三次多项式插值得到，图一的补全完成；Step3：对图4.2的操作重复前面两步即可完成；Step4：图4.3、图4.

12、4首先重复第一步，但由于图中缺失部分附近的等高线弯曲程度太大，导致灰度矩阵在读取时无法一一对应，所以在读出其灰度矩阵后要先对矩阵进行转置，然后重复第二步，图三图四的补全完成；Step5：图4.1-4.4的补全均完成，将其与原等高图在同一窗口显示，得到最终补全的等高图。3. 计算结果：经过上述的处理后，将最终补全得到的等高图保存为bmp，单色格式，图片如下（程序见附录）： 图5.1.1 三次样条插值结果 图5.1.2三次多项式插值结果（效果更好） 由多种方式插值得到的图像可知，相比于三次样条插值，三次多项式插值得到的完整等高线更为光滑，效果也更好。因此最终我们选择三次多项式插值的结果作为等高线的

13、补全图。5.1.3问题一结果的分析与验证对于补全后的等高图，我们可以采用分析点的差异度来判断曲线的补全是否合理。对于点的分析包括该点的位置和其对应的切向量，然后对两个点进行切向量和差异函数的计算，差异函数值越低则表示图像中的两个点越接近事实。5.2.1问题二模型的分析问题要求我们在问题一的基础上，根据补全的等高线图建立景区地表的三维地形图，并求出景区地表面积。通过第一问对于等高线的修补，结合等高线的高程信息，我们可以得到景区地表的离散坐标，即我们可以用一组有序数值阵列形式表示地面高程信息。但我们所得到的点比较稀疏，数量还远远不能拟合出地表曲面。因此，我们需要进行二维插值，进行二维插值模型。根据

14、已有点信息插值得到更多地表曲面上的点，然后MATLAB进行拟合曲面。得到拟合后的地表曲面后，意味着我们得到了地表曲面上所有点的坐标信息以及曲面的表达式，紧接着我们便建立三角形网模型，采用以直代取，将曲面面积用足够多的小三角形平面面积去表示，最后得到的三角形面积总和即为景区地表曲面面积。而我们欲得到景区地表三维模型，首先应根据像素与比例尺求得地图中各等高线的实际高度，我们首先将修补好的等高线图根据各等高线高程不同进行分割，进而采用MATLAB循环赋值的方法进行数据的前期处理（程序见附录）。等高线1-8如下图所示： 5.2.2问题二模型建立首先我们先采用MATLAB内部二维插值工具Curve Fi

15、tting Tool中的多种插值方式Cubic、v4，以及可以求出曲面表达式的polynomial等插值方法进行求解，求得的三维地形效果差不多，但是根据可以求出曲面表达式的polynoimal插值方法得到的曲面效果却不理想，因为polynoimal插值法x,y最高次为二元五次表示式，故我们用Custom Equation自拟表达式来提高x,y次数，发现随着次数的增高，效果并没有变得更理想。 图5.2.1 Cubic插值 图5.2.2 V4插值图5.2.3 二元五次表达式拟合曲面 图5.2.4 二元六次表达式拟合曲面 因此我们需要重新选取方法，考虑到求解表面积需要插值点坐标以及本问题特点，我们根

16、据部分等高线高程进行二维插值，而二维插值分为结点插值与散乱结点插值，因本问题中曲面上的点并不是全部落在网络节点上，因此我们采用散乱结点插值，即本问题可简化为： 已知个结点，求点处的插值。 对于上述问题。MATLAB中提供了插值函数gridata，其格式为：其中，x,y,z均为n维向量，指明所给数据点的横坐标，纵坐标和竖坐标；向量是给定的网格点的横坐标和纵坐标；返回值为网格处的函数值。为方向不同的向量，即一个是行向量，另一个是列向量。为保证准确性以及减小误差，我们采用混合插值的方式，进行立方差值和最近点插值的混合插值，通过编写程序可得到三维地形图（程序见附录）：我们首先得到了以50米为“地平线”

17、的景区地形三维图，即得到了海拔高度50米以上的景区地形图。图5.2.5 50米为地平线的景区三维图但题目要求我们建立整个景区的地形三维模型。对于海拔小于50米的部分，我们首先想到的是将其他部分都假设海拔高度为0，但这种想法会产生大量断崖，显然不符合实际。而实际情况应该为其余部分海拔高度逐渐上升，进而与50米处平缓连接。因此，我们考虑了其余区域，将其也列入插值区域，重新进行插值，并得到了重新插值后的三维地形图及其俯视图：图5.2.6 景区地形三维图 图5.2.7 景区地形三维图（俯视角） 由最终插值所得到的景区地形三维地形图及其俯视图可知，最终插值将海拔50米以下的区域考虑进去，弥补了第一次插值

18、的缺陷，建立了（）的整个景区的地形三维图。成功建立了景区的地形三维图后，我们利用已建立的模型以及数据进行求解景区地表表面积，我们首先想到的方法就是通过多项式拟合的方式求出曲面的表达式，进而利用曲面积分便可求得地表曲面的面积，即：便可求得地表曲面的面积。但是从表达式进行拟合曲面效果图来看，可以发现拟合效果并不理想，误差太大，因此我们用表达式进行曲面积分而求出景区表面积的方法是不可取的，所以我们建立三角形网模型来近似求得景区地表表面积。在本模型中我们先将景区地表划分为n个小三角形，以直代曲，求景区面积就近似等于这n个小三角形的面积之和。因为三点可以确定一个平面，而且等高图为512*512像素的图像

19、，容易想到将整个景区表面积划分为512*512个网格，选取其中一个网格示意图如下： 其中A,B,C,D分别代表上述网格的四个顶点。Matlab中的cross函数，用于返回两个向量的数量积，为了方便编程，所以采用向量的数量积来对三角形的面积进行求解，则面积为同理，可以求得面积为：那么一个网格的面积为最终得到景区地表面积的为：由此，对所有的小三角形进行遍历求和，便可近到整个景区地表面积S 5.2.3问题二模型的求解1. 计算思路：需要在第一问得到的完整等高线基础上得到景区地形三维模型，我们可以采用二维插值的方式进行求解，对于表面积的求解，我们可以采用对曲面进行曲面积分的方法或采用三角网模型，以曲代

20、直，进而得到景区地表表面积。2. 计算步骤：Step1：将第一问得到的完整等高线按照高程分离，得到8张分离后的等高线。Step2：利用混合插值的方法还原其三维地形图。Step3：将景区表面积划分三角网，以直代曲，通过matlab编程对所有三角形面积求和最终得到景区地表表面积。3. 计算结果： 由我们在之前的到的二维插值结果，编写遍历循环程序（程序见附录）对所有小三角形面积进行求和。最后得到整个景区地表表面积近似为： 5.2.4问题二结果的分析与验证1.三维地形图的结果分析：在构建景区的三维地形图时，采用了多种插值方法对图形进行绘制，采用MATLAB内部二维插值工具Curve Fitting T

21、ool中的多种插值方式Cubic、v4，以及可以求出曲面表达式的polynomial等插值方法进行求解，求得的三维地形大致符合要求，但是在光滑度等方面较差。通过比较我们采用混合插值，将多种方法结合，最终得到满意的三维地形图。2.景区地表表面积的结果分析：在问题二中计算面积时，我们并没有对等高线区域以外的部分进行条件限制，若假设等高线区域以外的点高度都是50米，利用相同的算法，计算可得其面积为将两种处理方式进行对比，计算其差异度由此可见，两种处理计算计算出的表面积的差异度很小。假设该景区全为平地，那么其表面积，因为构建出的景区地形三维图存在坡面，所以其地表表面积一定大于全为平地情况下的表面积，通

22、过比较：。结合上述两种方法的分析验证，我们可以确定计算出的地表表面积是正确的。5.3.1问题三问题分析 问题三在问题二的基础上，要求建立一条救火的最佳路径。我们假设救火员在空间任意两点之间的速度是相同的，因此，我们可以将问题转换为最短路径的问题。因为此救火路径为空间中位于地表曲面上的一条曲线，我们可先将位于固定曲面上两点的最短曲线转换为二维顶点之间最短路问题的数学规划模型。我们引入赋权网络图的概念，通过景区地表曲面在平面内的投影区域构造赋权有向图。5.3.2问题三模型建立构造赋权图:（1）以景区地表曲面在平面内的投影区域上划分网格中的结点为网络图的结点，记做。顶点集，这里表示各个结点。（2）以

23、景区地表曲面在平面内的投影区域上划分网格中的结点为网络图的边。记做。为边的集合，（3）以各条边在景区地表曲面上对应的曲线段的长度作为其相应边的权值，即为网络图中结点到的权值。我们设投影到顶点的空间坐标为，投影到顶点的空间坐标为，根据假设，我们近似以曲面上两点之间的线段长度来代替两点之间的曲线段长度，即：（4） 邻接矩阵是表示顶点之间相邻关系的矩阵，以邻边矩阵表示网络图，记做，因本模型建立的为赋权图，有： 我们设救火员起点A为,终点B为,我们引入0-1变量：我们的目标函数是求赋D权图上指定的两个顶点之间具有最小权的路。即最短路径上各条边的长度权值最小，即： 故最终我们建立的数学规划表达式为： 5

24、.3.3问题三模型的求解1.计算思路：问题三在问题二的基础上，要求建立一条救火的最佳路径。假设在空间任意两点之间的速度是相同的，因此将问题转换为最短路径的问题。因为此救火路径为空间中位于地表曲面上的一条曲线，我们可先将位于固定曲面上两点的最短曲线转换为二维顶点之间最短路问题的数学规划模型。我们引入赋权网络图的概念，通过景区地表曲面在平面内的投影区域构造赋权有向图。确定起点A和终点B利用Dijkstra算法确定行走路径。迪克斯拉特（Dijkstra）算法，作为求解最短路径一种成熟的算法，其基本思想是按距离从近到远为顺序，依次求得到G的各项点的最短路和距离，直至，算法结束。为避免重复并保留每一步的

25、计算信息，采用了标号算法。下面为该算法：Setp1:令=0,对，令。Setp2:对每个，用代替，这里用来表示顶点之间边的权值。计算，把达到这个最小值的一个顶点记为，令。Setp3:若，则停止；若，则用代替，转Setp2。算法结束时，从到各顶点的距离由的最后一次标号给出。在进入之前的标号叫T标号，进入之后的标号叫P标号。算法就是不断修改各项顶点的T标号，直至获得P标号。2.计算步骤： 问题三我们采用Dijkstra算法，但由于数据量过大，如果对于全局采用Dijkstra算法，运算量会很大，导致算法效率不高。因此，我们采用逐步优化的搜索方法，从而达到减少运算量的目的。逐步优化的步骤如下：Step1

26、：输入景区地形的数值模型，降低网格密度，取80m为单位网格长度Step2；在step1得出的稀疏网格中，利用Dijkstra算法求出最优路径，将此路径上的点作为控制点，转Step3。Step3：利用step2中的控制点，将每个80m*80m的控制节点分为64个10m*10m的网格，局部优化，即二次优化。Step4：将每个控制节点与其局部优化的节点结合，最终得出最优路径。3.计算结果：我们首先利用C+软件编程我们得到了在稀疏网格中节点的距离矩阵（程序见附录），再通过MATLAB软件编程我们得到了最佳救火路线在地形三维图上的表示，如下图所示（程序见附录）： 5.3.4问题三模型的改进在求解问题三时

27、，我们不考虑曲面坡度等其他因素对于救生员速度的影响，假定救火员在景区地表行走时速度恒定。但在实际情况中，由于受到坡度、地面平坦度和树木分布等因素，这种假设是不合理的。所以问题三转化为一个多目标决策问题，我们根据层次分析法的原理，可以建立筛选最优路线的层次结构模型。模型兼顾距离、坡度、地面平坦度、树木分布、火势蔓延等因素，确定各方案的优劣顺序，进行最优路线分析，构建AHP判断矩阵并计算权重，最终得到最优路线。结果表明，改进后的模型对于救火路线方案进行优选更加科学、方便。6、 模型的评价与推广6.1 模型的评价优点：1. 问题一、问题二绘制图形时采用了多种插值方法进行比较，最终得出最理的结果。2.

28、 问题三的模型简单易懂，并且考虑到效率问题，采取了二次优化，大大减少了计算量。缺点：1. 在问题三中没有考虑坡度等因素，在实际情况中其实是不允许的。2. 由于给出的等高图较为简单，所以在处理图像时很多都是人工进行操作，对于那些等高线密集的图像，这种处理方式是很不方便的。6.2 模型的推广 对于问题一中使用的三次多项式插值，通过利用一些特征节点，用多项式拟合的方法来产生平滑的插值曲线，其易操作且计算量不大，对于温度、高程、地下水位高度和污染浓度等问题同样适用。对于问题二，在使用插值的基础上，利用软件重建其三维地形图，由平面图得到其三维图，在机械三维、建筑三维、美术三维等方面都是极其重要的。问题三

29、对于生活生产中最短路径的问题中具有参考意义。7、 参考文献1数学建模算法与应用 司守奎，孙玺菁编著 国防工业出版社 2016.72MATLAB在数学建模中的应用（第二版）卓金武著 北京航空航天大学出版社 2014.93MATLAB教程 张志涌，杨祖樱等编著 北京航空航天大学出版社 2015.14基于等高线的三角形建模及真实感地形重建 翁巧琳,姜昱明 西安电子科技大学计算机学院，陕西 西安 5基于MSDAE的等髙线地图断线重连技术研究汪梁 西安电子科技大学计算机学院 TP3116DEM数据提取关键技术研究 乔利军 国防科技大学研究生院 2005.117彩色地形图中等高线断裂重建算法的研究 李晓清

30、 西安电子科技大学硕士学位论文 2014.12八、 附录附录清单：附录1：求解问题一的MATLAB程序附录2：求解问题二的MATLAB程序附录3：求解问题三的C+程序附录1 问题一中补全等高线MATLAB程序8.1.1三次样条插值程序%3图4图 无需转置 a=ones(512,512);imshow(a);hold onfor n=3:4imageName=strcat(num2str(n),'.bmp');A=imread(imageName);k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1;break

31、endendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'spline');plot(xi,yi,'k')axis(1,512,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test1.bmp')c=imread('test1.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test3.bmp')%1图2图 要转置a=on

32、es(512,512);imshow(a);hold onfor n=1:2imageName=strcat(num2str(n),'.bmp');A=imread(imageName);A=A'k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1;breakendendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'spline');plot(xi,yi,'k')axis(1,51

33、2,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test2.bmp')c=imread('test2.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test4.bmp')%3张图合并到一张图里并保存a=imread('quzaoyuantu.bmp')b=imread('test3.bmp')c=imread('test4.bmp')c=c'for i=1:512for j=1:512

34、if a(i,j)&b(i,j)&c(i,j)=1d(i,j)=1;elsed(i,j)=0;endendenddimshow(d)imwrite(d,'zuizhongtu.bmp')8.1.2三次多次式插值程序%3图4图 无需转置 a=ones(512,512);imshow(a);hold onfor n=3:4imageName=strcat(num2str(n),'.bmp');A=imread(imageName);k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1

35、;breakendendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'cubic');plot(xi,yi,'k')axis(1,512,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test1.bmp')c=imread('test1.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test3.bmp')%1图2图 要转

36、置a=ones(512,512);imshow(a);hold onfor n=1:2imageName=strcat(num2str(n),'.bmp');A=imread(imageName);A=A'k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1;breakendendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'cubic');plot(xi,yi,'k')axis(

37、1,512,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test2.bmp')c=imread('test2.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test4.bmp')%3张图合并到一张图里并保存a=imread('quzaoyuantu.bmp')b=imread('test3.bmp')c=imread('test4.bmp')c=c'for i=1:512for j=1

38、:512if a(i,j)&b(i,j)&c(i,j)=1d(i,j)=1;elsed(i,j)=0;endendenddimshow(d)imwrite(d,'zuizhongtu.bmp') 附录2 问题二中建立三维模型及求得地表表面积MATLAB程序clear C=; a=0; z=; for n=1:8imageName=strcat(num2str(n),'.bmp');A=imread(imageName);k=1;for i=1:512for j=1:512if A(i,j)=0B(1,k)=i;B(2,k)=j;k=k+1;end

39、endendB;C=C B;B=;x=C(1,:);y=C(2,:);D=size(C);b=D(:,2)-a;a=D(:,2);E=50*n*ones(1,b);z=z E;endxmm=minmax(x)ymm=minmax(y)xi=xmm(1):xmm(2);yi=ymm(1):ymm(2);zi1=griddata(x,y,z,xi,yi','cubic');zi2=griddata(x,y,z,xi,yi','nearest');zi=zi1;zi(isnan(zi1)=zi2(isnan(zi1)subplot(1,2,1),plo

40、t(x,y,'*')subplot(1,2,2),mesh(xi,yi,zi)s=0;for i=1:403 for j=1:367 a=i*10,j*10,zi(i,j); b=i*10,j*10+10,zi(i,j+1); c=i*10+10,j*10,zi(i+1,j); d=i*10+10,j*10+10,zi(i+1,j+1); s1=0.5*norm(cross(b-a,c-a); s2=0.5*norm(cross(b-d,c-d); s=s+s1+s2; endendsst=5120*5120-3680*4040+s求解结果：s = 1.5174e+07st =

41、 2.6521e+07附录3 求解问题三中最短路径的MATLAB程序8.3.1求解距离矩阵的C+程序：#include<iostream>#include<math.h>#include<iomanip>using namespace std;void main(void)int qd,zd;long double qz,qz1,qz2,qz3,qz4,qz5,qz6,qz7,qz8,qz9,qz10;int i,n,m;cout<<"请输入长度为10边的个数n&qu

42、ot;cin>>n;for(i=0;i<n;i+)        cin>>qdi>>zdi>>qz1i>>qz2i;qzi=pow(double(10*10+pow(abs(qz1i-qz2i),2),(double)1/2)+pow(double(10*10+pow(abs(qz2i-qz3i),2),(double)1/2)+pow(double(10*10+pow(abs(qz3i-qz4i),2),(double)1/2)+pow(doub

43、le(10*10+pow(abs(qz4i-qz5i),2),(double)1/2)+pow(double(10*10+pow(abs(qz5i-qz6i),2),(double)1/2)+pow(double(10*10+pow(abs(qz6i-qz7i),2),(double)1/2)+pow(double(10*10+pow(abs(qz7i-qz8i),2),(double)1/2)+pow(double(10*10+pow(abs(qz8i-qz9i),2),(double)1/2)+pow(double(10*10+pow(abs(qz9i-qz10i),2),(double)

44、1/2);cout<<"nnnnnnnnnnnnnnn计算完成nnnnnnnnnnnnnnn"for(i=0;i<n;i+)cout<<setw(10)<<qdi<<setw(10)<<zdi<<setw(10)<<qzi<<endl; 8.3.2求解救火路线的C+程序：#include <iostream>#include <limits>using namespace std;struct Nod

45、e  /定义表结点  int adjvex; /该边所指向的顶点的位置  int weight;/ 边的权值  Node *next; /下一条边的指针;struct HeadNode / 定义头结点    int nodeName; / 顶点信息    int inDegree; / 入

46、度    int d; /表示当前情况下起始顶点至该顶点的最短路径,初始化为无穷大    bool isKnown; /表示起始顶点至该顶点的最短路径是否已知,true表示已知，false表示未知    int parent; /表示最短路径的上一个顶点    Node *link; /指向第一条依附该顶点的边的指针;/G表示指向头结点数组的第一个结点的指针

47、/deNum表示结点个数/arcNum表示边的个数void createGraph(HeadNode *G, int nodeNum, int arcNum)   cout << "开始创建图(" << nodeNum << ", " << arcNum  << ")"

48、 << endl;  /初始化头结点  for (int i = 0; i < nodeNum; i+)     Gi.nodeName = i+1; /位置0上面存储的是结点v1,依次类推    Gi.inDegree = 0; /入度为0    G

49、i.link = NULL;    for (int j = 0; j < arcNum; j+)     int begin, end, weight;    cin >> begin >> end >> weight; &#

50、160;  / 创建新的结点插入链接表    Node *node = new Node;    node->adjvex = end - 1;    node->weight = weight;    +Gend-1.inDegree; /入度加1   &

51、#160;/插入链接表的第一个位置    node->next = Gbegin-1.link;    Gbegin-1.link = node;  void printGraph(HeadNode *G, int nodeNum)   for (int i = 0; i < nodeNum; i

52、+)     cout << "结点v" << Gi.nodeName << "的入度为"    cout << Gi.inDegree << ", 以它为起始顶点的边为: "    Node *node =&

53、#160;Gi.link;    while (node != NULL)       cout << "v" << Gnode->adjvex.nodeName << "(权:" << node->weight << ")"&

54、#160;<< "  "      node = node->next;        cout << endl;  /得到begin->end权重int getWeight(HeadNode *G, int begin, int end) &

55、#160; Node *node = Gbegin-1.link;  while (node)     if (node->adjvex = end - 1)       return node->weight;        node =

56、0;node->next;  /从start开始，计算其到每一个顶点的最短路径void Dijkstra(HeadNode *G, int nodeNum, int start)   /初始化所有结点  for (int i = 0; i < nodeNum; i+)     Gi.d = INT_MAX;

57、 /到每一个顶点的距离初始化为无穷大    Gi.isKnown = false; / 到每一个顶点的距离为未知数    Gstart-1.d = 0; /到其本身的距离为0  Gstart-1.parent = -1; /表示该结点是起始结点  while(true)    /= 如果所有的结点的最短距离都已知,&#

58、160;那么就跳出循环   int k;   bool ok = true; /表示是否全部ok   for (k = 0; k < nodeNum; k+)      /只要有一个顶点的最短路径未知,ok就设置为false     if (!Gk.isKnown

59、)        ok = false;       break;            if (ok) return;   /=   /= 搜索未知结点中d最小的,将其变为known   

60、/= 这里其实可以用最小堆来实现   int i;   int minIndex = -1;   for (i = 0; i < nodeNum; i+)      if (!Gi.isKnown)        if (m

61、inIndex = -1)         minIndex = i;       else if (GminIndex.d > Gi.d)         minIndex = i;     

62、      /=   cout << "当前选中的结点为: v" << (minIndex+1) << endl;     GminIndex.isKnown = true; /将其加入最短路径已知的顶点集     / 将以minIndex为起

63、始顶点的所有的d更新     Node *node = GminIndex.link;     while (node != NULL)        int begin = minIndex + 1;       int end

64、 = node->adjvex + 1;       int weight = getWeight(G, begin, end);       if (GminIndex.d + weight < Gend-1.d)