景区灭火数学建模模型

PAGEPAGE1西南交通大学2012年大学生数学建模竞赛题目：B题参赛队员1参赛队员2参赛队员3姓名甘罗谈盼明浩学号201038742010388720103878学院数学学院数学学院数学学院专业数学与应用数学数学与应用数学数学与应用数学电话134028869451878294300813688127427Email1522407022@576946951@173026319西南交通大学教务处西南交通大学实验室及设备管理处西南交通大学数学建模创新实践基地题目：景区灭火的数学模型摘要本文借助matlab的图像处理功能，运用插值拟合，有限元的思想，建立了消防员救火最佳路线和消防站合理选址的优化模型，再利用lingo软件求得最优解，由最优解确定救火最佳路线和消防站合理选址。附件所给等高线图为像素为512512的二维灰度图.利用matlab的imread函数读出等高线图的灰度矩阵.根据灰度矩阵求得灰度图中黑色点的坐标；补全等高线图转化为一元插值问题;通过对线性插值、三体样条插值、三次多项式插值结果进行分析比较，三次多项式插值效果更为理想。由问题一得出的完整等高线图可求出等高线上对应山体表面上点的三维坐标，利用matlab中的SurfaceFittingTool拟合得到景区的三维地形图和拟合后的曲面函数。用matlab算出山体表面的曲面积分，由像素与比例尺的值可得山体表面面积。从消防站A派遣消防员去B点灭火，最佳的灭火路线应使得消防员到达火灾蔓延区域边界的时间最短。将火势蔓延路线与消防员灭火路线合并为一条曲线，在曲线上均等插入n-1个点，得到的折线可以近似代替该曲线。折线上n-1各点的位置可以由n-1参变量唯一确定，则最短时间可由n-1个参变量唯一确定。求最佳的灭火路线转化为转化为以最短时间为目标函数，n-1个参变量的取值范围为约束条件的优化问题。利用lingo软件可求得最优解。根据实际情形，消防站选址的合理性由多方面因素决定。本文定义消防站的合理选址应使得火灾发生对森林公园造成损失的期望最小。消防员到达着火蔓延区域边界的时间由第三问的模型可得。再由假设4、5，可建立以火灾造成损失的期望最小为目标函数的优化模型。考虑到运算量过大，我们将景区分成了16个片区，取各个片区的中心代表该片区，近似估计最优选址点的位置。关键词：灰度矩阵一元插值曲面拟合有限元法优化模型等高线图的补全就转化为一元插值问题;对所给等高线图进行处理,得到只包含缺失部分附近黑色点的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中的SurfaceFittingTool将景区的三维空间曲面函数拟合出来，则可以得到景区的空间曲面图，如图2.2所示，图2.2同时SurfaceFittingTool得到拟合的多项式函数为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*x^2+p11*x*y+p02*y^2+p30*x^3+p21*x^2*y+p12*x*y^2+p03*y^3+p40*x^4+p31*x^3*y+p22*x^2*y^2+p13*x*y^3+p04*y^4+p50*x^5+p41*x^4*y+p32*x^3*y^2+p23*x^2*y^3+p14*x*y^4+p05*y^5;利用曲面积分公式建立表面积模型s=100用matlab编程即可将上式面积求出，程序如附件中程序2.2，算得结果为s=26112100。5.3问题三的模型建立与求解从景区消防站A派遣消防员去B点灭火，最佳的灭火路线应使得消防员到达火灾蔓延区域的时间最短。如果消防员灭火路线与火势蔓延经过的某条路线相交所需的时间最短，那么该路线即为消防员最佳的灭火路线。消防员灭火路线与火势蔓延经过的曲线合为一条曲线。在曲线上插入n-1个点，将相邻点连接，得到一条折线，当n-1很大时，可用该折线代替曲线；以消防站B在水平面上的投影B0为原点，B0B为z轴建立右手空间直角坐标系，连接A0与B0两点，将线段A0B0均等插入n个点P0i，分别以线段Mi为半径、B0B为中轴线做n个圆柱面，与山体表面相交曲线，依次在曲线上选取折点Pi（i=1,2,3…n-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,3…n-1 的函数，当n取定值时，t1i、t2i是关于（i=0,1,2,3…n-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.模型四从实际出发，对选址合理性的定义符合实际缺点matlab，lingo程序运行效率不高二．模型的推广1.模型一、模型二可用于分析统计数据，对位置样本给出估计2.模型三、四贴近现实，对工厂，医院等公共机构选址具有指导作用八．参考文献[1]姜启源谢金星叶俊数学模型北京高等教育出版社2011.1[2]韩忠庚数学建模方法及其应用北京高等教育出版社2009[3]薛长虹于凯大学数学实验——MATLAB应用篇成都西南交通大学出版社2003附录问题一附录只包含缺失部分附近黑色点的18张图片1-16：17-18：程序1.1,对第1-16张图片插值拟合的程序clearclca=ones(512,512);imshow(a);holdonform=1:16url=strcat('C:\MATLAB7\work\',num2str(m));imageurl=strcat(url,'.bmp');A=imread(imageurl);k=1;forj=1:512fori=1:512ifA(i,j)==0B(1,k)=j;B(2,k)=i;k=k+1;break;endendendB;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])holdonclearendholdoffF1=getframe;imwrite(F1.cdata,'test1.bmp')c=imread('test1.bmp')程序1.2,用matlab读取第17-18张图片对应的灰度矩阵转置的灰度图插值拟合clearclca=ones(512,512);imshow(a);holdonform=17:18url=strcat('C:\MATLAB7\work\',num2str(m));imageurl=strcat(url,'.bmp');A=imread(imageurl);A=A';k=1;forj=1:512fori=1:512ifA(i,j)==0B(1,k)=j;B(2,k)=i;k=k+1;break;endendendB;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]);holdonclearendholdoffF1=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';fori=1:512forj=1:512ifa(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=[];fors=1:8url=strcat('C:\MATLAB7\work\second\',num2str(s));imageurl=strcat(url,'.bmp');A=imread(imageurl);k=1;fori=1:512forj=1:512ifA(i,j)==0B(1,k)=i;B(2,k)=j;k=k+1;endendendB;C=[CB];B=[];x=C(1,:);y=C(2,:);D=size(C);b=D(:,2)-a;a=D(:,2);E=50*s*ones(1,b);z=[zE];endxyz程序2.2:计算山体曲面积分Symsxyp00=-553;p10=7.00