景區(qū)滅火數學建模模型.docx

西南交通大學2012年大學生數學建模競賽題目： B題 參賽隊員1參賽隊員2參賽隊員3姓名甘羅談盼明浩學號201038742010388720103878學院數學學院數學學院數學學院專業(yè)數學與應用數學數學與應用數學數學與應用數學電話134028869451878294300813688127427E173026319西南交通大學教務處西南交通大學實驗室及設備管理處西南交通大學數學建模創(chuàng)新實踐基地題目： 景區(qū)滅火的數學模型 摘要本文借助matlab的圖像處理功能，運用插值擬合，有限元的思想，建立了消防員救火最佳路線和消防站合理選址的優(yōu)化模型，再利用lingo軟件求得最優(yōu)解，由最優(yōu)解確定救火最佳路線和消防站合理選址。附件所給等高線圖為像素為512512的二維灰度圖.利用matlab的imread函數讀出等高線圖的灰度矩陣.根據灰度矩陣求得灰度圖中黑色點的坐標；補全等高線圖轉化為一元插值問題;通過對線性插值、三體樣條插值、三次多項式插值結果進行分析比較，三次多項式插值效果更為理想。由問題一得出的完整等高線圖可求出等高線上對應山體表面上點的三維坐標，利用matlab中的Surface Fitting Tool擬合得到景區(qū)的三維地形圖和擬合后的曲面函數。用matlab算出山體表面的曲面積分，由像素與比例尺的值可得山體表面面積。從消防站A派遣消防員去B點滅火，最佳的滅火路線應使得消防員到達火災蔓延區(qū)域邊界的時間最短。將火勢蔓延路線與消防員滅火路線合并為一條曲線，在曲線上均等插入n-1個點，得到的折線可以近似代替該曲線。折線上n-1各點的位置可以由n-1參變量唯一確定，則最短時間可由n-1個參變量唯一確定。求最佳的滅火路線轉化為轉化為以最短時間為目標函數，n-1個參變量的取值范圍為約束條件的優(yōu)化問題。利用lingo軟件可求得最優(yōu)解。根據實際情形，消防站選址的合理性由多方面因素決定。本文定義消防站的合理選址應使得火災發(fā)生對森林公園造成損失的期望最小。消防員到達著火蔓延區(qū)域邊界的時間由第三問的模型可得。再由假設4、5，可建立以火災造成損失的期望最小為目標函數的優(yōu)化模型。考慮到運算量過大，我們將景區(qū)分成了16個片區(qū)，取各個片區(qū)的中心代表該片區(qū)，近似估計最優(yōu)選址點的位置。關鍵詞：灰度矩陣 一元插值 曲面擬合 有限元法優(yōu)化模型一問題提出某國家級森林公園由于風景區(qū)植被豐富，擁有大量的國家級重點保護動植物，因此旅游管理部門在風景區(qū)的A點設置了景區(qū)消防站，當景區(qū)發(fā)生火災時能及時控制和消滅火情。圖1說明：該圖水平及豎直方向以10m為單位，山高以50m為單位。根據所給出的風景區(qū)地形等高線圖形提出了以下問題：1、由于人為原因，圖1所示的等高圖出現(xiàn)了局部破損的情況，根據所學知識利用適當的方法建立數學模型修補好該等高圖；2、在完成第一問的基礎上，結合數學模型建立該景區(qū)的三維地形圖，并估計該景區(qū)的地表面積；3、某天圖1所示的B點發(fā)生了火災，于是需要從景區(qū)消防站派遣消防員去B點滅火，建立模型確定最佳滅火路線；4、 如果需要對景區(qū)消防站進行重新選址，請建立模型確定合理的消防站地址。二問題分析2.1問題1的分析：根據所學的數學知識可以采用曲線插值的方法將等高圖中缺損的部分補充完全；2.2問題2的分析：要計算景區(qū)的地表面積，可以采用曲面積分的方法。這就要知道地表曲面的曲面方程或函數，而曲面的函數z=fx,y可以根據所給出的等高線圖上已知的點用matlab擬合出來，那么地表面積就可以用1+zy2+zy2dxdy求出；2.3問題3的分析：要確定最佳的滅火路線即確定火勢蔓延到某處與消防員到達這里的最短時間，而最短時間取決于火勢蔓延的速度和消防員行走的速度，那么就需要求出A、B兩點間的路線，設路線為L，則可以將問題轉化為優(yōu)化問題，再利用有限元的思想就可以求解出最佳路線；2.4問題4的分析：要使景區(qū)消防站更合理，義消防站的合理選址應使得火災發(fā)生對森林公園造成損失的期望最小。，那么可以將景區(qū)合理劃分成各個小片區(qū)，再求每個小片區(qū)的中心點建消防站對應的損失期望，取該值對應的片區(qū)的中心即為最佳選址。三問題假設1.假設消防員前往火災發(fā)生點的速度和火蔓延的速度只與坡度有關2.消防員在山體表面任意一點都可通行3.消防站在著火時馬上收到警報，同時消防員趕往火災地點4.景區(qū)山體表面上每一點發(fā)生火災的概率相等5.景區(qū)山體表面每一點發(fā)生火災造成的損失只消防員到達火災區(qū)的時間正相關四符號說明s:山的表面積A: 消防站B:火災發(fā)生點A0: A在水平面上的投影B0: B在水平面上的投影t:消防員到達火勢蔓延區(qū)域邊界的時間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點為消防站時火災發(fā)生造成損失的期望五模型建立與求解5.1 問題一的模型建立與求解附件所給等高線圖為像素為512512的二維灰度圖.利用matlab的imread函數可以讀出灰度圖灰度矩陣.根據灰度矩陣可以求得灰度圖中黑色點的坐標.那么等高線圖的補全就轉化為一元插值問題;對所給等高線圖進行處理,得到只包含缺失部分附近黑色點的18張圖片;將灰度矩陣的行與列分別對應坐標系中的縱坐標與橫坐標.當同一橫坐標上有多個數據點時,只保留其中一個.觀察可知,第1-16張圖片等高線缺失部分可直接由灰度矩陣得到的坐標點數據直接插值得到.對線性插值，三次抽樣插值，三次多項式插值結果比較，發(fā)現(xiàn)三次多項式插值更為合理。通過程序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將景區(qū)的三維空間曲面函數擬合出來，則可以得到景區(qū)的空間曲面圖，如圖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 問題三的模型建立與求解從景區(qū)消防站A派遣消防員去B點滅火，最佳的滅火路線應使得消防員到達火災蔓延區(qū)域的時間最短。如果消防員滅火路線與火勢蔓延經過的某條路線相交所需的時間最短，那么該路線即為消防員最佳的滅火路線。消防員滅火路線與火勢蔓延經過的曲線合為一條曲線。在曲線上插入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到達火勢蔓延區(qū)域則有可得m=m（,.）消防員t時刻在折線段PmPm+1到達火勢蔓延區(qū)域，關于t有當n取定值時，t可以表示為,.的函數，將min（t(,.)）作為目標函數，02作為約束條件求最優(yōu)解的問題；目標函數：min（t）= min（t(,.)）約束條件：1.2.3.t1i=，t2i= 4.02在、確定的情況下t1i、t2i可以表示為,.的函數，t=t(,.在02上求最優(yōu)解，利用lingo軟件，可解出,.，,.唯一確定最佳路線，對應t為消防員到達火勢蔓延區(qū)域的最短時間。5.4 問題四的模型建立與求解森林公園中某點發(fā)生火災將造成植被破壞，而植被破壞的面積與火持續(xù)燒的時間成正相關，火持續(xù)的時間主要與消防員到達的時間有關；火災造成的損失主要與植被破壞的面積有關。無論是從經濟還是生態(tài)考慮，合理的消防站選址應該使得火災發(fā)生造成的損失期望最小。根據假設4、5，火災發(fā)生造成的損失期望最小等價于消防員到達火勢蔓延區(qū)域邊界的時間t的期望最小。假設點Y：(x，y)為消防站，Y：（x,y）為山上任意一點；消防站為Y,Y發(fā)生火災時對應的t可由問題三的模型求得，記為t（x,y）消防站設為Y時對應的t的期望：E(Y)=求出山上各點的E(Y)，min（E(Y)）即為火災發(fā)生損失的期望的最小值；考慮到t（x,y）運算量過大，我們考慮將景區(qū)分成片區(qū)，取各個片區(qū)的中心代表該片區(qū)。例如我們可以等高線圖上的比例尺以128128為一個單位，將山體表面劃分為16個片區(qū)，取各個片區(qū)的中心代表該片區(qū)。求出E(Y)=求出山上各點的E(Y)，min（E(Y)）即為火災發(fā)生損失的期望的最小值；min（E(Y)）對應的（x,y）即為最佳選址地點。六評價與推廣一模型的優(yōu)缺點優(yōu)點1.模型建立過程邏輯性強、通俗易懂2.運用matlab的進行圖像處理，簡潔直觀3.問題三的模型化難為簡，利用有限元的思想以直代曲，將積分問題轉化為離散型優(yōu)化問題4.模型四從實際出發(fā)，對選址合理性的定義符合實際缺點1. matlab，lingo程序運行效率不高二模型的推廣1.模型一、模型二可用于分析統(tǒng)計數據，對位置樣本給出估計2.模型三、四貼近現(xiàn)實，對工廠，醫(yī)院等公共機構選址具有指導作用八參考文獻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=