景区灭火数学建模模型.docx

西南交通大学2012年大学生数学建模竞赛题目： B题 参赛队员1参赛队员2参赛队员3姓名甘罗谈盼明浩学号201038742010388720103878学院数学学院数学学院数学学院专业数学与应用数学数学与应用数学数学与应用数学电话134028869451878294300813688127427Email1522407022qq.com576946951qq.com173026319西南交通大学教务处西南交通大学实验室及设备管理处西南交通大学数学建模创新实践基地题目： 景区灭火的数学模型 摘要本文借助matlab的图像处理功能，运用插值拟合，有限元的思想，建立了消防员救火最佳路线和消防站合理选址的优化模型，再利用lingo软件求得最优解，由最优解确定救火最佳路线和消防站合理选址。附件所给等高线图为像素为512512的二维灰度图.利用matlab的imread函数读出等高线图的灰度矩阵.根据灰度矩阵求得灰度图中黑色点的坐标；补全等高线图转化为一元插值问题;通过对线性插值、三体样条插值、三次多项式插值结果进行分析比较，三次多项式插值效果更为理想。由问题一得出的完整等高线图可求出等高线上对应山体表面上点的三维坐标，利用matlab中的Surface Fitting Tool拟合得到景区的三维地形图和拟合后的曲面函数。用matlab算出山体表面的曲面积分，由像素与比例尺的值可得山体表面面积。从消防站A派遣消防员去B点灭火，最佳的灭火路线应使得消防员到达火灾蔓延区域边界的时间最短。将火势蔓延路线与消防员灭火路线合并为一条曲线，在曲线上均等插入n-1个点，得到的折线可以近似代替该曲线。折线上n-1各点的位置可以由n-1参变量唯一确定，则最短时间可由n-1个参变量唯一确定。求最佳的灭火路线转化为转化为以最短时间为目标函数，n-1个参变量的取值范围为约束条件的优化问题。利用lingo软件可求得最优解。根据实际情形，消防站选址的合理性由多方面因素决定。本文定义消防站的合理选址应使得火灾发生对森林公园造成损失的期望最小。消防员到达着火蔓延区域边界的时间由第三问的模型可得。再由假设4、5，可建立以火灾造成损失的期望最小为目标函数的优化模型。考虑到运算量过大，我们将景区分成了16个片区，取各个片区的中心代表该片区，近似估计最优选址点的位置。关键词：灰度矩阵 一元插值 曲面拟合 有限元法优化模型一问题提出某国家级森林公园由于风景区植被丰富，拥有大量的国家级重点保护动植物，因此旅游管理部门在风景区的A点设置了景区消防站，当景区发生火灾时能及时控制和消灭火情。图1说明：该图水平及竖直方向以10m为单位，山高以50m为单位。根据所给出的风景区地形等高线图形提出了以下问题：1、由于人为原因，图1所示的等高图出现了局部破损的情况，根据所学知识利用适当的方法建立数学模型修补好该等高图；2、在完成第一问的基础上，结合数学模型建立该景区的三维地形图，并估计该景区的地表面积；3、某天图1所示的B点发生了火灾，于是需要从景区消防站派遣消防员去B点灭火，建立模型确定最佳灭火路线；4、 如果需要对景区消防站进行重新选址，请建立模型确定合理的消防站地址。二问题分析2.1问题1的分析：根据所学的数学知识可以采用曲线插值的方法将等高图中缺损的部分补充完全；2.2问题2的分析：要计算景区的地表面积，可以采用曲面积分的方法。这就要知道地表曲面的曲面方程或函数，而曲面的函数z=fx,y可以根据所给出的等高线图上已知的点用matlab拟合出来，那么地表面积就可以用1+zy2+zy2dxdy求出；2.3问题3的分析：要确定最佳的灭火路线即确定火势蔓延到某处与消防员到达这里的最短时间，而最短时间取决于火势蔓延的速度和消防员行走的速度，那么就需要求出A、B两点间的路线，设路线为L，则可以将问题转化为优化问题，再利用有限元的思想就可以求解出最佳路线；2.4问题4的分析：要使景区消防站更合理，义消防站的合理选址应使得火灾发生对森林公园造成损失的期望最小。，那么可以将景区合理划分成各个小片区，再求每个小片区的中心点建消防站对应的损失期望，取该值对应的片区的中心即为最佳选址。三问题假设1.假设消防员前往火灾发生点的速度和火蔓延的速度只与坡度有关2.消防员在山体表面任意一点都可通行3.消防站在着火时马上收到警报，同时消防员赶往火灾地点4.景区山体表面上每一点发生火灾的概率相等5.景区山体表面每一点发生火灾造成的损失只消防员到达火灾区的时间正相关四符号说明s:山的表面积A: 消防站B:火灾发生点A0: A在水平面上的投影B0: B在水平面上的投影t:消防员到达火势蔓延区域边界的时间Pi:以A、B为端点的折线的n-1个折点（i=1,2,3n-1，P0、Pn分别表示B、A）P0i：Pi在水平面上的投影Mi: 线段A0 B0的第i个等分点（i=1,2,3n-1，M0、Mn分别表示B0、A0）：过B0和P0i的线段与x轴的夹角 ：PiPi+1与z轴的夹角v1：消防员速度v2：火势蔓延的速度E(Y)：以Y点为消防站时火灾发生造成损失的期望五模型建立与求解5.1 问题一的模型建立与求解附件所给等高线图为像素为512512的二维灰度图.利用matlab的imread函数可以读出灰度图灰度矩阵.根据灰度矩阵可以求得灰度图中黑色点的坐标.那么等高线图的补全就转化为一元插值问题;对所给等高线图进行处理,得到只包含缺失部分附近黑色点的18张图片;将灰度矩阵的行与列分别对应坐标系中的纵坐标与横坐标.当同一横坐标上有多个数据点时,只保留其中一个.观察可知,第1-16张图片等高线缺失部分可直接由灰度矩阵得到的坐标点数据直接插值得到.对线性插值，三次抽样插值，三次多项式插值结果比较，发现三次多项式插值更为合理。通过程序1.1将第1-16张图片插值得到后的图形显示在同一窗口中,如图1.1所示图1.2第17-18张图片等高线缺失部分与水平方向垂直；可先对灰度矩阵转置对应的灰度图插值拟合，程序1.2运行结果结果如图1.2:图1.2将原等高线，图1.1，图1.2对应灰度矩阵的转置的灰度图在同一图像窗口中显示，并将其存为bmp、单色格式。即得补全的等高线图。如图1.3图1.35.2 问题二的模型建立与求解利用5.1的模型得到的等高线图结果，首先将图片中完整的等高线分离，如图2.1，图2.1说明：左边为等高线1，右边为等高线3再用matlab软件中的imread函数将上述8张分离的等高线图读出灰度矩阵，然后将8个矩阵中的等高线上的点的位置取出并合成一个矩阵C，再取得等高线上的点的横坐标x、纵坐标y，同时计算出各个等高线上点的竖坐标z，具体matlab程序如附件中程序2.1所示，得到x,y,z后，利用matlab中的Surface Fitting Tool将景区的三维空间曲面函数拟合出来，则可以得到景区的空间曲面图，如图2.2所示，图2.2同时Surface Fitting Tool得到拟合的多项式函数为p00=-553 ;p10 =7.007;p01 = -0.2078 ;p20 = -0.03322 ;p11 = 0.02264;p02 = 0.03541;p30 = 8.848e-005 ;p21 = -0.0002331 ;p12 = -0.0001663 ;p03 = -0.0001483 ;p40 = -8.863e-008 ;p31 = 5.859e-007 ;p22 = 3.941e-007 ;p13 = 4.976e-007 ;p04 = 1.904e-007;p50 = -8.63e-011 ;p41 = -6.388e-011 ;p32 = -9.059e-010 ;p23 = -1.191e-010 ;p14 = -4.111e-010 ;p05 = -8.165e-011 ;z = p00 + p10*x + p01*y + p20*x2 + p11*x*y + p02*y2 + p30*x3 + p21*x2*y+ p12*x*y2 + p03*y3 + p40*x4 + p31*x3*y + p22*x2*y2 + p13*x*y3 + p04*y4 + p50*x5 + p41*x4*y + p32*x3*y2 + p23*x2*y3 + p14*x*y4 + p05*y5;利用曲面积分公式建立表面积模型s=100*1+zy2+zy2dxdy用matlab编程即可将上式面积求出，程序如附件中程序2.2，算得结果为s=26112100。5.3 问题三的模型建立与求解从景区消防站A派遣消防员去B点灭火，最佳的灭火路线应使得消防员到达火灾蔓延区域的时间最短。如果消防员灭火路线与火势蔓延经过的某条路线相交所需的时间最短，那么该路线即为消防员最佳的灭火路线。消防员灭火路线与火势蔓延经过的曲线合为一条曲线。在曲线上插入n-1个点，将相邻点连接，得到一条折线，当n-1很大时，可用该折线代替曲线；以消防站B在水平面上的投影B0为原点，B0 B为z轴建立右手空间直角坐标系，连接A0与B0两点，将线段A0B0均等插入n个点P0i，分别以线段Mi为半径、B0 B为中轴线做n个圆柱面，与山体表面相交曲线，依次在曲线上选取折点Pi（i=1,2,3n-1，P0、Pn分别表示B、A）, Pi在水平面上的投影点为P0i；B0P0i与x轴正向夹角为；Pi在水平面上的投影P0i与B在水平面上的投影B0点的距离| B0P0i|=P0i坐标为：(| B0P0i|,| B0P0i|,0)Pi坐标(xi,yi.,zi)为：(,f(,)相邻两折点之间的距离|PiPi+1|=PiPi+1与z轴的夹角=消防员经过线段PiPi+1所需时间 t1i=火势沿线段PiPi+1蔓延所需时间 t2i= 在、确定的情况下t1i、 t2i可以表示为（i=0,1,2,3n-1的函数，当n取定值时，t1i、 t2i是关于（i=0,1,2,3n-1）的函数设消防员在折线段PmPm+1到达火势蔓延区域则有可得m=m（,.）消防员t时刻在折线段PmPm+1到达火势蔓延区域，关于t有当n取定值时，t可以表示为,.的函数，将min（t(,.)）作为目标函数，02作为约束条件求最优解的问题；目标函数：min（t）= min（t(,.)）约束条件：1.2.3.t1i=，t2i= 4.02在、确定的情况下t1i、t2i可以表示为,.的函数，t=t(,.在02上求最优解，利用lingo软件，可解出,.，,.唯一确定最佳路线，对应t为消防员到达火势蔓延区域的最短时间。5.4 问题四的模型建立与求解森林公园中某点发生火灾将造成植被破坏，而植被破坏的面积与火持续烧的时间成正相关，火持续的时间主要与消防员到达的时间有关；火灾造成的损失主要与植被破坏的面积有关。无论是从经济还是生态考虑，合理的消防站选址应该使得火灾发生造成的损失期望最小。根据假设4、5，火灾发生造成的损失期望最小等价于消防员到达火势蔓延区域边界的时间t的期望最小。假设点Y：(x，y)为消防站，Y：（x,y）为山上任意一点；消防站为Y,Y发生火灾时对应的t可由问题三的模型求得，记为t（x,y）消防站设为Y时对应的t的期望：E(Y)=求出山上各点的E(Y)，min（E(Y)）即为火灾发生损失的期望的最小值；考虑到t（x,y）运算量过大，我们考虑将景区分成片区，取各个片区的中心代表该片区。例如我们可以等高线图上的比例尺以128128为一个单位，将山体表面划分为16个片区，取各个片区的中心代表该片区。求出E(Y)=求出山上各点的E(Y)，min（E(Y)）即为火灾发生损失的期望的最小值；min（E(Y)）对应的（x,y）即为最佳选址地点。六评价与推广一模型的优缺点优点1.模型建立过程逻辑性强、通俗易懂2.运用matlab的进行图像处理，简洁直观3.问题三的模型化难为简，利用有限元的思想以直代曲，将积分问题转化为离散型优化问题4.模型四从实际出发，对选址合理性的定义符合实际缺点1. matlab，lingo程序运行效率不高二模型的推广1.模型一、模型二可用于分析统计数据，对位置样本给出估计2.模型三、四贴近现实，对工厂，医院等公共机构选址具有指导作用八参考文献1 姜启源 谢金星 叶俊 数学模型 北京 高等教育出版社 2011.12韩忠庚 数学建模方法及其应用 北京 高等教育出版社 20093 薛长虹于凯大学数学实验MATLAB应用篇成都西南交通大学出版社2003附录问题一附录只包含缺失部分附近黑色点的18张图片1-16：17-18：程序1.1,对第1-16张图片插值拟合的程序clearclca=ones(512,512);imshow(a);hold onfor m=1:16url=strcat('C:MATLAB7work',num2str(m); imageurl=strcat(url,'.bmp'); A=imread(imageurl) ;k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i; k=k+1;break;endendend B; 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')程序1.2,用matlab读取第17-18张图片对应的灰度矩阵转置的灰度图插值拟合clearclca=ones(512,512);imshow(a);hold onfor m=17:18 url=strcat('C:MATLAB7work',num2str(m); imageurl=strcat(url,'.bmp'); A=imread(imageurl) ;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;break;endendend B; 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,'test2.bmp')c=imread('test2.bmp')程序1.3 得到补全的等高线图的程序clearclca=imread('原图.bmp')b=imread('test1.bmp')c=imread('test2.bmp')c=c'for i=1:512for j=1:512 if a(i,j)&b(i,j)&c(i,j)=1d(i,j)=1;elsed(i,j)=0;endendenddimshow(d)F1=getframe;imwrite(F1.cdata,'test3.bmp')c=imread('test3.bmp')问题二附件:八张海拔高度相等的等高线图:程序2.1：由等高线图包含的数据二元插值拟合山体曲面Clc;Clear; C=;a=0;z=;for s=1:8url=strcat(' C:MATLAB7worksecond',num2str(s);imageurl=strcat(url,'.bmp');A=imread(imageurl) ;k=1;for i=1:512for j=1:512if A(i,j)=0B(1,k)=i;B(2,k)=j;k=k+1;endendendB;C=C B;B=;x=C(1,:);y=C(2,:);D=size(C);b=D(:,2)-a;a=D(:,2);E=50*s*ones(1,b);z=z E;endxyz程序2.2:计算山体曲面积分Syms x yp00=-553 ;p10 =7.007;p01 = -0.2078 ;p20 = -0.03322 ;p11 = 0.02264;p02 = 0.03541;p30 = 8.848e-005 ;p21 = -0.0002331 ;p12 = -0.0001663 ;p03 = -0.0001483 ;p40 = -8.863e-008 ;p31 = 5.859e-007 ;p22 = 3.941e-007 ;p13 = 4.976e-007 ;p04 = 1.904e-007;p50 = -8.63e-011 ;p41 = -6.388e-011 ;p32 = -9.059e-010 ;p23 = -1.191e-010 ;p14 = -4.111e-010 ;p05 = -8.165e-011 ;f = p00 + p10*x + p01*y + p20*x2 + p11*x*y + p02*y2 + p30*x3 + p21*x2*y+ p12*x*y2 + p03*y3 + p40*x4 + p31*x3*y + p22*x2*y2 + p13*x*y3 + p04*y4 + p50*x5 + p41*x4*y + p32*x3*y2 + p23*x2*y3 + p14*x*y4 + p05*y5;s=100*int(int('sqrt(1+diff(f,x)2+diff(f,y)2) ',x,1,512),y,1,512)16